From a7c44c0d78bd6be66b6795ed339751657b1c3a3d 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/33] 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 925e7a50b..05d016c05 100644
--- a/inst/private/python_header.py
+++ b/inst/private/python_header.py
@@ -152,7 +152,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
@@ -161,6 +161,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 b618e5fdd0bb1080d913d965ed4040e6f1f6bb6a 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/33] 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 925e7a50b..859bc650a 100644
--- a/inst/private/python_header.py
+++ b/inst/private/python_header.py
@@ -235,6 +235,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 d1bd19328..698c93a81 100644
--- a/inst/private/python_ipc_native.m
+++ b/inst/private/python_ipc_native.m
@@ -111,6 +111,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 53059ac82043d4a472a997aab07072a1b34e91f3 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/33] @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 54c827199578e49cc24eae99a6202d6c6861d325 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/33] @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 ce779776e..f2f999a0b 100644
--- a/inst/@sym/private/mat_rclist_asgn.m
+++ b/inst/@sym/private/mat_rclist_asgn.m
@@ -1,6 +1,7 @@
 %% Copyright (C) 2014, 2016-2017, 2019, 2022 Colin B. Macdonald
 %% Copyright (C) 2020 Mike Miller
 %% Copyright (C) 2020 Fernando Alvarruiz
+%% Copyright (C) 2022 Alex Vong
 %%
 %% This file is part of OctSymPy.
 %%
@@ -59,40 +60,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 f1b210fdf07f92ce03edee5fa6f2564d343b03be 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/33] @sym/{horzcat,vertcat}: Fix typo.

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

diff --git a/inst/@sym/horzcat.m b/inst/@sym/horzcat.m
index cac1fbd65..0f562c2ee 100644
--- a/inst/@sym/horzcat.m
+++ b/inst/@sym/horzcat.m
@@ -1,4 +1,5 @@
 %% Copyright (C) 2014-2017, 2019 Colin B. Macdonald
+%% Copyright (C) 2022 Alex Vong
 %%
 %% This file is part of OctSymPy.
 %%
@@ -21,7 +22,7 @@
 %% @defop  Method   @@sym {horzcat} {(@var{x}, @var{y}, @dots{})}
 %% @defopx Operator @@sym {[@var{x}, @var{y}, @dots{}]} {}
 %% @defopx Operator @@sym {[@var{x} @var{y} @dots{}]} {}
-%% Horizontally concatentate symbolic arrays.
+%% Horizontally concatenate symbolic arrays.
 %%
 %% Example:
 %% @example
diff --git a/inst/@sym/vertcat.m b/inst/@sym/vertcat.m
index 09a8ea787..713f283bc 100644
--- a/inst/@sym/vertcat.m
+++ b/inst/@sym/vertcat.m
@@ -1,4 +1,5 @@
 %% Copyright (C) 2014-2017, 2019 Colin B. Macdonald
+%% Copyright (C) 2022 Alex Vong
 %%
 %% This file is part of OctSymPy.
 %%
@@ -20,7 +21,7 @@
 %% @documentencoding UTF-8
 %% @defop  Method   @@sym {vertcat} {(@var{x}, @var{y}, @dots{})}
 %% @defopx Operator @@sym {[@var{x}; @var{y}; @dots{}]} {}
-%% Vertically concatentate symbolic arrays.
+%% Vertically concatenate symbolic arrays.
 %%
 %% Example:
 %% @example

From d3b154731803c95460af08722e7a62a718d20870 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/33] @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 713f283bc..23dacb960 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
   % compatibilty
-  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 136184d6df506bc9945d8f341e9b46691895c5b4 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/33] @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 0d9360275c10a2b6f2cd57e4c56c84b64ddc8d13 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/33] @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 0f562c2ee..bbb583913 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
-  % compatibilty
-  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 fb044e387ddbe5c25d419928de0100d332b3388c 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/33] @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 f2f999a0b..852ff38dc 100644
--- a/inst/@sym/private/mat_rclist_asgn.m
+++ b/inst/@sym/private/mat_rclist_asgn.m
@@ -60,8 +60,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 65da8ab60c3a7b46a285697e20a1425ced570e6b 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/33] @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 56bd71902ca0faf5fa88c1bc7fbeba8bf6f99ce7 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/33] 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 c6c9ee3b4..a76cea0da 100644
--- a/inst/private/python_header.py
+++ b/inst/private/python_header.py
@@ -241,7 +241,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 30af4ef043df5fb295be63d58caa921e6d8cc8ef 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 12/33] 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 8b10a44f37569d3c79395b4e6d705a6252c9f36a 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 13/33] 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 0e0ea4bdb..cdefda7ee 100644
--- a/inst/@sym/logical.m
+++ b/inst/@sym/logical.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.
 %%
