Direct Delta Mush Skinning and Variants Eq. 1

v_i = ∑_j w_i,j M_j u_i

where

w ∈ ℝ^(n×m)
M_j ∈ ℝ^(4×4)
u_i ∈ ℝ^4

/*
v_i = ∑_j w_i,j M_j u_i

where

w ∈ ℝ^(n×m)
M_j ∈ ℝ^(4×4)
u_i ∈ ℝ^4
*/
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <iostream>
#include <set>

struct delta_mush_1ResultType {
    std::vector<Eigen::Matrix<double, 4, 1>> v;
    delta_mush_1ResultType(const std::vector<Eigen::Matrix<double, 4, 1>> & v)
    : v(v)
    {}
};

delta_mush_1ResultType delta_mush_1(
    const Eigen::MatrixXd & w,
    const std::vector<Eigen::Matrix<double, 4, 4>> & M,
    const std::vector<Eigen::Matrix<double, 4, 1>> & u)
{
    const long n = w.rows();
    const long m = w.cols();
    const long dim_0 = M.size();
    const long dim_1 = u.size();
    assert( w.rows() == n );
    assert( dim_0 == m );
    assert( dim_1 == n );

    std::vector<Eigen::Matrix<double, 4, 1>> v(dim_1);
    for( int i=1; i<=dim_1; i++){
        Eigen::MatrixXd sum_0 = Eigen::MatrixXd::Zero(4, 1);
        for(int j=1; j<=M.size(); j++){
            sum_0 += w(i-1, j-1) * M.at(j-1) * u.at(i-1);
        }
        v.at(i-1) = sum_0;
    }

    return delta_mush_1ResultType(v);
}


void generateRandomData(Eigen::MatrixXd & w,
    std::vector<Eigen::Matrix<double, 4, 4>> & M,
    std::vector<Eigen::Matrix<double, 4, 1>> & u)
{
    const int n = rand()%10;
    const int dim_1 = n;
    const int m = rand()%10;
    const int dim_0 = m;
    w = Eigen::MatrixXd::Random(n, m);
    M.resize(dim_0);
    for(int i=0; i<dim_0; i++){
        M[i] = Eigen::MatrixXd::Random(4, 4);
    }
    u.resize(dim_1);
    for(int i=0; i<dim_1; i++){
        u[i] = Eigen::VectorXd::Random(4);
    }
}


int main(int argc, char *argv[])
{
    srand((int)time(NULL));
    Eigen::MatrixXd w;
    std::vector<Eigen::Matrix<double, 4, 4>> M;
    std::vector<Eigen::Matrix<double, 4, 1>> u;
    generateRandomData(w, M, u);
    delta_mush_1ResultType func_value = delta_mush_1(w, M, u);
    std::cout<<"vector return value:"<<std::endl;
    for(int i=0; i<func_value.v.size(); i++){
        std::cout<<"i:"<<i<<", value:\n"<<func_value.v.at(i)<<std::endl;
    }
    return 0;
}

"""
v_i = ∑_j w_i,j M_j u_i

where

w ∈ ℝ^(n×m)
M_j ∈ ℝ^(4×4)
u_i ∈ ℝ^4
"""
import numpy as np
import scipy
import scipy.linalg
from scipy import sparse
from scipy.integrate import quad
from scipy.optimize import minimize


class delta_mush_1ResultType:
    def __init__( self, v):
        self.v = v


def delta_mush_1(w, M, u):
    w = np.asarray(w, dtype=np.float64)
    M = np.asarray(M, dtype=np.float64)
    u = np.asarray(u, dtype=np.float64)

    n = w.shape[0]
    m = w.shape[1]
    dim_0 = M.shape[0]
    dim_1 = u.shape[0]
    assert w.shape == (n, m)
    assert M.shape == (dim_0, 4, 4)
    assert u.shape == (dim_1, 4, )
    assert dim_0 == m 
    assert dim_1 == n 

    v = np.zeros((dim_1, 4, ))
    for i in range(1, dim_1+1):
        sum_0 = np.zeros((4, ))
        for j in range(1, len(M)+1):
            sum_0 += w[i-1, j-1] * M[j-1] @ u[i-1]
        v[i-1] = sum_0
    return delta_mush_1ResultType(v)


def generateRandomData():
    n = np.random.randint(10)
    dim_1 = n
    m = np.random.randint(10)
    dim_0 = m
    w = np.random.randn(n, m)
    M = np.random.randn(dim_0, 4, 4)
    u = np.random.randn(dim_1, 4, )
    return w, M, u


if __name__ == '__main__':
    w, M, u = generateRandomData()
    print("w:", w)
    print("M:", M)
    print("u:", u)
    func_value = delta_mush_1(w, M, u)
    print("return value: ", func_value.v)

function output = delta_mush_1(w, M, u)
% output = delta_mush_1(w, M, u)
%
%    v_i = ∑_j w_i,j M_j u_i
%    
%    where
%    
%    w ∈ ℝ^(n×m)
%    M_j ∈ ℝ^(4×4)
%    u_i ∈ ℝ^4
    if nargin==0
        warning('generating random input data');
        [w, M, u] = generateRandomData();
    end
    function [w, M, u] = generateRandomData()
        n = randi(10);
        dim_1 = n;
        m = randi(10);
        dim_0 = m;
        w = randn(n, m);
        M = randn(dim_0,4,4);
        u = randn(dim_1,4);
    end

    n = size(w, 1);
    m = size(w, 2);
    dim_0 = size(M, 1);
    dim_1 = size(u, 1);
    assert( isequal(size(w), [n, m]) );
    assert( isequal(size(M), [dim_0, 4, 4]) );
    assert( isequal(size(u), [dim_1, 4]) );
    assert( dim_0 == m );
    assert( dim_1 == n );

    v = zeros(dim_1, 4);
    for i = 1:dim_1
        sum_0 = zeros(4,1);
        for j = 1:size(M, 1)
            sum_0 = sum_0 + w(i, j) * squeeze(M(j,:,:)) * u(i,:)';
        end
        v(i,:) = (sum_0)';
    end
    output.v = v;
end

\documentclass[12pt]{article}
\usepackage{mathdots}
\usepackage[bb=boondox]{mathalfa}
\usepackage{mathtools}
\usepackage{amssymb}
\usepackage{libertine}
\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator*{\argmin}{arg\,min}
\usepackage[paperheight=8in,paperwidth=4in,margin=.3in,heightrounded]{geometry}
\let\originalleft\left
\let\originalright\right
\renewcommand{\left}{\mathopen{}\mathclose\bgroup\originalleft}
\renewcommand{\right}{\aftergroup\egroup\originalright}
\begin{document}

\begin{center}
\resizebox{\textwidth}{!} 
{
\begin{minipage}[c]{\textwidth}
\begin{align*}
\mathit{v}_{ \mathit{i} } & = \sum_\mathit{j} \mathit{w}_{\mathit{i}, \mathit{j}}\mathit{M}_{ \mathit{j} }\mathit{u}_{ \mathit{i} } \\
\intertext{where} 
\mathit{w} & \in \mathbb{R}^{ \mathit{n} \times \mathit{m} } \\
\mathit{M}_{\mathit{j}} & \in \mathbb{R}^{ 4 \times 4 } \\
\mathit{u}_{\mathit{i}} & \in \mathbb{R}^{ 4} \\
\\
\end{align*}
\end{minipage}
}
\end{center}

\end{document}