Convex Optimization page 208

`I(X;Y)` = ∑_i ∑_j x_j p_i,j log₂(p_i,j/∑_k x_k p_i,k)

where

x ∈ ℝ^n
p ∈ ℝ^(m×n)

/*
`I(X;Y)` = ∑_i ∑_j x_j p_i,j log₂(p_i,j/∑_k x_k p_i,k)

where

x ∈ ℝ^n
p ∈ ℝ^(m×n)
*/
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <iostream>
#include <set>

struct convex_optimization_208ResultType {
    double I_left_parenthesis_X_semicolon_Y_right_parenthesis;
    convex_optimization_208ResultType(const double & I_left_parenthesis_X_semicolon_Y_right_parenthesis)
    : I_left_parenthesis_X_semicolon_Y_right_parenthesis(I_left_parenthesis_X_semicolon_Y_right_parenthesis)
    {}
};

convex_optimization_208ResultType convex_optimization_208(
    const Eigen::VectorXd & x,
    const Eigen::MatrixXd & p)
{
    const long n = x.size();
    const long m = p.rows();
    assert( p.cols() == n );

    double sum_0 = 0;
    for(int i=1; i<=p.rows(); i++){
        double sum_1 = 0;
        for(int j=1; j<=x.size(); j++){
            double sum_2 = 0;
            for(int k=1; k<=p.cols(); k++){
                sum_2 += x[k-1] * p(i-1, k-1);
            }
            sum_1 += x[j-1] * p(i-1, j-1) * (log10(p(i-1, j-1) / double(sum_2)) / log10(2));
        }
        sum_0 += sum_1;
    }
    double I_left_parenthesis_X_semicolon_Y_right_parenthesis = sum_0;

    return convex_optimization_208ResultType(I_left_parenthesis_X_semicolon_Y_right_parenthesis);
}


void generateRandomData(Eigen::VectorXd & x,
    Eigen::MatrixXd & p)
{
    const int n = rand()%10;
    const int m = rand()%10;
    x = Eigen::VectorXd::Random(n);
    p = Eigen::MatrixXd::Random(m, n);
}


int main(int argc, char *argv[])
{
    srand((int)time(NULL));
    Eigen::VectorXd x;
    Eigen::MatrixXd p;
    generateRandomData(x, p);
    convex_optimization_208ResultType func_value = convex_optimization_208(x, p);
    std::cout<<"return value:\n"<<func_value.I_left_parenthesis_X_semicolon_Y_right_parenthesis<<std::endl;
    return 0;
}

"""
`I(X;Y)` = ∑_i ∑_j x_j p_i,j log₂(p_i,j/∑_k x_k p_i,k)

where

x ∈ ℝ^n
p ∈ ℝ^(m×n)
"""
import numpy as np
import scipy
import scipy.linalg
from scipy import sparse
from scipy.integrate import quad
from scipy.optimize import minimize


class convex_optimization_208ResultType:
    def __init__( self, I_left_parenthesis_X_semicolon_Y_right_parenthesis):
        self.I_left_parenthesis_X_semicolon_Y_right_parenthesis = I_left_parenthesis_X_semicolon_Y_right_parenthesis


def convex_optimization_208(x, p):
    x = np.asarray(x, dtype=np.float64)
    p = np.asarray(p, dtype=np.float64)

    n = x.shape[0]
    m = p.shape[0]
    assert x.shape == (n,)
    assert p.shape == (m, n)

    sum_0 = 0
    for i in range(1, p.shape[0]+1):
        sum_1 = 0
        for j in range(1, len(x)+1):
            sum_2 = 0
            for k in range(1, p.shape[1]+1):
                sum_2 += x[k-1] * p[i-1, k-1]
            sum_1 += x[j-1] * p[i-1, j-1] * np.log2(p[i-1, j-1] / sum_2)
        sum_0 += sum_1
    I_left_parenthesis_X_semicolon_Y_right_parenthesis = sum_0
    return convex_optimization_208ResultType(I_left_parenthesis_X_semicolon_Y_right_parenthesis)


def generateRandomData():
    n = np.random.randint(10)
    m = np.random.randint(10)
    x = np.random.randn(n)
    p = np.random.randn(m, n)
    return x, p


if __name__ == '__main__':
    x, p = generateRandomData()
    print("x:", x)
    print("p:", p)
    func_value = convex_optimization_208(x, p)
    print("return value: ", func_value.I_left_parenthesis_X_semicolon_Y_right_parenthesis)

function output = convex_optimization_208(x, p)
% output = convex_optimization_208(x, p)
%
%    `I(X;Y)` = ∑_i ∑_j x_j p_i,j log₂(p_i,j/∑_k x_k p_i,k)
%    
%    where
%    
%    x ∈ ℝ^n
%    p ∈ ℝ^(m×n)
    if nargin==0
        warning('generating random input data');
        [x, p] = generateRandomData();
    end
    function [x, p] = generateRandomData()
        n = randi(10);
        m = randi(10);
        x = randn(n,1);
        p = randn(m, n);
    end

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

    n = size(x, 1);
    m = size(p, 1);
    assert( numel(x) == n );
    assert( isequal(size(p), [m, n]) );

    sum_0 = 0;
    for i = 1:size(p,1)
        sum_1 = 0;
        for j = 1:size(x,1)
            sum_2 = 0;
            for k = 1:size(p,2)
                sum_2 = sum_2 + x(k) * p(i, k);
            end
            sum_1 = sum_1 + x(j) * p(i, j) * log2(p(i, j) / sum_2);
        end
        sum_0 = sum_0 + sum_1;
    end
    I_X_semicolon_Y = sum_0;
    output.I_X_semicolon_Y = I_X_semicolon_Y;
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{I(X;Y)} & = \sum_\mathit{i} \sum_\mathit{j} \mathit{x}_{ \mathit{j} }\mathit{p}_{\mathit{i}, \mathit{j}} \log_2{ \frac{\mathit{p}_{\mathit{i}, \mathit{j}}}{\sum_\mathit{k} \mathit{x}_{ \mathit{k} }\mathit{p}_{\mathit{i}, \mathit{k}}} } \\
\intertext{where} 
\mathit{x} & \in \mathbb{R}^{ \mathit{n}} \\
\mathit{p} & \in \mathbb{R}^{ \mathit{m} \times \mathit{n} } \\
\\
\end{align*}
\end{minipage}
}
\end{center}

\end{document}