Geometry Processing: Parameterization

L_i,j = { w_i,j if (i,j) ∈ E
          0 otherwise
L_i,i = -∑_(ℓ for ℓ ≠ i) L_i,ℓ

where
L ∈ ℝ^(n×n)
w ∈ ℝ^(n×n): edge weight matrix
E ∈ {ℤ²} index: edges

/*
L_i,j = { w_i,j if (i,j) ∈ E
          0 otherwise
L_i,i = -∑_(ℓ for ℓ ≠ i) L_i,ℓ

where
L ∈ ℝ^(n×n)
w ∈ ℝ^(n×n): edge weight matrix
E ∈ {ℤ²} index: edges
*/
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <iostream>
#include <set>

struct course_parameterizationResultType {
    Eigen::SparseMatrix<double> L;
    course_parameterizationResultType(const Eigen::SparseMatrix<double> & L)
    : L(L)
    {}
};

/**
 * course_parameterization
 *
 * @param w  edge weight matrix
 * @param E  edges
 * @return L
 */
course_parameterizationResultType course_parameterization(
    const Eigen::MatrixXd & w,
    const std::set<std::tuple< int, int > > & E)
{
    const long n = w.cols();
    assert( w.rows() == n );

    Eigen::SparseMatrix<double> L(n, n);
    std::vector<Eigen::Triplet<double> > tripletList_L;
    for( int i=1; i<=n; i++){
        for( int j=1; j<=n; j++){
            if(E.find(std::tuple< int, int >(i-1, j-1)) != E.end()){
                tripletList_L.push_back(Eigen::Triplet<double>(i-1, j-1, w(i-1, j-1)));
            }
        }
    }
    L.setFromTriplets(tripletList_L.begin(), tripletList_L.end());

    for( int i=1; i<=n; i++){
        double sum_0 = 0;
        for(int ℓ=1; ℓ<=L.cols(); ℓ++){
            if(ℓ != i){
                sum_0 += L.coeff(i-1, ℓ-1);
            }
        }
        tripletList_L.push_back(Eigen::Triplet<double>(i-1, i-1, -sum_0));
    }
    L.setFromTriplets(tripletList_L.begin(), tripletList_L.end());

    return course_parameterizationResultType(L);
}


void generateRandomData(Eigen::MatrixXd & w,
    std::set<std::tuple< int, int > > & E)
{
    const int n = rand()%10;
    w = Eigen::MatrixXd::Random(n, n);
    const int dim_1 = rand()%10;
    for(int i=0; i<dim_1; i++){
        E.insert(std::make_tuple(rand()%10, rand()%10));
    }
}


int main(int argc, char *argv[])
{
    srand((int)time(NULL));
    Eigen::MatrixXd w;
    std::set<std::tuple< int, int > > E;
    generateRandomData(w, E);
    course_parameterizationResultType func_value = course_parameterization(w, E);
    std::cout<<"return value:\n"<<func_value.L<<std::endl;
    return 0;
}

"""
L_i,j = { w_i,j if (i,j) ∈ E
          0 otherwise
L_i,i = -∑_(ℓ for ℓ ≠ i) L_i,ℓ

where
L ∈ ℝ^(n×n)
w ∈ ℝ^(n×n): edge weight matrix
E ∈ {ℤ²} index: edges
"""
import numpy as np
import scipy
import scipy.linalg
from scipy import sparse
from scipy.integrate import quad
from scipy.optimize import minimize


class course_parameterizationResultType:
    def __init__( self, L):
        self.L = L


def course_parameterization(w, E):
    """
    :param :w : edge weight matrix
    :param :E : edges
    """
    w = np.asarray(w, dtype=np.float64)
    E = frozenset(E)

    n = w.shape[1]
    assert w.shape == (n, n)

    Lij_0 = []
    Lvals_0 = []
    for i in range(1, n+1):
        for j in range(1, n+1):
            if (i-1, j-1) in E:
                Lij_0.append((i-1, j-1))
                Lvals_0.append(w[i-1, j-1])
    sparse_0 = scipy.sparse.coo_matrix((Lvals_0, np.asarray(Lij_0).T), shape=(n, n))
    L = sparse_0
    for i in range(1, n+1):
        sum_0 = 0
        for ℓ in range(1, L.shape[1]+1):
            if(ℓ != i):
                sum_0 += L.tocsr()[i-1, ℓ-1]
        Lij_0.append((i - 1, i - 1))
        Lvals_0.append(-sum_0)
    L = scipy.sparse.coo_matrix((Lvals_0, np.asarray(Lij_0).T), shape=(n, n))

    return course_parameterizationResultType(L)


def generateRandomData():
    n = np.random.randint(10)
    w = np.random.randn(n, n)
    E = []
    dim_0 = np.random.randint(1, 10)
    for i in range(dim_0):
        E.append((np.random.randint(10), np.random.randint(10)))
    return w, E


if __name__ == '__main__':
    w, E = generateRandomData()
    print("w:", w)
    print("E:", E)
    func_value = course_parameterization(w, E)
    print("return value: ", func_value.L)

function output = course_parameterization(w, E)
% output = course_parameterization(w, E)
%
%    L_i,j = { w_i,j if (i,j) ∈ E
%              0 otherwise
%    L_i,i = -∑_(ℓ for ℓ ≠ i) L_i,ℓ
%    
%    where
%    L ∈ ℝ^(n×n)
%    w ∈ ℝ^(n×n): edge weight matrix
%    E ∈ {ℤ²} index: edges
    if nargin==0
        warning('generating random input data');
        [w, E] = generateRandomData();
    end
    function [w, E] = generateRandomData()
        n = randi(10);
        w = randn(n, n);
        E = [];
        dim_2 = randi(10);
        for i = 1:dim_2 
            E = [E;randi(10), randi(10)];
        end
    end

    n = size(w, 2);
    assert( isequal(size(w), [n, n]) );
    assert(size(E,2) == 2)

    Lij_0 = zeros(2,0);
    Lvals_0 = zeros(1,0);
    for i = 1:n
        for j = 1:n
            if ismember([i, j],E,'rows')
                Lij_0(1:2,end+1) = [i;j];
                Lvals_0(end+1) = w(i, j);
            end
        end
    end
    sparse_0 = sparse(Lij_0(1,:),Lij_0(2,:),Lvals_0,n,n);
    L = sparse_0;
    for i = 1:n
        sum_0 = 0;
        for ell = 1:size(L,2)
            if ell ~= i
              sum_0 = sum_0 + L(i, ell);
            end
        end
        Lij_0(1:2,end+1) = [i;i];
        Lvals_0(end+1) = -sum_0;
    end
    L = sparse(Lij_0(1,:),Lij_0(2,:),Lvals_0,n,n);
    output.L = L;
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{L}_{\mathit{i}, \mathit{j}} & = \begin{cases} \mathit{w}_{\mathit{i}, \mathit{j}} & \text{if }  \left( \mathit{i}, \mathit{j} \right) \in \mathit{E}  \\ 0 & \text{otherwise} \end{cases} \\
\mathit{L}_{\mathit{i}, \mathit{i}} & = -\sum_{\mathit{ℓ} ≠ \mathit{i}} \mathit{L}_{\mathit{i}, \mathit{ℓ}} \\
\intertext{where} 
\mathit{L} & \in \mathbb{R}^{ \mathit{n} \times \mathit{n} } \\
\mathit{w} & \in \mathbb{R}^{ \mathit{n} \times \mathit{n} } \text{ edge weight matrix} \\
\mathit{E} & \in \{\mathbb{Z}^{ 2 }\} \mathit{{ index}} \text{ edges} \\
\\
\end{align*}
\end{minipage}
}
\end{center}

\end{document}