Scientific programming and reproducible computation

State the model before the program

Scientific programming expresses a model or analysis as operations a computer can perform. Start with the question and quantities. Our original teaching model is a rectangle with side lengths L and W and area A = L W. L and W are lengths in metres, and A is an area in square metres. The values are exact stipulated inputs to a mathematical model, not measurements of an assessed physical object.

Choose the local domain 0 ≤ L ≤ 10 m and 0 ≤ W ≤ 10 m. Zero is allowed as a degenerate mathematical boundary with zero area. The upper limit is a teaching-program restriction, not a universal physical limit on rectangles. The program must reject inputs outside this contract rather than silently assign them an unintended meaning.

Separate the physical quantities from their numerical representations. A computer value such as 2.5 carries no metre unit by itself in the program below. Our interface requires the numerical metre values and returns the numerical square-metre value. Names and documentation communicate that contract, but cannot establish that a caller actually converted the input correctly.

Convert input units explicitly

Supply L = 250 cm and W = 125 cm. Since 1 m = 100 cm, their numerical metre values are 250/100 = 2.5 and 125/100 = 1.25. The quantity equation gives A = (2.5 m)(1.25 m) = 3.125 m². Multiplying the original centimetre values instead gives 31250 cm², which is the same area because 1 m² = 10000 cm².

Passing the raw numbers 250 and 125 to an interface that expects numerical metre values would represent different lengths, outside our local domain. The guarded function below rejects them. In an unguarded multiplication, attaching m² to their product 31250 would overstate this stipulated area by a factor of 10000. The bug would be in the quantity interpretation, even if the multiplication itself ran normally.

Keep the input unit, conversion and output unit alongside the run. Convert each input once at a clearly identified stage. A check that a value is a number or lies between 0 and 10 cannot detect every unit mistake: for example, the bare value 2 could mean a numerical metre value or an unconverted centimetre value. Correct units require a reliable input record as well as code checks.

Define a small explicit function

Use these original Python 3 definitions in order: def reject_length(): raise ValueError("length-domain"); def valid_length(value): return (type(value) is int or type(value) is float) and 0 <= value <= 10; def area_m2(length_m, width_m): return length_m * width_m if valid_length(length_m) and valid_length(width_m) else reject_length(). Each definition is a complete single line. Python permits a simple function body on the same line after the colon. def defines a function, the parenthesized names are its inputs, and return supplies its result.

The predicate valid_length returns a truth value indicating whether its input meets our numerical contract. (type(value) is int or type(value) is float) deliberately accepts only the exact built-in integer and floating-point types; it excludes truth values such as True, text, complex numbers and subclasses. Here is compares the identity of type objects, and or permits either listed type. This is a local teaching choice, not a rule that scientific programs should reject every other numeric type. The range check also rejects negative values, values above 10, infinities and a floating-point not-a-number value. Python's and tests the type condition first and does not evaluate the range condition when that type test fails.

In area_m2, the conditional expression checks the input contract before evaluating its selected branch. Multiplication occurs only if both checks pass; if the first check fails, and skips the second. A rejected input calls reject_length, which raises ValueError("length-domain") rather than returning an area. The variables encode numerical metre values, and the returned number is interpreted in m². These built-in arithmetic operations use the interpreter's numeric rules: many decimal floating-point inputs are approximate. The selected main-case values 2.5, 1.25 and 3.125 have exact finite binary representations; that fact does not make every accepted floating-point product exact.

def reject_length(): raise ValueError("length-domain")
def valid_length(value): return (type(value) is int or type(value) is float) and 0 <= value <= 10
def area_m2(length_m, width_m): return length_m * width_m if valid_length(length_m) and valid_length(width_m) else reject_length()

Derive tests before trusting output

Derive an expected answer independently of the code being tested. For the supplied main case, ordinary quantity arithmetic gives 3.125 m², so area_m2(2.5, 1.25) should return the numerical value 3.125. Also check the swapped inputs, the allowed zero boundary and the upper boundary. These selected values are exactly represented in the checked interpreter; an equality check is appropriate for these particular expected results.

Different tests address different faults. Swapping the two inputs checks symmetry for that case; zero checks a boundary. A deliberately wrong rule A_bad = 2 L W still passes both of those checks while returning 6.25 m² for the main case. Therefore symmetry and a zero result alone do not establish the rectangle formula. A known-answer test distinguishes this wrong implementation.

