/*
Copyright (C) 1994-1998 Evgenii Rudnyi, http://Evgenii.Rudnyi.Ru/

Implementation of the algorithm described in  
E.B. Rudnyi. Statistical model of systematic errors: linear error model. Chemometrics and Intelligent Laboratory Systems. 1996, V. 34, N 1, p. 41-54. 

This software is a copyrighted work licensed under the terms, described
in the file "FREE_LICENSE". 
*/

#include <new.h>
#include <iomanip.h>
#include <ctype.h>
#include <math.h>
#include <stdlib.h>
#include "sumsqr.h"
#include "imsl.h"

set<residual, less<residual> > res;
set<residual, less<residual> >::iterator ires;
list<series> ser;
list<series>::iterator iser;

fl_print print = results;
fl_show show = used;
fl_series print_series = no_series;
fl_start start = variances;
int iss = 0;
ofstream of;

double sr2t[nvar];
int Nst[nvar];
double gat[nvar];
int Ngat[nvar];
double gbt[nvar];
int Ngbt[nvar];

void (*old_handler)();

void my_handler()
{
  warning("memory is out - you are going to have problems");
  set_new_handler(old_handler);
}

static const char* NOSER = "\nNothing to do - there is no series \n";

void free_mem()
{
}

int main(int argc, char *argv[])
{
  try
  {
  //  memory test;
    ifstream df, sf;
    old_handler = set_new_handler(my_handler);
    if (!open_streams(argc, argv, df, sf))
      return 1;

    for (int ii = 0; ii < nvar; ++ii)
    {
      sr2t[ii] = 1.;
      Nst[ii] = 0;
      gat[ii] = 0.;
      Ngat[ii] = 0;
      gbt[ii] = 0.;
      Ngbt[ii] = 0;
    }
    double zxssq_machine_eps = 12;
    init_res();
    read_series(df);
    init_par();
    if (ser.empty())
    {
      cout << NOSER;
      free_mem();
      return 1;
    }

  // parameters for ZXSSQ
    int n;
    int m;
    int i, j;
    int nsig = 5;
    double eps = 1e-7;
    double delta = 0;
    int maxfn = 100;
    int iopt = 1;
    int ixjac;
    int ier;
    int infer;
    double parm[4];
    double ssq;
//    int level = 0, levold;
//    uerset(level,levold);
//    ugetio(3, &df, &of);

  // parameters for main loop
    int iter;
    int max_iter = 10;
    double eps_L = 1e-5;
    double eps_v = 1e-5;
    double difL;
    double difv;

  // other stuff
    double SS;
    double L = 0.;
    double Lold = 0.;

  // start read setfile
    do
    {
      init_est();
      difL = eps_L + 1.;
      difv = eps_v + 1.;
      of << endl << endl << "-------------Start of calculation-------------"
         << endl;
      read_setfile(sf);
      anal_series();
      of << "initial values of parameters" << endl;
      for (j = 0; j < nparam; j++)
      {
        par[j].sx = 0.;
        print_par(par[j]);
      }
      of << endl;

    // calculations
      iter = 1;
      n = 0;
      for (j = 0; j < nparam; j++)
        if (par[j].atr) n++;
      m = 0;
      iser = ser.begin();
      while (iser != ser.end())
      {
        if (!(*iser).hide)
          m += (*iser).Ns();
        iser++;
      }
      if (m && n)
      {
        ixjac = m;
        double* x = new double[n];
        double* f = new double[m];
        double * xjac = new double [long(n)*long(m)];
        double* xjtj = new double[n*(n+1)/2];
        double* work = new double[5*n + 2*m + (n+1)*n/2];
        if (!x || !f || !xjac || !xjtj || !work)
        {
          cout << endl << "not enough memory to run ZXSSQ - aborted" << endl;
          free_mem();
          return 1;
        }
        i = 0;
        for (j = 0; j < nparam; j++)
          if (par[j].atr)
          {
            x[i] = par[j].x;
            i++;
          }
        if (start == variances)
        {
          init_est();
          iser = ser.begin();
          while(iser != ser.end())
          {
            series& t = *iser;
            if (!t.hide)
            set_sums_1(t);
            iser++;
          }
          cout << endl << "   var_comp       ";
          eval_var_comp();
          eval_L(L, SS);
          of << "value of L is " << setw(15) << L;
          of << "  and of SS is " << setw(15)
             << SS << endl << endl;
        }

        while ((iter < 3) || (iter <= max_iter && difL > eps_L && difv > eps_v))
        {
          iss = 0;
          cout << endl;
          cout << "iteration " << setw(4) << iter << "      ZXSSQ       ";
          if (print != results)
            of << endl << "***** iteration" << setw(4) << iter
               << "   *******" << endl;
          zxssq1_(ss, &m, &n, &nsig, &eps, &delta, &maxfn, &iopt,
              parm, x, &ssq, f, xjac, &ixjac, xjtj, work, &infer, &ier, 
							&zxssq_machine_eps);
          ss(x, &m, &n, f);
          if (print != results)
          {
            print_zxssq_1(infer, ier, work);
            print_ss_par_1();
            eval_L(L, SS);
            of << "value of L is " << setw(15)
               << L;
            of << "  and of SS is " << setw(15)
               << SS << endl << endl;
          }

          cout << "   var_comp       ";
          eval_var_comp();
          eval_L(L, SS);
          difL = fabs((L - Lold)/L);
          Lold = L;
          difv = 0.;
          iser = ser.begin();
          while(iser != ser.end())
          {
            series& t = *iser;
            if (!t.hide)
            {
              if (t.fl_sr2 != series::fixed)
                difv = max(difv, fabs((t.sr2 - t.sr2old)/t.sr2));
              if (t.fl_ga != series::fixed)
              {
                if (t.ga)
                  difv = max(difv, fabs((t.ga - t.gaold)/t.ga));
                else
                  difv = max(difv, fabs(t.ga - t.gaold));
              }
              if (t.fl_gb != series::fixed)
              {
                if (t.gb)
                  difv = max(difv, fabs((t.gb - t.gbold)/t.gb));
                else
                  difv = max(difv, fabs(t.gb - t.gbold));
              }
              t.sr2old = t.sr2;
              t.gaold = t.ga;
              t.gbold = t.gb;
            }
            iser++;
          }
          iter++;
        }
        of << endl << "Process ended after " << setw(4) << iter - 1
           << " iterations" << endl;
        of << endl << "Last ZXSSQ report" << endl;
        print_zxssq_1(infer, ier, work);
        of << endl;
        of << "difL " << setw(wdth) << difL;
        of << ", difv " << setw(wdth) << difv << endl;
        of << "value of L is " << setw(15) << L;
        of << "  and of SS is " << setw(15) << SS
           << endl << endl;
        cout << endl;
    // estimate confidence limits
        int ijob = 1;
        linv3p_(xjtj, work, &ijob, &n, &ier);
        if (ier)
        {
          warning("can't invert (X'*X) - skipped confidence intervals");
          of << "Final values of parameters" << endl;
          for (j = 0; j < nparam; j++)
            print_par(par[j]);
          of << endl;
        }
        else
        {
          of << endl << "Parameter +/- standard deviation" << endl;
          i = 1;   //Fortran i
          for (j = 0; j < nparam; j++)
          {
            if (!par[j].atr)
              print_par(par[j]);
            else
            {
              par[j].sx = sqrt(xjtj[i*(i+1)/2 - 1]);
              print_par(par[j], par[j].sx);
              i++;
            }
          }
          of << endl << "Correlation matrix - only for free parameters" << endl;
          of << setw(wdth) << 1 << endl;
          for (i = 2; i <= n; i++)      // Fortran i
          {
            for (j = 1 ; j <= i - 1; j++)  // Fortran j
              of << setw(wdth)
                 << xjtj[i*(i-1)/2+j-1]/sqrt(xjtj[i*(i+1)/2-1]*xjtj[j*(j+1)/2-1]);
              of << setw(wdth) << 1 << endl;
          }
        }
        delete x;
        delete f;
        delete xjac;
        delete xjtj;
        delete work;
      }
      else
      {
        warning("nothing to calculate - skipped main loop");
        of << "values of parameters" << endl;
        fl_show old = show;
        if (show == free_par) show = used;
        for (j = 0; j < nparam; j++)
          print_par(par[j]);
        show = old;
        of << endl;
        if (m)
        {
          init_est();
          iser = ser.begin();
          while(iser != ser.end())
          {
            series& t = *iser;
            if (!t.hide)
            set_sums_1(t);
            iser++;
          }
          cout << "   var_comp       ";
          eval_var_comp();
          cout << endl;
          eval_L(L, SS);
          of << "value of L is " << setw(15) << L;
          of << "  and of SS is " << setw(15)
             << SS << endl << endl;
        }
      }

    // print all the points
      if (print == results)
      {
        init_est();
        print_all_points();
        post_analysis();
      }
    // end calculations
    }
    while (sf >> ws, sf.peek() != EOF);
  // end read setfile

//    uerset(levold,levold);
    free_mem();
  }
//  catch (gError &t)
//  {
//    cout << t.message << endl;
//  }
  catch (...)
  {
    cout << endl << endl << "Unknown exception has been thrown" << endl;
  }
  return 0;
}