pure-lang-svn Mailing List for Pure (Page 2)
Status: Beta
Brought to you by:
agraef
Menu
▾
▴
pure-lang-svn — SVN commits.
You can subscribe to this list here.
| 2008 |
Jan
|
Feb
|
Mar
|
Apr
(5) |
May
(141) |
Jun
(184) |
Jul
(97) |
Aug
(232) |
Sep
(196) |
Oct
|
Nov
|
Dec
|
|---|
|
From: <ag...@us...> - 2008-09-25 09:46:15
|
Revision: 860
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=860&view=rev
Author: agraef
Date: 2008-09-25 09:46:09 +0000 (Thu, 25 Sep 2008)
Log Message:
-----------
Bugfixes.
Modified Paths:
--------------
pure/trunk/runtime.cc
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-25 09:37:07 UTC (rev 859)
+++ pure/trunk/runtime.cc 2008-09-25 09:46:09 UTC (rev 860)
@@ -4470,13 +4470,14 @@
case EXPR::MATRIX: {
gsl_matrix_symbolic *m = (gsl_matrix_symbolic*)x->data.mat.p;
if (n1*n2!=m->size1*m->size2) return 0;
- if (m->tda == m->size2) {
+ if (m->tda == (m->size2>0?m->size2:1)) {
// No copying necessary, just create a new view of this matrix.
gsl_matrix_symbolic *m1 =
(gsl_matrix_symbolic*)malloc(sizeof(gsl_matrix_symbolic));
assert(m1);
*m1 = *m;
m1->size1 = n1; m1->tda = m1->size2 = n2;
+ if (m1->tda == 0) m1->tda = 1;
p = m1;
} else {
// Create a new matrix containing the same elements.
@@ -4484,6 +4485,7 @@
if (y) {
gsl_matrix_symbolic *m = (gsl_matrix_symbolic*)y->data.mat.p;
m->size1 = n1; m->tda = m->size2 = n2;
+ if (m->tda == 0) m->tda = 1;
}
return y;
}
@@ -4493,12 +4495,13 @@
case EXPR::DMATRIX: {
gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
if (n1*n2!=m->size1*m->size2) return 0;
- if (m->tda == m->size2) {
+ if (m->tda == (m->size2>0?m->size2:1)) {
// No copying necessary, just create a new view of this matrix.
gsl_matrix *m1 = (gsl_matrix*)malloc(sizeof(gsl_matrix));
assert(m1);
*m1 = *m;
m1->size1 = n1; m1->tda = m1->size2 = n2;
+ if (m1->tda == 0) m1->tda = 1;
p = m1;
} else {
// Create a new matrix containing the same elements.
@@ -4506,6 +4509,7 @@
if (y) {
gsl_matrix *m = (gsl_matrix*)y->data.mat.p;
m->size1 = n1; m->tda = m->size2 = n2;
+ if (m->tda == 0) m->tda = 1;
}
return y;
}
@@ -4514,13 +4518,14 @@
case EXPR::CMATRIX: {
gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
if (n1*n2!=m->size1*m->size2) return 0;
- if (m->tda == m->size2) {
+ if (m->tda == (m->size2>0?m->size2:1)) {
// No copying necessary, just create a new view of this matrix.
gsl_matrix_complex *m1 =
(gsl_matrix_complex*)malloc(sizeof(gsl_matrix_complex));
assert(m1);
*m1 = *m;
m1->size1 = n1; m1->tda = m1->size2 = n2;
+ if (m1->tda == 0) m1->tda = 1;
p = m1;
} else {
// Create a new matrix containing the same elements.
@@ -4528,6 +4533,7 @@
if (y) {
gsl_matrix_complex *m = (gsl_matrix_complex*)y->data.mat.p;
m->size1 = n1; m->tda = m->size2 = n2;
+ if (m->tda == 0) m->tda = 1;
}
return y;
}
@@ -4536,13 +4542,14 @@
case EXPR::IMATRIX: {
gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
if (n1*n2!=m->size1*m->size2) return 0;
- if (m->tda == m->size2) {
+ if (m->tda == (m->size2>0?m->size2:1)) {
// No copying necessary, just create a new view of this matrix.
gsl_matrix_int *m1 =
(gsl_matrix_int*)malloc(sizeof(gsl_matrix_int));
assert(m1);
*m1 = *m;
m1->size1 = n1; m1->tda = m1->size2 = n2;
+ if (m1->tda == 0) m1->tda = 1;
p = m1;
} else {
// Create a new matrix containing the same elements.
@@ -4550,6 +4557,7 @@
if (y) {
gsl_matrix_int *m = (gsl_matrix_int*)y->data.mat.p;
m->size1 = n1; m->tda = m->size2 = n2;
+ if (m->tda == 0) m->tda = 1;
}
return y;
}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-25 09:37:16
|
Revision: 859
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=859&view=rev
Author: agraef
Date: 2008-09-25 09:37:07 +0000 (Thu, 25 Sep 2008)
Log Message:
-----------
Comment change.
Modified Paths:
--------------
pure/trunk/examples/gauss.pure
Modified: pure/trunk/examples/gauss.pure
===================================================================
--- pure/trunk/examples/gauss.pure 2008-09-25 09:36:44 UTC (rev 858)
+++ pure/trunk/examples/gauss.pure 2008-09-25 09:37:07 UTC (rev 859)
@@ -17,12 +17,13 @@
/* Here is a brief rundown of the algorithm: First we find the pivot element
in column j of the matrix. (We're doing partial pivoting here, i.e., we
- only look for the pivot in column j, starting at row i.) If the pivot is
- zero then we're done (the pivot column is already zeroed out). Otherwise,
- we bring it into the pivot position (swapping row i and the pivot row),
- divide the picot row by the pivot, and subtract suitable multiples of the
- pivot row to eliminate the pivot column in all subsequent rows. Finally we
- update i and p accordingly and return the result. */
+ only look for the pivot in column j, starting at row i. That's usually good
+ enough to achieve numerical stability.) If the pivot is zero then we're
+ done (the pivot column is already zeroed out). Otherwise, we bring it into
+ the pivot position (swapping row i and the pivot row), divide the pivot row
+ by the pivot, and subtract suitable multiples of the pivot row to eliminate
+ the elements of the pivot column in all subsequent rows. Finally we update
+ i and p accordingly and return the result. */
step (p,i,x) j
= if max_x>0 then
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-25 09:36:52
|
Revision: 858
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=858&view=rev
Author: agraef
Date: 2008-09-25 09:36:44 +0000 (Thu, 25 Sep 2008)
Log Message:
-----------
Update documentation.
Modified Paths:
--------------
pure/trunk/pure.1.in
Modified: pure/trunk/pure.1.in
===================================================================
--- pure/trunk/pure.1.in 2008-09-25 09:19:21 UTC (rev 857)
+++ pure/trunk/pure.1.in 2008-09-25 09:36:44 UTC (rev 858)
@@ -1239,9 +1239,9 @@
vectors aren't the same size then you'll get an `out_of_bounds' exception with
the definition above.)
.PP
-The matrix product now boils down to a simple matrix comprehension which just
-multiplies all rows of x with all columns of y (the rows and cols functions
-are prelude operations found in matrices.pure):
+The general matrix product now boils down to a simple matrix comprehension
+which just multiplies all rows of x with all columns of y (the rows and cols
+functions are prelude operations found in matrices.pure):
.sp
.nf
> x::matrix * y::matrix = {u*v | u = rows x; v = cols y};
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-25 09:22:02
|
Revision: 857
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=857&view=rev
Author: agraef
Date: 2008-09-25 09:19:21 +0000 (Thu, 25 Sep 2008)
Log Message:
-----------
Update documentation.
Modified Paths:
--------------
pure/trunk/pure.1.in
Modified: pure/trunk/pure.1.in
===================================================================
--- pure/trunk/pure.1.in 2008-09-25 09:06:41 UTC (rev 856)
+++ pure/trunk/pure.1.in 2008-09-25 09:19:21 UTC (rev 857)
@@ -1189,7 +1189,7 @@
a amounts to mapping the function \ex->a*x to x, which can be done as follows:
.sp
.nf
-> a::int * x::matrix = map (\ex->a*x) x;
+> a * x::matrix = map (\ex->a*x) x \fBif\fP not matrixp a;
> 2*{1,2,3;4,5,6};
{2,4,6;8,10,12}
.fi
@@ -1254,27 +1254,20 @@
equations. The algorithm brings a matrix into ``row echelon'' form, a
generalization of triangular matrices. The resulting system can then be solved
quite easily using back substitution. Here is a Pure implementation of the
-algorithm (please refer to any good textbook on numeric mathematics for a
-closer description of the algorithm):
+algorithm:
.sp
.nf
gauss_elimination x::matrix = p,x
\fBwhen\fP n,m = dim x; p,_,x = foldl step (0..n-1,0,x) (0..m-1) \fBend\fP;
-.fi
-.PP
-The actual pivoting and elimination step is a bit involved. x is our matrix, i
-the current row index, j the current column index, and p keeps track of the
-current permutation of the row indices performed during pivoting. The
-algorithm returns the updated matrix x, row index i and row permutation p.
-.sp
-.nf
+
+// One pivoting and elimination step in column j of the matrix:
step (p,i,x) j
= \fBif\fP max_x>0 \fBthen\fP
// updated row permutation and index:
transp i max_i p, i+1,
{// the top rows of the matrix remain unchanged:
x!!(0..i-1,0..m-1);
- // the pivot row, divided by the pivot:
+ // the pivot row, divided by the pivot element:
{x!(i,l)/x!(i,j) | l=0..m-1};
// subtract suitable multiples of the pivot row:
{x!(k,l)-x!(k,j)*x!(i,l)/x!(i,j) | k=i+1..n-1; l=0..m-1}}
@@ -1288,9 +1281,28 @@
\fBend\fP;
.fi
.PP
-We also need the following little helper functions to swap two rows of a
-matrix (this is used in the pivoting step above) and to apply a transposition
-to a permutation (represented as a list):
+The real meat is in the pivoting and elimination step ('step' function) which
+is iterated over all columns of the input matrix. In each step, x is the
+current matrix, i the current row index, j the current column index, and p
+keeps track of the current permutation of the row indices performed during
+pivoting. The algorithm returns the updated matrix x, row index i and row
+permutation p.
+.PP
+Please refer to any good textbook on numeric mathematics for a closer
+description of the algorithm. But here is a brief rundown of what happens in
+each elimination step: First we find the pivot element in column j of the
+matrix. (We're doing partial pivoting here, i.e., we only look for the element
+with the largest absolute value in column j, starting at row i. That's usually
+good enough to achieve numerical stability.) If the pivot is zero then we're
+done (the rest of the pivot column is already zeroed out). Otherwise, we bring
+it into the pivot position (swapping row i and the pivot row), divide the
+pivot row by the pivot, and subtract suitable multiples of the pivot row to
+eliminate the elements of the pivot column in all subsequent rows. Finally we
+update i and p accordingly and return the result.
+.PP
+In order to complete the implementation, we still need the following little
+helper functions to swap two rows of a matrix (this is used in the pivoting
+step) and to apply a transposition to a permutation (represented as a list):
.sp
.nf
swap x i j = x!!(transp i j (0..n-1),0..m-1) \fBwhen\fP n,m = dim x \fBend\fP;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-25 09:09:09
|
Revision: 856
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=856&view=rev
Author: agraef
Date: 2008-09-25 09:06:41 +0000 (Thu, 25 Sep 2008)
Log Message:
-----------
Comment change.
Modified Paths:
--------------
pure/trunk/examples/gauss.pure
Modified: pure/trunk/examples/gauss.pure
===================================================================
--- pure/trunk/examples/gauss.pure 2008-09-25 09:03:51 UTC (rev 855)
+++ pure/trunk/examples/gauss.pure 2008-09-25 09:06:41 UTC (rev 856)
@@ -30,7 +30,7 @@
transp i max_i p, i+1,
{// the top rows of the matrix remain unchanged:
x!!(0..i-1,0..m-1);
- // the pivot row, divided by the pivot:
+ // the pivot row, divided by the pivot element:
{x!(i,l)/x!(i,j) | l=0..m-1};
// subtract suitable multiples of the pivot row:
{x!(k,l)-x!(k,j)*x!(i,l)/x!(i,j) | k=i+1..n-1; l=0..m-1}}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-25 09:04:52
|
Revision: 855
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=855&view=rev
Author: agraef
Date: 2008-09-25 09:03:51 +0000 (Thu, 25 Sep 2008)
Log Message:
-----------
Bugfixes.
Modified Paths:
--------------
pure/trunk/lib/matrices.pure
Modified: pure/trunk/lib/matrices.pure
===================================================================
--- pure/trunk/lib/matrices.pure 2008-09-25 07:40:27 UTC (rev 854)
+++ pure/trunk/lib/matrices.pure 2008-09-25 09:03:51 UTC (rev 855)
@@ -266,14 +266,15 @@
= matrix_redim x n m if n>=0 && m>=0 && n*m==#x;
/* You can also redim a matrix to a given row size. In this case the row size
- must be positive and divide the total size of the matrix, */
+ must divide the total size of the matrix, */
redim x::matrix m::int = redim x (#x div m,m) if m>0 && #x mod m==0;
+ = x if m==0 && #x==0;
/* Convert a matrix to a row or column vector. */
rowvector x::matrix = redim x (1,#x);
-colvector x::matrix = redim x 1;
+colvector x::matrix = redim x (#x,1);
/* Transpose a matrix. */
@@ -331,19 +332,25 @@
takewhile p x::matrix = colcat (takewhile p (list x));
tail x::matrix = x!!(1..#x-1) if not null x;
-zip x::matrix y::matrix = flip redim (dim x!1) $
+private zipdim zip3dim;
+zipdim x::matrix y::matrix
+ = min (dim x!0) (dim y!0),dim x!1;
+zip3dim x::matrix y::matrix z::matrix
+ = min (dim x!0) (min (dim y!0) (dim z!0)),dim x!1;
+
+zip x::matrix y::matrix = flip redim (zipdim x y) $
colcat (zip (list x) (list y))
if dim x!1==dim y!1;
zip3 x::matrix y::matrix z::matrix
- = flip redim (dim x!1) $
+ = flip redim (zip3dim x y z) $
colcat (zip3 (list x) (list y) (list z))
if dim x!1==dim y!1 && dim x!1==dim z!1;
zipwith f x::matrix y::matrix
- = flip redim (dim x!1) $
+ = flip redim (zipdim x y) $
colcat (zipwith f (list x) (list y))
if dim x!1==dim y!1;
zipwith3 f x::matrix y::matrix z::matrix
- = flip redim (dim x!1) $
+ = flip redim (zip3dim x y z) $
colcat (zipwith3 f (list x) (list y) (list z))
if dim x!1==dim y!1 && dim x!1==dim z!1;
dowith f x::matrix y::matrix
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-25 07:40:34
|
Revision: 854
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=854&view=rev
Author: agraef
Date: 2008-09-25 07:40:27 +0000 (Thu, 25 Sep 2008)
Log Message:
-----------
Bugfixes.
Modified Paths:
--------------
pure/trunk/etc/pure-mode.el.in
Modified: pure/trunk/etc/pure-mode.el.in
===================================================================
--- pure/trunk/etc/pure-mode.el.in 2008-09-25 07:31:35 UTC (rev 853)
+++ pure/trunk/etc/pure-mode.el.in 2008-09-25 07:40:27 UTC (rev 854)
@@ -126,7 +126,7 @@
:group 'pure)
(defcustom pure-msg-regexp
- "^[ \t]*\\(\\([^:\n]+\\):\\([0-9]+\\)\\(\\.[0-9]+\\)?\\): "
+ "^[ \t]*\\(\\([^:\n]+\\):\\([0-9]+\\)\\(\\.[0-9]+\\(-[0-9]+\\)?\\)?\\): "
"*Regexp to match error and warning messages with source line references in
the Pure eval buffer. Expression 1 denotes the whole source line info,
expression 2 the file name and expression 3 the corresponding line number."
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-25 07:31:37
|
Revision: 853
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=853&view=rev
Author: agraef
Date: 2008-09-25 07:31:35 +0000 (Thu, 25 Sep 2008)
Log Message:
-----------
Bugfixes.
Modified Paths:
--------------
pure/trunk/lib/matrices.pure
Modified: pure/trunk/lib/matrices.pure
===================================================================
--- pure/trunk/lib/matrices.pure 2008-09-25 01:02:54 UTC (rev 852)
+++ pure/trunk/lib/matrices.pure 2008-09-25 07:31:35 UTC (rev 853)
@@ -112,12 +112,14 @@
of int values. As with list slicing, index ranges may be non-contiguous
and/or non-monotonous. However, the case of contiguous and monotonous
ranges is optimized by making good use of the 'submat' operation below. We
- also add some rules to make slicing work with ranges drawn from a matrix
- rather than a list. */
+ also add some convenience rules to handle matrix ranges as well "mixed"
+ ranges (ns,ms) where either ns or ms is a singleton. */
x!!ns::matrix = x!!list ns;
-x!!(ns::matrix,ms::matrix)
- = x!!(list ns,list ms);
+x!!(ns::matrix,ms) = x!!(list ns,ms);
+x!!(ns,ms::matrix) = x!!(ns,list ms);
+x!!(ns::int,ms) = x!!([ns],ms);
+x!!(ns,ms::int) = x!!(ns,[ms]);
x::matrix!!(ns,ms) = case ns,ms of
ns@(n:_),ms@(m:_) = submat x (n,m) (#ns,#ms)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-25 01:03:04
|
Revision: 852
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=852&view=rev
Author: agraef
Date: 2008-09-25 01:02:54 +0000 (Thu, 25 Sep 2008)
Log Message:
-----------
Updated NEWS.
Modified Paths:
--------------
pure/trunk/NEWS
Modified: pure/trunk/NEWS
===================================================================
--- pure/trunk/NEWS 2008-09-25 00:55:11 UTC (rev 851)
+++ pure/trunk/NEWS 2008-09-25 01:02:54 UTC (rev 852)
@@ -2,59 +2,38 @@
** Pure 0.7 (in progress)
Basic GSL (GNU Scientific Library) matrix support has been implemented,
-including matrix comprehensions. Here's some preliminary documentation (at the
-time of this writing, the manual still needs to be updated). For more
-information on GSL please refer to http://www.gnu.org/software/gsl.
+including matrix comprehensions. Here's a brief overview of the new features.
+For more information on GSL please refer to http://www.gnu.org/software/gsl.
GSL double, complex and integer matrices can be created with the new {x,y;u,v}
syntax, which works more or less like Octave/MATLAB matrices, but using curly
braces instead of brackets. Also, indices are zero-based (rather than
-one-based) to be consistent with Pure's list and tuple structures.
+one-based) to be consistent with Pure's list and tuple structures. Here are
+some simple examples of matrices:
-Here are some simple examples of matrices:
-
- {1,2,3} a row vector consisting of machine ints
- {1.0;2.0;3.0} a column vector of double values (';' separates rows)
- {1,2;3,4} a 2x2 int matrix
-- {1L,y+1;foo,bar} a "symbolic" matrix (see explanation below)
+- {1L,y+1;foo,bar} a symbolic matrix
-Note that the {...} syntax can be used only on the right-hand side of a
-definition to construct new matrices. Matrix patterns are not supported.
-However, the new 'matrix' type tag can be used in patterns to match matrix
-values, e.g.:
-
-foo a::matrix = {a!(i,j)+1 | i=0..n-1; j=0..m-1} when n,m = dim a end;
-
-(Additional predicates to check for the different subtypes of matrices are
-available in the prelude.)
-
-As mentioned above, there's also a special built-in type of "symbolic"
-matrices. These work like the numeric matrices, but can contain an arbitrary
-mixture of Pure values, and also work if GSL support isn't available.
-
-GSL matrices are always homogeneous, i.e., they only contain values from one
+Symbolic matrices work like the numeric matrices, but can contain an arbitrary
+mixture of Pure values, and also work if GSL support isn't available. GSL
+matrices are always homogeneous, i.e., they only contain values from one
numeric type, which in the Pure implementation can be machine ints, double and
-complex double values. (The latter are represented using Pure's infix notation
-for complex values defined in math.pure; this requires math.pure to be
-loaded. The same notation is also used to print complex matrices. Note that,
-for performance reasons, the expression printer doesn't use the __show__
-function to print individual matrix elements. It is possible to override the
-default print representation of matrix values as a whole, though.)
+complex double values. If a matrix contains values of different types, or Pure
+values which cannot be stored in a GSL matrix, then a symbolic matrix is
+created instead (this also includes the case of bigints, which are considered
+as symbolic values as far as matrix construction is concerned). If the
+interpreter was built without GSL support then symbolic matrices are the only
+kind of matrices supported by the interpreter.
-If a matrix contains values of different types, or Pure values which cannot be
-stored in a GSL matrix, then a symbolic matrix is created instead (this also
-includes the case of bigints, which are considered as symbolic values as far
-as matrix construction is concerned). If the interpreter was built without GSL
-support then symbolic matrices are the only kind of matrices supported by the
-interpreter.
-
Pure also provides so-called matrix comprehensions as a convenient means to
create matrices from a template expression (which can denote either a scalar
or a submatrix), drawing values from lists and (optionally) filtering the
elements with predicates. These work pretty much like list comprehensions, but
return matrices instead of lists. Generator clauses in matrix comprehensions
-alternate between row and column generation so that most common mathematical
-abbreviations carry over quite easily. E.g., here's a simple example which
+alternate between row and column generation so that customary mathematical
+notation carries over quite easily. E.g., here's a simple example which
illustrates how you can define a function which returns a square identity
matrix of a given dimension:
@@ -70,7 +49,7 @@
and inspect Pure matrix values from external C modules. Moreover, the
'pointer' function in the prelude can be used to convert matrices to Pure
pointer values, and marshalling of GSL matrices in the C interface is also
-available.
+provided.
Here is a brief synopsis of some important operations which are implemented in
the prelude (you can find the definitions of these and a bunch of other matrix
@@ -84,7 +63,9 @@
- x!!is, x!!(is,js) slicing (is and js are lists of machine ints)
- x==y, x!=y matrix comparisons
-Adding other operations by interfacing to GSL should be a piece of cake.
+Adding other operations by interfacing to GSL should be a piece of cake. In
+fact, we plan to provide a comprehensive Pure interface to GSL as a separate
+library in the future.
Other user-visible changes prompted by the introduction of the matrix syntax
are listed below:
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-25 00:55:21
|
Revision: 851
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=851&view=rev
Author: agraef
Date: 2008-09-25 00:55:11 +0000 (Thu, 25 Sep 2008)
Log Message:
-----------
Add Gaussian elimination example.
Modified Paths:
--------------
pure/trunk/NEWS
Added Paths:
-----------
pure/trunk/examples/gauss.pure
Modified: pure/trunk/NEWS
===================================================================
--- pure/trunk/NEWS 2008-09-25 00:54:30 UTC (rev 850)
+++ pure/trunk/NEWS 2008-09-25 00:55:11 UTC (rev 851)
@@ -78,7 +78,7 @@
- #x total number of elements
- dim x number of rows and columns (as a pair)
-- x' transposed matrix
+- x' transpose of a matrix
- x!i ith matrix element in row-major order (zero-based)
- x!(i,j) jth element in ith row of the matrix (zero-based)
- x!!is, x!!(is,js) slicing (is and js are lists of machine ints)
Added: pure/trunk/examples/gauss.pure
===================================================================
--- pure/trunk/examples/gauss.pure (rev 0)
+++ pure/trunk/examples/gauss.pure 2008-09-25 00:55:11 UTC (rev 851)
@@ -0,0 +1,69 @@
+
+/* Pure matrix example: Gaussian elimination algorithm. 2008-09-25 AG */
+
+/* The Gaussian elimination algorithm brings a matrix into row echelon form,
+ which makes it possible to solve the original system of linear equations
+ using back substitution. Our version of the algorithm just returns the row
+ echelon form of the matrix, along with the corresponding permutation of the
+ row indices performed during pivoting. */
+
+gauss_elimination x::matrix = p,x
+when n,m = dim x; p,_,x = foldl step (0..n-1,0,x) (0..m-1) end;
+
+/* The elimination step. x is our matrix, i the current row index, j the
+ current column index, and p keeps track of the current permutation of the
+ row indices performed during pivoting. The algorithm returns the updated
+ matrix x, row index i and row permutation p. */
+
+/* Here is a brief rundown of the algorithm: First we find the pivot element
+ in column j of the matrix. (We're doing partial pivoting here, i.e., we
+ only look for the pivot in column j, starting at row i.) If the pivot is
+ zero then we're done (the pivot column is already zeroed out). Otherwise,
+ we bring it into the pivot position (swapping row i and the pivot row),
+ divide the picot row by the pivot, and subtract suitable multiples of the
+ pivot row to eliminate the pivot column in all subsequent rows. Finally we
+ update i and p accordingly and return the result. */
+
+step (p,i,x) j
+= if max_x>0 then
+ // updated row permutation and index:
+ transp i max_i p, i+1,
+ {// the top rows of the matrix remain unchanged:
+ x!!(0..i-1,0..m-1);
+ // the pivot row, divided by the pivot:
+ {x!(i,l)/x!(i,j) | l=0..m-1};
+ // subtract suitable multiples of the pivot row:
+ {x!(k,l)-x!(k,j)*x!(i,l)/x!(i,j) | k=i+1..n-1; l=0..m-1}}
+ else p,i,x
+when
+ n,m = dim x; max_i, max_x = pivot i (col x j);
+ x = if max_x>0 then swap x i max_i else x;
+end with
+ pivot i x = foldl max (0,0) [j,abs (x!j)|j=i..#x-1];
+ max (i,x) (j,y) = if x<y then j,y else i,x;
+end;
+
+/* Swap rows i and j of the matrix x. This operation is used in the partial
+ pivoting step. (This is not really needed because we could just index the
+ rows indirectly through the row permutation readily available in our
+ implementation of the algorithm, but we omitted this here for clarity.) */
+
+swap x i j = x!!(transp i j (0..n-1),0..m-1) when n,m = dim x end;
+
+/* A little helper function to apply a transposition to a permutation. */
+
+transp i j p = [p!tr k | k=0..#p-1]
+with tr k = if k==i then j else if k==j then i else k end;
+
+/* For convenience, print a double matrix in "short" format a la Octave. */
+
+using system;
+__show__ x::matrix
+= strcat [printd j (x!(i,j))|i=0..n-1; j=0..m-1] + "\n"
+with printd 0 = sprintf "\n%10.5f"; printd _ = sprintf "%10.5f" end
+when n,m = dim x end if dmatrixp x;
+
+/* Example: */
+
+let x = dmatrix {2,1,-1,8; -3,-1,2,-11; -2,1,2,-3};
+x; gauss_elimination x;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-25 00:54:41
|
Revision: 850
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=850&view=rev
Author: agraef
Date: 2008-09-25 00:54:30 +0000 (Thu, 25 Sep 2008)
Log Message:
-----------
Update documentation.
Modified Paths:
--------------
pure/trunk/pure.1.in
Modified: pure/trunk/pure.1.in
===================================================================
--- pure/trunk/pure.1.in 2008-09-24 13:08:16 UTC (rev 849)
+++ pure/trunk/pure.1.in 2008-09-25 00:54:30 UTC (rev 850)
@@ -336,6 +336,8 @@
.B false
and any non-zero value
.BR true ).
+Pure also provides some built-in support for lists and matrices, although most
+of the corresponding operations are actually defined in the prelude.
.SS Expression Syntax
Expressions consist of the following elements:
.TP
@@ -430,15 +432,48 @@
tuple 1,2, the integer 3, and another tuple 4,5. Likewise, [(1,2,3)] is list
with a single element, the tuple 1,2,3.
.TP
-\fBList comprehensions:\fP [x,y; x = 1..n; y = 1..m; x<y]
-Pure also has list comprehensions which generate lists from an expression and
-one or more ``generator'' and ``filter'' clauses (the former bind a pattern to
-values drawn from a list, the latter are just predicates determining which
-generated elements should actually be added to the output list). List
-comprehensions are in fact syntactic sugar for a combination of nested
-lambdas, conditional expressions and ``catmaps'' (a list operation which
-combines list concatenation and mapping a function over a list, defined in the
-prelude), but they are often much easier to write. Some examples of list
+\fBMatrices:\fP {1.0,2.0,3.0}, {1,2;3,4}, {1L,y+1;foo,bar}
+Pure also offers matrices, a kind of arrays, as a built-in data structure
+which provides efficient storage and element access. These work more or less
+like their Octave/MATLAB equivalents, but using curly braces instead of
+brackets. As indicated, commas are used to separate the columns of a matrix,
+semicolons for its rows. In fact, the {...} construct is rather general,
+allowing you to construct new matrices from individual elements and/or
+submatrices, provided that all dimensions match up. E.g., {{1;3},{2;4}} is
+another way to write a 2x2 matrix in ``column-major'' form (however,
+internally all matrices are stored in C's row-major format).
+.sp
+If the interpreter was built with support for the GNU Scientific Library (GSL)
+then both numeric and symbolic matrices are available. The former are thin
+wrappers around GSL's homogeneous arrays of double, complex double or
+(machine) int matrices, while the latter can contain any mixture of Pure
+expressions. Pure will pick the appropriate type for the data at hand. If a
+matrix contains values of different types, or Pure values which cannot be
+stored in a numeric matrix, then a symbolic matrix is created instead (this
+also includes the case of bigints, which are considered as symbolic values as
+far as matrix construction is concerned). If the interpreter was built without
+GSL support then symbolic matrices are the only kind of matrices supported by
+the interpreter.
+.sp
+More information about matrices and corresponding examples can be found in the
+EXAMPLES section below.
+.TP
+\fBComprehensions:\fP [x,y | x=1..n; y=1..m; x<y], {i!=j | i=1..n; j=1..m}
+Pure provides the usual comprehension syntax as a convenient means to
+construct both list and matrix values from a ``template'' expression and one
+or more ``generator'' and ``filter'' clauses (the former bind a pattern to
+values drawn from a list or matrix, the latter are just predicates determining
+which generated elements should actually be added to the result). Both list
+and matrix comprehensions are in fact syntactic sugar for a combination of
+nested lambdas, conditional expressions and ``catmaps'' (a collection of
+operations which combine list or matrix construction and mapping a function
+over a list or matrix, defined in the prelude), but they are often much easier
+to write.
+.sp
+Matrix comprehensions work pretty much like list comprehensions, but produce
+matrices instead of lists. Generator clauses in matrix comprehensions
+alternate between row and column generation so that most common mathematical
+abbreviations carry over quite easily. Examples of both kinds of
comprehensions can be found in the EXAMPLES section below.
.TP
\fBFunction applications:\fP foo x y z
@@ -712,9 +747,113 @@
symbols needed in an evaluation
.I before
entering the expression to be evaluated.
+.SH RULE SYNTAX
+Basically, the same rule syntax is used in all kinds of global and local
+definitions. However, some constructs (specifically, \fBwhen\fP, \fBlet\fP,
+\fBconst\fP and \fBdef\fP) use a restricted rule syntax where no guards or
+multiple left-hand and right-hand sides are permitted. When matching against a
+function or macro call, or the subject term in a \fBcase\fP expression, the
+rules are always considered in the order in which they are written, and the
+first matching rule (whose guard evaluates to a nonzero value, if applicable)
+is picked. (Again, the \fBwhen\fP construct is treated differently, because
+each rule is actually a separate definition.)
+.PP
+In any case, the left-hand side pattern (which, as already mentioned, is
+always a simple expression) must not contain repeated variables (i.e., rules
+must be ``left-linear''), except for the anonymous variable `_' which matches
+an arbitrary value without binding a variable symbol.
+.PP
+A left-hand side variable (including the anonymous variable) may be followed
+by one of the special type tags \fB::int\fP, \fB::bigint\fP, \fB::double\fP,
+\fB::string\fP, \fB::matrix\fP, \fB::pointer\fP, to indicate that it can only
+match a constant value of the corresponding built-in type. (This is useful if
+you want to write rules matching \fIany\fP object of one of these types; note
+that there is no way to write out all ``constructors'' for the built-in types,
+as there are infinitely many.)
+.PP
+Pure also supports Haskell-style ``as'' patterns of the form
+.IB variable @ pattern
+which binds the given variable to the expression matched by the subpattern
+.I pattern
+(in addition to the variables bound by
+.I pattern
+itself). This is convenient if the value matched by the subpattern is to be
+used on the right-hand side of an equation. Syntactically, ``as'' patterns are
+primary expressions; if the subpattern is not a primary expression, it must be
+parenthesized. For instance, the following function duplicates the head
+element of a list:
+.sp
+.nf
+foo xs@(x:_) = x:xs;
+.fi
+.PP
+The left-hand side of a rule can be omitted if it is the same as for the
+previous rule. This provides a convenient means to write out a collection of
+equations for the same left-hand side which discriminates over different
+conditions:
+.sp
+.nf
+\fIlhs\fR = \fIrhs\fP \fBif\fP \fIguard\fP;
+ = \fIrhs\fP \fBif\fP \fIguard\fP;
+ ...
+ = \fIrhs\fP \fBotherwise\fP;
+.fi
+.PP
+For instance:
+.sp
+.nf
+fact n = n*fact (n-1) \fBif\fP n>0;
+ = 1 \fBotherwise\fP;
+.fi
+.PP
+Pure also allows a collection of rules with different left-hand sides but the
+same right-hand side(s) to be abbreviated as follows:
+.sp
+.nf
+\fIlhs\fP |
+ ...
+\fIlhs\fP = \fIrhs\fP;
+.fi
+.PP
+This is useful if you need different specializations of the same rule which
+use different type tags on the left-hand side variables. For instance:
+.sp
+.nf
+fact n::int |
+fact n::double |
+fact n = n*fact(n-1) \fBif\fP n>0;
+ = 1 \fBotherwise\fP;
+.fi
+.PP
+In fact, the left-hand sides don't have to be related at all, so that you can
+also write something like:
+.sp
+.nf
+foo x | bar y = x*y;
+.fi
+.PP
+However, this is most useful when using an ``as'' pattern to bind a common
+variable to a parameter value
+.I after
+checking that it matches one of several possible argument patterns (which is
+slightly more efficient than using an equivalent type-checking guard). E.g.,
+the following definition binds the xs variable to the parameter of foo, if it
+is either the empty list or a list starting with an integer:
+.sp
+.nf
+foo xs@[] | foo xs@(_::int:_) = ... xs ...;
+.fi
+.PP
+The same construct also works in
+.B case
+expressions, which is convenient if different cases should be mapped to the
+same value, e.g.:
+.sp
+.nf
+\fBcase\fP ans \fBof\fP "y" | "Y" = 1; _ = 0; \fBend\fP;
+.fi
.SH EXAMPLES
-Here are a few examples of simple Pure programs (see the following section for
-a closer discussion of the rule syntax).
+Here are a few examples of simple Pure programs.
.PP
The factorial:
.sp
@@ -1032,110 +1171,159 @@
the number of elements printed until memory is exhausted. Calling `do' on a
fresh instance of the stream of primes allows `do' to get rid of each `cons'
cell after having printed the corresponding stream element.)
-.SH RULE SYNTAX
-Basically, the same rule syntax is used in all kinds of global and local
-definitions. However, some constructs (specifically, \fBwhen\fP, \fBlet\fP,
-\fBconst\fP and \fBdef\fP) use a restricted rule syntax where no guards or
-multiple left-hand and right-hand sides are permitted. When matching against a
-function or macro call, or the subject term in a \fBcase\fP expression, the
-rules are always considered in the order in which they are written, and the
-first matching rule (whose guard evaluates to a nonzero value, if applicable)
-is picked. (Again, the \fBwhen\fP construct is treated differently, because
-each rule is actually a separate definition.)
+.SS Matrix Computations
+Pure offers a number of basic matrix operations, such as matrix construction,
+indexing, slicing, as well as getting the size and dimensions of a matrix
+(these are briefly described in the STANDARD LIBRARY section below). However,
+it does
+.I not
+supply built-in support for matrix arithmetic and other linear algebra
+algorithms. The idea is that these can and should be provided through separate
+libraries, such as a GSL interface (which will hopefully be available in the
+near future).
.PP
-In any case, the left-hand side pattern (which, as already mentioned, is
-always a simple expression) must not contain repeated variables (i.e., rules
-must be ``left-linear''), except for the anonymous variable `_' which matches
-an arbitrary value without binding a variable symbol.
+But Pure's facilities for matrix and list processing also make it easy to roll
+your own, if desired. First, the prelude provides matrix versions of the
+common list operations like map, fold, zip etc., which provide a way to
+implement common matrix operations. E.g., multiplying a matrix x with a scalar
+a amounts to mapping the function \ex->a*x to x, which can be done as follows:
+.sp
+.nf
+> a::int * x::matrix = map (\ex->a*x) x;
+> 2*{1,2,3;4,5,6};
+{2,4,6;8,10,12}
+.fi
.PP
-A left-hand side variable (including the anonymous variable) may be followed
-by one of the special type tags \fB::int\fP, \fB::bigint\fP, \fB::double\fP,
-\fB::string\fP, to indicate that it can only match a constant value of the
-corresponding built-in type. (This is useful if you want to write rules
-matching \fIany\fP object of one of these types; note that there is no way to
-write out all ``constructors'' for the built-in types, as there are infinitely
-many.)
+Likewise, matrix addition and other element-wise operations can be realized
+using zipwith, which combines corresponding elements of two matrices using a
+given binary function:
+.sp
+.nf
+> x::matrix + y::matrix = zipwith (+) x y;
+> {1,2,3;4,5,6}+{1,2,1;3,2,3};
+{2,4,4;7,7,9}
+.fi
.PP
-Pure also supports Haskell-style ``as'' patterns of the form
-.IB variable @ pattern
-which binds the given variable to the expression matched by the subpattern
-.I pattern
-(in addition to the variables bound by
-.I pattern
-itself). This is convenient if the value matched by the subpattern is to be
-used on the right-hand side of an equation. Syntactically, ``as'' patterns are
-primary expressions; if the subpattern is not a primary expression, it must be
-parenthesized. For instance, the following function duplicates the head
-element of a list:
+Second, matrix comprehensions make it easy to express a variety of algorithms
+which would be implemented using `for' loops in conventional programming
+languages. To illustrate the use of matrix comprehensions, here is how we can
+define an operation to create a square identity matrix of a given dimension:
.sp
.nf
-foo xs@(x:_) = x:xs;
+> eye n = {i==j | i = 1..n; j = 1..n};
+> eye 3;
+{1,0,0;0,1,0;0,0,1}
.fi
.PP
-The left-hand side of a rule can be omitted if it is the same as for the
-previous rule. This provides a convenient means to write out a collection of
-equations for the same left-hand side which discriminates over different
-conditions:
+Note that the i==j term is just a Pure idiom for the Kronecker symbol. Another
+point worth mentioning here is that the generator clauses of matrix
+comprehensions alternate between row and column generation
+automatically. (More precisely, the last generator, which varies most quickly,
+always yields a row, the next-to-last one a column of these row vectors, and
+so on.) This makes matrix comprehensions resemble customary mathematical
+notation very closely.
+.PP
+As a slightly more comprehensive example (no pun intended!), here is a
+definition of matrix multiplication in Pure. Let's start out with the simple
+case of the ``dot'' product of two vectors:
.sp
.nf
-\fIlhs\fR = \fIrhs\fP \fBif\fP \fIguard\fP;
- = \fIrhs\fP \fBif\fP \fIguard\fP;
- ...
- = \fIrhs\fP \fBotherwise\fP;
+> x::matrix * y::matrix = sum [x!i*y!i | i=0..#x-1]
+> \fBif\fP vectorp x && vectorp y;
+> sum = foldl (+) 0;
+> {1,2,3}*{1,0,1};
+4
.fi
.PP
-For instance:
+(For the sake of simplicity, this doesn't do much error checking; if the two
+vectors aren't the same size then you'll get an `out_of_bounds' exception with
+the definition above.)
+.PP
+The matrix product now boils down to a simple matrix comprehension which just
+multiplies all rows of x with all columns of y (the rows and cols functions
+are prelude operations found in matrices.pure):
.sp
.nf
-fact n = n*fact (n-1) \fBif\fP n>0;
- = 1 \fBotherwise\fP;
+> x::matrix * y::matrix = {u*v | u = rows x; v = cols y};
+> {0,1;1,0;1,1}*{1,2,3;4,5,6};
+{4,5,6;1,2,3;5,7,9}
.fi
.PP
-Pure also allows a collection of rules with different left-hand sides but the
-same right-hand side(s) to be abbreviated as follows:
+Well, that was easy. So let's take a look at a more challenging example,
+Gaussian elimination, which can be used to solve systems of linear
+equations. The algorithm brings a matrix into ``row echelon'' form, a
+generalization of triangular matrices. The resulting system can then be solved
+quite easily using back substitution. Here is a Pure implementation of the
+algorithm (please refer to any good textbook on numeric mathematics for a
+closer description of the algorithm):
.sp
.nf
-\fIlhs\fP |
- ...
-\fIlhs\fP = \fIrhs\fP;
+gauss_elimination x::matrix = p,x
+\fBwhen\fP n,m = dim x; p,_,x = foldl step (0..n-1,0,x) (0..m-1) \fBend\fP;
.fi
.PP
-This is useful if you need different specializations of the same rule which
-use different type tags on the left-hand side variables. For instance:
+The actual pivoting and elimination step is a bit involved. x is our matrix, i
+the current row index, j the current column index, and p keeps track of the
+current permutation of the row indices performed during pivoting. The
+algorithm returns the updated matrix x, row index i and row permutation p.
.sp
.nf
-fact n::int |
-fact n::double |
-fact n = n*fact(n-1) \fBif\fP n>0;
- = 1 \fBotherwise\fP;
+step (p,i,x) j
+= \fBif\fP max_x>0 \fBthen\fP
+ // updated row permutation and index:
+ transp i max_i p, i+1,
+ {// the top rows of the matrix remain unchanged:
+ x!!(0..i-1,0..m-1);
+ // the pivot row, divided by the pivot:
+ {x!(i,l)/x!(i,j) | l=0..m-1};
+ // subtract suitable multiples of the pivot row:
+ {x!(k,l)-x!(k,j)*x!(i,l)/x!(i,j) | k=i+1..n-1; l=0..m-1}}
+ \fBelse\fP p,i,x
+\fBwhen\fP
+ n,m = dim x; max_i, max_x = pivot i (col x j);
+ x = \fBif\fP max_x>0 \fBthen\fP swap x i max_i \fBelse\fP x;
+\fBend\fP \fBwith\fP
+ pivot i x = foldl max (0,0) [j,abs (x!j)|j=i..#x-1];
+ max (i,x) (j,y) = \fBif\fP x<y \fBthen\fP j,y \fBelse\fP i,x;
+\fBend\fP;
.fi
.PP
-In fact, the left-hand sides don't have to be related at all, so that you can
-also write something like:
+We also need the following little helper functions to swap two rows of a
+matrix (this is used in the pivoting step above) and to apply a transposition
+to a permutation (represented as a list):
.sp
.nf
-foo x | bar y = x*y;
+swap x i j = x!!(transp i j (0..n-1),0..m-1) \fBwhen\fP n,m = dim x \fBend\fP;
+transp i j p = [p!tr k | k=0..#p-1]
+\fBwith\fP tr k = \fBif\fP k==i \fBthen\fP j \fBelse\fP \fBif\fP k==j \fBthen\fP i \fBelse\fP k \fBend\fP;
.fi
.PP
-However, this is most useful when using an ``as'' pattern to bind a common
-variable to a parameter value
-.I after
-checking that it matches one of several possible argument patterns (which is
-slightly more efficient than using an equivalent type-checking guard). E.g.,
-the following definition binds the xs variable to the parameter of foo, if it
-is either the empty list or a list starting with an integer:
+Finally, let us define a convenient print representation of double matrices a
+la Octave (the meaning of the __show__ function is explained in the CAVEATS
+and NOTES section):
.sp
.nf
-foo xs@[] | foo xs@(_::int:_) = ... xs ...;
+\fBusing\fP system;
+__show__ x::matrix
+= strcat [printd j (x!(i,j))|i=0..n-1; j=0..m-1] + "\en"
+\fBwith\fP printd 0 = sprintf "\en%10.5f"; printd _ = sprintf "%10.5f" \fBend\fP
+\fBwhen\fP n,m = dim x \fBend\fP \fBif\fP dmatrixp x;
.fi
.PP
-The same construct also works in
-.B case
-expressions, which is convenient if different cases should be mapped to the
-same value, e.g.:
+Example:
.sp
.nf
-\fBcase\fP ans \fBof\fP "y" | "Y" = 1; _ = 0; \fBend\fP;
+> \fBlet\fP x = dmatrix {2,1,-1,8; -3,-1,2,-11; -2,1,2,-3};
+> x; gauss_elimination x;
+
+ 2.00000 1.00000 -1.00000 8.00000
+ -3.00000 -1.00000 2.00000 -11.00000
+ -2.00000 1.00000 2.00000 -3.00000
+
+[1,2,0],
+ 1.00000 0.33333 -0.66667 3.66667
+ 0.00000 1.00000 0.40000 2.60000
+ 0.00000 0.00000 1.00000 -1.00000
.fi
.SH MACROS
Macros are a special type of functions to be executed as a kind of
@@ -1688,18 +1876,19 @@
unsigned integers as well (if necessary, you can use a bigint to pass positive
values which are too big to fit into a machine int). Also note that when an
unsigned integer is returned by a C routine which is too big to fit into the
-corresponding signed integer type, it will become negative. In this case,
-depending on the target type, you can use the ubyte, ushort, uint and ulong
-functions provided by the prelude to convert the result back to an unsigned
-quantity.
+corresponding signed integer type, it will ``wrap around'' and become
+negative. In this case, depending on the target type, you can use the ubyte,
+ushort, uint and ulong functions provided by the prelude to convert the result
+back to an unsigned quantity.
.PP
Concerning the pointer types, char* is for string arguments and return values
which need translation between Pure's internal utf-8 representation and the
system encoding, while void* is for any generic kind of pointer (including
strings, which are \fInot\fP translated when passed/returned as void*). Any
-other kind of pointer (except expr*, see below) is effectively treated as
-void* right now, although in a future version the interpreter may keep track
-of the type names for the purpose of checking parameter types.
+other kind of pointer (except expr* and the GSL matrix pointer types, which
+are discussed below) is effectively treated as void* right now, although in a
+future version the interpreter may keep track of the type names for the
+purpose of checking parameter types.
.PP
The expr* pointer type is special; it indicates a Pure expression parameter or
return value which is just passed through unchanged. All other types of values
@@ -1708,6 +1897,32 @@
Pure). All of this is handled by the runtime system in a transparent way, of
course.
.PP
+The matrix pointer types dmatrix*, cmatrix* and imatrix* can be used to pass
+double, complex double and int matrices to GSL functions taking pointers to
+the corresponding GSL types (gsl_matrix, gsl_matrix_complex and
+gsl_matrix_int) as arguments or returning them as results. Note that there is
+no marshalling of Pure's symbolic matrix type, as these aren't supported by
+GSL anyway. Also note that matrices are always passed by reference. If you
+need to pass a matrix as an output parameter of a GSL matrix routine, you can
+either create a zero matrix or a copy of an existing matrix. The prelude
+provides various operations for that purpose (in particular, see the dmatrix,
+cmatrix, imatrix and pack functions in matrices.pure). For instance, here is
+how you can quickly wrap up GSL's double matrix addition function in a way
+that preserves value semantics:
+.sp
+.nf
+> \fBextern\fP int gsl_matrix_add(dmatrix*, dmatrix*);
+> x::matrix + y::matrix = gsl_matrix_add x y $$ x \fBwhen\fP x = pack x \fBend\fP;
+> \fBlet\fP x = dmatrix {1,2,3}; \fBlet\fP y = dmatrix {2,3,2}; x; y; x+y;
+{1.0,2.0,3.0}
+{2.0,3.0,2.0}
+{3.0,5.0,5.0}
+.fi
+.PP
+Most GSL matrix routines can be wrapped in this fashion quite easily. A
+ready-made GSL interface providing access to all of GSL's numeric functions
+will be provided in the future.
+.PP
As already mentioned, it is possible to augment an external C function with
ordinary Pure equations, but in this case you have to make sure that the
.B extern
@@ -1795,17 +2010,81 @@
.B using
clause. The prelude offers the necessary functions to work with the built-in
types (including arithmetic and logical operations) and to do most kind of
-list processing you can find in ML- and Haskell-like languages. Please refer
-to the
+list processing you can find in ML- and Haskell-like languages. It also
+provides a collection of basic string and matrix operations. Please refer to
+the
.B prelude.pure
-file for details on the provided operations. Common container data structures
-like sets and dictionaries are also available, see
-.BR set.pure ,
+file (as well as the modules included there, specifically
+.BR primitives.pure ,
+.B matrices.pure
+and
+.BR strings.pure )
+for details on the provided operations. Here is a very brief summary of some
+of the prelude operations which, besides the usual arithmetic and logical
+operators, are probably used most frequently:
+.TP
+x+y
+This is also used to denote list concatenation.
+.TP
+x:y
+This is the list-consing operation. x becomes the head of the list, y its tail.
+.TP
+x..y
+Constructs arithmetic sequences. x:y..z can be used to denote sequences with
+arbitrary stepsize y-x. Infinite sequences can be constructed using an
+infinite bound (i.e., inf or -inf). E.g., 1:3..inf denotes the stream of all
+positive odd (machine) integers.
+.TP
+#x
+The size (number of elements) of the list, tuple or matrix x. In addition, dim
+x yields the dimensions (number of rows and columns) of a matrix.
+.TP
+x'
+The transpose of a matrix.
+.TP
+x!y
+This is the list, tuple and matrix indexing operation. Note that all indices
+in Pure are zero-based, thus x!0 and x!(#x-1) are the first and last element
+of a list, tuple or matrix, respectively. In the case of matrices, the
+subscript may also be a pair of row and column indices, such as x!(1,2).
+.TP
+x!!ys
+This is Pure's list, tuple and matrix ``slicing'' operation, which returns the
+list, tuple or matrix of all x!y while y runs through the (list or matrix) ys.
+Thus, e.g., x!!(i..j) returns all the elements between i and j (inclusive).
+Indices which fall outside the valid index range are quietly discarded. In
+fact, the index range ys may contain any number of indices (also duplicates),
+in any order. Thus x![0|i=1..n] returns the first element of x n times, and,
+if ys is a permutation of the range 0..#x-1, then x!!ys yields the
+corresponding permutation of the elements of x. In the case of matrices the
+index range may also contain two-dimensional subscripts, or the index range
+itself may be specified as a pair of row/column index lists such as
+x!!(i..j,k..l).
+.PP
+The prelude also offers support operations for the implementation of list and
+matrix comprehensions, as well as the customary list operations like head,
+tail, drop, take, filter, map, foldl, foldr, scanl, scanr, zip, unzip, etc.,
+which make list programming so much fun in modern FPLs. In Pure, these also
+work on strings as well as matrices, although, for reasons of efficiency,
+these data structures are internally represented as different kinds of array
+data structures.
+.PP
+Besides the prelude, Pure's standard library also comprises a growing number
+of additional library modules which we can only mention in passing here. In
+particular, the
+.B math.pure
+module provides additional mathematical functions as well as Pure's complex
+and rational number data types. Common container data structures like sets and
+dictionaries are implemented in the
+.B set.pure
+and
.B dict.pure
-etc. Moreover, the (beginnings of a) system interface can be found in the
+modules, among others. Moreover, the (beginnings of a) system interface can be
+found in the
.B system.pure
-module. In particular, this module also includes operations to do basic
-I/O. More stuff will likely be provided in future releases.
+module. In particular, this module also provides operations to do basic
+C-style I/O, including printf and scanf. More stuff will likely be provided in
+future releases.
.SH INTERACTIVE USAGE
In interactive mode, the interpreter reads definitions and expressions and
processes them as usual. If the
@@ -2288,12 +2567,24 @@
operations with side effects (it does allow you to call any C function after
all), but with a few exceptions the standard library operations are free of
those. Just stay away from operations marked ``IMPURE'' in the library sources
-(most notably, eval and catch/throw) and avoid the system module, then your
-program will behave according to the semantics of term rewriting.
+(most notably, eval, catch/throw, references, sentries and direct pointer
+manipulations) and avoid the system module, then your program will behave
+according to the semantics of term rewriting.
.PP
The short answer is that I simply liked the name, and there wasn't any
programming language named ``Pure'' yet (quite a feat nowadays), so there's
one now. :)
+.SS Backward Compatibility
+Pure 0.7 introduced built-in matrix structures, which called for some minor
+changes in the syntax of comprehensions and arithmetic
+sequences. Specifically, the template expression and generator/filter clauses
+of a comprehension are now separated with '|'. (For the time being, the old
+[x; ...] list comprehension syntax is still supported, but the compiler will
+warn you about such constructs and flag them as deprecated.) Moreover,
+arithmetic sequences with arbitrary stepsize are now written x:y..z instead of
+x,y..z, and the `..' operator now has a lower precedence than the `,'
+operator. This makes writing matrix slices like x!!(i..j,k..l) much more
+convenient.
.SS Debugging
There's no symbolic debugger yet. So
.BR printf (3)
@@ -2377,6 +2668,24 @@
1:2:3:...
.fi
.PP
+Another case which needs special consideration are numeric matrices. For
+efficiency, the expression printer will always use the default representation
+for these, unless you override the representation of the matrix as a
+whole. E.g., the following rule for double matrices mimics Octave's default
+output format (for the sake of simplicity, this isn't perfect, but you get the
+idea):
+.sp
+.nf
+> __show__ x::matrix =
+> strcat [printd j (x!(i,j))|i=0..n-1; j=0..m-1] + "\en"
+> \fBwith\fP printd 0 = sprintf "\en%10.5f"; printd _ = sprintf "%10.5f" \fBend\fP
+> \fBwhen\fP n,m = dim x \fBend\fP \fBif\fP dmatrixp x;
+> {1.0,1/2;1/3,4.0};
+
+ 1.00000 0.50000
+ 0.33333 4.00000
+.fi
+.PP
Finally, by just purging the definition of the __show__ function you can
easily go back to the standard print syntax:
.sp
@@ -2826,15 +3135,25 @@
Albert Graef <Dr....@t-...>, Dept. of Computer Music, Johannes
Gutenberg University of Mainz, Germany.
.SH SEE ALSO
+(All software listed here is freely available, usually under the GNU Public
+License.)
.TP
.B Aardappel
Another functional programming language based on term rewriting,
\fIhttp://wouter.fov120.com/aardappel\fP.
.TP
.B Alice ML
-A version of ML (see below) with ``futures'',
-\fIhttp://www.ps.uni-sb.de/alice\fP.
+A version of ML (see below) from which Pure borrows its model of lazy
+evaluation, \fIhttp://www.ps.uni-sb.de/alice\fP.
.TP
+.B GNU Octave
+A popular high-level language for numeric applications and free MATLAB
+replacement, \fIhttp://www.gnu.org/software/octave\fP.
+.TP
+.B GNU Scientific Library
+A free software library for numeric applications, required for Pure's
+numeric matrix support, \fIhttp://www.gnu.org/software/gsl\fP.
+.TP
.B Haskell
A popular non-strict FPL, \fIhttp://www.haskell.org\fP.
.TP
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-24 13:08:22
|
Revision: 849
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=849&view=rev
Author: agraef
Date: 2008-09-24 13:08:16 +0000 (Wed, 24 Sep 2008)
Log Message:
-----------
Fix some minor glitches, comment changes.
Modified Paths:
--------------
pure/trunk/lib/matrices.pure
pure/trunk/lib/prelude.pure
Modified: pure/trunk/lib/matrices.pure
===================================================================
--- pure/trunk/lib/matrices.pure 2008-09-24 07:57:03 UTC (rev 848)
+++ pure/trunk/lib/matrices.pure 2008-09-24 13:08:16 UTC (rev 849)
@@ -130,7 +130,8 @@
nth x n = catch (cst {}) (row x n);
mth x m = catch (cst {}) (col x m);
end;
-x::matrix!!ns = if packed x then rowvector x!!([0],ns)
+x::matrix!!ns = if all intp ns && packed x
+ then rowvector x!!([0],ns)
else colcatmap (nth x) ns with
nth x n = catch (cst {}) {x!n};
end;
@@ -243,7 +244,10 @@
/* Pack a matrix. This creates a copy of the matrix which has the data in
contiguous storage. It also frees up extra memory if the matrix was created
as a slice from a bigger matrix (see 'submat' above). The 'packed'
- predicate can be used to verify whether a matrix is already packed. */
+ predicate can be used to verify whether a matrix is already packed. Note
+ that even if a matrix is already packed, 'pack' will make a copy of it
+ anyway, so this routine also provides a quick way to copy a matrix, e.g.,
+ if you want to pass it as an input/output parameter to a GSL routine. */
pack x::matrix = colcat [x,{}];
packed x::matrix = stride x==dim x!1;
Modified: pure/trunk/lib/prelude.pure
===================================================================
--- pure/trunk/lib/prelude.pure 2008-09-24 07:57:03 UTC (rev 848)
+++ pure/trunk/lib/prelude.pure 2008-09-24 13:08:16 UTC (rev 849)
@@ -243,13 +243,14 @@
stream () = [];
stream xs@(_,_) = stream (list xs);
-/* Slicing. xs!!ns returns the list of xs!n for all members n of the index
- list ns which are in the valid index range. This is a generic definition
- which will work with any kind of container data structure which defines (!)
- in such a manner that it throws an exception when the index is out of
- bounds. */
+/* Slicing. xs!!ns returns the list (or tuple) of xs!n for all members n of
+ the index list ns which are in the valid index range. This is a generic
+ definition which will work with any kind of container data structure which
+ defines (!) in such a manner that it throws an exception when the index is
+ out of bounds. */
-xs!!ns = catmap (nth xs) ns
+xs!!ns = if tuplep xs then tuple ys else ys
+ when ys = catmap (nth xs) ns end
with nth xs n = catch (cst []) [xs!n] end;
/* Arithmetic sequences. */
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-24 07:57:14
|
Revision: 848
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=848&view=rev
Author: agraef
Date: 2008-09-24 07:57:03 +0000 (Wed, 24 Sep 2008)
Log Message:
-----------
Bugfixes.
Modified Paths:
--------------
pure/trunk/lib/system.pure
Modified: pure/trunk/lib/system.pure
===================================================================
--- pure/trunk/lib/system.pure 2008-09-24 01:11:44 UTC (rev 847)
+++ pure/trunk/lib/system.pure 2008-09-24 07:57:03 UTC (rev 848)
@@ -259,6 +259,8 @@
res = pure_fprintf fp s;
count = if res>=0 then count+res else throw (printf_error res);
end;
+ do_fprintf fp (count,[]) (printf_format_spec t s) =
+ throw (printf_value_error s ());
do_fprintf fp (count,_) _ = throw (this_cant_happen count);
end;
@@ -328,6 +330,8 @@
u = if res>=0 then u + cstring buf
else free buf $$ throw (printf_error res);
end;
+ do_sprintf (u,[]) (printf_format_spec t s) =
+ throw (printf_value_error s ());
do_sprintf (u,_) _ = throw (this_cant_happen u);
check_buf buf = throw malloc_error if null buf;
= buf otherwise;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-24 01:11:51
|
Revision: 847
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=847&view=rev
Author: agraef
Date: 2008-09-24 01:11:44 +0000 (Wed, 24 Sep 2008)
Log Message:
-----------
Updated NEWS.
Modified Paths:
--------------
pure/trunk/NEWS
Modified: pure/trunk/NEWS
===================================================================
--- pure/trunk/NEWS 2008-09-24 00:46:01 UTC (rev 846)
+++ pure/trunk/NEWS 2008-09-24 01:11:44 UTC (rev 847)
@@ -72,9 +72,9 @@
pointer values, and marshalling of GSL matrices in the C interface is also
available.
-Here is a brief synopsis of some important operations which are already
-implemented in the prelude (you can find the definitions of these and the
-other matrix operations in primitives.pure):
+Here is a brief synopsis of some important operations which are implemented in
+the prelude (you can find the definitions of these and a bunch of other matrix
+operations in matrices.pure):
- #x total number of elements
- dim x number of rows and columns (as a pair)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-24 00:46:08
|
Revision: 846
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=846&view=rev
Author: agraef
Date: 2008-09-24 00:46:01 +0000 (Wed, 24 Sep 2008)
Log Message:
-----------
Naming of matrix-pointer functions, code refactoring, fix some smaller glitches.
Modified Paths:
--------------
pure/trunk/lib/matrices.pure
Modified: pure/trunk/lib/matrices.pure
===================================================================
--- pure/trunk/lib/matrices.pure 2008-09-23 23:56:53 UTC (rev 845)
+++ pure/trunk/lib/matrices.pure 2008-09-24 00:46:01 UTC (rev 846)
@@ -62,18 +62,6 @@
x::matrix != y::matrix = not x == y;
-/* Convenience functions to create zero matrices with the given dimensions
- (either a pair denoting the number of rows and columns, or just the row
- size in order to create a row vector). */
-
-dmatrix (n::int,m::int) = dmatrix_from_array_dup (pointer 0) (n,m);
-cmatrix (n::int,m::int) = cmatrix_from_array_dup (pointer 0) (n,m);
-imatrix (n::int,m::int) = imatrix_from_array_dup (pointer 0) (n,m);
-
-dmatrix n::int = dmatrix (1,n);
-cmatrix n::int = cmatrix (1,n);
-imatrix n::int = imatrix (1,n);
-
/* Size of a matrix (#x) and its dimensions (dim x=n,m where n is the number
of rows, m the number of columns). Note that Pure supports empty matrices,
thus the total size may well be zero, meaning that either the number of
@@ -118,7 +106,7 @@
when n::int,m::int = dim x end);
= throw out_of_bounds otherwise;
-/* Matrix slices (x!!ns). As with simple indexing, elements can addressed
+/* Matrix slices (x!!ns). As with simple indexing, elements can be addressed
using either singleton (row-major) indices or index pair (row,column). The
former is specified as a list of int values, the latter as a pair of lists
of int values. As with list slicing, index ranges may be non-contiguous
@@ -197,26 +185,6 @@
submat x::matrix (i::int,j::int) (n::int,m::int)
= matrix_slice x i j (i+n-1) (j+m-1);
-/* Change the dimensions of a matrix without changing its size. The total
- number of elements must match that of the input matrix. Reuses the
- underlying storage of the input matrix if possible. */
-
-private matrix_redim;
-extern expr* matrix_redim(expr* x, int n, int m);
-
-redim x::matrix (n::int,m::int)
- = matrix_redim x n m if n>=0 && m>=0 && n*m==#x;
-
-/* You can also redim a matrix to a given row size. In this case the row size
- must be positive and divide the total size of the matrix, */
-
-redim x::matrix m::int = redim x (#x div m,m) if m>0 && #x mod m==0;
-
-/* Convenience functions for converting a matrix to a row or column vector. */
-
-rowvector x::matrix = redim x (#x);
-colvector x::matrix = redim x 1;
-
/* Construct matrices from lists of rows and columns. These take either
scalars or submatrices as inputs; corresponding dimensions must match.
rowcat combines submatrices vertically, like {x;y}; colcat combines them
@@ -242,6 +210,65 @@
def void (rowcatmap f x) = do f x;
def void (colcatmap f x) = do f x;
+/* Convenience functions to create zero matrices with the given dimensions
+ (either a pair denoting the number of rows and columns, or just the row
+ size in order to create a row vector). */
+
+dmatrix (n::int,m::int) = dmatrix_dup (pointer 0,n,m);
+cmatrix (n::int,m::int) = cmatrix_dup (pointer 0,n,m);
+imatrix (n::int,m::int) = imatrix_dup (pointer 0,n,m);
+
+dmatrix n::int = dmatrix (1,n);
+cmatrix n::int = cmatrix (1,n);
+imatrix n::int = imatrix (1,n);
+
+/* Matrix conversions. These convert between different types of numeric
+ matrices. You can also extract the real and imaginary parts of a (complex)
+ matrix. */
+
+private matrix_double matrix_complex matrix_int;
+extern expr* matrix_double(expr *x), expr* matrix_complex(expr *x),
+ expr* matrix_int(expr *x);
+
+dmatrix x::matrix = matrix_double x if imatrixp x || dmatrixp x;
+imatrix x::matrix = matrix_int x if imatrixp x || dmatrixp x;
+cmatrix x::matrix = matrix_complex x if nmatrixp x;
+
+private matrix_re matrix_im;
+extern expr* matrix_re(expr *x), expr* matrix_im(expr *x);
+
+re x::matrix = matrix_re x if nmatrixp x;
+im x::matrix = matrix_im x if nmatrixp x;
+
+/* Pack a matrix. This creates a copy of the matrix which has the data in
+ contiguous storage. It also frees up extra memory if the matrix was created
+ as a slice from a bigger matrix (see 'submat' above). The 'packed'
+ predicate can be used to verify whether a matrix is already packed. */
+
+pack x::matrix = colcat [x,{}];
+packed x::matrix = stride x==dim x!1;
+
+/* Change the dimensions of a matrix without changing its size. The total
+ number of elements must match that of the input matrix. Reuses the
+ underlying storage of the input matrix if possible (i.e., if the matrix is
+ packed). */
+
+private matrix_redim;
+extern expr* matrix_redim(expr* x, int n, int m);
+
+redim x::matrix (n::int,m::int)
+ = matrix_redim x n m if n>=0 && m>=0 && n*m==#x;
+
+/* You can also redim a matrix to a given row size. In this case the row size
+ must be positive and divide the total size of the matrix, */
+
+redim x::matrix m::int = redim x (#x div m,m) if m>0 && #x mod m==0;
+
+/* Convert a matrix to a row or column vector. */
+
+rowvector x::matrix = redim x (1,#x);
+colvector x::matrix = redim x 1;
+
/* Transpose a matrix. */
private matrix_transpose;
@@ -328,43 +355,16 @@
flip redim (dim x) (colcat w)
when u,v,w = unzip3 (list x) end;
-/* Matrix conversions. These convert between different types of numeric
- matrices. You can also extract the real and imaginary parts of a (complex)
- matrix. */
-
-private matrix_double matrix_complex matrix_int;
-extern expr* matrix_double(expr *x), expr* matrix_complex(expr *x),
- expr* matrix_int(expr *x);
-
-dmatrix x::matrix = matrix_double x if imatrixp x || dmatrixp x;
-imatrix x::matrix = matrix_int x if imatrixp x || dmatrixp x;
-cmatrix x::matrix = matrix_complex x if nmatrixp x;
-
-private matrix_re matrix_im;
-extern expr* matrix_re(expr *x), expr* matrix_im(expr *x);
-
-re x::matrix = matrix_re x if nmatrixp x;
-im x::matrix = matrix_im x if nmatrixp x;
-
-/* 'Pack' a matrix. This creates a copy of the matrix which has the data in
- contiguous storage. It also frees up extra memory if the matrix was created
- as a slice from a bigger matrix. The 'packed' predicate can be used to
- verify whether a matrix is already packed. */
-
-pack x::matrix = colcat [x,{}];
-packed x::matrix = stride x==dim x!1;
-
/* Low-level operations for converting between matrices and raw pointers.
These are typically used to shovel around massive amounts of numeric data
between Pure and external C routines, when performance and throughput is an
important consideration (e.g., graphics, video and audio applications). The
usual caveats apply. */
-/* Convert a matrix to a pointer. Returns a pointer pointing directly to the
- underlying C array. You have to be careful when passing such a pointer to C
- functions if the underlying data is non-contiguous; when in doubt, first
- use the 'pack' function from above to place the data in contiguous
- storage. */
+/* Convert a matrix to a pointer, which points directly to the underlying C
+ array. You have to be careful when passing such a pointer to C functions if
+ the underlying data is non-contiguous; when in doubt, first use the 'pack'
+ function from above to place the data in contiguous storage. */
private pure_pointerval;
extern expr* pure_pointerval(expr*);
@@ -383,19 +383,19 @@
extern expr* matrix_from_complex_array(void* p, int n, int m);
extern expr* matrix_from_int_array(void* p, int n, int m);
-dmatrix_from_array p::pointer (n::int,m::int)
+dmatrix (p::pointer,n::int,m::int)
= matrix_from_double_array p n m;
-cmatrix_from_array p::pointer (n::int,m::int)
+cmatrix (p::pointer,n::int,m::int)
= matrix_from_complex_array p n m;
-imatrix_from_array p::pointer (n::int,m::int)
+imatrix (p::pointer,n::int,m::int)
= matrix_from_int_array p n m;
-dmatrix_from_array p::pointer n::int
- = dmatrix_from_array p (1,n);
-cmatrix_from_array p::pointer n::int
- = cmatrix_from_array p (1,n);
-imatrix_from_array p::pointer n::int
- = imatrix_from_array p (1,n);
+dmatrix (p::pointer,n::int)
+ = dmatrix (p,1,n);
+cmatrix (p::pointer,n::int)
+ = cmatrix (p,1,n);
+imatrix (p::pointer,n::int)
+ = imatrix (p,1,n);
/* Create a numeric matrix from a pointer, copying the data to fresh memory.
The source pointer p may also be NULL, in which case the new matrix is
@@ -408,23 +408,23 @@
extern expr* matrix_from_complex_array_dup(void* p, int n, int m);
extern expr* matrix_from_int_array_dup(void* p, int n, int m);
-dmatrix_from_array_dup p::pointer (n::int,m::int)
+dmatrix_dup (p::pointer,n::int,m::int)
= matrix_from_double_array_dup p n m;
-cmatrix_from_array_dup p::pointer (n::int,m::int)
+cmatrix_dup (p::pointer,n::int,m::int)
= matrix_from_complex_array_dup p n m;
-imatrix_from_array_dup p::pointer (n::int,m::int)
+imatrix_dup (p::pointer,n::int,m::int)
= matrix_from_int_array_dup p n m;
-dmatrix_from_array_dup p::pointer n::int
- = dmatrix_from_array_dup p (1,n);
-cmatrix_from_array_dup p::pointer n::int
- = cmatrix_from_array_dup p (1,n);
-imatrix_from_array_dup p::pointer n::int
- = imatrix_from_array_dup p (1,n);
+dmatrix_dup (p::pointer,n::int)
+ = dmatrix_dup (p,1,n);
+cmatrix_dup (p::pointer,n::int)
+ = cmatrix_dup (p,1,n);
+imatrix_dup (p::pointer,n::int)
+ = imatrix_dup (p,1,n);
-/* Some helper routines for dealing with alternate base types. These work like
- the routines above, but initialize the data from float, complex float,
- short and byte arrays, respectively. */
+/* Some additional functions for alternate base types. These work like the
+ routines above, but initialize the data from float, complex float, short
+ and byte arrays, respectively. */
private matrix_from_float_array_dup;
private matrix_from_complex_float_array_dup;
@@ -435,20 +435,20 @@
extern expr* matrix_from_short_array_dup(void* p, int n, int m);
extern expr* matrix_from_byte_array_dup(void* p, int n, int m);
-fmatrix_from_array_dup p::pointer (n::int,m::int)
+float_dmatrix_dup (p::pointer,n::int,m::int)
= matrix_from_float_array_dup p n m;
-cfmatrix_from_array_dup p::pointer (n::int,m::int)
+float_cmatrix_dup (p::pointer,n::int,m::int)
= matrix_from_complex_float_array_dup p n m;
-smatrix_from_array_dup p::pointer (n::int,m::int)
+short_imatrix_dup (p::pointer,n::int,m::int)
= matrix_from_short_array_dup p n m;
-bmatrix_from_array_dup p::pointer (n::int,m::int)
+byte_imatrix_dup (p::pointer,n::int,m::int)
= matrix_from_byte_array_dup p n m;
-fmatrix_from_array_dup p::pointer n::int
- = fmatrix_from_array_dup p (1,n);
-cfmatrix_from_array_dup p::pointer n::int
- = cfmatrix_from_array_dup p (1,n);
-smatrix_from_array_dup p::pointer n::int
- = smatrix_from_array_dup p (1,n);
-bmatrix_from_array_dup p::pointer n::int
- = bmatrix_from_array_dup p (1,n);
+float_dmatrix_dup (p::pointer,n::int)
+ = float_dmatrix_dup (p,1,n);
+float_cmatrix_dup (p::pointer,n::int)
+ = float_cmatrix_dup (p,1,n);
+short_imatrix_dup (p::pointer,n::int)
+ = short_imatrix_dup (p,1,n);
+byte_imatrix_dup (p::pointer,n::int)
+ = byte_imatrix_dup (p,1,n);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-23 23:57:04
|
Revision: 845
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=845&view=rev
Author: agraef
Date: 2008-09-23 23:56:53 +0000 (Tue, 23 Sep 2008)
Log Message:
-----------
Add some more matrix-pointer operations for alternative base types.
Modified Paths:
--------------
pure/trunk/lib/matrices.pure
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/lib/matrices.pure
===================================================================
--- pure/trunk/lib/matrices.pure 2008-09-23 23:25:21 UTC (rev 844)
+++ pure/trunk/lib/matrices.pure 2008-09-23 23:56:53 UTC (rev 845)
@@ -421,3 +421,34 @@
= cmatrix_from_array_dup p (1,n);
imatrix_from_array_dup p::pointer n::int
= imatrix_from_array_dup p (1,n);
+
+/* Some helper routines for dealing with alternate base types. These work like
+ the routines above, but initialize the data from float, complex float,
+ short and byte arrays, respectively. */
+
+private matrix_from_float_array_dup;
+private matrix_from_complex_float_array_dup;
+private matrix_from_short_array_dup;
+private matrix_from_byte_array_dup;
+extern expr* matrix_from_float_array_dup(void* p, int n, int m);
+extern expr* matrix_from_complex_float_array_dup(void* p, int n, int m);
+extern expr* matrix_from_short_array_dup(void* p, int n, int m);
+extern expr* matrix_from_byte_array_dup(void* p, int n, int m);
+
+fmatrix_from_array_dup p::pointer (n::int,m::int)
+ = matrix_from_float_array_dup p n m;
+cfmatrix_from_array_dup p::pointer (n::int,m::int)
+ = matrix_from_complex_float_array_dup p n m;
+smatrix_from_array_dup p::pointer (n::int,m::int)
+ = matrix_from_short_array_dup p n m;
+bmatrix_from_array_dup p::pointer (n::int,m::int)
+ = matrix_from_byte_array_dup p n m;
+
+fmatrix_from_array_dup p::pointer n::int
+ = fmatrix_from_array_dup p (1,n);
+cfmatrix_from_array_dup p::pointer n::int
+ = cfmatrix_from_array_dup p (1,n);
+smatrix_from_array_dup p::pointer n::int
+ = smatrix_from_array_dup p (1,n);
+bmatrix_from_array_dup p::pointer n::int
+ = bmatrix_from_array_dup p (1,n);
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-23 23:25:21 UTC (rev 844)
+++ pure/trunk/runtime.cc 2008-09-23 23:56:53 UTC (rev 845)
@@ -5191,6 +5191,152 @@
#endif
}
+extern "C"
+pure_expr *matrix_from_float_array_dup(void *p, uint32_t n1, uint32_t n2)
+{
+#ifdef HAVE_GSL
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_double_matrix(create_double_matrix(n1, n2));
+ if (!p)
+ p = calloc(n1*n2, sizeof(double));
+ else {
+ void *q = malloc(n1*n2*sizeof(double));
+ float *p1 = (float*)p; double *q1 = (double*)q;
+ for (size_t i = 0; i < n1*n2; i++) q1[i] = (double)p1[i];
+ p = q;
+ }
+ if (!p) return 0;
+ gsl_matrix_view v = gsl_matrix_view_array((double*)p, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix *m = (gsl_matrix*)malloc(sizeof(gsl_matrix));
+ gsl_block *b = (gsl_block*)malloc(sizeof(gsl_block));
+ assert(m && b && v.matrix.data);
+ *m = v.matrix;
+ b->size = n1*n2;
+ b->data = m->data;
+ m->block = b;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::DMATRIX;
+ x->data.mat.p = m;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_from_complex_float_array_dup(void *p, uint32_t n1, uint32_t n2)
+{
+#ifdef HAVE_GSL
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_complex_matrix(create_complex_matrix(n1, n2));
+ if (!p)
+ p = calloc(2*n1*n2, sizeof(double));
+ else {
+ void *q = malloc(2*n1*n2*sizeof(double));
+ float *p1 = (float*)p; double *q1 = (double*)q;
+ for (size_t i = 0; i < 2*n1*n2; i++) q1[i] = (double)p1[i];
+ p = q;
+ }
+ if (!p) return 0;
+ gsl_matrix_complex_view v =
+ gsl_matrix_complex_view_array((double*)p, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix_complex *m =
+ (gsl_matrix_complex*)malloc(sizeof(gsl_matrix_complex));
+ gsl_block_complex *b = (gsl_block_complex*)malloc(sizeof(gsl_block_complex));
+ assert(m && b && v.matrix.data);
+ *m = v.matrix;
+ b->size = n1*n2;
+ b->data = m->data;
+ m->block = b;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::CMATRIX;
+ x->data.mat.p = m;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_from_short_array_dup(void *p, uint32_t n1, uint32_t n2)
+{
+#ifdef HAVE_GSL
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_int_matrix(create_int_matrix(n1, n2));
+ if (!p)
+ p = calloc(n1*n2, sizeof(int));
+ else {
+ void *q = malloc(n1*n2*sizeof(int));
+ short *p1 = (short*)p; int *q1 = (int*)q;
+ for (size_t i = 0; i < n1*n2; i++) q1[i] = (int)p1[i];
+ p = q;
+ }
+ if (!p) return 0;
+ gsl_matrix_int_view v = gsl_matrix_int_view_array((int*)p, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix_int *m = (gsl_matrix_int*)malloc(sizeof(gsl_matrix_int));
+ gsl_block_int *b = (gsl_block_int*)malloc(sizeof(gsl_block_int));
+ assert(m && b && v.matrix.data);
+ *m = v.matrix;
+ b->size = n1*n2;
+ b->data = m->data;
+ m->block = b;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::IMATRIX;
+ x->data.mat.p = m;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_from_byte_array_dup(void *p, uint32_t n1, uint32_t n2)
+{
+#ifdef HAVE_GSL
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_int_matrix(create_int_matrix(n1, n2));
+ if (!p)
+ p = calloc(n1*n2, sizeof(int));
+ else {
+ void *q = malloc(n1*n2*sizeof(int));
+ int8_t *p1 = (int8_t*)p; int *q1 = (int*)q;
+ for (size_t i = 0; i < n1*n2; i++) q1[i] = (int)p1[i];
+ p = q;
+ }
+ if (!p) return 0;
+ gsl_matrix_int_view v = gsl_matrix_int_view_array((int*)p, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix_int *m = (gsl_matrix_int*)malloc(sizeof(gsl_matrix_int));
+ gsl_block_int *b = (gsl_block_int*)malloc(sizeof(gsl_block_int));
+ assert(m && b && v.matrix.data);
+ *m = v.matrix;
+ b->size = n1*n2;
+ b->data = m->data;
+ m->block = b;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::IMATRIX;
+ x->data.mat.p = m;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
+#else
+ return 0;
+#endif
+}
+
static uint32_t mpz_hash(const mpz_t z)
{
uint32_t h = 0;
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-23 23:25:21 UTC (rev 844)
+++ pure/trunk/runtime.h 2008-09-23 23:56:53 UTC (rev 845)
@@ -746,6 +746,7 @@
/* The following routines work like the above, but copy the data to newly
allocated memory, so the original data can be freed after the call.
+ Additional routines are provided for various alternate data sizes.
Moreover, the source pointer p may also be NULL in which case the new
matrix is filled with zeros instead. */
@@ -753,6 +754,11 @@
pure_expr *matrix_from_complex_array_dup(void *p, uint32_t n, uint32_t m);
pure_expr *matrix_from_int_array_dup(void *p, uint32_t n, uint32_t m);
+pure_expr *matrix_from_float_array_dup(void *p, uint32_t n, uint32_t m);
+pure_expr *matrix_from_complex_float_array_dup(void *p, uint32_t n, uint32_t m);
+pure_expr *matrix_from_short_array_dup(void *p, uint32_t n, uint32_t m);
+pure_expr *matrix_from_byte_array_dup(void *p, uint32_t n, uint32_t m);
+
/* Compute a 32 bit hash code of a Pure expression. This makes it possible to
use arbitary Pure values as keys in a hash table. */
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-23 23:25:32
|
Revision: 844
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=844&view=rev
Author: agraef
Date: 2008-09-23 23:25:21 +0000 (Tue, 23 Sep 2008)
Log Message:
-----------
Overhaul of matrix-pointer operations.
Modified Paths:
--------------
pure/trunk/lib/matrices.pure
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/lib/matrices.pure
===================================================================
--- pure/trunk/lib/matrices.pure 2008-09-23 22:03:18 UTC (rev 843)
+++ pure/trunk/lib/matrices.pure 2008-09-23 23:25:21 UTC (rev 844)
@@ -43,17 +43,36 @@
rowvectorp x = matrixp x && dim x!0==1;
colvectorp x = matrixp x && dim x!1==1;
-/* Convenience functions to create numeric matrices with the given dimensions
+/* Matrix comparisons. */
+
+x::matrix == y::matrix = x === y
+ if nmatrixp x && matrix_type x == matrix_type y;
+// mixed numeric cases
+ = cmatrix x === y if nmatrixp x && cmatrixp y;
+ = x === cmatrix y if cmatrixp x && nmatrixp y;
+ = dmatrix x === y if imatrixp x && dmatrixp y;
+ = x === dmatrix y if dmatrixp x && imatrixp y;
+// comparisons with symbolic matrices
+ = 0 if dim x != dim y;
+ = compare 0 with
+ compare i::int = 1 if i>=n;
+ = 0 if x!i != y!i;
+ = compare (i+1);
+ end when n::int = #x end;
+
+x::matrix != y::matrix = not x == y;
+
+/* Convenience functions to create zero matrices with the given dimensions
(either a pair denoting the number of rows and columns, or just the row
- size in order to create a row vector), and all elements zero. */
+ size in order to create a row vector). */
-dmatrix (n::int,m::int) = cdmatrix (pointer 0) (n,m);
-cmatrix (n::int,m::int) = ccmatrix (pointer 0) (n,m);
-imatrix (n::int,m::int) = cimatrix (pointer 0) (n,m);
+dmatrix (n::int,m::int) = dmatrix_from_array_dup (pointer 0) (n,m);
+cmatrix (n::int,m::int) = cmatrix_from_array_dup (pointer 0) (n,m);
+imatrix (n::int,m::int) = imatrix_from_array_dup (pointer 0) (n,m);
-dmatrix n::int = cdmatrix (pointer 0) n;
-cmatrix n::int = ccmatrix (pointer 0) n;
-imatrix n::int = cimatrix (pointer 0) n;
+dmatrix n::int = dmatrix (1,n);
+cmatrix n::int = cmatrix (1,n);
+imatrix n::int = imatrix (1,n);
/* Size of a matrix (#x) and its dimensions (dim x=n,m where n is the number
of rows, m the number of columns). Note that Pure supports empty matrices,
@@ -328,34 +347,34 @@
im x::matrix = matrix_im x if nmatrixp x;
/* 'Pack' a matrix. This creates a copy of the matrix which has the data in
- contiguous storage. If possible, it also frees up extra memory if the
- matrix was created as a slice from a bigger matrix. The 'packed' predicate
- can be used to verify whether a matrix is already packed. */
+ contiguous storage. It also frees up extra memory if the matrix was created
+ as a slice from a bigger matrix. The 'packed' predicate can be used to
+ verify whether a matrix is already packed. */
pack x::matrix = colcat [x,{}];
packed x::matrix = stride x==dim x!1;
-/* Convert a matrix to a pointer. This returns a pointer to the underlying C
- array, which allows you to change the data in-place by shelling out to C
- functions, so beware. You also need to be careful when passing such a
- pointer to C functions if the underlying data is non-contiguous; when in
- doubt use the 'pack' function from above to place the data in contiguous
- storage first. */
+/* Low-level operations for converting between matrices and raw pointers.
+ These are typically used to shovel around massive amounts of numeric data
+ between Pure and external C routines, when performance and throughput is an
+ important consideration (e.g., graphics, video and audio applications). The
+ usual caveats apply. */
+/* Convert a matrix to a pointer. Returns a pointer pointing directly to the
+ underlying C array. You have to be careful when passing such a pointer to C
+ functions if the underlying data is non-contiguous; when in doubt, first
+ use the 'pack' function from above to place the data in contiguous
+ storage. */
+
private pure_pointerval;
extern expr* pure_pointerval(expr*);
pointer x::matrix = pure_pointerval x;
-/* Create a numeric matrix from a pointer. The pointer to be converted must
- point to a malloc'ed, properly initialized memory area big enough to
- accommodate the requested dimensions. The pointer may also be NULL in which
- case a suitable pointer is allocated automatically. This memory is taken
- over by Pure and will be freed automatically when the matrix object is
- garbage-collected. CAVEAT: These are low-level operations and hence should
- be used with care. They are useful, in particular, if you need to handle
- massive amounts of matrix or vector data generated or processed by external
- C functions, such as routines dealing with graphics and sound data. */
+/* Create a numeric matrix from a pointer, without copying the data. The
+ caller must ensure that the pointer points to properly initialized memory
+ big enough to accommodate the requested dimensions, which persists for the
+ entire lifetime of the matrix object. */
private matrix_from_double_array;
private matrix_from_complex_array;
@@ -364,35 +383,41 @@
extern expr* matrix_from_complex_array(void* p, int n, int m);
extern expr* matrix_from_int_array(void* p, int n, int m);
-cdmatrix p::pointer (n::int,m::int)
+dmatrix_from_array p::pointer (n::int,m::int)
= matrix_from_double_array p n m;
-ccmatrix p::pointer (n::int,m::int)
+cmatrix_from_array p::pointer (n::int,m::int)
= matrix_from_complex_array p n m;
-cimatrix p::pointer (n::int,m::int)
+imatrix_from_array p::pointer (n::int,m::int)
= matrix_from_int_array p n m;
-cdmatrix p::pointer n::int
- = cdmatrix p (1,n);
-ccmatrix p::pointer n::int
- = ccmatrix p (1,n);
-cimatrix p::pointer n::int
- = cimatrix p (1,n);
+dmatrix_from_array p::pointer n::int
+ = dmatrix_from_array p (1,n);
+cmatrix_from_array p::pointer n::int
+ = cmatrix_from_array p (1,n);
+imatrix_from_array p::pointer n::int
+ = imatrix_from_array p (1,n);
-/* Matrix comparisons. */
+/* Create a numeric matrix from a pointer, copying the data to fresh memory.
+ The source pointer p may also be NULL, in which case the new matrix is
+ filled with zeros instead. */
-x::matrix == y::matrix = x === y
- if nmatrixp x && matrix_type x == matrix_type y;
-// mixed numeric cases
- = cmatrix x === y if nmatrixp x && cmatrixp y;
- = x === cmatrix y if cmatrixp x && nmatrixp y;
- = dmatrix x === y if imatrixp x && dmatrixp y;
- = x === dmatrix y if dmatrixp x && imatrixp y;
-// comparisons with symbolic matrices
- = 0 if dim x != dim y;
- = compare 0 with
- compare i::int = 1 if i>=n;
- = 0 if x!i != y!i;
- = compare (i+1);
- end when n::int = #x end;
+private matrix_from_double_array_dup;
+private matrix_from_complex_array_dup;
+private matrix_from_int_array_dup;
+extern expr* matrix_from_double_array_dup(void* p, int n, int m);
+extern expr* matrix_from_complex_array_dup(void* p, int n, int m);
+extern expr* matrix_from_int_array_dup(void* p, int n, int m);
-x::matrix != y::matrix = not x == y;
+dmatrix_from_array_dup p::pointer (n::int,m::int)
+ = matrix_from_double_array_dup p n m;
+cmatrix_from_array_dup p::pointer (n::int,m::int)
+ = matrix_from_complex_array_dup p n m;
+imatrix_from_array_dup p::pointer (n::int,m::int)
+ = matrix_from_int_array_dup p n m;
+
+dmatrix_from_array_dup p::pointer n::int
+ = dmatrix_from_array_dup p (1,n);
+cmatrix_from_array_dup p::pointer n::int
+ = cmatrix_from_array_dup p (1,n);
+imatrix_from_array_dup p::pointer n::int
+ = imatrix_from_array_dup p (1,n);
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-23 22:03:18 UTC (rev 843)
+++ pure/trunk/runtime.cc 2008-09-23 23:25:21 UTC (rev 844)
@@ -388,19 +388,19 @@
#ifdef HAVE_GSL
case EXPR::DMATRIX: {
gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
- m->owner = owner;
+ m->owner = owner && m->block;
gsl_matrix_free(m);
break;
}
case EXPR::CMATRIX: {
gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
- m->owner = owner;
+ m->owner = owner && m->block;
gsl_matrix_complex_free(m);
break;
}
case EXPR::IMATRIX: {
gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
- m->owner = owner;
+ m->owner = owner && m->block;
gsl_matrix_int_free(m);
break;
}
@@ -5016,14 +5016,94 @@
#ifdef HAVE_GSL
if (n1 == 0 || n2 == 0) // empty matrix
return pure_double_matrix(create_double_matrix(n1, n2));
- if (!p) p = calloc(n1*n2, sizeof(double));
if (!p) return 0;
gsl_matrix_view v = gsl_matrix_view_array((double*)p, n1, n2);
// take a copy of the view matrix
gsl_matrix *m = (gsl_matrix*)malloc(sizeof(gsl_matrix));
- gsl_block *b = (gsl_block*)malloc(sizeof(gsl_block));
assert(m && v.matrix.data);
*m = v.matrix;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::DMATRIX;
+ x->data.mat.p = m;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_from_complex_array(void *p, uint32_t n1, uint32_t n2)
+{
+#ifdef HAVE_GSL
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_complex_matrix(create_complex_matrix(n1, n2));
+ if (!p) return 0;
+ gsl_matrix_complex_view v =
+ gsl_matrix_complex_view_array((double*)p, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix_complex *m =
+ (gsl_matrix_complex*)malloc(sizeof(gsl_matrix_complex));
+ assert(m && v.matrix.data);
+ *m = v.matrix;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::CMATRIX;
+ x->data.mat.p = m;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_from_int_array(void *p, uint32_t n1, uint32_t n2)
+{
+#ifdef HAVE_GSL
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_int_matrix(create_int_matrix(n1, n2));
+ if (!p) return 0;
+ gsl_matrix_int_view v = gsl_matrix_int_view_array((int*)p, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix_int *m = (gsl_matrix_int*)malloc(sizeof(gsl_matrix_int));
+ assert(m && v.matrix.data);
+ *m = v.matrix;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::IMATRIX;
+ x->data.mat.p = m;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_from_double_array_dup(void *p, uint32_t n1, uint32_t n2)
+{
+#ifdef HAVE_GSL
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_double_matrix(create_double_matrix(n1, n2));
+ if (!p)
+ p = calloc(n1*n2, sizeof(double));
+ else {
+ void *q = malloc(n1*n2*sizeof(double));
+ memcpy(q, p, n1*n2*sizeof(double));
+ p = q;
+ }
+ if (!p) return 0;
+ gsl_matrix_view v = gsl_matrix_view_array((double*)p, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix *m = (gsl_matrix*)malloc(sizeof(gsl_matrix));
+ gsl_block *b = (gsl_block*)malloc(sizeof(gsl_block));
+ assert(m && b && v.matrix.data);
+ *m = v.matrix;
b->size = n1*n2;
b->data = m->data;
m->block = b;
@@ -5040,12 +5120,18 @@
}
extern "C"
-pure_expr *matrix_from_complex_array(void *p, uint32_t n1, uint32_t n2)
+pure_expr *matrix_from_complex_array_dup(void *p, uint32_t n1, uint32_t n2)
{
#ifdef HAVE_GSL
if (n1 == 0 || n2 == 0) // empty matrix
return pure_complex_matrix(create_complex_matrix(n1, n2));
- if (!p) p = calloc(2*n1*n2, sizeof(double));
+ if (!p)
+ p = calloc(2*n1*n2, sizeof(double));
+ else {
+ void *q = malloc(2*n1*n2*sizeof(double));
+ memcpy(q, p, 2*n1*n2*sizeof(double));
+ p = q;
+ }
if (!p) return 0;
gsl_matrix_complex_view v =
gsl_matrix_complex_view_array((double*)p, n1, n2);
@@ -5053,7 +5139,7 @@
gsl_matrix_complex *m =
(gsl_matrix_complex*)malloc(sizeof(gsl_matrix_complex));
gsl_block_complex *b = (gsl_block_complex*)malloc(sizeof(gsl_block_complex));
- assert(m && v.matrix.data);
+ assert(m && b && v.matrix.data);
*m = v.matrix;
b->size = n1*n2;
b->data = m->data;
@@ -5071,18 +5157,24 @@
}
extern "C"
-pure_expr *matrix_from_int_array(void *p, uint32_t n1, uint32_t n2)
+pure_expr *matrix_from_int_array_dup(void *p, uint32_t n1, uint32_t n2)
{
#ifdef HAVE_GSL
if (n1 == 0 || n2 == 0) // empty matrix
return pure_int_matrix(create_int_matrix(n1, n2));
- if (!p) p = calloc(n1*n2, sizeof(int));
+ if (!p)
+ p = calloc(n1*n2, sizeof(int));
+ else {
+ void *q = malloc(n1*n2*sizeof(int));
+ memcpy(q, p, n1*n2*sizeof(int));
+ p = q;
+ }
if (!p) return 0;
gsl_matrix_int_view v = gsl_matrix_int_view_array((int*)p, n1, n2);
// take a copy of the view matrix
gsl_matrix_int *m = (gsl_matrix_int*)malloc(sizeof(gsl_matrix_int));
gsl_block_int *b = (gsl_block_int*)malloc(sizeof(gsl_block_int));
- assert(m && v.matrix.data);
+ assert(m && b && v.matrix.data);
*m = v.matrix;
b->size = n1*n2;
b->data = m->data;
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-23 22:03:18 UTC (rev 843)
+++ pure/trunk/runtime.h 2008-09-23 23:25:21 UTC (rev 844)
@@ -734,17 +734,25 @@
pure_expr *matrix_im(pure_expr *x);
/* Create a matrix object of the given dimensions which uses the given pointer
- p as its underlying C array. There are no checks whatsoever, so the caller
- is responsible for making sure that the memory has the right size and is
- properly initialized. p must point to a malloc'ed memory area, which is
- taken over by Pure and will be freed automatically when the matrix is
- garbage-collected. p may also be NULL in which case a suitable pointer is
- allocated automatically. */
+ p as its underlying storage. There are no checks whatsoever and the data is
+ *not* copied, so the caller is responsible for making sure that the memory
+ has the right size, is properly initialized and is not freed while the
+ matrix is still in use. The memory is *not* freed when the matrix is
+ garbage-collected but remains in the ownership of the caller. */
pure_expr *matrix_from_double_array(void *p, uint32_t n, uint32_t m);
pure_expr *matrix_from_complex_array(void *p, uint32_t n, uint32_t m);
pure_expr *matrix_from_int_array(void *p, uint32_t n, uint32_t m);
+/* The following routines work like the above, but copy the data to newly
+ allocated memory, so the original data can be freed after the call.
+ Moreover, the source pointer p may also be NULL in which case the new
+ matrix is filled with zeros instead. */
+
+pure_expr *matrix_from_double_array_dup(void *p, uint32_t n, uint32_t m);
+pure_expr *matrix_from_complex_array_dup(void *p, uint32_t n, uint32_t m);
+pure_expr *matrix_from_int_array_dup(void *p, uint32_t n, uint32_t m);
+
/* Compute a 32 bit hash code of a Pure expression. This makes it possible to
use arbitary Pure values as keys in a hash table. */
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-23 22:03:28
|
Revision: 843
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=843&view=rev
Author: agraef
Date: 2008-09-23 22:03:18 +0000 (Tue, 23 Sep 2008)
Log Message:
-----------
Bugfixes.
Modified Paths:
--------------
pure/trunk/lib/matrices.pure
Modified: pure/trunk/lib/matrices.pure
===================================================================
--- pure/trunk/lib/matrices.pure 2008-09-23 21:56:51 UTC (rev 842)
+++ pure/trunk/lib/matrices.pure 2008-09-23 22:03:18 UTC (rev 843)
@@ -342,6 +342,9 @@
doubt use the 'pack' function from above to place the data in contiguous
storage first. */
+private pure_pointerval;
+extern expr* pure_pointerval(expr*);
+
pointer x::matrix = pure_pointerval x;
/* Create a numeric matrix from a pointer. The pointer to be converted must
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-23 21:56:57
|
Revision: 842
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=842&view=rev
Author: agraef
Date: 2008-09-23 21:56:51 +0000 (Tue, 23 Sep 2008)
Log Message:
-----------
Add some missing memory operations.
Modified Paths:
--------------
pure/trunk/lib/primitives.pure
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-23 18:56:29 UTC (rev 841)
+++ pure/trunk/lib/primitives.pure 2008-09-23 21:56:51 UTC (rev 842)
@@ -396,30 +396,40 @@
/* Direct memory accesses. Use with care ... or else! */
-private pointer_get_byte pointer_get_int pointer_get_double
+private pointer_get_byte pointer_get_short pointer_get_int
+ pointer_get_float pointer_get_double
pointer_get_string pointer_get_pointer;
extern int pointer_get_byte(void *ptr);
+extern int pointer_get_short(void *ptr);
extern int pointer_get_int(void *ptr);
+extern double pointer_get_float(void *ptr);
extern double pointer_get_double(void *ptr);
extern char *pointer_get_string(void *ptr);
extern void *pointer_get_pointer(void *ptr);
get_byte x::pointer = pointer_get_byte x;
+get_short x::pointer = pointer_get_short x;
get_int x::pointer = pointer_get_int x;
+get_float x::pointer = pointer_get_float x;
get_double x::pointer = pointer_get_double x;
get_string x::pointer = pointer_get_string x;
get_pointer x::pointer = pointer_get_pointer x;
-private pointer_put_byte pointer_put_int pointer_put_double
+private pointer_put_byte pointer_put_short pointer_put_int
+ pointer_put_float pointer_put_double
pointer_put_string pointer_put_pointer;
extern void pointer_put_byte(void *ptr, int x); // IMPURE!
+extern void pointer_put_short(void *ptr, int x); // IMPURE!
extern void pointer_put_int(void *ptr, int x); // IMPURE!
+extern void pointer_put_float(void *ptr, double x); // IMPURE!
extern void pointer_put_double(void *ptr, double x); // IMPURE!
extern void pointer_put_string(void *ptr, char *x); // IMPURE!
extern void pointer_put_pointer(void *ptr, void *x); // IMPURE!
put_byte x::pointer y::int = pointer_put_byte x y;
+put_short x::pointer y::int = pointer_put_short x y;
put_int x::pointer y::int = pointer_put_int x y;
+put_float x::pointer y::double = pointer_put_float x y;
put_double x::pointer y::double = pointer_put_double x y;
put_string x::pointer y::string = pointer_put_string x y;
put_pointer x::pointer y::string = pointer_put_pointer x y;
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-23 18:56:29 UTC (rev 841)
+++ pure/trunk/runtime.cc 2008-09-23 21:56:51 UTC (rev 842)
@@ -5301,11 +5301,18 @@
extern "C"
int32_t pointer_get_byte(void *ptr)
{
- uint8_t *p = (uint8_t*)ptr;
+ int8_t *p = (int8_t*)ptr;
return *p;
}
extern "C"
+int32_t pointer_get_short(void *ptr)
+{
+ int16_t *p = (int16_t*)ptr;
+ return *p;
+}
+
+extern "C"
int32_t pointer_get_int(void *ptr)
{
int32_t *p = (int32_t*)ptr;
@@ -5313,6 +5320,13 @@
}
extern "C"
+double pointer_get_float(void *ptr)
+{
+ float *p = (float*)ptr;
+ return *p;
+}
+
+extern "C"
double pointer_get_double(void *ptr)
{
double *p = (double*)ptr;
@@ -5343,11 +5357,18 @@
extern "C"
void pointer_put_byte(void *ptr, int32_t x)
{
- uint8_t *p = (uint8_t*)ptr;
+ int8_t *p = (int8_t*)ptr;
*p = x;
}
extern "C"
+void pointer_put_short(void *ptr, int32_t x)
+{
+ int16_t *p = (int16_t*)ptr;
+ *p = x;
+}
+
+extern "C"
void pointer_put_int(void *ptr, int32_t x)
{
int32_t *p = (int32_t*)ptr;
@@ -5355,6 +5376,13 @@
}
extern "C"
+void pointer_put_float(void *ptr, double x)
+{
+ float *p = (float*)ptr;
+ *p = x;
+}
+
+extern "C"
void pointer_put_double(void *ptr, double x)
{
double *p = (double*)ptr;
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-23 18:56:29 UTC (rev 841)
+++ pure/trunk/runtime.h 2008-09-23 21:56:51 UTC (rev 842)
@@ -767,14 +767,18 @@
you'll have to use the memory management routines above to do that. */
int32_t pointer_get_byte(void *ptr);
+int32_t pointer_get_short(void *ptr);
int32_t pointer_get_int(void *ptr);
+double pointer_get_float(void *ptr);
double pointer_get_double(void *ptr);
char *pointer_get_string(void *ptr);
void *pointer_get_pointer(void *ptr);
pure_expr *pointer_get_expr(void *ptr);
void pointer_put_byte(void *ptr, int32_t x);
+void pointer_put_short(void *ptr, int32_t x);
void pointer_put_int(void *ptr, int32_t x);
+void pointer_put_float(void *ptr, double x);
void pointer_put_double(void *ptr, double x);
void pointer_put_string(void *ptr, const char *x);
void pointer_put_pointer(void *ptr, void *x);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-23 18:59:59
|
Revision: 841
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=841&view=rev
Author: agraef
Date: 2008-09-23 18:56:29 +0000 (Tue, 23 Sep 2008)
Log Message:
-----------
Add missing string conversion routines.
Modified Paths:
--------------
pure/trunk/lib/strings.pure
Modified: pure/trunk/lib/strings.pure
===================================================================
--- pure/trunk/lib/strings.pure 2008-09-23 13:01:30 UTC (rev 840)
+++ pure/trunk/lib/strings.pure 2008-09-23 18:56:29 UTC (rev 841)
@@ -157,6 +157,14 @@
stream s::string = stream (chars s);
tuple s::string = tuple (chars s);
+string [] = "";
+string xs@(_::string:_) = strcat xs if all stringp xs;
+
+/* Conversions between strings and matrices. */
+
+matrix s::string = matrix (chars s);
+string x::matrix = string (list x) if all stringp x;
+
/* Define the customary list operations on strings, so that these can mostly
be used as if they were lists. */
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-23 13:01:40
|
Revision: 840
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=840&view=rev
Author: agraef
Date: 2008-09-23 13:01:30 +0000 (Tue, 23 Sep 2008)
Log Message:
-----------
Comment change.
Modified Paths:
--------------
pure/trunk/lib/matrices.pure
Modified: pure/trunk/lib/matrices.pure
===================================================================
--- pure/trunk/lib/matrices.pure 2008-09-23 12:56:30 UTC (rev 839)
+++ pure/trunk/lib/matrices.pure 2008-09-23 13:01:30 UTC (rev 840)
@@ -352,7 +352,7 @@
garbage-collected. CAVEAT: These are low-level operations and hence should
be used with care. They are useful, in particular, if you need to handle
massive amounts of matrix or vector data generated or processed by external
- C software, such as graphics and sound libraries. */
+ C functions, such as routines dealing with graphics and sound data. */
private matrix_from_double_array;
private matrix_from_complex_array;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-23 12:56:57
|
Revision: 839
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=839&view=rev
Author: agraef
Date: 2008-09-23 12:56:30 +0000 (Tue, 23 Sep 2008)
Log Message:
-----------
Add more matrix creation and pointer->matrix conversion functions.
Modified Paths:
--------------
pure/trunk/lib/matrices.pure
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/lib/matrices.pure
===================================================================
--- pure/trunk/lib/matrices.pure 2008-09-23 10:58:30 UTC (rev 838)
+++ pure/trunk/lib/matrices.pure 2008-09-23 12:56:30 UTC (rev 839)
@@ -43,6 +43,18 @@
rowvectorp x = matrixp x && dim x!0==1;
colvectorp x = matrixp x && dim x!1==1;
+/* Convenience functions to create numeric matrices with the given dimensions
+ (either a pair denoting the number of rows and columns, or just the row
+ size in order to create a row vector), and all elements zero. */
+
+dmatrix (n::int,m::int) = cdmatrix (pointer 0) (n,m);
+cmatrix (n::int,m::int) = ccmatrix (pointer 0) (n,m);
+imatrix (n::int,m::int) = cimatrix (pointer 0) (n,m);
+
+dmatrix n::int = cdmatrix (pointer 0) n;
+cmatrix n::int = ccmatrix (pointer 0) n;
+imatrix n::int = cimatrix (pointer 0) n;
+
/* Size of a matrix (#x) and its dimensions (dim x=n,m where n is the number
of rows, m the number of columns). Note that Pure supports empty matrices,
thus the total size may well be zero, meaning that either the number of
@@ -332,6 +344,37 @@
pointer x::matrix = pure_pointerval x;
+/* Create a numeric matrix from a pointer. The pointer to be converted must
+ point to a malloc'ed, properly initialized memory area big enough to
+ accommodate the requested dimensions. The pointer may also be NULL in which
+ case a suitable pointer is allocated automatically. This memory is taken
+ over by Pure and will be freed automatically when the matrix object is
+ garbage-collected. CAVEAT: These are low-level operations and hence should
+ be used with care. They are useful, in particular, if you need to handle
+ massive amounts of matrix or vector data generated or processed by external
+ C software, such as graphics and sound libraries. */
+
+private matrix_from_double_array;
+private matrix_from_complex_array;
+private matrix_from_int_array;
+extern expr* matrix_from_double_array(void* p, int n, int m);
+extern expr* matrix_from_complex_array(void* p, int n, int m);
+extern expr* matrix_from_int_array(void* p, int n, int m);
+
+cdmatrix p::pointer (n::int,m::int)
+ = matrix_from_double_array p n m;
+ccmatrix p::pointer (n::int,m::int)
+ = matrix_from_complex_array p n m;
+cimatrix p::pointer (n::int,m::int)
+ = matrix_from_int_array p n m;
+
+cdmatrix p::pointer n::int
+ = cdmatrix p (1,n);
+ccmatrix p::pointer n::int
+ = ccmatrix p (1,n);
+cimatrix p::pointer n::int
+ = cimatrix p (1,n);
+
/* Matrix comparisons. */
x::matrix == y::matrix = x === y
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-23 10:58:30 UTC (rev 838)
+++ pure/trunk/runtime.cc 2008-09-23 12:56:30 UTC (rev 839)
@@ -105,7 +105,7 @@
s.tda = m->tda;
s.block = m->block;
s.owner = 0;
- view.matrix = s;
+ view.matrix = s;
return view;
}
}
@@ -5010,6 +5010,95 @@
#endif
}
+extern "C"
+pure_expr *matrix_from_double_array(void *p, uint32_t n1, uint32_t n2)
+{
+#ifdef HAVE_GSL
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_double_matrix(create_double_matrix(n1, n2));
+ if (!p) p = calloc(n1*n2, sizeof(double));
+ if (!p) return 0;
+ gsl_matrix_view v = gsl_matrix_view_array((double*)p, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix *m = (gsl_matrix*)malloc(sizeof(gsl_matrix));
+ gsl_block *b = (gsl_block*)malloc(sizeof(gsl_block));
+ assert(m && v.matrix.data);
+ *m = v.matrix;
+ b->size = n1*n2;
+ b->data = m->data;
+ m->block = b;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::DMATRIX;
+ x->data.mat.p = m;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_from_complex_array(void *p, uint32_t n1, uint32_t n2)
+{
+#ifdef HAVE_GSL
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_complex_matrix(create_complex_matrix(n1, n2));
+ if (!p) p = calloc(2*n1*n2, sizeof(double));
+ if (!p) return 0;
+ gsl_matrix_complex_view v =
+ gsl_matrix_complex_view_array((double*)p, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix_complex *m =
+ (gsl_matrix_complex*)malloc(sizeof(gsl_matrix_complex));
+ gsl_block_complex *b = (gsl_block_complex*)malloc(sizeof(gsl_block_complex));
+ assert(m && v.matrix.data);
+ *m = v.matrix;
+ b->size = n1*n2;
+ b->data = m->data;
+ m->block = b;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::CMATRIX;
+ x->data.mat.p = m;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_from_int_array(void *p, uint32_t n1, uint32_t n2)
+{
+#ifdef HAVE_GSL
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_int_matrix(create_int_matrix(n1, n2));
+ if (!p) p = calloc(n1*n2, sizeof(int));
+ if (!p) return 0;
+ gsl_matrix_int_view v = gsl_matrix_int_view_array((int*)p, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix_int *m = (gsl_matrix_int*)malloc(sizeof(gsl_matrix_int));
+ gsl_block_int *b = (gsl_block_int*)malloc(sizeof(gsl_block_int));
+ assert(m && v.matrix.data);
+ *m = v.matrix;
+ b->size = n1*n2;
+ b->data = m->data;
+ m->block = b;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::IMATRIX;
+ x->data.mat.p = m;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
+#else
+ return 0;
+#endif
+}
+
static uint32_t mpz_hash(const mpz_t z)
{
uint32_t h = 0;
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-23 10:58:30 UTC (rev 838)
+++ pure/trunk/runtime.h 2008-09-23 12:56:30 UTC (rev 839)
@@ -733,6 +733,18 @@
pure_expr *matrix_re(pure_expr *x);
pure_expr *matrix_im(pure_expr *x);
+/* Create a matrix object of the given dimensions which uses the given pointer
+ p as its underlying C array. There are no checks whatsoever, so the caller
+ is responsible for making sure that the memory has the right size and is
+ properly initialized. p must point to a malloc'ed memory area, which is
+ taken over by Pure and will be freed automatically when the matrix is
+ garbage-collected. p may also be NULL in which case a suitable pointer is
+ allocated automatically. */
+
+pure_expr *matrix_from_double_array(void *p, uint32_t n, uint32_t m);
+pure_expr *matrix_from_complex_array(void *p, uint32_t n, uint32_t m);
+pure_expr *matrix_from_int_array(void *p, uint32_t n, uint32_t m);
+
/* Compute a 32 bit hash code of a Pure expression. This makes it possible to
use arbitary Pure values as keys in a hash table. */
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-23 10:58:40
|
Revision: 838
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=838&view=rev
Author: agraef
Date: 2008-09-23 10:58:30 +0000 (Tue, 23 Sep 2008)
Log Message:
-----------
Make slicing work with matrix ranges.
Modified Paths:
--------------
pure/trunk/lib/matrices.pure
Modified: pure/trunk/lib/matrices.pure
===================================================================
--- pure/trunk/lib/matrices.pure 2008-09-23 10:48:11 UTC (rev 837)
+++ pure/trunk/lib/matrices.pure 2008-09-23 10:58:30 UTC (rev 838)
@@ -92,8 +92,14 @@
former is specified as a list of int values, the latter as a pair of lists
of int values. As with list slicing, index ranges may be non-contiguous
and/or non-monotonous. However, the case of contiguous and monotonous
- ranges is optimized by making good use of the 'submat' operation below. */
+ ranges is optimized by making good use of the 'submat' operation below. We
+ also add some rules to make slicing work with ranges drawn from a matrix
+ rather than a list. */
+x!!ns::matrix = x!!list ns;
+x!!(ns::matrix,ms::matrix)
+ = x!!(list ns,list ms);
+
x::matrix!!(ns,ms) = case ns,ms of
ns@(n:_),ms@(m:_) = submat x (n,m) (#ns,#ms)
if cont ns && cont ms;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-23 10:48:19
|
Revision: 837
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=837&view=rev
Author: agraef
Date: 2008-09-23 10:48:11 +0000 (Tue, 23 Sep 2008)
Log Message:
-----------
Add optimization rules for "void" matrix comprehensions.
Modified Paths:
--------------
pure/trunk/lib/matrices.pure
Modified: pure/trunk/lib/matrices.pure
===================================================================
--- pure/trunk/lib/matrices.pure 2008-09-23 10:31:58 UTC (rev 836)
+++ pure/trunk/lib/matrices.pure 2008-09-23 10:48:11 UTC (rev 837)
@@ -190,7 +190,8 @@
extern expr* matrix_rows(expr *x) = rowcat;
extern expr* matrix_columns(expr *x) = colcat;
-/* Combinations of rowcat/colcat and map. */
+/* Combinations of rowcat/colcat and map. These are used, in particular, for
+ implementing matrix comprehensions. */
rowcatmap f [] = {};
rowcatmap f xs@(_:_) = rowcat (map f xs);
@@ -198,6 +199,12 @@
colcatmap f [] = {};
colcatmap f xs@(_:_) = colcat (map f xs);
+/* Optimization rules for "void" matrix comprehensions (cf. the catmap
+ optimization rules at the beginning of prelude.pure). */
+
+def void (rowcatmap f x) = do f x;
+def void (colcatmap f x) = do f x;
+
/* Transpose a matrix. */
private matrix_transpose;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-23 10:32:03
|
Revision: 836
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=836&view=rev
Author: agraef
Date: 2008-09-23 10:31:58 +0000 (Tue, 23 Sep 2008)
Log Message:
-----------
Moved the matrix operations from primitives.pure to their own module.
Modified Paths:
--------------
pure/trunk/ChangeLog
pure/trunk/lib/prelude.pure
pure/trunk/lib/primitives.pure
Added Paths:
-----------
pure/trunk/lib/matrices.pure
Modified: pure/trunk/ChangeLog
===================================================================
--- pure/trunk/ChangeLog 2008-09-23 09:15:34 UTC (rev 835)
+++ pure/trunk/ChangeLog 2008-09-23 10:31:58 UTC (rev 836)
@@ -1,5 +1,8 @@
2008-09-23 Albert Graef <Dr....@t-...>
+ * lib/matrices.pure: Moved the matrix operations from
+ primitives.pure to their own module.
+
* lib/primitives.pure: Added a bunch of new matrix operations. In
particular, list operations like filter and map now work on
matrices, too.
Added: pure/trunk/lib/matrices.pure
===================================================================
--- pure/trunk/lib/matrices.pure (rev 0)
+++ pure/trunk/lib/matrices.pure 2008-09-23 10:31:58 UTC (rev 836)
@@ -0,0 +1,339 @@
+
+/* Basic matrix functions. */
+
+/* Copyright (c) 2008 by Albert Graef <Dr....@t-...>.
+
+ This file is part of the Pure programming language and system.
+
+ Pure is free software: you can redistribute it and/or modify it under the
+ terms of the GNU General Public License as published by the Free Software
+ Foundation, either version 3 of the License, or (at your option) any later
+ version.
+
+ Pure is distributed in the hope that it will be useful, but WITHOUT ANY
+ WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+ details.
+
+ You should have received a copy of the GNU General Public License along
+ with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+/* Additional matrix predicates. The numeric matrix types will only be
+ available if the interpreter was built with support for the GNU Scientific
+ Library (GSL). This is the case iff the global gsl_version variable is
+ predefined by the interpreter. */
+
+extern int matrix_type(expr* x);
+
+dmatrixp x = case x of _::matrix = matrix_type x==1; _ = 0 end;
+cmatrixp x = case x of _::matrix = matrix_type x==2; _ = 0 end;
+imatrixp x = case x of _::matrix = matrix_type x==3; _ = 0 end;
+
+/* The nmatrixp predicate checks for any kind of numeric (double, complex or
+ int) matrix, smatrix for symbolic matrices. */
+
+nmatrixp x = case x of _::matrix = matrix_type x>=1; _ = 0 end;
+smatrixp x = case x of _::matrix = matrix_type x==0; _ = 0 end;
+
+/* Pure represents row and column vectors as matrices with 1 row or column,
+ respectively. The following predicates allow you to check for these special
+ kinds of matrices. */
+
+vectorp x = matrixp x && (n==1 || m==1 when n::int,m::int = dim x end);
+rowvectorp x = matrixp x && dim x!0==1;
+colvectorp x = matrixp x && dim x!1==1;
+
+/* Size of a matrix (#x) and its dimensions (dim x=n,m where n is the number
+ of rows, m the number of columns). Note that Pure supports empty matrices,
+ thus the total size may well be zero, meaning that either the number of
+ rows or the number of columns is zero, or both. The null predicate checks
+ for empty matrices. */
+
+private matrix_size matrix_dim;
+extern int matrix_size(expr *x), expr* matrix_dim(expr *x);
+
+#x::matrix = matrix_size x;
+dim x::matrix = matrix_dim x;
+null x::matrix = #x==0;
+
+/* The stride of a matrix denotes the real row size of the underlying C array,
+ see the description of the 'pack' function below for further details.
+ There's little use for this value in Pure, but it may be needed when
+ interfacing to C. */
+
+private matrix_stride;
+extern int matrix_stride(expr *x);
+
+stride x::matrix = matrix_stride x;
+
+/* Indexing. Note that in difference to Octave and MATLAB, all indices are
+ zero-based, and you *must* use Pure's indexing operator '!' to retrieve an
+ element. You can either get an element by specifying its row-major index in
+ the range 0..#x-1, or with a two-dimensional subscript of the form of a
+ pair (i,j), where i denotes the row and j the column index. Both operations
+ take only constant time, and an 'out_of_bounds' exception is thrown if an
+ index falls outside the valid range. */
+
+private matrix_elem_at matrix_elem_at2;
+extern expr* matrix_elem_at(expr* x, int i);
+extern expr* matrix_elem_at2(expr* x, int i, int j);
+
+x::matrix!i::int = matrix_elem_at x i if i>=0 && i<#x;
+ = throw out_of_bounds otherwise;
+
+x::matrix!(i::int,j::int)
+ = matrix_elem_at2 x i j
+ if (i>=0 && i<n && j>=0 && j<m
+ when n::int,m::int = dim x end);
+ = throw out_of_bounds otherwise;
+
+/* Matrix slices (x!!ns). As with simple indexing, elements can addressed
+ using either singleton (row-major) indices or index pair (row,column). The
+ former is specified as a list of int values, the latter as a pair of lists
+ of int values. As with list slicing, index ranges may be non-contiguous
+ and/or non-monotonous. However, the case of contiguous and monotonous
+ ranges is optimized by making good use of the 'submat' operation below. */
+
+x::matrix!!(ns,ms) = case ns,ms of
+ ns@(n:_),ms@(m:_) = submat x (n,m) (#ns,#ms)
+ if cont ns && cont ms;
+ _ = colcatmap (mth (rowcatmap (nth x) ns)) ms;
+ end with
+ cont [n] = 1;
+ cont (n::int:ns@(m::int:_)) = cont ns if m==n+1;
+ cont _ = 0 otherwise;
+ nth x n = catch (cst {}) (row x n);
+ mth x m = catch (cst {}) (col x m);
+ end;
+x::matrix!!ns = if packed x then rowvector x!!([0],ns)
+ else colcatmap (nth x) ns with
+ nth x n = catch (cst {}) {x!n};
+ end;
+
+/* Extract rows and columns from a matrix. */
+
+private matrix_slice;
+extern expr* matrix_slice(expr* x, int i1, int j1, int i2, int j2);
+
+row x::matrix i::int = if i>=0 && i<n then matrix_slice x i 0 i (m-1)
+ else throw out_of_bounds
+ when n::int,m::int = dim x end;
+
+col x::matrix j::int = if j>=0 && j<m then matrix_slice x 0 j (n-1) j
+ else throw out_of_bounds
+ when n::int,m::int = dim x end;
+
+rows x::matrix = map (row x) (0..n-1) when n::int,_ = dim x end;
+
+cols x::matrix = map (col x) (0..m-1) when _,m::int = dim x end;
+
+/* Convert a matrix to a list and vice versa. list x converts a matrix x to a
+ flat list of its elements, while list2 converts it to a list of lists.
+ Conversely, matrix xs converts a list of lists or matrices specifying the
+ rows of the matrix to the corresponding matrix; otherwise, the result is a
+ row vector. NOTE: The matrix function may throw a 'bad_matrix_value x' in
+ case of dimension mismatch, where x denotes the offending submatrix. */
+
+list x::matrix = [x!i|i=0..#x-1];
+list2 x::matrix = [[x!(i,j)|j=0..m-1]|i=0..n-1]
+ when n::int,m::int = dim x end;
+
+matrix [] = {};
+matrix xs@(x:_) = rowcatmap colcat xs if all listp xs;
+ = rowcat xs if any matrixp xs;
+ = colcat xs otherwise;
+
+/* Extract (sub-,super-) diagonals from a matrix. Sub- and super-diagonals for
+ k=0 return the main diagonal. Indices for sub- and super-diagonals can also
+ be negative, in which case the corresponding super- or sub-diagonal is
+ returned instead. In each case the result is a row vector. */
+
+private matrix_diag matrix_subdiag matrix_supdiag;
+extern expr* matrix_diag(expr* x) = diag;
+extern expr* matrix_subdiag(expr* x, int k) = subdiag;
+extern expr* matrix_supdiag(expr* x, int k) = supdiag;
+
+/* Extract a submatrix of a given size at a given offset. The result shares
+ the underlying storage with the input matrix (i.e., matrix elements are
+ *not* copied) and so this is a comparatively cheap operation. */
+
+submat x::matrix (i::int,j::int) (n::int,m::int)
+ = matrix_slice x i j (i+n-1) (j+m-1);
+
+/* Change the dimensions of a matrix without changing its size. The total
+ number of elements must match that of the input matrix. Reuses the
+ underlying storage of the input matrix if possible. */
+
+private matrix_redim;
+extern expr* matrix_redim(expr* x, int n, int m);
+
+redim x::matrix (n::int,m::int)
+ = matrix_redim x n m if n>=0 && m>=0 && n*m==#x;
+
+/* You can also redim a matrix to a given row size. In this case the row size
+ must be positive and divide the total size of the matrix, */
+
+redim x::matrix m::int = redim x (#x div m,m) if m>0 && #x mod m==0;
+
+/* Convenience functions for converting a matrix to a row or column vector. */
+
+rowvector x::matrix = redim x (#x);
+colvector x::matrix = redim x 1;
+
+/* Construct matrices from lists of rows and columns. These take either
+ scalars or submatrices as inputs; corresponding dimensions must match.
+ rowcat combines submatrices vertically, like {x;y}; colcat combines them
+ horizontally, like {x,y}. NOTE: Like the built-in matrix constructs, these
+ operations may throw a 'bad_matrix_value x' in case of dimension mismatch,
+ where x denotes the offending submatrix. */
+
+extern expr* matrix_rows(expr *x) = rowcat;
+extern expr* matrix_columns(expr *x) = colcat;
+
+/* Combinations of rowcat/colcat and map. */
+
+rowcatmap f [] = {};
+rowcatmap f xs@(_:_) = rowcat (map f xs);
+
+colcatmap f [] = {};
+colcatmap f xs@(_:_) = colcat (map f xs);
+
+/* Transpose a matrix. */
+
+private matrix_transpose;
+extern expr* matrix_transpose(expr *x);
+
+x::matrix' = matrix_transpose x;
+
+/* Reverse a matrix. 'rowrev' reverses the rows, 'colrev' the columns,
+ 'reverse' both dimensions. */
+
+rowrev x::matrix = rowcat (reverse (rows x));
+colrev x::matrix = colcat (reverse (cols x));
+reverse x::matrix = rowrev (colrev x);
+
+/* catmap et al on matrices. This allows list and matrix comprehensions to
+ draw values from matrices as well as from lists, treating the matrix as a
+ flat list of its elements. */
+
+catmap f x::matrix = catmap f (list x);
+rowcatmap f x::matrix = rowcat (map f (list x));
+colcatmap f x::matrix = colcat (map f (list x));
+
+/* Implementations of the other customary list operations, so that these can
+ be used on matrices, too. These operations treat the matrix essentially as
+ if it was a flat list of its elements. Aggregate results are then converted
+ back to matrices with the appropriate dimensions. (This depends on the
+ particular operation; functions like map and zip keep the dimensions of the
+ input matrix intact, while other functions like filter, take or scanl
+ always return a flat row vector. Also note that the zip-style operations
+ require that the row sizes of all arguments match.) */
+
+cycle x::matrix = cycle (list x);
+cyclen n::int x::matrix = cyclen n (list x) if not null x;
+
+all p x::matrix = all p (list x);
+any p x::matrix = any p (list x);
+do f x::matrix = do f (list x);
+drop k::int x::matrix = x!!(k..#x-1);
+dropwhile p x::matrix = colcat (dropwhile p (list x));
+filter p x::matrix = colcat (filter p (list x));
+foldl f a x::matrix = foldl f a (list x);
+foldl1 f x::matrix = foldl1 f (list x);
+foldr f a x::matrix = foldr f a (list x);
+foldr1 f x::matrix = foldr1 f (list x);
+head x::matrix = x!0 if not null x;
+init x::matrix = x!!(0..#x-2) if not null x;
+last x::matrix = x!(#x-1) if not null x;
+map f x::matrix = flip redim (dim x) $ colcat (map f (list x));
+scanl f a x::matrix = colcat (scanl f a (list x));
+scanl1 f x::matrix = colcat (scanl1 f (list x));
+scanr f a x::matrix = colcat (scanr f a (list x));
+scanr1 f x::matrix = colcat (scanr1 f (list x));
+take k::int x::matrix = x!!(0..k-1);
+takewhile p x::matrix = colcat (takewhile p (list x));
+tail x::matrix = x!!(1..#x-1) if not null x;
+
+zip x::matrix y::matrix = flip redim (dim x!1) $
+ colcat (zip (list x) (list y))
+ if dim x!1==dim y!1;
+zip3 x::matrix y::matrix z::matrix
+ = flip redim (dim x!1) $
+ colcat (zip3 (list x) (list y) (list z))
+ if dim x!1==dim y!1 && dim x!1==dim z!1;
+zipwith f x::matrix y::matrix
+ = flip redim (dim x!1) $
+ colcat (zipwith f (list x) (list y))
+ if dim x!1==dim y!1;
+zipwith3 f x::matrix y::matrix z::matrix
+ = flip redim (dim x!1) $
+ colcat (zipwith3 f (list x) (list y) (list z))
+ if dim x!1==dim y!1 && dim x!1==dim z!1;
+dowith f x::matrix y::matrix
+ = dowith f (list x) (list y)
+ if dim x!1==dim y!1;
+dowith3 f x::matrix y::matrix z::matrix
+ = dowith3 f (list x) (list y) (list z)
+ if dim x!1==dim y!1 && dim x!1==dim z!1;
+
+unzip x::matrix = flip redim (dim x) (colcat u),
+ flip redim (dim x) (colcat v)
+ when u,v = unzip (list x) end;
+unzip3 x::matrix = flip redim (dim x) (colcat u),
+ flip redim (dim x) (colcat v),
+ flip redim (dim x) (colcat w)
+ when u,v,w = unzip3 (list x) end;
+
+/* Matrix conversions. These convert between different types of numeric
+ matrices. You can also extract the real and imaginary parts of a (complex)
+ matrix. */
+
+private matrix_double matrix_complex matrix_int;
+extern expr* matrix_double(expr *x), expr* matrix_complex(expr *x),
+ expr* matrix_int(expr *x);
+
+dmatrix x::matrix = matrix_double x if imatrixp x || dmatrixp x;
+imatrix x::matrix = matrix_int x if imatrixp x || dmatrixp x;
+cmatrix x::matrix = matrix_complex x if nmatrixp x;
+
+private matrix_re matrix_im;
+extern expr* matrix_re(expr *x), expr* matrix_im(expr *x);
+
+re x::matrix = matrix_re x if nmatrixp x;
+im x::matrix = matrix_im x if nmatrixp x;
+
+/* 'Pack' a matrix. This creates a copy of the matrix which has the data in
+ contiguous storage. If possible, it also frees up extra memory if the
+ matrix was created as a slice from a bigger matrix. The 'packed' predicate
+ can be used to verify whether a matrix is already packed. */
+
+pack x::matrix = colcat [x,{}];
+packed x::matrix = stride x==dim x!1;
+
+/* Convert a matrix to a pointer. This returns a pointer to the underlying C
+ array, which allows you to change the data in-place by shelling out to C
+ functions, so beware. You also need to be careful when passing such a
+ pointer to C functions if the underlying data is non-contiguous; when in
+ doubt use the 'pack' function from above to place the data in contiguous
+ storage first. */
+
+pointer x::matrix = pure_pointerval x;
+
+/* Matrix comparisons. */
+
+x::matrix == y::matrix = x === y
+ if nmatrixp x && matrix_type x == matrix_type y;
+// mixed numeric cases
+ = cmatrix x === y if nmatrixp x && cmatrixp y;
+ = x === cmatrix y if cmatrixp x && nmatrixp y;
+ = dmatrix x === y if imatrixp x && dmatrixp y;
+ = x === dmatrix y if dmatrixp x && imatrixp y;
+// comparisons with symbolic matrices
+ = 0 if dim x != dim y;
+ = compare 0 with
+ compare i::int = 1 if i>=n;
+ = 0 if x!i != y!i;
+ = compare (i+1);
+ end when n::int = #x end;
+
+x::matrix != y::matrix = not x == y;
Modified: pure/trunk/lib/prelude.pure
===================================================================
--- pure/trunk/lib/prelude.pure 2008-09-23 09:15:34 UTC (rev 835)
+++ pure/trunk/lib/prelude.pure 2008-09-23 10:31:58 UTC (rev 836)
@@ -74,7 +74,7 @@
Note that the math and system modules are *not* included here, so you have
to do that yourself if your program requires any of those operations. */
-using primitives, strings;
+using primitives, matrices, strings;
/* Basic combinators. */
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-23 09:15:34 UTC (rev 835)
+++ pure/trunk/lib/primitives.pure 2008-09-23 10:31:58 UTC (rev 836)
@@ -41,33 +41,8 @@
doublep x = case x of _::double = 1; _ = 0 end;
stringp x = case x of _::string = 1; _ = 0 end;
pointerp x = case x of _::pointer = 1; _ = 0 end;
-
-/* Matrix predicates. The numeric matrix types will only be available if the
- interpreter was built with support for the GNU Scientific Library (GSL).
- This is the case iff the global gsl_version variable is predefined by the
- interpreter. */
-
-extern int matrix_type(expr* x);
-
matrixp x = case x of _::matrix = 1; _ = 0 end;
-dmatrixp x = case x of _::matrix = matrix_type x==1; _ = 0 end;
-cmatrixp x = case x of _::matrix = matrix_type x==2; _ = 0 end;
-imatrixp x = case x of _::matrix = matrix_type x==3; _ = 0 end;
-/* The nmatrixp predicate checks for any kind of numeric (double, complex or
- int) matrix, smatrix for symbolic matrices. */
-
-nmatrixp x = case x of _::matrix = matrix_type x>=1; _ = 0 end;
-smatrixp x = case x of _::matrix = matrix_type x==0; _ = 0 end;
-
-/* Pure represents row and column vectors as matrices with 1 row or column,
- respectively. The following predicates allow you to check for these special
- kinds of matrices. */
-
-vectorp x = matrixp x && (n==1 || m==1 when n::int,m::int = dim x end);
-rowvectorp x = matrixp x && dim x!0==1;
-colvectorp x = matrixp x && dim x!1==1;
-
/* Predicates to check for function objects, global (unbound) variables,
function applications, proper lists, list nodes and tuples. */
@@ -117,8 +92,7 @@
pointer x::int |
pointer x::bigint |
pointer x::double |
-pointer x::string |
-pointer x::matrix = pure_pointerval x;
+pointer x::string = pure_pointerval x;
/* Convert signed (8/16/32/64) bit integers to the corresponding unsigned
quantities. These functions behave as if the value was "cast" to the
@@ -410,266 +384,6 @@
x::pointer==y::pointer = bigint x == bigint y;
x::pointer!=y::pointer = bigint x != bigint y;
-/* Basic matrix operations: size, dimensions and indexing. */
-
-private matrix_size matrix_dim matrix_stride;
-extern int matrix_size(expr *x), int matrix_stride(expr *x);
-extern expr* matrix_dim(expr *x);
-
-null x::matrix = #x==0;
-#x::matrix = matrix_size x;
-dim x::matrix = matrix_dim x;
-stride x::matrix = matrix_stride x;
-
-private matrix_elem_at matrix_elem_at2;
-extern expr* matrix_elem_at(expr* x, int i);
-extern expr* matrix_elem_at2(expr* x, int i, int j);
-
-x::matrix!i::int = matrix_elem_at x i if i>=0 && i<#x;
- = throw out_of_bounds otherwise;
-
-x::matrix!(i::int,j::int)
- = matrix_elem_at2 x i j
- if (i>=0 && i<n && j>=0 && j<m
- when n::int,m::int = dim x end);
- = throw out_of_bounds otherwise;
-
-/* Matrix slices. */
-
-x::matrix!!(ns,ms) = case ns,ms of
- // optimize the case of contiguous slices
- ns@(n:_),ms@(m:_) = submat x (n,m) (#ns,#ms)
- if cont ns && cont ms;
- _ = colcatmap (mth (rowcatmap (nth x) ns)) ms;
- end with
- cont [n] = 1;
- cont (n::int:ns@(m::int:_)) = cont ns if m==n+1;
- cont _ = 0 otherwise;
- nth x n = catch (cst {}) (row x n);
- mth x m = catch (cst {}) (col x m);
- end;
-x::matrix!!ns = if packed x then rowvector x!!([0],ns)
- else colcatmap (nth x) ns with
- nth x n = catch (cst {}) {x!n};
- end;
-
-/* Extract rows and columns from a matrix. */
-
-private matrix_slice;
-extern expr* matrix_slice(expr* x, int i1, int j1, int i2, int j2);
-
-row x::matrix i::int = if i>=0 && i<n then matrix_slice x i 0 i (m-1)
- else throw out_of_bounds
- when n::int,m::int = dim x end;
-
-col x::matrix j::int = if j>=0 && j<m then matrix_slice x 0 j (n-1) j
- else throw out_of_bounds
- when n::int,m::int = dim x end;
-
-rows x::matrix = map (row x) (0..n-1) when n::int,_ = dim x end;
-
-cols x::matrix = map (col x) (0..m-1) when _,m::int = dim x end;
-
-/* Convert a matrix to a list and vice versa. list x converts a matrix x to a
- flat list of its elements, while list2 converts it to a list of lists.
- Conversely, matrix xs converts a list of lists or matrices specifying the
- rows of the matrix to the corresponding matrix; otherwise, the result is a
- row vector. NOTE: The matrix function may throw a 'bad_matrix_value x' in
- case of dimension mismatch, where x denotes the offending submatrix. */
-
-list x::matrix = [x!i|i=0..#x-1];
-list2 x::matrix = [[x!(i,j)|j=0..m-1]|i=0..n-1]
- when n::int,m::int = dim x end;
-
-matrix [] = {};
-matrix xs@(x:_) = rowcatmap colcat xs if all listp xs;
- = rowcat xs if any matrixp xs;
- = colcat xs otherwise;
-
-/* Extract (sub-,super-) diagonals from a matrix. Sub- and super-diagonals for
- k=0 return the main diagonal. Indices for sub- and super-diagonals can also
- be negative, in which case the corresponding super- or sub-diagonal is
- returned instead. In each case the result is a row vector. */
-
-private matrix_diag matrix_subdiag matrix_supdiag;
-extern expr* matrix_diag(expr* x) = diag;
-extern expr* matrix_subdiag(expr* x, int k) = subdiag;
-extern expr* matrix_supdiag(expr* x, int k) = supdiag;
-
-/* Extract a submatrix of a given size at a given offset. The result shares
- the underlying storage with the input matrix (i.e., matrix elements are
- *not* copied) and so this is a comparatively cheap operation. */
-
-submat x::matrix (i::int,j::int) (n::int,m::int)
- = matrix_slice x i j (i+n-1) (j+m-1);
-
-/* Change the dimensions of a matrix without changing its size. The total
- number of elements must match that of the input matrix. Reuses the
- underlying storage of the input matrix if possible. */
-
-private matrix_redim;
-extern expr* matrix_redim(expr* x, int n, int m);
-
-redim x::matrix (n::int,m::int)
- = matrix_redim x n m if n>=0 && m>=0 && n*m==#x;
-
-/* You can also redim a matrix to a given positive row size. In this case the
- row size must divide the total size of the matrix, */
-
-redim x::matrix m::int = redim x (#x div m,m) if m>0 && #x mod m==0;
-
-/* Convenience functions for converting a matrix to a row or column vector. */
-
-rowvector x::matrix = redim x (#x);
-colvector x::matrix = redim x 1;
-
-/* Construct matrices from lists of rows and columns. These take either
- scalars or submatrices as inputs; corresponding dimensions must match.
- rowcat combines submatrices vertically, like {x;y}; colcat combines them
- horizontally, like {x,y}. NOTE: Like the built-in matrix constructs, these
- operations may throw a 'bad_matrix_value x' in case of dimension mismatch,
- where x denotes the offending submatrix. */
-
-extern expr* matrix_rows(expr *x) = rowcat;
-extern expr* matrix_columns(expr *x) = colcat;
-
-/* Combinations of rowcat/colcat and map. */
-
-rowcatmap f [] = {};
-rowcatmap f xs@(_:_) = rowcat (map f xs);
-
-colcatmap f [] = {};
-colcatmap f xs@(_:_) = colcat (map f xs);
-
-/* 'Pack' a matrix. This creates a copy of the matrix which has the data in
- contiguous storage. If possible, it also frees up extra memory if the
- matrix was created as a slice from a bigger matrix. The 'packed' predicate
- can be used to verify whether a matrix is already packed. */
-
-pack x::matrix = colcat [x,{}];
-packed x::matrix = stride x==dim x!1;
-
-/* Transpose a matrix. */
-
-private matrix_transpose;
-extern expr* matrix_transpose(expr *x);
-
-x::matrix' = matrix_transpose x;
-
-/* Reverse a matrix. rowrev reverses the rows, colrev the columns, reverse
- both dimensions. */
-
-rowrev x::matrix = rowcat (reverse (rows x));
-colrev x::matrix = colcat (reverse (cols x));
-reverse x::matrix = rowrev (colrev x);
-
-/* catmap et al on matrices. This allows list and matrix comprehensions to
- draw values from matrices just as well as from lists, treating the matrix
- as a flat list of its elements. */
-
-catmap f x::matrix = catmap f (list x);
-rowcatmap f x::matrix = rowcat (map f (list x));
-colcatmap f x::matrix = colcat (map f (list x));
-
-/* Implementations of the other customary list operations, so that these can
- be used on matrices, too. These operations treat the matrix essentially as
- if it was a flat list of its elements. Aggregate results are then converted
- back to matrices with the appropriate dimensions. (This depends on the
- operation; operations like map and zip keep the dimensions of the input
- matrix intact, while other functions like filter, take or scanl always
- return a flat row vector. Also note that the zip-style operations require
- that the row sizes of all arguments match.) */
-
-cycle x::matrix = cycle (list x);
-cyclen n::int x::matrix = cyclen n (list x) if not null x;
-
-all p x::matrix = all p (list x);
-any p x::matrix = any p (list x);
-do f x::matrix = do f (list x);
-drop k::int x::matrix = x!!(k..#x-1);
-dropwhile p x::matrix = colcat (dropwhile p (list x));
-filter p x::matrix = colcat (filter p (list x));
-foldl f a x::matrix = foldl f a (list x);
-foldl1 f x::matrix = foldl1 f (list x);
-foldr f a x::matrix = foldr f a (list x);
-foldr1 f x::matrix = foldr1 f (list x);
-head x::matrix = x!0 if not null x;
-init x::matrix = x!!(0..#x-2) if not null x;
-last x::matrix = x!(#x-1) if not null x;
-map f x::matrix = flip redim (dim x) $ colcat (map f (list x));
-scanl f a x::matrix = colcat (scanl f a (list x));
-scanl1 f x::matrix = colcat (scanl1 f (list x));
-scanr f a x::matrix = colcat (scanr f a (list x));
-scanr1 f x::matrix = colcat (scanr1 f (list x));
-take k::int x::matrix = x!!(0..k-1);
-takewhile p x::matrix = colcat (takewhile p (list x));
-tail x::matrix = x!!(1..#x-1) if not null x;
-
-zip x::matrix y::matrix = flip redim (dim x!1) $
- colcat (zip (list x) (list y))
- if dim x!1==dim y!1;
-zip3 x::matrix y::matrix z::matrix
- = flip redim (dim x!1) $
- colcat (zip3 (list x) (list y) (list z))
- if dim x!1==dim y!1 && dim x!1==dim z!1;
-zipwith f x::matrix y::matrix
- = flip redim (dim x!1) $
- colcat (zipwith f (list x) (list y))
- if dim x!1==dim y!1;
-zipwith3 f x::matrix y::matrix z::matrix
- = flip redim (dim x!1) $
- colcat (zipwith3 f (list x) (list y) (list z))
- if dim x!1==dim y!1 && dim x!1==dim z!1;
-dowith f x::matrix y::matrix
- = dowith f (list x) (list y)
- if dim x!1==dim y!1;
-dowith3 f x::matrix y::matrix z::matrix
- = dowith3 f (list x) (list y) (list z)
- if dim x!1==dim y!1 && dim x!1==dim z!1;
-
-unzip x::matrix = flip redim (dim x) (colcat u),
- flip redim (dim x) (colcat v)
- when u,v = unzip (list x) end;
-unzip3 x::matrix = flip redim (dim x) (colcat u),
- flip redim (dim x) (colcat v),
- flip redim (dim x) (colcat w)
- when u,v,w = unzip3 (list x) end;
-
-/* Matrix conversions. */
-
-private matrix_double matrix_complex matrix_int;
-extern expr* matrix_double(expr *x), expr* matrix_complex(expr *x),
- expr* matrix_int(expr *x);
-
-dmatrix x::matrix = matrix_double x if imatrixp x || dmatrixp x;
-imatrix x::matrix = matrix_int x if imatrixp x || dmatrixp x;
-cmatrix x::matrix = matrix_complex x if nmatrixp x;
-
-private matrix_re matrix_im;
-extern expr* matrix_re(expr *x), expr* matrix_im(expr *x);
-
-re x::matrix = matrix_re x if nmatrixp x;
-im x::matrix = matrix_im x if nmatrixp x;
-
-/* Matrix comparisons. */
-
-x::matrix == y::matrix = x === y
- if nmatrixp x && matrix_type x == matrix_type y;
-// mixed numeric cases
- = cmatrix x === y if nmatrixp x && cmatrixp y;
- = x === cmatrix y if cmatrixp x && nmatrixp y;
- = dmatrix x === y if imatrixp x && dmatrixp y;
- = x === dmatrix y if dmatrixp x && imatrixp y;
-// comparisons with symbolic matrices
- = 0 if dim x != dim y;
- = compare 0 with
- compare i::int = 1 if i>=n;
- = 0 if x!i != y!i;
- = compare (i+1);
- end when n::int = #x end;
-
-x::matrix != y::matrix = not x == y;
-
/* IEEE floating point infinities and NaNs. Place these after the definitions
of the built-in operators so that the double arithmetic works. */
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
