//=============================================================//
//     Numerical Analysis Package in Java                      //
//     Non linear least square by Levenberg-marquart method    //
//                                                             //
//     Author : Dr. Turhan Coban                               //
//     N106 WSEC, University of Nebraska- Lincoln              //
//     phone : 1-402-472 3793                                  //
//     email : MCOBAN2@unl.edu                                 //
//             MCOBAN2@netscape.net                            //
//=============================================================//
import java.io.*;

abstract class f_xa
{
  abstract double func(double x,double a[]);
}



public class marquart
{
   public static double X2(f_xa y,double xi[], double yi[], double a[]) 
   {
     int i;
     double sum = 0;
     double dif;
     double yf;
     for (i = 0; i < xi.length; i++) 
     {
       yf=y.func(xi[i], a);
	   dif = yi[i] - yf;
       //System.out.println("i="+i+"xi="+xi[i]+"yi="+yi[i]+"yf="+yf+"dif="+dif);
       sum += dif * dif;
     }
     return sum;
   }

   public static double alpha(f_xak dy,int k, int l, double xi[], double yi[], double a[]) {
     int i;
     double sum = 0;
     double dif;
     for (i = 0; i < xi.length; i++) {
       dif = dy.func(xi[i], a, k) * dy.func(xi[i], a, l);
       sum += dif;
     }
     return sum;
   }

   public static double beta(f_xa y,f_xak dy,int k, double xi[], double yi[], double a[]) {
     int i;
     double sum = 0;
     double dif;
     for (i = 0; i < xi.length; i++) {
       dif = (yi[i] - y.func(xi[i], a)) * dy.func(xi[i], a, k);
       sum += dif;
     }
     return sum;
   }
    
   public static double[] B(f_xa y,f_xak dy, double xi[], double yi[], double a[])
   {
	int n=a.length;
	double b[]=new double[n];   
    for(int k=0;k<n;k++)
       b[k]=beta(y,dy,k,xi,yi,a);
    return b;
   }
   
   public static double alphap(f_xak dy,int k, int l, double xi[], double yi[], double a[], double lambda) {
     if (k == l)return alpha(dy,k, l, xi, yi, a) * (1.0 + lambda);
     else return alpha(dy,k, l, xi, yi, a);
   }
   
   public static double[][] A(f_xak dy, double xi[], double yi[], double a[],double lambda)
   {
   int n=a.length;
   double c[][]=new double[n][n];
   for(int k=0;k<n;k++)
      for(int l=0;l<n;l++)
         c[k][l]=alphap(dy,k,l,xi, yi,a,lambda);
   return c;
   }

   public static double[] da(f_xa y,f_xak dy,double xi[], double yi[], double a[],double lambda) 
   {
     int M = a.length;
     double AA[][]=A(dy,xi,yi,a,lambda);
     double BB[]=B(y,dy,xi,yi,a);
     double dd[]=new double[M];
     dd=Matrix.AXB(AA,BB);
     return dd;
   }
   
   public static double[] aa(f_xa y,f_xak dy,double xi[], double yi[], double a[])
   {
   int M = a.length;
     double tot=0;
     double tolerance=1e-4;
     double k2,k2a;
     double lambda=0.001;
     double dd[]=new double[M];
     double b[]=new double[M];
     k2=X2(y,xi,yi,a);
     do
     {
     double AA[][]=A(dy,xi,yi,a,lambda);
     double BB[]=B(y,dy,xi,yi,a);
     dd=Matrix.AXB(AA,BB);
     b=Matrix.add(a,dd);
     k2a=X2(y,xi,yi,b);
     tot=0;
     for(int i=0;i<M;i++) tot+=dd[i];
     tot=Math.abs(tot);
     if(k2a>=k2) lambda*=10.478967546;
     else {lambda/=9.587975;a=Matrix.add(a,dd);k2=k2a;}
     System.out.println("k2="+k2+"k2a="+k2a+"tot="+tot);   
     } while(tot>tolerance);
     return a;
   }
 }
