Skip to content

cds2.modeling

Symbolic mathematics: expression trees, differentiation, LaTeX export, equation solving.

Symbolic mathematics: expression trees, differentiation, LaTeX export.

Expression

Base class of the symbolic expression tree.

Source code in src/cds2/modeling.py 
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
class Expression:
    """Base class of the symbolic expression tree."""

    def __add__(self, other: object) -> Expression:
        return Add(self, _wrap(other))

    def __radd__(self, other: object) -> Expression:
        return Add(_wrap(other), self)

    def __sub__(self, other: object) -> Expression:
        return Subtract(self, _wrap(other))

    def __rsub__(self, other: object) -> Expression:
        return Subtract(_wrap(other), self)

    def __mul__(self, other: object) -> Expression:
        return Multiply(self, _wrap(other))

    def __rmul__(self, other: object) -> Expression:
        return Multiply(_wrap(other), self)

    def __truediv__(self, other: object) -> Expression:
        return Divide(self, _wrap(other))

    def __rtruediv__(self, other: object) -> Expression:
        return Divide(_wrap(other), self)

    def __pow__(self, exponent: object) -> Expression:
        return Power(self, _wrap(exponent))

    def __rpow__(self, base: object) -> Expression:
        return Power(_wrap(base), self)

    def __neg__(self) -> Expression:
        return Negate(self)

    def variables(self) -> set[str]:
        raise NotImplementedError

    def _latex(self) -> str:
        raise NotImplementedError

Constant

Bases: Expression

A fixed numeric literal.

Source code in src/cds2/modeling.py 
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
class Constant(Expression):
    """A fixed numeric literal."""

    def __init__(self, value: float) -> None:
        self.value = float(value)

    def variables(self) -> set[str]:
        return set()

    def _latex(self) -> str:
        text = f"{self.value:g}"
        return text

Variable

Bases: Expression

A named free variable.

Source code in src/cds2/modeling.py 
106
107
108
109
110
111
112
113
114
115
116
117
118
119
class Variable(Expression):
    """A named free variable."""

    def __init__(self, name: str) -> None:
        if not name.isidentifier():
            msg = f"variable name must be an identifier: {name!r}"
            raise ValueError(msg)
        self.name = name

    def variables(self) -> set[str]:
        return {self.name}

    def _latex(self) -> str:
        return self.name

BinaryOp

Bases: Expression

Base for binary operators.

Source code in src/cds2/modeling.py 
127
128
129
130
131
132
133
134
135
136
137
138
class BinaryOp(Expression):
    """Base for binary operators."""

    operator_symbol = "?"
    latex_operator = r"\operatorname{?}"

    def __init__(self, left: Expression, right: Expression) -> None:
        self.left = left
        self.right = right

    def variables(self) -> set[str]:
        return self.left.variables() | self.right.variables()

Negate

Bases: Expression

Unary negation.

Source code in src/cds2/modeling.py 
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
class Negate(Expression):
    """Unary negation."""

    def __init__(self, operand: Expression) -> None:
        self.operand = operand

    def variables(self) -> set[str]:
        return self.operand.variables()

    def _latex(self) -> str:
        inner = (
            f"({self.operand._latex()})"
            if isinstance(self.operand, (Add, Subtract))
            else self.operand._latex()
        )
        return f"-{inner}"

MathModel

Equation-centric modelling helper built on the expression tree.

Source code in src/cds2/modeling.py 
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
class MathModel:
    """Equation-centric modelling helper built on the expression tree."""

    def __init__(self, equation_lhs: Expression, equation_rhs: Expression) -> None:
        self.equation = simplify(Subtract(equation_lhs, equation_rhs))

    @classmethod
    def from_formula(cls, lhs: Expression, rhs: Expression) -> MathModel:
        return cls(lhs, rhs)

    def residual(self, bindings: dict[str, float]) -> float:
        """Evaluate LHS - RHS under complete variable bindings."""
        substituted = substitute(self.equation, dict(bindings))
        return evaluate(substituted, {})

    def solve_equation(
        self,
        unknown: str,
        target_value: float,
        known: dict[str, float],
        initial_guess: float = 1.0,
    ) -> float:
        """Solve for ``unknown`` so that LHS equals ``target_value`` (Newton)."""

        def root_function(value: float) -> float:
            bindings = dict(known)
            bindings[unknown] = float(value)
            try:
                substituted = substitute(self.equation, bindings)
                return evaluate(substituted, {})
            except ZeroDivisionError:
                return float("inf")

        guess = initial_guess
        for _ in range(60):
            function_value = root_function(guess) - target_value
            step = 1e-7 * max(1.0, abs(guess))
            derivative_value = (root_function(guess + step) - root_function(guess - step)) / (
                2 * step
            )
            if abs(derivative_value) < 1e-14:
                break
            next_guess = guess - function_value / derivative_value
            if abs(next_guess - guess) < 1e-12 * max(1.0, abs(guess)):
                guess = next_guess
                break
            guess = next_guess
        return float(guess)

    def fit_parameters(
        self,
        parameters: dict[str, float],
        observations: list[tuple[dict[str, float], float]],
        learning_rate: float = 1e-3,
        epochs: int = 500,
    ) -> dict[str, float]:
        """Grid-free least-squares refinement of free parameters by GD."""
        current = dict(parameters)
        for _ in range(epochs):
            gradient: dict[str, float] = {name: 0.0 for name in current}
            for known, target in observations:
                predicted_bindings = {**known, **current}
                substituted = substitute(self.equation, predicted_bindings)
                error = evaluate(substituted, {}) - target
                for name in current:
                    probe = dict(predicted_bindings)
                    probe[name] += 1e-6
                    bumped = evaluate(substitute(self.equation, probe), {})
                    gradient[name] += 2.0 * error * (bumped - evaluate(substituted, {})) / 1e-6
            for name in current:
                current[name] -= learning_rate * gradient[name] / len(observations)
        return current

