Skip to content

Latest commit

 

History

History
181 lines (138 loc) · 5.69 KB

scipy.org

File metadata and controls

181 lines (138 loc) · 5.69 KB
Raw

Scipy : high-level scientific computing

This section was adapted from the Scipy lecture notes, which was written by Gaël Varoquaux, Adrien Chauve, Andre Espaze, Emmanuelle Gouillart, and Ralf Gommers.

License: Creative Commons Attribution 4.0 International License (CC-by)

Adapted for py4cl by B.Dudson (2019).

The following examples assume that you have already loaded py4cl

(ql:quickload :py4cl)

File input/output: scipy.io

Matlab files: Loading and saving

Import the io module:

(py4cl:import-module "numpy" :as "np")
(py4cl:import-module "scipy.io" :as "spio")
(py4cl:import-function "dict") ; Use to make dictionaries (hash tables)

;; Create a 3x3 matrix and save to file.mat as "a"
(spio:savemat "file.mat" (dict :a (np:ones '(3 3))))

;; Load file.mat returning a hash table
(let ((data (spio:loadmat "file.mat")))
  (gethash "a" data))

; => #2A((1.0 1.0 1.0) (1.0 1.0 1.0) (1.0 1.0 1.0))

Linear algebra operations: scipy.linalg

The scipy.linalg module provides standard linear algebra operations, relying on an underlying efficient implementation (BLAS, LAPACK).

(py4cl:import-module "scipy.linalg" :as "linalg")

The scipy.linalg.det function computes the determinant of a square matrix:

(linalg:det #2A((1 2) (3 4))) ; => -2.0
(linalg:det (np:ones '(3 4))) ; => PY4CL:PYTHON-ERROR "expected square matrix"

The scipy.linalg.inv function computes the inverse of a square matrix:

(linalg:inv #2A((1 2)
                (3 4)))
; => #2A((-2.0 1.0)
;        (1.5 -0.5))
(let* ((arr #2A((1 2)
                (3 4)))
       (iarr (linalg:inv arr)))
  (np:allclose (np:dot arr iarr) (np:eye 2))) ; => T

Finally computing the inverse of a singular matrix (its determinant is zero) will raise a condition:

(linalg:inv #2A((3 2) (6 4))) ; => PY4CL:PYTHON-ERROR "singular matrix"

More advanced operations are available, for example singular-value decomposition (SVD):

(let ((arr (py4cl:python-eval
              (np:reshape (np:arange 9) '(3 3))
              "+"
              (np:diag '(1 0 1)))))
  (destructuring-bind (uarr spec vharr) (linalg:svd arr)
    spec))  ; => #(14.889826 0.45294237 0.29654968)

The original matrix can be re-composed by matrix multiplication of the outputs of svd with np.dot:

(let ((arr (py4cl:python-eval
              (np:reshape (np:arange 9) '(3 3))
              "+"
              (np:diag '(1 0 1)))))
  (destructuring-bind (uarr spec vharr) (linalg:svd arr)
    (let* ((sarr (np:diag spec))
           (svd_mat (py4cl:chain uarr (dot sarr) (dot vharr))))
      (np:allclose svd_mat arr)))) ; => T

SVD is commonly used in statistics and signal processing. Many other standard decompositions (QR, LU, Cholesky, Schur), as well as solvers for linear systems, are available in scipy.linalg.

Interpolation: scipy.interpolate

scipy.interpolate is useful for fitting a function from experimental data and thus evaluating points where no measure exists. The module is based on the FITPACK Fortran subroutines.

By imagining experimental data close to a sine function, scipy.interpolate.interp1d can build a linear interpolation function:

(py4cl:import-function "interp1d" :from "scipy.interpolate")

(py4cl:import-module "numpy" :as "np")
(py4cl:import-module "numpy.random" :as "nprandom")
(py4cl:import-module "operator" :as "op") ; For arithmetic
(py4cl:remote-objects*  ; Keep data in python except final result
  (let* ((measured_time (np:linspace 0 1 10))
         (noise (op:mul
                  (op:sub
                    (op:mul (nprandom:random 10) 2)
                    1.0)
                  1e-1))
         (measures (op:add 
                     (np:sin (op:mul (* 2 pi) measured_time))
                     noise)))

    ;; scipy.interpolate.interp1d can build a linear interpolation function:
    (let ((linear_interp (interp1d measured_time measures))
          (interpolation_time (np:linspace 0 1 50)))
      ;;Then the result can be evaluated at the time of interest:
      (py4cl:python-call linear_interp interpolation_time))))

Note that arithmetic with NumPy arrays is currently best done in a remote-objects environment. A 1D array returned to lisp is converted to a list when sent to python. This list then doesn’t support arithmetic, even though the original 1D array did support arithmetic.