Moved input-reading outside of solver

Author: CodeByDrescher

IDAWin/CMakeLists.txt

@@ -1,18 +1,25 @@
 project(IDAWin)
 
 set (SRC_FILES 
-	OdeResultSet.cpp
-	StoppedByUserException.cpp	
-	VCellCVodeSolver.cpp
-	VCellIDASolver.cpp
-	VCellSundialsSolver.cpp
+		OdeResultSet.cpp
+		StoppedByUserException.cpp
+		VCellCVodeSolver.cpp
+		VCellIDASolver.cpp
+		VCellSundialsSolver.cpp
+		SteadyStateDriver.cpp
+		VCellSolverFactory.cpp
+		VCellSolverInput.cpp
 )
 set (HEADER_FILES 
-	OdeResultSet.h
-	StoppedByUserException.h	
-	VCellCVodeSolver.h
-	VCellIDASolver.h
-	VCellSundialsSolver.h
+		OdeResultSet.h
+		StoppedByUserException.h
+		VCellSolverFactory.h
+		VCellSolver.h
+		VCellSolverInput.h
+		VCellSundialsSolver.h
+		VCellCVodeSolver.h
+		VCellIDASolver.h
+        SteadyStateDriver.h
 )
 
 set (EXE_SRC_FILES

IDAWin/SteadyStateDriver.cpp

@@ -0,0 +1,4 @@
+//
+// Created by Logan Drescher on 6/25/25.
+//
+#include "SteadyStateDriver.h"

IDAWin/SteadyStateDriver.h

@@ -0,0 +1,8 @@
+//
+// Created by Logan Drescher on 6/25/25.
+//
+
+#ifndef STEADYSTATEDRIVER_H
+#define STEADYSTATEDRIVER_H
+
+#endif //STEADYSTATEDRIVER_H

IDAWin/SundialsSolverStandalone.cpp

@@ -13,6 +13,8 @@
 #include "StoppedByUserException.h"
 #include <VCELL/GitDescribe.h>
 #include <argparse/argparse.hpp>
+
+#include "VCellSolverFactory.h"
 #define CVODE_SOLVER "CVODE"
 #define IDA_SOLVER "IDA"
 
@@ -86,90 +88,47 @@ int parseAndRunWithArgParse(int argc, char *argv[]) {
 }
 
 void activateSolver(std::ifstream& inputFileStream, FILE* outputFile, int taskID) {
-	std::string solver;
-
-	while (!inputFileStream.eof()) { // Note break statement if "SOLVER" encountered
-		std::string nextToken;
-		inputFileStream >> nextToken;
-		if (nextToken.empty()) continue;
-		if (nextToken[0] == '#') getline(inputFileStream, nextToken);
-		else if (nextToken == "JMS_PARAM_BEGIN") {
-			loadJMSInfo(inputFileStream, taskID);
-			#ifdef USE_MESSAGING
-			SimulationMessaging::getInstVar()->start(); // start the thread
-			#endif
-		} else if (nextToken == "SOLVER") {
-			inputFileStream >> solver;
-			break;
-		}
-	}
-	#ifdef USE_MESSAGING
-	// should only happen during testing for solver compiled with messaging but run locally.
-	if (SimulationMessaging::getInstVar() == nullptr) { SimulationMessaging::create(); }
-	#endif
-
-	if (solver.empty()) { throw "Solver not defined "; }
-	VCellSundialsSolver *vss = nullptr;
-	if (solver == IDA_SOLVER) {
-		vss = new VCellIDASolver();
-	} else if (solver == CVODE_SOLVER) {
-		vss = new VCellCVodeSolver();
-	} else {
-		std::stringstream ss;
-		ss << "Solver " << solver << " not defined!";
-		throw ss.str();
-	}
-	vss->readInput(inputFileStream);
-	vss->solve(nullptr, true, outputFile, VCellSundialsSolver::checkStopRequested);
-
-	delete vss;
-}
-
-void loadJMSInfo(std::istream &ifsInput, int taskID) {
-	#ifndef USE_MESSAGING
-		return; // Only useful for messaging; let's not waste time!
-	#else
-
-	if (taskID < 0) {
-		SimulationMessaging::create();
-		return; // No need to do any parsing
-	}
-	std::string broker;
-	std::string smqUserName;
-	std::string password;
-	std::string qName;
-	std::string topicName;
-	std::string vCellUsername;
-	int simKey, jobIndex;
-
-	while (!ifsInput.eof()) {
-		std::string nextToken;
-		ifsInput >> nextToken;
-		if (nextToken.empty()) continue;
-		if (nextToken[0] == '#') {
-			// getline(ifsInput, nextToken); // Is this ignoring because of a comment?
-			ifsInput.ignore('\n');
-			continue;
-		}
-		if (nextToken == "JMS_PARAM_END") { ifsInput.ignore(EOF); } else if (
-			nextToken == "JMS_BROKER") { ifsInput >> broker; } else if (
-			nextToken == "JMS_USER") { ifsInput >> smqUserName >> password; } else if (
-			nextToken == "JMS_QUEUE") { ifsInput >> qName; } else if (
-			nextToken == "JMS_TOPIC") { ifsInput >> topicName; } else if (nextToken == "VCELL_USER") {
-			ifsInput >> vCellUsername;
-		} else if (nextToken == "SIMULATION_KEY") {
-			ifsInput >> simKey;
-			continue;
-		} else if (nextToken == "JOB_INDEX") {
-			ifsInput >> jobIndex;
-			continue;
-		}
-	}
-
-	SimulationMessaging::create(broker.c_str(), smqUserName.c_str(),
-	                            password.c_str(), qName.c_str(), topicName.c_str(),
-	                            vCellUsername.c_str(), simKey, jobIndex, taskID);
-	#endif
+	// std::string solver;
+	//
+	// while (!inputFileStream.eof()) { // Note break statement if "SOLVER" encountered
+	// 	std::string nextToken;
+	// 	inputFileStream >> nextToken;
+	// 	if (nextToken.empty()) continue;
+	// 	if (nextToken[0] == '#') getline(inputFileStream, nextToken);
+	// 	else if (nextToken == "JMS_PARAM_BEGIN") {
+	// 		loadJMSInfo(inputFileStream, taskID);
+	// 		#ifdef USE_MESSAGING
+	// 		SimulationMessaging::getInstVar()->start(); // start the thread
+	// 		#endif
+	// 	} else if (nextToken == "SOLVER") {
+	// 		inputFileStream >> solver;
+	// 		break;
+	// 	}
+	// }
+	// #ifdef USE_MESSAGING
+	// // should only happen during testing for solver compiled with messaging but run locally.
+	// if (SimulationMessaging::getInstVar() == nullptr) { SimulationMessaging::create(); }
+	// #endif
+	//
+	// if (solver.empty()) { throw "Solver not defined "; }
+	// VCellSundialsSolver *vss = nullptr;
+	//
+	// if (solver == IDA_SOLVER) {
+	// 	vss = new VCellIDASolver();
+	// } else if (solver == CVODE_SOLVER) {
+	// 	vss = new VCellCVodeSolver();
+	// } else {
+	// 	std::stringstream ss;
+	// 	ss << "Solver " << solver << " not defined!";
+	// 	throw ss.str();
+	// }
+	//
+	// vss->solve(nullptr, true, outputFile, VCellSundialsSolver::checkStopRequested);
+	// delete vss;
+
+
+	VCellSolver* targetSolver = VCellSolverFactory::produceVCellSolver(inputFileStream, taskID);
+	targetSolver->solve(nullptr, true, outputFile, VCellSundialsSolver::checkStopRequested);
 }
 
 void errExit(int returnCode, std::string &errorMsg) {

IDAWin/VCellCVodeSolver.cpp

@@ -9,20 +9,16 @@
 #include <sys/timeb.h>
 #include <sstream>
 #include <string>
-
-#ifdef USE_MESSAGING
-#include <VCELL/SimulationMessaging.h>
-#endif
-
 #include <memory.h>
-
 #include <cvode/cvode.h>             /* prototypes for CVODE fcts. and consts. */
 #include <nvector/nvector_serial.h>  /* serial N_Vector types, fcts., and macros */
 #include <cvode/cvode_dense.h>       /* prototype for CVDense */
 //#include <cvode/cvode_spgmr.h>       /* prototype for CVSPGMR */
 #include <sundials/sundials_dense.h> /* definitions DenseMat DENSE_ELEM */
 #include <sundials/sundials_types.h> /* definition of type realtype */
-
+#ifdef USE_MESSAGING
+#include <VCELL/SimulationMessaging.h>
+#endif
 #define ToleranceType CV_SS
 
 /**
@@ -151,13 +147,13 @@ VCellCVodeSolver::VCellCVodeSolver() : VCellSundialsSolver() {
 VCellCVodeSolver::~VCellCVodeSolver() {
 	CVodeFree(&solver);
 
-	for (int i = 0; i < NEQ; i ++) {
+	for (int i = 0; i < NUM_EQUATIONS; i ++) {
 		delete rateExpressions[i];
 	}
 	delete[] rateExpressions;
 }
 
-/*
+/**----------------------------------------------------
 Input format:
 	STARTING_TIME 0.0
 	ENDING_TIME 0.1
@@ -172,15 +168,15 @@ Input format:
 		 RATE ((20.0 * x_o * D_B0) - (50.0 * x_i));
 	ODE x_o INIT 0.0;
 		 RATE ( - ((20.0 * x_o * D_B0) - (50.0 * x_i)) + (1505000.0 * (3.322259136212625E-4 - (3.322259136212625E-4 * x_o) - (3.322259136212625E-4 * x_i))) - (100.0 * x_o));
-*/
+--------------------------------------------------------------*/
 void VCellCVodeSolver::readEquations(std::istream& inputstream) {
 	try {
 		std::string name;
 		std::string exp;
 
-		rateExpressions = new Expression*[NEQ];		
+		rateExpressions = new Expression*[NUM_EQUATIONS];
 
-		for (int i = 0; i < NEQ; i ++) {
+		for (int i = 0; i < NUM_EQUATIONS; i ++) {
 			// ODE
 			inputstream >> name >> variableNames[i];			
 
@@ -209,11 +205,35 @@ void VCellCVodeSolver::readEquations(std::istream& inputstream) {
 	}
 }
 
+void VCellCVodeSolver::readEquations(VCellSolverInputBreakdown& inputBreakdown) {
+	try {
+		for (int i = 0; i < NUM_EQUATIONS; i ++) {
+			variableNames[i] = inputBreakdown.modelSettings.VARIABLE_NAMES[i];
+			try {
+				initialConditionExpressions[i] = new Expression(inputBreakdown.modelSettings.INITIAL_CONDITION_EXPRESSIONS[i]);
+			} catch (VCell::Exception& ex) {
+				throw VCell::Exception(std::string("Initial condition expression for [") + variableNames[i] + "] " + ex.getMessage());
+			}
+		}
+
+		rateExpressions = new Expression*[inputBreakdown.modelSettings.RATE_EXPRESSIONS.size()];
+		for (int i = 0; i < inputBreakdown.modelSettings.RATE_EXPRESSIONS.size(); i++) {
+			rateExpressions[i] = new Expression(inputBreakdown.modelSettings.RATE_EXPRESSIONS[i]);
+		}
+	} catch (char* ex) {
+		throw VCell::Exception(std::string("VCellCVodeSolver::readInput() : ") + ex);
+	} catch (VCell::Exception& ex) {
+		throw VCell::Exception(std::string("VCellCVodeSolver::readInput() : ") + ex.getMessage());
+	} catch (...) {
+		throw "VCellCVodeSolver::readInput() : caught unknown exception";
+	}
+}
+
 void VCellCVodeSolver::initialize() {
 	VCellSundialsSolver::initialize();
 
 	// rate can be function of variables, parameters and discontinuities.
-	for (int i = 0; i < NEQ; i ++) {
+	for (int i = 0; i < NUM_EQUATIONS; i ++) {
 		rateExpressions[i]->bindExpression(defaultSymbolTable);
 	}
 }
@@ -222,7 +242,7 @@ int VCellCVodeSolver::RHS (realtype t, N_Vector y, N_Vector r) {
 	try {
 		updateTandVariableValues(t, y);
 		double* r_data = NV_DATA_S(r);
-		for (int i = 0; i < NEQ; i ++) {
+		for (int i = 0; i < NUM_EQUATIONS; i ++) {
 			r_data[i] = rateExpressions[i]->evaluateVector(values);
 		}
 		recoverableErrMsg = "";
@@ -264,17 +284,17 @@ void VCellCVodeSolver::solve(double* paramValues, bool bPrintProgress, FILE* out
 	this->odeResultSet->clearData();
 
 	// copy parameter values to the end of values, these will stay the same during solving
-	memset(values, 0, (NEQ + 1) * sizeof(double));
-	memcpy(values + 1 + NEQ, paramValues, NPARAM * sizeof(double));
-	memset(values + 1 + NEQ + NPARAM, 0, numDiscontinuities * sizeof(double));
+	memset(values, 0, (NUM_EQUATIONS + 1) * sizeof(double));
+	memcpy(values + 1 + NUM_EQUATIONS, paramValues, NUM_PARAMETERS * sizeof(double));
+	memset(values + 1 + NUM_EQUATIONS + NUM_PARAMETERS, 0, numDiscontinuities * sizeof(double));
 
 	initCVode(paramValues);
 	cvodeSolve(bPrintProgress, outputFile, checkStopRequested);
 }
 
 void VCellCVodeSolver::initCVode(double* paramValues) {
 	//Initialize y, variable portion of values
-	for (int i = 0; i < NEQ; i ++) {
+	for (int i = 0; i < NUM_EQUATIONS; i ++) {
 		NV_Ith_S(y, i) = initialConditionExpressions[i]->evaluateVector(paramValues);
 	}
 
@@ -305,7 +325,7 @@ void VCellCVodeSolver::reInit(double t) {
 
 		flag = CVodeSetFdata(solver, this);
 		checkCVodeFlag(flag);
-		flag = CVDense(solver, NEQ);
+		flag = CVDense(solver, NUM_EQUATIONS);
 		//flag = CVSpgmr(solver, PREC_NONE, 0);
 		checkCVodeFlag(flag);
 
@@ -317,8 +337,8 @@ void VCellCVodeSolver::reInit(double t) {
 }
 
 bool VCellCVodeSolver::fixInitialDiscontinuities(double t) {
-	double* oldy = new double[NEQ];
-	memcpy(oldy, NV_DATA_S(y), NEQ * sizeof(realtype));
+	double* oldy = new double[NUM_EQUATIONS];
+	memcpy(oldy, NV_DATA_S(y), NUM_EQUATIONS * sizeof(realtype));
 
 	double epsilon = std::max(1e-15, ENDING_TIME * 1e-10);
 	double currentTime = t;	
@@ -341,11 +361,11 @@ bool VCellCVodeSolver::fixInitialDiscontinuities(double t) {
 		}
 	}	
 	if (bInitChanged) {
-		memcpy(values + 1 + NEQ + NPARAM, discontinuityValues, numDiscontinuities * sizeof(double));
+		memcpy(values + 1 + NUM_EQUATIONS + NUM_PARAMETERS, discontinuityValues, numDiscontinuities * sizeof(double));
 	}
 
 	//revert y
-	memcpy(NV_DATA_S(y), oldy, NEQ * sizeof(realtype));
+	memcpy(NV_DATA_S(y), oldy, NUM_EQUATIONS * sizeof(realtype));
 	reInit(t);
 
 	delete[] oldy;
@@ -427,7 +447,7 @@ void VCellCVodeSolver::cvodeSolve(bool bPrintProgress, FILE* outputFile, void (*
 
 				if (returnCode == CV_ROOT_RETURN || iterationCount % keepEvery == 0 || Time >= ENDING_TIME){
 					outputCount++;
-					if (outputCount * (NEQ + 1) * bytesPerSample > MaxFileSizeBytes){ 
+					if (outputCount * (NUM_EQUATIONS + 1) * bytesPerSample > MaxFileSizeBytes){
 						/* if more than one gigabyte, then fail */
 						const std::string msg{"output exceeded maximum " + std::to_string(MaxFileSizeBytes) + " bytes"};
 						throw VCell::Exception(msg.c_str());
@@ -487,5 +507,5 @@ void VCellCVodeSolver::cvodeSolve(bool bPrintProgress, FILE* outputFile, void (*
 
 void VCellCVodeSolver::updateTandVariableValues(realtype t, N_Vector y) {
 	values[0] = t;
-	memcpy(values + 1, NV_DATA_S(y), NEQ * sizeof(realtype));
+	memcpy(values + 1, NV_DATA_S(y), NUM_EQUATIONS * sizeof(realtype));
 }

IDAWin/VCellCVodeSolver.h

@@ -14,8 +14,10 @@ class VCellCVodeSolver : public VCellSundialsSolver {
 
 protected:
 	void readEquations(std::istream& inputstream) override;
+	void readEquations(VCellSolverInputBreakdown& inputBreakdown) override;
 	void initialize() override;
 	std::string getSolverName() override { return "CVODE"; }
+	VCellSolverTypes getSolverType() override { return VCellSolverTypes::CVODE; }
 
 private:
 	VCell::Expression** rateExpressions;