Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

draft: Use a SymPy Array not a Matrix for non-Expr #1194

Open
wants to merge 40 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
40 commits
Select commit Hold shift + click to select a range
880b600
Use a SymPy Array not a Matrix for non-Expr
cbm755 Oct 18, 2021
626c1a6
Python header: Add new function 'make_matrix_or_array'.
Aug 28, 2022
cdcc561
@sym/private/mat_rclist_access.m: Test Array-compatible code.
Aug 22, 2022
c18f830
@sym/private/mat_rclist_asgn.m: Test Array-compatible code.
Aug 24, 2022
10d26f2
@sym/{horzcat,vertcat}: Fix typo.
Aug 30, 2022
5f171e7
@sym/vertcat: Make function Array-compatible.
Aug 29, 2022
34c1f52
@sym/transpose: Make function Array-compatible.
Aug 30, 2022
6bfddcd
@sym/horzcat: Make function Array-compatible by simplifying it.
Aug 30, 2022
ce23c11
@sym/private/mat_rclist_{access,asgn}: Remove unused variables.
Aug 31, 2022
f95648e
@sym/repmat: Re-implement function without using 'pycall_sympy__'.
Sep 1, 2022
4a74fe2
elementwise_op: use 2d loop for Array
cbm755 Sep 2, 2022
0525456
private/python_ipc_native: Disable 'dbg_no_array' to keep in sync...
Sep 4, 2022
883326b
uniop_bool_helper: Array support in "bool" case
cbm755 Sep 2, 2022
5cb250e
logical: remove 1-d indexing for Array support
cbm755 Sep 4, 2022
d645f33
Minor fix to test: ctranspose of logical is tested elsewhere
cbm755 Sep 4, 2022
e3e0abc
Oops MatrixBase
cbm755 Sep 4, 2022
7386556
isAlways: avoid 1-d matrix/Array indexing
cbm755 Sep 4, 2022
af8caaa
levels=1 sufficient with tolist() from 2dArray/Matrix
cbm755 Sep 6, 2022
2ec6f02
subs: use make_matrix_or_array and use 2-d indexing
cbm755 Sep 3, 2022
978d8f0
Use S.Zero not integer 0
cbm755 Sep 5, 2022
535ecae
Array: fix a failing test in subasgn
cbm755 Sep 5, 2022
6cc099b
@sym/vertcat: Use sympy function 'flatten'.
Sep 6, 2022
f51d8ad
@sym/private/mat_rclist_asgn: Use sympy function 'flatten'.
Sep 6, 2022
84962e8
@sym/transpose: Use sympy function 'transpose'.
Sep 6, 2022
fd7efde
@sym/private/mat_replace: Make matrix mutable before modifying.
Sep 6, 2022
c154b6d
make_matrix_or_array: Generalise to accept an iterable of iterables.
Sep 6, 2022
25e0ebc
isconstant: fix behaviour with Array and add tests for errors
cbm755 Sep 6, 2022
60d5c12
find: work with Array
cbm755 Sep 6, 2022
6b02735
Fix some unused code
cbm755 Sep 6, 2022
d0ab2a3
elementwise_op: call make_array_or_matrix
cbm755 Sep 7, 2022
ea977c3
Rough gating of Array on older SymPy
cbm755 Sep 7, 2022
a90f041
Python header: Rename 'make_matrix_or_array' -> 'make_2d_sym'.
Aug 31, 2022
27b64eb
make_2d_sym: support flat input and shape
cbm755 Sep 13, 2022
174cbdc
Workaround nsolve with immutable x vars
cbm755 Sep 13, 2022
ee92b66
subs: after Array port was creating incorrect empties
cbm755 Sep 13, 2022
6f06abb
vertcat: tests for empty cases
cbm755 Sep 13, 2022
edcaec2
vertcat: use a flat list and shape for correct empty
cbm755 Sep 13, 2022
90af895
horzcat/vertcat: add tests for horzcat as well
cbm755 Sep 13, 2022
c9252ee
Port helpers from list-of-lists to flat+shape
cbm755 Sep 13, 2022
47834bd
isntance doesn't need special None handling
cbm755 Nov 12, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions inst/@sym/find.m
Original file line number Diff line number Diff line change
Expand Up @@ -132,8 +132,8 @@
' return r'
'#'
'x, = _ins'
'if x is not None and x.is_Matrix:'
' x = [a for a in x.T]' % note transpose
'if isinstance(x, (MatrixBase, NDimArray)):'
' x = [a for a in flatten(transpose(x).tolist(), levels=1)]' % note transpose
'else:'
' x = [x,]'
'return [scalar2tf(a) for a in x],' };
Expand Down
53 changes: 33 additions & 20 deletions inst/@sym/horzcat.m
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
%% SPDX-License-Identifier: GPL-3.0-or-later
%% Copyright (C) 2014-2017, 2019, 2024 Colin B. Macdonald
%% Copyright (C) 2014-2017, 2019, 2022, 2024 Colin B. Macdonald
%% Copyright (C) 2022 Alex Vong
%%
%% This file is part of OctSymPy.
%%
Expand Down Expand Up @@ -45,25 +46,8 @@

