// file: $isip/class/pr/RegressionDecisionTreeNode/rdtnod_06.cc
// version: $Id: rdtnod_06.cc 9920 2004-12-23 21:07:46Z gao $
//

// isip include files
//
#include "RegressionDecisionTreeNode.h"
#include <MixtureModel.h>

// method: nodeScore
//
// arguments:
//  Vector<StatisticalModel>& stat_models: (input) vector of statistical
// 
// return: logical error status
//
// Calculate the node score
//
bool8 RegressionDecisionTreeNode::nodeScore(Vector<StatisticalModel>& stat_models_a) {

  // local variables
  //
  VectorFloat mean;
  float64 num_occ = 0.0;
  cluster_score_d = 0.0;

  // if in debug mode, print the average mean for this node
  //
  if (debug_level_d >= Integral::ALL) {  
    average_mean_d.debug(L"average_mean_d=");
  }
  
  // loop all models in the node
  //
  for (bool8 more_index = gaussian_models_d.gotoFirst(); more_index;
       more_index = gaussian_models_d.gotoNext()) {

    // get the first datapoint in triple
    //
    RDataPoint* datapoint = gaussian_models_d.getCurr();
    
    // get the statistical model index
    //
    int32 sm_index = (int32)datapoint->first();
    
    // get the gaussian model index
    //
    int32 gm_index = (int32)datapoint->second();

    num_occ = stat_models_a(sm_index).getOccupancy();      
        
    // if the model to adapt is a single Gaussian model, compute its
    // contribution to the global adaptation
    //
    if (stat_models_a(sm_index).getType() == StatisticalModel::GAUSSIAN_MODEL) {
      GaussianModel& gm = stat_models_a(sm_index).getGaussianModel();
      
      gm.getMean(mean);
    }
    
    // if the model to adapt is a mixture of Gaussian models, compute
    // the contribution of mixture component match the index to the
    // global adaptation
    //
    else if (stat_models_a(sm_index).getType() == StatisticalModel::MIXTURE_MODEL) {

      // get the linked list of mixture components
      //
      SingleLinkedList<StatisticalModel>& models =
	const_cast<StatisticalModel&>(stat_models_a(sm_index)).
	getMixtureModel().getModels();

      // check whether each mixture component is Gaussian model and
      // compute its contribution to the global adaptation
      //
      bool8 more = models.gotoFirst();
      int32 local_gaussian_index = 0;
      
      while (more) {

	// get the current mixture component
	//
	StatisticalModel* model = models.getCurr();

	// check whether it is Gaussian
	//
	if (model->getType() == StatisticalModel::GAUSSIAN_MODEL) {
	  GaussianModel& gm = model->getGaussianModel();
	  
	  // compute the contribution of this Gaussian model to the
	  // global adaptation
	  //
	  if (local_gaussian_index == gm_index) {
	    gm.getMean(mean);
	    break;
	  }
	  else
	    local_gaussian_index++;
	}
	else {
	  return Error::handle(name(), L"nodeScore",
	         RegressionDecisionTreeNode::ERR_ADAPT_NO_GAUSSIAN,
	    	       __FILE__, __LINE__);
	}
	
	// move to the next mixture component
	//
	more = models.gotoNext();
	  
      } // end of loop over all mixtures
    } // end of else if mixture model

    // in debug mode, print out the mean value and score for this node
    //
    if (debug_level_d >= Integral::ALL) {  
      mean.debug(L"mean=");
      Double(num_occ).debug(L"num_occ");
      Double(mean.distance(average_mean_d) * num_occ).debug(L"score");
    }
    
    cluster_score_d += mean.distance(average_mean_d) * num_occ;
    
  } // end of all model loop

  // exit gracefully
  //
  return true;
}

// method: updateDistribution
//
// arguments:
//  Vector<StatisticalModel>& stat_models: (input) vector of statistical model
//
// return: logical error status
//
// update the distribution for this node
//
bool8 RegressionDecisionTreeNode::updateDistribution(Vector<StatisticalModel>& stat_models_a) {

  if (!calcClusterDistribution(stat_models_a)) {
    return Error::handle(name(), L"updateDistribution",
			 ERR, __FILE__, __LINE__);
  }

  if (!nodeScore(stat_models_a)) {
    return Error::handle(name(), L"updateDistribution",
			 ERR, __FILE__, __LINE__);
  }
   
  // exit gracefully
  //
  return true;
}

