.gitattributes

+# Set the deifault behavior, in case people don't have core.autocrlf set.
+*        text=auto
+
+# Explicitly declare text files you want to always be normalized and converted
+# to native line endings on checkout.
+
+# C/C++
+*.c      text
+*.cpp    text
+*.cxx    text
+*.cc     text
+*.h      text
+*.hpp    text
+
+# Python
+*.py     text
+
+# Object Pascal
+*.pas    text
+
+# Markdown
+*.md     text
+
+# CMake
+*.cmake  text
+
+# Data
+*.xml    text
+*.dat    text
+
+# Misc Windows files
+*.def    text
+*.txt    text
+*.bat    text
+
+# Doxygen
+*.doxy   text
+
+# Installer scripts (Inno?)
+*.iss    text
+
+# RoadRunner custom
+*.rrtest text

.travis.yml

@@ -25,13 +25,15 @@ before_install:
   # install dependencies
   - sudo apt-get -y install cmake swig python-dev python-numpy llvm-3.5-dev
 
-  - mkdir ~/install
-
   # install roadrunner dependencies via apt:
-  - sudo sh -c 'echo "\n\ndeb http://cdn.rawgit.com/sys-bio/roadrunner-debs/master ./ \n\n" >> /etc/apt/sources.list'
-  - sudo cat /etc/apt/sources.list
-  - sudo apt-get update -qq
-  - sudo apt-get -y --force-yes install libroadrunner-deps-dev
+  - mkdir ~/debs
+  - cd ~/debs
+  - wget http://cdn.rawgit.com/sys-bio/roadrunner-debs/b8288d2f3e9623ad0f13ff8a1aa97fabf481413a/libroadrunner-deps_1.5_amd64.deb
+  - sudo dpkg -i libroadrunner-deps_1.5_amd64.deb
+  - sudo cp /usr/local/include/rr-libstruct/* /usr/local/include
+  - sudo cp /usr/local/include/clapack/* /usr/local/include
+
+  - mkdir ~/install
 
   - echo "TRAVIS_BUILD_DIR $TRAVIS_BUILD_DIR"
   # build roadrunner:

CMakeLists.txt

@@ -5,8 +5,8 @@ project(rr)
 # source files expose library version information
 
 set(ROADRUNNER_VERSION_MAJOR 1)
-set(ROADRUNNER_VERSION_MINOR 4)
-set(ROADRUNNER_VERSION_PATCH 8)
+set(ROADRUNNER_VERSION_MINOR 5)
+set(ROADRUNNER_VERSION_PATCH 1)
 
 set(ROADRUNNER_VERSION "${ROADRUNNER_VERSION_MAJOR}.${ROADRUNNER_VERSION_MINOR}.${ROADRUNNER_VERSION_PATCH}")
 
@@ -208,7 +208,7 @@ else()
 endif()
 
 find_library(LIBSBML_LIBRARY sbml libsbml HINTS ${THIRD_PARTY_INSTALL_FOLDER}/lib)
-find_library(LIBSBML_STATIC_LIBRARY sbml-static libsbml-static HINTS ${THIRD_PARTY_INSTALL_FOLDER}/lib)
+find_library(LIBSBML_STATIC_LIBRARY sbml-static libsbml-static HINTS ${THIRD_PARTY_INSTALL_FOLDER}/lib ${THIRD_PARTY_INSTALL_FOLDER}/lib64)
 find_path(LIBSBML_INCLUDE_DIR sbml/SBMLTypes.h HINTS ${THIRD_PARTY_INSTALL_FOLDER}/include)
 
 # find libSBML deps
@@ -241,6 +241,7 @@ if (${BUILD_LLVM})
         ${LIBRARY_OUTPUT_PATH}
         ${THIRD_PARTY_FOLDER}/dependencies/libsbml/lib
         ${THIRD_PARTY_INSTALL_FOLDER}/lib
+        ${THIRD_PARTY_INSTALL_FOLDER}/lib64
         ${LLVM_LIBRARY_DIRS}
     )
 else()
@@ -248,6 +249,7 @@ else()
         ${LIBRARY_OUTPUT_PATH}
         ${THIRD_PARTY_FOLDER}/dependencies/libsbml/lib
         ${THIRD_PARTY_INSTALL_FOLDER}/lib
+        ${THIRD_PARTY_INSTALL_FOLDER}/lib64
     )
 endif (${BUILD_LLVM})
 

LICENSE.txt

 Licensing and Distribution Terms for RoadRunner
- * Copyright (C) 2012-2013
+ * Copyright (C) 2012-2018
  *   University of Washington, Seattle, WA, USA
  *
  * Licensed under the Apache License, Version 2.0 (the "License");

README.md

@@ -2,7 +2,7 @@
 [![GitHub version](https://badge.fury.io/gh/sys-bio%2Froadrunner.svg)](http://badge.fury.io/gh/sys-bio%2Froadrunner)
 [![Build Status](https://travis-ci.org/sys-bio/roadrunner.svg?branch=develop)](https://travis-ci.org/sys-bio/roadrunner)
 
-Copyright 2013-2016
+Copyright 2013-2018
 
 E. T. Somogyi <sup>1</sup>, J. K. Medley <sup>3</sup>, M. T. Karlsson <sup>2</sup>, M. Swat <sup>1</sup>, M. Galdzicki <sup>3</sup>, K. Choi <sup>3</sup>, W. Copeland <sup>3</sup> and H. M. Sauro <sup>3</sup>
 

debian/changelog

-libroadrunner (1.4.8-1) UNRELEASED; urgency=medium
+libroadrunner (1.5.1-1) UNRELEASED; urgency=medium
 
   * Initial release (Closes: #<bug>)
 
- -- Andreas Tille <tille@debian.org>  Tue, 31 Jan 2017 09:50:22 +0100
+ -- Andreas Tille <tille@debian.org>  Mon, 29 Oct 2018 08:45:58 +0100

debian/compat

-10
+11

debian/control

@@ -2,8 +2,9 @@ Source: libroadrunner
 Maintainer: Debian Med Packaging Team <debian-med-packaging@lists.alioth.debian.org>
 Uploaders: Andreas Tille <tille@debian.org>
 Section: science
+Testsuite: autopkgtest-pkg-python
 Priority: optional
-Build-Depends: debhelper (>= 10),
+Build-Depends: debhelper (>= 11~),
                llvm-dev,
                cmake,
                libncurses5-dev,
@@ -15,11 +16,11 @@ Build-Depends: debhelper (>= 10),
                swig,
                libsbml5-dev,
                libxml2-dev,
-               libsundials-serial-dev,
+               libsundials-dev,
                libpoco-dev
-Standards-Version: 3.9.8
-Vcs-Browser: https://anonscm.debian.org/cgit/debian-med/libroadrunner.git
-Vcs-Git: https://anonscm.debian.org/git/debian-med/libroadrunner.git
+Standards-Version: 4.2.1
+Vcs-Browser: https://salsa.debian.org/med-team/libroadrunner
+Vcs-Git: https://salsa.debian.org/med-team/libroadrunner.git
 Homepage: http://libroadrunner.org/
 
 Package: libroadrunner
@@ -38,8 +39,9 @@ Description: simulation engine for systems and synthetic biology
  simulation support.
 
 Package: libroadrunner-dev
-Section: libdevel
 Architecture: any
+Section: libdevel
+Testsuite: autopkgtest-pkg-python
 Depends: libroadrunner (= ${binary:Version})
 Description: simulation engine for systems and synthetic biology (development)
  LibRoadRunner is a high performance and portable simulation engine for
@@ -55,8 +57,9 @@ Description: simulation engine for systems and synthetic biology (development)
  This package contains the header files for libroadrunner.
 
 Package: python-roadrunner
-Section: python
 Architecture: any
+Section: python
+Testsuite: autopkgtest-pkg-python
 Depends: libroadrunner (= ${binary:Version}),
          python-numpy
 Description: Python bindings for libroadrunner

packaging/deb/libroadrunner/debian/changelog

+libroadrunner (1.4.15) unstable; urgency=low
+
+  * Sync with Github
+
+ -- Kyle Medley <0u812github@gmail.com>  Wed, 23 Mar 2016 12:00:00 -0700
+
+libroadrunner (1.4.4) unstable; urgency=low
+
+  * Add raterules to ss calc
+
+ -- Kyle Medley <0u812github@gmail.com>  Wed, 23 Mar 2016 12:00:00 -0700
+
 libroadrunner (1.4.2) unstable; urgency=low
 
   * Bufix release

packaging/deb/libroadrunner/debian/rules

@@ -17,4 +17,4 @@ CFLAGS   = -Wno-error
 CPPFLAGS = -Wno-error
 CXXFLAGS = -Wno-error
 
-DEB_CMAKE_NORMAL_ARGS += -DLLVM_CONFIG_EXECUTABLE=/usr/lib/llvm-3.5/bin/llvm-config -DTHIRD_PARTY_INSTALL_FOLDER=/home/user/devel/install/roadrunner-deps -DBUILD_PYTHON=TRUE
+DEB_CMAKE_NORMAL_ARGS += -DLLVM_CONFIG_EXECUTABLE=/usr/lib/llvm-3.5/bin/llvm-config -DTHIRD_PARTY_INSTALL_FOLDER=/usr -DBUILD_PYTHON=TRUE

source/CMakeLists.txt

@@ -278,8 +278,10 @@ if(RR_BUILD_SHARED_CORE)
             LIBRARY DESTINATION lib COMPONENT rr_core
             ARCHIVE DESTINATION lib COMPONENT rr_core
             )
+        if(${BUILD_PYTHON})
 			install(TARGETS roadrunner DESTINATION ${PYTHON_PACKAGE_DEST_DIR})
 		endif()
+    endif()
 
 
     if(BUILD_LLVM)

source/CVODEIntegrator.cpp

@@ -20,6 +20,8 @@
 #include <assert.h>
 #include <Poco/Logger.h>
 
+#include <iostream>
+
 #define CVODE_INT_TYPECODE 0x7799ff00
 
 using namespace std;
@@ -94,8 +96,8 @@ namespace rr
         Solver::resetSettings();
 
         // Set default integrator settings.
-        addSetting("relative_tolerance", 1e-6, "Relative Tolerance", "Specifies the scalar relative tolerance (double).", "(double) CVODE calculates a vector of error weights which is used in all error and convergence tests. The weighted RMS norm for the relative tolerance should not become smaller than this value.");
-        addSetting("absolute_tolerance", 1e-15, "Absolute Tolerance", "Specifies the scalar absolute tolerance (double).", "(double) CVODE calculates a vector of error weights which is used in all error and convergence tests. The weighted RMS norm for the absolute tolerance should not become smaller than this value.");
+        addSetting("relative_tolerance", Config::getDouble(Config::CVODE_MIN_RELATIVE), "Relative Tolerance", "Specifies the scalar relative tolerance (double).", "(double) CVODE calculates a vector of error weights which is used in all error and convergence tests. The weighted RMS norm for the relative tolerance should not become smaller than this value.");
+        addSetting("absolute_tolerance", Config::getDouble(Config::CVODE_MIN_ABSOLUTE), "Absolute Tolerance", "Specifies the scalar absolute tolerance (double).", "(double) CVODE calculates a vector of error weights which is used in all error and convergence tests. The weighted RMS norm for the absolute tolerance should not become smaller than this value.");
         addSetting("stiff",              true, "Stiff", "Specifies whether the integrator attempts to solve stiff equations. (bool)", "(bool) Specifies whether the integrator attempts to solve stiff equations. Ensure the integrator can solver stiff differential equations by setting this value to true.");
         addSetting("maximum_bdf_order",  mDefaultMaxBDFOrder, "Maximum BDF Order", "Specifies the maximum order for Backward Differentiation Formula integration. (int)", "(int) Specifies the maximum order for Backward Differentiation Formula integration. This integration method is used for stiff problems. Default value is 5.");
         addSetting("maximum_adams_order",mDefaultMaxAdamsOrder, "Maximum Adams Order", "Specifies the maximum order for Adams-Moulton intergration. (int)", "(int) Specifies the maximum order for Adams-Moulton intergration. This integration method is used for non-stiff problems. Default value is 12.");
@@ -198,13 +200,6 @@ namespace rr
         // MULTIPLE STEPS
         bVal = Config::getBool(Config::SIMULATEOPTIONS_MULTI_STEP);
         Integrator::setValue("multiple_steps", bVal);
-
-        // ABSOLUTE TOLERANCE
-//         Integrator::setValue("absolute_tolerance", Config::getDouble(Config::SIMULATEOPTIONS_ABSOLUTE));
-//         Integrator::setValue("relative_tolerance", Config::getDouble(Config::SIMULATEOPTIONS_RELATIVE));
-
-        CVODEIntegrator::setValue("absolute_tolerance", Config::getDouble(Config::CVODE_MIN_ABSOLUTE));
-        CVODEIntegrator::setValue("relative_tolerance", Config::getDouble(Config::CVODE_MIN_RELATIVE));
     }
 
 	void CVODEIntegrator::loadSBMLSettings(std::string const& filename)
@@ -734,6 +729,8 @@ namespace rr
 			return;
 		}
 
+		lastEventTime = 0.0;
+
 		// apply any events that trigger before or at time 0.
 		// important NOT to set model time before we check get
 		// the initial event state, initially time is < 0.

source/EulerIntegrator.h

@@ -11,9 +11,11 @@
 #include "Integrator.h"
 #include "rrExecutableModel.h"
 #include "rrRoadRunnerOptions.h"
+
 #include <string>
 #include <stdexcept>
 #include <sstream>
+#include <iostream>
 
 namespace rr
 {
@@ -70,22 +72,28 @@ namespace rr
         * @param o: a reference to a SimulatOptions object where the configuration
         * parameters will be read from.
         */
