- return steady state jacobian

js/copasi.d.ts

@@ -96,6 +96,14 @@ export default class COPASI {
     readonly globalParameterValues : number[];
     readonly localParameterNames : string[];
     readonly localParameterValues : number[];
+    readonly jacobian: object;
+    readonly jacobian2D: number[][];
+    readonly eigenValues: number[][];
+    readonly reducedJacobian: object;
+    readonly reducedJacobian2D: number[][];
+    readonly eigenValuesReduced2D: number[][];
+
+
 }
 
 export {COPASI};

js/copasi.js

@@ -366,6 +366,48 @@ class COPASI {
     get selectedValues() {
         return this._vectorToArray(this.Module.getSelectedValues());
     }
+
+    /**
+     * @property {object} jacobian the jacobian as object
+     */
+    get jacobian() {
+        return JSON.parse(this.Module.getJacbian());
+    }
+
+    /**
+     * @property {number[][]} jacobian2D returns the jacobian as 2D array
+     */
+    get jacobian2D() {
+        return this._vector2dToArray(this.Module.getJacobian2D());
+    }
+
+    /**
+     * @property { number[][]} eigenValues2D returns the eigenvalues as 2D array
+     */
+    get eigenValues2D() {
+        return this._vector2dToArray(this.Module.getEigenValues2D());
+    }
+
+    /**
+     * @property {object} reducedJacobian the reduced jacobian as object
+     */
+    get reducedJacobian() {
+        return JSON.parse(this.Module.getJacobianReduced());
+    }
+
+    /**
+     * @property {number[][]} reducedJacobian2D returns the reduced jacobian as 2D array
+     */
+    get reducedJacobian2D() {
+        return this._vector2dToArray(this.Module.getJacobianReduced2D());
+    }
+
+    /**
+     * @property { number[][]} eigenValuesReduced2D returns the eigenvalues of the reduced Jacobian as 2D array
+     */
+    get eigenValuesReduced2D() {
+        return this._vector2dToArray(this.Module.getEigenValuesReduced2D());
+    }
 }
 
 // if module is defined, export the COPASI class

src/copasijs.cpp

@@ -7,6 +7,9 @@ using namespace emscripten;
 
 #include "copasijs.h"
 
+#include <copasi/core/CDataContainer.h>
+#include <copasi/steadystate/CEigen.h>
+
 #ifdef __cplusplus
 #define EXTERN extern "C"
 #else
@@ -1007,12 +1010,138 @@ std::vector<std::vector<double>> getSimulationResults2D()
     return results;
 }
 
+std::vector<std::vector<double>> convertCEigen(const CEigen& eValues)
+{
+    std::vector<std::vector<double>> values;
+
+    auto& imag = eValues.getI();
+    auto& real = eValues.getR();
+
+    for (size_t i = 0; i < imag.size(); ++i)
+    {
+        std::vector<double> row;
+        row.push_back(real[i]);
+        row.push_back(imag[i]);
+        values.push_back(row);
+    }
+
+    return values;
+}
+
+std::vector<std::vector<double>> convertCArray( CArrayInterface* pArray)
+{
+    std::vector<std::vector<double>> values;
+    if (pArray == NULL || pArray->dimensionality() != 2)
+        return values;
+
+    auto arrayDim = pArray->size();
+
+    for (size_t i = 0; i < arrayDim[0]; ++i)
+    {
+        std::vector<double> row;
+        for (size_t j = 0; j < arrayDim[1]; ++j)
+        {
+            CArrayInterface::index_type idx = {i, j};
+            row.push_back((*pArray)[idx]);
+        }
+        values.push_back(row);
+    }
+
+    return values;
+}
+
+ordered_json convertDataArray(const CDataArray* pArray)
+{
+    ordered_json result;
+
+    if (pArray == NULL || pArray->dimensionality() != 2)
+        return result;
+
+    std::vector<std::string> columns = pArray->getAnnotationsString(1);
+    std::vector<std::string> rows = pArray->getAnnotationsString(0);
+    std::vector<std::vector<double>> values = convertCArray(const_cast<CArrayInterface*>( pArray->getArray()));
+
+    result["columns"] = columns;
+    result["rows"] = rows;
+    result["values"] = values;
+
+    return result;
+}
+
+std::string getJacobian()
+{
+    ensureModel();
+
+    auto &task = dynamic_cast<CSteadyStateTask &>((*pDataModel->getTaskList())["Steady-State"]);
+    auto* pMatrix = task.getJacobianAnnotated();
+
+    return convertDataArray(pMatrix).dump(2);
+
+}
+
+std::vector<std::vector<double>> getJacobian2D()
+{
+    ensureModel();
+
+    auto &task = dynamic_cast<CSteadyStateTask &>((*pDataModel->getTaskList())["Steady-State"]);
+    auto* pMatrix = task.getJacobianAnnotated();
+
+    return convertCArray(pMatrix->getArray());
+}
+
+std::vector<std::vector<double>> getEigenValues2D()
+{
+    ensureModel();
+
+    auto &task = dynamic_cast<CSteadyStateTask &>((*pDataModel->getTaskList())["Steady-State"]);
+    const auto& eValues = task.getEigenValues();
+
+    return convertCEigen(eValues);
+}
+
+std::string getJacobianReduced()
+{
+    ensureModel();
+
+    auto &task = dynamic_cast<CSteadyStateTask &>((*pDataModel->getTaskList())["Steady-State"]);
+    auto* pMatrix = task.getJacobianXAnnotated();
+
+    return convertDataArray(pMatrix).dump(2);
+
+}
+
+std::vector<std::vector<double>> getJacobianReduced2D()
+{
+    ensureModel();
+
+    auto &task = dynamic_cast<CSteadyStateTask &>((*pDataModel->getTaskList())["Steady-State"]);
+    auto* pMatrix = task.getJacobianXAnnotated();
+
+    return convertCArray(pMatrix->getArray());
+}
+
+std::vector<std::vector<double>> getEigenValuesReduced2D()
+{
+    ensureModel();
+
+    auto &task = dynamic_cast<CSteadyStateTask &>((*pDataModel->getTaskList())["Steady-State"]);
+    const auto& eValues = task.getEigenValuesReduced();
+
+    return convertCEigen(eValues);
+}
+
+
 double steadyState()
 {
     ensureModel();
 
     auto &task = dynamic_cast<CSteadyStateTask &>((*pDataModel->getTaskList())["Steady-State"]);
+
     task.setUpdateModel(true);
+    auto *problem = dynamic_cast<CSteadyStateProblem *>(task.getProblem());
+    bool request = true;
+    problem->setStabilityAnalysisRequested(request);
+
 
     if (!task.initialize(CCopasiTask::OUTPUT_UI, pDataModel, NULL))
         return std::numeric_limits<double>::quiet_NaN();
@@ -1431,5 +1560,12 @@ EMSCRIPTEN_BINDINGS(copasi_binding)
     emscripten::function("getSelectedValues", &getSelectionValues);
     emscripten::function("setSelectionList", &setSelectionList);
 
+    emscripten::function("getJacobian", &getJacobian);
+    emscripten::function("getJacobian2D", &getJacobian2D);
+    emscripten::function("getEigenValues2D", &getEigenValues2D);
+    emscripten::function("getJacobianReduced", &getJacobianReduced);
+    emscripten::function("getJacobianReduced2D", &getJacobianReduced2D);
+    emscripten::function("getEigenValuesReduced2D", &getEigenValuesReduced2D);
+
 }
 #endif

