<|start_filename|>cli.c<|end_filename|>
#include <ctype.h>
#include <unistd.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <assert.h>
#include <zlib.h>
#include "psmc.h"
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)

#define DEFAULT_PATTERN "4+5*3+4"
#define PSMC_VERSION "0.6.5-r67"

static char conv_table[256] = {
	 2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,
	 2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,
	 2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,
	 0, 1, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,
	 2, 0, 2, 0,  2, 2, 2, 0,  2, 2, 2, 1,  2, 1, 2, 2,
	 2, 2, 1, 1,  0, 2, 2, 1,  2, 1, 2, 2,  2, 2, 2, 2,
	 2, 0, 2, 0,  2, 2, 2, 0,  2, 2, 2, 1,  2, 1, 2, 2,
	 2, 2, 1, 1,  0, 2, 2, 1,  2, 1, 2, 2,  2, 2, 2, 2,
	 2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,
	 2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,
	 2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,
	 2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,
	 2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,
	 2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,
	 2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,
	 2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2,  2, 2, 2, 2
};

static void usage(psmc_par_t *p)
{
	fprintf(stderr, "\n");
	fprintf(stderr, "Program: psmc (Pairwise SMC Model)\n");
	fprintf(stderr, "Version: %s\n", PSMC_VERSION);
	fprintf(stderr, "Contact: <http://hengli.uservoice.com/>\n\n");
	fprintf(stderr, "Usage:   psmc [options] input.txt\n\n");
	fprintf(stderr, "Options: -p STR      pattern of parameters [%s]\n", DEFAULT_PATTERN);
	fprintf(stderr, "         -t FLOAT    maximum 2N0 coalescent time [%lg]\n", p->max_t);
	fprintf(stderr, "         -N INT      maximum number of iterations [%d]\n", p->n_iters);
	fprintf(stderr, "         -r FLOAT    initial theta/rho ratio [%lg]\n", p->tr_ratio);
	fprintf(stderr, "         -c FILE     CpG counts generated by cntcpg [null]\n");
//	fprintf(stderr, "         -l FLOAT    time skipping factor [0.1]\n"); // maybe not to change this
	fprintf(stderr, "         -o FILE     output file [stdout]\n");
	fprintf(stderr, "         -i FILE     input parameter file [null]\n");
//	fprintf(stderr, "         -I FLOAT    amplitude of random initialization for EM iteration [%lg]\n", p->ran_init);
	fprintf(stderr, "         -T FLOAT    initial divergence time; -1 to disable [-1]\n");
	fprintf(stderr, "         -b          bootstrap (input be preprocessed with split_psmcfa)\n");
	fprintf(stderr, "         -S          simulate sequence\n");
	fprintf(stderr, "         -d          perform decoding\n");
	fprintf(stderr, "         -D          print full posterior probabilities\n");
	fprintf(stderr, "\n");
	psmc_delete_par(p);
	exit(1);
}

// parse a pattern like "4+5*3+4"
// the returned array holds which parameters are linked together
// number of parameters and number of free parameters will be also returned

int *psmc_parse_pattern(const char *pattern, int *n_free, int *n_pars)
{
	char *q, *p, *tmp;
	int top = 0, *stack = (int*)malloc(sizeof(int) * 0x100);
	int *pars_map, k, l, i;
	p = q = tmp = strdup(pattern);
	k = 1;
	while (1) {
		assert(isdigit(*p) || *p == '*' || *p == '+' || *p == '\0'); // allowed characters
		if (*p == '+' || *p == '\0') {
			int is_end = (*p == 0)? 1 : 0;
			*p++ = '\0';
			l = atoi(q); q = p;
			for (i = 0; i < k; ++i) {
				stack[top++] = l;
				assert(top <= 0xff);
			}
			k = 1;
			if (is_end) break;
		} else if (*p == '*') {
			*p = '\0';
			k = atoi(q); // number of repeats
			*p++ = '*'; q = p;
		} else ++p;
	}
	for (k = l = 0; k != top; ++k) l += stack[k];
	*n_pars = l - 1; *n_free = top;
	pars_map = (int*)malloc(sizeof(int) * (*n_pars + 1));
	for (k = i = 0; k != top; ++k)
		for (l = 0; l < stack[k]; ++l)
			pars_map[i++] = k;
	free(tmp); free(stack);
	return pars_map;
}

