/*
Copyright (c) 2007 Takahiro Shinozaki

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
*/


// This is a demo program of CV-EM and Ag-EM algorithms and is written in 
// MaTX language. It trains GMM using the algorithms and evaluates it. The 
// details of the algorithms are found at
// http://www.furui.cs.titech.ac.jp/~shinot
//
// To compile this program, MaTX is required. MaTX is a freeware whose 
// language syntax is somewhat similar to matlab. Assuming that MaTX is
// installed, run the following command to make an executable binary:
// % matc gmmcvagem.mm
// where % is the command line prompt.
//
// If MaTX is not installed to the standard place, you might need to specify
// matc library and include directories. For example,
// % matc gmmcvagem.mm  -LMaTX/lib/ -IMaTX/include/
// assuming that MaTX library and include directories are at MaTX/lib/ 
// and MaTX/include/.
//
// After the compilation, this program is run by
// ./gmmcvagem
// See the usage printed by the command for more details.

Integer maxint;
Integer warnlv;
Matrix llresults;

// Return Gaussian likelihood
Func Real gcprob(m, v, x)
	Real m, v;
	Real x;
{
	Real prob;
	prob = 1/sqrt(2*PI*v)*exp( -(x-m)^2/(2*v) );
	if (0<warnlv && (prob < 0.0 || !isfinite(prob))) {
		fprintf(2, "WARN: GC prob is out of range. mean=%e var=%e obs=%e\n", m, v, x);
	}
	return prob;
}

// Get sufficient statistics
Func List getss(model, data)
	Matrix model;
	Matrix data;
{
	Integer nmix;
	Integer nsamp;
	Matrix mixPr;
	Integer t,m;
	Real mean, var, weight;
	Real obsPr, lobsPr, totalll, dimProb;
	Matrix dval;
	Matrix ss;
	Integer dim, d;
	Matrix dtmp;

	nmix = Cols(model);
	nsamp = Cols(data);
	dim = Rows(data);
	ss = Z(2*dim+1, nmix);
	dtmp = Z(dim, 1);
	totalll=0.0;
	for (t=1; t<=nsamp; t++) {
		// compute the posterior probability distribution of the hidden variable "mix".
		// i.e. get the observation counts.
		mixPr = Z(1,nmix);
		dval = data(1, t, dtmp);
		obsPr = 0.0;
		for (m=1; m<=nmix; m++) {
			weight = model(1, m);
			dimProb = 1.0;
			for (d=1; d<=dim; d++) {
				mean = model(2*d, m);
				var = model(2*d+1, m);
				dimProb = dimProb * gcprob(mean, var, dval(d));
			}
			mixPr(m) = dimProb * weight;
			if (mixPr(m) < 0.0 || !isfinite(mixPr(m))) {
				fprintf(2, "ERROR: Observation and mix %d joint Prob is out of range : %e\n", m, mixPr(m));
				exit(1);
			}
			obsPr = obsPr + mixPr(m);
		}
		if (obsPr <= 0.0 || !isfinite(obsPr)) {
			if (0<warnlv) {
				fprintf(2, "WARN: Observation prob is out of range : %e\n", obsPr);
			}
			continue;
		}
		for (m=1; m<=nmix; m++) {
			mixPr(m) = mixPr(m)/obsPr;
		}
		// accumulate the statistics
		for (m=1; m<=nmix; m++) {
			ss(1,m) = ss(1,m) + mixPr(m);
			for (d=1; d<=dim; d++) {
				ss(2*d,m) = ss(2*d,m) + dval(d)*mixPr(m);
				ss(2*d+1,m) = ss(2*d+1,m) + dval(d)*dval(d)*mixPr(m);
			}
		}
		lobsPr = log(obsPr);
		if (!isfinite(lobsPr)) {
			fprintf(2, "ERROR: Log likelihood is out of range %e. (likelihood = %e)\n", lobsPr, obsPr);
			exit(1);
		}
		totalll = totalll + lobsPr; // total observation likelihood of the data
		if (!isfinite(totalll)) {
			fprintf(2, "ERROR: Total log likelihood is out of range %e\n", totalll);
			exit(1);
		}
	}
	return {ss, totalll};
}