Test rejected inputs as well as accepted ones: area_m2(-1, 1), area_m2(10.5, 1) and area_m2(True, 1) should raise the stated exception. Selected tests provide evidence about the exercised cases, not a proof of all program behavior or physical validity. Reference recommends tests against simplified expectations and preserving a discovered bug as a regression test, meaning a check rerun after changes to detect its return.

Make the execution order visible

For a complete main-case run, put the three definitions above into a file named area_run.py, each on its own line, then append these three lines in order: length_m = 250 / 100, width_m = 125 / 100 and print(area_m2(length_m, width_m)). Run the saved file using python area_run.py with the intended Python 3 interpreter. The printed number is 3.125; the interface supplies its m² interpretation. No external package, random draw or measured dataset is needed for this supplied run.

A function is one part of an execution recipe. The conversions, definitions and call must belong to the intended saved version. An interactive session can retain names from earlier work; executing a later cell may then use an older value or function. A fresh run of the complete saved recipe checks that the displayed result does not depend on such leftover state. It does not prove that every interactive workflow is irreproducible.

Record the interpreter version and relevant environment separately from the source file. An environment can include software versions, operating system, hardware and configuration; The Turing Way explains why these conditions matter for computational reruns. The article's displayed syntax was checked using Python 3.12.14, while the directly inspected language documentation reports Python 3.13.16. These are distinct version records, not a claim that all versions or machines produce identical behavior.

def reject_length(): raise ValueError("length-domain")
def valid_length(value): return (type(value) is int or type(value) is float) and 0 <= value <= 10
def area_m2(length_m, width_m): return length_m * width_m if valid_length(length_m) and valid_length(width_m) else reject_length()
length_m = 250 / 100
width_m = 125 / 100
print(area_m2(length_m, width_m))

Declare a numerical comparison rule

First say which outputs are being compared. Exact numerical equality, identical output text and identical stored bits are different criteria. A joint Japanese–French numerical-computing workshop discusses a distinction between strong bit-level and weaker numerical reproducibility in its own context. Here we define a separate deterministic comparison rule for finite area values, without adopting that source's probabilistic method or assuming agreement across every platform.

Let r be a reference area and z a candidate area, both in m². Choose an absolute threshold α = 0.02 m² and a dimensionless relative coefficient ρ = 0.01. Our local criterion is |z − r| ≤ max(α, ρ max(|z|, |r|)). The absolute difference and both quantities inside the outer maximum have units m². The absolute term gives a usable comparison near zero. These thresholds are stipulated choices for this exercise, not measurement uncertainties, confidence intervals or universal accuracy requirements.

For r = 3.125 m² and the separately supplied candidate z = 3.15 m², the exact stipulated decimal values give a difference of 0.025 m² and threshold 0.0315 m², so the criterion passes. For z = 3.20 m², the difference is 0.075 m² and threshold 0.032 m², so it fails. These candidates illustrate a checker; they are not claimed outputs of the main-case function. A passing comparison means agreement under this declared rule, not exact equality or proof that either result is physically correct.

Preserve a reviewable run recipe

Keep the exact program revision, input numbers with their units, conversions, execution command and order, interpreter and relevant dependency versions, configuration, output interpretation and test results. For this small run, explicitly record that it uses Python built-ins and no external packages or random numbers. Larger analyses may require additional inputs and environment details; their requirements depend on the computation. References support recording computational provenance and environment conditions.

Version control records identifiable source revisions and changes so that a prior program can be retrieved and compared. A saved revision helps answer which code produced a result. It does not by itself capture every untracked input or runtime setting, ensure that the recorded command was used, or establish code correctness. Pair the revision with the run recipe and retain the relevant outputs and checks.

When a change fixes the doubled-area bug, keep the independently expected main-case answer and rerun it with the boundary and rejection checks. Explain what changed and why the old output was wrong. If a dependency, input conversion or comparison criterion changes, record that change too. Another reader should be able to distinguish a faithful rerun of the old procedure from a revised scientific calculation.

Separate rerunning from validity

Three questions remain separate: can the stated procedure be rerun, does its code implement the intended model, and is that model adequate for the intended physical situation? A repeated result addresses the first question under the recorded conditions. Independent answers and code review address implementation. Physical evidence and the model's assumptions address adequacy. Repeating the doubled-area rule can reproduce the same wrong 6.25 m² indefinitely.

