TheBoussinesqModel  3.2.1
 All Data Structures Files Functions Variables Typedefs Macros Pages
datamanipulation.c
Go to the documentation of this file.
1 #include "turtle.h"
2 #include "t_datamanipulation.h"
3 #include "t_alloc.h"
4 //#include "t_statistics.h"
5 #include "write_dem.h"
6 
7 /*-------------------------------------------------------------------------------*/
8 
9 
10 
11 DOUBLEVECTOR *vectorize_doublematrix(DOUBLEMATRIX *input)
12 
13 {
14 
15  DOUBLEVECTOR *U;
16 
17  U=(DOUBLEVECTOR *)malloc(sizeof(DOUBLEVECTOR));
18  if(!U){
19  t_error("This vector cannot be allocated");
20  }
21 
22  U->nl=1;
23  U->nh=(input->nrh)*(input->nch);
24  U->isdynamic=1;
25  if(input->co==NULL || input->isdynamic!=1 ){
26 
27  t_error("This matrix is not properly allocated");
28 
29  } else {
30 
31  U->co=input->co[1];
32 
33  }
34 
35  /* print_doublevector_elements(U,10); */
36 
37 
38  free(input);
39 
40  return U;
41 
42 
43 }
44 
45 
46 /*-------------------------------------------------------------------------------*/
47 
48 
49 
50 LONGVECTOR *vectorize_longmatrix(LONGMATRIX *input)
51 
52 {
53 
54  LONGVECTOR *U;
55 
56  U=(LONGVECTOR *)malloc(sizeof(LONGVECTOR));
57  if(!U){
58  t_error("This vector cannot be allocated");
59  }
60 
61  U->nl=1;
62  U->nh=(input->nrh)*(input->nch);
63  U->isdynamic=1;
64  if(input->co==NULL || input->isdynamic!=1 ){
65 
66  t_error("This matrix is not properly allocated");
67 
68  } else {
69 
70  U->co=input->co[1];
71 
72  }
73 
74  /* print_doublevector_elements(U,10); */
75 
76 
77  free(input);
78 
79  return U;
80 
81 
82 }
83 /*-------------------------------------------------------------------------------*/
84 
85 
86 
87 FLOATVECTOR *vectorize_floatmatrix(FLOATMATRIX *input)
88 
89 {
90 
91  FLOATVECTOR *U;
92 
93  U=(FLOATVECTOR *)malloc(sizeof(FLOATVECTOR));
94  if(!U){
95  t_error("This vector cannot be allocated");
96  }
97 
98  U->nl=1;
99  U->nh=(input->nrh)*(input->nch);
100  U->isdynamic=1;
101  if(input->co==NULL || input->isdynamic!=1 ){
102 
103  t_error("This matrix is not properly allocated");
104 
105  } else {
106 
107  U->co=input->co[1];
108 
109  }
110 
111  /* print_floatvector_elements(U,10); */
112 
113 
114  free(input);
115 
116  return U;
117 
118 
119 }
120 
121 
122 /*-------------------------------------------------------------------------------*/
123 SHORTVECTOR *vectorize_shortmatrix(SHORTMATRIX *input)
124 
125 {
126 
127  SHORTVECTOR *U;
128 
129  U=(SHORTVECTOR *)malloc(sizeof(SHORTVECTOR));
130  if(!U){
131  t_error("This vector cannot be allocated");
132  }
133 
134  U->nl=1;
135  U->nh=(input->nrh)*(input->nch);
136  U->isdynamic=1;
137  if(input->co==NULL || input->isdynamic!=1 ){
138 
139  t_error("This matrix is not properly allocated");
140 
141  } else {
142 
143  U->co=input->co[1];
144 
145  }
146 
147  /* print_floatvector_elements(U,10); */
148 
149 
150  free(input);
151 
152  return U;
153 
154 
155 }
156 
157 
158 /*-------------------------------------------------------------------------------*/
159 
160 void sortreal(DOUBLEVECTOR *ra)
161 
162 {
163 
164  long n;
165  unsigned int i,ir,j,l;
166  double rra;
167 
168  n=ra->nh;
169 
170  if (n <2) return;
171  l=(n>>1)+1;
172  ir=n;
173 
174  for(;;){
175  if (l >1) {
176  rra=ra->co[--l];
177  } else {
178  rra=ra->co[ir];
179  ra->co[ir]=ra->co[1];
180  if (--ir ==1){
181  ra->co[1]=rra;
182  break;
183  }
184  }
185  i=l;
186  j=l+l;
187  while (j<=ir){
188  if(j<ir && ra->co[j] < ra->co[j+1]) j++;
189  if(rra < ra->co[j]){
190  ra->co[i]=ra->co[j];
191  i=j;
192  j <<=1;
193  } else j=ir+1;
194  }
195  ra->co[i]=rra;
196  }
197 
198 }
199 
200 /*-------------------------------------------------------------------------------*/
201 void sort2realvectors(DOUBLEVECTOR *ra,DOUBLEVECTOR *rb)
202 
203 {
204 
205  long n;
206  unsigned int i,ir,j,l;
207  double rra,rrb;
208 
209  n=ra->nh;
210 
211  if (n <2) return;
212  l=(n>>1)+1;
213  ir=n;
214 
215  for(;;){
216  if (l >1) {
217  rra=ra->co[--l];
218  rrb=rb->co[l];
219  } else {
220  rra=ra->co[ir];
221  ra->co[ir]=ra->co[1];
222  rrb=rb->co[ir];
223  rb->co[ir]=rb->co[1];
224 
225  if (--ir ==1){
226  ra->co[1]=rra;
227  rb->co[1]=rrb;
228  break;
229  }
230  }
231  i=l;
232  j=l+l;
233  while (j<=ir){
234  if(j<ir && ra->co[j] < ra->co[j+1]) j++;
235  if(rra < ra->co[j]){
236  ra->co[i]=ra->co[j];
237  rb->co[i]=rb->co[j];
238  i=j;
239  j <<=1;
240  } else j=ir+1;
241  }
242  ra->co[i]=rra;
243  rb->co[i]=rrb;
244 
245  }
246 
247 
248 }
249 
250 /*--------------------------------------------------------------------------------*/
251 
252 void sort2floatvectors(FLOATVECTOR *ra,FLOATVECTOR *rb)
253 
254 {
255 
256  long n;
257  unsigned int i,ir,j,l;
258  float rra;
259  float rrb;
260 
261  n=ra->nh;
262 
263  if (n <2) return;
264  l=(n>>1)+1;
265  ir=n;
266  for(;;){
267  if (l >1) {
268  rra=ra->co[--l];
269  rrb=rb->co[l];
270  } else {
271  rra=ra->co[ir];
272  ra->co[ir]=ra->co[1];
273  rrb=rb->co[ir];
274  rb->co[ir]=rb->co[1];
275 
276  if (--ir ==1){
277  ra->co[1]=rra;
278  rb->co[1]=rrb;
279  break;
280  }
281  }
282  i=l;
283  j=l+l;
284  while (j<=ir){
285  if(j<ir && ra->co[j] < ra->co[j+1]) j++;
286  if(rra < ra->co[j]){
287  ra->co[i]=ra->co[j];
288  rb->co[i]=rb->co[j];
289  i=j;
290  j <<=1;
291  } else j=ir+1;
292  }
293  ra->co[i]=rra;
294  rb->co[i]=rrb;
295 
296  }
297 
298 }
299 
300 /*-------------------------------------------------------------------------------*/
301 void sort2vectors(LONGVECTOR *ra,FLOATVECTOR *rb)
302 
303 {
304 
305  long n;
306  unsigned int i,ir,j,l;
307  long rra;
308  float rrb;
309 
310  n=ra->nh;
311 
312  if (n <2) return;
313  l=(n>>1)+1;
314  ir=n;
315 
316  for(;;){
317  if (l >1) {
318  rra=ra->co[--l];
319  rrb=rb->co[l];
320  } else {
321  rra=ra->co[ir];
322  ra->co[ir]=ra->co[1];
323  rrb=rb->co[ir];
324  rb->co[ir]=rb->co[1];
325 
326  if (--ir ==1){
327  ra->co[1]=rra;
328  rb->co[1]=rrb;
329  break;
330  }
331  }
332  i=l;
333  j=l+l;
334  while (j<=ir){
335  if(j<ir && ra->co[j] < ra->co[j+1]) j++;
336  if(rra < ra->co[j]){
337  ra->co[i]=ra->co[j];
338  rb->co[i]=rb->co[j];
339  i=j;
340  j <<=1;
341  } else j=ir+1;
342  }
343  ra->co[i]=rra;
344  rb->co[i]=rrb;
345 
346  }
347 
348 
349 }
350 
351 
352 
353 /*-------------------------------------------------------------------------------*/
354 /*
355 Inputs: 1) a sorted vector containing the data to be binned; 2) the number of bins
356 (instead of fixing the size of each bin is usually convenient to select a fixed number
357 of bins ). If the number of bins is set to 0, is simply counted the number of elements with the
358 same abscissa; 3) the minimum number of elements required in each bin (this implies
359 that if a bin has not enough elements the numbers of elements of two adjacent bins are summed).
360 In the case of the exponential binning, the j-th bin extends from base^(j*delta)
361 to base^((j+1)*delta) where base is the fourth input field and delta the extension of
362 each bin in the logarithm axis.
363 
364 
365 Return: a matrix of double: the first column contains the number of elements
366 in each mean, the second coulumns the mean abscissa of the data in the bin, the third
367 the highest limit of each bin interval */
368 
369 DOUBLEMATRIX *simplehystogram(DOUBLEVECTOR *U,long N,long mn)
370 
371 {
372 
373  long j,count,count1;
374  double delta,min,max,tmp;
375  DOUBLEMATRIX *indx=NULL;
376  XYZ *bins,*head,*pre;
377  /* char ch; */
378 
379  if(N <0 ) t_error("A negative number of bins is not allowed");
380 
381 
382  head=new_xyz();
383  head->y=U->co[1];
384  pre=head;
385  bins=head;
386  count1=1;
387  count=2;
388  while(count <=U->nh){
389 
390  while(U->co[count]==U->co[count-1] && count <=U->nh){
391  count++;
392  }
393 
394  bins->next=new_xyz();
395  bins=bins->next;
396  bins->x=count-count1;
397  count1=count;
398  head->x++;
399  bins->y=U->co[count-1];
400  count++;
401 
402  }
403 
404 
405  bins=head->next;
406 
407  while(bins!=NULL){
408  count=bins->x;
409  while(bins->x < mn && bins->next!=NULL){
410  bins->y=(bins->y*bins->x+bins->next->y*bins->next->x)/(bins->x+bins->next->x);
411  bins->x+=bins->next->x;
412  bins->z=0;
413  delete_xyz(head,bins);
414  (head->x)--;
415  }
416 
417  if(bins->x< mn) {
418 
419  if(bins!=head){
420  pre->x+=bins->x;
421  delete_xyz(head,pre);
422  }else {
423  printf("\nWarning::Obtaining just one bin\n");
424  }
425 
426  }
427 
428  pre=bins;
429  bins=bins->next;
430  }
431 
432  /* print_xyz_elements(head,10); */
433 
434  if(N<1){
435 
436  indx=new_doublematrix(head->x,3);
437  xyz_into_doublematrix(head,indx);
438  delete_xyz_list(head);
439 
440  }else{
441 
442  if(head->x <N) printf("\nWarning::Some bin are remaining empty\n");
443  /* cnt=count_xyz_elements(head); */
444  indx=new_doublematrix(N,3);
445  delta=( U->co[U->nh]-head->next->y)/(N-1);
446  max=U->co[U->nh]+0.5*delta;
447  min=head->next->y-0.5*delta;
448  head->y=head->z=min;
449  bins=head->next;
450  for(j=1;j<=N;j++){
451  indx->co[j][1]=0;
452  indx->co[j][2]=0;
453  indx->co[j][3]=min+j*delta;
454  tmp=0;
455  if(bins!=NULL){
456  while(bins->y<=min+j*delta){
457  tmp+=bins->x;
458  indx->co[j][2]=(indx->co[j][1]*indx->co[j][2]+bins->y*bins->x)/tmp;
459  indx->co[j][1]=tmp;
460  bins=bins->next;
461  if(bins==NULL) break;
462  }
463  }
464 
465  }
466 
467  delete_xyz_list(head);
468 
469  }
470 
471  /* print_doublematrix_elements(indx,10); */
472 
473  return indx;
474 
475 }
476 
477 /*----------------------------------------------------------------------*/
478 
479 DOUBLEMATRIX *exponentialhystogram(DOUBLEVECTOR *U,long N,long mn,long base)
480 
481 {
482 
483  long j,count,count1;
484  double logdelta,min,max,tmp,logbase;
485  DOUBLEMATRIX *indx=NULL;
486  XYZ *bins,*head,*pre;
487 
488 
489  if(N <0 ) t_error("A negative number of bins is not allowed");
490 
491  logbase=log(base);
492  head=new_xyz();
493  head->y=U->co[1];
494  pre=head;
495  bins=head;
496  count1=1;
497  count=2;
498  while(count <=U->nh){
499 
500  while(U->co[count]==U->co[count-1] && count <=U->nh){
501  count++;
502  }
503 
504  bins->next=new_xyz();
505  bins=bins->next;
506  bins->x=count-count1;
507  count1=count;
508  head->x++;
509  bins->y=U->co[count-1];
510  count++;
511 
512  }
513 
514 
515  bins=head->next;
516 
517  while(bins!=NULL){
518  count=bins->x;
519  while(bins->x < mn && bins->next!=NULL){
520  bins->y=(bins->y*bins->x+bins->next->y*bins->next->x)/(bins->x+bins->next->x);
521  bins->x+=bins->next->x;
522  bins->z=0;
523  delete_xyz(head,bins);
524  (head->x)--;
525  }
526 
527  if(bins->x< mn) {
528 
529  if(bins!=head){
530  pre->x+=bins->x;
531  delete_xyz(head,pre);
532  }else {
533  printf("\nWarning::Obtaining just one bin\n");
534  }
535 
536  }
537 
538  pre=bins;
539  bins=bins->next;
540  }
541 
542  /* print_xyz_elements(head,10); */
543 
544 
545  if(N<1){
546 
547  indx=new_doublematrix(head->x,3);
548  xyz_into_doublematrix(head,indx);
549  /* print_doublematrix_elements(indx,10); */
550  delete_xyz_list(head);
551 
552  }else{
553 
554  if(head->x <N) printf("\nWarning::Some bin is remaining empty\n");
555  /* cnt=count_xyz_elements(head); */
556  indx=new_doublematrix(N,3);
557  logdelta=(log( U->co[U->nh])-log(head->next->y))/(logbase*(N-1));
558  max=log(U->co[U->nh])/log(base)+0.5*logdelta;
559  min=log(head->next->y)/logbase-0.5*logdelta;
560  head->y=head->z=min;
561  bins=head->next;
562 
563  for(j=1;j<=N;j++){
564  indx->co[j][1]=0;
565  indx->co[j][2]=0;
566  indx->co[j][3]=pow(base,min+j*logdelta);
567  tmp=0;
568  if(bins!=NULL){
569  while(log(bins->y)/logbase<=min+j*logdelta){
570  tmp+=bins->x;
571  indx->co[j][2]=(indx->co[j][1]*indx->co[j][2]+bins->y*bins->x)/tmp;
572  indx->co[j][1]=tmp;
573  bins=bins->next;
574  if(bins==NULL) break;
575  }
576  }
577  }
578 
579  delete_xyz_list(head);
580 
581  }
582 
583  /* print_doublematrix_elements(indx,10); */
584 
585  return indx;
586 
587 }
588 
589 
590 /*----------------------------------------------------------------------*/
591 void realpair_into_doublematrix(REALPAIR * head,DOUBLEMATRIX *indx )
592 {
593 
594  long i;
595  REALPAIR *tmp;
596 
597  if(head==NULL || indx==NULL){
598 
599  t_error("Inputs are not properly allocated");
600 
601  } else if(head->x!=indx->nrh) {
602 
603  t_error("Inputs do not have the same dimensions");
604  }
605 
606  tmp=head->next;
607 
608  for(i=1;i<=head->x;i++){
609  if(tmp==NULL) t_error("Something was wrong here");
610  indx->co[i][1]=tmp->x;
611  indx->co[i][2]=tmp->y;
612  tmp=tmp->next;
613  }
614 
615 }
616 
617 /*----------------------------------------------------------------------*/
618 
619 
620 void xyz_into_doublematrix(XYZ * head,DOUBLEMATRIX *indx )
621 {
622 
623  long i;
624  XYZ *tmp;
625 
626  if(head==NULL || indx==NULL){
627 
628  t_error("Inputs are not properly allocated");
629 
630  } else if(head->x!=indx->nrh) {
631 
632  t_error("Inputs do not have the same dimensions");
633  }
634 
635  tmp=head->next;
636 
637  for(i=1;i<=head->x;i++){
638  if(tmp==NULL) t_error("Something was wrong here");
639  indx->co[i][1]=tmp->x;
640  indx->co[i][2]=tmp->y;
641  indx->co[i][3]=tmp->z;
642  tmp=tmp->next;
643  }
644 
645 }
646 
647 /*---------------------------------------------------------------------------*/
648 void initialize_longvector(LONGVECTOR *L,long sign)
649 
650 {
651 
652  long i;
653 
654  if(L!=NULL){
655  if(L->isdynamic==1){
656  for(i=1;i<=L->nh;i++){
657  L->co[i]=sign;
658  }
659  }else{
660  t_error("This vector was no properly allocated");
661  }
662  }else {
663  t_error("A null vector was addressed");
664  }
665 }
666 
667 /*---------------------------------------------------------------------------*/
668 void initialize_shortvector(SHORTVECTOR *L,short sign)
669 
670 {
671 
672  long i;
673 
674  if(L!=NULL){
675  if(L->isdynamic==1){
676  for(i=1;i<=L->nh;i++){
677  L->co[i]=sign;
678  }
679  }else{
680  t_error("This vector was no properly allocated");
681  }
682  }else {
683  t_error("A null vector was addressed");
684  }
685 }
686 
687 /*---------------------------------------------------------------------------*/
688 void initialize_doublevector(DOUBLEVECTOR *L,double sign)
689 
690 {
691 
692  long i;
693 
694  if(L!=NULL){
695  if(L->isdynamic==1){
696  for(i=1;i<=L->nh;i++){
697  L->co[i]=sign;
698  }
699  }else{
700  t_error("This vector was no properly allocated");
701  }
702  }else {
703  t_error("A null vector was addressed");
704  }
705 }
706 
707 /*---------------------------------------------------------------------------*/
708 void initialize_floatvector(FLOATVECTOR *L,float sign)
709 
710 {
711 
712  long i;
713 
714  if(L!=NULL){
715  if(L->isdynamic==1){
716  for(i=1;i<=L->nh;i++){
717  L->co[i]=sign;
718  }
719  }else{
720  t_error("This vector was no properly allocated");
721  }
722  }else {
723  t_error("A null vector was addressed");
724  }
725 }
726 
727 /*---------------------------------------------------------------------------*/
728 void initialize_shortmatrix(SHORTMATRIX *L,short sign)
729 
730 {
731 
732  long i,j;
733 
734  if(L!=NULL){
735  if(L->isdynamic==1){
736  for(i=1;i<=L->nrh;i++){
737  for(j=1;j<=L->nch;j++){
738  L->co[i][j]=sign;
739  }
740  }
741  }else{
742  t_error("This matrix was no properly allocated");
743  }
744  }else {
745  t_error("A null matrix was addressed");
746  }
747 }
748 
749 
750 /*---------------------------------------------------------------------------*/
751 void initialize_longmatrix(LONGMATRIX *L,long sign)
752 
753 {
754 
755  long i,j;
756 
757  if(L!=NULL){
758  if(L->isdynamic==1){
759  for(i=1;i<=L->nrh;i++){
760  for(j=1;j<=L->nch;j++){
761  L->co[i][j]=sign;
762  }
763  }
764  }else{
765  t_error("This matrix was no properly allocated");
766  }
767  }else {
768  t_error("A null matrix was addressed");
769  }
770 }
771 
772 
773 /*---------------------------------------------------------------------------*/
774 void initialize_floatmatrix(FLOATMATRIX *L,float sign)
775 
776 {
777 
778  long i,j;
779 
780  if(L!=NULL){
781  if(L->isdynamic==1){
782  for(i=1;i<=L->nrh;i++){
783  for(j=1;j<=L->nch;j++){
784  L->co[i][j]=sign;
785  }
786  }
787  }else{
788  t_error("This matrix was no properly allocated");
789  }
790  }else {
791  t_error("A null matrix was addressed");
792  }
793 }
794 
795 /*---------------------------------------------------------------------------*/
796 void initialize_doublematrix(DOUBLEMATRIX *L,double sign)
797 
798 {
799 
800  long i,j;
801 
802  if(L!=NULL){
803  if(L->isdynamic==1){
804  for(i=1;i<=L->nrh;i++){
805  for(j=1;j<=L->nch;j++){
806  L->co[i][j]=sign;
807  }
808  }
809  }else{
810  t_error("This matrix was no properly allocated");
811  }
812  }else {
813  t_error("A null matrix was addressed");
814  }
815 }
816 
817 
818 /*--------------------------------------------------------------------------*/
819 
820 void copy_shortmatrix(SHORTMATRIX *origin,SHORTMATRIX *destination)
821 
822 {
823 
824  long i,j;
825 
826  if(origin==NULL || destination==NULL || origin->co==NULL || destination->co==NULL){
827 
828  t_error("A matrix was not allocated");
829 
830  } else if(origin->isdynamic!=1 || destination->isdynamic!=1 || origin->nrh <1 || destination->nrh <1 || origin->nch <1 || destination->nch <1 ){
831 
832  t_error("A matrix was not allocated properly");
833 
834  }else if( origin->nrh != destination->nrh || origin->nch != destination->nch ){
835 
836  t_error("The matrixes do not have the same dimensions");
837 
838  }
839 
840  for(i=1;i<=origin->nrh;i++){
841  for(j=1;j<=origin->nch;j++){
842 
843  destination->co[i][j]=origin->co[i][j];
844 
845  }
846  }
847 
848 }
849 
850 /*--------------------------------------------------------------------------*/
851 
852 void copy_intmatrix(INTMATRIX *origin,INTMATRIX *destination)
853 
854 {
855 
856  long i,j;
857 
858  if(origin==NULL || destination==NULL || origin->co==NULL || destination->co==NULL){
859 
860  t_error("A matrix was not allocated");
861 
862  } else if(origin->isdynamic!=1 || destination->isdynamic!=1 || origin->nrh <1 || destination->nrh <1 || origin->nch <1 || destination->nch <1 ){
863 
864  t_error("A matrix was not allocated properly");
865 
866  }else if( origin->nrh != destination->nrh || origin->nch != destination->nch ){
867 
868  t_error("The matrixes do not have the same dimensions");
869 
870  }
871 
872  for(i=1;i<=origin->nrh;i++){
873  for(j=1;j<=origin->nch;j++){
874 
875  destination->co[i][j]=origin->co[i][j];
876 
877  }
878  }
879 
880 }
881 
882 
883 
884 
885 
886 /*--------------------------------------------------------------------------*/
887 
888 void copy_longmatrix(LONGMATRIX *origin,LONGMATRIX *destination)
889 
890 {
891 
892  long i,j;
893 
894  if(origin==NULL || destination==NULL || origin->co==NULL || destination->co==NULL){
895 
896  t_error("A matrix was not allocated");
897 
898  } else if(origin->isdynamic!=1 || destination->isdynamic!=1 || origin->nrh <1 || destination->nrh <1 || origin->nch <1 || destination->nch <1 ){
899 
900  t_error("A matrix was not allocated properly");
901 
902  }else if( origin->nrh != destination->nrh || origin->nch != destination->nch ){
903 
904  t_error("The matrixes do not have the same dimensions");
905 
906  }
907 
908  for(i=1;i<=origin->nrh;i++){
909  for(j=1;j<=origin->nch;j++){
910 
911  destination->co[i][j]=origin->co[i][j];
912 
913  }
914  }
915 
916 }
917 
918 
919 
920 /*--------------------------------------------------------------------------*/
921 
922 void copy_floatmatrix(FLOATMATRIX *origin,FLOATMATRIX *destination)
923 
924 {
925 
926  long i,j;
927 
928  if(origin==NULL || destination==NULL || origin->co==NULL || destination->co==NULL){
929 
930  t_error("A matrix was not allocated");
931 
932  } else if(origin->isdynamic!=1 || destination->isdynamic!=1 || origin->nrh <1 || destination->nrh <1 || origin->nch <1 || destination->nch <1 ){
933 
934  t_error("A matrix was not allocated properly");
935 
936  }else if( origin->nrh != destination->nrh || origin->nch != destination->nch ){
937 
938  t_error("The matrixes do not have the same dimensions");
939 
940  }
941 
942  for(i=1;i<=origin->nrh;i++){
943  for(j=1;j<=origin->nch;j++){
944 
945  destination->co[i][j]=origin->co[i][j];
946 
947  }
948  }
949 
950 }
951 
952 
953 /*--------------------------------------------------------------------------*/
954 
955 void copy_doublematrix(DOUBLEMATRIX *origin,DOUBLEMATRIX *destination)
956 
957 {
958 
959  long i,j;
960 
961  if(origin==NULL || destination==NULL || origin->co==NULL || destination->co==NULL){
962 
963  t_error("A matrix was not allocated");
964 
965  } else if(origin->isdynamic!=1 || destination->isdynamic!=1 || origin->nrh <1 || destination->nrh <1 || origin->nch <1 || destination->nch <1 ){
966 
967  t_error("A matrix was not allocated properly");
968 
969  }else if( origin->nrh != destination->nrh || origin->nch != destination->nch ){
970 
971  t_error("The matrixes do not have the same dimensions");
972 
973  }
974 
975  for(i=1;i<=origin->nrh;i++){
976  for(j=1;j<=origin->nch;j++){
977 
978  destination->co[i][j]=origin->co[i][j];
979 
980  }
981  }
982 
983 }
984 
985 
986 
987 
988 /*--------------------------------------------------------------------------*/
989 
990 void copy_shortvector(SHORTVECTOR *origin,SHORTVECTOR *destination)
991 
992 {
993 
994  long i;
995 
996  if(origin==NULL || destination==NULL || origin->co==NULL || destination->co==NULL){
997 
998  t_error("A vector was not allocated");
999 
1000  } else if(origin->isdynamic!=1 || destination->isdynamic!=1 || origin->nh <1 || destination->nh <1 ){
1001 
1002  t_error("A vector was not allocated properly");
1003 
1004  }else if( origin->nh != destination->nh ){
1005 
1006  t_error("The vector do not have the same dimensions");
1007 
1008  }
1009 
1010  for(i=1;i<=origin->nh;i++){
1011 
1012 
1013  destination->co[i]=origin->co[i];
1014 
1015 
1016  }
1017 
1018 }
1019 
1020 /*--------------------------------------------------------------------------*/
1021 
1022 void copy_intvector(INTVECTOR *origin,INTVECTOR *destination)
1023 
1024 {
1025 
1026  long i;
1027 
1028  if(origin==NULL || destination==NULL || origin->co==NULL || destination->co==NULL){
1029 
1030  t_error("A vector was not allocated");
1031 
1032  } else if(origin->isdynamic!=1 || destination->isdynamic!=1 || origin->nh <1 || destination->nh <1 ){
1033 
1034  t_error("A vector was not allocated properly");
1035 
1036  }else if( origin->nh != destination->nh ){
1037 
1038  t_error("The vector do not have the same dimensions");
1039 
1040  }
1041 
1042  for(i=1;i<=origin->nh;i++){
1043 
1044 
1045  destination->co[i]=origin->co[i];
1046 
1047 
1048  }
1049 
1050 }
1051 
1052 
1053 /*--------------------------------------------------------------------------*/
1054 
1055 void copy_longvector(LONGVECTOR *origin,LONGVECTOR *destination)
1056 
1057 {
1058 
1059  long i;
1060 
1061  if(origin==NULL || destination==NULL || origin->co==NULL || destination->co==NULL){
1062 
1063  t_error("A vector was not allocated");
1064 
1065  } else if(origin->isdynamic!=1 || destination->isdynamic!=1 || origin->nh <1 || destination->nh <1 ){
1066 
1067  t_error("A vector was not allocated properly");
1068 
1069  }else if( origin->nh != destination->nh ){
1070 
1071  t_error("The vector do not have the same dimensions");
1072 
1073  }
1074 
1075  for(i=1;i<=origin->nh;i++){
1076 
1077 
1078  destination->co[i]=origin->co[i];
1079 
1080 
1081  }
1082 
1083 }
1084 
1085 
1086 /*--------------------------------------------------------------------------*/
1087 
1088 void copy_floatvector(FLOATVECTOR *origin,FLOATVECTOR *destination)
1089 
1090 {
1091 
1092  long i;
1093 
1094  if(origin==NULL || destination==NULL || origin->co==NULL || destination->co==NULL){
1095 
1096  t_error("A vector was not allocated");
1097 
1098  } else if(origin->isdynamic!=1 || destination->isdynamic!=1 || origin->nh <1 || destination->nh <1 ){
1099 
1100  t_error("A vector was not allocated properly");
1101 
1102  }else if( origin->nh != destination->nh ){
1103 
1104  t_error("The vector do not have the same dimensions");
1105 
1106  }
1107 
1108  for(i=1;i<=origin->nh;i++){
1109 
1110 
1111  destination->co[i]=origin->co[i];
1112 
1113 
1114  }
1115 
1116 }
1117 
1118 /*--------------------------------------------------------------------------*/
1119 
1120 void copy_doublevector(DOUBLEVECTOR *origin,DOUBLEVECTOR *destination)
1121 
1122 {
1123 
1124  long i;
1125 
1126  if(origin==NULL || destination==NULL || origin->co==NULL || destination->co==NULL){
1127 
1128  t_error("A vector was not allocated");
1129 
1130  } else if(origin->isdynamic!=1 || destination->isdynamic!=1 || origin->nh <1 || destination->nh <1 ){
1131 
1132  t_error("A vector was not allocated properly");
1133 
1134  }else if( origin->nh != destination->nh ){
1135 
1136  t_error("The vector do not have the same dimensions");
1137 
1138  }
1139 
1140  for(i=1;i<=origin->nh;i++){
1141 
1142 
1143  destination->co[i]=origin->co[i];
1144 
1145 
1146  }
1147 
1148 }
1149 
1150 
1151 /*--------------------------------------------------------------------------*/
1152 void floatmatrix_element_multiplication(FLOATMATRIX* dist,FLOATMATRIX *ca)
1153 
1154 {
1155 
1156  long i,j;
1157 
1158  if(dist==NULL || ca==NULL || dist->co==NULL || ca->co==NULL){
1159  t_error("Matrixes were not allocated");
1160  } else if(dist->isdynamic!=1 || dist->isdynamic!=1 || dist->nrh <2 ){
1161  t_error("Matrixes were not properly allocated");
1162  } else if(dist->nrh !=ca->nrh || dist->nch!=ca->nch) {
1163  t_error("These matrixes do not have the same dimensions");
1164  }
1165 
1166  for(i=1;i<=dist->nrh;i++){
1167  for(j=1;j<=dist->nch;j++){
1168  dist->co[i][j]*=ca->co[i][j];
1169  }
1170  }
1171 
1172 
1173 }
1174 
1175 /*--------------------------------------------------------------------------*/
1176 void doublematrix_element_multiplication(DOUBLEMATRIX* dist,DOUBLEMATRIX *ca)
1177 
1178 {
1179 
1180  long i,j;
1181 
1182  if(dist==NULL || ca==NULL || dist->co==NULL || ca->co==NULL){
1183  t_error("Matrixes were not allocated");
1184  } else if(dist->isdynamic!=1 || dist->isdynamic!=1 || dist->nrh <2 ){
1185  t_error("Matrixes were not properly allocated");
1186  } else if(dist->nrh !=ca->nrh || dist->nch!=ca->nch) {
1187  t_error("These matrixes do not have the same dimensions");
1188  }
1189 
1190  for(i=1;i<=dist->nrh;i++){
1191  for(j=1;j<=dist->nch;j++){
1192  dist->co[i][j]*=ca->co[i][j];
1193  }
1194  }
1195 
1196 
1197 }
1198 
1199 /*--------------------------------------------------------------------------*/
1200 void longmatrix_element_multiplication(LONGMATRIX* dist,LONGMATRIX *ca)
1201 
1202 {
1203 
1204  long i,j;
1205 
1206  if(dist==NULL || ca==NULL || dist->co==NULL || ca->co==NULL){
1207  t_error("Matrixes were not allocated");
1208  } else if(dist->isdynamic!=1 || dist->isdynamic!=1 || dist->nrh <2 ){
1209  t_error("Matrixes were not properly allocated");
1210  } else if(dist->nrh !=ca->nrh || dist->nch!=ca->nch) {
1211  t_error("These matrixes do not have the same dimensions");
1212  }
1213 
1214  for(i=1;i<=dist->nrh;i++){
1215  for(j=1;j<=dist->nch;j++){
1216  dist->co[i][j]*=ca->co[i][j];
1217  }
1218  }
1219 
1220 
1221 }
1222 
1223 
1224 /*--------------------------------------------------------------------------*/
1225 void shortmatrix_element_multiplication(SHORTMATRIX* dist,SHORTMATRIX *ca)
1226 
1227 {
1228 
1229  long i,j;
1230 
1231  if(dist==NULL || ca==NULL || dist->co==NULL || ca->co==NULL){
1232  t_error("Matrixes were not allocated");
1233  } else if(dist->isdynamic!=1 || dist->isdynamic!=1 || dist->nrh <2 ){
1234  t_error("Matrixes were not properly allocated");
1235  } else if(dist->nrh !=ca->nrh || dist->nch!=ca->nch) {
1236  t_error("These matrixes do not have the same dimensions");
1237  }
1238 
1239  for(i=1;i<=dist->nrh;i++){
1240  for(j=1;j<=dist->nch;j++){
1241  dist->co[i][j]*=ca->co[i][j];
1242  }
1243  }
1244 
1245 
1246 }
1247 
1248 /*--------------------------------------------------------------------------*/
1249 void shortvector_element_multiplication(SHORTVECTOR* dist,SHORTVECTOR *ca)
1250 
1251 {
1252 
1253  long i;
1254 
1255  if(dist==NULL || ca==NULL || dist->co==NULL || ca->co==NULL){
1256  t_error("Vectors were not allocated");
1257  } else if(dist->isdynamic!=1 || dist->isdynamic!=1 || dist->nh <2 ){
1258  t_error("Vectors were not properly allocated");
1259  } else if(dist->nh !=ca->nh ) {
1260  t_error("These vectors do not have the same dimensions");
1261  }
1262 
1263  for(i=1;i<=dist->nh;i++){
1264  dist->co[i]*=ca->co[i];
1265  }
1266 
1267 
1268 }
1269 
1270 /*--------------------------------------------------------------------------*/
1271 void longvector_element_multiplication(LONGVECTOR* dist,LONGVECTOR *ca)
1272 
1273 {
1274 
1275  long i;
1276 
1277  if(dist==NULL || ca==NULL || dist->co==NULL || ca->co==NULL){
1278  t_error("Vectors were not allocated");
1279  } else if(dist->isdynamic!=1 || dist->isdynamic!=1 || dist->nh <2 ){
1280  t_error("Vectors were not properly allocated");
1281  } else if(dist->nh !=ca->nh ) {
1282  t_error("These vectors do not have the same dimensions");
1283  }
1284 
1285  for(i=1;i<=dist->nh;i++){
1286  dist->co[i]*=ca->co[i];
1287  }
1288 
1289 
1290 }
1291 
1292 /*--------------------------------------------------------------------------*/
1293 void floatvector_element_multiplication(FLOATVECTOR* dist,FLOATVECTOR *ca)
1294 
1295 {
1296 
1297  long i;
1298 
1299  if(dist==NULL || ca==NULL || dist->co==NULL || ca->co==NULL){
1300  t_error("Vectors were not allocated");
1301  } else if(dist->isdynamic!=1 || dist->isdynamic!=1 || dist->nh <2 ){
1302  t_error("Vectors were not properly allocated");
1303  } else if(dist->nh !=ca->nh ) {
1304  t_error("These vectors do not have the same dimensions");
1305  }
1306 
1307  for(i=1;i<=dist->nh;i++){
1308  dist->co[i]*=ca->co[i];
1309  }
1310 
1311 
1312 }
1313 
1314 /*--------------------------------------------------------------------------*/
1315 void doublevector_element_multiplication(DOUBLEVECTOR* dist,DOUBLEVECTOR *ca)
1316 
1317 {
1318 
1319  long i;
1320 
1321  if(dist==NULL || ca==NULL || dist->co==NULL || ca->co==NULL){
1322  t_error("Vectors were not allocated");
1323  } else if(dist->isdynamic!=1 || dist->isdynamic!=1 || dist->nh <2 ){
1324  t_error("Vectors were not properly allocated");
1325  } else if(dist->nh !=ca->nh ) {
1326  t_error("These vectors do not have the same dimensions");
1327  }
1328 
1329  for(i=1;i<=dist->nh;i++){
1330  dist->co[i]*=ca->co[i];
1331  }
1332 
1333 
1334 }
1335 /*-----------------------------------------------------------------------*/
1336 
1337 
1338 void clean_floatmatrix(FLOATMATRIX *iv,FLOATMATRIX *ov,FLOATVECTOR *U,FLOATVECTOR *V)
1339 
1340 
1341 {
1342 
1343  long i,j;
1344 
1345  if(ov->nrh!=iv->nrh || ov->nch!=iv->nch){
1346 
1347  t_error("Matrixes do not have proper dimensions");
1348 
1349  }
1350 
1351  if(V->co[1] >0){
1352  for(i=1;i<=iv->nrh;i++){
1353  for(j=1;j<=iv->nch;j++){
1354  if(ov->co[i][j] >= V->co[2]){
1355  iv->co[i][j]=U->co[2];
1356  }
1357 
1358 
1359  }
1360  }
1361  }else if(V->co[1] <0){
1362 
1363  for(i=1;i<=iv->nrh;i++){
1364  for(j=1;j<=iv->nch;j++){
1365  if(ov->co[i][j] <= V->co[2]){
1366  iv->co[i][j]=U->co[2];
1367  }
1368 
1369  }
1370  }
1371 
1372  }else{
1373 
1374  for(i=1;i<=iv->nrh;i++){
1375  for(j=1;j<=iv->nch;j++){
1376  if(ov->co[i][j] == V->co[2]){
1377  iv->co[i][j]=U->co[2];
1378  }
1379 
1380  }
1381  }
1382 
1383 
1384  }
1385 
1386 }
1387 
1388 /*---------------------------------------------------------------------------*/
1389 
1390 DOUBLEBIN *split(DOUBLEVECTOR *U,long N,FLOATVECTOR *novalue)
1391 
1392 
1393 {
1394 
1395  double delta,min,max;
1396  long i,count,count1,minposition,maxposition,bin_vuoti,num_max;
1397  DOUBLEBIN *l;
1398  LONGPAIR *head,*bins;
1399  /* char ch; */
1400 
1401  minposition=1;
1402  maxposition=U->nh;
1403 
1404  if(N <=1){
1405 
1406 
1407  head=new_longpair();
1408  head->i=0;
1409  bins=head;
1410  count1=1;
1411  count=2;
1412  printf("This options is memory expensive\n");
1413  printf("ENTER THE MAXIMUM NUMBER OF ALLOWED BIN\n");
1414  scanf("%ld",&num_max);
1415 
1416  while(count <= U->nh){
1417  while(U->co[count]==U->co[count-1] && count <=U->nh){
1418  count++;
1419  }
1420  /* printf("%f -> %d-%d=%d\n",U->co[count-1],count,count1,bins->i);
1421  scanf("%c",&ch);
1422  */
1423 
1424  bins->next=new_longpair();
1425  bins=bins->next;
1426  bins->i=count-count1;
1427  bins->j=1;
1428  count1=count-1;
1429  head->i++;
1430  count1=count;
1431  count++;
1432  if(head->i>num_max)t_error("Il numero di bin eccede il numero massimo che hai consentito");
1433 
1434  }
1435 
1436  }else if(N >1){
1437  if(novalue->co[1]==0){
1438  minposition=1;
1439  maxposition=U->nh;
1440  }else if(novalue->co[1]==-1){
1441  minposition=2;
1442  /* printf("*****%d\n",novalue); */
1443  max=U->co[U->nh];
1444  while(U->co[minposition]<=novalue->co[2]){
1445  minposition++;
1446  }
1447  min=U->co[minposition];
1448  /*printf("***%f %f\n",min,U->co[minposition-1]);
1449  scanf(" %c",&ch);*/
1450  maxposition=U->nh;
1451  }else{
1452  maxposition=U->nh-1;
1453  while(U->co[maxposition]>=novalue->co[2]){
1454  /* printf("%d\n",maxposition);
1455  scanf("%c",&ch);
1456  */ maxposition--;
1457  }
1458  max=U->co[maxposition];
1459  minposition=1;
1460  min=U->co[1];
1461  }
1462 
1463  printf("THE MINIMUM VALUE IS %f\n",min);
1464  printf("THE MAXIMUM VALUE IS %f\n",max);
1465  delta=(max-min)/(N-1);
1466  printf("DELTA IS%f\n",delta);
1467  /* scanf("%c",&ch); */
1468  head=new_longpair();
1469  head->i=0;
1470  bins=head;
1471  count1=minposition;
1472  count=minposition+1;
1473  bin_vuoti=0;
1474 
1475  while(count <= maxposition){
1476  if(U->co[count] < min+0.5*delta){
1477  while(U->co[count] < min+0.5*delta && count<=maxposition){
1478  count++;
1479  }
1480 
1481  bins->next=new_longpair();
1482  bins=bins->next;
1483  bins->i=count-count1;
1484  /* printf("%f -> %d-%d=%d\n",U->co[count-1],count,count1,bins->i);
1485  scanf("%c",&ch);
1486  */
1487  bins->j=1;
1488  count1=count-1;
1489  (head->i)++;
1490  count1=count;
1491  count++;
1492 
1493  }else{
1494  bin_vuoti++;
1495  printf("\nWarning:: An empty bin was encountered\n");
1496  }
1497  min+=delta;
1498  }
1499  if(bin_vuoti!=0){
1500  printf("\nThere was look for (%ld) empty bins\n",bin_vuoti);
1501  }
1502  }
1503 
1504  l=(DOUBLEBIN *)malloc(sizeof(DOUBLEBIN));
1505  if (!l) t_error("allocation failure in new_doublebin()");
1506  l->isdynamic=isDynamic;
1507  if(head->i < 2 ){
1508  t_error("Something wrong happened in binning");
1509  }else{
1510  l->index=new_longvector(head->i);
1511  (l->index)->nl=1;
1512  (l->index)->nh=head->i;
1513  bins=head->next;
1514  for(i=1;i<=head->i;i++){
1515  (l->index)->co[i]=bins->i;
1516  bins=bins->next;
1517  if(bins==NULL) break;
1518  }
1519  l->co=(double **)malloc((size_t) ((l->index->nh+NR_END)*sizeof(double *)));
1520  if (!l->co) t_error("allocation failure in new_doublebin()");
1521  l+=NR_END-1;
1522 
1523  l->co[1]=U->co+minposition-1;
1524  if(!(l->co[1])) t_error("Failure in splitting this vector");
1525  l->co[1]+=NR_END;
1526  l->co[1]-=1;
1527  bins=head->next;
1528  for(i=2;i<=l->index->nh; i++){
1529  l->co[i]=l->co[i-1]+bins->i;
1530  bins=bins->next;
1531  if(bins==NULL) break;
1532  }
1533 
1534  }
1535 
1536 
1537  U->isdynamic=-1;
1538  U->nl=-1;
1539  U->nh=1;
1540 
1541  free(U);
1542 
1543  delete_longpair_list(head);
1544 
1545  return l;
1546 
1547 
1548 }
1549 
1550 
1551 /*---------------------------------------------------------------------------*/
1552 
1553 DOUBLEBIN *esponentialsplit(DOUBLEVECTOR *U,long N,double base,FLOATVECTOR* novalue)
1554 
1555 
1556 {
1557 
1558  double min,max,logdelta,logbase,logmin;
1559  long i,count,count1,minposition,maxposition,bin_vuoti,num_max;
1560  DOUBLEBIN *l;
1561  LONGPAIR *head,*bins;
1562  /*char ch;
1563 
1564  minposition=1;
1565  maxposition=U->nh;
1566  printf("%d\n",N);
1567  //scanf("%c",&ch);*/
1568  logbase=log10(base);
1569 
1570  if(N <=1){
1571  printf("This options is memory expensive\n");
1572  printf("ENTER THE MAXIMUM NUMBER OF ALLOWED BIN\n");
1573  scanf("%ld",&num_max);
1574  head=new_longpair();
1575  head->i=0;
1576  bins=head;
1577  count1=1;
1578  count=2;
1579 
1580  while(count <= U->nh){
1581  while(U->co[count]==U->co[count-1] && count <=U->nh){
1582  count++;
1583  }
1584 
1585  bins->next=new_longpair();
1586  bins=bins->next;
1587  bins->i=count-count1;
1588  /* printf("%f -> %d-%d=%d\n",U->co[count-1],count,count1,bins->i);
1589  scanf("%c",&ch);
1590  */
1591  bins->j=1;
1592  count1=count-1;
1593  head->i++;
1594  count1=count;
1595  count++;
1596  if(head->i>num_max)t_error("Il numero di bin eccede il numero massimo che hai consentito");
1597  }
1598  }else if(N >1){
1599  if(novalue->co[1]==0){
1600  minposition=1;
1601  if(U->co[minposition] <=0){
1602  while(U->co[minposition]<=0){
1603  minposition++;
1604  if( minposition > U->nh) t_error("No valid data for esponential binning");
1605  }
1606  }
1607  maxposition=U->nh;
1608  }else if(novalue->co[1]==-1){
1609  minposition=2;
1610  /* printf("*****%d\n",novalue); */
1611  max=U->co[U->nh];
1612  while(U->co[minposition]<=novalue->co[2]){
1613  minposition++;
1614  }
1615 
1616  min=U->co[minposition];
1617  maxposition=U->nh;
1618  }else{
1619  maxposition=U->nh-1;
1620  while(U->co[maxposition]>=novalue->co[2]){
1621  /* printf("%d\n",maxposition);
1622  scanf("%c",&ch);
1623  */ maxposition--;
1624  }
1625  max=U->co[maxposition];
1626  minposition=1;
1627  min=U->co[1];
1628  }
1629 
1630  printf("THE MINIMUM VALUE IS %f\n",min);
1631  printf("THE MAXIMUM VALUE IS %f\n",max);
1632  logdelta=(log10(max)-log10(min))/(logbase*(N-1));
1633  printf("LOG_10 DELTA IS%f\n",logdelta);
1634  /* scanf("%c",&ch); */
1635  head=new_longpair();
1636  head->i=0;
1637  bins=head;
1638  count1=minposition;
1639  count=minposition+1;
1640  logmin=log10(min);
1641  bin_vuoti=0;
1642 
1643  while(count <= maxposition){
1644  if(log10(U->co[count]) < logmin+0.5*logdelta){
1645  while(log10(U->co[count]) < logmin+0.5*logdelta && count<=maxposition){
1646  count++;
1647  }
1648 
1649  bins->next=new_longpair();
1650  bins=bins->next;
1651  bins->i=count-count1;
1652  /* printf("%f -> %d-%d=%d\n",U->co[count-1],count,count1,bins->i);
1653  scanf("%c",&ch);
1654  */
1655  bins->j=1;
1656  count1=count-1;
1657  (head->i)++;
1658  count1=count;
1659  count++;
1660  }else{
1661  bin_vuoti++;
1662  printf("\nWarning:: An empty bin was encountered\n");
1663  }
1664  logmin+=logdelta;
1665  }
1666  if(bin_vuoti!=0){
1667  printf("\nThere was look for (%ld) empty bins\n",bin_vuoti);
1668  }
1669  }
1670 
1671  l=(DOUBLEBIN *)malloc(sizeof(DOUBLEBIN));
1672  if (!l) t_error("allocation failure in new_doublebin()");
1673  l->isdynamic=isDynamic;
1674  if(head->i < 2 ){
1675  t_error("Something wrong happened in binning");
1676  }else{
1677  l->index=new_longvector(head->i);
1678  (l->index)->nl=1;
1679  (l->index)->nh=head->i;
1680  bins=head->next;
1681  for(i=1;i<=head->i;i++){
1682  (l->index)->co[i]=bins->i;
1683  bins=bins->next;
1684  if(bins==NULL) break;
1685  }
1686  l->co=(double **)malloc((size_t) ((l->index->nh+NR_END)*sizeof(double *)));
1687  if (!l->co) t_error("allocation failure in new_doublebin()");
1688  l+=NR_END-1;
1689 
1690  l->co[1]=U->co+minposition-1;
1691  if(!(l->co[1])) t_error("Failure in splitting this vector");
1692  l->co[1]+=NR_END;
1693  l->co[1]-=1;
1694  bins=head->next;
1695  for(i=2;i<=l->index->nh; i++){
1696  l->co[i]=l->co[i-1]+bins->i;
1697  bins=bins->next;
1698  if(bins==NULL) break;
1699  }
1700 
1701  }
1702 
1703 
1704  U->isdynamic=-1;
1705  U->nl=-1;
1706  U->nh=1;
1707 
1708  free(U);
1709 
1710  delete_longpair_list(head);
1711 
1712  return l;
1713 
1714 
1715 }
1716 
1717 /*---------------------------------------------------------------------------*/
1718 
1719 double split2realvectors(DOUBLEVECTOR *U,DOUBLEVECTOR *T,DOUBLEBIN *l,DOUBLEBIN *m, long N,long num_max,FLOATVECTOR *novalue)
1720 
1721 
1722 {
1723 
1724  double delta,min,max;
1725  long i,count,count1,minposition,maxposition,bin_vuoti;
1726  LONGPAIR *head,*bins;
1727 
1728 
1729  minposition=1;
1730  maxposition=U->nh;
1731  /*printf("%d\n",N);*/
1732  /*scanf("%c",&ch);*/
1733  if(N <=1){
1734  head=new_longpair();
1735  head->i=0;
1736  bins=head;
1737  count1=1;
1738  count=2;
1739  while(count <= U->nh){
1740  while(U->co[count]==U->co[count-1] && count <=U->nh){
1741  count++;
1742  }
1743 
1744  bins->next=new_longpair();
1745  bins=bins->next;
1746  bins->i=count-count1;
1747  /*printf("%f -> %d-%d=%d\n",U->co[count-1],count,count1,bins->i);
1748  scanf("%c",&ch);*/
1749 
1750  bins->j=1;
1751  count1=count-1;
1752  head->i++;
1753  count1=count;
1754  count++;
1755  if(head->i>num_max)t_error("Il numero di bin eccede il numero massimo che hai consentito");
1756 
1757  }
1758 
1759  }else if(N >1){
1760  if(novalue->co[1]==0){
1761  minposition=1;
1762  maxposition=U->nh;
1763  }else if(novalue->co[1]==-1){
1764  minposition=2;
1765  /*printf("*****%d\n",novalue);*/
1766  /*scanf("%hd",&ch);*/
1767  max=U->co[U->nh];
1768  while(U->co[minposition]<=novalue->co[2]){
1769  /*printf("%d\n",minposition);
1770  scanf("%c",&ch);*/
1771  minposition++;
1772  }
1773  min=U->co[minposition];
1774  maxposition=U->nh;
1775  }else{
1776  maxposition=U->nh-1;
1777  while(U->co[maxposition]>=novalue->co[2]){
1778  /*printf("%d\n",maxposition);
1779  scanf("%c",&ch);*/
1780  maxposition--;
1781  }
1782  max=U->co[maxposition];
1783  minposition=1;
1784  min=U->co[1];
1785  }
1786 
1787  printf("THE MINIMUM VALUE IS %f\n",min);
1788  printf("THE MAXIMUM VALUE IS %f\n",max);
1789  delta=(max-min)/(N-1);
1790  printf("DELTA IS%f\n",delta);
1791  /* scanf("%c",&ch); */
1792  head=new_longpair();
1793  head->i=0;
1794  bins=head;
1795  count1=minposition;
1796  count=minposition+1;
1797  bin_vuoti=0;
1798  while(count <= maxposition){
1799  if(U->co[count] < min +0.5*delta){
1800  while(U->co[count] < min+0.5*delta && count<=maxposition){
1801  count++;
1802  }
1803  bins->next=new_longpair();
1804  bins=bins->next;
1805  bins->i=count-count1;
1806  /*printf("%f -> %d-%d=%d\n",U->co[count-1],count,count1,bins->i);
1807  scanf("%c",&ch);*/
1808 
1809  bins->j=1;
1810  count1=count-1;
1811  (head->i)++;
1812  count1=count;
1813  count++;
1814 
1815  }else{
1816  bin_vuoti++;
1817  printf("\nWarning:: An empty bin was encountered\n");
1818  }
1819  min+=delta;
1820  }
1821  if(bin_vuoti!=0){
1822  printf("\nThere was look for (%ld) empty bins\n",bin_vuoti);
1823  }
1824 
1825  }
1826  if(head->i < 2 ){
1827  t_error("Something wrong happened in binning");
1828  }else{
1829  m->index=new_longvector(head->i);
1830  l->index=new_longvector(head->i);
1831  (l->index)->nl=1;
1832  (l->index)->nh=head->i;
1833  (m->index)->nl=1;
1834  (m->index)->nh=head->i;
1835 
1836  bins=head->next;
1837  for(i=1;i<=head->i;i++){
1838  (l->index)->co[i]=bins->i;
1839  (m->index)->co[i]=bins->i;
1840  bins=bins->next;
1841  if(bins==NULL) break;
1842  }
1843  l->co=(double **)malloc((size_t) ((l->index->nh+NR_END)*sizeof(double *)));
1844  m->co=(double **)malloc((size_t) ((m->index->nh+NR_END)*sizeof(double *)));
1845 
1846  if (!l->co || !m->co) t_error("allocation failure in new_doublebin()");
1847  l+=NR_END-1;
1848  l->co[1]=U->co+minposition-1;
1849  m+=NR_END-1;
1850  m->co[1]=T->co+minposition-1;
1851 
1852  if(!(l->co[1]) || !(m->co[1])) t_error("Failure in splitting this vector");
1853  l->co[1]+=NR_END;
1854  l->co[1]-=1;
1855  m->co[1]+=NR_END;
1856  m->co[1]-=1;
1857 
1858  bins=head->next;
1859  for(i=2;i<=l->index->nh; i++){
1860  l->co[i]=l->co[i-1]+bins->i;
1861  m->co[i]=m->co[i-1]+bins->i;
1862  bins=bins->next;
1863  if(bins==NULL) break;
1864  }
1865 
1866  }
1867  U->isdynamic=-1;
1868  U->nl=-1;
1869  U->nh=1;
1870 
1871  free(U);
1872  T->isdynamic=-1;
1873  T->nl=-1;
1874  T->nh=1;
1875 
1876  free(T);
1877 
1878  delete_longpair_list(head);
1879 
1880  if(N < 2) delta=0;
1881  return delta;
1882 
1883 }
1884 
1885 
1886 
1887 /*---------------------------------------------------------------------------*/
1888 
1889 double esponentialsplit2realvectors(DOUBLEVECTOR *U,DOUBLEVECTOR *W, DOUBLEBIN* l,DOUBLEBIN* m,long N,long num_max,double base,FLOATVECTOR *novalue)
1890 
1891 
1892 {
1893 
1894  double min,max,logdelta,logbase,logmin;
1895  long i,j,count,count1,minposition,maxposition,bin_vuoti;
1896  LONGPAIR *head,*bins;
1897  short mode;
1898  SHORTVECTOR *segna;
1899  minposition=1;
1900  maxposition=U->nh;
1901  logbase=log10(base);
1902 
1903  if(N <=1){
1904  head=new_longpair();
1905  head->i=0;
1906  bins=head;
1907  count1=1;
1908  count=2;
1909 
1910  while(count <= U->nh){
1911  while(U->co[count]==U->co[count-1] && count <=U->nh){
1912  count++;
1913  }
1914 
1915  bins->next=new_longpair();
1916  bins=bins->next;
1917  bins->i=count-count1;
1918  /*printf("%f -> %d-%d=%d\n",U->co[count-1],count,count1,bins->i);
1919  scanf("%c",&ch);*/
1920 
1921  bins->j=1;
1922  count1=count-1;
1923  head->i++;
1924  count1=count;
1925  count++;
1926  if(head->i>num_max)t_error("Il numero di bin eccede il numero massimo che hai consentito");
1927  }
1928  }else if(N >1){
1929  if(novalue->co[1]==0){
1930  minposition=1;
1931  if(U->co[minposition] <=0){
1932  while(U->co[minposition]<=0){
1933  minposition++;
1934  if( minposition > U->nh) t_error("No valid data for esponential binning");
1935  }
1936  }
1937  maxposition=U->nh;
1938  }else if(novalue->co[1]==-1){
1939  minposition=2;
1940  /* printf("*****%d\n",novalue->co[1]); */
1941  max=U->co[U->nh];
1942  while(U->co[minposition]<=novalue->co[2]){
1943  minposition++;
1944  }
1945  min=U->co[minposition];
1946  maxposition=U->nh;
1947  }else{
1948  maxposition=U->nh-1;
1949  while(U->co[maxposition]>=novalue->co[2]){
1950  /* printf("%d\n",maxposition);
1951  scanf("%c",&ch);
1952  */ maxposition--;
1953  }
1954  max=U->co[maxposition];
1955  minposition=1;
1956  min=U->co[1];
1957  }
1958  printf("THE MINIMUM VALUE IS \t%f\n",min);
1959  printf("THE MAXIMUM VALUE IS \t%f\n",max);
1960  logdelta=(log10(max)-log10(min))/(logbase*(N-1));
1961  printf("LOG_10 DELTA IS \t%f\t",logdelta);
1962  /*scanf("%c",&ch); */
1963  head=new_longpair();
1964  head->i=0;
1965  bins=head;
1966  count1=minposition;
1967  count=minposition+1;
1968  logmin=log10(min);
1969  segna=new_shortvector(N);
1970  j=0;
1971  bin_vuoti=0;
1972  while(count <= maxposition){
1973  mode=1;
1974  j++;
1975  segna->co[j]=0;
1976  if(log10(U->co[count])<=logmin+0.5*logdelta){
1977  while(log10(U->co[count]) <= logmin+0.5*logdelta && count<=maxposition){
1978  count++;
1979  }
1980  bins->next=new_longpair();
1981  bins=bins->next;
1982  bins->i=count-count1;
1983  /*printf("%d %f -> %d-%d=%d\n",j,U->co[count-1],count,count1,bins->i);
1984  bins->j=1;*/
1985  count1=count-1;
1986  (head->i)++;
1987  count1=count;
1988  count++;
1989  }else{
1990  bin_vuoti++;
1991  printf("\nWarning:: An empty bin was encountered\n");
1992  }
1993  logmin+=logdelta;
1994  }
1995  if(bin_vuoti!=0){
1996  printf("\nThere was look for (%ld) empty bins\n",bin_vuoti);
1997  }
1998  }
1999  if(head->i < 2 ){
2000  t_error("Something wrong happened in binning");
2001  }else{
2002  m->index=new_longvector(head->i);
2003  l->index=new_longvector(head->i);
2004  (l->index)->nl=1;
2005  (l->index)->nh=head->i;
2006  (m->index)->nl=1;
2007  (m->index)->nh=head->i;
2008  bins=head->next;
2009  for(i=1;i<=head->i;i++){
2010  (l->index)->co[i]=bins->i;
2011  (m->index)->co[i]=bins->i;
2012  bins=bins->next;
2013  if(bins==NULL) break;
2014  }
2015  l->co=(double **)malloc((size_t) ((l->index->nh+NR_END)*sizeof(double *)));
2016  m->co=(double **)malloc((size_t) ((m->index->nh+NR_END)*sizeof(double *)));
2017 
2018  if (!l->co || !m->co) t_error("allocation failure in new_doublebin()");
2019  l+=NR_END-1;
2020  l->co[1]=U->co+minposition-1;
2021  m+=NR_END-1;
2022  m->co[1]=W->co+minposition-1;
2023 
2024  if(!(l->co[1]) || !(m->co[1])) t_error("Failure in splitting this vector");
2025  l->co[1]+=NR_END;
2026  l->co[1]-=1;
2027  m->co[1]+=NR_END;
2028  m->co[1]-=1;
2029 
2030  bins=head->next;
2031  for(i=2;i<=l->index->nh; i++){
2032  l->co[i]=l->co[i-1]+bins->i;
2033  m->co[i]=m->co[i-1]+bins->i;
2034  bins=bins->next;
2035  if(bins==NULL) break;
2036  }
2037 
2038  }
2039 
2040  U->isdynamic=-1;
2041  U->nl=-1;
2042  U->nh=1;
2043 
2044  free(U);
2045  W->isdynamic=-1;
2046  W->nl=-1;
2047  W->nh=1;
2048 
2049  free(W);
2050 
2051  delete_longpair_list(head);
2052 
2053  return logdelta;
2054 
2055 }
2056 
2057 
2058 
2059 
2060 DOUBLEMATRIX *shrink_doublematrix(DOUBLEMATRIX *data,FLOATVECTOR *Q)
2061 
2062 {
2063 
2064  long i,count=0;
2065 
2066 
2067  DOUBLEMATRIX *newm;
2068 
2069  if(Q->co[1]<0){
2070 
2071  for(i=1;i<=data->nrh;i++){
2072 
2073  if(data->co[i][2] > Q->co[2]){
2074  count++;
2075  }
2076  }
2077  }else if(Q->co[1]>0){
2078 
2079  for(i=1;i<=data->nrh;i++){
2080 
2081  if(data->co[i][2] < Q->co[2]){
2082  count++;
2083  }
2084  }
2085 
2086  } else{
2087 
2088  for(i=1;i<=data->nrh;i++){
2089 
2090  if(data->co[i][2] != Q->co[2]){
2091  count++;
2092  }
2093  }
2094 
2095 
2096  }
2097 
2098 
2099  newm=new_doublematrix(count,2);
2100 
2101  count=1;
2102 
2103  if(Q->co[1]<0){
2104 
2105  for(i=1;i<=data->nrh;i++){
2106 
2107  if(data->co[i][2] > Q->co[2]){
2108  newm->co[count][1]=data->co[i][1];
2109  newm->co[count][2]=data->co[i][2];
2110 
2111  count++;
2112  }
2113  }
2114  }else if(Q->co[1]>0){
2115 
2116  for(i=1;i<=data->nrh;i++){
2117 
2118  if(data->co[i][2] < Q->co[2]){
2119  newm->co[count][1]=data->co[i][1];
2120  newm->co[count][2]=data->co[i][2];
2121  count++;
2122  }
2123  }
2124 
2125  } else{
2126 
2127  for(i=1;i<=data->nrh;i++){
2128 
2129  if(data->co[i][2] != Q->co[2]){
2130  newm->co[count][1]=data->co[i][1];
2131  newm->co[count][2]=data->co[i][2];
2132 
2133  count++;
2134  }
2135  }
2136 
2137 
2138  }
2139 
2140 
2141  return newm;
2142 
2143 }
2144 
2145 /* interpolating_function */
2146 
2147 DOUBLEMATRIX *interpolating_function(DOUBLEMATRIX *cleandata)
2148 
2149 {
2150  long i;
2151  DOUBLEMATRIX * W;
2152 
2153  W=new_doublematrix(cleandata->nrh,2);
2154 
2155  for(i=1;i<=cleandata->nrh-1;i++){
2156  W->co[i][1]=(cleandata->co[i+1][2]-cleandata->co[i][2])/(cleandata->co[i+1][1]-cleandata->co[i][1]);
2157  W->co[i][2]=cleandata->co[i][2]-W->co[i][1]*cleandata->co[i][1];
2158  }
2159 
2160  return W;
2161 }
2162 
2163 /* interpolate */
2164 
2165 /*
2166 
2167 double interpolate(double x,DOUBLEMATRIX *cleandata,DOUBLEMATRIX *W)
2168 
2169 {
2170 long count;
2171 double y=0;
2172 
2173 count=W->co[W->nrh][1];
2174 
2175 if(x< cleandata->co[1][1] || x>cleandata->co[cleandata->nrh][1] ){
2176 t_error("Out of bound");
2177 } else if ( x < cleandata->co[count][1]){
2178 while(x < cleandata->co[count][1]){
2179 count--;
2180 }
2181 W->co[W->nrh][1]=count;
2182 y=W->co[count][1]*x+W->co[count][2];
2183 }else{
2184 while(x > cleandata->co[count][1]){
2185 count++;
2186 }
2187 count--;
2188 W->co[W->nrh][1]=count;
2189 y=W->co[count][1]*x+W->co[count][2];
2190 }
2191 
2192 return y;
2193 }
2194 
2195 */
2196 
2197 /*--------------------------------------------------------------------------*/
2198 
2199 double interpolate(double x,DOUBLEMATRIX *cleandata,DOUBLEMATRIX *W)
2200 
2201 {
2202  long count;
2203  double y=0;
2204 
2205  count=W->nrh;
2206 
2207  if(x< cleandata->co[1][1] || x>cleandata->co[cleandata->nrh][1] ){
2208  t_error("Out of bound");
2209  } else if ( x <= cleandata->co[count][1]){
2210  while(x <= cleandata->co[count][1]){
2211  count--;
2212  if(count<=0){
2213  count++;
2214  break;
2215  }
2216  }
2217 
2218  y=W->co[count][1]*x+W->co[count][2];
2219  }
2220 
2221  return y;
2222 }
2223 
2224 /*--------------------------------------------------------------------------*/
2225 
2226 
2227 /* interpolate_floatmatrix */
2228 
2229 void interpolate_floatmatrix(FLOATMATRIX *matrice, float dt, float istante, FLOATVECTOR *vettore)
2230 
2231 
2238 {
2239  long i1,j,time;
2240  double resto;
2241 
2242  time=floor(istante/dt);
2243 
2244  if((istante/dt-time)==0){
2245  i1=time;
2246  }else{
2247  i1=time+1;
2248  }
2249 
2250  if(i1<1) t_error("tempo troppo basso");
2251  if(i1>=matrice->nrh) t_error("tempo troppo alto");
2252  resto=(istante-(i1-1)*dt)/dt;
2253  for(j=1;j<=matrice->nch;j++){
2254  vettore->co[j]=matrice->co[i1][j]+(matrice->co[i1+1][j]-matrice->co[i1][j])*resto;
2255  }
2256 }
2257 
2258 /*--------------------------------------------------------------------------*/
2259 
2260 /* interpolate_doublematrix */
2261 
2262 void interpolate_doublematrix(DOUBLEMATRIX *matrice, float dt, float istante, DOUBLEVECTOR *vettore)
2263 
2264 
2272 {
2273  long i1,j,time;
2274  double resto;
2275 
2276  time=floor(istante/dt);
2277 
2278  if((istante/dt-time)==0){
2279  i1=time;
2280  }else{
2281  i1=time+1;
2282  }
2283 
2284  if(i1<1){
2285  printf("istant=%f\n",istante);
2286  t_error("tempo troppo basso");
2287  }
2288  if(i1>=matrice->nrh){
2289  printf("istant=%f,imax=%f\n",istante,(matrice->nrh-1)*dt);
2290  t_error("tempo troppo alto");
2291  }
2292  resto=(istante-(i1-1)*dt)/dt;
2293  for(j=1;j<=matrice->nch;j++){
2294  vettore->co[j]=matrice->co[i1][j]+(matrice->co[i1+1][j]-matrice->co[i1][j])*resto;
2295  }
2296 }
2297 
2298 /*--------------------------------------------------------------------------*/
2299 double mean_doublematrix_column(DOUBLEMATRIX* net,long column)
2300 {
2301 
2302  double mean;
2303  long i;
2304 
2305  mean=0;
2306 
2307  for(i=1;i<=net->nrh;i++){
2308  mean+=net->co[i][column];
2309  }
2310 
2311  return mean/=net->nrh;
2312 
2313 }
2314 /*--------------------------------------------------------------------------*/
2315 
2316 double variance_doublematrix_column(DOUBLEMATRIX* net,long column,double mean)
2317 
2318 {
2319 
2320  double variance;
2321  long i;
2322 
2323  variance=0;
2324 
2325  for(i=1;i<=net->nrh;i++){
2326  variance+=(net->co[i][column]-mean)*(net->co[i][column]-mean);
2327  }
2328 
2329  return variance/=net->nrh;
2330 
2331 }
2332 
2333 /*--------------------------------------------------------------------------*/
2334 
2335 double approximate_2_multiple(double number,double div)
2336 
2337 {
2338 
2339  return number-fabs(fmod(number,div));
2340 
2341 }
2342 
2343 /*-------- ricampiona ------------------------------------------------------------------*/
2344 
2345 DOUBLEMATRIX *ricampiona(DOUBLEMATRIX *cleandata, float ti2, float dt2, long n2, float dt1)
2346 
2347 {
2348  long i,j,k,nint;
2349  DOUBLEMATRIX *result, *coeff, *xy, *interval, *trapezio;
2350  long count;
2351  float perc;
2352 
2353  FLOATVECTOR *unused;
2354  unused=new_floatvector(1);
2355  initialize_floatvector(unused,0);
2356 
2357 
2358 
2359  /* matrix with x e y data */
2360  xy=new_doublematrix(cleandata->nrh,2);
2361  initialize_doublematrix(xy,0);
2362 
2363  /* assign xy x data */
2364  for(i=1;i<=cleandata->nrh;i++){
2365  xy->co[i][1]=cleandata->co[i][1];
2366  }
2367  /*doublematrix_control(xy,"xy.txt","xy",PRINT); */
2368 
2369  /* matrix with results */
2370  result=new_doublematrix(n2,cleandata->nch);
2371  initialize_doublematrix(result,0);
2372 
2373  /* doublematrix_dem(result,unused, unused, "result.txt","doublematrix of result",NOPRINT);
2374  doublematrix_dem(cleandata,unused, unused, "result.txt","doublematrix of cleandata",NOPRINT); */
2375  printf("ti2=%f,dt2=%f,n2=%ld,dt1=%f \n",ti2,dt2,n2,dt1);
2376 
2377  /* assign results x- data */
2378  for(i=1;i<=n2;i++){
2379  /* result->co[i][1]=ti2+i*dt2-dt1; */
2380  result->co[i][1]=ti2+dt2/2.+(i-1)*dt2;
2381  }
2382  /*doublematrix_control(result,"result.txt","result",PRINT); */
2383 
2384  /* matrix with intervals */
2385  interval=new_doublematrix(n2+1,2);
2386  initialize_doublematrix(interval,0);
2387 
2388  /* assign intervals x data: points where interpolate data*/
2389  /* interval->co[1][1]=ti2; */
2390  for(i=1;i<=n2+1;i++){
2391  /* interval->co[i][1]=ti2+(i-1)*dt2-dt1; */
2392  interval->co[i][1]=ti2+(i-1)*dt2;
2393 
2394  }
2395  /* doublematrix_control(interval,"prova.interval.txt","interval",PRINT); */
2396 
2397  /* matrix with point for each interval */
2398  trapezio=new_doublematrix(cleandata->nrh/n2+10,2);
2399  initialize_doublematrix(trapezio,0);
2400 
2401  /*doublematrix_control(trapezio,"prova.trapezio.txt","trapezio",PRINT);*/
2402 
2403  perc=0;
2404  count=0;
2405 
2406  /* make calculations for each cols of cleandata */
2407  for(j=2;j<=cleandata->nch;j++){
2408 
2409  /* assign xy y data */
2410  for(i=1;i<=cleandata->nrh;i++){
2411  xy->co[i][2]=cleandata->co[i][j];
2412  }
2413 
2414  /* creates a matrix with regession par */
2415  coeff=interpolating_function(xy);
2416  /*doublematrix_control(coeff,"prova.coeff.txt","coeff",PRINT);
2417  doublematrix_control(xy,"prova.xy.txt","xy",PRINT); */
2418 
2419  /* works for each interval */
2420  for(i=1;i<=n2;i++){
2421 
2422  /* writes controls */
2423 
2424  count++;
2425  if(count==100){
2426  perc=(i+(j-1)*n2)/((float)(cleandata->nch*n2))*100.;
2427  count=0;
2428  printf("computation %4.1f perc\n",perc);
2429  }
2430 
2431  /* count number of xy in the interval and assign x-values to trapezio */
2432  nint=1;
2433  trapezio->co[1][1]=interval->co[i][1];
2434  for(k=1;k<=xy->nrh;k++){
2435  if((xy->co[k][1]>interval->co[i][1])&&(xy->co[k][1]<interval->co[i+1][1])){
2436  nint+=1;
2437  trapezio->co[nint][1]=xy->co[k][1];
2438  }
2439  }
2440  nint+=1;
2441  trapezio->co[nint][1]=interval->co[i+1][1];
2442 
2443  /* assign y-values to trapezio */
2444  for(k=1;k<=nint;k++){
2445  trapezio->co[k][2]=interpolate(trapezio->co[k][1],xy,coeff);
2446  }
2447  /* correct bag of interpolate */
2448  if(trapezio->co[1][1]==xy->co[1][1]){
2449  trapezio->co[1][2]=xy->co[1][2];
2450  }
2451 
2452  /* calculate mean value to trapezio */
2453  result->co[i][j]=mean_function(trapezio, nint);
2454 
2455  /*printf("mean[%d][%d]=%f",i,j,result->co[i][j]);
2456  doublematrix_control(trapezio,"prova.trapezio.txt","trapezio",PRINT); */
2457  }
2458  }
2459 
2460  return result;
2461 }
2462 
2463 
2464 /*-------- quickinterpolate ------------------------------------------------------------------*/
2465 
2466 DOUBLEMATRIX *quickinterpolate(DOUBLEMATRIX *cleandata, short nint)
2467 
2468  /* input: cleandata: double matrix with original data
2469  nint: number of intervals to aggegate
2470  return: double matrix with the interpolated data with (cleandata->nrh)/nint rows */
2471 
2472 {
2473 
2474  long i,j,h,k,row2,row1;
2475  DOUBLEMATRIX *result;
2476 
2477  row1=cleandata->nrh;
2478  row2=row1/nint;
2479 
2480  result=new_doublematrix(row2,cleandata->nch);
2481  initialize_doublematrix(result,0);
2482 
2483  for(j=1;j<=result->nch;j++){
2484  for(i=1;i<=result->nrh;i++){
2485  for(k=1;k<=nint;k++){
2486  h=(i-1)*nint+k;
2487  result->co[i][j]+=cleandata->co[h][j];
2488  }
2489  result->co[i][j]/=(float)nint;
2490  }
2491  }
2492 
2493  return result;
2494 }
2495 
2496 
2497 /*-------- mean_function ------------------------------------------------------------------*/
2498 
2499 double mean_function(DOUBLEMATRIX *data, long n)
2500 
2501 {
2502 
2503  double mean;
2504  long i;
2505 
2506  mean=0;
2507 
2508  for(i=1;i<n;i++){
2509  mean+=(data->co[i+1][2]+data->co[i][2])*
2510  (data->co[i+1][1]-data->co[i][1]);
2511  }
2512 
2513 
2514  return mean/=(2.*(data->co[n][1]-data->co[1][1]));
2515 
2516 }
2517 
2518 /*-------- cleandata_matrix ------------------------------------------------------------------*/
2519 
2520 
2521 DOUBLEMATRIX *cleandata_matrix(DOUBLEMATRIX *data, FLOATVECTOR *V, SHORTMATRIX *control)
2522 
2523 {
2524  DOUBLEMATRIX *cleandata;
2525  double nv, x_i, x_f;
2526  long i,j, n_i, n_f;
2527 
2528  nv=V->co[2];
2529  cleandata=new_doublematrix(data->nrh,data->nch);
2530  initialize_shortmatrix(control,0);
2531  copy_doublematrix(data,cleandata);
2532 
2533  for(j=1;j<=data->nch;j++){
2534  /* if nv in first row...*/
2535  i=1;
2536  n_i=0;
2537  n_f=0;
2538  x_i=0;
2539  x_f=0;
2540  while(i<data->nrh){
2541  while(data->co[i][j]!=nv && i<data->nrh) i++;
2542  x_i=data->co[i-1][j];
2543  n_i=i-1;
2544  if(i==data->nrh) n_i=0;
2545  if(n_i!=0){
2546  while(data->co[i][j]==nv && i<data->nrh) i++;
2547  x_f=data->co[i][j];
2548  n_f=i;
2549  /* call fill_data in datamanipulation.c */
2550  fill_data(cleandata,control,j,n_i,n_f,x_i,x_f);
2551 
2552 
2553  }
2554  }
2555 
2556  /* if nv in last row...*/
2557  }
2558 
2559  return cleandata;
2560 }
2561 
2562 
2563 /*-------- fill_data ------------------------------------------------------------------*/
2564 
2565 void fill_data(DOUBLEMATRIX *cleandata,SHORTMATRIX *control,long j,long n_i,long n_f,double x_i,double x_f)
2566 
2567 {
2568 
2569  double m;
2570  long i,n;
2571 
2572  n=n_f-n_i;
2573  m=(x_f-x_i)/n;
2574 
2575  for(i=1;i<n;i++){
2576  cleandata->co[n_i+i][j]=x_i+m*i;
2577  control->co[n_i+i][j]=1;
2578  }
2579 
2580 }
2581 
2582 /*-------- aggregate ------------------------------------------------------------------*/
2583 
2584 DOUBLEMATRIX *aggregate(DOUBLEMATRIX *data, long col, float nv)
2585 
2586  /* calculates the mean aggregating all data wich have the same value in the column col
2587  input: data: double matrix with original data
2588  long col: the column with,
2589  float nv: novalue
2590  return: double matrix aggregated data as mean values;
2591  in the last column you have the number of aggregated elements */
2592 
2593 {
2594  long i,j,nrows,count,num;
2595  DOUBLEMATRIX *result;
2596  FLOATVECTOR *U;
2597  U=new_floatvector(1);
2598  doublematrix_dem(data,U, U,"data.tmp", "data",NOPRINT);
2599 
2600  /* count rows */
2601  nrows=1;
2602  for(i=2;i<=data->nrh;i++){
2603  if(data->co[i-1][col]!=data->co[i][col]) nrows++;
2604  }
2605  result=new_doublematrix(nrows,data->nch+1);
2606 
2607  /* first row */
2608  for(j=1;j<=data->nch;j++){
2609  if(data->co[1][j]!=nv){
2610  result->co[1][j]=data->co[1][j];
2611  }else{
2612  result->co[1][j]=nv;
2613  }
2614  }
2615  result->co[1][data->nch+1]=1;
2616  /* makes mean */
2617  count=1;
2618  num=1;
2619  for(i=2;i<=data->nrh;i++){
2620  if(data->co[i-1][col]!=data->co[i][col]) {
2621  for(j=1;j<=data->nch;j++){
2622  if(result->co[count][j]!=nv){
2623  result->co[count][j]/=num;
2624  }
2625  if(data->co[i][j]!=nv){
2626  result->co[count+1][j]=data->co[i][j];
2627  }else{
2628  result->co[count+1][j]=nv;
2629  }
2630  }
2631  result->co[count][data->nch+1]=num;
2632  count++;
2633  num=1;
2634 
2635  }else{
2636  for(j=1;j<=data->nch;j++){
2637  result->co[count][j]+=data->co[i][j];
2638  }
2639  num++;
2640  }
2641  }
2642  /* last row */
2643  if(result->co[count][j]!=nv){
2644  for(j=1;j<=data->nch;j++){
2645  result->co[count][j]/=num;
2646  }
2647  }
2648  result->co[count][data->nch+1]=num;
2649 
2650  doublematrix_dem(result,U, U,"result.tmp", "result",NOPRINT);
2651  return result;
2652 }
2653