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

Software for paper 
E. B. Rudnyi. Statistical model of systematic errors: An assessment of the Ba-Cu and Cu-Y phase diagram. Chemometrics and Intelligent Laboratory Systems, 1997, v. 36, p. 213-227.

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

#include <fstream.h>
#include <math.h>
#include <float.h>
#include <iomanip.h>
#include <sumsqr.h>
#include <vector.h>
#include <map.h>
#include "bin_td.h"

int irun = 0;

void set_hide()
{
  iser = ser.begin();
  while(iser != ser.end())
  {
    series& t = *iser;
    double sr2, gb, ga;
    if (t.hide)
    {
      set_sums_1(t);
      if (t.Ns() > 1)
      {
        if (t.Ps())
        {
          if (t.Ns() == 2)
          {
            sr2 = t.sumi/2.;
            gb = 0.;
            ga = 0.;
          }
          else
          {
            sr2 = (t.sumi - t.suma*t.suma/t.Ns() - t.sumb*t.sumb/t.Ps())
                               /(t.Ns() - 2);
            gb = (t.sumb*t.sumb/t.Ps()/t.sr2 - 1)/t.Ps();
            ga = (t.suma*t.suma/t.Ns()/t.sr2 - 1)/t.Ns();
          }
        }
        else
        {
          sr2 = (t.sumi - t.suma*t.suma/t.Ns())/(t.Ns() - 1);
          gb = 0.;
          ga = (t.suma*t.suma/t.Ns()/t.sr2 - 1)/t.Ns();
        }
      }
      else
      {
        sr2 = t.sumi;
        ga = 0.;
        gb = 0.;
      }
      if (ga < 0.) t.ga = 0.;
      if (gb < 0.) t.gb = 0.;
      if (t.fl_sr2 == series::own)
        t.sr2 = sr2;
      else if (t.fl_sr2 == series::same)
        t.sr2 = sr2t[t.isr2];
      if (t.fl_ga == series::own)
        t.ga = ga;
      else if (t.fl_ga == series::same)
        t.ga = gat[t.igat];
      if (t.fl_gb == series::own)
        t.gb = gb;
      else if (t.fl_ga == series::same)
        t.gb = gbt[t.igbt];
    }
    iser++;
  }
}

class err
{
public:
  cString ID;
  double err_a;
  double err_b;
  err() : ID(10), err_a(0.), err_b(0.) {}
};

class vec_err : public vector <err>
{
public:
  double sr;
};

vector <vec_err> errors;

void deviates()
{
//  ofstream of1((file_name + "ax1").c_str(), ios::out | ios::app);
  ofstream of1((file_name + "ax1").c_str());
  of1.precision(prcsn);
  of1.setf(ios::showpoint);
  iser = ser.begin();
  while(iser != ser.end())
  {
    series& t = *iser;
    of1 << setw(11 - t.ID.stringlen() + 2) << 'x' << t.ID.substr(2, 100);
    of1 << setw(11 - t.ID.stringlen() + 4) << t.ID.substr(2, 100);
    iser++;
  }
  of1 << endl;
  int end = 0;
  int ii, jj;
  for (int i = 0; !end; ++i)
  {
    end = 1;
    iser = ser.begin();
    while(iser != ser.end())
    {
      series& t = *iser;
      if (i < t.Ntot())
      {
        end = 0;
        if (t.Nvar() > 1)
          of1 << setw(11) << t.val(i, 1);
        else
          of1 << setw(11) << patch_x(t.ID, i);
//        of1 << setw(11) << (*t.f)(t(i))/sqrt(t.sr2);
//        of1 << setw(11) << (*t.f)(t(i));
        for (ii = 0; ii < errors.size(); ++ii)
          for (jj = 0; jj < errors[ii].size(); ++jj)
            if (errors[ii][jj].ID == t.ID)
              goto found;
        cerr << "deviates - series not found" << endl;
        exit(1);
      found:
        of1 << setw(11) << (*t.f)(t(i))/errors[ii].sr;
      }
      else
      {
        of1 << setw(11) << "miss";
        of1 << setw(11) << "miss";
      }
      iser++;
    }
    of1 << endl;
  }
}

void read_errors(char* name)
{
  errors.erase(errors.begin(), errors.end());
  ifstream in(name);
  err t;
  char buffer[100];
  int i = 0;
  while(in)
  {
    in.getline(buffer, 100);
    errors.push_back(vec_err());
    errors[i].sr = atof(buffer);
    while (1)
    {
      in.getline(buffer, 100);
      if (!strlen(buffer))
        break;
      t.ID = buffer;
      errors[i].push_back(t);
    }
    i++;
  }
}

