sys-bio / roadrunner

libRoadRunner: A high-performance SBML simulator
http://libroadrunner.org/
Other
39 stars 24 forks source link

Calculating steady state changes the selection vector #1148

Closed luciansmith closed 10 months ago

luciansmith commented 11 months ago

The selection vector needs to be left alone and not changed. Here's an example of the problem:

import tellurium as te
import roadrunner

r = te.loada("""
// Created by libAntimony v2.13.4
model *DB_for_loglinear()

  // Compartments and Species:
  species x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14;
  species x15;

  // Reactions:
  _J0: x2 + x4 -> x0 + x6; kf0*x2*x4;
  _J1: x0 + x6 -> x2 + x4; kb0*x0*x6;
  _J2: x2 + x7 -> x3 + x6; kf1*x2*x7;
  _J3: x3 + x6 -> x2 + x7; kb1*x3*x6;
  _J4: x3 + x5 -> x1 + x7; kf2*x3*x5;
  _J5: x1 + x7 -> x3 + x5; kb2*x1*x7;
  _J6: x0 + x10 -> x2 + x8; kf3*x0*x10;
  _J7: x2 + x8 -> x0 + x10; kb3*x2*x8;
  _J8: x1 + x11 -> x3 + x9; kf4*x1*x11;
  _J9: x3 + x9 -> x1 + x11; kb4*x3*x9;
  _J10: x3 + x10 -> x2 + x11; kf5*x3*x10;
  _J11: x2 + x11 -> x3 + x10; kb5*x2*x11;
  _J12: x4 + x14 -> x6 + x12; kf6*x4*x14;
  _J13: x6 + x12 -> x4 + x14; kb6*x6*x12;
  _J14: x5 + x15 -> x7 + x13; kf7*x5*x15;
  _J15: x7 + x13 -> x5 + x15; kb7*x7*x13;
  _J16: x7 + x14 -> x6 + x15; kf8*x7*x14;
  _J17: x6 + x15 -> x7 + x14; kb8*x6*x15;

  // Species initializations:
  x0 = 5;
  x1 = 10;
  x2 = 5;
  x3 = 11;
  x4 = 8;
  x5 = 9;
  x6 = 8;
  x7 = 10;
  x8 = 3;
  x9 = 6;
  x10 = 4;
  x11 = 7;
  x12 = 5;
  x13 = 3;
  x14 = 10;
  x15 = 12;

  // Variable initializations:
  kf0 = 1;
  kb0 = 1;
  kf1 = 1;
  kb1 = 1;
  kf2 = 1;
  kb2 = 1;
  kf3 = 1;
  kb3 = 1;
  kf4 = 1;
  kb4 = 1;
  kf5 = 1;
  kb5 = 1;
  kf6 = 1;
  kb6 = 1;
  kf7 = 1;
  kb7 = 1;
  kf8 = 1;
  kb8 = 1;

  // Other declarations:
  const kf0, kb0, kf1, kb1, kf2, kb2, kf3, kb3, kf4, kb4, kf5, kb5, kf6, kb6;
  const kf7, kb7, kf8, kb8;
end

""")

print(r.steadyStateSelections)
# r.steadyStateSelections = ['x0']
print(r.getSteadyStateValues())
print(r.steadyStateSelections)
# r.steadyStateSelections = ['x0']
print(r.getSteadyStateValues()) 

The SteadyStateSelections vector is completely re-ordered when 'getSteadyStateValues' is called, and the resulting values correspond to the new order, not the old. This is particularly problematic when the steady state selection vector was set by the user, as in the commented-out bits: it was set to x0, but gives the result for x2!