View Raw SPL
/*****************************************************************************
* *
* COLPAIRWISE.SPL Copyright (C) 2025 DSP Development Corporation *
* All Rights Reserved *
* *
* Author: Randy Race *
* *
* Synopsis: Processes all pairs of columns with a specified function *
* *
* Revisions: 7 Nov 2025 RRR Creation *
* *
*****************************************************************************/
#if @HELP_COLPAIRWISE
COLPAIRWISE
Purpose: Applies a specified function on all combinations of column pairs
from one or two multi-column input series.
Syntax: COLPAIRWISE("funname", x, y, opt1, opt2, ..., opt2)
"funname" - A string, the name of the processing function.
x - A multi-column series.
y - Optional. A multi-column series, Defaults to X.
optn - Optional. Anything, zero or more optional
arguments supported by "funname".
Returns: A multi-column series. If "F" is the function name and X is the
input series, the columns of the multi-column result are:
C11, C12, ..., C1N, C21, C22, ..., C2N, ..., CN1, CN2, ..., CNN
where Cab is the result of F(col(X, a), col(X, b))
For two input series X and Y, the resulting columns are:
C11, C12, ..., C1N, C21, C22, ..., C2N, ..., CN1, CN2, ..., CNN
where Cab is the result of F(col(X, a), col(Y, b))
Example:
W1: ravel(1..12, 4)
W2: colpairwise("plus", w1)
W1 == {{1, 5, 9},
{2, 6, 10},
{3, 7, 11},
{4, 8, 12}}
W2 == {{2, 6, 10, 6, 10, 14, 10, 14, 18},
{4, 8, 12, 8, 12, 16, 12, 16, 20},
{6, 10, 14, 10, 14, 18, 14, 18, 22},
{8, 12, 16, 12, 16, 20, 16, 20, 24}}
The result contains a total of 9 columns. The first 3 columns
are equivalent to:
col(w1,1)+col(w1,1), col(w1,1)+col(w1,2), col(w1,2)+col(w1,3)
The next three columns are equivalent to:
col(w1,2)+col(w1,1), col(w1,2)+col(w1,2), col(w1,2)+col(w1,3)
And the final three columns are equivalent to:
col(w1,3)+col(w1,1), col(w1,3)+col(w1,2), col(w1,3)+col(w1,3)
Example:
W1: ravel(1..12, 4)
W2: w1
W3: colpairwise("plus", w1, w2)
W1 == {{1, 5, 9},
{2, 6, 10},
{3, 7, 11},
{4, 8, 12}}
W2 == {{1, 5, 9},
{2, 6, 10},
{3, 7, 11},
{4, 8, 12}}
W3 == {{2, 6, 10, 6, 10, 14, 10, 14, 18},
{4, 8, 12, 8, 12, 16, 12, 16, 20},
{6, 10, 14, 10, 14, 18, 14, 18, 22},
{8, 12, 16, 12, 16, 20, 16, 20, 24}}
Same as above, except two identical input series are provided.
Example:
W1: ravel(1..12, 4)
W2: ravel(1..12, 3);
W3: colpairwise("plus", w1, w2)
W4: colpairwise("plus", w2, w1)
W1 == {{1, 5, 9},
{2, 6, 10},
{3, 7, 11},
{4, 8, 12}}
W2 == {{1, 4, 7, 10},
{2, 5, 8, 11},
{3, 6, 9, 12}}
W3 == {{2, 5, 8, 11, 6, 9, 12, 15, 10, 13, 16, 19},
4, 7, 10, 13, 8, 11, 14, 17, 12, 15, 18, 21},
6, 9, 12, 15, 10, 13, 16, 19, 14, 17, 20, 23},
4, 4, 4, 4, 8, 8, 8, 8, 12, 12, 12, 12}}
W4 == {{2, 6, 10, 5, 9, 13, 8, 12, 16, 11, 15, 19},
{4, 8, 12, 7, 11, 15, 10, 14, 18, 13, 17, 21},
{6, 10, 14, 9, 13, 17, 12, 16, 20, 15, 19, 23},
{4, 8, 12, 4, 8, 12, 4, 8, 12, 4, 8, 12}}
W1 is a 4x3 array.
W2 is a 3x4 array.
W3 is a 4x12 array.
W4 is a 4x12 array.
W3 adds each column of W1 to each column of W2. Result is a 4×12
matrix.
W4 performs the same operation, but with inputs reversed. The
result is also 4×12, but with a different column order.
The number of output rows is the maximum of the number of rows
of the two input series. During processing, the values of short
columns are assumed to be zero.
The number of output columns is the product of the number of
columns of the two input series.
Example:
W1: ravel(1..12, 4)
W2: ravel(1..12, 3);
W3: colpairwise("times", w1, w2)
W4: colpairwise("times", w2, w1)
W1 == {{1, 5, 9},
{2, 6, 10},
{3, 7, 11},
{4, 8, 12}}
W2 == {{1, 4, 7, 10},
{2, 5, 8, 11},
{3, 6, 9, 12}}
W3 == {{1, 4, 7, 10, 5, 20, 35, 50, 9, 36, 63, 90},
4, 10, 16, 22, 12, 30, 48, 66, 20, 50, 80, 110},
9, 18, 27, 36, 21, 42, 63, 84, 33, 66, 99, 132},
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}}
W4 == {{1, 5, 9, 4, 20, 36, 7, 35, 63, 10, 50, 90},
{4, 12, 20, 10, 30, 50, 16, 48, 80, 22, 66, 110},
{9, 21, 33, 18, 42, 66, 27, 63, 99, 36, 84, 132},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}}
Similar to above, except the pairwise product is performed. Because
short rows are extended with zeros, the last rows of W3 and W4
are all zeros.
Example:
W1: gsweep(10240, 1/10240, 10, 1000);ravel(w0, 5120);stripchart
W2: w1*w1
W3: cpsd(w1, w2, 1024, 512, 2048, "mag");stripchart;loglog
W4: colpairwise("cpsd", w1, w2, 1024, 512, 2048, "mag");stripchart;loglog
W1 contains two columns of a swept sinewave.
W2 squares the series in W1.
W3 computes the cross-power spectral density of the columns of
W1 and W2. Optional CPSD arguments are specified. The result is the
CPSD of W1 column 1 and W2 column 1, and the CPSD of W1 column
2 and W2 column 2 for a total of 2 columns.
W4 computes the CPSD of all column pairs using the same optional
arguments. The result is 4 columns:
CPSD(col(w1,1),col(w2,1))
CPSD(col(w1,1),col(w2,2))
CPSD(col(w1,2),col(w2,1))
CPSD(col(w1,2),col(w2,2))
Note that column 1 of W4 matches column 1 of W3 and column 4 of
W4 matches column 2 of W3.
Remarks:
COLPAIRWISE applies a specified processing function to every
combination of column pairs from one or two multi-column series.
Useful for generating pairwise transformations, comparisons, or
interactions across columns.
The output column order is row-major: for each column in X, all
pairings with columns in Y are evaluated.
The function "funname" must accept two series inputs and return
a single series output.
For a single input series X, the number of output columns is
numcols(X) * numcols(X). The number of output rows is numrows(X).
For two input series X and Y, the number of output columns is
numcols(X) * numcols(Y). The number of output rows is
max(numrows(X), numrows(Y)).
Optional arguments are passed through unchanged to "funname".
See Also:
Iterate
Loop
#endif
/* process each pair of columns with funname */
colpairwise(funname = "", argv)
{
local k, j, nc1, nc2, nc3, cols, cmd;
local s1 = {}, s2 = {}, matout = {};
/* optional args */
k = 0;
/* function name */
if (not(isstring(funname)) || (strlen(funname) == 0))
{
error(sprintf("%s - Function Name Required", __FUNC__));
}
/* get input args */
loop (j = 1..argvc)
{
if (isarray(getargv(j)))
{
if (isempty(s1))
{
s1 = getargv(j);
}
else if (isempty(s2))
{
s2 = getargv(j);
}
else
{
/* optional series arg */
k++;
setlocal(sprintf("arg_%d", k), getargv(j));
}
}
else
{
/* optional arg */
k++;
setlocal(sprintf("arg_%d", k), getargv(j));
}
}
if (isempty(s1))
{
return({});
}
if (isempty(s2))
{
/* default to process pairs in s1 */
s2 = refseries(s1);
}
/* number of output columns */
nc1 = numcols(s2);
nc2 = numcols(s1);
nc3 = nc1 * nc2;
/* build function call */
cmd = "feval(funname, s2, col(s1, j)";
/* add any optional arguments */
loop (j = 1..k)
{
cmd += sprintf(", arg_%d", j);
}
/* closing ( */
cmd += ")";
/* iterate through all pairs */
loop (j = 1..nc2)
{
matout = ravel(matout, eval(cmd));
}
return(matout);
}