diff --git a/CHANGES b/CHANGES new file mode 100644 --- /dev/null +++ b/CHANGES @@ -0,0 +1,2 @@ +0.1.0.1: + renamed from hstatistics 0.1.0.5 diff --git a/INSTALL b/INSTALL new file mode 100644 --- /dev/null +++ b/INSTALL @@ -0,0 +1,35 @@ +----------------------------------------------- + A simple signal processing library for Haskell +----------------------------------------------- + +INSTALLATION + +Recommended method (ok in Ubuntu/Debian systems): + $ cabal install hsignal + +INSTALLATION ON WINDOWS ---------------------------------------- + +1) Install a recent ghc (e.g. ghc-6.10.3) + +2) Install cabal-install. A binary for windows can be obtained from: + + http://www.haskell.org/cabal/release/cabal-install-0.6.2/cabal.exe + + Put it somewhere in the path, for instance in c:\ghc\ghc-6.10.3\bin + +3) Download and uncompress hmatrix-x.y.z.tar.gz from Hackage: + + http://hackage.haskell.org/cgi-bin/hackage-scripts/package/hmatrix + +4) Open a terminal, cd to the hmatrix folder, and run + + > cabal install + +5) Download and uncompress hsignal-x.y.z.tar.gz from Hackage: + + http://hackage.haskell.org/cgi-bin/hackage-scripts/package/hsignal + +6) Open a terminal, cd to the hsignal folder, and run + + > cabal install + diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,2 @@ +Copyright Alberto Ruiz 2006-2007 +GPL license diff --git a/README b/README new file mode 100644 --- /dev/null +++ b/README @@ -0,0 +1,12 @@ +----------------------------------------- + A simple scientific library for Haskell +----------------------------------------- + +INSTALLATION + +See the INSTALL file. + +ACKNOWLEDGEMENTS ----------------------------------------------------- + +I thank ALberto Ruiz for hmatrix, especially for setup files which I +have shamelessly copied and modified diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,5 @@ +#! /usr/bin/env runhaskell + +> import Distribution.Simple + +> main = defaultMainWithHooks autoconfUserHooks diff --git a/configure b/configure new file mode 100644 --- /dev/null +++ b/configure @@ -0,0 +1,3 @@ +#! /bin/sh + +runhaskell configure.hs $* diff --git a/configure.hs b/configure.hs new file mode 100644 --- /dev/null +++ b/configure.hs @@ -0,0 +1,145 @@ +#! /usr/bin/env runhaskell +{- configure.hs for hmatrix-gsl-stats, 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 hmatrix-gsl-stats.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 hmatrix-gsl-stats --configure-option=link:lib1,lib2,lib3,etc." + writeFile "hmatrix-gsl-stats.buildinfo" ("buildable: False\n") + Just ops -> do + putStrLn " OK" + g <- checkCommand $ gsl112 bInfo buildInfo + writeFile "hmatrix-gsl-stats.buildinfo" $ "extra-libraries: " ++ + ops ++ "\n" ++ + if g + then "" + else "cc-options: -DGSL110\n" diff --git a/hmatrix-gsl-stats.cabal b/hmatrix-gsl-stats.cabal new file mode 100644 --- /dev/null +++ b/hmatrix-gsl-stats.cabal @@ -0,0 +1,81 @@ +Name: hmatrix-gsl-stats +Version: 0.1.0.1 +License: GPL +License-file: LICENSE +Copyright: (c) A.V.H. McPhail 2010 +Author: Vivian McPhail +Maintainer: haskell.vivian.mcphail gmail com +Stability: provisional +Homepage: http://code.haskell.org/hmatrix-gsl-stats +Synopsis: GSL Statistics interface +Description: Purely functional interface for statistics based on hmatrix and GSL +Category: Math +tested-with: GHC ==6.12.1 + +cabal-version: >=1.2 + +build-type: Custom + +extra-source-files: configure configure.hs README INSTALL CHANGES +extra-tmp-files: hmatrix-gsl-stats.buildinfo + +flag mkl + description: Link with Intel's MKL optimized libraries. + default: False + +library + + Build-Depends: base >= 3 && < 5, + storable-complex, hmatrix >= 0.9.3 + + Extensions: ForeignFunctionInterface + + 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 + 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 + + ghc-prof-options: -auto + + ghc-options: -Wall -fno-warn-missing-signatures + -fno-warn-orphans + -fno-warn-unused-binds + + if flag(mkl) + if arch(x86_64) + extra-libraries: gsl mkl_lapack mkl_intel_lp64 mkl_sequential mkl_core + else + extra-libraries: gsl mkl_lapack mkl_intel mkl_sequential mkl_core + + if os(OSX) + extra-lib-dirs: /opt/local/lib/ + include-dirs: /opt/local/include/ + extra-libraries: gsl + frameworks: Accelerate + + if os(windows) + extra-libraries: gsl-0 blas lapack + +-- The extra-libraries required for GSL +-- should now be automatically detected by configure(.hs) + + extra-libraries: + extra-lib-dirs: + + source-repository head + type: darcs + location: http://code.haskell.org/hmatrix-gsl-stats + diff --git a/lib/Numeric/GSL/Distribution/Common.hs b/lib/Numeric/GSL/Distribution/Common.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/GSL/Distribution/Common.hs @@ -0,0 +1,29 @@ +{-# 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 new file mode 100644 --- /dev/null +++ b/lib/Numeric/GSL/Distribution/Continuous.hs @@ -0,0 +1,299 @@ +{-# 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 new file mode 100644 --- /dev/null +++ b/lib/Numeric/GSL/Distribution/Discrete.hs @@ -0,0 +1,210 @@ +{-# 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 new file mode 100644 --- /dev/null +++ b/lib/Numeric/GSL/Fitting/Linear.hs @@ -0,0 +1,181 @@ +{-# 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 new file mode 100644 --- /dev/null +++ b/lib/Numeric/GSL/Histogram.hs @@ -0,0 +1,468 @@ +{-# 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 new file mode 100644 --- /dev/null +++ b/lib/Numeric/GSL/Histogram2D.hs @@ -0,0 +1,534 @@ +{-# 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 new file mode 100644 --- /dev/null +++ b/lib/Numeric/GSL/Permutation.hs @@ -0,0 +1,349 @@ +{-# 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 new file mode 100644 --- /dev/null +++ b/lib/Numeric/GSL/Sort.hs @@ -0,0 +1,34 @@ +{-# 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 new file mode 100644 --- /dev/null +++ b/lib/Numeric/GSL/Statistics.hs @@ -0,0 +1,308 @@ +{-# 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 new file mode 100644 --- /dev/null +++ b/lib/Numeric/GSL/distribution-aux.c @@ -0,0 +1,749 @@ +#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 new file mode 100644 --- /dev/null +++ b/lib/Numeric/GSL/fitting-aux.c @@ -0,0 +1,87 @@ +#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 new file mode 100644 --- /dev/null +++ b/lib/Numeric/GSL/histogram-aux.c @@ -0,0 +1,151 @@ +#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 new file mode 100644 --- /dev/null +++ b/lib/Numeric/GSL/permutation-aux.c @@ -0,0 +1,80 @@ +#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 new file mode 100644 --- /dev/null +++ b/lib/Numeric/GSL/sort-aux.c @@ -0,0 +1,10 @@ +#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 new file mode 100644 --- /dev/null +++ b/lib/Numeric/GSL/statistics-aux.c @@ -0,0 +1,224 @@ +#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; +} + + +