//
//     framework.java
//
// M.Lampton UCB SSL 2003:   mlampton@SSL.berkeley.edu  
//
//
// 
// One systematic parameter + 3*NFILT filter parms; NFILT=9 filters
// 
//  >> assume that filter centerLam and widthLam are *given*
//   not determined by joint stellar photometry
//
//  >> assume that filter grasp is totally unknown, except for
//  its concordance with stellar photometry.
//
// Estimates filter parm errors from the one systematic error 
// plus a collection of photometry errors, e.g. standards
//
//   ** this example: all photometry errors are equal, 1%rms
//   ** this example: systematic error 1%rms linear in wavelength
//
//
// Method: Parameters => Jacobian => Fisher matrix => covariance matrix
//
// Datum 0 is one SYSCAL observation of overall slope coefficient.
//   this datum does not involve any of the filters. 
//   this datum is obtained outside the system of filter determination.
//   its SNR is set by SYSSNR = 100.0 and SIGMA=1.0
//
// Data 1...10 are ten photometry stars seen through filter 1.
// data 11..20 the same photometry stars in filter 2. 
// data 21..30 the same photometry stars in filter 3. etc
//
// Parameter 0 is overall systematic slope parameter
// Parameters 1,2,3 are filter1 throughput, lambdaPeak, width. 
// parameters 4,5,6 are filter2 throughput, lambdaPeak, width.
// parameters 7,8,9 are filter3. etc
//
// 
//
// This code evaluates fluxes, Jacobian, Fisher matrix, covariance matrix.
//
//
///////////// A Stellar Spectral Flux Library, 1150-25000A ////////////////////////
//
// A.J. Pickles, PASP 110, 863, 1998
// http://www.ifa.hawaii.edu/~pickles/AJP/hilib.html
//
//******************************************************************************
// HILIB Stellar Flux library contents, procedures and summary. HILIB consists of:
// UVILIB (1150-10620/5A, 131 spectra),
// UVKLIB (1150-25000/5A, 66 with ir spectra, rest with JHK continuum),
// UVMLIB (1150-52000/5A,  16 spectra).
//******************************************************************************
//
// The first column is wavelength in Angstrom, always 1150 - 25000A in steps of 5A. 
// The second column is F(lambda) for the spectrum, normalised to unity at 5556A. 




import java.io.*;   
import java.util.*;
import java.text.DecimalFormat;


public class framework
{
    static PrintWriter pw; 
    static FileOutputStream fos; 

    ////////// photometry and spectra //////////////////////

    static final int    NFILTS =      9;                // how many filters
    static final int    NPPF   =      3;                // how many parms per filter 
    static final int    NSTARS =      10;               // how many stellar spectra
    static final int    NSYS   =      1;                // how many systematic parms

    static final int    NPTS   =  NSYS + NSTARS*NFILTS; // total data points
    static final int    NPARMS =  NSYS + NPPF*NFILTS;   // total parameters
    static final int    NADJ   =  10;                   // how many parms are active

    static final int[]  which  = {0,1,4,7,10,13,16,19,22,25}; // sys + 9 throughputs

    static final double SYSSNR =  100;                  // systematic slope SNR
    static final double PHOSNR =  100;                  // photom SNR, each star, each band 
    static final double SIGMA = 1.0;                    // std deviation everywhere

    static final int    NWAVES =   2501;                // how many wavelengths per star
    static final double DELTAP =   1e-6;                // for two sided derivatives

    static double[][]   spectra= new double[NSTARS][NWAVES];  // photons/nm each star
    static double[][]   coefs  = new double[NSTARS][NFILTS];  // makes SNRs = 100.0
    static double[]     error  = new double[NPTS];      // the error vector

    //////////////////// create the true parameters ////////////////////////////

