hmatrix 0.15.0.1 → 0.15.2.0
raw patch · 10 files changed
+400/−228 lines, 10 filesdep ~base
Dependency ranges changed: base
Files
- CHANGELOG +188/−0
- CHANGES.md +0/−177
- INSTALL.md +1/−2
- THANKS.md +13/−0
- hmatrix.cabal +2/−2
- lib/Data/Packed/Development.hs +1/−0
- lib/Numeric/GSL/Integration.hs +50/−1
- lib/Numeric/GSL/Minimization.hs +34/−0
- lib/Numeric/GSL/gsl-aux.c +67/−0
- lib/Numeric/LinearAlgebra/Algorithms.hs +44/−46
+ CHANGELOG view
@@ -0,0 +1,188 @@+0.15.2.0+--------++ * general pinvTol and improved pinv++0.15.1.0+--------++ * One-dimensional minimization++ * Doubly-adaptive quadrature for difficult integrands++0.15.0.0+--------++ * Data.Packed.Foreign (additional FFI helpers)++ * NFData instance of Matrix++ * Unidimensional root finding++ * In Numeric.LinearAlgebra.Util:+ pairwise2D, rowOuters, null1, null1sym, size, unitary, mt, (¦), (?), (¿)++ * diagBlock++ * meanCov moved to Container++0.14.1.0+--------++ * In Numeric.LinearAlgebra.Util:+ convolution: corr, conv, corr2, conv2, separable, corrMin+ kronecker: vec, vech, dup, vtrans++0.14.0.0+--------++ * integration over infinite intervals++ * msadams and msbdf methods for ode++ * Numeric.LinearAlgebra.Util++ * (<\>) extended to multiple right-hand sides++ * orth++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.+
− CHANGES.md
@@ -1,177 +0,0 @@-0.15.0.0------------ Data.Packed.Foreign (additional FFI helpers)--- NFData instance of Matrix--- Unidimensional root finding--- In Numeric.LinearAlgebra.Util:- pairwise2D, rowOuters, null1, null1sym, size, unitary, mt, (¦), (?), (¿)--- diagBlock--- meanCov moved to Container--0.14.1.0------------ In Numeric.LinearAlgebra.Util:- convolution: corr, conv, corr2, conv2, separable, corrMin- kronecker: vec, vech, dup, vtrans--0.14.0.0------------ integration over infinite intervals--- msadams and msbdf methods for ode--- Numeric.LinearAlgebra.Util--- (<\>) extended to multiple right-hand sides--- orth--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.md view
@@ -1,4 +1,3 @@- # [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**.)@@ -24,7 +23,7 @@ GSL must be installed via MacPorts: - $ sudo port install gsl-devel +universal+ $ sudo port install gsl +universal $ cabal install hmatrix (Contributed by Heinrich Apfelmus and Torsten Kemps-Benedix).
THANKS.md view
@@ -128,3 +128,16 @@ - Greg Horn fixed the bus error on ghci 64-bit. +- Kristof Bastiaensen added bindings for one-dimensional minimization.++- Matthew Peddie added bindings for gsl_integrate_cquad doubly-adaptive quadrature+ for difficult integrands.++- Ben Gamari exposed matrixFromVector for Development.++- greg94301 reported tolerance issues in the tests.++- Clemens Lang updated the MacPort installation instructions.++- Henning Thielemann reported the pinv inefficient implementation.+
hmatrix.cabal view
@@ -1,5 +1,5 @@ Name: hmatrix-Version: 0.15.0.1+Version: 0.15.2.0 License: GPL License-file: LICENSE Author: Alberto Ruiz@@ -27,7 +27,7 @@ build-type: Custom -extra-source-files: Config.hs THANKS.md INSTALL.md CHANGES.md+extra-source-files: Config.hs THANKS.md INSTALL.md CHANGELOG extra-source-files: examples/deriv.hs examples/integrate.hs
lib/Data/Packed/Development.hs view
@@ -21,6 +21,7 @@ app1, app2, app3, app4, app5, app6, app7, app8, app9, app10, MatrixOrder(..), orderOf, cmat, fmat,+ matrixFromVector, unsafeFromForeignPtr, unsafeToForeignPtr, check, (//),
lib/Numeric/GSL/Integration.hs view
@@ -20,7 +20,8 @@ integrateQAGS, integrateQAGI, integrateQAGIU,- integrateQAGIL+ integrateQAGIL,+ integrateCQUAD ) where import Foreign.C.Types@@ -196,3 +197,51 @@ c_integrate_qagil :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt ++--------------------------------------------------------------------+{- | Numerical integration using /gsl_integration_cquad/ (quadrature+for general integrands). From the GSL manual:++@CQUAD is a new doubly-adaptive general-purpose quadrature routine+which can handle most types of singularities, non-numerical function+values such as Inf or NaN, as well as some divergent integrals. It+generally requires more function evaluations than the integration+routines in QUADPACK, yet fails less often for difficult integrands.@++For example:++@\> let quad = integrateCQUAD 1E-12 1000 +\> let f a x = exp(-a * x^2)+\> quad (f 0.5) 2 5+(5.7025405463957006e-2,9.678874441303705e-16,95)@++Unlike other quadrature methods, integrateCQUAD also returns the+number of function evaluations required.++-}++integrateCQUAD :: Double -- ^ precision (e.g. 1E-9)+ -> Int -- ^ size of auxiliary workspace (e.g. 1000)+ -> (Double -> Double) -- ^ function to be integrated on the interval (a, b)+ -> Double -- ^ a+ -> Double -- ^ b+ -> (Double, Double, Int) -- ^ result of the integration, error and number of function evaluations performed+integrateCQUAD prec n f a b = unsafePerformIO $ do+ r <- malloc+ e <- malloc+ neval <- malloc+ fp <- mkfun (\x _ -> f x)+ c_integrate_cquad fp a b prec (fromIntegral n) r e neval // check "integrate_cquad"+ vr <- peek r+ ve <- peek e+ vneval <- peek neval+ let result = (vr,ve,vneval)+ free r+ free e+ free neval+ freeHaskellFunPtr fp+ return result++foreign import ccall safe "gsl-aux.h integrate_cquad" + c_integrate_cquad :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> Double -> CInt+ -> Ptr Double -> Ptr Double -> Ptr Int -> IO CInt
lib/Numeric/GSL/Minimization.hs view
@@ -53,6 +53,7 @@ module Numeric.GSL.Minimization ( minimize, minimizeV, MinimizeMethod(..), minimizeD, minimizeVD, MinimizeMethodD(..),+ uniMinimize, UniMinimizeMethod(..), minimizeNMSimplex, minimizeConjugateGradient,@@ -80,6 +81,39 @@ minimizeVectorBFGS2 step tol eps maxit f g xi = minimizeD VectorBFGS2 eps maxit step tol f g xi -------------------------------------------------------------------------++data UniMinimizeMethod = GoldenSection+ | BrentMini+ | QuadGolden+ deriving (Enum, Eq, Show, Bounded)++-- | Onedimensional minimization.++uniMinimize :: UniMinimizeMethod -- ^ The method used.+ -> Double -- ^ desired precision of the solution+ -> Int -- ^ maximum number of iterations allowed+ -> (Double -> Double) -- ^ function to minimize+ -> Double -- ^ guess for the location of the minimum+ -> Double -- ^ lower bound of search interval+ -> Double -- ^ upper bound of search interval+ -> (Double, Matrix Double) -- ^ solution and optimization path++uniMinimize method epsrel maxit fun xmin xl xu = uniMinimizeGen (fi (fromEnum method)) fun xmin xl xu epsrel maxit++uniMinimizeGen m f xmin xl xu epsrel maxit = unsafePerformIO $ do+ fp <- mkDoublefun f+ rawpath <- createMIO maxit 4+ (c_uniMinize m fp epsrel (fi maxit) xmin xl xu)+ "uniMinimize"+ let it = round (rawpath @@> (maxit-1,0))+ path = takeRows it rawpath+ [sol] = toLists $ dropRows (it-1) path+ freeHaskellFunPtr fp+ return (sol !! 1, path)+++foreign import ccall safe "uniMinimize"+ c_uniMinize:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> Double -> TM data MinimizeMethod = NMSimplex | NMSimplex2
lib/Numeric/GSL/gsl-aux.c view
@@ -29,6 +29,7 @@ #include <gsl/gsl_poly.h> #include <gsl/gsl_multimin.h> #include <gsl/gsl_multiroots.h>+#include <gsl/gsl_min.h> #include <gsl/gsl_complex_math.h> #include <gsl/gsl_rng.h> #include <gsl/gsl_randist.h>@@ -802,7 +803,20 @@ OK } +int integrate_cquad(double f(double,void*), double a, double b, double prec,+ int w, double *result, double* error, int *neval) {+ DEBUGMSG("integrate_cquad");+ gsl_integration_cquad_workspace * wk = gsl_integration_cquad_workspace_alloc (w);+ gsl_function F;+ F.function = f;+ F.params = NULL;+ int res = gsl_integration_cquad (&F, a, b, 0, prec, wk, result, error, neval); + CHECK(res,res);+ gsl_integration_cquad_workspace_free (wk); + OK+} + int polySolve(KRVEC(a), CVEC(z)) { DEBUGMSG("polySolve"); REQUIRES(an>1,BAD_SIZE);@@ -891,6 +905,59 @@ free(p); return res; }++double only_f_aux_root(double x, void *pars);+int uniMinimize(int method, double f(double),+ double epsrel, int maxit, double min,+ double xl, double xu, RMAT(sol)) {+ REQUIRES(solr == maxit && solc == 4,BAD_SIZE);+ DEBUGMSG("minimize_only_f");+ gsl_function my_func;+ my_func.function = only_f_aux_root;+ my_func.params = f;+ size_t iter = 0;+ int status;+ const gsl_min_fminimizer_type *T;+ gsl_min_fminimizer *s;+ // Starting point+ switch(method) {+ case 0 : {T = gsl_min_fminimizer_goldensection; break; }+ case 1 : {T = gsl_min_fminimizer_brent; break; }+ case 2 : {T = gsl_min_fminimizer_quad_golden; break; }+ default: ERROR(BAD_CODE);+ }+ s = gsl_min_fminimizer_alloc (T);+ gsl_min_fminimizer_set (s, &my_func, min, xl, xu);+ do {+ double current_min, current_lo, current_hi;+ status = gsl_min_fminimizer_iterate (s);+ current_min = gsl_min_fminimizer_x_minimum (s);+ current_lo = gsl_min_fminimizer_x_lower (s);+ current_hi = gsl_min_fminimizer_x_upper (s);+ solp[iter*solc] = iter + 1;+ solp[iter*solc+1] = current_min;+ solp[iter*solc+2] = current_lo;+ solp[iter*solc+3] = current_hi;+ iter++;+ if (status) /* check if solver is stuck */+ break;+ + status =+ gsl_min_test_interval (current_lo, current_hi, 0, epsrel);+ }+ while (status == GSL_CONTINUE && iter < maxit);+ int i;+ for (i=iter; i<solr; i++) {+ solp[i*solc+0] = iter;+ solp[i*solc+1]=0.;+ solp[i*solc+2]=0.;+ solp[i*solc+3]=0.;+ }+ gsl_min_fminimizer_free(s);+ OK+}++ // this version returns info about intermediate steps int minimize(int method, double f(int, double*), double tolsize, int maxit,
lib/Numeric/LinearAlgebra/Algorithms.hs view
@@ -31,7 +31,7 @@ cholSolve, linearSolveLS, linearSolveSVD,- inv, pinv,+ inv, pinv, pinvTol, det, invlndet, rank, rcond, -- * Matrix factorizations@@ -73,7 +73,6 @@ -- * Util haussholder, unpackQR, unpackHess,- pinvTol, ranksv ) where @@ -95,7 +94,9 @@ Container Vector t, Container Matrix t, Normed Matrix t,- Normed Vector t) => Field t where+ Normed Vector t,+ Floating t,+ RealOf t ~ Double) => Field t where svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t) thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t) sv' :: Matrix t -> Vector Double@@ -330,9 +331,9 @@ -- | Joint computation of inverse and logarithm of determinant of a square matrix.-invlndet :: (Floating t, Field t)+invlndet :: Field t => Matrix t- -> (Matrix t, (t, t)) -- ^ (inverse, (log abs det, sign or phase of det)) + -> (Matrix t, (t, t)) -- ^ (inverse, (log abs det, sign or phase of det)) invlndet m | square m = (im,(ladm,sdm)) | otherwise = error $ "invlndet of nonsquare "++ shSize m ++ " matrix" where@@ -363,10 +364,43 @@ inv m | square m = m `linearSolve` ident (rows m) | otherwise = error $ "inv of nonsquare "++ shSize m ++ " matrix" --- | Pseudoinverse of a general matrix.++-- | Pseudoinverse of a general matrix with default tolerance ('pinvTol' 1, similar to GNU-Octave). pinv :: Field t => Matrix t -> Matrix t-pinv m = linearSolveSVD m (ident (rows m))+pinv = pinvTol 1 +{- | @pinvTol r@ computes the pseudoinverse of a matrix with tolerance @tol=r*g*eps*(max rows cols)@, where g is the greatest singular value.++@\> let m = 'fromLists' [[1,0, 0]+ ,[0,1, 0]+ ,[0,0,1e-10]]+\ --+\> 'pinv' m+1. 0. 0.+0. 1. 0.+0. 0. 10000000000.+\ --+\> pinvTol 1E8 m+1. 0. 0.+0. 1. 0.+0. 0. 1.@++-}++pinvTol :: Field t => Double -> Matrix t -> Matrix t+pinvTol t m = conj v' `mXm` diag s' `mXm` ctrans u' where+ (u,s,v) = thinSVD m+ sl@(g:_) = toList s+ s' = real . fromList . map rec $ sl+ rec x = if x <= g*tol then x else 1/x+ tol = (fromIntegral (max r c) * g * t * eps)+ r = rows m+ c = cols m+ d = dim s+ u' = takeColumns d u+ v' = takeColumns d v++ -- | Numeric rank of a matrix from the SVD decomposition. rankSVD :: Element t => Double -- ^ numeric zero (e.g. 1*'eps')@@ -439,39 +473,6 @@ ------------------------------------------------------------------------ -{- Pseudoinverse of a real matrix with the desired tolerance, expressed as a-multiplicative factor of the default tolerance used by GNU-Octave (see 'pinv').--@\> let m = 'fromLists' [[1,0, 0]- ,[0,1, 0]- ,[0,0,1e-10]]-\ ---\> 'pinv' m -1. 0. 0.-0. 1. 0.-0. 0. 10000000000.-\ ---\> pinvTol 1E8 m-1. 0. 0.-0. 1. 0.-0. 0. 1.@---}---pinvTol :: Double -> Matrix Double -> Matrix Double-pinvTol t m = v' `mXm` diag s' `mXm` trans u' where- (u,s,v) = thinSVDRd m- sl@(g:_) = toList s- s' = fromList . map rec $ sl- rec x = if x < g*tol then 1 else 1/x- tol = (fromIntegral (max r c) * g * t * eps)- r = rows m- c = cols m- d = dim s- u' = takeColumns d u- v' = takeColumns d v------------------------------------------------------------------------ -- many thanks, quickcheck! haussholder :: (Field a) => a -> Vector a -> Matrix a@@ -545,7 +546,7 @@ matFunc :: (Complex Double -> Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) matFunc f m = case diagonalize m of Just (l,v) -> v `mXm` diag (mapVector f l) `mXm` inv v- Nothing -> error "Sorry, matFunc requires a diagonalizable matrix" + Nothing -> error "Sorry, matFunc requires a diagonalizable matrix" -------------------------------------------------------------- @@ -556,6 +557,7 @@ c = fact (p+q) * fact (p+q+1) fact n = product [1..n] +epslist :: [(Int,Double)] epslist = [ (fromIntegral k, golubeps k k) | k <- [1..]] geps delta = head [ k | (k,g) <- epslist, g<delta]@@ -567,11 +569,7 @@ 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 :: Field 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)