{- | Demonstration of Gaussian process classification using the  
     demonstration problem from

     www.gaussianprocess.org/gpml/code/matlab/doc/

     This demo uses the EP approximation approach.

     For details of the algorithms involved see www.gaussianprocesses.org. 
     For a detailed explanation of the following code see the HasGP user 
     manual.

     Copyright (C) 2011 Sean Holden. sbh11@cl.cam.ac.uk.
-}
{- This file is part of HasGP.

   HasGP is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   HasGP is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with HasGP.  If not, see <http://www.gnu.org/licenses/>.
-}
module HasGP.Demos.ClassificationDemo2 where

import Numeric.LinearAlgebra
import Numeric.GSL.Minimization

import HasGP.Types.MainTypes
import HasGP.Support.Linear
import HasGP.Classification.EP.ClassificationEP
import HasGP.Covariance.SquaredExpARD
import HasGP.Covariance.Basic

-- | This function defines when iteration stops.
stopEP::EPConvergenceTest
stopEP s1 s2 = ((count s2) == 100) || 
             ((eValue s1) > (eValue s2)) ||
             (abs ((eValue s1) - (eValue s2)) < 0.001)

demo = do 
  putStrLn "Loading the training data..."
  
  inputs <- loadMatrix "gpml-classifier-x.txt"
  targets <- fscanfVector "gpml-classifier-y.txt" 120
  points <- loadMatrix "gpml-classifier-test.txt"

  putStrLn "Learning and predicting: EP + hyperparameter optimization..."

  let cov = SquaredExponentialARD (log 1.0) (constant (log 1.0) 2)
  let c = covarianceMatrix cov inputs
  
  let f = (\v -> gpClassifierEPLogEvidenceVec inputs targets cov  
                 generateRandomSiteOrder stopEP v)
  let ev = fst . f
  let gev = snd . f
  let (solution, path) = 
       minimizeVD ConjugatePR 0.0001 50 1 0.0001 ev gev 
                      (constant (log 1) 3)
                      
  putStrLn $ "Solution: " ++ (show $ mapVector exp solution)
  putStrLn $ "Path: "
  putStrLn $ show path

  let cov' = SquaredExponentialARD (solution @> 0) 
                             (fromList [(solution @> 1), (solution @> 2)])
  let c' = covarianceMatrix cov' inputs
  let (epValue, epState) = 
         gpClassifierEPLearn c' targets generateRandomSiteOrder stopEP
  let classify = 
          gpClassifierEPPredict (siteState epValue) inputs targets c' cov'  
  let newOuts = classify points
  
  fprintfVector "gpml-hasgp-outputs.txt" "%g" newOuts

  putStrLn $ "Done"
  

  return ()