X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=STEER%2FAliCheb3D.cxx;h=6d93b1f305307c53c1066d35a5fec6061477daf2;hb=f7a1cc68313147ec921d4c82df1890abe00e4032;hp=1c24a5305021756f359eddbd647b902807bc5332;hpb=99adacae21d3c29dd0d759fbe4e816cfd7d557b2;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliCheb3D.cxx b/STEER/AliCheb3D.cxx index 1c24a530502..6d93b1f3053 100644 --- a/STEER/AliCheb3D.cxx +++ b/STEER/AliCheb3D.cxx @@ -13,186 +13,103 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* $Id$ */ - -// Author: ruben.shahoyan@cern.ch 09/09/2006 -// -//////////////////////////////////////////////////////////////////////////////// -// // -// AliCheb3D produces the interpolation of the user 3D->NDimOut arbitrary // -// function supplied in "void (*fcn)(float* inp,float* out)" format // -// either in a separate macro file or as a function pointer. // -// Only coefficients needed to guarantee the requested precision are kept. // -// // -// The user-callable methods are: // -// To create the interpolation use: // -// AliCheb3D(const char* funName, // name of the file with user function // -// or // -// AliCheb3D(void (*ptr)(float*,float*),// pointer on the user function // -// Int_t DimOut, // dimensionality of the function's output // -// Float_t *bmin, // lower 3D bounds of interpolation domain // -// Float_t *bmax, // upper 3D bounds of interpolation domain // -// Int_t *npoints, // number of points in each of 3 input // -// // dimension, defining the interpolation grid // -// Float_t prec=1E-6); // requested max.absolute difference between // -// // the interpolation and any point on grid // -// // -// To test obtained parameterization use the method // -// TH1* TestRMS(int idim,int npoints = 1000,TH1* histo=0); // -// it will compare the user output of the user function and interpolation // -// for idim-th output dimension and fill the difference in the supplied // -// histogram. If no histogram is supplied, it will be created. // -// // -// To save the interpolation data: // -// SaveData(const char* filename, Bool_t append ) // -// write text file with data. If append is kTRUE and the output file already // -// exists, data will be added in the end of the file. // -// Alternatively, SaveData(FILE* stream) will write the data to // -// already existing stream. // -// // -// To read back already stored interpolation use either the constructor // -// AliCheb3D(const char* inpFile); // -// or the default constructor AliCheb3D() followed by // -// AliCheb3D::LoadData(const char* inpFile); // -// // -// To compute the interpolation use Eval(float* par,float *res) method, with // -// par being 3D vector of arguments (inside the validity region) and res is // -// the array of DimOut elements for the output. // -// // -// If only one component (say, idim-th) of the output is needed, use faster // -// Float_t Eval(Float_t *par,int idim) method. // -// // -// void Print(option="") will print the name, the ranges of validity and // -// the absolute precision of the parameterization. Option "l" will also print // -// the information about the number of coefficients for each output // -// dimension. // -// // -// NOTE: during the evaluation no check is done for parameter vector being // -// outside the interpolation region. If there is such a risk, use // -// Bool_t IsInside(float *par) method. Chebyshev parameterization is not // -// good for extrapolation! // -// // -// For the properties of Chebyshev parameterization see: // -// H.Wind, CERN EP Internal Report, 81-12/Rev. // -// // -//////////////////////////////////////////////////////////////////////////////// - #include #include -#include #include +#include +#include +#include +#include +#include #include "AliCheb3D.h" -#include "AliLog.h" - - +#include "AliCheb3DCalc.h" ClassImp(AliCheb3D) -AliCheb3D::AliCheb3D(): - TNamed("", ""), - fDimOut(0), - fPrec(0.), - fChebCalc(), - fMaxCoefs(0), - fResTmp(0), - fGrid(0), - fUsrFunName(), - fUsrMacro(0) -{ - // Default constructor - Init0(); -} - -AliCheb3D::AliCheb3D(const char* inputFile): - TNamed("", ""), - fDimOut(0), - fPrec(0.), - fChebCalc(), - fMaxCoefs(0), - fResTmp(0), - fGrid(0), - fUsrFunName(), - fUsrMacro(0) +//__________________________________________________________________________________________ +AliCheb3D::AliCheb3D() : + fDimOut(0), + fPrec(0), + fChebCalc(1), + fMaxCoefs(0), + fResTmp(0), + fGrid(0), + fUsrFunName(""), + fUsrMacro(0) { - // Default constructor - Init0(); - LoadData(inputFile); + for (int i=3;i--;) fBMin[i] = fBMax[i] = fBScale[i] = fBOffset[i] = 0; } - - -AliCheb3D::AliCheb3D(FILE* stream): - TNamed("", ""), - fDimOut(0), - fPrec(0.), - fChebCalc(), - fMaxCoefs(0), - fResTmp(0), - fGrid(0), - fUsrFunName(), - fUsrMacro(0) +//__________________________________________________________________________________________ +AliCheb3D::AliCheb3D(const AliCheb3D& src) : + TNamed(src), + fDimOut(src.fDimOut), + fPrec(src.fPrec), + fChebCalc(1), + fMaxCoefs(src.fMaxCoefs), + fResTmp(0), + fGrid(0), + fUsrFunName(src.fUsrFunName), + fUsrMacro(0) { - // Default constructor - Init0(); - LoadData(stream); + // read coefs from text file + for (int i=3;i--;) { + fBMin[i] = src.fBMin[i]; + fBMax[i] = src.fBMax[i]; + fBScale[i] = src.fBScale[i]; + fBOffset[i] = src.fBOffset[i]; + fNPoints[i] = src.fNPoints[i]; + } + for (int i=0;i%dD function. Precision: %e\n",GetName(),fDimOut,fPrec); printf("Region of validity: [%+.5e:%+.5e] [%+.5e:%+.5e] [%+.5e:%+.5e]\n",fBMin[0],fBMax[0],fBMin[1],fBMax[1],fBMin[2],fBMax[2]); TString opts = opt; opts.ToLower(); @@ -262,22 +287,7 @@ void AliCheb3D::Print(Option_t* opt) const } //__________________________________________________________________________________________ -void AliCheb3D::Init0() -{ - // Initialisation - for (int i=3;i--;) fBMin[i] = fBMax[i] = fBScale[i] = fBOffset[i] = 0; - fMaxCoefs = 0; - fGrid = 0; - fResTmp = 0; - fUsrFunName = ""; - fUsrMacro = 0; -#ifdef _INC_CREATION_ALICHEB3D_ - gUsrFunAliCheb3D = 0; -#endif -} - -//__________________________________________________________________________________________ -void AliCheb3D::PrepareBoundaries(Float_t *bmin,Float_t *bmax) +void AliCheb3D::PrepareBoundaries(const Float_t *bmin, const Float_t *bmax) { // Set and check boundaries defined by user, prepare coefficients for their conversion to [-1:1] interval // @@ -295,8 +305,20 @@ void AliCheb3D::PrepareBoundaries(Float_t *bmin,Float_t *bmax) // } + //__________________________________________________________________________________________ #ifdef _INC_CREATION_ALICHEB3D_ + +// Pointer on user function (faster altrnative to TMethodCall) +void (*gUsrFunAliCheb3D) (float* ,float* ); + +void AliCheb3D::EvalUsrFunction() +{ + // call user supplied function + if (gUsrFunAliCheb3D) gUsrFunAliCheb3D(fArgsTmp,fResTmp); + else fUsrMacro->Execute(); +} + void AliCheb3D::SetUsrFunction(const char* name) { // load user macro with function definition and compile it @@ -324,6 +346,8 @@ void AliCheb3D::SetUsrFunction(const char* name) #ifdef _INC_CREATION_ALICHEB3D_ void AliCheb3D::SetUsrFunction(void (*ptr)(float*,float*)) { + // assign user training function + // if (fUsrMacro) delete fUsrMacro; fUsrMacro = 0; fUsrFunName = ""; @@ -333,7 +357,10 @@ void AliCheb3D::SetUsrFunction(void (*ptr)(float*,float*)) //__________________________________________________________________________________________ #ifdef _INC_CREATION_ALICHEB3D_ -void AliCheb3D::EvalUsrFunction(Float_t *x, Float_t *res) { +void AliCheb3D::EvalUsrFunction(const Float_t *x, Float_t *res) +{ + // evaluate user function value + // for (int i=3;i--;) fArgsTmp[i] = x[i]; if (gUsrFunAliCheb3D) gUsrFunAliCheb3D(fArgsTmp,fResTmp); else fUsrMacro->Execute(); @@ -343,7 +370,7 @@ void AliCheb3D::EvalUsrFunction(Float_t *x, Float_t *res) { //__________________________________________________________________________________________ #ifdef _INC_CREATION_ALICHEB3D_ -Int_t AliCheb3D::CalcChebCoefs(Float_t *funval,int np, Float_t *outCoefs, Float_t prec) +Int_t AliCheb3D::CalcChebCoefs(const Float_t *funval,int np, Float_t *outCoefs, Float_t prec) { // Calculate Chebyshev coeffs using precomputed function values at np roots. // If prec>0, estimate the highest coeff number providing the needed precision @@ -383,12 +410,14 @@ void AliCheb3D::DefineGrid(Int_t* npoints) for (int id=3;id--;) { fNPoints[id] = npoints[id]; if (fNPoints[id]=fracStep) { + frac = fr; + printf("\b\b\b\b\b\b\b\b\b\b\b"); + printf("%05.2f%% Done",fr*100); + fflush(stdout); + } + // } int nc = CalcChebCoefs(fvals,fNPoints[0], tmpCoef1D, fPrec); for (int id0=fNPoints[0];id0--;) tmpCoef2D[id1 + id0*fNPoints[1]] = tmpCoef1D[id0]; @@ -474,8 +517,7 @@ Int_t AliCheb3D::ChebFit(int dmOut) for (int id2=fNPoints[2];id2--;) { int id = id2 + fNPoints[2]*(id1+id0*fNPoints[1]); Float_t cfa = TMath::Abs(tmpCoef3D[id]); - if (cfa < RTiny) {tmpCoef3D[id] = 0; continue;} // neglect coeefs below the threshold - + if (cfa < rTiny) {tmpCoef3D[id] = 0; continue;} // neglect coefs below the threshold resid += cfa; if (resid0 && tmpCoefSurf[(id1-1)+id0*fNPoints[1]]==0) id1--; tmpCols[id0] = id1; } // find max significant row - for (int id0=NRows;id0--;) {if (tmpCols[id0]>0) break; NRows--;} + for (int id0=nRows;id0--;) {if (tmpCols[id0]>0) break; nRows--;} // find max significant column and fill the permanent storage for the max sigificant column of each row - cheb->InitRows(NRows); // create needed arrays; - int *NColsAtRow = cheb->GetNColsAtRow(); - int *ColAtRowBg = cheb->GetColAtRowBg(); - int NCols = 0; + cheb->InitRows(nRows); // create needed arrays; + int *nColsAtRow = cheb->GetNColsAtRow(); + int *colAtRowBg = cheb->GetColAtRowBg(); + int nCols = 0; int NElemBound2D = 0; - for (int id0=0;id0InitCols(NCols); + cheb->InitCols(nCols); delete[] tmpCols; // // create the 2D matrix defining the boundary of significance for 3D coeffs.matrix // and count the number of siginifacnt coefficients // cheb->InitElemBound2D(NElemBound2D); - int *CoefBound2D0 = cheb->GetCoefBound2D0(); - int *CoefBound2D1 = cheb->GetCoefBound2D1(); + int *coefBound2D0 = cheb->GetCoefBound2D0(); + int *coefBound2D1 = cheb->GetCoefBound2D1(); fMaxCoefs = 0; // redefine number of coeffs - for (int id0=0;id0InitCoefs(fMaxCoefs); - Float_t *Coefs = cheb->GetCoefs(); + Float_t *coefs = cheb->GetCoefs(); int count = 0; - for (int id0=0;id0SaveData(stream); fprintf(stream,"#\nEND %s\n#\n",GetName()); @@ -614,7 +658,8 @@ void AliCheb3D::SaveData(FILE* stream) const //_______________________________________________ void AliCheb3D::LoadData(const char* inpFile) { - // Load data from input file + // load coefficients data from txt file + // TString strf = inpFile; gSystem->ExpandPathName(strf); FILE* stream = fopen(strf.Data(),"r"); @@ -626,7 +671,8 @@ void AliCheb3D::LoadData(const char* inpFile) //_______________________________________________ void AliCheb3D::LoadData(FILE* stream) { - // Load data from input stream stream + // load coefficients data from stream + // if (!stream) {Error("LoadData","No stream provided.\nStop"); exit(1);} TString buffs; Clear(); @@ -667,12 +713,12 @@ void AliCheb3D::LoadData(FILE* stream) } //_______________________________________________ -void AliCheb3D::SetDimOut(int d) +void AliCheb3D::SetDimOut(const int d) { - // Set the dimension of the output array + // init output dimensions fDimOut = d; if (fResTmp) delete fResTmp; - fResTmp = new Float_t[fDimOut]; // RRR + fResTmp = new Float_t[fDimOut]; fChebCalc.Delete(); for (int i=0;i2) {printf("Maximum 3 dimensions are supported\n"); return;} fBMin[id] += dif; fBMax[id] += dif; @@ -712,3 +759,82 @@ TH1* AliCheb3D::TestRMS(int idim,int npoints,TH1* histo) // } #endif + +//_______________________________________________ +#ifdef _INC_CREATION_ALICHEB3D_ +void AliCheb3D::EstimateNPoints(float Prec, int gridBC[3][3]) +{ + // Estimate number of points to generate a training data + // + const float sclA[9] = {0.1, 0.5, 0.9, 0.1, 0.5, 0.9, 0.1, 0.5, 0.9} ; + const float sclB[9] = {0.1, 0.1, 0.1, 0.5, 0.5, 0.5, 0.9, 0.9, 0.9} ; + const float sclDim[2] = {0.01,0.99}; + const int compDim[3][2] = { {1,2}, {2,0}, {0,1} }; + static float xyz[3]; + // + for (int i=3;i--;)for (int j=3;j--;) gridBC[i][j] = -1; + // + for (int idim=0;idim<3;idim++) { + float dimMN = fBMin[idim] + sclDim[0]*(fBMax[idim]-fBMin[idim]); + float dimMX = fBMin[idim] + sclDim[1]*(fBMax[idim]-fBMin[idim]); + // + for (int it=0;it<9;it++) { // test in 9 points + int id1 = compDim[idim][0]; // 1st fixed dim + int id2 = compDim[idim][1]; // 2nd fixed dim + xyz[ id1 ] = fBMin[id1] + sclA[it]*( fBMax[id1]-fBMin[id1] ); + xyz[ id2 ] = fBMin[id2] + sclB[it]*( fBMax[id2]-fBMin[id2] ); + // + int* npt = GetNCNeeded(xyz,idim, dimMN,dimMX, Prec); // npoints for Bx,By,Bz + for (int ib=0;ib<3;ib++) if (npt[ib]>gridBC[ib][idim]) gridBC[ib][idim] = npt[ib]+2; + // + } + } +} + + +int* AliCheb3D::GetNCNeeded(float xyz[3],int DimVar, float mn,float mx, float prec) +{ + // estimate needed number of chebyshev coefs for given function desctiption in DimVar dimension + // The values for two other dimensions must be set beforehand + // + static int curNC[3]; + static int retNC[3]; + const int kMaxPoint = 400; + float* gridVal = new float[3*kMaxPoint]; + float* coefs = new float[3*kMaxPoint]; + // + float scale = mx-mn; + float offs = mn + scale/2.0; + scale = 2./scale; + // + int curNP; + int maxNC=-1; + int maxNCPrev=-1; + for (int i=0;i<3;i++) retNC[i] = -1; + for (int i=0;i<3;i++) fArgsTmp[i] = xyz[i]; + // + for (curNP=5; curNP3 && (maxNC-maxNCPrev)<1 ) break; + maxNCPrev = maxNC; + // + } + delete[] gridVal; + delete[] coefs; + return retNC; + // +} +#endif