Fourier, minimization, polySolve

7 files changed, 377 insertions, 6 deletions

diff --git a/HSSL.cabal b/HSSL.cabal
index 21d98fe..2bd30bc 100644
--- a/HSSL.cabal
+++ b/HSSL.cabal
@@ -24,9 +24,11 @@ Exposed-modules:    Data.Packed.Internal, Data.Packed.Internal.Common,
                     Data.Packed.Internal.Tensor,
                     LAPACK,
                     GSL.Vector,
-                    GSL.Differentiation,
-                    GSL.Integration,
-                    GSL.Special
+                    GSL.Differentiation, GSL.Integration,
+                    GSL.Special,
+                    GSL.Fourier,
+                    GSL.Polynomials,
+                    GSL.Minimization
 Other-modules:
 C-sources:          lib/Data/Packed/Internal/aux.c,
                     lib/LAPACK/lapack-aux.c,
diff --git a/examples/tests.hs b/examples/tests.hs
index 3b1d878..22c4674 100644
--- a/examples/tests.hs
+++ b/examples/tests.hs
@@ -12,11 +12,14 @@ import GSL.Vector
 import GSL.Integration
 import GSL.Differentiation
 import GSL.Special
+import GSL.Fourier
+import GSL.Polynomials
 import LAPACK
 import Test.QuickCheck
 import Test.HUnit
 import Complex
 
+
 {-
 -- Bravo por quickCheck!    
 
@@ -244,6 +247,29 @@ pinvSVDR m = linearSolveSVDR Nothing m (ident (rows m))
 
 pinvSVDC m = linearSolveSVDC Nothing m (ident (rows m))
 
+--------------------------------------------------------------------
+
+polyEval cs x = foldr (\c ac->ac*x+c) 0 cs
+
+polySolveTest' p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p))
+    where l1 |~~| l2 = eps > aproxL magnitude l1 l2
+
+polySolveTest = assertBool "polySolve" (polySolveTest' [1,2,3,4])
+
+---------------------------------------------------------------------
+
+quad f a b = fst $ integrateQAGS 1E-9 100 f a b  
+
+-- A multiple integral can be easily defined using partial application
+quad2 f a b g1 g2 = quad h a b
+    where h x = quad (f x) (g1 x) (g2 x)
+
+volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y)) 
+                        0 r (const 0) (\x->sqrt (r*r-x*x))
+
+integrateTest = assertBool "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < eps)
+
+
 ---------------------------------------------------------------------
 
 arit1 u = vectorMapValR PowVS 2 (vectorMapR Sin u)
@@ -271,6 +297,8 @@ exponentialTest = do
 tests = TestList
     [ TestCase $ besselTest
     , TestCase $ exponentialTest
+    , TestCase $ polySolveTest
+    , TestCase $ integrateTest
     ]
 
 ----------------------------------------------------------------------
@@ -315,9 +343,12 @@ main = do
     quickCheck arit2
     putStrLn "--------- GSL ------"
     runTestTT tests
+    quickCheck $ \v -> ifft (fft v) ~~ v
 
 kk = (2><2)
  [  1.0, 0.0
  , -1.5, 1.0 ::Double]
 
-v = 11 # [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0::Double]
\ No newline at end of file
+v = 11 # [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0::Double]
+
+pol =  [14.125,-7.666666666666667,-14.3,-13.0,-7.0,-9.6,4.666666666666666,13.0,0.5]
\ No newline at end of file
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 2c0acdf..fccf8bb 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -281,3 +281,11 @@ fromColumns m = trans . fromRows $ m
 -- | Creates a list of vectors from the columns of a matrix
 toColumns :: Field t => Matrix t -> [Vector t]
 toColumns m = toRows . trans $ m
