@@ -1,9 +1,7 @@

 findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max.time=-1, max.iter=-1, num.trials=1, eps.try=NULL, num.threads=1, output.if.not.converged=FALSE, filter.reads=TRUE) {

 	## Intercept user input

-	use.states <- 0:6	# set fixed

 	IDcheck <- ID  #trigger error if not defined

-	if (check.nonnegative.integer.vector(use.states)!=0) stop("argument 'use.states' expects a vector of non-negative integers")

 	if (check.positive(eps)!=0) stop("argument 'eps' expects a positive numeric")

 	if (check.integer(max.time)!=0) stop("argument 'max.time' expects an integer")

 	if (check.integer(max.iter)!=0) stop("argument 'max.iter' expects an integer")

@@ -14,21 +12,19 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max.time=-1, m

 	if (check.positive.integer(num.threads)!=0) stop("argument 'num.threads' expects a positive integer")

 	if (check.logical(output.if.not.converged)!=0) stop("argument 'output.if.not.converged' expects a logical (TRUE or FALSE)")

-

 	war <- NULL

 	if (is.null(eps.try)) eps.try <- eps

 	## Assign variables

-# 	state.labels # assigned globally outside this function

-	use.state.labels <- state.labels[use.states+1]

-	numstates <- length(use.states)

+# 	state.labels # assigned in global.R

+	numstates <- length(state.labels)

 	numbins <- length(binned.data)

 	reads <- mcols(binned.data)$reads

 	iniproc <- which(init==c("standard","random","empiric")) # transform to int

 	# Check if there are reads in the data, otherwise HMM will blow up

 	if (!any(reads!=0)) {

-		stop("All reads in data are zero. No univariate HMM done.")

+		stop("All reads in data are zero. No HMM done.")

 	}

 	# Filter high reads out, makes HMM faster

@@ -37,9 +33,10 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max.time=-1, m

 		mask <- reads > read.cutoff

 		reads[mask] <- read.cutoff

 		numfiltered <- length(which(mask))

-		if (numfiltered > 0) {

-			warning(paste("There are very high read counts in your data (probably artificial). Replaced read counts > ",read.cutoff," (99.99% quantile) by ",read.cutoff," in ",numfiltered," bins. Set option 'filter.reads=FALSE' to disable this filtering.", sep=""))

-		}

+		cat(paste0("Replaced read counts > ",read.cutoff," (99.99% quantile) by ",read.cutoff," in ",numfiltered," bins. Set option 'filter.reads=FALSE' to disable this filtering.\n"))

+# 		if (numfiltered > 0) {

+# 			warning(paste("There are very high read counts in your data (probably artificial). Replaced read counts > ",read.cutoff," (99.99% quantile) by ",read.cutoff," in ",numfiltered," bins. Set option 'filter.reads=FALSE' to disable this filtering.", sep=""))

+# 		}

 	}

@@ -51,7 +48,6 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max.time=-1, m

 			reads = as.integer(reads), # int* O

 			num.bins = as.integer(numbins), # int* T

 			num.states = as.integer(numstates), # int* N

-			use.states = as.integer(use.states), # int* statelabels

 			size = double(length=numstates), # double* size

 			prob = double(length=numstates), # double* prob

 			num.iterations = as.integer(max.iter), #  int* maxiter

@@ -62,6 +58,7 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max.time=-1, m

 			proba = double(length=numstates), # double* proba

 			loglik = double(length=1), # double* loglik

 			weights = double(length=numstates), # double* weights

+			distr.type = as.integer(state.distributions), # int* distr_type

 			ini.proc = as.integer(iniproc), # int* iniproc

 			size.initial = double(length=numstates), # double* initial_size

 			prob.initial = double(length=numstates), # double* initial_prob

@@ -73,22 +70,7 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max.time=-1, m

 			read.cutoff = as.integer(read.cutoff) # int* read_cutoff

 		)

-		names(hmm$weights) <- use.state.labels

 		hmm$eps <- eps.try

-		hmm$A <- matrix(hmm$A, ncol=hmm$num.states, byrow=TRUE)

-		rownames(hmm$A) <- use.state.labels

-		colnames(hmm$A) <- use.state.labels

