laydi/fluents/lib/engines.py

"""Module contain algorithms for low-rank models.

There is almost no typechecking of any kind here, just focus on speed
"""

import math
import warnings
from scipy.linalg import svd,inv
from scipy import dot,empty,eye,newaxis,zeros,sqrt,diag,\
     apply_along_axis,mean,ones,randn,empty_like,outer,r_,c_,\
     rand,sum,cumsum,matrix, expand_dims,minimum,where,arange,inner,tile
has_sym = True
has_arpack = True
try:
    from symeig import symeig
except:
    has_sym = False
try:
    from scipy.sandbox import arpack
except:
    has_arpack = False


def pca(a, aopt,scale='scores',mode='normal',center_axis=0):
    """ Principal Component Analysis.

    Performs PCA on given matrix and returns results in a dictionary.

    :Parameters:
        a : array
    Data measurement matrix, (samples x variables)
        aopt : int
    Number of components to use, aopt<=min(samples, variables)

    :Returns:
    results : dict
        keys -- values,  T -- scores, P -- loadings, E -- residuals,
        lev --leverages, ssq -- sum of squares, expvar -- cumulative
        explained variance, aopt -- number of components used
    
    :OtherParam eters:
    mode : str
        Amount of info retained, ('fast', 'normal', 'detailed')
    center_axis : int
        Center along given axis. If neg.: no centering (-inf,..., matrix modes)
  
    :SeeAlso:
        - pcr : other blm
        - pls : other blm
        - lpls : other blm

    Notes
    -----
    Uses kernel speed-up if m>>n or m<<n.

    If residuals turn rank deficient, a lower number of component than given
    in input will be used. The number of components used is given in
    results-dict.
    
    Examples
    --------
    
    >>> import scipy,engines
    >>> a=scipy.asarray([[1,2,3],[2,4,5]])
    >>> dat=engines.pca(a, 2)
    >>> dat['expvarx']
    array([0.,99.8561562,  100.])

    """
    m, n = a.shape
    assert(aopt<=min(m,n))
    if center_axis>=0:
        a = a - expand_dims(a.mean(center_axis), center_axis)
    if m>(n+100) or n>(m+100):
        u, s, v = esvd(a, amax=None) # fixme:amax option need to work with expl.var
    else:
        u, s, vt = svd(a, 0)
        v = vt.T
    e = s**2
    tol = 1e-10
    eff_rank = sum(s>s[0]*tol)
    aopt = minimum(aopt, eff_rank)
    T = u*s
    s = s[:aopt]
    T = T[:,:aopt]
    P = v[:,:aopt]
    
    if scale=='loads':
        T = T/s
        P = P*s

    if mode == 'fast':
        return {'T':T, 'P':P, 'aopt':aopt}
    
    if mode=='detailed':
        E = empty((aopt, m, n))
        ssq = []
        lev = []
        for ai in range(aopt):
            E[ai,:,:] = a - dot(T[:,:ai+1], P[:,:ai+1].T)
            ssq.append([(E[ai,:,:]**2).mean(0), (E[ai,:,:]**2).mean(1)])
            if scale=='loads':
                lev.append([((s*T)**2).sum(1), (P**2).sum(1)])
            else:
                lev.append([(T**2).sum(1), ((s*P)**2).sum(1)])
    else:
        # residuals
        E = a - dot(T, P.T)
        #E = a
        SEP = E**2
        ssq = [SEP.sum(0), SEP.sum(1)]
        # leverages
        if scale=='loads':
            lev = [(1./m)+(T**2).sum(1), (1./n)+((P/s)**2).sum(1)]
        else:
            lev = [(1./m)+((T/s)**2).sum(1), (1./n)+(P**2).sum(1)]
        # variances
    expvarx = r_[0, 100*e.cumsum()/e.sum()][:aopt+1]
    
    return {'T':T, 'P':P, 'E':E, 'expvarx':expvarx, 'levx':lev, 'ssqx':ssq, 'aopt':aopt, 'eigvals': e[:aopt,newaxis]}