// read the input

void psmc_read_seq(const char *fn, psmc_par_t *pp)
{
	kseq_t *ks;
	int i, j;
	gzFile fp;
	fp = strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(fileno(stdin), "r");
	ks = kseq_init(fp);
	while (kseq_read(ks) >= 0) {
		char *s;
		int L_e = 0, n_e = 0, L = 0;
		if ((pp->n_seqs&0xff) == 0)
			pp->seqs = (psmc_seq_t*)realloc(pp->seqs, sizeof(psmc_seq_t) * (pp->n_seqs + 0x100));
		pp->seqs[pp->n_seqs].name = strdup(ks->name.s);
		s = pp->seqs[pp->n_seqs].seq = (char*)calloc(ks->seq.l + 1, 1);
		L = 0;
		for (j = 0; j != ks->seq.l; ++j) {
			char c = conv_table[(int)ks->seq.s[j]];
			s[L++] = c;
			if (c < 2) {
				++L_e;
				if (c == 1) ++n_e;
			}
		}
		pp->seqs[pp->n_seqs].L_e = L_e;
		pp->seqs[pp->n_seqs].n_e = n_e;
		pp->seqs[pp->n_seqs++].L = L;
	}
	kseq_destroy(ks);
	gzclose(fp);
	// count the sum of lengths and number of segsites
	pp->sum_n = pp->sum_L = 0;
	for (i = 0; i != pp->n_seqs; ++i) {
		pp->sum_n += pp->seqs[i].n_e;
		pp->sum_L += pp->seqs[i].L_e;
	}
}

// parse the command-line options

psmc_par_t *psmc_new_par()
{
	psmc_par_t *par;
	par = (psmc_par_t*)calloc(1, sizeof(psmc_par_t));
	par->n_iters = 30;
	par->fpout = stdout;
	par->max_t = 15.0;
	par->tr_ratio = 4.0;
	par->alpha = 0.1;
	par->dt0 = -1.0;
	return par;
}

void psmc_delete_par(psmc_par_t *par)
{
	int i;
	if (par == 0) return;
	free(par->par_map); free(par->inp_pa);
	free(par->pre_fn); free(par->pattern);
	fclose(par->fpout);
	for (i = 0; i != par->n_seqs; ++i) {
		free(par->seqs[i].name);
		free(par->seqs[i].seq);
	}
	if (par->fpcnt) fclose(par->fpcnt);
	free(par->seqs);
	free(par);
}

