diff --git a/antsSCCANObject.hxx b/antsSCCANObject.hxx index 7eaab10..e5cb793 100644 --- a/antsSCCANObject.hxx +++ b/antsSCCANObject.hxx @@ -764,7 +764,9 @@ TRealType antsSCCANObject 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 vnl_diag_matrix 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(sub_covmat); + sub_trace=vnl_trace(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 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; diff --git a/scripts/sccan_tests b/scripts/sccan_tests index 061742f..5077a20 100755 --- a/scripts/sccan_tests +++ b/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 `