def pcr(a, b, aopt, scale='scores',mode='normal',center_axis=0):
    """ Principal Component Regression.

    Performs PCR on given matrix and returns results in a dictionary.

    :Parameters:
        a : array
    Data measurement matrix, (samples x variables)
        b : array
    Data response matrix, (samples x responses)
        aopt : int
    Number of components to use, aopt<=min(samples, variables)

    :Returns:
    results : dict
        keys -- values,  T -- scores, P -- loadings, E -- residuals,
        levx -- leverages, ssqx -- sum of squares, expvarx -- cumulative
        explained variance, aopt -- number of components used
    
    :OtherParameters:
    mode : str
        Amount of info retained, ('fast', 'normal', 'detailed')
    center_axis : int
        Center along given axis. If neg.: no centering (-inf,..., matrix modes)
  
    :SeeAlso:
        - pca : other blm
        - pls : other blm
        - lpls : other blm
 
    Notes
    -----

    Uses kernel speed-up if m>>n or m<<n.

    If residuals turn rank deficient, a lower number of component than given
    in input will be used. The number of components used is given in results-dict. 
    

    Examples
    --------
    
    >>> import scipy,engines
    >>> a=scipy.asarray([[1,2,3],[2,4,5]])
    >>> b=scipy.asarray([[1,1],[2,3]])
    >>> dat=engines.pcr(a, 2)
    >>> dat['expvarx']
    array([0.,99.8561562,  100.])

    """
    k, l = m_shape(b)
    if center_axis>=0:
        b = b - expand_dims(b.mean(center_axis), center_axis)
    dat = pca(a, aopt=aopt, scale=scale, mode=mode, center_axis=center_axis)
    T = dat['T']
    weights = apply_along_axis(vnorm, 0, T)**2
    if scale=='loads':
        Q = dot(b.T, T*weights)
    else:
        Q = dot(b.T, T/weights)

    if mode=='fast':
        dat.update({'Q':Q})
        return dat
    if mode=='detailed':
        F = empty((aopt, k, l))
        for i in range(aopt):
            F[i,:,:] = b - dot(T[:,:i+1], Q[:,:i+1].T)
    else:
        F = b - dot(T, Q.T)
    expvary = r_[0,  100*((T**2).sum(0)*(Q**2).sum(0)/(b**2).sum()).cumsum()[:aopt]]
    #fixme: Y-var leverages
    dat.update({'Q':Q, 'F':F, 'expvary':expvary})
    return dat

