|
| 1 | +// Generate src/data/presets.json: one preset per realizable onset x offset |
| 2 | +// dynamotype class. For each pair we search our actual bifurcation curves for |
| 3 | +// control points that produce a seizure the model classifies as that exact |
| 4 | +// pair, choosing the simplest path type (arc -> circle -> piecewise) that works. |
| 5 | +// |
| 6 | +// Run from this directory: node make_presets.mjs |
| 7 | +import { simulate } from "../src/model.js"; |
| 8 | +import MAP from "../src/data/curves.json" with { type: "json" }; |
| 9 | +import { writeFileSync } from "node:fs"; |
| 10 | + |
| 11 | +const R = 0.4, TSTEP = 0.02; |
| 12 | +const ICTAL = [-0.043, -0.209, -0.338]; // a vigorous seizure-core direction (model frame) |
| 13 | + |
| 14 | +/* ---- geometry (mirrors App.jsx) ---- */ |
| 15 | +const onR = v => { const m = Math.hypot(...v) || 1e-9; return [v[0]/m*R, v[1]/m*R, v[2]/m*R]; }; |
| 16 | +const dot = (a,b) => a[0]*b[0]+a[1]*b[1]+a[2]*b[2]; |
| 17 | +const cross = (a,b) => [a[1]*b[2]-a[2]*b[1], a[2]*b[0]-a[0]*b[2], a[0]*b[1]-a[1]*b[0]]; |
| 18 | +const segN = (ang,k) => Math.max(2, Math.round(ang/(k*TSTEP))); |
| 19 | +function geoArc(A,B,k){ A=onR(A);B=onR(B); const th=Math.acos(Math.max(-1,Math.min(1,dot(A,B)/(R*R)))); if(th<1e-4)return[A.slice()]; const s=Math.sin(th),n=segN(th,k),o=[]; for(let i=0;i<=n;i++){const t=i/n,a=Math.sin((1-t)*th)/s,b=Math.sin(t*th)/s;o.push([a*A[0]+b*B[0],a*A[1]+b*B[1],a*A[2]+b*B[2]]);} return o; } |
| 20 | +function circleArc(A,V,B,k){ A=onR(A);V=onR(V);B=onR(B); let nr=cross([V[0]-A[0],V[1]-A[1],V[2]-A[2]],[B[0]-A[0],B[1]-A[1],B[2]-A[2]]); const nl=Math.hypot(...nr); if(nl<1e-7)return geoArc(A,B,k); nr=[nr[0]/nl,nr[1]/nl,nr[2]/nl]; const d=dot(nr,A),C=[nr[0]*d,nr[1]*d,nr[2]*d],r=Math.sqrt(Math.max(1e-9,R*R-d*d)); let E=[A[0]-C[0],A[1]-C[1],A[2]-C[2]]; const el=Math.hypot(...E); E=[E[0]/el,E[1]/el,E[2]/el]; const F=cross(nr,E); const ang=p=>{const q=[p[0]-C[0],p[1]-C[1],p[2]-C[2]];return Math.atan2(dot(q,F),dot(q,E));}; const TAU=2*Math.PI,wrap=x=>{x%=TAU;return x<0?x+TAU:x;}; const aV=wrap(ang(V)),aB=wrap(ang(B)); let dir,total; if(aV<=aB){dir=1;total=aB;}else{dir=-1;total=TAU-aB;} const n=segN(total,k),o=[]; for(let i=0;i<=n;i++){const th=dir*total*i/n,c=Math.cos(th),s=Math.sin(th);o.push([C[0]+r*(c*E[0]+s*F[0]),C[1]+r*(c*E[1]+s*F[1]),C[2]+r*(c*E[2]+s*F[2])]);} return o; } |
| 21 | +function piece(P,k){ let o=[]; for(let i=0;i<P.length-1;i++){const s=geoArc(P[i],P[i+1],k);o=o.concat(i?s.slice(1):s);} return o; } |
| 22 | + |
| 23 | +/* ---- classify (mirrors App.jsx) ---- */ |
| 24 | +const minD = (p,cv) => { let m=1e9; for(const q of cv.pts){const d=Math.hypot(p[0]-q[0],p[1]-q[1],p[2]-q[2]); if(d<m)m=d;} return m; }; |
| 25 | +function classify(point, role){ |
| 26 | + let best=1e9, hit=null; |
| 27 | + for(const cv of MAP.curves){ if(cv.role!=="both"&&cv.role!==role)continue; const dm=minD(point,cv); if(dm<best){best=dm;hit=cv;} } |
| 28 | + if(hit && role==="onset" && (hit.dyno==="SN"||hit.dyno==="SupH")){ |
| 29 | + for(const d of ["SNIC","SubH"]){ const cv=MAP.curves.find(c=>c.dyno===d); if(cv&&minD(point,cv)<0.03)return cv; } |
| 30 | + } |
| 31 | + return hit; |
| 32 | +} |
| 33 | +const ONSET = new Set(["SN","SNIC","SupH","SubH"]); |
| 34 | +const OFFSET = new Set(["SNIC","SH","SupH","FLC"]); |
| 35 | +const ptsOf = dyno => [].concat(...MAP.curves.filter(c=>c.dyno===dyno).map(c=>c.pts)); |
| 36 | + |
| 37 | +const niceMargins = m => m && m.onset>0.08 && m.onset<0.45 && m.offset>0.55 && m.offset<0.94; |
| 38 | + |
| 39 | +// score a candidate: returns {ok, score} — prefer clean margins, seizure present, correct class |
| 40 | +function evalPath(path, oD, fD){ |
| 41 | + const o = simulate(path, { fs:200, dur:16, freq:10, noise:0, drift:0, seed:7 }); |
| 42 | + if(!o.markers) return null; |
| 43 | + const onPt = path[0], offPt = path[path.length-1]; |
| 44 | + if(classify(onPt,"onset").dyno !== oD) return null; |
| 45 | + if(classify(offPt,"offset").dyno !== fD) return null; |
| 46 | + // score: closeness of margins to ideal (on 0.2 / off 0.8) |
| 47 | + const s = Math.abs(o.markers.onset-0.2) + Math.abs(o.markers.offset-0.8) + (niceMargins(o.markers)?0:0.5); |
| 48 | + return { score:s, markers:o.markers }; |
| 49 | +} |
| 50 | + |
| 51 | +const FRACS = [0.15,0.3,0.45,0.6,0.75,0.9]; |
| 52 | +const VIA_CANDS = [ICTAL, [0.0,-0.3,-0.1], [-0.1,-0.25,0.1], [0.1,-0.28,-0.05], [-0.2,-0.2,-0.2]]; |
| 53 | +const K = 0.02; |
| 54 | + |
| 55 | +// path type the DfD tutorial associates with each family |
| 56 | +function preferType(oD, fD){ |
| 57 | + if(oD === "SN" || oD === "SubH") return "arc"; // hysteresis-loop family |
| 58 | + if(oD === "SupH" && fD === "SupH") return "piecewise"; // piecewise family |
| 59 | + return "circle"; // slow-wave family |
| 60 | +} |
| 61 | + |
| 62 | +function searchType(type, A, B, oD, fD){ |
| 63 | + const pick = (arr,f)=>arr[Math.floor(f*(arr.length-1))]; |
| 64 | + let best=null; |
| 65 | + if(type==="arc"){ |
| 66 | + for(const fa of FRACS) for(const fb of FRACS){ const on=pick(A,fa),off=pick(B,fb); |
| 67 | + const r=evalPath(geoArc(on,off,K),oD,fD); if(r&&(!best||r.score<best.score)) best={type,pts:[on,off],...r}; } |
| 68 | + } else if(type==="circle"){ |
| 69 | + for(const fa of [0.2,0.4,0.6,0.8]) for(const fb of [0.2,0.4,0.6,0.8]) for(const via of VIA_CANDS){ const on=pick(A,fa),off=pick(B,fb); |
| 70 | + const r=evalPath(circleArc(on,via,off,K),oD,fD); if(r&&(!best||r.score<best.score)) best={type,pts:[on,via,off],...r}; } |
| 71 | + } else { |
| 72 | + for(const fa of [0.3,0.6]) for(const fb of [0.4,0.7]) for(const via of VIA_CANDS){ const on=pick(A,fa),off=pick(B,fb); |
| 73 | + const v2=[(via[0]+off[0])/2,(via[1]+off[1])/2,(via[2]+off[2])/2]; |
| 74 | + const r=evalPath(piece([on,via,v2,off],K),oD,fD); if(r&&(!best||r.score<best.score)) best={type,pts:[on,via,v2,off],...r}; } |
| 75 | + } |
| 76 | + return best; |
| 77 | +} |
| 78 | + |
| 79 | +function findPreset(oD, fD){ |
| 80 | + const A = ptsOf(oD).filter(p=>classify(p,"onset").dyno===oD); |
| 81 | + const B = ptsOf(fD).filter(p=>classify(p,"offset").dyno===fD); |
| 82 | + if(!A.length||!B.length) return null; |
| 83 | + const pref = preferType(oD, fD); |
| 84 | + const order = [pref, ...["arc","circle","piecewise"].filter(t=>t!==pref)]; |
| 85 | + const found = {}; |
| 86 | + for(const t of order){ const r = searchType(t, A, B, oD, fD); if(r) found[t]=r; } |
| 87 | + // use the preferred type if it verifies reasonably cleanly; else the best available |
| 88 | + if(found[pref] && found[pref].score < 0.8) return found[pref]; |
| 89 | + let best=null; for(const t in found) if(!best||found[t].score<best.score) best=found[t]; |
| 90 | + return best; |
| 91 | +} |
| 92 | + |
| 93 | +const round4 = p => p.map(v=>Math.round(v*1e4)/1e4); |
| 94 | +const presets = []; |
| 95 | +for(const oD of ["SN","SNIC","SupH","SubH"]) for(const fD of ["SH","FLC","SNIC","SupH"]){ |
| 96 | + const r = findPreset(oD, fD); |
| 97 | + const name = `${oD} / ${fD}`; |
| 98 | + if(r){ |
| 99 | + presets.push({ name, onset:oD, offset:fD, type:r.type, pts:r.pts.map(round4) }); |
| 100 | + console.log(`${name.padEnd(14)} -> ${r.type.padEnd(9)} on${(r.markers.onset*100)|0}%/off${(r.markers.offset*100)|0}%`); |
| 101 | + } else { |
| 102 | + console.log(`${name.padEnd(14)} -> (no realizable path found)`); |
| 103 | + } |
| 104 | +} |
| 105 | +writeFileSync(new URL("../src/data/presets.json", import.meta.url), JSON.stringify(presets)); |
| 106 | +console.log(`\nwrote presets.json: ${presets.length} classes`); |
0 commit comments