psmc_par_t *psmc_parse_cli(int argc, char *argv[])
{
	int c, is_bootstrap = 0;
	psmc_par_t *par;
	par = psmc_new_par();
	while ((c = getopt(argc, argv, "i:t:l:r:N:p:o:dI:c:bSDT:")) >= 0) {
		switch (c) {
		case 'S': par->flag |= PSMC_F_SIMU; break;
		case 'd': par->flag |= PSMC_F_DECODE; break;
		case 'D': par->flag |= PSMC_F_FULLDEC; break;
		case 'T': par->dt0 = atof(optarg); if (par->dt0 >= 0) par->flag |= PSMC_F_DIVERG; break;
		case 't': par->max_t = atof(optarg); break;
		case 'l': par->alpha = atof(optarg); break;
		case 'r': par->tr_ratio = atof(optarg); break;
		case 'N': par->n_iters = atoi(optarg); break;
		case 'i': par->pre_fn = strdup(optarg); break;
		case 'o': par->fpout = fopen(optarg, "w"); assert(par->fpout); break;
		case 'c': par->fpcnt = fopen(optarg, "rb"); assert(par->fpcnt); break;
		case 'p': par->pattern = strdup(optarg); break;
		case 'I': par->ran_init = atof(optarg); break;
		case 'b': is_bootstrap = 1; break;
		default:  usage(par);
		}
	}
	if (par->flag & PSMC_F_FULLDEC) par->flag |= PSMC_F_DECODE;
	if (optind == argc) usage(par);
	fprintf(par->fpout, "CC\n");
	fprintf(par->fpout, "CC\tBrief Description of the file format:\n");
	fprintf(par->fpout, "CC\t  CC  comments\n");
	fprintf(par->fpout, "CC\t  MM  useful-messages\n");
	fprintf(par->fpout, "CC\t  RD  round-of-iterations\n");
	fprintf(par->fpout, "CC\t  LL  \\log[P(sequence)]\n");
	fprintf(par->fpout, "CC\t  QD  Q-before-opt Q-after-opt\n");
	fprintf(par->fpout, "CC\t  TR  \\theta_0 \\rho_0\n");
	fprintf(par->fpout, "CC\t  RS  k t_k \\lambda_k \\pi_k \\sum_{l\\not=k}A_{kl} A_{kk}\n");
	fprintf(par->fpout, "CC\t  DC  begin end best-k t_k+\\Delta_k max-prob\n");
	fprintf(par->fpout, "CC\n");
	if (par->pre_fn) { // load parameters from input
		psmc_read_param(par);
	} else {
		if (par->pattern == 0) par->pattern = strdup(DEFAULT_PATTERN);
		par->par_map = psmc_parse_pattern(par->pattern, &par->n_free, &par->n);
	}
	fprintf(par->fpout, "MM\tVersion: %s\n", PSMC_VERSION);
	fprintf(par->fpout, "MM\tpattern:%s, n:%d, n_free_lambdas:%d\n", par->pattern, par->n, par->n_free);
	fprintf(par->fpout, "MM\tn_iterations:%d, skip:1, max_t:%g, theta/rho:%g\n",
			par->n_iters, par->max_t, par->tr_ratio);
	fprintf(par->fpout, "MM\tis_decoding:%d\n", (par->flag&PSMC_F_DECODE)? 1 : 0);
	psmc_read_seq(argv[optind], par);
	if (is_bootstrap) psmc_resamp(par); // resampling if required
	fprintf(par->fpout, "MM\tn_seqs:%d, sum_L:%lld, sum_n:%d\n", par->n_seqs, (long long)par->sum_L, par->sum_n);
	fflush(par->fpout);
	return par;
}


<|start_filename|>psmc.h<|end_filename|>
#ifndef PSMC_H_LH3
#define PSMC_H_LH3

#include <sys/types.h>
#include <stdio.h>
#include "khmm.h"

#ifndef FLOAT
#define FLOAT double
#endif

#define PSMC_N_PARAMS 3

#define PSMC_T_INF 1000.0

#define PSMC_F_DECODE   0x1
#define PSMC_F_FULLDEC  0x2
#define PSMC_F_SIMU     0x4
#define PSMC_F_DIVERG   0x8

typedef struct
{
	int L, L_e, n_e;
	char *seq, *name;
} psmc_seq_t;

typedef struct
{
	int n; // $n$ in psmc.tex. number of intervals equals to $n+1$
	int n_free; // number of free lambdas
	int *par_map; // parameter groups

	int n_seqs; // number of sequences
	psmc_seq_t *seqs; // sequences of segregating sites
	int64_t sum_L;
	int sum_n;

	FILE *fpout, *fpcnt;

	FLOAT *inp_pa; // parameters from the input
	FLOAT *inp_ti;

	int flag;
	int n_iters; // number of iterations
	char *pre_fn; // previous results
	char *pattern; // pattern
	FLOAT max_t; // maximum time
	FLOAT tr_ratio; // theta/rho ratio
	FLOAT alpha;
	FLOAT ran_init;
	FLOAT dt0;
} psmc_par_t;

typedef struct
{
	hmm_par_t *hp;
	FLOAT *t; // time boundaries
	int n_params; // equals to (psmc_par_t::n_free + 2)
	FLOAT *params; // all free parameters: \theta_0, \rho_0 and free \lambda_k
	FLOAT *sigma, *post_sigma;

	FLOAT C_pi, C_sigma; // normalization factor
	FLOAT lk; // log-likelihood
	FLOAT lk_lo; // log-likelihood of the left-out sequence
	FLOAT Q0, Q1; // EM-Q before and after maximization
} psmc_data_t;

