diff --git a/CHANGES b/CHANGES deleted file mode 100644 --- a/CHANGES +++ /dev/null @@ -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. - diff --git a/CHANGES.md b/CHANGES.md new file mode 100644 --- /dev/null +++ b/CHANGES.md @@ -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. + diff --git a/INSTALL b/INSTALL deleted file mode 100644 --- a/INSTALL +++ /dev/null @@ -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. - diff --git a/INSTALL.md b/INSTALL.md new file mode 100644 --- /dev/null +++ b/INSTALL.md @@ -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 + + diff --git a/Setup.lhs b/Setup.lhs --- a/Setup.lhs +++ b/Setup.lhs @@ -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() diff --git a/THANKS b/THANKS deleted file mode 100644 --- a/THANKS +++ /dev/null @@ -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. - diff --git a/THANKS.md b/THANKS.md new file mode 100644 --- /dev/null +++ b/THANKS.md @@ -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. + diff --git a/examples/bool.hs b/examples/bool.hs --- a/examples/bool.hs +++ b/examples/bool.hs @@ -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 diff --git a/examples/parallel.hs b/examples/parallel.hs --- a/examples/parallel.hs +++ b/examples/parallel.hs @@ -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 + diff --git a/examples/tests.hs b/examples/tests.hs deleted file mode 100644 --- a/examples/tests.hs +++ /dev/null @@ -1,3 +0,0 @@ -import Numeric.LinearAlgebra.Tests - -main = runTests 20 diff --git a/hmatrix.cabal b/hmatrix.cabal --- a/hmatrix.cabal +++ b/hmatrix.cabal @@ -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 diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs --- a/lib/Data/Packed/Internal/Matrix.hs +++ b/lib/Data/Packed/Internal/Matrix.hs @@ -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) diff --git a/lib/Data/Packed/Internal/Signatures.hs b/lib/Data/Packed/Internal/Signatures.hs --- a/lib/Data/Packed/Internal/Signatures.hs +++ b/lib/Data/Packed/Internal/Signatures.hs @@ -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 -- diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs --- a/lib/Data/Packed/Internal/Vector.hs +++ b/lib/Data/Packed/Internal/Vector.hs @@ -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 diff --git a/lib/Numeric/Container.hs b/lib/Numeric/Container.hs --- a/lib/Numeric/Container.hs +++ b/lib/Numeric/Container.hs @@ -59,11 +59,7 @@ readMatrix, fscanfVector, fprintfVector, freadVector, fwriteVector, -- * Experimental - build', konst', - -- * Deprecated - (.*),(*/),(<|>),(<->), - vectorMax,vectorMin, - vectorMaxIndex, vectorMinIndex + build', konst' ) where import Data.Packed diff --git a/lib/Numeric/ContainerBoot.hs b/lib/Numeric/ContainerBoot.hs --- a/lib/Numeric/ContainerBoot.hs +++ b/lib/Numeric/ContainerBoot.hs @@ -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 diff --git a/lib/Numeric/GSL/Differentiation.hs b/lib/Numeric/GSL/Differentiation.hs --- a/lib/Numeric/GSL/Differentiation.hs +++ b/lib/Numeric/GSL/Differentiation.hs @@ -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) diff --git a/lib/Numeric/GSL/Fitting.hs b/lib/Numeric/GSL/Fitting.hs --- a/lib/Numeric/GSL/Fitting.hs +++ b/lib/Numeric/GSL/Fitting.hs @@ -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) ------------------------------------------------------------------------- diff --git a/lib/Numeric/GSL/Fourier.hs b/lib/Numeric/GSL/Fourier.hs --- a/lib/Numeric/GSL/Fourier.hs +++ b/lib/Numeric/GSL/Fourier.hs @@ -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 diff --git a/lib/Numeric/GSL/Integration.hs b/lib/Numeric/GSL/Integration.hs --- a/lib/Numeric/GSL/Integration.hs +++ b/lib/Numeric/GSL/Integration.hs @@ -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) diff --git a/lib/Numeric/GSL/Internal.hs b/lib/Numeric/GSL/Internal.hs --- a/lib/Numeric/GSL/Internal.hs +++ b/lib/Numeric/GSL/Internal.hs @@ -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) diff --git a/lib/Numeric/GSL/Minimization.hs b/lib/Numeric/GSL/Minimization.hs --- a/lib/Numeric/GSL/Minimization.hs +++ b/lib/Numeric/GSL/Minimization.hs @@ -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) ------------------------------------------------------------------------ diff --git a/lib/Numeric/GSL/ODE.hs b/lib/Numeric/GSL/ODE.hs --- a/lib/Numeric/GSL/ODE.hs +++ b/lib/Numeric/GSL/ODE.hs @@ -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) ------------------------------------------------------------------------- diff --git a/lib/Numeric/GSL/Polynomials.hs b/lib/Numeric/GSL/Polynomials.hs --- a/lib/Numeric/GSL/Polynomials.hs +++ b/lib/Numeric/GSL/Polynomials.hs @@ -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 diff --git a/lib/Numeric/GSL/Root.hs b/lib/Numeric/GSL/Root.hs --- a/lib/Numeric/GSL/Root.hs +++ b/lib/Numeric/GSL/Root.hs @@ -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) ------------------------------------------------------------------------- diff --git a/lib/Numeric/GSL/Vector.hs b/lib/Numeric/GSL/Vector.hs --- a/lib/Numeric/GSL/Vector.hs +++ b/lib/Numeric/GSL/Vector.hs @@ -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 diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs --- a/lib/Numeric/LinearAlgebra/Algorithms.hs +++ b/lib/Numeric/LinearAlgebra/Algorithms.hs @@ -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 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 -------------------------------------------------------------- diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs --- a/lib/Numeric/LinearAlgebra/LAPACK.hs +++ b/lib/Numeric/LinearAlgebra/LAPACK.hs @@ -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) diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs deleted file mode 100644 --- a/lib/Numeric/LinearAlgebra/Tests.hs +++ /dev/null @@ -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 diff --git a/lib/Numeric/LinearAlgebra/Tests/Instances.hs b/lib/Numeric/LinearAlgebra/Tests/Instances.hs deleted file mode 100644 --- a/lib/Numeric/LinearAlgebra/Tests/Instances.hs +++ /dev/null @@ -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>< 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> 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> 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 - diff --git a/lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h b/lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h deleted file mode 100644 --- a/lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h +++ /dev/null @@ -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 -