residual(bindings)

Evaluate LHS - RHS under complete variable bindings.

Source code in src/cds2/modeling.py 
557
558
559
560
def residual(self, bindings: dict[str, float]) -> float:
    """Evaluate LHS - RHS under complete variable bindings."""
    substituted = substitute(self.equation, dict(bindings))
    return evaluate(substituted, {})

solve_equation(unknown, target_value, known, initial_guess=1.0)

Solve for unknown so that LHS equals target_value (Newton).

Source code in src/cds2/modeling.py 
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
def solve_equation(
    self,
    unknown: str,
    target_value: float,
    known: dict[str, float],
    initial_guess: float = 1.0,
) -> float:
    """Solve for ``unknown`` so that LHS equals ``target_value`` (Newton)."""

    def root_function(value: float) -> float:
        bindings = dict(known)
        bindings[unknown] = float(value)
        try:
            substituted = substitute(self.equation, bindings)
            return evaluate(substituted, {})
        except ZeroDivisionError:
            return float("inf")

    guess = initial_guess
    for _ in range(60):
        function_value = root_function(guess) - target_value
        step = 1e-7 * max(1.0, abs(guess))
        derivative_value = (root_function(guess + step) - root_function(guess - step)) / (
            2 * step
        )
        if abs(derivative_value) < 1e-14:
            break
        next_guess = guess - function_value / derivative_value
        if abs(next_guess - guess) < 1e-12 * max(1.0, abs(guess)):
            guess = next_guess
            break
        guess = next_guess
    return float(guess)

fit_parameters(parameters, observations, learning_rate=0.001, epochs=500)

Grid-free least-squares refinement of free parameters by GD.

Source code in src/cds2/modeling.py 
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
def fit_parameters(
    self,
    parameters: dict[str, float],
    observations: list[tuple[dict[str, float], float]],
    learning_rate: float = 1e-3,
    epochs: int = 500,
) -> dict[str, float]:
    """Grid-free least-squares refinement of free parameters by GD."""
    current = dict(parameters)
    for _ in range(epochs):
        gradient: dict[str, float] = {name: 0.0 for name in current}
        for known, target in observations:
            predicted_bindings = {**known, **current}
            substituted = substitute(self.equation, predicted_bindings)
            error = evaluate(substituted, {}) - target
            for name in current:
                probe = dict(predicted_bindings)
                probe[name] += 1e-6
                bumped = evaluate(substitute(self.equation, probe), {})
                gradient[name] += 2.0 * error * (bumped - evaluate(substituted, {})) / 1e-6
        for name in current:
            current[name] -= learning_rate * gradient[name] / len(observations)
    return current

symbol(name)

Shorthand factory for a named :class:Variable.

Source code in src/cds2/modeling.py 
122
123
124
def symbol(name: str) -> Variable:
    """Shorthand factory for a named :class:`Variable`."""
    return Variable(name)

diff(expression, variable)

Symbolic partial derivative of expression w.r.t. variable.