+
+
+-- | Reads a matrix position.
+(@@>) :: Field t => Matrix t -> (Int,Int) -> t
+infixl 9 @@>
+m@M {rows = r, cols = c} @@> (i,j)
+    | i<0 || i>=r || j<0 || j>=c = error "matrix indexing out of range"
+    | otherwise   = cdat m `at` (i*c+j)
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index 6ea76c0..ec5744d 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -14,7 +14,7 @@
 
 module Data.Packed.Matrix (
     Matrix(rows,cols), Field,
-    toLists, (><), (>|<),
+    toLists, (><), (>|<), (@@>),
     trans,
     reshape,
     fromRows, toRows, fromColumns, toColumns,
@@ -22,7 +22,9 @@ module Data.Packed.Matrix (
     flipud, fliprl,
     liftMatrix, liftMatrix2,
     multiply,
-    subMatrix, diag, takeDiag, diagRect, ident
+    subMatrix,
+    takeRows, dropRows, takeColumns, dropColumns,
+    diag, takeDiag, diagRect, ident
 ) where
 
 import Data.Packed.Internal
@@ -69,3 +71,20 @@ r >|< c = f where
         | otherwise    = error $ "inconsistent list size = "
                                  ++show (dim v) ++"in ("++show r++"><"++show c++")"
         where v = fromList l
+
+----------------------------------------------------------------
+
+-- | Creates a matrix with the first n rows of another matrix
+takeRows :: Field t => Int -> Matrix t -> Matrix t
+takeRows n mat = subMatrix (0,0) (n, cols mat) mat
+-- | Creates a copy of a matrix without the first n rows
+dropRows :: Field t => Int -> Matrix t -> Matrix t
+dropRows n mat = subMatrix (n,0) (rows mat - n, cols mat) mat
+-- |Creates a matrix with the first n columns of another matrix
+takeColumns :: Field t => Int -> Matrix t -> Matrix t
+takeColumns n mat = subMatrix (0,0) (rows mat, n) mat
+-- | Creates a copy of a matrix without the first n columns
+dropColumns :: Field t => Int -> Matrix t -> Matrix t
+dropColumns n mat = subMatrix (0,n) (rows mat, cols mat - n) mat
+
+----------------------------------------------------------------
diff --git a/lib/GSL/Fourier.hs b/lib/GSL/Fourier.hs
new file mode 100644
index 0000000..c8c79f0
--- /dev/null
+++ b/lib/GSL/Fourier.hs
@@ -0,0 +1,46 @@
+{-# OPTIONS_GHC -fglasgow-exts #-}
+-----------------------------------------------------------------------------
+{- |
+Module      : 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.
+
+<http://www.gnu.org/software/gsl/manual/html_node/Fast-Fourier-Transforms.html#Fast-Fourier-Transforms>
+
+-}
+-----------------------------------------------------------------------------
+module GSL.Fourier (
+    fft,
+    ifft
+) where
+
+import Data.Packed.Internal
+import Complex
+import Foreign
+
+genfft code v = unsafePerformIO $ do
+    r <- createVector (dim v)
+    c_fft code // vec v // vec r // check "fft" [v]
+    return r
+
+foreign import ccall "gsl-aux.h fft" c_fft ::  Int -> Complex Double :> Complex Double :>  IO Int
+
+
+{- | Fast 1D Fourier transform of a 'Vector' @(@'Complex' 'Double'@)@ using /gsl_fft_complex_forward/. It uses the same scaling conventions as GNU Octave.
+
+@> fft ('GSL.Matrix.fromList' [1,2,3,4])
+vector (4) [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/lib/GSL/Minimization.hs b/lib/GSL/Minimization.hs
new file mode 100644
index 0000000..bc228a7
--- /dev/null
+++ b/lib/GSL/Minimization.hs
@@ -0,0 +1,211 @@
+{-# OPTIONS_GHC -fglasgow-exts #-}
+-----------------------------------------------------------------------------
+{- |
+Module      :  GSL.Minimization
+Copyright   :  (c) Alberto Ruiz 2006
+License     :  GPL-style
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  uses ffi
+
+Minimization of a multidimensional function Minimization of a multidimensional function using some of the algorithms described in:
+
+<http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Minimization.html>
+
+-}
+-----------------------------------------------------------------------------
+module GSL.Minimization (
+    minimizeConjugateGradient,
+    minimizeNMSimplex
+) where
+
+
+import Data.Packed.Internal
+import Data.Packed.Matrix
+import Foreign
+import Complex
+
+
+-------------------------------------------------------------------------
+
+{- | The method of Nelder and Mead, implemented by /gsl_multimin_fminimizer_nmsimplex/. The gradient of the function is not required. This is the example in the GSL manual:
+
+@minimize f xi = minimizeNMSimplex f xi (replicate (length xi) 1) 1e-2 100
+\ 
+f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30
+\ 
+main = do
+    let (s,p) = minimize f [5,7]
+    print s
+    print p
+\ 
+\> main
+[0.9920430849306285,1.9969168063253164]
+0. 512.500    1.082 6.500    5.
+ 1. 290.625    1.372 5.250    4.
+ 2. 290.625    1.372 5.250    4.
+ 3. 252.500    1.372 5.500    1.
+ 4. 101.406    1.823 2.625 3.500
+ 5. 101.406    1.823 2.625 3.500
+ 6.     60.    1.823    0.    3.
+ 7.  42.275    1.303 2.094 1.875
+ 8.  42.275    1.303 2.094 1.875
+ 9.  35.684    1.026 0.258 1.906
+10.  35.664    0.804 0.588 2.445
+11.  30.680    0.467 1.258 2.025
+12.  30.680    0.356 1.258 2.025
+13.  30.539    0.285 1.093 1.849
+14.  30.137    0.168 0.883 2.004
+15.  30.137    0.123 0.883 2.004
+16.  30.090    0.100 0.958 2.060
+17.  30.005 6.051e-2 1.022 2.004
+18.  30.005 4.249e-2 1.022 2.004
+19.  30.005 4.249e-2 1.022 2.004
+20.  30.005 2.742e-2 1.022 2.004
+21.  30.005 2.119e-2 1.022 2.004
+22.  30.001 1.530e-2 0.992 1.997
+23.  30.001 1.259e-2 0.992 1.997
+24.  30.001 7.663e-3 0.992 1.997@
+
+The path to the solution can be graphically shown by means of:
+
+@'GSL.Plot.mplot' $ drop 3 ('toColumns' p)@
+
+-}
+minimizeNMSimplex :: ([Double] -> Double) -- ^ function to minimize
+          -> [Double]            -- ^ starting point
+          -> [Double]            -- ^ sizes of the initial search box
+          -> Double              -- ^ desired precision of the solution
+          -> Int                 -- ^ maximum number of iterations allowed
+          -> ([Double], Matrix Double)
+          -- ^ solution vector, and the optimization trajectory followed by the algorithm
+minimizeNMSimplex f xi sz tol maxit = unsafePerformIO $ do
+    let xiv = fromList xi
+        szv = fromList sz
+        n   = dim xiv
+    fp <- mkVecfun (iv (f.toList))
+    rawpath <- createMIO maxit (n+3)
+                         (c_minimizeNMSimplex fp tol maxit // vec xiv // vec szv) 
+                         "minimizeNMSimplex" [xiv,szv]
+    let it = round (rawpath @@> (maxit-1,0))
+        path = takeRows it rawpath
+        [sol] = toLists $ dropRows (it-1) path
+    freeHaskellFunPtr fp
+    return (drop 3 sol, path)
+
+
+foreign import ccall "gsl-aux.h minimize"
+    c_minimizeNMSimplex:: FunPtr (Int -> Ptr Double -> Double) -> Double -> Int
+                       -> Double :> Double :> Double ::> IO Int
+
+----------------------------------------------------------------------------------
+
+{- | The Fletcher-Reeves conjugate gradient algorithm /gsl_multimin_fminimizer_conjugate_fr/. This is the example in the GSL manual:
+
+@minimize = minimizeConjugateGradient 1E-2 1E-4 1E-3 30
+f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30
+\ 
+df [x,y] = [20*(x-1), 40*(y-2)]
+\  
+main = do
+    let (s,p) = minimize f df [5,7]
+    print s
+    print p
+\ 
+\> main
+[1.0,2.0]
+ 0. 687.848 4.996 6.991
+ 1. 683.555 4.989 6.972
+ 2. 675.013 4.974 6.935
+ 3. 658.108 4.944 6.861
+ 4. 625.013 4.885 6.712
+ 5. 561.684 4.766 6.415
+ 6. 446.467 4.528 5.821
+ 7. 261.794 4.053 4.632
+ 8.  75.498 3.102 2.255
+ 9.  67.037 2.852 1.630
+10.  45.316 2.191 1.762
+11.  30.186 0.869 2.026
+12.     30.    1.    2.@
+
+The path to the solution can be graphically shown by means of:
+
+@'GSL.Plot.mplot' $ drop 2 ('toColumns' p)@
+
+-}     
+minimizeConjugateGradient ::
+       Double        -- ^ initial step size
+    -> Double        -- ^ minimization parameter   
+    -> Double        -- ^ desired precision of the solution (gradient test)
+    -> Int           -- ^ maximum number of iterations allowed
+    -> ([Double] -> Double) -- ^ function to minimize
+    -> ([Double] -> [Double])      -- ^ gradient
+    -> [Double]             -- ^ starting point
+    -> ([Double], Matrix Double)        -- ^ solution vector, and the optimization trajectory followed by the algorithm
+minimizeConjugateGradient istep minimpar tol maxit f df xi = unsafePerformIO $ do
+    let xiv = fromList xi
+        n = dim xiv
+        f' = f . toList
+        df' = (fromList . df . toList)
+    fp <- mkVecfun (iv f')
+    dfp <- mkVecVecfun (aux_vTov df')
+    rawpath <- createMIO maxit (n+2)
+                         (c_minimizeConjugateGradient fp dfp istep minimpar tol maxit // vec xiv)
+                         "minimizeDerivV" [xiv]
+    let it = round (rawpath @@> (maxit-1,0))
+        path = takeRows it rawpath
+        sol = toList $ cdat $ dropColumns 2 $ dropRows (it-1) path
+    freeHaskellFunPtr fp
+    freeHaskellFunPtr dfp
+    return (sol,path)
+
+
+foreign import ccall "gsl-aux.h minimizeWithDeriv"
+    c_minimizeConjugateGradient :: FunPtr (Int -> Ptr Double -> Double)
+                                -> FunPtr (Int -> Ptr Double -> Ptr Double -> IO ())
+                                -> Double -> Double -> Double -> Int 
+                                -> Double :> Double ::> IO Int
+
+---------------------------------------------------------------------
+iv :: (Vector Double -> Double) -> (Int -> Ptr Double -> Double)
+iv f n p = f (createV n copy "iv" []) where
+    copy n q = do 
+        copyArray q p n
+        return 0
+
+-- | conversion of Haskell functions into function pointers that can be used in the C side
+foreign import ccall "wrapper"
+    mkVecfun :: (Int -> Ptr Double -> Double)
+             -> IO( FunPtr (Int -> Ptr Double -> Double))
+
+-- | another required conversion
+foreign import ccall "wrapper"
+    mkVecVecfun :: (Int -> Ptr Double -> Ptr Double -> IO ())
+                -> IO (FunPtr (Int -> Ptr Double -> Ptr Double->IO()))
+
+aux_vTov :: (Vector Double -> Vector Double) -> (Int -> Ptr Double -> Ptr Double -> IO())
+aux_vTov f n p r = g where
+    V {fptr = pr, ptr = p} = f x
+    x = createV n copy "aux_vTov" []
+    copy n q = do 
+        copyArray q p n
+        return 0
+    g = withForeignPtr pr $ \_ -> copyArray r p n
+
+--------------------------------------------------------------------
+
+createV n fun msg ptrs = unsafePerformIO $ do
+    r <- createVector n
+    fun // vec r // check msg ptrs
+    return r
+
+createM r c fun msg ptrs = unsafePerformIO $ do
+    r <- createMatrix RowMajor r c
+    fun // mat cdat r // check msg ptrs
+    return r
+
+createMIO r c fun msg ptrs = do
+    r <- createMatrix RowMajor r c
+    fun // mat cdat r // check msg ptrs
+    return r
diff --git a/lib/GSL/Polynomials.hs b/lib/GSL/Polynomials.hs
new file mode 100644
index 0000000..d7db9a3
--- /dev/null
+++ b/lib/GSL/Polynomials.hs
@@ -0,0 +1,54 @@
+{-# OPTIONS_GHC -fglasgow-exts #-}
+-----------------------------------------------------------------------------
+{- |
+Module      :  GSL.Polynomials
+Copyright   :  (c) Alberto Ruiz 2006
+License     :  GPL-style
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  uses ffi
+
+Polynomials.
+
+<http://www.gnu.org/software/gsl/manual/html_node/General-Polynomial-Equations.html#General-Polynomial-Equations>
+
+-}
+-----------------------------------------------------------------------------
+module GSL.Polynomials (
+    polySolve
+) where
+
+import Data.Packed.Internal
+import Complex
+import Foreign
+
+{- | 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]
+[(-1.9999999999999998) :+ 0.0,
+ 1.0 :+ 1.732050807568877,
+ 1.0 :+ (-1.732050807568877)]@
+
+The example in the GSL manual: To find the roots of x^5 -1 = 0:
+
+@\> polySolve [-1, 0, 0, 0, 0, 1]
+[(-0.8090169943749475) :+ 0.5877852522924731,
+(-0.8090169943749475) :+ (-0.5877852522924731),
+0.30901699437494734 :+ 0.9510565162951536,
+0.30901699437494734 :+ (-0.9510565162951536),
+1.0 :+ 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)
+    c_polySolve // vec v // vec r // check "polySolve" [v]
+    return r
+             | otherwise = error "polySolve on a polynomial of degree zero"
+
+foreign import ccall "gsl-aux.h polySolve" c_polySolve:: Double :> Complex Double :> IO Int