Allink  v0.1
MatematicaSign.cpp
1 //#include "../include/Matematica.h"
2 #include <Matematica.h>
3 
4 #ifdef __GSL__
5 #include <gsl/gsl_randist.h>
6 double Matematica::Casuale(){
7  return gsl_rng_uniform(rng);
8 }
9 bool Matematica::InizializzaGaussiano(double Scarto,int N){
10  return 0;
11 }
12 double Matematica::Gaussiano(double Media,double Scarto){//Gia inizializzato
13  return Media + gsl_ran_gaussian(rng,Scarto);
14 }
15 #else
16 bool Matematica::InizializzaGaussiano(double Scarto,int N){
17  zigset(86947731);
18  return 0;
19 }
23 double Matematica::Gaussiano(double Media,double Scarto){
24  double Rand = genrand_real1() ;
25  double Length = -log( Rand ) ;
26  double Theta = 2. * M_PI * genrand_real2() ;
27  double kdX = sqrt( 2. * Length ) * cos( Theta ) ;
28  return Scarto * kdX + Media ;
29  return RandomGaussian(Media,Scarto);
30 }
32  return genrand_real2();
33  return RandomUniform();
34  int IA=16807,IM=2147483647,IQ=127773,IR=2836,NDIV=1+(IM-1)/NTAB,kcas,jcas;
35  double AM,EPS,RNMX;
36  EPS=1.7*pow(10.,-7.);
37  RNMX=1.-EPS;
38  AM=1./IM;
39  static int iv[NTAB];
40  static int iy=10;
41  static int seme = 2351653;
42  // printf("%d,%d\n",seme,iy);
43  if(seme<=0 || iy==0){
44  seme = -seme>1 ? -seme : 1;
45  for(jcas=NTAB+8;jcas>=1;jcas--){
46  kcas=seme/IQ;
47  seme=IA*(seme-kcas*IQ)-IR*kcas;
48  if(seme<=0) seme += IM;
49  if(jcas<=NTAB) iv[jcas]=seme;
50  }
51  iy = iv[1];
52  }
53  kcas=seme/IQ;
54  seme=IA*(seme-kcas*IQ)-IR*kcas;
55  if(seme<=0) seme+=IM;
56  jcas=1+iy/NDIV;
57  iy=iv[jcas];
58  iv[jcas]=seme;
59  double casuale=AM*iy < RNMX ? AM*iy : RNMX;
60  // fprintf(CONTROLLA,"%f\n",casuale);
61  return casuale;
62 }
63 #endif
64 double Matematica::RandDiscrProb(double *Prob,int NBin){
65  double Ran = Casuale();
66  double InvNBin = 1./(double)NBin;
67  int j = 0;
68  for(int i=0;i<NBin-1;i++){
69  if(Ran > Prob[i] && Ran <= Prob[i+1]){
70  j = i;
71  break;
72  }
73  }
74  double xi = Ran - floor(Ran*NBin)*InvNBin;
75  double m = (Prob[j+1] - Prob[j])*NBin;
76  double y = j*InvNBin + m*xi;
77  return y;
78 }
79 #ifdef USE_FFTW
80 void Matematica::Spettro(double *st,double *sw,int NMax){//N=int^2
81  //printf("using fftw\n");
82  fftw_complex *uscita = (fftw_complex *) fftw_malloc( sizeof(fftw_complex)*(NMax) );
83  fftw_complex *entrata = (fftw_complex *) fftw_malloc( sizeof(fftw_complex)*(NMax) );
84  fftw_plan p;
85  for(int i=0;i<NMax;i++){
86  entrata[0][i] = st[i];
87  }
88  //for(int i=0;i<NMax;i++)printf("%d %lf \n",i,st[i]);
89  //p = fftw_plan_dft_r2c_1d(NMax,st,uscita,FFTW_FORWARD);
90  p = fftw_plan_dft_1d(NMax,entrata,uscita,FFTW_FORWARD,FFTW_ESTIMATE);
91  fftw_execute(p);
92  for(int i=0;i<NMax;i++){
93  sw[i] = QUAD(uscita[i][0])+QUAD(uscita[i][1]);
94  }
95  //for(int i=0;i<NMax;i++) printf("%d %lf %lf\n",i,st[i],sw[i]);
96  fftw_destroy_plan(p);
97  fftw_free(entrata);
98  fftw_free(uscita);
99 }
100 void Matematica::Spettro2d(double *st,double *sw,int NMax){//N=int^2
101  //printf("using fftw\n");
102  fftw_complex *uscita = (fftw_complex *) fftw_malloc( sizeof(fftw_complex)*(NMax*NMax) );
103  fftw_plan p = fftw_plan_dft_r2c_2d(NMax,NMax,st,uscita,FFTW_MEASURE);
104  fftw_execute(p);
105  for(int i=0;i<NMax;i++){
106  sw[i] = 0.;
107  for(int j=0;j<NMax;j++){
108  sw[i*NMax+j] = SQR(uscita[i*NMax+j][0])+SQR(uscita[i*NMax+j][0]);
109  }
110  }
111  fftw_destroy_plan(p);
112  fftw_free(uscita);
113 }
114 void Matematica::Spettro2d(double *st,double **sw,int NMax){//N=int^2
115  //printf("using fftw\n");
116  double NMaxInv = 1./(double)NMax;
117  int NHalf = (int)(NMax/2.);
118  fftw_plan p;
119  //fftw_complex *uscita = (fftw_complex *) fftw_malloc( sizeof(fftw_complex)*(NMax*NMax) );
120  fftw_complex *uscita = (fftw_complex *)calloc(NMax*NMax,sizeof(fftw_complex));
121  //for(int i=0;i<NMax;i++)for(int j=0;j<NMax;j++)printf("%d %d %lf \n",i,j,st[i*NMax+j]);
122  p = fftw_plan_dft_r2c_2d(NMax,NMax,st,uscita,FFTW_MEASURE);
123  fftw_execute(p);
124  for(int i=0;i<NMax;i++){
125  for(int j=0;j<NMax;j++){
126  sw[i][j] = QUAD(uscita[0][i*NMax+j])+QUAD(uscita[1][i*NMax+j]);
127  sw[i][j] *= QUAD(NMaxInv);
128  //printf("%d %d %lf %lf %lf\n",i,j,st[i*NMax+j],uscita[0][i*NMax+j],uscita[1][i*NMax+j]);
129  }
130  }
131  fftw_destroy_plan(p);
132  //fftw_free(uscita);
133  free(uscita);
134 }
135 #else
136 void Matematica::Spettro(double *st,double *sw,int NMax){//N=int^2
137  double Re1=0.,Im1=0.,Re2=0.,Im2=0.;
138  double FMass=(double)NMax/DUE_PI;
139  double dNMax = 1./(double)NMax;
140  for(int j=0;j<NMax;j++){
141  Re1=0.;Re2=0.;Im1=0.;Im2=0.;
142  for(int i=0;i<NMax;i++){
143  Re1+=st[i]*cos(DUE_PI*dNMax*(i)*j);
144  Im1-=st[i]*sin(DUE_PI*dNMax*(i)*j);
145  }
146  sw[j]= sqrt(QUAD((Re1+Re2))*dNMax+QUAD((Im1+Im2))*dNMax);
147  }
148  // double Parseval1 = 0.;
149  // double Parseval2 = 0.;
150  // for(int j=0;j<NMax;j++){
151  // Parseval1 += QUAD(st[j]);
152  // Parseval2 += sw[j];
153  // }
154  // printf("Parseval %lf=%lf\n",Parseval1,Parseval2);
155 }
156 void Matematica::Spettro2d(double *st,double *sw,int NMax){
157  printf("dft\n");
158  double dNMax = 1./(double)NMax;
159  int NHalf = (int)(NMax/2.);
160  for(int kx=-NMax/2;kx<NMax/2;kx++){
161  double qx = kx*dNMax;
162  for(int ky=-NMax/2;ky<NMax/2;ky++){
163  double qy = ky*dNMax;
164  double Re2=0.,Im2=0.;
165  double Re1=0.,Im1=0.;
166  for(int lx=0;lx<NMax;lx++){
167  for(int ly=0;ly<NMax;ly++){
168  double Arg = 2.*M_PI*(lx*kx + ly*ky);
169  double cy = cos(Arg);
170  double sy = sin(Arg);
171  Re1 += st[lx*NMax+ly]*cy;
172  Im1 += st[lx*NMax+ly]*sy;
173  }
174  }
175  int kkx = kx + NMax/2;
176  int kky = ky + NMax/2;
177  sw[kkx*NMax+kky] = SQR(Re1*dNMax) + SQR(Im1*dNMax);
178  }
179  }
180 }
181 // void Matematica::Spettro2d(double *st,double *sw,int NMax){
182 // double dNMax = 1./(double)NMax;
183 // int NHalf = (int)(NMax/2.);
184 // for(int kx=0;kx<NMax;kx++){
185 // for(int ky=0;ky<NMax;ky++){
186 // double Re2=0.,Im2=0.;
187 // for(int lx=0;lx<NMax;lx++){
188 // double cx = cos(kx*lx*dNMax*DUE_PI);
189 // double sx = sin(kx*lx*dNMax*DUE_PI);
190 // double Re1=0.,Im1=0.;
191 // for(int ly=0;ly<NMax;ly++){
192 // double cy = cos(ky*ly*dNMax*DUE_PI);
193 // double sy = sin(ky*ly*dNMax*DUE_PI);
194 // Re1 += st[lx*NMax + ly]*cy;
195 // Im1 += st[lx*NMax + ly]*sy;
196 // }
197 // Re2 += cx*Re1 - sx*Im1;
198 // Im2 -= sx*Re1 + cx*Im1;
199 // }
200 // sw[kx*NMax+ky] = SQR(Re2*dNMax) + SQR(Im2*dNMax);
201 // }
202 // }
203 // }
204 void Matematica::Spettro2d(double *st,double **sw,int NMax){
205  printf("No fftw\n");
206  double dNMax = 1./(double)NMax;
207  int NHalf = (int)(NMax/2.);
208  for(int l=0;l<NMax;l++){
209  for(int k=0;k<NMax;k++){
210  double Re1=0.,Im1=0.,Re2=0.,Im2=0.;
211  for(int n=0;n<NMax;n++){
212  double Cosn = cos(DUE_PI*l*n*dNMax);
213  double Sinn = sin(DUE_PI*l*n*dNMax);
214  for(int m=0;m<NMax;m++){
215  double Arg = DUE_PI*((l-NHalf)*n*dNMax + (k-NHalf)*m*dNMax);
216  // Re1 += cos(Arg)*st[n*NMax + m];
217  // Im2 += sin(Arg)*st[n*NMax + m];
218  Re1 += Cosn*cos(DUE_PI*k*m*dNMax)*st[n*NMax + m];
219  Re2 -= Sinn*sin(DUE_PI*k*m*dNMax)*st[n*NMax + m];
220  Im1 += Cosn*sin(DUE_PI*k*m*dNMax)*st[n*NMax + m];
221  Im2 += Sinn*cos(DUE_PI*k*m*dNMax)*st[n*NMax + m];
222  }
223  }
224  sw[l][k] = QUAD((Re1+Re2)*dNMax*dNMax) + QUAD((Im1+Im2)*dNMax*dNMax);
225  }
226  }
227  // double Parseval1 = 0.;
228  // double Parseval2 = 0.;
229  // for(int v=0;v<NMax;v++){
230  // for(int vv=0;vv<NMax;vv++){
231  // Parseval1 += QUAD(sw[v][vv]);
232  // Parseval2 += QUAD(st[v*NMax+vv]);
233  // }
234  // }
235  // printf("Parseval %lf=%lf\n",Parseval1,Parseval2);
236 }
237 #endif
238 void Matematica::Radice(double *st,double *sw,int N){
239  for(int j=0;j<N;j++){
240  if(st[j] > 0.)
241  sw[j] = sqrt(st[j]);
242  else
243  sw[j] = sqrt(-st[j]);
244  }
245 }
246 void Matematica::Autocor(double *st,double *sAutocor,int NMax){
247  for(int i=0;i<NMax;i++){
248  sAutocor[i] = 0.;
249  for(int j=0,k=0;j<NMax;j++){
250  k = i-j;
251  if(k<0)
252  k = NMax - i + j;
253  sAutocor[i] += (double)(st[k]*st[j]);
254  }
255  sAutocor[i] /= (double)(NMax);
256  }
257 }
258 void Matematica::MediaMobile(double *st,int NMax,double *sw,int Parti){
259  if(Parti <= 1){
260  return;
261  }
262  double PartiInv = 1./(double)Parti;
263  int Nw = (int) (NMax/(double)Parti);
264  for(int i=0;i<Nw;i++)
265  sw[i] = 0.;
266  for(int i=0,j=0;i<NMax;i++){
267  sw[j] += st[i];
268  if( (i%Parti)==0 ){
269  sw[j] *= PartiInv;
270  j++;
271  }
272  }
273 }
274 int Matematica::MediaMobile(double *st,int NMax,double *sw,double *sErr,int NParti){
275  if(NParti <= 0 ){
276  return NParti;
277  }
278  int Nw = (int) (NMax/(double)NParti);
279  double InvNParti = 1./(double)NParti;
280  for(int i=0;i<Nw;i++){
281  sw[i] = 0.;
282  sErr[i] = 0.;
283  }
284  for(int i=0,j=0;i<NMax;i++){
285  sw[j] += st[i];
286  sErr[j] += SQR((st[i]));
287  if( ((i+1)%NParti)==0 ){
288  sw[j] *= InvNParti;
289  sErr[j] = sqrt( (sErr[j] - SQR(sw[j])*NParti)*InvNParti );
290  j++;
291  if(j==Nw) return Nw;
292  }
293  }
294  return Nw;
295 }
296 int Matematica::CorrelaDuePunti(double *st,int NMax,double *sw,int Punti){
297  int PuntiMass = NMax - Punti;
298  for(int i=0;i<PuntiMass;i++){
299  sw[i] = (st[i] + st[i+Punti]) / 2.;
300  }
301  return PuntiMass;
302 }
303 void Matematica::Autosimilarita(double *st,int NMax,double *sw,int Potenze){
304  for(int j=0;j<Potenze;j++){
305  sw[j] = 0.;
306  }
307  for(int i=0;i<NMax;i++){
308  double Temp = 1;
309  for(int j=0;j<Potenze;j++){
310  Temp *= st[i];
311  sw[j] += Temp;
312  }
313  }
314  for(int j = 0;j<Potenze;j++){
315  sw[j] = log10(POS(sw[j]))/log10((double)NMax);
316  // printf("sw[%d] %f\n",j,sw[j]);
317  }
318 }
319 int Matematica::Normalizza(double *st,double *sw,int NMax){
320  double yMass=st[0];
321  double yMin = st[0];
322  double Media=0.;
323  int Zeri=0;
324  int Campioni=0;
325  for(int i = 1;i<NMax;i++){
326  if(yMin > st[i])
327  yMin = st[i];
328  if(yMass< st[i])
329  yMass = st[i];
330  Media += st[i];
331  }
332  Media /= (double)NMax;
333  for(int i=0;i<NMax;i++){
334  sw[i] = (st[i]-yMin)/(yMass-yMin);
335  if(st[i] > Media && st[i-1] < Media)
336  Zeri++;
337  }
338  Campioni = (int)(NMax/(double)Zeri);
339  return Campioni;
340 }
341 int Matematica::Normalizza(double *st,int NMax){
342  double yMass=st[0];
343  double yMin = st[0];
344  double Media=0.;
345  int Zeri=0;
346  int Campioni=0;
347  for(int i = 1;i<NMax;i++){
348  if(yMin > st[i])
349  yMin = st[i];
350  if(yMass< st[i])
351  yMass = st[i];
352  Media += st[i];
353  }
354  Media /= (double)NMax;
355  for(int i=0;i<NMax;i++){
356  st[i] = (st[i]-yMin)/(yMass-yMin);
357  if(st[i] > Media && st[i-1] < Media)
358  Zeri++;
359  }
360  Campioni = (int)(NMax/(double)Zeri);
361  return Campioni;
362 }
363 int Matematica::NormalizeArea(double *st,int NMax){
364  double Area = 0.;
365  for(int i = 0;i<NMax;i++){
366  Area += st[i];
367  }
368  for(int i=0;i<NMax;i++){
369  st[i] /= Area;
370  }
371  return 0;
372 }
373 MOMENTI Matematica::Distribuzione(const double *st,int NMax){
374  MOMENTI m1; m1.Uno=0.; m1.Due=0.; m1.Tre=0.;m1.Delta=0.;m1.Num=0;
375  m1.Min = st[0];
376  m1.Max = st[0];
377  for(int i=0;i<NMax;i++){
378  m1.Uno += st[i];
379  if(st[i] < m1.Min) m1.Min = st[i];
380  if(st[i] > m1.Max) m1.Max = st[i];
381  }
382  m1.Uno /= (double)(NMax);
383  m1.Delta = (m1.Max-m1.Min);
384  for(int i=0;i<NMax;i++){
385  m1.Due+=QUAD((st[i]-m1.Uno));
386  m1.Tre+=QUAD((st[i]-m1.Uno))*(st[i]-m1.Uno);
387  m1.Num++;
388  }
389  m1.Due=sqrt(m1.Due/(double)(NMax-1));
390  m1.Tre = pow(m1.Tre / (double)(NMax - 2),.33333);
391  return m1;
392 }
393 MOMENTI Matematica::DistrErr(const double *st,int NMax,double *Distr,double *Err,int NBin,double *Confine,int IfNorm){
394  MOMENTI m1 = Distribuzione(st,NMax,Distr,NBin,Confine,IfNorm);
395  for(int i=0;i<NBin;i++)
396  Err[i] = 0.;
397  double Add = 1./(double)NMax;
398  for(int i=0;i<NMax;i++){
399  int j = (int)((st[i]-Confine[0])/(Confine[1]-Confine[0])*NBin);
400  if(j >= NBin || j < 0) continue;
401  Err[j] += Add;
402  }
403  for(int i=0;i<NBin;i++){
404  Err[i] *= Distr[i];
405  }
406  return m1;
407 }
408 MOMENTI Matematica::Distribuzione(const double *st,int NMax,double *Distr,int NBin,double *Confine,int IfNorm){
409  MOMENTI m1; m1.Uno=0.; m1.Due=0.; m1.Tre=0.;m1.Delta=0.;m1.Num=0;
410  m1.Min = Confine[0];
411  m1.Max = Confine[1];
412  for(int i=0;i<NBin;i++) Distr[i]=0.;
413  for(int i=0;i<NMax;i++) m1.Uno += st[i];
414  m1.Uno/=(double)(NMax);
415  m1.Delta = (m1.Max - m1.Min)/(double)NBin;
416  double Add = !IfNorm ? 1. : 1./(double)NMax;
417  for(int i=0;i<NMax;i++){
418  m1.Due += QUAD(st[i]-m1.Uno);
419  m1.Tre += m1.Due*(st[i]-m1.Uno);
420  int j = (int)((st[i]-Confine[0])/(Confine[1]-Confine[0])*NBin);
421  if(j >= NBin || j < 0) continue;
422  Distr[j] += Add;
423  m1.Num++;
424  }
425  m1.yMin = Distr[0];
426  m1.yMax = Distr[0];
427  for(int v=0;v<NBin;v++){
428  if(m1.yMin > Distr[v]) m1.yMin = Distr[v];
429  if(m1.yMax < Distr[v]) m1.yMax = Distr[v];
430  }
431  m1.Due=sqrt(m1.Due/((double)(NMax-1)));
432  m1.Tre = pow(m1.Tre,.33333) / (double)(NMax - 2);
433  // printf("[-] Massimo %g m1.Minimo %g Media %g [-] \n[-] Scarto %g Terzo %g Deltax/3 %g [-]\n",m1.Massimo,m1.Minimo,m1.Uno,m1.Due,m1.Tre,(m1.Massimo-m1.Minimo)/6.);
434  return m1;
435 }
436 MOMENTI Matematica::Distribuzione(const double *st,int NMax,double *Distr,int NBin,int IfNorm){
437  MOMENTI m1; m1.Uno=0.; m1.Due=0.; m1.Tre=0.;m1.Delta=0.;m1.Num=0;
438  m1.Min = st[0];
439  m1.Max = st[0];
440  double Sum2 = 0.;
441  double Sum3 = 0.;
442  for(int i=0;i<NBin;i++) Distr[i] = 0.;
443  for(int i=0;i<NMax;i++){
444  if(st[i] < m1.Min) m1.Min = st[i];
445  if(st[i] > m1.Max) m1.Max = st[i];
446  m1.Uno += st[i];
447  Sum2 += QUAD(st[i]);
448  Sum3 += st[i]*st[i]*st[i];
449  }
450  double Add = !IfNorm ? 1. : 1./(double)NMax;
451  m1.Delta = (m1.Max-m1.Min)/(double)NBin;
452  for(int i=0;i<NMax;i++){
453  int j = (int)((st[i]-m1.Min)/(m1.Max-m1.Min)*NBin);
454  if(j >= NBin) continue;
455  Distr[j] += Add;
456  m1.Num++;
457  }
458  m1.yMin = Distr[0];
459  m1.yMax = Distr[0];
460  for(int v=0;v<NBin;v++){
461  if(m1.yMin > Distr[v]) m1.yMin = Distr[v];
462  if(m1.yMax < Distr[v]) m1.yMax = Distr[v];
463  }
464  m1.Uno/=(double)(NMax);
465  m1.Due = sqrt((Sum2 - QUAD(m1.Uno)*NMax)/((double)NMax-1.));
466  m1.Tre = Sum3-3.*Sum2*m1.Uno+3.*NMax*QUAD(m1.Uno)-NMax*m1.Uno;
467  m1.Tre = pow( m1.Tre , 1./3.)/(double)(NMax-1);
468  // printf("[-] Massimo %g m1.Minimo %g Media %g [-] \n[-] Scarto %g Terzo %g Deltax/3 %g [-]\n",m1.Massimo,m1.Minimo,m1.Uno,m1.Due,m1.Tre,(m1.Massimo-m1.Minimo)/6.);
469  return m1;
470 }
471 void Matematica::DistrSample(double *Px,double *Py,int NMax,double **Distr,int NBin,const int NSample,int IfNorm,double *yBorder){
472  double xBorder[2];
473  double yDelta;
474  double Norm[NSample];
475  xBorder[1] = Px[0];
476  xBorder[0] = Px[0];
477  for(int n=0;n<NMax;n++){
478  if(xBorder[1] < Px[n]) xBorder[1] = Px[n];
479  if(xBorder[0] > Px[n]) xBorder[0] = Px[n];
480  }
481  double Dx = (xBorder[1] - xBorder[0]) > 0. ? 1./(xBorder[1] - xBorder[0]) : 1.;
482  if(1==0){//search for borders
483  for(int s=0;s<NSample;s++){
484  yBorder[0] = Py[s];
485  yBorder[1] = Py[s];
486  }
487  for(int n=0;n<NMax;n++){
488  int sx = (int)((Px[n]-xBorder[0])*Dx*NSample);
489  if(sx < 0 || sx >= NSample) continue;
490  if(yBorder[0] > Py[n]) yBorder[0] = Py[n];
491  if(yBorder[1] < Py[n]) yBorder[1] = Py[n];
492  }
493  }
494  yDelta = yBorder[1] - yBorder[0] > 0. ? 1./(yBorder[1] - yBorder[0]) : 1.;
495  for(int n=0;n<NMax;n++){
496  int sx = (int)((Px[n]-xBorder[0])*Dx*NSample);
497  if(sx < 0 || sx >= NSample) continue;
498  int by = (int)((Py[n]-yBorder[0])*yDelta*NBin);
499  if(by < 0 || by >= NBin) continue;
500  Distr[sx][by] += 1.;
501  }
502  if(IfNorm){
503  for(int sx=0;sx<NSample;sx++){
504  Norm[sx] = 0.;
505  for(int by=0;by<NBin;by++){
506  Norm[sx] += Distr[sx][by];
507  }
508  }
509  for(int sx=0;sx<NSample;sx++){
510  Norm[sx] = Norm[sx] > 1. ? 1./Norm[sx] : 1.;
511  for(int by=0;by<NBin;by++){
512  Distr[sx][by] *= Norm[sx];
513  }
514  }
515  }
516 }
517 MOMENTI Matematica::DistribuzioneGauss(const double *st,int NMax,double *Distr,double *dInt,int NBin,int IfNorm){
518  MOMENTI m1;
519  m1 = Distribuzione(st,NMax,Distr,NBin,IfNorm);
520  m1.Chi = 0.;
521  for(int v=0;v<NBin;v++){
522  dInt[v] = NMax*Gauss( m1.Uno , m1.Due , (v)*m1.Delta + m1.Min)*m1.Delta;
523  m1.Chi += QUAD( (Distr[v] - dInt[v])/m1.Due);
524  }
525  m1.Chi /= NBin-2;
526  return m1;
527 }
528 MOMENTI Matematica::DistribuzioneMaxwell(const double *st,int NMax,double *Distr,double *dInt,int NBin,int IfNorm){
529  MOMENTI m1;
530  m1 = Distribuzione(st,NMax,Distr,NBin,IfNorm);
531  m1.Chi = 0.;
532  for(int i=0;i<NBin;i++){
533  dInt[i] = 2*DUE_PI*QUAD((i*m1.Delta + m1.Min))*NMax*Gauss( m1.Uno , m1.Due , i*m1.Delta + m1.Min)*m1.Delta;
534  m1.Chi += QUAD( Distr[i] - dInt[i]);
535  }
536  m1.Chi /= NBin-2;
537  return m1;
538 }
539 MOMENTI Matematica::WeightAverage(const double *sx,const double *sy,int NMax){
540  MOMENTI m1; m1.Uno=0.; m1.Due=0.; m1.Tre=0.;m1.Delta=0.;m1.Num=0;
541  m1.Min = sx[0];
542  m1.Max = sx[0];
543  double TotWei = 0.;
544  for(int i=0;i<NMax;i++){
545  if(sx[i] < m1.Min) m1.Min = sx[i];
546  if(sx[i] > m1.Max) m1.Max = sx[i];
547  m1.Uno += sx[i]*sy[i];
548  TotWei += sy[i];
549  m1.Due += sy[i] > 0. ? 1./sy[i] : 0.;
550  }
551  m1.Uno /= TotWei;
552  m1.Due = sqrt(1./m1.Due);
553  m1.Delta = (m1.Max-m1.Min);
554  return m1;
555 }
556 void Matematica::WeightHisto(double **hist,double *Border,int NBin,int NHisto,double tolerance,double *OrPos,double *kSpring){
557  if(1==0) {
558  cout<<"Usage:\n\n";
559  cout<<"do_wham_c <TEMPERATURE> <MIN> <MAX> <BINS> <TOLERANCE> <OUTPUT> <INPUT1> <INPUT2> ....\n\n";
560  cout<<"TEMPERATURE\t\tTemperature that simulations were performed at\n";
561  cout<<"MIN\t\t\tMinimum value used in histogram. All histograms must have been computed with matching value.\n";
562  cout<<"MAX\t\t\tMaximum value used in histogram. All histograms must have been computed with matching value.\n";
563  cout<<"NBIN\t\t\tNumber of NBin in histogram. All histograms must have been computed with matching value.\n";
564  cout<<"TOLERANCE\t\tTolerance value used to determine convergence. Value of 0.00001 seems to work fine.\n";
565  cout<<"OUTPUT\t\t\tName of output file.\n";
566  cout<<"INPUT\t\t\tNames of input files. These are histogram files produced from the .pdo files using make_histo.pl.\n";
567  cout<<endl;
568  }
569  double temperature = 1.;
570  double min = Border[0];
571  double max = Border[1];
572  double RT = temperature * 8.314472e-3;
573  ofstream output("WeightHistoAnal.dat");
574  double* F = new double[NHisto];
575  double* old_F = new double[NHisto];
576  double* result = new double[NBin];
577  double* N = new double[NHisto];
578  double* energy = new double[NBin];
579  double * numerator = new double[NBin];
580  double ** expU = new double *[NHisto];
581  double ** NexpU = new double *[NHisto];
582  //filling the exponential and normalizing factor
583  for(int h = 0;h < NHisto;++h) {
584  F[h] = 1.0;
585  old_F[h] = 0.0;
586  N[h] = 0.0;
587  expU[h] = new double[NBin];
588  NexpU[h] = new double[NBin];
589  for(int i = 0; i < NBin; ++i) {
590  N[h] += hist[h][i];
591  double pos = (i + 0.5)/NBin * (max - min) + min;
592  expU[h][i] = exp(-(0.5 * kSpring[h] * (pos - OrPos[h]) * (pos - OrPos[h]))/RT);
593  }
594  for(int i=0; i < NBin; ++i) {
595  NexpU[h][i] = N[h] * expU[h][i];
596  }
597  }
598  for(int bin = 0; bin < NBin; ++bin) {
599  result[bin] = 0.0;
600  }
601  double d_max = 0.;
602  for(int bin = 0; bin < NBin; ++bin) {
603  numerator[bin] = 0.0;
604  for(int file = 0; file < NHisto; ++file) {
605  numerator[bin] += hist[file][bin];
606  }
607  }
608  cerr<<"Done loading histograms. Doing WHAM."<<endl;
609  int count = 0;
610  int IfContinue = 1;
611  do {
612  d_max = 0.;
613  for(int bin=0; bin < NBin; ++bin) {
614  double denom = 0.0;
615  for(int file = 0; file < NHisto; ++file) {
616  denom += NexpU[file][bin] / F[file];
617  }
618  //printf("%lf %lf %lf\n",denom,numerator[bin],result[bin]);
619  if(denom > 1.e9) IfContinue = 0;
620  denom = denom <= 0.0000001 ? 1. : denom;
621  result[bin] = numerator[bin] / denom;
622  }
623  double norm = 0.0;
624  for(int bin=0; bin < NBin; ++bin) norm += result[bin];
625  norm = norm > 0. ? norm : 1.;
626  for(int bin=0; bin < NBin; ++bin) result[bin] /= norm;
627  for(int file = 0; file < NHisto; ++file) {
628  double Z = 0.;
629  for (int bin = 0; bin < NBin; ++bin){
630  Z += result[bin] * expU[file][bin];
631  }
632  if(isnan(Z)){IfContinue = 0;break;}
633  if(isinf(Z)){IfContinue = 0;break;}
634  if(Z <= 1.e-9){IfContinue = 0;break;}
635  double dZ = (max - min)/(double)NBin;
636  Z = Z * dZ;
637  old_F[file] = F[file];
638  //printf("%d %d %lf %lf\n",count,file,F[file],Z);
639  F[file] = Z;
640  }
641  if(!IfContinue) break;
642  double temp = F[0];
643  for(int file = 0; file < NHisto; ++file) {
644  temp = RT * fabs(log(F[file]) - log(old_F[file]));
645  if(temp > d_max)
646  d_max = temp;
647  }
648  if(!(count % 1000))
649  cerr<<"Step: "<<count<<"\tTolerance: "<<d_max<<endl;
650  ++count;
651  if(count == 2) break;
652  } while(d_max > tolerance);
653  for(int bin = 0; bin < NBin; ++bin) {
654  energy[bin] = 0.0;
655  }
656  for(int bin = 0; bin < NBin; ++bin) {
657  double res = result[bin] > 0. ? log(result[bin]) : 0.;
658  energy[bin] = -RT*res;
659  double pos = (double)bin / NBin * (max - min) + min;
660  output << pos << " " << result[bin] << " " << energy[bin] << endl;
661  }
662  output.close();
663  delete [] F;
664  delete [] old_F;
665  delete [] result;
666  delete [] N;
667  delete [] energy;
668  delete [] numerator;
669  for(int current = 0;current < NHisto;++current) {
670  delete [] expU[current];
671  delete [] NexpU[current];
672  }
673  delete [] expU;
674  delete [] NexpU;
675 }
676 void Matematica::Sort(double *Sign,int NMax){
677  for(int i=0;i<NMax;i++){
678  for(int j=i;j>0;j--){
679  if(Sign[j] < Sign[j-1])
680  Swap(j,j-1,Sign);
681  else
682  break;
683  }
684  }
685 }
686 void Matematica::Swap(int i, int j,double *Sign){
687  double Temp = Sign[j];
688  Sign[j] = Sign[i];
689  Sign[i] = Temp;
690 }
691 void Matematica::Swap(double *s,int si,double *t,int ti,const int NDim){
692  double Temp[NDim];
693  for(int d=0;d<NDim;d++){
694  Temp[d] = s[si*NDim+d];
695  s[si*NDim+d] = t[ti*NDim+d];
696  t[ti*NDim+d] = Temp[d];
697  }
698 }
699 void Matematica::Sort(int *Sign,int NMax){
700  for(int i=0;i<NMax;i++){
701  for(int j=i;j>0;j--){
702  if(Sign[j] < Sign[j-1])
703  Swap(j,j-1,Sign);
704  else
705  break;
706  }
707  }
708 }
709 void Matematica::Swap(int i, int j,int *Sign){
710  int Temp = Sign[j];
711  Sign[j] = Sign[i];
712  Sign[i] = Temp;
713 }
714 void Matematica::ConvWeight(double *st,int NMax,double *sw,int *WIndex,int NWeight){
715  double *st1 = (double *)calloc(NMax,sizeof(double));
716  for(int n=0;n<NMax;n++){
717  st1[n] = st[n];
718  }
719  for(int n=0;n<NMax;n++){
720  st[n] = 0;
721  if(n+WIndex[0] < 0 || n+WIndex[NWeight-1] >= NMax){
722  st[n] = st1[n];
723  continue;
724  }
725  for(int w=0;w<NWeight;w++){
726  st[n] += sw[w]*st1[n+WIndex[w]];
727  }
728  }
729  free(st1);
730 }
731 void Matematica::FillWeightGauss(double *st,int *WIndex,int NWeight,double CutOff,double Sigma){
732  int Half = NWeight/2;
733  double Norm = 0.;
734  for(int w=0;w<NWeight;w++){
735  WIndex[w] = w-Half;
736  double x = CutOff*(w-Half)/(double)NWeight;
737  double r2 = SQR(x);
738  double Gauss = exp(-r2*.5/SQR(Sigma));
739  st[w] = Gauss;
740  Norm += Gauss;
741  }
742  for(int w=0;w<NWeight;w++){
743  st[w] /= Norm;
744  }
745 }
746 void Matematica::ExecCommand(double *st,double *sw,int NMax,char *cmd){
747  for(int i=0;i<NMax;i++){
748  sw[i] = ExecFormula(st[i],1.,cmd);
749  }
750 }
751 double Matematica::ExecFormula(double x,double y,char *cmd){
752  if(strlen(cmd) == 0) return x;
753  const int NOp = 20;
754  int PosOp[NOp];
755  int PosPar[NOp];
756  int PosNum[2*NOp];
757  PosNum[0] = 0;
758  int nOp = 1;
759  int Comm[NOp];
760  double Numb;
761  char *cPart = cmd;
762  for(int i=0;i<NOp;i++){
763  cPart = strpbrk(cPart+1,"*/^+-");
764  if(cPart == NULL) break;
765  Comm[i] = (int)*cPart;
766  PosOp[i] = cPart-cmd+1;
767  PosNum[nOp] = PosOp[i];
768  nOp++;
769  }
770  Sort(PosNum,nOp);
771  PosNum[nOp++] = strlen(cmd);
772  double Res = 0.;
773  for(int i=0;i<nOp-1;i++){
774  if( *(cmd+PosNum[i]) == 'x'){
775  Numb = x;
776  }
777  if( *(cmd+PosNum[i]) == 'y'){
778  Numb = y;
779  }
780  else
781  sscanf(cmd+PosNum[i],"%lf",&Numb);
782  if(i == 0){
783  Res = Numb;
784  continue;
785  }
786  switch(Comm[i-1]){
787  case '+':
788  Res += Numb;
789  break;
790  case '-':
791  Res -= Numb;
792  break;
793  case '/':
794  Res /= Numb;
795  break;
796  case '*':
797  Res *= Numb;
798  break;
799  }
800  //printf("%d) %s %d %lf %lf %lf %c %lf\n",i,cmd+PosNum[i],PosNum[i],Numb,x,y,Comm[i],Res);
801  }
802  return Res;
803  // cPart = cmd;
804  // for(int i=0;i<NOp;i++){
805  // cPart = strpbrk(cPart+1,"()");
806  // if(cPart == NULL) break;
807  // PosPar[i] = cPart-cmd+1;
808  // PosNum[nOp] = PosOp[i];
809  // nOp++;
810  // }
811  //
812 }
813 double Matematica::ExecFormula(double **st,int n,char *cmd){
814  if(strlen(cmd) == 0) return st[0][n];
815  const int NOp = 20;
816  int PosOp[NOp];
817  int PosPar[NOp];
818  int PosNum[2*NOp];
819  PosNum[0] = 0;
820  int nOp = 1;
821  int Comm[NOp];
822  double Numb;
823  char *cPart = cmd;
824  for(int i=0;i<NOp;i++){
825  cPart = strpbrk(cPart+1,"*/^+-");
826  if(cPart == NULL) break;
827  Comm[i] = (int)*cPart;
828  PosOp[i] = cPart-cmd+1;
829  PosNum[nOp] = PosOp[i];
830  nOp++;
831  }
832  Sort(PosNum,nOp);
833  PosNum[nOp++] = strlen(cmd);
834  double Res = 0.;
835  int NVar = 0;
836  for(int i=0;i<nOp-1;i++){
837  if( *(cmd+PosNum[i]) == 'x'){
838  sscanf(cmd+PosNum[i]+1,"%d",&NVar);
839  Numb = st[NVar][n];
840  }
841  else
842  sscanf(cmd+PosNum[i],"%lf",&Numb);
843  //printf("%d) %s %d %lf %lf %c\n",i,cmd+PosNum[i],PosNum[i],Numb,x,Comm[i]);
844  if(i == 0){
845  Res = Numb;
846  continue;
847  }
848  switch(Comm[i-1]){
849  case '+':
850  Res += Numb;
851  break;
852  case '-':
853  Res -= Numb;
854  break;
855  case '/':
856  Res /= Numb;
857  break;
858  case '*':
859  Res *= Numb;
860  break;
861  }
862  }
863  return Res;
864 }
double RandDiscrProb(double *Prob, int NBin)
Random number following a discrete probability.
void Radice(double *st, double *sw, int N)
Compute the root of the signal.
void FillWeightGauss(double *st, int *WIndex, int NWeight, double CutOff, double Sigma)
Fill the weight array with a gaussian fuction.
int Zeri(double a, double b, double *Radici, int NRadici)
Find the.
double Min
Minimum value.
MOMENTI DistribuzioneMaxwell(const double *st, int NMass, double *Intervalli, double *dInt, int Valori, int IfNorm)
Look for the Maxwellian distribution.
bool InizializzaGaussiano(double Scarto, int N)
Initialize the Gaussian number generator.
void ConvWeight(double *st, int NMax, double *sw, int *WIndex, int NWeight)
Convolute with a weight.
double Due
Second moment.
double Tre
Third moment.
double yMin
Minimum of the ordinate.
MOMENTI Distribuzione(const double *st, int NMass)
Moments of a signal.
MOMENTI WeightAverage(const double *sx, const double *sy, int NMax)
Calculate the weighted average.
double yMax
Maximum of the ordinate.
Moments of a distribution.
void Sort(double *Sign, int NMass)
Sort.
void WeightHisto(double **hist, double *Border, int NBin, int NHisto, double tolerance, double *OrPos, double *kSpring)
Weighted histogram analysis.
MOMENTI DistrErr(const double *st, int NMass, double *Intervalli, double *Err, int Valori, double *Confine, int IfNorm)
Moments and histogram of a signal between two values.
double Casuale()
Random uniform number.
int NormalizeArea(double *st, int NMass)
Normalize.
void Swap(int i, int j, double *Sign)
Swap to indices.
void DistrSample(double *Px, double *Py, int NMax, double **Distr, int NBin, const int NSample, int IfNorm, double *xBound)
Compare the distribution of a sample of data.
MOMENTI DistribuzioneGauss(const double *st, int NMass, double *Intervalli, double *dInt, int Valori, int IfNorm)
Look for the Gaussian distribution.
int CorrelaDuePunti(double *st, int NMass, double *sw, int Punti)
Two points correlation.
void Spettro2d(double *st, double *sw, int NMass)
Compute the 2d spectrum of.
double Gauss(double Media, double Scarto, double x)
Gaussian.
double Gaussiano(double Media, double Scarto)
Gaussian random number.
double Max
Maximum value.
double F(double TD, double T)
Boh.
double Uno
First moment.
double ExecFormula(double x, double y, char *cmd)
Execute a formula.
double Chi
Chi square.
double Norm(double *st, int NMass)
Norm of an array.
int Num
Number of points considered.
void Autocor(bool *st, double *sAuto, int N)
Compute the autocorrelation of a boolean signal.
void ExecCommand(double *st, double *st1, int NMass, char *cmd)
Execute a command defined in string.
int Normalizza(double *st, int NMass)
Normalize.
double Delta
Step interval.
void Spettro(double *st, double *sw, int NMass)
Compute the spectrum.
void Autosimilarita(double *st, int NMass, double *sw, int Valori)
Self similarity.
void MediaMobile(double *st, int NMass, double *sw, int Parti)
Running average.