// Estimate a GMM from a bank of SS
Func Matrix getmodel(bss, dim)
	Array bss;
	Integer dim;
{
	Integer nmix, nbank;
	Matrix ss, tmpss;
	Integer i, m, d;
	Real totcount;
	Matrix model;

	nmix = Cols(bss);
	nbank = Rows(bss)/(2*dim+1);

	// accumulate SS over banks
	ss = Z(2*dim+1, nmix);
	for (i=1; i<=nbank; i++) {
		tmpss = bss(i, 1, ss);
		ss = ss + tmpss;
	}

	// estimate model parameters
	totcount = 0.0;
	for (m=1; m<=nmix; m++) {
		if (ss(1,m) <= 0 || !isfinite(ss(1,m))) {
			if (0<warnlv) {
				fprintf(2, "WARN: mix %d observation count is out of range in model estimation.\n", m);
				fprintf(2, "      SS values: count %e, sum %e, sqsum %e\n", ss(1,m), ss(2,m), ss(3,m));
			}
			// add small counts to the mixture to avoid numeric error
			ss(1, m) = 1.0E-20; ss(2, m) = 0.0; ss(3, m) = 0.0;
		}
		totcount = totcount + ss(1,m);
	}
	if (0<warnlv && (totcount <= 0 || !isfinite(totcount))) {
		fprintf(2, "WARN: total observation count is out of range in model estimation: %e\n", totcount);
	}
	model = Z(2*dim+1, nmix);
	for (m=1; m<=nmix; m++) {
		model(1, m) = ss(1, m)/totcount; // mix weight
		for (d=1; d<=dim; d++) {
			model(2*d, m) = ss(2*d, m)/ss(1, m); // GC mean
			model(2*d+1, m) = ss(2*d+1, m)/ss(1, m) - model(2*d, m)*model(2*d, m); // GC var
		}
		// flooring
		if (model(1, m) < 1.0E-5) { // weight
			model(1, m) = 1.0E-5; 
		}
		for (d=1; d<=dim; d++) {
			if (model(2*d+1, m) < 1.0E-5) {  // var
				model(2*d+1, m) = 1.0E-5; 
			}
		}
	}

	return model;
}

// Evaluate and return test set log likelihood
Func Real evalm(model, test)
	Matrix model;
	Matrix test;
{
	Integer nmix, nsamp;
	Integer t, m;
	Real prob, ll, fll, lltmp;
	Integer num;
	Integer dim, d;
	Real dimProb;

	nmix = Cols(model);
	dim = Rows(test);
	nsamp = Cols(test);
	ll=0.0; num=0;
	for (t=1; t<=Cols(test); t++) {
		prob=0.0;
		for (m=1; m<=nmix; m++) {
			dimProb=1.0;
			for (d=1; d<=dim; d++) {
				dimProb = dimProb * gcprob(model(2*d, m), model(2*d+1, m), test(d, t));
			}
			prob = prob + dimProb * model(1, m);
		}
		lltmp = log(prob);
		if (!isfinite(lltmp)) {
			if (0<warnlv) {
				fprintf(2, "WARN: observation prob is out of range : prob %e ll %e \n", prob, ll);
			}
		} else {
			ll = ll + lltmp;
			num++;
		}
	}
	fll = ll/num;
	return fll;
}


// Sample integers from [1 m] for n times approximately without replacement
Func Matrix SampIntsAWOrep(m, n)
  Integer m;
  Integer n;
{
  Matrix s;
  Index idx;
  Matrix sorted;

  if (m < n) {
    fprintf(2, "ERROR SampIntsAWOrep: input is out of range m=%d n=%d.\n", m, n);
    exit(1);
  }
  if (m<1000 || m/50 < n) {
    // strict sampling without replacement
    {sorted, idx} = sort(rand(1, m));
    return Matrix(idx(1:n));
  } else {
    // approximate with sampling with replacement
    return Matrix(ceil(rand(1,n)*m));
  }
}

Func Integer isdiv(a, b)
	Integer a;
	Integer b;
{
	return ((a/b)*b==a);
}

