From bd0e7e725725f8375915113a77a726149a9ec6b0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek.necada@aalto.fi>
Date: Tue, 9 Aug 2016 15:40:36 +0300
Subject: [PATCH] ND matrix tensor multiplication

Former-commit-id: b24a08cb61ccac2464bfe9d904fa973cab7596b1
---
 qpms/qpms_p.py | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/qpms/qpms_p.py b/qpms/qpms_p.py
index 2205d9c..c0690d9 100644
--- a/qpms/qpms_p.py
+++ b/qpms/qpms_p.py
@@ -1050,6 +1050,16 @@ def apply_matrix_left(matrix, tensor, axis):
     tmp = np.tensordot(matrix, tensor, axes=(-1,axis))
     return np.moveaxis(tmp, 0, axis)
 
+def apply_ndmatrix_left(matrix,tensor,axes):
+    """
+    Generalized apply_matrix_left, the matrix can have more (2N) abstract dimensions,
+    like M[i,j,k,...z,i,j,k,...,z]. N axes have to be specified in a tuple, corresponding
+    to the axes 0,1,...N-1 of the matrix
+    """
+    N = len(axes)
+    matrix = np.tensordot(matrix, tensor, axes=([-N+axn for axn in range(N)],axes))
+    matrix = np.moveaxis(matrix, range(N), axes)
+    return matrix
 
 ####################
 # Array simulations