-        EulerIntegrator(ExecutableModel *m) {
+        EulerIntegrator(ExecutableModel *m)
+			:	eventStatus(std::vector<unsigned char>(m->getNumEvents(),false)),
+				previousEventStatus(std::vector<unsigned char>(m->getNumEvents(), false)) {
+            resetSettings();
+
             model = m;
             exampleParameter1 = 3.14;
             exampleParameter2 = "hello";
+			Log(Logger::LOG_WARNING) << "Euler integrator is inaccurate";
+			//std::cerr << "Number of event triggers: " << m->getEventTriggers(0, 0, 0) << std::endl;
 
             if(model) {
                 // calling the getStateVector with a NULL argument returns
                 // the size of teh state vector.
                 stateVectorSize = model->getStateVector(NULL);
                 rateBuffer = new double[stateVectorSize];
-                stateBuffer1 = new double[stateVectorSize];
-                stateBuffer2 = new double[stateVectorSize];
+                stateBufferBegin = new double[stateVectorSize];
+                stateBufferEnd = new double[stateVectorSize];
             } else {
                 rateBuffer = NULL;
-                stateBuffer1 = NULL;
-                stateBuffer2 = NULL;
+                stateBufferBegin = NULL;
+                stateBufferEnd = NULL;
             }
         }
 
@@ -94,8 +102,8 @@ namespace rr
         */
         virtual ~EulerIntegrator() {
             delete[] rateBuffer;
-            delete[] stateBuffer1;
-            delete[] stateBuffer2;
+            delete[] stateBufferBegin;
+            delete[] stateBufferEnd;
         };
 
         /**
@@ -108,7 +116,13 @@ namespace rr
         * @return the end time.
         */
         virtual double integrate(double t0, double h) {
-            if(model) {
+			int internal_steps = getValueAsInt("subdivision_steps");
+			if (model == (rr::ExecutableModel*)NULL) return 0;
+
+            double finalTimeEnd;
+
+			h /= internal_steps;
+			for (int subdiv = 0; subdiv < internal_steps; ++subdiv) {
 				// evaluate and copy the rate of change of the state vector
 				// rate into the local buffer. If the 2nd argument is NULL,
 				// the current model state is used to evaluate the
@@ -116,28 +130,61 @@ namespace rr
 				model->getStateVectorRate(t0, NULL, rateBuffer);
 
 				// copy the current state vector into a local buffer
-                model->getStateVector(stateBuffer1);
+				model->getStateVector(stateBufferBegin);
 
 				// perform the Euler integration step, i.e.
 				// y_{n+1} = y_{n} + h * y'_{n}
 				for (int i = 0; i < stateVectorSize; ++i) {
-                    stateBuffer2[i] = stateBuffer1[i] + h * rateBuffer[i];
+					stateBufferEnd[i] = stateBufferBegin[i] + h * rateBuffer[i];
 				}
 
 				// set the model state to the newly calculated state
-                model->setStateVector(stateBuffer2);
+				model->setStateVector(stateBufferEnd);
 
 				// update the model time to the new time
-                model->setTime(t0 + h);
+				double timeEnd = t0 + h;
+				model->setTime(timeEnd);
 
 				// if we have a client, notify them that we have taken
 				// a time step
-                if (listener)
-                {
-                    listener->onTimeStep(this, model, t0 + h);
+				if (listener) {
+					listener->onTimeStep(this, model, timeEnd);
+				}
+
+				// events
+				bool triggered = false;
+
+				model->getEventTriggers(eventStatus.size(), NULL, eventStatus.size() ? &eventStatus[0] : NULL);
+				for (int k_ = 0; k_ < eventStatus.size(); ++k_) {
+					if (eventStatus.at(k_)) {
+						triggered = true;
+						//std::cerr << "Triggered" << std::endl;
 					}
 				}
-            return t0 + h;
+
+				if (triggered) {
+					// applyEvents takes the list of events which were previously triggered
+					//std::cerr << "An event was triggered at " << t0 << std::endl;
+					applyEvents(timeEnd, previousEventStatus);
+				}
+
+				if (eventStatus.size()) {
+					previousEventStatus = eventStatus;
+				}
+
+                finalTimeEnd = timeEnd;
+			}
+
+            return finalTimeEnd;
+        }
+
+		void applyEvents(double timeEnd, std::vector<unsigned char> &previousEventStatus) {
+			//std::cerr << "Size of previous events: " << previousEventStatus.size() << std::endl;
+			// If we pass in the events including the ones just triggered, they won't be applied, so use previousEventStatus
+			model->applyEvents(timeEnd, previousEventStatus.size() == 0 ? NULL : &previousEventStatus[0], stateBufferEnd, NULL);
+			// AHu: jk I think that model->applyEvents does update the mode's state vector
+			// The previous statement loaded the result into the final stateBufferEnd, so now update the model's state vector
+			//model->setStateVector(stateBufferEnd);
 		}
 
         /**
@@ -234,21 +281,9 @@ namespace rr
          * @brief Get the hint for this integrator
          */
         static std::string getEulerHint() {
-//             return "An elementary (my dear Watson) Euler integrator";
             return "A simple Euler integrator";
         }
 
-        /**
-         * @author JKM
-         * @brief Reset all integrator settings to their respective default values
-         */
-        void resetSettings()
-        {
-            Solver::resetSettings();
-
-            // Euler integrator has no settings
-        }
-
         /**
          * @author JKM
          * @brief Always deterministic for Euler
@@ -270,7 +305,9 @@ namespace rr
             }
 
             if(key == "exampleParameter2") {
+				// Ahu: Why is this cast here, and is this a static or dynamic cast?
                 exampleParameter2 = (string)value;
+				return;
             }
 
             // they did not give a valid key, so throw an exception.
@@ -326,33 +363,14 @@ namespace rr
             return keys;
         }
 
-        /**
-        * list of keys that this integrator supports.
-        *
-        * This method is called by the IntegratorFactory to build a list of
-        * all the options that all the integrators support.
-        */
-        static const Dictionary* getIntegratorOptions() {
-            static BasicDictionary dict;
-
-            // fill out the required options for defining this integrator
-            dict.setItem("integrator", "euler");
-            dict.setItem("integrator.description",
-                    "Forward Euler Integrator: an example of creating an integrator");
-            dict.setItem("integrator.hint", "Forward Euler Integrator");
-
-            // this uses two example parameters, add them here
-            dict.setItem("exampleParameter1", 3.14);
-            dict.setItem("exampleParameter1.description",
-                    "an example parameter that does nothing");
-            dict.setItem("exampleParameter1.hint", "a hint for the example parameter");
-
-            dict.setItem("exampleParameter2", string("hello"));
-            dict.setItem("exampleParameter2.description",
-                    "another parameter which does nothing");
-            dict.setItem("exampleParameter2.hint", "a hint for the other parameter");
+		void resetSettings() {
+			Solver::resetSettings();
 
-            return &dict;
+			// Set default integrator settings.
+            addSetting("subdivision_steps", 1,
+                "Subdivision Steps",
+                "The number of subdivisions of the Euler step size (int).",
+                "(int) For each point, up to this many extra steps will be taken as smaller steps within each step, although their values are not saved");
         }
 
     private:
@@ -375,13 +393,16 @@ namespace rr
         * two buffers to store the state vector rate, and
         * new state vector
         */
-        double *rateBuffer, *stateBuffer1, *stateBuffer2;
+        double *rateBuffer, *stateBufferBegin, *stateBufferEnd;
 
         /**
         * size of state vector
         */
         int stateVectorSize;
 
+		std::vector<unsigned char> eventStatus;
+		std::vector<unsigned char> previousEventStatus;
+
         /**
         * Clients may register a listener to listen for time steps, or
         * sbml events. Time step events are much more usefull for variable

source/IntegratorRegistration.cpp

@@ -35,7 +35,7 @@ namespace rr
         IntegratorFactory::getInstance().registerIntegrator(new GillespieIntegratorRegistrar());
         IntegratorFactory::getInstance().registerIntegrator(new RK4IntegratorRegistrar());
         IntegratorFactory::getInstance().registerIntegrator(new RK45IntegratorRegistrar());
-//         IntegratorFactory::getInstance().registerIntegrator(new EulerIntegratorRegistrar());
+        IntegratorFactory::getInstance().registerIntegrator(new EulerIntegratorRegistrar());
     }
 
     void IntegratorRegistrationMgr::Register() {

source/NLEQSolver.cpp

@@ -3,7 +3,7 @@
 #include "rrExecutableModel.h"
 #include "rrStringUtils.h"
 #include "rrUtils.h"
-#include "nleq/nleq1.h"
+#include "nleq/nleq2.h"
 #include "rrLogger.h"
 #include "rrUtils.h"
 #include "rrException.h"
@@ -39,9 +39,19 @@ void NLEQSolver::loadConfigSettings()
     SteadyStateSolver::loadConfigSettings();
     // Load settings specific to solver integrator
 
-	NLEQSolver::setValue("maximum_iterations", Config::getInt(Config::STEADYSTATE_MAXIMUM_NUM_STEPS));
+    NLEQSolver::setValue("allow_presimulation", Config::getBool(Config::STEADYSTATE_PRESIMULATION));
+    NLEQSolver::setValue("presimulation_tolerance", Config::getDouble(Config::STEADYSTATE_PRESIMULATION_TOL));
+    NLEQSolver::setValue("presimulation_maximum_steps", Config::getInt(Config::STEADYSTATE_PRESIMULATION_MAX_STEPS));
+    NLEQSolver::setValue("presimulation_time", Config::getDouble(Config::STEADYSTATE_PRESIMULATION_TIME));
+    NLEQSolver::setValue("allow_approx", Config::getBool(Config::STEADYSTATE_APPROX));
+    NLEQSolver::setValue("approx_tolerance", Config::getDouble(Config::STEADYSTATE_APPROX_TOL));
+    NLEQSolver::setValue("approx_maximum_steps", Config::getInt(Config::STEADYSTATE_APPROX_MAX_STEPS));
+    NLEQSolver::setValue("approx_time", Config::getDouble(Config::STEADYSTATE_APPROX_TIME));
     NLEQSolver::setValue("relative_tolerance", Config::getDouble(Config::STEADYSTATE_RELATIVE));
+    NLEQSolver::setValue("maximum_iterations", Config::getInt(Config::STEADYSTATE_MAXIMUM_NUM_STEPS));
     NLEQSolver::setValue("minimum_damping", Config::getDouble(Config::STEADYSTATE_MINIMUM_DAMPING));
+    NLEQSolver::setValue("broyden_method", Config::getInt(Config::STEADYSTATE_BROYDEN));
+    NLEQSolver::setValue("linearity", Config::getInt(Config::STEADYSTATE_LINEARITY));
 }
 
 void NLEQSolver::resetSettings()
@@ -49,9 +59,19 @@ void NLEQSolver::resetSettings()
     Solver::resetSettings();
 
     // Set default settings.
+    addSetting("allow_presimulation", false, "Allow Presimulation", "Flag for starting steady state analysis with simulation (bool).", "(bool) This flag does not affect the usage of approximation routine when the default steaty state solver fails");
+    addSetting("presimulation_tolerance", 1e-6, "Presimulation Tolerance", "Tolerance for presimulation before steady state analysis (double).", "(double) Absolute tolerance used by presimulation routine. Only used when allow_presimulation is True");
+    addSetting("presimulation_maximum_steps", 10000, "Presimulation Maximum Steps", "Maximum number of steps that can be taken for presimulation before steady state analysis (int).", "(int) Takes priority over presimulation_time. Only used when allow_presimulation is True");
+    addSetting("presimulation_time", 10000, "Presimulation Time", "End time for presimulation steady state analysis (double).", "(double) presimulation_maximum_steps takes priority. Only used when allow_presimulation is True");
+    addSetting("allow_approx", true, "Allow Approximiation", "Flag for using steady state approximation routine when steady state solver fails (bool).", "(bool) Approximation routine will run only when the default solver fails to fine a solution. This flag does not affect usage of approximation routine for pre-simulation");
+    addSetting("approx_tolerance", 1e-12, "Approximation Tolerance", "Tolerance for steady state approximation routine (double).", "(double) Absolute tolerance used by steady state approximation routine. Only used when steady state approximation routine is used");
+    addSetting("approx_maximum_steps", 10000, "Approximation Maximum Steps", "Maximum number of steps that can be taken for steady state approximation routine (int).", "(int) Takes priority over approx_time. Only used when steady state approximation routine is used");
+    addSetting("approx_time", 10000, "Approximation Time", "End time for steady state approximation routine (double).", "(double) approx_maximum_steps takes priority. Only used when steady state approximation routine is used");
+    addSetting("relative_tolerance", 1e-16, "Relative Tolerance", "Specifies the relative tolerance (double).", "(double) Relative tolerance used by the solver");
     addSetting("maximum_iterations", 100, "Maximum Iterations", "The maximum number of iterations the solver is allowed to use (int)", "(int) Iteration caps off at the maximum, regardless of whether a solution has been reached");
     addSetting("minimum_damping", 1e-4, "Minimum Damping", "The minimum damping factor (double).", "(double) Minumum damping factor used by the algorithm");
-    addSetting("relative_tolerance", 1e-16, "Relative Tolerance", "Specifies the relative tolerance (double).", "(double) Relative tolerance used by the solver");
+    addSetting("broyden_method", 0, "Broyden Method", "Switches on Broyden method (int)", "(int) Broyden method is a quasi-Newton approximation for rank-1 updates");
+    addSetting("linearity", 3, "Problem Linearity", "Specifies linearity of the problem (int).", "(int) 1 is for linear problem and 4 is for extremly nonlinear problem");
     NLEQSolver::loadConfigSettings();
 }
 
@@ -60,7 +80,7 @@ std::string NLEQSolver::getName() const {
 }
 
 std::string NLEQSolver::getNLEQName() {
-	return "nleq";
+	return "nleq2";
 }
 
 std::string NLEQSolver::getDescription() const {
@@ -68,7 +88,7 @@ std::string NLEQSolver::getDescription() const {
 }
 
 std::string NLEQSolver::getNLEQDescription() {
-	return "NLEQ is a non-linear equation solver which uses a global Newton "
+	return "NLEQ2 is a non-linear equation solver which uses a global Newton "
      "method with adaptive damping strategies (see http://www.zib.de/weimann/NewtonLib/index.html)";
 }
 

source/RK45Integrator.cpp

@@ -21,8 +21,9 @@ namespace rr
 
     RK45Integrator::RK45Integrator(ExecutableModel *m)
     {
-        Log(Logger::LOG_NOTICE) << "Creating Runge-Kutta Fehlberg integrator";
-        stateVectorSize = hCurrent = hmin = hmax = 0;
+        Log(Logger::LOG_NOTICE) << "Creating Runge-Kutta-Fehlberg integrator";
+        resetSettings();
+        stateVectorSize = hmin = hmax = 0;
         k1 = k2 = k3 = k4 = k5 = k6 = err = y = ytmp = NULL;
         syncWithModel(m);
     }
@@ -55,11 +56,10 @@ namespace rr
             err = new double[stateVectorSize];
             y = new double[stateVectorSize];
             ytmp = new double[stateVectorSize];
-            hCurrent = 0.;
             hmin = getValueAsDouble("minimum_time_step"); // TODO: replace hmin with getValueAsDouble("minimum_time_step")
             hmax = getValueAsDouble("maximum_time_step"); // TODO: replace hmax with getValueAsDouble("maximum_time_step")
         } else {
-            stateVectorSize = hCurrent = hmin = hmax = 0;
+            stateVectorSize = hmin = hmax = 0;
             k1 = k2 = k3 = k4 = k5 = k6 = err = y = ytmp = NULL;
         }
     }
@@ -77,19 +77,9 @@ namespace rr
         delete []ytmp;
     }
 
-    double RK45Integrator::integrate(double t, double tEnd)
+    double RK45Integrator::integrate(double t, double tDiff)
     {
-        double h = getValueAsDouble("initial_time_step");
-        if (hCurrent != 0.)
-          h = hCurrent;
-        // NOTE: @kirichoi, please implement rk45 here
-//         static const double epsilon = 1e-3;
-//         double tf = 0;
-        bool singleStep;
-
-        assert(h > hmin && "h must be > hmin");
-//         tf = t + h;
-        singleStep = false;
+        double h = getValueAsDouble("maximum_time_step");
 
         if (!model) {
             throw std::runtime_error("RK45Integrator::integrate: No model");
@@ -98,7 +88,6 @@ namespace rr
         Log(Logger::LOG_DEBUG) <<
                 "RK45Integrator::integrate(" << t << ", " << h << ")";
 
-        // blas daxpy: y -> y + \alpha x
         integer n = stateVectorSize;
         integer inc = 1;
         integer i;
@@ -188,10 +177,31 @@ namespace rr
 
           Log(Logger::LOG_DEBUG) <<
 	    "RK45 step: t = " << t << ", error = " << error << ", epsilon = " << getValueAsDouble("epsilon") << ", h = " << h;
-          if (error <= getValueAsDouble("epsilon")) {
+          if (q <= 0.1) {
+            h = 0.1*h;
+          } else if (q >= 4) {
+            h = 4*h;
+          } else {
+            h = q*h;
+          }
 
-            Log(Logger::LOG_DEBUG) <<
-                  "RK45: Update state vector";
+          if (h > hmax) { 
+              h = hmax; 
+          }
+
+          if (t > tDiff + t) {
+              return tDiff + t;
+          }
+          else if (h > tDiff) {
+              h = tDiff;
+          } else if (h < hmin) {
+              h = hmin;
+              //throw std::runtime_error("RK45Integrator::integrate: Stepsize became smaller than specified minimum.");
+          }
+
+        } while ( error > getValueAsDouble("epsilon"));
+
+        Log(Logger::LOG_DEBUG) << "RK45: Update state vector";
 
         // y = y + (1408*h/2565) k_3
         alpha = 1408.*h / 2565;
@@ -201,7 +211,6 @@ namespace rr
         alpha = (2197. / 4104)*h;
         daxpy_(&n, &alpha, k4, &inc, y, &inc);
 
-
         // y = (y + (1408/2565)*h k_3 + (2197/4104)*h k_4) - (1/5)*h k_5
         alpha = (-1. / 5)*h;
         daxpy_(&n, &alpha, k5, &inc, y, &inc);
@@ -212,32 +221,15 @@ namespace rr
 
         model->setTime(t + h);
         model->setStateVector(y);
-            t = t + h;
 
         for (int i = 0; i<stateVectorSize; ++i) {
             Log(Logger::LOG_DEBUG) << "  " << y[i];
         }
-          }
-
-          if (q <= 0.1) {
-            h = 0.1*h;
-          } else if (q > 4) {
-            h = 4*h;
-          } else {
-            h = q*h;
-          }
-
-          if (h > hmax) { h = hmax; }
-
-        } while ( error > getValueAsDouble("epsilon") && h > hmin );
-
-        hCurrent = h;
-        if (tEnd - t < hCurrent)
-          hCurrent = tEnd - t;
 
         Log(Logger::LOG_DEBUG) <<
                 "RK45: end of step";
-        return t;
+
+        return t + h;
     }
 
     void RK45Integrator::testRootsAtInitialTime()
@@ -285,19 +277,6 @@ namespace rr
         return IntegratorListenerPtr();
     }
 
-    std::string RK45Integrator::toString() const
-    {
-        return toRepr();
-    }
-
-    std::string RK45Integrator::toRepr() const
-    {
-        std::stringstream ss;
-        ss << "< roadrunner.RK45Integrator() { 'this' : "
-                << (void*)this << " }>";
-        return ss.str();
-    }
-
     std::string RK45Integrator::getName() const {
         return RK45Integrator::getRK45Name();
     }
@@ -311,7 +290,9 @@ namespace rr
     }
 
     std::string RK45Integrator::getRK45Description() {
-        return "Kiri, please fill in";
+        return "Runge-Kutta-Fehlberg methods are a family of algorithms for solving "
+			"ODEs. This integrator is based on explicit Runge-Kutta (4,5) formula, automatically "
+			"determining adaptive stepsize.";
     }
 
     std::string RK45Integrator::getHint() const {
@@ -339,11 +320,10 @@ namespace rr
     {
         Solver::resetSettings();
 
-        addSetting("variable_step_size", false, "Variable Step Size", "Perform a variable time step simulation. (bool)", "(bool) Enabling this setting will allow the integrator to adapt the size of each time step. This will result in a non-uniform time column.");
-        addSetting("initial_time_step",  0.5, "Initial Time Step", "Specifies the initial time step size. (double)", "(double) Specifies the initial time step size. If inappropriate, CVODE will attempt to estimate a better initial time step.");
+        addSetting("variable_step_size", true, "Variable Step Size", "Perform a variable time step simulation. (bool)", "(bool) Enabling this setting will allow the integrator to adapt the size of each time step. This will result in a non-uniform time column.");
         addSetting("minimum_time_step",  1e-12, "Minimum Time Step", "Specifies the minimum absolute value of step size allowed. (double)", "(double) The minimum absolute value of step size allowed.");
         addSetting("maximum_time_step",  1.0, "Maximum Time Step", "Specifies the maximum absolute value of step size allowed. (double)", "(double) The maximum absolute value of step size allowed.");
-        addSetting("epsilon",  1e-12, "Maximum error", "TODO: fill in. (double)", "(double) TODO: fill in.");
+        addSetting("epsilon",  1e-12, "Maximum error tolerance", "Specifies the maximum error tolerance allowed. (double)", "(double) The maximum error tolerance allowed.");
     }
 
 } /* namespace rr */

source/RK45Integrator.h

@@ -73,16 +73,6 @@ namespace rr
          */
         virtual void restart(double t0);
 
-        /**
-         * @brief Get a description of this object, compatable with python __str__
-         */
-        virtual std::string toString() const;
-
-        /**
-         * @brief Get a short descriptions of this object, compatable with python __repr__.
-         */
-        virtual std::string toRepr() const;
-
         // ** Meta Info ********************************************************
 
         /**