(* ::Package:: *)

(* ::Subsection::Closed:: *)
(*pre*)


(* :Title: ReactionKinetics *)

(* :Authors: Attila Laszlo Nagy, David Papp, Janos TOTH *)

(* :Summary:
	This is still a package intended to help doing fine reaction kinetics...
*)

(* :Discussion:
	...
*)

(* :Context: ReactionKinetics1` *)

(* :Package Version: 1.0 *)

(* :Copyright: Copyright 2011 by ... *)

(* :Acknowledgements: 
	Janos TOTH
*)

(* :Keywords:
    reaction kinetics, reaction decomposition, linear diophantine equations, stochastic simulation, Petri nets, stochastic Petri nets
*)

(* :History:
	Version 0.1: David Papp, May 2005.
*)

(* :Sources:
	Papp, D., Vizvari, B.: Effective solution of linear Diophantine equation systems
                            with an application in chemistry, RUTCOR Research Reports, RRR 28-2004, September 2004
*)

(* :Limitations:
	LPPartitions seems to work improperly sometimes due to bugs in LinearProgramming and related functions.
*)

(* :Warnings:
	* Needs Mathematica 6.0 or higher.
	* No other known restrictions, no known compatiblity problems!!!
 *)

(* :Mathematica Version: 8.0 *)

(*
   :Long-term `To-Do' list (beside the obvious and omnipresent minor upgrades needed)
	 1. Finish integrated Mathematica Help file
	 2. Waiting for a Mma version with bugfree Minimize (see MyMinimize)
*)


(*
	ReactionKinetics1.m is a Mathematica 8.0 package that helps you 
	to investigate, simulate, etc. (bio)chemical reactions and models..
    Copyright (C) 2010  Janos TOTH, David Papp, Attila Laszlo Nagy

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*)


(*
	''You acknowledge that Software is not designed, licensed or intended for
	use in the design, construction, operation or maintenance of any nuclear
	facility.'' (Java SDK demo applet)
*)


(* ::Subsection::Closed:: *)
(**begin*)


BeginPackage["ReactionKinetics1`",
	{"PlotLegends`",
	 "Units`",
	 "PhysicalConstants`",
	 "MultivariateStatistics`"
	 (*"MathSBML`"*)
	}]


If[$VersionNumber<6, 
	Off[DiscreteMath`Combinatorica`Star]; 
];

If[$VersionNumber<6, 
	Needs["Graphics`"];
	Needs["Graphics`Legend`"];
	Needs["LinearAlgebra`MatrixManipulation`"];
	Needs["Utilities`"];
]; 

If[$VersionNumber >= 6,
	Needs["XML`"];
	Needs["JLink`"];
	Needs["GUIKit`"];
];  (* XML` is loaded automatically in versions 4 & 5 *)

If[$VersionNumber>5,
	Needs["DifferentialEquations`NDSolveProblems`"];
	Needs["DifferentialEquations`NDSolveUtilities`"]; 
	Needs["FunctionApproximations`"];
	Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
	If[$VersionNumber<6, 
		Needs["DiscreteMath`GraphPlot`"];
	]
];

If[$VersionNumber > 7,
	Off[General::compat];
	Needs["Combinatorica`"];
	Needs["GraphUtilities`"];
	On[General::compat];
	Needs["EquationTrekker`"];
]; 

If[$VersionNumber<6, 
	On[DiscreteMath`Combinatorica`Star];
];


Unprotect[
	$ReactionKinetics1,
	$ReactionKineticsPackageLoaded,
	$ReactionKineticsVersionNumber,
	$Maybe,
	(**)
	Global`aExplodator,(*Explodator coeffs*)
	Global`bExplodator,(*Explodator coeffs*)
	Global`fOregonator,(*Explodator coeffs*)
	Models,(**)
	Reactions,
	BioModels,
	BioReactions,
	GetProblem,
	(**)
	Global`c,(*reaction rate equation: func symbol*)
	ComplexColours,(**)
	ExternalSpecies,
	FormattedOutput,
	GeneratedRateCoefficient,
	PlotFunction,
	StronglyConnectedComponentsColours,
	ToCanonicalForm,(**)
	ToReversible,
	ReactionsData,(**)
	ReversibleQ,
	WeaklyReversibleQ,
	ShowFHJGraph,
	ShowVolpertGraph,
	Atoms,(**)
	ToAtomMatrix,
	AtomConservingQ,
	AtomsQ,
	FromAtomMatrix,
	DetailedBalanced,
	(**)
	Global`t, (*time*)
	Global`k, (*reaction rate symbol*)
	Conditions, (**)
	Positivity,
	TimeLimit,
	RightHandSide, (**)
	DeterministicModel,
	Concentrations,
	StationaryPoints,
	MassConservation,
	EigenvaluesJacobianMatrix,
	AbsoluteRobustness, (**)
	(**)
	VolpertIndexing, (**)
	(**)
	Global`\[DoubleStruckCapitalE],(*symbol of expected value process*)
	Global`g,(*generating func symbol*)
	Global`\[CapitalPi],(*stat dist symbol*)
	Global`P,(*abs probability func symbol*)
	Global`x,(*state vector*)
	Global`z,(*gen func arg symbol*)
	CombinatorialMoments, (**)
	FixedStepSize,
	MaxIteration,
	Method, (* StiffNessSwitching, slow- fast variables *)
	Species,
	Verbose,
	Volume, (**)
	MassConservationQ, (**)
	ChemicalMasterEquation,
	StationaryProbabilityDensityEquation,
	ProbabilityGeneratingFunctionEquation,
	SolveProbabilityGeneratingFunctionEquation,
	ExpectedvalueProcess,
	DependencyGraph,
	ShowDependencyGraph,
	Simulation, (**)
	SimulationPlot,
	SimulationPlot2D,
	SimulationPlot3D,
	(**)
	OpenReactionKineticsPalette,
	OpenReactionKineticsNamesPalette,
	ReactionRatesNotebook
];


(*
	Circular,
	ContejeanDevie,
	CoveringDecompositionSet,
	DecompositionQ,
	Decompositions,
	ElementaryReactionQ,
	ElementaryReactions,
	FastSelection,
	Filter,
	FromStep,
	GeneralizedDecompositionQ,
	GreedySelection,
	Indexable,
	Kinetic,
	LPBased,
	MinimalDecompositions,
	NonZeros,
	ObjectiveFunction,
	Obligatory,
	Omittable,
	OptimalPartitions,
	OriginalSelection,
	Preprocess,
	PrintDecomposition,
	ReactionQ,
	SelectMinimalDecompositions,
	ToStep,
	VectorFromPoints
*)


(* ================================================================ *)
(* usage messages for the exported functions and the context itself *)
(* ================================================================ *)


$ReactionKinetics1::usage = 
"\!\(\*\nStyleBox[\"ReactionsKinetics1.m\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\) \
is a package that is no doubt the nicest package ever created to help you doing fine reaction kinetics.";

$ReactionKineticsVersionNumber::usage = "\!\(\*\nStyleBox[\"$ReactionKineticsVersion\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\) \
gives the version number and the date of the program you are using.";

$ReactionKineticsPackageLoaded::usage = "\!\(\*\nStyleBox[\"$ReactionKineticsPackageLoaded\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\) \
returns true if the current version of package ReactionKinetics is loaded.";


$ReactionKineticsVersionNumber = "1.0 [01-01-2011]";


$Maybe::usage = "$Maybe is a special symbol returned by certain functions when they cannot decide whether the answer is yes or no.";


Global`aExplodator::usage = "An unknown coefficient in the model Oregonator";
Global`bExplodator::usage = "An unknown coefficient in the model Oregonator";
Global`fOregonator::usage = "An unknown coefficient in the model Oregonator";
Global`t::usage = "Default symbol used in the kinetic differential equation for independent variable, partial differential equations, etc.";
Global`k::usage = "Default symbol of rate constants used in kinetic differential equation.";
Global`c::usage = "Default symbol for the unknown function of kinetic differential equation.";
Global`\[DoubleStruckCapitalE]::usage = "Default function symbol of the expected value process.";
Global`g::usage = "Default function symbol of the equation of probability generating function.";
Global`\[CapitalPi]::usage = "Default symbol of stationary probability density function.";
Global`P::usage = "Default symbol of the absolute probability function.";
Global`x::usage = "Default symbol of state.";
Global`z::usage = "Default argument symbol of the probability generating function.";


ExternalSpecies::usage = "ExternalSpecies is an option for several functions involving ReactionsData, ShowFHJGraph, Concentrations, Simulation, etc. \
The option describes ...";


(*options*)
CombinatorialMoments::usage = "...";
ComplexColours::usage = "...";
Conditions::usage = "...";
FormattedOutput::usage = "...";
GeneratedRateCoefficient::usage = "...";
PlotFunction::usage = "...";
Positivity::usage = "...";
Species::usage = "...";
StronglyConnectedComponentsColours::usage = "...";
TimeLimit::usage = "...";


(* unprotect any system functions for which definitions will be made *)
(* esetleg az All-hoz es Method-hoz kellene hozzaadni az itteni hasznalatat! *)
(* protected = Unprotect[] *)


(* ::Subsection::Closed:: *)
(**main*)


(* =============================================== *)
(* begin the private context (implementation part) *)
(* =============================================== *)

Begin["`Private`"]


(* ::Subsubsection::Closed:: *)
(**Auxiliary functions*)


(* ============================================================== *)
(* definition of auxiliary functions and local (static) variables *)
(* ============================================================== *)

solvedLPs = 0;
currentbest = Infinity;

