In [ ]:

    

import numpy as np
from menpo.shape import UndirectedGraph
from scipy.sparse import lil_matrix
from menpo.visualize import print_dynamic, progress_bar_str


    

In [ ]:

    

def _covariance_matrix_inverse(cov_mat, n_appearance_parameters):
    if n_appearance_parameters is None:
        return np.linalg.inv(cov_mat)
    else:
        s, v, d = np.linalg.svd(cov_mat)
        s = s[:, :n_appearance_parameters]
        v = v[:n_appearance_parameters]
        d = d[:n_appearance_parameters, :]
        return s.dot(np.diag(1/v)).dot(d)


    

In [ ]:

    

graph = UndirectedGraph(np.array([[0, 1], [0, 3], [2, 3], [0, 2]]))

patch_len = 2

n_appearance_parameters = None

all_patches_array = np.random.rand(patch_len * graph.n_vertices, 200)

# initialize block sparse covariance matrix
all_cov = lil_matrix((graph.n_vertices * patch_len,
                      graph.n_vertices * patch_len))

# compute covariance matrix for each edge
for e in range(graph.n_edges):
    # edge vertices
    v1 = np.min(graph.adjacency_array[e, :])
    v2 = np.max(graph.adjacency_array[e, :])

    # find indices in target covariance matrix
    v1_from = v1 * patch_len
    v1_to = (v1 + 1) * patch_len
    v2_from = v2 * patch_len
    v2_to = (v2 + 1) * patch_len

    # extract data
    edge_data = np.concatenate((all_patches_array[v1_from:v1_to, :],
                                all_patches_array[v2_from:v2_to, :]))

    # compute covariance
    cov_mat = _covariance_matrix_inverse(np.cov(edge_data), n_appearance_parameters)

    # v1, v2
    all_cov[v1_from:v1_to, v2_from:v2_to] += cov_mat[:patch_len, patch_len::]
    
    # v2, v1
    all_cov[v2_from:v2_to, v1_from:v1_to] += cov_mat[patch_len::, :patch_len]

    # v1, v1
    all_cov[v1_from:v1_to, v1_from:v1_to] += cov_mat[:patch_len, :patch_len]
    
    # v2, v2
    all_cov[v2_from:v2_to, v2_from:v2_to] += cov_mat[patch_len::, patch_len::]

all_cov = np.array(all_cov.todense())


    

In [ ]:

    

a = np.random.rand(patch_len * graph.n_vertices, 1)
print a.T.dot(all_cov).dot(a)


    

In [ ]:

    

res = 0
for e in range(graph.n_edges):
    # edge vertices
    v1 = np.min(graph.adjacency_array[e, :])
    v2 = np.max(graph.adjacency_array[e, :])
    
    # find indices in target covariance matrix
    v1_from = v1 * patch_len
    v1_to = (v1 + 1) * patch_len
    v2_from = v2 * patch_len
    v2_to = (v2 + 1) * patch_len
    
    # extract data
    edge_data = np.concatenate((all_patches_array[v1_from:v1_to, :],
                                all_patches_array[v2_from:v2_to, :]))
    
    # compute covariance
    cov_mat = _covariance_matrix_inverse(np.cov(edge_data), n_appearance_parameters)
    a_tmp = np.concatenate((a[v1_from:v1_to],
                            a[v2_from:v2_to]))
    
    res += a_tmp.T.dot(cov_mat).dot(a_tmp)
print res


    

In [ ]:

    

print a.todense()


    

In [ ]:

    

all_patches_array[5:10, :].shape


    

In [ ]:

    

import numpy as np
print np.cov(np.random.rand(6, 100))


    

In [ ]: