BVB Source Codes

jsfeat Show jsfeat_linalg.js Source code

Return Download jsfeat: download jsfeat_linalg.js Source code - Download jsfeat Source code - Type:.js
  1. /**
  2.  * @author Eugene Zatepyakin / http://inspirit.ru/
  3.  *
  4.  */
  5.  
  6. (function(global) {
  7.     "use strict";
  8.     //
  9.  
  10.     var linalg = (function() {
  11.  
  12.         var swap = function(A, i0, i1, t) {
  13.             t = A[i0];
  14.             A[i0] = A[i1];
  15.             A[i1] = t;
  16.         }
  17.  
  18.         var hypot = function(a, b) {
  19.             a = Math.abs(a);
  20.             b = Math.abs(b);
  21.             if( a > b ) {
  22.                 b /= a;
  23.                 return a*Math.sqrt(1.0 + b*b);
  24.             }
  25.             if( b > 0 ) {
  26.                 a /= b;
  27.                 return b*Math.sqrt(1.0 + a*a);
  28.             }
  29.             return 0.0;
  30.         }
  31.  
  32.         var JacobiImpl = function(A, astep, W, V, vstep, n) {
  33.             var eps = jsfeat.EPSILON;
  34.             var i=0,j=0,k=0,m=0,l=0,idx=0,_in=0,_in2=0;
  35.             var iters=0,max_iter=n*n*30;
  36.             var mv=0.0,val=0.0,p=0.0,y=0.0,t=0.0,s=0.0,c=0.0,a0=0.0,b0=0.0;
  37.  
  38.             var indR_buff = jsfeat.cache.get_buffer(n<<2);
  39.             var indC_buff = jsfeat.cache.get_buffer(n<<2);
  40.             var indR = indR_buff.i32;
  41.             var indC = indC_buff.i32;
  42.  
  43.             if(V) {
  44.                 for(; i < n; i++) {
  45.                     k = i*vstep;
  46.                     for(j = 0; j < n; j++) {
  47.                         V[k + j] = 0.0;
  48.                     }
  49.                     V[k + i] = 1.0;
  50.                 }
  51.             }
  52.  
  53.             for(k = 0; k < n; k++) {
  54.                 W[k] = A[(astep + 1)*k];
  55.                 if(k < n - 1) {
  56.                     for(m = k+1, mv = Math.abs(A[astep*k + m]), i = k+2; i < n; i++) {
  57.                         val = Math.abs(A[astep*k+i]);
  58.                         if(mv < val)
  59.                             mv = val, m = i;
  60.                     }
  61.                     indR[k] = m;
  62.                 }
  63.                 if(k > 0) {
  64.                     for(m = 0, mv = Math.abs(A[k]), i = 1; i < k; i++) {
  65.                         val = Math.abs(A[astep*i+k]);
  66.                         if(mv < val)
  67.                             mv = val, m = i;
  68.                     }
  69.                     indC[k] = m;
  70.                 }
  71.             }
  72.  
  73.             if(n > 1) for( ; iters < max_iter; iters++) {
  74.                 // find index (k,l) of pivot p
  75.                 for(k = 0, mv = Math.abs(A[indR[0]]), i = 1; i < n-1; i++) {
  76.                     val = Math.abs(A[astep*i + indR[i]]);
  77.                     if( mv < val )
  78.                         mv = val, k = i;
  79.                 }
  80.                 l = indR[k];
  81.                 for(i = 1; i < n; i++) {
  82.                     val = Math.abs(A[astep*indC[i] + i]);
  83.                     if( mv < val )
  84.                         mv = val, k = indC[i], l = i;
  85.                 }
  86.                
  87.                 p = A[astep*k + l];
  88.  
  89.                 if(Math.abs(p) <= eps) break;
  90.  
  91.                 y = (W[l] - W[k])*0.5;
  92.                 t = Math.abs(y) + hypot(p, y);
  93.                 s = hypot(p, t);
  94.                 c = t/s;
  95.                 s = p/s; t = (p/t)*p;
  96.                 if(y < 0)
  97.                     s = -s, t = -t;
  98.                 A[astep*k + l] = 0;
  99.                
  100.                 W[k] -= t;
  101.                 W[l] += t;
  102.                
  103.                 // rotate rows and columns k and l
  104.                 for (i = 0; i < k; i++) {
  105.                     _in = (astep * i + k);
  106.                     _in2 = (astep * i + l);
  107.                     a0 = A[_in];
  108.                     b0 = A[_in2];
  109.                     A[_in] = a0 * c - b0 * s;
  110.                     A[_in2] = a0 * s + b0 * c;
  111.                 }
  112.                 for (i = (k + 1); i < l; i++) {
  113.                     _in = (astep * k + i);
  114.                     _in2 = (astep * i + l);
  115.                     a0 = A[_in];
  116.                     b0 = A[_in2];
  117.                     A[_in] = a0 * c - b0 * s;
  118.                     A[_in2] = a0 * s + b0 * c;
  119.                 }
  120.                 i = l + 1;
  121.                 _in = (astep * k + i);
  122.                 _in2 = (astep * l + i);
  123.                 for (; i < n; i++, _in++, _in2++) {
  124.                     a0 = A[_in];
  125.                     b0 = A[_in2];
  126.                     A[_in] = a0 * c - b0 * s;
  127.                     A[_in2] = a0 * s + b0 * c;
  128.                 }
  129.                
  130.                 // rotate eigenvectors
  131.                 if (V) {
  132.                     _in = vstep * k;
  133.                     _in2 = vstep * l;
  134.                     for (i = 0; i < n; i++, _in++, _in2++) {
  135.                         a0 = V[_in];
  136.                         b0 = V[_in2];
  137.                         V[_in] = a0 * c - b0 * s;
  138.                         V[_in2] = a0 * s + b0 * c;
  139.                     }
  140.                 }
  141.                
  142.                 for(j = 0; j < 2; j++) {
  143.                     idx = j == 0 ? k : l;
  144.                     if(idx < n - 1) {
  145.                         for(m = idx+1, mv = Math.abs(A[astep*idx + m]), i = idx+2; i < n; i++) {
  146.                             val = Math.abs(A[astep*idx+i]);
  147.                             if( mv < val )
  148.                                 mv = val, m = i;
  149.                         }
  150.                         indR[idx] = m;
  151.                     }
  152.                     if(idx > 0) {
  153.                         for(m = 0, mv = Math.abs(A[idx]), i = 1; i < idx; i++) {
  154.                             val = Math.abs(A[astep*i+idx]);
  155.                             if( mv < val )
  156.                                 mv = val, m = i;
  157.                         }
  158.                         indC[idx] = m;
  159.                     }
  160.                 }
  161.             }
  162.  
  163.             // sort eigenvalues & eigenvectors
  164.             for(k = 0; k < n-1; k++) {
  165.                 m = k;
  166.                 for(i = k+1; i < n; i++) {
  167.                     if(W[m] < W[i])
  168.                         m = i;
  169.                 }
  170.                 if(k != m) {
  171.                     swap(W, m, k, mv);
  172.                     if(V) {
  173.                         for(i = 0; i < n; i++) {
  174.                             swap(V, vstep*m + i, vstep*k + i, mv);
  175.                         }
  176.                     }
  177.                 }
  178.             }
  179.  
  180.  
  181.             jsfeat.cache.put_buffer(indR_buff);
  182.             jsfeat.cache.put_buffer(indC_buff);
  183.         }
  184.  
  185.         var JacobiSVDImpl = function(At, astep, _W, Vt, vstep, m, n, n1) {
  186.             var eps = jsfeat.EPSILON * 2.0;
  187.             var minval = jsfeat.FLT_MIN;
  188.             var i=0,j=0,k=0,iter=0,max_iter=Math.max(m, 30);
  189.             var Ai=0,Aj=0,Vi=0,Vj=0,changed=0;
  190.             var c=0.0, s=0.0, t=0.0;
  191.             var t0=0.0,t1=0.0,sd=0.0,beta=0.0,gamma=0.0,delta=0.0,a=0.0,p=0.0,b=0.0;
  192.             var seed = 0x1234;
  193.             var val=0.0,val0=0.0,asum=0.0;
  194.  
  195.             var W_buff = jsfeat.cache.get_buffer(n<<3);
  196.             var W = W_buff.f64;
  197.            
  198.             for(; i < n; i++) {
  199.                 for(k = 0, sd = 0; k < m; k++) {
  200.                     t = At[i*astep + k];
  201.                     sd += t*t;
  202.                 }
  203.                 W[i] = sd;
  204.                
  205.                 if(Vt) {
  206.                     for(k = 0; k < n; k++) {
  207.                         Vt[i*vstep + k] = 0;
  208.                     }
  209.                     Vt[i*vstep + i] = 1;
  210.                 }
  211.             }
  212.            
  213.             for(; iter < max_iter; iter++) {
  214.                 changed = 0;
  215.                
  216.                 for(i = 0; i < n-1; i++) {
  217.                     for(j = i+1; j < n; j++) {
  218.                         Ai = (i*astep)|0, Aj = (j*astep)|0;
  219.                         a = W[i], p = 0, b = W[j];
  220.                        
  221.                         k = 2;
  222.                         p += At[Ai]*At[Aj];
  223.                         p += At[Ai+1]*At[Aj+1];
  224.  
  225.                         for(; k < m; k++)
  226.                             p += At[Ai+k]*At[Aj+k];
  227.                        
  228.                         if(Math.abs(p) <= eps*Math.sqrt(a*b)) continue;
  229.                        
  230.                         p *= 2.0;
  231.                         beta = a - b, gamma = hypot(p, beta);
  232.                         if( beta < 0 ) {
  233.                             delta = (gamma - beta)*0.5;
  234.                             s = Math.sqrt(delta/gamma);
  235.                             c = (p/(gamma*s*2.0));
  236.                         } else {
  237.                             c = Math.sqrt((gamma + beta)/(gamma*2.0));
  238.                             s = (p/(gamma*c*2.0));
  239.                         }
  240.                        
  241.                         a=0.0, b=0.0;
  242.                        
  243.                         k = 2; // unroll
  244.                         t0 = c*At[Ai] + s*At[Aj];
  245.                         t1 = -s*At[Ai] + c*At[Aj];
  246.                         At[Ai] = t0; At[Aj] = t1;
  247.                         a += t0*t0; b += t1*t1;
  248.  
  249.                         t0 = c*At[Ai+1] + s*At[Aj+1];
  250.                         t1 = -s*At[Ai+1] + c*At[Aj+1];
  251.                         At[Ai+1] = t0; At[Aj+1] = t1;
  252.                         a += t0*t0; b += t1*t1;
  253.  
  254.                         for( ; k < m; k++ )
  255.                         {
  256.                             t0 = c*At[Ai+k] + s*At[Aj+k];
  257.                             t1 = -s*At[Ai+k] + c*At[Aj+k];
  258.                             At[Ai+k] = t0; At[Aj+k] = t1;
  259.                            
  260.                             a += t0*t0; b += t1*t1;
  261.                         }
  262.                        
  263.                         W[i] = a; W[j] = b;
  264.                        
  265.                         changed = 1;
  266.                        
  267.                         if(Vt) {
  268.                             Vi = (i*vstep)|0, Vj = (j*vstep)|0;
  269.  
  270.                             k = 2;
  271.                             t0 = c*Vt[Vi] + s*Vt[Vj];
  272.                             t1 = -s*Vt[Vi] + c*Vt[Vj];
  273.                             Vt[Vi] = t0; Vt[Vj] = t1;
  274.  
  275.                             t0 = c*Vt[Vi+1] + s*Vt[Vj+1];
  276.                             t1 = -s*Vt[Vi+1] + c*Vt[Vj+1];
  277.                             Vt[Vi+1] = t0; Vt[Vj+1] = t1;
  278.  
  279.                             for(; k < n; k++) {
  280.                                 t0 = c*Vt[Vi+k] + s*Vt[Vj+k];
  281.                                 t1 = -s*Vt[Vi+k] + c*Vt[Vj+k];
  282.                                 Vt[Vi+k] = t0; Vt[Vj+k] = t1;
  283.                             }
  284.                         }
  285.                     }
  286.                 }
  287.                 if(changed == 0) break;
  288.             }
  289.            
  290.             for(i = 0; i < n; i++) {
  291.                 for(k = 0, sd = 0; k < m; k++) {
  292.                     t = At[i*astep + k];
  293.                     sd += t*t;
  294.                 }
  295.                 W[i] = Math.sqrt(sd);
  296.             }
  297.            
  298.             for(i = 0; i < n-1; i++) {
  299.                 j = i;
  300.                 for(k = i+1; k < n; k++) {
  301.                     if(W[j] < W[k])
  302.                         j = k;
  303.                 }
  304.                 if(i != j) {
  305.                     swap(W, i, j, sd);
  306.                     if(Vt) {
  307.                         for(k = 0; k < m; k++) {
  308.                             swap(At, i*astep + k, j*astep + k, t);
  309.                         }
  310.                        
  311.                         for(k = 0; k < n; k++) {
  312.                             swap(Vt, i*vstep + k, j*vstep + k, t);
  313.                         }
  314.                     }
  315.                 }
  316.             }
  317.            
  318.             for(i = 0; i < n; i++) {
  319.                 _W[i] = W[i];
  320.             }
  321.            
  322.             if(!Vt) {
  323.                 jsfeat.cache.put_buffer(W_buff);
  324.                 return;
  325.             }
  326.  
  327.             for(i = 0; i < n1; i++) {
  328.  
  329.                 sd = i < n ? W[i] : 0;
  330.                
  331.                 while(sd <= minval) {
  332.                     // if we got a zero singular value, then in order to get the corresponding left singular vector
  333.                     // we generate a random vector, project it to the previously computed left singular vectors,
  334.                     // subtract the projection and normalize the difference.
  335.                     val0 = (1.0/m);
  336.                     for(k = 0; k < m; k++) {
  337.                         seed = (seed * 214013 + 2531011);
  338.                         val = (((seed >> 16) & 0x7fff) & 256) != 0 ? val0 : -val0;
  339.                         At[i*astep + k] = val;
  340.                     }
  341.                     for(iter = 0; iter < 2; iter++) {
  342.                         for(j = 0; j < i; j++) {
  343.                             sd = 0;
  344.                             for(k = 0; k < m; k++) {
  345.                                 sd += At[i*astep + k]*At[j*astep + k];
  346.                             }
  347.                             asum = 0.0;
  348.                             for(k = 0; k < m; k++) {
  349.                                 t = (At[i*astep + k] - sd*At[j*astep + k]);
  350.                                 At[i*astep + k] = t;
  351.                                 asum += Math.abs(t);
  352.                             }
  353.                             asum = asum ? 1.0/asum : 0;
  354.                             for(k = 0; k < m; k++) {
  355.                                 At[i*astep + k] *= asum;
  356.                             }
  357.                         }
  358.                     }
  359.                     sd = 0;
  360.                     for(k = 0; k < m; k++) {
  361.                         t = At[i*astep + k];
  362.                         sd += t*t;
  363.                     }
  364.                     sd = Math.sqrt(sd);
  365.                 }
  366.                
  367.                 s = (1.0/sd);
  368.                 for(k = 0; k < m; k++) {
  369.                     At[i*astep + k] *= s;
  370.                 }
  371.             }
  372.  
  373.             jsfeat.cache.put_buffer(W_buff);
  374.         }
  375.        
  376.         return {
  377.  
  378.             lu_solve: function(A, B) {
  379.                 var i=0,j=0,k=0,p=1,astep=A.cols;
  380.                 var ad=A.data, bd=B.data;
  381.                 var t,alpha,d,s;
  382.  
  383.                 for(i = 0; i < astep; i++) {
  384.                     k = i;                    
  385.                     for(j = i+1; j < astep; j++) {
  386.                         if(Math.abs(ad[j*astep + i]) > Math.abs(ad[k*astep+i])) {
  387.                             k = j;
  388.                         }
  389.                     }
  390.                    
  391.                     if(Math.abs(ad[k*astep+i]) < jsfeat.EPSILON) {
  392.                         return 0; // FAILED
  393.                     }
  394.                    
  395.                     if(k != i) {
  396.                         for(j = i; j < astep; j++ ) {
  397.                             swap(ad, i*astep+j, k*astep+j, t);
  398.                         }
  399.                        
  400.                         swap(bd, i, k, t);
  401.                         p = -p;
  402.                     }
  403.                    
  404.                     d = -1.0/ad[i*astep+i];
  405.                    
  406.                     for(j = i+1; j < astep; j++) {
  407.                         alpha = ad[j*astep+i]*d;
  408.                        
  409.                         for(k = i+1; k < astep; k++) {
  410.                             ad[j*astep+k] += alpha*ad[i*astep+k];
  411.                         }
  412.                        
  413.                         bd[j] += alpha*bd[i];
  414.                     }
  415.                    
  416.                     ad[i*astep+i] = -d;
  417.                 }
  418.                
  419.                 for(i = astep-1; i >= 0; i--) {
  420.                     s = bd[i];
  421.                     for(k = i+1; k < astep; k++) {
  422.                         s -= ad[i*astep+k]*bd[k];
  423.                     }
  424.                     bd[i] = s*ad[i*astep+i];
  425.                 }
  426.  
  427.                 return 1; // OK
  428.             },
  429.  
  430.             cholesky_solve: function(A, B) {
  431.                 var col=0,row=0,col2=0,cs=0,rs=0,i=0,j=0;
  432.                 var size = A.cols;
  433.                 var ad=A.data, bd=B.data;
  434.                 var val,inv_diag;
  435.  
  436.                 for (col = 0; col < size; col++) {
  437.                     inv_diag = 1.0;
  438.                     cs = (col * size);
  439.                     rs = cs;
  440.                     for (row = col; row < size; row++)
  441.                     {
  442.                         // correct for the parts of cholesky already computed
  443.                         val = ad[(rs+col)];
  444.                         for (col2 = 0; col2 < col; col2++) {
  445.                             val -= ad[(col2*size+col)] * ad[(rs+col2)];
  446.                         }
  447.                         if (row == col) {
  448.                             // this is the diagonal element so don't divide
  449.                             ad[(rs+col)] = val;
  450.                             if(val == 0) {
  451.                                 return 0;
  452.                             }
  453.                             inv_diag = 1.0 / val;
  454.                         } else {
  455.                             // cache the value without division in the upper half
  456.                             ad[(cs+row)] = val;
  457.                             // divide my the diagonal element for all others
  458.                             ad[(rs+col)] = val * inv_diag;
  459.                         }
  460.                         rs = (rs + size);
  461.                     }
  462.                 }
  463.  
  464.                 // first backsub through L
  465.                 cs = 0;
  466.                 for (i = 0; i < size; i++) {
  467.                     val = bd[i];
  468.                     for (j = 0; j < i; j++) {
  469.                         val -= ad[(cs+j)] * bd[j];
  470.                     }
  471.                     bd[i] = val;
  472.                     cs = (cs + size);
  473.                 }
  474.                 // backsub through diagonal
  475.                 cs = 0;
  476.                 for (i = 0; i < size; i++) {
  477.                     bd[i] /= ad[(cs + i)];
  478.                     cs = (cs + size);
  479.                 }
  480.                 // backsub through L Transpose
  481.                 i = (size-1);
  482.                 for (; i >= 0; i--) {
  483.                     val = bd[i];
  484.                     j = (i + 1);
  485.                     cs = (j * size);
  486.                     for (; j < size; j++) {
  487.                         val -= ad[(cs + i)] * bd[j];
  488.                         cs = (cs + size);
  489.                     }
  490.                     bd[i] = val;
  491.                 }
  492.  
  493.                 return 1;
  494.             },
  495.  
  496.             svd_decompose: function(A, W, U, V, options) {
  497.                 if (typeof options === "undefined") { options = 0; };
  498.                 var at=0,i=0,j=0,_m=A.rows,_n=A.cols,m=_m,n=_n;
  499.                 var dt = A.type | jsfeat.C1_t; // we only work with single channel
  500.  
  501.                 if(m < n) {
  502.                     at = 1;
  503.                     i = m;
  504.                     m = n;
  505.                     n = i;
  506.                 }
  507.  
  508.                 var a_buff = jsfeat.cache.get_buffer((m*m)<<3);
  509.                 var w_buff = jsfeat.cache.get_buffer(n<<3);
  510.                 var v_buff = jsfeat.cache.get_buffer((n*n)<<3);
  511.  
  512.                 var a_mt = new jsfeat.matrix_t(m, m, dt, a_buff.data);
  513.                 var w_mt = new jsfeat.matrix_t(1, n, dt, w_buff.data);
  514.                 var v_mt = new jsfeat.matrix_t(n, n, dt, v_buff.data);
  515.  
  516.                 if(at == 0) {
  517.                     // transpose
  518.                     jsfeat.matmath.transpose(a_mt, A);
  519.                 } else {
  520.                     for(i = 0; i < _n*_m; i++) {
  521.                         a_mt.data[i] = A.data[i];
  522.                     }
  523.                     for(; i < n*m; i++) {
  524.                         a_mt.data[i] = 0;
  525.                     }
  526.                 }
  527.  
  528.                 JacobiSVDImpl(a_mt.data, m, w_mt.data, v_mt.data, n, m, n, m);
  529.  
  530.                 if(W) {
  531.                     for(i=0; i < n; i++) {
  532.                         W.data[i] = w_mt.data[i];
  533.                     }
  534.                     for(; i < _n; i++) {
  535.                         W.data[i] = 0;
  536.                     }
  537.                 }
  538.  
  539.                 if (at == 0) {
  540.                     if(U && (options & jsfeat.SVD_U_T)) {
  541.                         i = m*m;
  542.                         while(--i >= 0) {
  543.                             U.data[i] = a_mt.data[i];
  544.                         }
  545.                     } else if(U) {
  546.                         jsfeat.matmath.transpose(U, a_mt);
  547.                     }
  548.  
  549.                     if(V && (options & jsfeat.SVD_V_T)) {
  550.                         i = n*n;
  551.                         while(--i >= 0) {
  552.                             V.data[i] = v_mt.data[i];
  553.                         }
  554.                     } else if(V) {
  555.                         jsfeat.matmath.transpose(V, v_mt);
  556.                     }
  557.                 } else {
  558.                     if(U && (options & jsfeat.SVD_U_T)) {
  559.                         i = n*n;
  560.                         while(--i >= 0) {
  561.                             U.data[i] = v_mt.data[i];
  562.                         }
  563.                     } else if(U) {
  564.                         jsfeat.matmath.transpose(U, v_mt);
  565.                     }
  566.  
  567.                     if(V && (options & jsfeat.SVD_V_T)) {
  568.                         i = m*m;
  569.                         while(--i >= 0) {
  570.                             V.data[i] = a_mt.data[i];
  571.                         }
  572.                     } else if(V) {
  573.                         jsfeat.matmath.transpose(V, a_mt);
  574.                     }
  575.                 }
  576.  
  577.                 jsfeat.cache.put_buffer(a_buff);
  578.                 jsfeat.cache.put_buffer(w_buff);
  579.                 jsfeat.cache.put_buffer(v_buff);
  580.  
  581.             },
  582.  
  583.             svd_solve: function(A, X, B) {
  584.                 var i=0,j=0,k=0;
  585.                 var pu=0,pv=0;
  586.                 var nrows=A.rows,ncols=A.cols;
  587.                 var sum=0.0,xsum=0.0,tol=0.0;
  588.                 var dt = A.type | jsfeat.C1_t;
  589.  
  590.                 var u_buff = jsfeat.cache.get_buffer((nrows*nrows)<<3);
  591.                 var w_buff = jsfeat.cache.get_buffer(ncols<<3);
  592.                 var v_buff = jsfeat.cache.get_buffer((ncols*ncols)<<3);
  593.  
  594.                 var u_mt = new jsfeat.matrix_t(nrows, nrows, dt, u_buff.data);
  595.                 var w_mt = new jsfeat.matrix_t(1, ncols, dt, w_buff.data);
  596.                 var v_mt = new jsfeat.matrix_t(ncols, ncols, dt, v_buff.data);
  597.  
  598.                 var bd = B.data, ud = u_mt.data, wd = w_mt.data, vd = v_mt.data;
  599.  
  600.                 this.svd_decompose(A, w_mt, u_mt, v_mt, 0);
  601.  
  602.                 tol = jsfeat.EPSILON * wd[0] * ncols;
  603.  
  604.                 for (; i < ncols; i++, pv += ncols) {
  605.                     xsum = 0.0;
  606.                     for(j = 0; j < ncols; j++) {
  607.                         if(wd[j] > tol) {
  608.                             for(k = 0, sum = 0.0, pu = 0; k < nrows; k++, pu += ncols) {
  609.                                 sum += ud[pu + j] * bd[k];
  610.                             }
  611.                             xsum += sum * vd[pv + j] / wd[j];
  612.                         }
  613.                     }
  614.                     X.data[i] = xsum;
  615.                 }
  616.  
  617.                 jsfeat.cache.put_buffer(u_buff);
  618.                 jsfeat.cache.put_buffer(w_buff);
  619.                 jsfeat.cache.put_buffer(v_buff);
  620.             },
  621.  
  622.             svd_invert: function(Ai, A) {
  623.                 var i=0,j=0,k=0;
  624.                 var pu=0,pv=0,pa=0;
  625.                 var nrows=A.rows,ncols=A.cols;
  626.                 var sum=0.0,tol=0.0;
  627.                 var dt = A.type | jsfeat.C1_t;
  628.  
  629.                 var u_buff = jsfeat.cache.get_buffer((nrows*nrows)<<3);
  630.                 var w_buff = jsfeat.cache.get_buffer(ncols<<3);
  631.                 var v_buff = jsfeat.cache.get_buffer((ncols*ncols)<<3);
  632.  
  633.                 var u_mt = new jsfeat.matrix_t(nrows, nrows, dt, u_buff.data);
  634.                 var w_mt = new jsfeat.matrix_t(1, ncols, dt, w_buff.data);
  635.                 var v_mt = new jsfeat.matrix_t(ncols, ncols, dt, v_buff.data);
  636.  
  637.                 var id = Ai.data, ud = u_mt.data, wd = w_mt.data, vd = v_mt.data;
  638.  
  639.                 this.svd_decompose(A, w_mt, u_mt, v_mt, 0);
  640.  
  641.                 tol = jsfeat.EPSILON * wd[0] * ncols;
  642.  
  643.                 for (; i < ncols; i++, pv += ncols) {
  644.                     for (j = 0, pu = 0; j < nrows; j++, pa++) {
  645.                         for (k = 0, sum = 0.0; k < ncols; k++, pu++) {
  646.                             if (wd[k] > tol) sum += vd[pv + k] * ud[pu] / wd[k];
  647.                         }
  648.                         id[pa] = sum;
  649.                     }
  650.                 }
  651.  
  652.                 jsfeat.cache.put_buffer(u_buff);
  653.                 jsfeat.cache.put_buffer(w_buff);
  654.                 jsfeat.cache.put_buffer(v_buff);
  655.             },
  656.  
  657.             eigenVV: function(A, vects, vals) {
  658.                 var n=A.cols,i=n*n;
  659.                 var dt = A.type | jsfeat.C1_t;
  660.  
  661.                 var a_buff = jsfeat.cache.get_buffer((n*n)<<3);
  662.                 var w_buff = jsfeat.cache.get_buffer(n<<3);
  663.                 var a_mt = new jsfeat.matrix_t(n, n, dt, a_buff.data);
  664.                 var w_mt = new jsfeat.matrix_t(1, n, dt, w_buff.data);
  665.  
  666.                 while(--i >= 0) {
  667.                     a_mt.data[i] = A.data[i];
  668.                 }
  669.  
  670.                 JacobiImpl(a_mt.data, n, w_mt.data, vects ? vects.data : null, n, n);
  671.  
  672.                 if(vals) {
  673.                     while(--n >= 0) {
  674.                         vals.data[n] = w_mt.data[n];
  675.                     }
  676.                 }
  677.  
  678.                 jsfeat.cache.put_buffer(a_buff);
  679.                 jsfeat.cache.put_buffer(w_buff);
  680.             }
  681.  
  682.         };
  683.  
  684.     })();
  685.  
  686.     global.linalg = linalg;
  687.  
  688. })(jsfeat);