    static final double[] truep = {0.0,                  // systematic slope
                                   1.0, 0.42, 0.095,     // filter at 0.42 um 
                                   1.0, 0.48, 0.11,
                                   1.0, 0.56, 0.13,
                                   1.0, 0.64, 0.15,
                                   1.0, 0.73, 0.17,
                                   1.0, 0.85, 0.20,
                                   1.0, 0.97, 0.23,
                                   1.0, 1.12, 0.26,
                                   1.0, 1.29, 0.39};      // filter at 1.4 um
    


    static final String picklenames_three[] = { 
          "UKO5V.DAT", 
          "UKK0III.DAT", 
          "UKM7III.DAT" };

    static final String picklenames_ten[] = {
          "UKO5V.DAT", 
          "UKB9V.DAT", 
          "UKA7V.DAT",
          "UKF5V.DAT", 
          "UKK0III.DAT", 
          "UKK5III.DAT", 
          "UKK7V.DAT",
          "UKM2V.DAT", 
          "UKM5III.DAT",
          "UKM7III.DAT" };


    static double T10[] = {3000,4000,5000,6000,8000,10000,15000,20000,40000,80000};
    static double T4[]  = {3000, 5000, 10000, 80000}; 
    static final String plancknames_four[] = {"3000",  "5000", "10000", "80000"};


    static void setupspectrum_Planck(int istar, double t )
    // sets up one spectrum, temperature t kelvins.
    // Planck: 2*h*c2/lam^5 * 1/(exp(hc/(kt*lam))-1)  joules/m wavel
    // Planck: 2*c/lam^4 * 1/(exp(hc/(kt*lam))-1)    photons/m wavel
    // but here we normalize the overall coefficient to unity at 1000nm
    // Note:  hc/k = 0.0144 meter-kelvins = 14400 micron kelvins
    {
           for (int nm=0; nm<NWAVES; nm++)
           {
              double microns = 0.001 * nm;  
              if (nm<150)
                spectra[istar][nm] = 0.0;
              else
                spectra[istar][nm] = Math.pow(microns,-4.0)/(Math.exp(14400.0/(t*microns))-1.0);
           }

           // now normalize the spectra at 1000nm, to one photon/nm
           double s1000 = spectra[istar][1000];
           for (int nm=1; nm<NWAVES; nm++)
             spectra[istar][nm] *= 1.0/s1000;
    }


    static void setupspectra_Planck( )
    // sets up NSTARS spectra, using temperature array Tn; here n=4
    {
        for (int i=0; i<NSTARS; i++)
          setupspectrum_Planck(i, T4[i]); 
    }    

    static void setupspectra_Pickles( ) throws IOException
    // http://www.ifa.hawaii.edu/~pickles/AJP/hilib.html "UVKLIB"
    // for stars with given flambda, need photons ~ flambda * lambda
    {
        for (int i=0; i<NSTARS; i++)
        {
            getstar(picklenames_ten[i], spectra[i]); 

            // convert from ergs/nm to photons/nm
            for (int nm=1; nm<NWAVES; nm++)
              spectra[i][nm] *= nm; 

            // get peak intensity
            double smax = 0.0; 
            for (int nm =1; nm<NWAVES; nm++)
               if (spectra[i][nm] > smax)
                 smax = spectra[i][nm]; 

            // normalize at peak
            for (int nm=1; nm<NWAVES; nm++)
               spectra[i][nm] *= 0.1/smax;  
        }
    }


    static int getstar(String fname, double flux[]) throws IOException
    {
       int nrecords=0; 
       File myfile = new File(fname); 
       if (myfile.exists() && myfile.canRead())   // LJ p.308
         try
         {
            FileReader fr = new FileReader(myfile); 
            BufferedReader br = new BufferedReader(fr); 
            String line;
            while ((line = br.readLine()) != null)
            {
               nrecords++;           
               int nfield = 0;
               int i=0; 
               double d=0.0; 
               StringTokenizer st = new StringTokenizer(line); 
               while (st.hasMoreTokens())
               {
                  String nt = st.nextToken(); 
                  try
                  {
                     d = Double.parseDouble(nt); 
                  }
                  catch (NumberFormatException e)
                  {
                     d = 0.0;
                  }
                  if (nfield == 0)
                    i = (int) (0.1*d + 0.5); 
                  if ((nfield == 1) && (i>0) && (i<2500))
                    flux[i] = d;  
                  nfield++;
               }
            }
         }
         catch (FileNotFoundException e)
         {
             System.out.println(fname + " is unavailable"); 
             beep(); 
         }
         return nrecords;
    }


    static void prenormalize()
    // initializes the array of coefficients to unity before initial photometry
    {
        for (int istar=0; istar<NSTARS; istar++)
          for (int jfilt=0; jfilt<NFILTS; jfilt++)
             coefs[istar][jfilt] = 1.0; 
    }

    static void renormalize()
    // modifies the coefficients to give PHOSNR after initial photometry
    {
        for (int istar=0; istar<NSTARS; istar++)
          for (int jfilt=0; jfilt<NFILTS; jfilt++)
          { 
             int idat = NSYS + istar + jfilt*NSTARS; 
             coefs[istar][jfilt] *= PHOSNR / f(idat, truep);
          }
             
    }


    static double chuckb(double microns)
    // Lampton's take on Chuck Bower's B filter function
    // only here I want a single point per call
    // peak = 1.000 is at 0.42 microns
    // integral chuckb dlam = 0.095139 um = 0.22652 * peakLambda
    // HM at 0.3900 and 0.4845 um; FWHM = 0.0945 um.
    {
        double nm = 1000.0*microns; 
        if (nm < 360.0)
          return 0.0; 
        if (nm > 560.0)
          return 0.0; 
        if (nm < 420.0)
          return 1./(1. + Math.exp(-0.17*(nm-390.0))) + 0.006*(nm-390.0)/30.0; 
        double cosfun = Math.cos(1.57079633*(nm-420.0)/140.0); 
        return Math.pow(cosfun, 2.4); 
    }


    static double filter(double microns, double p[])
    //  chuckb filter form, made tunable: 3 filter parms
    //  M.Lampton UCB SSL 2003
    //       p[0] = integrated filter throughput
    //       p[1] = filter peak microns;
    //       p[2] = filter band width, nominally 0.2262*peakmicrons. 
    {
        double arg = 0.42 + (microns - p[1]) * (0.0945/p[2]); 
        return chuckb(arg) * (p[0]/p[2]); 
    } 

    static void setupNineFilters()
    // sets up parameters for nine photometric filters
    {
        for (int i=0; i<9; i++)
        {
           double wavel = 0.42 * Math.pow(1.15, (double) i); 
           double width = 0.234 * wavel; 
           int index = 1 + 3*i; 
           truep[index] = 1.0; 
           truep[index+1] = wavel; 
           truep[index+2] = width; 
        }
    } 
                

    static double systematic(double microns, double parms[])
    // systematic linear wavelength-dependent calibration error
    // uses only parms[0] for overall slope; that is, NSYS=1
    {
        return 1.0 + microns * parms[0]; 
    }


    static double f(int idat, double p[])
    // noiseless photometry for one filt, one star, given parameters p[]
    {
        if (idat < NSYS)
          return 1.0;

        // determine which star, which filter, is this datum
        int istar = (idat - NSYS) % NSTARS; 
        int jfilt = (idat - NSYS) / NSTARS;  

        // get this filter's parms from our parm vector...
        double tp[] = new double[NPPF]; 
        for (int j=0; j<NPPF; j++)
          tp[j] = p[NSYS + j + jfilt*NPPF]; 

        double sum = 0.0; 
        for (int nm=0; nm<NWAVES; nm++)
        {
            double microns = 0.001 * nm; 
            sum += spectra[istar][nm] * systematic(microns, p) * filter(microns, tp); 
        }
        return sum * coefs[istar][jfilt]; 
    }


