gblfit.py

Go to the documentation of this file.

'''
Track fit with general broken lines.
 
Created on Jul 27, 2011
 
@author: kleinwrt
'''
 
 
 
import numpy as np
import math
from gblnum import BorderedBandMatrix
from mille import MilleRecord
 
 
 
class GblMeasurement(object):
 
  
  def __init__(self, aMeasurement, minPrecision=0.):
    
    self.__measurement = aMeasurement
    
    self.__measDim = aMeasurement[1].shape[0]
    
    self.__measMinPrec = minPrecision
    
    self.__measTransformation = None    
    
    self.__localDerivatives = None
    
    self.__globalLabels = None
    
    self.__globalDerivatives = None
       
    # full precision matrix?
    if (aMeasurement[2].ndim == 2):
      # need to diagonalize
      eigenVal, eigenVec = np.linalg.eigh(aMeasurement[2])
      self.__measTransformation = eigenVec.T
      # transform measurement
      if (aMeasurement[0] is None):
        self.__measurement[0] = self.__measTransformation
      else:
        self.__measurement[0] = np.dot(self.__measTransformation, aMeasurement[0])
      self.__measurement[1] = np.dot(self.__measTransformation, aMeasurement[1])
      self.__measurement[2] = eigenVal
 
  
  def hasMeasurement(self):
 
    return (self.__measurement is not None)
 
  
  def getMeasurement(self):
    return self.__measurement
 
  
  def getMeasDim(self):
    return self.__measDim
  
  
  def getMeasMinPrec(self):
    return self.__measMinPrec
 
  
  def addLocals(self, derivatives):
    if (self.__measDim > 0):
      if (self.__measTransformation is None):
        self.__localDerivatives = derivatives
      else:
        self.__localDerivatives = np.dot(self.__measTransformation, derivatives)        
 
  
  def addGlobals(self, labels, derivatives):
    if (self.__measDim > 0):
      self.__globalLabels = labels
      if (self.__measTransformation is None):
        self.__globalDerivatives = derivatives
      else:
        self.__globalDerivatives = np.dot(self.__measTransformation, derivatives) 
 
  
  def getNumLocals(self):
    if (self.__localDerivatives is not None):
      return self.__localDerivatives.shape[1]
    else:
      return 0
 
  
  def getLocalDerivatives(self):
    return self.__localDerivatives
 
  
  def getGlobalLabels(self):
    return self.__globalLabels
 
  
  def getGlobalDerivatives(self):
    return self.__globalDerivatives
 
  
  def printMeasurement(self):
    print "  measurement ", self.__measDim, len(self.__localDerivatives), len(self.__globalDerivatives)
 
#------------------------------------------------------------------------------ 
 
 
 
class GblPoint(object):
 
  
  def __init__(self, aJacobian):
    
    self.__label = 0
    
    self.__offset = 0
    
    self.__type = 0
    
    self.__p2pJacobian = aJacobian
    
    self.__jacobians = [ [], [] ]
    
    self.__measurements = []
    
    self.__scatDim = 0
    
    self.__scatterer = None
      
# for extension to retrieval of residuals, pulls    
#    self.__dataMeas = [0, 0]
# for extension to retrieval of residuals, pulls    
#    self.__dataScat = [0, 0]
 
  
  def addMeasurement(self, aMeasurement, minPrecision=0.):
    self.__measurements.append(GblMeasurement(aMeasurement, minPrecision))
    
  
  def hasMeasurement(self):
    return (self.__measurements <> [])
 
  
  def getMeasurements(self):
    return self.__measurements
 
  
  def addScatterer(self, aScatterer):
    self.__scatDim = aScatterer[0].shape[0]  # 2 (thin) or 4 (thick)
    self.__scatterer = [ None ] + aScatterer
    if (aScatterer[1].ndim == 2):  # full precision matrix, need to diagonalize
      eigenVal, eigenVec = np.linalg.eigh(aScatterer[1])
      scatTransformation = eigenVec.T
