From 880b600eee5210cd3093218b30ef04cb8a7b326f Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Sun, 17 Oct 2021 22:05:56 -0700
Subject: [PATCH 01/40] Use a SymPy Array not a Matrix for non-Expr

It has mostly the same semantics as Matrix but can contain more general
things.  There might be some issues about different shapes of empties
though...  WIP on a fix for #1052.
---
 inst/@sym/private/elementwise_op.m | 9 ++++++---
 inst/private/python_header.py      | 5 ++++-
 2 files changed, 10 insertions(+), 4 deletions(-)

diff --git a/inst/@sym/private/elementwise_op.m b/inst/@sym/private/elementwise_op.m
index fe68ca95e..d03e0e94a 100644
--- a/inst/@sym/private/elementwise_op.m
+++ b/inst/@sym/private/elementwise_op.m
@@ -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.
@@ -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:'
@@ -97,7 +97,10 @@
           '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)' ];
+          'if all(isinstance(x, Expr) for x in elements):'
+          '    return Matrix(*q.shape, elements)'
+          'dbout(f"elementwise_op: returning an Array not a Matrix")'
+          'return Array(elements, shape=q.shape)' ];
 
   z = pycall_sympy__ (cmd, varargin{:});
 
diff --git a/inst/private/python_header.py b/inst/private/python_header.py
index 3506bf56a..f73e2e33b 100644
--- a/inst/private/python_header.py
+++ b/inst/private/python_header.py
@@ -155,7 +155,7 @@ def octoutput(x, et):
             f.text = str(OCTCODE_BOOL)
             f = ET.SubElement(a, "f")
             f.text = str(x)
-        elif x is None or isinstance(x, (sp.Basic, sp.MatrixBase)):
+        elif x is None or isinstance(x, (sp.Basic, sp.MatrixBase, sp.NDimArray)):
             # FIXME: is it weird to pretend None is a SymPy object?
             if isinstance(x, (sp.Matrix, sp.ImmutableMatrix)):
                 _d = x.shape
@@ -164,6 +164,9 @@ def octoutput(x, et):
                 _d = [float(r) if (isinstance(r, sp.Basic) and r.is_Integer)
                       else float('nan') if isinstance(r, sp.Basic)
                       else r for r in x.shape]
+            elif isinstance(x, sp.NDimArray):
+                _d = x.shape
+                dbout(f"I am here with an array with shape {_d}")
             elif x is None:
                 _d = (1,1)
             else:

From 626c1a6f0bc09bdc4ceac6e6137638854023535e Mon Sep 17 00:00:00 2001
From: Alex Vong <alexvong1995 AT protonmail DOT com>
Date: Sun, 28 Aug 2022 03:06:53 +0000
Subject: [PATCH 02/40] Python header: Add new function 'make_matrix_or_array'.

Fixes #1211.

* inst/private/python_header.py: Add it.
* inst/private/python_ipc_native.m: Add it.
---
 inst/private/python_header.py    | 24 ++++++++++++++++++++++++
 inst/private/python_ipc_native.m | 18 ++++++++++++++++++
 2 files changed, 42 insertions(+)

diff --git a/inst/private/python_header.py b/inst/private/python_header.py
index f73e2e33b..3882f572d 100644
--- a/inst/private/python_header.py
+++ b/inst/private/python_header.py
@@ -241,6 +241,30 @@ def octoutput(x, et):
 except:
     echo_exception_stdout("in python_header defining fcns block 4")
     raise
+
+
+try:
+    def make_matrix_or_array(ls_of_ls, dbg_no_array=True):
+        # should be kept in sync with the same function
+        # defined in inst/private/python_ipc_native.m
+        # FIXME: dbg_no_array is currently used for debugging,
+        # remove it once sympy drops non-Expr supp in Matrix
+        """
+        Given a list of list of syms LS_OF_LS
+        If all elements of LS_OF_LS are Expr,
+        construct the corresponding Matrix.
+        Otherwise, construct the corresponding 2D array.
+        """
+        elts = itertools.chain.from_iterable(ls_of_ls)
+        if (dbg_no_array
+            or all(isinstance(elt, Expr) for elt in elts)):
+            return Matrix(ls_of_ls)
+        else:
+            dbout(f"make_matrix_or_array: making 2D Array...")
+            return Array(ls_of_ls)
+except:
+    echo_exception_stdout("in python_header defining fcns block 5")
+    raise
 # end of python header, now couple blank lines
 
 
diff --git a/inst/private/python_ipc_native.m b/inst/private/python_ipc_native.m
index 5beb23cca..17a3f8822 100644
--- a/inst/private/python_ipc_native.m
+++ b/inst/private/python_ipc_native.m
@@ -112,6 +112,24 @@
                     '    # should be kept in sync with the same function'
                     '    # defined in inst/private/python_header.py'
                     '    sys.stderr.write("pydebug: " + str(l) + "\n")'
+                    'def make_matrix_or_array(ls_of_ls, dbg_no_array=True):'
+                    '    # should be kept in sync with the same function'
+                    '    # defined in inst/private/python_header.py'
+                    '    # FIXME: dbg_no_array is currently used for debugging,'
+                    '    # remove it once sympy drops non-Expr supp in Matrix'
+                    '    """'
+                    '    Given a list of list of syms LS_OF_LS'
+                    '    If all elements of LS_OF_LS are Expr,'
+                    '    construct the corresponding Matrix.'
+                    '    Otherwise, construct the corresponding 2D array.'
+                    '    """'
+                    '    elts = itertools.chain.from_iterable(ls_of_ls)'
+                    '    if (dbg_no_array'
+                    '        or all(isinstance(elt, Expr) for elt in elts)):'
+                    '        return Matrix(ls_of_ls)'
+                    '    else:'
+                    '        dbout(f"make_matrix_or_array: making 2D Array...")'
+                    '        return Array(ls_of_ls)'
                   }, newl))
     have_headers = true;
   end

From cdcc5616650a3653ea3a3066e5293473fa1d01e5 Mon Sep 17 00:00:00 2001
From: Alex Vong <alexvong1995 AT protonmail DOT com>
Date: Mon, 22 Aug 2022 12:23:59 +0000
Subject: [PATCH 03/40] @sym/private/mat_rclist_access.m: Test Array-compatible
 code.

See #1194 for more information.

* inst/@sym/private/mat_rclist_access.m: Test it.
---
 inst/@sym/private/mat_rclist_access.m | 16 +++++++---------
 1 file changed, 7 insertions(+), 9 deletions(-)

diff --git a/inst/@sym/private/mat_rclist_access.m b/inst/@sym/private/mat_rclist_access.m
index 46593030b..460ec99af 100644
--- a/inst/@sym/private/mat_rclist_access.m
+++ b/inst/@sym/private/mat_rclist_access.m
@@ -1,4 +1,5 @@
 %% Copyright (C) 2014, 2016, 2019, 2022 Colin B. Macdonald
+%% Copyright (C) 2022 Alex Vong
 %%
 %% This file is part of OctSymPy.
 %%
@@ -32,15 +33,12 @@
     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 = {'dbg_no_array = True'
+         '(A, rr, cc) = _ins'
+         'AA = A.tolist() if isinstance(A, (MatrixBase, NDimArray)) else [[A]]'
+         'MM = [[AA[i][j]] for i, j in zip(rr, cc)]'
+         'M = make_matrix_or_array(MM)'
+         'return M,'};
 
   rr = num2cell(int32(r-1));
   cc = num2cell(int32(c-1));

From c18f83089af326a6999a9ad04a62c6dcca59450f Mon Sep 17 00:00:00 2001
From: Alex Vong <alexvong1995 AT protonmail DOT com>
Date: Wed, 24 Aug 2022 04:29:56 +0000
Subject: [PATCH 04/40] @sym/private/mat_rclist_asgn.m: Test Array-compatible
 code.

See #1194 for more information.

* inst/@sym/private/mat_rclist_asgn.m: Test it.
---
 inst/@sym/private/mat_rclist_asgn.m | 61 ++++++++++++++---------------
 1 file changed, 30 insertions(+), 31 deletions(-)

diff --git a/inst/@sym/private/mat_rclist_asgn.m b/inst/@sym/private/mat_rclist_asgn.m
index 636cdb40c..be419ecec 100644
--- a/inst/@sym/private/mat_rclist_asgn.m
+++ b/inst/@sym/private/mat_rclist_asgn.m
@@ -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.
 %%
