laydi/fluents/lib/validation.py

from scipy import ones,mean,sqrt,dot,newaxis,zeros,sum,empty,\
     apply_along_axis,eye, kron
from scipy.linalg import triu,inv,svd,norm

from select_generators import w_pls_gen,w_pls_gen_jk,pls_gen,pca_gen,diag_pert
from engines import w_simpls,pls, bridge,pca
from pylab import *

def w_pls_cv_val(X, Y, amax, n_blocks=None, algo='simpls'):
    """RMSEP calc for pls with wide X.
    """
    k, l = Y.shape
    PRESS = zeros((l, amax+1), dtype='f')
    # X,Y are centered
    if n_blocks==None:
        n_blocks = Y.shape[0]
    V = w_pls_gen(dot(X, X.T), Y, n_blocks=n_blocks, center=True)
    for Din, Doi, Yin, Yout in V:
        ym = -sum(Yout, 0)[newaxis]/(1.0*Yin.shape[0])
        Yin = Yin - ym
        PRESS[:,0] = PRESS[:,0] + ((Yout - ym)**2).sum(0)
        if algo=='simpls':
            dat = w_simpls(Din, Yin, amax)
            Q,U,H = dat['Q'], dat['U'], dat['H']
            That = dot(Doi, dot(U, inv(triu(dot(H.T,U))) ))
        else:
            "Other algo-support comming soon"
            raise NotImplementedError
        #Yhat = empty((amax, k, l),dtype='<f8')
        Yhat = []
        for j in range(l):
            TQ = dot(That, triu(dot(Q[j,:][:,newaxis], ones((1,amax)))) )
            E = Yout[:,j][:,newaxis] - TQ
            E = E + sum(E, 0)/Din.shape[0]
            PRESS[j,1:] = PRESS[j,1:] + sum(E**2, 0)
    #Yhat = Y - dot(That,Q.T)
    return sqrt(PRESS/Y.shape[0])

def pls_val(X, Y, amax=2, n_blocks=10,algo='pls'):
    """ Validation results of pls model. 
    """
    
    k, l = Y.shape
    PRESS = zeros((l, amax+1), dtype='<f8')
    EE = zeros((amax, k, l), dtype='<f8')
    Yhat = zeros((amax, k, l), dtype='<f8')
    # X,Y are centered
    V = pls_gen(X, Y, n_blocks=n_blocks, center=True, index_out=True)
    for Xin, Xout, Yin, Yout, out in V:
        ym = -sum(Yout,0)[newaxis]/Yin.shape[0]
        Yin = (Yin - ym)
        PRESS[:,0] = PRESS[:,0] + ((Yout - ym)**2).sum(0)
        if algo=='pls':
            dat = pls(Xin, Yin, amax, mode='normal')
        elif algo=='bridge':
            dat = simpls(Xin, Yin, amax, mode='normal')
        
        for a in range(amax):
            Ba = dat['B'][a,:,:]
            Yhat[a,out[:],:] = dot(Xout, Ba)
            E = Yout -  dot(Xout, Ba)
            EE[a,out,:] = E
            PRESS[:,a+1] = PRESS[:,a+1] + sum(E**2,0)

    return sqrt(PRESS/(k-1.)), EE, Yhat

def pca_alter_val(a, amax, n_sets=10,method='diag'):
    """Pca validation by altering elements in X.
    """
    # todo: it is just as easy to do jk-estimates her as well
    V = diag_pert(a, n_sets, center=True, index_out=True)
    sep = empty((n_sets, amax), dtype='f')
    for i, (xi, ind) in enumerate(V):
        dat_i = pca(xi, amax, mode='detailed')
        Ti,Pi = dat_i['T'],dat_i['P']
        for j in xrange(amax):
            Xhat = dot(Ti[:,:j+1], Pi[:,:j+1].T)
            a_sub = a.ravel().take(ind)
            EE = a_sub - Xhat.ravel().take(ind)
            tot = (a_sub**2).sum()
            sep[i,j] = (EE**2).sum()/tot
    return sqrt(sep.mean(0))
    #return sep

def pca_cv_val(X, amax, n_sets):
    """ Cross validation of pca using random sets crossval.
    """
    m, n = X.shape
    xtot = (X**2).sum()
    V = pca_gen(X, n_sets=7, center=True, index_out=True)
    E = empty((amax, m, n), dtype='f')
    for xi,xout,ind in V:
        dat_i = pca(xi, amax, mode='detailed')
        Pi = dat_i['P']
        for a in xrange(amax):
            Pia = Pi[:,:a+1]
            E[a][ind,:] = (X[ind,:] - dot(xout, dot(Pia,Pia.T) ))**2

    sep = []
    for a in xrange(amax):
        sep.append(E[a].sum()/xtot)
    return sqrt(sep.mean(0))

def pls_jkW(a, b, amax, n_blocks=None, algo='pls', use_pack=True):
    """ Returns CV-segments of paramter W for wide X.

    todo: add support for T,Q and B
    """
    if n_blocks == None:
        n_blocks = b.shape[0]

    WW = empty((n_blocks, a.shape[1], amax), dtype='f')

    if use_pack:
        u, s, inflater = svd(a, full_matrices=0)
        a = u*s
    V = pls_gen(a, b, n_blocks=n_blocks)
    for nn,(a_in, a_out, b_in, b_out) in enumerate(V):
        if algo=='pls':
            dat = pls(a_in, b_in, amax, 'loads', 'fast')
        elif algo=='bridge':
            dat = bridge(a_in, b_in, amax, 'loads', 'fast')
        W = dat['W']
        if use_pack:
            W = dot(inflater.T, W)
        WW[nn,:,:] = W
        
    return WW

def pca_jkP(a, aopt, n_blocks=None):
    """ Returns CV-segments of paramter P.
    todo: add support for T
    fixme: more efficient to add this in validation loop
    """
    if n_blocks == None:
        n_blocks = a.shape[0]

    PP = empty((n_blocks, a.shape[1], aopt), dtype='f')
    V = pca_gen(a, n_sets=n_blocks, center=True)
    for nn,(a_in, a_out) in enumerate(V):  
        dat = pca(a_in, aopt, mode='fast')
        P = dat['P']
        PP[nn,:,:] = P
        
    return PP