Skip to content

Commit

Permalink
Checkpointing work. Refs #7653.
Browse files Browse the repository at this point in the history
  • Loading branch information
@wdzhou
wdzhou committed Aug 27, 2013
1 parent f496e9e commit 05bedc7
Show file tree
Hide file tree
Showing 2 changed files with 265 additions and 13 deletions.
127 changes: 126 additions & 1 deletion Code/Mantid/Framework/CurveFitting/test/NeutronBk2BkExpConvPVoigtTest.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,7 +103,7 @@ class NeutronBk2BkExpConvPVoigtTest : public CxxTest::TestSuite
//----------------------------------------------------------------------------------------------
/** Calculate peak positions: data is from Fullprof's sample: arg_si
*/
void test_calculatePeakShape()
void TODO_test_calculatePeakShape()
{
NeutronBk2BkExpConvPVoigt func;
func.initialize();
Expand Down Expand Up @@ -156,6 +156,131 @@ class NeutronBk2BkExpConvPVoigtTest : public CxxTest::TestSuite
*/
}

//----------------------------------------------------------------------------------------------
/** Calculate peak positions: data is from Fullprof's sample: arg_si
*/
void TODO_test_calculateVulcanPeakPositions()
{
// (2, 2, 0)
NeutronBk2BkExpConvPVoigt func;
func.initialize();

func.setParameter("Dtt1", 16370.650);
func.setParameter("Dtt2", 0.100);
func.setParameter("Zero", 0.000);
func.setParameter("LatticeConstant", 5.431363); // Silicon

func.setMillerIndex(3, 3, 1);
func.calculateParameters(false);
double dh1 = func.getPeakParameter("d_h");
double tofh1 = func.centre();

cout << "Peak [111]: d_h = " << dh1 << ", TOF_h = " << tofh1 << ".\n";

TS_ASSERT_DELTA(tofh1, 23421.7207, 0.01);

/*
// (2, 2, 0)
NeutronBk2BkExpConvPVoigt func220;
func220.initialize();
func220.setParameter("Dtt1", 16370.650);
func220.setParameter("Dtt2", -1.540);
func220.setParameter("Zero", -9.227);
func220.setParameter("LatticeConstant", 5.431363); // Silicon
func220.setMillerIndex(2, 2, 0);
func220.calculateParameters(false);
double tofh2 = func220.centre();
TS_ASSERT_DELTA(tofh2, 14342.8350, 0.01);
// (3,1,1)
NeutronBk2BkExpConvPVoigt func311;
func311.initialize();
func311.setParameter("Dtt1", 7476.910);
func311.setParameter("Dtt2", -1.540);
func311.setParameter("Zero", -9.227);
func311.setParameter("LatticeConstant", 5.431363); // Silicon
func311.setMillerIndex(3, 1, 1);
func311.calculateParameters(false);
double tofh3 = func311.centre();
TS_ASSERT_DELTA(tofh3, 12230.9648, 0.01);
// (2, 2, 2)
NeutronBk2BkExpConvPVoigt func222;
func222.initialize();
func222.setParameter("Dtt1", 7476.910);
func222.setParameter("Dtt2", -1.540);
func222.setParameter("Zero", -9.227);
func222.setParameter("LatticeConstant", 5.431363); // Silicon
func222.setMillerIndex(2, 2, 2);
func222.calculateParameters(false);
double tofh4 = func222.centre();
TS_ASSERT_DELTA(tofh4, 11710.0332, 0.01);
*/

return;
}

