00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #ifdef SEGMENTED
00024 #if defined(MPW) || defined(MW_CW)
00025 #pragma segment Rotfac
00026 #endif
00027 #endif
00028
00029 #include <math.h>
00030 #include "varimax.h"
00031
00032 #ifndef ROTMAXITER
00033 #define ROTMAXITER 100
00034 #endif
00035
00036 #ifndef ROTEPSILON
00037 #define ROTEPSILON 1e-5
00038 #endif
00039
00040
00041
00042
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
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
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 }
00109
00110 return (crit);
00111
00112 }
00113
00114
00115
00116
00117
00118
00119
00120 int varmx(float *aload,int nv, int nf, int method, float *params,
00121 float *fnorm,
00122 int itmax, float eps, int )
00123
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;
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
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 }
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 }
00183 }
00184 aloadk += nv;
00185 }
00186 aloadj += nv;
00187 }
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 }
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 }
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 }
00235 }
00236 aloadk += nv;
00237 }
00238 aloadj += nv;
00239 }
00240 }
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257 return (ict < itmax);
00258 }
00259
00260
00261
00262
00263
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
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 }
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 }
00327 return (reply);
00328 }
00329