single-python.c

     
   1  //! @file single-python.c
   2  //! @author J. Marcel van der Veer
   3  //!
   4  //! @section Copyright
   5  //!
   6  //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
   7  //! Copyright 2001-2023 J. Marcel van der Veer [algol68g@xs4all.nl].
   8  //!
   9  //! @section License
  10  //!
  11  //! This program is free software; you can redistribute it and/or modify it 
  12  //! under the terms of the GNU General Public License as published by the 
  13  //! Free Software Foundation; either version 3 of the License, or 
  14  //! (at your option) any later version.
  15  //!
  16  //! This program is distributed in the hope that it will be useful, but 
  17  //! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
  18  //! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for 
  19  //! more details. You should have received a copy of the GNU General Public 
  20  //! License along with this program. If not, see [http://www.gnu.org/licenses/].
  21  
  22  //! @section Synopsis
  23  //!
  24  //! REAL vector and matrix routines, Python look-a-likes.
  25  
  26  #include "a68g.h"
  27  
  28  #if defined (HAVE_GSL)
  29  
  30  #include "a68g-genie.h"
  31  #include "a68g-prelude.h"
  32  #include "a68g-prelude-gsl.h"
  33  #include "a68g-torrix.h"
  34  
  35  //! Print REAL vector and matrix using GSL.
  36  
  37  #define FMT "%12.4g"
  38  #define COLS 6
  39  
  40  void print_vector (gsl_vector *A, unt w)
  41  {
  42    unt N = SIZE (A);
  43    fprintf (stdout, "[%d]\n", N);
  44    if (w <= 0 || N <= 2 * w) {
  45      for (unt i = 0, k = 0; i < N; i++, k++) {
  46        if (k > 0 && k % COLS == 0) {
  47          fprintf (stdout, "\n");
  48          k = 0;
  49        }
  50        fprintf (stdout, FMT, gsl_vector_get (A, i));
  51      }
  52    } else {
  53      for (unt i = 0, k = 0; i < w; i++, k++) {
  54        if (k > 0 && k % COLS == 0) {
  55          fprintf (stdout, "\n");
  56          k = 0;
  57        }
  58        fprintf (stdout, FMT, gsl_vector_get (A, i));
  59      }
  60      fprintf (stdout, " ... ");
  61      for (unt i = N - w, k = 0; i < N; i++, k++) {
  62        if (k > 0 && k % COLS == 0) {
  63          fprintf (stdout, "\n");
  64          k = 0;
  65        }
  66        fprintf (stdout, FMT, gsl_vector_get (A, i));
  67      }
  68    }
  69    fprintf (stdout, "\n");
  70  }
  71  
  72  void print_row (gsl_matrix *A, unt m, unt w)
  73  {
  74    unt N = SIZE2 (A);
  75    if (w <= 0 || N <= 2 * w) {
  76      for (unt i = 0, k = 0; i < N; i++, k++) {
  77        if (k > 0 && k % COLS == 0) {
  78          fprintf (stdout, "\n");
  79          k = 0;
  80        }
  81        fprintf (stdout, FMT, gsl_matrix_get (A, m, i));
  82      }
  83    } else {
  84      for (unt i = 0, k = 0; i < w; i++, k++) {
  85        if (k > 0 && k % COLS == 0) {
  86          fprintf (stdout, "\n");
  87          k = 0;
  88        }
  89        fprintf (stdout, FMT, gsl_matrix_get (A, m, i));
  90      }
  91      fprintf (stdout, " ... ");
  92      for (unt i = N - w, k = 0; i < N; i++, k++) {
  93        if (k > 0 && k % COLS == 0) {
  94          fprintf (stdout, "\n");
  95          k = 0;
  96        }
  97        fprintf (stdout, FMT, gsl_matrix_get (A, m, i));
  98      }
  99    }
 100    fprintf (stdout, "\n");
 101  }
 102  
 103  void print_matrix (gsl_matrix *A, unt w)
 104  {
 105    unt M = SIZE1 (A), N = SIZE2 (A);
 106    fprintf (stdout, "[%d, %d]\n", M, N);
 107    if (w <= 0 || M <= 2 * w) {
 108      for (unt i = 0; i < M; i++) {
 109        print_row (A, i, w);
 110      }
 111    } else {
 112      for (unt i = 0; i < w; i++) {
 113        print_row (A, i, w);
 114      }
 115      fprintf (stdout, " ...\n");
 116      for (unt i = M - w; i < M; i++) {
 117        print_row (A, i, w);
 118      }
 119    }
 120    fprintf (stdout, "\n");
 121  }
 122  
 123  //! @brief PROC print vector = ([] REAL v, INT width) VOID
 124  
 125  void genie_print_vector (NODE_T *p)
 126  {
 127    A68_INT width;
 128    POP_OBJECT (p, &width, A68_INT);
 129    gsl_vector *V = pop_vector (p, A68_TRUE);
 130    print_vector (V, VALUE (&width));
 131    gsl_vector_free (V);
 132  }
 133  
 134  //! @brief PROC print matrix = ([, ] REAL v, INT width) VOID
 135  
 136  void genie_print_matrix (NODE_T *p)
 137  {
 138    A68_INT width;
 139    POP_OBJECT (p, &width, A68_INT);
 140    gsl_matrix *M = pop_matrix (p, A68_TRUE);
 141    print_matrix (M, VALUE (&width));
 142    gsl_matrix_free (M);
 143  }
 144  
 145  //! @brief Convert VECTOR to [] REAL.
 146  
 147  A68_ROW vector_to_row (NODE_T * p, gsl_vector * v)
 148  {
 149    unt len = (unt) (SIZE (v));
 150    A68_ROW desc, row; A68_ARRAY arr; A68_TUPLE tup;
 151    NEW_ROW_1D (desc, row, arr, tup, M_ROW_REAL, M_REAL, len);
 152    BYTE_T *base = DEREF (BYTE_T, &ARRAY (&arr));
 153    int idx = VECTOR_OFFSET (&arr, &tup);
 154    unt inc = SPAN (&tup) * ELEM_SIZE (&arr);
 155    for (unt k = 0; k < len; k++, idx += inc) {
 156      A68_REAL *x = (A68_REAL *) (base + idx);
 157      STATUS (x) = INIT_MASK;
 158      VALUE (x) = gsl_vector_get (v, (size_t) k);
 159      CHECK_REAL (p, VALUE (x));
 160    }
 161    return desc;
 162  }
 163  
 164  //! @brief Convert MATRIX to [, ] REAL.
 165  
 166  A68_ROW matrix_to_row (NODE_T * p, gsl_matrix * a)
 167  {
 168    int len1 = (int) (SIZE1 (a)), len2 = (int) (SIZE2 (a));
 169    A68_REF desc; A68_ROW row; A68_ARRAY arr; A68_TUPLE tup1, tup2;
 170    desc = heap_generator (p, M_ROW_ROW_REAL, DESCRIPTOR_SIZE (2));
 171    row = heap_generator (p, M_ROW_ROW_REAL, len1 * len2 * SIZE (M_REAL));
 172    DIM (&arr) = 2;
 173    MOID (&arr) = M_REAL;
 174    ELEM_SIZE (&arr) = SIZE (M_REAL);
 175    SLICE_OFFSET (&arr) = FIELD_OFFSET (&arr) = 0;
 176    ARRAY (&arr) = row;
 177    LWB (&tup1) = 1; UPB (&tup1) = len1; SPAN (&tup1) = 1;
 178    SHIFT (&tup1) = LWB (&tup1); K (&tup1) = 0;
 179    LWB (&tup2) = 1; UPB (&tup2) = len2; SPAN (&tup2) = ROW_SIZE (&tup1);
 180    SHIFT (&tup2) = LWB (&tup2) * SPAN (&tup2); K (&tup2) = 0;
 181    PUT_DESCRIPTOR2 (arr, tup1, tup2, &desc);
 182    BYTE_T *base = DEREF (BYTE_T, &ARRAY (&arr));
 183    int idx1 = MATRIX_OFFSET (&arr, &tup1, &tup2);
 184    int inc1 = SPAN (&tup1) * ELEM_SIZE (&arr), inc2 = SPAN (&tup2) * ELEM_SIZE (&arr);
 185    for (int k1 = 0; k1 < len1; k1++, idx1 += inc1) {
 186      for (int k2 = 0, idx2 = idx1; k2 < len2; k2++, idx2 += inc2) {
 187        A68_REAL *x = (A68_REAL *) (base + idx2);
 188        STATUS (x) = INIT_MASK;
 189        VALUE (x) = gsl_matrix_get (a, (size_t) k1, (size_t) k2);
 190        CHECK_REAL (p, VALUE (x));
 191      }
 192    }
 193    return desc;
 194  }
 195  
 196  //! @brief OP NORM = ([, ] REAL) REAL
 197  
 198  REAL_T matrix_norm (gsl_matrix *a)
 199  {
 200  #if (A68_LEVEL >= 3)
 201    DOUBLE_T sum = 0.0;
 202    unt M = SIZE1 (a), N = SIZE2 (a);
 203    for (unt i = 0; i < M; i++) {
 204      for (unt j = 0; j < N; j++) {
 205        DOUBLE_T z = gsl_matrix_get (a, i, j);
 206        sum += z * z;
 207      }
 208    }
 209    return ((REAL_T) sqrtq (sum));
 210  #else
 211    REAL_T sum = 0.0;
 212    unt M = SIZE1 (a), N = SIZE2 (a);
 213    for (unt i = 0; i < M; i++) {
 214      for (unt j = 0; j < N; j++) {
 215        REAL_T z = gsl_matrix_get (a, i, j);
 216        sum += z * z;
 217      }
 218    }
 219    return (sqrt (sum));
 220  #endif
 221  }
 222  
 223  void genie_matrix_norm (NODE_T * p)
 224  {
 225    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
 226    torrix_error_node = p;
 227    gsl_matrix *a = pop_matrix (p, A68_TRUE);
 228    PUSH_VALUE (p, matrix_norm (a), A68_REAL);
 229    gsl_matrix_free (a);
 230    (void) gsl_set_error_handler (save_handler);
 231  }
 232  
 233  //! @brief OP HCAT = ([, ] REAL, [, ] REAL) [, ] REAL
 234  
 235  gsl_matrix *matrix_hcat (NODE_T *p, gsl_matrix *u, gsl_matrix *v)
 236  {
 237    unt Mv = SIZE1 (v), Nv = SIZE2 (v);
 238    unt Mu, Nu;
 239    if (u == NO_REAL_MATRIX) {
 240      Mu = Nu = 0;
 241    } else {
 242      Mu = SIZE1 (u), Nu = SIZE2 (u);
 243    }
 244    if (Mu == 0 && Nu == 0) {
 245      Mu = Mv;
 246    } else {
 247      MATH_RTE (p, Mu != Mv, M_INT, "number of rows does not match");
 248    }
 249    gsl_matrix *w = gsl_matrix_calloc (Mu, Nu + Nv);
 250    for (unt i = 0; i < Mu; i++) {
 251      unt k = 0;
 252      for (unt j = 0; j < Nu; j++) {
 253        gsl_matrix_set (w, i, k++, gsl_matrix_get (u, i, j));
 254      }
 255      for (unt j = 0; j < Nv; j++) {
 256        gsl_matrix_set (w, i, k++, gsl_matrix_get (v, i, j));
 257      }
 258    }
 259    return (w);
 260  }
 261  
 262  void genie_matrix_hcat (NODE_T *p)
 263  {
 264  // Yield [u v].
 265    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
 266    torrix_error_node = p;
 267    gsl_matrix *v = pop_matrix (p, A68_TRUE);
 268    gsl_matrix *u = pop_matrix (p, A68_TRUE);
 269    gsl_matrix *w = matrix_hcat (p, u, v);
 270    push_matrix (p, w);
 271    gsl_matrix_free (u);
 272    gsl_matrix_free (v);
 273    gsl_matrix_free (w);
 274    (void) gsl_set_error_handler (save_handler);
 275  }
 276  
 277  //! @brief OP VCAT = ([, ] REAL, [, ] REAL) [, ] REAL
 278  
 279  gsl_matrix *matrix_vcat (NODE_T *p, gsl_matrix *u, gsl_matrix *v)
 280  {
 281    unt Mv = SIZE1 (v), Nv = SIZE2 (v);
 282    unt Mu, Nu;
 283    if (u == NO_REAL_MATRIX) {
 284      Mu = Nu = 0;
 285    } else {
 286      Mu = SIZE1 (u), Nu = SIZE2 (u);
 287    }
 288    if (Mu == 0 && Nu == 0) {
 289      Nu = Nv;
 290    } else {  
 291      MATH_RTE (p, Nu != Nv, M_INT, "number of columns does not match");
 292    }
 293    gsl_matrix *w = gsl_matrix_calloc (Mu + Mv, Nu);
 294    for (unt j = 0; j < Nu; j++) {
 295      unt k = 0;
 296      for (unt i = 0; i < Mu; i++) {
 297        gsl_matrix_set (w, k++, j, gsl_matrix_get (u, i, j));
 298      }
 299      for (unt i = 0; i < Mv; i++) {
 300        gsl_matrix_set (w, k++, j, gsl_matrix_get (v, i, j));
 301      }
 302    }
 303    return (w);
 304  }
 305  
 306  void genie_matrix_vcat (NODE_T *p)
 307  {
 308  // Yield [u; v].
 309    gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
 310    torrix_error_node = p;
 311    gsl_matrix *v = pop_matrix (p, A68_TRUE);
 312    gsl_matrix *u = pop_matrix (p, A68_TRUE);
 313    gsl_matrix *w = matrix_vcat (p, u, v);
 314    push_matrix (p, w);
 315    gsl_matrix_free (u);
 316    gsl_matrix_free (v);
 317    gsl_matrix_free (w);
 318    (void) gsl_set_error_handler (save_handler);
 319  }
 320  
 321  #endif