// Get the number of choosing K out of N.
Func Integer combination(N, K)
	Integer N;
	Integer K;
{
  Integer i, j;
  Integer c;
  Matrix n, d;
  Integer orgK;


  orgK=K;
  if (N-K < K) {
    K = N-K;
  }

  for (i=1; i<=K; i++) {
    n(i) = N-i+1;
  }
  for (i=1; i<=K; i++) {
    d(i) = K-i+1;
  }

  // reduce
  for (i=1; i<=K; i++) {
    for (j=1; j<=K; j++) {
      for (c=Integer(d(j)); 1<c; c--) {
	if (isdiv(Integer(n(i)), c) && isdiv(Integer(d(j)), c)) {
          n(i) = n(i)/c;
          d(j) = d(j)/c;
        }
      }
    }
  }
  c = 1;
  for (i=1; i<=K; i++) {
    if (maxint/(Integer(n(i))+1) < c) {
      fprintf(2, "Error: combination C(N,K) overflow N=%d K=%d\n", N, orgK);
      exit(1);
    }
    c = c * Integer(n(i));
  }
  return c;
}

/*
This version of the combination function is the same as above but easily overflows
Func Integer combination(N, K)
	Integer N;
	Integer K;
{
  Integer i;
  Integer c;

  if (N-K < K) {
    K = N-K;
  }
  c = 1;
  for (i=N; N-K<i; i--) {
    c = c*i;
  }
  for (i=1; i<=K; i++) {
    c = c/i;
  }
  return c;
}
*/

// This function lists all the combination of choosing K out of N in order.
// This function is just for showing how the combinations are sorted and not 
// actually used in the main routine.
Func void listcombinations(N, K)
  Integer N;
  Integer K;
{
  Integer i, j;
  Matrix cmb;

  // Initialize
  cmb = Z(1, K);
  for (i=1; i<=K; i++) {
    cmb(i) = i;
  }
  cmb(K) = cmb(K)-1;
  // list all the combinations
  j = 1;
  while (1) { 
    i=K;
    while (1<=i) {
      if (cmb(i)<N-(K-i)) {
	break;
      }
      i--;
    }
    if (i<1) {
      break;
    } else {
      cmb(i) = cmb(i)+1;
      for (i++; i<=K; i++) {
	cmb(i) = cmb(i-1)+1;
      }
    }
    printf("cmb%d:", j++);
    for (i=1; i<=K; i++) {
      printf(" %d", Integer(cmb(i)));
    }
    printf("\n");
  }
}

// Get idx-th element of the combination list defined by listcombinations(N, K)
Func Matrix
idx2cmb(N, K, idx, cmb)
	Integer N;
	Integer K;
	Integer idx; // 1 <= idx <= C(N,K)
	Matrix cmb;
{
  Integer i, j;
  Integer l, c, sum;

  if (idx < 1 || combination(N, K) < idx) {
    fprintf(2, "idx is out of range idx %d [%d %d]\n", idx, 1, combination(N, K));
    exit(1);
  }
  l = 1;
  for (i=1; i<=K; i++) {
    sum = 0;
    for (j=l; j<=N-K+i; j++) {
      c = combination(N-l, K-i);
      l = j+1;
      if (idx <= sum + c) {
	break; 
      }
      sum = sum + c;
    }
    cmb(i) = j;
    idx = idx - sum;
  }
  return cmb;
}

// Sample C(N,K) for m times (approximately) without replacement
Func Matrix
mCnk(N, K, m)	// return m by K matrix
	Integer N;
	Integer K;
	Integer m;
{
	Integer numc;
	Matrix samples;
	Integer i;
	Matrix retv;
	Integer idx;
	Matrix cmb;

	numc = combination(N, K);
	samples = SampIntsAWOrep(numc, m);
	for (i=1; i<=m; i++) {
		idx = Integer(samples(i));
		cmb = Z(1, K);
		idx2cmb(N, K, idx, cmb);
		retv(i, :) = cmb;
	}
	return retv;
}

