Remez Algorithm - Source Code
This page displays the source code of key functions and methods used in the Remez algorithm implementation.
initializeCFRemezNodesRational(targetFunc, a, b, numNodes, numerDeg, denomDeg) {
try {
if (!Number.isFinite(a) || !Number.isFinite(b) || a >= b || numNodes < 2) {
throw new Error(`Invalid inputs: a=${a}, b=${b}, numNodes=${numNodes}`);
}
if (numerDeg < 0 || denomDeg < 0) {
throw new Error(`Invalid degrees: n=${numerDeg}, d=${denomDeg}`);
}
const n = numerDeg;
const d = denomDeg;
const totalDeg = n + d;
if (!Number.isFinite(totalDeg) || totalDeg < 0 || totalDeg > 1000) {
throw new Error(`Invalid total degree: ${totalDeg} (n=${n}, d=${d})`);
}
const numChebCoeffs = Math.min(2 * totalDeg + 200, 1500);
const numSamples = Math.max(40 * (totalDeg + 1), 5000);
console.log(`CF Rational: n=${n}, d=${d}, chebCoeffs=${numChebCoeffs}, samples=${numSamples}`);
const chebCoeffs = new Array(numChebCoeffs).fill(0);
for (let k = 0; k < numSamples; k++) {
const t = Math.cos(Math.PI * k / (numSamples - 1));
const x = (b - a) * (t + 1) / 2 + a;
const y = targetFunc(x);
if (!Number.isFinite(y)) {
console.warn(`Non-finite function value at x=${x}`);
continue;
}
for (let j = 0; j < numChebCoeffs && j < numSamples; j++) {
const Tj = this.chebyshevPolynomial(j, t);
chebCoeffs[j] += y * Tj;
}
}
for (let j = 0; j < chebCoeffs.length; j++) {
chebCoeffs[j] *= (j === 0) ? 1.0 / numSamples : 2.0 / numSamples;
if (Math.abs(chebCoeffs[j]) < 1e-14) chebCoeffs[j] = 0;
}
console.log('Chebyshev coefficients:', chebCoeffs.slice(0, Math.min(10, chebCoeffs.length)));
const cfResult = this.solveRationalCFSystem(chebCoeffs, n, d);
if (!cfResult.success) {
throw new Error('CF rational system solver failed');
}
const { numerCoeffs, denomCoeffs } = cfResult;
const cfRationalApprox = (x) => {
const t = 2 * (x - a) / (b - a) - 1;
let numer = 0;
for (let k = 0; k <= n && k < numerCoeffs.length; k++) {
numer += numerCoeffs[k] * this.chebyshevPolynomial(k, t);
}
let denom = 0;
for (let k = 0; k <= d && k < denomCoeffs.length; k++) {
denom += denomCoeffs[k] * this.chebyshevPolynomial(k, t);
}
if (Math.abs(denom) < 1e-12) {
console.warn(`Near-zero denominator at x=${x}, t=${t}`);
return numer / (Math.sign(denom) * 1e-12 || 1e-12);
}
return numer / denom;
};
console.log('CF rational - numerator coeffs:', numerCoeffs);
console.log('CF rational - denominator coeffs:', denomCoeffs);
const testError = Math.abs(targetFunc((a+b)/2) - cfRationalApprox((a+b)/2));
console.log(`CF test error at midpoint: ${testError.toExponential(3)}`);
const denseSamples = Math.max(1000, 100 * numNodes);
const extrema = this.findAlternatingExtremaDetailed(
targetFunc,
cfRationalApprox,
a, b,
numNodes,
denseSamples
);
const nodes = extrema.map(e => e.x);
const maxError = Math.max(...extrema.map(e => Math.abs(e.y)));
console.log(`CF rational: found ${nodes.length} extrema, maxError=${maxError}`);
if (nodes.length < numNodes) {
console.warn(`CF found only ${nodes.length} extrema, need ${numNodes}`);
return { nodes: nodes, success: true, error: 'Insufficient extrema' };
}
return {
numerCoeffs: numerCoeffs,
denomCoeffs: denomCoeffs,
nodes: nodes,
cfApprox: cfRationalApprox,
maxError: maxError,
success: true,
isRational: true,
degrees: { n, d }
};
} catch (error) {
console.error('CF rational approximation failed:', error);
return { nodes: [], success: false, error: error.message };
}
}initializeL2RemezNodes(targetFunc, a, b, numNodes) {
try {
const l2A = remMath.polynomialL2Approximation(targetFunc, numNodes-2, a, b, { silent: this.autoMode });
if (!l2A || !l2A.lsApprox) {
throw new Error('L2 approximation returned invalid result');
}
const lsApprox = l2A.lsApprox;
let extrema = this.findAlternatingExtrema(targetFunc, l2A.lsApprox, a, b, numNodes);
const nodes = extrema.map(e => e.x);
return {
coeffs: l2A.coeffs,
nodes: nodes,
lsApprox: l2A.lsApprox,
success: true
};
} catch (error) {
console.warn('L2 approximation failed:', error.message);
return { nodes: [], success: false, error: error.message };
}
}initializeDiscreteRemezNodes(targetFunc, a, b, numNodes, gridDensity = 100) {
const degree = numNodes - 2;
const numGridPoints = Math.max(1000, gridDensity * numNodes);
const gridPoints = Array.from({ length: numGridPoints }, (_, i) =>
a + i * (b - a) / (numGridPoints - 1)
);
const targetValues = gridPoints.map(x => targetFunc(x));
const discreteResult = this.solveDiscreteRemez(gridPoints, targetValues, degree);
if (!discreteResult.success) {
console.warn('Discrete Remez failed:', discreteResult.error);
return { nodes: [], success: false, error: discreteResult.error };
}
const extrema = this.findDiscreteExtremaFromSolution(
gridPoints,
targetValues,
discreteResult.approxValues,
numNodes
);
const nodes = extrema.map(e => e.x);
return {
coeffs: discreteResult.coeffs,
nodes: nodes,
discreteApprox: discreteResult.approxFunc,
maxError: Math.max(...extrema.map(e => Math.abs(e.error))),
success: true
};
}findReferenceSimple(errorFunction, coefficients, currentReference) {
const coarseSteps = parseInt(document.getElementById('bruteSteps').value);
const refinementSteps = 10;
const refinementIterations = 5;
const eps = 0.000001/n;
let newReference = [];
for (let k = 0; k < currentReference.length; k++) {
const isFirst = k == 0;
const isLast = k == currentReference.length - 1;
// If first/last point lies on the boundary, re-evaluate it there without searching
if (isFirst && currentReference[k].x <= aP.intervalStart + eps) {
const x = aP.intervalStart;
newReference.push({x, y: errorFunction(coefficients, x), originalIndex: k});
continue;
}
if (isLast && currentReference[k].x >= aP.intervalEnd - eps) {
const x = aP.intervalEnd;
newReference.push({x, y: errorFunction(coefficients, x), originalIndex: k});
continue;
}
// Free points: search within the interval between neighbors
const left = isFirst ? aP.intervalStart + eps : currentReference[k-1].x + eps;
const right = isLast ? aP.intervalEnd - eps : currentReference[k+1].x - eps;
if (left >= right || !isFinite(left) || !isFinite(right)) {
console.warn('Current reference: ', currentReference);
console.warn(`Invalid interval at k=${k}: [${left}, ${right}]`);
newReference.push({...currentReference[k], originalIndex: k});
continue;
}
let extremum = currentReference[k];
let delta = (right - left) / coarseSteps;
for (let step = 0; step <= coarseSteps; step++) {
const currentX = left + step * delta;
if (currentX > right) break;
const currentY = errorFunction(coefficients, currentX);
if (!isFinite(currentY)) continue;
if ((currentReference[k].y > 0 && currentY > extremum.y) ||
(currentReference[k].y < 0 && currentY < extremum.y)) {
extremum = {x: currentX, y: currentY};
}
}
for (let i = 0; i < refinementIterations; i++) {
const refinementLeft = Math.max(left, extremum.x - delta);
const refinementRight = Math.min(right, extremum.x + delta);
delta = (refinementRight - refinementLeft) / refinementSteps;
if (delta <= 0 || !isFinite(delta)) break;
for (let step = 0; step <= refinementSteps; step++) {
const currentX = refinementLeft + step * delta;
if (currentX > refinementRight) break;
const currentY = errorFunction(coefficients, currentX);
if (!isFinite(currentY)) continue;
if ((currentReference[k].y > 0 && currentY > extremum.y) ||
(currentReference[k].y < 0 && currentY < extremum.y)) {
extremum = {x: currentX, y: currentY};
}
}
}
newReference.push({...extremum, originalIndex: k});
}
return newReference;
}findReferenceParabolic(calculateError, optimizationCoefficients, currentReferencePoints) {
const maxInterpolationIterations = parseInt(document.getElementById('parabolicSteps').value);
const minimumPointSeparation = 0.01;
const convergenceThreshold = 1e-10;
const epsilon = 1e-12;
const goldenRatio = 0.618033988749;
// Enhanced parabolic interpolation with better numerical stability
function interpolateParabola(x0, x1, x2, y0, y1, y2) {
// Use more numerically stable form
const h1 = x1 - x0;
const h2 = x2 - x1;
const d1 = (y1 - y0) / h1;
const d2 = (y2 - y1) / h2;
const denominator = h1 + h2;
if (Math.abs(denominator) < epsilon || Math.abs(d2 - d1) < epsilon) {
return (x0 + x2) / 2; // Fallback to midpoint
}
const a = (d2 - d1) / denominator;
if (Math.abs(a) < epsilon) {
return (x0 + x2) / 2; // Nearly linear, use midpoint
}
// Compute vertex of parabola
const vertex = x1 - d1 / (2 * a) - h1 / 2;
return vertex;
}
// Adaptive bracket width calculation
function calculateInitialBracket(pointIndex, currentPoints) {
const baseWidth = 0.1; // Default bracket width
let adaptiveWidth = baseWidth;
if (pointIndex > 0) {
const leftSpacing = currentPoints[pointIndex].x - currentPoints[pointIndex - 1].x;
adaptiveWidth = Math.min(adaptiveWidth, leftSpacing * 0.3);
}
if (pointIndex < currentPoints.length - 1) {
const rightSpacing = currentPoints[pointIndex + 1].x - currentPoints[pointIndex].x;
adaptiveWidth = Math.min(adaptiveWidth, rightSpacing * 0.3);
}
return Math.max(adaptiveWidth, minimumPointSeparation * 2);
}
// Enhanced caching with LRU eviction
class ErrorCache {
constructor(maxSize = 1000) {
this.cache = new Map();
this.maxSize = maxSize;
}
get(x) {
const key = x.toFixed(12); // Round to avoid floating point issues
if (this.cache.has(key)) {
const value = this.cache.get(key);
this.cache.delete(key);
this.cache.set(key, value); // Move to end (most recent)
return value;
}
const value = calculateError(optimizationCoefficients, x);
if (this.cache.size >= this.maxSize) {
const firstKey = this.cache.keys().next().value;
this.cache.delete(firstKey);
}
this.cache.set(key, value);
return value;
}
}
const errorCache = new ErrorCache();
// Main optimization loop
let optimizedReferencePoints = [];
for (let pointIndex = 0; pointIndex < currentReferencePoints.length; pointIndex++) {
const currentPoint = currentReferencePoints[pointIndex];
// Calculate boundaries
const minAllowedX = (pointIndex === 0) ?
aP.intervalStart :
optimizedReferencePoints[pointIndex - 1].x + minimumPointSeparation;
const maxAllowedX = (pointIndex === currentReferencePoints.length - 1) ?
aP.intervalEnd :
currentReferencePoints[pointIndex + 1].x - minimumPointSeparation;
// Initialize with adaptive bracket
const bracketWidth = calculateInitialBracket(pointIndex, currentReferencePoints);
let leftX = Math.max(minAllowedX, currentPoint.x - bracketWidth);
let midX = currentPoint.x;
let rightX = Math.min(maxAllowedX, currentPoint.x + bracketWidth);
// Ensure proper ordering and spacing
if (rightX - leftX < 2 * minimumPointSeparation) {
const center = (leftX + rightX) / 2;
leftX = Math.max(minAllowedX, center - minimumPointSeparation);
rightX = Math.min(maxAllowedX, center + minimumPointSeparation);
midX = center;
}
let bestX = midX;
let bestY = errorCache.get(midX);
let iterationCount = 0;
let consecutiveSmallSteps = 0;
let lastStepSize = Infinity;
// Target direction: maximize absolute error
const targetSign = Math.sign(currentPoint.y);
while (iterationCount < maxInterpolationIterations) {
const leftY = errorCache.get(leftX);
const midY = errorCache.get(midX);
const rightY = errorCache.get(rightX);
// Try parabolic interpolation
let newX = interpolateParabola(leftX, midX, rightX,
leftY * targetSign, midY * targetSign, rightY * targetSign);
// Clamp to valid range
newX = Math.max(minAllowedX, Math.min(maxAllowedX, newX));
// Avoid getting stuck at boundaries
if (newX === minAllowedX || newX === maxAllowedX) {
newX = (leftX + rightX) / 2;
}
const stepSize = Math.abs(newX - bestX);
// Check convergence
if (stepSize < convergenceThreshold) {
consecutiveSmallSteps++;
if (consecutiveSmallSteps >= 3) break;
} else {
consecutiveSmallSteps = 0;
}
// Detect oscillation or stagnation
if (stepSize > lastStepSize * 0.9 && iterationCount > 3) {
// Switch to golden section if parabolic isn't converging
const range = rightX - leftX;
if (newX > midX) {
newX = midX + goldenRatio * (rightX - midX);
} else {
newX = midX - goldenRatio * (midX - leftX);
}
newX = Math.max(minAllowedX, Math.min(maxAllowedX, newX));
}
const newY = errorCache.get(newX);
// Update best point if we found improvement
if (Math.abs(newY) > Math.abs(bestY)) {
bestX = newX;
bestY = newY;
}
// Update bracket based on function behavior
if (newX < midX) {
if (Math.abs(newY) > Math.abs(midY)) {
rightX = midX;
midX = newX;
} else {
leftX = newX;
}
} else if (newX > midX) {
if (Math.abs(newY) > Math.abs(midY)) {
leftX = midX;
midX = newX;
} else {
rightX = newX;
}
}
// Ensure bracket doesn't collapse
if (rightX - leftX < 2 * convergenceThreshold) {
const expansion = Math.max(convergenceThreshold * 10, (maxAllowedX - minAllowedX) * 0.01);
const center = (leftX + rightX) / 2;
leftX = Math.max(minAllowedX, center - expansion);
rightX = Math.min(maxAllowedX, center + expansion);
midX = center;
}
lastStepSize = stepSize;
iterationCount++;
}
// Final adjustment for minimum separation
if (optimizedReferencePoints.length > 0) {
const minRequiredX = optimizedReferencePoints[optimizedReferencePoints.length - 1].x + minimumPointSeparation;
if (bestX < minRequiredX) {
bestX = minRequiredX;
bestY = errorCache.get(bestX);
}
}
optimizedReferencePoints.push({
x: bestX,
y: bestY,
iterations: iterationCount,
converged: consecutiveSmallSteps >= 3 || iterationCount < maxInterpolationIterations
});
}
console.log("Optimized Reference Points:", optimizedReferencePoints);
return optimizedReferencePoints;
}