#     transform measurement
      self.__scatterer[0] = scatTransformation
      self.__scatterer[1] = np.dot(scatTransformation, aScatterer[0])
      self.__scatterer[2] = eigenVal
 
  
  def getScatDim(self):
    return self.__scatDim
 
  
  def getScatterer(self):
    return self.__scatterer
 
  
  def getReducedScatterer(self):
    # get scatterer input back from diagonalization
    scatRes = np.dot(self.__scatterer[0].T, self.__scatterer[1])
    # positional covariance from diagonalization
    posCov = np.dot(np.dot(self.__scatterer[0][:, 2:].T, np.diag(1. / self.__scatterer[2])), self.__scatterer[0][:, 2:])
    
    scatPrec = np.zeros((4, 4))
    scatPrec[2:, 2:] = np.linalg.inv(posCov)
    # diagonalize again
    eigenVal, eigenVec = np.linalg.eigh(scatPrec)
    scatTransformation = eigenVec.T
    return [scatTransformation, np.dot(scatTransformation, scatRes), eigenVal]
 
  
  def addLocals(self, derivatives):
    if self.__measurements <> []:
      self.__measurements[-1].addLocals(derivatives)
 
  
  def addGlobals(self, labels, derivatives):
    if self.__measurements <> []:
      self.__measurements[-1].addGlobals(labels, derivatives)
 
  
  def setLabel(self, aLabel):
    self.__label = aLabel
 
  
  def getLabel(self):
    return self.__label
 
  
  def setOffset(self, anOffset):
    self.__offset = anOffset
 
  
  def getOffset(self):
    return self.__offset
  
  
  def setType(self, aType):
    self.__type = aType
 
  
  def isFirst(self):
    return (self.__type < 0)
 
  
  def isLast(self):
    return (self.__type > 0)
 
#  def setDataMeas(self, aIndex, aData):
# for extension to retrieval of residuals, pulls    
#    self.__dataMeas[aIndex] = aData
 
#  def setDataScat(self, aIndex, aData):
# for extension to retrieval of residuals, pulls    
#    self.__dataScat[aIndex] = aData
 
  
  def getP2pJacobian(self):
    return self.__p2pJacobian         
 
  
  def addPrevJacobian(self, aJacobian):
    self.__jacobians[0] = np.linalg.inv(aJacobian)
 
  
  def addNextJacobian(self, aJacobian):
    self.__jacobians[1] = aJacobian
 
  
  def getDerivatives(self, index):
    aJacobian = self.__jacobians[index]
    matJ = aJacobian[3:5, 3:5]  # J
    matS = aJacobian[3:5, 1:3]  # S 
    vecD = aJacobian[3:5, 0:1]  # d
    if (index < 1):
      matS = -matS
    matW = np.linalg.inv(matS)  # W = +/- S^-1
    return matW, np.dot(matW, matJ), np.dot(matW, vecD)  # W, W*J, W*d
 
  
  def printPoint(self):
    print " point ", self.__label, self.__offset, self.__scatDim, len(self.__measurements)
 
#------------------------------------------------------------------------------ 
 
 
 
class GblData(object):
 
  
  def __init__(self, aLabel=0, aType=0, aValue=0., aPrec=0., aMeas=0):
    
    self.__label = aLabel
    
    self.__type = aType
    
    self.__index = aMeas
    
    self.__value = aValue
    
    self.__precision = aPrec
    
    self.__dwMethod = 0
    
    self.__downWeight = 1.
    
    self.__prediction = 0.
    
    self.__parameters = []
    
    self.__derivatives = []
    
    self.__globalLabels = []
    
    self.__globalDerivatives = []
 
#    self.__predictionVar = 0.
  
  
  def addDerivatives(self, iRow, labDer, matDer, \
                     derLocal=None, labGlobal=None, derGlobal=None):
    if (derLocal is not None):
      for i in range(derLocal.shape[1]):  # local derivatives
        if (derLocal[iRow, i] != 0.):
          self.__derivatives.append(derLocal[iRow, i])
          self.__parameters.append(i + 1)
        
    for i in range(len(labDer)):  # curvature, offset derivatives
      if (labDer[i] != 0 and matDer[iRow, i] != 0.):
        self.__derivatives.append(matDer[iRow , i])
        self.__parameters.append(labDer[i])
 
    if (derGlobal is not None): 
      for i in range(derGlobal.shape[1]):  # global derivatives
        if (derGlobal[iRow, i] != 0.):
          self.__globalLabels.append(labGlobal[iRow, i])
          self.__globalDerivatives.append(derGlobal[iRow, i])
 
  
  def addExtDerivatives(self, indexExt, derExt):
    for i in range(len(derExt)):  # external derivatives
      if (derExt[i] != 0.):
        self.__derivatives.append(derExt[i])
        self.__parameters.append(indexExt[i])  
 
  
  def getMatrices(self): 
    aVector = np.array([ self.__derivatives ])
    aMatrix = np.dot(aVector.T, aVector)
    aValue = self.__value
    aWeight = self.__precision * self.__downWeight
    return self.__parameters, aValue * aVector * aWeight, aMatrix * aWeight
 
  
  def setPrediction(self, aVector):
    self.__prediction = 0.
    for i in range(len(self.__parameters)):
      self.__prediction += self.__derivatives[i] * aVector[ self.__parameters[i] - 1 ]
 
