/* ========================================================================== *
 *  Evolutionary cascade-correlation.
 *
 *  Author: Firas Safadi (fsafadi@ulg.ac.be)
 * ========================================================================== */

#include "CascadeCorrelation.h"
#include "GeneticAlgorithm.h"
#include "CrossEntropy.h"

#define DEFAULT_OPTIMIZER GeneticAlgorithmEAID

#ifndef DOUBLE_PRECISION
#define FABS fabsf
#define CUBLAS_GEMM cublasSgemm
#define CUBLAS_GEMV cublasSgemv
#define CUBLAS_COPY cublasScopy
#else
#define FABS fabs
#define CUBLAS_GEMM cublasDgemm
#define CUBLAS_GEMV cublasDgemv
#define CUBLAS_COPY cublasDcopy
#endif

#define ACTFUN(x) ((x) / (1 + FABS(x)))

CCSimulator::CCSimulator(CascadeCorrelation* cc, CCSimFunc sim)
{
  if (!cc || !sim)
    error("Simulator cannot be null");

  m_pCascadeCorrelation = cc;
  m_fSim = sim;
}

CCSimulator::~CCSimulator()
{

}

void CCSimulator::sim(EvolutionaryAlgorithm* ea) {
  (m_pCascadeCorrelation->*m_fSim)(ea);
}

CascadeCorrelation::CascadeCorrelation()
{
  m_pUnitEA = NULL;
  m_pOutputEA = NULL;

  m_pUnitSimulator = new CCSimulator(this, &CascadeCorrelation::simUnit);
  m_pOutputSimulator = new CCSimulator(this, &CascadeCorrelation::simOutput);

  m_iSamples = 0;
  m_iInputDim = 0;
  m_iOutputDim = 0;

  m_dPropTestSamples = 0.0;
  m_iTrainingSamples = 0;
  m_iTestSamples = 0;

  m_iMaxHiddenLayers = 0;
  m_iHiddenLayers = 0;

  m_dBias = 1.0;

  m_dErrorGoal = 0.0;
  m_dTrainingError = 0.0;
  m_dTestError = 0.0;

  m_bNetworkBuilt = false;
  m_bNetworkTrained = false;
  m_bDataLoaded = false;

  m_bRandomizeTestSamples = true;
  m_bSilentMode = false;

  m_iTopPopulation = 0;
  m_iTopIndividual = 0;
  m_dTopFitness = 0.0;

  m_iOptimalLayerCount = 0;

  m_iOverfitCounter = 0;
  m_iOverfitStop = 0;

  m_iUnitEAPopCount = 0;
  m_iUnitEAPopSize = 0;

  m_iOutputEAPopCount = 0;
  m_iOutputEAPopSize = 0;

  m_eUnitEAID = DEFAULT_OPTIMIZER;
  m_eOutputEAID = DEFAULT_OPTIMIZER;

  cublas_call(cublasCreate(&m_pHandle));
}

CascadeCorrelation::~CascadeCorrelation()
{
  if (m_bNetworkBuilt) {
    delete m_pUnitEA;
    delete m_pOutputEA;
  }

  delete m_pUnitSimulator;
  delete m_pOutputSimulator;

  cublas_call(cublasDestroy(m_pHandle));
}

void CascadeCorrelation::setNetworkParameters(nat samples, real testProp, nat inputDim, nat outputDim, real bias, real errorGoal, nat maxHiddenLayers, nat overfitStop, bool randomTestSamples)
{
  if (m_bNetworkBuilt)
    error("Cannot change existing network parameters");

  m_iSamples = samples;
  m_dPropTestSamples = testProp;

  m_iInputDim = inputDim;
  m_iOutputDim = outputDim;

  m_dBias = bias;

  m_dErrorGoal = errorGoal;

  m_iMaxHiddenLayers = maxHiddenLayers;

  m_iOverfitStop = overfitStop;

  m_bRandomizeTestSamples = randomTestSamples;

  m_iTestSamples = static_cast<nat>(m_iSamples*m_dPropTestSamples);
  m_iTrainingSamples = m_iSamples - m_iTestSamples;

  validateNetworkParameters();
}

