forked from: forked from: flash on 2014-3-18
forked from forked from: flash on 2014-3-18 (diff: 164)
ActionScript3 source code
/**
* Copyright imo_ ( http://wonderfl.net/user/imo_ )
* MIT License ( http://www.opensource.org/licenses/mit-license.php )
* Downloaded from: http://wonderfl.net/c/crXb
*/
// forked from imo_'s forked from: flash on 2014-3-18
// forked from imo_'s flash on 2014-3-18
package {
import flash.utils.getTimer;
import flash.text.TextFormat;
import flash.text.TextField;
import flash.display.Sprite;
import com.bit101.components.*;
public class GenInvMat extends Sprite {
private var ta:TextArea;
public function GenInvMat() {
Style.embedFonts = false;
Style.fontName = "Courier New";
Style.fontSize = 12;
ta = new TextArea(this, 0, 0, "");
ta.editable = false;
ta.width = 465;
ta.height = 465;
calc();
}
private function join(...nums:Array):String {
var len:int = nums.length;
var txt:String = "| ";
for (var i:int = 0; i < len; i++) {
if (i != 0) txt += " ";
txt += Number(nums[i]).toFixed(4);
}
txt += " |";
return txt;
}
private function log(txt:Object):void {
if (txt is Mat) {
var m:Mat = Mat(txt);
log(join(m.e00, m.e01, m.e02, m.e03));
log(join(m.e10, m.e11, m.e12, m.e13));
log(join(m.e20, m.e21, m.e22, m.e23));
log(join(m.e30, m.e31, m.e32, m.e33));
return;
}
ta.text += txt + "\n";
}
private function calc():void {
log("Moore-Penrose Pseudoinverse Calculation Test");
log("");
var mat:Mat = new Mat(); // a matrix whose det = 0
mat.e00 = 2; mat.e01 = 1; mat.e02 = 1; mat.e03 = 0;
mat.e10 = 1; mat.e11 = 2; mat.e12 = 0; mat.e13 = 1;
mat.e20 = 1; mat.e21 = 0; mat.e22 = 2; mat.e23 = 1;
mat.e30 = 0; mat.e31 = 1; mat.e32 = 1; mat.e33 = 2;
/*
mat.e00 = 1; mat.e01 = 1; mat.e02 = 1; mat.e03 = 1;
mat.e10 = 1; mat.e11 = 1; mat.e12 = 1; mat.e13 = 1;
mat.e20 = 1; mat.e21 = 1; mat.e22 = 1; mat.e23 = 1;
mat.e30 = 1; mat.e31 = 1; mat.e32 = 1; mat.e33 = 1;
*/
var invMat:Mat = calcMain(mat);
/*
// benchmark
var i:int;
var st:int;
var en:int;
st = getTimer();
for (i = 0; i < 10000000; i++) compute(mat, invMat);
en = getTimer();
log("time: " + (en - st) + "ms");
st = getTimer();
for (i = 0; i < 10000000; i++) compute(mat, invMat);
en = getTimer();
log("time: " + (en - st) + "ms");
st = getTimer();
for (i = 0; i < 10000000; i++) compute(mat, invMat);
en = getTimer();
log("time: " + (en - st) + "ms");
*/
var tmp:Mat = new Mat();
log("M: ");
log(mat);
log("");
log("M+: ");
log(invMat);
log("");
log("MM+ (= (MM+)*): ");
tmp.mul(mat, invMat);
log(tmp);
log("");
log("M+M (= (M+M)*): ");
tmp.mul(invMat, mat);
log(tmp);
log("");
log("MM+M (= M): ");
tmp.mul(mat, invMat);
tmp.mul(tmp, mat);
log(tmp);
log("");
log("M+MM+ (= M+): ");
tmp.mul(invMat, mat);
tmp.mul(tmp, invMat);
log(tmp);
log("");
}
private function calcMain(m:Mat):Mat {
log("calculating...");
var inv:Mat = new Mat();
compute(m, inv);
log("done.");
log("");
return inv;
}
private function compute(m:Mat, inv:Mat):void {
// --- variables ---
var len:Number;
var d0:Number;
var d1:Number;
var d2:Number;
var i00:Number;
var i01:Number;
var i02:Number;
var i03:Number;
var i10:Number;
var i11:Number;
var i12:Number;
var i13:Number;
var i20:Number;
var i21:Number;
var i22:Number;
var i23:Number;
var i30:Number;
var i31:Number;
var i32:Number;
var i33:Number;
var e00:Number;
var e01:Number;
var e02:Number;
var e03:Number;
var e10:Number;
var e11:Number;
var e12:Number;
var e13:Number;
var e20:Number;
var e21:Number;
var e22:Number;
var e23:Number;
var e30:Number;
var e31:Number;
var e32:Number;
var e33:Number;
// --- input ---
e00 = m.e00;
e01 = m.e01;
e02 = m.e02;
e03 = m.e03;
e10 = m.e10;
e11 = m.e11;
e12 = m.e12;
e13 = m.e13;
e20 = m.e20;
e21 = m.e21;
e22 = m.e22;
e23 = m.e23;
e30 = m.e30;
e31 = m.e31;
e32 = m.e32;
e33 = m.e33;
// --- compute ---
// <oimo's method>
// ---------------------------------------------
// For an m×2 matrix M,
//
// O if M*M = O
// M+ = inv(M*M)M* if M*M ≠ O, det|M*M| ≠ 0
// M* / tr(M*M) if M*M ≠ O, det|M*M| = 0
// ---------------------------------------------
// calc A2+
if ((d0 = e00 * e00 + e10 * e10 + e20 * e20 + e30 * e30) + (d1 = e01 * e01 + e11 * e11 + e21 * e21 + e31 * e31) > 1e-10) { // A2*A2 ≠ O
d2 = e00 * e01 + e10 * e11 + e20 * e21 + e30 * e31;
if ((len = d0 * d1 - d2 * d2) > 1e-10 || len < -1e-10) { // det|A2*A2| ≠ 0
// A2+ = inv(A2*A2)A2*
len = 1 / len;
d0 *= len;
d1 *= len;
d2 *= len;
i00 = d1 * e00 - d2 * e01;
i01 = d1 * e10 - d2 * e11;
i02 = d1 * e20 - d2 * e21;
i03 = d1 * e30 - d2 * e31;
i10 = d0 * e01 - d2 * e00;
i11 = d0 * e11 - d2 * e10;
i12 = d0 * e21 - d2 * e20;
i13 = d0 * e31 - d2 * e30;
} else { // det|A2*A2| = 0
// A2+ = A2* / tr(A2*A2)
len = 1 / (d0 + d1);
i00 = e00 * len;
i01 = e10 * len;
i02 = e20 * len;
i03 = e30 * len;
i10 = e01 * len;
i11 = e11 * len;
i12 = e21 * len;
i13 = e31 * len;
}
} else { // A2*A2 = O
// A2+ = O
i00 = 0;
i01 = 0;
i02 = 0;
i03 = 0;
i10 = 0;
i11 = 0;
i12 = 0;
i13 = 0;
}
// </oimo's method>
// <optimized greville's method>
// calc A3+
d0 = i00 * e02 + i01 * e12 + i02 * e22 + i03 * e32;
d1 = i10 * e02 + i11 * e12 + i12 * e22 + i13 * e32;
i20 = e02 - e00 * d0 - e01 * d1;
i21 = e12 - e10 * d0 - e11 * d1;
i22 = e22 - e20 * d0 - e21 * d1;
i23 = e32 - e30 * d0 - e31 * d1;
if ((len = i20 * i20 + i21 * i21 + i22 * i22 + i23 * i23) > 1e-10) {
len = 1 / len;
i20 *= len;
i21 *= len;
i22 *= len;
i23 *= len;
} else {
len = 1 / (1 + d0 * d0 + d1 * d1);
i20 = (d0 * i00 + d1 * i10) * len;
i21 = (d0 * i01 + d1 * i11) * len;
i22 = (d0 * i02 + d1 * i12) * len;
i23 = (d0 * i03 + d1 * i13) * len;
}
i00 -= d0 * i20;
i01 -= d0 * i21;
i02 -= d0 * i22;
i03 -= d0 * i23;
i10 -= d1 * i20;
i11 -= d1 * i21;
i12 -= d1 * i22;
i13 -= d1 * i23;
// calc A4+
d0 = i00 * e03 + i01 * e13 + i02 * e23 + i03 * e33;
d1 = i10 * e03 + i11 * e13 + i12 * e23 + i13 * e33;
d2 = i20 * e03 + i21 * e13 + i22 * e23 + i23 * e33;
i30 = e03 - e00 * d0 - e01 * d1 - e02 * d2;
i31 = e13 - e10 * d0 - e11 * d1 - e12 * d2;
i32 = e23 - e20 * d0 - e21 * d1 - e22 * d2;
i33 = e33 - e30 * d0 - e31 * d1 - e32 * d2;
if ((len = i30 * i30 + i31 * i31 + i32 * i32 + i33 * i33) > 1e-10) {
len = 1 / len;
i30 *= len;
i31 *= len;
i32 *= len;
i33 *= len;
} else {
len = 1 / (1 + d0 * d0 + d1 * d1 + d2 * d2);
i30 = (d0 * i00 + d1 * i10 + d2 * i20) * len;
i31 = (d0 * i01 + d1 * i11 + d2 * i21) * len;
i32 = (d0 * i02 + d1 * i12 + d2 * i22) * len;
i33 = (d0 * i03 + d1 * i13 + d2 * i23) * len;
}
i00 -= d0 * i30;
i01 -= d0 * i31;
i02 -= d0 * i32;
i03 -= d0 * i33;
i10 -= d1 * i30;
i11 -= d1 * i31;
i12 -= d1 * i32;
i13 -= d1 * i33;
i20 -= d2 * i30;
i21 -= d2 * i31;
i22 -= d2 * i32;
i23 -= d2 * i33;
// </optimized greville's method>
// --- output ---
inv.e00 = i00;
inv.e01 = i01;
inv.e02 = i02;
inv.e03 = i03;
inv.e10 = i10;
inv.e11 = i11;
inv.e12 = i12;
inv.e13 = i13;
inv.e20 = i20;
inv.e21 = i21;
inv.e22 = i22;
inv.e23 = i23;
inv.e30 = i30;
inv.e31 = i31;
inv.e32 = i32;
inv.e33 = i33;
}
}
}
class Mat {
public var e00:Number;
public var e01:Number;
public var e02:Number;
public var e03:Number;
public var e10:Number;
public var e11:Number;
public var e12:Number;
public var e13:Number;
public var e20:Number;
public var e21:Number;
public var e22:Number;
public var e23:Number;
public var e30:Number;
public var e31:Number;
public var e32:Number;
public var e33:Number;
public function Mat() {
e00 = 1; e01 = 0; e02 = 0; e03 = 0;
e10 = 0; e11 = 1; e12 = 0; e13 = 0;
e20 = 0; e21 = 0; e22 = 1; e23 = 0;
e30 = 0; e31 = 0; e32 = 0; e33 = 1;
}
public function mul(m1:Mat, m2:Mat):void {
var t00:Number = m1.e00 * m2.e00 + m1.e01 * m2.e10 + m1.e02 * m2.e20 + m1.e03 * m2.e30;
var t01:Number = m1.e00 * m2.e01 + m1.e01 * m2.e11 + m1.e02 * m2.e21 + m1.e03 * m2.e31;
var t02:Number = m1.e00 * m2.e02 + m1.e01 * m2.e12 + m1.e02 * m2.e22 + m1.e03 * m2.e32;
var t03:Number = m1.e00 * m2.e03 + m1.e01 * m2.e13 + m1.e02 * m2.e23 + m1.e03 * m2.e33;
var t10:Number = m1.e10 * m2.e00 + m1.e11 * m2.e10 + m1.e12 * m2.e20 + m1.e13 * m2.e30;
var t11:Number = m1.e10 * m2.e01 + m1.e11 * m2.e11 + m1.e12 * m2.e21 + m1.e13 * m2.e31;
var t12:Number = m1.e10 * m2.e02 + m1.e11 * m2.e12 + m1.e12 * m2.e22 + m1.e13 * m2.e32;
var t13:Number = m1.e10 * m2.e03 + m1.e11 * m2.e13 + m1.e12 * m2.e23 + m1.e13 * m2.e33;
var t20:Number = m1.e20 * m2.e00 + m1.e21 * m2.e10 + m1.e22 * m2.e20 + m1.e23 * m2.e30;
var t21:Number = m1.e20 * m2.e01 + m1.e21 * m2.e11 + m1.e22 * m2.e21 + m1.e23 * m2.e31;
var t22:Number = m1.e20 * m2.e02 + m1.e21 * m2.e12 + m1.e22 * m2.e22 + m1.e23 * m2.e32;
var t23:Number = m1.e20 * m2.e03 + m1.e21 * m2.e13 + m1.e22 * m2.e23 + m1.e23 * m2.e33;
var t30:Number = m1.e30 * m2.e00 + m1.e31 * m2.e10 + m1.e32 * m2.e20 + m1.e33 * m2.e30;
var t31:Number = m1.e30 * m2.e01 + m1.e31 * m2.e11 + m1.e32 * m2.e21 + m1.e33 * m2.e31;
var t32:Number = m1.e30 * m2.e02 + m1.e31 * m2.e12 + m1.e32 * m2.e22 + m1.e33 * m2.e32;
var t33:Number = m1.e30 * m2.e03 + m1.e31 * m2.e13 + m1.e32 * m2.e23 + m1.e33 * m2.e33;
e00 = t00;
e01 = t01;
e02 = t02;
e03 = t03;
e10 = t10;
e11 = t11;
e12 = t12;
e13 = t13;
e20 = t20;
e21 = t21;
e22 = t22;
e23 = t23;
e30 = t30;
e31 = t31;
e32 = t32;
e33 = t33;
}
}
