1'use strict'; 2 3var CSynth, msgfixlog, springs, spearson, G, format, msgfix, sleep, log, numInstances, distxyz, col3, throws, 4Eigenvalues, VEC3, uniforms, geneOverrides, copyFrom, inworker, Worker, currentGenes, applyMatop, height, width, GO, framenum, glsl, S, array2Table, getstats,$; 5 6var setViewports, genedefs, mutate, slots, vps, setObjUniforms, renderObjsInner, mainvp, V, rot4toGenes, refmain, setAllLots, msgfixerrorlog, 7 clamp, U, addscript, makeDraggable, W, basescale; 8 9// get positions for key or ready made positions 10CSynth.pos = function(inputDef) { 11 const cc = CSynth.current; 12 let use; 13 // Allow for all sorts of inputs 14 if (inputDef === 'model') { // use current model as used by dynamics 15 throws('cannot use "model" to compute positions'); 16 } else if (typeof inputDef[0] === 'number') { // ready made distances (correct format assumed), Array.isArray(inputDef) 17 throws('cannot use distance array to compute positions'); 18 } else if (Array.isArray(inputDef) && 'z' in inputDef[0]) { // array of points 19 use = {coords: inputDef, reason: 'points as input'}; 20 } else if (!isNaN(parseInt(inputDef))) { // implicit index into xyzs 21 use = cc.xyzs[inputDef]; 22 } else if (typeof(inputDef) === 'string' && inputDef[0] === 'x') { // explict index into xyzs 23 const ssss = + inputDef.substring(1); 24 use = cc.xyzs[ssss]; 25 } else if (inputDef.isXyz) { // xyz structure 26 use = inputDef; 27 } else if (inputDef.startsWith('cur')) { // current positions 28 use = {coords: springs.getpos(), reason: 'current positions'}; 29 } else if (typeof(inputDef) === 'string' && inputDef[0] === 'c') { // explicit index into contacts 30 throws('cannot use contact numnber to compute positions'); 31 } else if (inputDef.isContact) { // contact structure 32 throws('cannot use contacts to compute positions'); 33 } else { 34 throw new Error('Invalid input to CSynth.dists', inputDef) 35 } 36 return use.coords; 37} 38 39/** compute dists for positions, return [distances, reason], dres can have array ready for distances */ 40CSynth.dists = function csynthdists(inputDef, statsres = CSynth.statsres, dres=undefined) { 41 const cc = CSynth.current; 42 let use; 43 44 // Allow for all sorts of inputs 45 // If this resolves to contacts, return associated distances 46 // If already array of numbers assume ready made distances 47 // If it resolves to points, compute distances 48 if (inputDef === 'model') { // use current model as used by dynamics 49 const sss = cc.selectedSpringSource; 50 const ncc = cc.contacts.length; 51 if (sss < ncc) use = cc.contacts[sss]; 52 else use = cc.xyzs[sss-ncc]; 53 } else if (typeof inputDef[0] === 'number') { // ready made distances (correct format assumed), Array.isArray(inputDef) fails on typed arrays 54 return [inputDef, 'dists as input']; 55 } else if (inputDef.length === 2 && typeof inputDef[0][0] === 'number' && typeof inputDef[1] === 'string') { // ready made dist pair 56 return inputDef; 57 } else if (Array.isArray(inputDef) && 'z' in inputDef[0]) { // array of points 58 use = {coords: inputDef, reason: 'points as input'}; 59 } else if (!isNaN(parseInt(inputDef))) { // implicit index into xyzs 60 use = cc.xyzs[inputDef]; 61 } else if (inputDef.startsWith('cur')) { // current positions 62 // special caching for current, handle and return 63 if (CSynth.dists.framenum === framenum && CSynth.dists.stepsSoFar === springs.getStepsSoFar()) 64 return CSynth.dists.framedists; 65 CSynth.dists.framenum = framenum; 66 CSynth.dists.stepsSoFar = springs.getStepsSoFar(); 67 return CSynth.dists.framedists = [gendists(springs.getpos()), 'current positions'] 68 } else if (typeof(inputDef) === 'string' && inputDef[0] === 'x') { // explict index into xyzs 69 const ssss = + inputDef.substring(1); 70 use = cc.xyzs[ssss]; 71 } else if (typeof(inputDef) === 'string' && inputDef.startsWith('cap')) { // explict index into contacts 72 const ssss = + inputDef.substring(3); 73 return [CSynth.captures[ssss].dists[0], 'capture' + ssss]; 74 } else if (inputDef.isXyz) { // xyz structure 75 use = inputDef; 76 } else if (inputDef.isContact) { // contact structure 77 use = inputDef; 78 } else { 79 throw new Error('Invalid input to CSynth.dists', inputDef) 80 } 81 let key = statsres + inputDef; 82 if (use.isContact) // check for changed details that effect contact distances 83 key += 'contactforcesc,patchval,pushapartforce,powBaseDist,pushapartpow,m_k,m_alph
83a,representativeContact,wrongfade' 84 .split(',').map(k=>GO[k]).join(','); 85 86 let result = use.dists && use.dists[key]; 87 if (result) return result; 88 89 if (use.isContact) { 90 result = [CSynth.contactToDistArray(use), use.shortname]; 91 } else { 92 result = [gendists(use.coords), use.shortname || use.reason]; 93 } 94 if (use.dists) { use.dists = {}; use.dists[key] = result; } // do not keep more than one 95 return result; 96 97 function gendists(s) { 98 // resample to correct resolution if necessary 99 if (s.length !== statsres) 100 s = CSynth.sample(s, statsres); 101 102 // and generate distances 103 const n = s.length; 104 const d = dres || new Float32Array(n * (n-1)/2); 105 let k = 0; 106 for (let i=0; i<n; i++) 107 for (let j=i+1; j<n; j++) 108 d[k++] = distxyz(s[i], s[j]); // do not use .distanceTo, may not be THREE vectors 109 return d; 110 } 111} 112 113/** compute correlation between two distance maps 114flags may be any (string) combination of rmse wrmse pearson spearman to choose stats to compute 115and nonormalize to prevent normalization 116dres1 and dres2 are optional arrays to hold distances for s1 and s2 117*/ 118CSynth.correl = function csynthcorrel(s1, s2, flags='all', statsres=undefined, dres1=undefined, dres2=undefined) { 119 flags = ' ' + flags + ' '; 120 if (flags.indexOf(' all ') !== -1) flags += ' rmse wrmse pearson spearman '; // bucket '; 121 let [d1, reason1] = CSynth.dists(s1, statsres, dres1); 122 let [d2, reason2] = CSynth.dists(s2, statsres, dres2); 123 if (d1.length !== d2.length) {msgfixlog('rmse', 'cannot compare rmse, wrong lengths', d1.length, d2.length)} 124 125 if (flags.indexOf(' filter1 ') !== -1) { 126 const mm = d1.reduce((c,x) => Math.max(c,x), 0); // why does d1.reduce(Math.max) not work 127 // const mm = Math.max(...d1); // blows it up 128 let d1a = [], d2a = []; 129 for (let i = 0; i < d1.length; i++) { 130 if (d1[i] !== mm) { 131 d1a.push(d1[i]); 132 d2a.push(d2[i]); 133 } 134 } 135 d1 = d1a; d2 = d2a; 136 } 137 138 const r = {}; 139 if (flags.indexOf( ' nonormalize ') === -1) { 140 const ss1 = d1.reduce((c,x) => c + x, 0); 141 const ss2 = d2.reduce((c,x) => c + x, 0); 142 r.scale = ss1/ss2; 143 if (CSynth.correlNormBase) { 144 //log('correl Normalize to base', ss1, ss2, CSynth.correlNormBase); 145 const scale1 = CSynth.correlNormBase/ss1; 146 const scale2 = CSynth.correlNormBase/ss2; 147 d1 = d1.map(x => x * scale1); 148 d2 = d2.map(x => x * scale2); 149 } else { 150 //log('correl Normalize to first', ss1, ss2); 151 d2 = d2.map(x => x * r.scale); 152 } 153 } 154 155 if (flags.indexOf( ' rmse ') !== -1) r.rmse = CSynth.rmse(d1, d2); 156 if (flags.indexOf( ' wrmse ') !== -1) r.wrmse = CSynth.wrmse(d1, d2); 157 if (flags.indexOf( ' pearson ') !== -1) r.pearson = spearson.correlation.pearson(d1, d2); // same as spearson.correl 158 if (flags.indexOf( ' oldpearson ') !== -1) r.oldpearson = spearson.correlation.oldpearson(d1, d2); 159 if (flags.indexOf( ' spearman ') !== -1) r.spearman = spearson.correlation.spearman(d1, d2, 'quick'); 160 if (flags.indexOf( ' oldspearman ') !== -1) r.oldpearman = spearson.correlation.spearman(d1, d2, true); 161 if (flags.indexOf( ' bucket ') !== -1) r.bucket = spearson.correlation.bucket(d1, d2); 162 // const sp2 = spearmanCorrelation(d1,d2); 163 return r; 164} 165CSynth.correlNormBase = 0; // 1000000; beware changes when using workers, <<< 166 167/** root mean square error */ 168CSynth.rmse = function csynthrmse(d1, d2) { 169 const n = d1.length; 170 let ss = 0; 171 if (d1[0].x === undefined) { 172 for (let i=0; i<n; i++) 173 ss += (d1[i]-d2[i]) ** 2; 174 } else { 175 for (let i=0; i<n; i++) { 176 const v1 = d1[i], v2 = d2[i]; 177 ss += (v1.x-v2.x) ** 2 + (v1.y-v2.y) ** 2 + (v1.z-v2.z) ** 2; 178 } 179 } 180 return Math.sqrt(ss / n); 181} 182 183/** weighted root mean square error */ 184CSynth.wrmse = function csynthwrmse (d1, d2) { 185 const n = d1.length; 186 let ss = 0, sw = 0; 187 for (let i=0; i<n; i++) { 188 const w = 1 / Math.max(d1[i], d2[i]); 189 ss += (w * (d1[i]-d2[i])) ** 2; 190 sw += w * w; // was 1, before that w ??? 191 } 192 return Math.sqrt(ss / sw); 193} 194 195 196// from https://stackoverflow.com/questions/15886527/javascript-library-for-pearson-and-or-spearman-correlations 197// used as cross-check.
198// This version is only valid where all contact/distance values are distinct. 199function spearmanCorrelation(p1, p2){ 200 const N=p1.length; 201 const order=[]; 202 let sum=0; 203 204 for(let i=0;i<N;i++){ 205 order.push([p1[i], p2[i]]); 206 } 207 208 order.sort(function(a,b){ 209 return a[0]-b[0] 210 }); 211 212 for(let i=0;i<N;i++){ 213 order[i].push(i+1); 214 } 215 216 order.sort(function(a,b){ 217 return a[1]-b[1] 218 }); 219 220 for(let i=0;i<N;i++){ 221 order[i].push(i+1); 222 } 223 for(let i=0;i<N;i++){ 224 sum+=Math.pow((order[i][2])-(order[i][3]), 2); 225 } 226 227 const r=1-(6*sum/(N*(N*N-1))); 228 229 return r; 230} 231 232/** get simple values for computing contact<=>dist relationship for current settings */ 233CSynth._getContactMultPow = function() { 234 let {contactforcesc, pushapartforce, powBaseDist, pushapartpow, m_k, m_alpha, representativeContact, patchval} = U; 235 let contactMult, contactPow, path; 236 // const buff = U.contactbuff.source.data.data; 237 if (contactforcesc !== 0) { 238 contactMult = contactforcesc / pushapartforce * (powBaseDist ** pushapartpow); 239 contactPow = (1 / (pushapartpow - 1)); 240 path = 'csy'; 241 } else { 242 contactMult = m_k ** (1/-m_alpha) / representativeContact; 243 contactPow = -m_alpha; 244 path = 'lor'; 245 } 246 return [contactMult, contactPow, path]; 247} 248 249/** compute equivalent wish dist for given contacts value with current settings */ 250CSynth.contact2Dist = CSynth.contact2dist = function(contact, opts = U) { 251 const [contactMult, contactPow] = CSynth._getContactMultPow(); 252 if (contact < 0) contact = U.patchval; 253 // else if (contact < minv) contact = minv; 254 const d = (contact * contactMult) ** contactPow; // regular distance 255 return d; 256} 257 258/** compute equivalent wish dist for given contacts value with current settings */ 259CSynth.dist2contact = function(dist, opts = U) { 260 const [contactMult, contactPow] = CSynth._getContactMultPow(); 261 return dist ** (1/contactPow) / contactMult; 262} 263 264 265 266/** compute complete array of target/wish dists for contacts, 267 * input may be number, contact or array 268 * if contactforcesc is nonzero assume CSynth implicit distances 269 * else assume v **-alpha style 270 * see CSynth.alignModels in csynth.js for some workings to deduce formula below 271 * and also compare with nval() in matrix.js 272 */ 273CSynth.contactToDistArray = function(c, alpha = G.m_alpha, k = G.m_k ) { 274 const cc = CSynth.current; 275 if (typeof(c) === 'number') c = cc.contacts[c]; 276 const minv = c.minv ? c.minv : 0.00001; // should usually be minv 277 278 if (c.isContact) { 279 CSynth.contactsToTexture(c); // in case not already done 280 c = c.textureData; 281 } 282 if (!(c instanceof Float32Array)) { 283 log('CSynth.contactToDistArray input not array', c); // not serious may be fixed in a frame or so 284 return; 285 } 286 287 288 const basen = Math.sqrt(c.length); // particles for real data 289 const n = CSynth.statsres || basen; //sample # particles 290 const ds = new Array(n * (n-1)/2); 291 // const pow = (a,b) => a**b; 292 293 // use overridden distances, and cache them for efficiency 294 const {contactforcesc, patchval, pushapartforce, powBaseDist, pushapartpow, wrongfade} = GO; 295 const representativeContact = cc.representativeContact; 296 const [contactMult, contactPow, path] = CSynth._getContactMultPow(); 297 if (n < 4) log('contactToDistArray', path, contactMult, contactPow); // debug 298 299 if (n === basen) { // original quick version 300 let p = 0; 301 for (let i=0; i<n; i++) { 302 for (let j=i+1; j<n; j++) { 303 let contact = c[i*n + j]; // / representativeContact; 304 if (contact < 0) contact = patchval; 305 else if (contact < minv) contact = minv; 306 const d = (contact * contactMult) ** contactPow; // regular distance 307 ds[p++] = d; 308 } 309 } 310 } else { // version allowing sample 311 const ss = CSynth.sample(new Array(basen).fill(0).map((x,i) => i), n); // sample inddices 312 let p = 0; 313 for (let i=0; i<n; i++) { 314 const ii = ss[i]; 315 for (let j=i+1; j<n; j++) { 316 let contact = c[ii*basen + ss[j]]; // / representativeContact; 317 if (contact < 0) contact = patchval; 318 else if (contact < minv) contact = minv; 319 const d = (contact * contactMult) ** contactPow; // regular distance 320 ds[p++] = d; 321 } 322 } 323 } 324 return ds; 325} 326 327/** get and cache mean distances --- WARNING no check on model parameters changing cached values */ 328CSynth.meandists = function(contact) { 329 const dists = CSynth.contactToDistArray(contact); 330 const n = contact.numInstances; 331 const c = contact.textureData; 332 let sv=0, swd = 0, sd = 0, p = 0; 333 for (let i=0; i<n; i++) { // complicated as c is square, dists triangular 334 for (let j=i+1; j<n; j++) { 335 const d = dists[p]; // MAY GET OUT OF DATE 336 const v = c[i*n + j]; 337 if (v >= 0) { 338 sd += d; 339 swd += d * v; 340 sv += v; 341 } 342 p++; 343 } 344 }
345 const alpha = -1 / (G.pushapartpow - 1); 346 const r = { 347 wmeandist: swd / sv, meandist: sd / p, meanndv: sv / p, 348 wmeandistI: sv / swd, meandistI: p / sd, 349 wmeandistC: (swd / sv) ** (-1/alpha), meandistC: (sd / p) ** (-1/alpha) 350 } ; 351 copyFrom(contact, r); // << may get out of date 352 return r; 353} 354 355// compare the weights, assume two IF inputs, for debug 356var yaml; 357CSynth.compareWeights = function() { 358 const cc = CSynth.current, ccc0 = cc.contacts[0], ccc1 = cc.contacts[1]; 359 CSynth.meandists(ccc0); 360 CSynth.meandists(ccc1); 361 const r = {}; 362 CSynth.representativeSources.forEach(s => { 363 let a = ccc0[s]; 364 let b = ccc1[s]; 365 366 const rat = b/a; 367 r[s] = {s,a,b,rat} 368 }); 369 log(yaml.safeDump(r, {skipInvalid: true})); 370 return r; 371} 372 373/** show all stats relative to current model, older version without workers */ 374CSynth.allstatsSlow = async function allstatsSlow(statsres) { 375 msgfix.show('correlations'); 376 msgfixlog('correlations', 'pending'); 377 const cc = CSynth.current; 378 let [refd, refreason] = CSynth.dists('model', statsres); 379 const msg = ['vs ' + refreason + '. "Y" to recompute.']; 380 // msg.push('current: ' + format(CSynth.correl(refd, 'current'), 6)); 381 for (let i=0; i < cc.xyzs.length; i++) { 382 msgfix('correlations', msg.join('<br>'), '<br>...'); 383 await sleep(0); 384 const xyz = cc.xyzs[i]; 385 const mm = log(xyz.shortname + ': ' + format(CSynth.correl(refd, i), 6)); 386 msg.push(mm); 387 } 388 msgfixlog('correlations', msg.join('<br>'), '<br>done'); 389} 390CSynth.workers = []; 391CSynth.statsFields = Object.keys({scale: 1.196170, rmse: 6.396037, wrmse: 0.238274, pearson: 0.721054, spearman: 0.929094}); 392 393/** show all stats relative to current model using workers */ 394CSynth.allstats = async function allstats(statsres = CSynth.statsres) { 395 msgfix.show('correlations'); 396 msgfixlog('correlations', 'pending'); 397 const cc = CSynth.current; 398 let [refd, refreason] = CSynth.dists('model', statsres); 399 const msg = ['vs ' + refreason + '. "Y" to recompute.', format(CSynth.statsFields, 8)]; 400 // msg.push('current: ' + format(CSynth.correl(refd, 'current'), 6)); 401 for (let i=0; i < cc.xyzs.length; i++) { 402 const xyz = cc.xyzs[i]; 403 let w = CSynth.workers[i]; 404 if (!w) { 405 w = CSynth.workers[i] = new Worker('csynthWorker.js'); 406 } 407 w.onmessage = e => { 408 msg[i+2] = xyz.shortname + ': ' + format(Object.values(e.data[1]), 6); 409 msgfix('correlations', msg.join('<br>')); 410 }; 411 w.postMessage(['correlCD', refd, xyz.coords, 'all', statsres]); 412 msg.push(xyz.shortname + ': ...'); 413 await sleep(0); 414 } 415} 416 417CSynth.clearWorkers = function() { 418 CSynth.workers.forEach(w => w.postMessage(['close'])); 419 CSynth.workers = []; 420} 421 422 423CSynth.curstats = function(statsres) { 424 msgfix.show('cstats'); 425 let [refd, refreason] = CSynth.dists('model', statsres); 426 427 msgfix('cstats', 'current vs ' + refreason + '<br>' + format(CSynth.correl('c0', 'current'), 6)); 428} 429 430/** run n tests and collect stats against distance derived from contact map 0 (c0) */ 431CSynth.statsn = function(opts) { 432 let {n, w, base} = Object.assign({n: 5, w: 2000, base: 'c0'}, opts); 433 const basedists = CSynth.dists(base); 434 const r = []; 435 for (let i = 1; i <= n; i++) { 436 CSynth.randpos(undefined, i); 437 springs.step(w); 438 r.push(`csynth${i}: ${format(Object.values(CSynth.correl(basedists, 'current')), 6)}`); 439 msgfix('statsn', 'csynth against ref IF<br>', r.join('<br>')); 440 } 441} 442 443/** clear the captured data */ 444CSynth.clearCaptures = function() { 445 return CSynth.captures = []; // {sr: [], dists: [], positions: []}; 446}; 447CSynth.captures = CSynth.clearCaptures(); 448 449/** capture data into a capture set and show stats so far */ 450CSynth.capture = function({set = CSynth.captures, correl = false} = {}) { 451 // const {sr, dists, positions} = set; 452 const i = set.length; 453 const sr = []; 454 const positions = springs.getpos(); 455 const dists = CSynth.dists('current'); // collect the distances from the simulated positions 456 const genes = Object.assign({}, G); 457 set.push({sr, positions, dists, genes}) 458 if (correl) { 459 for (let j = 0; j < i; j++) { // compare run i with previous runs j 460 sr[j] = CSynth.correl(dists, set[j].dists, correl); // perform correlation 461 } 462 CSynth.displayCaptureStats(set); 463 } 464} 465 466/** goto a particular capture, w.i.p. */ 467CSynth.gotoCapture = function(n = 0, set = CSynth.captures) { 468 // const {sr, dists, positions} = set; 469 const s = set[n]; 470 if (!s) return msgfixlog('cannot go to capture item ', n); 471 springs.setpos(s.positions);
472 Object.assign(currentGenes, s.genes); 473 setObjUniforms(currentGenes, uniforms); 474} 475 476/** display stats information about a set (probably captures) */ 477CSynth.displayCaptureStats = function(set = CSynth.captures, name = 'captures', msg = 'stats for captured data' ) { 478 // const {sr, dists, positions} = CSynth.captures; 479 const r = []; 480 if (set.length < 2) return 'not enough capture data'; 481 const tto = ['', '']; 482 tto.push(...Object.keys(set[1].sr[0])) 483 r.push(tto) 484 for (let i = 0; i < set.length; i++) { 485 for (let j = 0; j < i; j++) { 486 const tt = [i,j] 487 tt.push(...Object.values(set[i].sr[j]).map(c => c.toFixed(3))) 488 489 // let tt = `csynth ${i} ${j} ${format(Object.values(set[i].sr[j]), 6)}` 490 //if (j === 0) tt = `'<green>${tt}</green>` 491 r.push(tt); // and format 492 if (j === 0) r.push( r.pop().map(c => '<green>' + c + '</green>')); 493 } 494 } 495 msgfix(name, msg + array2Table(r)); 496 //msgfix(name, msg + '<br>' + r.join('<br>')); 497} 498 499/** generate an Array from a definition object (def) and options 500 * def is an object; 501 * each field value is either 502 * an object {from, to, [step]}; the field name will take appropriate values (to inclusive) 503 * an array; the field name will take values from each element of the array 504 * other: the field name will take the given value 505 * The result is an array of objects, with the same keys as the definition object 506 * but with each element of the array taking one of the candidate values. 507 * 508 * To consider: make genArray a generator 509 */ 510CSynth.genArray = function(def, popts = {}) { 511 const opts = Object.assign({ 512 keys: Object.keys(def), // the object keys still to iterate 513 vals: {}, // the object to hold each result in turn 514 res: [] // the array to accumulate the result 515 }, popts); 516 const {keys, vals, fun, res} = opts; 517 const key = keys[0]; 518 let litem = def[key]; 519 520 const xopts = Object.assign({}, opts); 521 xopts.clear = false; 522 xopts.keys = keys.slice(1); 523 if (litem.from !== undefined) { // 524 const x = []; 525 for(let v=litem.from; v <= litem.to; v += litem.step || 1) x.push(v); 526 litem = x; 527 } 528 if (!Array.isArray(litem)) litem = [litem]; 529 for (const v of litem) { 530 xopts.vals[key] = v; 531 if (keys.length === 1) { 532 res.push(Object.assign({}, xopts.vals)); 533 } else { 534 CSynth.genArray(def, xopts); 535 } 536 } 537 return res; 538} 539 540var Maestro; 541// capture a set of stats given ab 542CSynth.capture8 = function(w = 5000, opts={}) { 543 CSynth.captures = CSynth.clearCaptures(); 544 const arr = CSynth.genArray({contactforce: [100,400], pushapartpow: [0, -1], springpow: [0, -2] }); 545 for (const item of arr) { 546 Object.assign(currentGenes, item); 547 // Maestro.trigger('prespringstep'); 548 setObjUniforms(currentGenes, uniforms); 549 CSynth.randpos(undefined, 2); // generate random positions from seed i 550 springs.step(w); // run w simulation steps 551 if (opts.autoscale) CSynth.autoscale(); 552 // framenum++; 553 // await S.frame(); 554 CSynth.capture(); 555 } 556} 557 558var GX, setBackgroundColor, vpborder, onframe; 559// reminder of contruction of paperS1X 560CSynth.paperS1X = function() { 561 // http://localhost:8800/csynth.html?startscript=C75X1200/load_data.js 562 V.gui.visible = false; 563 GX.setValue('matrix/visible', false); 564 setBackgroundColor(1); 565 vpborder = 0; 566 CSynth.capture8(5000, {autoscale: true}); 567 onframe(CSynth.draw8); 568} 569 570/** run n simulations of w simulation steps against different seeds, 571 * and collect stats against each other */ 572CSynth.comparen = function(n = 5, w = 10000) { 573 const sr = [], r = []; 574 const dists = []; // collect the distance maps 575 for (let i = 1; i <= n; i++) { 576 sr[i] = []; 577 CSynth.randpos(undefined, i); // generate random positions from seed i 578 springs.step(w); // run w simulation steps 579 dists[i] = CSynth.dists('current'); // collect the distances from the simulated positions 580 for (let j = 1; j < i; j++) { // compare run i with previous runs j 581 sr[i][j] = CSynth.correl(dists[i], dists[j]); // perform correlation 582 r.push(`csynth ${i} ${j} ${format(Object.values(sr[i][j]), 6)}`); // and format 583 } 584 msgfix('statsij', `CSynth against each other <br>${r.join('<br>')}`); // display output 585 } 586 return sr; // return structured output for further analysis (if wanted) 587} 588 589/** compare n LorDG models */ 590CSynth.comparex = function(n = 5) { 591 const sr = [], r = []; 592 const dists = []; // collect the distance maps 593 for (let i = 1; i <= n; i++) { 594 sr[i] = []; 595 dists[i] = CSynth.dists('x' + (i-1)); // collect the distances from the simulated positions 596 for (let j = 1; j < i; j++) { // compare run i with previous runs j 597 sr[i][j] = CSynth.correl(dists[i], dists[j]); // perform correlation 598 r.push(`LorDG ${i} ${j} ${format(Object.values(sr[i][j]), 6)}`); // and format 599 } 600 msgfix('statslorij', `LorDG against each other <br>${r.join('<br>')}`); // display output 601 } 602 return sr; // return structured output for further analysis (if wanted) 603} 604 605CSynth.maxRealtimeStatsSize = 1000; 606 607/** cache and compute referenceDistances if necessary */ 608CSynth.referenceDistances = function refdist() { 609 const cc = CSynth.current; 610 if (!cc || !cc.contacts || !cc.contacts[0]) return; 611 const ccc0 = cc.contacts[0];
612 let check=[]; 613 if (!CSynth._referenceDistances || cc.contacts[0].referenceDistances !== CSynth._referenceDistances 614 || G.m_alpha !== check[0] || G.m_k !== check[1]) { 615 cc.contacts[0].referenceDistances = CSynth._referenceDistances = CSynth.contactToDistArray(0, G.m_alpha, G.m_k); 616 check = [G.m_alpha, G.m_k]; 617 } 618 return CSynth._referenceDistances; 619} 620 621CSynth.startDynamicStats = function(flag = 'all') { 622 const cc = CSynth.current; 623 624 for (let i=0; i < cc.contacts.length; i++) 625 msgfix('!stats cur c' + i, ()=>CSynth.correl('c' + i, 'current', flag)); 626 for (let i=0; i < cc.xyzs.length; i++) 627 msgfix('!stats cur x' + i, ()=>CSynth.correl('x' + i, 'current', flag)); 628} 629 630CSynth.fixedStats = function(flag = 'all') { 631 const cc = CSynth.current; 632 let cl = cc.contacts.length, tl = cl + cc.xyzs.length; 633 for (let i = 0; i < tl; i++) { 634 const ii = i < cl ? 'c' + i : 'x' + (i-cl); 635 for (let j = i+1; j < tl; j++) { 636 const jj = j < cl ? 'c' + j : 'x' + (j-cl); 637 msgfix('stats ' + ii + ' ' + jj, CSynth.correl(ii,jj)); 638 } 639 } 640} 641 642var gpuStats; 643// get 'volume' of array of points, current points by default 644// return eigenvalues and 'volume' 645CSynth.stats = function(p = springs.getpos()) { 646 // let x=0, y=0, z=0, xx=0, yy=0, zz=0, xy=0, yz=0, zx=0; 647 let max = VEC3(-Number.MAX_VALUE, -Number.MAX_VALUE, -Number.MAX_VALUE); 648 let min = VEC3(Number.MAX_VALUE, Number.MAX_VALUE, Number.MAX_VALUE); 649 const n = p.length; 650 651 // compute raw sums 652 let sx=0, sy=0, sz=0, sxx=0, syy=0, szz=0, sxy=0, syz=0, szx=0; 653 for (let i=0; i < n; i++) { 654 const v = p[i]; 655 sx += v.x; 656 sy += v.y; 657 sz += v.z; 658 max = max.max(v); 659 min = min.min(v); 660 sxx += v.x * v.x; 661 syy += v.y * v.y; 662 szz += v.z * v.z; 663 sxy += v.x * v.y; 664 syz += v.y * v.z; 665 szx += v.z * v.x; 666 } 667 const r =CSynth._stats2({sx, sy, sz, sxx, syy, szz, sxy, syz, szx, n}) 668 r.min = min; r.max = max; 669 r.rx = max.x - min.x; r.ry = max.y - min.y; r.rz = max.z - min.z; 670 r.rmax = Math.max(r.rx, r.ry, r.rz); 671 return r; 672} 673 674/** shared function for anayzing cumulative statistics */ 675CSynth._stats2 = function({sx, sy, sz, sxx, syy, szz, sxy, syz, szx, n}) { 676 // compare raw sum values with equivalent values from gpuStats 677 // debug 678 //const r = gpuStats(); 679 //const [bx, by, bz, bxx, byy, bzz, byz, bzx, bxy] = 680 // [r[0][0], r[1][0], r[2][0], r[0][1], r[1][1], r[2][1], r[0][2], r[1][2], r[2][2]]; 681 //log('x,y,z', sx, bx, sy,by, sz,bz); 682 //log(sxx, bxx, syy,byy, szz,bzz, sxy,bxy, syz,byz, szx,bzx); 683 //log(f(sxx- bxx)+'', f(syy-byy)+'', f(szz-bzz)+'', f(sxy-bxy)+'', f(syz-byz)+'', f(szx-bzx)+''); 684 685 // 'compensated' sums for more convenient stats 686 const xx = sxx - sx*sx/n; 687 const yy = syy - sy*sy/n; 688 const zz = szz - sz*sz/n; 689 const xy = sxy - sx*sy/n; 690 const yz = syz - sy*sz/n; 691 const zx = szx - sz*sx/n; 692 const x = sx/n; 693 const y = sy/n; 694 const z = sz/n; 695 696 // now derive the various stats from the raw sum figures 697 const mat3 = [xx/n, xy/n, zx/n, xy/n, yy/n, yz/n, zx/n, yz/n, zz/n]; 698 const trace = xx/n + yy/n + zz/n; 699 let eigv = Eigenvalues.eigenvalues(mat3); 700 const trace2 = eigv[0] + eigv[1] + eigv[2]; // for crosscheck 701 // log('traces', trace, trace2, f(trace2-trace)+''); 702 // const hessenberg = Eigenvalues.hessenberg(mat3); 703 const centroid = VEC3(x, y, z); 704 let radii = eigv.map(v=> Math.sqrt(v)); 705 const volume = 4/3 * Math.PI * radii[0] * radii[1] * radii[2]; 706 707 708 // compuute eigenvectors using Cayley–Hamilton theorem 709 // https://en.wikipedia.org/wiki/Eigenvalue_algorithm#Hessenberg_and_tridiagonal_matrices 710 // https://en.wikipedia.org/wiki/Cayley%E2%80%93Hamilton_theorem 711 const {full, d1,d2,d3,s1,s2,s3} = CSynth.stats.mats; // CSynth.stats.mats contains useful working objects 712 full.elements = mat3; 713 d1.copy(full); d1.elements[0] -= eigv[0]; d1.elements[4] -= eigv[0]; d1.elements[8] -= eigv[0]; 714 d2.copy(full); d2.elements[0] -= eigv[1]; d2.elements[4] -= eigv[1]; d2.elements[8] -= eigv[1]; 715 d3.copy(full); d3.elements[0] -= eigv[2]; d3.elements[4] -= eigv[2]; d3.elements[8] -= eigv[2]; 716 s1.multiplyMatrices(d2,d3); 717 s2.multiplyMatrices(d3,d1); 718 s3.multiplyMatrices(d1,d2); 719 let eigenvectors = [ 720 VEC3(s1.elements).normalize().multiplyScalar(radii[0]), 721 VEC3(s2.elements).normalize().multiplyScalar(radii[1]), 722 VEC3(s3.elements).normalize().multiplyScalar(radii[2]) 723 ] 724 725 // order eigenvalues, eigenvectors and radii 726 const oo = eigv.map((v,i) => [v,i]).sort((a,b)=>b[0] - a[0]).map(v=>v[1]); 727 eigv = oo.map(v=> eigv[v]); 728 eigenvectors = oo.map(v=> eigenvectors[v]); 729 radii = oo.map(v=> radii[v]); 730 731 732 return { eigenvalues: eigv, radii, volume, centroid, trace, eigenvectors} ; 733} 734var THREE, ColorKeywords; 735if (!inworker) 736CSynth.stats.mats = { 737 full: new THREE.Matrix3(), d1: new THREE.Matrix3(), d2: new THREE.Matrix3(), d3: new THREE.Matrix3(), 738 s1: new THREE.Matrix3(), s2: new THREE.Matrix3(), s3: new THREE.Matrix3() 739}; 740 741/** show the eigenvectors and centroid, d is a stats object */ 742CSynth.showEigen = function(doOrient, d = CSynth.stats()) { 743 const lineGeometry = new THREE.BufferGeometry(), positions=[], colors=[]; 744 const vert = /*glsl*/` 745 precision highp float; 746 ${CSynth.CommonShaderCode()} 747 // attribute vec3 position; 748 attribute vec3 color; 749 varying vec3 vcol; 750 // const float xx = 1.; 751 void main() { 752 vec4 p= pposToWorld(vec4(position, 1)); 753 p.w = 1.; 754 gl_Position = logdepth(projectionMatrix * viewMatrix * p); // model part already in partposWSorld 755 vcol = color; 756 } 757 `; 758 const frag = /*glsl*/` 759 precision highp float; 760 varying vec3 vcol; 761 void main() { 762 gl_FragColor = vec4(vcol, 1); 763 } 764 `; 765 // const lineMaterial = new THREE.LineBasicMaterial( { linewidth:1, color: 0xffffff, vertexColors: THREE.VertexColors } ); 766 const lineMaterial = new THREE.RawShaderMaterial({ 767 fragmentShader: frag, vertexShader: vert, uniforms: uniforms 768 }); 769 770 if (CSynth.eigenLine) V.rawscene.remove(CSynth.eigenLine); 771 772 const line = CSynth.eigenLine = new THREE.LineSegments(lineGeometry, lineMaterial); 773 V.rawscene.add(line);
774 const kcols = ['red','green','blue'].map(c=>new THREE.Color(c)); 775 const pushv = (v) => positions.push(v.x, v.y, v.z); 776 const pushc = (v) => colors.push(v.r, v.g, v.b); 777 for (let i = 0; i < 3; i++) { 778 pushv(d.centroid.clone().sub(d.eigenvectors[i])); 779 pushv(d.centroid.clone().add(d.eigenvectors[i])); 780 pushc(kcols[i]); 781 pushc(kcols[i]); 782 } 783 lineGeometry.setAttribute('position', new THREE.BufferAttribute(new Float32Array(positions), 3)); 784 lineGeometry.setAttribute('color', new THREE.BufferAttribute(new Float32Array(colors), 3)); 785 786 if (doOrient) 787 CSynth.statsOrient(d, {angled: doOrient === 'angled', locate: false, scale: false}); 788 789} 790 791// not really CSynth ... 792CSynth.statsOrient = function statsOrient(stats = skelstats(), {angled = false, locate = true, scale = true } = {}) { 793 const e = stats.eigenvectors; 794 for (let i=0; i<2; i++) e[i].normalize(); 795 const t = currentGenes._rot4_ele; 796 e[2].crossVectors(e[0], e[1]); // third eigenvector sometimes give -ve determinant 797 798 const v = stats.eigenvalues.map(x=>x**0.5); 799 const cen = stats.centroid.multiplyScalar(locate ? 1/basescale : 0); 800 const vsize = Math.max(v[0], v[1] * width / height); 801 if (scale) G._uScale = basescale/vsize/1.5; 802 803 const x = e[0], y = e[1], z = e[2]; 804 G._rot4_ele = [x.x, x.y, x.z, cen.x, y.x, y.y, y.z, cen.y, z.x, z.y, z.z, cen.z, 0,0,0,1]; 805 806 if (angled) 807 applyMatop(0,1, Math.atan2(height, width)); 808 else 809 rot4toGenes(); 810 refmain(); 811} 812 813/** orient with centroid to middle,particle # px to the right and py to the top */ 814CSynth.orient = function(px = CSynth.orient.right, py = CSynth.orient.up) { 815 if (typeof px === 'string') { 816 const p = CSynth.picks['g-ribbon']; 817 if (p < 0) return msgfix('orient', 'no particle selected for RIGHT/UP'); 818 CSynth.orient[px] = p.partid; 819 msgfix('orient', `particles are ${CSynth.orient.right} RIGHT and ${CSynth.orient.up} UP.`) 820 return; 821 } 822 if (px === undefined) px = Math.round(numInstances/2); 823 if (py === undefined) py = numInstances - 1; 824 825 const p = springs.getpos(); 826 const c = p.reduce((cc,v) => cc.add(v), VEC3()).multiplyScalar(1 / numInstances); 827 const x = p[px].sub(c).normalize(); 828 const y1 = c.clone().sub(p[py]); 829 const z = y1.cross(x).normalize(); 830 const y = z.clone().cross(x); 831 G._rot4_ele = [x.x, x.y, x.z, 0, y.x, y.y, y.z, 0, z.x, z.y, z.z, 0, 0,0,0,1]; 832} 833 834/** tilt so x axis aligns with screen diagonal */ 835CSynth.tilt = function() {applyMatop(0,1, Math.atan2(height, width));} 836 837 838/** work out forces for different force styles at given distance with current settings: 839This uses G values without override, so computes forces whether or not they are enabled. 840If no length is provided, use representativeDistance = wish length from Lorentz model based on representativeContact. 841If no contact provided, use representativeContact. 842+ve attraction force for long distance */ 843CSynth.forces = function(len, contact = CSynth.current.representativeContact) { 844 845 const cc = CSynth.current; 846 const {pow, sqrt, max} = Math; 847 const {representativeContact} = cc; 848 // const g = {}; copyFrom(g, G); copyFrom(g, geneOverrides); 849 const {pushapartforce, powBaseDist, pushapartpow, xyzforce, m_k, m_force, m_alpha, m_c, contactthreshold} = G; 850 const r = {} 851 const rc = cc.representativeContact; 852 853 const targlen = m_k * pow(contact/representativeContact, -m_alpha) 854 855 // lor forces 856 const rcontact = contact/representativeContact; 857 const d = m_k * pow(rcontact, -m_alpha); // target distance 858 if (len === undefined) len = d; 859 const dd = d - len; 860 const dem = m_c * m_c + dd*dd; 861 r.lorforce = m_force * rcontact * -2 * m_c * m_c * dd / (dem*dem); 862 863 // csynth forces, note r._ values just for debug 864 const ccontact = max(0, contact - contactthreshold); 865 const contactforcesc = G.contactforce * 1e-6 / CSynth.current.representativeContact; 866 const contactforce = r._contactforce = contactforcesc * ccontact * len; 867 868 const globalpushapart = r._globalpushapart = pushapartforce * pow(len/powBaseDist, pushapartpow); 869 r.csyforce = contactforce - globalpushapart; 870 871 let dlen = len - targlen; 872 // sqrt a bit arbitrary here, but helped ensure correct refolding of e.g. different lor examples 873 r.xyzforce = dlen * xyzforce / sqrt(targlen); 874 875 return r; 876} 877 878/** work out forces at representativeDistance to check they all stabalize at 0,
879and at k (default 1.1) times to find effective force from different forces */ 880CSynth.forcetest = function(k = 1.1, errlevel = 1e-10, len = G.m_k) { 881 const r0 = CSynth.forces(len); 882 // log(r0); 883 const abs = Math.abs; 884 for (let f in r0) 885 if (f[0] !== '_' && Math.abs(r0[f]) > errlevel) 886 console.error(`Force ${f} not balanced ${r0[f]} at representative length ${len}.`); 887 const r = CSynth.forces(len*k); 888 // log(r); 889 return r; 890} 891 892/** align force strengths to put all into same range; e.g. same value for reference distance * 1.1 */ 893CSynth.alignForces = function(use = 'lor') { 894 895 for (const gn in CSynth.defvals) G[gn] = G[gn] || CSynth.defvals[gn]; // set 0 values to defval value, this should render then necxt few lines unnecessary 896 const v = G.contactforce * G.pushapartforce * G.m_force * G.xyzforce; 897 if (v === 0 || isNaN(v)) { 898 console.error('cannot align forces while one is 0 or NaN'); 899 return; 900 } 901 const r = CSynth.forcetest(1.01, 1e-15); // list of current changes 902 switch (use.replace('force','')) { 903 case 'lor': 904 G.contactforce *= r.lorforce / r.csyforce; 905 G.pushapartforce *= r.lorforce / r.csyforce; 906 G.xyzforce *= r.lorforce / r.xyzforce; 907 break; 908 case 'csy': 909 case 'old': 910 G.m_force *= r.csyforce / r.lorforce; 911 G.xyzforce *= r.csyforce / r.xyzforce; 912 break; 913 case 'xyz': 914 G.contactforce *= r.xyzforce / r.csyforce; 915 G.pushapartforce *= r.xyzforce / r.csyforce; 916 G.m_force *= r.xyzforce / r.lorforce; 917 break; 918 default: 919 console.error('bad parameter to CSynth.alignForces', use); 920 } 921 922} 923 924/** align models to given model (? that is, align Lorenz m_ values with our contactforce etc values */ 925CSynth.alignModels = function(type = 'auto') { 926 // Lorentz, wish distance = m_k / if ** m_alpha 927 // CSynth, implicit wish distance (if * contactforce/pushapartforce) ** (1/(pushapartpow-1)) 928 // = (if * contactforce/pushapartforce) ** -alpha 929 // = (contactforce/pushapartforce) ** -alpha / if ** m_alpha 930 /***** 931 Workings for CSynth implicit wish distance (allowing for ) 932 933 contactforcesc = contactforce * 1e-6 / representativeContact 934 f = contactforcesc * contact * len 935 f = pushapartforce * pow(len/powBaseDist, pushapartpow) 936 = pushapartforce / pow(powBaseDist, pushapartpow) * pow(len, pushapartpow) 937 938 contactforcesc * contact / pushapartforce * pow(powBaseDist, pushapartpow) = pow(len, pushapartpow-1) 939 len = pow(contactforcesc * contact / pushapartforce * pow(powBaseDist, pushapartpow), 1/(pushapartpow-1)) 940 */ 941 942 const oal = CSynth.lastAlign; 943 // const nal = [G.pushapartpow, G.contactforce, G.m_alpha, G.m_k]; 944 945 const types = CSynth.forcetypes; // defined CSynth.switchSpringSettings 946 947 let forceChanged = true; 948 if (type === 'auto') { 949 // track the types and find the last one to have changed 950 //if (types.every(x => oal[x[1]] === G[x[0]] )) 951 // return; // nothing changed 952 types.forEach( (p,i) => {if (oal[i][1] !== G[p[0]]) type = p[1];} ); 953 if (type === 'auto') return; // nothing changed 954 forceChanged = oal.some( x => x[0].endsWith('force') && x[1] !== G[x[0]]) 955 } 956 957 958 const x = {}; 959 // both models want to be 'pure' of extra forces; NO leave to the calling environment 960 // x.backboneforce = 0; 961 // x.pushapartlocalforce = 0; 962 963 // -alpha = 1 / (pushapartpow-1) 964 if (type === 'lor') { 965 x.pushapartpow = 1 - 1 / G.m_alpha; 966 x.contactforce = 1e6 * G.pushapartforce / G.powBaseDist * (G.m_k / G.powBaseDist) ** (-1/G.m_alpha) 967 //if (CSynth.springSettings.contact) 968 // copyFrom(CSynth.springSettings.contact, x); 969 //else 970 // copyFrom(G, x); 971 } else if (type === 'csy') { 972 x.m_alpha = -1 / (G.pushapartpow - 1); 973 x.m_k = G.powBaseDist * Math.pow(G.contactforce/1e6/G.pushapartforce * G.powBaseDist, -x.m_alph
973a); 974 975 // // check debug, at distance m_k what are the csy forces 976 // copyFrom(G, x); 977 // var len = x.m_k 978 // var repv = CSynth.current.representativeContact 979 // var backboneScale = 0.00001 // ??? 980 // var contactforcesc = G.contactforce * 1e-6 / repv; // === G.contactforcesc 981 // var cforce = repv * contactforcesc * len 982 // var pforce = G.pushapartforce * Math.pow(Math.max(len,backboneScale)/G.powBaseDist, G.pushapartpow) 983 // log('forces len', CSynth.forces(len)) 984 // log('forces len', CSynth.forces(len)) 985 // var xx = {cforce, pforce} 986 // log('test', xx) 987 // xx=xx 988 // //if (CSynth.springSettings.contact Lor) 989 // // copyFrom(CSynth.springSettings.contact Lor, x); 990 // //else 991 // // copyFrom(G, x); 992 } else if (type === 'xyz') { // 993 994 } else { 995 console.error(`Unexpected type '${type}' in CSynth.alignModels.`); 996 } 997 copyFrom(G, x); 998 999 CSynth.alignForces(forceChanged ? type : 'old'); 1000 1001 // G.m_k = G.powBaseDist * Math.pow(G.contactforce/1e6/G.pushapartforce * G.powBaseDist, -G.m_alpha) 1002 1003 // make x apply to both current and pending (if any) 1004 1005 // const q = {}; // for debug logging 1006 // q.m_alpha = -1 / (G.pushapartpow - 1) 1007 // q.m_k = G.powBaseDist * Math.pow(G.contactforce/1e6/G.pushapartforce * G.powBaseDist, -G.m_alpha) 1008 // log (`diffs, m_alpha ${q.m_alpha} ${G.m_alpha} ... ${q.m_k} ${G.m_k}`); 1009 CSynth.lastAlign = types.map(p => [p[0], G[p[0]]]); 1010} 1011CSynth.lastAlign = []; 1012 1013CSynth.saveRuns = async function CSynth_compareRuns({nruns=10, fid, time1=4000, time2=1000}) { 1014 if (!fid) fid = '>csynth_' + (new Date()).toJSON().replace(/:/g, ".") + '_'; 1015 G.stepsPerStep = 40; 1016 for (let run=0; run<nruns; run++) { 1017 CSynth.randpos(6); 1018 G.damp=1; 1019 await sleep(time1) 1020 G.damp=0.99; 1021 await sleep(time2) 1022 springs.save(fid + ('00000' + run).substr(-6) + '.xyz'); 1023 msgfix('runs', run+1, nruns); 1024 if (CSynth.compareRuns.break) { 1025 CSynth.compareRuns.break = false; 1026 break; 1027 } 1028 } 1029} 1030 1031 1032CSynth.compareRuns = async function CSynth_compareRuns(nruns=10) { 1033 G.stepsPerStep = 40; 1034 let runs = CSynth.compareRuns.runs = []; 1035 for (let run=0; run<nruns; run++) { 1036 CSynth.randpos(6); 1037 G.damp=1; 1038 await sleep(4000) 1039 G.damp=0.99; 1040 await sleep(1000) 1041 const a = {coords: springs.getpos(), numInstances} 1042 // CSynth.xyzToTexture(a) 1043 runs.push(a) 1044 msgfix('runs', run, nruns); 1045 if (CSynth.compareRuns.break) { 1046 CSynth.compareRuns.break = false; 1047 break; 1048 } 1049 } 1050 return CSynth.statsForRuns(runs); 1051} 1052 1053CSynth.statsForRuns = function CSynth_statsForRuns(runs = CSynth.compareRuns.runs) { 1054 const stats = new Array(numInstances) 1055 const nruns = runs.length; 1056 const pp = new Array(nruns); 1057 for (let p=0; p < numInstances; p++) { 1058 for (let r=0; r < nruns; r++) { 1059 pp[r] = runs[r].coords[p]; 1060 } 1061 stats[p] = CSynth.stats(pp); 1062 } 1063 1064 const radv = stats.map(v=>v.radii[0] + v.radii[1] + v.radii[2] ) 1065 const res = {runs, stats, radv}; 1066 CSynth.compareRuns.res = res; 1067 msgfix('runs', 'complete: stats ready'); 1068 return res; 1069 1070 // CSynth.xyzsExact(runs[1]) 1071 1072 /** 1073 const r = []; 1074 for (let run=0; run<nruns; run++) { 1075 r[run] = []; 1076 for (let run1=0; run1<run; run1++) { 1077 const rmse = CSynth.rmse(runs[run].textureData, runs[run1].textureData); 1078 r[run][run1] = r[run1][run] = rmse 1079 } 1080 } 1081 return {runs, r} 1082 **/ 1083} 1084// let rundata = CSynth.compareRuns(5) 1085 1086var Plane; 1087/** 1088Make planes to illustrate the inside packing of the 60 capsid chain groups. 1089*/ 1090CSynth.planesForInside = function(glmol = CSynth.glmol["All_capsid_proteins.pdb"], maxd=120, chains=undefined, pgroup = CSynth.rawgroup) { 1091 if (typeof glmol === 'string') glmol = CSynth.glmol[glmol]; 1092 const planes = CSynth.planes = []; 1093 const atoms = glmol.atoms; 1094 let ccc = []; 1095 let lastChaingroup, lastcol, lastResi; 1096 glmol.colorByChaingroup(); 1097 const bychain = typeof maxd === 'string'; 1098 for (let i=0; i<atoms.length; i++) { 1099 const atom = atoms[i]; 1100 if (!atom) continue; 1101 if (chains && chains.indexOf(atom.chain) === -1) continue; 1102 if (atom.chaingroup !== lastChaingroup || atom.resi < lastResi) 1103 makePlane(lastcol); 1104 lastChaingroup = atom.chaingroup; 1105 lastResi = atom.resi; 1106 lastcol = atom.color; 1107 if (bychain ? atom.chain === maxd : atom.originDistance < maxd) 1108 ccc.push(atom); 1109 } 1110 makePlane(); 1111 return Plane.drawSet(planes, 'planesInside'); 1112 1113 function makePlane(col) { 1114 if (ccc.length === 0) return; 1115 const stats = CSynth.stats(ccc); 1116 const ax3 = stats.eigenvectors[2].normalize(); 1117 const dot = ax3.dot(stats.centroid); 1118 if (dot < 0) 1119 ax3.multiplyScalar(-1); 1120 planes.push(new Plane(ax3, stats.centroid, col3().set(col))); 1121 ccc = []; 1122 } 1123 1124} 1125 1126CSynth.manyPlanes = function(show = false, fid=undefined, chains=undefined) { 1127 if (fid) CSynth.planesForInside(fid, undefined, chains).visible = show; 1128 if (CSynth.current.fixedPointsData) CSynth.planesForFixed().visible = show; 1129 Plane.icosa().visible = show; 1130 Plane.dodeca().visible = show; 1131 Plane.varying().visible = show; 1132} 1133 1134// sample array to get n samples, a can be an array of anything 1135CSynth.sample = function CSynth_sample(a, n) { 1136 const r = []; 1137 const l = a.length; 1138 for (let i = 0.5; i < n; i++) {
1139 const k = Math.floor(i*l / n); 1140 r.push(a[k]); 1141 } 1142 return r; 1143} 1144 1145// sample array to get n averaged samples, a should be array of numbers 1146CSynth.average = function CSynth_sample(a, n) { 1147 const r = []; 1148 const l = a.length; 1149 let j = 0; 1150 for (let i = 1; i <= n; i++) { // as we are going to compute END of range 1151 let k = Math.floor(i*l / n); 1152 let nn = k - j; 1153 if (nn === 0) { 1154 r.push(a[k]) 1155 } else { 1156 let s = 0; 1157 for( ;j < k; j++) s += a[j] ; 1158 r.push(s/nn); 1159 } 1160 } 1161 return r; 1162} 1163 1164 1165// look at the gaps between fixed points in backbone dustance and geometric distance 1166// if flag, only use 1167CSynth.fixedCheck = function(flag) { 1168 const cc = CSynth.current; 1169 const cos = cc.fixedPointsData.coords; 1170 let xc = cc.extraContacts; 1171 if (flag) xc = xc.filter( e => e.flag === flag); 1172 const r = []; 1173 for (let i = 1; i < xc.length; i++) { 1174 const xa = xc[i-1]; 1175 const xb = xc[i]; 1176 const bpa = xa.bp; 1177 const bpb = xb.bp; 1178 const bpdist = bpb - bpa; 1179 const fposa = xa.fpos; 1180 const fposb = xb.fpos; 1181 const pa = cos[fposa]; 1182 const pb = cos[fposb]; 1183 const gdist = pa.distanceTo(pb); 1184 r.push({bpdist, gdist, rel: gdist/bpdist, bpa, bpb, flags: xa.flag + xb.flag}); 1185 } 1186 log(r); 1187 return r; 1188} 1189 1190/** compute average dist for given backbone dist, idef can be coordinates or code (see CSynth.pos) */ 1191CSynth.avgDist = function(idef = 'cur', bd = 1) { 1192 const d = CSynth.pos(idef); 1193 let sd = 0; 1194 for (let i = bd; i < d.length; i++) { 1195 sd += d[i].distanceTo(d[i-bd]); 1196 } 1197 return sd / (d.length - bd); 1198} 1199 1200/** get volume, very course grid based estimate */ 1201CSynth.volume = function({res = 1, p = springs.getpos()} = {}) { 1202 const st = CSynth.stats(p); 1203 const lx = Math.floor(st.min.x / res), hx = Math.floor(st.max.x / res), nx = hx-lx+1; 1204 const ly = Math.floor(st.min.y / res), hy = Math.floor(st.max.y / res), ny = hy-ly+1; 1205 const lz = Math.floor(st.min.z / res), hz = Math.floor(st.max.z / res), nz = hz-lz+1; 1206 const n = new Uint16Array(nx*ny*nz); 1207 p.forEach(v => { 1208 const ix = Math.floor(v.x/res) - lx, iy = Math.floor(v.y/res) - ly, iz = Math.floor(v.z/res) - lz; 1209 //if (ix < 0 || ix > nx || iy < 0 || iy > ny || iz < 0 || iz > nz) 1210 // log('err'); 1211 n[ix + iy * nx + iz * nx * ny]++; 1212 }); 1213 const v = n.reduce((c,x) => c + (x > 0), 0); 1214 log('vol', v, 'of', n.length); 1215 return v; 1216} 1217 1218/** work in progress, mutation for CSynth parameters */ 1219CSynth.mutate = async function(w = 2000, use = 'pushapartpow springpow'.split(' ')) { 1220 if (vps[0] + vps[1] === 0) setViewports([3,3]); 1221 // V.gui.visible = false; 1222 for (let gn in genedefs) {genedefs[gn].free = +use.includes(gn)} 1223 mutate(); 1224 for (const slot of slots) { 1225 if (slot === undefined) continue; 1226 const dispobj = slot.dispobj; 1227 setObjUniforms(dispobj.genes, uniforms); 1228 CSynth.randpos(undefined, 1); 1229 springs.step(w); 1230 CSynth.autoscale(); 1231 dispobj.render(1); 1232 await S.frame(); 1233 } 1234} 1235 1236// w.i.p. display 8 objects 1237var alwaysNewframe, dustbinvp, refall; 1238CSynth.draw8 = function() { 1239 // the prerenderObj callback will ensure correct spring positions and genes 1240 // before the slot/dispobj is rerendered 1241 if (CSynth.prerenderObj) Maestro.remove('prerenderObj', CSynth.prerenderObj) 1242 Maestro.on('prerenderObj', (evt,a,b,c) => { 1243 const dispobj = evt.dispobj; 1244 const vn = dispobj.vn; 1245 CSynth.gotoCapture(vn); 1246 slots[vn].dispobj.genes = currentGenes; 1247 }); 1248 // alwaysNewframe = 0; // 1249 if (vps[0] + vps[1] === 0) setViewports([4,2]); 1250 refall(); // ensure all objects are repainted 1251 dustbinvp = 99; 1252} 1253 1254var skelbuffer, readWebGlFloatDirect; 1255/** make stats from skeleton */ 1256function skelstats() { 1257 const d = readWebGlFloatDirect(skelbuffer, {width:U.skelnum+1 + 2*U.skelends, height: U.horncount}) 1258 let maxx = -Number.MAX_VALUE, maxy = -Number.MAX_VALUE, maxz = -Number.MAX_VALUE; 1259 let minx = Number.MAX_VALUE, miny = Number.MAX_VALUE, minz = Number.MAX_VALUE; 1260 const mxx = Math.max; 1261 const mnn = Math.min; 1262 const n = d.length; 1263 1264 // compute raw sums, assuming the w is radius used as weight 1265 let sx=0, sy=0, sz=0, sxx=0, syy=0, szz=0, sxy=0, syz=0, szx=0, sw = 0; 1266 for (let i=0; i < n; ) { 1267 const x = d[i++], y = d[i++], z = d[i++], w = d[i++]; 1268 sx += x * w; 1269 sy += y * w; 1270 sz += z * w; 1271 maxx = mxx(maxx, x); maxy = mxx(maxy, y); maxz = mxx(maxz, z); 1272 minx = mnn(maxx, x); miny = mnn(maxy, y); minz = mnn(maxz, z); 1273 sxx += x * x * w; 1274 syy += y * y * w; 1275 szz += z * z * w; 1276 sxy += x * y * w; 1277 syz += y * z * w; 1278 szx += z * x * w; 1279 sw += w; 1280 } 1281 const r = CSynth._stats2({sx, sy, sz, sxx, syy, szz, sxy, syz, szx, n:sw}) 1282 Object.assign(r, {minx, miny, minz, maxx, maxy, maxz}) 1283 r.rx = maxx - minx; r.ry = maxy - miny; r.rz = maxz - minz; 1284 r.rmax = Math.max(r.rx, r.ry, r.rz); 1285 return r; 1286} 1287// CSynth.showEigen('angleed', skelstats()); 1288 1289/** align the current conformation (assumed distances for now) with the given fixed positions 1290using pullsprings */ 1291CSynth.alignConformation = async function({maxpull = 1, minpull = 0.001, degrade = 0.9} = {}) { 1292 console.time('align') 1293 const cc = CSynth.current; 1294 if (!CSynth.alignmentTarget) return msgfixerrorlog('align request with no target set'); 1295 1296 setAllLots('force', 1e-26); // not sure why 0 fails 1297 G.springmaxvel = 9999; // make high so we can get really fast convergence 1298 G.xyzforce = 0.01; // strongest we can get away with for Lorentz initial test 1299 G.pullskelforce = 0; // relevant ??? 1300 springs.finishFix(); // in case in a fixed position already 1301 1302 CSynth.alignmentTarget.forEach((v,i,a) => springs.addpull(i, v.x, v.y, v.z, 1)); 1303 let p = maxpull, i = 0; 1304 while (p > minpull) { 1305 G.pullspringforce = p; 1306 springs.step(1); 1307 if (i++ % 40 === 0) { 1308 await S.frame(); 1309 log(i, 'pull', p) 1310 } 1311 p *= degrade; 1312 } 1313 G.pullspringforce = 0; 1314 log(i, 'alignment complete with pull', p); 1315 console.timeEnd('align') 1316 1317} 1318 1319/** test alignment between two xyz's, by default move current to previously set target */ 1320CSynth.alignConformationtest = async function(opts = {}) { 1321 1322 GX.setValue('simulationsettings/autoalign', false); 1323 1324 const cc = CSynth.current; 1325 const src = opts.src ?? 'current'; 1326 if (typeof src === 'number') { 1327 const sname = cc.xyzs[src].shortname 1328 GX.getgui('modes/' + sname + '\npositions').press(); 1329 await S.frame(); 1330 GX.getgui('modes/' + sname + '\ndists').press(); 1331 await S.frame(); 1332 } else { 1333 log('capture current position as distance texture') 1334 const texture = CSynth.xyzsToTexture(springs.getpos()); 1335 U.distbuff = texture; 1336 } 1337 1338 await CSynth.alignConformation(opts); 1339} 1340 1341/** set the alignment target, but do not use it yet */ 1342CSynth.setAlignmentTarget = function(targ = 'current') { 1343 const cc = CSynth.current; 1344 if (typeof targ === 'number') 1345 CSynth.alignmentTarget = cc.xyzs[targ].coords; 1346 else 1347 CSynth.alignmentTarget = springs.getpos(); 1348} 1349 1350CSynth.alignConformationNow = async function(opts = {}) { 1351 if (!CSynth.alignmentTarget) return msgfixerrorlog('align request with no target set'); 1352 const p2 = springs.getpos(); 1353 CSynth.alignTransform(CSynth.alignmentTarget, p2, {apply: true}); 1354 CSynth.usedata(p2); 1355} 1356 1357/** align by transform two array vectors p1 and p2 */ 1358CSynth.alignTransform = function(p1, p2, { 1359 loop = 1000, apply = false, 1360 sample = 1, tol = 1e-3, vdamprate = 0.8, 1361 t = [0,0,0,0,0,0], mirror = false, 1362 imax = 6, loglow = 2, lograte = 50 1363 } = {}) { 1364 console.time('aligntr') 1365 const cc = CSynth.current; 1366 if (typeof p1 === 'number') p1 = cc.xyzs[p1].coords; 1367 if (typeof p2 === 'number') p2 = cc.xyzs[p2].coords; 1368 let mmirror = new THREE.Matrix4().makeScale(-1, 1, 1); 1369 let m = new THREE.Matrix4(); 1370 const dp = 0.1, dr = 0.1;
1371 let del = [dp, dp, dp, dr, dr, dr]; 1372 const tv = new THREE.Vector3(); 1373 function makem() { 1374 m.makeRotationFromEuler({x:t[3], y:t[4], z:t[5], order: 'XYZ'}).setPosition(t[0], t[1], t[2]); 1375 if (mirror) m.multiply(mmirror); 1376 return m; 1377 } 1378 let l; 1379 1380 const dd = (p) => p1[p].distanceToSquared(tv.copy(p2[p]).applyMatrix4(m)); // error or particle 1381 // const ddd = () => { makem(); return p1.reduce( (c,v,p) => c + dd(p), 0); }; // overall error 1382 const ddd = (ss = sample) => { makem(); let c = 0; for (let p=0; p < p1.length; p += ss) c += dd(p); return c; }; // overall error 1383 let dampv = ddd(); // used for convergence test 1384 let err = () => Math.min(dampv, Math.abs(ddd()-dampv) / dampv); 1385 1386 // perform one test interation on one variable; could do several on given i so only need fewer ddd() tests 1387 function test(i, d = del[i]) { 1388 const t1 = t[i]; 1389 const v1 = ddd(); 1390 const t2 = t[i] = t1 + d; 1391 const v2 = ddd(); 1392 let v; 1393 if (v1 > v2) { 1394 del[i] *= 1.1; 1395 v = v2; 1396 } else { 1397 del[i] *= -0.9; 1398 t[i] = t1; 1399 v = v1; 1400 } 1401 if (l < loglow) log('test', {l, i, d, t1, t2, v1, v2, v, dampv, err: err()}) 1402 if (l%lograte === 0 && i === 0) log({l, del, t, v , dampv, err: err().toString()}); 1403 return v; 1404 } 1405 const test6 = _=> { let v; for (let i = 0; i < imax; i++) v = test(i); return v; }; // perform single iteration on all variables 1406 1407 const drate = vdamprate; 1408 const vv = new Float32Array(loop); 1409 // W.xvv = vv; 1410 for (l = 0; l < loop; l++) { 1411 const v = test6(); // perform all iterations on all variables 1412 vv[l] = v; 1413 dampv = drate * dampv + (1-drate) * v; 1414 if (l > 20 && err() < tol) 1415 break; 1416 } 1417 log('finish at', {l, vsamp: ddd(), vfull: ddd(1), dampv, err: err(), del, t}); 1418 1419 if (apply) { 1420 for (let p = 0; p < p2.length; p++) p2[p].applyMatrix4(m); 1421 const f = cc.xyzs.find(x => x.coords === p2); 1422 if (f) 1423 GX.getgui('modes/' + f.shortname + '\npositions').press(); 1424 } 1425 log('mat', m.toString()); 1426 log('t', t.toString()) 1427 console.timeEnd('aligntr') 1428 return m; 1429} 1430 1431springs.velstats = function() { 1432 const a = springs.getpos() 1433 springs.step(1) 1434 const b = springs.getpos() 1435 const v = a.map((x,i) => x.distanceTo(b[i])); 1436 return getstats(v); 1437} 1438 1439/** compute running stats of half width w, eg point i will use values for i-w to i+w, 2*w+1 items */ 1440CSynth.runningStats = function(w = 10, p = springs.getpos()) { 1441 // const stats1 = {sx:0, sy: 0, sz: 0, sxx: 0, syy:0, szz:0, sxy: 0, syz:0, szx: 0, n:0 } 1442 const ccx = 1, ccy = 1, ccz = 1; // ?? not used in a signficant way 1443 let sx, sy, sz, sxx, syy, szz, sxy, syz, szx, n=0; 1444 1445 const zero = () => { sx = 0, sy= 0, sz= 0, sxx= 0, syy=0, szz=0, sxy= 0, syz=0, szx= 0, n=0; } 1446 zero(); 1447 const r = []; 1448 const add = pi => { 1449 if (!pi) return; 1450 const {x, y, z} = pi; 1451 sx += x; sy += y; sz += z; 1452 sxx += x*x; syy += y*y; szz += z*z; 1453 sxy += x*y; syz += y*z; szx += z*x; 1454 n++; 1455 } 1456 const sub = pi => { 1457 if (!pi) return; 1458 const {x, y, z} = pi; 1459 sx -= x; sy -= y; sz -= z; 1460 sxx -= x*x; syy -= y*y; szz -= z*z; 1461 sxy -= x*y; syz -= y*z; szx -= z*x; 1462 n--; 1463 } 1464 // for (let i = 0; i < w; i++) { 1465 // add(p[i]); 1466 // } 1467 for (let i = 0; i < p.length+w; i++) { 1468 add(p[i]); 1469 if (p[i-w]) { 1470 const 1471 xx = (sxx - sx * sx / n) / n * ccx * ccx, 1472 yy = (syy - sy * sy / n) / n * ccy * ccy, 1473 zz = (szz - sz * sz / n) / n * ccz * ccz, 1474 xy = (sxy - sx * sy / n) / n * ccx * ccy, 1475 yz = (syz - sy * sz / n) / n * ccy * ccz, 1476 zx = (szx - sz * sx / n) / n * ccz * ccx, 1477 cx = sx / n, // * ccx; 1478 cy = sy / n, // * ccy; 1479 cz = sz / n, // * ccz; 1480 tr = (xx + yy + zz) ** 0.5; 1481 1482 r[i-w] = {i: i-w, x:cx, y:cy,z:cz, tr, xx, yy, zz, xy, yz, zx, n} // , sx, sy, sz, sxx, syy, szz, sxy, syz, szx, n, } 1483 } 1484 sub(p[i - 2*w]); 1485 } 1486 1487 const rl = r[0], rh = r[p.length-1]; 1488 r.forEach((x,i,a) => { 1489 x.dw = distxyz(a[i-w] ?? rl, a[i+w] ?? rh); // inter block dist 1490 // x.tr1 = a[i-w]?.tr; // trace to leftinter block dist 1491 }); 1492 1493 return r; 1494} 1495 1496/** generate bed file automatically, w is running width */ 1497CSynth.genbed = function(w = 10, thresh = 15, peakwidth = 5) { 1498 const rs = CSynth.runningStats(w); 1499 const r = rs.map(x => x.dw); 1500 const tr = rs.map(x => x.tr); // /20 1501 //const tr1 = rs.map(x => x.tr1/20); 1502 1503 /** max value from i-n to i+n */ 1504 const mm = (a, i, n) => a.slice(i-n, i+n).reduce((c,v) => Math.max(c,v), -Infinity); 1505 1506 CSynth.clearMarkers(); 1507 r.forEach((x,i,a) => { 1508 if (x > thresh && x === mm(r, i, peakwidth) && mm(tr, i, peakwidth) < x) 1509 CSynth.setMarker(-1, i*5000) 1510 }); 1511 CSynth.markers2Bed(); 1512 1513 CSynth.plot([{data: r, label: 'dw'}, {data:tr, label: 'tr'}]); 1514} 1515 1516 1517var Chart; 1518CSynth.plotev = function plotev(e) { 1519 const chart = CSynth.chart; 1520 // coordinates of click relative to canvas 1521 const { x, y }
1521 = Chart.helpers.getRelativePosition(e, chart); 1522 // can also use const x = e.native.offsetX, y = e.native.offsetY; 1523 1524 // get values relative to chart axes 1525 const dataX = chart.scales.x.getValueForPixel(x); 1526 // const dataY = chart.scales.y.getValueForPixel(y); 1527 CSynth.setMarker(15, CSynth.bp4particle(dataX), 'chart'); 1528 1529 // TODO 8 is to allow for standard margin+padding when dragable 1530 W.plothr.style.left = (x+8) + 'px'; // e.clientX + 'px' 1531 // log({x, dataX, off: offx(e)}) 1532} 1533 1534// xx=(data,label) => ({data:CSynth.makess(data).map(x=>Math.log10(x)), label}); CSynth.plot([xx(ccc1.unitData, 'raw'), xx(ccc1.textureData, 'csynorm'),xx(ccc0.unitData, 'xnorm'), xx(ccc0.textureData, 'xnorm,csynorm')] ) 1535CSynth.plot = function(rlabel, {xmin, xmax, ymin = -0.02, ymax = 0.02, type = 'line', lineWidth = 1, pointRadius = 0} = {}) { 1536 if (!Chart) { 1537 log('adding Chart') 1538 addscript("https://cdn.jsdelivr.net/npm/chart.js"); 1539 setTimeout(() => CSynth.plot(rlabel), 1000); 1540 return; 1541 } 1542 if (Array.isArray(rlabel) && typeof rlabel[0] === 'number') rlabel = [{data: rlabel, label: '?'}] 1543 if (!Array.isArray(rlabel)) rlabel = [rlabel]; 1544 1545 let div = CSynth.plotdiv, canvas = CSynth.plotcanvas; 1546 if (!div) { 1547 div = CSynth.plotdiv = document.createElement('div'); 1548 div.innerHTML = `<hr id="plothr" style=" 1549 width:1px; 1550 height:100%; 1551 left: 80%; 1552 display: block; 1553 position: absolute; 1554 background: white; 1555 opacity: 0.3; 1556 border: none; 1557 ">` 1558 div.style = 'position:fixed; top:0px; right:0px; z-index:9999; width:600px; height: 400px; background: rgba(0,0,0,0.9); border: 1px green solid' 1559 document.body.appendChild(div); 1560 canvas = CSynth.plotcanvas = document.createElement('canvas'); 1561 div.appendChild(canvas); 1562 div.addEventListener('mousemove', CSynth.plotev); 1563 makeDraggable(div, {button: 1, movecallback: canvsize, upcallback: sizechange}); 1564 } 1565 if (CSynth.chart) { CSynth.chart.destroy(); CSynth.chart = undefined; } 1566 function canvsize() { // size so you can see what is happening 1567 canvas.style.width = div.style.width; 1568 canvas.style.height = div.style.height; 1569 } 1570 function sizechange() { // proper redraw for decent quality 1571 const h = +div.style.height.pre('px'), w = +div.style.width.pre('px'); 1572 if (canvas.width === w && canvas.height === h) return; 1573 CSynth.chart.options.aspectRatio = w/h; 1574 CSynth.chart.resize(w, h); 1575 } 1576 1577 // let chart; 1578 const datasets = []; 1579 for (let i=0; i<rlabel.length; i++) 1580 datasets.push(Object.assign( 1581 { borderWidth: lineWidth, // borderWidth is lineWidth 1582 pointRadius}, rlabel[i])) 1583 1584 const cfg = { 1585 type, 1586 data: { 1587 datasets, 1588 labels: new Array(rlabel[0].data.length).fill(''), // labels needed otherwise it collapses x, to give a vertical line 1589 }, 1590 options: { 1591 responsive: false, // we'll handle in makeDraggable 1592 aspectRatio: CSynth.plotdiv.clientWidth / CSynth.plotdiv.clientHeight, 1593 onClick: CSynth.plotev, 1594 onMousemove: CSynth.plotev, 1595 onMouseMove: CSynth.plotev, 1596 animation: false, 1597 //scales: {yAxes: [{ticks: {max: 440.2}, {ticks: {max: 0.0002}]}} 1598 scales: { 1599 y: {max: ymax, min:ymin}, 1600 xAxes: {type: 'linear', position: 'bottom', min: xmin, max: xmax} 1601 } 1602 } 1603 // labels: r, 1604 } 1605 1606 try { 1607 CSynth.chart = new Chart(CSynth.plotcanvas, cfg); 1608 if (!CSynth.chart) {console.error('failed to make chart'); return; } 1609 } catch(e) { 1610 console.error('new CSynth chart call failed', e.message); 1611 document.body.removeChild(CSynth.plotdiv); 1612 CSynth.plotdiv = 0; // force refresh, cChart.js confused ??? and has done some allocation work 1613 return; 1614 } 1615 sizechange(); 1616 // chart corrupts the size, so reset it 1617 1618 //await sleep(500) 1619 //CSynth.plotdiv.style.width = '600px'; // = 'position:fixed; top:0px; right:0px; z-index:9999; width:600px; height: 400px;' 1620 //CSynth.plotdiv.style.height = '400px'; 1621 //log(cfg) 1622} 1623 1624/** compute medial filter for near diagonal elements */ 1625CSynth.medial = function({c = U.contactbuff.source.data.data, h = 250, hstep = 5, w = 15, wstep = 1, perc = 0.6, maxr = Infinity, avoid = 1, plot = true, deadends = 0} = {}) { 1626 console.time('medial'); 1627 const n = Math.round(c.length ** 0.5); 1628 const res = [] 1629 for (const neg of [false, true]) { 1630 const r = new Float32Array(n); // to collect result 1631 const dr = new Float32Array(n); // to collect result 1632 const xy = new Float32Array(h * (2*w + 1)); // to collect contributing elements, reuse each i 1633 for (let i=deadends; i < n-deadends; i++) { // for each particle 1634 if (c[i*n + i] < 0) {r[i] = r[i-1] ?? 0; continue; } // to handle blank regions 1635 let p = 0; 1636 for (let x = i-w; x <= i+w; x += wstep) { // for neighbours 1637 if (x < 0 || x >= n) continue; 1638 let [st, et] = neg ? [x+avoid, x+h] : [x-h, x-avoid] 1639 for (let y = st; y < et; y += hstep) { 1640 if (y < 0 || y >= n || Math.abs(x-y) <= avoid) continue; 1641 const v = c[y*n + x]; 1642 if (v < 0) continue; // eg -999 1643 xy[p++] = v; 1644 } 1645 } 1646 const ss = xy.subarray(0,p).sort((x,y) => x-y); 1647 r[i] = p === 0 ? 0 : Math.min(maxr, xy[Math.floor(p * perc)]); 1648 } 1649 1650 const rl = r[0], rh = r.slice(-1)[0]; 1651 const negi = neg ? -1 : 1; 1652 const negs = neg ? 'start' : 'end' 1653 for (let i = deadends+w; i < n - deadends-w; i++) { 1654 dr[i] = ((r[i+w]??rh) - (r[i-w]??rl)) * negi 1655 } 1656 res.push({data: r, label: 'medial ' + negs}, {data: dr, label: 'dmedial ' + negs}); 1657 } // neg 1658 console.timeEnd('medial'); 1659 if (plot) CSynth.plot(res) 1660 return res; 1661} 1662 1663var ccc0; 1664/** test that missing values always go in complete rows/cols */ 1665CSynth.testMissingRows = function(a = ccc0.rawData, m = -999, n = numInstances) { 1666 const ra=[], rb=[]; 1667 for (let i=0; i<n;i++) { 1668 for (let j=0; j<n;j++) { 1669 const v = a[i+j*n]; 1670 const vi = a[i+i*n]; 1671 const vj = a[j+j*n]; 1672 if (v < 0 && vi >= 0 && vj >= 0) 1673 ra.push([i,j]); 1674 if (v >= 0 && (vi < 0 ||vj < 0)) 1675 rb.push([i,j]); 1676 } 1677 } 1678 return [ra, rb]; 1679}
Line numbers count LF bytes from the start of the resource, as the search results do. Vendor segments are library code the classifier recognised; they are stored but not indexed. Bytes are shown as Latin1 characters, one per byte.