void CascadeCorrelation::setUnitEAParameters(nat popCount, nat popSize, EAID id)
{
  m_iUnitEAPopCount = popCount;
  m_iUnitEAPopSize = popSize;
  m_eUnitEAID = id;
}

void CascadeCorrelation::setOutputEAParameters(nat popCount, nat popSize, EAID id)
{
  m_iOutputEAPopCount = popCount;
  m_iOutputEAPopSize = popSize;
  m_eOutputEAID = id;
}

void CascadeCorrelation::validateNetworkParameters()
{
  if (m_iSamples == 0 || m_iInputDim == 0 || m_iOutputDim == 0 || m_iMaxHiddenLayers == 0 || m_iOverfitStop == 0)
    error("Invalid network parameters");

  if (m_dPropTestSamples < 0.0 || m_dPropTestSamples > 1.0)
    error("Test samples proportion cannot be lower than 0.0 or greater than 1.0");

  if (m_iTrainingSamples == 0 || m_iTestSamples == 0)
    error("At least one training sample and one test sample are required");

  if (m_iOutputDim != 1)
    error("Only one output variable is allowed (regression mode)");
}

void CascadeCorrelation::buildNetwork()
{
  if (m_bNetworkBuilt)
    error("Cannot build network: already exists");

  if (m_pUnitEA || m_pOutputEA)
    error("Cannot build network: stray optimizer");

  validateNetworkParameters();

  //Create optimizers
  m_pUnitEA = makeEA(m_iUnitEAPopCount, m_iUnitEAPopSize, m_eUnitEAID);
  m_pOutputEA = makeEA(m_iOutputEAPopCount, m_iOutputEAPopSize, m_eOutputEAID);

  //Set the simulators
  m_pUnitEA->setSimulator(m_pUnitSimulator);
  m_pOutputEA->setSimulator(m_pOutputSimulator);

  //Configure genome pregeneration
  m_pUnitEA->disablePregen();
  m_pOutputEA->enablePregen();

  //The cache stores bias, input vector and hidden layer outputs for each sample
  m_pTrainingCache.resize(m_iTrainingSamples, maxCacheWidth());
  m_pTestCache.resize(m_iTestSamples, maxCacheWidth());

  //Target vectors
  m_pTrainingTargets.resize(1, m_iTrainingSamples);
  m_pTestTargets.resize(1, m_iTestSamples);

  //Training error is used as a buffer for host-device data transfers
  m_pTrainingError.resize(1, m_iTrainingSamples);

  //Network weights
  m_pUnitWeights.resize(m_iMaxHiddenLayers, maxCacheWidth());
  m_pOutputWeights.resize(1, maxCacheWidth());

  //Simulation matrices are used to store the result of the matrix multiplication at the heart of the population evaluation process
  m_pTrainingSimulation.resize(m_iTrainingSamples, MAX(m_pUnitEA->totalInd(), m_pOutputEA->totalInd()));
  m_pTestSimulation.resize(1, m_iTestSamples);

  //History matrices are used to keep track of previous network states
  m_pOutputWeightsHistory.resize(1 + m_iMaxHiddenLayers, maxCacheWidth());
  m_pErrorHistory.resize(2, 1 + m_iMaxHiddenLayers);

  //Keeps track of original sample numbers
  m_pTestSetIDs.resize(1, m_iTestSamples);

  //Spawn everything
  m_pTrainingCache.spawn(Matrix::HostAndDevice);
  m_pTestCache.spawn(Matrix::HostAndDevice);

  m_pTrainingTargets.spawn(Matrix::HostAndDevice);
  m_pTestTargets.spawn(Matrix::HostAndDevice);

  m_pTrainingError.spawn(Matrix::HostAndDevice);

  m_pUnitWeights.spawn(Matrix::HostAndDevice);
  m_pOutputWeights.spawn(Matrix::HostAndDevice);

  m_pTrainingSimulation.spawn(Matrix::DeviceOnly);
  m_pTestSimulation.spawn(Matrix::HostAndDevice);

  m_pOutputWeightsHistory.spawn(Matrix::HostOnly);
  m_pErrorHistory.spawn(Matrix::HostOnly);

  m_pTestSetIDs.spawn(Natrix::HostOnly);

  //Initialize matrices to avoid data persistence between successive runs
  m_pTrainingCache.init(0.0);
  m_pTestCache.init(0.0);

  m_pTrainingTargets.init(0.0);
  m_pTestTargets.init(0.0);

  m_pTrainingError.init(0.0);

  m_pUnitWeights.init(0.0);
  m_pOutputWeights.init(0.0);

  m_pTrainingSimulation.initDevice(0.0);
  m_pTestSimulation.init(0.0);

  m_pOutputWeightsHistory.initHost(0.0);
  m_pErrorHistory.initHost(0.0);

  m_pTestSetIDs.initHost(0);

  m_bNetworkBuilt = true;
}