@@ -100,8 +100,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 53fbee24bb605f3915c07a1e0e2ba1566afa65fb 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 14/33] 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 cdefda7ee..b105e3273 100644
--- a/inst/@sym/logical.m
+++ b/inst/@sym/logical.m
@@ -178,7 +178,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 d36285cd0d964aed088583832d510f27603a2714 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 15/33] 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 b105e3273..eed490c6c 100644
--- a/inst/@sym/logical.m
+++ b/inst/@sym/logical.m
@@ -100,7 +100,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 1381757e56ba0afed82a2cd86bf8709b7f05a3b5 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 16/33] 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 bda13094f..ec69f3a1a 100644
--- a/inst/@sym/isAlways.m
+++ b/inst/@sym/isAlways.m
@@ -131,8 +131,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())]'  % note tranpose
     'else:'
     '    r = [x,]'
     'r = [simplify_tfn(a) for a in r]'

From 37724d624c6952669709c85e058f3e0fc72b604c 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 17/33] 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 698c93a81..6805e47ae 100644
--- a/inst/private/python_ipc_native.m
+++ b/inst/private/python_ipc_native.m
@@ -111,7 +111,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 a1d3608adc06af89bd3ba69fc555047a118ee087 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/33] 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 ec69f3a1a..6c253027e 100644
--- a/inst/@sym/isAlways.m
+++ b/inst/@sym/isAlways.m
@@ -132,7 +132,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 = [simplify_tfn(a) for a in r]'
diff --git a/inst/@sym/logical.m b/inst/@sym/logical.m
index eed490c6c..2682f548d 100644
--- a/inst/@sym/logical.m
+++ b/inst/@sym/logical.m
@@ -101,7 +101,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 79d0526044bec4984a13074452c0cd026a5963c0 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/33] 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 543344090..c54cfeca5 100644
--- a/inst/@sym/subs.m
+++ b/inst/@sym/subs.m
@@ -256,14 +256,14 @@
     'sizes = {(a.shape if a.is_Matrix else (1, 1)) for a in yy}'
     'sizes.discard((1, 1))'
     'assert len(sizes) == 1, "all substitions 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 ef9744d1d63f0a517b121ea7e8c37b9a745fb55c 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/33] 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 852ff38dc..1932efa87 100644
--- a/inst/@sym/private/mat_rclist_asgn.m
+++ b/inst/@sym/private/mat_rclist_asgn.m
@@ -80,7 +80,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 551081bb9ebe463eec0ea18f5367b0e5a2701524 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/33] 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 b24005e2c4a5089244c0703e65c92a5f10cd46fd 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 22/33] 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 693e055d0eaa7230310831a3d25e01c537fc8754 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 23/33] 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 9ce0a2657128be4e3f946645e8641fa5d8d3fc70 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 24/33] 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 bd795e8a6bd2547e24e6c6df62e4f2952a0fb3b6 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 25/33] @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 23dacb960..ea2b5c32c 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 1f8163e0e2c948fd55641facac83064c9e69dfba 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 26/33] @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 852ff38dc..0d818c1f1 100644
--- a/inst/@sym/private/mat_rclist_asgn.m
+++ b/inst/@sym/private/mat_rclist_asgn.m
@@ -72,8 +72,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 c9f6f0f9b1d0edaee5c29ce7d9e4ae2ece8708d4 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 27/33] @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 321ee1709f011ea568a243d5f8722fa6e3a14db9 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 28/33] @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 137aef88e..8e3a19d35 100644
--- a/inst/@sym/private/mat_replace.m
+++ b/inst/@sym/private/mat_replace.m
@@ -2,6 +2,7 @@
 %% Copyright (C) 2016 Lagu
 %% Copyright (C) 2016 Abhinav Tripathi
 %% Copyright (C) 2020 Fernando Alvarruiz
+%% Copyright (C) 2022 Alex Vong
 %%
 %% This file is part of OctSymPy.
 %%
