import cern.jet.math.Functions; //导入依赖的package包/类
/**
 * Normalizes matrix of p(z|u) such that \forall_u: \sum_z p(z|u) = 1.
 *
 * @param pu_z normalized matrix of p(z|u)
 */
@Override
protected void normalizePuz(DoubleMatrix2D pu_z) {
    for (int u = 0; u < pu_z.rows(); u++) {
        DoubleMatrix1D tmp = pu_z.viewRow(u);
        double norm = tmp.aggregate(Functions.plus, Functions.identity);
        if (norm != 0.0) {
            tmp.assign(Functions.mult(1 / norm));
        }
    }
}
 

import cern.jet.math.Functions; //导入依赖的package包/类
/**
 * Normalizes matrix of p(i|z) such that \forall_z: \sum_i p(i|z) = 1.
 *
 * @param piz normalized matrix of p(i|z)
 */
@Override
protected void normalizePiz(DoubleMatrix2D piz) {
    for (int i = 0; i < piz.columns(); i++) {
        DoubleMatrix1D tmp = piz.viewColumn(i);
        double norm = tmp.aggregate(Functions.plus, Functions.identity);
        if (norm != 0.0) {
            tmp.assign(Functions.mult(1 / norm));
        }
    }
}
 

import cern.jet.math.Functions; //导入依赖的package包/类
/**
 * Update the initial state probabilities. 
 * Calculations done in log-space.
 */
private void updatePi() {
  //set up a matrix whose columns will hold the gamma(z_0) vector of each
  //analysis region
  DoubleMatrix2D logGamma_zeros = new DenseDoubleMatrix2D(hmm_numStates, analysisRegions.size());
  for (int i = 0; i < analysisRegions.size(); i++) {
    HMMAnalysisRegion region = analysisRegions.get(i);
    logGamma_zeros.viewColumn(i).assign(region.getLogGamma().viewRow(0));      
  }
  
  //the sums along each of the rows will be the numerators for the update eqs  
  DoubleMatrix1D hmm_LogPi_new = new DenseDoubleMatrix1D(hmm_numStates);
  for (int k = 0; k < hmm_numStates; k++) {
    hmm_LogPi_new.setQuick(k, Numerical.log_sum_exp(logGamma_zeros.viewRow(k)));
  }
  
  //summing all of those numerators will give the denominator for the update eq
  double logGammaSum = Numerical.log_sum_exp(hmm_LogPi_new);
  hmm_LogPi_new.assign(Functions.minus(logGammaSum));
  
  hmm_LogPi = hmm_LogPi_new;
}
 

import cern.jet.math.Functions; //导入依赖的package包/类
/**
 * 
 * @param n => (n_max x N) Matrix where each row is state of reservoir for each time step
 * @param y => (n_max x L) Matrix where each row is state of output for each time step
 * @param alpha The idenity coefficient. This can be found using cross-validation
 * @return
 */
public static DoubleMatrix2D ridge(DoubleMatrix2D n, DoubleMatrix2D y, double alpha){
	DoubleFactory2D factory = DoubleFactory2D.dense;

	//force same type, this is giving a cast error. Bug on their part?
	DoubleMatrix2D _n = factory.make(n.toArray());
	DoubleMatrix2D _y = factory.make(y.toArray());

	
	//aI
	DoubleMatrix2D scaledIdentity = factory.identity(_n.rows()).assign(Functions.mult(alpha));
	
	
	Algebra alg = new Algebra();
	// ((n x n^T) + aI)^-1
	// (n_max x N) x (N x n_max)
	DoubleMatrix2D inverse = alg.inverse(n.zMult(_n.viewDice(), null).assign(scaledIdentity, Functions.plus));
	//System.out.println("-1=" + inverse);

	//X = ((n x n^T) + aI)^-1 x n
	// (n_max x n_max) x 
	DoubleMatrix2D X = inverse.zMult(_n, null);
	//System.out.println("X=" + X);
	//((n x n^T) + aI)^-1 x n x y
	DoubleMatrix2D weights = X.viewDice().zMult(_y, null);
	//X.viewDice()
	return weights;	
}
 