void CascadeCorrelation::loadDataSet(char* path)
{
  if (!m_bNetworkBuilt)
    error("Cannot load dataset: network not built");

  if (m_bDataLoaded)
    error("Cannot load dataset: already loaded");

  File data(path);

  bool* inTestSet = new bool[m_iSamples];

  //Flag each sample for either training or testing
  std::fill(inTestSet, inTestSet + m_iTrainingSamples, false);
  std::fill(inTestSet + m_iTrainingSamples, inTestSet + m_iSamples, true);

  if (m_bRandomizeTestSamples)
    std::random_shuffle(inTestSet, inTestSet + m_iSamples);

  nat trainingCount = 0;
  nat testCount = 0;

  //Read samples from file
  for (nat i = 0; i < m_iSamples; i++) {
    real* cache;
    real* target;
    
    if (inTestSet[i]) {
      m_pTestSetIDs(0, testCount) = i;
      cache = m_pTestCache.hr(testCount);
      target = m_pTestTargets.hp() + testCount;
      testCount++;
    }
    else {
      cache = m_pTrainingCache.hr(trainingCount);
      target = m_pTrainingTargets.hp() + trainingCount;
      trainingCount++;
    }

    //Write bias
    cache[0] = m_dBias;

    //Read input vector
    for (nat j = 1; j <= m_iInputDim; j++)
      data >> cache[j];

    //Read output value
    data >> *target;
  }

  data.close();

  delete[] inTestSet;

  //Copy samples to device
  m_pTrainingCache.syncDeviceWithHost();
  m_pTestCache.syncDeviceWithHost();

  m_pTrainingTargets.syncDeviceWithHost();
  m_pTestTargets.syncDeviceWithHost();

  m_bDataLoaded = true;
}

void CascadeCorrelation::loadDataSet(Matrix& dataset, bool* inTestSet)
{
  if (!m_bNetworkBuilt)
    error("Cannot load dataset: network not built");

  if (m_bDataLoaded)
    error("Cannot load dataset: already loaded");

  bool userPartition = true;

  if (inTestSet == NULL) {
    inTestSet = new bool[m_iSamples];

    //Flag each sample for either training or testing
    std::fill(inTestSet, inTestSet + m_iTrainingSamples, false);
    std::fill(inTestSet + m_iTrainingSamples, inTestSet + m_iSamples, true);

    if (m_bRandomizeTestSamples)
      std::random_shuffle(inTestSet, inTestSet + m_iSamples);

    userPartition = false;
  }

  nat trainingCount = 0;
  nat testCount = 0;

  //Copy samples from matrix
  for (nat i = 0; i < m_iSamples; i++) {
    real* cache;
    real* target;

    if (inTestSet[i]) {
      m_pTestSetIDs(0, testCount) = i;
      cache = m_pTestCache.hr(testCount);
      target = m_pTestTargets.hp() + testCount;
      testCount++;
    }
    else {
      cache = m_pTrainingCache.hr(trainingCount);
      target = m_pTrainingTargets.hp() + trainingCount;
      trainingCount++;
    }

    //Write bias
    cache[0] = m_dBias;

    //Copy input vector
    memcpy(cache + 1, dataset.hr(i), m_iInputDim*dataset.cellSize());

    //Copy output value
    *target = dataset(i, m_iInputDim);
  }

  if (!userPartition)
    delete[] inTestSet;

  //Copy samples to device
  m_pTrainingCache.syncDeviceWithHost();
  m_pTestCache.syncDeviceWithHost();

  m_pTrainingTargets.syncDeviceWithHost();
  m_pTestTargets.syncDeviceWithHost();

  m_bDataLoaded = true;
}

