diff options
Diffstat (limited to 'packages/base')
-rw-r--r-- | packages/base/CHANGELOG | 13 | ||||
-rw-r--r-- | packages/base/THANKS.md | 30 | ||||
-rw-r--r-- | packages/base/hmatrix.cabal | 20 | ||||
-rw-r--r-- | packages/base/src/C/vector-aux.c | 114 | ||||
-rw-r--r-- | packages/base/src/Data/Packed/Internal/Matrix.hs | 3 | ||||
-rw-r--r-- | packages/base/src/Data/Packed/Matrix.hs | 6 | ||||
-rw-r--r-- | packages/base/src/Data/Packed/Numeric.hs | 24 | ||||
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/Algorithms.hs | 9 | ||||
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/Data.hs | 4 | ||||
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/HMatrix.hs | 34 | ||||
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/Static.hs | 10 | ||||
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs | 2 | ||||
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/Util.hs | 44 | ||||
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs | 3 |
14 files changed, 271 insertions, 45 deletions
diff --git a/packages/base/CHANGELOG b/packages/base/CHANGELOG index 35ccdbc..c137285 100644 --- a/packages/base/CHANGELOG +++ b/packages/base/CHANGELOG | |||
@@ -1,3 +1,10 @@ | |||
1 | 0.16.1.0 | ||
2 | -------- | ||
3 | |||
4 | * Added (|||) and (===) for "besides" and "above" | ||
5 | |||
6 | * rowOuters | ||
7 | |||
1 | 0.16.0.0 | 8 | 0.16.0.0 |
2 | -------- | 9 | -------- |
3 | 10 | ||
@@ -11,9 +18,11 @@ | |||
11 | Numeric.LinearAlgebra.Data | 18 | Numeric.LinearAlgebra.Data |
12 | Numeric.LinearAlgebra.Devel | 19 | Numeric.LinearAlgebra.Devel |
13 | 20 | ||
14 | The documentation is now hidden for Data.Packed.*, Numeric.Container, and | 21 | For normal usage we only need to import Numeric.LinearAlgebra.HMatrix. |
22 | |||
23 | (The documentation is now hidden for Data.Packed.*, Numeric.Container, and | ||
15 | the other Numeric.LinearAlgebra.* modules, but they continue to be exposed | 24 | the other Numeric.LinearAlgebra.* modules, but they continue to be exposed |
16 | for backwards compatibility. | 25 | for backwards compatibility.) |
17 | 26 | ||
18 | * Added support for empty arrays, extending automatic conformability | 27 | * Added support for empty arrays, extending automatic conformability |
19 | (very useful for construction of block matrices). | 28 | (very useful for construction of block matrices). |
diff --git a/packages/base/THANKS.md b/packages/base/THANKS.md index b1417a6..fdbbe14 100644 --- a/packages/base/THANKS.md +++ b/packages/base/THANKS.md | |||
@@ -103,7 +103,7 @@ module reorganization, monadic mapVectorM, and many other improvements. | |||
103 | deprecation warnings, used more explicit imports, and updated to ghc-7.4. | 103 | deprecation warnings, used more explicit imports, and updated to ghc-7.4. |
104 | 104 | ||
105 | - Tom Nielsen discovered a problem in Config.hs, exposed by link problems | 105 | - Tom Nielsen discovered a problem in Config.hs, exposed by link problems |
106 | in Ubuntu 11.10 beta. | 106 | in Ubuntu 11.10 beta, and fixed the link options on freebsd. |
107 | 107 | ||
108 | - Daniel Fischer reported some Haddock markup errors. | 108 | - Daniel Fischer reported some Haddock markup errors. |
109 | 109 | ||
@@ -120,7 +120,8 @@ module reorganization, monadic mapVectorM, and many other improvements. | |||
120 | 120 | ||
121 | - Takano Akio fixed off-by-one errors in gsl-aux.c producing segfaults. | 121 | - Takano Akio fixed off-by-one errors in gsl-aux.c producing segfaults. |
122 | 122 | ||
123 | - Alex Lang implemented uniRoot and uniRootJ for one-dimensional root-finding. | 123 | - Alex Lang implemented uniRoot and uniRootJ for one-dimensional root-finding, and |
124 | fixed asRow and asColumn for empty vectors. | ||
124 | 125 | ||
125 | - Mike Ledger contributed alternative FFI helpers for matrix interoperation with C | 126 | - Mike Ledger contributed alternative FFI helpers for matrix interoperation with C |
126 | 127 | ||
@@ -139,7 +140,8 @@ module reorganization, monadic mapVectorM, and many other improvements. | |||
139 | 140 | ||
140 | - Clemens Lang updated the MacPort installation instructions. | 141 | - Clemens Lang updated the MacPort installation instructions. |
141 | 142 | ||
142 | - Henning Thielemann reported the pinv inefficient implementation. | 143 | - Henning Thielemann reported the pinv inefficient implementation and the need of |
144 | pkgconfig-depends. | ||
143 | 145 | ||
144 | - bdoering reported the problem of zero absolute tolerance in the integration functions. | 146 | - bdoering reported the problem of zero absolute tolerance in the integration functions. |
145 | 147 | ||
@@ -155,3 +157,25 @@ module reorganization, monadic mapVectorM, and many other improvements. | |||
155 | 157 | ||
156 | - Samium Gromoff reported a build failure caused by a size_t - int mismatch. | 158 | - Samium Gromoff reported a build failure caused by a size_t - int mismatch. |
157 | 159 | ||
160 | - Denis Laxalde separated the gsl tests from the base ones. | ||
161 | |||
162 | - Dominic Steinitz (idontgetoutmuch) reported a bug in the static diagonal creation functions. | ||
163 | |||
164 | - Dylan Thurston reported an error in the glpk documentation and ambiguity in | ||
165 | the description of linearSolve. | ||
166 | |||
167 | - Adrian Victor Crisciu developed an installation method for platforms which | ||
168 | don't provide shared lapack libraries. | ||
169 | |||
170 | - Ian Ross reported the max/minIndex bug. | ||
171 | |||
172 | - Niklas Hambüchen improved the documentation. | ||
173 | |||
174 | - "erdeszt" optimized "conv" using a direct vector reverse. | ||
175 | |||
176 | - John Shahbazian added support for openBLAS. | ||
177 | |||
178 | - "yongqli" reported the bug in randomVector (rand() is not thread-safe and drand48_r() is not portable). | ||
179 | |||
180 | - Kiwamu Ishikura improved randomVector for OSX | ||
181 | |||
diff --git a/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal index ba719a1..8ba3e06 100644 --- a/packages/base/hmatrix.cabal +++ b/packages/base/hmatrix.cabal | |||
@@ -1,5 +1,5 @@ | |||
1 | Name: hmatrix | 1 | Name: hmatrix |
2 | Version: 0.16.0.3 | 2 | Version: 0.16.1.3 |
3 | License: BSD3 | 3 | License: BSD3 |
4 | License-file: LICENSE | 4 | License-file: LICENSE |
5 | Author: Alberto Ruiz | 5 | Author: Alberto Ruiz |
@@ -14,6 +14,12 @@ Description: Linear algebra based on BLAS and LAPACK. | |||
14 | ["Numeric.LinearAlgebra.HMatrix"] Starting point and recommended import module for most applications. | 14 | ["Numeric.LinearAlgebra.HMatrix"] Starting point and recommended import module for most applications. |
15 | . | 15 | . |
16 | ["Numeric.LinearAlgebra.Static"] Experimental alternative interface. | 16 | ["Numeric.LinearAlgebra.Static"] Experimental alternative interface. |
17 | . | ||
18 | ["Numeric.LinearAlgebra.Devel"] Tools for extending the library. | ||
19 | . | ||
20 | (Other modules are exposed with hidden documentation for backwards compatibility.) | ||
21 | . | ||
22 | Code examples: <http://dis.um.es/~alberto/hmatrix/hmatrix.html> | ||
17 | 23 | ||
18 | Category: Math | 24 | Category: Math |
19 | tested-with: GHC==7.8 | 25 | tested-with: GHC==7.8 |
@@ -26,6 +32,10 @@ extra-source-files: THANKS.md CHANGELOG | |||
26 | 32 | ||
27 | extra-source-files: src/C/lapack-aux.h | 33 | extra-source-files: src/C/lapack-aux.h |
28 | 34 | ||
35 | flag openblas | ||
36 | description: Link with OpenBLAS (https://github.com/xianyi/OpenBLAS) optimized libraries. | ||
37 | default: False | ||
38 | |||
29 | library | 39 | library |
30 | 40 | ||
31 | Build-Depends: base >= 4 && < 5, | 41 | Build-Depends: base >= 4 && < 5, |
@@ -94,7 +104,11 @@ library | |||
94 | 104 | ||
95 | cpp-options: -DBINARY | 105 | cpp-options: -DBINARY |
96 | 106 | ||
97 | extra-libraries: blas lapack | 107 | if flag(openblas) |
108 | extra-lib-dirs: /usr/lib/openblas/lib | ||
109 | extra-libraries: openblas | ||
110 | else | ||
111 | extra-libraries: blas lapack | ||
98 | 112 | ||
99 | if os(OSX) | 113 | if os(OSX) |
100 | extra-lib-dirs: /opt/local/lib/ | 114 | extra-lib-dirs: /opt/local/lib/ |
@@ -108,7 +122,7 @@ library | |||
108 | if os(freebsd) | 122 | if os(freebsd) |
109 | extra-lib-dirs: /usr/local/lib | 123 | extra-lib-dirs: /usr/local/lib |
110 | include-dirs: /usr/local/include | 124 | include-dirs: /usr/local/include |
111 | extra-libraries: blas lapack | 125 | extra-libraries: blas lapack gfortran |
112 | 126 | ||
113 | if os(windows) | 127 | if os(windows) |
114 | extra-libraries: blas lapack | 128 | extra-libraries: blas lapack |
diff --git a/packages/base/src/C/vector-aux.c b/packages/base/src/C/vector-aux.c index 2f47c8f..dda47cb 100644 --- a/packages/base/src/C/vector-aux.c +++ b/packages/base/src/C/vector-aux.c | |||
@@ -172,7 +172,7 @@ double vector_min(KDVEC(x)) { | |||
172 | double vector_max_index(KDVEC(x)) { | 172 | double vector_max_index(KDVEC(x)) { |
173 | int k, r = 0; | 173 | int k, r = 0; |
174 | for (k = 1; k<xn; k++) { | 174 | for (k = 1; k<xn; k++) { |
175 | if(xp[k]>xp[0]) { | 175 | if(xp[k]>xp[r]) { |
176 | r = k; | 176 | r = k; |
177 | } | 177 | } |
178 | } | 178 | } |
@@ -182,7 +182,7 @@ double vector_max_index(KDVEC(x)) { | |||
182 | double vector_min_index(KDVEC(x)) { | 182 | double vector_min_index(KDVEC(x)) { |
183 | int k, r = 0; | 183 | int k, r = 0; |
184 | for (k = 1; k<xn; k++) { | 184 | for (k = 1; k<xn; k++) { |
185 | if(xp[k]<xp[0]) { | 185 | if(xp[k]<xp[r]) { |
186 | r = k; | 186 | r = k; |
187 | } | 187 | } |
188 | } | 188 | } |
@@ -237,7 +237,7 @@ float vector_min_f(KFVEC(x)) { | |||
237 | float vector_max_index_f(KFVEC(x)) { | 237 | float vector_max_index_f(KFVEC(x)) { |
238 | int k, r = 0; | 238 | int k, r = 0; |
239 | for (k = 1; k<xn; k++) { | 239 | for (k = 1; k<xn; k++) { |
240 | if(xp[k]>xp[0]) { | 240 | if(xp[k]>xp[r]) { |
241 | r = k; | 241 | r = k; |
242 | } | 242 | } |
243 | } | 243 | } |
@@ -247,7 +247,7 @@ float vector_max_index_f(KFVEC(x)) { | |||
247 | float vector_min_index_f(KFVEC(x)) { | 247 | float vector_min_index_f(KFVEC(x)) { |
248 | int k, r = 0; | 248 | int k, r = 0; |
249 | for (k = 1; k<xn; k++) { | 249 | for (k = 1; k<xn; k++) { |
250 | if(xp[k]<xp[0]) { | 250 | if(xp[k]<xp[r]) { |
251 | r = k; | 251 | r = k; |
252 | } | 252 | } |
253 | } | 253 | } |
@@ -700,17 +700,86 @@ int saveMatrix(char * file, char * format, KDMAT(a)){ | |||
700 | 700 | ||
701 | //////////////////////////////////////////////////////////////////////////////// | 701 | //////////////////////////////////////////////////////////////////////////////// |
702 | 702 | ||
703 | #ifdef __APPLE__ | ||
704 | |||
705 | #pragma message "randomVector is not thread-safe in OSX" | ||
706 | |||
707 | inline double urandom() { | ||
708 | const long max_random = 2147483647; // 2**31 - 1 | ||
709 | return (double)random() / (double)max_random; | ||
710 | } | ||
711 | |||
712 | double gaussrand(int *phase, double *pV1, double *pV2, double *pS) | ||
713 | { | ||
714 | double V1=*pV1, V2=*pV2, S=*pS; | ||
715 | double X; | ||
716 | |||
717 | if(*phase == 0) { | ||
718 | do { | ||
719 | double U1 = urandom(); | ||
720 | double U2 = urandom(); | ||
721 | |||
722 | V1 = 2 * U1 - 1; | ||
723 | V2 = 2 * U2 - 1; | ||
724 | S = V1 * V1 + V2 * V2; | ||
725 | } while(S >= 1 || S == 0); | ||
726 | |||
727 | X = V1 * sqrt(-2 * log(S) / S); | ||
728 | } else | ||
729 | X = V2 * sqrt(-2 * log(S) / S); | ||
730 | |||
731 | *phase = 1 - *phase; | ||
732 | *pV1=V1; *pV2=V2; *pS=S; | ||
733 | |||
734 | return X; | ||
735 | |||
736 | } | ||
737 | |||
738 | int random_vector(unsigned int seed, int code, DVEC(r)) { | ||
739 | int phase = 0; | ||
740 | double V1,V2,S; | ||
741 | |||
742 | srandom(seed); | ||
743 | |||
744 | int k; | ||
745 | switch (code) { | ||
746 | case 0: { // uniform | ||
747 | for (k=0; k<rn; k++) { | ||
748 | rp[k] = urandom(); | ||
749 | } | ||
750 | OK | ||
751 | } | ||
752 | case 1: { // gaussian | ||
753 | for (k=0; k<rn; k++) { | ||
754 | rp[k] = gaussrand(&phase,&V1,&V2,&S); | ||
755 | } | ||
756 | OK | ||
757 | } | ||
758 | |||
759 | default: ERROR(BAD_CODE); | ||
760 | } | ||
761 | } | ||
762 | |||
763 | #else | ||
764 | |||
765 | inline double urandom(struct random_data * buffer) { | ||
766 | int32_t res; | ||
767 | random_r(buffer,&res); | ||
768 | return (double)res/RAND_MAX; | ||
769 | } | ||
770 | |||
771 | |||
703 | // http://c-faq.com/lib/gaussian.html | 772 | // http://c-faq.com/lib/gaussian.html |
704 | double gaussrand() | 773 | double gaussrand(struct random_data *buffer, |
774 | int *phase, double *pV1, double *pV2, double *pS) | ||
705 | { | 775 | { |
706 | static double V1, V2, S; | 776 | double V1=*pV1, V2=*pV2, S=*pS; |
707 | static int phase = 0; | ||
708 | double X; | 777 | double X; |
709 | 778 | ||
710 | if(phase == 0) { | 779 | if(*phase == 0) { |
711 | do { | 780 | do { |
712 | double U1 = (double)rand() / RAND_MAX; | 781 | double U1 = urandom(buffer); |
713 | double U2 = (double)rand() / RAND_MAX; | 782 | double U2 = urandom(buffer); |
714 | 783 | ||
715 | V1 = 2 * U1 - 1; | 784 | V1 = 2 * U1 - 1; |
716 | V2 = 2 * U2 - 1; | 785 | V2 = 2 * U2 - 1; |
@@ -721,24 +790,37 @@ double gaussrand() | |||
721 | } else | 790 | } else |
722 | X = V2 * sqrt(-2 * log(S) / S); | 791 | X = V2 * sqrt(-2 * log(S) / S); |
723 | 792 | ||
724 | phase = 1 - phase; | 793 | *phase = 1 - *phase; |
794 | *pV1=V1; *pV2=V2; *pS=S; | ||
725 | 795 | ||
726 | return X; | 796 | return X; |
797 | |||
727 | } | 798 | } |
728 | 799 | ||
729 | int random_vector(int seed, int code, DVEC(r)) { | 800 | int random_vector(unsigned int seed, int code, DVEC(r)) { |
730 | srand(seed); | 801 | struct random_data buffer; |
802 | char random_state[128]; | ||
803 | memset(&buffer, 0, sizeof(struct random_data)); | ||
804 | memset(random_state, 0, sizeof(random_state)); | ||
805 | |||
806 | initstate_r(seed,random_state,sizeof(random_state),&buffer); | ||
807 | // setstate_r(random_state,&buffer); | ||
808 | // srandom_r(seed,&buffer); | ||
809 | |||
810 | int phase = 0; | ||
811 | double V1,V2,S; | ||
812 | |||
731 | int k; | 813 | int k; |
732 | switch (code) { | 814 | switch (code) { |
733 | case 0: { // uniform | 815 | case 0: { // uniform |
734 | for (k=0; k<rn; k++) { | 816 | for (k=0; k<rn; k++) { |
735 | rp[k] = (double)rand()/RAND_MAX; | 817 | rp[k] = urandom(&buffer); |
736 | } | 818 | } |
737 | OK | 819 | OK |
738 | } | 820 | } |
739 | case 1: { // gaussian | 821 | case 1: { // gaussian |
740 | for (k=0; k<rn; k++) { | 822 | for (k=0; k<rn; k++) { |
741 | rp[k] = gaussrand(); | 823 | rp[k] = gaussrand(&buffer,&phase,&V1,&V2,&S); |
742 | } | 824 | } |
743 | OK | 825 | OK |
744 | } | 826 | } |
@@ -747,6 +829,8 @@ int random_vector(int seed, int code, DVEC(r)) { | |||
747 | } | 829 | } |
748 | } | 830 | } |
749 | 831 | ||
832 | #endif | ||
833 | |||
750 | //////////////////////////////////////////////////////////////////////////////// | 834 | //////////////////////////////////////////////////////////////////////////////// |
751 | 835 | ||
752 | int smXv(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) { | 836 | int smXv(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) { |
diff --git a/packages/base/src/Data/Packed/Internal/Matrix.hs b/packages/base/src/Data/Packed/Internal/Matrix.hs index 91a9466..150b978 100644 --- a/packages/base/src/Data/Packed/Internal/Matrix.hs +++ b/packages/base/src/Data/Packed/Internal/Matrix.hs | |||
@@ -253,7 +253,8 @@ compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2 | |||
253 | operations for selected element types. | 253 | operations for selected element types. |
254 | It provides unoptimised defaults for any 'Storable' type, | 254 | It provides unoptimised defaults for any 'Storable' type, |
255 | so you can create instances simply as: | 255 | so you can create instances simply as: |
256 | @instance Element Foo@. | 256 | |
257 | >instance Element Foo | ||
257 | -} | 258 | -} |
258 | class (Storable a) => Element a where | 259 | class (Storable a) => Element a where |
259 | subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position | 260 | subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position |
diff --git a/packages/base/src/Data/Packed/Matrix.hs b/packages/base/src/Data/Packed/Matrix.hs index 6445ce4..70b9232 100644 --- a/packages/base/src/Data/Packed/Matrix.hs +++ b/packages/base/src/Data/Packed/Matrix.hs | |||
@@ -284,7 +284,7 @@ fromLists = fromRows . map fromList | |||
284 | -- [ 1.0, 2.0, 3.0, 4.0, 5.0 ] | 284 | -- [ 1.0, 2.0, 3.0, 4.0, 5.0 ] |
285 | -- | 285 | -- |
286 | asRow :: Storable a => Vector a -> Matrix a | 286 | asRow :: Storable a => Vector a -> Matrix a |
287 | asRow v = reshape (dim v) v | 287 | asRow = trans . asColumn |
288 | 288 | ||
289 | -- | creates a 1-column matrix from a vector | 289 | -- | creates a 1-column matrix from a vector |
290 | -- | 290 | -- |
@@ -297,7 +297,7 @@ asRow v = reshape (dim v) v | |||
297 | -- , 5.0 ] | 297 | -- , 5.0 ] |
298 | -- | 298 | -- |
299 | asColumn :: Storable a => Vector a -> Matrix a | 299 | asColumn :: Storable a => Vector a -> Matrix a |
300 | asColumn = trans . asRow | 300 | asColumn v = reshape 1 v |
301 | 301 | ||
302 | 302 | ||
303 | 303 | ||
@@ -476,7 +476,7 @@ mapMatrixWithIndexM g m = liftM (reshape c) . mapVectorWithIndexM (mk c g) . fla | |||
476 | 476 | ||
477 | {- | | 477 | {- | |
478 | 478 | ||
479 | >>> mapMatrixWithIndex (\\(i,j) v -> 100*v + 10*fromIntegral i + fromIntegral j) (ident 3:: Matrix Double) | 479 | >>> mapMatrixWithIndex (\(i,j) v -> 100*v + 10*fromIntegral i + fromIntegral j) (ident 3:: Matrix Double) |
480 | (3><3) | 480 | (3><3) |
481 | [ 100.0, 1.0, 2.0 | 481 | [ 100.0, 1.0, 2.0 |
482 | , 10.0, 111.0, 12.0 | 482 | , 10.0, 111.0, 12.0 |
diff --git a/packages/base/src/Data/Packed/Numeric.hs b/packages/base/src/Data/Packed/Numeric.hs index 7aa53f1..6027f43 100644 --- a/packages/base/src/Data/Packed/Numeric.hs +++ b/packages/base/src/Data/Packed/Numeric.hs | |||
@@ -160,7 +160,29 @@ instance Mul Vector Matrix Vector where | |||
160 | 160 | ||
161 | -------------------------------------------------------------------------------- | 161 | -------------------------------------------------------------------------------- |
162 | 162 | ||
163 | -- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD) | 163 | {- | Least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD) |
164 | |||
165 | @ | ||
166 | a = (3><2) | ||
167 | [ 1.0, 2.0 | ||
168 | , 2.0, 4.0 | ||
169 | , 2.0, -1.0 ] | ||
170 | @ | ||
171 | |||
172 | @ | ||
173 | v = vector [13.0,27.0,1.0] | ||
174 | @ | ||
175 | |||
176 | >>> let x = a <\> v | ||
177 | >>> x | ||
178 | fromList [3.0799999999999996,5.159999999999999] | ||
179 | |||
180 | >>> a #> x | ||
181 | fromList [13.399999999999999,26.799999999999997,1.0] | ||
182 | |||
183 | It also admits multiple right-hand sides stored as columns in a matrix. | ||
184 | |||
185 | -} | ||
164 | infixl 7 <\> | 186 | infixl 7 <\> |
165 | (<\>) :: (LSDiv c, Field t) => Matrix t -> c t -> c t | 187 | (<\>) :: (LSDiv c, Field t) => Matrix t -> c t -> c t |
166 | (<\>) = linSolve | 188 | (<\>) = linSolve |
diff --git a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs index 25700bc..02ac6a0 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs | |||
@@ -809,7 +809,7 @@ expGolub m = iterate msq f !! j | |||
809 | -------------------------------------------------------------- | 809 | -------------------------------------------------------------- |
810 | 810 | ||
811 | {- | Matrix square root. Currently it uses a simple iterative algorithm described in Wikipedia. | 811 | {- | Matrix square root. Currently it uses a simple iterative algorithm described in Wikipedia. |
812 | It only works with invertible matrices that have a real solution. For diagonalizable matrices you can try @matFunc sqrt@. | 812 | It only works with invertible matrices that have a real solution. |
813 | 813 | ||
814 | @m = (2><2) [4,9 | 814 | @m = (2><2) [4,9 |
815 | ,0,4] :: Matrix Double@ | 815 | ,0,4] :: Matrix Double@ |
@@ -819,6 +819,13 @@ It only works with invertible matrices that have a real solution. For diagonaliz | |||
819 | [ 2.0, 2.25 | 819 | [ 2.0, 2.25 |
820 | , 0.0, 2.0 ] | 820 | , 0.0, 2.0 ] |
821 | 821 | ||
822 | For diagonalizable matrices you can try 'matFunc' @sqrt@: | ||
823 | |||
824 | >>> matFunc sqrt ((2><2) [1,0,0,-1]) | ||
825 | (2><2) | ||
826 | [ 1.0 :+ 0.0, 0.0 :+ 0.0 | ||
827 | , 0.0 :+ 0.0, 0.0 :+ 1.0 ] | ||
828 | |||
822 | -} | 829 | -} |
823 | sqrtm :: Field t => Matrix t -> Matrix t | 830 | sqrtm :: Field t => Matrix t -> Matrix t |
824 | sqrtm = sqrtmInv | 831 | sqrtm = sqrtmInv |
diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs index b1a31fc..6dea407 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Data.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs | |||
@@ -40,7 +40,7 @@ module Numeric.LinearAlgebra.Data( | |||
40 | takeRows, dropRows, takeColumns, dropColumns, subMatrix, (?), (¿), fliprl, flipud, | 40 | takeRows, dropRows, takeColumns, dropColumns, subMatrix, (?), (¿), fliprl, flipud, |
41 | 41 | ||
42 | -- * Block matrix | 42 | -- * Block matrix |
43 | fromBlocks, (¦), (——), diagBlock, repmat, toBlocks, toBlocksEvery, | 43 | fromBlocks, (|||), (===), diagBlock, repmat, toBlocks, toBlocksEvery, |
44 | 44 | ||
45 | -- * Mapping functions | 45 | -- * Mapping functions |
46 | conj, cmap, step, cond, | 46 | conj, cmap, step, cond, |
@@ -66,7 +66,7 @@ module Numeric.LinearAlgebra.Data( | |||
66 | arctan2, | 66 | arctan2, |
67 | rows, cols, | 67 | rows, cols, |
68 | separable, | 68 | separable, |
69 | 69 | (¦),(——), | |
70 | module Data.Complex, | 70 | module Data.Complex, |
71 | 71 | ||
72 | Vector, Matrix, GMatrix, nRows, nCols | 72 | Vector, Matrix, GMatrix, nRows, nCols |
diff --git a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs index d2cae6c..677f9ee 100644 --- a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs +++ b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs | |||
@@ -134,7 +134,7 @@ module Numeric.LinearAlgebra.HMatrix ( | |||
134 | Seed, RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample, | 134 | Seed, RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample, |
135 | 135 | ||
136 | -- * Misc | 136 | -- * Misc |
137 | meanCov, peps, relativeError, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv, | 137 | meanCov, rowOuters, peps, relativeError, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv, |
138 | ℝ,ℂ,iC, | 138 | ℝ,ℂ,iC, |
139 | -- * Auxiliary classes | 139 | -- * Auxiliary classes |
140 | Element, Container, Product, Numeric, LSDiv, | 140 | Element, Container, Product, Numeric, LSDiv, |
@@ -194,7 +194,37 @@ mul :: Numeric t => Matrix t -> Matrix t -> Matrix t | |||
194 | mul = mXm | 194 | mul = mXm |
195 | 195 | ||
196 | 196 | ||
197 | -- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, returning Nothing for a singular system. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'. | 197 | {- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, returning Nothing for a singular system. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'. |
198 | |||
199 | @ | ||
200 | a = (2><2) | ||
201 | [ 1.0, 2.0 | ||
202 | , 3.0, 5.0 ] | ||
203 | @ | ||
204 | |||
205 | @ | ||
206 | b = (2><3) | ||
207 | [ 6.0, 1.0, 10.0 | ||
208 | , 15.0, 3.0, 26.0 ] | ||
209 | @ | ||
210 | |||
211 | >>> linearSolve a b | ||
212 | Just (2><3) | ||
213 | [ -1.4802973661668753e-15, 0.9999999999999997, 1.999999999999997 | ||
214 | , 3.000000000000001, 1.6653345369377348e-16, 4.000000000000002 ] | ||
215 | |||
216 | >>> let Just x = it | ||
217 | >>> disp 5 x | ||
218 | 2x3 | ||
219 | -0.00000 1.00000 2.00000 | ||
220 | 3.00000 0.00000 4.00000 | ||
221 | |||
222 | >>> a <> x | ||
223 | (2><3) | ||
224 | [ 6.0, 1.0, 10.0 | ||
225 | , 15.0, 3.0, 26.0 ] | ||
226 | |||
227 | -} | ||
198 | linearSolve m b = A.mbLinearSolve m b | 228 | linearSolve m b = A.mbLinearSolve m b |
199 | 229 | ||
200 | -- | return an orthonormal basis of the null space of a matrix. See also 'nullspaceSVD'. | 230 | -- | return an orthonormal basis of the null space of a matrix. See also 'nullspaceSVD'. |
diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs index 5749c40..037396d 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Static.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs | |||
@@ -512,7 +512,10 @@ crossC (extract -> x) (extract -> y) = mkC (LA.fromList [z1, z2, z3]) | |||
512 | -------------------------------------------------------------------------------- | 512 | -------------------------------------------------------------------------------- |
513 | 513 | ||
514 | diagRectR :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => ℝ -> R k -> L m n | 514 | diagRectR :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => ℝ -> R k -> L m n |
515 | diagRectR x v = r | 515 | diagRectR x v |
516 | | m' == 1 = mkL (LA.diagRect x ev m' n') | ||
517 | | m'*n' > 0 = r | ||
518 | | otherwise = matrix [] | ||
516 | where | 519 | where |
517 | r = mkL (asRow (vjoin [scalar x, ev, zeros])) | 520 | r = mkL (asRow (vjoin [scalar x, ev, zeros])) |
518 | ev = extract v | 521 | ev = extract v |
@@ -521,7 +524,10 @@ diagRectR x v = r | |||
521 | 524 | ||
522 | 525 | ||
523 | diagRectC :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => ℂ -> C k -> M m n | 526 | diagRectC :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => ℂ -> C k -> M m n |
524 | diagRectC x v = r | 527 | diagRectC x v |
528 | | m' == 1 = mkM (LA.diagRect x ev m' n') | ||
529 | | m'*n' > 0 = r | ||
530 | | otherwise = fromList [] | ||
525 | where | 531 | where |
526 | r = mkM (asRow (vjoin [scalar x, ev, zeros])) | 532 | r = mkM (asRow (vjoin [scalar x, ev, zeros])) |
527 | ev = extract v | 533 | ev = extract v |
diff --git a/packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs b/packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs index 339ef7d..ec02cf6 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs | |||
@@ -150,7 +150,7 @@ gmat st xs' | |||
150 | (xs,rest) = splitAt (m'*n') xs' | 150 | (xs,rest) = splitAt (m'*n') xs' |
151 | v = LA.fromList xs | 151 | v = LA.fromList xs |
152 | x = reshape n' v | 152 | x = reshape n' v |
153 | ok = rem (LA.size v) n' == 0 && LA.size x == (m',n') && null rest | 153 | ok = null rest && ((n' == 0 && dim v == 0) || n'> 0 && (rem (LA.size v) n' == 0) && LA.size x == (m',n')) |
154 | m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int | 154 | m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int |
155 | n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int | 155 | n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int |
156 | abort info = error $ st ++" "++show m' ++ " " ++ show n'++" can't be created from elements " ++ info | 156 | abort info = error $ st ++" "++show m' ++ " " ++ show n'++" can't be created from elements " ++ info |
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util.hs b/packages/base/src/Numeric/LinearAlgebra/Util.hs index 6bb9d15..043aa21 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Util.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Util.hs | |||
@@ -33,7 +33,7 @@ module Numeric.LinearAlgebra.Util( | |||
33 | diagl, | 33 | diagl, |
34 | row, | 34 | row, |
35 | col, | 35 | col, |
36 | (&), (¦), (——), (#), | 36 | (&), (¦), (|||), (——), (===), (#), |
37 | (?), (¿), | 37 | (?), (¿), |
38 | Indexable(..), size, | 38 | Indexable(..), size, |
39 | Numeric, | 39 | Numeric, |
@@ -157,25 +157,34 @@ a & b = vjoin [a,b] | |||
157 | 157 | ||
158 | {- | horizontal concatenation of real matrices | 158 | {- | horizontal concatenation of real matrices |
159 | 159 | ||
160 | (unicode 0x00a6, broken bar) | 160 | >>> ident 3 ||| konst 7 (3,4) |
161 | |||
162 | >>> ident 3 ¦ konst 7 (3,4) | ||
163 | (3><7) | 161 | (3><7) |
164 | [ 1.0, 0.0, 0.0, 7.0, 7.0, 7.0, 7.0 | 162 | [ 1.0, 0.0, 0.0, 7.0, 7.0, 7.0, 7.0 |
165 | , 0.0, 1.0, 0.0, 7.0, 7.0, 7.0, 7.0 | 163 | , 0.0, 1.0, 0.0, 7.0, 7.0, 7.0, 7.0 |
166 | , 0.0, 0.0, 1.0, 7.0, 7.0, 7.0, 7.0 ] | 164 | , 0.0, 0.0, 1.0, 7.0, 7.0, 7.0, 7.0 ] |
167 | 165 | ||
168 | -} | 166 | -} |
167 | infixl 3 ||| | ||
168 | (|||) :: Matrix Double -> Matrix Double -> Matrix Double | ||
169 | a ||| b = fromBlocks [[a,b]] | ||
170 | |||
171 | -- | a synonym for ('|||') (unicode 0x00a6, broken bar) | ||
169 | infixl 3 ¦ | 172 | infixl 3 ¦ |
170 | (¦) :: Matrix Double -> Matrix Double -> Matrix Double | 173 | (¦) :: Matrix Double -> Matrix Double -> Matrix Double |
171 | a ¦ b = fromBlocks [[a,b]] | 174 | (¦) = (|||) |
175 | |||
172 | 176 | ||
173 | -- | vertical concatenation of real matrices | 177 | -- | vertical concatenation of real matrices |
174 | -- | 178 | -- |
175 | -- (unicode 0x2014, em dash) | 179 | (===) :: Matrix Double -> Matrix Double -> Matrix Double |
180 | infixl 2 === | ||
181 | a === b = fromBlocks [[a],[b]] | ||
182 | |||
183 | -- | a synonym for ('===') (unicode 0x2014, em dash) | ||
176 | (——) :: Matrix Double -> Matrix Double -> Matrix Double | 184 | (——) :: Matrix Double -> Matrix Double -> Matrix Double |
177 | infixl 2 —— | 185 | infixl 2 —— |
178 | a —— b = fromBlocks [[a],[b]] | 186 | (——) = (===) |
187 | |||
179 | 188 | ||
180 | (#) :: Matrix Double -> Matrix Double -> Matrix Double | 189 | (#) :: Matrix Double -> Matrix Double -> Matrix Double |
181 | infixl 2 # | 190 | infixl 2 # |
@@ -356,7 +365,26 @@ pairwiseD2 x y | ok = x2 `outer` oy + ox `outer` y2 - 2* x <> trans y | |||
356 | 365 | ||
357 | -------------------------------------------------------------------------------- | 366 | -------------------------------------------------------------------------------- |
358 | 367 | ||
359 | -- | outer products of rows | 368 | {- | outer products of rows |
369 | |||
370 | >>> a | ||
371 | (3><2) | ||
372 | [ 1.0, 2.0 | ||
373 | , 10.0, 20.0 | ||
374 | , 100.0, 200.0 ] | ||
375 | >>> b | ||
376 | (3><3) | ||
377 | [ 1.0, 2.0, 3.0 | ||
378 | , 4.0, 5.0, 6.0 | ||
379 | , 7.0, 8.0, 9.0 ] | ||
380 | |||
381 | >>> rowOuters a (b ||| 1) | ||
382 | (3><8) | ||
383 | [ 1.0, 2.0, 3.0, 1.0, 2.0, 4.0, 6.0, 2.0 | ||
384 | , 40.0, 50.0, 60.0, 10.0, 80.0, 100.0, 120.0, 20.0 | ||
385 | , 700.0, 800.0, 900.0, 100.0, 1400.0, 1600.0, 1800.0, 200.0 ] | ||
386 | |||
387 | -} | ||
360 | rowOuters :: Matrix Double -> Matrix Double -> Matrix Double | 388 | rowOuters :: Matrix Double -> Matrix Double -> Matrix Double |
361 | rowOuters a b = a' * b' | 389 | rowOuters a b = a' * b' |
362 | where | 390 | where |
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs b/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs index c8c7536..c9e75de 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs | |||
@@ -16,6 +16,7 @@ module Numeric.LinearAlgebra.Util.Convolution( | |||
16 | corr2, conv2, separable | 16 | corr2, conv2, separable |
17 | ) where | 17 | ) where |
18 | 18 | ||
19 | import qualified Data.Vector.Storable as SV | ||
19 | import Data.Packed.Numeric | 20 | import Data.Packed.Numeric |
20 | 21 | ||
21 | 22 | ||
@@ -51,7 +52,7 @@ conv ker v | |||
51 | | dim ker == 0 = konst 0 (dim v) | 52 | | dim ker == 0 = konst 0 (dim v) |
52 | | otherwise = corr ker' v' | 53 | | otherwise = corr ker' v' |
53 | where | 54 | where |
54 | ker' = (flatten.fliprl.asRow) ker | 55 | ker' = SV.reverse ker |
55 | v' = vjoin [z,v,z] | 56 | v' = vjoin [z,v,z] |
56 | z = konst 0 (dim ker -1) | 57 | z = konst 0 (dim ker -1) |
57 | 58 | ||