Week 1: Coding Period Begins

Hi everyone,

The official coding period has now begun, and I've diving deeper into the LSODA internals. Continuing from the last blog post, I focused primarily on refactoring helper methods within the LSODA package, while also setting the foundation for a robust and modular design. 

After reading some of the core literature on LSODA, I decided to examine the libLSODA codebase—the C library of LSODA translated from its original Fortran implementation in ODEPACK. My first step was to outline the overall program flow of the LSODA algorithm, which helped me visualize where and how different components fit together.  

Upon exploring the main lsoda function, I chose to start by implementing and refining the helper functions, as they form the building blocks of the algorithm. I thoroughly reviewed the SBSCL LSODA implementation so far and identified areas for improvement.

Here's a more detailed explanation of my work (PR), followed by what I did this week.

1. Made the SM1 array a final variable in LSODACommon, as it is not meant to change during execution. This also allowed the removal of the now-unnecessary setSM1() function.
   
   
2. The LSODAFunction wasn’t yet implemented in SBSCL. Its design was crucial as it determines how the algorithm will integrate with the SBML parsing pipeline.
   To address this, I carefully reviewed both the requirements of LSODAFunction and the design of the SBML integration mechanism within SBSCL.
   Based on that, I decided to use the existing DESystem interface in place of LSODAFunction, as it already provides the required functionalities—such as getDimension() and computeDerivatives(). DESystem itself implements FirstOrderDifferentialEquations from Apache Commons Math, making it a logical and robust fit.
   
   
3. Introduced a new odeSystem variable (of type DESystem) in the LSODAContext class, and added its setter and getter methods.
   
   
4. Refactored all methods previously dependent on LSODAFunction to work with DESystem, including: prja and correction in LSODAIntegrator class and stoda in LSODAStepper class. Because LSODA internally uses 1-based indexing, but computeDerivatives() assumes 0-based indexing, I implemented a utility method called findDerivatives() to handle this translation. It accepts a 1-based y array and returns a 1-based derivative array using temporary conversions internally. 
   
   
5. I also reviewed and refactored the following methods:
   ddot, dscal, daxpy, scaleh, idamax, intdy, dgefa, dgesl, solsy, vmnorm, fnorm, ewset, cfode, corfailure, resetCoeff, and endStoda, addressing subtle bugs and inconsistencies.
   
   
6. Following all this refactoring, I worked on the unit tests as suggested by mentors in last meeting. I ensured that all existing unit tests in LSODAIntegratorTest pass successfully. 
   
   
7. I added new unit tests where I verified the correctness of linear system solutions using dgefa (LU decomposition) and dgesl (solving Ax = B). 
   
   
8. I refactored tests for prja and correction, stubbing the DESystem interface where needed. I calculated expected results manually and validated them against the method outputs.
   
   
9. The LSODAIntegratorTest class now includes 111 passing tests.
   

Meeting with Mentors

• We had our weekly meeting on 5 June 2025, where we continued our discussion on the PR. I shared my progress in the week and communicated my thought process behind the changes.
• Dr. Drager, based on his discussions with Dr. Funahashi, emphasized the importance of rigrous testing. We discussed strategies to cover extreme input values (very large and small), edge cases, exception handling, NaN and invalid inputs etc.
• Arthur will be sharing additional reading materials and resources related to the algorithm. He will also review the PR and let me know of any changes. 
• Dr. Drager told me to post a meeting reminder one day ahead of time every week. 

Next Up

Moving forward, I will be working on the implementation of orderswitch() and methodswitch() functions. These are essential for enabling LSODA to switch between Adams and BDF methods dynamically, based on stiffness detection.

Until next time,
Bye!