@@ -60,40 +61,38 @@
   %    AA[0, 0] = A
   % Also usefil: .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 = {'dbg_no_array = True'
+         '(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]]'
+         'bb_elts = itertools.chain.from_iterable(bb)'
+         'entries = dict(zip(zip(rr, cc), bb_elts))'
+         '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 0'
+         '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_matrix_or_array(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'

From 10d26f2f7425ab563e57c92d3db577fa1d065b5a Mon Sep 17 00:00:00 2001
From: Alex Vong <alexvong1995 AT protonmail DOT com>
Date: Tue, 30 Aug 2022 03:28:36 +0000
Subject: [PATCH 05/40] @sym/{horzcat,vertcat}: Fix typo.

* inst/@sym/{horzcat,vertcat}.m: Fix them.
---
 inst/@sym/horzcat.m | 1 +
 inst/@sym/vertcat.m | 1 +
 2 files changed, 2 insertions(+)

diff --git a/inst/@sym/horzcat.m b/inst/@sym/horzcat.m
index a31e2bb40..1a3eac869 100644
--- a/inst/@sym/horzcat.m
+++ b/inst/@sym/horzcat.m
@@ -1,5 +1,6 @@
 %% SPDX-License-Identifier: GPL-3.0-or-later
 %% Copyright (C) 2014-2017, 2019, 2024 Colin B. Macdonald
+%% Copyright (C) 2022 Alex Vong
 %%
 %% This file is part of OctSymPy.
 %%
diff --git a/inst/@sym/vertcat.m b/inst/@sym/vertcat.m
index 3a6907fa6..dbf0ba597 100644
--- a/inst/@sym/vertcat.m
+++ b/inst/@sym/vertcat.m
@@ -1,5 +1,6 @@
 %% SPDX-License-Identifier: GPL-3.0-or-later
 %% Copyright (C) 2014-2017, 2019, 2024 Colin B. Macdonald
+%% Copyright (C) 2022 Alex Vong
 %%
 %% This file is part of OctSymPy.
 %%

From 5f171e78efbcf83adf768cc30a193f80a4c710f4 Mon Sep 17 00:00:00 2001
From: Alex Vong <alexvong1995 AT protonmail DOT com>
Date: Mon, 29 Aug 2022 16:30:39 +0000
Subject: [PATCH 06/40] @sym/vertcat: Make function Array-compatible.

* inst/@sym/vertcat.m: Implement it.
---
 inst/@sym/vertcat.m | 36 +++++++++++++++++++-----------------
 1 file changed, 19 insertions(+), 17 deletions(-)

diff --git a/inst/@sym/vertcat.m b/inst/@sym/vertcat.m
index dbf0ba597..d225499b1 100644
--- a/inst/@sym/vertcat.m
+++ b/inst/@sym/vertcat.m
@@ -45,23 +45,25 @@
 
   % 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.vstack(*_proc),'
-          };
-
-  for i = 1:nargin
-    varargin{i} = sym(varargin{i});
-  end
-  h = pycall_sympy__ (cmd, varargin{:});
+  cmd = {'def is_matrix_or_array(x):'
+         '    return isinstance(x, (MatrixBase, NDimArray))'
+         'def number_of_columns(x):'
+         '    return x.shape[1] if is_matrix_or_array(x) else 1'
+         'def all_equal(*ls):'
+         '    return True if ls == [] else all(ls[0] == x for x in ls[1:])'
+         'def as_list_of_list(x):'
+         '    return x.tolist() if is_matrix_or_array(x) else [[x]]'
+         'args = [x for x in _ins if x != zeros(0, 0)] # remove 0x0 matrices'
+         'ncols = [number_of_columns(x) for x in args]'
+         'if not all_equal(*ncols):'
+         '    msg = "vertcat: all inputs must have the same number of columns"'
+         '    raise ShapeError(msg)'
+         'CCC = [as_list_of_list(x) for x in args]'
+         'CC = list(itertools.chain.from_iterable(CCC))'
+         'return make_matrix_or_array(CC)'};
+
+  args = cellfun (@sym, varargin, 'UniformOutput', false);
+  h = pycall_sympy__ (cmd, args{:});
 
 end
 

From 34c1f525d847a6e15d899370f802a5a9991b409d Mon Sep 17 00:00:00 2001
From: Alex Vong <alexvong1995 AT protonmail DOT com>
Date: Tue, 30 Aug 2022 15:03:53 +0000
Subject: [PATCH 07/40] @sym/transpose: Make function Array-compatible.

This is a prerequisite for making @sym/horzcat Array-compatible.

Partially fixes <https://github.com/cbm755/octsympy/issues/1215>.

* inst/@sym/transpose.m: Implement it.
---
 inst/@sym/transpose.m | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/inst/@sym/transpose.m b/inst/@sym/transpose.m
index 214766b35..8ae9beb31 100644
--- a/inst/@sym/transpose.m
+++ b/inst/@sym/transpose.m
@@ -1,4 +1,5 @@
 %% Copyright (C) 2014-2016, 2019 Colin B. Macdonald
+%% Copyright (C) 2022 Alex Vong
 %%
 %% This file is part of OctSymPy.
 %%
@@ -66,8 +67,11 @@
   end
 
   cmd = { 'x = _ins[0]'
-          'if x.is_Matrix:'
+          'if isinstance(x, MatrixBase):'
           '    return x.T'
+          'elif isinstance(x, NDimArray):'
+          '    xx = x.tolist()'
+          '    return make_matrix_or_array(list(zip(*xx)))'
           'else:'
           '    return x' };
 
@@ -92,3 +96,8 @@
 %!test
 %! A = [1 2] + 1i;
 %! assert(isequal( sym(A).' , sym(A.') ))
+
+%!test
+%! % see https://github.com/cbm755/octsympy/issues/1215
+%! none = pycall_sympy__ ('return None');
+%! assert (isequal (none.', none));

From 6bfddcd0642ecc1d71136adec10c2c706a1cd660 Mon Sep 17 00:00:00 2001
From: Alex Vong <alexvong1995 AT protonmail DOT com>
Date: Tue, 30 Aug 2022 13:44:02 +0000
Subject: [PATCH 08/40] @sym/horzcat: Make function Array-compatible by
 simplifying it.

* inst/@sym/horzcat.m: Simplify it.
---
 inst/@sym/horzcat.m | 21 ++-------------------
 1 file changed, 2 insertions(+), 19 deletions(-)

diff --git a/inst/@sym/horzcat.m b/inst/@sym/horzcat.m
index 1a3eac869..fed82f1cf 100644
--- a/inst/@sym/horzcat.m
+++ b/inst/@sym/horzcat.m
@@ -46,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{:}).';
 
 end
 

From ce23c11f0e4e9592060d75ac5cfd7a7006d6cd4d Mon Sep 17 00:00:00 2001
From: Alex Vong <alexvong1995 AT protonmail DOT com>
Date: Wed, 31 Aug 2022 13:43:44 +0000
Subject: [PATCH 09/40] @sym/private/mat_rclist_{access,asgn}: Remove unused
 variables.

* inst/@sym/private/mat_rclist_{access,asgn}.m: Remove them.
---
 inst/@sym/private/mat_rclist_access.m | 3 +--
 inst/@sym/private/mat_rclist_asgn.m   | 3 +--
 2 files changed, 2 insertions(+), 4 deletions(-)

diff --git a/inst/@sym/private/mat_rclist_access.m b/inst/@sym/private/mat_rclist_access.m
index 460ec99af..aea30b506 100644
--- a/inst/@sym/private/mat_rclist_access.m
+++ b/inst/@sym/private/mat_rclist_access.m
@@ -33,8 +33,7 @@
     error('this routine is for a list of rows and cols');
   end
 
-  cmd = {'dbg_no_array = True'
-         '(A, rr, cc) = _ins'
+  cmd = {'(A, rr, cc) = _ins'
          'AA = A.tolist() if isinstance(A, (MatrixBase, NDimArray)) else [[A]]'
          'MM = [[AA[i][j]] for i, j in zip(rr, cc)]'
          'M = make_matrix_or_array(MM)'
diff --git a/inst/@sym/private/mat_rclist_asgn.m b/inst/@sym/private/mat_rclist_asgn.m
index be419ecec..d73cd0893 100644
--- a/inst/@sym/private/mat_rclist_asgn.m
+++ b/inst/@sym/private/mat_rclist_asgn.m
@@ -61,8 +61,7 @@
   %    AA[0, 0] = A
   % Also usefil: .copyin_matrix
 
-  cmd = {'dbg_no_array = True'
-         '(A, rr, cc, b) = _ins'
+  cmd = {'(A, rr, cc, b) = _ins'
          'assert A == [] or not isinstance(A, list), "unexpectedly non-empty list: report bug!"'
          'if A == []:'
          '    AA = []'

From f95648e8d44b448d299f01bb0975b0d02e2ef871 Mon Sep 17 00:00:00 2001
From: Alex Vong <alexvong1995 AT protonmail DOT com>
Date: Thu, 1 Sep 2022 23:59:30 +0000
Subject: [PATCH 10/40] @sym/repmat: Re-implement function without using
 'pycall_sympy__'.

The function is now compatible with non-Matrix 2D sym (e.g. Array,
TableForm). Also, it now outputs empty matrices in the same way as
upstream Octave (fixes #1218).

* inst/@sym/repmat.m: Re-implement it and add tests accordingly.
---
 inst/@sym/repmat.m | 34 +++++++++++++++++++++++-----------
 1 file changed, 23 insertions(+), 11 deletions(-)

diff --git a/inst/@sym/repmat.m b/inst/@sym/repmat.m
index af16d0d8f..7a8d521fe 100644
--- a/inst/@sym/repmat.m
+++ b/inst/@sym/repmat.m
@@ -1,4 +1,5 @@
 %% Copyright (C) 2014, 2016, 2019 Colin B. Macdonald
+%% Copyright (C) 2022 Alex Vong
 %%
 %% This file is part of OctSymPy.
 %%
@@ -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
 
@@ -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]))));

From 4a74fe217e9010c60031dfcb422a3a5f810db4ef Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Fri, 2 Sep 2022 11:09:47 -0700
Subject: [PATCH 11/40] elementwise_op: use 2d loop for Array

Array linear indexing is not the same as Matrix linear indexing.
Fixes #1222.
---
 inst/@sym/private/elementwise_op.m | 6 ++++--
 inst/private/python_header.py      | 2 +-
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/inst/@sym/private/elementwise_op.m b/inst/@sym/private/elementwise_op.m
index d03e0e94a..441cc6d85 100644
--- a/inst/@sym/private/elementwise_op.m
+++ b/inst/@sym/private/elementwise_op.m
@@ -95,8 +95,10 @@
           % 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]))'
