hmatrix-gsl-stats (empty) → 0.1.0.1
raw patch · 23 files changed
+3998/−0 lines, 23 filesdep +basedep +hmatrixdep +storable-complexbuild-type:Customsetup-changed
Dependencies added: base, hmatrix, storable-complex
Files
- CHANGES +2/−0
- INSTALL +35/−0
- LICENSE +2/−0
- README +12/−0
- Setup.lhs +5/−0
- configure +3/−0
- configure.hs +145/−0
- hmatrix-gsl-stats.cabal +81/−0
- lib/Numeric/GSL/Distribution/Common.hs +29/−0
- lib/Numeric/GSL/Distribution/Continuous.hs +299/−0
- lib/Numeric/GSL/Distribution/Discrete.hs +210/−0
- lib/Numeric/GSL/Fitting/Linear.hs +181/−0
- lib/Numeric/GSL/Histogram.hs +468/−0
- lib/Numeric/GSL/Histogram2D.hs +534/−0
- lib/Numeric/GSL/Permutation.hs +349/−0
- lib/Numeric/GSL/Sort.hs +34/−0
- lib/Numeric/GSL/Statistics.hs +308/−0
- lib/Numeric/GSL/distribution-aux.c +749/−0
- lib/Numeric/GSL/fitting-aux.c +87/−0
- lib/Numeric/GSL/histogram-aux.c +151/−0
- lib/Numeric/GSL/permutation-aux.c +80/−0
- lib/Numeric/GSL/sort-aux.c +10/−0
- lib/Numeric/GSL/statistics-aux.c +224/−0
+ CHANGES view
@@ -0,0 +1,2 @@+0.1.0.1: + renamed from hstatistics 0.1.0.5
+ INSTALL view
@@ -0,0 +1,35 @@+-----------------------------------------------+ A simple signal processing library for Haskell+-----------------------------------------------++INSTALLATION++Recommended method (ok in Ubuntu/Debian systems):+ $ cabal install hsignal++INSTALLATION ON WINDOWS ----------------------------------------++1) Install a recent ghc (e.g. ghc-6.10.3)++2) Install cabal-install. A binary for windows can be obtained from:++ http://www.haskell.org/cabal/release/cabal-install-0.6.2/cabal.exe++ Put it somewhere in the path, for instance in c:\ghc\ghc-6.10.3\bin++3) Download and uncompress hmatrix-x.y.z.tar.gz from Hackage:++ http://hackage.haskell.org/cgi-bin/hackage-scripts/package/hmatrix++4) Open a terminal, cd to the hmatrix folder, and run++ > cabal install++5) Download and uncompress hsignal-x.y.z.tar.gz from Hackage:++ http://hackage.haskell.org/cgi-bin/hackage-scripts/package/hsignal++6) Open a terminal, cd to the hsignal folder, and run++ > cabal install+
+ LICENSE view
@@ -0,0 +1,2 @@+Copyright Alberto Ruiz 2006-2007+GPL license
+ README view
@@ -0,0 +1,12 @@+-----------------------------------------+ A simple scientific library for Haskell+-----------------------------------------++INSTALLATION++See the INSTALL file.++ACKNOWLEDGEMENTS -----------------------------------------------------++I thank ALberto Ruiz for hmatrix, especially for setup files which I+have shamelessly copied and modified
+ Setup.lhs view
@@ -0,0 +1,5 @@+#! /usr/bin/env runhaskell++> import Distribution.Simple++> main = defaultMainWithHooks autoconfUserHooks
+ configure view
@@ -0,0 +1,3 @@+#! /bin/sh++runhaskell configure.hs $*
+ configure.hs view
@@ -0,0 +1,145 @@+#! /usr/bin/env runhaskell+{- configure.hs for hmatrix-gsl-stats, copied from hmatrix+ ------------------------++ GSL and LAPACK may require auxiliary libraries which depend on OS,+ distribution, and implementation. This script tries to to find out+ the correct link command for your system.+ Suggestions and contributions are welcome.++ By default we try to link -lgsl -llapack. This works in ubuntu/debian,+ both with and without ATLAS.+ If this fails we try different sets of additional libraries which are+ known to work in some systems.++ The desired libraries can also be explicitly given by the user using cabal+ flags (e.g., -fmkl, -faccelerate) or --configure-option=link:lib1,lib2,lib3,...++-}++import System+import System.Directory(createDirectoryIfMissing,getCurrentDirectory,getDirectoryContents)+import Data.List(isPrefixOf, intercalate)+import Distribution.Simple.LocalBuildInfo+import Distribution.Simple.Configure+import Distribution.PackageDescription++-- possible additional dependencies for the desired libs (by default gsl lapack)++opts = [ "" -- Ubuntu/Debian+ , "blas"+ , "blas cblas"+ , "cblas"+ , "gslcblas"+ , "blas gslcblas"+ , "f77blas"+ , "f77blas cblas atlas gcc_s" -- Arch Linux (older version of atlas-lapack)+ , "blas gslcblas gfortran" -- Arch Linux with normal blas and lapack+ ]++-- compile a simple program with symbols from GSL and LAPACK with the given libs+testprog bInfo buildInfo libs fmks =+ "echo \"#include <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 hmatrix-gsl-stats.cabal (using a cabal flags)+ -- and from a posible --configure-option=link:lib1,lib2,lib3+ -- by default the desired libs are gsl lapack.++ let pref = if null (words (base ++ " " ++ auxpref)) then "gsl lapack" else auxpref+ fullOpts = map ((pref++" ")++) opts++ -- create the build directory (used for tmp files) if necessary+ createDirectoryIfMissing True $ buildDir bInfo+ + r <- try bInfo buildInfo base fwks fullOpts+ case r of+ Nothing -> do+ putStrLn " FAIL"+ g <- checkCommand $ gsl bInfo buildInfo+ if g+ then putStrLn " *** Sorry, I can't link LAPACK."+ else putStrLn " *** Sorry, I can't link GSL."+ putStrLn " *** Please make sure that the appropriate -dev packages are installed."+ putStrLn " *** You can also specify the required libraries using"+ putStrLn " *** cabal install hmatrix-gsl-stats --configure-option=link:lib1,lib2,lib3,etc."+ writeFile "hmatrix-gsl-stats.buildinfo" ("buildable: False\n")+ Just ops -> do+ putStrLn " OK"+ g <- checkCommand $ gsl112 bInfo buildInfo+ writeFile "hmatrix-gsl-stats.buildinfo" $ "extra-libraries: " +++ ops ++ "\n" +++ if g+ then ""+ else "cc-options: -DGSL110\n"
+ hmatrix-gsl-stats.cabal view
@@ -0,0 +1,81 @@+Name: hmatrix-gsl-stats+Version: 0.1.0.1+License: GPL+License-file: LICENSE+Copyright: (c) A.V.H. McPhail 2010+Author: Vivian McPhail+Maintainer: haskell.vivian.mcphail <at> gmail <dot> com+Stability: provisional+Homepage: http://code.haskell.org/hmatrix-gsl-stats+Synopsis: GSL Statistics interface+Description: Purely functional interface for statistics based on hmatrix and GSL+Category: Math+tested-with: GHC ==6.12.1++cabal-version: >=1.2++build-type: Custom++extra-source-files: configure configure.hs README INSTALL CHANGES+extra-tmp-files: hmatrix-gsl-stats.buildinfo++flag mkl+ description: Link with Intel's MKL optimized libraries.+ default: False++library++ Build-Depends: base >= 3 && < 5,+ storable-complex, hmatrix >= 0.9.3++ Extensions: ForeignFunctionInterface++ hs-source-dirs: lib+ Exposed-modules: Numeric.GSL.Sort+ Numeric.GSL.Statistics+ Numeric.GSL.Histogram+ Numeric.GSL.Histogram2D+ Numeric.GSL.Permutation+ Numeric.GSL.Distribution.Continuous+ Numeric.GSL.Distribution.Discrete+ Numeric.GSL.Distribution.Common+ Numeric.GSL.Fitting.Linear+ other-modules: + C-sources: lib/Numeric/GSL/statistics-aux.c+ lib/Numeric/GSL/sort-aux.c+ lib/Numeric/GSL/histogram-aux.c+ lib/Numeric/GSL/permutation-aux.c+ lib/Numeric/GSL/distribution-aux.c+ lib/Numeric/GSL/fitting-aux.c++ ghc-prof-options: -auto++ ghc-options: -Wall -fno-warn-missing-signatures+ -fno-warn-orphans+ -fno-warn-unused-binds++ if flag(mkl)+ if arch(x86_64)+ extra-libraries: gsl mkl_lapack mkl_intel_lp64 mkl_sequential mkl_core+ else+ extra-libraries: gsl mkl_lapack mkl_intel mkl_sequential mkl_core++ if os(OSX)+ extra-lib-dirs: /opt/local/lib/+ include-dirs: /opt/local/include/+ extra-libraries: gsl+ frameworks: Accelerate++ if os(windows)+ extra-libraries: gsl-0 blas lapack++-- The extra-libraries required for GSL+-- should now be automatically detected by configure(.hs)++ extra-libraries:+ extra-lib-dirs:++ source-repository head+ type: darcs+ location: http://code.haskell.org/hmatrix-gsl-stats+
+ lib/Numeric/GSL/Distribution/Common.hs view
@@ -0,0 +1,29 @@+{-# OPTIONS_GHC -fglasgow-exts #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.GSL.Distribution.Common+-- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- License : GPL-style+--+-- Maintainer : haskell.vivian.mcphail <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 view
@@ -0,0 +1,299 @@+{-# OPTIONS_GHC -fglasgow-exts #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.GSL.Distribution.Continuous+-- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- License : GPL-style+--+-- Maintainer : haskell.vivian.mcphail <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 view
@@ -0,0 +1,210 @@+{-# OPTIONS_GHC -fglasgow-exts #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.GSL.Distribution.Discrete+-- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- License : GPL-style+--+-- Maintainer : haskell.vivian.mcphail <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 view
@@ -0,0 +1,181 @@+{-# OPTIONS_GHC -fglasgow-exts #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.GSL.Fitting.Linear+-- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- License : GPL-style+--+-- Maintainer : haskell.vivian.mcphail <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 view
@@ -0,0 +1,468 @@+{-# OPTIONS_GHC -fglasgow-exts #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.GSL.Histogram+-- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- License : GPL-style+--+-- Maintainer : haskell.vivian.mcphail <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 view
@@ -0,0 +1,534 @@+{-# OPTIONS_GHC -fglasgow-exts #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.GSL.Histogram2D+-- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- License : GPL-style+--+-- Maintainer : haskell.vivian.mcphail <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 view
@@ -0,0 +1,349 @@+{-# OPTIONS_GHC -fglasgow-exts #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.GSL.Permutation+-- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- License : GPL-style+--+-- Maintainer : haskell.vivian.mcphail <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 view
@@ -0,0 +1,34 @@+{-# OPTIONS_GHC -fglasgow-exts #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.GSL.Sort+-- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- License : GPL-style+--+-- Maintainer : haskell.vivian.mcphail <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 view
@@ -0,0 +1,308 @@+{-# OPTIONS_GHC -fglasgow-exts #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.GSL.Statistics+-- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- License : GPL-style+--+-- Maintainer : haskell.vivian.mcphail <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 view
@@ -0,0 +1,749 @@+#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 view
@@ -0,0 +1,87 @@+#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 view
@@ -0,0 +1,151 @@+#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 view
@@ -0,0 +1,80 @@+#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 view
@@ -0,0 +1,10 @@+#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 view
@@ -0,0 +1,224 @@+#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;+}+++