(* ::Package:: *)

(* ::Subsection::Closed:: *)
(*prre*)


(* :Title: ReactionKinetics *)

(* :Authors: Attila Laszlo NAGY, David PAPP and 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: 

*)

(* :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`",
	{
	 (*"PhysicalConstants`", (*ezzel egyutt a Units` is bejon*)*)
	 (*"Utilities`FilterOptions`",*)
	 "MultivariateStatistics`"
	}]


If[$VersionNumber<6, 
	Off[DiscreteMath`Combinatorica`Star]; 
	Needs["Graphics`"];
	Needs["Graphics`Legend`"];
	Needs["LinearAlgebra`MatrixManipulation`"];
	Needs["Utilities`"];
	(*Needs["DiscreteMath`GraphPlot`"];*)
	On[DiscreteMath`Combinatorica`Star];
]; 

(* 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`"];
	Needs["PlotLegends`"];
	Needs["MathSBML`"];
	Needs["EquationTrekker`"];
];
*)

If[$VersionNumber < 8,
	Off[General::compat];
	Needs["Combinatorica`"];
	Needs["GraphUtilities`"];
    On[General::compat];
]; 


Unprotect[
	$ReactionKinetics1,
	$ReactionKineticsPackageLoaded,
	$ReactionKineticsVersionNumber,
	$Maybe,
	(**)
	Global`aExplodator,(*Explodator coeffs*)
	Global`bExplodator,(*Explodator coeffs*)
	Global`fOregonator,(*Explodator coeffs*)
	(**)
	ReactionImport,
	(**)
	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`T,(*temperature*)
(*	Global`x,(*state vector*)*)
	Global`X, (*species*)
	(**)
	Global`z,(*gen func arg symbol*)
	CombinatorialMoments, (**)
	FixedStepSize,
	MaxIteration,
	Method, (* StiffNessSwitching, slow- fast variables *)
	Species,
	Verbose,
	Volume, (**)
	MassConservationQ, (**)
	ChemicalMasterEquation,
	StationaryProbabilityDensityEquation,
	ProbabilityGeneratingFunctionEquation,
	SolveProbabilityGeneratingFunctionEquation,
	MomentEquations,
	DependencyGraph,
	ShowDependencyGraph,
	Simulation, (**)
	SimulationPlot,
	SimulationPlot2D,
	SimulationPlot3D,
	CompareModels,
	(**)
	OpenReactionKineticsPalette,
	OpenReactionKineticsNamesPalette,
	ReactionRatesNotebook,

	(**)
	CHEMKINExport,
	CHEMKINImport, 
	FromStoichiometry, 
	DeleteAutocatalysis,
	ReversibleFHJRepresentation,
	Memo,
	EdgeLabeled,
	Highlight,
	Bipartitized,
	Indexed, 
	Numbered,
	SubgraphHighlight,
	Side,
	FilterReactions,
	MinFHJWeaklyConnectedComponents, 
	MinFHJStronglyConnectedComponents, 
	MaxFHJWeaklyConnectedComponents, 
	MaxFHJStronglyConnectedComponents,
	VolpertSpecies
];


Unprotect[
	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,
(**)
	VolpertIndexingWithMatrices,
	ComponentwiseSmallerEqualQ, 
	MinimalQ, 
	MyMinimize,
	MyMaximize,
	ZeroVectorQ, 
	MyFloor, 
	MyCeiling, 
	SMaux, 
	CDPartitions, 
	UnrestrictedSolutions, 
	FindAllSolutionsWithPreprocessing, 
	LPPartitions, 
	FindAllSolutions
];


(* ================================================================ *)
(* 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 [06-01-2013]";


$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*)


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

MyFilterOptions[f_, opts___]:= 
(*	If[$VersionNumber<6,
			Utilities`FilterOptions`FilterOptions[f, opts],*)
			Sequence @@ FilterRules[Flatten[{opts}], Options[f]];


(*filteroptions[f_, opts___] := Sequence @@ FilterRules[{opts}, Options@f];*)


(* ::Subsection::Closed:: *)
(*kinetic 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",
	"Birth-Death process",
	"Briggs-Rauscher",
	"Brusselator",
	"Chapman cycle",
	"Colquhoun",
	"Decomposition",
	"DeYoung-Keizer",
	"Dimerization",
	"Eigen",
	"\[CapitalEAcute]rdi-Ropolyi",
	"Explodator",
	"Feinberg-Horn I",
	"Feinberg-Horn II",
	"FKN mechanism",
	"Goldbeter",
	"Hudson-R\[ODoubleDot]ssler",
	"Huxel",
	"Hyver",
	"Inflow",
	"Ivanova",
	"Kliemann",
	"Leonard-Reichl",
	"Lotka",
	"Lotka-Volterra",
(*	"Lotka-Volterra Cyclic",*)
	"Michaelis-Menten",
	"Mole",
	"Oregonator",
	"Outflow",
	"Robertson",
	"Schl\[ODoubleDot]gl",
	"Triangle reaction",
(*	"Triangle reaction Reversible",*)
	"Vaiman",
	"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,+]"]
	   }*)
	
	(* Birth-Death process *)
{"X"->2"X", "X"->"A"},

	(* 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"*)

	(*Chapman cycle *)
{"O2"->"O"+"O", "O2"+"O"+"M"->"O3"+"M", "O3"->"O2"+"O", "O"+"O3"->"O2"+"O2"},
(*externals: "M"*)

	(* 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"},

	(* Goldbeter *)
{"M"\[Equilibrium]"0"\[Equilibrium]"C","0"\[Equilibrium]"X","C"\[RightArrow]"0"},

	(* Hudson-R\[ODoubleDot]ssler *) 
{"P"->"A","Q"->"B","A"+"B"->2"B","B"->"R","A"\[LeftRightArrow]"C"}, 
(*externals: "P","Q","R"*)
	
	(* Huxel, special Lotka-Volterra *)
{"X"\[Equilibrium]2"X", "X"+"Y"->2"Y","Y"->"0"},

	(* 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"},

	(* Mole *)
{"X" + "Y" \[RightArrow] 2 "X" + 2 "Y", "X" \[LeftRightArrow] "0" \[LeftRightArrow] "Y"},

	(* 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"},
*)

	(* Vaiman *)
{"O2" + 2 "K" -> 2 "KO", "H2S" + "K" -> "KH2S", "O2" + 2 "KH2S" -> 2 "H2O" + "K" + "KS2", 
    "H2S" + "KO" -> "H2O" + "KS", "KO" + "KH2S" -> "H2O" + "KS" + "K", 2 "KS" -> "K" + "KS2", "KS2" -> "S2" + "K"},

	(* 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::Closed:: *)
(**fitting to standards*)


(*
MyMemberQ[x_,s_] := MemberQ[x,s] || Or@@(If[StringQ[#], StringCount[#,s]==1, False] &/@ x);

MyTotal[x_] := 
	Module[{ pos },

		Total[

			If[ DigitQ[StringTake[#,1]],

				pos = Position[DigitQ /@ Characters[#], False, 1, 1][[1,1]]-1;

				ToExpression[StringTake[#,pos]]*StringDrop[#,pos],

				#

			] & /@ x

		]

	];

FromStringToStepWithPars[x_]:=
	Module[{ll, arrhs},

		ll = x /.{Longest[a__?StringQ],b__} :> Length[{a}]; (*If[StringQ[x[[2]]], ll=3;, ll=1;];*)

		If[Length[x]-ll < 3, ll = ll - 1;];

		arrhs = Drop[x, ll];

		Append[
			Flatten[
					StringSplit[
							StringReplace[StringJoin[x[[1;;ll]]],
							{y_~~"+AR":>y,y_~~"(+AR)":>y,y_~~"+HE":>y,y_~~"(+HE)":>y,y_~~"(+M)":>y,y_~~"+M":>y,y_~~"(+Ar)":>y,y_~~"(+O2)":>y,y_~~"(+H2O)":>y,y_~~"(+N2)":>y}
							],
					{"=>"->0,"="->1,"<=>"->1,"+"}
					] /. {{a__,0,b__}:>{RightArrow[MyTotal[{a}],MyTotal[{b}]]},{a__,1,b__}:>{Equilibrium@@(MyTotal/@{{a},{b}})}}
																						(*Equilibrium@@(MyTotal/@Sort[Sort/@{{a},{b}}])*)												
			],
			(arrhs /. Thread[{"!","!me"}->Sequence[]] /. y_?StringQ:>ToExpression[StringReplace[y,{"!":>"",z_~~"+":>z~~"*10^"}]])[[1;;3]]]
];
(*Remarks: felkialtojeles reakciok ma'r nincsenek benne, ill. Konnov fele felkialtojeles seb.egy.hatok kiveve: ez az [[1;;3]]*)
(*a vegen: elobb a felkialtojel kezelese, utana ha meg maradt benne string egyszeruen delete, Konnov-ba meg sok a szemet*)
(*egyszeru program arra, hogy ha zarojelben van megadva egy anyag reagenskent akkor azt vegyuk ki, de ne kezzel!!*)
(*1.000+136 ---> "!"->""*)
*)


(*
(* MEMO IMPORT FUNCTIONS *)

MyMemoImport[filename_] := MyMemoImport[filename] = 
		Check[Import[filename], Return[$Failed];, {Import::nffil,Import::chtype}];

MyImport[filename_,data_] := MyImport[filename,data] =
	Module[{ m, mm = filename, x, mf, arr },

		m = MyMemoImport[mm];

		If[m===$Failed, Message[CHEMKINImport::"bimp"]; Abort[];];

		If[ Position[m, "ELEM"] =!= {}, x = "ELEM"; , 
				If[Position[m, "ELEMENTS"] =!= {}, x = "ELEMENTS";, x = "" ]];

		If[ x =!= "",
			mf = Flatten[m];

			Switch[data,
				"chemkinelements", ( mf /. {a___, x, Shortest[b___],"END",c___} :> (Flatten[{b}] /. "!" -> Sequence[] ) ),

				"chemkinspecies", ( mf /. {a___,"SPECIES",Shortest[b___],"END",c___} :> (Flatten[{b}] /. "!" -> Sequence[] ) ),
					(*/. {a___,{b___,"SPECIES",c___},Shortest[d___],{e___,"END",f___},g___}:>{d}*)
				"chemkinarrhenius" | "chemkinreactions",
					(arr = FromStringToStepWithPars /@ 
						((m /. (
									{a___,{b___,"REACTIONS",c___},Shortest[d___],{e___,"END",f___},g___} :>  
									Select[{d}, MyMemberQ[#,"="]&&FreeQ[#,"DUPLICATE"|"LOW"|"TROE"|"REV"|"END"|"Deleted"]&&StringFreeQ[First[#],"!"]&]
								))
									/. {a__, "!", b___} :> {a} 
					 ));
					If[data === "chemkinarrhenius", arr, DeleteDuplicates[First /@ arr]]
					
			],
		 
			Switch[data,
				"chemkinarrhenius" | "chemkinreactions",
					(arr = FromStringToStepWithPars /@ 
								Select[m /. {} -> Sequence[], StringQ[First[#]] && Not[StringFreeQ[ToString[First[#]],"="]]&];
					 );
					If[data === "chemkinarrhenius", arr, DeleteDuplicates[First /@ arr]],
				_,
				Message[CHEMKINImport::misdat,data];
				Return[{}]
			]
		]

	];
*)


MyTotal[x_] := 
	Module[{ pos },

		Total[
			If[ DigitQ[StringTake[#,1]],

				pos = Position[DigitQ /@ Characters[#], False, 1, 1][[1,1]]-1;

				ToExpression[StringTake[#,pos]]*StringDrop[#,pos],

				#
			] & /@ x
		]
	];

EqualDependentQ := StringCount[#,"="] === 1&;

SS[x_] := 
	Module[{ r, read },
		read = (If[
					MemberQ[
						r=Quiet[
							ReadList[
								StringToStream[
									StringReplace[#, "."~~Longest[a__?DigitQ]~~x1_?(#=!="E"&&DigitQ[#]&)~~"+"~~b_?DigitQ :> "."~~a~~x1~~"E+"~~b]
								],
								Number]
						   ], $Failed],
						#,
					    Sequence@@r
				]) &/@ x;

				read /. {xx___,Longest[yy__?StringQ],Longest[zz__?NumberQ],uu___} :> 
					{
						StringSplit[
							StringReplace[
								StringJoin[yy],
								{"(+"~~Shortest[___]~~")":>""}
							],
							{"=>"->0,"="->1,"<=>"->1,"+"}
						] /. {{a__, 0, b__} :> {RightArrow[MyTotal[{a}],MyTotal[{b}]]}, {a__, 1, b__} :> {Equilibrium@@(MyTotal/@{{a},{b}})}}, {zz}
					}
];


(* MEMO IMPORT FUNCTIONS *)

MyMemoImport[filename_] := MyMemoImport[filename] = 
		Check[Import[filename,"Text"], Return[$Failed];, {Import::nffil,Import::chtype}];


MyImport[filename_, data_] := MyImport[filename, data] =
	Module[{ m, mm = filename, x, y, z, arr },

		m = MyMemoImport[mm]; (*text-kent beolvasva*)

		If[m===$Failed, Message[CHEMKINImport::"bimp"]; Return[$Failed];];

		m = StringReplace[m,"!"~~Shortest[___]~~"\n"->"\n"]; (*szemet-kommentek kitorolve*)

		If[ StringCases[m, "ELEMENTS"] =!= {}, x = "ELEMENTS"; , 
				If[StringCases[m, "ELEMS"] =!= {}, x = "ELEMS";, x = ""; ]];

		If[ StringCases[m, "SPECIES"] =!= {}, y = "SPECIES"; , y = ""; ];

		If[ StringCases[m, "REACTIONS"] =!= {}, z = "REACTIONS"; , z = ""; ];

		If[ z === "", 
				Message[CHEMKINImport::misdat,"REACTIONS"];		
				Return[{}],

				arr = Cases[
						StringSplit[
							First[StringCases[m,z~~Shortest[a___]~~"END":>a]],
								"\n"],
						x_?EqualDependentQ :> (SS[StringSplit[x]] /. Thread[{"M","HE","AR"}:>0] /. {"BASE",___}:>Sequence[])
					  ];
		];

		Switch[ {x =!= "",y =!= ""},

			{True, True},   Switch[data,
								"chemkinelements", Flatten[StringSplit[StringCases[m, x~~Shortest[a___]~~"END":>a], {Whitespace,"\n"}]],

								"chemkinspecies", Flatten[StringSplit[StringCases[m, y~~Shortest[a___]~~"END":>a], {Whitespace,"\n"}]],

								"chemkinarrhenius" | "chemkinreactions",
									If[data === "chemkinarrhenius", {#[[1,1]], #[[2]]}& /@ arr, DeleteDuplicates[Flatten[First /@ arr]]]
								],
				
			{True, False},  
							Switch[data,
								"chemkinelements", Flatten[StringSplit[StringCases[m, x~~Shortest[a___]~~"END":>a], {Whitespace,"\n"}]],

								"chemkinarrhenius" | "chemkinreactions",
									If[data === "chemkinarrhenius", {#[[1,1]], #[[2]]}& /@ arr, DeleteDuplicates[Flatten[First /@ arr]]],
								_,
									Message[CHEMKINImport::misdat,"SPECIES"];  
									Return[{}]
							],

			{False, True},   
							Switch[data,
								"chemkinspecies", Flatten[StringSplit[StringCases[m, y~~Shortest[a___]~~"END":>a], {Whitespace,"\n"}]],

								"chemkinarrhenius" | "chemkinreactions",
									If[data === "chemkinarrhenius", {#[[1,1]], #[[2]]}& /@ arr, DeleteDuplicates[Flatten[First /@ arr]]],
								_, 
									Message[CHEMKINImport::misdat,"ELEMENTS"]; 
									Return[{}]
							],

			{False, False}, 
							Switch[data,
								"chemkinarrhenius" | "chemkinreactions",
									If[data === "chemkinarrhenius", {#[[1,1]], #[[2]]}& /@ arr, DeleteDuplicates[Flatten[First /@ arr]]],
								_,
									Message[CHEMKINImport::misdat,"ELEMENTS"];  
									Message[CHEMKINImport::misdat,"SPECIES"];
									Return[{}]
							]
		]
	];


PropertiesCHEMKINImport = 
	{
		"chemkinelements","chemkinspecies","chemkinreactions","chemkinarrhenius"
	};


Options[CHEMKINImport] := {Memo -> False};


CHEMKINImport::usage = "CHEMKINImport[\"filename\"] imports all the relevant informations relating to the reaction mechanism given in the file. \
For further options and specific data see ...[data].";


CHEMKINImport::illarg = "The argument of CHEMKINImport needs a string with a file name of which extension is .dat!";
CHEMKINImport::"badname" = "At least one of the arguments '`1`' is a non-identified property. For the possibilities use the command \
\!\(\*\nStyleBox[\"CHEMKINImport\",\nFontFamily->\"Courier\",\nFontColor->RGBColor[0, 0, 1]]\)[\"Properties\"].";
CHEMKINImport::bimp = "Some problem occured around the import procedure. For possible further issues, see the function Import."; 
CHEMKINImport::misdat = "Missing DATA: the CHEMKIN file to be imported does not contain data for '`1`'.";
CHEMKINImport::badarg = "Illegal argument of function CHEMKINImport.";

(*SyntaxInformation[CHEMKINImport]={"ArgumentsPattern"->{__}};*)

CHEMKINImport[filename_?StringQ]["Properties"] := PropertiesCHEMKINImport;

CHEMKINImport[filename_?StringQ][] := Thread[PropertiesCHEMKINImport->CHEMKINImport[filename]["All"]];

CHEMKINImport[filename_?StringQ]["All"] := CHEMKINImport[filename][PropertiesCHEMKINImport];

CHEMKINImport[filename_?StringQ][data__?StringQ] := 
	(
	If[ Last[StringSplit[StringTrim[filename],"."]] =!= "dat", 
		Message[CHEMKINImport::"illarg"];
		Return[$Failed];
	];
	If[Nand @@ Map[MemberQ[PropertiesCHEMKINImport,#]&,{data}],
		Message[CHEMKINImport::"badname",data];
		Return[$Failed];
	];
	(*Catch[*)
		Part[Map[(*Check[*)MyImport[filename, #](*,Throw[] ...]*)&, {data}], Length[{data}] /. x_/;x>1 :> (1;;x)
		]
	(*]*)
	);

CHEMKINImport[filename_?StringQ][data__] := 
	CHEMKINImport[filename][Sequence @@ (ToString /@ Flatten[{data}])];

(*CHEMKINImport[files__?StringQ] := CHEMKINImport /@ Flatten[{files}];*)

Format[CHEMKINImport[filename_?StringQ]] := 
	DynamicModule[{x, chi},
					x = Dynamic[Refresh[Round[Clock[Infinity]],UpdateInterval->1]];
					Monitor[chi = CHEMKINImport[filename][];,
									Column[{Row[{"CHEMKINImport is now importing and 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[];
					chi(*FormatString[{reactions},{externals}]*)
	];

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


(**********************************)


CHEMKINExport::usage = "CHEMKINExport[\"filename\",\"reactions\"] exports the reaction mechanism given in the second argument in CHEMKIN format. \
For further options and specific data see ...";


CHEMKINExport::illarg = "The first argument of CHEMKINExport needs a string with a file name (optionally with path) of which extension is .dat.";
CHEMKINExport::fexist = "The given file(name) in which to be exported the reaction mechanism have already existed and will be overwritten.";
CHEMKINExport::wrgformat = "The given list of reactions has wrong fromat. It is recommended to be a list of reaction steps with their arrhenius parameters.";
CHEMKINExport::badarg = "Illegal argument of function CHEMKINExport.";

CHEMKINExport[filename_?StringQ, reactions_List] := 
	Module[{filestr, path, i, r, str, ns, n, sq, step, rstep, arrh, w, nform, b = 0},
		
		If[ Last[StringSplit[StringTrim[filename],"."]] =!= "dat", 
			Message[CHEMKINExport::"illarg"];
			Return[$Failed];
		];

		If[FileNameDepth[filename] <= 1 || StringMatchQ[filename, "*:/*"], 
			path = filename;, path = FileNameJoin[{Directory[],filename}];
		];

		If[FileExistsQ[path], Message[CHEMKINExport::"fexist"]; ]; (*Do you want to delete or rename DeleteFile[path]*)
		
		filestr = Check[OpenWrite[path], Return[$Failed];, {General::"aofil", OpenWrite::"noopen"}];

		If[ filestr === $Failed, Return[$Failed]; ];

		r = Length[reactions];		

		Do[
			step = reactions[[i]];			

			If[Head[step]===List,
				If[Length[step]===2,
					If[Length[Last[step]]=!=3,
						Message[CHEMKINExport::wrgformat];
						Return[$Failed];
					];,
					Message[CHEMKINExport::wrgformat];
					Return[$Failed];
				];,
				Message[CHEMKINExport::wrgformat];
				Return[$Failed];
			];

			rstep = First[step];
			arrh = Last[step];

			If[
				Length[
					Join@@(Position[rstep,#]&/@
								{Rule,RightArrow,LongRightArrow,Equilibrium,ReverseEquilibrium,LeftRightArrow,LongLeftRightArrow,LongLeftArrow,LeftArrow}
							)] =!= 1,
					Message[CHEMKINExport::wrgformat];
					Return[$Failed];,
					If[
						MemberQ[{LongLeftArrow, LeftArrow}, Head[rstep]],
							w = ToString[rstep[[2]]->rstep[[1]]];,
							w = ToString[rstep];
					];
					(*ns = If[StringQ[#], NumericQ[ToExpression[#,InputForm]], NumericQ[#]]& /@ arrh;*)
					ns = NumericQ /@ arrh;
					If[And @@ ns,
						(*
						sq = MapIndexed[{#1, First[#2]}&, StringQ /@ arrh];
						str = StringToStream[""];

						n = Map[
							If[First[#], 
									Close[str]; 
									str = StringToStream[arrh[[Last[#]]]]; 
									N[Read[str, Number]], 
									N[arrh[[Last[#]]]]
							]&, sq
						];
						Close[str];
						*)

						nform = Riffle[StringReplace[ToString[FortranForm[#]],"e"~~x_:>If[x==="-","E"~~x,"E+"~~x]] &/@ arrh,"\t"];

						(*nform = ScientificForm[#,{Infinity,10},NumberFormat->(StringJoin[{#1,"E",If[#3==="","+0",If[StringTake[#3,1]==="-",#3,"+"~~#3]],"\t"}]&)] &/@ n;*)
						(*nform = NumberForm[#,NumberFormat->(StringJoin[{#1,"E+",If[#3==="","0",#3],"\t"}]&)] &/@ n;*)
						WriteString[filestr, ToString[i]<>"\n"];
						WriteString[filestr,
							StringReplace[w, {" ":>"","->"->"=>","\[RightArrow]"->"=>","\[LongRightArrow]"->"=>","\[Equilibrium]":>"<=>","\[ReverseEquilibrium]":>"<=>","\[LeftRightArrow]":>"<=>","\[LongLeftRightArrow]":>"<=>"}]<>"\t\t"<>StringJoin[nform]<>"\n"]; (*StringJoin[ToString /@ nform]*)
						b = i;,

						Message[CHEMKINExport::wrgformat];
						Return[$Failed];
					];
					
			];,
			{i, 1, r}
		];

		Close[filestr];

		If[ b === r, 
			Print["Data have been written successfully to the file (path): " <> filename];,
			Return[$Failed];
		];

	];

(* {FileNameSetter[Dynamic[f]],Dynamic[f]} *)

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


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


(* ::Subsubsection::Closed:: *)
(*ReactionsData and some related functions*)


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[]
			]
			]
];


complextospecies[cplxs_] := DeleteDuplicates[
								((Flatten[{cplxs} /. Plus:>List] 
										/. r_*A_?StringQ :> A /. A_?StringQ*r_ :> A)
										/. r_*A_?SpeciesQ :> A /. A_?SpeciesQ*r_ :> A)
										/. A_?StringQ :> A /. A_?SpeciesQ :> A /. r_*"0":>"0"
							];
									(*collect, symbolic coefficient, multip coeff*)
									(*Union*)

complextocoefflist[cplxs_] := Block[{spec},
								spec = complextospecies[cplxs];
								Thread[List[spec, Coefficient[cplxs, spec]]]
							 ];

subgraphedges[graph_,vtceslist_] := 
			If[ Length[#] === 1, {{}, #}, Join[{EdgeList[Subgraph[graph, #]] /. DirectedEdge->Rule}, {#}]] & /@ vtceslist;



Rdata[{reactions__},externals__] := Rdata[{reactions},externals] = 
	Module[
			{ reacs,
			  m, builtin, exs, exsrule, 
			  steps, steprarrow, steprule, fhjgraph, complexes, vgggraph, vgraph, vggraphwzerocomplex, sp, 
			  indexed, a, b, e, nloc, lloc, sloc,
			  lsp, lrs, alpha, beta, gamma, fhjweakvertices, fhjweakcomponents, 
			  fhjstrongvertices, fhjstrongcomponents,
			  fhjterminalstrongvertices, fhjterminalstrongcomponents
			},

			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*)

				exsrule = Thread[exs -> 0];
				(**)
				(* MAIN *)

				steps = ReactionsToList[reacs] /. (0->"0"); (*"0" - zero complex*)

				steprarrow = DeleteDuplicates[steps /. exsrule /. (0->"0")];

				steprule = steprarrow /. RightArrow -> Rule;

				(*Print[steprule];*)
				fhjgraph = Graph[steprule];

				complexes = DeleteDuplicates[Flatten[steprarrow /. RightArrow -> List]]; (*Sort*)
				(*VertexList[fhjgraph];*)

				vgggraph = Flatten[(Thread /@ {#->complextocoefflist[Last[#]], complextocoefflist[First[#]]->#}) &/@ steprarrow];
				vgraph = vgggraph /. Rule[A_RightArrow,B_]:>{Rule[A, First[B]], Last[B]} /. Rule[A_,B_RightArrow]:>{Rule[First[A], B], Last[A]};
				(*Print[vgraph];*)

				sp = complextospecies[complexes] /. "0" -> Sequence[]; (*komplexek csak belso anyagfajta erejeig egyertelmuek*)

				indexed = Flatten[MapIndexed[#1->First[#2]&, #]& /@ {sp, steprarrow}];
				(*Print[indexed];*)

				{lsp, lrs} = Length /@ {sp, steprule};

				vggraphwzerocomplex = vgraph /. {Rule["0",y_RightArrow],z_} :> Sequence[] /. {Rule[x_RightArrow,"0"],z_} :> Sequence[];

				a = Cases[vggraphwzerocomplex, {Rule[x_,y_RightArrow],z_}:> Rule[{x,y} /. indexed, z]];
				b = Cases[vggraphwzerocomplex, {Rule[x_RightArrow,y_],z_}:> Rule[{y,x} /. indexed, z]];
				(*Print[{Join[{a,b}],lsp,lrs}];*)
				alpha = SparseArray[a, {lsp, lrs}];
				beta = SparseArray[b, {lsp, lrs}];
				gamma = beta-alpha;
				(*SparseArray[Join[a,b] //. {x___,Rule[{y_,z_},k_],u___,Rule[{y_,z_},l_],v___}:>{x,Rule[{y,z},k+l],u,v}, {lsp,lrs}];*)
				(*Print[gamma];*)

				fhjweakvertices = ConnectedComponents[UndirectedGraph[fhjgraph]];
				fhjweakcomponents = subgraphedges[fhjgraph, fhjweakvertices];

				fhjstrongvertices = ConnectedComponents[fhjgraph];
				fhjstrongcomponents = subgraphedges[fhjgraph, fhjstrongvertices];

(*				volpertgraph = Graph[First /@ vgraph];
				volpertbasis = Flatten[Map[If[
											VertexInDegree[volpertgraph,#]===0, 
											#, {}
										    ]&, 
											Append[sp,"0"]
									      ] /. "0" :> {}
							         ];
*)
(*				fhjstronglength = Length[fhjstrongvertices];
				fhjstronggraph = Graph[
									DeleteDuplicates[
										steprarrow /. 
											Sort[Flatten[MapThread[Thread[#1 -> #2]&, {fhjstrongvertices, Range[fhjstronglength]}]], 
													Length[First[#1]]>Length[First[#2]]&
												]
									] /. (RightArrow[y_,y_] :> RightArrow[y,"\[Infinity]"]) /. RightArrow -> Rule];
				(*Print[fhjstronggraph];*)
				volpertbasis = Map[If[
									VertexInDegree[fhjstronggraph,#]===0, 
									Append[fhjstrongcomponents[[#]],Flatten[(First /@ complextocoefflist[#])& /@ fhjstrongcomponents[[#,2]]]], {{}}
									]&,
									(VertexList[fhjstronggraph] /. "\[Infinity]":>Sequence[])
								] /. {{}} :> Sequence[];
*)		

				fhjterminalstrongvertices = Map[
												Complement[
													Union @@ Map[Cases[steprarrow, RightArrow[#,y_]:>y]&, #], #
												]==={} /. True-># /. False->Sequence[]&, fhjstrongvertices
											];
				fhjterminalstrongcomponents = subgraphedges[fhjgraph, fhjterminalstrongvertices];

				Dispatch[{

					"species" -> sp,

					"M" -> lsp,

					"externalspecies" -> (e = exs /. "0" -> Sequence[]),
					(*e = If[MemberQ[complexes,"0"], exs, Rest[exs]]*)

					"E" -> Length[e],

					"reactionsteps" -> steps, (*steprule*)

					"R" -> lrs,

					"complexes" -> complexes,

					"fhjgraph" -> steprule, 
					(*FromOrderedPairs[steprule /. Rule[c1_,c2_] :> {Position[complexes, c1][[1, 1]], Position[complexes, c2][[1, 1]]}]*)

					"fhjweaklyconnectedcomponents" -> fhjweakcomponents,
					(*(fhjweak = WeakComponents[steprule])*)

					"fhjstronglyconnectedcomponents" -> fhjstrongcomponents,
					(*(fhjstrong = StrongComponents[steprule])*)

					"fhjterminalstronglyconnectedcomponents" -> fhjterminalstrongcomponents,

					(*"volpertbasis" -> {},*)

					"n" -> (nloc = Length[complexes]),

					"l" -> (lloc = Length[fhjweakvertices]),

					"s" -> (sloc = MatrixRank[gamma]),
					(*MatrixRank[gamma = First[Differences[ab = (Normal/@{a,b})]]]*)

					"deficiency" -> "\[Delta]=n-l-s="<>ToString[nloc]<>"-"<>ToString[lloc]<>"-"<>ToString[sloc]<>"="<>ToString[nloc-lloc-sloc],

					"\[Alpha]" -> alpha,

					"\[Beta]" -> beta,

					"\[Gamma]" -> gamma,

					"reactionorders" -> (Total /@ Transpose[alpha]),

					"variables" -> (Subscript[Global`c, #] & /@ sp),

					"volpertgraph" -> {vgraph, Drop[indexed, lsp]}

				}]
			]
	];