#  def setPredictionVariance(self, aMatrix):
#    '''Calculate variance of prediction for data from fit.'''
# var(residual) = 1./precision**2 - var(prediction)
#    aBlockMatrix = aMatrix.getBlockMatrix(self.__parameters)
#    self.__predictionVar = np.dot(self.__derivatives.T, \
#                             np.dot(aBlockMatrix, self.__derivatives))
 
  
  def getGlobalLabels(self):
    return self.__globalLabels
 
  
  def setDownWeighting(self, aMethod): 
    self.__dwMethod = aMethod
    scaledResidual = abs(self.__value - self.__prediction) * math.sqrt(self.__precision)   
    if (aMethod == 1):  # Tukey
      if (scaledResidual < 4.6851):
        aWeight = (1.0 - (scaledResidual / 4.6851) ** 2) ** 2
      else:
        aWeight = 0.
    elif (aMethod == 2):  # Huber
      if (scaledResidual < 1.345):
        aWeight = 1.
      else:
        aWeight = 1.345 / scaledResidual
    elif (aMethod == 3):  # Cauchy
      aWeight = 1.0 / (1.0 + (scaledResidual / 2.3849) ** 2)      
    self.__downWeight = aWeight
    return aWeight
  
  
  def getChi2(self):
    scaledResidual = abs(self.__value - self.__prediction) * math.sqrt(self.__precision)
    Chi2 = scaledResidual ** 2
    if (self.__dwMethod == 1):  # Tukey
      if (scaledResidual < 4.6851):
        Chi2 = 4.6851 ** 2 / 3. * (1. - (1. - (scaledResidual / 4.6851) ** 2) ** 3)
      else: 
        Chi2 = 4.6851 ** 2 / 3.
    elif (self.__dwMethod == 2):  # Huber
      if (scaledResidual > 1.345):
        Chi2 = 1.345 * (2.*scaledResidual - 1.345) 
    elif (self.__dwMethod == 3):  # Cauchy
      Chi2 = math.log(1. + (scaledResidual / 2.3849) ** 2) * 2.3849 ** 2     
    return Chi2
  
  
  def getLabel(self):
    return self.__label
 
  
  def getType(self):
    return self.__type
 
  
  def getIndex(self):
    return self.__index
      
  
  def getResidual(self):
    return self.__value - self.__prediction, 1.0 / self.__precision, self.__downWeight, self.__parameters, self.__derivatives
 
  
  def toRecord(self):
    return self.__value, self.__precision, self.__parameters, self.__derivatives, \
            self.__globalLabels, self.__globalDerivatives
 
  
  def fromRecord(self, dataList):
    self.__value, self.__precision, self.__parameters, self.__derivatives, \
            self.__globalLabels, self.__globalDerivatives = dataList
 
  
  def analyzeData(self, maxBand):
    maxPar = self.__parameters[-1]
    maxBor = 0
    for i in self.__parameters:
      if (i < maxPar - maxBand):
        maxBor = i
    return maxPar, maxBor
 
  
  def printData(self):
    print " measurement at label ", self.__label, " with type ", self.__type, " : ", self.__value, self.__precision
    print " param ", self.__parameters
    print " deriv ", self.__derivatives
    print " global labels ", self.__globalLabels
    print " global deriv  ", self.__globalDerivatives
 
#------------------------------------------------------------------------------ 
 
 
 
 
 
class GblTrajectory(object):
 
  
  def __init__(self, hasCurv=True, aDim=[0, 1]):
    
    self.__numPoints = 0
    
    self.__numOffsets = 0
    
    self.__numCurvature = (1 if hasCurv else 0)
    
    self.__numParameters = 0
    
    self.__numLocals = 0
    
    self.__externalPoint = 0
    
    self.__dimensions = aDim
    
    self.__points = [] 
    
    self.__data = []
    
    self.__externalSeed = None
    
    self.__measDataIndex = []
    
    self.__scatDataIndex = []
    
    self.__skippedMeasLabel = 0
    
  
  def addPoint(self, point):
    self.__numPoints += 1
    label = self.__numPoints
    point.setLabel(label)
    self.__points.append(point)
    for m in point.getMeasurements():
      self.__numLocals = max(self.__numLocals, m.getNumLocals())
    return label
 
  
  def getNumPoints(self):
    return self.__numPoints
  
  
  def getNumbers(self):
    return self.__numLocals, self.__numCurvature, self.__numParameters
 
  
  def addExternalSeed(self, aLabel, aSeed):
    self.__externalPoint = aLabel
    self.__externalSeed = aSeed
 
  
  def printPoints(self):
    print "GblPoints"
    for p in self.__points:
      p.printPoint() 
 
  
  def printData(self):
    print "GblData blocks"
    for d in self.__data:
      d.printData() 
 
  
  def getData(self):
    return self.__data
 
  
  def milleOut(self, aFile, doublePrec=False):
    rec = MilleRecord(doublePrec)
