simple-squiggle

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

csCounts.js (3581B)

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