downloadjsfeat_linalg.js Source code - Download jsfeat Source code
Related Source Codes/Software:
flakes - Flakes is an Admin Template Framework. A combinati... 2017-04-15
capstone - Capstone disassembly/disassembler framework: Core ... 2017-04-15
nginx-resources - A collection of resources covering Nginx, Nginx + ... 2017-04-15
utron - A lightweight MVC framework for Go(Golang) 2017-04-15
cfssl - CFSSL: Cloudflare's PKI and TLS toolkit ... 2017-04-15
MLeaksFinder - Find memory leaks in your iOS app at develop time. 2017-04-16
qt - Qt binding for Go (Golang) which supports Windows ... 2017-04-16
rainloop-webmail - Simple, modern & fast web-based email client ... 2017-04-16
pencil - Multiplatform GUI Prototyping/Wireframing 2017-04-16
x64dbg - An open-source x64/x32 debugger for windows. ... 2017-04-16
Toucan - Fabulous Image Processing in Swift 2017-04-23
CoffeeScriptRedux - 2017-04-23
breakpoint - Really simple media queries in Sa 2017-04-23
libsvm - 2017-04-22
grr - GRR Rapid Response: remote live forensics for inci... 2017-04-22
grit - **Grit is no longer maintained. Check out libgit2/... 2017-04-22
guard-livereload - Guard::LiveReload automatically reload your browse... 2017-04-22
Begin-Latex-in-minutes - Brief Intro to LaTeX for beginners that helps you ... 2017-04-22
wicked - Use wicked to turn your controller into a wizard ... 2017-04-22
flexboxfroggy - A game for learning CSS flexbox ... 2017-04-22

 Back to top