+	  'assert len(q.shape) == 2, "non-2D arrays/tensors not yet supported"'
+          'for i in range(0, q.shape[0]):'
+          '    for j in range(0, q.shape[1]):'
+          '        elements.append(_op(*[k[i, j] if isinstance(k, (MatrixBase, NDimArray)) else k for k in _ins]))'
           'if all(isinstance(x, Expr) for x in elements):'
           '    return Matrix(*q.shape, elements)'
           'dbout(f"elementwise_op: returning an Array not a Matrix")'
diff --git a/inst/private/python_header.py b/inst/private/python_header.py
index 3882f572d..baa18a391 100644
--- a/inst/private/python_header.py
+++ b/inst/private/python_header.py
@@ -244,7 +244,7 @@ def octoutput(x, et):
 
 
 try:
-    def make_matrix_or_array(ls_of_ls, dbg_no_array=True):
+    def make_matrix_or_array(ls_of_ls, dbg_no_array=False):
         # should be kept in sync with the same function
         # defined in inst/private/python_ipc_native.m
         # FIXME: dbg_no_array is currently used for debugging,

From 052545646ffb3bf2052c65a3be8462b3b6712c15 Mon Sep 17 00:00:00 2001
From: Alex Vong <alexvong1995 AT protonmail DOT com>
Date: Sun, 4 Sep 2022 20:01:25 +0000
Subject: [PATCH 12/40] private/python_ipc_native: Disable 'dbg_no_array' to
 keep in sync...

...with 'inst/private/python_header.py'.

Follow-up to 56bd71902ca0faf5fa88c1bc7fbeba8bf6f99ce7.

* inst/private/python_ipc_native.m: Disable it.
---
 inst/private/python_ipc_native.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/inst/private/python_ipc_native.m b/inst/private/python_ipc_native.m
index 17a3f8822..453307230 100644
--- a/inst/private/python_ipc_native.m
+++ b/inst/private/python_ipc_native.m
@@ -112,7 +112,7 @@
                     '    # should be kept in sync with the same function'
                     '    # defined in inst/private/python_header.py'
                     '    sys.stderr.write("pydebug: " + str(l) + "\n")'
-                    'def make_matrix_or_array(ls_of_ls, dbg_no_array=True):'
+                    'def make_matrix_or_array(ls_of_ls, dbg_no_array=False):'
                     '    # should be kept in sync with the same function'
                     '    # defined in inst/private/python_header.py'
                     '    # FIXME: dbg_no_array is currently used for debugging,'

From 883326b532d969f1a15592fd7359abf0bc533239 Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Fri, 2 Sep 2022 13:16:12 -0700
Subject: [PATCH 13/40] uniop_bool_helper: Array support in "bool" case

This routine is supposed to return a logical Octave array.  Do linear
indexing on the flattened tolist output which will work on both Array
and Matrix.  Fixes #1221.
---
 inst/@sym/private/uniop_bool_helper.m | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/inst/@sym/private/uniop_bool_helper.m b/inst/@sym/private/uniop_bool_helper.m
index c9d6ab201..03cc064d8 100644
--- a/inst/@sym/private/uniop_bool_helper.m
+++ b/inst/@sym/private/uniop_bool_helper.m
@@ -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.
 %%
@@ -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())],'
               'return bool(sf(x, *pp))' ];
 
       r = pycall_sympy__ (cmd, x, varargin{:});

From 5cb250e8ee39edf1f43b79970260f68534c2cf3b Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Sat, 3 Sep 2022 20:07:11 -0700
Subject: [PATCH 14/40] logical: remove 1-d indexing for Array support

---
 inst/@sym/logical.m | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/inst/@sym/logical.m b/inst/@sym/logical.m
index d664defd8..95106f4e5 100644
--- a/inst/@sym/logical.m
+++ b/inst/@sym/logical.m
@@ -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.
 %%
@@ -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, (Matrix, NDimArray)):'
+    '    r = [a for a in flatten(transpose(x).tolist())]'  % note tranpose
     'else:'
     '    r = [x,]'
     'r = [scalar2tfn(a) for a in r]'

From d645f33359e3bb4904288ae58b3122312d7059a5 Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Sat, 3 Sep 2022 20:13:00 -0700
Subject: [PATCH 15/40] Minor fix to test: ctranspose of logical is tested
 elsewhere

In fact it is failing elsewhere but the point of *this* test is to
change the orientation of a vector.
---
 inst/@sym/logical.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/inst/@sym/logical.m b/inst/@sym/logical.m
index 95106f4e5..bb86c24a9 100644
--- a/inst/@sym/logical.m
+++ b/inst/@sym/logical.m
@@ -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]))

From e3e0abc76f9065c0011e3652e509a2b184467cd9 Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Sat, 3 Sep 2022 20:35:51 -0700
Subject: [PATCH 16/40] Oops MatrixBase

---
 inst/@sym/logical.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/inst/@sym/logical.m b/inst/@sym/logical.m
