//+------------------------------------------------------------------+ //| dataanalysis.mqh | //| Copyright 2003-2022 Sergey Bochkanov (ALGLIB project) | //| Copyright 2012-2026, MetaQuotes Ltd. | //| www.mql5.com | //+------------------------------------------------------------------+ //| Implementation of ALGLIB library in MetaQuotes Language 5 | //| | //| The features of the library include: | //| - Linear algebra (direct algorithms, EVD, SVD) | //| - Solving systems of linear and non-linear equations | //| - Interpolation | //| - Optimization | //| - FFT (Fast Fourier Transform) | //| - Numerical integration | //| - Linear and nonlinear least-squares fitting | //| - Ordinary differential equations | //| - Computation of special functions | //| - Descriptive statistics and hypothesis testing | //| - Data analysis - classification, regression | //| - Implementing linear algebra algorithms, interpolation, etc. | //| in high-precision arithmetic (using MPFR) | //| | //| This file is free software; you can redistribute it and/or | //| modify it under the terms of the GNU General Public License as | //| published by the Free Software Foundation (www.fsf.org); either | //| version 2 of the License, or (at your option) any later version. | //| | //| This program 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 General Public License for more details. | //+------------------------------------------------------------------+ #include "ap.mqh" #include "optimization.mqh" #include "statistics.mqh" #include "solvers.mqh" //+------------------------------------------------------------------+ //| Auxiliary class for CBdSS | //+------------------------------------------------------------------+ struct CCVReport { double m_RelCLSError; double m_AvgCE; double m_RMSError; double m_AvgError; double m_AvgRelError; //--- CCVReport(void) { ZeroMemory(this); } ~CCVReport(void) {} //--- void Copy(const CCVReport &obj); //--- overloading void operator=(const CCVReport &obj) { Copy(obj); } }; //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ void CCVReport::Copy(const CCVReport &obj) { m_RelCLSError=obj.m_RelCLSError; m_AvgCE=obj.m_AvgCE; m_RMSError=obj.m_RMSError; m_AvgError=obj.m_AvgError; m_AvgRelError=obj.m_AvgRelError; } //+------------------------------------------------------------------+ //| Data analysis | //+------------------------------------------------------------------+ class CBdSS { public: static void DSErrAllocate(const int nclasses,double &buf[]); static void DSErrAllocate(const int nclasses,CRowDouble &buf); static void DSErrAccumulate(double &buf[],double &y[],double &desiredy[]); static void DSErrAccumulate(CRowDouble &buf,CRowDouble &y,CRowDouble &desiredy); static void DSErrFinish(double &buf[]); static void DSErrFinish(CRowDouble &buf); static void DSNormalize(CMatrixDouble &xy,const int npoints,const int nvars,int &info,double &means[],double &sigmas[]); static void DSNormalize(CMatrixDouble &xy,const int npoints,const int nvars,int &info,CRowDouble &means,CRowDouble &sigmas); static void DSNormalizeC(CMatrixDouble &xy,const int npoints,const int nvars,int &info,double &means[],double &sigmas[]); static void DSNormalizeC(CMatrixDouble &xy,const int npoints,const int nvars,int &info,CRowDouble &means,CRowDouble &sigmas); static double DSGetMeanMindIstance(CMatrixDouble &xy,const int npoints,const int nvars); static void DSTie(double &a[],const int n,int &ties[],int &tiecount,int &p1[],int &p2[]); static void DSTie(CRowDouble &a,const int n,CRowInt &ties,int &tiecount,CRowInt &p1,CRowInt &p2); static void DSTieFastI(double &a[],int &b[],const int n,int &ties[],int &tiecount,double &bufr[],int &bufi[]); static void DSTieFastI(CRowDouble &a,CRowInt &b,const int n,CRowInt &ties,int &tiecount,CRowDouble &bufr,CRowInt &bufi); static void DSOptimalSplit2(double &ca[],int &cc[],const int n,int &info,double &threshold,double &pal,double &pbl,double &par,double &pbr,double &cve); static void DSOptimalSplit2(CRowDouble &ca,CRowInt &cc,const int n,int &info,double &threshold,double &pal,double &pbl,double &par,double &pbr,double &cve); static void DSOptimalSplit2Fast(double &a[],int &c[],int &tiesbuf[],int &cntbuf[],double &bufr[],int &bufi[],const int n,const int nc,double alpha,int &info,double &threshold,double &rms,double &cvrms); static void DSOptimalSplit2Fast(CRowDouble &a,CRowInt &c,CRowInt &tiesbuf,CRowInt &cntbuf,CRowDouble &bufr,CRowInt &bufi,const int n,const int nc,double alpha,int &info,double &threshold,double &rms,double &cvrms); static void DSSplitK(double &ca[],int &cc[],const int n,const int nc,int kmax,int &info,double &thresholds[],int &ni,double &cve); static void DSSplitK(CRowDouble &ca,CRowInt &cc,const int n,const int nc,int kmax,int &info,CRowDouble &thresholds,int &ni,double &cve); static void DSOptimalSplitK(double &ca[],int &cc[],const int n,const int nc,int kmax,int &info,double &thresholds[],int &ni,double &cve); static void DSOptimalSplitK(CRowDouble &ca,CRowInt &cc,const int n,const int nc,int kmax,int &info,CRowDouble &thresholds,int &ni,double &cve); private: static double XLnY(const double x,const double y); static double GetCV(CRowInt &cnt,const int nc); static void TieAddC(CRowInt &c,CRowInt &ties,const int ntie,const int nc,CRowInt &cnt); static void TieSubC(CRowInt &c,CRowInt &ties,const int ntie,const int nc,CRowInt &cnt); }; //+------------------------------------------------------------------+ //| This set of routines (DSErrAllocate, DSErrAccumulate, | //| DSErrFinish) calculates different error functions (classification| //| error, cross-entropy, rms, avg, avg.rel errors). | //| 1. DSErrAllocate prepares buffer. | //| 2. DSErrAccumulate accumulates individual errors: | //| * Y contains predicted output (posterior probabilities for | //| classification) | //| * DesiredY contains desired output (class number for | //| classification) | //| 3. DSErrFinish outputs results: | //| * Buf[0] contains relative classification error (zero for | //| regression tasks) | //| * Buf[1] contains avg. cross-entropy (zero for regression | //| tasks) | //| * Buf[2] contains rms error (regression, classification) | //| * Buf[3] contains average error (regression, classification) | //| * Buf[4] contains average relative error (regression, | //| classification) | //| NOTES(1): | //| "NClasses>0" means that we have classification task. | //| "NClasses<0" means regression task with -NClasses real | //| outputs. | //| NOTES(2): | //| rms. avg, avg.rel errors for classification tasks are | //| interpreted as errors in posterior probabilities with | //| respect to probabilities given by training/test set. | //+------------------------------------------------------------------+ void CBdSS::DSErrAllocate(const int nclasses,double &buf[]) { CRowDouble Buf; DSErrAllocate(nclasses,Buf); Buf.ToArray(buf); } //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ void CBdSS::DSErrAllocate(const int nclasses,CRowDouble &buf) { //--- allocation buf=vector::Zeros(8); //--- initialization buf.Set(5,(double)nclasses); } //+------------------------------------------------------------------+ //| See DSErrAllocate for comments on this routine. | //+------------------------------------------------------------------+ void CBdSS::DSErrAccumulate(double &buf[],double &y[],double &desiredy[]) { CRowDouble Buf=buf; CRowDouble Y=y; CRowDouble DesiredY=desiredy; DSErrAccumulate(Buf,Y,DesiredY); Buf.ToArray(buf); } //+------------------------------------------------------------------+ //| See DSErrAllocate for comments on this routine. | //+------------------------------------------------------------------+ void CBdSS::DSErrAccumulate(CRowDouble &buf,CRowDouble &y,CRowDouble &desiredy) { //--- create variables int offs=5; int nclasses=(int)MathRound(buf[offs]); int nout=0; int mmax=0; int rmax=0; int j=0; double v=0; //--- initialization vector ev=vector::Zeros(MathAbs(nclasses)); //--- check if(nclasses>0) { //--- Classification rmax=(int)MathRound(desiredy[0]); //--- initialization mmax=(int)y.ArgMax(); //--- check if(mmax!=rmax) buf.Add(0,1); //--- check if(y[rmax]>0.0) buf.Add(1,-MathLog(y[rmax])); else buf.Add(1,MathLog(CMath::m_maxrealnumber)); //--- calculation if(rmax>=0 && rmax err=y.ToVector()-ev; buf.Add(2,MathPow(err,2.0).Sum()); buf.Add(3,MathAbs(err).Sum()); //--- change value buf.Add(offs+1,1); } else { //--- Regression nout=-nclasses; //--- initialization rmax=(int)desiredy.ArgMax(); //--- initialization mmax=(int)y.ArgMax(); //--- check if(mmax!=rmax) buf.Add(0,1); //--- calculation ev=desiredy.ToVector(); vector err=y.ToVector()-ev; buf.Add(2,MathPow(err,2.0).Sum()); buf.Add(3,MathAbs(err).Sum()); for(j=0; j tmp=xy.ToMatrix(); tmp.Resize(npoints,nvars); //--- change value info=1; //--- Standartization means=tmp.Mean(0); sigmas=tmp.Std(0); //--- calculation for(int j=0; j tmp; vector tmp2; //--- Test parameters if(npoints<=0 || nvars<1) return(0); //--- Process tmp=vector::Full(npoints,CMath::m_maxrealnumber); //--- allocation for(int i=0; i1) | //| * -1, incorrect pararemets were passed (N<=0). | //| * 1, OK | //| Threshold- partiton boundary. Left part contains values | //| which are strictly less than Threshold. Right | //| part contains values which are greater than or | //| equal to Threshold. | //| PAL, PBL- probabilities P(0|v=Threshold) and | //| P(1|v>=Threshold) | //| CVE - cross-validation estimate of cross-entropy | //+------------------------------------------------------------------+ void CBdSS::DSOptimalSplit2(double &ca[],int &cc[],const int n, int &info,double &threshold,double &pal, double &pbl,double &par,double &pbr, double &cve) { CRowDouble A=ca; CRowInt C=cc; DSOptimalSplit2(A,C,n,info,threshold,pal,pbl,par,pbr,cve); } //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ void CBdSS::DSOptimalSplit2(CRowDouble &ca,CRowInt &cc,const int n, int &info,double &threshold,double &pal, double &pbl,double &par,double &pbr, double &cve) { //--- create variables int i=0; int t=0; double s=0; int tiecount=0; int k=0; int koptimal=0; double pak=0; double pbk=0; double cvoptimal=0; double cv=0; //--- creating arrays CRowInt ties; CRowInt p1; CRowInt p2; //--- copy CRowDouble a=ca; CRowInt c=cc; //--- initialization info=0; threshold=0; pal=0; pbl=0; par=0; pbr=0; cve=0; //--- Test for errors in inputs if(n<=0) { info=-1; return; } for(i=0; i 1 //--- NOTE: we assume that P.Get(i,j) equals to 0 or 1, //--- intermediate values are not allowed. pal=0; pbl=0; par=0; pbr=0; for(i=0; i0"=OK, "<0"=bad) | //| RMS training set RMS error | //| CVRMS leave-one-out RMS error | //| Note: | //| content of all arrays is changed by subroutine; | //| it doesn't allocate temporaries. | //+------------------------------------------------------------------+ void CBdSS::DSOptimalSplit2Fast(double &a[],int &c[],int &tiesbuf[], int &cntbuf[],double &bufr[],int &bufi[], const int n,const int nc,double alpha, int &info,double &threshold, double &rms,double &cvrms) { CRowDouble A=a; CRowInt C=c; CRowInt TiesBuf=tiesbuf; CRowInt CntBuf=cntbuf; CRowDouble BufR=bufr; CRowInt BufI=bufi; DSOptimalSplit2Fast(A,C,TiesBuf,CntBuf,BufR,BufI,n,nc,alpha,info,threshold,rms,cvrms); A.ToArray(a); C.ToArray(c); TiesBuf.ToArray(tiesbuf); CntBuf.ToArray(cntbuf); BufR.ToArray(bufr); BufI.ToArray(bufi); } //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ void CBdSS::DSOptimalSplit2Fast(CRowDouble &a,CRowInt &c,CRowInt &tiesbuf, CRowInt &cntbuf,CRowDouble &bufr,CRowInt &bufi, const int n,const int nc,double alpha, int &info,double &threshold, double &rms,double &cvrms) { //--- create variables int i=0; int k=0; int cl=0; int tiecount=0; double cbest=0; double cc=0; int koptimal=0; int sl=0; int sr=0; double v=0; double w=0; double x=0; //--- initialization info=0; threshold=0; rms=0; cvrms=0; //--- Test for errors in inputs if(n<=0 || nc<2) { info=-1; return; } for(i=0; i=nc) { info=-2; return; } } //--- change value info=1; //--- Tie DSTieFastI(a,c,n,tiesbuf,tiecount,bufr,bufi); //--- Special case: number of ties is 1. if(tiecount==1) { info=-3; return; } //--- General case,number of ties > 1 cntbuf.Fill(0,0,2*nc); for(i=0; i1) { w=cntbuf[i]; cvrms+=w*CMath::Sqr((w-1.0)/(sl-1.0)-1.0); cvrms+=(sl-w)*CMath::Sqr(w/(sl-1.0)); } else { w=cntbuf[i]; cvrms+=w*CMath::Sqr(1.0/(double)nc-1.0); cvrms+=(sl-w)*CMath::Sqr(1.0/(double)nc); } //--- check if(sr>1) { w=cntbuf[nc+i]; cvrms+=w*CMath::Sqr((w-1.0)/(sr-1.0)-1.0); cvrms+=(sr-w)*CMath::Sqr(w/(sr-1.0)); } else { w=cntbuf[nc+i]; cvrms+=w*CMath::Sqr(1.0/(double)nc-1.0); cvrms+=(sr-w)*CMath::Sqr(1.0/(double)nc); } } //--- change value cvrms=MathSqrt(cvrms/(double)(nc*n)); } } //--- Calculate threshold. //--- Code is a bit complicated because there can be such //--- numbers that 0.5(A+B) equals to A or B (if A-B=epsilon) threshold=0.5*(a[tiesbuf[koptimal]]+a[tiesbuf[koptimal+1]]); //--- check if(threshold<=a[tiesbuf[koptimal]]) threshold=a[tiesbuf[koptimal+1]]; } //+------------------------------------------------------------------+ //| Automatic non-optimal discretization, internal subroutine. | //+------------------------------------------------------------------+ void CBdSS::DSSplitK(double &ca[],int &cc[],const int n,const int nc, int kmax,int &info,double &thresholds[],int &ni, double &cve) { CRowDouble A=ca; CRowInt C=cc; CRowDouble Thresholds=thresholds; DSSplitK(A,C,n,nc,kmax,info,Thresholds,ni,cve); Thresholds.ToArray(thresholds); } //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ void CBdSS::DSSplitK(CRowDouble &ca,CRowInt &cc,const int n,const int nc, int kmax,int &info,CRowDouble &thresholds,int &ni, double &cve) { //--- create variables int i=0; int j=0; int j1=0; int k=0; int tiecount=0; double v2=0; int bestk=0; double bestcve=0; double curcve=0; //--- creating arrays CRowInt ties; CRowInt p1; CRowInt p2; CRowInt cnt; CRowInt bestsizes; CRowInt cursizes; CRowDouble a=ca; CRowInt c=cc; a.Resize(n); c.Resize(n); //--- initialization info=0; ni=0; cve=0; //--- Test for errors in inputs if((n<=0 || nc<2) || kmax<2) { info=-1; return; } for(i=0; i=nc) { info=-2; return; } } //--- change value info=1; //--- Tie DSTie(a,n,ties,tiecount,p1,p2); //--- swap for(i=0; i0,__FUNCTION__+": internal error #1!")) return; //--- change values bestk=2; bestsizes.Set(0,ties[j]); bestsizes.Set(1,n-j); bestcve=0; //--- calculation cnt.Fill(0); for(i=0; i=nc) { info=-2; return; } } //--- change value info=1; //--- Tie DSTie(a,n,ties,tiecount,p1,p2); //--- swap for(i=0; i=0,__FUNCTION__+": internal error #1!")) return; //--- check if(koptimal==0) { //--- Special case: best partition is one big interval. //--- Even 2-partition is not better. //--- This is possible when dealing with "weak" predictor variables. //--- Make binary split as close to the median as possible. v2=CMath::m_maxrealnumber; j=-1; for(i=1; i<=tiecount-1; i++) { //--- check if(MathAbs(ties[i]-0.5*(n-1))0,__FUNCTION__+": internal error #2!")) return; //--- allocation thresholds.Resize(1); //--- change values thresholds.Set(0,0.5*(a[ties[j-1]]+a[ties[j]])); ni=2; cve=0; //--- calculation cnt.Fill(0); for(i=0; i=1; k--) { thresholds.Set(k-1,0.5*(a[ties[jl-1]]+a[ties[jl]])); jr=jl-1; jl=splits.Get(k-1,jl-1); } } } //+------------------------------------------------------------------+ //| Internal function | //+------------------------------------------------------------------+ double CBdSS::XLnY(const double x,const double y) { //--- check if(x==0.0) return(0.0); //--- return result return(x*MathLog(y)); } //+------------------------------------------------------------------+ //| Internal function, | //| returns number of samples of class I in Cnt[I] | //+------------------------------------------------------------------+ double CBdSS::GetCV(CRowInt &cnt,const int nc) { //--- create variables double result=0; int i=0; double s=0; //--- calculation s=0; for(i=0; i=1 | //| NVars - number of independent variables, NVars>=1 | //| NClasses - indicates type of the problem being solved: | //| * NClasses>=2 means that classification problem is | //| solved (last column of the dataset stores class | //| number) | //| * NClasses=1 means that regression problem is | //| solved (last column of the dataset stores | //| variable value) | //| OUTPUT PARAMETERS: | //| S - decision forest builder | //+------------------------------------------------------------------+ void CDForest::DFBuilderSetDataset(CDecisionForestBuilder &s, CMatrixDouble &xy, int npoints, int nvars, int nclasses) { //--- create variables int i=0; int j=0; //--- Check parameters if(!CAp::Assert(npoints>=1,__FUNCTION__": npoints<1")) return; if(!CAp::Assert(nvars>=1,__FUNCTION__": nvars<1")) return; if(!CAp::Assert(nclasses>=1,__FUNCTION__": nclasses<1")) return; if(!CAp::Assert(xy.Rows()>=npoints,__FUNCTION__": rows(xy)nvars,__FUNCTION__": cols(xy)1) for(i=0; i=0 && j1) { CApServ::IVectorSetLengthAtLeast(s.m_DSIVal,npoints); for(i=0; i=1, number of trees to train | //| OUTPUT PARAMETERS: | //| D - decision forest. You can compress this forest to | //| more compact 16-bit representation with | //| DFBinaryCompression() | //| Rep - report, see below for information on its fields. | //| == report information produced by forest construction function = | //| Decision forest training report includes following information: | //| * training set errors | //| * out-of-bag estimates of errors | //| * variable importance ratings | //| Following fields are used to store information: | //| * training set errors are stored in rep.RelCLSError, rep.AvgCE,| //| rep.RMSError, rep.AvgError and rep.AvgRelError | //| * out-of-bag estimates of errors are stored in | //| rep.oobrelclserror, rep.oobavgce, rep.oobrmserror, | //| rep.oobavgerror and rep.oobavgrelerror | //| Variable importance reports, if requested by | //| DFBuilderSetImportanceGini(), DFBuilderSetImportanceTrnGini() or | //| DFBuilderSetImportancePermutation() call, are stored in: | //| * rep.varimportances field stores importance ratings | //| * rep.topvars stores variable indexes ordered from the most | //| important to less important ones | //| You can find more information about report fields in: | //| * comments on CDFReport structure | //| * comments on DFBuilderSetImportanceGini function | //| * comments on DFBuilderSetImportanceTrnGini function | //| * comments on DFBuilderDetImportancePermutation function | //+------------------------------------------------------------------+ void CDForest::DFBuilderBuildRandomForest(CDecisionForestBuilder &s, int ntrees, CDecisionForest &df, CDFReport &rep) { //--- create variables int nvars=0; int nclasses=0; int npoints=0; int trnsize=0; int maxtreesize=0; int sessionseed=0; CDFVoteBuf buf; //--- check if(!CAp::Assert(ntrees>=1,__FUNCTION__": ntrees<1")) return; CleanReport(s,rep); npoints=s.m_NPoints; nvars=s.m_NVars; nclasses=s.m_NClasses; //--- Set up progress counter s.m_RDFProgress=0; s.m_RDFTotal=ntrees*npoints; if(s.m_RDFImportance==m_NeedPermutation) s.m_RDFTotal+=ntrees*npoints; //--- Quick exit for empty dataset if(s.m_DSType==-1 || npoints==0) { if(!CAp::Assert(m_LeafNodeWdth==2,__FUNCTION__": integrity check failed")) return; df.m_ForestFormat=m_DFUncompressedV0; df.m_NVars=s.m_NVars; df.m_NClasses=s.m_NClasses; df.m_NTrees=1; df.m_BufSize=1+m_LeafNodeWdth; df.m_Trees.Resize(1+m_LeafNodeWdth); df.m_Trees.Set(0,1+m_LeafNodeWdth); df.m_Trees.Set(1,-1.0); df.m_Trees.Set(2,0.0); DFCreateBuffer(df,df.m_Buffer); return; } if(!CAp::Assert(npoints>0,__FUNCTION__": integrity check failed")) return; //--- Analyze dataset statistics, perform preprocessing AnalyzeAndPreprocessDataset(s); //--- Prepare "work", "vote" and "tree" pools and other settings trnsize=(int)MathRound(npoints*s.m_RDFRatio); trnsize=MathMax(trnsize,1); trnsize=MathMin(trnsize,npoints); maxtreesize=1+m_InnerNodeWidth*(trnsize-1)+m_LeafNodeWdth*trnsize; //--- Allocation ResetLastError(); if(!CAp::Assert(ArrayResize(s.m_WorkBuf,ntrees)==ntrees,StringFormat("%s: resize buffer failed (%d)",__FUNCTION__,GetLastError()))) return; if(!CAp::Assert(ArrayResize(s.m_VoteBuf,ntrees)==ntrees,StringFormat("%s: resize buffer failed (%d)",__FUNCTION__,GetLastError()))) return; if(!CAp::Assert(ArrayResize(s.m_TreeBuf,ntrees)==ntrees,StringFormat("%s: resize buffer failed (%d)",__FUNCTION__,GetLastError()))) return; s.m_WorkBuf[0].m_varpool.Resize(nvars); s.m_WorkBuf[0].m_trnset.Resize(trnsize); s.m_WorkBuf[0].m_oobset.Resize(npoints-trnsize); s.m_WorkBuf[0].m_tmp0i.Resize(npoints); s.m_WorkBuf[0].m_tmp1i.Resize(npoints); s.m_WorkBuf[0].m_tmp0r.Resize(npoints); s.m_WorkBuf[0].m_tmp1r.Resize(npoints); s.m_WorkBuf[0].m_tmp2r.Resize(npoints); s.m_WorkBuf[0].m_tmp3r.Resize(npoints); s.m_WorkBuf[0].m_trnlabelsi.Resize(npoints); s.m_WorkBuf[0].m_trnlabelsr.Resize(npoints); s.m_WorkBuf[0].m_ooblabelsi.Resize(npoints); s.m_WorkBuf[0].m_ooblabelsr.Resize(npoints); s.m_WorkBuf[0].m_curvals.Resize(npoints); s.m_WorkBuf[0].m_bestvals.Resize(npoints); s.m_WorkBuf[0].m_classpriors.Resize(nclasses); s.m_WorkBuf[0].m_classtotals0.Resize(nclasses); s.m_WorkBuf[0].m_classtotals1.Resize(nclasses); s.m_WorkBuf[0].m_classtotals01.Resize(2*nclasses); s.m_WorkBuf[0].m_treebuf.Resize(maxtreesize); s.m_WorkBuf[0].m_trnsize=trnsize; s.m_WorkBuf[0].m_oobsize=npoints-trnsize; s.m_VoteBuf[0].m_trntotals=vector::Zeros(npoints*nclasses); s.m_VoteBuf[0].m_oobtotals=s.m_VoteBuf[0].m_trntotals; s.m_VoteBuf[0].m_trncounts.Resize(npoints); s.m_VoteBuf[0].m_oobcounts.Resize(npoints); s.m_VoteBuf[0].m_trncounts.Fill(0); s.m_VoteBuf[0].m_oobcounts.Fill(0); s.m_VoteBuf[0].m_giniimportances=vector::Zeros(nvars); for(int treeIdx=1; treeIdx1) { //--- classification-specific code int k=(int)MathRound(xy.Get(i,df.m_NVars)); //--- check if(y[k]!=0.0) result-=MathLog(y[k]); else result-=MathLog(CMath::m_minrealnumber); } } //--- return result return(result/npoints); } //+------------------------------------------------------------------+ //| RMS error on the test set | //| INPUT PARAMETERS: | //| DF - decision forest model | //| XY - test set | //| NPoints - test set size | //| RESULT: | //| root mean square error. | //| Its meaning for regression task is obvious. As for | //| classification task,RMS error means error when estimating | //| posterior probabilities. | //+------------------------------------------------------------------+ double CDForest::DFRMSError(CDecisionForest &df,CMatrixDouble &xy, const int npoints) { //--- create variables double result=0; CRowDouble x; CRowDouble y; for(int i=0; i1) { //--- classification-specific code int k=(int)MathRound(xy.Get(i,df.m_NVars)); y.Add(k,-1); result+=MathPow(y.ToVector()+0,2.0).Sum(); } else { //--- regression-specific code result+=CMath::Sqr(y[0]-xy.Get(i,df.m_NVars)); } } //--- return result return(MathSqrt(result/(npoints*df.m_NClasses))); } //+------------------------------------------------------------------+ //| Average error on the test set | //| INPUT PARAMETERS: | //| DF - decision forest model | //| XY - test set | //| NPoints - test set size | //| RESULT: | //| Its meaning for regression task is obvious. As for | //| classification task, it means average error when estimating | //| posterior probabilities. | //+------------------------------------------------------------------+ double CDForest::DFAvgError(CDecisionForest &df,CMatrixDouble &xy, const int npoints) { //--- create variables double result=0; CRowDouble x; CRowDouble y; for(int i=0; i1) { //--- classification-specific code int k=(int)MathRound(xy.Get(i,df.m_NVars)); y.Add(k,-1); result+=(y.Abs()+0).Sum(); } else { //--- regression-specific code result+=MathAbs(y[0]-xy.Get(i,df.m_NVars)); } } //--- return result return(result/(npoints*df.m_NClasses)); } //+------------------------------------------------------------------+ //| Average relative error on the test set | //| INPUT PARAMETERS: | //| DF - decision forest model | //| XY - test set | //| NPoints - test set size | //| RESULT: | //| Its meaning for regression task is obvious. As for | //| classification task, it means average relative error when | //| estimating posterior probability of belonging to the correct | //| class. | //+------------------------------------------------------------------+ double CDForest::DFAvgRelError(CDecisionForest &df,CMatrixDouble &xy, const int npoints) { //--- create variables double result=0; int relcnt=0; //--- creating arrays CRowDouble x; CRowDouble y; //--- initialization for(int i=0; i1) { //--- classification-specific code int k=(int)MathRound(xy.Get(i,df.m_NVars)); if(k<0 || k>=df.m_NVars) continue; result+=MathAbs(y[k]-1); relcnt++; } else { //--- regression-specific code if(xy.Get(i,df.m_NVars)!=0.0) { result+=MathAbs((y[0]-xy.Get(i,df.m_NVars))/xy.Get(i,df.m_NVars)); relcnt ++; } } } //--- check if(relcnt>0) result=result/relcnt; //--- return result return(result); } //+------------------------------------------------------------------+ //| Copying of DecisionForest strucure | //| INPUT PARAMETERS: | //| DF1 - original | //| OUTPUT PARAMETERS: | //| DF2 - copy | //+------------------------------------------------------------------+ void CDForest::DFCopy(CDecisionForest &df1,CDecisionForest &df2) { if(df1.m_ForestFormat==m_DFUncompressedV0) { df2.m_ForestFormat=df1.m_ForestFormat; df2.m_NVars=df1.m_NVars; df2.m_NClasses=df1.m_NClasses; df2.m_NTrees=df1.m_NTrees; df2.m_BufSize=df1.m_BufSize ; df2.m_Trees=df1.m_Trees; DFCreateBuffer(df2,df2.m_Buffer); return; } if(df1.m_ForestFormat==m_DFCompressedV0) { df2.m_ForestFormat=df1.m_ForestFormat; df2.m_NVars=df1.m_NVars; df2.m_NClasses=df1.m_NClasses; df2.m_NTrees=df1.m_NTrees; df2.m_BufSize=df1.m_BufSize ; df2.m_UseMantissa8=df1.m_UseMantissa8; df2.m_Trees8=df1.m_Trees8; DFCreateBuffer(df2,df2.m_Buffer); return; } CAp::Assert(false,__FUNCTION__": unexpected forest format"); } //+------------------------------------------------------------------+ //| Serializer: allocation | //+------------------------------------------------------------------+ void CDForest::DFAlloc(CSerializer &s,CDecisionForest &forest) { //--- preparation to serialize if(forest.m_ForestFormat==m_DFUncompressedV0) { s.Alloc_Entry(); s.Alloc_Entry(); s.Alloc_Entry(); s.Alloc_Entry(); s.Alloc_Entry(); s.Alloc_Entry(); CApServ::AllocRealArray(s,forest.m_Trees,forest.m_BufSize); return; } if(forest.m_ForestFormat==m_DFCompressedV0) { s.Alloc_Entry(); s.Alloc_Entry(); s.Alloc_Entry(); s.Alloc_Entry(); s.Alloc_Entry(); s.Alloc_Entry(); CApServ::AllocIntegerArray(s,forest.m_Trees8,forest.m_BufSize); return; } CAp::Assert(false,__FUNCTION__": unexpected forest format"); } //+------------------------------------------------------------------+ //| Serializer: serialization | //+------------------------------------------------------------------+ void CDForest::DFSerialize(CSerializer &s,CDecisionForest &forest) { //--- serializetion if(forest.m_ForestFormat==m_DFUncompressedV0) { s.Serialize_Int(CSCodes::GetRDFSerializationCode()); s.Serialize_Int(m_DFUncompressedV0); s.Serialize_Int(forest.m_NVars); s.Serialize_Int(forest.m_NClasses); s.Serialize_Int(forest.m_NTrees); s.Serialize_Int(forest.m_BufSize); CApServ::SerializeRealArray(s,forest.m_Trees,forest.m_BufSize); return; } if(forest.m_ForestFormat==m_DFCompressedV0) { s.Serialize_Int(CSCodes::GetRDFSerializationCode()); s.Serialize_Int(forest.m_ForestFormat); s.Serialize_Bool(forest.m_UseMantissa8); s.Serialize_Int(forest.m_NVars); s.Serialize_Int(forest.m_NClasses); s.Serialize_Int(forest.m_NTrees); CApServ::SerializeIntegerArray(s,forest.m_Trees8,forest.m_Trees8.Size()); return; } CAp::Assert(false,__FUNCTION__": unexpected forest format"); } //+------------------------------------------------------------------+ //| Serializer: unserialization | //+------------------------------------------------------------------+ void CDForest::DFUnserialize(CSerializer &s,CDecisionForest &forest) { bool processed=false; //--- check correctness of header int i0=s.Unserialize_Int(); if(!CAp::Assert(i0==CSCodes::GetRDFSerializationCode(),__FUNCTION__": stream header corrupted")) return; //--- Read forest int forestformat=s.Unserialize_Int(); if(forestformat==m_DFUncompressedV0) { //--- Unserialize data forest.m_ForestFormat=forestformat; forest.m_NVars=s.Unserialize_Int(); forest.m_NClasses=s.Unserialize_Int(); forest.m_NTrees=s.Unserialize_Int(); forest.m_BufSize=s.Unserialize_Int(); CApServ::UnserializeRealArray(s,forest.m_Trees); processed=true; } if(forestformat==m_DFCompressedV0) { //--- Unserialize data forest.m_ForestFormat=forestformat; forest.m_UseMantissa8=s.Unserialize_Bool(); forest.m_NVars=s.Unserialize_Int(); forest.m_NClasses=s.Unserialize_Int(); forest.m_NTrees=s.Unserialize_Int(); CApServ::UnserializeIntegerArray(s,forest.m_Trees8); processed=true; } if(!CAp::Assert(processed,__FUNCTION__": unexpected forest format")) return; //--- Prepare buffer DFCreateBuffer(forest,forest.m_Buffer); } //+------------------------------------------------------------------+ //| This subroutine builds random decision forest. | //| ---- DEPRECATED VERSION! USE DECISION FOREST BUILDER OBJECT ---- | //+------------------------------------------------------------------+ void CDForest::DFBuildRandomDecisionForest(CMatrixDouble &xy, int npoints, int nvars, int nclasses, int ntrees, double r, int &info, CDecisionForest &df, CDFReport &rep) { info=0; //--- check if(r<=0.0 || r>1.0) { info=-1; return; } int samplesize=MathMax((int)MathRound(r*npoints),1); DFBuildInternal(xy,npoints,nvars,nclasses,ntrees,samplesize,MathMax(nvars/2,1),m_DFUseStrongSplits+m_DFUseEVS,info,df,rep); } //+------------------------------------------------------------------+ //| This subroutine builds random decision forest. | //| ---- DEPRECATED VERSION! USE DECISION FOREST BUILDER OBJECT ---- | //+------------------------------------------------------------------+ void CDForest::DFBuildRandomDecisionForestX1(CMatrixDouble &xy, int npoints, int nvars, int nclasses, int ntrees, int nrndvars, double r, int &info, CDecisionForest &df, CDFReport &rep) { info=0; //--- check if(r<=0.0 || r>1.0) { info=-1; return; } if(nrndvars<=0 || nrndvars>nvars) { info=-1; return; } int samplesize=MathMax((int)MathRound(r*npoints),1); DFBuildInternal(xy,npoints,nvars,nclasses,ntrees,samplesize,nrndvars,m_DFUseStrongSplits+m_DFUseEVS,info,df,rep); } //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ void CDForest::DFBuildInternal(CMatrixDouble &xy, int npoints, int nvars, int nclasses, int ntrees, int samplesize, int nfeatures, int flags, int &info, CDecisionForest &df, CDFReport &rep) { CDecisionForestBuilder builder; info=0; //--- Test for inputs if(npoints<1 || samplesize<1 || samplesize>npoints || nvars<1 || nclasses<1 || ntrees<1 || nfeatures<1) { info=-1; return; } if(nclasses>1) for(int i=0; i<=npoints-1; i++) if(MathRound(xy.Get(i,nvars))<0 || MathRound(xy.Get(i,nvars))>=nclasses) { info=-2; return; } info=1; DFBuilderCreate(builder); DFBuilderSetDataset(builder,xy,npoints,nvars,nclasses); DFBuilderSetSubsampleRatio(builder,(double)samplesize/(double)npoints); DFBuilderSetRndVars(builder,nfeatures); DFBuilderBuildRandomForest(builder,ntrees,df,rep); } //+------------------------------------------------------------------+ //| Builds a range of random trees [TreeIdx0, TreeIdx1) using | //| decision forest algorithm. Tree index is used to seed per - tree | //| RNG. | //+------------------------------------------------------------------+ void CDForest::BuildRandomTree(CDecisionForestBuilder &s,int treeidx0, int treeidx1) { //--- create variables int i=0; int j=0; int npoints=0; int nvars=0; int nclasses=0; CHighQualityRandState rs; int treesize=0; int varstoselect=0; int workingsetsize=0; double meanloss=0; //--- Prepare npoints=s.m_NPoints; nvars=s.m_NVars; nclasses=s.m_NClasses; //--- if(treeidx0>treeidx1) CApServ::Swap(treeidx0,treeidx1); if(treeidx0==treeidx1) treeidx1++; for(int treeidx=treeidx0; treeidx0) { CHighQualityRand::HQRndSeed(s.m_RDFGlobalSeed,1+treeidx,rs); } else { CHighQualityRand::HQRndSeed(CMath::RandomInteger(30000),1+treeidx,rs); } //--- Prepare everything for tree construction. if(!CAp::Assert(s.m_WorkBuf[treeidx].m_trnsize>=1,__FUNCTION__": integrity check failed (34636)")) return; if(!CAp::Assert(s.m_WorkBuf[treeidx].m_oobsize>=0,__FUNCTION__": integrity check failed (45745)")) return; if(!CAp::Assert(s.m_WorkBuf[treeidx].m_trnsize+s.m_WorkBuf[treeidx].m_oobsize==npoints,__FUNCTION__": integrity check failed (89415)")) return; workingsetsize=-1; s.m_WorkBuf[treeidx].m_varpoolsize=0; for(i=0; i=0,__FUNCTION__": integrity check failed (73f5)")) return; for(i=0; i1) s.m_WorkBuf[treeidx].m_trnlabelsi.Set(i,s.m_DSIVal[s.m_WorkBuf[treeidx].m_tmp0i[i]]); else s.m_WorkBuf[treeidx].m_trnlabelsr.Set(i,s.m_DSRVal[s.m_WorkBuf[treeidx].m_tmp0i[i]]); if(s.m_NeedIOBMatrix) s.m_IOBMatrix.Set(treeidx,s.m_WorkBuf[treeidx].m_trnset[i],(int)true); } for(i=0; i1) s.m_WorkBuf[treeidx].m_ooblabelsi.Set(i,s.m_DSIVal[j]); else s.m_WorkBuf[treeidx].m_ooblabelsr.Set(i,s.m_DSRVal[j]); } varstoselect=(int)MathRound(MathSqrt(nvars)); if(s.m_RDFVars>0.0) varstoselect=(int)MathRound(s.m_RDFVars); else if(s.m_RDFVars<0.0) { varstoselect=(int)MathRound(-(nvars*s.m_RDFVars)); } varstoselect=MathMax(varstoselect,1); varstoselect=MathMin(varstoselect,nvars); //--- Perform recurrent construction if(s.m_RDFImportance==m_NeedTrnGini) meanloss=MeanNRMS2(nclasses,s.m_WorkBuf[treeidx].m_trnlabelsi,s.m_WorkBuf[treeidx].m_trnlabelsr,0,s.m_WorkBuf[treeidx].m_trnsize,s.m_WorkBuf[treeidx].m_trnlabelsi,s.m_WorkBuf[treeidx].m_trnlabelsr,0,s.m_WorkBuf[treeidx].m_trnsize,s.m_WorkBuf[treeidx].m_tmpnrms2); else meanloss=MeanNRMS2(nclasses,s.m_WorkBuf[treeidx].m_trnlabelsi,s.m_WorkBuf[treeidx].m_trnlabelsr,0,s.m_WorkBuf[treeidx].m_trnsize,s.m_WorkBuf[treeidx].m_ooblabelsi,s.m_WorkBuf[treeidx].m_ooblabelsr,0,s.m_WorkBuf[treeidx].m_oobsize,s.m_WorkBuf[treeidx].m_tmpnrms2); treesize=1; BuildRandomTreeRec(s,s.m_WorkBuf[treeidx],workingsetsize,varstoselect,s.m_WorkBuf[treeidx].m_treebuf,s.m_VoteBuf[treeidx],rs,0,s.m_WorkBuf[treeidx].m_trnsize,0,s.m_WorkBuf[treeidx].m_oobsize,meanloss,meanloss,treesize); s.m_WorkBuf[treeidx].m_treebuf.Set(0,treesize); //--- Store tree s.m_TreeBuf[treeidx].m_treebuf=s.m_WorkBuf[treeidx].m_treebuf; s.m_TreeBuf[treeidx].m_treeidx=treeidx; } } //+------------------------------------------------------------------+ //| Recurrent tree construction function using caller - allocated | //| buffers and caller - initialized RNG. | //| Following iterms are processed: | //| * items [Idx0, Idx1) of WorkBuf.TrnSet | //| * items [OOBIdx0, OOBIdx1) of WorkBuf.OOBSet | //| TreeSize on input must be 1(header element of the tree), on | //| output it contains size of the tree. | //| OOBLoss on input must contain value of MeanNRMS2(...) computed | //| for entire dataset. | //| Variables from #0 to #WorkingSet - 1 from WorkBuf.VarPool are | //| used(for block algorithm: blocks, not vars) | //+------------------------------------------------------------------+ void CDForest::BuildRandomTreeRec(CDecisionForestBuilder &s, CDFWorkBuf &workbuf, int workingset, int varstoselect, CRowDouble &treebuf, CDFVoteBuf &votebuf, CHighQualityRandState &rs, int idx0, int idx1, int oobidx0, int oobidx1, double meanloss, double topmostmeanloss, int &treesize) { //---create variables int npoints=0; int nclasses=0; int i=0; int j=0; int j0=0; double v=0; bool labelsaresame; int offs=0; int varbest=0; double splitbest=0; int i1=0; int i2=0; int idxtrn=0; int idxoob=0; double meanloss0=0; double meanloss1=0; //--- check if(!CAp::Assert(s.m_DSType==0,__FUNCTION__": not supported skbdgfsi!")) return; if(!CAp::Assert(idx01) { labelsaresame=true; workbuf.m_classpriors.Fill(0); j0=workbuf.m_trnlabelsi[idx0]; for(i=idx0; i1) v=workbuf.m_trnlabelsi[idx0+CHighQualityRand::HQRndUniformI(rs,idx1-idx0)]; else v=workbuf.m_trnlabelsr[idx0+CHighQualityRand::HQRndUniformI(rs,idx1-idx0)]; OutputLeaf(s,workbuf,treebuf,votebuf,idx0,idx1,oobidx0,oobidx1,treesize,v); return; } //--- Good split WAS found, we can perform it: //--- * first, we split training set //--- * then, we similarly split OOB set if(!CAp::Assert(s.m_DSType==0,__FUNCTION__": not supported 54bfdh")) return; offs=npoints*varbest; i1=idx0; i2=idx1-1; while(i1<=i2) { //--- Reorder indexes so that left partition is in [Idx0..I1), //--- and right partition is in [I2+1..Idx1) if(workbuf.m_bestvals[i1]=splitbest) { i2--; continue; } workbuf.m_trnset.Swap(i1,i2); if(nclasses>1) workbuf.m_trnlabelsi.Swap(i1,i2); else workbuf.m_trnlabelsr.Swap(i1,i2); i1++; i2--; } if(!CAp::Assert(i1==i2+1,__FUNCTION__": integrity check failed (45rds3)")) return; idxtrn=i1; if(oobidx0=splitbest) { i2--; continue; } workbuf.m_oobset.Swap(i1,i2); if(nclasses>1) workbuf.m_ooblabelsi.Swap(i1,i2); else workbuf.m_ooblabelsr.Swap(i1,i2); i1++; i2--; } if(!CAp::Assert(i1==i2+1,__FUNCTION__": integrity check failed (643fs3)")) return; idxoob=i1; } else idxoob=oobidx0; //--- Compute estimates of NRMS2 loss over TRN or OOB subsets, update Gini importances if(s.m_RDFImportance==m_NeedTrnGini) { meanloss0=MeanNRMS2(nclasses,workbuf.m_trnlabelsi,workbuf.m_trnlabelsr,idx0,idxtrn,workbuf.m_trnlabelsi,workbuf.m_trnlabelsr,idx0,idxtrn,workbuf.m_tmpnrms2); meanloss1=MeanNRMS2(nclasses,workbuf.m_trnlabelsi,workbuf.m_trnlabelsr,idxtrn,idx1,workbuf.m_trnlabelsi,workbuf.m_trnlabelsr,idxtrn,idx1,workbuf.m_tmpnrms2); } else { meanloss0=MeanNRMS2(nclasses,workbuf.m_trnlabelsi,workbuf.m_trnlabelsr,idx0,idxtrn,workbuf.m_ooblabelsi,workbuf.m_ooblabelsr,oobidx0,idxoob,workbuf.m_tmpnrms2); meanloss1=MeanNRMS2(nclasses,workbuf.m_trnlabelsi,workbuf.m_trnlabelsr,idxtrn,idx1,workbuf.m_ooblabelsi,workbuf.m_ooblabelsr,idxoob,oobidx1,workbuf.m_tmpnrms2); } votebuf.m_giniimportances.Add(varbest,(meanloss-(meanloss0+meanloss1))/(topmostmeanloss+1.0e-20)); //--- Generate tree node and subtrees (recursively) treebuf.Set(treesize,varbest); treebuf.Set(treesize+1,splitbest); i=treesize; treesize=treesize+m_InnerNodeWidth; BuildRandomTreeRec(s,workbuf,workingset,varstoselect,treebuf,votebuf,rs,idx0,idxtrn,oobidx0,idxoob,meanloss0,topmostmeanloss,treesize); treebuf.Set(i+2,treesize); BuildRandomTreeRec(s,workbuf,workingset,varstoselect,treebuf,votebuf,rs,idxtrn,idx1,idxoob,oobidx1,meanloss1,topmostmeanloss,treesize); } //+------------------------------------------------------------------+ //| Estimates permutation variable importance ratings for a range of | //| dataset points. | //| Initial call to this function should span entire range of the | //| dataset, [Idx0, Idx1) = [0, NPoints), because function performs | //| initialization of some internal structures when called with these| //| arguments. | //+------------------------------------------------------------------+ void CDForest::EstimateVariableImportance(CDecisionForestBuilder &s, int sessionseed, CDecisionForest &df, int ntrees, CDFReport &rep) { //--- create variables int npoints=0; int nvars=0; int nclasses=0; int nperm=0; int i=0; int j=0; int k=0; CRowDouble tmpr0; CRowDouble tmpr1; CRowInt tmpi0; vector losses; CDFPermimpBuf permseed; double nopermloss=0; double totalpermloss=0; CHighQualityRandState varimprs; npoints=s.m_NPoints; nvars=s.m_NVars; nclasses=s.m_NClasses; //--- No importance rating if(s.m_RDFImportance==0) return; //--- Gini importance if(s.m_RDFImportance==m_NeedTrnGini || s.m_RDFImportance==m_NeedOOBGini) { //--- Merge OOB Gini importances computed during tree generation int total=ArraySize(s.m_VoteBuf); losses=s.m_VoteBuf[0].m_giniimportances.ToVector(); for(int vote=1; vote=ntrees && s.m_IOBMatrix.Cols()>=npoints,__FUNCTION__": integrity check failed (IOB)")) return; //--- Generate packed representation of the shuffle which is applied to all variables //--- Ideally we want to apply different permutations to different variables, //--- i.e. we have to generate and store NPoints*NVars random numbers. //--- However due to performance and memory restrictions we prefer to use compact //--- representation: //---*we store one "reference" permutation P_ref in VarImpShuffle2[0:NPoints-1] //--- * a permutation P_j applied to variable J is obtained by circularly shifting //--- elements in P_ref by VarImpShuffle2[NPoints+J] CHighQualityRand::HQRndSeed(sessionseed,1117,varimprs); s.m_VarImpShuffle2.Resize(npoints+nvars); for(i=0; i::Zeros(nperm); permseed.m_yv.Resize(nperm*nclasses); permseed.m_xraw.Resize(nvars); permseed.m_xdist.Resize(nvars); permseed.m_xcur.Resize(nvars); permseed.m_targety.Resize(nclasses); permseed.m_startnodes.Resize(nvars); permseed.m_y.Resize(nclasses); //--- Recursively split subset EstimatePermutationImportances(s,df,ntrees,permseed,0,npoints); //--- Merge results losses=permseed.m_losses+vector::Full(nperm,1.0e-20); //--- Compute importances nopermloss=losses[nvars+1]; totalpermloss=losses[nvars]; losses= losses/totalpermloss-nopermloss/totalpermloss; losses.Clip(0,1); losses.Resize(nvars); rep.m_varimportances=losses; //--- Compute topvars[] array tmpr0=losses*(-1); rep.m_topvars.Resize(nvars); for(j=0; j=0 && idx0<=idx1 && idx1<=npoints,__FUNCTION__": integrity check failed (idx)")) return; if(!CAp::Assert(s.m_IOBMatrix.Rows()>=ntrees && s.m_IOBMatrix.Cols()>=npoints,__FUNCTION__": integrity check failed (IOB)")) return; //--- main loop nperm=nvars+2; //--- Process range of points [idx0,idx1) for(i=idx0; i1) { permimpbuf.m_targety=vector::Zeros(nclasses); permimpbuf.m_targety.Set(s.m_DSIVal[i],1); } else permimpbuf.m_targety.Set(0,s.m_DSRVal[i]); //--- Process all trees, for each tree compute NPerm losses corresponding //--- to various permutations of variable values permimpbuf.m_yv=vector::Zeros(nperm*nclasses); oobcounts=0; treeroot=0; for(k=0; k1) { j=(int)MathRound(prediction); permimpbuf.m_yv.Add(varidx*nclasses+j,1); } else permimpbuf.m_yv.Add(varidx,prediction); //--- Save loss for all variables being perturbed (XDist). //--- This loss is used as a reference loss when we compute R-squared. varidx=nvars; permimpbuf.m_y=vector::Zeros(nclasses); DFProcessInternalUncompressed(df,treeroot,treeroot+1,permimpbuf.m_xdist,permimpbuf.m_y); permimpbuf.m_yv.Add(varidx*nclasses+j,permimpbuf.m_y[j]); //--- Compute losses for variable #VarIdx being perturbed. Quite an often decision //--- process does not actually depend on the variable #VarIdx (path from the tree //--- root does not include splits on this variable). In such cases we perform //--- quick exit from the loop with precomputed value. permimpbuf.m_xcur=permimpbuf.m_xraw; for(varidx=0; varidx=0) { //--- Path from tree root to the final leaf involves split on variable #VarIdx. //--- Restart computation from the position first split on #VarIdx. if(!CAp::Assert(df.m_ForestFormat==m_DFUncompressedV0,__FUNCTION__": integrity check failed (ff)")) return; permimpbuf.m_xcur.Set(varidx,permimpbuf.m_xdist[varidx]); nodeoffs=permimpbuf.m_startnodes[varidx]; while(true) { if(df.m_Trees[nodeoffs]==-1.0) { if(nclasses>1) { j=(int)MathRound(df.m_Trees[nodeoffs+1]); permimpbuf.m_yv.Add(varidx*nclasses+j,1.0); } else permimpbuf.m_yv.Add(varidx,df.m_Trees[nodeoffs+1]); break; } j=(int)MathRound(df.m_Trees[nodeoffs]); if(permimpbuf.m_xcur[j]1) { j=(int)MathRound(prediction); permimpbuf.m_yv.Add(varidx*nclasses+j,1.0); } else permimpbuf.m_yv.Add(varidx,prediction); } } //--- update OOB counter oobcounts++; } treeroot+=(int)MathRound(df.m_Trees[treeroot]); } //--- Now YV[] stores NPerm versions of the forest output for various permutations of variable values. //--- Update losses. for(j=0; j::Zeros(s.m_NVars); for(int i=0; i<=s.m_NVars-1; i++) rep.m_topvars.Set(i,i); } //+------------------------------------------------------------------+ //| This function returns NRMS2 loss(sum of squared residuals) for a | //| constant - output model: | //| * model output is a mean over TRN set being passed (for | //| classification problems - NClasses - dimensional vector | //| of class probabilities) | //| * model is evaluated over TST set being passed, with L2 loss | //| being returned | //| Input parameters: | //| NClasses - ">1" for classification, "=1" for regression | //| TrnLabelsI - training set labels, class indexes (for | //| NClasses > 1) | //| TrnLabelsR - training set output values (for NClasses = 1) | //| TrnIdx0, | //| TrnIdx1 - a range [Idx0, Idx1) of elements in LabelsI/R | //| is considered | //| TstLabelsI - training set labels, class indexes (for | //| NClasses > 1) | //| TstLabelsR - training set output values (for NClasses = 1) | //| TstIdx0, | //| TstIdx1 - a range [Idx0, Idx1) of elements in LabelsI/R | //| is considered | //| TmpI - temporary array, reallocated as needed | //| Result: sum of squared residuals; for NClasses >= 2 it coincides | //| with Gini impurity times(Idx1 - Idx0) | //| Following fields of WorkBuf are used as temporaries: | //| * TmpMeanNRMS2 | //+------------------------------------------------------------------+ double CDForest::MeanNRMS2(int nclasses,CRowInt &trnlabelsi, CRowDouble &trnlabelsr,int trnidx0, int trnidx1,CRowInt &tstlabelsi, CRowDouble &tstlabelsr,int tstidx0, int tstidx1,CRowInt &tmpi) { //--- create variables double result=0; int i=0; int k=0; int ntrn=trnidx1-trnidx0; int ntst=tstidx1-tstidx0; double v=0; double vv=0; double invntrn=0; double pitrn=0; double nitst=0; //--- check if(!CAp::Assert(trnidx0<=trnidx1,__FUNCTION__": integrity check failed (8754)")) return(result); if(!CAp::Assert(tstidx0<=tstidx1,__FUNCTION__": integrity check failed (8754)")) return(result); //--- Quick exit if(ntrn==0 || ntst==0) return(result); invntrn=1.0/(double)ntrn; if(nclasses>1) { //--- Classification problem tmpi.Resize(2*nclasses); tmpi.Fill(0); for(i=trnidx0; i0 && (varbest<0 || currms<=errbest)) { errbest=currms; varbest=varcur; splitbest=split; for(i=idx0; i1) { //--- Classification problem workbuf.m_classtotals0.Fill(0); sl=0; for(i=idx0; i1) { for(i=0; i=20) advanceby=MathMax(2,(int)MathRound(n*0.05)); info=-1; threshold=0; e=CMath::m_maxrealnumber; //--- Random split if(s.m_RDFSplitStrength==0) { //--- Evaluate minimum, maximum and randomly selected values vmin=x[0]; vmax=x[0]; for(i=1; ivmax) vmax=v; } if(vmin==vmax) return; v=x[CHighQualityRand::HQRndUniformI(rs,n)]; if(v==vmin) v=vmax; //--- Calculate RMS error associated with the split workbuf.m_classtotals0.Fill(0,0,nclasses); n0=0; for(i=0; i0 && n00 && n00) e=MathSqrt(e/(nclasses*n)); break; default: CAp::Assert(false,__FUNCTION__": ClassifierSplit(),critical error"); break; } } //+------------------------------------------------------------------+ //| Regression model split | //+------------------------------------------------------------------+ void CDForest::RegressionSplit(CDecisionForestBuilder &s, CDFWorkBuf &workbuf, CRowDouble &x, CRowDouble &y, int n, int &info, double &threshold, double &e, CRowDouble &sortrbuf, CRowDouble &sortrbuf2) { //--- create variables int i=0; double vmin=0; double vmax=0; double bnd01=0; double bnd12=0; double bnd23=0; int total0=0; int total1=0; int total2=0; int total3=0; int cnt0=0; int cnt1=0; int cnt2=0; int cnt3=0; int n0=0; int advanceby=0; double v=0; double v0=0; double v1=0; double rms=0; int n0prev=0; int k0=0; int k1=0; info=0; threshold=0; e=0; advanceby=1; if(n>=20) advanceby=MathMax(2,(int)MathRound(n*0.05)); //--- Sort data //--- Quick check for degeneracy CTSort::TagSortFastR(x,y,sortrbuf,sortrbuf2,n); v=0.5*(x[0]+x[n-1]); if(!(x[0]vmax) vmax=v; } bnd12=0.5*(vmin+vmax); bnd01=0.5*(vmin+bnd12); bnd23=0.5*(vmax+bnd12); total0=0; total1=0; total2=0; total3=0; for(i=0; i0) e=MathSqrt(e/(4*n)); } //+------------------------------------------------------------------+ //| Returns split: either deterministic split at the middle | //| of [A, B], or randomly chosen split. | //| It is guaranteed that A < Split <= B. | //+------------------------------------------------------------------+ double CDForest::GetSplit(CDecisionForestBuilder &s,double a, double b,CHighQualityRandState &rs) { double result=0.5*(a+b); //--- check if(result<=a) result=b; //--- return result return(result); } //+------------------------------------------------------------------+ //| Outputs leaf to the tree | //| Following items of TRN and OOB sets are updated in the voting | //| buffer: | //| * items [Idx0, Idx1) of WorkBuf.TrnSet | //| * items [OOBIdx0, OOBIdx1) of WorkBuf.OOBSet | //+------------------------------------------------------------------+ void CDForest::OutputLeaf(CDecisionForestBuilder &s, CDFWorkBuf &workbuf, CRowDouble &treebuf, CDFVoteBuf &votebuf, int idx0, int idx1, int oobidx0, int oobidx1, int &treesize, double leafval) { //--- create variables int nclasses=s.m_NClasses; int leafvali=0; int i=0; int j=0; if(nclasses==1) { //--- Store split to the tree treebuf.Set(treesize,-1); treebuf.Set(treesize+1,leafval); //--- Update training and OOB voting stats for(i=idx0; i0) CHighQualityRand::HQRndSeed(s.m_RDFGlobalSeed,3532,rs); else CHighQualityRand::HQRndSeed(CMath::RandomInteger(30000),3532,rs); //--- Generic processing if(!CAp::Assert(npoints>=1,__FUNCTION__": integrity check failed")) return; s.m_DSMin.Resize(nvars); s.m_DSMax.Resize(nvars); CApServ::BVectorSetLengthAtLeast(s.m_DSBinary,nvars); for(i=0; iv1) v1=v; } s.m_DSMin.Set(i,v0); s.m_DSMax.Set(i,v1); if(!CAp::Assert(v0<=v1,__FUNCTION__": strange integrity check failure")) return; isbinary=true; for(j=0; j= 1, number of trees to train | //| OUTPUT PARAMETERS: | //| DF - decision forest | //| Rep - report | //+------------------------------------------------------------------+ void CDForest::MergeTrees(CDecisionForestBuilder &s, CDecisionForest &df) { //--- create variables int i=0; int cursize=0; int offs=0; CRowInt treesizes; CRowInt treeoffsets; //--- initialize df.m_ForestFormat=m_DFUncompressedV0; df.m_NVars=s.m_NVars; df.m_NClasses=s.m_NClasses; df.m_BufSize=0; //--- Determine trees count df.m_NTrees=ArraySize(s.m_TreeBuf); if(!CAp::Assert(df.m_NTrees>0,__FUNCTION__": integrity check failed,zero trees count")) return; //--- Determine individual tree sizes and total buffer size treesizes.Resize(df.m_NTrees); treesizes.Fill(-1); for(int idx=0; idx=0 && s.m_TreeBuf[idx].m_treeidx0,__FUNCTION__": integrity check failed (wrong TreeSize)")) return; df.m_BufSize+=size; treesizes.Set(s.m_TreeBuf[idx].m_treeidx,size); } //--- Determine offsets for individual trees in output buffer treeoffsets.Resize(df.m_NTrees); treeoffsets.Set(0,0); for(i=1; i0,__FUNCTION__": integrity check failed")) return; if(!CAp::Assert(nvars>0,__FUNCTION__": integrity check failed")) return; if(!CAp::Assert(nclasses>0,__FUNCTION__": integrity check failed")) return; //--- Prepare vote buffer buf.m_trntotals=vector::Zeros(npoints*nclasses); buf.m_oobtotals=vector::Zeros(npoints*nclasses); buf.m_trncounts.Resize(npoints); buf.m_oobcounts.Resize(npoints); buf.m_trncounts.Fill(0); buf.m_oobcounts.Fill(0); //--- Merge voting arrays for(int idx=0; idx1) { //--- classification-specific code k=s.m_DSIVal[i]; for(j=0; jbuf.m_trntotals[i*nclasses+k1]) k1=j; } if(k1!=k) rep.m_RelCLSError++; k1=0; for(j=1; jbuf.m_oobtotals[i*nclasses+k1]) k1=j; } if(k1!=k) rep.m_oobrelclserror++; } else { //--- regression-specific code v=buf.m_trntotals[i]-s.m_DSRVal[i]; rep.m_RMSError+=CMath::Sqr(v); rep.m_AvgError+=MathAbs(v); if(s.m_DSRVal[i]!=0.0) { rep.m_AvgRelError+=MathAbs(v/s.m_DSRVal[i]); avgrelcnt++; } v=buf.m_oobtotals[i]-s.m_DSRVal[i]; rep.m_oobrmserror+=CMath::Sqr(v); rep.m_oobavgerror+=MathAbs(v); if(s.m_DSRVal[i]!=0.0) { rep.m_oobavgrelerror+=MathAbs(v/s.m_DSRVal[i]); oobavgrelcnt++; } } } rep.m_RelCLSError/=npoints; rep.m_RMSError=MathSqrt(rep.m_RMSError/(npoints*nclasses)); rep.m_AvgError/=(npoints*nclasses); rep.m_AvgRelError/=CApServ::Coalesce(avgrelcnt,1); rep.m_oobrelclserror/=npoints; rep.m_oobrmserror=MathSqrt(rep.m_oobrmserror/(npoints*nclasses)); rep.m_oobavgerror/=(npoints*nclasses); rep.m_oobavgrelerror/=CApServ::Coalesce(oobavgrelcnt,1); } //+------------------------------------------------------------------+ //| This function performs binary compression of decision forest, | //| using either 8-bit mantissa (a bit more compact representation) | //| or 16-bit mantissa for splits and regression outputs. | //| Forest is compressed in - place. | //| Return value is a compression factor. | //+------------------------------------------------------------------+ double CDForest::BinaryCompression(CDecisionForest &df, bool usemantissa8) { //--- create variables double result=0; int size8=0; int size8i=0; int offssrc=0; int offsdst=0; int maxrawtreesize=0; CRowInt dummyi; CRowInt compressedsizes; //--- Quick exit if already compressed if(df.m_ForestFormat==m_DFCompressedV0) return(1.0); //--- Check that source format is supported if(!CAp::Assert(df.m_ForestFormat==m_DFUncompressedV0,__FUNCTION__": unexpected forest format")) return(result); //--- Compute sizes of uncompressed and compressed trees. for(int i=0; i=VAL then BRANCH0 else BRANCH1" //--- * stream value used for splitting //--- * stream children #0 and #1 StreamUint(buf,dstoffs,varidx+df.m_NVars); StreamFloat(buf,usemantissa8,dstoffs,splitval); StreamUint(buf,dstoffs,child1size); CompressRec(df,usemantissa8,treeroot,treeroot+jmponbranch,compressedsizes,buf,dstoffs); CompressRec(df,usemantissa8,treeroot,treepos+m_InnerNodeWidth,compressedsizes,buf,dstoffs); } } //--- Integrity check at the end if(!CAp::Assert(dstoffs-dstoffsold==compressedsizes[treepos-treeroot],__FUNCTION__": integrity check failed (compressed size at leaf)")) return; } //+------------------------------------------------------------------+ //| This function returns exact number of bytes required to store | //| compressed unsigned integer number(negative arguments result in | //| assertion being generated). | //+------------------------------------------------------------------+ int CDForest::ComputeCompressedUintSize(int v) { //--- check if(!CAp::Assert(v>=0)) return(0); //--- create variables int result=1; //--- while(v>=128) { v=v/128; result=result+1; } //--- return result return(result); } //+------------------------------------------------------------------+ //| This function stores compressed unsigned integer number (negative| //| arguments result in assertion being generated) to byte array at | //| location Offs and increments Offs by number of bytes being stored| //+------------------------------------------------------------------+ void CDForest::StreamUint(CRowInt &buf, int &offs, int v) { int v0=0; //--- check if(!CAp::Assert(v>=0)) return; //--- while(true) { //--- Save 7 least significant bits of V, use 8th bit as a flag which //--- tells us whether subsequent 7-bit packages will be sent. v0=v%128; if(v>=128) v0=v0+128; buf.Set(offs,(uint)v0); offs++; v/=128; if(v==0) break; } } //+------------------------------------------------------------------+ //| This function reads compressed unsigned integer number from byte | //| array starting at location Offs and increments Offs by number of | //| bytes being read. | //+------------------------------------------------------------------+ int CDForest::UnstreamUint(CRowInt &buf,int &offs) { //--- create variables int result=0; int v0=0; int p=1; //--- while(true) { //--- Rad 7 bits of V, use 8th bit as a flag which tells us whether //--- subsequent 7-bit packages will be received. v0=buf[offs]; offs=offs+1; result+=v0%128*p; if(v0<128) break; p*=128; } //--- return result return(result); } //+------------------------------------------------------------------+ //| This function stores compressed floating point number to byte | //| array at location Offs and increments Offs by number of bytes | //| being stored. Either 8-bit mantissa or 16-bit mantissa is used. | //| The exponent is always 7 bits of exponent + sign. Values which | //| do not fit into exponent range are truncated to fit. | //+------------------------------------------------------------------+ void CDForest::StreamFloat(CRowInt &buf, bool usemantissa8, int &offs, double v) { //--- create variables int signbit=0; int e=0; int m=0; double twopow30=0; double twopowm30=0; double twopow10=0; double twopowm10=0; //--- check if(!CAp::Assert(MathIsValidNumber(v),__FUNCTION__": V is not finite number")) return; //--- Special case: zero if(v==0.0) { if(usemantissa8) { buf.Set(offs,0); buf.Set(offs+1,0); offs+=2; } else { buf.Set(offs,0); buf.Set(offs+1,0); buf.Set(offs+2,0); offs+=3; } return; } //--- Handle sign signbit=0; if(v<0.0) { v=-v; signbit=128; } //--- Compute exponent twopow30=1073741824; twopow10=1024; twopowm30=1.0/twopow30; twopowm10=1.0/twopow10; e=0; while(v>=twopow30) { v*=twopowm30; e+=30; } while(v>=twopow10) { v*=twopowm10; e+=10; } while(v>=1.0) { v*=0.5; e++; } while(v=0.5 && v<1.0,__FUNCTION__": integrity check failed")) return; //--- Handle exponent underflow/overflow if(e<-63) { signbit=0; e=0; v=0; } if(e>63) { e=63; v=1.0; } //--- Save to stream if(usemantissa8) { m=(int)MathRound(v*256); if(m==256) { m=m/2; e=MathMin(e+1,63); } buf.Set(offs,e+64+signbit); buf.Set(offs+1,m); offs+=2; } else { m=(int)MathRound(v*65536); if(m==65536) { m=m/2; e=MathMin(e+1,63); } buf.Set(offs,e+64+signbit); buf.Set(offs+1,m%256); buf.Set(offs+2,m/256); offs+=3; } } //+------------------------------------------------------------------+ //| This function reads compressed floating point number from the | //| byte array starting from location Offs and increments Offs by | //| number of bytes being read. | //| Either 8-bit mantissa or 16-bit mantissa is used. The exponent | //| is always 7 bits of exponent + sign. Values which do not fit | //| into exponent range are truncated to fit. | //+------------------------------------------------------------------+ double CDForest::UnstreamFloat(CRowInt &buf, bool usemantissa8, int &offs) { //--- create variables double result=0; int e=0; double v=0; double inv256=1.0/256.0; //--- Read from stream if(usemantissa8) { e=buf[offs]; v=buf[offs+1]*inv256; offs+=2; } else { e=buf[offs]; v=(buf[offs+1]*inv256+buf[offs+2])*inv256; offs+=3 ; } //--- Decode if(e>128) { v=-v; e=e-128; } e=e-64; result=MathPow(2.0,e)*v; //--- return result return(result); } //+------------------------------------------------------------------+ //| Classification error | //+------------------------------------------------------------------+ int CDForest::DFClsError(CDecisionForest &df,CMatrixDouble &xy, const int npoints) { //--- create variables int result=0; int k=0; int tmpi=0; int i_=0; //--- creating arrays CRowDouble x; CRowDouble y; //--- check if(df.m_NClasses<=1) return(0); //--- initialization result=0; for(int i=0; iy[tmpi]) tmpi=j; } //--- check if(tmpi!=k) result++; } //--- return result return(result); } //+------------------------------------------------------------------+ //| Internal subroutine for processing one decision tree stored in | //| uncompressed format starting at SubtreeRoot (this index points | //| to the header of the tree, not its first node). First node being | //| processed is located at NodeOffs. | //+------------------------------------------------------------------+ void CDForest::DFProcessInternalUncompressed(CDecisionForest &df, int subtreeroot, int nodeoffs, CRowDouble &x, CRowDouble &y) { int idx=0; //--- check if(!CAp::Assert(df.m_ForestFormat==m_DFUncompressedV0,__FUNCTION__+": unexpected forest format")) return; //--- Navigate through the tree while(true) { if(df.m_Trees[nodeoffs]==(-1.0)) { if(df.m_NClasses==1) y.Add(0,df.m_Trees[nodeoffs+1]); else { idx=(int)MathRound(df.m_Trees[nodeoffs+1]); y.Add(idx,1); } break; } if(x[(int)MathRound(df.m_Trees[nodeoffs])]=splitval) offs+=jmplen; } else { //--- The split rule is "if VAR>=VAL then BRANCH0 else BRANCH1" varidx-=df.m_NVars; if(x[varidx]=K | //| NVars - number of variables, NVars>=1 | //| K - desired number of clusters, K>=1 | //| Restarts - number of restarts, Restarts>=1 | //| OUTPUT PARAMETERS: | //| Info - return code: | //| * -3, if task is degenerate (number of | //| distinct points is less than K) | //| * -1, if incorrect | //| NPoints/NFeatures/K/Restarts was passed| //| * 1, if subroutine finished successfully | //| C - array[0..NVars-1,0..K-1].matrix whose columns| //| store cluster's centers | //| XYC - array[NPoints], which contains cluster | //| indexes | //+------------------------------------------------------------------+ void CKMeans::KMeansGenerate(CMatrixDouble &xy,const int npoints, const int nvars,const int k, const int restarts,int &info, CMatrixDouble &c,int &xyc[]) { //--- create variables int i=0; int j=0; double e=0; double ebest=0; double v=0; int cclosest=0; bool waschanges; bool zerosizeclusters; int pass=0; int i_=0; double dclosest=0; //--- creating arrays int xycbest[]; double x[]; double tmp[]; double d2[]; double p[]; int csizes[]; bool cbusy[]; double work[]; //--- create matrix CMatrixDouble ct; CMatrixDouble ctbest; //--- initialization info=0; //--- Test parameters if(npoints=0 | //| NVars - number of independent variables, NVars>=1 | //| NClasses - number of classes, NClasses>=2 | //| OUTPUT PARAMETERS: | //| Info - return code: | //| * -4, if internal EVD subroutine hasn't | //| converged | //| * -2, if there is a point with class number | //| outside of [0..NClasses-1]. | //| * -1, if incorrect parameters was passed | //| (NPoints<0, NVars<1, NClasses<2) | //| * 1, if task has been solved | //| * 2, if there was a multicollinearity in | //| training set, but task has been solved.| //| W - linear combination coefficients, | //| array[0..NVars-1] | //+------------------------------------------------------------------+ void CLDA::FisherLDA(CMatrixDouble &xy,const int npoints, const int nvars,const int nclasses, int &info,double &w[]) { CRowDouble W; FisherLDA(xy,npoints,nvars,nclasses,info,W); W.ToArray(w); } //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ void CLDA::FisherLDA(CMatrixDouble &xy,const int npoints, const int nvars,const int nclasses, int &info,CRowDouble &w) { //--- create matrix CMatrixDouble w2; //--- initialization info=0; //--- function call FisherLDAN(xy,npoints,nvars,nclasses,info,w2); //--- check and copy if(info>0) w=w2.Col(0)+0; } //+------------------------------------------------------------------+ //| N-dimensional multiclass Fisher LDA | //| Subroutine finds coefficients of linear combinations which | //| optimally separates | //| training set on classes. It returns N-dimensional basis whose | //| vector are sorted | //| by quality of training set separation (in descending order). | //| INPUT PARAMETERS: | //| XY - training set, array[0..NPoints-1,0..NVars]. | //| First NVars columns store values of | //| independent variables, next column stores | //| number of class (from 0 to NClasses-1) which | //| dataset element belongs to. Fractional values| //| are rounded to nearest integer. | //| NPoints - training set size, NPoints>=0 | //| NVars - number of independent variables, NVars>=1 | //| NClasses - number of classes, NClasses>=2 | //| OUTPUT PARAMETERS: | //| Info - return code: | //| * -4, if internal EVD subroutine hasn't | //| converged | //| * -2, if there is a point with class number | //| outside of [0..NClasses-1]. | //| * -1, if incorrect parameters was passed | //| (NPoints<0, NVars<1, NClasses<2) | //| * 1, if task has been solved | //| * 2, if there was a multicollinearity in | //| training set, but task has been solved.| //| W - basis, array[0..NVars-1,0..NVars-1] | //| columns of matrix stores basis vectors, | //| sorted by quality of training set separation | //| (in descending order) | //+------------------------------------------------------------------+ void CLDA::FisherLDAN(CMatrixDouble &xy,const int npoints,const int nvars, const int nclasses,int &info,CMatrixDouble &w) { //--- create variables int i=0; int j=0; int k=0; int m=0; double v=0; int i_=0; //--- creating arrays int c[]; vector mu; int nc[]; vector tf; CRowDouble d; CRowDouble d2; CRowDouble work; //--- create matrix CMatrixDouble muc; CMatrixDouble sw; CMatrixDouble st; CMatrixDouble z; CMatrixDouble z2; CMatrixDouble tm; CMatrixDouble sbroot; CMatrixDouble a; CMatrixDouble xyc; CMatrixDouble xyproj; CMatrixDouble wproj; //--- initialization info=0; //--- Test data if((npoints<0 || nvars<1) || nclasses<2) { info=-1; return; } mu=xy.Col(nvars); mu.Resize(npoints); //--- check if((int)MathRound(mu.Min())<0 || (int)MathRound(mu.Max())>=nclasses) { info=-2; return; } //--- change value info=1; //--- Special case: NPoints<=1 //--- Degenerate task. if(npoints<=1) { info=2; //--- initialization w=matrix::Identity(nvars,nvars); //--- exit the function return; } //--- Prepare temporaries tf.Resize(nvars); work.Resize(MathMax(nvars,npoints)+1); //--- Convert class labels from reals to integers (just for convenience) ArrayResizeAL(c,npoints); for(i=0; i::Zeros(nclasses,nvars+1); ArrayResize(nc,nclasses); ArrayInitialize(nc,0); //--- calculation for(i=0; i::Identity(nvars,nvars); //--- exit the function return; } //--- Special case: degenerate ST matrix,multicollinearity found. //--- Since we know ST eigenvalues/vectors we can translate task to //--- non-degenerate form. //--- Let WG is orthogonal basis of the non zero variance subspace //--- of the ST and let WZ is orthogonal basis of the zero variance //--- subspace. //--- Projection on WG allows us to use LDA on reduced M-dimensional //--- subspace,N-M vectors of WZ allows us to update reduced LDA //--- factors to full N-dimensional subspace. m=0; for(k=0; kNVars+1 | //| NVars - number of independent variables | //| OUTPUT PARAMETERS: | //| Info - return code: | //| * -255, in case of unknown internal error | //| * -4, if internal SVD subroutine haven't | //| converged | //| * -1, if incorrect parameters was passed | //| (NPoints::Ones(npoints); //--- function call LRBuildS(xy,s,npoints,nvars,info,lm,ar); //--- check if(info<0) return; //--- calculation sigma2=CMath::Sqr(ar.m_RMSError)*npoints/(npoints-nvars-1); ar.m_c*=sigma2; } //+------------------------------------------------------------------+ //| Linear regression | //| Variant of LRBuild which uses vector of standatd deviations | //| (errors in function values). | //| INPUT PARAMETERS: | //| XY - training set, array [0..NPoints-1,0..NVars]: | //| * NVars columns - independent variables | //| * last column - dependent variable | //| S - standard deviations (errors in function | //| values) array[0..NPoints-1], S[i]>0. | //| NPoints - training set size, NPoints>NVars+1 | //| NVars - number of independent variables | //| OUTPUT PARAMETERS: | //| Info - return code: | //| * -255, in case of unknown internal error | //| * -4, if internal SVD subroutine haven't | //| converged | //| * -1, if incorrect parameters was passed | //| (NPoints means; vector sigmas=vector::Zeros(nvars); CRowDouble x; //--- create array CMatrixDouble xyi; //--- Copy data,add one more column (constant term) xyi=xy; xyi.Resize(npoints,nvars+2); xyi.Col(nvars+1,vector::Ones(npoints)); xyi.SwapCols(nvars,nvars+1); //--- Standartization means=xyi.Mean(0); means.Resize(nvars); for(int i=0; i::Zeros(npoints)); xyi.SwapCols(nvars,nvars+1); //--- Standartization: unusual scaling c.Resize(nvars); for(int j=0; jMathSqrt(variance)) { //--- variation is relatively small,it is better to //--- bring mean value to 1 c.Set(j,mean); } else { //--- variation is large,it is better to bring variance to 1 if(variance==0.0) c.Set(j,1.0); else c.Set(j,MathSqrt(variance)); } xyi.Col(j,x/c[j]+0); } //--- Internal processing LRInternal(xyi,s,npoints,nvars+1,info,lm,ar); //--- check if(info<0) return; //--- Un-standartization offs=(int)MathRound(lm.m_w[3]); for(int j=0; j::Ones(npoints); //--- initialization info=0; //--- function call LRBuildZS(xy,s,npoints,nvars,info,lm,ar); //--- check if(info<0) return; //--- calculation double sigma2=CMath::Sqr(ar.m_RMSError)*npoints/(npoints-nvars-1); ar.m_c*=sigma2; } //+------------------------------------------------------------------+ //| Unpacks coefficients of linear model. | //| INPUT PARAMETERS: | //| LM - linear model in ALGLIB format | //| OUTPUT PARAMETERS: | //| V - coefficients,array[0..NVars] | //| constant term (intercept) is stored in the | //| V[NVars]. | //| NVars - number of independent variables (one less | //| than number of coefficients) | //+------------------------------------------------------------------+ void CLinReg::LRUnpack(CLinearModel &lm,double &v[],int &nvars) { CRowDouble V; LRUnpack(lm,V,nvars); V.ToArray(v); } //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ void CLinReg::LRUnpack(CLinearModel &lm,CRowDouble &v,int &nvars) { //--- create variables int offs=0; int i1_=0; //--- initialization nvars=0; //--- check if(!CAp::Assert((int)MathRound(lm.m_w[1])==m_lrvnum,__FUNCTION__+": Incorrect LINREG version!")) { v.Resize(0); return; } //--- change values nvars=(int)MathRound(lm.m_w[2]); offs=(int)MathRound(lm.m_w[3]); //--- allocation v.Resize(nvars+1); //--- calculation i1_=offs; for(int i_=0; i_<=nvars; i_++) v.Set(i_,lm.m_w[i_+i1_]); } //+------------------------------------------------------------------+ //| "Packs" coefficients and creates linear model in ALGLIB format | //| (LRUnpack reversed). | //| INPUT PARAMETERS: | //| V - coefficients, array[0..NVars] | //| NVars - number of independent variables | //| OUTPUT PAREMETERS: | //| LM - linear model. | //+------------------------------------------------------------------+ void CLinReg::LRPack(double &v[],const int nvars,CLinearModel &lm) { CRowDouble V=v; LRPack(V,nvars,lm); } //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ void CLinReg::LRPack(CRowDouble &v,const int nvars,CLinearModel &lm) { //--- create variables int offs=0; int i1_=0; //--- allocation lm.m_w.Resize(5+nvars); //--- change values offs=4; lm.m_w.Set(0,4+nvars+1); lm.m_w.Set(1,m_lrvnum); lm.m_w.Set(2,nvars); lm.m_w.Set(3,offs); //--- calculation i1_=-offs; for(int i_=offs; i_<=offs+nvars; i_++) lm.m_w.Set(i_,v[i_+i1_]); } //+------------------------------------------------------------------+ //| Procesing | //| INPUT PARAMETERS: | //| LM - linear model | //| X - input vector, array[0..NVars-1]. | //| Result: | //| value of linear model regression estimate | //+------------------------------------------------------------------+ double CLinReg::LRProcess(CLinearModel &lm,double &x[]) { CRowDouble X=x; return(LRProcess(lm,X)); } //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ double CLinReg::LRProcess(CLinearModel &lm,CRowDouble &x) { //--- create variables double v=0; int offs=0; int nvars=0; int i1_=0; //--- check if(!CAp::Assert((int)MathRound(lm.m_w[1])==m_lrvnum,__FUNCTION__+": Incorrect LINREG version!")) return(EMPTY_VALUE); //--- change values nvars=(int)MathRound(lm.m_w[2]); offs=(int)MathRound(lm.m_w[3]); i1_=offs; v=0.0; //--- calculation for(int i_=0; i_ X=xy.Col(0); vector Y=xy.Col(1); vector S=s.ToVector(); S.Resize(n); X.Resize(n); Y.Resize(n); vector SS=MathPow(S*S,-1); ss=SS.Sum(); sx=(X*SS).Sum(); sy=(Y*SS).Sum(); sxx=(X*X*SS).Sum(); //--- Test for condition number t=MathSqrt(4*CMath::Sqr(sx)+CMath::Sqr(ss-sxx)); e1=0.5*(ss+sxx+t); e2=0.5*(ss+sxx-t); //--- check if(MathMin(e1,e2)<=1000*CMath::m_machineepsilon*MathMax(e1,e2)) { info=-3; return; } //--- Calculate A,B vector T=(X-sx/ss)/S; b=((T*Y)/S).Sum(); stt=(T*T).Sum(); b=b/stt; a=(sy-sx*b)/ss; //--- Calculate goodness-of-fit if(n>2) { chi2=(MathPow((Y-X*b-a)/S,2)).Sum(); //--- function call p=CIncGammaF::IncompleteGammaC((double)(n-2)/2.0,chi2/2.0); } else p=1; //--- Calculate other parameters vara=(1+CMath::Sqr(sx)/(ss*stt))/ss; varb=1/stt; covab=-(sx/(ss*stt)); corrab=covab/MathSqrt(vara*varb); } //+------------------------------------------------------------------+ //| Class method | //+------------------------------------------------------------------+ void CLinReg::LRLine(CMatrixDouble &xy,const int n,int &info, double &a,double &b) { //--- initialization info=0; a=0; b=0; //--- check if(n<2) { info=-1; return; } //--- create variables double vara=0; double varb=0; double covab=0; double corrab=0; double p=0; //--- create array CRowDouble s=vector::Ones(n); //--- function call LRLines(xy,s,n,info,a,b,vara,varb,covab,corrab,p); } //+------------------------------------------------------------------+ //| Internal linear regression subroutine | //+------------------------------------------------------------------+ void CLinReg::LRInternal(CMatrixDouble &xy,CRowDouble &s,const int npoints, const int nvars,int &info,CLinearModel &lm, CLRReport &ar) { //--- Check for errors in data if(npoints::Zeros(nvars,nvars); return; } if(sv[nvars-1]<=(epstol*CMath::m_machineepsilon*sv[0])) { //--- Degenerate case, non-zero design matrix. //--- We can leave it and solve task in SVD least squares fashion. //--- Solution and covariance matrix will be obtained correctly, //--- but CV error estimates - will not. It is better to reduce //--- it to non-degenerate task and to obtain correct CV estimates. for(k=nvars; k>=1; k--) { if(sv[k-1]>(epstol*CMath::m_machineepsilon*sv[0])) { //--- Reduce xym.Resize(npoints,k+1); for(i=0; i(epstol*CMath::m_machineepsilon*sv[0])) svi.Set(i,1/sv[i]); else svi.Set(i,0); } t=vector::Zeros(nvars); for(i=0; i(1-epstol*CMath::m_machineepsilon)) { ar.m_cvdefects.Set(ar.m_ncvdefects,i); ar.m_ncvdefects++; continue; } r=s[i]*(r/s[i]-b[i]*p)/(1-p); ar.m_cvrmserror+=CMath::Sqr(r-xy.Get(i,nvars)); ar.m_cvavgerror+=MathAbs(r-xy.Get(i,nvars)); if(xy.Get(i,nvars)!=0.0) { ar.m_cvavgrelerror+=MathAbs((r-xy.Get(i,nvars))/xy.Get(i,nvars)); nacv++; } ncv++; } if(ncv==0) { //--- Something strange: ALL ui are degenerate. //--- Unexpected... info=-255; return; } ar.m_RMSError=MathSqrt(ar.m_RMSError/npoints); ar.m_AvgError=ar.m_AvgError/npoints; if(na!=0) ar.m_AvgRelError=ar.m_AvgRelError/na; ar.m_cvrmserror=MathSqrt(ar.m_cvrmserror/ncv); ar.m_cvavgerror=ar.m_cvavgerror/ncv; if(nacv!=0) ar.m_cvavgrelerror=ar.m_cvavgrelerror/nacv; } //+------------------------------------------------------------------+ //| Model's errors: | //| * RelCLSError - fraction of misclassified cases. | //| * AvgCE - acerage cross-entropy | //| * RMSError - root-mean-square error | //| * AvgError - average error | //| * AvgRelError - average relative error | //| NOTE 1: RelCLSError/AvgCE are zero on regression problems. | //| NOTE 2: on classification problems RMSError/AvgError/AvgRelError | //| contain errors in prediction of posterior probabilities | //+------------------------------------------------------------------+ struct CModelErrors { double m_RelCLSError; double m_AvgCE; double m_RMSError; double m_AvgError; double m_AvgRelError; //--- CModelErrors(void) { ZeroMemory(this); } ~CModelErrors(void) {} //--- void Copy(const CModelErrors &obj); //--- overloading void operator=(const CModelErrors &obj) { Copy(obj); } }; //+------------------------------------------------------------------+ //| Copy | //+------------------------------------------------------------------+ void CModelErrors::Copy(const CModelErrors &obj) { m_RelCLSError=obj.m_RelCLSError; m_AvgCE=obj.m_AvgCE; m_RMSError=obj.m_RMSError; m_AvgError=obj.m_AvgError; m_AvgRelError=obj.m_AvgRelError; } //+------------------------------------------------------------------+ //|This structure is used to store MLP error and gradient. | //+------------------------------------------------------------------+ struct CSMLPGrad { double m_F; CRowDouble m_G; //--- CSMLPGrad(void) { m_F=0; } ~CSMLPGrad(void) {} //--- void Copy(const CSMLPGrad &obj) { m_F=obj.m_F; m_G=obj.m_G; } //--- overloading void operator=(const CSMLPGrad &obj) { Copy(obj); } }; //+------------------------------------------------------------------+ //| Auxiliary class for CMLPBase | //+------------------------------------------------------------------+ class CMultilayerPerceptron { public: //--- variables int m_hlnetworktype; int m_hlnormtype; CModelErrors m_err; CSMLPGrad m_grad; //--- arrays CRowInt m_dummyidx; CRowInt m_hllayersizes; CRowInt m_hlconnections; CRowInt m_hlneurons; CRowInt m_integerbuf; CRowInt m_structinfo; CRowDouble m_weights; CRowDouble m_columnmeans; CRowDouble m_columnsigmas; CRowDouble m_neurons; CRowDouble m_dfdnet; CRowDouble m_derror; CRowDouble m_x; CRowDouble m_y; CRowDouble m_xyrow; CRowDouble m_nwbuf; CRowDouble m_rndbuf; //--- matrix CMatrixDouble m_xy; CMatrixDouble m_dummydxy; CSparseMatrix m_dummysxy; //--- constructor, destructor CMultilayerPerceptron(void) { m_hlnetworktype=0; m_hlnormtype=0; } ~CMultilayerPerceptron(void) {} //--- copy void Copy(const CMultilayerPerceptron &obj); //--- overloading void operator=(const CMultilayerPerceptron &obj) { Copy(obj); } }; //+------------------------------------------------------------------+ //| Copy | //+------------------------------------------------------------------+ void CMultilayerPerceptron::Copy(const CMultilayerPerceptron &obj) { //--- copy variables m_hlnetworktype=obj.m_hlnetworktype; m_hlnormtype=obj.m_hlnormtype; m_err=obj.m_err; m_grad=obj.m_grad; //--- copy arrays m_dummyidx=obj.m_dummyidx; m_hllayersizes=obj.m_hllayersizes; m_hlconnections=obj.m_hlconnections; m_hlneurons=obj.m_hlneurons; m_integerbuf=obj.m_integerbuf; m_structinfo=obj.m_structinfo; m_weights=obj.m_weights; m_columnmeans=obj.m_columnmeans; m_columnsigmas=obj.m_columnsigmas; m_neurons=obj.m_neurons; m_dfdnet=obj.m_dfdnet; m_derror=obj.m_derror; m_x=obj.m_x; m_y=obj.m_y; m_xyrow=obj.m_xyrow; m_nwbuf=obj.m_nwbuf; m_rndbuf=obj.m_rndbuf; //--- copy matrix m_xy=obj.m_xy; m_dummydxy=obj.m_dummydxy; m_dummysxy=obj.m_dummysxy; } //+------------------------------------------------------------------+ //| This class is a shell for class CMultilayerPerceptron | //+------------------------------------------------------------------+ class CMultilayerPerceptronShell { private: CMultilayerPerceptron m_innerobj; public: //--- constructors, destructor CMultilayerPerceptronShell(void) {} CMultilayerPerceptronShell(CMultilayerPerceptron &obj) { m_innerobj.Copy(obj); } ~CMultilayerPerceptronShell(void) {} //--- method CMultilayerPerceptron *GetInnerObj(void) { return(GetPointer(m_innerobj)); } }; //+------------------------------------------------------------------+ //| Multilayer perceptron class | //+------------------------------------------------------------------+ class CMLPBase { public: //--- variables static const int m_mlpvnum; static const int m_mlpfirstversion; static const int m_nfieldwidth; static const int m_hlconm_nfieldwidth; static const int m_hlm_nfieldwidth; static const int m_gradbasecasecost; static const int m_microbatchsize; //--- public methods static int MLPGradSplitCost(void); static int MLPGradSplitSize(void); static void MLPCreate0(const int nin,const int nout,CMultilayerPerceptron &network); static void MLPCreate1(const int nin,const int nhid,const int nout,CMultilayerPerceptron &network); static void MLPCreate2(const int nin,const int nhid1,const int nhid2,const int nout,CMultilayerPerceptron &network); static void MLPCreateB0(const int nin,const int nout,const double b,double d,CMultilayerPerceptron &network); static void MLPCreateB1(const int nin,const int nhid,const int nout,const double b,double d,CMultilayerPerceptron &network); static void MLPCreateB2(const int nin,const int nhid1,const int nhid2,const int nout,const double b,double d,CMultilayerPerceptron &network); static void MLPCreateR0(const int nin,const int nout,const double a,const double b,CMultilayerPerceptron &network); static void MLPCreateR1(const int nin,const int nhid,const int nout,const double a,const double b,CMultilayerPerceptron &network); static void MLPCreateR2(const int nin,const int nhid1,const int nhid2,const int nout,const double a,const double b,CMultilayerPerceptron &network); static void MLPCreateC0(const int nin,const int nout,CMultilayerPerceptron &network); static void MLPCreateC1(const int nin,const int nhid,const int nout,CMultilayerPerceptron &network); static void MLPCreateC2(const int nin,const int nhid1,const int nhid2,const int nout,CMultilayerPerceptron &network); static void MLPCopy(const CMultilayerPerceptron &network1,CMultilayerPerceptron &network2); static bool MLPSameArchitecture(CMultilayerPerceptron &network1,CMultilayerPerceptron &network2); static void MLPCopyTunableParameters(CMultilayerPerceptron &network1,CMultilayerPerceptron &network2); static void MLPExportTunableParameters(CMultilayerPerceptron &network,CRowDouble &p,int &pcount); static void MLPImportTunableParameters(CMultilayerPerceptron &network,CRowDouble &p); static void MLPSerializeOld(CMultilayerPerceptron &network,double &ra[],int &rlen); static void MLPSerializeOld(CMultilayerPerceptron &network,CRowDouble &ra,int &rlen); static void MLPUnserializeOld(double &ra[],CMultilayerPerceptron &network); static void MLPUnserializeOld(CRowDouble &ra,CMultilayerPerceptron &network); static void MLPRandomize(CMultilayerPerceptron &network); static void MLPRandomizeFull(CMultilayerPerceptron &network); static void MLPInitPreprocessor(CMultilayerPerceptron &network,CMatrixDouble &xy,const int ssize); static void MLPInitPreprocessorSparse(CMultilayerPerceptron &network,CSparseMatrix &xy,int ssize); static void MLPInitPreprocessorSubset(CMultilayerPerceptron &network,CMatrixDouble &xy,int setsize,CRowInt &idx,int subsetsize); static void MLPInitPreprocessorSparseSubset(CMultilayerPerceptron &network,CSparseMatrix &xy,int setsize,CRowInt &idx,int subsetsize); static void MLPProperties(CMultilayerPerceptron &network,int &nin,int &nout,int &wcount); static int MLPNTotal(CMultilayerPerceptron &network); static int MLPGetInputsCount(CMultilayerPerceptron &network); static int MLPGetOutputsCount(CMultilayerPerceptron &network); static int MLPGetWeightsCount(CMultilayerPerceptron &network); static bool MLPIsSoftMax(CMultilayerPerceptron &network); static int MLPGetLayersCount(CMultilayerPerceptron &network); static int MLPGetLayerSize(CMultilayerPerceptron &network,const int k); static void MLPGetInputScaling(CMultilayerPerceptron &network,const int i,double &mean,double &sigma); static void MLPGetOutputScaling(CMultilayerPerceptron &network,const int i,double &mean,double &sigma); static void MLPGetNeuronInfo(CMultilayerPerceptron &network,const int k,const int i,int &fkind,double &threshold); static double MLPGetWeight(CMultilayerPerceptron &network,const int k0,const int i0,const int k1,const int i1); static void MLPSetInputScaling(CMultilayerPerceptron &network,const int i,const double mean,double sigma); static void MLPSetOutputScaling(CMultilayerPerceptron &network,const int i,const double mean,double sigma); static void MLPSetNeuronInfo(CMultilayerPerceptron &network,const int k,const int i,const int fkind,const double threshold); static void MLPSetWeight(CMultilayerPerceptron &network,const int k0,const int i0,const int k1,const int i1,const double w); static void MLPActivationFunction(double net,const int k,double &f,double &df,double &d2f); static void MLPProcess(CMultilayerPerceptron &network,double &x[],double &y[]); static void MLPProcess(CMultilayerPerceptron &network,CRowDouble &x,CRowDouble &y); static void MLPProcessI(CMultilayerPerceptron &network,double &x[],double &y[]); static void MLPProcessI(CMultilayerPerceptron &network,CRowDouble &x,CRowDouble &y); static double MLPError(CMultilayerPerceptron &network,CMatrixDouble &xy,const int npoints); static double MLPErrorSparse(CMultilayerPerceptron &network,CSparseMatrix &xy,int npoints); static double MLPErrorN(CMultilayerPerceptron &network,CMatrixDouble &xy,const int ssize); static int MLPClsError(CMultilayerPerceptron &network,CMatrixDouble &xy,const int ssize); static double MLPRelClsError(CMultilayerPerceptron &network,CMatrixDouble &xy,const int npoints); static double MLPRelClsErrorSparse(CMultilayerPerceptron &network,CSparseMatrix &xy,int npoints); static double MLPAvgCE(CMultilayerPerceptron &network,CMatrixDouble &xy,const int npoints); static double MLPAvgCESparse(CMultilayerPerceptron &network,CSparseMatrix &xy,int npoints); static double MLPRMSError(CMultilayerPerceptron &network,CMatrixDouble &xy,const int npoints); static double MLPRMSErrorSparse(CMultilayerPerceptron &network,CSparseMatrix &xy,int npoints); static double MLPAvgError(CMultilayerPerceptron &network,CMatrixDouble &xy,const int npoints); static double MLPAvgErrorSparse(CMultilayerPerceptron &network,CSparseMatrix &xy,int npoints); static double MLPAvgRelError(CMultilayerPerceptron &network,CMatrixDouble &xy,const int npoints); static double MLPAvgRelErrorSparse(CMultilayerPerceptron &network,CSparseMatrix &xy,int npoints); static void MLPGrad(CMultilayerPerceptron &network,double &x[],double &desiredy[],double &e,double &grad[]); static void MLPGrad(CMultilayerPerceptron &network,CRowDouble &x,CRowDouble &desiredy,double &e,CRowDouble &grad); static void MLPGradN(CMultilayerPerceptron &network,double &x[],double &desiredy[],double &e,double &grad[]); static void MLPGradN(CMultilayerPerceptron &network,CRowDouble &x,CRowDouble &desiredy,double &e,CRowDouble &grad); static void MLPGradBatch(CMultilayerPerceptron &network,CMatrixDouble &xy,const int ssize,double &e,double &grad[]); static void MLPGradBatch(CMultilayerPerceptron &network,CMatrixDouble &xy,const int ssize,double &e,CRowDouble &grad); static void MLPGradBatchSparse(CMultilayerPerceptron &network,CSparseMatrix &xy,int ssize,double &e,CRowDouble &grad); static void MLPGradBatchSubset(CMultilayerPerceptron &network,CMatrixDouble &xy,int setsize,CRowInt &idx,int subsetsize,double &e,CRowDouble &grad); static void MLPGradBatchSparseSubset(CMultilayerPerceptron &network,CSparseMatrix &xy,int setsize,CRowInt &idx,int subsetsize,double &e,CRowDouble &grad); static void MLPGradBatchX(CMultilayerPerceptron &network,CMatrixDouble &densexy,CSparseMatrix &sparsexy,int datasetsize,int datasettype,CRowInt &idx,int subset0,int subset1,int subsettype,double &e,CRowDouble &grad); static void MLPGradNBatch(CMultilayerPerceptron &network,CMatrixDouble &xy,const int ssize,double &e,double &grad[]); static void MLPGradNBatch(CMultilayerPerceptron &network,CMatrixDouble &xy,const int ssize,double &e,CRowDouble &grad); static void MLPHessianNBatch(CMultilayerPerceptron &network,CMatrixDouble &xy,const int ssize,double &e,double &grad[],CMatrixDouble &h); static void MLPHessianNBatch(CMultilayerPerceptron &network,CMatrixDouble &xy,const int ssize,double &e,CRowDouble &grad,CMatrixDouble &h); static void MLPHessianBatch(CMultilayerPerceptron &network,CMatrixDouble &xy,const int ssize,double &e,double &grad[],CMatrixDouble &h); static void MLPHessianBatch(CMultilayerPerceptron &network,CMatrixDouble &xy,const int ssize,double &e,CRowDouble &grad,CMatrixDouble &h); static void MLPInternalProcessVector(int &structinfo[],double &weights[],double &columnmeans[],double &columnsigmas[],double &neurons[],double &dfdnet[],double &x[],double &y[]); static void MLPInternalProcessVector(CRowInt &structinfo,CRowDouble &weights,CRowDouble &columnmeans,CRowDouble &columnsigmas,CRowDouble &neurons,CRowDouble &dfdnet,CRowDouble &x,CRowDouble &y); static void MLPAlloc(CSerializer &s,CMultilayerPerceptron &network); static void MLPSerialize(CSerializer &s,CMultilayerPerceptron &network); static void MLPUnserialize(CSerializer &s,CMultilayerPerceptron &network); static void MLPAllErrorsSubset(CMultilayerPerceptron &network,CMatrixDouble &xy,int setsize,CRowInt &subset,int subsetsize,CModelErrors &rep); static void MLPAllErrorsSparseSubset(CMultilayerPerceptron &network,CSparseMatrix &xy,int setsize,CRowInt &subset,int subsetsize,CModelErrors &rep); static double MLPErrorSubset(CMultilayerPerceptron &network,CMatrixDouble &xy,int setsize,CRowInt &subset,int subsetsize); static double MLPErrorSparseSubset(CMultilayerPerceptron &network,CSparseMatrix &xy,int setsize,CRowInt &subset,int subsetsize); static void MLPAllErrorsX(CMultilayerPerceptron &network,CMatrixDouble &densexy,CSparseMatrix &sparsexy,int datasetsize,int datasettype,CRowInt &idx,int subset0,int subset1,int subsettype,CModelErrors &rep); private: static void AddInputLayer(const int ncount,CRowInt &lsizes,CRowInt <ypes,CRowInt &lconnfirst,CRowInt &lconnlast,int &lastproc); static void AddBiasedSummatorLayer(const int ncount,CRowInt &lsizes,CRowInt <ypes,CRowInt &lconnfirst,CRowInt &lconnlast,int &lastproc); static void AddActivationLayer(const int functype,CRowInt &lsizes,CRowInt <ypes,CRowInt &lconnfirst,CRowInt &lconnlast,int &lastproc); static void AddZeroLayer(CRowInt &lsizes,CRowInt <ypes,CRowInt &lconnfirst,CRowInt &lconnlast,int &lastproc); static void HLAddInputLayer(CMultilayerPerceptron &network,int &connidx,int &neuroidx,int &structinfoidx,int nin); static void HLAddOutputLayer(CMultilayerPerceptron &network,int &connidx,int &neuroidx,int &structinfoidx,int &weightsidx,const int k,const int nprev,const int nout,const bool iscls,const bool islinearout); static void HLAddHiddenLayer(CMultilayerPerceptron &network,int &connidx,int &neuroidx,int &structinfoidx,int &weightsidx,const int k,const int nprev,const int ncur); static void FillHighLevelInformation(CMultilayerPerceptron &network,const int nin,const int nhid1,const int nhid2,const int nout,const bool iscls,const bool islinearout); static void MLPCreate(const int nin,const int nout,CRowInt &lsizes,CRowInt <ypes,CRowInt &lconnfirst,CRowInt &lconnlast,const int layerscount,const bool isclsnet,CMultilayerPerceptron &network); static void MLPHessianBatchInternal(CMultilayerPerceptron &network,CMatrixDouble &xy,const int ssize,const bool naturalerr,double &e,CRowDouble &grad,CMatrixDouble &h); static void MLPInternalCalculateGradient(CMultilayerPerceptron &network,CRowDouble &neurons,CRowDouble &weights,CRowDouble &derror,CRowDouble &grad,const bool naturalerrorfunc); static void MLPChunkedGradient(CMultilayerPerceptron &network,CMatrixDouble &xy,const int cstart,const int csize,CRowDouble &batch4buf,CRowDouble &hpcbuf,double &e,const bool naturalerrorfunc); static void MLPChunkedProcess(CMultilayerPerceptron &network,CMatrixDouble &xy,int cstart,int csize,CRowDouble &batch4buf,CRowDouble &hpcbuf); static double SafeCrossEntropy(const double t,const double z); static void RandomizeBackwardPass(CMultilayerPerceptron &network,int neuronidx,double v); }; //+------------------------------------------------------------------+ //| Initialize constants | //+------------------------------------------------------------------+ const int CMLPBase::m_mlpvnum=7; const int CMLPBase::m_mlpfirstversion=0; const int CMLPBase::m_nfieldwidth=4; const int CMLPBase::m_hlconm_nfieldwidth=5; const int CMLPBase::m_hlm_nfieldwidth=4; const int CMLPBase::m_gradbasecasecost=50000; const int CMLPBase::m_microbatchsize=64; //+------------------------------------------------------------------+ //| This function returns number of weights updates which is required| //| for gradient calculation problem to be splitted. | //+------------------------------------------------------------------+ int CMLPBase::MLPGradSplitCost(void) { return(m_gradbasecasecost); } //+------------------------------------------------------------------+ //| This function returns number of elements in subset of dataset | //| which is required for gradient calculation problem to be splitted| //+------------------------------------------------------------------+ int CMLPBase::MLPGradSplitSize(void) { return(m_microbatchsize); } //+------------------------------------------------------------------+ //| Creates neural network with NIn inputs, NOut outputs, | //| without hidden layers, with linear output layer. Network weights | //| are filled with small random values. | //+------------------------------------------------------------------+ void CMLPBase::MLPCreate0(const int nin,const int nout, CMultilayerPerceptron &network) { //--- create variables int layerscount=0; int lastproc=0; //--- creating arrays CRowInt lsizes; CRowInt ltypes; CRowInt lconnfirst; CRowInt lconnlast; //--- initialization layerscount=4; //--- Allocate arrays lsizes.Resize(layerscount); ltypes.Resize(layerscount); lconnfirst.Resize(layerscount); lconnlast.Resize(layerscount); //--- Layers AddInputLayer(nin,lsizes,ltypes,lconnfirst,lconnlast,lastproc); //--- function call AddBiasedSummatorLayer(nout,lsizes,ltypes,lconnfirst,lconnlast,lastproc); //--- function call AddActivationLayer(-5,lsizes,ltypes,lconnfirst,lconnlast,lastproc); //--- Create MLPCreate(nin,nout,lsizes,ltypes,lconnfirst,lconnlast,layerscount,false,network); //--- function call FillHighLevelInformation(network,nin,0,0,nout,false,true); } //+------------------------------------------------------------------+ //| Same as MLPCreate0, but with one hidden layer (NHid neurons) with| //| non-linear activation function. Output layer is linear. | //+------------------------------------------------------------------+ void CMLPBase::MLPCreate1(const int nin,const int nhid,const int nout, CMultilayerPerceptron &network) { //--- create variables int layerscount=0; int lastproc=0; //--- creating arrays CRowInt lsizes; CRowInt ltypes; CRowInt lconnfirst; CRowInt lconnlast; //--- create variables layerscount=7; //--- Allocate arrays lsizes.Resize(layerscount); ltypes.Resize(layerscount); lconnfirst.Resize(layerscount); lconnlast.Resize(layerscount); //--- Layers AddInputLayer(nin,lsizes,ltypes,lconnfirst,lconnlast,lastproc); //--- function call AddBiasedSummatorLayer(nhid,lsizes,ltypes,lconnfirst,lconnlast,lastproc); //--- function call AddActivationLayer(1,lsizes,ltypes,lconnfirst,lconnlast,lastproc); //--- function call AddBiasedSummatorLayer(nout,lsizes,ltypes,lconnfirst,lconnlast,lastproc); //--- function call AddActivationLayer(-5,lsizes,ltypes,lconnfirst,lconnlast,lastproc); //--- Create MLPCreate(nin,nout,lsizes,ltypes,lconnfirst,lconnlast,layerscount,false,network); //--- function call FillHighLevelInformation(network,nin,nhid,0,nout,false,true); } //+------------------------------------------------------------------+ //| Same as MLPCreate0,but with two hidden layers (NHid1 and NHid2 | //| neurons) with non-linear activation function. Output layer is | //| linear. | //| $ALL | //+------------------------------------------------------------------+ void CMLPBase::MLPCreate2(const int nin,const int nhid1,const int nhid2, const int nout,CMultilayerPerceptron &network) { //--- create variables int layerscount=0; int lastproc=0; //--- creating arrays CRowInt lsizes; CRowInt ltypes; CRowInt lconnfirst; CRowInt lconnlast; //--- initialization layerscount=10; //--- Allocate arrays lsizes.Resize(layerscount); ltypes.Resize(layerscount); lconnfirst.Resize(layerscount); lconnlast.Resize(layerscount); //--- Layers AddInputLayer(nin,lsizes,ltypes,lconnfirst,lconnlast,lastproc); //--- function call AddBiasedSummatorLayer(nhid1,lsizes,ltypes,lconnfirst,lconnlast,lastproc); //--- function call AddActivationLayer(1,lsizes,ltypes,lconnfirst,lconnlast,lastproc); //--- function call AddBiasedSummatorLayer(nhid2,lsizes,ltypes,lconnfirst,lconnlast,lastproc); //--- function call AddActivationLayer(1,lsizes,ltypes,lconnfirst,lconnlast,lastproc); //--- function call AddBiasedSummatorLayer(nout,lsizes,ltypes,lconnfirst,lconnlast,lastproc); //--- function call AddActivationLayer(-5,lsizes,ltypes,lconnfirst,lconnlast,lastproc); //--- Create MLPCreate(nin,nout,lsizes,ltypes,lconnfirst,lconnlast,layerscount,false,network); //--- function call FillHighLevelInformation(network,nin,nhid1,nhid2,nout,false,true); } //+------------------------------------------------------------------+ //| Creates neural network with NIn inputs, NOut outputs, without | //| hidden layers with non-linear output layer. Network weights are | //| filled with small random values. | //| Activation function of the output layer takes values: | //| (B, +INF), if D>=0 | //| or | //| (-INF, B), if D<0. | //+------------------------------------------------------------------+ void CMLPBase::MLPCreateB0(const int nin,const int nout,const double b, double d,CMultilayerPerceptron &network) { //--- create variables int layerscount=0; int lastproc=0; //--- creating arrays CRowInt lsizes; CRowInt ltypes; CRowInt lconnfirst; CRowInt lconnlast; //--- initialization layerscount=4; //--- check if(d>=0.0) d=1; else d=-1; //--- Allocate arrays lsizes.Resize(layerscount); ltypes.Resize(layerscount); lconnfirst.Resize(layerscount); lconnlast.Resize(layerscount); //--- Layers AddInputLayer(nin,lsizes,ltypes,lconnfirst,lconnlast,lastproc); //--- function call AddBiasedSummatorLayer(nout,lsizes,ltypes,lconnfirst,lconnlast,lastproc); //--- function call AddActivationLayer(3,lsizes,ltypes,lconnfirst,lconnlast,lastproc); //--- Create MLPCreate(nin,nout,lsizes,ltypes,lconnfirst,lconnlast,layerscount,false,network); //--- function call FillHighLevelInformation(network,nin,0,0,nout,false,false); //--- Turn on ouputs shift/scaling. for(int i=nin; i