share/amatrix/amatrix.mac

   1 /* copyright 2007 by Robert Dodier
   2  * I release this file under the terms of the GNU General Public License.
   3  */
   4
   5 /* operations on amatrix objects */
   6
   7 matchdeclare ([A%, A1%, A2%], amatrixp);
   8 matchdeclare ([X%, X1%, X2%], lambda ([e], not amatrixp (e)));
   9
  10 /* troublesome ??: tellsimp (is (A%), amatrixmap (is, A%)); */
  11
  12 tellsimpafter (A% < X%, amatrixmap (buildq ([X%], lambda ([a%], a% < X%)), A%));
  13 tellsimpafter (A% <= X%, amatrixmap (buildq ([X%], lambda ([a%], a% <= X%)), A%));
  14 tellsimpafter (A% = X%, amatrixmap (buildq ([X%], lambda ([a%], a% = X%)), A%));
  15 tellsimpafter (equal (A%, X%), amatrixmap (buildq ([X%], lambda ([a%], equal (a%, X%))), A%));
  16 tellsimpafter (notequal (A%, X%), amatrixmap (buildq ([X%], lambda ([a%], notequal (a%, X%))), A%));
  17 tellsimpafter (A% # X%, amatrixmap (buildq ([X%], lambda ([a%], a% # X%)), A%));
  18 tellsimpafter (A% >= X%, amatrixmap (buildq ([X%], lambda ([a%], a% >= X%)), A%));
  19 tellsimpafter (A% > X%, amatrixmap (buildq ([X%], lambda ([a%], a% > X%)), A%));
  20
  21 tellsimpafter (X% < A%, amatrixmap (buildq ([X%], lambda ([a%], X% < a%)), A%));
  22 tellsimpafter (X% <= A%, amatrixmap (buildq ([X%], lambda ([a%], X% <= a%)), A%));
  23 tellsimpafter (X% = A%, amatrixmap (buildq ([X%], lambda ([a%], X% = a%)), A%));
  24 tellsimpafter (equal (X%, A%), amatrixmap (buildq ([X%], lambda ([a%], equal (X%, a%))), A%));
  25 tellsimpafter (notequal (X%, A%), amatrixmap (buildq ([X%], lambda ([a%], notequal (X%, a%))), A%));
  26 tellsimpafter (X% # A%, amatrixmap (buildq ([X%], lambda ([a%], X% # a%)), A%));
  27 tellsimpafter (X% >= A%, amatrixmap (buildq ([X%], lambda ([a%], X% >= a%)), A%));
  28 tellsimpafter (X% > A%, amatrixmap (buildq ([X%], lambda ([a%], X% > a%)), A%));
  29
  30 tellsimpafter (A1% < A2%, amatrixmap ("<", A1%, A2%));
  31 tellsimpafter (A1% <= A2%, amatrixmap ("<=", A1%, A2%));
  32 tellsimpafter (A1% = A2%, amatrixmap ("=", A1%, A2%));
  33 tellsimpafter (equal (A1%, A2%), amatrixmap ('equal, A1%, A2%));
  34 tellsimpafter (notequal (A1%, A2%), amatrixmap ('notequal, A1%, A2%));
  35 tellsimpafter (A1% # A2%, amatrixmap ("#", A1%, A2%));
  36 tellsimpafter (A1% >= A2%, amatrixmap (">=", A1%, A2%));
  37 tellsimpafter (A1% > A2%, amatrixmap (">", A1%, A2%));
  38
  39 tellsimpafter (not A%, amatrixmap ("not", A%));
  40 tellsimpafter (A1% and A2%, amatrixmap (lambda ([a1%, a2%], a1% and a2%), A1%, A2%));
  41 tellsimpafter (A1% or A2%, amatrixmap (lambda ([a1%, a2%], a1% or a2%), A1%, A2%));
  42
  43 /* implementation */
  44
  45 defstruct (amatrix (nr, r0, rinc, nc, c0, cinc, storage));
  46
  47 elements (m) := create_list (get_element (m, i, j), i, 1, m@nr, j, 1, m@nc);
  48
  49 make_matrix  (nr, nc, [a]) := block
  50   ([aa : ?gensym()],
  51     ?putprop (aa, 1, '?refcount),
  52
  53     if a = []
  54     then a : make_array (fixnum, nr*nc)
  55     /* ELSE STORAGE ARRAY IS SPECIFIED; SHOULD VERIFY HERE THAT
  56      * (1) IT IS 1-DIMENSIONAL, AND (2) IT IS BIG ENOUGH
  57      */
  58     else a : a[1],
  59
  60     ?putprop (aa, a, '?storage_array),
  61     new (amatrix (nr, 0, 1, nc, 0, nr, aa)));
  62
  63 /* DECLARING A% SYMBOLP CAUSES TROUBLE IF STORAGE GENSYM IS EVALUATED ...
  64  * CONSIDER ATTACHING ARRAY TO GENSYM AS A PROPERTY INSTEAD OF A VALUE (SIGH)
  65  */
  66 matchdeclare
  67    (a%, symbolp,
  68     [LI%, LI1%, LI2%], integer_listp,
  69     [MB%, MB1%, MB2%], putative_boolean_amatrixp,
  70     [ii%, jj%, i0%, j0%, u%, v%, x%, y%], integerp);
  71
  72 integer_listp (e) := not atom(e) and op(e) = "[" and every (integerp, e);
  73 putative_boolean_amatrixp (e) :=
  74   not atom(e)
  75     and op(e) = amatrix
  76     and every (lambda ([x], not atom(x) or booleanp(x)), elements (e));
  77 booleanp (e) := is (e = true or e = false);
  78
  79 /* USEFUL FOR DEBUGGING
  80 announce_rules_firing : true;
  81  */
  82
  83 incprop (aa, bb) := ?putprop (aa, ?get (aa, bb) + 1, bb);
  84
  85 decprop (aa, bb) := ?putprop (aa, ?get (aa, bb) - 1, bb);
  86
  87 simp : false;
  88
  89 get_element (m%, ii%, jj%) :=
  90     get_element_internal (m%@nr, m%@r0, m%@rinc, m%@nc, m%@c0, m%@cinc, m%@storage, ii%, jj%);
  91
  92 get_element_internal (u%, i0%, x%, v%, j0%, y%, a%, ii%, jj%) := block
  93
  94   ([a% : ?get (a%, '?storage_array),
  95     i% : (i0% + ii% - 1)*x% + (j0% + jj% - 1)*y%],
  96
  97     if ii% < 1 or jj% < 1 or ii% > u% or jj% > v%
  98     then error ("Matrix index out of bounds")
  99     else ?aref (a%, i%));
 100
 101 /* one integer */
 102
 103 tellsimpafter
 104     (amatrix (1, i0%, x%, v%, j0%, y%, a%) [ii%],
 105      get_element_internal (1, i0%, x%, v%, j0%, y%, a%, 1, ii%));
 106
 107 tellsimpafter
 108     (amatrix (u%, i0%, x%, 1, j0%, y%, a%) [ii%],
 109      get_element_internal (u%, i0%, x%, 1, j0%, y%, a%, ii%, 1));
 110
 111 /* integers both */
 112
 113 tellsimpafter
 114     (amatrix (u%, i0%, x%, v%, j0%, y%, a%) [ii%, jj%],
 115      get_element_internal (u%, i0%, x%, v%, j0%, y%, a%, ii%, jj%));
 116
 117 /* one 'all */
 118
 119 tellsimpafter
 120     (amatrix (1, i0%, x%, v%, j0%, y%, a%) [all],
 121     (incprop (a%, '?refcount),
 122      amatrix (1, i0%, x%, v%, j0%, y%, a%)));
 123
 124 tellsimpafter
 125     (amatrix (u%, i0%, x%, 1, j0%, y%, a%) [all],
 126     (incprop (a%, '?refcount),
 127      amatrix (u%, i0%, x%, 1, j0%, y%, a%)));
 128
 129 /* 'all and an integer */
 130
 131 tellsimpafter
 132     (amatrix (u%, i0%, x%, v%, j0%, y%, a%) [all, jj%],
 133     (incprop (a%, '?refcount),
 134      amatrix (u%, i0%, x%, 1, j0% + jj% - 1, y%, a%)));
 135
 136 tellsimpafter
 137     (amatrix (u%, i0%, x%, v%, j0%, y%, a%) [ii%, all],
 138     (incprop (a%, '?refcount),
 139      amatrix (1, i0% + ii% - 1, x%, v%, j0%, y%, a%)));
 140
 141 /* 'all both */
 142
 143 tellsimpafter
 144     (amatrix (u%, i0%, x%, v%, j0%, y%, a%) [all, all],
 145     (incprop (a%, '?refcount),
 146      amatrix (u%, i0%, x%, v%, j0%, y%, a%)));
 147
 148 /* one list of integers */
 149
 150 tellsimpafter
 151     (amatrix (1, i0%, x%, v%, j0%, y%, a%) [LI%],
 152      submatrix_by_indices (1, i0%, x%, v%, j0%, y%, a%, [1], LI%));
 153
 154 tellsimpafter
 155     (amatrix (u%, i0%, x%, 1, j0%, y%, a%) [LI%],
 156      submatrix_by_indices (u%, i0%, x%, 1, j0%, y%, a%, LI%, [1]));
 157
 158 /* list of integers and an integer */
 159
 160 tellsimpafter
 161     (amatrix (u%, i0%, x%, v%, j0%, y%, a%) [ii%, LI%],
 162      submatrix_by_indices (u%, i0%, x%, v%, j0%, y%, a%, [ii%], LI%));
 163
 164 tellsimpafter
 165     (amatrix (u%, i0%, x%, v%, j0%, y%, a%) [LI%, jj%],
 166      submatrix_by_indices (u%, i0%, x%, v%, j0%, y%, a%, LI%, [jj%]));
 167
 168 /* list of integers and 'all */
 169
 170 tellsimpafter
 171     (amatrix (u%, i0%, x%, v%, j0%, y%, a%) [all, LI%],
 172      submatrix_by_indices (u%, i0%, x%, v%, j0%, y%, a%, makelist (i, i, 1, u%), LI%));
 173
 174 tellsimpafter
 175     (amatrix (u%, i0%, x%, v%, j0%, y%, a%) [LI%, all],
 176      submatrix_by_indices (u%, i0%, x%, v%, j0%, y%, a%, LI%, makelist (i, i, 1, v%)));
 177
 178 /* lists of integers both */
 179
 180 tellsimpafter
 181     (amatrix (u%, i0%, x%, v%, j0%, y%, a%) [LI1%, LI2%],
 182      submatrix_by_indices (u%, i0%, x%, v%, j0%, y%, a%, LI1%, LI2%));
 183
 184 /* one amatrix of booleans */
 185
 186 tellsimpafter
 187     (amatrix (1, i0%, x%, v%, j0%, y%, a%) [MB%],
 188      submatrix_by_indices_and_flags (1, i0%, x%, v%, j0%, y%, a%, [1], MB%));
 189
 190 tellsimpafter
 191     (amatrix (u%, i0%, x%, 1, j0%, y%, a%) [MB%],
 192      submatrix_by_flags_and_indices (u%, i0%, x%, 1, j0%, y%, a%, MB%, [1]));
 193
 194 /* amatrix of booleans and an integer */
 195
 196 tellsimpafter
 197     (amatrix (u%, i0%, x%, v%, j0%, y%, a%) [ii%, MB%],
 198      submatrix_by_indices_and_flags (u%, i0%, x%, v%, j0%, y%, a%, [ii%], MB%));
 199
 200 tellsimpafter
 201     (amatrix (u%, i0%, x%, v%, j0%, y%, a%) [MB%, jj%],
 202      submatrix_by_flags_and_indices (u%, i0%, x%, v%, j0%, y%, a%, MB%, [jj%]));
 203
 204 /* amatrix of booleans and 'all */
 205
 206 tellsimpafter
 207     (amatrix (u%, i0%, x%, v%, j0%, y%, a%) [all, MB%],
 208      submatrix_by_all_and_flags (u%, i0%, x%, v%, j0%, y%, a%, MB%));
 209
 210 tellsimpafter
 211     (amatrix (u%, i0%, x%, v%, j0%, y%, a%) [MB%, all],
 212      submatrix_by_flags_and_all (u%, i0%, x%, v%, j0%, y%, a%, MB%));
 213
 214 /* amatrix of booleans and a list of integers */
 215
 216 tellsimpafter
 217     (amatrix (u%, i0%, x%, v%, j0%, y%, a%) [LI%, MB%],
 218      submatrix_by_indices_and_flags (u%, i0%, x%, v%, j0%, y%, a%, LI%, MB%));
 219
 220 tellsimpafter
 221     (amatrix (u%, i0%, x%, v%, j0%, y%, a%) [MB%, LI%],
 222      submatrix_by_flags_and_indices (u%, i0%, x%, v%, j0%, y%, a%, MB%, LI%));
 223
 224 /* amatrix of booleans both */
 225
 226 tellsimpafter
 227     (amatrix (u%, i0%, x%, v%, j0%, y%, a%) [MB1%, MB2%],
 228      submatrix_by_flags_and_flags (u%, i0%, x%, v%, j0%, y%, a%, MB1%, MB2%));
 229
 230 simp : true;
 231
 232 submatrix_by_indices (u%, i0%, x%, v%, j0%, y%, a%, LI1%, LI2%) := block
 233    ([M0 : amatrix (u%, i0%, x%, v%, j0%, y%, a%),
 234      M : make_matrix (length (LI1%), length (LI2%)),
 235      ?\*afterflag : false],
 236     /* Do this the simple way.
 237      * Probably some special cases could be faster.
 238      */
 239     for i:1 thru nrows(M)
 240         do for j:1 thru ncols(M)
 241             do M[i, j] : get_element (M0, LI1%[i], LI2%[j]),
 242     M);
 243
 244 submatrix_by_indices_and_flags (u%, i0%, x%, v%, j0%, y%, a%, LI%, MB%) := block
 245    ([LI2% : sublist_indices (elements (MB%), lambda ([x], is (is(x) = true)))],
 246     submatrix_by_indices (u%, i0%, x%, v%, j0%, y%, a%, LI%, LI2%));
 247
 248 submatrix_by_flags_and_indices (u%, i0%, x%, v%, j0%, y%, a%, MB%, LI%) := block
 249    ([LI1% : sublist_indices (elements (MB%), lambda ([x], is (is(x) = true)))],
 250     submatrix_by_indices (u%, i0%, x%, v%, j0%, y%, a%, LI1%, LI%));
 251
 252 submatrix_by_all_and_flags (u%, i0%, x%, v%, j0%, y%, a%, MB%) := block
 253   ([LI1% : makelist (i, i, 1, u%),
 254     LI2% : sublist_indices (elements (MB%), lambda ([x], is (is(x) = true)))],
 255     submatrix_by_indices (u%, i0%, x%, v%, j0%, y%, a%, LI1%, LI2%));
 256
 257 submatrix_by_flags_and_all (u%, i0%, x%, v%, j0%, y%, a%, MB%) := block
 258   ([LI1% : sublist_indices (elements (MB%), lambda ([x], is (is(x) = true)))],
 259     LI2% : makelist (i, i, 1, v%),
 260     submatrix_by_indices (u%, i0%, x%, v%, j0%, y%, a%, LI1%, LI2%));
 261
 262 submatrix_by_flags_and_flags (u%, i0%, x%, v%, j0%, y%, a%, MB1%, MB2%) := block
 263   ([LI1% : sublist_indices (elements (MB1%), lambda ([x], is (is(x) = true))),
 264     LI2% : sublist_indices (elements (MB2%), lambda ([x], is (is(x) = true)))],
 265     submatrix_by_indices (u%, i0%, x%, v%, j0%, y%, a%, LI1%, LI2%));
 266
 267 t(m):=
 268     (incprop (m@storage, '?refcount),
 269      amatrix
 270         (m@nc, m@c0,
 271          if m@rinc = 1 then m@nr else 1,
 272          m@nr, m@r0,
 273          if m@cinc = 1 then m@nc else 1,
 274          m@storage));
 275
 276 compute_index1 (aa, i, j) :=
 277     if i < 1 or j < 1 or i > aa@nr or j > aa@nc
 278     then error ("Matrix index out of bounds")
 279     else 1 + (aa@r0 + i - 1)*aa@rinc + (aa@c0 + j - 1)*aa@cinc;
 280
 281 compute_index0 (aa, i, j) :=
 282     if i < 0 or j < 0 or i >= aa@nr or j >= aa@nc
 283     then error ("Matrix index out of bounds")
 284     else (aa@r0 + i)*aa@rinc + (aa@c0 + j)*aa@cinc;
 285
 286 /* REVISIT THE FOLLOWING !! IT IS BROKEN AS IT STANDS !! */
 287 copy_amatrix (m) :=
 288     if atom(m) or op(m) # 'amatrix
 289     then copy (m)
 290     else
 291        (incprop (m@storage, '?refcount),
 292         amatrix (m@nr, m@r0, m@rinc, m@nc, m@c0, m@cinc, m@storage));
 293         /* copy(m) copies everything but the storage array -- that's just what we want */
 294         /* comment out until storage is actually an array
 295         block ([m2 : copy(m)], m2@refcount : m2@refcount + 1, m2);
 296          */
 297
 298 copy_array (a) := if ?arrayp (a) then ?copy\-seq (a) else ?copy\-tree (a);
 299
 300 if not ?get ('amatrix, 'present) then load ("amatrix.lisp");
 301
 302 /* SHOULD TRY TO MERGE THIS WITH EXISTING MEAN FUNCTION IN SHARE PACKAGE DESCRIPTIVE */
 303 amatrix_mean (x) := block
 304   ([m : nrows(x),
 305     n : ncols(x)],
 306     if n=1
 307         then 1/m * sum (x [i], i, 1, m)
 308         else makelist (mean (x [all, j]), j, 1, n));
 309
 310 matchdeclare (A%, amatrixp);
 311 tellsimp (mean (A%), amatrix_mean (A%));
 312
 313 /* FOLLOWING NEEDED MOSTLY FOR TESTING ALTHOUGH IT WOULD BE OK TO LEAVE IT IN */
 314
 315 random_matrix (nr, nc) := block
 316   ([M : make_matrix (nr, nc), A],
 317     A : ?get (M@storage, '?storage_array),
 318     for i from 0 thru nr*nc-1
 319         do A[i] : floor (random (256.0)) / 256.0,
 320     M);
 321
 322 new2old_matrix (M) :=
 323   if nrows(M) = 0 or ncols(M) = 0
 324     then matrix()
 325     else genmatrix (lambda ([i, j], get_element (M, i, j)), nrows(M), ncols(M));
 326
 327
 328 old2new_matrix (M) := block ([nr, nc, M2],
 329     nr : length (M),
 330     nc : if nr = 0 then 0 else length (M [1]),
 331     M2 : make_matrix (nr, nc),
 332     for i thru nr
 333         do for j thru nc
 334             do M2 [i, j] : M [i, j],
 335     M2);
 336
 337