Polygon Mesh Processing page 42

given
α_T : ℝ
n_T : ℝ³

`n(v)` = ( ∑_(T for T ∈ `N₁`_v) α_T n_T )/‖ ∑_(T for T ∈ `N₁`_v) α_T n_T ‖

where
v ∈ ℤ
`N₁`_i ∈ {ℤ}

/*
given
α_T : ℝ
n_T : ℝ³

`n(v)` = ( ∑_(T for T ∈ `N₁`_v) α_T n_T )/‖ ∑_(T for T ∈ `N₁`_v) α_T n_T ‖

where
v ∈ ℤ
`N₁`_i ∈ {ℤ}
*/
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <iostream>
#include <set>

struct pmp_42ResultType {
    Eigen::Matrix<double, 3, 1> n_left_parenthesis_v_right_parenthesis;
    pmp_42ResultType(const Eigen::Matrix<double, 3, 1> & n_left_parenthesis_v_right_parenthesis)
    : n_left_parenthesis_v_right_parenthesis(n_left_parenthesis_v_right_parenthesis)
    {}
};

/**
 * pmp_42
 *
 * @param α  ℝ
 * @param n  ℝ³
 * @return n_left_parenthesis_v_right_parenthesis
 */
pmp_42ResultType pmp_42(
    const std::vector<double> & α,
    const std::vector<Eigen::Matrix<double, 3, 1>> & n,
    const int & v,
    const std::vector<std::set<int >> & N₁)
{
    const long dim_0 = α.size();
    const long dim_1 = N₁.size();
    assert( n.size() == dim_0 );

    Eigen::MatrixXd sum_0 = Eigen::MatrixXd::Zero(3, 1);
    for(int T=1; T<=n.size(); T++){
        std::set<int > set_0 = N₁.at(v-1);
        if(set_0.find(int(T)) != set_0.end()){
            sum_0 += α.at(T-1) * n.at(T-1);
        }
    }
    Eigen::MatrixXd sum_1 = Eigen::MatrixXd::Zero(3, 1);
    for(int T=1; T<=n.size(); T++){
        std::set<int > set_1 = N₁.at(v-1);
        if(set_1.find(int(T)) != set_1.end()){
            sum_1 += α.at(T-1) * n.at(T-1);
        }
    }
    Eigen::Matrix<double, 3, 1> n_left_parenthesis_v_right_parenthesis = (sum_0) / double((sum_1).lpNorm<2>());

    return pmp_42ResultType(n_left_parenthesis_v_right_parenthesis);
}


void generateRandomData(std::vector<double> & α,
    std::vector<Eigen::Matrix<double, 3, 1>> & n,
    int & v,
    std::vector<std::set<int >> & N₁)
{
    v = rand() % 10;
    const int dim_0 = rand()%10;
    const int dim_1 = rand()%10;
    α.resize(dim_0);
    for(int i=0; i<dim_0; i++){
        α[i] = rand() % 10;
    }
    n.resize(dim_0);
    for(int i=0; i<dim_0; i++){
        n[i] = Eigen::VectorXd::Random(3);
    }
    N₁.resize(dim_1);
    for(int i=0; i<dim_1; i++){
        const int dim_3 = rand()%10;
        for(int j=0; j<dim_3; j++){
            N₁[i].insert(rand()%10);
        }
    }
}


int main(int argc, char *argv[])
{
    srand((int)time(NULL));
    std::vector<double> α;
    std::vector<Eigen::Matrix<double, 3, 1>> n;
    int v;
    std::vector<std::set<int >> N₁;
    generateRandomData(α, n, v, N₁);
    pmp_42ResultType func_value = pmp_42(α, n, v, N₁);
    std::cout<<"return value:\n"<<func_value.n_left_parenthesis_v_right_parenthesis<<std::endl;
    return 0;
}

"""
given
α_T : ℝ
n_T : ℝ³

`n(v)` = ( ∑_(T for T ∈ `N₁`_v) α_T n_T )/‖ ∑_(T for T ∈ `N₁`_v) α_T n_T ‖

where
v ∈ ℤ
`N₁`_i ∈ {ℤ}
"""
import numpy as np
import scipy
import scipy.linalg
from scipy import sparse
from scipy.integrate import quad
from scipy.optimize import minimize


class pmp_42ResultType:
    def __init__( self, n_left_parenthesis_v_right_parenthesis):
        self.n_left_parenthesis_v_right_parenthesis = n_left_parenthesis_v_right_parenthesis