index bb86c24a9..42a8dea3b 100644
--- a/inst/@sym/logical.m
+++ b/inst/@sym/logical.m
@@ -102,7 +102,7 @@
 
   cmd = vertcat(cmd, {
     '(x, unknown) = _ins'
-    'if isinstance(x, (Matrix, NDimArray)):'
+    'if isinstance(x, (MatrixBase, NDimArray)):'
     '    r = [a for a in flatten(transpose(x).tolist())]'  % note tranpose
     'else:'
     '    r = [x,]'

From 7386556193a68fac2bea4a7d5bfc6daf8f77208e Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Sat, 3 Sep 2022 20:36:05 -0700
Subject: [PATCH 17/40] isAlways: avoid 1-d matrix/Array indexing

---
 inst/@sym/isAlways.m | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/inst/@sym/isAlways.m b/inst/@sym/isAlways.m
index 08989b5e9..146a759ed 100644
--- a/inst/@sym/isAlways.m
+++ b/inst/@sym/isAlways.m
@@ -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 x is not None and isinstance(x, (MatrixBase, NDimArray)):'
+    '    r = [a for a in flatten(transpose(x).tolist())]'  % note tranpose
     'else:'
     '    r = [x,]'
     'r = [simplify_true_false_none(a) for a in r]'

From af8caaa66d7b8e74019fc9d6a57616181caab277 Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Mon, 5 Sep 2022 17:19:10 -0700
Subject: [PATCH 18/40] levels=1 sufficient with tolist() from 2dArray/Matrix

---
 inst/@sym/isAlways.m                  | 2 +-
 inst/@sym/logical.m                   | 2 +-
 inst/@sym/private/uniop_bool_helper.m | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/inst/@sym/isAlways.m b/inst/@sym/isAlways.m
index 146a759ed..b27d5db6a 100644
--- a/inst/@sym/isAlways.m
+++ b/inst/@sym/isAlways.m
@@ -129,7 +129,7 @@
   cmd = vertcat(cmd, {
     '(x, map_unknown_to) = _ins'
     'if x is not None and isinstance(x, (MatrixBase, NDimArray)):'
-    '    r = [a for a in flatten(transpose(x).tolist())]'  % note tranpose
+    '    r = [a for a in flatten(transpose(x).tolist(), levels=1)]'  % note tranpose
     'else:'
     '    r = [x,]'
     'r = [simplify_true_false_none(a) for a in r]'
diff --git a/inst/@sym/logical.m b/inst/@sym/logical.m
index 42a8dea3b..1fa56eabf 100644
--- a/inst/@sym/logical.m
+++ b/inst/@sym/logical.m
@@ -103,7 +103,7 @@
   cmd = vertcat(cmd, {
     '(x, unknown) = _ins'
     'if isinstance(x, (MatrixBase, NDimArray)):'
-    '    r = [a for a in flatten(transpose(x).tolist())]'  % note tranpose
+    '    r = [a for a in flatten(transpose(x).tolist(), levels=1)]'  % note tranpose
     'else:'
     '    r = [x,]'
     'r = [scalar2tfn(a) for a in r]'
diff --git a/inst/@sym/private/uniop_bool_helper.m b/inst/@sym/private/uniop_bool_helper.m
index 03cc064d8..bb81e0207 100644
--- a/inst/@sym/private/uniop_bool_helper.m
+++ b/inst/@sym/private/uniop_bool_helper.m
@@ -53,7 +53,7 @@
               'pp = _ins[1:]'
               'if isinstance(x, (MatrixBase, NDimArray)):'
               '    # bool will map None to False'
-              '    return [bool(sf(a, *pp)) for a in flatten(transpose(x).tolist())],'
+              '    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{:});

From 2ec6f026e3dde3d56f872b172b321fa32a6c81b3 Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Sat, 3 Sep 2022 16:28:24 -0700
Subject: [PATCH 19/40] subs: use make_matrix_or_array and use 2-d indexing

Fixes Issue #1226.  Related to Issue #1222.
---
 inst/@sym/subs.m | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/inst/@sym/subs.m b/inst/@sym/subs.m
index 66d786163..132291199 100644
--- a/inst/@sym/subs.m
+++ b/inst/@sym/subs.m
@@ -255,14 +255,14 @@
     'sizes = {(a.shape if a.is_Matrix else (1, 1)) for a in yy}'
     'sizes.discard((1, 1))'
     'assert len(sizes) == 1, "all substitutions must be same size or scalar"'
-    'size = sizes.pop()'
-    'numel = prod(size)'
-    'g = [0]*numel'
-    'for i in range(len(g)):'
-    '    yyy = [y[i] if y.is_Matrix else y for y in yy]'
-    '    sublist = list(zip(xx, yyy))'
-    '    g[i] = f.subs(sublist, simultaneous=True).doit()'
-    'return Matrix(*size, g)'
+    'm, n = sizes.pop()'
+    'g = [[0]*n for i in range(m)]'
+    'for i in range(m):'
+    '    for j in range(n):'
+    '        yyy = [y[i, j] if y.is_Matrix else y for y in yy]'
+    '        sublist = list(zip(xx, yyy))'
+    '        g[i][j] = f.subs(sublist, simultaneous=True).doit()'
+    'return make_matrix_or_array(g)'
   };
 
   g = pycall_sympy__ (cmd, f, in, out);

From 978d8f04f7ea306e1639d5fbb7d672fe3f2555c9 Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Sun, 4 Sep 2022 21:53:14 -0700
Subject: [PATCH 20/40] Use S.Zero not integer 0

Fixes #1228.
---
 inst/@sym/private/mat_rclist_asgn.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/inst/@sym/private/mat_rclist_asgn.m b/inst/@sym/private/mat_rclist_asgn.m
index d73cd0893..a3841bf70 100644
--- a/inst/@sym/private/mat_rclist_asgn.m
+++ b/inst/@sym/private/mat_rclist_asgn.m
@@ -81,7 +81,7 @@
          '    elif i < nrows_A and j < ncols_A:'
          '        return AA[i][j]'
          '    else:'
-         '        return 0'
+         '        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)]'

From 535ecaed82b8c0f2ce158a20f5b5f75b4b27b791 Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Sun, 4 Sep 2022 22:07:10 -0700
Subject: [PATCH 21/40] Array: fix a failing test in subasgn

---
 inst/@sym/reshape.m | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/inst/@sym/reshape.m b/inst/@sym/reshape.m
index abbdfdb28..b80879193 100644
--- a/inst/@sym/reshape.m
+++ b/inst/@sym/reshape.m
@@ -1,4 +1,4 @@
-%% Copyright (C) 2014, 2016, 2019 Colin B. Macdonald
+%% Copyright (C) 2014, 2016, 2019, 2022 Colin B. Macdonald
 %%
 %% This file is part of OctSymPy.
 %%
@@ -75,9 +75,10 @@
 
   cmd = {
       '(A, n, m) = _ins'
-      'if A is not None and A.is_Matrix:'
-      '    #sympy is row-based'
-      '    return A.T.reshape(m,n).T'
+      'if isinstance(A, (MatrixBase, NDimArray)):'
+      '    #sympy is row-based, but Array does not have .T'
+      '    #return A.T.reshape(m,n).T'
+      '    return transpose(transpose(A).reshape(m, n))'
       'else:'
       '    if n != 1 or m != 1:'
       '        raise ValueError("cannot reshape scalar to non-1x1 size")'

From 6cc099bdc81e54d671d8f15faed99c5cc20f6ea1 Mon Sep 17 00:00:00 2001
From: Alex Vong <alexvong1995 AT protonmail DOT com>
Date: Tue, 6 Sep 2022 15:48:33 +0000
Subject: [PATCH 22/40] @sym/vertcat: Use sympy function 'flatten'.

Follow-up to d3b154731803c95460af08722e7a62a718d20870.

* inst/@sym/vertcat.m: Use it.
---
 inst/@sym/vertcat.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/inst/@sym/vertcat.m b/inst/@sym/vertcat.m
index d225499b1..23a66087c 100644
--- a/inst/@sym/vertcat.m
+++ b/inst/@sym/vertcat.m
@@ -59,7 +59,7 @@
          '    msg = "vertcat: all inputs must have the same number of columns"'
          '    raise ShapeError(msg)'
          'CCC = [as_list_of_list(x) for x in args]'
-         'CC = list(itertools.chain.from_iterable(CCC))'
+         'CC = flatten(CCC, levels=1)'
          'return make_matrix_or_array(CC)'};
 
   args = cellfun (@sym, varargin, 'UniformOutput', false);

From f51d8ad4aa0f0415d1ad437ecd340094fcd65114 Mon Sep 17 00:00:00 2001
From: Alex Vong <alexvong1995 AT protonmail DOT com>
Date: Tue, 6 Sep 2022 19:01:33 +0000
Subject: [PATCH 23/40] @sym/private/mat_rclist_asgn: Use sympy function
 'flatten'.

* inst/@sym/private/mat_rclist_asgn.m: Use it.
---
 inst/@sym/private/mat_rclist_asgn.m | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/inst/@sym/private/mat_rclist_asgn.m b/inst/@sym/private/mat_rclist_asgn.m
index a3841bf70..836714519 100644
--- a/inst/@sym/private/mat_rclist_asgn.m
+++ b/inst/@sym/private/mat_rclist_asgn.m
@@ -73,8 +73,7 @@
          '    AA = [[A]]'
          '    (nrows_A, ncols_A) = (1, 1)'
          'bb = b.tolist() if isinstance(b, (MatrixBase, NDimArray)) else [[b]]'
-         'bb_elts = itertools.chain.from_iterable(bb)'
-         'entries = dict(zip(zip(rr, cc), bb_elts))'
+         'entries = dict(zip(zip(rr, cc), flatten(bb, levels=1)))'
          'def entry(i, j):'
          '    if (i, j) in entries:'
          '        return entries[i, j]'

From 84962e80abc5c1bf85ed5cb07eafab2779c2abf7 Mon Sep 17 00:00:00 2001
From: Alex Vong <alexvong1995 AT protonmail DOT com>
Date: Tue, 6 Sep 2022 14:46:06 +0000
Subject: [PATCH 24/40] @sym/transpose: Use sympy function 'transpose'.

Follow-up to 136184d6df506bc9945d8f341e9b46691895c5b4.

* inst/@sym/transpose.m: Use it.
---
 inst/@sym/transpose.m | 12 ++++--------
 1 file changed, 4 insertions(+), 8 deletions(-)

diff --git a/inst/@sym/transpose.m b/inst/@sym/transpose.m
index 8ae9beb31..161554803 100644
--- a/inst/@sym/transpose.m
+++ b/inst/@sym/transpose.m
@@ -66,14 +66,10 @@
     print_usage ();
   end
 
-  cmd = { 'x = _ins[0]'
-          'if isinstance(x, MatrixBase):'
-          '    return x.T'
-          'elif isinstance(x, NDimArray):'
-          '    xx = x.tolist()'
-          '    return make_matrix_or_array(list(zip(*xx)))'
-          'else:'
-          '    return x' };
+  cmd = {'def is_matrix_or_array(x):'
+         '    return isinstance(x, (MatrixBase, NDimArray))'
+         'x, = _ins'
+         'return transpose(x) if is_matrix_or_array(x) else x'};
 
   z = pycall_sympy__ (cmd, x);
 

From fd7efde08a1e60524609128ba258e9a7d30f4f9f Mon Sep 17 00:00:00 2001
From: Alex Vong <alexvong1995 AT protonmail DOT com>
Date: Tue, 6 Sep 2022 20:52:34 +0000
Subject: [PATCH 25/40] @sym/private/mat_replace: Make matrix mutable before
 modifying.

We did not run into issues in the past because many sympy functions
return a mutable matrix. However, 'transpose' used in
c9f6f0f9b1d0edaee5c29ce7d9e4ae2ece8708d4 returns an immutable
matrix. This change should avoid similar problems in the future.

Follow-up to c9f6f0f9b1d0edaee5c29ce7d9e4ae2ece8708d4.

* inst/@sym/private/mat_replace.m: Make it mutable.
---
 inst/@sym/private/mat_replace.m | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/inst/@sym/private/mat_replace.m b/inst/@sym/private/mat_replace.m
index 9919a0fda..931c9e9d6 100644
--- a/inst/@sym/private/mat_replace.m
+++ b/inst/@sym/private/mat_replace.m
@@ -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.
 %%
@@ -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,'
@@ -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,'

From c154b6d57bcad060d48ceaa7c2a39bd59abdcd1b Mon Sep 17 00:00:00 2001
From: Alex Vong <alexvong1995 AT protonmail DOT com>
Date: Tue, 6 Sep 2022 17:49:46 +0000
Subject: [PATCH 26/40] make_matrix_or_array: Generalise to accept an iterable
 of iterables.

* inst/private/{python_header.py,python_ipc_native.m}: Generalise
function 'make_matrix_or_array', use sympy function 'flatten' in
'make_matrix_or_array'.
---
 inst/private/python_header.py    | 9 +++++----
 inst/private/python_ipc_native.m | 9 +++++----
 2 files changed, 10 insertions(+), 8 deletions(-)

diff --git a/inst/private/python_header.py b/inst/private/python_header.py
index baa18a391..ad1554262 100644
--- a/inst/private/python_header.py
+++ b/inst/private/python_header.py
@@ -244,18 +244,19 @@ def octoutput(x, et):
 
 
 try:
-    def make_matrix_or_array(ls_of_ls, dbg_no_array=False):
+    def make_matrix_or_array(it_of_it, dbg_no_array=False):
         # should be kept in sync with the same function
         # defined in inst/private/python_ipc_native.m
         # FIXME: dbg_no_array is currently used for debugging,
         # remove it once sympy drops non-Expr supp in Matrix
         """
-        Given a list of list of syms LS_OF_LS
-        If all elements of LS_OF_LS are Expr,
+        Given an iterable of iterables of syms IT_OF_IT
+        If all elements of IT_OF_IT are Expr,
         construct the corresponding Matrix.
         Otherwise, construct the corresponding 2D array.
         """
-        elts = itertools.chain.from_iterable(ls_of_ls)
+        ls_of_ls = [[elt for elt in it] for it in it_of_it]
+        elts = flatten(ls_of_ls, levels=1)
         if (dbg_no_array
             or all(isinstance(elt, Expr) for elt in elts)):
             return Matrix(ls_of_ls)
diff --git a/inst/private/python_ipc_native.m b/inst/private/python_ipc_native.m
index 453307230..31a26f684 100644
--- a/inst/private/python_ipc_native.m
+++ b/inst/private/python_ipc_native.m
@@ -112,18 +112,19 @@
                     '    # should be kept in sync with the same function'
                     '    # defined in inst/private/python_header.py'
                     '    sys.stderr.write("pydebug: " + str(l) + "\n")'
-                    'def make_matrix_or_array(ls_of_ls, dbg_no_array=False):'
+                    'def make_matrix_or_array(it_of_it, dbg_no_array=False):'
                     '    # should be kept in sync with the same function'
                     '    # defined in inst/private/python_header.py'
                     '    # FIXME: dbg_no_array is currently used for debugging,'
                     '    # remove it once sympy drops non-Expr supp in Matrix'
                     '    """'
-                    '    Given a list of list of syms LS_OF_LS'
-                    '    If all elements of LS_OF_LS are Expr,'
+                    '    Given an iterable of iterables of syms IT_OF_IT'
+                    '    If all elements of IT_OF_IT are Expr,'
                     '    construct the corresponding Matrix.'
                     '    Otherwise, construct the corresponding 2D array.'
                     '    """'
-                    '    elts = itertools.chain.from_iterable(ls_of_ls)'
+                    '    ls_of_ls = [[elt for elt in it] for it in it_of_it]'
+                    '    elts = flatten(ls_of_ls, levels=1)'
                     '    if (dbg_no_array'
                     '        or all(isinstance(elt, Expr) for elt in elts)):'
                     '        return Matrix(ls_of_ls)'

From 25e0ebc8dced147ab5ed309c781b80b6e16410d2 Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Mon, 5 Sep 2022 17:15:22 -0700
Subject: [PATCH 27/40] isconstant: fix behaviour with Array and add tests for
 errors

---
 inst/@sym/isconstant.m | 13 ++++++++++++-
 1 file changed, 12 insertions(+), 1 deletion(-)

diff --git a/inst/@sym/isconstant.m b/inst/@sym/isconstant.m
index fc580a4a6..ccd900f45 100644
--- a/inst/@sym/isconstant.m
+++ b/inst/@sym/isconstant.m
@@ -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));
@@ -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);

From 60d5c1213fc99d7184822694a9c878388870f6a8 Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Mon, 5 Sep 2022 21:22:08 -0700
Subject: [PATCH 28/40] find: work with Array

---
 inst/@sym/find.m | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/inst/@sym/find.m b/inst/@sym/find.m
index fde8a9fa4..af7b58850 100644
--- a/inst/@sym/find.m
+++ b/inst/@sym/find.m
@@ -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],' };

From 6b0273529b2ab52f4bbc80fa16b31a3f1321c405 Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Mon, 5 Sep 2022 22:13:08 -0700
Subject: [PATCH 29/40] Fix some unused code

This was broken with syntax errors since 2016...  Apply fix for Array
but leave a comment that this is unused and untested.  It is referred to
in isprime but commented out.
---
 inst/@sym/private/uniop_bool_helper.m | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/inst/@sym/private/uniop_bool_helper.m b/inst/@sym/private/uniop_bool_helper.m
index bb81e0207..b2f3eb927 100644
--- a/inst/@sym/private/uniop_bool_helper.m
+++ b/inst/@sym/private/uniop_bool_helper.m
@@ -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)' ];
 

