simple-squiggle

A restricted subset of Squiggle
Log | Files | Refs | README

csCounts.js (3396B)

      1 import { factory } from '../../../utils/factory.js';
      2 import { csLeaf } from './csLeaf.js';
      3 var name = 'csCounts';
      4 var dependencies = ['transpose'];
      5 export var createCsCounts = /* #__PURE__ */factory(name, dependencies, _ref => {
      6   var {
      7     transpose
      8   } = _ref;
      9 
     10   /**
     11    * Computes the column counts using the upper triangular part of A.
     12    * It transposes A internally, none of the input parameters are modified.
     13    *
     14    * @param {Matrix} a           The sparse matrix A
     15    *
     16    * @param {Matrix} ata         Count the columns of A'A instead
     17    *
     18    * @return                     An array of size n of the column counts or null on error
     19    *
     20    * Reference: http://faculty.cse.tamu.edu/davis/publications.html
     21    */
     22   return function (a, parent, post, ata) {
     23     // check inputs
     24     if (!a || !parent || !post) {
     25       return null;
     26     } // a matrix arrays
     27 
     28 
     29     var asize = a._size; // rows and columns
     30 
     31     var m = asize[0];
     32     var n = asize[1]; // variables
     33 
     34     var i, j, k, J, p, p0, p1; // workspace size
     35 
     36     var s = 4 * n + (ata ? n + m + 1 : 0); // allocate workspace
     37 
     38     var w = []; // (s)
     39 
     40     var ancestor = 0; // first n entries
     41 
     42     var maxfirst = n; // next n entries
     43 
     44     var prevleaf = 2 * n; // next n entries
     45 
     46     var first = 3 * n; // next n entries
     47 
     48     var head = 4 * n; // next n + 1 entries (used when ata is true)
     49 
     50     var next = 5 * n + 1; // last entries in workspace
     51     // clear workspace w[0..s-1]
     52 
     53     for (k = 0; k < s; k++) {
     54       w[k] = -1;
     55     } // allocate result
     56 
     57 
     58     var colcount = []; // (n)
     59     // AT = A'
     60 
     61     var at = transpose(a); // at arrays
     62 
     63     var tindex = at._index;
     64     var tptr = at._ptr; // find w[first + j]
     65 
     66     for (k = 0; k < n; k++) {
     67       j = post[k]; // colcount[j]=1 if j is a leaf
     68 
     69       colcount[j] = w[first + j] === -1 ? 1 : 0;
     70 
     71       for (; j !== -1 && w[first + j] === -1; j = parent[j]) {
     72         w[first + j] = k;
     73       }
     74     } // initialize ata if needed
     75 
     76 
     77     if (ata) {
     78       // invert post
     79       for (k = 0; k < n; k++) {
     80         w[post[k]] = k;
     81       } // loop rows (columns in AT)
     82 
     83 
     84       for (i = 0; i < m; i++) {
     85         // values in column i of AT
     86         for (k = n, p0 = tptr[i], p1 = tptr[i + 1], p = p0; p < p1; p++) {
     87           k = Math.min(k, w[tindex[p]]);
     88         } // place row i in linked list k
     89 
     90 
     91         w[next + i] = w[head + k];
     92         w[head + k] = i;
     93       }
     94     } // each node in its own set
     95 
     96 
     97     for (i = 0; i < n; i++) {
     98       w[ancestor + i] = i;
     99     }
    100 
    101     for (k = 0; k < n; k++) {
    102       // j is the kth node in postordered etree
    103       j = post[k]; // check j is not a root
    104 
    105       if (parent[j] !== -1) {
    106         colcount[parent[j]]--;
    107       } // J=j for LL'=A case
    108 
    109 
    110       for (J = ata ? w[head + k] : j; J !== -1; J = ata ? w[next + J] : -1) {
    111         for (p = tptr[J]; p < tptr[J + 1]; p++) {
    112           i = tindex[p];
    113           var r = csLeaf(i, j, w, first, maxfirst, prevleaf, ancestor); // check A(i,j) is in skeleton
    114 
    115           if (r.jleaf >= 1) {
    116             colcount[j]++;
    117           } // check account for overlap in q
    118 
    119 
    120           if (r.jleaf === 2) {
    121             colcount[r.q]--;
    122           }
    123         }
    124       }
    125 
    126       if (parent[j] !== -1) {
    127         w[ancestor + j] = parent[j];
    128       }
    129     } // sum up colcount's of each child
    130 
    131 
    132     for (j = 0; j < n; j++) {
    133       if (parent[j] !== -1) {
    134         colcount[parent[j]] += colcount[j];
    135       }
    136     }
    137 
    138     return colcount;
    139   };
    140 });