src/copasijs.h

@@ -351,6 +351,53 @@ void setValue(const std::string& nameOrId, double value);
 /// sets the value by display name
 void setValueByName(const std::string& name, double value);
 
+/// @brief returns the Jacobian at steady state as JSON string
+/// @return the Jacobian as JSON string in the following format
+///
+/// ```json
+/// {
+///     "rows": ["X", "Y", ...],
+///     "columns": ["X", "Y", ...],
+///     "values": [
+///         [1.0, 0.0, ...],
+///         [0.0, 1.0, ...],
+///         ...
+///     ]
+/// }
+/// ```
+std::string getJacobian();
+
+/// @brief returns the Jacobian at steady state as 2D double vector
+std::vector<std::vector<double>> getJacobian2D();
+
+/// @brief returns the eigenvalues of the Jacobian at steady state as 2D double vector
+std::vector<std::vector<double>> getEigenValues2D();
+
+/// @brief returns the reduced Jacobian at steady state as JSON string
+/// @return the Jacobian as JSON string in the following format
+///
+/// ```json
+/// {
+///     "rows": ["X", "Y", ...],
+///     "columns": ["X", "Y", ...],
+///     "values": [
+///         [1.0, 0.0, ...],
+///         [0.0, 1.0, ...],
+///         ...
+///     ]
+/// }
+/// ```
+std::string getJacobianReduced();
+
+/// @brief returns the reduced Jacobian at steady state as 2D double vector
+std::vector<std::vector<double>> getJacobianReduced2D();
+
+/// @brief returns the reduced eigenvalues of the Jacobian at steady state as 2D double vector
+std::vector<std::vector<double>> getEigenValuesReduced2D();
+
+
+
+
 #pragma region  // internal calls Internal
 
 /// @brief frees a pointer allocated by the COPASI library

test/copasijs_test.cpp

@@ -328,3 +328,41 @@ TEST_CASE("Load COVID Model", "[copasijs][slow]")
     CAPTURE(json["titles"][0].size());
     REQUIRE(json["columns"][0].size() == json["recorded_steps"].get<int>());
 }
+
+
+
+TEST_CASE("Steadystate test", "[copasijs][steadystate]") 
+{
+    Instance instance;
+    std::string model = loadFromFile(getTestFile("../example_files/brusselator.cps"));
+    REQUIRE(!model.empty());
+    REQUIRE(model != "Error loading model");
+
+    auto selectionList = getSelectionList();
+    REQUIRE_THAT(selectionList, Catch::Matchers::Equals(std::vector<std::string>{"Time", "X", "Y"}));
+    auto selectionValues = getSelectionValues();
+    REQUIRE(selectionValues.size() == 3);
+    REQUIRE_THAT(selectionValues, Catch::Matchers::Approx(std::vector<double>{0, 3, 3}));
+
+    // run to steady state
+    double closeNess = steadyState();
+    REQUIRE(closeNess < 0.1);
+
+    // look at steady state values
+    selectionValues = getSelectionValues();
+    REQUIRE(selectionValues.size() == 3);
+    REQUIRE_THAT(selectionValues, Catch::Matchers::Approx(std::vector<double>{0, 0.5, 6}));
+
+    // get the jacobian
+    auto jacobian = getJacobian();
+    auto jac2d = getJacobian2D();    
+    auto eigenValues = getEigenValues2D();
+
+    REQUIRE(jac2d.size() == 2);
+
+    jacobian = getJacobianReduced();
+    jac2d = getJacobianReduced2D();
+    eigenValues = getEigenValuesReduced2D();
+
+    REQUIRE(jac2d.size() == 2);
+}