source: branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Eigen/HeuristicLab.Eigen.cpp @ 9674

// This is the main DLL file.

#include "stdafx.h"

#include <iostream>
#include "HeuristicLab.Eigen.h"
#include "Eigen/Dense"
#include "Eigen/Core"
using namespace Eigen;
using namespace std;

namespace HeuristicLabEigen {
  void MyEigen::Solve(double a[], double x[], double alpha[], double invL[], double sqrSigmaNoise, int n, double* negLogLikelihood, int* computationInfo) {
    MatrixXd A = Map<Matrix<double, Dynamic, Dynamic,ColMajor>>(a, n, n);
    LLT<MatrixXd> chol = A.llt();
    *computationInfo = chol.info();
    if(chol.info() != Success) return;

    VectorXd res = chol.solve(Map<VectorXd>(x, n)) / sqrSigmaNoise;

    for(int i=0;i<n;i++) {
      alpha[i] = res(i);
    }

    *negLogLikelihood = 0.5 * Map<VectorXd>(x, n).dot(res) 
      + chol.matrixU().toDenseMatrix().diagonal<0>().array().log().sum() 
      + (n / 2.0) * log(2.0 * M_PI * sqrSigmaNoise);

    MatrixXd invLRes = chol.solve(MatrixXd::Identity(n,n)) / sqrSigmaNoise - res * res.transpose();
    for(int i=0;i<n;i++) {
      for(int j=0;j<n;j++) {
        invL[j + i * n] = invLRes(i,j);
      }
    }
  }



  //double diagSum = Enumerable.Range(0, n).Select(i => Math.Log(chol[i * i])).Sum();

  //// solve for alpha
  //ILArray<double> ym = ILMath.array(y.Zip(m, (a, b) => a - b).ToArray());

  //MatrixProperties lowerTriProps = MatrixProperties.LowerTriangular;
  //MatrixProperties upperTriProps = MatrixProperties.UpperTriangular;
  //ILArray<double> alpha;
  //try {
  //  ILArray<double> t1 = ILMath.linsolve(chol.T, ym, ref lowerTriProps);
  //  alpha = ILMath.linsolve(chol, t1, ref upperTriProps) / sqrSigmaNoise;
  //}
  //catch (Exception e) {
  //  throw new ArgumentException("covariance matrix is not positive definite", e);
  //}
  //if (alpha == null || alpha.IsEmpty) throw new ArgumentException();
  //this.alpha = new double[alpha.Length];
  //alpha.ExportValues(ref this.alpha);

  //negativeLogLikelihood = 0.5 * (double)ILMath.multiply(ym.T, alpha) + diagSum + (n / 2.0) * Math.Log(2.0 * Math.PI * sqrSigmaNoise);

  //// derivatives
  //int nAllowedVariables = x.GetLength(1);

  //ILArray<double> lCopy;
  //try {
  //  ILArray<double> t2 = ILMath.linsolve(chol.T, ILMath.eye<double>(n, n), ref lowerTriProps);
  //  lCopy = ILMath.linsolve(chol, t2, ref upperTriProps) / sqrSigmaNoise - ILMath.multiply(alpha, alpha.T);
  //}
  //catch (Exception e) {
  //  throw new ArgumentException("covariance matrix is not positive definite", e);
  //}

  //double noiseGradient = sqrSigmaNoise * ILMath.sumall(ILMath.diag(lCopy)).GetValue(0, 0);

  //double[] meanGradients = new double[meanFunction.GetNumberOfParameters(nAllowedVariables)];
  //for (int k = 0; k < meanGradients.Length; k++) {
  //  ILArray<double> meanGrad =
  //    ILMath.array(
  //      Enumerable.Range(0, alpha.Length)
  //        .Select(r => mean.Gradient(x, r, k))
  //        .ToArray());
  //  meanGradients[k] = -(double)ILMath.multiply(meanGrad.T, alpha);
  //}

}