#ifndef cvodeobj_h
#define cvodeobj_h

#include "nrnmpi.h"
#include "nrnneosm.h"
//#include "shared/nvector_serial.h"
#include "shared/nvector.h"
#include "membfunc.h"

// when USE_NVSERIAL is 0 then N_VNew_NrnThread is used even for single thread
// Then, (with USELONGDOUBLE in nvector_nrnthread.c), one should obtain
// quantitative double precision results independent of number of threads.
// Note that USE_NVSERIAL 1 should give double precision identity for
// the serial version and this version with nrn_nthread = 1
#define USE_NVSERIAL 1

class NetCvode;
class Daspk;
class TQItem;
class TQueue;
class PreSynList;
struct BAMech;
struct NrnThread;
class PlayRecList;
class PlayRecord;
class STEList;
class HTList;

class CvMembList {
public:
	CvMembList();
	virtual ~CvMembList();
	CvMembList* next;
	Memb_list* ml;
	int index;
};

class BAMechList {
public:
	BAMechList(BAMechList** first);
	BAMechList* next;
	BAMech* bam;
	Memb_list* ml;
	int* indices; // if nil then all in ml
	int cnt; // number of indices if non nil
	static void destruct(BAMechList** first);
	static void alloc(BAMechList* first);
};

#define CTD(i) ctd_[((nctd_ > 1) ? (i) : 0)]
class CvodeThreadData {
public:
	CvodeThreadData();
	virtual ~CvodeThreadData();
	void delete_memb_list(CvMembList*);

	int no_cap_count_; // number of nodes with no capacitance
	int no_cap_child_count_;
	Node** no_cap_node_;
	Node** no_cap_child_; // connected to nodes that have no capacitance
	CvMembList* cv_memb_list_;
	CvMembList* cmlcap_;
	CvMembList* cmlext_; // used only by daspk
	CvMembList* no_cap_memb_; // used only by cvode, point processes in the no cap nodes
	BAMechList* before_breakpoint_;
	BAMechList* after_solve_;
	BAMechList* before_step_;
	int rootnodecount_;
	int v_node_count_;
	Node** v_node_;
	Node** v_parent_;
	PreSynList* psl_th_; // with a threshold
	HTList* watch_list_;
	STEList* ste_list_;
	double** pv_;
	double** pvdot_;
	int nvoffset_;
	int nvsize_;
	int neq_v_; //for daspk
	PlayRecList* record_;
	PlayRecList* play_;
};

class Cvode {
public:
	Cvode(NetCvode*);
	Cvode();
	virtual ~Cvode();
	
	virtual int handle_step(NetCvode*, double);
	virtual int init(double t);
	virtual int advance_tn();
	virtual int interpolate(double t);
	virtual double tn() { return tn_;} // furthest time of advance
	virtual double t0() { return t0_;} // beginning of last real step
	void init_prepare();