import cern.jet.math.Functions; //导入依赖的package包/类
public static double log_sum_exp(DoubleMatrix1D exponents) {
  double maxExp = exponents.aggregate(Functions.max, Functions.identity);
  if (maxExp == Double.NEGATIVE_INFINITY) {
    return Double.NEGATIVE_INFINITY;
  }
  else {
    double expSum = exponents.aggregate(Functions.plus, Functions.chain(Functions.exp, Functions.minus(maxExp)));
    return maxExp + Math.log(expSum);
  }
}
 

import cern.jet.math.Functions; //导入依赖的package包/类
public Matrix plus(double value) {
	Matrix result = new ColtSparseDoubleMatrix2D((SparseDoubleMatrix2D) matrix.copy().assign(Functions.plus(value)));
	MapMatrix<String, Object> a = getMetaData();
	if (a != null) {
		result.setMetaData(a.clone());
	}
	return result;
}
 

import cern.jet.math.Functions; //导入依赖的package包/类
public Matrix times(double value) {
	Matrix result = new ColtSparseDoubleMatrix2D((SparseDoubleMatrix2D) matrix.copy().assign(Functions.mult(value)));
	MapMatrix<String, Object> a = getMetaData();
	if (a != null) {
		result.setMetaData(a.clone());
	}
	return result;
}
 

import cern.jet.math.Functions; //导入依赖的package包/类
public Matrix divide(double value) {
	Matrix result = new ColtSparseDoubleMatrix2D((SparseDoubleMatrix2D) matrix.copy().assign(Functions.div(value)));
	MapMatrix<String, Object> a = getMetaData();
	if (a != null) {
		result.setMetaData(a.clone());
	}
	return result;
}
 

import cern.jet.math.Functions; //导入依赖的package包/类
public Matrix minus(double value) {
	Matrix result = new ColtSparseDoubleMatrix2D((SparseDoubleMatrix2D) matrix.copy()
			.assign(Functions.minus(value)));
	MapMatrix<String, Object> a = getMetaData();
	if (a != null) {
		result.setMetaData(a.clone());
	}
	return result;
}
 

import cern.jet.math.Functions; //导入依赖的package包/类
public Matrix plus(double value) {
	Matrix result = new ColtDenseDoubleMatrix2D((DenseDoubleMatrix2D) matrix.copy().assign(Functions.plus(value)));
	MapMatrix<String, Object> a = getMetaData();
	if (a != null) {
		result.setMetaData(a.clone());
	}
	return result;
}
 

import cern.jet.math.Functions; //导入依赖的package包/类
public Matrix minus(double value) {
	Matrix result = new ColtDenseDoubleMatrix2D((DenseDoubleMatrix2D) matrix.copy().assign(Functions.minus(value)));
	MapMatrix<String, Object> a = getMetaData();
	if (a != null) {
		result.setMetaData(a.clone());
	}
	return result;
}
 

import cern.jet.math.Functions; //导入依赖的package包/类
public Matrix times(double value) {
	Matrix result = new ColtDenseDoubleMatrix2D((DenseDoubleMatrix2D) matrix.copy().assign(Functions.mult(value)));
	MapMatrix<String, Object> a = getMetaData();
	if (a != null) {
		result.setMetaData(a.clone());
	}
	return result;
}
 

import cern.jet.math.Functions; //导入依赖的package包/类
public Matrix divide(double value) {
	Matrix result = new ColtDenseDoubleMatrix2D((DenseDoubleMatrix2D) matrix.copy().assign(Functions.div(value)));
	MapMatrix<String, Object> a = getMetaData();
	if (a != null) {
		result.setMetaData(a.clone());
	}
	return result;
}
 

