diff --git a/src/Matrix/FspMatrixBase.cpp b/src/Matrix/FspMatrixBase.cpp index 0dd8dba..3765b9a 100644 --- a/src/Matrix/FspMatrixBase.cpp +++ b/src/Matrix/FspMatrixBase.cpp @@ -38,6 +38,8 @@ int FspMatrixBase::Action(PetscReal t, Vec x, Vec y) { ierr = VecSet(y, 0.0); CHKERRQ(ierr); + if (has_values_ == PETSC_FALSE) return 0; + if (!tv_reactions_.empty()) { ierr = t_fun_(t, num_reactions_, time_coefficients_.memptr(), t_fun_args_); PACMENSLCHKERRQ(ierr); @@ -244,6 +246,7 @@ PacmenslErrorCode FspMatrixBase::GenerateValues(const StateSetBase &fsp, ierr = MatAssemblyEnd(ti_mat_, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); } + has_values_ = PETSC_TRUE; return 0; } @@ -266,7 +269,7 @@ int FspMatrixBase::Destroy() { ierr = VecDestroy(work_.mem()); CHKERRQ(ierr); } - + has_values_ = PETSC_FALSE; return 0; } diff --git a/src/Matrix/FspMatrixBase.h b/src/Matrix/FspMatrixBase.h index 84af916..ac23b24 100644 --- a/src/Matrix/FspMatrixBase.h +++ b/src/Matrix/FspMatrixBase.h @@ -58,7 +58,7 @@ class FspMatrixBase { * @param comm Communicator context for the processes that own the matrix. */ explicit FspMatrixBase(MPI_Comm comm); - + /** * @brief Populate the fields of the FSP matrix object. * This method is __collective__, meaning that it must be called by all processes that own this object. @@ -69,13 +69,13 @@ class FspMatrixBase { virtual PacmenslErrorCode GenerateValues(const StateSetBase &fsp, const Model &model); - + /** * @brief Populate the fields of the FSP matrix object. * This method is __collective__, meaning that it must be called by all processes that own this object. * @param fsp (in) The FSP-truncated state space from which the transition rate matrix is built. * @param SM (in) The stoichiometry matrix of the reaction network. - * @param time_varying (in) List of reactions whose propensities are time-varying. + * @param time_varying (in) List of reactions whose propensities are time-varying. If empty, all reaction propensities are time-independent. * @param new_prop_t (in) Function pointer to evaluate the vector of time-varying coefficients $c_r(t,\theta)$ at time $t$. * @param new_prop_x (in) Function pointer to evaluate the time-independent component of the propensities. * @param enable_reactions (in) Vector of reactions that are active. If left empty, we assume all reactions are active. Propensity functions of inactive reactions are not added to the data structure. @@ -92,11 +92,15 @@ class FspMatrixBase { const std::vector &enable_reactions, void *prop_t_args, void *prop_x_args); - + PacmenslErrorCode SetTimeFun(TcoefFun new_t_fun, void *new_t_fun_args); - + + /** + * @brief Clear internal state of the object. + * @return Error code: 0 (success) or -1 (failure). + */ virtual int Destroy(); - + /** * @brief Compute y = A*x. * This method is __collective__, meaning that it must be called by all processes that own this object. @@ -106,7 +110,7 @@ class FspMatrixBase { * @return Error code: 0 (success) or -1 (failure) */ virtual PacmenslErrorCode Action(PetscReal t, Vec x, Vec y); - + /** * @brief Generate a PETSc Matrix object with the same sparsity structure as the FSP matrix. This serves as the Jacobian matrix required by ODE solvers. * This method is __collective__, meaning that it must be called by all processes that own this object. @@ -114,7 +118,7 @@ class FspMatrixBase { * @return Error code. 0 (success) or -1 (failure) */ virtual PacmenslErrorCode CreateRHSJacobian(Mat *A); - + /** * @brief Populate the values of a PETSc Mat object with the entries of the FSP matrix at any given time. * This method is __collective__, meaning that it must be called by all processes that own this object. @@ -124,7 +128,7 @@ class FspMatrixBase { * @return Error code: 0 (success) or -1 (failure) */ virtual PacmenslErrorCode ComputeRHSJacobian(PetscReal t, Mat A); - + /** * @brief Get an estimate for the number of floating-point operations needed by the calling process per matrix-vector multiplication (i.e., per call to the \ref Action() method). * This method is __local__, meaning that it can be called independently by the local process. @@ -132,14 +136,14 @@ class FspMatrixBase { * @return Error code: 0 (success) or -1 (failure) */ virtual PacmenslErrorCode GetLocalMVFlops(PetscInt *nflops); - + /** * @brief Get the number of rows owned by the calling process. * This method is __local__, meaning that it can be called independently by the local process. * @return Number of rows owned by the calling process. */ int GetNumLocalRows() const { return num_rows_local_; }; - + /** * @brief Object destructor. */ @@ -148,26 +152,26 @@ class FspMatrixBase { MPI_Comm comm_ = MPI_COMM_NULL; ///< Communicator context int rank_; ///< Rank of local process int comm_size_; ///< Total number of processes that own this object - + Int num_reactions_ = 0; ///< Number of reactions Int num_rows_global_ = 0; ///< Global number of rows Int num_rows_local_ = 0; ///< Number of matrix rows owned by this process - + std::vector enable_reactions_ = std::vector(); std::vector tv_reactions_ = std::vector(); std::vector ti_reactions_ = std::vector(); - + // Local data of the matrix std::vector> tv_mats_; ///< List of PETSc Mat objects for storing the $A_r$ factors for reactions $r$ with time-varying propensities Petsc ti_mat_; ///< PETsc Mat object to store the time-independent term of the matrix sum decomposition - + // Data for computing the matrix action Petsc work_; ///< Work vector for computing operator times vector - + TcoefFun t_fun_ = nullptr; ///< Pointer to external function that evaluates time-varying factors of reaction propensities. void *t_fun_args_ = nullptr; ///< Pointer for extra data needed in \ref t_fun_ evaluation. arma::Row time_coefficients_; ///< Vector to store the time-varying factors of reaction propensities - + /** * @brief Populate the data members with information about the inter-process layout of the object. * This method is __collective__, meaning that it must be called by all processes that own this object. @@ -176,7 +180,7 @@ class FspMatrixBase { * Modified fields: \ref num_rows_local_, \ref num_rows_global_, \ref work_. */ virtual int DetermineLayout_(const StateSetBase &fsp); - + // arrays for counting nonzero entries on the diagonal and off-diagonal blocks arma::Mat dblock_nz_, oblock_nz_; arma::Col ti_dblock_nz_, ti_oblock_nz_; @@ -185,6 +189,8 @@ class FspMatrixBase { // array of matrix values arma::Mat offdiag_vals_; arma::Mat diag_vals_; + + PetscBool has_values_ = PETSC_FALSE; ///< whether the matrix has values. }; } diff --git a/src/Matrix/FspMatrixConstrained.cpp b/src/Matrix/FspMatrixConstrained.cpp index 2ca5e41..41c98d3 100644 --- a/src/Matrix/FspMatrixConstrained.cpp +++ b/src/Matrix/FspMatrixConstrained.cpp @@ -25,111 +25,97 @@ SOFTWARE. #include "FspMatrixConstrained.h" namespace pacmensl { -FspMatrixConstrained::FspMatrixConstrained(MPI_Comm comm) : FspMatrixBase(comm) -{ +FspMatrixConstrained::FspMatrixConstrained(MPI_Comm comm) : FspMatrixBase(comm) { } +int FspMatrixConstrained::Action(PetscReal t, Vec x, Vec y) { + int ierr; -int FspMatrixConstrained::Action(PetscReal t,Vec x,Vec y) -{ - int ierr; - - ierr = FspMatrixBase::Action(t,x,y); - PACMENSLCHKERRQ(ierr); - ierr = VecGetLocalVectorRead(x,xx); - CHKERRQ(ierr); - ierr = VecSet(sink_entries_,0.0); - CHKERRQ(ierr); - if (!tv_reactions_.empty()) - { - for (int i {0}; i < tv_reactions_.size(); ++i) - { - ierr = MatMult(tv_sinks_mat_[i],xx,sink_tmp); - CHKERRQ(ierr); - ierr = VecAXPY(sink_entries_,time_coefficients_[tv_reactions_[i]],sink_tmp); - CHKERRQ(ierr); - } - } + ierr = FspMatrixBase::Action(t, x, y); + PACMENSLCHKERRQ(ierr); + if (has_values_ == PETSC_FALSE) return 0; - if (ti_sinks_mat_ != nullptr) - { - ierr = MatMult(ti_sinks_mat_,xx,sink_tmp); + ierr = VecGetLocalVectorRead(x, xx); CHKERRQ(ierr); - ierr = VecAXPY(sink_entries_,1.0,sink_tmp); + ierr = VecSet(sink_entries_, 0.0); CHKERRQ(ierr); - } - ierr = VecScatterBegin(sink_scatter_ctx_,sink_entries_,y,ADD_VALUES,SCATTER_FORWARD); - CHKERRQ(ierr); - ierr = VecScatterEnd(sink_scatter_ctx_,sink_entries_,y,ADD_VALUES,SCATTER_FORWARD); - CHKERRQ(ierr); - ierr = VecRestoreLocalVectorRead(x,xx); - CHKERRQ(ierr); - return 0; -} + if (!tv_reactions_.empty()) { + for (int i{0}; i < tv_reactions_.size(); ++i) { + ierr = MatMult(tv_sinks_mat_[i], xx, sink_tmp); + CHKERRQ(ierr); + ierr = VecAXPY(sink_entries_, time_coefficients_[tv_reactions_[i]], sink_tmp); + CHKERRQ(ierr); + } + } -int FspMatrixConstrained::Destroy() -{ - PetscErrorCode ierr; - FspMatrixBase::Destroy(); - if (sink_entries_ != nullptr) - { - ierr = VecDestroy(&sink_entries_); - CHKERRQ(ierr); - } - if (sink_tmp != nullptr) - { - ierr = VecDestroy(&sink_tmp); - CHKERRQ(ierr); - } - if (!tv_sinks_mat_.empty()) - { - for (int i{0}; i < tv_sinks_mat_.size(); ++i) - { - if (tv_sinks_mat_[i]) - { - ierr = MatDestroy(&tv_sinks_mat_[i]); + if (ti_sinks_mat_ != nullptr) { + ierr = MatMult(ti_sinks_mat_, xx, sink_tmp); + CHKERRQ(ierr); + ierr = VecAXPY(sink_entries_, 1.0, sink_tmp); CHKERRQ(ierr); - } } - tv_sinks_mat_.clear(); - } - if (xx != nullptr) - { - ierr = VecDestroy(&xx); + ierr = VecScatterBegin(sink_scatter_ctx_, sink_entries_, y, ADD_VALUES, SCATTER_FORWARD); CHKERRQ(ierr); - } - if (ti_sinks_mat_ != nullptr) - { - ierr = MatDestroy(&ti_sinks_mat_); + ierr = VecScatterEnd(sink_scatter_ctx_, sink_entries_, y, ADD_VALUES, SCATTER_FORWARD); CHKERRQ(ierr); - } - if (sink_scatter_ctx_ != nullptr) - { - ierr = VecScatterDestroy(&sink_scatter_ctx_); + ierr = VecRestoreLocalVectorRead(x, xx); CHKERRQ(ierr); - } - sinkmat_nnz.clear(); - sinkmat_inz.clear(); - sinkmat_entries.clear(); + return 0; +} + +int FspMatrixConstrained::Destroy() { + PetscErrorCode ierr; + FspMatrixBase::Destroy(); + if (sink_entries_ != nullptr) { + ierr = VecDestroy(&sink_entries_); + CHKERRQ(ierr); + } + if (sink_tmp != nullptr) { + ierr = VecDestroy(&sink_tmp); + CHKERRQ(ierr); + } + if (!tv_sinks_mat_.empty()) { + for (int i{0}; i < tv_sinks_mat_.size(); ++i) { + if (tv_sinks_mat_[i]) { + ierr = MatDestroy(&tv_sinks_mat_[i]); + CHKERRQ(ierr); + } + } + tv_sinks_mat_.clear(); + } + if (xx != nullptr) { + ierr = VecDestroy(&xx); + CHKERRQ(ierr); + } + if (ti_sinks_mat_ != nullptr) { + ierr = MatDestroy(&ti_sinks_mat_); + CHKERRQ(ierr); + } + if (sink_scatter_ctx_ != nullptr) { + ierr = VecScatterDestroy(&sink_scatter_ctx_); + CHKERRQ(ierr); + } + sinkmat_nnz.clear(); + sinkmat_inz.clear(); + sinkmat_entries.clear(); - return 0; + return 0; } -FspMatrixConstrained::~FspMatrixConstrained() -{ - Destroy(); +FspMatrixConstrained::~FspMatrixConstrained() { + Destroy(); } PacmenslErrorCode FspMatrixConstrained::GenerateValues(const StateSetBase &fsp, const Model &model) { - return FspMatrixConstrained::GenerateValues(fsp, - model.stoichiometry_matrix_, - model.tv_reactions_, - model.prop_t_, - model.prop_x_, - std::vector(), - model.prop_t_args_, - model.prop_x_args_ - ); + return FspMatrixConstrained::GenerateValues(fsp, + model.stoichiometry_matrix_, + model.tv_reactions_, + model.prop_t_, + model.prop_x_, + std::vector(), + model.prop_t_args_, + model.prop_x_args_ + ); } PacmenslErrorCode FspMatrixConstrained::GenerateValues(const StateSetBase &state_set, @@ -139,383 +125,343 @@ PacmenslErrorCode FspMatrixConstrained::GenerateValues(const StateSetBase &state const PropFun &prop, const std::vector &enable_reactions, void *prop_t_args, - void *prop_args) -{ - PetscErrorCode ierr{0}; - - PetscInt own_start,own_end; - - auto *constrained_fss_ptr = dynamic_cast(&state_set); - if (!constrained_fss_ptr) ierr = -1; - PACMENSLCHKERRQ(ierr); - - sinks_rank_ = comm_size_ - 1; // rank of the processor that holds sink states_ - try - { - num_constraints_ = constrained_fss_ptr->GetNumConstraints(); - } catch (std::runtime_error &err) - { - ierr = -1; + void *prop_args) { + PetscErrorCode ierr{0}; + + PetscInt own_start, own_end; + + auto *constrained_fss_ptr = dynamic_cast(&state_set); + if (!constrained_fss_ptr) ierr = -1; PACMENSLCHKERRQ(ierr); - } - - // Generate the entries corresponding to usual states_ - ierr = FspMatrixBase::GenerateValues(state_set, - SM, - time_vayring, - new_prop_t, - prop, - enable_reactions, - prop_t_args, - prop_args); - PACMENSLCHKERRQ(ierr); - - // Generate the extra blocks corresponding to sink states_ - const arma::Mat &state_list = constrained_fss_ptr->GetStatesRef(); - int n_local_states = constrained_fss_ptr->GetNumLocalStates(); - int n_constraints = constrained_fss_ptr->GetNumConstraints(); - arma::Mat can_reach_my_state; - - arma::Mat reachable_from_X(state_list.n_rows,state_list.n_cols); - // Workspace for checking constraints - arma::Mat constraints_satisfied(n_local_states,n_constraints); - - sinkmat_nnz.resize(n_constraints,num_reactions_); - sinkmat_inz.resize(n_constraints * num_reactions_); - sinkmat_entries.resize(n_constraints * num_reactions_); - tv_sinks_mat_.resize(tv_reactions_.size()); - for (auto i_reaction : enable_reactions_) - { - // Count nnz for rows that represent sink states_ - can_reach_my_state = state_list + arma::repmat(SM.col(i_reaction),1,state_list.n_cols); - ierr = constrained_fss_ptr->CheckConstraints(n_local_states,can_reach_my_state.colptr(0), - constraints_satisfied.colptr(0)); + + sinks_rank_ = comm_size_ - 1; // rank of the processor that holds sink states_ + try { + num_constraints_ = constrained_fss_ptr->GetNumConstraints(); + } catch (std::runtime_error &err) { + ierr = -1; + PACMENSLCHKERRQ(ierr); + } + + // Generate the entries corresponding to usual states_ + ierr = FspMatrixBase::GenerateValues(state_set, + SM, + time_vayring, + new_prop_t, + prop, + enable_reactions, + prop_t_args, + prop_args); PACMENSLCHKERRQ(ierr); - for (int i_constr = 0; i_constr < n_constraints; ++i_constr) - { - sinkmat_nnz(i_constr,i_reaction) = n_local_states - arma::sum(constraints_satisfied.col(i_constr)); - sinkmat_inz.at(n_constraints * i_reaction + i_constr).set_size(sinkmat_nnz(i_constr,i_reaction)); - sinkmat_entries.at(n_constraints * i_reaction + i_constr).set_size(sinkmat_nnz(i_constr,i_reaction)); + // Generate the extra blocks corresponding to sink states_ + const arma::Mat &state_list = constrained_fss_ptr->GetStatesRef(); + int n_local_states = constrained_fss_ptr->GetNumLocalStates(); + int n_constraints = constrained_fss_ptr->GetNumConstraints(); + arma::Mat can_reach_my_state; + + arma::Mat reachable_from_X(state_list.n_rows, state_list.n_cols); + // Workspace for checking constraints + arma::Mat constraints_satisfied(n_local_states, n_constraints); + + sinkmat_nnz.resize(n_constraints, num_reactions_); + sinkmat_inz.resize(n_constraints * num_reactions_); + sinkmat_entries.resize(n_constraints * num_reactions_); + tv_sinks_mat_.resize(tv_reactions_.size()); + for (auto i_reaction : enable_reactions_) { + // Count nnz for rows that represent sink states_ + can_reach_my_state = state_list + arma::repmat(SM.col(i_reaction), 1, state_list.n_cols); + ierr = constrained_fss_ptr->CheckConstraints(n_local_states, can_reach_my_state.colptr(0), + constraints_satisfied.colptr(0)); + PACMENSLCHKERRQ(ierr); + + for (int i_constr = 0; i_constr < n_constraints; ++i_constr) { + sinkmat_nnz(i_constr, i_reaction) = n_local_states - arma::sum(constraints_satisfied.col(i_constr)); + sinkmat_inz.at(n_constraints * i_reaction + i_constr).set_size(sinkmat_nnz(i_constr, i_reaction)); + sinkmat_entries.at(n_constraints * i_reaction + i_constr).set_size(sinkmat_nnz(i_constr, i_reaction)); + } + // Store the column indices and values of the nonzero entries on the sink rows + for (int i_constr = 0; i_constr < n_constraints; ++i_constr) { + int count = 0; + for (int i_state = 0; i_state < n_local_states; ++i_state) { + if (constraints_satisfied(i_state, i_constr) == 0) { + sinkmat_inz.at(n_constraints * i_reaction + i_constr).at(count) = i_state; + prop(i_reaction, state_list.n_rows, 1, state_list.colptr(i_state), + &sinkmat_entries.at(n_constraints * i_reaction + i_constr)[count], prop_args); + count += 1; + } + } + } + } + + // Fill values for the time varying matrix + for (int i{0}; i < tv_reactions_.size(); ++i) { + int i_reaction = tv_reactions_[i]; + ierr = MatCreate(PETSC_COMM_SELF, &tv_sinks_mat_[i]); + CHKERRQ(ierr); + ierr = MatSetType(tv_sinks_mat_[i], MATSEQAIJ); + CHKERRQ(ierr); + ierr = MatSetSizes(tv_sinks_mat_[i], n_constraints, num_rows_local_, n_constraints, num_rows_local_); + CHKERRQ(ierr); + ierr = MatSeqAIJSetPreallocation(tv_sinks_mat_[i], NULL, sinkmat_nnz.colptr(i_reaction)); + CHKERRQ(ierr); + for (auto i_constr{0}; i_constr < n_constraints; i_constr++) { + ierr = MatSetValues(tv_sinks_mat_[i], 1, &i_constr, sinkmat_nnz(i_constr, i_reaction), + sinkmat_inz.at(i_reaction * n_constraints + i_constr).memptr(), + sinkmat_entries.at(i_reaction * n_constraints + i_constr).memptr(), ADD_VALUES); + CHKERRQ(ierr); + } + ierr = MatAssemblyBegin(tv_sinks_mat_[i], MAT_FINAL_ASSEMBLY); + CHKERRQ(ierr); + ierr = MatAssemblyEnd(tv_sinks_mat_[i], MAT_FINAL_ASSEMBLY); + CHKERRQ(ierr); } - // Store the column indices and values of the nonzero entries on the sink rows - for (int i_constr = 0; i_constr < n_constraints; ++i_constr) - { - int count = 0; - for (int i_state = 0; i_state < n_local_states; ++i_state) - { - if (constraints_satisfied(i_state,i_constr) == 0) - { - sinkmat_inz.at(n_constraints * i_reaction + i_constr).at(count) = i_state; - prop(i_reaction,state_list.n_rows,1,state_list.colptr(i_state), - &sinkmat_entries.at(n_constraints * i_reaction + i_constr)[count],prop_args); - count += 1; + + // Fill values for the time-invariant matrix + if (!ti_reactions_.empty()) { + ierr = MatCreate(PETSC_COMM_SELF, &ti_sinks_mat_); + CHKERRQ(ierr); + ierr = MatSetType(ti_sinks_mat_, MATSELL); + CHKERRQ(ierr); + ierr = MatSetSizes(ti_sinks_mat_, n_constraints, num_rows_local_, n_constraints, num_rows_local_); + CHKERRQ(ierr); + ierr = MatSeqSELLSetPreallocation(ti_sinks_mat_, n_constraints * num_rows_local_, NULL); + CHKERRQ(ierr); + for (auto i_reaction: ti_reactions_) { + for (auto i_constr{0}; i_constr < n_constraints; i_constr++) { + ierr = MatSetValues(ti_sinks_mat_, 1, &i_constr, sinkmat_nnz(i_constr, i_reaction), + sinkmat_inz.at(i_reaction * n_constraints + i_constr).memptr(), + sinkmat_entries.at(i_reaction * n_constraints + i_constr).memptr(), ADD_VALUES); + CHKERRQ(ierr); + } } - } + ierr = MatAssemblyBegin(ti_sinks_mat_, MAT_FINAL_ASSEMBLY); + CHKERRQ(ierr); + ierr = MatAssemblyEnd(ti_sinks_mat_, MAT_FINAL_ASSEMBLY); + CHKERRQ(ierr); } - } - // Fill values for the time varying matrix - for (int i{0}; i < tv_reactions_.size(); ++i) - { - int i_reaction = tv_reactions_[i]; - ierr = MatCreate(PETSC_COMM_SELF,&tv_sinks_mat_[i]); - CHKERRQ(ierr); - ierr = MatSetType(tv_sinks_mat_[i],MATSEQAIJ); - CHKERRQ(ierr); - ierr = MatSetSizes(tv_sinks_mat_[i],n_constraints,num_rows_local_,n_constraints,num_rows_local_); - CHKERRQ(ierr); - ierr = MatSeqAIJSetPreallocation(tv_sinks_mat_[i],NULL,sinkmat_nnz.colptr(i_reaction)); + // Update sinkmat_inz to use global indices + ierr = VecGetOwnershipRange(work_, &own_start, &own_end); CHKERRQ(ierr); - for (auto i_constr{0}; i_constr < n_constraints; i_constr++) - { - ierr = MatSetValues(tv_sinks_mat_[i],1,&i_constr,sinkmat_nnz(i_constr,i_reaction), - sinkmat_inz.at(i_reaction * n_constraints + i_constr).memptr(), - sinkmat_entries.at(i_reaction * n_constraints + i_constr).memptr(),ADD_VALUES); - CHKERRQ(ierr); + for (auto i_reaction:enable_reactions_) { + for (auto i_constr{0}; i_constr < n_constraints; i_constr++) { + sinkmat_inz.at(i_reaction * n_constraints + i_constr) += own_start; + } } - ierr = MatAssemblyBegin(tv_sinks_mat_[i],MAT_FINAL_ASSEMBLY); + + // Local vectors for computing sink entries + ierr = VecCreateSeq(PETSC_COMM_SELF, n_constraints, &sink_entries_); CHKERRQ(ierr); - ierr = MatAssemblyEnd(tv_sinks_mat_[i],MAT_FINAL_ASSEMBLY); + ierr = VecSetUp(sink_entries_); CHKERRQ(ierr); - } - // Fill values for the time-invariant matrix - if (!ti_reactions_.empty()) - { - ierr = MatCreate(PETSC_COMM_SELF,&ti_sinks_mat_); + ierr = VecDuplicate(sink_entries_, &sink_tmp); CHKERRQ(ierr); - ierr = MatSetType(ti_sinks_mat_,MATSELL); + ierr = VecSetUp(sink_tmp); CHKERRQ(ierr); - ierr = MatSetSizes(ti_sinks_mat_,n_constraints,num_rows_local_,n_constraints,num_rows_local_); + + ierr = VecCreateSeq(PETSC_COMM_SELF, num_rows_local_, &xx); CHKERRQ(ierr); - ierr = MatSeqSELLSetPreallocation(ti_sinks_mat_,n_constraints * num_rows_local_,NULL); + ierr = VecSetUp(xx); CHKERRQ(ierr); - for (auto i_reaction: ti_reactions_) - { - for (auto i_constr{0}; i_constr < n_constraints; i_constr++) - { - ierr = MatSetValues(ti_sinks_mat_,1,&i_constr,sinkmat_nnz(i_constr,i_reaction), - sinkmat_inz.at(i_reaction * n_constraints + i_constr).memptr(), - sinkmat_entries.at(i_reaction * n_constraints + i_constr).memptr(),ADD_VALUES); - CHKERRQ(ierr); - } + + // Scatter context for adding sink values + int *sink_global_indices = new int[n_constraints]; + for (int i{0}; i < n_constraints; ++i) { + sink_global_indices[i] = constrained_fss_ptr->GetNumGlobalStates() + i; } - ierr = MatAssemblyBegin(ti_sinks_mat_,MAT_FINAL_ASSEMBLY); + IS sink_is; + ierr = ISCreateGeneral(comm_, n_constraints, sink_global_indices, PETSC_COPY_VALUES, &sink_is); CHKERRQ(ierr); - ierr = MatAssemblyEnd(ti_sinks_mat_,MAT_FINAL_ASSEMBLY); + ierr = VecScatterCreate(sink_entries_, NULL, work_, sink_is, &sink_scatter_ctx_); CHKERRQ(ierr); - } - - // Update sinkmat_inz to use global indices - ierr = VecGetOwnershipRange(work_,&own_start,&own_end); - CHKERRQ(ierr); - for (auto i_reaction:enable_reactions_) - { - for (auto i_constr{0}; i_constr < n_constraints; i_constr++) - { - sinkmat_inz.at(i_reaction * n_constraints + i_constr) += own_start; - } - } - - // Local vectors for computing sink entries - ierr = VecCreateSeq(PETSC_COMM_SELF,n_constraints,&sink_entries_); - CHKERRQ(ierr); - ierr = VecSetUp(sink_entries_); - CHKERRQ(ierr); - - ierr = VecDuplicate(sink_entries_,&sink_tmp); - CHKERRQ(ierr); - ierr = VecSetUp(sink_tmp); - CHKERRQ(ierr); - - ierr = VecCreateSeq(PETSC_COMM_SELF,num_rows_local_,&xx); - CHKERRQ(ierr); - ierr = VecSetUp(xx); - CHKERRQ(ierr); - - // Scatter context for adding sink values - int *sink_global_indices = new int[n_constraints]; - for (int i{0}; i < n_constraints; ++i) - { - sink_global_indices[i] = constrained_fss_ptr->GetNumGlobalStates() + i; - } - IS sink_is; - ierr = ISCreateGeneral(comm_,n_constraints,sink_global_indices,PETSC_COPY_VALUES,&sink_is); - CHKERRQ(ierr); - ierr = VecScatterCreate(sink_entries_,NULL,work_,sink_is,&sink_scatter_ctx_); - CHKERRQ(ierr); - ierr = ISDestroy(&sink_is); - CHKERRQ(ierr); - delete[] sink_global_indices; - return 0; + ierr = ISDestroy(&sink_is); + CHKERRQ(ierr); + delete[] sink_global_indices; + return 0; } -PacmenslErrorCode FspMatrixConstrained::DetermineLayout_(const StateSetBase &fsp) -{ - PetscErrorCode ierr; - - num_rows_local_ = fsp.GetNumLocalStates(); - if (rank_ == sinks_rank_) num_rows_local_ += num_constraints_; - - // Generate matrix layout from Fsp's layout - ierr = VecCreate(comm_,work_.mem()); - CHKERRQ(ierr); - ierr = VecSetFromOptions(work_); - CHKERRQ(ierr); - ierr = VecSetSizes(work_,num_rows_local_,PETSC_DECIDE); - CHKERRQ(ierr); - ierr = VecSetUp(work_); - CHKERRQ(ierr); - ierr = VecGetSize(work_,&num_rows_global_); - CHKERRQ(ierr); - return 0; +PacmenslErrorCode FspMatrixConstrained::DetermineLayout_(const StateSetBase &fsp) { + PetscErrorCode ierr; + + num_rows_local_ = fsp.GetNumLocalStates(); + if (rank_ == sinks_rank_) num_rows_local_ += num_constraints_; + + // Generate matrix layout from Fsp's layout + ierr = VecCreate(comm_, work_.mem()); + CHKERRQ(ierr); + ierr = VecSetFromOptions(work_); + CHKERRQ(ierr); + ierr = VecSetSizes(work_, num_rows_local_, PETSC_DECIDE); + CHKERRQ(ierr); + ierr = VecSetUp(work_); + CHKERRQ(ierr); + ierr = VecGetSize(work_, &num_rows_global_); + CHKERRQ(ierr); + return 0; } -int FspMatrixConstrained::CreateRHSJacobian(Mat *A) -{ - int ierr; - - PetscInt own_start,own_end,num_local_states,itmp; - arma::Col d_nz,o_nz,tmp; - - ierr = MatCreate(comm_,A); - CHKERRQ(ierr); - ierr = MatSetFromOptions(*A); - CHKERRQ(ierr); - ierr = MatSetSizes(*A,num_rows_local_,num_rows_local_,num_rows_global_,num_rows_global_); - CHKERRQ(ierr); - ierr = MatSetUp(*A); - CHKERRQ(ierr); - - num_local_states = offdiag_vals_.n_rows; - d_nz.set_size(num_rows_local_); - o_nz.set_size(num_rows_local_); - d_nz.fill(1); - o_nz.fill(0); - for (auto ir: tv_reactions_) - { - for (int i{0}; i < num_local_states; ++i) - { - d_nz(i) += (dblock_nz_(i,ir) - 1); - o_nz(i) += oblock_nz_(i,ir); - } - } - - if (!ti_reactions_.empty()) - { - for (int i{0}; i < num_local_states; ++i) - { - d_nz(i) += (ti_dblock_nz_(i) - 1); - o_nz(i) += ti_oblock_nz_(i); +int FspMatrixConstrained::CreateRHSJacobian(Mat *A) { + int ierr; + + PetscInt own_start, own_end, num_local_states, itmp; + arma::Col d_nz, o_nz, tmp; + + ierr = MatCreate(comm_, A); + CHKERRQ(ierr); + ierr = MatSetFromOptions(*A); + CHKERRQ(ierr); + ierr = MatSetSizes(*A, num_rows_local_, num_rows_local_, num_rows_global_, num_rows_global_); + CHKERRQ(ierr); + ierr = MatSetUp(*A); + CHKERRQ(ierr); + + num_local_states = offdiag_vals_.n_rows; + d_nz.set_size(num_rows_local_); + o_nz.set_size(num_rows_local_); + d_nz.fill(1); + o_nz.fill(0); + for (auto ir: tv_reactions_) { + for (int i{0}; i < num_local_states; ++i) { + d_nz(i) += (dblock_nz_(i, ir) - 1); + o_nz(i) += oblock_nz_(i, ir); + } } - } - - - // Find nnz for rows of the sink states - tmp.set_size(num_constraints_); - tmp.zeros(); - if (rank_ == sinks_rank_) - { - for (auto ir: enable_reactions_) - { - for (int i = 0; i < num_constraints_; ++i) - { - d_nz(num_local_states + i) += sinkmat_nnz(i,ir); - } + + if (!ti_reactions_.empty()) { + for (int i{0}; i < num_local_states; ++i) { + d_nz(i) += (ti_dblock_nz_(i) - 1); + o_nz(i) += ti_oblock_nz_(i); + } } - } else - { - for (auto ir: enable_reactions_) - { - tmp += sinkmat_nnz.col(ir); + + + // Find nnz for rows of the sink states + tmp.set_size(num_constraints_); + tmp.zeros(); + if (rank_ == sinks_rank_) { + for (auto ir: enable_reactions_) { + for (int i = 0; i < num_constraints_; ++i) { + d_nz(num_local_states + i) += sinkmat_nnz(i, ir); + } + } + } else { + for (auto ir: enable_reactions_) { + tmp += sinkmat_nnz.col(ir); + } } - } - ierr = - MPI_Reduce(&tmp[0],o_nz.memptr() + num_local_states,num_constraints_,MPIU_INT,MPIU_SUM,sinks_rank_,comm_); - CHKERRMPI(ierr); + ierr = + MPI_Reduce(&tmp[0], o_nz.memptr() + num_local_states, num_constraints_, MPIU_INT, MPIU_SUM, sinks_rank_, comm_); + CHKERRMPI(ierr); - if (rank_ == sinks_rank_){ - for (int i = 0; i < num_constraints_; ++i) - { - d_nz(num_local_states + i) = std::min(d_nz(num_local_states+i), num_rows_local_); - o_nz(num_local_states + i) = std::min(o_nz(num_local_states+i), num_rows_global_ - num_rows_local_); + if (rank_ == sinks_rank_) { + for (int i = 0; i < num_constraints_; ++i) { + d_nz(num_local_states + i) = std::min(d_nz(num_local_states + i), num_rows_local_); + o_nz(num_local_states + i) = std::min(o_nz(num_local_states + i), num_rows_global_ - num_rows_local_); + } } - } - ierr = MatMPIAIJSetPreallocation(*A,PETSC_NULL,&d_nz[0],PETSC_NULL,&o_nz[0]); - CHKERRQ(ierr); + ierr = MatMPIAIJSetPreallocation(*A, PETSC_NULL, &d_nz[0], PETSC_NULL, &o_nz[0]); + CHKERRQ(ierr); - ierr = VecGetOwnershipRange(work_,&own_start,&own_end); - CHKERRQ(ierr); + ierr = VecGetOwnershipRange(work_, &own_start, &own_end); + CHKERRQ(ierr); - for (auto ir: enable_reactions_) - { - for (int i{0}; i < num_local_states; ++i) - { - ierr = MatSetValue(*A,own_start + i,offdiag_col_idxs_(i,ir),0.0,INSERT_VALUES); - CHKERRQ(ierr); + for (auto ir: enable_reactions_) { + for (int i{0}; i < num_local_states; ++i) { + ierr = MatSetValue(*A, own_start + i, offdiag_col_idxs_(i, ir), 0.0, INSERT_VALUES); + CHKERRQ(ierr); + } } - } - for (int i{0}; i < num_rows_local_; ++i) - { - ierr = MatSetValue(*A,own_start + i,own_start + i,0.0,INSERT_VALUES); - CHKERRQ(ierr); - } - - for (auto i_reaction: enable_reactions_) - { - for (int i_constr{0}; i_constr < num_constraints_; i_constr++) - { - for (int j{0}; j < sinkmat_nnz(i_constr,i_reaction); j++) - { - itmp = num_rows_global_ - num_constraints_ + i_constr; - ierr = MatSetValue(*A,itmp,*(sinkmat_inz.at(i_reaction * num_constraints_ + i_constr).memptr() + j), - 0.0,INSERT_VALUES); + for (int i{0}; i < num_rows_local_; ++i) { + ierr = MatSetValue(*A, own_start + i, own_start + i, 0.0, INSERT_VALUES); CHKERRQ(ierr); - } } - } - MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY); + for (auto i_reaction: enable_reactions_) { + for (int i_constr{0}; i_constr < num_constraints_; i_constr++) { + for (int j{0}; j < sinkmat_nnz(i_constr, i_reaction); j++) { + itmp = num_rows_global_ - num_constraints_ + i_constr; + ierr = MatSetValue(*A, itmp, *(sinkmat_inz.at(i_reaction * num_constraints_ + i_constr).memptr() + j), + 0.0, INSERT_VALUES); + CHKERRQ(ierr); + } + } + } - return 0; + MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY); + + return 0; } -int FspMatrixConstrained::ComputeRHSJacobian(PetscReal t,Mat A) -{ - int ierr = 0; - PetscInt itmp,jtmp; - PetscInt own_start,own_end; +int FspMatrixConstrained::ComputeRHSJacobian(PetscReal t, Mat A) { + int ierr = 0; + PetscInt itmp, jtmp; + PetscInt own_start, own_end; - ierr = FspMatrixBase::ComputeRHSJacobian(t,A); - PACMENSLCHKERRQ(ierr); + ierr = FspMatrixBase::ComputeRHSJacobian(t, A); + PACMENSLCHKERRQ(ierr); - ierr = VecGetOwnershipRange(work_,&own_start,&own_end); - CHKERRQ(ierr); + ierr = VecGetOwnershipRange(work_, &own_start, &own_end); + CHKERRQ(ierr); - PetscReal atmp; - if (!tv_reactions_.empty()) - { - for (int i_reaction: tv_reactions_) - { - for (int i_constr{0}; i_constr < num_constraints_; i_constr++) - { + PetscReal atmp; + if (!tv_reactions_.empty()) { + for (int i_reaction: tv_reactions_) { + for (int i_constr{0}; i_constr < num_constraints_; i_constr++) { - for (int j{0}; j < sinkmat_nnz(i_constr,i_reaction); j++) - { - itmp = num_rows_global_ - num_constraints_ + i_constr; - atmp = time_coefficients_(i_reaction) * sinkmat_entries.at(i_reaction * num_constraints_ + i_constr)(j); - jtmp = sinkmat_inz.at(i_reaction * num_constraints_ + i_constr)(j); + for (int j{0}; j < sinkmat_nnz(i_constr, i_reaction); j++) { + itmp = num_rows_global_ - num_constraints_ + i_constr; + atmp = time_coefficients_(i_reaction) + * sinkmat_entries.at(i_reaction * num_constraints_ + i_constr)(j); + jtmp = sinkmat_inz.at(i_reaction * num_constraints_ + i_constr)(j); - ierr = MatSetValue(A,itmp,jtmp,atmp,ADD_VALUES); - CHKERRQ(ierr); + ierr = MatSetValue(A, itmp, jtmp, atmp, ADD_VALUES); + CHKERRQ(ierr); + } + } } - } } - } - - if (!ti_reactions_.empty()) - { - for (auto i_reaction: ti_reactions_) - { - for (auto i_constr{0}; i_constr < num_constraints_; i_constr++) - { - itmp = num_rows_global_ - num_constraints_ + i_constr; - ierr = MatSetValues(A,1,&itmp,sinkmat_nnz(i_constr,i_reaction), - sinkmat_inz.at(i_reaction * num_constraints_ + i_constr).memptr(), - sinkmat_entries.at(i_reaction * num_constraints_ + i_constr).memptr(), - ADD_VALUES); - CHKERRQ(ierr); - } + + if (!ti_reactions_.empty()) { + for (auto i_reaction: ti_reactions_) { + for (auto i_constr{0}; i_constr < num_constraints_; i_constr++) { + itmp = num_rows_global_ - num_constraints_ + i_constr; + ierr = MatSetValues(A, 1, &itmp, sinkmat_nnz(i_constr, i_reaction), + sinkmat_inz.at(i_reaction * num_constraints_ + i_constr).memptr(), + sinkmat_entries.at(i_reaction * num_constraints_ + i_constr).memptr(), + ADD_VALUES); + CHKERRQ(ierr); + } + } } - } - MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); - return 0; + MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); + return 0; } -PacmenslErrorCode FspMatrixConstrained::GetLocalMVFlops(PetscInt *nflops) -{ - PetscInt ierr; - MatInfo minfo; +PacmenslErrorCode FspMatrixConstrained::GetLocalMVFlops(PetscInt *nflops) { + PetscInt ierr; + MatInfo minfo; - ierr = FspMatrixBase::GetLocalMVFlops(nflops); PACMENSLCHKERRQ(ierr); + ierr = FspMatrixBase::GetLocalMVFlops(nflops); + PACMENSLCHKERRQ(ierr); - ierr = MatGetInfo(ti_sinks_mat_, MAT_LOCAL, &minfo); CHKERRQ(ierr); - *nflops += 2*((PetscInt) minfo.nz_used); + ierr = MatGetInfo(ti_sinks_mat_, MAT_LOCAL, &minfo); + CHKERRQ(ierr); + *nflops += 2 * (( PetscInt ) minfo.nz_used); - for (int j{0}; j < tv_sinks_mat_.size(); ++j){ - ierr = MatGetInfo(tv_sinks_mat_[j], MAT_LOCAL, &minfo); CHKERRQ(ierr); - *nflops += 2*((PetscInt) minfo.nz_used) + num_constraints_; - } + for (int j{0}; j < tv_sinks_mat_.size(); ++j) { + ierr = MatGetInfo(tv_sinks_mat_[j], MAT_LOCAL, &minfo); + CHKERRQ(ierr); + *nflops += 2 * (( PetscInt ) minfo.nz_used) + num_constraints_; + } - return 0; + return 0; } - } \ No newline at end of file diff --git a/src/Models/Model.h b/src/Models/Model.h index 2db01bf..d5cdbfb 100644 --- a/src/Models/Model.h +++ b/src/Models/Model.h @@ -62,13 +62,12 @@ using TcoefFun = std::function stoichiometry_matrix_; ///< Stoichiometry matrix of the reaction network, each column corresponds to a reaction - TcoefFun prop_t_; ///< Callable to evaluate the time-dependent coefficients of the propensity functions. \n - ///< This function must have syntax `int(double t, int num_coefs, double *outputs, void *args)` - void *prop_t_args_; ///< Pointer to extra data (if there is any) needed for the excution of \ref prop_t_ - PropFun prop_x_; ///< Function to evaluate the state-dependent coefficients of the propensity functions.\n - ///< This function must have syntax `int (const int reaction, const int num_species, const int num_states, const int *states, double *outputs, void *args)` - void *prop_x_args_; ///< Pointer to extra data (if there is any) needed for the excution of \ref prop_x_ + arma::Mat + stoichiometry_matrix_; ///< Stoichiometry matrix of the reaction network, each column corresponds to a reaction + TcoefFun prop_t_; ///< Callable to evaluate the time-dependent coefficients of the propensity functions. \n + void *prop_t_args_; ///< Pointer to extra data (if there is any) needed for the excution of \ref prop_t_ + PropFun prop_x_; ///< Function to evaluate the state-dependent coefficients of the propensity functions. + void *prop_x_args_; ///< Pointer to extra data (if there is any) needed for the excution of \ref prop_x_ std::vector tv_reactions_; ///< List of reactions whose propensities are time-varying /** diff --git a/src/Models/SensModel.cpp b/src/Models/SensModel.cpp index 1a90901..6d40191 100644 --- a/src/Models/SensModel.cpp +++ b/src/Models/SensModel.cpp @@ -24,118 +24,34 @@ SOFTWARE. #include "SensModel.h" -pacmensl::SensModel::SensModel(const arma::Mat &stoichiometry_matrix, +pacmensl::SensModel::SensModel(const int num_parameters, + const arma::Mat &stoichiometry_matrix, const std::vector &tv_reactions, - const TcoefFun &prop_t, - const PropFun &prop_x, - const std::vector &dprop_t, - const std::vector &dprop_x, - const std::vector &dprop_ic, - const std::vector &dprop_rowptr, + const pacmensl::TcoefFun &prop_t, + const pacmensl::PropFun &prop_x, + const pacmensl::DTcoefFun &dprop_t, + const std::vector> &dprop_t_sp, + const pacmensl::DPropFun &dprop_x, + const std::vector> &dprop_x_sp, void *prop_t_args, void *prop_x_args, - const std::vector &dprop_t_args, - const std::vector &dprop_x_args) -{ - num_reactions_ = stoichiometry_matrix.n_cols; - num_parameters_ = dprop_x.size(); + void *dprop_t_args, + void *dprop_x_args) { + num_parameters_ = num_parameters; + num_reactions_ = stoichiometry_matrix.n_cols; stoichiometry_matrix_ = stoichiometry_matrix; - prop_t_ = prop_t; - prop_x_ = prop_x; - dprop_t_ = dprop_t; - dprop_x_ = dprop_x; - dpropensity_ic_ = dprop_ic; - dpropensity_rowptr_ = dprop_rowptr; - if (dpropensity_rowptr_.size() == num_parameters_) dpropensity_rowptr_.push_back(dpropensity_ic_.size()); - + prop_t_ = prop_t; + prop_x_ = prop_x; + dprop_t_ = dprop_t; + dprop_t_sp_ = dprop_t_sp; + dprop_x_ = dprop_x; + dprop_x_sp_ = dprop_x_sp; + prop_t_args_ = prop_t_args; prop_x_args_ = prop_x_args; - if (dprop_t_args.empty()) { - dprop_t_args_ = std::vector(num_parameters_, nullptr); - } else { - dprop_t_args_ = dprop_t_args; - } - if (dprop_x_args.empty()) { - dprop_x_args_ = std::vector(num_parameters_, nullptr); - } else { - dprop_x_args_ = dprop_x_args; - } - + dprop_t_args_ = dprop_t_args; + dprop_x_args_ = dprop_x_args; + tv_reactions_ = tv_reactions; -} - -pacmensl::SensModel::SensModel(const pacmensl::SensModel &model) { - num_reactions_ = model.stoichiometry_matrix_.n_cols; - num_parameters_ = model.dprop_x_.size(); - stoichiometry_matrix_ = model.stoichiometry_matrix_; - prop_t_ = model.prop_t_; - prop_x_ = model.prop_x_; - prop_t_args_ = model.prop_t_args_; - prop_x_args_ = model.prop_x_args_; - dprop_t_ = model.dprop_t_; - dprop_x_ = model.dprop_x_; - dprop_t_args_ = model.dprop_t_args_; - dprop_x_args_ = model.dprop_x_args_; - dpropensity_ic_ = model.dpropensity_ic_; - dpropensity_rowptr_ = model.dpropensity_rowptr_; - tv_reactions_ = model.tv_reactions_; -} - -pacmensl::SensModel &pacmensl::SensModel::operator=(const pacmensl::SensModel &model) noexcept { - if (this == &model) return *this; - - stoichiometry_matrix_.clear(); - dprop_t_.clear(); - dprop_x_.clear(); - dpropensity_ic_.clear(); - dpropensity_rowptr_.clear(); - dprop_t_args_.clear(); - dprop_x_args_.clear(); - - num_reactions_ = model.stoichiometry_matrix_.n_cols; - num_parameters_ = model.dprop_x_.size(); - stoichiometry_matrix_ = model.stoichiometry_matrix_; - prop_t_ = model.prop_t_; - prop_x_ = model.prop_x_; - prop_t_args_ = model.prop_t_args_; - prop_x_args_ = model.prop_x_args_; - dprop_t_ = model.dprop_t_; - dprop_x_ = model.dprop_x_; - dprop_t_args_ = model.dprop_t_args_; - dprop_x_args_ = model.dprop_x_args_; - dpropensity_ic_ = model.dpropensity_ic_; - dpropensity_rowptr_ = model.dpropensity_rowptr_; - - tv_reactions_ = model.tv_reactions_; - return *this; -} - -pacmensl::SensModel &pacmensl::SensModel::operator=(pacmensl::SensModel &&model) noexcept { - if (this == &model) return *this; - - stoichiometry_matrix_.clear(); - dprop_x_args_.clear(); - dprop_x_.clear(); - dprop_t_.clear(); - dprop_t_args_.clear(); - dpropensity_ic_.clear(); - dpropensity_rowptr_.clear(); - - - num_reactions_ = std::move(model.stoichiometry_matrix_.n_cols); - num_parameters_ = std::move(model.dprop_x_.size()); - stoichiometry_matrix_ = std::move(model.stoichiometry_matrix_); - prop_t_ = std::move(model.prop_t_); - prop_x_ = std::move(model.prop_x_); - prop_t_args_ = std::move(model.prop_t_args_); - prop_x_args_ = std::move(model.prop_x_args_); - dprop_t_ = std::move(model.dprop_t_); - dprop_x_ = std::move(model.dprop_x_); - dprop_t_args_ = std::move(model.dprop_t_args_); - dprop_x_args_ = std::move(model.dprop_x_args_); - dpropensity_ic_ = model.dpropensity_ic_; - dpropensity_rowptr_ = model.dpropensity_rowptr_; - tv_reactions_ = std::move(model.tv_reactions_); - - return *this; + } diff --git a/src/Models/SensModel.h b/src/Models/SensModel.h index 1c3e366..89fd567 100644 --- a/src/Models/SensModel.h +++ b/src/Models/SensModel.h @@ -29,52 +29,71 @@ SOFTWARE. #include "Sys.h" #include "Model.h" +/** + * @file SensModel.h + */ + namespace pacmensl { +using DTcoefFun = std::function; + +using DPropFun = std::function; + /* * Class for representing a stochastic chemical reaction network model with sensitivity information. */ -class SensModel -{ +class SensModel { public: - int num_reactions_ = 0; - int num_parameters_ = 0; - arma::Mat stoichiometry_matrix_; - // Which reactions have time-varying propensities? - std::vector tv_reactions_; - // propensities - PropFun prop_x_; - void *prop_x_args_; - TcoefFun prop_t_; - void *prop_t_args_; - // derivatives of propensities - std::vector dprop_x_; - std::vector dprop_t_; - std::vector dprop_x_args_; - std::vector dprop_t_args_; - std::vector dpropensity_ic_; - std::vector dpropensity_rowptr_; + int num_reactions_ = 0; ///< number of reactions + int num_parameters_ = 0; ///< number of parameters + arma::Mat stoichiometry_matrix_; ///< stoichiometry matrix, reactions are ordered column-wise - SensModel() {}; + std::vector tv_reactions_; ///< container for reaction indices with time-varying propensities? - explicit SensModel(const arma::Mat &stoichiometry_matrix, - const std::vector &tv_reactions, - const TcoefFun &prop_t, - const PropFun &prop_x, - const std::vector &dprop_t, - const std::vector &dprop_x, - const std::vector &dprop_ic = std::vector(), - const std::vector &dprop_rowptr = std::vector(), - void *prop_t_args = nullptr, - void *prop_x_args = nullptr, - const std::vector &dprop_t_args = std::vector(), - const std::vector &dprop_x_args = std::vector()); + // propensities + PropFun prop_x_; ///< function to evaluate state-dependent factors + void *prop_x_args_; ///< pointer to extra arguments needed by \ref prop_x_ + TcoefFun prop_t_; ///< function to evaluate time-dependent factors + void *prop_t_args_; ///< pointer to extra arguments needed by \ref prop_t_ - SensModel(const SensModel &model); + // derivatives of propensities + DTcoefFun dprop_t_; + void *dprop_t_args_; + std::vector> dprop_t_sp_; - SensModel &operator=(const SensModel &model) noexcept; + DPropFun dprop_x_; + void *dprop_x_args_; + std::vector> dprop_x_sp_; + + SensModel() {}; - SensModel &operator=(SensModel &&model) noexcept; + explicit SensModel( + const int num_parameters, + const arma::Mat &stoichiometry_matrix, + const std::vector &tv_reactions, + const TcoefFun &prop_t, + const PropFun &prop_x, + const DTcoefFun &dprop_t, + const std::vector> &dprop_t_sp, + const DPropFun &dprop_x, + const std::vector> &dprop_x_sp = std::vector>(), + void *prop_t_args = nullptr, + void *prop_x_args = nullptr, + void *dprop_t_args = nullptr, + void *dprop_x_args = nullptr); }; } diff --git a/src/OdeSolver/CvodeFsp.h b/src/OdeSolver/CvodeFsp.h index 49c597d..2259ac2 100644 --- a/src/OdeSolver/CvodeFsp.h +++ b/src/OdeSolver/CvodeFsp.h @@ -35,13 +35,12 @@ SOFTWARE. #include "StateSetConstrained.h" #include "Sys.h" - namespace pacmensl { class CvodeFsp : public OdeSolverBase { public: explicit CvodeFsp(MPI_Comm _comm, int lmm = CV_BDF); - PacmenslErrorCode SetUp() override ; + PacmenslErrorCode SetUp() override; PetscInt Solve() override; diff --git a/src/SensFsp/SensFspMatrix.h b/src/SensFsp/SensFspMatrix.h index 4e5157c..78e3346 100644 --- a/src/SensFsp/SensFspMatrix.h +++ b/src/SensFsp/SensFspMatrix.h @@ -28,6 +28,7 @@ SOFTWARE. #include "SensModel.h" #include "FspMatrixBase.h" #include "FspMatrixConstrained.h" +#include "PetscWrap.h" namespace pacmensl { @@ -35,14 +36,13 @@ namespace pacmensl { * @brief Templated class for the time-dependent matrix on the right hand side of the finite state projection-based forward sensitivity * equation. * @details We currently assume that the CME matrix could be decomposed into the form - * \f$ A(t) = \sum_{r=1}^{M}{c_r(t, \theta)A_r} \f$ - * where c_r(t,\theta) are scalar-valued functions that depend on the time variable and parameters, while the matrices \f$ A_r \f$ are constant. + * \f$ A(t) = \sum_{r=1}^{M}{c_r(t, \theta)A_r(\theta)} \f$ + * where c_r(t,\theta) are scalar-valued functions that depend on the time variable and parameters, while the matrices \f$ A_r(\theta) \f$ depend only on parameters. * As a consequence, the derivative of \f$ A(t) \f$ with respect to the \f$ j \f$-th parameter component is assumed to have the form - * \f$ \partial_{j}A(t) = \sum_{r=1}^{M}{\partial_j{c_r}(t,\theta)}A_r \f$ + * \f$ \partial_{j}A(t) = \sum_{r=1}^{M}{\partial_j{c_r}(t,\theta)}A_r(\theta) \f$ */ template -class SensFspMatrix -{ +class SensFspMatrix { public: NOT_COPYABLE_NOT_MOVABLE(SensFspMatrix); @@ -56,20 +56,20 @@ class SensFspMatrix int GetNumLocalRows() const; protected: MPI_Comm comm_ = MPI_COMM_NULL; - int rank_ = 0; + int rank_ = 0; - int num_parameters_ = 0; ///< Number of sensitivity parameters - FspMatrixT A_; ///< The FSP-truncated infinitesimal generator - std::vector dA_; ///< The derivatives of A with respect to sensitivity parameters + int num_parameters_ = 0; ///< Number of sensitivity parameters + FspMatrixT A_; ///< The FSP-truncated infinitesimal generator + std::vector dcxA_; ///< The derivatives of A with respect to sensitivity parameters + std::vector cxdA_; + Petsc work_; ///< Scratch vector for performing matrix-vector multiplication }; } template -pacmensl::SensFspMatrix::SensFspMatrix(MPI_Comm comm) : A_(comm) -{ +pacmensl::SensFspMatrix::SensFspMatrix(MPI_Comm comm) : A_(comm) { int ierr; - if (comm == MPI_COMM_NULL) - { + if (comm == MPI_COMM_NULL) { std::cout << "Null pointer detected.\n"; } comm_ = comm; @@ -78,93 +78,136 @@ pacmensl::SensFspMatrix::SensFspMatrix(MPI_Comm comm) : A_(comm) } template -PacmenslErrorCode pacmensl::SensFspMatrix::Destroy() -{ +PacmenslErrorCode pacmensl::SensFspMatrix::Destroy() { PacmenslErrorCode ierr; ierr = A_.Destroy(); PACMENSLCHKERRQ(ierr); - for (int i{0}; i < num_parameters_; ++i) - { - ierr = dA_[i].Destroy(); + for (int i{0}; i < num_parameters_; ++i) { + ierr = dcxA_[i].Destroy(); + PACMENSLCHKERRQ(ierr); + ierr = cxdA_[i].Destroy(); PACMENSLCHKERRQ(ierr); } - dA_.clear(); + dcxA_.clear(); + cxdA_.clear(); return 0; } template -pacmensl::SensFspMatrix::~SensFspMatrix() -{ +pacmensl::SensFspMatrix::~SensFspMatrix() { Destroy(); comm_ = MPI_COMM_NULL; } template -PacmenslErrorCode pacmensl::SensFspMatrix::GenerateValues(const pacmensl::StateSetBase &state_set, - const pacmensl::SensModel &model) -{ +PacmenslErrorCode +pacmensl::SensFspMatrix::GenerateValues(const pacmensl::StateSetBase &state_set, + const pacmensl::SensModel &model) { PacmenslErrorCode ierr; - ierr = A_.GenerateValues(state_set, model.stoichiometry_matrix_, model.tv_reactions_, model.prop_t_, model.prop_x_, - std::vector(), model.prop_t_args_, model.prop_x_args_); + ierr = A_.GenerateValues(state_set, + model.stoichiometry_matrix_, + model.tv_reactions_, + model.prop_t_, + model.prop_x_, + std::vector(), + model.prop_t_args_, + model.prop_x_args_); PACMENSLCHKERRQ(ierr); + + ierr = VecCreate(comm_, work_.mem()); + CHKERRQ(ierr); + ierr = VecSetSizes(work_, A_.GetNumLocalRows(), PETSC_DECIDE); + CHKERRQ(ierr); + ierr = VecSetUp(work_); + CHKERRQ(ierr); + num_parameters_ = model.num_parameters_; - dA_.reserve(num_parameters_); - if (!model.dpropensity_ic_.empty()) - { - for (int i{0}; i < num_parameters_; ++i) - { - auto first = model.dpropensity_ic_.begin() + model.dpropensity_rowptr_[i]; - auto last = model.dpropensity_ic_.begin() + model.dpropensity_rowptr_[i + 1]; - std::vector enable_reactions(first, last); - dA_.emplace_back(FspMatrixT(comm_)); - ierr = dA_[i].GenerateValues(state_set, - model.stoichiometry_matrix_, - model.tv_reactions_, - model.dprop_t_[i], - model.dprop_x_[i], - enable_reactions, - model.dprop_t_args_[i], - model.dprop_x_args_[i]); - PACMENSLCHKERRQ(ierr); + + // Generate the first part of the derivative of A wrt parameters, consisting of the terms dc*A + dcxA_.reserve(num_parameters_); + if (!model.dprop_t_sp_.empty()) { + for (int i{0}; i < num_parameters_; ++i) { + dcxA_.emplace_back(FspMatrixT(comm_)); + if (!model.dprop_t_sp_[i].empty()) { + auto dprop_t = std::bind(model.dprop_t_, i, + std::placeholders::_1, + std::placeholders::_2, + std::placeholders::_3, + std::placeholders::_4); + ierr = dcxA_[i].GenerateValues(state_set, + model.stoichiometry_matrix_, + model.tv_reactions_, + dprop_t, + model.prop_x_, + model.dprop_t_sp_[i], + model.dprop_t_args_, + model.prop_x_args_); + PACMENSLCHKERRQ(ierr); + } + } + } else { + for (int i{0}; i < num_parameters_; ++i) { + dcxA_.emplace_back(FspMatrixT(comm_)); + } + } + + // Generate the second part of the derivative of A wrt parameters, consisting of the terms c*dA + cxdA_.reserve(num_parameters_); + if (!model.dprop_x_sp_.empty()) { + for (int i{0}; i < num_parameters_; ++i) { + cxdA_.emplace_back(FspMatrixT(comm_)); + if (!model.dprop_x_sp_[i].empty()) { + auto dprop_x = std::bind(model.dprop_x_, i, + std::placeholders::_1, + std::placeholders::_2, + std::placeholders::_3, + std::placeholders::_4, + std::placeholders::_5, + std::placeholders::_6); + + ierr = cxdA_[i].GenerateValues(state_set, + model.stoichiometry_matrix_, + model.tv_reactions_, + model.prop_t_, + dprop_x, + model.dprop_x_sp_[i], + model.prop_t_args_, + model.dprop_x_args_); + PACMENSLCHKERRQ(ierr); + } } - } else - { - for (int i{0}; i < num_parameters_; ++i) - { - dA_.emplace_back(FspMatrixT(comm_)); - ierr = dA_[i].GenerateValues(state_set, - model.stoichiometry_matrix_, - model.tv_reactions_, - model.dprop_t_[i], - model.dprop_x_[i], - std::vector(), - model.dprop_t_args_[i], - model.dprop_x_args_[i]); - PACMENSLCHKERRQ(ierr); + } else { + for (int i{0}; i < num_parameters_; ++i) { + cxdA_.emplace_back(FspMatrixT(comm_)); } } return 0; } template -PacmenslErrorCode pacmensl::SensFspMatrix::Action(PetscReal t, Vec x, Vec y) -{ +PacmenslErrorCode +pacmensl::SensFspMatrix::Action(PetscReal t, Vec x, Vec y) { return A_.Action(t, x, y); } template -PacmenslErrorCode pacmensl::SensFspMatrix::SensAction(int i_par, PetscReal t, Vec x, Vec y) -{ - if (i_par < 0 || i_par >= num_parameters_) - { +PacmenslErrorCode +pacmensl::SensFspMatrix::SensAction(int i_par, PetscReal t, Vec x, Vec y) { + int ierr; + if (i_par < 0 || i_par >= num_parameters_) { return -1; } - return dA_[i_par].Action(t, x, y); + ierr = dcxA_[i_par].Action(t, x, y); + PACMENSLCHKERRQ(ierr); + ierr = cxdA_[i_par].Action(t, x, work_); + PACMENSLCHKERRQ(ierr); + ierr = VecAXPY(y, 1.0, work_); + CHKERRQ(ierr); + return 0; } template -int pacmensl::SensFspMatrix::GetNumLocalRows() const -{ +int pacmensl::SensFspMatrix::GetNumLocalRows() const { return A_.GetNumLocalRows(); } diff --git a/src/SensFsp/SensFspSolverMultiSinks.cpp b/src/SensFsp/SensFspSolverMultiSinks.cpp index 5ec4d2b..3b0c706 100644 --- a/src/SensFsp/SensFspSolverMultiSinks.cpp +++ b/src/SensFsp/SensFspSolverMultiSinks.cpp @@ -125,10 +125,6 @@ PacmenslErrorCode pacmensl::SensFspSolverMultiSinks::SetUp() // Make sure all the necessary parameters have been set try { - if (model_.prop_t_ == nullptr) - { - throw std::runtime_error("Temporal signal was not set before calling FspSolver.SetUp()."); - } if (model_.prop_x_ == nullptr) { throw std::runtime_error("Propensity was not set before calling FspSolver.SetUp()."); @@ -137,10 +133,6 @@ PacmenslErrorCode pacmensl::SensFspSolverMultiSinks::SetUp() { throw std::runtime_error("Empty stoichiometry matrix cannot be used for FspSolver."); } - if (model_.dprop_x_.empty()) - { - throw std::runtime_error("Empty sensitivity information.\n"); - } if (init_states_.n_elem == 0 || init_probs_.n_elem == 0) { throw std::runtime_error("Initial states and/or probabilities were not set before calling FspSolver.SetUp()."); diff --git a/tests/test_sensfsp_solver.cpp b/tests/test_sensfsp_solver.cpp index 4d65c49..fc30ee3 100644 --- a/tests/test_sensfsp_solver.cpp +++ b/tests/test_sensfsp_solver.cpp @@ -1 +1 @@ -/* MIT License Copyright (c) 2020 Huy Vo Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ static char help[] = "Test interface to CVODE for solving the CME of the toggle model.\n\n"; #include #include "SensFspSolverMultiSinks.h" #include "pacmensl_test_env.h" namespace toggle_cme { /* Stoichiometric matrix of the toggle switch model */ arma::Mat SM{{1, 1, -1, 0, 0, 0}, {0, 0, 0, 1, 1, -1}}; const int nReaction = 6; /* Parameters for the propensity functions */ const double ayx{2.6e-3}, axy{6.1e-3}, nyx{3.0e0}, nxy{2.1e0}, kx0{2.2e-3}, kx{1.7e-2}, dx{3.8e-4}, ky0{6.8e-5}, ky{ 1.6e-2}, dy{3.8e-4}; // Function to constraint the shape of the Fsp void lhs_constr(PetscInt num_species, PetscInt num_constrs, PetscInt num_states, PetscInt *states, int *vals, void *args) { for (int i{0}; i < num_states; ++i) { vals[i * num_constrs] = states[num_species * i]; vals[i * num_constrs + 1] = states[num_species * i + 1]; vals[i * num_constrs + 2] = states[num_species * i] * states[num_species * i + 1]; } } arma::Row rhs_constr{200, 200, 2000}; arma::Row expansion_factors{0.2, 0.2, 0.2}; // propensity function for toggle int propensity(const int reaction, const int num_species, const int num_states, const PetscInt *X, double *outputs, void *args) { int (*X_view)[2] = ( int (*)[2] ) X; switch (reaction) { case 0:for (int i{0}; i < num_states; ++i) { outputs[i] = 1.0; } break; case 1: for (int i{0}; i < num_states; ++i) { outputs[i] = 1.0 / (1.0 + ayx * pow(PetscReal(X_view[i][1]), nyx)); } break; case 2:for (int i{0}; i < num_states; ++i) { outputs[i] = PetscReal(X_view[i][0]); } break; case 3:for (int i{0}; i < num_states; ++i) { outputs[i] = 1.0; } break; case 4: for (int i{0}; i < num_states; ++i) { outputs[i] = 1.0 / (1.0 + axy * pow(PetscReal(X_view[i][0]), nxy)); } break; case 5:for (int i{0}; i < num_states; ++i) { outputs[i] = PetscReal(X_view[i][1]); } break; default:return -1; } return 0; } int t_fun(PetscReal t, int n_coefs, double *outputs, void *args) { outputs[0] = kx0; outputs[1] = kx; outputs[2] = dx; outputs[3] = ky0; outputs[4] = ky; outputs[5] = dy; return 0; } int get_sens_t_fun(int i, PetscReal t, int n_coefs, double *outputs, void *args) { outputs[i] = 1.0; return 0; } } using namespace pacmensl; class SensFspTest : public ::testing::Test { protected: SensFspTest() {} void SetUp() override { int n_par = 6; t_final = 100.0; fsp_tol = 1.0e-10; X0 = X0.t(); dp0 = std::vector>(n_par, arma::Col({0.0})); std::vector d_t_fun(n_par); for (int i{0}; i < n_par; ++i) { d_t_fun[i] = std::bind(toggle_cme::get_sens_t_fun, i, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4); } toggle_model = SensModel(toggle_cme::SM, std::vector({0,1,2,3,4,5}), toggle_cme::t_fun, toggle_cme::propensity, d_t_fun, std::vector(n_par, toggle_cme::propensity), std::vector({0, 1, 2, 3, 4, 5}), std::vector({0, 1, 2, 3, 4, 5, 6})); } void TearDown() override { } PetscReal t_final, fsp_tol; arma::Mat X0{0, 0}; arma::Col p0 = {1.0}; std::vector> dp0; SensModel toggle_model; arma::Row fsp_size = {1, 1}; arma::Row expansion_factors = {0.25, 0.25}; }; TEST_F(SensFspTest, test_wrong_call_sequence_detection) { int ierr; SensFspSolverMultiSinks fsp(PETSC_COMM_WORLD); ierr = fsp.SetUp(); ASSERT_EQ(ierr, -1); } // //TEST_F(SensFspTest, test_handling_t_fun_error) //{ // int ierr; // // SensFspSolverMultiSinks fsp(PETSC_COMM_WORLD); // std::vector tspan = // arma::conv_to>::from(arma::linspace>( // 0.0, // t_final, // 3)); // SensModel bad_model = toggle_model; // bad_model.prop_t_ = [&](double t, int n, double *vals, void *args) { // return -1; // }; // // ierr = fsp.SetModel(bad_model); // ASSERT_FALSE(ierr); // ierr = fsp.SetInitialBounds(fsp_size); // ASSERT_FALSE(ierr); // ierr = fsp.SetExpansionFactors(expansion_factors); // ASSERT_FALSE(ierr); // ierr = fsp.SetVerbosity(0); // ASSERT_FALSE(ierr); // ierr = fsp.SetInitialDistribution(X0, p0, dp0); // ASSERT_FALSE(ierr); // // ierr = fsp.SetUp(); // ASSERT_FALSE(ierr); // ASSERT_THROW(fsp.Solve(t_final, fsp_tol), std::runtime_error); // fsp.ClearState(); //} TEST_F(SensFspTest, toggle_sens_solve_with_cvode) { PetscInt ierr; PetscReal stmp; SensDiscreteDistribution p_final_bdf; std::vector p_snapshots_bdf; std::vector tspan; Vec q; tspan = arma::conv_to>::from(arma::linspace>(0.0, t_final, 3)); SensFspSolverMultiSinks fsp(PETSC_COMM_WORLD); ierr = fsp.SetModel(toggle_model); ASSERT_FALSE(ierr); ierr = fsp.SetInitialBounds(fsp_size); ASSERT_FALSE(ierr); ierr = fsp.SetExpansionFactors(expansion_factors); ASSERT_FALSE(ierr); ierr = fsp.SetInitialDistribution(X0, p0, dp0); ASSERT_FALSE(ierr); ierr = fsp.SetUp(); ASSERT_FALSE(ierr); p_final_bdf = fsp.Solve(t_final, fsp_tol); fsp.ClearState(); ierr = VecSum(p_final_bdf.p_, &stmp); ASSERT_FALSE(ierr); ASSERT_GE(stmp, 1.0-fsp_tol); for (int i{0}; i < 6; ++i) { ierr = VecSum(p_final_bdf.dp_[i], &stmp); ASSERT_FALSE(ierr); ASSERT_LE(abs(stmp), 1.0e-6); } // arma::Mat FIM; // ComputeFIM(p_final_bdf, FIM); // std::cout << FIM; } class SensFspPoissonTest : public ::testing::Test { protected: SensFspPoissonTest() {} void SetUp() override { auto propensity = [&](int reaction, int num_species, int num_states, const int *state, PetscReal *output, void *args) { for (int i{0}; i < num_states; ++i) { output[i] = 1.0; } return 0; }; auto t_fun = [&](double t, int num_coefs, double *outputs, void *args) { outputs[0] = lambda; return 0; }; auto d_t_fun = [&](double t, int num_coefs, double *outputs, void *args) { outputs[0] = 1.0; return 0; }; poisson_model = SensModel(stoich_matrix, std::vector({0}), t_fun, propensity, std::vector({d_t_fun}), std::vector({propensity})); } void TearDown() override { } ~SensFspPoissonTest() {} SensModel poisson_model; PetscReal lambda = 2.0; arma::Mat stoich_matrix = {1}; arma::Mat x0 = {0}; arma::Col p0 = {1.0}; arma::Col s0 = {0.0}; arma::Row fsp_size = {5}; arma::Row expansion_factors = {0.1}; PetscReal t_final{1.0}, fsp_tol{1.0e-7}; }; TEST_F(SensFspPoissonTest, test_poisson_analytic) { PetscInt ierr; PetscReal stmp; SensDiscreteDistribution p_final; SensFspSolverMultiSinks fsp(PETSC_COMM_WORLD); ierr = fsp.SetModel(poisson_model); ASSERT_FALSE(ierr); ierr = fsp.SetInitialBounds(fsp_size); ASSERT_FALSE(ierr); ierr = fsp.SetExpansionFactors(expansion_factors); ASSERT_FALSE(ierr); ierr = fsp.SetInitialDistribution(x0, p0, std::vector>({s0})); ASSERT_FALSE(ierr); ierr = fsp.SetUp(); ASSERT_FALSE(ierr); p_final = fsp.Solve(t_final, fsp_tol); fsp.ClearState(); // Check that the solution is close to Poisson stmp = 0.0; PetscReal *p_dat; int num_states; p_final.GetProbView(num_states, p_dat); PetscReal pdf; int n; for (int i = 0; i < num_states; ++i) { n = p_final.states_(0, i); pdf = exp(-lambda * t_final) * pow(lambda * t_final, double(n)) / tgamma(n + 1); stmp += abs(p_dat[i] - pdf); } p_final.RestoreProbView(p_dat); MPI_Allreduce(&stmp, MPI_IN_PLACE, 1, MPIU_REAL, MPIU_SUM, PETSC_COMM_WORLD); ASSERT_LE(stmp, fsp_tol); PetscReal *s_dat; p_final.GetSensView(0, num_states, s_dat); for (int i = 0; i < num_states; ++i) { n = p_final.states_(0, i); pdf = -t_final*exp(-lambda*t_final)*pow(lambda*t_final, double(n))/tgamma(n+1) + exp(-lambda*t_final)*t_final*pow(lambda*t_final, double(n-1))/tgamma(n); stmp += abs(s_dat[i] - pdf); } p_final.RestoreSensView(0, s_dat); MPI_Allreduce(&stmp, MPI_IN_PLACE, 1, MPIU_REAL, MPIU_SUM, PETSC_COMM_WORLD); PetscPrintf(PETSC_COMM_WORLD, "Sensitivity error = %.2e \n", stmp); ASSERT_LE(stmp, 1.0e-6); } \ No newline at end of file +/* MIT License Copyright (c) 2020 Huy Vo Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ static char help[] = "Test the sensitivity FSP implementation.\n\n"; #include #include "SensFspSolverMultiSinks.h" #include "pacmensl_test_env.h" using namespace pacmensl; /* * ===================================================================================================================== * TEST FIXTURE: Toggle-switch with time-varying propensities. Only the time-varying factors of the propensity functions are * parameter-dependent. */ namespace toggle_cme { /* Stoichiometric matrix of the toggle switch model */ arma::Mat SM{{1, 1, -1, 0, 0, 0}, {0, 0, 0, 1, 1, -1}}; const int nReaction = 6; /* Parameters for the propensity functions */ const double ayx{2.6e-3}, axy{6.1e-3}, nyx{3.0e0}, nxy{2.1e0}, kx0{2.2e-3}, kx{1.7e-2}, dx{3.8e-4}, ky0{6.8e-5}, ky{1.6e-2}, dy{3.8e-4}; // Function to constraint the shape of the Fsp void lhs_constr(PetscInt num_species, PetscInt num_constrs, PetscInt num_states, PetscInt *states, int *vals, void *args) { for (int i{0}; i < num_states; ++i) { vals[i * num_constrs] = states[num_species * i]; vals[i * num_constrs + 1] = states[num_species * i + 1]; vals[i * num_constrs + 2] = states[num_species * i] * states[num_species * i + 1]; } } arma::Row rhs_constr{200, 200, 2000}; arma::Row expansion_factors{0.2, 0.2, 0.2}; // propensity function for toggle int propensity(const int reaction, const int num_species, const int num_states, const PetscInt *X, double *outputs, void *args) { int (*X_view)[2] = ( int (*)[2] ) X; switch (reaction) { case 0: for (int i{0}; i < num_states; ++i) { outputs[i] = 1.0; } break; case 1: for (int i{0}; i < num_states; ++i) { outputs[i] = 1.0 / (1.0 + ayx * pow(PetscReal(X_view[i][1]), nyx)); } break; case 2: for (int i{0}; i < num_states; ++i) { outputs[i] = PetscReal(X_view[i][0]); } break; case 3: for (int i{0}; i < num_states; ++i) { outputs[i] = 1.0; } break; case 4: for (int i{0}; i < num_states; ++i) { outputs[i] = 1.0 / (1.0 + axy * pow(PetscReal(X_view[i][0]), nxy)); } break; case 5: for (int i{0}; i < num_states; ++i) { outputs[i] = PetscReal(X_view[i][1]); } break; default:return -1; } return 0; } int t_fun(PetscReal t, int n_coefs, double *outputs, void *args) { outputs[0] = kx0; outputs[1] = kx; outputs[2] = dx; outputs[3] = ky0; outputs[4] = ky; outputs[5] = dy; return 0; } int dt_fun(int parameter_idx, PetscReal t, int n_coefs, double *outputs, void *args) { outputs[parameter_idx] = 1.0; return 0; } } class SensFspToggleTest : public ::testing::Test { protected: SensFspToggleTest() {} void SetUp() override { int n_par = 6; t_final = 100.0; fsp_tol = 1.0e-10; X0 = X0.t(); dp0 = std::vector>(n_par, arma::Col({0.0})); toggle_model = SensModel(6, toggle_cme::SM, std::vector({0, 1, 2, 3, 4, 5}), toggle_cme::t_fun, toggle_cme::propensity, toggle_cme::dt_fun, {{0}, {1}, {2}, {3}, {4}, {5}}, nullptr, {}); } void TearDown() override { } PetscReal t_final, fsp_tol; arma::Mat X0{0, 0}; arma::Col p0 = {1.0}; std::vector> dp0; SensModel toggle_model; arma::Row fsp_size = {1, 1}; arma::Row expansion_factors = {0.25, 0.25}; }; TEST_F(SensFspToggleTest, toggle_sens_solve_with_cvode) { PetscInt ierr; PetscReal stmp; SensDiscreteDistribution p_final_bdf; std::vector p_snapshots_bdf; std::vector tspan; Vec q; tspan = arma::conv_to>::from(arma::linspace>(0.0, t_final, 3)); SensFspSolverMultiSinks fsp(PETSC_COMM_WORLD); ierr = fsp.SetModel(toggle_model); ASSERT_FALSE(ierr); ierr = fsp.SetInitialBounds(fsp_size); ASSERT_FALSE(ierr); ierr = fsp.SetExpansionFactors(expansion_factors); ASSERT_FALSE(ierr); ierr = fsp.SetInitialDistribution(X0, p0, dp0); ASSERT_FALSE(ierr); ierr = fsp.SetUp(); ASSERT_FALSE(ierr); p_final_bdf = fsp.Solve(t_final, fsp_tol); fsp.ClearState(); ierr = VecSum(p_final_bdf.p_, &stmp); ASSERT_FALSE(ierr); ASSERT_GE(stmp, 1.0 - fsp_tol); for (int i{0}; i < 6; ++i) { ierr = VecSum(p_final_bdf.dp_[i], &stmp); ASSERT_FALSE(ierr); ASSERT_LE(abs(stmp), 1.0e-6); } } /* * ===================================================================================================================== * TEST FIXTURE: Pure birth model with constant-birth rate. */ class SensFspPoissonTest : public ::testing::Test { protected: SensFspPoissonTest() {} void SetUp() override { auto propensity = [&](int reaction, int num_species, int num_states, const int *state, PetscReal *output, void *args) { for (int i{0}; i < num_states; ++i) { output[i] = 1.0; } return 0; }; auto t_fun = [&](double t, int num_coefs, double *outputs, void *args) { outputs[0] = lambda; return 0; }; auto d_t_fun = [&](int par_idx, double t, int num_coefs, double *outputs, void *args) { outputs[0] = 1.0; return 0; }; poisson_model = SensModel(1, stoich_matrix, std::vector({0}), t_fun, propensity, d_t_fun, {{0}}, nullptr); } void TearDown() override { } ~SensFspPoissonTest() {} SensModel poisson_model; PetscReal lambda = 2.0; arma::Mat stoich_matrix = {1}; arma::Mat x0 = {0}; arma::Col p0 = {1.0}; arma::Col s0 = {0.0}; arma::Row fsp_size = {5}; arma::Row expansion_factors = {0.1}; PetscReal t_final{1.0}, fsp_tol{1.0e-7}; }; TEST_F(SensFspPoissonTest, test_poisson_analytic) { PetscInt ierr; PetscReal stmp; SensDiscreteDistribution p_final; SensFspSolverMultiSinks fsp(PETSC_COMM_WORLD); ierr = fsp.SetModel(poisson_model); ASSERT_FALSE(ierr); ierr = fsp.SetInitialBounds(fsp_size); ASSERT_FALSE(ierr); ierr = fsp.SetExpansionFactors(expansion_factors); ASSERT_FALSE(ierr); ierr = fsp.SetInitialDistribution(x0, p0, std::vector>({s0})); ASSERT_FALSE(ierr); ierr = fsp.SetUp(); ASSERT_FALSE(ierr); p_final = fsp.Solve(t_final, fsp_tol); fsp.ClearState(); // Check that the solution is close to Poisson stmp = 0.0; PetscReal *p_dat; int num_states; p_final.GetProbView(num_states, p_dat); PetscReal pdf; int n; for (int i = 0; i < num_states; ++i) { n = p_final.states_(0, i); pdf = exp(-lambda * t_final) * pow(lambda * t_final, double(n)) / tgamma(n + 1); stmp += abs(p_dat[i] - pdf); } p_final.RestoreProbView(p_dat); MPI_Allreduce(&stmp, MPI_IN_PLACE, 1, MPIU_REAL, MPIU_SUM, PETSC_COMM_WORLD); ASSERT_LE(stmp, fsp_tol); PetscReal *s_dat; p_final.GetSensView(0, num_states, s_dat); for (int i = 0; i < num_states; ++i) { n = p_final.states_(0, i); pdf = -t_final * exp(-lambda * t_final) * pow(lambda * t_final, double(n)) / tgamma(n + 1) + exp(-lambda * t_final) * t_final * pow(lambda * t_final, double(n - 1)) / tgamma(n); stmp += abs(s_dat[i] - pdf); } p_final.RestoreSensView(0, s_dat); MPI_Allreduce(&stmp, MPI_IN_PLACE, 1, MPIU_REAL, MPIU_SUM, PETSC_COMM_WORLD); ASSERT_LE(stmp, 1.0e-6); } /* * ===================================================================================================================== * TEST FIXTURE: Bursting-gene (i.e. telegraph) model with time-invariant propensities. */ namespace telegraph_cme { /* Stoichiometric matrix of the telegraph model */ arma::Mat SM{{-1, 1, 0, 0}, {1, -1, 0, 0}, {0, 0, 1, -1}}; /* Parameters for the propensity functions */ const double k01{1.0e-2}, k10{1.0e-1}, kr{10.0}, gamma{1.0}; // propensity function for telegraph model int propensity(const int reaction, const int num_species, const int num_states, const PetscInt *X, double *outputs, void *args) { int (*X_view)[3] = ( int (*)[3] ) X; switch (reaction) { case 0: for (int i{0}; i < num_states; ++i) { outputs[i] = k01 * X_view[i][0]; } break; case 1: for (int i{0}; i < num_states; ++i) { outputs[i] = k10 * X_view[i][1]; } break; case 2: for (int i{0}; i < num_states; ++i) { outputs[i] = kr * X_view[i][1]; } break; case 3: for (int i{0}; i < num_states; ++i) { outputs[i] = gamma * X_view[i][2]; } break; default:return -1; } return 0; } // Derivatives of the propensity functions of the telegraph model int propensity_derivatives( const int parameter_idx, const int reaction, const int num_species, const int num_states, const PetscInt *X, double *outputs, void *args) { int (*X_view)[3] = ( int (*)[3] ) X; switch (parameter_idx){ case 0: if (reaction == 0){ for (int i{0}; i < num_states; ++i) { outputs[i] = X_view[i][0]; } } break; case 1: if (reaction == 1){ for (int i{0}; i < num_states; ++i) { outputs[i] = X_view[i][1]; } } break; case 2: for (int i{0}; i < num_states; ++i) { outputs[i] = X_view[i][1]; } break; case 3: for (int i{0}; i < num_states; ++i) { outputs[i] = X_view[i][2]; } break; default: return -1; } return 0; } } class SensFspTelegraphTest : public ::testing::Test { protected: SensFspTelegraphTest() {} void SetUp() override { int n_par = 6; t_final = 100.0; fsp_tol = 1.0e-10; X0 = X0.t(); dp0 = std::vector>(n_par, arma::Col({0.0})); telegraph_model = SensModel(4, telegraph_cme::SM, {}, nullptr, telegraph_cme::propensity, nullptr, {}, telegraph_cme::propensity_derivatives, {{0}, {1}, {2}, {3}}); } void TearDown() override { } PetscReal t_final, fsp_tol; arma::Mat X0{1, 0, 0}; arma::Col p0 = {1.0}; std::vector> dp0; SensModel telegraph_model; arma::Row fsp_size = {2, 2, 1}; arma::Row expansion_factors = {0.25, 0.25, 0.25}; }; TEST_F(SensFspTelegraphTest, telegraph_sens_solve_with_cvode) { PetscInt ierr; PetscReal stmp; SensDiscreteDistribution p_final_bdf; std::vector p_snapshots_bdf; SensFspSolverMultiSinks fsp(PETSC_COMM_WORLD); ierr = fsp.SetModel(telegraph_model); ASSERT_FALSE(ierr); ierr = fsp.SetInitialBounds(fsp_size); ASSERT_FALSE(ierr); ierr = fsp.SetExpansionFactors(expansion_factors); ASSERT_FALSE(ierr); ierr = fsp.SetInitialDistribution(X0, p0, dp0); ASSERT_FALSE(ierr); ierr = fsp.SetUp(); ASSERT_FALSE(ierr); p_final_bdf = fsp.Solve(t_final, fsp_tol); fsp.ClearState(); ierr = VecSum(p_final_bdf.p_, &stmp); ASSERT_FALSE(ierr); ASSERT_GE(stmp, 1.0 - fsp_tol); for (int i{0}; i < telegraph_model.num_parameters_; ++i) { ierr = VecSum(p_final_bdf.dp_[i], &stmp); ASSERT_FALSE(ierr); ASSERT_LE(abs(stmp), 1.0e-6); } } \ No newline at end of file diff --git a/tests/test_sensmat.cpp b/tests/test_sensmat.cpp index f6801aa..a79b6b3 100644 --- a/tests/test_sensmat.cpp +++ b/tests/test_sensmat.cpp @@ -43,183 +43,188 @@ class SensMatrixTest : public ::testing::Test { SensMatrixTest() {} void SetUp() override { - fsp_size = arma::Row({12}); - t_fun = [&](double t, int num_coefs, double *outputs, void *args) { - outputs[0] = rate_right; - outputs[1] = rate_left; - return 0; - }; - - d_t_fun1 = [&](double t, int num_coefs, double *outputs, void *args) { - outputs[0] = 1.0; - return 0; - }; - - d_t_fun2 = [&](double t, int num_coefs, double *outputs, void *args) { - outputs[1] = 1.0; - return 0; - }; - - propensity = [&](const int reaction, const int num_species, const int num_states, const int *X, - double *outputs, void *args) { - switch (reaction) { - case 0: - for (int i{0}; i < num_states; ++i) { - outputs[i] = 1.0; - } - break; - case 1: - for (int i{0}; i < num_states; ++i) { - outputs[i] = (X[i] > 0); - } - break; - default:return -1; - } - return 0; - }; - - - int rank; - MPI_Comm_rank(PETSC_COMM_WORLD, &rank); - - int ierr; - // Read options for state_set_ - char opt[100]; - PetscBool opt_set; - PartitioningType fsp_par_type = PartitioningType::GRAPH; - ierr = PetscOptionsGetString(NULL, PETSC_NULL, "-fsp_partitioning_type", opt, 100, &opt_set); - ASSERT_FALSE(ierr); - - arma::Mat X0(1, 1); - X0.fill(0); - X0(0) = 0; - state_set = std::make_shared(PETSC_COMM_WORLD); - ierr = state_set->SetStoichiometryMatrix(stoichiometry); - ASSERT_FALSE(ierr); - ierr = state_set->SetShapeBounds(fsp_size); - ASSERT_FALSE(ierr); - ierr = state_set->SetUp(); - ASSERT_FALSE(ierr); - ierr = state_set->AddStates(X0); - ASSERT_FALSE(ierr); - ierr = state_set->Expand(); - ASSERT_FALSE(ierr); - - smodel = pacmensl::SensModel(stoichiometry, - std::vector({0,1}), - t_fun, - propensity, - std::vector({d_t_fun1, d_t_fun2}), - std::vector({propensity, propensity}), - std::vector({0,1}), - std::vector({0,1,2}) - ); + fsp_size = arma::Row({12}); + + t_fun = [&](double t, int num_coefs, double *outputs, void *args) { + outputs[0] = rate_right; + outputs[1] = rate_left; + return 0; + }; + + dt_fun = [&](int par_idx, double t, int num_coefs, double *outputs, void *args) { + switch (par_idx) { + case 0:outputs[0] = 1.0; + break; + case 1:outputs[1] = 1.0; + break; + default:break; + } + return 0; + }; + + std::vector> dt_sp = {{0}, {1}}; + + propensity = [&](const int reaction, const int num_species, const int num_states, const int *X, + double *outputs, void *args) { + switch (reaction) { + case 0: + for (int i{0}; i < num_states; ++i) { + outputs[i] = 1.0; + } + break; + case 1: + for (int i{0}; i < num_states; ++i) { + outputs[i] = (X[i] > 0); + } + break; + default:return -1; + } + return 0; + }; + + int rank; + MPI_Comm_rank(PETSC_COMM_WORLD, &rank); + + int ierr; + // Read options for state_set_ + char opt[100]; + PetscBool opt_set; + PartitioningType fsp_par_type = PartitioningType::GRAPH; + ierr = PetscOptionsGetString(NULL, PETSC_NULL, "-fsp_partitioning_type", opt, 100, &opt_set); + ASSERT_FALSE(ierr); + + arma::Mat X0(1, 1); + X0.fill(0); + X0(0) = 0; + state_set = std::make_shared(PETSC_COMM_WORLD); + ierr = state_set->SetStoichiometryMatrix(stoichiometry); + ASSERT_FALSE(ierr); + ierr = state_set->SetShapeBounds(fsp_size); + ASSERT_FALSE(ierr); + ierr = state_set->SetUp(); + ASSERT_FALSE(ierr); + ierr = state_set->AddStates(X0); + ASSERT_FALSE(ierr); + ierr = state_set->Expand(); + ASSERT_FALSE(ierr); + + smodel = pacmensl::SensModel(2, + stoichiometry, + std::vector({0, 1}), + t_fun, + propensity, + dt_fun, + dt_sp, + nullptr, + {} + ); }; void TearDown() override { }; std::shared_ptr state_set; - arma::Row fsp_size; - const double rate_right = 2.0, rate_left = 3.0; - const arma::Mat stoichiometry{1, -1}; - pacmensl::TcoefFun t_fun; - pacmensl::TcoefFun d_t_fun1, d_t_fun2; - pacmensl::PropFun propensity; + arma::Row fsp_size; + const double rate_right = 2.0, rate_left = 3.0; + const arma::Mat stoichiometry{1, -1}; + pacmensl::TcoefFun t_fun; + pacmensl::DTcoefFun dt_fun; + pacmensl::PropFun propensity; + pacmensl::DPropFun dpropensity; pacmensl::SensModel smodel; }; TEST_F(SensMatrixTest, mat_base_generation) { - int ierr; - - double Q_sum; - - pacmensl::SensFspMatrix A(PETSC_COMM_WORLD); - ierr = A.GenerateValues(*state_set, smodel); - ASSERT_FALSE(ierr); - - Vec P, Q; - ierr = VecCreate(PETSC_COMM_WORLD, &P); - ASSERT_FALSE(ierr); - ierr = VecSetSizes(P, state_set->GetNumLocalStates(), PETSC_DECIDE); - ASSERT_FALSE(ierr); - ierr = VecSetFromOptions(P); - ASSERT_FALSE(ierr); - ierr = VecSet(P, 1.0); - ASSERT_FALSE(ierr); - ierr = VecSetUp(P); - ASSERT_FALSE(ierr); - ierr = VecDuplicate(P, &Q); - ASSERT_FALSE(ierr); - ierr = VecSetUp(Q); - ASSERT_FALSE(ierr); - - ierr = A.Action(0.0, P, Q); - ASSERT_FALSE(ierr); - - ierr = VecSum(Q, &Q_sum); - ASSERT_FALSE(ierr); - ASSERT_DOUBLE_EQ(Q_sum, -1.0 * rate_right); - - for (int i_par{0}; i_par < 1; ++i_par) { - PetscReal dqsum = (i_par == 0) ? -1.0 : 0.0; - ierr = A.SensAction(i_par, 0.0, P, Q); + int ierr; + + double Q_sum; + + pacmensl::SensFspMatrix A(PETSC_COMM_WORLD); + ierr = A.GenerateValues(*state_set, smodel); ASSERT_FALSE(ierr); - ierr = VecSum(Q, &Q_sum); + + Vec P, Q; + ierr = VecCreate(PETSC_COMM_WORLD, &P); + ASSERT_FALSE(ierr); + ierr = VecSetSizes(P, state_set->GetNumLocalStates(), PETSC_DECIDE); + ASSERT_FALSE(ierr); + ierr = VecSetFromOptions(P); + ASSERT_FALSE(ierr); + ierr = VecSet(P, 1.0); + ASSERT_FALSE(ierr); + ierr = VecSetUp(P); + ASSERT_FALSE(ierr); + ierr = VecDuplicate(P, &Q); + ASSERT_FALSE(ierr); + ierr = VecSetUp(Q); + ASSERT_FALSE(ierr); + + ierr = A.Action(0.0, P, Q); ASSERT_FALSE(ierr); - ASSERT_DOUBLE_EQ(Q_sum, dqsum); - } - ierr = VecDestroy(&P); - ASSERT_FALSE(ierr); - ierr = VecDestroy(&Q); - ASSERT_FALSE(ierr); + ierr = VecSum(Q, &Q_sum); + ASSERT_FALSE(ierr); + ASSERT_DOUBLE_EQ(Q_sum, -1.0 * rate_right); + + for (int i_par{0}; i_par < 1; ++i_par) { + PetscReal dqsum = (i_par == 0) ? -1.0 : 0.0; + ierr = A.SensAction(i_par, 0.0, P, Q); + ASSERT_FALSE(ierr); + ierr = VecSum(Q, &Q_sum); + ASSERT_FALSE(ierr); + ASSERT_DOUBLE_EQ(Q_sum, dqsum); + } + + ierr = VecDestroy(&P); + ASSERT_FALSE(ierr); + ierr = VecDestroy(&Q); + ASSERT_FALSE(ierr); } TEST_F(SensMatrixTest, mat_constr_generation) { - int ierr; - - double Q_sum; + int ierr; - pacmensl::SensFspMatrix A(PETSC_COMM_WORLD); + double Q_sum; - ierr = A.GenerateValues(*state_set, smodel); - ASSERT_FALSE(ierr); + pacmensl::SensFspMatrix A(PETSC_COMM_WORLD); - Vec P, Q; - ierr = VecCreate(PETSC_COMM_WORLD, &P); - ASSERT_FALSE(ierr); - ierr = VecSetSizes(P, A.GetNumLocalRows(), PETSC_DECIDE); - ASSERT_FALSE(ierr); - ierr = VecSetFromOptions(P); - ASSERT_FALSE(ierr); - ierr = VecSet(P, 1.0); - ASSERT_FALSE(ierr); - ierr = VecSetUp(P); - ASSERT_FALSE(ierr); - ierr = VecDuplicate(P, &Q); - ASSERT_FALSE(ierr); - ierr = VecSetUp(Q); - ASSERT_FALSE(ierr); - - ierr = A.Action(0.0, P, Q); - ASSERT_FALSE(ierr); + ierr = A.GenerateValues(*state_set, smodel); + ASSERT_FALSE(ierr); - ierr = VecSum(Q, &Q_sum); - ASSERT_FALSE(ierr); - ASSERT_DOUBLE_EQ(Q_sum, 0.0); + Vec P, Q; + ierr = VecCreate(PETSC_COMM_WORLD, &P); + ASSERT_FALSE(ierr); + ierr = VecSetSizes(P, A.GetNumLocalRows(), PETSC_DECIDE); + ASSERT_FALSE(ierr); + ierr = VecSetFromOptions(P); + ASSERT_FALSE(ierr); + ierr = VecSet(P, 1.0); + ASSERT_FALSE(ierr); + ierr = VecSetUp(P); + ASSERT_FALSE(ierr); + ierr = VecDuplicate(P, &Q); + ASSERT_FALSE(ierr); + ierr = VecSetUp(Q); + ASSERT_FALSE(ierr); - for (int i_par{0}; i_par < 2; ++i_par) { - ierr = A.SensAction(i_par, 0.0, P, Q); + ierr = A.Action(0.0, P, Q); ASSERT_FALSE(ierr); + ierr = VecSum(Q, &Q_sum); ASSERT_FALSE(ierr); ASSERT_DOUBLE_EQ(Q_sum, 0.0); - } - ierr = VecDestroy(&P); - ASSERT_FALSE(ierr); - ierr = VecDestroy(&Q); - ASSERT_FALSE(ierr); + for (int i_par{0}; i_par < 2; ++i_par) { + ierr = A.SensAction(i_par, 0.0, P, Q); + ASSERT_FALSE(ierr); + ierr = VecSum(Q, &Q_sum); + ASSERT_FALSE(ierr); + ASSERT_DOUBLE_EQ(Q_sum, 0.0); + } + + ierr = VecDestroy(&P); + ASSERT_FALSE(ierr); + ierr = VecDestroy(&Q); + ASSERT_FALSE(ierr); }