@@ -154,7 +155,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,'
@@ -171,7 +173,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 59176a2d17418c5288ae07f5052c2abba9eb1a1e 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 29/33] 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 a76cea0da..a0d5e1cc0 100644
--- a/inst/private/python_header.py
+++ b/inst/private/python_header.py
@@ -241,18 +241,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 6805e47ae..58ce5d277 100644
--- a/inst/private/python_ipc_native.m
+++ b/inst/private/python_ipc_native.m
@@ -111,18 +111,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 5ea9b3039eb16701e3cc63bc0eb9be41138ecd70 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/33] 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 3c0e95c1118946ef7be3350ecbd609552970052c 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/33] 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 a0d5e1cc0..1578c40de 100644
--- a/inst/private/python_header.py
+++ b/inst/private/python_header.py
@@ -254,6 +254,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 58ce5d277..168d7027b 100644
--- a/inst/private/python_ipc_native.m
+++ b/inst/private/python_ipc_native.m
@@ -124,6 +124,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 e74ce50c3a52c3756b8b5d90885b0bb61a30f70f 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/33] 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 5bf107687..bf88254b7 100644
--- a/inst/@sym/private/mat_rclist_asgn.m
+++ b/inst/@sym/private/mat_rclist_asgn.m
@@ -83,7 +83,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 c54cfeca5..aa6c9e093 100644
--- a/inst/@sym/subs.m
+++ b/inst/@sym/subs.m
@@ -263,7 +263,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 ea2b5c32c..3a57250a2 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 1578c40de..c93c83097 100644
--- a/inst/private/python_header.py
+++ b/inst/private/python_header.py
@@ -241,27 +241,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 168d7027b..5bfa94e5a 100644
--- a/inst/private/python_ipc_native.m
+++ b/inst/private/python_ipc_native.m
@@ -111,27 +111,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 d4f54273696b960232ad9f805f57818bf23c8136 Mon Sep 17 00:00:00 2001
From: Alex Vong <alexvong1995 AT protonmail DOT com>
Date: Wed, 31 Aug 2022 18:28:15 +0000
Subject: [PATCH 33/33] Python header: Add more 2D sym functions.

Add more 2D sym functions to abstract over the idea of 2D sym 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.

WIP: This is incomplete. More functions need to use these 2D sym
functions. Not sure if this is needed anymore as TableForm was shown
to be lacking various features such as indexing and taking
transpose. No other candidates are being considered at the moment.

* inst/private/{python_header.py,python_ipc_native.m}: Add new 2D sym
functions 'is_2d_sym', 'is_matrix', 'is_non_matrix_2d_sym',
'list_from_2d_sym' and 'shape_of_2d_sym'.
* inst/@sym/private/mat_rclist_{access,asgn}.m: Use them.
* inst/@sym/{transpose.m,vertcat.m}: Use them.
---
 inst/@sym/private/mat_rclist_access.m |  2 +-
 inst/@sym/private/mat_rclist_asgn.m   |  8 ++++----
 inst/@sym/transpose.m                 |  6 ++----
 inst/@sym/vertcat.m                   |  8 +++-----
 inst/private/python_header.py         | 24 ++++++++++++++++++++++--
 inst/private/python_ipc_native.m      | 24 ++++++++++++++++++++++--
 6 files changed, 54 insertions(+), 18 deletions(-)

diff --git a/inst/@sym/private/mat_rclist_access.m b/inst/@sym/private/mat_rclist_access.m
index c46036a31..8662ac578 100644
--- a/inst/@sym/private/mat_rclist_access.m
+++ b/inst/@sym/private/mat_rclist_access.m
@@ -34,7 +34,7 @@
   end
 
   cmd = {'(A, rr, cc) = _ins'
-         'AA = A.tolist() if isinstance(A, (MatrixBase, NDimArray)) else [[A]]'
+         'AA = list_from_2d_sym(A) if is_2d_sym(A) else [[A]]'
          'MM = [[AA[i][j]] for i, j in zip(rr, cc)]'
          'M = make_2d_sym(MM)'
          'return M,'};
diff --git a/inst/@sym/private/mat_rclist_asgn.m b/inst/@sym/private/mat_rclist_asgn.m
index bf88254b7..91b60762b 100644
--- a/inst/@sym/private/mat_rclist_asgn.m
+++ b/inst/@sym/private/mat_rclist_asgn.m
@@ -65,13 +65,13 @@
          'if A == []:'
          '    AA = []'
          '    (nrows_A, ncols_A) = (0, 0)'
-         'elif isinstance(A, (MatrixBase, NDimArray)):'
-         '    AA = A.tolist()'
-         '    (nrows_A, ncols_A) = A.shape'
+         'elif is_2d_sym(A):'
+         '    AA = list_from_2d_sym(A)'
+         '    (nrows_A, ncols_A) = shape_of_2d_sym(A)'
          'else:'
          '    AA = [[A]]'
          '    (nrows_A, ncols_A) = (1, 1)'
-         'bb = b.tolist() if isinstance(b, (MatrixBase, NDimArray)) else [[b]]'
+         'bb = list_from_2d_sym(b) if is_2d_sym(b) else [[b]]'
          'entries = dict(zip(zip(rr, cc), flatten(bb, levels=1)))'
          'def entry(i, j):'
          '    if (i, j) in entries:'
