Joe says  
Hello everyone, I have to fit a surface to data. I'm thinking that a polynomial in x and y with cross terms would suffice (i.e. nothing like exponentials, logs). Would least squares be the easiest way to implement this? I've already pretty much finished an algorithm that does this, but it seems to take forever. This is most likely due to the fact that the input file is over 4 million lines long. 

Total Topic Karma: 0   More by this Author 
Bryan says 


Now I am no graduate student. Hell I am just starting college this fall. But that DOES seem like it would be a BIT long.  
 Author's History  22 June, 2007 
Joe says 


The input file is long or the code would be? The input file is that long because I'm fitting a surface to a 2048x2048 CCD plate. So thats 2048*2048 pixels and each pixel has its own value. It adds up. If you meant the code would be long, it isn't really. The algorithm I made already generates a matrix in one loop. Then takes two more loops to get it in row reduced echelon form. 

 Author's History  22 June, 2007 

deth says 


Depends alot on the Algorithm you wrote. Calculate the time complexity of your algorithm and try to find a way to help your computer, like reducing number of variables and so on. Does your input file really needs to be that long? 

 Author's History  27 June, 2007 
Joe says 


I think that my algorithm is pretty streamlined as far as I can tell. And yes my input file has to be that long. There is no getting around that. These are 2048x2048 images and each pixel has a value. I did, however, find the use of the O2 (that's an "oh") with gcc. When doing a fourth order fit it cut the time from five minutes to under two! 

 Author's History  27 June, 2007 
Valerio says 


I would appreciate if you post some line of pseudocode or else. I'm working on a thesis on Algorithms so I find problems like that fascinating :)  
 Author's History  27 June, 2007 
Joe says 


I can give you the whole thing if you want. I apologize for lack of comments. [code=default] #include <stdio.h> #include <math.h> double tothepower(double base, double exp); int main(int argc, char *argv[]) { double val[4], exp[4]; int i, j, k, kk,flag,midflag,ii,fit; char line[300]; FILE *in, *out; fit=4; double xtx[(fit+1)*(fit+1)][(fit+1)*(fit+1)+1]; in=fopen(argv[1],"r"); k=0; \\ Set everything to zero. for(i=0;i<fit*fit;i++) { for(j=0;j<fit*fit+1;j++) { xtx[i][j]=0; } } xtx[i][j]=0; while (fgets(line,50,in) != NULL) { val[0]=1; sscanf(line,"%lf %lf %lf",&val[1],&val[2],&val[3]); \\ Initialize the array. for (i=0;i<fit*fit;i++) { for(j=0;j<fit*fit;j++) { exp[0]=(ii%fit)/fit; exp[1]=i%fit; exp[2]=(jj%fit)/fit; exp[3]=j%fit; if (exp[0]+exp[1]<fit) { if(exp[2]+exp[3]<fit) { xtx[i][j]=xtx[i][j]+tothepower(val[1],exp[0])*tothepower(val[2],exp[1])*tothepower(val[1],exp[2])*tothepower(val[2],exp[3]); } } } } for (i=0;i<fit*fit;i++) { exp[0]=(ii%fit)/fit; exp[1]=i%fit; if(exp[0]+exp[1]<fit) { xtx[i][fit*fit]=xtx[i][fit*fit]+tothepower(val[1],exp[0])*tothepower(val[2],exp[1])*val[3]; } } } \\ Make matrix upper triangular for(ii=0;ii<fit*fit;ii++) { if (xtx[ii][ii]!=0) { for(j=fit*fit;j>=ii;j) { xtx[ii][j]=xtx[ii][j]/xtx[ii][ii]; for(i=ii+1;i<fit*fit;i++) { xtx[i][j]=xtx[i][j]xtx[i][ii]*xtx[ii][j]; } } } } \\ Make matrix identity for(ii=fit*fit2;ii>=0;ii) { for(i=ii;i>=0;i) { for(j=fit*fit;j>=ii+1;j) { xtx[i][j]=xtx[i][j]xtx[ii+1][j]*xtx[i][ii+1]; } } } out=fopen(argv[2],"w"); \\ Output best fit parameters. for(i=0;i<fit*fit;i++) { fprintf(out,"%25.23e\n",xtx[i][fit*fit]); } } double tothepower(double base, double exp) { int carpow=1; double carry=1; while (carpow<=exp) { carry=carry*base; carpow++; } return carry; } 

 Author's History  27 June, 2007 
deth says 


This code looks so similar to Matlab code! In what language is it C++ ? 

 Author's History  30 July, 2007 