import cern.jet.math.Functions; //导入依赖的package包/类
/**
  * Compute the stock portfolio return standard deviation (sigma, a number).
  * 
  * A. Compute the stock portfolio covariance matrix:
  *   1. Calculate log-returns: aR(t) = log(p(t)/p(t-1)).
  *   2. Estimate the log-returns: maR(t) = d*aR(t) + (1-d)*maR(t-1) where d is memory.
  *   3. Estimate the covariance matrix:
  *      Cov(t) = d*(outer(aR(t)-maR(t),aR(t)-maR(t)) + (1-d)*Cov(t-1).
  * 
  * B. Compute the deviation of portfolio returns:
  *   1. Variance sigma**2 of the stock portfolio : sigma**2 = w**T.Cov.w 
  *   2. Assume, for simplicity, that loans have same variance.
  *   3. Return sqrt(sigma**2 + (1-sum(w))**2*s**2).
  */
static double CalculatePortfolioStandardDeviation(
   int numberOfFirmsToInvestIn,
   List<Double> firmMarketValuesNow,
   List<Double> firmMarketValuesLast,
   double[] firmStockWeights,
   double covarianceAdaptationRate,
   DoubleMatrix1D meanLogStockReturns,
   DoubleMatrix2D meanStockCovarianceMatrix
   ) {
   Algebra algebra = new Algebra();
   DenseDoubleMatrix1D logOfStockReturns = new DenseDoubleMatrix1D(numberOfFirmsToInvestIn);
   for (int i = 0; i < numberOfFirmsToInvestIn; ++i) {
      double presentOverLastValue = (firmMarketValuesNow.get(i) / firmMarketValuesLast.get(i));
      if (presentOverLastValue > 0.) 
         logOfStockReturns.set(i, Math.log(presentOverLastValue));
      else logOfStockReturns.set(i, 0.);
   }
   
   meanLogStockReturns.assign(
      logOfStockReturns.copy().assign(Functions.mult(covarianceAdaptationRate))
      .assign(meanLogStockReturns.copy().assign(Functions.mult(1. - covarianceAdaptationRate)), 
         Functions.plus));
   
   logOfStockReturns.assign(meanLogStockReturns, Functions.minus);
   DoubleMatrix2D newCovMatrix = 
      DoubleFactory2D.dense.make(numberOfFirmsToInvestIn, numberOfFirmsToInvestIn);
   algebra.multOuter(logOfStockReturns, logOfStockReturns, newCovMatrix);
   
   meanStockCovarianceMatrix.assign(
      newCovMatrix.copy().assign(Functions.mult(covarianceAdaptationRate))
      .assign(meanStockCovarianceMatrix.copy().assign(
         Functions.mult(1. - covarianceAdaptationRate)), Functions.plus));
   
   DoubleMatrix1D stockWeightsVector =
      new DenseDoubleMatrix1D(firmStockWeights);
   final double
      loanWeight = (1. - stockWeightsVector.zSum()),
      sigmaSquared = algebra.mult(stockWeightsVector,
         algebra.mult(meanStockCovarianceMatrix, stockWeightsVector));
   
   return Math.sqrt(sigmaSquared * (1. + loanWeight * loanWeight));
}
 

import cern.jet.math.Functions; //导入依赖的package包/类
private void checkLikelihood(DoubleMatrix2D logAlpha, DoubleMatrix2D logBeta) {
	double logLikelihood = Numerical.log_sum_exp(logAlpha.viewRow(0).copy().assign(logBeta.viewRow(0), Functions.plus));
	TreeMap<Integer, Double> mismatches = new TreeMap<Integer, Double>();
	for (int m = 1; m < logAlpha.rows(); m++) {
		double currentLL = Numerical.log_sum_exp(logAlpha.viewRow(m).copy().assign(logBeta.viewRow(m), Functions.plus));
		for (int j = 0; j < hmm_numStates; j++) {
			if (isStateFixedDuration[j] && (stateDurations[j] > 1)) {
				int start = Math.max(0, m - stateDurations[j] + 1);
				for (int n = start; n < m; n++) {
					currentLL = Numerical.log_add(currentLL, (logAlpha.getQuick(n, j) + logBeta.getQuick(n,j)));
				}
			}
		}
		if (Math.abs(logLikelihood - currentLL) > (double)1E-6) {
			mismatches.put(m, currentLL);
		}
	}
	
	if (!mismatches.isEmpty()) {
		StringBuffer errormsg = new StringBuffer();
		errormsg.append("Log-Likelihood calculations from forward-backward algorithm are NaN or inconsistent\n" + "0: " + logLikelihood);
		for (int pos : mismatches.keySet()) {
			errormsg.append(pos + ": " + mismatches.get(pos) + "\n");
		}
		logger.error(errormsg.toString());
		System.exit(1);
	}
}
 