void CascadeCorrelation::trainNetwork()
{
  if (!m_bNetworkBuilt)
    error("Cannot train network: not built");

  if (!m_bDataLoaded)
    error("Cannot train network: no data");

  if (m_bNetworkTrained)
    error("Cannot train network: already trained");

  if (!m_bSilentMode) {
    std::cout << std::endl;
    std::cout << "Output EA populations: " << m_pOutputEA->popCount() << std::endl;
    std::cout << "Output EA population size: " << m_pOutputEA->popSize() << std::endl;
    std::cout << "Output EA threads per block: " << m_pOutputEA->threadsPerBlock() << std::endl;

    std::cout << std::endl;
    std::cout << "Unit EA populations: " << m_pUnitEA->popCount() << std::endl;
    std::cout << "Unit EA population size: " << m_pUnitEA->popSize() << std::endl;
    std::cout << "Unit EA threads per block: " << m_pUnitEA->threadsPerBlock() << std::endl;

    std::cout << std::endl;
    std::cout << "Training network [error goal = " << m_dErrorGoal << ", max hidden layers = " << m_iMaxHiddenLayers << ", overfit stop = " << m_iOverfitStop << "]" << std::endl;
    std::cout << std::endl;
  }

  //Hatch output population
  m_pOutputEA->hatchSpecies();

  //Train network once before adding hidden layers
  trainOutputWeights();

  //Add hidden layer and train network again until error goal is reached
  while (m_dTrainingError > m_dErrorGoal && m_iHiddenLayers < m_iMaxHiddenLayers && m_iOverfitCounter < m_iOverfitStop && !_kbhit()) {

    //Add hidden layer
    addHiddenUnit();

    //Add a new gene to output EA
    m_pOutputEA->incGenSize();

    //Train network
    trainOutputWeights();
  }

  //Copy network to host
  m_pOutputWeights.syncHostWithDevice();
  m_pUnitWeights.syncHostWithDevice();

  if (!m_bSilentMode) {
    std::cout << std::endl;
    std::cout << "Training stopped [" << (m_dTrainingError > m_dErrorGoal ? m_iHiddenLayers < m_iMaxHiddenLayers ? m_iOverfitCounter < m_iOverfitStop ?
      "interrupted by user" : "overfit stop triggered" : "max hidden layers reached" : "error goal reached") << "]" << std::endl;

    std::cout << std::endl;
    std::cout << "Top network performance achieved with " << m_iOptimalLayerCount << " hidden layers (MSE = " << m_pErrorHistory(1, m_iOptimalLayerCount) << ")" << std::endl;
  }

  m_bNetworkTrained = true;
}

__global__ void sim_output_K(real* sim, nat simLdt, real* tar, real* fit, nat* board, nat samples, nat totalInd, nat popSize)
{
  nat x = threadIdx.x + blockIdx.x*blockDim.x;

  if (x < totalInd) {
    board[x] = x & (popSize - 1);

    sim += x;

    real* end = sim + samples*simLdt;

    real f = 0.0;

    while (sim < end) {
      f += FABS(*sim - *tar);
      sim += simLdt;
      tar++;
    }

    fit[x] = 1.0/(1.0 + f);
  }
}

void CascadeCorrelation::simOutput(EvolutionaryAlgorithm* ea)
{
  real alpha = 1.0, beta = 0.0;

  //Simulate individual results for all training samples
  cublas_call(CUBLAS_GEMM(m_pHandle, CUBLAS_OP_N, CUBLAS_OP_N,
    ea->totalInd(), m_iTrainingSamples, cacheWidth(),
    &alpha,
    ea->getCurrPop().dp(), ea->getCurrPop().ldt(),
    m_pTrainingCache.dp(), m_pTrainingCache.ldt(),
    &beta,
    m_pTrainingSimulation.dp(), m_pTrainingSimulation.ldt()));

  //Compute individual fitness
  AsyncKernelLaunch(sim_output_K, ea->blockCount(), ea->threadsPerBlock(),
    m_pTrainingSimulation.dp(), m_pTrainingSimulation.ldt(), m_pTrainingTargets.dp(), ea->getCurrFit().dp(), ea->getCurrBoard().dp(), m_iTrainingSamples, ea->totalInd(), ea->popSize());
}

