Wrong reaction order in an autocatalytic reaction

[caprilesport]

Please complete the following tasks

• I have searched the documentation and help resources before reporting.
• I have searched the discussions to avoid duplicates.
• I have searched the open and closed issues to avoid duplicates.

Describe the bug
Consider the autocatalytic reaction "2B -> B + C":

import numpy as np

import overreact as rx

scheme = rx.parse_reactions("2 B -> B + C")

stoichiometry = np.asarray(scheme.A)
kinetic_orders = np.where(stoichiometry > 0, 0, -stoichiometry)
reaction_order = np.sum(kinetic_orders[:, 0])

print("reaction (should be 2nd order in B):", scheme.reactions[0])
print("computed net stoichiometry (A):", stoichiometry[:, 0])
print("reaction order:", reaction_order)

# actually computing the schema to show it's first order:
dydt = rx.get_dydt(scheme, k=np.array([1.0]))
rate_b1 = dydt(0.0, np.array([1.0, 0.0]))[1]
rate_b2 = dydt(0.0, np.array([2.0, 0.0]))[1]

print("rate at [B] = 1:", float(rate_b1))
print("rate at [B] = 2:", float(rate_b2))

Output:

reaction (should be 2nd order in B): 2 B -> B + C
computed net stoichiometry (A): [-1.  1.]
reaction order: 1.0
An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.
rate at [B] = 1: 1.0
rate at [B] = 2: 2.0

To Reproduce
Steps to reproduce the behavior:

1. Open a terminal.
2. Save the above in a script with a proper overreact environment.
3. Execute and see the inferred reaction order in B is 1 instead of 2.

Expected behavior
Should be 2nd order in B.

Origin

I think _add_reaction() when taken together with _dydt() is the faulty point.

 def _add_reaction(reactants, products, is_half_equilibrium, transition):
      """Local helper function with side-effects."""
      # TODO(schneiderfelipe): what if reaction is defined, then redefined
      # as equilibrium, or vice-versa?
      if (
          transition is not None
          and (reactants, products, is_half_equilibrium, transition) in reactions
      ):
          return

      # found new reaction
      reactions[(reactants, products, is_half_equilibrium, transition)] = None

      A_vector = np.zeros(len(compounds))
      for coefficient, reactant in reactants:
          A_vector[compounds[reactant]] = -coefficient
      B_vector = A_vector

      if transition is not None:
          B_vector = A_vector.copy()

          # it's assumed that
          #   1. there's a singe transition compound, and
          #   2. its coefficient equals one
          B_vector[compounds[transition[-1][-1]]] = 1

      for coefficient, product in products:
          A_vector[compounds[product]] += coefficient

      if (
          is_half_equilibrium
          and (products, reactants, is_half_equilibrium, transition) in reactions
      ):
          B_vector = np.zeros(len(compounds))

      A.append(A_vector)
      B.append(B_vector)

The parser sotres the reactants:

A_vector[B] = -2

Then it adds the products:

A_vector[B] += 1
A_vector[C] += 1

Which is fine for the net concentration changes, but when we get to _dydt(), more specifically, in the line below:

M = jnp.where(A > 0, 0, -A).T

we compute M based on A, which results in a first order reaction.

This has appeared when trying to implement the Robertson problem as a test in #817.