#include <iostream.h>
#include <iomanip.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <strstream.h>
#include "field.h"
#include "matrix.h"
#define ALPHA 0.02        // alpha
#define NEURONSIZE 12     // trading days to learn from * 4
#define MAINLOOP 1000     // number of times to learn training data
#define PERCENT 0.75      // allways train on last 1 - PRECENT
                          // but do not retrain if it is a bad
                          // prediction
#define INNERLOOP 255     // training data (the rest is data to compare
                          // to predictions)
#define DIFFERENCE 0.0001 // how far off a prediction can be
double logsig(double n);
double purelin(double n);
double hardlim(double d);
double poslin(double n);
 
using namespace std;
 
template<class F>
ostream &operator<<( ostream &os, const Matrix<F>&  m)
{
  Matrix<char*> s( m.size().first, m.size().second );
  int *maxLen = new int[m.size().second];
     //  any entry in matrix column j.
  int colSep = 2;
  for ( size_t j = 0 ; j < m.size().second ; ++j )
  {
    maxLen[j] = 0;
    for ( size_t i = 0 ; i < m.size().first ; ++i )
    {
      ostrstream oss;
      oss.flags( os.flags());
      oss.precision( os.precision());
      oss.fill( os.fill());
      oss << m(i,j) << ends;
      int len = strlen(oss.str());
      s(i,j) = strcpy( new char[len+1], oss.str());
      if ( len > maxLen[j] )
        maxLen[j] = len;
    }
  }
  for ( size_t i = 0 ; i < m.size().first ; ++i )
  {
    os << "|"; // MTS: changed output
    for ( size_t j = 0 ; j < m.size().second ; ++j )
    {
      os << setw( maxLen[j] + (j?colSep:1) ) << s(i,j);
      delete[] s(i,j);
    }
    os << " |" << endl;
  }
  delete[] maxLen;
  return os;
}
 