// method: calcClusterDistribution
//
// arguments:
//  Vector<StatisticalModel>& stat_models: (input) vector of statistical model
//
// return: logical error status
//
// Calculate the cluster distribution, which includes pooled mean and
// pooled variance.
//
// reference Eq 7
// J. Zhao, et al, "Tutorial for Decision Tree-Based
// State-Tying For Acoustic Modeling", pp. 6, June, 1999.
//
// note that in this implementation we assume the mean as a row
// vector while the reference equation assumes the mean to be a column
// vector.
//
// |C| = determinant(C)
//
// C = (a / sum_num_occ) - b
//
// where a = sigma ( num_occ(s) ( cov(s) + transpose(mean(s)) * mean(s) ) )
//             s
//
//       sum_num_occ = sigma ( state_num_occ(s) )
//                       s
//
//       b = transpose(c) * c
//
//       c  = (sigma ( num_occ(s) * mean(s) )) / sum_num_occ ;
//               s
//
//       where s = statistical models at this node
// 
bool8 RegressionDecisionTreeNode::calcClusterDistribution(Vector<StatisticalModel>& stat_models_a) {

  VectorFloat mean, temp_mean, mean_sum;
  VectorFloat weights;
  MatrixFloat covar, cov_sum;
  
  float32 weight = 0.0;
  int32 i = 0;

  // local variables
  //
  float32 sum_num_occ;
  MatrixFloat a;
  MatrixFloat b;
  VectorFloat c;
  MatrixFloat pooled_covar;
  bool8 diagonal = true;
  bool8 res = true;
    
  // local variables
  //
  float64 datapoint_num_occ = 0.0;
  MatrixFloat mprod1;
  VectorFloat vprod1;
  MatrixFloat mprod2;
  MatrixFloat msum1;
    
  // loop all models in the node
  //
  for (bool8 more_index = gaussian_models_d.gotoFirst(); more_index;
       more_index = gaussian_models_d.gotoNext()) {

    // get the first datapoint in triple
    //
    RDataPoint* datapoint = gaussian_models_d.getCurr();
    
    // get the statistical model index
    //
    int32 sm_index = (int32)datapoint->first();
    
    // get the gaussian model index
    //
    int32 gm_index = (int32)datapoint->second();
	    
    // if the model to adapt is a single Gaussian model, compute its
    // contribution to the global adaptation
    //
    if (stat_models_a(sm_index).getType() == StatisticalModel::GAUSSIAN_MODEL) {
      GaussianModel& gm = stat_models_a(sm_index).getGaussianModel();
      
      gm.getMean(mean);
      gm.getCovariance(covar);      
      datapoint_num_occ = gm.getOccupancy();
      weight = 1.0;
    }
    
    // if the model to adapt is a mixture of Gaussian models, compute
    // the contribution of mixture component match the index to the
    // global adaptation
    //
    else if (stat_models_a(sm_index).getType() == StatisticalModel::MIXTURE_MODEL) {

      // get the linked list of mixture components
      //
      SingleLinkedList<StatisticalModel>& models =
	const_cast<StatisticalModel&>(stat_models_a(sm_index)).
	getMixtureModel().getModels();

      res =
      	const_cast<StatisticalModel&>(stat_models_a(sm_index)).
	getMixtureModel().getWeights(weights);
      
      // check whether each mixture component is Gaussian model and
      // compute its contribution to the global adaptation
      //
      bool8 more = models.gotoFirst();
      int32 local_gaussian_index = 0;
      
      while (more) {

	// get the current mixture component
	//
	StatisticalModel* model = models.getCurr();

	// check whether it is Gaussian
	//
	if (model->getType() == StatisticalModel::GAUSSIAN_MODEL) {
	  GaussianModel& gm = model->getGaussianModel();

	  // compute the contribution of this Gaussian model to the
	  // global adaptation
	  //
	  if (local_gaussian_index == gm_index) {
	    gm.getMean(mean);
	    gm.getCovariance(covar);

	    datapoint_num_occ = gm.getOccupancy();
	    weight = weights(gm_index);
	    break;
	  }
	  else
	    local_gaussian_index++;
	}
	else {
	  return Error::handle(name(), L"calcClusterDistribution",
		 RegressionDecisionTreeNode::ERR_ADAPT_NO_GAUSSIAN,
			       __FILE__, __LINE__);
	}
	
	// move to the next mixture component
	//
	more = models.gotoNext();
	  
      } // end of loop over all mixtures
    } // end of else if mixture model


    // check if the covariance is non-diagonal
    //
    if (!covar.isDiagonal()) {
      diagonal = false;
    }

    // intermediate computations to compute matrices "a" and "c"
    //
    res = mprod1.outerProduct(mean, mean);
    res = msum1.add(covar, mprod1);
    res = msum1.mult(datapoint_num_occ);
    a.setDimensions(msum1);
    res = a.add(msum1);
    res = vprod1.assign(mean);
    res = vprod1.mult(datapoint_num_occ);
    c.setLength(vprod1.length());
    res = c.add(vprod1);

    temp_mean.assign(mean);

    //    temp_mean.mult(weight);
    temp_mean.mult(datapoint_num_occ);
    
    if (i == 0) {
      mean_sum.assign(temp_mean);
      i++;
    }
    else {      
      mean_sum.add(temp_mean);
      i++;      
    }
       
  } // end of all model loop

  // get the sum of occupancies at this node
  //
  res = computeSumOccupancy(stat_models_a, sum_num_occ);

  //  compute the pooled covariance matrix and its determinant
  //
  res = a.div(sum_num_occ);
  res = c.div(sum_num_occ);
  res = b.outerProduct(c, c);
  res = pooled_covar.sub(a, b);
  
  // set the pooled covariance to diagonal if all the datapoints had
  // diagonal covariance. this is an assumptiom though the equations
  // don't generate a diagonal even if all the input data points have
  // diagonal covariances
  //
  if (diagonal) {

    // local variables
    //
    MatrixFloat temp;

    temp.setDimensions(pooled_covar);
    temp.setDiagonal(pooled_covar);
    average_covariance_d.assign(temp);
  }
  
  // now calculate mean and covariance for this cluster
  //
  average_mean_d.assign(mean_sum);
  average_mean_d.div(sum_num_occ);  
  
  cluster_accumulate_d = sum_num_occ;

  if (debug_level_d >= Integral::ALL) {
    average_mean_d.debug(L"average_mean_d");
    average_covariance_d.debug(L"average_covariance_d");
    cluster_accumulate_d.debug(L"cluster_accumulate_d");
  }
  
  // exit gracefully
  //
  return res;
}


