Free cookie consent management tool by TermsFeed Policy Generator

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

Last change on this file since 10662 was 9562, checked in by gkronber, 12 years ago

#1967 worked on Gaussian process evolution.

File size: 3.0 KB
Line 
1// This is the main DLL file.
2
3#include "stdafx.h"
4
5#include <iostream>
6#include "HeuristicLab.Eigen.h"
7#include "Eigen/Dense"
8#include "Eigen/Core"
9using namespace Eigen;
10using namespace std;
11
12namespace HeuristicLabEigen {
13  void MyEigen::Solve(double a[], double x[], double alpha[], double invL[], double sqrSigmaNoise, int n, double* negLogLikelihood, int* computationInfo) {
14    MatrixXd A = Map<Matrix<double, Dynamic, Dynamic,ColMajor>>(a, n, n);
15    LLT<MatrixXd> chol = A.llt();
16    *computationInfo = chol.info();
17    if(chol.info() != Success) return;
18
19    VectorXd res = chol.solve(Map<VectorXd>(x, n)) / sqrSigmaNoise;
20
21    for(int i=0;i<n;i++) {
22      alpha[i] = res(i);
23    }
24
25    *negLogLikelihood = 0.5 * Map<VectorXd>(x, n).dot(res)
26      + chol.matrixU().toDenseMatrix().diagonal<0>().array().log().sum()
27      + (n / 2.0) * log(2.0 * M_PI * sqrSigmaNoise);
28
29    MatrixXd invLRes = chol.solve(MatrixXd::Identity(n,n)) / sqrSigmaNoise - res * res.transpose();
30    for(int i=0;i<n;i++) {
31      for(int j=0;j<n;j++) {
32        invL[j + i * n] = invLRes(i,j);
33      }
34    }
35  }
36
37
38
39  //double diagSum = Enumerable.Range(0, n).Select(i => Math.Log(chol[i * i])).Sum();
40
41  //// solve for alpha
42  //ILArray<double> ym = ILMath.array(y.Zip(m, (a, b) => a - b).ToArray());
43
44  //MatrixProperties lowerTriProps = MatrixProperties.LowerTriangular;
45  //MatrixProperties upperTriProps = MatrixProperties.UpperTriangular;
46  //ILArray<double> alpha;
47  //try {
48  //  ILArray<double> t1 = ILMath.linsolve(chol.T, ym, ref lowerTriProps);
49  //  alpha = ILMath.linsolve(chol, t1, ref upperTriProps) / sqrSigmaNoise;
50  //}
51  //catch (Exception e) {
52  //  throw new ArgumentException("covariance matrix is not positive definite", e);
53  //}
54  //if (alpha == null || alpha.IsEmpty) throw new ArgumentException();
55  //this.alpha = new double[alpha.Length];
56  //alpha.ExportValues(ref this.alpha);
57
58  //negativeLogLikelihood = 0.5 * (double)ILMath.multiply(ym.T, alpha) + diagSum + (n / 2.0) * Math.Log(2.0 * Math.PI * sqrSigmaNoise);
59
60  //// derivatives
61  //int nAllowedVariables = x.GetLength(1);
62
63  //ILArray<double> lCopy;
64  //try {
65  //  ILArray<double> t2 = ILMath.linsolve(chol.T, ILMath.eye<double>(n, n), ref lowerTriProps);
66  //  lCopy = ILMath.linsolve(chol, t2, ref upperTriProps) / sqrSigmaNoise - ILMath.multiply(alpha, alpha.T);
67  //}
68  //catch (Exception e) {
69  //  throw new ArgumentException("covariance matrix is not positive definite", e);
70  //}
71
72  //double noiseGradient = sqrSigmaNoise * ILMath.sumall(ILMath.diag(lCopy)).GetValue(0, 0);
73
74  //double[] meanGradients = new double[meanFunction.GetNumberOfParameters(nAllowedVariables)];
75  //for (int k = 0; k < meanGradients.Length; k++) {
76  //  ILArray<double> meanGrad =
77  //    ILMath.array(
78  //      Enumerable.Range(0, alpha.Length)
79  //        .Select(r => mean.Gradient(x, r, k))
80  //        .ToArray());
81  //  meanGradients[k] = -(double)ILMath.multiply(meanGrad.T, alpha);
82  //}
83
84}
85
Note: See TracBrowser for help on using the repository browser.