$$\phantom{******** If you see this on the webpage then the browser could not locate *********}$$
$$\phantom{******** jsMath/easy/load.js or the variable root is set wrong in this file *********}$$
$$\newcommand{\vector}[1]{\left[\begin{array}{c} #1 \end{array}\right]}$$ $$\newenvironment{matrix}{\left[\begin{array}{cccccccccc}} {\end{array}\right]}$$ $$\newcommand{\A}{{\cal A}}$$ $$\newcommand{\W}{{\cal W}}$$

/* Copyright 2008-2014 Research Foundation State University of New York   */

/* This file is part of QUB Express.                                      */

/* QUB Express is free software; you can redistribute it and/or modify    */
/* it under the terms of the GNU General Public License as published by   */
/* the Free Software Foundation, either version 3 of the License, or      */
/* (at your option) any later version.                                    */

/* QUB Express 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,     */
/* named LICENSE.txt, in the QUB Express program directory.  If not, see  */
/* <http://www.gnu.org/licenses/>.                                        */

/*

*/

#ifdef _WIN32
#include <windows.h>
#else
#include <stdlib.h>
#define BOOL int
#define TRUE 1
#define FALSE 0
#endif

#include <string.h>
#include <iostream>
#include <fstream>

#include <string.h>
#include <set>

#include <boost/math/special_functions/next.hpp>
#define TDEAD_DELTA_FLOATS 5

#include "max_inter_ll.h"
#include "../qubfast/max_ll_util.h"
#include "../qubfast/ublas_matrixutil.h"

using std::cout;
using std::cerr;
using std::endl;
using std::set;


extern "C" MAXILL_API int inter_ll_arr(
						 int Ns, int *clazz, double *P0_in,  double *K0_in, double *K1_in, double *K2_in,
						 int *Ligand_in, int *Voltage_in, int *Pressure_in,
						 double *Constants, double tdead,
						 int Nseg, int *dwellCounts, int **classeses, float **durationses,
						 double *ll)
{
  std::vector<double*> K0, K1, K2;
  std::vector<int*> Lig, Vol, Pres;
  for (int i=0; i<Ns; ++i) {
    K0.push_back( K0_in + i*Ns );
    K1.push_back( K1_in + i*Ns );
    K2.push_back( K2_in + i*Ns );
    Lig.push_back( Ligand_in + i*Ns );
    Vol.push_back( Voltage_in + i*Ns );
    Pres.push_back( Pressure_in + i*Ns );
  }
  
  return inter_ll(Ns, clazz, P0_in, &(K0[0]), &(K1[0]), &(K2[0]), &(Lig[0]), &(Vol[0]), &(Pres[0]),
		  Constants, tdead, Nseg, dwellCounts, classeses, durationses, ll);
}


extern "C" MAXILL_API int inter_ll(
					     int Ns, int *clazz, double *P0_in,  double **K0_in, double **K1_in, double **K2_in,
					     int **Ligand_in, int **Voltage_in, int **Pressure_in,
					     double *Constants, double tdead,
					     int Nseg, int *dwellCounts, int **classeses, float **durationses,
					     double *ll)
{
  //std::ofstream logger;
  //logger.open("s:\\maxill_log.txt");
  
  matrix<double> K0;         matrix_from_arrays(K0, K0_in, Ns, Ns);
  matrix<double> K1;         matrix_from_arrays(K1, K1_in, Ns, Ns);
  matrix<double> K2;         matrix_from_arrays(K2, K1_in, Ns, Ns);
  matrix<int>    Ligand;     matrix_from_arrays(Ligand, Ligand_in, Ns, Ns);
  matrix<int>    Voltage;    matrix_from_arrays(Voltage, Voltage_in, Ns, Ns);
  matrix<int>    Pressure;    matrix_from_arrays(Pressure, Voltage_in, Ns, Ns);
  
  matrix<double> Q;          BuildQ_p(Ns, Q, K0, K1, K2, Ligand, Voltage, Pressure, Constants);
  matrix<double> eQ;         QtoQe(Q, clazz, tdead, eQ);
  matrix<double> P0;
  if ( P0_in ) 
    {
      matrix_from_array(P0, P0_in, 1, Ns);
      double sum = 0.0;
      for (int i=0; i<Ns; ++i)
	sum += P0_in[i];
      for (int i=0; i<Ns; ++i)
	P0(0,i) = (sum > 1e-8 ) ? (P0(0,i) / sum) : (1.0 / Ns);
    }
  else
    QtoPe(eQ, P0);
  int            Na;
  
  //logger << "K0: " << K0 << endl;
  //logger << "K1: " << K1 << endl;
  //logger << "K2: " << K2 << endl;
  //logger << "L:  " << Ligand << endl;
  //logger << "V:  " << Voltage << endl;
  //logger << "P:  " << Pressure << endl;
  //logger << "Q:  " << Q << endl;
  //logger << "eQ: " << eQ << endl;
  //logger << "P0: " << P0 << endl;
  
  // max class
  int maxcls = clazz[0];
  for (int i=1; i<Ns; ++i)
    if ( clazz[i] > maxcls )
      maxcls = clazz[i];
  std::vector< std::vector<int> > clazz_states(maxcls+1);
  for (int i=0; i<Ns; ++i)
    clazz_states[clazz[i]].push_back(i);
  
  int maxcls_data = -1;
  for (int i=0; i<Nseg; ++i) {
    int *end = classeses[i] + dwellCounts[i];
    for (int *cl = classeses[i]; cl != end; ++cl) 
      if ( *cl > maxcls_data )
	maxcls_data = *cl;
  }
  if ( maxcls_data > maxcls ) {
    cerr << "Not enough classes in the model." << endl;
    return -8;
  }
  
  // we eigendecompose each submatrix eQaa for quick exponentiation in the inner loop:
  std::vector< std::vector<double> > eigvals_aa(maxcls+1);
  std::vector< matrix<double> > eigvecs_aa(maxcls+1);
  std::vector< matrix<double> > eiginvs_aa(maxcls+1);
  vector<double> eig_imag(Ns);
  for (int a=0; a<=maxcls; ++a) {
    Na = (int) clazz_states[a].size();
    if ( Na ) {
      int *aa = &(clazz_states[a][0]);
      matrix<double> eQaa;   matrix_to_partition(eQ, Na, Na, aa, aa, eQaa);
      eigvals_aa[a].resize(Na);
      eigvecs_aa[a].resize(Na, Na);
      eiginvs_aa[a].resize(Na, Na);
      if ( ! eigen_vals_vecs(eQaa, &(eigvals_aa[a][0]), &(eig_imag[0]), eigvecs_aa[a], eiginvs_aa[a], cerr) ) {
	//output << "No convergence in eigen calculations" << endl;
	return -2;
      }
      for (int i=0; i<Na; ++i) {
	if ( fabs(eig_imag[i]) > 1e-5 ) {
	  //output << "Complex eigenvalues found: " << eig_imag[i] << endl;
	  return -2;
	}
	if ( eigvals_aa[a][i] > 1e-5 ) {
	  //output << "Positive eigenvalues found: " << eigvals_aa[a][i] << endl;
	  return -2;
	}
      }
    }
  }
  // Then for efficiency, allocate space for expm_eigen temps
  std::vector< matrix<double> > tmp_nxn(Ns+1), tmp_nxn_diag(Ns+1);
  for (int i=1; i<=Ns; ++i) {
    tmp_nxn[i].resize(i, i);
    tmp_nxn_diag[i].resize(i, i);
    for (int j=0; j<i; ++j)
      for (int k=0; k<i; ++k)
	tmp_nxn_diag[i](j,k) = 0.0;
  }
  // and a place for the final exponentiated matrices:
  std::vector< matrix<double> > eQ_aa(maxcls+1);
  for (int a=0; a<=maxcls; ++a) {
    Na = (int) clazz_states[a].size();
    if ( Na )
      eQ_aa[a].resize(Na, Na);
  }
  
  // partitions eQab are used for transition probability (per second)
  std::vector< std::vector< matrix<double> > > eQab;
  eQab.resize(maxcls+1);
  for (int a=0; a<=maxcls; ++a) {
    eQab[a].resize(maxcls+1);
    Na = (int) clazz_states[a].size();
    if ( Na ) {
      for (int b=0; b<=maxcls; ++b) {
	int Nb = (int) clazz_states[b].size();
	if ( (b != a) && Nb ) {
	  matrix<double>& subQ = eQab[a][b];
	  subQ.resize(Na, Nb);
	  matrix_to_partition(eQ, Na, Nb, &(clazz_states[a][0]), &(clazz_states[b][0]), subQ);
	}
      }
    }
  }
  
  // alpha (temp) rows per class, since they're different sizes
  // actually a pair apiece to avoid temporaries alpha_a <- alpha_a * At
  std::vector< vector<double> > alpha_pre(maxcls+1), alpha_post(maxcls+1);
  for (int a=0; a<=maxcls; ++a) {
    alpha_pre[a].resize( clazz_states[a].size() ); // 1, clazz...size
    alpha_post[a].resize( clazz_states[a].size() ); // 1, clazz...size
  }
  
  // scale accumulator
  *ll = 0.0;
  
  for (int iseg = 0; iseg<Nseg; ++iseg) {
    int Nd = dwellCounts[iseg];
    int *classes = classeses[iseg];
    float *durations = durationses[iseg];
    int a = classes[0];
    Na = (int) clazz_states[a].size();
    int b = a;
    int Nb = Na;
    
    // normalized P0 among possible start states
    double tot = 0.0;
    for (int i=0; i<Na; ++i) {
      alpha_pre[a](i) = P0(0,clazz_states[a][i]);
      tot += alpha_pre[a](i);
    }
    if ( tot ) {
      for (int i=0; i<Na; ++i)
	alpha_pre[a](i) /= tot;
    } else {
      for (int i=0; i<Na; ++i)
	alpha_pre[a](i) = 1.0 / Na;
    }
    *ll += alpha_logscale(alpha_pre[a]);
    
    for (int d=1; d<Nd; ++d) {
      a = b;
      Na = Nb;
      b = classes[d];
      Nb = (int) clazz_states[b].size();
      
      // * exp(Qaa * dur) [scale?]   * eQab [scale]
      expm_eig(&(eigvals_aa[a][0]), eigvecs_aa[a], eiginvs_aa[a], durations[d-1],
	       tmp_nxn_diag[Na], tmp_nxn[Na], eQ_aa[a]);
      noalias(alpha_post[a]) = prod(alpha_pre[a], eQ_aa[a]);
      // *ll += alpha_logscale(alpha_post[a]);
      
      noalias(alpha_pre[b]) = prod(alpha_post[a], eQab[a][b]);
      *ll += alpha_logscale(alpha_pre[b]);
    }
    
    a = b;
    Na = Nb;
    expm_eig(&(eigvals_aa[a][0]), eigvecs_aa[a], eiginvs_aa[a], durations[Nd-1],
	     tmp_nxn_diag[Na], tmp_nxn[Na], eQ_aa[a]);
    noalias(alpha_post[a]) = prod(alpha_pre[a], eQ_aa[a]);
    // *ll += alpha_logscale(alpha_post[a]);
    
    std::vector<int> not_a;
    for (int i=0; i<Ns; ++i)
      if ( clazz[i] != a )
	not_a.push_back(i);
    int Nz = (int) not_a.size();
    matrix<double> eQaz(Na, Nz);
    matrix_to_partition(eQ, Na, Nz, &(clazz_states[a][0]), &(not_a[0]), eQaz);
    vector<double> alpha_z(Nz);
    noalias(alpha_z) = prod(alpha_post[a], eQaz);
    *ll += alpha_logscale(alpha_z);
  }
  
  return 0; // success
}




extern "C" MAXILL_API void mil_prep_events( int *ndwt, int *idwt, float *tdwt, double tdead_ms )
{
  float td = float(tdead_ms * 1.0e-3);
  int ndwells = *ndwt;
  int i; // outgoing index
  int j; // incoming index
  int cls;
  float tm;
  
  // first, drop initial too-short events
  for ( j=0; j<ndwells; j++ ) {
    tm = tdwt[j] * 1.0e-3f;
    if ( (boost::math::float_distance(tm, td) < TDEAD_DELTA_FLOATS) || (tm > td) )
      break;
  }
  // then add the time-values of other too-shorts to the prev. dwell
  for ( i=0; j<ndwells; j++ ) {
    cls = idwt[j];
    tm = tdwt[j] * 1.0e-3f;
    
    if ( (i > 0) && (idwt[i-1] == cls) )
      tdwt[i-1] = (float)(tdwt[i-1] + tm);
    else if ( (i > 0) && ((boost::math::float_distance(tm, td) >= TDEAD_DELTA_FLOATS) && (tm <= td)) )
      tdwt[i-1] = (float)(tdwt[i-1] + tm);
    else {
      idwt[i] = cls;
      tdwt[i] = (float) (tm - td);
      i++;
    }
  }
  *ndwt = i;
}

extern "C" MAXILL_API void mil_prep_segments( int Nseg, int *dwellCounts, int **classeses, float **durationses, double tdead_ms )
{
  for (int i=0; i<Nseg; ++i)
    mil_prep_events(dwellCounts+i, classeses[i], durationses[i], tdead_ms);
}