#   data measurements and kinks        
    for aData in self.__data: 
      rec.addData(aData.toRecord())
                    
    rec.writeRecord(aFile)
 
  
  def milleIn(self, aFile):
    rec = MilleRecord()
    rec.readRecord(aFile)
    mPar = 0
    mBor = 0
    mBand = 4 * len(self.__dimensions) - 1  # max band width
    while (rec.moreData()):
      aTag = rec.specialDataTag()
      if (aTag < 0):
# get data
        aData = GblData()
        aData.fromRecord(rec.getData())
        self.__data.append(aData)
        nPar, nBor = aData.analyzeData(mBand)
        mPar = max(mPar, nPar)
        mBor = max(mBor, nBor)
        
    self.__numParameters = mPar
    self.__numLocals = mBor - self.__numCurvature
 
  
  def __getJacobian(self, aLabel):
    aDim = self.__dimensions
    nDim = len(aDim)
    anIndex = abs(aLabel) - 1
#   check consistency of (index, direction)    
    if (aLabel > 0):
      nJacobian = 1
      if (anIndex >= self.__numPoints - 1):
        anIndex = self.__numPoints - 1
        nJacobian = 0
    else:
      nJacobian = 0
      if (anIndex <= 0):
        anIndex = 0
        nJacobian = 1
# Jacobian broken lines (q/p,..,u_i,u_i+1..) to local (q/p,u',u) parameters   
    nCurv = self.__numCurvature
    nLocals = self.__numLocals   
    nBorder = nCurv + nLocals
    nParBRL = nBorder + 2 * nDim
    nParLoc = nLocals + 5
    aJacobian = np.zeros((nParLoc, nParBRL))
    aPoint = self.__points[anIndex]
    anIndex = []
    labDer, matDer = self.__getFitToLocalJacobian(aPoint, 5, nJacobian)
#   from local parameters
    for i in range(nLocals):
      aJacobian[i + 5, i] = 1.0;
      anIndex.append(i + 1);
  
#   from trajectory parameters
    iCol = nLocals;
    for i in range(5):
      if (labDer[i] > 0):
        anIndex.append(labDer[i]);
        for j in range(5):
          aJacobian[j, iCol] = matDer[j, i];
        iCol += 1
 
    return anIndex, aJacobian
 
  
  def __getFitToLocalJacobian(self, aPoint, measDim, nJacobian=1):
    aDim = self.__dimensions
    nDim = len(aDim)
    nCurv = self.__numCurvature
    nLocals = self.__numLocals      
    nOffset = aPoint.getOffset()
    scatDim = aPoint.getScatDim()
    anIndex = [0, 0, 0, 0, 0]
    aJacobian = np.zeros((measDim, 5))
    labOffset = measDim - 2
    labSlope = measDim - 4
    labCurv = measDim - 5
    
    if (nOffset < 0):  # need interpolation
      prevW, prevWJ, prevWd = aPoint.getDerivatives(0)  # W-, W- * J-, W- * d-
      nextW, nextWJ, nextWd = aPoint.getDerivatives(1)  # W+, W+ * J+, W+ * d+
      sumWJ = prevWJ + nextWJ
      matN = np.linalg.inv(sumWJ)  # N = (W- * J- + W+ * J+)^-1 
#     local offset
      if (labOffset >= 0):
#       derivatives for u_int      
        prevNW = np.dot(matN, prevW)  # N * W-
        nextNW = np.dot(matN, nextW)  # N * W+
        prevNd = np.dot(matN, prevWd)  # N * W- * d-
        nextNd = np.dot(matN, nextWd)  # N * W+ * d+
        iOff = nDim * (-nOffset - 1) + nLocals + nCurv + 1  # first offset ('i' in u_i)
        if (nCurv > 0):
          aJacobian[labOffset:measDim, 0:1] = -prevNd - nextNd  # from curvature
          anIndex[0] = nLocals + 1
        aJacobian[labOffset:measDim, 1:3] = prevNW
        aJacobian[labOffset:measDim, 3:5] = nextNW
        for i in range(nDim):
          anIndex[1 + aDim[i]] = iOff + i
          anIndex[3 + aDim[i]] = iOff + nDim + i
