Java tutorial
/******************************************************************************* * Copyright (c) 2014, 2015 IBM Corporation * * 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. *******************************************************************************/ package hulo.localization.models.obs; import hulo.localization.kernels.KernelFunction; import hulo.localization.likelihood.LikelihoodModel; import hulo.localization.thread.ExecutorServiceHolder; import hulo.localization.utils.JSONUtils; import hulo.localization.utils.StatUtils; import java.util.concurrent.ExecutionException; import java.util.concurrent.ExecutorService; import java.util.concurrent.Future; import org.apache.commons.math3.linear.LUDecomposition; import org.apache.commons.math3.linear.MatrixUtils; import org.apache.commons.math3.linear.RealMatrix; import org.apache.wink.json4j.JSONArray; import org.apache.wink.json4j.JSONException; import org.apache.wink.json4j.JSONObject; public class GaussianProcess { int optConstVar = 0; LikelihoodModel likelihoodModel = null; public void setLikelihoodModel(LikelihoodModel likelihoodModel) { this.likelihoodModel = likelihoodModel; } double sigmaN = 1; double[][] K; KernelFunction kernel; double[][] X; double[][] Y; double[][] dY; double[][] mX; double[][] Ky = null; double[][] invKy = null; double[][] invKyDY = null; double detKy; int matrices_are_computed = 0; double[] Kstar_temp; double[] Ks_invKy_temp; public GaussianProcess() { } public GaussianProcess(GaussianProcess gp) { this.optConstVar = gp.optConstVar; this.likelihoodModel = gp.likelihoodModel; this.kernel = gp.kernel; this.K = gp.K; this.sigmaN = gp.sigmaN; this.X = gp.X; this.Y = gp.Y; this.dY = gp.dY; this.mX = gp.mX; this.invKyDY = gp.invKyDY; this.matrices_are_computed = gp.matrices_are_computed; this.Kstar_temp = gp.Kstar_temp; this.Ks_invKy_temp = gp.Ks_invKy_temp; } public JSONObject toJSON() { JSONObject jobj = new JSONObject(); if (invKyDY == null) { return jobj; } JSONArray jsonInvKyDY = JSONUtils.toJSONArray(invKyDY); try { jobj.put("invKyDY", jsonInvKyDY); } catch (JSONException e) { e.printStackTrace(); } return jobj; } public void fromJSON(JSONObject jobj) { try { JSONArray jsonInvKyDY = jobj.getJSONArray("invKyDY"); invKyDY = JSONUtils.array2DFromJSONArray(jsonInvKyDY); } catch (JSONException e) { e.printStackTrace(); } } public void setKernel(KernelFunction kernel) { this.kernel = kernel; } public KernelFunction getKernel() { return this.kernel; } public void setSigmaN(double sigmaN) { this.sigmaN = sigmaN; } protected double[][] getX() { return X; } protected double[][] getY() { return Y; } protected double[][] getdY() { return dY; } protected double[][] computeGramMatrix(double[][] X) { int ns = X.length; double[][] K = new double[ns][ns]; for (int i = 0; i < ns; i++) { K[i][i] = kernel.computeKernel(X[i], X[i]); for (int j = i + 1; j < ns; j++) { double buff = kernel.computeKernel(X[i], X[j]); K[i][j] = buff; K[j][i] = buff; } } return K; } public GaussianProcess fit(double[][] X, double[][] Y) { int ns = X.length; this.X = X; this.Y = Y; // Compute Gram Matrix (Symmetric matrix) K = computeGramMatrix(X); RealMatrix KMat = MatrixUtils.createRealMatrix(K); RealMatrix sigma_n2I = MatrixUtils.createRealIdentityMatrix(ns).scalarMultiply(sigmaN * sigmaN); RealMatrix Ky = KMat.add(sigma_n2I); this.Ky = Ky.getData(); LUDecomposition LUDecomp = new LUDecomposition(Ky); detKy = LUDecomp.getDeterminant(); RealMatrix invKyMat = LUDecomp.getSolver().getInverse(); invKy = invKyMat.getData(); // Precompute dY = Y - mX int ny = Y[0].length; this.mX = new double[ns][ny]; this.dY = new double[ns][ny]; for (int i = 0; i < ns; i++) { for (int j = 0; j < ny; j++) { mX[i][j] = meanFunc(X[i], j); dY[i][j] = Y[i][j] - mX[i][j]; } } if (optConstVar == 1) { invKyDY = computeInvKyDY(invKy, dY); } return this; } double[][] getInvKyDY() { if (invKyDY == null) { invKyDY = computeInvKyDY(invKy, dY); } else { if (optConstVar == 1) { return invKyDY; } else { invKyDY = computeInvKyDY(invKy, dY); } } return this.invKyDY; } public void setOptConstVar(int optConstVar) { this.optConstVar = optConstVar; System.out.println("optConstVar=" + optConstVar + ": use a constant variance"); if (optConstVar == 1 && invKy != null && dY != null) { invKyDY = computeInvKyDY(invKy, dY); } } public int getOptConstVar() { return optConstVar; } double[][] computeInvKyDY(double[][] invKy, double[][] dY) { int ns = dY.length; int ny = dY[0].length; double[][] invKyDY = new double[ns][ny]; for (int i = 0; i < ns; i++) { for (int j = 0; j < ny; j++) { invKyDY[i][j] = 0; double tmp = 0; for (int k = 0; k < ns; k++) { tmp += invKy[i][k] * dY[k][j]; } invKyDY[i][j] = tmp; } } return invKyDY; } public double[][] predictYgivenX(double[][] X) { int nin = X.length; double[][] Y = getY(); int ny = Y[0].length; double[][] Ypred = new double[nin][ny]; for (int i = 0; i < nin; i++) { Ypred[i] = predictygivenx(X[i]); } return Ypred; } protected double meanFunc(double[] x, int index_y) { return meanFunc(x); } protected double meanFunc(double[] x) { return 0; } public double[] predictygivenx(double[] x) { return predictfgivenx(x); } public double[] predictfgivenx(double[] x) { if (optConstVar == 0) { return predictfgivenxActual(x); } else if (optConstVar == 1) { return predictfgivenxConstVar(x); } return null; } private void predict_y_var_givenx(double[] x, double[] fpred, double[] variance_y) { double[][] Y = getY(); double[][] dY = getdY(); int ns = Y.length; int ny = Y[0].length; double[] Kstar = computeKstar(x); double[] Ks_invKy = computeKs_invKy(Kstar, invKy); double kstar = kernel.computeKernel(x, x); double fi = 0; for (int i = 0; i < ny; i++) { fi = meanFunc(x, i); for (int k = 0; k < ns; k++) { fi += Ks_invKy[k] * dY[k][i]; } fpred[i] = fi; } double variance_f_i = kstar; for (int j = 0; j < ns; j++) { variance_f_i += -Ks_invKy[j] * Kstar[j]; } double sigmaN2 = sigmaN * sigmaN; for (int i = 0; i < ny; i++) { variance_y[i] = variance_f_i + sigmaN2; } } private double[] predictfgivenxActual(double[] x) { double[][] Y = getY(); double[][] dY = getdY(); int ns = Y.length; double[] Kstar; if (matrices_are_computed == 1) { Kstar = Kstar_temp; } else { Kstar = computeKstar(x); Kstar_temp = Kstar; } int ny = Y[0].length; double[] Ks_invKpSig; if (matrices_are_computed == 1) { Ks_invKpSig = Ks_invKy_temp; } else { Ks_invKpSig = computeKs_invKy(Kstar, invKy); Ks_invKy_temp = Ks_invKpSig; } double[] fpred = new double[ny]; double fi = 0; for (int i = 0; i < ny; i++) { fi = meanFunc(x, i); for (int k = 0; k < ns; k++) { fi += Ks_invKpSig[k] * dY[k][i]; } fpred[i] = fi; } return fpred; } private double[] predictfgivenxConstVar(double[] x) { double[][] Y = getY(); double[][] invKyDY = getInvKyDY(); int ns = Y.length; double[] Kstar = computeKstar(x); int ny = Y[0].length; double[] fpred = new double[ny]; double fi = 0; for (int i = 0; i < ny; i++) { fi = meanFunc(x, i); for (int k = 0; k < ns; k++) { fi += Kstar[k] * invKyDY[k][i]; } fpred[i] = fi; } return fpred; } private double[] computeKstar(double[] x) { int ns = X.length; double[] Kstar = new double[ns]; for (int i = 0; i < ns; i++) { Kstar[i] = kernel.computeKernel(x, X[i]); } return Kstar; } private double[] computeKs_invKy(double[] Kstar, double[][] invKy) { int ns = Kstar.length; double[] Ks_invKy = new double[ns]; for (int j = 0; j < ns; j++) { Ks_invKy[j] = 0.0; for (int k = 0; k < ns; k++) { double buff = Kstar[k] * invKy[k][j]; Ks_invKy[j] += buff; } } return Ks_invKy; } protected double[] variancefgivenx(double[] x) { if (optConstVar == 0) { return variancefgivenxActual(x); } else if (optConstVar == 1) { return variancefgivenxConstVar(x); } return null; } private double[] variancefgivenxActual(double[] x) { double[][] Y = getY(); int ns = Y.length; int ny = Y[0].length; double kstar = kernel.computeKernel(x, x); double[] Kstar; if (matrices_are_computed == 1) { Kstar = Kstar_temp; } else { Kstar = computeKstar(x); } double[] Ks_invKy; if (matrices_are_computed == 1) { Ks_invKy = Ks_invKy_temp; } else { Ks_invKy = computeKs_invKy(Kstar, invKy); } double variance_f_i = kstar; for (int j = 0; j < ns; j++) { variance_f_i += -Ks_invKy[j] * Kstar[j]; } double[] variance_f = new double[ny]; for (int i = 0; i < ny; i++) { variance_f[i] = variance_f_i; } Kstar_temp = Kstar; Ks_invKy_temp = Ks_invKy; return variance_f; } private double[] variancefgivenxConstVar(double[] x) { double[][] Y = getY(); int ny = Y[0].length; double kstar = kernel.computeKernel(x, x); double[] variance_f = new double[ny]; for (int i = 0; i < ny; i++) { variance_f[i] = kstar; } return variance_f; } public double[] varianceygivenx(double[] x) { double[] varfgiveny = variancefgivenx(x); int ny = varfgiveny.length; double varN = sigmaN * sigmaN; for (int i = 0; i < ny; i++) { varfgiveny[i] += varN; } return varfgiveny; } public double logProbaygivenx(double[] x, double[] y) { double[] mean_y = new double[y.length]; double[] var_y = new double[y.length]; if (optConstVar == 0) { predict_y_var_givenx(x, mean_y, var_y); } if (optConstVar == 1) { mean_y = predictfgivenx(x); var_y = varianceygivenx(x); } int ny = mean_y.length; double sumLogLL = 0; double sigma = 0; for (int j = 0; j < ny; j++) { sigma = Math.sqrt(var_y[j]); sumLogLL += logProba(y[j], mean_y[j], sigma); } return sumLogLL; } private double logProba(double y, double mean_y, double sigma) { if (likelihoodModel == null) { return StatUtils.logProbaNormal(y, mean_y, sigma); } else { return likelihoodModel.logProba(y, mean_y, sigma); } } public double[] logProbaygivenX(double[][] X, double[] y) { return logProbaygivenXMultiThread(X, y); } public double[] logProbaygivenXSingleThread(double[][] X, double[] y) { int n = X.length; final double logpro[] = new double[n]; for (int i = 0; i < n; i++) { logpro[i] = logProbaygivenx(X[i], y); } return logpro; } public double[] logProbaygivenXMultiThread(final double[][] X, final double[] y) { ExecutorService ex = ExecutorServiceHolder.getExecutorService(); int n = X.length; final double logpro[] = new double[n]; Future<?>[] futures = new Future<?>[n]; for (int i = 0; i < n; i++) { final int idx = i; futures[i] = ex.submit(new Runnable() { public void run() { logpro[idx] = logProbaygivenx(X[idx], y); } }); } for (int i = 0; i < n; i++) { try { futures[i].get(); } catch (InterruptedException | ExecutionException e) { e.printStackTrace(); } } return logpro; } public double looMSE() { double[][] Y = getY(); double[][] dY = getdY(); int ns = X.length; int ny = Y[0].length; RealMatrix invKy = MatrixUtils.createRealMatrix(this.invKy); RealMatrix K = MatrixUtils.createRealMatrix(this.K); RealMatrix Hmat = invKy.multiply(K); double[][] H = Hmat.getData(); double sum = 0; double count = 0; for (int j = 0; j < ny; j++) { for (int i = 0; i < ns; i++) { double preddY = 0; for (int k = 0; k < ns; k++) { preddY += H[i][k] * dY[k][j]; } double diff = (dY[i][j] - preddY) / (1.0 - H[i][i]); sum += (diff * diff); count += 1; } } sum /= count; return sum; } public double looPredLogLikelihood() { double[][] Y = getY(); double[][] dY = getdY(); int ns = X.length; int ny = Y[0].length; RealMatrix Ky = MatrixUtils.createRealMatrix(this.Ky); RealMatrix invKy = new LUDecomposition(Ky).getSolver().getInverse(); RealMatrix dYmat = MatrixUtils.createRealMatrix(dY); double[] LOOPredLL = new double[ny]; for (int j = 0; j < ny; j++) { RealMatrix dy = MatrixUtils.createColumnRealMatrix(dYmat.getColumn(j)); RealMatrix invKdy = invKy.multiply(dy); double sum = 0; for (int i = 0; i < ns; i++) { double mu_i = dYmat.getEntry(i, j) - invKdy.getEntry(i, 0) / invKy.getEntry(i, i); double sigma_i2 = 1 / invKy.getEntry(i, i); double logLL = StatUtils.logProbaNormal(dYmat.getEntry(i, j), mu_i, Math.sqrt(sigma_i2)); sum += logLL; } LOOPredLL[j] = sum; } double sumLOOPredLL = 0; for (int j = 0; j < ny; j++) { sumLOOPredLL += LOOPredLL[j]; } return sumLOOPredLL; } }