def pls(a, b, aopt=2, scale='scores', mode='normal', center_axis=-1, ab=None):
    """Partial Least Squares Regression.

    Performs PLS on given matrix and returns results in a dictionary.

    :Parameters:
        a : array
    Data measurement matrix, (samples x variables)
        b : array
    Data response matrix, (samples x responses)
        aopt : int
    Number of components to use, aopt<=min(samples, variables)

    :Returns:
    results : dict
        keys -- values,  T -- scores, P -- loadings, E -- residuals,
        levx -- leverages, ssqx -- sum of squares, expvarx -- cumulative
        explained variance of descriptors, expvary -- cumulative explained
        variance of responses, aopt -- number of components used
    
    :OtherParameters:
    mode : str
        Amount of info retained, ('fast', 'normal', 'detailed')
    center_axis : int
        Center along given axis. If neg.: no centering (-inf,..., matrix modes)
    
    :SeeAlso:
        - pca : other blm
        - pcr : other blm
        - lpls : other blm
 
    Notes
    -----

    Uses kernel speed-up if m>>n or m<<n.

    If residuals turn rank deficient, a lower number of component than given
    in input will be used. The number of components used is given in results-dict.

    Examples
    --------
    
    >>> import scipy,engines
    >>> a=scipy.asarray([[1,2,3],[2,4,5]])
    >>> b=scipy.asarray([[1,1],[2,3]])
    >>> dat=engines.pls(a, b, 2)
    >>> dat['expvarx']
    array([0.,99.8561562,  100.])

    """
    
    m, n = m_shape(a)
    if ab!=None:
        mm, l = m_shape(ab)
        assert(m==mm)
    else:
         k, l = m_shape(b)

    if center_axis>=0:
        a = a - expand_dims(a.mean(center_axis), center_axis)
        b = b - expand_dims(b.mean(center_axis), center_axis)
    
    W = empty((n, aopt))
    P = empty((n, aopt))
    R = empty((n, aopt))
    Q = empty((l, aopt))
    T = empty((m, aopt))
    B = empty((aopt, n, l))
    tt = empty((aopt,))

    if ab==None:
        ab = dot(a.T, b)
    for i in range(aopt):
        if ab.shape[1]==1: #pls 1
            w = ab.reshape(n, l)
            w = w/vnorm(w)
        elif n<l: # more yvars than xvars
            if has_sym:
                s, w = symeig(dot(ab, ab.T),range=[n,n],overwrite=True)
            else:
                w, s, vh = svd(dot(ab, ab.T))
            w = w[:,:1]
        else: # standard wide xdata
            if has_sym:
                s, q = symeig(dot(ab.T, ab),range=[l,l],overwrite=True)
            else:
                q, s, vh = svd(dot(ab.T, ab))
                q = q[:,:1]
            w = dot(ab, q)
            w = w/vnorm(w)
        r = w.copy()
        if i>0:
            for j in range(0, i, 1):
                r = r - dot(P[:,j].T, w)*R[:,j][:,newaxis]
        
        t = dot(a, r)
        tt[i] = tti = dot(t.T, t).ravel()
        p  = dot(a.T, t)/tti
        q = dot(r.T, ab).T/tti
        ab = ab - dot(p, q.T)*tti
        T[:,i] = t.ravel()
        W[:,i] = w.ravel()

        if mode=='fast' and i==aopt-1:
            if scale=='loads':
                tnorm = sqrt(tt)
                T = T/tnorm
                W = W*tnorm
            return {'T':T, 'W':W}

        P[:,i] = p.ravel()
        R[:,i] = r.ravel()
        Q[:,i] = q.ravel()
        #B[i] = dot(R[:,:i+1], Q[:,:i+1].T)


        
    qnorm = apply_along_axis(vnorm, 0, Q)
    tnorm = sqrt(tt)
    pp = (P**2).sum(0)
    if mode=='detailed':
        E = empty((aopt, m, n))
        F = empty((aopt, k, l))
        ssqx, ssqy = [], []
        leverage = empty((aopt, m))
        h2x = [] #hotellings T^2
        h2y = []
        for ai in range(aopt):
            E[ai,:,:] = a - dot(T[:,:ai+1], P[:,:ai+1].T)
            F[i-1] = b - dot(T[:,:i], Q[:,:i].T)
            ssqx.append([(E[ai,:,:]**2).mean(0), (E[ai,:,:]**2).mean(1)])
            ssqy.append([(F[ai,:,:]**2).mean(0), (F[ai,:,:]**2).mean(1)])
            leverage[ai,:] = 1./m + ((T[:,:ai+1]/tnorm[:ai+1])**2).sum(1)
            h2y.append(1./k + ((Q[:,:ai+1]/qnorm[:ai+1])**2).sum(1))
    else:
        # residuals
        E = a - dot(T, P.T)
        F = b - dot(T, Q.T)
        sepx = E**2
        ssqx = [sepx.sum(0), sepx.sum(1)]
        sepy = F**2
        ssqy = [sepy.sum(0), sepy.sum(1)]
        # leverage
        leverage = 1./m + ((T/tnorm)**2).sum(1)
        h2x = []
        h2y = []
    # variances
    tp= tt*pp
    tq = tt*qnorm*qnorm
    expvarx = r_[0, 100*tp/(a*a).sum()]
    expvary = r_[0, 100*tq/(b*b).sum()]
    
    if scale=='loads':
        T = T/tnorm
        W = W*tnorm
        Q = Q*tnorm
        P = P*tnorm

    return {'Q':Q, 'P':P, 'T':T, 'W':W, 'R':R, 'E':E, 'F':F,
            'expvarx':expvarx, 'expvary':expvary, 'ssqx':ssqx, 'ssqy':ssqy,
            'leverage':leverage, 'h2':h2x}