#     local slope 
      if (labSlope >= 0):
#       derivatives for u'_int
        prevWPN = np.dot(nextWJ, prevNW)  # W+ * J+ * N * W-
        nextWPN = np.dot(prevWJ, nextNW)  # W- * J- * N * W+
        prevWNd = np.dot(nextWJ, prevNd)  # W+ * J+ * N * W- * d-
        nextWNd = np.dot(prevWJ, nextNd)  # W- * J- * N * W+ * d+
        if (nCurv > 0):
          aJacobian[labSlope:labOffset, 0:1] = prevWNd - nextWNd  # from curvature
        aJacobian[labSlope:labOffset, 1:3] = -prevWPN
        aJacobian[labSlope:labOffset, 3:5] = nextWPN
 
    else:  # at point               
#     anIndex must be sorted
#     forward : iOff2 = iOff1 + nDim, index1 = 1, index2 = 3
#     backward: iOff2 = iOff1 - nDim, index1 = 3, index2 = 1
      # adjust for thick scatterer
      if (scatDim == 4):
        nOffset += 1
      iOff1 = nDim * nOffset + nCurv + nLocals + 1  # first offset ('i' in u_i)
      index1 = 3 - 2 * nJacobian  # index of first offset
      iOff2 = iOff1 + nDim * (nJacobian * 2 - 1)  # second offset ('i' in u_i)
      index2 = 1 + 2 * nJacobian  # index of second offset
#     local offset
      if (labOffset >= 0):
        aJacobian[labOffset, index1] = 1.0  # from 1st Offset
        aJacobian[labOffset + 1, index1 + 1] = 1.0
        for i in range(nDim):
          anIndex[index1 + aDim[i]] = iOff1 + i
#     local slope and curvature
      if (labSlope >= 0):
        matW, matWJ, vecWd = aPoint.getDerivatives(nJacobian)  # W, W * J, W * d
        sign = 2 * nJacobian - 1
        if (nCurv > 0):
          aJacobian[labSlope:labOffset, 0:1] = -sign * vecWd  # from curvature
          anIndex[0] = nLocals + 1
        aJacobian[labSlope:labOffset, index1:index1 + 2] = -sign * matWJ
        aJacobian[labSlope:labOffset, index2:index2 + 2] = sign * matW
        for i in range(nDim):
          anIndex[index2 + aDim[i]] = iOff2 + i  
          
#   local curvature
    if (labCurv >= 0):
      if (nCurv > 0):
        aJacobian[labCurv, labCurv] = 1.0  
                    
    return anIndex, aJacobian    
 
  
  def __getFitToKinkJacobian(self, aPoint): 
    aDim = self.__dimensions
    nDim = len(aDim)
    nCurv = self.__numCurvature
    nLocals = self.__numLocals        
    nOffset = aPoint.getOffset() 
    anIndex = [0, 0, 0, 0, 0, 0, 0]
    aJacobian = np.zeros((2, 7))    
 
    prevW, prevWJ, prevWd = aPoint.getDerivatives(0)  # W-, W- * J-, W- * d-
    nextW, nextWJ, nextWd = aPoint.getDerivatives(1)  # W+, W+ * J+, W+ * d+
    sumWJ = prevWJ + nextWJ  # W- * J- + W+ * J+
    sumWd = prevWd + nextWd  # W+ * d+ + W- * d-
    iOff = (nOffset - 1) * nDim + nCurv + nLocals + 1  # first offset ('i' in u_i)
 
#   local offset
    if (nCurv > 0):
      aJacobian[:, 0:1] = -sumWd  # from curvature
      anIndex[0] = nLocals + 1      
    aJacobian[:, 1:3] = prevW  # from 1st Offset
    aJacobian[:, 3:5] = -sumWJ  # from 2nd Offset
    aJacobian[:, 5:7] = nextW  # from 3rd Offset
    for i in range(nDim):
      anIndex[1 + aDim[i]] = iOff + i
      anIndex[3 + aDim[i]] = iOff + nDim + i
      anIndex[5 + aDim[i]] = iOff + nDim * 2 + i
        
    return anIndex, aJacobian    
 
  
  def __getFitToStepJacobian(self, aPoint):
    aDim = self.__dimensions
    nDim = len(aDim)
    nCurv = self.__numCurvature
    nLocals = self.__numLocals
    nOffset = aPoint.getOffset()
    anIndex = [0, 0, 0, 0]
    aJacobian = np.zeros((4, 4))
 
    iOff = (nOffset - 1) * nDim + nCurv + nLocals + 1  # first offset ('i' in u_i)
 
