# -*- coding: utf-8 -*-
"""
Created on Thu Jun 29 11:39:00 2017

@author: durand

ACP probabiliste :
Calcule l'inverse de la matrice de variance-covariance de la loi de W|X,Z, ...
Modele : X = mu + WZ + sigma epsilon
reecrit X = mu + UW  + sigma epsilon
V = U * U.T + sigma**2 I
"""

import scipy, scipy.linalg

def CalculU(Z,D,sigma2=0):
    """
    Calcule U
    """
    n, M = Z.shape
    U = scipy.matrix(scipy.zeros((n*D,D*M)))
    for i in range(n):
        for j in range(D):
            # print(i*D+j)
            # print(range((j*M),(j*M+M)))
            U[i*D+j,(j*M):(j*M+M)] = Z[i,:]
    return(U)
    
#def InvV(Z,D,sigma2=0):
#    """
#    Calcule V^-1    
#    """
#    
    
    
# Exemple
n = 3
M = 4
D = 5
vals = scipy.randn(n*M)
Z = scipy.matrix(vals)
Z = Z.reshape(n,M)
U = CalculU(Z,D,0)
V = U*U.T
VM1 = scipy.linalg.inv(V)
print(VM1)
Sigma = Z*Z.T
print(scipy.linalg.inv(Sigma))
UTVM1 = U.T * VM1
print(UTVM1)
UTVM1U = UTVM1 * U
print(UTVM1U)
UTVM1U[0:M,0:M]
UTVM1U[M:2*M,M:2*M]
Z.T * scipy.linalg.inv(Sigma) * Z

# On remarque que la matrice de variance-covariance de la loi de W|X,Z est
# diagonale par blocs, les blocs diagonaux etant egaux a 
# Z.T * (Z * Z.T + sigma^2 * I)^{-1} * Z