Bug 102785

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;
}