EMAN2
varimax.cpp
Go to the documentation of this file.
00001 // this file is taken from MacAnova by Wei Zhang 1/4/2007. MacAnova is
00002 // also a GPL licensed software, so it is OK to use its source code here
00003 // if things change in the future, we will need re-implement this part
00004 
00005 
00006 /*
00007 *(C)* The file is part of the source distribution of MacAnova
00008 *(C)* version 4.12 or later
00009 *(C)*
00010 *(C)* Copyright (c) 2001 by Gary Oehlert and Christopher Bingham
00011 *(C)* unless indicated otherwise
00012 *(C)*
00013 *(C)* You may give out copies of this software;  for conditions see the
00014 *(C)* file COPYING included with this distribution
00015 *(C)*
00016 *(C)* This file is distributed WITHOUT ANY WARANTEE; without even
00017 *(C)* the implied warantee of MERCHANTABILITY or FITNESS FOR
00018 *(C)* A PARTICULAR PURPOSE
00019 *(C)*
00020 */
00021 
00022 
00023 #ifdef SEGMENTED
00024 #if defined(MPW) || defined(MW_CW)
00025 #pragma segment Rotfac
00026 #endif /*MPW||MW_CW*/
00027 #endif /*SEGMENTED*/
00028 
00029 #include <math.h>
00030 #include "varimax.h"
00031 
00032 #ifndef ROTMAXITER
00033 #define ROTMAXITER  100
00034 #endif /*ROTMAXITER*/
00035 
00036 #ifndef ROTEPSILON
00037 #define ROTEPSILON 1e-5
00038 #endif /*ROTEPSILON*/
00039 
00040 /*
00041   010614 reorganized code and added capability for quartimax and
00042          orthomax rotation
00043 */
00044 methodName     Methods[NMETHODS] =
00045 {
00046         {"Varimax",    3},
00047         {"Quartimax",  5},
00048         {"Equimax",    4},
00049         {"Orthomax",   5},
00050         {"Oblimin",    7}
00051 };
00052 
00053 /*
00054   Code acquired from Doug Hawkins, University of Minnesota, 12/93
00055 
00056   Comments from Fortran original
00057         Routine to do a varimax rotation on the real array aload(nv,nnf).
00058         If nnf is positive, the routine feels free to reverse and reorder
00059         the factors;this is suppressed in nnf is entered as the negative
00060         of its actual value.  This suppression is desirable when doing
00061         a q-mode analysis.
00062 
00063   Translated to C by C. Bingham with the following modifications
00064      Argument nnf is now nf and is assumed positive and additional argument
00065      fnorm replaces the local variable fnorm.  If fnorm == (float *) 0,
00066          reversing and reordering the factors is suppressed
00067 
00068          Also argument itmax has been added to control the maximum number of
00069          iterations and argument eps provides a convergence limit (originally
00070          hard wired as 1e-4).
00071          
00072          Information on the iteration is printed if verbose != 0
00073 
00074          varmx() returns a non-zero value if and only if fewer than itmax
00075          iterations are required.
00076 
00077          960919 Modified code to make it easier to add new rotation methods.
00078             except for slgihtly different error messages, it should have no
00079             effect on what it does.
00080      980303 Changed check before computation of rotaton angle to avoid
00081             atan2 domain error
00082 
00083      010614 added argument lambda to implement orthomax; lambda = 1 is
00084             varimax; lambda = 0 is quartimax
00085             Also computation of the criterion moved to separate function
00086             compcrit and names of certain variables were changed
00087 */
00088 
00089 float compcrit(float *loadings, int nv, int nf, float lambda)
00090 {
00091         float       crit = 0.0, *loadingj = loadings;
00092         float       fnv = (float) nv;
00093         int         i, j;
00094         
00095         for (j = 0; j < nf ;j++)
00096         {
00097                 float       s2 = 0.0;
00098 
00099                 for (i = 0;i < nv ;i++)
00100                 {
00101                         float        sq = loadingj[i]*loadingj[i];
00102 
00103                         s2 += sq;
00104                         crit += sq*sq;
00105                 }
00106                 crit -= lambda*s2*s2/fnv;
00107                 loadingj += nv;
00108         } /*for (j = 0;j < nf ;j++)*/
00109 
00110         return (crit);
00111         
00112 } /*compcrit()*/
00113 
00114 /*
00115   010614 added arguments method and param and modified code so that it
00116          finds optimal orthomax rotation with parameter lambda = params[0]
00117          lambda == 1 <==> variamx
00118          lambda == 2 <==> quartimax
00119 */
00120 int varmx(float *aload,int nv, int nf, int method, float *params,
00121                                   float *fnorm,
00122                                   int itmax, float eps, int )
00123 /*float aload[nv][1];*/
00124 {
00125         float         crit, startCrit, fnv = (float) nv;
00126         float        *aloadj, *aloadk;
00127         float         denominator, numerator, angl, trot;
00128         float         eps1 = eps, eps2 = eps;
00129         float         lambda = 0.0;     //avoid use lambda uninitialized
00130         int           inoim = 0, ict = 0, irot = 0;
00131         int           i, j, k, iflip, nf1 = nf - 1;
00132 
00133         if (method <= IORTHOMAX)
00134         {
00135                 lambda = params[0];
00136         }
00137 
00138         startCrit = crit = compcrit(aload, nv, nf, lambda);
00139 
00140         do /*while (inoim < 2 && ict < itmax && iflip);*/
00141         {
00142                 float      oldCrit = crit;
00143 
00144                 iflip = 0;
00145                 aloadj = aload;
00146 
00147                 for (j = 0;j < nf1 ;j++)
00148                 {
00149                         aloadk = aloadj + nv;
00150                         for (k = j + 1;k < nf ;k++)
00151                         {
00152                                 float      a = 0.0, b = 0.0, c = 0.0, d = 0.0, s = 0.0;
00153                                 
00154                                 for (i = 0;i < nv ;i++)
00155                                 {
00156                                         float    c2 = aloadj[i]*aloadj[i] - aloadk[i]*aloadk[i];
00157                                         float    s2 = 2.0f*aloadj[i]*aloadk[i];
00158 
00159                                         a += c2;
00160                                         b += s2;
00161                                         c += c2*c2 - s2*s2;
00162                                         d += c2*s2;
00163                                 } /*for (i = 0;i < nv ;i++)*/
00164 
00165                                 denominator = fnv*c + lambda*(b*b - a*a);
00166                                 numerator = 2.0f*(fnv*d - lambda*a*b);
00167 
00168                                 if (fabs(numerator) > eps1*fabs(denominator))
00169                                 {
00170                                         iflip = 1;
00171                                         irot++;
00172                                         angl = 0.25f*atan2(numerator,denominator);
00173                                         
00174                                         c = cos(angl);
00175                                         s = sin(angl);
00176                                         for (i = 0;i < nv ;i++)
00177                                         {
00178                                                 float   t = c*aloadj[i] + s*aloadk[i];
00179 
00180                                                 aloadk[i] = -s*aloadj[i] + c*aloadk[i];
00181                                                 aloadj[i] = t;
00182                                         } /*for (i = 0;i < nv ;i++)*/
00183                                 } /*if (fabs(numerator) >= eps1*fabs(denominator))*/
00184                                 aloadk += nv;
00185                         } /*for (k = j + 1;k < nf ;k++)*/
00186                         aloadj += nv;
00187                 } /*for (j = 0;j < nf1 ;j++)*/
00188                 ict++;
00189 
00190                 crit = compcrit(aload, nv, nf, lambda);
00191                 trot = (crit > 0.0f) ? (crit - oldCrit)/crit : 0.0f;
00192                 inoim++;
00193                 if (trot > eps2)
00194                 {
00195                         inoim = 0;
00196                 }
00197         } while (inoim < 2 && ict < itmax && iflip);
00198 
00199         if (fnorm != (float *) 0)
00200         {
00201                 aloadj = aload;
00202                 for (j = 0;j < nf ;j++)
00203                 {
00204                         float     ssj = 0, sj = 0;
00205                         
00206                         for (i = 0;i < nv ;i++)
00207                         {
00208                                 sj += aloadj[i];
00209                                 ssj += aloadj[i]*aloadj[i];
00210                         } /*for (i = 0;i < nv ;i++)*/
00211                         fnorm[j] = ssj;
00212                         if (sj <= 0.0)
00213                         {
00214                                 for (i = 0;i < nv ;i++)
00215                                 {
00216                                         aloadj[i] = -aloadj[i];
00217                                 }
00218                         } /*if (sj <= 0.0)*/
00219 
00220                         aloadk = aload;
00221                         for (k = 0;k < j ;k++)
00222                         {
00223                                 if (fnorm[k] < fnorm[j])
00224                                 {
00225                                         float       t = fnorm[k];
00226 
00227                                         fnorm[k] = fnorm[j];
00228                                         fnorm[j] = t;
00229                                         for (i = 0;i < nv ;i++)
00230                                         {
00231                                                 t = aloadj[i];
00232                                                 aloadj[i] = aloadk[i];
00233                                                 aloadk[i] = t;
00234                                         } /*for (i = 0;i < nv ;i++)*/
00235                                 } /*if (fnorm[k] < fnorm[j])*/
00236                                 aloadk += nv;
00237                         } /*for (k = 0;k < j ;k++)*/
00238                         aloadj += nv;
00239                 } /*for (j = 0;j < nf ;j++)*/
00240         } /*if (fnorm != (float *) 0)*/
00241 
00242         /*
00243         if (verbose)
00244         {
00245                 char   *outstr = OUTSTR;
00246                 
00247                 outstr += formatChar(outstr, Methods[method].name, CHARASIS);
00248                 outstr += formatChar(outstr, " starting criterion = ", CHARASIS);
00249                 outstr += formatDouble(outstr, startCrit, DODEFAULT | TRIMLEFT);
00250                 outstr += formatChar(outstr, ", final criterion = ", CHARASIS);
00251                 outstr += formatDouble(outstr, crit, DODEFAULT | TRIMLEFT);
00252                 putOUTSTR();
00253                 sprintf(OUTSTR, "%ld iterations and %ld rotations", ict, irot);
00254                 putOUTSTR();
00255         } if (verbose)*/
00256 
00257         return (ict < itmax);
00258 } /*varmx()*/
00259 
00260 /*
00261   011125 added argument knorm. 
00262          knorm != (float *) 0 signals Kaiser normalization with knorm
00263          providing scratch for row norms
00264 */
00265 int doRotation(int method, float *aload, int nv, int nf, 
00266                                            float *params, float *fnorm, float *knorm,
00267                                            int itmax, float eps, int verbose)
00268 {
00269         int        reply = -1;
00270 
00271         if (method <= IORTHOMAX)
00272         {
00273                 if (method < IORTHOMAX)
00274                 {
00275                         if (method == IVARIMAX)
00276                         {
00277                                 params[0] = 1.0;
00278                         }
00279                         else if (method == IQUARTIMAX)
00280                         {
00281                                 params[0] = 0.0;
00282                         }
00283                         else
00284                         {
00285                                 /* equimax */
00286                                 params[0] = .5f*(float) nf;
00287                         }
00288                 }
00289                 if (knorm != (float *) 0)
00290                 {
00291                         int         i, j, k;
00292                         float       s;
00293 
00294                         for (i = 0; i < nv; i++)
00295                         {
00296                                 k = i;
00297                                 s = 0.0;
00298                                 for (j = 0; j < nf; j++, k += nv)
00299                                 {
00300                                         s += aload[k]*aload[k];
00301                                 }
00302                                 knorm[i] = s = (s > 0) ? sqrt(s) : 1.0f;
00303 
00304                                 k = i;
00305                                 for (j = 0; j < nf; j++, k += nv)
00306                                 {
00307                                         aload[k] /= s;
00308                                 }
00309                         }
00310                 } /*if (fcopy != (float *) 0)*/
00311                 
00312                 reply = varmx(aload, nv, nf, method, params, fnorm, itmax,
00313                                           eps, verbose);
00314                 if (knorm != (float *) 0)
00315                 {
00316                         int      i, j, k = 0;
00317                         
00318                         for (j = 0; j < nf; j++)
00319                         {
00320                                 for (i = 0; i < nv; i++, k++)
00321                                 {
00322                                         aload[k] *= knorm[i];
00323                                 }
00324                         }
00325                 }
00326         } /*if (method <= IORTHOMAX)*/
00327         return (reply);
00328 } /*doRotation()*/
00329