diff --git a/inst/@sym/transpose.m b/inst/@sym/transpose.m
index 161554803..316c14035 100644
--- a/inst/@sym/transpose.m
+++ b/inst/@sym/transpose.m
@@ -66,10 +66,8 @@
     print_usage ();
   end
 
-  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'};
+  cmd = {'x, = _ins'
+         'return transpose(x) if is_2d_sym(x) else x'};
 
   z = pycall_sympy__ (cmd, x);
 
diff --git a/inst/@sym/vertcat.m b/inst/@sym/vertcat.m
index 3a57250a2..a9c5d3ce4 100644
--- a/inst/@sym/vertcat.m
+++ b/inst/@sym/vertcat.m
@@ -45,14 +45,12 @@
 
   % special case for 0x0 but other empties should be checked for
   % compatibilty
-  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'
+  cmd = {'def number_of_columns(x):'
+         '    return shape_of_2d_sym(x)[1] if is_2d_sym(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]]'
+         '    return list_from_2d_sym(x) if is_2d_sym(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):'
diff --git a/inst/private/python_header.py b/inst/private/python_header.py
index c93c83097..4c3ae2fd6 100644
--- a/inst/private/python_header.py
+++ b/inst/private/python_header.py
@@ -241,9 +241,10 @@ def octoutput(x, et):
 
 
 try:
+    # begin: 2D sym funcs
+    # 2D sym funcs defined in inst/private/python_ipc_native.m
+    # and inst/private/python_header.py should be kept in sync
     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_matrix_only is used for debugging, remove
         # it once sympy drops non-Expr support in Matrix
         """
@@ -262,7 +263,26 @@ def make_2d_sym(it_of_it, dbg_matrix_only=False):
             return Matrix(ls_of_ls)
         else:
             dbout(f"make_2d_sym: constructing 2D sym...")
+            # FIXME: should we use Array or TableForm?
             return Array(ls_of_ls)
+    def is_2d_sym(x):
+        types = (MatrixBase, NDimArray, TableForm)
+        return isinstance(x, types)
+    def is_matrix(x):
+        return isinstance(x, MatrixBase)
+    def is_non_matrix_2d_sym(x):
+        return isinstance(x, (NDimArray, TableForm))
+    def list_from_2d_sym(X):
+        if isinstance(X, TableForm):
+            return [[x for x in tup] for tup in X._lines]
+        else:
+            return X.tolist()
+    def shape_of_2d_sym(X):
+        if isinstance(X, TableForm):
+            return (X._h, X._w)
+        else:
+            return X.shape
+    # end: 2D sym funcs
 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 5bfa94e5a..652f17e81 100644
--- a/inst/private/python_ipc_native.m
+++ b/inst/private/python_ipc_native.m
@@ -111,9 +111,10 @@
                     '    # should be kept in sync with the same function'
                     '    # defined in inst/private/python_header.py'
                     '    sys.stderr.write("pydebug: " + str(l) + "\n")'
+                    '# begin: 2D sym funcs'
+                    '# 2D sym funcs defined in inst/private/python_ipc_native.m'
+                    '# and inst/private/python_header.py should be kept in sync'
                     '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_matrix_only is used for debugging, remove'
                     '    # it once sympy drops non-Expr support in Matrix'
                     '    """'
@@ -132,7 +133,26 @@
                     '        return Matrix(ls_of_ls)'
                     '    else:'
                     '        dbout(f"make_2d_sym: constructing 2D sym...")'
+                    '        # FIXME: should we use Array or TableForm?'
                     '        return Array(ls_of_ls)'
+                    'def is_2d_sym(x):'
+                    '    types = (MatrixBase, NDimArray, TableForm)'
+                    '    return isinstance(x, types)'
+                    'def is_matrix(x):'
+                    '    return isinstance(x, MatrixBase)'
+                    'def is_non_matrix_2d_sym(x):'
+                    '    return isinstance(x, (NDimArray, TableForm))'
+                    'def list_from_2d_sym(X):'
+                    '    if isinstance(X, TableForm):'
+                    '        return [[x for x in tup] for tup in X._lines]'
+                    '    else:'
+                    '        return X.tolist()'
+                    'def shape_of_2d_sym(X):'
+                    '    if isinstance(X, TableForm):'
+                    '        return (X._h, X._w)'
+                    '    else:'
+                    '        return X.shape'
+                    '# end: 2D sym funcs'
                   }, newl))
     have_headers = true;
   end