Atlas Refinement with Bounded Packing Efficiency Eq. 3

`G_σ(s^k_i)` = ∑_j l_j exp(-dist(`bᵢ`, b_j)²/(2σ²)) `s^k`_j

where
l_j ∈ ℝ : the length of bj
dist: ℝ², ℝ² → ℝ : measures the geodesic distance 
`bᵢ` ∈ ℝ²
b_j ∈ ℝ²
σ ∈ ℝ
`s^k`_j ∈ ℝ² : direction vector

/*
`G_σ(s^k_i)` = ∑_j l_j exp(-dist(`bᵢ`, b_j)²/(2σ²)) `s^k`_j

where
l_j ∈ ℝ : the length of bj
dist: ℝ², ℝ² → ℝ : measures the geodesic distance 
`bᵢ` ∈ ℝ²
b_j ∈ ℝ²
σ ∈ ℝ
`s^k`_j ∈ ℝ² : direction vector
*/
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <iostream>
#include <set>

struct atlas_refinement_3ResultType {
    Eigen::Matrix<double, 2, 1> G_σ_left_parenthesis_s_circumflex_accent_k_i_right_parenthesis;
    atlas_refinement_3ResultType(const Eigen::Matrix<double, 2, 1> & G_σ_left_parenthesis_s_circumflex_accent_k_i_right_parenthesis)
    : G_σ_left_parenthesis_s_circumflex_accent_k_i_right_parenthesis(G_σ_left_parenthesis_s_circumflex_accent_k_i_right_parenthesis)
    {}
};

/**
 * atlas_refinement_3
 *
 * @param l  the length of bj
 * @param dist  ℝ², ℝ² → ℝ : measures the geodesic distance 
 * @param s_circumflex_accent_k  direction vector
 * @return G_σ_left_parenthesis_s_circumflex_accent_k_i_right_parenthesis
 */
atlas_refinement_3ResultType atlas_refinement_3(
    const std::vector<double> & l,
    const std::function<double(Eigen::Matrix<double, 2, 1>, Eigen::Matrix<double, 2, 1>)> & dist,
    const Eigen::Matrix<double, 2, 1> & bᵢ,
    const std::vector<Eigen::Matrix<double, 2, 1>> & b,
    const double & σ,
    const std::vector<Eigen::Matrix<double, 2, 1>> & s_circumflex_accent_k)
{
    const long dim_0 = l.size();
    assert( b.size() == dim_0 );
    assert( s_circumflex_accent_k.size() == dim_0 );

    Eigen::MatrixXd sum_0 = Eigen::MatrixXd::Zero(2, 1);
    for(int j=1; j<=s_circumflex_accent_k.size(); j++){
        sum_0 += l.at(j-1) * exp(-pow(dist(bᵢ, b.at(j-1)), 2) / double((2 * pow(σ, 2)))) * s_circumflex_accent_k.at(j-1);
    }
    Eigen::Matrix<double, 2, 1> G_σ_left_parenthesis_s_circumflex_accent_k_i_right_parenthesis = sum_0;

    return atlas_refinement_3ResultType(G_σ_left_parenthesis_s_circumflex_accent_k_i_right_parenthesis);
}


void generateRandomData(std::vector<double> & l,
    std::function<double(Eigen::Matrix<double, 2, 1>, Eigen::Matrix<double, 2, 1>)> & dist,
    Eigen::Matrix<double, 2, 1> & bᵢ,
    std::vector<Eigen::Matrix<double, 2, 1>> & b,
    double & σ,
    std::vector<Eigen::Matrix<double, 2, 1>> & s_circumflex_accent_k)
{
    σ = rand() % 10;
    const int dim_0 = rand()%10;
    l.resize(dim_0);
    for(int i=0; i<dim_0; i++){
        l[i] = rand() % 10;
    }
    dist = [](Eigen::Matrix<double, 2, 1>, Eigen::Matrix<double, 2, 1>)->double{
        return rand() % 10;
    };
    bᵢ = Eigen::VectorXd::Random(2);
    b.resize(dim_0);
    for(int i=0; i<dim_0; i++){
        b[i] = Eigen::VectorXd::Random(2);
    }
    s_circumflex_accent_k.resize(dim_0);
    for(int i=0; i<dim_0; i++){
        s_circumflex_accent_k[i] = Eigen::VectorXd::Random(2);
    }
}


