Hi There. I'm trying to fit a function Acos(wt + theta) to my data that I have binned and need to find the best values for the variable parameters A, w and theta so that the cos function fits to my data. However I've tried making a loop for this but it crashes every time I run it. I'm not sure weather this is because my programme uses a ridiculous amount of memory as I try each combination or I'm doing something wrong with my code? I am pretty new to programming so I haven't done this before so have no idea if there is a more efficient method of doing it?

#include <stdio.h>
#include<math.h>
#define N_HITS 1800
#define N_TIME_BINS 78
   int main()
     {
     int i, time_bins[N_TIME_BINS],j,k;
	 float event,time_error[N_TIME_BINS];
     float hit[N_HITS],value;
     FILE *fp, *fw, *fy;
     fp = fopen("damadata.txt" , "r");
     fw = fopen("bindata.txt", "w");
	 fy = fopen("chidata.txt", "w");
            
 for (i=0; i<N_HITS ; i++)  /* enter loop to set array elements to zero*/
     {
       hit[i] = 0;
     }
 for (i=0; i<N_TIME_BINS ; i++)
     {
       time_bins[i] = 0;
     }
            
 if (fp != NULL)
     {                                 
        for(i=0; i<N_HITS; i++)                      /*Enter loop to read all lines of data into hit array*/
          {
          fscanf(fp,"%f\n", &event);                 /*Each line is read into a seperate array element*/
          hit[i]=event; 
          }                                 
     for (i = 0; i < N_HITS; i++)                  /*Now cycle through array containing data elements*/
         {                                          /*Open Loop 1*/                          
            for (j = 0; j < N_TIME_BINS; j++)
            {
               value = hit[i];                       /*Value now equal to element i of array*/                     /*enter loop to bin data into 77 bins of group of 60*/
            if (value >= (j*60) && value < ((j+1)*60))    /*Tests if the value lies within a particular boundary*/
                {
                time_bins[j] = (time_bins[j] + 1);                       /*If it is add 1 element into bin j*/
                }                                  
            else
                { 
                time_bins[j] = (time_bins[j]);                       /*If not element j remains the same*/
                }
           }                                             /*Close Loop 1*/
        }                                             /*Close Loop 2*/
    for (j = 0; j <N_TIME_BINS; j++)                       /*Loop to print the number of data points in each bin*/
          {
		  k=time_bins[j];
          time_error[j]=sqrt(k);
          fprintf(fw, "%d\t%d\t%f\n" , (j*60), k, time_error[j]);       /*Prints the bin number and the number of elements in this bin*/
	      }
   }
   else {
         printf("Error: The file you specified could not be found!");
        }
	int theta;
	double w,a,chi_squ_tot,chi_squ_end[10000000],chi,chi_squ,ideal[N_TIME_BINS];
	j=0;
	chi_squ_tot=0;
	for (a=7; a<14; a=a+0.1){                                        /* Varies my 'a' variable through set of sensible values */
		for(w=0.002;w<0.001; w=w+0.005){                             /* Varies my 'w' variable through set of sensible values */
			for(theta=150; theta<250; theta=theta+1){                 /* Varies my 'a' variable through set of sensible values */
               for(i=0; i<N_TIME_BINS; i++){
					ideal[i]=0;                                       
                   ideal[i]=a*cos((w*(i*60))+theta)+11.11;           /*Calcualtes the flux at t given this set of variables*/
				   chi= (time_bins[i]-ideal[i])/time_error[i];       /*Set of setps to calculate chi squared to evaluate fit*/
					chi_squ=chi*chi; 
					chi_squ_tot=chi_squ_tot + chi_squ;

				}
				chi_squ_end[j]=chi_squ_tot;                          /*Passes final value of chi squared to an external variable*/
				j=j+1;
				chi_squ_tot=0;
			}
		}
	}
	for (j=0; j<10000000; j++){
	fprintf(fy, "%d\t%f\n" , j,chi_squ_end[j] );                     /*Prints all values of chi squared so can see which produces best fit*/
		}
   
          fclose(fp);
          fclose(fw);
	      fclose(fy);
     return(0);
}

Recommended Answers

All 10 Replies

Without generating the data myself I'm kind of limited in what I can do to help.
However, chi_squ_end[10000000] bugs me. A lot. And being a double type, on a 32-bit system, that thing's 80 Million bytes huge.
May I humbly suggest moving outside and before the main block. This takes it outside the "stack" and places it on the heap. Honestly, I have no idea how well it'll work.
But after doing work on small memory devices and discovering the "joys" of sprintf in confined spaces, I've learned to become deathly wary of putting big items on the stack.

What does this do and why is it doing it that way?

for(w=0.002;w<0.001; w=w+0.005)

>> However, chi_squ_end[10000000] bugs me. A lot.

Ditto. I'm not sure whether this is complete, but at least the way you have it now, I see no reason to have it as an array. It seems like you can simply substitute a single variable for it. As far as I can tell, change line 78 to this and move it after line 71

fprintf(fy, "%d\t%f\n" , j,chi_squ_end );

delete line 77, change line 71 to this...

chi_squ_end=chi_squ_tot;

and change the declaration on line 57 to a double rather than an array, and you have the exact same code, correct? Again, this only deals with the code as it is now. Perhaps you actually need an array and are adding to the program, but as the program stands now, there are no calculations based on this array and all you do is print, so I see no reason to store anything but the most recent value. It can be discarded after printing.

Thanks for your replies, I have lectures now but I'll give this a go and report back. Thanks again

OK, so I've got rid of the huge array, but now I am getting no output written to my chi.txt file when I should be getting the values of a, w, theta and chi squared for each combination and I have no idea why!?!?

#include <stdio.h>
#include<math.h>
#define N_HITS 1800
#define N_TIME_BINS 78
   int main()
     {
     int i, time_bins[N_TIME_BINS],j,k;
	 float event,time_error[N_TIME_BINS];
     float hit[N_HITS],value;
     FILE *fp, *fw, *fy;
     fp = fopen("damadata.txt" , "r");
     fw = fopen("bindata.txt", "w");
	 fy = fopen("chidata.txt", "w");
            
 for (i=0; i<N_HITS ; i++)  /* enter loop to set array elements to zero*/
     {
       hit[i] = 0;
     }
 for (i=0; i<N_TIME_BINS ; i++)
     {
       time_bins[i] = 0;
     }
 if (fy == NULL) {
         printf("Error: The file you specified could not be found!");
        } 
 if (fp != NULL)
     {                                 
        for(i=0; i<N_HITS; i++)                      /*Enter loop to read all lines of data into hit array*/
          {
          fscanf(fp,"%f\n", &event);                 /*Each line is read into a seperate array element*/
          hit[i]=event; 
          }                                 
     for (i = 0; i < N_HITS; i++)                  /*Now cycle through array containing data elements*/
         {                                          /*Open Loop 1*/                          
            for (j = 0; j < N_TIME_BINS; j++)
            {
               value = hit[i];                       /*Value now equal to element i of array*/                     /*enter loop to bin data into 77 bins of group of 60*/
            if (value >= (j*60) && value < ((j+1)*60))    /*Tests if the value lies within a particular boundary*/
                {
                time_bins[j] = (time_bins[j] + 1);                       /*If it is add 1 element into bin j*/
                }                                  
            else
                { 
                time_bins[j] = (time_bins[j]);                       /*If not element j remains the same*/
                }
           }                                             /*Close Loop 1*/
        }                                             /*Close Loop 2*/
    for (j = 0; j <N_TIME_BINS; j++)                       /*Loop to print the number of data points in each bin*/
          {
		  k=time_bins[j];
          time_error[j]=sqrt(k);
          fprintf(fw, "%d\t%d\t%f\n" , (j*60), k, time_error[j]);       /*Prints the bin number and the number of elements in this bin*/
	      }
   }
   else {
         printf("Error: The file you specified could not be found!");
        }
	int theta;
	double w,a,chi_squ_tot,chi,chi_squ,ideal;
	j=0;
	chi_squ_tot=0;
	for (a=7; a<14; a=a+0.1){                                        /* Varies my 'a' variable through set of sensible values */
		for(w=0.002;w<0.001; w=w+0.005){                             /* Varies my 'w' variable through set of sensible values */
			for(theta=150; theta<250; theta=theta+1){                 /* Varies my 'a' variable through set of sensible values */
               for(i=0; i<N_TIME_BINS; i++){                                      
                   ideal=a*cos((w*(i*60))+theta)+11.11;           /*Calcualtes the flux at t given this set of variables*/
				   chi= (time_bins[i]-ideal)/time_error[i];       /*Set of setps to calculate chi squared to evaluate fit*/
					chi_squ=chi*chi; 
					chi_squ_tot=chi_squ_tot + chi_squ;

				}
				fprintf(fy, "%f\t%f\t%d\t%f\n" , a, w, theta, chi_squ_tot );
				chi_squ_tot=0;
			}
		}
	}
   
          fclose(fp);
          fclose(fw);
	      fclose(fy);
     return(0);
}