// Train and return a GMM
Func Matrix trainGMM(initm, train, test, niter, nlot, nsspm, nmodel, trtype)
	Matrix initm;
	Matrix train, test;
	Integer niter;
	Integer nlot;
	Integer nsspm;   // number of subsets used to estimate a model
	Integer nmodel;  // number of parallel models
	Integer trtype;  // 0 for conventional EM, 1 for CV-EM, 2 for Ag-EM
{
	Matrix lottmp;
	Matrix sstmp;
	Matrix model, bmodel;
	Integer lotsize;
	Matrix bss, cvbss, agbss;
	Integer nmix;
	Integer i, j, k, l;
	Real trll, tsll;
	Integer dim;
	Matrix msstable;
	Index idx;
	Matrix sorted;
	Matrix ssret;
	Real llret;

	nmix = Cols(initm);
	dim = Rows(train);
	lotsize = Cols(train)/nlot; //  fraction is discarded.
	lottmp = Z(dim, lotsize);
	sstmp = Z(2*dim+1, nmix);

	// Evaluate model performance for the 0-th iteration
	model = initm;
	trll = evalm(model, train);
	tsll = evalm(model, test);
	/* printf("Iter %d train ll %e test ll %e\n", 0, trll, tsll); */
	llresults(2, 1) = llresults(2, 1) + trll; // training set likelihood by general model
	llresults(3, 1) = llresults(3, 1) + tsll; // test set likelihood by general model

	bmodel = Z(nmodel, 1, initm); // bank of models
	for (i=1; i<=nmodel; i++) {
		bmodel(i, 1, initm) = initm; // initialize the models with the initial model parameters
	}

	// EM, CV-EM, or Ag-EM iteration
	for (i=1; i<=niter; i++) {
		/* printf("Iter %d\n", i); */
		// E-step
		bss = Z(nlot, 1, sstmp);    // bank of SSs
		if (trtype==0) {
			// p-EM
			for (j=1; j<=nlot; j++) {
				// Although getss is run sequentially for each training subset in this sample code, 
				// it can be run in parallel and that is the purpose of parallel EM to shorten
				// the turn-around time.
				{ssret, llret} = getss(model, train(1, j, lottmp));
				bss(j, 1, sstmp) = ssret; // sufficient statistics
				llresults(1, i) = llresults(1, i) + llret; // E-step training set likelihood
			}
		} else if (trtype==1) {
			// CV-EM
			for (j=1; j<=nlot; j++) {
				// process j-th data subset using j-th CV model
				{ssret, llret} = getss(bmodel(j, 1, initm), train(1, j, lottmp));
				bss(j, 1, sstmp) = ssret; // sufficient statistics
				llresults(1, i) = llresults(1, i) + llret; // E-step training set likelihood
			}
		} else {
			// Ag-EM
			for (j=1; j<=nlot; j++) {
				// process j-th data subset using multiple (nmodel) models
				for (k=1; k<=nmodel; k++) {
					{ssret, llret} = getss(bmodel(k, 1, initm), train(1, j, lottmp));
					bss(j, 1, sstmp) = bss(j, 1, sstmp) + ssret;
					llresults(1, i) = llresults(1, i) + llret; // E-step training set likelihood
				}
				bss(j, 1, sstmp) = bss(j, 1, sstmp) / nmodel; 
			}
		}
		// M-step
		if (trtype==0) {
			// p-EM
			model = getmodel(bss, dim);
			for (j=1; j<=nmodel; j++) {
				bmodel(j, 1, initm) = model; // just copy
			}
		} else if (trtype==1) {
			// CV-EM (nmodel-fold CV)
			for (j=1; j<=nmodel; j++) {
				// exclude j-th subset
				cvbss=Z(nsspm, 1, sstmp);
				l=1;
				for (k=1; k<=nlot; k++) {
					if (k != j) {
						cvbss(l, 1, sstmp) = bss(k, 1, sstmp);
						l = l+1;
					}
				}
				// estimate j-th model
				bmodel(j, 1, initm) = getmodel(cvbss, dim);
			}
			// M-step for general model
			model = getmodel(bss, dim);
		} else if (trtype==2) {
			// Ag-EM
			msstable = mCnk(nlot, nsspm, nmodel); // model set sampling table
			for (j=1; j<=nmodel; j++) {
				agbss=Z(nsspm, 1, sstmp);
				for (k=1; k<=nsspm; k++) {
					agbss(k, 1, sstmp) = bss(Integer(msstable(j, k)), 1, sstmp);
				}
				// estimate j-th model
				bmodel(j, 1, initm) = getmodel(agbss, dim);
			}
			// M-step for general model
			model = getmodel(bss, dim);
		}

		// Evaluate model performance using the integrated model
		trll = evalm(model, train);
		tsll = evalm(model, test);
		/* printf("Iter %d train ll %e test ll %e\n", i, trll, tsll); */
		llresults(2, i+1) = llresults(2, i+1) + trll; // training set likelihood by general model
		llresults(3, i+1) = llresults(3, i+1) + tsll; // test set likelihood by general model
	}
	return model;
}