int main(int argc, char *argv[])
{
    srand((int)time(NULL));
    std::vector<double> l;
    std::function<double(Eigen::Matrix<double, 2, 1>, Eigen::Matrix<double, 2, 1>)> dist;
    Eigen::Matrix<double, 2, 1> bᵢ;
    std::vector<Eigen::Matrix<double, 2, 1>> b;
    double σ;
    std::vector<Eigen::Matrix<double, 2, 1>> s_circumflex_accent_k;
    generateRandomData(l, dist, bᵢ, b, σ, s_circumflex_accent_k);
    atlas_refinement_3ResultType func_value = atlas_refinement_3(l, dist, bᵢ, b, σ, s_circumflex_accent_k);
    std::cout<<"return value:\n"<<func_value.G_σ_left_parenthesis_s_circumflex_accent_k_i_right_parenthesis<<std::endl;
    return 0;
}

"""
`G_σ(s^k_i)` = ∑_j l_j exp(-dist(`bᵢ`, b_j)²/(2σ²)) `s^k`_j

where
l_j ∈ ℝ : the length of bj
dist: ℝ², ℝ² → ℝ : measures the geodesic distance 
`bᵢ` ∈ ℝ²
b_j ∈ ℝ²
σ ∈ ℝ
`s^k`_j ∈ ℝ² : direction vector
"""
import numpy as np
import scipy
import scipy.linalg
from scipy import sparse
from scipy.integrate import quad
from scipy.optimize import minimize


class atlas_refinement_3ResultType:
    def __init__( self, G_σ_left_parenthesis_s_circumflex_accent_k_i_right_parenthesis):
        self.G_σ_left_parenthesis_s_circumflex_accent_k_i_right_parenthesis = G_σ_left_parenthesis_s_circumflex_accent_k_i_right_parenthesis


def atlas_refinement_3(l, dist, bᵢ, b, σ, s_circumflex_accent_k):
    """
    :param :l : the length of bj
    :param :dist : ℝ², ℝ² → ℝ : measures the geodesic distance 
    :param :s_circumflex_accent_k : direction vector
    """
    l = np.asarray(l, dtype=np.float64)
    bᵢ = np.asarray(bᵢ, dtype=np.float64)
    b = np.asarray(b, dtype=np.float64)
    s_circumflex_accent_k = np.asarray(s_circumflex_accent_k, dtype=np.float64)

    dim_0 = l.shape[0]
    assert l.shape == (dim_0,)
    assert bᵢ.shape == (2,)
    assert b.shape == (dim_0, 2, )
    assert np.ndim(σ) == 0
    assert s_circumflex_accent_k.shape == (dim_0, 2, )

    sum_0 = np.zeros((2, ))
    for j in range(1, len(s_circumflex_accent_k)+1):
        sum_0 += l[j-1] * np.exp(-np.power(dist(bᵢ, b[j-1]), 2) / (2 * np.power(σ, 2))) * s_circumflex_accent_k[j-1]
    G_σ_left_parenthesis_s_circumflex_accent_k_i_right_parenthesis = sum_0
    return atlas_refinement_3ResultType(G_σ_left_parenthesis_s_circumflex_accent_k_i_right_parenthesis)


def generateRandomData():
    σ = np.random.randn()
    dim_0 = np.random.randint(10)
    l = np.random.randn(dim_0)
    def dist(p0, p1):
        return np.random.randn()
    bᵢ = np.random.randn(2)
    b = np.random.randn(dim_0, 2, )
    s_circumflex_accent_k = np.random.randn(dim_0, 2, )
    return l, dist, bᵢ, b, σ, s_circumflex_accent_k