void eb_vs_ea()
{
  ofstream of2((file_name + "ax2").c_str());
  of2.precision(prcsn);
  of2.setf(ios::showpoint);
  int i, j;
  for (i = 0; i < errors.size(); ++i)
    of2 << " ID" << i + 1 << " err_a" << i + 1 << " err_b" << i + 1;
  of2 << endl;
  iser = ser.begin();
  while(iser != ser.end())
  {
    series& t = *iser;
    for (i = 0; i < errors.size(); ++i)
      for (j = 0; j < errors[i].size(); ++j)
        if (errors[i][j].ID == t.ID)
          goto found;
    iser++;
    continue;
  found:
    if (t.hide)
      errors[i][j].ID = "*";
    else
      errors[i][j].ID = "";
    errors[i][j].ID += t.ID.substr(2, 100);
    errors[i][j].err_a = t.suma/t.Ns()/errors[i].sr;
    if (t.Ps())
      errors[i][j].err_b =
            t.sumb/t.Ps()*sqrt(t.Ps()/t.Ns())/errors[i].sr;
    else
      errors[i][j].err_b = 0.;
    iser++;
  }
  int end = 0;
  for (j = 0; !end; ++j)
  {
    end = 1;
    for (i = 0; i < errors.size(); ++i)
      if (j < errors[i].size())
      {
        end = 0;
        of2 << setw(wdth) << errors[i][j].ID
            << setw(wdth) << errors[i][j].err_a
            << setw(wdth) << errors[i][j].err_b;
      }
      else
        of2 << setw(wdth) << "miss"
            << setw(wdth) << "miss"
            << setw(wdth) << "miss";
   of2 << endl;
  }
}

typedef map<cString, vector<double>, less<cString> > maptype;
maptype extra_x;

void read_extra_x(char *name)
{
  extra_x.erase(extra_x.begin(), extra_x.end());
  ifstream in(name);
  cString ID;
  vector<double> x;
  double value;
  char buffer[100];
  while(in)
  {
    ID = "";
    x.erase(x.begin(), x.end());
    in.getline(buffer, 100);
    ID = buffer;
    while (1)
    {
      in.getline(buffer, 100);
      if (!strlen(buffer))
        break;
      x.push_back(atof(buffer));
    }
    extra_x.insert(pair<const cString, vector<double> >(ID, x));
  }
}

double patch_x(cString &ID, size_t i)
{
  vector<double> &x = extra_x[ID];
  if (i > x.size())
    return 0.;
  else
    return x[i];
}

void test_convex_L()
{
  of << "testing Gibbs energy of the melt" << endl;
  for (double T = 300.; T <= 3000.01; T+=50.)
  {
    int okay = 1;
    for (double x = 0.02; x <= 0.98; x+= 0.02)
      if (Dmixconv(T, x) <= 0.)
      {
        if (okay)
        {
          of << setw(wdth) << T
             << setw(wdth) << x;
          okay = 0;
        }
        else if (x > 0.97)
          of << setw(wdth) << x << endl;
      }
      else
      {
        if (!okay)
        {
          of << setw(wdth) << x << endl;
          okay = 1;
        }
      }
  }
  of << endl;
}

void test_convex_solids()
{
  of << "testing Gibbs energies of the solids" << endl;
  double Gl, Gm, Gr;
  double xl, xm, xr;
  int i, j, k;
  for (double T = 100.; T <= 2000.01; T+=50.)
  {
    vector <linephase1> lp1(nlp);
    for (i = 0; i < nlp; ++i)
    {
      lp1[i].name = lp[i].name;
      lp1[i].xs = lp[i].xs;
      lp1[i].H = lp[i].H;
      lp1[i].S = lp[i].S;
      lp1[i].C = lp[i].C;
    }
    int okay = 1;
    for (i = -1, j = 0, k= 1; k <= lp1.size(); )
    {
      if (i == -1)
      {
        Gl = 0;
        xl = e1.xs;
      }
      else
      {
        Gl = lp1[i].G(T);
        xl = lp1[i].xs;
      }
      if (k == lp1.size())
      {
        Gr = 0;
        xr = e2.xs;
      }
      else
      {
        Gr = lp1[k].G(T);
        xr = lp1[k].xs;
      }
      Gm = lp1[j].G(T);
      xm = lp1[j].xs;
      if (Gm > (Gr - Gl)/(xr - xl)*(xm - xl) + Gl)
      {
        okay = 0;
        of << lp1[j].name << " unstable "
           << setw(wdth) << T
           << setw(wdth) << xl
           << setw(wdth) << xm
           << setw(wdth) << xr
           << setw(wdth) << Gl
           << setw(wdth) << Gm
           << setw(wdth) << Gr
           << endl;
        lp1.erase(lp1.begin() + j);
        if (i > -1)
        {
          --i;
          --j;
          --k;
        }
      }
      else
      {
        if (!okay)
          of << lp1[j].name << " stable "
             << setw(wdth) << T
             << setw(wdth) << xl
             << setw(wdth) << xm
             << setw(wdth) << xr
             << setw(wdth) << Gl
             << setw(wdth) << Gm
             << setw(wdth) << Gr
             << endl;
        ++i;
        ++j;
        ++k;
      }
    }
  }
}

/*
int main()
{
  memory test;
  read_errors("bacu_e.dat");
  for (int i = 0; i < errors.size(); ++i)
  {
    cout << errors[i].sr << " ";
    for (int j = 0; j < errors[i].size(); ++j)
      cout << " " << errors[i][j].ID;
    cout << endl;
  }
  errors.erase(errors.begin(), errors.end());
  return 0;
}
*/
/*
int main()
{
  memory test;
  maptype::iterator i;
  read_extra_x("bacu_x.dat");
  for (i = extra_x.begin(); i != extra_x.end(); ++i)
  {
    cout << (*i).first << " ";
    for (int j = 0; j < (*i).second.size(); ++j)
      cout << (*i).second[j] << " ";
    cout << endl;
  }
  extra_x.erase(extra_x.begin(), extra_x.end());
  return 0;
}
*/