//----------------------------------------------------------------------------------------------
/** Calculate peak positions: data is from Fullprof's sample: arg_si
*/
void test_calculateVulcanProfile()
{
NeutronBk2BkExpConvPVoigt func;
func.initialize();

func.setParameter("Dtt1", 16370.650);
func.setParameter("Dtt2", 0.100);
func.setParameter("Zero", 0.000);
func.setParameter("LatticeConstant", 5.431363); // Silicon

func.setParameter("Alph0", 1.000000);
func.setParameter("Alph1", 0.000000);
func.setParameter("Beta0", 0.109036);
func.setParameter("Beta1", 0.009834);
func.setParameter("Sig0", sqrt(0.000));
func.setParameter("Sig1", sqrt(1119.230));
func.setParameter("Sig2", sqrt(91.127));
func.setParameter("Gam0", 0.000);
func.setParameter("Gam1", 2.604);
func.setParameter("Gam2", 0.000);

// Peak 220
func.setMillerIndex(2, 2, 0);
func.calculateParameters(false);

// Peak centre
double tofh1 = func.centre();
TS_ASSERT_DELTA(tofh1, 23421.7207, 0.01);

// Peak shape
func.setParameter("Height", 1.0);

double fwhm = func.fwhm();
TS_ASSERT_DELTA(fwhm, 47.049, 0.001);

cout << "Peak 220: TOF_h = " << tofh1 << ", FWHM = " << fwhm << ".\n";

vector<double> vecX;
double tof = tofh1 - 10*fwhm;
while (tof < tofh1 + 10*fwhm)
{
vecX.push_back(tof);
tof += fwhm*0.1;
}
vector<double> vecY(vecX.size(), 0.0);
func.function(vecY, vecX);
for (size_t i = 0; i < vecX.size(); ++i)
cout << vecX[i] << "\t\t" << vecY[i] << "\n";

return;
}


};


