From 8e607c7b1a7ea353a04c448dd7cc8fd28844cfa8 Mon Sep 17 00:00:00 2001
From: flatberg <flatberg@pvv.ntnu.no>
Date: Wed, 31 Jan 2007 11:57:04 +0000
Subject: [PATCH] bugfixed scaling issues

---
 fluents/lib/cx_stats.py | 36 ++++++++++++++++--------------------
 1 file changed, 16 insertions(+), 20 deletions(-)

diff --git a/fluents/lib/cx_stats.py b/fluents/lib/cx_stats.py
index 31d0fe4..078df5f 100644
--- a/fluents/lib/cx_stats.py
+++ b/fluents/lib/cx_stats.py
@@ -9,7 +9,8 @@ from select_generators import shuffle_1d
 from engines import *
 import time
 
-def hotelling(P, Pfull, p_center='med', cov_center='med',
+
+def hotelling(Pcv, P, p_center='med', cov_center='med',
               alpha=0.3, crot=True, strict=False, metric=None):
     """Returns regularized hotelling T^2.
     
@@ -20,35 +21,35 @@ def hotelling(P, Pfull, p_center='med', cov_center='med',
     alpha -- regularisation
     crot -- rotate submodels toward full?
     strict -- only rotate 90 degree ?
-    metric -- inverse metric matrix (if P and Pfull from metric pca/pls)
+    metric -- inverse metric matrix (if Pcv and P from metric pca/pls)
     
     """
-    m, n = Pfull.shape
+    m, n = P.shape
     if metric==None:
         metric = eye(m, dtype='<f8')
-    Pfull = dot(metric.T, asarray(Pfull))
-    n_sets,n,amax = P.shape
+    P = dot(metric.T, asarray(P))
+    n_sets, n, amax = Pcv.shape
     # allocate
     T_sq = empty((n, ),dtype='f')
     Cov_i = zeros((n, amax, amax),dtype='f')
     
     # rotate sub_models to full model
     if crot:
-        for i,Pi in enumerate(P):
+        for i, Pi in enumerate(Pcv):
             Pi = dot(metric.T, Pi)
-            P[i] = procrustes(Pfull, Pi, strict=strict)
+            Pcv[i] = procrustes(P, Pi, strict=strict)
 
     # center of pnull
     if p_center=='med':
-        P_ctr = median(P, 0)
+        P_ctr = median(Pcv, 0)
     elif p_center=='mean':
         # fixme: mean is unstable
-        P_ctr = mean(P, 0)
+        P_ctr = mean(Pcv, 0)
     else: #use full
-        P_ctr = Pfull
+        P_ctr = P
 
     for i in xrange(n):
-        Pi = P[:,i,:] # (n_sets x amax) 
+        Pi = Pcv[:,i,:] # (n_sets x amax) 
         Pi_ctr = P_ctr[i,:] # (1 x amax)
         Pim = (Pi - Pi_ctr[newaxis])*sqrt(n_sets-1)
         Cov_i[i] = (1./n_sets)*dot(Pim.T, Pim)
@@ -131,8 +132,6 @@ def pls_qvals(a, b, aopt=None, alpha=.3,
     crot -- bool, use rotations of sub models?
     strict -- bool, use stict (rot/flips only) rotations?
     metric -- bool, use row metric?
-
-    
     """
     
     m, n = a.shape
@@ -144,8 +143,8 @@ def pls_qvals(a, b, aopt=None, alpha=.3,
         dat = bridge(a, b, aopt, 'loads', 'fast')
     else:
         dat = pls(a, b, aopt, 'loads', 'fast')
-    W = pls_jkW(a, b, aopt, n_blocks=None, algo=algo)
-    tsq_full = hotelling(W, dat['W'], p_center=p_center,
+    Wcv = pls_jkW(a, b, aopt, n_blocks=None, algo=algo)
+    tsq_full = hotelling(Wcv, dat['W'], p_center=p_center,
                          alpha=alpha, crot=crot, strict=strict,
                          cov_center=cov_center, metric=metric)
     t0 = time.time()
@@ -156,8 +155,8 @@ def pls_qvals(a, b, aopt=None, alpha=.3,
             dat = bridge(a, b_shuff, aopt, 'loads','fast')
         else:
             dat = pls(a, b, aopt, 'loads', 'fast')
-        W = pls_jkW(a, b_shuff, aopt, n_blocks=None, algo=algo)
-        TSQ[:,i] = hotelling(W, dat['W'],p_center=p_center,
+        Wcv = pls_jkW(a, b_shuff, aopt, n_blocks=None, algo=algo)
+        TSQ[:,i] = hotelling(Wcv, dat['W'], p_center=p_center,
                              alpha=alpha, crot=crot, strict=strict,
                              cov_center=cov_center, metric=metric)
         print time.time() - t1
@@ -194,8 +193,6 @@ def ensure_strict(C, only_flips=True):
         raise ValueError,"Implement this!"
     return Cm*S
 
-
-
 def leverage(aopt=1,*args):
     """Returns leverages
     input : aopt, number of components to base leverage calculations on
@@ -213,7 +210,6 @@ def leverage(aopt=1,*args):
         lev.append(lev_u)
     return lev
 
-
 def variances(a,t,p):
     """Returns explained variance and ind. var from blm-params.
     input: