/*----------------------------------------------------------------------------*
\Module{Concavity Test}

This program takes as its input (through stdin) lines of data points,
comma delimited.  It assumes that these data points are points at which 
all first derivatives are zero.  Data points must be formatted as:

M1,M2,M3,...,B1,B2,B3,...

The neuron numbers are those specified in the config file.  Run this
with the same config file that was used to run the simulation and
everything will work out correctly.  It then calculates the Hessian
matrix, and finds whether it is definite or indefinite.  An indefinite
matrix means that the point is a saddle, and a definite one means it
is an extrema.  This is then output to stdout.  

 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "neural.h"
#include "protos.h"
#include "concavity.h"

int main(int argc, char *argv[]) {
  int i,j;
  Org* o;
  // Load the config file
  char *configfile;
  configfile = (char *) 0;
  
  for (i = 1; i <  argc; i++)
    if (argv[i][0] != '-')
      configfile = argv[1];

  if (!configfile)
  {
    printf("usage: concavity [options] inputfile\n");
    return 1;
  }

  if (!ReadConfiguration(configfile) |
      !ReadArgConfiguration(argc, argv))
    return 1;

  DisplayConfiguration();

  // Initialize the random seed
  RandStart(c.seed);

  // Figure out how many mutable synapses were present in the config file
  int nsynapses = 0;
  for (i=0; i<c.M->n; i++)
    if (c.M->m[i].mut == TRUE)
      nsynapses++;
  for (i=0; i<c.B->n; i++)
    if (c.B->m[i].mut == TRUE)
      nsynapses++;

  // Read the input from stdin, stopping at the newline character or EOF
  #define MAXLINE 1024 // No more than 1024 chars on a line
  char c;
  char *ch;
  char buff[MAXLINE];
  real coordinates[nsynapses];
  // Continue looping until we have reached the end of the stream.
  while (c != EOF) {
    for (i=0; i<MAXLINE; i++) {
      c = getc(stdin); // Using getc instead of gets to avoid a buffer overflow
      if (c == EOF || c == '\n') {
        buff[i] = '\0';
        break;
      }
      if (i == MAXLINE-1)
        printf("Error, input line too long.\n"), exit(5);
      buff[i] = c;
    }
    // Exit the program successfully on an empty line.
    if (strlen(buff) == 0) return 0;


    // We assume here that the number of columns in the input data is
    // equal to the number of synapses.  All coordinates from the
    // input line are split by the "," character and stored in
    // `coordinates`.
    ch = strtok(buff, ",");
    for (i=0; i<nsynapses; i++) {
      if (ch == NULL)
        printf("Error, insufficient coordinates.\n"), exit(7);
      coordinates[i] = atof(ch);
      ch = strtok(NULL, ",");
    }

    // Now we create the Hessian matrix for the input coordinates.
    Matrix* hessian = MatCreate(nsynapses, nsynapses);
    hessian->n = 0;
    createHessian(coordinates, nsynapses, hessian);
    MatDisplay(hessian, "Hessian");

  }
  return 0;
}

/* 
ENTRY CONDITIONS: 
    `coords` holds a list of coordinates that is `ncoords` long
    `hess` is a pointer to a matrix that has the dimensions `ncoords`
        by `ncoords`
EXIT CONDITIONS:
    `hess` will be filled with an approximated Hessian matrix at that point
    Returns 0 on success or 1 on failure
*/
#define _H 0.5
int createHessian(real coords[], int ncoords, Matrix *hess) {
  // No one portion of the Hessian depends on values used to calculate
  // another portion, so we do not need to cache values.  Hence, we
  // can calculate the matrix element by element.  It is symmetrical,
  // so we only need to calculate the elements on or above the
  // diagonal.  See 4FE1DC0C.n for a more complete discussion.
  int i,j;
  Org *o;
  // First do the diagonals.  These are just the second derivatives.
  // NOTE: This would be a good place to verify that the first
  // derivatives are 0.
  for (i=0; i<ncoords; i++) {
    // Calculate fitness at the point
    real fi = orgFitness(coords, ncoords);
    // Calculate the fitness + and - _H
    coords[i] += _H;
    real fiplus = orgFitness(coords, ncoords);
    coords[i] -= 2*_H;
    real fiminus = orgFitness(coords, ncoords);
    // Reset to the original value for the next iteration.
    coords[i] += _H;
    // Now calculate the Hessian value and store it
    real hessval = (fiplus + fiminus - 2*fi)/(2*_H*_H);
    MatStore(hess, i, i, hessval, TRUE);
    printf("Stored value %i,%i\n", i, i);
    printf("First derivative is: %.8f", (fiplus-fiminus)/(2*_H));
  }
  // Now, non-diagonals
  for (i=0; i<ncoords; i++) {
    for (j=i+1; j<ncoords; j++) {
      // Calculate all necessary perturbations
      coords[i] += _H;
      coords[j] += _H;
      real fiplusjplus = orgFitness(coords, ncoords);
      coords[j] -= 2*_H;
      real fiplusjminus = orgFitness(coords, ncoords);
      coords[i] -= 2*_H;
      real fiminusjminus = orgFitness(coords, ncoords);
      coords[j] += 2*_H;
      real fiminusjplus = orgFitness(coords, ncoords);
      // Reset to the original value for the next iteration.
      coords[j] -= _H;
      coords[i] += _H;
      // Now calculate the Hessian value and store it
      real hessval = (fiplusjplus + fiminusjminus - fiplusjminus - fiminusjplus)/(4*_H*_H);
      MatStore(hess, i, j, hessval, TRUE);
      MatStore(hess, j, i, hessval, TRUE); // Symmetric matrix
      printf("Stored value %i,%i\n", i, j);
    }
  }
  calc_concavity_with_octave(hess, "Hesstest");
  // Success
  return 0;
}


