Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

varimax.cpp File Reference

#include <math.h>
#include "varimax.h"

Include dependency graph for varimax.cpp:

Include dependency graph

Go to the source code of this file.

Defines

#define ROTMAXITER   100
#define ROTEPSILON   1e-5

Functions

float compcrit (float *loadings, int nv, int nf, float lambda)
int varmx (float *aload, int nv, int nf, int method, float *params, float *fnorm, int itmax, float eps, int)
int doRotation (int method, float *aload, int nv, int nf, float *params, float *fnorm, float *knorm, int itmax, float eps, int verbose)

Variables

methodName Methods [NMETHODS]


Define Documentation

#define ROTEPSILON   1e-5
 

Definition at line 37 of file varimax.cpp.

#define ROTMAXITER   100
 

Definition at line 33 of file varimax.cpp.


Function Documentation

float compcrit float *  loadings,
int  nv,
int  nf,
float  lambda
 

Definition at line 89 of file varimax.cpp.

Referenced by varmx().

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()*/

int doRotation int  method,
float *  aload,
int  nv,
int  nf,
float *  params,
float *  fnorm,
float *  knorm,
int  itmax,
float  eps,
int  verbose
 

Definition at line 265 of file varimax.cpp.

References sqrt(), and varmx().

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()*/

int varmx float *  aload,
int  nv,
int  nf,
int  method,
float *  params,
float *  fnorm,
int  itmax,
float  eps,
int 
 

Definition at line 120 of file varimax.cpp.

References b, compcrit(), and t.

Referenced by EMAN::varimax::analyze(), and doRotation().

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()*/


Variable Documentation

methodName Methods[NMETHODS]
 

Initial value:

{
        {"Varimax",    3},
        {"Quartimax",  5},
        {"Equimax",    4},
        {"Orthomax",   5},
        {"Oblimin",    7}
}

Definition at line 44 of file varimax.cpp.


Generated on Tue Jun 11 13:47:34 2013 for EMAN2 by  doxygen 1.3.9.1