diff --git a/CHANGELOG b/CHANGELOG new file mode 100644 --- /dev/null +++ b/CHANGELOG @@ -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. + diff --git a/CHANGES.md b/CHANGES.md deleted file mode 100644 --- a/CHANGES.md +++ /dev/null @@ -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. - diff --git a/INSTALL.md b/INSTALL.md --- a/INSTALL.md +++ b/INSTALL.md @@ -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). diff --git a/THANKS.md b/THANKS.md --- a/THANKS.md +++ b/THANKS.md @@ -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. + diff --git a/hmatrix.cabal b/hmatrix.cabal --- a/hmatrix.cabal +++ b/hmatrix.cabal @@ -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 diff --git a/lib/Data/Packed/Development.hs b/lib/Data/Packed/Development.hs --- a/lib/Data/Packed/Development.hs +++ b/lib/Data/Packed/Development.hs @@ -21,6 +21,7 @@ app1, app2, app3, app4, app5, app6, app7, app8, app9, app10, MatrixOrder(..), orderOf, cmat, fmat, + matrixFromVector, unsafeFromForeignPtr, unsafeToForeignPtr, check, (//), 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,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 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 @@ -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 diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c --- a/lib/Numeric/GSL/gsl-aux.c +++ b/lib/Numeric/GSL/gsl-aux.c @@ -29,6 +29,7 @@ #include #include #include +#include #include #include #include @@ -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 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 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)