[Bug]: Simplifications cause heavisides to evaluate as booleans

[aabills]

PyBaMM Version

develop

Python Version

3.11

Describe the bug

Simplifications cause Heaviside functions to evaluate as booleans when not multiplied by something else. This makes it not possible to have a zoh as described in #3446 with a value of 1

Steps to Reproduce

import pybamm
import numpy as np
import matplotlib.pyplot as plt
pybamm.set_logging_level("DEBUG")

model = pybamm.lithium_ion.DFN()
param = model.default_parameter_values

# Works
param.update({"Current function [A]": pybamm.Scalar(1.000001) * (pybamm.t<10.0) + pybamm.Scalar(-1.000001) * (pybamm.t>=10.0)})
solver = pybamm.IDAKLUSolver()
sim = pybamm.Simulation(model, solver=solver, parameter_values=param)
sol = sim.solve(np.array([0.0, 10.0, 20.0]))

# Errors
param.update({"Current function [A]": pybamm.Scalar(1) * (pybamm.t<10.0) + pybamm.Scalar(-1) * (pybamm.t>=10.0)})
solver = pybamm.IDAKLUSolver()
sim = pybamm.Simulation(model, solver=solver, parameter_values=param)
sol = sim.solve(np.array([0.0, 10.0, 20.0]))

### Relevant log output

```shell
2025-06-25 14:13:28.785 - [DEBUG] discretisation.process_dict(738): Discretise 'Electrode current density [A.m-2]'
Traceback (most recent call last):
  File "/Users/abills/pybamm/pybamm/src/pybamm/discretisations/discretisation.py", line 769, in process_symbol
    return self._discretised_symbols[symbol]
           ~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^
KeyError: Concatenation(-0x2b7cf2d7ff82eb9f, concatenation, children=['-58.56620185738528 * grad(Negative electrode potential [V])', 'broadcast(0.0)', '-5.856620185738528 * grad(Positive electrode potential [V])'], domains={'primary': ['negative electrode', 'separator', 'positive electrode'], 'secondary': ['current collector']})

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/Users/abills/pybamm/pybamm/test.py", line 19, in <module>
    sol = sim.solve(np.array([0.0, 10.0, 20.0]))
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/abills/pybamm/pybamm/src/pybamm/simulation.py", line 474, in solve
    self.build(initial_soc=initial_soc, inputs=inputs)
  File "/Users/abills/pybamm/pybamm/src/pybamm/simulation.py", line 334, in build
    self._built_model = self._disc.process_model(
                        ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/abills/pybamm/pybamm/src/pybamm/discretisations/discretisation.py", line 213, in process_model
    model_disc.variables = self.process_dict(pre_processed_variables)
                           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/abills/pybamm/pybamm/src/pybamm/discretisations/discretisation.py", line 739, in process_dict
    processed_eqn = self.process_symbol(eqn)
                    ^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/abills/pybamm/pybamm/src/pybamm/discretisations/discretisation.py", line 771, in process_symbol
    discretised_symbol = self._process_symbol(symbol)
                         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/abills/pybamm/pybamm/src/pybamm/discretisations/discretisation.py", line 983, in _process_symbol
    new_symbol = spatial_method.concatenation(new_children)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/abills/pybamm/pybamm/src/pybamm/spatial_methods/finite_volume.py", line 1234, in concatenation
    child_size = child.size
                 ^^^^^^^^^^
  File "/Users/abills/miniconda3/envs/pybamm/lib/python3.11/functools.py", line 1001, in __get__
    val = self.func(instance)
          ^^^^^^^^^^^^^^^^^^^
  File "/Users/abills/pybamm/pybamm/src/pybamm/expression_tree/symbol.py", line 1019, in size
    return np.prod(self.shape)
                   ^^^^^^^^^^
  File "/Users/abills/miniconda3/envs/pybamm/lib/python3.11/functools.py", line 1001, in __get__
    val = self.func(instance)
          ^^^^^^^^^^^^^^^^^^^
  File "/Users/abills/pybamm/pybamm/src/pybamm/expression_tree/symbol.py", line 1030, in shape
    evaluated_self = self.evaluate(0, y, y, inputs="shape test")
                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/abills/pybamm/pybamm/src/pybamm/expression_tree/binary_operators.py", line 162, in evaluate
    right = self.right.evaluate(t, y, y_dot, inputs)
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/abills/pybamm/pybamm/src/pybamm/expression_tree/binary_operators.py", line 162, in evaluate
    right = self.right.evaluate(t, y, y_dot, inputs)
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/abills/pybamm/pybamm/src/pybamm/expression_tree/binary_operators.py", line 162, in evaluate
    right = self.right.evaluate(t, y, y_dot, inputs)
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/abills/pybamm/pybamm/src/pybamm/expression_tree/binary_operators.py", line 163, in evaluate
    return self._binary_evaluate(left, right)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/abills/pybamm/pybamm/src/pybamm/expression_tree/binary_operators.py", line 314, in _binary_evaluate
    return left - right
           ~~~~~^~~~~~~
TypeError: numpy boolean subtract, the `-` operator, is not supported, use the bitwise_xor, the `^` operator, or the logical_xor function instead.

[rtimms]

need this as part of fixing #4893 . hoping to wrap up today but incase you need a workaround sooner, you need this in binary_operators.py lines 1123-1129

    # anything multiplied by a scalar one returns itself
    if pybamm.is_scalar_one(left) and not isinstance(right, _Heaviside):
        return right

    # anything multiplied by a scalar negative one returns negative itself
    if pybamm.is_scalar_minus_one(left) and not isinstance(right, _Heaviside):
        return -right