/*
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 <math.h>
#include <float.h>
#include <iomanip.h>
#include <imsl.h>
#include <sumsqr.h>
#include "bin_td.h"

inline double fifdsign(double a, double b)
{
  return b < 0. ? -a: a;
}

double R = 8.31441;

int iii = 1;
double xtmp;
elem* subst1;
elem* subst2;
double Tmin = 300.;
double Tmax = 4000.;
double xmin = 1e-10;
double xmax = 1. - 1e-10;
cString log_name;
ofstream ofpd;


/*************************************************************
from DMIXG_RK.RED
ON FORT;


DmixG;

      ANS=-(((2.*X-1.)**3*C3+(2.*X-1.)**2*C2+(2.*X-1.)*C1+
     . C0)*(X-1.)*LOG(T)*T*X+(LOG(-(X-1.))*R*T+A0*X+B0*T*X)
     . *(X-1.)+(A1+B1*T)*(2.*X-1.)*(X-1.)*X+(A2+B2*T)*(2.*X
     . -1.)**2*(X-1.)*X+(A3+B3*T)*(2.*X-1.)**3*(X-1.)*X-LOG
     . (X)*R*T*X)

DmixS;

      ANS=((2.*X-1.)**3*C3+(2.*X-1.)**2*C2+(2.*X-1.)*C1+C0)
     . *(X-1.)*LOG(T)*X+(LOG(-(X-1.))*R+B0*X+C0*X)*(X-1.)+(
     . B1+C1)*(2.*X-1.)*(X-1.)*X+(B2+C2)*(2.*X-1.)**2*(X-1.
     . )*X+(B3+C3)*(2.*X-1.)**3*(X-1.)*X-LOG(X)*R*X

DmixH;

      ANS=-((A0-C0*T)+(A1-C1*T)*(2.*X-1.)+(A2-C2*T)*(2.*X-
     . 1.)**2+(A3-C3*T)*(2.*X-1.)**3)*(X-1.)*X

Dmixg1;

      ANS=((8.*X-7.)*(2.*X-1.)**2*C3+(6.*X-5.)*(2.*X-1.)*C2
     . +(4.*X-3.)*C1+C0)*LOG(T)*T*X**2+(A1+B1*T)*(4.*X-3.)*
     . X**2+(A2+B2*T)*(6.*X-5.)*(2.*X-1.)*X**2+(A3+B3*T)*(
     . 8.*X-7.)*(2.*X-1.)**2*X**2+LOG(-(X-1.))*R*T+A0*X**2+
     . B0*T*X**2

Dmixg2;

      ANS=((8.*X-1.)*(2.*X-1.)**2*C3+(6.*X-1.)*(2.*X-1.)*C2
     . +(4.*X-1.)*C1+C0)*(X-1.)**2*LOG(T)*T+(A0+B0*T)*(X-1.
     . )**2+(A1+B1*T)*(4.*X-1.)*(X-1.)**2+(A2+B2*T)*(6.*X-
     . 1.)*(2.*X-1.)*(X-1.)**2+(A3+B3*T)*(8.*X-1.)*(2.*X-1.)
     . **2*(X-1.)**2+LOG(X)*R*T

Dmixs1;

      ANS=-(((8.*X-7.)*(2.*X-1.)**2*C3+(6.*X-5.)*(2.*X-1.)*
     . C2+(4.*X-3.)*C1+C0)*LOG(T)*X**2+(B1+C1)*(4.*X-3.)*X
     . **2+(B2+C2)*(6.*X-5.)*(2.*X-1.)*X**2+(B3+C3)*(8.*X-
     . 7.)*(2.*X-1.)**2*X**2+LOG(-(X-1.))*R+B0*X**2+C0*X**2)

Dmixs2;

      ANS=-(((8.*X-1.)*(2.*X-1.)**2*C3+(6.*X-1.)*(2.*X-1.)*
     . C2+(4.*X-1.)*C1+C0)*(X-1.)**2*LOG(T)+(B0+C0)*(X-1.)
     . **2+(B1+C1)*(4.*X-1.)*(X-1.)**2+(B2+C2)*(6.*X-1.)*(
     . 2.*X-1.)*(X-1.)**2+(B3+C3)*(8.*X-1.)*(2.*X-1.)**2*(X-
     . 1.)**2+LOG(X)*R)

Dmixh1;

      ANS=((A1-C1*T)*(4.*X-3.)+(A2-C2*T)*(6.*X-5.)*(2.*X-1.
     . )+(A3-C3*T)*(8.*X-7.)*(2.*X-1.)**2+A0-C0*T)*X**2

Dmixh2;

      ANS=((A0-C0*T)+(A1-C1*T)*(4.*X-1.)+(A2-C2*T)*(6.*X-1.
     . )*(2.*X-1.)+(A3-C3*T)*(8.*X-1.)*(2.*X-1.)**2)*(X-1.)
     . **2

Dmixconv;

      ANS=-(2.*((40.*X**2-40.*X+7.)*(2.*X-1.)*C3+(24.*X**2-
     . 24.*X+5.)*C2+3.*(2.*X-1.)*C1+C0)*(X-1.)*LOG(T)*T*X+2.
     . *(A0+B0*T)*(X-1.)*X+6.*(A1+B1*T)*(2.*X-1.)*(X-1.)*X+
     . 2.*(A2+B2*T)*(24.*X**2-24.*X+5.)*(X-1.)*X+2.*(A3+B3*T
     . )*(40.*X**2-40.*X+7.)*(2.*X-1.)*(X-1.)*X+R*T)/((X-1.
     . )*X)
****************************************************************/

