gold 2/27/2026.
Here are some simple snippets for numerical methods. The goal is to use Tcl's minimalism as a learning tool. Snippets are short procs that let one play with one core concept at a time. All snippets are Playground V9 safe. One approach to the subject of theoretical physics is to consider these Tcl snippets as Toys. Some snippets here are listed as Toys. These Tcl procs are tiny entry points into physics. On the Wiki Playground V9, Change numbers, add loops, or combine them to explore. Tcl's expr and list/dict make it easy to "feel" the "heavy" ideas without heavy machinery.
The TCL Snippets illustrate ideal mathematical behavior only and do not perform full simulation, actual measurements, or state vector evolution. The tool only visualizes ideal math structure, whereas no state vector simulation, probabilities, or actual measurement outcomes are derived. This tool for visualization does not simulate actual measurement outcomes or state vector evolution during operations. These are idealized protocols for tutorial purposes. Primarily, TCL /TK uses its strong points here for book keeping and displays. The example tool is not a full emulator. Meaning, limited scope for tutorial purposes.
In debugging the calculations, some of the printout values reflect roughly 17-digit precision output from a typical double-precision computation. It's not "true exact" beyond 5 significant figures. Extra significant figures are used to check the calculations from other computer set-ups, not necessarily to infer accuracy of data measurements here. Typically, the slight differences in decimal places on far right of decimal point are normal floating-point behavior in Tcl's expr.
The page presents Tcl code snippets. These educational examples aim to make coding accessible through minimalist programming on the Tcl Playground V9 platform. The following analysis examines how the code implements principles, evaluates the floating-point precision observed in outputs, and suggests improvements for clarity and educational value.
This executive summary introduces the Whittaker-Eilers smoother, a signal-smoothing method grounded in penalised least squares (PLS), a mathematical technique that balances fitting accuracy against smoothness. The summary draws on Bowell's guide "The Perfect Way to Smooth Your Noisy Data," the pybaselines documentation, and Paul H. C. Eilers's foundational 2003 paper, "A Perfect Smoother," published in Analytical Chemistry. The goal is to explain what the method does, how to configure key parameters, and where practical limitations arise.
Edward Whittaker originally proposed the underlying framework roughly 80 years before Eilers formalised it in 2003. The smoother works by finding an output series that stays close to the original noisy measurements while also being as smooth as possible. The method balances two competing forces: fidelity to the original data and smoothness of the output. A concrete example is smoothing global temperature anomaly records from 1880 to 2022, where year-to-year noise obscures the long-term warming trend. Applying the smoother reveals the underlying trajectory clearly. Unlike the Savitzky-Golay (S-G) filter, which processes small windows of data sequentially, the Whittaker-Eilers smoother operates on the entire dataset at once, making the approach both globally consistent and computationally fast.
Two parameters govern behaviour. The first is lambda (λ), a scaling value that controls how much smoothness is penalised relative to data fidelity. A small lambda, such as 10, produces an output that closely tracks every fluctuation in the raw signal. A large lambda, such as 10,000,000, forces the output toward a near-straight line. Selecting the right lambda requires either visual inspection or an automated method called leave-one-out cross validation (LOOCV), which tests multiple lambda values and identifies the one producing the lowest prediction error. The second parameter is order (d), which determines how differences between adjacent points are measured. Order two is the default and suits most applications. Higher orders penalise curvature more aggressively and are appropriate for signals with complex, rapidly changing shapes.
Three specific drawbacks deserve attention. The first concerns serially correlated data, meaning data where each measurement is statistically related to the one before it, as is common in satellite time series or financial prices. Standard LOOCV assumes measurement errors are independent; serial correlation violates this assumption and can cause the method to under-smooth the data. A practical fix is to sample every fifth or tenth data point before running cross validation, which removes much of the correlation. The second drawback is lambda sensitivity to dataset size. A lambda value that fits a dataset of 100 points may need to increase by several orders of magnitude to fit the same underlying signal with 10,000 points. Whattoknow Practitioners who reuse a lambda value across datasets of different lengths may produce inconsistent results. The third drawback is that the basic formulation assumes evenly spaced data; unevenly spaced inputs require additional configuration before the smoother can be applied reliably.
The Whittaker-Eilers smoother offers speed, built-in interpolation for missing data, and a single intuitive parameter in lambda, making the method well suited to large or gappy datasets. Practitioners should use cross validation to select lambda rather than guessing, adjust expectations when data is serially correlated, and rescale lambda when dataset size changes substantially
Gist.
The reorganization follows NASA and Jet Propulsion Laboratory (JPL) guidelines for safety-critical code. Procedures stay under 25 lines, include at least two assertions for validation, use descriptive variable names.
gold 2/28/2026.
Maintain single-purpose TCL procedures under 30 lines each. Test TCL invariants after every transformation. Use descriptive names embedding physics meaning in the TCL code. Convert tabs to spaces uniformly. Print variable states at computation boundaries. These steps transform debugging into verification process reliably. Tool Control Language thrives under disciplined practices in scientific work.
The program demonstrates math patterns.
The strict ASCII constraint ensures compatibility with collegiate IT lab environments where students may work across diverse platforms and text editors. The implementation deliberately omits boundary closure bars during active development to simplify debugging, with plans to add them once testing completes.
Note. These Snippets on Theoretical Physics are a set, not stand alones. Recommend read all of the set.
This is a draft, still debugging on Playground V9. convert to strict 7-bit ASCII for Playground V9.
# Whittaker Eilers Digital Smoothing Filter - VERSION V2
# Tcl/Tk 8.6+ 7-bit ASCII safe. NASA/JPL defensive programming style.
# NASA/JPL defensive programming style.
# Compatible with Tcl/Tk (Tool Command Language / Toolkit) 8.6+
# Written for Windows 11 on ActiveState Tcl.
# Working under strict 7-bit ASCII encoding.
# Optimized for collegiate information technology lab environments.
# TCL club, 02/27/2026
# This preserves constant and linear signals without amplitude distortion.
#
# Tcl/Tk (Tool Command Language / Toolkit) 8.6+ 7-bit ASCII safe.
# NASA/JPL defensive programming style.
# Compatible with Windows 11 on ActiveState Tcl.
#
# WHAT THIS PROGRAM DOES:
# Applies the Whittaker-Eilers smoothing algorithm to a list of noisy
# numeric values. The smoother solves one linear system over the entire
# dataset at once, unlike the Savitzky-Golay filter which slides a window
# point by point.
#
# MATHEMATICAL PRINCIPLE:
# Find smoothed values z that minimize:
# S = sum( (y_i - z_i)^2 ) + lambda * sum( (Delta^d z_i)^2 )
# The first term penalizes departure from the raw data y.
# The second term penalizes roughness in the smoothed result z.
# Delta^d means the d-th order finite difference operator.
# lambda (lambdaPenalty) controls the trade-off:
# lambda = 0 -> output equals input exactly (no smoothing at all)
# lambda large -> output approaches a polynomial of degree (d-1)
#
# The minimizer satisfies the linear system:
# ( I + lambda * D_transposed * D ) * z = y
# where D is the finite difference matrix of order d,
# I is the identity matrix, and D_transposed is the transpose of D.
#
# KEY PROPERTIES (verified by autotests below):
# 1. A constant input signal passes through unchanged for any lambda.
# 2. A linear ramp passes through unchanged for any lambda with order >= 2.
# 3. Setting lambda = 0 always returns the input unchanged.
# 4. A single peak is smoothed and spread symmetrically for order 2.
#
# DIFFERENCE FROM SAVITZKY-GOLAY:
# Savitzky-Golay fits local polynomial windows and trims edge points.
# Whittaker-Eilers keeps all points and handles edges without special cases.
#
# PARAMETERS:
# inputDataList - Tcl list of numeric values to smooth (must be evenly spaced)
# lambdaPenalty - Non-negative smoothing strength (0.0 to millions)
# differenceOrder - Degree of difference penalty (1 = first differences,
# 2 = second differences, 2 is standard default)
#
# LIMITATIONS:
# Uses dense Gaussian elimination: O(n^3) operations.
# Suitable for lab datasets up to a few hundred points.
# For very large datasets a banded or sparse solver would be faster.
# Assumes data are evenly spaced in time or space.
#
# REFERENCES:
# Eilers, P.H.C. (2003). "A Perfect Smoother."
# Analytical Chemistry, 75(14), 3631-3636.
# Whittaker, E.T. (1923). "On a New Method of Graduation."
# Proceedings of the Edinburgh Mathematical Society, 41, 63-75.
# ===========================================================================
console show
# ---------------------------------------------------------------------------
# assertConditionIsTrue
# Purpose: Defensive assertion that halts execution with a clear message
# if the given condition is false. Follows NASA/JPL Rule 7: all code
# must check return values and assert pre/post conditions explicitly.
# ---------------------------------------------------------------------------
proc assertConditionIsTrue {conditionBoolean failureMessage} {
if {!$conditionBoolean} {
error "ASSERTION FAILED: $failureMessage"
}
}
# ---------------------------------------------------------------------------
# computeBinomialCoefficient
# Purpose: Compute the binomial coefficient C(topValue, bottomValue),
# also written "topValue choose bottomValue".
# Formula: topValue! / ( bottomValue! * (topValue - bottomValue)! )
# Used to build the coefficients of finite difference rows.
#
# Example: C(2,0)=1, C(2,1)=2, C(2,2)=1 (these are Pascal triangle row 2)
#
# Implementation uses the multiplicative formula to stay in integers:
# C(n,k) = product from step=0 to k-1 of (n - step) / (step + 1)
# Integer division is exact at each step due to combinatorial properties.
# ---------------------------------------------------------------------------
proc computeBinomialCoefficient {topValue bottomValue} {
assertConditionIsTrue \
[expr {$bottomValue >= 0 && $bottomValue <= $topValue}] \
"binomial coefficient requires 0 <= bottomValue <= topValue"
set runningProduct 1
for {set stepIndex 0} {$stepIndex < $bottomValue} {incr stepIndex} {
set runningProduct [expr {$runningProduct * ($topValue - $stepIndex) / ($stepIndex + 1)}]
}
return $runningProduct
}
# ---------------------------------------------------------------------------
# buildFiniteDifferenceMatrix
# Purpose: Build the finite difference matrix D of the requested order.
# The matrix has (dataLength - differenceOrder) rows and dataLength columns.
# Row i contains the coefficients for computing the i-th finite difference.
#
# General coefficient rule for row i, column position (i + coeffStep):
# coefficient = (-1)^(differenceOrder - coeffStep) * C(differenceOrder, coeffStep)
#
# Order 1 example, row i: [ -1, 1, 0, ... ] starting at column i
# Computes: y[i+1] - y[i]
# Order 2 example, row i: [ 1, -2, 1, 0, ... ] starting at column i
# Computes: y[i] - 2*y[i+1] + y[i+2]
#
# Second-order differences of a constant or linear ramp equal zero.
# That is why constants and ramps pass through the smoother unchanged.
# ---------------------------------------------------------------------------
proc buildFiniteDifferenceMatrix {dataLength differenceOrder} {
assertConditionIsTrue \
[expr {$dataLength > $differenceOrder}] \
"dataLength must exceed differenceOrder to form a valid difference matrix"
assertConditionIsTrue \
[expr {$differenceOrder >= 1}] \
"differenceOrder must be at least 1"
set numberOfDifferenceRows [expr {$dataLength - $differenceOrder}]
set numberOfMatrixCols $dataLength
# Precompute the coefficient values for a single row.
# There are (differenceOrder + 1) nonzero values per row.
set rowCoefficientPattern {}
for {set coeffStep 0} {$coeffStep <= $differenceOrder} {incr coeffStep} {
set binomialValue [computeBinomialCoefficient $differenceOrder $coeffStep]
set signExponent [expr {$differenceOrder - $coeffStep}]
# (-1)^signExponent: positive when signExponent is even, negative when odd
if {$signExponent % 2 == 0} {
set signMultiplier 1
} else {
set signMultiplier -1
}
lappend rowCoefficientPattern [expr {$signMultiplier * $binomialValue}]
}
# Initialize the full matrix as all zeros.
set differenceMatrix {}
for {set rowIndex 0} {$rowIndex < $numberOfDifferenceRows} {incr rowIndex} {
lappend differenceMatrix [lrepeat $numberOfMatrixCols 0]
}
# Place the coefficient pattern into each row starting at the diagonal.
for {set rowIndex 0} {$rowIndex < $numberOfDifferenceRows} {incr rowIndex} {
for {set coeffStep 0} {$coeffStep <= $differenceOrder} {incr coeffStep} {
set targetColIndex [expr {$rowIndex + $coeffStep}]
set coeffValue [lindex $rowCoefficientPattern $coeffStep]
lset differenceMatrix $rowIndex $targetColIndex $coeffValue
}
}
assertConditionIsTrue \
[expr {[llength $differenceMatrix] == $numberOfDifferenceRows}] \
"difference matrix row count does not match expected numberOfDifferenceRows"
return $differenceMatrix
}
# ---------------------------------------------------------------------------
# multiplyMatrixTransposeBySelf
# Purpose: Compute M_transposed * M for a matrix M with given dimensions.
# inputMatrix has inputRowCount rows and inputColCount columns.
# The result is a square matrix of size inputColCount x inputColCount.
#
# Formula: result[colA][colB] = sum over all rows of M[row][colA] * M[row][colB]
#
# The result is always symmetric: result[colA][colB] = result[colB][colA].
# The loop exploits symmetry by computing only the upper triangle and
# then copying to the lower triangle, halving the work.
#
# In Whittaker-Eilers usage: inputMatrix is D, result is D_transposed * D.
# D_transposed * D captures the roughness structure of the data.
# ---------------------------------------------------------------------------
proc multiplyMatrixTransposeBySelf {inputMatrix inputRowCount inputColCount} {
assertConditionIsTrue \
[expr {$inputRowCount > 0 && $inputColCount > 0}] \
"matrix must have positive row and column counts"
# Initialize result as a square zero matrix of size inputColCount.
set resultMatrix {}
for {set rowIndex 0} {$rowIndex < $inputColCount} {incr rowIndex} {
lappend resultMatrix [lrepeat $inputColCount 0.0]
}
# Compute upper triangle (including diagonal) and mirror to lower triangle.
for {set colA 0} {$colA < $inputColCount} {incr colA} {
for {set colB $colA} {$colB < $inputColCount} {incr colB} {
set dotProductAccumulator 0.0
for {set sharedRowIndex 0} {$sharedRowIndex < $inputRowCount} {incr sharedRowIndex} {
set valueFromColA [lindex [lindex $inputMatrix $sharedRowIndex] $colA]
set valueFromColB [lindex [lindex $inputMatrix $sharedRowIndex] $colB]
set dotProductAccumulator \
[expr {$dotProductAccumulator + $valueFromColA * $valueFromColB}]
}
lset resultMatrix $colA $colB $dotProductAccumulator
lset resultMatrix $colB $colA $dotProductAccumulator
}
}
return $resultMatrix
}
# ---------------------------------------------------------------------------
# buildPenalizedSystemMatrix
# Purpose: Construct the matrix A = I + lambda * D_transposed * D.
# Solving the linear system A * smoothedValues = inputValues
# gives the Whittaker-Eilers optimal smoothed signal.
#
# I is the identity matrix (ones on diagonal, zeros elsewhere).
# D is the finite difference matrix of the requested order.
# lambda is the smoothing strength parameter.
#
# When lambda = 0: A = I, so solution = input (no smoothing).
# When lambda is large: solution approaches a polynomial of degree (order-1).
# ---------------------------------------------------------------------------
proc buildPenalizedSystemMatrix {dataLength lambdaPenalty differenceOrder} {
set differenceMatrix [buildFiniteDifferenceMatrix $dataLength $differenceOrder]
set numberOfDifferenceRows [expr {$dataLength - $differenceOrder}]
set penaltyMatrix [multiplyMatrixTransposeBySelf \
$differenceMatrix $numberOfDifferenceRows $dataLength]
# Build A = I + lambda * penaltyMatrix element by element.
set systemMatrix {}
for {set rowIndex 0} {$rowIndex < $dataLength} {incr rowIndex} {
set newRow {}
for {set colIndex 0} {$colIndex < $dataLength} {incr colIndex} {
set penaltyEntry [lindex [lindex $penaltyMatrix $rowIndex] $colIndex]
# Identity matrix contributes 1.0 on the diagonal, 0.0 elsewhere.
if {$rowIndex == $colIndex} {
set identityEntry 1.0
} else {
set identityEntry 0.0
}
lappend newRow [expr {$identityEntry + $lambdaPenalty * $penaltyEntry}]
}
lappend systemMatrix $newRow
}
return $systemMatrix
}
# ---------------------------------------------------------------------------
# solveLinearSystemByGaussianElimination
# Purpose: Solve the linear system coefficientMatrix * solutionVector = rhsVector
# using Gaussian elimination with partial pivoting.
#
# Partial pivoting: at each step, swap rows to place the largest available
# coefficient on the diagonal before eliminating. This avoids division by
# very small numbers and improves numerical stability.
#
# The procedure works on internal copies and does not modify the inputs.
# Returns the solution as a Tcl list of systemDimension floating-point values.
#
# Complexity: O(n^3) for n = systemDimension. Acceptable for lab datasets.
# ---------------------------------------------------------------------------
proc solveLinearSystemByGaussianElimination \
{coefficientMatrix rhsVector systemDimension} {
assertConditionIsTrue \
[expr {$systemDimension >= 1}] \
"systemDimension must be at least 1"
assertConditionIsTrue \
[expr {[llength $rhsVector] == $systemDimension}] \
"rhsVector length must equal systemDimension"
# Working copies so the caller's data is not changed.
set workingMatrix $coefficientMatrix
set workingRhs $rhsVector
# ---- Forward Elimination with Partial Pivoting -------------------------
for {set pivotRowIndex 0} {$pivotRowIndex < $systemDimension} {incr pivotRowIndex} {
# Search for the row with the largest absolute value in the pivot column.
set largestAbsoluteValue \
[expr {abs([lindex [lindex $workingMatrix $pivotRowIndex] $pivotRowIndex])}]
set largestAbsoluteRowIndex $pivotRowIndex
for {set searchRowIndex [expr {$pivotRowIndex + 1}]} \
{$searchRowIndex < $systemDimension} \
{incr searchRowIndex} {
set candidateAbsoluteValue \
[expr {abs([lindex [lindex $workingMatrix $searchRowIndex] $pivotRowIndex])}]
if {$candidateAbsoluteValue > $largestAbsoluteValue} {
set largestAbsoluteValue $candidateAbsoluteValue
set largestAbsoluteRowIndex $searchRowIndex
}
}
# Swap the current pivot row with the row holding the largest value.
if {$largestAbsoluteRowIndex != $pivotRowIndex} {
set savedRow [lindex $workingMatrix $pivotRowIndex]
lset workingMatrix $pivotRowIndex \
[lindex $workingMatrix $largestAbsoluteRowIndex]
lset workingMatrix $largestAbsoluteRowIndex $savedRow
set savedRhsValue [lindex $workingRhs $pivotRowIndex]
lset workingRhs $pivotRowIndex \
[lindex $workingRhs $largestAbsoluteRowIndex]
lset workingRhs $largestAbsoluteRowIndex $savedRhsValue
}
set pivotDiagonalValue \
[lindex [lindex $workingMatrix $pivotRowIndex] $pivotRowIndex]
assertConditionIsTrue \
[expr {abs($pivotDiagonalValue) > 1.0e-14}] \
"singular or near-singular matrix at pivot row $pivotRowIndex - check lambda and data"
# Eliminate all rows below the current pivot row.
for {set targetRowIndex [expr {$pivotRowIndex + 1}]} \
{$targetRowIndex < $systemDimension} \
{incr targetRowIndex} {
set eliminationMultiplier \
[expr {[lindex [lindex $workingMatrix $targetRowIndex] $pivotRowIndex] \
/ $pivotDiagonalValue}]
for {set colIndex $pivotRowIndex} \
{$colIndex < $systemDimension} \
{incr colIndex} {
set targetEntry \
[lindex [lindex $workingMatrix $targetRowIndex] $colIndex]
set pivotEntry \
[lindex [lindex $workingMatrix $pivotRowIndex] $colIndex]
lset workingMatrix $targetRowIndex $colIndex \
[expr {$targetEntry - $eliminationMultiplier * $pivotEntry}]
}
set targetRhsValue [lindex $workingRhs $targetRowIndex]
set pivotRhsValue [lindex $workingRhs $pivotRowIndex]
lset workingRhs $targetRowIndex \
[expr {$targetRhsValue - $eliminationMultiplier * $pivotRhsValue}]
}
}
# ---- Back Substitution ------------------------------------------------
# The matrix is now upper triangular. Solve from the last row upward.
set solutionVector [lrepeat $systemDimension 0.0]
for {set solveRowIndex [expr {$systemDimension - 1}]} \
{$solveRowIndex >= 0} \
{incr solveRowIndex -1} {
set accumulatedRhsValue [lindex $workingRhs $solveRowIndex]
for {set alreadySolvedCol [expr {$solveRowIndex + 1}]} \
{$alreadySolvedCol < $systemDimension} \
{incr alreadySolvedCol} {
set matrixEntry \
[lindex [lindex $workingMatrix $solveRowIndex] $alreadySolvedCol]
set knownSolutionValue [lindex $solutionVector $alreadySolvedCol]
set accumulatedRhsValue \
[expr {$accumulatedRhsValue - $matrixEntry * $knownSolutionValue}]
}
set diagonalEntry \
[lindex [lindex $workingMatrix $solveRowIndex] $solveRowIndex]
lset solutionVector $solveRowIndex \
[expr {$accumulatedRhsValue / $diagonalEntry}]
}
return $solutionVector
}
# ---------------------------------------------------------------------------
# smoothDataWithWhittakerEilers
# Purpose: Main entry point. Apply Whittaker-Eilers smoothing to a numeric
# data list and return a smoothed list of the same length.
#
# Arguments:
# inputDataList - list of numeric values (evenly spaced measurements)
# lambdaPenalty - smoothing strength: 0.0 = no change, 1e6 = very smooth
# differenceOrder - roughness order: 1 = first differences, 2 = second (default)
#
# Returns:
# A Tcl list of smoothed floating-point values, same length as inputDataList.
# ---------------------------------------------------------------------------
proc smoothDataWithWhittakerEilers {inputDataList lambdaPenalty differenceOrder} {
set dataLength [llength $inputDataList]
assertConditionIsTrue \
[expr {$dataLength > $differenceOrder}] \
"data list length must exceed differenceOrder"
assertConditionIsTrue \
[expr {$lambdaPenalty >= 0.0}] \
"lambdaPenalty must be non-negative (zero means no smoothing)"
assertConditionIsTrue \
[expr {$differenceOrder >= 1}] \
"differenceOrder must be at least 1"
set systemMatrix [buildPenalizedSystemMatrix \
$dataLength $lambdaPenalty $differenceOrder]
# Convert input list to a floating-point right-hand-side vector.
set rhsVector {}
foreach rawDataValue $inputDataList {
lappend rhsVector [expr {double($rawDataValue)}]
}
set smoothedValues [solveLinearSystemByGaussianElimination \
$systemMatrix $rhsVector $dataLength]
assertConditionIsTrue \
[expr {[llength $smoothedValues] == $dataLength}] \
"smoother must return exactly as many values as the input"
return $smoothedValues
}
# ===========================================================================
# TOLERANCE AND COMPARISON HELPERS
# ===========================================================================
# ---------------------------------------------------------------------------
# numericAbsoluteDifference
# Purpose: Safe absolute-value of (valueA minus valueB) without the Tcl
# abs() function to remain explicit and portable.
# ---------------------------------------------------------------------------
proc numericAbsoluteDifference {valueA valueB} {
set rawDifference [expr {$valueA - $valueB}]
if {$rawDifference < 0.0} {
return [expr {0.0 - $rawDifference}]
}
return $rawDifference
}
# ---------------------------------------------------------------------------
# numericApproximatelyEqualWithMargin
# Purpose: Return 1 if two floating-point values agree within a combined
# relative and absolute tolerance, 0 otherwise.
# Tolerance = absoluteMargin + relativeMargin * max(1.0, |expectedValue|)
# This avoids false failures when expected values are very large or very small.
# ---------------------------------------------------------------------------
proc numericApproximatelyEqualWithMargin \
{expectedValue actualValue relativeMargin absoluteMargin} {
set absoluteDifference [numericAbsoluteDifference $expectedValue $actualValue]
if {abs($expectedValue) < 1.0} {
set referenceMagnitude 1.0
} else {
set referenceMagnitude [expr {abs($expectedValue)}]
}
set allowedTolerance \
[expr {$absoluteMargin + $relativeMargin * $referenceMagnitude}]
if {$absoluteDifference <= $allowedTolerance} {
return 1
}
return 0
}
# ---------------------------------------------------------------------------
# listApproximatelyEqualWithMargin
# Purpose: Element-wise tolerance test across two numeric lists.
# Returns 1 if every corresponding pair passes numericApproximatelyEqualWithMargin.
# Returns 0 if lists differ in length or any element pair fails.
# ---------------------------------------------------------------------------
proc listApproximatelyEqualWithMargin \
{expectedList actualList relativeMargin absoluteMargin} {
if {[llength $expectedList] != [llength $actualList]} {
return 0
}
foreach expectedValue $expectedList actualValue $actualList {
set elementPasses [numericApproximatelyEqualWithMargin \
$expectedValue $actualValue $relativeMargin $absoluteMargin]
if {!$elementPasses} {
return 0
}
}
return 1
}
# ---------------------------------------------------------------------------
# runSingleWhittakerEilersAutotest
# Purpose: Execute one named autotest, print inputs and outputs, and report
# PASS or FAIL based on the tolerance comparison.
# ---------------------------------------------------------------------------
proc runSingleWhittakerEilersAutotest \
{testName inputDataList lambdaPenalty differenceOrder \
expectedOutputList relativeMargin absoluteMargin} {
puts "---- AUTOTEST: $testName ----"
puts "input: $inputDataList"
puts "lambda: $lambdaPenalty"
puts "order: $differenceOrder"
puts "expected: $expectedOutputList"
set actualOutputList [smoothDataWithWhittakerEilers \
$inputDataList $lambdaPenalty $differenceOrder]
puts "actual: $actualOutputList"
set passFlag [listApproximatelyEqualWithMargin \
$expectedOutputList $actualOutputList $relativeMargin $absoluteMargin]
if {$passFlag} {
puts "RESULT: PASS 1"
} else {
puts "RESULT: FAIL 0"
}
puts "-----------------------------"
}
# ===========================================================================
# AUTOTEST SUITE
# Tolerance: 1% relative + 0.0001 absolute (tighter than SG suite because
# the solver produces near-exact answers for these structured inputs).
#
# MATHEMATICAL BASIS FOR EXPECTED VALUES:
#
# Property A (constant preservation):
# D^2 of a constant vector = 0.
# Therefore D_transposed * D * constant = 0.
# Therefore (I + lambda * D_transposed * D) * constant = constant.
# Output equals input regardless of lambda.
#
# Property B (linear ramp preservation with order 2):
# D^2 of a linear ramp = 0 (second differences of a straight line are zero).
# Same reasoning as Property A: output equals input regardless of lambda.
#
# Property C (lambda = 0 gives identity):
# System matrix = I + 0 * D_transposed * D = I.
# Solution of I * z = y is z = y.
# Output equals input regardless of data shape.
#
# Property D (single peak, n=5, lambda=1, order=2):
# Hand-verified by solving the 5x5 linear system algebraically.
# Input {0 0 10 0 0} -> output {5/12, 5/2, 25/6, 5/2, 5/12}
# = {0.41667, 2.5, 4.16667, 2.5, 0.41667}
# Symmetry confirmed: z[0]=z[4], z[1]=z[3], peak reduced from 10 to 4.167.
#
# Property E (single peak, n=3, lambda=1, order=2):
# Minimum-size dataset for order-2 differences (needs at least 3 points).
# Hand-solved 3x3 system: input {0 1 0} -> output {2/7, 3/7, 2/7}
# = {0.28571, 0.42857, 0.28571}
# ===========================================================================
set autotestRelativeTolerance 0.01
set autotestAbsoluteTolerance 0.0001
# ---------------------------------------------------------------------------
# Autotest 1: CONSTANT INPUT
# Keywords: constant signal, any lambda, output unchanged, Property A
# Input: ten identical values of 5.0
# Lambda: 1000.0 (strong smoothing)
# Order: 2 (standard second-difference penalty)
# Expected: ten identical values of 5.0
# Reason: Second differences of a constant are zero; penalty has no effect.
# ---------------------------------------------------------------------------
runSingleWhittakerEilersAutotest "constant_input_lambda1000_order2" \
{5 5 5 5 5 5 5 5 5 5} \
1000.0 \
2 \
{5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0} \
$autotestRelativeTolerance $autotestAbsoluteTolerance
# ---------------------------------------------------------------------------
# Autotest 2: LINEAR RAMP INPUT
# Keywords: linear ramp, strong lambda, output unchanged, Property B
# Input: 1 through 10 (arithmetic progression, step 1)
# Lambda: 500.0 (strong smoothing)
# Order: 2 (second-difference penalty)
# Expected: 1 through 10 unchanged
# Reason: Second differences of a straight line are zero; penalty has no effect.
# ---------------------------------------------------------------------------
runSingleWhittakerEilersAutotest "linear_ramp_lambda500_order2" \
{1 2 3 4 5 6 7 8 9 10} \
500.0 \
2 \
{1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0} \
$autotestRelativeTolerance $autotestAbsoluteTolerance
# ---------------------------------------------------------------------------
# Autotest 3: ZERO LAMBDA (IDENTITY BEHAVIOR)
# Keywords: lambda zero, identity, arbitrary input, no smoothing, Property C
# Input: irregular values {3 7 1 9 4 6 2 8}
# Lambda: 0.0 (no penalty applied)
# Order: 2
# Expected: output equals input exactly
# Reason: System matrix = I when lambda = 0; solution of I*z = y is z = y.
# ---------------------------------------------------------------------------
runSingleWhittakerEilersAutotest "zero_lambda_identity_order2" \
{3 7 1 9 4 6 2 8} \
0.0 \
2 \
{3.0 7.0 1.0 9.0 4.0 6.0 2.0 8.0} \
$autotestRelativeTolerance $autotestAbsoluteTolerance
# ---------------------------------------------------------------------------
# Autotest 4: SINGLE PEAK, N=5, LAMBDA=1, ORDER=2
# Keywords: single peak, symmetric smoothing, hand-verified algebra, Property D
# Input: {0 0 10 0 0} -- spike of height 10 at center position
# Lambda: 1.0 (light smoothing, balances fidelity and roughness equally)
# Order: 2
# Expected: {5/12, 5/2, 25/6, 5/2, 5/12}
# = {0.41667, 2.5, 4.16667, 2.5, 0.41667}
# Reason: Symmetric input produces symmetric output.
# Peak reduced from 10 to 25/6 = 4.167 and energy spread to neighbors.
# The 5x5 linear system was solved algebraically:
# Let a = z[0] = z[4], b = z[1] = z[3], c = z[2].
# Row 0: 2a - 2b + c = 0
# Row 1: -2a + 7b - 4c = 0 (using symmetry z[3]=b)
# Row 2: 2a - 8b + 7c = 10
# Solving: b = 6a, c = 10a, then 24a = 10, so a = 5/12.
# ---------------------------------------------------------------------------
runSingleWhittakerEilersAutotest "single_peak_n5_lambda1_order2" \
{0 0 10 0 0} \
1.0 \
2 \
{0.41666667 2.5 4.16666667 2.5 0.41666667} \
$autotestRelativeTolerance $autotestAbsoluteTolerance
# ---------------------------------------------------------------------------
# Autotest 5: SINGLE PEAK, N=3, LAMBDA=1, ORDER=2 (MINIMUM VALID SIZE)
# Keywords: minimum dataset, n=3, hand-solved 3x3 system, Property E
# Input: {0 1 0} -- unit spike at center of three-point dataset
# Lambda: 1.0
# Order: 2 (requires n > 2; this is the smallest valid input for order 2)
# Expected: {2/7, 3/7, 2/7} = {0.28571, 0.42857, 0.28571}
# Reason: D_transposed * D for n=3 is a 3x3 matrix with structure:
# [ 1 -2 1 ]
# [-2 4 -2 ]
# [ 1 -2 1 ]
# System A = I + D_transposed*D:
# [ 2 -2 1 ]
# [-2 5 -2 ]
# [ 1 -2 2 ]
# Solving with y = {0, 1, 0} gives z[0]=z[2]=2/7, z[1]=3/7.
# Peak reduced from 1.0 to 3/7 = 0.4286 with symmetric wing spread.
# ---------------------------------------------------------------------------
runSingleWhittakerEilersAutotest "single_peak_n3_lambda1_order2_minimum_size" \
{0 1 0} \
1.0 \
2 \
{0.28571429 0.42857143 0.28571429} \
$autotestRelativeTolerance $autotestAbsoluteTolerance
# End of deck
---- AUTOTEST: constant_input_lambda1000_order2 ---- input: 5 5 5 5 5 5 5 5 5 5 lambda: 1000.0 order: 2 expected: 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 actual: 5.000000000000498 5.000000000000405 5.000000000000313 5.000000000000223 5.000000000000135 5.000000000000051 4.999999999999969 4.999999999999888 4.999999999999807 4.999999999999726 RESULT: PASS 1 ----------------------------- ---- AUTOTEST: linear_ramp_lambda500_order2 ---- input: 1 2 3 4 5 6 7 8 9 10 lambda: 500.0 order: 2 expected: 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 actual: 1.00000000000012 2.0000000000000893 3.0000000000000577 4.000000000000025 4.99999999999999 5.999999999999954 6.9999999999999165 7.999999999999878 8.99999999999984 9.999999999999803 RESULT: PASS 1 ----------------------------- ---- AUTOTEST: zero_lambda_identity_order2 ---- input: 3 7 1 9 4 6 2 8 lambda: 0.0 order: 2 expected: 3.0 7.0 1.0 9.0 4.0 6.0 2.0 8.0 actual: 3.0 7.0 1.0 9.0 4.0 6.0 2.0 8.0 RESULT: PASS 1 ----------------------------- ---- AUTOTEST: single_peak_n5_lambda1_order2 ---- input: 0 0 10 0 0 lambda: 1.0 order: 2 expected: 0.41666667 2.5 4.16666667 2.5 0.41666667 actual: 0.4166666666666665 2.4999999999999996 4.166666666666666 2.4999999999999996 0.4166666666666666 RESULT: PASS 1 ----------------------------- ---- AUTOTEST: single_peak_n3_lambda1_order2_minimum_size ---- input: 0 1 0 lambda: 1.0 order: 2 expected: 0.28571429 0.42857143 0.28571429 actual: 0.2857142857142857 0.42857142857142855 0.2857142857142857 RESULT: PASS 1 ----------------------------- (bin) 1 %
Note. Program or algorithm is working on Playground V9. But troubles in grabbing the extra long output file from Playground here. This output is from a reduced deck here.
The Playground V9 session was dropping the proc definitions because they were entered separately from the autotest calls.
"-----------------------------"
> }
(tcl) 12 % set autotestRelativeTolerance 0.01
0.01
(tcl) 13 % set autotestAbsoluteTolerance 0.0001
0.0001
(tcl) 14 % runSingleWhittakerEilersAutotest "constant_input_lambda1000_order2" \
> {5 5 5 5 5 5 5 5 5 5} \
> 1000.0 2 \
> {5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0} \
> $autotestRelativeTolerance $autotestAbsoluteTolerance
---- AUTOTEST: constant_input_lambda1000_order2 ----
input: 5 5 5 5 5 5 5 5 5 5
lambda: 1000.0
order: 2
expected: 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0
actual: 5.000000000000498 5.000000000000405 5.000000000000313 5.000000000000223 5.000000000000135 5.000000000000051 4.999999999999969 4.999999999999888 4.999999999999807 4.999999999999726
RESULT: PASS 1
-----------------------------
(tcl) 15 % runSingleWhittakerEilersAutotest "linear_ramp_lambda500_order2" \
> {1 2 3 4 5 6 7 8 9 10} \
> 500.0 2 \
> {1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0} \
> $autotestRelativeTolerance $autotestAbsoluteTolerance
---- AUTOTEST: linear_ramp_lambda500_order2 ----
input: 1 2 3 4 5 6 7 8 9 10
lambda: 500.0
order: 2
expected: 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0
actual: 1.00000000000012 2.0000000000000893 3.0000000000000577 4.000000000000025 4.99999999999999 5.999999999999954 6.9999999999999165 7.999999999999878 8.99999999999984 9.999999999999803
RESULT: PASS 1
-----------------------------
(tcl) 16 % runSingleWhittakerEilersAutotest "zero_lambda_identity_order2" \
> {3 7 1 9 4 6 2 8} \
> 0.0 2 \
> {3.0 7.0 1.0 9.0 4.0 6.0 2.0 8.0} \
> $autotestRelativeTolerance $autotestAbsoluteTolerance
---- AUTOTEST: zero_lambda_identity_order2 ----
input: 3 7 1 9 4 6 2 8
lambda: 0.0
order: 2
expected: 3.0 7.0 1.0 9.0 4.0 6.0 2.0 8.0
actual: 3.0 7.0 1.0 9.0 4.0 6.0 2.0 8.0
RESULT: PASS 1
-----------------------------
(tcl) 17 % runSingleWhittakerEilersAutotest "single_peak_n5_lambda1_order2" \
> {0 0 10 0 0} \
> 1.0 2 \
> {0.41666667 2.5 4.16666667 2.5 0.41666667} \
> $autotestRelativeTolerance $autotestAbsoluteTolerance
---- AUTOTEST: single_peak_n5_lambda1_order2 ----
input: 0 0 10 0 0
lambda: 1.0
order: 2
expected: 0.41666667 2.5 4.16666667 2.5 0.41666667
actual: 0.4166666666666665 2.4999999999999996 4.166666666666666 2.4999999999999996 0.4166666666666666
RESULT: PASS 1
-----------------------------
(tcl) 18 % runSingleWhittakerEilersAutotest "single_peak_n3_lambda1_order2_minimum_size" \
> {0 1 0} \
> 1.0 2 \
> {0.28571429 0.42857143 0.28571429} \
> $autotestRelativeTolerance $autotestAbsoluteTolerance
---- AUTOTEST: single_peak_n3_lambda1_order2_minimum_size ----
input: 0 1 0
lambda: 1.0
order: 2
expected: 0.28571429 0.42857143 0.28571429
actual: 0.2857142857142857 0.42857142857142855 0.2857142857142857
RESULT: PASS 1
-----------------------------
(tcl) 19 %
This is a draft, still debugging on Playground V9. convert to strict 7-bit ASCII for Playground V9.
This page is under development. Comments are welcome, but please load any comments in the comments section at the bottom of the page. Please include your wiki MONIKER and date in your comment with the same courtesy that I will give you. Aside from your courtesy, your wiki MONIKER and date as a signature and minimal good faith of any internet post are the rules of this TCL-WIKI. Its very hard to reply reasonably without some background of the correspondent on his WIKI bio page. Thanks, gold 5Jan2026
gold 2/25/2026. Added categories, so can find message in Wiki.
gold 2/3/2025. Testing, encountered initial difficulty in saving work? Long code blocks with or unmatched wiki markup can sometimes confuse the Tcl Wiki formatting engine, especially if fences are not balanced or a line begins with markup it treats specially.
gold 2/14/2026. Added Automatic Dump of Examples, Using ActiveState.
gold 2/14/2026. convert to strict 7-bit ASCII for Playground V9. reporting error at bottom. program should run to completion with automatic test suite.
gold 2/14/2026. convert to strict 7-bit ASCII for Playground V9. variables need to be human readable and very explanatory. avoid variables with single letter names. Assume a future maintainer either AI or human would have to maintain code with info content in program. the program is working the numbers correctly . so minimal changes.
Please place any comments here with your wiki MONIKER and date, Thanks.gold 1/30/2026
Note. Testing computer methods and computer programs, maybe wrong numbers.
| Category Numerical Analysis | Category Toys | Category Calculator | Category Mathematics | Category Example | Toys and Games | Category Games | Category Application | Category GUI |