function h = horzcat(varargin)

% special case for 0x0 but other empties should be checked for
% compatibility
cmd = {
'_proc = []'
'for i in _ins:'
' if i is None or not i.is_Matrix:'
' _proc.append(sp.Matrix([[i]]))'
' else:'
' if i.shape == (0, 0):'
' pass'
' else:'
' _proc.append(i)'
'return sp.MatrixBase.hstack(*_proc),'
};

for i = 1:nargin
varargin{i} = sym(varargin{i});
end
h = pycall_sympy__ (cmd, varargin{:});
args = cellfun (@transpose, varargin, 'UniformOutput', false);
h = vertcat (args{:}).';
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

while I think of it: I approve of this kind of "surface decrease" (having fewer functions that call pycall_sympy__) BUT it is important to be aware that there is a downside: namely more "roundtrips" between Python and Octave. This is most noticeable with ipc-system but even with popen2, there are additional copies being made as everything is serialized back-and-forth over the pipe.

Its analogous to HTTP where all objects keep getting converted back-and-forth to JSON.

No action required but its something to keep in mind: perfect software design here might blow-up the latency and IO costs.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Interesting thoughts!

Normally when working with a single language, the classic solution is to use an (optimizing) compiler to optimize away all the redundant conversions. Unfortunately, this doesn't work for us because of language barrier: a compiler can't optimize across language boundary :(


Some of my thoughts:

First, pythonic may have already solved our problem! While I am not familiar with the internals of pythonic, perhaps using pythonic's py.... instead of pycall_sympy__ can avoid the needless conversions when going through a pipe. Of course, we should benchmark to see if it's indeed the case.

I believe this was the original future subproject but we changed it to the current (more urgent) one.


Another way is to create a octsmypy python library which contains one python function for each m-file function.

For example, say we want to implement m-file function @sym/transpose, @sym/vertcat and @sym/horzcat. Then our octsmypy python library should contain them as well. More precisely, we can do:

from octsympy.sym import transpose, vertcat, horzcat

To implement horzcat in terms of vertcat and transpose, we do it inside the octsmypy python library:

def horzcat(*XX):
    return transpose(vertcat(*[transpose(X) for X in XX]))

without the redundant conversions.

Finally, the m-file functions @sym/transpose, @sym/vertcat and @sym/horzcat are now thin wrappers over the corresponding octsympy python library functions, such as:

function h = vertcat (varargin)
    args = ...
    h = pycall_sympy__ ('from return octsympy.sym.vertcat(*_ins)', args{:});
end

and so on...

Perhaps this can be mentioned in the final report as possible future projects (?)

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree Pythonic offers possibilities. I think it currently still does conversion: a 'sym' is really just a struct of the 'srepr' (string representation) from SymPy, indep of IPC mechanism. But your 'xsym' idea would change that.

Agree about mentioing these things ad future projects.

Another thing to mention: CI pythonic runs are twice as slow as Popen2. I susect there migit be some 'low hanging fruit' in the Pythonic project.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And just for the record: another alternative that requires very little effort is (a) being mindful of this issue and aiming for 1 (or a small number) of language crosses per function (b) balancing that against technical debt of code duplication etc. This is the approach we've been using so far.


end

Expand Down Expand Up @@ -133,8 +117,37 @@
%! q = sym(ones(3, 0));
%! w = horzcat(v, q);

%!error <ShapeError>
%! z30 = sym (zeros (3, 0));
%! z40 = sym (zeros (4, 0));
%! % Note Issue #1238: unclear error message:
%! % [z30 z40];
%! % but this gives the ShapeError
%! horzcat(z30, z40);

%!test
%! % special case for the 0x0 empty: no error
%! z00 = sym (zeros (0, 0));
%! z30 = sym (zeros (3, 0));
%! [z00 z30];

%!test
%! % issue #700
%! A = sym ([1 2]);
%! B = simplify (A);
%! assert (isequal ([B A], [A B]))

%!test
%! % issue #1236, correct empty sizes
%! syms x
%! z00 = sym (zeros (0, 0));
%! z30 = sym (zeros (3, 0));
%! z03 = sym (zeros (0, 3));
%! z04 = sym (zeros (0, 4));
%! assert (size ([z00 z00]), [0 0])
%! assert (size ([z00 z03]), [0 3])
%! assert (size ([z03 z03]), [0 6])
%! assert (size ([z03 z04]), [0 7])
%! assert (size ([z03 z00 z04]), [0 7])
%! assert (size ([z30 z30]), [3 0])
%! assert (size ([z30 z30 z30]), [3 0])
4 changes: 2 additions & 2 deletions inst/@sym/isAlways.m
Original file line number Diff line number Diff line change
Expand Up @@ -128,8 +128,8 @@

cmd = vertcat(cmd, {
'(x, map_unknown_to) = _ins'
'if x is not None and x.is_Matrix:'
' r = [a for a in x.T]' % note transpose
'if isinstance(x, (MatrixBase, NDimArray)):'
' r = [a for a in flatten(transpose(x).tolist(), levels=1)]' % note tranpose

Check failure on line 132 in inst/@sym/isAlways.m

View workflow job for this annotation

GitHub Actions / codespell

tranpose ==> transpose

Check failure on line 132 in inst/@sym/isAlways.m

View workflow job for this annotation

GitHub Actions / codespell

tranpose ==> transpose
'else:'
' r = [x,]'
'r = [simplify_true_false_none(a) for a in r]'
Expand Down
13 changes: 12 additions & 1 deletion inst/@sym/isconstant.m
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@
function z = isconstant(x)

cmd = { '(x,) = _ins'
'if x is not None and x.is_Matrix:'
'if isinstance(x, (MatrixBase, NDimArray)):'
' return x.applyfunc(lambda a: a.is_constant()),'
'return x.is_constant(),' };
z = pycall_sympy__ (cmd, sym(x));
Expand All @@ -68,3 +68,14 @@
%! A = [x 2; 3 x];
%! B = [false true; true false];
%! assert (isequal (isconstant (A), B))

%!error
%! % semantically not an error, but is using SymPy's Array
%! t = sym(true);
%! A = [t t];
%! isconstant (A);

%!error
%! syms x
%! A = [x == 10; x <= 10];
%! isconstant (A);
8 changes: 4 additions & 4 deletions inst/@sym/logical.m
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
%% SPDX-License-Identifier: GPL-3.0-or-later
%% Copyright (C) 2014-2016, 2019, 2024 Colin B. Macdonald
%% Copyright (C) 2014-2016, 2019, 2022, 2024 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
Expand Down Expand Up @@ -102,8 +102,8 @@

cmd = vertcat(cmd, {
'(x, unknown) = _ins'
'if x is not None and x.is_Matrix:'
' r = [a for a in x.T]' % note transpose
'if isinstance(x, (MatrixBase, NDimArray)):'
' r = [a for a in flatten(transpose(x).tolist(), levels=1)]' % note tranpose

Check failure on line 106 in inst/@sym/logical.m

View workflow job for this annotation

GitHub Actions / codespell

tranpose ==> transpose

Check failure on line 106 in inst/@sym/logical.m

View workflow job for this annotation

GitHub Actions / codespell

tranpose ==> transpose
'else:'
' r = [x,]'
'r = [scalar2tfn(a) for a in r]'
Expand Down Expand Up @@ -180,7 +180,7 @@
%! w = logical(e);
%! assert (islogical (w))
%! assert (isequal (w, [true false true]))
%! e = e';
%! e = e.';
%! w = logical(e);
%! assert (islogical (w))
%! assert (isequal (w, [true; false; true]))
Expand Down
15 changes: 9 additions & 6 deletions inst/@sym/private/elementwise_op.m
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
%% Copyright (C) 2014, 2016, 2018-2019, 2022 Colin B. Macdonald
%% Copyright (C) 2014, 2016, 2018-2019, 2021-2022 Colin B. Macdonald
%% Copyright (C) 2016 Lagu
%%
%% This file is part of OctSymPy.
Expand Down Expand Up @@ -84,7 +84,7 @@
% Make sure all matrices in the input are the same size, and set q to one of them
'q = None'
'for A in _ins:'
' if isinstance(A, MatrixBase):'
' if isinstance(A, (MatrixBase, NDimArray)):'
' if q is None:'
' q = A'
' else:'
Expand All @@ -94,10 +94,13 @@
' return _op(*_ins)'
% at least one input was a matrix:
'# dbout(f"at least one matrix param, shape={q.shape}")'
'elements = []'
'for i in range(0, len(q)):'
' elements.append(_op(*[k[i] if isinstance(k, MatrixBase) else k for k in _ins]))'
'return Matrix(*q.shape, elements)' ];
'assert len(q.shape) == 2, "non-2D arrays/tensors not yet supported"'
'm, n = q.shape'
'g = []'
'for i in range(m):'
' for j in range(n):'
' g.append(_op(*[k[i, j] if isinstance(k, (MatrixBase, NDimArray)) else k for k in _ins]))'
'return _make_2d_sym(g, shape=q.shape)' ];

z = pycall_sympy__ (cmd, varargin{:});

Expand Down
15 changes: 6 additions & 9 deletions inst/@sym/private/mat_rclist_access.m
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
%% Copyright (C) 2014, 2016, 2019, 2022 Colin B. Macdonald
%% Copyright (C) 2022 Alex Vong
%%
%% This file is part of OctSymPy.
%%
Expand Down Expand Up @@ -32,15 +33,11 @@
error('this routine is for a list of rows and cols');
end

cmd = { '(A, rr, cc) = _ins'
'if A is None or not A.is_Matrix:'
' A = sp.Matrix([A])'
'n = len(rr)'
'M = [[0] for i in range(n)]'
'for i in range(0, n):'
' M[i][0] = A[rr[i],cc[i]]'
'M = sp.Matrix(M)'
'return M,' };
cmd = {'(A, rr, cc) = _ins'
'AA = A.tolist() if isinstance(A, (MatrixBase, NDimArray)) else [[A]]'
'M = [AA[i][j] for i, j in zip(rr, cc)]'
'return _make_2d_sym(M, shape=(len(rr), 1)),'
};

rr = num2cell(int32(r-1));
cc = num2cell(int32(c-1));
Expand Down
61 changes: 29 additions & 32 deletions inst/@sym/private/mat_rclist_asgn.m
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
%% Copyright (C) 2014, 2016-2017, 2019, 2022, 2024 Colin B. Macdonald
%% Copyright (C) 2020 Mike Miller
%% Copyright (C) 2020 Fernando Alvarruiz
%% Copyright (C) 2022 Alex Vong
%%
%% This file is part of OctSymPy.
%%
Expand Down Expand Up @@ -58,42 +59,38 @@
% Easy trick to copy A into larger matrix AA:
% AA = sp.Matrix.zeros(n, m)
% AA[0, 0] = A
% Also usefil: .copyin_matrix
% Also useful: .copyin_matrix

cmd = { '(A, r, c, B) = _ins'
'# B linear access fix, transpose for sympy row-based'
'if B is None or not B.is_Matrix:'
' B = sp.Matrix([[B]])'
'BT = B.T'
'# make a resized copy of A, and copy existing stuff in'
'if isinstance(A, list):'
' assert len(A) == 0, "unexpectedly non-empty list: report bug!"'
' n = max(max(r) + 1, 1)'
' m = max(max(c) + 1, 1)'
' AA = [[0]*m for i in range(n)]'
'elif A is None or not isinstance(A, MatrixBase):'
' # we have non-matrix, put in top-left'
' n = max(max(r) + 1, 1)'
' m = max(max(c) + 1, 1)'
' AA = [[0]*m for i in range(n)]'
' AA[0][0] = A'
'else:'
' # build bigger matrix'
' n = max(max(r) + 1, A.rows)'
' m = max(max(c) + 1, A.cols)'
' AA = [[0]*m for i in range(n)]'
' # copy current matrix in'
' for i in range(A.rows):'
' for j in range(A.cols):'
' AA[i][j] = A[i, j]'
'# now insert the new bits from B'
'for i, (r, c) in enumerate(zip(r, c)):'
' AA[r][c] = BT[i]'
'return sp.Matrix(AA),' };
cmd = {'(A, rr, cc, b) = _ins'
'assert A == [] or not isinstance(A, list), "unexpectedly non-empty list: report bug!"'
'if A == []:'
' AA = []'
' (nrows_A, ncols_A) = (0, 0)'
'elif isinstance(A, (MatrixBase, NDimArray)):'
' AA = A.tolist()'
' (nrows_A, ncols_A) = A.shape'
'else:'
' AA = [[A]]'
' (nrows_A, ncols_A) = (1, 1)'
'bb = b.tolist() if isinstance(b, (MatrixBase, NDimArray)) else [[b]]'
'entries = dict(zip(zip(rr, cc), flatten(bb, levels=1)))'
'def entry(i, j):'
' if (i, j) in entries:'
' return entries[i, j]'
' elif i < nrows_A and j < ncols_A:'
' return AA[i][j]'
' else:'
' return S.Zero'
'n = max(max(rr) + 1, nrows_A)'
'm = max(max(cc) + 1, ncols_A)'
'MM = [[entry(i, j) for j in range(m)] for i in range(n)]'
'M = make_2d_sym(MM)'
'return M,'};

rr = num2cell(int32(r-1));
cc = num2cell(int32(c-1));
z = pycall_sympy__ (cmd, A, rr, cc, B);
b = vec (B); # B is accessed with linear indexing, as a column vector
z = pycall_sympy__ (cmd, A, rr, cc, b);

% a simpler earlier version, but only for scalar r,c
%cmd = { '(A, r, c, b) = _ins'
Expand Down
7 changes: 5 additions & 2 deletions inst/@sym/private/mat_replace.m
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
%% Copyright (C) 2016 Abhinav Tripathi
%% Copyright (C) 2017 NVS Abhilash
%% Copyright (C) 2020 Fernando Alvarruiz
%% Copyright (C) 2022 Alex Vong
%%
%% This file is part of OctSymPy.
%%
Expand Down Expand Up @@ -161,7 +162,8 @@
if isscalar (A)
z = sym(zeros (1, 0));
else
cmd = { 'A, subs = _ins'
cmd = { 'AA, subs = _ins'
'A = AA.as_mutable()'
'if isinstance(subs, Integer):'
' A.col_del(subs - 1)'
' return A,'
Expand All @@ -178,7 +180,8 @@
% no test coverage: not sure how to hit this
z = sym(zeros (0, 1));
else
cmd = { 'A, subs = _ins'
cmd = { 'AA, subs = _ins'
'A = AA.as_mutable()'
'if isinstance(subs, Integer):'
' A.row_del(subs - 1)'
' return A,'
Expand Down
9 changes: 5 additions & 4 deletions inst/@sym/private/uniop_bool_helper.m
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
%% Copyright (C) 2016, 2019 Colin B. Macdonald
%% Copyright (C) 2016, 2019, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
Expand Down Expand Up @@ -51,9 +51,9 @@
cmd = [ cmd
'x = _ins[0]'
'pp = _ins[1:]'
'if x is not None and x.is_Matrix:'
'if isinstance(x, (MatrixBase, NDimArray)):'
' # bool will map None to False'
' return [bool(sf(a, *pp)) for a in x.T],'
' return [bool(sf(a, *pp)) for a in flatten(transpose(x).tolist(), levels=1)],'
'return bool(sf(x, *pp))' ];

r = pycall_sympy__ (cmd, x, varargin{:});
Expand All @@ -65,10 +65,11 @@
case 'sym'
warning('FIXME: not working for scalars')

% currently unused, and certainly not tested
cmd = [ cmd
'x = _ins[0]'
'pp = _ins[1:]'
'if x if not None and x.is_Matrix:'
'if isinstance(x, (MatrixBase, NDimArray)):'
' return x.applyfunc(sf, *pp)'
'return sf(x, *pp)' ];

Expand Down
34 changes: 23 additions & 11 deletions inst/@sym/repmat.m
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
%% Copyright (C) 2014, 2016, 2019 Colin B. Macdonald
%% Copyright (C) 2022 Alex Vong
%%
%% This file is part of OctSymPy.
%%
Expand Down Expand Up @@ -51,18 +52,21 @@
print_usage ();
end

cmd = { '(A, n, m) = _ins'
'if n == 0 or m == 0:'
' return sp.Matrix(n, m, [])'
'if A is None or not A.is_Matrix:'
' A = sp.Matrix([A])'
'L = [A]*m'
'B = sp.Matrix.hstack(*L)'
'L = [B]*n'
'B = sp.Matrix.vstack(*L)'
'return B' };
if n == 0 || m == 0
[nrows_A, ncols_A] = size (A);
B = sym (zeros (n * nrows_A, m * ncols_A));
return
end

%% duplicate A horizontally m times
MM = cellfun (@(varargin) A, cell (1, m), 'UniformOutput', false);
M = horzcat (MM{:});

%% now duplicate that result vertically n times
NN = cellfun (@(varargin) M, cell (1, n), 'UniformOutput', false);
N = vertcat (NN{:});

B = pycall_sympy__ (cmd, sym(A), int32(n), int32(m));
B = N;

end

Expand Down Expand Up @@ -95,3 +99,11 @@
%! assert (isequal (size(A), [0 3]))
%! A = repmat(sym(pi), [2 0]);
%! assert (isequal (size(A), [2 0]))

%!test
%! % even more empties, see also https://github.com/cbm755/octsympy/issues/1218
%! A = randi (16, 5, 7);
%! assert (isa (repmat (sym (A), [0 3]), 'sym'));
%! assert (isa (repmat (sym (A), [2 0]), 'sym'));
%! assert (isequal (sym (repmat (A, [0 3])), (repmat (sym (A), [0 3]))));
%! assert (isequal (sym (repmat (A, [2 0])), (repmat (sym (A), [2 0]))));
Loading
Loading