def w_simpls(aat, b, aopt):
    """ Simpls for wide matrices.
    Fast pls for crossval, used in calc rmsep for wide X
    There is no P or W.  T is normalised
    """
    bb = b.copy()
    m, m = aat.shape
    U = empty((m, aopt)) # W
    T = empty((m, aopt))
    H = empty((m, aopt)) # R
    PROJ = empty((m, aopt)) # P?

    for i in range(aopt):
        q, s, vh = svd(dot(dot(b.T, aat), b), full_matrices=0)
        u = dot(b, q[:,:1]) #y-factor scores
        U[:,i] = u.ravel()
        t = dot(aat, u)
        t = t/vnorm(t)
        T[:,i] = t.ravel()
        h = dot(aat, t) #score-weights
        H[:,i] = h.ravel()
        PROJ[:,:i+1] = dot(T[:,:i+1], inv(dot(T[:,:i+1].T, H[:,:i+1])) )
        if i<aopt:
            b = b - dot(PROJ[:,:i+1], dot(H[:,:i+1].T,b) )
    C = dot(bb.T, T)

    return {'T':T, 'U':U, 'Q':C, 'H':H}

def w_pls(aat, b, aopt):
    """ Pls for wide matrices.
    Fast pls for crossval, used in calc rmsep for wide X
    There is no P or W.  T is normalised

    aat = centered kernel matrix
    b = centered y
    """
    bb = b.copy()
    k, l = m_shape(b)
    m, m = m_shape(aat)
    U = empty((m, aopt)) # W
    T = empty((m, aopt))
    R = empty((m, aopt)) # R
    PROJ = empty((m, aopt)) # P?

    for i in range(aopt):
        if has_sym:
            s, q = symeig(dot(dot(b.T, aat), b), range=(l,l),overwrite=True)
        else:
            q, s, vh = svd(dot(dot(b.T, aat), b), full_matrices=0)
            q = q[:,:1]
        u = dot(b , q) #y-factor scores
        U[:,i] = u.ravel()
        t = dot(aat, u)
        
        t = t/vnorm(t)
        T[:,i] = t.ravel()
        r = dot(aat, t)#score-weights
        #r = r/vnorm(r)
        R[:,i] = r.ravel()
        PROJ[:,: i+1] = dot(T[:,:i+1], inv(dot(T[:,:i+1].T, R[:,:i+1])) )
        if i<aopt:
            b = b - dot(PROJ[:,:i+1], dot(R[:,:i+1].T,  b) )
    C = dot(bb.T, T)

    return {'T':T, 'U':U, 'Q':C, 'R':R}

def bridge(a, b, aopt, scale='scores', mode='normal', r=0):
    """Undeflated Ridged svd(X'Y)
    """
    m, n = m_shape(a)
    k, l = m_shape(b)
    u, s, vt = svd(b, full_matrices=0)
    g0 = dot(u*s, u.T)
    g = (1 - r)*g0 + r*eye(m)
    ag = dot(a.T, g)
    u, s, vt = svd(ag, full_matrices=0)
    W = u[:,:aopt]
    K = vt[:aopt,:].T
    T = dot(a, W)
    tnorm = apply_along_axis(vnorm, 0, T) # norm of T-columns

    if mode == 'fast':
        if scale=='loads':
            T = T/tnorm
            W = W*tnorm
        return {'T':T, 'W':W}

    U = dot(g0, K) #fixme check this 
    Q = dot(b.T, dot(T, inv(dot(T.T, T)) ))
    B = zeros((aopt, n, l), dtype='f')
    for i in range(aopt):
        B[i] = dot(W[:,:i+1], Q[:,:i+1].T)

    if mode == 'detailed':
        E = empty((aopt, m, n))
        F = empty((aopt, k, l))
        for i in range(aopt):
            E[i] = a - dot(T[:,:i+1], W[:,:i+1].T)
            F[i] = b - dot(a, B[i])
    else: #normal
        F = b - dot(a, B[-1])
        E = a - dot(T, W.T)

    if scale=='loads':
        T = T/tnorm
        W = W*tnorm
        Q = Q*tnorm

    return {'B':B, 'W':W, 'T':T, 'Q':Q, 'E':E, 'F':F, 'U':U, 'P':W}


def nipals_lpls(X, Y, Z, a_max, alpha=.7, mean_ctr=[2, 0, 1], scale='scores', verbose=False):
    """ L-shaped Partial Least Sqaures Regression by the nipals algorithm.

    (X!Z)->Y
    :input:
        X : data matrix (m, n)
        Y : data matrix (m, l)
        Z : data matrix (n, o)

    :output:
      T : X-scores
      W : X-weights/Z-weights
      P : X-loadings
      Q : Y-loadings
      U : X-Y relation
      L : Z-scores
      K : Z-loads
      B : Regression coefficients X->Y
      b0: Regression coefficient intercept
      evx : X-explained variance
      evy : Y-explained variance
      evz : Z-explained variance
      mnx : X location
      mny : Y location
      mnz : Z location

    :Notes:
    
    """
    if mean_ctr!=None:
        xctr, yctr, zctr = mean_ctr
        X, mnX = center(X, xctr)
        Y, mnY = center(Y, yctr)
        Z, mnZ = center(Z, zctr)

    varX = (X**2).sum()
    varY = (Y**2).sum()
    varZ = (Z**2).sum()
    
    m, n = X.shape
    k, l = Y.shape
    u, o = Z.shape

    # initialize 
    U = empty((k, a_max))
    Q = empty((l, a_max))
    T = empty((m, a_max))
    W = empty((n, a_max))
    P = empty((n, a_max))
    K = empty((o, a_max))
    L = empty((u, a_max))
    B = empty((a_max, n, l))
    #b0 = empty((a_max, 1, l))
    var_x = empty((a_max,))
    var_y = empty((a_max,))
    var_z = empty((a_max,))

    MAX_ITER = 250
    LIM = 1e-1
    for a in range(a_max):
        if verbose:
            print "\nWorking on comp. %s" %a
        u = Y[:,:1]
        diff = 1
        niter = 0
        while (diff>LIM and niter<MAX_ITER):
            niter += 1
            u1 = u.copy()
            w = dot(X.T, u)
            w = w/sqrt(dot(w.T, w))
            #w = w/dot(w.T, w)
            l = dot(Z, w)
            k = dot(Z.T, l)
            k = k/sqrt(dot(k.T, k))
            #k = k/dot(k.T, k)
            w = alpha*k + (1-alpha)*w
            #print sqrt(dot(w.T, w))
            w = w/sqrt(dot(w.T, w))
            t = dot(X, w)
            c = dot(Y.T, t)
            c = c/sqrt(dot(c.T, c))
            u = dot(Y, c)
            diff = dot((u-u1).T, (u-u1))
        if verbose:
            print "Converged after %s iterations" %niter
            print "Error: %.2E" %diff
        tt = dot(t.T, t)
        p = dot(X.T, t)/tt
        q = dot(Y.T, t)/tt
        l = dot(Z, w)
        
        U[:,a] = u.ravel()
        W[:,a] = w.ravel()
        P[:,a] = p.ravel()
        T[:,a] = t.ravel()
        Q[:,a] = q.ravel()
        L[:,a] = l.ravel()
        K[:,a] = k.ravel()

        X = X - dot(t, p.T)
        Y = Y - dot(t, q.T)
        Z = (Z.T - dot(w, l.T)).T

        var_x[a] = pow(X, 2).sum()
        var_y[a] = pow(Y, 2).sum()
        var_z[a] = pow(Z, 2).sum()
    
        B[a] = dot(dot(W[:,:a+1], inv(dot(P[:,:a+1].T, W[:,:a+1]))), Q[:,:a+1].T)
        #b0[a] = mnY - dot(mnX, B[a])
        
    
    # variance explained
    evx = 100.0*(1 - var_x/varX)
    evy = 100.0*(1 - var_y/varY)
    evz = 100.0*(1 - var_z/varZ)
    if scale=='loads':
        tnorm = apply_along_axis(vnorm, 0, T)
        T = T/tnorm
        W = W*tnorm
        Q = Q*tnorm
        knorm = apply_along_axis(vnorm, 0, K)
        L = L*knorm
        K = K/knorm
    
    return {'T':T, 'W':W, 'P':P, 'Q':Q, 'U':U, 'L':L, 'K':K, 'B':B, 'evx':evx, 'evy':evy, 'evz':evz,'mnx': mnX, 'mny': mnY, 'mnz': mnZ}    



def nipals_pls(X, Y, a_max, alpha=.7, ax_center=0, mode='normal', scale='scores', verbose=False):
    """Partial Least Sqaures Regression by the nipals algorithm.

    (X!Z)->Y
    :input:
        X : data matrix (m, n)
        Y : data matrix (m, l)

    :output:
      T : X-scores
      W : X-weights/Z-weights
      P : X-loadings
      Q : Y-loadings
      U : X-Y relation
      B : Regression coefficients X->Y
      b0: Regression coefficient intercept
      evx : X-explained variance
      evy : Y-explained variance
      evz : Z-explained variance

    :Notes:
    
    """
    if ax_center>=0:
        mn_x = expand_dims(X.mean(ax_center), ax_center)
        mn_y = expand_dims(Y.mean(ax_center), ax_center)
        X = X - mn_x
        Y = Y - mn_y

    varX = pow(X, 2).sum()
    varY = pow(Y, 2).sum()
    
    m, n = X.shape
    k, l = Y.shape

    # initialize 
    U = empty((k, a_max))
    Q = empty((l, a_max))
    T = empty((m, a_max))
    W = empty((n, a_max))
    P = empty((n, a_max))
    B = empty((a_max, n, l))
    b0 = empty((a_max, m, l))
    var_x = empty((a_max,))
    var_y = empty((a_max,))

    t1 = X[:,:1]
    for a in range(a_max):
        if verbose:
            print "\n Working on comp. %s" %a
        u = Y[:,:1]
        diff = 1
        MAX_ITER = 100
        lim = 1e-16
        niter = 0
        while (diff>lim and niter<MAX_ITER):
            niter += 1
            #u1 = u.copy()
            w = dot(X.T, u)
            w = w/sqrt(dot(w.T, w))
            #l = dot(Z, w)
            #k = dot(Z.T, l)
            #k = k/sqrt(dot(k.T, k))
            #w = alpha*k + (1-alpha)*w
            #w = w/sqrt(dot(w.T, w))
            t = dot(X, w)
            q = dot(Y.T, t)
            q = q/sqrt(dot(q.T, q))
            u = dot(Y, q)
            diff = vnorm(t1 - t)
            t1 = t.copy()
        if verbose:
            print "Converged after %s iterations" %niter
        #tt = dot(t.T, t)
        #p = dot(X.T, t)/tt
        #q = dot(Y.T, t)/tt
        #l = dot(Z, w)
        p = dot(X.T, t)/dot(t.T, t)
        p_norm = vnorm(p)
        t = t*p_norm
        w = w*p_norm
        p = p/p_norm
        
        U[:,a] = u.ravel()
        W[:,a] = w.ravel()
        P[:,a] = p.ravel()
        T[:,a] = t.ravel()
        Q[:,a] = q.ravel()

        X = X - dot(t, p.T)
        Y = Y - dot(t, q.T)

        var_x[a] = pow(X, 2).sum()
        var_y[a] = pow(Y, 2).sum()
    
        B[a] = dot(dot(W[:,:a+1], inv(dot(P[:,:a+1].T, W[:,:a+1]))), Q[:,:a+1].T)
        b0[a] =  mn_y - dot(mn_x, B[a])
    
    # variance explained
    evx = 100.0*(1 - var_x/varX)
    evy = 100.0*(1 - var_y/varY)
    
    if scale=='loads':
        tnorm = apply_along_axis(vnorm, 0, T)
        T = T/tnorm
        W = W*tnorm
        Q = Q*tnorm
    
    return {'T':T, 'W':W, 'P':P, 'Q':Q, 'U':U, 'B':B, 'b0':b0, 'evx':evx, 'evy':evy,
            'mnx': mnX, 'mny': mnY, 'xc': X, 'yc': Y}


########### Helper routines #########

def m_shape(array):
    return matrix(array).shape

def esvd(data, amax=None):
    """SVD with the option of economy sized calculation
    Calculate subspaces of X'X or XX' depending on the shape
    of the matrix.

    Good for extreme fat or thin matrices

    :notes:
    Numpy supports this by setting full_matrices=0
    """
    m, n = data.shape
    if m>=n:
        kernel = dot(data.T, data)
        if has_arpack:
            if amax==None:
                amax = n
                s, v = arpack.eigen_symmetric(kernel,k=amax, which='LM',
                                              maxiter=200,tol=1e-5)
        if has_sym:
            if amax==None:
                amax = n
                pcrange = None
            else:
                pcrange = [n-amax, n]
            s, v = symeig(kernel, range=pcrange, overwrite=True)
            s = s[::-1].real
            v = v[:,::-1].real
        else:
            u, s, vt = svd(kernel)
            v = vt.T
        s = sqrt(s)
        u = dot(data, v)/s
    else:
        kernel = dot(data, data.T)
        if has_sym:
            if amax==None:
                amax = m
                pcrange = None
            else:
                pcrange = [m-amax, m]
            s, u = symeig(kernel, range=pcrange, overwrite=True)
            s = s[::-1]
            u = u[:,::-1]
        else:
            u, s, vt = svd(kernel)
        s = sqrt(s)
        v = dot(data.T, u)/s
    # some use of symeig returns the 0 imaginary part
    return u.real, s.real, v.real

