Millepede-II V04-17-03
tinypede.py
Go to the documentation of this file.
1#!/usr/bin/env python3
2
3
24
25import numpy as np
26from scipy.stats import chi2
27import array
28import math
29import time
30
31
32
75class MilleRecord(object):
76
77
79 def __init__(self, doublePrec=False):
80
81 self.__doublePrecision = doublePrec
82
83 self.__position = 1
84
85 self.__numData = 0
86
87 self.__recLen = 0
88
89 self.__iMeas = 0
90
91 self.__iErr = 0
92
93 self.__inder = array.array('i')
94
95 self.__glder = array.array('d' if doublePrec else 'f')
96
97
101 def addData(self, dataList):
102 if (self.__numData == 0): # first word is error counter
103 self.__inder.append(0)
104 self.__glder.append(0.)
105 self.__numData += 1
106
107 aMeas, aPrec, indLocal, derLocal, labGlobal, derGlobal = dataList
108 self.__inder.append(0)
109 self.__glder.append(aMeas)
110 self.__inder.fromlist(indLocal)
111 self.__glder.fromlist(derLocal)
112 self.__inder.append(0)
113 self.__glder.append(1.0 / math.sqrt(aPrec)) # convert to error
114 self.__inder.fromlist(labGlobal)
115 self.__glder.fromlist(derGlobal)
116
117
121 def getData(self):
122 aMeas = self.__glder[self.__iMeas]
123 indLocal = []
124 derLocal = []
125 for i in range(self.__iMeas + 1, self.__iErr):
126 indLocal.append(self.__inder[i])
127 derLocal.append(self.__glder[i])
128 aPrec = 1.0 / self.__glder[self.__iErr] ** 2 # convert to precision
129 indGlobal = []
130 derGlobal = []
131 for i in range(self.__iErr + 1, self.__position):
132 indGlobal.append(self.__inder[i])
133 derGlobal.append(self.__glder[i])
134 return aMeas, aPrec, indLocal, derLocal, indGlobal, derGlobal
135
136
137 def printRecord(self):
138 print(" MilleRecord, len: ", len(self.__inder))
139 print(self.__inder)
140 print(self.__glder)
141
142
146 def writeRecord(self, aFile):
147 header = array.array('i') # header with number of words
148 header.append(-len(self.__inder) * 2 if self.__doublePrecision else len(self.__inder) * 2)
149 header.tofile(aFile)
150 self.__glder.tofile(aFile)
151 self.__inder.tofile(aFile)
152
153
158 def readRecord(self, aFile, aType):
159 if aType == 'F':
160 # FOrtran record info
161 lenf = array.array('i')
162 lenf.fromfile(aFile, 1)
163 header = array.array('i') # header with number of words
164 header.fromfile(aFile, 1)
165 self.__recLen = abs(header[0] // 2)
166 if header[0] < 0:
167 self.__glder = array.array('d')
168 else:
169 self.__glder = array.array('f')
170 self.__glder.fromfile(aFile, self.__recLen)
171 self.__inder.fromfile(aFile, self.__recLen)
172 if aType == 'F':
173 # FOrtran record info
174 lenf = array.array('i')
175 lenf.fromfile(aFile, 1)
176
177
181 def moreData(self):
182 if (self.__position < self.__recLen):
183 while (self.__position < self.__recLen and self.__inder[self.__position] != 0):
184 self.__position += 1
185 self.__iMeas = self.__position
186 self.__position += 1
187 while (self.__position < self.__recLen and self.__inder[self.__position] != 0):
188 self.__position += 1
189 self.__iErr = self.__position
190 self.__position += 1
191 # special data?
192 if (self.__iMeas + 1 == self.__iErr and self.__glder[self.__iErr] < 0):
193 self.__position += int(-self.__glder[self.__iErr])
194 else:
195 while (self.__position < self.__recLen and self.__inder[self.__position] != 0):
196 self.__position += 1
197 self.__numData += 1
198 return True
199 else:
200 return False
201
202
206 def specialDataTag(self):
207 aTag = -1
208 if (self.__iMeas + 1 == self.__iErr and self.__glder[self.__iErr] < 0):
209 aTag = int(-self.__glder[self.__iErr] * 10. + 0.5) % 10
210 return aTag
211
212
213
221class Pede(object):
222
223
225 def __init__(self):
226
228
229 self.__numPar = 0
230
232
233 self.__parIndices = {}
234
235 self.__matrix = None
236
237 self.__vector = None
238
239 self.__sumChi2 = 0.
240
241 self.__sumNdf = 0
242
243 self.__numRec = 0
244
246
247 self.__numBad = 0
248
249 self.__numCons = 0
250
251 self.__listCons = []
252
253 #
254 print()
255 print(" TinyPede - a simple PEDE implementation in python")
256 print()
257
258
267 def addBinaryFile(self, binaryFileName, maxRec=1000000, fileType='C'):
268 self.__binaryFiles.append((binaryFileName, maxRec, fileType))
269
270
279 def defineParameters(self, entriesCut=25, fixedPar=[]):
280 for fileName, maxRecord, fileType in self.__binaryFiles:
281 # open binary file
282 if fileType == 'F':
283 print(" Reading Fortran binary file ", fileName)
284 else:
285 print(" Reading binary file ", fileName)
286 binaryFile = open(fileName, "rb")
287 nRec = 0
288 try:
289
290 while nRec < maxRecord:
291 nRec += 1
292 # read record
293 rec = MilleRecord()
294 rec.readRecord(binaryFile, fileType)
295 nData = 0
296 while (rec.moreData()):
297 aTag = rec.specialDataTag()
298 if (aTag >= 0):
299 continue
300 # get data
301 nData += 1
302 indGlobal = rec.getData()[4]
303 for label in indGlobal:
304 if label not in self.__parCounters:
305 self.__parCounters[label] = 0
306 self.__parCounters[label] += 1
307
308 except EOFError:
309 pass
310
311 print(" records found ", nRec)
312 binaryFile.close()
313 # fix parameters
314 for label in fixedPar:
315 if label in self.__parCounters:
316 self.__parCounters[label] *= -1
317 # active parameters
318 for l in sorted(self.__parCounters):
319 if self.__parCounters[l] >= entriesCut:
320 self.__parIndices[l] = self.__numPar
321 self.__numPar += 1
322
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)
327 print()
328
329
337 def addConstraints(self, textFile):
338 if self.__numPar == 0:
339 print(" no variable parameters defined ")
340 return
341 # read text file
342 print(" Reading constraint file ", textFile)
343 nline = 0
344 nCon = 0
345 pairs = [] # (index, coefficient)
346 tf = open(textFile)
347 for line in iter(tf):
348 nline += 1
349 fields = line.split()
350 if len(fields) < 2:
351 continue
352 if fields[0] == 'Constraint':
353 # next constraint
354 nCon += 1
355 # previous constraint had variable parameters (pairs)?
356 if len(pairs) > 0:
357 self.__numCons += 1
358 self.__listCons.append(pairs)
359 pairs = []
360 continue
361 if fields[0] == '*':
362 continue
363 # label and coefficient
364 label = int(fields[0])
365 value = float(fields[1])
366 if label in self.__parIndices:
367 # variable parameter (label)
368 pairs.append((self.__parIndices[label], value))
369
370 tf.close()
371 # last constraint
372 if len(pairs) > 0:
373 self.__numCons += 1
374 self.__listCons.append(pairs)
375
376 print(" number of all constraints ", nCon)
377 print(" number of accepted constr. ", self.__numCons)
378 #print self.__listCons
379
380
387 def construct(self, chi2Factor=50.):
388 # create global matrix and vector
389 # size: number of parameters + number of Lagrange multipliers for constraints
390 # no initial non-zero values for global parameters (simplifying local fit, Lagrange multipliers)
391 self.__vector = np.zeros(self.__numPar + self.__numCons)
392 self.__matrix = np.zeros((self.__numPar + self.__numCons, self.__numPar + self.__numCons))
393 print()
394 print(" Constructing linear equations system")
395 print(" local chi2 scaling factor ", chi2Factor)
396 for fileName, maxRecord, fileType in self.__binaryFiles:
397 # open binary file
398 binaryFile = open(fileName, "rb")
399 nRec = 0
400 try:
401 while nRec < maxRecord:
402 # read record
403 rec = MilleRecord()
404 rec.readRecord(binaryFile, fileType)
405 nRec += 1
406 dataBlocks = []
407 # number of local parameters
408 numLoc = 0
409 # number of measurements
410 numMeas = 0
411
412 while (rec.moreData()):
413 aTag = rec.specialDataTag()
414 if (aTag >= 0):
415 continue
416 # get data
417 dataBlocks.append(rec.getData())
418 numLoc = max(numLoc, dataBlocks[-1][2][-1])
419 numMeas += 1
420
421 #print(" numLoc ", nRec, numLoc, numMeas)
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))
428 # loop over data blocks
429 for iData, data in enumerate(dataBlocks):
430 #aRes = data.getResidual()[0]
431 aMeas, aPrec, indLocal, derLocal, indGlobal, derGlobal = data
432 # update local-local matrix
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
438 # update local-global matrix
439 for l, der in zip(indGlobal, derGlobal):
440 if l not in self.__parIndices:
441 continue
442 i = self.__parIndices[l]
443 gloDer[iData, i] = der
444 gloDerW[iData, i] = der * aPrec
445
446 # local fit
447 matLocLoc = np.dot(locDerW.T, locDer)
448 vecLoc = np.dot(locDerW.T, vecMeas)
449 try:
450 locCov = np.linalg.inv(matLocLoc)
451 except np.linalg.LinAlgError as e:
452 print(" local fit failed:", e)
453 self.__numBad += 1
454 continue
455 locSol = np.dot(locCov, vecLoc)
456 # chi2, ndf
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)
460 # check rank
461 nRank = np.linalg.matrix_rank(matLocLoc)
462 if nRank < numLoc:
463 self.__numBad += 1
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]) # global correlation (squared)
467 print(" locPar ", i + 1, matLocLoc[i, i], locCov[i, i], gcor2)
468 continue
469 # reject (at 3 sigma)
470 if prob < 0.0027:
471 self.__numRejects += 1
472 continue
473 # accepted
474 self.__sumNdf += locNdf
475 self.__sumChi2 += locChi2
476 # update global matrix, vector (parameter part)
477 self.__matrix[:self.__numPar,:self.__numPar] += np.dot(gloDerW.T, gloDer)
478 self.__vector[:self.__numPar] += np.dot(gloDerW.T, vecMeas)
479 # account for correlations (via local parameters)
480 matGloLoc = np.dot(gloDerW.T, locDer)
481 self.__matrix[:self.__numPar,:self.__numPar] -= np.dot(matGloLoc, np.dot(locCov, matGloLoc.T))
482 self.__vector[:self.__numPar] -= np.dot(matGloLoc, locSol)
483
484 except EOFError:
485 pass
486
487 binaryFile.close()
488 self.__numRec += nRec
489
490 # add Lagrange multipliers
491 indexCons = self.__numPar
492 for pairs in self.__listCons:
493 for idx, val in pairs:
494 self.__matrix[indexCons, idx] = val
495 self.__matrix[idx, indexCons] = val
496 indexCons += 1
497
498 # record statistics
499 print(" total number of records ", self.__numRec)
500 print(" number of rejected records ", self.__numRejects)
501 print(" number of bad records ", self.__numBad)
502 # matrix information
503 matSize = self.__numPar + self.__numCons
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), '%')
507
508
509 def dump(self):
510 print(" Global matrix, non-zero elements")
511 for i in range(self.__numPar + self.__numCons):
512 for j in range(i, self.__numPar + self.__numCons):
513 if self.__matrix[i, j] != 0.:
514 print(" ", i, j, self.__matrix[i, j])
515
516
520 def solve(self):
521 print()
522 print(" Solving linear equations system")
523 # solve equation system
524 try:
525 aCov = np.linalg.inv(self.__matrix)
526 except np.linalg.LinAlgError as e:
527 print(" matrix inversion failed:", e)
528 print(" >>> improve input data <<<")
529 return
530 aSol = np.dot(aCov, self.__vector)
531 print(" sum(ndf) ", self.__sumNdf - self.__numPar)
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)
536 print(" final sum(chi2)/sum(ndf) ", (self.__sumChi2 - deltaChi2) / (self.__sumNdf - self.__numPar))
537 print()
538 print(" Solution")
539 print(" parameter label, #entries, correction, error")
540 for l in sorted(self.__parCounters):
541 if l in self.__parIndices:
542 i = self.__parIndices[l]
543 print(" ", l, self.__parCounters[l], aSol[i], math.sqrt(aCov[i, i]) if aCov[i, i] > 0. else -math.sqrt(-aCov[i, i]))
544 else:
545 print(" ", l, self.__parCounters[l], " fixed")
546
547
548if __name__ == '__main__':
549 print(time.asctime())
550 #
551 # Some steering parameters
552 #
553 # maximum number of records to read (for each file)
554 maxRec = 10000
555 # minimum number of entries (appearances in binary files) for variable global parameters
556 entries = 25
557
558 #
559 # get pede instance
560 #
561 p = Pede()
562
563 #
564 # MP2 test case (from pede -t)
565 #
566 # add (Fortran) binary file(s)
567 p.addBinaryFile("mp2tst.bin", maxRec, 'F')
568 # define parameters
569 p.defineParameters(entries, fixedPar=[])
570 # add constraints
571 p.addConstraints("mp2con.txt")
572 # construct linear equation system
573 p.construct()
574 #p.dump()
575 # solve linear equation system
576 p.solve()
577
578 print()
579 print(time.asctime())
Millepede-II (binary) record.
Definition: tinypede.py:75
def readRecord(self, aFile, aType)
Read record from file.
Definition: tinypede.py:158
__glder
array with values, errors and derivatives; (float32 or float64)
Definition: tinypede.py:95
def printRecord(self)
Print record.
Definition: tinypede.py:137
__inder
array with markers (0) and labels; array(int32)
Definition: tinypede.py:93
__iMeas
position of value in current data block; int
Definition: tinypede.py:89
__recLen
record length; int
Definition: tinypede.py:87
__iErr
position of error in current data block; int
Definition: tinypede.py:91
__numData
number of data blocks in record; int
Definition: tinypede.py:85
def moreData(self)
Locate next data block.
Definition: tinypede.py:181
def getData(self)
Get data block from current position in record.
Definition: tinypede.py:121
def __init__(self, doublePrec=False)
Create MP-II binary record.
Definition: tinypede.py:79
def writeRecord(self, aFile)
Write record to file.
Definition: tinypede.py:146
def addData(self, dataList)
Add data block to (end of) record.
Definition: tinypede.py:101
def specialDataTag(self)
Get special data tag from block.
Definition: tinypede.py:206
__position
position in record, usually start of next data block; int
Definition: tinypede.py:83
__doublePrecision
flag for storage in as double values
Definition: tinypede.py:81
__parIndices
parameter indices
Definition: tinypede.py:233
def solve(self)
Solve linear equation system.
Definition: tinypede.py:520
def construct(self, chi2Factor=50.)
Construct linear equation system.
Definition: tinypede.py:387
def defineParameters(self, entriesCut=25, fixedPar=[])
Define Parameters.
Definition: tinypede.py:279
__numRejects
number of rejects (Chi2,ndf)
Definition: tinypede.py:245
__sumChi2
sum(chi2)
Definition: tinypede.py:239
def addConstraints(self, textFile)
Add constraints.
Definition: tinypede.py:337
def addBinaryFile(self, binaryFileName, maxRec=1000000, fileType='C')
Add binary file.
Definition: tinypede.py:267
__numCons
rank of global matrix
Definition: tinypede.py:249
__sumNdf
sum(ndf)
Definition: tinypede.py:241
__vector
vector
Definition: tinypede.py:237
__binaryFiles
binary files
Definition: tinypede.py:227
__numRec
number of records
Definition: tinypede.py:243
__parCounters
parameter counters
Definition: tinypede.py:231
__numBad
bad local fits (rank deficit)
Definition: tinypede.py:247
__listCons
list constraints
Definition: tinypede.py:251
def __init__(self)
Constructor.
Definition: tinypede.py:225
__matrix
matrix
Definition: tinypede.py:235
def dump(self)
Dump global matrix.
Definition: tinypede.py:509
__numPar
number of parameters
Definition: tinypede.py:229