diff --git a/configure b/configure deleted file mode 100644 --- a/configure +++ /dev/null @@ -1,3 +0,0 @@ -#! /bin/sh - -runhaskell configure.hs $* diff --git a/configure.hs b/configure.hs deleted file mode 100644 --- a/configure.hs +++ /dev/null @@ -1,145 +0,0 @@ -#! /usr/bin/env runhaskell -{- configure.hs for hstatistics, copied from hmatrix - ------------------------ - - GSL and LAPACK may require auxiliary libraries which depend on OS, - distribution, and implementation. This script tries to to find out - the correct link command for your system. - Suggestions and contributions are welcome. - - By default we try to link -lgsl -llapack. This works in ubuntu/debian, - both with and without ATLAS. - If this fails we try different sets of additional libraries which are - known to work in some systems. - - The desired libraries can also be explicitly given by the user using cabal - flags (e.g., -fmkl, -faccelerate) or --configure-option=link:lib1,lib2,lib3,... - --} - -import System -import System.Directory(createDirectoryIfMissing,getCurrentDirectory,getDirectoryContents) -import Data.List(isPrefixOf, intercalate) -import Distribution.Simple.LocalBuildInfo -import Distribution.Simple.Configure -import Distribution.PackageDescription - --- possible additional dependencies for the desired libs (by default gsl lapack) - -opts = [ "" -- Ubuntu/Debian - , "blas" - , "blas cblas" - , "cblas" - , "gslcblas" - , "blas gslcblas" - , "f77blas" - , "f77blas cblas atlas gcc_s" -- Arch Linux (older version of atlas-lapack) - , "blas gslcblas gfortran" -- Arch Linux with normal blas and lapack - ] - --- compile a simple program with symbols from GSL and LAPACK with the given libs -testprog bInfo buildInfo libs fmks = - "echo \"#include \nint main(){zgesvd_(); gsl_sf_gamma(5);}\"" - ++" > " ++ (buildDir bInfo) ++ "/dummy.c; gcc " - ++ (join $ ccOptions buildInfo) ++ " " - ++ (join $ cppOptions buildInfo) ++ " " - ++ (join $ map ("-I"++) $ includeDirs buildInfo) ++ " " - ++ (buildDir bInfo) ++ "/dummy.c -o " - ++ (buildDir bInfo) ++ "/dummy " - ++ (join $ map ("-L"++) $ extraLibDirs buildInfo) ++ " " - ++ (prepend "-l" $ libs) ++ " " - ++ (prepend "-framework " fmks) ++ " > /dev/null 2> /dev/null" - -join = intercalate " " -prepend x = unwords . map (x++) . words - -check bInfo buildInfo libs fmks = (ExitSuccess ==) `fmap` system (testprog bInfo buildInfo libs fmks) - --- simple test for GSL -gsl bInfo buildInfo = "echo \"#include \nint main(){gsl_sf_gamma(5);}\"" - ++" > " ++ (buildDir bInfo) ++ "/dummy.c; gcc " - ++ (join $ ccOptions buildInfo) ++ " " - ++ (join $ cppOptions buildInfo) ++ " " - ++ (join $ map ("-I"++) $ includeDirs buildInfo) ++ " " - ++ (buildDir bInfo) ++ "/dummy.c -o " - ++ (buildDir bInfo) ++ "/dummy " - ++ (join $ map ("-L"++) $ extraLibDirs buildInfo) ++ " -lgsl -lgslcblas" - ++ " > /dev/null 2> /dev/null" - --- test for gsl >= 1.12 -gsl112 bInfo buildInfo = - "echo \"#include \nint main(){gsl_sf_exprel_n_CF_e(1,1,0);}\"" - ++" > " ++ (buildDir bInfo) ++ "/dummy.c; gcc " - ++ (buildDir bInfo) ++ "/dummy.c " - ++ (join $ ccOptions buildInfo) ++ " " - ++ (join $ cppOptions buildInfo) ++ " " - ++ (join $ map ("-I"++) $ includeDirs buildInfo) - ++" -o " ++ (buildDir bInfo) ++ "/dummy " - ++ (join $ map ("-L"++) $ extraLibDirs buildInfo) ++ " -lgsl -lgslcblas" - ++ " > /dev/null 2> /dev/null" - - -checkCommand c = (ExitSuccess ==) `fmap` system c - --- test different configurations until the first one works -try _ _ _ _ [] = return Nothing -try l i b f (opt:rest) = do - ok <- check l i (b ++ " " ++ opt) f - if ok then return (Just opt) - else try l i b f rest - --- read --configure-option=link:lib1,lib2,lib3,etc -linkop = "link:" -getUserLink = concatMap (g . drop (length linkop)) . filter (isPrefixOf linkop) - where g = map cs - cs ',' = ' ' - cs x = x - -main = do - dir <- getCurrentDirectory - putStrLn $ "Current directory: " ++ dir - --- files <- getDirectoryContents dir --- mapM_ putStrLn files - - putStr "Checking foreign libraries..." - - args <- getArgs - Just bInfo <- maybeGetPersistBuildConfig "dist" - - let Just lib = library . localPkgDescr $ bInfo - buildInfo = libBuildInfo lib - base = unwords . extraLibs $ buildInfo - fwks = unwords . frameworks $ buildInfo - auxpref = getUserLink args - - -- We extract the desired libs from hstatistics.cabal (using a cabal flags) - -- and from a posible --configure-option=link:lib1,lib2,lib3 - -- by default the desired libs are gsl lapack. - - let pref = if null (words (base ++ " " ++ auxpref)) then "gsl lapack" else auxpref - fullOpts = map ((pref++" ")++) opts - - -- create the build directory (used for tmp files) if necessary - createDirectoryIfMissing True $ buildDir bInfo - - r <- try bInfo buildInfo base fwks fullOpts - case r of - Nothing -> do - putStrLn " FAIL" - g <- checkCommand $ gsl bInfo buildInfo - if g - then putStrLn " *** Sorry, I can't link LAPACK." - else putStrLn " *** Sorry, I can't link GSL." - putStrLn " *** Please make sure that the appropriate -dev packages are installed." - putStrLn " *** You can also specify the required libraries using" - putStrLn " *** cabal install hstatistics --configure-option=link:lib1,lib2,lib3,etc." - writeFile "hstatistics.buildinfo" ("buildable: False\n") - Just ops -> do - putStrLn " OK" - g <- checkCommand $ gsl112 bInfo buildInfo - writeFile "hstatistics.buildinfo" $ "extra-libraries: " ++ - ops ++ "\n" ++ - if g - then "" - else "cc-options: -DGSL110\n" diff --git a/hstatistics.cabal b/hstatistics.cabal --- a/hstatistics.cabal +++ b/hstatistics.cabal @@ -1,5 +1,5 @@ Name: hstatistics -Version: 0.1.0.5 +Version: 0.2.0.1 License: GPL License-file: LICENSE Copyright: (c) A.V.H. McPhail 2010 @@ -8,59 +8,34 @@ Stability: provisional Homepage: http://code.haskell.org/hstatistics Synopsis: Statistics -Description: Purely functional interface for statistics based on hmatrix and GSL +Description: Purely functional interface for statistics based on hmatrix and hmatrix-gsl-stats Category: Math tested-with: GHC ==6.12.1 cabal-version: >=1.2 -build-type: Custom +build-type: Simple -extra-source-files: configure configure.hs README INSTALL CHANGES +extra-source-files: README INSTALL CHANGES extra-tmp-files: hstatistics.buildinfo library Build-Depends: base >= 3 && < 5, - storable-complex, hmatrix >= 0.9.3 + hmatrix >= 0.9.3, hmatrix-gsl-stats >= 0.1.0.1 - Extensions: ForeignFunctionInterface + Extensions: hs-source-dirs: lib - Exposed-modules: Numeric.GSL.Sort - Numeric.GSL.Statistics - Numeric.GSL.Histogram - Numeric.GSL.Histogram2D - Numeric.GSL.Permutation - Numeric.GSL.Distribution.Continuous - Numeric.GSL.Distribution.Discrete - Numeric.GSL.Distribution.Common - Numeric.GSL.Fitting.Linear + Exposed-modules: Numeric.Statistics.Shannon other-modules: - C-sources: lib/Numeric/GSL/statistics-aux.c - lib/Numeric/GSL/sort-aux.c - lib/Numeric/GSL/histogram-aux.c - lib/Numeric/GSL/permutation-aux.c - lib/Numeric/GSL/distribution-aux.c - lib/Numeric/GSL/fitting-aux.c + C-sources: ghc-prof-options: -auto ghc-options: -Wall -fno-warn-missing-signatures -fno-warn-orphans -fno-warn-unused-binds - - if os(OSX) - extra-lib-dirs: /opt/local/lib/ - include-dirs: /opt/local/include/ - extra-libraries: gsl - frameworks: Accelerate - --- The extra-libraries required for GSL --- should now be automatically detected by configure(.hs) - - extra-libraries: - extra-lib-dirs: source-repository head type: darcs diff --git a/lib/Numeric/GSL/Distribution/Common.hs b/lib/Numeric/GSL/Distribution/Common.hs deleted file mode 100644 --- a/lib/Numeric/GSL/Distribution/Common.hs +++ /dev/null @@ -1,29 +0,0 @@ -{-# OPTIONS_GHC -fglasgow-exts #-} ------------------------------------------------------------------------------ --- | --- Module : Numeric.GSL.Distribution.Common --- Copyright : (c) Alexander Vivian Hugh McPhail 2010 --- License : GPL-style --- --- Maintainer : haskell.vivian.mcphail gmail com --- Stability : provisional --- Portability : uses ffi --- --- GSL common data types for distributions --- ------------------------------------------------------------------------------ - -module Numeric.GSL.Distribution.Common ( - DistFunc(..) - ) where - ------------------------------------------------------------------------------ - -data DistFunc = Density -- ^ pdf - | Lower -- ^ lower cdf - | Upper -- ^ upper cdf - | LowInv -- ^ lower inverse cdf - | UppInv -- ^ upper inverse cdf - deriving(Enum,Eq) - ------------------------------------------------------------------------------ diff --git a/lib/Numeric/GSL/Distribution/Continuous.hs b/lib/Numeric/GSL/Distribution/Continuous.hs deleted file mode 100644 --- a/lib/Numeric/GSL/Distribution/Continuous.hs +++ /dev/null @@ -1,299 +0,0 @@ -{-# OPTIONS_GHC -fglasgow-exts #-} ------------------------------------------------------------------------------ --- | --- Module : Numeric.GSL.Distribution.Continuous --- Copyright : (c) Alexander Vivian Hugh McPhail 2010 --- License : GPL-style --- --- Maintainer : haskell.vivian.mcphail gmail com --- Stability : provisional --- Portability : uses ffi --- --- GSL continuous random distribution functions --- ------------------------------------------------------------------------------ - -module Numeric.GSL.Distribution.Continuous ( - ZeroParamDist(..), OneParamDist(..) - , TwoParamDist(..), ThreeParamDist(..) - , MultiParamDist(..) - , BivariateDist(..) - , DistFunc(..) - , random_0p, density_0p - , random_1p, density_1p - , random_2p, density_2p - , random_3p, density_3p - , random_mp, density_mp - , random_biv, density_biv - , spherical_vector - ) where - ------------------------------------------------------------------------------ - -import Data.Packed.Vector ---import Data.Packed.Matrix hiding(toLists) -import Data.Packed.Development - ---import Numeric.LinearAlgebra.Linear - ---import Control.Monad(when) - -import Foreign hiding(shift) ---import Foreign.ForeignPtr ---import Foreign.Marshal.Alloc(alloca) ---import Foreign.C.Types(CInt,CChar) -import Foreign.C.Types(CInt) ---import Foreign.C.String(newCString,peekCString) - ---import GHC.ForeignPtr (mallocPlainForeignPtrBytes) - ---import GHC.Base ---import GHC.IOBase - ---import Prelude hiding(reverse) - -import Numeric.GSL.Distribution.Common - ------------------------------------------------------------------------------ - -data ZeroParamDist = Landau - deriving Enum - -data OneParamDist = Gaussian -- ^ standard deviation - | Exponential -- ^ mean - | Laplace -- ^ width - | Cauchy -- ^ scale - | Rayleigh -- ^ standard deviation - | ChiSq -- ^ degrees of freedom - | TDist -- ^ degrees of freedom - | Logistic -- ^ scale - deriving Enum - -data TwoParamDist = GaussianTail -- ^ limit, standard deviation - | ExpPower -- ^ scale, exponent - | RayleighTail -- ^ lower limit, standard deviation - | Levy -- ^ scale, exponent - | Gamma -- ^ par1, par2 - | Uniform -- ^ lower, upper - | Lognormal -- ^ offset, standard deviation - | FDist -- ^ degrees of freedom, degrees of freedom - | Beta -- ^ parameter a, parameter b - | Pareto -- ^ exponent, scale - | Weibull -- ^ scale, exponent - | GumbellI -- ^ A, B - | GumbellII -- ^ A, B - deriving Enum - -data ThreeParamDist = LevySkew -- ^ scale, exponent, skewness - deriving Enum - -data MultiParamDist = Dirichlet -- ^ size, alpha - deriving Enum - -data BivariateDist = BiGaussian -- ^ standard deviation, standard deviation, correlation coefficient - deriving Enum - ------------------------------------------------------------------------------ - -fromei x = fromIntegral (fromEnum x) :: CInt - ------------------------------------------------------------------------------ - --- | draw a sample from a zero parameter distribution -random_0p :: ZeroParamDist -- ^ distribution type - -> Int -- ^ random seed - -> Double -- ^ result -random_0p d s = unsafePerformIO $ distribution_random_zero_param (fromIntegral s) (fromei d) - --- | probability of a variate take a value outside the argument -density_0p :: ZeroParamDist -- ^ density type - -> DistFunc -- ^ distribution function type - -> Double -- ^ value - -> Double -- ^ result -density_0p d f x = unsafePerformIO $ do - case d of - Landau -> density_only f d x - where density_only f' d' x' = if f' /= Density - then error "distribution has no CDF" - else distribution_dist_zero_param (fromei f') (fromei d') x' - -foreign import ccall "distribution-aux.h random0" distribution_random_zero_param :: CInt -> CInt -> IO Double -foreign import ccall "distribution-aux.h random0_dist" distribution_dist_zero_param :: CInt -> CInt -> Double -> IO Double - ------------------------------------------------------------------------------ - --- | draw a sample from a one parameter distribution -random_1p :: OneParamDist -- ^ distribution type - -> Int -- ^ random seed - -> Double -- ^ parameter - -> Double -- ^ result -random_1p d s p = unsafePerformIO $ - alloca $ \r -> do - check "random1p" $ distribution_random_one_param (fromIntegral s) (fromei d) p r - r' <- peek r - return r' - --- | probability of a variate take a value outside the argument -density_1p :: OneParamDist -- ^ density type - -> DistFunc -- ^ distribution function type - -> Double -- ^ parameter - -> Double -- ^ value - -> Double -- ^ result -density_1p d f p x = unsafePerformIO $ distribution_dist_one_param (fromei f) (fromei d) x p - -foreign import ccall "distribution-aux.h random1" distribution_random_one_param :: CInt -> CInt -> Double -> Ptr Double -> IO CInt -foreign import ccall "distribution-aux.h random1_dist" distribution_dist_one_param :: CInt -> CInt -> Double -> Double -> IO Double - ------------------------------------------------------------------------------ - --- | draw a sample from a two parameter distribution -random_2p :: TwoParamDist -- ^ distribution type - -> Int -- ^ random seed - -> Double -- ^ parameter 1 - -> Double -- ^ parameter 2 - -> Double -- ^ result -random_2p d s p1 p2 = unsafePerformIO $ - alloca $ \r -> do - check "random2p" $ distribution_random_two_param (fromIntegral s) (fromei d) p1 p2 r - r' <- peek r - return r' - --- | probability of a variate take a value outside the argument -density_2p :: TwoParamDist -- ^ density type - -> DistFunc -- ^ distribution function type - -> Double -- ^ parameter 1 - -> Double -- ^ parameter 2 - -> Double -- ^ value - -> Double -- ^ result -density_2p d f p1 p2 x = unsafePerformIO $ do - case d of - GaussianTail -> density_only f d p1 p2 x - ExpPower -> no_inverse f d p1 p2 x - RayleighTail -> density_only f d p1 p2 x - Levy -> error "no PDF or CDF for Levy" - _ -> distribution_dist_two_param (fromei f) (fromei d) x p1 p2 - where density_only f' d' p1' p2' x' = if f' /= Density - then error "distribution has no CDF" - else distribution_dist_two_param (fromei f') (fromei d') x' p1' p2' - no_inverse f' d' p1' p2' x' = if (f' == LowInv || f' == UppInv) - then error "distribution has no inverse CDF" - else distribution_dist_two_param (fromei f') (fromei d') x' p1' p2' - -foreign import ccall "distribution-aux.h random2" distribution_random_two_param :: CInt -> CInt -> Double -> Double -> Ptr Double -> IO CInt -foreign import ccall "distribution-aux.h random2_dist" distribution_dist_two_param :: CInt -> CInt -> Double -> Double -> Double -> IO Double - ------------------------------------------------------------------------------ - --- | draw a sample from a three parameter distribution -random_3p :: ThreeParamDist -- ^ distribution type - -> Int -- ^ random seed - -> Double -- ^ parameter 1 - -> Double -- ^ parameter 2 - -> Double -- ^ parameter 3 - -> Double -- ^ result -random_3p d s p1 p2 p3 = unsafePerformIO $ - alloca $ \r -> do - check "random_3p" $ distribution_random_three_param (fromIntegral s) (fromei d) p1 p2 p3 r - r' <- peek r - return r' - --- | probability of a variate take a value outside the argument -density_3p :: ThreeParamDist -- ^ density type - -> DistFunc -- ^ distribution function type - -> Double -- ^ parameter 1 - -> Double -- ^ parameter 2 - -> Double -- ^ parameter 3 - -> Double -- ^ value - -> Double -- ^ result -density_3p d f p1 p2 p3 x = unsafePerformIO $ do - case d of - LevySkew -> error "Levy skew has no PDF or CDF" - where density_only f' d' p1' p2' p3' x' = if f' /= Density - then error "distribution has no CDF" - else distribution_dist_three_param (fromei f') (fromei d') x' p1' p2' p3' - -foreign import ccall "distribution-aux.h random3" distribution_random_three_param :: CInt -> CInt -> Double -> Double -> Double -> Ptr Double -> IO CInt -foreign import ccall "distribution-aux.h random3_dist" distribution_dist_three_param :: CInt -> CInt -> Double -> Double -> Double -> Double -> IO Double - ------------------------------------------------------------------------------ - --- | draw a sample from a three parameter distribution -random_mp :: MultiParamDist -- ^ distribution type - -> Int -- ^ random seed - -> Vector Double -- ^ parameters - -> Vector Double -- ^ result -random_mp d s p = unsafePerformIO $ do - r <- createVector $ dim p - app2 (distribution_random_multi_param (fromIntegral s) (fromei d)) vec p vec r "random_mp" - return r - --- | probability of a variate take a value outside the argument -density_mp :: MultiParamDist -- ^ density type - -> DistFunc -- ^ distribution function type - -> Vector Double -- ^ parameters - -> Vector Double -- ^ values - -> Double -- ^ result -density_mp d f p q = unsafePerformIO $ do - case d of - Dirichlet -> density_only f d p q - where density_only f' d' p' q' = if f' /= Density - then error "distribution has no CDF" - else alloca $ \r -> do - app2 (distribution_dist_multi_param (fromei f') (fromei d') r) vec p vec q "density_mp" - r' <- peek r - return r' - -foreign import ccall "distribution-aux.h random_mp" distribution_random_multi_param :: CInt -> CInt -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt -foreign import ccall "distribution-aux.h random_mm_dist" distribution_dist_multi_param :: CInt -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt - ------------------------------------------------------------------------------ - --- | draw a sample from a bivariate distribution -random_biv :: BivariateDist -- ^ distribution type - -> Int -- ^ random seed - -> Double -- ^ parameter 1 - -> Double -- ^ parameter 2 - -> Double -- ^ parameter 3 - -> (Double,Double) -- ^ result -random_biv d s p1 p2 p3 = unsafePerformIO $ - alloca $ \r1 -> - alloca $ \r2 -> do - distribution_random_bivariate (fromIntegral s) (fromei d) p1 p2 p3 r1 r2 - r1' <- peek r1 - r2' <- peek r2 - return (r1',r2') - --- | probability of a variate take a value outside the argument -density_biv :: BivariateDist -- ^ density type - -> DistFunc -- ^ distribution function type - -> Double -- ^ parameter 1 - -> Double -- ^ parameter 2 - -> Double -- ^ parameter 3 - -> (Double,Double) -- ^ value - -> Double -- ^ result -density_biv d f p1 p2 p3 (x,y) = unsafePerformIO $ do - case d of - BiGaussian -> density_only f d p1 p2 p3 x y - where density_only f' d' p1' p2' p3' x' y' = if f' /= Density - then error "distribution has no CDF" - else distribution_dist_bivariate (fromei f') (fromei d') x' y' p1' p2' p3' - -foreign import ccall "distribution-aux.h random_biv" distribution_random_bivariate :: CInt -> CInt -> Double -> Double -> Double -> Ptr Double -> Ptr Double -> IO () -foreign import ccall "distribution-aux.h random_biv_dist" distribution_dist_bivariate :: CInt -> CInt -> Double -> Double -> Double -> Double -> Double -> IO Double - ------------------------------------------------------------------------------ - --- | returns a normalised random direction vector from a multivariate gaussian distribution -spherical_vector :: Int -- ^ seed - -> Int -- ^ vector size - -> Vector Double -- result -spherical_vector s vs = unsafePerformIO $ do - r <- createVector vs - app1 (distribution_spherical_vector (fromIntegral s)) vec r "spherical_vector" - return r - -foreign import ccall "distribution-aux.h random_vector" distribution_spherical_vector :: CInt -> CInt -> Ptr Double -> IO CInt - ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ diff --git a/lib/Numeric/GSL/Distribution/Discrete.hs b/lib/Numeric/GSL/Distribution/Discrete.hs deleted file mode 100644 --- a/lib/Numeric/GSL/Distribution/Discrete.hs +++ /dev/null @@ -1,210 +0,0 @@ -{-# OPTIONS_GHC -fglasgow-exts #-} ------------------------------------------------------------------------------ --- | --- Module : Numeric.GSL.Distribution.Discrete --- Copyright : (c) Alexander Vivian Hugh McPhail 2010 --- License : GPL-style --- --- Maintainer : haskell.vivian.mcphail gmail com --- Stability : provisional --- Portability : uses ffi --- --- GSL discrete random distribution functions --- ------------------------------------------------------------------------------ - -module Numeric.GSL.Distribution.Discrete ( - OneParamDist(..) - , TwoParamDist(..), ThreeParamDist(..) - , MultiParamDist(..) - , DistFunc(..) - , random_1p, density_1p - , random_2p, density_2p - , random_3p, density_3p - , random_mp, density_mp - ) where - ------------------------------------------------------------------------------ - -import Data.Packed.Vector ---import Data.Packed.Matrix hiding(toLists) -import Data.Packed.Development - ---import Numeric.LinearAlgebra.Linear - ---import Control.Monad(when) - -import Foreign hiding(shift) ---import Foreign.ForeignPtr ---import Foreign.Marshal.Alloc(alloca) ---import Foreign.C.Types(CInt,CUInt,CChar) -import Foreign.C.Types(CInt,CUInt) ---import Foreign.C.String(newCString,peekCString) - ---import GHC.ForeignPtr (mallocPlainForeignPtrBytes) - ---import GHC.Base ---import GHC.IOBase - ---import Prelude hiding(reverse) - -import Numeric.GSL.Distribution.Common - ------------------------------------------------------------------------------ - -data OneParamDist = Poisson -- ^ mean - | Bernoulli -- ^ probability - | Geometric -- ^ probability - | Logarithmic -- ^ probability - deriving Enum - -data TwoParamDist = Binomial -- ^ probability, successes - | NegBinomial -- ^ probability, successes - | Pascal -- ^ probability, n - deriving Enum - -data ThreeParamDist = HyperGeometric -- ^ number type 1, number type 2, samples - deriving Enum - -data MultiParamDist = Multinomial -- ^ trials, probabilities - deriving Enum - ------------------------------------------------------------------------------ - -fromei x = fromIntegral (fromEnum x) :: CInt - ------------------------------------------------------------------------------ - --- | draw a sample from a one parameter distribution -random_1p :: OneParamDist -- ^ distribution type - -> Int -- ^ random seed - -> Double -- ^ parameter - -> Int -- ^ result -random_1p d s p = unsafePerformIO $ - alloca $ \r -> do - check "random1p" $ distribution_discrete_one_param (fromIntegral s) (fromei d) p r - r' <- peek r - return $ fromIntegral r' - --- | probability of a variate take a value outside the argument -density_1p :: OneParamDist -- ^ density type - -> DistFunc -- ^ distribution function type - -> Double -- ^ parameter - -> Int -- ^ value - -> Double -- ^ result -density_1p d f p x = unsafePerformIO $ do - case d of - Poisson -> no_inverse f d p x - Geometric -> no_inverse f d p x - Logarithmic -> pdf_only f d p x - Bernoulli -> pdf_only f d p x - _ -> distribution_dist_one_param (fromei f) (fromei d) (fromIntegral x) p - where pdf_only f' d' p' x' = if f' /= Density - then error "no CDF" - else distribution_dist_one_param (fromei f) (fromei d) (fromIntegral x) p - no_inverse f' d' p' x' = if (f' == LowInv || f' == UppInv) - then error "No inverse CDF" - else distribution_dist_one_param (fromei f) (fromei d) (fromIntegral x) p - -foreign import ccall "distribution-aux.h discrete1" distribution_discrete_one_param :: CInt -> CInt -> Double -> Ptr CUInt -> IO CInt -foreign import ccall "distribution-aux.h discrete1_dist" distribution_dist_one_param :: CInt -> CInt -> CUInt -> Double -> IO Double - ------------------------------------------------------------------------------ - --- | draw a sample from a two parameter distribution -random_2p :: TwoParamDist -- ^ distribution type - -> Int -- ^ random seed - -> Double -- ^ parameter 1 - -> Int -- ^ parameter 2 - -> Int -- ^ result -random_2p d s p1 p2 = unsafePerformIO $ - alloca $ \r -> do - check "random2p" $ distribution_discrete_two_param (fromIntegral s) (fromei d) p1 (fromIntegral p2) r - r' <- peek r - return $ fromIntegral r' - --- | probability of a variate take a value outside the argument -density_2p :: TwoParamDist -- ^ density type - -> DistFunc -- ^ distribution function type - -> Double -- ^ parameter 1 - -> Int -- ^ parameter 2 - -> Int -- ^ value - -> Double -- ^ result -density_2p d f p1 p2 x = unsafePerformIO $ do - case d of - _ -> no_inverse f d p1 p2 x - where no_inverse f' d' p1' p2' x' = if (f' == LowInv || f' == UppInv) - then error "distribution has no inverse CDF" - else distribution_dist_two_param (fromei f') (fromei d') (fromIntegral x') p1' (fromIntegral p2') - -foreign import ccall "distribution-aux.h discrete2" distribution_discrete_two_param :: CInt -> CInt -> Double -> CUInt -> Ptr CUInt -> IO CInt -foreign import ccall "distribution-aux.h discrete2_dist" distribution_dist_two_param :: CInt -> CInt -> CUInt -> Double -> CUInt -> IO Double - ------------------------------------------------------------------------------ - --- | draw a sample from a three parameter distribution -random_3p :: ThreeParamDist -- ^ distribution type - -> Int -- ^ random seed - -> Int -- ^ parameter 1 - -> Int -- ^ parameter 2 - -> Int -- ^ parameter 3 - -> Int -- ^ result -random_3p d s p1 p2 p3 = unsafePerformIO $ - alloca $ \r -> do - check "random_3p" $ distribution_discrete_three_param (fromIntegral s) (fromei d) (fromIntegral p1) (fromIntegral p2) (fromIntegral p3) r - r' <- peek r - return $ fromIntegral r' - --- | probability of a variate take a value outside the argument -density_3p :: ThreeParamDist -- ^ density type - -> DistFunc -- ^ distribution function type - -> Int -- ^ parameter 1 - -> Int -- ^ parameter 2 - -> Int -- ^ parameter 3 - -> Int -- ^ value - -> Double -- ^ result -density_3p d f p1 p2 p3 x = unsafePerformIO $ do - case d of - HyperGeometric -> no_inverse f d p1 p2 p3 x - where no_inverse f' d' p1' p2' p3' x' = if (f' == LowInv || f' == UppInv) - then error "No inverse CDF" - else distribution_dist_three_param (fromei f') (fromei d') (fromIntegral x') (fromIntegral p1') (fromIntegral p1') (fromIntegral p1') - -foreign import ccall "distribution-aux.h discrete3" distribution_discrete_three_param :: CInt -> CInt -> CUInt -> CUInt -> CUInt -> Ptr CUInt -> IO CInt -foreign import ccall "distribution-aux.h discrete3_dist" distribution_dist_three_param :: CInt -> CInt -> CUInt -> CUInt -> CUInt -> CUInt -> IO Double - ------------------------------------------------------------------------------ - --- | draw a sample from a three parameter distribution -random_mp :: MultiParamDist -- ^ distribution type - -> Int -- ^ random seed - -> Int -- ^ trials - -> Vector Double -- ^ parameters - -> Vector Int -- ^ result -random_mp d s t p = unsafePerformIO $ do - r <- createVector $ dim p - app2 (distribution_discrete_multi_param (fromIntegral s) (fromei d) (fromIntegral t)) vec p vec r "random_mp" - return $ mapVector (\x -> (fromIntegral x) :: Int) r - --- | probability of a variate take a value outside the argument -density_mp :: MultiParamDist -- ^ density type - -> DistFunc -- ^ distribution function type - -> Vector Double -- ^ parameters - -> Vector Int -- ^ values - -> Double -- ^ result -density_mp d f p q = unsafePerformIO $ do - case d of - Multinomial -> density_only f d p q - where density_only f' d' p' q' = if f' /= Density - then error "distribution has no CDF" - else alloca $ \r -> do - app2 (distribution_dist_multi_param (fromei f') (fromei d') r) vec p vec (mapVector (\x -> (fromIntegral x) :: CUInt) q) "density_mp" - r' <- peek r - return r' - -foreign import ccall "distribution-aux.h discrete_mp" distribution_discrete_multi_param :: CInt -> CInt -> CUInt -> CInt -> Ptr Double -> CInt -> Ptr CUInt -> IO CInt -foreign import ccall "distribution-aux.h discrete_mm_dist" distribution_dist_multi_param :: CInt -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr CUInt -> IO CInt - ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ diff --git a/lib/Numeric/GSL/Fitting/Linear.hs b/lib/Numeric/GSL/Fitting/Linear.hs deleted file mode 100644 --- a/lib/Numeric/GSL/Fitting/Linear.hs +++ /dev/null @@ -1,181 +0,0 @@ -{-# OPTIONS_GHC -fglasgow-exts #-} ------------------------------------------------------------------------------ --- | --- Module : Numeric.GSL.Fitting.Linear --- Copyright : (c) Alexander Vivian Hugh McPhail 2010 --- License : GPL-style --- --- Maintainer : haskell.vivian.mcphail gmail com --- Stability : provisional --- Portability : uses ffi --- --- GSL linear regression functions --- ------------------------------------------------------------------------------ - -module Numeric.GSL.Fitting.Linear ( - linear, linear_w, linear_est, - multifit, multifit_w, multifit_est - ) where - ------------------------------------------------------------------------------ - -import Data.Packed.Vector -import Data.Packed.Matrix -import Data.Packed.Development - ---import Numeric.LinearAlgebra.Linear - ---import Control.Monad(when) - -import Foreign ---import Foreign.ForeignPtr ---import Foreign.Marshal.Alloc(alloca) -import Foreign.C.Types(CInt) ---import Foreign.C.String(newCString,peekCString) - ---import GHC.ForeignPtr (mallocPlainForeignPtrBytes) - ---import GHC.Base ---import GHC.IOBase - ---import Prelude hiding(reverse) - ------------------------------------------------------------------------------ - --- | fits the model Y = C X -linear :: Vector Double -- ^ x data - -> Vector Double -- ^ y data - -> (Double,Double,Double,Double,Double,Double) -- ^ (c_0,c_1,cov_00,cov_01,cov_11,chi_sq) -linear x y = unsafePerformIO $ do - alloca $ \c0 -> - alloca $ \c1 -> - alloca $ \chi_sq -> - alloca $ \cov00 -> - alloca $ \cov01 -> - alloca $ \cov11 -> do - app2 (fitting_linear c0 c1 chi_sq cov00 cov01 cov11) vec x vec y "linear" - c0' <- peek c0 - c1' <- peek c1 - cov00' <- peek cov00 - cov01' <- peek cov01 - cov11' <- peek cov11 - chi_sq' <- peek chi_sq - return (c0',c1',cov00',cov01',cov11',chi_sq') - - ------------------------------------------------------------------------------ - -foreign import ccall "fitting-aux.h linear" fitting_linear :: Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt - ------------------------------------------------------------------------------ - --- | fits the model Y = C X, with x data weighted -linear_w :: Vector Double -- ^ x data - -> Vector Double -- ^ x weights - -> Vector Double -- ^ y data - -> (Double,Double,Double,Double,Double,Double) -- ^ (c_0,c_1,cov_00,cov_01,cov_11,chi_sq) -linear_w x w y = unsafePerformIO $ do - alloca $ \c0 -> - alloca $ \c1 -> - alloca $ \chi_sq -> - alloca $ \cov00 -> - alloca $ \cov01 -> - alloca $ \cov11 -> do - app3 (fitting_linear_w c0 c1 chi_sq cov00 cov01 cov11) vec x vec w vec y "linear_w" - c0' <- peek c0 - c1' <- peek c1 - cov00' <- peek cov00 - cov01' <- peek cov01 - cov11' <- peek cov11 - chi_sq' <- peek chi_sq - return (c0',c1',cov00',cov01',cov11',chi_sq') - ------------------------------------------------------------------------------ - -foreign import ccall "fitting-aux.h linear_weighted" fitting_linear_w :: Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt - ------------------------------------------------------------------------------ - --- | computes the fitted function and standard deviation at the input point -linear_est :: Double -- ^ x data point - -> Double -- ^ c0 - -> Double -- ^ c1 - -> Double -- ^ cov00 - -> Double -- ^ cov01 - -> Double -- ^ cov11 - -> (Double,Double) -- ^ (y,error) -linear_est x c0 c1 cov00 cov01 cov11 = unsafePerformIO $ do - alloca $ \y -> - alloca $ \e -> do - check "linear_est" $ fitting_linear_est x c0 c1 cov00 cov01 cov11 y e - y' <- peek y - e' <- peek e - return (y',e') - ------------------------------------------------------------------------------ - -foreign import ccall "fitting-aux.h linear_estimate" fitting_linear_est :: Double -> Double -> Double -> Double -> Double -> Double -> Ptr Double -> Ptr Double -> IO CInt - ------------------------------------------------------------------------------ - --- | fit the model Y = C X, with design matrix X --- | X is a design matrix X_{ij} = x_j(i) with i observations and p predictors --- | a polynomial would be X_{ij} = x_i^j --- | a fourier series would be X_{ij} = sin (\omega_j x_i) -multifit :: Matrix Double -- ^ design matrix (X) - -> Vector Double -- ^ observations - -> (Vector Double,Matrix Double,Double) -- ^ (coefficients,covariance,chi_sq) -multifit x y = unsafePerformIO $ do - let ys = dim y - cov <- createMatrix RowMajor ys ys - p <- createVector ys - alloca$ \chi_sq -> do - app4 (fitting_multifit chi_sq) mat x vec y vec p mat cov "multifit" - chi_sq' <- peek chi_sq - return (p,cov,chi_sq') - ------------------------------------------------------------------------------ - -foreign import ccall "fitting_aux.h multifit" fitting_multifit :: Ptr Double -> CInt -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt - ------------------------------------------------------------------------------ - --- | fit the model Y = C X, with design matrix X, and x weighted -multifit_w :: Matrix Double -- ^ design matrix (X) - -> Vector Double -- ^ weights - -> Vector Double -- ^ observations - -> (Vector Double,Matrix Double,Double) -- ^ (coefficients,covariance,chi_sq) -multifit_w x w y = unsafePerformIO $ do - let ys = dim y - cov <- createMatrix RowMajor ys ys - p <- createVector ys - alloca$ \chi_sq -> do - app5 (fitting_multifit_w chi_sq) mat x vec w vec y vec p mat cov "multifit" - chi_sq' <- peek chi_sq - return (p,cov,chi_sq') - ------------------------------------------------------------------------------ - -foreign import ccall "fitting_aux.h multifit_weighted" fitting_multifit_w :: Ptr Double -> CInt -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt - ------------------------------------------------------------------------------ - --- | computes the fitted function and standard deviation at the input point -multifit_est :: Vector Double -- ^ input point - -> Vector Double -- ^ the coefficients - -> Matrix Double -- ^ the covariance matrix - -> (Double,Double) -- ^ (y,y_error_ -multifit_est x c cov = unsafePerformIO $ do - alloca $ \y -> - alloca $ \e -> do - app3 (fitting_multifit_est y e) vec x vec c mat cov "multifit_estimate" - y' <- peek y - e' <- peek e - return (y',e') - ------------------------------------------------------------------------------ - -foreign import ccall "fitting_aux.h multifit_estimate" fitting_multifit_est :: Ptr Double -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt - ------------------------------------------------------------------------------ diff --git a/lib/Numeric/GSL/Histogram.hs b/lib/Numeric/GSL/Histogram.hs deleted file mode 100644 --- a/lib/Numeric/GSL/Histogram.hs +++ /dev/null @@ -1,468 +0,0 @@ -{-# OPTIONS_GHC -fglasgow-exts #-} ------------------------------------------------------------------------------ --- | --- Module : Numeric.GSL.Histogram --- Copyright : (c) Alexander Vivian Hugh McPhail 2010 --- License : GPL-style --- --- Maintainer : haskell.vivian.mcphail gmail com --- Stability : provisional --- Portability : uses ffi --- --- GSL histogram functions --- ------------------------------------------------------------------------------ - -module Numeric.GSL.Histogram ( - Histogram - , fromRanges, fromLimits - , addList, addVector, addListWeighted, addVectorWeighted - , toVectors - , getBin, getRange - , getMax, getMin, getBins - , find - , maxVal, maxBin, minVal, minBin - , mean, stddev, sum - , equalBins - , add, subtract, multiply, divide, shift, scale - , fwriteHistogram, freadHistogram, fprintfHistogram, fscanfHistogram - -- - , HistogramPDF - , fromHistogram - , sample - ) where - ------------------------------------------------------------------------------ - -import Data.Packed.Vector ---import Data.Packed.Matrix hiding(toLists) -import Data.Packed.Development - ---import Numeric.LinearAlgebra.Linear - ---import Control.Monad - -import Foreign hiding(shift) ---import Foreign.ForeignPtr ---import Foreign.Marshal.Alloc(alloca) -import Foreign.C.Types(CInt,CChar) ---import Foreign.C.String(newCString,peekCString) -import Foreign.C.String(newCString) ---import Control.Monad(when) - ---import GHC.ForeignPtr (mallocPlainForeignPtrBytes) - ---import GHC.Base ---import GHC.IOBase - -import Prelude hiding(subtract,sum) - ------------------------------------------------------------------------------ - -data Hist -type HistHandle = Ptr Hist --- | A histogram structure -data Histogram = H { hdim :: {-# UNPACK #-} !Int -- ^ number of bins - , hist :: {-# UNPACK #-} !(ForeignPtr Hist) } - -data PDF -type PDFHandle = Ptr PDF --- | A histogram-derived cumulative distribution function (CDF) -data HistogramPDF = P { pdf :: {-# UNPACK #-} !(ForeignPtr PDF)} - ------------------------------------------------------------------------------ - -instance Eq Histogram where - (==) = equalBins -{- -instance Num Histogram where - (+) = add - (-) = subtract - negate = flip scale (-1.0) - (*) = multiply - signum = error "can't signumm Histogram" - abs = error "can't abs Histogram" - fromInteger x = fromLimits (fromInteger x) (0,1) - -instance Fractional Histogram where - fromRational x = fromLimits (round x) (0,fromRational x) - (/) = divide --} ------------------------------------------------------------------------------ - -foreign import ccall "gsl-histogram.h gsl_histogram_alloc" histogram_new :: CInt -> IO HistHandle -foreign import ccall "gsl-histogram.h &gsl_histogram_free" histogram_free :: FunPtr (HistHandle -> IO ()) - ------------------------------------------------------------------------------ - --- | create a histogram with n bins from ranges (x0->x1),(x1->x2)..(xn->xn+1) -fromRangesIO :: Vector Double -> IO Histogram -fromRangesIO v = do - let sz = fromIntegral $ dim v - 1 - h <- histogram_new sz - h' <- newForeignPtr histogram_free h - app1 (\d p -> withForeignPtr h' $ \f -> histogram_set_ranges f (fromIntegral d) p) vec v "fromRanges" - return $ H (fromIntegral sz) h' - --- | create a histogram with n bins and lower and upper limits -fromLimitsIO :: Int -> (Double,Double) -> IO Histogram -fromLimitsIO n (l,u) = do - h <- histogram_new (fromIntegral n) - h' <- newForeignPtr histogram_free h - check "set_ranges_uniform" $ withForeignPtr h' (\f -> histogram_set_ranges_uniform f l u) - return $ H n h' - -foreign import ccall "gsl-histogram.h gsl_histogram_set_ranges" histogram_set_ranges :: HistHandle -> CInt -> Ptr Double -> IO CInt -foreign import ccall "gsl-histogram.h gsl_histogram_set_ranges_uniform" histogram_set_ranges_uniform :: HistHandle -> Double -> Double -> IO CInt - ------------------------------------------------------------------------------ - --- | create a histogram with n bins from ranges (x0->x1),(x1->x2)..(xn->xn+1) and increment from a vector -fromRanges :: Vector Double -- ^ the ranges - -> Vector Double -- ^ the data - -> Histogram -- ^ result -fromRanges r d = unsafePerformIO $ do - h <- fromRangesIO r - incrementVectorIO h d - return h - --- | create a histogram with n bins and lower and upper limits and increment from a vector -fromLimits :: Int -- ^ bins - -> (Double,Double) -- ^ lower and upper limits - -> Vector Double -- ^ the data - -> Histogram -- ^ result -fromLimits n r d = unsafePerformIO $ do - h <- fromLimitsIO n r - incrementVectorIO h d - return h - ------------------------------------------------------------------------------ - --- | extract the ranges and bins -toVectors :: Histogram -> (Vector Double,Vector Double) -- ^ (ranges,bins) -toVectors (H b h) = unsafePerformIO $ do - rs <- createVector (b+1) - bs <- createVector b - app2 (\s1 p1 s2 p2 -> withForeignPtr h $ \f -> histogram_to_vectors f s1 p1 s2 p2) vec rs vec bs "toVectors" - return (rs,bs) - -foreign import ccall "gsl-histogram.h to_vectors" histogram_to_vectors :: HistHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt - ------------------------------------------------------------------------------ - --- | create a copy of a histogram -cloneHistogram :: Histogram -> IO Histogram -cloneHistogram (H b h) = do - h' <- withForeignPtr h histogram_clone - h'' <- newForeignPtr histogram_free h' - return $ H b h'' - -foreign import ccall "gsl-histogram.h gsl_histogram_clone" histogram_clone :: HistHandle -> IO HistHandle - ------------------------------------------------------------------------------ - --- | adds 1.0 to the correct bin for each element of the list -addList :: Histogram -> [Double] -> Histogram -addList h xs = unsafePerformIO $ do - h' <- cloneHistogram h - incrementListIO h' xs - return h' - --- | adds 1.0 to the correct bin for each element of the vector -addVector :: Histogram -> Vector Double -> Histogram -addVector h v = unsafePerformIO $ do - h' <- cloneHistogram h - incrementVectorIO h' v - return h' - --- | adds the appropriate weight for each element of the list -addListWeighted :: Histogram -> [(Double,Double)] -> Histogram -addListWeighted h xs = unsafePerformIO $ do - h' <- cloneHistogram h - accumulateListIO h' xs - return h' - --- | adds the appropriate weight for each element of the list -addVectorWeighted :: Histogram -> Vector Double -> Vector Double -> Histogram -addVectorWeighted h v w = unsafePerformIO $ do - h' <- cloneHistogram h - accumulateVectorIO h' v w - return h' - --- | add 1.0 to the correct bin, fails silently if the value is outside the range -incrementIO :: Histogram -> Double -> IO () -incrementIO (H _ h) x = withForeignPtr h (\f -> check "increment" $ histogram_increment f x) - --- | add 1.0 to the correct bin for each element of the vector, fails silently if the value is outside the range -incrementListIO :: Histogram -> [Double] -> IO () -incrementListIO (H _ h) xs = withForeignPtr h (\f -> mapM_ (histogram_increment f) xs) - --- | add 1.0 to the correct bin for each element of the vector, fails silently if the value is outside the range -incrementVectorIO :: Histogram -> Vector Double -> IO () -incrementVectorIO (H _ h) v = do - app1 (\s p -> withForeignPtr h (\f -> histogram_increment_vector f s p)) vec v "incrementVector" - return () - --- | adds the weight (second Double) to the bin appropriate for the value (first Double) -accumulateIO :: Histogram -> Double -> Double -> IO () -accumulateIO (H _ h) x w = withForeignPtr h (\f -> check "accumulate" $ histogram_accumulate f x w) - --- | add the weight (snd) to the correct bin for each (fst) element of the list, fails silently if the value is outside the range -accumulateListIO :: Histogram -> [(Double,Double)] -> IO () -accumulateListIO (H _ h) xs = do - withForeignPtr h (\f -> mapM_ (\(x,w) -> histogram_accumulate f x w) xs) - return () - --- | add the weight (second vector) to the correct bin for each element of the first vector, fails silently if the value is outside the range -accumulateVectorIO :: Histogram -> Vector Double -> Vector Double -> IO () -accumulateVectorIO (H _ h) v w = do - app2 (\s1 p1 s2 p2 -> withForeignPtr h (\f -> histogram_accumulate_vector f s1 p1 s2 p2)) vec v vec w "accumulateVector" - return () - --- | returns the contents of the i-th bin -getBin :: Histogram -> Int -> Double -getBin (H _ h) b = unsafePerformIO $ do - withForeignPtr h (\f -> histogram_get f (fromIntegral b)) - --- | returns the upper and lower limits of the i-th bin -getRange :: Histogram -> Int -> (Double,Double) -getRange (H _ h) b = unsafePerformIO $ do - alloca $ \l -> - alloca $ \u -> do - check "get_range" $ withForeignPtr h (\f -> histogram_get_range f (fromIntegral b) l u) - l' <- peek l - u' <- peek u - return (l',u') - --- | the maximum upper range limit -getMax :: Histogram -> Double -getMax (H _ h) = unsafePerformIO $ withForeignPtr h histogram_max - --- | the minimum lower range limit -getMin :: Histogram -> Double -getMin (H _ h) = unsafePerformIO $ withForeignPtr h histogram_min - --- | the number of bins -getBins :: Histogram -> Int -getBins (H b _) = b - --- | reset all the bins to zero -reset :: Histogram -> IO () -reset (H _ h) = withForeignPtr h histogram_reset - -foreign import ccall "gsl-histogram.h gsl_histogram_increment" histogram_increment :: HistHandle -> Double -> IO CInt -foreign import ccall "gsl-histogram.h gsl_histogram_accumulate" histogram_accumulate :: HistHandle -> Double -> Double -> IO CInt -foreign import ccall "gsl-histogram.h gsl_histogram_get" histogram_get :: HistHandle -> CInt -> IO Double -foreign import ccall "gsl-histogram.h gsl_histogram_get_range" histogram_get_range :: HistHandle -> CInt -> Ptr Double -> Ptr Double -> IO CInt -foreign import ccall "gsl-histogram.h gsl_histogram_max" histogram_max :: HistHandle -> IO Double -foreign import ccall "gsl-histogram.h gsl_histogram_min" histogram_min :: HistHandle -> IO Double -foreign import ccall "gsl-histogram.h gsl_histogram_bins" histogram_bins :: HistHandle -> IO Int -foreign import ccall "gsl-histogram.h gsl_histogram_reset" histogram_reset :: HistHandle -> IO () - -foreign import ccall "histogram-aux.h increment_vector" histogram_increment_vector :: HistHandle -> CInt -> Ptr Double -> IO CInt -foreign import ccall "histogram-aux.h accumulate_vector" histogram_accumulate_vector :: HistHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt - ------------------------------------------------------------------------------ - --- | find the bin corresponding to the value -find :: Histogram -> Double -> Maybe Int -find (H _ h) x = unsafePerformIO $ do - alloca $ \b -> do - err <- withForeignPtr h (\f -> histogram_find f x b) - if err == 0 - then do - b' <- peek b - return $ Just $ fromIntegral b' - else return Nothing - -foreign import ccall "gsl-histogram.h gsl_histogram_find" histogram_find :: HistHandle -> Double -> Ptr CInt -> IO CInt - ------------------------------------------------------------------------------ - --- | the maximum value contained in the bins -maxVal :: Histogram -> Double -maxVal (H _ h) = unsafePerformIO $ withForeignPtr h histogram_max_val - --- | the index of the bin containing the maximum value -maxBin :: Histogram -> Int -maxBin (H _ h) = unsafePerformIO $ do - i <- withForeignPtr h histogram_max_bin - return $ fromIntegral i - --- | the minimum value contained in the bins -minVal :: Histogram -> Double -minVal (H _ h) = unsafePerformIO $ withForeignPtr h histogram_min_val - --- | the index of the bin containing the minimum value -minBin :: Histogram -> Int -minBin (H _ h) = unsafePerformIO $ do - i <- withForeignPtr h histogram_min_bin - return $ fromIntegral i - --- | the mean of the values, accuracy limited by bin width -mean :: Histogram -> Double -mean (H _ h) = unsafePerformIO $ withForeignPtr h histogram_mean - --- | the standard deviation of the values, accuracy limited by bin width -stddev :: Histogram -> Double -stddev (H _ h) = unsafePerformIO $ withForeignPtr h histogram_sigma - --- | the sum of the values, accuracy limited by bin width -sum :: Histogram -> Double -sum (H _ h) = unsafePerformIO $ withForeignPtr h histogram_sum - -foreign import ccall "gsl-histogram.h gsl_histogram_max_val" histogram_max_val :: HistHandle -> IO Double -foreign import ccall "gsl-histogram.h gsl_histogram_max_bin" histogram_max_bin :: HistHandle -> IO CInt -foreign import ccall "gsl-histogram.h gsl_histogram_min_val" histogram_min_val :: HistHandle -> IO Double -foreign import ccall "gsl-histogram.h gsl_histogram_min_bin" histogram_min_bin :: HistHandle -> IO CInt -foreign import ccall "gsl-histogram.h gsl_histogram_mean" histogram_mean :: HistHandle -> IO Double -foreign import ccall "gsl-histogram.h gsl_histogram_sigma" histogram_sigma :: HistHandle -> IO Double -foreign import ccall "gsl-histogram.h gsl_histogram_sum" histogram_sum :: HistHandle -> IO Double - ------------------------------------------------------------------------------ - --- | returns True of all the individual bin ranges of the two histograms are identical -equalBins :: Histogram -> Histogram -> Bool -equalBins (H _ h1) (H _ h2) = unsafePerformIO $ do - i <- withForeignPtr h1 $ \p1 -> do - withForeignPtr h2 $ \p2 -> histogram_equal_bins p1 p2 - if (fromIntegral i) == (1 :: Int) - then return True - else return False - --- | adds the contents of the bins of the second histogram to the first -add :: Histogram -> Histogram -> Histogram -add d (H _ s) = unsafePerformIO $ do - h@(H _ d') <- cloneHistogram d - check "add" $ - withForeignPtr d' $ \dp -> - withForeignPtr s $ \sp -> histogram_add dp sp - return h - --- | subtracts the contents of the bins of the second histogram from the first -subtract :: Histogram -> Histogram -> Histogram -subtract d (H _ s) = unsafePerformIO $ do - h@(H _ d') <- cloneHistogram d - check "subtract" $ - withForeignPtr d' $ \dp -> - withForeignPtr s $ \sp -> histogram_sub dp sp - return h - --- | multiplies the contents of the bins of the second histogram by the first -multiply :: Histogram -> Histogram -> Histogram -multiply d (H _ s) = unsafePerformIO $ do - h@(H _ d') <- cloneHistogram d - check "multiply" $ - withForeignPtr d' $ \dp -> - withForeignPtr s $ \sp -> histogram_mul dp sp - return h - --- | divides the contents of the bins of the first histogram by the second -divide :: Histogram -> Histogram -> Histogram -divide d (H _ s) = unsafePerformIO $ do - h@(H _ d') <- cloneHistogram d - check "divide" $ - withForeignPtr d' $ \dp -> - withForeignPtr s $ \sp -> histogram_div dp sp - return h - --- | multiplies the contents of the bins by a constant -scale :: Histogram -> Double -> Histogram -scale d s = unsafePerformIO $ do - h@(H _ d') <- cloneHistogram d - check "scale" $ - withForeignPtr d' $ (\f -> histogram_scale f s) - return h - --- | adds a constant to the contents of the bins -shift :: Histogram -> Double -> Histogram -shift d s = unsafePerformIO $ do - h@(H _ d') <- cloneHistogram d - check "shift" $ - withForeignPtr d' $ (\f -> histogram_shift f s) - return h - -foreign import ccall "gsl-histogram.h gsl_histogram_equal_bins_p" histogram_equal_bins :: HistHandle -> HistHandle -> IO CInt -foreign import ccall "gsl-histogram.h gsl_histogram_add" histogram_add :: HistHandle -> HistHandle -> IO CInt -foreign import ccall "gsl-histogram.h gsl_histogram_sub" histogram_sub :: HistHandle -> HistHandle -> IO CInt -foreign import ccall "gsl-histogram.h gsl_histogram_mul" histogram_mul :: HistHandle -> HistHandle -> IO CInt -foreign import ccall "gsl-histogram.h gsl_histogram_div" histogram_div :: HistHandle -> HistHandle -> IO CInt -foreign import ccall "gsl-histogram.h gsl_histogram_scale" histogram_scale :: HistHandle -> Double -> IO CInt -foreign import ccall "gsl-histogram.h gsl_histogram_shift" histogram_shift :: HistHandle -> Double -> IO CInt - ------------------------------------------------------------------------------ - --- | write a histogram in the native binary format (may not be portable) -fwriteHistogram :: FilePath -> Histogram -> IO () -fwriteHistogram fn (H _ h) = do - cn <- newCString fn - check "fwriteHistogram" $ - withForeignPtr h $ histogram_fwrite cn - free cn - --- | read a histogram in the native binary format, number of bins must be known -freadHistogram :: FilePath -> Int -> IO Histogram -freadHistogram fn b = do - h <- histogram_new (fromIntegral b) - h' <- newForeignPtr histogram_free h - cn <- newCString fn - check "freadHistogram" $ - withForeignPtr h' $ histogram_fread cn - return $ H b h' - --- | saves the histogram with the given formats (%f,%e,%g) for ranges and bins --- each line comprises: range[i] range[i+1] bin[i] -fprintfHistogram :: FilePath -> String -> String -> Histogram -> IO () -fprintfHistogram fn fr fb (H _ h) = do - cn <- newCString fn - cr <- newCString fr - cb <- newCString fb - check "fprintfHistogram" $ - withForeignPtr h $ histogram_fprintf cn cr cb - free cn - free cr - free cb - return () - --- | reads formatted data as written by fprintf, the number of bins must be known in advance -fscanfHistogram :: FilePath -> Int -> IO Histogram -fscanfHistogram fn b = do - h <- histogram_new (fromIntegral b) - h' <- newForeignPtr histogram_free h - cn <- newCString fn - check "fscanfHistogram" $ - withForeignPtr h' $ histogram_fscanf cn - return $ H b h' - -foreign import ccall "histogram-aux.h hist_fwrite" histogram_fwrite :: Ptr CChar -> HistHandle -> IO CInt -foreign import ccall "histogram-aux.h hist_fread" histogram_fread :: Ptr CChar -> HistHandle -> IO CInt -foreign import ccall "histogram-aux.h hist_fprintf" histogram_fprintf :: Ptr CChar -> Ptr CChar -> Ptr CChar -> HistHandle -> IO CInt -foreign import ccall "histogram-aux.h hist_fscanf" histogram_fscanf :: Ptr CChar -> HistHandle -> IO CInt ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ - -foreign import ccall "gsl-histogram.h gsl_histogram_pdf_alloc" histogram_pdf_new :: CInt -> IO PDFHandle -foreign import ccall "gsl-histogram.h &gsl_histogram_pdf_free" histogram_pdf_free :: FunPtr (PDFHandle -> IO ()) - ------------------------------------------------------------------------------ - --- | create a histogram PDF from a histogram -fromHistogram :: Histogram -> HistogramPDF -fromHistogram (H b h) = unsafePerformIO $ do - p <- histogram_pdf_new $ fromIntegral b - p' <- newForeignPtr histogram_pdf_free p - withForeignPtr p' $ \p'' -> - withForeignPtr h $ \h' -> check "pdf_init" $ histogram_pdf_init p'' h' - return $ P p' - --- | given a randomm from the uniform distribution [0,1], draw a random sample from the PDF -sample :: HistogramPDF -> Double -sample (P p) = unsafePerformIO $ withForeignPtr p $ \p' -> histogram_pdf_sample p' - -foreign import ccall "gsl-histogram.h gsl_histogram_pdf_init" histogram_pdf_init :: PDFHandle -> HistHandle -> IO CInt -foreign import ccall "gsl-histogram.h gsl_histogram_pdf_sample" histogram_pdf_sample :: PDFHandle -> IO Double - ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ diff --git a/lib/Numeric/GSL/Histogram2D.hs b/lib/Numeric/GSL/Histogram2D.hs deleted file mode 100644 --- a/lib/Numeric/GSL/Histogram2D.hs +++ /dev/null @@ -1,534 +0,0 @@ -{-# OPTIONS_GHC -fglasgow-exts #-} ------------------------------------------------------------------------------ --- | --- Module : Numeric.GSL.Histogram2D --- Copyright : (c) Alexander Vivian Hugh McPhail 2010 --- License : GPL-style --- --- Maintainer : haskell.vivian.mcphail gmail com --- Stability : provisional --- Portability : uses ffi --- --- GSL 2D histogram functions --- ------------------------------------------------------------------------------ - -module Numeric.GSL.Histogram2D ( - Histogram2D - , fromRanges, fromLimits - , addList, addVector, addListWeighted, addVectorWeighted - , toMatrix - , getBin, getXRange, getYRange - , getXMax, getYMax, getXMin, getYMin, getXBins, getYBins - , reset - , find - , maxVal, maxBin, minVal, minBin - , xmean, ymean, xstddev, ystddev, covariance, sum - , equalBins - , add, subtract, multiply, divide, shift, scale - , fwriteHistogram2D, freadHistogram2D, fprintfHistogram2D, fscanfHistogram2D - -- - , Histogram2DPDF - , fromHistogram2D - , sample - ) where - ------------------------------------------------------------------------------ - -import Data.Packed.Vector -import Data.Packed.Matrix -import Data.Packed.Development - ---import Numeric.LinearAlgebra.Linear - ---import Control.Monad ---import Control.Monad(when) - -import Foreign hiding(shift) ---import Foreign.ForeignPtr ---import Foreign.Marshal.Alloc(alloca) -import Foreign.C.Types(CInt,CChar) ---import Foreign.C.String(newCString,peekCString) -import Foreign.C.String(newCString) - ---import GHC.ForeignPtr (mallocPlainForeignPtrBytes) - ---import GHC.Base ---import GHC.IOBase - -import Prelude hiding(subtract,sum) - ------------------------------------------------------------------------------ - -data Hist2D -type Hist2DHandle = Ptr Hist2D --- | A histogram structure -data Histogram2D = H { hxdim :: {-# UNPACK #-} !Int -- ^ number of bins - , hydim :: {-# UNPACK #-} !Int -- ^ number of bins - , hist :: {-# UNPACK #-} !(ForeignPtr Hist2D) } - -data PDF -type PDFHandle = Ptr PDF --- | A histogram-derived cumulative distribution function (CDF) -data Histogram2DPDF = P { pdf :: {-# UNPACK #-} !(ForeignPtr PDF)} - ------------------------------------------------------------------------------ - -instance Eq Histogram2D where - (==) = equalBins -{- -instance Num Histogram2D where - (+) = add - (-) = subtract - negate = flip scale (-1.0) - (*) = multiply - signum = error "can't signumm Histogram2D" - abs = error "can't abs Histogram2D" - fromInteger x = fromLimits (fromInteger x) (0,1) - -instance Fractional Histogram2D where - fromRational x = fromLimits (round x) (0,fromRational x) - (/) = divide --} ------------------------------------------------------------------------------ - -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_alloc" histogram2d_new :: CInt -> CInt -> IO Hist2DHandle -foreign import ccall "gsl-histogram2d.h &gsl_histogram2d_free" histogram2d_free :: FunPtr (Hist2DHandle -> IO ()) - ------------------------------------------------------------------------------ - --- | create a histogram with n bins from ranges (x0->x1),(x1->x2),..,(xn->xn+1) and y0,..,yn+1 -fromRangesIO :: Vector Double -> Vector Double -> IO Histogram2D -fromRangesIO v w = do - let sx = fromIntegral $ dim v - 1 - let sy = fromIntegral $ dim w - 1 - h <- histogram2d_new sx sy - h' <- newForeignPtr histogram2d_free h - app2 (\xs xp ys yp -> withForeignPtr h' (\f -> histogram2d_set_ranges f xp xs yp ys)) vec v vec w "fromRanges" - return $ H (fromIntegral sx) (fromIntegral sy) h' - --- | create a histogram with n bins and lower and upper limits -fromLimitsIO :: Int -> Int -- ^ number of bins - -> (Double,Double) -- ^ xmin, xmax - -> (Double,Double) -- ^ ymin, ymax - -> IO Histogram2D -fromLimitsIO nx ny (lx,ux) (uy,ly) = do - h <- histogram2d_new (fromIntegral nx) (fromIntegral ny) - h' <- newForeignPtr histogram2d_free h - check "set_ranges_uniform" $ withForeignPtr h' (\f -> histogram2d_set_ranges_uniform f lx ux ly uy) - return $ H nx ny h' - -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_set_ranges" histogram2d_set_ranges :: Hist2DHandle -> Ptr Double -> CInt -> Ptr Double -> CInt -> IO CInt -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_set_ranges_uniform" histogram2d_set_ranges_uniform :: Hist2DHandle -> Double -> Double -> Double -> Double -> IO CInt - ------------------------------------------------------------------------------ - --- | create a histogram with n bins from ranges (x0->x1),(x1->x2)..(xn->xn+1) and increment from a vector -fromRanges :: Vector Double -- ^ the x ranges - -> Vector Double -- ^ the y ranges - -> [(Double,Double)] -- ^ the data - -> Histogram2D -- ^ result -fromRanges rx ry d = unsafePerformIO $ do - h <- fromRangesIO rx ry - incrementListIO h d - return h - --- | create a histogram with n bins and lower and upper limits and increment from a vector -fromLimits :: Int -> Int -- ^ bins - -> (Double,Double) -- ^ x lower and upper limits - -> (Double,Double) -- ^ y lower and upper limits - -> [(Double,Double)] -- ^ the data - -> Histogram2D -- ^ result -fromLimits nx ny rx ry d = unsafePerformIO $ do - h <- fromLimitsIO nx ny rx ry - incrementListIO h d - return h - ------------------------------------------------------------------------------ - --- | extract the ranges and bins -toMatrix :: Histogram2D -> (Vector Double,Vector Double,Matrix Double) -- ^ (ranges,bins) -toMatrix (H bx by h) = unsafePerformIO $ do - rx <- createVector (bx+1) - ry <- createVector (by+1) - bs <- createMatrix RowMajor bx by - app3 (\s1 p1 s2 p2 sx sy p -> withForeignPtr h $ \f -> histogram_to_matrix f s1 p1 s2 p2 sx sy p) vec rx vec ry mat bs "toMatrix" - return (rx,ry,bs) - -foreign import ccall "histogram-aux.h to_matrix" histogram_to_matrix :: Hist2DHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt - ------------------------------------------------------------------------------ - --- | create a copy of a histogram -cloneHistogram2D :: Histogram2D -> IO Histogram2D -cloneHistogram2D (H nx ny h) = do - h' <- withForeignPtr h histogram2d_clone - h'' <- newForeignPtr histogram2d_free h' - return $ H nx ny h'' - -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_clone" histogram2d_clone :: Hist2DHandle -> IO Hist2DHandle - ------------------------------------------------------------------------------ - --- | add 1.0 to the correct bin for each element of the list, fails silently if the value is outside the range -addList :: Histogram2D -> [(Double,Double)] -> Histogram2D -addList h xs = unsafePerformIO $ do - h' <- cloneHistogram2D h - incrementListIO h' xs - return h' - --- | add 1.0 to the correct bin for each element of the vector pair, fails silently if the value is outside the range -addVector :: Histogram2D -> Vector Double -> Vector Double -> Histogram2D -addVector h x y = unsafePerformIO $ do - h' <- cloneHistogram2D h - incrementVectorIO h' x y - return h' - --- add the appropriate weight for each element of the list, fails silently if the value is outside the range -addListWeighted :: Histogram2D -> [(Double,Double,Double)] -> Histogram2D -addListWeighted h xs = unsafePerformIO $ do - h' <- cloneHistogram2D h - accumulateListIO h' xs - return h' - --- add the appropriate weight for each element of the vector pair, fails silently if the value is outside the range -addVectorWeighted :: Histogram2D -> Vector Double -> Vector Double -> Vector Double -> Histogram2D -addVectorWeighted h x y w = unsafePerformIO $ do - h' <- cloneHistogram2D h - accumulateVectorIO h' x y w - return h' - - ------------------------------------------------------------------------------ - --- | add 1.0 to the correct bin, fails silently if the value is outside the range -incrementIO :: Histogram2D -> Double -> Double -> IO () -incrementIO (H _ _ h) x y = do - check "increment" $ withForeignPtr h (\f -> histogram2d_increment f x y) - return () - --- | add 1.0 to the correct bin for each element of the vector pair, fails silently if the value is outside the range -incrementVectorIO :: Histogram2D -> Vector Double -> Vector Double -> IO () -incrementVectorIO (H _ _ h) x y = do - app2 (\xs xp ys yp -> withForeignPtr h (\f -> histogram2d_increment_matrix f xs xp ys yp)) vec x vec y "incrementVector" - return () - --- | add 1.0 to the correct bin for each element of the list, fails silently if the value is outside the range -incrementListIO :: Histogram2D -> [(Double,Double)] -> IO () -incrementListIO (H _ _ h) zs = withForeignPtr h (\f -> mapM_ (\(x,y) -> histogram2d_increment f x y) zs) - --- | Adds the weight (third Double) to the bin appropriate for the value (first two Doubles) -accumulateIO :: Histogram2D -> Double -> Double -> Double -> IO () -accumulateIO (H _ _ h) x y w = do - check "accumulate" $ withForeignPtr h (\f -> histogram2d_accumulate f x y w) - return () - --- | add the weight (third) to the correct bin for each vector pair element, fails silently if the value is outside the range -accumulateVectorIO :: Histogram2D -> Vector Double -> Vector Double -> Vector Double -> IO () -accumulateVectorIO (H _ _ h) x y w = do - app3 (\xs xp ys yp ws wp -> withForeignPtr h (\f -> histogram2d_accumulate_matrix f xs xp ys yp ws wp)) vec x vec y vec w "accumulateVector" - return () - --- | add the weight (snd) to the correct bin for each (fst) element of the list, fails silently if the value is outside the range -accumulateListIO :: Histogram2D -> [(Double,Double,Double)] -> IO () -accumulateListIO (H _ _ h) zs = withForeignPtr h (\f -> mapM_ (\(x,y,w) -> histogram2d_accumulate f x y w) zs) - --- | returns the contents of the i-th bin -getBin :: Histogram2D -> Int -> Int -> Double -getBin (H _ _ h) bx by = unsafePerformIO $ do - withForeignPtr h (\f -> histogram2d_get f (fromIntegral bx) (fromIntegral by)) - --- | returns the upper and lower limits in the first dimension of the i-th bin -getXRange :: Histogram2D -> Int -> (Double,Double) -getXRange (H _ _ h) b = unsafePerformIO $ do - alloca $ \l -> - alloca $ \u -> do - check "get_xrange" $ withForeignPtr h (\f -> histogram2d_get_xrange f (fromIntegral b) l u) - l' <- peek l - u' <- peek u - return (l',u') - --- | returns the upper and lower limits in the second dimension of the i-th bin -getYRange :: Histogram2D -> Int -> (Double,Double) -getYRange (H _ _ h) b = unsafePerformIO $ do - alloca $ \l -> - alloca $ \u -> do - check "get_yrange" $ withForeignPtr h (\f -> histogram2d_get_yrange f (fromIntegral b) l u) - l' <- peek l - u' <- peek u - return (l',u') - --- | the maximum upper range limit in the first dimension -getXMax :: Histogram2D -> Double -getXMax (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_xmax - --- | the minimum lower range limit in the first dimension -getXMin :: Histogram2D -> Double -getXMin (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_xmin - --- | the number of binsin the first dimension -getXBins :: Histogram2D -> Int -getXBins (H bx _ _) = bx - --- | the maximum upper range limit in the first dimension -getYMax :: Histogram2D -> Double -getYMax (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_ymax - --- | the minimum lower range limit in the first dimension -getYMin :: Histogram2D -> Double -getYMin (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_ymin - --- | the number of binsin the first dimension -getYBins :: Histogram2D -> Int -getYBins (H _ by _) = by - --- | reset all the bins to zero -reset :: Histogram2D -> IO () -reset (H _ _ h) = withForeignPtr h histogram2d_reset - -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_increment" histogram2d_increment :: Hist2DHandle -> Double -> Double -> IO CInt -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_accumulate" histogram2d_accumulate :: Hist2DHandle -> Double -> Double -> Double -> IO CInt -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_get" histogram2d_get :: Hist2DHandle -> CInt -> CInt -> IO Double -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_get_xrange" histogram2d_get_xrange :: Hist2DHandle -> CInt -> Ptr Double -> Ptr Double -> IO CInt -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_get_yrange" histogram2d_get_yrange :: Hist2DHandle -> CInt -> Ptr Double -> Ptr Double -> IO CInt -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_xmax" histogram2d_xmax :: Hist2DHandle -> IO Double -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_xmin" histogram2d_xmin :: Hist2DHandle -> IO Double -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_nx" histogram2d_xn :: Hist2DHandle -> IO Int -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_ymax" histogram2d_ymax :: Hist2DHandle -> IO Double -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_ymin" histogram2d_ymin :: Hist2DHandle -> IO Double -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_ny" histogram2d_yn :: Hist2DHandle -> IO Int -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_reset" histogram2d_reset :: Hist2DHandle -> IO () - -foreign import ccall "histogram-aux.h increment_matrix" histogram2d_increment_matrix :: Hist2DHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt -foreign import ccall "histogram-aux.h accumulate_matrix" histogram2d_accumulate_matrix :: Hist2DHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt - ------------------------------------------------------------------------------ - --- | find the bin corresponding to the value -find :: Histogram2D -> (Double,Double) -> Maybe (Int,Int) -find (H _ _ h) (x,y) = unsafePerformIO $ do - alloca $ \bx -> - alloca $ \by -> do - err <- withForeignPtr h (\f -> histogram2d_find f x y bx by) - if err == 0 - then do - bx' <- peek bx - by' <- peek by - return $ Just $ (fromIntegral bx',fromIntegral by') - else return Nothing - -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_find" histogram2d_find :: Hist2DHandle -> Double -> Double -> Ptr CInt -> Ptr CInt -> IO CInt - ------------------------------------------------------------------------------ - --- | the maximum value contained in the bins -maxVal :: Histogram2D -> Double -maxVal (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_max_val - --- | the index of the bin containing the maximum value -maxBin :: Histogram2D -> (Int,Int) -maxBin (H _ _ h) = unsafePerformIO $ do - alloca $ \bx -> - alloca $ \by -> do - withForeignPtr h (\f -> histogram2d_max_bin f bx by) - bx' <- peek bx - by' <- peek by - return $ (fromIntegral bx',fromIntegral by') - --- | the minimum value contained in the bins -minVal :: Histogram2D -> Double -minVal (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_min_val - --- | the index of the bin containing the minimum value -minBin :: Histogram2D -> (Int,Int) -minBin (H _ _ h) = unsafePerformIO $ do - alloca $ \bx -> - alloca $ \by -> do - withForeignPtr h (\f -> histogram2d_min_bin f bx by) - bx' <- peek bx - by' <- peek by - return $ (fromIntegral bx',fromIntegral by') - --- | the mean of the values in the first dimension, accuracy limited by bin width -xmean :: Histogram2D -> Double -xmean (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_xmean - --- | the mean of the values in the second dimension, accuracy limited by bin width -ymean :: Histogram2D -> Double -ymean (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_ymean - --- | the standard deviation of the values in thee first dimension, accuracy limited by bin width -xstddev :: Histogram2D -> Double -xstddev (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_xsigma - --- | the standard deviation of the values in thee first dimension, accuracy limited by bin width -ystddev :: Histogram2D -> Double -ystddev (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_ysigma - --- | the covariance of the first and second dimensions -covariance :: Histogram2D -> Double -covariance (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_cov - --- | the sum of the values, accuracy limited by bin width -sum :: Histogram2D -> Double -sum (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_sum - -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_max_val" histogram2d_max_val :: Hist2DHandle -> IO Double -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_max_bin" histogram2d_max_bin :: Hist2DHandle -> Ptr CInt -> Ptr CInt -> IO () -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_min_val" histogram2d_min_val :: Hist2DHandle -> IO Double -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_min_bin" histogram2d_min_bin :: Hist2DHandle -> Ptr CInt -> Ptr CInt -> IO () -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_xmean" histogram2d_xmean :: Hist2DHandle -> IO Double -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_ymean" histogram2d_ymean :: Hist2DHandle -> IO Double -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_xsigma" histogram2d_xsigma :: Hist2DHandle -> IO Double -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_ysigma" histogram2d_ysigma :: Hist2DHandle -> IO Double -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_cov" histogram2d_cov :: Hist2DHandle -> IO Double -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_sum" histogram2d_sum :: Hist2DHandle -> IO Double - ------------------------------------------------------------------------------ - --- | returns True of all the individual bin ranges of the two histograms are identical -equalBins :: Histogram2D -> Histogram2D -> Bool -equalBins (H _ _ h1) (H _ _ h2) = unsafePerformIO $ do - i <- withForeignPtr h1 $ \p1 -> do - withForeignPtr h2 $ \p2 -> histogram2d_equal_bins p1 p2 - if (fromIntegral i) == (1 :: Int) - then return True - else return False - --- | adds the contents of the bins of the second histogram to the first -add :: Histogram2D -> Histogram2D -> Histogram2D -add d (H _ _ s) = unsafePerformIO $ do - h@(H _ _ d') <- cloneHistogram2D d - check "add" $ - withForeignPtr d' $ \dp -> - withForeignPtr s $ \sp -> histogram2d_add dp sp - return h - --- | subtracts the contents of the bins of the second histogram from the first -subtract :: Histogram2D -> Histogram2D -> Histogram2D -subtract d (H _ _ s) = unsafePerformIO $ do - h@(H _ _ d') <- cloneHistogram2D d - check "subtract" $ - withForeignPtr d' $ \dp -> - withForeignPtr s $ \sp -> histogram2d_sub dp sp - return h - --- | multiplies the contents of the bins of the second histogram by the first -multiply :: Histogram2D -> Histogram2D -> Histogram2D -multiply d (H _ _ s) = unsafePerformIO $ do - h@(H _ _ d') <- cloneHistogram2D d - check "multiply" $ - withForeignPtr d' $ \dp -> - withForeignPtr s $ \sp -> histogram2d_mul dp sp - return h - --- | divides the contents of the bins of the first histogram by the second -divide :: Histogram2D -> Histogram2D -> Histogram2D -divide d (H _ _ s) = unsafePerformIO $ do - h@(H _ _ d') <- cloneHistogram2D d - check "divide" $ - withForeignPtr d' $ \dp -> - withForeignPtr s $ \sp -> histogram2d_div dp sp - return h - --- | multiplies the contents of the bins by a constant -scale :: Histogram2D -> Double -> Histogram2D -scale d s = unsafePerformIO $ do - h@(H _ _ d') <- cloneHistogram2D d - check "scale" $ - withForeignPtr d' $ (\f -> histogram2d_scale f s) - return h - --- | adds a constant to the contents of the bins -shift :: Histogram2D -> Double -> Histogram2D -shift d s = unsafePerformIO $ do - h@(H _ _ d') <- cloneHistogram2D d - check "shift" $ - withForeignPtr d' $ (\f -> histogram2d_shift f s) - return h - -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_equal_bins_p" histogram2d_equal_bins :: Hist2DHandle -> Hist2DHandle -> IO CInt -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_add" histogram2d_add :: Hist2DHandle -> Hist2DHandle -> IO CInt -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_sub" histogram2d_sub :: Hist2DHandle -> Hist2DHandle -> IO CInt -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_mul" histogram2d_mul :: Hist2DHandle -> Hist2DHandle -> IO CInt -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_div" histogram2d_div :: Hist2DHandle -> Hist2DHandle -> IO CInt -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_scale" histogram2d_scale :: Hist2DHandle -> Double -> IO CInt -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_shift" histogram2d_shift :: Hist2DHandle -> Double -> IO CInt - ------------------------------------------------------------------------------ - --- | write a histogram in the native binary format (may not be portable) -fwriteHistogram2D :: FilePath -> Histogram2D -> IO () -fwriteHistogram2D fn (H _ _ h) = do - cn <- newCString fn - check "fwriteHistogram2d2D" $ - withForeignPtr h $ histogram2d_fwrite cn - free cn - --- | read a histogram in the native binary format, number of bins must be known -freadHistogram2D :: FilePath -> Int -> Int -> IO Histogram2D -freadHistogram2D fn bx by = do - h <- histogram2d_new (fromIntegral bx) (fromIntegral by) - h' <- newForeignPtr histogram2d_free h - cn <- newCString fn - check "freadHistogram2d2D" $ - withForeignPtr h' $ histogram2d_fread cn - return $ H bx by h' - --- | saves the histogram with the given formats (%f,%e,%g) for ranges and bins --- each line comprises: xrange[i] xrange[i+1] xrange[j] xrange[j+1] bin(i,j) -fprintfHistogram2D :: FilePath -> String -> String -> Histogram2D -> IO () -fprintfHistogram2D fn fr fb (H _ _ h) = do - cn <- newCString fn - cr <- newCString fr - cb <- newCString fb - check "fprintfHistogram2d2D" $ - withForeignPtr h $ histogram2d_fprintf cn cr cb - free cn - free cr - free cb - return () - --- | reads formatted data as written by fprintf, the number of bins must be known in advance -fscanfHistogram2D :: FilePath -> Int -> Int -> IO Histogram2D -fscanfHistogram2D fn bx by = do - h <- histogram2d_new (fromIntegral bx) (fromIntegral by) - h' <- newForeignPtr histogram2d_free h - cn <- newCString fn - check "fscanfHistogram2d2D" $ - withForeignPtr h' $ histogram2d_fscanf cn - return $ H bx by h' - -foreign import ccall "histogram-aux.h hist2d_fwrite" histogram2d_fwrite :: Ptr CChar -> Hist2DHandle -> IO CInt -foreign import ccall "histogram-aux.h hist2d_fread" histogram2d_fread :: Ptr CChar -> Hist2DHandle -> IO CInt -foreign import ccall "histogram-aux.h hist2d_fprintf" histogram2d_fprintf :: Ptr CChar -> Ptr CChar -> Ptr CChar -> Hist2DHandle -> IO CInt -foreign import ccall "histogram-aux.h hist2d_fscanf" histogram2d_fscanf :: Ptr CChar -> Hist2DHandle -> IO CInt ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ - -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_pdf_alloc" histogram2d_pdf_new :: CInt -> CInt -> IO PDFHandle -foreign import ccall "gsl-histogram2d.h &gsl_histogram2d_pdf_free" histogram2d_pdf_free :: FunPtr (PDFHandle -> IO ()) - ------------------------------------------------------------------------------ - --- | create a histogram PDF from a histogram -fromHistogram2D :: Histogram2D -> Histogram2DPDF -fromHistogram2D (H bx by h) = unsafePerformIO $ do - p <- histogram2d_pdf_new (fromIntegral bx) (fromIntegral by) - p' <- newForeignPtr histogram2d_pdf_free p - withForeignPtr p' $ \p'' -> - withForeignPtr h $ \h' -> do - check "pdf_init" $ histogram2d_pdf_init p'' h' - return $ P p' - --- | given a randomm from the uniform distribution [0,1], draw a random sample from the PDF -sample :: Histogram2DPDF -> Double -sample (P p) = unsafePerformIO $ withForeignPtr p $ \p' -> histogram2d_pdf_sample p' - -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_pdf_init" histogram2d_pdf_init :: PDFHandle -> Hist2DHandle -> IO CInt -foreign import ccall "gsl-histogram2d.h gsl_histogram2d_pdf_sample" histogram2d_pdf_sample :: PDFHandle -> IO Double - ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ diff --git a/lib/Numeric/GSL/Permutation.hs b/lib/Numeric/GSL/Permutation.hs deleted file mode 100644 --- a/lib/Numeric/GSL/Permutation.hs +++ /dev/null @@ -1,349 +0,0 @@ -{-# OPTIONS_GHC -fglasgow-exts #-} ------------------------------------------------------------------------------ --- | --- Module : Numeric.GSL.Permutation --- Copyright : (c) Alexander Vivian Hugh McPhail 2010 --- License : GPL-style --- --- Maintainer : haskell.vivian.mcphail gmail com --- Stability : provisional --- Portability : uses ffi --- --- GSL permutation functions --- ------------------------------------------------------------------------------ - -module Numeric.GSL.Permutation ( - Permutation, CanPerm - , random_permute - , get, swap, swapList - , size - , valid - , reverse - , inverse - , next, prev - , permute, inverse_permute, mul - , fwritePermutation, freadPermutation, fprintfPermutation, fscanfPermutation - , canonical, linear - , inversions - , cyclesLinear, cyclesCanonical - ) where - ------------------------------------------------------------------------------ - -import Data.Packed.Vector ---import Data.Packed.Matrix hiding(toLists) -import Data.Packed.Development - ---import Numeric.LinearAlgebra.Linear - ---import Control.Monad(when) - -import Foreign hiding(shift) ---import Foreign.ForeignPtr ---import Foreign.Marshal.Alloc(alloca) -import Foreign.C.Types(CInt,CChar) ---import Foreign.C.String(newCString,peekCString) -import Foreign.C.String(newCString) - ---import GHC.ForeignPtr (mallocPlainForeignPtrBytes) - ---import GHC.Base ---import GHC.IOBase - -import Prelude hiding(reverse) - ------------------------------------------------------------------------------ - -data Perm -type PermHandle = Ptr Perm --- | A permutation structure -data Permutation = P { pdim :: {-# UNPACK #-} !Int -- ^ number of bins - , perm :: {-# UNPACK #-} !(ForeignPtr Perm) } --- | A canonical permutation structure -data CanPerm = CP { cdim :: {-# UNPACK #-} !Int -- ^ number of bins - , canperm :: {-# UNPACK #-} !(ForeignPtr Perm) } - ------------------------------------------------------------------------------ -{- -instance Eq Permutation where - (==) = equalBins - -instance Num Permutation where - (+) = add - (-) = subtract - negate = flip scale (-1.0) - (*) = multiply - signum = error "can't signumm Permutation" - abs = error "can't abs Permutation" - fromInteger x = fromLimits (fromInteger x) (0,1) - -instance Fractional Permutation where - fromRational x = fromLimits (round x) (0,fromRational x) - (/) = divide --} ------------------------------------------------------------------------------ - -foreign import ccall "gsl-permutation.h gsl_permutation_alloc" permutation_new :: CInt -> IO PermHandle -foreign import ccall "gsl-permutation.h gsl_permutation_calloc" permutation_init :: CInt -> IO PermHandle -foreign import ccall "gsl-permutation.h &gsl_permutation_free" permutation_free :: FunPtr (PermHandle -> IO ()) - ------------------------------------------------------------------------------ - -nullPermutation :: Int -> IO Permutation -nullPermutation n = do - p <- permutation_init (fromIntegral n) - p' <- newForeignPtr permutation_free p - return (P n p') - -clonePermutation :: Permutation -> IO Permutation -clonePermutation (P n s) = do - d <- permutation_new (fromIntegral n) - d' <- newForeignPtr permutation_free d - check "clonePermutation" $ - withForeignPtr s $ \s' -> - withForeignPtr d' $ \d'' -> - permutation_clone d'' s' - return (P n d') - -foreign import ccall "gsl-permutation.h gsl_permutation_memcpy" permutation_clone :: PermHandle -> PermHandle -> IO CInt - ------------------------------------------------------------------------------ - --- | generate a random permutation -random_permute :: Int -- ^ seed - -> Int -- ^ size - -> Permutation -random_permute s n = unsafePerformIO $ do - (P _ p) <- nullPermutation n - check "random_permute" $ - withForeignPtr p $ \p' -> permutation_random_permute (fromIntegral s) p' - return (P n p) - -foreign import ccall "permutation-aux.h random_permute" permutation_random_permute :: CInt -> PermHandle -> IO CInt - ------------------------------------------------------------------------------ - --- | returns the value of the i-th element of the permutation -get :: Permutation -> Int -> Int -get (P _ p) i = unsafePerformIO $ do - j <- withForeignPtr p $ \p' -> permutation_get p' (fromIntegral i) - return $ fromIntegral j - --- | swaps the i-th and j-th elements -swapIO :: Permutation -> Int -> Int -> IO () -swapIO (P _ p) i j = do - check "swap" $ - withForeignPtr p $ \p' -> permutation_swap p' (fromIntegral i) (fromIntegral j) - --- | swaps the i-th and j-th elements -swap :: Permutation -> Int -> Int -> Permutation -swap p i j = unsafePerformIO $ do - p' <- clonePermutation p - swapIO p' i j - return p' - --- | swaps pairs of elements -swapList :: Permutation -> [(Int,Int)] -> Permutation -swapList p xs = unsafePerformIO $ do - (P n p') <- clonePermutation p - withForeignPtr p' $ \p'' -> mapM_ (\(i,j) -> permutation_swap p'' (fromIntegral i) (fromIntegral j)) xs - return (P n p') - -foreign import ccall "gsl-permutation.h gsl_permutation_get" permutation_get :: PermHandle -> CInt -> IO CInt -foreign import ccall "gsl-permutation.h gsl_permutation_swap" permutation_swap :: PermHandle -> CInt -> CInt -> IO CInt - ------------------------------------------------------------------------------ - --- | get the length of the permutation -size :: Permutation -> Int -size (P s _) = s - --- | checks that the permutation is valid -valid :: Permutation -> Bool -valid (P _ p) = unsafePerformIO $ do - v <- withForeignPtr p $ \p' -> permutation_valid p' - if v == 0 - then return False - else return True - -foreign import ccall "gsl-permutation.h gsl_permutation_valid" permutation_valid :: PermHandle -> IO CInt - ------------------------------------------------------------------------------ - --- | reverse the elements of the permutation -reverseIO :: Permutation -> IO () -reverseIO (P _ p) = do - check "reverseIO" $ - withForeignPtr p $ \p' -> permutation_reverse p' - --- | reverse the elements of the permutation -reverse :: Permutation -> Permutation -reverse p = unsafePerformIO $ do - p' <- clonePermutation p - reverseIO p' - return p' - --- | computes the inverse of the permutation -inverse :: Permutation -> Permutation -inverse (P n p) = unsafePerformIO $ do - d <- permutation_new (fromIntegral n) - d' <- newForeignPtr permutation_free d - check "inverse" $ - withForeignPtr d' $ \d'' -> - withForeignPtr p $ \p' -> permutation_inverse d'' p' - return (P n d') - --- | advances the permutation to the next in lexicographic order, if there is one -next :: Permutation -> IO Bool -next (P _ p) = do - err <- withForeignPtr p $ \p' -> permutation_next p' - if err == 0 - then return True - else return False - --- | steps the permutation back to the previous in lexicographic order, if there is one -prev :: Permutation -> IO Bool -prev (P _ p) = do - err <- withForeignPtr p $ \p' -> permutation_prev p' - if err == 0 - then return True - else return False - - -foreign import ccall "gsl-permutation.h gsl_permutation_reverse" permutation_reverse :: PermHandle -> IO CInt -foreign import ccall "gsl-permutation.h gsl_permutation_inverse" permutation_inverse :: PermHandle -> PermHandle -> IO CInt -foreign import ccall "gsl-permutation.h gsl_permutation_next" permutation_next :: PermHandle -> IO CInt -foreign import ccall "gsl-permutation.h gsl_permutation_prev" permutation_prev :: PermHandle -> IO CInt - ------------------------------------------------------------------------------ - --- | apply the permutation to a vector -permute :: Permutation -> Vector Double -> Vector Double -permute (P n p) v = unsafePerformIO $ do - r <- createVector n - app2 (\vs vp rs rp -> withForeignPtr p $ \p' -> permutation_permute p' vs vp rs rp) vec v vec r "permute" - return r - --- | apply the inverse permutation to a vector -inverse_permute :: Permutation -> Vector Double -> Vector Double -inverse_permute (P n p) v = unsafePerformIO $ do - r <- createVector n - app2 (\vs vp rs rp -> withForeignPtr p $ \p' -> permutation_permute_inverse p' vs vp rs rp) vec v vec r "permute" - return r - --- | multiply two permutations, P = PA * PB -mul :: Permutation -> Permutation -> Permutation -mul (P n p1) (P _ p2) = unsafePerformIO $ do - p <- permutation_new (fromIntegral n) - p' <- newForeignPtr permutation_free p - check "mul" $ - withForeignPtr p' $ \p'' -> - withForeignPtr p1 $ \p1' -> - withForeignPtr p2 $ \p2' -> permutation_mul p'' p1' p2' - return (P n p') - -foreign import ccall "permutation-aux.h permute" permutation_permute :: PermHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt -foreign import ccall "permutation-aux.h permute_inverse" permutation_permute_inverse :: PermHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt -foreign import ccall "gsl-permutation.h gsl_permutation_mul" permutation_mul :: PermHandle -> PermHandle -> PermHandle -> IO CInt - ------------------------------------------------------------------------------ - --- | write a permutation in the native binary format (may not be portable) -fwritePermutation :: FilePath -> Permutation -> IO () -fwritePermutation fn (P _ p) = do - cn <- newCString fn - check "fwritePermutation" $ - withForeignPtr p $ permutation_fwrite cn - free cn - --- | read a permutation in the native binary format, length must be known -freadPermutation :: FilePath -> Int -> IO Permutation -freadPermutation fn b = do - h <- permutation_new (fromIntegral b) - h' <- newForeignPtr permutation_free h - cn <- newCString fn - check "freadPermutation" $ - withForeignPtr h' $ permutation_fread cn - return $ P b h' - --- | saves the permutation with the given format -fprintfPermutation :: FilePath -> String -> Permutation -> IO () -fprintfPermutation fn fr (P _ p) = do - cn <- newCString fn - cr <- newCString fr - check "fprintfPermutation" $ - withForeignPtr p $ permutation_fprintf cn cr - free cn - free cr - return () - --- | reads formatted data as written by fprintf, the number of bins must be known in advance -fscanfPermutation :: FilePath -> Int -> IO Permutation -fscanfPermutation fn b = do - h <- permutation_new (fromIntegral b) - h' <- newForeignPtr permutation_free h - cn <- newCString fn - check "fscanfPermutation" $ - withForeignPtr h' $ permutation_fscanf cn - return $ P b h' - -foreign import ccall "permutation-aux.h perm_fwrite" permutation_fwrite :: Ptr CChar -> PermHandle -> IO CInt -foreign import ccall "permutation-aux.h perm_fread" permutation_fread :: Ptr CChar -> PermHandle -> IO CInt -foreign import ccall "permutation-aux.h perm_fprintf" permutation_fprintf :: Ptr CChar -> Ptr CChar -> PermHandle -> IO CInt -foreign import ccall "permutation-aux.h perm_fscanf" permutation_fscanf :: Ptr CChar -> PermHandle -> IO CInt - ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ - --- | compute the canonical form -canonical :: Permutation -> CanPerm -canonical (P n p) = unsafePerformIO $ do - q <- permutation_new (fromIntegral n) - q' <- newForeignPtr permutation_free q - check "canonical" $ - withForeignPtr p $ \p' -> - withForeignPtr q' $ \q'' -> - permutation_linear_to_canonical q'' p' - return (CP n q') - --- | convert from canonical to linear -linear :: CanPerm -> Permutation -linear (CP n p) = unsafePerformIO $ do - q <- permutation_new (fromIntegral n) - q' <- newForeignPtr permutation_free q - check "linear" $ - withForeignPtr p $ \p' -> - withForeignPtr q' $ \q'' -> - permutation_canonical_to_linear q'' p' - return (P n q') - --- | a count of the inversions -inversions :: Permutation -> Int -inversions (P _ p) = unsafePerformIO $ do - i <- withForeignPtr p $ \p' -> permutation_inversions p' - return $ fromIntegral i - --- | a count of the cycles of a permutation in linear form -cyclesLinear :: Permutation -> Int -cyclesLinear (P _ p) = unsafePerformIO $ do - i <- withForeignPtr p $ \p' -> permutation_linear_cycles p' - return $ fromIntegral i - --- | a count of the cycles of a permutation in canonical form -cyclesCanonical :: CanPerm -> Int -cyclesCanonical (CP _ p) = unsafePerformIO $ do - i <- withForeignPtr p $ \p' -> permutation_canonical_cycles p' - return $ fromIntegral i - - -foreign import ccall "gsl-permutation.h gsl_permutation_linear_to_canonical" permutation_linear_to_canonical :: PermHandle -> PermHandle -> IO CInt -foreign import ccall "gsl-permutation.h gsl_permutation_canonical_to_linear" permutation_canonical_to_linear :: PermHandle -> PermHandle -> IO CInt -foreign import ccall "gsl-permutation.h gsl_permutation_inversions" permutation_inversions :: PermHandle -> IO CInt -foreign import ccall "gsl-permutation.h gsl_permutation_linear_cycles" permutation_linear_cycles :: PermHandle -> IO CInt -foreign import ccall "gsl-permutation.h gsl_permutation_canonical_cycles" permutation_canonical_cycles :: PermHandle -> IO CInt - ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ diff --git a/lib/Numeric/GSL/Sort.hs b/lib/Numeric/GSL/Sort.hs deleted file mode 100644 --- a/lib/Numeric/GSL/Sort.hs +++ /dev/null @@ -1,34 +0,0 @@ -{-# OPTIONS_GHC -fglasgow-exts #-} ------------------------------------------------------------------------------ --- | --- Module : Numeric.GSL.Sort --- Copyright : (c) Alexander Vivian Hugh McPhail 2010 --- License : GPL-style --- --- Maintainer : haskell.vivian.mcphail gmail com --- Stability : provisional --- Portability : uses ffi --- --- GSL sorting functions --- ------------------------------------------------------------------------------ - -module Numeric.GSL.Sort ( - sort - ) where - - -import Data.Packed.Vector -import Data.Packed.Development - -import Foreign -import Foreign.C.Types(CInt) - --- | sort the elements of a vector into ascending order -sort :: Vector Double -> Vector Double -sort v = unsafePerformIO $ do - r <- createVector (dim v) - app2 sort_sort vec v vec r "sort" - return r - -foreign import ccall "sort-aux.h sort" sort_sort :: CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt diff --git a/lib/Numeric/GSL/Statistics.hs b/lib/Numeric/GSL/Statistics.hs deleted file mode 100644 --- a/lib/Numeric/GSL/Statistics.hs +++ /dev/null @@ -1,308 +0,0 @@ -{-# OPTIONS_GHC -fglasgow-exts #-} ------------------------------------------------------------------------------ --- | --- Module : Numeric.GSL.Statistics --- Copyright : (c) Alexander Vivian Hugh McPhail 2010 --- License : GPL-style --- --- Maintainer : haskell.vivian.mcphail gmail com --- Stability : provisional --- Portability : uses ffi --- --- GSL statistics functions --- ------------------------------------------------------------------------------ - -module Numeric.GSL.Statistics ( - mean - , variance,variance_m,variance_pm - , stddev,stddev_m,stddev_pm - , tot_sumsq,tot_sumsq_m - , absdev, absdev_m - , skew, skew_m_sd - , kurtosis, kurtosis_m_sd - -- - , mean_w - , variance_w,variance_w_m,variance_w_pm - , stddev_w,stddev_w_m,stddev_w_pm - , tot_sumsq_w,tot_sumsq_w_m - , absdev_w, absdev_w_m - , skew_w, skew_w_m_sd - , kurtosis_w, kurtosis_w_m_sd - -- - , lag1auto - , covariance, covariance_m - , correlation - -- - , median, quantile - ) where - ------------------------------------------------------------------------------ - -import Data.Packed.Vector ---import Data.Packed(Container(..)) - -import Data.Packed.Development - ---import Numeric.GSL.Vector ---import Numeric.LinearAlgebra.Instances() ---import Numeric.LinearAlgebra.Linear(Linear(..)) - -import Foreign -import Foreign.C.Types(CInt) ---import Foreign.Marshal.Alloc(alloca) - ------------------------------------------------------------------------------ - -type PD = Ptr Double - ------------------------------------------------------------------------------ - -getD1 f s v = unsafePerformIO $ do - alloca $ \r -> do - app1 (f r) vec v s - r' <- peek r - return r' - -getD2 f s v w = unsafePerformIO $ do - alloca $ \r -> do - app2 (f r) vec v vec w s - r' <- peek r - return r' - ------------------------------------------------------------------------------ - --- | the mean of the elements of a vector -mean :: Vector Double -> Double -mean = getD1 statistics_mean "mean" - --- | the sample variance -variance :: Vector Double -> Double -variance = getD1 statistics_variance "variance" - --- | the sample variance given the precomputed mean -variance_m :: Double -> Vector Double -> Double -variance_m m = getD1 (statistics_variance_m m) "variance_m" - --- | the population variance given the a priori mean -variance_pm :: Double -> Vector Double -> Double -variance_pm m = getD1 (statistics_var_with_fixed_m m) "variance_pm" - --- | the sample standard deviation -stddev :: Vector Double -> Double -stddev = getD1 statistics_stddev "stddev" - --- | the sample standard deviation given the precomputed mean -stddev_m :: Double -> Vector Double -> Double -stddev_m m = getD1 (statistics_stddev_m m) "stddev_m" - --- | the population standard deviation given the a priori mean -stddev_pm :: Double -> Vector Double -> Double -stddev_pm m = getD1 (statistics_stddev_with_fixed_m m) "stddev_pm" - --- | the total sum of squares about the mean -tot_sumsq :: Vector Double -> Double -tot_sumsq = getD1 statistics_tot_sumsq "tot_sumsq" - --- | the total sum of squares about the precomputed mean -tot_sumsq_m :: Double -> Vector Double -> Double -tot_sumsq_m m = getD1 (statistics_tot_sumsq_m m) "totsumsq_m" - ------------------------------------------------------------------------------ - -foreign import ccall "statistics-aux.h mean" statistics_mean :: PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h variance" statistics_variance :: PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h variance_m" statistics_variance_m :: Double -> PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h var_with_fixed_m" statistics_var_with_fixed_m :: Double -> PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h stddev" statistics_stddev :: PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h stddev_m" statistics_stddev_m :: Double -> PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h stddev_with_fixed_m" statistics_stddev_with_fixed_m :: Double -> PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h tot_sumsq" statistics_tot_sumsq :: PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h tot_sumsq_m" statistics_tot_sumsq_m :: Double -> PD -> CInt -> PD -> IO CInt - ------------------------------------------------------------------------------ - --- | the absolute deviation from the mean -absdev :: Vector Double -> Double -absdev = getD1 statistics_absdev "absdev" - --- | the absolute deviation from the precomputed mean -absdev_m :: Double -> Vector Double -> Double -absdev_m m = getD1 (statistics_absdev_m m) "absdev_m" - ------------------------------------------------------------------------------ - -foreign import ccall "statistics-aux.h absdev" statistics_absdev :: PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h absdev_m" statistics_absdev_m :: Double -> PD -> CInt -> PD -> IO CInt - ------------------------------------------------------------------------------ - --- | the skewness of the data (asymmetry of tails) -skew :: Vector Double -> Double -skew = getD1 statistics_skew "skew" - --- | the skewness of the data (asymmetry of tails) with precomputed mean and sd -skew_m_sd :: Double -> Double -> Vector Double -> Double -skew_m_sd m sd = getD1 (statistics_skew_m_sd m sd) "skew_m_sd" - --- | the kurtosis of the data (sharpness of peak relative to width) -kurtosis :: Vector Double -> Double -kurtosis = getD1 statistics_kurtosis "kurtosis" - --- | the kurtosis of the data (sharpness of peak relative to width) with precomputed mean and sd -kurtosis_m_sd :: Double -> Double -> Vector Double -> Double -kurtosis_m_sd m sd = getD1 (statistics_kurtosis_m_sd m sd) "kurtosis_m_sd" - ------------------------------------------------------------------------------ - -foreign import ccall "statistics-aux.h skew" statistics_skew :: PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h skew_m_sd" statistics_skew_m_sd :: Double -> Double -> PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h kurtosis" statistics_kurtosis :: PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h kurtosis_m_sd" statistics_kurtosis_m_sd :: Double -> Double -> PD -> CInt -> PD -> IO CInt - ------------------------------------------------------------------------------ - --- | the lag-1 autocorrelation of the data -lag1auto :: Vector Double -> Double -lag1auto = getD1 statistics_lag1auto "lag1auto" - ------------------------------------------------------------------------------ - -foreign import ccall "statistics-aux.h lag1_autocorrelation" statistics_lag1auto :: PD -> CInt -> PD -> IO CInt - ----------------------------------------------------------------------------- - --- | the covariance of two datasets of the same length -covariance :: Vector Double -> Vector Double -> Double -covariance = getD2 statistics_covariance "covariance" - --- | the covariance of two datasets of the same length -covariance_m :: Double -> Double -> Vector Double -> Vector Double -> Double -covariance_m m1 m2 = getD2 (statistics_covariance_m m1 m2) "covariance_m" - ------------------------------------------------------------------------------ - -foreign import ccall "statistics-aux.h covariance" statistics_covariance :: PD -> CInt -> PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h covariance_m" statistics_covariance_m :: Double -> Double -> PD -> CInt -> PD -> CInt -> PD -> IO CInt - ------------------------------------------------------------------------------ - --- | the Pearson correlation of two datasets of the same length -correlation :: Vector Double -> Vector Double -> Double -correlation = getD2 statistics_correlation "correlation" - ------------------------------------------------------------------------------ - -foreign import ccall "statistics-aux.h correlation" statistics_correlation :: PD -> CInt -> PD -> CInt -> PD -> IO CInt - ------------------------------------------------------------------------------ - --- | the weighted mean of the elements of a vector -mean_w :: Vector Double -- ^ weights - -> Vector Double -- ^ dataset - -> Double -mean_w = getD2 statistics_w_mean "w_mean" - --- | the weighted sample variance -variance_w :: Vector Double -> Vector Double -> Double -variance_w = getD2 statistics_w_variance "w_variance" - --- | the weighted sample variance given the precomputed mean -variance_w_m :: Double -> Vector Double -> Vector Double -> Double -variance_w_m m = getD2 (statistics_w_variance_m m) "w_variance_m" - --- | the weighted population variance given the a priori mean -variance_w_pm :: Double -> Vector Double -> Vector Double -> Double -variance_w_pm m = getD2 (statistics_w_var_with_fixed_m m) "w_variance_pm" - --- | the weighted sample standard deviation -stddev_w :: Vector Double -> Vector Double -> Double -stddev_w = getD2 statistics_w_stddev "w_stddev" - --- | the weighted sample standard deviation given the precomputed mean -stddev_w_m :: Double -> Vector Double -> Vector Double -> Double -stddev_w_m m = getD2 (statistics_w_stddev_m m) "w_stddev_m" - --- | the weighted population standard deviation given the a priori mean -stddev_w_pm :: Double -> Vector Double -> Vector Double -> Double -stddev_w_pm m = getD2 (statistics_w_stddev_with_fixed_m m) "w_stddev_pm" - --- | the weighted total sum of squares about the mean -tot_sumsq_w :: Vector Double -> Vector Double -> Double -tot_sumsq_w = getD2 statistics_w_tot_sumsq "w_tot_sumsq" - --- | the weighted total sum of squares about the precomputed mean -tot_sumsq_w_m :: Double -> Vector Double -> Vector Double -> Double -tot_sumsq_w_m m = getD2 (statistics_w_tot_sumsq_m m) "w_totsumsq_m" - ------------------------------------------------------------------------------ - -foreign import ccall "statistics-aux.h w_mean" statistics_w_mean :: PD -> CInt -> PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h w_variance" statistics_w_variance :: PD -> CInt -> PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h w_variance_m" statistics_w_variance_m :: Double -> PD -> CInt -> PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h w_var_with_fixed_m" statistics_w_var_with_fixed_m :: Double -> PD -> CInt -> PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h w_stddev" statistics_w_stddev :: PD -> CInt -> PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h w_stddev_m" statistics_w_stddev_m :: Double -> PD -> CInt -> PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h w_stddev_with_fixed_m" statistics_w_stddev_with_fixed_m :: Double -> PD -> CInt -> PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h w_tot_sumsq" statistics_w_tot_sumsq :: PD -> CInt -> PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h w_tot_sumsq_m" statistics_w_tot_sumsq_m :: Double -> PD -> CInt -> PD -> CInt -> PD -> IO CInt - ------------------------------------------------------------------------------ - --- | the weighted absolute deviation from the mean -absdev_w :: Vector Double -> Vector Double -> Double -absdev_w = getD2 statistics_w_absdev "w_absdev" - --- | the weighted absolute deviation from the precomputed mean -absdev_w_m :: Double -> Vector Double -> Vector Double -> Double -absdev_w_m m = getD2 (statistics_w_absdev_m m) "w_absdev_m" - ------------------------------------------------------------------------------ - -foreign import ccall "statistics-aux.h w_absdev" statistics_w_absdev :: PD -> CInt -> PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h w_absdev_m" statistics_w_absdev_m :: Double -> PD -> CInt -> PD -> CInt -> PD -> IO CInt - ------------------------------------------------------------------------------ - --- | the weighted skewness of the data (asymmetry of tails) -skew_w :: Vector Double -> Vector Double -> Double -skew_w = getD2 statistics_w_skew "w_skew" - --- | the weighted skewness of the data (asymmetry of tails) with precomputed mean and sd -skew_w_m_sd :: Double -> Double -> Vector Double -> Vector Double -> Double -skew_w_m_sd m sd = getD2 (statistics_w_skew_m_sd m sd) "w_skew_m_sd" - --- | the weighted kurtosis of the data (sharpness of peak relative to width) -kurtosis_w :: Vector Double -> Vector Double -> Double -kurtosis_w = getD2 statistics_w_kurtosis "w_kurtosis" - --- | the weighted kurtosis of the data (sharpness of peak relative to width) with precomputed mean and sd -kurtosis_w_m_sd :: Double -> Double -> Vector Double -> Vector Double -> Double -kurtosis_w_m_sd m sd = getD2 (statistics_w_kurtosis_m_sd m sd) "w_kurtosis_m_sd" - ------------------------------------------------------------------------------ - -foreign import ccall "statistics-aux.h w_skew" statistics_w_skew :: PD -> CInt -> PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h w_skew_m_sd" statistics_w_skew_m_sd :: Double -> Double -> PD -> CInt -> PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h w_kurtosis" statistics_w_kurtosis :: PD -> CInt -> PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h w_kurtosis_m_sd" statistics_w_kurtosis_m_sd :: Double -> Double -> PD -> CInt -> PD -> CInt -> PD -> IO CInt - ------------------------------------------------------------------------------ - --- | the median value of the dataset, which must be sorted -median :: Vector Double -> Double -median = getD1 statistics_median "median" - --- | the quantile value of the dataset, which must be sorted -quantile :: Double -- ^ the desired quantile from [0..1] - -> Vector Double -- ^ the dataset - -> Double -quantile f = getD1 (statistics_quantile f) "quantile" - ------------------------------------------------------------------------------ - -foreign import ccall "statistics-aux.h median" statistics_median :: PD -> CInt -> PD -> IO CInt -foreign import ccall "statistics-aux.h quantile" statistics_quantile :: Double -> PD -> CInt -> PD -> IO CInt - ------------------------------------------------------------------------------ - diff --git a/lib/Numeric/GSL/distribution-aux.c b/lib/Numeric/GSL/distribution-aux.c deleted file mode 100644 --- a/lib/Numeric/GSL/distribution-aux.c +++ /dev/null @@ -1,749 +0,0 @@ -#include -#include -#include - -int random1(int s, int type, double par, double* r) -{ - const gsl_rng_type * T; - gsl_rng * rng; - - gsl_rng_env_setup(); - T = gsl_rng_default; - rng = gsl_rng_alloc(T); - gsl_rng_set(rng,s); - - switch(type) { - case 0: { (*r) = gsl_ran_gaussian(rng,par); break; } - case 1: { (*r) = gsl_ran_exponential(rng,par); break; } - case 2: { (*r) = gsl_ran_laplace(rng,par); break; } - case 3: { (*r) = gsl_ran_cauchy(rng,par); break; } - case 4: { (*r) = gsl_ran_rayleigh(rng,par); break; } - case 5: { (*r) = gsl_ran_chisq(rng,par); break; } - case 6: { (*r) = gsl_ran_tdist(rng,par); break; } - case 7: { (*r) = gsl_ran_logistic(rng,par); break; } - } - - gsl_rng_free(rng); - - return 0; -} - -double random1_pdf(int type, double x, double par) -{ - switch (type) { - case 0: return gsl_ran_gaussian_pdf(x,par); - case 1: return gsl_ran_exponential_pdf(x,par); - case 2: return gsl_ran_laplace_pdf(x,par); - case 3: return gsl_ran_cauchy_pdf(x,par); - case 4: return gsl_ran_rayleigh_pdf(x,par); - case 5: return gsl_ran_chisq_pdf(x,par); - case 6: return gsl_ran_tdist_pdf(x,par); - case 7: return gsl_ran_logistic_pdf(x,par); - } -} - -double random1_cdf_lower(int type, double x, double par) -{ - switch (type) { - case 0: return gsl_cdf_gaussian_P(x,par); - case 1: return gsl_cdf_exponential_P(x,par); - case 2: return gsl_cdf_laplace_P(x,par); - case 3: return gsl_cdf_cauchy_P(x,par); - case 4: return gsl_cdf_rayleigh_P(x,par); - case 5: return gsl_cdf_chisq_P(x,par); - case 6: return gsl_cdf_tdist_P(x,par); - case 7: return gsl_cdf_logistic_P(x,par); - } -} - -double random1_cdf_upper(int type, double x, double par) -{ - switch (type) { - case 0: return gsl_cdf_gaussian_Q(x,par); - case 1: return gsl_cdf_exponential_Q(x,par); - case 2: return gsl_cdf_laplace_Q(x,par); - case 3: return gsl_cdf_cauchy_Q(x,par); - case 4: return gsl_cdf_rayleigh_Q(x,par); - case 5: return gsl_cdf_chisq_Q(x,par); - case 6: return gsl_cdf_tdist_Q(x,par); - case 7: return gsl_cdf_logistic_Q(x,par); - } -} - -double random1_cdf_lower_inv(int type, double X, double par) -{ - switch (type) { - case 0: return gsl_cdf_gaussian_Pinv(X,par); - case 1: return gsl_cdf_exponential_Pinv(X,par); - case 2: return gsl_cdf_laplace_Pinv(X,par); - case 3: return gsl_cdf_cauchy_Pinv(X,par); - case 4: return gsl_cdf_rayleigh_Pinv(X,par); - case 5: return gsl_cdf_chisq_Pinv(X,par); - case 6: return gsl_cdf_tdist_Pinv(X,par); - case 7: return gsl_cdf_logistic_Pinv(X,par); - } -} - -double random1_cdf_upper_inv(int type, double X, double par) -{ - switch (type) { - case 0: return gsl_cdf_gaussian_Qinv(X,par); - case 1: return gsl_cdf_exponential_Qinv(X,par); - case 2: return gsl_cdf_laplace_Qinv(X,par); - case 3: return gsl_cdf_cauchy_Qinv(X,par); - case 4: return gsl_cdf_rayleigh_Qinv(X,par); - case 5: return gsl_cdf_chisq_Qinv(X,par); - case 6: return gsl_cdf_tdist_Qinv(X,par); - case 7: return gsl_cdf_logistic_Qinv(X,par); - } -} - -double random1_dist(int df, int type, double x, double par) -{ - switch(df) { - case 0: return random1_pdf(type,x,par); - case 1: return random1_cdf_lower(type,x,par); - case 2: return random1_cdf_upper(type,x,par); - case 3: return random1_cdf_lower_inv(type,x,par); - case 4: return random1_cdf_upper_inv(type,x,par); - } -} - -////////////////////////////////////////////////////////////////////// - -int random2(int s, int type, double par1, double par2, double* r) -{ - const gsl_rng_type * T; - gsl_rng * rng; - - gsl_rng_env_setup(); - T = gsl_rng_default; - rng = gsl_rng_alloc(T); - gsl_rng_set(rng,s); - - switch(type) { - case 0: { (*r) = gsl_ran_gaussian_tail(rng,par1,par2); break; } - case 1: { (*r) = gsl_ran_exppow(rng,par1,par2); break; } - case 2: { (*r) = gsl_ran_rayleigh_tail(rng,par1,par2); break; } - case 3: { (*r) = gsl_ran_levy(rng,par1,par2); break; } - case 4: { (*r) = gsl_ran_gamma(rng,par1,par2); break; } - case 5: { (*r) = gsl_ran_flat(rng,par1,par2); break; } - case 6: { (*r) = gsl_ran_lognormal(rng,par1,par2); break; } - case 7: { (*r) = gsl_ran_fdist(rng,par1,par2); break; } - case 8: { (*r) = gsl_ran_beta(rng,par1,par2); break; } - case 9: { (*r) = gsl_ran_pareto(rng,par1,par2); break; } - case 10: { (*r) = gsl_ran_weibull(rng,par1,par2); break; } - case 11: { (*r) = gsl_ran_gumbel1(rng,par1,par2); break; } - case 12: { (*r) = gsl_ran_gumbel2(rng,par1,par2); break; } - } - - gsl_rng_free(rng); - - return 0; -} - -double random2_pdf(int type, double x, double par1, double par2) -{ - switch (type) { - case 0: return gsl_ran_gaussian_tail_pdf(x,par1,par2); - case 1: return gsl_ran_exppow_pdf(x,par1,par2); - case 2: return gsl_ran_rayleigh_tail_pdf(x,par1,par2); - case 4: return gsl_ran_gamma_pdf(x,par1,par2); - case 5: return gsl_ran_flat_pdf(x,par1,par2); - case 6: return gsl_ran_lognormal_pdf(x,par1,par2); - case 7: return gsl_ran_fdist_pdf(x,par1,par2); - case 8: return gsl_ran_beta_pdf(x,par1,par2); - case 9: return gsl_ran_pareto_pdf(x,par1,par2); - case 10: return gsl_ran_weibull_pdf(x,par1,par2); - case 11: return gsl_ran_gumbel1_pdf(x,par1,par2); - case 12: return gsl_ran_gumbel2_pdf(x,par1,par2); - } -} - -double random2_cdf_lower(int type, double x, double par1, double par2) -{ - switch (type) { - case 1: return gsl_cdf_exppow_P(x,par1,par2); - case 4: return gsl_cdf_gamma_P(x,par1,par2); - case 5: return gsl_cdf_flat_P(x,par1,par2); - case 6: return gsl_cdf_lognormal_P(x,par1,par2); - case 7: return gsl_cdf_fdist_P(x,par1,par2); - case 8: return gsl_cdf_beta_P(x,par1,par2); - case 9: return gsl_cdf_pareto_P(x,par1,par2); - case 10: return gsl_cdf_weibull_P(x,par1,par2); - case 11: return gsl_cdf_gumbel1_P(x,par1,par2); - case 12: return gsl_cdf_gumbel2_P(x,par1,par2); - } -} - -double random2_cdf_upper(int type, double x, double par1, double par2) -{ - switch (type) { - case 1: return gsl_cdf_exppow_Q(x,par1,par2); - case 4: return gsl_cdf_gamma_Q(x,par1,par2); - case 5: return gsl_cdf_flat_Q(x,par1,par2); - case 6: return gsl_cdf_lognormal_Q(x,par1,par2); - case 7: return gsl_cdf_fdist_Q(x,par1,par2); - case 8: return gsl_cdf_beta_Q(x,par1,par2); - case 9: return gsl_cdf_pareto_Q(x,par1,par2); - case 10: return gsl_cdf_weibull_Q(x,par1,par2); - case 11: return gsl_cdf_gumbel1_Q(x,par1,par2); - case 12: return gsl_cdf_gumbel2_Q(x,par1,par2); - } -} - -double random2_cdf_lower_inv(int type, double X, double par1, double par2) -{ - switch (type) { - case 4: return gsl_cdf_gamma_Pinv(X,par1,par2); - case 5: return gsl_cdf_flat_Pinv(X,par1,par2); - case 6: return gsl_cdf_lognormal_Pinv(X,par1,par2); - case 7: return gsl_cdf_fdist_Pinv(X,par1,par2); - case 8: return gsl_cdf_beta_Pinv(X,par1,par2); - case 9: return gsl_cdf_pareto_Pinv(X,par1,par2); - case 10: return gsl_cdf_weibull_Pinv(X,par1,par2); - case 11: return gsl_cdf_gumbel1_Pinv(X,par1,par2); - case 12: return gsl_cdf_gumbel2_Pinv(X,par1,par2); - } -} - -double random2_cdf_upper_inv(int type, double X, double par1, double par2) -{ - switch (type) { - case 4: return gsl_cdf_gamma_Qinv(X,par1,par2); - case 5: return gsl_cdf_flat_Qinv(X,par1,par2); - case 6: return gsl_cdf_lognormal_Qinv(X,par1,par2); - case 7: return gsl_cdf_fdist_Qinv(X,par1,par2); - case 8: return gsl_cdf_beta_Qinv(X,par1,par2); - case 9: return gsl_cdf_pareto_Qinv(X,par1,par2); - case 10: return gsl_cdf_weibull_Qinv(X,par1,par2); - case 11: return gsl_cdf_gumbel1_Qinv(X,par1,par2); - case 12: return gsl_cdf_gumbel2_Qinv(X,par1,par2); - } -} - -double random2_dist(int df, int type, double x, double par1, double par2) -{ - switch(df) { - case 0: return random2_pdf(type,x,par1,par2); - case 1: return random2_cdf_lower(type,x,par1,par2); - case 2: return random2_cdf_upper(type,x,par1,par2); - case 3: return random2_cdf_lower_inv(type,x,par1,par2); - case 4: return random2_cdf_upper_inv(type,x,par1,par2); - } -} - -////////////////////////////////////////////////////////////////////// - -int random3(int s, int type, double par1, double par2, double par3, double* r) -{ - const gsl_rng_type * T; - gsl_rng * rng; - - gsl_rng_env_setup(); - T = gsl_rng_default; - rng = gsl_rng_alloc(T); - gsl_rng_set(rng,s); - - switch(type) { - case 0: { (*r) = gsl_ran_levy_skew(rng,par1,par2,par3); break; } - } - - gsl_rng_free(rng); - - return 0; -} - -double random3_pdf(int type, double x, double par1, double par2, double par3) -{ - switch (type) { - default : return 0; - } -} - -double random3_cdf_lower(int type, double x, double par1, double par2, double par3) -{ - switch (type) { - default : return 0; - } -} - -double random3_cdf_upper(int type, double x, double par1, double par2, double par3) -{ - switch (type) { - default : return 0; - } -} - -double random3_cdf_lower_inv(int type, double X, double par1, double par2, double par3) -{ - switch (type) { - default : return 0; - } -} - -double random3_cdf_upper_inv(int type, double X, double par1, double par2, double par3) -{ - switch (type) { - default : return 0; - } -} - -double random3_dist(int df, int type, double x, double par1, double par2, double par3) -{ - switch(df) { - case 0: return random3_pdf(type,x,par1,par2,par3); - case 1: return random3_cdf_lower(type,x,par1,par2,par3); - case 2: return random3_cdf_upper(type,x,par1,par2,par3); - case 3: return random3_cdf_lower_inv(type,x,par1,par2,par3); - case 4: return random3_cdf_upper_inv(type,x,par1,par2,par3); - } -} - -////////////////////////////////////////////////////////////////////// - -int random_biv(int s, int type, double par1, double par2, double par3, double* r1, double* r2) -{ - const gsl_rng_type * T; - gsl_rng * rng; - - gsl_rng_env_setup(); - T = gsl_rng_default; - rng = gsl_rng_alloc(T); - gsl_rng_set(rng,s); - - switch(type) { - case 0: { gsl_ran_bivariate_gaussian(rng,par1,par2,par3,r1,r2); break; } - } - - gsl_rng_free(rng); - - return 0; -} - -double random_biv_pdf(int type, double x, double y, double par1, double par2, double par3) -{ - switch (type) { - case 0: return gsl_ran_bivariate_gaussian_pdf(x,y,par1,par2,par3); - } -} - -double random_biv_cdf_lower(int type, double x, double y, double par1, double par2, double par3) -{ - switch (type) { - default : return 0; - } -} - -double random_biv_cdf_upper(int type, double x, double y, double par1, double par2, double par3) -{ - switch (type) { - default : return 0; - } -} - -double random_biv_cdf_lower_inv(int type, double x, double y, double par1, double par2, double par3) -{ - switch (type) { - default : return 0; - } -} - -double random_biv_cdf_upper_inv(int type, double x, double y, double par1, double par2, double par3) -{ - switch (type) { - default : return 0; - } -} - -double random_biv_dist(int df, int type, double x, double y, double par1, double par2, double par3) -{ - switch(df) { - case 0: return random_biv_pdf(type,x,y,par1,par2,par3); - case 1: return random_biv_cdf_lower(type,x,y,par1,par2,par3); - case 2: return random_biv_cdf_upper(type,x,y,par1,par2,par3); - case 3: return random_biv_cdf_lower_inv(type,x,y,par1,par2,par3); - case 4: return random_biv_cdf_upper_inv(type,x,y,par1,par2,par3); - } -} - -////////////////////////////////////////////////////////////////////// - -int random_mp(int s, int type, int ps, const double* p, int rs, double* r) -{ - if (ps != rs) return 2000; // BAD_SIZE - - const gsl_rng_type * T; - gsl_rng * rng; - - gsl_rng_env_setup(); - T = gsl_rng_default; - rng = gsl_rng_alloc(T); - gsl_rng_set(rng,s); - - switch(type) { - case 0: { gsl_ran_dirichlet(rng,ps,p,r); break; } - } - - gsl_rng_free(rng); - - return 0; -} - -double random_mp_pdf(int type, int ps, const double* p, int qs, const double* q) -{ - switch (type) { - case 0: return gsl_ran_dirichlet_pdf(ps,p,q); - } - return 0; -} - -double random_mp_cdf_lower(int type, int ps, const double* p, int qs, const double* q) -{ - switch (type) { - default : return 0; - } -} - -double random_mp_cdf_upper(int type, int ps, const double* p, int qs, const double* q) -{ - switch (type) { - default : return 0; - } -} - -double random_mp_cdf_lower_inv(int type, int ps, const double* p, int qs, const double* q) -{ - switch (type) { - default : return 0; - } -} - -double random_mp_cdf_upper_inv(int type, int ps, const double* p, int qs, const double* q) -{ - switch (type) { - default : return 0; - } -} - - int random_mp_dist(int df, int type, double* r, int ps, const double* p, int qs, const double* q) -{ - switch(df) { - case 0: { (*r) = random_mp_pdf(type,ps,p,qs,q); break; } - case 1: { (*r) = random_mp_cdf_lower(type,ps,p,qs,q); break; } - case 2: { (*r) = random_mp_cdf_upper(type,ps,p,qs,q); break; } - case 3: { (*r) = random_mp_cdf_lower_inv(type,ps,p,qs,q); break; } - case 4: { (*r) = random_mp_cdf_upper_inv(type,ps,p,qs,q); break; } - } - return 0; -} - -////////////////////////////////////////////////////////////////////// - -int random_vector(int s, int rs, double* r) -{ - const gsl_rng_type * T; - gsl_rng * rng; - - gsl_rng_env_setup(); - T = gsl_rng_default; - rng = gsl_rng_alloc(T); - gsl_rng_set(rng,s); - - gsl_ran_dir_nd(rng,rs,r); - - gsl_rng_free(rng); - - return 0; -} - -////////////////////////////////////////////////////////////////////// - -int discrete1(int s, int type, double par, unsigned int* r) -{ - const gsl_rng_type * T; - gsl_rng * rng; - - gsl_rng_env_setup(); - T = gsl_rng_default; - rng = gsl_rng_alloc(T); - gsl_rng_set(rng,s); - - switch(type) { - case 0: { (*r) = gsl_ran_poisson(rng,par); break; } - case 1: { (*r) = gsl_ran_bernoulli(rng,par); break; } - case 2: { (*r) = gsl_ran_laplace(rng,par); break; } - case 3: { (*r) = gsl_ran_logarithmic(rng,par); break; } - } - - gsl_rng_free(rng); - - return 0; -} - -double discrete1_pdf(int type, unsigned int x, double par) -{ - switch (type) { - case 0: return gsl_ran_poisson_pdf(x,par); - case 1: return gsl_ran_bernoulli_pdf(x,par); - case 2: return gsl_ran_geometric_pdf(x,par); - case 3: return gsl_ran_logarithmic_pdf(x,par); - } -} - -double discrete1_cdf_lower(int type, unsigned int x, double par) -{ - switch (type) { - case 0: return gsl_cdf_poisson_P(x,par); - case 2: return gsl_cdf_geometric_P(x,par); - } -} - -double discrete1_cdf_upper(int type, unsigned int x, double par) -{ - switch (type) { - case 0: return gsl_cdf_poisson_Q(x,par); - case 2: return gsl_cdf_geometric_Q(x,par); - } -} - -double discrete1_cdf_lower_inv(int type, unsigned int X, double par) -{ - switch (type) { - default : return 0; - } -} - -double discrete1_cdf_upper_inv(int type, unsigned int X, double par) -{ - switch (type) { - default : return 0; - } -} - -double discrete1_dist(int df, int type, unsigned int x, double par) -{ - switch(df) { - case 0: return discrete1_pdf(type,x,par); - case 1: return discrete1_cdf_lower(type,x,par); - case 2: return discrete1_cdf_upper(type,x,par); - case 3: return discrete1_cdf_lower_inv(type,x,par); - case 4: return discrete1_cdf_upper_inv(type,x,par); - } -} - -////////////////////////////////////////////////////////////////////// - -int discrete2(int s, int type, double par1, unsigned int par2, unsigned int* r) -{ - const gsl_rng_type * T; - gsl_rng * rng; - - gsl_rng_env_setup(); - T = gsl_rng_default; - rng = gsl_rng_alloc(T); - gsl_rng_set(rng,s); - - switch(type) { - case 0: { (*r) = gsl_ran_binomial(rng,par1,par2); break; } - case 1: { (*r) = gsl_ran_negative_binomial(rng,par1,par2*1.0); break; } - case 2: { (*r) = gsl_ran_pascal(rng,par1,par2); break; } - } - - gsl_rng_free(rng); - - return 0; -} - -double discrete2_pdf(int type, unsigned int x, double par1, unsigned int par2) -{ - switch (type) { - case 0: return gsl_ran_binomial_pdf(x,par1,par2); - case 1: return gsl_ran_negative_binomial_pdf(x,par1,par2*1.0); - case 2: return gsl_ran_pascal_pdf(x,par1,par2); - } -} - -double discrete2_cdf_lower(int type, unsigned int x, double par1, unsigned int par2) -{ - switch (type) { - case 0: return gsl_cdf_binomial_P(x,par1,par2); - case 1: return gsl_cdf_negative_binomial_P(x,par1,par2*1.0); - case 2: return gsl_cdf_pascal_P(x,par1,par2*1.0); - } -} - -double discrete2_cdf_upper(int type, unsigned int x, double par1, unsigned int par2) -{ - switch (type) { - case 0: return gsl_cdf_binomial_Q(x,par1,par2); - case 1: return gsl_cdf_negative_binomial_Q(x,par1,par2*1.0); - case 2: return gsl_cdf_pascal_Q(x,par1,par2*1.0); - } -} - -double discrete2_cdf_lower_inv(int type, unsigned int X, double par1, unsigned int par2) -{ - switch (type) { - default : return 0; - } -} - -double discrete2_cdf_upper_inv(int type, unsigned int X, double par1, unsigned int par2) -{ - switch (type) { - default : return 0; - } -} - -double discrete2_dist(int df, int type, unsigned int x, double par1, unsigned int par2) -{ - switch(df) { - case 0: return discrete2_pdf(type,x,par1,par2); - case 1: return discrete2_cdf_lower(type,x,par1,par2); - case 2: return discrete2_cdf_upper(type,x,par1,par2); - case 3: return discrete2_cdf_lower_inv(type,x,par1,par2); - case 4: return discrete2_cdf_upper_inv(type,x,par1,par2); - } -} - -////////////////////////////////////////////////////////////////////// - -int discrete3(int s, int type, unsigned int par1, unsigned int par2, unsigned int par3, double* r) -{ - const gsl_rng_type * T; - gsl_rng * rng; - - gsl_rng_env_setup(); - T = gsl_rng_default; - rng = gsl_rng_alloc(T); - gsl_rng_set(rng,s); - - switch(type) { - case 0: { (*r) = gsl_ran_hypergeometric(rng,par1,par2,par3); break; } - } - - gsl_rng_free(rng); - - return 0; -} - -double discrete3_pdf(int type, double x, unsigned int par1, unsigned int par2, unsigned int par3) -{ - switch (type) { - case 0: return gsl_ran_hypergeometric_pdf(x,par1,par2,par3); - } -} - -double discrete3_cdf_lower(int type, double x, unsigned int par1, unsigned int par2, unsigned int par3) -{ - switch (type) { - case 0: return gsl_cdf_hypergeometric_P(x,par1,par2,par3); - } -} - -double discrete3_cdf_upper(int type, double x, unsigned int par1, unsigned int par2, unsigned int par3) -{ - switch (type) { - case 0: return gsl_cdf_hypergeometric_Q(x,par1,par2,par3); - } -} - -double discrete3_cdf_lower_inv(int type, double X, unsigned int par1, unsigned int par2, unsigned int par3) -{ - switch (type) { - default : return 0; - } -} - -double discrete3_cdf_upper_inv(int type, double X, unsigned int par1, unsigned int par2, unsigned int par3) -{ - switch (type) { - default : return 0; - } -} - -double discrete3_dist(int df, int type, double x, unsigned int par1, unsigned int par2, unsigned int par3) -{ - switch(df) { - case 0: return discrete3_pdf(type,x,par1,par2,par3); - case 1: return discrete3_cdf_lower(type,x,par1,par2,par3); - case 2: return discrete3_cdf_upper(type,x,par1,par2,par3); - case 3: return discrete3_cdf_lower_inv(type,x,par1,par2,par3); - case 4: return discrete3_cdf_upper_inv(type,x,par1,par2,par3); - } -} - -////////////////////////////////////////////////////////////////////// - -int discrete_mp(int s, int type, unsigned int n, int ps, const double* p, int rs, unsigned int* r) -{ - if (ps != rs) return 2000; // BAD_SIZE - - const gsl_rng_type * T; - gsl_rng * rng; - - gsl_rng_env_setup(); - T = gsl_rng_default; - rng = gsl_rng_alloc(T); - gsl_rng_set(rng,s); - - switch(type) { - case 0: { gsl_ran_multinomial(rng,ps,n,p,r); break; } - } - - gsl_rng_free(rng); - - return 0; -} - -double discrete_mp_pdf(int type, int ps, const double* p, int qs, const unsigned int* q) -{ - switch (type) { - case 0: return gsl_ran_multinomial_pdf(ps,p,q); - } - return 0; -} - -double discrete_mp_cdf_lower(int type, int ps, const double* p, int qs, const unsigned int* q) -{ - switch (type) { - default : return 0; - } -} - -double discrete_mp_cdf_upper(int type, int ps, const double* p, int qs, const unsigned int* q) -{ - switch (type) { - default : return 0; - } -} - -double discrete_mp_cdf_lower_inv(int type, int ps, const double* p, int qs, const unsigned int* q) -{ - switch (type) { - default : return 0; - } -} - -double discrete_mp_cdf_upper_inv(int type, int ps, const double* p, int qs, const unsigned int* q) -{ - switch (type) { - default : return 0; - } -} - - int discrete_mp_dist(int df, int type, double* r, int ps, const double* p, int qs, const unsigned int* q) -{ - switch(df) { - case 0: { (*r) = discrete_mp_pdf(type,ps,p,qs,q); break; } - case 1: { (*r) = discrete_mp_cdf_lower(type,ps,p,qs,q); break; } - case 2: { (*r) = discrete_mp_cdf_upper(type,ps,p,qs,q); break; } - case 3: { (*r) = discrete_mp_cdf_lower_inv(type,ps,p,qs,q); break; } - case 4: { (*r) = discrete_mp_cdf_upper_inv(type,ps,p,qs,q); break; } - } - return 0; -} - -////////////////////////////////////////////////////////////////////// - diff --git a/lib/Numeric/GSL/fitting-aux.c b/lib/Numeric/GSL/fitting-aux.c deleted file mode 100644 --- a/lib/Numeric/GSL/fitting-aux.c +++ /dev/null @@ -1,87 +0,0 @@ -#include -#include - -#include -#include - -int linear(double* c0, double* c1, double* chi_sq, - double* cov00, double* cov01, double* cov11, - int xs, const double* x, int ys, const double* y) -{ - if (xs != ys) return 2000; //BAD_SIZE - - return gsl_fit_linear(x,1,y,1,xs,c0,c1,cov00,cov01,cov11,chi_sq); -} - -int linear_weighted(double* c0, double* c1, double* chi_sq, - double* cov00, double* cov01, double* cov11, - int xs, const double* x, - int ws, const double * w, - int ys, const double* y) -{ - if (xs != ys || xs != ws) return 2000; //BAD_SIZE - - return gsl_fit_wlinear(x,1,w,1,y,1,xs,c0,c1,cov00,cov01,cov11,chi_sq); -} - - -int linear_estimate(double x, double c0, double c1, - double cov00, double cov01, double cov11, - double* y, double* e) -{ - return gsl_fit_linear_est(x,c0,c1,cov00,cov01,cov11,y,e); -} - -int multifit(double* chi_sq, - int xrs, int xcs, const double* x, - int ys, const double* y, - int cs, double* c, - int covrs, int covcs, double* cov) -{ - gsl_multifit_linear_workspace* wsp = gsl_multifit_linear_alloc(xrs,xcs); - - gsl_matrix_const_view X = gsl_matrix_const_view_array(x,xrs,xcs); - gsl_vector_const_view Y = gsl_vector_const_view_array(y,ys); - gsl_vector_view C = gsl_vector_view_array(c,cs); - gsl_matrix_view V = gsl_matrix_view_array(cov,covrs,covcs); - - int err = gsl_multifit_linear(&X.matrix,&Y.vector,&C.vector,&V.matrix,chi_sq,wsp); - - gsl_multifit_linear_free(wsp); - - return err; -} - -int multifit_weighted(double* chi_sq, - int xrs, int xcs, const double* x, - int ws, const double* w, - int ys, const double* y, - int cs, double* c, - int covrs, int covcs, double* cov) -{ - gsl_multifit_linear_workspace* wsp = gsl_multifit_linear_alloc(xrs,xcs); - - gsl_matrix_const_view X = gsl_matrix_const_view_array(x,xrs,xcs); - gsl_vector_const_view W = gsl_vector_const_view_array(w,ws); - gsl_vector_const_view Y = gsl_vector_const_view_array(y,ys); - gsl_vector_view C = gsl_vector_view_array(c,cs); - gsl_matrix_view V = gsl_matrix_view_array(cov,covrs,covcs); - - int err = gsl_multifit_wlinear(&X.matrix,&W.vector,&Y.vector,&C.vector,&V.matrix,chi_sq,wsp); - - gsl_multifit_linear_free(wsp); - - return err; -} - -int multifit_estimate(double* y, double* e, - int xs, const double* x, - int cs, double* c, - int covrs, int covcs, double* cov) -{ - gsl_vector_const_view X = gsl_vector_const_view_array(x,xs); - gsl_vector_const_view C = gsl_vector_const_view_array(c,cs); - gsl_matrix_const_view V = gsl_matrix_const_view_array(cov,covrs,covcs); - - return gsl_multifit_linear_est(&X.vector,&C.vector,&V.matrix,y,e); -} diff --git a/lib/Numeric/GSL/histogram-aux.c b/lib/Numeric/GSL/histogram-aux.c deleted file mode 100644 --- a/lib/Numeric/GSL/histogram-aux.c +++ /dev/null @@ -1,151 +0,0 @@ -#include -#include - -#include -#include - -#include - - -int to_vectors(gsl_histogram * H, int rs, double* r, int bs, double* b) -{ - int i; - for (i = 0; i < rs; i++) - r[i] = H->range[i]; - for (i = 0; i < bs; i++) - b[i] = H->bin[i]; - return 0; -} - -int increment_vector(gsl_histogram* H, int vs, const double* v) -{ - int i; - for (i = 0; i < vs; i++) - gsl_histogram_increment(H,v[i]); - return 0; -} - -int accumulate_vector(gsl_histogram* H, int vs, const double* v, int ws, const double* w) -{ - if (vs != ws) return 2000; // BAD_SIZE - int i; - for (i = 0; i < vs; i++) - gsl_histogram_accumulate(H,v[i],w[i]); - return 0; -} - -int hist_fwrite(const char* filename, const gsl_histogram* h) -{ - int err; - FILE* f = fopen(filename,"w"); - if (!f) return 2003; // BAD_FILE - err = gsl_histogram_fwrite(f,h); - fclose(f); - return err; -} - -int hist_fread(const char* filename, gsl_histogram* h) -{ - int err; - FILE* f = fopen(filename,"r"); - if (!f) return 2003; // BAD_FILE - err = gsl_histogram_fread(f,h); - fclose(f); - return err; -} - -int hist_fprintf(const char* filename, const char* rfmt, const char* bfmt, const gsl_histogram* h) -{ - int err; - FILE* f = fopen(filename,"w"); - if (!f) return 2003; // BAD_FILE - err = gsl_histogram_fprintf(f,h,rfmt,bfmt); - fclose(f); - return err; -} - -int hist_fscanf(const char* filename, gsl_histogram* h) -{ - int err; - FILE* f = fopen(filename,"r"); - if (!f) return 2003; //BAD_FILE - err = gsl_histogram_fscanf(f,h); - fclose(f); - return err; -} - -////////////////////// - -int to_matrix(gsl_histogram2d * H, int rxs, double* rx, int rys, double* ry, int bx, int by, double* b) -{ - int bz = (rxs-1)*(rys-1); - if (bx*by != bz) return 2000; // BAD_SIZE - int i,j; - for (i = 0; i < rxs; i++) - rx[i] = H->xrange[i]; - for (i = 0; i < rys; i++) - ry[i] = H->yrange[i]; - for (i = 0; i < bz; i++) - b[i] = H->bin[i]; - return 0; -} - -int increment_matrix(gsl_histogram2d* H, int xs, const double* x, int ys, const double* y) -{ - if (xs != ys) return 2000; // BAD_SIZE - int i; - for (i = 0; i < xs; i++) - gsl_histogram2d_increment(H,x[i],y[i]); - return 0; -} - -int accumulate_matrix(gsl_histogram2d* H, int xs, const double* x, int ys, const double* y, int ws, const double* w) -{ - if (xs != ys) return 2000; // BAD_SIZE - if (xs != ws) return 2000; // BAD_SIZE - int i; - for (i = 0; i < xs; i++) - gsl_histogram2d_accumulate(H,x[i],y[i],w[i]); - return 0; -} - -int hist2d_fwrite(const char* filename, const gsl_histogram2d* h) -{ - int err; - FILE* f = fopen(filename,"w"); - if (!f) return 2003; // BAD_FILE - err = gsl_histogram2d_fwrite(f,h); - fclose(f); - return err; -} - -int hist2d_fread(const char* filename, gsl_histogram2d* h) -{ - int err; - FILE* f = fopen(filename,"r"); - if (!f) return 2003; // BAD_FILE - err = gsl_histogram2d_fread(f,h); - fclose(f); - return err; -} - -int hist2d_fprintf(const char* filename, const char* rfmt, const char* bfmt, const gsl_histogram2d* h) -{ - int err; - FILE* f = fopen(filename,"w"); - if (!f) return 2003; // BAD_FILE - err = gsl_histogram2d_fprintf(f,h,rfmt,bfmt); - fclose(f); - return err; -} - -int hist2d_fscanf(const char* filename, gsl_histogram2d* h) -{ - int err; - FILE* f = fopen(filename,"r"); - if (!f) return 2003; //BAD_FILE - err = gsl_histogram2d_fscanf(f,h); - fclose(f); - return err; -} - diff --git a/lib/Numeric/GSL/permutation-aux.c b/lib/Numeric/GSL/permutation-aux.c deleted file mode 100644 --- a/lib/Numeric/GSL/permutation-aux.c +++ /dev/null @@ -1,80 +0,0 @@ -#include -#include - -#include -#include - -int random_permute(int s, gsl_permutation* p) -{ - const gsl_rng_type * T; - gsl_rng * r; - - gsl_rng_env_setup(); - T = gsl_rng_default; - r = gsl_rng_alloc(T); - gsl_rng_set(r,s); - - gsl_ran_shuffle(r, p->data, p->size, 1); - - gsl_rng_free(r); - - return 0; -} - -int permute(const gsl_permutation* p, int vs, const double* v, int rs, double* r) -{ - int i; - for (i = 0; i < vs; i++) - r[i] = v[i]; - return gsl_permute(p->data,r,1,vs); -} - -int permute_inverse(const gsl_permutation* p, int vs, const double* v, int rs, double* r) -{ - int i; - for (i = 0; i < vs; i++) - r[i] = v[i]; - return gsl_permute_inverse(p->data,r,1,vs); -} - -int perm_fwrite(const char* filename, const gsl_permutation* p) -{ - int err; - FILE* f = fopen(filename,"w"); - if (!f) return 2003; // BAD_FILE - err = gsl_permutation_fwrite(f,p); - fclose(f); - return err; -} - -int perm_fread(const char* filename, gsl_permutation* p) -{ - int err; - FILE* f = fopen(filename,"r"); - if (!f) return 2003; // BAD_FILE - err = gsl_permutation_fread(f,p); - fclose(f); - return err; -} - -int perm_fprintf(const char* filename, const char* fmt, const gsl_permutation* p) -{ - int err; - FILE* f = fopen(filename,"w"); - if (!f) return 2003; // BAD_FILE - err = gsl_permutation_fprintf(f,p,fmt); - fclose(f); - return err; -} - -int perm_fscanf(const char* filename, gsl_permutation* p) -{ - int err; - FILE* f = fopen(filename,"r"); - if (!f) return 2003; //BAD_FILE - err = gsl_permutation_fscanf(f,p); - fclose(f); - return err; -} - - diff --git a/lib/Numeric/GSL/sort-aux.c b/lib/Numeric/GSL/sort-aux.c deleted file mode 100644 --- a/lib/Numeric/GSL/sort-aux.c +++ /dev/null @@ -1,10 +0,0 @@ -#include -#include - -int sort(int vs, const double* v, int rs, double* r) -{ - if (vs != rs) return 2000; // BAD_SIZE - int i; - for (i = 0; i < vs; i++) r[i] = v[i]; - gsl_sort(r,1,rs); -} diff --git a/lib/Numeric/GSL/statistics-aux.c b/lib/Numeric/GSL/statistics-aux.c deleted file mode 100644 --- a/lib/Numeric/GSL/statistics-aux.c +++ /dev/null @@ -1,224 +0,0 @@ -#include - -inline int mean(double* r, int vs, const double* v) -{ - (*r) = gsl_stats_mean(v,1,vs); - return 0; -} - -inline int variance(double* r, int vs, const double* v) -{ - (*r) = gsl_stats_variance(v,1,vs); - return 0; -} - -inline int variance_m(double m, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_variance_m(v,1,vs,m); - return 0; -} - -inline int stddev(double* r, int vs, const double* v) -{ - (*r) = gsl_stats_sd(v,1,vs); - return 0; -} - -inline int stddev_m(double m, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_sd_m(v,1,vs,m); - return 0; -} - -inline int tot_sumsq(double* r, int vs, const double* v) -{ - (*r) = gsl_stats_tss(v,1,vs); - return 0; -} - -inline int tot_sumsq_m(double m, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_tss_m(v,1,vs,m); - return 0; -} - -inline int var_with_fixed_m(double m, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_variance_with_fixed_mean(v,1,vs,m); - return 0; -} - -inline int stddev_with_fixed_m(double m, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_sd_with_fixed_mean(v,1,vs,m); - return 0; -} - -inline int absdev(double m, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_absdev(v,1,vs); - return 0; -} - -inline int absdev_m(double m, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_absdev_m(v,1,vs,m); - return 0; -} - -inline int skew(double m, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_skew(v,1,vs); - return 0; -} - -inline int skew_m_sd(double m, double sd, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_skew_m_sd(v,1,vs,m,sd); - return 0; -} - -inline int kurtosis(double m, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_kurtosis(v,1,vs); - return 0; -} - -inline int kurtosis_m_sd(double m, double sd, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_kurtosis_m_sd(v,1,vs,m,sd); - return 0; -} - -inline int lag1_autocorrelation(double* r, int vs, const double* v) -{ - (*r) = gsl_stats_lag1_autocorrelation(v,1,vs); - return 0; -} - -inline int covariance(double* r, int vs, const double* v, int ws, const double* w) -{ - if (vs != ws) return 2000; // BAD_LENGTH - (*r) = gsl_stats_covariance(v,1,w,1,vs); - return 0; -} - -inline int covariance_m(double mv, double mw, double* r, int vs, const double* v, int ws, const double* w) -{ - if (vs != ws) return 2000; // BAD_LENGTH - (*r) = gsl_stats_covariance_m(v,1,w,1,vs,mv,mw); - return 0; -} - -inline int correlation(double* r, int vs, const double* v, int ws, const double* w) -{ - if (vs != ws) return 2000; // BAD_LENGTH - (*r) = gsl_stats_correlation(v,1,w,1,vs); - return 0; -} - - -inline int w_mean(int ws, const double* w, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_wmean(w,1,v,1,vs); - return 0; -} - -inline int w_variance(int ws, const double* w, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_wvariance(w,1,v,1,vs); - return 0; -} - -inline int w_variance_m(double m, int ws, const double* w, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_wvariance_m(w,1,v,1,vs,m); - return 0; -} - -inline int w_stddev(int ws, const double* w, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_wsd(w,1,v,1,vs); - return 0; -} - -inline int w_stddev_m(double m, int ws, const double* w, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_wsd_m(w,1,v,1,vs,m); - return 0; -} - -inline int w_tot_sumsq(int ws, const double* w, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_wtss(w,1,v,1,vs); - return 0; -} - -inline int w_tot_sumsq_m(double m, int ws, const double* w, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_wtss_m(w,1,v,1,vs,m); - return 0; -} - -inline int w_var_with_fixed_m(double m, int ws, const double* w, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_wvariance_with_fixed_mean(w,1,v,1,vs,m); - return 0; -} - -inline int w_stddev_with_fixed_m(double m, int ws, const double* w, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_wsd_with_fixed_mean(w,1,v,1,vs,m); - return 0; -} - -inline int w_absdev(double m, int ws, const double* w, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_wabsdev(w,1,v,1,vs); - return 0; -} - -inline int w_absdev_m(double m, int ws, const double* w, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_wabsdev_m(w,1,v,1,vs,m); - return 0; -} - -inline int w_skew(double m, int ws, const double* w, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_wskew(w,1,v,1,vs); - return 0; -} - -inline int w_skew_m_sd(double m, double sd, int ws, const double* w, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_wskew_m_sd(w,1,v,1,vs,m,sd); - return 0; -} - -inline int w_kurtosis(double m, int ws, const double* w, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_wkurtosis(w,1,v,1,vs); - return 0; -} - -inline int w_kurtosis_m_sd(double m, double sd, int ws, const double* w, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_wkurtosis_m_sd(w,1,v,1,vs,m,sd); - return 0; -} - -inline int median(double* r, int vs, const double* v) -{ - (*r) = gsl_stats_median_from_sorted_data(v,1,vs); - return 0; -} - -inline int quantile(double f, double* r, int vs, const double* v) -{ - (*r) = gsl_stats_quantile_from_sorted_data(v,1,vs,f); - return 0; -} - - - diff --git a/lib/Numeric/Statistics/Shannon.hs b/lib/Numeric/Statistics/Shannon.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/Statistics/Shannon.hs @@ -0,0 +1,38 @@ +{-# OPTIONS_GHC -fglasgow-exts #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Statistics.Shannon +-- Copyright : (c) Alexander Vivian Hugh McPhail 2010 +-- License : GPL-style +-- +-- Maintainer : haskell.vivian.mcphail gmail com +-- Stability : provisional +-- Portability : portable +-- +-- Shannon entropy +-- +----------------------------------------------------------------------------- + +module Numeric.Statistics.Shannon ( + entropy + ) where + + +import Data.Packed.Vector + +import Numeric.GSL.Histogram hiding(sum) + +import Numeric.LinearAlgebra.Algorithms +import Numeric.LinearAlgebra.Interface() + +import Prelude hiding (sum) + +sum x = dot x (constant 1 (dim x)) + +-- | the entropy \sum p_i l\ln{p_i} of a sequence +entropy :: Histogram -- the underlying distribution + -> Vector Double -- the sequence + -> Double -- the entropy +entropy p x = let ps = mapVector (\y -> let Just y' = find p y in getBin p y') x + in sum (ps * log ps) +