/*
ENTRY CONDITIONS:
    `coords` is an array of M values in the order to be loaded
        directly into an organism, followed by a similar array of B
        values.  It must be `ncoords` long.
EXIT CONDITIONS:
    Returns the fitness of the organism.
*/

real orgFitness(real coords[], int ncoords) {
  Org *o;
  real fitness;
  int i,j;
  o = OrgCreate(1);
  OrgInit(o);

  // We assume here for simplicity's sake that the coordinate inputs
  // are in the same order as in `o`.  (i.e. o->ploid[0].M->m[i].r
  // is the same as coords[i].)
  int k=0;
  for (i=0; i<o->ploid[0].M->n; i++) 
    if (o->ploid[0].M->m[i].mut)
      o->ploid[0].M->m[i].r = coords[k], k++;
  for (i=0; i<o->ploid[0].B->n; i++)
    if (o->ploid[0].B->m[i].mut)
      o->ploid[0].B->m[i].r = coords[k], k++;

  if (ncoords != k)
    printf("Error, check to make sure you are using the right config file.\n"), exit(8);
      
  // OrgEval takes us through all of the iterations necessary for
  // the organism, and puts the results in o->wdist, o->edist, and
  // o->dist.
  OrgEval(o);

  fitness = o->dist;
  printf("Fitness is %f", fitness);
  OrgDestroy(o);
  return fitness;
}

/* calc_concavity_with_octave(Matrix M)
ENTRY CONDITIONS: 
    `M` contains the hessian matrix at a point.
    `ident` contains some sort of identifier string that allows this
        point to be distinguished from others when many are output.
        This will just be a string containing the coordinates in most
        cases.
EXIT CONDITIONS:
    Prints the Matlab/Octave code necessary to calculate the concavity
        of the point.  If all of the eigenvalues are real and
        positive, it is a minimum.  If they are all real and negative,
        it is a maximum.  Otherwise, it is a saddle.
*/
void calc_concavity_with_octave(Matrix *M, char ident[]) {
  printf("printf(\"%s:\");\n", ident);
  printf("M = [");
  int i,j;
  for (i=0; i<M->di; i++) {
    for (j=0; j<M->dj; j++)
      printf(" %.7f", (double)MatLookup(M, i, j));
    printf("\n");
  }
  printf("]\n");
  printf("evs = eig(M);\n"
         "if ! isreal(evs)\n"
         "    printf(\"Saddle point\\n\");\n"
         "elseif min(evs) > 0\n"
         "    printf(\"Minimum\\n\");\n"
         "elseif max(evs) < 0\n"
         "    printf(\"Maximum\\n\");\n"
         "else\n"
         "    printf(\"Saddle point\\n\");\n"
         "endif;\n\n");
}