From d0ab2a388584d17819a0d392e47e14f5f5b2ff64 Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Tue, 6 Sep 2022 17:28:09 -0700
Subject: [PATCH 30/40] elementwise_op: call make_array_or_matrix

Instead of doing it ourselves, prepare a list of lists and call the
centralized helper function.
---
 inst/@sym/private/elementwise_op.m | 16 +++++++---------
 1 file changed, 7 insertions(+), 9 deletions(-)

diff --git a/inst/@sym/private/elementwise_op.m b/inst/@sym/private/elementwise_op.m
index 441cc6d85..fb324005f 100644
--- a/inst/@sym/private/elementwise_op.m
+++ b/inst/@sym/private/elementwise_op.m
@@ -94,15 +94,13 @@
           '    return _op(*_ins)'
           % at least one input was a matrix:
           '# dbout(f"at least one matrix param, shape={q.shape}")'
-          'elements = []'
-	  'assert len(q.shape) == 2, "non-2D arrays/tensors not yet supported"'
-          'for i in range(0, q.shape[0]):'
-          '    for j in range(0, q.shape[1]):'
-          '        elements.append(_op(*[k[i, j] if isinstance(k, (MatrixBase, NDimArray)) else k for k in _ins]))'
-          'if all(isinstance(x, Expr) for x in elements):'
-          '    return Matrix(*q.shape, elements)'
-          'dbout(f"elementwise_op: returning an Array not a Matrix")'
-          'return Array(elements, shape=q.shape)' ];
+          'assert len(q.shape) == 2, "non-2D arrays/tensors not yet supported"'
+          'm, n = q.shape'
+          'g = [[0]*n for i in range(m)]'
+          'for i in range(m):'
+          '    for j in range(n):'
+          '        g[i][j] = _op(*[k[i, j] if isinstance(k, (MatrixBase, NDimArray)) else k for k in _ins])'
+          'return make_matrix_or_array(g)' ];
 
   z = pycall_sympy__ (cmd, varargin{:});
 

From ea977c3a238ee93adbadcbd93bcfcd46fdcb630b Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Tue, 6 Sep 2022 17:48:37 -0700
Subject: [PATCH 31/40] Rough gating of Array on older SymPy

