search

NOAA-PMEL / Ferret

The Ferret program from NOAA/PMEL
https://ferret.pmel.noaa.gov/Ferret/
The Unlicense
55 stars 20 forks source link

Functions for Matrix computations #1780

Closed karlmsmith closed 6 years ago

karlmsmith commented 6 years ago

Reported by @AnsleyManke on 16 Feb 2017 00:55 UTC As one of our GFDL deliverables for Spring 2017 we had this item:

''Improved support for matrix algebra. Ferret currently uses only array-based (element-wise) data handling. Functionality could be added to work with 2-D array data as a matrix, such as transpose and matrix multiplication.''

There may be more to offer at a later time, maybe via calls to Python libraries. However as a first step, add some Ferret functions, e.g.

let AB = matrix_product(a,b)

We already have TRANSPOSE functions in all permutations of XYZT. If we add more functions to include the E and F dimensions then all of these would just need to be listed with other MATRIX functions in the documentation.

yes? sh func *trans*
TRANSPOSE_XT(VAR)
    transposes X and T axes of given variable
    VAR: variable to transpose in X and T
TRANSPOSE_XY(VAR)
    transposes x and y axes of given variable
    VAR: variable to transpose in X and Y
TRANSPOSE_XZ(VAR)
    transposes x and z axes of given variable
    VAR: variable to transpose in X and Z
TRANSPOSE_YT(VAR)
    transposes Y and T axes of given variable
    VAR: variable to transpose in Y and T
TRANSPOSE_YZ(VAR)
    transposes Y and Z axes of given variable
    VAR: variable transposed in Y and Z
TRANSPOSE_ZT(VAR)
    transposes Z and T axes of given variable
    VAR: variable to transpose in Z and T

Migrated-From: http://dunkel.pmel.noaa.gov/trac/ferret/ticket/2508

karlmsmith commented 6 years ago

Comment by @karlmsmith on 16 Feb 2017 17:55 UTC Looking at the transpose functions make me think the matrix_product may need either similar variants, or an '/ALONG=' to indicate the axis that they have in common (if it is not obvious from the grids) - or at least the same number of points along the axis.

karlmsmith commented 6 years ago

Comment by steven.c.hankin on 16 Feb 2017 18:04 UTC I second your comment.

It would make sense to decide upon a syntax (function names and args) and motivate it with an example script or two in which the matrix math is used in realistic calculations.

karlmsmith commented 6 years ago

Comment by @AnsleyManke on 16 Feb 2017 18:39 UTC Two things could be done for implementation:

The notes from last year's meeting didn't give us too much idea what the use cases or motivation might be. -thoughts, anyone? I think it came out of a discussion about extensions via Python library calls, and definitely things in the way of matrix-solving functions should happen that way rather than via Ferret-side functions.

karlmsmith commented 6 years ago

Comment by @karlmsmith on 16 Feb 2017 20:05 UTC This could be implemented as a PyEF (PyFerret external function written in Python).

karlmsmith commented 6 years ago

Comment by @AnsleyManke on 21 Mar 2017 18:28 UTC Comments about Matrix Multiply from email exchanges,

I was mulling over how the mechanics of matrix multiplies would look in the context of Ferret's 6-dimensional grids, and realized quickly that you were right. There is only a single "inner" axis (or the inner product), and all of the other (non-inner) axes act equivalently. Which means that a general matrix multiply function can be as simple as

MATMUL(a, b, ax)

where

    a and b are typical Ferret 6d gridded variables.

        a and b must share the inner axis, "ax".  On all other dimensions they must either share the same axis, or be normal to that axis.
        The inner product axis, ax, will disappear from the result.

    ax is the axis of the inner product

        With no change at all to Ferret syntax ax can be expressed as, say "Y[gy=a]".  Alternatively an integer 1, ...6.
        (As a somewhat separate matter, it might be nice to develop a syntax that can pass an axis name to a Ferret function without any actual evaluation.  Say new pre-defined variables that are in fact constants (1, ...,6) with names like x_ax, y_ax, or whatever.)

Simplest example,

If Axy is a variable in the XY plane, and Byz is in the YZ plane then

    MATMUL(Axy, Byz, Y[gy=axy])

is the matrix product -- lying in the XZ plane

Higher dimension example,

    If Axyt is a variable on an XYT grid, and Byz is in the YZ plane then

        MATMUL(Axyt, Byz, Y[gy=axy])

    lies in the XZT space.  It can be thought of as the same matrix product as the previous -- lying in the XZ plane -- but evaluated over a range of T indices.
karlmsmith commented 6 years ago

Comment by @AnsleyManke on 21 Mar 2017 18:39 UTC for transpose operations see new ticket,#2519

karlmsmith commented 6 years ago

Comment by @AnsleyManke on 23 Mar 2017 19:32 UTC The description in comment 5 is logically clear, but can't be implemented using the external functions framework. The function definition needs to know exactly how the grid of the result depends on the arguments. It can't be set to keep all the axes from the arguments except for one which will be determined when the user calls the function. This is the same reason we have the proliferation of functions in all directions where if the user wants to SORT or COMPRESS or whatever in X or Z, those are different function calls.

So, I suppose it's six functions: inner_product_x, inner_product_y, etc, and then maybe a script

yes? go define_inner_product a b "y" result_var

which would set up the Y[gy=a] and call the right function to define the result.

karlmsmith commented 6 years ago

Comment by steven.c.hankin on 24 Mar 2017 14:25 UTC Yea, I did play fast and loose with the syntax in comment 5. But I think there is a way that inner products can be a single function. As long as the argument that indicates which axis is to be used as the inner product is a constant, that is explicitly given in the function call, then the result grid has no "data dependencies". I think there is already code to support this in the GCF_ routines. If so, this should work:

MATMUL(Axyt, Byz, 2)

The last argument would be 1, 2, 3, ..., 6 for the 6 axes.

To make this more elegant a nice addition would be the introduction of some pre-defined constants into the Ferret syntax. With that addition, the previous function call might look like (say)

MATMUL(Axyt, Byz, _Yaxis)

In addition to predefined constants _Xaxis, ....,_Faxis (as 1,...,6), _PI would be a nice addition, too. For their implementation, predefined constants would probably want to be a new "cat_" in Ferret, so they could be readily detected as constants in the GCF routines. Predefined constants would not be a lot of work, and would have some nice payoffs.

karlmsmith commented 6 years ago

Modified by @AnsleyManke on 24 Mar 2017 22:58 UTC

karlmsmith commented 6 years ago

Comment by @AnsleyManke on 12 Apr 2017 22:26 UTC I've introduced a new function type, "direction-changing functions" which are a subset of grid-changing functions. They have designated arguments which are direction arguments. The grid inheritance is initially set to inherit axes from some function arguments, but then the result grid may be changed further when the function is evaluated, so the grid setup needs to always define a new result grid.

The INNERPRODUCT(VAR1,VAR2,IDIM) function is implemented this way, to return the inner product computed along the "idim" direction.

Also as a test, a TRANSPOSE(VAR,DIM1,DIM2) function, which calls one of the existing transpose_xy, transpose_yz, etc. functions depending on the values of arguments 2 and 3.

I have not implemented pre-defined constants as suggested in comment 8. The functions are called with integers pointing to the axis direction, 1 means X, 2 means Y, etc. I will open a new ticket on that and on further Matrix operations such as an SVD.

karlmsmith commented 6 years ago

Comment by @AndrewWittenberg on 24 Apr 2017 22:56 UTC Thanks -- this is great! Suggestions:

1) To make scripts more readable, shorten the function name to just "DOT()".

2) If ARG3 (axis direction) isn't present and there's only 1 shared axis, then just use that shared axis. If there's no shared axis or more than one, then raise an error.

3) If ARG3 has more than 1 axis, then do a grand inner product of ARG1 & ARG2 along all of ARG3's axes at once. All axes of ARG3 must be present in both ARG1 & ARG2 (otherwise raise an error). So dot(XYZT,XYZT,XY) would have shape ZT.

With those changes, here are what some use cases would look like:

1) Scalar dot product:

yes? let a = {1,2,3,4}
yes? let b = {-1,6,1,0}

yes? say `dot(a,b)`
14

2) Matrix multiplication:

yes? let a = reshape({1,2,3, -1,0,2}, x[gx=1:3:1]+y[gy=1:2:1])
yes? let b = reshape({4,-2,0,5, 1,3,2,-2, 3,1,0,-1}, x[gx=1:3:1]+z[gz=1:4:1])
yes? list a
             VARIABLE : RESHAPE({1,2,3, -1,0,2}, X[GX=1:3:1]+Y[GY=1:2:1])
             SUBSET   : 3 by 2 points (X-Y)
             1      2      3
             1      2      3
 1   / 1:  1.000  2.000  3.000
 2   / 2: -1.000  0.000  2.000
yes? list b
             VARIABLE : RESHAPE({4,-2,0,5, 1,3,2,-2, 3,1,0,-1}, X[GX=1:3:1]+Z[GZ=1:4:1])
             SUBSET   : 3 by 4 points (X-Z)
             1      2      3
             1      2      3
 1   / 1:  4.000 -2.000  0.000
 2   / 2:  5.000  1.000  3.000
 3   / 3:  2.000 -2.000  3.000
 4   / 4:  1.000  0.000 -1.000
yes? list dot(a,b)
             VARIABLE : DOT(A,B)
             SUBSET   : 2 by 4 points (Y-Z)
             1      2
             1      2
 1   / 1:   0.00  -4.00
 2   / 2:  16.00   1.00
 3   / 3:   7.00   4.00
 4   / 4:  -2.00  -3.00

3) Signal synthesis:

yes? let twopi = 2*3.14159265
yes? let tseries = sin(twopi*t[gt=0:5:.05]/esequence({1,2,4}))
yes? let patterns = sin(twopi*x[gx=0:4:.05]/esequence({3,.5,1}))
yes? shade dot(tseries, patterns)
yes? shade dot(tseries[m=1], patterns[m=1])
yes? shade dot(tseries[m=1:2], patterns[m=1:2])

4) Harmonic analysis:

yes? let twopi = 2*3.14159265
yes? let tval = t[gt=0:`16-.1`:.1]
yes? let tseries = sin(twopi*tval) + 2*cos(twopi*tval/2) + 3*sin(twopi*tval/4)
yes? let period = fsequence({1,2,4,8})
yes? let elements = fcat(sin(twopi*tval/period), cos(twopi*tval/period))
yes? let amplitudes = 2*dot(tseries,elements)/tseries[t=@ngd]
yes? list amplitudes
             VARIABLE : 2*INNERPRODUCT(TSERIES,ELEMENTS,4)/TSERIES[T=@NGD]
             SUBSET   : 8 points (F)
 1   / 1:  1.000
 2   / 2:  0.000
 3   / 3:  3.000
 4   / 4:  0.000
 5   / 5: -0.000
 6   / 6:  2.000
 7   / 7: -0.000
 8   / 8: -0.000

Related to these would be projecting signals onto EOFs or singular vectors.