hstatistics 0.1.0.5 → 0.2.0.1
raw patch · 19 files changed
+46/−3894 lines, 19 filesdep +hmatrix-gsl-statsdep −storable-complex
Dependencies added: hmatrix-gsl-stats
Dependencies removed: storable-complex
Files
- configure +0/−3
- configure.hs +0/−145
- hstatistics.cabal +8/−33
- lib/Numeric/GSL/Distribution/Common.hs +0/−29
- lib/Numeric/GSL/Distribution/Continuous.hs +0/−299
- lib/Numeric/GSL/Distribution/Discrete.hs +0/−210
- lib/Numeric/GSL/Fitting/Linear.hs +0/−181
- lib/Numeric/GSL/Histogram.hs +0/−468
- lib/Numeric/GSL/Histogram2D.hs +0/−534
- lib/Numeric/GSL/Permutation.hs +0/−349
- lib/Numeric/GSL/Sort.hs +0/−34
- lib/Numeric/GSL/Statistics.hs +0/−308
- lib/Numeric/GSL/distribution-aux.c +0/−749
- lib/Numeric/GSL/fitting-aux.c +0/−87
- lib/Numeric/GSL/histogram-aux.c +0/−151
- lib/Numeric/GSL/permutation-aux.c +0/−80
- lib/Numeric/GSL/sort-aux.c +0/−10
- lib/Numeric/GSL/statistics-aux.c +0/−224
- lib/Numeric/Statistics/Shannon.hs +38/−0
− configure
@@ -1,3 +0,0 @@-#! /bin/sh--runhaskell configure.hs $*
− configure.hs
@@ -1,145 +0,0 @@-#! /usr/bin/env runhaskell-{- configure.hs for hstatistics, copied from hmatrix- -------------------------- GSL and LAPACK may require auxiliary libraries which depend on OS,- distribution, and implementation. This script tries to to find out- the correct link command for your system.- Suggestions and contributions are welcome.-- By default we try to link -lgsl -llapack. This works in ubuntu/debian,- both with and without ATLAS.- If this fails we try different sets of additional libraries which are- known to work in some systems.-- The desired libraries can also be explicitly given by the user using cabal- flags (e.g., -fmkl, -faccelerate) or --configure-option=link:lib1,lib2,lib3,...---}--import System-import System.Directory(createDirectoryIfMissing,getCurrentDirectory,getDirectoryContents)-import Data.List(isPrefixOf, intercalate)-import Distribution.Simple.LocalBuildInfo-import Distribution.Simple.Configure-import Distribution.PackageDescription---- possible additional dependencies for the desired libs (by default gsl lapack)--opts = [ "" -- Ubuntu/Debian- , "blas"- , "blas cblas"- , "cblas"- , "gslcblas"- , "blas gslcblas"- , "f77blas"- , "f77blas cblas atlas gcc_s" -- Arch Linux (older version of atlas-lapack)- , "blas gslcblas gfortran" -- Arch Linux with normal blas and lapack- ]---- compile a simple program with symbols from GSL and LAPACK with the given libs-testprog bInfo buildInfo libs fmks =- "echo \"#include <gsl/gsl_sf_gamma.h>\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 <gsl/gsl_sf_gamma.h>\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 <gsl/gsl_sf_exp.h>\nint main(){gsl_sf_exprel_n_CF_e(1,1,0);}\""- ++" > " ++ (buildDir bInfo) ++ "/dummy.c; gcc " - ++ (buildDir bInfo) ++ "/dummy.c "- ++ (join $ ccOptions buildInfo) ++ " "- ++ (join $ cppOptions buildInfo) ++ " "- ++ (join $ map ("-I"++) $ includeDirs buildInfo)- ++" -o " ++ (buildDir bInfo) ++ "/dummy "- ++ (join $ map ("-L"++) $ extraLibDirs buildInfo) ++ " -lgsl -lgslcblas"- ++ " > /dev/null 2> /dev/null"---checkCommand c = (ExitSuccess ==) `fmap` system c---- test different configurations until the first one works-try _ _ _ _ [] = return Nothing-try l i b f (opt:rest) = do- ok <- check l i (b ++ " " ++ opt) f- if ok then return (Just opt)- else try l i b f rest---- read --configure-option=link:lib1,lib2,lib3,etc-linkop = "link:"-getUserLink = concatMap (g . drop (length linkop)) . filter (isPrefixOf linkop)- where g = map cs- cs ',' = ' '- cs x = x--main = do- dir <- getCurrentDirectory- putStrLn $ "Current directory: " ++ dir---- files <- getDirectoryContents dir--- mapM_ putStrLn files-- putStr "Checking foreign libraries..."-- args <- getArgs- Just bInfo <- maybeGetPersistBuildConfig "dist"-- let Just lib = library . localPkgDescr $ bInfo- buildInfo = libBuildInfo lib- base = unwords . extraLibs $ buildInfo- fwks = unwords . frameworks $ buildInfo- auxpref = getUserLink args-- -- We extract the desired libs from hstatistics.cabal (using a cabal flags)- -- and from a posible --configure-option=link:lib1,lib2,lib3- -- by default the desired libs are gsl lapack.-- let pref = if null (words (base ++ " " ++ auxpref)) then "gsl lapack" else auxpref- fullOpts = map ((pref++" ")++) opts-- -- create the build directory (used for tmp files) if necessary- createDirectoryIfMissing True $ buildDir bInfo- - r <- try bInfo buildInfo base fwks fullOpts- case r of- Nothing -> do- putStrLn " FAIL"- g <- checkCommand $ gsl bInfo buildInfo- if g- then putStrLn " *** Sorry, I can't link LAPACK."- else putStrLn " *** Sorry, I can't link GSL."- putStrLn " *** Please make sure that the appropriate -dev packages are installed."- putStrLn " *** You can also specify the required libraries using"- putStrLn " *** cabal install hstatistics --configure-option=link:lib1,lib2,lib3,etc."- writeFile "hstatistics.buildinfo" ("buildable: False\n")- Just ops -> do- putStrLn " OK"- g <- checkCommand $ gsl112 bInfo buildInfo- writeFile "hstatistics.buildinfo" $ "extra-libraries: " ++- ops ++ "\n" ++- if g- then ""- else "cc-options: -DGSL110\n"
hstatistics.cabal view
@@ -1,5 +1,5 @@ Name: hstatistics-Version: 0.1.0.5+Version: 0.2.0.1 License: GPL License-file: LICENSE Copyright: (c) A.V.H. McPhail 2010@@ -8,59 +8,34 @@ Stability: provisional Homepage: http://code.haskell.org/hstatistics Synopsis: Statistics-Description: Purely functional interface for statistics based on hmatrix and GSL+Description: Purely functional interface for statistics based on hmatrix and hmatrix-gsl-stats Category: Math tested-with: GHC ==6.12.1 cabal-version: >=1.2 -build-type: Custom+build-type: Simple -extra-source-files: configure configure.hs README INSTALL CHANGES+extra-source-files: README INSTALL CHANGES extra-tmp-files: hstatistics.buildinfo library Build-Depends: base >= 3 && < 5,- storable-complex, hmatrix >= 0.9.3+ hmatrix >= 0.9.3, hmatrix-gsl-stats >= 0.1.0.1 - Extensions: ForeignFunctionInterface+ Extensions: hs-source-dirs: lib- Exposed-modules: Numeric.GSL.Sort- Numeric.GSL.Statistics- Numeric.GSL.Histogram- Numeric.GSL.Histogram2D- Numeric.GSL.Permutation- Numeric.GSL.Distribution.Continuous- Numeric.GSL.Distribution.Discrete- Numeric.GSL.Distribution.Common- Numeric.GSL.Fitting.Linear+ Exposed-modules: Numeric.Statistics.Shannon other-modules: - C-sources: lib/Numeric/GSL/statistics-aux.c- lib/Numeric/GSL/sort-aux.c- lib/Numeric/GSL/histogram-aux.c- lib/Numeric/GSL/permutation-aux.c- lib/Numeric/GSL/distribution-aux.c- lib/Numeric/GSL/fitting-aux.c+ C-sources: ghc-prof-options: -auto ghc-options: -Wall -fno-warn-missing-signatures -fno-warn-orphans -fno-warn-unused-binds-- if os(OSX)- extra-lib-dirs: /opt/local/lib/- include-dirs: /opt/local/include/- extra-libraries: gsl- frameworks: Accelerate---- The extra-libraries required for GSL--- should now be automatically detected by configure(.hs)-- extra-libraries:- extra-lib-dirs: source-repository head type: darcs
− lib/Numeric/GSL/Distribution/Common.hs
@@ -1,29 +0,0 @@-{-# OPTIONS_GHC -fglasgow-exts #-}--------------------------------------------------------------------------------- |--- Module : Numeric.GSL.Distribution.Common--- Copyright : (c) Alexander Vivian Hugh McPhail 2010--- License : GPL-style------ Maintainer : haskell.vivian.mcphail <at> gmail <dot> 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)-------------------------------------------------------------------------------
− lib/Numeric/GSL/Distribution/Continuous.hs
@@ -1,299 +0,0 @@-{-# OPTIONS_GHC -fglasgow-exts #-}--------------------------------------------------------------------------------- |--- Module : Numeric.GSL.Distribution.Continuous--- Copyright : (c) Alexander Vivian Hugh McPhail 2010--- License : GPL-style------ Maintainer : haskell.vivian.mcphail <at> gmail <dot> 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-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
− lib/Numeric/GSL/Distribution/Discrete.hs
@@ -1,210 +0,0 @@-{-# OPTIONS_GHC -fglasgow-exts #-}--------------------------------------------------------------------------------- |--- Module : Numeric.GSL.Distribution.Discrete--- Copyright : (c) Alexander Vivian Hugh McPhail 2010--- License : GPL-style------ Maintainer : haskell.vivian.mcphail <at> gmail <dot> 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-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
− lib/Numeric/GSL/Fitting/Linear.hs
@@ -1,181 +0,0 @@-{-# OPTIONS_GHC -fglasgow-exts #-}--------------------------------------------------------------------------------- |--- Module : Numeric.GSL.Fitting.Linear--- Copyright : (c) Alexander Vivian Hugh McPhail 2010--- License : GPL-style------ Maintainer : haskell.vivian.mcphail <at> gmail <dot> 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-------------------------------------------------------------------------------
− lib/Numeric/GSL/Histogram.hs
@@ -1,468 +0,0 @@-{-# OPTIONS_GHC -fglasgow-exts #-}--------------------------------------------------------------------------------- |--- Module : Numeric.GSL.Histogram--- Copyright : (c) Alexander Vivian Hugh McPhail 2010--- License : GPL-style------ Maintainer : haskell.vivian.mcphail <at> gmail <dot> 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-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
− lib/Numeric/GSL/Histogram2D.hs
@@ -1,534 +0,0 @@-{-# OPTIONS_GHC -fglasgow-exts #-}--------------------------------------------------------------------------------- |--- Module : Numeric.GSL.Histogram2D--- Copyright : (c) Alexander Vivian Hugh McPhail 2010--- License : GPL-style------ Maintainer : haskell.vivian.mcphail <at> gmail <dot> 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-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
− lib/Numeric/GSL/Permutation.hs
@@ -1,349 +0,0 @@-{-# OPTIONS_GHC -fglasgow-exts #-}--------------------------------------------------------------------------------- |--- Module : Numeric.GSL.Permutation--- Copyright : (c) Alexander Vivian Hugh McPhail 2010--- License : GPL-style------ Maintainer : haskell.vivian.mcphail <at> gmail <dot> 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-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
− lib/Numeric/GSL/Sort.hs
@@ -1,34 +0,0 @@-{-# OPTIONS_GHC -fglasgow-exts #-}--------------------------------------------------------------------------------- |--- Module : Numeric.GSL.Sort--- Copyright : (c) Alexander Vivian Hugh McPhail 2010--- License : GPL-style------ Maintainer : haskell.vivian.mcphail <at> gmail <dot> 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
− lib/Numeric/GSL/Statistics.hs
@@ -1,308 +0,0 @@-{-# OPTIONS_GHC -fglasgow-exts #-}--------------------------------------------------------------------------------- |--- Module : Numeric.GSL.Statistics--- Copyright : (c) Alexander Vivian Hugh McPhail 2010--- License : GPL-style------ Maintainer : haskell.vivian.mcphail <at> gmail <dot> 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--------------------------------------------------------------------------------
− lib/Numeric/GSL/distribution-aux.c
@@ -1,749 +0,0 @@-#include <gsl/gsl_rng.h>-#include <gsl/gsl_randist.h>-#include <gsl/gsl_cdf.h>--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;-}--//////////////////////////////////////////////////////////////////////-
− lib/Numeric/GSL/fitting-aux.c
@@ -1,87 +0,0 @@-#include <gsl/gsl_fit.h>-#include <gsl/gsl_multifit.h>--#include <gsl/gsl_vector.h>-#include <gsl/gsl_matrix.h>--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);-}
− lib/Numeric/GSL/histogram-aux.c
@@ -1,151 +0,0 @@-#include <gsl/gsl_math.h>-#include <gsl/gsl_vector.h>--#include <gsl/gsl_histogram.h>-#include <gsl/gsl_histogram2d.h>--#include <stdio.h>---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;-}-
− lib/Numeric/GSL/permutation-aux.c
@@ -1,80 +0,0 @@-#include <gsl/gsl_permutation.h>-#include <gsl/gsl_permute.h>- -#include <gsl/gsl_rng.h>-#include <gsl/gsl_randist.h>--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;-}--
− lib/Numeric/GSL/sort-aux.c
@@ -1,10 +0,0 @@-#include <gsl/gsl_vector.h>-#include <gsl/gsl_sort.h>--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);-}
− lib/Numeric/GSL/statistics-aux.c
@@ -1,224 +0,0 @@-#include <gsl/gsl_statistics_double.h>--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;-}---
+ lib/Numeric/Statistics/Shannon.hs view
@@ -0,0 +1,38 @@+{-# OPTIONS_GHC -fglasgow-exts #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.Statistics.Shannon+-- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- License : GPL-style+--+-- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com+-- Stability : provisional+-- Portability : portable+--+-- Shannon entropy+--+-----------------------------------------------------------------------------++module Numeric.Statistics.Shannon (+ entropy+ ) where+++import Data.Packed.Vector++import Numeric.GSL.Histogram hiding(sum)++import Numeric.LinearAlgebra.Algorithms+import Numeric.LinearAlgebra.Interface()++import Prelude hiding (sum)++sum x = dot x (constant 1 (dim x))++-- | the entropy \sum p_i l\ln{p_i} of a sequence+entropy :: Histogram -- the underlying distribution+ -> Vector Double -- the sequence+ -> Double -- the entropy+entropy p x = let ps = mapVector (\y -> let Just y' = find p y in getBin p y') x+ in sum (ps * log ps) +