#ifndef _IMPULSE_H_
#define _IMPULSE_H_

#include <cstring>
#include <string>
#include <vector>
#include <map>
using namespace std;

// General
#define N_AMODEL_PARMS  (6) // (h0,h1,h2,t1,t2,b)    
#define N_BMODEL_PARMS  (4) // (h0,h1,t1,b)	    
#define N_CMODEL_PARMS  (7) // (h0,h1,h2,t1,t2,b1,b2)
#define N_DMODEL_PARMS (10) // (h0,h1,h2,h3,h4,t1,t2,t3,t4,b);
#define N_EMODEL_PARMS  (6) // (h0,h1,h2,t1,t2,b)    
#define N_FMODEL_PARMS  (7) // (h0,h1,h2,t1,t2,b1,b2)
#define N_GMODEL_PARMS N_AMODEL_PARMS
#define N_HMODEL_PARMS N_CMODEL_PARMS
#define N_IMODEL_PARMS N_AMODEL_PARMS
#define N_JMODEL_PARMS N_BMODEL_PARMS

#define EXPRESS_DIR    "/Data/Expression/Yeast/"
#define BASE_DIR       "/Projects/TimeCourses/Runs/"
#define AUX_DATA_DIR   "/Projects/TimeCourses/Runs/SingleGenesFit/AuxData/"
#define GENES_FILENAME  "genes.lst"

#define MODEL_GH_NUM_EXTRA_ZEROS (3)

#define N_REPEATS (2500)


// Handle priors
float pdf_gauss(float *x, float *m,float *c,int n,int mode);
float deriv_pdf_gauss(float *x, int i, float *m,float *s, int n,int mode);
float log_pdf_gauss(float *x, float *m,float *s,int n,int mode);
float calc_emp_prior(float *x, float **m, float **s, float *weights,
		 int n_clusters, int dim, int mode);
float calc_hand_prior(float *x, int dim);

void randperm (int *perm, int n);
void print_theta(float*theta,char model);
void sigmoid (float *local_xs, float *local_ys, int n,
              float x_shift, float slope, float y_scale, float y_shift);
void calc_ys_Amodel(float *local_xs, float* local_ys, float *theta, int local_n_points);
void calc_ys_Bmodel(float *local_xs, float* local_ys, float *theta, int local_n_points);
void calc_ys_Cmodel(float *local_xs, float* local_ys, float *theta, int local_n_points);
void calc_ys_Dmodel(float *local_xs, float* local_ys, float *theta, int local_n_points);
void calc_ys_Emodel(float *local_xs, float* local_ys, float *theta, int local_n_points);
void calc_ys_Fmodel(float *local_xs, float* local_ys, float *theta, int local_n_points);

float deriv_pdf_gauss(float *x, int i_x, float *means,float *covars, int n,int mode);
float log_pdf_gauss(float *x, float *means,float *covars,int n,int mode);
float pdf_gauss(float *x, float *means,float *covars,int n,int mode);
float calc_emp_prior(float *x, float **means, float **covar, float *weights,
                     int n_clusters, int dim, int mode);
float calc_hand_prior(float *x, int n);
float calc_deriv_hand_prior(float *x, int n);
float calc_deriv_model_log_emp_prior(float *x, int i_x, float **means,
				     float **covar, float *weights, int n_clusters,
				     int dim, int mode);

void calc_deriv_Amodel(float *theta, float *d_thetas);
void calc_deriv_Bmodel(float *theta, float *d_thetas);
void calc_deriv_Cmodel(float *theta, float *d_thetas);
void calc_deriv_Dmodel(float *theta, float *d_thetas);
void calc_deriv_Emodel(float *theta, float *d_thetas);
void calc_deriv_Fmodel(float *theta, float *d_thetas);

float calc_error_Amodel(float*theta);
float calc_error_Bmodel(float*theta);
float calc_error_Cmodel(float*theta);
float calc_error_Dmodel(float*theta);
float calc_error_Emodel(float*theta);
float calc_error_Fmodel(float*theta);

double cg_fit_single_nr(char model, float *theta, int zero_transform, int n_repeats);
double cg_fit_single_nr_Gmodel(float *best_theta, int zero_transform, int n_repeats);
double cg_fit_single_nr_Bmodel(float *best_theta,int zero_transform,int n_repeats);
double cg_fit_single_nr_Bmodel(float *best_theta,int zero_transform,int n_repeats,int n_extra_zeros);
double cg_fit_single_nr_Cmodel(float *best_theta,int zero_transform,int n_repeats);
double cg_fit_single_nr_Dmodel(float *best_theta,int zero_transform,int n_repeats);
double cg_fit_single_nr_Emodel(float *best_theta,int zero_transform,int n_repeats);

int count_valid(bool ** mask, float ** express, int i_gene, vector<int> segs_inds);
int set_n_model_parms(char model);
void set_xs_ys(float ** express, bool  ** mask, int i_gene, vector<float> segs_lags,
               vector<int> segs_inds, int zero_transform, int num_extra_zeros);
void alloc_xs_ys(int size);
double estimate_model(float * theta, float** express, bool  ** mask, int i_gene,
                      vector<float> segs_lags, vector<int> segs_inds, ofstream& outstream,
                      int zero_transform, char model);
void estimate_pvals(float *  theta, float ** express, bool  ** mask, vector<float> segs_lags,
		      vector<int> segs_inds, int i_gene, int n_shuffles, int n_valid,
		      double err, ofstream &outstream, int zero_transform, char model);
void estimate_stability(float *  theta, float ** express, bool  ** mask, vector<float> segs_lags,
                        vector<int> segs_inds, int i_gene, int n_shuffles, int n_valid, double err,
                        ofstream &outstream, int zero_transform, char model);
void estimate_cvall(float * theta, float** express, bool  ** mask, int i_gene,
                    vector<float> segs_lags, vector<int> segs_inds, int n_valid,
                    ofstream &outstream, int zero_transform, char model);
void concat_to_float_mat(float **&mat,float **mat2,int n_lines, int n_fields, int n_fields2);
void concat_to_bool_mat(bool **&mat,bool **mat2,int n_lines, int n_fields, int n_fields2);
void recalc_data(float ** & express, bool ** & mask, int & n_genes, int & n_fields);



#endif /* _IMPULSE_H_ */