-		hmm$distributions <- cbind(size=hmm$size, prob=hmm$prob, mu=fmean(hmm$size,hmm$prob), variance=fvariance(hmm$size,hmm$prob))

-		rownames(hmm$distributions) <- use.state.labels

-		# Treat 'null-mixed' separately

-			hmm$distributions[2,'mu'] <- (1-hmm$distributions[2,'prob'])/hmm$distributions[2,'prob']

-			hmm$distributions[2,'variance'] <- hmm$distributions[2,'mu']/hmm$distributions[2,'prob']

-			hmm$distributions[2,'size'] <- NA

-		hmm$A.initial <- matrix(hmm$A.initial, ncol=hmm$num.states, byrow=TRUE)

-		rownames(hmm$A.initial) <- use.state.labels

-		colnames(hmm$A.initial) <- use.state.labels

-		hmm$distributions.initial <- cbind(size=hmm$size.initial, prob=hmm$prob.initial, mu=fmean(hmm$size.initial,hmm$prob.initial), variance=fvariance(hmm$size.initial,hmm$prob.initial))

-		rownames(hmm$distributions.initial) <- use.state.labels

 		if (num.trials > 1) {

 			if (hmm$loglik.delta > hmm$eps) {

 				warning("HMM did not converge in trial run ",i_try,"!\n")

@@ -111,7 +93,6 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max.time=-1, m

 			reads = as.integer(reads), # int* O

 			num.bins = as.integer(numbins), # int* T

 			num.states = as.integer(numstates), # int* N

-			use.states = as.integer(use.states), # int* statelabels

 			size = double(length=numstates), # double* size

 			prob = double(length=numstates), # double* prob

 			num.iterations = as.integer(max.iter), #  int* maxiter

@@ -122,9 +103,10 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max.time=-1, m

 			proba = double(length=numstates), # double* proba

 			loglik = double(length=1), # double* loglik

 			weights = double(length=numstates), # double* weights

+			distr.type = as.integer(state.distributions), # int* distr_type

 			ini.proc = as.integer(iniproc), # int* iniproc

-			size.initial = as.vector(hmm$distributions[,'size']), # double* initial_size

-			prob.initial = as.vector(hmm$distributions[,'prob']), # double* initial_prob

+			size.initial = as.vector(hmm$size), # double* initial_size

+			prob.initial = as.vector(hmm$prob), # double* initial_prob

 			A.initial = as.vector(hmm$A), # double* initial_A

 			proba.initial = as.vector(hmm$proba), # double* initial_proba

 			use.initial.params = as.logical(1), # bool* use_initial_params

@@ -136,27 +118,30 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max.time=-1, m

 	# Add useful entries

 	hmm$ID <- ID

-	names(hmm$weights) <- use.state.labels

+	names(hmm$weights) <- state.labels

 	hmm$coordinates <- data.frame(as.character(seqnames(binned.data)), start(ranges(binned.data)), end(ranges(binned.data)))

 	names(hmm$coordinates) <- coordinate.names

 	hmm$seqlengths <- seqlengths(binned.data)

 	class(hmm) <- class.aneufinder.hmm

-	hmm$states <- factor(use.state.labels, levels=use.state.labels)[hmm$states+1]

+	hmm$states <- factor(state.labels, levels=state.labels)[hmm$states+1]

 	hmm$eps <- eps

 	hmm$A <- matrix(hmm$A, ncol=hmm$num.states, byrow=TRUE)

-	rownames(hmm$A) <- use.state.labels

-	colnames(hmm$A) <- use.state.labels

-	hmm$distributions <- cbind(size=hmm$size, prob=hmm$prob, mu=fmean(hmm$size,hmm$prob), variance=fvariance(hmm$size,hmm$prob))

-	rownames(hmm$distributions) <- use.state.labels

+	rownames(hmm$A) <- state.labels

+	colnames(hmm$A) <- state.labels

+	hmm$distributions <- data.frame(type=state.distributions, size=hmm$size, prob=hmm$prob, mu=fmean(hmm$size,hmm$prob), variance=fvariance(hmm$size,hmm$prob))

+	rownames(hmm$distributions) <- state.labels

 	# Treat 'null-mixed' separately

-		hmm$distributions[2,'mu'] <- (1-hmm$distributions[2,'prob'])/hmm$distributions[2,'prob']

-		hmm$distributions[2,'variance'] <- hmm$distributions[2,'mu']/hmm$distributions[2,'prob']

-		hmm$distributions[2,'size'] <- NA

+	if ('null-mixed' %in% state.labels) {

+		hmm$distributions['null-mixed','mu'] <- (1-hmm$distributions['null-mixed','prob'])/hmm$distributions['null-mixed','prob']

+		hmm$distributions['null-mixed','variance'] <- hmm$distributions['null-mixed','mu']/hmm$distributions['null-mixed','prob']

+		hmm$distributions['null-mixed','size'] <- NA

+	}

 	hmm$A.initial <- matrix(hmm$A.initial, ncol=hmm$num.states, byrow=TRUE)

-	rownames(hmm$A.initial) <- use.state.labels

-	colnames(hmm$A.initial) <- use.state.labels

-	hmm$distributions.initial <- cbind(size=hmm$size.initial, prob=hmm$prob.initial, mu=fmean(hmm$size.initial,hmm$prob.initial), variance=fvariance(hmm$size.initial,hmm$prob.initial))

-	rownames(hmm$distributions.initial) <- use.state.labels

+	rownames(hmm$A.initial) <- state.labels

+	colnames(hmm$A.initial) <- state.labels

+	hmm$distributions.initial <- data.frame(type=state.distributions, size=hmm$size.initial, prob=hmm$prob.initial, mu=fmean(hmm$size.initial,hmm$prob.initial), variance=fvariance(hmm$size.initial,hmm$prob.initial))

+	rownames(hmm$distributions.initial) <- state.labels

+	hmm$distributions.initial['nullsomy',2:5] <- c(0,1,0,0)

 	hmm$filter.reads <- filter.reads

 	# Delete redundant entries

@@ -166,7 +151,7 @@ findCNVs <- function(binned.data, ID, eps=0.001, init="standard", max.time=-1, m

 	hmm$prob.initial <- NULL

 	hmm$use.initial.params <- NULL

 	hmm$read.cutoff <- NULL

-	hmm$use.states <- NULL

+	hmm$distr.type <- NULL

 	# Issue warnings

 	if (num.trials == 1) {

@@ -2,6 +2,7 @@

 # Some global variables that can be used in all functions

 # =======================================================

 state.labels <- c("nullsomy","null-mixed","monosomy","disomy","trisomy","tetrasomy","multisomy")

+state.distributions <- factor(c('delta','dgeom','dnbinom','dnbinom','dnbinom','dnbinom','dnbinom'), levels=c('delta','dgeom','dnbinom'))

 coordinate.names <- c("chrom","start","end")

 binned.data.names <- c(coordinate.names,"reads")

 class.aneufinder.hmm <- "aneufinder.hmm"

@@ -70,13 +70,17 @@ plot.distribution <- function(model, state=NULL, chrom=NULL, start=NULL, end=NUL

 	x <- 0:rightxlim

 	distributions <- list(x)

-	# zero-inflation

-	distributions[[length(distributions)+1]] <- c(weights[1],rep(0,length(x)-1))

-	# geometric

-	distributions[[length(distributions)+1]] <- weights[2] * dgeom(x, model$distributions[2,'prob'])

-	# negative binomials

-	for (istate in 3:numstates) {

-		distributions[[length(distributions)+1]] <- weights[istate] * dnbinom(x, model$distributions[istate,'size'], model$distributions[istate,'prob'])

+	for (istate in 1:nrow(model$distributions)) {

+		if (model$distributions[istate,'type']=='delta') {

+			# zero-inflation

+			distributions[[length(distributions)+1]] <- c(weights[istate],rep(0,length(x)-1))

+		} else if (model$distributions[istate,'type']=='dgeom') {

+			# geometric

+			distributions[[length(distributions)+1]] <- weights[istate] * dgeom(x, model$distributions[istate,'prob'])

+		} else if (model$distributions[istate,'type']=='dnbinom') {

+			# negative binomials

+			distributions[[length(distributions)+1]] <- weights[istate] * dnbinom(x, model$distributions[istate,'size'], model$distributions[istate,'prob'])

+		}

 	}

 	distributions <- as.data.frame(distributions)

 	names(distributions) <- c("x",state.labels)

@@ -7,14 +7,14 @@

 // 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* statelabels, double* size, double* prob, int* maxiter, int* maxtime, double* eps, int* states, double* A, double* proba, double* loglik, double* weights, int* iniproc, double* initial_size, double* initial_prob, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* read_cutoff)

+void R_univariate_hmm(int* O, int* T, int* N, double* size, double* prob, int* maxiter, int* maxtime, double* eps, int* states, double* A, double* proba, double* loglik, double* weights, int* distr_type, int* iniproc, double* initial_size, double* initial_prob, double* initial_A, double* initial_proba, bool* use_initial_params, int* num_threads, int* error, int* read_cutoff)

 {

 	// Define logging level

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

 // 	Output2FILE::Stream() = pFile;

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

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

+//  	FILELog::ReportingLevel() = FILELog::FromString("DEBUG1");

 	// Parallelization settings

 	omp_set_num_threads(*num_threads);

@@ -76,6 +76,7 @@ void R_univariate_hmm(int* O, int* T, int* N, int* statelabels, double* size, do

 	Rprintf("data mean = %g, data variance = %g\n", mean, variance);		

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

+	// This loop assumes that the negative binomial states come last and are consecutive

 	double imean, ivariance;

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

 	{

@@ -91,29 +92,41 @@ void R_univariate_hmm(int* O, int* T, int* N, int* statelabels, double* size, do

 			if (*iniproc == 1)

 			{

 				// Simple initialization based on data mean, assumed to be the disomic mean

-				imean = mean/2 * statelabels[i_state];

-				ivariance = variance/2 * statelabels[i_state];

+				if (distr_type[i_state] == 1) { }

+				else if (distr_type[i_state] == 2) { }

+				else if (distr_type[i_state] == 3)

+				{

+					for (int ii_state=i_state; ii_state<*N; ii_state++)

+					{

+						imean = mean/2 * (ii_state-i_state+1);

+						ivariance = imean * 5;

+// 						ivariance = variance/2 * (i_state-1);

+						// Calculate r and p from mean and variance

+						initial_size[ii_state] = pow(imean,2)/(ivariance-imean);

+						initial_prob[ii_state] = imean/ivariance;

+					}

+					break;

+				}

 			}

-			// Calculate r and p from mean and variance

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

-			initial_prob[i_state] = imean/ivariance;

-

 		}

+	}

-		if (i_state == 0)

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

+	{

+		if (distr_type[i_state] == 1)

 		{

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

+			FILE_LOG(logDEBUG1) << "Using delta distribution for state " << i_state;

 			ZeroInflation *d = new ZeroInflation(O, *T); // delete is done inside ~ScaleHMM()

 			hmm->densityFunctions.push_back(d);

 		}

-		else if (i_state == 1)

+		else if (distr_type[i_state] == 2)

 		{

 			FILE_LOG(logDEBUG1) << "Using geometric distribution for state " << i_state;

 			Geometric *d = new Geometric(O, *T, 0.9); // delete is done inside ~ScaleHMM()

 			hmm->densityFunctions.push_back(d);

 		}

-		else if (i_state >= 2)

+		else if (distr_type[i_state] == 3)

 		{

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

 			NegativeBinomial *d = new NegativeBinomial(O, *T, initial_size[i_state], initial_prob[i_state]); // delete is done inside ~ScaleHMM()

@@ -166,7 +179,7 @@ void R_univariate_hmm(int* O, int* T, int* N, int* statelabels, double* size, do

 		{

 			posterior_per_t[iN] = hmm->get_posterior(iN, t);

 		}

-		states[t] = statelabels[argMax(posterior_per_t, *N)];

+		states[t] = argMax(posterior_per_t, *N);

 	}

 	FILE_LOG(logDEBUG1) << "Return parameters";

@@ -325,8 +325,10 @@ void NegativeBinomial::update(double* weights)

 void NegativeBinomial::update_constrained(double** weights, int fromState, int toState)

 {

 	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

+	FILE_LOG(logDEBUG1) << "r = "<<this->size << ", p = "<<this->prob;

 	double eps = 1e-4, kmax;

 	double numerator, denominator, rhere, dr, Fr, dFrdr, DigammaR, DigammaRplusDR;

+	double logp = log(this->prob);

 	// Update prob (p)

 	numerator=denominator=0.0;

 // 	clock_t time, dtime;

@@ -339,8 +341,7 @@ void NegativeBinomial::update_constrained(double** weights, int fromState, int t

 			denominator+=weights[i+fromState][t]*(this->size*(i+1)+this->obs[t]);

 		}

 	}

-	this->prob = numerator/denominator; // Update of size (r) is now done with updated prob

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

+	this->prob = numerator/denominator; // Update of size (r) is now done with old prob

 // 	dtime = clock() - time;

 // 	FILE_LOG(logDEBUG1) << "updateP(): "<<dtime<< " clicks";

 	// Update of size (r) with Newton Method

@@ -370,13 +371,13 @@ void NegativeBinomial::update_constrained(double** weights, int fromState, int t

 				{

 					if(this->obs[t]==0)

 					{

-						Fr+=weights[i+fromState][t]*logp;

+						Fr+=weights[i+fromState][t]*(i+1)*logp;

 						//dFrdr+=0;

 					}

 					if(this->obs[t]!=0)

 					{

-						Fr+=weights[i+fromState][t]*(logp-DigammaR+DigammaRplusX[(int)obs[t]]);

-						dFrdr+=weights[i+fromState][t]/((i+1)*dr)*(DigammaR-DigammaRplusDR+DigammaRplusDRplusX[(int)obs[t]]-DigammaRplusX[(int)obs[t]]);

+						Fr+=weights[i+fromState][t]*(i+1)*(logp-DigammaR+DigammaRplusX[(int)obs[t]]);

+						dFrdr+=weights[i+fromState][t]/dr*(i+1)*(DigammaR-DigammaRplusDR+DigammaRplusDRplusX[(int)obs[t]]-DigammaRplusX[(int)obs[t]]);

 					}

 				}

 				if(fabs(Fr)<eps)

@@ -405,13 +406,13 @@ void NegativeBinomial::update_constrained(double** weights, int fromState, int t

 					DigammaRplusDRplusX = digamma((i+1)*(rhere+dr)+this->obs[t]); // boost::math::digamma<>(rhere+dr+this->obs[ti]);

 					if(this->obs[t]==0)

 					{

-						Fr+=weights[i+fromState][t]*logp;

+						Fr+=weights[i+fromState][t]*(i+1)*logp;

 						//dFrdr+=0;

 					}

 					if(this->obs[t]!=0)

 					{

-						Fr+=weights[i+fromState][t]*(logp-DigammaR+DigammaRplusX);

-						dFrdr+=weights[i+fromState][t]/((i+1)*dr)*(DigammaR-DigammaRplusDR+DigammaRplusDRplusX-DigammaRplusX);

+						Fr+=weights[i+fromState][t]*(i+1)*(logp-DigammaR+DigammaRplusX);

+						dFrdr+=weights[i+fromState][t]/dr*(i+1)*(DigammaR-DigammaRplusDR+DigammaRplusDRplusX-DigammaRplusX);

 					}

 				}

 			}

@@ -425,6 +426,8 @@ void NegativeBinomial::update_constrained(double** weights, int fromState, int t

 	}

 	this->size = rhere;

 	FILE_LOG(logDEBUG1) << "r = "<<this->size << ", p = "<<this->prob;

+	this->mean = this->fmean(this->size, this->prob);

+	this->variance = this->fvariance(this->size, this->prob);

 // 	dtime = clock() - time;

 // 	FILE_LOG(logDEBUG1) << "updateR(): "<<dtime<< " clicks";

@@ -699,7 +702,7 @@ void Geometric::calc_densities(double* dens)

 	}

 } 

-void Geometric::update(double* weight)

+void Geometric::update(double* weights)

 {

 	FILE_LOG(logDEBUG2) << __PRETTY_FUNCTION__;

 	double numerator, denominator;

@@ -707,8 +710,8 @@ void Geometric::update(double* weight)

 	numerator=denominator=0.0;

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

 	{

-		numerator+=weight[t];

-		denominator+=weight[t]*(1+this->obs[t]);

+		numerator+=weights[t];

+		denominator+=weights[t]*(1+this->obs[t]);

 	}

 	this->prob = numerator/denominator;

 	FILE_LOG(logDEBUG1) << "p = "<<this->prob;

@@ -256,45 +256,39 @@ void ScaleHMM::baumWelch(int* maxiter, int* maxtime, double* eps)

 		clock_t clocktime = clock(), dtime;

-// 		// Update distribution of state 1 (null-mixed)

-// 		this->densityFunctions[1]->update(this->gamma[1]);

-// 		// Update distribution of state 2 (monosomic)

-// 		int iN=2;

-// 		this->densityFunctions[iN]->update_constrained(this->gamma, iN, this->N);

-// 		double mean1 = this->densityFunctions[iN]->get_mean();

-// 		double variance1 = this->densityFunctions[iN]->get_variance();

-// 		// Set others as multiples

+// 		// Update all distributions independantly

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

 // 		{

-// 			if (iN!=1 and iN!=2)

-// 			{

-// 				this->densityFunctions[iN]->set_mean(mean1 * (iN-1));

-// 				this->densityFunctions[iN]->set_variance(variance1 * (iN-1));

-// 			}

-// 			FILE_LOG(logDEBUG1) << "mean(state="<<iN<<") = " << this->densityFunctions[iN]->get_mean();

+// 			this->densityFunctions[iN]->update(this->gamma[iN]);

 // 		}

-		

-		// Update distribution of state 3 (disomic)

-		this->densityFunctions[3]->update(this->gamma[3]);

-		double mean2 = this->densityFunctions[3]->get_mean();

-		double variance2 = this->densityFunctions[3]->get_variance();

-		// Update distribution of state 1 (null-mixed)

-		this->densityFunctions[1]->update(this->gamma[1]);

-// 		// Set mean of state 1 (null-mixed) to half of monosomic

-// 		this->densityFunctions[1]->set_mean(mean2/2 / 2);

-// 		this->densityFunctions[1]->set_variance(variance2/2 / 2);

-		// Set the others as multiples of disomic

+

+		// Update distribution of state 'null-mixed' and 'monosomic', set others as multiples of 'monosomic'

+		// This loop assumes that the negative binomial states come last and are consecutive

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

 		{

-			if (iN!=1 and iN!=3)

+			if (this->densityFunctions[iN]->get_name() == ZERO_INFLATION) {}

+			if (this->densityFunctions[iN]->get_name() == GEOMETRIC)

 			{

-				this->densityFunctions[iN]->set_mean(mean2/2 * (iN-1));

-				this->densityFunctions[iN]->set_variance(variance2/2 * (iN-1));

+				this->densityFunctions[iN]->update(this->gamma[iN]);

+			}

+			if (this->densityFunctions[iN]->get_name() == NEGATIVE_BINOMIAL)

+			{

+				FILE_LOG(logDEBUG1) << "mean(state="<<iN<<") = " << this->densityFunctions[iN]->get_mean() << ", var(state="<<iN<<") = " << this->densityFunctions[iN]->get_variance();

+				this->densityFunctions[iN]->update_constrained(this->gamma, iN, this->N);

+				double mean1 = this->densityFunctions[iN]->get_mean();

+				double variance1 = this->densityFunctions[iN]->get_variance();

+				FILE_LOG(logDEBUG1) << "mean(state="<<iN<<") = " << this->densityFunctions[iN]->get_mean() << ", var(state="<<iN<<") = " << this->densityFunctions[iN]->get_variance();

+				// Set others as multiples

+				for (int jN=iN+1; jN<this->N; jN++)

+				{

+					this->densityFunctions[jN]->set_mean(mean1 * (jN-iN+1));

+					this->densityFunctions[jN]->set_variance(variance1 * (jN-iN+1));

+				FILE_LOG(logDEBUG1) << "mean(state="<<jN<<") = " << this->densityFunctions[jN]->get_mean() << ", var(state="<<jN<<") = " << this->densityFunctions[jN]->get_variance();

+				}

+				break;

 			}

 		}

-

-

-

+			

 		dtime = clock() - clocktime;

 	 	FILE_LOG(logDEBUG) << "updating distributions: " << dtime << " clicks";

 		R_CheckUserInterrupt();