26from scipy.stats
import chi2
95 self.
__glder = array.array(
'd' if doublePrec
else 'f')
107 aMeas, aPrec, indLocal, derLocal, labGlobal, derGlobal = dataList
110 self.
__inder.fromlist(indLocal)
111 self.
__glder.fromlist(derLocal)
113 self.
__glder.append(1.0 / math.sqrt(aPrec))
114 self.
__inder.fromlist(labGlobal)
115 self.
__glder.fromlist(derGlobal)
126 indLocal.append(self.
__inder[i])
127 derLocal.append(self.
__glder[i])
132 indGlobal.append(self.
__inder[i])
133 derGlobal.append(self.
__glder[i])
134 return aMeas, aPrec, indLocal, derLocal, indGlobal, derGlobal
138 print(
" MilleRecord, len: ", len(self.
__inder))
147 header = array.array(
'i')
161 lenf = array.array(
'i')
162 lenf.fromfile(aFile, 1)
163 header = array.array(
'i')
164 header.fromfile(aFile, 1)
167 self.
__glder = array.array(
'd')
169 self.
__glder = array.array(
'f')
174 lenf = array.array(
'i')
175 lenf.fromfile(aFile, 1)
255 print(
" TinyPede - a simple PEDE implementation in python")
268 self.
__binaryFiles.append((binaryFileName, maxRec, fileType))
283 print(
" Reading Fortran binary file ", fileName)
285 print(
" Reading binary file ", fileName)
286 binaryFile = open(fileName,
"rb")
290 while nRec < maxRecord:
294 rec.readRecord(binaryFile, fileType)
296 while (rec.moreData()):
297 aTag = rec.specialDataTag()
302 indGlobal = rec.getData()[4]
303 for label
in indGlobal:
311 print(
" records found ", nRec)
314 for label
in fixedPar:
323 print(
" Found in binary files")
324 print(
" total number of parameters ", len(self.
__parCounters))
325 print(
" minimum number of entries ", entriesCut)
326 print(
" number of free parameters ", self.
__numPar)
339 print(
" no variable parameters defined ")
342 print(
" Reading constraint file ", textFile)
347 for line
in iter(tf):
349 fields = line.split()
352 if fields[0] ==
'Constraint':
364 label =
int(fields[0])
365 value = float(fields[1])
376 print(
" number of all constraints ", nCon)
377 print(
" number of accepted constr. ", self.
__numCons)
394 print(
" Constructing linear equations system")
395 print(
" local chi2 scaling factor ", chi2Factor)
398 binaryFile = open(fileName,
"rb")
401 while nRec < maxRecord:
404 rec.readRecord(binaryFile, fileType)
412 while (rec.moreData()):
413 aTag = rec.specialDataTag()
417 dataBlocks.append(rec.getData())
418 numLoc = max(numLoc, dataBlocks[-1][2][-1])
422 vecMeas = np.zeros(numMeas)
423 vecMeasW = np.zeros(numMeas)
424 locDer = np.zeros((numMeas, numLoc))
425 locDerW = np.zeros((numMeas, numLoc))
426 gloDer = np.zeros((numMeas, self.
__numPar))
427 gloDerW = np.zeros((numMeas, self.
__numPar))
429 for iData, data
in enumerate(dataBlocks):
431 aMeas, aPrec, indLocal, derLocal, indGlobal, derGlobal = data
433 for i, der
in zip(indLocal, derLocal):
434 locDer[iData, i - 1] = der
435 locDerW[iData, i - 1] = der * aPrec
436 vecMeas[iData] = aMeas
437 vecMeasW[iData] = aMeas * aPrec
439 for l, der
in zip(indGlobal, derGlobal):
443 gloDer[iData, i] = der
444 gloDerW[iData, i] = der * aPrec
447 matLocLoc = np.dot(locDerW.T, locDer)
448 vecLoc = np.dot(locDerW.T, vecMeas)
450 locCov = np.linalg.inv(matLocLoc)
451 except np.linalg.LinAlgError
as e:
452 print(
" local fit failed:", e)
455 locSol = np.dot(locCov, vecLoc)
457 locNdf = numMeas - numLoc
458 locChi2 = np.dot(vecMeas - np.dot(locDer, locSol), vecMeasW - np.dot(locDerW, locSol))
459 prob = 1. - chi2.cdf(locChi2 / chi2Factor, locNdf)
461 nRank = np.linalg.matrix_rank(matLocLoc)
464 print(
" ??? bad local fit ", nRec, locNdf, locChi2, numLoc, nRank)
465 for i
in range(numLoc):
466 gcor2 = 1.0 - 1.0 / (matLocLoc[i, i] * locCov[i, i])
467 print(
" locPar ", i + 1, matLocLoc[i, i], locCov[i, i], gcor2)
480 matGloLoc = np.dot(gloDerW.T, locDer)
493 for idx, val
in pairs:
499 print(
" total number of records ", self.
__numRec)
500 print(
" number of rejected records ", self.
__numRejects)
501 print(
" number of bad records ", self.
__numBad)
504 nonZero = np.count_nonzero(self.
__matrix)
505 print(
" size of global matrix ", matSize)
506 print(
" fraction of elements <> 0. ",
int(100.*nonZero / matSize ** 2),
'%')
510 print(
" Global matrix, non-zero elements")
514 print(
" ", i, j, self.
__matrix[i, j])
522 print(
" Solving linear equations system")
526 except np.linalg.LinAlgError
as e:
527 print(
" matrix inversion failed:", e)
528 print(
" >>> improve input data <<<")
532 print(
" initial sum(chi2) ", self.
__sumChi2)
533 deltaChi2 = np.dot(self.
__vector, aSol)
534 print(
" chi2 reduction by fit ", deltaChi2)
535 print(
" final sum(chi2) ", self.
__sumChi2 - deltaChi2)
539 print(
" parameter label, #entries, correction, error")
543 print(
" ", l, self.
__parCounters[l], aSol[i], math.sqrt(aCov[i, i])
if aCov[i, i] > 0.
else -math.sqrt(-aCov[i, i]))
548if __name__ ==
'__main__':
549 print(time.asctime())
567 p.addBinaryFile(
"mp2tst.bin", maxRec,
'F')
569 p.defineParameters(entries, fixedPar=[])
571 p.addConstraints(
"mp2con.txt")
579 print(time.asctime())
Millepede-II (binary) record.
def readRecord(self, aFile, aType)
Read record from file.
__glder
array with values, errors and derivatives; (float32 or float64)
def printRecord(self)
Print record.
__inder
array with markers (0) and labels; array(int32)
__iMeas
position of value in current data block; int
__recLen
record length; int
__iErr
position of error in current data block; int
__numData
number of data blocks in record; int
def moreData(self)
Locate next data block.
def getData(self)
Get data block from current position in record.
def __init__(self, doublePrec=False)
Create MP-II binary record.
def writeRecord(self, aFile)
Write record to file.
def addData(self, dataList)
Add data block to (end of) record.
def specialDataTag(self)
Get special data tag from block.
__position
position in record, usually start of next data block; int
__doublePrecision
flag for storage in as double values
__parIndices
parameter indices
def solve(self)
Solve linear equation system.
def construct(self, chi2Factor=50.)
Construct linear equation system.
def defineParameters(self, entriesCut=25, fixedPar=[])
Define Parameters.
__numRejects
number of rejects (Chi2,ndf)
def addConstraints(self, textFile)
Add constraints.
def addBinaryFile(self, binaryFileName, maxRec=1000000, fileType='C')
Add binary file.
__numCons
rank of global matrix
__binaryFiles
binary files
__numRec
number of records
__parCounters
parameter counters
__numBad
bad local fits (rank deficit)
__listCons
list constraints
def __init__(self)
Constructor.
def dump(self)
Dump global matrix.
__numPar
number of parameters