Func void usage()
{
	String myname;
	myname = "gmmcvagem";
	printf("NAME\n");
	printf("  gmmcvagem - a sample program of CV-EM and Ag-EM algorithms\n");
	printf("\n");
	printf("SYNOPSIS\n");
	printf("  %s em niter nlot lsize\n", myname);
	printf("  %s cv niter nlot lsize\n", myname);
	printf("  %s ag niter nlot lsize nsel enss\n", myname);
	printf("\n");
	printf("DESCRIPTION\n");
	printf("  Train and evaluate GMM using parallel-EM (p-EM), Cross-validation EM(CV-EM),\n");
	printf("  and Aggregated-EM(Ag-EM) algorithms.\n");
	printf("\n");
	printf("  Arguments\n");
	printf("        (em|cv|ag) : \"em\" for p-EM, \"cv\" for CV-EM, \"ag\" for Ag-EM.\n");
	printf("	niter : Number of EM/CV-EM/Ag-EM iterations.\n");
	printf("	nlot  : Number of training subsets.\n");
	printf("	lsize : Number of training samples per subset.\n");
	printf("		Total training set size = nlot x lsize.\n");
	printf("	nsel  : Number of subset selection to estimate a model.\n");
	printf("		(nsel should be around 0.6*nlot for good results)\n");
	printf("	enss  : Ensemble size.\n");
	printf("\n");
	printf("EXAMPLES\n");
	printf("  %s em 10 20 1\n", myname);
	printf("  %s cv 10 20 1\n", myname);
	printf("  %s ag 10 20 1 12 8\n", myname);
	printf("\n");
	printf("AUTHOR\n");
	printf("  Written by Takahiro Shinozaki\n");
	printf("  URL http://www.furui.cs.titech.ac.jp/~shinot\n");
	printf("\n");
	exit(1);
}

