{- | Demonstration of Gaussian process classification using the  
     1-dimensional problem from Rasmussen and Williams' book.

     This demo compares the Laplace and EP approximation approaches.

     For details of the algorithms involved see www.gaussianprocesses.org.

     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.ClassificationDemo1 where 

import Numeric.LinearAlgebra
import Numeric.GSL.Minimization

import HasGP.Data.RWData1
import HasGP.Types.MainTypes
import HasGP.Support.Linear
import HasGP.Classification.Laplace.ClassificationLaplace as L
import HasGP.Classification.EP.ClassificationEP as EP
import HasGP.Covariance.SquaredExp
import HasGP.Covariance.Basic
import HasGP.Likelihood.LogPhi
import HasGP.Likelihood.Basic 

-- | This function defines when iteration stops for the Laplace version.
stopLaplace::LaplaceConvergenceTest
stopLaplace s1 s2 = ((L.count s2) == 100) || 
             (lengthV ((fValue s1) - (fValue s2)) < 0.001)

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

demo = do 
  putStrLn "Generating the training and testing data..."

  let (inputs, t) = simpleClassificationData 0
  let targets = mapVector (\x -> if (x==1) then 1 else -1) t
  saveMatrix "training-inputs.txt" "%g" inputs
  fprintfVector "training-targets.txt" "%g" targets

  let points = asColumn $ linspace 50 (-9.0,5.0) 
  saveMatrix "test-inputs.txt" "%g" points

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

  let cov1 = SquaredExponential (log 1.0) (log 1.0)
  let c1 = covarianceMatrix cov1 inputs
  
  let f1 = (\v -> gpCLaplaceLogEvidenceVec inputs targets cov1 LogPhi 
                 stopLaplace v)
  let ev1 = fst . f1
  let gev1 = snd . f1
  let (solution1, path1) = 
       minimizeVD ConjugatePR 0.0001 50 0.01 0.1 ev1 gev1 
                      (constant (log 1) 2)
                      
  putStrLn $ "Solution: " ++ (show $ mapVector exp solution1)
  putStrLn $ "Path: "
  putStrLn $ show path1

  let cov1' = SquaredExponential (solution1 @> 0) (solution1 @> 1)
  let c1' = covarianceMatrix cov1' inputs
  let result1 = gpCLaplaceLearn c1' targets LogPhi stopLaplace
  let classify1 = 
          gpCLaplacePredict' (fValue result1) inputs targets c1' cov1' LogPhi 
  let newOuts1 = fromList $ map convertToP_CG $ classify1 $ points
  
  fprintfVector "test-outputs-laplace.txt" "%g" newOuts1

  putStrLn $ "Done"

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

  let cov2 = SquaredExponential (log 1.0) (log 1.0)
  let c2 = covarianceMatrix cov2 inputs
  
  let f2 = (\v -> gpClassifierEPLogEvidenceVec inputs targets cov2  
                 generateRandomSiteOrder stopEP v)
  let ev2 = fst . f2
  let gev2 = snd . f2
  let (solution2, path2) = 
       minimizeVD ConjugatePR 0.0001 50 0.01 0.1 ev2 gev2 
                      (constant (log 1) 2)
                      
  putStrLn $ "Solution: " ++ (show $ mapVector exp solution2)
  putStrLn $ "Path: "
  putStrLn $ show path2

  let cov2' = SquaredExponential (solution2 @> 0) (solution2 @> 1)
  let c2' = covarianceMatrix cov2' inputs
  let (epValue, epState) = 
         gpClassifierEPLearn c2' targets generateRandomSiteOrder stopEP
  let classify2 = 
          gpClassifierEPPredict (siteState epValue) inputs targets c2' cov2'  
  let newOuts2 = classify2 points
  
  fprintfVector "test-outputs-ep.txt" "%g" newOuts2

  putStrLn $ "Done"

  return ()