#   step
    aJacobian[2:,:2] = -np.eye(2)  # from 2nd Offset
    aJacobian[2:, 2:] = np.eye(2)  # from 3rd Offset
 
    for i in range(nDim):
      anIndex[aDim[i]] = iOff + nDim + i
      anIndex[2 + aDim[i]] = iOff + nDim * 2 + i
 
    return anIndex, aJacobian
 
  
  def __getFitToKinkAndStepJacobian(self, aPoint):
    aDim = self.__dimensions
    nDim = len(aDim)
    nCurv = self.__numCurvature
    nLocals = self.__numLocals
    nOffset = aPoint.getOffset()
    anIndex = [0, 0, 0, 0, 0, 0, 0, 0, 0]
    aJacobian = np.zeros((4, 9))
 
    prevW, prevWJ, prevWd = aPoint.getDerivatives(0)  # W-, W- * J-, W- * d-
    nextW, nextWJ, nextWd = aPoint.getDerivatives(1)  # W+, W+ * J+, W+ * d+
    iOff = (nOffset - 1) * nDim + nCurv + nLocals + 1  # first offset ('i' in u_i)
 
#   kink
    if (nCurv > 0):
      aJacobian[:2, 0:1] = -(prevWd + nextWd)  # from curvature
      anIndex[0] = nLocals + 1
    aJacobian[:2, 1:3] = prevW  # from 1st Offset
    aJacobian[:2, 3:5] = -prevWJ  # from 2nd Offset
    aJacobian[:2, 5:7] = -nextWJ  # from 3rd Offset
    aJacobian[:2, 7:9] = nextW  # from 4th Offset
#   step
    aJacobian[2:, 3:5] = -np.eye(2)  # from 2nd Offset
    aJacobian[2:, 5:7] = np.eye(2)  # from 3rd Offset
 
    for i in range(nDim):
      anIndex[1 + aDim[i]] = iOff + i
      anIndex[3 + aDim[i]] = iOff + nDim + i
      anIndex[5 + aDim[i]] = iOff + nDim * 2 + i
      anIndex[7 + aDim[i]] = iOff + nDim * 3 + i
 
    return anIndex, aJacobian
 
  def __getResAndErr(self, aData, used=True):
    aResidual, aMeasVar, aDownWeight, indLocal, derLocal = self.__data[aData].getResidual()
    aVec = np.array(derLocal)  # compressed vector
    aMat = self.__matrix.getBlockMatrix(indLocal)  # compressed matrix     
    aFitVar = np.dot(aVec, np.dot(aMat, aVec.T))  # variance from track fit
    aFitVar *= aDownWeight  # account for down-weighting (of measurement in fit)
    aMeasError = math.sqrt(aMeasVar)  # error of measurement
    if used:
      aResError = math.sqrt(aMeasVar - aFitVar) if aFitVar < aMeasVar else 0.  # error of biased residual
    else:
      aResError = math.sqrt(aMeasVar + aFitVar)  # error of unbiased residual
    return aResidual, aMeasError, aResError, aDownWeight 
  
  
  def getResults(self, aLabel):
    anIndex, aJacobian = self.__getJacobian(aLabel)
    nParBRL = len(anIndex)
    aVec = np.empty(nParBRL)
    for i in range(nParBRL):
      aVec[i] = self.__vector[anIndex[i] - 1]  # compressed vector
    aMat = self.__matrix.getBlockMatrix(anIndex)  # compressed matrix
    locPar = np.dot(aJacobian, aVec) 
    locCov = np.dot(aJacobian, np.dot(aMat, aJacobian.T))
    return  locPar, locCov 
  
  
  def getMeasResults(self, aLabel):
    firstData = self.__measDataIndex[aLabel - 1]  # first data block with measurement
    numData = self.__measDataIndex[aLabel] - firstData  # number of data blocks
    if numData <= 0:
      return 0, None, None, None, None
    #
    aResiduals = np.empty(numData)
    aMeasErr = np.empty(numData)
    aResErr = np.empty(numData)
    aDownWeight = np.empty(numData)
    for i in range(numData):
      aResiduals[i], aMeasErr[i], aResErr[i], aDownWeight[i] = self.__getResAndErr(firstData + i, (aLabel <> self.__skippedMeasLabel))
    return numData, aResiduals, aMeasErr, aResErr, aDownWeight
 
  
  def getScatResults(self, aLabel):
    firstData = self.__scatDataIndex[aLabel - 1]  # first data block with measurement
    numData = self.__scatDataIndex[aLabel] - firstData  # number of data blocks
    if numData <= 0:
      return 0, None, None, None, None
    #
    aResiduals = np.empty(numData)
    aMeasErr = np.empty(numData)
    aResErr = np.empty(numData)
    aDownWeight = np.empty(numData)
    for i in range(numData):
      aResiduals[i], aMeasErr[i], aResErr[i], aDownWeight[i] = self.__getResAndErr(firstData + i)
    return numData, aResiduals, aMeasErr, aResErr, aDownWeight
 
  
  def getResidual(self, iData):
    return self.__getResAndErr(iData)   
         
  
  def fit(self, optionList="", aLabel=0):
 
    
    def defineOffsets():
