Version: 1.0.0 (using KDE KDE 3.3.2) Installed from: Debian testing/unstable Packages Compiler: gcc 3.3.5 OS: Linux Here's the plugin source code. ----------- #include <stdlib.h> #include <math.h> #include <stdio.h> #define RADEG 57.29578 #define N_P 3 #define MAX_ITER 10000 int vary_param[N_P] = {1,1,1}; extern "C" int xpol(const double *const inArrays[], const int inArrayLens[], const double is[], double *outArrays[], int outArrayLens[], double outScalars[]); double fit_function (double x, double*p) { return p[0]*(1.0 - (1-p[1])*sin((x-p[2])/RADEG)*sin((x-p[2])/RADEG)); } double Chi2(double*x,double*y,double*p,int npts) { double term; double sum=0.0; int i; for (i=0;i<npts;i++) { term = y[i] - fit_function(x[i],p); sum += term*term; } return sum; } double maximum(double* x,int npts) { int i; double result = -1e99; for (i=0;i<npts;i++) { if (x[i]>result) result = x[i]; } return result; } double minimum(double*x,int npts) { int i; double result = 1e99; for (i=0;i<npts;i++) { if (x[i]<result) result = x[i]; } return result; } int mindex(double*x,int npts) { int i; int result=0; double min=1e99; for (i=0;i<npts;i++) { if (x[i]<min) {min = x[i];result=i;} } return result; } int maxdex(double*x,int npts) { int i; int result=0; double max=0.0; for (i=0;i<npts;i++) { if (x[i]>max) {max = x[i];result=i;} } return result; } double guess_phase(double*x,double*y,int npts){ int min,max,i; min = mindex(y,npts); max = maxdex(y,npts); i = max; return (double)x[i]; } void Search_BruteFit(double*x,double*y,int npts,double*p, double*error,double*covar,int skip) { double *diff; int curpar; double *input,*guess_apr; int done,really_done; int m,n_iter,i=0,j; double a1,a2,a3; double x1,x2,x3,x4; double last_x; diff = (double*)malloc(N_P*sizeof(double)); input = (double*)malloc(N_P*sizeof(double)); guess_apr = (double*)malloc(N_P*sizeof(double)); really_done = 0; for (i=0;i<N_P;i++) { guess_apr[i] = p[i]; diff[i] = p[i]/100.0; } if (!skip) { while (!really_done) { n_iter=0; curpar=0; n_iter++; while ((curpar<N_P)) { // figure out whether to go up or down in parameter space done = 0; x2=Chi2(x,y,guess_apr,npts); guess_apr[curpar]+=diff[curpar]; x3=Chi2(x,y,guess_apr,npts); guess_apr[curpar]-=2.0*diff[curpar]; x1=Chi2(x,y,guess_apr,npts); guess_apr[curpar]+=diff[curpar]; // which way is the direction of decreasing chi^2 ? if (x1<x2) { diff[curpar]*=-1.0;} if (x3<x2) { diff[curpar]*=1.0;} if ((x1<x2)||(x3<x2)){ really_done=0; // keep incrementing parameter until chi^2 increases again do { guess_apr[curpar]+=diff[curpar]; last_x=x2; x2 = Chi2(x,y,guess_apr,npts); if (n_iter>MAX_ITER) {done=1; really_done=1;} } while ((x2<last_x)&&(!done)); } else { really_done=1;} if (!done) { guess_apr[curpar]-=diff[curpar]; if (diff[curpar]<0.0) {diff[curpar]*=-1.0;} // use parabola method to find real minimum in this parameter a1 = guess_apr[curpar]-diff[curpar]; a2 = guess_apr[curpar]; a3 = guess_apr[curpar]+diff[curpar]; for (m=0;m<N_P;m++) { input[m]=guess_apr[m];} input[curpar]=a1; x1=Chi2(x,y,input,npts); input[curpar]=a2; x2=Chi2(x,y,input,npts); input[curpar]=a3; x3=Chi2(x,y,input,npts); guess_apr[curpar]=a3-diff[curpar]*((x3-x2)/(x1-2*x2+x3)+0.5); error[curpar]=2.0/(x1-2*x2+x3); if (error[curpar]>=0.0) { error[curpar]=diff[curpar]*sqrt(error[curpar]); } else { error[curpar]=99.; } } else { error[curpar]=99.; p[curpar]=99.; really_done=1; } curpar++; while (!vary_param[curpar] && (curpar<N_P)) curpar++; if ((!really_done)&&(curpar==N_P)) curpar=0; if (n_iter>MAX_ITER) {really_done=1;error[0]=99.;} } } for (m=0;m<N_P;m++) { p[m]=guess_apr[m]; if (diff[m]<0.0) diff[m]*=-1.0; } // compute co-variance for (m=0;m<N_P;m++){ i = m+1; if (m==N_P) i=0; for (j=0;j<N_P;j++){ input[j]=p[j]; } input[i]+=diff[i]; input[m]+=diff[m]; x1 = Chi2(x,y,input,npts); input[m]-=2.*diff[m]; x2 = Chi2(x,y,input,npts); input[i]-=2.*diff[i]; x4 = Chi2(x,y,input,npts); input[m]+=2.*diff[m]; x3 = Chi2(x,y,input,npts); covar[m]=(x1-x2-x3+x4)/diff[m]/diff[i]/4.0; covar[m]=1.0/covar[m]; } } else { for (i=0;i<N_P;i++) p[i]=guess_apr[i]; } free(guess_apr); free(input); free(diff); } int xpol(const double *const inArrays[], const int inArrayLens[], const double is[], double *outArrays[], int outArrayLens[], double outScalars[]) { double *p; double *error; double *covar; double *x,*y; int npts; int i; npts = inArrayLens[0]; outArrayLens[0] = inArrayLens[0]; outArrays[0]=(double*)realloc(outArrays[0], npts*sizeof(double)); p = (double*)malloc(N_P*sizeof(double)); error = (double*)malloc(N_P*sizeof(double)); covar = (double*)malloc(N_P*sizeof(double)); x = (double*)malloc(npts*sizeof(double)); y = (double*)malloc(npts*sizeof(double)); for (i=0;i<npts;i++){ x[i] = inArrays[0][i]; y[i] = inArrays[1][i]; } // initial guess p[0] = maximum(y,npts); p[1] = minimum(y,npts)/p[0]; p[2] = guess_phase(x,y,npts); printf("guess:%g\n",p[1]); Search_BruteFit(x,y,npts,p,error,covar,0); outScalars[0] = p[1]; printf("result:%g\n",p[1]); for (i=0;i<inArrayLens[0]; i++) outArrays[0][i] = fit_function(x[i],p); free(x); free(y); free(p); free(error); free(covar); printf("at end of program..\n"); return 0; }
ah, I figured out the problem - I'm using mallocs, but only reallocs are allowed. Should have read the documentation a little more closely :) thanks!
On Tuesday 29 March 2005 13:51, Brendan Crill wrote: > ------- ah, I figured out the problem - I'm using mallocs, but only > reallocs are allowed. Should have read the documentation a little more > closely :) thanks! _______________________________________________ Note: you can only use realloc() on memory passed in and out of Kst, but you can use malloc or any other allocator for your own memory.