void CascadeCorrelation::trainOutputWeights()
{
  //Evolve individuals
  m_pOutputEA->evolveSpecies();

  //Find top performer
  m_pOutputEA->topPerf(m_iTopPopulation, m_iTopIndividual, m_dTopFitness);

  //Copy top individual to output weights
  cublas_call(CUBLAS_COPY(m_pHandle, cacheWidth(), m_pOutputEA->getCurrPop().dp() + topOutputInd(), m_pOutputEA->getCurrPop().ldt(), m_pOutputWeights.dp(), 1));

  //Compute residual error in the network
  computeNetworkError();

  //Record current ouput weights
  cuda_call(cudaMemcpy(m_pOutputWeightsHistory.hr(m_iHiddenLayers), m_pOutputWeights.dp(), cacheWidth()*m_pOutputWeights.cellSize(), cudaMemcpyDeviceToHost));

  //Record current error
  m_pErrorHistory(0, m_iHiddenLayers) = m_dTrainingError;
  m_pErrorHistory(1, m_iHiddenLayers) = m_dTestError;

  //Detect overfitting
  if (m_iHiddenLayers == 0 || m_dTestError < m_pErrorHistory(1, m_iOptimalLayerCount)) {
    m_iOptimalLayerCount = m_iHiddenLayers;
    m_iOverfitCounter = 0;
  }
  else
    m_iOverfitCounter++;

  //Display progress
  if (!m_bSilentMode)
    std::cout << "Hidden layers:    " << std::setw(6) << m_iHiddenLayers << "    Training MSE:    " << std::setw(10) << m_dTrainingError << "    Test MSE:    " << std::setw(10) << m_dTestError << std::endl;
}

void CascadeCorrelation::computeNetworkError()
{
  //Compute error on training set
  computeTrainingError();

  //Compute error on test set
  computeTestError();
}

void CascadeCorrelation::computeTrainingError()
{
  //Copy results of top individual to host
  cuda_call(cudaMemcpy2D(m_pTrainingError.hp(), m_pTrainingError.cellSize(), m_pTrainingSimulation.dp() + topOutputInd(), m_pTrainingSimulation.pitch(), m_pTrainingSimulation.cellSize(), m_iTrainingSamples, cudaMemcpyDeviceToHost));

  real* err = m_pTrainingError.hp();
  real* tar = m_pTrainingTargets.hp();

  m_dTrainingError = 0.0;

  //Compute the error for each sample as well as the MSE
  for (nat i = 0; i < m_iTrainingSamples; i++) {
    real dif = err[i] - tar[i];
    err[i] = dif;
    m_dTrainingError += dif*dif;
  }

  m_dTrainingError /= m_iTrainingSamples;

  //Copy back error vector to device for next hidden unit simulation
  m_pTrainingError.syncDeviceWithHost();
}

void CascadeCorrelation::computeTestError()
{
  if (m_iTestSamples == 0)
    return;

  real alpha = 1.0, beta = 0.0;

  //Simulate results of top individual for all test samples
  cublas_call(CUBLAS_GEMV(m_pHandle, CUBLAS_OP_T, cacheWidth(), m_iTestSamples, &alpha, m_pTestCache.dp(), m_pTestCache.ldt(), m_pOutputWeights.dp(), 1, &beta, m_pTestSimulation.dp(), 1));

  //Copy results to host
  m_pTestSimulation.syncHostWithDevice();

  real* err = m_pTestSimulation.hp();
  real* tar = m_pTestTargets.hp();

  m_dTestError = 0.0;

  //Compute the error for each sample as well as the MSE
  for (nat i = 0; i < m_iTestSamples; i++) {
    real dif = err[i] - tar[i];
    m_dTestError += dif*dif;
  }

  m_dTestError /= m_iTestSamples;
}

__global__ void sim_unit_K(real* sim, nat simLdt, real* err, real* fit, nat* board, nat samples, nat totalInd, nat popSize)
{
  nat x = threadIdx.x + blockIdx.x*blockDim.x;

  if (x < totalInd) {
    board[x] = x & (popSize - 1);

    sim += x;

    real* end = sim + samples*simLdt;

    real ve_sum = 0.0;
    real v_sum = 0.0;
    real e_sum = 0.0;

    while (sim < end) {
      real v = ACTFUN(*sim);
      real e = *err;

      ve_sum += v*e;
      v_sum += v;
      e_sum += e;

      sim += simLdt;
      err++;
    }

    fit[x] = FABS(ve_sum - v_sum*(e_sum / samples));
  }
}

void CascadeCorrelation::simUnit(EvolutionaryAlgorithm* ea)
{
  real alpha = 1.0, beta = 0.0;

  //Simulate individual results for all training samples
  cublas_call(CUBLAS_GEMM(m_pHandle, CUBLAS_OP_N, CUBLAS_OP_N,
    ea->totalInd(), m_iTrainingSamples, cacheWidth(),
    &alpha,
    ea->getCurrPop().dp(), ea->getCurrPop().ldt(),
    m_pTrainingCache.dp(), m_pTrainingCache.ldt(),
    &beta,
    m_pTrainingSimulation.dp(), m_pTrainingSimulation.ldt()));

  //Compute individual fitness
  AsyncKernelLaunch(sim_unit_K, ea->blockCount(), ea->threadsPerBlock(),
    m_pTrainingSimulation.dp(), m_pTrainingSimulation.ldt(), m_pTrainingError.dp(), ea->getCurrFit().dp(), ea->getCurrBoard().dp(), m_iTrainingSamples, ea->totalInd(), ea->popSize());
}

void CascadeCorrelation::addHiddenUnit()
{
  //Hatch a new unit population
  m_pUnitEA->hatchSpecies(cacheWidth());

  //Evolve individuals
  m_pUnitEA->evolveSpecies();

  //Find top performer
  nat topPop;
  nat topInd;
  real topFit;

  m_pOutputEA->topPerf(topPop, topInd, topFit);

  nat topUnitInd = m_pUnitEA->popSize()*topPop + topInd;

  //Copy top individual to new unit weights
  cublas_call(CUBLAS_COPY(m_pHandle, cacheWidth(), m_pUnitEA->getCurrPop().dp() + topUnitInd, m_pUnitEA->getCurrPop().ldt(), m_pUnitWeights.dr(m_iHiddenLayers), 1));

  //Simulate results of top individual for all test samples
  real alpha = 1.0, beta = 0.0;

  cublas_call(CUBLAS_GEMV(m_pHandle, CUBLAS_OP_T, cacheWidth(), m_iTestSamples, &alpha, m_pTestCache.dp(), m_pTestCache.ldt(), m_pUnitWeights.dr(m_iHiddenLayers), 1, &beta, m_pTestSimulation.dp(), 1));

  //Copy results of top individual to host
  cuda_call(cudaMemcpy2D(m_pTrainingError.hp(), m_pTrainingError.cellSize(), m_pTrainingSimulation.dp() + topUnitInd, m_pTrainingSimulation.pitch(), m_pTrainingSimulation.cellSize(), m_iTrainingSamples, cudaMemcpyDeviceToHost));
  
  m_pTestSimulation.syncHostWithDevice();

  //Compute unit output for all samples
  real* trainingOutput = m_pTrainingError.hp();
  real* testOutput = m_pTestSimulation.hp();

  for (nat i = 0; i < m_iTrainingSamples; i++)
    trainingOutput[i] = ACTFUN(trainingOutput[i]);

  for (nat i = 0; i < m_iTestSamples; i++)
    testOutput[i] = ACTFUN(testOutput[i]);

  //Copy back data to device cache
  cuda_call(cudaMemcpy2D(m_pTrainingCache.dp() + cacheWidth(), m_pTrainingCache.pitch(), m_pTrainingError.hp(), m_pTrainingError.cellSize(), m_pTrainingError.cellSize(), m_iTrainingSamples, cudaMemcpyHostToDevice));

  cuda_call(cudaMemcpy2D(m_pTestCache.dp() + cacheWidth(), m_pTestCache.pitch(), m_pTestSimulation.hp(), m_pTestSimulation.cellSize(), m_pTestSimulation.cellSize(), m_iTestSamples, cudaMemcpyHostToDevice));

  //Increase hidden layers count
  m_iHiddenLayers++;
}

