/* ---------------------------------------------------------------------------

Copyright (C) 1990, 1991, 1993, 1994, 2002.

A. Ronald Gallant
Post Office Box 659
Chapel Hill NC 27514-0659 
USA   

Permission to use, copy, modify, and distribute this software and its 
documentation for any purpose and without fee is hereby granted, provided 
that the above copyright notice appears in all copies and that both that 
copyright notice and this permission notice appear in supporting 
documentation.  

This software is provided "as is" without any expressed or implied warranty.

----------------------------------------------------------------------------*/

#include "model.h"

model::model() { }

model::~model() { }

void model::e(INTEGER t)
{
  REAL peak, inter, base;
  REAL t1, t2, t3, t4, t5, t6, t7, t8;
  REAL r1, r2, r3;
  REAL x1, x2, x3;
  REAL y1, y2, y3;

  t1 = s.theta[1];
  t2 = s.theta[2];
  t3 = s.theta[3];
  t4 = s.theta[4];
  t5 = s.theta[5];
  t6 = s.theta[6];
  t7 = s.theta[7];
  t8 = s.theta[8];
  
  y1 = data(1,t);
  y2 = data(2,t);
  y3 = data(3,t);

  r1 = data(4,t);
  r2 = data(5,t);
  r3 = data(6,t);

  x1 = r1 - y3;
  x2 = r2 - y3;
  x3 = r3 - y3;

  peak  = t1 + t2*x1 + t3*x2 + t4*x3;  
  inter = t5 + t3*x1 + t6*x2 + t7*x3;
  base  = -1 + t4*x1 + t7*x2 + t8*x3;

  e_t[1] = y1 - log(peak/base);
  e_t[2] = y2 - log(inter/base);
}

void model::dele(INTEGER t) 
{
  REAL peak, inter, base;
  REAL t1, t2, t3, t4, t5, t6, t7, t8;
  REAL r1, r2, r3;
  REAL x1, x2, x3;
  REAL y1, y2, y3;

  t1 = s.theta[1];
  t2 = s.theta[2];
  t3 = s.theta[3];
  t4 = s.theta[4];
  t5 = s.theta[5];
  t6 = s.theta[6];
  t7 = s.theta[7];
  t8 = s.theta[8];

  y1 = data(1,t);
  y2 = data(2,t);
  y3 = data(3,t);

  r1 = data(4,t);
  r2 = data(5,t);
  r3 = data(6,t);

  x1 = r1 - y3;
  x2 = r2 - y3;
  x3 = r3 - y3;

  peak  = t1 + t2*x1 + t3*x2 + t4*x3;  
  inter = t5 + t3*x1 + t6*x2 + t7*x3;
  base  = -1 + t4*x1 + t7*x2 + t8*x3;

//e_t[1] = y1 - log(peak/base);
//e_t[2] = y2 - log(inter/base);


  dele_t(1,1) = -    (1.0/peak);  
  dele_t(1,2) = - (1.0/peak)*x1;  
  dele_t(1,3) = - (1.0/peak)*x2;
  dele_t(1,4) = - (1.0/peak)*x3   + (1.0/base)*x1;  
  dele_t(1,5) =               0;
  dele_t(1,6) =               0;  
  dele_t(1,7) =                     (1.0/base)*x2;
  dele_t(1,8) =                     (1.0/base)*x3;

  dele_t(2,1) =                0;  
  dele_t(2,2) =                0;  
  dele_t(2,3) = - (1.0/inter)*x1;
  dele_t(2,4) =                     (1.0/base)*x1;  
  dele_t(2,5) = -    (1.0/inter);
  dele_t(2,6) = - (1.0/inter)*x2;  
  dele_t(2,7) = - (1.0/inter)*x3  + (1.0/base)*x2;
  dele_t(2,8) =                     (1.0/base)*x3;
}