---
 inst/private/python_header.py    | 3 +++
 inst/private/python_ipc_native.m | 3 +++
 2 files changed, 6 insertions(+)

diff --git a/inst/private/python_header.py b/inst/private/python_header.py
index ad1554262..df855dd15 100644
--- a/inst/private/python_header.py
+++ b/inst/private/python_header.py
@@ -257,6 +257,9 @@ def make_matrix_or_array(it_of_it, dbg_no_array=False):
         """
         ls_of_ls = [[elt for elt in it] for it in it_of_it]
         elts = flatten(ls_of_ls, levels=1)
+        if Version(spver) <= Version("1.11.1"):
+            # never use Array on older SymPy
+            dbg_no_array = True
         if (dbg_no_array
             or all(isinstance(elt, Expr) for elt in elts)):
             return Matrix(ls_of_ls)
diff --git a/inst/private/python_ipc_native.m b/inst/private/python_ipc_native.m
index 31a26f684..956537b52 100644
--- a/inst/private/python_ipc_native.m
+++ b/inst/private/python_ipc_native.m
@@ -125,6 +125,9 @@
                     '    """'
                     '    ls_of_ls = [[elt for elt in it] for it in it_of_it]'
                     '    elts = flatten(ls_of_ls, levels=1)'
+                    '    if Version(spver) <= Version("1.11.1"):'
+                    '        # never use Array on older SymPy'
+                    '        dbg_no_array = True'
                     '    if (dbg_no_array'
                     '        or all(isinstance(elt, Expr) for elt in elts)):'
                     '        return Matrix(ls_of_ls)'

From a90f041e663170a323fbcb8b0c04cfc782731b04 Mon Sep 17 00:00:00 2001
From: Alex Vong <alexvong1995 AT protonmail DOT com>
Date: Wed, 31 Aug 2022 13:40:56 +0000
Subject: [PATCH 32/40] Python header: Rename 'make_matrix_or_array' ->
 'make_2d_sym'.

Generalises 'make_matrix_or_array' so that we can choose to represent
non-Matrix 2D sym to be something other than an Array in the
future. For example, we could use TableForm.

* inst/private/{python_header.py,python_ipc_native.m}: Rename function
'make_matrix_or_array' -> 'make_2d_sym', variable 'dbg_no_array' ->
'dbg_matrix_only' and adjust accordingly.
* inst/@sym/private/mat_rclist_{access,asgn}.m: Adjust accordingly.
* inst/@sym/vertcat.m: Adjust accordingly.
---
 inst/@sym/private/elementwise_op.m    |  2 +-
 inst/@sym/private/mat_rclist_access.m |  2 +-
 inst/@sym/private/mat_rclist_asgn.m   |  2 +-
 inst/@sym/subs.m                      |  2 +-
 inst/@sym/vertcat.m                   |  2 +-
 inst/private/python_header.py         | 20 ++++++++++----------
 inst/private/python_ipc_native.m      | 20 ++++++++++----------
 7 files changed, 25 insertions(+), 25 deletions(-)

diff --git a/inst/@sym/private/elementwise_op.m b/inst/@sym/private/elementwise_op.m
index fb324005f..996bd43d5 100644
--- a/inst/@sym/private/elementwise_op.m
+++ b/inst/@sym/private/elementwise_op.m
@@ -100,7 +100,7 @@
           'for i in range(m):'
           '    for j in range(n):'
           '        g[i][j] = _op(*[k[i, j] if isinstance(k, (MatrixBase, NDimArray)) else k for k in _ins])'
-          'return make_matrix_or_array(g)' ];
+          'return make_2d_sym(g)' ];
 
   z = pycall_sympy__ (cmd, varargin{:});
 
diff --git a/inst/@sym/private/mat_rclist_access.m b/inst/@sym/private/mat_rclist_access.m
index aea30b506..c46036a31 100644
--- a/inst/@sym/private/mat_rclist_access.m
+++ b/inst/@sym/private/mat_rclist_access.m
@@ -36,7 +36,7 @@
   cmd = {'(A, rr, cc) = _ins'
          'AA = A.tolist() if isinstance(A, (MatrixBase, NDimArray)) else [[A]]'
          'MM = [[AA[i][j]] for i, j in zip(rr, cc)]'
-         'M = make_matrix_or_array(MM)'
+         'M = make_2d_sym(MM)'
          'return M,'};
 
   rr = num2cell(int32(r-1));
diff --git a/inst/@sym/private/mat_rclist_asgn.m b/inst/@sym/private/mat_rclist_asgn.m
index 836714519..cc5fb0760 100644
--- a/inst/@sym/private/mat_rclist_asgn.m
+++ b/inst/@sym/private/mat_rclist_asgn.m
@@ -84,7 +84,7 @@
          '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_matrix_or_array(MM)'
+         'M = make_2d_sym(MM)'
          'return M,'};
 
   rr = num2cell(int32(r-1));
diff --git a/inst/@sym/subs.m b/inst/@sym/subs.m
index 132291199..a612cf8c7 100644
--- a/inst/@sym/subs.m
+++ b/inst/@sym/subs.m
@@ -262,7 +262,7 @@
     '        yyy = [y[i, j] if y.is_Matrix else y for y in yy]'
     '        sublist = list(zip(xx, yyy))'
     '        g[i][j] = f.subs(sublist, simultaneous=True).doit()'
-    'return make_matrix_or_array(g)'
+    'return make_2d_sym(g)'
   };
 
   g = pycall_sympy__ (cmd, f, in, out);
diff --git a/inst/@sym/vertcat.m b/inst/@sym/vertcat.m
index 23a66087c..1e49a5116 100644
--- a/inst/@sym/vertcat.m
+++ b/inst/@sym/vertcat.m
@@ -60,7 +60,7 @@
          '    raise ShapeError(msg)'
          'CCC = [as_list_of_list(x) for x in args]'
          'CC = flatten(CCC, levels=1)'
-         'return make_matrix_or_array(CC)'};
+         'return make_2d_sym(CC)'};
 
   args = cellfun (@sym, varargin, 'UniformOutput', false);
   h = pycall_sympy__ (cmd, args{:});
diff --git a/inst/private/python_header.py b/inst/private/python_header.py
index df855dd15..c4013c01b 100644
--- a/inst/private/python_header.py
+++ b/inst/private/python_header.py
@@ -244,27 +244,27 @@ def octoutput(x, et):
 
 
 try:
-    def make_matrix_or_array(it_of_it, dbg_no_array=False):
+    def make_2d_sym(it_of_it, dbg_matrix_only=False):
         # should be kept in sync with the same function
         # defined in inst/private/python_ipc_native.m
-        # FIXME: dbg_no_array is currently used for debugging,
-        # remove it once sympy drops non-Expr supp in Matrix
+        # FIXME: dbg_matrix_only is used for debugging, remove
+        # it once sympy drops non-Expr support in Matrix
         """
-        Given an iterable of iterables of syms IT_OF_IT
-        If all elements of IT_OF_IT are Expr,
-        construct the corresponding Matrix.
-        Otherwise, construct the corresponding 2D array.
+        Given an iterable of iterables of syms IT_OF_IT.
+        If all elements of IT_OF_IT are Expr, construct the
+        corresponding Matrix. Otherwise, construct the
+        corresponding non-Matrix 2D sym.
         """
         ls_of_ls = [[elt for elt in it] for it in it_of_it]
         elts = flatten(ls_of_ls, levels=1)
         if Version(spver) <= Version("1.11.1"):
             # never use Array on older SymPy