import cern.jet.math.Functions; //导入依赖的package包/类
public static double matrixMax(DoubleMatrix1D mat) {
  return mat.aggregate(Functions.max, Functions.identity);
}
 

import cern.jet.math.Functions; //导入依赖的package包/类
@Override
public double calculate(final double[] parameters) {
    // 20 fitness values
    DoubleMatrix1D f = DoubleFactory1D.dense.make(20);

    // 19 parameters are estimated, passed in parameters[]
    for (int i = 0; i < parameters.length; i++) f.setQuick(i, parameters[i]);

    // One is fixed
    f.setQuick(19, Constants.FITNESS_FIXED_FOR_RELATIVE);

    double mean = f.zSum() / 20.0;
    f.assign(Functions.minus(mean));
    return -coefficient * f.aggregate(Functions.plus, Functions.pow(2.0));
}
 

import cern.jet.math.Functions; //导入依赖的package包/类
@Override
public double calculate(final double[] fitness) {
    DoubleMatrix1D theta = DoubleFactory1D.dense.make(fitness);
    theta.assign(Functions.exp);
    theta.assign(Functions.div(1 + theta.zSum()));


    // original:
    return ((alpha - 1) * theta.aggregate(Functions.plus, Functions.log))
            + ((alpha - 1) * Math.log(1 - theta.zSum()))
            + Math.log(Math.abs(Algebra.DEFAULT.det(jacobian(theta))));

}
 

import cern.jet.math.Functions; //导入依赖的package包/类
private double calculateMdR(final double[] fitness) {
    DoubleMatrix1D theta = DoubleFactory1D.dense.make(20, 0.0);
    for (int i = 0; i < fitness.length; i++) theta.setQuick(i, fitness[i]);

    theta.assign(Functions.exp);
    theta.assign(Functions.div(theta.zSum()));

    return alpha * theta.aggregate(Functions.plus, Functions.log);
}
 

import cern.jet.math.Functions; //导入依赖的package包/类
public static DoubleMatrix2D sprand(int m, int n, double density){
	int nnzwanted = (int)Math.round(m * n * Math.min(density, 1));
	DoubleFactory2D factory = DoubleFactory2D.sparse;
	DoubleMatrix2D i = factory.random(nnzwanted, 1).assign(customRandomizer).assign(Functions.mult(m)).assign(Functions.floor);//.assign(Functions.plus(1));
	DoubleMatrix2D j = factory.random(nnzwanted, 1).assign(customRandomizer).assign(Functions.mult(n)).assign(Functions.floor);//.assign(Functions.plus(1));

	DoubleMatrix2D rows=null;
	DoubleMatrix2D cols=null;
	//Combine i j
	DoubleMatrix2D ij = unique(factory.appendColumns(i, j));
	if (ij.cardinality() != 0){
		//separate into i and j again
		rows = ij.viewSelection(null, new int[]{0});
		cols= ij.viewSelection(null, new int[]{1});
	}
	
	
	DoubleMatrix2D rands = factory.random(rows.rows(), 1);
	
	double [][] R = new double [m][n];
	for (int x=0; x<rows.rows(); x++ ){
	//	for (int y=0; y<n; y++){
			R[(int)rows.get(x, 0)][(int)cols.get(x, 0)] = rands.get(x, 0);
	//	}
	}
	
	DoubleMatrix2D y = new SparseDoubleMatrix2D(R);
	return y;
}