def vnorm(x):
    # assume column arrays (or vectors)
    return math.sqrt(dot(x.T, x))

def center(a, axis):
    # 0 = col center, 1 = row center, 2 = double center
    # -1 = nothing

    # check if we have a vector
    is_vec = len(a.shape)==1
    if not is_vec:
        is_vec = a.shape[0]==1 or a.shape[1]==1
    if is_vec:
        if axis==2:
            warnings.warn("Double centering of vecor ignored, using ordinary centering")
        if axis==-1:
            mn = 0
        else:
            mn = a.mean()
        return a - mn, mn
    # !!!fixme: use broadcasting
    if axis==-1:
        mn = zeros((1,a.shape[1],))
        #mn = tile(mn, (a.shape[0], 1))
    elif axis==0:
        mn = a.mean(0)[newaxis]
        #mn = tile(mn, (a.shape[0], 1)) 
    elif axis==1:
        mn = a.mean(1)[:,newaxis]
        #mn = tile(mn, (1, a.shape[1]))
    elif axis==2:
        mn = a.mean(0)[newaxis] + a.mean(1)[:,newaxis] - a.mean()
        return a - mn , a.mean(0)[newaxis]
    else:
        raise IOError("input error: axis must be in [-1,0,1,2]")

    return a - mn, mn

def scale(a, axis):
    if axis==-1:
        sc = zeros((a.shape[1],))
    elif axis==0:
        sc = a.std(0)
    elif axis==1:
        sc = a.std(1)[:,newaxis]
    else:
        raise IOError("input error: axis must be in [-1,0,1]")

    return a - sc, sc



## #PCA CALCS

## %  Calculate Q limit using unused eigenvalues
## temp = diag(s);
## if n < m
##   emod = temp(lv+1:n,:);
## else
##   emod = temp(lv+1:m,:);
## end
## th1 = sum(emod);
## th2 = sum(emod.^2);
## th3 = sum(emod.^3);
## h0 = 1 - ((2*th1*th3)/(3*th2^2));
## if h0 <= 0.0
## h0 = .0001;
## disp('  ')
## disp('Warning:  Distribution of unused eigenvalues indicates that')
## disp('          you should probably retain more PCs in the model.')
## end
## q = th1*(((1.65*sqrt(2*th2*h0^2)/th1) + 1 + th2*h0*(h0-1)/th1^2)^(1/h0));
## disp('  ')
## disp('The 95% Q limit is')
## disp(q)
## if plots >= 1
##   lim = [q q];
##   plot(scl,res,scllim,lim,'--b')
##   str = sprintf('Process Residual Q with 95 Percent Limit Based on %g PC Model',lv);
##   title(str)
##   xlabel('Sample Number')
##   ylabel('Residual')
##   pause
## end
## %  Calculate T^2 limit using ftest routine
## if lv > 1
##   if m > 300
##     tsq = (lv*(m-1)/(m-lv))*ftest(.95,300,lv,2);
##   else
##     tsq = (lv*(m-1)/(m-lv))*ftest(.95,m-lv,lv,2);
##   end
##   disp('  ')
##   disp('The 95% T^2 limit is')
##   disp(tsq)
## %  Calculate the value of T^2 by normalizing the scores to
## %  unit variance and summing them up
##   if plots >= 1.0
##     temp2 = scores*inv(diag(ssq(1:lv,2).^.5));
##     tsqvals = sum((temp2.^2)');
##     tlim = [tsq tsq];
##     plot(scl,tsqvals,scllim,tlim,'--b')
##     str = sprintf('Value of T^2 with 95 Percent Limit Based on %g PC Model',lv);
##     title(str)
##     xlabel('Sample Number')
##     ylabel('Value of T^2')
##   end
## else
##   disp('T^2 not calculated when number of latent variables = 1')
##   tsq = 1.96^2;
## end