    static double dfdp(int row, int col, double p[])
    // returns one element of the Jacobian; 
    // two-sided derivative, dPhotom/dParameter
    {
        if ((row==0) && (col==0))
          return 100.0; 
        if (row==0)
          return 0.0; 
        double pminus[] = new double[NPARMS]; 
        double pplus[] = new double[NPARMS];
        for (int j=0; j<NPARMS; j++)
          pminus[j] = pplus[j] = p[j]; 

        int jparm = which[col]; 
        pminus[jparm] -= DELTAP;
        pplus[jparm] += DELTAP; 
        return (f(row, pplus) - f(row, pminus))/(2.0*DELTAP); 
    }


    static void makejacobian(double p[], double jacobian[][])
    // builds the entire jacobian matrix
    {
        for (int row=0; row<NPTS; row++)
          for (int col=0; col<NADJ; col++)
            jacobian[row][col] = 0.0; 

        jacobian[0][0] = SYSSNR;   // the systematic SNR

        for (int row=1; row<NPTS; row++)
          for (int col=0; col<NADJ; col++)
            jacobian[row][col] = dfdp(row, col, p); 
    }
                        
        
    static double gaussj( double a[][], int N )
    // inverts the double array a[N][N] by Gauss-Jordan method
    {
        double det = 1.0, big, save;
        int i,j,k,L;
        int[] ik = new int[100];
        int[] jk = new int[100];
        for (k=0; k<N; k++)
        {
            big = 0.0;
            for (i=k; i<N; i++)
              for (j=k; j<N; j++)          // find biggest element
                if (Math.abs(big) <= Math.abs(a[i][j]))
                {
                    big = a[i][j];
                    ik[k] = i;
                    jk[k] = j;
                }
            if (big == 0.0) return 0.0;
            i = ik[k];
            if (i>k)
              for (j=0; j<N; j++)          // exchange rows
              {
                  save = a[k][j];
                  a[k][j] = a[i][j];
                  a[i][j] = -save;
              }
            j = jk[k];
            if (j>k)
              for (i=0; i<N; i++)
              {
                  save = a[i][k];
                  a[i][k] = a[i][j];
                  a[i][j] = -save;
              }
            for (i=0; i<N; i++)            // build the inverse
              if (i != k)
                a[i][k] = -a[i][k]/big;
            for (i=0; i<N; i++)
              for (j=0; j<N; j++)
                if ((i != k) && (j != k))
                  a[i][j] += a[i][k]*a[k][j];
            for (j=0; j<N; j++)
              if (j != k)
                a[k][j] /= big;
            a[k][k] = 1.0/big;
            det *= big;                    // bomb point
        }                                  // end k loop
        for (L=0; L<N; L++)
        {
            k = N-L-1;
            j = ik[k];
            if (j>k)
              for (i=0; i<N; i++)
              {
                  save = a[i][k];
                  a[i][k] = -a[i][j];
                  a[i][j] = save;
              }
            i = jk[k];
            if (i>k)
              for (j=0; j<N; j++)
              {
                  save = a[k][j];
                  a[k][j] = -a[i][j];
                  a[i][j] = save;
              }
        }
        return det;
    }


    static void getfisher(double jacobian[][], double error[], double fisher[][])
    // supplies Fisher matrix to a caller-supplied array[NADJ][NADJ]
    // Depends only on the Jacobian and on the errors of each datum. 
    // Method: Fjk = sumi[Jij * Jik * Wi]
    {
        for (int j=0; j<NADJ; j++)        
          for (int k=0; k<NADJ; k++)
          {
             fisher[j][k] = 0.0;
             for (int i=0; i<NPTS; i++)
               fisher[j][k] += jacobian[i][j]*jacobian[i][k]/SQR(error[i]);
          }
    }