PropertiesReactionsData = 
	{
		"species","M","externalspecies","E","reactionsteps","R","complexes",
		"fhjgraph","fhjweaklyconnectedcomponents","fhjstronglyconnectedcomponents","fhjterminalstronglyconnectedcomponents",
		"n","l","s","deficiency","\[Alpha]","\[Beta]","\[Gamma]","reactionorders","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 fhjweaklyconnectedcomponents:..\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[{},externals:(_?VectorQ|OptionsPattern[])][] := {};

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[{},externals:(_?VectorQ|OptionsPattern[])]["Properties"] := PropertiesReactionsData;

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

ReactionsData[{},externals:(_?VectorQ|OptionsPattern[])]["All"] := {};

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

ReactionsData[{},externals:(_?VectorQ|OptionsPattern[])][data__?StringQ] := 
	(
	If[Nand @@ Map[MemberQ[PropertiesReactionsData,#]&,{data}],
		Message[ReactionsData::"badname",data];
	];
	Return[{}]
	);

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[{},externals:(_?VectorQ|OptionsPattern[])][data__] := ReactionsData[{}][Sequence @@ ToString /@ Flatten[{data}]];

ReactionsData[{reactions__},externals:(_?VectorQ|OptionsPattern[])][data__] := 
		Block[{d, pos, rdata, sol = {data}},

				d = Flatten[sol];

				rdata = ReactionsData[{reactions},externals][Sequence @@ ToString /@ d];
				(*Check[, Return[$Failed], {ReactionsData::wrreac,ReactionsData::badarg,ReactionsData::badname}]];*)

				If[Length[d]===1, rdata = {rdata}];				

				pos = Flatten[Position[sol, #]& /@ d, 1];

				Do[Part[sol, Sequence@@pos[[n]]] = rdata[[n]], {n, Length[d]}];
				
				If[ Length[sol] === 1,
					Flatten[sol, 1],
					sol
				]
		];

(*
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])


(* ::Subsubsection::Closed:: *)
(*ToCanonicalForm, ToReversible, ReversibleFHJRepresentation, ReversibleQ, WeaklyReversibleQ*)


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)


ReversibleFHJRepresentation::usage = "ReversibleFHJRepresentation[reactions,externalspecies] results the Feinberg--Horn--Jackson \
graph representation of equilibriums and rules of the given chemical mechanism.";

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


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

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

ReversibleFHJRepresentation[{reactions__},externals:(_?VectorQ|OptionsPattern[])] :=
	ReplaceRepeated[
		Check[ReactionsData[{reactions},externals]["fhjgraph"],Abort[];,{ReactionsData::wrreac,ReactionsData::badarg}],
		{x___, Rule[y_,z_], t___, Rule[z_,y_], u___} :> {x, Equilibrium[y,z], t, u}
	];

(*avagy: Complement[EdgeList[g], EdgeList[ReverseGraph[g]]], where g = FHJ Graph object*)

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


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

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


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

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

ReversibleQ[{reactions__},externals:(_?VectorQ|OptionsPattern[])] := 
		FreeQ[ReversibleFHJRepresentation[{reactions},externals],_Rule];
		(*SymmetricQ[Check[ReactionsData[{reactions}]["fhjgraph"],Abort[];,{ReactionsData::wrreac,ReactionsData::badarg}]];*)

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


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

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


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

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

(*ez lenyegeben arra megy, hogy egyetlen erosen osszefuggo komponensbol all-e az fhj graf*)
WeaklyReversibleQ[{reactions__},externals:(_?VectorQ|OptionsPattern[])] := 
		(Length[Check[ReactionsData[{reactions},externals]["fhjstronglyconnectedcomponents"],Abort[];,{ReactionsData::wrreac,ReactionsData::badarg}]] === 1);
		(*SymmetricQ[TransitiveClosure[Check[ReactionsData[{reactions}]["fhjgraph"],Abort[];,{ReactionsData::wrreac,ReactionsData::badarg}]]];*)

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


(* ::Subsubsection::Closed:: *)
(**FilterReactions, FromStoichiometry, DeleteAutocatalysis*)


FilterReactions::usage = "FilterReactions[]... \
....";


FilterReactions::badsp = "The given species are not contained in the given mechanism.";
FilterReactions::badopt = "The given option Side has optional values only: Left, Right, or All.";
FilterReactions::badarg = "Illegal argument of function FilterReactions.";

Options[FilterReactions]:={ExternalSpecies -> {}, Side -> All};

FilterReactions[{reactions__},specs_?VectorQ,opts : OptionsPattern[]]:=
	Module[{ options, rdata, fhj, sp, spp, side },

		options = Flatten[{opts, Options[FilterReactions]}];

		rdata = Check[ReactionsData[{reactions},ExternalSpecies /. options]["fhjgraph","species"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];
		If[rdata === $Failed, Return[$Failed];];
		
		{fhj, sp} = rdata;
			
		If[ (spp = Intersection[specs,sp]) =!= {},

			If[ MemberQ[{Left,Right,All}, (side = Side /. options)],
					
				Switch[side,
					Left, 
						Flatten[If[Intersection[complextospecies[First[#]],specs]=!={}, #, {}] &/@ fhj],
					Right, 
						Flatten[If[Intersection[complextospecies[Last[#]],specs]=!={}, #, {}] &/@ fhj],
					All, 
						Flatten[If[Intersection[complextospecies[# /. Rule -> Plus],specs]=!={}, #, {}] &/@ fhj]
				],
				Message[FilterReactions::"badopt"]; 
				Return[$Failed]
			],
			Message[FilterReactions::"badsp"];
			Return[$Failed]
		]
		
	];

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


FromStoichiometry::usage = "FromStoichiometry... \
....";


FromStoichiometry::dimmx = "The type of the two matrix are not the same, must have the same size.";
FromStoichiometry::vars = "The number of species (third argument) is not in correspondence with the row dimensions of alpha (first argument) and beta (second argument).";
FromStoichiometry::badarg = "Illegal argument of function FromStoichiometry.";

FromStoichiometry[alpha_?MatrixQ,beta_?MatrixQ] := 
	FromStoichiometry[alpha, beta, Array[Global`X[#]&,Length[alpha]]];

FromStoichiometry[alpha_?MatrixQ,beta_?MatrixQ,specs_?VectorQ]:=
	Module[{ m, r },
		{m, r} = Dimensions[alpha];

		If[ {m, r} === Dimensions[beta],

			If[ m === Length[specs],

				Thread[specs.alpha -> specs.beta], (*Rightarrow?*)

				 Message[FromStoichiometry::"vars"];
				 Return[$Failed];
			],

		 Message[FromStoichiometry::"dimmx"];
		 Return[$Failed];
		]
	];

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


SetAttributes[{PositivePart,NegativePart}, Listable];

PositivePart[x_] := Max[x,0];

NegativePart[x_] := -Min[x,0];


DeleteAutocatalysis::usage = "DeleteAutocatalysis... \
....";


DeleteAutocatalysis::dimmx = "The type of the two matrix are not the same, must have the same size.";
DeleteAutocatalysis::vars = "The number of species (third argument) is not in correspondence with the row dimensions of alpha (first argument) and beta (second argument).";
DeleteAutocatalysis::badarg = "Illegal argument of function DeleteAutocatalysis.";

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

DeleteAutocatalysis[{reactions__}, OptionsPattern[]] := 
	DeleteAutocatalysis @@ Check[ReactionsData[{reactions},OptionValue[ExternalSpecies]]["\[Alpha]","\[Beta]","species"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];

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

DeleteAutocatalysis[alpha_?MatrixQ,beta_?MatrixQ] := 
	DeleteAutocatalysis[alpha, beta, Array[Global`X[#]&,Length[alpha]]];

DeleteAutocatalysis[alpha_?MatrixQ,beta_?MatrixQ,specs_?VectorQ]:=
	Module[{ gamma, m, r, reacs, revreacs },
		{m, r} = Dimensions[alpha];

		If[ {m, r} === Dimensions[beta],

			If[ m === Length[specs],
			gamma = beta-alpha;

			reacs = DeleteDuplicates[Thread[specs.NegativePart[gamma] -> specs.PositivePart[gamma]]];
			(*Thread[Rule@@(species.#[gamma]&/@{NegativePart,PositivePart})]]*)
			revreacs = Reverse /@ reacs;

			Join[Complement[reacs,revreacs] /. Rule -> RightArrow, Union[Sort /@ Intersection[reacs, revreacs]] /. Rule -> Equilibrium],
		
			 Message[DeleteAutocatalysis::"vars"];
			 Return[$Failed];
			],

		 Message[DeleteAutocatalysis::"dimmx"];
		 Return[$Failed];
		]
	];

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


(* ::Subsubsection::Closed:: *)
(*MinFHJWeaklyConnectedComponents, MinFHJStronglyConnectedComponents, MaxFHJWeaklyConnectedComponents, MaxFHJStronglyConnectedComponents*)


MinFHJWeaklyConnectedComponents::usage = "MinFHJWeaklyConnectedComponents... \
....";


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

MinFHJWeaklyConnectedComponents[{reactions__}, OptionsPattern[]] := 
	Module[{fhjwcc, lengthvtxes, minvtxes},

		fhjwcc = Check[ReactionsData[{reactions},OptionValue[ExternalSpecies]]["fhjweaklyconnectedcomponents"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];

		lengthvtxes = Length[Last[#]] &/@ fhjwcc;

		minvtxes = Min[lengthvtxes];

		(*fhjwcc[[Flatten[Position[lengthvtxes, minvtxes]]]]*)
		Select[fhjwcc, Length[Last[#]]===minvtxes &]

	];

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


MinFHJStronglyConnectedComponents::usage = "MinFHJStronglyConnectedComponents... \
....";


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

MinFHJStronglyConnectedComponents[{reactions__}, OptionsPattern[]] := 
	Module[{fhjscc, lengthvtxes, minvtxes},

		fhjscc = Check[ReactionsData[{reactions},OptionValue[ExternalSpecies]]["fhjstronglyconnectedcomponents"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];

		lengthvtxes = Length[Last[#]] &/@ fhjscc;

		minvtxes = Min[lengthvtxes];

		(*fhjscc[[Flatten[Position[lengthvtxes, minvtxes]]]]*)
		Select[fhjscc, Length[Last[#]]===minvtxes &]

	];

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


MaxFHJWeaklyConnectedComponents::usage = "MaxFHJWeaklyConnectedComponents... \
....";


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

MaxFHJWeaklyConnectedComponents[{reactions__}, OptionsPattern[]] := 
	Module[{fhjwcc, lengthvtxes, maxvtxes},

		fhjwcc = Check[ReactionsData[{reactions},OptionValue[ExternalSpecies]]["fhjweaklyconnectedcomponents"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];

		lengthvtxes = Length[Last[#]] &/@ fhjwcc;

		maxvtxes = Max[lengthvtxes];

		(*fhjwcc[[Flatten[Position[lengthvtxes, maxvtxes]]]]*)
		Select[fhjwcc, Length[Last[#]]===maxvtxes &]

	];

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


MaxFHJStronglyConnectedComponents::usage = "MaxFHJStronglyConnectedComponents... \
....";


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

MaxFHJStronglyConnectedComponents[{reactions__}, OptionsPattern[]] := 
	Module[{fhjscc, lengthvtxes, maxvtxes},

		fhjscc = Check[ReactionsData[{reactions},OptionValue[ExternalSpecies]]["fhjstronglyconnectedcomponents"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];

		lengthvtxes = Length[Last[#]] &/@ fhjscc;

		maxvtxes = Max[lengthvtxes];

		(*fhjscc[[Flatten[Position[lengthvtxes, maxvtxes]]]]*)
		Select[fhjscc, Length[Last[#]]===maxvtxes &]

	];

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


(* ::Subsubsection::Closed:: *)
(*ShowFHJGraph*)


GraphPlotFunctionQ := MemberQ[{"Graph","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 Graph, GraphPlot, GraphPlot3D, LayeredGraphPlot 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]:={ExternalSpecies -> {}, ComplexColours -> {}, PlotFunction -> "GraphPlot", StronglyConnectedComponentsColours -> {}, Numbered -> False};

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

ShowFHJGraph[{reactions__}, rates_List:{}, opts___?OptionQ] := 
	Module[{ cc, externals, rdata, exs, 
			 n, r, fhj, cxs, fhjscc, graphfunc, nmd,
			 nopt, c, sc, rscc, rs 
		   },

			cc = FilterRules[{opts},Options[ShowFHJGraph]];
			externals = FilterRules[{opts},ExternalSpecies];

			rdata = Check[ReactionsData[{reactions},Flatten[externals]]["n","R","fhjgraph","complexes","fhjstronglyconnectedcomponents"], Return[$Failed], {ReactionsData::wrreac,ReactionsData::badarg}];
			If[rdata === $Failed, Return[$Failed]];

			{n, r, fhj, cxs, fhjscc} = rdata;

			graphfunc = (PlotFunction /. cc) /. PlotFunction->"GraphPlot";			
			If[GraphPlotFunctionQ[graphfunc],
				graphfunc = ToExpression[graphfunc],
				Message[ShowFHJGraph::"funcarg",graphfunc];
				Return[$Failed]				
			];

			fhjscc = Last /@ fhjscc; 
			If[((Numbered /. cc) /. Numbered -> False)===True, 
				cxs = Thread[Rule[cxs,Range[n]]];
				fhj = fhj /. cxs;
				fhjscc = fhjscc /. cxs;
				cxs = Last /@ cxs;
			];

			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 = fhj /. Flatten[Thread[First[#]->Thread[Style[First[#], Plain, Last[#]]]] &/@ Thread[{fhjscc,sc}]];
						If[rates==={},
							graphfunc[rscc,nopt],
							If[r===Length[rates],
								graphfunc[Thread[{rscc,rates}],nopt],
								Message[ShowFHJGraph::"srates"];
								Return[$Failed]
							]
						],
					Message[ShowFHJGraph::"sccomps"];
					$Failed
				],

				If[c =!= {} && c =!= ComplexColours,

					If[VectorQ[c] && n===Length[c],
						rs = fhj /. (First[#]->Style[First[#], Plain, Last[#]] &/@ Thread[{cxs,c}]);
						If[rates==={},
							graphfunc[rs,nopt],
							If[r===Length[rates],
								graphfunc[Thread[{rs,rates}],nopt],
								Message[ShowFHJGraph::"srates"];
								Return[$Failed]
							]
						],
						Message[ShowFHJGraph::"ccols"];
						$Failed
					],

					If[rates === {},
						GraphPlot[fhj,nopt],
						If[r===Length[rates],
								graphfunc[Thread[{fhj,rates}],nopt],
								Message[ShowFHJGraph::"srates"];
								Return[$Failed]
							]
					]
				]
			]
	];

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


(* ::Subsubsection::Closed:: *)
(*Volpertindexing, ShowVolpertGraph, VolpertSpecies*)


(*
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[{reactions}, initspecies, options] gives the Volpert indices \
of the species and reactions with respect to the initspecies given as \
the second argument. VolpertIndexing[] gives a pretty printed table of the results, \
with the reactions denoted by the strings given in the list tags.";


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::verbtrumis = "The given species and/or reactionsteps, whose indices are required, are not contained in \
the current reaction mechanism or the option Verbose is prevent to be printed. \
Let the option Verbose be True in order to get more info.";

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 -> {}, Verbose -> False };

VolpertIndexing[{reactions__}, initspecies_?VectorQ, specsreacs_List : {}, OptionsPattern[]]:=
	Module[{
			rdata, m, r, sp, reacs, vggraph, vgraph, w, w0, indexed, complx, initcomplx,
			indexedr, indexedsp, a, b, v, req, z1, z2, t1, t2, vindeces
		   },

		rdata = Check[ReactionsData[{reactions},OptionValue[ExternalSpecies]]["M","R","species","fhjgraph","volpertgraph"],Return[$Failed],
							     {ReactionsData::wrreac,ReactionsData::badarg}];

		If[rdata === $Failed, Return[$Failed], {m, r, sp, reacs, vggraph} = rdata; ];
 
		complx = Transpose[reacs /. Rule->List];
		initcomplx = First[complx];

		If[MemberQ[Flatten[complx],"0"], sp = Append[sp,"0"]; m = m+1; ];
		If[MemberQ[initcomplx,"0"], w0 = {m};, w0 = {};];

		w = Flatten[Position[sp, #] &/@ DeleteDuplicates[initspecies]];
		If[w =!= {} || initspecies === {},
	
			reacs = reacs /. Rule->RightArrow;
			indexedr = Last[vggraph];
			indexedsp = MapIndexed[#1->First[#2]&, sp];
			indexed = Join[indexedr,indexedsp];
			vgraph = First[vggraph];

			a = Cases[vgraph, {Rule[x_,y_RightArrow],z_}:> Rule[{x,y} /. indexed, z]];
			b = Cases[vgraph, {Rule[x_RightArrow,y_],z_}:> Rule[{y,x} /. indexed, z]];

			v = VolpertIndexing[SparseArray[a, {m, r}], SparseArray[b, {m, r}], DeleteDuplicates[Join[w,w0]]];
			vindeces = MapThread[Thread[Rule[#1,#2]]&, {{sp, reacs}, v}];,

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

		];

		If[OptionValue[Verbose],

			req = Flatten[specsreacs]; (*ToCanonicalForm*)		
			
			If[req === {}, 
						z1 = sp; 
						z2 = reacs;,
						z1 = Intersection[sp, req];
						z2 = Intersection[reacs, req];
						If[ z1 === {}, 
								If[ z2 === {},
										Message[VolpertIndexing::"verbtrumis"]; 
										Return[$Failed], 
										z1 = sp;
								];,
								If[ z2 === {}, z2 = reacs; ];
						];
			];
			z2 = z2 /. Rule->RightArrow;
			t1 = Sort[Transpose[{z1, z1 /. First[vindeces]}], #1[[2]]<#2[[2]]&];
			t2 = Sort[Transpose[{z2, z2 /. Last[vindeces]}], #1[[2]]<#2[[2]]&];

			{
				Grid[Join[{{"Species","Indices"}},t1],Dividers->{{False,True},{False,True}}],
				Grid[Join[{{"Reaction steps","Indices"}},t2],Dividers->{{False,True},{False,True}}]
			}
			(*
			mz = Length[z1];
			rz = Length[z2];
			z1 = Join[#, Array[""&, Max[rz-mz, 0]]] & /@ {z1, z1 /. First[vindeces]};
			z2 = Join[#, Array[""&, Max[mz-rz, 0]]] & /@ {z2, z2 /. Last[vindeces]};
			Grid[Join[{{"Species", "Indices", "Reactionsteps", "Indices"}}, Transpose[Join[z1, z2]]],
							Dividers -> {{False,False,True}, {False,True}}]
			*)
			,
			vindeces
		]
	];

(*
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[___] := (Message[VolpertIndexing::"badarg"]; $Failed)


GraphPlotFunction2Q := MemberQ[{"Graph","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", EdgeLabeled -> False, 
								Highlight -> {}, SubgraphHighlight -> False, Bipartitized -> False, Indexed -> False, Numbered -> False};

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

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

ShowVolpertGraph[{reactions__},opts___?OptionQ] := 
	Module[{
			fullopts, graphfunc, exps, vgp, vindices, vtempind, 
			vg, vind, complx, numbd, ind, initsp, highlight, edgelabeled, 
			nopts, g, bipd, fhj, sp, m, r
		   },

		fullopts = Join[{opts},Options[ShowVolpertGraph]];		
		graphfunc = PlotFunction /. FilterRules[fullopts, PlotFunction];

		If[GraphPlotFunction2Q[graphfunc],

			exps = FilterRules[{opts}, ExternalSpecies];
			vgp = Check[ReactionsData[{reactions},exps]["volpertgraph"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];
			
			If[ vgp === $Failed, Return[$Failed];, {vg, vind} = vgp;];

			highlight = Flatten[{Highlight /. FilterRules[fullopts, Highlight]}]; (*ToCanonicalForm*)

			ind = Indexed /. FilterRules[fullopts, Indexed];
			sp = ReactionsData[{reactions},exps]["species"];
			initsp = Intersection[Flatten[{ind}],sp];
			If[ ind =!= False, 
				vindices = VolpertIndexing[{reactions},initsp] /. Rule[x_,y_] :> Rule[x,{x,y}];
				(*Print[{sp,initsp,Flatten[{ind}]}];*)
				If[ind =!= True && initsp === {} && Flatten[{ind}] =!= {}, Message[VolpertIndexing::"specerr",Flatten[{ind}]];];
				vg = Cases[vg, {Rule[x_,y_],z_} :> If[Head[x]===RightArrow,
														{Rule[x /. Last[vindices], y /. First[vindices]],z},
														{Rule[x /. First[vindices], y /. Last[vindices]],z}] 
					 ];
				(*
				vg = vg /. x_RightArrow :> {x /. vtempind, Last[x /. Last[vindices]]};
				vg = vg /. First[vindices] /. (Reverse /@ vtempind);
				*)
				vtempind = vind /. Rule[x_,y_] :> Rule[x,Subscript[xyz,y]];
				highlight = highlight /. x_RightArrow :> {x /. vtempind, Last[x /. Last[vindices]]};
				highlight = highlight /. First[vindices] /. (Reverse /@ vtempind);
				(*If[Head[#]===RightArrow, {#, # /. Last[vindices]}, {#, # /. First[vindices]}] &/@ highlight;*)
			];
			(*VertexReplace[, Table[i -> Subscript[v, i], {i, 4}]]*)

			numbd = (Numbered /. FilterRules[fullopts, Numbered]) === True; 
			If[ numbd, 
					vg = vg /. vind; 
					highlight = highlight /. vind;
			];

			edgelabeled = (EdgeLabeled /. FilterRules[fullopts, EdgeLabeled]) === True;

			If[ highlight =!= {},
				
				If[ edgelabeled, vg = Labeled@@# &/@ vg;, vg = First /@ vg; ];
				(*ez itt erzeketlen a plotfunction-ra, illetve a bipartitized-ra*)
				g = Graph[vg];
				If[ (SubgraphHighlight /. FilterRules[fullopts,SubgraphHighlight]) === True, highlight = Subgraph[g,highlight];]; 
				HighlightGraph[g, highlight, FilterRules[{opts},Options[HighlightGraph]]],
				
				nopts = FilterRules[{opts},Options[ToExpression[graphfunc]]];		
				Switch[graphfunc,
						"Graph",
							If[ edgelabeled, vg = Labeled@@# &/@ vg;, vg = First /@ vg; ]; 
							(*ez meg itt erzeketlen a bipartitized-ra*)
							Graph[vg, nopts],
						_, 
							bipd = (Bipartitized /. FilterRules[fullopts, Bipartitized]) === True;
							nopts = Join[(FilterRules[nopts, EdgeLabeling] /. {}->{EdgeLabeling->False}), nopts];
							If[ edgelabeled, nopts = Join[{EdgeLabeling -> True}, nopts]; ];
							If[ bipd,
								fhj = ReactionsData[{reactions},exps]["fhjgraph"] /. Rule->RightArrow;
								{m, r} = ReactionsData[{reactions},exps]["M","R"];
								complx = ReactionsData[{reactions},exps]["complexes"];
								If[MemberQ[complx,"0"], sp = Append[sp, "0"]; m = m+1;];
								If[ ind =!= False, 
									sp = sp /. Flatten[vindices];
									fhj = fhj /. Flatten[vindices];
								];
								If[ numbd, fhj = fhj /. vind;];
								(*Print[{vg,sp,fhj,m,r,nopts}];*)

								GraphPlot[vg,
									EdgeRenderingFunction -> (
										If[
											MemberQ[{RightArrow,Integer},Head[If[VectorQ[#2[[1]]],#2[[1,1]],#2[[1]]]]], 
											{Orange, Arrow[#1,0.15], If[edgelabeled, Inset[Cases[vg,{Rule@@#2,x_}:>x][[1]], Mean[#1], Background->White], Sequence[]]}, 
											{Blue, Arrow[#1,0.15], If[edgelabeled, Inset[Cases[vg,{Rule@@#2,x_}:>x][[1]], Mean[#1], Background->White], Sequence[]]}
										]&),
									VertexCoordinateRules -> 
										Join[
											Table[(sp[[i]])->{0,i},{i,1,m,1}],Table[(fhj[[i]])->{Ceiling[Min[{m,r}]/2],i},{i,1,r,1}]
										], nopts
								],
								ToExpression[graphfunc][vg, nopts]
							]
					]
			  ],
			Message[ShowVolpertGraph::"funcarg",graphfunc];
			Return[$Failed]
		]
	];

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


(*SubsetQ[x_,y_] := (Intersection[x,y] == x); ez miert nem jo?*)
SubsetQ[x_,y_] := And @@ (MemberQ[x,#] &/@ y);


VolpertSpecies::usage = "VolpertSpecies[{reactions},option] finds the set of possible species \
that obtains finite Volpert indexing.";


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

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

VolpertSpecies[{reactions__}, OptionsPattern[]] :=
	Module[{
			fhj, fhjsp, initcompl, graph, strongcomp, fhjstronglength, 
			fhjstronggraph, volpertbasis, vertxs
		   },

		fhj = Check[ReactionsData[{reactions},OptionValue[ExternalSpecies]]["fhjgraph"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];
		If[fhj === $Failed, Return[$Failed];];

		fhjsp = DeleteDuplicates[(complextospecies /@ #) &/@ fhj] /. Rule->RightArrow;
		initcompl = Union[First /@ fhjsp];
		(*Print[fhjsp];
		Print[initcompl];*)

		graph = DeleteDuplicates[
					Flatten[
						Outer[
							{If[SubsetQ[#2,First[#1]],DirectedEdge[#2, #1],{}],If[SubsetQ[Flatten[#1 /. RightArrow->List],#2],DirectedEdge[#1, #2],{}]}&,
							fhjsp, initcompl,1
						]
					]
				];
		(*Print[Graph[graph,VertexLabels->"Name"]];*)
		strongcomp = ConnectedComponents[Graph[graph]];
		(*Print[strongcomp];*)

		fhjstronglength = Length[strongcomp];
		vertxs= Range[fhjstronglength];
		fhjstronggraph = Graph[
							DeleteDuplicates[
								graph /. 
									Sort[Flatten[MapThread[Thread[ #1 -> #2 ]&, {strongcomp, vertxs}]], 
										Length[First[#1]]>Length[First[#2]]&
									]
							] /. (DirectedEdge[y_,y_] :> DirectedEdge[y,"\[Infinity]"])];

		volpertbasis = Map[If[
							VertexInDegree[fhjstronggraph,#]===0, 
							strongcomp[[#]] /. x_RightArrow :> Sequence[],
							{{}}
							]&, vertxs
		   			] /. {{}} :> Sequence[];
		
		If[ Length[volpertbasis] === 1, 
				Sort[First[volpertbasis]], 
				Sort[volpertbasis] 
		]
		
	];

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


(* ::Subsubsection::Closed:: *)
(**MassConservationQ, Atoms, ToAtomMatrix, AtomConservingQ, FromAtomMatrix*)


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)


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.";

Options[ToAtomMatrix] := {FormattedOutput->False};

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

ToAtomMatrix[molecules_?VectorQ,opts___?OptionQ]:=
	Module[{ ats, l, toa, ind, atomm, fout, grid },

		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};

		fout = (FormattedOutput /. Flatten[{opts}]) /. FormattedOutput -> False;

		If[fout,

			grid = FilterRules[{opts},Options[Grid]];
			Grid[MapThread[Prepend,{Prepend@@toa,Prepend[molecules,"Molecules\[Backslash]Atoms"]},1],grid],

			ind = Sort[MapIndexed[{#1,First[#2]}&,Most[Last[toa]]]];
			atomm = Transpose[First[toa]];

			Join[ {Join[First/@ind,{"charge"}]}, {Transpose[Join[atomm[[#]]& /@ (Last/@ind), {-Last[atomm]}]]} ]
		]
	];

(* Prepend[MapThread[Prepend,{atommatrix, molecules},1],Join[{"molecules \\ atoms"},atoms]] *)

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 = ToAtomMatrix[ls];
		(*Transpose[a].atommatrix == Transpose[b].atommatrix*)

		(#===ConstantArray[0,Length[First[atommatrix]]]) &/@ (Transpose[alpha-beta].Last[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::Closed:: *)
(*DetailedBalanced*)


USet[a_Integer,b_Integer,s_List] :=
	Module[{sa = FSet[a,s], sb = FSet[b,s], set = s},
	      If[ set[[sa,2]] < set[[sb,2]], {sa,sb} = {sb,sa} ];
	        set[[sa]] = {sa, Max[ set[[sa,2]], set[[sb,2]]+1 ]};
	        set[[sb]] = {sa, set[[sb,2]]};
	        set
	];

FSet[n_Integer,s_List] := 
        Block[{$RecursionLimit = Infinity}, 
              If [n == s[[n,1]], n, FSet[s[[n,1]],s]]
        ];

MST[e_List,n_Integer]:=
	Module[{sorte, i, s},
			sorte = Sort[e, (#1[[2]] <= #2[[2]])&];
			s = Table[{i,1},{i,n}];
			Select[sorte,(If[FSet[#[[1]],s]!=FSet[#[[2]], s], s = USet[#[[1]],#[[2]], s]; True,False])&]
	];


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}];

			(*If[rr === $Failed, Return[$Failed]];*)

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

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

		]
	];

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

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

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

			{cxes, m, r, nn, ll, fhj, gamma, deficiency} = 
					Check[ReactionsData[{reactions},exs]["complexes","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[StringReplace[deficiency,Longest[__]~~"="->""]];
							(*ToExpression[Last[StringSplit[deficiency,"="]]];*)

			fhj = fhj /. Sort[Thread[cxes->Range[nn]], Length[#1[[1]]/.Plus->List]>Length[#2[[1]]/.Plus->List]&];

			fhjunlist = Union[Sort /@ (fhj /. Rule->UndirectedEdge)];
			fhjun = Graph[fhjunlist];
			rsteplistind = fhj /. Rule->List;
			shadow = Array[Subscript[shw,#]&,r];
			rrates = Thread[rsteplistind -> shadow];
			(**)
			fhjun0 = fhjunlist /. UndirectedEdge->List;
			mintreeun = MST[fhjun0,nn];
			mintree = Graph[UndirectedEdge@@# &/@ mintreeun];
			(**)
			rev = Thread[shadow -> rates];
			(**)
			circconds := (
							If[Length[c0 = Complement[fhjun0, mintreeun]] =!= r/2-nn+ll,
								Message[DetailedBalanced::"misscirc"];
							];
							circs = Map[Partition[Append[FindShortestPath[mintree,#[[1]],#[[2]]],#[[1]]],2,1]&, c0];

							Map[
								(Times @@ (# /. rrates) == Times @@ (Map[Reverse,#] /. rrates)) /. rev
								&, circs
							]
						 );
			(**)
			spanforconds := (
							spanreacs = Flatten[Position[rsteplistind, #]& /@ mintreeun];
							If[Length[mm = NullSpace[Part[gamma, All, spanreacs]]] =!= deficiency,
								Message[DetailedBalanced::"missspf"];
							];

							Map[
								Numerator[Together[
											Times @@ ((mintreeun /. rrates)^#) + Times @@ ((Map[Reverse,mintreeun] /. rrates)^#)
										 ]
								] /. Plus -> Equal /. rev				
								&, mm
							]
							);

			If[ deficiency === 0,

				If[AcyclicGraphQ[fhjun],

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

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

				If[AcyclicGraphQ[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[{reactions__}, rates_?VectorQ, variables_?VectorQ, OptionsPattern[]] := 
	RightHandSide[Sequence@@Check[ReactionsData[{reactions},OptionValue[ExternalSpecies]]["\[Alpha]","\[Beta]"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}],variables,rates];

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

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

RightHandSide[alpha_?MatrixQ,beta_?MatrixQ,variables_?VectorQ,ratefunc_Function] := 
	Dot[beta-alpha, Check[ratefunc @@ variables, Return[$Failed];, {Function::"fpct",Function::"slotn",General::"stop"}]];

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 number of rate constants (or triplets) is unequal to the number of given reactionsteps.";
DeterministicModel::enthalpies = "The given number of enthalpies is 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.";
DeterministicModel::steady = "The given vector for steady state concentrations has wrong size.";

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];

(*without Arrhenius parameters*)
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
			]
	];

(*with Arrhenius parameters, KELL: vars-os + homerseklet mas jelolessel,V mas jelolessel*)
DeterministicModel[{reactions__},arrhenius_?MatrixQ,enthalpies_?VectorQ,C0_,V_,Ta_,tres_,p1_,p2_,s_Symbol:Global`t,OptionsPattern[]] :=
	Module[{ kkk, maf, alpha, beta, vars, reacs = Check[ReactionsData[{reactions},OptionValue[ExternalSpecies]]["R"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}]},
			If[reacs===Length[arrhenius],
				If[reacs===Length[enthalpies],
					{alpha, beta, vars} = ReactionsData[{reactions},OptionValue[ExternalSpecies]]["\[Alpha]","\[Beta]","variables"];
					vars = Through[vars[s]];
					If[Length[C0]===Length[vars],
			(*kiemelni*)	kkk[{a_,n_,e_}] := a Global`T[s]^n Exp[- e /( Global`T[s] MolarGasConstant Kelvin Mole/Joule )];
					maf = MassActionKineticFunction[alpha,vars,kkk/@arrhenius];
					{Join[Thread[D[vars,s] == (Dot[beta-alpha,maf] + (C0-vars)/tres)], {D[Global`T[s],s] == - Dot[maf,enthalpies]/p1 - (1/tres + p2/p1)(Global`T[s]-Ta)/V}], Join[vars,{Global`T[s]}]},
					Message[DeterministicModel::"steady"];
					$Failed
					],
					Message[DeterministicModel::"enthalpies"];
					$Failed
				],
				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 *)
(*generally*)
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
			]
	];

(*with Arrhenius parameters*)
Concentrations[{reactions__},arrhenius_?MatrixQ,enthalpies_?VectorQ,C0_?VectorQ,V_?Positive,Ta_?Positive,tres_?Positive,p1_?Positive,p2_?Positive,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},arrhenius,enthalpies,C0,V,Ta,tres,p1,p2,s,exs],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg,DeterministicModel::"rates",DeterministicModel::"enthalpies",DeterministicModel::"steady"}];
			If[Length[init]===(ReactionsData[{reactions},exs]["M"]+1),
				{ 
					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] gives us the complex balanced equilibriums.";


Options[StationaryPoints] := {ExternalSpecies->{}, Positivity->False, Conditions->{}, Method->"Automatic"}; (*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::unrecogmethod = "The system did not recognize the given `1` method, try Automatic of LinearFirstIntegrals.";

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,Superscript[#,"*"]& /@ 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,Superscript[#,"*"]& /@ vars,opts]
	];

StationaryPoints[{reactions__},rates_?VectorQ,init_?VectorQ,opts___?OptionQ] := 
	StationaryPoints[{reactions},rates,init,Check[Superscript[#,"*"]& /@ 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, gamma, m, r, method, y, z, statp, balances, rateconst, inimass, nonneg, conds, sol, fless, null, linfirstin },
			
			ropts = Sequence@@FilterRules[{opts},Options[Reduce]];
			
			{alpha, beta, gamma, m, r} = 
					Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["\[Alpha]","\[Beta]","\[Gamma]","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],
			
				Switch[(method = (Method/.{opts}) /. Method->"Automatic"),

					"Automatic",

					z = Subscript[y,ToString[#]]& /@ Range[r];

					{statp,balances,rateconst,inimass,nonneg} = 
						Thread/@{RightHandSide[alpha,beta,vars,rates] == 0, (beta-alpha).z == vars-init, rates > 0, Cases[init,Except[0]]>0,fless[0,vars]};

					conds = Join[balances,rateconst,inimass,nonneg];

					If[(sol = Reduce[Join[statp,conds,(Conditions/.{opts})/.Conditions->{}],Join[vars,z],Reals,Backsubstitution->True,ropts])===False,
						Message[StationaryPoints::noeql];
						Return[$Failed],

						Join[{vars},{Sort[Sort[MyToRules[Select[#/.And->List,Not[MemberQ[z,First[#]]]&]]]&/@Flatten[{sol/.Thread[rateconst->True]/.Thread[inimass->True]/.Or->List}]]}]
					],

					"LinearFirstIntegrals",

					null = NullSpace[Transpose[gamma]] /. {}->{ConstantArray[0,m]};

					{statp,linfirstin,rateconst,inimass,nonneg} = 
						Thread/@{RightHandSide[alpha,beta,vars,rates] == 0, null.(vars-init) == 0, rates > 0, Cases[init,Except[0]]>0,fless[0,vars]};

					conds = Join[rateconst, inimass, nonneg];

					sol := Reduce[Join[linfirstin, statp, rateconst, inimass, nonneg, (Conditions/.{opts})/.Conditions->{}],vars, Reals, Backsubstitution->True, ropts];

					If[sol === False,
						Message[StationaryPoints::noeql];
						Return[$Failed],

						Join[{vars},{Sort[Sort[MyToRules[# /. And->List]]& /@ Flatten[{sol/.Thread[conds->True]/.Or->List}]]}]
					],

					_,

					Message[StationaryPoints::unrecogmethod,method];
					Return[$Failed]
				],

				Message[StationaryPoints::arglth,rates,init,vars];
				Return[$Failed]
			]
	];

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


MassConservation::usage = "MassConservation[{reactions__},init,vars,opts] computes the linear first integrals of formal mechanism {reactions}, \ 
and gives us also the solution of its solution(s).";


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, gamma, m, ns, ropts, conds, sol },

			{alpha, beta, gamma, m} = Check[ReactionsData[{reactions},FilterRules[{opts},ExternalSpecies]]["\[Alpha]","\[Beta]","\[Gamma]","M"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];

			If[m===Length[init] && m===Length[vars],

				ns = NullSpace[Transpose[gamma]];

				ropts = Sequence@@FilterRules[{opts},Options[Reduce]];

				conds = Flatten[Thread /@ {Cases[init,Except[0]]>0, vars >= 0}];

				sol := Reduce[Join[Thread[ns.(vars - init) == 0],conds],vars,Reals,Backsubstitution->True,ropts];

				If[ns === {} (*|| sol === False*), 

					Message[MassConservation::nomass];
					$Failed,

					(*Join[{Thread[ns.vars == ns.init]},{vars},{Sort[Sort[MyToRules[#/.And->List]]&/@Flatten[{sol/.Thread[conds->True]/.Or->List}]]}]*)
					Thread[ns.vars == ns.init]

				],

				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;

			Eigensystem[D[RightHandSide[alpha,beta,newvars,rates],{newvars}]],

			Message[EigenvaluesJacobianMatrix::arglth,rates,vars];
			$Failed
		]
	];

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


(* ::Subsubsection::Closed:: *)
(*AbsoluteConcentrationsRobustness*)


AbsoluteRobustness::usage = "AbsoluteRobustness[] checks the Shinar-Feinberg conditions...";


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[{reactions}],

			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[{reactions}],

				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[{reactions}],
						False,
							Message[AbsoluteRobustness::"maybe","There are no positive stationary points"]; 
							Return[{reactions}],
						True,
							Union[acr /. r_*A_?StringQ :> A /. A_?StringQ*r_ :> A /. r_*A_?SpeciesQ :> A /. A_?SpeciesQ*r_ :> {A}](*!!!!*)
				]	
			]
		]
	];

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


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


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

INTENSITY = Compile[{{x,_Integer,1},{alphatr,_Integer,2}},

				Times @@ FactorialPower[x,#] &/@ alphatr,

				(*RuntimeAttributes->{Listable},*) 
				Parallelization->Automatic,
				RuntimeOptions->{"WarningMessages"->False,"CatchMachineIntegerOverflow"->False} (*"EvaluateSymbolically"*)

			];

(*
MassActionKineticIntensity[alphatr_,variables_,rates_] := rates * ((Times @@ Power[variables, #]) &/@ alphatr);
*)

SetAttributes[HornerForm,Listable];

(*MyDerivative[d_, f_, x_] := Derivative[Sequence @@ d][f][Sequence @@ x];*)

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

avogadroconstant = 6.02214179`*^23;

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


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


ChemicalMasterEquation::usage = "...";


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

ChemicalMasterEquation::badarg = "Illegal argument of function ChemicalMasterEquation.";
ChemicalMasterEquation::rvarg = "Argument `1` or `2` has wrong size. \
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, alphatr, gammatrind, n},
			
		{alpha, gamma, m, r} = Check[ReactionsData[{reactions},OptionValue[ExternalSpecies]]["\[Alpha]","\[Gamma]","M","R"],Return[$Failed];,{ReactionsData::wrreac,ReactionsData::badarg}];
			
			alphatr = Transpose[alpha];

			gammatrind = MapIndexed[{#1, First[#2]}&, Transpose[gamma]];

			If[ Length[rates] === r && Length[vars] === m+1,

				{
					D[h @@ vars, First[vars]] == 

					- (h @@ vars)*Total[FunctionExpand[Intensity[rates, Rest[vars], alphatr]]]

					+ Total[Map[
						h@@(vars-Prepend[First[#],0])*
							FunctionExpand[rates[[Last[#]]]*(Times@@FactorialPower[Rest[vars]-First[#],alphatr[[Last[#]]]])] &, gammatrind
							]
					  ]
				},
				
				Message[ChemicalMasterEquation::"rvarg", rates, vars];
				Return[$Failed];

			  ]

	];

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, 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[#]==0/.(h[Global`t,q___]:>h[q])) &/@ ChemicalMasterEquation[{reactions}, rates, Prepend[vars,Global`t], h, ExternalSpecies->exs],

				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:: *)
(**MomentEquations*)


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]*)


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


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

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

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

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

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

(* arg. csere*)
MomentEquations[{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]]]) =!= {},
						MomentEquations[{reactions}, #, rates, e, var, ExternalSpecies->exs, CombinatorialMoments->combmom]& /@ ((UnitVector[m,#]& /@ sppos)*Last[a1]),
						Message[MomentEquations::"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[MomentEquations[{reactions},#,rates,e,var,ExternalSpecies->exs,CombinatorialMoments->True]]& /@ Rest[Tuples[Prepend[Range[Max[a]],0],l]]]
						],
						Message[MomentEquations::"aarg",a];
						Return[$Failed];
					],
					_,
						Message[MomentEquations::"nonnegarg",a];
						Return[$Failed];
			]
	];

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

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


(* 
SolveMomentEquations

*)


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


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[rates,X0,alphatr]; (*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[Part[rates,newrs],Last[next],newrs/.alphatrules]]];
				(*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[rates,vars,alphatr];
		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[rates,vars,alphatr];
		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[Thread[vars-Last[#]==Total[rhs*gammatr]],vars,Reals]];*)
				sol = Flatten[FindInstance[Join[nonneg,Thread[vars-Last[#]==Total[rhs*gammatr]]],vars,Reals]];
				(*Print[sol];*)
				(*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::intarg = "The volume is finite, so the initial numbers of species are considered to be nonnegative integers.";
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, coi, 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[v===Infinity && (Or @@ (Not[IntegerQ[#]]& /@ X0)), Message[Simulation::"intarg"]; Return[$Failed];];
		If[Positive[v],(*NumberQ[v]*)
			cou[vol_,atr_] := If[vol===Infinity, rates, ChangeOfUnits[vol,rates,atr]];
			coi[vol_,atr_] := If[vol===Infinity, atr, Round[atr*avogadroconstant*vol]];,
			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],coi[v,X0],maxtime,maxit,verbose],
					"FirstReactionStep",
						FirstReactionStepMethod[sp,alphatr,gammatr,cou[v,alphatr],coi[v,X0],maxtime,maxit,verbose],
					"ModifiedFirstReactionStep", 
						Return[$Failed];,
					"NextReactionStep", 
						NextReactionStepMethod[sp,m,r,alphatr,gammatr,cou[v,alphatr],coi[v,X0],maxtime,maxit,verbose],
					"ExplicitTau-Leaping",
						ExplicitTauLeapingMethod[sp,alphatr,gammatr,cou[v,alphatr],coi[v,X0],maxtime,maxit,fss,verbose],
					"ImplicitTau-Leaping",
						ImplicitTauLeapingMethod[sp,alphatr,gammatr,cou[v,alphatr],vars,coi[v,X0],maxtime,maxit,fss,verbose],
					"TrapezoidalTau-Leaping",
						TrapezoidalTauLeapingMethod[sp,alphatr,gammatr,cou[v,alphatr],vars,coi[v,X0],maxtime,maxit,fss,verbose],
					"K-Leaping",
						error = fss/2;
						KLeapingMethod[sp,alphatr,gammatr,cou[v,alphatr],coi[v,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::intarg,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::intarg,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::intarg,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:: *)
(*CompareMethods*)


(* ::Subsubsection::Closed:: *)
(*CompareModels*)


CompareModels::usage = "...";


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

CompareModels[{reactions__}, rates_?VectorQ, c0_?VectorQ, maxtime_?NumericQ, volume_?NumericQ, opts___?OptionQ]:=
	Module[{
			simopt = FilterRules[{opts},Join[Options[SimulationPlot],Options[Simulation]]],
			plotopt = FilterRules[{opts},Options[Plot]],
			concopt = FilterRules[{opts},Options[Concentrations]], plotopts
		   },
		Show[
			SimulationPlot[{reactions},rates,c0,maxtime,Join[{Volume->volume},simopt,plotopt]],
			Plot[Evaluate[(ReplaceAll@@Concentrations[{reactions},rates,c0,{0,maxtime},concopt])*avogadroconstant*volume],
			{Global`t,0,maxtime},Evaluate[Complement[plotopt,simopt,concopt]]]
		]
	];

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


(* ::Subsubsection::Closed:: *)
(**new methods, Langevin equation, cont. approx.*)


(*
based on int. approx.
Langevin

*)


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


(* ::Subsubsection::Closed:: *)
(*aux functions, 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).. }] &;


(* ::Subsubsection::Closed:: *)
(*decomposition functions*)


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, MyFilterOptions[Decompositions, opts]}]],
						Insert[#, -2, {{i}}] & /@ Decompositions[2 molec, others, maxlen, Flatten[{Verbose -> False, MyFilterOptions[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, MyFilterOptions[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)
(*  ElementaryStep, ez mi? nem ElementaryReactions?  *)


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 = MyFilterOptions[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)


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, nem constraints? *)
		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, MyFilterOptions[LPPartitions, opts]],
				ContejeanDevie, CDPartitions[overall, steps, maxlen, MyFilterOptions[CDPartitions, opts]]
				(* , SubsetEnumeration, SEPartitions[overall, steps, maxlen, MyFilterOptions[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 -> 24]}, (*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 -> "ReactionKinetics 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["\"fhjweaklyconnectedcomponents\""]},
			{ButtonBox["\"fhjstronglyconnectedcomponents\""]},
			{ButtonBox["\"fhjterminalstronglyconnectedcomponents\""]},
			{ButtonBox["\"reactionorders\""]},
			{ButtonBox["\"reactionsteps\""]},
			{ButtonBox["\"species\""]},
			{ButtonBox["\"variables\""]},
			{ButtonBox["\"volpertgraph\""]}
		(*  {ButtonBox["Models"]},
			{ButtonBox["Reactions"]}  *)
		},RowSpacings ->.9, ColumnSpacings ->1.]],
        FontFamily -> "TimesNewRoman", FontWeight -> "Normal", FontSize -> 22]}, (*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 -> "ReactionKinetics 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["dpapp@iems.northwestern.edu","mailto:dpapp@iems.northwestern.edu"], Style[", or "],
			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["ReactionKinetics Palette"]] === 0,
			OpenReactionKineticsPalette[];
	];
	If[ Length[Notebooks["ReactionKinetics 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`T,
(*	Global`x,*)
	Global`X,(**)
	Global`z,
	Bipartitized,
	CombinatorialMoments,
	ComplexColours,
	Conditions,
	EdgeLabeled,
	ExternalSpecies,
	FixedStepSize,
	FormattedOutput,
	GeneratedRateCoefficient,
	Highlight,
	Indexed, 
	MaxIteration,
	Memo,
	Method,
	Numbered,
	PlotFunction,
	Positivity,
	Side,
	Species,
	StronglyConnectedComponentsColours,
	SubgraphHighlight,
	TimeLimit,
	Verbose,
	Volume
];


SetAttributes[
	{
		ReactionImport,
		(**)
		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,
		MomentEquations,
		DependencyGraph,
		ShowDependencyGraph,
		Simulation,
		SimulationPlot,
		SimulationPlot2D,
		SimulationPlot3D,
		CompareModels,
		(**)
		OpenReactionKineticsPalette,
		OpenReactionKineticsNamesPalette,
		ReactionRatesNotebook,

		(**)
		CHEMKINExport,
		CHEMKINImport, 
		FilterReactions,
		FromStoichiometry, 
		DeleteAutocatalysis,
		ReversibleFHJRepresentation,
		MinFHJWeaklyConnectedComponents, 
		MinFHJStronglyConnectedComponents, 
		MaxFHJWeaklyConnectedComponents, 
		MaxFHJStronglyConnectedComponents,
		VolpertSpecies

	},
	{Protected,ReadProtected}
];


SetAttributes[
	{
	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,
(**)
	VolpertIndexingWithMatrices,
	ComponentwiseSmallerEqualQ, 
	MinimalQ, 
	MyMinimize,
	MyMaximize,
	ZeroVectorQ, 
	MyFloor, 
	MyCeiling, 
	SMaux, 
	CDPartitions, 
	UnrestrictedSolutions, 
	FindAllSolutionsWithPreprocessing, 
	LPPartitions, 
	FindAllSolutions
    },
	{Protected,ReadProtected}
];


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


EndPackage[ ];
