`generate_steady_states` fails when you edit a model's optimal percentage twice

[hariszaf]

model = dingo.MetabolicNetwork.from_json("ext_data/e_coli_core.json")

model.set_opt_percentage(90)
sampler = dingo.PolytopeSampler(model)
sampler.generate_steady_states()

model.set_opt_percentage(20)
sampler = dingo.PolytopeSampler(model)
sampler.generate_steady_states()

would return:

phase 1: number of correlated samples = 500, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 4548.22
phase 2: number of correlated samples = 500, effective sample size = 123, ratio of the maximum singilar value over the minimum singular value = 3.07681
phase 3: number of correlated samples = 500, effective sample size = 148, ratio of the maximum singilar value over the minimum singular value = 2.8596
phase 4: number of correlated samples = 1900, effective sample size = 754
[5]total ess 1028: number of correlated samples = 3400


[5]maximum marginal PSRF: 1.03027


UnboundLocalError                         Traceback (most recent call last)
[/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/dev.ipynb](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/dev.ipynb) Cell 26 line 8
      [6](vscode-notebook-cell:/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/dev.ipynb#X36sZmlsZQ%3D%3D?line=5) model.set_opt_percentage(20)
      [7](vscode-notebook-cell:/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/dev.ipynb#X36sZmlsZQ%3D%3D?line=6) sampler = dingo.PolytopeSampler(model)
----> [8](vscode-notebook-cell:/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/dev.ipynb#X36sZmlsZQ%3D%3D?line=7) sampler.generate_steady_states()

File [~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:161](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:161), in PolytopeSampler.generate_steady_states(self, ess, psrf, parallel_mmcs, num_threads)
    [149](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:149) def generate_steady_states(
    [150](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:150)     self, ess=1000, psrf=False, parallel_mmcs=False, num_threads=1
    [151](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:151) ):
    [152](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:152)     """A member function to sample steady states.
    [153](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:153) 
    [154](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:154)     Keyword arguments:
   (...)
    [158](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:158)     num_threads -- the number of threads to use for parallel mmcs
    [159](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:159)     """
--> [161](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:161)     self.get_polytope()
    [163](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:163)     P = HPolytope(self._A, self._b)
    [165](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:165)     if self._parameters["fast_computations"]:

File [~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:82](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:82), in PolytopeSampler.get_polytope(self)
     [66](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:66) """A member function to derive the corresponding full dimensional polytope
     [67](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:67) and a isometric linear transformation that maps the latter to the initial space.
     [68](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:68) """
     [70](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:70) if (
     [71](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:71)     self._A == []
     [72](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:72)     or self._b == []
   (...)
     [76](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:76)     or self._T_shift == []
     [77](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:77) ):
     [79](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:79)     (
     [80](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:80)         max_biomass_flux_vector,
     [81](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:81)         max_biomass_objective,
---> [82](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:82)     ) = self._metabolic_network.fba()
     [84](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:84)     if (
     [85](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:85)         self._parameters["fast_computations"]
     [86](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:86)         and self._parameters["remove_redundant_facets"]
     [87](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:87)     ):
     [89](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:89)         A, b, Aeq, beq = fast_remove_redundant_facets(
     [90](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:90)             self._metabolic_network.lb,
     [91](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:91)             self._metabolic_network.ub,
   (...)
     [94](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:94)             self._parameters["opt_percentage"],
     [95](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/PolytopeSampler.py:95)         )

File [~/github_repos/GeomScale/dingo/dingo/MetabolicNetwork.py:117](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/MetabolicNetwork.py:117), in MetabolicNetwork.fba(self)
    [114](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/MetabolicNetwork.py:114) """A member function to apply the FBA method on the metabolic network."""
    [116](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/MetabolicNetwork.py:116) if self._parameters["fast_computations"]:
--> [117](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/MetabolicNetwork.py:117)     return fast_fba(self._lb, self._ub, self._S, self._biomass_function)
    [118](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/MetabolicNetwork.py:118) else:
    [119](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/MetabolicNetwork.py:119)     return slow_fba(self._lb, self._ub, self._S, self._biomass_function)

File [~/github_repos/GeomScale/dingo/dingo/gurobi_based_implementations.py:112](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/gurobi_based_implementations.py:112), in fast_fba(lb, ub, S, c)
    [109](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/gurobi_based_implementations.py:109)     v = model.getVars()
    [111](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/gurobi_based_implementations.py:111) for i in range(n):
--> [112](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/gurobi_based_implementations.py:112)     optimum_sol.append(v[i].x)
    [114](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/gurobi_based_implementations.py:114) optimum_sol = np.asarray(optimum_sol)
    [116](https://file+.vscode-resource.vscode-cdn.net/home/luna.kuleuven.be/u0156635/github_repos/GeomScale/dingo/~/github_repos/GeomScale/dingo/dingo/gurobi_based_implementations.py:116) return optimum_sol, optimum_value

UnboundLocalError: local variable 'v' referenced before assignment

It seems the error is related to the https://github.com/GeomScale/dingo/blob/aaae4ae63f432da45eb0f4363a92767ba537a074/dingo/gurobi_based_implementations.py#L109
and the status == GRB.OPTIMAL that is not the case in the second edit.
where the status is INF_OR_UNBD for some reason (https://www.gurobi.com/documentation/current/refman/optimization_status_codes.html)

[AdityaTaggar05]

Hi,
I tried to understand what could be causing the problem, but the thing is, it's confusing.
First of all, I didn't set up gurobi on my system, and maybe that is why I am not getting the same error.
Secondly, when changing the optimal percentage value, the error varies.

Using optimal percentage = 10/20 %:

model = dingo.MetabolicNetwork.from_json("ext_data/e_coli_core.json")

model.set_opt_percentage(10)
sampler = dingo.PolytopeSampler(model)
sampler.generate_steady_states()

model.set_opt_percentage(20)
sampler = dingo.PolytopeSampler(model)
sampler.generate_steady_states()

Output 1

phase 1: number of correlated samples = 500, effective sample size = 10, ratio of the maximum singilar value over the minimum singular value = 1260.43
phase 2: number of correlated samples = 500, effective sample size = 82, ratio of the maximum singilar value over the minimum singular value = 2.04456
phase 3: number of correlated samples = 2200, effective sample size = 912
[5]total ess 1004: number of correlated samples = 3200
[5]maximum marginal PSRF: 1.00247

phase 1: number of correlated samples = 600, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 3599.46
phase 2: number of correlated samples = 600, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 5630.56
phase 3: number of correlated samples = 600, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 3455.94
phase 4: number of correlated samples = 600, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 3331.84
phase 5: number of correlated samples = 600, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 4463.4
phase 6: number of correlated samples = 600, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 3094.86
phase 7: number of correlated samples = 600, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 3511.64
phase 8: number of correlated samples = 600, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 4926.28
phase 9: number of correlated samples = 600, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 4621.42
phase 10: number of correlated samples = 600, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 4499.42
phase 11: number of correlated samples = 600, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 4172.76
phase 12: number of correlated samples = 600, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 4811.6
phase 13: number of correlated samples = 600, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 3620.84
phase 14: number of correlated samples = 600, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 3519.4
phase 15: number of correlated samples = 600, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 3785.91
phase 16: number of correlated samples = 600, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 3675.98
phase 17: number of correlated samples = 600, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 4776.99
phase 18: number of correlated samples = 600, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 4373.06
phase 19: number of correlated samples = 600, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 5028.11
phase 20: number of correlated samples = 600, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 4325.08
phase 21: number of correlated samples = 2800, effective sample size = 18
phase 22: number of correlated samples = 2800, effective sample size = 17
phase 23: number of correlated samples = 2800, effective sample size = 17
phase 24: number of correlated samples = 2800, effective sample size = 19
phase 25: number of correlated samples = 2800, effective sample size = 17
phase 26: number of correlated samples = 2800, effective sample size = 18
phase 27: number of correlated samples = 2800, effective sample size = 17
phase 28: number of correlated samples = 2800, effective sample size = 17
phase 29: number of correlated samples = 2800, effective sample size = 18
phase 30: number of correlated samples = 2800, effective sample size = 18
phase 31: number of correlated samples = 2800, effective sample size = 18
phase 32: number of correlated samples = 2800, effective sample size = 17
phase 33: number of correlated samples = 2800, effective sample size = 17
phase 34: number of correlated samples = 2800, effective sample size = 18
phase 35: number of correlated samples = 2800, effective sample size = 18
phase 36: number of correlated samples = 2800, effective sample size = 18
phase 37: number of correlated samples = 2800, effective sample size = 17
phase 38: number of correlated samples = 2800, effective sample size = 18
phase 39: number of correlated samples = 2800, effective sample size = 18
phase 40: number of correlated samples = 2800, effective sample size = 18
phase 41: number of correlated samples = 2800, effective sample size = 17
phase 42: number of correlated samples = 2800, effective sample size = 18
phase 43: number of correlated samples = 2800, effective sample size = 17
phase 44: number of correlated samples = 2800, effective sample size = 18
phase 45: number of correlated samples = 2800, effective sample size = 17
phase 46: number of correlated samples = 2800, effective sample size = 18
phase 47: number of correlated samples = 2800, effective sample size = 18
phase 48: number of correlated samples = 2800, effective sample size = 18
phase 49: number of correlated samples = 2800, effective sample size = 18
phase 50: number of correlated samples = 2800, effective sample size = 18
phase 51: number of correlated samples = 2800, effective sample size = 18
phase 52: number of correlated samples = 2800, effective sample size = 18
phase 53: number of correlated samples = 2800, effective sample size = 18
phase 54: number of correlated samples = 2800, effective sample size = 17
phase 55: number of correlated samples = 2800, effective sample size = 17
phase 56: number of correlated samples = 2800, effective sample size = 17
phase 57: number of correlated samples = 2800, effective sample size = 18
phase 58: number of correlated samples = 2800, effective sample size = 17
phase 59: number of correlated samples = 2800, effective sample size = 18
phase 60: number of correlated samples = 2800, effective sample size = 17
phase 61: number of correlated samples = 2800, effective sample size = 18
phase 62: number of correlated samples = 2800, effective sample size = 18
phase 63: number of correlated samples = 2800, effective sample size = 17
phase 64: number of correlated samples = 2800, effective sample size = 18
phase 65: number of correlated samples = 2800, effective sample size = 16
phase 66: number of correlated samples = 2800, effective sample size = 17
phase 67: number of correlated samples = 2800, effective sample size = 18
phase 68: number of correlated samples = 2800, effective sample size = 19
phase 69: number of correlated samples = 2800, effective sample size = 18
phase 70: number of correlated samples = 2800, effective sample size = 17
phase 71: number of correlated samples = 2800, effective sample size = 18
phase 72: number of correlated samples = 2800, effective sample size = 18
phase 73: number of correlated samples = 2800, effective sample size = 17
phase 74: number of correlated samples = 1300, effective sample size = 8
[5]total ess 1002: number of correlated samples = 161700
[5]maximum marginal PSRF: 1.01695


Segmentation fault (core dumped)

As pointed out by @rkstu, a fix to the code runs the code but causes a Segmentation fault at the end. On my end, even on not making any changes results in the same.

Using optimal percentage = 90/20 %:

Output 2

phase 1: number of correlated samples = 500, effective sample size = 3, ratio of the maximum singilar value over the minimum singular value = 3031.23
phase 2: number of correlated samples = 500, effective sample size = 115, ratio of the maximum singilar value over the minimum singular value = 5.55305
phase 3: number of correlated samples = 500, effective sample size = 157, ratio of the maximum singilar value over the minimum singular value = 2.59314
phase 4: number of correlated samples = 1700, effective sample size = 741
[5]total ess 1016: number of correlated samples = 3200
[5]maximum marginal PSRF: 1.0461


Traceback (most recent call last):
File "/mnt/d/Python/GSoC 2025/dingo/dingo/pyoptinterface_based_impl.py", line 378, in remove_redundant_facets
	model_iter.optimize()
File "/home/aditya05/.cache/pypoetry/virtualenvs/dingo-HlYi8y3j-py3.8/lib/python3.8/site-packages/pyoptinterface/_src/highs.py", line 438, in optimize
	super().optimize()
MemoryError: std::bad_alloc

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "test.py", line 11, in <module>
	sampler.generate_steady_states()
File "/mnt/d/Python/GSoC 2025/dingo/dingo/PolytopeSampler.py", line 139, in generate_steady_states
	self.get_polytope()
File "/mnt/d/Python/GSoC 2025/dingo/dingo/PolytopeSampler.py", line 72, in get_polytope
	A, b, Aeq, beq = remove_redundant_facets(
File "/mnt/d/Python/GSoC 2025/dingo/dingo/pyoptinterface_based_impl.py", line 492, in remove_redundant_facets
	except poi.TerminationStatusCode.NUMERICAL_ERROR as e:
TypeError: catching classes that do not inherit from BaseException is not allowed

The most likely reason for this error would be that using 90% as the optimal percentage and then generate_steady_states() followed by any percentage and generate_steady_states() takes up most of the available RAM (~85-95%).
Another evidence of this would be that using del model, sampler and then reinitializing the mode and sampler and then running generate_steady_states() works perfectly fine.

Using optimal percentage = 20/10/5 % (i.e. running generate_steady_states() 3 times) results in the following error in the end

double free or corruption (!prev)
Aborted (core dumped)

Thirdly, I ran the method generate_steady_states() twice without changing the optimal percentage. The first time, it ran only 3-4 phases. The second time, it ran >70 phases. Why does it change so much?