#ifdef __cplusplus
extern "C" {
#endif

	// functions for parsing command-line options and inputs
	int *psmc_parse_pattern(const char *pattern, int *n_free, int *n_pars);
	void psmc_read_param(psmc_par_t *pp);
	psmc_par_t *psmc_parse_cli(int argc, char *argv[]);
	void psmc_delete_par(psmc_par_t *pp);
	void psmc_read_segsites(const char *fn, psmc_par_t *pp);
	void psmc_resamp(psmc_par_t *pp);
	// new/delete psmc_data_t
	psmc_data_t *psmc_new_data(psmc_par_t *pp);
	void psmc_delete_data(psmc_data_t *pd);
	// calculate HMM parameters given population parameters
	void psmc_update_hmm(const psmc_par_t *pp, psmc_data_t *pd);
	void psmc_avg_t(const psmc_par_t *pp, const psmc_data_t *pd, double *avg_t);
	// Baum-Welch algorithm
	FLOAT psmc_em(psmc_par_t *pp, psmc_data_t *data);
	// calculate and print decoding
	void psmc_decode(const psmc_par_t *pp, const psmc_data_t *pd);
	// print parameters
	void psmc_print_data(const psmc_par_t *pp, const psmc_data_t *pd);
	// simulate sequences
	void psmc_simulate(const psmc_par_t *pp, const psmc_data_t *pd);

	void psmc_update_intv(int n, FLOAT t[], FLOAT max_t, FLOAT alpha, FLOAT *inp_ti);

#ifdef __cplusplus
}
#endif

#endif


<|start_filename|>core.c<|end_filename|>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "psmc.h"

void psmc_update_intv(int n, FLOAT t[], FLOAT max_t, FLOAT alpha, FLOAT *inp_ti)
{
	int k;
	if (inp_ti == 0) {
		FLOAT beta;
		beta = log(1.0 + max_t / alpha) / n; // beta controls the sizes of intervals
		for (k = 0; k < n; ++k)
			t[k] = alpha * (exp(beta * k) - 1);
		t[n] = max_t; t[n+1] = PSMC_T_INF; // the infinity: exp(PSMC_T_INF) > 1e310 = inf
	} else {
		memcpy(t, inp_ti, (n+1) * sizeof(FLOAT));
		t[n+1] = PSMC_T_INF;
	}
}