# set labels for previous/next offsets
#     first point is offset    
      self.__points[0].setOffset(0)
      self.__points[0].setType(-1)
      nOffsets = 1
#     intermediate scatterers are offsets    
      for aPoint in self.__points[1:-1]:
        scatDim = aPoint.getScatDim()
        if (scatDim > 0):
          aPoint.setOffset(nOffsets)
          nOffsets += 1
          # thick scatterer ?
          if (scatDim == 4):
            nOffsets += 1
        else: 
          aPoint.setOffset(-nOffsets)
#     last point is offset    
      self.__points[-1].setOffset(nOffsets)
      self.__points[-1].setType(1)
      if (self.__points[-1].getScatDim() == 4):
        # thick scatterer
        nOffsets += 1
      self.__numOffsets = nOffsets + 1
      self.__numParameters = self.__numOffsets * len(self.__dimensions) \
                           +self.__numCurvature + self.__numLocals       
 
    
    def calcJacobians():
      scatJacobian = np.empty((5, 5)) 
#     forward propagation (all)
      lastPoint = 0;
      numStep = 0;
      for iPoint in range(1, self.__numPoints):
        if (numStep == 0):
          scatJacobian = self.__points[iPoint].getP2pJacobian()
        else: 
          scatJacobian = np.dot(self.__points[iPoint].getP2pJacobian(), scatJacobian)
        numStep += 1
        self.__points[iPoint].addPrevJacobian(scatJacobian)  # aPoint -> previous scatterer
        if (self.__points[iPoint].getOffset() >= 0):
          self.__points[lastPoint].addNextJacobian(scatJacobian)  # lastPoint -> next scatterer
          numStep = 0;
          lastPoint = iPoint
#     backward propagation (without scatterers)
      for iPoint in range(self.__numPoints - 1, 0, -1):
        if (self.__points[iPoint].getOffset() >= 0):
          scatJacobian = self.__points[iPoint].getP2pJacobian()
          continue  # skip offsets
        self.__points[iPoint].addNextJacobian(scatJacobian)  # iPoint -> next scatterer
        scatJacobian = np.dot(scatJacobian, self.__points[iPoint].getP2pJacobian())
 
    
    def prepare():
      aDim = self.__dimensions
# measurements
      self.__measDataIndex.append(len(self.__data))  # offset
      for aPoint in self.__points:
        if (aPoint.hasMeasurement()):
          nLabel = aPoint.getLabel()
          measIndex = 0
          for m in aPoint.getMeasurements():
            measIndex += 1
            measDim = m.getMeasDim()
            measPrecision = m.getMeasMinPrec()
            localDer = m.getLocalDerivatives()
            globalLab = m.getGlobalLabels()
            globalDer = m.getGlobalDerivatives()
            matP, aMeas, aPrec = m.getMeasurement()
            nJacobian = 1 if aPoint.getOffset() < self.__numOffsets - 1 else 0  # last point needs backward propagation
            labDer, matDer = self.__getFitToLocalJacobian(aPoint, measDim, nJacobian)
            matPDer = matDer if matP is None else np.dot(matP, matDer)
            for i in range(measDim):
              if (aPrec[i] > measPrecision):
                aData = GblData(nLabel, 1, aMeas[i], aPrec[i], measIndex)
                aData.addDerivatives(i, labDer, matPDer, localDer, \
                                     globalLab, globalDer)
                self.__data.append(aData)
        self.__measDataIndex.append(len(self.__data))
