hmatrix 0.12.0.2 → 0.13.0.0
raw patch · 32 files changed
+431/−1727 lines, 32 filesdep −HUnitdep −QuickCheckdep −randomdep ~arraydep ~basedep ~vectorsetup-changedPVP ok
version bump matches the API change (PVP)
Dependencies removed: HUnit, QuickCheck, random
Dependency ranges changed: array, base, vector
API changes (from Hackage documentation)
- Numeric.Container: (*/) :: Container c e => c e -> e -> c e
- Numeric.Container: (.*) :: Container c e => e -> c e -> c e
- Numeric.Container: (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t
- Numeric.Container: (<|>) :: (Element t, Joinable a b) => a t -> b t -> Matrix t
- Numeric.Container: vectorMax :: (Container Vector t, Element t) => Vector t -> t
- Numeric.Container: vectorMaxIndex :: Vector Double -> Int
- Numeric.Container: vectorMin :: (Container Vector t, Element t) => Vector t -> t
- Numeric.Container: vectorMinIndex :: Vector Double -> Int
- Numeric.LinearAlgebra.Tests: findNaN :: Int -> Bool
- Numeric.LinearAlgebra.Tests: instance Monad (State s)
- Numeric.LinearAlgebra.Tests: instance Monad m => Monad (MaybeT m)
- Numeric.LinearAlgebra.Tests: qCheck :: Testable prop => Int -> prop -> IO ()
- Numeric.LinearAlgebra.Tests: runBenchmarks :: IO ()
- Numeric.LinearAlgebra.Tests: runTests :: Int -> IO ()
Files
- CHANGES +0/−136
- CHANGES.md +141/−0
- INSTALL +0/−40
- INSTALL.md +110/−0
- Setup.lhs +1/−3
- THANKS +0/−107
- THANKS.md +107/−0
- examples/bool.hs +3/−0
- examples/parallel.hs +3/−2
- examples/tests.hs +0/−3
- hmatrix.cabal +9/−23
- lib/Data/Packed/Internal/Matrix.hs +1/−1
- lib/Data/Packed/Internal/Signatures.hs +3/−3
- lib/Data/Packed/Internal/Vector.hs +18/−2
- lib/Numeric/Container.hs +1/−5
- lib/Numeric/ContainerBoot.hs +4/−105
- lib/Numeric/GSL/Differentiation.hs +1/−1
- lib/Numeric/GSL/Fitting.hs +1/−1
- lib/Numeric/GSL/Fourier.hs +1/−1
- lib/Numeric/GSL/Integration.hs +1/−1
- lib/Numeric/GSL/Internal.hs +1/−1
- lib/Numeric/GSL/Minimization.hs +1/−1
- lib/Numeric/GSL/ODE.hs +1/−1
- lib/Numeric/GSL/Polynomials.hs +5/−1
- lib/Numeric/GSL/Root.hs +2/−3
- lib/Numeric/GSL/Vector.hs +1/−1
- lib/Numeric/LinearAlgebra/Algorithms.hs +14/−7
- lib/Numeric/LinearAlgebra/LAPACK.hs +1/−1
- lib/Numeric/LinearAlgebra/Tests.hs +0/−723
- lib/Numeric/LinearAlgebra/Tests/Instances.hs +0/−249
- lib/Numeric/LinearAlgebra/Tests/Properties.hs +0/−272
- lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h +0/−33
− CHANGES
@@ -1,136 +0,0 @@-0.11.2.0-========--- geigSH' (symmetric generalized eigensystem)--- mapVectorWithIndex---0.11.1.0-========--- exported Mul--- mapMatrixWithIndex{,M,M_}--0.11.0.0-========--- flag -fvector default = True--- invlndet (inverse and log of determinant)--- step, cond--- find--- assoc, accum--0.10.0.0-========--- Module reorganization--- Support for Float and Complex Float elements (excluding LAPACK computations)--- Binary instances for Vector and Matrix--- optimiseMult--- mapVectorM, mapVectorWithIndexM, unzipVectorWith, and related functions.--- diagRect admits diagonal vectors of any length without producing an error,- and takes an additional argument for the off-diagonal elements.--- different signatures in some functions--0.9.3.0-=======--- flag -fvector to optionally use Data.Vector.Storable.Vector- without any conversion.--- Simpler module structure.--- toBlocks, toBlocksEvery--- cholSolve, mbCholSH--- GSL Nonlinear Least-Squares fitting using Levenberg-Marquardt.--- GSL special functions moved to separate package hmatrix-special.--- Added offset of Vector, allowing fast, noncopy subVector (slice).- Vector is now identical to Roman Leshchinskiy's Data.Vector.Storable.Vector,- so we can convert from/to them in O(1).--- Removed Data.Packed.Convert, see examples/vector.hs--0.8.3.0-=======--- odeSolve--- Matrix arithmetic automatically replicates matrix with single row/column--- latexFormat, dispcf--0.8.2.0-=======--- fromRows/fromColumns now automatically expand vectors of dim 1- to match the common dimension.- fromBlocks also replicates single row/column matrices.- Previously all dimensions had to be exactly the same.--- display utilities: dispf, disps, vecdisp--- scalar--- minimizeV, minimizeVD, using Vector instead of lists.--0.8.1.0-=======--- runBenchmarks--0.8.0.0-=======--- singularValues, fullSVD, thinSVD, compactSVD, leftSV, rightSV- and complete interface to [d|z]gesdd.- Algorithms based on the SVD of large matrices can now be- significantly faster.--- eigenvalues, eigenvaluesSH--- linearSolveLS, rq--0.7.2.0-=======--- ranksv--0.7.1.0-=======--- buildVector/buildMatrix--- removed NFData instances--0.6.0.0-=======--- added randomVector, gaussianSample, uniformSample, meanCov--- added rankSVD, nullspaceSVD--- rank, nullspacePrec, and economy svd defined in terms of ranksvd.--- economy svd now admits zero rank matrices and return a "degenerate- rank 1" decomposition with zero singular value.--- added NFData instances for Matrix and Vector.--- liftVector, liftVector2 replaced by mapVector, zipVector.-
+ CHANGES.md view
@@ -0,0 +1,141 @@+0.13.0.0+--------++- tests moved to new package hmatrix-tests++0.11.2.0+--------++- geigSH' (symmetric generalized eigensystem)++- mapVectorWithIndex+++0.11.1.0+--------++- exported Mul++- mapMatrixWithIndex{,M,M_}++0.11.0.0+--------++- flag -fvector default = True++- invlndet (inverse and log of determinant)++- step, cond++- find++- assoc, accum++0.10.0.0+--------++- Module reorganization++- Support for Float and Complex Float elements (excluding LAPACK computations)++- Binary instances for Vector and Matrix++- optimiseMult++- mapVectorM, mapVectorWithIndexM, unzipVectorWith, and related functions.++- diagRect admits diagonal vectors of any length without producing an error,+ and takes an additional argument for the off-diagonal elements.++- different signatures in some functions++0.9.3.0+--------++- flag -fvector to optionally use Data.Vector.Storable.Vector+ without any conversion.++- Simpler module structure.++- toBlocks, toBlocksEvery++- cholSolve, mbCholSH++- GSL Nonlinear Least-Squares fitting using Levenberg-Marquardt.++- GSL special functions moved to separate package hmatrix-special.++- Added offset of Vector, allowing fast, noncopy subVector (slice).+ Vector is now identical to Roman Leshchinskiy's Data.Vector.Storable.Vector,+ so we can convert from/to them in O(1).++- Removed Data.Packed.Convert, see examples/vector.hs++0.8.3.0+--------++- odeSolve++- Matrix arithmetic automatically replicates matrix with single row/column++- latexFormat, dispcf++0.8.2.0+--------++- fromRows/fromColumns now automatically expand vectors of dim 1+ to match the common dimension.+ fromBlocks also replicates single row/column matrices.+ Previously all dimensions had to be exactly the same.++- display utilities: dispf, disps, vecdisp++- scalar++- minimizeV, minimizeVD, using Vector instead of lists.++0.8.1.0+--------++- runBenchmarks++0.8.0.0+--------++- singularValues, fullSVD, thinSVD, compactSVD, leftSV, rightSV+ and complete interface to [d|z]gesdd.+ Algorithms based on the SVD of large matrices can now be+ significantly faster.++- eigenvalues, eigenvaluesSH++- linearSolveLS, rq++0.7.2.0+--------++- ranksv++0.7.1.0+--------++- buildVector/buildMatrix++- removed NFData instances++0.6.0.0+--------++- added randomVector, gaussianSample, uniformSample, meanCov++- added rankSVD, nullspaceSVD++- rank, nullspacePrec, and economy svd defined in terms of ranksvd.++- economy svd now admits zero rank matrices and return a "degenerate+ rank 1" decomposition with zero singular value.++- added NFData instances for Matrix and Vector.++- liftVector, liftVector2 replaced by mapVector, zipVector.+
− INSTALL
@@ -1,40 +0,0 @@-INSTALLATION--Recommended method (tested on Ubuntu/Debian systems):- $ sudo apt-get install libgsl0-dev liblapack-dev- $ cabal install hmatrix--Detailed installation instructions:- http://code.haskell.org/hmatrix/install.html--INSTALLATION ON WINDOWS--1) Install the Haskell Platform (tested on 2011.2.0.1)-- > cabal update--2) Download and unzip the following file into a stable folder %GSL%-- http://perception.inf.um.es/hmatrix/gsl-lapack-windows.zip--3.a) In a msys shell installation should be fully automatic:-- $ cabal install hmatrix --extra-lib-dir=${GSL} --extra-include-dir=${GSL}--3.b) Alternatively, in a normal windows cmd:-- > cabal unpack hmatrix-- Edit hmatrix.cabal, in line 28 change build-type to "Simple", and then-- > cabal install --extra-lib-dir=%GSL% --extra-include-dir=%GSL%-- It may be necessary to put the dlls in the search path.--4) If everything is ok we can run the tests:-- > ghci- Prelude> Numeric.LinearAlgebra.Tests.runTests 20--NOTE: The examples using graphics do not yet work in windows.-
+ INSTALL.md view
@@ -0,0 +1,110 @@++# [hmatrix][hmatrix2] installation++This package requires the [Glasgow Haskell Compiler](http://www.haskell.org/ghc/index.html) ghc >= 6.10, and [cabal-install](http://www.haskell.org/haskellwiki/Cabal-Install), conveniently available in the [Haskell Platform](http://hackage.haskell.org/platform), and the development packages for [GSL](http://www.gnu.org/software/gsl) and BLAS/[LAPACK](http://www.netlib.org/lapack). (The graphical functions also require **gnuplot** and **imagemagick**.)++[hmatrix]: http://code.haskell.org/hmatrix+[hmatrix2]: http://perception.inf.um.es/hmatrix+++## Linux ##################################################+++Ubuntu/Debian:++ $ sudo apt-get install libgsl0-dev liblapack-dev+ $ cabal install hmatrix++Arch Linux: If the automatic installation from Hackage fails, install atlas-lapack and gsl, unpack the source, change the build-type to Simple in hmatrix.cabal (line 28) and add extra-libraries: gsl lapack (line 194). ++Other distributions may require additional libraries. They can be given in a **--configure-option**.++## Mac OS/X ###############################################+++GSL must be installed via MacPorts:++ $ sudo port install gsl-devel +universal+ $ cabal install hmatrix++(Contributed by Heinrich Apfelmus and Torsten Kemps-Benedix).++## Windows ###############################################++We use this [GSL binary](http://www.miscdebris.net/blog/2009/04/20/mingw-345-binaries-of-gnu-scientific-library-112-for-use-with-mingw-and-visual-c/), and blas/lapack dlls built with g77 (contributed by Gilberto Camara). All required files are in [gsl-lapack-windows.zip][winpack].+++1) Install the Haskell Platform (tested on 2011.2.0.1)++ > cabal update++2) Download and unzip the following file into a stable folder %GSL%++ http://perception.inf.um.es/hmatrix/gsl-lapack-windows.zip++3.a) In a msys shell the installation should be fully automatic:++ $ cabal install hmatrix --extra-lib-dir=${GSL} --extra-include-dir=${GSL}++3.b) Alternatively, in a normal windows cmd:++ > cabal unpack hmatrix++ Edit hmatrix.cabal, in line 28 change build-type to "Simple", and then++ > cabal install --extra-lib-dir=%GSL% --extra-include-dir=%GSL%++ It may be necessary to put the dlls in the search path.+++NOTE: The examples using graphics do not yet work in windows.++[install]: http://code.haskell.org/hmatrix/INSTALL+[install2]: http://patch-tag.com/r/aruiz/hmatrix/snapshot/current/content/pretty/INSTALL+[winpack]: http://perception.inf.um.es/hmatrix/gsl-lapack-windows.zip+++## Tests ###############################################++After installation we must verify that the library works as expected:++ $ cabal install hmatrix-tests --enable-tests+ $ ghci+ > Numeric.LinearAlgebra.Tests.runTests 20+ OK, passed 100 tests. + OK, passed 100 tests.+ ... etc...++If you get any failure please run lapack's own tests to confirm that your version is not broken. For instance, in ubuntu 9.04, **libatlas-sse2** does not work (see this [bug report](https://bugs.launchpad.net/ubuntu/+source/atlas/+bug/368478)). If your lapack library is ok but hmatrix's tests fail please send a bug report!+++## Optimized BLAS/LAPACK ##########################################++I have successfully tested ATLAS and MKL on Linux.++### [ATLAS](http://math-atlas.sourceforge.net/) ####################++In Ubuntu >= 9.04 we need:++ $ sudo apt-get install libatlas-base-dev++In older Ubuntu/Debian versions we needed:++ $ sudo apt-get install refblas3-dev lapack3-dev atlas3-base-dev++We may use a version (sse2, 3dnow, etc.) optimized for the machine.++### Intel's MKL ###############################################++There is a free noncommercial download available from Intel's website. To use it I have added the following lines in my .bashrc configuration file:++ export LD_LIBRARY_PATH=/path/to/mkl/lib/arch+ export LIBRARY_PATH=/path/to/mkl/lib/arch++where arch = 32 or em64t.++The library must be installed with the -fmkl flag:++ $ cabal install hmatrix -fmkl++
Setup.lhs view
@@ -8,13 +8,11 @@ > import System.Process(system) > import Config(config) -> main = defaultMainWithHooks simpleUserHooks { confHook = c, runTests = t }+> main = defaultMainWithHooks simpleUserHooks { confHook = c } > c x y = do > binfo <- confHook simpleUserHooks x y > pbi <- config binfo > let pkg_descr = localPkgDescr binfo > return $ binfo { localPkgDescr = updatePackageDescription pbi pkg_descr }--> t _ _ _ _ = system ( "runhaskell examples/tests.hs") >> return()
− THANKS
@@ -1,107 +0,0 @@-I thank Don Stewart, Henning Thielemann, Bulat Ziganshin, Heinrich Apfelmus,-and all the people in the Haskell mailing lists for their help.--I am particularly grateful to Vivian McPhail for his excellent-contributions: improved configure.hs, Binary instances for-Vector and Matrix, support for Float and Complex Float elements,-module reorganization, monadic mapVectorM, and many other improvements.--- Nico Mahlo discovered a bug in the eigendecomposition wrapper.--- Frederik Eaton discovered a bug in the design of the wrappers.--- Eric Kidd has created a wiki page explaining the installation on MacOS X:- http://www.haskell.org/haskellwiki/GSLHaskell_on_MacOS_X--- Fawzi Mohamed discovered a portability bug in the lapack wrappers.--- Pedro E. López de Teruel fixed the interface to lapack.--- Antti Siira discovered a bug in the plotting functions.--- Paulo Tanimoto helped to fix the configuration of the required libraries.- He also discovered the segfault of minimize.hs in ghci.--- Xiao-Yong Jin reported a bug on x86_64 caused by the assumptions in f2c.h,- which are wrong for this architecture.--- Jason Schroeder reported an error in the documentation.--- Bulat Ziganshin gave invaluable help for the ST monad interface to- in-place modifications.--- Don Stewart fixed the implementation of the internal data structures- to achieve excellent, C-like performance in Haskell functions which- explicitly work with the elements of vectors and matrices.--- Dylan Alex Simon improved the numeric instances to allow optimized- implementations of signum and abs on Vectors.--- Pedro E. López de Teruel discovered the need of asm("finit") to- avoid the wrong NaNs produced by foreign functions.--- Reiner Pope added support for luSolve, based on (d|z)getrs.--- Simon Beaumont reported the need of QuickCheck<2 and the invalid- asm("finit") on ppc. He also contributed the configuration options- for the accelerate framework on OS X.--- Daniel Schüssler added compatibility with QuickCheck 2 as well- as QuickCheck 1 using the C preprocessor. He also added some- implementations for the new "shrink" method of class Arbitrary.--- Tracy Wadleigh improved the definitions of (|>) and (><), which now- apply an appropriate 'take' to the given lists so that they may be- safely used on lists that are too long (or infinite).--- Chris Waterson improved the configure.hs program for OS/X.--- Erik de Castro Lopo added buildVector and buildMatrix, which take a- size parameter(s) and a function that maps vector/matrix indices- to the values at that position.--- Jean-Francois Tremblay discovered an error in the tutorial.--- Gilberto Camara contributed improved blas and lapack dlls for Windows.--- Heinrich Apfelmus fixed hmatrix.cabal for OS/X. He also tested the package- on PPC discovering a problem in zgesdd.--- Felipe Lessa tested the performance of GSL special function bindings- and contributed the cabal flag "safe-cheap".--- Ozgur Akgun suggested better symbols for the Bound constructors in the- Linear Programming package.--- Tim Sears reported the zgesdd problem also in intel mac.--- Max Suica simplified the installation on Windows and improved the instructions.--- John Billings first reported an incompatibility with QuickCheck>=2.1.1--- Alexey Khudyakov cleaned up PRAGMAS and fixed some hlint suggestions.--- Torsten Kemps-Benedix reported an installation problem in OS/X.--- Stefan Kersten fixed hmatrix.cabal for 64-bit ghc-7 in OS/X--- Sacha Sokoloski reported an installation problem on Arch Linux and- helped with the configuration.--- Carter Schonwald helped with the configuration for Homebrew OS X and- found a tolerance problem in test "1E5 rots". He also discovered- a bug in the signature of cmap.--- Duncan Coutts reported a problem with configure.hs and contributed- a solution and a simplified Setup.lhs.--- Mark Wright fixed the import of vector >= 0.8.--- Bas van Dijk fixed the import of vector >= 0.8, got rid of some- deprecation warnings and used more explicit imports.--- Tom Nielsen discovered a problem in Config.hs, exposed by link problems- in Ubuntu 11.10 beta.--- Daniel Fischer reported some Haddock markup errors.-
+ THANKS.md view
@@ -0,0 +1,107 @@+I thank Don Stewart, Henning Thielemann, Bulat Ziganshin, Heinrich Apfelmus,+and all the people in the Haskell mailing lists for their help.++I am particularly grateful to Vivian McPhail for his excellent+contributions: improved configure.hs, Binary instances for+Vector and Matrix, support for Float and Complex Float elements,+module reorganization, monadic mapVectorM, and many other improvements.++- Nico Mahlo discovered a bug in the eigendecomposition wrapper.++- Frederik Eaton discovered a bug in the design of the wrappers.++- Eric Kidd has created a wiki page explaining the installation on MacOS X:+ http://www.haskell.org/haskellwiki/GSLHaskell_on_MacOS_X++- Fawzi Mohamed discovered a portability bug in the lapack wrappers.++- Pedro E. López de Teruel fixed the interface to lapack.++- Antti Siira discovered a bug in the plotting functions.++- Paulo Tanimoto helped to fix the configuration of the required libraries.+ He also discovered the segfault of minimize.hs in ghci.++- Xiao-Yong Jin reported a bug on x86_64 caused by the assumptions in f2c.h,+ which are wrong for this architecture.++- Jason Schroeder reported an error in the documentation.++- Bulat Ziganshin gave invaluable help for the ST monad interface to+ in-place modifications.++- Don Stewart fixed the implementation of the internal data structures+ to achieve excellent, C-like performance in Haskell functions which+ explicitly work with the elements of vectors and matrices.++- Dylan Alex Simon improved the numeric instances to allow optimized+ implementations of signum and abs on Vectors.++- Pedro E. López de Teruel discovered the need of asm("finit") to+ avoid the wrong NaNs produced by foreign functions.++- Reiner Pope added support for luSolve, based on (d|z)getrs.++- Simon Beaumont reported the need of QuickCheck<2 and the invalid+ asm("finit") on ppc. He also contributed the configuration options+ for the accelerate framework on OS X.++- Daniel Schüssler added compatibility with QuickCheck 2 as well+ as QuickCheck 1 using the C preprocessor. He also added some+ implementations for the new "shrink" method of class Arbitrary.++- Tracy Wadleigh improved the definitions of (|>) and (><), which now+ apply an appropriate 'take' to the given lists so that they may be+ safely used on lists that are too long (or infinite).++- Chris Waterson improved the configure.hs program for OS/X.++- Erik de Castro Lopo added buildVector and buildMatrix, which take a+ size parameter(s) and a function that maps vector/matrix indices+ to the values at that position.++- Jean-Francois Tremblay discovered an error in the tutorial.++- Gilberto Camara contributed improved blas and lapack dlls for Windows.++- Heinrich Apfelmus fixed hmatrix.cabal for OS/X. He also tested the package+ on PPC discovering a problem in zgesdd.++- Felipe Lessa tested the performance of GSL special function bindings+ and contributed the cabal flag "safe-cheap".++- Ozgur Akgun suggested better symbols for the Bound constructors in the+ Linear Programming package.++- Tim Sears reported the zgesdd problem also in intel mac.++- Max Suica simplified the installation on Windows and improved the instructions.++- John Billings first reported an incompatibility with QuickCheck>=2.1.1++- Alexey Khudyakov cleaned up PRAGMAS and fixed some hlint suggestions.++- Torsten Kemps-Benedix reported an installation problem in OS/X.++- Stefan Kersten fixed hmatrix.cabal for 64-bit ghc-7 in OS/X++- Sacha Sokoloski reported an installation problem on Arch Linux and+ helped with the configuration.++- Carter Schonwald helped with the configuration for Homebrew OS X and+ found a tolerance problem in test "1E5 rots". He also discovered+ a bug in the signature of cmap.++- Duncan Coutts reported a problem with configure.hs and contributed+ a solution and a simplified Setup.lhs.++- Mark Wright fixed the import of vector >= 0.8.++- Bas van Dijk fixed the import of vector >= 0.8, got rid of some+ deprecation warnings, used more explicit imports, and updated to ghc-7.4.++- Tom Nielsen discovered a problem in Config.hs, exposed by link problems+ in Ubuntu 11.10 beta.++- Daniel Fischer reported some Haddock markup errors.+
examples/bool.hs view
@@ -22,6 +22,9 @@ taut x = minElement x == 1 +minEvery a b = cond a b a a b+maxEvery a b = cond a b b b a+ -- examples clip a b x = cond y b y y b where y = cond x a a x x
examples/parallel.hs view
@@ -1,4 +1,5 @@--- $ runhaskell parallel.hs 2000+-- $ ghc --make -O -rtsopts -threaded parallel.hs+-- $ ./parallel 3000 +RTS -N4 -s -A200M import System.Environment(getArgs) import Numeric.LinearAlgebra@@ -7,7 +8,6 @@ inParallel = parMap rwhnf id - -- matrix product decomposed into p parallel subtasks parMul p x y = fromBlocks [ inParallel ( map (x <>) ys ) ] where [ys] = toBlocksEvery (rows y) (cols y `div` p) y@@ -25,3 +25,4 @@ act t1 <- getClockTime print $ tdSec $ normalizeTimeDiff $ diffClockTimes t1 t0+
− examples/tests.hs
@@ -1,3 +0,0 @@-import Numeric.LinearAlgebra.Tests--main = runTests 20
hmatrix.cabal view
@@ -1,5 +1,5 @@ Name: hmatrix-Version: 0.12.0.2+Version: 0.13.0.0 License: GPL License-file: LICENSE Author: Alberto Ruiz@@ -21,17 +21,15 @@ . - "Numeric.LinearAlgebra": everything + instances of standard Haskell numeric classes Category: Math-tested-with: GHC ==6.10.4, GHC ==6.12.1, GHC ==6.12.3, GHC ==7.0.1, GHC==7.0.2, GHC==7.2.1+tested-with: GHC ==7.4 -cabal-version: >=1.6+cabal-version: >=1.8 build-type: Custom -extra-source-files: lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h- Config.hs THANKS INSTALL CHANGES+extra-source-files: Config.hs THANKS.md INSTALL.md CHANGES.md -extra-source-files: examples/tests.hs- examples/deriv.hs+extra-source-files: examples/deriv.hs examples/integrate.hs examples/minimize.hs examples/root.hs@@ -60,10 +58,6 @@ extra-source-files: lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h, lib/Numeric/LinearAlgebra/LAPACK/clapack.h -flag tests- description: Build tests- default: True- flag dd description: svd = zgesdd default: True@@ -152,17 +146,11 @@ Build-Depends: binary cpp-options: -DBINARY - if flag(tests)- Build-Depends: QuickCheck, HUnit, random- exposed-modules: Numeric.LinearAlgebra.Tests- other-modules: Numeric.LinearAlgebra.Tests.Instances,- Numeric.LinearAlgebra.Tests.Properties- -- ghc-prof-options: -auto - ghc-options: -Wall -fno-warn-missing-signatures- -fno-warn-orphans- -fno-warn-unused-binds+ ghc-options: -Wall -fno-warn-missing-signatures+ -fno-warn-orphans+ -fno-warn-unused-binds if flag(unsafe) cpp-options: -DUNSAFE@@ -220,7 +208,5 @@ type: git location: https://github.com/AlbertoRuiz/hmatrix --- Test-Suite tests--- type: exitcode-stdio-1.0--- main-is: examples/tests.hs+-- The tests are in package hmatrix-tests
lib/Data/Packed/Internal/Matrix.hs view
@@ -44,7 +44,7 @@ import Foreign.Ptr(Ptr, castPtr) import Foreign.Storable(Storable, peekElemOff, pokeElemOff, poke, sizeOf) import Data.Complex(Complex)-import Foreign.C.Types(CInt, CChar)+import Foreign.C.Types import Foreign.C.String(newCString) import System.IO.Unsafe(unsafePerformIO)
lib/Data/Packed/Internal/Signatures.hs view
@@ -14,9 +14,9 @@ module Data.Packed.Internal.Signatures where -import Foreign-import Data.Complex-import Foreign.C.Types+import Foreign.Ptr(Ptr)+import Data.Complex(Complex)+import Foreign.C.Types(CInt) type PF = Ptr Float -- type PD = Ptr Double --
lib/Data/Packed/Internal/Vector.hs view
@@ -23,6 +23,7 @@ createVector, vec, asComplex, asReal, float2DoubleV, double2FloatV, stepF, stepD, condF, condD,+ conjugateQ, conjugateC, fwriteVector, freadVector, fprintfVector, fscanfVector, cloneVector, unsafeToForeignPtr,@@ -38,7 +39,7 @@ import Foreign.Ptr(Ptr) import Foreign.Storable(Storable, peekElemOff, pokeElemOff, sizeOf) import Foreign.C.String-import Foreign.C.Types(CInt,CChar)+import Foreign.C.Types import Data.Complex import Control.Monad(when) import System.IO.Unsafe(unsafePerformIO)@@ -344,7 +345,22 @@ foreign import ccall "condF" c_condF :: CInt -> PF -> CInt -> PF -> CInt -> PF -> TFFF foreign import ccall "condD" c_condD :: CInt -> PD -> CInt -> PD -> CInt -> PD -> TVVV -----------------------------------------------------------------+--------------------------------------------------------------------------------++conjugateAux fun x = unsafePerformIO $ do+ v <- createVector (dim x)+ app2 fun vec x vec v "conjugateAux"+ return v++conjugateQ :: Vector (Complex Float) -> Vector (Complex Float)+conjugateQ = conjugateAux c_conjugateQ+foreign import ccall "conjugateQ" c_conjugateQ :: TQVQV++conjugateC :: Vector (Complex Double) -> Vector (Complex Double)+conjugateC = conjugateAux c_conjugateC+foreign import ccall "conjugateC" c_conjugateC :: TCVCV++-------------------------------------------------------------------------------- cloneVector :: Storable t => Vector t -> IO (Vector t) cloneVector v = do
lib/Numeric/Container.hs view
@@ -59,11 +59,7 @@ readMatrix, fscanfVector, fprintfVector, freadVector, fwriteVector, -- * Experimental- build', konst',- -- * Deprecated- (.*),(*/),(<|>),(<->),- vectorMax,vectorMin,- vectorMaxIndex, vectorMinIndex+ build', konst' ) where import Data.Packed
lib/Numeric/ContainerBoot.hs view
@@ -1,3 +1,4 @@+{-# LANGUAGE CPP #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-}@@ -37,11 +38,7 @@ IndexOf, module Data.Complex, -- * Experimental- build', konst',- -- * Deprecated- (.*),(*/),(<|>),(<->),- vectorMax,vectorMin,- vectorMaxIndex, vectorMinIndex+ build', konst' ) where import Data.Packed@@ -49,22 +46,19 @@ import Numeric.Conversion import Data.Packed.Internal import Numeric.GSL.Vector- import Data.Complex import Control.Monad(ap) import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ) -import System.IO.Unsafe- ------------------------------------------------------------------- -type family IndexOf c+type family IndexOf (c :: * -> *) type instance IndexOf Vector = Int type instance IndexOf Matrix = (Int,Int) -type family ArgOf c a+type family ArgOf (c :: * -> *) a type instance ArgOf Vector a = a -> a type instance ArgOf Matrix a = a -> a -> a@@ -512,101 +506,6 @@ type instance ElementOf (Matrix a) = a --------------------------------------------------------------conjugateAux fun x = unsafePerformIO $ do- v <- createVector (dim x)- app2 fun vec x vec v "conjugateAux"- return v--conjugateQ :: Vector (Complex Float) -> Vector (Complex Float)-conjugateQ = conjugateAux c_conjugateQ-foreign import ccall "conjugateQ" c_conjugateQ :: TQVQV--conjugateC :: Vector (Complex Double) -> Vector (Complex Double)-conjugateC = conjugateAux c_conjugateC-foreign import ccall "conjugateC" c_conjugateC :: TCVCV--------------------------------------------------------{-# DEPRECATED (.*) "use scale a x or scalar a * x" #-}---- -- | @x .* a = scale x a@--- (.*) :: (Linear c a) => a -> c a -> c a-infixl 7 .*-a .* x = scale a x--------------------------------------------------------{-# DEPRECATED (*/) "use scale (recip a) x or x / scalar a" #-}---- -- | @a *\/ x = scale (recip x) a@--- (*/) :: (Linear c a) => c a -> a -> c a-infixl 7 */-v */ x = scale (recip x) v-----------------------------------------------------{-# DEPRECATED (<|>) "define operator a & b = fromBlocks[[a,b]] and use asRow/asColumn to join vectors" #-}-{-# DEPRECATED (<->) "define operator a // b = fromBlocks[[a],[b]] and use asRow/asColumn to join vectors" #-}--class Joinable a b where- joinH :: Element t => a t -> b t -> Matrix t- joinV :: Element t => a t -> b t -> Matrix t--instance Joinable Matrix Matrix where- joinH m1 m2 = fromBlocks [[m1,m2]]- joinV m1 m2 = fromBlocks [[m1],[m2]]--instance Joinable Matrix Vector where- joinH m v = joinH m (asColumn v)- joinV m v = joinV m (asRow v)--instance Joinable Vector Matrix where- joinH v m = joinH (asColumn v) m- joinV v m = joinV (asRow v) m--infixl 4 <|>-infixl 3 <->--{-- - | Horizontal concatenation of matrices and vectors:--@> (ident 3 \<-\> 3 * ident 3) \<|\> fromList [1..6.0]-(6><4)- [ 1.0, 0.0, 0.0, 1.0- , 0.0, 1.0, 0.0, 2.0- , 0.0, 0.0, 1.0, 3.0- , 3.0, 0.0, 0.0, 4.0- , 0.0, 3.0, 0.0, 5.0- , 0.0, 0.0, 3.0, 6.0 ]@--}--- (<|>) :: (Element t, Joinable a b) => a t -> b t -> Matrix t-a <|> b = joinH a b---- -- | Vertical concatenation of matrices and vectors.--- (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t-a <-> b = joinV a b-----------------------------------------------------------------------{-# DEPRECATED vectorMin "use minElement" #-}-vectorMin :: (Container Vector t, Element t) => Vector t -> t-vectorMin = minElement--{-# DEPRECATED vectorMax "use maxElement" #-}-vectorMax :: (Container Vector t, Element t) => Vector t -> t-vectorMax = maxElement---{-# DEPRECATED vectorMaxIndex "use minIndex" #-}-vectorMaxIndex :: Vector Double -> Int-vectorMaxIndex = round . toScalarR MaxIdx--{-# DEPRECATED vectorMinIndex "use maxIndex" #-}-vectorMinIndex :: Vector Double -> Int-vectorMinIndex = round . toScalarR MinIdx------------------------------------------------------- class Build f where build' :: BoundsOf f -> f -> ContainerOf f
lib/Numeric/GSL/Differentiation.hs view
@@ -22,7 +22,7 @@ derivBackward ) where -import Foreign.C.Types(CInt)+import Foreign.C.Types import Foreign.Marshal.Alloc(malloc, free) import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr) import Foreign.Storable(peek)
lib/Numeric/GSL/Fitting.hs view
@@ -54,7 +54,7 @@ import Numeric.GSL.Internal import Foreign.Ptr(FunPtr, freeHaskellFunPtr)-import Foreign.C.Types(CInt)+import Foreign.C.Types import System.IO.Unsafe(unsafePerformIO) -------------------------------------------------------------------------
lib/Numeric/GSL/Fourier.hs view
@@ -22,7 +22,7 @@ import Data.Packed.Internal import Data.Complex-import Foreign.C.Types(CInt)+import Foreign.C.Types import System.IO.Unsafe (unsafePerformIO) genfft code v = unsafePerformIO $ do
lib/Numeric/GSL/Integration.hs view
@@ -20,7 +20,7 @@ integrateQAGS ) where -import Foreign.C.Types(CInt)+import Foreign.C.Types import Foreign.Marshal.Alloc(malloc, free) import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr) import Foreign.Storable(peek)
lib/Numeric/GSL/Internal.hs view
@@ -16,7 +16,7 @@ import Foreign.Marshal.Array(copyArray) import Foreign.Ptr(Ptr, FunPtr)-import Foreign.C.Types(CInt)+import Foreign.C.Types import System.IO.Unsafe(unsafePerformIO) iv :: (Vector Double -> Double) -> (CInt -> Ptr Double -> Double)
lib/Numeric/GSL/Minimization.hs view
@@ -65,7 +65,7 @@ import Numeric.GSL.Internal import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)-import Foreign.C.Types(CInt)+import Foreign.C.Types import System.IO.Unsafe(unsafePerformIO) ------------------------------------------------------------------------
lib/Numeric/GSL/ODE.hs view
@@ -36,7 +36,7 @@ import Numeric.GSL.Internal import Foreign.Ptr(FunPtr, nullFunPtr, freeHaskellFunPtr)-import Foreign.C.Types(CInt)+import Foreign.C.Types import System.IO.Unsafe(unsafePerformIO) -------------------------------------------------------------------------
lib/Numeric/GSL/Polynomials.hs view
@@ -1,4 +1,4 @@-{-# LANGUAGE ForeignFunctionInterface #-}+{-# LANGUAGE CPP, ForeignFunctionInterface #-} ----------------------------------------------------------------------------- {- | Module : Numeric.GSL.Polynomials@@ -22,6 +22,10 @@ import Data.Packed.Internal import Data.Complex import System.IO.Unsafe (unsafePerformIO)++#if __GLASGOW_HASKELL__ >= 704+import Foreign.C.Types (CInt(..))+#endif {- | Solution of general polynomial equations, using /gsl_poly_complex_solve/. For example, the three solutions of x^3 + 8 = 0
lib/Numeric/GSL/Root.hs view
@@ -52,9 +52,8 @@ import Data.Packed.Internal import Data.Packed.Matrix import Numeric.GSL.Internal--import Foreign.Ptr(FunPtr, freeHaskellFunPtr)-import Foreign.C.Types(CInt)+import Foreign.Ptr(FunPtr, freeHaskellFunPtr)+import Foreign.C.Types import System.IO.Unsafe(unsafePerformIO) -------------------------------------------------------------------------
lib/Numeric/GSL/Vector.hs view
@@ -31,7 +31,7 @@ import Foreign.Marshal.Alloc(free) import Foreign.Marshal.Array(newArray) import Foreign.Ptr(Ptr)-import Foreign.C.Types(CInt)+import Foreign.C.Types import System.IO.Unsafe(unsafePerformIO) fromei x = fromIntegral (fromEnum x) :: CInt
lib/Numeric/LinearAlgebra/Algorithms.hs view
@@ -3,6 +3,7 @@ {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE TypeFamilies #-}+ ----------------------------------------------------------------------------- {- | Module : Numeric.LinearAlgebra.Algorithms@@ -82,7 +83,7 @@ import Numeric.LinearAlgebra.LAPACK as LAPACK import Data.List(foldl1') import Data.Array-import Numeric.ContainerBoot hiding ((.*),(*/))+import Numeric.ContainerBoot {- | Class used to define generic linear algebra computations for both real and complex matrices. Only double precision is supported in this version (we can@@ -567,6 +568,18 @@ geps delta = head [ k | (k,g) <- epslist, g<delta] ++{- | Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan,+ based on a scaled Pade approximation.+-}+expm :: Field t => Matrix t -> Matrix t+expm = expGolub++expGolub :: ( Fractional t, Element t, Field t+ , Normed Matrix t+ , RealFrac (RealOf t)+ , Floating (RealOf t)+ ) => Matrix t -> Matrix t expGolub m = iterate msq f !! j where j = max 0 $ floor $ logBase 2 $ pnorm Infinity m a = m */ fromIntegral ((2::Int)^j)@@ -586,12 +599,6 @@ v */ x = scale (recip x) v (.*) = scale (|+|) = add--{- | Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan,- based on a scaled Pade approximation.--}-expm :: Field t => Matrix t -> Matrix t-expm = expGolub --------------------------------------------------------------
lib/Numeric/LinearAlgebra/LAPACK.hs view
@@ -47,7 +47,7 @@ import Numeric.GSL.Vector(vectorMapValR, FunCodeSV(Scale)) import Foreign.Ptr(nullPtr)-import Foreign.C.Types(CInt)+import Foreign.C.Types import Control.Monad(when) import System.IO.Unsafe(unsafePerformIO)
− lib/Numeric/LinearAlgebra/Tests.hs
@@ -1,723 +0,0 @@-{-# LANGUAGE CPP #-}-{-# OPTIONS_GHC -fno-warn-unused-imports -fno-warn-incomplete-patterns #-}-------------------------------------------------------------------------------{- |-Module : Numeric.LinearAlgebra.Tests-Copyright : (c) Alberto Ruiz 2007-9-License : GPL-style--Maintainer : Alberto Ruiz (aruiz at um dot es)-Stability : provisional-Portability : portable--Some tests.---}--module Numeric.LinearAlgebra.Tests(--- module Numeric.LinearAlgebra.Tests.Instances,--- module Numeric.LinearAlgebra.Tests.Properties,- qCheck, runTests, runBenchmarks, findNaN---, runBigTests-) where--import Data.Packed.Random-import Numeric.LinearAlgebra-import Numeric.LinearAlgebra.LAPACK-import Numeric.LinearAlgebra.Tests.Instances-import Numeric.LinearAlgebra.Tests.Properties-import Test.HUnit hiding ((~:),test,Testable,State)-import System.Info-import Data.List(foldl1')-import Numeric.GSL-import Prelude hiding ((^))-import qualified Prelude-import System.CPUTime-import Text.Printf-import Data.Packed.Development(unsafeFromForeignPtr,unsafeToForeignPtr)-import Control.Arrow((***))-import Debug.Trace--#include "Tests/quickCheckCompat.h"--debug x = trace (show x) x--a ^ b = a Prelude.^ (b :: Int)--utest str b = TestCase $ assertBool str b--a ~~ b = fromList a |~| fromList b--feye n = flipud (ident n) :: Matrix Double---------------------------------------------------------------detTest1 = det m == 26- && det mc == 38 :+ (-3)- && det (feye 2) == -1- where- m = (3><3) - [ 1, 2, 3- , 4, 5, 7- , 2, 8, 4 :: Double- ]- mc = (3><3)- [ 1, 2, 3- , 4, 5, 7- , 2, 8, i- ]--detTest2 = inv1 |~| inv2 && [det1] ~~ [det2]- where- m = complex (feye 6)- inv1 = inv m- det1 = det m- (inv2,(lda,sa)) = invlndet m- det2 = sa * exp lda------------------------------------------------------------------------polyEval cs x = foldr (\c ac->ac*x+c) 0 cs--polySolveProp p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p))-------------------------------------------------------------------------quad f a b = fst $ integrateQAGS 1E-9 100 f a b---- A multiple integral can be easily defined using partial application-quad2 f a b g1 g2 = quad h a b- where h x = quad (f x) (g1 x) (g2 x)--volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y)) - 0 r (const 0) (\x->sqrt (r*r-x*x))-------------------------------------------------------------------------derivTest = abs (d (\x-> x * d (\y-> x+y) 1) 1 - 1) < 1E-10- where d f x = fst $ derivCentral 0.01 f x--------------------------------------------------------------------------- besselTest = utest "bessel_J0_e" ( abs (r-expected) < e )--- where (r,e) = bessel_J0_e 5.0--- expected = -0.17759677131433830434739701---- exponentialTest = utest "exp_e10_e" ( abs (v*10^e - expected) < 4E-2 )--- where (v,e,_err) = exp_e10_e 30.0--- expected = exp 30.0-------------------------------------------------------------------------nd1 = (3><3) [ 1/2, 1/4, 1/4- , 0/1, 1/2, 1/4- , 1/2, 1/4, 1/2 :: Double]--nd2 = (2><2) [1, 0, 1, 1:: Complex Double]--expmTest1 = expm nd1 :~14~: (3><3)- [ 1.762110887278176- , 0.478085470590435- , 0.478085470590435- , 0.104719410945666- , 1.709751181805343- , 0.425725765117601- , 0.851451530235203- , 0.530445176063267- , 1.814470592751009 ]--expmTest2 = expm nd2 :~15~: (2><2)- [ 2.718281828459045- , 0.000000000000000- , 2.718281828459045- , 2.718281828459045 ]-------------------------------------------------------------------------minimizationTest = TestList- [ utest "minimization conjugatefr" (minim1 f df [5,7] ~~ [1,2])- , utest "minimization nmsimplex2" (minim2 f [5,7] `elem` [24,25])- ]- where f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30- df [x,y] = [20*(x-1), 40*(y-2)]- minim1 g dg ini = fst $ minimizeD ConjugateFR 1E-3 30 1E-2 1E-4 g dg ini- minim2 g ini = rows $ snd $ minimize NMSimplex2 1E-2 30 [1,1] g ini-------------------------------------------------------------------------rootFindingTest = TestList [ utest "root Hybrids" (fst sol1 ~~ [1,1])- , utest "root Newton" (rows (snd sol2) == 2)- ]- where sol1 = root Hybrids 1E-7 30 (rosenbrock 1 10) [-10,-5]- sol2 = rootJ Newton 1E-7 30 (rosenbrock 1 10) (jacobian 1 10) [-10,-5]- rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ]- jacobian a b [x,_y] = [ [-a , 0]- , [-2*b*x, b] ]-------------------------------------------------------------------------odeTest = utest "ode" (last (toLists sol) ~~ [-1.7588880332411019, 8.364348908711941e-2])- where sol = odeSolveV RK8pd 1E-6 1E-6 0 (l2v $ vanderpol 10) Nothing (fromList [1,0]) ts- ts = linspace 101 (0,100)- l2v f = \t -> fromList . f t . toList- vanderpol mu _t [x,y] = [y, -x + mu * y * (1-x^2) ]-------------------------------------------------------------------------fittingTest = utest "levmar" (ok1 && ok2)- where- xs = map return [0 .. 39]- sigma = 0.1- ys = map return $ toList $ fromList (map (head . expModel [5,0.1,1]) xs)- + scalar sigma * (randomVector 0 Gaussian 40)- dats = zip xs (zip ys (repeat sigma))- dat = zip xs ys-- expModel [a,lambda,b] [t] = [a * exp (-lambda * t) + b]- expModelDer [a,lambda,_b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]]-- sols = fst $ fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dats [1,0,0]- sol = fst $ fitModel 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0]-- ok1 = and (zipWith f sols [5,0.1,1]) where f (x,d) r = abs (x-r)<2*d- ok2 = norm2 (fromList (map fst sols) - fromList sol) < 1E-5---------------------------------------------------------mbCholTest = utest "mbCholTest" (ok1 && ok2) where- m1 = (2><2) [2,5,5,8 :: Double]- m2 = (2><2) [3,5,5,9 :: Complex Double]- ok1 = mbCholSH m1 == Nothing- ok2 = mbCholSH m2 == Just (chol m2)-------------------------------------------------------------------------randomTestGaussian = c :~1~: snd (meanCov dat) where- a = (3><3) [1,2,3,- 2,4,0,- -2,2,1]- m = 3 |> [1,2,3]- c = a <> trans a- dat = gaussianSample 7 (10^6) m c--randomTestUniform = c :~1~: snd (meanCov dat) where- c = diag $ 3 |> map ((/12).(^2)) [1,2,3]- dat = uniformSample 7 (10^6) [(0,1),(1,3),(3,6)]-------------------------------------------------------------------------rot :: Double -> Matrix Double-rot a = (3><3) [ c,0,s- , 0,1,0- ,-s,0,c ]- where c = cos a- s = sin a--rotTest = fun (10^5) :~11~: rot 5E4- where fun n = foldl1' (<>) (map rot angles)- where angles = toList $ linspace n (0,1)-------------------------------------------------------------------------- vector <= 0.6.0.2 bug discovered by Patrick Perry--- http://trac.haskell.org/vector/ticket/31--offsetTest = y == y' where- x = fromList [0..3 :: Double]- y = subVector 1 3 x- (f,o,n) = unsafeToForeignPtr y- y' = unsafeFromForeignPtr f o n-------------------------------------------------------------------------normsVTest = TestList [- utest "normv2CD" $ norm2PropC v- , utest "normv2CF" $ norm2PropC (single v)-#ifndef NONORMVTEST- , utest "normv2D" $ norm2PropR x- , utest "normv2F" $ norm2PropR (single x)-#endif- , utest "normv1CD" $ norm1 v == 8- , utest "normv1CF" $ norm1 (single v) == 8- , utest "normv1D" $ norm1 x == 6- , utest "normv1F" $ norm1 (single x) == 6-- , utest "normvInfCD" $ normInf v == 5- , utest "normvInfCF" $ normInf (single v) == 5- , utest "normvInfD" $ normInf x == 3- , utest "normvInfF" $ normInf (single x) == 3-- ] where v = fromList [1,-2,3:+4] :: Vector (Complex Double)- x = fromList [1,2,-3] :: Vector Double-#ifndef NONORMVTEST- norm2PropR a = norm2 a =~= sqrt (dot a a)-#endif- norm2PropC a = norm2 a =~= realPart (sqrt (dot a (conj a)))- a =~= b = fromList [a] |~| fromList [b]--normsMTest = TestList [- utest "norm2mCD" $ pnorm PNorm2 v =~= 8.86164970498005- , utest "norm2mCF" $ pnorm PNorm2 (single v) =~= 8.86164970498005- , utest "norm2mD" $ pnorm PNorm2 x =~= 5.96667765076216- , utest "norm2mF" $ pnorm PNorm2 (single x) =~= 5.96667765076216-- , utest "norm1mCD" $ pnorm PNorm1 v == 9- , utest "norm1mCF" $ pnorm PNorm1 (single v) == 9- , utest "norm1mD" $ pnorm PNorm1 x == 7- , utest "norm1mF" $ pnorm PNorm1 (single x) == 7-- , utest "normmInfCD" $ pnorm Infinity v == 12- , utest "normmInfCF" $ pnorm Infinity (single v) == 12- , utest "normmInfD" $ pnorm Infinity x == 8- , utest "normmInfF" $ pnorm Infinity (single x) == 8-- , utest "normmFroCD" $ pnorm Frobenius v =~= 8.88819441731559- , utest "normmFroCF" $ pnorm Frobenius (single v) =~~= 8.88819441731559- , utest "normmFroD" $ pnorm Frobenius x =~= 6.24499799839840- , utest "normmFroF" $ pnorm Frobenius (single x) =~~= 6.24499799839840-- ] where v = (2><2) [1,-2*i,3:+4,7] :: Matrix (Complex Double)- x = (2><2) [1,2,-3,5] :: Matrix Double- a =~= b = fromList [a] :~10~: fromList [b]- a =~~= b = fromList [a] :~5~: fromList [b]-------------------------------------------------------------------------sumprodTest = TestList [- utest "sumCD" $ sumElements z == 6- , utest "sumCF" $ sumElements (single z) == 6- , utest "sumD" $ sumElements v == 6- , utest "sumF" $ sumElements (single v) == 6-- , utest "prodCD" $ prodProp z- , utest "prodCF" $ prodProp (single z)- , utest "prodD" $ prodProp v- , utest "prodF" $ prodProp (single v)- ] where v = fromList [1,2,3] :: Vector Double- z = fromList [1,2-i,3+i]- prodProp x = prodElements x == product (toList x)-------------------------------------------------------------------------chainTest = utest "chain" $ foldl1' (<>) ms |~| optimiseMult ms where- ms = [ diag (fromList [1,2,3 :: Double])- , konst 3 (3,5)- , (5><10) [1 .. ]- , konst 5 (10,2)- ]-------------------------------------------------------------------------conjuTest m = mapVector conjugate (flatten (trans m)) == flatten (ctrans m)-------------------------------------------------------------------------newtype State s a = State { runState :: s -> (a,s) }--instance Monad (State s) where- return a = State $ \s -> (a,s)- m >>= f = State $ \s -> let (a,s') = runState m s- in runState (f a) s'--state_get :: State s s-state_get = State $ \s -> (s,s)--state_put :: s -> State s ()-state_put s = State $ \_ -> ((),s)--evalState :: State s a -> s -> a-evalState m s = let (a,s') = runState m s- in seq s' a--newtype MaybeT m a = MaybeT { runMaybeT :: m (Maybe a) }--instance Monad m => Monad (MaybeT m) where- return a = MaybeT $ return $ Just a- m >>= f = MaybeT $ do- res <- runMaybeT m- case res of- Nothing -> return Nothing- Just r -> runMaybeT (f r)- fail _ = MaybeT $ return Nothing--lift_maybe m = MaybeT $ do- res <- m- return $ Just res---- | apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs---successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool-successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ stp (subVector 1 (dim v - 1) v))) (v @> 0)- where stp e = do- ep <- lift_maybe $ state_get- if t e ep- then lift_maybe $ state_put e- else (fail "successive_ test failed")---- | operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input---successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b-successive f v = evalState (mapVectorM stp (subVector 1 (dim v - 1) v)) (v @> 0)- where stp e = do- ep <- state_get- state_put e- return $ f ep e---succTest = utest "successive" $- successive_ (>) (fromList [1 :: Double,2,3,4]) == True- && successive_ (>) (fromList [1 :: Double,3,2,4]) == False- && successive (+) (fromList [1..10 :: Double]) == 9 |> [3,5,7,9,11,13,15,17,19]-------------------------------------------------------------------------findAssocTest = utest "findAssoc" ok- where- ok = m1 == m2- m1 = assoc (6,6) 7 $ zip (find (>0) (ident 5 :: Matrix Float)) [10 ..] :: Matrix Double- m2 = diagRect 7 (fromList[10..14]) 6 6-------------------------------------------------------------------------condTest = utest "cond" ok- where- ok = step v * v == cond v 0 0 0 v- v = fromList [-7 .. 7 ] :: Vector Float-------------------------------------------------------------------------conformTest = utest "conform" ok- where- ok = 1 + row [1,2,3] + col [10,20,30,40] + (4><3) [1..]- == (4><3) [13,15,17- ,26,28,30- ,39,41,43- ,52,54,56]- row = asRow . fromList- col = asColumn . fromList :: [Double] -> Matrix Double-------------------------------------------------------------------------accumTest = utest "accum" ok- where- x = ident 3 :: Matrix Double- ok = accum x (+) [((1,2),7), ((2,2),3)]- == (3><3) [1,0,0- ,0,1,7- ,0,0,4]- &&- toList (flatten x) == [1,0,0,0,1,0,0,0,1] --------------------------------------------------------------------------- | All tests must pass with a maximum dimension of about 20--- (some tests may fail with bigger sizes due to precision loss).-runTests :: Int -- ^ maximum dimension- -> IO ()-runTests n = do- setErrorHandlerOff- let test p = qCheck n p- putStrLn "------ mult Double"- test (multProp1 10 . rConsist)- test (multProp1 10 . cConsist)- test (multProp2 10 . rConsist)- test (multProp2 10 . cConsist)- putStrLn "------ mult Float"- test (multProp1 6 . (single *** single) . rConsist)- test (multProp1 6 . (single *** single) . cConsist)- test (multProp2 6 . (single *** single) . rConsist)- test (multProp2 6 . (single *** single) . cConsist)- putStrLn "------ sub-trans"- test (subProp . rM)- test (subProp . cM)- putStrLn "------ ctrans"- test (conjuTest . cM)- test (conjuTest . zM)- putStrLn "------ lu"- test (luProp . rM)- test (luProp . cM)- putStrLn "------ inv (linearSolve)"- test (invProp . rSqWC)- test (invProp . cSqWC)- putStrLn "------ luSolve"- test (linearSolveProp (luSolve.luPacked) . rSqWC)- test (linearSolveProp (luSolve.luPacked) . cSqWC)- putStrLn "------ cholSolve"- test (linearSolveProp (cholSolve.chol) . rPosDef)- test (linearSolveProp (cholSolve.chol) . cPosDef)- putStrLn "------ luSolveLS"- test (linearSolveProp linearSolveLS . rSqWC)- test (linearSolveProp linearSolveLS . cSqWC)- test (linearSolveProp2 linearSolveLS . rConsist)- test (linearSolveProp2 linearSolveLS . cConsist)- putStrLn "------ pinv (linearSolveSVD)"- test (pinvProp . rM)- test (pinvProp . cM)- putStrLn "------ det"- test (detProp . rSqWC)- test (detProp . cSqWC)- putStrLn "------ svd"- test (svdProp1 . rM)- test (svdProp1 . cM)- test (svdProp1a svdR)- test (svdProp1a svdC)- test (svdProp1a svdRd)- test (svdProp1b svdR)- test (svdProp1b svdC)- test (svdProp1b svdRd)- test (svdProp2 thinSVDR)- test (svdProp2 thinSVDC)- test (svdProp2 thinSVDRd)- test (svdProp2 thinSVDCd)- test (svdProp3 . rM)- test (svdProp3 . cM)- test (svdProp4 . rM)- test (svdProp4 . cM)- test (svdProp5a)- test (svdProp5b)- test (svdProp6a)- test (svdProp6b)- test (svdProp7 . rM)- test (svdProp7 . cM)- putStrLn "------ svdCd"-#ifdef NOZGESDD- putStrLn "Omitted"-#else- test (svdProp1a svdCd)- test (svdProp1b svdCd)-#endif- putStrLn "------ eig"- test (eigSHProp . rHer)- test (eigSHProp . cHer)- test (eigProp . rSq)- test (eigProp . cSq)- test (eigSHProp2 . rHer)- test (eigSHProp2 . cHer)- test (eigProp2 . rSq)- test (eigProp2 . cSq)- putStrLn "------ nullSpace"- test (nullspaceProp . rM)- test (nullspaceProp . cM)- putStrLn "------ qr"- test (qrProp . rM)- test (qrProp . cM)- test (rqProp . rM)- test (rqProp . cM)- test (rqProp1 . cM)- test (rqProp2 . cM)- test (rqProp3 . cM)- putStrLn "------ hess"- test (hessProp . rSq)- test (hessProp . cSq)- putStrLn "------ schur"- test (schurProp2 . rSq)- test (schurProp1 . cSq)- putStrLn "------ chol"- test (cholProp . rPosDef)- test (cholProp . cPosDef)- test (exactProp . rPosDef)- test (exactProp . cPosDef)- putStrLn "------ expm"- test (expmDiagProp . complex. rSqWC)- test (expmDiagProp . cSqWC)- putStrLn "------ fft"- test (\v -> ifft (fft v) |~| v)- putStrLn "------ vector operations - Double"- test (\u -> sin u ^ 2 + cos u ^ 2 |~| (1::RM))- test $ (\u -> sin u ^ 2 + cos u ^ 2 |~| (1::CM)) . liftMatrix makeUnitary- test (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM))- test (\u -> cos u * tan u |~| sin (u::RM))- test $ (\u -> cos u * tan u |~| sin (u::CM)) . liftMatrix makeUnitary- putStrLn "------ vector operations - Float"- test (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::FM))- test $ (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::ZM)) . liftMatrix makeUnitary- test (\u -> sin u ** 2 + cos u ** 2 |~~| (1::FM))- test (\u -> cos u * tan u |~~| sin (u::FM))- test $ (\u -> cos u * tan u |~~| sin (u::ZM)) . liftMatrix makeUnitary- putStrLn "------ read . show"- test (\m -> (m::RM) == read (show m))- test (\m -> (m::CM) == read (show m))- test (\m -> toRows (m::RM) == read (show (toRows m)))- test (\m -> toRows (m::CM) == read (show (toRows m)))- test (\m -> (m::FM) == read (show m))- test (\m -> (m::ZM) == read (show m))- test (\m -> toRows (m::FM) == read (show (toRows m)))- test (\m -> toRows (m::ZM) == read (show (toRows m)))- putStrLn "------ some unit tests"- _ <- runTestTT $ TestList- [ utest "1E5 rots" rotTest- , utest "det1" detTest1- , utest "invlndet" detTest2- , utest "expm1" (expmTest1)- , utest "expm2" (expmTest2)- , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM)- , utest "arith2" $ ((scalar (1+i) * ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( scalar (140*i-51) :: CM)- , utest "arith3" $ exp (scalar i * ones(10,10)*pi) + 1 |~| 0- , utest "<\\>" $ (3><2) [2,0,0,3,1,1::Double] <\> 3|>[4,9,5] |~| 2|>[2,3]--- , utest "gamma" (gamma 5 == 24.0)--- , besselTest--- , exponentialTest- , utest "deriv" derivTest- , utest "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < 1E-8)- , utest "polySolve" (polySolveProp [1,2,3,4])- , minimizationTest- , rootFindingTest- , utest "randomGaussian" randomTestGaussian- , utest "randomUniform" randomTestUniform- , utest "buildVector/Matrix" $- complex (10 |> [0::Double ..]) == buildVector 10 fromIntegral- && ident 5 == buildMatrix 5 5 (\(r,c) -> if r==c then 1::Double else 0)- , utest "rank" $ rank ((2><3)[1,0,0,1,6*eps,0]) == 1- && rank ((2><3)[1,0,0,1,7*eps,0]) == 2- , utest "block" $ fromBlocks [[ident 3,0],[0,ident 4]] == (ident 7 :: CM)- , odeTest- , fittingTest- , mbCholTest- , utest "offset" offsetTest- , normsVTest- , normsMTest- , sumprodTest- , chainTest- , succTest- , findAssocTest- , condTest- , conformTest- , accumTest- ]- return ()----- single precision approximate equality-infixl 4 |~~|-a |~~| b = a :~6~: b--makeUnitary v | realPart n > 1 = v / scalar n- | otherwise = v- where n = sqrt (conj v <.> v)---- -- | Some additional tests on big matrices. They take a few minutes.--- runBigTests :: IO ()--- runBigTests = undefined---- testcase for nonempty fpu stack-findNaN :: Int -> Bool-findNaN n = all (bugProp . eye) (take n $ cycle [1..20])- where eye m = ident m :: Matrix ( Double)-------------------------------------------------------------------------------------- | Performance measurements.-runBenchmarks :: IO ()-runBenchmarks = do- --cholBench- solveBench- subBench- multBench- svdBench- eigBench- putStrLn ""------------------------------------time msg act = do- putStr (msg++" ")- t0 <- getCPUTime- act `seq` putStr " "- t1 <- getCPUTime- printf "%6.2f s CPU\n" $ (fromIntegral (t1 - t0) / (10^12 :: Double)) :: IO ()- return ()------------------------------------manymult n = foldl1' (<>) (map rot2 angles) where- angles = toList $ linspace n (0,1)- rot2 :: Double -> Matrix Double- rot2 a = (3><3) [ c,0,s- , 0,1,0- ,-s,0,c ]- where c = cos a- s = sin a--multb n = foldl1' (<>) (replicate (10^6) (ident n :: Matrix Double))------------------------------------subBench = do- putStrLn ""- let g = foldl1' (.) (replicate (10^5) (\v -> subVector 1 (dim v -1) v))- time "0.1M subVector " (g (constant 1 (1+10^5) :: Vector Double) @> 0)- let f = foldl1' (.) (replicate (10^5) (fromRows.toRows))- time "subVector-join 3" (f (ident 3 :: Matrix Double) @@>(0,0))- time "subVector-join 10" (f (ident 10 :: Matrix Double) @@>(0,0))------------------------------------multBench = do- let a = ident 1000 :: Matrix Double- let b = ident 2000 :: Matrix Double- a `seq` b `seq` putStrLn ""- time "product of 1M different 3x3 matrices" (manymult (10^6))- putStrLn ""- time "product of 1M constant 1x1 matrices" (multb 1)- time "product of 1M constant 3x3 matrices" (multb 3)- --time "product of 1M constant 5x5 matrices" (multb 5)- time "product of 1M const. 10x10 matrices" (multb 10)- --time "product of 1M const. 15x15 matrices" (multb 15)- time "product of 1M const. 20x20 matrices" (multb 20)- --time "product of 1M const. 25x25 matrices" (multb 25)- putStrLn ""- time "product (1000 x 1000)<>(1000 x 1000)" (a<>a)- time "product (2000 x 2000)<>(2000 x 2000)" (b<>b)------------------------------------eigBench = do- let m = reshape 1000 (randomVector 777 Uniform (1000*1000))- s = m + trans m- m `seq` s `seq` putStrLn ""- time "eigenvalues symmetric 1000x1000" (eigenvaluesSH' m)- time "eigenvectors symmetric 1000x1000" (snd $ eigSH' m)- time "eigenvalues general 1000x1000" (eigenvalues m)- time "eigenvectors general 1000x1000" (snd $ eig m)------------------------------------svdBench = do- let a = reshape 500 (randomVector 777 Uniform (3000*500))- b = reshape 1000 (randomVector 777 Uniform (1000*1000))- fv (_,_,v) = v@@>(0,0)- a `seq` b `seq` putStrLn ""- time "singular values 3000x500" (singularValues a)- time "thin svd 3000x500" (fv $ thinSVD a)- time "full svd 3000x500" (fv $ svd a)- time "singular values 1000x1000" (singularValues b)- time "full svd 1000x1000" (fv $ svd b)------------------------------------solveBenchN n = do- let x = uniformSample 777 (2*n) (replicate n (-1,1))- a = trans x <> x- b = asColumn $ randomVector 666 Uniform n- a `seq` b `seq` putStrLn ""- time ("svd solve " ++ show n) (linearSolveSVD a b)- time (" ls solve " ++ show n) (linearSolveLS a b)- time (" solve " ++ show n) (linearSolve a b)- time ("cholSolve " ++ show n) (cholSolve (chol a) b)--solveBench = do- solveBenchN 500- solveBenchN 1000- -- solveBenchN 1500------------------------------------cholBenchN n = do- let x = uniformSample 777 (2*n) (replicate n (-1,1))- a = trans x <> x- a `seq` putStrLn ""- time ("chol " ++ show n) (chol a)--cholBench = do- cholBenchN 1200- cholBenchN 600- cholBenchN 300--- cholBenchN 150--- cholBenchN 50
− lib/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -1,249 +0,0 @@-{-# LANGUAGE FlexibleContexts, UndecidableInstances, CPP, FlexibleInstances #-}-{-# OPTIONS_GHC -fno-warn-unused-imports #-}-------------------------------------------------------------------------------{- |-Module : Numeric.LinearAlgebra.Tests.Instances-Copyright : (c) Alberto Ruiz 2008-License : GPL-style--Maintainer : Alberto Ruiz (aruiz at um dot es)-Stability : provisional-Portability : portable--Arbitrary instances for vectors, matrices.---}--module Numeric.LinearAlgebra.Tests.Instances(- Sq(..), rSq,cSq,- Rot(..), rRot,cRot,- Her(..), rHer,cHer,- WC(..), rWC,cWC,- SqWC(..), rSqWC, cSqWC,- PosDef(..), rPosDef, cPosDef,- Consistent(..), rConsist, cConsist,- RM,CM, rM,cM,- FM,ZM, fM,zM-) where--import System.Random--import Numeric.LinearAlgebra-import Control.Monad(replicateM)-#include "quickCheckCompat.h"--#if MIN_VERSION_QuickCheck(2,0,0)-shrinkListElementwise :: (Arbitrary a) => [a] -> [[a]]-shrinkListElementwise [] = []-shrinkListElementwise (x:xs) = [ y:xs | y <- shrink x ]- ++ [ x:ys | ys <- shrinkListElementwise xs ]--shrinkPair :: (Arbitrary a, Arbitrary b) => (a,b) -> [(a,b)]-shrinkPair (a,b) = [ (a,x) | x <- shrink b ] ++ [ (x,b) | x <- shrink a ]-#endif--#if MIN_VERSION_QuickCheck(2,1,1)-#else-instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where- arbitrary = do- re <- arbitrary- im <- arbitrary- return (re :+ im)--#if MIN_VERSION_QuickCheck(2,0,0)- shrink (re :+ im) = - [ u :+ v | (u,v) <- shrinkPair (re,im) ]-#else- -- this has been moved to the 'Coarbitrary' class in QuickCheck 2- coarbitrary = undefined -#endif--#endif--chooseDim = sized $ \m -> choose (1,max 1 m)--instance (Field a, Arbitrary a) => Arbitrary (Vector a) where - arbitrary = do m <- chooseDim- l <- vector m- return $ fromList l--#if MIN_VERSION_QuickCheck(2,0,0)- -- shrink any one of the components- shrink = map fromList . shrinkListElementwise . toList--#else- coarbitrary = undefined-#endif--instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where - arbitrary = do- m <- chooseDim- n <- chooseDim- l <- vector (m*n)- return $ (m><n) l--#if MIN_VERSION_QuickCheck(2,0,0)- -- shrink any one of the components- shrink a = map (rows a >< cols a)- . shrinkListElementwise- . concat . toLists - $ a-#else- coarbitrary = undefined-#endif----- a square matrix-newtype (Sq a) = Sq (Matrix a) deriving Show-instance (Element a, Arbitrary a) => Arbitrary (Sq a) where- arbitrary = do- n <- chooseDim- l <- vector (n*n)- return $ Sq $ (n><n) l--#if MIN_VERSION_QuickCheck(2,0,0)- shrink (Sq a) = [ Sq b | b <- shrink a ]-#else- coarbitrary = undefined-#endif----- a unitary matrix-newtype (Rot a) = Rot (Matrix a) deriving Show-instance (Field a, Arbitrary a) => Arbitrary (Rot a) where- arbitrary = do- Sq m <- arbitrary- let (q,_) = qr m- return (Rot q)--#if MIN_VERSION_QuickCheck(2,0,0)-#else- coarbitrary = undefined-#endif----- a complex hermitian or real symmetric matrix-newtype (Her a) = Her (Matrix a) deriving Show-instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where- arbitrary = do- Sq m <- arbitrary- let m' = m/2- return $ Her (m' + ctrans m')--#if MIN_VERSION_QuickCheck(2,0,0)-#else- coarbitrary = undefined-#endif--class (Field a, Arbitrary a, Element (RealOf a), Random (RealOf a)) => ArbitraryField a-instance ArbitraryField Double-instance ArbitraryField (Complex Double)----- a well-conditioned general matrix (the singular values are between 1 and 100)-newtype (WC a) = WC (Matrix a) deriving Show-instance (ArbitraryField a) => Arbitrary (WC a) where- arbitrary = do- m <- arbitrary- let (u,_,v) = svd m- r = rows m- c = cols m- n = min r c- sv' <- replicateM n (choose (1,100))- let s = diagRect 0 (fromList sv') r c- return $ WC (u <> real s <> trans v)--#if MIN_VERSION_QuickCheck(2,0,0)-#else- coarbitrary = undefined-#endif----- a well-conditioned square matrix (the singular values are between 1 and 100)-newtype (SqWC a) = SqWC (Matrix a) deriving Show-instance (ArbitraryField a) => Arbitrary (SqWC a) where- arbitrary = do- Sq m <- arbitrary- let (u,_,v) = svd m- n = rows m- sv' <- replicateM n (choose (1,100))- let s = diag (fromList sv')- return $ SqWC (u <> real s <> trans v)--#if MIN_VERSION_QuickCheck(2,0,0)-#else- coarbitrary = undefined-#endif----- a positive definite square matrix (the eigenvalues are between 0 and 100)-newtype (PosDef a) = PosDef (Matrix a) deriving Show-instance (ArbitraryField a, Num (Vector a)) - => Arbitrary (PosDef a) where- arbitrary = do- Her m <- arbitrary- let (_,v) = eigSH m- n = rows m- l <- replicateM n (choose (0,100))- let s = diag (fromList l)- p = v <> real s <> ctrans v- return $ PosDef (0.5 * p + 0.5 * ctrans p)--#if MIN_VERSION_QuickCheck(2,0,0)-#else- coarbitrary = undefined-#endif----- a pair of matrices that can be multiplied-newtype (Consistent a) = Consistent (Matrix a, Matrix a) deriving Show-instance (Field a, Arbitrary a) => Arbitrary (Consistent a) where- arbitrary = do- n <- chooseDim- k <- chooseDim- m <- chooseDim- la <- vector (n*k)- lb <- vector (k*m)- return $ Consistent ((n><k) la, (k><m) lb)--#if MIN_VERSION_QuickCheck(2,0,0)- shrink (Consistent (x,y)) = [ Consistent (u,v) | (u,v) <- shrinkPair (x,y) ]-#else- coarbitrary = undefined-#endif----type RM = Matrix Double-type CM = Matrix (Complex Double)-type FM = Matrix Float-type ZM = Matrix (Complex Float)---rM m = m :: RM-cM m = m :: CM-fM m = m :: FM-zM m = m :: ZM---rHer (Her m) = m :: RM-cHer (Her m) = m :: CM--rRot (Rot m) = m :: RM-cRot (Rot m) = m :: CM--rSq (Sq m) = m :: RM-cSq (Sq m) = m :: CM--rWC (WC m) = m :: RM-cWC (WC m) = m :: CM--rSqWC (SqWC m) = m :: RM-cSqWC (SqWC m) = m :: CM--rPosDef (PosDef m) = m :: RM-cPosDef (PosDef m) = m :: CM--rConsist (Consistent (a,b)) = (a,b::RM)-cConsist (Consistent (a,b)) = (a,b::CM)-
− lib/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -1,272 +0,0 @@-{-# LANGUAGE CPP, FlexibleContexts #-}-{-# OPTIONS_GHC -fno-warn-unused-imports #-}-------------------------------------------------------------------------------{- |-Module : Numeric.LinearAlgebra.Tests.Properties-Copyright : (c) Alberto Ruiz 2008-License : GPL-style--Maintainer : Alberto Ruiz (aruiz at um dot es)-Stability : provisional-Portability : portable--Testing properties.---}--module Numeric.LinearAlgebra.Tests.Properties (- dist, (|~|), (~:), Aprox((:~)),- zeros, ones,- square,- unitary,- hermitian,- wellCond,- positiveDefinite,- upperTriang,- upperHessenberg,- luProp,- invProp,- pinvProp,- detProp,- nullspaceProp,- bugProp,- svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4,- svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7,- eigProp, eigSHProp, eigProp2, eigSHProp2,- qrProp, rqProp, rqProp1, rqProp2, rqProp3,- hessProp,- schurProp1, schurProp2,- cholProp, exactProp,- expmDiagProp,- multProp1, multProp2,- subProp,- linearSolveProp, linearSolveProp2-) where--import Numeric.LinearAlgebra --hiding (real,complex)-import Numeric.LinearAlgebra.LAPACK-import Debug.Trace-#include "quickCheckCompat.h"-----real x = real'' x---complex x = complex'' x--debug x = trace (show x) x---- relative error-dist :: (Normed c t, Num (c t)) => c t -> c t -> Double-dist a b = realToFrac r- where norm = pnorm Infinity- na = norm a- nb = norm b- nab = norm (a-b)- mx = max na nb- mn = min na nb- r = if mn < peps- then mx- else nab/mx--infixl 4 |~|-a |~| b = a :~10~: b---a |~| b = dist a b < 10^^(-10)--data Aprox a = (:~) a Int--- (~:) :: (Normed a, Num a) => Aprox a -> a -> Bool-a :~n~: b = dist a b < 10^^(-n)----------------------------------------------------------square m = rows m == cols m---- orthonormal columns-orthonormal m = ctrans m <> m |~| ident (cols m)--unitary m = square m && orthonormal m--hermitian m = square m && m |~| ctrans m--wellCond m = rcond m > 1/100--positiveDefinite m = minimum (toList e) > 0- where (e,_v) = eigSH m--upperTriang m = rows m == 1 || down == z- where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m))- z = constant 0 (dim down)--upperHessenberg m = rows m < 3 || down == z- where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m))- z = constant 0 (dim down)--zeros (r,c) = reshape c (constant 0 (r*c))--ones (r,c) = zeros (r,c) + 1---------------------------------------------------------luProp m = m |~| p <> l <> u && f (det p) |~| f s- where (l,u,p,s) = lu m- f x = fromList [x]--invProp m = m <> inv m |~| ident (rows m)--pinvProp m = m <> p <> m |~| m- && p <> m <> p |~| p- && hermitian (m<>p)- && hermitian (p<>m)- where p = pinv m--detProp m = s d1 |~| s d2- where d1 = det m- d2 = det' * det q- det' = product $ toList $ takeDiag r- (q,r) = qr m- s x = fromList [x]--nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c)- && orthonormal (fromColumns nl))- where nl = nullspacePrec 1 m- n = fromColumns nl- r = rows m- c = cols m - rank m------------------------------------------------------------------------ testcase for nonempty fpu stack--- uncommenting unitary' signature eliminates the problem-bugProp m = m |~| u <> real d <> trans v && unitary' u && unitary' v- where (u,d,v) = fullSVD m- -- unitary' :: (Num (Vector t), Field t) => Matrix t -> Bool- unitary' a = unitary a------------------------------------------------------------------------ fullSVD-svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v- where (u,d,v) = fullSVD m--svdProp1a svdfun m = m |~| u <> real d <> trans v && unitary u && unitary v where- (u,s,v) = svdfun m- d = diagRect 0 s (rows m) (cols m)--svdProp1b svdfun m = unitary u && unitary v where- (u,_,v) = svdfun m---- thinSVD-svdProp2 thinSVDfun m = m |~| u <> diag (real s) <> trans v && orthonormal u && orthonormal v && dim s == min (rows m) (cols m)- where (u,s,v) = thinSVDfun m---- compactSVD-svdProp3 m = (m |~| u <> real (diag s) <> trans v- && orthonormal u && orthonormal v)- where (u,s,v) = compactSVD m--svdProp4 m' = m |~| u <> real (diag s) <> trans v- && orthonormal u && orthonormal v- && (dim s == r || r == 0 && dim s == 1)- where (u,s,v) = compactSVD m- m = fromBlocks [[m'],[m']]- r = rank m'--svdProp5a m = all (s1|~|) [s2,s3,s4,s5,s6] where- s1 = svR m- s2 = svRd m- (_,s3,_) = svdR m- (_,s4,_) = svdRd m- (_,s5,_) = thinSVDR m- (_,s6,_) = thinSVDRd m--svdProp5b m = all (s1|~|) [s2,s3,s4,s5,s6] where- s1 = svC m- s2 = svCd m- (_,s3,_) = svdC m- (_,s4,_) = svdCd m- (_,s5,_) = thinSVDC m- (_,s6,_) = thinSVDCd m--svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'- where (u,s,v) = svdR m- (s',v') = rightSVR m- (u',s'') = leftSVR m--svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'- where (u,s,v) = svdC m- (s',v') = rightSVC m- (u',s'') = leftSVC m--svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s'''- where (u,s,v) = svd m- (s',v') = rightSV m- (u',_s'') = leftSV m- s''' = singularValues m----------------------------------------------------------------------eigProp m = complex m <> v |~| v <> diag s- where (s, v) = eig m--eigSHProp m = m <> v |~| v <> real (diag s)- && unitary v- && m |~| v <> real (diag s) <> ctrans v- where (s, v) = eigSH m--eigProp2 m = fst (eig m) |~| eigenvalues m--eigSHProp2 m = fst (eigSH m) |~| eigenvaluesSH m----------------------------------------------------------------------qrProp m = q <> r |~| m && unitary q && upperTriang r- where (q,r) = qr m--rqProp m = r <> q |~| m && unitary q && upperTriang' r- where (r,q) = rq m--rqProp1 m = r <> q |~| m- where (r,q) = rq m--rqProp2 m = unitary q- where (r,q) = rq m--rqProp3 m = upperTriang' r- where (r,q) = rq m--upperTriang' r = upptr (rows r) (cols r) * r |~| r- where upptr f c = buildMatrix f c $ \(r',c') -> if r'-t > c' then 0 else 1- where t = f-c--hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h- where (p,h) = hess m--schurProp1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s- where (u,s) = schur m--schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme- where (u,s) = schur m--cholProp m = m |~| ctrans c <> c && upperTriang c- where c = chol m--exactProp m = chol m == chol (m+0)--expmDiagProp m = expm (logm m) :~ 7 ~: complex m- where logm = matFunc log---- reference multiply-mulH a b = fromLists [[ doth ai bj | bj <- toColumns b] | ai <- toRows a ]- where doth u v = sum $ zipWith (*) (toList u) (toList v)--multProp1 p (a,b) = (a <> b) :~p~: (mulH a b)--multProp2 p (a,b) = (ctrans (a <> b)) :~p~: (ctrans b <> ctrans a)--linearSolveProp f m = f m m |~| ident (rows m)--linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b)- where q = min (rows a) (cols a)- b = a <> x- wc = rank a == q--subProp m = m == (trans . fromColumns . toRows) m-
− lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h
@@ -1,33 +0,0 @@-#ifndef MIN_VERSION_QuickCheck-#define MIN_VERSION_QuickCheck(A,B,C) 1-#endif--#if MIN_VERSION_QuickCheck(2,0,0)-import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector- ,sized,classify,Testable,Property-- ,quickCheckWith,maxSize,stdArgs,shrink)--#else-import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector- ,sized,classify,Testable,Property-- ,check,configSize,defaultConfig,trivial)-#endif----#if MIN_VERSION_QuickCheck(2,0,0)-trivial :: Testable a => Bool -> a -> Property-trivial = (`classify` "trivial")-#else-#endif----- define qCheck, which used to be in Tests.hs-#if MIN_VERSION_QuickCheck(2,0,0)-qCheck n = quickCheckWith stdArgs {maxSize = n}-#else-qCheck n = check defaultConfig {configSize = const n}-#endif-