Source code in src/cds2/modeling.py 
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
def diff(expression: Expression, variable: str) -> Expression:
    """Symbolic partial derivative of ``expression`` w.r.t. ``variable``."""
    if isinstance(expression, Constant):
        return Constant(0.0)
    if isinstance(expression, Variable):
        return Constant(1.0 if expression.name == variable else 0.0)
    if isinstance(expression, Negate):
        return Negate(diff(expression.operand, variable))
    if isinstance(expression, Add):
        return Add(diff(expression.left, variable), diff(expression.right, variable))
    if isinstance(expression, Subtract):
        return Subtract(diff(expression.left, variable), diff(expression.right, variable))
    if isinstance(expression, Multiply):
        return Add(
            Multiply(diff(expression.left, variable), expression.right),
            Multiply(expression.left, diff(expression.right, variable)),
        )
    if isinstance(expression, Divide):
        numerator = Subtract(
            Multiply(diff(expression.left, variable), expression.right),
            Multiply(expression.left, diff(expression.right, variable)),
        )
        return Divide(numerator, Power(expression.right, Constant(2.0)))
    if isinstance(expression, Power):
        if isinstance(expression.right, Constant):
            reduced_power = Power(expression.left, Constant(expression.right.value - 1.0))
            return Multiply(
                Multiply(expression.right, reduced_power),
                diff(expression.left, variable),
            )
        if isinstance(expression.left, Constant):
            log_term = Multiply(Constant(math.log(expression.left.value)), expression)
            return Multiply(log_term, diff(expression.right, variable))
        raise NotImplementedError("d/dx u(x)^v(x): take logarithms or provide a numeric base")
    raise NotImplementedError(f"differentiation not defined for {type(expression).__name__}")

integrate(expression, variable)

Symbolic indefinite integral of a polynomial expression w.r.t. variable.

Handles sums, constant factors, powers with numeric exponents (except x**-1) and plain variables; anything else raises NotImplementedError. The integration constant is omitted.

Source code in src/cds2/modeling.py 
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
def integrate(expression: Expression, variable: str) -> Expression:
    """Symbolic indefinite integral of a polynomial expression w.r.t. ``variable``.

    Handles sums, constant factors, powers with numeric exponents (except
    ``x**-1``) and plain variables; anything else raises ``NotImplementedError``.
    The integration constant is omitted.
    """
    if isinstance(expression, Constant):
        return Multiply(expression, Variable(variable))
    if isinstance(expression, Variable):
        if expression.name == variable:
            return Divide(Power(Variable(variable), Constant(2.0)), Constant(2.0))
        return Multiply(expression, Variable(variable))
    if isinstance(expression, Negate):
        return Negate(integrate(expression.operand, variable))
    if isinstance(expression, Add):
        return Add(integrate(expression.left, variable), integrate(expression.right, variable))
    if isinstance(expression, Subtract):
        return Subtract(integrate(expression.left, variable), integrate(expression.right, variable))
    if isinstance(expression, Multiply):
        left_free = variable in expression.left.variables()
        right_free = variable in expression.right.variables()
        if not left_free and not right_free:
            return Multiply(expression, Variable(variable))
        if not left_free:
            return Multiply(expression.left, integrate(expression.right, variable))
        if not right_free:
            return Multiply(integrate(expression.left, variable), expression.right)
        raise NotImplementedError("integration of x*x style products needs expansion first")
    if isinstance(expression, Power):
        base = expression.left
        exponent = expression.right
        if isinstance(base, Variable) and base.name == variable and isinstance(exponent, Constant):
            if math.isclose(exponent.value, -1.0):
                raise NotImplementedError("integral of 1/x is out of scope")
            new_power = Power(base, Constant(exponent.value + 1.0))
            return Divide(new_power, Constant(exponent.value + 1.0))
        if variable not in expression.variables():
            return Multiply(expression, Variable(variable))
        raise NotImplementedError("unsupported power pattern in integral")
    if isinstance(expression, Divide):
        numerator_free = variable in expression.left.variables()
        denominator_free = variable in expression.right.variables()
        if not denominator_free:
            return Divide(integrate(expression.left, variable), expression.right)
        if not numerator_free and isinstance(expression.right, Variable):
            raise NotImplementedError("integral of 1/x is out of scope")
        raise NotImplementedError("unsupported division pattern in integral")
    raise NotImplementedError(f"integration not defined for {type(expression).__name__}")

polynomial_coefficients(expression, variable)

Coefficients (ascending powers) of the polynomial in variable.

Raises ValueError when the expression is not a polynomial in the requested variable.

Source code in src/cds2/modeling.py 
299
300
301
302
303
304
305
def polynomial_coefficients(expression: Expression, variable: str) -> FloatArray:
    """Coefficients (ascending powers) of the polynomial in ``variable``.

    Raises ``ValueError`` when the expression is not a polynomial in the
    requested variable.
    """
    return np.asarray(_poly_coefficients(expression, variable), dtype=float)

solve_polynomial(lhs, rhs, unknown, real_only=False)

Roots of lhs = rhs treated as a polynomial in unknown.

Returns complex roots (ascending powers input); with real_only only the real roots are returned, sorted ascending.

