// ntrees.cc File to demonstrate to calculate number of trees of each depth // W.B.Langdon cwi.nl // based on ntrees.cc r1.2 #define main_version "12 March 1999 \$Revision: 1.23 \$" //WBL 29 Dec 21 Convert to gcc version 4.8.5 and gcc version 9.3.1 //WBL 20 Aug 02 Add code from ntrees.cc r1.3 to make self contained //WBL 5 Apr 99 Add command line choice and make cache a proper class // Add calculation of mean shape //WBL 16 Mar 99 re-definition of depth so root depth=1 (not inside count) //WBL 14 Mar 99 Add gradient //WBL 12 Mar 99 Add calculation of percentiles //WBL 14 Dec 98 //compile g++ -o ntreeshape -O3 ntreeshape.cc //Running on SGI ./ntreeshape 1,141,2 1,100 // Includes main function /*Software bit rot help given by https://stackoverflow.com/questions/13049978/fatal-error-iostream-h-no-such-file-or-directory https://stackoverflow.com/questions/15185801/cout-was-not-declared-in-this-scope https://stackoverflow.com/questions/3223732/how-to-instruct-gcc-to-stop-after-5-errors https://stackoverflow.com/questions/27705689/dealing-with-ambiguous-declarations-in-c */ #include using namespace std; #include #include #include #include #include #include #include #define BOOL int #define TRUE 1 #define FALSE 0 //dont change! #define max_arity 4 double* log_fact; int log_fact_size = 0; void new_log_fact(const int x) { log_fact = new double[x+1]; log_fact_size = x+1; double f = 0; log_fact[0] = 0; for(int i=1;i=0 && s<=size && d>=0 && d<=depth) { //only binary trees address = addressS(s)*stride+d; return TRUE; } return FALSE; }; public: Cache(const int s, const int d): size(s), depth(d), stride(d+1) { const int datasize = (addressS(size)+1)*stride; data = new datatype[datasize]; memset(data,invalid,sizeof(datatype)*datasize); }; ~Cache() { delete[] data; }; inline BOOL cached(const int s, const int d, double& out) const { int a; if(address(s,d,a) && data[a].flag != INVALID ) { out = data[a].value; return TRUE; } else return FALSE; } inline void set(const int s, const int d, const double in) { int a; if(address(s,d,a)) data[a].value = in; } }; Cache* cache; //number of binary trees of size and max depth (root node is depth=0,size=1) double count(const int size, const int depth /*, const char* ds*/) { double ans; if(size<=0||depth<0)//||size<=2*depth) ans = 0; else if(depth==0) ans = (size==1)? 1 : 0; else { double d; if(cache->cached(size,depth,d)) return d; double sum1 = 0; double sum2 = 0; {double cmd1=1;//dummy initial value for(int m=2*depth-1;(m0);m+=2) { cmd1 = count(m,depth-1 /*,"1a"*/); for(int s=0;s=(2*depth-1))&&c1>0;m+=2){ c1 = count(m, depth-1 /*,"2a"*/); c2 = count(size-m-1,depth-1 /*,"2b"*/); sum2 += c1*c2; }} ans = 2*sum1+sum2; } //cout<<"count("<set(size,depth,ans); return ans; }//end count //as count but do quick sanity check first (root node is depth=1,size=1) double count2(const int i, const int d /*, const char* ds*/) { return ((i%2==1)&&(i0) cout<threshold) set(d,v); } void peakset(const int d, const double v) { if(v>value) set(d,v); } void cset(const int d, const double v) { if(loc==d) value = v; } void print() { pr2(loc,value," "); cout<<"\t"; } };//endclass pair class mm { pair min, pl, peak, pu, max; double sum; double estimatedmean; //for improved numerical stability double msum, msum2; public: mm(const int size) { int binary[5] = {0,0,size/2,0,0}; const int arityc[5] = {1,0,1,0,0}; const double total = exp(tree_count(size,binary,arityc)); pl.thresh(total*0.05); pu.thresh(total*0.95); min.loc = 1+int(log(size)/log(2)); //(root node is depth=1,size=1) pu.loc = min.loc; peak.loc = min.loc; pl.loc = min.loc; max.loc = (size+1)/2; max.value= exp((max.loc-1)*log(2)); sum = 0; const double pi = 3.1415927; estimatedmean = 2*sqrt(pi*(size-1)/2);//sedgewick p256 T5.8 ignore O(n**0.25) msum = 0; msum2 = 0; } void update(const int d, const double t) { const double dev = d-estimatedmean; msum += dev*t; msum2 += dev*dev*t; sum += t; min.cset(d,t); pl.update(d,sum,t); peak.peakset(d,t); max.cset(d,t); pu.update(d,sum,t); } void print(const int i, const char* s) { cout<0)? sqrt(var) : 0; cout<<"\t"<0&&lastt==0) pr(i-incr,d,lastt,"df"); if(t>0||lastt>0) pr(i , d, t,"df"); lastt=t; } cout<0) cout<max_z) max_z = dfd; //acurate enough for log,log scaling if(dfs>max_z) max_z = dfs; } } double llmaxz = log10(log10(max_z)); cout<0) { const double dfd = count2(i, d+1) - f0; const double dfs = count2(i+2,d) - f0; const double theta = atan2(dfs,dfd); const double sint = sin(theta); const double cost = cos(theta); //const double z1 = fabs(dfd*sint+dfs*cost); //const double z = (z1