Free cookie consent management tool by TermsFeed Policy Generator

Ignore:
Timestamp:
07/08/16 14:40:02 (8 years ago)
Author:
gkronber
Message:

#2434: merged trunk changes r12934:14026 from trunk to branch

Location:
branches/crossvalidation-2434
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • branches/crossvalidation-2434

  • branches/crossvalidation-2434/HeuristicLab.Algorithms.DataAnalysis

  • branches/crossvalidation-2434/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessModel.cs

    r12820 r14029  
    3434  [StorableClass]
    3535  [Item("GaussianProcessModel", "Represents a Gaussian process posterior.")]
    36   public sealed class GaussianProcessModel : NamedItem, IGaussianProcessModel {
     36  public sealed class GaussianProcessModel : RegressionModel, IGaussianProcessModel {
     37    public override IEnumerable<string> VariablesUsedForPrediction {
     38      get { return allowedInputVariables; }
     39    }
     40
    3741    [Storable]
    3842    private double negativeLogLikelihood;
     
    6165      get { return meanFunction; }
    6266    }
    63     [Storable]
    64     private string targetVariable;
    65     public string TargetVariable {
    66       get { return targetVariable; }
    67     }
     67
    6868    [Storable]
    6969    private string[] allowedInputVariables;
     
    124124      this.meanFunction = cloner.Clone(original.meanFunction);
    125125      this.covarianceFunction = cloner.Clone(original.covarianceFunction);
    126       this.inputScaling = cloner.Clone(original.inputScaling);
     126      if (original.inputScaling != null)
     127        this.inputScaling = cloner.Clone(original.inputScaling);
    127128      this.trainingDataset = cloner.Clone(original.trainingDataset);
    128129      this.negativeLogLikelihood = original.negativeLogLikelihood;
    129       this.targetVariable = original.targetVariable;
    130130      this.sqrSigmaNoise = original.sqrSigmaNoise;
    131131      if (original.meanParameter != null) {
     
    144144    }
    145145    public GaussianProcessModel(IDataset ds, string targetVariable, IEnumerable<string> allowedInputVariables, IEnumerable<int> rows,
    146       IEnumerable<double> hyp, IMeanFunction meanFunction, ICovarianceFunction covarianceFunction)
    147       : base() {
     146      IEnumerable<double> hyp, IMeanFunction meanFunction, ICovarianceFunction covarianceFunction,
     147      bool scaleInputs = true)
     148      : base(targetVariable) {
    148149      this.name = ItemName;
    149150      this.description = ItemDescription;
    150151      this.meanFunction = (IMeanFunction)meanFunction.Clone();
    151152      this.covarianceFunction = (ICovarianceFunction)covarianceFunction.Clone();
    152       this.targetVariable = targetVariable;
    153153      this.allowedInputVariables = allowedInputVariables.ToArray();
    154154
     
    163163                                             .ToArray();
    164164      sqrSigmaNoise = Math.Exp(2.0 * hyp.Last());
    165       CalculateModel(ds, rows);
    166     }
    167 
    168     private void CalculateModel(IDataset ds, IEnumerable<int> rows) {
     165      try {
     166        CalculateModel(ds, rows, scaleInputs);
     167      }
     168      catch (alglib.alglibexception ae) {
     169        // wrap exception so that calling code doesn't have to know about alglib implementation
     170        throw new ArgumentException("There was a problem in the calculation of the Gaussian process model", ae);
     171      }
     172    }
     173
     174    private void CalculateModel(IDataset ds, IEnumerable<int> rows, bool scaleInputs = true) {
    169175      this.trainingDataset = (IDataset)ds.Clone();
    170176      this.trainingRows = rows.ToArray();
    171       this.inputScaling = new Scaling(trainingDataset, allowedInputVariables, rows);
    172       this.x = CalculateX(trainingDataset, allowedInputVariables, rows, inputScaling);
    173       var y = ds.GetDoubleValues(targetVariable, rows);
     177      this.inputScaling = scaleInputs ? new Scaling(ds, allowedInputVariables, rows) : null;
     178
     179      x = GetData(ds, this.allowedInputVariables, this.trainingRows, this.inputScaling);
     180
     181      IEnumerable<double> y;
     182      y = ds.GetDoubleValues(TargetVariable, rows);
    174183
    175184      int n = x.GetLength(0);
    176185
     186      var columns = Enumerable.Range(0, x.GetLength(1)).ToArray();
    177187      // calculate cholesky decomposed (lower triangular) covariance matrix
    178       var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, x.GetLength(1)));
     188      var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, columns);
    179189      this.l = CalculateL(x, cov, sqrSigmaNoise);
    180190
    181191      // calculate mean
    182       var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, Enumerable.Range(0, x.GetLength(1)));
     192      var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, columns);
    183193      double[] m = Enumerable.Range(0, x.GetLength(0))
    184194        .Select(r => mean.Mean(x, r))
    185195        .ToArray();
    186 
    187 
    188196
    189197      // calculate sum of diagonal elements for likelihood
     
    219227      double[] meanGradients = new double[meanFunction.GetNumberOfParameters(nAllowedVariables)];
    220228      for (int k = 0; k < meanGradients.Length; k++) {
    221         var meanGrad = Enumerable.Range(0, alpha.Length)
    222         .Select(r => mean.Gradient(x, r, k));
     229        var meanGrad = new double[alpha.Length];
     230        for (int g = 0; g < meanGrad.Length; g++)
     231          meanGrad[g] = mean.Gradient(x, g, k);
    223232        meanGradients[k] = -Util.ScalarProd(meanGrad, alpha);
    224233      }
     
    228237        for (int i = 0; i < n; i++) {
    229238          for (int j = 0; j < i; j++) {
    230             var g = cov.CovarianceGradient(x, i, j).ToArray();
     239            var g = cov.CovarianceGradient(x, i, j);
    231240            for (int k = 0; k < covGradients.Length; k++) {
    232241              covGradients[k] += lCopy[i, j] * g[k];
     
    234243          }
    235244
    236           var gDiag = cov.CovarianceGradient(x, i, i).ToArray();
     245          var gDiag = cov.CovarianceGradient(x, i, i);
    237246          for (int k = 0; k < covGradients.Length; k++) {
    238247            // diag
     
    249258    }
    250259
    251     private static double[,] CalculateX(IDataset ds, IEnumerable<string> allowedInputVariables, IEnumerable<int> rows, Scaling inputScaling) {
    252       return AlglibUtil.PrepareAndScaleInputMatrix(ds, allowedInputVariables, rows, inputScaling);
     260    private static double[,] GetData(IDataset ds, IEnumerable<string> allowedInputs, IEnumerable<int> rows, Scaling scaling) {
     261      if (scaling != null) {
     262        return AlglibUtil.PrepareAndScaleInputMatrix(ds, allowedInputs, rows, scaling);
     263      } else {
     264        return AlglibUtil.PrepareInputMatrix(ds, allowedInputs, rows);
     265      }
    253266    }
    254267
     
    286299
    287300    #region IRegressionModel Members
    288     public IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {
     301    public override IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {
    289302      return GetEstimatedValuesHelper(dataset, rows);
    290303    }
    291     public GaussianProcessRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
     304    public override IRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
    292305      return new GaussianProcessRegressionSolution(this, new RegressionProblemData(problemData));
    293306    }
    294     IRegressionSolution IRegressionModel.CreateRegressionSolution(IRegressionProblemData problemData) {
    295       return CreateRegressionSolution(problemData);
    296     }
    297307    #endregion
    298308
    299309
    300310    private IEnumerable<double> GetEstimatedValuesHelper(IDataset dataset, IEnumerable<int> rows) {
    301       if (x == null) {
    302         this.x = CalculateX(trainingDataset, allowedInputVariables, trainingRows, inputScaling);
    303       }
    304       int n = x.GetLength(0);
    305 
    306       var newX = AlglibUtil.PrepareAndScaleInputMatrix(dataset, allowedInputVariables, rows, inputScaling);
    307       int newN = newX.GetLength(0);
    308 
    309       var Ks = new double[newN, n];
    310       var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, Enumerable.Range(0, newX.GetLength(1)));
    311       var ms = Enumerable.Range(0, newX.GetLength(0))
    312       .Select(r => mean.Mean(newX, r))
    313       .ToArray();
    314       var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, newX.GetLength(1)));
    315       for (int i = 0; i < newN; i++) {
    316         for (int j = 0; j < n; j++) {
    317           Ks[i, j] = cov.CrossCovariance(x, newX, j, i);
    318         }
    319       }
    320 
    321       return Enumerable.Range(0, newN)
    322         .Select(i => ms[i] + Util.ScalarProd(Util.GetRow(Ks, i), alpha));
     311      try {
     312        if (x == null) {
     313          x = GetData(trainingDataset, allowedInputVariables, trainingRows, inputScaling);
     314        }
     315        int n = x.GetLength(0);
     316
     317        double[,] newX = GetData(dataset, allowedInputVariables, rows, inputScaling);
     318        int newN = newX.GetLength(0);
     319
     320        var Ks = new double[newN][];
     321        var columns = Enumerable.Range(0, newX.GetLength(1)).ToArray();
     322        var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, columns);
     323        var ms = Enumerable.Range(0, newX.GetLength(0))
     324        .Select(r => mean.Mean(newX, r))
     325        .ToArray();
     326        var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, columns);
     327        for (int i = 0; i < newN; i++) {
     328          Ks[i] = new double[n];
     329          for (int j = 0; j < n; j++) {
     330            Ks[i][j] = cov.CrossCovariance(x, newX, j, i);
     331          }
     332        }
     333
     334        return Enumerable.Range(0, newN)
     335          .Select(i => ms[i] + Util.ScalarProd(Ks[i], alpha));
     336      }
     337      catch (alglib.alglibexception ae) {
     338        // wrap exception so that calling code doesn't have to know about alglib implementation
     339        throw new ArgumentException("There was a problem in the calculation of the Gaussian process model", ae);
     340      }
    323341    }
    324342
    325343    public IEnumerable<double> GetEstimatedVariance(IDataset dataset, IEnumerable<int> rows) {
    326       if (x == null) {
    327         this.x = CalculateX(trainingDataset, allowedInputVariables, trainingRows, inputScaling);
    328       }
    329       int n = x.GetLength(0);
    330 
    331       var newX = AlglibUtil.PrepareAndScaleInputMatrix(dataset, allowedInputVariables, rows, inputScaling);
    332       int newN = newX.GetLength(0);
    333 
    334       var kss = new double[newN];
    335       double[,] sWKs = new double[n, newN];
    336       var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, x.GetLength(1)));
    337 
    338       if (l == null) {
    339         l = CalculateL(x, cov, sqrSigmaNoise);
    340       }
    341 
    342       // for stddev
    343       for (int i = 0; i < newN; i++)
    344         kss[i] = cov.Covariance(newX, i, i);
    345 
    346       for (int i = 0; i < newN; i++) {
    347         for (int j = 0; j < n; j++) {
    348           sWKs[j, i] = cov.CrossCovariance(x, newX, j, i) / Math.Sqrt(sqrSigmaNoise);
    349         }
    350       }
    351 
    352       // for stddev
    353       alglib.ablas.rmatrixlefttrsm(n, newN, l, 0, 0, false, false, 0, ref sWKs, 0, 0);
    354 
    355       for (int i = 0; i < newN; i++) {
    356         var sumV = Util.ScalarProd(Util.GetCol(sWKs, i), Util.GetCol(sWKs, i));
    357         kss[i] -= sumV;
    358         if (kss[i] < 0) kss[i] = 0;
    359       }
    360       return kss;
    361     }
     344      try {
     345        if (x == null) {
     346          x = GetData(trainingDataset, allowedInputVariables, trainingRows, inputScaling);
     347        }
     348        int n = x.GetLength(0);
     349
     350        var newX = GetData(dataset, allowedInputVariables, rows, inputScaling);
     351        int newN = newX.GetLength(0);
     352
     353        var kss = new double[newN];
     354        double[,] sWKs = new double[n, newN];
     355        var columns = Enumerable.Range(0, newX.GetLength(1)).ToArray();
     356        var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, columns);
     357
     358        if (l == null) {
     359          l = CalculateL(x, cov, sqrSigmaNoise);
     360        }
     361
     362        // for stddev
     363        for (int i = 0; i < newN; i++)
     364          kss[i] = cov.Covariance(newX, i, i);
     365
     366        for (int i = 0; i < newN; i++) {
     367          for (int j = 0; j < n; j++) {
     368            sWKs[j, i] = cov.CrossCovariance(x, newX, j, i) / Math.Sqrt(sqrSigmaNoise);
     369          }
     370        }
     371
     372        // for stddev
     373        alglib.ablas.rmatrixlefttrsm(n, newN, l, 0, 0, false, false, 0, ref sWKs, 0, 0);
     374
     375        for (int i = 0; i < newN; i++) {
     376          var col = Util.GetCol(sWKs, i).ToArray();
     377          var sumV = Util.ScalarProd(col, col);
     378          kss[i] += sqrSigmaNoise; // kss is V(f), add noise variance of predictive distibution to get V(y)
     379          kss[i] -= sumV;
     380          if (kss[i] < 0) kss[i] = 0;
     381        }
     382        return kss;
     383      }
     384      catch (alglib.alglibexception ae) {
     385        // wrap exception so that calling code doesn't have to know about alglib implementation
     386        throw new ArgumentException("There was a problem in the calculation of the Gaussian process model", ae);
     387      }
     388    }
     389
    362390  }
    363391}
Note: See TracChangeset for help on using the changeset viewer.