diff --git a/.gitignore b/.gitignore index a89409e03670800ce83e4d919efdb1f9cfb85955..62910df9a4ef3a51a5a3f462a08c0dade1a7c94f 100644 --- a/.gitignore +++ b/.gitignore @@ -18,3 +18,5 @@ *.swp .root_hist /example/temp.root +*.d +*.so diff --git a/example/TNArray/test_TNArray.C b/example/TNArray/test_TNArray.C new file mode 100644 index 0000000000000000000000000000000000000000..c68a0b7c0a570620884bb930395dc85f128a0bce --- /dev/null +++ b/example/TNArray/test_TNArray.C @@ -0,0 +1,40 @@ +#include "TNArray.h" +#include "TF2.h" +#include "TMath.h" +#include "TRandom.h" +#include "TNtupleD.h" +void test_TNArray() { + art::TNArray *array = new art::TNArray("table","gaus2d","a:r:value"); + array->Add("angle",0.,90.,181); + array->Add("radius",0.,10,501); + double tstep = array->GetVar(0).GetStep(); + double rstep = array->GetVar(1).GetStep(); + + for (int it = 0, nt = array->GetVar(0).GetNumVals(); it < nt; ++it) { + for (int ir = 0, nr = array->GetVar(1).GetNumVals(); ir < nr; ++ir) { + double r = ir * rstep; + double adeg = it * tstep; + double a = adeg * TMath::DegToRad(); + double x = r * TMath::Cos(a); + double y = r * TMath::Sin(a); + array->Fill(adeg,r, TMath::Gaus(x,0.,1.,1)*TMath::Gaus(y,0.,1.,1) ); + } + } + array->Load(); + + TNtupleD *tuple = new TNtupleD("tuple","tuple","x:y:r:a:tval:tval2:ival:diff:diff2"); + + for (int ix = 0; ix < 100; ++ix) { + for (int iy = 0; iy < 100; ++iy) { + double x = gRandom->Uniform(0.,5.); + double y = gRandom->Uniform(0.,5.); + double r = TMath::Sqrt(x*x + y*y); + double a = TMath::ATan(y/x) * TMath::RadToDeg(); + double xval [] = {a,r}; + double tval = array->Eval(xval); + double tval2 = array->Eval2(xval); + double ival = TMath::Gaus(x,0.,1.,1)*TMath::Gaus(y,0.,1.,1); + tuple->Fill(x,y,r,a,tval,tval2,ival,tval-ival,tval2-ival); + } + } +} diff --git a/sources/cont/TNArray.cc b/sources/cont/TNArray.cc index 00fdd5a91b589c9dd1b0811b674d57d393074334..1cfad72af24a922d29dae540dc1b124e7de4bcd9 100644 --- a/sources/cont/TNArray.cc +++ b/sources/cont/TNArray.cc @@ -3,13 +3,16 @@ * @brief * * @date Created : 2016-01-29 14:16:43 JST - * Last Modified : 2018-05-07 17:41:40 JST (ota) + * Last Modified : 2020-05-07 18:43:03 JST (ota) * @author Shinsuke OTA * * (C) 2016 Shinsuke OTA */ #include "TNArray.h" +#include "TMatrixD.h" +#include "TVectorD.h" +#include "TDecompChol.h" using art::TNArray; @@ -106,6 +109,170 @@ Double_t TNArray::Eval(Double_t *x) return retval; } + +void TNArray::CalcBilinearWeights(double dx, std::vector& weight) { + weight.resize(2); + weight[0] = 1 - dx; + weight[1] = dx; +} + +void TNArray::CalcBicubicWeights(double dx, std::vector& weight, double param) { + weight.resize(4); + double d = 1 + dx; + weight[0] = param * (d * d * d - 5 * d * d + 8 * d - 4); + d = dx; + weight[1] = (param + 2) * d * d * d - (param + 3) * d * d + 1; + d = 1 - dx; + weight[2] = (param + 2) * d * d * d - (param + 3) * d * d + 1; + d= 2 - dx; + weight[3] = param * (d * d * d - 5 * d * d + 8 * d - 4); +} + +double TNArray::WeightedSum(int ipar, double w, int idxbase, + const std::vector< std::vector > &weight, + const std::vector< std::vector > &indexes) +{ + double sum = 0.; + for (int i = 0, n = weight[ipar].size(); i < n; ++i) { + double this_w = w * weight[ipar][i]; + int idx = idxbase + indexes[ipar][i] * fTotalNumVals[ipar]; + if (ipar != 0) { + sum += WeightedSum(ipar - 1,this_w,idx, weight, indexes); + } else { +// printf("ipar = %d, idx = %10d, w = %10.5f, wprod = %10.5f\n",ipar,idx,weight[ipar][i],w); + + sum += this_w * fValues[idx]; +// printf("ipar = %d, idx = %10d, w = %10.5f, wprod = %10.5f, value = %f\n",ipar,idx,weight[ipar][i],w,fValues[idx]); + + } +// printf("ipar = %2d, ip = %2d, w = %10.5f, sum = %10.5f\n",ipar,i, weight[ipar][i], sum); + } + return sum; +} + + +Double_t TNArray::Eval2(Double_t *x) +{ +// Int_t baseIndex = 0; +// Double_t retval = 0.; + std::vector idx(fNumVars); + std::vector deriv(fNumVars); + std::vector< std::vector< double > > weight(fNumVars); + std::vector< std::vector< int > > indexes(fNumVars); + + for (Int_t i = 0; i!=fNumVars; ++i) { + if (x[i] < fVars[i].GetMin()) { + return TMath::QuietNaN(); + } else if (fVars[i].GetMax() < x[i]) { + return TMath::QuietNaN(); + } + idx[i] = fVars[i].IndexI(x[i]); + deriv[i] = fVars[i].Derivative(x[i]); +// printf("ipar = %d, idx = %d, deriv = %f\n",i,idx[i],deriv[i]); + + + if (0 < idx[i] && idx[i] + 2 < fVars[i].GetNumVals()) { + CalcBicubicWeights(deriv[i],weight[i]); + indexes[i].push_back(idx[i] - 1); + indexes[i].push_back(idx[i]); + indexes[i].push_back(idx[i] + 1); + indexes[i].push_back(idx[i] + 2); + } else if (0 == idx[i] || idx[i] + 1 < fVars[i].GetNumVals()) { + CalcBilinearWeights(deriv[i],weight[i]); + indexes[i].push_back(idx[i]); + indexes[i].push_back(idx[i] + 1); + } else if (idx[i] + 1 == fVars[i].GetNumVals()) { + weight[i].push_back(1.); + indexes[i].push_back(idx[i]); + } + } + + return WeightedSum(fNumVars - 1, 1., 0, weight, indexes); + +#if 0 + printf("baseIndex = %d\n",baseIndex); + Double_t reference = retval = fValues[baseIndex]; + + for (Int_t i = 0; i!=fNumVars; ++i) { + Double_t diff = 0; + Int_t nval = 1; + int idxdiff = 1; + const int npt = 8; + double node[npt]; + double node2[npt]; + double node3[npt]; + double node4[npt]; + double val[npt]; + + node[0] = 0; + node2[0] = 0; + val[0] = 0.; + + + + while (nval < npt) { + if (idx[i] + idxdiff < fVars[i].GetNumVals()) { + node[nval] = idxdiff; + node2[nval] = idxdiff * idxdiff; + node3[nval] = idxdiff * idxdiff * idxdiff; + node4[nval] = idxdiff * idxdiff * idxdiff * idxdiff; + val[nval] = fValues[baseIndex + idxdiff * fTotalNumVals[i]] - reference; + nval++; + if (nval >= npt) break; + } + if (idx[i] - idxdiff > 0) { + node[nval] = - idxdiff; + node2[nval] = idxdiff * idxdiff; + node3[nval] = idxdiff * idxdiff * idxdiff; + node4[nval] = idxdiff * idxdiff * idxdiff * idxdiff; + val[nval] = fValues[baseIndex - idxdiff * fTotalNumVals[i]] - reference; + nval++; + } + idxdiff++; + } + TVectorD x; x.Use(npt,node); + TVectorD x2; x2.Use(npt,node2); + TVectorD x3; x3.Use(npt,node3); + TVectorD x4; x4.Use(npt,node4); + TVectorD y; y.Use(npt,val); + TMatrixD A(npt,5); // 4th order polynomial; + TMatrixDColumn(A,0) = 1.0; + TMatrixDColumn(A,1) = x; + TMatrixDColumn(A,2) = x2; + TMatrixDColumn(A,3) = x3; + TMatrixDColumn(A,4) = x4; + const TVectorD c_norm = NormalEqn(A,y); + const double deriv = fVars[i].Derivative(x[i]); + + retval += c_norm[0] + deriv * (c_norm[1] + deriv * (c_norm[2] + deriv * (c_norm[3] + deriv * (c_norm[4])))); + +#if 0 + for (int i = 0; i < npt; ++i) { + printf("[%d] node = %5.1f, node2 = %5.1f, val = %10.5f\n",i,node[i],node2[i],val[i]); + } + A.Print(); + y.Print(); + c_norm.Print(); +#endif + +#if 0 + printf("diff[%d] = %f %f %f @ idx=%d\n",i,diff,fVars[i].Derivative(x[i]),fValues[baseIndex+fTotalNumVals[i]],baseIndex+fTotalNumVals[i]); + printf(" min = %g\n",fVars[i].GetMin()); + printf(" x-min = %g step = %g\n",x[i]-fVars[i].GetMin(),fVars[i].GetStep()); + printf(" x-min = %g\n",((x[i]-fVars[i].GetMin())*(1./fVars[i].GetStep()))); + printf(" x-min = %d\n",(Int_t)((x[i]-fVars[i].GetMin())*(1./fVars[i].GetStep()))); + printf(" IndexI = %d\n",fVars[i].IndexI(x[i])); + +#endif + } + return retval; +#endif +} + + + + + void TNArray::Load() { Int_t nEntries = GetEntries(); @@ -144,8 +311,10 @@ void TNArray::Load() printf("index = %d used for entry %d\n",index,iEntry); for (Int_t iVar = 0; iVar != nVars; ++iVar) { printf("iVar = %d vars[%d] = %.10g IndexI = %d\n",iVar,iVar,vars[iVar],fVars[iVar].IndexI(vars[iVar])); - printf(" %.20e / %.20f = %.20e => %.20e\n",vars[iVar],fVars[iVar].GetStep(),vars[iVar]/fVars[iVar].GetStep(), - TMath::Floor(((Float_t)vars[iVar])/((Float_t)fVars[iVar].GetStep()))); + printf(" %.20e / %.20f = %.20e => %.20e\n",(vars[iVar]-fVars[iVar].GetMin()),fVars[iVar].GetStep(),(vars[iVar]-fVars[iVar].GetMin())/fVars[iVar].GetStep(), + TMath::Floor((vars[iVar]-fVars[iVar].GetMin())/(fVars[iVar].GetStep()))); + printf("index = %d\n",fVars[iVar].IndexI(vars[iVar])); + } return; } @@ -169,6 +338,8 @@ void TNArray::Load() } + + TNArray::Variable::Variable(const char* name, Double_t step, Double_t min, Double_t max, Int_t num) : TNamed(name,name), fStep(step), fMin(min), fMax(max), fNumVals(num) { diff --git a/sources/cont/TNArray.h b/sources/cont/TNArray.h index 2e7d697d13047abe58dca0af63faf3988ce0bf2f..d6e88d9299b9654a09e1f91a1d1f7c2f43539c20 100644 --- a/sources/cont/TNArray.h +++ b/sources/cont/TNArray.h @@ -3,7 +3,7 @@ * @brief n-dimension array * * @date Created : 2016-01-29 11:34:04 JST - * Last Modified : 2018-04-23 23:00:03 JST (ota) + * Last Modified : 2020-05-07 18:15:50 JST (ota) * @author Shinsuke OTA * * (C) 2016 Shinsuke OTA @@ -35,11 +35,18 @@ public: virtual void Add(const char* name, Double_t min, Double_t max, Int_t num); virtual Double_t Eval(Double_t *x); + virtual Double_t Eval2(Double_t *x); virtual const Variable& GetVar(Int_t i) { return fVars[i]; } virtual Double_t Value(Int_t i) const { return fValues[i]; } + virtual double WeightedSum(int ipar, double w, int idx, const std::vector< std::vector< double > >& weights, + const std::vector< std::vector< int > > & indexes); + virtual void CalcBicubicWeights(double dx, std::vector& weight, double param = -0.5); + virtual void CalcBilinearWeights(double dx, std::vector& weight); + + protected: private: @@ -59,9 +66,12 @@ public: // Int_t IndexI(Double_t x) { return TMath::Nint((x-fMin)*(1./fStep)); } // Int_t IndexI(Double_t x) { return (Int_t)TMath::Floor(((Float_t)(x-fMin))*((Float_t)(1./fStep))); } - Int_t IndexI(Double_t x) { return TMath::FloorNint((x-fMin)/fStep + 0.5); } + Int_t IndexI(Double_t x) { return TMath::FloorNint(TMath::Floor((x-fMin)/fStep*10 + 0.5)/10); } Bool_t CheckBounce(Double_t x) { return (fMin <= x ) && (x <= fMax); } - Double_t Derivative(Double_t x) { return (x-fMin)*(1/fStep) - TMath::Floor((x-fMin)*(1/fStep) + 0.5); } + Double_t Derivative(Double_t x) { + double ret = (x-fMin)*(1/fStep) - TMath::Floor((TMath::Floor((x-fMin)/fStep * 10 + 0.5))/10.); + return ret > 0. ? ret : 0.; + } virtual Int_t GetNumVals() const { return fNumVals; } virtual Double_t GetMin() const{ return fMin; }