#                aPoint.setDataMeas(i, len(self.__data)) 
# pseudo measurements from kinks
      self.__scatDataIndex.append(len(self.__data))  # offset
      self.__scatDataIndex.append(len(self.__data))  # first point
      for aPoint in self.__points[1:]:
        scatDim = aPoint.getScatDim()
        if (scatDim > 0):
          nLabel = aPoint.getLabel()        
          if (scatDim == 4):
          # thick scatterer, last point?
            if (aPoint.isLast()):
              # steps (only)
              matT, aMeas, aPrec = aPoint.getReducedScatterer()
              labDer, matDer = self.__getFitToStepJacobian(aPoint)
            else:
              # kinks+steps
              matT, aMeas, aPrec = aPoint.getScatterer()
              labDer, matDer = self.__getFitToKinkAndStepJacobian(aPoint)
          else:
          # thin scatterer, last point?
            if (aPoint.isLast()):
              continue
            else:
              # kinks
              matT, aMeas, aPrec = aPoint.getScatterer()
              labDer, matDer = self.__getFitToKinkJacobian(aPoint)
          matTDer = matDer if matT is None else np.dot(matT, matDer)
          for i in aDim:
            if (aPrec[i] > 0.):
              aData = GblData(nLabel, 2, aMeas[i], aPrec[i])
              aData.addDerivatives(i, labDer, matTDer)
              self.__data.append(aData)
          # thick scatterer?
          if (scatDim == 4):
            for i in aDim:
              if (aPrec[i + 2] > 0.):
                aData = GblData(nLabel, 2, aMeas[i + 2], aPrec[i + 2])
                aData.addDerivatives(i + 2, labDer, matTDer)
                self.__data.append(aData)
        self.__scatDataIndex.append(len(self.__data))
#              aPoint.setDataScat(i, len(self.__data)) 
      self.__scatDataIndex.append(len(self.__data))  # last point
#     external seed
      if (self.__externalPoint != 0):
        externalIndex, aJacobian = self.__getJacobian(self.__externalPoint)
        eigenVal, eigenVec = np.linalg.eigh(self.__externalSeed)
        aMatrix = np.dot(eigenVec.T, aJacobian)
        for i in range(len(eigenVec)):
          if (eigenVal[i] > 0.):
            externalDerivatives = []
            for j in range(len(externalIndex)):
              externalDerivatives.append(aMatrix[i, j])
            aData = GblData(self.__externalPoint, 3, 0., eigenVal[i])
            aData.addExtDerivatives(externalIndex, externalDerivatives)
            self.__data.append(aData)
        self.__measDataIndex.append(len(self.__data))
 
    
    def buildLinearEquationSystem():
      nBorder = self.__numCurvature + self.__numLocals
      self.__matrix = BorderedBandMatrix(self.__numParameters, nBorder)
      self.__vector = np.zeros(self.__numParameters)
      for aData in self.__data:
        # skipped (internal) measurement?
        if aData.getLabel() == self.__skippedMeasLabel and aData.getType() == 1:
          continue 
        index, aVector, aMatrix = aData.getMatrices()
        for i in range(len(index)):
          self.__vector[ index[i] - 1 ] += aVector[0, i]  # update vector
        self.__matrix.addBlockMatrix(index, aMatrix)  # update matrix
 
    
    def downWeight(aMethod):
      aLoss = 0.
      for aData in self.__data: 
        aLoss += (1. - aData.setDownWeighting(aMethod))
      return aLoss
 
    
    def predict():
      for aData in self.__data: 
        aData.setPrediction(self.__vector)
    
    if (self.__data == []):  # generate data from points   
      defineOffsets()    
      calcJacobians()                      
      prepare()
 
    # skip measurements from point?
    self.__skippedMeasLabel = aLabel
 
    buildLinearEquationSystem()  # create linear equations system from data
    #
    try:
      aMethod = 0
      lostWeight = 0.
      self.__vector = self.__matrix.solveAndInvertBorderedBand(self.__vector)
      predict()
 
      for o in optionList:  # down weighting iterations    
        try:
          aMethod = "THC".index(o.upper()) + 1
          lostWeight = downWeight(aMethod)
          buildLinearEquationSystem()
          self.__vector = self.__matrix.solveAndInvertBorderedBand(self.__vector)
          predict()
        except ValueError:
          pass                  
             
      Ndf = -self.__numParameters 
      Chi2 = 0.
      for aData in self.__data: 
        # skipped (internal) measurement?
        if aData.getLabel() == self.__skippedMeasLabel and aData.getType() == 1:
          continue 
        Chi2 += aData.getChi2()
        Ndf += 1
      Chi2 /= [1.0, 0.8737, 0.9326, 0.8228 ][aMethod]  
      return Chi2, Ndf, lostWeight #, self.__matrix.getBandCond()
    
    except (ZeroDivisionError, np.linalg.linalg.LinAlgError):
      return  0., -1, 0. #, 0.