Func void main(argc, argv)
	Integer argc;
	List argv;
{
	Integer nmix, niter, nlot, nsspm, nmodel, lsize, tssize, trtype, dim;
	Matrix model, mother;
	Integer m;
	Real norm;
	Real rnd;
	Matrix train, test;
	Integer i, d;
	Integer nexp, t;

	// Initialize global variables
	maxint = 2147483647;
	warnlv = 0;

	// Initialize variables
	nmix = 8;	// Number of Gaussin mixtures
	dim = 4;	// feature dimension
	tssize = 1000;	// test set size
	nexp = 20;	// number of experiments

	printf("gmmcvagem ver.1.0\n");
	printf("Copyright (c) 2007 Takahiro Shinozaki\n");
	printf("\n");
	if (argc > 1) {
		if (argv(2, String) == "em") {
			// EM
			if (argc != 5) {usage();}
			trtype = 0;
			printf("Parallel-EM training\n");
		} else if (argv(2, String) == "cv") {
			// CV-EM
			if (argc != 5) {usage();}
			trtype = 1;
			printf("CV-EM training\n");
		} else if (argv(2, String) == "ag") {
			// Ag-EM
			if (argc != 7) {usage();}
			trtype = 2;
			printf("Ag-EM training\n");
		} else {
			usage();
		}
	} else {
		usage();
	}

	niter = Integer(argv(3, String));
	nlot = Integer(argv(4, String));
	lsize = Integer(argv(5, String));
	if (trtype==0) {
		nmodel = nlot;
		nsspm = nlot;
	} else if (trtype==1) {
		nmodel = nlot;
		nsspm = nlot-1;
	} else if (trtype==2) {
		nsspm = Integer(argv(6, String));
		nmodel = Integer(argv(7, String));
		if (nlot < nsspm) {
			fprintf(2, "nlot(%d) must be smaller than nsspm(%d)\n", nlot, nsspm);
			exit(1);
		}
	}

	printf("train set size = %d (subset size = %d, num of subsets = %d)\n", nlot*lsize, lsize, nlot);
	printf("test set size = %d\n", tssize);
	if (trtype == 2) {
		printf("num subset selection = %d, ensemble size = %d\n", nsspm, nmodel);
	}
	printf("num experiments = %d (Averaged results are output)\n", nexp);

	llresults = Z(3, niter+1);  // 1st row: E-step LL, 2nd: train set LL by global model, 3rd: test set LL by global model
	for (t=1; t<=nexp; t++) { // Repeat experiments with different data set
		// Initialize random seed
		randn("seed", t);
		rand("seed", t);
		// Initialize mother distribution
		mother=Z(2*dim+2, nmix);
		norm = 0.0;
		for (m=1; m<=nmix; m++) {
			mother(1, m) = rand()+0.01;  // weight
			norm = norm + mother(1, m);
			for (d=1; d<=dim; d++) {
				mother(2*d, m) = (rand()-0.5)*50.0;  // mean
				mother(2*d+1, m) = rand()*1.0+0.01;  // var
			}
		}
		for (m=1; m<=nmix; m++) {
			mother(1, m) = mother(1, m)/norm;
		}
		mother(2*dim+2, 1) = mother(1, 1);
		for (m=2; m<=nmix; m++) {
			mother(2*dim+2, m) = mother(2*dim+2, m-1) + mother(1, m);
		}

		// Sample train & test data
		train = Z(dim, nlot*lsize);
		for (i=1; i<=nlot*lsize; i++) {
			// choose mix index according to its probability
			rnd = rand();
			m = 1;
			while (mother(2*dim+2, m)<rnd) {
				m++;
			}
			if (nmix<m) {
				fprintf(2, "ERROR: Assertion error nmix<m %d %d\n", nmix, m);
				exit(1);
			}
			for (d=1; d<=dim; d++) {
				train(d, i) = randn() * sqrt(mother(2*d+1, m)) + mother(2*d, m); // x = z * sigma + mean
			}
		}
		test = Z(dim, tssize);
		for (i=1; i<=tssize; i++) {
			// choose mix index according to its probability
			rnd = rand();
			m = 1;
			while (mother(2*dim+2, m)<rnd) {
				m++;
			}
			if (nmix<m) {
				fprintf(2, "ERROR: Assertion error nmix<m %d %d\n", nmix, m);
				exit(1);
			}
			for (d=1; d<=dim; d++) {
				test(d, i) = randn() * sqrt(mother(2*d+1, m)) + mother(2*d, m); // x = z * sigma + mean
			}
		}

		// Initialize model parameter with global mean and var
		model=Z(2*dim+1, nmix);
		norm = 0.0;
		for (m=1; m<=nmix; m++) {
			model(1, m) = 1.0;  // weight
			for (d=1; d<=dim; d++) {
				model(2*d, m) = mean(train)+(rand()-0.5)*2*std(train);  // mean
				model(2*d+1, m) = std(train)*std(train);  // var
			}
			norm = norm + model(1, m);
		}
		for (m=1; m<=nmix; m++) {
			model(1, m) = model(1, m)/norm;
		}

		// train and evaluate GMM
		trainGMM(model, train, test, niter, nlot, nsspm, nmodel, trtype);
	}
	if (trtype==2) {
		// If Ag-EM, normalize the E-step LL by the number of models as the same
		// training subset is multiply processed.
		for (i=1; i<=niter; i++) {
			llresults(1, i) = llresults(1, i) / nmodel;
		}
	}
	for (i=1; i<=niter; i++) {
		// Normalize by the number of training set samples
		llresults(1, i) = llresults(1, i) / (nlot*lsize);
	}

	printf("\n");
	printf("--Training and evaluation results--\n");
	printf("# 'iter'      is the training iterations.\n");
	printf("# 'E-step-LL' is the per frame E-step log likelihood.\n");
	printf("#             (The likelihood by the first E-step is shown at 0-th iteration)\n");
	printf("# 'train-LL'  is the per frame training set log likelihood by a general model.\n");
	printf("#             (The likelihood by an initial model is shown at 0-th iteration)\n");
	printf("# 'test-LL'   is the per frame test set log likelihood by a general model.\n");
	printf("#             (The likelihood by an initial model is shown at 0-th iteration)\n");
	printf("iter    E-step-LL   train-LL   test-LL\n");
	for (i=1; i<=niter; i++) {
		printf("%-5d %10.3f %10.3f %10.3f\n", i-1, llresults(1, i)/nexp, llresults(2, i)/nexp, llresults(3, i)/nexp);
	}
	i = niter+1;
	printf("%-5d    N/A     %10.3f %10.3f\n", i-1, llresults(2, i)/nexp, llresults(3, i)/nexp);
	printf("Done\n");
}