ENH: hack to estimate variance explained

antsSCCANObject.hxx

@@ -764,7 +764,9 @@ TRealType antsSCCANObject<TInputImage, TRealType>
   unsigned long mind=this->m_MatrixP.rows();
   if ( mind > this->m_MatrixP.cols() ) mind=this->m_MatrixP.cols();
   double avgdifffromevec=0;
-  double evalsum=0;
+  double eval_sum=0;
+  double trace_sum=0;
+  double sub_trace=0;
   // we estimate variance explained by  \sum_i eigenvalue_i / trace(A) ...
   double vex=0;  
   for ( unsigned int i=0; i < n_vecs ; i++ ) 
@@ -773,15 +775,19 @@ TRealType antsSCCANObject<TInputImage, TRealType>
     vnl_diag_matrix<TRealType> indicator(this->m_MatrixP.cols(),1);
     for ( unsigned int j=0; j< u.size(); j++) if ( fabs(u(j)) < this->m_Epsilon ) indicator(j,j)=0; 
     MatrixType pmod=this->m_MatrixP*indicator; 
+    if ( sub_trace <= 0 ) {
     MatrixType sub_covmat=pmod*pmod.transpose(); 
-    double sub_trace=vnl_trace<double>(sub_covmat); 
+    sub_trace=vnl_trace<double>(sub_covmat);
+    } 
     VectorType proj=(pmod*u);
     VectorType m=pmod.transpose()*proj;
     double eigenvalue_i=1.e-5;
     double unorm=u.two_norm();
     if (unorm > this->m_Epsilon) eigenvalue_i=m.two_norm()/unorm;
-    if ( i < mind-1 ) evalsum+=eigenvalue_i;
-    vex+=eigenvalue_i/sub_trace;
+    if ( i < mind-1 ) eval_sum+=eigenvalue_i;
+    trace_sum+=sub_trace;
+    vex+=(eigenvalue_i/sub_trace)/(double)n_vecs;
+    //    std::cout << eigenvalue_i << " vs trace " << sub_trace << " trace " << trace <<  " vex " << eigenvalue_i/sub_trace << std::endl;
     VectorType diff=u*eigenvalue_i-m;
     /** this is a good metric of the difference from the eigenvalue 
      *  because, over iterations, we minimize \|  u*eigenvalue_i-m \|^2  
@@ -791,7 +797,6 @@ TRealType antsSCCANObject<TInputImage, TRealType>
     avgdifffromevec+=d;
     this->m_CanonicalCorrelations[i]=eigenvalue_i;
   }
-  if ( vex > 1 ) vex=0;  
   std::cout <<" variance explained: " << vex << std::endl;
   //  std::cout <<" eval diff " <<   avgdifffromevec/n_vecs << std::endl;
   return avgdifffromevec/(TRealType)n_vecs; 

scripts/sccan_tests

@@ -22,8 +22,8 @@ params=" -e 0 -n 2 -i 20 -r 0 -p 0  " # most common parameters
 t1=`${PATHTOSCCAN}/sccan --scca two-view[ ${PATHTODATA}myview1.csv , ${PATHTODATA}myview_mismatch.csv , na , na , 0.5 , 0.5 ] -o ${OUT}.nii.gz $params `
 test1=$?
 t2=`${PATHTOSCCAN}/sccan --sparse-svd [ ${PATHTODATA}myview3.csv , na , -0.15 ] -o ${OUT}.nii.gz -n 15 -i 500 `
-#${PATHTOSCCAN}/sccan --sparse-svd [ ${PATHTODATA}myview2.csv , na , 0.15 ] -o ${OUT}.nii.gz -n 28 -i 111
-#exit
+${PATHTOSCCAN}/sccan --sparse-svd [ ${PATHTODATA}myview1.csv , na , -.05 ] -o ${OUT}.nii.gz -n 28 -i 111
+exit
 test2=$?
 paramsb=" -e 0 -n 5 -i 50 -r 0 -p 0  --partial-scca-option PQminusR " # most common parameters 
 t3=`${PATHTOSCCAN}/sccan --scca partial[ ${PATHTODATA}myview2.csv , ${PATHTODATA}myview3.csv , ${PATHTODATA}myview4.csv , na , na , na , -0.4 , -0.9 , 1] -o ${OUT}.nii.gz $paramsb `