/*
# Copyright (C) 2005 Andrew E Firth, University of Otago, Dunedin, 
# New Zealand, aef(at)sanger.otago.ac.nz
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License (version 2) as
# published by the Free Software Foundation.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
# 02110-1301, USA.
*/

#include <iostream>
#include <cmath>
#include <fstream>
#include <cstdlib>
#include <iomanip>
#include <cstring>
using namespace std;
#define maxnseqs 1000

int main()
{
  int i, k, nslines, ndlines;
  double obs[6], mcs[6][maxnseqs], mcd[6][maxnseqs], num[6], sum[6], ssq[6];
  double mean_s[6], sigma_s[6], mean_d[6], sigma_d[6];
  double zs, zd, prob[6], pi, lpr[6];

  pi = 3.14159265;

  // Read in observations data
  ifstream obsfile("calcprob.obs.dat"); 
  if (!obsfile) {
    cerr << "Aborting: can't find file 'calcprob.obs.dat'.\n";
    exit(EXIT_FAILURE);
  }
  obsfile >> obs[0];
  obsfile >> obs[1];
  obsfile >> obs[2];
  obsfile >> obs[3];
  obsfile >> obs[4];
  obsfile >> obs[5];
  obsfile.close();

  // Read in simulations data (single)
  ifstream mcsfile("calcprob.mc.0.dat"); 
  if (!mcsfile) {
    cerr << "Aborting: can't find file 'calcprob.mc.0.dat'.\n";
    exit(EXIT_FAILURE);
  }
  nslines = -1;
  while (mcsfile.ignore(1000, '\n')) {
    ++nslines;
  }
  mcsfile.clear();
  mcsfile.seekg(0);
  for (k = 0; k < nslines; ++k) {
    mcsfile >> mcs[0][k];
    mcsfile >> mcs[1][k];
    mcsfile >> mcs[2][k];
    mcsfile >> mcs[3][k];
    mcsfile >> mcs[4][k];
    mcsfile >> mcs[5][k];
    mcsfile.ignore(1000, '\n');
  }
  mcsfile.close();

  // Read in simulations data (double)
  ifstream mcdfile("calcprob.mc.1.dat"); 
  if (!mcdfile) {
    cerr << "Aborting: can't find file 'calcprob.mc.1.dat'.\n";
    exit(EXIT_FAILURE);
  }
  ndlines = -1;
  while (mcdfile.ignore(1000, '\n')) {
    ++ndlines;
  }
  mcdfile.clear();
  mcdfile.seekg(0);
  for (k = 0; k < ndlines; ++k) {
    mcdfile >> mcd[0][k];
    mcdfile >> mcd[1][k];
    mcdfile >> mcd[2][k];
    mcdfile >> mcd[3][k];
    mcdfile >> mcd[4][k];
    mcdfile >> mcd[5][k];
    mcdfile.ignore(1000, '\n');
  }
  mcdfile.close();

  /*
  // Transform data (N123, NsynNnon, but not MLOGD)
  for (i = 1; i < 6; ++i) {
    obs[i] = asin(sqrt(obs[i]));
    for (k = 0; k < nslines; ++k) {
      mcs[i][k] = asin(sqrt(mcs[i][k]));
    }
    for (k = 0; k < ndlines; ++k) {
      mcd[i][k] = asin(sqrt(mcd[i][k]));
    }
  }
  */

  // Calculate mean and standard deviation for simulations
  for (i = 0; i < 6; ++i) {
    num[i] = sum[i] = ssq[i] = 0.;
    for (k = 0; k < nslines; ++k) {
      num[i] += 1.;
      sum[i] += mcs[i][k];
      ssq[i] += mcs[i][k]*mcs[i][k];
    }
    if (num[i] > 0.) {
      mean_s[i] = sum[i]/num[i];
    } else {
      mean_s[i] = 0.;
    }
    if (num[i] > 1.) {
      sigma_s[i] = sqrt((ssq[i]-sum[i]*sum[i]/num[i])/(num[i]-1.));
    } else {
      sigma_s[i] = 0.;
    }
  }
  for (i = 0; i < 6; ++i) {
    num[i] = sum[i] = ssq[i] = 0.;
    for (k = 0; k < ndlines; ++k) {
      num[i] += 1.;
      sum[i] += mcd[i][k];
      ssq[i] += mcd[i][k]*mcd[i][k];
    }
    if (num[i] > 0.) {
      mean_d[i] = sum[i]/num[i];
    } else {
      mean_d[i] = 0.;
    }
    if (num[i] > 1.) {
      sigma_d[i] = sqrt((ssq[i]-sum[i]*sum[i]/num[i])/(num[i]-1.));
    } else {
      sigma_d[i] = 0.;
    }
  }

  // Calculate z = (x - mu) / sigma, and log(P(z|d)/P(z|s))
  for (i = 0; i < 6; ++i) {
    if (sigma_s[i] <= 0. || sigma_d[i] <= 0.) {
      lpr[i] = 99999.;
    } else {
      if (sigma_s[i] > 0. && sigma_d[i] > 0.) {
	zs = (obs[i] - mean_s[i]) / sigma_s[i];
	zd = (obs[i] - mean_d[i]) / sigma_d[i];
	lpr[i] = 0.5*((zs*zs)-(zd*zd)) + log(sigma_s[i]/sigma_d[i]);
      } else {
	zs = zd = 0.;  // actually equals Inf unless obs[i] = mean_[sd][i]
	lpr[i] = 0.;   // actually +/- Inf, but at least if you return 0 it
	// tells the user that the results are unusable
      } 
    }
  }

  // Output log(P(z|d)/P(z|s)) values
  ofstream outfile("calcprob.out.dat"); 
  if (!outfile) {
    cerr << "Aborting: can't open file 'calcprob.out.dat'.\n";
    exit(EXIT_FAILURE);
  }
  for (i = 0; i < 6; ++i) {
    outfile << lpr[i] << " ";
  }
  outfile << "\n";

  return(0);
}