double DmixG(double& T, double& x)
{
  return (x <= 0. || x >= 1.) ? 0. :
  -((pow(2.*x-1.,3)*C3+pow(2.*x-1.,2)*C2+(2.*x-1.)*C1+C0)*(x-1.)*log(T)*T*x
  +(log(-(x-1.))*R*T+A0*x+B0*T*x)*(x-1.)+(A1+B1*T)*(2.*x-1.)*(x-1.)*x
  +(A2+B2*T)*pow(2.*x-1.,2)*(x-1.)*x
  +(A3+B3*T)*pow(2.*x-1.,3)*(x-1.)*x-log(x)*R*T*x);
}

double DmixS(double& T, double& x)
{
  return
  (pow(2.*x-1.,3)*C3+pow(2.*x-1.,2)*C2+(2.*x-1.)*C1+C0)*(x-1.)*log(T)*x
  +(log(-(x-1.))*R+B0*x+C0*x)*(x-1.)+(B1+C1)*(2.*x-1.)*(x-1.)*x
  +(B2+C2)*pow(2.*x-1.,2)*(x-1.)*x+(B3+C3)*pow(2.*x-1.,3)*(x-1.)*x-log(x)*R*x;
}

double DmixH(double& T, double& x)
{
  return
  -((A0-C0*T)+(A1-C1*T)*(2.*x-1.)+(A2-C2*T)*pow(2.*x-1.,2)
  +(A3-C3*T)*pow(2.*x-1.,3))*(x-1.)*x;
}

double Dmixg1(double& T, double& x)
{
  return
  ((8.*x-7.)*pow(2.*x-1.,2)*C3+(6.*x-5.)*(2.*x-1.)*C2
  +(4.*x-3.)*C1+C0)*log(T)*T*pow(x,2)+(A1+B1*T)*(4.*x-3.)*pow(x,2)
  +(A2+B2*T)*(6.*x-5.)*(2.*x-1.)*pow(x,2)
  +(A3+B3*T)*(8.*x-7.)*pow(2.*x-1.,2)*pow(x,2)+log(-(x-1.))*R*T
  +A0*pow(x,2)+B0*T*pow(x,2);
}

double Dmixg2(double& T, double& x)
{
  return
  ((8.*x-1.)*pow(2.*x-1.,2)*C3+(6.*x-1.)*(2.*x-1.)*C2+(4.*x-1.)*C1+C0)
  *pow(x-1.,2)*log(T)*T+(A0+B0*T)*pow(x-1.,2)+(A1+B1*T)*(4.*x-1.)*pow(x-1.,2)
  +(A2+B2*T)*(6.*x-1.)*(2.*x-1.)*pow(x-1.,2)
  +(A3+B3*T)*(8.*x-1.)*pow(2.*x-1.,2)*pow(x-1.,2)+log(x)*R*T;
}

