/*
Copyright (C) 1995 Evgenii Rudnyi, http://Evgenii.Rudnyi.Ru

Software for paper 
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 <sys/time.h>
#include <math.h>
#include <float.h>
#include <iomanip.h>
#include <fstream.h>
#include <strstream>
#include <string.h>
#include <stdlib.h>
#include <string>
#include "imsl.h"

using namespace::std;

static float** dev = NULL;
static float** x = NULL;
static float* deva = NULL;
static float* devb = NULL;
static int* Ni = NULL;
static float* Pi = NULL;
static float* xav = NULL;
static int N = 0;
static char key[10];

void clear()
{
	int i;
  for (i = 0; i < N; i++)
    delete [] dev[i];
  delete [] dev;
  delete [] deva;
  delete [] Ni;
  if (strcmp(key, "line") == 0)
  {
    for (i = 0; i < N; i++)
      delete [] x[i];
    delete [] x;
    delete [] devb;
    delete [] Pi;
    delete [] xav;
  }
}

void my_handler()
{
  cout << "not enough memory" << endl;
  clear();
  exit(1);
}

int main(int argc, char* argv[])
{
  double seed;
  int Ni_max = 0;
  float si = 1.;
  float sa = 0.;
  float sb = 0.;
  float a = 0.;
  float b = 0.;
  const int buf_sz = 1024*5;
  char oneline[buf_sz];
  unsigned int i, j;
  ifstream seedf("generate.rnd", ios::in);
  ifstream setf;
  ofstream linearf;
  string filename;

  set_new_handler(my_handler);
  if (argc < 2)
  {
    cout <<
"Evgenii Rudnyi 1994, (c) All rights reserved \n\
version 1.05 (28 June 1994) \n\n\
to generete normal distributed random values to check the program LINEAR \n\
see file README for more information \n\n\
    program setfile [outfile sga sgb] \n\n\
see files ONEWAY.CFG and LINE.CFG as examples of the setfile \n\
see discussion of sga and sgb values in the file README \n";
    return 1;
  }
  filename = argv[1];
  if (filename.find(".") == string::npos)
    filename += ".cfg";
  setf.open(filename.c_str(), ios::in);
  if (!setf)
  {
    cout << "file " << filename << " does not exist" << endl;
    return 1;
  }
  if (!seedf)
  {
    cout << "random value generated by gettime" << endl;
//    seed = biostime(0, 0);
		timeval tim;
// This is working on my Mac. You may need change it.		
		gettimeofday(&tim, (struct timezone*)0);
		seed = double(tim.tv_sec) + tim.tv_usec*1e-6;
    ofstream seedf("generate.rnd");
 		if (seed < 1.)
			seed = 1.;
		if (seed > 2147483647.)
			seed = 2147483640.;
		seedf << seed << endl; 
  }
  else
  {
    seedf >> seed;
    seedf.close();
  }
  if (argc > 2)
  {
    filename = argv[2];
    if (filename.find(".") == string::npos)
    {
      linearf.open((filename + ".dat").c_str());
    }
    else
    {
      linearf.open((filename.substr(0, filename.find_last_of(".") + 1) + "dat").c_str());
    }
  }
  linearf.precision(3);
  if (argc > 3)
    sa = atof(argv[3]);
  if (argc > 4)
    sb = atof(argv[4]);

  setf.getline(oneline, buf_sz);
  {
    istrstream in(oneline);
    in.width(10);
    in >> key;
  }
  if (strcmp(key, "oneway") == 0)
  {
    setf.getline(oneline, buf_sz);
    a = atof(oneline);
    setf.getline(oneline, buf_sz);
    si = atof(oneline);
    sa *= si;
    while (setf)
    {
      setf.getline(oneline, buf_sz);
      i = atoi(oneline);
      if (i)
      {
        if (i > Ni_max) Ni_max = i;
        N++;
      }
      else
        break;
    }
    if (!N)
    {
      cout << "number of series is zero - nothing to do" << endl;
      return 1;
    }
    dev = new float*[N];
    deva = new float[N];
    Ni = new int[N];
    setf.clear();
    setf.seekg(0);
    setf.getline(oneline, buf_sz);
    setf.getline(oneline, buf_sz);
    setf.getline(oneline, buf_sz);
    for (i = 0; i < N; i++)
    {
      setf.getline(oneline, buf_sz);
      Ni[i] = atoi(oneline);
    }
    ggnml_(&seed, &N, deva);
    cout << "systematic errors" << endl;
    for (i = 0; i < N; i++)
    {
      deva[i] *= sa;
      cout << setw(5) << int(deva[i]);
    }
    cout << endl << endl;
    for (i = 0; i < N; i++)
    {
      cout << "series " << setw(2) << i + 1 << endl;
      linearf << "series_" << i + 1 << ", justa, y,";
      dev[i] = new float[Ni[i]];
      ggnml_(&seed, &Ni[i], dev[i]);
      for (j = 0; j < Ni[i]; j++)
      {
        dev[i][j] = a + dev[i][j]*si + deva[i];
        cout << setw(5) << (int)dev[i][j];
        if (j % 5 == 0)
          linearf << endl;
        linearf << dev[i][j];
        if (j == Ni[i] - 1)
          linearf << ';';
        else
          linearf << ",  ";
      }
      cout << endl;
      linearf << endl;
    }
  }
  else if (strcmp(key, "line") == 0)
  {
    setf.getline(oneline, buf_sz);
    a = atof(oneline);
    setf.getline(oneline, buf_sz);
    b = atof(oneline);
    setf.getline(oneline, buf_sz);
    si = atof(oneline);
    sa *= si;
    sb *= si;
    while (setf)
    {
      setf.getline(oneline, buf_sz);
      istrstream in(oneline);
      float tmp;
      int i = -1;
      do
      {
        i++;
        tmp = HUGE_VAL;
        in >> tmp;
      }
      while (tmp != HUGE_VAL);
      if (i)
      {
        if (i > Ni_max) Ni_max = i;
        N++;
      }
      else
        break;
    }
    if (!N)
    {
      cout << "number of series is zero - nothing to do" << endl;
      return 1;
    }
    dev = new float*[N];
    x = new float*[N];
    deva = new float[N];
    devb = new float[N];
    Ni = new int[N];
    Pi = new float[N];
    xav = new float[N];
    setf.clear();
    setf.seekg(0);
    setf.getline(oneline, buf_sz);
    setf.getline(oneline, buf_sz);
    setf.getline(oneline, buf_sz);
    setf.getline(oneline, buf_sz);
    for (i = 0; i < N; i++)
    {
      setf.getline(oneline, buf_sz);
      istrstream in(oneline);
      float tmp;
      int j = -1;
      do
      {
        j++;
        tmp = HUGE_VAL;
        in >> tmp;
      }
      while (tmp != HUGE_VAL);
      Ni[i] = j;
      Pi[i] = 0.;
      xav[i] = 0.;
      dev[i] = new float[j];
      x[i] = new float[j];
      in.clear();
      in.seekg(0);
      for (j = 0; j < Ni[i]; j++)
      {
        in >> tmp;
        x[i][j] = tmp;
        xav[i] += tmp;
        Pi[i] += tmp*tmp;
      }
      xav[i] /= Ni[i];
      Pi[i] -= Ni[i]*(xav[i]*xav[i]);
    }
    ggnml_(&seed, &N, deva);
    ggnml_(&seed, &N, devb);
    cout << "systematic errors" << endl;
    for (i = 0; i < N; i++)
    {
      deva[i] *= sa;
      cout << setw(5) << int(deva[i]);
    }
    cout << endl;
    for (i = 0; i < N; i++)
    {
      devb[i] *= sb;
      cout << setw(5) << int(devb[i]);
    }
    cout << endl << endl;
    for (i = 0; i < N; i++)
    {
      cout << "series " << setw(2) << i + 1 << endl;
      linearf << "series_" << i + 1 << ", line, y x,";
      ggnml_(&seed, &Ni[i], dev[i]);
      for (j = 0; j < Ni[i]; j++)
      {
        dev[i][j] = a + b*x[i][j] + dev[i][j]*si + deva[i]
                                  + devb[i]*(x[i][j] - xav[i]);
        cout << setw(5) << int(x[i][j]);
        if (j % 3 == 0)
          linearf << endl;
        linearf << dev[i][j] << "  " << x[i][j];
        if (j == Ni[i] - 1)
          linearf << ';';
        else
          linearf << ",  ";
      }
      cout << endl;
      linearf << endl;
      for (j = 0; j < Ni[i]; j++)
        cout << setw(5) << int(dev[i][j]);
      cout << endl;
    }
  }
  else
  {
    cout << "don't know how to handle the case - " << key << endl;
    return 1;
  }
  clear();
  {
    ofstream seedf("generate.rnd");
		seedf << seed << endl;
  }
  return 0;
}