Thanks!

Line 63. Seems like this loop will never execute. w is initialized as .002. Loop condition is that w is less than .001, so it's false before the code inside executes even once, which means you can delete the loop and it'll have no effect.

for(w=0.002;w<0.001; w=w+0.005)

ah yea, cheers!

Is there a limit to the number of loops you can have inside eachother as my bit of code

for (a=7; a<14.5; a=a+0.5){                                        /* Varies my 'a' variable through set of sensible values */
		for(w=0.01;w<0.02; w=w+0.002); {                           /* Varies my 'w' variable through set of sensible values */
			for(theta=150; theta<201; theta=theta+1){                 /* Varies my 'a' variable through set of sensible values */
               for(i=0; i<N_TIME_BINS; i++){                                      
                   ideal[i]=(a*cos((w*(i*60))+theta))+11.11;           /*Calcualtes the flux at t given this set of variables*/
				   chi[i]= ((time_bins[i]-ideal[i])/time_error[i]);       /*Set of setps to calculate chi squared to evaluate fit*/
				   chi_squ[i]=chi[i]*chi[i]; 
				   chi_squ_tot=chi_squ_tot + chi_squ[i];
				   ideal[i]=0;
                   chi_squ[i]=0;
		           chi[i]=0;

				}
				re_chi=((chi_squ_tot)/(N_TIME_BINS-4));
				fprintf(fy, "%f\t%f\t%d\t%f\t%f\n" ,a,w,theta, chi_squ_tot,re_chi );
				chi_squ_tot=0;
			}
		  }
		}

Will cycle through all the a and theta values fine but it takes the w value as constant even though I want to change this variable as well. When I print the data to a file I have values for A fro 7-14, theta from 150->200 for each A value but the only w value I get is 0.02, it doesn't build up from 0.01

>> Is there a limit to the number of loops you can have inside eachother as my bit of code

No.


>> but it takes the w value as constant

Actually, it doesn't. It loops through. More on that later.


>> the only w value I get is 0.02

Note that this value(.02) is the upper end of your w loop. Try changing the loop from .02 to .03 and see what happens. Pattern?

Now try changing the loop in line 2 to this...

for(int w=0.01;w<0.02; w=w+0.002)

What happens? Add that all up and what it tells you is that the loop goes through all the way, but the code that the loop executes isn't what you want it to be. First check your brackets.

They seem to match, except actually they don't. Look closely, really closely, at line 2.

Do the lines just differ by the 'int' deceleration? I thought I had to use float as I want using integer values for w?

Whoops. My bad. Yes, float or double, not int. I didn't notice that and didn't intend that and wasn't paying attention to that because it's actually irrelevant to my point(though obviously not irrelevant!). The problem is a scoping issue and the scoping issue comes right here...

for(w=0.01;w<0.02; w=w+0.002); {

The semicolon screws everything up. The example was supposed to be this...

for(float w=0.01;w<0.02; w=w+0.002); {

or double, whichever you prefer. You'd notice that your printouts would NOT print 0.02.


I have to proofread better. Here I am giving an example trying to lead you to the problem and instead it leads you in the wrong direction. Sorry about that. Get rid of the semicolon.

Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, networking, learning, and sharing knowledge.