if __name__ == '__main__':
    l, dist, bᵢ, b, σ, s_circumflex_accent_k = generateRandomData()
    print("l:", l)
    print("dist:", dist)
    print("bᵢ:", bᵢ)
    print("b:", b)
    print("σ:", σ)
    print("s_circumflex_accent_k:", s_circumflex_accent_k)
    func_value = atlas_refinement_3(l, dist, bᵢ, b, σ, s_circumflex_accent_k)
    print("return value: ", func_value.G_σ_left_parenthesis_s_circumflex_accent_k_i_right_parenthesis)

function output = atlas_refinement_3(l, dist, b_i, b, sigma, s_k)
% output = atlas_refinement_3(l, dist, `bᵢ`, b, σ, `s^k`)
%
%    `G_σ(s^k_i)` = ∑_j l_j exp(-dist(`bᵢ`, b_j)²/(2σ²)) `s^k`_j
%    
%    where
%    l_j ∈ ℝ : the length of bj
%    dist: ℝ², ℝ² → ℝ : measures the geodesic distance 
%    `bᵢ` ∈ ℝ²
%    b_j ∈ ℝ²
%    σ ∈ ℝ
%    `s^k`_j ∈ ℝ² : direction vector
    if nargin==0
        warning('generating random input data');
        [l, dist, b_i, b, sigma, s_k] = generateRandomData();
    end
    function [l, dist, b_i, b, sigma, s_k] = generateRandomData()
        sigma = randn();
        dim_0 = randi(10);
        l = randn(dim_0,1);
        dist = @distFunc;
        rseed = randi(2^32);
        function tmp =  distFunc(p0, p1)
            rng(rseed);
            tmp = randn();
        end

        b_i = randn(2,1);
        b = randn(dim_0,2);
        s_k = randn(dim_0,2);
    end

    l = reshape(l,[],1);
    b_i = reshape(b_i,[],1);

    dim_0 = size(l, 1);
    assert( size(l,1) == dim_0 );
    assert( numel(b_i) == 2 );
    assert( isequal(size(b), [dim_0, 2]) );
    assert(numel(sigma) == 1);
    assert( isequal(size(s_k), [dim_0, 2]) );

    sum_0 = zeros(2,1);
    for j = 1:size(s_k, 1)
        sum_0 = sum_0 + l(j) * exp(-dist(b_i, b(j,:)').^2 / (2 * sigma.^2)) * s_k(j,:)';
    end
    G_sigma_s_k_i = sum_0;
    output.G_sigma_s_k_i = G_sigma_s_k_i;
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*}
\textit{G\_σ(s\^k\_i)} & = \sum_\mathit{j} \mathit{l}_{ \mathit{j} }exp\left( -\frac{{\mathit{dist}\left( \textit{b\textsubscript{i}},\mathit{b}_{ \mathit{j} } \right)}^{2}}{2{\mathit{σ}}^{2}} \right)\textit{s\^k}_{ \mathit{j} } \\
\intertext{where} 
\mathit{l}_{\mathit{j}} & \in \mathbb{R} \text{ the length of bj} \\
\mathit{dist} & \in \mathbb{R}^{ 2},\mathbb{R}^{ 2}\rightarrow \mathbb{R} \text{ measures the geodesic distance } \\
\textit{b\textsubscript{i}} & \in \mathbb{R}^{ 2} \\
\mathit{b}_{\mathit{j}} & \in \mathbb{R}^{ 2} \\
\mathit{σ} & \in \mathbb{R} \\
\textit{s\^k}_{\mathit{j}} & \in \mathbb{R}^{ 2} \text{ direction vector} \\
\\
\end{align*}
\end{minipage}
}
\end{center}

\end{document}