Expand Down
151 changes: 139 additions & 12 deletions Code/Mantid/Framework/PythonInterface/plugins/algorithms/LoadFullprofFile.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,7 @@ def PyInit(self):
"""
self.setWikiSummary("""Load file generated by Fullprof.""")

self.declareProperty(FileProperty("Filename","", FileAction.Load, ['.hkl', '.prf']),
self.declareProperty(FileProperty("Filename","", FileAction.Load, ['.hkl', '.prf', '.dat']),
"Name of [http://www.ill.eu/sites/fullprof/ Fullprof] .hkl or .prf file.")

#self.declareProperty("Bank", 1, "Bank ID for output if there are more than one bank in .irf file.")
Expand All @@ -76,6 +76,9 @@ def PyExec(self):
elif fpfilename.lower().endswith(".prf") is True:
# (.prf) file
self._tableWS, self._dataWS= self._loadFullprofPrfFile(fpfilename)
elif fpfilename.lower().endswith(".dat") is True:
# (.dat) file: Fullprof data file
self._tableWS, self._dataWS = self._loadFullprofDataFile(fpfilename)
else:
raise NotImplementedError("File %s is neither .hkl nor .prf. It is not supported." % (fpfilename))

Expand Down Expand Up @@ -231,7 +234,9 @@ def _loadFullprofPrfFile(self, prffilename):
def _parseFullprofPrfFile(self, filename):
""" Parse Fullprof .prf file to a information dictionary and a data set (list of list)
"""
# 1. Import file
import re

# Import .prf file
try:
pfile = open(filename, "r")
except IOError:
Expand All @@ -247,7 +252,8 @@ def _parseFullprofPrfFile(self, filename):
infodict = {}
dataset = []

# 2. Parse information line
# Parse information lines
# Line 0: header
infoline = lines[0]
if infoline.count("CELL:") == 1:
terms = infoline.split("CELL:")[1].split()
Expand All @@ -268,24 +274,39 @@ def _parseFullprofPrfFile(self, filename):
# spacegroup = terms[1].strip()
# infodict["SpaceGroup"] = spacegroup

# 3. Parse data
# Find data line header
firstline = -1
for i in xrange(1, len(lines)):
if lines[i].startswith("T.O.F"):
if lines[i].count("Yobs-Ycal") > 0:
firstline = i+1
dataheader = lines[i].strip()
break

if firstline < 0:
raise NotImplementedError("File format is incorrect. Unable to locate data title line")

# FIXME - Need to parse header line: T.O.F. Yobs Ycal Yobs-Ycal Backg Bragg ...
# to determine how the data line look alike (==5 or >= 5)
# Parse header line: T.O.F. Yobs Ycal Yobs-Ycal Backg Bragg ...
# to determine how the data line look alike (==5 or >= 5)
headerterms = dataheader.split()
dataperline = 5
# TOF., ... h k l ...
reflectionperline = len(headerterms)-5+3


# Parse data
count = 0
for i in xrange(firstline, len(lines)):
terms = lines[i].split()
if (len(terms) == 5 and lines[i].count(")") == 0) or (len(terms) > 5 and terms[2].count("(") == 0):
# it is a data line
line = lines[i].strip()
if len(line) == 0: # empty line
continue

if line.count(")") == 0 and line.count("(") == 0:
# A regular data line
terms = line.split()
if len(terms) != 5:
self.log().warning("Pure data line %d (%s) has irregular number of data points" % (i, line))
continue

x = float(terms[0])
yobs = float(terms[1])
ycal = float(terms[2])
Expand All @@ -295,8 +316,38 @@ def _parseFullprofPrfFile(self, filename):
dataset.append([x, yobs, ycal, ydif, ybak])
count += 1

else:
break
elif line.count(")") == 1 and line.count("(") == 1:
# A line can be either pure reflection line or a combined data/reflection line
# remove '(' and ')'
newline = re.sub('[()]', ' ', line)

terms = newline.strip().split()
if len(terms) < 9:
# Pure reflection line
tofh = float(terms[0])
hklstr = line.split(")")[1].split(")")[0].strip()
infodict[hklstr] = tofh

else:
# Mixed line: least number of items: data(5) + TOF+hkl = 9

x = float(terms[0])
yobs = float(terms[1])
ycal = float(terms[2])
ydif = float(terms[3])
ybak = float(terms[4])

dataset.append([x, yobs, ycal, ydif, ybak])
count += 1

raise NotImplementedError("Need a sample line of this use case.")
hklstr = line.split(")")[1].split(")")[0].strip()
infodict[hklstr] = tofh
# ENDIFELSE (terms)
else:
self.log().warning("%d-th line (%s) is not well-defined." % (i, line))

# ENDIF-ELIF-ELSE (line.count)
# ENDFOR

print "Data set counter = ", count
Expand All @@ -311,6 +362,82 @@ def _makeEmptyDataWorkspace(self):

return dataws

def _loadFullprofDataFile(self, datafilename):
""" Parse a Fullprof (multiple) column file
"""
# Import file
datafile = open(datafilename, "r")
rawlines = datafile.readlines()
datafile.close()

# Parse head
iline = 0
parseheader = True
title = ""
while iline < len(rawlines) and parseheader is True:
line = rawlines[iline].strip()
if len(line) > 0:
if line.count("BANK") == 0:
# header
title += line + ", "
else:
# line information
terms = line.split()
if terms[0] != 'BANK':
raise NotImplementedError("First word must be 'BANK', but not %s" % (terms[0]))
bankid = int(terms[1])
numdata = int(terms[2])
numlines = int(terms[3])

parseheader = False
# ENDIF
# ENDIF
iline += 1
# ENDWHILE (iline)

# Data vectors
vecx = []
vecy = []
vece = []

for i in xrange(iline, len(rawlines)):
line = rawlines[i].strip()
if len(line) == 0:
continue

terms = line.split()
numitems = len(terms)
if numitems % 3 != 0:
print "%d-th line '%s' is not a data line" % (i, line)
continue

numpts = numitems/3
for j in xrange(numpts):
x = float(terms[j*3])
y = float(terms[j*3+1])
e = float(terms[j*3+2])

vecx.append(x)
vecy.append(y)
vece.append(e)
# ENDFOR
# ENDFOR (i)

# Check
self.log().notice("Expected to read %d data points; Exactly read %d data points. " % (numdata*numlines, len(vecx)))

# Create output workspaces
tablews = WorkspaceFactory.createTable()

# Create the data workspace
datasize = len(vecx)
dataws = WorkspaceFactory.create("Workspace2D", 1, datasize, datasize)
for i in xrange(datasize):
dataws.dataX(0)[i] = vecx[i]
dataws.dataY(0)[i] = vecy[i]
dataws.dataE(0)[i] = vece[i]

return (tablews, dataws)

# Register algorithm with Mantid
#registerAlgorithm(LoadFullprofFile)
Expand Down

0 comments on commit 05bedc7

Please sign in to comment.