double Dmixs1(double& T, double& x)
{
  return
  -(((8.*x-7.)*pow(2.*x-1.,2)*C3+(6.*x-5.)*(2.*x-1.)*C2
  +(4.*x-3.)*C1+C0)*log(T)*pow(x,2)+(B1+C1)*(4.*x-3.)*pow(x,2)
  +(B2+C2)*(6.*x-5.)*(2.*x-1.)*pow(x,2)
  +(B3+C3)*(8.*x-7.)*pow(2.*x-1.,2)*pow(x,2)+log(-(x-1.))*R
  +B0*pow(x,2)+C0*pow(x,2));
}

double Dmixs2(double& T, double& x)
{
  return
  -(((8.*x-1.)*pow(2.*x-1.,2)*C3+(6.*x-1.)*(2.*x-1.)*C2
  +(4.*x-1.)*C1+C0)*pow(x-1.,2)*log(T)+(B0+C0)*pow(x-1.,2)
  +(B1+C1)*(4.*x-1.)*pow(x-1.,2)+(B2+C2)*(6.*x-1.)*(2.*x-1.)*pow(x-1.,2)
  +(B3+C3)*(8.*x-1.)*pow(2.*x-1.,2)*pow(x-1.,2)+log(x)*R);
}

double Dmixh1(double& T, double& x)
{
  return
  ((A1-C1*T)*(4.*x-3.)+(A2-C2*T)*(6.*x-5.)*(2.*x-1.)
  +(A3-C3*T)*(8.*x-7.)*pow(2.*x-1.,2)+A0-C0*T)*pow(x,2);
}

double Dmixh2(double& T, double& x)
{
  return
  ((A0-C0*T)+(A1-C1*T)*(4.*x-1.)+(A2-C2*T)*(6.*x-1.)*(2.*x-1.)
  +(A3-C3*T)*(8.*x-1.)*pow(2.*x-1.,2))*pow(x-1.,2);
}

double Dmixconv(double& T, double& x)
{
  return
  -(2.*((40.*pow(x,2)-40.*x+7.)*(2.*x-1.)*C3+(24.*pow(x,2)-24.*x+5.)*C2
  +3.*(2.*x-1.)*C1+C0)*(x-1.)*log(T)*T*x+2.*(A0+B0*T)*(x-1.)*x
  +6.*(A1+B1*T)*(2.*x-1.)*(x-1.)*x
  +2.*(A2+B2*T)*(24.*pow(x,2)-24.*x+5.)*(x-1.)*x
  +2.*(A3+B3*T)*(40.*pow(x,2)-40.*x+7.)*(2.*x-1.)*(x-1.)*x+R*T)/((x-1.)*x);
}

double elem::h(double& T)
{
  if (xs == 0.)
    return h1s(T);
  else if (xs == 1.)
    return h2s(T);
  else
  {
    of << "something wrong in elem::h" << endl;
    exit(1);
    return 0.;
  }
}

double elem::dif(double T, double x)
{
  if (xs == 0.)
    return g1l(T, x) - g1s(T);
  else if (xs == 1.)
    return g2l(T, x) - g2s(T);
  else
  {
    of << "something wrong in elem::dif" << endl;
    exit(1);
    return 0.;
  }
}

double linephase::dif(double T, double x)
{
  return (1. - xs)*g1l(T, x) + xs*g2l(T, x) - g(T);
}

void init_est()
{
  for (int i = 0; i < nnon; i++)
    non[i].est = false;
}

double difT(double *T)
{
  double dif = subst1->dif(*T);
    if (print == all_iter)
    of << subst1->name << "_L_solve "
       << setw(wdth) << T
       << setw(wdth) << dif << endl;
  return dif;
}

double T_solve(double *reg)
{
  double& T = reg[0];
  double& x = reg[1];
  double eps = 0., eps2 = 0., eta = 0.;
  double Tc = T;
  xtmp = x;
  int nsig = 10, n = 1, itmax = 70, ier;
  zreal2_(difT, &eps, &eps2, &eta, &nsig, &n, &Tc, &itmax, &ier);
  if (ier)
  {
    nsig = 8;
    itmax = 70;
    Tc = T;
    zreal1_(difT, &eps, &eps2, &eta, &nsig, &n, &Tc, &itmax, &ier);
    if (ier)
    {
//      if (print == all_iter)
      of << "WARNING - can't find solution for " << subst1->name << "_L"
         << endl;
      double delta = 10*sqrt(DBL_EPSILON);
      delta = delta*(T + delta);
			double x = T - delta;
      double f1 = difT(&x);
			x = T + delta;
      double f2 = difT(&x);
  //  y = y1 + (y2 - y1)/(x2 - x1)*(x - x1)
  //  y = 0  ->  x = x1 - y1*(x2 - x1)/(y2 - y1)
      return (f1 + f2)*delta/(f2 - f1);
    }
  }
  return T - Tc;
}