psmc_data_t *psmc_new_data(psmc_par_t *pp)
{
	psmc_data_t *pd;
	int k, n = pp->n;
	pd = (psmc_data_t*)calloc(1, sizeof(psmc_data_t));
	pd->n_params = pp->n_free + PSMC_N_PARAMS + ((pp->flag & PSMC_F_DIVERG)? 1 : 0); // one addition parameter for the divergence model
	pd->hp = hmm_new_par(2, n + 1);
	// initialize
	pd->sigma = (FLOAT*)calloc(n+1, sizeof(FLOAT));
	pd->post_sigma = (FLOAT*)calloc(n+1, sizeof(FLOAT));
	pd->t = (FLOAT*)malloc(sizeof(FLOAT) * (n + 2)); // $t_0,\ldots,t_{n+1}$
	pd->params = (FLOAT*)calloc(pd->n_params, sizeof(FLOAT)); // free lambdas + theta + rho
	// initialize t[] and params[]
	if (pp->inp_pa) { // pameters are loaded from a file
		memcpy(pd->params, pp->inp_pa, sizeof(FLOAT) * pd->n_params); // FIXME: not working for the divergence model
	} else {
		FLOAT theta;
		// initialize psmc_data_t::params[]
		theta = -log(1.0 - (FLOAT)pp->sum_n / pp->sum_L);
		pd->params[0] = theta; // \theta_0
		pd->params[1] = theta / pp->tr_ratio; // \rho_0
		pd->params[2] = pp->max_t;
		for (k = PSMC_N_PARAMS; k != pp->n_free + PSMC_N_PARAMS; ++k) {
			// \lambda_k
			pd->params[k] = 1.0 + (drand48() * 2.0 - 1.0) * pp->ran_init;
			if (pd->params[k] < 0.1) pd->params[k] = 0.1;
		}
		if (pp->flag & PSMC_F_DIVERG) pd->params[pd->n_params - 1] = pp->dt0;
	}
	psmc_update_hmm(pp, pd);
	return pd;
}
void psmc_delete_data(psmc_data_t *pd)
{
	if (pd == 0) return;
	free(pd->sigma); free(pd->post_sigma);
	free(pd->t); free(pd->params);
	hmm_delete_par(pd->hp);
	free(pd);
}
void psmc_update_hmm(const psmc_par_t *pp, psmc_data_t *pd) // calculate the a_{kl} and e_k(b)
{
	FLOAT *q, tmp, sum_t, max_t, *alpha, *beta, *q_aux, *lambda, theta, rho, *t, *tau, dt = 0;
	hmm_par_t *hp = pd->hp;
	int k, l, n = pp->n;
	t = pd->t;
	lambda = (FLOAT*)malloc(sizeof(FLOAT) * (n + 1)); // \lambda_k
	alpha = (FLOAT*)malloc(sizeof(FLOAT) * (n + 2)); // \alpha_k
	beta = (FLOAT*)malloc(sizeof(FLOAT) * (n + 1)); // \beta_k
	q_aux = (FLOAT*)malloc(sizeof(FLOAT) * n); // for acceleration
	q = (FLOAT*)malloc(sizeof(FLOAT) * (n + 1)); // q_{kl}
	tau = (FLOAT*)malloc(sizeof(FLOAT) * (n + 1)); // \tau_k
	// calculate population parameters: \theta_0, \rho_0, \lambda_k and max_t
	theta = pd->params[0]; rho = pd->params[1]; max_t = pd->params[2];
	for (k = 0; k <= n; ++k)
		lambda[k] = pd->params[pp->par_map[k] + PSMC_N_PARAMS];
	psmc_update_intv(pp->n, pd->t, max_t, pp->alpha, pp->inp_ti);
	// set the divergence time parameter if necessary
	if (pp->flag & PSMC_F_DIVERG) {
		dt = pd->params[pd->n_params - 1];
		if (dt < 0) dt = 0;
	}
	// calculate \tau_k
	for (k = 0; k <= n; ++k) tau[k] = t[k+1] - t[k];
	// calculate \alpha
	for (k = 1, alpha[0] = 1.0; k <= n; ++k)
		alpha[k] = alpha[k-1] * exp(-tau[k-1] / lambda[k-1]);
	alpha[k] = 0.0;
	// calculate \beta
	for (k = 1, beta[0] = 0.0; k <= n; ++k)
		beta[k] = beta[k-1] + lambda[k-1] * (1.0 / alpha[k] - 1.0 / alpha[k-1]);
	// calculate q_aux
	for (l = 0; l < n; ++l)
		q_aux[l] = (alpha[l] - alpha[l+1]) * (beta[l] - lambda[l] / alpha[l]) + tau[l];
	// calculate C_pi and C_sigma
	for (l = 0, pd->C_pi = 0.0; l <= n; ++l)
		pd->C_pi += lambda[l] * (alpha[l] - alpha[l+1]);
	pd->C_sigma = 1.0 / (pd->C_pi * rho) + 0.5;
	// calculate all the rest
	for (k = 0, sum_t = 0.0; k <= n; ++k) {
		FLOAT *aa, avg_t, ak1, lak, pik, cpik;
		ak1 = alpha[k] - alpha[k+1]; lak = lambda[k]; // just for convenient
		// calculate $\pi_k$, $\sigma_k$ and Lak
		cpik = ak1 * (sum_t + lak) - alpha[k+1] * tau[k];
		pik = cpik / pd->C_pi;
		pd->sigma[k] = (ak1 / (pd->C_pi * rho) + pik / 2.0) / pd->C_sigma;
		// calculate avg_t, the average time point where mutation happens
		avg_t = - log(1.0 - pik / (pd->C_sigma*pd->sigma[k])) / rho;
		if (isnan(avg_t) || avg_t < sum_t || avg_t > sum_t + tau[k]) // in case something bad happens
			avg_t = sum_t + (lak - tau[k] * alpha[k+1] / (alpha[k] - alpha[k+1]));
		// calculate q_{kl}
		tmp = ak1 / cpik;
		for (l = 0; l < k; ++l) q[l] = tmp * q_aux[l]; // q_{kl}, l<k
		q[l++] = (ak1 * ak1 * (beta[k] - lak/alpha[k]) + 2*lak*ak1 - 2*alpha[k+1]*tau[k]) / cpik; // q_{kk}
		if (k < n) {
			tmp = q_aux[k] / cpik;
			for (; l <= n; ++l) q[l] = (alpha[l] - alpha[l+1]) * tmp; // q_{kl}, l>k
		}
		// calculate p_{kl} and e_k(b)
		tmp = pik / (pd->C_sigma * pd->sigma[k]);
		for (aa = hp->a[k], l = 0; l <= n; ++l) aa[l] = tmp * q[l];
		aa[k] = tmp * q[k] + (1.0 - tmp);
		hp->a0[k] = pd->sigma[k];
		hp->e[0][k] = exp(-theta * (avg_t + dt));
		hp->e[1][k] = 1.0 - hp->e[0][k];
		// update sum_lt
		sum_t += tau[k];
		// for (l = 0, tmp = 0.0; l <= n; ++l) tmp += q[l]; fprintf(stderr, "%d\t%lf\n", k, tmp); // for testing only
	}
	// for (l = 0, tmp = 0.0; l <= n; ++l) tmp += hp->a0[l]; fprintf(stderr, "%lf\n", tmp); // for testing only
	// free
	free(q); free(alpha); free(beta); free(q_aux); free(lambda); free(tau);
}
// compute the average time for each interval; a strip-down version of psmc_update_hmm()
void psmc_avg_t(const psmc_par_t *pp, const psmc_data_t *pd, double *avg_t)
{
	FLOAT sum_t, *alpha, *lambda, rho = pd->params[1], *tau, dt = 0;
	int k, n = pp->n;
	lambda = (FLOAT*)alloca(sizeof(FLOAT) * (n + 1)); // \lambda_k
	alpha = (FLOAT*)alloca(sizeof(FLOAT) * (n + 2)); // \alpha_k
	tau = (FLOAT*)alloca(sizeof(FLOAT) * (n + 1)); // \tau_k
	for (k = 0; k <= n; ++k)
		lambda[k] = pd->params[pp->par_map[k] + PSMC_N_PARAMS];
	if (pp->flag & PSMC_F_DIVERG) {
		dt = pd->params[pd->n_params - 1];
		if (dt < 0) dt = 0;
	}
	for (k = 0; k <= n; ++k) tau[k] = pd->t[k+1] - pd->t[k];
	for (k = 1, alpha[0] = 1.0; k <= n; ++k)
		alpha[k] = alpha[k-1] * exp(-tau[k-1] / lambda[k-1]);
	alpha[k] = 0.0;
	for (k = 0, sum_t = 0.0; k <= n; ++k) {
		FLOAT ak1, lak, pik;
		ak1 = alpha[k] - alpha[k+1]; lak = lambda[k]; // just for convenient
		pik = (ak1 * (sum_t + lak) - alpha[k+1] * tau[k]) / pd->C_pi;
		avg_t[k] = - log(1.0 - pik / (pd->C_sigma*pd->sigma[k])) / rho;
		if (isnan(avg_t[k]) || avg_t[k] < sum_t || avg_t[k] > sum_t + tau[k]) // in case something bad happens
			avg_t[k] = sum_t + (lak - tau[k] * alpha[k+1] / (alpha[k] - alpha[k+1]));
		avg_t[k] += dt;
		sum_t += tau[k];
	}
}
