Add comprehensive biohydrogenation debugging and CAS dependencies

Major changes:
- Add biohydrogenation example with extensive debugging
- Add CAS dependencies (Oscar, Singular, Groebner, Nemo, DynamicPolynomials)
- Extensive debugging for identifiability analysis
- Use Symbolics.substitute explicitly throughout
- Add debug output for polynomial system construction
- Track derivative order decisions and rank analysis

src/core/parameter_estimation.jl:
- Add comprehensive debug output for identifiability analysis
- Track derivative order selection (PE formula vs heuristic)
- Debug Jacobian construction and nullspace analysis
- Output full polynomial system for debugging
- Add UNTESTED FIX comments for OrderedDict usage
- Add EXPERIMENTAL FIX for derivative order calculation

src/examples/biohydrogenation/:
- Complete biohydrogenation example with saved systems
- Algebraic diagnostics summary
- Data files and results

Dependencies added for CAS diagnostics:
- Oscar, Singular, Groebner for algebraic analysis
- DynamicPolynomials, MultivariatePolynomials for polynomial manipulation
- Nemo for number theory support

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>

Author: orebas

Project.toml

@@ -10,30 +10,36 @@ CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GaussianProcesses = "891a1506-143c-57d2-908e-e1f8e92e6de9"
+Groebner = "0b43b601-686d-58a3-8a1c-6623616c7cd4"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LeastSquaresOptim = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 NLSolversBase = "d41bc354-129a-5804-8e4c-c37616107c6c"
+Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Oscar = "f1435218-dba5-11e9-1e4d-f1a5fab5fc13"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PolynomialRoots = "3a141323-8675-5d76-9d11-e1df1406c778"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 RS = "01933e51-4c7d-4354-b809-e991ab46e7ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RationalUnivariateRepresentation = "17246082-4fa7-4287-87d4-6a5617fca855"
+Singular = "bcd08a7b-43d2-5ff7-b6d4-c458787f915c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
@@ -48,6 +54,7 @@ PEtab = "48d54b35-e43e-4a66-a5a1-dde6b987cf69"
 ODEParameterEstimationPEtabExt = ["PEtab"]
 
 [compat]
+Oscar = "1.2.2"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/ODEParameterEstimation.jl

@@ -28,6 +28,11 @@ using NonlinearSolve
 using PolynomialRoots
 using Suppressor
 using Logging
+using Singular
+using Groebner
+using DynamicPolynomials
+using Nemo#using GLPK
+using Oscar
 
 using AbstractAlgebra
 using RationalUnivariateRepresentation
@@ -129,7 +134,7 @@ export substr_test, global_unident_test, sum_test, trivial_unident
 		# data_sample = load("/home/ad7760/parameter_estimation_tests/data/julia/lotka-volterra_0.jld2", "data")
 
 		estimation_problem = sample_problem_data(pep, datasize = datasize, time_interval = time_interval, noise_level = 0.0)
-		res = analyze_parameter_estimation_problem(estimation_problem, nooutput = true, system_solver = solve_with_nlopt)
+		res = analyze_parameter_estimation_problem(estimation_problem, nooutput = true, system_solver = solve_with_nlopt, shooting_points = 1)
 		#analysis_result, besterror = 
 		#	analyze_estimation_result(estimation_problem, res, nooutput = true)