-            dbg_no_array = True
-        if (dbg_no_array
+            dbg_matrix_only = True
+        if (dbg_matrix_only
             or all(isinstance(elt, Expr) for elt in elts)):
             return Matrix(ls_of_ls)
         else:
-            dbout(f"make_matrix_or_array: making 2D Array...")
+            dbout(f"make_2d_sym: constructing 2D sym...")
             return Array(ls_of_ls)
 except:
     echo_exception_stdout("in python_header defining fcns block 5")
diff --git a/inst/private/python_ipc_native.m b/inst/private/python_ipc_native.m
index 956537b52..685a48519 100644
--- a/inst/private/python_ipc_native.m
+++ b/inst/private/python_ipc_native.m
@@ -112,27 +112,27 @@
                     '    # should be kept in sync with the same function'
                     '    # defined in inst/private/python_header.py'
                     '    sys.stderr.write("pydebug: " + str(l) + "\n")'
-                    'def make_matrix_or_array(it_of_it, dbg_no_array=False):'
+                    'def make_2d_sym(it_of_it, dbg_matrix_only=False):'
                     '    # should be kept in sync with the same function'
                     '    # defined in inst/private/python_header.py'
-                    '    # FIXME: dbg_no_array is currently used for debugging,'
-                    '    # remove it once sympy drops non-Expr supp in Matrix'
+                    '    # FIXME: dbg_matrix_only is used for debugging, remove'
+                    '    # it once sympy drops non-Expr support in Matrix'
                     '    """'
-                    '    Given an iterable of iterables of syms IT_OF_IT'
-                    '    If all elements of IT_OF_IT are Expr,'
-                    '    construct the corresponding Matrix.'
-                    '    Otherwise, construct the corresponding 2D array.'
+                    '    Given an iterable of iterables of syms IT_OF_IT.'
+                    '    If all elements of IT_OF_IT are Expr, construct the'
+                    '    corresponding Matrix. Otherwise, construct the'
+                    '    corresponding non-Matrix 2D sym.'
                     '    """'
                     '    ls_of_ls = [[elt for elt in it] for it in it_of_it]'
                     '    elts = flatten(ls_of_ls, levels=1)'
                     '    if Version(spver) <= Version("1.11.1"):'
                     '        # never use Array on older SymPy'
-                    '        dbg_no_array = True'
-                    '    if (dbg_no_array'
+                    '        dbg_matrix_only = True'
+                    '    if (dbg_matrix_only'
                     '        or all(isinstance(elt, Expr) for elt in elts)):'
                     '        return Matrix(ls_of_ls)'
                     '    else:'
-                    '        dbout(f"make_matrix_or_array: making 2D Array...")'
+                    '        dbout(f"make_2d_sym: constructing 2D sym...")'
                     '        return Array(ls_of_ls)'
                   }, newl))
     have_headers = true;

From 27b64eb12918cbf240420164974fb3eaabb6d734 Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Tue, 13 Sep 2022 00:23:36 -0700
Subject: [PATCH 33/40] make_2d_sym: support flat input and shape

Fixes #1236.  For now, just add a new underscore version that takes a
shape kwarg.
---
 inst/@sym/private/elementwise_op.m |  6 +++---
 inst/private/python_header.py      | 15 +++++++++++++++
 inst/private/python_ipc_native.m   | 15 +++++++++++++++
 3 files changed, 33 insertions(+), 3 deletions(-)

diff --git a/inst/@sym/private/elementwise_op.m b/inst/@sym/private/elementwise_op.m
index 996bd43d5..74c998442 100644
--- a/inst/@sym/private/elementwise_op.m
+++ b/inst/@sym/private/elementwise_op.m
@@ -96,11 +96,11 @@
           '# dbout(f"at least one matrix param, shape={q.shape}")'
           'assert len(q.shape) == 2, "non-2D arrays/tensors not yet supported"'
           'm, n = q.shape'
-          'g = [[0]*n for i in range(m)]'
+          'g = []'
           'for i in range(m):'
           '    for j in range(n):'
-          '        g[i][j] = _op(*[k[i, j] if isinstance(k, (MatrixBase, NDimArray)) else k for k in _ins])'
-          'return make_2d_sym(g)' ];
+          '        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{:});
 
diff --git a/inst/private/python_header.py b/inst/private/python_header.py
index c4013c01b..2c2651900 100644
--- a/inst/private/python_header.py
+++ b/inst/private/python_header.py
@@ -266,6 +266,21 @@ def make_2d_sym(it_of_it, dbg_matrix_only=False):
         else:
             dbout(f"make_2d_sym: constructing 2D sym...")
             return Array(ls_of_ls)
+    def _make_2d_sym(flat, shape, dbg_matrix_only=False):
+        """
+        If all elements of FLAT are Expr, construct the
+        corresponding Matrix. Otherwise, construct the
+        corresponding non-Matrix 2D sym.
+        """
+        flat = list(flat)
+        if Version(spver) <= Version("1.11.1"):
+            # never use Array on older SymPy
+            dbg_matrix_only = True
+        if (dbg_matrix_only
+            or all(isinstance(elt, Expr) for elt in flat)):
+            return Matrix(*shape, flat)
+        dbout(f"make_2d_sym: constructing 2D sym...")
+        return Array(flat, shape)
 except:
     echo_exception_stdout("in python_header defining fcns block 5")
     raise
diff --git a/inst/private/python_ipc_native.m b/inst/private/python_ipc_native.m
index 685a48519..b62c47e02 100644
--- a/inst/private/python_ipc_native.m
+++ b/inst/private/python_ipc_native.m
@@ -134,6 +134,21 @@
                     '    else:'
                     '        dbout(f"make_2d_sym: constructing 2D sym...")'
                     '        return Array(ls_of_ls)'
+                    'def _make_2d_sym(flat, shape, dbg_matrix_only=False):'
+                    '    """'
+                    '    If all elements of FLAT are Expr, construct the'
+                    '    corresponding Matrix. Otherwise, construct the'
+                    '    corresponding non-Matrix 2D sym.'
+                    '    """'
+                    '    flat = list(flat)'
+                    '    if Version(spver) <= Version("1.11.1"):'
+                    '        # never use Array on older SymPy'
+                    '        dbg_matrix_only = True'
+                    '    if (dbg_matrix_only'
+                    '        or all(isinstance(elt, Expr) for elt in flat)):'
+                    '        return Matrix(*shape, flat)'
+                    '    dbout(f"make_2d_sym: constructing 2D sym...")'
+                    '    return Array(flat, shape)'
                   }, newl))
     have_headers = true;
   end

From 174cbdcfc0ab305720fc8c89c87601511cfe154a Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Tue, 13 Sep 2022 00:51:49 -0700
Subject: [PATCH 34/40] Workaround nsolve with immutable x vars

---
 inst/vpasolve.m | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/inst/vpasolve.m b/inst/vpasolve.m
index 6fb80e780..2b0a7308b 100644
--- a/inst/vpasolve.m
+++ b/inst/vpasolve.m
@@ -105,6 +105,9 @@
       '(e, x, x0, n) = _ins'
       'import mpmath'
       'mpmath.mp.dps = n'
+      'if isinstance(x, MatrixBase):'
+      '    # TODO: find/file upstream issue'
+      '    x = x.as_mutable()'
       'r = nsolve(e, x, x0)'
       'return r' };
   r = pycall_sympy__ (cmd, sym(e), x, x0, n);

From ee92b6619b700b7ad3e324ba59fa8a5f910fe9f6 Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Tue, 13 Sep 2022 09:25:09 -0700
Subject: [PATCH 35/40] subs: after Array port was creating incorrect empties

Fixes #1236.  Or at least the part of #1236 that pertains to `subs`.
Added new tests.
---
 inst/@sym/subs.m | 16 +++++++++++++---
 1 file changed, 13 insertions(+), 3 deletions(-)

diff --git a/inst/@sym/subs.m b/inst/@sym/subs.m
index a612cf8c7..056ba60a8 100644
--- a/inst/@sym/subs.m
+++ b/inst/@sym/subs.m
@@ -256,13 +256,13 @@
     'sizes.discard((1, 1))'
     'assert len(sizes) == 1, "all substitutions must be same size or scalar"'
     'm, n = sizes.pop()'
-    'g = [[0]*n for i in range(m)]'
+    'g = []'
     'for i in range(m):'
     '    for j in range(n):'
     '        yyy = [y[i, j] if y.is_Matrix else y for y in yy]'
     '        sublist = list(zip(xx, yyy))'
-    '        g[i][j] = f.subs(sublist, simultaneous=True).doit()'
-    'return make_2d_sym(g)'
+    '        g.append(f.subs(sublist, simultaneous=True).doit())'
+    'return _make_2d_sym(g, shape=(m, n))'
   };
 
   g = pycall_sympy__ (cmd, f, in, out);
@@ -459,3 +459,13 @@
 %! g = subs (f);
 %! syms x
 %! assert (isequal (g, 6*x*y))
+
+%!test
+%! % issue #1236, correct empty sizes
+%! syms x
+%! g = subs (x, x, zeros (0, 0));
+%! assert (size (g), [0 0])
+%! g = subs (x, x, zeros (0, 3));
+%! assert (size (g), [0 3])
+%! g = subs (x, x, zeros (3, 0));
+%! assert (size (g), [3 0])

From 6f06abb91a9c9a0b3a10e27250e2c690fff29857 Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Tue, 13 Sep 2022 09:49:36 -0700
Subject: [PATCH 36/40] vertcat: tests for empty cases

---
 inst/@sym/vertcat.m | 42 +++++++++++++++++++++++++++++++++++-------
 1 file changed, 35 insertions(+), 7 deletions(-)

diff --git a/inst/@sym/vertcat.m b/inst/@sym/vertcat.m
index 1e49a5116..a93cc4be1 100644
--- a/inst/@sym/vertcat.m
+++ b/inst/@sym/vertcat.m
@@ -1,5 +1,5 @@
 %% 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.
@@ -123,12 +123,6 @@
 %! a = [v; []; []];
 %! assert (isequal (a, v))
 