// method: createTransform
//
// arguments:
//  Vector<StatisticalModel>& stat_models: (input) vector of statistical model
//
// return: logical error status
//
// Create the transform
//
bool8 RegressionDecisionTreeNode::createTransform(Vector<StatisticalModel>& stat_models_a) {

  // get the number of features
  //
  int32 num_feat = 0;

  // if the model to adapt is Gaussian model, get feature dimension
  //
  if (stat_models_a(0).getType()
      == StatisticalModel::GAUSSIAN_MODEL) {
    
    VectorFloat mean;
    stat_models_a(0).getGaussianModel().getMean(mean);
    num_feat = mean.length();
  }

  // if the model to adapt is mixture model, get feature dimension of
  // the first mixture component
  //
  else if (stat_models_a(0).getType()
	   == StatisticalModel::MIXTURE_MODEL) {
    
    VectorFloat mean;
    StatisticalModel* sm = stat_models_a(0).getMixtureModel()
      .getModels().getFirst();
    sm->getMean(mean);
    num_feat = mean.length();
  }
  else {
    return Error::handle(name(), L"createTransform",
	   RegressionDecisionTreeNode::ERR_ADAPT_NO_GAUSSIAN,
			 __FILE__, __LINE__);
  }

  // define Z matrix
  //
  MatrixFloat z_glob(num_feat, num_feat + 1);

  // define G matrix
  //
  Vector<MatrixDouble> g_glob(num_feat);
  for (int32 i = 0; i < num_feat; i++) {
    g_glob(i).setDimensions(num_feat + 1, num_feat + 1);
  }

  float32 sum_num_occ;
  
  // get the sum of occupancies at this node
  //
  computeSumOccupancy(stat_models_a, sum_num_occ);

  if (debug_level_d >= Integral::ALL) {
    Float(sum_num_occ).debug(L"sum_num_occ on this node");
  }
  
  if (sum_num_occ < 1000.0) {
    if (debug_level_d >= Integral::ALL) {
      Float(sum_num_occ).debug(L"Not enough occupation acount for creating transformation");
    }
    transform_flag_d = false;
    return false;
  }
  
  transform_flag_d = true;
  
  // loop over all gaussian, but not silence model and
  // compute their contributions to the global adaptation
  //
  for (bool8 more_index = gaussian_models_d.gotoFirst(); more_index;
       more_index = gaussian_models_d.gotoNext()) {

    // get the first datapoint in triple
    //
    RDataPoint* datapoint = gaussian_models_d.getCurr();

    // get the statistical model index
    //
    int32 sm_index = (int32)datapoint->first();

    // get the gaussian model index
    //
    int32 gm_index = (int32)datapoint->second();

    // if the model to adapt is a single Gaussian model, compute its
    // contribution to the global adaptation
    //
    if (stat_models_a(sm_index).getType() == StatisticalModel::GAUSSIAN_MODEL) {
      GaussianModel& gm = stat_models_a(sm_index).getGaussianModel();

      // compute the contribution of this Gaussian model to the global
      // adaptation
      //
      adaptPart(g_glob, z_glob, gm);
    }
    
    // if the model to adapt is a mixture of Gaussian models, compute
    // the contribution of mixture component match the index to the
    // global adaptation
    //
    else if (stat_models_a(sm_index).getType() == StatisticalModel::MIXTURE_MODEL) {

      // get the linked list of mixture components
      //
      SingleLinkedList<StatisticalModel>& models =
	const_cast<StatisticalModel&>(stat_models_a(sm_index)).
	getMixtureModel().getModels();

      // check whether each mixture component is Gaussian model and
      // compute its contribution to the global adaptation
      //
      bool8 more = models.gotoFirst();
      int32 local_gaussian_index = 0;
      
      while (more) {

	// get the current mixture component
	//
	StatisticalModel* model = models.getCurr();

	// check whether it is Gaussian
	//
	if (model->getType() == StatisticalModel::GAUSSIAN_MODEL) {
	  GaussianModel& gm = model->getGaussianModel();

	  // compute the contribution of this Gaussian model to the
	  // global adaptation
	  //
	  if (local_gaussian_index == gm_index) {
	    adaptPart(g_glob, z_glob, gm);
	    break;
	  }
	  else
	    local_gaussian_index++;
	}
	else {
	  return Error::handle(name(), L"createTransform",
		 RegressionDecisionTreeNode::ERR_ADAPT_NO_GAUSSIAN,
			       __FILE__, __LINE__);
	}
	
	// move to the next mixture component
	//
	more = models.gotoNext();
	  
      } // end of loop over all mixtures
    } // end of else if mixture model

  } // end of all model loop
  
  // define transformation matrix W (one row for each Gaussian
  // dimension)
  //
  w_transform_d.setLength(num_feat);
  
  for (int32 i = 0; i < num_feat; i++) {
    
    // since transformation is applied to the extended mean
    //
    w_transform_d(i).setLength(num_feat + 1);
  }

  // compute the inverse of G
  //
  Vector<MatrixDouble> g_inv(num_feat);

  for (int32 i = 0; i < num_feat; i++) {

    float64 det = g_glob(i).determinant();
    Double Det = det;

    // if in debug mode, print the determinant of G
    //
    if (debug_level_d >= Integral::ALL) {
      String out;
      out.concat(i);
      out.concat(L". det:");
      out.concat(Det);
      Console::put(out);
    }

    // compute the threshold for a singularity of G matrix
    //
    float64 threshold = Integral::pow(0.1, num_feat);

    // if the matrix is singular with a given threshold, compute
    // inverse using SVD
    //
    if (g_glob(i).isSingular(threshold)) {

      MatrixDouble uu;
      MatrixDouble ww;
      MatrixDouble vv;
      g_glob(i).decompositionSVD(uu, ww, vv);

      // compute the inverse of ww
      //
      for (int32 j = 0; j < num_feat; j++) {
	Double diag_element = ww(j, j);
	
	// replace diagonal elements that are close to zero by zero
	//
	if (diag_element.almostEqual(0)) {
	  ww.setValue(j, j, 0);
	}
	
	// otherwise compute inverted value
	//
	else {
	  ww.setValue(j, j, 1 / ww(j, j));
	}
      }

      // finally compute the inverse of G
      //
      g_inv(i).mult(vv, ww);
      uu.transpose();
      g_inv(i).mult(uu);
    }

    // if matrix is not singular with a given threshold, compute
    // inverse by a standard method
    //
    else {
      g_inv(i).inverse(g_glob(i));
    }
  }

  // compute the rows of Z matrix
  //
  Vector<VectorFloat> z_r(num_feat);
  for (int32 i = 0; i < num_feat; i++) {
    z_glob.getRow(z_r(i), i);
  }

  // convert G_inv matrices to float32 by looping over all elements
  //
  Vector<MatrixFloat> g_inv_f(num_feat);
  for (int32 i = 0; i < num_feat; i++) {
    g_inv_f(i).setDimensions(num_feat + 1, num_feat + 1);
    for (int32 j = 0; j < num_feat + 1; j++) {
      for (int32 k = 0; k < num_feat + 1; k++) {
	g_inv_f(i).setValue(j, k, g_inv(i)(j, k));
      }
    }
  }
  
  // compute rows of the transformation matrix W
  //
  for (int32 i = 0; i < num_feat; i++) {
    
    g_inv_f(i).multv(w_transform_d(i), z_r(i));

    // if in debug mode, print the rows of the transformation matrix W
    //
    if (debug_level_d >= Integral::ALL) {
      String comment(L"w(i), i: ");
      comment.concat(i);
      w_transform_d(i).debug(comment);
    }
  }

  w_transform_d.debug(L"transform");
  
  // exit gracefully
  //
  return true;  
}