EvolutionaryAlgorithm* CascadeCorrelation::makeEA(nat popCount, nat popSize, EAID id)
{
  EvolutionaryAlgorithm* ea = NULL;

  switch (id) {
  case GeneticAlgorithmEAID:
    ea = new GeneticAlgorithm(popCount, popSize, cacheWidth(), maxCacheWidth());
    break;
  case CrossEntropyEAID:
    break;
  default:
    error("Snake, answer me. Snake? SNAAAAKE!");
  }

  return ea;
}

void CascadeCorrelation::writeNetwork(char* path) const
{
  if (!m_bNetworkTrained)
    error("Cannot write network: not trained");

  File f(path, std::ios_base::out);

  //Write hidden layer weights
  for (nat i = 0; i < m_iOptimalLayerCount; i++) {
    for (nat j = 0; j < 1 + m_iInputDim + i; j++)
      f << m_pUnitWeights(i, j) << " ";
    f << File::endl;
  }

  //Write output weights
  for (nat i = 0; i < 1 + m_iInputDim + m_iOptimalLayerCount; i++)
    f << m_pOutputWeightsHistory(m_iOptimalLayerCount, i) << " ";
  f << File::endl;

  f.close();
}

void CascadeCorrelation::testNetwork() const
{
  if (!m_bNetworkTrained)
    error("Cannot test network: not trained");

  std::cout << std::endl;

  //Compute absolute value of the targets' mean
  real mean = 0.0;
  
  for (nat x = 0; x < m_iTrainingSamples; x++)
    mean += FABS(m_pTrainingTargets(0, x));

  for (nat x = 0; x < m_iTestSamples; x++)
    mean += FABS(m_pTestTargets(0, x));

  mean = mean / m_iSamples;

  //Create output vector
  real* output = new real[1 + m_iInputDim + m_iOptimalLayerCount];

  //Iterate over test samples
  for (nat x = 0; x < m_iTestSamples; x++) {

    //Copy bias and inputs
    memcpy(output, m_pTestCache.hr(x), (1 + m_iInputDim)*m_pTestCache.cellSize());

    //Compute output of each hidden unit
    for (nat i = 0; i < m_iOptimalLayerCount; i++) {
      
      real u = 0.0;

      for (nat j = 0; j < 1 + m_iInputDim + i; j++)
        u += output[j] * m_pUnitWeights(i, j);

      output[1 + m_iInputDim + i] = ACTFUN(u);
    }

    //Compute network output
    real out = 0.0;

    for (nat i = 0; i < 1 + m_iInputDim + m_iOptimalLayerCount; i++)
      out += output[i] * m_pOutputWeightsHistory(m_iOptimalLayerCount, i);

    //Display results
    set_color(ODD(x) ? DARKGRAY : GRAY);
    std::cout << "Testing sample #  ";

    set_color(ODD(x) ? DARKPINK : PINK);
    std::cout << std::setw(6) << m_pTestSetIDs(0, x);

    set_color(ODD(x) ? DARKGRAY : GRAY);
    std::cout << "    Network output =    ";

    set_color(ODD(x) ? DARKYELLOW : YELLOW);
    std::cout << std::setw(10) << out;

    set_color(ODD(x) ? DARKGRAY : GRAY);
    std::cout << "    Target value =    ";

    set_color(ODD(x) ? DARKYELLOW : YELLOW);
    std::cout << std::setw(10) << m_pTestTargets(0, x);

    set_color(ODD(x) ? DARKGRAY : GRAY);
    std::cout << "    Accuracy =    ";

    nat accuracy = static_cast<nat>(100.0 * (1.0 - MIN(mean, FABS(out - m_pTestTargets(0, x))) / mean));

    if (accuracy < 60)
      set_color(ODD(x) ? DARKRED : RED);
    else
      set_color(ODD(x) ? DARKGREEN : GREEN);

    std::cout << std::setw(4) << accuracy << "%";
    std::cout << std::endl;
  }

  set_color(GRAY);

  delete[] output;
}

void CascadeCorrelation::topPerf(nat& hiddenLayers, real& trainingError, real& testError)
{
  if (!m_bNetworkTrained)
    error("Cannot get top performance record: network not trained");

  hiddenLayers = m_iOptimalLayerCount;
  trainingError = m_pErrorHistory(0, m_iOptimalLayerCount);
  testError = m_pErrorHistory(1, m_iOptimalLayerCount);
}