int model::initialize()
{
  if (s.p != 8 || s.M != 2) {cerr << "Error, model::initialize.\n"; exit(1);}


  ifstream electa("electa.dat");
  if (!electa) {
    cerr << "Error, model::model, Cannot open electa.dat\n";
    exit(1);
  }

  ifstream electc1("electc1.dat");
  if (!electc1) {
    cerr << "Error, model::model, Cannot open electc1.dat\n";
    exit(1);
  }

  ifstream electc2("electc2.dat");
  if (!electc2) {
    cerr << "Error, model::model, Cannot open electc2.dat\n";
    exit(1);
  }

  double t1, treat, base, inter, peak, expend;
  double t2, famsize, income, sqfeet, heatlos, range, wash, dry, cac, wac;
  double t3, single, duplex, mobile, hwh, frez, ref;

  data.resize(20,224);

  int missing;

  s.n = 224;

  INTEGER t0 = 1;

  for (INTEGER t=1; t<=224; t++) {

    missing = 0;

    electa >>  t1;
    electa >>  treat;
    electa >>  base;
    electa >>  inter;
    electa >>  peak;
    electa >>  expend;

    electc1 >> t2;
    electc1 >> famsize; if(famsize == -9999.) missing = 1;  //used for d8 d10
    electc1 >> income;  if(income == -9999.) missing = 1;   //used for d0
    electc1 >> sqfeet;  if(sqfeet == -9999.) missing = 1;   //used for d3
    electc1 >> heatlos; if(heatlos == -9999.) missing = 1;  //used for d6
    electc1 >> range;   if(range == -9999.) missing = 1;    //used for d13
    electc1 >> wash;    if(wash == -9999.) missing = 1;     //used for d9
    electc1 >> dry;     if(dry == -9999.) missing = 1;      //used for d10
    electc1 >> cac;     if(cac == -9999.) missing = 1;      //used for d6
    electc1 >> wac;     if(wac == -9999.) missing = 1;      //used for d7

    electc2 >> t3;
    electc2 >> single; // if(single == -9999.) missing = 1; //not used
    electc2 >> duplex;  if(duplex == -9999.) missing = 1;   //used for d4
    electc2 >> mobile;  if(mobile == -9999.) missing = 1;   //used for d5
    electc2 >> hwh;     if(hwh == -9999.) missing = 1;      //used for d8 d9
    electc2 >> frez;    if(frez == -9999.) missing = 1;     //used for d12
    electc2 >> ref;     if(ref == -9999.) missing = 1;      //used for d11
   
    if (t1 != t2 || t1 != t3) {
       cerr << "Error, model::initialize, t, t1, t2, t3 = " 
            << t << "  " << t1 << "  " << t2 << "  " << t3 << ".\n"; 
       exit(1);
     }

    double p1, p2, p3;

    p1 = p2 = p3 = 0;

    if (treat == 1) {  p1=3.90; p2=2.86; p3=1.06; }
    if (treat == 2) {  p1=3.90; p2=2.86; p3=1.78; }
    if (treat == 3) {  p1=3.90; p2=3.90; p3=1.06; }
    if (treat == 4) {  p1=3.90; p2=3.90; p3=1.78; }
    if (treat == 5) {  p1=5.06; p2=3.34; p3=1.37; }
    if (treat == 6) {  p1=6.56; p2=2.86; p3=1.06; }
    if (treat == 7) {  p1=6.56; p2=2.86; p3=1.78; }
    if (treat == 8) {  p1=6.56; p2=3.90; p3=1.06; }
    if (treat == 9) {  p1=6.56; p2=3.90; p3=1.78; }


    double y1  = log(peak/base);
    double y2  = log(inter/base);
    double y3  = log(expend);

    double r1  = log(p1);
    double r2  = log(p2);
    double r3  = log(p3);

    double d0, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11, d12, d13;

    if (missing) {
      cout << "Missing observation at t = " << t << ".\n";
      d0= d1= d2= d3= d4= d5= d6= d7= d8= d9= d10= d11= d12= d13= -9999.;
    }
    else {
      d0  = 1;
      d1  = log(10*p1+6*p2+8*p3);
      d2  = log(income);
      d3  = log(sqfeet);
      d4  = duplex;
      d5  = mobile;
      d6  = cac*log(heatlos);
      d7  = 0; if (wac>0) d7=log(wac);
      d8  = 0; if (hwh>0) d8=log(famsize+1);
      d9  = 0; if (hwh>0 && wash>0) d9=1;
      d10 = 0; if (dry>0) d10=log(famsize+1);
      d11 = 0; if (ref>0) d11=log(ref);
      d12 = 0; if (frez>0) d12=log(frez);
      d13 = range;
    }

    if (strcmp(s.method,"SUR") == 0  || !missing ) {
      data(1,t0)  = (REAL)y1;
      data(2,t0)  = (REAL)y2;
      data(3,t0)  = (REAL)y3;
      data(4,t0)  = (REAL)r1;
      data(5,t0)  = (REAL)r2;
      data(6,t0)  = (REAL)r3;
      data(7,t0)  = (REAL)d0;
      data(8,t0)  = (REAL)d1;
      data(9,t0)  = (REAL)d2;
      data(10,t0) = (REAL)d3;
      data(11,t0) = (REAL)d4;
      data(12,t0) = (REAL)d5;
      data(13,t0) = (REAL)d6;
      data(14,t0) = (REAL)d7;
      data(15,t0) = (REAL)d8;
      data(16,t0) = (REAL)d9;
      data(17,t0) = (REAL)d10;
      data(18,t0) = (REAL)d11;
      data(19,t0) = (REAL)d12;
      data(20,t0) = (REAL)d13;
      t0++;
    }
    else {
      s.n--;
    }
  }

  s.starting="tmp.dat";
   s.to(s.starting);

//  if (strcmp(s.detail,"full") == 0) cout << starbox("/data//_") << data;
                
  return 1;
}

void model::Z(INTEGER t)
{
  INTEGER i;

  if (s.K < 4 || s.K > 17) {cerr << "Error, model::Z.\n"; exit(1);}
  
  for (i = 1; i<=3;   i++)    Z_t[i] = data(i+3,t); 
  for (i = 4; i<=s.K; i++)    Z_t[i] = data(i+3,t); 

}

int model::terminate() { return 0; }