// method: adaptPart
//
// arguments:
//  Vector<MatrixDouble>& g: (output) vector of G matrices
//  MatrixFloat& z: (output) Z matrix
//  GaussianModel& gm: (input) Gaussian model with its statistics
//
// return: a bool8 value indicating status
//
// this method computes contribution of one Gaussian model to
// cumulative adaptation matrices G and Z
//
bool8 RegressionDecisionTreeNode::adaptPart(Vector<MatrixDouble>& g_a,
					      MatrixFloat& z_a,
					      GaussianModel& gm_a) {

  // print the access count in debug mode
  //
  if (debug_level_d >= Integral::ALL) {
    String out(L"access count:");
    out.concat(gm_a.getAccessCount());
    Console::put(out);
  }

  // if model was accessed, compute its contribution to the global
  // adaptation and add it to the global adaptation matrices
  //
  if (gm_a.getAccessCount() > 0) {
  
    // form the extended mean
    //
    VectorFloat extended_mean;
    gm_a.getMean(extended_mean);

    int32 num_feat = extended_mean.length();
    extended_mean.setLength(num_feat + 1);
    for (int32 i = num_feat; i > 0; i--) {
      extended_mean(i) = extended_mean(i - 1);
    }
    extended_mean(0) = 1;

    // inverse the covariance matrix
    //
    MatrixFloat inv_cov;
    gm_a.getCovariance(inv_cov);
    inv_cov.inverse();
      
    // convert mean accumulator to float32
    //
    VectorDouble mean_accum;
    gm_a.getMeanAccumulator(mean_accum);
    VectorFloat mean_accum_f(num_feat);
    for (int32 i = 0; i < num_feat; i++) {
      mean_accum_f(i) = (Float)mean_accum(i);
    }

    // compute contribution of this model to the matrix Z
    //
    VectorFloat tmp;
    inv_cov.multv(tmp, mean_accum_f);

    MatrixFloat z_r;
    z_r.outerProduct(tmp, extended_mean);
    z_a.add(z_r);

    // compute V matrix
    //
    MatrixFloat v(inv_cov);
    Double occup_accum = gm_a.getOccupancy();
    v.mult((float32)occup_accum);

    // compute D matrix
    //
    MatrixFloat d;
    d.outerProduct(extended_mean, extended_mean);

    // check for the diagonal covariance before computing G matrix
    //
    if (!inv_cov.isDiagonal()) {
      return Error::handle(name(),
			   L"adapt - only diagonal matrices supported",
			   Error::ARG, __FILE__, __LINE__);
    }
    
    // compute contribution of this model to the matrix G
    //
    MatrixFloat g_r;
    for (int32 i = 0; i < num_feat; i++) {
      g_r.assign(d);
      g_r.mult(v(i, i));
      g_a(i).add(g_r);
    }

    // exit gracefully
    //
    return true;
  }

  // else the model was not accessed
  //
  else {
    
    // exit ungracefully
    //
    return false;
  }
}


