@@ -5,3 +5,4 @@

 *.swp

 *.swo

 *~

+*_

@@ -1,6 +1,6 @@

-Package: anafinder

+Package: aneufinder

 Type: Package

-Title: Anaploidy analysis of Seq-data using a Hidden Markov Model

+Title: Aneuploidy analysis of Seq-data using a Hidden Markov Model

 Version: 0.1

 Date: 2014-05-01

 Author: Aaron Taudt, David Porubsky

@@ -1,11 +1,11 @@

-useDynLib(chromstar)

+useDynLib(aneufinder)

 # Export all names

 # exportPattern(".")

 export(

 bed2bin,

 bam2bin,

-bedGraph2bin,

+bedGraph2bin

 )

 # Import all packages listed as Imports or Depends

@@ -1,16 +1,14 @@

-#include "loghmm.h"

 #include "scalehmm.h"

-

 // ---------------------------------------------------------------

 // void R_univariate_hmm()

 // This function takes parameters from R, creates a univariate HMM object, creates the distributions, runs the Baum-Welch and returns the result to R.

 // ---------------------------------------------------------------

 extern "C" {

-void R_univariate_hmm(int* O, int* T, int* N, int* densityNames, double* r, double* p, int* maxiter, int* maxtime, double* eps, double* post, double* A, double* proba, double* loglik, double* softweights, double* initial_r, double* initial_p, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads) {

+void R_univariate_hmm(int* O, int* T, int* N, double* means, double* variances, int* maxiter, int* maxtime, double* eps, double* post, double* A, double* proba, double* loglik, double* weights, double* initial_means, double* initial_variances, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads) {

 	// Define logging level

-// 	FILE* pFile = fopen("chromStar.log", "w");

+// 	FILE* pFile = fopen("aneufinder.log", "w");

 // 	Output2FILE::Stream() = pFile;

 //  	FILELog::ReportingLevel() = FILELog::FromString("INFO");

  	FILELog::ReportingLevel() = FILELog::FromString("ITERATION");

@@ -42,11 +40,10 @@ void R_univariate_hmm(int* O, int* T, int* N, int* densityNames, double* r, doub

 	// Create the HMM

 	FILE_LOG(logDEBUG1) << "Creating a univariate HMM";

-// 	LogHMM* model = new LogHMM(*T, *N);

-	ScaleHMM* model = new ScaleHMM(*T, *N);

+	ScaleHMM* hmm = new ScaleHMM(*T, *N);

 	// Initialize the transition probabilities and proba

-	model->initialize_transition_probs(initial_A, *use_initial_params);

-	model->initialize_proba(initial_proba, *use_initial_params);

+	hmm->initialize_transition_probs(initial_A, *use_initial_params);

+	hmm->initialize_proba(initial_proba, *use_initial_params);

 	// Calculate mean and variance of data

 	double mean = 0, variance = 0;

@@ -62,386 +59,70 @@ void R_univariate_hmm(int* O, int* T, int* N, int* densityNames, double* r, doub

 	variance = variance / *T;

 	FILE_LOG(logINFO) << "data mean = " << mean << ", data variance = " << variance;		

-	// Go through all states of the model and assign the density functions

+	// Go through all states of the hmm and assign the density functions

 	srand (clock());

 	int rand1, rand2;

 	double imean, ivariance;

-	//int i_count = 0;

-	for (int i_state=0; i_state<model->N; i_state++)

+	for (int istate=0; istate<*N; istate++)

 	{

 		if (*use_initial_params) {

-			FILE_LOG(logINFO) << "Using given parameters for r and p";

-			imean = (1-initial_p[i_state])*initial_r[i_state] / initial_p[i_state];

-			ivariance = imean / initial_p[i_state];

+			FILE_LOG(logINFO) << "Using given parameters for mean and variance";

+// 			TODO

 			FILE_LOG(logDEBUG3) << "imean = " << imean;

 			FILE_LOG(logDEBUG3) << "ivariance = " << ivariance;

 		} else {

-

-// 			// Disturb mean and variance for use as randomized initial parameters

-// 			FILE_LOG(logINFO) << "Using random initialization for r and p";

-// 			rand1 = rand();

-// 			rand2 = rand();

-// 			imean = (double)rand1/(double)RAND_MAX * 10*mean;

-// // 			ivariance = imean + (double)rand2/(double)RAND_MAX * 10*variance; // variance has to be greater than mean, otherwise r will be negative

-// 			ivariance = imean + (double)rand2/(double)RAND_MAX * 20*imean; // variance has to be greater than mean, otherwise r will be negative

-// 			FILE_LOG(logDEBUG3) << "RAND_MAX = " << RAND_MAX;

-// 			FILE_LOG(logDEBUG3) << "rand1 = " << rand1;

-// 			FILE_LOG(logDEBUG3) << "rand2 = " << rand2;

-// 			FILE_LOG(logDEBUG3) << "imean = " << imean;

-// 			FILE_LOG(logDEBUG3) << "ivariance = " << ivariance;

-

-// 	 		// Empirical initialization

-// 	 		if (i_state == 1) {

-// 				FILE_LOG(logINFO) << "Initializing r and p empirically for state 1";

-// 	 			imean = mean/2;

-// 	 			ivariance = imean*2;

-// 	 		} else if (i_state == 2) {

-// 				FILE_LOG(logINFO) << "Initializing r and p empirically for state 2";

-// 	 			imean = mean*2;

-// 	 			ivariance = imean*2;

-// 	 		} 

-

-			// Simple initialization, seems to give the fastest convergence

-	 		if (i_state == 1) {

-				FILE_LOG(logINFO) << "Initializing r and p for state 1";

-	 			imean = mean;

-	 			ivariance = variance;

-	 		} else if (i_state == 2) {

-				FILE_LOG(logINFO) << "Initializing r and p for state 2";

-	 			imean = mean+1;

-	 			ivariance = variance;

-	 		} 

-			

-			// Calculate r and p from mean and variance

-			initial_r[i_state] = pow(imean,2)/(ivariance-imean);

-			initial_p[i_state] = imean/ivariance;

+			// Simple initialization

+			imean = mean * pow(2, istate-1);

+			ivariance = variance * pow(2, istate-1);

 		}

-		if (densityNames[i_state] == NB)

-		{

-			FILE_LOG(logDEBUG1) << "Using negative binomial for state " << i_state;

-			NegativeBinomial *d = new NegativeBinomial(O, model->T, initial_r[i_state], initial_p[i_state]); // delete is done inside ~LogHMM()

-			model->densityFunctions.push_back(d);

-		}

-		else if (densityNames[i_state] == ZI)

-		{

-			FILE_LOG(logDEBUG1) << "Using only zeros for state " << i_state;

-			OnlyZeros *d = new OnlyZeros(O, model->T); // delete is done inside ~LogHMM()

-			model->densityFunctions.push_back(d);

-		}

-		else

-		{

-			FILE_LOG(logWARNING) << "Density not specified, using default negative binomial for state " << i_state;

-			NegativeBinomial *d = new NegativeBinomial(O, model->T, initial_r[i_state], initial_p[i_state]);

-			model->densityFunctions.push_back(d);

-		}

+		FILE_LOG(logDEBUG1) << "Using gaussian normal for state " << istate;

+		Normal *d = new Normal(O, *T, imean, ivariance); // delete is done inside ~ScaleHMM()

+		hmm->densityFunctions.push_back(d);

+

 	}

 	// Do the Baum-Welch to estimate the parameters

 	FILE_LOG(logDEBUG1) << "Starting Baum-Welch estimation";

-	model->baumWelch(maxiter, maxtime, eps);

+	hmm->baumWelch(maxiter, maxtime, eps);

 	FILE_LOG(logDEBUG1) << "Finished with Baum-Welch estimation";

 	// Compute the posteriors and save results directly to the R pointer

-	double** posteriors = allocDoubleMatrix(model->N, model->T);

-	model->get_posteriors(posteriors);

+	double** posteriors = allocDoubleMatrix(*N, *T);

+	hmm->get_posteriors(posteriors);

 	FILE_LOG(logDEBUG1) << "Recode posteriors into column representation";

 	#pragma omp parallel for

-	for (int iN=0; iN<model->N; iN++)

+	for (int iN=0; iN<*N; iN++)

 	{

-		for (int t=0; t<model->T; t++)

+		for (int t=0; t<*T; t++)

 		{

-			post[t + iN * model->T] = posteriors[iN][t];

+			post[t + iN * (*T)] = posteriors[iN][t];

 		}

 	}

-	freeDoubleMatrix(posteriors, model->N);

+	freeDoubleMatrix(posteriors, *N);

 	FILE_LOG(logDEBUG1) << "Return parameters";

 	// also return the estimated transition matrix and the initial probs

-	for (int i=0; i<model->N; i++)

+	for (int i=0; i<*N; i++)

 	{

-		proba[i] = model->get_proba(i);

-		for (int j=0; j<model->N; j++)

+		proba[i] = hmm->get_proba(i);

+		for (int j=0; j<*N; j++)

 		{

-          A[i*model->N + j] = model->A[i][j];

+			A[i * (*N) + j] = hmm->get_A(i,j);

 		}

 	}

 	// copy the estimated distribution params

-	for (int i=0; i<model->N; i++)

+	for (int i=0; i<*N; i++)

 	{

-			if (model->densityFunctions[i]->getType() == NB) 

-			{

-					NegativeBinomial* d = (NegativeBinomial*)(model->densityFunctions[i]);

-					r[i] = d->getR();

-					p[i] = d->getP();

-			}

+		means[i] = hmm->densityFunctions[i]->get_mean();

+		variances[i] = hmm->densityFunctions[i]->get_variance();

 	}

-	*loglik = model->logP;

-	model->calc_weights(softweights);

+	*loglik = hmm->get_logP();

+	hmm->calc_weights(weights);

-	FILE_LOG(logDEBUG1) << "Deleting the model";

-	delete model;

+	FILE_LOG(logDEBUG1) << "Deleting the hmm";

+	delete hmm;

 }

 } // extern C

-// ---------------------------------------------------------------

-// void R_multivariate_hmm()

-// This function takes parameters from R, creates a multivariate HMM object, creates the distributions, runs the Baum-Welch and returns the result to R.

-// ---------------------------------------------------------------

-extern "C" {

-void R_multivariate_hmm(int* O, int* T, int* N, int *Nmod, int* states, double* r, double* p, double* w, double* cor_matrix_inv, double* det, int* maxiter, int* maxtime, double* eps, double* post, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads){

-

-	// Define logging level {"ERROR", "WARNING", "INFO", "ITERATION", "DEBUG", "DEBUG1", "DEBUG2", "DEBUG3", "DEBUG4"}

-// 	FILE* pFile = fopen("chromStar.log", "w");

-// 	Output2FILE::Stream() = pFile;

- 	FILELog::ReportingLevel() = FILELog::FromString("ITERATION");

-//  	FILELog::ReportingLevel() = FILELog::FromString("DEBUG");

-

-	// Parallelization settings

-	omp_set_num_threads(*num_threads);

-

-	// Print some information

-	FILE_LOG(logINFO) << "number of states = " << *N;

-	FILE_LOG(logINFO) << "number of bins = " << *T;

-	if (*maxiter < 0)

-	{

-		FILE_LOG(logINFO) << "maximum number of iterations = none";

-	} else {

-		FILE_LOG(logINFO) << "maximum number of iterations = " << *maxiter;

-	}

-	if (*maxtime < 0)

-	{

-		FILE_LOG(logINFO) << "maximum running time = none";

-	} else {

-		FILE_LOG(logINFO) << "maximum running time = " << *maxtime << " sec";

-	}

-	FILE_LOG(logINFO) << "epsilon = " << *eps;

-	FILE_LOG(logINFO) << "number of modifications = " << *Nmod;

-

-	// Recode the observation vector to matrix representation

-	clock_t clocktime = clock(), dtime;

-	int** multiO = allocIntMatrix(*Nmod, *T);

-	for (int imod=0; imod<*Nmod; imod++)

-	{

-		for (int t=0; t<*T; t++)

-		{

-			multiO[imod][t] = O[imod*(*T)+t];

-		}

-	}

-	dtime = clock() - clocktime;

-	FILE_LOG(logDEBUG1) << "recoding observation vector to matrix representation: " << dtime << " clicks";

-

-	// Create the HMM

-	FILE_LOG(logDEBUG1) << "Creating the multivariate HMM";

-// 	LogHMM* model = new LogHMM(*T, *N, *Nmod);

-	ScaleHMM* model = new ScaleHMM(*T, *N, *Nmod);

-	// Initialize the transition probabilities and proba

-	model->initialize_transition_probs(initial_A, *use_initial_params);

-	model->initialize_proba(initial_proba, *use_initial_params);

-	

-	// Print logproba and A

-// 	for (int iN=0; iN<*N; iN++)

-// 	{

-// 		FILE_LOG(logINFO) << "proba["<<iN<<"] = " <<exp(model->logproba[iN]);

-// 		for (int jN=0; jN<*N; jN++)

-// 		{

-// 			FILE_LOG(logINFO) << "A["<<iN<<"]["<<jN<<"] = " << model->A[iN][jN];

-// 		}

-// 	}

-

-	// Prepare the binary_states (univariate) vector: binary_states[N][Nmod], e.g., binary_states[iN][imod] tells me at state states[iN], modification imod is non-enriched (0) or enriched (1)

-	FILE_LOG(logDEBUG1) << "Preparing the binary_states vector";

-	bool **binary_states = allocBoolMatrix(model->N, model->Nmod);

-	for(int iN=0; iN < model->N; iN++) //for each comb state considered

-	{

-		for(int imod=0; imod < model->Nmod; imod++) //for each modification of this comb state

-		{

-			binary_states[iN][imod] = states[iN]&(int)pow(2,model->Nmod-imod-1);//if =0->hidden state states[iN] has modification imod non enriched; if !=0->enriched

-			if (binary_states[iN][imod] != 0 )

-				binary_states[iN][imod] = 1;

-		}

-	}

-

-	/* initialize the distributions */

-	FILE_LOG(logDEBUG1) << "Initializing the distributions";

-	for (int iN=0; iN<model->N; iN++) //for each combinatorial state

-	{

-		vector <Density*> tempMarginals;            

-		for (int imod=0; imod < model->Nmod; imod++) //for each modification

-		{

-			Density *d;

-			if (binary_states[iN][imod]) //construct the marginal density function for modification imod being enriched

-			{

-				d = new NegativeBinomial(multiO[imod], model->T, r[2*imod+1], p[2*imod+1]); // delete is done inside ~MVCopulaApproximation()

-			}

-			else //construct the density function for modification imod being non-enriched

-			{

-				d = new ZiNB(multiO[imod], model->T, r[2*imod], p[2*imod], w[imod]); // delete is done inside ~MVCopulaApproximation()

-			}

-			tempMarginals.push_back(d);

-		}

-		//MVCopulaApproximation *tempMVdens = new MVCopulaApproximation(O, tempMarginals, &(cor_matrix_inv[iN*Nmod*Nmod]), det[iN]);

-		FILE_LOG(logDEBUG1) << "Calling MVCopulaApproximation for state " << iN;

-		MVCopulaApproximation *tempMVdens = new MVCopulaApproximation(multiO, model->T, tempMarginals, &(cor_matrix_inv[iN*model->Nmod*model->Nmod]), det[iN]); // delete is done inside ~LogHMM()

-		model->densityFunctions.push_back(tempMVdens);

-	}

-	

-	// Estimate the parameters

-	FILE_LOG(logDEBUG1) << "Starting Baum-Welch estimation";

-	model->baumWelch(maxiter, maxtime, eps);

-	FILE_LOG(logDEBUG1) << "Finished with Baum-Welch estimation";

-	

-	// Compute the posteriors and save results directly to the R pointer

-	double** posteriors = allocDoubleMatrix(model->N, model->T);

-	model->get_posteriors(posteriors);

-	FILE_LOG(logDEBUG1) << "Recode posteriors into column representation";

-	for (int iN=0; iN<model->N; iN++)

-	{

-		for (int t=0; t<model->T; t++)

-		{

-			post[t + iN * model->T] = posteriors[iN][t];

-		}

-	}

-	freeDoubleMatrix(posteriors, model->N);

-	

-	FILE_LOG(logDEBUG1) << "Return parameters";

-	// also return the estimated transition matrix and the initial probs

-	for (int i=0; i<model->N; i++)

-	{

-		proba[i] = model->get_proba(i);

-		for (int j=0; j<model->N; j++)

-		{

-				A[i*model->N + j] = model->A[i][j];

-		}

-	}

-	*loglik = model->logP;

-

-	FILE_LOG(logDEBUG1) << "Deleting the model";

-	delete model;

-	freeIntMatrix(multiO, *Nmod);

-	freeBoolMatrix(binary_states, *N);

-}

-} // extern C

-

-// ---------------------------------------------------------------

-// void R_multivariate_hmm_productBernoulli()

-// This function takes parameters from R, creates a multivariate HMM object, creates the distributions, runs the Baum-Welch and returns the result to R.

-// ---------------------------------------------------------------

-extern "C" {//observation is now the posterior for being UNMODIFIED

-void R_multivariate_hmm_productBernoulli(double* O, int* T, int* N, int *Nmod, int* states, int* maxiter, int* maxtime, double* eps, double* post, double* A, double* proba, double* loglik, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads){

-

-	// Define logging level

-// 	FILELog::ReportingLevel() = FILELog::FromString("DEBUG3");

-	FILELog::ReportingLevel() = FILELog::FromString("ITERATION");

-

-	// Parallelization settings

-	omp_set_num_threads(*num_threads);

-

-	// Print some information

-	FILE_LOG(logINFO) << "number of states = " << *N;

-	FILE_LOG(logINFO) << "number of bins = " << *T;

-	if (*maxiter < 0)

-	{

-		FILE_LOG(logINFO) << "maximum number of iterations = none";

-	} else {

-		FILE_LOG(logINFO) << "maximum number of iterations = " << *maxiter;

-	}

-	if (*maxtime < 0)

-	{

-		FILE_LOG(logINFO) << "maximum running time = none";

-	} else {

-		FILE_LOG(logINFO) << "maximum running time = " << *maxtime << " sec";

-	}

-	FILE_LOG(logINFO) << "epsilon = " << *eps;

-	FILE_LOG(logINFO) << "number of modifications = " << *Nmod;

-

-	// Recode the observation vector to matrix representation

-	clock_t clocktime = clock(), dtime;

-	double** multiO = allocDoubleMatrix(*Nmod, *T);

-	for (int imod=0; imod<*Nmod; imod++)

-	{

-		for (int t=0; t<*T; t++)

-		{

-			multiO[imod][t] = O[imod*(*T)+t];

-		}

-	}

-	dtime = clock() - clocktime;

-	FILE_LOG(logDEBUG1) << "recoding observation and probability vectors to matrix representation: " << dtime << " clicks";

-

-	// Create the HMM

-	FILE_LOG(logDEBUG1) << "Creating the multivariate HMM";

-	LogHMM* model = new LogHMM(*T, *N, *Nmod);

-

-	// Initialize the transition probabilities and proba

-	model->initialize_transition_probs(initial_A, *use_initial_params);

-	model->initialize_proba(initial_proba, *use_initial_params);

-	

-	// Print logproba and A

-// 	for (int iN=0; iN<*N; iN++)

-// 	{

-// 		FILE_LOG(logINFO) << "proba["<<iN<<"] = " <<exp(model->logproba[iN]);

-// 		for (int jN=0; jN<*N; jN++)

-// 		{

-// 			FILE_LOG(logINFO) << "A["<<iN<<"]["<<jN<<"] = " << model->A[iN][jN];

-// 		}

-// 	}

-

-	// Prepare the binary_states (univariate) vector: binary_states[N][Nmod], e.g., binary_states[iN][imod] tells me at state states[iN], modification imod is non-enriched (0) or enriched (1)

-	FILE_LOG(logDEBUG1) << "Preparing the binary_states vector";

-	bool **binary_states = allocBoolMatrix(model->N, model->Nmod);

-	for(int iN=0; iN < model->N; iN++) //for each comb state considered

-	{

-		for(int imod=0; imod < model->Nmod; imod++) //for each modification of this comb state

-		{

-			binary_states[iN][imod] = states[iN]&(int)pow(2,model->Nmod-imod-1);//if =0->hidden state states[iN] has modification imod non enriched; if !=0->enriched

-			if (binary_states[iN][imod] != 0 )

-				binary_states[iN][imod] = 1;

-		}

-	}

-

-	/* initialize the distributions */

-	FILE_LOG(logDEBUG1) << "Initializing the distributions";

-	for (int iN=0; iN<model->N; iN++) //for each combinatorial state

-	{

-		FILE_LOG(logDEBUG1) << "Calling BernoulliProduct";

-		BernoulliProduct *tempBP = new BernoulliProduct(multiO, binary_states[iN], *T, *Nmod); // delete is done inside ~LogHMM()

-		model->densityFunctions.push_back(tempBP);

-	}

-

-	// Estimate the parameters

-	FILE_LOG(logDEBUG1) << "Starting Baum-Welch estimation";

-	model->baumWelch(maxiter, maxtime, eps);

-	FILE_LOG(logDEBUG1) << "Finished with Baum-Welch estimation";

-	

-    // Compute the posteriors and save results directly to the R pointer

-    double** posteriors = allocDoubleMatrix(model->N, model->T);

-    model->get_posteriors(posteriors);

-	FILE_LOG(logDEBUG1) << "Recode posteriors into column representation";

-	for (int iN=0; iN<model->N; iN++)

-	{

-		for (int t=0; t<model->T; t++)

-		{

-			post[t + iN * model->T] = posteriors[iN][t];

-		}

-	}

-    freeDoubleMatrix(posteriors, model->N);

-	

-	FILE_LOG(logDEBUG1) << "Return parameters";

-	// also return the estimated transition matrix and the initial probs

-	for (int i=0; i<model->N; i++)

-	{

-		proba[i] = model->get_proba(i);

-		for (int j=0; j<model->N; j++)

-		{

-				A[i*model->N + j] = model->A[i][j];

-		}

-	}

-	*loglik = model->logP;

-

-	FILE_LOG(logDEBUG1) << "Deleting the model";

-	delete model;

-	freeDoubleMatrix(multiO, *Nmod);

-	freeBoolMatrix(binary_states, *N);

-}

-} // extern C

@@ -1,1128 +1,86 @@

 #include "densities.h"

-#include "utility.h"

-// ------------------------------------------------------------

-// Zero inflated negative binomial

-// ------------------------------------------------------------

+// ============================================================

+// Normal (Gaussian) density

+// ============================================================

-ZiNB::ZiNB()

-{

-	this->lxfactorials = NULL;

-}

-

-ZiNB::ZiNB(int* observations, int T, double r, double p, double w)

+// Constructor and Destructor ---------------------------------

+Normal::Normal(int* observations, int T, double mean, double variance)

 {

 	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

-	this->O = observations;

+	this->obs = observations;

 	this->T = T;

-	this->p = p;

-	this->r = r;

-	this->w = w;

-	this->lxfactorials = NULL;

-	if (this->O != NULL)

-	{

-		this->max_O = intMax(observations, T);

-		this->lxfactorials = (double*) calloc(max_O+1, sizeof(double));

-		this->lxfactorials[0] = 0.0;	// Not necessary, already 0 because of calloc

-		this->lxfactorials[1] = 0.0;

-		for (int j=2; j<=max_O; j++)

-		{

-			this->lxfactorials[j] = this->lxfactorials[j-1] + log(j);

-		}

-	}

+	this->mean = mean;

+	this->var = variance;

+	this->sd = sqrt(variance);

 }

-ZiNB::~ZiNB()

+Normal::~Normal()

 {

 	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

-	if (this->O != NULL)

-	{

-		free(this->lxfactorials);

-	}

 }

-void ZiNB::calc_logdensities(double* logdens)

-{

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

-	double logp = log(this->p);

-	double log1minusp = log(1-this->p);

-	double lGammaR,lGammaRplusX,lxfactorial;

-	lGammaR=lgamma(this->r);

-	// Select strategy for computing gammas

-	if (this->max_O <= this->T)

-	{

-		FILE_LOG(logDEBUG2) << "Precomputing gammas in " << __func__ << " for every O[t], because max(O)<=T";

-		double lGammaRplusX[this->max_O+1];

-		for (int j=0; j<=this->max_O; j++)

-		{

-			lGammaRplusX[j] = lgamma(this->r + j);

-		}

-		for (int t=0; t<this->T; t++)

-		{

-			lxfactorial = this->lxfactorials[(int) this->O[t]];

-			if (O[t] == 0)

-			{

-				logdens[t] = log( this->w + (1-this->w) * exp( lGammaRplusX[(int) this->O[t]] - lGammaR - lxfactorial + this->r * logp + this->O[t] * log1minusp ) );

-			}

-			else

-			{

-				logdens[t] = log(1-this->w) + lGammaRplusX[(int) this->O[t]] - lGammaR - lxfactorial + this->r * logp + this->O[t] * log1minusp;

-			}

-			if (isnan(logdens[t]))

-			{

-				FILE_LOG(logWARNING) << __PRETTY_FUNCTION__;

-				FILE_LOG(logWARNING) << "logdens["<<t<<"] = "<< logdens[t];

-				exit(1);

-			}

-		}

-	}

-	else

-	{

-		FILE_LOG(logDEBUG2) << "Computing gammas in " << __func__ << " for every t, because max(O)>T";

-		for (int t=0; t<this->T; t++)

-		{

-			lGammaRplusX = lgamma(this->r + this->O[t]);

-			lxfactorial = this->lxfactorials[(int) this->O[t]];

-			if (O[t] == 0)

-			{

-				logdens[t] = log( this->w + (1-this->w) * exp( lGammaRplusX - lGammaR - lxfactorial + this->r * logp + this->O[t] * log1minusp ) );

-			}

-			else

-			{

-				logdens[t] = log(1-this->w) + lGammaRplusX - lGammaR - lxfactorial + this->r * logp + this->O[t] * log1minusp;

-			}

-			if (isnan(logdens[t]))

-			{

-				FILE_LOG(logWARNING) << __PRETTY_FUNCTION__;

-				FILE_LOG(logWARNING) << "logdens["<<t<<"] = "<< logdens[t];

-				exit(1);

-			}

-		}

-	}

-}

-

-void ZiNB::calc_densities(double* dens)

-{

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

-	double logp = log(this->p);

-	double log1minusp = log(1-this->p);

-	double lGammaR,lGammaRplusX,lxfactorial;

-	lGammaR=lgamma(this->r);

-	// Select strategy for computing gammas

-	if (this->max_O <= this->T)

-	{

-		FILE_LOG(logDEBUG2) << "Precomputing gammas in " << __func__ << " for every O[t], because max(O)<=T";

-		double lGammaRplusX[this->max_O+1];

-		for (int j=0; j<=this->max_O; j++)

-		{

-			lGammaRplusX[j] = lgamma(this->r + j);

-		}

-		for (int t=0; t<this->T; t++)

-		{

-			lxfactorial = this->lxfactorials[(int) this->O[t]];

-			if (O[t] == 0)

-			{

-				dens[t] = this->w + (1-this->w) * exp( lGammaRplusX[(int) this->O[t]] - lGammaR - lxfactorial + this->r * logp + this->O[t] * log1minusp );

-			}

-			else

-			{

-				dens[t] = (1-this->w) * exp( lGammaRplusX[(int) this->O[t]] - lGammaR - lxfactorial + this->r * logp + this->O[t] * log1minusp );

-			}

-			if (isnan(dens[t]))

-			{

-				FILE_LOG(logWARNING) << __PRETTY_FUNCTION__;

-				FILE_LOG(logWARNING) << "dens["<<t<<"] = "<< dens[t];

-				exit(1);

-			}

-		}

-	}

-	else

-	{

-		FILE_LOG(logDEBUG2) << "Computing gammas in " << __func__ << " for every t, because max(O)>T";

-		for (int t=0; t<this->T; t++)

-		{

-			lGammaRplusX = lgamma(this->r + this->O[t]);

-			lxfactorial = this->lxfactorials[(int) this->O[t]];

-			if (O[t] == 0)

-			{

-				dens[t] = this->w + (1-this->w) * exp( lGammaRplusX - lGammaR - lxfactorial + this->r * logp + this->O[t] * log1minusp );

-			}

-			else

-			{

-				dens[t] = (1-this->w) * exp( lGammaRplusX - lGammaR - lxfactorial + this->r * logp + this->O[t] * log1minusp );

-			}

-			if (isnan(dens[t]))

-			{

-				FILE_LOG(logWARNING) << __PRETTY_FUNCTION__;

-				FILE_LOG(logWARNING) << "dens["<<t<<"] = "<< dens[t];

-				exit(1);

-			}

-		}

-	}

-}

-

-void ZiNB::calc_CDFs(double* CDF)

-{

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

-	double log1minusp = log(1-this->p);

-	double lGammaR = lgamma(this->r);

-	double lppowerr = this->r * log(this->p);

-	// No selection strategy here, because we must precompute CDF to deal with 1s by shifting them

-// 	// Select strategy for computing gammas

-// 	if (this->max_O <= this->T)

-//	{

-		FILE_LOG(logDEBUG2) << "Precomputing gammas in " << __func__ << " for every O[t], because max(O)<=T";

-		double lGamma1plusRplusX[this->max_O+1], lGamma2plusX[this->max_O+1], lHyper[this->max_O+1], lppowert[this->max_O+1];

-		double precomputed_CDF[this->max_O+1];

-		for (int j=0; j<=this->max_O; j++)

-		{

-			lGamma1plusRplusX[j] = lgamma(1 + this->r + j);

-			lGamma2plusX[j] = lgamma(2 + j);

-			lHyper[j] = log(gsl_sf_hyperg_2F1(1, 1 + this->r + j, 2 + j, 1-this->p));

-			lppowert[j] = (1+j) * log1minusp;

-			precomputed_CDF[j] = 1 - exp( log(1-this->w) + lppowerr + lppowert[j] + lGamma1plusRplusX[j] + lHyper[j] - lGammaR - lGamma2plusX[j] );

-			if (precomputed_CDF[j] == 1)

-			{

-				FILE_LOG(logDEBUG4) << "CDF = 1 for O[t] = "<<j<< ", shifting to value of O[t] = "<<j-1;

-				precomputed_CDF[j] = precomputed_CDF[j-1]; 

-			}

-		}

-		for (int t=0; t<this->T; t++)

-		{

-			CDF[t] = precomputed_CDF[(int)O[t]];

-			if (isnan(CDF[t]))

-			{

-				FILE_LOG(logWARNING) << __PRETTY_FUNCTION__;

-				FILE_LOG(logWARNING) << "CDF["<<t<<"] = "<< CDF[t];

-				exit(1);

-			}

-		}

-// 	}

-// 	else

-// 	{

-// 		FILE_LOG(logDEBUG2) << "Computing gammas in " << __func__ << " for every t, because max(O)>T";

-// 		double lGamma1plusRplusX, lGamma2plusX, lHyper, lppowert;

-// 		for (int t=0; t<this->T; t++)

-//		{

-// 			lGamma1plusRplusX = lgamma(1 + this->r + this->O[t]);

-// 			lGamma2plusX = lgamma(2 + this->O[t]);

-// 			lHyper = log(gsl_sf_hyperg_2F1(1, 1 + this->r + this->O[t], 2 + this->O[t], 1-this->p));

-// 			lppowert = (1+this->O[t]) * log1minusp;

-// 			CDF[t] = 1 - exp( log(1-this->w) + lppowerr + lppowert + lGamma1plusRplusX + lHyper - lGammaR - lGamma2plusX ); //TODO: Check formula for log

-// 			if(CDF[t] == 0)

-// 			{

-// 				FILE_LOG(logWARNING) << "CDF["<<t<<"] = "<< CDF[t]; //TODO: Check if this works

-// // 				cout<<"CAUTION!!!! current ="<<current<<endl; current = this->cdf->at(i-1);

-// 			}

-// 			if (isnan(CDF[t]))

-// 			{

-// 				FILE_LOG(logWARNING) << __PRETTY_FUNCTION__;

-// 				FILE_LOG(logWARNING) << "CDF["<<t<<"] = "<< CDF[t];

-// 				exit(1);

-// 			}

-// 		}

-// 	}

-}

-

-void ZiNB::calc_logCDFs(double* logCDF)

-{

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

-	double log1minusp = log(1-this->p);

-	double lGamma1plusRplusX, lHyper, lGammaR, lGamma2plusX, lppowert, lppowerr;

-	lGammaR = lgamma(this->r);

-	lppowerr = this->r * log(this->p);

-	// Select strategy for computing gammas

-	if (this->max_O <= this->T)

-	{

-		FILE_LOG(logDEBUG2) << "Precomputing gammas in " << __func__ << " for every O[t], because max(O)<=T";

-		double lGamma1plusRplusX[this->max_O+1], lGamma2plusX[this->max_O+1], lHyper[this->max_O+1], lppowert[this->max_O+1];

-		for (int j=0; j<=this->max_O; j++)

-		{

-			lGamma1plusRplusX[j] = lgamma(1 + this->r + j);

-			lGamma2plusX[j] = lgamma(2 + j);

-			lHyper[j] = log(gsl_sf_hyperg_2F1(1, 1 + this->r + j, 2 + j, 1-this->p));

-			lppowert[j] = (1+j) * log1minusp;

-		}

-		for (int t=0; t<this->T; t++)

-		{

-			logCDF[t] = log(1 - exp( log(1-this->w) + lppowerr + lppowert[(int)O[t]] + lGamma1plusRplusX[(int)O[t]] + lHyper[(int)O[t]] - lGammaR - lGamma2plusX[(int)O[t]] ));

-			if(logCDF[t] == 0)

-			{

-				FILE_LOG(logWARNING) << "logCDF["<<t<<"] = 0";

-// 				cout<<"CAUTION!!!! current ="<<current<<endl; current = this->cdf->at(i-1);

-			}

-			if (isnan(logCDF[t]))

-			{

-				FILE_LOG(logWARNING) << __PRETTY_FUNCTION__;

-				FILE_LOG(logWARNING) << "logCDF["<<t<<"] = "<< logCDF[t];

-				exit(1);

-			}

-		}

-	}

-	else

-	{

-		FILE_LOG(logDEBUG2) << "Computing gammas in " << __func__ << " for every t, because max(O)>T";

-		for (int t=0; t<this->T; t++)

-		{

-			lGamma1plusRplusX = lgamma(1 + this->r + this->O[t]);

-			lGamma2plusX = lgamma(2 + this->O[t]);

-			lHyper = log(gsl_sf_hyperg_2F1(1, 1 + this->r + this->O[t], 2 + this->O[t], 1-this->p));

-			lppowert = (1+this->O[t]) * log1minusp;

-			logCDF[t] = log(1 - exp( log(1-this->w) + lppowerr + lppowert + lGamma1plusRplusX + lHyper - lGammaR - lGamma2plusX )); //TODO: Check formula for log

-			if(logCDF[t] == 0)

-			{

-				FILE_LOG(logWARNING) << "logCDF["<<t<<"] = 0"; //TODO: Check if this works

-// 				cout<<"CAUTION!!!! current ="<<current<<endl; current = this->cdf->at(i-1);

-			}

-			if (isnan(logCDF[t]))

-			{

-				FILE_LOG(logWARNING) << __PRETTY_FUNCTION__;

-				FILE_LOG(logWARNING) << "logCDF["<<t<<"] = "<< logCDF[t];

-				exit(1);

-			}

-		}

-	}

-}

-

-double ZiNB::getMean()

-{

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

-	return (1-this->w)*this->r*(1-this->p)/this->p;

-}

-

-double ZiNB::getVariance()

-{

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

-	return (1-this->w)*this->r*(1-this->p)/this->p/this->p; //TODO: Is this correct?

-}

-

-DensityName ZiNB::getType()

-{

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

-	return(ZINB);

-}

-

-void ZiNB::setR (double newR)

-{

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

-  this->r = newR;

-}

-

-double ZiNB::getR()

-{

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

-  return(this->r);

-}

-

-void ZiNB::setP (double newP)

-{

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

-  this->p = newP;

-}

-

-double ZiNB::getP()

-{

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

-  return(this->p);

-}

-

-void ZiNB::setW (double newW)

-{

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

-    this->w = newW;

-}

-

-double ZiNB::getW()

-{

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

-	return(this->w);

-}

-

-void ZiNB::copy(Density* other)

-{

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

-	ZiNB* o = (ZiNB*)other;

-	this->p = o->p;

-	this->r = o->r;

-	this->O = o->O;

-	this->w = o->w;

-}

-

-double ZiNB::getLogDensityAt10Variance()

-{

-	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

-	double logp = log(this->p);

-	double log1minusp = log(1-this->p);

-	double lGammaR,lGammaRplusX,lxfactorial;

-	double logdens;

-	// Calculate variance

-	double mean = 0, variance = 0;

-	for(int t=0; t<this->T; t++)

-	{

-		mean += O[t];

-	}

-	mean = mean / this->T;

-	for(int t=0; t<this->T; t++)

-	{

-		variance += pow(O[t] - mean, 2);

-	}

-	variance = variance / this->T;

-	// Calculate logdensity

-	int x = 10*((int)variance+1);

-	lGammaR=lgamma(this->r);

-	lGammaRplusX = lgamma(this->r + x);

-	lxfactorial = this->lxfactorials[x];

-	if (x == 0)

-	{

-		logdens = log( this->w + (1-this->w) * exp( lGammaRplusX - lGammaR - lxfactorial + this->r * logp + x * log1minusp ) );

-	}

-	else

-	{

-		logdens = log(1-this->w) + lGammaRplusX - lGammaR - lxfactorial + this->r * logp + x * log1minusp;

-	}

-	if (isnan(logdens))

-	{

-		FILE_LOG(logWARNING) << __PRETTY_FUNCTION__;

-		FILE_LOG(logWARNING) << "logdens = "<< logdens;

-		exit(1);

-	}