Commit 31e726d7 authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

Python3 version (gblpy3) added

parent 9c2c5893
README ident
gbltst.py ident
gblsit.py ident
gbltst*.py ident
gblsit*.py ident
example*.cpp ident
LAST_COMMIT export-subst
......
......@@ -7,7 +7,7 @@ PROJECT(GBL)
# project version
SET( ${PROJECT_NAME}_VERSION_MAJOR 2 )
SET( ${PROJECT_NAME}_VERSION_MINOR 4 )
SET( ${PROJECT_NAME}_VERSION_PATCH 0 )
SET( ${PROJECT_NAME}_VERSION_PATCH 1 )
# make life easier and simply use the ilcsoft default settings
# load default ilcsoft settings (install prefix, build type, rpath, etc.)
......
$Rev$
$Date$
$Id$
C. Kleinwort, DESY
===========================================================
Helmholtz Terascale Alliance
......@@ -15,5 +15,4 @@ V. Blobel, NIM A, 566 (2006), pp. 14-17
===========================================================
Implementation in Python. Requires at least Python version
2.5 and the 'numpy' module. Documentation generated by
epydoc is available in doc (see doc/index.html) subdirectory.
2.5 and the 'numpy' module.
This diff is collapsed.
$Id$
C. Kleinwort, DESY
===========================================================
Helmholtz Terascale Alliance
General Broken Lines Analysis Centre
Statistics Tools Group
===========================================================
Track refitting with broken lines in 3D building on:
A new fast track-fit algorithm based on broken lines
V. Blobel, NIM A, 566 (2006), pp. 14-17
===========================================================
Implementation in Python3. Requires the 'numpy' module.
This diff is collapsed.
'''
Algebra for linear equation system with bordered band matrix.
Created on Jul 27, 2011
@author: kleinwrt
'''
## \file
# Bordered band matrix.
#
# \author Claus Kleinwort, DESY, 2021 (Claus.Kleinwort@desy.de)
#
# \copyright
# Copyright (c) 2011 - 2016 Deutsches Elektronen-Synchroton,
# Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY \n\n
# This library is free software; you can redistribute it and/or modify
# it under the terms of the GNU Library General Public License as
# published by the Free Software Foundation; either version 2 of the
# License, or (at your option) any later version. \n\n
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Library General Public License for more details. \n\n
# You should have received a copy of the GNU Library General Public
# License along with this program (see the file COPYING.LIB for more
# details); if not, write to the Free Software Foundation, Inc.,
# 675 Mass Ave, Cambridge, MA 02139, USA.
import numpy as np
## (Symmetric) Bordered Band Matrix.
#
# Separate storage of border, mixed and band parts.
#
#\verbatim
# Example for matrix size=8 with border size and band width of two:
#
# +- -+
# | B11 B12 M13 M14 M15 M16 M17 M18 |
# | B12 B22 M23 M24 M25 M26 M27 M28 |
# | M13 M23 C33 C34 C35 0. 0. 0. |
# | M14 M24 C34 C44 C45 C46 0. 0. |
# | M15 M25 C35 C45 C55 C56 C57 0. |
# | M16 M26 0. C46 C56 C66 C67 C68 |
# | M17 M27 0. 0. C57 C67 C77 C78 |
# | M18 M28 0. 0. 0. C68 C78 C88 |
# +- -+
#
# Is stored as:
#
# +- -+ +- -+
# | B11 B12 | | M13 M14 M15 M16 M17 M18 |
# | B12 B22 | | M23 M24 M25 M26 M27 M28 |
# +- -+ +- -+
#
# +- -+
# | C33 C44 C55 C66 C77 C88 |
# | C34 C45 C56 C67 C78 0. |
# | C35 C46 C57 C68 0. 0. |
# +- -+
#\endverbatim
#
class BorderedBandMatrix(object):
## Create new BBmatrix.
#
# @param nSize size of matrix; int
# @param nBorder size of border (default: 1, 'curvature'); int
# @param nBand (maximal) band width (5); int
#
def __init__(self, nSize, nBorder=1, nBand=5):
nSizeBand = nSize - nBorder
## size of matrix; int
self.__numSize = nSize
## size of border; int
self.__numBorder = nBorder
## size of border; int
self.__numBand = 0 # actual band width
## size of band part of matrix; int
self.__numCol = nSizeBand
## border part B; matrix(float)
self.__border = np.zeros((nBorder, nBorder))
## mixed part M; matrix(float)
self.__mixed = np.zeros((nBorder, nSizeBand))
## band part C; matrix(float)
self.__band = np.zeros((nBand + 1, nSizeBand))
# print " new BBM ", self.__border__.shape, self.__mixed__.shape, self.__band__.shape
## Add (compressed) block to BBmatrix:
#
# BBmatrix(aIndex(i),aIndex(j)) += aMatrix(i,j)
#
# @param aIndex list of indices; list(int)
# @param aMatrix (compressed) matrix; matrix(float)
#
def addBlockMatrix(self, aIndex, aMatrix):
nBorder = self.__numBorder
for i in range(len(aIndex)):
iIndex = aIndex[i] - 1
for j in range(i + 1):
jIndex = aIndex[j] - 1
if (iIndex < nBorder):
self.__border[iIndex, jIndex] += aMatrix[i, j]
if (iIndex != jIndex):
self.__border[jIndex, iIndex] += aMatrix[i, j]
elif (jIndex < nBorder):
self.__mixed[jIndex, iIndex - nBorder] += aMatrix[i, j]
else:
nBand = iIndex - jIndex
self.__band[nBand, jIndex - nBorder] += aMatrix[i, j]
self.__numBand = max(self.__numBand, nBand)
return self.__numBand
## Retrieve (compressed) block from BBmatrix:
#
# aMatrix(i,j) = BBmatrix(aIndex(i),aIndex(j))
#
# @param aIndex list of indices; list(int)
# @return (compressed) matrix; matrix(float)
#
def getBlockMatrix(self, aIndex):
nBorder = self.__numBorder
nBand = self.__numBand
nSize = len(aIndex)
aMatrix = np.empty((nSize, nSize))
for i in range(nSize):
iIndex = aIndex[i] - 1
for j in range(i + 1):
jIndex = aIndex[j] - 1
iMax = max(iIndex, jIndex)
iMin = min(iIndex, jIndex)
if (iMax < nBorder):
aMatrix[i, j] = self.__border[iIndex, jIndex]
elif (iMin < nBorder):
aMatrix[i, j] = self.__mixed[iMin, iMax - nBorder]
else:
iBand = iMax - iMin
if iBand <= nBand:
aMatrix[i, j] = self.__band[iBand, iMin - nBorder]
else:
aMatrix[i, j] = 0.
aMatrix[j, i] = aMatrix[i, j]
return aMatrix
## Print BBmatrix.
#
def printMatrix(self):
print(" block part ")
nRow = self.__numBorder
for i in range(nRow):
print(" row ", i, self.__border[i])
print(" mixed part ")
for i in range(nRow):
print(" row ", i, self.__mixed[i])
nRow = self.__numBand + 1
print(" band part ")
for i in range(nRow):
print(" diagonal ", i, self.__band[i])
## Solve linear equation A*x=b system with BBmatrix A, calculate BB part of inverse of A.
#
# @param aRightHandSide right hand side 'b' of linear equation system; vector(float)
# @return solution; vector(float)
# @exception ZeroDivisionError Band matrix is not positive definite
# @note BBmatrix is replaced by BB part of it's inverse
#
def solveAndInvertBorderedBand(self, aRightHandSide):
#============================================================================
## from Dbandmatrix.F (MillePede-II by V. Blobel, Univ. Hamburg)
#============================================================================
##(root free) Cholesky decomposition of band part: C=LDL^T
#
# @note band part (C) is replaced by its decomposition (D,L)
#
def decomposeBand():
nRow = self.__numBand + 1
nCol = self.__numCol
auxVec = np.copy(self.__band[0]) * 16.0 # save diagonal elements
for i in range(nCol):
if ((self.__band[0, i] + auxVec[i]) != self.__band[0, i]):
self.__band[0, i] = 1.0 / self.__band[0, i]
else:
self.__band[0, i] = 0.0 # singular
for j in range(min(nRow, nCol - i) - 1):
rxw = self.__band[j + 1, i] * self.__band[0, i]
for k in range(min(nRow, nCol - i) - j - 1):
self.__band[k, i + j + 1] -= self.__band[k + j + 1, i] * rxw
self.__band[j + 1, i] = rxw
## Solve linear equation system for band part.
#
# @param aRightHandSide right hand side; aRightHandSide vector(float)
# @return solution; vector(float)
#
def solveBand(aRightHandSide):
nRow = self.__numBand + 1
nCol = self.__numCol
aSolution = np.copy(aRightHandSide)
for i in range(nCol): # forward substitution
for j in range(min(nRow, nCol - i) - 1):
aSolution[j + i + 1] -= self.__band[j + 1, i] * aSolution[i]
for i in range(nCol - 1, -1, -1): # backward substitution
rxw = self.__band[0, i] * aSolution[i]
for j in range(min(nRow, nCol - i) - 1):
rxw -= self.__band[j + 1, i] * aSolution[j + i + 1]
aSolution[i] = rxw
return aSolution
## Invert band part.
#
# @return band part; matrix(float)
#
def invertBand():
nRow = self.__numBand + 1
nCol = self.__numCol
inverseBand = np.zeros((nRow, nCol))
for i in range(nCol - 1, -1, -1):
rxw = self.__band[0, i]
for j in range(i, max(0, i - nRow + 1) - 1, -1):
for k in range(j + 1, min(nCol, j + nRow)):
rxw -= inverseBand[abs(i - k), min(i, k)] * self.__band[k - j, j]
inverseBand[i - j, j] = rxw
rxw = 0.0
return inverseBand
## Calculate band part of A*V*A^T.
#
# @param anArray matrix A; matrix(float)
# @param aSymArray symmetric matrix V; matrix(float)
# @return band part; matrix(float)
#
def bandOfAVAT(anArray, aSymArray):
nCol = self.__numCol
nBand = self.__numBand
aBand = np.empty((nBand + 1, nCol))
for i in range(nCol):
for j in range(max(0, i - nBand), i + 1):
aBand[i - j, j] = np.dot(anArray[i], np.dot(aSymArray, anArray[j]))
return aBand
# setup
nBorder = self.__numBorder
# nRow = self.__numBand +1
nCol = self.__numCol
aSolution = np.empty(nBorder + nCol)
# decompose band
decomposeBand()
if ((self.__band[0] <= 0.0).any()):
raise ZeroDivisionError("Band matrix not positive definite")
if (nBorder > 0):
auxMat = np.empty((nBorder, nCol))
# solve for mixed part
for i in range(nBorder):
auxMat[i] = solveBand(self.__mixed[i])
auxMatT = auxMat.T
# solve for border
auxVec = aRightHandSide[:nBorder] - np.dot(auxMat, aRightHandSide[nBorder:])
invBorder = np.linalg.inv(self.__border - np.dot(self.__mixed, auxMatT))
aSolution[:nBorder] = np.dot(invBorder, auxVec)
# solve for band part
aSolution[nBorder:] = solveBand(aRightHandSide[nBorder:]) \
- np.dot(auxMatT, aSolution[:nBorder])
# parts of inverse
self.__border = invBorder
self.__mixed = np.dot(-invBorder, auxMat)
self.__band = invertBand() + bandOfAVAT(auxMatT, invBorder)
else:
# solve for band part (only)
aSolution = solveBand(aRightHandSide)
self.__band = invertBand()
return aSolution
'''
Input/output of MP-II binary records.
Created on Aug 1, 2011
@author: kleinwrt
'''
## \file
# Millepede-II (binary) record
#
# \author Claus Kleinwort, DESY, 2021 (Claus.Kleinwort@desy.de)
#
# \copyright
# Copyright (c) 2011 - 2018 Deutsches Elektronen-Synchroton,
# Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY \n\n
# This library is free software; you can redistribute it and/or modify
# it under the terms of the GNU Library General Public License as
# published by the Free Software Foundation; either version 2 of the
# License, or (at your option) any later version. \n\n
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Library General Public License for more details. \n\n
# You should have received a copy of the GNU Library General Public
# License along with this program (see the file COPYING.LIB for more
# details); if not, write to the Free Software Foundation, Inc.,
# 675 Mass Ave, Cambridge, MA 02139, USA.
import array, math
## Millepede-II (binary) record.
#
# Containing information for local (track) and global fit.
#
# The data blocks are collected in two arrays, a real array
# (containing float or double values) and integer array, of same length.
# A positive record length indicate _float_ and a negative one _double_ values.
# The content of the record is:
#
#\verbatim
# real array integer array
# 0 0.0 error count (this record)
# 1 RMEAS, measured value 0 __iMeas -+
# 2 local derivative index of local derivative |
# 3 local derivative index of local derivative |
# 4 ... | block
# SIGMA, error (>0) 0 __iErr |
# global derivative label of global derivative |
# global derivative label of global derivative |
# ... -+
# RMEAS, measured value 0 __position
# local derivative index of local derivative
# local derivative index of local derivative
# ...
# SIGMA, error 0
# global derivative label of global derivative
# global derivative label of global derivative
# ...
# global derivative label of global derivative __recLen
#\endverbatim
#
# Special data block (other/debug information).
# Contains no local derivatives and (error) SIGMA is negative (-Number of SPecial data words).
#
#\verbatim
# real array integer array
# 0.0 0 __iMeas -+
# -float(NSP) 0 __iErr |
# special data special data | special block (2+NSP words)
# special data special data |
# ... -+
#\endverbatim
#
class MilleRecord(object):
## Create MP-II binary record.
#
def __init__(self, doublePrec=False):
## flag for storage in as *double* values
self.__doublePrecision = doublePrec
## position in record, usually start of next data block; int
self.__position = 1
## number of data blocks in record; int
self.__numData = 0
## record length; int
self.__recLen = 0
## position of value in current data block; int
self.__iMeas = 0
## position of error in current data block; int
self.__iErr = 0
## array with markers (0) and labels; array(int32)
self.__inder = array.array('i')
## array with values, errors and derivatives; (float32 or float64)
self.__glder = array.array('d' if doublePrec else 'f')
## Add data block to (end of) record.
#
# @param dataList list with measurement, error, labels and derivatives; list
#
def addData(self, dataList):
if (self.__numData == 0): # first word is error counter
self.__inder.append(0)
self.__glder.append(0.)
self.__numData += 1
aMeas, aPrec, indLocal, derLocal, labGlobal, derGlobal = dataList
self.__inder.append(0)
self.__glder.append(aMeas)
self.__inder.fromlist(indLocal)
self.__glder.fromlist(derLocal)
self.__inder.append(0)
self.__glder.append(1.0 / math.sqrt(aPrec)) # convert to error
self.__inder.fromlist(labGlobal)
self.__glder.fromlist(derGlobal)
## Get data block from current position in record.
#
# @return list with measurement, error, labels and derivatives; list
#
def getData(self):
aMeas = self.__glder[self.__iMeas]
indLocal = []
derLocal = []
for i in range(self.__iMeas + 1, self.__iErr):
indLocal.append(self.__inder[i])
derLocal.append(self.__glder[i])
aPrec = 1.0 / self.__glder[self.__iErr] ** 2 # convert to precision
indGlobal = []
derGlobal = []
for i in range(self.__iErr + 1, self.__position):
indGlobal.append(self.__inder[i])
derGlobal.append(self.__glder[i])
return aMeas, aPrec, indLocal, derLocal, indGlobal, derGlobal
## Print record.
def printRecord(self):
print(" MilleRecord, len: ", len(self.__inder))
print(self.__inder)
print(self.__glder)
## Write record to file.
#
# @param aFile (binary) file
#
def writeRecord(self, aFile):
header = array.array('i') # header with number of words
header.append(-len(self.__inder) * 2 if self.__doublePrecision else len(self.__inder) * 2)
header.tofile(aFile)
self.__glder.tofile(aFile)
self.__inder.tofile(aFile)
## Read record from file.
#
# @param aFile (binary) file
#
def readRecord(self, aFile):
header = array.array('i') # header with number of words
header.fromfile(aFile, 1)
self.__recLen = abs(header[0] / 2)
if header[0] < 0:
self.__glder = array.array('d')
self.__glder.fromfile(aFile, self.__recLen)
self.__inder.fromfile(aFile, self.__recLen)
## Locate next data block.
#
# @return next block exists; bool
#
def moreData(self):
if (self.__position < self.__recLen):
while (self.__position < self.__recLen and self.__inder[self.__position] != 0):
self.__position += 1
self.__iMeas = self.__position
self.__position += 1
while (self.__position < self.__recLen and self.__inder[self.__position] != 0):
self.__position += 1
self.__iErr = self.__position
self.__position += 1
# special data?
if (self.__iMeas + 1 == self.__iErr and self.__glder[self.__iErr] < 0):
self.__position += int(-self.__glder[self.__iErr])
else:
while (self.__position < self.__recLen and self.__inder[self.__position] != 0):
self.__position += 1
self.__numData += 1
return True
else:
return False
## Get special data tag from block.
#
# @return tag or -1 for ordinary data block; int
#
def specialDataTag(self):
aTag = -1
if (self.__iMeas + 1 == self.__iErr and self.__glder[self.__iErr] < 0):
aTag = int(-self.__glder[self.__iErr] * 10. + 0.5) % 10
return aTag
This diff is collapsed.
'''
Simple Test Program for General Broken Lines.
Created on Jul 27, 2011
@author: kleinwrt
'''
## \file
# (technical) examples
#
# \author Claus Kleinwort, DESY, 2011 (Claus.Kleinwort@desy.de)
#
# \copyright
# Copyright (c) 2011 - 2021 Deutsches Elektronen-Synchroton,
# Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY \n\n
# This library is free software; you can redistribute it and/or modify
# it under the terms of the GNU Library General Public License as
# published by the Free Software Foundation; either version 2 of the
# License, or (at your option) any later version. \n\n
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Library General Public License for more details. \n\n
# You should have received a copy of the GNU Library General Public
# License along with this program (see the file COPYING.LIB for more
# details); if not, write to the Free Software Foundation, Inc.,
# 675 Mass Ave, Cambridge, MA 02139, USA.
import numpy as np
import math
import time
from gblpy3.gblfit import GblPoint, GblTrajectory
#
## Create points on initial trajectory, create trajectory from points,
# fit and write trajectory to MP-II binary file,
# get track parameter corrections and covariance matrix at points.
#
# Equidistant measurement layers and thin scatterers, propagation
# with simple jacobian (quadratic in arc length differences).
# Curvilinear system (T,U,V) or measurement system (I,J,K) as local
# coordinate system and (q/p, slopes, offsets) as local track parameters.
#
# This example simulates and refits tracks in a system of planar detectors
# with 2D measurements in a constant magnet field in Z direction using
# the curvilinear system as local system. The true track parameters are
# randomly smeared with respect to a (constant and straight) reference
# trajectory with direction (lambda, phi) and are used (only) for the
# on-the-fly simulation of the measurements and scatterers. The predictions
# from the reference trajectory are therefore always zero and the residuals
# needed (by addMeasurement) are equal to the measurements.
#
def example1():
#
np.random.seed(47117)
nTry = 1000 #: number of tries
nLayer = 5 #: number of detector layers
print(" Gbltst $Id$ ", nTry, nLayer)
start = time.process_time()
# track direction
sinLambda = 0.3
cosLambda = math.sqrt(1.0 - sinLambda ** 2)
sinPhi = 0.
cosPhi = math.sqrt(1.0 - sinPhi ** 2)
# Curvilinear system: track direction T, U = Z x T / |Z x T|, V = T x U
tuvDir = np.array([[cosLambda * cosPhi, cosLambda * sinPhi, sinLambda], \
[-sinPhi, cosPhi, 0.], \
[-sinLambda * cosPhi, -sinLambda * sinPhi, cosLambda]])
# measurement resolution
measErr = np.array([ 0.001, 0.001]) # 10 mu
measPrec = 1.0 / measErr ** 2
# scattering error