    static int invert(double given[][], double result[][])
    // inverts a square real matrix using Gauss-Jordan method
    // caller supplies both given[][] and result[][]
    {
        for(int i=0; i<NADJ; i++)
          for(int j=0; j<NADJ; j++)
            result[i][j] = given[i][j];

        double det = gaussj(result, NADJ);
        if (Math.abs(det) < 1e-99)
          return -1;
        return 0; 
    }



    /////////////////////////////// main //////////////////////////

    public static void main (String args[]) throws IOException
    {

        setupspectra_Pickles(); 
        // setupspectra_Planck(); 

        prenormalize(); 

        setupNineFilters(); 

        //////////////// show the chosen parameters //////////////

        System.out.println(); 
        System.out.println("Parameters:");
        System.out.println(" SystematicSlope="+fwd(truep[0], 12, 6)); 
        for (int i=0; i<NFILTS; i++)
          System.out.println(" Filter"+fwi(i,2)+":  Thru="+fwd(truep[1+3*i],12,6)
             +"  Peak="+fwd(truep[2+3*i],12,6)+"  Width="+fwd(truep[3+3*i],12,6)); 


        ///////////////// show the raw photometry //////
 
        System.out.println(); 
        System.out.println("Integrated signal from each star, each band:"); 
        for (int istar=0; istar<NSTARS; istar++)
        {
            System.out.print(fws(picklenames_ten[istar],12)); 
            for (int jfilt=0; jfilt<NFILTS; jfilt++)
            {
                int idat = NSYS + istar + jfilt*NSTARS;  
                System.out.print(fwd(f(idat,truep),12,3));  // line 305
            }
            System.out.println(); 
        } 


        renormalize(); 

        /////// show the photometry normalized ////////
 
        System.out.println(); 
        System.out.println("Normalized signal from each star, each band:"); 
        for (int istar=0; istar<NSTARS; istar++)
        {
            System.out.print(fws(picklenames_ten[istar],12)); 
            for (int jfilt=0; jfilt<NFILTS; jfilt++)
            {
                int idat = NSYS + istar + jfilt*NSTARS;  
                System.out.print(fwe(f(idat,truep)));
            }
            System.out.println(); 
        } 



        ////////////////// create and show the Jacobian ///////////

        double[][] jacobian = new double[NPTS][NADJ];

        makejacobian(truep, jacobian); 

        System.out.println(); 
        System.out.println("Jacobian matrix at true parms:");
        System.out.println("                Sys"+
             "  Thru Peak Width   Thru Peak Width   Thru Peak Width"); 

        for (int i=0; i<NPTS; i++)
        {
            if (i==0)
              System.out.print(fws("SYSTEMATIC", 12)); 
            else
              System.out.print(fws(picklenames_ten[(i-1) % NSTARS], 12)); 
            for (int j=0; j<NADJ; j++)
              System.out.print(fwd(jacobian[i][j],8,0)); 
            System.out.println(); 
        }

        try 
        {
            fos = new FileOutputStream("JACOBIAN.TXT"); 
            pw = new PrintWriter(fos, true); 
            for (int i=0; i<NPTS; i++)
            {
                if (i==0)
                  pw.print(fws("SYSTEMATIC", 12)); 
                else
                  pw.print(fws(picklenames_ten[(i-1) % NSTARS], 12)); 
                for (int j=0; j<NADJ; j++)
                  pw.print(fwd(jacobian[i][j],8,0)); 
                pw.println(); 
            }
        }
        catch (Exception e) {}
        fos.close(); 


        ////////// Generate the Fisher matrix ////////////////////////

        for (int row=0; row<NPTS; row++)
          error[row] = SIGMA; 

        double[][] fisher = new double[NADJ][NADJ];

        getfisher(jacobian, error, fisher); 

        System.out.println(); 
        System.out.println("Fisher matrix at true parms:");
        for (int i=0; i<NADJ; i++)
        {
           for (int j=0; j<NADJ; j++)
             System.out.print(fwd(fisher[i][j],8,0)); 
           System.out.println(""); 
        }

        try 
        {
            fos = new FileOutputStream("FISHER.TXT"); 
            pw = new PrintWriter(fos, true); 
            for (int i=0; i<NADJ; i++)
            {
                for (int j=0; j<NADJ; j++)
                  pw.print(fwd(fisher[i][j],8,0)); 
                pw.println(); 
            }
        }
        catch (Exception e) {}
        fos.close(); 


        ///////// Generate the covariance matrix ////////////////////

        double[][] covar  = new double[NADJ][NADJ];   
 
        invert(fisher, covar);

        System.out.println(); 
        System.out.println("Covariance matrix at true parms:");
        for (int i=0; i<NADJ; i++)
        {
           for (int j=0; j<NADJ; j++)
             System.out.print(fwe(covar[i][j]));         
           System.out.println(""); 
        }

        System.out.println(); 
        System.out.println("RMS absolute parameter errors = sqrt(cov[i,i])...");
        for (int i=0; i<NADJ; i++)
          System.out.print(fwd(Math.sqrt(covar[i][i]), 8,4)); 
        System.out.println(); 

        try 
        {
            fos = new FileOutputStream("COVAR.TXT"); 
            pw = new PrintWriter(fos, true); 
            for (int i=0; i<NADJ; i++)
            {
                for (int j=0; j<NADJ; j++)
                  pw.print(fwe(covar[i][j])); 
                pw.println(""); 
            }
            pw.println(""); 
            pw.println("RMS absolute parm errors = sqrt(covar[i,i]..."); 
            for (int i=0; i<NADJ; i++)
              pw.print(fwd(Math.sqrt(covar[i][i]), 8, 4)); 
            pw.println(); 
        }
        catch (Exception e) {}
        fos.close(); 



        /////////////////// shut down //////////////////////////////

        System.out.println(); 
        System.out.println("That's all folks");
    }



