{- | HasGP Gaussian Process Library. This module contains assorted functions that support GP calculations and are specifically related to linear algebra. 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.Support.Linear where import Data.Packed.ST import Control.Monad (mapM_,zipWithM_) import Control.Monad.ST import Numeric.LinearAlgebra import HasGP.Types.MainTypes import HasGP.Support.Functions as F -- | Sum the elements in a vector. sumVector :: DVector -> Double sumVector = foldVector (+) 0.0 -- | Sum of elements in a vector, divided by an Int. sumVectorDiv :: Int -> DVector -> Double sumVectorDiv d v = (sumVector v)/(fromIntegral d) -- | Length of a vector. lengthV :: (Normed a b) => a b -> RealOf b lengthV = pnorm PNorm2 -- | Generate a vector equal to the first column of a matrix. toVector :: Matrix Double -> Vector Double toVector x = head $ toColumns x -- | Replace the element at a specified position in a vector. -- NOTE: hmatrix numbers from 0, which is odd. This numbers from 1. -- The result is returned by overwriting v. This is implemented -- via runSTVector because the increase in efficiency is HUGE. replaceInVector :: DVector -> Int -> Double -> DVector replaceInVector v i n | (1 <= i) && (i <= (dim v)) = runSTVector $ do v2 <- thawVector v writeVector v2 (i-1) n return v2 | otherwise = error "Index out of range in replaceInVector" -- | Efficiently pre multiply by a diagonal matrix (passed as a vector) preMultiply :: DVector -> DMatrix -> DMatrix preMultiply v m = fromRows $ zipWith scale (toList v) (toRows m) -- | Efficiently post multiply by a diagonal matrix (passed as a vector) postMultiply :: DMatrix -> DVector -> DMatrix postMultiply m v = fromColumns $ zipWith scale (toList v) (toColumns m) -- | Compute x^T A x when A is diagonal. The second argument is the -- diagonal of A. xAxDiag :: DVector -> DVector -> Double xAxDiag x a | (d == dim a) = a <.> (x * x) | otherwise = error "Incorrect dimensions in xAxDiag" where d = dim x -- | Compute the diagonal only of the product of two square matrices abDiagOnly :: DMatrix -> DMatrix -> DVector abDiagOnly a b = fromList $ zipWith (<.>) (toRows a) (toColumns b) -- | Compute ABA where A is diagonal. The first argument is the diagonal of A. abaDiagDiag :: DVector -> DMatrix -> DMatrix abaDiagDiag a b = (d><d) (zipWith (*) bL aA) where d = dim a aL = toList a aA = [(a1 * a2) | a1 <- aL, a2 <- aL] bL = toList $ join $ toRows b -- | Compute aBa where a is a vector and B is a matrix abaVV :: DVector -> DMatrix -> Double abaVV a b = (flatten ((asRow a) <> b)) <.> a