ComponentwiseSmallerEqualQ[v1_?VectorQ, v2_?VectorQ] := NonPositive[Max[v1 - v2]];
MinimalQ[v_?VectorQ, vlist_List] :=
	If[Scan[If[ComponentwiseSmallerEqualQ[#, v], Return[False]] &, vlist] === False, False, True];

MyFloor[x_, prec_] := Floor[x/prec]*prec;
MyCeiling[x_, prec_] := Ceiling[x/prec]*prec;

MyMinimize[var_, conlist_, vars_] :=
	Module[{c, b, A, lp},
		(* EZ LENNE A NORMALIS MEGOLDAS *)
		(* First[Minimize[var, Thread[conlist >= 0], vars]] *)
		(* c = Normal[Last[CoefficientArrays[var, vars]]]; *)
		c = Coefficient[var, #] & /@ vars;
		{b, A} = Normal[CoefficientArrays[conlist, vars]];
		lp = LinearProgramming[c, A, -b, -Infinity, Method -> "Simplex"]; (* Automatic, Simplex, RevisedSimplex, InteriorPoint *)
		If[Head[lp] === LinearProgramming, "no solution", c.lp]
	];

MyMaximize[var_, conlist_, vars_] := -MyMinimize[-var, conlist, vars];

ZeroVectorQ[v_List] := Max[Abs[v]] === 0;


VectorFromPoints::usage = "...";
VectorFromPoints::badarg = "Illegal argument of function VectorFromPoints (should \
VectorFromPoints[len_Integer, pts_List, vals_List, < pad_ > ]).";

VectorFromPoints[len_Integer, pts_List, vals_List, pad_:0] :=
	Module[ {u = Table[pad, {len}]}, 
			(u[[ First[#] ]] = Last[#]) & /@ Transpose[{pts, vals}];
			 u
	];

VectorFromPoints[___] := (Message[VectorFromPoints::"badarg"]; Return[$Failed];)


(*
FromStep::usage = "FromStep[step, species] converts the reaction step to the standard vector notation. \
The components correspond to the species given as the second argument, the order of components is the
order the species are given.";
FromStep::badarg = "Illegal argument of function FromStep (should be FromStep[step_?ReactionQ, species_List]).";

FromStep[step_, species_] :=
	Subtract @@ (VectorFromPoints[Length[species],
					Sequence @@ Transpose[#[Flatten[List[#]] & /@ (step /. Plus -> List) /.
				  	coeff_. spec_?(MemberQ[species, #]&) :> {Position[species, spec][[1,1]], coeff}]]] & /@ {Last, First}
	);

FromStep[___] := (Message[FromStep::badarg]; Return[$Failed];)

ToStep::usage = "ToStep[vector, species] pretty prints the reaction assigned \
to vector using species as the name of species.";
ToStep::badarg = "Illegal argument of function ToStep (should be ToStep[vector_?VectorQ, species_List]).";

ToStep[vector_?ReactionQ, species_List] :=
	First[
    		Inner[	Rule,
			{Abs[#2[#1, #3@Positive]].#4},
			{#2[#1, #3@Negative].#4},
			List
		] & [vector, ReplaceAll, _?# -> 0 &, species]
	];

ToStep[___] := (Message[ToStep::"badarg"]; Return[$Failed];)
*)


(*
ReactionQ::usage = "ReactionQ[expr] gives True if expr is a reaction, and False \
otherwise.";
ReactionQ = MatchQ[#, {__Integer}] &;

DecompositionQ::usage = "DecompositionQ[expr] gives True if expr is a decomposition, \
and False otherwise.";
DecompositionQ = MatchQ[#, {__Integer?NonNegative}] &;

ElementaryReactionQ::usage = "ElementaryReactionQ[expr] gives True if expr is an \
elementary reaction (i.e. it starts from at most two reactant species), and False \
otherwise.";
ElementaryReactionQ = ReactionQ[#] && Total[Select[#, Negative]] >= -2 &;

GeneralizedDecompositionQ::usage = "GeneralizedDecompositionQ[expr] gives True if expr is a \
generalized decomposition, and False otherwise.";
GeneralizedDecompositionQ = MatchQ[#, { (_Integer?NonNegative | _Rational?NonNegative).. }] &;
*)


(* ::Subsection:: *)
(*reactionkinetic models*)


Models::usage = "\!\(\*\nStyleBox[\"Models\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\) \
is a list of names of widespread mechanisms of \
reaction kinetics such as the Lotka-Volterra model, the Brusselator, the Explodator etc.";

Models := 
{
	"Autocatalator",
	"Belousov-Zhabotinsky",
	"Briggs-Rauscher",
	"Brusselator",
	"Colquhoun",
	"Decomposition",
	"DeYoung-Keizer",
	"Dimerization",
	"Eigen",
	"\[CapitalEAcute]rdi-Ropolyi",
	"Explodator",
	"Feinberg-Horn I",
	"Feinberg-Horn II",
	"FKN mechanism",
	"Hudson-R\[ODoubleDot]ssler",
	"Hyver",
	"Inflow",
	"Ivanova",
	"Kliemann",
	"Leonard-Reichl",
	"Lotka",
	"Lotka-Volterra",
(*	"Lotka-Volterra Cyclic",*)
	"Michaelis-Menten",
	"Oregonator",
	"Outflow",
	"Robertson",
	"Schl\[ODoubleDot]gl",
	"Triangle reaction",
(*	"Triangle reaction Reversible",*)
	"Verhulst",
	"Volpert",
	"Wegscheider",
	"Willamowski-R\[ODoubleDot]ssler"
};


Reactions::usage = "\!\(\*\nStyleBox[\"Reactions\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\) \
is a list of widespread mechanisms given as lists of reversible and irreversible reaction steps \
of reaction kinetics such as the Lotka-Volterra model, the Brusselator, the Explodator etc.";

Reactions :=
Thread[
	Rule[Models,
		{
	(* Autocatalator *) 
{"A"->"X","X"->"Y","X"+2"Y"->3"Y","Y"->"P"}, 
(*externals: "A", "P"*)

	(* Belousov-Zhabotinsky *) 
{"X"+"Y"+"H"->2"V",
"Y"+"A"+2"H"->"X"+"V",
2"X"->"V",
1/2"X"+"A"+"H"->"X"+"Z",
"X"+"Z"->1/2"X",
"Z"+"M"->"Q",
"Z"+"V"->"Y",
"V"->"Y",
"X"->"0",
"Y"->"0",
"Z"->"0",
"V"->"0"}, 
(* Bazsa, Gy. (szerk.): Nem.lin.din., 1992. pp. 105--107.
	externals: "A" (BrO3-), "H" (H+), M (CH2 (COOH)2), "Q" (???), "0" (kiaramlas a reaktorbol)*)
	(*{ X->Tooltip[X,"Subscript[HBrO,2]"], Y->Tooltip[Y,"Superscript[Br,-]"], Z->Tooltip[Z,"Z"],
		V->Tooltip[V,"Superscript[Ce,4+]"], M->Tooltip[M,"Subscript[CH, 2](COOHSubscript[), 2]"],
		A->Tooltip[A,"\!\*SubsuperscriptBox[\(BrO\), \(3\), \(-\)]"], H->Tooltip[H,"Superscript[H,+]"]
	   }*)

	(* Briggs-Rauscher *) 
{"IO3 -"+2"H2O2"+"CH2(CO2H)2"+"H +"->"ICH(CO2H)2"+2"O2"+3"H2O",
"IO3 -"+2"H2O2"+"H +"->"HIO"+2"O2"+2"H2O",
"HIO"+"CH2(CO2H)2"->"ICH(CO2H)2"+"H2O",
"IO3 -"+"I -"+2"H +"->"HIO2"+"HIO",
"HIO2"+"I -"+"H +"->2"HIO",
"HIO"+"H2O2"->"I -"+"O2"+"H +"+"H2O",
"IO3 -"+"HIO2"+"H +"->2"IO2"+"H2O",
"IO2"+"Mn 2+"+"H2O"->"HIO2"+"Mn(OH) 2+",
(*!!!*)"Mn 2+"+"H2O2"->"Mn 2+"+"H2O"+"HO2",
2"HO2"->"H2O2"+"O2",
2"HIO2"->"IO3 -"+"HIO"+"H +",
"I -"+"HIO"+"H +"->"I2"+"H2O",
"I2"+"CH2(CO2H)2"->"ICH(CO2H)2"+"H +"+"I -"},
(*http://www.chem.leeds.ac.uk/delights/texts/expt_11.html*)

	(* Brusselator *) 
{"X"\[LeftArrow]"A","B"+"X"->"Y"+"D","Y"+2"X"->3"X","X"->"E"}, 
(*externals: "A","B","D","E"*)

	(* Colquhoun *)
{"AR"\[Equilibrium]"A"+"R"\[Equilibrium]"RA","A"+"AR"\[Equilibrium]"ARA"\[Equilibrium]"A"+"RA"},

	(* Decomposition *)
{"C" -> "A"+"B"},

	(*DeYoung-Keizer "*)
{"S(000)"\[Equilibrium]"S(001)"\[Equilibrium]"S(100)"\[Equilibrium]"S(010)","S(100)"\[Equilibrium]"S(101)"\[Equilibrium]"S(110)","S(001)"\[Equilibrium]"S(101)"\[Equilibrium]"S(011)",
"S(011)"\[Equilibrium]"S(010)"\[Equilibrium]"S(111)","S(110)"\[Equilibrium]"S(010)"\[Equilibrium]"S(111)","S(111)"\[Equilibrium]"S(101)"},

	(* Dimerization*)
{"A"+"B" -> "C"},

	(* Eigen *) 
{"X"->0,"X"->2"X"},

	(* \[CapitalEAcute]rdi-Ropolyi *)
{"R"+4"T"\[Equilibrium]"T4R1"\[Equilibrium]"T4R2"\[Equilibrium]"T4R3"},
	
	(* Explodator *) 
{"A"+"X"->(1+Global`aExplodator)"X","X"+"Y"->"Z","Z"->(1+Global`bExplodator)"Y","Y"->"P"},
(*externals: "A","P", aExplodator=bExplodator=1*)

	(* Feinberg-Horn I *)
{3"X"->"X"+2"Y"->3"Y"->2"X"+"Y"->3"X"},

	(* Feinberg-Horn II *) 
{"X"->"Y","X"+"Z"->"T"->"Y"+"U"->"X"+"Z"},

	(* FKN mechanism *)
{"HBrO2"+"Br -"->2"HOBr","Br -"+"BrO3 -"->"HBrO2"+"HOBr",2"HBrO2"->"BrO3 -"+"HOBr","HOBr"+"Br -"->"Br2"+"H2O","BrO3 -"+"HBrO2"->2"BrO2"+"H2O",
"HBrO2"+"H2O"->"BrO3 -","Br2"+"MA"->"BrMA"+"Br -","MA"+2"H2O"->C+2"CO2","BrMA"+2"H2O"->"Br -"+"HCO2H"+2"CO2"},

	(* Hudson-R\[ODoubleDot]ssler *) 
{"P"->"A","Q"->"B","A"+"B"->2"B","B"->"R","A"\[LeftRightArrow]"C"}, 
(*externals: "P","Q","R"*)

	(* Hyver *) 
{},

	(* Inflow *)
{"0"->"X"},
(*externals: "A"*)

	(* Ivanova *)
{"X"+"Y"->2"Y","Y"+"Z"->2"Z","X"+"Z"->2"X"},

	(* Kliemann *)
{"X"\[Equilibrium]"Y",2"X"\[Equilibrium]"X"+"Y"\[Equilibrium]2"Y"},

	(* Leonard-Reichl *) 
{2"X"\[LeftArrowRightArrow]"Y"},

	(* Lotka *) 
{"A"->"X","X"+"Y"->2"Y","Y"->"P"},

	(* Lotka-Volterra *)
{"A"+"X"->2"X","X"+"Y"->2"Y","B"\[LeftArrow]"Y"},
(*	(* Lotka-Volterra Cyclic *)
{"A"+"B"->2"A","B"+"C"->2"B","C"+"A"->2"C"},*)

	(* Michaelis-Menten *)
{"E"+"S"\[Equilibrium]"C"->"E"+"P"},

	(* Oregonator *) 
{"A"+"Y"->"X","X"+"Y"->"P","B"+"X"->2"X"+"Z",2"X"->"Q","Z"->Global`fOregonator*"Y"},

	(* Outflow *)
{"X"->"0"},
(*externals: "A"*)

	(* Robertson *)
{"A"->"B", 2"B"->"C"+"B", "B"+"C"->"A"+"C"},

	(* Schl\[ODoubleDot]gl *) 
{0\[Equilibrium]"X",2"X"\[Equilibrium]3"X"},

	(* Triangle reaction *)
{"A"->"B"\[RightArrow]"C"->"A"},
(*	(* Triangle reaction Reversible *)
{"A"\[ReverseEquilibrium]"B"\[LeftRightArrow]"C"\[RightArrowLeftArrow]"A"},
*)

	(* Verhulst *)
{"X"\[Equilibrium]2"X"},

	(* Volpert *)
{"X"+"Y"->"T","Y"+"Z"->"U"},

	(* Wegscheider *)
{"A"\[Equilibrium]"B",2"A"\[Equilibrium]"A"+"B"},

	(* Willamowski-R\[ODoubleDot]ssler *)
{"A"+"X"\[ReverseEquilibrium]2"X","X"+"Y"\[ReverseEquilibrium]2"Y","Q"+"Y"\[ReverseEquilibrium]"B","X"+"Z"\[ReverseEquilibrium]"C","P"+"Z"\[ReverseEquilibrium]2"Z"}
(*http://delfin.klte.hu/~gasparv/doktori/kaosz.html*)
		}
		]
];


BioModels::usage = "\!\(\*\nStyleBox[\"BioModels\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\) \
is a list of names of widespread used biochemical models such as Szent-Gy\[ODoubleDot]rgyi-Krebs cycle etc.";

BioModels := 
{
	"Calvin Cycle (Photosynthesis)",
	"Glycolysis",
	"Glyoxylate Cycle",
	"SzentGy\[ODoubleDot]rgyi-Krebs Cycle"
};


BioReactions::usage = "\!\(\*\nStyleBox[\"BioReactions\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\) \
is a list of widespread biochemical reactions given as lists of reversible and irreversible reaction steps \
such as the Calvin cycle or the Szent-Gy\[ODoubleDot]rgyi-Krebs cycle etc.";

BioReactions=
Thread[
	Rule[BioModels,
		{
	(* Calvin Cycle (Photosynthesis) *)
{"DG3P"->"GP","DF1,6BP"->"GP"+"DG3P","DF1,6BP"+"H2O"->"DF6P"+"Pi","S7P"+"DG3P"->"DR5P"+"DX5P","DR5P"->"DRL5P","DRL5P"->"DX5P",
"ATP"+"DRL5P"->"ADP"+"DRL1,5BP","3PDG"+2"H+"->"DRL1,5BP"+"CO2"+"H2O" ,"ATP"+"3PDG"->"ADP"+"3PDGP","DG3P"+"Pi"+"NADP+"->"3PDGP"+ "NADPH" + "H+"}/.
	{"DG3P"->Tooltip["DG3P","D-glyceraldehyde 3-phosphate"],
	"GP"->Tooltip["GP","glycerone phosphate"],
	"DRL5P"->Tooltip["DRL5P","D-ribulose 5-phosphate"],
	"DR5P"->Tooltip["DR5P","D-ribose 5-phosphate"],
	"DX5P"->Tooltip["DX5P","D-xylulose 5-phosphate"],
	"S7P"->Tooltip["S7P","sedoheptulose 7-phosphate"],
	"DF6P"->Tooltip["DF6P","D-fructose 1,6-bisphosphate"],
	"DG3P"->Tooltip["DG3P","D-glyceraldehyde 3-phosphate"],
	"DF1,6BP"->Tooltip["DF1,6BP","D-fructose 1,6-bisphosphate"],
	"3PDG"->Tooltip["3PDG","3-phospho-D-glycerate"],
	"Pi"->Tooltip["Pi","phosphate"],
	"DRL1,5BP"->Tooltip["DRL1,5BP","D-ribulose 1,5-bisphosphate"],
	"3PDGP"->Tooltip["3PDGP","3-phospho-D-glyceroyl phosphate"],
	"H2O"->Tooltip["H2O","dihydrogen monoxide"],
	"H+"->Tooltip["H+","hydrogen ion"],
	"CO2"->Tooltip["CO2","carbon dioxide"],
	"ATP"->Tooltip["ATP","adenosine-5'-triphosphate"],
	"ADP"->Tooltip["ADP","adenosine diphosphate"],
	"NADP+"->Tooltip["NADP+","nicotinamide adenine dinucleotide phosphate"],
	"NADPH"->Tooltip["NADPH","NADPH"](*???*)},

	(* Glycolysis *)
{"Glc"+"ATP"-> "G6P"+"ADP"+"H+","G6P"\[Equilibrium]"F6P","F6P"+"ATP"->"ADP"+"F1,6BP","F1,6BP"\[RightArrowLeftArrow]"GADP"+"DHAP","DHAP"\[Equilibrium]"GADP",
"GADP"+"NAD+"+"Pi"-> "NADH"+"H+"+"1,3BPG","1,3BPG"->"GADP","1,3BPG"+"ADP"->"ATP"+"3PG","3PG"\[ReverseEquilibrium]"2PG",
"2PG"->"PEP"+"H2O","PEP"->"2PG","PEP"+"ADP"-> "ATP"+"Pyr"}/.
	{"Glc"->Tooltip["Glc","D-glucose"],
	"G6P"->Tooltip["G6P","\[Alpha]-D-glucose-6-phosphate"],
	"F6P"->Tooltip["F6P","\[Beta]-D-fructose 6-phosphate"],
	"Pi"->Tooltip["Pi","phosphate"],
	"F1,6BP"->Tooltip["F1,6BP","\[Beta]-D-fructose 1,6-bisphosphate"],
	"GADP"->Tooltip["GADP","D-glyceraldehyde 3-phosphate"],
	"DHAP"->Tooltip["DHAP","dihydroxyacetone phosphate"],
	"1,3BPG"->Tooltip["1,3BPG","D-1,3-bisphosphoglycerate"],
	"3PG"->Tooltip["3PG","3-phosphoglycerate"],
	"2PG"->Tooltip["2PG","2-phosphoglycerate"],
	"PEP"->Tooltip["PEP","phosphoenolpyruvate"],
	"Pyr"->Tooltip["Pyr","pyruvate"],
	"H+"->Tooltip["H+","hydrogen ion"],
	"H2O"->Tooltip["H2O","dihydrogen monoxide"],
	"ATP"->Tooltip["ATP","adenosine-5'-triphosphate"],
	"ADP"->Tooltip["ADP","adenosine diphosphate"],
	"NAD+"->Tooltip["NAD+","nicotinamide adenine dinucleotide"],
	"NADH"->Tooltip["NADH","NADH"](*???*)},

	(* Glyoxylate Cycle *)
{"SM"+"NAD+"->"Oxc" + "NADH" + "H+","ACoA"+"H2O"+"Oxc"->"C"+"HS-CoA"+"H+","ACoA"+"H2O"+"G"->"SM"+"HS-CoA"+"H+","IsoC"->"S"+"G","C"->"cA"+"H2O","cA"+"H2O"->"IsoC"}/.
	{"Oxc"-> Tooltip["Oxc","oxalocetate"],
	"C"->Tooltip["C","citrate"],
	"ACoA"-> Tooltip["ACoA","acetyl CoA"],
	"SM"->Tooltip["SM","S-malate"],
	"S"-> Tooltip["S","succinate"],
	"IsoC"->Tooltip["IsoC","isocitrate"],
	"cA"->Tooltip["cA","cis-aconitate"],
	"G"->Tooltip["G","glyoxylate"],
	"H+"->Tooltip["H+","hydrogen ion"],
	"H2O"->Tooltip["H2O","dihydrogen monoxide"],
	"NAD+"->Tooltip["NAD+","nicotinamide adenine dinucleotide"],
	"HS-CoA"->Tooltip["HS-CoA","coenzyme A"],
	"NADH"->Tooltip["NADH","NADH"](*???*)},

	(* SzentGy\[ODoubleDot]rgyi-Krebs Cycle *)
{"Oc"+"ACoA"+"H2O"->"C"+"CSH","C"->"cA"+"H2O"->"IsoC","IsoC"+"NAD+"-> "Os"+"NADH"+"H+","Os"-> "\[Alpha]K"+"CO2","\[Alpha]K"+"NAD+"+"CSH"->"SCoA"+"NADH"+"H+"+"CO2",
"SCoA"+"GDP"+"Pi"->"S"+"CSH"+"GTP","S"+"Q"->"F"+"QH2","F"+"H2O"->"LM","LM"+"NAD+"->"Oc"+"NADH"+"H+"}/.
	{"Oc"->Tooltip["Oc","oxalocetate"],
	"ACoA"->Tooltip["ACoA","acetyl CoA"],
	"C"->Tooltip["C","citrate"],
	"CSH"->Tooltip["CSH","CoA-SH"],
	"cA"->Tooltip["cA","cis-aconitate"],
	"IsoC"->Tooltip["IsoC","isocitrate"],
	"Os"->Tooltip["Os","oxalosuccinate"],
	"\[Alpha]K"->Tooltip["\[Alpha]K","\[Alpha]-ketoglutarate"],
	"SCoA"->Tooltip["SCoA","succinyl-CoA"],
	"S"->Tooltip["S","succinate"],
	"Q"->Tooltip["Q","ubiquinone"],
	"F"->Tooltip["F","fumarate"],
	"QH2"-> Tooltip["QH2","ubiquinol"],
	"LM"->Tooltip["LM","L-malate"],
	"Pi"->Tooltip["Pi","phosphate"],
	"H2O"->Tooltip["H2O","dihydrogen monoxide"],
	"H+"->Tooltip["H+","hydrogen ion"],
	"CO2"->Tooltip["CO2","carbon dioxide"],
	"NAD+"->Tooltip["NAD+","nicotinamide adenine dinucleotide"],
	"GDP"->Tooltip["GDP","guanosine diphosphate"],
	"GTP"->Tooltip["GTP","guanosine-5'-triphosphate"],
	"NADH"->Tooltip["NADH","NADH"](*???*)}
		}
		]
];


ModelsSumPrivate := Join[BioModels, Models];

ReactionsSumPrivate := Join[BioReactions, Reactions];


GetProblem::usage=
"\!\(\*\nStyleBox[\"GetProblem\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\)\
\!\(\*\nStyleBox[\"[\",\nFontFamily->\"Courier\"]\)\
\"\!\(\*\nStyleBox[\"Models\",\nFontFamily->\"Courier\",\nFontSize->10]\)\"\
\!\(\*\nStyleBox[\"]\",\nFontFamily->\"Courier\"]\) gives all available built-in models.\n\
\!\(\*\nStyleBox[\"GetProblem\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\)\
\!\(\*\nStyleBox[\"[\",\nFontFamily->\"Courier\"]\)\
\!\(\*\nStyleBox[\"expr\",\nFontFamily->\"Courier\",\nFontSlant->\"Italic\",\nFontSize->10]\)\
\!\(\*\nStyleBox[\"]\",\nFontFamily->\"Courier\"]\) \
gives the mechanism with the name \!\(\*StyleBox[\"expr\",\nFontSlant->\"Italic\",\nFontSize->10]\), \
if \!\(\*\nStyleBox[\"expr\",\nFontSlant->\"Italic\",\nFontSize->10]\) \
is a string and is one of the built-in models, it fails otherwise.";

(*See for examples the list \!\(\*\nStyleBox[\" Models \",\n FontFamily->\" Courier \",\n FontColor->RGBColor[0, 0, 1]]\).*)


GetProblem::badarg = "Illegal argument of function GetProblem.";
GetProblem::nvmod = "Argument '`1`' is not a valid built-in (bio)model. For the existing models use \
\!\(\*\nStyleBox[\"GetProblem\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\)\
\!\(\*\nStyleBox[\"[\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 0]]\)\
\"\!\(\*\nStyleBox[\"Models\",\nFontFamily->\"Courier\",\nFontSize->10]\)\"\
\!\(\*\nStyleBox[\"]\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 0]]\).";

SyntaxInformation[GetProblem]={"ArgumentsPattern"->{__}};

GetProblem["Models"] := ModelsSumPrivate; (*Models*)

GetProblem[x_?StringQ] := 
	If[
		MemberQ[ModelsSumPrivate,x], (*Models*)
			x /. ReactionsSumPrivate, (*Reactions*)
			Message[GetProblem::"nvmod",x];
			$Failed
	  ];

SetAttributes[GetProblem,Listable];

GetProblem[x__?StringQ] := GetProblem[{x}];

GetProblem[___] := (Message[GetProblem::"badarg"]; $Failed)


(* ::Subsection:: *)
(**stoichiometry*)


(* ::Subsubsection::Closed:: *)
(**ReactionsData aux. funcs, ToCanonicalForm, ToReversible, ReactionsData, ReversibleQ, WeaklyReversibleQ*)


SpeciesQ := UpperCaseQ[StringTake[ToString[#],1]]&;

ExternalsQ := OptionQ[{#}]&&Union[First/@Flatten[{#}]]==={ExternalSpecies}&;

NotExternalsQ := Not[ExternalsQ[#]]&;

AllExternals := Flatten[#/.Rule->List/.ExternalSpecies->Sequence[]]&;


SymbolQ := Head[#]===Symbol&;

PlusQ := Head[#]===Plus&;

TimesQ := Head[#]===Times&;

SubscriptQ := Head[#]===Subscript&;

SuperscriptQ := Head[#]===Superscript&;

RightArrowQ := (Head[#]===RightArrow&&Length[#]===2&&
				And@@Map[StringQ[#]||SymbolQ[#]||PlusQ[#]||TimesQ[#]||SubscriptQ[#]||SuperscriptQ[#]&,#/.RightArrow->List])&;

LeftArrowQ := (Head[#]===LeftArrow&&Length[#]===2&&
				And@@Map[StringQ[#]||SymbolQ[#]||PlusQ[#]||TimesQ[#]||SubscriptQ[#]||SuperscriptQ[#]&,#/.LeftArrow->List])&;

EquilibriumQ := (Head[#]===Equilibrium&&Length[#]===2&&
				And@@Map[StringQ[#]||SymbolQ[#]||PlusQ[#]||TimesQ[#]||SubscriptQ[#]||SuperscriptQ[#]&,#/.Equilibrium->List])&;

ReverseEquilibriumQ := (Head[#]===ReverseEquilibrium && Length[#]===2 && 
						And@@Map[StringQ[#]||SymbolQ[#]||PlusQ[#]||TimesQ[#]||SubscriptQ[#]||SuperscriptQ[#]&,#/.ReverseEquilibrium->List])&;


ReactionsToList[{reactions__?RightArrowQ}]:= DeleteDuplicates[{reactions}];

ReactionsToList[{reactions__?LeftArrowQ}] := DeleteDuplicates[Reverse /@ {reactions} /. LeftArrow->RightArrow];

ReactionsToList[{reactions__?EquilibriumQ}]:=
		DeleteDuplicates[Sequence@@({#,Reverse[#]}/.Equilibrium->RightArrow)&/@{reactions}];

ReactionsToList[{reactions__?ReverseEquilibriumQ}]:=
		DeleteDuplicates[Sequence@@({Reverse[#],#}/.ReverseEquilibrium->RightArrow)&/@{reactions}];

ReactionsToList[{reactions__}]:=
	Module[{form, sp, a},
			form = DeleteDuplicates[ToCanonicalForm[reactions]];
			DeleteDuplicates[
			Flatten[
				Map[
					Switch[#,
						_?RightArrowQ,
							#,
						_?LeftArrowQ,
							Reverse[#] /. LeftArrow->RightArrow,
						_?EquilibriumQ,
							{#, Reverse[#]} /. Equilibrium->RightArrow,
						_?ReverseEquilibrium,
							{Reverse[#], #} /. ReverseEquilibrium->RightArrow,
						_,
							(sp = Flatten[ToBoxes[#,StandardForm]/.RowBox->List/.Thread[{"(",")","{","}",","}:>Sequence[]]];(*MakeBoxes*)
							a = Sort[Flatten[Position[sp,#]&/@{"\[Equilibrium]","\[RightArrow]","\[ReverseEquilibrium]","\[LeftArrow]"}]];
							(Partition[
									Riffle[
										ReleaseHold[
											MakeExpression[RowBox[Take[sp,#+{1,-1}]],StandardForm]
										] &/@ Partition[Join[{0},a,{Length[sp]+1}],2,1],
									Part[sp, a]],
							3,2] /. {
										{x_,"\[Equilibrium]",y_}->Sequence[RightArrow[x,y],RightArrow[y,x]],
										{x_,"\[ReverseEquilibrium]",y_}->Sequence[RightArrow[y,x],RightArrow[x,y]],
										{x_,"\[LeftArrow]",y_}->RightArrow[y,x],
										{x_,"\[RightArrow]",y_}->RightArrow[x,y]
									 }
							)
							)
					]&,
					form
				] /. RightArrow[A_,A_]->Sequence[]
			]
			]
];
(*DeleteDuplicates*)


Rdata[{reactions__},externals__] := Rdata[{reactions},externals] = 
	Module[
			{ reacs,
			  m, builtin, 
			  exs, coefflist, steps, steprule, complexes, vg, sp, indexed, vgraph,
			  lsp, lrs, e, fhjweak, fhjstrong, ab, nloc, lloc, sloc, gamma, aleph, a, b
			},
			reacs = Flatten[{reactions}];
			If[(m=Cases[reacs,_String])=!={},
				builtin = Check[(#->GetProblem[#])&/@ m, Return[$Failed];, GetProblem::"nvmod"];(*!*)
				Rdata[DeleteDuplicates[Join[Flatten[reacs/.builtin],reacs/.Thread[m->Sequence[]]]],externals]
				(*DeleteDuplicates*),
				exs = DeleteDuplicates[Prepend[AllExternals[{externals}],"0"]/.(0->"0")]; (*internalspecies megadas*)
				coefflist[cplxs_] := 
						Block[{cp = cplxs /. Thread[exs->0], cpa, spec},
							cpa = If[Head[cp]===Plus,
										Rest[Level[cp,1,Heads->True]],
										{cp}
								  ];
							spec = Union[((Flatten[cpa/.(0->{})] (*collect, symbolic coefficient, multip. coeff.*)
											/. r_*A_?StringQ :> A /. A_?StringQ*r_ :> A)
											/. r_*A_?SpeciesQ :> A /. A_?SpeciesQ*r_ :> A)
											/. A_?StringQ :> A /. A_?SpeciesQ :> A];
							Thread[List[spec,Coefficient[cplxs,spec]]]
						];
				steps = ReactionsToList[reacs] /.(0->"0"); (*"0" - zero complex*)
				steprule = steps /. RightArrow->Rule;
				complexes = Sort[VertexList[steprule]];
				vg = Flatten[(Thread /@ {#->coefflist[Last[#]],coefflist[First[#]]->#}) &/@ steps]
							/. Rule[A_RightArrow,B_]:>Prepend[Reverse[B],A->First[B]]
							/. Rule[A_,B_RightArrow]:>Prepend[Reverse[A],First[A]->B];
				(*Print[vg];*)
				sp = Union[Last /@ vg]; (*Union[Map[# /. {A_,r_}:>A &, coefflist[complexes]]];*)
				vgraph = Most /@ vg;
				indexed = Flatten[MapIndexed[#1->First[#2]&,#]&/@{sp,steps}];
				{lsp, lrs} = Length /@ {sp, steps};
				a = SparseArray[Cases[vgraph,{Rule[x_,y_RightArrow],z_}:> Rule[{x,y}/.indexed,z]],{lsp,lrs}];
				b = SparseArray[Cases[vgraph,{Rule[x_RightArrow,y_],z_}:> Rule[{y,x}/.indexed,z]],{lsp,lrs}];
				(*Print[{a,b}];*)
				Dispatch[{
					"species" -> sp,
					(*"reactantspecies" -> ...,
					"productspecies" -> ...,*)
					"M" -> lsp,
					"externalspecies" -> (e = If[MemberQ[complexes,"0"], exs, Rest[exs]]),
					"E" -> Length[e],
					"reactionsteps" -> steps,
					"R" -> lrs,
					"complexes" -> complexes,
					(*"reactantcomplexes" -> ...,
					"productcomplexes" -> ...,*)
					"fhjgraph" -> FromOrderedPairs[steps /. RightArrow[c1_,c2_] :> {Position[complexes, c1][[1, 1]], Position[complexes, c2][[1, 1]]}],
					"fhjconnectedcomponents" -> (fhjweak = WeakComponents[steprule]),
					"fhjstronglyconnectedcomponents" -> (fhjstrong = StrongComponents[steprule]),
					"fhjterminalstronglyconnectedcomponents"->
						Map[Complement[Union@@Map[Cases[steps,RightArrow[#,y_]:>y]&,#],#]==={}/.True->#/.False->Sequence[]&,fhjstrong],
					"n" -> (nloc = Length[complexes]),
					"l" -> (lloc = Length[fhjweak]),
					"s" -> (sloc = MatrixRank[gamma = First[Differences[ab = (Normal/@{a,b})]]]),
					(*MatrixRank[First[Differences[ab=Transpose /@ Transpose[Outer[Coefficient,steprule /. Rule->List,sp]]]]]*)
					(*option for SPARSEARRAY; for symbolic coefficient!*)
					"deficiency" -> "\[Delta]=n-l-s="<>ToString[nloc]<>"-"<>ToString[lloc]<>"-"<>ToString[sloc]<>"="<>ToString[nloc-lloc-sloc],
					"\[Alpha]" -> (aleph = Part[ab,1]),
					"\[Alpha] sparsearray" -> a,
					"\[Beta]" -> Part[ab,2],
					"\[Beta] sparsearray" -> b,
					"\[Gamma]" -> gamma,
					(*"\[Gamma] sparsearray" -> ...*)
					"reactionorder" -> Max[Total /@ Transpose[aleph]],
					"variables" -> (Subscript[Global`c, #] & /@ sp),
					"volpertgraph" -> vgraph	(* /. {A_,1}->A *)
					(*Flatten[vg /. Rule[A_RightArrow,B_List]:>{A->First[B],Last[B]} /. Rule[A_List,B_RightArrow]:>{First[A]->B,Last[A]},2]*)
				}]
			]
	];


ToCanonicalForm::usage = "\!\(\*\nStyleBox[\"ToCanonicalForm\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\)\
\!\(\*\nStyleBox[\"[\",\nFontFamily->\"Courier\"]\)reactions\!\(\*\nStyleBox[\"]\",\nFontFamily->\"Courier\"]\)\n\
	gives the canonical form of a given reaction network involving only the rightarrow, leftarrow and equilibrium.";

ToCanonicalForm::badarg = "Illegal argument of function ToCanonicalForm.";


SyntaxInformation[ToCanonicalForm]={"ArgumentsPattern"->{__}};

ToCanonicalForm[reactions__]:=
	Flatten[{reactions} /. Reactions] /. 
					  Join[Thread[{LeftRightArrow,DoubleLeftRightArrow,LongLeftRightArrow,DoubleLongLeftRightArrow,RightArrowLeftArrow}->Equilibrium],
									{LeftArrowRightArrow->ReverseEquilibrium},
									Thread[{Rule,ShortRightArrow,DoubleRightArrow,LongRightArrow,DoubleLongRightArrow}->RightArrow],
									Thread[{ShortLeftArrow,DoubleLeftArrow,LongLeftArrow,DoubleLongLeftArrow}->LeftArrow]];

ToCanonicalForm[___] := (Message[ToCanonicalForm::"badarg"]; $Failed)


ToReversible::usage = "\!\(\*\nStyleBox[\"ToReversible\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\)\
\!\(\*\nStyleBox[\"[\",\nFontFamily->\"Courier\"]\)reactions\!\(\*\nStyleBox[\"]\",\nFontFamily->\"Courier\"]\)\n\
	gives the reaction network associated to reactions in which every reaction step is reversible.";

ToReversible::badarg = "Illegal argument of function ToReversible.";


SyntaxInformation[ToReversible]={"ArgumentsPattern"->{__}};

ToReversible[reactions__] := ToCanonicalForm[reactions]/.RightArrow->Equilibrium/.LeftArrow->Equilibrium;

ToReversible[___] := (Message[ToReversible::"badarg"]; $Failed)


PropertiesReactionsData = 
	{
		"species","M","externalspecies","E","reactionsteps","R","complexes",
		"fhjgraph","fhjconnectedcomponents","fhjstronglyconnectedcomponents","fhjterminalstronglyconnectedcomponents",
		"n","l","s","deficiency","\[Alpha]","\[Alpha] sparsearray","\[Beta]","\[Beta] sparsearray","\[Gamma]","reactionorder","variables","volpertgraph"
	};


ReactionsForm[reactionlist_] := 
	If[Length[Flatten[reactionlist]]<=2,ToString[reactionlist,StandardForm],"{"<>ToString[First[reactionlist],StandardForm]<>",...,"<>ToString[Last[reactionlist],StandardForm]<>"}"];


ReactionsData::usage = "\!\(\*\nStyleBox[\"ReactionsData\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\)\
\!\(\*\nStyleBox[\"[\",\nFontFamily->\"Courier\"]\){reactions},externalspecies\!\(\*\nStyleBox[\"]\",\nFontFamily->\"Courier\"]\) \
yields a list of the most important characteristics of the chemical mechanism given by the list of \
reactions. With optional argument externalspecies, one can treats the species \
given as the function's second argument as external species having constant \
concentration.\n \
\!\(\*\nStyleBox[\"ReactionsData\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\)\
\!\(\*\nStyleBox[\"[\",\nFontFamily->\"Courier\"]\){reactions},externalspecies\!\(\*\nStyleBox[\"]\",\nFontFamily->\"Courier\"]\)\
\!\(\*\nStyleBox[\"[\",\nFontFamily->\"Courier\"]\)data\!\(\*\nStyleBox[\"]\",\nFontFamily->\"Courier\"]\)\n\
\t Data name: \n\t species:\n\t M:.\n\t externalspecies:\n\t E:.\n\t reactionsteps:\n\t R:...\n\t complexes:\n\t \
fhjgraph:...\n\t fhjconnectedcomponents:..\n\t fhjstronglyconnectedcomponents:.\n\t fhjterminalstronglyconnectedcomponents:..\n\t \
n\n\t l\n\t s\n\t deficiency:..\n\t \[Alpha]:..\n\t \[Alpha] sparsearray:..\n\t \[Beta]:..\n\t \[Beta] sparsearray:..\n\t \[Gamma]:...\n\t reactionorder:..\n\t variables:..\n\t \
volpertgraph:..\n\
Or Use Usage...";


ReactionsData::badarg = "Illegal argument of function ReactionsData.";
ReactionsData::badname = "At least one of the arguments '`1`' is a non-identified property. For the possibilities use the command \
\!\(\*\nStyleBox[\"ReactionsData\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\)[\"Properties\"].";
ReactionsData::wrreac = "Argument `1` may not have a correct form. The list of reactions must contain the operations Plus, Times, String, Symbol and many kinds of arrows only\ 
in a reasonable form. For help check the palette \!\(\*\nStyleBox[\"OpenReactionKineticsPalette\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\).";


SyntaxInformation[ReactionsData]={"ArgumentsPattern"->{__,OptionsPattern[]}};

Options[ReactionsData] = {ExternalSpecies->{}};

ReactionsData["Properties"] := PropertiesReactionsData;

ReactionsData[{reactions__},externals:(_?VectorQ|OptionsPattern[])][] := 
	Check[Rdata[{reactions},{externals}],
					Message[ReactionsData::"wrreac",ReactionsForm[{reactions}]];
					Return[$Failed]
	];
(*
ReactionsData[Longest[reactions__?NotExternalsQ],Shortest[externals:(___?ExternalsQ)]][] := 
	Check[Rdata[{reactions},{externals}],
					Message[ReactionsData::"wrreac",ReactionsForm[{reactions}]];
					Return[$Failed]
	];
*)
ReactionsData[{reactions__},externals:(_?VectorQ|OptionsPattern[])]["Properties"] := PropertiesReactionsData;
(*ReactionsData[Longest[reactions__?NotExternalsQ],Shortest[externals:(___?ExternalsQ)]]["Properties"] := PropertiesReactionsData;*)

ReactionsData[{reactions__},externals:(_?VectorQ|OptionsPattern[])]["All"] := ReactionsData[{reactions},externals][PropertiesReactionsData];
(*ReactionsData[Longest[reactions__?NotExternalsQ],Shortest[externals:(___?ExternalsQ)]]["All"] := ReactionsData[{reactions},externals][PropertiesReactionsData];*)

ReactionsData[{reactions__},externals:(_?VectorQ|OptionsPattern[])][data__?StringQ] := 
	(
	If[Nand @@ Map[MemberQ[PropertiesReactionsData,#]&,{data}],
		Message[ReactionsData::"badname",data];
	];
	Catch[
		Part[
			ReplaceAll[{data},
					Check[Rdata[{reactions},{externals}],
							Throw[
								Message[ReactionsData::"wrreac",ReactionsForm[{reactions}]];
								Return[$Failed]
							]
					]
			], Length[{data}] /. x_/;x>1 :> (1;;x)
		]
	]
	);

ReactionsData[{reactions__},externals:(_?VectorQ|OptionsPattern[])][data__] := 
	ReactionsData[{reactions},externals][Sequence @@ (ToString /@ Flatten[{data}])];
(*{externals}!*)

(*
ReactionsData[Longest[reactions__?NotExternalsQ],Shortest[externals:(___?ExternalsQ)]][data__?StringQ] := ReactionsData[{reactions},externals][data];
ReactionsData[Longest[reactions__?NotExternalsQ],Shortest[externals:(___?ExternalsQ)]][{data__?StringQ}] := ReactionsData[{reactions},externals][data];

FormatString[{reactions__},{externals___}] := 
	"The \!\(\*\nStyleBox[\"ReactionsData\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\) object is generated for the reaction(s)/model(s) "<>
	ReactionsForm[{reactions}] <> " with external species " <> ToString[AllExternals[{externals}],StandardForm] <> ", beside 0 or \"0\".\n" <> 
	"For more help check the function \!\(\*\nStyleBox[\"ReactionsData\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\), "<>
	"the palette of names: \!\(\*\nStyleBox[\"OpenReactionKineticsNamesPalette\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\) "<>
	"or the list of properties using \!\(\*\nStyleBox[\"ReactionsData\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\)[\"Properties\"].";
*)
Format[ReactionsData[{reactions__},externals:(_?VectorQ|OptionsPattern[])]] := 
	DynamicModule[{x, rd},
					x = Dynamic[Refresh[Round[Clock[Infinity]],UpdateInterval->1]];
					Monitor[rd = ReactionsData[{reactions},externals][];,
									Column[{Row[{"ReactionsData is now calculating ",Dynamic[Mod[First[x],5]/.{0->"",1->".",2->"..",3->"...",4->"...."}]}],
											ProgressIndicator[Dynamic[Clock[Infinity]],Indeterminate,ImageMargins->1,BaselinePosition->Center],
											Row[{x,". ","seconds passed"}]
									}]	
					];
					FinishDynamic[];
					rd(*FormatString[{reactions},{externals}]*)
	];
(*
Format[ReactionsData[Longest[reactions__?NotExternalsQ],Shortest[externals:(___?ExternalsQ)]]] := 
	(ReactionsData[{reactions},externals][];
	FormatString[{reactions},{externals}])
*)
ReactionsData[___][___] := (Message[ReactionsData::"badarg"]; Return[$Failed])


ReversibleQ::usage = "ReversibleQ[reactions] yields True if and only if the given formal \
chemical mechanism is reversible.";


ReversibleQ::badarg = "Illegal argument of function ReversibleQ.";


SyntaxInformation[ReversibleQ]={"ArgumentsPattern"->{__}};

ReversibleQ[reactions__] := 
		SymmetricQ[Check[ReactionsData[{reactions}]["fhjgraph"],Abort[];,{ReactionsData::wrreac,ReactionsData::badarg}]];
(*Abort[];*)

ReversibleQ[___] := (Message[ReversibleQ::"badarg"]; $Failed)


WeaklyReversibleQ::usage = "WeaklyReversibleQ[reactions] yields True if and only if the formal \
chemical mechanism reactions is weakly reversible.";


WeaklyReversibleQ::badarg = "Illegal argument of function WeaklyReversibleQ.";


SyntaxInformation[WeaklyReversibleQ]={"ArgumentsPattern"->{__}};

WeaklyReversibleQ[reactions__] := 
		SymmetricQ[TransitiveClosure[Check[ReactionsData[{reactions}]["fhjgraph"],Return[$Failed],{ReactionsData::wrreac,ReactionsData::badarg}]]];

WeaklyReversibleQ[___] := (Message[WeaklyReversibleQ::"badarg"]; $Failed)


(* ::Subsubsection::Closed:: *)
(*ShowFHJGraph, ShowVolpertGraph*)


(*
FHJComponents::usage = "FHJComponents[{reactions}] gives the strongly connected components of the Feinberg-Horn-Jackson \
graph of the formal chemical mechanism reactions.";

FHJStrongComponents::usage = "FHJStrongComponents[{reactions}] gives the connected components of the (undirected) \
Feinberg-Horn-Jackson graph of the formal chemical mechanism reactions.";
*)


GraphPlotFunctionQ := MemberQ[{"GraphPlot","LayeredGraphPlot","GraphPlot3D","TreePlot"},#]&;


ShowFHJGraph::usage = "ShowFHJGraph[{reactions},rates,options] visualizes the Feinberg-Horn-Jackson graph of the mechanism given as steps in \
reactions accepting also a set of external species which are to be deleted and it allows the use of options of GraphPlot. Options...";


ShowFHJGraph::badarg = "Illegal argument of function ShowFHJGraph.";
ShowFHJGraph::funcarg = "Argument `1` is not a valid graph plot function. Try GraphPlot, LayeredGraphPlot, GraphPlot3D or TreePlot.";
ShowFHJGraph::ccols = "The given colours may have not in the correct vector form or have wrong size.";
ShowFHJGraph::srates = "The given list of rate constants may have wrong size.";
ShowFHJGraph::sccomps = "The given list of colours for the strongly connected components may have wrong form or length or it is also given complex colours.";

Options[ShowFHJGraph]:={ComplexColours->{}, PlotFunction->"GraphPlot", StronglyConnectedComponentsColours->{}};

(*SyntaxInformation[ShowFHJGraph]={"ArgumentsPattern"->{{__},___}};*)

ShowFHJGraph[{reactions__}, rates_List:{}, opts___?OptionQ] := 
	Module[{ cc, n, r, rsteps, cxs, fhjscc, nopt, graphfunc, c, sc, rscc, rs },
			cc = FilterRules[{opts},Options[ShowFHJGraph]];
			{n, r, rsteps, cxs, fhjscc} = 
					Check[ReactionsData[{reactions}]["n","R","reactionsteps","complexes","fhjstronglyconnectedcomponents"] /. RightArrow->Rule, Return[$Failed], {ReactionsData::wrreac,ReactionsData::badarg}];
			graphfunc = (PlotFunction /. cc) /. PlotFunction->"GraphPlot";
			If[GraphPlotFunctionQ[graphfunc],
				graphfunc = ToExpression[graphfunc],
				Message[ShowFHJGraph::"funcarg",graphfunc];
				Return[$Failed](*Abort[];*)				
			];
			nopt = FilterRules[{opts},Options[graphfunc]];
			c = ComplexColours /. cc;
			sc = StronglyConnectedComponentsColours /. cc;
			If[sc=!={} && sc =!= StronglyConnectedComponentsColours,
				If[(VectorQ[sc]&&Length[fhjscc] === Length[sc])&&(c === {} || c === ComplexColours),
					rscc = rsteps /. Flatten[Thread[First[#]->Thread[Style[First[#],Italic,Last[#]]]] &/@ Thread[{fhjscc,sc}]];
						If[rates==={},
							graphfunc[rscc,nopt],
							If[r===Length[rates],
								graphfunc[Thread[{rscc,rates}],nopt],
								Message[ShowFHJGraph::"srates"];
								$Failed
							]
						],
					Message[ShowFHJGraph::"sccomps"];
					$Failed
				],
				If[c=!={} && c =!= ComplexColours,
					If[VectorQ[c] && n===Length[c],
						rs = rsteps /. (First[#]->Style[First[#],Italic,Last[#]] &/@ Thread[{cxs,c}]);
						If[rates==={},
							graphfunc[rs,nopt],
							If[r===Length[rates],
								graphfunc[Thread[{rs,rates}],nopt],
								Message[ShowFHJGraph::"srates"];
								$Failed
							]
						],
						Message[ShowFHJGraph::"ccols"];
						$Failed
					],
					If[rates==={},
						GraphPlot[rsteps,nopt],
						If[r===Length[rates],
								graphfunc[Thread[{rsteps,rates}],nopt],
								Message[ShowFHJGraph::"srates"];
								$Failed
							]
					]
				]
			]
	];

ShowFHJGraph[___] := (Message[ShowFHJGraph::"badarg"]; $Failed)


(*FHJ - graph $Failed -> Return[$Failed]!!!!*)


GraphPlotFunction2Q := MemberQ[{"GraphPlot","LayeredGraphPlot","GraphPlot3D","TreePlot","ShowGraph","ShowLabeledGraph"},#]&;


ShowVolpertGraph::usage = "ShowVolpertGraph[{reactions},options] visualizes the Volpert graph of the \
formal chemical mechanism reactions. Options...";


Options[ShowVolpertGraph]:={ExternalSpecies->{}, PlotFunction->"GraphPlot"};

ShowVolpertGraph::badarg = "Illegal argument of function ShowVolpertGraph.";
ShowVolpertGraph::funcarg = "Argument `1` is not a valid graph plot function. Try GraphPlot, LayeredGraphPlot, GraphPlot3D, TreePlot, ShowGraph or ShowLabeledGraph.";

(*SyntaxInformation[ShowFHJGraph]={"ArgumentsPattern"->{{__},___}};*)

ShowVolpertGraph[{reactions__},opts___?OptionQ] := 
	Module[{graphfunc, vg, nopts},
		graphfunc = (PlotFunction /. FilterRules[{opts},PlotFunction])/.PlotFunction->"GraphPlot";
		If[GraphPlotFunction2Q[graphfunc],
			graphfunc = ToExpression[graphfunc];
			vg = Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["volpertgraph"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];
			nopts = FilterRules[{opts},Options[graphfunc]];
			Switch[graphfunc,
					ShowGraph,ShowGraph[ToCombinatoricaGraph[vg],nopts],
					ShowLabeledGraph,ShowLabeledGraph[ToCombinatoricaGraph[vg],nopts],
					_,graphfunc[vg, nopts]
			],
			Message[ShowVolpertGraph::"funcarg",graphfunc];
			Return[$Failed]
		]
	];

ShowVolpertGraph[___] := (Message[ShowVolpertGraph::"badarg"]; $Failed)


(* ::Subsubsection::Closed:: *)
(*Atoms, ToAtomMatrix, FromAtomMatrix...*)


Atoms::usage = "\!\(\*\nStyleBox[\"Atoms\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\) \
is the list of current names of currently known atoms from Ac through Uuh to Y.";

Atoms = 
	{
	  "Ubn","Ubu","Uub","Uue","Uuh","Uuo","Uup","Uuq","Uus","Uut",
	  "Ac","Ag","Al","Am","Ar","As","At","Au","Ba","Be","Bh","Bi","Bk","Br","Ca","Cd","Ce","Cf","Cl","Cm","Co","Cr","Cs","Cu",
	  "Db","Ds","Dy","Er","Es","Eu","Fe","Fm","Fr","Ga","Gd","Ge","He","Hf","Hg","Ho","Hs","In","Ir","Kr","La","Li","Lr","Lu",
	  "Md","Mg","Mn","Mo","Mt","Na","Nb","Nd","Ne","Ni","No","Np","Os","Pa","Pb","Pd","Pm","Po","Pr","Pt","Pu","Ra","Rb","Re","Rf","Rg","Rh","Rn","Ru",
	  "Sb","Sc","Se","Sg","Si","Sm","Sn","Sr","Ta","Tb","Tc","Te","Th","Ti","Tl","Tm","Xe","Yb","Zn","Zr",
	  "B","C","F","H","I","K","N","O","P","S","U","V","W","Y"
	};
(* sorrend! *)


NestedMoleculeQ[s_?StringQ] := StringFreeQ[NestWhile[StringReplace[#,("("~~x___~~")"/;StringCount[x,{"(",")"}]===0)->""]&,s,UnsameQ[##]&,2],{"(",")"}];

MoleculeToAtoms[s_?StringQ,ats_?VectorQ]:=
	StringCases[s, 
			Join[" "~~r1___?DigitQ~~#:>ToExpression[#<>If[r1=="","1",r1]]*"charge"&/@{"+","-"},
			#~~r2___?DigitQ:>ToExpression[If[r2 == "","1",r2]]*#&/@ats,
			{"("~~x___?NestedMoleculeQ~~")"~~r3__?DigitQ:>Sequence@@Times[ToExpression[If[r3=="","1",r3]],MoleculeToAtoms[x,ats]]}],
				Overlaps->False];
(*tobbszorosen beagyazott molekulakra*)


ToAtomMatrix::usage = "ToAtomMatrix[molecules,options] constructs the atommatrix using the molecular structure of the species. Options...";
(*Checks if the atoms are conserved in the reaction steps and list those steps where they are not*)

ToAtomMatrix::badarg = "Illegal argument of function ToAtomMatrix.";

SyntaxInformation[ToAtomMatrix]={"ArgumentsPattern"->{_,OptionsPattern[]}};

ToAtomMatrix[molecules_?VectorQ,OptionsPattern[{FormattedOutput->False}]]:=
	Module[{ ats, l, toa },
		ats := Append[Reverse[Sort[Union@@StringCases[molecules,Atoms],StringLength[#1]<StringLength[#2]&]],"charge"];
		l := Length[ats];
		toa := {Map[Total[MoleculeToAtoms[#,ats]]&,molecules]/.Thread[ats->(UnitVector[l,#]&/@Range[l])],ats};
		If[OptionValue[FormattedOutput],
			Column[{Row[Riffle[Most[Last[toa]]," "]],
					Grid[MapThread[Prepend,{Prepend@@toa,Prepend[molecules,"Molecules\[Backslash]Atoms"]},1],Dividers->{{False,True},{False,True}}]}],
			toa
		]
	];

ToAtomMatrix[___] := (Message[ToAtomMatrix::"badarg"]; $Failed)


AtomConservingQ::usage = "AtomConservingQ[{reactions},externalspecies] gives True if the number of atoms are conserved in all \
the steps of reactions given in terms of species with atomic structure.";

AtomConservingQ::badarg = "Illegal argument of function AtomConservingQ.";

Options[AtomConservingQ]:={ExternalSpecies->{}};

SyntaxInformation[AtomConservingQ]={"ArgumentsPattern"->{{__},___}};

AtomConservingQ[{reactions__},externals:(_?VectorQ|OptionsPattern[])] :=
	Module[{ls, alpha, beta, atommatrix},
		{ls, alpha, beta} = Check[ReactionsData[{reactions},externals]["species","\[Alpha]","\[Beta]"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];
		atommatrix = First[ToAtomMatrix[ls]];
		(*Transpose[a].atommatrix == Transpose[b].atommatrix*)
		(#===ConstantArray[0,Length[First[atommatrix]]]) &/@ (Transpose[alpha-beta].atommatrix)
	];

AtomConservingQ[___] := (Message[AtomConservingQ::"badarg"]; $Failed)


AtomsQ::usage = "\!\(\*\nStyleBox[\"AtomsQ\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\)\
\!\(\*\nStyleBox[\"[\",\nFontFamily->\"Courier\"]\)\!\(\*\nStyleBox[\"list\",\nFontFamily->\"Courier\",\nFontSlant->\"Italic\",\nFontSize->10]\)\
\!\(\*\nStyleBox[\"]\",\nFontFamily->\"Courier\"]\) \
gives True if every elements of \!\(\*\nStyleBox[\"list\",\nFontFamily->\"Courier\",\nFontSlant->\"Italic\",\nFontSize->10]\) are the name of an atom, \
that is member of the list \!\(\*\nStyleBox[\"Atoms\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\) and gives False otherwise.";

AtomsQ::badarg = "Illegal argument of function AtomsQ.";

SyntaxInformation[AtomsQ] = {"ArgumentsPattern"->{_}};

AtomsQ[v_?VectorQ] := And@@(MemberQ[Atoms,#]&/@v);

AtomsQ[___] := (Message[AtomsQ::"badarg"]; $Failed)


FromAtomMatrix::usage = "FromAtomMatrix[atommatrix,atoms] constructs the molecular formula of the species from the atommatrix, \
if the members of the list atoms are members of the list Atoms. Otherwise...";

FromAtomMatrix::badarg = "Illegal argument of function FromAtomMatrix.";

SyntaxInformation[FromAtomMatrix]={"ArgumentsPattern"->{_,_}};

FromAtomMatrix[atommatrix_?MatrixQ,atoms_?AtomsQ] := 
	Inner[
		Switch[{#1,#2},{0,_},"",{_,"charge"},StringJoin[" ",If[Abs[#1]===1,"",Abs[#1]],If[#1<0,"-","+"]],{1,_},#2,_,StringJoin[#2,ToString[#1]]]&,
		atommatrix,
		Append[atoms,"charge"],
		StringJoin
	];

FromAtomMatrix[___] := (Message[FromAtomMatrix::"badarg"]; $Failed)


(* ::Subsubsection:: *)
(**DetailedBalanced*)


DetailedBalanced::usage = "DetailedBalanced[{reactions}, rates, opts] gives necessary and sufficient conditions \
to the given formal mechanism be detailed balanced. Options: ExternalSpecies->{}, GeneratedRateCoefficient->k, ...";


DetailedBalanced::badarg = "Illegal argument of function DetailedBalanced.";
DetailedBalanced::genparerr = "The option GeneratedRateCoefficient must be a symbol type variable.";
DetailedBalanced::notrev = "The given formal mechanism is not reversible, so we have no condition to be detailed balanced.";
DetailedBalanced::defeqz = "The given formal mechanism has deficiency zero, we are given only the circuit condition(s).";
DetailedBalanced::nocycle = "The FHJ-graph of the given formal mechanism has no cycle, so we are given only the spanning forest condition(s).";
DetailedBalanced::nocycdefz = "The given formal mechanism has deficiency zero as well as the FHJ-graph has no cycle, \
so the reaction is unconditionally detailed balanced.";
DetailedBalanced::rates = "The list of reactionsteps is not the same type with the given list of the reaction rates: `1`.";
DetailedBalanced::misscirc = "We may have missing circuit condition(s).";
DetailedBalanced::missspf = "We may have missing spanning forest condition(s).";

Options[DetailedBalanced] := { ExternalSpecies -> {}, GeneratedRateCoefficient -> Global`k };

SyntaxInformation[DetailedBalanced]={"ArgumentsPattern"->{{__},___}};

DetailedBalanced[{reactions__}, opts : OptionsPattern[]] :=
	Block[ { genpar, rr },

		genpar = GeneratedRateCoefficient /. Flatten[{opts, Options[DetailedBalanced]}];

		If[ Head[genpar] === Symbol,

			rr = Check[ReactionsData[{reactions}, ExternalSpecies /. Flatten[{opts, Options[DetailedBalanced]}]]["R"], Return[$Failed],
				 {ReactionsData::wrreac,ReactionsData::badarg}];

			DetailedBalanced[{reactions}, Array[Subscript[genpar,#]&,rr], opts],

			Message[DetailedBalanced::"genparerr"];
			Return[$Failed];

		]
	];

DetailedBalanced[{reactions__}, rates_?VectorQ, opts : OptionsPattern[]] := 
	Module[{ exs, m, r, nn, ll, fhj, gamma, deficiency, fhjun, rsteplistind, rrates,
			 fhjun0, mintree, mintreeun, c0, circs, spanreacs, mm, ccc, circconds, spanforconds
		   },

			If[Not[ReversibleQ[reactions]],
				Message[DetailedBalanced::"notrev"];
				Return[$Failed];
			];

			exs = ExternalSpecies /. Flatten[{opts, Options[DetailedBalanced]}];

			{m, r, nn, ll, fhj, gamma, deficiency} = 
					Check[ReactionsData[{reactions},exs]["M","R","n","l","fhjgraph","\[Gamma]","deficiency"], Return[$Failed], 
						 {ReactionsData::wrreac,ReactionsData::badarg}];

			If[Length[rates] =!= r,
				Message[DetailedBalanced::"rates",rates];
				Return[$Failed];
			];
			
			deficiency = ToExpression[StringTake[deficiency,-1]];

			fhjun = MakeUndirected[fhj];
			rsteplistind = ToOrderedPairs[fhj];
			rrates = Thread[rsteplistind -> rates];
			(**)
			fhjun0 = ToUnorderedPairs[fhjun];
			mintree = MinimumSpanningTree[fhjun];
			mintreeun = ToUnorderedPairs[mintree];
			(**)
			circconds := (
							If[Length[c0 = Complement[fhjun0, mintreeun]] =!= r/2-nn+ll,
								Message[DetailedBalanced::"misscirc"];
							];
							circs = Map[Partition[FindCycle[AddEdge[mintree,#]],2,1]&, c0];
							Map[Times @@ (#/.rrates) == Times @@ (Map[Reverse,#] /. rrates)&, circs]
						 );
			(**)
			spanforconds := (
							spanreacs = Flatten[Position[rsteplistind, #]& /@ mintreeun];
							If[Length[mm = NullSpace[Part[gamma, All, spanreacs]]] =!= deficiency,
								Message[DetailedBalanced::"missspf"];
							];
							Map[
								If[ And @@ (NumericQ /@ (ccc = {Times @@ ((mintreeun /. rrates)^#),Times @@ ((Map[Reverse,mintreeun] /. rrates)^#)})),
									ccc /. List -> Equal,
									Numerator[Together[Total[ccc]]] /. Plus -> Equal
								]&, mm
							]
							);

			If[ deficiency === 0,

				If[AcyclicQ[fhjun],

					Message[DetailedBalanced::"nocycdefz"];
					Return[True];,

					Message[DetailedBalanced::"defeqz"];
					circconds
				],

				If[AcyclicQ[fhjun],

					Message[DetailedBalanced::"nocycle"];
					spanforconds,

					Join[circconds,spanforconds]	
				]
			]
	];

DetailedBalanced[___] := (Message[DetailedBalanced::"badarg"]; Return[$Failed];)


(* ::Subsubsection::Closed:: *)
(**old stoichiometry (From David)*)


(*
ToCanonical[chains_List] :=
	DeleteDuplicates[
		Flatten[
		Flatten[
			chains //. (op1:(RightArrow | LeftArrow | LeftRightArrow | Rule | Equilibrium | ReverseEquilibrium))[x_, y_, z__] -> op1[op1[x, y],z]
			       //. (op1:(RightArrow | LeftArrow | LeftRightArrow | Rule | Equilibrium | ReverseEquilibrium))[
							(op2:(RightArrow | LeftArrow | LeftRightArrow | Rule | Equilibrium | ReverseEquilibrium))[x_, y_], z_
						] -> {op2[x, y], op1[y, z]}
		] /. Rule[x_, y_] -> RightArrow[x, y]
          /. Equilibrium[x_, y_] -> LeftRightArrow[x, y]
          /. ReverseEquilibrium[x_, y_] -> LeftRightArrow[x, y]
          /.	{
					LeftArrow[x_, y_]      -> RightArrow[y, x],
					LeftRightArrow[x_, y_] -> Sequence[RightArrow[x, y], RightArrow[y, x]]
				}
		]
	];

ToStoichiometry::usage = "ToStoichiometry[reactionlist] yields a list of the \
most important characteristics of the chemical mechanism given by the list of \
reaction list. ToStoichiometry[reactionlist, externalspecies] treats the species \
given as the function's second argument as external species having constant \
concentration.";

ToStoichiometry[reactionlist_List, externalspecies_List : {}] :=
	Module[{
			ReactionSteps,
			Species, Externals,
			Cmplxes, ComplexVectors,
			Reactants, ReactantVectors,
			Products, ProductVectors,
			alfa, beta,
			ToVariableName
		   },

		Externals     = Union[externalspecies, {"0"}];
		ReactionSteps = ToCanonical[reactionlist] //. A_ /; MemberQ[Externals, A] -> 0 /. 0 -> "0";
		Cmplxes       = Union[Cases[ReactionSteps, _, {2}]];
		Species       = Sort[Complement[
			Flatten[ReactionSteps //. {RightArrow[x_, y_] -> {x, y}, n_Integer x_ -> x, x_ + y_ -> {x, y}}],
			Externals
		]];
		Reactants = First /@ ReactionSteps;
		Products  = Last  /@ ReactionSteps;

		{ComplexVectors, ReactantVectors, ProductVectors} =
			(Outer[Coefficient[#1, #2] &, #1, Species, 1] &) /@ {Cmplxes, Reactants, Products};

		{alfa, beta} = Transpose /@ {ReactantVectors, ProductVectors};

		{
			Global`species -> Species,
			Global`complexes -> Cmplxes,
			Global`M -> Length[Species],
			Global`reactionSteps -> ReactionSteps,
			Global`R -> Length[ReactionSteps],
			Global`\[Alpha] -> alfa,
			Global`\[Beta] -> beta,
			Global`\[Gamma] -> beta-alfa,
			Global`variables -> (Subscript[Global`c, #] & /@ Species) 
		}
	];

FHJ[tst_List] :=
	FromOrderedPairs[Global`reactionSteps /. tst /. RightArrow[c1_, c2_] :>
		{Position[Global`complexes /. tst, c1][[1, 1]], Position[Global`complexes /. tst, c2][[1, 1]]}
	]

FHJComponents::usage = "FHJComponents[mech] gives the strongly connected components of the Feinberg-Horn-Jackson \
graph of the formal chemical mechanism mech.";

FHJComponents[tst_List] :=
	Module[{i},
		ConnectedComponents[MakeUndirected[FHJ[tst]]] /. i_Integer :> (Global`complexes /. tst)[[i]]
    ]

FHJStrongComponents::usage = "FHJStrongComponents[mech] gives the connected components of the (undirected) \
Feinberg-Horn-Jackson graph of the formal chemical mechanism mech.";

FHJStrongComponents[tst_List] :=
	Module[{i},
		StronglyConnectedComponents[FHJ[tst]] /. i_Integer :> (Global`complexes /. tst)[[i]]
    ]

RightHandSide::usage = "RightHandSide[mech, rrates] gives the right hand side of the \
kinetic differential equations assigned to the chemical mechanism mech. (The argument \
mech is a list yielded by the fuction ToStoichiometry.) The second argument contains \
the rate constants or rate functions of the reaction steps.";

RightHandSide[tst_List, reactionrates_?VectorQ] :=
	Module[{w, i, p, r, massactionrates = reactionrates Times @@ (Global`variables ^ Global`\[Alpha] /. tst)},
                w = Table[
	               If[FreeQ[reactionrates[[i]], Alternatives @@ (Prepend[Global`variables /. tst, t])],
	                  massactionrates[[i]],
		          reactionrates[[i]]
	               ],
                       {i, Global`R /. tst}
                ];
                Table[
		    Sum[(Global`\[Gamma] /. tst)[[m, r]] w[[r]], {r, 1, Global`R /. tst}],
		{m, 1, Global`M /. tst}]
	];

NonZeros::usage = "...";

Options[Kinetic] = {NonZeros -> {}};

Kinetic[tst_List, reactionrates_List, opts___?OptionQ] :=
    Module[{eqs, w, i, massactionrates = reactionrates Times @@ (Global`variables ^ Global`\[Alpha] /. tst)},
		w = Table[
				If[FreeQ[reactionrates[[i]], Alternatives @@ (Prepend[Global`variables /. tst, t])],
					massactionrates[[i]],
					reactionrates[[i]]
				],
				{i, Global`R /. tst}
			];

		eqs = Flatten[{	Thread[reactionrates != 0],
						Thread[(NonZeros /. Flatten[{opts, Options[Kinetic]}]) != 0],
						Thread[(Global`\[Gamma] /. tst).w == 0] }];
		Union[{ToRules[Reduce[eqs, Global`variables /. tst]]}]
	]

Concentrations[tst_List, reactionrates_?VectorQ, initialvalues_?VectorQ, opts___?OptionQ] :=
	Module[ {vars = Global`variables /. tst, attempt},

		attempt = DSolve[
				Join[
					Thread[ #'[Global`t]& /@ vars ==
							(RightHandSide[tst, reactionrates] /. (# -> #[Global`t] & /@ vars)) ], (* equations *)
					Thread[ (#[0]& /@ vars == initialvalues) ]                                      (* initial values *)
				],
				#[Global`t]& /@ vars,
				Global`t,
				FilterOptions[DSolve, opts]
			];
		If[Head[attempt]===List, First[attempt], attempt]
	]

Concentrations[tst_List, reactionrates_?VectorQ, initialvalues_?VectorQ, t0_?NumericQ, t1_?NumericQ, opts___?OptionQ] :=
	Module[ {vars = Global`variables /. tst, attempt},
		attempt = NDSolve[
				Join[
					Thread[ #'[Global`t]& /@ vars ==
							(RightHandSide[tst, reactionrates] /. (# -> #[Global`t] & /@ vars)) ], (* equations *)
					Thread[ (#[0]& /@ vars == initialvalues) ]                                     (* initial values *)
				],
				#[Global`t]& /@ vars,
				{Global`t, t0, t1},
				FilterOptions[NDSolve, opts]
			];
		If[Head[attempt]===List, First[attempt], attempt]
	];
*)


(* ::Subsection::Closed:: *)
(**deterministic approach*)


VariablesQ[x_] := And@@((Head[#]===Symbol)&/@x);

MassActionKineticFunction[alpha_,variables_,rates_] := rates*((Times @@ Power[variables,#]) &/@ Transpose[alpha]);

(*SpeciesQ := UpperCaseQ[StringTake[ToString[#],1]]&;*)


(* ::Subsubsection::Closed:: *)
(*RightHandSide, DeterministicModel*)


RightHandSide::usage = "RightHandSide[{reactions}, reactionrates, options] gives the right hand side of the \
kinetic differential equations assigned to the chemical mechanism reactions. The second argument contains \
the rate constants or rate functions of the reaction steps.";


RightHandSide::badarg = "Illegal argument of function RightHandSide.";

Options[RightHandSide] := {ExternalSpecies->{}};

RightHandSide[{reactions__}, rates_?VectorQ, OptionsPattern[]] := 
	RightHandSide[Sequence@@Check[ReactionsData[{reactions},OptionValue[ExternalSpecies]]["\[Alpha]","\[Beta]","variables"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}],rates];

RightHandSide[alpha_?MatrixQ,beta_?MatrixQ,variables_?VectorQ,rates_?VectorQ] := Dot[beta-alpha,MassActionKineticFunction[alpha,variables,rates]];

RightHandSide[___] := (Message[RightHandSide::"badarg"]; $Failed)


DeterministicModel::usage = "DeterministicModel[{reactions}, reactionrates, options] gives the symbolic form of the \
kinetic differential equations assigned to the chemical mechanism reactions... Several options...";


DeterministicModel::badarg = "Illegal argument of function DeterministicModel.";
DeterministicModel::rates = "The size of rate constants unequal to the number of given reactionsteps.";
DeterministicModel::vars = "The type of vector of variables must equal to the number of internal species.";
DeterministicModel::func = "The given rate function does not have an appropriate domain or image dimension.";

Options[DeterministicModel] := {ExternalSpecies->{}};

(*SyntaxInformation*)

DeterministicModel[{reactions__},s_Symbol:Global`t, opts : OptionsPattern[]] :=
	DeterministicModel[{reactions},Array[Subscript[Global`k,#]&,
		Check[ReactionsData[{reactions},OptionValue[ExternalSpecies]]["R"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}]],s, opts];

DeterministicModel[{reactions__},rates_?VectorQ,s_Symbol:Global`t,OptionsPattern[]] :=
	Module[{ alpha, beta, vars, reacs = Check[ReactionsData[{reactions},OptionValue[ExternalSpecies]]["R"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}]},
			If[reacs===Length[rates],
				{alpha, beta, vars} = ReactionsData[{reactions},OptionValue[ExternalSpecies]]["\[Alpha]","\[Beta]","variables"];
				vars = Through[vars[s]];
				{Thread[D[vars,s] == RightHandSide[alpha,beta,vars,rates]], vars},
				Message[DeterministicModel::"rates"];
				$Failed
			]
	];

(*?VariablesQ; Abort[]*)
DeterministicModel[{reactions__},rates_?VectorQ,vars_?VectorQ,s_Symbol:Global`t, opts : OptionsPattern[]] :=
	Module[{ eq, vnew, vold },
			{eq, vold} = 
				Check[DeterministicModel[{reactions},rates,s,opts],Abort[];,{DeterministicModel::"rates",ReactionsData::wrreac,ReactionsData::badarg}];
			If[Length[vars]===ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["M"],
				vnew = Through[vars[s]];
				{eq /. Thread[vold->vnew] /. Thread[D[vold,s]->D[vnew,s]], vnew},
				Message[DeterministicModel::"vars"];
				$Failed
			]
	];

(*vars van de rates nincs eset kezelese!*)

DeterministicModel[{reactions__},ratefunc_Function,s_Symbol:Global`t,OptionsPattern[]] := 
	Module[ {r, alpha, beta, vars, nonlin},
			{r, alpha, beta, vars} = Check[ReactionsData[{reactions},OptionValue[ExternalSpecies]]["R","\[Alpha]","\[Beta]","variables"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];
			vars = Through[vars[s]];
			nonlin = Check[ratefunc @@ vars, Message[DeterministicModel::"func"]; Return[$Failed];, {Function::"fpct",Function::"slotn",General::"stop"}];
			(*Check[ratefunc @@ Most[vars], Message[DeterministicModel::"func"]; Abort[], {Function::"fpct",Function::"slotn",General::"stop"}];*)
			If[r===Length[nonlin],
				{Thread[D[vars,s] ==  Dot[beta-alpha,nonlin]], vars},
				Message[DeterministicModel::"func"];
				$Failed
			]
	];

(*?VariablesQ*)
DeterministicModel[{reactions__},ratefunc_Function,vars_?VectorQ,s_Symbol:Global`t,OptionsPattern[]] := 
	Module[ {eq, vold, vnew},
			{eq, vold} = Check[DeterministicModel[{reactions},ratefunc,s,ExternalSpecies->OptionValue[ExternalSpecies]], Return[$Failed];,
							{ReactionsData::wrreac,ReactionsData::badarg,Function::"fpct",Function::"slotn",General::"stop",DeterministicModel::"func"}];
			If[Length[vars]===ReactionsData[{reactions},OptionValue[ExternalSpecies]]["M"],
				vnew = Through[vars[s]];
				{eq /. Thread[vold->vnew] /. Thread[D[vold,s]->D[vnew,s]], vnew},
				Message[DeterministicModel::"vars"];
				$Failed
			]
	];

DeterministicModel[___] := (Message[DeterministicModel::"badarg"]; $Failed)


(* ::Subsubsection::Closed:: *)
(*Concentrations*)


Concentrations::usage = "Concentrations[{reactions}, reactionrates, initialvalues, options] calculates the concentration-time \
functions of each species of the mechanism mech by solving the kinetic differential equaitons. The first \
argument is the list that describes the mechanism, the second is the list of reaction rates. (See also ReactionRatesNotebook.) \
The third argument sets the initial concentrations. Options... \
This way Concentrations uses DSolve, accepts symbolic parameters, and computes the symbolic solution.\n\
Concentrations[{reactions}, reactionrates, initialvalues, {t0, t1}, options] uses NDSolve to compute a numerical solution on the \
interval [t0, t1].";


Options[Concentrations] := {ExternalSpecies->{}};

Concentrations::badarg = "Illegal argument of function Concentrations.";
Concentrations::init = "The initialization vector does not have an appropiate size.";

(* Solving the ODE symbolically *)

Concentrations[{reactions__},rates_?VectorQ,init_?VectorQ,s_Symbol:Global`t,opts___?OptionQ] := 
	Concentrations[{reactions},rates,init,0,s,opts];

Concentrations[{reactions__},rates_?VectorQ,init_?VectorQ,vars_?VariablesQ,s_Symbol:Global`t,opts___?OptionQ] := 
	Concentrations[{reactions},rates,init,0,vars,s,opts];

Concentrations[{reactions__},rates_?VectorQ,init_?VectorQ,t0_?NumericQ(*Symbol*),s_Symbol:Global`t,opts___?OptionQ] :=
	Module[{kineq, vars, exs},
		exs = Sequence@@FilterRules[{opts},Options[Concentrations]];
		{kineq, vars} = Check[DeterministicModel[{reactions},rates,s,exs],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg,DeterministicModel::"rates"}];
		If[Length[init]===ReactionsData[{reactions},exs]["M"],
			{
				vars,
				First[DSolve[Join[kineq, Thread[(vars /. s -> t0) == init]], vars, s, FilterRules[{opts},Options[DSolve]]]]
			},
			Message[Concentrations::"init"];
			$Failed
		]
	];

(*?VariablesQ*)
Concentrations[{reactions__},rates_?VectorQ,init_?VectorQ,t0_?NumericQ(*Symbol*),vars_?VectorQ,s_Symbol:Global`t,opts___?OptionQ] :=
	Module[{kineq, newvars, exs},
		exs = Sequence@@FilterRules[{opts},Options[Concentrations]];
		{kineq, newvars} = Check[DeterministicModel[{reactions},rates,vars,s,exs],Return[$Failed];,
									{ReactionsData::wrreac,ReactionsData::badarg,DeterministicModel::"rates",DeterministicModel::"vars"}];
		If[Length[init]===ReactionsData[{reactions},exs]["M"],
			{
				newvars,
				First[DSolve[Join[kineq, Thread[(newvars /. s -> t0) == init]], newvars, s, FilterRules[{opts},Options[DSolve]]]]
			},
			Message[Concentrations::"init"];
			$Failed
		]
	];

Concentrations[{reactions__},ratefunc_Function,init_?VectorQ,t0_?NumericQ,s_Symbol:Global`t,opts___?OptionQ] :=
	Module[{kineq, vars, exs},
		exs = Sequence@@FilterRules[{opts},Options[Concentrations]];
		{kineq, vars} = Check[DeterministicModel[{reactions},ratefunc,s,exs],Return[$Failed];,
									{ReactionsData::wrreac,ReactionsData::badarg,Function::"fpct",Function::"slotn",General::"stop",DeterministicModel::"func"}];
		If[Length[init]===ReactionsData[{reactions},exs]["M"],
			{
				vars,
				First[DSolve[Join[kineq, Thread[(vars /. s -> t0) == init]], vars, s, FilterRules[{opts},Options[DSolve]]]]
			},
			Message[Concentrations::"init"];
			$Failed
		]
	];

(*?VariablesQ*)
Concentrations[{reactions__},ratefunc_Function,init_?VectorQ,t0_?NumericQ,vars_?VectorQ,s_Symbol:Global`t,opts___?OptionQ] :=
	Module[{kineq, newvars, exs},
		exs = Sequence@@FilterRules[{opts},Options[Concentrations]];
		{kineq, newvars} = Check[DeterministicModel[{reactions},ratefunc,vars,s,exs],Return[$Failed];,
									{ReactionsData::wrreac,ReactionsData::badarg,Function::"fpct",Function::"slotn",General::"stop",DeterministicModel::"func"}];
		If[Length[init]===ReactionsData[{reactions},exs]["M"],
			{
				newvars,
				First[DSolve[Join[kineq, Thread[(newvars /. s -> t0) == init]], newvars, s, FilterRules[{opts},Options[DSolve]]]]
			},
			Message[Concentrations::"init"];
			$Failed
		]
	];

(* Solving the ODE numerically *)

Concentrations[{reactions__},rates_?VectorQ,init_?VectorQ,{t0_?NumericQ,t1_?NumericQ},s_Symbol:Global`t,opts___?OptionQ] :=
	Module[{kineq, vars, exs},
			exs = Sequence@@FilterRules[{opts},Options[Concentrations]];
			{kineq, vars} = Check[DeterministicModel[{reactions},rates,s,exs],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg,DeterministicModel::"rates"}];
			If[Length[init]===ReactionsData[{reactions},exs]["M"],
				{ 
					vars,
					First[(*vars/.*)NDSolve[Join[kineq,Thread[(vars /. s -> t0) == init]], vars, {s, t0, t1}, FilterRules[{opts},Options[NDSolve]]]]
				},
				Message[Concentrations::"init"];
				$Failed
			]
	];

Concentrations[{reactions__},rates_?VectorQ,init_?VectorQ,{t0_?NumericQ,t1_?NumericQ},vars_?VariablesQ,s_Symbol:Global`t,opts___?OptionQ] :=
	Module[{kineq, newvars, exs},
			exs = Sequence@@FilterRules[{opts},Options[Concentrations]];
			{kineq, newvars} = Check[DeterministicModel[{reactions},rates,vars,s,exs],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg,DeterministicModel::"rates",DeterministicModel::"vars"}];
			If[Length[init]===ReactionsData[{reactions},exs]["M"],
				{ 
					newvars,
					First[NDSolve[Join[kineq,Thread[(newvars /. s -> t0) == init]], newvars, {s, t0, t1}, FilterRules[{opts},Options[NDSolve]]]]
				},
				Message[Concentrations::"init"];
				$Failed
			]
	];

(*vars-os valtozat, lent is, fent is*)

Concentrations[{reactions__},ratefunc_Function,init_?VectorQ,{t0_?NumericQ,t1_?NumericQ},s_Symbol:Global`t,opts___?OptionQ] :=
	Module[{kineq, vars, exs},
		exs = Sequence@@FilterRules[{opts},Options[Concentrations]];
		{kineq, vars} = Check[DeterministicModel[{reactions},ratefunc,s,exs],Return[$Failed];,
									{ReactionsData::wrreac,ReactionsData::badarg,Function::"fpct",Function::"slotn",General::"stop",DeterministicModel::"func"}];
		If[Length[init]===ReactionsData[{reactions},exs]["M"],
			{
				vars,
				First[(*vars/.*)NDSolve[Join[kineq,Thread[(vars /. s -> t0) == init]], vars, {s, t0, t1}, FilterRules[{opts},Options[NDSolve]]]]
			},
			Message[Concentrations::"init"];
			$Failed
		]
	];

Concentrations[___] := (Message[Concentrations::"badarg"]; $Failed)


(* ::Subsubsection::Closed:: *)
(*StationaryPoints, MassConservation, EigenvaluesJacobianMatrix*)


MyToRules[list_] := Flatten[If[Head[#]===Equal,ToRules[#],#]&/@list];


StationaryPoints::usage = "StationaryPoints[{reactions},rates,init]...";


Options[StationaryPoints] := {ExternalSpecies->{}, Positivity->False, Conditions->{}}; (*Variables->*)

StationaryPoints::badarg = "Illegal argument of function StationaryPoints.";
StationaryPoints::arglth = "One of the vector arguments `1`, `2` or `3` seem to have wrong sizes.";
StationaryPoints::noeql = "The current reaction simplex has no stationary points inside.";

StationaryPoints[{reactions__},opts___?OptionQ] := 
	Module[{r, vars},
		{r, vars} = Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["R","variables"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];
		StationaryPoints[{reactions},Array[Subscript[Global`k,#]&,r],Superscript[#,0]&/@vars,vars,opts]
	];

StationaryPoints[{reactions__},rates_?VectorQ,opts___?OptionQ] := 
	Module[{vars},
		vars = Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["variables"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];
		StationaryPoints[{reactions},rates,Superscript[#,0]&/@vars,vars,opts]
	];

StationaryPoints[{reactions__},rates_?VectorQ,init_?VectorQ,opts___?OptionQ] := 
	StationaryPoints[{reactions},rates,init,Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["variables"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}],opts];

StationaryPoints[{reactions__},rates_?VectorQ,init_?VectorQ,vars_?VectorQ,opts___?OptionQ]:=
	Module[{ equilibriums, ropts, alpha, beta, m, r, v, y, z, statp, balances, rateconst, inimass, nonneg, conds, sol, fless },
			ropts = Sequence@@FilterRules[{opts},Options[Reduce]];
			{alpha, beta, m, r} = 
					Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["\[Alpha]","\[Beta]","M","R"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];
			fless = If[(Positivity/.{opts})/.Positivity->False,Less,LessEqual];
			If[r===Length[rates]&&m===Length[init]&&m===Length[vars],
				v = Superscript[#,"*"]& /@ vars;
				z = Subscript[y,ToString[#]]&/@Range[r];
				{statp,balances,rateconst,inimass,nonneg} = 
					Thread/@{RightHandSide[alpha,beta,v,rates] == 0,(beta-alpha).z==v-init,rates > 0,Cases[init,Except[0]]>0,fless[0,v]};
				conds = Join[balances,rateconst,inimass,nonneg];
				If[(sol=Reduce[Join[statp,conds,(Conditions/.{opts})/.Conditions->{}],Join[v,z],Reals,Backsubstitution->True,ropts])===False,
					Message[StationaryPoints::noeql];
					$Failed,
					Join[{v},{Sort[Sort[MyToRules[Select[#/.And->List,Not[MemberQ[z,First[#]]]&]]]&/@Flatten[{sol/.Thread[rateconst->True]/.Thread[inimass->True]/.Or->List}]]}]
				],
				Message[StationaryPoints::arglth,rates,init,vars];
				$Failed
			]
	];

(*vars-os StatPoints...*)

StationaryPoints[___] := (Message[StationaryPoints::"badarg"]; $Failed)


MassConservation::usage = "...";


Options[MassConservation] := {ExternalSpecies->{}};

MassConservation::badarg = "Illegal argument of function MassConservation.";
MassConservation::arglth = "One of the vector arguments `1` or `2` seem to have wrong sizes. Both should have the size of internal species.";
MassConservation::nomass = "There is no mass to be conserved.";

MassConservation[{reactions__},init_?VectorQ,opts___?OptionQ] := 
	MassConservation[{reactions},init,Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["variables"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}],opts]

MassConservation[{reactions__},init_?VectorQ,vars_?VectorQ,opts___?OptionQ] :=
	Module[ {alpha, beta, m, ns, ropts, inimass, sol},
			{alpha, beta, m} = Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["\[Alpha]","\[Beta]","M"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];
			If[m===Length[init]&&m===Length[vars],
				ns = NullSpace[Transpose[beta-alpha]];
				ropts = Sequence@@FilterRules[{opts},Options[Reduce]];
				(*{mass, inimass, nonneg} = Thread /@ {ns.(vars - init) == 0, Cases[init,Except[0]]>0, vars >= 0 };*)
				inimass = Cases[init,Except[0]]>0;
				If[ns === {} || (sol=Reduce[Join@@(Thread /@ {ns.(vars - init) == 0, inimass, vars >= 0 }),vars,Reals,Backsubstitution->True,ropts])===False, 
					Message[MassConservation::nomass];
					$Failed,
					Join[{vars},{Sort[Sort[MyToRules[#/.And->List]]&/@Flatten[{sol/.Thread[inimass->True]/.Or->List}]]}]
				],
				Message[MassConservation::"arglth",init,vars];
				$Failed
			]
	];

MassConservation[___] := (Message[MassConservation::"badarg"]; $Failed)


EigenvaluesJacobianMatrix::usage = "...";


EigenvaluesJacobianMatrix::badarg = "Illegal argument of function EigenvaluesJacobianMatrix.";
EigenvaluesJacobianMatrix::arglth = "One of the arguments `1` or `2` seem to have wrong sizes. First one should have the size of reaction steps \
while the latter one of internal species.";

Options[EigenvaluesJacobianMatrix] := {ExternalSpecies->{}}; (*Variables->*)

EigenvaluesJacobianMatrix[{reactions__}, opts : OptionsPattern[]] := 
	Module[{r, vars},
		{r, vars} = Check[ReactionsData[{reactions},opts]["R","variables"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];
		EigenvaluesJacobianMatrix[{reactions},Array[Subscript[Global`k,#]&,r],vars,opts]
	];

EigenvaluesJacobianMatrix[{reactions__},rates_?VectorQ,opts : OptionsPattern[]] := 
	EigenvaluesJacobianMatrix[{reactions},rates,
			Check[ReactionsData[{reactions},opts]["variables"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}],opts];

EigenvaluesJacobianMatrix[{reactions__},rates_?VectorQ,vars_?VectorQ,OptionsPattern[]] := 
	Module[{alpha, beta, m, r, newvars},
		{alpha, beta, m, r} = Check[ReactionsData[{reactions},OptionValue[ExternalSpecies]]["\[Alpha]","\[Beta]","M","R"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];
		If[r===Length[rates]&&m===Length[vars],
			newvars = Superscript[#,"*"]& /@ vars;
			Eigenvalues[D[RightHandSide[alpha,beta,newvars,rates],{newvars}]],
			Message[EigenvaluesJacobianMatrix::arglth,rates,vars];
			$Failed
		]
	];

EigenvaluesJacobianMatrix[___] := (Message[EigenvaluesJacobianMatrix::"badarg"]; $Failed)


(* ::Subsubsection::Closed:: *)
(*AbsoluteRobustness*)


AbsoluteRobustness::usage = "complex...";


AbsoluteRobustness::badarg = "Illegal argument of function AbsoluteRobustness.";
AbsoluteRobustness::maybe = "`1`.";
AbsoluteRobustness::timeexc = "Time spent on calculations of positive stationary points exceeded `1` seconds, \
without being known whether it exists a positive stationary point or not in spite of the deficiency is one and \
there exist nonterminal complexes that differ only one species. Let the option TimeLimit be Infinity or try the function StationaryPoints.";

Options[AbsoluteRobustness] := { TimeLimit -> 60, ExternalSpecies -> {} };

AbsoluteRobustness[{reactions__},opts___?OptionQ] := 
	Module[{exs, time, delta, fhjstrong, fhjterminal, nonterms, acr},
		exs = Sequence @@ FilterRules[{opts},ExternalSpecies];
		time = (TimeLimit /. FilterRules[{opts},TimeLimit]) /. TimeLimit->60;
		{delta, fhjstrong, fhjterminal} = 
			Check[ReactionsData[{reactions},exs]["deficiency","fhjstronglyconnectedcomponents","fhjterminalstronglyconnectedcomponents"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];
		If[ Last[StringSplit[delta,"="]]=!="1",
			Message[AbsoluteRobustness::"maybe","The deficiency is different from 1"];
			Return[$Maybe],
			nonterms = Flatten[fhjstrong/.Thread[fhjterminal->Sequence[]]] /. "0"->0 /. Thread[exs->0];
			acr = Flatten[(nonterms-# &/@ nonterms) /. Plus[A_,B_]->0 /. 0->{}]; (*/. r_*A_:>A;*)
			If[ acr === {},
				Message[AbsoluteRobustness::"maybe","In spite of the deficiency is one there are no nonterminal complexes that differ only one species"];
				Return[$Maybe],
				Switch[TimeConstrained[Check[Length[StationaryPoints[{reactions},exs,Positivity->True,Sequence@@FilterRules[{opts},Conditions]]]>=2,False,{StationaryPoints::noeql}],time,"Abort"],
						"Abort",
							Message[AbsoluteRobustness::"timeexc",time];
							Return[$Maybe],
						False,
							Message[AbsoluteRobustness::"maybe","There are no positive stationary points"]; 
							Return[$Maybe],
						True,
							Union[acr /. r_*A_?StringQ :> A /. A_?StringQ*r_ :> A /. r_*A_?SpeciesQ :> A /. A_?SpeciesQ*r_ :> {A}](*!!!!*)
				]	
			]
		]
	];

AbsoluteRobustness[___][___] := (Message[AbsoluteRobustness::"badarg"]; $Failed)


(* ::Subsubsection::Closed:: *)
(**Volpert-indexing*)


VolpertIndexing::usage = "VolpertIndexing[{reactions}, initspecies, options] gives the Volpert indices \
of the species and reactions in the decomposition with respect to the initspecies given as \
the second argument.";


VolpertIndexing::badarg = "Illegal argument of function VolpertIndexing (see ?VolpertIndexing).";
VolpertIndexing::specerr = "None of the given species `1` is in the list of internal species of given mechanism.";

VolpertIndexing[alpha_?MatrixQ,beta_?MatrixQ, m0_?VectorQ] :=
	Module[
		{iterationlist = FixedPointList[
				{#[[1]](1 - Sign[beta.(1 - #[[2]])]),#[[2]]} & @ {#[[1]], #[[2]] Sign[#[[1]].alpha]} &,
						{
						 MapAt[# - 1 &, Array[1 &, Length[alpha]], Transpose[{m0}]],
						 Array[1 &, Length[First[alpha]]]
						}
						]
		},
		MapAt[# - 1 &, (Last[iterationlist] /. (1 -> Infinity)) + Total[iterationlist], 2]
	];

Options[VolpertIndexing] := {ExternalSpecies->{}};

VolpertIndexing[{reactions__}, initspecies_?VectorQ, OptionsPattern[]]:=
	Module[{alpha, beta, sp, w},

		{alpha, beta, sp} = Check[ReactionsData[{reactions},OptionValue[ExternalSpecies]]["\[Alpha]","\[Beta]","species"],Return[$Failed];,
							     {ReactionsData::wrreac,ReactionsData::badarg}];

		If[(w = Union[Flatten[Position[sp,#]&/@initspecies]]) =!= {},

			VolpertIndexing[alpha, beta, w],

			Message[VolpertIndexing::"specerr",initspecies];
			Return[$Failed]

		]

	];

VolpertIndexing[___][___] := (Message[VolpertIndexing::"badarg"]; $Failed)


(*
VolpertIndexingWithMatrices[A_?MatrixQ, B_?MatrixQ, M0_List] :=
	Module[
		{iterationlist = FixedPointList[{
						 #[[1]](1 - Sign[B.(1 - #[[2]])]),
						 #[[2]]
						} & @ {#[[1]], #[[2]] Sign[#[[1]].A]} &,

						{
						 MapAt[# - 1 &, Array[1 &, Length[A]], Transpose[{M0}]],
						 Array[1 &, Length[First[A]]]
						}
						]
		},
		MapAt[# - 1 &, (Last[iterationlist] /. (1 -> Infinity)) + Total[iterationlist], 2]
	]

Indexable::usage = "Indexable[decomp, atomm, init] yields True if and only if the decomposition with
respect to the given atommatrix is Volpert indexable, that is, all species have finite index.";
Indexable::badarg = "Illegal argument of function Indexable (see ?Indexable).";

Indexable[decomp_?DecompositionQ, atomm_?MatrixQ, initial_List] :=
 	Max[Last[VolpertIndexing[decomp, atomm, initial]]] < Infinity

Indexable[reactions : {__?ReactionQ}, initial_List] :=
	Max[Last[VolpertIndexing[reactions, initial]]] < Infinity

Indexable[___] := (Message[Indexable::"badarg"]; $Failed)

VolpertIndexing::usage = "VolpertIndexing[decomp, atomm, init] gives the Volpert indices \
of the species and reactions in the decomposition with respect to the atommatrix given as \
the second argument. The species given in the list init are the initial species. \
VolpertIndexing[decomp, atomm, init, tags] gives a pretty printed table of the results, \
with the reactions denoted by the strings given in the list tags.";

VolpertIndexing[decomp_?DecompositionQ, atomm_?MatrixQ, initial_List] :=
	VolpertIndexing[ atomm[[ Flatten[Position[decomp, _?Positive]] ]], initial ]

VolpertIndexing[reactions : {__?ReactionQ}, initial_List] :=
	VolpertIndexingWithMatrices[
		Transpose[UnitStep[-reactions - 1]],
		Transpose[UnitStep[ reactions - 1]],
		initial
	]

VolpertIndexing[decomp_?DecompositionQ, atomm_?MatrixQ, initial_List, species_List, reactionsigns_List] :=
	VolpertIndexing[ atomm[[ Flatten[Position[decomp, _?Positive]] ]], initial, species, reactionsigns]

VolpertIndexing[reactions : {__?ReactionQ}, initial_List, species_List, reactionsigns_List] :=
	Block[ {speciesindices, reactionindices},

		{speciesindices, reactionindices} = VolpertIndexingWithMatrices[
										Transpose[UnitStep[-reactions - 1]],
										Transpose[UnitStep[ reactions - 1]],
										initial
						    ];
		Print[
			TableForm[
				Transpose[{species, speciesindices}],
				TableHeadings -> {{}, {"Species", "Indices"}}
			]
		];
		Print[
			TableForm[
				Transpose[{reactionsigns, reactionindices}],
				TableHeadings -> {{}, {"Steps", "Indices"}}
			]
		];
	];
*)


(* ::Subsection:: *)
(**stochastic approach*)


(*MyDerivative[d_, f_, x_] := Derivative[Sequence @@ d][f][Sequence @@ x];*)
(*MassActionKineticFunction[alpha_,variables_,rates_] := rates*((Times @@ Power[variables,#]) &/@ Transpose[alpha]);*)


Intensity[alphatr_,X_,rates_] := (*Intensity[alphatr,X,rates] =*) rates*(Times@@FactorialPower[X,#] &/@ alphatr);

INTENSITY = Compile[{{rates,_Real,1},{x,_Integer,1},{alphatr,_Integer,2}},(rates*(Times@@FactorialPower[x,#] &/@ alphatr)),
	RuntimeAttributes->{Listable}, Parallelization->True, RuntimeOptions->{"CatchMachineUnderflow"->False,"CatchMachineIntegerOverflow"->False}]; (*Intensity*)
(*RuntimeOptions, "EvaluateSymbolically", RuntimeAttributes, COMPILE*)

SetAttributes[HornerForm,Listable];

IntensityDerivative[rates_,alphatr_,gamma_,vars_]:=
	HornerForm[D[rates*First[Outer[(FactorialPower[#1,#2]/.List->Times)&,{#},alphatr,1]]&@vars,{vars,1}].gamma]/.Thread[vars->#]&;

ChangeOfUnits[V_,rates_,alphatr_] := rates*((AvogadroConstant*Mole)*V)^(1 - Total /@ alphatr);


(*MyFilterRules[rule_,opts___] := (First[rule]/.FilterRules[{opts},First[rule]]) /. rule;*)

MyFilterOptions[f_, opt___]:= 
	If[$VersionNumber<6,
			Utilities`FilterOptions`FilterOptions[f, opt],
			Sequence@@FilterRules[Flatten[{opt}],Options[f]]
	];


(* ::Subsubsection::Closed:: *)
(*MassConservationQ*)


MassConservationQ::usage = "...";


Options[MassConservationQ] := { ExternalSpecies->{} };

MassConservationQ::badarg = "Illegal argument of function MassConservationQ.";
MassConservationQ::masscons = "The given formal reaction mechanism `1` does not conserve the mass.";

(*SyntaxInformation[MassConservationQ]={"ArgumentsPattern"->{{__},___}};*)

MassConservationQ[{reactions__},opts___?OptionQ] := 
	Module[{gamma, vars, ropts, sol},
			{gamma, vars} = Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["\[Gamma]","variables"],Return[$Failed];,
								 {ReactionsData::wrreac,ReactionsData::badarg}];
			ropts = MyFilterOptions[Reduce,opts];
			If[(sol=Reduce[Flatten[Thread/@{Transpose[gamma].vars==0,vars>0}],vars,Reals,Backsubstitution->True,ropts])=!=False,
				sol,
				Message[MassConservationQ::"masscons",{reactions}];
				Return[$Failed];
			]
	];

(* opts-ra pl. {"A"} es hibauzenet, Externalspecies *)

MassConservationQ[___] := (Message[MassConservationQ::"badarg"]; $Failed)


(* ::Subsubsection::Closed:: *)
(*ChemicalMasterEquation*)


ChemicalMasterEquation::usage = "...";


Options[ChemicalMasterEquation] := { ExternalSpecies->{} };

ChemicalMasterEquation::badarg = "Illegal argument of function ChemicalMasterEquation.";
ChemicalMasterEquation::rvarg = "Argument `1` or `2` have wrong sizes. \
First one should have the length of the reaction steps while the latter one of internal species including time parameter.";

ChemicalMasterEquation[{reactions__}, h_Symbol:Global`P, opts : OptionsPattern[]] := 
	ChemicalMasterEquation[{reactions},
		Array[Subscript[Global`k,#]&,Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["R"],Return[$Failed];,
										  {ReactionsData::wrreac,ReactionsData::badarg}]
			 ], h, opts];

ChemicalMasterEquation[{reactions__}, rates_?VectorQ, h_Symbol:Global`P, opts : OptionsPattern[]] := 
	ChemicalMasterEquation[{reactions},rates,
		Prepend[Array[Subscript[Global`x, #]&,Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["M"],Return[$Failed];,
												   {ReactionsData::wrreac,ReactionsData::badarg}]
					 ],Global`t
		], h, opts];

ChemicalMasterEquation[{reactions__}, rates_?VectorQ, vars_?VectorQ, h_Symbol:Global`P, OptionsPattern[]] := 
	Module[{alpha, gamma, m, r, s},
			{alpha, gamma, m, r} = Check[ReactionsData[{reactions},OptionValue[ExternalSpecies]]["\[Alpha]","\[Gamma]","M","R"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];
			{alpha, gamma} = Transpose /@ {alpha, gamma}; (*MyTranspose*)
			If[ Length[rates] === r && Length[vars] === m+1,
				D[h @@ vars,First[vars]] == 
					- (h @@ vars)*Total[FunctionExpand[Intensity[alpha, Rest[vars], rates]]] (*Expand!!!*)
					+ Total[Map[h@@(vars-Prepend[#,0])*Total[FunctionExpand[Intensity[alpha,Rest[vars]-#,rates]]]&,gamma]],
				Message[ChemicalMasterEquation::"rvarg", rates, vars];
				Return[$Failed];
			  ]
	];

(*rates vars csere*)

ChemicalMasterEquation[___] := (Message[ChemicalMasterEquation::"badarg"]; $Failed)


(* ::Subsubsection::Closed:: *)
(*StationaryProbabilityDensityEquation*)


StationaryProbabilityDensityEquation::usage = "...";


Options[StationaryProbabilityDensityEquation] := { ExternalSpecies->{} };

StationaryProbabilityDensityEquation::badarg = "Illegal argument of function ChemicalMasterEquation.";
StationaryProbabilityDensityEquation::rvarg = "Argument `1` or `2` have wrong sizes. \
First one should have the length of the reaction steps while the latter one of internal species.";

StationaryProbabilityDensityEquation[{reactions__}, h_Symbol:Global`\[CapitalPi], opts : OptionsPattern[]] := 
	StationaryProbabilityDensityEquation[{reactions},
		Array[Subscript[Global`k,#]&,Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["R"],Return[$Failed];,
										  {ReactionsData::wrreac,ReactionsData::badarg}]
			 ], h, opts];

StationaryProbabilityDensityEquation[{reactions__}, rates_?VectorQ, h_Symbol:Global`\[CapitalPi], opts : OptionsPattern[]] := 
	StationaryProbabilityDensityEquation[{reactions},rates,
		Array[Subscript[Global`x, #]&,Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["M"],Return[$Failed];,
										   {ReactionsData::wrreac,ReactionsData::badarg}]
			 ], h, opts];

StationaryProbabilityDensityEquation[{reactions__}, rates_?VectorQ, vars_?VectorQ, h_Symbol:Global`\[CapitalPi], OptionsPattern[]] := 
	Module[{alpha, gamma, m, r, s, exs = OptionValue[ExternalSpecies]},
			{m, r} = Check[ReactionsData[{reactions},exs]["M","R"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];
			If[ Length[rates] === r && Length[vars] === m,
				Last[ChemicalMasterEquation[{reactions},rates,Prepend[vars,Global`t],h,ExternalSpecies->exs]]==0/.(h[Global`t,q___]:>h[q]),
				Message[StationaryProbabilityDensityEquation::"rvarg",rates,vars];
				$Failed
			  ]
	];

StationaryProbabilityDensityEquation[___] := (Message[StationaryProbabilityDensityEquation::"badarg"]; $Failed)


(*
Solve Stationary Dist...
RSolve,
GeneratingFunction, etc.
*)


(* ::Subsubsection::Closed:: *)
(*ProbabilityGeneratingFunctionEquation, SolveProbabilityGeneratingFunctionEquation*)


ProbabilityGeneratingFunctionEquation::usage = "ProbabilityGeneratingFunctionEquation[{reactions},rates,vars,options] ...";


ProbabilityGeneratingFunctionEquation::badarg = "Illegal argument of function ProbabilityGeneratingFunctionEquation.";
ProbabilityGeneratingFunctionEquation::kvarg = "Argument `1` or `2` have wrong sizes. \
First one should have the length of the reaction steps while the latter one of internal species including time parameter.";

Options[ProbabilityGeneratingFunctionEquation] := { ExternalSpecies->{} };

(*SyntaxInformation*)

ProbabilityGeneratingFunctionEquation[{reactions__}, h_Symbol:Global`g, opts : OptionsPattern[]] := 
	ProbabilityGeneratingFunctionEquation[{reactions},
		Array[Subscript[Global`k,#]&,Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["R"],Return[$Failed];,
										  {ReactionsData::wrreac,ReactionsData::badarg}]
			 ], h, opts];

ProbabilityGeneratingFunctionEquation[{reactions__}, rates_?VectorQ, h_Symbol:Global`g, opts : OptionsPattern[]] := 
	ProbabilityGeneratingFunctionEquation[{reactions},rates,
		Prepend[Array[Subscript[Global`z, #]&,Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["M"],Return[$Failed];,
												   {ReactionsData::wrreac,ReactionsData::badarg}]
					 ], Global`t
			   ], h, opts];

ProbabilityGeneratingFunctionEquation[{reactions__}, rates_?VectorQ, vars_?VectorQ, h_Symbol:Global`g, OptionsPattern[]] := 
	Module[{alpha, beta, m, r, tv = First[vars], zv = Rest[vars]},
			{alpha, beta, m, r} = Check[ReactionsData[{reactions},OptionValue[ExternalSpecies]]["\[Alpha]","\[Beta]","M","R"],Return[$Failed];,
									   {ReactionsData::wrreac,ReactionsData::badarg}];
			If[ Length[rates] === r && Length[vars] === m+1,
				{alpha, beta} = Transpose /@ {alpha, beta};
				D[h @@ vars,tv]==Total[rates*MapThread[(Times@@(zv^#1)-Times@@(zv^#2))*Derivative[0,Sequence@@#2][h]@@vars&,{beta,alpha}]],
				Message[ProbabilityGeneratingFunctionEquation::"kvarg",rates,vars]; 
				Return[$Failed];
			  ]
	];

(*rates vars csere*)

ProbabilityGeneratingFunctionEquation[___] := (Message[ProbabilityGeneratingFunctionEquation::"badarg"]; $Failed)


SolveProbabilityGeneratingFunctionEquation::usage = "SolveProbabilityGeneratingFunctionEquation[{reactions},rates,init,vars,options] argument vars is optional.. \
\n Options: ExternalSpecies->{}, and Method->\"Built-in\" (default). Other methods: \"MatrixExponential\" and \"Characteristics\".";


SolveProbabilityGeneratingFunctionEquation::badarg = "Illegal argument of function ProbabilityGeneratingFunctionEquation.";
SolveProbabilityGeneratingFunctionEquation::initkvarg = "Argument `1`, `2` or `3` have wrong sizes. \
First one should have the length of the reaction steps, the second one of internal species including time parameter, \
while the latter one of internal species.";
SolveProbabilityGeneratingFunctionEquation::methoderr = "Argument `1` is not a valid optional method for this function. \
Try the methods \"Built-in\", \"Characteristics\" or \"MatrixExponential\".";
SolveProbabilityGeneratingFunctionEquation::orderlone = "The reactionorder of the given formal reaction mechanism is zero or 1. \
You can also try the method \"Characteristics\", if the following does not lead to appropiate result.";
SolveProbabilityGeneratingFunctionEquation::ordererr = "The method `1` cannot be applied to the given formal reaction mechanism, since its \
reactionorder is greater than 1. Try the method \"Built-in\".";
SolveProbabilityGeneratingFunctionEquation::betamaxlone = "The given formal reaction mechanism has no product species \
with stoichiometry coefficient greater than one. Thus, you can also apply the method \"MatrixExponential\", which can be more efficient.";
SolveProbabilityGeneratingFunctionEquation::betamaxerr = "The given formal reaction mechanism has at least one product species \
with stoichiometry coefficient greater than one. Thus, you cannot apply method `1`. Try one of the following methods: \"Characteristics\", \"Built-in\".";

Options[SolveProbabilityGeneratingFunctionEquation] := { ExternalSpecies->{}, Method->"Built-in" }; (*MatrixExponential, Characteristics*)

(* solve the pde symbolically *)

SolveProbabilityGeneratingFunctionEquation[{reactions__}, init_?VectorQ, h_Symbol:Global`g, opts___?OptionQ] := 
	SolveProbabilityGeneratingFunctionEquation[{reactions},
		Array[Subscript[Global`k,#]&,Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["R"],Return[$Failed];,
										  {ReactionsData::"wrreac",ReactionsData::"badarg"}]
			 ], init, h, opts];
(*
SolveProbabilityGeneratingFunctionEquation[{reactions__}, rates_?VectorQ, h_Symbol:Global`g, opts___?OptionQ] := 
	SolveProbabilityGeneratingFunctionEquation[{reactions},rates,
		Prepend[Array[Subscript[Global`z, #]&,Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["M"],Return[$Failed];,
												   {ReactionsData::wrreac,ReactionsData::badarg}]
					 ], Global`t
			   ], h, opts];
*)
SolveProbabilityGeneratingFunctionEquation[{reactions__}, rates_?VectorQ, init_?VectorQ, h_Symbol:Global`g, opts___?OptionQ]:=
	SolveProbabilityGeneratingFunctionEquation[{reactions}, rates, init, 
		Prepend[Array[Subscript[Global`z, #]&,Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["M"],Return[$Failed];,
												   {ReactionsData::wrreac,ReactionsData::badarg}]
					 ], Global`t
			   ], h, opts];
(*
SolveProbabilityGeneratingFunctionEquation[{reactions__}, rates_?VectorQ, vars_?VectorQ, h_Symbol:Global`g, opts___?OptionQ] := 
	Module[ {
			 pgfe = Check[ProbabilityGeneratingFunctionEquation[{reactions},rates,vars,h,FilterRules[{opts},ExternalSpecies]],Return[$Failed];,
						{ProbabilityGeneratingFunctionEquation::"kvarg",ProbabilityGeneratingFunctionEquation::"badarg"}]
			},
			{h @@ vars, Flatten[DSolve[pgfe,h @@ vars,vars,FilterRules[{opts},Options[DSolve]]]]}
	];
*)
SolveProbabilityGeneratingFunctionEquation[{reactions__}, rates_?VectorQ, init_?VectorQ, vars_?VectorQ, h_Symbol:Global`g, opts___?OptionQ] := 
	Module[{ 
			 exs = FilterRules[{opts},ExternalSpecies], tv = First[vars], zv = Rest[vars], rd, m, method, 
			 alphatr, betatr, gammatr, order, maxbeta, InhomogeneousMultip,
			 zp, matrixexp, u, s, zarg, z0, C, eqsol, pgfe
		   },
			rd = Check[ReactionsData[{reactions},exs],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];
			If[(m = rd["M"]) === Length[init] && Length[rates] === rd["R"] && Length[vars] === m+1,
				method = Method /. Flatten[{opts,Options[SolveProbabilityGeneratingFunctionEquation]}]; (*mindig az elso*)
				{alphatr, betatr, gammatr} = Transpose /@ rd["\[Alpha]","\[Beta]","\[Gamma]"];
				order = rd["reactionorder"];
				maxbeta = Max[betatr];
				InhomogeneousMultip[atr_,btr_,k_,x_,w_]:=
							Block[{
									zr = Flatten[Position[Total /@ Sign[atr], 0]] (*reactionstep order!*)
								  },
									Exp[Total[k[[zr]]*((Times@@(x^#)& /@ btr[[zr]])-1)]*w] (*inhomogeneous part*)
							];
				Switch[method,
					"Built-in",
						If[order <= 1,
							Message[SolveProbabilityGeneratingFunctionEquation::orderlone];
						];
						pgfe = (D[h @@ vars,tv]==Total[rates*MapThread[(Times@@(zv^#1)-Times@@(zv^#2))*Derivative[0,Sequence@@#2][h]@@vars&,{betatr,alphatr}]]);
						{h @@ vars, Flatten[DSolve[{pgfe,((h @@ vars)/.{tv->0})==Times@@(Rest[vars]^init)},h @@ vars,vars,FilterRules[{opts},Options[DSolve]]]]},
					"Characteristics",
						If[order > 1,
							Message[SolveProbabilityGeneratingFunctionEquation::ordererr,method];
							Return[$Failed];,
							If[maxbeta <= 1,
								Message[SolveProbabilityGeneratingFunctionEquation::betamaxlone];
							]
						];
						zarg = Through[zv[tv]];
						z0 = Array[C[#]&,m];
						eqsol = Flatten[
									DSolve[Join[Thread /@ {D[zarg,tv] == Total[-rates*MapThread[Times@@(zarg^#1)-Times@@(zarg^#2)&,{betatr,alphatr}]*alphatr],(zarg /. tv->0) == z0}], zarg, tv]] /. Rule->Equal;
						(*Quiet[multip*Times@@((z0/.Sort[Flatten[Solve[eqsol,z0]]/.Thread[zarg->zv]])^init),{Solve::"ifun"}]);*)
						InhomogeneousMultip[alphatr,betatr,rates,zv,tv]*Times@@((z0/.Sort[Flatten[Solve[eqsol,z0]]/.Thread[zarg->zv]])^init)
						(*{h @@ vars, h @@ vars -> ...}*)
						(*Print[eqsol]*),
					"MatrixExponential",
						If[order > 1,
							Message[SolveProbabilityGeneratingFunctionEquation::ordererr,method];
							Return[$Failed];,
							If[maxbeta>1,
								Message[SolveProbabilityGeneratingFunctionEquation::betamaxerr,method];
								Return[$Failed];,
								Message[SolveProbabilityGeneratingFunctionEquation::orderlone];
							]
						];
						zp = Flatten[Position[Total /@ Sign[betatr], 0]];
						matrixexp[u_] := matrixexp[u] = MatrixExp[Total[MapThread[KroneckerProduct,{alphatr, gammatr*rates}]]*u];
						InhomogeneousMultip[alphatr,betatr,rates,zv,tv]*Inner[Power,matrixexp[tv].zv+Total[rates[[zp]]*Integrate[Transpose[matrixexp[s]][[zp]],{s,0,tv}]],init,Times]
						(*{h @@ vars, h @@ vars -> ...}*)
						(*Exp[Total[rates[[zeros]]*(MapThread[Times,zv^betatr[[zeros]]]-1)]*tv]*),
					_,
						Message[SolveProbabilityGeneratingFunctionEquation::"methoderr",method];
						Return[$Failed];						
				],
				Message[SolveProbabilityGeneratingFunctionEquation::"initkvarg",rates,init,vars]; 
				Return[$Failed];
			]
	];

(*rates, init, vars ; init nelkul DSolve mukodhet*)

SolveProbabilityGeneratingFunctionEquation[___] := (Message[SolveProbabilityGeneratingFunctionEquation::"badarg"]; $Failed)


(* solve the pde numerically... *)


(* ::Subsubsection::Closed:: *)
(*ExpectedvalueProcess*)


NonNegativeIntegerNotZeroVectorQ := Max[#]>0 && MatchQ[#, {__Integer?NonNegative}]&;

MomentswithSpeciesQ := MatchQ[#,{_,_Integer?NonNegative}]&;

ExpectedSubscript[e_,newvars_,ylist_,vars_,var_] := 
	Block[{ coeffrules },
		coeffrules = CoefficientRules[Expand[FunctionExpand[Times@@FactorialPower[vars,ylist]]],vars];
		(coeffrules/.Rule[a_,b_]:>b).Map[Subscript[e,ToString[#,StandardForm]][var]&,coeffrules/.Rule[a_,b_]:>Inner[Power,newvars,a,Times]]
	];

(*Subscript[e,ToString[FunctionExpand[Times@@FactorialPower[newvars,{y}]]]][var]*)


ExpectedvalueProcess::usage = "ExpectedvalueProcess[{reactions},a,rates,...] ... Options: CombinatorialMoments, ExternalSpecies";


ExpectedvalueProcess::badarg = "Illegal argument of function ExpectedvalueProcess.";
ExpectedvalueProcess::aarg = "Argument `1` have wrong size. It should have the size of internal species.";
ExpectedvalueProcess::nonnegarg = "Argument `1` must be a not identically zero vector having nonnegative integer entries.";
ExpectedvalueProcess::speciesarg = "Argument `1` contains only invalid internal species.";

Options[ExpectedvalueProcess] := { CombinatorialMoments->False, ExternalSpecies->{} };

(*SyntaxInformation ...*)

ExpectedvalueProcess[{reactions__}, a_, opts : OptionsPattern[]] :=
	ExpectedvalueProcess[{reactions},a,
		Array[Subscript[Global`k,#]&,Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["R"],Return[$Failed];,
										  {ReactionsData::"wrreac",ReactionsData::"badarg"}]
			 ], Global`\[DoubleStruckCapitalE], Global`t, opts];

ExpectedvalueProcess[{reactions__}, a_, rates_?VectorQ, opts : OptionsPattern[]] :=
	ExpectedvalueProcess[{reactions}, a, rates, Global`\[DoubleStruckCapitalE], Global`t, opts];

(*{e_Symbol, var_Symbol} or e_Symbol:Global`E..., optional*)
ExpectedvalueProcess[{reactions__}, a_, e_Symbol, var_Symbol, opts : OptionsPattern[]] :=
	ExpectedvalueProcess[{reactions},a,
		Array[Subscript[Global`k,#]&,Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["R"],Return[$Failed];,
										  {ReactionsData::"wrreac",ReactionsData::"badarg"}]
			 ], e, var, opts];

(* arg. csere*)

ExpectedvalueProcess[{reactions__}, a_, rates_?VectorQ, e_Symbol, var_Symbol, OptionsPattern[]] := 
	Module[{ 
			species, m, l, sppos, a1, pgfe, vars, newvars, exp, es,
			exs = OptionValue[ExternalSpecies], combmom = OptionValue[CombinatorialMoments] 
		   },
			{species, m} = Check[ReactionsData[{reactions},exs]["species","M"],Return[$Failed];,{ReactionsData::"wrreac",ReactionsData::"badarg"}];
			Switch[a, 
				{__?MomentswithSpeciesQ},
					If[(sppos = Flatten[Position[species,#]& /@ First[a1 = Transpose[a]]]) =!= {},
						ExpectedvalueProcess[{reactions}, #, rates, e, var, ExternalSpecies->exs, CombinatorialMoments->combmom]& /@ ((UnitVector[m,#]& /@ sppos)*Last[a1]),
						Message[ExpectedvalueProcess::"speciesarg",a];
						Return[$Failed];
					],
				_?NonNegativeIntegerNotZeroVectorQ,
					l = Length[a];
					If[l === m,
						If[combmom,
							pgfe = Check[ProbabilityGeneratingFunctionEquation[{reactions},rates,ExternalSpecies->exs],Return[$Failed];,
										{ProbabilityGeneratingFunctionEquation::kvarg,ProbabilityGeneratingFunctionEquation::badarg}];
							vars = First[pgfe] /. Derivative[x__][y_][z__]:>{z};
							newvars = ToString /@ species;
							D[pgfe,Sequence@@Thread[List[Rest[vars],a]]]/.Thread[Rest[vars]->1]
												/. Derivative[x_,y__][f_][z__]:>D[ExpectedSubscript[e,newvars,{y},Rest[vars],var],{var,x}]
												/. Global`g[Global`t,Sequence@@Array[1&,l]]->1 /. func_[Global`t]:>func[var],
							D[Subscript[e,ToString[Inner[Power,species,a,Times],StandardForm]][var],var] == 
								Total[(Times@@StirlingS2[a,#])*Last[ExpectedvalueProcess[{reactions},#,rates,e,var,ExternalSpecies->exs,CombinatorialMoments->True]]& /@ Rest[Tuples[Prepend[Range[Max[a]],0],l]]]
						],
						Message[ExpectedvalueProcess::"aarg",a];
						Return[$Failed];
					],
					_,
						Message[ExpectedvalueProcess::"nonnegarg",a];
						Return[$Failed];
			]
	];

(*  
es = Union[Cases[exp,Subscript[h,_][var],3],{Integrate[First[exp],var]}];
{Collect[exp,es],es},
Expand...
*)

ExpectedvalueProcess[___] := (Message[ExpectedvalueProcess::"badarg"]; $Failed)


(* SolveExpectedvalueProcess 

*)


(* ::Subsection::Closed:: *)
(*stochastic simulation algos*)


MonitorMSG[prop_,pill_,species_] := 
		Column[{
				Row[{ProgressIndicator[prop]," " <> ToString[Ceiling[100*prop]]<>"%"}],
				"Current time: " <> ToString[First[pill],TraditionalForm] <> " units.",
				Grid[Prepend[Thread[List[ToString /@ species, Last[pill]]], {"Species","Current state"}],
					 Dividers->{{False,True}, {False,True}}]
			   }
		];


(* ::Subsubsection::Closed:: *)
(*DependencyGraph*)


(* Memo *)
DependencyGraphPrivate[alphatrules_, gammatrules_, r_] := DependencyGraphPrivate[alphatrules, gammatrules, r] = 
	Map[(Min[(First[#]/.gammatrules)+(Last[#]/.alphatrules)]<Infinity)/.{True->Rule[#,1],False->Sequence[]}&,Tuples[Range[r],2]];
(* SparseArray, Normal! AdjacencyMatrix *)

(* 
Module[{res = {}, a = alphatrindeces/.(0->Infinity), g = gammatrindeces/.(0->Infinity), i, j},
			Do[ If[Min[(i/.g)+(j/.a)]<Infinity,
					AppendTo[res,{i,j}->1]
				],{i, r}, {j, r}
			];
			res
	];
*)
(* 
DG = Compile[{{alpha, _Integer, 2},{gamma, _Integer, 2},{r, _Integer}},
		Module[{res = {}, ran = r, a = alpha, g = gamma, i, j},...
*)
(*
DependencyGraph[alphatr_,gammatr_,m_,r_] := 
	Map[#->Boole[Length[Union[Flatten[(Position[#,0]&)/@{Part[gammatr,First[#]],Part[alphatr,Last[#]]}]]]<m]&,Tuples[Range[r],2]];
*)


DependencyGraph::usage = "DependencyGraph[{reactions},externals] gives the dependency graph of the given formal reaction mechanism with \
(optional) externalspecies: externals. More help...";


DependencyGraph::badarg = "Illegal argument of function DependencyGraph.";

Options[DependencyGraph] := { ExternalSpecies->{} };

(*SyntaxInformation*)

(*errors?*)
DependencyGraph[{reactions__}, externals : (_?VectorQ | OptionsPattern[])] :=
	Module[{m, r, rsteps, alpha, gamma },
			{m, r, rsteps, alpha, gamma} = Check[ReactionsData[{reactions},externals]["M","R","reactionsteps","\[Alpha]","\[Gamma]"],Return[$Failed];,
												{ReactionsData::wrreac,ReactionsData::badarg}];
			{alpha, gamma} = Transpose /@ {alpha,gamma};
			((Length[Union[Flatten[(Position[#,0]&)/@{Part[gamma,First[#]],Part[alpha,Last[#]]}]]] < m) /. {True->Rule @@ Part[rsteps,#],False->Sequence[]})& /@ Tuples[Range[r],2]
			(*Flatten[Array[{#1,#2}&,{R,R}],1]*)
	];

DependencyGraph[__] := (Message[DependencyGraph::"badarg"]; $Failed);


ShowDependencyGraph::usage = "ShowDependencyGraph[{reactions},options] visualizes the dependency graph of the given formal reaction mechanism with \
(optional) externalspecies: externals. Options: PlotFunction, ExternalSpecies.. More help...";


ShowDependencyGraph::badarg = "Illegal argument of function ShowDependencyGraph.";

Options[ShowDependencyGraph] := { ExternalSpecies->{}, PlotFunction->"GraphPlot" };(*plotfunction!*)

(*SyntaxInformation*)

ShowDependencyGraph[{reactions__},opts___?OptionQ] :=
	GraphPlot[DependencyGraph[{reactions},FilterRules[{opts},ExternalSpecies]],FilterRules[{opts},Options[GraphPlot]]];

ShowDependencyGraph[__] := (Message[ShowDependencyGraph::"badarg"]; $Failed);


(* ::Subsubsection::Closed:: *)
(*Direct- and First Reaction Methods*)


(*$IterationLimit, $RecursionLimit*)
(* t0 *)
DirectMethod[species_,alphatr_,gammatr_,rates_,X0_,maxtime_,maxiteration_,verbose_] :=
	Module[{it, itotal, nwl, pill, cfint},
		nwl := NestWhileList[
				(pill:=#;(*sorrend, lehet a vegere*)
				# + {RandomReal[ExponentialDistribution[itotal]],First[RandomChoice[it->gammatr,1]]})&,
				{0,X0}, (First[#]<=maxtime && (itotal = Total[it = INTENSITY[rates,Last[#],alphatr]])>0)&, 1, maxiteration
			   ];
		If[verbose,
			Monitor[nwl,
				MonitorMSG[First[pill]/maxtime,pill,species]
			],
			nwl
		]
	];

MyExponentialDistribution[0.`] := Infinity;
MyExponentialDistribution[0] := Infinity;
(*if Intensity is used, then MyExponentialDistribution[0] := Infinity; machine prec.!*)
MyExponentialDistribution[x_] := RandomReal[ExponentialDistribution[x]];

SetAttributes[MyExponentialDistribution,Listable];

FirstReactionStepMethod[species_,alphatr_,gammatr_,rates_,X0_,maxtime_,maxiteration_,verbose_] := 
	Module[{itd, it, min, pill, nwl},
		nwl := NestWhileList[
				(pill := #;(*sorrend, lehet a vegere*)
				itd = MyExponentialDistribution[it]; (*Parallelize*)
				min = Min[itd];
				# + {min, min/.Thread[itd->gammatr]})&, (*Through[Join[List,Part[gammatr,First[Position[itd,#]]]&][Min[itd]]]*)
				{0,X0}, (First[#]<=maxtime && Total[it = INTENSITY[rates,Last[#],alphatr]]>0)&, 1, maxiteration
			   ];
		If[verbose,
			Monitor[nwl,
				MonitorMSG[First[pill]/maxtime,pill,species]
			],
			nwl
		]
	];


(* ::Subsubsection::Closed:: *)
(*NextReactionStep Method (!)*)


Indexing[m_] := MapIndexed[First[#2]->#1 &, m];

Sorting[v_] := Sort[v, #1[[2]] < #2[[2]] &];


(*a_new ?? in point c) and t<-tau????, dg sima graph, nem kell sparsearray, it -indexed?, t_ 1, t_ 2, milyen idok?*)

NextReactionStepMethod[species_,m_,r_,alphatr_,gammatr_,rates_,X0_,maxtime_,maxiteration_,verbose_] :=
	Module[{
			alphatrules, gammatrules, dg, it, tau, nwl, pill, next, 
			tmu1, tmu2, newrs, newits, store, oldprop, newprop, up
		   },
		{alphatrules, gammatrules} = Dispatch[Indexing[#]]& /@ {alphatr, gammatr};
		dg = DependencyGraphPrivate[alphatrules/.(0->Infinity),gammatrules/.(0->Infinity),r]; (*Inflow! if[{},...]*)
		it = Intensity[alphatr,X0,rates]; (*INTENSITY[rates,X0,alphatr], {} problem*)
		(*Print[{it,dg}];*)
		{tau, it} = {Sorting[Indexing[MyExponentialDistribution[it]]], Indexing[it]};
		store = Cases[it, Rule[x_,0] :> Rule[x,{0,0}]];
		nwl := NestWhileList[
				(pill := #;(*vegere*)
				{tmu1, tmu2} = Reverse[First[tau]/.Rule->List];
				next = # + {tmu1, tmu2/.gammatrules};
				newrs = Cases[dg,Rule[{tmu2,x_},1]:>x]; (*lehet olyan, hogy ez ures? igen: Inflow: {}*)
				(*Print[newrs];*)
				newits = Thread[Rule[newrs,Intensity[newrs/.alphatrules,Last[next],Part[rates,newrs]]]];
				(*INTENSITY[Part[rates,newrs],Last[next],newrs/.alphatrules]*)
				tau = Sorting[
						tau /. (Rule[#,x_] -> 
									Switch[ {newprop = #/.newits, oldprop = #/.it},
										{0, 0}, Rule[#,Infinity],
										{0, _}, AppendTo[store, #->{oldprop, #/.tau}];
												Rule[#,Infinity],
										{_, 0}, up = # /. store;
												store /. (Rule[#,x_]->Sequence[]);
(*Michaelis-Menten, kell meg, hogy First[up] ===0 !!!!!*)If[(#/.tau) === Infinity, Rule[#,tmu1], Rule[#, First[up]*((#/.tau) - Last[up])/newprop + tmu1]],
										{_,_},  Rule[#, oldprop*((#/.tau) - tmu1)/newprop + tmu1]
									]& /@ (newrs /. tmu2 -> Sequence[])) 
							/. Rule[tmu2, x_] -> Rule[tmu2, MyExponentialDistribution[tmu2/.newits] + tmu1]
					  ];
				it = it /. Thread[Thread[Rule[newrs,x_]]->newits];
				next )&,
				{0,X0}, (First[#]<=maxtime && Total[Last /@ it]>0)&, 1, maxiteration
			   ];
		If[verbose,
			Monitor[nwl,
				MonitorMSG[First[pill]/maxtime,pill,species]
			],
			nwl
		]
	];


(* ::Subsubsection::Closed:: *)
(*Tau - Leaping Methods*)


MyPoissonDistribution[0.`] := 0;
MyPoissonDistribution[0] := 0;
MyPoissonDistribution[x_] := RandomInteger[PoissonDistribution[x]];

SetAttributes[MyPoissonDistribution,Listable];


ExplicitTauLeapingMethod[species_,alphatr_,gammatr_,rates_,X0_,maxtime_,maxiteration_,stepsize_,verbose_] :=
	Module[{ nwl, pill, it, tau = stepsize },
		nwl := NestWhileList[
				(pill := #;(*sorrend, lehet a vegere*)
				# + {tau, Total[MyPoissonDistribution[it*tau]*gammatr]})&, (*Parallelize, ido most nem is lenyeges*)
				{0,X0}, (First[#]<=maxtime && Min[Last[#]]>=0 && Total[it=INTENSITY[rates,Last[#],alphatr]]>0)&, 1, maxiteration
			   ];
		If[verbose,
			Monitor[nwl,
				MonitorMSG[First[pill]/maxtime,pill,species]
			],
			nwl
		]
	];


(*condition positivity*)


ImplicitTauLeapingMethod[species_,alphatr_,gammatr_,rates_,vars_,X0_,maxtime_,maxiteration_,stepsize_,verbose_] :=
	Module[{nwl, pill, it, tau = stepsize, poi, rhs, sol, nonneg, intvars},
		nonneg = Thread[vars>=0];
		intvars = Intensity[alphatr,vars,rates];
		nwl := NestWhileList[
				(pill := #;(*sorrend, lehet a vegere*)
				poi = MyPoissonDistribution[it*tau];
				rhs = poi-it*tau+intvars*tau;(*tau*???*)
				sol = Flatten[FindInstance[Join[nonneg,Thread[vars-Last[#]==Total[rhs*gammatr]]],vars,Reals]];(*ha nincs megoldas, osszefuzes mas modszerekkel*)
				# + {tau, Total[Round[rhs/.sol]*gammatr]})&,
				{0,X0}, (First[#]<=maxtime && Min[Last[#]]>=0 && Total[it=INTENSITY[rates,Last[#],alphatr]]>0)&, 1, maxiteration
			   ];
		If[verbose,
			Monitor[nwl,
				MonitorMSG[First[pill]/maxtime,pill,species]
			],
			nwl
		]
	];


TrapezoidalTauLeapingMethod[species_,alphatr_,gammatr_,rates_,vars_,X0_,maxtime_,maxiteration_,stepsize_,verbose_] :=
	Module[{nwl, pill, it, tau = stepsize, poi, rhs, sol, nonneg, intvars},
		nonneg = Thread[vars>=0];
		intvars = Intensity[alphatr,vars,rates];
		nwl := NestWhileList[
				(pill := #;(*sorrend, lehet a vegere*)
				poi = MyPoissonDistribution[it*tau];
				rhs = poi-it*tau/2+intvars*tau/2;(*alt. int kozelites*)
				sol = Flatten[FindInstance[Join[nonneg,Thread[vars-Last[#]==Total[rhs*gammatr]]],vars,Reals]];(*gamma.rhs, mi van ha nincs megoldas, osszefuzes mas modszerekkel*)
				# + {tau, Total[Round[rhs/.sol]*gammatr]})&,
				{0,X0}, (First[#]<=maxtime && Min[Last[#]]>=0 && Total[it=INTENSITY[rates,Last[#],alphatr]]>0)&, 1, maxiteration
			   ];
		If[verbose,
			Monitor[nwl,
				MonitorMSG[First[pill]/maxtime,pill,species]
			],
			nwl
		]
	];


(* ::Subsubsection::Closed:: *)
(*K - Leaping Method*)


(*condition positivity*)


KLeapingMethod[species_,alphatr_,gammatr_,rates_,X0_,maxtime_,maxiteration_,error_,verbose_] :=
	Module[{gamma, nwl, pill, it, itotal, old, K, monitormsg, sz},
		gamma = Transpose[gammatr]; (**)
		nwl := NestWhileList[
				(pill := #;
				sz := First[pill]/maxtime;
				monitormsg := MonitorMSG[sz,pill,species];
							(* Column[{Row[{ProgressIndicator[sz]," " <> ToString[Ceiling[100*sz]]<>"%"}],
									"Current time: " <> ToString[First[pill],TraditionalForm] <> " units.\nCurrent state: " <> ToString[Last[pill]] <> "."}];*)
				old = MapIndexed[{First[#2],#1}&, Last[#]];
				If[(K = Round[Max[Min[(Max[error*Last[#],1] / Max[Abs[Part[gamma,First[#]]]])& /@ old],1]]) === 1,
					# + {RandomReal[ExponentialDistribution[itotal]], First[RandomChoice[it->gammatr,1]]},
					# + {RandomReal[GammaDistribution[K,1/itotal]], gamma.RandomInteger[MultinomialDistribution[K,it/itotal]]}
				])&,
				{0,X0}, (First[#]<=maxtime && (itotal=Total[it=INTENSITY[rates,Last[#],alphatr]])>0)&, 1, maxiteration
			   ];
		If[verbose,
			Monitor[nwl,monitormsg], (**)
			nwl
		]
	];


(* ::Subsubsection::Closed:: *)
(*Simulation*)


Simulation::usage = "Simulation[{reactions}, rates, init, maxtime, options] simulates the \
mechanism reactions with reaction rates given in the list rates. The initial state is determined \
by init, the stopping time is specified by maxtime. The result is given in the format {time, state}, \
where time is a positive real, when a reaction has occured, and state is the list of the integers, the corresponding state.\n\
Option name\t Default value\n ExternalSpecies\t {},\n MaxIteration\t Infinity,\n Method\t \"Direct\",\n \
Verbose\t True,\n Volume\t False.\n\
other options... more help needed. Methods: \"Direct\", \"FirstReactionStep\", \"ModifiedFirstReactionStep\", \
\"NextReactionStep\", \"ExplicitTau-Leaping\", \"ImplicitTau-Leaping\", \"TrapezoidalTau-Leaping\", \"K-Leaping\".";

(*Options*)
FixedStepSize::usage = "FixedStepSize ...";
MaxIteration::usage = "MaxIteration ...";
(*
Method::usage = "Methods: \"Direct\", \"FirstReactionStep\", \"ModifiedFirstReactionStep\", \
\"NextReactionStep\", \"ExplicitTau-Leaping\", \"ImplicitTau-Leaping\", \"TrapezoidalTau-Leaping\", \"K-Leaping\"."
*)
Verbose::usage = "Verbose is an option of multiple functions such as Concentrations, \
CoveringDecompositionSet etc. that makes the functions give more detailed information \
during the calculations. (The results are unaffected.) Default value may vary. Enable \
it using Verbose -> True (default).";
Volume::usage = "Volume ...";


Options[Simulation] = { ExternalSpecies->{}, FixedStepSize -> 0.012, MaxIteration -> Infinity, Method -> "Direct", Verbose -> True, Volume -> Infinity }; 
(*Error, StartingStepSize, Methods*)

Simulation::badarg = "Illegal argument of function Simulation.";
Simulation::method = "The given method '`1`' have not recognized yet. Try one of the following: \"Direct\", \"FirstReactionStep\", ...";
Simulation::time = "Argument `1` is not positive. The maximum time parameter is expected to be a positive real number.";
Simulation::rsarg = "Argument '`1`' or '`2`' have wrong sizes or bad entries. First one is expected to be the size of reaction steps \
while the latter one of internal species. \
Each of the vectors must have nonnegative real and integer entries, respectively.";
Simulation::volarg = "Argument `1` is not positive. The volume is expected to be a positive real number.";

SyntaxInformation[Simulation]={"ArgumentsPattern"->{{__},{__},{__},_,OptionsPattern[]}};

Simulation[{reactions__}, rates_?VectorQ, X0_?VectorQ, maxtime_?NumericQ, opts : OptionsPattern[]] := 
	Module[{
			exs, fss, maxit, method, verbose, v, cou, rd, res, alpha, beta, gamma, 
			m, r, sp, vars, alphatr, betatr, gammatr, x, error
			},
		{exs, fss, maxit, method, verbose, v} = {ExternalSpecies, FixedStepSize, MaxIteration, Method, Verbose, Volume} /. Flatten[{opts, Options[Simulation]}];
		If[Positive[v],(*NumberQ[v]*)
			cou[vol_,atr_] := If[vol===Infinity, rates, ChangeOfUnits[vol,rates,atr]];,
			Message[Simulation::"volarg",v];
			Return[$Failed];
		];
		rd := Check[ReactionsData[{reactions},Flatten[{exs}]]["\[Alpha]","\[Beta]","\[Gamma]","M","R","species","variables"], Return[$Failed];,
				   {ReactionsData::wrreac,ReactionsData::badarg}]; 
		If[verbose,
			x = Dynamic[Refresh[Round[Clock[Infinity]],UpdateInterval->1]];
			Monitor[res = rd,
						Column[{Row[{"ReactionsData is now calculating ",Dynamic[Mod[First[x],5]/.{0->"",1->".",2->"..",3->"...",4->"...."}]}],
								ProgressIndicator[Dynamic[Clock[Infinity]],Indeterminate,ImageMargins->1,BaselinePosition->Center],
								Row[{x,". ","seconds passed"}],
								Row[{"Simulation process will be started soon."}]
						}]	
			];
			FinishDynamic[];,
			res = rd;
		];
		{alpha, beta, gamma, m, r, sp, vars} = res;
		{alphatr, betatr, gammatr} = Transpose /@ {alpha, beta, gamma};
		(*calc := Check[ReactionsData[Equilibrium@@#&/@Partition[Array[Subscript[Global`A,#]&,1000],3,2]][],Return[$Failed];,
					   {ReactionsData::wrreac,ReactionsData::badarg}];*)
		If[r===Length[rates] && m===Length[X0] && Min[X0]>=0 && Min[rates]>=0,
			If[maxtime>0,
				Switch[method,
					"Direct",
						DirectMethod[sp,alphatr,gammatr,cou[v,alphatr],X0,maxtime,maxit,verbose],
					"FirstReactionStep",
						FirstReactionStepMethod[sp,alphatr,gammatr,cou[v,alphatr],X0,maxtime,maxit,verbose],
					"ModifiedFirstReactionStep", 
						Return[$Failed];,
					"NextReactionStep", 
						NextReactionStepMethod[sp,m,r,alphatr,gammatr,cou[v,alphatr],X0,maxtime,maxit,verbose],
					"ExplicitTau-Leaping",
						ExplicitTauLeapingMethod[sp,alphatr,gammatr,cou[v,alphatr],X0,maxtime,maxit,fss,verbose],
					"ImplicitTau-Leaping",
						ImplicitTauLeapingMethod[sp,alphatr,gammatr,cou[v,alphatr],vars,X0,maxtime,maxit,fss,verbose],
					"TrapezoidalTau-Leaping",
						TrapezoidalTauLeapingMethod[sp,alphatr,gammatr,cou[v,alphatr],vars,X0,maxtime,maxit,fss,verbose],
					"K-Leaping",
						error = fss/2;
						KLeapingMethod[sp,alphatr,gammatr,cou[v,alphatr],X0,maxtime,maxit,error,verbose],
					_, 
						Message[Simulation::"method",method]; 
						Return[$Failed];
				],
				Message[Simulation::"time",maxtime]; 
				Return[$Failed];
			],
			Message[Simulation::"rsarg",rates,X0]; 
			Return[$Failed];
		]
	];

Simulation[___] := (Message[Simulation::"badarg"]; $Failed)


(* ::Subsubsection::Closed:: *)
(*SimulationPlot, 2D, 3D*)


SimulationPlot::usage = "SimulationPlot[{reactions}, rates, init, maxtime, options] simulates the \
mechanism reactions and plots the result. The initial state is determined \
by init, the stopping time is specified by maxtime.\n\
Option name\t Default value\n ExternalSpecies\t {},\n PlotFunction\t \"ListLinePlot\",\n Species\t All,\n\
and other options: see Simulation... more help needed.";


listplotfunctions = {"ListPointPlot3D", "ListLinePlot", "ListPlot", "ListLogPlot"}; (* .... DiscretePlot, PerformanceGoal->"Speed"*)


Options[SimulationPlot] = { ExternalSpecies->{}, PlotFunction -> "ListLinePlot", Species -> All };

SimulationPlot::badarg = "Illegal argument of function SimulationPlot.";
SimulationPlot::species = "Wrong format or none of the given species '`1`' is in the list of internal species of reaction network {`2`,...}.";
SimulationPlot::plotfunc = "The given symbol '`1`' is not a list plotting function. Try one of the following: \
\"ListPointPlot3D\", \"ListLinePlot\", \"ListPlot\" or \"ListLogPlot\".";

(*SyntaxInformation[SimulationPlot] = {"ArgumentsPattern"->{{__},{__},{__},_,___}};*)

SimulationPlot[{reactions__}, rates_?VectorQ, X0_?VectorQ, maxtime_?NumericQ, opts___?OptionQ] := 
	Module[{ ops, plot, sp, givensp, proc, w },
			ops = Flatten[{opts, Options[SimulationPlot]}];
			If[MemberQ[listplotfunctions,(plot = PlotFunction /. ops)],
				(*Monitor*)
				plot = ToExpression[plot]; (*Flatten[{ExternalSpecies/.ops}]*)
				sp = Check[ReactionsData[{reactions},ExternalSpecies/.ops]["species"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];
				givensp = (Species /. ops) /. All->sp;
				If[(w = Flatten[Position[sp,#] &/@ Union[Flatten[{givensp}]]]) =!= {},
					proc = Check[Simulation[{reactions},rates,X0,maxtime,FilterRules[{opts},Options[Simulation]]],Return[$Failed];,
								{Simulation::volarg,Simulation::rsarg,Simulation::time,Simulation::method,Simulation::badarg}];
					(*Transpose[Map[Thread[List[#]]&,proc]][[w]], maskent kiszedve az anyagokat!!!*) (* Outer[Flatten[#2][[{1,#1+1}]]&,w,proc,1] *)
					plot[Transpose[Map[Thread[List@@#]&,proc]][[w]],FilterRules[{opts},Options[plot]]], 
					(*Show, PlotLegend*)
					Message[SimulationPlot::"species",givensp,First[{reactions}]];
					Return[$Failed];
				],
				Message[SimulationPlot::"plotfunc",plot];
				Return[$Failed];
			]
	];

SimulationPlot[___] := (Message[SimulationPlot::"badarg"]; $Failed)


SimulationPlot2D::usage = "SimulationPlot2D[{reactions}, rates, init, maxtime, options] simulates the \
mechanism reactions and plots the result in dimension 2. The initial state is determined \
by init, the stopping time is specified by maxtime.\n\
Option name\t Default value\n ExternalSpecies\t {},\n PlotFunction\t \"ListLinePlot\",\n Species\t All,\n\
and other options: see Simulation... more help needed.";


SpecSubsets[l_,w_] := Part[#,w] &/@ l;

listplot2dfunctions = {"ListPlot","ListLinePlot","ParametricPlot???"}; (**)


Options[SimulationPlot2D] = { ExternalSpecies->{}, PlotFunction->"ListPlot", Species -> All };

SimulationPlot2D::badarg = "Illegal argument of function SimulationPlot2D.";
SimulationPlot2D::species = "Wrong format or none of the given list of pair of species '`1`' is in the list of internal species of reaction network {`2`,...}.";
SimulationPlot2D::plotfunc = "The given symbol '`1`' is not a valid 2D plotting function. Try one of the following: \
\"ListPlot\"...";

(*SyntaxInformation[SimulationPlot2D] = {"ArgumentsPattern"->{{__},{__},{__},_,___}};*)

SimulationPlot2D[{reactions__}, rates_?VectorQ, X0_?VectorQ, maxtime_?NumericQ, opts___?OptionQ] := 
	Module[{ ops, plot, plotopts, sp, givensp, proc, w },
			ops = Flatten[{opts, Options[SimulationPlot2D]}];
			If[MemberQ[listplot2dfunctions,(plot = PlotFunction /. ops)],
				(*Monitor*)
				plot = ToExpression[plot];
				sp = Check[ReactionsData[{reactions},ExternalSpecies/.ops]["species"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];
				givensp = (Species /. ops) /. All->Subsets[sp,{2}];
				If[MatrixQ[givensp] && Last[Dimensions[givensp]]===2 &&(w = Flatten[Position[sp,#]&/@#] &/@ givensp) =!= {},
					proc = Check[Last /@ Simulation[{reactions},rates,X0,maxtime,FilterRules[{opts},Options[Simulation]]],Return[$Failed];,
								{Simulation::volarg,Simulation::rsarg,Simulation::time,Simulation::method,Simulation::badarg}];
					(*Print[w];*)
					Show[
						plot[SpecSubsets[proc,#]& /@ w,FilterRules[{opts},Options[plot]]]
					],
					(*Show, PlotLegend*)
					Message[SimulationPlot2D::"species",givensp,First[{reactions}]];
					Return[$Failed];
				],
				Message[SimulationPlot2D::"plotfunc",plot];
				Return[$Failed];
			]
	];

SimulationPlot2D[___] := (Message[SimulationPlot2D::"badarg"]; $Failed)


SimulationPlot3D::usage = "SimulationPlot3D[{reactions}, rates, init, maxtime, options] simulates the \
mechanism reactions and plots the result in dimension 3. The initial state is determined \
by init, the stopping time is specified by maxtime.\n\
Option name\t Default value\n ExternalSpecies\t {},\n PlotFunction\t \"ListLinePlot\",\n Species\t All,\n\
and other options: see Simulation... more help needed.";


listplot3dfunctions = {"ListPlot3D","ListPointPlot3D","ParametricPlot3D???"}; (**)


Options[SimulationPlot3D] = { ExternalSpecies->{}, PlotFunction->"ListPointPlot3D", Species -> All };

SimulationPlot3D::badarg = "Illegal argument of function SimulationPlot3D.";
SimulationPlot3D::species = "Wrong format or none of the given list of species '`1`' is in the list of internal species of reaction network {`2`,...}.";
SimulationPlot3D::plotfunc = "The given symbol '`1`' is not a valid 3D plotting function. Try one of the following: \
\"ListPlot3D\"...";

(*SyntaxInformation[SimulationPlot3D] = {"ArgumentsPattern"->{{__},{__},{__},_,___}};*)

SimulationPlot3D[{reactions__}, rates_?VectorQ, X0_?VectorQ, maxtime_?NumericQ, opts___?OptionQ]:=
	Module[{ ops, plot, plotopts, sp, givensp, proc, w },
			ops = Flatten[{opts, Options[SimulationPlot3D]}];
			If[MemberQ[listplot3dfunctions,(plot = PlotFunction /. ops)],
				plot = ToExpression[plot];
				(*Monitor*)
				sp = Check[ReactionsData[{reactions},ExternalSpecies/.ops]["species"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];
				givensp = (Species /. ops) /. All->Subsets[sp,{3}];
				If[MatrixQ[givensp] && Last[Dimensions[givensp]]===3 &&(w = Flatten[Position[sp,#]&/@#] &/@ givensp) =!= {},
					proc = Check[Last /@ Simulation[{reactions},rates,X0,maxtime,FilterRules[{opts},Options[Simulation]]],Return[$Failed];,
								{Simulation::volarg,Simulation::rsarg,Simulation::time,Simulation::method,Simulation::badarg}];
					Show[
						plot[SpecSubsets[proc,#]& /@ w,FilterRules[{opts},Options[plot]]]
					],
					Message[SimulationPlot3D::"species",givensp,First[{reactions}]];
					Return[$Failed];
				],
				Message[SimulationPlot3D::"plotfunc",plot];
				Return[$Failed];
			]
	];

SimulationPlot3D[___] := (Message[SimulationPlot3D::"badarg"]; $Failed)


(* ::Subsubsection::Closed:: *)
(*comparemethods*)


(*
CompareMethods

*)


(* ::Subsubsection::Closed:: *)
(*New methods, Langevin equation, etc.*)


(*
based on int. approx.

*)


(* ::Subsection:: *)
(**decompositions*)


(* ============== *)
(* Decompositions *)
(* ============== *)


Circular::usage = "...";

Circular[overall_?ReactionQ, steps : {__?ReactionQ}] :=
	Block[{$Messages},
    	Check[
			LinearProgramming[
				Array[-1 &, Length[steps]],
				Transpose[steps],
				{#, 0} & /@ overall
			],
      		Return[True],
      		LinearProgramming::"lpsub"
      	]; False
    ]	

CoveringDecompositionSet::usage = "CoveringDecompositionSet[overall, reactions] gives some of the decompositions \
of the reaction denoted by overall into the reaction steps denoted by steps. See also: \
Options[CoveringDecompositionSet].";
CoveringDecompositionSet::badarg = "Illegal argument of function CoveringDecompositionSet \
(should be CoveringDecompositionSet[overall_?ReactionQ, reactions_?MatrixQ]).";

ObjectiveFunction::usage = "...";
FastSelection::usage = "ObjectiveFunction \[Rule] FastSelection is an option for \
CoveringDecompositionSet to specify how the consecutive decompositions should be searched for.";
GreedySelection::usage = "ObjectiveFunction \[Rule] GreedySelection is an option for \
CoveringDecompositionSet to specify how the consecutive decompositions should be searched for.";
OriginalSelection::usage = "ObjectiveFunction \[Rule] OriginalSelection is an option for \
CoveringDecompositionSet to specify how the consecutive decompositions should be searched for. \
(The naming is purely historical and has no mathematical justification.)";

Options[CoveringDecompositionSet] = {ObjectiveFunction -> GreedySelection, Verbose -> True}

CoveringDecompositionSet[overall_?ReactionQ, reactions_?MatrixQ, opts___?OptionQ] :=
	Block[{$Messages,
			verbose = Verbose /. Flatten[{opts, Options[CoveringDecompositionSet]}],
			objective = ObjectiveFunction /. Flatten[{opts, Options[CoveringDecompositionSet]}],
			markedspecies = Array[1 &, Length[reactions]],
			decompositions = {},
			nextdecomposition,
			coefficientvector,
			mainf,
			monitormsg = "Working..."
	      },

	coefficientvector :=
		Switch[objective,
				ReactionKinetics`GreedySelection, 1 - markedspecies,
				ReactionKinetics`FastSelection, Array[0 &, Length[reactions]],
				ReactionKinetics`OriginalSelection, markedspecies
			   ];

	mainf := While[True,
(*markedspecies ?????? marked helyett*)		If[ Max[marked] === 0, Break[] ];
		nextdecomposition =
			Check[
				LinearProgramming[
					coefficientvector,
					Prepend[Transpose[reactions],markedspecies],
					Prepend[{#,0} & /@ overall, {1,1}]
				],
          			Break[],
          			LinearProgramming::"lpsnf"
          		];
		PrependTo[decompositions, nextdecomposition];
		markedspecies *= (1-Sign[nextdecomposition]);
		monitormsg = "Found " <> ToString[Length[decompositions]] <>
					  " decomposition(s) covering " <> ToString[Total[1 - markedspecies]] <>
					  " of the " <> ToString[Length[reactions]] <> " reactions."
	];
	If[verbose, Monitor[mainf, monitormsg], mainf];

	If[verbose, Print["Found " <> ToString[Length[decompositions]] <> " decomposition(s) covering " <>
		ToString[Total[1 - markedspecies]] <> " of the " <> ToString[Length[reactions]] <> " reactions."
	]];
	decompositions
]

ElementaryReactions::usage = "ElementaryReactions[amx, maxp] computes every possible \
elementary steps among the species described by the atomic matrix amx having at most maxp products. The \
second argument is optional, its default value is Infinity.";

Options[ElementaryReactions] = {Verbose -> True}

ElementaryReactions[atommatrix_?MatrixQ, maxlen : (_Integer?Positive | Infinity) : Infinity, opts___?OptionQ] :=
	Block[{$Messages, verbose = Verbose /. Flatten[{opts, Options[ElementaryReactions]}],
		molec, others, reac, i, j, fda = First[Dimensions[atommatrix]], monitorf, monitormsg = "Working...", casesdone},
		monitorf := Join[
			Flatten[
				Table[
					molec = atommatrix[[i]];
					casesdone = 2*i;
	    			others = Delete[atommatrix, {{i}}];
	    			{
						Insert[#, -1, {{i}}] & /@ Decompositions[  molec, others, maxlen, Flatten[{Verbose -> False, FilterOptions[Decompositions, opts]}]],
						Insert[#, -2, {{i}}] & /@ Decompositions[2 molec, others, maxlen, Flatten[{Verbose -> False, FilterOptions[Decompositions, opts]}]]
    				},
					{i, fda}
				],
				2
			],
			Flatten[
				Table[
    				molec = atommatrix[[i]] + atommatrix[[j]];
					casesdone = -1 + fda - i/2 + fda*i - i^2/2 + j; 
	    			others = Delete[atommatrix, {{i}, {j}}];
					Insert[#, -1, {{i}, {j-1}}] & /@ Decompositions[molec, others, maxlen, Flatten[{Verbose -> False, FilterOptions[Decompositions, opts]}]]
    				,
	    			{i, fda-1}, {j, i+1, fda}
    			],
    			2
    		]
    	];
		If[verbose, Monitor[monitorf, ToString[IntegerPart[100.*casesdone/((fda*(3 + fda))/2)]]<>"% done"], monitorf]
    ]

ElementaryReactions[___] := (Message[ElementaryStep::badarg]; $Failed)

Obligatory::usage = "Obligatory[reactions, overall] selects the reactions \
which are provably contained in each decomposition of overall into reactions.";
Obligatory::badarg = "Illegal argument of function Obligatory (see ?Obligatory)";

Obligatory[reactions : {__?ReactionQ}, overall_?ReactionQ] :=
	Transpose[{
		Flatten[Position[#, _?Positive]],
		Select[#, Positive]
	}] & [Array[
				LinearProgramming[
					Table[If[i===#, 1, 0], {i,Length[reactions]}],
					Transpose[reactions],
					{#, 0} & /@ overall
				][[#]] &,
				Length[reactions]
			]]

Obligatory[___] := (Message[Obligatory::"badarg"]; $Failed)

Omittable::usage = "...";
Omittable::badarg = "Illegal argument of function Omittable (should \
be Omittable[overall_?ReactionQ, reactions_?MatrixQ]).";

Omittable[overall_?ReactionQ, reactions_?MatrixQ, opts___?OptionQ] :=
	Block[ {optCDS = FilterOptions[CoveringDecompositionSet, opts]},
	Complement[
		Range[Length[reactions]],
		Union[Last /@ Position[CoveringDecompositionSet[overall, reactions, optCDS], _?Positive]]
	]
]

Omittable[___] := (Message[Omittable::"badarg"]; $Failed)

PrintDecomposition::usage = "PrintDecomposition[decomp, reactions] pretty prints \
the decomposition assigned to the vector decomp using the chemical equations \
given in reactions.\nPrintDecomposition[decomp, atomm, species] does the same, but \
automatically generates the chemical equations from the atomic matrix and the \
list of species. It does not work with generalized decompositions.";
PrintDecomposition::badarg = "Illegal argument of function PrintDecomposition (see ?PrintDecomposition).";

PrintDecomposition[decomp_?DecompositionQ, reactionsmatrix_?MatrixQ, species_List] :=
	PrintDecomposition[decomp, ToStep[#, species] & /@ reactionsmatrix]

PrintDecomposition[decomp_?DecompositionQ, reactions_List] :=
	GridBox[
		Transpose[{
			"(" <> ToString[#] <> ")" & /@ #,
			decomp[[#]],
			reactions[[#]]
		}] & [Flatten[Position[decomp, _?Positive]]],

		ColumnAlignments -> {Right, Right, Left},
		ColumnSpacings -> {2, 1.2}
	] // DisplayForm
(*------- DisplayForm[*]*)
PrintDecomposition[___] := (Message[PrintDecomposition::"badarg"]; $Failed)

SMaux[minimals_List, {}] := minimals
SMaux[minimals_List, {next_List, rest___List}] :=
	SMaux[Join[minimals, Select[next, MinimalQ[#, minimals] &]], {rest}]

SelectMinimalDecompositions::usage = "SelectMinimalDecompositions[vlist] selects those vectors \
from vlist, which are minimal with respect to the \"componentwise less or equal\" partial order. \
Note, that using this function to select minimal decompositions will not necessarily give truly \
minimal decompositions, unless all decompositions up to a certain length are given in vlist.";
SelectMinimalDecompositions::badarg = "Illegal argument of function SelectMinimalDecompositions \
(should be SelectMinimalDecompositions[vlist : {__?DecompositionQ}]).";

SelectMinimalDecompositions[vlist : {__?DecompositionQ}] :=
	SMaux[First[#], Rest[#]] &[
		Map[Last,
			Split[Sort[{Total[#], #} & /@ vlist, First[#1] < First[#2] &], First[#1] === First[#2] &],
			{2}]
		]

SelectMinimalDecompositions[___] := (Message[SelectMinimalDecompositions::"badarg"]; $Failed)


(* ============== *)
(* Decompositions *)
(* ============== *)

Options[CDPartitions] = {Verbose -> True, Filter -> MinimalDecompositions}

CDPartitions[overall_?ReactionQ, steps : {__?ReactionQ}, stepno_, opts___?OptionQ] :=
	Module[
		{verbose, filter,
		 vecs = Append[steps, -overall],
		 im = IdentityMatrix[Length[steps] + 1],
		 solutions = {},
		 stack, product, tmpsol, stepstogo,
		 chooseable, properoriented,
		 x, verboseStepnum = 0, verboseMaxdepth = 0,
		 monitormsg = "Working...", monitorf},

		{verbose, filter} = {Verbose, Filter} /. Flatten[{opts, Options[Decompositions]}];
		If[filter === All, Message[Decompositions::cdwall]; Return[$Failed]];

		stack = Array[
					{steps[[#]] - overall,
					 im[[#]] + im[[Length[vecs]]],
					 stepno,
					 Range[#]} &,
					Length[steps]
				];

		monitorf := While[Length[stack] =!= 0,

			{product, tmpsol, stepstogo, chooseable} = First[stack];
			++verboseStepnum;
			monitormsg = "Running... " <> ToString[Length[solutions]] <>
								  " decomposition(s) found in " <> ToString[verboseStepnum] <>
								  " steps.";
			verboseMaxdepth = Max[verboseMaxdepth, Total[tmpsol]];

			Which[
				(*Total[tmpsol] > stepno ||*) stepstogo === 0 || ! MinimalQ[tmpsol, solutions], stack = Rest[stack],
				ZeroVectorQ[product], PrependTo[solutions, tmpsol];
									  stack = Rest[stack],
				True, properoriented = Select[chooseable, product.vecs[[#]] <= -(product.product)/stepstogo && (* this is the improvement *)
														  product.vecs[[#]] < 0 &&
														  MinimalQ[tmpsol + im[[#]], solutions] &];
					  stack = Join[
								Map[Function[x, {product + vecs[[x]],
            									 tmpsol + im[[x]],
												 stepstogo - 1,
												 Complement[chooseable, Select[properoriented, # > x &]]
												}],
        							properoriented
								],
								Rest[stack]];
			];
		];
		If[verbose, Monitor[monitorf, monitormsg], monitorf];

		If[verbose,
			Print[Length[solutions], " decomposition(s) found in ", verboseStepnum,
			 " steps.\nMaximum search depth was ", verboseMaxdepth, ", the longest decomposition consists of ",
        	 Max[Total /@ solutions] - 1, " step(s)."
			]
		];
		Most /@ solutions
	]

(* ----------------------------------------------------------------------------------------------------- *)

UnrestrictedSolutions[m_?MatrixQ, b_List] :=
	Module[{ls, ns},
		Check[
			ls = LinearSolve[m, b],
			Return[{}],
			LinearSolve::"nosol"
		];
    	If [(ns = NullSpace[m]) === {}, ls, ls + Transpose[ns].(x /@ Range[Length[ns]])]
    ]

Options[LPPartitions] = {Filter -> All, Preprocess -> False, Verbose -> True}

LPPartitions[overall_?ReactionQ, steps : {__?ReactionQ}, maxlen_, opts___?OptionQ] :=
	Block[{verbose, filter, preprocess, generalsol, vars, solvedLPs = 0, constraints},

		{verbose, filter, preprocess} = {Verbose, Filter, Preprocess} /. Flatten[{opts, Options[LPPartitions]}];

		If[maxlen === Infinity,
			If[Circular[overall, steps] && filter === All, Message[Decompositions::unbnd]; Return[$Failed]];
			If[Circular[overall, steps], Message[Decompositions::wrmet]; Return[$Failed]]
		];

		If[ (generalsol = UnrestrictedSolutions[Transpose[steps], overall]) === {}, Return[{}] ];

		vars = Union[Cases[generalsol, x[_], Infinity]];
		If[vars === {}, Return[
							If[And @@ (IntegerQ[#] && NonNegative[#] & /@ generalsol) && Total[generalsol] <= maxlen,
          					   {generalsol},
							   {}
							]]
    	];

    	constraints = If[maxlen =!= Infinity, Append[generalsol, maxlen - Total[generalsol]], generalsol];

    	(* Check whether the polyhedron of the relaxed problem is empty *)
    	(* Saves a lot of breath and time if it is, good for nothing if not *)
    	If[MyMinimize[x[1], constraints, vars] === "no solution",
    		If[verbose,	Print["0 solution(s) found with 1 solved LPs."]];
			Return[{}]
		];

		constraintsols = (* FindAllSolutionsWithPreprocessing[constraints, vars, verbose]; *)
			If[preprocess,
				FindAllSolutionsWithPreprocessing[constraints, vars, verbose],
				FindAllSolutions[constraints, vars]
			];

		If[verbose,
			Print[ToString[Length[constraintsols]] <> " solution(s) found with " <>
				  ToString[solvedLPs + If[preprocess, 2 Length[vars], 0]] <> " solved LPs."]
		];
    	generalsol /. Thread[vars -> #] & /@ constraintsols

    ]

(* Az uj valtozathoz *)
FindAllSolutionsWithPreprocessing[conlist_List, vars_List, verbose_] :=
	Module[{bounds, bnd, bnd2, filtered = {}, f1 = 0, f2 = 0, newf2 = True, unfiltered, newconlist = conlist,
	filteredsolutions, i, lvar = Length[vars], steprec, monitorf, monitormsg = "Working..."},

		(* If[verbose, Print["Preprocessing..."]]; *)

	monitorf = While[newf2 && Length[filtered] < lvar,
		newf2 = False;
		For[i = 1, i <= lvar, i++,
			monitormsg = ToString[i-1] <> "/" <> ToString[lvar] <> " variables examined, " <>
											ToString[f1] <> " + " <> ToString[f2] <> " variables filtered by Rule #1 and #2.";
			If[! And @@ (IntegerQ /@ Select[newconlist, NumberQ]), Return[{}]];
			If[MemberQ[filtered, i], Continue[]];

			(* steprec = GCD @@ (Cases[newconlist, a_. + b_. x[i] /; NumberQ[a] :> b, {1}]);
			If[steprec == 0, Print[{i, conlist, newconlist}]]; *)

			solvedLPs += 2;
			Check[
(*x ????*)				bnd = #[x[i], newconlist, vars] & /@ {MyMinimize, MyMaximize},
				Return[{}],
				LinearProgramming::lpsnf
				(* ?? Miert nem lep ki soha? Akkor a kovetkezo mar nem is kene *)
			];
			If[ !MatchQ[bnd, {_?NumberQ, _?NumberQ}], Return[{}] ];

			If[ SameQ @@ bnd,

				(* Print["Variable ", i, " is FILTERED by Rule no. 1."]; *)
				AppendTo[filtered, i];
				(* AppendTo[bounds, bnd]; *)
				bounds[i] = bnd;
				newconlist = newconlist /. x[i] -> First[bnd];
				++f1;
				If[! And @@ (IntegerQ /@ Select[newconlist, NumberQ]), Return[{}]],

				If[ Length[
						bnd2 = Intersection @@
								(Cases[newconlist, a_. + b_. x[i] /; NumberQ[a] :> Range[
																					NestWhile[# + 1/b &, -a/b, If[Sign[b]>0, # < First[bnd], # > Last[bnd]] &],
																					If[Sign[b]>0, Last[bnd], First[bnd]],
																					1/b
																				   ],
								 {1}])
					] === 1,

					(* Print["Variable ", i, " is FILTERED by Rule no. 2."]; *)
					(*Print[newconlist];
					Print[bnd];
					Print[bnd2];*)

				 	(* Ceiling[First[bnd] steprec] + 1 > steprec Last[bnd] *)

					bnd2 = First[bnd2];
					AppendTo[filtered, i];
					(* bnd2 = Ceiling[First[bnd] steprec] / steprec; *)
					(* AppendTo[bounds, {bnd2, bnd2}]; *)
					bounds[i] = {bnd2, bnd2};
					newconlist = newconlist /. x[i] -> bnd2;
					(* Print[newconlist]; *)
					++f2;
					If[! And @@ (IntegerQ /@ Select[newconlist, NumberQ]), Return[{}]];
					newf2 = True,

						(*Print["Variable ", i, " is NOT filtered."]; *)
						(* AppendTo[bounds, bnd]; *)
						bounds[i] = bnd;
				]
			]
		]
	];

	If[verbose, Monitor[monitorf, monitormsg], monitorf];

		unfiltered = Complement[Range[lvar], filtered];

		If[verbose, Print[Length[filtered], " out of ", lvar, " variable(s) filtered out in preprocessing."]
		];

		If[unfiltered === {},
			Return[ {Array[First[bounds[#]] &, {lvar}]} ],

			(* Print[newconlist /. ReactionKinetics`Private`x -> Global`x];
			Print[First /@ bounds /@ unfiltered];
			Print[Last /@ bounds /@ unfiltered]; *)
			filteredsolutions = FindAllSolutions[
									newconlist,
									x /@ unfiltered,
									(* bounds[[unfiltered, 2]] *)
									First /@ bounds /@ unfiltered,
									Last /@ bounds /@ unfiltered
								];
			(* Print["fsols: ", filteredsolutions]; *)
		];

		(* Return *)
		Array[If[MemberQ[unfiltered, #], 0, (* bounds[[#, 1]] *) First[bounds[#]]] &, {lvar}] +
		VectorFromPoints[lvar, unfiltered, #] & /@ filteredsolutions
	]

(* UJ VALTOZAT, ALSO-FELSO KORLATOKKAL, LESZAMLALASSAL *)
FindAllSolutions[conlist_List, vars : {var_, restvars___}, lowerbounds_List, upperbounds_List] :=
	Module[{lower, upper, newcons, a, b, step, offset},

		{step, offset} = First[Sort[
			Cases[conlist, a_. + b_. var /; NumberQ[a] :> {Abs[1/b], -a/b}, {1}],
			First[#1] > First[#2] &
		]];

		(*	step   = 1 / GCD @@ (Cases[conlist, a_. + b_. var /; NumberQ[a] :> b, {1}]);
			offset = First[Cases[conlist, a_. + b_. var /; NumberQ[a] :> -a/b, {1}, 1]]; *)

		(* Print[{var, step, offset}]; *)
		(* Print[conlist /. ReactionKinetics`Private`x -> Global`x]; *)
		(* If[vars =!= Union[Cases[conlist,x[_],Infinity]], Print[vars]]; *)
		(* Print[solvedLPs]; *)
		(* Print[{Length[vars], step, offset}]; *)

		(* solvedLPs += 2;
		Check[
			lower = MyCeiling[MyMinimize[var, conlist, vars] - offset, step] + offset;
			upper = MyFloor[MyMaximize[var, conlist, vars] - offset, step] + offset,
			Return[{}],
			LinearProgramming::lpsnf
		];*)

		lower = MyCeiling[First[lowerbounds] - offset, step] + offset;
		upper = MyFloor[First[upperbounds] - offset, step] + offset;
		While[Or @@ Negative[conlist /. var -> lower] && (lower <= upper), lower += step];
		While[Or @@ Negative[conlist /. var -> upper] && (lower <= upper), upper -= step];

		(* Return *)
		If[lower > upper,
			{},

			If[Length[vars] === 1,

				(* Print[Length[Select[Range[lower, upper, step], And @@ IntegerQ /@ (conlist /. var -> #) &]]]; *)
				Transpose[{Select[Range[lower, upper, step], And @@ IntegerQ /@ (conlist /. var -> #) &]}],

				Join @@ Table[
					newcons = conlist /. var -> i;
					If[ Not[And @@ IntegerQ /@ Select[newcons, NumberQ]],
					    {},
						Prepend[#, i] & /@ FindAllSolutions[
												DeleteCases[newcons, _Integer],
												{restvars},
												Rest[lowerbounds],
												Rest[upperbounds]
										   ]
					], {i, lower, upper, step}
				] (* Join@@Table*)
			] (* If length *)
		] (* If lower *)
	] (* Module *)

(* REGEBBI, NINCS ELOFElDOLGOZAS *)
FindAllSolutions[conlist_List, vars : {var_, restvars___}] :=
	Module[{lower, upper, newcons, a, b, step, offset},

		{step, offset} = First[Sort[
			Cases[conlist, a_. + b_. var /; NumberQ[a] :> {Abs[1/b], -a/b}, {1}],
			First[#1] > First[#2] &
		]];

		(* Print[{var, step, offset}]; *)
		(* Print[conlist /. ReactionKinetics`Private`x -> Global`x]; *)
		(* If[vars =!= Union[Cases[conlist,x[_],Infinity]], Print[vars]]; *)
		(* Print[solvedLPs]; *)
		(* Print[{Length[vars], step, offset}]; *)

		solvedLPs += 2;
		Check[
			lower = MyCeiling[MyMinimize[var, conlist, vars] - offset, step] + offset;
			upper = MyFloor[MyMaximize[var, conlist, vars] - offset, step] + offset,
			Return[{}],
			LinearProgramming::lpsnf
		];
		(*
		lower = MyCeiling[First[lowerbounds] - offset, step] + offset;
		upper = MyFloor[First[upperbounds] - offset, step] + offset;
		While[Or @@ Negative[conlist /. var -> lower] && (lower <= upper), lower += step];
		While[Or @@ Negative[conlist /. var -> upper] && (lower <= upper), upper -= step];
		*)
		(* Return *)
		If[lower > upper,
			{},
			If[Length[vars] === 1,
				(* Print[Length[Select[Range[lower, upper, step], And @@ IntegerQ /@ (conlist /. var -> #) &]]]; *)
				Transpose[{Select[Range[lower, upper, step], And @@ IntegerQ /@ (conlist /. var -> #) &]}],
				Join @@ Table[
					newcons = conlist /. var -> i;
					If[ Not[And @@ IntegerQ /@ Select[newcons, NumberQ]],
					    {},
						Prepend[#, i] & /@ FindAllSolutions[
												DeleteCases[newcons, _Integer],
												{restvars}
										   ]
					], {i, lower, upper, step}
				] (* Join@@Table*)
			] (* If length *)
		] (* If lower *)
	] (* Module *)

(* ----------------------------------------------------------------------------------------------------- *)

Decompositions::usage = "Decompositions[overall, steps] gives the decompositions of the reaction \
denoted by overall into the reaction steps denoted by steps. See also: Options[Decompositions].";

Decompositions::badarg = "Illegal argument of function Decompositions (should \
be\n\tDecompositions[overall_?ReactionQ, steps : {__?ReactionQ}, maxlen : \
(_Integer?Positive | Infinity) : Infinity]). See ??Decompositions for further \
instructions.";

Decompositions::cdwall = "Method -> ContejeanDevie can yield the minimal \
solutions only. Choose another method, or use the option Filter -> \
MinimalDecompositions or MinimalSupportDecompositions.";

Decompositions::unbnd = "The number of the decompositions is unbounded. Give \
upper bound for solution length, or use the option Filter -> \
MinimalDecompositions or MinimalSupportDecompositions and select the \
appropriate method.";

Decompositions::wrmet = "Incompatible method and filter option (perhaps due to the circularity of\
the network of steps). (Ez fel\[EAcute]r egy bad command or file name-vel, egyel\[OTilde]re... [??])";

(*Methods of Decompositions*)
ContejeanDevie::usage = "...";
LPBased::usage = "...";
(*Options of Decompositions*)
Filter::usage = "...";
Method::usage = "...";
Preprocess::usage = "...";
Verbose::usage = "Verbose is an option of multiple functions such as Concentrations, \
CoveringDecompositionSet etc. that makes the functions give more detailed information \
during the calculations. (The results are unaffected.) Default value may vary. Enable \
it using Verbose -> True.";
(*Filter*)
OptimalPartitions::usage = "...";
MinimalDecompositions::usage = "...";
(*All*)

Options[Decompositions] = {Filter -> MinimalDecompositions, Method -> ContejeanDevie, Preprocess -> True, Verbose -> True}

Decompositions[overall_?ReactionQ, steps : {__?ReactionQ},
			   maxlen : (_Integer?Positive | Infinity) : Infinity, opts___?OptionQ] :=

	Block[{method, $RecursionLimit = Infinity},

		{method} = {Method} /. Flatten[{opts, Options[Decompositions]}];

		Switch[method,
				LPBased, LPPartitions[overall, steps, maxlen, FilterOptions[LPPartitions, opts]],
				ContejeanDevie, CDPartitions[overall, steps, maxlen, FilterOptions[CDPartitions, opts]]
				(* , SubsetEnumeration, SEPartitions[overall, steps, maxlen, FilterOptions[SEPartitions, opts]] *)
		]
    ];

Decompositions[___] := (Message[Decompositions::"badarg"]; Return[$Failed];)


(* ::Subsection::Closed:: *)
(*palettes, reaction rates notebook*)


OpenReactionKineticsPalette::usage = "\!\(\*\nStyleBox[\"OpenReactionKinetics\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\)\
\!\(\*\nStyleBox[\"[\",\nFontFamily->\"Courier\"]\)\!\(\*\nStyleBox[\"]\",\nFontFamily->\"Courier\"]\) \
provides a palette involving the basic symbols and notations of reaction kinetics.";

OpenReactionKineticsNamesPalette::usage = "\!\(\*\nStyleBox[\"OpenReactionKineticsNamesPalette\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\)\
\!\(\*\nStyleBox[\"[\",\nFontFamily->\"Courier\"]\)\!\(\*\nStyleBox[\"]\",\nFontFamily->\"Courier\"]\) \
provides a palette involving the basic names and notions of reaction kinetics.";


OpenReactionKineticsPalette::badarg = "Illegal argument of function OpenReactionKineticsPalette. \
It should not be given any arguments.";

OpenReactionKineticsNamesPalette::badarg = "Illegal argument of function OpenReactionKineticsNamesPalette. \
It should not be given any arguments.";


SyntaxInformation[OpenReactionKineticsPalette] = { "ArgumentsPattern"->{} };

OpenReactionKineticsPalette[] :=
    NotebookPut[Notebook[{
        Cell[BoxData[GridBox[
		{
		 {ButtonBox["\[Alpha]"], ButtonBox["\[Beta]"], ButtonBox["\[Gamma]"], ButtonBox["\[Delta]"], ButtonBox["\[Theta]"], ButtonBox["\[Lambda]"]},
		 {ButtonBox["\"M\""], ButtonBox["\"E\""], ButtonBox["\"R\""], ButtonBox["(\[Placeholder])"], ButtonBox["{\[Placeholder]}"], ButtonBox["[\[Placeholder]]"]},
		 {ButtonBox["\[ShortRightArrow]"], ButtonBox["\[RightArrow]"], ButtonBox["\[LongRightArrow]"], ButtonBox["\[ShortLeftArrow]"],ButtonBox["\[LeftArrow]"], ButtonBox["\[LongLeftArrow]"]},
		 {ButtonBox["\[LeftRightArrow]"], ButtonBox["\[LongLeftRightArrow]"], ButtonBox["\[Equilibrium]"], ButtonBox["\[ReverseEquilibrium]"], ButtonBox["\[RightArrowLeftArrow]"], ButtonBox["\[LeftArrowRightArrow]"]}, 
		 {ButtonBox["\[DoubleRightArrow]"], ButtonBox["\[DoubleLongRightArrow]"], ButtonBox["\[DoubleLeftArrow]"], ButtonBox["\[DoubleLongLeftArrow]"], ButtonBox["\[DoubleLeftRightArrow]"], ButtonBox["\[DoubleLongLeftRightArrow]"]},
			(*ButtonBox["\[Placeholder]\!\(\*OverscriptBox[\"\[RightArrow]\", \"\[Placeholder]\"]\)\[Placeholder]"]},*)
		 {ButtonBox[SubscriptBox["\[Placeholder]", "\[Placeholder]"]], ButtonBox[SuperscriptBox["\[Placeholder]", "\[Placeholder]"]], ButtonBox[SubsuperscriptBox["\[Placeholder]", "\[Placeholder]", "\[Placeholder]"]],
		  ButtonBox["\!\(\*SubscriptBox[\"\[InvisiblePrefixScriptBase]\", \"\[Placeholder]\"]\)\[Placeholder]"], ButtonBox["\!\(\*SuperscriptBox[\"\[InvisiblePrefixScriptBase]\", \"\[Placeholder]\"]\)\[Placeholder]"], ButtonBox["\!\(\*SubsuperscriptBox[\"\", \"\[Placeholder]\", \"\[Placeholder]\"]\)\[Placeholder]"]}
		}, 
		RowSpacings -> 1.5, ColumnSpacings -> 1.5]],
        FontFamily -> "TimesNewRoman", FontWeight -> "Normal", FontSize -> 18]}, (*Bold*)
		Enabled->True,
        CellMargins -> {{0, 0}, {0, 0}},
		WindowMargins->{{Automatic,0},{Automatic,0}},
        Editable -> True, (*False*)
        ShowCellBracket -> False,
        WindowElements -> {},
        WindowFrame -> "Palette",
        WindowFrameElements -> "CloseBox",
        WindowSize -> {Fit, Fit},
        WindowToolbars -> {},
        WindowClickSelect -> False,
        WindowTitle -> "ReactionKinetics1 Palette",
		WindowFloating->True]
];

OpenReactionKineticsPalette[__] := (Message[OpenReactionKineticsPalette::"badarg"]; $Failed);

SyntaxInformation[OpenReactionKineticsNamesPalette] = { "ArgumentsPattern"->{} };

OpenReactionKineticsNamesPalette[] := 
	NotebookPut[Notebook[{
        Cell[BoxData[GridBox[
		{
			{ButtonBox["\"atommatrix\""]},
			{ButtonBox["\"atoms\""]},
			{ButtonBox["\"complexes\""]},
			{ButtonBox["\"deficiency\""]},
			{ButtonBox["\"externalspecies\""]},
			{ButtonBox["\"fhjgraph\""]},
			{ButtonBox["\"fhjconnectedcomponents\""]},
			{ButtonBox["\"fhjstronglyconnectedcomponents\""]},
			{ButtonBox["\"reactionorder\""]},
			{ButtonBox["\"reactionsteps\""]},
			{ButtonBox["\"species\""]},
			{ButtonBox["\"variables\""]},
			{ButtonBox["\"volpertgraph\""]}
		(*  {ButtonBox["Atoms"]},
			{ButtonBox["Models"]},
			{ButtonBox["Reactions"]}  *)
		},RowSpacings ->.9, ColumnSpacings ->1.]],
        FontFamily -> "TimesNewRoman", FontWeight -> "Normal", FontSize -> 18]}, (*Bold*)
		Enabled -> True,
        CellMargins -> {{0, 0}, {0, 0}},
		WindowMargins->{{Automatic,0},{110,Automatic}},
        Editable -> False,
        ShowCellBracket -> False,
        WindowElements -> {},
        WindowFrame -> "Palette",
        WindowFrameElements -> "CloseBox",
        WindowSize -> {Fit, Fit},
        WindowToolbars -> {},
        WindowClickSelect -> False,
        WindowTitle -> "ReactionKinetics1 Names",
		WindowFloating->True,WindowElements->{"VerticalScrollBar"}]]

OpenReactionKineticsNamesPalette[__] := (Message[OpenReactionKineticsNamesPalette::"badarg"]; $Failed);


ReactionRatesNotebook::usage = 
"\!\(\*\nStyleBox[\"ReactionRatesNotebook\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\)\
\!\(\*\nStyleBox[\"[\",\nFontFamily->\"Courier\"]\){reactions}\
\!\(\*\nStyleBox[\",\",\nFontFamily->\"Courier\"]\)externalspecies\
\!\(\*\nStyleBox[\"]\",\nFontFamily->\"Courier\"]\) \
opens a new notebook with a command that assigns reaction rates to the given reaction network. \
The result is a list of reaction rates or functions in the same order as they have \
to be written in every command that uses reaction rates as parameter \
(This is the order the reactionsteps are enumerated in the reaction network).";


ReactionRatesNotebook::"badarg" = "Illegal argument of function ReactionRatesNotebook.";


ReactionKinetics1`Private`myReactionRatesCounter = 0;

Options[ReactionRatesNotebook] := {ExternalSpecies->{}};

SyntaxInformation[ReactionRatesNotebook]={"ArgumentsPattern"->{__,OptionsPattern[]}};

ReactionRatesNotebook[{reactions__}] := ReactionRatesNotebook[{reactions},{}];

ReactionRatesNotebook[{reactions__},externals:(_?VectorQ|OptionsPattern[])] :=
	(++ReactionKinetics1`Private`myReactionRatesCounter;
	NotebookPut[Notebook[
		{
		Cell[TextData[{
			"Fill in the following table in every placeholder with the appropriate reaction rate constants or functions as well.\nEvaluate the cell when finished."}], 
			"Text", FontSize -> 15, FontColor -> RGBColor[0,0, 0], FontWeight -> Bold
			],
		Cell[BoxData[RowBox[{RowBox[{"myReactionRates" <> ToString[ReactionKinetics1`Private`myReactionRatesCounter], "=",
				RowBox[{"Last", "/@",GridBox[
					Thread[List @@ Through[List[ToBoxes/@First[#]&,ToBoxes/@Array[Placeholder[{##}]&,Last[#]]&][Check[ReactionsData[{reactions},externals]["reactionsteps","R"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}]]]],
						RowLines -> True, RowSpacings -> 2.5, ColumnSpacings -> 3, ColumnAlignments -> {Center, Right}]
				}]}]}]], "Input"
			]
		},
		WindowSize -> {720, 500},
		WindowTitle -> "Reaction Rates Notebook " <> ToString[ReactionKinetics1`Private`myReactionRatesCounter]
		]
	];)

ReactionRatesNotebook[___] := (Message[ReactionRatesNotebook::"badarg"]; $Failed)


(* ::Subsection::Closed:: *)
(*welcoming words*)


If[Length[Names["FLAGS`ECHOLOAD"]]>0,
	ReactionKinetics1`Private`msgflag=ToExpression["FLAGS`ECHOLOAD"],
	ReactionKinetics1`Private`msgflag=True;
];

If[ReactionKinetics1`Private`msgflag,
	Print[Style["\!\(\*
		StyleBox[\"ReactionKinetics\",\nFontFamily->\"Courier New\"]\) Version "<> $ReactionKineticsVersionNumber <>" using Mathematica Version "<>
		$Version<>" (Version "<>ToString[$VersionNumber]<>", Release "<>ToString[$ReleaseNumber]<>")"<>
		" loaded "<>DateString[{"Day"," ","MonthName", " ", "Year", " at ", "Hour24",":", "Minute", " ", "TimeZone"}]
		<>"\nGNU General Public License (GPLv3) Terms Apply. ",FontWeight-> Bold,FontColor->Black]];
];

If[ReactionKinetics1`Private`msgflag,
	Print[Style["Please report reaction kinetics issues or problems at\n",FontWeight-> Bold,FontColor->Black], 
			Hyperlink["nagyal@math.bme.hu", "mailto:nagyal@math.bme.hu"],Style[", "],
			Hyperlink["papp.david@gmail.com","mailto:papp.david@gmail.com"], Style[", "],
			Hyperlink["jtoth@math.bme.hu","mailto:jtoth@math.bme.hu"]]
];

FLAGS`ECHOLOAD = ReactionKinetics1`Private`msgflag;
Remove[ReactionKinetics1`Private`msgflag];


(* ::Subsection::Closed:: *)
(*end*)


End[ ];


If[ $ReactionKineticsPackageLoaded =!= True, 
	If[ Length[Notebooks["ReactionKinetics1 Palette"]] === 0,
			OpenReactionKineticsPalette[];
	];
	If[ Length[Notebooks["ReactionKinetics1 Names"]] === 0, 
			OpenReactionKineticsNamesPalette[];
	];
];

$ReactionKineticsPackageLoaded = True;


SetAttributes[
	{
		$ReactionKinetics1,
		$ReactionKineticsPackageLoaded,
		$ReactionKineticsVersionNumber,
		$Maybe
	},
	{Protected,Locked}
];


Protect[
	Global`aExplodator,
	Global`bExplodator,
	Global`fOregonator,
	Global`c,
	Global`t,
	Global`k,
	Global`\[DoubleStruckCapitalE],
	Global`g,
	Global`\[CapitalPi],
	Global`P,
	Global`x,
	Global`z,
	CombinatorialMoments,
	ComplexColours,
	Conditions,
	ExternalSpecies,
	FixedStepSize,
	FormattedOutput,
	GeneratedRateCoefficient,
	MaxIteration,
	Method,
	PlotFunction,
	Positivity,
	Species,
	StronglyConnectedComponentsColours,
	TimeLimit,
	Verbose,
	Volume
];


SetAttributes[
	{
		Models,
		Reactions,
		BioModels,
		BioReactions,
		GetProblem,
		(**)
		ToCanonicalForm,
		ToReversible,
		ReactionsData,
		ReversibleQ,
		WeaklyReversibleQ,
		ShowFHJGraph,
		ShowVolpertGraph,
		Atoms,
		ToAtomMatrix,
		AtomConservingQ,
		AtomsQ,
		FromAtomMatrix,
		DetailedBalanced,
		(**)
		RightHandSide,
		DeterministicModel,
		Concentrations,
		StationaryPoints,
		MassConservation,
		EigenvaluesJacobianMatrix,
		AbsoluteRobustness,
		(**)
		VolpertIndexing,
		(**)
		MassConservationQ,
		ChemicalMasterEquation,
		StationaryProbabilityDensityEquation,
		ProbabilityGeneratingFunctionEquation,
		SolveProbabilityGeneratingFunctionEquation,
		ExpectedvalueProcess,
		DependencyGraph,
		ShowDependencyGraph,
		Simulation,
		SimulationPlot,
		SimulationPlot2D,
		SimulationPlot3D,
		(**)
		OpenReactionKineticsPalette,
		OpenReactionKineticsNamesPalette,
		ReactionRatesNotebook
	},
	{Protected,ReadProtected}
];


(*$Epilog := NotebookClose/@Join[Notebooks["ReactionKinetics1 Palette"],Notebooks["ReactionKinetics1 Names"]];*)


EndPackage[ ];