// ***** Start of Main Program ******
main()
{
  int n,y, z, i = 0, j, testdata, updated;
  char wholeline[128], * tempch;
  double stockarray[1024][5];
 
  cerr << "NN Project" << endl;
 
  cin >> wholeline;
  while(cin >> wholeline) 
  {
    j=0;
    tempch = strtok(wholeline, ",");
    while(tempch = strtok(NULL, ","))
      stockarray[i][j++] = atof(tempch);
    i++;
  }

  testdata = i - INNERLOOP;

  Matrix<double> open_a1(NEURONSIZE, 1);
  Matrix<double> high_a1(NEURONSIZE, 1);
  Matrix<double> low_a1(NEURONSIZE, 1);
  Matrix<double> close_a1(NEURONSIZE, 1);
  Matrix<double> open_w1(NEURONSIZE, NEURONSIZE);
  Matrix<double> high_w1(NEURONSIZE, NEURONSIZE);
  Matrix<double> low_w1(NEURONSIZE, NEURONSIZE);
  Matrix<double> close_w1(NEURONSIZE, NEURONSIZE);
  Matrix<double> p1(NEURONSIZE, 1);
  Matrix<double> p1_day(NEURONSIZE, 1);
  Matrix<double> open_b1(NEURONSIZE, 1);
  Matrix<double> high_b1(NEURONSIZE, 1);
  Matrix<double> low_b1(NEURONSIZE, 1);
  Matrix<double> close_b1(NEURONSIZE, 1);
  Matrix<double> open_a2(1, 1);
  Matrix<double> high_a2(1, 1);
  Matrix<double> low_a2(1, 1);
  Matrix<double> close_a2(1, 1);
  Matrix<double> open_a2_day(1, 1);
  Matrix<double> high_a2_day(1, 1);
  Matrix<double> low_a2_day(1, 1);
  Matrix<double> close_a2_day(1, 1);
  Matrix<double> open_w2(1, NEURONSIZE);
  Matrix<double> high_w2(1, NEURONSIZE);
  Matrix<double> low_w2(1, NEURONSIZE);
  Matrix<double> close_w2(1, NEURONSIZE);
  Matrix<double> open_b2(1, 1);
  Matrix<double> high_b2(1, 1);
  Matrix<double> low_b2(1, 1);
  Matrix<double> close_b2(1, 1);
  Matrix<double> open_t(1, 1);
  Matrix<double> high_t(1, 1);
  Matrix<double> low_t(1, 1);
  Matrix<double> close_t(1, 1);
  Matrix<double> open_err(1, 1);
  Matrix<double> high_err(1, 1);
  Matrix<double> low_err(1, 1);
  Matrix<double> close_err(1, 1);
  Matrix<double> open_s1(NEURONSIZE, 1);
  Matrix<double> high_s1(NEURONSIZE, 1);
  Matrix<double> low_s1(NEURONSIZE, 1);
  Matrix<double> close_s1(NEURONSIZE, 1);
  Matrix<double> open_s2(1, 1);
  Matrix<double> high_s2(1, 1);
  Matrix<double> low_s2(1, 1);
  Matrix<double> close_s2(1, 1);
  Matrix<double> open_f(NEURONSIZE, NEURONSIZE);
  Matrix<double> high_f(NEURONSIZE, NEURONSIZE);
  Matrix<double> low_f(NEURONSIZE, NEURONSIZE);
  Matrix<double> close_f(NEURONSIZE, NEURONSIZE);
  Matrix<double> open_diff(1, 1);
  Matrix<double> high_diff(1, 1);
  Matrix<double> low_diff(1, 1);
  Matrix<double> close_diff(1, 1);

  for(i = 0; i < NEURONSIZE; i++)
    for(j = 0; j < NEURONSIZE; j++)
    {
      open_w1(i, j)  = rand() % 10;
      high_w1(i, j)  = rand() % 10;
      low_w1(i, j)   = rand() % 10;
      close_w1(i, j) = rand() % 10;
    }
 
  for(i = 0; i < NEURONSIZE; i++)
  {
    open_b1(i, 0)  = 1;
    high_b1(i, 0)  = 1;
    low_b1(i, 0)   = 1;
    close_b1(i, 0) = 1;
  }
 
  for(i = 0; i < NEURONSIZE; i++)
  {
    open_w2(0, i)  = rand() % 10;
    high_w2(0, i)  = rand() % 10;
    low_w2(0, i)   = rand() % 10;
    close_w2(0, i) = rand() % 10;
  }

  open_b2(0,0)  = 1;
  high_b2(0,0)  = 1;
  low_b2(0,0)   = 1;
  close_b2(0,0) = 1;
 
  for(y = 0; y < MAINLOOP; y++)
  {
    cerr << "Learning loop " << y + 1 <<  " of " << MAINLOOP << endl; 
    for(z = NEURONSIZE / 4; z < (INNERLOOP - 1); z++)
    {
      if ((z + 1 - NEURONSIZE / 4) % 25 == 0)
        cerr << "Data loop " << z + 1 - NEURONSIZE / 4 
             <<  " of " << INNERLOOP - NEURONSIZE / 4 - 1 << endl; 
      do { 
        updated = 0;
        for(i = 0; i < NEURONSIZE; i += 4)
        {
          p1(i, 0)     = stockarray[z - (i/4)][0] 
                         - stockarray[z - (i/4) - 1][0];
          p1(i + 1, 0) = stockarray[z - (i/4)][1]
                         - stockarray[z - (i/4) - 1][1];
          p1(i + 2, 0) = stockarray[z - (i/4)][2] 
                         - stockarray[z - (i/4) - 1][2];
          p1(i + 3, 0) = stockarray[z - (i/4)][3]
                         - stockarray[z - (i/4) - 1][3];
        }

        open_a1 = open_w1 * p1 + open_b1;
        high_a1 = high_w1 * p1 + high_b1;
        low_a1 = low_w1 * p1 + low_b1;
        close_a1 = close_w1 * p1 + close_b1;

        for(i = 0; i < NEURONSIZE; i++) {
          open_a1(i, 0) = logsig(open_a1(i, 0));
          high_a1(i, 0) = logsig(high_a1(i, 0));
          low_a1(i, 0) = logsig(low_a1(i, 0));
          close_a1(i, 0) = logsig(close_a1(i, 0));
        }

        // purelin well thats easy...
        open_a2 = open_w2 * open_a1 + open_b2;
        high_a2 = high_w2 * high_a1 + high_b2;
        low_a2 = low_w2 * low_a1 + low_b2;
        close_a2 = close_w2 * close_a1 + close_b2;

        open_t(0, 0)  = stockarray[z + 1][0] - stockarray[z][0];
        open_err      = open_t - open_a2;
        high_t(0, 0)  = stockarray[z + 1][1] - stockarray[z][1];
        high_err      = high_t - high_a2;
        low_t(0, 0)   = stockarray[z + 1][2] - stockarray[z][2];
        low_err       = low_t - low_a2;
        close_t(0, 0) = stockarray[z + 1][3] - stockarray[z][3];
        close_err     = close_t - close_a2;

        for(i = 0; i < NEURONSIZE; i++)
          for(j = 0; j < NEURONSIZE; j++)
          {
            open_f(i, j)  = 0.0;
            high_f(i, j)  = 0.0;
            low_f(i, j)   = 0.0;
            close_f(i, j) = 0.0;
          }

        for(i = 0; i < NEURONSIZE; i++)
        {
          open_f(i, i)  = (1 - open_a1(i, 0)) * open_a1(i, 0);
          high_f(i, i)  = (1 - high_a1(i, 0)) * high_a1(i, 0);
          low_f(i, i)   = (1 - low_a1(i, 0)) * low_a1(i, 0);
          close_f(i, i) = (1 - close_a1(i, 0)) * close_a1(i, 0);
        } 

        open_diff(0,0) = open_err(0,0);
        if (open_diff(0,0) < 0)
          open_diff(0,0) = -1 * open_diff(0,0);
        if (open_diff(0,0) > DIFFERENCE || y >= MAINLOOP * PERCENT)
        {
          open_s2 = -2 * open_err;
          open_s1 = open_f * open_w2.transpose() * open_s2;
          open_w2 = open_w2 - ALPHA * open_s2 * open_a1.transpose();
          open_b2 = open_b2 - ALPHA * open_s2;
          open_w1 = open_w1 - ALPHA * open_s1 * p1.transpose();
          open_b1 = open_b1 - ALPHA * open_s1;
          updated = 1;
        }
        high_diff(0,0) = high_err(0,0);
        if (high_diff(0,0) < 0)
          high_diff(0,0) = -1 * high_diff(0,0);
        if (high_diff(0,0) > DIFFERENCE || y >= MAINLOOP * PERCENT)
        {
          high_s2 = -2 * high_err;
          high_s1 = high_f * high_w2.transpose() * high_s2;
          high_w2 = high_w2 - ALPHA * high_s2 * high_a1.transpose();
          high_b2 = high_b2 - ALPHA * high_s2;
          high_w1 = high_w1 - ALPHA * high_s1 * p1.transpose();
          high_b1 = high_b1 - ALPHA * high_s1;
          updated = 1;
        }
        low_diff(0,0) = low_err(0,0);
        if (low_diff(0,0) < 0)
          low_diff(0,0) = -1 * low_diff(0,0);
        if (low_diff(0,0) > DIFFERENCE || y >= MAINLOOP * PERCENT)
        {
          low_s2 = -2 * low_err;
          low_s1 = low_f * low_w2.transpose() * low_s2;
          low_w2 = low_w2 - ALPHA * low_s2 * low_a1.transpose();
          low_b2 = low_b2 - ALPHA * low_s2;
          low_w1 = low_w1 - ALPHA * low_s1 * p1.transpose();
          low_b1 = low_b1 - ALPHA * low_s1;
          updated = 1;
        }
        close_diff(0,0) = close_err(0,0);
        if (close_diff(0,0) < 0)
          close_diff(0,0) = -1 * close_diff(0,0);
        if (close_diff(0,0) > DIFFERENCE || y >= MAINLOOP * PERCENT)
        {
          close_s2 = -2 * close_err;
          close_s1 = close_f * close_w2.transpose() * close_s2;
          close_w2 = close_w2 - ALPHA * close_s2 * close_a1.transpose();
          close_b2 = close_b2 - ALPHA * close_s2;
          close_w1 = close_w1 - ALPHA * close_s1 * p1.transpose();
          close_b1 = close_b1 - ALPHA * close_s1;
          updated = 1;
        } 
      } while(updated && y < MAINLOOP * PERCENT);
    }
  }

  for(i = 0; i < NEURONSIZE; i+=4)
  {
    p1(i, 0)     = stockarray[INNERLOOP - (i/4) - 1][0] 
                   - stockarray[INNERLOOP - (i/4) - 2][0];
    p1(i + 1, 0) = stockarray[INNERLOOP - (i/4) - 1][1]
                   - stockarray[INNERLOOP - (i/4) - 2][1];
    p1(i + 2, 0) = stockarray[INNERLOOP - (i/4) - 1][2] 
                   - stockarray[INNERLOOP - (i/4) - 2][2];
    p1(i + 3, 0) = stockarray[INNERLOOP - (i/4) - 1][3] 
                   - stockarray[INNERLOOP - (i/4) - 2][3];
    p1_day(i, 0) = stockarray[INNERLOOP - (i/4) - 1][0] 
                   - stockarray[INNERLOOP - (i/4) - 2][0];
    p1_day(i + 1, 0) = stockarray[INNERLOOP - (i/4) - 1][1]
                       - stockarray[INNERLOOP - (i/4) - 2][1];
    p1_day(i + 2, 0) = stockarray[INNERLOOP - (i/4) - 1][2] 
                       - stockarray[INNERLOOP - (i/4) - 2][2];
    p1_day(i + 3, 0) = stockarray[INNERLOOP - (i/4) - 1][3] 
                       - stockarray[INNERLOOP - (i/4) - 2][3];
  }

  double lastopen = stockarray[INNERLOOP - 1][0];
  double lasthigh = stockarray[INNERLOOP - 1][1];
  double lastlow = stockarray[INNERLOOP - 1][2];
  double lastclose = stockarray[INNERLOOP -1][3];
  double openday, highday, lowday, closeday;
  
  cout << "pnum,"
       << "popen,popenday,open,popenerr,popendayerr,"
       << "phigh,phighday,high,phigherr,phighdayerr,"
       << "plow,plowday,low,plowerr,plowdayerr,"
       << "pclose,pcloseday,close,pcloseerr,pclosedayerr\n";

  for (z = 0; z < testdata; z++)
  {
    open_a1 = open_w1 * p1 + open_b1;
    for(i = 0; i < NEURONSIZE; i++)
      open_a1(i, 0) = logsig(open_a1(i, 0));
    open_a2 = open_w2 * open_a1 + open_b2;
    open_a1 = open_w1 * p1_day + open_b1;
    for(i = 0; i < NEURONSIZE; i++)
      open_a1(i, 0) = logsig(open_a1(i, 0));
    open_a2_day = open_w2 * open_a1 + open_b2;

    high_a1 = high_w1 * p1 + high_b1;
    for(i = 0; i < NEURONSIZE; i++)
      high_a1(i, 0) = logsig(high_a1(i, 0));
    high_a2 = high_w2 * high_a1 + high_b2;
    high_a1 = high_w1 * p1_day + high_b1;
    for(i = 0; i < NEURONSIZE; i++)
      high_a1(i, 0) = logsig(high_a1(i, 0));
    high_a2_day = high_w2 * high_a1 + high_b2;

    low_a1 = low_w1 * p1 + low_b1;
    for(i = 0; i < NEURONSIZE; i++)
      low_a1(i, 0) = logsig(low_a1(i, 0));
    low_a2 = low_w2 * low_a1 + low_b2;
    low_a1 = low_w1 * p1_day + low_b1;
    for(i = 0; i < NEURONSIZE; i++)
      low_a1(i, 0) = logsig(low_a1(i, 0));
    low_a2_day = low_w2 * low_a1 + low_b2;

    close_a1 = close_w1 * p1 + close_b1;
    for(i = 0; i < NEURONSIZE; i++)
      close_a1(i, 0) = logsig(close_a1(i, 0));
    close_a2 = close_w2 * close_a1 + close_b2;
    close_a1 = close_w1 * p1_day + close_b1;
    for(i = 0; i < NEURONSIZE; i++)
      close_a1(i, 0) = logsig(close_a1(i, 0));
    close_a2_day = close_w2 * close_a1 + close_b2;

    lastopen += open_a2(0, 0);
    lasthigh += high_a2(0, 0);
    lastlow += low_a2(0, 0);
    lastclose += close_a2(0, 0);
    openday = stockarray[INNERLOOP + z - 1][0] + open_a2_day(0,0);
    highday = stockarray[INNERLOOP + z - 1][1] + high_a2_day(0,0);
    lowday = stockarray[INNERLOOP + z - 1][2] + low_a2_day(0,0);
    closeday = stockarray[INNERLOOP + z - 1][3] + close_a2_day(0,0);

    cout << z << "," << lastopen  
         << "," << openday
         << "," << stockarray[INNERLOOP + z][0]
         << "," << lastopen - stockarray[INNERLOOP + z][0]
         << "," << openday - stockarray[INNERLOOP + z][0]
         << "," << lasthigh 
         << "," << highday
         << "," << stockarray[INNERLOOP + z][1] 
         << "," << lasthigh - stockarray[INNERLOOP + z][1]
         << "," << highday - stockarray[INNERLOOP + z][1]
         << "," << lastlow 
         << "," << lowday
         << "," << stockarray[INNERLOOP + z][2]   
         << "," << lastlow - stockarray[INNERLOOP + z][2]
         << "," << lowday - stockarray[INNERLOOP + z][2]
         << "," << lastclose 
         << "," << closeday
         << "," << stockarray[INNERLOOP + z][3]
         << "," << lastclose - stockarray[INNERLOOP + z][3]
         << "," << closeday - stockarray[INNERLOOP + z][3] << endl;

    for (i = NEURONSIZE - 1; i > 3; i--)
    {
      p1(i, 0) = p1(i - 4, 0);
      p1_day(i, 0) = p1(i - 4, 0);
    }

    p1(0, 0) = open_a2(0, 0);
    p1(1, 0) = high_a2(0, 0);
    p1(2, 0) = low_a2(0, 0);
    p1(3, 0) = close_a2(0, 0);
    p1_day(0, 0) = stockarray[INNERLOOP + z][0] 
                   - stockarray[INNERLOOP + z - 1][0];
    p1_day(1, 0) = stockarray[INNERLOOP + z][1] 
                   - stockarray[INNERLOOP + z - 1][1];
    p1_day(2, 0) = stockarray[INNERLOOP + z][2] 
                   - stockarray[INNERLOOP + z - 1][2];
    p1_day(3, 0) = stockarray[INNERLOOP + z][3] 
                   - stockarray[INNERLOOP + z - 1][3];
  }
}
 
double logsig(double n)
{
  return(1.0/(1.0 + exp(n * -1.0)));
}