From 1925c123d7d8184a1d2ddc0a413e0fd2776e1083 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Thu, 8 May 2014 08:48:12 +0200 Subject: empty hmatrix-base --- packages/hmatrix/src/Numeric/Chain.hs | 140 ++ packages/hmatrix/src/Numeric/Container.hs | 303 ++++ packages/hmatrix/src/Numeric/ContainerBoot.hs | 611 ++++++++ packages/hmatrix/src/Numeric/Conversion.hs | 91 ++ packages/hmatrix/src/Numeric/GSL.hs | 43 + .../hmatrix/src/Numeric/GSL/Differentiation.hs | 87 ++ packages/hmatrix/src/Numeric/GSL/Fitting.hs | 179 +++ packages/hmatrix/src/Numeric/GSL/Fourier.hs | 47 + packages/hmatrix/src/Numeric/GSL/Integration.hs | 250 ++++ packages/hmatrix/src/Numeric/GSL/Internal.hs | 76 + packages/hmatrix/src/Numeric/GSL/Minimization.hs | 227 +++ packages/hmatrix/src/Numeric/GSL/ODE.hs | 138 ++ packages/hmatrix/src/Numeric/GSL/Polynomials.hs | 58 + packages/hmatrix/src/Numeric/GSL/Root.hs | 199 +++ packages/hmatrix/src/Numeric/GSL/Vector.hs | 328 +++++ packages/hmatrix/src/Numeric/GSL/gsl-aux.c | 1541 ++++++++++++++++++++ packages/hmatrix/src/Numeric/GSL/gsl-ode.c | 182 +++ packages/hmatrix/src/Numeric/HMatrix.hs | 136 ++ packages/hmatrix/src/Numeric/HMatrix/Data.hs | 69 + packages/hmatrix/src/Numeric/HMatrix/Devel.hs | 69 + packages/hmatrix/src/Numeric/IO.hs | 165 +++ packages/hmatrix/src/Numeric/LinearAlgebra.hs | 30 + .../src/Numeric/LinearAlgebra/Algorithms.hs | 746 ++++++++++ .../hmatrix/src/Numeric/LinearAlgebra/LAPACK.hs | 555 +++++++ .../src/Numeric/LinearAlgebra/LAPACK/lapack-aux.c | 1489 +++++++++++++++++++ .../src/Numeric/LinearAlgebra/LAPACK/lapack-aux.h | 60 + packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs | 295 ++++ .../src/Numeric/LinearAlgebra/Util/Convolution.hs | 114 ++ packages/hmatrix/src/Numeric/Matrix.hs | 98 ++ packages/hmatrix/src/Numeric/Vector.hs | 158 ++ 30 files changed, 8484 insertions(+) create mode 100644 packages/hmatrix/src/Numeric/Chain.hs create mode 100644 packages/hmatrix/src/Numeric/Container.hs create mode 100644 packages/hmatrix/src/Numeric/ContainerBoot.hs create mode 100644 packages/hmatrix/src/Numeric/Conversion.hs create mode 100644 packages/hmatrix/src/Numeric/GSL.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/Differentiation.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/Fitting.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/Fourier.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/Integration.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/Internal.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/Minimization.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/ODE.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/Polynomials.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/Root.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/Vector.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/gsl-aux.c create mode 100644 packages/hmatrix/src/Numeric/GSL/gsl-ode.c create mode 100644 packages/hmatrix/src/Numeric/HMatrix.hs create mode 100644 packages/hmatrix/src/Numeric/HMatrix/Data.hs create mode 100644 packages/hmatrix/src/Numeric/HMatrix/Devel.hs create mode 100644 packages/hmatrix/src/Numeric/IO.hs create mode 100644 packages/hmatrix/src/Numeric/LinearAlgebra.hs create mode 100644 packages/hmatrix/src/Numeric/LinearAlgebra/Algorithms.hs create mode 100644 packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK.hs create mode 100644 packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.c create mode 100644 packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.h create mode 100644 packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs create mode 100644 packages/hmatrix/src/Numeric/LinearAlgebra/Util/Convolution.hs create mode 100644 packages/hmatrix/src/Numeric/Matrix.hs create mode 100644 packages/hmatrix/src/Numeric/Vector.hs (limited to 'packages/hmatrix/src/Numeric') diff --git a/packages/hmatrix/src/Numeric/Chain.hs b/packages/hmatrix/src/Numeric/Chain.hs new file mode 100644 index 0000000..e1ab7da --- /dev/null +++ b/packages/hmatrix/src/Numeric/Chain.hs @@ -0,0 +1,140 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Chain +-- Copyright : (c) Vivian McPhail 2010 +-- License : GPL-style +-- +-- Maintainer : Vivian McPhail gmail.com> +-- Stability : provisional +-- Portability : portable +-- +-- optimisation of association order for chains of matrix multiplication +-- +----------------------------------------------------------------------------- + +module Numeric.Chain ( + optimiseMult, + ) where + +import Data.Maybe + +import Data.Packed.Matrix +import Numeric.ContainerBoot + +import qualified Data.Array.IArray as A + +----------------------------------------------------------------------------- +{- | + Provide optimal association order for a chain of matrix multiplications + and apply the multiplications. + + The algorithm is the well-known O(n\^3) dynamic programming algorithm + that builds a pyramid of optimal associations. + +> m1, m2, m3, m4 :: Matrix Double +> m1 = (10><15) [1..] +> m2 = (15><20) [1..] +> m3 = (20><5) [1..] +> m4 = (5><10) [1..] + +> >>> optimiseMult [m1,m2,m3,m4] + +will perform @((m1 `multiply` (m2 `multiply` m3)) `multiply` m4)@ + +The naive left-to-right multiplication would take @4500@ scalar multiplications +whereas the optimised version performs @2750@ scalar multiplications. The complexity +in this case is 32 (= 4^3/2) * (2 comparisons, 3 scalar multiplications, 3 scalar additions, +5 lookups, 2 updates) + a constant (= three table allocations) +-} +optimiseMult :: Product t => [Matrix t] -> Matrix t +optimiseMult = chain + +----------------------------------------------------------------------------- + +type Matrices a = A.Array Int (Matrix a) +type Sizes = A.Array Int (Int,Int) +type Cost = A.Array Int (A.Array Int (Maybe Int)) +type Indexes = A.Array Int (A.Array Int (Maybe ((Int,Int),(Int,Int)))) + +update :: A.Array Int (A.Array Int a) -> (Int,Int) -> a -> A.Array Int (A.Array Int a) +update a (r,c) e = a A.// [(r,(a A.! r) A.// [(c,e)])] + +newWorkSpaceCost :: Int -> A.Array Int (A.Array Int (Maybe Int)) +newWorkSpaceCost n = A.array (1,n) $ map (\i -> (i, subArray i)) [1..n] + where subArray i = A.listArray (1,i) (repeat Nothing) + +newWorkSpaceIndexes :: Int -> A.Array Int (A.Array Int (Maybe ((Int,Int),(Int,Int)))) +newWorkSpaceIndexes n = A.array (1,n) $ map (\i -> (i, subArray i)) [1..n] + where subArray i = A.listArray (1,i) (repeat Nothing) + +matricesToSizes :: [Matrix a] -> Sizes +matricesToSizes ms = A.listArray (1,length ms) $ map (\m -> (rows m,cols m)) ms + +chain :: Product a => [Matrix a] -> Matrix a +chain [] = error "chain: zero matrices to multiply" +chain [m] = m +chain [ml,mr] = ml `multiply` mr +chain ms = let ln = length ms + ma = A.listArray (1,ln) ms + mz = matricesToSizes ms + i = chain_cost mz + in chain_paren (ln,ln) i ma + +chain_cost :: Sizes -> Indexes +chain_cost mz = let (_,u) = A.bounds mz + cost = newWorkSpaceCost u + ixes = newWorkSpaceIndexes u + (_,_,i) = foldl chain_cost' (mz,cost,ixes) (order u) + in i + +chain_cost' :: (Sizes,Cost,Indexes) -> (Int,Int) -> (Sizes,Cost,Indexes) +chain_cost' sci@(mz,cost,ixes) (r,c) + | c == 1 = let cost' = update cost (r,c) (Just 0) + ixes' = update ixes (r,c) (Just ((r,c),(r,c))) + in (mz,cost',ixes') + | otherwise = minimum_cost sci (r,c) + +minimum_cost :: (Sizes,Cost,Indexes) -> (Int,Int) -> (Sizes,Cost,Indexes) +minimum_cost sci fu = foldl (smaller_cost fu) sci (fulcrum_order fu) + +smaller_cost :: (Int,Int) -> (Sizes,Cost,Indexes) -> ((Int,Int),(Int,Int)) -> (Sizes,Cost,Indexes) +smaller_cost (r,c) (mz,cost,ixes) ix@((lr,lc),(rr,rc)) = let op_cost = fromJust ((cost A.! lr) A.! lc) + + fromJust ((cost A.! rr) A.! rc) + + fst (mz A.! (lr-lc+1)) + * snd (mz A.! lc) + * snd (mz A.! rr) + cost' = (cost A.! r) A.! c + in case cost' of + Nothing -> let cost'' = update cost (r,c) (Just op_cost) + ixes'' = update ixes (r,c) (Just ix) + in (mz,cost'',ixes'') + Just ct -> if op_cost < ct then + let cost'' = update cost (r,c) (Just op_cost) + ixes'' = update ixes (r,c) (Just ix) + in (mz,cost'',ixes'') + else (mz,cost,ixes) + + +fulcrum_order (r,c) = let fs' = zip (repeat r) [1..(c-1)] + in map (partner (r,c)) fs' + +partner (r,c) (a,b) = ((r-b, c-b), (a,b)) + +order 0 = [] +order n = order (n-1) ++ zip (repeat n) [1..n] + +chain_paren :: Product a => (Int,Int) -> Indexes -> Matrices a -> Matrix a +chain_paren (r,c) ixes ma = let ((lr,lc),(rr,rc)) = fromJust $ (ixes A.! r) A.! c + in if lr == rr && lc == rc then (ma A.! lr) + else (chain_paren (lr,lc) ixes ma) `multiply` (chain_paren (rr,rc) ixes ma) + +-------------------------------------------------------------------------- + +{- TESTS -} + +-- optimal association is ((m1*(m2*m3))*m4) +m1, m2, m3, m4 :: Matrix Double +m1 = (10><15) [1..] +m2 = (15><20) [1..] +m3 = (20><5) [1..] +m4 = (5><10) [1..] \ No newline at end of file diff --git a/packages/hmatrix/src/Numeric/Container.hs b/packages/hmatrix/src/Numeric/Container.hs new file mode 100644 index 0000000..7e46147 --- /dev/null +++ b/packages/hmatrix/src/Numeric/Container.hs @@ -0,0 +1,303 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE UndecidableInstances #-} + +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Container +-- Copyright : (c) Alberto Ruiz 2010-14 +-- License : GPL +-- +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- Portability : portable +-- +-- Basic numeric operations on 'Vector' and 'Matrix', including conversion routines. +-- +-- The 'Container' class is used to define optimized generic functions which work +-- on 'Vector' and 'Matrix' with real or complex elements. +-- +-- Some of these functions are also available in the instances of the standard +-- numeric Haskell classes provided by "Numeric.LinearAlgebra". +-- +----------------------------------------------------------------------------- +{-# OPTIONS_HADDOCK hide #-} + +module Numeric.Container ( + -- * Basic functions + module Data.Packed, + konst, build, + constant, linspace, + diag, ident, + ctrans, + -- * Generic operations + Container(..), + -- * Matrix product + Product(..), udot, + Mul(..), + Contraction(..), + optimiseMult, + mXm,mXv,vXm,LSDiv(..), cdot, (·), dot, (<.>), + outer, kronecker, + -- * Random numbers + RandDist(..), + randomVector, + gaussianSample, + uniformSample, + meanCov, + -- * Element conversion + Convert(..), + Complexable(), + RealElement(), + + RealOf, ComplexOf, SingleOf, DoubleOf, + + IndexOf, + module Data.Complex, + -- * IO + dispf, disps, dispcf, vecdisp, latexFormat, format, + loadMatrix, saveMatrix, fromFile, fileDimensions, + readMatrix, + fscanfVector, fprintfVector, freadVector, fwriteVector, +) where + +import Data.Packed +import Data.Packed.Internal(constantD) +import Numeric.ContainerBoot +import Numeric.Chain +import Numeric.IO +import Data.Complex +import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD) +import Data.Packed.Random + +------------------------------------------------------------------ + +{- | creates a vector with a given number of equal components: + +@> constant 2 7 +7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@ +-} +constant :: Element a => a -> Int -> Vector a +-- constant x n = runSTVector (newVector x n) +constant = constantD-- about 2x faster + +{- | Creates a real vector containing a range of values: + +>>> linspace 5 (-3,7::Double) +fromList [-3.0,-0.5,2.0,4.5,7.0]@ + +>>> linspace 5 (8,2+i) :: Vector (Complex Double) +fromList [8.0 :+ 0.0,6.5 :+ 0.25,5.0 :+ 0.5,3.5 :+ 0.75,2.0 :+ 1.0] + +Logarithmic spacing can be defined as follows: + +@logspace n (a,b) = 10 ** linspace n (a,b)@ +-} +linspace :: (Container Vector e) => Int -> (e, e) -> Vector e +linspace 0 (a,b) = fromList[(a+b)/2] +linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n-1] + where s = (b-a)/fromIntegral (n-1) + +-- | dot product: @cdot u v = 'udot' ('conj' u) v@ +cdot :: (Container Vector t, Product t) => Vector t -> Vector t -> t +cdot u v = udot (conj u) v + +-------------------------------------------------------- + +class Contraction a b c | a b -> c, c -> a b + where + infixr 7 × + {- | Matrix-matrix product, matrix-vector product, and unconjugated dot product + +(unicode 0x00d7, multiplication sign) + +Examples: + +>>> let a = (3><4) [1..] :: Matrix Double +>>> let v = fromList [1,0,2,-1] :: Vector Double +>>> let u = fromList [1,2,3] :: Vector Double + +>>> a +(3><4) + [ 1.0, 2.0, 3.0, 4.0 + , 5.0, 6.0, 7.0, 8.0 + , 9.0, 10.0, 11.0, 12.0 ] + +matrix × matrix: + +>>> disp 2 (a × trans a) +3x3 + 30 70 110 + 70 174 278 +110 278 446 + +matrix × vector: + +>>> a × v +fromList [3.0,11.0,19.0] + +unconjugated dot product: + +>>> fromList [1,i] × fromList[2*i+1,3] +1.0 :+ 5.0 + +(×) is right associative, so we can write: + +>>> u × a × v +82.0 :: Double + +-} + (×) :: a -> b -> c + +instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where + (×) = mXv + +instance Product t => Contraction (Matrix t) (Matrix t) (Matrix t) where + (×) = mXm + +instance Contraction (Vector Double) (Vector Double) Double where + (×) = udot + +instance Contraction (Vector Float) (Vector Float) Float where + (×) = udot + +instance Contraction (Vector (Complex Double)) (Vector (Complex Double)) (Complex Double) where + (×) = udot + +instance Contraction (Vector (Complex Float)) (Vector (Complex Float)) (Complex Float) where + (×) = udot + + +-- | alternative function for the matrix product (×) +mmul :: Contraction a b c => a -> b -> c +mmul = (×) + +-------------------------------------------------------------------------------- + +class Mul a b c | a b -> c where + infixl 7 <> + -- | Matrix-matrix, matrix-vector, and vector-matrix products. + (<>) :: Product t => a t -> b t -> c t + +instance Mul Matrix Matrix Matrix where + (<>) = mXm + +instance Mul Matrix Vector Vector where + (<>) m v = flatten $ m <> asColumn v + +instance Mul Vector Matrix Vector where + (<>) v m = flatten $ asRow v <> m + +-------------------------------------------------------------------------------- + +class LSDiv c where + infixl 7 <\> + -- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD) + (<\>) :: Field t => Matrix t -> c t -> c t + +instance LSDiv Vector where + m <\> v = flatten (linearSolveSVD m (reshape 1 v)) + +instance LSDiv Matrix where + (<\>) = linearSolveSVD + +-------------------------------------------------------- + +{- | Dot product : @u · v = 'cdot' u v@ + + (unicode 0x00b7, middle dot, Alt-Gr .) + +>>> fromList [1,i] · fromList[2*i+1,3] +1.0 :+ (-1.0) + +-} +(·) :: (Container Vector t, Product t) => Vector t -> Vector t -> t +infixl 7 · +u · v = cdot u v + +-------------------------------------------------------------------------------- + +-- bidirectional type inference +class Konst e d c | d -> c, c -> d + where + -- | + -- >>> konst 7 3 :: Vector Float + -- fromList [7.0,7.0,7.0] + -- + -- >>> konst i (3::Int,4::Int) + -- (3><4) + -- [ 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 + -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 + -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 ] + -- + konst :: e -> d -> c e + +instance Container Vector e => Konst e Int Vector + where + konst = konst' + +instance Container Vector e => Konst e (Int,Int) Matrix + where + konst = konst' + +-------------------------------------------------------------------------------- + +class Build d f c e | d -> c, c -> d, f -> e, f -> d, f -> c, c e -> f, d e -> f + where + -- | + -- >>> build 5 (**2) :: Vector Double + -- fromList [0.0,1.0,4.0,9.0,16.0] + -- + -- Hilbert matrix of order N: + -- + -- >>> let hilb n = build (n,n) (\i j -> 1/(i+j+1)) :: Matrix Double + -- >>> putStr . dispf 2 $ hilb 3 + -- 3x3 + -- 1.00 0.50 0.33 + -- 0.50 0.33 0.25 + -- 0.33 0.25 0.20 + -- + build :: d -> f -> c e + +instance Container Vector e => Build Int (e -> e) Vector e + where + build = build' + +instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e + where + build = build' + +-------------------------------------------------------------------------------- + +{- | Compute mean vector and covariance matrix of the rows of a matrix. + +>>> meanCov $ gaussianSample 666 1000 (fromList[4,5]) (diagl[2,3]) +(fromList [4.010341078059521,5.0197204699640405], +(2><2) + [ 1.9862461923890056, -1.0127225830525157e-2 + , -1.0127225830525157e-2, 3.0373954915729318 ]) + +-} +meanCov :: Matrix Double -> (Vector Double, Matrix Double) +meanCov x = (med,cov) where + r = rows x + k = 1 / fromIntegral r + med = konst k r `vXm` x + meds = konst 1 r `outer` med + xc = x `sub` meds + cov = scale (recip (fromIntegral (r-1))) (trans xc `mXm` xc) + +-------------------------------------------------------------------------------- + +{-# DEPRECATED dot "use udot" #-} +dot :: Product e => Vector e -> Vector e -> e +dot = udot + +-- | contraction operator, equivalent to (x) +infixr 7 <.> +(<.>) :: Contraction a b c => a -> b -> c +(<.>) = (×) + diff --git a/packages/hmatrix/src/Numeric/ContainerBoot.hs b/packages/hmatrix/src/Numeric/ContainerBoot.hs new file mode 100644 index 0000000..ea4262c --- /dev/null +++ b/packages/hmatrix/src/Numeric/ContainerBoot.hs @@ -0,0 +1,611 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} + +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.ContainerBoot +-- Copyright : (c) Alberto Ruiz 2010 +-- License : GPL-style +-- +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- Portability : portable +-- +-- Module to avoid cyclyc dependencies. +-- +----------------------------------------------------------------------------- + +module Numeric.ContainerBoot ( + -- * Basic functions + ident, diag, ctrans, + -- * Generic operations + Container(..), + -- * Matrix product and related functions + Product(..), udot, + mXm,mXv,vXm, + outer, kronecker, + -- * Element conversion + Convert(..), + Complexable(), + RealElement(), + + RealOf, ComplexOf, SingleOf, DoubleOf, + + IndexOf, + module Data.Complex +) where + +import Data.Packed +import Data.Packed.ST as ST +import Numeric.Conversion +import Data.Packed.Internal +import Numeric.GSL.Vector +import Data.Complex +import Control.Applicative((<*>)) + +import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ) + +------------------------------------------------------------------- + +type family IndexOf (c :: * -> *) + +type instance IndexOf Vector = Int +type instance IndexOf Matrix = (Int,Int) + +type family ArgOf (c :: * -> *) a + +type instance ArgOf Vector a = a -> a +type instance ArgOf Matrix a = a -> a -> a + +------------------------------------------------------------------- + +-- | Basic element-by-element functions for numeric containers +class (Complexable c, Fractional e, Element e) => Container c e where + -- | create a structure with a single element + -- + -- >>> let v = fromList [1..3::Double] + -- >>> v / scalar (norm2 v) + -- fromList [0.2672612419124244,0.5345224838248488,0.8017837257372732] + -- + scalar :: e -> c e + -- | complex conjugate + conj :: c e -> c e + scale :: e -> c e -> c e + -- | scale the element by element reciprocal of the object: + -- + -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@ + scaleRecip :: e -> c e -> c e + addConstant :: e -> c e -> c e + add :: c e -> c e -> c e + sub :: c e -> c e -> c e + -- | element by element multiplication + mul :: c e -> c e -> c e + -- | element by element division + divide :: c e -> c e -> c e + equal :: c e -> c e -> Bool + -- + -- element by element inverse tangent + arctan2 :: c e -> c e -> c e + -- + -- | cannot implement instance Functor because of Element class constraint + cmap :: (Element b) => (e -> b) -> c e -> c b + -- | constant structure of given size + konst' :: e -> IndexOf c -> c e + -- | create a structure using a function + -- + -- Hilbert matrix of order N: + -- + -- @hilb n = build' (n,n) (\\i j -> 1/(i+j+1))@ + build' :: IndexOf c -> (ArgOf c e) -> c e + -- | indexing function + atIndex :: c e -> IndexOf c -> e + -- | index of min element + minIndex :: c e -> IndexOf c + -- | index of max element + maxIndex :: c e -> IndexOf c + -- | value of min element + minElement :: c e -> e + -- | value of max element + maxElement :: c e -> e + -- the C functions sumX/prodX are twice as fast as using foldVector + -- | the sum of elements (faster than using @fold@) + sumElements :: c e -> e + -- | the product of elements (faster than using @fold@) + prodElements :: c e -> e + + -- | A more efficient implementation of @cmap (\\x -> if x>0 then 1 else 0)@ + -- + -- >>> step $ linspace 5 (-1,1::Double) + -- 5 |> [0.0,0.0,0.0,1.0,1.0] + -- + + step :: RealElement e => c e -> c e + + -- | Element by element version of @case compare a b of {LT -> l; EQ -> e; GT -> g}@. + -- + -- Arguments with any dimension = 1 are automatically expanded: + -- + -- >>> cond ((1><4)[1..]) ((3><1)[1..]) 0 100 ((3><4)[1..]) :: Matrix Double + -- (3><4) + -- [ 100.0, 2.0, 3.0, 4.0 + -- , 0.0, 100.0, 7.0, 8.0 + -- , 0.0, 0.0, 100.0, 12.0 ] + -- + + cond :: RealElement e + => c e -- ^ a + -> c e -- ^ b + -> c e -- ^ l + -> c e -- ^ e + -> c e -- ^ g + -> c e -- ^ result + + -- | Find index of elements which satisfy a predicate + -- + -- >>> find (>0) (ident 3 :: Matrix Double) + -- [(0,0),(1,1),(2,2)] + -- + + find :: (e -> Bool) -> c e -> [IndexOf c] + + -- | Create a structure from an association list + -- + -- >>> assoc 5 0 [(3,7),(1,4)] :: Vector Double + -- fromList [0.0,4.0,0.0,7.0,0.0] + -- + -- >>> assoc (2,3) 0 [((0,2),7),((1,0),2*i-3)] :: Matrix (Complex Double) + -- (2><3) + -- [ 0.0 :+ 0.0, 0.0 :+ 0.0, 7.0 :+ 0.0 + -- , (-3.0) :+ 2.0, 0.0 :+ 0.0, 0.0 :+ 0.0 ] + -- + assoc :: IndexOf c -- ^ size + -> e -- ^ default value + -> [(IndexOf c, e)] -- ^ association list + -> c e -- ^ result + + -- | Modify a structure using an update function + -- + -- >>> accum (ident 5) (+) [((1,1),5),((0,3),3)] :: Matrix Double + -- (5><5) + -- [ 1.0, 0.0, 0.0, 3.0, 0.0 + -- , 0.0, 6.0, 0.0, 0.0, 0.0 + -- , 0.0, 0.0, 1.0, 0.0, 0.0 + -- , 0.0, 0.0, 0.0, 1.0, 0.0 + -- , 0.0, 0.0, 0.0, 0.0, 1.0 ] + -- + -- computation of histogram: + -- + -- >>> accum (konst 0 7) (+) (map (flip (,) 1) [4,5,4,1,5,2,5]) :: Vector Double + -- fromList [0.0,1.0,1.0,0.0,2.0,3.0,0.0] + -- + + accum :: c e -- ^ initial structure + -> (e -> e -> e) -- ^ update function + -> [(IndexOf c, e)] -- ^ association list + -> c e -- ^ result + +-------------------------------------------------------------------------- + +instance Container Vector Float where + scale = vectorMapValF Scale + scaleRecip = vectorMapValF Recip + addConstant = vectorMapValF AddConstant + add = vectorZipF Add + sub = vectorZipF Sub + mul = vectorZipF Mul + divide = vectorZipF Div + equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0 + arctan2 = vectorZipF ATan2 + scalar x = fromList [x] + konst' = constantD + build' = buildV + conj = id + cmap = mapVector + atIndex = (@>) + minIndex = emptyErrorV "minIndex" (round . toScalarF MinIdx) + maxIndex = emptyErrorV "maxIndex" (round . toScalarF MaxIdx) + minElement = emptyErrorV "minElement" (toScalarF Min) + maxElement = emptyErrorV "maxElement" (toScalarF Max) + sumElements = sumF + prodElements = prodF + step = stepF + find = findV + assoc = assocV + accum = accumV + cond = condV condF + +instance Container Vector Double where + scale = vectorMapValR Scale + scaleRecip = vectorMapValR Recip + addConstant = vectorMapValR AddConstant + add = vectorZipR Add + sub = vectorZipR Sub + mul = vectorZipR Mul + divide = vectorZipR Div + equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0 + arctan2 = vectorZipR ATan2 + scalar x = fromList [x] + konst' = constantD + build' = buildV + conj = id + cmap = mapVector + atIndex = (@>) + minIndex = emptyErrorV "minIndex" (round . toScalarR MinIdx) + maxIndex = emptyErrorV "maxIndex" (round . toScalarR MaxIdx) + minElement = emptyErrorV "minElement" (toScalarR Min) + maxElement = emptyErrorV "maxElement" (toScalarR Max) + sumElements = sumR + prodElements = prodR + step = stepD + find = findV + assoc = assocV + accum = accumV + cond = condV condD + +instance Container Vector (Complex Double) where + scale = vectorMapValC Scale + scaleRecip = vectorMapValC Recip + addConstant = vectorMapValC AddConstant + add = vectorZipC Add + sub = vectorZipC Sub + mul = vectorZipC Mul + divide = vectorZipC Div + equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 + arctan2 = vectorZipC ATan2 + scalar x = fromList [x] + konst' = constantD + build' = buildV + conj = conjugateC + cmap = mapVector + atIndex = (@>) + minIndex = emptyErrorV "minIndex" (minIndex . fst . fromComplex . (mul <*> conj)) + maxIndex = emptyErrorV "maxIndex" (maxIndex . fst . fromComplex . (mul <*> conj)) + minElement = emptyErrorV "minElement" (atIndex <*> minIndex) + maxElement = emptyErrorV "maxElement" (atIndex <*> maxIndex) + sumElements = sumC + prodElements = prodC + step = undefined -- cannot match + find = findV + assoc = assocV + accum = accumV + cond = undefined -- cannot match + +instance Container Vector (Complex Float) where + scale = vectorMapValQ Scale + scaleRecip = vectorMapValQ Recip + addConstant = vectorMapValQ AddConstant + add = vectorZipQ Add + sub = vectorZipQ Sub + mul = vectorZipQ Mul + divide = vectorZipQ Div + equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 + arctan2 = vectorZipQ ATan2 + scalar x = fromList [x] + konst' = constantD + build' = buildV + conj = conjugateQ + cmap = mapVector + atIndex = (@>) + minIndex = emptyErrorV "minIndex" (minIndex . fst . fromComplex . (mul <*> conj)) + maxIndex = emptyErrorV "maxIndex" (maxIndex . fst . fromComplex . (mul <*> conj)) + minElement = emptyErrorV "minElement" (atIndex <*> minIndex) + maxElement = emptyErrorV "maxElement" (atIndex <*> maxIndex) + sumElements = sumQ + prodElements = prodQ + step = undefined -- cannot match + find = findV + assoc = assocV + accum = accumV + cond = undefined -- cannot match + +--------------------------------------------------------------- + +instance (Container Vector a) => Container Matrix a where + scale x = liftMatrix (scale x) + scaleRecip x = liftMatrix (scaleRecip x) + addConstant x = liftMatrix (addConstant x) + add = liftMatrix2 add + sub = liftMatrix2 sub + mul = liftMatrix2 mul + divide = liftMatrix2 divide + equal a b = cols a == cols b && flatten a `equal` flatten b + arctan2 = liftMatrix2 arctan2 + scalar x = (1><1) [x] + konst' v (r,c) = matrixFromVector RowMajor r c (konst' v (r*c)) + build' = buildM + conj = liftMatrix conj + cmap f = liftMatrix (mapVector f) + atIndex = (@@>) + minIndex = emptyErrorM "minIndex of Matrix" $ + \m -> divMod (minIndex $ flatten m) (cols m) + maxIndex = emptyErrorM "maxIndex of Matrix" $ + \m -> divMod (maxIndex $ flatten m) (cols m) + minElement = emptyErrorM "minElement of Matrix" (atIndex <*> minIndex) + maxElement = emptyErrorM "maxElement of Matrix" (atIndex <*> maxIndex) + sumElements = sumElements . flatten + prodElements = prodElements . flatten + step = liftMatrix step + find = findM + assoc = assocM + accum = accumM + cond = condM + + +emptyErrorV msg f v = + if dim v > 0 + then f v + else error $ msg ++ " of Vector with dim = 0" + +emptyErrorM msg f m = + if rows m > 0 && cols m > 0 + then f m + else error $ msg++" "++shSize m + +---------------------------------------------------- + +-- | Matrix product and related functions +class (Num e, Element e) => Product e where + -- | matrix product + multiply :: Matrix e -> Matrix e -> Matrix e + -- | sum of absolute value of elements (differs in complex case from @norm1@) + absSum :: Vector e -> RealOf e + -- | sum of absolute value of elements + norm1 :: Vector e -> RealOf e + -- | euclidean norm + norm2 :: Vector e -> RealOf e + -- | element of maximum magnitude + normInf :: Vector e -> RealOf e + +instance Product Float where + norm2 = emptyVal (toScalarF Norm2) + absSum = emptyVal (toScalarF AbsSum) + norm1 = emptyVal (toScalarF AbsSum) + normInf = emptyVal (maxElement . vectorMapF Abs) + multiply = emptyMul multiplyF + +instance Product Double where + norm2 = emptyVal (toScalarR Norm2) + absSum = emptyVal (toScalarR AbsSum) + norm1 = emptyVal (toScalarR AbsSum) + normInf = emptyVal (maxElement . vectorMapR Abs) + multiply = emptyMul multiplyR + +instance Product (Complex Float) where + norm2 = emptyVal (toScalarQ Norm2) + absSum = emptyVal (toScalarQ AbsSum) + norm1 = emptyVal (sumElements . fst . fromComplex . vectorMapQ Abs) + normInf = emptyVal (maxElement . fst . fromComplex . vectorMapQ Abs) + multiply = emptyMul multiplyQ + +instance Product (Complex Double) where + norm2 = emptyVal (toScalarC Norm2) + absSum = emptyVal (toScalarC AbsSum) + norm1 = emptyVal (sumElements . fst . fromComplex . vectorMapC Abs) + normInf = emptyVal (maxElement . fst . fromComplex . vectorMapC Abs) + multiply = emptyMul multiplyC + +emptyMul m a b + | x1 == 0 && x2 == 0 || r == 0 || c == 0 = konst' 0 (r,c) + | otherwise = m a b + where + r = rows a + x1 = cols a + x2 = rows b + c = cols b + +emptyVal f v = + if dim v > 0 + then f v + else 0 + +-- FIXME remove unused C wrappers +-- | (unconjugated) dot product +udot :: Product e => Vector e -> Vector e -> e +udot u v + | dim u == dim v = val (asRow u `multiply` asColumn v) + | otherwise = error $ "different dimensions "++show (dim u)++" and "++show (dim v)++" in dot product" + where + val m | dim u > 0 = m@@>(0,0) + | otherwise = 0 + +---------------------------------------------------------- + +-- synonym for matrix product +mXm :: Product t => Matrix t -> Matrix t -> Matrix t +mXm = multiply + +-- matrix - vector product +mXv :: Product t => Matrix t -> Vector t -> Vector t +mXv m v = flatten $ m `mXm` (asColumn v) + +-- vector - matrix product +vXm :: Product t => Vector t -> Matrix t -> Vector t +vXm v m = flatten $ (asRow v) `mXm` m + +{- | Outer product of two vectors. + +>>> fromList [1,2,3] `outer` fromList [5,2,3] +(3><3) + [ 5.0, 2.0, 3.0 + , 10.0, 4.0, 6.0 + , 15.0, 6.0, 9.0 ] + +-} +outer :: (Product t) => Vector t -> Vector t -> Matrix t +outer u v = asColumn u `multiply` asRow v + +{- | Kronecker product of two matrices. + +@m1=(2><3) + [ 1.0, 2.0, 0.0 + , 0.0, -1.0, 3.0 ] +m2=(4><3) + [ 1.0, 2.0, 3.0 + , 4.0, 5.0, 6.0 + , 7.0, 8.0, 9.0 + , 10.0, 11.0, 12.0 ]@ + +>>> kronecker m1 m2 +(8><9) + [ 1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 0.0, 0.0, 0.0 + , 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 0.0, 0.0, 0.0 + , 7.0, 8.0, 9.0, 14.0, 16.0, 18.0, 0.0, 0.0, 0.0 + , 10.0, 11.0, 12.0, 20.0, 22.0, 24.0, 0.0, 0.0, 0.0 + , 0.0, 0.0, 0.0, -1.0, -2.0, -3.0, 3.0, 6.0, 9.0 + , 0.0, 0.0, 0.0, -4.0, -5.0, -6.0, 12.0, 15.0, 18.0 + , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0 + , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ] + +-} +kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t +kronecker a b = fromBlocks + . splitEvery (cols a) + . map (reshape (cols b)) + . toRows + $ flatten a `outer` flatten b + +------------------------------------------------------------------- + + +class Convert t where + real :: Container c t => c (RealOf t) -> c t + complex :: Container c t => c t -> c (ComplexOf t) + single :: Container c t => c t -> c (SingleOf t) + double :: Container c t => c t -> c (DoubleOf t) + toComplex :: (Container c t, RealElement t) => (c t, c t) -> c (Complex t) + fromComplex :: (Container c t, RealElement t) => c (Complex t) -> (c t, c t) + + +instance Convert Double where + real = id + complex = comp' + single = single' + double = id + toComplex = toComplex' + fromComplex = fromComplex' + +instance Convert Float where + real = id + complex = comp' + single = id + double = double' + toComplex = toComplex' + fromComplex = fromComplex' + +instance Convert (Complex Double) where + real = comp' + complex = id + single = single' + double = id + toComplex = toComplex' + fromComplex = fromComplex' + +instance Convert (Complex Float) where + real = comp' + complex = id + single = id + double = double' + toComplex = toComplex' + fromComplex = fromComplex' + +------------------------------------------------------------------- + +type family RealOf x + +type instance RealOf Double = Double +type instance RealOf (Complex Double) = Double + +type instance RealOf Float = Float +type instance RealOf (Complex Float) = Float + +type family ComplexOf x + +type instance ComplexOf Double = Complex Double +type instance ComplexOf (Complex Double) = Complex Double + +type instance ComplexOf Float = Complex Float +type instance ComplexOf (Complex Float) = Complex Float + +type family SingleOf x + +type instance SingleOf Double = Float +type instance SingleOf Float = Float + +type instance SingleOf (Complex a) = Complex (SingleOf a) + +type family DoubleOf x + +type instance DoubleOf Double = Double +type instance DoubleOf Float = Double + +type instance DoubleOf (Complex a) = Complex (DoubleOf a) + +type family ElementOf c + +type instance ElementOf (Vector a) = a +type instance ElementOf (Matrix a) = a + +------------------------------------------------------------ + +buildM (rc,cc) f = fromLists [ [f r c | c <- cs] | r <- rs ] + where rs = map fromIntegral [0 .. (rc-1)] + cs = map fromIntegral [0 .. (cc-1)] + +buildV n f = fromList [f k | k <- ks] + where ks = map fromIntegral [0 .. (n-1)] + +-------------------------------------------------------- +-- | conjugate transpose +ctrans :: (Container Vector e, Element e) => Matrix e -> Matrix e +ctrans = liftMatrix conj . trans + +-- | Creates a square matrix with a given diagonal. +diag :: (Num a, Element a) => Vector a -> Matrix a +diag v = diagRect 0 v n n where n = dim v + +-- | creates the identity matrix of given dimension +ident :: (Num a, Element a) => Int -> Matrix a +ident n = diag (constantD 1 n) + +-------------------------------------------------------- + +findV p x = foldVectorWithIndex g [] x where + g k z l = if p z then k:l else l + +findM p x = map ((`divMod` cols x)) $ findV p (flatten x) + +assocV n z xs = ST.runSTVector $ do + v <- ST.newVector z n + mapM_ (\(k,x) -> ST.writeVector v k x) xs + return v + +assocM (r,c) z xs = ST.runSTMatrix $ do + m <- ST.newMatrix z r c + mapM_ (\((i,j),x) -> ST.writeMatrix m i j x) xs + return m + +accumV v0 f xs = ST.runSTVector $ do + v <- ST.thawVector v0 + mapM_ (\(k,x) -> ST.modifyVector v k (f x)) xs + return v + +accumM m0 f xs = ST.runSTMatrix $ do + m <- ST.thawMatrix m0 + mapM_ (\((i,j),x) -> ST.modifyMatrix m i j (f x)) xs + return m + +---------------------------------------------------------------------- + +condM a b l e t = matrixFromVector RowMajor (rows a'') (cols a'') $ cond a' b' l' e' t' + where + args@(a'':_) = conformMs [a,b,l,e,t] + [a', b', l', e', t'] = map flatten args + +condV f a b l e t = f a' b' l' e' t' + where + [a', b', l', e', t'] = conformVs [a,b,l,e,t] + diff --git a/packages/hmatrix/src/Numeric/Conversion.hs b/packages/hmatrix/src/Numeric/Conversion.hs new file mode 100644 index 0000000..8941451 --- /dev/null +++ b/packages/hmatrix/src/Numeric/Conversion.hs @@ -0,0 +1,91 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE UndecidableInstances #-} + +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Conversion +-- Copyright : (c) Alberto Ruiz 2010 +-- License : GPL-style +-- +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- Portability : portable +-- +-- Conversion routines +-- +----------------------------------------------------------------------------- + +module Numeric.Conversion ( + Complexable(..), RealElement, + module Data.Complex +) where + +import Data.Packed.Internal.Vector +import Data.Packed.Internal.Matrix +import Data.Complex +import Control.Arrow((***)) + +------------------------------------------------------------------- + +-- | Supported single-double precision type pairs +class (Element s, Element d) => Precision s d | s -> d, d -> s where + double2FloatG :: Vector d -> Vector s + float2DoubleG :: Vector s -> Vector d + +instance Precision Float Double where + double2FloatG = double2FloatV + float2DoubleG = float2DoubleV + +instance Precision (Complex Float) (Complex Double) where + double2FloatG = asComplex . double2FloatV . asReal + float2DoubleG = asComplex . float2DoubleV . asReal + +-- | Supported real types +class (Element t, Element (Complex t), RealFloat t +-- , RealOf t ~ t, RealOf (Complex t) ~ t + ) + => RealElement t + +instance RealElement Double +instance RealElement Float + + +-- | Structures that may contain complex numbers +class Complexable c where + toComplex' :: (RealElement e) => (c e, c e) -> c (Complex e) + fromComplex' :: (RealElement e) => c (Complex e) -> (c e, c e) + comp' :: (RealElement e) => c e -> c (Complex e) + single' :: Precision a b => c b -> c a + double' :: Precision a b => c a -> c b + + +instance Complexable Vector where + toComplex' = toComplexV + fromComplex' = fromComplexV + comp' v = toComplex' (v,constantD 0 (dim v)) + single' = double2FloatG + double' = float2DoubleG + + +-- | creates a complex vector from vectors with real and imaginary parts +toComplexV :: (RealElement a) => (Vector a, Vector a) -> Vector (Complex a) +toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i] + +-- | the inverse of 'toComplex' +fromComplexV :: (RealElement a) => Vector (Complex a) -> (Vector a, Vector a) +fromComplexV z = (r,i) where + [r,i] = toColumns $ reshape 2 $ asReal z + + +instance Complexable Matrix where + toComplex' = uncurry $ liftMatrix2 $ curry toComplex' + fromComplex' z = (reshape c *** reshape c) . fromComplex' . flatten $ z + where c = cols z + comp' = liftMatrix comp' + single' = liftMatrix single' + double' = liftMatrix double' + diff --git a/packages/hmatrix/src/Numeric/GSL.hs b/packages/hmatrix/src/Numeric/GSL.hs new file mode 100644 index 0000000..5f39a3e --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL.hs @@ -0,0 +1,43 @@ +{- | + +Module : Numeric.GSL +Copyright : (c) Alberto Ruiz 2006-7 +License : GPL-style + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional +Portability : uses -fffi and -fglasgow-exts + +This module reexports all available GSL functions. + +The GSL special functions are in the separate package hmatrix-special. + +-} + +module Numeric.GSL ( + module Numeric.GSL.Integration +, module Numeric.GSL.Differentiation +, module Numeric.GSL.Fourier +, module Numeric.GSL.Polynomials +, module Numeric.GSL.Minimization +, module Numeric.GSL.Root +, module Numeric.GSL.ODE +, module Numeric.GSL.Fitting +, module Data.Complex +, setErrorHandlerOff +) where + +import Numeric.GSL.Integration +import Numeric.GSL.Differentiation +import Numeric.GSL.Fourier +import Numeric.GSL.Polynomials +import Numeric.GSL.Minimization +import Numeric.GSL.Root +import Numeric.GSL.ODE +import Numeric.GSL.Fitting +import Data.Complex + + +-- | This action removes the GSL default error handler (which aborts the program), so that +-- GSL errors can be handled by Haskell (using Control.Exception) and ghci doesn't abort. +foreign import ccall unsafe "GSL/gsl-aux.h no_abort_on_error" setErrorHandlerOff :: IO () diff --git a/packages/hmatrix/src/Numeric/GSL/Differentiation.hs b/packages/hmatrix/src/Numeric/GSL/Differentiation.hs new file mode 100644 index 0000000..93c5007 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Differentiation.hs @@ -0,0 +1,87 @@ +{-# OPTIONS #-} +----------------------------------------------------------------------------- +{- | +Module : Numeric.GSL.Differentiation +Copyright : (c) Alberto Ruiz 2006 +License : GPL-style + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional +Portability : uses ffi + +Numerical differentiation. + + + +From the GSL manual: \"The functions described in this chapter compute numerical derivatives by finite differencing. An adaptive algorithm is used to find the best choice of finite difference and to estimate the error in the derivative.\" +-} +----------------------------------------------------------------------------- +module Numeric.GSL.Differentiation ( + derivCentral, + derivForward, + derivBackward +) where + +import Foreign.C.Types +import Foreign.Marshal.Alloc(malloc, free) +import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr) +import Foreign.Storable(peek) +import Data.Packed.Internal(check,(//)) +import System.IO.Unsafe(unsafePerformIO) + +derivGen :: + CInt -- ^ type: 0 central, 1 forward, 2 backward + -> Double -- ^ initial step size + -> (Double -> Double) -- ^ function + -> Double -- ^ point where the derivative is taken + -> (Double, Double) -- ^ result and error +derivGen c h f x = unsafePerformIO $ do + r <- malloc + e <- malloc + fp <- mkfun (\y _ -> f y) + c_deriv c fp x h r e // check "deriv" + vr <- peek r + ve <- peek e + let result = (vr,ve) + free r + free e + freeHaskellFunPtr fp + return result + +foreign import ccall safe "gsl-aux.h deriv" + c_deriv :: CInt -> FunPtr (Double -> Ptr () -> Double) -> Double -> Double + -> Ptr Double -> Ptr Double -> IO CInt + + +{- | Adaptive central difference algorithm, /gsl_deriv_central/. For example: + +>>> let deriv = derivCentral 0.01 +>>> deriv sin (pi/4) +(0.7071067812000676,1.0600063101654055e-10) +>>> cos (pi/4) +0.7071067811865476 + +-} +derivCentral :: Double -- ^ initial step size + -> (Double -> Double) -- ^ function + -> Double -- ^ point where the derivative is taken + -> (Double, Double) -- ^ result and absolute error +derivCentral = derivGen 0 + +-- | Adaptive forward difference algorithm, /gsl_deriv_forward/. The function is evaluated only at points greater than x, and never at x itself. The derivative is returned in result and an estimate of its absolute error is returned in abserr. This function should be used if f(x) has a discontinuity at x, or is undefined for values less than x. A backward derivative can be obtained using a negative step. +derivForward :: Double -- ^ initial step size + -> (Double -> Double) -- ^ function + -> Double -- ^ point where the derivative is taken + -> (Double, Double) -- ^ result and absolute error +derivForward = derivGen 1 + +-- | Adaptive backward difference algorithm, /gsl_deriv_backward/. +derivBackward ::Double -- ^ initial step size + -> (Double -> Double) -- ^ function + -> Double -- ^ point where the derivative is taken + -> (Double, Double) -- ^ result and absolute error +derivBackward = derivGen 2 + +{- | conversion of Haskell functions into function pointers that can be used in the C side +-} +foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double)) diff --git a/packages/hmatrix/src/Numeric/GSL/Fitting.hs b/packages/hmatrix/src/Numeric/GSL/Fitting.hs new file mode 100644 index 0000000..c4f3a91 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Fitting.hs @@ -0,0 +1,179 @@ +{- | +Module : Numeric.GSL.Fitting +Copyright : (c) Alberto Ruiz 2010 +License : GPL + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional +Portability : uses ffi + +Nonlinear Least-Squares Fitting + + + +The example program in the GSL manual (see examples/fitting.hs): + +@ +dat = [ + ([0.0],([6.0133918608118675],0.1)), + ([1.0],([5.5153769909966535],0.1)), + ([2.0],([5.261094606015287],0.1)), + ... + ([39.0],([1.0619821710802808],0.1))] + +expModel [a,lambda,b] [t] = [a * exp (-lambda * t) + b] + +expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]] + +(sol,path) = fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0] +@ + +>>> path +(6><5) + [ 1.0, 76.45780563978782, 1.6465931240727802, 1.8147715267618197e-2, 0.6465931240727797 + , 2.0, 37.683816318260355, 2.858760367632973, 8.092094813253975e-2, 1.4479636296208662 + , 3.0, 9.5807893736187, 4.948995119561291, 0.11942927999921617, 1.0945766509238248 + , 4.0, 5.630494933603935, 5.021755718065913, 0.10287787128056883, 1.0338835440862608 + , 5.0, 5.443976278682909, 5.045204331329302, 0.10405523433131504, 1.019416067207375 + , 6.0, 5.4439736648994685, 5.045357818922331, 0.10404905846029407, 1.0192487112786812 ] +>>> sol +[(5.045357818922331,6.027976702418132e-2), +(0.10404905846029407,3.157045047172834e-3), +(1.0192487112786812,3.782067731353722e-2)] + +-} +----------------------------------------------------------------------------- + +module Numeric.GSL.Fitting ( + -- * Levenberg-Marquardt + nlFitting, FittingMethod(..), + -- * Utilities + fitModelScaled, fitModel +) where + +import Data.Packed.Internal +import Numeric.LinearAlgebra +import Numeric.GSL.Internal + +import Foreign.Ptr(FunPtr, freeHaskellFunPtr) +import Foreign.C.Types +import System.IO.Unsafe(unsafePerformIO) + +------------------------------------------------------------------------- + +data FittingMethod = LevenbergMarquardtScaled -- ^ Interface to gsl_multifit_fdfsolver_lmsder. This is a robust and efficient version of the Levenberg-Marquardt algorithm as implemented in the scaled lmder routine in minpack. Minpack was written by Jorge J. More, Burton S. Garbow and Kenneth E. Hillstrom. + | LevenbergMarquardt -- ^ This is an unscaled version of the lmder algorithm. The elements of the diagonal scaling matrix D are set to 1. This algorithm may be useful in circumstances where the scaled version of lmder converges too slowly, or the function is already scaled appropriately. + deriving (Enum,Eq,Show,Bounded) + + +-- | Nonlinear multidimensional least-squares fitting. +nlFitting :: FittingMethod + -> Double -- ^ absolute tolerance + -> Double -- ^ relative tolerance + -> Int -- ^ maximum number of iterations allowed + -> (Vector Double -> Vector Double) -- ^ function to be minimized + -> (Vector Double -> Matrix Double) -- ^ Jacobian + -> Vector Double -- ^ starting point + -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path + +nlFitting method epsabs epsrel maxit fun jac xinit = nlFitGen (fi (fromEnum method)) fun jac xinit epsabs epsrel maxit + +nlFitGen m f jac xiv epsabs epsrel maxit = unsafePerformIO $ do + let p = dim xiv + n = dim (f xiv) + fp <- mkVecVecfun (aux_vTov (checkdim1 n p . f)) + jp <- mkVecMatfun (aux_vTom (checkdim2 n p . jac)) + rawpath <- createMatrix RowMajor maxit (2+p) + app2 (c_nlfit m fp jp epsabs epsrel (fi maxit) (fi n)) vec xiv mat rawpath "c_nlfit" + let it = round (rawpath @@> (maxit-1,0)) + path = takeRows it rawpath + [sol] = toRows $ dropRows (it-1) path + freeHaskellFunPtr fp + freeHaskellFunPtr jp + return (subVector 2 p sol, path) + +foreign import ccall safe "nlfit" + c_nlfit:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> Double -> CInt -> CInt -> TVM + +------------------------------------------------------- + +checkdim1 n _p v + | dim v == n = v + | otherwise = error $ "Error: "++ show n + ++ " components expected in the result of the function supplied to nlFitting" + +checkdim2 n p m + | rows m == n && cols m == p = m + | otherwise = error $ "Error: "++ show n ++ "x" ++ show p + ++ " Jacobian expected in nlFitting" + +------------------------------------------------------------ + +err (model,deriv) dat vsol = zip sol errs where + sol = toList vsol + c = max 1 (chi/sqrt (fromIntegral dof)) + dof = length dat - (rows cov) + chi = norm2 (fromList $ cost (resMs model) dat sol) + js = fromLists $ jacobian (resDs deriv) dat sol + cov = inv $ trans js <> js + errs = toList $ scalar c * sqrt (takeDiag cov) + + + +-- | Higher level interface to 'nlFitting' 'LevenbergMarquardtScaled'. The optimization function and +-- Jacobian are automatically built from a model f vs x = y and its derivatives, and a list of +-- instances (x, (y,sigma)) to be fitted. + +fitModelScaled + :: Double -- ^ absolute tolerance + -> Double -- ^ relative tolerance + -> Int -- ^ maximum number of iterations allowed + -> ([Double] -> x -> [Double], [Double] -> x -> [[Double]]) -- ^ (model, derivatives) + -> [(x, ([Double], Double))] -- ^ instances + -> [Double] -- ^ starting point + -> ([(Double, Double)], Matrix Double) -- ^ (solution, error) and optimization path +fitModelScaled epsabs epsrel maxit (model,deriv) dt xin = (err (model,deriv) dt sol, path) where + (sol,path) = nlFitting LevenbergMarquardtScaled epsabs epsrel maxit + (fromList . cost (resMs model) dt . toList) + (fromLists . jacobian (resDs deriv) dt . toList) + (fromList xin) + + + +-- | Higher level interface to 'nlFitting' 'LevenbergMarquardt'. The optimization function and +-- Jacobian are automatically built from a model f vs x = y and its derivatives, and a list of +-- instances (x,y) to be fitted. + +fitModel :: Double -- ^ absolute tolerance + -> Double -- ^ relative tolerance + -> Int -- ^ maximum number of iterations allowed + -> ([Double] -> x -> [Double], [Double] -> x -> [[Double]]) -- ^ (model, derivatives) + -> [(x, [Double])] -- ^ instances + -> [Double] -- ^ starting point + -> ([Double], Matrix Double) -- ^ solution and optimization path +fitModel epsabs epsrel maxit (model,deriv) dt xin = (toList sol, path) where + (sol,path) = nlFitting LevenbergMarquardt epsabs epsrel maxit + (fromList . cost (resM model) dt . toList) + (fromLists . jacobian (resD deriv) dt . toList) + (fromList xin) + +cost model ds vs = concatMap (model vs) ds + +jacobian modelDer ds vs = concatMap (modelDer vs) ds + +-- | Model-to-residual for association pairs with sigma, to be used with 'fitModel'. +resMs :: ([Double] -> x -> [Double]) -> [Double] -> (x, ([Double], Double)) -> [Double] +resMs m v = \(x,(ys,s)) -> zipWith (g s) (m v x) ys where g s a b = (a-b)/s + +-- | Associated derivative for 'resMs'. +resDs :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, ([Double], Double)) -> [[Double]] +resDs m v = \(x,(_,s)) -> map (map (/s)) (m v x) + +-- | Model-to-residual for association pairs, to be used with 'fitModel'. It is equivalent +-- to 'resMs' with all sigmas = 1. +resM :: ([Double] -> x -> [Double]) -> [Double] -> (x, [Double]) -> [Double] +resM m v = \(x,ys) -> zipWith g (m v x) ys where g a b = a-b + +-- | Associated derivative for 'resM'. +resD :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, [Double]) -> [[Double]] +resD m v = \(x,_) -> m v x diff --git a/packages/hmatrix/src/Numeric/GSL/Fourier.hs b/packages/hmatrix/src/Numeric/GSL/Fourier.hs new file mode 100644 index 0000000..86aedd6 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Fourier.hs @@ -0,0 +1,47 @@ +{-# LANGUAGE ForeignFunctionInterface #-} +----------------------------------------------------------------------------- +{- | +Module : Numeric.GSL.Fourier +Copyright : (c) Alberto Ruiz 2006 +License : GPL-style + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional +Portability : uses ffi + +Fourier Transform. + + + +-} +----------------------------------------------------------------------------- +module Numeric.GSL.Fourier ( + fft, + ifft +) where + +import Data.Packed.Internal +import Data.Complex +import Foreign.C.Types +import System.IO.Unsafe (unsafePerformIO) + +genfft code v = unsafePerformIO $ do + r <- createVector (dim v) + app2 (c_fft code) vec v vec r "fft" + return r + +foreign import ccall unsafe "gsl-aux.h fft" c_fft :: CInt -> TCVCV + + +{- | Fast 1D Fourier transform of a 'Vector' @(@'Complex' 'Double'@)@ using /gsl_fft_complex_forward/. It uses the same scaling conventions as GNU Octave. + +>>> fft (fromList [1,2,3,4]) +fromList [10.0 :+ 0.0,(-2.0) :+ 2.0,(-2.0) :+ 0.0,(-2.0) :+ (-2.0)] + +-} +fft :: Vector (Complex Double) -> Vector (Complex Double) +fft = genfft 0 + +-- | The inverse of 'fft', using /gsl_fft_complex_inverse/. +ifft :: Vector (Complex Double) -> Vector (Complex Double) +ifft = genfft 1 diff --git a/packages/hmatrix/src/Numeric/GSL/Integration.hs b/packages/hmatrix/src/Numeric/GSL/Integration.hs new file mode 100644 index 0000000..5f0a415 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Integration.hs @@ -0,0 +1,250 @@ +{-# OPTIONS #-} +----------------------------------------------------------------------------- +{- | +Module : Numeric.GSL.Integration +Copyright : (c) Alberto Ruiz 2006 +License : GPL-style + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional +Portability : uses ffi + +Numerical integration routines. + + +-} +----------------------------------------------------------------------------- + +module Numeric.GSL.Integration ( + integrateQNG, + integrateQAGS, + integrateQAGI, + integrateQAGIU, + integrateQAGIL, + integrateCQUAD +) where + +import Foreign.C.Types +import Foreign.Marshal.Alloc(malloc, free) +import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr) +import Foreign.Storable(peek) +import Data.Packed.Internal(check,(//)) +import System.IO.Unsafe(unsafePerformIO) + +eps = 1e-12 + +{- | conversion of Haskell functions into function pointers that can be used in the C side +-} +foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double)) + +-------------------------------------------------------------------- +{- | Numerical integration using /gsl_integration_qags/ (adaptive integration with singularities). For example: + +>>> let quad = integrateQAGS 1E-9 1000 +>>> let f a x = x**(-0.5) * log (a*x) +>>> quad (f 1) 0 1 +(-3.999999999999974,4.871658632055187e-13) + +-} + +integrateQAGS :: Double -- ^ precision (e.g. 1E-9) + -> Int -- ^ size of auxiliary workspace (e.g. 1000) + -> (Double -> Double) -- ^ function to be integrated on the interval (a,b) + -> Double -- ^ a + -> Double -- ^ b + -> (Double, Double) -- ^ result of the integration and error +integrateQAGS prec n f a b = unsafePerformIO $ do + r <- malloc + e <- malloc + fp <- mkfun (\x _ -> f x) + c_integrate_qags fp a b eps prec (fromIntegral n) r e // check "integrate_qags" + vr <- peek r + ve <- peek e + let result = (vr,ve) + free r + free e + freeHaskellFunPtr fp + return result + +foreign import ccall safe "integrate_qags" c_integrate_qags + :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double + -> Double -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt + +----------------------------------------------------------------- +{- | Numerical integration using /gsl_integration_qng/ (useful for fast integration of smooth functions). For example: + +>>> let quad = integrateQNG 1E-6 +>>> quad (\x -> 4/(1+x*x)) 0 1 +(3.141592653589793,3.487868498008632e-14) + +-} + +integrateQNG :: Double -- ^ precision (e.g. 1E-9) + -> (Double -> Double) -- ^ function to be integrated on the interval (a,b) + -> Double -- ^ a + -> Double -- ^ b + -> (Double, Double) -- ^ result of the integration and error +integrateQNG prec f a b = unsafePerformIO $ do + r <- malloc + e <- malloc + fp <- mkfun (\x _ -> f x) + c_integrate_qng fp a b eps prec r e // check "integrate_qng" + vr <- peek r + ve <- peek e + let result = (vr,ve) + free r + free e + freeHaskellFunPtr fp + return result + +foreign import ccall safe "integrate_qng" c_integrate_qng + :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double + -> Double -> Double -> Ptr Double -> Ptr Double -> IO CInt + +-------------------------------------------------------------------- +{- | Numerical integration using /gsl_integration_qagi/ (integration over the infinite integral -Inf..Inf using QAGS). +For example: + +>>> let quad = integrateQAGI 1E-9 1000 +>>> let f a x = exp(-a * x^2) +>>> quad (f 0.5) +(2.5066282746310002,6.229215880648858e-11) + +-} + +integrateQAGI :: Double -- ^ precision (e.g. 1E-9) + -> Int -- ^ size of auxiliary workspace (e.g. 1000) + -> (Double -> Double) -- ^ function to be integrated on the interval (-Inf,Inf) + -> (Double, Double) -- ^ result of the integration and error +integrateQAGI prec n f = unsafePerformIO $ do + r <- malloc + e <- malloc + fp <- mkfun (\x _ -> f x) + c_integrate_qagi fp eps prec (fromIntegral n) r e // check "integrate_qagi" + vr <- peek r + ve <- peek e + let result = (vr,ve) + free r + free e + freeHaskellFunPtr fp + return result + +foreign import ccall safe "integrate_qagi" c_integrate_qagi + :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double + -> CInt -> Ptr Double -> Ptr Double -> IO CInt + +-------------------------------------------------------------------- +{- | Numerical integration using /gsl_integration_qagiu/ (integration over the semi-infinite integral a..Inf). +For example: + +>>> let quad = integrateQAGIU 1E-9 1000 +>>> let f a x = exp(-a * x^2) +>>> quad (f 0.5) 0 +(1.2533141373155001,3.114607940324429e-11) + +-} + +integrateQAGIU :: Double -- ^ precision (e.g. 1E-9) + -> Int -- ^ size of auxiliary workspace (e.g. 1000) + -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf) + -> Double -- ^ a + -> (Double, Double) -- ^ result of the integration and error +integrateQAGIU prec n f a = unsafePerformIO $ do + r <- malloc + e <- malloc + fp <- mkfun (\x _ -> f x) + c_integrate_qagiu fp a eps prec (fromIntegral n) r e // check "integrate_qagiu" + vr <- peek r + ve <- peek e + let result = (vr,ve) + free r + free e + freeHaskellFunPtr fp + return result + +foreign import ccall safe "integrate_qagiu" c_integrate_qagiu + :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double + -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt + +-------------------------------------------------------------------- +{- | Numerical integration using /gsl_integration_qagil/ (integration over the semi-infinite integral -Inf..b). +For example: + +>>> let quad = integrateQAGIL 1E-9 1000 +>>> let f a x = exp(-a * x^2) +>>> quad (f 0.5) 0 +(1.2533141373155001,3.114607940324429e-11) + +-} + +integrateQAGIL :: Double -- ^ precision (e.g. 1E-9) + -> Int -- ^ size of auxiliary workspace (e.g. 1000) + -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf) + -> Double -- ^ b + -> (Double, Double) -- ^ result of the integration and error +integrateQAGIL prec n f b = unsafePerformIO $ do + r <- malloc + e <- malloc + fp <- mkfun (\x _ -> f x) + c_integrate_qagil fp b eps prec (fromIntegral n) r e // check "integrate_qagil" + vr <- peek r + ve <- peek e + let result = (vr,ve) + free r + free e + freeHaskellFunPtr fp + return result + +foreign import ccall safe "gsl-aux.h integrate_qagil" c_integrate_qagil + :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double + -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt + + +-------------------------------------------------------------------- +{- | Numerical integration using /gsl_integration_cquad/ (quadrature +for general integrands). From the GSL manual: + +@CQUAD is a new doubly-adaptive general-purpose quadrature routine +which can handle most types of singularities, non-numerical function +values such as Inf or NaN, as well as some divergent integrals. It +generally requires more function evaluations than the integration +routines in QUADPACK, yet fails less often for difficult integrands.@ + +For example: + +>>> let quad = integrateCQUAD 1E-12 1000 +>>> let f a x = exp(-a * x^2) +>>> quad (f 0.5) 2 5 +(5.7025405463957006e-2,9.678874441303705e-16,95) + +Unlike other quadrature methods, integrateCQUAD also returns the +number of function evaluations required. + +-} + +integrateCQUAD :: Double -- ^ precision (e.g. 1E-9) + -> Int -- ^ size of auxiliary workspace (e.g. 1000) + -> (Double -> Double) -- ^ function to be integrated on the interval (a, b) + -> Double -- ^ a + -> Double -- ^ b + -> (Double, Double, Int) -- ^ result of the integration, error and number of function evaluations performed +integrateCQUAD prec n f a b = unsafePerformIO $ do + r <- malloc + e <- malloc + neval <- malloc + fp <- mkfun (\x _ -> f x) + c_integrate_cquad fp a b eps prec (fromIntegral n) r e neval // check "integrate_cquad" + vr <- peek r + ve <- peek e + vneval <- peek neval + let result = (vr,ve,vneval) + free r + free e + free neval + freeHaskellFunPtr fp + return result + +foreign import ccall safe "integrate_cquad" c_integrate_cquad + :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double + -> Double -> Double -> CInt -> Ptr Double -> Ptr Double -> Ptr Int -> IO CInt + diff --git a/packages/hmatrix/src/Numeric/GSL/Internal.hs b/packages/hmatrix/src/Numeric/GSL/Internal.hs new file mode 100644 index 0000000..69a9750 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Internal.hs @@ -0,0 +1,76 @@ +-- Module : Numeric.GSL.Internal +-- Copyright : (c) Alberto Ruiz 2009 +-- License : GPL +-- +-- Maintainer : Alberto Ruiz (aruiz at um dot es) +-- Stability : provisional +-- Portability : uses ffi +-- +-- Auxiliary functions. +-- +-- #hide + +module Numeric.GSL.Internal where + +import Data.Packed.Internal + +import Foreign.Marshal.Array(copyArray) +import Foreign.Ptr(Ptr, FunPtr) +import Foreign.C.Types +import System.IO.Unsafe(unsafePerformIO) + +iv :: (Vector Double -> Double) -> (CInt -> Ptr Double -> Double) +iv f n p = f (createV (fromIntegral n) copy "iv") where + copy n' q = do + copyArray q p (fromIntegral n') + return 0 + +-- | conversion of Haskell functions into function pointers that can be used in the C side +foreign import ccall safe "wrapper" + mkVecfun :: (CInt -> Ptr Double -> Double) + -> IO( FunPtr (CInt -> Ptr Double -> Double)) + +foreign import ccall safe "wrapper" + mkVecVecfun :: TVV -> IO (FunPtr TVV) + +foreign import ccall safe "wrapper" + mkDoubleVecVecfun :: (Double -> TVV) -> IO (FunPtr (Double -> TVV)) + +foreign import ccall safe "wrapper" + mkDoublefun :: (Double -> Double) -> IO (FunPtr (Double -> Double)) + +aux_vTov :: (Vector Double -> Vector Double) -> TVV +aux_vTov f n p nr r = g where + v = f x + x = createV (fromIntegral n) copy "aux_vTov" + copy n' q = do + copyArray q p (fromIntegral n') + return 0 + g = do unsafeWith v $ \p' -> copyArray r p' (fromIntegral nr) + return 0 + +foreign import ccall safe "wrapper" + mkVecMatfun :: TVM -> IO (FunPtr TVM) + +foreign import ccall safe "wrapper" + mkDoubleVecMatfun :: (Double -> TVM) -> IO (FunPtr (Double -> TVM)) + +aux_vTom :: (Vector Double -> Matrix Double) -> TVM +aux_vTom f n p rr cr r = g where + v = flatten $ f x + x = createV (fromIntegral n) copy "aux_vTov" + copy n' q = do + copyArray q p (fromIntegral n') + return 0 + g = do unsafeWith v $ \p' -> copyArray r p' (fromIntegral $ rr*cr) + return 0 + +createV n fun msg = unsafePerformIO $ do + r <- createVector n + app1 fun vec r msg + return r + +createMIO r c fun msg = do + res <- createMatrix RowMajor r c + app1 fun mat res msg + return res diff --git a/packages/hmatrix/src/Numeric/GSL/Minimization.hs b/packages/hmatrix/src/Numeric/GSL/Minimization.hs new file mode 100644 index 0000000..1879dab --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Minimization.hs @@ -0,0 +1,227 @@ +{-# LANGUAGE ForeignFunctionInterface #-} +----------------------------------------------------------------------------- +{- | +Module : Numeric.GSL.Minimization +Copyright : (c) Alberto Ruiz 2006-9 +License : GPL-style + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional +Portability : uses ffi + +Minimization of a multidimensional function using some of the algorithms described in: + + + +The example in the GSL manual: + +@ +f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30 + +main = do + let (s,p) = minimize NMSimplex2 1E-2 30 [1,1] f [5,7] + print s + print p +@ + +>>> main +[0.9920430849306288,1.9969168063253182] + 0.000 512.500 1.130 6.500 5.000 + 1.000 290.625 1.409 5.250 4.000 + 2.000 290.625 1.409 5.250 4.000 + 3.000 252.500 1.409 5.500 1.000 + ... +22.000 30.001 0.013 0.992 1.997 +23.000 30.001 0.008 0.992 1.997 + +The path to the solution can be graphically shown by means of: + +@'Graphics.Plot.mplot' $ drop 3 ('toColumns' p)@ + +Taken from the GSL manual: + +The vector Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm is a quasi-Newton method which builds up an approximation to the second derivatives of the function f using the difference between successive gradient vectors. By combining the first and second derivatives the algorithm is able to take Newton-type steps towards the function minimum, assuming quadratic behavior in that region. + +The bfgs2 version of this minimizer is the most efficient version available, and is a faithful implementation of the line minimization scheme described in Fletcher's Practical Methods of Optimization, Algorithms 2.6.2 and 2.6.4. It supercedes the original bfgs routine and requires substantially fewer function and gradient evaluations. The user-supplied tolerance tol corresponds to the parameter \sigma used by Fletcher. A value of 0.1 is recommended for typical use (larger values correspond to less accurate line searches). + +The nmsimplex2 version is a new O(N) implementation of the earlier O(N^2) nmsimplex minimiser. It calculates the size of simplex as the rms distance of each vertex from the center rather than the mean distance, which has the advantage of allowing a linear update. + +-} + +----------------------------------------------------------------------------- +module Numeric.GSL.Minimization ( + minimize, minimizeV, MinimizeMethod(..), + minimizeD, minimizeVD, MinimizeMethodD(..), + uniMinimize, UniMinimizeMethod(..), + + minimizeNMSimplex, + minimizeConjugateGradient, + minimizeVectorBFGS2 +) where + + +import Data.Packed.Internal +import Data.Packed.Matrix +import Numeric.GSL.Internal + +import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr) +import Foreign.C.Types +import System.IO.Unsafe(unsafePerformIO) + +------------------------------------------------------------------------ + +{-# DEPRECATED minimizeNMSimplex "use minimize NMSimplex2 eps maxit sizes f xi" #-} +minimizeNMSimplex f xi szs eps maxit = minimize NMSimplex eps maxit szs f xi + +{-# DEPRECATED minimizeConjugateGradient "use minimizeD ConjugateFR eps maxit step tol f g xi" #-} +minimizeConjugateGradient step tol eps maxit f g xi = minimizeD ConjugateFR eps maxit step tol f g xi + +{-# DEPRECATED minimizeVectorBFGS2 "use minimizeD VectorBFGS2 eps maxit step tol f g xi" #-} +minimizeVectorBFGS2 step tol eps maxit f g xi = minimizeD VectorBFGS2 eps maxit step tol f g xi + +------------------------------------------------------------------------- + +data UniMinimizeMethod = GoldenSection + | BrentMini + | QuadGolden + deriving (Enum, Eq, Show, Bounded) + +-- | Onedimensional minimization. + +uniMinimize :: UniMinimizeMethod -- ^ The method used. + -> Double -- ^ desired precision of the solution + -> Int -- ^ maximum number of iterations allowed + -> (Double -> Double) -- ^ function to minimize + -> Double -- ^ guess for the location of the minimum + -> Double -- ^ lower bound of search interval + -> Double -- ^ upper bound of search interval + -> (Double, Matrix Double) -- ^ solution and optimization path + +uniMinimize method epsrel maxit fun xmin xl xu = uniMinimizeGen (fi (fromEnum method)) fun xmin xl xu epsrel maxit + +uniMinimizeGen m f xmin xl xu epsrel maxit = unsafePerformIO $ do + fp <- mkDoublefun f + rawpath <- createMIO maxit 4 + (c_uniMinize m fp epsrel (fi maxit) xmin xl xu) + "uniMinimize" + let it = round (rawpath @@> (maxit-1,0)) + path = takeRows it rawpath + [sol] = toLists $ dropRows (it-1) path + freeHaskellFunPtr fp + return (sol !! 1, path) + + +foreign import ccall safe "uniMinimize" + c_uniMinize:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> Double -> TM + +data MinimizeMethod = NMSimplex + | NMSimplex2 + deriving (Enum,Eq,Show,Bounded) + +-- | Minimization without derivatives +minimize :: MinimizeMethod + -> Double -- ^ desired precision of the solution (size test) + -> Int -- ^ maximum number of iterations allowed + -> [Double] -- ^ sizes of the initial search box + -> ([Double] -> Double) -- ^ function to minimize + -> [Double] -- ^ starting point + -> ([Double], Matrix Double) -- ^ solution vector and optimization path + +-- | Minimization without derivatives (vector version) +minimizeV :: MinimizeMethod + -> Double -- ^ desired precision of the solution (size test) + -> Int -- ^ maximum number of iterations allowed + -> Vector Double -- ^ sizes of the initial search box + -> (Vector Double -> Double) -- ^ function to minimize + -> Vector Double -- ^ starting point + -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path + +minimize method eps maxit sz f xi = v2l $ minimizeV method eps maxit (fromList sz) (f.toList) (fromList xi) + where v2l (v,m) = (toList v, m) + +ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2 + +minimizeV method eps maxit szv f xiv = unsafePerformIO $ do + let n = dim xiv + fp <- mkVecfun (iv f) + rawpath <- ww2 vec xiv vec szv $ \xiv' szv' -> + createMIO maxit (n+3) + (c_minimize (fi (fromEnum method)) fp eps (fi maxit) // xiv' // szv') + "minimize" + let it = round (rawpath @@> (maxit-1,0)) + path = takeRows it rawpath + sol = cdat $ dropColumns 3 $ dropRows (it-1) path + freeHaskellFunPtr fp + return (sol, path) + + +foreign import ccall safe "gsl-aux.h minimize" + c_minimize:: CInt -> FunPtr (CInt -> Ptr Double -> Double) -> Double -> CInt -> TVVM + +---------------------------------------------------------------------------------- + + +data MinimizeMethodD = ConjugateFR + | ConjugatePR + | VectorBFGS + | VectorBFGS2 + | SteepestDescent + deriving (Enum,Eq,Show,Bounded) + +-- | Minimization with derivatives. +minimizeD :: MinimizeMethodD + -> Double -- ^ desired precision of the solution (gradient test) + -> Int -- ^ maximum number of iterations allowed + -> Double -- ^ size of the first trial step + -> Double -- ^ tol (precise meaning depends on method) + -> ([Double] -> Double) -- ^ function to minimize + -> ([Double] -> [Double]) -- ^ gradient + -> [Double] -- ^ starting point + -> ([Double], Matrix Double) -- ^ solution vector and optimization path + +-- | Minimization with derivatives (vector version) +minimizeVD :: MinimizeMethodD + -> Double -- ^ desired precision of the solution (gradient test) + -> Int -- ^ maximum number of iterations allowed + -> Double -- ^ size of the first trial step + -> Double -- ^ tol (precise meaning depends on method) + -> (Vector Double -> Double) -- ^ function to minimize + -> (Vector Double -> Vector Double) -- ^ gradient + -> Vector Double -- ^ starting point + -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path + +minimizeD method eps maxit istep tol f df xi = v2l $ minimizeVD + method eps maxit istep tol (f.toList) (fromList.df.toList) (fromList xi) + where v2l (v,m) = (toList v, m) + + +minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do + let n = dim xiv + f' = f + df' = (checkdim1 n . df) + fp <- mkVecfun (iv f') + dfp <- mkVecVecfun (aux_vTov df') + rawpath <- vec xiv $ \xiv' -> + createMIO maxit (n+2) + (c_minimizeD (fi (fromEnum method)) fp dfp istep tol eps (fi maxit) // xiv') + "minimizeD" + let it = round (rawpath @@> (maxit-1,0)) + path = takeRows it rawpath + sol = cdat $ dropColumns 2 $ dropRows (it-1) path + freeHaskellFunPtr fp + freeHaskellFunPtr dfp + return (sol,path) + +foreign import ccall safe "gsl-aux.h minimizeD" + c_minimizeD :: CInt + -> FunPtr (CInt -> Ptr Double -> Double) + -> FunPtr TVV + -> Double -> Double -> Double -> CInt + -> TVM + +--------------------------------------------------------------------- + +checkdim1 n v + | dim v == n = v + | otherwise = error $ "Error: "++ show n + ++ " components expected in the result of the gradient supplied to minimizeD" diff --git a/packages/hmatrix/src/Numeric/GSL/ODE.hs b/packages/hmatrix/src/Numeric/GSL/ODE.hs new file mode 100644 index 0000000..9a29085 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/ODE.hs @@ -0,0 +1,138 @@ +{- | +Module : Numeric.GSL.ODE +Copyright : (c) Alberto Ruiz 2010 +License : GPL + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional +Portability : uses ffi + +Solution of ordinary differential equation (ODE) initial value problems. + + + +A simple example: + +@ +import Numeric.GSL.ODE +import Numeric.LinearAlgebra +import Numeric.LinearAlgebra.Util(mplot) + +xdot t [x,v] = [v, -0.95*x - 0.1*v] + +ts = linspace 100 (0,20 :: Double) + +sol = odeSolve xdot [10,0] ts + +main = mplot (ts : toColumns sol) +@ + +-} +----------------------------------------------------------------------------- + +module Numeric.GSL.ODE ( + odeSolve, odeSolveV, ODEMethod(..), Jacobian +) where + +import Data.Packed.Internal +import Numeric.GSL.Internal + +import Foreign.Ptr(FunPtr, nullFunPtr, freeHaskellFunPtr) +import Foreign.C.Types +import System.IO.Unsafe(unsafePerformIO) + +------------------------------------------------------------------------- + +type Jacobian = Double -> Vector Double -> Matrix Double + +-- | Stepping functions +data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method. + | RK4 -- ^ 4th order (classical) Runge-Kutta. The error estimate is obtained by halving the step-size. For more efficient estimate of the error, use the embedded methods. + | RKf45 -- ^ Embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good general-purpose integrator. + | RKck -- ^ Embedded Runge-Kutta Cash-Karp (4, 5) method. + | RK8pd -- ^ Embedded Runge-Kutta Prince-Dormand (8,9) method. + | RK2imp Jacobian -- ^ Implicit 2nd order Runge-Kutta at Gaussian points. + | RK4imp Jacobian -- ^ Implicit 4th order Runge-Kutta at Gaussian points. + | BSimp Jacobian -- ^ Implicit Bulirsch-Stoer method of Bader and Deuflhard. The method is generally suitable for stiff problems. + | RK1imp Jacobian -- ^ Implicit Gaussian first order Runge-Kutta. Also known as implicit Euler or backward Euler method. Error estimation is carried out by the step doubling method. + | MSAdams -- ^ A variable-coefficient linear multistep Adams method in Nordsieck form. This stepper uses explicit Adams-Bashforth (predictor) and implicit Adams-Moulton (corrector) methods in P(EC)^m functional iteration mode. Method order varies dynamically between 1 and 12. + | MSBDF Jacobian -- ^ A variable-coefficient linear multistep backward differentiation formula (BDF) method in Nordsieck form. This stepper uses the explicit BDF formula as predictor and implicit BDF formula as corrector. A modified Newton iteration method is used to solve the system of non-linear equations. Method order varies dynamically between 1 and 5. The method is generally suitable for stiff problems. + + +-- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists. +odeSolve + :: (Double -> [Double] -> [Double]) -- ^ xdot(t,x) + -> [Double] -- ^ initial conditions + -> Vector Double -- ^ desired solution times + -> Matrix Double -- ^ solution +odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts + where hi = (ts@>1 - ts@>0)/100 + epsAbs = 1.49012e-08 + epsRel = 1.49012e-08 + l2v f = \t -> fromList . f t . toList + +-- | Evolution of the system with adaptive step-size control. +odeSolveV + :: ODEMethod + -> Double -- ^ initial step size + -> Double -- ^ absolute tolerance for the state vector + -> Double -- ^ relative tolerance for the state vector + -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x) + -> Vector Double -- ^ initial conditions + -> Vector Double -- ^ desired solution times + -> Matrix Double -- ^ solution +odeSolveV RK2 = odeSolveV' 0 Nothing +odeSolveV RK4 = odeSolveV' 1 Nothing +odeSolveV RKf45 = odeSolveV' 2 Nothing +odeSolveV RKck = odeSolveV' 3 Nothing +odeSolveV RK8pd = odeSolveV' 4 Nothing +odeSolveV (RK2imp jac) = odeSolveV' 5 (Just jac) +odeSolveV (RK4imp jac) = odeSolveV' 6 (Just jac) +odeSolveV (BSimp jac) = odeSolveV' 7 (Just jac) +odeSolveV (RK1imp jac) = odeSolveV' 8 (Just jac) +odeSolveV MSAdams = odeSolveV' 9 Nothing +odeSolveV (MSBDF jac) = odeSolveV' 10 (Just jac) + + +odeSolveV' + :: CInt + -> Maybe (Double -> Vector Double -> Matrix Double) -- ^ optional jacobian + -> Double -- ^ initial step size + -> Double -- ^ absolute tolerance for the state vector + -> Double -- ^ relative tolerance for the state vector + -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x) + -> Vector Double -- ^ initial conditions + -> Vector Double -- ^ desired solution times + -> Matrix Double -- ^ solution +odeSolveV' method mbjac h epsAbs epsRel f xiv ts = unsafePerformIO $ do + let n = dim xiv + fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t)) + jp <- case mbjac of + Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t)) + Nothing -> return nullFunPtr + sol <- vec xiv $ \xiv' -> + vec (checkTimes ts) $ \ts' -> + createMIO (dim ts) n + (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' ) + "ode" + freeHaskellFunPtr fp + return sol + +foreign import ccall safe "ode" + ode_c :: CInt -> Double -> Double -> Double -> FunPtr (Double -> TVV) -> FunPtr (Double -> TVM) -> TVVM + +------------------------------------------------------- + +checkdim1 n v + | dim v == n = v + | otherwise = error $ "Error: "++ show n + ++ " components expected in the result of the function supplied to odeSolve" + +checkdim2 n m + | rows m == n && cols m == n = m + | otherwise = error $ "Error: "++ show n ++ "x" ++ show n + ++ " Jacobian expected in odeSolve" + +checkTimes ts | dim ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts + | otherwise = error "odeSolve requires increasing times" + where ts' = toList ts diff --git a/packages/hmatrix/src/Numeric/GSL/Polynomials.hs b/packages/hmatrix/src/Numeric/GSL/Polynomials.hs new file mode 100644 index 0000000..290c615 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Polynomials.hs @@ -0,0 +1,58 @@ +{-# LANGUAGE CPP, ForeignFunctionInterface #-} +----------------------------------------------------------------------------- +{- | +Module : Numeric.GSL.Polynomials +Copyright : (c) Alberto Ruiz 2006 +License : GPL-style + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional +Portability : uses ffi + +Polynomials. + + + +-} +----------------------------------------------------------------------------- +module Numeric.GSL.Polynomials ( + polySolve +) where + +import Data.Packed.Internal +import Data.Complex +import System.IO.Unsafe (unsafePerformIO) + +#if __GLASGOW_HASKELL__ >= 704 +import Foreign.C.Types (CInt(..)) +#endif + +{- | Solution of general polynomial equations, using /gsl_poly_complex_solve/. + +For example, the three solutions of x^3 + 8 = 0 + +>>> polySolve [8,0,0,1] +[(-2.0) :+ 0.0,1.0 :+ 1.7320508075688776,1.0 :+ (-1.7320508075688776)] + + +The example in the GSL manual: To find the roots of x^5 -1 = 0: + +>>> polySolve [-1, 0, 0, 0, 0, 1] +[(-0.8090169943749472) :+ 0.5877852522924731, +(-0.8090169943749472) :+ (-0.5877852522924731), +0.30901699437494756 :+ 0.9510565162951535, +0.30901699437494756 :+ (-0.9510565162951535), +1.0000000000000002 :+ 0.0] + +-} +polySolve :: [Double] -> [Complex Double] +polySolve = toList . polySolve' . fromList + +polySolve' :: Vector Double -> Vector (Complex Double) +polySolve' v | dim v > 1 = unsafePerformIO $ do + r <- createVector (dim v-1) + app2 c_polySolve vec v vec r "polySolve" + return r + | otherwise = error "polySolve on a polynomial of degree zero" + +foreign import ccall unsafe "gsl-aux.h polySolve" c_polySolve:: TVCV diff --git a/packages/hmatrix/src/Numeric/GSL/Root.hs b/packages/hmatrix/src/Numeric/GSL/Root.hs new file mode 100644 index 0000000..9d561c4 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Root.hs @@ -0,0 +1,199 @@ +{- | +Module : Numeric.GSL.Root +Copyright : (c) Alberto Ruiz 2009 +License : GPL + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional +Portability : uses ffi + +Multidimensional root finding. + + + +The example in the GSL manual: + +>>> let rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ] +>>> let (sol,path) = root Hybrids 1E-7 30 (rosenbrock 1 10) [-10,-5] +>>> sol +[1.0,1.0] +>>> disp 3 path +11x5 + 1.000 -10.000 -5.000 11.000 -1050.000 + 2.000 -3.976 24.827 4.976 90.203 + 3.000 -3.976 24.827 4.976 90.203 + 4.000 -3.976 24.827 4.976 90.203 + 5.000 -1.274 -5.680 2.274 -73.018 + 6.000 -1.274 -5.680 2.274 -73.018 + 7.000 0.249 0.298 0.751 2.359 + 8.000 0.249 0.298 0.751 2.359 + 9.000 1.000 0.878 -0.000 -1.218 +10.000 1.000 0.989 -0.000 -0.108 +11.000 1.000 1.000 0.000 0.000 + +-} +----------------------------------------------------------------------------- + +module Numeric.GSL.Root ( + uniRoot, UniRootMethod(..), + uniRootJ, UniRootMethodJ(..), + root, RootMethod(..), + rootJ, RootMethodJ(..), +) where + +import Data.Packed.Internal +import Data.Packed.Matrix +import Numeric.GSL.Internal +import Foreign.Ptr(FunPtr, freeHaskellFunPtr) +import Foreign.C.Types +import System.IO.Unsafe(unsafePerformIO) + +------------------------------------------------------------------------- + +data UniRootMethod = Bisection + | FalsePos + | Brent + deriving (Enum, Eq, Show, Bounded) + +uniRoot :: UniRootMethod + -> Double + -> Int + -> (Double -> Double) + -> Double + -> Double + -> (Double, Matrix Double) +uniRoot method epsrel maxit fun xl xu = uniRootGen (fi (fromEnum method)) fun xl xu epsrel maxit + +uniRootGen m f xl xu epsrel maxit = unsafePerformIO $ do + fp <- mkDoublefun f + rawpath <- createMIO maxit 4 + (c_root m fp epsrel (fi maxit) xl xu) + "root" + let it = round (rawpath @@> (maxit-1,0)) + path = takeRows it rawpath + [sol] = toLists $ dropRows (it-1) path + freeHaskellFunPtr fp + return (sol !! 1, path) + +foreign import ccall safe "root" + c_root:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> TM + +------------------------------------------------------------------------- +data UniRootMethodJ = UNewton + | Secant + | Steffenson + deriving (Enum, Eq, Show, Bounded) + +uniRootJ :: UniRootMethodJ + -> Double + -> Int + -> (Double -> Double) + -> (Double -> Double) + -> Double + -> (Double, Matrix Double) +uniRootJ method epsrel maxit fun dfun x = uniRootJGen (fi (fromEnum method)) fun + dfun x epsrel maxit + +uniRootJGen m f df x epsrel maxit = unsafePerformIO $ do + fp <- mkDoublefun f + dfp <- mkDoublefun df + rawpath <- createMIO maxit 2 + (c_rootj m fp dfp epsrel (fi maxit) x) + "rootj" + let it = round (rawpath @@> (maxit-1,0)) + path = takeRows it rawpath + [sol] = toLists $ dropRows (it-1) path + freeHaskellFunPtr fp + return (sol !! 1, path) + +foreign import ccall safe "rootj" + c_rootj :: CInt -> FunPtr (Double -> Double) -> FunPtr (Double -> Double) + -> Double -> CInt -> Double -> TM + +------------------------------------------------------------------------- + +data RootMethod = Hybrids + | Hybrid + | DNewton + | Broyden + deriving (Enum,Eq,Show,Bounded) + +-- | Nonlinear multidimensional root finding using algorithms that do not require +-- any derivative information to be supplied by the user. +-- Any derivatives needed are approximated by finite differences. +root :: RootMethod + -> Double -- ^ maximum residual + -> Int -- ^ maximum number of iterations allowed + -> ([Double] -> [Double]) -- ^ function to minimize + -> [Double] -- ^ starting point + -> ([Double], Matrix Double) -- ^ solution vector and optimization path + +root method epsabs maxit fun xinit = rootGen (fi (fromEnum method)) fun xinit epsabs maxit + +rootGen m f xi epsabs maxit = unsafePerformIO $ do + let xiv = fromList xi + n = dim xiv + fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList)) + rawpath <- vec xiv $ \xiv' -> + createMIO maxit (2*n+1) + (c_multiroot m fp epsabs (fi maxit) // xiv') + "multiroot" + let it = round (rawpath @@> (maxit-1,0)) + path = takeRows it rawpath + [sol] = toLists $ dropRows (it-1) path + freeHaskellFunPtr fp + return (take n $ drop 1 sol, path) + + +foreign import ccall safe "multiroot" + c_multiroot:: CInt -> FunPtr TVV -> Double -> CInt -> TVM + +------------------------------------------------------------------------- + +data RootMethodJ = HybridsJ + | HybridJ + | Newton + | GNewton + deriving (Enum,Eq,Show,Bounded) + +-- | Nonlinear multidimensional root finding using both the function and its derivatives. +rootJ :: RootMethodJ + -> Double -- ^ maximum residual + -> Int -- ^ maximum number of iterations allowed + -> ([Double] -> [Double]) -- ^ function to minimize + -> ([Double] -> [[Double]]) -- ^ Jacobian + -> [Double] -- ^ starting point + -> ([Double], Matrix Double) -- ^ solution vector and optimization path + +rootJ method epsabs maxit fun jac xinit = rootJGen (fi (fromEnum method)) fun jac xinit epsabs maxit + +rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do + let xiv = fromList xi + n = dim xiv + fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList)) + jp <- mkVecMatfun (aux_vTom (checkdim2 n . fromLists . jac . toList)) + rawpath <- vec xiv $ \xiv' -> + createMIO maxit (2*n+1) + (c_multirootj m fp jp epsabs (fi maxit) // xiv') + "multiroot" + let it = round (rawpath @@> (maxit-1,0)) + path = takeRows it rawpath + [sol] = toLists $ dropRows (it-1) path + freeHaskellFunPtr fp + freeHaskellFunPtr jp + return (take n $ drop 1 sol, path) + +foreign import ccall safe "multirootj" + c_multirootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM + +------------------------------------------------------- + +checkdim1 n v + | dim v == n = v + | otherwise = error $ "Error: "++ show n + ++ " components expected in the result of the function supplied to root" + +checkdim2 n m + | rows m == n && cols m == n = m + | otherwise = error $ "Error: "++ show n ++ "x" ++ show n + ++ " Jacobian expected in rootJ" diff --git a/packages/hmatrix/src/Numeric/GSL/Vector.hs b/packages/hmatrix/src/Numeric/GSL/Vector.hs new file mode 100644 index 0000000..6204b8e --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Vector.hs @@ -0,0 +1,328 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.GSL.Vector +-- Copyright : (c) Alberto Ruiz 2007 +-- License : GPL-style +-- +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- Portability : portable (uses FFI) +-- +-- Low level interface to vector operations. +-- +----------------------------------------------------------------------------- + +module Numeric.GSL.Vector ( + sumF, sumR, sumQ, sumC, + prodF, prodR, prodQ, prodC, + dotF, dotR, dotQ, dotC, + FunCodeS(..), toScalarR, toScalarF, toScalarC, toScalarQ, + FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ, + FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF, vectorMapValQ, + FunCodeVV(..), vectorZipR, vectorZipC, vectorZipF, vectorZipQ, + RandDist(..), randomVector +) where + +import Data.Packed.Internal.Common +import Data.Packed.Internal.Signatures +import Data.Packed.Internal.Vector + +import Data.Complex +import Foreign.Marshal.Alloc(free) +import Foreign.Marshal.Array(newArray) +import Foreign.Ptr(Ptr) +import Foreign.C.Types +import System.IO.Unsafe(unsafePerformIO) +import Control.Monad(when) + +fromei x = fromIntegral (fromEnum x) :: CInt + +data FunCodeV = Sin + | Cos + | Tan + | Abs + | ASin + | ACos + | ATan + | Sinh + | Cosh + | Tanh + | ASinh + | ACosh + | ATanh + | Exp + | Log + | Sign + | Sqrt + deriving Enum + +data FunCodeSV = Scale + | Recip + | AddConstant + | Negate + | PowSV + | PowVS + deriving Enum + +data FunCodeVV = Add + | Sub + | Mul + | Div + | Pow + | ATan2 + deriving Enum + +data FunCodeS = Norm2 + | AbsSum + | MaxIdx + | Max + | MinIdx + | Min + deriving Enum + +------------------------------------------------------------------ + +-- | sum of elements +sumF :: Vector Float -> Float +sumF x = unsafePerformIO $ do + r <- createVector 1 + app2 c_sumF vec x vec r "sumF" + return $ r @> 0 + +-- | sum of elements +sumR :: Vector Double -> Double +sumR x = unsafePerformIO $ do + r <- createVector 1 + app2 c_sumR vec x vec r "sumR" + return $ r @> 0 + +-- | sum of elements +sumQ :: Vector (Complex Float) -> Complex Float +sumQ x = unsafePerformIO $ do + r <- createVector 1 + app2 c_sumQ vec x vec r "sumQ" + return $ r @> 0 + +-- | sum of elements +sumC :: Vector (Complex Double) -> Complex Double +sumC x = unsafePerformIO $ do + r <- createVector 1 + app2 c_sumC vec x vec r "sumC" + return $ r @> 0 + +foreign import ccall unsafe "gsl-aux.h sumF" c_sumF :: TFF +foreign import ccall unsafe "gsl-aux.h sumR" c_sumR :: TVV +foreign import ccall unsafe "gsl-aux.h sumQ" c_sumQ :: TQVQV +foreign import ccall unsafe "gsl-aux.h sumC" c_sumC :: TCVCV + +-- | product of elements +prodF :: Vector Float -> Float +prodF x = unsafePerformIO $ do + r <- createVector 1 + app2 c_prodF vec x vec r "prodF" + return $ r @> 0 + +-- | product of elements +prodR :: Vector Double -> Double +prodR x = unsafePerformIO $ do + r <- createVector 1 + app2 c_prodR vec x vec r "prodR" + return $ r @> 0 + +-- | product of elements +prodQ :: Vector (Complex Float) -> Complex Float +prodQ x = unsafePerformIO $ do + r <- createVector 1 + app2 c_prodQ vec x vec r "prodQ" + return $ r @> 0 + +-- | product of elements +prodC :: Vector (Complex Double) -> Complex Double +prodC x = unsafePerformIO $ do + r <- createVector 1 + app2 c_prodC vec x vec r "prodC" + return $ r @> 0 + +foreign import ccall unsafe "gsl-aux.h prodF" c_prodF :: TFF +foreign import ccall unsafe "gsl-aux.h prodR" c_prodR :: TVV +foreign import ccall unsafe "gsl-aux.h prodQ" c_prodQ :: TQVQV +foreign import ccall unsafe "gsl-aux.h prodC" c_prodC :: TCVCV + +-- | dot product +dotF :: Vector Float -> Vector Float -> Float +dotF x y = unsafePerformIO $ do + r <- createVector 1 + app3 c_dotF vec x vec y vec r "dotF" + return $ r @> 0 + +-- | dot product +dotR :: Vector Double -> Vector Double -> Double +dotR x y = unsafePerformIO $ do + r <- createVector 1 + app3 c_dotR vec x vec y vec r "dotR" + return $ r @> 0 + +-- | dot product +dotQ :: Vector (Complex Float) -> Vector (Complex Float) -> Complex Float +dotQ x y = unsafePerformIO $ do + r <- createVector 1 + app3 c_dotQ vec x vec y vec r "dotQ" + return $ r @> 0 + +-- | dot product +dotC :: Vector (Complex Double) -> Vector (Complex Double) -> Complex Double +dotC x y = unsafePerformIO $ do + r <- createVector 1 + app3 c_dotC vec x vec y vec r "dotC" + return $ r @> 0 + +foreign import ccall unsafe "gsl-aux.h dotF" c_dotF :: TFFF +foreign import ccall unsafe "gsl-aux.h dotR" c_dotR :: TVVV +foreign import ccall unsafe "gsl-aux.h dotQ" c_dotQ :: TQVQVQV +foreign import ccall unsafe "gsl-aux.h dotC" c_dotC :: TCVCVCV + +------------------------------------------------------------------ + +toScalarAux fun code v = unsafePerformIO $ do + r <- createVector 1 + app2 (fun (fromei code)) vec v vec r "toScalarAux" + return (r `at` 0) + +vectorMapAux fun code v = unsafePerformIO $ do + r <- createVector (dim v) + app2 (fun (fromei code)) vec v vec r "vectorMapAux" + return r + +vectorMapValAux fun code val v = unsafePerformIO $ do + r <- createVector (dim v) + pval <- newArray [val] + app2 (fun (fromei code) pval) vec v vec r "vectorMapValAux" + free pval + return r + +vectorZipAux fun code u v = unsafePerformIO $ do + r <- createVector (dim u) + when (dim u > 0) $ app3 (fun (fromei code)) vec u vec v vec r "vectorZipAux" + return r + +--------------------------------------------------------------------- + +-- | obtains different functions of a vector: norm1, norm2, max, min, posmax, posmin, etc. +toScalarR :: FunCodeS -> Vector Double -> Double +toScalarR oper = toScalarAux c_toScalarR (fromei oper) + +foreign import ccall unsafe "gsl-aux.h toScalarR" c_toScalarR :: CInt -> TVV + +-- | obtains different functions of a vector: norm1, norm2, max, min, posmax, posmin, etc. +toScalarF :: FunCodeS -> Vector Float -> Float +toScalarF oper = toScalarAux c_toScalarF (fromei oper) + +foreign import ccall unsafe "gsl-aux.h toScalarF" c_toScalarF :: CInt -> TFF + +-- | obtains different functions of a vector: only norm1, norm2 +toScalarC :: FunCodeS -> Vector (Complex Double) -> Double +toScalarC oper = toScalarAux c_toScalarC (fromei oper) + +foreign import ccall unsafe "gsl-aux.h toScalarC" c_toScalarC :: CInt -> TCVV + +-- | obtains different functions of a vector: only norm1, norm2 +toScalarQ :: FunCodeS -> Vector (Complex Float) -> Float +toScalarQ oper = toScalarAux c_toScalarQ (fromei oper) + +foreign import ccall unsafe "gsl-aux.h toScalarQ" c_toScalarQ :: CInt -> TQVF + +------------------------------------------------------------------ + +-- | map of real vectors with given function +vectorMapR :: FunCodeV -> Vector Double -> Vector Double +vectorMapR = vectorMapAux c_vectorMapR + +foreign import ccall unsafe "gsl-aux.h mapR" c_vectorMapR :: CInt -> TVV + +-- | map of complex vectors with given function +vectorMapC :: FunCodeV -> Vector (Complex Double) -> Vector (Complex Double) +vectorMapC oper = vectorMapAux c_vectorMapC (fromei oper) + +foreign import ccall unsafe "gsl-aux.h mapC" c_vectorMapC :: CInt -> TCVCV + +-- | map of real vectors with given function +vectorMapF :: FunCodeV -> Vector Float -> Vector Float +vectorMapF = vectorMapAux c_vectorMapF + +foreign import ccall unsafe "gsl-aux.h mapF" c_vectorMapF :: CInt -> TFF + +-- | map of real vectors with given function +vectorMapQ :: FunCodeV -> Vector (Complex Float) -> Vector (Complex Float) +vectorMapQ = vectorMapAux c_vectorMapQ + +foreign import ccall unsafe "gsl-aux.h mapQ" c_vectorMapQ :: CInt -> TQVQV + +------------------------------------------------------------------- + +-- | map of real vectors with given function +vectorMapValR :: FunCodeSV -> Double -> Vector Double -> Vector Double +vectorMapValR oper = vectorMapValAux c_vectorMapValR (fromei oper) + +foreign import ccall unsafe "gsl-aux.h mapValR" c_vectorMapValR :: CInt -> Ptr Double -> TVV + +-- | map of complex vectors with given function +vectorMapValC :: FunCodeSV -> Complex Double -> Vector (Complex Double) -> Vector (Complex Double) +vectorMapValC = vectorMapValAux c_vectorMapValC + +foreign import ccall unsafe "gsl-aux.h mapValC" c_vectorMapValC :: CInt -> Ptr (Complex Double) -> TCVCV + +-- | map of real vectors with given function +vectorMapValF :: FunCodeSV -> Float -> Vector Float -> Vector Float +vectorMapValF oper = vectorMapValAux c_vectorMapValF (fromei oper) + +foreign import ccall unsafe "gsl-aux.h mapValF" c_vectorMapValF :: CInt -> Ptr Float -> TFF + +-- | map of complex vectors with given function +vectorMapValQ :: FunCodeSV -> Complex Float -> Vector (Complex Float) -> Vector (Complex Float) +vectorMapValQ oper = vectorMapValAux c_vectorMapValQ (fromei oper) + +foreign import ccall unsafe "gsl-aux.h mapValQ" c_vectorMapValQ :: CInt -> Ptr (Complex Float) -> TQVQV + +------------------------------------------------------------------- + +-- | elementwise operation on real vectors +vectorZipR :: FunCodeVV -> Vector Double -> Vector Double -> Vector Double +vectorZipR = vectorZipAux c_vectorZipR + +foreign import ccall unsafe "gsl-aux.h zipR" c_vectorZipR :: CInt -> TVVV + +-- | elementwise operation on complex vectors +vectorZipC :: FunCodeVV -> Vector (Complex Double) -> Vector (Complex Double) -> Vector (Complex Double) +vectorZipC = vectorZipAux c_vectorZipC + +foreign import ccall unsafe "gsl-aux.h zipC" c_vectorZipC :: CInt -> TCVCVCV + +-- | elementwise operation on real vectors +vectorZipF :: FunCodeVV -> Vector Float -> Vector Float -> Vector Float +vectorZipF = vectorZipAux c_vectorZipF + +foreign import ccall unsafe "gsl-aux.h zipF" c_vectorZipF :: CInt -> TFFF + +-- | elementwise operation on complex vectors +vectorZipQ :: FunCodeVV -> Vector (Complex Float) -> Vector (Complex Float) -> Vector (Complex Float) +vectorZipQ = vectorZipAux c_vectorZipQ + +foreign import ccall unsafe "gsl-aux.h zipQ" c_vectorZipQ :: CInt -> TQVQVQV + +----------------------------------------------------------------------- + +data RandDist = Uniform -- ^ uniform distribution in [0,1) + | Gaussian -- ^ normal distribution with mean zero and standard deviation one + deriving Enum + +-- | Obtains a vector of pseudorandom elements from the the mt19937 generator in GSL, with a given seed. Use randomIO to get a random seed. +randomVector :: Int -- ^ seed + -> RandDist -- ^ distribution + -> Int -- ^ vector size + -> Vector Double +randomVector seed dist n = unsafePerformIO $ do + r <- createVector n + app1 (c_random_vector (fi seed) ((fi.fromEnum) dist)) vec r "randomVector" + return r + +foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV diff --git a/packages/hmatrix/src/Numeric/GSL/gsl-aux.c b/packages/hmatrix/src/Numeric/GSL/gsl-aux.c new file mode 100644 index 0000000..410d157 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/gsl-aux.c @@ -0,0 +1,1541 @@ +#include + +#define RVEC(A) int A##n, double*A##p +#define RMAT(A) int A##r, int A##c, double* A##p +#define KRVEC(A) int A##n, const double*A##p +#define KRMAT(A) int A##r, int A##c, const double* A##p + +#define CVEC(A) int A##n, gsl_complex*A##p +#define CMAT(A) int A##r, int A##c, gsl_complex* A##p +#define KCVEC(A) int A##n, const gsl_complex*A##p +#define KCMAT(A) int A##r, int A##c, const gsl_complex* A##p + +#define FVEC(A) int A##n, float*A##p +#define FMAT(A) int A##r, int A##c, float* A##p +#define KFVEC(A) int A##n, const float*A##p +#define KFMAT(A) int A##r, int A##c, const float* A##p + +#define QVEC(A) int A##n, gsl_complex_float*A##p +#define QMAT(A) int A##r, int A##c, gsl_complex_float* A##p +#define KQVEC(A) int A##n, const gsl_complex_float*A##p +#define KQMAT(A) int A##r, int A##c, const gsl_complex_float* A##p + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#define MACRO(B) do {B} while (0) +#define ERROR(CODE) MACRO(return CODE;) +#define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);}) +#define OK return 0; + +#define MIN(A,B) ((A)<(B)?(A):(B)) +#define MAX(A,B) ((A)>(B)?(A):(B)) + +#ifdef DBG +#define DEBUGMSG(M) printf("*** calling aux C function: %s\n",M); +#else +#define DEBUGMSG(M) +#endif + +#define CHECK(RES,CODE) MACRO(if(RES) return CODE;) + +#ifdef DBG +#define DEBUGMAT(MSG,X) printf(MSG" = \n"); gsl_matrix_fprintf(stdout,X,"%f"); printf("\n"); +#else +#define DEBUGMAT(MSG,X) +#endif + +#ifdef DBG +#define DEBUGVEC(MSG,X) printf(MSG" = \n"); gsl_vector_fprintf(stdout,X,"%f"); printf("\n"); +#else +#define DEBUGVEC(MSG,X) +#endif + +#define DVVIEW(A) gsl_vector_view A = gsl_vector_view_array(A##p,A##n) +#define DMVIEW(A) gsl_matrix_view A = gsl_matrix_view_array(A##p,A##r,A##c) +#define CVVIEW(A) gsl_vector_complex_view A = gsl_vector_complex_view_array((double*)A##p,A##n) +#define CMVIEW(A) gsl_matrix_complex_view A = gsl_matrix_complex_view_array((double*)A##p,A##r,A##c) +#define KDVVIEW(A) gsl_vector_const_view A = gsl_vector_const_view_array(A##p,A##n) +#define KDMVIEW(A) gsl_matrix_const_view A = gsl_matrix_const_view_array(A##p,A##r,A##c) +#define KCVVIEW(A) gsl_vector_complex_const_view A = gsl_vector_complex_const_view_array((double*)A##p,A##n) +#define KCMVIEW(A) gsl_matrix_complex_const_view A = gsl_matrix_complex_const_view_array((double*)A##p,A##r,A##c) + +#define FVVIEW(A) gsl_vector_float_view A = gsl_vector_float_view_array(A##p,A##n) +#define FMVIEW(A) gsl_matrix_float_view A = gsl_matrix_float_view_array(A##p,A##r,A##c) +#define QVVIEW(A) gsl_vector_complex_float_view A = gsl_vector_float_complex_view_array((float*)A##p,A##n) +#define QMVIEW(A) gsl_matrix_complex_float_view A = gsl_matrix_float_complex_view_array((float*)A##p,A##r,A##c) +#define KFVVIEW(A) gsl_vector_float_const_view A = gsl_vector_float_const_view_array(A##p,A##n) +#define KFMVIEW(A) gsl_matrix_float_const_view A = gsl_matrix_float_const_view_array(A##p,A##r,A##c) +#define KQVVIEW(A) gsl_vector_complex_float_const_view A = gsl_vector_complex_float_const_view_array((float*)A##p,A##n) +#define KQMVIEW(A) gsl_matrix_complex_float_const_view A = gsl_matrix_complex_float_const_view_array((float*)A##p,A##r,A##c) + +#define V(a) (&a.vector) +#define M(a) (&a.matrix) + +#define GCVEC(A) int A##n, gsl_complex*A##p +#define KGCVEC(A) int A##n, const gsl_complex*A##p + +#define GQVEC(A) int A##n, gsl_complex_float*A##p +#define KGQVEC(A) int A##n, const gsl_complex_float*A##p + +#define BAD_SIZE 2000 +#define BAD_CODE 2001 +#define MEM 2002 +#define BAD_FILE 2003 + + +void no_abort_on_error() { + gsl_set_error_handler_off(); +} + + +int sumF(KFVEC(x),FVEC(r)) { + DEBUGMSG("sumF"); + REQUIRES(rn==1,BAD_SIZE); + int i; + float res = 0; + for (i = 0; i < xn; i++) res += xp[i]; + rp[0] = res; + OK +} + +int sumR(KRVEC(x),RVEC(r)) { + DEBUGMSG("sumR"); + REQUIRES(rn==1,BAD_SIZE); + int i; + double res = 0; + for (i = 0; i < xn; i++) res += xp[i]; + rp[0] = res; + OK +} + +int sumQ(KQVEC(x),QVEC(r)) { + DEBUGMSG("sumQ"); + REQUIRES(rn==1,BAD_SIZE); + int i; + gsl_complex_float res; + res.dat[0] = 0; + res.dat[1] = 0; + for (i = 0; i < xn; i++) { + res.dat[0] += xp[i].dat[0]; + res.dat[1] += xp[i].dat[1]; + } + rp[0] = res; + OK +} + +int sumC(KCVEC(x),CVEC(r)) { + DEBUGMSG("sumC"); + REQUIRES(rn==1,BAD_SIZE); + int i; + gsl_complex res; + res.dat[0] = 0; + res.dat[1] = 0; + for (i = 0; i < xn; i++) { + res.dat[0] += xp[i].dat[0]; + res.dat[1] += xp[i].dat[1]; + } + rp[0] = res; + OK +} + +int prodF(KFVEC(x),FVEC(r)) { + DEBUGMSG("prodF"); + REQUIRES(rn==1,BAD_SIZE); + int i; + float res = 1; + for (i = 0; i < xn; i++) res *= xp[i]; + rp[0] = res; + OK +} + +int prodR(KRVEC(x),RVEC(r)) { + DEBUGMSG("prodR"); + REQUIRES(rn==1,BAD_SIZE); + int i; + double res = 1; + for (i = 0; i < xn; i++) res *= xp[i]; + rp[0] = res; + OK +} + +int prodQ(KQVEC(x),QVEC(r)) { + DEBUGMSG("prodQ"); + REQUIRES(rn==1,BAD_SIZE); + int i; + gsl_complex_float res; + float temp; + res.dat[0] = 1; + res.dat[1] = 0; + for (i = 0; i < xn; i++) { + temp = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1]; + res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0]; + res.dat[0] = temp; + } + rp[0] = res; + OK +} + +int prodC(KCVEC(x),CVEC(r)) { + DEBUGMSG("prodC"); + REQUIRES(rn==1,BAD_SIZE); + int i; + gsl_complex res; + double temp; + res.dat[0] = 1; + res.dat[1] = 0; + for (i = 0; i < xn; i++) { + temp = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1]; + res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0]; + res.dat[0] = temp; + } + rp[0] = res; + OK +} + +int dotF(KFVEC(x), KFVEC(y), FVEC(r)) { + DEBUGMSG("dotF"); + REQUIRES(xn==yn,BAD_SIZE); + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("dotF"); + KFVVIEW(x); + KFVVIEW(y); + gsl_blas_sdot(V(x),V(y),rp); + OK +} + +int dotR(KRVEC(x), KRVEC(y), RVEC(r)) { + DEBUGMSG("dotR"); + REQUIRES(xn==yn,BAD_SIZE); + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("dotR"); + KDVVIEW(x); + KDVVIEW(y); + gsl_blas_ddot(V(x),V(y),rp); + OK +} + +int dotQ(KQVEC(x), KQVEC(y), QVEC(r)) { + DEBUGMSG("dotQ"); + REQUIRES(xn==yn,BAD_SIZE); + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("dotQ"); + KQVVIEW(x); + KQVVIEW(y); + gsl_blas_cdotu(V(x),V(y),rp); + OK +} + +int dotC(KCVEC(x), KCVEC(y), CVEC(r)) { + DEBUGMSG("dotC"); + REQUIRES(xn==yn,BAD_SIZE); + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("dotC"); + KCVVIEW(x); + KCVVIEW(y); + gsl_blas_zdotu(V(x),V(y),rp); + OK +} + +int toScalarR(int code, KRVEC(x), RVEC(r)) { + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("toScalarR"); + KDVVIEW(x); + double res; + switch(code) { + case 0: { res = gsl_blas_dnrm2(V(x)); break; } + case 1: { res = gsl_blas_dasum(V(x)); break; } + case 2: { res = gsl_vector_max_index(V(x)); break; } + case 3: { res = gsl_vector_max(V(x)); break; } + case 4: { res = gsl_vector_min_index(V(x)); break; } + case 5: { res = gsl_vector_min(V(x)); break; } + default: ERROR(BAD_CODE); + } + rp[0] = res; + OK +} + +int toScalarF(int code, KFVEC(x), FVEC(r)) { + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("toScalarF"); + KFVVIEW(x); + float res; + switch(code) { + case 0: { res = gsl_blas_snrm2(V(x)); break; } + case 1: { res = gsl_blas_sasum(V(x)); break; } + case 2: { res = gsl_vector_float_max_index(V(x)); break; } + case 3: { res = gsl_vector_float_max(V(x)); break; } + case 4: { res = gsl_vector_float_min_index(V(x)); break; } + case 5: { res = gsl_vector_float_min(V(x)); break; } + default: ERROR(BAD_CODE); + } + rp[0] = res; + OK +} + + +int toScalarC(int code, KCVEC(x), RVEC(r)) { + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("toScalarC"); + KCVVIEW(x); + double res; + switch(code) { + case 0: { res = gsl_blas_dznrm2(V(x)); break; } + case 1: { res = gsl_blas_dzasum(V(x)); break; } + default: ERROR(BAD_CODE); + } + rp[0] = res; + OK +} + +int toScalarQ(int code, KQVEC(x), FVEC(r)) { + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("toScalarQ"); + KQVVIEW(x); + float res; + switch(code) { + case 0: { res = gsl_blas_scnrm2(V(x)); break; } + case 1: { res = gsl_blas_scasum(V(x)); break; } + default: ERROR(BAD_CODE); + } + rp[0] = res; + OK +} + + +inline double sign(double x) { + if(x>0) { + return +1.0; + } else if (x<0) { + return -1.0; + } else { + return 0.0; + } +} + +inline float float_sign(float x) { + if(x>0) { + return +1.0; + } else if (x<0) { + return -1.0; + } else { + return 0.0; + } +} + +inline gsl_complex complex_abs(gsl_complex z) { + gsl_complex r; + r.dat[0] = gsl_complex_abs(z); + r.dat[1] = 0; + return r; +} + +inline gsl_complex complex_signum(gsl_complex z) { + gsl_complex r; + double mag; + if (z.dat[0] == 0 && z.dat[1] == 0) { + r.dat[0] = 0; + r.dat[1] = 0; + } else { + mag = gsl_complex_abs(z); + r.dat[0] = z.dat[0]/mag; + r.dat[1] = z.dat[1]/mag; + } + return r; +} + +#define OP(C,F) case C: { for(k=0;k1,BAD_SIZE); + gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc (an); + int res = gsl_poly_complex_solve ((double*)ap, an, w, (double*)zp); + CHECK(res,res); + gsl_poly_complex_workspace_free (w); + OK; +} + +int vector_fscanf(char*filename, RVEC(a)) { + DEBUGMSG("gsl_vector_fscanf"); + DVVIEW(a); + FILE * f = fopen(filename,"r"); + CHECK(!f,BAD_FILE); + int res = gsl_vector_fscanf(f,V(a)); + CHECK(res,res); + fclose (f); + OK +} + +int vector_fprintf(char*filename, char*fmt, RVEC(a)) { + DEBUGMSG("gsl_vector_fprintf"); + DVVIEW(a); + FILE * f = fopen(filename,"w"); + CHECK(!f,BAD_FILE); + int res = gsl_vector_fprintf(f,V(a),fmt); + CHECK(res,res); + fclose (f); + OK +} + +int vector_fread(char*filename, RVEC(a)) { + DEBUGMSG("gsl_vector_fread"); + DVVIEW(a); + FILE * f = fopen(filename,"r"); + CHECK(!f,BAD_FILE); + int res = gsl_vector_fread(f,V(a)); + CHECK(res,res); + fclose (f); + OK +} + +int vector_fwrite(char*filename, RVEC(a)) { + DEBUGMSG("gsl_vector_fwrite"); + DVVIEW(a); + FILE * f = fopen(filename,"w"); + CHECK(!f,BAD_FILE); + int res = gsl_vector_fwrite(f,V(a)); + CHECK(res,res); + fclose (f); + OK +} + +int matrix_fprintf(char*filename, char*fmt, int ro, RMAT(m)) { + DEBUGMSG("matrix_fprintf"); + FILE * f = fopen(filename,"w"); + CHECK(!f,BAD_FILE); + int i,j,sr,sc; + if (ro==1) { sr = mc; sc = 1;} else { sr = 1; sc = mr;} + #define AT(M,r,c) (M##p[(r)*sr+(c)*sc]) + for (i=0; isize,sizeof(double)); + int k; + for(k=0;ksize;k++) { + p[k] = gsl_vector_get(x,k); + } + double res = f(x->size,p); + free(p); + return res; +} + +double only_f_aux_root(double x, void *pars); +int uniMinimize(int method, double f(double), + double epsrel, int maxit, double min, + double xl, double xu, RMAT(sol)) { + REQUIRES(solr == maxit && solc == 4,BAD_SIZE); + DEBUGMSG("minimize_only_f"); + gsl_function my_func; + my_func.function = only_f_aux_root; + my_func.params = f; + size_t iter = 0; + int status; + const gsl_min_fminimizer_type *T; + gsl_min_fminimizer *s; + // Starting point + switch(method) { + case 0 : {T = gsl_min_fminimizer_goldensection; break; } + case 1 : {T = gsl_min_fminimizer_brent; break; } + case 2 : {T = gsl_min_fminimizer_quad_golden; break; } + default: ERROR(BAD_CODE); + } + s = gsl_min_fminimizer_alloc (T); + gsl_min_fminimizer_set (s, &my_func, min, xl, xu); + do { + double current_min, current_lo, current_hi; + status = gsl_min_fminimizer_iterate (s); + current_min = gsl_min_fminimizer_x_minimum (s); + current_lo = gsl_min_fminimizer_x_lower (s); + current_hi = gsl_min_fminimizer_x_upper (s); + solp[iter*solc] = iter + 1; + solp[iter*solc+1] = current_min; + solp[iter*solc+2] = current_lo; + solp[iter*solc+3] = current_hi; + iter++; + if (status) /* check if solver is stuck */ + break; + + status = + gsl_min_test_interval (current_lo, current_hi, 0, epsrel); + } + while (status == GSL_CONTINUE && iter < maxit); + int i; + for (i=iter; ifval; + solp[iter*solc+2] = size; + + int k; + for(k=0;kx,k); + } + iter++; + if (status) break; + status = gsl_multimin_test_size (size, tolsize); + } while (status == GSL_CONTINUE && iter < maxit); + int i,j; + for (i=iter; isize,sizeof(double)); + int k; + for(k=0;ksize;k++) { + p[k] = gsl_vector_get(x,k); + } + double res = fdf->f(x->size,p); + free(p); + return res; +} + + +void df_aux_min(const gsl_vector * x, void * pars, gsl_vector * g) { + Tfdf * fdf = ((Tfdf*) pars); + double* p = (double*)calloc(x->size,sizeof(double)); + double* q = (double*)calloc(g->size,sizeof(double)); + int k; + for(k=0;ksize;k++) { + p[k] = gsl_vector_get(x,k); + } + + fdf->df(x->size,p,g->size,q); + + for(k=0;ksize;k++) { + gsl_vector_set(g,k,q[k]); + } + free(p); + free(q); +} + +void fdf_aux_min(const gsl_vector * x, void * pars, double * f, gsl_vector * g) { + *f = f_aux_min(x,pars); + df_aux_min(x,pars,g); +} + + +int minimizeD(int method, double f(int, double*), int df(int, double*, int, double*), + double initstep, double minimpar, double tolgrad, int maxit, + KRVEC(xi), RMAT(sol)) { + REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE); + DEBUGMSG("minimizeWithDeriv (conjugate_fr)"); + gsl_multimin_function_fdf my_func; + // extract function from pars + my_func.f = f_aux_min; + my_func.df = df_aux_min; + my_func.fdf = fdf_aux_min; + my_func.n = xin; + Tfdf stfdf; + stfdf.f = f; + stfdf.df = df; + my_func.params = &stfdf; + size_t iter = 0; + int status; + const gsl_multimin_fdfminimizer_type *T; + gsl_multimin_fdfminimizer *s = NULL; + // Starting point + KDVVIEW(xi); + // conjugate gradient fr + switch(method) { + case 0 : {T = gsl_multimin_fdfminimizer_conjugate_fr; break; } + case 1 : {T = gsl_multimin_fdfminimizer_conjugate_pr; break; } + case 2 : {T = gsl_multimin_fdfminimizer_vector_bfgs; break; } + case 3 : {T = gsl_multimin_fdfminimizer_vector_bfgs2; break; } + case 4 : {T = gsl_multimin_fdfminimizer_steepest_descent; break; } + default: ERROR(BAD_CODE); + } + s = gsl_multimin_fdfminimizer_alloc (T, my_func.n); + gsl_multimin_fdfminimizer_set (s, &my_func, V(xi), initstep, minimpar); + do { + status = gsl_multimin_fdfminimizer_iterate (s); + solp[iter*solc+0] = iter+1; + solp[iter*solc+1] = s->f; + int k; + for(k=0;kx,k); + } + iter++; + if (status) break; + status = gsl_multimin_test_gradient (s->gradient, tolgrad); + } while (status == GSL_CONTINUE && iter < maxit); + int i,j; + for (i=iter; if)(x); +} + +double jf_aux_uni(double x, void * pars) { + uniTfjf * fjf = ((uniTfjf*) pars); + return (fjf->jf)(x); +} + +void fjf_aux_uni(double x, void * pars, double * f, double * g) { + *f = f_aux_uni(x,pars); + *g = jf_aux_uni(x,pars); +} + +int rootj(int method, double f(double), + double df(double), + double epsrel, int maxit, + double x, RMAT(sol)) { + REQUIRES(solr == maxit && solc == 2,BAD_SIZE); + DEBUGMSG("root_fjf"); + gsl_function_fdf my_func; + // extract function from pars + my_func.f = f_aux_uni; + my_func.df = jf_aux_uni; + my_func.fdf = fjf_aux_uni; + uniTfjf stfjf; + stfjf.f = f; + stfjf.jf = df; + my_func.params = &stfjf; + size_t iter = 0; + int status; + const gsl_root_fdfsolver_type *T; + gsl_root_fdfsolver *s; + // Starting point + switch(method) { + case 0 : {T = gsl_root_fdfsolver_newton;; break; } + case 1 : {T = gsl_root_fdfsolver_secant; break; } + case 2 : {T = gsl_root_fdfsolver_steffenson; break; } + default: ERROR(BAD_CODE); + } + s = gsl_root_fdfsolver_alloc (T); + + gsl_root_fdfsolver_set (s, &my_func, x); + + do { + double x0; + status = gsl_root_fdfsolver_iterate (s); + x0 = x; + x = gsl_root_fdfsolver_root(s); + solp[iter*solc+0] = iter+1; + solp[iter*solc+1] = x; + + iter++; + if (status) /* check if solver is stuck */ + break; + + status = + gsl_root_test_delta (x, x0, 0, epsrel); + } + while (status == GSL_CONTINUE && iter < maxit); + + int i; + for (i=iter; isize,sizeof(double)); + double* q = (double*)calloc(y->size,sizeof(double)); + int k; + for(k=0;ksize;k++) { + p[k] = gsl_vector_get(x,k); + } + f(x->size,p,y->size,q); + for(k=0;ksize;k++) { + gsl_vector_set(y,k,q[k]); + } + free(p); + free(q); + return 0; //hmmm +} + +int multiroot(int method, void f(int, double*, int, double*), + double epsabs, int maxit, + KRVEC(xi), RMAT(sol)) { + REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE); + DEBUGMSG("root_only_f"); + gsl_multiroot_function my_func; + // extract function from pars + my_func.f = only_f_aux_multiroot; + my_func.n = xin; + my_func.params = f; + size_t iter = 0; + int status; + const gsl_multiroot_fsolver_type *T; + gsl_multiroot_fsolver *s; + // Starting point + KDVVIEW(xi); + switch(method) { + case 0 : {T = gsl_multiroot_fsolver_hybrids;; break; } + case 1 : {T = gsl_multiroot_fsolver_hybrid; break; } + case 2 : {T = gsl_multiroot_fsolver_dnewton; break; } + case 3 : {T = gsl_multiroot_fsolver_broyden; break; } + default: ERROR(BAD_CODE); + } + s = gsl_multiroot_fsolver_alloc (T, my_func.n); + gsl_multiroot_fsolver_set (s, &my_func, V(xi)); + + do { + status = gsl_multiroot_fsolver_iterate (s); + + solp[iter*solc+0] = iter+1; + + int k; + for(k=0;kx,k); + } + for(k=xin;k<2*xin;k++) { + solp[iter*solc+k+1] = gsl_vector_get(s->f,k-xin); + } + + iter++; + if (status) /* check if solver is stuck */ + break; + + status = + gsl_multiroot_test_residual (s->f, epsabs); + } + while (status == GSL_CONTINUE && iter < maxit); + + int i,j; + for (i=iter; isize,sizeof(double)); + double* q = (double*)calloc(y->size,sizeof(double)); + int k; + for(k=0;ksize;k++) { + p[k] = gsl_vector_get(x,k); + } + (fjf->f)(x->size,p,y->size,q); + for(k=0;ksize;k++) { + gsl_vector_set(y,k,q[k]); + } + free(p); + free(q); + return 0; +} + +int jf_aux(const gsl_vector * x, void * pars, gsl_matrix * jac) { + Tfjf * fjf = ((Tfjf*) pars); + double* p = (double*)calloc(x->size,sizeof(double)); + double* q = (double*)calloc((jac->size1)*(jac->size2),sizeof(double)); + int i,j,k; + for(k=0;ksize;k++) { + p[k] = gsl_vector_get(x,k); + } + + (fjf->jf)(x->size,p,jac->size1,jac->size2,q); + + k=0; + for(i=0;isize1;i++) { + for(j=0;jsize2;j++){ + gsl_matrix_set(jac,i,j,q[k++]); + } + } + free(p); + free(q); + return 0; +} + +int fjf_aux(const gsl_vector * x, void * pars, gsl_vector * f, gsl_matrix * g) { + f_aux(x,pars,f); + jf_aux(x,pars,g); + return 0; +} + +int multirootj(int method, int f(int, double*, int, double*), + int jac(int, double*, int, int, double*), + double epsabs, int maxit, + KRVEC(xi), RMAT(sol)) { + REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE); + DEBUGMSG("root_fjf"); + gsl_multiroot_function_fdf my_func; + // extract function from pars + my_func.f = f_aux; + my_func.df = jf_aux; + my_func.fdf = fjf_aux; + my_func.n = xin; + Tfjf stfjf; + stfjf.f = f; + stfjf.jf = jac; + my_func.params = &stfjf; + size_t iter = 0; + int status; + const gsl_multiroot_fdfsolver_type *T; + gsl_multiroot_fdfsolver *s; + // Starting point + KDVVIEW(xi); + switch(method) { + case 0 : {T = gsl_multiroot_fdfsolver_hybridsj;; break; } + case 1 : {T = gsl_multiroot_fdfsolver_hybridj; break; } + case 2 : {T = gsl_multiroot_fdfsolver_newton; break; } + case 3 : {T = gsl_multiroot_fdfsolver_gnewton; break; } + default: ERROR(BAD_CODE); + } + s = gsl_multiroot_fdfsolver_alloc (T, my_func.n); + + gsl_multiroot_fdfsolver_set (s, &my_func, V(xi)); + + do { + status = gsl_multiroot_fdfsolver_iterate (s); + + solp[iter*solc+0] = iter+1; + + int k; + for(k=0;kx,k); + } + for(k=xin;k<2*xin;k++) { + solp[iter*solc+k+1] = gsl_vector_get(s->f,k-xin); + } + + iter++; + if (status) /* check if solver is stuck */ + break; + + status = + gsl_multiroot_test_residual (s->f, epsabs); + } + while (status == GSL_CONTINUE && iter < maxit); + + int i,j; + for (i=iter; if); + + int k; + for(k=0;kx,k); + } + + iter++; + if (status) /* check if solver is stuck */ + break; + + status = gsl_multifit_test_delta (s->dx, s->x, epsabs, epsrel); + } + while (status == GSL_CONTINUE && iter < maxit); + + int i,j; + for (i=iter; iJ, 0.0, M(cov)); + + gsl_multifit_fdfsolver_free (s); + OK +} + + +////////////////////////////////////////////////////// + + +#define RAN(C,F) case C: { for(k=0;k + +typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode; + +int odefunc (double t, const double y[], double f[], void *params) { + Tode * P = (Tode*) params; + (P->f)(t,P->n,y,P->n,f); + return GSL_SUCCESS; +} + +int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) { + Tode * P = ((Tode*) params); + (P->j)(t,P->n,y,P->n,P->n,dfdy); + int j; + for (j=0; j< P->n; j++) + dfdt[j] = 0.0; + return GSL_SUCCESS; +} + + +int ode(int method, double h, double eps_abs, double eps_rel, + int f(double, int, const double*, int, double*), + int jac(double, int, const double*, int, int, double*), + KRVEC(xi), KRVEC(ts), RMAT(sol)) { + + const gsl_odeiv_step_type * T; + + switch(method) { + case 0 : {T = gsl_odeiv_step_rk2; break; } + case 1 : {T = gsl_odeiv_step_rk4; break; } + case 2 : {T = gsl_odeiv_step_rkf45; break; } + case 3 : {T = gsl_odeiv_step_rkck; break; } + case 4 : {T = gsl_odeiv_step_rk8pd; break; } + case 5 : {T = gsl_odeiv_step_rk2imp; break; } + case 6 : {T = gsl_odeiv_step_rk4imp; break; } + case 7 : {T = gsl_odeiv_step_bsimp; break; } + case 8 : { printf("Sorry: ODE rk1imp not available in this GSL version\n"); exit(0); } + case 9 : { printf("Sorry: ODE msadams not available in this GSL version\n"); exit(0); } + case 10: { printf("Sorry: ODE msbdf not available in this GSL version\n"); exit(0); } + default: ERROR(BAD_CODE); + } + + gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin); + gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel); + gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin); + + Tode P; + P.f = f; + P.j = jac; + P.n = xin; + + gsl_odeiv_system sys = {odefunc, odejac, xin, &P}; + + double t = tsp[0]; + + double* y = (double*)calloc(xin,sizeof(double)); + int i,j; + for(i=0; i< xin; i++) { + y[i] = xip[i]; + solp[i] = xip[i]; + } + + for (i = 1; i < tsn ; i++) + { + double ti = tsp[i]; + while (t < ti) + { + gsl_odeiv_evolve_apply (e, c, s, + &sys, + &t, ti, &h, + y); + // if (h < hmin) h = hmin; + } + for(j=0; j + +typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode; + +int odefunc (double t, const double y[], double f[], void *params) { + Tode * P = (Tode*) params; + (P->f)(t,P->n,y,P->n,f); + return GSL_SUCCESS; +} + +int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) { + Tode * P = ((Tode*) params); + (P->j)(t,P->n,y,P->n,P->n,dfdy); + int j; + for (j=0; j< P->n; j++) + dfdt[j] = 0.0; + return GSL_SUCCESS; +} + + +int ode(int method, double h, double eps_abs, double eps_rel, + int f(double, int, const double*, int, double*), + int jac(double, int, const double*, int, int, double*), + KRVEC(xi), KRVEC(ts), RMAT(sol)) { + + const gsl_odeiv2_step_type * T; + + switch(method) { + case 0 : {T = gsl_odeiv2_step_rk2; break; } + case 1 : {T = gsl_odeiv2_step_rk4; break; } + case 2 : {T = gsl_odeiv2_step_rkf45; break; } + case 3 : {T = gsl_odeiv2_step_rkck; break; } + case 4 : {T = gsl_odeiv2_step_rk8pd; break; } + case 5 : {T = gsl_odeiv2_step_rk2imp; break; } + case 6 : {T = gsl_odeiv2_step_rk4imp; break; } + case 7 : {T = gsl_odeiv2_step_bsimp; break; } + case 8 : {T = gsl_odeiv2_step_rk1imp; break; } + case 9 : {T = gsl_odeiv2_step_msadams; break; } + case 10: {T = gsl_odeiv2_step_msbdf; break; } + default: ERROR(BAD_CODE); + } + + Tode P; + P.f = f; + P.j = jac; + P.n = xin; + + gsl_odeiv2_system sys = {odefunc, odejac, xin, &P}; + + gsl_odeiv2_driver * d = + gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel); + + double t = tsp[0]; + + double* y = (double*)calloc(xin,sizeof(double)); + int i,j; + int status=0; + for(i=0; i< xin; i++) { + y[i] = xip[i]; + solp[i] = xip[i]; + } + + for (i = 1; i < tsn ; i++) + { + double ti = tsp[i]; + + status = gsl_odeiv2_driver_apply (d, &t, ti, y); + + if (status != GSL_SUCCESS) { + printf ("error in ode, return value=%d\n", status); + break; + } + +// printf ("%.5e %.5e %.5e\n", t, y[0], y[1]); + + for(j=0; j>> fromList [1,2,3] * fromList [3,0,-2 :: Double] + -- fromList [3.0,0.0,-6.0] + -- + -- >>> (3><3) [1..9] * ident 3 :: Matrix Double + -- (3><3) + -- [ 1.0, 0.0, 0.0 + -- , 0.0, 5.0, 0.0 + -- , 0.0, 0.0, 9.0 ] + -- + -- In arithmetic operations single-element vectors and matrices + -- (created from numeric literals or using 'scalar') automatically + -- expand to match the dimensions of the other operand: + -- + -- >>> 5 + 2*ident 3 :: Matrix Double + -- (3><3) + -- [ 7.0, 5.0, 5.0 + -- , 5.0, 7.0, 5.0 + -- , 5.0, 5.0, 7.0 ] + -- + + -- * Products + (×), + + -- | The matrix product is also implemented in the "Data.Monoid" instance for Matrix, where + -- single-element matrices (created from numeric literals or using 'scalar') + -- are used for scaling. + -- + -- >>> let m = (2><3)[1..] :: Matrix Double + -- >>> m <> 2 <> diagl[0.5,1,0] + -- (2><3) + -- [ 1.0, 4.0, 0.0 + -- , 4.0, 10.0, 0.0 ] + -- + -- mconcat uses 'optimiseMult' to get the optimal association order. + + (·), outer, kronecker, cross, + scale, + sumElements, prodElements, absSum, + + -- * Linear Systems + (<\>), + linearSolve, + linearSolveLS, + linearSolveSVD, + luSolve, + cholSolve, + + -- * Inverse and pseudoinverse + inv, pinv, pinvTol, + + -- * Determinant and rank + rcond, rank, ranksv, + det, invlndet, + + -- * Singular value decomposition + svd, + fullSVD, + thinSVD, + compactSVD, + singularValues, + leftSV, rightSV, + + -- * Eigensystems + eig, eigSH, eigSH', + eigenvalues, eigenvaluesSH, eigenvaluesSH', + geigSH', + + -- * QR + qr, rq, qrRaw, qrgr, + + -- * Cholesky + chol, cholSH, mbCholSH, + + -- * Hessenberg + hess, + + -- * Schur + schur, + + -- * LU + lu, luPacked, + + -- * Matrix functions + expm, + sqrtm, + matFunc, + + -- * Nullspace + nullspacePrec, + nullVector, + nullspaceSVD, + null1, null1sym, + + orth, + + -- * Norms + norm1, norm2, normInf, pnorm, NormType(..), + + -- * Correlation and Convolution + corr, conv, corrMin, corr2, conv2, + + -- * Random arrays + rand, randn, RandDist(..), randomVector, gaussianSample, uniformSample, + + -- * Misc + meanCov, peps, relativeError, haussholder, optimiseMult, udot, cdot, (<.>) +) where + +import Numeric.HMatrix.Data + +import Numeric.Matrix() +import Numeric.Vector() +import Numeric.Container +import Numeric.LinearAlgebra.Algorithms +import Numeric.LinearAlgebra.Util + diff --git a/packages/hmatrix/src/Numeric/HMatrix/Data.hs b/packages/hmatrix/src/Numeric/HMatrix/Data.hs new file mode 100644 index 0000000..568dc05 --- /dev/null +++ b/packages/hmatrix/src/Numeric/HMatrix/Data.hs @@ -0,0 +1,69 @@ +-------------------------------------------------------------------------------- +{- | +Module : Numeric.HMatrix.Data +Copyright : (c) Alberto Ruiz 2014 +License : GPL + +Maintainer : Alberto Ruiz +Stability : provisional + +Basic data processing. + +-} +-------------------------------------------------------------------------------- + +module Numeric.HMatrix.Data( + + -- * Vector + -- | 1D arrays are storable vectors from the vector package. + + Vector, (|>), dim, (@>), + + -- * Matrix + Matrix, (><), size, (@@>), trans, ctrans, + + -- * Construction + scalar, konst, build, assoc, accum, linspace, -- ones, zeros, + + -- * Diagonal + ident, diag, diagl, diagRect, takeDiag, + + -- * Data manipulation + fromList, toList, subVector, takesV, vjoin, + flatten, reshape, asRow, asColumn, row, col, + fromRows, toRows, fromColumns, toColumns, fromLists, toLists, fromArray2D, + takeRows, dropRows, takeColumns, dropColumns, subMatrix, (?), (¿), fliprl, flipud, + + -- * Block matrix + fromBlocks, (¦), (——), diagBlock, repmat, toBlocks, toBlocksEvery, + + -- * Mapping functions + conj, cmap, step, cond, + + -- * Find elements + find, maxIndex, minIndex, maxElement, minElement, atIndex, + + -- * IO + disp, dispf, disps, dispcf, latexFormat, format, + loadMatrix, saveMatrix, fromFile, fileDimensions, + readMatrix, + fscanfVector, fprintfVector, freadVector, fwriteVector, + +-- * Conversion + Convert(..), + + -- * Misc + arctan2, + rows, cols, + separable, + + module Data.Complex + +) where + +import Data.Packed.Vector +import Data.Packed.Matrix +import Numeric.Container +import Numeric.LinearAlgebra.Util +import Data.Complex + diff --git a/packages/hmatrix/src/Numeric/HMatrix/Devel.hs b/packages/hmatrix/src/Numeric/HMatrix/Devel.hs new file mode 100644 index 0000000..b921f44 --- /dev/null +++ b/packages/hmatrix/src/Numeric/HMatrix/Devel.hs @@ -0,0 +1,69 @@ +-------------------------------------------------------------------------------- +{- | +Module : Numeric.HMatrix.Devel +Copyright : (c) Alberto Ruiz 2014 +License : GPL + +Maintainer : Alberto Ruiz +Stability : provisional + +The library can be easily extended using the tools in this module. + +-} +-------------------------------------------------------------------------------- + +module Numeric.HMatrix.Devel( + -- * FFI helpers + -- | Sample usage, to upload a perspective matrix to a shader. + -- + -- @ glUniformMatrix4fv 0 1 (fromIntegral gl_TRUE) \`appMatrix\` perspective 0.01 100 (pi\/2) (4\/3) + -- @ + module Data.Packed.Foreign, + + -- * FFI tools + -- | Illustrative usage examples can be found + -- in the @examples\/devel@ folder included in the package. + module Data.Packed.Development, + + -- * ST + -- | In-place manipulation inside the ST monad. + -- See examples\/inplace.hs in the distribution. + + -- ** Mutable Vectors + STVector, newVector, thawVector, freezeVector, runSTVector, + readVector, writeVector, modifyVector, liftSTVector, + -- ** Mutable Matrices + STMatrix, newMatrix, thawMatrix, freezeMatrix, runSTMatrix, + readMatrix, writeMatrix, modifyMatrix, liftSTMatrix, + -- ** Unsafe functions + newUndefinedVector, + unsafeReadVector, unsafeWriteVector, + unsafeThawVector, unsafeFreezeVector, + newUndefinedMatrix, + unsafeReadMatrix, unsafeWriteMatrix, + unsafeThawMatrix, unsafeFreezeMatrix, + + -- * Special maps and zips + mapVectorWithIndex, zipVector, zipVectorWith, unzipVector, unzipVectorWith, + mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_, + foldLoop, foldVector, foldVectorG, foldVectorWithIndex, + mapMatrixWithIndex, mapMatrixWithIndexM, mapMatrixWithIndexM_, + liftMatrix, liftMatrix2, liftMatrix2Auto, + + -- * Auxiliary classes + Element, Container, Product, Contraction, LSDiv, + Complexable(), RealElement(), + RealOf, ComplexOf, SingleOf, DoubleOf, + IndexOf, + Field, Normed +) where + +import Data.Packed.Foreign +import Data.Packed.Development +import Data.Packed.ST +import Numeric.Container(Container,Contraction,LSDiv,Product, + Complexable(),RealElement(), + RealOf, ComplexOf, SingleOf, DoubleOf, IndexOf) +import Data.Packed +import Numeric.LinearAlgebra.Algorithms(Field,Normed) + diff --git a/packages/hmatrix/src/Numeric/IO.hs b/packages/hmatrix/src/Numeric/IO.hs new file mode 100644 index 0000000..836f352 --- /dev/null +++ b/packages/hmatrix/src/Numeric/IO.hs @@ -0,0 +1,165 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.IO +-- Copyright : (c) Alberto Ruiz 2010 +-- License : GPL +-- +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- Portability : portable +-- +-- Display, formatting and IO functions for numeric 'Vector' and 'Matrix' +-- +----------------------------------------------------------------------------- + +module Numeric.IO ( + dispf, disps, dispcf, vecdisp, latexFormat, format, + loadMatrix, saveMatrix, fromFile, fileDimensions, + readMatrix, fromArray2D, + fscanfVector, fprintfVector, freadVector, fwriteVector +) where + +import Data.Packed +import Data.Packed.Internal +import System.Process(readProcess) +import Text.Printf(printf) +import Data.List(intersperse) +import Data.Complex + +{- | Creates a string from a matrix given a separator and a function to show each entry. Using +this function the user can easily define any desired display function: + +@import Text.Printf(printf)@ + +@disp = putStr . format \" \" (printf \"%.2f\")@ + +-} +format :: (Element t) => String -> (t -> String) -> Matrix t -> String +format sep f m = table sep . map (map f) . toLists $ m + +{- | Show a matrix with \"autoscaling\" and a given number of decimal places. + +>>> putStr . disps 2 $ 120 * (3><4) [1..] +3x4 E3 + 0.12 0.24 0.36 0.48 + 0.60 0.72 0.84 0.96 + 1.08 1.20 1.32 1.44 + +-} +disps :: Int -> Matrix Double -> String +disps d x = sdims x ++ " " ++ formatScaled d x + +{- | Show a matrix with a given number of decimal places. + +>>> dispf 2 (1/3 + ident 3) +"3x3\n1.33 0.33 0.33\n0.33 1.33 0.33\n0.33 0.33 1.33\n" + +>>> putStr . dispf 2 $ (3><4)[1,1.5..] +3x4 +1.00 1.50 2.00 2.50 +3.00 3.50 4.00 4.50 +5.00 5.50 6.00 6.50 + +>>> putStr . unlines . tail . lines . dispf 2 . asRow $ linspace 10 (0,1) +0.00 0.11 0.22 0.33 0.44 0.56 0.67 0.78 0.89 1.00 + +-} +dispf :: Int -> Matrix Double -> String +dispf d x = sdims x ++ "\n" ++ formatFixed (if isInt x then 0 else d) x + +sdims x = show (rows x) ++ "x" ++ show (cols x) + +formatFixed d x = format " " (printf ("%."++show d++"f")) $ x + +isInt = all lookslikeInt . toList . flatten + +formatScaled dec t = "E"++show o++"\n" ++ ss + where ss = format " " (printf fmt. g) t + g x | o >= 0 = x/10^(o::Int) + | otherwise = x*10^(-o) + o | rows t == 0 || cols t == 0 = 0 + | otherwise = floor $ maximum $ map (logBase 10 . abs) $ toList $ flatten t + fmt = '%':show (dec+3) ++ '.':show dec ++"f" + +{- | Show a vector using a function for showing matrices. + +>>> putStr . vecdisp (dispf 2) $ linspace 10 (0,1) +10 |> 0.00 0.11 0.22 0.33 0.44 0.56 0.67 0.78 0.89 1.00 + +-} +vecdisp :: (Element t) => (Matrix t -> String) -> Vector t -> String +vecdisp f v + = ((show (dim v) ++ " |> ") ++) . (++"\n") + . unwords . lines . tail . dropWhile (not . (`elem` " \n")) + . f . trans . reshape 1 + $ v + +{- | Tool to display matrices with latex syntax. + +>>> latexFormat "bmatrix" (dispf 2 $ ident 2) +"\\begin{bmatrix}\n1 & 0\n\\\\\n0 & 1\n\\end{bmatrix}" + +-} +latexFormat :: String -- ^ type of braces: \"matrix\", \"bmatrix\", \"pmatrix\", etc. + -> String -- ^ Formatted matrix, with elements separated by spaces and newlines + -> String +latexFormat del tab = "\\begin{"++del++"}\n" ++ f tab ++ "\\end{"++del++"}" + where f = unlines . intersperse "\\\\" . map unwords . map (intersperse " & " . words) . tail . lines + +-- | Pretty print a complex number with at most n decimal digits. +showComplex :: Int -> Complex Double -> String +showComplex d (a:+b) + | isZero a && isZero b = "0" + | isZero b = sa + | isZero a && isOne b = s2++"i" + | isZero a = sb++"i" + | isOne b = sa++s3++"i" + | otherwise = sa++s1++sb++"i" + where + sa = shcr d a + sb = shcr d b + s1 = if b<0 then "" else "+" + s2 = if b<0 then "-" else "" + s3 = if b<0 then "-" else "+" + +shcr d a | lookslikeInt a = printf "%.0f" a + | otherwise = printf ("%."++show d++"f") a + + +lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx + where shx = show x + +isZero x = show x `elem` ["0.0","-0.0"] +isOne x = show x `elem` ["1.0","-1.0"] + +-- | Pretty print a complex matrix with at most n decimal digits. +dispcf :: Int -> Matrix (Complex Double) -> String +dispcf d m = sdims m ++ "\n" ++ format " " (showComplex d) m + +-------------------------------------------------------------------- + +-- | reads a matrix from a string containing a table of numbers. +readMatrix :: String -> Matrix Double +readMatrix = fromLists . map (map read). map words . filter (not.null) . lines + +{- | obtains the number of rows and columns in an ASCII data file + (provisionally using unix's wc). +-} +fileDimensions :: FilePath -> IO (Int,Int) +fileDimensions fname = do + wcres <- readProcess "wc" ["-w",fname] "" + contents <- readFile fname + let tot = read . head . words $ wcres + c = length . head . dropWhile null . map words . lines $ contents + if tot > 0 + then return (tot `div` c, c) + else return (0,0) + +-- | Loads a matrix from an ASCII file formatted as a 2D table. +loadMatrix :: FilePath -> IO (Matrix Double) +loadMatrix file = fromFile file =<< fileDimensions file + +-- | Loads a matrix from an ASCII file (the number of rows and columns must be known in advance). +fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double) +fromFile filename (r,c) = reshape c `fmap` fscanfVector filename (r*c) + diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra.hs b/packages/hmatrix/src/Numeric/LinearAlgebra.hs new file mode 100644 index 0000000..1db860c --- /dev/null +++ b/packages/hmatrix/src/Numeric/LinearAlgebra.hs @@ -0,0 +1,30 @@ +----------------------------------------------------------------------------- +{- | +Module : Numeric.LinearAlgebra +Copyright : (c) Alberto Ruiz 2006-10 +License : GPL-style + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional +Portability : uses ffi + +This module reexports all normally required functions for Linear Algebra applications. + +It also provides instances of standard classes 'Show', 'Read', 'Eq', +'Num', 'Fractional', and 'Floating' for 'Vector' and 'Matrix'. +In arithmetic operations one-component vectors and matrices automatically +expand to match the dimensions of the other operand. + +-} +----------------------------------------------------------------------------- +{-# OPTIONS_HADDOCK hide #-} + +module Numeric.LinearAlgebra ( + module Numeric.Container, + module Numeric.LinearAlgebra.Algorithms +) where + +import Numeric.Container +import Numeric.LinearAlgebra.Algorithms +import Numeric.Matrix() +import Numeric.Vector() diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/Algorithms.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/Algorithms.hs new file mode 100644 index 0000000..8c4b610 --- /dev/null +++ b/packages/hmatrix/src/Numeric/LinearAlgebra/Algorithms.hs @@ -0,0 +1,746 @@ +{-# LANGUAGE FlexibleContexts, FlexibleInstances #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE TypeFamilies #-} + +----------------------------------------------------------------------------- +{- | +Module : Numeric.LinearAlgebra.Algorithms +Copyright : (c) Alberto Ruiz 2006-9 +License : GPL-style + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional +Portability : uses ffi + +High level generic interface to common matrix computations. + +Specific functions for particular base types can also be explicitly +imported from "Numeric.LinearAlgebra.LAPACK". + +-} +----------------------------------------------------------------------------- +{-# OPTIONS_HADDOCK hide #-} + + +module Numeric.LinearAlgebra.Algorithms ( +-- * Supported types + Field(), +-- * Linear Systems + linearSolve, + luSolve, + cholSolve, + linearSolveLS, + linearSolveSVD, + inv, pinv, pinvTol, + det, invlndet, + rank, rcond, +-- * Matrix factorizations +-- ** Singular value decomposition + svd, + fullSVD, + thinSVD, + compactSVD, + singularValues, + leftSV, rightSV, +-- ** Eigensystems + eig, eigSH, eigSH', + eigenvalues, eigenvaluesSH, eigenvaluesSH', + geigSH', +-- ** QR + qr, rq, qrRaw, qrgr, +-- ** Cholesky + chol, cholSH, mbCholSH, +-- ** Hessenberg + hess, +-- ** Schur + schur, +-- ** LU + lu, luPacked, +-- * Matrix functions + expm, + sqrtm, + matFunc, +-- * Nullspace + nullspacePrec, + nullVector, + nullspaceSVD, + orth, +-- * Norms + Normed(..), NormType(..), + relativeError, +-- * Misc + eps, peps, i, +-- * Util + haussholder, + unpackQR, unpackHess, + ranksv +) where + + +import Data.Packed.Internal hiding ((//)) +import Data.Packed.Matrix +import Numeric.LinearAlgebra.LAPACK as LAPACK +import Data.List(foldl1') +import Data.Array +import Numeric.ContainerBoot + + +{- | Generic linear algebra functions for double precision real and complex matrices. + +(Single precision data can be converted using 'single' and 'double'). + +-} +class (Product t, + Convert t, + Container Vector t, + Container Matrix t, + Normed Matrix t, + Normed Vector t, + Floating t, + RealOf t ~ Double) => Field t where + svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t) + thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t) + sv' :: Matrix t -> Vector Double + luPacked' :: Matrix t -> (Matrix t, [Int]) + luSolve' :: (Matrix t, [Int]) -> Matrix t -> Matrix t + linearSolve' :: Matrix t -> Matrix t -> Matrix t + cholSolve' :: Matrix t -> Matrix t -> Matrix t + linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t + linearSolveLS' :: Matrix t -> Matrix t -> Matrix t + eig' :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) + eigSH'' :: Matrix t -> (Vector Double, Matrix t) + eigOnly :: Matrix t -> Vector (Complex Double) + eigOnlySH :: Matrix t -> Vector Double + cholSH' :: Matrix t -> Matrix t + mbCholSH' :: Matrix t -> Maybe (Matrix t) + qr' :: Matrix t -> (Matrix t, Vector t) + qrgr' :: Int -> (Matrix t, Vector t) -> Matrix t + hess' :: Matrix t -> (Matrix t, Matrix t) + schur' :: Matrix t -> (Matrix t, Matrix t) + + +instance Field Double where + svd' = svdRd + thinSVD' = thinSVDRd + sv' = svR + luPacked' = luR + luSolve' (l_u,perm) = lusR l_u perm + linearSolve' = linearSolveR -- (luSolve . luPacked) ?? + cholSolve' = cholSolveR + linearSolveLS' = linearSolveLSR + linearSolveSVD' = linearSolveSVDR Nothing + eig' = eigR + eigSH'' = eigS + eigOnly = eigOnlyR + eigOnlySH = eigOnlyS + cholSH' = cholS + mbCholSH' = mbCholS + qr' = qrR + qrgr' = qrgrR + hess' = unpackHess hessR + schur' = schurR + +instance Field (Complex Double) where +#ifdef NOZGESDD + svd' = svdC + thinSVD' = thinSVDC +#else + svd' = svdCd + thinSVD' = thinSVDCd +#endif + sv' = svC + luPacked' = luC + luSolve' (l_u,perm) = lusC l_u perm + linearSolve' = linearSolveC + cholSolve' = cholSolveC + linearSolveLS' = linearSolveLSC + linearSolveSVD' = linearSolveSVDC Nothing + eig' = eigC + eigOnly = eigOnlyC + eigSH'' = eigH + eigOnlySH = eigOnlyH + cholSH' = cholH + mbCholSH' = mbCholH + qr' = qrC + qrgr' = qrgrC + hess' = unpackHess hessC + schur' = schurC + +-------------------------------------------------------------- + +square m = rows m == cols m + +vertical m = rows m >= cols m + +exactHermitian m = m `equal` ctrans m + +-------------------------------------------------------------- + +-- | Full singular value decomposition. +svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) +svd = {-# SCC "svd" #-} svd' + +-- | A version of 'svd' which returns only the @min (rows m) (cols m)@ singular vectors of @m@. +-- +-- If @(u,s,v) = thinSVD m@ then @m == u \<> diag s \<> trans v@. +thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) +thinSVD = {-# SCC "thinSVD" #-} thinSVD' + +-- | Singular values only. +singularValues :: Field t => Matrix t -> Vector Double +singularValues = {-# SCC "singularValues" #-} sv' + +-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values. +-- +-- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> trans v@. +fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t) +fullSVD m = (u,d,v) where + (u,s,v) = svd m + d = diagRect 0 s r c + r = rows m + c = cols m + +-- | Similar to 'thinSVD', returning only the nonzero singular values and the corresponding singular vectors. +compactSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) +compactSVD m = (u', subVector 0 d s, v') where + (u,s,v) = thinSVD m + d = rankSVD (1*eps) m s `max` 1 + u' = takeColumns d u + v' = takeColumns d v + + +-- | Singular values and all right singular vectors. +rightSV :: Field t => Matrix t -> (Vector Double, Matrix t) +rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v) + | otherwise = let (_,s,v) = svd m in (s,v) + +-- | Singular values and all left singular vectors. +leftSV :: Field t => Matrix t -> (Matrix t, Vector Double) +leftSV m | vertical m = let (u,s,_) = svd m in (u,s) + | otherwise = let (u,s,_) = thinSVD m in (u,s) + + +-------------------------------------------------------------- + +-- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'. +luPacked :: Field t => Matrix t -> (Matrix t, [Int]) +luPacked = {-# SCC "luPacked" #-} luPacked' + +-- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization obtained by 'luPacked'. +luSolve :: Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t +luSolve = {-# SCC "luSolve" #-} luSolve' + +-- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'. +-- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system. +linearSolve :: Field t => Matrix t -> Matrix t -> Matrix t +linearSolve = {-# SCC "linearSolve" #-} linearSolve' + +-- | Solve a symmetric or Hermitian positive definite linear system using a precomputed Cholesky decomposition obtained by 'chol'. +cholSolve :: Field t => Matrix t -> Matrix t -> Matrix t +cholSolve = {-# SCC "cholSolve" #-} cholSolve' + +-- | Minimum norm solution of a general linear least squares problem Ax=B using the SVD. Admits rank-deficient systems but it is slower than 'linearSolveLS'. The effective rank of A is determined by treating as zero those singular valures which are less than 'eps' times the largest singular value. +linearSolveSVD :: Field t => Matrix t -> Matrix t -> Matrix t +linearSolveSVD = {-# SCC "linearSolveSVD" #-} linearSolveSVD' + + +-- | Least squared error solution of an overconstrained linear system, or the minimum norm solution of an underconstrained system. For rank-deficient systems use 'linearSolveSVD'. +linearSolveLS :: Field t => Matrix t -> Matrix t -> Matrix t +linearSolveLS = {-# SCC "linearSolveLS" #-} linearSolveLS' + +-------------------------------------------------------------- + +-- | Eigenvalues and eigenvectors of a general square matrix. +-- +-- If @(s,v) = eig m@ then @m \<> v == v \<> diag s@ +eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) +eig = {-# SCC "eig" #-} eig' + +-- | Eigenvalues of a general square matrix. +eigenvalues :: Field t => Matrix t -> Vector (Complex Double) +eigenvalues = {-# SCC "eigenvalues" #-} eigOnly + +-- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part. +eigSH' :: Field t => Matrix t -> (Vector Double, Matrix t) +eigSH' = {-# SCC "eigSH'" #-} eigSH'' + +-- | Similar to 'eigenvaluesSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part. +eigenvaluesSH' :: Field t => Matrix t -> Vector Double +eigenvaluesSH' = {-# SCC "eigenvaluesSH'" #-} eigOnlySH + +-- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix. +-- +-- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@ +eigSH :: Field t => Matrix t -> (Vector Double, Matrix t) +eigSH m | exactHermitian m = eigSH' m + | otherwise = error "eigSH requires complex hermitian or real symmetric matrix" + +-- | Eigenvalues of a complex hermitian or real symmetric matrix. +eigenvaluesSH :: Field t => Matrix t -> Vector Double +eigenvaluesSH m | exactHermitian m = eigenvaluesSH' m + | otherwise = error "eigenvaluesSH requires complex hermitian or real symmetric matrix" + +-------------------------------------------------------------- + +-- | QR factorization. +-- +-- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular. +qr :: Field t => Matrix t -> (Matrix t, Matrix t) +qr = {-# SCC "qr" #-} unpackQR . qr' + +qrRaw m = qr' m + +{- | generate a matrix with k orthogonal columns from the output of qrRaw +-} +qrgr n (a,t) + | dim t > min (cols a) (rows a) || n < 0 || n > dim t = error "qrgr expects k <= min(rows,cols)" + | otherwise = qrgr' n (a,t) + +-- | RQ factorization. +-- +-- If @(r,q) = rq m@ then @m == r \<> q@, where q is unitary and r is upper triangular. +rq :: Field t => Matrix t -> (Matrix t, Matrix t) +rq m = {-# SCC "rq" #-} (r,q) where + (q',r') = qr $ trans $ rev1 m + r = rev2 (trans r') + q = rev2 (trans q') + rev1 = flipud . fliprl + rev2 = fliprl . flipud + +-- | Hessenberg factorization. +-- +-- If @(p,h) = hess m@ then @m == p \<> h \<> ctrans p@, where p is unitary +-- and h is in upper Hessenberg form (it has zero entries below the first subdiagonal). +hess :: Field t => Matrix t -> (Matrix t, Matrix t) +hess = hess' + +-- | Schur factorization. +-- +-- If @(u,s) = schur m@ then @m == u \<> s \<> ctrans u@, where u is unitary +-- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is +-- upper triangular in 2x2 blocks. +-- +-- \"Anything that the Jordan decomposition can do, the Schur decomposition +-- can do better!\" (Van Loan) +schur :: Field t => Matrix t -> (Matrix t, Matrix t) +schur = schur' + + +-- | Similar to 'cholSH', but instead of an error (e.g., caused by a matrix not positive definite) it returns 'Nothing'. +mbCholSH :: Field t => Matrix t -> Maybe (Matrix t) +mbCholSH = {-# SCC "mbCholSH" #-} mbCholSH' + +-- | Similar to 'chol', without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part. +cholSH :: Field t => Matrix t -> Matrix t +cholSH = {-# SCC "cholSH" #-} cholSH' + +-- | Cholesky factorization of a positive definite hermitian or symmetric matrix. +-- +-- If @c = chol m@ then @c@ is upper triangular and @m == ctrans c \<> c@. +chol :: Field t => Matrix t -> Matrix t +chol m | exactHermitian m = cholSH m + | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix" + + +-- | Joint computation of inverse and logarithm of determinant of a square matrix. +invlndet :: Field t + => Matrix t + -> (Matrix t, (t, t)) -- ^ (inverse, (log abs det, sign or phase of det)) +invlndet m | square m = (im,(ladm,sdm)) + | otherwise = error $ "invlndet of nonsquare "++ shSize m ++ " matrix" + where + lp@(lup,perm) = luPacked m + s = signlp (rows m) perm + dg = toList $ takeDiag $ lup + ladm = sum $ map (log.abs) dg + sdm = s* product (map signum dg) + im = luSolve lp (ident (rows m)) + + +-- | Determinant of a square matrix. To avoid possible overflow or underflow use 'invlndet'. +det :: Field t => Matrix t -> t +det m | square m = {-# SCC "det" #-} s * (product $ toList $ takeDiag $ lup) + | otherwise = error $ "det of nonsquare "++ shSize m ++ " matrix" + where (lup,perm) = luPacked m + s = signlp (rows m) perm + +-- | Explicit LU factorization of a general matrix. +-- +-- If @(l,u,p,s) = lu m@ then @m == p \<> l \<> u@, where l is lower triangular, +-- u is upper triangular, p is a permutation matrix and s is the signature of the permutation. +lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t) +lu = luFact . luPacked + +-- | Inverse of a square matrix. See also 'invlndet'. +inv :: Field t => Matrix t -> Matrix t +inv m | square m = m `linearSolve` ident (rows m) + | otherwise = error $ "inv of nonsquare "++ shSize m ++ " matrix" + + +-- | Pseudoinverse of a general matrix with default tolerance ('pinvTol' 1, similar to GNU-Octave). +pinv :: Field t => Matrix t -> Matrix t +pinv = pinvTol 1 + +{- | @pinvTol r@ computes the pseudoinverse of a matrix with tolerance @tol=r*g*eps*(max rows cols)@, where g is the greatest singular value. + +@ +m = (3><3) [ 1, 0, 0 + , 0, 1, 0 + , 0, 0, 1e-10] :: Matrix Double +@ + +>>> pinv m +1. 0. 0. +0. 1. 0. +0. 0. 10000000000. + +>>> pinvTol 1E8 m +1. 0. 0. +0. 1. 0. +0. 0. 1. + +-} + +pinvTol :: Field t => Double -> Matrix t -> Matrix t +pinvTol t m = conj v' `mXm` diag s' `mXm` ctrans u' where + (u,s,v) = thinSVD m + sl@(g:_) = toList s + s' = real . fromList . map rec $ sl + rec x = if x <= g*tol then x else 1/x + tol = (fromIntegral (max r c) * g * t * eps) + r = rows m + c = cols m + d = dim s + u' = takeColumns d u + v' = takeColumns d v + + +-- | Numeric rank of a matrix from the SVD decomposition. +rankSVD :: Element t + => Double -- ^ numeric zero (e.g. 1*'eps') + -> Matrix t -- ^ input matrix m + -> Vector Double -- ^ 'sv' of m + -> Int -- ^ rank of m +rankSVD teps m s = ranksv teps (max (rows m) (cols m)) (toList s) + +-- | Numeric rank of a matrix from its singular values. +ranksv :: Double -- ^ numeric zero (e.g. 1*'eps') + -> Int -- ^ maximum dimension of the matrix + -> [Double] -- ^ singular values + -> Int -- ^ rank of m +ranksv teps maxdim s = k where + g = maximum s + tol = fromIntegral maxdim * g * teps + s' = filter (>tol) s + k = if g > teps then length s' else 0 + +-- | The machine precision of a Double: @eps = 2.22044604925031e-16@ (the value used by GNU-Octave). +eps :: Double +eps = 2.22044604925031e-16 + + +-- | 1 + 0.5*peps == 1, 1 + 0.6*peps /= 1 +peps :: RealFloat x => x +peps = x where x = 2.0 ** fromIntegral (1 - floatDigits x) + + +-- | The imaginary unit: @i = 0.0 :+ 1.0@ +i :: Complex Double +i = 0:+1 + +----------------------------------------------------------------------- + +-- | The nullspace of a matrix from its precomputed SVD decomposition. +nullspaceSVD :: Field t + => Either Double Int -- ^ Left \"numeric\" zero (eg. 1*'eps'), + -- or Right \"theoretical\" matrix rank. + -> Matrix t -- ^ input matrix m + -> (Vector Double, Matrix t) -- ^ 'rightSV' of m + -> [Vector t] -- ^ list of unitary vectors spanning the nullspace +nullspaceSVD hint a (s,v) = vs where + tol = case hint of + Left t -> t + _ -> eps + k = case hint of + Right t -> t + _ -> rankSVD tol a s + vs = drop k $ toRows $ ctrans v + + +-- | The nullspace of a matrix. See also 'nullspaceSVD'. +nullspacePrec :: Field t + => Double -- ^ relative tolerance in 'eps' units (e.g., use 3 to get 3*'eps') + -> Matrix t -- ^ input matrix + -> [Vector t] -- ^ list of unitary vectors spanning the nullspace +nullspacePrec t m = nullspaceSVD (Left (t*eps)) m (rightSV m) + +-- | The nullspace of a matrix, assumed to be one-dimensional, with machine precision. +nullVector :: Field t => Matrix t -> Vector t +nullVector = last . nullspacePrec 1 + +orth :: Field t => Matrix t -> [Vector t] +-- ^ Return an orthonormal basis of the range space of a matrix +orth m = take r $ toColumns u + where + (u,s,_) = compactSVD m + r = ranksv eps (max (rows m) (cols m)) (toList s) + +------------------------------------------------------------------------ + +-- many thanks, quickcheck! + +haussholder :: (Field a) => a -> Vector a -> Matrix a +haussholder tau v = ident (dim v) `sub` (tau `scale` (w `mXm` ctrans w)) + where w = asColumn v + + +zh k v = fromList $ replicate (k-1) 0 ++ (1:drop k xs) + where xs = toList v + +zt 0 v = v +zt k v = vjoin [subVector 0 (dim v - k) v, konst' 0 k] + + +unpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t) +unpackQR (pq, tau) = {-# SCC "unpackQR" #-} (q,r) + where cs = toColumns pq + m = rows pq + n = cols pq + mn = min m n + r = fromColumns $ zipWith zt ([m-1, m-2 .. 1] ++ repeat 0) cs + vs = zipWith zh [1..mn] cs + hs = zipWith haussholder (toList tau) vs + q = foldl1' mXm hs + +unpackHess :: (Field t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t) +unpackHess hf m + | rows m == 1 = ((1><1)[1],m) + | otherwise = (uH . hf) m + +uH (pq, tau) = (p,h) + where cs = toColumns pq + m = rows pq + n = cols pq + mn = min m n + h = fromColumns $ zipWith zt ([m-2, m-3 .. 1] ++ repeat 0) cs + vs = zipWith zh [2..mn] cs + hs = zipWith haussholder (toList tau) vs + p = foldl1' mXm hs + +-------------------------------------------------------------------------- + +-- | Reciprocal of the 2-norm condition number of a matrix, computed from the singular values. +rcond :: Field t => Matrix t -> Double +rcond m = last s / head s + where s = toList (singularValues m) + +-- | Number of linearly independent rows or columns. +rank :: Field t => Matrix t -> Int +rank m = rankSVD eps m (singularValues m) + +{- +expm' m = case diagonalize (complex m) of + Just (l,v) -> v `mXm` diag (exp l) `mXm` inv v + Nothing -> error "Sorry, expm not yet implemented for non-diagonalizable matrices" + where exp = vectorMapC Exp +-} + +diagonalize m = if rank v == n + then Just (l,v) + else Nothing + where n = rows m + (l,v) = if exactHermitian m + then let (l',v') = eigSH m in (real l', v') + else eig m + +-- | Generic matrix functions for diagonalizable matrices. For instance: +-- +-- @logm = matFunc log@ +-- +matFunc :: (Complex Double -> Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) +matFunc f m = case diagonalize m of + Just (l,v) -> v `mXm` diag (mapVector f l) `mXm` inv v + Nothing -> error "Sorry, matFunc requires a diagonalizable matrix" + +-------------------------------------------------------------- + +golubeps :: Integer -> Integer -> Double +golubeps p q = a * fromIntegral b / fromIntegral c where + a = 2^^(3-p-q) + b = fact p * fact q + c = fact (p+q) * fact (p+q+1) + fact n = product [1..n] + +epslist :: [(Int,Double)] +epslist = [ (fromIntegral k, golubeps k k) | k <- [1..]] + +geps delta = head [ k | (k,g) <- epslist, g Matrix t -> Matrix t +expm = expGolub + +expGolub :: Field t => Matrix t -> Matrix t +expGolub m = iterate msq f !! j + where j = max 0 $ floor $ logBase 2 $ pnorm Infinity m + a = m */ fromIntegral ((2::Int)^j) + q = geps eps -- 7 steps + eye = ident (rows m) + work (k,c,x,n,d) = (k',c',x',n',d') + where k' = k+1 + c' = c * fromIntegral (q-k+1) / fromIntegral ((2*q-k+1)*k) + x' = a <> x + n' = n |+| (c' .* x') + d' = d |+| (((-1)^k * c') .* x') + (_,_,_,nf,df) = iterate work (1,1,eye,eye,eye) !! q + f = linearSolve df nf + msq x = x <> x + + (<>) = multiply + v */ x = scale (recip x) v + (.*) = scale + (|+|) = add + +-------------------------------------------------------------- + +{- | Matrix square root. Currently it uses a simple iterative algorithm described in Wikipedia. +It only works with invertible matrices that have a real solution. For diagonalizable matrices you can try @matFunc sqrt@. + +@m = (2><2) [4,9 + ,0,4] :: Matrix Double@ + +>>> sqrtm m +(2><2) + [ 2.0, 2.25 + , 0.0, 2.0 ] + +-} +sqrtm :: Field t => Matrix t -> Matrix t +sqrtm = sqrtmInv + +sqrtmInv x = fst $ fixedPoint $ iterate f (x, ident (rows x)) + where fixedPoint (a:b:rest) | pnorm PNorm1 (fst a |-| fst b) < peps = a + | otherwise = fixedPoint (b:rest) + fixedPoint _ = error "fixedpoint with impossible inputs" + f (y,z) = (0.5 .* (y |+| inv z), + 0.5 .* (inv y |+| z)) + (.*) = scale + (|+|) = add + (|-|) = sub + +------------------------------------------------------------------ + +signlp r vals = foldl f 1 (zip [0..r-1] vals) + where f s (a,b) | a /= b = -s + | otherwise = s + +swap (arr,s) (a,b) | a /= b = (arr // [(a, arr!b),(b,arr!a)],-s) + | otherwise = (arr,s) + +fixPerm r vals = (fromColumns $ elems res, sign) + where v = [0..r-1] + s = toColumns (ident r) + (res,sign) = foldl swap (listArray (0,r-1) s, 1) (zip v vals) + +triang r c h v = (r>=h then v else 1 - v + +luFact (l_u,perm) | r <= c = (l ,u ,p, s) + | otherwise = (l',u',p, s) + where + r = rows l_u + c = cols l_u + tu = triang r c 0 1 + tl = triang r c 0 0 + l = takeColumns r (l_u |*| tl) |+| diagRect 0 (konst' 1 r) r r + u = l_u |*| tu + (p,s) = fixPerm r perm + l' = (l_u |*| tl) |+| diagRect 0 (konst' 1 c) r c + u' = takeRows c (l_u |*| tu) + (|+|) = add + (|*|) = mul + +--------------------------------------------------------------------------- + +data NormType = Infinity | PNorm1 | PNorm2 | Frobenius + +class (RealFloat (RealOf t)) => Normed c t where + pnorm :: NormType -> c t -> RealOf t + +instance Normed Vector Double where + pnorm PNorm1 = norm1 + pnorm PNorm2 = norm2 + pnorm Infinity = normInf + pnorm Frobenius = norm2 + +instance Normed Vector (Complex Double) where + pnorm PNorm1 = norm1 + pnorm PNorm2 = norm2 + pnorm Infinity = normInf + pnorm Frobenius = pnorm PNorm2 + +instance Normed Vector Float where + pnorm PNorm1 = norm1 + pnorm PNorm2 = norm2 + pnorm Infinity = normInf + pnorm Frobenius = pnorm PNorm2 + +instance Normed Vector (Complex Float) where + pnorm PNorm1 = norm1 + pnorm PNorm2 = norm2 + pnorm Infinity = normInf + pnorm Frobenius = pnorm PNorm2 + + +instance Normed Matrix Double where + pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns + pnorm PNorm2 = (@>0) . singularValues + pnorm Infinity = pnorm PNorm1 . trans + pnorm Frobenius = pnorm PNorm2 . flatten + +instance Normed Matrix (Complex Double) where + pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns + pnorm PNorm2 = (@>0) . singularValues + pnorm Infinity = pnorm PNorm1 . trans + pnorm Frobenius = pnorm PNorm2 . flatten + +instance Normed Matrix Float where + pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns + pnorm PNorm2 = realToFrac . (@>0) . singularValues . double + pnorm Infinity = pnorm PNorm1 . trans + pnorm Frobenius = pnorm PNorm2 . flatten + +instance Normed Matrix (Complex Float) where + pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns + pnorm PNorm2 = realToFrac . (@>0) . singularValues . double + pnorm Infinity = pnorm PNorm1 . trans + pnorm Frobenius = pnorm PNorm2 . flatten + +-- | Approximate number of common digits in the maximum element. +relativeError :: (Normed c t, Container c t) => c t -> c t -> Int +relativeError x y = dig (norm (x `sub` y) / norm x) + where norm = pnorm Infinity + dig r = round $ -logBase 10 (realToFrac r :: Double) + +---------------------------------------------------------------------- + +-- | Generalized symmetric positive definite eigensystem Av = lBv, +-- for A and B symmetric, B positive definite (conditions not checked). +geigSH' :: Field t + => Matrix t -- ^ A + -> Matrix t -- ^ B + -> (Vector Double, Matrix t) +geigSH' a b = (l,v') + where + u = cholSH b + iu = inv u + c = ctrans iu <> a <> iu + (l,v) = eigSH' c + v' = iu <> v + (<>) = mXm + diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK.hs new file mode 100644 index 0000000..11394a6 --- /dev/null +++ b/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK.hs @@ -0,0 +1,555 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.LinearAlgebra.LAPACK +-- Copyright : (c) Alberto Ruiz 2006-7 +-- License : GPL-style +-- +-- Maintainer : Alberto Ruiz (aruiz at um dot es) +-- Stability : provisional +-- Portability : portable (uses FFI) +-- +-- Functional interface to selected LAPACK functions (). +-- +----------------------------------------------------------------------------- +{-# OPTIONS_HADDOCK hide #-} + +module Numeric.LinearAlgebra.LAPACK ( + -- * Matrix product + multiplyR, multiplyC, multiplyF, multiplyQ, + -- * Linear systems + linearSolveR, linearSolveC, + lusR, lusC, + cholSolveR, cholSolveC, + linearSolveLSR, linearSolveLSC, + linearSolveSVDR, linearSolveSVDC, + -- * SVD + svR, svRd, svC, svCd, + svdR, svdRd, svdC, svdCd, + thinSVDR, thinSVDRd, thinSVDC, thinSVDCd, + rightSVR, rightSVC, leftSVR, leftSVC, + -- * Eigensystems + eigR, eigC, eigS, eigS', eigH, eigH', + eigOnlyR, eigOnlyC, eigOnlyS, eigOnlyH, + -- * LU + luR, luC, + -- * Cholesky + cholS, cholH, mbCholS, mbCholH, + -- * QR + qrR, qrC, qrgrR, qrgrC, + -- * Hessenberg + hessR, hessC, + -- * Schur + schurR, schurC +) where + +import Data.Packed.Internal +import Data.Packed.Matrix +import Numeric.Conversion +import Numeric.GSL.Vector(vectorMapValR, FunCodeSV(Scale)) + +import Foreign.Ptr(nullPtr) +import Foreign.C.Types +import Control.Monad(when) +import System.IO.Unsafe(unsafePerformIO) + +----------------------------------------------------------------------------------- + +foreign import ccall unsafe "multiplyR" dgemmc :: CInt -> CInt -> TMMM +foreign import ccall unsafe "multiplyC" zgemmc :: CInt -> CInt -> TCMCMCM +foreign import ccall unsafe "multiplyF" sgemmc :: CInt -> CInt -> TFMFMFM +foreign import ccall unsafe "multiplyQ" cgemmc :: CInt -> CInt -> TQMQMQM + +isT Matrix{order = ColumnMajor} = 0 +isT Matrix{order = RowMajor} = 1 + +tt x@Matrix{order = ColumnMajor} = x +tt x@Matrix{order = RowMajor} = trans x + +multiplyAux f st a b = unsafePerformIO $ do + when (cols a /= rows b) $ error $ "inconsistent dimensions in matrix product "++ + show (rows a,cols a) ++ " x " ++ show (rows b, cols b) + s <- createMatrix ColumnMajor (rows a) (cols b) + app3 (f (isT a) (isT b)) mat (tt a) mat (tt b) mat s st + return s + +-- | Matrix product based on BLAS's /dgemm/. +multiplyR :: Matrix Double -> Matrix Double -> Matrix Double +multiplyR a b = {-# SCC "multiplyR" #-} multiplyAux dgemmc "dgemmc" a b + +-- | Matrix product based on BLAS's /zgemm/. +multiplyC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) +multiplyC a b = multiplyAux zgemmc "zgemmc" a b + +-- | Matrix product based on BLAS's /sgemm/. +multiplyF :: Matrix Float -> Matrix Float -> Matrix Float +multiplyF a b = multiplyAux sgemmc "sgemmc" a b + +-- | Matrix product based on BLAS's /cgemm/. +multiplyQ :: Matrix (Complex Float) -> Matrix (Complex Float) -> Matrix (Complex Float) +multiplyQ a b = multiplyAux cgemmc "cgemmc" a b + +----------------------------------------------------------------------------- +foreign import ccall unsafe "svd_l_R" dgesvd :: TMMVM +foreign import ccall unsafe "svd_l_C" zgesvd :: TCMCMVCM +foreign import ccall unsafe "svd_l_Rdd" dgesdd :: TMMVM +foreign import ccall unsafe "svd_l_Cdd" zgesdd :: TCMCMVCM + +-- | Full SVD of a real matrix using LAPACK's /dgesvd/. +svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) +svdR = svdAux dgesvd "svdR" . fmat + +-- | Full SVD of a real matrix using LAPACK's /dgesdd/. +svdRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) +svdRd = svdAux dgesdd "svdRdd" . fmat + +-- | Full SVD of a complex matrix using LAPACK's /zgesvd/. +svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) +svdC = svdAux zgesvd "svdC" . fmat + +-- | Full SVD of a complex matrix using LAPACK's /zgesdd/. +svdCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) +svdCd = svdAux zgesdd "svdCdd" . fmat + +svdAux f st x = unsafePerformIO $ do + u <- createMatrix ColumnMajor r r + s <- createVector (min r c) + v <- createMatrix ColumnMajor c c + app4 f mat x mat u vec s mat v st + return (u,s,trans v) + where r = rows x + c = cols x + + +-- | Thin SVD of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'S\'. +thinSVDR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) +thinSVDR = thinSVDAux dgesvd "thinSVDR" . fmat + +-- | Thin SVD of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'S\'. +thinSVDC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) +thinSVDC = thinSVDAux zgesvd "thinSVDC" . fmat + +-- | Thin SVD of a real matrix, using LAPACK's /dgesdd/ with jobz == \'S\'. +thinSVDRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) +thinSVDRd = thinSVDAux dgesdd "thinSVDRdd" . fmat + +-- | Thin SVD of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'S\'. +thinSVDCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) +thinSVDCd = thinSVDAux zgesdd "thinSVDCdd" . fmat + +thinSVDAux f st x = unsafePerformIO $ do + u <- createMatrix ColumnMajor r q + s <- createVector q + v <- createMatrix ColumnMajor q c + app4 f mat x mat u vec s mat v st + return (u,s,trans v) + where r = rows x + c = cols x + q = min r c + + +-- | Singular values of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'N\'. +svR :: Matrix Double -> Vector Double +svR = svAux dgesvd "svR" . fmat + +-- | Singular values of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'N\'. +svC :: Matrix (Complex Double) -> Vector Double +svC = svAux zgesvd "svC" . fmat + +-- | Singular values of a real matrix, using LAPACK's /dgesdd/ with jobz == \'N\'. +svRd :: Matrix Double -> Vector Double +svRd = svAux dgesdd "svRd" . fmat + +-- | Singular values of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'N\'. +svCd :: Matrix (Complex Double) -> Vector Double +svCd = svAux zgesdd "svCd" . fmat + +svAux f st x = unsafePerformIO $ do + s <- createVector q + app2 g mat x vec s st + return s + where r = rows x + c = cols x + q = min r c + g ra ca pa nb pb = f ra ca pa 0 0 nullPtr nb pb 0 0 nullPtr + + +-- | Singular values and all right singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'N\' and jobvt == \'A\'. +rightSVR :: Matrix Double -> (Vector Double, Matrix Double) +rightSVR = rightSVAux dgesvd "rightSVR" . fmat + +-- | Singular values and all right singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'N\' and jobvt == \'A\'. +rightSVC :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) +rightSVC = rightSVAux zgesvd "rightSVC" . fmat + +rightSVAux f st x = unsafePerformIO $ do + s <- createVector q + v <- createMatrix ColumnMajor c c + app3 g mat x vec s mat v st + return (s,trans v) + where r = rows x + c = cols x + q = min r c + g ra ca pa = f ra ca pa 0 0 nullPtr + + +-- | Singular values and all left singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'A\' and jobvt == \'N\'. +leftSVR :: Matrix Double -> (Matrix Double, Vector Double) +leftSVR = leftSVAux dgesvd "leftSVR" . fmat + +-- | Singular values and all left singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'A\' and jobvt == \'N\'. +leftSVC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double) +leftSVC = leftSVAux zgesvd "leftSVC" . fmat + +leftSVAux f st x = unsafePerformIO $ do + u <- createMatrix ColumnMajor r r + s <- createVector q + app3 g mat x mat u vec s st + return (u,s) + where r = rows x + c = cols x + q = min r c + g ra ca pa ru cu pu nb pb = f ra ca pa ru cu pu nb pb 0 0 nullPtr + +----------------------------------------------------------------------------- + +foreign import ccall unsafe "eig_l_R" dgeev :: TMMCVM +foreign import ccall unsafe "eig_l_C" zgeev :: TCMCMCVCM +foreign import ccall unsafe "eig_l_S" dsyev :: CInt -> TMVM +foreign import ccall unsafe "eig_l_H" zheev :: CInt -> TCMVCM + +eigAux f st m = unsafePerformIO $ do + l <- createVector r + v <- createMatrix ColumnMajor r r + app3 g mat m vec l mat v st + return (l,v) + where r = rows m + g ra ca pa = f ra ca pa 0 0 nullPtr + + +-- | Eigenvalues and right eigenvectors of a general complex matrix, using LAPACK's /zgeev/. +-- The eigenvectors are the columns of v. The eigenvalues are not sorted. +eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double)) +eigC = eigAux zgeev "eigC" . fmat + +eigOnlyAux f st m = unsafePerformIO $ do + l <- createVector r + app2 g mat m vec l st + return l + where r = rows m + g ra ca pa nl pl = f ra ca pa 0 0 nullPtr nl pl 0 0 nullPtr + +-- | Eigenvalues of a general complex matrix, using LAPACK's /zgeev/ with jobz == \'N\'. +-- The eigenvalues are not sorted. +eigOnlyC :: Matrix (Complex Double) -> Vector (Complex Double) +eigOnlyC = eigOnlyAux zgeev "eigOnlyC" . fmat + +-- | Eigenvalues and right eigenvectors of a general real matrix, using LAPACK's /dgeev/. +-- The eigenvectors are the columns of v. The eigenvalues are not sorted. +eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double)) +eigR m = (s', v'') + where (s,v) = eigRaux (fmat m) + s' = fixeig1 s + v' = toRows $ trans v + v'' = fromColumns $ fixeig (toList s') v' + +eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double) +eigRaux m = unsafePerformIO $ do + l <- createVector r + v <- createMatrix ColumnMajor r r + app3 g mat m vec l mat v "eigR" + return (l,v) + where r = rows m + g ra ca pa = dgeev ra ca pa 0 0 nullPtr + +fixeig1 s = toComplex' (subVector 0 r (asReal s), subVector r r (asReal s)) + where r = dim s + +fixeig [] _ = [] +fixeig [_] [v] = [comp' v] +fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs) + | r1 == r2 && i1 == (-i2) = toComplex' (v1,v2) : toComplex' (v1,scale (-1) v2) : fixeig r vs + | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs) + where scale = vectorMapValR Scale +fixeig _ _ = error "fixeig with impossible inputs" + + +-- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'. +-- The eigenvalues are not sorted. +eigOnlyR :: Matrix Double -> Vector (Complex Double) +eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR" . fmat + + +----------------------------------------------------------------------------- + +eigSHAux f st m = unsafePerformIO $ do + l <- createVector r + v <- createMatrix ColumnMajor r r + app3 f mat m vec l mat v st + return (l,v) + where r = rows m + +-- | Eigenvalues and right eigenvectors of a symmetric real matrix, using LAPACK's /dsyev/. +-- The eigenvectors are the columns of v. +-- The eigenvalues are sorted in descending order (use 'eigS'' for ascending order). +eigS :: Matrix Double -> (Vector Double, Matrix Double) +eigS m = (s', fliprl v) + where (s,v) = eigS' (fmat m) + s' = fromList . reverse . toList $ s + +-- | 'eigS' in ascending order +eigS' :: Matrix Double -> (Vector Double, Matrix Double) +eigS' = eigSHAux (dsyev 1) "eigS'" . fmat + +-- | Eigenvalues and right eigenvectors of a hermitian complex matrix, using LAPACK's /zheev/. +-- The eigenvectors are the columns of v. +-- The eigenvalues are sorted in descending order (use 'eigH'' for ascending order). +eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) +eigH m = (s', fliprl v) + where (s,v) = eigH' (fmat m) + s' = fromList . reverse . toList $ s + +-- | 'eigH' in ascending order +eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) +eigH' = eigSHAux (zheev 1) "eigH'" . fmat + + +-- | Eigenvalues of a symmetric real matrix, using LAPACK's /dsyev/ with jobz == \'N\'. +-- The eigenvalues are sorted in descending order. +eigOnlyS :: Matrix Double -> Vector Double +eigOnlyS = vrev . fst. eigSHAux (dsyev 0) "eigS'" . fmat + +-- | Eigenvalues of a hermitian complex matrix, using LAPACK's /zheev/ with jobz == \'N\'. +-- The eigenvalues are sorted in descending order. +eigOnlyH :: Matrix (Complex Double) -> Vector Double +eigOnlyH = vrev . fst. eigSHAux (zheev 1) "eigH'" . fmat + +vrev = flatten . flipud . reshape 1 + +----------------------------------------------------------------------------- +foreign import ccall unsafe "linearSolveR_l" dgesv :: TMMM +foreign import ccall unsafe "linearSolveC_l" zgesv :: TCMCMCM +foreign import ccall unsafe "cholSolveR_l" dpotrs :: TMMM +foreign import ccall unsafe "cholSolveC_l" zpotrs :: TCMCMCM + +linearSolveSQAux f st a b + | n1==n2 && n1==r = unsafePerformIO $ do + s <- createMatrix ColumnMajor r c + app3 f mat a mat b mat s st + return s + | otherwise = error $ st ++ " of nonsquare matrix" + where n1 = rows a + n2 = cols a + r = rows b + c = cols b + +-- | Solve a real linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /dgesv/. For underconstrained or overconstrained systems use 'linearSolveLSR' or 'linearSolveSVDR'. See also 'lusR'. +linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double +linearSolveR a b = linearSolveSQAux dgesv "linearSolveR" (fmat a) (fmat b) + +-- | Solve a complex linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /zgesv/. For underconstrained or overconstrained systems use 'linearSolveLSC' or 'linearSolveSVDC'. See also 'lusC'. +linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) +linearSolveC a b = linearSolveSQAux zgesv "linearSolveC" (fmat a) (fmat b) + + +-- | Solves a symmetric positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholS'. +cholSolveR :: Matrix Double -> Matrix Double -> Matrix Double +cholSolveR a b = linearSolveSQAux dpotrs "cholSolveR" (fmat a) (fmat b) + +-- | Solves a Hermitian positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholH'. +cholSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) +cholSolveC a b = linearSolveSQAux zpotrs "cholSolveC" (fmat a) (fmat b) + +----------------------------------------------------------------------------------- +foreign import ccall unsafe "linearSolveLSR_l" dgels :: TMMM +foreign import ccall unsafe "linearSolveLSC_l" zgels :: TCMCMCM +foreign import ccall unsafe "linearSolveSVDR_l" dgelss :: Double -> TMMM +foreign import ccall unsafe "linearSolveSVDC_l" zgelss :: Double -> TCMCMCM + +linearSolveAux f st a b = unsafePerformIO $ do + r <- createMatrix ColumnMajor (max m n) nrhs + app3 f mat a mat b mat r st + return r + where m = rows a + n = cols a + nrhs = cols b + +-- | Least squared error solution of an overconstrained real linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /dgels/. For rank-deficient systems use 'linearSolveSVDR'. +linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double +linearSolveLSR a b = subMatrix (0,0) (cols a, cols b) $ + linearSolveAux dgels "linearSolverLSR" (fmat a) (fmat b) + +-- | Least squared error solution of an overconstrained complex linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /zgels/. For rank-deficient systems use 'linearSolveSVDC'. +linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) +linearSolveLSC a b = subMatrix (0,0) (cols a, cols b) $ + linearSolveAux zgels "linearSolveLSC" (fmat a) (fmat b) + +-- | Minimum norm solution of a general real linear least squares problem Ax=B using the SVD, based on LAPACK's /dgelss/. Admits rank-deficient systems but it is slower than 'linearSolveLSR'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used. +linearSolveSVDR :: Maybe Double -- ^ rcond + -> Matrix Double -- ^ coefficient matrix + -> Matrix Double -- ^ right hand sides (as columns) + -> Matrix Double -- ^ solution vectors (as columns) +linearSolveSVDR (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ + linearSolveAux (dgelss rcond) "linearSolveSVDR" (fmat a) (fmat b) +linearSolveSVDR Nothing a b = linearSolveSVDR (Just (-1)) (fmat a) (fmat b) + +-- | Minimum norm solution of a general complex linear least squares problem Ax=B using the SVD, based on LAPACK's /zgelss/. Admits rank-deficient systems but it is slower than 'linearSolveLSC'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used. +linearSolveSVDC :: Maybe Double -- ^ rcond + -> Matrix (Complex Double) -- ^ coefficient matrix + -> Matrix (Complex Double) -- ^ right hand sides (as columns) + -> Matrix (Complex Double) -- ^ solution vectors (as columns) +linearSolveSVDC (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ + linearSolveAux (zgelss rcond) "linearSolveSVDC" (fmat a) (fmat b) +linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) (fmat a) (fmat b) + +----------------------------------------------------------------------------------- +foreign import ccall unsafe "chol_l_H" zpotrf :: TCMCM +foreign import ccall unsafe "chol_l_S" dpotrf :: TMM + +cholAux f st a = do + r <- createMatrix ColumnMajor n n + app2 f mat a mat r st + return r + where n = rows a + +-- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/. +cholH :: Matrix (Complex Double) -> Matrix (Complex Double) +cholH = unsafePerformIO . cholAux zpotrf "cholH" . fmat + +-- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/. +cholS :: Matrix Double -> Matrix Double +cholS = unsafePerformIO . cholAux dpotrf "cholS" . fmat + +-- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/ ('Maybe' version). +mbCholH :: Matrix (Complex Double) -> Maybe (Matrix (Complex Double)) +mbCholH = unsafePerformIO . mbCatch . cholAux zpotrf "cholH" . fmat + +-- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/ ('Maybe' version). +mbCholS :: Matrix Double -> Maybe (Matrix Double) +mbCholS = unsafePerformIO . mbCatch . cholAux dpotrf "cholS" . fmat + +----------------------------------------------------------------------------------- +foreign import ccall unsafe "qr_l_R" dgeqr2 :: TMVM +foreign import ccall unsafe "qr_l_C" zgeqr2 :: TCMCVCM + +-- | QR factorization of a real matrix, using LAPACK's /dgeqr2/. +qrR :: Matrix Double -> (Matrix Double, Vector Double) +qrR = qrAux dgeqr2 "qrR" . fmat + +-- | QR factorization of a complex matrix, using LAPACK's /zgeqr2/. +qrC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double)) +qrC = qrAux zgeqr2 "qrC" . fmat + +qrAux f st a = unsafePerformIO $ do + r <- createMatrix ColumnMajor m n + tau <- createVector mn + app3 f mat a vec tau mat r st + return (r,tau) + where + m = rows a + n = cols a + mn = min m n + +foreign import ccall unsafe "c_dorgqr" dorgqr :: TMVM +foreign import ccall unsafe "c_zungqr" zungqr :: TCMCVCM + +-- | build rotation from reflectors +qrgrR :: Int -> (Matrix Double, Vector Double) -> Matrix Double +qrgrR = qrgrAux dorgqr "qrgrR" +-- | build rotation from reflectors +qrgrC :: Int -> (Matrix (Complex Double), Vector (Complex Double)) -> Matrix (Complex Double) +qrgrC = qrgrAux zungqr "qrgrC" + +qrgrAux f st n (a, tau) = unsafePerformIO $ do + res <- createMatrix ColumnMajor (rows a) n + app3 f mat (fmat a) vec (subVector 0 n tau') mat res st + return res + where + tau' = vjoin [tau, constantD 0 n] + +----------------------------------------------------------------------------------- +foreign import ccall unsafe "hess_l_R" dgehrd :: TMVM +foreign import ccall unsafe "hess_l_C" zgehrd :: TCMCVCM + +-- | Hessenberg factorization of a square real matrix, using LAPACK's /dgehrd/. +hessR :: Matrix Double -> (Matrix Double, Vector Double) +hessR = hessAux dgehrd "hessR" . fmat + +-- | Hessenberg factorization of a square complex matrix, using LAPACK's /zgehrd/. +hessC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double)) +hessC = hessAux zgehrd "hessC" . fmat + +hessAux f st a = unsafePerformIO $ do + r <- createMatrix ColumnMajor m n + tau <- createVector (mn-1) + app3 f mat a vec tau mat r st + return (r,tau) + where m = rows a + n = cols a + mn = min m n + +----------------------------------------------------------------------------------- +foreign import ccall unsafe "schur_l_R" dgees :: TMMM +foreign import ccall unsafe "schur_l_C" zgees :: TCMCMCM + +-- | Schur factorization of a square real matrix, using LAPACK's /dgees/. +schurR :: Matrix Double -> (Matrix Double, Matrix Double) +schurR = schurAux dgees "schurR" . fmat + +-- | Schur factorization of a square complex matrix, using LAPACK's /zgees/. +schurC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double)) +schurC = schurAux zgees "schurC" . fmat + +schurAux f st a = unsafePerformIO $ do + u <- createMatrix ColumnMajor n n + s <- createMatrix ColumnMajor n n + app3 f mat a mat u mat s st + return (u,s) + where n = rows a + +----------------------------------------------------------------------------------- +foreign import ccall unsafe "lu_l_R" dgetrf :: TMVM +foreign import ccall unsafe "lu_l_C" zgetrf :: TCMVCM + +-- | LU factorization of a general real matrix, using LAPACK's /dgetrf/. +luR :: Matrix Double -> (Matrix Double, [Int]) +luR = luAux dgetrf "luR" . fmat + +-- | LU factorization of a general complex matrix, using LAPACK's /zgetrf/. +luC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int]) +luC = luAux zgetrf "luC" . fmat + +luAux f st a = unsafePerformIO $ do + lu <- createMatrix ColumnMajor n m + piv <- createVector (min n m) + app3 f mat a vec piv mat lu st + return (lu, map (pred.round) (toList piv)) + where n = rows a + m = cols a + +----------------------------------------------------------------------------------- +type TW a = CInt -> PD -> a +type TQ a = CInt -> CInt -> PC -> a + +foreign import ccall unsafe "luS_l_R" dgetrs :: TMVMM +foreign import ccall unsafe "luS_l_C" zgetrs :: TQ (TW (TQ (TQ (IO CInt)))) + +-- | Solve a real linear system from a precomputed LU decomposition ('luR'), using LAPACK's /dgetrs/. +lusR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double +lusR a piv b = lusAux dgetrs "lusR" (fmat a) piv (fmat b) + +-- | Solve a real linear system from a precomputed LU decomposition ('luC'), using LAPACK's /zgetrs/. +lusC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double) +lusC a piv b = lusAux zgetrs "lusC" (fmat a) piv (fmat b) + +lusAux f st a piv b + | n1==n2 && n2==n =unsafePerformIO $ do + x <- createMatrix ColumnMajor n m + app4 f mat a vec piv' mat b mat x st + return x + | otherwise = error $ st ++ " on LU factorization of nonsquare matrix" + where n1 = rows a + n2 = cols a + n = rows b + m = cols b + piv' = fromList (map (fromIntegral.succ) piv) :: Vector Double + diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.c new file mode 100644 index 0000000..e5e45ef --- /dev/null +++ b/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.c @@ -0,0 +1,1489 @@ +#include +#include +#include +#include +#include +#include "lapack-aux.h" + +#define MACRO(B) do {B} while (0) +#define ERROR(CODE) MACRO(return CODE;) +#define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);}) + +#define MIN(A,B) ((A)<(B)?(A):(B)) +#define MAX(A,B) ((A)>(B)?(A):(B)) + +// #define DBGL + +#ifdef DBGL +#define DEBUGMSG(M) printf("\nLAPACK "M"\n"); +#else +#define DEBUGMSG(M) +#endif + +#define OK return 0; + +// #ifdef DBGL +// #define DEBUGMSG(M) printf("LAPACK Wrapper "M"\n: "); size_t t0 = time(NULL); +// #define OK MACRO(printf("%ld s\n",time(0)-t0); return 0;); +// #else +// #define DEBUGMSG(M) +// #define OK return 0; +// #endif + +#define TRACEMAT(M) {int q; printf(" %d x %d: ",M##r,M##c); \ + for(q=0;q=1 && ar==ac && ar==br,BAD_SIZE); + DEBUGMSG("linearSolveR_l"); + double*AC = (double*)malloc(n*n*sizeof(double)); + memcpy(AC,ap,n*n*sizeof(double)); + memcpy(xp,bp,n*nhrs*sizeof(double)); + integer * ipiv = (integer*)malloc(n*sizeof(integer)); + integer res; + dgesv_ (&n,&nhrs, + AC, &n, + ipiv, + xp, &n, + &res); + if(res>0) { + return SINGULAR; + } + CHECK(res,res); + free(ipiv); + free(AC); + OK +} + +//////////////////// general complex linear system //////////// + +/* Subroutine */ int zgesv_(integer *n, integer *nrhs, doublecomplex *a, + integer *lda, integer *ipiv, doublecomplex *b, integer *ldb, integer * + info); + +int linearSolveC_l(KCMAT(a),KCMAT(b),CMAT(x)) { + integer n = ar; + integer nhrs = bc; + REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); + DEBUGMSG("linearSolveC_l"); + doublecomplex*AC = (doublecomplex*)malloc(n*n*sizeof(doublecomplex)); + memcpy(AC,ap,n*n*sizeof(doublecomplex)); + memcpy(xp,bp,n*nhrs*sizeof(doublecomplex)); + integer * ipiv = (integer*)malloc(n*sizeof(integer)); + integer res; + zgesv_ (&n,&nhrs, + AC, &n, + ipiv, + xp, &n, + &res); + if(res>0) { + return SINGULAR; + } + CHECK(res,res); + free(ipiv); + free(AC); + OK +} + +//////// symmetric positive definite real linear system using Cholesky //////////// + +/* Subroutine */ int dpotrs_(char *uplo, integer *n, integer *nrhs, + doublereal *a, integer *lda, doublereal *b, integer *ldb, integer * + info); + +int cholSolveR_l(KDMAT(a),KDMAT(b),DMAT(x)) { + integer n = ar; + integer nhrs = bc; + REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); + DEBUGMSG("cholSolveR_l"); + memcpy(xp,bp,n*nhrs*sizeof(double)); + integer res; + dpotrs_ ("U", + &n,&nhrs, + (double*)ap, &n, + xp, &n, + &res); + CHECK(res,res); + OK +} + +//////// Hermitian positive definite real linear system using Cholesky //////////// + +/* Subroutine */ int zpotrs_(char *uplo, integer *n, integer *nrhs, + doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, + integer *info); + +int cholSolveC_l(KCMAT(a),KCMAT(b),CMAT(x)) { + integer n = ar; + integer nhrs = bc; + REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); + DEBUGMSG("cholSolveC_l"); + memcpy(xp,bp,n*nhrs*sizeof(doublecomplex)); + integer res; + zpotrs_ ("U", + &n,&nhrs, + (doublecomplex*)ap, &n, + xp, &n, + &res); + CHECK(res,res); + OK +} + +//////////////////// least squares real linear system //////////// + +/* Subroutine */ int dgels_(char *trans, integer *m, integer *n, integer * + nrhs, doublereal *a, integer *lda, doublereal *b, integer *ldb, + doublereal *work, integer *lwork, integer *info); + +int linearSolveLSR_l(KDMAT(a),KDMAT(b),DMAT(x)) { + integer m = ar; + integer n = ac; + integer nrhs = bc; + integer ldb = xr; + REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); + DEBUGMSG("linearSolveLSR_l"); + double*AC = (double*)malloc(m*n*sizeof(double)); + memcpy(AC,ap,m*n*sizeof(double)); + if (m>=n) { + memcpy(xp,bp,m*nrhs*sizeof(double)); + } else { + int k; + for(k = 0; k0) { + return SINGULAR; + } + CHECK(res,res); + free(work); + free(AC); + OK +} + +//////////////////// least squares complex linear system //////////// + +/* Subroutine */ int zgels_(char *trans, integer *m, integer *n, integer * + nrhs, doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, + doublecomplex *work, integer *lwork, integer *info); + +int linearSolveLSC_l(KCMAT(a),KCMAT(b),CMAT(x)) { + integer m = ar; + integer n = ac; + integer nrhs = bc; + integer ldb = xr; + REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); + DEBUGMSG("linearSolveLSC_l"); + doublecomplex*AC = (doublecomplex*)malloc(m*n*sizeof(doublecomplex)); + memcpy(AC,ap,m*n*sizeof(doublecomplex)); + if (m>=n) { + memcpy(xp,bp,m*nrhs*sizeof(doublecomplex)); + } else { + int k; + for(k = 0; k0) { + return SINGULAR; + } + CHECK(res,res); + free(work); + free(AC); + OK +} + +//////////////////// least squares real linear system using SVD //////////// + +/* Subroutine */ int dgelss_(integer *m, integer *n, integer *nrhs, + doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal * + s, doublereal *rcond, integer *rank, doublereal *work, integer *lwork, + integer *info); + +int linearSolveSVDR_l(double rcond,KDMAT(a),KDMAT(b),DMAT(x)) { + integer m = ar; + integer n = ac; + integer nrhs = bc; + integer ldb = xr; + REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); + DEBUGMSG("linearSolveSVDR_l"); + double*AC = (double*)malloc(m*n*sizeof(double)); + double*S = (double*)malloc(MIN(m,n)*sizeof(double)); + memcpy(AC,ap,m*n*sizeof(double)); + if (m>=n) { + memcpy(xp,bp,m*nrhs*sizeof(double)); + } else { + int k; + for(k = 0; k0) { + return NOCONVER; + } + CHECK(res,res); + free(work); + free(S); + free(AC); + OK +} + +//////////////////// least squares complex linear system using SVD //////////// + +// not in clapack.h + +int zgelss_(integer *m, integer *n, integer *nhrs, + doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, doublereal *s, + doublereal *rcond, integer* rank, + doublecomplex *work, integer* lwork, doublereal* rwork, + integer *info); + +int linearSolveSVDC_l(double rcond, KCMAT(a),KCMAT(b),CMAT(x)) { + integer m = ar; + integer n = ac; + integer nrhs = bc; + integer ldb = xr; + REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); + DEBUGMSG("linearSolveSVDC_l"); + doublecomplex*AC = (doublecomplex*)malloc(m*n*sizeof(doublecomplex)); + double*S = (double*)malloc(MIN(m,n)*sizeof(double)); + double*RWORK = (double*)malloc(5*MIN(m,n)*sizeof(double)); + memcpy(AC,ap,m*n*sizeof(doublecomplex)); + if (m>=n) { + memcpy(xp,bp,m*nrhs*sizeof(doublecomplex)); + } else { + int k; + for(k = 0; k0) { + return NOCONVER; + } + CHECK(res,res); + free(work); + free(RWORK); + free(S); + free(AC); + OK +} + +//////////////////// Cholesky factorization ///////////////////////// + +/* Subroutine */ int zpotrf_(char *uplo, integer *n, doublecomplex *a, + integer *lda, integer *info); + +int chol_l_H(KCMAT(a),CMAT(l)) { + integer n = ar; + REQUIRES(n>=1 && ac == n && lr==n && lc==n,BAD_SIZE); + DEBUGMSG("chol_l_H"); + memcpy(lp,ap,n*n*sizeof(doublecomplex)); + char uplo = 'U'; + integer res; + zpotrf_ (&uplo,&n,lp,&n,&res); + CHECK(res>0,NODEFPOS); + CHECK(res,res); + doublecomplex zero = {0.,0.}; + int r,c; + for (r=0; r=1 && ac == n && lr==n && lc==n,BAD_SIZE); + DEBUGMSG("chol_l_S"); + memcpy(lp,ap,n*n*sizeof(double)); + char uplo = 'U'; + integer res; + dpotrf_ (&uplo,&n,lp,&n,&res); + CHECK(res>0,NODEFPOS); + CHECK(res,res); + int r,c; + for (r=0; r=1 && n >=1 && rr== m && rc == n && taun == mn, BAD_SIZE); + DEBUGMSG("qr_l_R"); + double *WORK = (double*)malloc(n*sizeof(double)); + CHECK(!WORK,MEM); + memcpy(rp,ap,m*n*sizeof(double)); + integer res; + dgeqr2_ (&m,&n,rp,&m,taup,WORK,&res); + CHECK(res,res); + free(WORK); + OK +} + +/* Subroutine */ int zgeqr2_(integer *m, integer *n, doublecomplex *a, + integer *lda, doublecomplex *tau, doublecomplex *work, integer *info); + +int qr_l_C(KCMAT(a), CVEC(tau), CMAT(r)) { + integer m = ar; + integer n = ac; + integer mn = MIN(m,n); + REQUIRES(m>=1 && n >=1 && rr== m && rc == n && taun == mn, BAD_SIZE); + DEBUGMSG("qr_l_C"); + doublecomplex *WORK = (doublecomplex*)malloc(n*sizeof(doublecomplex)); + CHECK(!WORK,MEM); + memcpy(rp,ap,m*n*sizeof(doublecomplex)); + integer res; + zgeqr2_ (&m,&n,rp,&m,taup,WORK,&res); + CHECK(res,res); + free(WORK); + OK +} + +/* Subroutine */ int dorgqr_(integer *m, integer *n, integer *k, doublereal * + a, integer *lda, doublereal *tau, doublereal *work, integer *lwork, + integer *info); + +int c_dorgqr(KDMAT(a), KDVEC(tau), DMAT(r)) { + integer m = ar; + integer n = MIN(ac,ar); + integer k = taun; + DEBUGMSG("c_dorgqr"); + integer lwork = 8*n; // FIXME + double *WORK = (double*)malloc(lwork*sizeof(double)); + CHECK(!WORK,MEM); + memcpy(rp,ap,m*k*sizeof(double)); + integer res; + dorgqr_ (&m,&n,&k,rp,&m,(double*)taup,WORK,&lwork,&res); + CHECK(res,res); + free(WORK); + OK +} + +/* Subroutine */ int zungqr_(integer *m, integer *n, integer *k, + doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex * + work, integer *lwork, integer *info); + +int c_zungqr(KCMAT(a), KCVEC(tau), CMAT(r)) { + integer m = ar; + integer n = MIN(ac,ar); + integer k = taun; + DEBUGMSG("z_ungqr"); + integer lwork = 8*n; // FIXME + doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); + CHECK(!WORK,MEM); + memcpy(rp,ap,m*k*sizeof(doublecomplex)); + integer res; + zungqr_ (&m,&n,&k,rp,&m,(doublecomplex*)taup,WORK,&lwork,&res); + CHECK(res,res); + free(WORK); + OK +} + + +//////////////////// Hessenberg factorization ///////////////////////// + +/* Subroutine */ int dgehrd_(integer *n, integer *ilo, integer *ihi, + doublereal *a, integer *lda, doublereal *tau, doublereal *work, + integer *lwork, integer *info); + +int hess_l_R(KDMAT(a), DVEC(tau), DMAT(r)) { + integer m = ar; + integer n = ac; + integer mn = MIN(m,n); + REQUIRES(m>=1 && n == m && rr== m && rc == n && taun == mn-1, BAD_SIZE); + DEBUGMSG("hess_l_R"); + integer lwork = 5*n; // fixme + double *WORK = (double*)malloc(lwork*sizeof(double)); + CHECK(!WORK,MEM); + memcpy(rp,ap,m*n*sizeof(double)); + integer res; + integer one = 1; + dgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res); + CHECK(res,res); + free(WORK); + OK +} + + +/* Subroutine */ int zgehrd_(integer *n, integer *ilo, integer *ihi, + doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex * + work, integer *lwork, integer *info); + +int hess_l_C(KCMAT(a), CVEC(tau), CMAT(r)) { + integer m = ar; + integer n = ac; + integer mn = MIN(m,n); + REQUIRES(m>=1 && n == m && rr== m && rc == n && taun == mn-1, BAD_SIZE); + DEBUGMSG("hess_l_C"); + integer lwork = 5*n; // fixme + doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); + CHECK(!WORK,MEM); + memcpy(rp,ap,m*n*sizeof(doublecomplex)); + integer res; + integer one = 1; + zgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res); + CHECK(res,res); + free(WORK); + OK +} + +//////////////////// Schur factorization ///////////////////////// + +/* Subroutine */ int dgees_(char *jobvs, char *sort, L_fp select, integer *n, + doublereal *a, integer *lda, integer *sdim, doublereal *wr, + doublereal *wi, doublereal *vs, integer *ldvs, doublereal *work, + integer *lwork, logical *bwork, integer *info); + +int schur_l_R(KDMAT(a), DMAT(u), DMAT(s)) { + integer m = ar; + integer n = ac; + REQUIRES(m>=1 && n==m && ur==n && uc==n && sr==n && sc==n, BAD_SIZE); + DEBUGMSG("schur_l_R"); + //int k; + //printf("---------------------------\n"); + //printf("%p: ",ap); for(k=0;k0) { + return NOCONVER; + } + CHECK(res,res); + free(WR); + free(WI); + free(BWORK); + free(WORK); + OK +} + + +/* Subroutine */ int zgees_(char *jobvs, char *sort, L_fp select, integer *n, + doublecomplex *a, integer *lda, integer *sdim, doublecomplex *w, + doublecomplex *vs, integer *ldvs, doublecomplex *work, integer *lwork, + doublereal *rwork, logical *bwork, integer *info); + +int schur_l_C(KCMAT(a), CMAT(u), CMAT(s)) { + integer m = ar; + integer n = ac; + REQUIRES(m>=1 && n==m && ur==n && uc==n && sr==n && sc==n, BAD_SIZE); + DEBUGMSG("schur_l_C"); + memcpy(sp,ap,n*n*sizeof(doublecomplex)); + integer lwork = 6*n; // fixme + doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); + doublecomplex *W = (doublecomplex*)malloc(n*sizeof(doublecomplex)); + // W not really required in this call + logical *BWORK = (logical*)malloc(n*sizeof(logical)); + double *RWORK = (double*)malloc(n*sizeof(double)); + integer res; + integer sdim; + zgees_ ("V","N",NULL,&n,sp,&n,&sdim,W, + up,&n, + WORK,&lwork,RWORK,BWORK,&res); + if(res>0) { + return NOCONVER; + } + CHECK(res,res); + free(W); + free(BWORK); + free(WORK); + OK +} + +//////////////////// LU factorization ///////////////////////// + +/* Subroutine */ int dgetrf_(integer *m, integer *n, doublereal *a, integer * + lda, integer *ipiv, integer *info); + +int lu_l_R(KDMAT(a), DVEC(ipiv), DMAT(r)) { + integer m = ar; + integer n = ac; + integer mn = MIN(m,n); + REQUIRES(m>=1 && n >=1 && ipivn == mn, BAD_SIZE); + DEBUGMSG("lu_l_R"); + integer* auxipiv = (integer*)malloc(mn*sizeof(integer)); + memcpy(rp,ap,m*n*sizeof(double)); + integer res; + dgetrf_ (&m,&n,rp,&m,auxipiv,&res); + if(res>0) { + res = 0; // fixme + } + CHECK(res,res); + int k; + for (k=0; k=1 && n >=1 && ipivn == mn, BAD_SIZE); + DEBUGMSG("lu_l_C"); + integer* auxipiv = (integer*)malloc(mn*sizeof(integer)); + memcpy(rp,ap,m*n*sizeof(doublecomplex)); + integer res; + zgetrf_ (&m,&n,rp,&m,auxipiv,&res); + if(res>0) { + res = 0; // fixme + } + CHECK(res,res); + int k; + for (k=0; k0; + } + OK +} + +int stepD(DVEC(x),DVEC(y)) { + DEBUGMSG("stepD") + int k; + for(k=0;k0; + } + OK +} + +//////////////////// cond ///////////////////////// + +int condF(FVEC(x),FVEC(y),FVEC(lt),FVEC(eq),FVEC(gt),FVEC(r)) { + REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE); + DEBUGMSG("condF") + int k; + for(k=0;kyp[k]?gtp[k]:eqp[k]); + } + OK +} + +int condD(DVEC(x),DVEC(y),DVEC(lt),DVEC(eq),DVEC(gt),DVEC(r)) { + REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE); + DEBUGMSG("condD") + int k; + for(k=0;kyp[k]?gtp[k]:eqp[k]); + } + OK +} + diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.h b/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.h new file mode 100644 index 0000000..a3f1899 --- /dev/null +++ b/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.h @@ -0,0 +1,60 @@ +/* + * We have copied the definitions in f2c.h required + * to compile clapack.h, modified to support both + * 32 and 64 bit + + http://opengrok.creo.hu/dragonfly/xref/src/contrib/gcc-3.4/libf2c/readme.netlib + http://www.ibm.com/developerworks/library/l-port64.html + */ + +#ifdef _LP64 +typedef int integer; +typedef unsigned int uinteger; +typedef int logical; +typedef long longint; /* system-dependent */ +typedef unsigned long ulongint; /* system-dependent */ +#else +typedef long int integer; +typedef unsigned long int uinteger; +typedef long int logical; +typedef long long longint; /* system-dependent */ +typedef unsigned long long ulongint; /* system-dependent */ +#endif + +typedef char *address; +typedef short int shortint; +typedef float real; +typedef double doublereal; +typedef struct { real r, i; } complex; +typedef struct { doublereal r, i; } doublecomplex; +typedef short int shortlogical; +typedef char logical1; +typedef char integer1; + +typedef logical (*L_fp)(); +typedef short ftnlen; + +/********************************************************/ + +#define FVEC(A) int A##n, float*A##p +#define DVEC(A) int A##n, double*A##p +#define QVEC(A) int A##n, complex*A##p +#define CVEC(A) int A##n, doublecomplex*A##p +#define PVEC(A) int A##n, void* A##p, int A##s +#define FMAT(A) int A##r, int A##c, float* A##p +#define DMAT(A) int A##r, int A##c, double* A##p +#define QMAT(A) int A##r, int A##c, complex* A##p +#define CMAT(A) int A##r, int A##c, doublecomplex* A##p +#define PMAT(A) int A##r, int A##c, void* A##p, int A##s + +#define KFVEC(A) int A##n, const float*A##p +#define KDVEC(A) int A##n, const double*A##p +#define KQVEC(A) int A##n, const complex*A##p +#define KCVEC(A) int A##n, const doublecomplex*A##p +#define KPVEC(A) int A##n, const void* A##p, int A##s +#define KFMAT(A) int A##r, int A##c, const float* A##p +#define KDMAT(A) int A##r, int A##c, const double* A##p +#define KQMAT(A) int A##r, int A##c, const complex* A##p +#define KCMAT(A) int A##r, int A##c, const doublecomplex* A##p +#define KPMAT(A) int A##r, int A##c, const void* A##p, int A##s + diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs new file mode 100644 index 0000000..7d134bf --- /dev/null +++ b/packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs @@ -0,0 +1,295 @@ +{-# LANGUAGE FlexibleContexts #-} +----------------------------------------------------------------------------- +{- | +Module : Numeric.LinearAlgebra.Util +Copyright : (c) Alberto Ruiz 2013 +License : GPL + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional + +-} +----------------------------------------------------------------------------- +{-# OPTIONS_HADDOCK hide #-} + +module Numeric.LinearAlgebra.Util( + + -- * Convenience functions + size, disp, + zeros, ones, + diagl, + row, + col, + (&), (¦), (——), (#), + (?), (¿), + rand, randn, + cross, + norm, + unitary, + mt, + pairwiseD2, + rowOuters, + null1, + null1sym, + -- * Convolution + -- ** 1D + corr, conv, corrMin, + -- ** 2D + corr2, conv2, separable, + -- * Tools for the Kronecker product + -- + -- | (see A. Fusiello, A matter of notation: Several uses of the Kronecker product in + -- 3d computer vision, Pattern Recognition Letters 28 (15) (2007) 2127-2132) + + -- + -- | @`vec` (a \<> x \<> b) == ('trans' b ` 'kronecker' ` a) \<> 'vec' x@ + vec, + vech, + dup, + vtrans, + -- * Plot + mplot, + plot, parametricPlot, + splot, mesh, meshdom, + matrixToPGM, imshow, + gnuplotX, gnuplotpdf, gnuplotWin +) where + +import Numeric.Container +import Numeric.LinearAlgebra.Algorithms hiding (i) +import Numeric.Matrix() +import Numeric.Vector() + +import System.Random(randomIO) +import Numeric.LinearAlgebra.Util.Convolution +import Graphics.Plot + + +{- | print a real matrix with given number of digits after the decimal point + +>>> disp 5 $ ident 2 / 3 +2x2 +0.33333 0.00000 +0.00000 0.33333 + +-} +disp :: Int -> Matrix Double -> IO () + +disp n = putStrLn . dispf n + +-- | pseudorandom matrix with uniform elements between 0 and 1 +randm :: RandDist + -> Int -- ^ rows + -> Int -- ^ columns + -> IO (Matrix Double) +randm d r c = do + seed <- randomIO + return (reshape c $ randomVector seed d (r*c)) + +-- | pseudorandom matrix with uniform elements between 0 and 1 +rand :: Int -> Int -> IO (Matrix Double) +rand = randm Uniform + +{- | pseudorandom matrix with normal elements + +>>> x <- randn 3 5 +>>> disp 3 x +3x5 +0.386 -1.141 0.491 -0.510 1.512 +0.069 -0.919 1.022 -0.181 0.745 +0.313 -0.670 -0.097 -1.575 -0.583 + +-} +randn :: Int -> Int -> IO (Matrix Double) +randn = randm Gaussian + +{- | create a real diagonal matrix from a list + +>>> diagl [1,2,3] +(3><3) + [ 1.0, 0.0, 0.0 + , 0.0, 2.0, 0.0 + , 0.0, 0.0, 3.0 ] + +-} +diagl :: [Double] -> Matrix Double +diagl = diag . fromList + +-- | a real matrix of zeros +zeros :: Int -- ^ rows + -> Int -- ^ columns + -> Matrix Double +zeros r c = konst 0 (r,c) + +-- | a real matrix of ones +ones :: Int -- ^ rows + -> Int -- ^ columns + -> Matrix Double +ones r c = konst 1 (r,c) + +-- | concatenation of real vectors +infixl 3 & +(&) :: Vector Double -> Vector Double -> Vector Double +a & b = vjoin [a,b] + +{- | horizontal concatenation of real matrices + + (unicode 0x00a6, broken bar) + +>>> ident 3 ¦ konst 7 (3,4) +(3><7) + [ 1.0, 0.0, 0.0, 7.0, 7.0, 7.0, 7.0 + , 0.0, 1.0, 0.0, 7.0, 7.0, 7.0, 7.0 + , 0.0, 0.0, 1.0, 7.0, 7.0, 7.0, 7.0 ] + +-} +infixl 3 ¦ +(¦) :: Matrix Double -> Matrix Double -> Matrix Double +a ¦ b = fromBlocks [[a,b]] + +-- | vertical concatenation of real matrices +-- +-- (unicode 0x2014, em dash) +(——) :: Matrix Double -> Matrix Double -> Matrix Double +infixl 2 —— +a —— b = fromBlocks [[a],[b]] + +(#) :: Matrix Double -> Matrix Double -> Matrix Double +infixl 2 # +a # b = fromBlocks [[a],[b]] + +-- | create a single row real matrix from a list +row :: [Double] -> Matrix Double +row = asRow . fromList + +-- | create a single column real matrix from a list +col :: [Double] -> Matrix Double +col = asColumn . fromList + +{- | extract rows + +>>> (20><4) [1..] ? [2,1,1] +(3><4) + [ 9.0, 10.0, 11.0, 12.0 + , 5.0, 6.0, 7.0, 8.0 + , 5.0, 6.0, 7.0, 8.0 ] + +-} +infixl 9 ? +(?) :: Element t => Matrix t -> [Int] -> Matrix t +(?) = flip extractRows + +{- | extract columns + +(unicode 0x00bf, inverted question mark, Alt-Gr ?) + +>>> (3><4) [1..] ¿ [3,0] +(3><2) + [ 4.0, 1.0 + , 8.0, 5.0 + , 12.0, 9.0 ] + +-} +infixl 9 ¿ +(¿) :: Element t => Matrix t -> [Int] -> Matrix t +(¿)= flip extractColumns + + +cross :: Vector Double -> Vector Double -> Vector Double +-- ^ cross product (for three-element real vectors) +cross x y | dim x == 3 && dim y == 3 = fromList [z1,z2,z3] + | otherwise = error $ "cross ("++show x++") ("++show y++")" + where + [x1,x2,x3] = toList x + [y1,y2,y3] = toList y + z1 = x2*y3-x3*y2 + z2 = x3*y1-x1*y3 + z3 = x1*y2-x2*y1 + +norm :: Vector Double -> Double +-- ^ 2-norm of real vector +norm = pnorm PNorm2 + + +-- | Obtains a vector in the same direction with 2-norm=1 +unitary :: Vector Double -> Vector Double +unitary v = v / scalar (norm v) + +-- | ('rows' &&& 'cols') +size :: Matrix t -> (Int, Int) +size m = (rows m, cols m) + +-- | trans . inv +mt :: Matrix Double -> Matrix Double +mt = trans . inv + +---------------------------------------------------------------------- + +-- | Matrix of pairwise squared distances of row vectors +-- (using the matrix product trick in blog.smola.org) +pairwiseD2 :: Matrix Double -> Matrix Double -> Matrix Double +pairwiseD2 x y | ok = x2 `outer` oy + ox `outer` y2 - 2* x <> trans y + | otherwise = error $ "pairwiseD2 with different number of columns: " + ++ show (size x) ++ ", " ++ show (size y) + where + ox = one (rows x) + oy = one (rows y) + oc = one (cols x) + one k = constant 1 k + x2 = x * x <> oc + y2 = y * y <> oc + ok = cols x == cols y + +-------------------------------------------------------------------------------- + +-- | outer products of rows +rowOuters :: Matrix Double -> Matrix Double -> Matrix Double +rowOuters a b = a' * b' + where + a' = kronecker a (ones 1 (cols b)) + b' = kronecker (ones 1 (cols a)) b + +-------------------------------------------------------------------------------- + +-- | solution of overconstrained homogeneous linear system +null1 :: Matrix Double -> Vector Double +null1 = last . toColumns . snd . rightSV + +-- | solution of overconstrained homogeneous symmetric linear system +null1sym :: Matrix Double -> Vector Double +null1sym = last . toColumns . snd . eigSH' + +-------------------------------------------------------------------------------- + +vec :: Element t => Matrix t -> Vector t +-- ^ stacking of columns +vec = flatten . trans + + +vech :: Element t => Matrix t -> Vector t +-- ^ half-vectorization (of the lower triangular part) +vech m = vjoin . zipWith f [0..] . toColumns $ m + where + f k v = subVector k (dim v - k) v + + +dup :: (Num t, Num (Vector t), Element t) => Int -> Matrix t +-- ^ duplication matrix (@'dup' k \<> 'vech' m == 'vec' m@, for symmetric m of 'dim' k) +dup k = trans $ fromRows $ map f es + where + rs = zip [0..] (toRows (ident (k^(2::Int)))) + es = [(i,j) | j <- [0..k-1], i <- [0..k-1], i>=j ] + f (i,j) | i == j = g (k*j + i) + | otherwise = g (k*j + i) + g (k*i + j) + g j = v + where + Just v = lookup j rs + + +vtrans :: Element t => Int -> Matrix t -> Matrix t +-- ^ generalized \"vector\" transposition: @'vtrans' 1 == 'trans'@, and @'vtrans' ('rows' m) m == 'asColumn' ('vec' m)@ +vtrans p m | r == 0 = fromBlocks . map (map asColumn . takesV (replicate q p)) . toColumns $ m + | otherwise = error $ "vtrans " ++ show p ++ " of matrix with " ++ show (rows m) ++ " rows" + where + (q,r) = divMod (rows m) p + diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/Util/Convolution.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/Util/Convolution.hs new file mode 100644 index 0000000..82de476 --- /dev/null +++ b/packages/hmatrix/src/Numeric/LinearAlgebra/Util/Convolution.hs @@ -0,0 +1,114 @@ +{-# LANGUAGE FlexibleContexts #-} +----------------------------------------------------------------------------- +{- | +Module : Numeric.LinearAlgebra.Util.Convolution +Copyright : (c) Alberto Ruiz 2012 +License : GPL + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional + +-} +----------------------------------------------------------------------------- + +module Numeric.LinearAlgebra.Util.Convolution( + corr, conv, corrMin, + corr2, conv2, separable +) where + +import Numeric.LinearAlgebra + + +vectSS :: Element t => Int -> Vector t -> Matrix t +vectSS n v = fromRows [ subVector k n v | k <- [0 .. dim v - n] ] + + +corr :: Product t => Vector t -- ^ kernel + -> Vector t -- ^ source + -> Vector t +{- ^ correlation + +>>> corr (fromList[1,2,3]) (fromList [1..10]) +fromList [14.0,20.0,26.0,32.0,38.0,44.0,50.0,56.0] + +-} +corr ker v | dim ker <= dim v = vectSS (dim ker) v <> ker + | otherwise = error $ "corr: dim kernel ("++show (dim ker)++") > dim vector ("++show (dim v)++")" + + +conv :: (Product t, Num t) => Vector t -> Vector t -> Vector t +{- ^ convolution ('corr' with reversed kernel and padded input, equivalent to polynomial product) + +>>> conv (fromList[1,1]) (fromList [-1,1]) +fromList [-1.0,0.0,1.0] + +-} +conv ker v = corr ker' v' + where + ker' = (flatten.fliprl.asRow) ker + v' | dim ker > 1 = vjoin [z,v,z] + | otherwise = v + z = constant 0 (dim ker -1) + +corrMin :: (Container Vector t, RealElement t, Product t) + => Vector t + -> Vector t + -> Vector t +-- ^ similar to 'corr', using 'min' instead of (*) +corrMin ker v = minEvery ss (asRow ker) <> ones + where + minEvery a b = cond a b a a b + ss = vectSS (dim ker) v + ones = konst 1 (dim ker) + + + +matSS :: Element t => Int -> Matrix t -> [Matrix t] +matSS dr m = map (reshape c) [ subVector (k*c) n v | k <- [0 .. r - dr] ] + where + v = flatten m + c = cols m + r = rows m + n = dr*c + + +corr2 :: Product a => Matrix a -> Matrix a -> Matrix a +-- ^ 2D correlation +corr2 ker mat = dims + . concatMap (map (udot ker' . flatten) . matSS c . trans) + . matSS r $ mat + where + r = rows ker + c = cols ker + ker' = flatten (trans ker) + rr = rows mat - r + 1 + rc = cols mat - c + 1 + dims | rr > 0 && rc > 0 = (rr >< rc) + | otherwise = error $ "corr2: dim kernel ("++sz ker++") > dim matrix ("++sz mat++")" + sz m = show (rows m)++"x"++show (cols m) + +conv2 :: (Num a, Product a, Container Vector a) => Matrix a -> Matrix a -> Matrix a +-- ^ 2D convolution +conv2 k m = corr2 (fliprl . flipud $ k) pm + where + pm | r == 0 && c == 0 = m + | r == 0 = fromBlocks [[z3,m,z3]] + | c == 0 = fromBlocks [[z2],[m],[z2]] + | otherwise = fromBlocks [[z1,z2,z1] + ,[z3, m,z3] + ,[z1,z2,z1]] + r = rows k - 1 + c = cols k - 1 + h = rows m + w = cols m + z1 = konst 0 (r,c) + z2 = konst 0 (r,w) + z3 = konst 0 (h,c) + +-- TODO: could be simplified using future empty arrays + + +separable :: Element t => (Vector t -> Vector t) -> Matrix t -> Matrix t +-- ^ matrix computation implemented as separated vector operations by rows and columns. +separable f = fromColumns . map f . toColumns . fromRows . map f . toRows + diff --git a/packages/hmatrix/src/Numeric/Matrix.hs b/packages/hmatrix/src/Numeric/Matrix.hs new file mode 100644 index 0000000..e285ff2 --- /dev/null +++ b/packages/hmatrix/src/Numeric/Matrix.hs @@ -0,0 +1,98 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} + +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Matrix +-- Copyright : (c) Alberto Ruiz 2010 +-- License : GPL-style +-- +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- Portability : portable +-- +-- Provides instances of standard classes 'Show', 'Read', 'Eq', +-- 'Num', 'Fractional', and 'Floating' for 'Matrix'. +-- +-- In arithmetic operations one-component +-- vectors and matrices automatically expand to match the dimensions of the other operand. + +----------------------------------------------------------------------------- + +module Numeric.Matrix ( + ) where + +------------------------------------------------------------------- + +import Numeric.Container +import qualified Data.Monoid as M +import Data.List(partition) + +------------------------------------------------------------------- + +instance Container Matrix a => Eq (Matrix a) where + (==) = equal + +instance (Container Matrix a, Num (Vector a)) => Num (Matrix a) where + (+) = liftMatrix2Auto (+) + (-) = liftMatrix2Auto (-) + negate = liftMatrix negate + (*) = liftMatrix2Auto (*) + signum = liftMatrix signum + abs = liftMatrix abs + fromInteger = (1><1) . return . fromInteger + +--------------------------------------------------- + +instance (Container Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where + fromRational n = (1><1) [fromRational n] + (/) = liftMatrix2Auto (/) + +--------------------------------------------------------- + +instance (Floating a, Container Vector a, Floating (Vector a), Fractional (Matrix a)) => Floating (Matrix a) where + sin = liftMatrix sin + cos = liftMatrix cos + tan = liftMatrix tan + asin = liftMatrix asin + acos = liftMatrix acos + atan = liftMatrix atan + sinh = liftMatrix sinh + cosh = liftMatrix cosh + tanh = liftMatrix tanh + asinh = liftMatrix asinh + acosh = liftMatrix acosh + atanh = liftMatrix atanh + exp = liftMatrix exp + log = liftMatrix log + (**) = liftMatrix2Auto (**) + sqrt = liftMatrix sqrt + pi = (1><1) [pi] + +-------------------------------------------------------------------------------- + +isScalar m = rows m == 1 && cols m == 1 + +adaptScalarM f1 f2 f3 x y + | isScalar x = f1 (x @@>(0,0) ) y + | isScalar y = f3 x (y @@>(0,0) ) + | otherwise = f2 x y + +instance (Container Vector t, Eq t, Num (Vector t), Product t) => M.Monoid (Matrix t) + where + mempty = 1 + mappend = adaptScalarM scale mXm (flip scale) + + mconcat xs = work (partition isScalar xs) + where + work (ss,[]) = product ss + work (ss,ms) = scale' (product ss) (optimiseMult ms) + scale' x m + | isScalar x && x00 == 1 = m + | otherwise = scale x00 m + where + x00 = x @@> (0,0) + diff --git a/packages/hmatrix/src/Numeric/Vector.hs b/packages/hmatrix/src/Numeric/Vector.hs new file mode 100644 index 0000000..3f480a0 --- /dev/null +++ b/packages/hmatrix/src/Numeric/Vector.hs @@ -0,0 +1,158 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Vector +-- Copyright : (c) Alberto Ruiz 2011 +-- License : GPL-style +-- +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- Portability : portable +-- +-- Provides instances of standard classes 'Show', 'Read', 'Eq', +-- 'Num', 'Fractional', and 'Floating' for 'Vector'. +-- +----------------------------------------------------------------------------- + +module Numeric.Vector () where + +import Numeric.GSL.Vector +import Numeric.Container + +------------------------------------------------------------------- + +adaptScalar f1 f2 f3 x y + | dim x == 1 = f1 (x@>0) y + | dim y == 1 = f3 x (y@>0) + | otherwise = f2 x y + +------------------------------------------------------------------ + +instance Num (Vector Float) where + (+) = adaptScalar addConstant add (flip addConstant) + negate = scale (-1) + (*) = adaptScalar scale mul (flip scale) + signum = vectorMapF Sign + abs = vectorMapF Abs + fromInteger = fromList . return . fromInteger + +instance Num (Vector Double) where + (+) = adaptScalar addConstant add (flip addConstant) + negate = scale (-1) + (*) = adaptScalar scale mul (flip scale) + signum = vectorMapR Sign + abs = vectorMapR Abs + fromInteger = fromList . return . fromInteger + +instance Num (Vector (Complex Double)) where + (+) = adaptScalar addConstant add (flip addConstant) + negate = scale (-1) + (*) = adaptScalar scale mul (flip scale) + signum = vectorMapC Sign + abs = vectorMapC Abs + fromInteger = fromList . return . fromInteger + +instance Num (Vector (Complex Float)) where + (+) = adaptScalar addConstant add (flip addConstant) + negate = scale (-1) + (*) = adaptScalar scale mul (flip scale) + signum = vectorMapQ Sign + abs = vectorMapQ Abs + fromInteger = fromList . return . fromInteger + +--------------------------------------------------- + +instance (Container Vector a, Num (Vector a)) => Fractional (Vector a) where + fromRational n = fromList [fromRational n] + (/) = adaptScalar f divide g where + r `f` v = scaleRecip r v + v `g` r = scale (recip r) v + +------------------------------------------------------- + +instance Floating (Vector Float) where + sin = vectorMapF Sin + cos = vectorMapF Cos + tan = vectorMapF Tan + asin = vectorMapF ASin + acos = vectorMapF ACos + atan = vectorMapF ATan + sinh = vectorMapF Sinh + cosh = vectorMapF Cosh + tanh = vectorMapF Tanh + asinh = vectorMapF ASinh + acosh = vectorMapF ACosh + atanh = vectorMapF ATanh + exp = vectorMapF Exp + log = vectorMapF Log + sqrt = vectorMapF Sqrt + (**) = adaptScalar (vectorMapValF PowSV) (vectorZipF Pow) (flip (vectorMapValF PowVS)) + pi = fromList [pi] + +------------------------------------------------------------- + +instance Floating (Vector Double) where + sin = vectorMapR Sin + cos = vectorMapR Cos + tan = vectorMapR Tan + asin = vectorMapR ASin + acos = vectorMapR ACos + atan = vectorMapR ATan + sinh = vectorMapR Sinh + cosh = vectorMapR Cosh + tanh = vectorMapR Tanh + asinh = vectorMapR ASinh + acosh = vectorMapR ACosh + atanh = vectorMapR ATanh + exp = vectorMapR Exp + log = vectorMapR Log + sqrt = vectorMapR Sqrt + (**) = adaptScalar (vectorMapValR PowSV) (vectorZipR Pow) (flip (vectorMapValR PowVS)) + pi = fromList [pi] + +------------------------------------------------------------- + +instance Floating (Vector (Complex Double)) where + sin = vectorMapC Sin + cos = vectorMapC Cos + tan = vectorMapC Tan + asin = vectorMapC ASin + acos = vectorMapC ACos + atan = vectorMapC ATan + sinh = vectorMapC Sinh + cosh = vectorMapC Cosh + tanh = vectorMapC Tanh + asinh = vectorMapC ASinh + acosh = vectorMapC ACosh + atanh = vectorMapC ATanh + exp = vectorMapC Exp + log = vectorMapC Log + sqrt = vectorMapC Sqrt + (**) = adaptScalar (vectorMapValC PowSV) (vectorZipC Pow) (flip (vectorMapValC PowVS)) + pi = fromList [pi] + +----------------------------------------------------------- + +instance Floating (Vector (Complex Float)) where + sin = vectorMapQ Sin + cos = vectorMapQ Cos + tan = vectorMapQ Tan + asin = vectorMapQ ASin + acos = vectorMapQ ACos + atan = vectorMapQ ATan + sinh = vectorMapQ Sinh + cosh = vectorMapQ Cosh + tanh = vectorMapQ Tanh + asinh = vectorMapQ ASinh + acosh = vectorMapQ ACosh + atanh = vectorMapQ ATanh + exp = vectorMapQ Exp + log = vectorMapQ Log + sqrt = vectorMapQ Sqrt + (**) = adaptScalar (vectorMapValQ PowSV) (vectorZipQ Pow) (flip (vectorMapValQ PowVS)) + pi = fromList [pi] + -- cgit v1.2.3