| 31 | | package org.libspark.flartoolkit.core { |
|---|
| 32 | | import org.libspark.flartoolkit.FLARException; |
|---|
| 33 | | import org.libspark.flartoolkit.utils.ArrayUtil; |
|---|
| 34 | | |
|---|
| 35 | | /** |
|---|
| 36 | | * ARMat構造体に対応するクラス typedef struct { double *m; int row; int clm; }ARMat; |
|---|
| 37 | | * |
|---|
| 38 | | */ |
|---|
| 39 | | public class FLARMat { |
|---|
| 40 | | |
|---|
| 41 | | /** |
|---|
| 42 | | * 配列サイズと行列サイズは必ずしも一致しないことに注意 返された配列のサイズを行列の大きさとして使わないこと！ |
|---|
| 43 | | * |
|---|
| 44 | | */ |
|---|
| 45 | | protected var m:Array; |
|---|
| 46 | | // double[][] |
|---|
| 47 | | |
|---|
| 48 | | private var clm:int; |
|---|
| 49 | | private var row:int; |
|---|
| 50 | | |
|---|
| 51 | | /** |
|---|
| 52 | | * デフォルトコンストラクタは機能しません。 |
|---|
| 53 | | * |
|---|
| 54 | | * @throws FLARException |
|---|
| 55 | | */ |
|---|
| 56 | | // protected function FLARMat() |
|---|
| 57 | | // { |
|---|
| 58 | | // throw new FLARException(); |
|---|
| 59 | | // } |
|---|
| 60 | | |
|---|
| 61 | | public function FLARMat(i_row:int, i_clm:int) { |
|---|
| 62 | | // m = new double[i_row][i_clm]; |
|---|
| 63 | | m = ArrayUtil.createJaggedArray(i_row, i_clm); |
|---|
| 64 | | clm = i_clm; |
|---|
| 65 | | row = i_row; |
|---|
| 66 | | } |
|---|
| 67 | | |
|---|
| 68 | | /** |
|---|
| 69 | | * i_row x i_clmサイズの行列を格納できるように行列サイズを変更します。 実行後、行列の各値は不定になります。 |
|---|
| 70 | | * |
|---|
| 71 | | * @param i_row |
|---|
| 72 | | * @param i_clm |
|---|
| 73 | | */ |
|---|
| 74 | | public function realloc(i_row:int, i_clm:int):void { |
|---|
| 75 | | if (i_row <= this.m.length && i_clm <= this.m[0].length) { |
|---|
| 76 | | // 十分な配列があれば何もしない。 |
|---|
| 77 | | } else { |
|---|
| 78 | | // 不十分なら取り直す。 |
|---|
| 79 | | // m = new double[i_row][i_clm]; |
|---|
| 80 | | m = ArrayUtil.createJaggedArray(i_row, i_clm); |
|---|
| 81 | | } |
|---|
| 82 | | this.clm = i_clm; |
|---|
| 83 | | this.row = i_row; |
|---|
| 84 | | } |
|---|
| 85 | | |
|---|
| 86 | | public function getClm():int { |
|---|
| 87 | | return clm; |
|---|
| 88 | | } |
|---|
| 89 | | |
|---|
| 90 | | public function getRow():int { |
|---|
| 91 | | return row; |
|---|
| 92 | | } |
|---|
| 93 | | |
|---|
| 94 | | /** |
|---|
| 95 | | * 行列をゼロクリアする。 |
|---|
| 96 | | */ |
|---|
| 97 | | public function zeroClear():void { |
|---|
| 98 | | var i:int; |
|---|
| 99 | | var i2:int; |
|---|
| 100 | | // For順変更OK |
|---|
| 101 | | for (i = row - 1;i >= 0; i--) { |
|---|
| 102 | | for (i2 = clm - 1;i2 >= 0; i2--) { |
|---|
| 103 | | m[i][i2] = 0.0; |
|---|
| 104 | | } |
|---|
| 105 | | } |
|---|
| 106 | | } |
|---|
| 107 | | |
|---|
| 108 | | /** |
|---|
| 109 | | * i_copy_fromの内容を自分自身にコピーします。 高さ・幅は同一で無いと失敗します。 |
|---|
| 110 | | * |
|---|
| 111 | | * @param i_copy_from |
|---|
| 112 | | */ |
|---|
| 113 | | public function copyFrom(i_copy_from:FLARMat):void { |
|---|
| 114 | | // サイズ確認 |
|---|
| 115 | | if (this.row != i_copy_from.row || this.clm != i_copy_from.clm) { |
|---|
| 116 | | throw new FLARException(); |
|---|
| 117 | | } |
|---|
| 118 | | // 値コピー |
|---|
| 119 | | var r:int; |
|---|
| 120 | | var c:int; |
|---|
| 121 | | for (r = this.row - 1;r >= 0; r--) { |
|---|
| 122 | | for (c = this.clm - 1;c >= 0; c--) { |
|---|
| 123 | | this.m[r][c] = i_copy_from.m[r][c]; |
|---|
| 124 | | } |
|---|
| 125 | | } |
|---|
| 126 | | } |
|---|
| 127 | | |
|---|
| 128 | | public function getArray():Array { |
|---|
| 129 | | return m; |
|---|
| 130 | | } |
|---|
| 131 | | |
|---|
| 132 | | // public void getRowVec(int i_row,FLARVec o_vec) |
|---|
| 133 | | // { |
|---|
| 134 | | // o_vec.set(this.m[i_row],this.clm); |
|---|
| 135 | | // } |
|---|
| 136 | | /** |
|---|
| 137 | | * aとbの積を自分自身に格納する。arMatrixMul()の代替品 |
|---|
| 138 | | * |
|---|
| 139 | | * @param a |
|---|
| 140 | | * @param b |
|---|
| 141 | | * @throws FLARException |
|---|
| 142 | | */ |
|---|
| 143 | | public function matrixMul(a:FLARMat, b:FLARMat):void { |
|---|
| 144 | | if (a.clm != b.row || this.row != a.row || this.clm != b.clm) { |
|---|
| 145 | | throw new FLARException(); |
|---|
| 146 | | } |
|---|
| 147 | | var w:Number; |
|---|
| 148 | | var r:int; |
|---|
| 149 | | var c:int; |
|---|
| 150 | | var i:int; |
|---|
| 151 | | var am:Array = a.m; |
|---|
| 152 | | var bm:Array = b.m; |
|---|
| 153 | | var dm:Array = this.m; |
|---|
| 154 | | // For順変更禁止 |
|---|
| 155 | | for (r = 0;r < this.row; r++) { |
|---|
| 156 | | for (c = 0;c < this.clm; c++) { |
|---|
| 157 | | w = 0.0; |
|---|
| 158 | | // dest.setARELEM0(r, c,0.0); |
|---|
| 159 | | for (i = 0;i < a.clm; i++) { |
|---|
| 160 | | w += am[r][i] * bm[i][c];// ARELEM0(dest, r, c) +=ARELEM0(a, r, i) * ARELEM0(b,i, c); |
|---|
| 161 | | } |
|---|
| 162 | | dm[r][c] = w; |
|---|
| 163 | | } |
|---|
| 164 | | } |
|---|
| 165 | | } |
|---|
| 166 | | |
|---|
| 167 | | private var wk_nos_matrixSelfInv:Array = new Array(50); |
|---|
| 168 | | |
|---|
| 169 | | // new int[50]; |
|---|
| 170 | | |
|---|
| 171 | | // private final static double matrixSelfInv_epsl=1.0e-10; |
|---|
| 172 | | /** |
|---|
| 173 | | * i_targetを逆行列に変換する。arMatrixSelfInv()と、arMatrixSelfInv_minv()関数を合成してあります。 |
|---|
| 174 | | * OPTIMIZE STEP[485->422] |
|---|
| 175 | | * |
|---|
| 176 | | * @param i_target |
|---|
| 177 | | * 逆行列にする行列 |
|---|
| 178 | | * @return 逆行列があればTRUE/無ければFALSE |
|---|
| 179 | | * |
|---|
| 180 | | * @throws FLARException |
|---|
| 181 | | */ |
|---|
| 182 | | public function matrixSelfInv():Boolean { |
|---|
| 183 | | var ap:Array = this.m; |
|---|
| 184 | | var dimen:int = this.row; |
|---|
| 185 | | var dimen_1:int = dimen - 1; |
|---|
| 186 | | var ap_n:Array; |
|---|
| 187 | | var ap_ip:Array; |
|---|
| 188 | | var ap_i:Array; |
|---|
| 189 | | var j:int; |
|---|
| 190 | | var ip:int; |
|---|
| 191 | | var nwork:int; |
|---|
| 192 | | var nos:Array = wk_nos_matrixSelfInv; |
|---|
| 193 | | // この関数で初期化される。 |
|---|
| 194 | | // double epsl; |
|---|
| 195 | | var p:Number; |
|---|
| 196 | | var pbuf:Number; |
|---|
| 197 | | var work:Number; |
|---|
| 198 | | |
|---|
| 199 | | /* check size */ |
|---|
| 200 | | switch (dimen) { |
|---|
| 201 | | case 0: |
|---|
| 202 | | throw new FLARException(); |
|---|
| 203 | | case 1: |
|---|
| 204 | | ap[0][0] = 1.0 / ap[0][0]; |
|---|
| 205 | | // *ap = 1.0 / (*ap); |
|---|
| 206 | | return true;/* 1 dimension */ |
|---|
| 207 | | } |
|---|
| 208 | | |
|---|
| 209 | | for (var n:int = 0;n < dimen; n++) { |
|---|
| 210 | | nos[n] = n; |
|---|
| 211 | | } |
|---|
| 212 | | |
|---|
| 213 | | /* |
|---|
| 214 | | * nyatla memo ipが定まらないで計算が行われる場合があるので挿入。 ループ内で0初期化していいかが判らない。 |
|---|
| 215 | | */ |
|---|
| 216 | | ip = 0; |
|---|
| 217 | | // For順変更禁止 |
|---|
| 218 | | var i:int; |
|---|
| 219 | | for (n = 0;n < dimen; n++) { |
|---|
| 220 | | ap_n = ap[n]; |
|---|
| 221 | | // wcp = ap + n * rowa; |
|---|
| 222 | | p = 0.0; |
|---|
| 223 | | for (i = n;i < dimen; i++) { |
|---|
| 224 | | // for(i = n, wap = wcp, p = |
|---|
| 225 | | // 0.0; i < dimen ; i++, wap += |
|---|
| 226 | | // rowa) |
|---|
| 227 | | if (p < (pbuf = Math.abs(ap[i][0]))) { |
|---|
| 228 | | p = pbuf; |
|---|
| 229 | | ip = i; |
|---|
| 230 | | } |
|---|
| 231 | | } |
|---|
| 232 | | // if (p <= matrixSelfInv_epsl){ |
|---|
| 233 | | if (p == 0.0) { |
|---|
| 234 | | return false; |
|---|
| 235 | | // throw new FLARException(); |
|---|
| 236 | | } |
|---|
| 237 | | |
|---|
| 238 | | nwork = nos[ip]; |
|---|
| 239 | | nos[ip] = nos[n]; |
|---|
| 240 | | nos[n] = nwork; |
|---|
| 241 | | |
|---|
| 242 | | ap_ip = ap[ip]; |
|---|
| 243 | | for (j = 0;j < dimen; j++) { |
|---|
| 244 | | // for(j = 0, wap = ap + ip * rowa, |
|---|
| 245 | | // wbp = wcp; j < dimen ; j++) { |
|---|
| 246 | | work = ap_ip[j]; |
|---|
| 247 | | // work = *wap; |
|---|
| 248 | | ap_ip[j] = ap_n[j]; |
|---|
| 249 | | ap_n[j] = work; |
|---|
| 250 | | } |
|---|
| 251 | | |
|---|
| 252 | | work = ap_n[0]; |
|---|
| 253 | | for (j = 0;j < dimen_1; j++) { |
|---|
| 254 | | // for(j = 1, wap = wcp, work = |
|---|
| 255 | | // *wcp; j < dimen ; j++, wap++) |
|---|
| 256 | | ap_n[j] = ap_n[j + 1] / work;// *wap = *(wap + 1) / work; |
|---|
| 257 | | } |
|---|
| 258 | | ap_n[j] = 1.0 / work; |
|---|
| 259 | | // *wap = 1.0 / work; |
|---|
| 260 | | for (i = 0;i < dimen; i++) { |
|---|
| 261 | | if (i != n) { |
|---|
| 262 | | ap_i = ap[i]; |
|---|
| 263 | | // wap = ap + i * rowa; |
|---|
| 264 | | |
|---|
| 265 | | work = ap_i[0]; |
|---|
| 266 | | for (j = 0; j < dimen_1; j++) { |
|---|
| 267 | | // for(j = 1, wbp = wcp,work = *wap;j < dimen ;j++, wap++, wbp++) |
|---|
| 268 | | ap_i[j] = ap_i[j + 1] - work * ap_n[j];// wap = *(wap +1) - work *(*wbp); |
|---|
| 269 | | } |
|---|
| 270 | | ap_i[j] = -work * ap_n[j];// *wap = -work * (*wbp); |
|---|
| 271 | | } |
|---|
| 272 | | } |
|---|
| 273 | | } |
|---|
| 274 | | |
|---|
| 275 | | for (n = 0;n < dimen; n++) { |
|---|
| 276 | | for (j = n;j < dimen; j++) { |
|---|
| 277 | | if (nos[j] == n) { |
|---|
| 278 | | break; |
|---|
| 279 | | } |
|---|
| 280 | | } |
|---|
| 281 | | nos[j] = nos[n]; |
|---|
| 282 | | for (i = 0; i < dimen; i++) { |
|---|
| 283 | | // for(i = 0, wap = ap + j, wbp |
|---|
| 284 | | // = ap + n; i < dimen ;i++, wap |
|---|
| 285 | | // += rowa, wbp += rowa) { |
|---|
| 286 | | ap_i = ap[i]; |
|---|
| 287 | | work = ap_i[j]; |
|---|
| 288 | | // work = *wap; |
|---|
| 289 | | ap_i[j] = ap_i[n]; |
|---|
| 290 | | // *wap = *wbp; |
|---|
| 291 | | ap_i[n] = work;// *wbp = work; |
|---|
| 292 | | } |
|---|
| 293 | | } |
|---|
| 294 | | return true; |
|---|
| 295 | | } |
|---|
| 296 | | |
|---|
| 297 | | /** |
|---|
| 298 | | * sourceの転置行列をdestに得る。arMatrixTrans()の代替品 |
|---|
| 299 | | * |
|---|
| 300 | | * @param dest |
|---|
| 301 | | * @param source |
|---|
| 302 | | * @return |
|---|
| 303 | | */ |
|---|
| 304 | | public static function matrixTrans(dest:FLARMat, source:FLARMat):void { |
|---|
| 305 | | if (dest.row != source.clm || dest.clm != source.row) { |
|---|
| 306 | | throw new FLARException(); |
|---|
| 307 | | } |
|---|
| 308 | | FLARException.trap("未チェックのパス"); |
|---|
| 309 | | // For順変更禁止 |
|---|
| 310 | | var r:int; |
|---|
| 311 | | var c:int; |
|---|
| 312 | | for (r = 0; r < dest.row; r++) { |
|---|
| 313 | | for (c = 0; c < dest.clm; c++) { |
|---|
| 314 | | dest.m[r][c] = source.m[c][r]; |
|---|
| 315 | | } |
|---|
| 316 | | } |
|---|
| 317 | | } |
|---|
| 318 | | |
|---|
| 319 | | /** |
|---|
| 320 | | * unitを単位行列に初期化する。arMatrixUnitの代替品 |
|---|
| 321 | | * |
|---|
| 322 | | * @param unit |
|---|
| 323 | | */ |
|---|
| 324 | | public static function matrixUnit(unit:FLARMat):void { |
|---|
| 325 | | if (unit.row != unit.clm) { |
|---|
| 326 | | throw new FLARException(); |
|---|
| 327 | | } |
|---|
| 328 | | FLARException.trap("未チェックのパス"); |
|---|
| 329 | | // For順変更禁止 |
|---|
| 330 | | var r:int; |
|---|
| 331 | | var c:int; |
|---|
| 332 | | for (r = 0; r < unit.getRow(); r++) { |
|---|
| 333 | | for (c = 0; c < unit.getClm(); c++) { |
|---|
| 334 | | if (r == c) { |
|---|
| 335 | | unit.m[r][c] = 1.0; |
|---|
| 336 | | } else { |
|---|
| 337 | | unit.m[r][c] = 0.0; |
|---|
| 338 | | } |
|---|
| 339 | | } |
|---|
| 340 | | } |
|---|
| 341 | | } |
|---|
| 342 | | |
|---|
| 343 | | /** |
|---|
| 344 | | * sourceの内容を自身に複製する。 Optimized 2008.04.19 |
|---|
| 345 | | * |
|---|
| 346 | | * @param i_source |
|---|
| 347 | | * @return |
|---|
| 348 | | */ |
|---|
| 349 | | public function matrixDup(i_source:FLARMat):void { |
|---|
| 350 | | // 自身の配列サイズを相手のそれより大きいことを保障する。 |
|---|
| 351 | | this.realloc(i_source.row, i_source.clm); |
|---|
| 352 | | // 内容を転写 |
|---|
| 353 | | var r:int; |
|---|
| 354 | | var c:int; |
|---|
| 355 | | var src_m:Array; |
|---|
| 356 | | var dest_m:Array; |
|---|
| 357 | | src_m = i_source.m; |
|---|
| 358 | | dest_m = this.m; |
|---|
| 359 | | // コピーはFor順を変えてもOK |
|---|
| 360 | | for (r = this.row - 1; r >= 0; r--) { |
|---|
| 361 | | for (c = this.clm - 1; c >= 0; c--) { |
|---|
| 362 | | dest_m[r][c] = src_m[r][c]; |
|---|
| 363 | | } |
|---|
| 364 | | } |
|---|
| 365 | | } |
|---|
| 366 | | |
|---|
| 367 | | public function matrixAllocDup():FLARMat { |
|---|
| 368 | | var result:FLARMat = new FLARMat(this.row, this.clm); |
|---|
| 369 | | // コピー |
|---|
| 370 | | var r:int; |
|---|
| 371 | | var c:int; |
|---|
| 372 | | var dest_m:Array; |
|---|
| 373 | | var src_m:Array; |
|---|
| 374 | | dest_m = result.m; |
|---|
| 375 | | src_m = this.m; |
|---|
| 376 | | // コピーはFor順を変えてもOK |
|---|
| 377 | | for (r = this.row - 1; r >= 0; r--) { |
|---|
| 378 | | for (c = this.clm - 1; c >= 0; c--) { |
|---|
| 379 | | dest_m[r][c] = src_m[r][c]; |
|---|
| 380 | | } |
|---|
| 381 | | } |
|---|
| 382 | | return result; |
|---|
| 383 | | } |
|---|
| 384 | | |
|---|
| 385 | | /** |
|---|
| 386 | | * arMatrixInv関数の代替品です。 destにsourceの逆行列を返します。 |
|---|
| 387 | | * |
|---|
| 388 | | * @param dest |
|---|
| 389 | | * @param source |
|---|
| 390 | | * @throws FLARException |
|---|
| 391 | | */ |
|---|
| 392 | | public static function matrixInv(dest:FLARMat, source:FLARMat):void { |
|---|
| 393 | | FLARException.trap("未チェックのパス"); |
|---|
| 394 | | dest.matrixDup(source); |
|---|
| 395 | | |
|---|
| 396 | | FLARException.trap("未チェックのパス"); |
|---|
| 397 | | dest.matrixSelfInv(); |
|---|
| 398 | | } |
|---|
| 399 | | |
|---|
| 400 | | public function matrixAllocInv():FLARMat { |
|---|
| 401 | | FLARException.trap("未チェックのパス"); |
|---|
| 402 | | var result:FLARMat = matrixAllocDup(); |
|---|
| 403 | | |
|---|
| 404 | | FLARException.trap("未チェックのパス"); |
|---|
| 405 | | result.matrixSelfInv(); |
|---|
| 406 | | return result; |
|---|
| 407 | | } |
|---|
| 408 | | |
|---|
| 409 | | /** |
|---|
| 410 | | * dim x dim の単位行列を作る。 |
|---|
| 411 | | * |
|---|
| 412 | | * @param dim |
|---|
| 413 | | * @return |
|---|
| 414 | | * @throws FLARException |
|---|
| 415 | | */ |
|---|
| 416 | | public static function matrixAllocUnit(dim:int):FLARMat { |
|---|
| 417 | | FLARException.trap("未チェックのパス"); |
|---|
| 418 | | var result:FLARMat = new FLARMat(dim, dim); |
|---|
| 419 | | FLARException.trap("未チェックのパス"); |
|---|
| 420 | | FLARMat.matrixUnit(result); |
|---|
| 421 | | return result; |
|---|
| 422 | | } |
|---|
| 423 | | |
|---|
| 424 | | /** |
|---|
| 425 | | * arMatrixDispの代替品 |
|---|
| 426 | | * |
|---|
| 427 | | * @param m |
|---|
| 428 | | * @return |
|---|
| 429 | | */ |
|---|
| 430 | | public function matrixDisp():int { |
|---|
| 431 | | FLARException.trap("未チェックのパス"); |
|---|
| 432 | | // System.out.println(" === matrix (" + row + "," + clm + ") ===");// printf(" ===matrix (%d,%d) ===\n", m->row, m->clm); |
|---|
| 433 | | // for (var r:int = 0; r < row; r++) {// for(int r = 0; r < m->row; r++) { |
|---|
| 434 | | // System.out.print(" |");// printf(" |"); |
|---|
| 435 | | // for (var c:int = 0; c < clm; c++) {// for(int c = 0; c < m->clm; c++) { |
|---|
| 436 | | // System.out.print(" " + m[r][c]);// printf(" %10g", ARELEM0(m, r, c)); |
|---|
| 437 | | // } |
|---|
| 438 | | // System.out.println(" |");// printf(" |\n"); |
|---|
| 439 | | // } |
|---|
| 440 | | // System.out.println(" ======================");// printf(" ======================\n"); |
|---|
| 441 | | return 0; |
|---|
| 442 | | } |
|---|
| 443 | | |
|---|
| 444 | | private static const PCA_EPS:Number = 1e-6; |
|---|
| 445 | | // #define EPS 1e-6 |
|---|
| 446 | | |
|---|
| 447 | | private static const PCA_MAX_ITER:int = 100; |
|---|
| 448 | | // #define MAX_ITER 100 |
|---|
| 449 | | |
|---|
| 450 | | private static const PCA_VZERO:Number = 1e-16; |
|---|
| 451 | | |
|---|
| 452 | | // #define VZERO 1e-16 |
|---|
| 453 | | |
|---|
| 454 | | /** |
|---|
| 455 | | * static int EX( ARMat *input, ARVec *mean )の代替関数 Optimize:STEP:[144->110] |
|---|
| 456 | | * |
|---|
| 457 | | * @param input |
|---|
| 458 | | * @param mean |
|---|
| 459 | | * @return |
|---|
| 460 | | * @throws FLARException |
|---|
| 461 | | */ |
|---|
| 462 | | private function PCA_EX(mean:FLARVec):void { |
|---|
| 463 | | var lrow:int; |
|---|
| 464 | | var lclm:int; |
|---|
| 465 | | var i:int; |
|---|
| 466 | | var i2:int; |
|---|
| 467 | | lrow = this.row; |
|---|
| 468 | | lclm = this.clm; |
|---|
| 469 | | var lm:Array = this.m; |
|---|
| 470 | | |
|---|
| 471 | | if (lrow <= 0 || lclm <= 0) { |
|---|
| 472 | | throw new FLARException(); |
|---|
| 473 | | } |
|---|
| 474 | | if (mean.getClm() != lclm) { |
|---|
| 475 | | throw new FLARException(); |
|---|
| 476 | | } |
|---|
| 477 | | // double[] mean_array=mean.getArray(); |
|---|
| 478 | | // mean.zeroClear(); |
|---|
| 479 | | const mean_array:Array = mean.getArray(); |
|---|
| 480 | | var w:Number; |
|---|
| 481 | | // For順変更禁止 |
|---|
| 482 | | for (i2 = 0; i2 < lclm; i2++) { |
|---|
| 483 | | w = 0.0; |
|---|
| 484 | | for (i = 0; i < lrow; i++) { |
|---|
| 485 | | // *(v++) += *(m++); |
|---|
| 486 | | w += lm[i][i2]; |
|---|
| 487 | | } |
|---|
| 488 | | mean_array[i2] = w / lrow;// mean->v[i] /= row; |
|---|
| 489 | | } |
|---|
| 490 | | } |
|---|
| 491 | | |
|---|
| 492 | | /** |
|---|
| 493 | | * static int CENTER( ARMat *inout, ARVec *mean )の代替関数 |
|---|
| 494 | | * |
|---|
| 495 | | * @param inout |
|---|
| 496 | | * @param mean |
|---|
| 497 | | * @return |
|---|
| 498 | | */ |
|---|
| 499 | | private static function PCA_CENTER(inout:FLARMat, mean:FLARVec):void { |
|---|
| 500 | | var v:Array; |
|---|
| 501 | | var row:int; |
|---|
| 502 | | var clm:int; |
|---|
| 503 | | |
|---|
| 504 | | row = inout.getRow(); |
|---|
| 505 | | clm = inout.getClm(); |
|---|
| 506 | | if (mean.getClm() != clm) { |
|---|
| 507 | | throw new FLARException(); |
|---|
| 508 | | } |
|---|
| 509 | | var im:Array = inout.m; |
|---|
| 510 | | var im_i:Array; |
|---|
| 511 | | var w0:Number; |
|---|
| 512 | | var w1:Number; |
|---|
| 513 | | v = mean.getArray(); |
|---|
| 514 | | // 特にパフォーマンスが劣化するclm=1と2ときだけ、別パスで処理します。 |
|---|
| 515 | | var i:int; |
|---|
| 516 | | var j:int; |
|---|
| 517 | | switch (clm) { |
|---|
| 518 | | case 1: |
|---|
| 519 | | w0 = v[0]; |
|---|
| 520 | | for (i = 0;i < row; i++) { |
|---|
| 521 | | im[i][0] -= w0; |
|---|
| 522 | | } |
|---|
| 523 | | break; |
|---|
| 524 | | case 2: |
|---|
| 525 | | w0 = v[0]; |
|---|
| 526 | | w1 = v[1]; |
|---|
| 527 | | for (i = 0;i < row; i++) { |
|---|
| 528 | | im_i = im[i]; |
|---|
| 529 | | im_i[0] -= w0; |
|---|
| 530 | | im_i[1] -= w1; |
|---|
| 531 | | } |
|---|
| 532 | | break; |
|---|
| 533 | | default: |
|---|
| 534 | | for (i = 0;i < row; i++) { |
|---|
| 535 | | im_i = im[i]; |
|---|
| 536 | | for (j = 0;j < clm; j++) { |
|---|
| 537 | | // *(m++) -= *(v++); |
|---|
| 538 | | im_i[j] -= v[j]; |
|---|
| 539 | | } |
|---|
| 540 | | } |
|---|
| 541 | | break; |
|---|
| 542 | | } |
|---|
| 543 | | return; |
|---|
| 544 | | } |
|---|
| 545 | | |
|---|
| 546 | | /** |
|---|
| 547 | | * int x_by_xt( ARMat *input, ARMat *output )の代替関数 |
|---|
| 548 | | * |
|---|
| 549 | | * @param input |
|---|
| 550 | | * @param output |
|---|
| 551 | | * @throws FLARException |
|---|
| 552 | | */ |
|---|
| 553 | | private static function PCA_x_by_xt(input:FLARMat, output:FLARMat):void { |
|---|
| 554 | | FLARException.trap("動作未チェック/配列化未チェック"); |
|---|
| 555 | | var row:int; |
|---|
| 556 | | var clm:int; |
|---|
| 557 | | // double[][] out; |
|---|
| 558 | | var in1:Array; |
|---|
| 559 | | var in2:Array; |
|---|
| 560 | | |
|---|
| 561 | | FLARException.trap("未チェックのパス"); |
|---|
| 562 | | row = input.row; |
|---|
| 563 | | clm = input.clm; |
|---|
| 564 | | FLARException.trap("未チェックのパス"); |
|---|
| 565 | | if (output.row != row || output.clm != row) { |
|---|
| 566 | | throw new FLARException(); |
|---|
| 567 | | } |
|---|
| 568 | | |
|---|
| 569 | | // out = output.getArray(); |
|---|
| 570 | | var i:int; |
|---|
| 571 | | var j:int; |
|---|
| 572 | | var k:int; |
|---|
| 573 | | for (i = 0; i < row; i++) { |
|---|
| 574 | | for (j = 0; j < row; j++) { |
|---|
| 575 | | if (j < i) { |
|---|
| 576 | | FLARException.trap("未チェックのパス"); |
|---|
| 577 | | output.m[i][j] = output.m[j][i];// *out = |
|---|
| 578 | | // output->m[j*row+i]; |
|---|
| 579 | | } else { |
|---|
| 580 | | FLARException.trap("未チェックのパス"); |
|---|
| 581 | | in1 = input.m[i]; |
|---|
| 582 | | // input.getRowArray(i);//in1 = &(input->m[clm*i]); |
|---|
| 583 | | in2 = input.m[j]; |
|---|
| 584 | | // input.getRowArray(j);//in2 = &(input->m[clm*j]); |
|---|
| 585 | | output.m[i][j] = 0; |
|---|
| 586 | | // *out = 0.0; |
|---|
| 587 | | for (k = 0; k < clm; k++) { |
|---|
| 588 | | output.m[i][j] += (in1[k] * in2[k]);// *out += *(in1++) |
|---|
| 589 | | // * *(in2++); |
|---|
| 590 | | } |
|---|
| 591 | | } |
|---|
| 592 | | // out.incPtr(); |
|---|
| 593 | | } |
|---|
| 594 | | } |
|---|
| 595 | | } |
|---|
| 596 | | |
|---|
| 597 | | /** |
|---|
| 598 | | * static int xt_by_x( ARMat *input, ARMat *output )の代替関数 |
|---|
| 599 | | * Optimize:2008.04.19 |
|---|
| 600 | | * |
|---|
| 601 | | * @param input |
|---|
| 602 | | * @param i_output |
|---|
| 603 | | * @throws FLARException |
|---|
| 604 | | */ |
|---|
| 605 | | private static function PCA_xt_by_x(input:FLARMat, i_output:FLARMat):void { |
|---|
| 606 | | var in_:Array; |
|---|
| 607 | | var row:int; |
|---|
| 608 | | var clm:int; |
|---|
| 609 | | |
|---|
| 610 | | row = input.row; |
|---|
| 611 | | clm = input.clm; |
|---|
| 612 | | if (i_output.row != clm || i_output.clm != clm) { |
|---|
| 613 | | throw new FLARException(); |
|---|
| 614 | | } |
|---|
| 615 | | |
|---|
| 616 | | var i:int; |
|---|
| 617 | | var k:int; |
|---|
| 618 | | var j:int; |
|---|
| 619 | | var out_m:Array = i_output.m; |
|---|
| 620 | | var w:Number; |
|---|
| 621 | | for (i = 0; i < clm; i++) { |
|---|
| 622 | | for (j = 0; j < clm; j++) { |
|---|
| 623 | | if (j < i) { |
|---|
| 624 | | out_m[i][j] = out_m[j][i];// *out = output->m[j*clm+i]; |
|---|
| 625 | | } else { |
|---|
| 626 | | w = 0.0; |
|---|
| 627 | | // *out = 0.0; |
|---|
| 628 | | for (k = 0; k < row; k++) { |
|---|
| 629 | | in_ = input.m[k]; |
|---|
| 630 | | // in=input.getRowArray(k); |
|---|
| 631 | | w += (in_[i] * in_[j]);// *out += *in1 * *in2; |
|---|
| 632 | | } |
|---|
| 633 | | out_m[i][j] = w; |
|---|
| 634 | | } |
|---|
| 635 | | } |
|---|
| 636 | | } |
|---|
| 637 | | } |
|---|
| 638 | | |
|---|
| 639 | | private const wk_PCA_QRM_ev:FLARVec = new FLARVec(1); |
|---|
| 640 | | |
|---|
| 641 | | /** |
|---|
| 642 | | * static int QRM( ARMat *a, ARVec *dv )の代替関数 |
|---|
| 643 | | * |
|---|
| 644 | | * @param a |
|---|
| 645 | | * @param dv |
|---|
| 646 | | * @throws FLARException |
|---|
| 647 | | */ |
|---|
| 648 | | private function PCA_QRM(dv:FLARVec):void { |
|---|
| 649 | | var w:Number; |
|---|
| 650 | | var t:Number; |
|---|
| 651 | | var s:Number; |
|---|
| 652 | | var x:Number; |
|---|
| 653 | | var y:Number; |
|---|
| 654 | | var c:Number; |
|---|
| 655 | | var dim:int; |
|---|
| 656 | | var iter:int; |
|---|
| 657 | | var dv_array:Array = dv.getArray(); |
|---|
| 658 | | |
|---|
| 659 | | dim = this.row; |
|---|
| 660 | | if (dim != this.clm || dim < 2) { |
|---|
| 661 | | throw new FLARException(); |
|---|
| 662 | | } |
|---|
| 663 | | if (dv.getClm() != dim) { |
|---|
| 664 | | throw new FLARException(); |
|---|
| 665 | | } |
|---|
| 666 | | |
|---|
| 667 | | var ev:FLARVec = this.wk_PCA_QRM_ev; |
|---|
| 668 | | ev.realloc(dim); |
|---|
| 669 | | var ev_array:Array = ev.getArray(); |
|---|
| 670 | | if (ev == null) { |
|---|
| 671 | | throw new FLARException(); |
|---|
| 672 | | } |
|---|
| 673 | | const L_m:Array = this.m; |
|---|
| 674 | | this.vecTridiagonalize(dv, ev, 1); |
|---|
| 675 | | |
|---|
| 676 | | ev_array[0] = 0.0; |
|---|
| 677 | | // ev->v[0] = 0.0; |
|---|
| 678 | | |
|---|
| 679 | | var h:int; |
|---|
| 680 | | var j:int; |
|---|
| 681 | | var k:int; |
|---|
| 682 | | var i:int; |
|---|
| 683 | | for (h = dim - 1; h > 0; h--) { |
|---|
| 684 | | j = h; |
|---|
| 685 | | while (j > 0 && Math.abs(ev_array[j]) > PCA_EPS * (Math.abs(dv_array[j - 1]) + Math.abs(dv_array[j]))) { |
|---|
| 686 | | // while(j>0 && fabs(ev->v[j]) >EPS*(fabs(dv->v[j-1])+fabs(dv->v[j]))) |
|---|
| 687 | | // j--; |
|---|
| 688 | | j--; |
|---|
| 689 | | } |
|---|
| 690 | | if (j == h) { |
|---|
| 691 | | continue; |
|---|
| 692 | | } |
|---|
| 693 | | iter = 0; |
|---|
| 694 | | do { |
|---|
| 695 | | iter++; |
|---|
| 696 | | if (iter > PCA_MAX_ITER) { |
|---|
| 697 | | break; |
|---|
| 698 | | } |
|---|
| 699 | | w = (dv_array[h - 1] - dv_array[h]) / 2; |
|---|
| 700 | | // w = (dv->v[h-1] -dv->v[h]) / 2;//ここ？ |
|---|
| 701 | | t = ev_array[h] * ev_array[h]; |
|---|
| 702 | | // t = ev->v[h] * ev->v[h]; |
|---|
| 703 | | s = Math.sqrt(w * w + t); |
|---|
| 704 | | if (w < 0) { |
|---|
| 705 | | s = -s; |
|---|
| 706 | | } |
|---|
| 707 | | x = dv_array[j] - dv_array[h] + t / (w + s); |
|---|
| 708 | | // x = dv->v[j] -dv->v[h] +t/(w+s); |
|---|
| 709 | | y = ev_array[j + 1]; |
|---|
| 710 | | // y = ev->v[j+1]; |
|---|
| 711 | | for (k = j; k < h; k++) { |
|---|
| 712 | | if (Math.abs(x) >= Math.abs(y)) { |
|---|
| 713 | | if (Math.abs(x) > PCA_VZERO) { |
|---|
| 714 | | t = -y / x; |
|---|
| 715 | | c = 1 / Math.sqrt(t * t + 1); |
|---|
| 716 | | s = t * c; |
|---|
| 717 | | } else { |
|---|
| 718 | | c = 1.0; |
|---|
| 719 | | s = 0.0; |
|---|
| 720 | | } |
|---|
| 721 | | } else { |
|---|
| 722 | | t = -x / y; |
|---|
| 723 | | s = 1.0 / Math.sqrt(t * t + 1); |
|---|
| 724 | | c = t * s; |
|---|
| 725 | | } |
|---|
| 726 | | w = dv_array[k] - dv_array[k + 1]; |
|---|
| 727 | | // w = dv->v[k] -dv->v[k+1]; |
|---|
| 728 | | t = (w * s + 2 * c * ev_array[k + 1]) * s; |
|---|
| 729 | | // t = (w * s +2 * c *ev->v[k+1]) *s; |
|---|
| 730 | | dv_array[k] -= t; |
|---|
| 731 | | // dv->v[k] -= t; |
|---|
| 732 | | dv_array[k + 1] += t; |
|---|
| 733 | | // dv->v[k+1] += t; |
|---|
| 734 | | if (k > j) { |
|---|
| 735 | | FLARException.trap("未チェックパス"); |
|---|
| 736 | | ev_array[k] = c * ev_array[k] - s * y;// ev->v[k]= c *ev->v[k]- s * y; |
|---|
| 737 | | } |
|---|
| 738 | | ev_array[k + 1] += s * (c * w - 2 * s * ev_array[k + 1]); |
|---|
| 739 | | // ev->v[k+1]+= s * (c* w- 2* s *ev->v[k+1]); |
|---|
| 740 | | |
|---|
| 741 | | for (i = 0; i < dim; i++) { |
|---|
| 742 | | x = L_m[k][i]; |
|---|
| 743 | | // x = a->m[k*dim+i]; |
|---|
| 744 | | y = L_m[k + 1][i]; |
|---|
| 745 | | // y = a->m[(k+1)*dim+i]; |
|---|
| 746 | | L_m[k][i] = c * x - s * y; |
|---|
| 747 | | // a->m[k*dim+i] = c * x - s* y; |
|---|
| 748 | | L_m[k + 1][i] = s * x + c * y;// a->m[(k+1)*dim+i] = s* x + c * y; |
|---|
| 749 | | } |
|---|
| 750 | | if (k < h - 1) { |
|---|
| 751 | | FLARException.trap("未チェックパス"); |
|---|
| 752 | | // { |
|---|
| 753 | | x = ev_array[k + 1]; |
|---|
| 754 | | // x = ev->v[k+1]; |
|---|
| 755 | | y = -s * ev_array[k + 2]; |
|---|
| 756 | | // y = -s * ev->v[k+2]; |
|---|
| 757 | | ev_array[k + 2] *= c;// ev->v[k+2] *= c; |
|---|
| 758 | | // } |
|---|
| 759 | | } |
|---|
| 760 | | } |
|---|
| 761 | | } while (Math.abs(ev_array[h]) > PCA_EPS * (Math.abs(dv_array[h - 1]) + Math.abs(dv_array[h]))); |
|---|
| 762 | | } |
|---|
| 763 | | for (k = 0;k < dim - 1; k++) { |
|---|
| 764 | | h = k; |
|---|
| 765 | | t = dv_array[h]; |
|---|
| 766 | | // t = dv->v[h]; |
|---|
| 767 | | for (i = k + 1;i < dim; i++) { |
|---|
| 768 | | if (dv_array[i] > t) { |
|---|
| 769 | | // if( dv->v[i] > t ) { |
|---|
| 770 | | h = i; |
|---|
| 771 | | t = dv_array[h];// t = dv->v[h]; |
|---|
| 772 | | } |
|---|
| 773 | | } |
|---|
| 774 | | dv_array[h] = dv_array[k]; |
|---|
| 775 | | // dv->v[h] = dv->v[k]; |
|---|
| 776 | | dv_array[k] = t; |
|---|
| 777 | | // dv->v[k] = t; |
|---|
| 778 | | this.flipRow(h, k); |
|---|
| 779 | | } |
|---|
| 780 | | } |
|---|
| 781 | | |
|---|
| 782 | | /** |
|---|
| 783 | | * i_row_1番目の行と、i_row_2番目の行を入れ替える。 |
|---|
| 784 | | * |
|---|
| 785 | | * @param i_row_1 |
|---|
| 786 | | * @param i_row_2 |
|---|
| 787 | | */ |
|---|
| 788 | | private function flipRow(i_row_1:int, i_row_2:int):void { |
|---|
| 789 | | var i:int; |
|---|
| 790 | | var w:Number; |
|---|
| 791 | | var r1:Array = this.m[i_row_1]; |
|---|
| 792 | | var r2:Array = this.m[i_row_2]; |
|---|
| 793 | | // For順変更OK |
|---|
| 794 | | for (i = clm - 1;i >= 0; i--) { |
|---|
| 795 | | w = r1[i]; |
|---|
| 796 | | r1[i] = r2[i]; |
|---|
| 797 | | r2[i] = w; |
|---|
| 798 | | } |
|---|
| 799 | | } |
|---|
| 800 | | |
|---|
| 801 | | /** |
|---|
| 802 | | * static int EV_create( ARMat *input, ARMat *u, ARMat *output, ARVec *ev |
|---|
| 803 | | * )の代替関数 |
|---|
| 804 | | * |
|---|
| 805 | | * @param input |
|---|
| 806 | | * @param u |
|---|
| 807 | | * @param output |
|---|
| 808 | | * @param ev |
|---|
| 809 | | * @throws FLARException |
|---|
| 810 | | */ |
|---|
| 811 | | private static function PCA_EV_create(input:FLARMat, u:FLARMat, output:FLARMat, ev:FLARVec):void { |
|---|
| 812 | | FLARException.trap("未チェックのパス"); |
|---|
| 813 | | var row:int; |
|---|
| 814 | | var clm:int; |
|---|
| 815 | | row = input.row; |
|---|
| 816 | | // row = input->row; |
|---|
| 817 | | clm = input.clm; |
|---|
| 818 | | // clm = input->clm; |
|---|
| 819 | | if (row <= 0 || clm <= 0) { |
|---|
| 820 | | throw new FLARException(); |
|---|
| 821 | | } |
|---|
| 822 | | if (u.row != row || u.clm != row) { |
|---|
| 823 | | // if( u->row != row || u->clm != |
|---|
| 824 | | // row ){ |
|---|
| 825 | | throw new FLARException(); |
|---|
| 826 | | } |
|---|
| 827 | | if (output.row != row || output.clm != clm) { |
|---|
| 828 | | // if( output->row != |
|---|
| 829 | | // row || output->clm != |
|---|
| 830 | | // clm ){ |
|---|
| 831 | | throw new FLARException(); |
|---|
| 832 | | } |
|---|
| 833 | | if (ev.getClm() != row) { |
|---|
| 834 | | // if( ev->clm != row ){ |
|---|
| 835 | | throw new FLARException(); |
|---|
| 836 | | } |
|---|
| 837 | | var m:Array; |
|---|
| 838 | | var in_:Array; |
|---|
| 839 | | var m1:Array; |
|---|
| 840 | | var ev_array:Array; |
|---|
| 841 | | var sum:Number; |
|---|
| 842 | | var work:Number; |
|---|
| 843 | | |
|---|
| 844 | | FLARException.trap("未チェックのパス"); |
|---|
| 845 | | m = output.m; |
|---|
| 846 | | // m = output->m; |
|---|
| 847 | | in_ = input.m; |
|---|
| 848 | | |
|---|
| 849 | | var i:int; |
|---|
| 850 | | var j:int; |
|---|
| 851 | | var k:int; |
|---|
| 852 | | ev_array = ev.getArray(); |
|---|
| 853 | | for (i = 0; i < row; i++) { |
|---|
| 854 | | FLARException.trap("未チェックのパス"); |
|---|
| 855 | | if (ev_array[i] < PCA_VZERO) { |
|---|
| 856 | | // if( ev->v[i] < VZERO ){ |
|---|
| 857 | | break; |
|---|
| 858 | | } |
|---|
| 859 | | FLARException.trap("未チェックのパス"); |
|---|
| 860 | | work = 1 / Math.sqrt(Math.abs(ev_array[i])); |
|---|
| 861 | | // work = 1 / |
|---|
| 862 | | // sqrt(fabs(ev->v[i])); |
|---|
| 863 | | for (j = 0; j < clm; j++) { |
|---|
| 864 | | sum = 0.0; |
|---|
| 865 | | m1 = u.m[i]; |
|---|
| 866 | | // m1 = &(u->m[i*row]); |
|---|
| 867 | | // m2=input.getPointer(j);//m2 = &(input->m[j]); |
|---|
| 868 | | for (k = 0; k < row; k++) { |
|---|
| 869 | | sum += m1[k] + in_[k][j];// sum += *m1 * *m2; |
|---|
| 870 | | // m1.incPtr(); //m1++; |
|---|
| 871 | | // m2.addPtr(clm);//m2 += clm; |
|---|
| 872 | | } |
|---|
| 873 | | m1[j] = sum * work;// *(m++) = sum * work; |
|---|
| 874 | | // {//*(m++) = sum * work; |
|---|
| 875 | | // m.set(sum * work); |
|---|
| 876 | | // m.incPtr();} |
|---|
| 877 | | } |
|---|
| 878 | | } |
|---|
| 879 | | for (; i < row; i++) { |
|---|
| 880 | | FLARException.trap("未チェックのパス"); |
|---|
| 881 | | ev_array[i] = 0.0; |
|---|
| 882 | | // ev->v[i] = 0.0; |
|---|
| 883 | | for (j = 0; j < clm; j++) { |
|---|
| 884 | | m[i][j] = 0.0; |
|---|
| 885 | | // m.set(0.0);//*(m++) = 0.0; |
|---|
| 886 | | // m.incPtr(); |
|---|
| 887 | | } |
|---|
| 888 | | } |
|---|
| 889 | | } |
|---|
| 890 | | |
|---|
| 891 | | private var wk_PCA_PCA_u:FLARMat = null; |
|---|
| 892 | | |
|---|
| 893 | | /** |
|---|
| 894 | | * static int PCA( ARMat *input, ARMat *output, ARVec *ev ) |
|---|
| 895 | | * |
|---|
| 896 | | * @param output |
|---|
| 897 | | * @param o_ev |
|---|
| 898 | | * @throws FLARException |
|---|
| 899 | | */ |
|---|
| 900 | | private function PCA_PCA(o_output:FLARMat, o_ev:FLARVec):void { |
|---|
| 901 | | |
|---|
| 902 | | var l_row:int; |
|---|
| 903 | | var l_clm:int; |
|---|
| 904 | | var min:int; |
|---|
| 905 | | var ev_array:Array = o_ev.getArray(); |
|---|
| 906 | | |
|---|
| 907 | | l_row = this.row; |
|---|
| 908 | | // row = input->row; |
|---|
| 909 | | l_clm = this.clm; |
|---|
| 910 | | // clm = input->clm; |
|---|
| 911 | | min = (l_clm < l_row) ? l_clm : l_row; |
|---|
| 912 | | if (l_row < 2 || l_clm < 2) { |
|---|
| 913 | | throw new FLARException(); |
|---|
| 914 | | } |
|---|
| 915 | | if (o_output.clm != this.clm) { |
|---|
| 916 | | // if( output->clm != input->clm ){ |
|---|
| 917 | | throw new FLARException(); |
|---|
| 918 | | } |
|---|
| 919 | | if (o_output.row != min) { |
|---|
| 920 | | // if( output->row != min ){ |
|---|
| 921 | | throw new FLARException(); |
|---|
| 922 | | } |
|---|
| 923 | | if (o_ev.getClm() != min) { |
|---|
| 924 | | // if( ev->clm != min ){ |
|---|
| 925 | | throw new FLARException(); |
|---|
| 926 | | } |
|---|
| 927 | | |
|---|
| 928 | | var u:FLARMat; |
|---|
| 929 | | // u =new FLARMat( min, min ); |
|---|
| 930 | | if (this.wk_PCA_PCA_u == null) { |
|---|
| 931 | | u = new FLARMat(min, min); |
|---|
| 932 | | this.wk_PCA_PCA_u = u; |
|---|
| 933 | | } else { |
|---|
| 934 | | u = this.wk_PCA_PCA_u; |
|---|
| 935 | | u.realloc(min, min); |
|---|
| 936 | | } |
|---|
| 937 | | |
|---|
| 938 | | if (l_row < l_clm) { |
|---|
| 939 | | FLARException.trap("未チェックのパス"); |
|---|
| 940 | | PCA_x_by_xt(this, u);// if(x_by_xt( input, u ) < 0 ) { |
|---|
| 941 | | } else { |
|---|
| 942 | | PCA_xt_by_x(this, u);// if(xt_by_x( input, u ) < 0 ) { |
|---|
| 943 | | } |
|---|
| 944 | | u.PCA_QRM(o_ev); |
|---|
| 945 | | |
|---|
| 946 | | var m1:Array; |
|---|
| 947 | | var m2:Array; |
|---|
| 948 | | if (l_row < l_clm) { |
|---|
| 949 | | FLARException.trap("未チェックのパス"); |
|---|
| 950 | | PCA_EV_create(this, u, o_output, o_ev); |
|---|
| 951 | | } else { |
|---|
| 952 | | m1 = u.m; |
|---|
| 953 | | // m1 = u->m; |
|---|
| 954 | | m2 = o_output.m; |
|---|
| 955 | | // m2 = output->m; |
|---|
| 956 | | var i:int; |
|---|
| 957 | | var j:int; |
|---|
| 958 | | for (i = 0;i < min; i++) { |
|---|
| 959 | | if (ev_array[i] < PCA_VZERO) { |
|---|
| 960 | | // if( ev->v[i] < VZERO ){ |
|---|
| 961 | | break; |
|---|
| 962 | | } |
|---|
| 963 | | for (j = 0;j < min; j++) { |
|---|
| 964 | | m2[i][j] = m1[i][j];// *(m2++) = *(m1++); |
|---|
| 965 | | } |
|---|
| 966 | | } |
|---|
| 967 | | for (;i < min; i++) { |
|---|
| 968 | | // for( ; i < min; i++){ |
|---|
| 969 | | // コードを見た限りあってそうだからコメントアウト(2008/03/26)FLARException.trap("未チェックのパス"); |
|---|
| 970 | | ev_array[i] = 0.0; |
|---|
| 971 | | // ev->v[i] = 0.0; |
|---|
| 972 | | for (j = 0;j < min; j++) { |
|---|
| 973 | | m2[i][j] = 0.0;// *(m2++) = 0.0; |
|---|
| 974 | | } |
|---|
| 975 | | } |
|---|
| 976 | | } |
|---|
| 977 | | } |
|---|
| 978 | | |
|---|
| 979 | | private var wk_work_matrixPCA:FLARMat = null; |
|---|
| 980 | | |
|---|
| 981 | | /** |
|---|
| 982 | | * int arMatrixPCA( ARMat *input, ARMat *evec, ARVec *ev, ARVec *mean ); |
|---|
| 983 | | * 関数の置き換え。input引数がthisになる。 Optimize:2008.04.19 |
|---|
| 984 | | * |
|---|
| 985 | | * @param o_evec |
|---|
| 986 | | * @param o_ev |
|---|
| 987 | | * |
|---|
| 988 | | * @param mean |
|---|
| 989 | | * @throws FLARException |
|---|
| 990 | | */ |
|---|
| 991 | | public function matrixPCA(o_evec:FLARMat, o_ev:FLARVec, mean:FLARVec):void { |
|---|
| 992 | | var srow:Number; |
|---|
| 993 | | var sum:Number; |
|---|
| 994 | | var l_row:int; |
|---|
| 995 | | var l_clm:int; |
|---|
| 996 | | var check:int; |
|---|
| 997 | | |
|---|
| 998 | | l_row = this.row; |
|---|
| 999 | | // row = input->row; |
|---|
| 1000 | | l_clm = this.clm; |
|---|
| 1001 | | // clm = input->clm; |
|---|
| 1002 | | check = (l_row < l_clm) ? l_row : l_clm; |
|---|
| 1003 | | if (l_row < 2 || l_clm < 2) { |
|---|
| 1004 | | throw new FLARException(); |
|---|
| 1005 | | } |
|---|
| 1006 | | if (o_evec.clm != l_clm || o_evec.row != check) { |
|---|
| 1007 | | // if( evec->clm != |
|---|
| 1008 | | // input->clm || |
|---|
| 1009 | | // evec->row != |
|---|
| 1010 | | // check ){ |
|---|
| 1011 | | throw new FLARException(); |
|---|
| 1012 | | } |
|---|
| 1013 | | if (o_ev.getClm() != check) { |
|---|
| 1014 | | // if( ev->clm != check ){ |
|---|
| 1015 | | throw new FLARException(); |
|---|
| 1016 | | } |
|---|
| 1017 | | if (mean.getClm() != l_clm) { |
|---|
| 1018 | | // if( mean->clm != input->clm ){ |
|---|
| 1019 | | throw new FLARException(); |
|---|
| 1020 | | } |
|---|
| 1021 | | |
|---|
| 1022 | | // 自分の内容をワークにコピー(高速化の為に、1度作ったインスタンスは使いまわす) |
|---|
| 1023 | | var work:FLARMat; |
|---|
| 1024 | | if (this.wk_work_matrixPCA == null) { |
|---|
| 1025 | | work = this.matrixAllocDup(); |
|---|
| 1026 | | this.wk_work_matrixPCA = work; |
|---|
| 1027 | | } else { |
|---|
| 1028 | | work = this.wk_work_matrixPCA; |
|---|
| 1029 | | work.matrixDup(this);// arMatrixAllocDup( input );work = |
|---|
| 1030 | | // arMatrixAllocDup( input ); |
|---|
| 1031 | | } |
|---|
| 1032 | | |
|---|
| 1033 | | srow = Math.sqrt(l_row); |
|---|
| 1034 | | work.PCA_EX(mean); |
|---|
| 1035 | | |
|---|
| 1036 | | PCA_CENTER(work, mean); |
|---|
| 1037 | | |
|---|
| 1038 | | var i:int; |
|---|
| 1039 | | var j:int; |
|---|
| 1040 | | // For順変更OK |
|---|
| 1041 | | for (i = 0;i < l_row; i++) { |
|---|
| 1042 | | for (j = 0;j < l_clm; j++) { |
|---|
| 1043 | | work.m[i][j] /= srow;// work->m[i] /= srow; |
|---|
| 1044 | | } |
|---|
| 1045 | | } |
|---|
| 1046 | | |
|---|
| 1047 | | work.PCA_PCA(o_evec, o_ev); |
|---|
| 1048 | | |
|---|
| 1049 | | sum = 0.0; |
|---|
| 1050 | | var ev_array:Array = o_ev.getArray(); |
|---|
| 1051 | | var ev_clm:int = o_ev.getClm(); |
|---|
| 1052 | | // For順変更禁止 |
|---|
| 1053 | | for (i = 0;i < ev_clm; i++) { |
|---|
| 1054 | | // for(int i = 0; i < ev->clm; i++ ){ |
|---|
| 1055 | | sum += ev_array[i];// sum += ev->v[i]; |
|---|
| 1056 | | } |
|---|
| 1057 | | // For順変更禁止 |
|---|
| 1058 | | for (i = 0;i < ev_clm; i++) { |
|---|
| 1059 | | // for(int i = 0; i < ev->clm; i++ ){ |
|---|
| 1060 | | ev_array[i] /= sum;// ev->v[i] /= sum; |
|---|
| 1061 | | } |
|---|
| 1062 | | } |
|---|
| 1063 | | |
|---|
| 1064 | | /* int arMatrixPCA2( ARMat *input, ARMat *evec, ARVec *ev ); */ |
|---|
| 1065 | | public static function arMatrixPCA2(input:FLARMat, evec:FLARMat, ev:FLARVec):void { |
|---|
| 1066 | | FLARException.trap("未チェックのパス"); |
|---|
| 1067 | | var work:FLARMat; |
|---|
| 1068 | | // double srow; // unreferenced |
|---|
| 1069 | | var sum:Number; |
|---|
| 1070 | | var row:int; |
|---|
| 1071 | | var clm:int; |
|---|
| 1072 | | var check:int; |
|---|
| 1073 | | |
|---|
| 1074 | | row = input.row; |
|---|
| 1075 | | // row = input->row; |
|---|
| 1076 | | clm = input.clm; |
|---|
| 1077 | | // clm = input->clm; |
|---|
| 1078 | | check = (row < clm) ? row : clm; |
|---|
| 1079 | | if (row < 2 || clm < 2) { |
|---|
| 1080 | | throw new FLARException(); |
|---|
| 1081 | | } |
|---|
| 1082 | | if (evec.getClm() != input.clm || evec.row != check) { |
|---|
| 1083 | | // if( evec->clm!= input->clm|| evec->row!= check ){ |
|---|
| 1084 | | throw new FLARException(); |
|---|
| 1085 | | } |
|---|
| 1086 | | if (ev.getClm() != check) { |
|---|
| 1087 | | // if( ev->clm != check ){ |
|---|
| 1088 | | throw new FLARException(); |
|---|
| 1089 | | } |
|---|
| 1090 | | |
|---|
| 1091 | | FLARException.trap("未チェックのパス"); |
|---|
| 1092 | | work = input.matrixAllocDup(); |
|---|
| 1093 | | |
|---|
| 1094 | | FLARException.trap("未チェックパス"); |
|---|
| 1095 | | work.PCA_PCA(evec, ev); |
|---|
| 1096 | | // rval = PCA( work, evec, ev ); |
|---|
| 1097 | | sum = 0.0; |
|---|
| 1098 | | var ev_array:Array = ev.getArray(); |
|---|
| 1099 | | |
|---|
| 1100 | | var i:int; |
|---|
| 1101 | | for (i = 0;i < ev.getClm(); i++) { |
|---|
| 1102 | | // for( i = 0; i < ev->clm; i++ |
|---|
| 1103 | | // ){ |
|---|
| 1104 | | FLARException.trap("未チェックパス"); |
|---|
| 1105 | | sum += ev_array[i];// sum += ev->v[i]; |
|---|
| 1106 | | } |
|---|
| 1107 | | for (i = 0;i < ev.getClm(); i++) { |
|---|
| 1108 | | // for(int i = 0; i < ev->clm;i++ ){ |
|---|
| 1109 | | FLARException.trap("未チェックパス"); |
|---|
| 1110 | | ev_array[i] /= sum;// ev->v[i] /= sum; |
|---|
| 1111 | | } |
|---|
| 1112 | | return; |
|---|
| 1113 | | } |
|---|
| 1114 | | |
|---|
| 1115 | | public static function matrixAllocMul(a:FLARMat, b:FLARMat):FLARMat { |
|---|
| 1116 | | FLARException.trap("未チェックのパス"); |
|---|
| 1117 | | var dest:FLARMat = new FLARMat(a.row, b.clm); |
|---|
| 1118 | | FLARException.trap("未チェックのパス"); |
|---|
| 1119 | | dest.matrixMul(a, b); |
|---|
| 1120 | | return dest; |
|---|
| 1121 | | } |
|---|
| 1122 | | |
|---|
| 1123 | | /* static double mdet(double *ap, int dimen, int rowa) */ |
|---|
| 1124 | | private static function Det_mdet(ap:Array, dimen:int, rowa:int):Number { |
|---|
| 1125 | | FLARException.trap("動作未チェック/配列化未チェック"); |
|---|
| 1126 | | var det:Number = 1.0; |
|---|
| 1127 | | var work:Number; |
|---|
| 1128 | | var is_:int = 0; |
|---|
| 1129 | | var mmax:int; |
|---|
| 1130 | | |
|---|
| 1131 | | var k:int; |
|---|
| 1132 | | var i:int; |
|---|
| 1133 | | var j:int; |
|---|
| 1134 | | for (k = 0;k < dimen - 1; k++) { |
|---|
| 1135 | | mmax = k; |
|---|
| 1136 | | for (i = k + 1;i < dimen; i++) { |
|---|
| 1137 | | // if (Math.abs(arMatrixDet_MATRIX_get(ap, i, k, rowa)) > |
|---|
| 1138 | | // Math.abs(arMatrixDet_MATRIX_get(ap, mmax, k, rowa))){ |
|---|
| 1139 | | if (Math.abs(ap[i][k]) > Math.abs(ap[mmax][k])) { |
|---|
| 1140 | | mmax = i; |
|---|
| 1141 | | } |
|---|
| 1142 | | } |
|---|
| 1143 | | if (mmax != k) { |
|---|
| 1144 | | for (j = k;j < dimen; j++) { |
|---|
| 1145 | | work = ap[k][j]; |
|---|
| 1146 | | // work = MATRIX(ap, k, j, rowa); |
|---|
| 1147 | | ap[k][j] = ap[mmax][j]; |
|---|
| 1148 | | // MATRIX(ap, k, j, rowa) =MATRIX(ap, mmax, j, rowa); |
|---|
| 1149 | | ap[mmax][j] = work;// MATRIX(ap, mmax, j, rowa) = work; |
|---|
| 1150 | | } |
|---|
| 1151 | | is_++; |
|---|
| 1152 | | } |
|---|
| 1153 | | for (i = k + 1;i < dimen; i++) { |
|---|
| 1154 | | work = ap[i][k] / ap[k][k]; |
|---|
| 1155 | | // work = arMatrixDet_MATRIX_get(ap,i, k, rowa) /arMatrixDet_MATRIX_get(ap, k, k,rowa); |
|---|
| 1156 | | for (j = k + 1;j < dimen; j++) { |
|---|
| 1157 | | // MATRIX(ap, i, j, rowa) -= work * MATRIX(ap, k, j, rowa); |
|---|
| 1158 | | ap[i][j] -= work * ap[k][j]; |
|---|
| 1159 | | } |
|---|
| 1160 | | } |
|---|
| 1161 | | } |
|---|
| 1162 | | for (i = 0;i < dimen; i++) { |
|---|
| 1163 | | det = ap[i][i];// det *= MATRIX(ap, i, i, rowa); |
|---|
| 1164 | | } |
|---|
| 1165 | | for (i = 0;i < is_; i++) { |
|---|
| 1166 | | det *= -1.0; |
|---|
| 1167 | | } |
|---|
| 1168 | | return det; |
|---|
| 1169 | | } |
|---|
| 1170 | | |
|---|
| 1171 | | /* double arMatrixDet(ARMat *m); */ |
|---|
| 1172 | | public static function arMatrixDet(m:FLARMat):Number { |
|---|
| 1173 | | FLARException.trap("動作未チェック/配列化未チェック"); |
|---|
| 1174 | | if (m.row != m.clm) { |
|---|
| 1175 | | return 0.0; |
|---|
| 1176 | | } |
|---|
| 1177 | | return Det_mdet(m.getArray(), m.row, m.clm);// return mdet(m->m, m->row,m->row); |
|---|
| 1178 | | } |
|---|
| 1179 | | |
|---|
| 1180 | | private const wk_vecTridiagonalize_vec:FLARVec = new FLARVec(0); |
|---|
| 1181 | | |
|---|
| 1182 | | private const wk_vecTridiagonalize_vec2:FLARVec = new FLARVec(0); |
|---|
| 1183 | | |
|---|
| 1184 | | /** |
|---|
| 1185 | | * arVecTridiagonalize関数の代替品 a,d,e間で演算をしてる。何をどうしているかはさっぱりさっぱり |
|---|
| 1186 | | * |
|---|
| 1187 | | * @param a |
|---|
| 1188 | | * @param d |
|---|
| 1189 | | * @param e |
|---|
| 1190 | | * @param i_e_start |
|---|
| 1191 | | * 演算開始列(よくわからないけどarVecTridiagonalizeの呼び出し元でなんかしてる) |
|---|
| 1192 | | * @return |
|---|
| 1193 | | * @throws FLARException |
|---|
| 1194 | | */ |
|---|
| 1195 | | private function vecTridiagonalize(d:FLARVec, e:FLARVec, i_e_start:int):void { |
|---|
| 1196 | | var vec:FLARVec = wk_vecTridiagonalize_vec; |
|---|
| 1197 | | // double[][] a_array=a.getArray(); |
|---|
| 1198 | | var s:Number; |
|---|
| 1199 | | var t:Number; |
|---|
| 1200 | | var p:Number; |
|---|
| 1201 | | var q:Number; |
|---|
| 1202 | | var dim:int; |
|---|
| 1203 | | |
|---|
| 1204 | | if (this.clm != this.row) { |
|---|
| 1205 | | // if(a.getClm()!=a.getRow()){ |
|---|
| 1206 | | throw new FLARException(); |
|---|
| 1207 | | } |
|---|
| 1208 | | if (this.clm != d.getClm()) { |
|---|
| 1209 | | // if(a.getClm() != d.clm){ |
|---|
| 1210 | | throw new FLARException(); |
|---|
| 1211 | | } |
|---|
| 1212 | | if (this.clm != e.getClm()) { |
|---|
| 1213 | | // if(a.getClm() != e.clm){ |
|---|
| 1214 | | throw new FLARException(); |
|---|
| 1215 | | } |
|---|
| 1216 | | dim = this.getClm(); |
|---|
| 1217 | | |
|---|
| 1218 | | var d_vec:Array = d.getArray(); |
|---|
| 1219 | | var e_vec:Array = e.getArray(); |
|---|
| 1220 | | var a_vec_k:Array; |
|---|
| 1221 | | |
|---|
| 1222 | | var k:int; |
|---|
| 1223 | | var i:int; |
|---|
| 1224 | | var j:int; |
|---|
| 1225 | | for (k = 0;k < dim - 2; k++) { |
|---|
| 1226 | | |
|---|
| 1227 | | a_vec_k = this.m[k]; |
|---|
| 1228 | | vec.setNewArray(a_vec_k, clm); |
|---|
| 1229 | | // vec=this.getRowVec(k);//double[] |
|---|
| 1230 | | // vec_array=vec.getArray(); |
|---|
| 1231 | | FLARException.trap("未チェックパス"); |
|---|
| 1232 | | d_vec[k] = a_vec_k[k]; |
|---|
| 1233 | | // d.v[k]=vec.v[k];//d.set(k,v.get(k)); |
|---|
| 1234 | | // //d->v[k] = v[k]; |
|---|
| 1235 | | |
|---|
| 1236 | | // wv1.clm = dim-k-1; |
|---|
| 1237 | | // wv1.v = &(v[k+1]); |
|---|
| 1238 | | FLARException.trap("未チェックパス"); |
|---|
| 1239 | | e_vec[k + i_e_start] = vec.vecHousehold(k + 1); |
|---|
| 1240 | | // e.v[k+i_e_start]=vec.vecHousehold(k+1);//e->v[k]= arVecHousehold(&wv1); |
|---|
| 1241 | | if (e_vec[k + i_e_start] == 0.0) { |
|---|
| 1242 | | // if(e.v[k+i_e_start]== 0.0){//if(e.v[k+i_e_start]== 0.0){ |
|---|
| 1243 | | continue; |
|---|
| 1244 | | } |
|---|
| 1245 | | |
|---|
| 1246 | | for (i = k + 1;i < dim; i++) { |
|---|
| 1247 | | s = 0.0; |
|---|
| 1248 | | for (j = k + 1;j < i; j++) { |
|---|
| 1249 | | FLARException.trap("未チェックのパス"); |
|---|
| 1250 | | s += this.m[j][i] * a_vec_k[j];// s += a_array[j][i] *vec.v[j];//s +=a.get(j*dim+i) *v.get(j);//s +=a->m[j*dim+i] * v[j]; |
|---|
| 1251 | | } |
|---|
| 1252 | | for (j = i;j < dim; j++) { |
|---|
| 1253 | | FLARException.trap("未チェックのパス"); |
|---|
| 1254 | | s += this.m[i][j] * a_vec_k[j];// s += a_array[i][j] *vec.v[j];//s +=a.get(i*dim+j) *v.get(j);//s +=a->m[i*dim+j] * v[j]; |
|---|
| 1255 | | } |
|---|
| 1256 | | FLARException.trap("未チェックのパス"); |
|---|
| 1257 | | d_vec[i] = s;// d.v[i]=s;//d->v[i] = s; |
|---|
| 1258 | | } |
|---|
| 1259 | | |
|---|
| 1260 | | // wv1.clm = wv2.clm = dim-k-1; |
|---|
| 1261 | | // wv1.v = &(v[k+1]); |
|---|
| 1262 | | // wv2.v = &(d->v[k+1]); |
|---|
| 1263 | | a_vec_k = this.m[k]; |
|---|
| 1264 | | vec.setNewArray(a_vec_k, clm); |
|---|
| 1265 | | // vec=this.getRowVec(k); |
|---|
| 1266 | | // vec_array=vec.getArray(); |
|---|
| 1267 | | FLARException.trap("未チェックパス"); |
|---|
| 1268 | | t = vec.vecInnerproduct(d, k + 1) / 2; |
|---|
| 1269 | | for (i = dim - 1;i > k; i--) { |
|---|
| 1270 | | FLARException.trap("未チェックパス"); |
|---|
| 1271 | | p = a_vec_k[i]; |
|---|
| 1272 | | // p = v.get(i);//p = v[i]; |
|---|
| 1273 | | d_vec[i] -= t * p; |
|---|
| 1274 | | q = d_vec[i]; |
|---|
| 1275 | | // d.v[i]-=t*p;q=d.v[i];//q = d->v[i] -= t*p； |
|---|
| 1276 | | for (j = i;j < dim; j++) { |
|---|
| 1277 | | FLARException.trap("未チェックパス"); |
|---|
| 1278 | | this.m[i][j] -= p * (d_vec[j] + q * a_vec_k[j]);// a.m[i][j]-=p*(d.v[j] +q*vec.v[j]);//a->m[i*dim+j] -=p*(d->v[j]) + q*v[j]; |
|---|
| 1279 | | } |
|---|
| 1280 | | } |
|---|
| 1281 | | } |
|---|
| 1282 | | |
|---|
| 1283 | | if (dim >= 2) { |
|---|
| 1284 | | d_vec[dim - 2] = this.m[dim - 2][dim - 2]; |
|---|
| 1285 | | // d.v[dim-2]=a.m[dim-2][dim-2];//d->v[dim-2]=a->m[(dim-2)*dim+(dim-2)]; |
|---|
| 1286 | | e_vec[dim - 2 + i_e_start] = this.m[dim - 2][dim - 1];// e.v[dim-2+i_e_start]=a.m[dim-2][dim-1];//e->v[dim-2] = a->m[(dim-2)*dim+(dim-1)]; |
|---|
| 1287 | | } |
|---|
| 1288 | | |
|---|
| 1289 | | if (dim >= 1) { |
|---|
| 1290 | | d_vec[dim - 1] = this.m[dim - 1][dim - 1];// d.v[dim-1]=a_array[dim-1][dim-1];//d->v[dim-1] =a->m[(dim-1)*dim+(dim-1)]; |
|---|
| 1291 | | } |
|---|
| 1292 | | var vec2:FLARVec = this.wk_vecTridiagonalize_vec2; |
|---|
| 1293 | | for (k = dim - 1;k >= 0; k--) { |
|---|
| 1294 | | a_vec_k = this.m[k]; |
|---|
| 1295 | | vec.setNewArray(a_vec_k, clm); |
|---|
| 1296 | | // vec=this.getRowVec(k);//v =a.getPointer(k*dim);//v = &(a->m[k*dim]); |
|---|
| 1297 | | if (k < dim - 2) { |
|---|
| 1298 | | for (i = k + 1;i < dim; i++) { |
|---|
| 1299 | | // wv1.clm = wv2.clm = dim-k-1; |
|---|
| 1300 | | // wv1.v = &(v[k+1]); |
|---|
| 1301 | | // wv2.v = &(a->m[i*dim+k+1]); |
|---|
| 1302 | | vec2.setNewArray(this.m[i], clm); |
|---|
| 1303 | | // vec2=this.getRowVec(i); |
|---|
| 1304 | | |
|---|
| 1305 | | t = vec.vecInnerproduct(vec2, k + 1); |
|---|
| 1306 | | for (j = k + 1;j < dim; j++) { |
|---|
| 1307 | | FLARException.trap("未チェックパス"); |
|---|
| 1308 | | this.m[i][j] -= t * a_vec_k[j];// a_array[i][j]-=t*vec.v[j];//a.subValue(i*dim+j,t*v.get(j));//a->m[i*dim+j]-= t * v[j]; |
|---|
| 1309 | | } |
|---|
| 1310 | | } |
|---|
| 1311 | | } |
|---|
| 1312 | | for (i = 0;i < dim; i++) { |
|---|
| 1313 | | a_vec_k[i] = 0.0;// v.set(i,0.0);//v[i] = 0.0; |
|---|
| 1314 | | } |
|---|
| 1315 | | a_vec_k[k] = 1;// v.set(k,1);//v[k] = 1; |
|---|
| 1316 | | } |
|---|
| 1317 | | return; |
|---|
| 1318 | | } |
|---|
| 1319 | | } |
|---|
