package ChaosDemos;
import java.util.*;

public class svdfitClass {

double[][] u;
double[][] v;
double[] w;
double[] xminus;
double[] yminus;
static double TOL=1.0e-5;
int ma;
int ndata; 

public svdfitClass(int inNdata, int inMa) {
    ndata=inNdata;
    ma=inMa;
    u=new double[ndata+1][ma+1];
    v=new double[ma+1][ma+1];
    w=new double[ma+1];
}    
                                                    
public double svdfit(int[] ix, double y[],  double inx[], double iny[],
                      double sig[],  double a[])
{
    int j,i;
    double wmax,tmp,thresh,sum;
    double chisq;

    double[] b=new double[ndata+1];
    double[] afunc=new double[ma+1];
    xminus=inx;
    yminus=iny;

    for (i=1;i<=ndata;i++) {
        funcs(ix[i],afunc,ma);
        tmp=1.0/sig[i];
        for (j=1;j<=ma;j++) u[i][j]=afunc[j]*tmp;
        b[i]=y[i]*tmp;
    }
    svdcmp(ndata,ma);
    wmax=0.0;
    for (j=1;j<=ma;j++)
        if (w[j] > wmax) wmax=w[j];
    thresh=TOL*wmax;
    for (j=1;j<=ma;j++)
        if (w[j] < thresh) w[j]=0.0;
    svbksb(ndata,ma,b,a);
    chisq=0.0;
    for (i=1;i<=ndata;i++) {
        funcs(ix[i],afunc,ma);
        for (sum=0.0,j=1;j<=ma;j++) sum += a[j]*afunc[j];
        chisq += Math.pow((y[i]-sum)/sig[i],2);
    }
    return chisq;
}

private void svdcmp(int m, int n) {

    int i,its,j,jj,k;
    //Got error messages if these two weren't initialized: don;t know why
    int l=0;
    int nm=0;
    boolean flag;
    double anorm,c,f,g,h,s,scale,x,y,z;

    double[] rv1=new double[n+1];
    g=scale=anorm=0.0;
    for (i=1;i<=n;i++) {
        l=i+1;
        rv1[i]=scale*g;
        g=s=scale=0.0;
        if (i <= m) {
            for (k=i;k<=m;k++) scale += Math.abs(u[k][i]);
            if (scale!=0.) {
                for (k=i;k<=m;k++) {
                    u[k][i] /= scale;
                    s += u[k][i]*u[k][i];
                }
                f=u[i][i];
                g = -SIGN(Math.sqrt(s),f);
                h=f*g-s;
                u[i][i]=f-g;
                for (j=l;j<=n;j++) {
                    for (s=0.0,k=i;k<=m;k++) s += u[k][i]*u[k][j];
                    f=s/h;
                    for (k=i;k<=m;k++) u[k][j] += f*u[k][i];
                }
                for (k=i;k<=m;k++) u[k][i] *= scale;
            }
        }
        w[i]=scale *g;
        g=s=scale=0.0;
        if (i <= m && i != n) {
            for (k=l;k<=n;k++) scale += Math.abs(u[i][k]);
            if (scale!= 0.) {
                for (k=l;k<=n;k++) {
                    u[i][k] /= scale;
                    s += u[i][k]*u[i][k];
                }
                f=u[i][l];
                g = -SIGN(Math.sqrt(s),f);
                h=f*g-s;
                u[i][l]=f-g;
                for (k=l;k<=n;k++) rv1[k]=u[i][k]/h;
                for (j=l;j<=m;j++) {
                    for (s=0.0,k=l;k<=n;k++) s += u[j][k]*u[i][k];
                    for (k=l;k<=n;k++) u[j][k] += s*rv1[k];
                }
                for (k=l;k<=n;k++) u[i][k] *= scale;
            }
        }
        anorm=Math.max(anorm,(Math.abs(w[i])+Math.abs(rv1[i])));
    }
    for (i=n;i>=1;i--) {
        if (i < n) {
            if (g!=0.) {
                for (j=l;j<=n;j++)
                    v[j][i]=(u[i][j]/u[i][l])/g;
                for (j=l;j<=n;j++) {
                    for (s=0.0,k=l;k<=n;k++) s += u[i][k]*v[k][j];
                    for (k=l;k<=n;k++) v[k][j] += s*v[k][i];
                }
            }
            for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;
        }
        v[i][i]=1.0;
        g=rv1[i];
        l=i;
    }
    for (i=Math.min(m,n);i>=1;i--) {
        l=i+1;
        g=w[i];
        for (j=l;j<=n;j++) u[i][j]=0.0;
        if (g!=0.) {
            g=1.0/g;
            for (j=l;j<=n;j++) {
                for (s=0.0,k=l;k<=m;k++) s += u[k][i]*u[k][j];
                f=(s/u[i][i])*g;
                for (k=i;k<=m;k++) u[k][j] += f*u[k][i];
            }
            for (j=i;j<=m;j++) u[j][i] *= g;
        } else for (j=i;j<=m;j++) u[j][i]=0.0;
        ++u[i][i];
    }
    for (k=n;k>=1;k--) {
        for (its=1;its<=30;its++) {
            flag=true;
            for (l=k;l>=1;l--) {
                nm=l-1;
                if ((Math.abs(rv1[l])+anorm) == anorm) {
                    flag=false;
                    break;
                }
                if ((Math.abs(w[nm])+anorm) == anorm) break;
            }
            if (flag) {
                c=0.0;
                s=1.0;
                for (i=l;i<=k;i++) {
                    f=s*rv1[i];
                    rv1[i]=c*rv1[i];
                    if ((Math.abs(f)+anorm) == anorm) break;
                    g=w[i];
                    h=pythag(f,g);
                    w[i]=h;
                    h=1.0/h;
                    c=g*h;
                    s = -f*h;
                    for (j=1;j<=m;j++) {
                        y=u[j][nm];
                        z=u[j][i];
                        u[j][nm]=y*c+z*s;
                        u[j][i]=z*c-y*s;
                    }
                }
            }
            z=w[k];
            if (l == k) {
                if (z < 0.0) {
                    w[k] = -z;
                    for (j=1;j<=n;j++) v[j][k] = -v[j][k];
                }
                break;
            }
            if (its == 30) System.out.println("no convergence in 30 svdcmp iterations");
            x=w[l];
            nm=k-1;
            y=w[nm];
            g=rv1[nm];
            h=rv1[k];
            f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
            g=pythag(f,1.0);
            f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
            c=s=1.0;
            for (j=l;j<=nm;j++) {
                i=j+1;
                g=rv1[i];
                y=w[i];
                h=s*g;
                g=c*g;
                z=pythag(f,h);
                rv1[j]=z;
                c=f/z;
                s=h/z;
                f=x*c+g*s;
                g = g*c-x*s;
                h=y*s;
                y *= c;
                for (jj=1;jj<=n;jj++) {
                    x=v[jj][j];
                    z=v[jj][i];
                    v[jj][j]=x*c+z*s;
                    v[jj][i]=z*c-x*s;
                }
                z=pythag(f,h);
                w[j]=z;
                if (z!=0.) {
                    z=1.0/z;
                    c=f*z;
                    s=h*z;
                }
                f=c*g+s*y;
                x=c*y-s*g;
                for (jj=1;jj<=m;jj++) {
                    y=u[jj][j];
                    z=u[jj][i];
                    u[jj][j]=y*c+z*s;
                    u[jj][i]=z*c-y*s;
                }
            }
            rv1[l]=0.0;
            rv1[k]=f;
            w[k]=x;
        }
    }
}

private void svbksb(int m, int n, double[] b, double[] x)
{
    int jj,j,i;
    double s;
    double[] tmp;

    tmp=new double[n+1];
    for (j=1;j<=n;j++) {
        s=0.0;
        if (w[j]!=0.) {
            for (i=1;i<=m;i++) s += u[i][j]*b[i];
            s /= w[j];
        }
        tmp[j]=s;
    }
    for (j=1;j<=n;j++) {
        s=0.0;
        for (jj=1;jj<=n;jj++) s += v[j][jj]*tmp[jj];
        x[j]=s;
    }
}


private double SIGN(double a, double b) {
    if(b>=0) return Math.abs(a);
    else return -Math.abs(a);
}

private double pythag(double a, double b)
{
    double absa,absb;
    absa=Math.abs(a);
    absb=Math.abs(b);
    if(absa > absb)
      return absa*Math.sqrt(1.+Math.pow((absb/absa),2));
    else
      if(absb==0.)
        return 0.;
      else
        return absb*Math.sqrt(1.+Math.pow((absa/absb),2));
    }

    private void funcs(int ix, double[] afunc, int ma)  {
        afunc[1]=xminus[ix];
        afunc[2]=yminus[ix];
        afunc[3]=1.;
    }
    
}