void xtoz(double x, double& z, double xmin, double xmax)
{
  z = log((x - xmin)/(xmax - x));
}

void ztox(double& x, double z, double xmin, double xmax)
{
  x = (xmax - xmin)/(1 + exp(-z)) + xmin;
}

void difTx(double* x, int *m, int *n, double* f)
{
  double T;
  double xc;
  ztox(T, x[0], Tmin, Tmax);
  ztox(xc, x[1], xmin, xmax);
  f[0] = subst1->dif(T, xc);
  f[1] = subst2->dif(T, xc);
  if (print == all_iter)
  {
    of << subst1->name << "_L_" << subst2->name << "_solve - ";
    of << setw(wdth) << T
       << setw(wdth) << xc
       << setw(wdth) << sqrt(pow(f[0], 2) + pow(f[1], 2))
       << endl;
  }
}

double Tx_solve(double* reg, int i)
{
  double& T = reg[0];
  if (!non[i].est)
  {
    double x[2];
    double x1[2];
    x1[0] = non[i].Tin;
    x1[1] = non[i].xin;
    xtoz(x1[0], x[0], Tmin, Tmax);
    xtoz(x1[1], x[1], xmin, xmax);
    int n = 2;
    int m = n;
    int nsig = 10;
    double eps = 0;
    double delta = 0;
    int maxfn = 300;
    int iopt = 1;
    int ixjac;
    int ier;
    int infer;
    double parm[4];
    double ssq;
    ixjac = m;
    double* f = new double[m];
    double* xjac = new double[n*m];
    double* xjtj = new double[n*(n+1)/2];
    double* work = new double[5*n + 2*m + (n+1)*n/2];
    double sig = 15;
    zxssq1_(difTx, &m, &n, &nsig, &eps, &delta, &maxfn, &iopt,
          parm, x, &ssq, f, xjac, &ixjac, xjtj, work, &infer, &ier, &sig);
    difTx(x, &m, &n, f);
    delete f;
    delete xjac;
    delete xjtj;
    delete work;
    ztox(non[i].T, x[0], Tmin, Tmax);
    ztox(non[i].x, x[1], xmin, xmax);
    if (ssq > 1e-4 || fabs(non[i].x - x1[1]) > 0.08)
    {
  //    if (print == all_iter)
//      of << "WARNING - can't find solution for " << subst1->name << "_L_"
//         << subst2->name << endl;
//      print_zxssq_1(infer, ier, work);
      non[i].conv = false;
      non[i].T1 = x1[0] - T_solve(x1);
      subst1 = subst2;
      non[i].T2 = x1[0] - T_solve(x1);
      non[i].T = (non[i].T1 + non[i].T2)/2.;
      non[i].x = x1[1];
    }
    else
      non[i].conv = true;
    non[i].est = true;
  }
  if (non[i].conv)
    return T - non[i].T;
  else
  {
    double del1 = T - non[i].T1;
    double del2 = T - non[i].T2;
    if (fifdsign(1, del1)* fifdsign(1, del2) > 0.)
      return (del1 + del2)*5.;
    else
      return (del2 - del1)*5.;
  }
}

/* old solution
    double eps = 0.;
    double fnorm;
    int nsig = 1, n = 2, itmax = 10, ier;
    double *wa, *par;
//    wa = new double[(n+1)*(3*n+8)];
//    zscnt(difTx1, nsig, n, itmax, par, x, fnorm, wa, ier);
//    delete wa;
    itmax = 70;
    nsig = 10;
    wa = new double[(n+2)*(n-1)/2 + 3*n];
    zsystm(difTx, eps, nsig, n, x, itmax, wa, par, ier);
    delete wa;
end old solution */