// method: containModel
//
// arguments:
//  int32 sm_index: (input) statistical model index
//  int32 gm_index: (input) Gaussian model index
//
// return: logical error status
//
// To test if containing the models specified index in this node
//
bool8 RegressionDecisionTreeNode::containModel(int32 sm_index_a, int32 gm_index_a) {

  // loop all models in the node
  //
  for (bool8 more_index = gaussian_models_d.gotoFirst(); more_index;
       more_index = gaussian_models_d.gotoNext()) {

    // get the first datapoint in triple
    //
    RDataPoint* datapoint = gaussian_models_d.getCurr();
    
    // get the statistical model index
    //
    int32 sm_index = (int32)datapoint->first();
    
    // get the gaussian model index
    //
    int32 gm_index = (int32)datapoint->second();

    if ((sm_index_a == sm_index) && (gm_index_a == gm_index)) {
      return true;
    }
  }
    
  // exit gracefully
  //
  return false;
}


// method: computeSumOccupancy
//
// arguments:
//  Vector<StatisticalModel>& stat_models: (input) vector of statistical model
//  float32& sum_num_occ: (output) the sum of occupancy 
//
// return: a bool8 value indicating status
//
// this method sums the occupancies of all the datapoints.
//
bool8 RegressionDecisionTreeNode::computeSumOccupancy(Vector<StatisticalModel>& stat_models_a, float32& sum_num_occ_a) {
  
  // local variables
  //
  float64 datapoint_num_occ = 0.0;  
  VectorFloat weights;
  sum_num_occ_a = (float64)0;
  float64 weight = 0;
  bool8 res = true;
  
  // loop all models in the node
  //
  for (bool8 more_index = gaussian_models_d.gotoFirst(); more_index;
       more_index = gaussian_models_d.gotoNext()) {

    // get the first datapoint in triple
    //
    RDataPoint* datapoint = gaussian_models_d.getCurr();
    
    // get the statistical model index
    //
    int32 sm_index = (int32)datapoint->first();

    // get the gaussian model index
    //
    int32 gm_index = (int32)datapoint->second();
    
    // if the model to adapt is a single Gaussian model, compute its
    // contribution to the global adaptation
    //
    if (stat_models_a(sm_index).getType() == StatisticalModel::GAUSSIAN_MODEL) {
      GaussianModel& gm = stat_models_a(sm_index).getGaussianModel();
      
      datapoint_num_occ = gm.getOccupancy();
      weight = 1.0;
    }
    
    // if the model to adapt is a mixture of Gaussian models, compute
    // the contribution of mixture component match the index to the
    // global adaptation
    //
    else if (stat_models_a(sm_index).getType() == StatisticalModel::MIXTURE_MODEL) {

      // get the linked list of mixture components
      //
      SingleLinkedList<StatisticalModel>& models =
	const_cast<StatisticalModel&>(stat_models_a(sm_index)).
	getMixtureModel().getModels();

      res =
      	const_cast<StatisticalModel&>(stat_models_a(sm_index)).
	getMixtureModel().getWeights(weights);
      
      // check whether each mixture component is Gaussian model and
      // compute its contribution to the global adaptation
      //
      bool8 more = models.gotoFirst();
      int32 local_gaussian_index = 0;
      
      while (more) {

	// get the current mixture component
	//
	StatisticalModel* model = models.getCurr();

	// check whether it is Gaussian
	//
	if (model->getType() == StatisticalModel::GAUSSIAN_MODEL) {
	  GaussianModel& gm = model->getGaussianModel();

	  // compute the contribution of this Gaussian model to the
	  // global adaptation
	  //
	  if (local_gaussian_index == gm_index) {
	    datapoint_num_occ = gm.getOccupancy();
	    weight = weights(gm_index);
	    break;
	  }
	  else
	    local_gaussian_index++;
	}
	else {
	  return Error::handle(name(), L"calcClusterDistribution",
		 RegressionDecisionTreeNode::ERR_ADAPT_NO_GAUSSIAN,
			       __FILE__, __LINE__);
	}
	
	// move to the next mixture component
	//
	more = models.gotoNext();
	  
      } // end of loop over all mixtures
    } // end of else if mixture model

    // add this occupancy to the sum
    //
    sum_num_occ_a += datapoint_num_occ;

    Double(sum_num_occ_a).debug(L"sum_num_occ_a");
    
  } // end of all model loop

  // exit gracefully
  //
  return res;
}