For the supplied rectangle, the mathematical formula and unit conversions are established locally. No physical measurement, uncertainty estimate, nonrectangular boundary or material behavior has been assessed. The small input guard catches specified contract violations, but cannot discover every mislabeled quantity or justify using a rectangle model for another shape. Likewise, the numerical comparison thresholds encode an acceptance decision rather than a scientific truth test.

The next map subject, Floating-point arithmetic and numerical stability, examines how numerical representation and algorithms affect computations. This lesson has checked selected exact cases and bounded program behavior, without promising universal floating-point equality. A clear model, explicit interface, independent checks and a recorded run recipe make a computation open to inspection; the relevant scientific conclusions still require their own evidence.

Worked examples

Try it yourself

1. For a new stipulated rectangle with L = 75 cm and W = 50 cm, give the numerical metre inputs and expected output of the displayed function, including the output's physical unit.

   Show the answer and reasoning

   Use 0.75 and 0.5; the expected numerical output is 0.375, interpreted as 0.375 m².

   Dividing both centimetre values by 100 gives the metre values. Their product is (0.75 m)(0.5 m) = 0.375 m². Both lie in the allowed interval, and these selected fractions have exact finite binary representations in the checked interpreter.

2. Under the exact displayed definitions, classify area_m2(0, 10), area_m2(10.5, 1) and area_m2(True, 1). Does rejecting a truth value mean every scientific program must reject all other numeric types?

   Show the answer and reasoning

   The first returns 0. The second and third raise ValueError("length-domain"). The exact-type restriction is this lesson's local contract, not a universal scientific-programming requirement.

   Zero and 10 are allowed boundaries. The value 10.5 exceeds the numerical range, and True has the built-in truth-value type rather than the exact int or float type required here. Different interfaces can deliberately accept other types if their behavior and quantities are defined.

3. The wrong rule A_bad = 2 L W returns 0 m² for a zero-length input and agrees with itself when its inputs are swapped. For L = 0.75 m and W = 0.5 m, give the correct and wrong outputs. Have the two passing property checks proved the rule correct?

   Show the answer and reasoning

   The correct area is 0.375 m² and the doubled rule gives 0.75 m². The passing zero and symmetry checks have not proved correctness.

   The original model independently gives A = L W. Multiplying it by 2 preserves zero and symmetry while changing the nonzero answer. A known-answer case detects this specified fault; the same checks do not exhaust every possible fault.

4. Apply the same finite-area comparison rule with reference r = 0 m², α = 0.02 m² and ρ = 0.01. Treat z = 0.015 m² and z = 0.03 m² as exact supplied candidates. Which passes?

   Show the answer and reasoning

   The 0.015 m² candidate passes; the 0.03 m² candidate fails.

   For the first, the difference is 0.015 m² and the relative term is 0.00015 m², so the threshold is 0.02 m². For the second, the difference is 0.03 m² and the relative term is 0.0003 m², leaving the same threshold 0.02 m². The absolute term makes the declared criterion usable at a zero reference; it does not set a measurement uncertainty.

5. A record stores only the number 3.125 and the name area_run.py. Name four kinds of missing information needed to inspect this particular computation. Does keeping a filename alone identify the exact program used?

   Show the answer and reasoning

   Missing information includes the exact source revision, input values with units and conversions, interpreter/environment details, and the execution command/order with output units and check results. A filename alone does not identify the exact source version.

   The same filename can contain different definitions or conversions. A reviewable recipe links the numerical output to the saved source, stipulated input interpretation, execution conditions and tests. Which additional environment or input details matter depends on the program.

6. Two fresh runs of the doubled-area program both return 6.25 for the main inputs. A file comparison also shows identical output text. What has been established, and what remains wrong or unassessed?

   Show the answer and reasoning

   The two recorded runs agree in their output text under those conditions. The doubled rule still gives the wrong area for the stated rectangle model, and no actual physical object or measurement uncertainty has been assessed.

   Repetition and identical text do not turn A_bad = 2 L W into A = L W. The independently derived main-case area is 3.125 m². Code implementation and physical adequacy need their own checks and evidence in addition to a rerun record.

What to carry forward

• State the physical quantities and model before defining numerical inputs and outputs.
• A bare number does not verify its unit; keep conversions and quantity interpretation explicit.
• A small function's domain and rejected inputs are part of its interface.
• Independent known answers and different selected checks reveal different faults; passing tests alone is not universal proof.
• Declare the equality or numerical-agreement criterion, including units and near-zero behavior.
• Preserve the exact source and run conditions, and distinguish reproducible output from correct implementation and physical adequacy.