	int solve(); // checks event_flag and init or advance_tn
	void statistics();
	double gam();
	double time() const { return t_; }
	void free_cvodemem();
	int order();
	void maxorder(int), minstep(double), maxstep(double);
public:
	double tn_, t0_, t_;
	boolean initialize_;
	boolean can_retreat_; // only true after an integration step
	// statistics
	void stat_init();
	int advance_calls_, interpolate_calls_, init_calls_;
	int f_calls_, mxb_calls_, jac_calls_, ts_inits_;
private:
	void alloc_cvode();
	void alloc_daspk();
	int cvode_init(double);
	int cvode_advance_tn();
	int cvode_interpolate(double);
	int daspk_init(double);
	int daspk_advance_tn();
	int daspk_interpolate(double);
public:
	N_Vector nvnew(long);
	int setup(N_Vector ypred, N_Vector fpred);
	int solvex_thread(double* b, double* y, NrnThread* nt);
	int solvex_thread_part1(double* b, NrnThread* nt);
	int solvex_thread_part2(NrnThread* nt);
	int solvex_thread_part3(double* b, NrnThread* nt);
	void fun_thread(double t, double* y, double* ydot, NrnThread* nt);
	void fun_thread_part1(double t, double* y, NrnThread* nt);
	void fun_thread_part2(NrnThread* nt);
	void fun_thread_part3(double* ydot, NrnThread* nt);
	boolean at_time(double, NrnThread*);
	void set_init_flag();
	void check_deliver(NrnThread* nt = 0);
	void evaluate_conditions(NrnThread* nt = 0);
	void ste_check();
	void states(double*);
	void dstates(double*);
	void error_weights(double*);
	void acor(double*);
	void fill(Cvode* standard);
	// following 7 crucial for local time step recording, also used by global
	void delete_prl();
	void record_add(PlayRecord*);
	void record_continuous();
	void record_continuous_thread(NrnThread*);
	void play_add(PlayRecord*);
	void play_continuous(double t);
	void play_continuous_thread(double t, NrnThread*);
	void do_ode(NrnThread*);
	void do_nonode(NrnThread* nt = 0);
	double* n_vector_data(N_Vector, int);
private:
	void cvode_constructor();
	boolean init_global();
	void init_eqn();
	void daspk_init_eqn();
	void matmeth();
	void nocap_v(NrnThread*);
	void nocap_v_part1(NrnThread*);
	void nocap_v_part2(NrnThread*);
	void nocap_v_part3(NrnThread*);
	void solvemem(NrnThread*);
	void atolvec_alloc(int);
	double h();
	N_Vector ewtvec();
	N_Vector acorvec();
	void new_no_cap_memb(CvodeThreadData&, NrnThread*);
	void before_after(BAMechList*, NrnThread*);
public:
	// daspk
	boolean use_daspk_;
	Daspk* daspk_;
	int res(double, double*, double*, double*, NrnThread*);
	int psol(double, double*, double*, double, NrnThread*);
	void daspk_scatter_y(N_Vector); // daspk solves vi,vx instead of vm,vx
	void daspk_gather_y(N_Vector);
	void daspk_scatter_y(double*, int);
	void daspk_gather_y(double*, int);
	void scatter_y(double*, int);
	void gather_y(N_Vector);
	void gather_y(double*, int);
	void scatter_ydot(double*, int);
	void gather_ydot(N_Vector);
	void gather_ydot(double*, int);
public:
	void activate_maxstate(boolean);
	void maxstate(double*);
	void maxstate(boolean, NrnThread* nt = 0);
	void maxacor(double*);
public:
	void* mem_;
	N_Vector y_;
	N_Vector atolnvec_;
	N_Vector maxstate_;
	N_Vector maxacor_;
public:
	boolean structure_change_;
#if USENEOSIM
	TQueue* neosim_self_events_;
#endif
public:
	CvodeThreadData* ctd_;
	NrnThread* nth_; // for lvardt
	int nctd_;
	long int* nthsizes_; // N_Vector_NrnThread uses this copy of ctd_[i].nvsize_
	NetCvode* ncv_;
	int neq_;
	int event_flag_;
	double next_at_time_;
	double tstop_;
	double tstop_begin_, tstop_end_;

private:
	void rhs(NrnThread*);
	void rhs_memb(CvMembList*, NrnThread*);
	void lhs(NrnThread*);
	void lhs_memb(CvMembList*, NrnThread*);
	void triang(NrnThread*);
	void bksub(NrnThread*);
private:
	// segregation of old vectorized information to per cell info
	friend class NetCvode;
	boolean is_owner(double*); // for play and record in local step context.
	boolean local_;
	void daspk_setup1_tree_matrix(); //unused
	void daspk_setup2_tree_matrix(); //unused
	TQItem* tqitem_;
private:
	int prior2init_;
#if PARANEURON
public:
	boolean use_partrans_;
	int global_neq_;
	int mpicomm_;
	int opmode_; // 1 advance, 2 interpolate, 3 init; for testing
#endif // PARANEURON
};

#endif