def pmp_42(α, n, v, N1):
    """
    :param :α : ℝ
    :param :n : ℝ³
    """
    α = np.asarray(α, dtype=np.float64)
    n = np.asarray(n, dtype=np.float64)

    dim_0 = α.shape[0]
    dim_1 = N1.shape[0]
    assert α.shape == (dim_0,)
    assert n.shape == (dim_0, 3, )
    assert np.ndim(v) == 0

    sum_0 = np.zeros((3, ))
    for T in range(1, len(n)+1):
        if((T) in N1[v-1]):
            sum_0 += α[T-1] * n[T-1]
    sum_1 = np.zeros((3, ))
    for T in range(1, len(n)+1):
        if((T) in N1[v-1]):
            sum_1 += α[T-1] * n[T-1]
    n_left_parenthesis_v_right_parenthesis = (sum_0) / np.linalg.norm(sum_1, 2)
    return pmp_42ResultType(n_left_parenthesis_v_right_parenthesis)


def generateRandomData():
    v = np.random.randint(10)
    dim_0 = np.random.randint(10)
    dim_1 = np.random.randint(10)
    α = np.random.randn(dim_0)
    n = np.random.randn(dim_0, 3, )
    N1 = []
    for i in range(dim_1):
        N1_tmp = []
        dim_2 = np.random.randint(1, 10)
        for j in range(dim_2):
            N1_tmp.append((np.random.randint(10)))
        N1.append(N1_tmp)
    N1 = np.asarray(N1)
    return α, n, v, N1


if __name__ == '__main__':
    α, n, v, N1 = generateRandomData()
    print("α:", α)
    print("n:", n)
    print("v:", v)
    print("N1:", N1)
    func_value = pmp_42(α, n, v, N1)
    print("return value: ", func_value.n_left_parenthesis_v_right_parenthesis)

function output = pmp_42(alpha, n, v, N1)
% output = pmp_42(α, n, v, `N₁`)
%
%    given
%    α_T : ℝ
%    n_T : ℝ³
%    
%    `n(v)` = ( ∑_(T for T ∈ `N₁`_v) α_T n_T )/‖ ∑_(T for T ∈ `N₁`_v) α_T n_T ‖
%    
%    where
%    v ∈ ℤ
%    `N₁`_i ∈ {ℤ}
    if nargin==0
        warning('generating random input data');
        [alpha, n, v, N1] = generateRandomData();
    end
    function [alpha, n, v, N1] = generateRandomData()
        v = randi(10);
        dim_0 = randi(10);
        dim_1 = randi(10);
        alpha = randn(dim_0,1);
        n = randn(dim_0,3);
        N1 = {};
        for i = 1:dim_1
            N1_tmp = [];
            dim_4 = randi(10);
            for j = 1:dim_4 
                N1_tmp = [N1_tmp;randi(10)];
            end
            N1 = [N1, N1_tmp];
        end
    end

    alpha = reshape(alpha,[],1);

    dim_0 = size(alpha, 1);
    dim_1 = size(N1, 1);
    assert( size(alpha,1) == dim_0 );
    assert( isequal(size(n), [dim_0, 3]) );
    assert(numel(v) == 1);

    sum_0 = zeros(3,1);
    for T = 1:size(n, 1)
        if ismember([T],N1{v},'rows')
          sum_0 = sum_0 + alpha(T) * n(T,:)';
        end
    end
    sum_1 = zeros(3,1);
    for T = 1:size(n, 1)
        if ismember([T],N1{v},'rows')
          sum_1 = sum_1 + alpha(T) * n(T,:)';
        end
    end
    n_v = (sum_0) / norm(sum_1, 2);
    output.n_v = n_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*}
\intertext{given} 
\mathit{α}_{\mathit{T}} & \in \mathbb{R} \\
\mathit{n}_{\mathit{T}} & \in \mathbb{R}^{ 3} \\
\\
\textit{n(v)} & = \frac{\sum_{\mathit{T} \in \textit{N₁}_{ \mathit{v} } } \mathit{α}_{ \mathit{T} }\mathit{n}_{ \mathit{T} }}{\left\|\sum_{\mathit{T} \in \textit{N₁}_{ \mathit{v} } } \mathit{α}_{ \mathit{T} }\mathit{n}_{ \mathit{T} }\right\|_2} \\
\intertext{where} 
\mathit{v} & \in \mathbb{Z} \\
\textit{N₁}_{\mathit{i}} & \in \{\mathbb{Z}\} \\
\\
\end{align*}
\end{minipage}
}
\end{center}

\end{document}