30:T454,flux_plot = { const marks = [ Plot.frame(), Plot.dot(spectrum.data, {x:"wlen", y:"flux", r:1.5, fill:"limegreen"}), Plot.line(coarseFit.fit, {x:"wlen", y:"linear", stroke:"red", strokeDasharray: [2,3]}), Plot.line(finalFit.iknots, {x:i=>Math.exp(kwave[i]), y:(_,i)=>finalFit.F[i], stroke:"black"}) ]; if(wlen_zoom >= 0) marks.push( Plot.dot(smoothedRegions[wlen_zoom], { x:"wlen", y:"flux", r:3, fill:"magenta" }), Plot.dot(prunedRegions[wlen_zoom], { x:"wlen", y:"flux", r:5, stroke:"magenta" }), Plot.line(prunedRegions[wlen_zoom], { x:"wlen", y:"flux", r:5, stroke:"magenta", strokeDasharray:[3,6] }) ); if(regions.merged !== null) marks.push( Plot.rectX(regions.merged, { x1:d=>Math.exp(kwave[d.ilo]), x2:d=>Math.exp(kwave[d.ihi]), fill:"red", fillOpacity:0.2, title:(r,i)=>`region ${i+1}/${regions.merged.length}: ${r.ihi-r.ilo+1} => ${prunedRegions[i].length}` }) ); return Plot.plot({ width: width, x: { domain: wlen_range, grid: true }, y: { domain: fzoom.map(spectrum.fscaled), grid: true }, clip: true, marks }); }31:T4c0,finalKnots = { // Combine the pruned knots from each region. function insert_knots(i1, i2) { const ninsert = Math.floor((i2 - i1 + 1) / oversampling); if(ninsert == 0) return [ ]; const spacing = Math.round((i2 - i1 + 1) / (ninsert + 1)); const iknots = d3.range(ninsert).map(j => i1 + (j + 1) * spacing); console.log(`Inserting ${iknots.length} knots ${iknots} to fit continuum over [${i1},${i2}]`); return iknots; } const iknots = [ 0 ]; let iprev = 0; prunedRegions.forEach((pruned,iregion) => { const {ilo,ihi} = regions.merged[iregion]; // Insert knots for continuum fitting if necessary. const continuum_knots = insert_knots(iprev, ilo); iknots.push(...continuum_knots); // Add knots for this region. const region_knots = prunedRegions[iregion].map(({iknot}) => iknot); console.log(`Adding ${region_knots.length} knots ${region_knots} for region ${iregion} [${ilo},${ihi}]`); iknots.push(...region_knots); iprev = ihi; }); // Insert knots after last region if necessary. const ilast = kwave.length - 1; const continuum_knots = insert_knots(iprev, ilast); iknots.push(...continuum_knots); iknots.push(ilast); return iknots; }32:T675,function getPrunedKnotsContinuous(iregion, {mu=2}={}) { const offset = index * arraySize; const {ilo,ihi} = regions.merged[iregion]; const Fknots = smoothedRegions[iregion].map(({flux}) => flux); const n = ihi - ilo + 1; const ndata = breaks[ihi] - breaks[ilo]; // Precompute the error matrix E[i1,i2] for i1 < i2 const E = new Array(n * (n - 1) / 2).fill(0); let idx = 0; for(let i1=0; i1 < n; ++i1) { const logw1 = kwave[ilo + i1]; const k1 = breaks[ilo + i1]; for(let i2 = i1+1; i2 < n; ++i2) { const logw2 = kwave[ilo + i2], dlogw = logw2 - logw1; const k2 = breaks[ilo + i2]; let chisq = 0; for(let k=k1; k 0) cost += min_cost[i1] + mu * ndata / n; if(cost < min_cost[i2]) { min_cost[i2] = cost; backptr[i2] = i1; } } } // Reconstruct the best subset of knots. const pruned = [ ]; let j = n - 1; while(j != -1) { pruned.push(j); j = backptr[j]; } pruned.reverse(); return pruned.map(i => ({ iknot:i+ilo, wlen:Math.exp(kwave[i+ilo]), flux:Fknots[i] })); }33:T751,function getRegions(fit, {inset=4, pad=3, chisq_cut=4, SG_window_size=19, SG_order=1}={}) { // Calculate smoothed chisq. const chisq = Float64Array.from({length:arraySize}, (_,i) => fit.fit[i].chisq); const chisq_smooth = savitzkyGolay(chisq, SG_window_size, SG_order); // Normalize the cut to the median smooth chisq const chisq_median = d3.median(chisq_smooth); chisq_cut *= chisq_median; // Identify knots bounding regions of smoothed chisq above the cut. const regions = [ ]; let currentRegion = null; for(let i = inset; i < nSegments - inset; ++i) { const k1 = breaks[i], k2 = breaks[i+1]; let above = false; for(let k=k1; k < k2; ++k) { if(chisq_smooth[k] < chisq_cut) continue; above = true; break; } if(above) { // chisq_smooth > chisq_cut somewhere between knots [i,i+1] if(currentRegion === null) { // start a new region currentRegion = { ilo:i-pad, ihi:i+pad }; } else { // extend the current region currentRegion.ihi = i+pad; } } else { // chisq_smooth <= chisq_cut everywhere between knots [i,i+1] if(currentRegion !== null) { // end the current region currentRegion.ihi = i+pad; regions.push(currentRegion); currentRegion = null; } } } // Merge any regions that overlap due to padding. const merged = regions.slice(0,1); // [ ] or [regions[0]] for(let i = 1; i < regions.length; ++i) { const latest = merged[merged.length-1]; if(regions[i].ilo - latest.ihi <= 1) { // merge region i into latest merged region latest.ihi = regions[i].ihi; console.log(`merged region ${i}/${regions.length}`); } else { // start a new isolated region merged.push(regions[i]); } } return { chisq_smooth, chisq_median, merged }; }34:T44e,function savitzkyGolay(y, windowSize, polyOrder) { if (windowSize % 2 === 0 || windowSize < 3) { throw new Error("Window size must be an odd number greater than or equal to 3."); } if (polyOrder >= windowSize) { throw new Error("Polynomial order must be less than the window size."); } const halfWindow = Math.floor(windowSize / 2); const len = y.length; let result = new Array(len).fill(0); // Generate Savitzky-Golay coefficients const coeffs = generateSGCoefficients(windowSize, polyOrder); // Apply convolution with SG coefficients for (let i = 0; i < len; i++) { let smoothedValue = 0; let normFactor = 0; for (let j = -halfWindow; j <= halfWindow; j++) { let idx = i + j; if (idx < 0) idx = 0; // Reflect at the edges if (idx >= len) idx = len - 1; // Reflect at the edges smoothedValue += y[idx] * coeffs[j + halfWindow]; normFactor += coeffs[j + halfWindow]; } result[i] = smoothedValue / normFactor; } return result; }35:T85c,function smoothedKnotFlux(iknot, windowSize, polyOrder) { // Calculate a locally smoothed value of the flux at kwave[iknot] using // linear smoothing in wavelength and accounting for ivar weights. const halfWindow = Math.floor(windowSize / 2); const k0 = breaks[iknot]; const k1 = Math.max(0, k0 - halfWindow); const k2 = Math.min(arraySize, k0 + halfWindow + 1); const offset = index * arraySize; const x0 = Math.exp(kwave[iknot]); const xWindow = Array.from(fullwave.slice(k1, k2)).map(x => x - x0); const yWindow = Array.from(flux.slice(offset + k1, offset + k2)); const wWindow = Array.from(ivar.slice(offset + k1, offset + k2)); if(wWindow.every(w => w === 0)) return NaN; // Build the weighted design (Vandermonde) matrix. // Each row is [sqrt(w) * 1, sqrt(w) * x, sqrt(w) * x^2, ..., sqrt(w) * x^(polyOrder)] let Xprime = []; for (let j = 0; j < xWindow.length; j++) { const sqrtw = Math.sqrt(wWindow[j]); const row = []; for (let p = 0; p <= polyOrder; p++) { row.push(sqrtw * Math.pow(xWindow[j], p)); } Xprime.push(row); } // Convert the design matrix to a math.js matrix. Xprime = math.matrix(Xprime); // Create the weighted y vector: each element is sqrt(weight)*y. let yprime = yWindow.map((val, idx) => Math.sqrt(wWindow[idx]) * val); yprime = math.matrix(yprime); // Compute A = Xprime^T * Xprime. const Xt = math.transpose(Xprime); const A = math.multiply(Xt, Xprime); // Compute b = Xprime^T * yprime. const b = math.multiply(Xt, yprime); // Solve for the polynomial coefficients: A * coeffs = b. // We use math.lusolve for robust linear solving. let coeffs; try { // math.lusolve returns a column vector (matrix) of coefficients. coeffs = math.lusolve(A, b); // Convert the column vector to a flat array. coeffs = coeffs.valueOf().map(item => item[0]); } catch (e) { console.error(`Matrix solve failed: ${e}`); return NaN; } // For the "increasing" Vandermonde matrix, evaluating the fitted polynomial at x=0 // yields the constant term, which is the first coefficient. return coeffs[0]; }36:T93c,function weightedSavitzkyGolay(x, y, weights, windowSize, polyOrder) { const n = y.length; const halfWindow = Math.floor(windowSize / 2); const ySmooth = new Array(n); for (let i = 0; i < n; i++) { // Define window boundaries ensuring we don't exceed array limits. const start = Math.max(0, i - halfWindow); const end = Math.min(n, i + halfWindow + 1); // Extract local window data. // Center the x values by subtracting x[i] so that the current point is at 0. const xWindow = x.slice(start, end).map(val => val - x[i]); const yWindow = y.slice(start, end); const wWindow = weights.slice(start, end); // If all weights in the window are zero, mark as missing. if (wWindow.every(w => w === 0)) { ySmooth[i] = NaN; continue; } // Build the weighted design (Vandermonde) matrix. // Each row is [sqrt(w) * 1, sqrt(w) * x, sqrt(w) * x^2, ..., sqrt(w) * x^(polyOrder)] let Xprime = []; for (let j = 0; j < xWindow.length; j++) { const sqrtw = Math.sqrt(wWindow[j]); const row = []; for (let p = 0; p <= polyOrder; p++) { row.push(sqrtw * Math.pow(xWindow[j], p)); } Xprime.push(row); } // Convert the design matrix to a math.js matrix. Xprime = math.matrix(Xprime); // Create the weighted y vector: each element is sqrt(weight)*y. let yprime = yWindow.map((val, idx) => Math.sqrt(wWindow[idx]) * val); yprime = math.matrix(yprime); // Compute A = Xprime^T * Xprime. const Xt = math.transpose(Xprime); const A = math.multiply(Xt, Xprime); // Compute b = Xprime^T * yprime. const b = math.multiply(Xt, yprime); // Solve for the polynomial coefficients: A * coeffs = b. // We use math.lusolve for robust linear solving. let coeffs; try { // math.lusolve returns a column vector (matrix) of coefficients. coeffs = math.lusolve(A, b); // Convert the column vector to a flat array. coeffs = coeffs.valueOf().map(item => item[0]); } catch (e) { console.error(`Matrix solve failed at index ${i}: ${e}`); ySmooth[i] = NaN; continue; } // For the "increasing" Vandermonde matrix, evaluating the fitted polynomial at x=0 // yields the constant term, which is the first coefficient. ySmooth[i] = coeffs[0]; } return ySmooth; }37:T928,function fitFixedKnotsContinuous(iknots, {ftype=Float64Array}={}) { // Solve the linear least-squares problem to find the optimum continuous piecewise linear // function passing through the specified knots. const {min,max,abs} = Math; const offset = index * arraySize; // Build the coefficient matrix A[i,j] = sum_k { w_k phi_i(x_k) phi_j(x_k) } // and constant vector b[i] = sum_k { w_k f_k phi_i(x_k) } // where phi_i(x) is piecewise linear through (knot[i-1],0) - (knot[i],1) - (knot[i+1],0). const n = iknots.length; const Adiag = new ftype(n).fill(0); const Aupper = new ftype(n).fill(0); const Alower = new ftype(n).fill(0); const b = new ftype(n).fill(0); for(let i = 0; i < n; ++i) { const logwi = kwave[iknots[i]]; // Accumulate A and b. if(i > 0) { // contribution from [j=i-1,i] const logwj = kwave[iknots[i-1]]; for(let k = breaks[iknots[i-1]]; k < breaks[iknots[i]]; ++k) { const t = (logwave[k] - logwj) / (logwi - logwj); const wgt = ivar[offset + k]; Adiag[i] += wgt * t ** 2; Alower[i] += wgt * t * (1 - t); b[i] += wgt * flux[offset + k] * t; } } if(i < n - 1) { // contribution from [i,j=i+1] const logwj = kwave[iknots[i+1]]; for(let k = breaks[iknots[i]]; k < breaks[iknots[i+1]]; ++k) { const t = (logwave[k] - logwi) / (logwj - logwi); const wgt = ivar[offset + k]; Adiag[i] += wgt * (1 - t) ** 2; Aupper[i] += wgt * t * (1 - t); b[i] += wgt * flux[offset + k] * (1 - t); } } } // Solve the tridiagonal system Ax=b const F = thomasAlgorithm(Alower, Adiag, Aupper, b); // Tabulate the fit for each wavelength. let logw_lo, logw_hi = kwave[iknots[0]], chisq_sum = 0; const fit = [ ]; for(let i = 0; i < n - 1; ++i) { logw_lo = logw_hi; logw_hi = kwave[iknots[i+1]]; for(let k = breaks[iknots[i]]; k < breaks[iknots[i+1]]; ++k) { const t = (logwave[k] - logw_lo) / (logw_hi - logw_lo); const linear = F[i] * (1 - t) + F[i+1] * t; const chisq = ivar[offset + k] * (flux[offset + k] - linear) ** 2; //if(!Number.isFinite(chisq)) console.log([k, fullwave[k], j, F[j], F[j+1]]); chisq_sum += chisq; fit.push({ wlen:fullwave[k], linear, chisq }); } } return { F, fit, iknots }; }38:T540,function thomasAlgorithm(a, b, c, d) { /** * Solves a tridiagonal system Ax = d using the Thomas algorithm. * * The tridiagonal matrix A is represented by three arrays: * - a: sub-diagonal (a[0] is unused, so it can be set to 0) * - b: main diagonal * - c: super-diagonal (c[n-1] is unused, so it can be set to 0) * * @param {number[]} a - sub-diagonal coefficients (length n) * @param {number[]} b - main diagonal coefficients (length n) * @param {number[]} c - super-diagonal coefficients (length n) * @param {number[]} d - right-hand side vector (length n) * @returns {number[]} x - solution vector (length n) */ const n = b.length; // Make copies of c and d to store modified coefficients const cp = new c.constructor(n); const dp = new d.constructor(n); // First row modifications cp[0] = c[0] / b[0]; dp[0] = d[0] / b[0]; // Forward sweep: modify the coefficients for (let i = 1; i < n; i++) { // Denominator for the current row const m = b[i] - a[i] * cp[i - 1]; // For the last row, cp[i] is not used so we set it to 0 cp[i] = i < n - 1 ? c[i] / m : 0; dp[i] = (d[i] - a[i] * dp[i - 1]) / m; } // Back substitution: solve for x const x = new d.constructor(n); x[n - 1] = dp[n - 1]; for (let i = n - 2; i >= 0; i--) { x[i] = dp[i] - cp[i] * x[i + 1]; } return x; }39:T590,function fitPrunedKnotsDiscontinuous(iknots, { mu=2 } = {}) { const n = iknots.length - 1; const backptr = new Array(n + 1).fill(-1); const min_cost = new Array(n + 1).fill(Infinity); min_cost[0] = 0; for(let i2=1; i2 < n+1; ++i2) { for(let i1=0; i1 < i2; ++i1) { const sfit = segmentFit(iknots[i1], iknots[i2]); const ndata = breaks[iknots[i2]] - breaks[iknots[i1]]; let cost = sfit.chisq / ndata; if(i1 > 0) cost += min_cost[i1] + mu / iknots.length; if(cost < min_cost[i2]) { min_cost[i2] = cost; backptr[i2] = i1; } } } const pruned = [ ]; let j = n; while(j != -1) { pruned.push(iknots[j]); j = backptr[j]; } pruned.reverse(); // Tabulate the fit for each wavelength. const offset = index * arraySize; let logw_lo, logw_hi = kwave[pruned[0]], chisq_sum = 0; const fit = [ ]; for(let i = 0; i < pruned.length - 1; ++i) { const sfit = segmentFit(pruned[i], pruned[i+1]); logw_lo = logw_hi; logw_hi = kwave[pruned[i+1]]; for(let k = breaks[pruned[i]]; k < breaks[pruned[i+1]]; ++k) { const linear = sfit.linear(logwave[k]); const chisq = ivar[offset + k] * (flux[offset + k] - linear) ** 2; if(!Number.isFinite(chisq)) console.log(k, fullwave[k], i, sfit); chisq_sum += chisq; fit.push({ wlen:fullwave[k], linear, chisq }); } } return { pruned, fit, chisq_sum }; }3a:[{"id":"d5b93090cd0dbf559f028b55a8fe545c02424e75f296b2077e2af49301204d6d942362594bed5697d48ea6bc182e897104c6b76da9ba16b55ba0776c22676be9","href":"https://static.observableusercontent.com/files/d5b93090cd0dbf559f028b55a8fe545c02424e75f296b2077e2af49301204d6d942362594bed5697d48ea6bc182e897104c6b76da9ba16b55ba0776c22676be9","name":"data.bin","type":"application/octet-stream","lastModified":1740750063490.22,"size":13570064},{"id":"ba0ff9eb8a17b01ea4c34f18023fa5440ae8c274880fad2a9b0bf407557c5f0f4485ac527d7c08da2982d2c96a3b96099f9229c8dca60e0f66534485b9f80279","href":"https://static.observableusercontent.com/files/ba0ff9eb8a17b01ea4c34f18023fa5440ae8c274880fad2a9b0bf407557c5f0f4485ac527d7c08da2982d2c96a3b96099f9229c8dca60e0f66534485b9f80279","name":"ivar.bin","type":"application/octet-stream","lastModified":1740750276589.341,"size":13570064}] 2e:{"title":"Piecewise Linear Fits to DESI Spectra","cells":[{"id":0,"value":"# Piecewise Linear Fits to DESI Spectra","pinned":false,"hidden":false,"mode":"md"},{"id":30,"value":"viewof index = Inputs.range([0, nSpectra-1], {label: \"Spectrum Index\", step: 1, value: 5})","pinned":false,"hidden":false,"mode":"ojs"},{"id":816,"value":"viewof wlen_zoom = {\n const n = regions.merged.length;\n return Inputs.button([\n [\"Full Range\", value => -1],\n [\"Next Region\", value => (value + 1) % n],\n [\"Previous Region\", value => value == -1 ? n - 1 : (value + n - 1) % n]\n ], {value: -1, label: \"Zoom\"});\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":167,"value":"viewof fzoom = interval(\n [0, 1],\n { step: 0.01, value: [spectrum.tlo, spectrum.thi], label:\"Flux zoom\" })","pinned":false,"hidden":false,"mode":"ojs"},{"id":880,"value":"viewof fluxWindowSize = Inputs.range([5, 99], {label: \"flux window size\", step: 2, value: 9 })","pinned":false,"hidden":false,"mode":"ojs"},{"id":885,"value":"viewof fluxPolyOrder = Inputs.radio([1,2,3], {label: \"flux poly order\", value: 3})","pinned":false,"hidden":false,"mode":"ojs"},{"id":969,"value":"viewof mu = Inputs.range([0.1, 10], {label: \"Mu Penalty\", step: 0.1, value:2})","pinned":false,"hidden":false,"mode":"ojs"},{"id":66,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":159,"value":"chisq_plot = {\n const [clo, chi] = d3.extent(coarseFit.fit, d => d.chisq);\n const cscaled = t => clo * (1 - t) + chi * t;\n return Plot.plot({\n width: width,\n x: { domain: wlen_range, grid: true },\n y: { domain: czoom.map(cscaled), grid: true, domain:[0,20] },\n clip: true,\n marks: [\n Plot.frame(),\n Plot.dot(coarseFit.fit, {x:\"wlen\", y:\"chisq\", r:1.5, fill:\"limegreen\"}),\n Plot.line(regions.chisq_smooth, {x:(_,i)=>fullwave[i], y:y=>y, stroke:\"red\"}),\n Plot.rectX(regions.merged, {\n x1:d=>Math.exp(kwave[d.ilo]), x2:d=>Math.exp(kwave[d.ihi]),\n fill:\"red\", fillOpacity:0.2,\n title:r=>`${r.ihi-r.ilo+1}` }),\n Plot.ruleY([regions.chisq_median, chisq_cut * regions.chisq_median], {\n y:d=>d, stroke: \"blue\", strokeDasharray: [3,3]\n }),\n ]\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1097,"value":"viewof pruneContinuous = Inputs.toggle({label: \"Continous Prune\", value: true})","pinned":false,"hidden":false,"mode":"ojs"},{"id":184,"value":"viewof czoom = interval(\n [0, 1],\n { step: 0.01, value: [0,1], label:\"ChiSq zoom\" })","pinned":false,"hidden":false,"mode":"ojs"},{"id":478,"value":"viewof ChisqWindowSize = Inputs.range([5, 99], {label: \"chisq window size\", step: 2, value: 19 })","pinned":false,"hidden":false,"mode":"ojs"},{"id":474,"value":"viewof ChisqPolyOrder = Inputs.radio([1,2,3], {label: \"chisq poly order\", value: 1})","pinned":false,"hidden":false,"mode":"ojs"},{"id":574,"value":"viewof chisq_cut = Inputs.range([2, 10], {label: \"Smoothed ChiSq Cut\", step: 0.2, value: 4})","pinned":false,"hidden":false,"mode":"ojs"},{"id":596,"value":"viewof pad = Inputs.radio([0,1,2,3], {label: \"Region Padding\", value: 3})","pinned":false,"hidden":false,"mode":"ojs"},{"id":978,"value":"### Global Fit with Pruned Knots","pinned":false,"hidden":false,"mode":"md"},{"id":986,"value":"finalFit = fitFixedKnotsContinuous(finalKnots)","pinned":false,"hidden":false,"mode":"ojs"},{"id":980,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":924,"value":"### Fit Non-Continuum Regions","pinned":false,"hidden":false,"mode":"md"},{"id":929,"value":"prunedRegions = {\n const method = pruneContinuous ? getPrunedKnotsContinuous : getPrunedKnotsDiscontinuous;\n return d3.range(smoothedRegions.length).map(iregion => method(iregion, {mu}));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1088,"value":"function getPrunedKnotsDiscontinuous(iregion, {mu=2}={}) {\n\n const {ilo,ihi} = regions.merged[iregion];\n const { pruned, fit } = fitPrunedKnotsDiscontinuous(d3.range(ilo, ihi+1), {mu});\n return pruned.map(i => ({ iknot:i, wlen:Math.exp(kwave[i]) }));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":926,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":849,"value":"### Smooth Non-Continuum Regions","pinned":false,"hidden":false,"mode":"md"},{"id":851,"value":"smoothedRegions = {\n const data = [ ];\n const offset = index * arraySize;\n return regions.merged.map(({ilo,ihi}) => {\n return d3.range(ilo, ihi + 1)\n .map(iknot => {\n const wlen = Math.exp(kwave[iknot]);\n const flux = smoothedKnotFlux(iknot, fluxWindowSize, fluxPolyOrder);\n return { iknot, wlen, flux };\n })\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":775,"value":"### Identify Non-Continuum Regions","pinned":false,"hidden":false,"mode":"md"},{"id":819,"value":"wlen_range = {\n if(wlen_zoom == -1 || regions.merged.length == 0) return [ fullwave[0], fullwave[fullwave.length - 1] ];\n const { ilo, ihi } = regions.merged[wlen_zoom];\n const pad = 8;\n return [ Math.max(0, ilo-pad), Math.min(nSegments, ihi+pad) ].map(i => Math.exp(kwave[i]));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":793,"value":"regions = getRegions(coarseFit, {chisq_cut, pad, SG_window_size:ChisqWindowSize, SG_order:ChisqPolyOrder})","pinned":false,"hidden":false,"mode":"ojs"},{"id":783,"value":"coarseFit = {\n const coarseknots = d3.range(nSegmentsCoarse).map(i => i * oversampling);\n coarseknots.push(nSegments);\n return fitFixedKnotsContinuous(coarseknots);\n} ","pinned":false,"hidden":false,"mode":"ojs"},{"id":778,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":446,"value":"### Savitsky-Golay Smoothing","pinned":false,"hidden":false,"mode":"md"},{"id":469,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":448,"value":"function generateSGCoefficients(windowSize, polyOrder) {\n const halfWindow = Math.floor(windowSize / 2);\n let A = [];\n\n // Construct the design matrix A\n for (let i = -halfWindow; i <= halfWindow; i++) {\n let row = [];\n for (let j = 0; j <= polyOrder; j++) {\n row.push(Math.pow(i, j));\n }\n A.push(row);\n }\n\n // Compute (A^T A)^(-1) A^T e_0 using math.js\n let AT = math.transpose(A);\n let ATA = math.multiply(AT, A);\n let ATA_inv = math.inv(ATA); // Inverse of (A^T A)\n let e0 = new Array(windowSize).fill(0);\n //e0[Math.floor(windowSize / 2)] = 1; // Select the central point\n //let coeffs = math.multiply(ATA_inv, math.multiply(AT, e0));\n let coeffs = math.multiply(ATA_inv, AT)[0];\n\n return coeffs;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":893,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":845,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":424,"value":"### Continuous Fit with Fixed Knots","pinned":false,"hidden":false,"mode":"md"},{"id":767,"value":"fitFixedKnotsContinuous([0,1,2,3,4,5])","pinned":false,"hidden":false,"mode":"ojs"},{"id":298,"value":"$37","pinned":false,"hidden":false,"mode":"ojs"},{"id":327,"value":"$38","pinned":false,"hidden":false,"mode":"ojs"},{"id":406,"value":"### Discontinuous Fit with Automatically Pruned Knots","pinned":false,"hidden":false,"mode":"md"},{"id":377,"value":"$39","pinned":true,"hidden":false,"mode":"ojs"},{"id":740,"value":"fitPrunedKnotsDiscontinuous(d3.range(100, 200))","pinned":false,"hidden":false,"mode":"ojs"},{"id":209,"value":"function segmentFit(i1, i2, eps=1e-10) {\n // Calculate the sums over knots [i1,i2].\n const Sw = cumSums.Sw[i2] - cumSums.Sw[i1];\n const Sx = cumSums.Sx[i2] - cumSums.Sx[i1];\n const Sy = cumSums.Sy[i2] - cumSums.Sy[i1];\n const Sxx = cumSums.Sxx[i2] - cumSums.Sxx[i1];\n const Sxy = cumSums.Sxy[i2] - cumSums.Sxy[i1];\n const Syy = cumSums.Syy[i2] - cumSums.Syy[i1];\n // Calculate denominator for slope calculation.\n const denom = Sw * Sxx - Sx * Sx;\n // Calculate linear coefficients (a,b) and chisq for this segment.\n let a, b, chisq;\n if(denom > eps) {\n b = (Sw * Sxy - Sx * Sy) / denom;\n a = (Sy - b * Sx) / Sw;\n chisq = Syy - a * Sy - b * Sxy;\n }\n else {\n a = Sy / Sw;\n b = 0;\n chisq = Syy - a * Sy;\n }\n if(chisq < -1e-3) console.log(i1, i2, a, b, chisq);\n const xmid = (kwave[i1] + kwave[i2]) / 2;\n const amid = a + b * xmid;\n const linear = (x) => amid + b * (x - xmid);\n const F1 = linear(kwave[i1]), F2 = linear(kwave[i2]);\n \n return { a, b, chisq, linear, F1, F2 };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":544,"value":"cumSums = calculateCumSums(d3.range(kwave.length))","pinned":false,"hidden":false,"mode":"ojs"},{"id":204,"value":"function calculateCumSums(iknots, { ftype=Float64Array } = {}) {\n const n = iknots.length;\n // Calculate weighted sums of data between each knot.\n const Sw = new ftype(n).fill(0);\n const Sx = new ftype(n).fill(0);\n const Sy = new ftype(n).fill(0);\n const Sxx = new ftype(n).fill(0);\n const Sxy = new ftype(n).fill(0);\n const Syy = new ftype(n).fill(0);\n let idx = 1;\n const offset = index * arraySize;\n for(let i = 0; i < arraySize; ++i) {\n const w = ivar[offset + i], y = flux[offset + i], x = logwave[i];\n if(x > kwave[idx]) idx++;\n Sw[idx] += w;\n Sx[idx] += w * x;\n Sy[idx] += w * y;\n Sxx[idx] += w * x * x;\n Sxy[idx] += w * x * y;\n Syy[idx] += w * y * y;\n }\n // Calculate cumulative sums.\n for(let j = 2; j < n; ++j) {\n Sw[j] += Sw[j-1];\n Sx[j] += Sx[j-1];\n Sy[j] += Sy[j-1];\n Sxx[j] += Sxx[j-1];\n Sxy[j] += Sxy[j-1];\n Syy[j] += Syy[j-1];\n }\n return { Sw, Sx, Sy, Sxx, Sxy, Syy };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":413,"value":"### Select Possible Knot Locations","pinned":false,"hidden":false,"mode":"md"},{"id":306,"value":"breaks = {\n // for each kwave[i], find min j such that logwave[j] >= kwave[i]\n const breaks = [ ];\n let j = 0;\n for(let i = 0; i < kwave.length; ++i) {\n while(logwave[j] < kwave[i]) ++j;\n breaks.push(j);\n }\n breaks[breaks.length-1]++;\n return breaks;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":120,"value":"kwave = {\n const lo = logwave[0], hi = logwave[logwave.length-1], delta = (hi - lo) / nSegments;\n return d3.range(nSegments + 1).map(i => lo + i * delta);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":527,"value":"spacingPlot = {\n const kratio = (fullwave[arraySize-1] / fullwave[0]) ** (1 / nSegments);\n const dpix = fullwave.map(wlen => wlen * (kratio - 1) / pixsize);\n return Plot.plot({\n x: { label: \"Wavelength [A]\" },\n y: { label: `Knot Spacing [${pixsize}A pixels]` },\n grid: true,\n marks: [\n Plot.lineY(dpix, {x: (_,i) => fullwave[i], y: d => d}),\n Plot.frame(),\n ]\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":503,"value":"nSegments = oversampling * nSegmentsCoarse","pinned":false,"hidden":false,"mode":"ojs"},{"id":500,"value":"oversampling = 256","pinned":false,"hidden":false,"mode":"ojs"},{"id":497,"value":"nSegmentsCoarse = 8","pinned":false,"hidden":false,"mode":"ojs"},{"id":415,"value":"### Input Data","pinned":false,"hidden":false,"mode":"md"},{"id":33,"value":"spectrum = {\n const offset = index * arraySize;\n const data = d3.range(arraySize).map(i => {\n const w = ivar[offset + i];\n return {\n wlen: fullwave[i],\n flux: w > 0 ? flux[offset + i] : NaN,\n ivar: w,\n dflux: w > 0 ? w ** -0.5 : NaN,\n }\n });\n let [flo, fhi] = d3.extent(data, d => d.flux);\n flo = Math.min(flo, 0); // always include flux=0 in the range\n const fscaled = t => flo * (1 - t) + fhi * t;\n const tlo = (Math.min(0, d3.quantile(data, 0.001, d => d.flux)) - flo) / (fhi - flo);\n const thi = (d3.quantile(data, 0.999, d => d.flux) - flo) / (fhi - flo);\n return {fscaled, data, tlo, thi};\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1074,"value":"clipped = {\n const x = [ ], y = [ ], ivar2 = [ ];\n for(const {wlen,flux,ivar} of spectrum.data) {\n x.push(wlen.toFixed(1)*1);\n y.push(ivar > 0 ? flux.toFixed(3)*1 : NaN);\n ivar2.push(ivar.toFixed(4)*1);\n }\n return { x, y, ivar:ivar2 };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":228,"value":"nSpectra = {\n if(flux.length % arraySize != 0) throw new Error('ArraySize error');\n return Math.round(flux.length / arraySize);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":51,"value":"logwave = fullwave.map(w => Math.log(w))","pinned":false,"hidden":false,"mode":"ojs"},{"id":37,"value":"fullwave = Float64Array.from({length:arraySize}, (_,i) => (3600.0 + pixsize * i))","pinned":false,"hidden":false,"mode":"ojs"},{"id":521,"value":"pixsize = 0.8 // Angstrom","pinned":false,"hidden":false,"mode":"ojs"},{"id":517,"value":"arraySize = 7781","pinned":false,"hidden":false,"mode":"ojs"},{"id":18,"value":"flux = await FileAttachment(\"data.bin\").arrayBuffer().then(buf => new Float32Array(buf))","pinned":false,"hidden":false,"mode":"ojs"},{"id":27,"value":"ivar = await FileAttachment(\"ivar.bin\").arrayBuffer().then(buf => new Float32Array(buf))","pinned":false,"hidden":false,"mode":"ojs"},{"id":443,"value":"### Imports","pinned":false,"hidden":false,"mode":"md"},{"id":450,"value":"math = require(\"mathjs@14.3.0\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":134,"value":"import {interval} from '@mootari/range-slider'","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W3a","stdlib":"1","titling":"auto"}