    /////////////////////////////////utilities ////////////////////

    static double SQR(double x) 
    {
        return x*x;
    }


    static double gauss()
    {
        double sum = 0.0; 
        for (int i=0; i<12; i++)
           sum += Math.random(); 
        return sum - 6.0;
    }


    static String fwi(int n, int w)
    // converts an int to a string with given width.
    {
        String s = Integer.toString(n); 
        while (s.length() < w)
           s = " " + s; 
        return s;
    }


    static String fwd(double x, int w, int d)
    // converts a double to a string with given width and decimals.
    {
        java.text.DecimalFormat df = new DecimalFormat();
        df.setMaximumFractionDigits(d); 
        df.setMinimumFractionDigits(d);
        df.setGroupingUsed(false);   
        String s = df.format(x); 
        while (s.length() < w)
          s = " " + s;  
        if (s.length() > w)
        {
            s = "";
            for (int i=0; i<w; i++)
              s = s + "-";
        }  
        return s; 
    }

    static String fwe(double x)
    // converts a double into an E notation string
    {
        java.text.NumberFormat ef = new DecimalFormat("0.0E0"); 
        String s = ef.format(x);    
        while (s.length() < 8)
          s = " " + s; 
        return s; 
    }

    static String fws(String s, int len)
    // enlarges a String so that it occupies "len" characters
    {
        while (s.length() < len)
          s = s + " ";
        return s; 
    }

    static double log10(double arg)
    {
        if (arg>0.0)
          return Math.log(arg)/2.302585092994046;
        return -999.9;
    }

    static double dex(double arg)
    {
        if (arg > 300.0)
           arg = 300.0; 
        if (arg < -300.0)
           arg = -300.0; 
        return Math.pow(10.0, arg); 
    }

    static void beep()
    {
        char c = 7; 
        String s = ""+c; 
        System.out.print(s); 
    }
}