Source code in src/cds2/modeling.py 
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
def solve_polynomial(
    lhs: Expression,
    rhs: Expression,
    unknown: str,
    real_only: bool = False,
) -> NDArray[np.float64] | NDArray[np.complex128]:
    """Roots of ``lhs = rhs`` treated as a polynomial in ``unknown``.

    Returns complex roots (ascending powers input); with ``real_only`` only
    the real roots are returned, sorted ascending.
    """
    combined = Subtract(lhs, rhs)
    coefficients = np.asarray(_poly_coefficients(combined, unknown), dtype=float)
    degree = coefficients.size - 1
    while degree > 0 and math.isclose(coefficients[degree], 0.0):
        degree -= 1
    if degree >= 1:
        # numpy expects descending powers; our tree walk produces ascending.
        return_array: NDArray[np.complex128] = np.roots(coefficients[: degree + 1][::-1]).astype(
            np.complex128
        )
    else:
        return_array = np.zeros(0, dtype=np.complex128)
    if real_only:
        reals = return_array[np.abs(return_array.imag) < 1e-9].real
        return np.sort(np.asarray(reals, dtype=np.float64))
    return return_array

simplify(expression, max_rounds=8)

Apply algebraic simplification rules until stable.

Source code in src/cds2/modeling.py 
394
395
396
397
398
399
400
401
402
def simplify(expression: Expression, max_rounds: int = 8) -> Expression:
    """Apply algebraic simplification rules until stable."""
    current = expression
    for _ in range(max_rounds):
        simplified = _simplify_once(current)
        if _signature(simplified) == _signature(current):
            return simplified
        current = simplified
    return current

to_string(expression)

Plain-text rendering used for structural comparison and printing.

Source code in src/cds2/modeling.py 
471
472
473
474
475
476
477
478
479
480
481
482
def to_string(expression: Expression) -> str:
    """Plain-text rendering used for structural comparison and printing."""
    if isinstance(expression, Constant):
        return f"#({expression.value:g})"
    if isinstance(expression, Variable):
        return expression.name
    if isinstance(expression, Negate):
        return f"(~{to_string(expression.operand)})"
    if isinstance(expression, BinaryOp):
        inner = f"{to_string(expression.left)} {expression.operator_symbol} {to_string(expression.right)}"
        return f"({inner})"
    raise NotImplementedError(type(expression).__name__)

substitute(expression, replacements)

Replace variables with values or other expressions.

Source code in src/cds2/modeling.py 
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
def substitute(expression: Expression, replacements: Mapping[str, object]) -> Expression:
    """Replace variables with values or other expressions."""

    def walk(node: Expression) -> Expression:
        if isinstance(node, Constant):
            return node
        if isinstance(node, Variable):
            if node.name in replacements:
                return _wrap(replacements[node.name])
            return node
        if isinstance(node, Negate):
            return Negate(walk(node.operand))
        if isinstance(node, BinaryOp):
            return _rebuild(type(node), walk(node.left), walk(node.right))
        raise NotImplementedError(type(node).__name__)

    return simplify(walk(expression))

evaluate(expression, environment)

Numerically evaluate the expression under a variable binding.

Source code in src/cds2/modeling.py 
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
def evaluate(expression: Expression, environment: dict[str, float]) -> float:
    """Numerically evaluate the expression under a variable binding."""

    def walk(node: Expression) -> float:
        if isinstance(node, Constant):
            return node.value
        if isinstance(node, Variable):
            if node.name not in environment:
                msg = f"unbound variable: {node.name}"
                raise KeyError(msg)
            return float(environment[node.name])
        if isinstance(node, Negate):
            return -walk(node.operand)
        if isinstance(node, Add):
            return float(walk(node.left) + walk(node.right))
        if isinstance(node, Subtract):
            return float(walk(node.left) - walk(node.right))
        if isinstance(node, Multiply):
            return float(walk(node.left) * walk(node.right))
        if isinstance(node, Divide):
            return float(walk(node.left) / walk(node.right))
        if isinstance(node, Power):
            return float(walk(node.left) ** walk(node.right))
        raise NotImplementedError(type(node).__name__)

    return walk(expression)

iter_postorder(expression)

Yield expression nodes depth-first, children before parents.

Source code in src/cds2/modeling.py 
532
533
534
535
536
537
538
539
def iter_postorder(expression: Expression) -> Iterator[Expression]:
    """Yield expression nodes depth-first, children before parents."""
    if isinstance(expression, Negate):
        yield from iter_postorder(expression.operand)
    elif isinstance(expression, BinaryOp):
        yield from iter_postorder(expression.left)
        yield from iter_postorder(expression.right)
    yield expression

to_latex(expression)

Render the expression as a LaTeX fragment.

Source code in src/cds2/modeling.py 
542
543
544
def to_latex(expression: Expression) -> str:
    """Render the expression as a LaTeX fragment."""
    return expression._latex()