GeneralBrokenLines V03-01-02
gblpy3
gblfit.py
Go to the documentation of this file.
1'''
2Track fit with general broken lines.
3
4Created on Jul 27, 2011
5
6@author: kleinwrt
7'''
8
9
29
30import numpy as np
31import math
32from .gblnum import BorderedBandMatrix
33from .mille import MilleRecord
34
35
36
44class GblMeasurement(object):
45
46
55 def __init__(self, aMeasurement, minPrecision=0.):
56
57 self.__measurement = aMeasurement
58
59 self.__measDim = aMeasurement[1].shape[0]
60
61 self.__measMinPrec = minPrecision
62
64
66
67 self.__globalLabels = None
68
70
71 # full precision matrix?
72 if (aMeasurement[2].ndim == 2):
73 # need to diagonalize
74 eigenVal, eigenVec = np.linalg.eigh(aMeasurement[2])
75 self.__measTransformation = eigenVec.T
76 # transform measurement
77 if (aMeasurement[0] is None):
79 else:
80 self.__measurement[0] = np.dot(self.__measTransformation, aMeasurement[0])
81 self.__measurement[1] = np.dot(self.__measTransformation, aMeasurement[1])
82 self.__measurement[2] = eigenVal
83
84
88 def hasMeasurement(self):
89
90 return (self.__measurement is not None)
91
92
96 def getMeasurement(self):
97 return self.__measurement
98
99
103 def getMeasDim(self):
104 return self.__measDim
105
106
110 def getMeasMinPrec(self):
111 return self.__measMinPrec
112
113
117 def addLocals(self, derivatives):
118 if (self.__measDim > 0):
119 if (self.__measTransformation is None):
120 self.__localDerivatives = derivatives
121 else:
122 self.__localDerivatives = np.dot(self.__measTransformation, derivatives)
123
124
129 def addGlobals(self, labels, derivatives):
130 if (self.__measDim > 0):
131 self.__globalLabels = labels
132 if (self.__measTransformation is None):
133 self.__globalDerivatives = derivatives
134 else:
135 self.__globalDerivatives = np.dot(self.__measTransformation, derivatives)
136
137
141 def getNumLocals(self):
142 if (self.__localDerivatives is not None):
143 return self.__localDerivatives.shape[1]
144 else:
145 return 0
146
147
152 return self.__localDerivatives
153
154
158 return self.__globalLabels
159
160
165 return self.__globalDerivatives
166
167
169 print(" measurement ", self.__measDim, len(self.__localDerivatives), len(self.__globalDerivatives))
170
171#------------------------------------------------------------------------------
172
173
174
182class GblPoint(object):
183
184
188 def __init__(self, aJacobian):
189
190 self.__label = 0
191
192 self.__offset = 0
193
194 self.__type = 0
195
196 self.__p2pJacobian = aJacobian
197
198 self.__jacobians = [ [], [] ]
199
201
202 self.__scatDim = 0
203
204 self.__scatterer = None
205
206# for extension to retrieval of residuals, pulls
207# self.__dataMeas = [0, 0]
208# for extension to retrieval of residuals, pulls
209# self.__dataScat = [0, 0]
210
211
220 def addMeasurement(self, aMeasurement, minPrecision=0.):
221 self.__measurements.append(GblMeasurement(aMeasurement, minPrecision))
222
223
227 def hasMeasurement(self):
228 return (self.__measurements != [])
229
230
235 return self.__measurements
236
237
284 def addScatterer(self, aScatterer):
285 self.__scatDim = aScatterer[0].shape[0] # 2 (thin) or 4 (thick)
286 self.__scatterer = [ None ] + aScatterer
287 if (aScatterer[1].ndim == 2): # full precision matrix, need to diagonalize
288 eigenVal, eigenVec = np.linalg.eigh(aScatterer[1])
289 scatTransformation = eigenVec.T
290# transform measurement
291 self.__scatterer[0] = scatTransformation
292 self.__scatterer[1] = np.dot(scatTransformation, aScatterer[0])
293 self.__scatterer[2] = eigenVal
294
295
299 def getScatDim(self):
300 return self.__scatDim
301
302
306 def getScatterer(self):
307 return self.__scatterer
308
309
316 # get scatterer input back from diagonalization
317 scatRes = np.dot(self.__scatterer[0].T, self.__scatterer[1])
318 # positional covariance from diagonalization
319 posCov = np.dot(np.dot(self.__scatterer[0][:, 2:].T, np.diag(1. / self.__scatterer[2])), self.__scatterer[0][:, 2:])
320
321 scatPrec = np.zeros((4, 4))
322 scatPrec[2:, 2:] = np.linalg.inv(posCov)
323 # diagonalize again
324 eigenVal, eigenVec = np.linalg.eigh(scatPrec)
325 scatTransformation = eigenVec.T
326 return [scatTransformation, np.dot(scatTransformation, scatRes), eigenVal]
327
328
332 def addLocals(self, derivatives):
333 if self.__measurements != []:
334 self.__measurements[-1].addLocals(derivatives)
335
336
341 def addGlobals(self, labels, derivatives):
342 if self.__measurements != []:
343 self.__measurements[-1].addGlobals(labels, derivatives)
344
345
349 def setLabel(self, aLabel):
350 self.__label = aLabel
351
352
356 def getLabel(self):
357 return self.__label
358
359
363 def setOffset(self, anOffset):
364 self.__offset = anOffset
365
366
370 def getOffset(self):
371 return self.__offset
372
373
377 def setType(self, aType):
378 self.__type = aType
379
380
384 def isFirst(self):
385 return (self.__type < 0)
386
387
391 def isLast(self):
392 return (self.__type > 0)
393
394# def setDataMeas(self, aIndex, aData):
395# for extension to retrieval of residuals, pulls
396# self.__dataMeas[aIndex] = aData
397
398# def setDataScat(self, aIndex, aData):
399# for extension to retrieval of residuals, pulls
400# self.__dataScat[aIndex] = aData
401
402
406 def getP2pJacobian(self):
407 return self.__p2pJacobian
408
409
413 def addPrevJacobian(self, aJacobian):
414 self.__jacobians[0] = np.linalg.inv(aJacobian)
415
416
420 def addNextJacobian(self, aJacobian):
421 self.__jacobians[1] = aJacobian
422
423
428 def getDerivatives(self, index):
429 aJacobian = self.__jacobians[index]
430 matJ = aJacobian[3:5, 3:5] # J
431 matS = aJacobian[3:5, 1:3] # S
432 vecD = aJacobian[3:5, 0:1] # d
433 if (index < 1):
434 matS = -matS
435 matW = np.linalg.inv(matS) # W = +/- S^-1
436 return matW, np.dot(matW, matJ), np.dot(matW, vecD) # W, W*J, W*d
437
438
439 def printPoint(self):
440 print(" point ", self.__label, self.__offset, self.__scatDim, len(self.__measurements))
441
442#------------------------------------------------------------------------------
443
444
445
449class GblData(object):
450
451
459 def __init__(self, aLabel=0, aType=0, aValue=0., aPrec=0., aMeas=0):
460
461 self.__label = aLabel
462
463 self.__type = aType
464
465 self.__index = aMeas
466
467 self.__value = aValue
468
469 self.__precision = aPrec
470
471 self.__dwMethod = 0
472
473 self.__downWeight = 1.
474
475 self.__prediction = 0.
476
477 self.__parameters = []
478
480
482
484
485# self.__predictionVar = 0.
486
487
496 def addDerivatives(self, iRow, labDer, matDer, \
497 derLocal=None, labGlobal=None, derGlobal=None):
498 if (derLocal is not None):
499 for i in range(derLocal.shape[1]): # local derivatives
500 if (derLocal[iRow, i] != 0.):
501 self.__derivatives.append(derLocal[iRow, i])
502 self.__parameters.append(i + 1)
503
504 for i in range(len(labDer)): # curvature, offset derivatives
505 if (labDer[i] != 0 and matDer[iRow, i] != 0.):
506 self.__derivatives.append(matDer[iRow, i])
507 self.__parameters.append(labDer[i])
508
509 if (derGlobal is not None):
510 for i in range(derGlobal.shape[1]): # global derivatives
511 if (derGlobal[iRow, i] != 0.):
512 self.__globalLabels.append(labGlobal[iRow, i])
513 self.__globalDerivatives.append(derGlobal[iRow, i])
514
515
520 def addExtDerivatives(self, indexExt, derExt):
521 for i in range(len(derExt)): # external derivatives
522 if (derExt[i] != 0.):
523 self.__derivatives.append(derExt[i])
524 self.__parameters.append(indexExt[i])
525
526
530 def getMatrices(self):
531 aVector = np.array([ self.__derivatives ])
532 aMatrix = np.dot(aVector.T, aVector)
533 aValue = self.__value
534 aWeight = self.__precision * self.__downWeight
535 return self.__parameters, aValue * aVector * aWeight, aMatrix * aWeight
536
537
541 def setPrediction(self, aVector):
542 self.__prediction = 0.
543 for i in range(len(self.__parameters)):
544 self.__prediction += self.__derivatives[i] * aVector[ self.__parameters[i] - 1 ]
545
546# def setPredictionVariance(self, aMatrix):
547# '''Calculate variance of prediction for data from fit.'''
548# var(residual) = 1./precision**2 - var(prediction)
549# aBlockMatrix = aMatrix.getBlockMatrix(self.__parameters)
550# self.__predictionVar = np.dot(self.__derivatives.T, \
551# np.dot(aBlockMatrix, self.__derivatives))
552
553
558 def setDownWeighting(self, aMethod):
559 self.__dwMethod = aMethod
560 scaledResidual = abs(self.__value - self.__prediction) * math.sqrt(self.__precision)
561 if (aMethod == 1): # Tukey
562 if (scaledResidual < 4.6851):
563 aWeight = (1.0 - (scaledResidual / 4.6851) ** 2) ** 2
564 else:
565 aWeight = 0.
566 elif (aMethod == 2): # Huber
567 if (scaledResidual < 1.345):
568 aWeight = 1.
569 else:
570 aWeight = 1.345 / scaledResidual
571 elif (aMethod == 3): # Cauchy
572 aWeight = 1.0 / (1.0 + (scaledResidual / 2.3849) ** 2)
573 self.__downWeight = aWeight
574 return aWeight
575
576
582 def getChi2(self):
583 scaledResidual = abs(self.__value - self.__prediction) * math.sqrt(self.__precision)
584 Chi2 = scaledResidual ** 2
585 if (self.__dwMethod == 1): # Tukey
586 if (scaledResidual < 4.6851):
587 Chi2 = 4.6851 ** 2 / 3. * (1. - (1. - (scaledResidual / 4.6851) ** 2) ** 3)
588 else:
589 Chi2 = 4.6851 ** 2 / 3.
590 elif (self.__dwMethod == 2): # Huber
591 if (scaledResidual > 1.345):
592 Chi2 = 1.345 * (2.*scaledResidual - 1.345)
593 elif (self.__dwMethod == 3): # Cauchy
594 Chi2 = math.log(1. + (scaledResidual / 2.3849) ** 2) * 2.3849 ** 2
595 return Chi2
596
597
601 def getLabel(self):
602 return self.__label
603
604
608 def getType(self):
609 return self.__type
610
611
615 def getIndex(self):
616 return self.__index
617
618
622 def getResidual(self):
623 return self.__value - self.__prediction, 1.0 / self.__precision, self.__downWeight, self.__parameters, self.__derivatives
624
625
629 def toRecord(self):
630 return self.__value, self.__precision, self.__parameters, self.__derivatives, \
632
633
637 def fromRecord(self, dataList):
638 self.__value, self.__precision, self.__parameters, self.__derivatives, \
639 self.__globalLabels, self.__globalDerivatives = dataList
640
641
647 def analyzeData(self, maxBand):
648 maxPar = self.__parameters[-1]
649 maxBor = 0
650 for i in self.__parameters:
651 if (i < maxPar - maxBand):
652 maxBor = i
653 return maxPar, maxBor
654
655
656 def printData(self):
657 print(" measurement at label ", self.__label, " with type ", self.__type, " : ", self.__value, self.__precision)
658 print(" param ", self.__parameters)
659 print(" deriv ", self.__derivatives)
660 print(" global labels ", self.__globalLabels)
661 print(" global deriv ", self.__globalDerivatives)
662
663#------------------------------------------------------------------------------
664
665
758
759
760
762class GblTrajectory(object):
763
764
769 def __init__(self, hasCurv=True, aDim=[0, 1]):
770
771 self.__numPoints = 0
772
774
775 self.__numCurvature = (1 if hasCurv else 0)
776
778
779 self.__numLocals = 0
780
782
783 self.__dimensions = aDim
784
785 self.__points = []
786
787 self.__data = []
788
789 self.__externalSeed = None
790
792
794
796
797
802 def addPoint(self, point):
803 self.__numPoints += 1
804 label = self.__numPoints
805 point.setLabel(label)
806 self.__points.append(point)
807 for m in point.getMeasurements():
808 self.__numLocals = max(self.__numLocals, m.getNumLocals())
809 return label
810
811
815 def getNumPoints(self):
816 return self.__numPoints
817
818
823 def addExternalSeed(self, aLabel, aSeed):
824 self.__externalPoint = aLabel
825 self.__externalSeed = aSeed
826
827
828 def printPoints(self):
829 print("GblPoints")
830 for p in self.__points:
831 p.printPoint()
832
833
834 def printData(self):
835 print("GblData blocks")
836 for d in self.__data:
837 d.printData()
838
839
843 def getData(self):
844 return self.__data
845
846
851 def milleOut(self, aFile, doublePrec=False):
852 rec = MilleRecord(doublePrec)
853# data measurements and kinks
854 for aData in self.__data:
855 rec.addData(aData.toRecord())
856
857 rec.writeRecord(aFile)
858
859
863 def milleIn(self, aFile):
864 rec = MilleRecord()
865 rec.readRecord(aFile)
866 mPar = 0
867 mBor = 0
868 mBand = 4 * len(self.__dimensions) - 1 # max band width
869 while (rec.moreData()):
870 aTag = rec.specialDataTag()
871 if (aTag < 0):
872# get data
873 aData = GblData()
874 aData.fromRecord(rec.getData())
875 self.__data.append(aData)
876 nPar, nBor = aData.analyzeData(mBand)
877 mPar = max(mPar, nPar)
878 mBor = max(mBor, nBor)
879
880 self.__numParameters = mPar
881 self.__numLocals = mBor - self.__numCurvature
882
883
888 def __getJacobian(self, aLabel):
889 aDim = self.__dimensions
890 nDim = len(aDim)
891 anIndex = abs(aLabel) - 1
892# check consistency of (index, direction)
893 if (aLabel > 0):
894 nJacobian = 1
895 if (anIndex >= self.__numPoints - 1):
896 anIndex = self.__numPoints - 1
897 nJacobian = 0
898 else:
899 nJacobian = 0
900 if (anIndex <= 0):
901 anIndex = 0
902 nJacobian = 1
903# Jacobian broken lines (q/p,..,u_i,u_i+1..) to local (q/p,u',u) parameters
904 nCurv = self.__numCurvature
905 nLocals = self.__numLocals
906 nBorder = nCurv + nLocals
907 nParBRL = nBorder + 2 * nDim
908 nParLoc = nLocals + 5
909 aJacobian = np.zeros((nParLoc, nParBRL))
910 aPoint = self.__points[anIndex]
911 anIndex = []
912 labDer, matDer = self.__getFitToLocalJacobian(aPoint, 5, nJacobian)
913# from local parameters
914 for i in range(nLocals):
915 aJacobian[i + 5, i] = 1.0;
916 anIndex.append(i + 1);
917
918# from trajectory parameters
919 iCol = nLocals;
920 for i in range(5):
921 if (labDer[i] > 0):
922 anIndex.append(labDer[i]);
923 for j in range(5):
924 aJacobian[j, iCol] = matDer[j, i];
925 iCol += 1
926
927 return anIndex, aJacobian
928
929
939 def __getFitToLocalJacobian(self, aPoint, measDim, nJacobian=1):
940 aDim = self.__dimensions
941 nDim = len(aDim)
942 nCurv = self.__numCurvature
943 nLocals = self.__numLocals
944 nOffset = aPoint.getOffset()
945 scatDim = aPoint.getScatDim()
946 anIndex = [0, 0, 0, 0, 0]
947 aJacobian = np.zeros((measDim, 5))
948 labOffset = measDim - 2
949 labSlope = measDim - 4
950 labCurv = measDim - 5
951
952 if (nOffset < 0): # need interpolation
953 prevW, prevWJ, prevWd = aPoint.getDerivatives(0) # W-, W- * J-, W- * d-
954 nextW, nextWJ, nextWd = aPoint.getDerivatives(1) # W+, W+ * J+, W+ * d+
955 sumWJ = prevWJ + nextWJ
956 matN = np.linalg.inv(sumWJ) # N = (W- * J- + W+ * J+)^-1
957# local offset
958 if (labOffset >= 0):
959# derivatives for u_int
960 prevNW = np.dot(matN, prevW) # N * W-
961 nextNW = np.dot(matN, nextW) # N * W+
962 prevNd = np.dot(matN, prevWd) # N * W- * d-
963 nextNd = np.dot(matN, nextWd) # N * W+ * d+
964 iOff = nDim * (-nOffset - 1) + nLocals + nCurv + 1 # first offset ('i' in u_i)
965 if (nCurv > 0):
966 aJacobian[labOffset:measDim, 0:1] = -prevNd - nextNd # from curvature
967 anIndex[0] = nLocals + 1
968 aJacobian[labOffset:measDim, 1:3] = prevNW
969 aJacobian[labOffset:measDim, 3:5] = nextNW
970 for i in range(nDim):
971 anIndex[1 + aDim[i]] = iOff + i
972 anIndex[3 + aDim[i]] = iOff + nDim + i
973# local slope
974 if (labSlope >= 0):
975# derivatives for u'_int
976 prevWPN = np.dot(nextWJ, prevNW) # W+ * J+ * N * W-
977 nextWPN = np.dot(prevWJ, nextNW) # W- * J- * N * W+
978 prevWNd = np.dot(nextWJ, prevNd) # W+ * J+ * N * W- * d-
979 nextWNd = np.dot(prevWJ, nextNd) # W- * J- * N * W+ * d+
980 if (nCurv > 0):
981 aJacobian[labSlope:labOffset, 0:1] = prevWNd - nextWNd # from curvature
982 aJacobian[labSlope:labOffset, 1:3] = -prevWPN
983 aJacobian[labSlope:labOffset, 3:5] = nextWPN
984
985 else: # at point
986# anIndex must be sorted
987# forward : iOff2 = iOff1 + nDim, index1 = 1, index2 = 3
988# backward: iOff2 = iOff1 - nDim, index1 = 3, index2 = 1
989 # adjust for thick scatterer
990 if (scatDim == 4):
991 nOffset += 1
992 iOff1 = nDim * nOffset + nCurv + nLocals + 1 # first offset ('i' in u_i)
993 index1 = 3 - 2 * nJacobian # index of first offset
994 iOff2 = iOff1 + nDim * (nJacobian * 2 - 1) # second offset ('i' in u_i)
995 index2 = 1 + 2 * nJacobian # index of second offset
996# local offset
997 if (labOffset >= 0):
998 aJacobian[labOffset, index1] = 1.0 # from 1st Offset
999 aJacobian[labOffset + 1, index1 + 1] = 1.0
1000 for i in range(nDim):
1001 anIndex[index1 + aDim[i]] = iOff1 + i
1002# local slope and curvature
1003 if (labSlope >= 0):
1004 matW, matWJ, vecWd = aPoint.getDerivatives(nJacobian) # W, W * J, W * d
1005 sign = 2 * nJacobian - 1
1006 if (nCurv > 0):
1007 aJacobian[labSlope:labOffset, 0:1] = -sign * vecWd # from curvature
1008 anIndex[0] = nLocals + 1
1009 aJacobian[labSlope:labOffset, index1:index1 + 2] = -sign * matWJ
1010 aJacobian[labSlope:labOffset, index2:index2 + 2] = sign * matW
1011 for i in range(nDim):
1012 anIndex[index2 + aDim[i]] = iOff2 + i
1013
1014# local curvature
1015 if (labCurv >= 0):
1016 if (nCurv > 0):
1017 aJacobian[labCurv, labCurv] = 1.0
1018
1019 return anIndex, aJacobian
1020
1021
1029 def __getFitToKinkJacobian(self, aPoint):
1030 aDim = self.__dimensions
1031 nDim = len(aDim)
1032 nCurv = self.__numCurvature
1033 nLocals = self.__numLocals
1034 nOffset = aPoint.getOffset()
1035 anIndex = [0, 0, 0, 0, 0, 0, 0]
1036 aJacobian = np.zeros((2, 7))
1037
1038 prevW, prevWJ, prevWd = aPoint.getDerivatives(0) # W-, W- * J-, W- * d-
1039 nextW, nextWJ, nextWd = aPoint.getDerivatives(1) # W+, W+ * J+, W+ * d+
1040 sumWJ = prevWJ + nextWJ # W- * J- + W+ * J+
1041 sumWd = prevWd + nextWd # W+ * d+ + W- * d-
1042 iOff = (nOffset - 1) * nDim + nCurv + nLocals + 1 # first offset ('i' in u_i)
1043
1044# local offset
1045 if (nCurv > 0):
1046 aJacobian[:, 0:1] = -sumWd # from curvature
1047 anIndex[0] = nLocals + 1
1048 aJacobian[:, 1:3] = prevW # from 1st Offset
1049 aJacobian[:, 3:5] = -sumWJ # from 2nd Offset
1050 aJacobian[:, 5:7] = nextW # from 3rd Offset
1051 for i in range(nDim):
1052 anIndex[1 + aDim[i]] = iOff + i
1053 anIndex[3 + aDim[i]] = iOff + nDim + i
1054 anIndex[5 + aDim[i]] = iOff + nDim * 2 + i
1055
1056 return anIndex, aJacobian
1057
1058
1068 def __getFitToStepJacobian(self, aPoint):
1069 aDim = self.__dimensions
1070 nDim = len(aDim)
1071 nCurv = self.__numCurvature
1072 nLocals = self.__numLocals
1073 nOffset = aPoint.getOffset()
1074 anIndex = [0, 0, 0, 0]
1075 aJacobian = np.zeros((4, 4))
1076
1077 iOff = (nOffset - 1) * nDim + nCurv + nLocals + 1 # first offset ('i' in u_i)
1078
1079# step
1080 aJacobian[2:,:2] = -np.eye(2) # from 2nd Offset
1081 aJacobian[2:, 2:] = np.eye(2) # from 3rd Offset
1082
1083 for i in range(nDim):
1084 anIndex[aDim[i]] = iOff + nDim + i
1085 anIndex[2 + aDim[i]] = iOff + nDim * 2 + i
1086
1087 return anIndex, aJacobian
1088
1089
1100 aDim = self.__dimensions
1101 nDim = len(aDim)
1102 nCurv = self.__numCurvature
1103 nLocals = self.__numLocals
1104 nOffset = aPoint.getOffset()
1105 anIndex = [0, 0, 0, 0, 0, 0, 0, 0, 0]
1106 aJacobian = np.zeros((4, 9))
1107
1108 prevW, prevWJ, prevWd = aPoint.getDerivatives(0) # W-, W- * J-, W- * d-
1109 nextW, nextWJ, nextWd = aPoint.getDerivatives(1) # W+, W+ * J+, W+ * d+
1110 iOff = (nOffset - 1) * nDim + nCurv + nLocals + 1 # first offset ('i' in u_i)
1111
1112# kink
1113 if (nCurv > 0):
1114 aJacobian[:2, 0:1] = -(prevWd + nextWd) # from curvature
1115 anIndex[0] = nLocals + 1
1116 aJacobian[:2, 1:3] = prevW # from 1st Offset
1117 aJacobian[:2, 3:5] = -prevWJ # from 2nd Offset
1118 aJacobian[:2, 5:7] = -nextWJ # from 3rd Offset
1119 aJacobian[:2, 7:9] = nextW # from 4th Offset
1120# step
1121 aJacobian[2:, 3:5] = -np.eye(2) # from 2nd Offset
1122 aJacobian[2:, 5:7] = np.eye(2) # from 3rd Offset
1123
1124 for i in range(nDim):
1125 anIndex[1 + aDim[i]] = iOff + i
1126 anIndex[3 + aDim[i]] = iOff + nDim + i
1127 anIndex[5 + aDim[i]] = iOff + nDim * 2 + i
1128 anIndex[7 + aDim[i]] = iOff + nDim * 3 + i
1129
1130 return anIndex, aJacobian
1131
1132 def __getResAndErr(self, aData, used=True):
1133 aResidual, aMeasVar, aDownWeight, indLocal, derLocal = self.__data[aData].getResidual()
1134 aVec = np.array(derLocal) # compressed vector
1135 aMat = self.__matrix.getBlockMatrix(indLocal) # compressed matrix
1136 aFitVar = np.dot(aVec, np.dot(aMat, aVec.T)) # variance from track fit
1137 aFitVar *= aDownWeight # account for down-weighting (of measurement in fit)
1138 aMeasError = math.sqrt(aMeasVar) # error of measurement
1139 if used:
1140 aResError = math.sqrt(aMeasVar - aFitVar) if aFitVar < aMeasVar else 0. # error of biased residual
1141 else:
1142 aResError = math.sqrt(aMeasVar + aFitVar) # error of unbiased residual
1143 return aResidual, aMeasError, aResError, aDownWeight
1144
1145
1154 def getResults(self, aLabel):
1155 anIndex, aJacobian = self.__getJacobian(aLabel)
1156 nParBRL = len(anIndex)
1157 aVec = np.empty(nParBRL)
1158 for i in range(nParBRL):
1159 aVec[i] = self.__vector[anIndex[i] - 1] # compressed vector
1160 aMat = self.__matrix.getBlockMatrix(anIndex) # compressed matrix
1161 locPar = np.dot(aJacobian, aVec)
1162 locCov = np.dot(aJacobian, np.dot(aMat, aJacobian.T))
1163 return locPar, locCov
1164
1165
1172 def getMeasResults(self, aLabel):
1173 firstData = self.__measDataIndex[aLabel - 1] # first data block with measurement
1174 numData = self.__measDataIndex[aLabel] - firstData # number of data blocks
1175 if numData <= 0:
1176 return 0, None, None, None, None
1177 #
1178 aResiduals = np.empty(numData)
1179 aMeasErr = np.empty(numData)
1180 aResErr = np.empty(numData)
1181 aDownWeight = np.empty(numData)
1182 for i in range(numData):
1183 aResiduals[i], aMeasErr[i], aResErr[i], aDownWeight[i] = self.__getResAndErr(firstData + i, (aLabel != self.__skippedMeasLabel))
1184 return numData, aResiduals, aMeasErr, aResErr, aDownWeight
1185
1186
1193 def getScatResults(self, aLabel):
1194 firstData = self.__scatDataIndex[aLabel - 1] # first data block with measurement
1195 numData = self.__scatDataIndex[aLabel] - firstData # number of data blocks
1196 if numData <= 0:
1197 return 0, None, None, None, None
1198 #
1199 aResiduals = np.empty(numData)
1200 aMeasErr = np.empty(numData)
1201 aResErr = np.empty(numData)
1202 aDownWeight = np.empty(numData)
1203 for i in range(numData):
1204 aResiduals[i], aMeasErr[i], aResErr[i], aDownWeight[i] = self.__getResAndErr(firstData + i)
1205 return numData, aResiduals, aMeasErr, aResErr, aDownWeight
1206
1207
1211 def getResidual(self, iData):
1212 return self.__getResAndErr(iData)
1213
1214
1220 def fit(self, optionList="", aLabel=0):
1221
1222
1223 def defineOffsets():
1224# set labels for previous/next offsets
1225# first point is offset
1226 self.__points[0].setOffset(0)
1227 self.__points[0].setType(-1)
1228 nOffsets = 1
1229# intermediate scatterers are offsets
1230 for aPoint in self.__points[1:-1]:
1231 scatDim = aPoint.getScatDim()
1232 if (scatDim > 0):
1233 aPoint.setOffset(nOffsets)
1234 nOffsets += 1
1235 # thick scatterer ?
1236 if (scatDim == 4):
1237 nOffsets += 1
1238 else:
1239 aPoint.setOffset(-nOffsets)
1240# last point is offset
1241 self.__points[-1].setOffset(nOffsets)
1242 self.__points[-1].setType(1)
1243 if (self.__points[-1].getScatDim() == 4):
1244 # thick scatterer
1245 nOffsets += 1
1246 self.__numOffsets = nOffsets + 1
1247 self.__numParameters = self.__numOffsets * len(self.__dimensions) \
1248 +self.__numCurvature + self.__numLocals
1249
1250
1251 def calcJacobians():
1252 scatJacobian = np.empty((5, 5))
1253# forward propagation (all)
1254 lastPoint = 0;
1255 numStep = 0;
1256 for iPoint in range(1, self.__numPoints):
1257 if (numStep == 0):
1258 scatJacobian = self.__points[iPoint].getP2pJacobian()
1259 else:
1260 scatJacobian = np.dot(self.__points[iPoint].getP2pJacobian(), scatJacobian)
1261 numStep += 1
1262 self.__points[iPoint].addPrevJacobian(scatJacobian) # aPoint -> previous scatterer
1263 if (self.__points[iPoint].getOffset() >= 0):
1264 self.__points[lastPoint].addNextJacobian(scatJacobian) # lastPoint -> next scatterer
1265 numStep = 0;
1266 lastPoint = iPoint
1267# backward propagation (without scatterers)
1268 for iPoint in range(self.__numPoints - 1, 0, -1):
1269 if (self.__points[iPoint].getOffset() >= 0):
1270 scatJacobian = self.__points[iPoint].getP2pJacobian()
1271 continue # skip offsets
1272 self.__points[iPoint].addNextJacobian(scatJacobian) # iPoint -> next scatterer
1273 scatJacobian = np.dot(scatJacobian, self.__points[iPoint].getP2pJacobian())
1274
1275
1276 def prepare():
1277 aDim = self.__dimensions
1278# measurements
1279 self.__measDataIndex.append(len(self.__data)) # offset
1280 for aPoint in self.__points:
1281 if (aPoint.hasMeasurement()):
1282 nLabel = aPoint.getLabel()
1283 measIndex = 0
1284 for m in aPoint.getMeasurements():
1285 measIndex += 1
1286 measDim = m.getMeasDim()
1287 measPrecision = m.getMeasMinPrec()
1288 localDer = m.getLocalDerivatives()
1289 globalLab = m.getGlobalLabels()
1290 globalDer = m.getGlobalDerivatives()
1291 matP, aMeas, aPrec = m.getMeasurement()
1292 nJacobian = 1 if aPoint.getOffset() < self.__numOffsets - 1 else 0 # last point needs backward propagation
1293 labDer, matDer = self.__getFitToLocalJacobian(aPoint, measDim, nJacobian)
1294 matPDer = matDer if matP is None else np.dot(matP, matDer)
1295 for i in range(measDim):
1296 if (aPrec[i] > measPrecision):
1297 aData = GblData(nLabel, 1, aMeas[i], aPrec[i], measIndex)
1298 aData.addDerivatives(i, labDer, matPDer, localDer, \
1299 globalLab, globalDer)
1300 self.__data.append(aData)
1301 self.__measDataIndex.append(len(self.__data))
1302# aPoint.setDataMeas(i, len(self.__data))
1303# pseudo measurements from kinks
1304 self.__scatDataIndex.append(len(self.__data)) # offset
1305 self.__scatDataIndex.append(len(self.__data)) # first point
1306 for aPoint in self.__points[1:]:
1307 scatDim = aPoint.getScatDim()
1308 if (scatDim > 0):
1309 nLabel = aPoint.getLabel()
1310 if (scatDim == 4):
1311 # thick scatterer, last point?
1312 if (aPoint.isLast()):
1313 # steps (only)
1314 matT, aMeas, aPrec = aPoint.getReducedScatterer()
1315 labDer, matDer = self.__getFitToStepJacobian(aPoint)
1316 else:
1317 # kinks+steps
1318 matT, aMeas, aPrec = aPoint.getScatterer()
1319 labDer, matDer = self.__getFitToKinkAndStepJacobian(aPoint)
1320 else:
1321 # thin scatterer, last point?
1322 if (aPoint.isLast()):
1323 continue
1324 else:
1325 # kinks
1326 matT, aMeas, aPrec = aPoint.getScatterer()
1327 labDer, matDer = self.__getFitToKinkJacobian(aPoint)
1328 matTDer = matDer if matT is None else np.dot(matT, matDer)
1329 for i in aDim:
1330 if (aPrec[i] > 0.):
1331 aData = GblData(nLabel, 2, aMeas[i], aPrec[i])
1332 aData.addDerivatives(i, labDer, matTDer)
1333 self.__data.append(aData)
1334 # thick scatterer?
1335 if (scatDim == 4):
1336 for i in aDim:
1337 if (aPrec[i + 2] > 0.):
1338 aData = GblData(nLabel, 2, aMeas[i + 2], aPrec[i + 2])
1339 aData.addDerivatives(i + 2, labDer, matTDer)
1340 self.__data.append(aData)
1341 self.__scatDataIndex.append(len(self.__data))
1342# aPoint.setDataScat(i, len(self.__data))
1343 self.__scatDataIndex.append(len(self.__data)) # last point
1344# external seed
1345 if (self.__externalPoint != 0):
1346 externalIndex, aJacobian = self.__getJacobian(self.__externalPoint)
1347 eigenVal, eigenVec = np.linalg.eigh(self.__externalSeed)
1348 aMatrix = np.dot(eigenVec.T, aJacobian)
1349 for i in range(len(eigenVec)):
1350 if (eigenVal[i] > 0.):
1351 externalDerivatives = []
1352 for j in range(len(externalIndex)):
1353 externalDerivatives.append(aMatrix[i, j])
1354 aData = GblData(self.__externalPoint, 3, 0., eigenVal[i])
1355 aData.addExtDerivatives(externalIndex, externalDerivatives)
1356 self.__data.append(aData)
1357 self.__measDataIndex.append(len(self.__data))
1358
1359
1360 def buildLinearEquationSystem():
1361 nBorder = self.__numCurvature + self.__numLocals
1363 self.__vector = np.zeros(self.__numParameters)
1364 for aData in self.__data:
1365 # skipped (internal) measurement?
1366 if aData.getLabel() == self.__skippedMeasLabel and aData.getType() == 1:
1367 continue
1368 index, aVector, aMatrix = aData.getMatrices()
1369 for i in range(len(index)):
1370 self.__vector[ index[i] - 1 ] += aVector[0, i] # update vector
1371 self.__matrix.addBlockMatrix(index, aMatrix) # update matrix
1372
1373
1378 def downWeight(aMethod):
1379 aLoss = 0.
1380 for aData in self.__data:
1381 aLoss += (1. - aData.setDownWeighting(aMethod))
1382 return aLoss
1383
1384
1385 def predict():
1386 for aData in self.__data:
1387 aData.setPrediction(self.__vector)
1388
1389 if (self.__data == []): # generate data from points
1390 defineOffsets()
1391 calcJacobians()
1392 prepare()
1393
1394 # skip measurements from point?
1395 self.__skippedMeasLabel = aLabel
1396
1397 buildLinearEquationSystem() # create linear equations system from data
1398#
1399 try:
1400 aMethod = 0
1401 lostWeight = 0.
1402 self.__vector = self.__matrix.solveAndInvertBorderedBand(self.__vector)
1403 predict()
1404
1405 for o in optionList: # down weighting iterations
1406 try:
1407 aMethod = "THC".index(o.upper()) + 1
1408 lostWeight = downWeight(aMethod)
1409 buildLinearEquationSystem()
1410 self.__vector = self.__matrix.solveAndInvertBorderedBand(self.__vector)
1411 predict()
1412 except ValueError:
1413 pass
1414
1415 Ndf = -self.__numParameters
1416 Chi2 = 0.
1417 for aData in self.__data:
1418 # skipped (internal) measurement?
1419 if aData.getLabel() == self.__skippedMeasLabel and aData.getType() == 1:
1420 continue
1421 Chi2 += aData.getChi2()
1422 Ndf += 1
1423 Chi2 /= [1.0, 0.8737, 0.9326, 0.8228 ][aMethod]
1424 return Chi2, Ndf, lostWeight
1425
1426 except (ZeroDivisionError, np.linalg.linalg.LinAlgError):
1427 return 0., -1, 0.
1428
Data (block) containing value, precision and derivatives for measurements, kinks and external seed.
Definition: gblfit.py:449
__dwMethod
down weighting method; int
Definition: gblfit.py:471
__globalDerivatives
derivatives (prediction vs global parameters); list(float)
Definition: gblfit.py:483
def toRecord(self)
Get data components (for copying to MP binaty record)
Definition: gblfit.py:629
__parameters
labels of fit parameters (with non zero derivative); list(int)
Definition: gblfit.py:477
def setDownWeighting(self, aMethod)
Outlier down weighting with M-estimators.
Definition: gblfit.py:558
def setPrediction(self, aVector)
Calculate prediction for data from fit.
Definition: gblfit.py:541
def analyzeData(self, maxBand)
Analyze labels of fit parameters to determine number of parameters and border size with given maximal...
Definition: gblfit.py:647
def getIndex(self)
Get index.
Definition: gblfit.py:615
def addExtDerivatives(self, indexExt, derExt)
Add derivatives to data (block) from external seed.
Definition: gblfit.py:520
__index
index for internal measurement
Definition: gblfit.py:465
def getMatrices(self)
Calculate compressed matrix and right hand side from data.
Definition: gblfit.py:530
def getType(self)
Get type.
Definition: gblfit.py:608
__value
value (residual or kink); float
Definition: gblfit.py:467
__downWeight
down weighting factor (M-estimators); float
Definition: gblfit.py:473
def printData(self)
Print data.
Definition: gblfit.py:656
def __init__(self, aLabel=0, aType=0, aValue=0., aPrec=0., aMeas=0)
Create new data.
Definition: gblfit.py:459
__globalLabels
labels of global parameters; list(int)
Definition: gblfit.py:481
__label
label of corresponding point; int
Definition: gblfit.py:461
__derivatives
derivatives (prediction vs fit parameters); list(float)
Definition: gblfit.py:479
def getLabel(self)
Get Label.
Definition: gblfit.py:601
def getChi2(self)
Calculate Chi2 (contribution) from data.
Definition: gblfit.py:582
def fromRecord(self, dataList)
Set data components (from MP binaty record)
Definition: gblfit.py:637
__precision
precision (diagonal element of inverse covariance matrix); float
Definition: gblfit.py:469
def addDerivatives(self, iRow, labDer, matDer, derLocal=None, labGlobal=None, derGlobal=None)
Add derivatives to data (block) from measurements and kinks.
Definition: gblfit.py:497
__type
type of data (0: none, 1: internal measurement, 2: internal kink, 3: external seed,...
Definition: gblfit.py:463
def getResidual(self)
Get data for residual (and errors).
Definition: gblfit.py:622
__prediction
prediction (for value from fit); float
Definition: gblfit.py:475
User supplied measurement at point on (initial) trajectory.
Definition: gblfit.py:44
def getMeasDim(self)
Retrieve dimension of measurement of a point.
Definition: gblfit.py:103
def addLocals(self, derivatives)
Add local derivatives.
Definition: gblfit.py:117
def getMeasurement(self)
Retrieve measurement of a point.
Definition: gblfit.py:96
__globalLabels
global labels; matrix(int)
Definition: gblfit.py:67
def getMeasMinPrec(self)
Retrieve minimal precision to accept measurement.
Definition: gblfit.py:110
__measTransformation
transformation (to eigen-vectors of precision matrix); matrix(float)
Definition: gblfit.py:63
def printMeasurement(self)
Print Measurement.
Definition: gblfit.py:168
__globalDerivatives
global derivatives; matrix(float)
Definition: gblfit.py:69
__measMinPrec
minimal precision to accept measurement <
Definition: gblfit.py:61
__measDim
dimension of measurement (2D, 4D or 5D); int
Definition: gblfit.py:59
__localDerivatives
local derivatives; matrix(float)
Definition: gblfit.py:65
def getNumLocals(self)
Get number of local derivatives.
Definition: gblfit.py:141
def __init__(self, aMeasurement, minPrecision=0.)
Create new measurement.
Definition: gblfit.py:55
def hasMeasurement(self)
Check point for a measurement.
Definition: gblfit.py:88
def getLocalDerivatives(self)
Get local derivatives.
Definition: gblfit.py:151
def addGlobals(self, labels, derivatives)
Add global derivatives.
Definition: gblfit.py:129
__measurement
measurement at point: projection (dm/du), residuals (to initial trajectory), precision; list(matrix(f...
Definition: gblfit.py:57
def getGlobalDerivatives(self)
Get global derivatives.
Definition: gblfit.py:164
def getGlobalLabels(self)
Get global labels.
Definition: gblfit.py:157
User supplied point on (initial) trajectory.
Definition: gblfit.py:182
def hasMeasurement(self)
Check point for a measurement.
Definition: gblfit.py:227
def addNextJacobian(self, aJacobian)
Add jacobian to next offset.
Definition: gblfit.py:420
__measurements
measurements at point; list(GblMeasurement)
Definition: gblfit.py:200
def printPoint(self)
Print point.
Definition: gblfit.py:439
def getMeasurements(self)
Retrieve measurements of a point.
Definition: gblfit.py:234
def setOffset(self, anOffset)
Define offset of a point and references to previous and next point with offset.
Definition: gblfit.py:363
def addScatterer(self, aScatterer)
Add a (thin or thick) scatterer to a point.
Definition: gblfit.py:284
def addPrevJacobian(self, aJacobian)
Add jacobian to previous offset.
Definition: gblfit.py:413
def __init__(self, aJacobian)
Create new point.
Definition: gblfit.py:188
def addLocals(self, derivatives)
Add local derivatives (to last measurement).
Definition: gblfit.py:332
__label
label for referencing point (0,1,..,number(points)-1); int
Definition: gblfit.py:190
def getLabel(self)
Retrieve label of a point.
Definition: gblfit.py:356
__scatterer
scatterer at point: (transformation or None,) initial kinks, precision (inverse covariance matrix); l...
Definition: gblfit.py:204
def addMeasurement(self, aMeasurement, minPrecision=0.)
Add a mesurement to a point.
Definition: gblfit.py:220
def getScatDim(self)
Get scatterer dimension.
Definition: gblfit.py:299
def setLabel(self, aLabel)
Define label of a point.
Definition: gblfit.py:349
def getScatterer(self)
Retrieve scatterer of a point.
Definition: gblfit.py:306
def addGlobals(self, labels, derivatives)
Add global derivatives (to last measurement).
Definition: gblfit.py:341
def getDerivatives(self, index)
Get derivatives for locally linearized track model (backward or forward propagation).
Definition: gblfit.py:428
__offset
>=0: offset number at point, <0: offset number at next point with offset; int
Definition: gblfit.py:192
def isFirst(self)
Is first point?
Definition: gblfit.py:384
__jacobians
jacobians for propagation to previous or next point with offsets; pair(matrix(float))
Definition: gblfit.py:198
def getP2pJacobian(self)
Retrieve jacobian of a point.
Definition: gblfit.py:406
__p2pJacobian
Point-to-point jacobian from previous point; matrix(float)
Definition: gblfit.py:196
def isLast(self)
Is last point?
Definition: gblfit.py:391
def getOffset(self)
Get offset of a point.
Definition: gblfit.py:370
__scatDim
scatterer dimension (valid 2(thin) or 4(thick))
Definition: gblfit.py:202
__type
type (-1: first, 0: inner, 1: last)
Definition: gblfit.py:194
def getReducedScatterer(self)
Retrieve (reduced) scatterer of a point.
Definition: gblfit.py:315
def setType(self, aType)
Define type of a point.
Definition: gblfit.py:377
General Broken Lines Trajectory.
Definition: gblfit.py:762
def __getFitToLocalJacobian(self, aPoint, measDim, nJacobian=1)
Definition: gblfit.py:939
def milleIn(self, aFile)
Read (data blocks of) trajectory from MP (binary) file.
Definition: gblfit.py:863
__numLocals
number of local parameters; int
Definition: gblfit.py:779
__numCurvature
'curvature' is fit parameter (=1); int
Definition: gblfit.py:775
def getNumPoints(self)
Get number of points on trajectory.
Definition: gblfit.py:815
__numOffsets
number of (points with) offsets on trajectory; int
Definition: gblfit.py:773
__dimensions
active components of offsets (both ([0,1]) or single ([0] or [1]); list(int)
Definition: gblfit.py:783
__externalSeed
external seed (for local, fit parameters); matrix(float)
Definition: gblfit.py:789
__measDataIndex
mapping points to data blocks from measurements; list(int)
Definition: gblfit.py:791
def getResidual(self, iData)
Get residuals from data of trajectory.
Definition: gblfit.py:1211
def __getJacobian(self, aLabel)
Get jacobian for transformation from fit to track parameters at point.
Definition: gblfit.py:888
def getMeasResults(self, aLabel)
Get residuals at point from measurement.
Definition: gblfit.py:1172
def fit(self, optionList="", aLabel=0)
Perform fit of trajectory.
Definition: gblfit.py:1220
def __getFitToKinkJacobian(self, aPoint)
Get jacobian for transformation from (trajectory) fit to kink parameters at point.
Definition: gblfit.py:1029
def getResults(self, aLabel)
Get results (corrections, covarinace matrix) at point in forward or backward direction.
Definition: gblfit.py:1154
def __getFitToStepJacobian(self, aPoint)
Get jacobian for transformation from (trajectory) fit to step parameters at point.
Definition: gblfit.py:1068
__scatDataIndex
mapping points to data blocks from scatterers; list(int)
Definition: gblfit.py:793
__matrix
Calculate Jacobians to previous/next scatterer from point to point ones.
Definition: gblfit.py:1362
__numParameters
number fit parameters; int
Definition: gblfit.py:777
__data
data (blocks) of trajectory; list(GblData)
Definition: gblfit.py:787
__points
points on trajectory; list(GblPoint)
Definition: gblfit.py:785
def addExternalSeed(self, aLabel, aSeed)
Add external seed to trajectory.
Definition: gblfit.py:823
__skippedMeasLabel
label of point with measurements skipped in fit (for unbiased residuals)
Definition: gblfit.py:795
def milleOut(self, aFile, doublePrec=False)
Write (data blocks of) trajectory to MP (binary) file (as float or double values).
Definition: gblfit.py:851
__numPoints
number of points on trajectory; int
Definition: gblfit.py:771
__externalPoint
label of point with external seed; int
Definition: gblfit.py:781
def __init__(self, hasCurv=True, aDim=[0, 1])
Create new trajectory.
Definition: gblfit.py:769
def printData(self)
Print data of trajectory.
Definition: gblfit.py:834
def getScatResults(self, aLabel)
Get residuals at point from scatterer.
Definition: gblfit.py:1193
def getData(self)
Get data of trajectory.
Definition: gblfit.py:843
def __getResAndErr(self, aData, used=True)
Definition: gblfit.py:1132
def printPoints(self)
Print points of trajectory.
Definition: gblfit.py:828
def addPoint(self, point)
Add point to trajectory.
Definition: gblfit.py:802
def __getFitToKinkAndStepJacobian(self, aPoint)
Get jacobian for transformation from (trajectory) fit to kink and step parameters at point.
Definition: gblfit.py:1099
(Symmetric) Bordered Band Matrix.
Definition: gblnum.py:65
Millepede-II (binary) record.
Definition: mille.py:76