-%!xtest
-%! % FIXME: is this Octave bug? worth worrying about
-%! syms x
-%! a = [x; [] []];
-%! assert (isequal (a, x))
-
 %!test
 %! % more empty vectors
 %! v = [sym(1) sym(2)];
@@ -140,6 +134,14 @@
 %! q = sym(ones(0, 3));
 %! w = vertcat(v, q);
 
+%!error <ShapeError>
+%! z03 = sym (zeros (0, 3));
+%! z04 = sym (zeros (0, 4));
+%! % Note Issue #1238: unclear error message:
+%! % [z03; z04];
+%! % but this gives the ShapeError
+%! vertcat(z03, z04);
+
 %!test
 %! % Octave 3.6 bug: should pass on 3.8.1 and matlab
 %! a = [sym(1) 2];
@@ -157,3 +159,29 @@
 %! A = sym ([1 2]);
 %! B = simplify (A);
 %! assert (isequal ([B; A], [A; B]))
+
+%!xtest
+%! % issue #1236, correct empty sizes
+%! syms x
+%! z00 = sym (zeros (0, 0));
+%! z30 = sym (zeros (3, 0));
+%! z40 = sym (zeros (4, 0));
+%! z03 = sym (zeros (0, 3));
+%! assert (size ([z00; z00]), [0 0])
+%! assert (size ([z00; z30]), [3 0])
+%! assert (size ([z30; z30]), [6 0])
+%! assert (size ([z30; z40]), [7 0])
+%! assert (size ([z30; z00; z30]), [6 0])
+%! assert (size ([z03; z03]), [0 3])
+%! assert (size ([z03; z03; z03]), [0 3])
+
+%!test
+%! % special case for the 0x0 empty: no error
+%! z00 = sym (zeros (0, 0));
+%! z03 = sym (zeros (0, 3));
+%! [z00; z03];
+
+%!test
+%! syms x
+%! a = [x; sym([]) []];
+%! assert (isequal (a, x))

From edcaec2acda169b8228f5f07bbf1a82323d9cf54 Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Tue, 13 Sep 2022 10:20:35 -0700
Subject: [PATCH 37/40] vertcat: use a flat list and shape for correct empty

Fixes #1236.
---
 inst/@sym/vertcat.m | 12 +++++++-----
 1 file changed, 7 insertions(+), 5 deletions(-)

diff --git a/inst/@sym/vertcat.m b/inst/@sym/vertcat.m
index a93cc4be1..766106f12 100644
--- a/inst/@sym/vertcat.m
+++ b/inst/@sym/vertcat.m
@@ -49,18 +49,20 @@
          '    return isinstance(x, (MatrixBase, NDimArray))'
          'def number_of_columns(x):'
          '    return x.shape[1] if is_matrix_or_array(x) else 1'
+         'def number_of_rows(x):'
+         '    return x.shape[0] if is_matrix_or_array(x) else 1'
          'def all_equal(*ls):'
          '    return True if ls == [] else all(ls[0] == x for x in ls[1:])'
-         'def as_list_of_list(x):'
-         '    return x.tolist() if is_matrix_or_array(x) else [[x]]'
          'args = [x for x in _ins if x != zeros(0, 0)] # remove 0x0 matrices'
          'ncols = [number_of_columns(x) for x in args]'
+         'nrows = [number_of_rows(x) for x in args]'
          'if not all_equal(*ncols):'
          '    msg = "vertcat: all inputs must have the same number of columns"'
          '    raise ShapeError(msg)'
-         'CCC = [as_list_of_list(x) for x in args]'
+         'ncol = 0 if not ncols else ncols[0]'
+         'CCC = [flatten(x, levels=1) if is_matrix_or_array(x) else [x] for x in args]'
          'CC = flatten(CCC, levels=1)'
-         'return make_2d_sym(CC)'};
+         'return _make_2d_sym(CC, shape=(sum(nrows), ncol))' };
 
   args = cellfun (@sym, varargin, 'UniformOutput', false);
   h = pycall_sympy__ (cmd, args{:});
@@ -160,7 +162,7 @@
 %! B = simplify (A);
 %! assert (isequal ([B; A], [A; B]))
 
-%!xtest
+%!test
 %! % issue #1236, correct empty sizes
 %! syms x
 %! z00 = sym (zeros (0, 0));

From 90af895ea4d952c45b5ff51e8b6b24addbcea910 Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Tue, 13 Sep 2022 10:30:47 -0700
Subject: [PATCH 38/40] horzcat/vertcat: add tests for horzcat as well

---
 inst/@sym/horzcat.m | 31 ++++++++++++++++++++++++++++++-
 inst/@sym/vertcat.m | 22 +++++++++++-----------
 2 files changed, 41 insertions(+), 12 deletions(-)

diff --git a/inst/@sym/horzcat.m b/inst/@sym/horzcat.m
index fed82f1cf..295ad7518 100644
--- a/inst/@sym/horzcat.m
+++ b/inst/@sym/horzcat.m
@@ -1,5 +1,5 @@
 %% 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.
@@ -117,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])
diff --git a/inst/@sym/vertcat.m b/inst/@sym/vertcat.m
index 766106f12..2f71a9ffc 100644
--- a/inst/@sym/vertcat.m
+++ b/inst/@sym/vertcat.m
@@ -144,6 +144,17 @@
 %! % but this gives the ShapeError
 %! vertcat(z03, z04);
 
+%!test
+%! % special case for the 0x0 empty: no error
+%! z00 = sym (zeros (0, 0));
+%! z03 = sym (zeros (0, 3));
+%! [z00; z03];
+
+%!test
+%! syms x
+%! a = [x; sym([]) []];
+%! assert (isequal (a, x))
+
 %!test
 %! % Octave 3.6 bug: should pass on 3.8.1 and matlab
 %! a = [sym(1) 2];
@@ -176,14 +187,3 @@
 %! assert (size ([z30; z00; z30]), [6 0])
 %! assert (size ([z03; z03]), [0 3])
 %! assert (size ([z03; z03; z03]), [0 3])
-
-%!test
-%! % special case for the 0x0 empty: no error
-%! z00 = sym (zeros (0, 0));
-%! z03 = sym (zeros (0, 3));
-%! [z00; z03];
-
-%!test
-%! syms x
-%! a = [x; sym([]) []];
-%! assert (isequal (a, x))

From c9252ee9df9e9bfecf9f1a3c0f3fb9e82571eab2 Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Tue, 13 Sep 2022 10:38:58 -0700
Subject: [PATCH 39/40] Port helpers from list-of-lists to flat+shape

Related to #1236.  I don't see a problem with the list-of-lists
approach, but these are the only routine using it so we can delete some
code if we port them to flat + shape.
---
 inst/@sym/private/mat_rclist_access.m | 6 +++---
 inst/@sym/private/mat_rclist_asgn.m   | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/inst/@sym/private/mat_rclist_access.m b/inst/@sym/private/mat_rclist_access.m
index c46036a31..92a7c06c2 100644
--- a/inst/@sym/private/mat_rclist_access.m
+++ b/inst/@sym/private/mat_rclist_access.m
@@ -35,9 +35,9 @@
 
   cmd = {'(A, rr, cc) = _ins'
          'AA = A.tolist() if isinstance(A, (MatrixBase, NDimArray)) else [[A]]'
-         'MM = [[AA[i][j]] for i, j in zip(rr, cc)]'
-         'M = make_2d_sym(MM)'
-         'return M,'};
+         '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));
diff --git a/inst/@sym/private/mat_rclist_asgn.m b/inst/@sym/private/mat_rclist_asgn.m
index cc5fb0760..1641229e3 100644
--- a/inst/@sym/private/mat_rclist_asgn.m
+++ b/inst/@sym/private/mat_rclist_asgn.m
@@ -59,7 +59,7 @@
   % 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, rr, cc, b) = _ins'
          'assert A == [] or not isinstance(A, list), "unexpectedly non-empty list: report bug!"'

From 47834bd12b2304e87b1c468929c0c38bf4b21faf Mon Sep 17 00:00:00 2001
From: "Colin B. Macdonald" <cbm@m.fsf.org>
Date: Tue, 12 Nov 2024 09:15:45 -0800
Subject: [PATCH 40/40] isntance doesn't need special None handling

---
 inst/@sym/isAlways.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/inst/@sym/isAlways.m b/inst/@sym/isAlways.m
index b27d5db6a..94cb79a0c 100644
--- a/inst/@sym/isAlways.m
+++ b/inst/@sym/isAlways.m
@@ -128,7 +128,7 @@
 
   cmd = vertcat(cmd, {
     '(x, map_unknown_to) = _ins'
-    'if x is not None and isinstance(x, (MatrixBase, NDimArray)):'
+    'if isinstance(x, (MatrixBase, NDimArray)):'
     '    r = [a for a in flatten(transpose(x).tolist(), levels=1)]'  % note tranpose
     'else:'
     '    r = [x,]'