text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2016 Varnerlab
# School of Chemical Engineering Purdue University
# W. Lafayette IN 46907 USA

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #
function DataFile(TSTART,TSTOP,Ts)
# ----------------------------------------------------------------------------------- #
# DataFile.jl was generated using the Kwatee code generation system.
# DataFile: Stores model parameters as key - value pairs in a Julia Dict()
# Username: jeffreyvarner
# Type: GRN-JULIA
# Version: 1.0
# Generation timestamp: 08-08-2016 20:01:41
#
# Input arguments:
# TSTART  - Time start
# TSTOP  - Time stop
# Ts - Time step
#
# Return arguments:
# data_dictionary  - Data dictionary instance (holds model parameters)
# ----------------------------------------------------------------------------------- #

# Load the stoichiometric matrix -
S = float(open(readdlm,"../Network.dat"));
(NSPECIES,NREACTIONS) = size(S);

# How many genes do we have in the model? -
number_of_genes = 23;

# Formulate the initial condition array -
initial_condition_array = Float64[];
push!(initial_condition_array,1.0);	#	1	gene_Smad
push!(initial_condition_array,1.0);	#	2	gene_Ecadherin
push!(initial_condition_array,1.0);	#	3	gene_Bcatenin
push!(initial_condition_array,1.0);	#	4	gene_PLCy
push!(initial_condition_array,1.0);	#	5	gene_PKC
push!(initial_condition_array,1.0);	#	6	gene_RAF
push!(initial_condition_array,1.0);	#	7	gene_RAS
push!(initial_condition_array,1.0);	#	8	gene_AP1_SP1
push!(initial_condition_array,1.0);	#	9	gene_GSK3
push!(initial_condition_array,1.0);	#	10	gene_tubulin
push!(initial_condition_array,1.0);	#	11	gene_SNAIL_SLUG
push!(initial_condition_array,1.0);	#	12	gene_Calcineurin
push!(initial_condition_array,1.0);	#	13	gene_APC_AXIN
push!(initial_condition_array,1.0);	#	14	gene_PI3K
push!(initial_condition_array,0.1);	#	15	mRNA_Smad
push!(initial_condition_array,0.1);	#	16	mRNA_Ecadherin
push!(initial_condition_array,0.1);	#	17	mRNA_Bcatenin
push!(initial_condition_array,0.1);	#	18	mRNA_PLCy
push!(initial_condition_array,0.1);	#	19	mRNA_PKC
push!(initial_condition_array,0.1);	#	20	mRNA_RAF
push!(initial_condition_array,0.1);	#	21	mRNA_RAS
push!(initial_condition_array,0.1);	#	22	mRNA_AP1_SP1
push!(initial_condition_array,0.1);	#	23	mRNA_Calcineurin
push!(initial_condition_array,0.0);	#	24	mRNA_SNAIL_SLUG
push!(initial_condition_array,0.1);	#	25	mRNA_tubulin
push!(initial_condition_array,0.1);	#	26	mRNA_APC_AXIN
push!(initial_condition_array,0.1);	#	27	mRNA_GSK3
push!(initial_condition_array,0.1);	#	28	mRNA_PI3K
push!(initial_condition_array,1.0);	#	29	PLCy
push!(initial_condition_array,0.0);	#	30	Active_PLCy
push!(initial_condition_array,1.0);	#	31	Calcineurin
push!(initial_condition_array,0.0);	#	32	Active_Calcineurin
push!(initial_condition_array,1.0);	#	33	NFATc
push!(initial_condition_array,0.0);	#	34	Active_NFATc
push!(initial_condition_array,1.0);	#	35	PKC
push!(initial_condition_array,0.0);	#	36	Active_PKC
push!(initial_condition_array,1.0);	#	37	RAF
push!(initial_condition_array,0.0);	#	38	RAF_P
push!(initial_condition_array,1.0);	#	39	Smad
push!(initial_condition_array,0.0);	#	40	pSmad
push!(initial_condition_array,1.0);	#	41	PI3K
push!(initial_condition_array,0.0);	#	42	PI3K_P
push!(initial_condition_array,1.0);	#	43	GSK3
push!(initial_condition_array,0.0);	#	44	GSK3_P
push!(initial_condition_array,1.0);	#	45	RAS
push!(initial_condition_array,0.0);	#	46	RAS_GTP
push!(initial_condition_array,1.0);	#	47	MAPK
push!(initial_condition_array,0.0);	#	48	pMAPK
push!(initial_condition_array,1.0);	#	49	AP1_SP1
push!(initial_condition_array,0.0);	#	50	AP1_SP1_P
push!(initial_condition_array,1.0);	#	51	Ecadherin
push!(initial_condition_array,1.0);	#	52	Bcatenin
push!(initial_condition_array,0.0);	#	53	Ecadherin_Bcatenin
push!(initial_condition_array,1.0);	#	54	APC_AXIN
push!(initial_condition_array,0.0);	#	55	APC_AXIN_Bcatenin
push!(initial_condition_array,1.0);	#	56	TCF4
push!(initial_condition_array,0.0);	#	57	Active_TCF4
push!(initial_condition_array,0.0);	#	58	LEF1
push!(initial_condition_array,0.0);	#	59	Active_LEF1
push!(initial_condition_array,100.0+10.0);	#	60	PPase
push!(initial_condition_array,1.0);	#	61	YREG1
push!(initial_condition_array,0.0);	#	62	SNAIL_SLUG
push!(initial_condition_array,0.0);	#	63	vimentin
push!(initial_condition_array,1.0);	#	64	tubulin
push!(initial_condition_array,270.0);	#	65	RNAP
push!(initial_condition_array,3599.9);	#	66	RIBOSOME
push!(initial_condition_array,0.0);	#	67	TGFB3_Ab
push!(initial_condition_array,0.0);	#	68	ss_virus
push!(initial_condition_array,0.0);	#	69	ecad_virus
push!(initial_condition_array,0.0);	#	70	U0126
push!(initial_condition_array,0.0);	#	71	vivit
push!(initial_condition_array,1.0);	#	72	gene_TCF4
push!(initial_condition_array,1.0);	#	73	gene_TGFB3
push!(initial_condition_array,1.0);	#	74	gene_YREG1
push!(initial_condition_array,1.0);	#	75	gene_LEF1
push!(initial_condition_array,1.0);	#	76	gene_vimentin
push!(initial_condition_array,1.0);	#	77	gene_R
push!(initial_condition_array,1.0);	#	78	gene_VEGFR
push!(initial_condition_array,1.0);	#	79	gene_NFATc
push!(initial_condition_array,1.0);	#	80	gene_MAPK
push!(initial_condition_array,0.1);	#	81	mRNA_TCF4
push!(initial_condition_array,0.0);	#	82	mRNA_TGFB3
push!(initial_condition_array,0.1);	#	83	mRNA_YREG1
push!(initial_condition_array,0.0);	#	84	mRNA_LEF1
push!(initial_condition_array,0.0);	#	85	mRNA_vimentin
push!(initial_condition_array,0.1);	#	86	mRNA_R
push!(initial_condition_array,0.1);	#	87	mRNA_VEGFR
push!(initial_condition_array,0.1);	#	88	mRNA_NFATc
push!(initial_condition_array,0.1);	#	89	mRNA_MAPK
push!(initial_condition_array,0.1);	#	90	R
push!(initial_condition_array,0.0);	#	91	TGFB12
push!(initial_condition_array,0.0);	#	92	TGFB3
push!(initial_condition_array,0.0);	#	93	VEGF
push!(initial_condition_array,0.1);	#	94	TGFB12_R
push!(initial_condition_array,0.1);	#	95	TGFB3_R
push!(initial_condition_array,0.1);	#	96	VEGFR
push!(initial_condition_array,0.0);	#	97	VEGF_VEGFR

# Formulate the time constant array -
time_constant_array = Float64[];
push!(time_constant_array,1.0);	#	1	 time constant: gene_Smad
push!(time_constant_array,1.0);	#	2	 time constant: gene_Ecadherin
push!(time_constant_array,1.0);	#	3	 time constant: gene_Bcatenin
push!(time_constant_array,1.0);	#	4	 time constant: gene_PLCy
push!(time_constant_array,1.0);	#	5	 time constant: gene_PKC
push!(time_constant_array,1.0);	#	6	 time constant: gene_RAF
push!(time_constant_array,1.0);	#	7	 time constant: gene_RAS
push!(time_constant_array,1.0);	#	8	 time constant: gene_AP1_SP1
push!(time_constant_array,1.0);	#	9	 time constant: gene_GSK3
push!(time_constant_array,1.0);	#	10	 time constant: gene_tubulin
push!(time_constant_array,1.0);	#	11	 time constant: gene_SNAIL_SLUG
push!(time_constant_array,1.0);	#	12	 time constant: gene_Calcineurin
push!(time_constant_array,1.0);	#	13	 time constant: gene_APC_AXIN
push!(time_constant_array,1.0);	#	14	 time constant: gene_PI3K
push!(time_constant_array,0.1);	#	15	 time constant: mRNA_Smad
push!(time_constant_array,0.1);	#	16	 time constant: mRNA_Ecadherin
push!(time_constant_array,0.1);	#	17	 time constant: mRNA_Bcatenin
push!(time_constant_array,0.1);	#	18	 time constant: mRNA_PLCy
push!(time_constant_array,0.1);	#	19	 time constant: mRNA_PKC
push!(time_constant_array,0.1);	#	20	 time constant: mRNA_RAF
push!(time_constant_array,0.1);	#	21	 time constant: mRNA_RAS
push!(time_constant_array,0.1);	#	22	 time constant: mRNA_AP1_SP1
push!(time_constant_array,0.1);	#	23	 time constant: mRNA_Calcineurin
push!(time_constant_array,0.1);	#	24	 time constant: mRNA_SNAIL_SLUG
push!(time_constant_array,0.1);	#	25	 time constant: mRNA_tubulin
push!(time_constant_array,0.1);	#	26	 time constant: mRNA_APC_AXIN
push!(time_constant_array,0.1);	#	27	 time constant: mRNA_GSK3
push!(time_constant_array,0.1);	#	28	 time constant: mRNA_PI3K
push!(time_constant_array,1.0);	#	29	 time constant: PLCy
push!(time_constant_array,1.0);	#	30	 time constant: Active_PLCy
push!(time_constant_array,1.0);	#	31	 time constant: Calcineurin
push!(time_constant_array,1.0);	#	32	 time constant: Active_Calcineurin
push!(time_constant_array,1.0);	#	33	 time constant: NFATc
push!(time_constant_array,1.0);	#	34	 time constant: Active_NFATc
push!(time_constant_array,1.0);	#	35	 time constant: PKC
push!(time_constant_array,1.0);	#	36	 time constant: Active_PKC
push!(time_constant_array,1.0);	#	37	 time constant: RAF
push!(time_constant_array,1.0);	#	38	 time constant: RAF_P
push!(time_constant_array,1.0);	#	39	 time constant: Smad
push!(time_constant_array,1.0);	#	40	 time constant: pSmad
push!(time_constant_array,1.0);	#	41	 time constant: PI3K
push!(time_constant_array,1.0);	#	42	 time constant: PI3K_P
push!(time_constant_array,1.0);	#	43	 time constant: GSK3
push!(time_constant_array,1.0);	#	44	 time constant: GSK3_P
push!(time_constant_array,1.0);	#	45	 time constant: RAS
push!(time_constant_array,1.0);	#	46	 time constant: RAS_GTP
push!(time_constant_array,1.0);	#	47	 time constant: MAPK
push!(time_constant_array,1.0);	#	48	 time constant: pMAPK
push!(time_constant_array,1.0);	#	49	 time constant: AP1_SP1
push!(time_constant_array,1.0);	#	50	 time constant: AP1_SP1_P
push!(time_constant_array,1.0);	#	51	 time constant: Ecadherin
push!(time_constant_array,1.0);	#	52	 time constant: Bcatenin
push!(time_constant_array,1.0);	#	53	 time constant: Ecadherin_Bcatenin
push!(time_constant_array,1.0);	#	54	 time constant: APC_AXIN
push!(time_constant_array,1.0);	#	55	 time constant: APC_AXIN_Bcatenin
push!(time_constant_array,1.0);	#	56	 time constant: TCF4
push!(time_constant_array,1.0);	#	57	 time constant: Active_TCF4
push!(time_constant_array,1.0);	#	58	 time constant: LEF1
push!(time_constant_array,1.0);	#	59	 time constant: Active_LEF1
push!(time_constant_array,1.0);	#	60	 time constant: PPase
push!(time_constant_array,1.0);	#	61	 time constant: YREG1
push!(time_constant_array,1.0);	#	62	 time constant: SNAIL_SLUG
push!(time_constant_array,1.0);	#	63	 time constant: vimentin
push!(time_constant_array,1.0);	#	64	 time constant: tubulin
push!(time_constant_array,1.0);	#	65	 time constant: RNAP
push!(time_constant_array,1.0);	#	66	 time constant: RIBOSOME
push!(time_constant_array,1.0);	#	67	 time constant: TGFB3_Ab
push!(time_constant_array,1.0);	#	68	 time constant: ss_virus
push!(time_constant_array,1.0);	#	69	 time constant: ecad_virus
push!(time_constant_array,1.0);	#	70	 time constant: U0126
push!(time_constant_array,1.0);	#	71	 time constant: vivit
push!(time_constant_array,1.0);	#	72	 time constant: gene_TCF4
push!(time_constant_array,1.0);	#	73	 time constant: gene_TGFB3
push!(time_constant_array,1.0);	#	74	 time constant: gene_YREG1
push!(time_constant_array,1.0);	#	75	 time constant: gene_LEF1
push!(time_constant_array,1.0);	#	76	 time constant: gene_vimentin
push!(time_constant_array,1.0);	#	77	 time constant: gene_R
push!(time_constant_array,1.0);	#	78	 time constant: gene_VEGFR
push!(time_constant_array,1.0);	#	79	 time constant: gene_NFATc
push!(time_constant_array,1.0);	#	80	 time constant: gene_MAPK
push!(time_constant_array,0.1);	#	81	 time constant: mRNA_TCF4
push!(time_constant_array,0.1);	#	82	 time constant: mRNA_TGFB3
push!(time_constant_array,0.1);	#	83	 time constant: mRNA_YREG1
push!(time_constant_array,0.1);	#	84	 time constant: mRNA_LEF1
push!(time_constant_array,0.1);	#	85	 time constant: mRNA_vimentin
push!(time_constant_array,0.1);	#	86	 time constant: mRNA_R
push!(time_constant_array,0.1);	#	87	 time constant: mRNA_VEGFR
push!(time_constant_array,0.1);	#	88	 time constant: mRNA_NFATc
push!(time_constant_array,0.1);	#	89	 time constant: mRNA_MAPK
push!(time_constant_array,1.0);	#	90	 time constant: R
push!(time_constant_array,1.0);	#	91	 time constant: TGFB12
push!(time_constant_array,1.0);	#	92	 time constant: TGFB3
push!(time_constant_array,1.0);	#	93	 time constant: VEGF
push!(time_constant_array,1.0);	#	94	 time constant: TGFB12_R
push!(time_constant_array,1.0);	#	95	 time constant: TGFB3_R
push!(time_constant_array,1.0);	#	96	 time constant: VEGFR
push!(time_constant_array,1.0);	#	97	 time constant: VEGF_VEGFR

# Formulate the rate constant array -
rate_constant_array = Float64[];
push!(rate_constant_array,1.0);	# 1	TGFB12_binding: TGFB12+R =([])=> TGFB12_R
push!(rate_constant_array,1.0);	# 2	TGFB12_binding: TGFB12_R =([])=> TGFB12+R (reverse)
push!(rate_constant_array,1.0);	# 3	TGFB3_binding: TGFB3+R =([])=> TGFB3_R
push!(rate_constant_array,1.0);	# 4	TGFB3_binding: TGFB3_R =([])=> TGFB3+R (reverse)
push!(rate_constant_array,1.0);	# 5	VEGF_binding: VEGF+VEGFR =([])=> VEGF_VEGFR
push!(rate_constant_array,1.0);	# 6	VEGF_binding: VEGF_VEGFR =([])=> VEGF+VEGFR (reverse)
push!(rate_constant_array,1.0);	# 7	VEGF_VEGFR_PLCy: PLCy =(VEGF_VEGFR)=> Active_PLCy
push!(rate_constant_array,1.0);	# 8	PLCy_activation_Calcineurin: Calcineurin =(Active_PLCy)=> Active_Calcineurin
push!(rate_constant_array,1.0);	# 9	Calcineurin_activation_NFATc: NFATc =(Active_Calcineurin)=> Active_NFATc
push!(rate_constant_array,1.0);	# 10	PLCy_activatation_PKC: PKC =(Active_PLCy)=> Active_PKC
push!(rate_constant_array,1.0);	# 11	PKC_activation_RAF: RAF =(Active_PKC)=> RAF_P
push!(rate_constant_array,1.0);	# 12	TGFB12_R_phosphorylation_Smad: Smad =(TGFB12_R)=> pSmad
push!(rate_constant_array,1.0);	# 13	TGFB3_R_phosphorylation_Smad: Smad =(TGFB3_R)=> pSmad
push!(rate_constant_array,1.0);	# 14	TGFB3_R_phosphorylation_PI3K: PI3K =(TGFB3_R)=> PI3K_P
push!(rate_constant_array,1.0);	# 15	PI3K_phosphorylation_GSK3: GSK3 =(PI3K_P)=> GSK3_P
push!(rate_constant_array,1.0);	# 16	TGFB12_R_activation_RAS: RAS =(TGFB12_R)=> RAS_GTP
push!(rate_constant_array,1.0);	# 17	RAS_GTP_phosphorylation_RAF: RAF =(RAS_GTP)=> RAF_P
push!(rate_constant_array,1.0);	# 18	RAF_P_phosphorylation_MAPK: MAPK =(RAF_P)=> pMAPK
push!(rate_constant_array,1.0);	# 19	MAPK_phosphorylation_Smad: Smad =(pMAPK)=> pSmad
push!(rate_constant_array,1.0);	# 20	MAPK_phosphorylation_AP1_SP1: AP1_SP1 =(pMAPK)=> AP1_SP1_P
push!(rate_constant_array,1.0);	# 21	Binding_Ecadherin_Bcatenin: Ecadherin+Bcatenin =([])=> Ecadherin_Bcatenin
push!(rate_constant_array,1.0);	# 22	Binding_Ecadherin_Bcatenin: Ecadherin_Bcatenin =([])=> Ecadherin+Bcatenin (reverse)
push!(rate_constant_array,1.0);	# 23	Binding_APC_AXIN_forw: APC_AXIN+Bcatenin =([])=> APC_AXIN_Bcatenin
push!(rate_constant_array,1.0);	# 24	Binding_APC_AXIN_revs: APC_AXIN_Bcatenin =([])=> APC_AXIN+Bcatenin
push!(rate_constant_array,1.0);	# 25	TCF4_Bcatenin_complex_formation: TCF4+Bcatenin =([])=> Active_TCF4
push!(rate_constant_array,1.0);	# 26	TCF4_Bcatenin_complex_formation: Active_TCF4 =([])=> TCF4+Bcatenin (reverse)
push!(rate_constant_array,1.0);	# 27	LEF1_Bcatenin_complex_formation: LEF1+Bcatenin =([])=> Active_LEF1
push!(rate_constant_array,1.0);	# 28	LEF1_Bcatenin_complex_formation: Active_LEF1 =([])=> LEF1+Bcatenin (reverse)
push!(rate_constant_array,0.9);	# 29	deactivation_pSmad: pSmad =(PPase)=> Smad
push!(rate_constant_array,1.0);	# 30	deactivation_RAF_P: RAF_P =(PPase)=> RAF
push!(rate_constant_array,1.0);	# 31	deactivation_RAS_GTP: RAS_GTP =([])=> RAS
push!(rate_constant_array,1.0);	# 32	deactivation_pMAPK: pMAPK =(PPase)=> MAPK
push!(rate_constant_array,1.0);	# 33	deactivation_AP1_SP1: AP1_SP1_P =(PPase)=> AP1_SP1
push!(rate_constant_array,25.0);	# 34	deactivation_Active_PLCy: Active_PLCy =([])=> PLCy
push!(rate_constant_array,1.0);	# 35	deactivation_Active_PKC: Active_PKC =([])=> PKC
push!(rate_constant_array,100.0);	# 36	deactivation_Active_Calcineurin: Active_Calcineurin =([])=> Calcineurin
push!(rate_constant_array,100.0);	# 37	deactivation_Active_NFATc: Active_NFATc =([])=> NFATc
push!(rate_constant_array,1.0);	# 38	deactivation_GSK3_P: GSK3_P =(PPase)=> GSK3
push!(rate_constant_array,0.1);	# 39	degradation_TGFB12_R: TGFB12_R =([])=> []
push!(rate_constant_array,0.1);	# 40	degradation_TGFB3_R: TGFB3_R =([])=> []
push!(rate_constant_array,0.1);	# 41	degradation_TGFB3: TGFB3 =([])=> []
push!(rate_constant_array,0.1);	# 42	VEGF_VEGFR_degradation: VEGF_VEGFR =([])=> []
push!(rate_constant_array,0.1);	# 43	Bcatenin_degradation: Bcatenin =([])=> []
push!(rate_constant_array,0.1);	# 44	Ecadherin_degrdation: Ecadherin =([])=> []
push!(rate_constant_array,0.1);	# 45	Smad_degradation: Smad =([])=> []
push!(rate_constant_array,0.1);	# 46	TCF4_degradation: TCF4 =([])=> []
push!(rate_constant_array,0.1);	# 47	Active_TCF4_degradation: Active_TCF4 =([])=> []
push!(rate_constant_array,0.1);	# 48	YREG1_degradation: YREG1 =([])=> []
push!(rate_constant_array,0.1);	# 49	SNAIL_SLUG_degradation: SNAIL_SLUG =([])=> []
push!(rate_constant_array,0.1);	# 50	LEF1_degradation: LEF1 =([])=> []
push!(rate_constant_array,0.1);	# 51	Active_LEF1_degradation: Active_LEF1 =([])=> []
push!(rate_constant_array,0.1);	# 52	vimentin_degradation: vimentin =([])=> []
push!(rate_constant_array,0.1);	# 53	VEGFR_degradation: VEGFR =([])=> []
push!(rate_constant_array,0.1);	# 54	NFATc_degradation: NFATc =([])=> []
push!(rate_constant_array,0.1);	# 55	Active_NFATc_degradation: Active_NFATc =([])=> []
push!(rate_constant_array,0.1);	# 56	degradation_MAPK: MAPK =([])=> []
push!(rate_constant_array,0.1);	# 57	degradation_R: R =([])=> []
push!(rate_constant_array,0.1);	# 58	Ecadherin_Bcatenin_degradation: Ecadherin_Bcatenin =([])=> []
push!(rate_constant_array,0.1);	# 59	AP1_SP1_degradation: AP1_SP1 =([])=> []
push!(rate_constant_array,0.1);	# 60	AP1_SP1_P_degradation: AP1_SP1_P =([])=> []
push!(rate_constant_array,0.1);	# 61	RAS_degradation: RAS =([])=> []
push!(rate_constant_array,0.1);	# 62	RAF_degradation: RAF =([])=> []
push!(rate_constant_array,0.1);	# 63	RAS_GTP_degradation: RAS_GTP =([])=> []
push!(rate_constant_array,0.1);	# 64	RAF_P_degradation: RAF_P =([])=> []
push!(rate_constant_array,0.1);	# 65	pMAPK_degradation: pMAPK =([])=> []
push!(rate_constant_array,0.1);	# 66	pSmad_degradation: pSmad =([])=> []
push!(rate_constant_array,0.1);	# 67	PLCy_degradation: PLCy =([])=> []
push!(rate_constant_array,0.1);	# 68	Active_PLCy_degradation: Active_PLCy =([])=> []
push!(rate_constant_array,0.1);	# 69	PKC_degradation: PKC =([])=> []
push!(rate_constant_array,0.1);	# 70	Active_PKC_degradation: Active_PKC =([])=> []
push!(rate_constant_array,0.1);	# 71	Calcineurin_degradation: Calcineurin =([])=> []
push!(rate_constant_array,0.1);	# 72	Active_Calcineurin_degradation: Active_Calcineurin =([])=> []
push!(rate_constant_array,0.1);	# 73	degradation_APC_AXIN: APC_AXIN =([])=> []
push!(rate_constant_array,0.1);	# 74	degradation_GSK3_P: GSK3_P =([])=> []
push!(rate_constant_array,0.1);	# 75	degradation_GSK3: GSK3 =([])=> []
push!(rate_constant_array,0.1);	# 76	degradation_PI3K: PI3K =([])=> []
push!(rate_constant_array,0.1);	# 77	degradation_APC_AXIN_Bcatenin: APC_AXIN_Bcatenin =([])=> []
push!(rate_constant_array,0.1);	# 78	tubulin_degradation: tubulin =([])=> []
push!(rate_constant_array,0.1);	# 79	generation_RNAP: RNAP =([])=> []
push!(rate_constant_array,1.0);	# 80	generation_RNAP: [] =([])=> RNAP (reverse)
push!(rate_constant_array,0.1);	# 81	generation_RIBOSOME: RIBOSOME =([])=> []
push!(rate_constant_array,1.0);	# 82	generation_RIBOSOME: [] =([])=> RIBOSOME (reverse)
push!(rate_constant_array,0.1);	# 83	generation_PPase: PPase =([])=> []
push!(rate_constant_array,1.0);	# 84	generation_PPase: [] =([])=> PPase (reverse)
push!(rate_constant_array,0.1);	# 85	addition_TGFB3_Ab: TGFB3_Ab =([])=> []
push!(rate_constant_array,1.0);	# 86	addition_TGFB3_Ab: [] =([])=> TGFB3_Ab (reverse)
push!(rate_constant_array,0.1);	# 87	addition_ss_virus: ss_virus =([])=> []
push!(rate_constant_array,1.0);	# 88	addition_ss_virus: [] =([])=> ss_virus (reverse)
push!(rate_constant_array,0.1);	# 89	addition_ecad_virus: ecad_virus =([])=> []
push!(rate_constant_array,1.0);	# 90	addition_ecad_virus: [] =([])=> ecad_virus (reverse)
push!(rate_constant_array,1.0);	# 91	addition_TGFb12: [] =([])=> TGFB12
push!(rate_constant_array,0.1);	# 92	addition_TGFb12: TGFB12 =([])=> [] (reverse)
push!(rate_constant_array,1.0);	# 93	addition_VEGF: [] =([])=> VEGF
push!(rate_constant_array,0.1);	# 94	addition_VEGF: VEGF =([])=> [] (reverse)
push!(rate_constant_array,1.0);	# 95	addition_TGFB3: [] =([])=> TGFB3
push!(rate_constant_array,0.1);	# 96	addition_TGFB3: TGFB3 =([])=> [] (reverse)
push!(rate_constant_array,0.1);	# 97	addition_U0126: U0126 =([])=> []
push!(rate_constant_array,1.0);	# 98	addition_U0126: [] =([])=> U0126 (reverse)
push!(rate_constant_array,0.1);	# 99	addition_vivit: vivit =([])=> []
push!(rate_constant_array,1.0);	# 100	addition_vivit: [] =([])=> vivit (reverse)
push!(rate_constant_array,0.1);	# 101	degradation_mRNA_Smad: mRNA_Smad = []
push!(rate_constant_array,0.1);	# 102	degradation_mRNA_Ecadherin: mRNA_Ecadherin = []
push!(rate_constant_array,0.1);	# 103	degradation_mRNA_Bcatenin: mRNA_Bcatenin = []
push!(rate_constant_array,0.1);	# 104	degradation_mRNA_TCF4: mRNA_TCF4 = []
push!(rate_constant_array,0.1);	# 105	degradation_mRNA_TGFB3: mRNA_TGFB3 = []
push!(rate_constant_array,0.1);	# 106	degradation_mRNA_YREG1: mRNA_YREG1 = []
push!(rate_constant_array,0.1);	# 107	degradation_mRNA_LEF1: mRNA_LEF1 = []
push!(rate_constant_array,0.1);	# 108	degradation_mRNA_vimentin: mRNA_vimentin = []
push!(rate_constant_array,0.1);	# 109	degradation_mRNA_VEGFR: mRNA_VEGFR = []
push!(rate_constant_array,0.1);	# 110	degradation_mRNA_NFATc: mRNA_NFATc = []
push!(rate_constant_array,0.1);	# 111	degradation_mRNA_MAPK: mRNA_MAPK = []
push!(rate_constant_array,0.1);	# 112	degradation_mRNA_R: mRNA_R = []
push!(rate_constant_array,0.1);	# 113	degradation_mRNA_PLCy: mRNA_PLCy = []
push!(rate_constant_array,0.1);	# 114	degradation_mRNA_PKC: mRNA_PKC = []
push!(rate_constant_array,0.1);	# 115	degradation_mRNA_RAF: mRNA_RAF = []
push!(rate_constant_array,0.1);	# 116	degradation_mRNA_RAS: mRNA_RAS = []
push!(rate_constant_array,0.1);	# 117	degradation_mRNA_AP1_SP1: mRNA_AP1_SP1 = []
push!(rate_constant_array,0.1);	# 118	degradation_mRNA_Calcineurin: mRNA_Calcineurin = []
push!(rate_constant_array,0.1);	# 119	degradation_mRNA_SNAIL_SLUG: mRNA_SNAIL_SLUG = []
push!(rate_constant_array,0.1);	# 120	degradation_mRNA_tubulin: mRNA_tubulin = []
push!(rate_constant_array,0.1);	# 121	degradation_mRNA_APC_AXIN: mRNA_APC_AXIN = []
push!(rate_constant_array,0.1);	# 122	degradation_mRNA_GSK3: mRNA_GSK3 = []
push!(rate_constant_array,0.1);	# 123	degradation_mRNA_PI3K: mRNA_PI3K = []
push!(rate_constant_array,1.0);	# 124	induction_gene_Smad: gene_Smad = mRNA_Smad
push!(rate_constant_array,1.0);	# 125	induction_gene_Ecadherin: gene_Ecadherin = mRNA_Ecadherin
push!(rate_constant_array,1.0);	# 126	induction_gene_Bcatenin: gene_Bcatenin = mRNA_Bcatenin
push!(rate_constant_array,1.0);	# 127	induction_gene_TCF4: gene_TCF4 = mRNA_TCF4
push!(rate_constant_array,1.0);	# 128	induction_gene_TGFB3: gene_TGFB3 = mRNA_TGFB3
push!(rate_constant_array,1.0);	# 129	induction_gene_YREG1: gene_YREG1 = mRNA_YREG1
push!(rate_constant_array,1.0);	# 130	induction_gene_LEF1: gene_LEF1 = mRNA_LEF1
push!(rate_constant_array,1.0);	# 131	induction_gene_R: gene_R = mRNA_R
push!(rate_constant_array,1.0);	# 132	induction_gene_vimentin: gene_vimentin = mRNA_vimentin
push!(rate_constant_array,1.0);	# 133	induction_gene_VEGFR: gene_VEGFR = mRNA_VEGFR
push!(rate_constant_array,1.0);	# 134	induction_gene_NFATc: gene_NFATc = mRNA_NFATc
push!(rate_constant_array,1.0);	# 135	induction_gene_PLCy: gene_PLCy = mRNA_PLCy
push!(rate_constant_array,1.0);	# 136	induction_gene_PKC: gene_PKC = mRNA_PKC
push!(rate_constant_array,1.0);	# 137	induction_gene_RAF: gene_RAF = mRNA_RAF
push!(rate_constant_array,1.0);	# 138	induction_gene_RAS: gene_RAS = mRNA_RAS
push!(rate_constant_array,1.0);	# 139	induction_gene_MAPK: gene_MAPK = mRNA_MAPK
push!(rate_constant_array,1.0);	# 140	induction_gene_AP1_SP1: gene_AP1_SP1 = mRNA_AP1_SP1
push!(rate_constant_array,1.0);	# 141	induction_gene_GSK3: gene_GSK3 = mRNA_GSK3
push!(rate_constant_array,1.0);	# 142	induction_gene_tubulin: gene_tubulin = mRNA_tubulin
push!(rate_constant_array,1.0);	# 143	induction_gene_SNAIL_SLUG: gene_SNAIL_SLUG = mRNA_SNAIL_SLUG
push!(rate_constant_array,1.0);	# 144	induction_gene_Calcineurin: gene_Calcineurin = mRNA_Calcineurin
push!(rate_constant_array,1.0);	# 145	induction_gene_APC_AXIN: gene_APC_AXIN = mRNA_APC_AXIN
push!(rate_constant_array,1.0);	# 146	induction_gene_PI3K: gene_PI3K = mRNA_PI3K
push!(rate_constant_array,1.0);	# 147	translation_mRNA_Smad: mRNA_Smad = Smad
push!(rate_constant_array,1.0);	# 148	translation_mRNA_Ecadherin: mRNA_Ecadherin = Ecadherin
push!(rate_constant_array,1.0);	# 149	translation_mRNA_Bcatenin: mRNA_Bcatenin = Bcatenin
push!(rate_constant_array,1.0);	# 150	translation_mRNA_TCF4: mRNA_TCF4 = TCF4
push!(rate_constant_array,1.0);	# 151	translation_mRNA_TGFB3: mRNA_TGFB3 = TGFB3
push!(rate_constant_array,1.0);	# 152	translation_mRNA_YREG1: mRNA_YREG1 = YREG1
push!(rate_constant_array,1.0);	# 153	translation_mRNA_LEF1: mRNA_LEF1 = LEF1
push!(rate_constant_array,1.0);	# 154	translation_mRNA_vimentin: mRNA_vimentin = vimentin
push!(rate_constant_array,1.0);	# 155	translation_mRNA_VEGFR: mRNA_VEGFR = VEGFR
push!(rate_constant_array,1.0);	# 156	translation_mRNA_NFATc: mRNA_NFATc = NFATc
push!(rate_constant_array,1.0);	# 157	translation_mRNA_MAPK: mRNA_MAPK = MAPK
push!(rate_constant_array,1.0);	# 158	translation_mRNA_R: mRNA_R = R
push!(rate_constant_array,1.0);	# 159	translation_mRNA_PLCy: mRNA_PLCy = PLCy
push!(rate_constant_array,1.0);	# 160	translation_mRNA_PKC: mRNA_PKC = PKC
push!(rate_constant_array,1.0);	# 161	translation_mRNA_RAF: mRNA_RAF = RAF
push!(rate_constant_array,1.0);	# 162	translation_mRNA_RAS: mRNA_RAS = RAS
push!(rate_constant_array,1.0);	# 163	translation_mRNA_AP1_SP1: mRNA_AP1_SP1 = AP1_SP1
push!(rate_constant_array,1.0);	# 164	translation_mRNA_Calcineurin: mRNA_Calcineurin = Calcineurin
push!(rate_constant_array,1.0);	# 165	translation_mRNA_SNAIL_SLUG: mRNA_SNAIL_SLUG = SNAIL_SLUG
push!(rate_constant_array,1.0);	# 166	translation_mRNA_tubulin: mRNA_tubulin = tubulin
push!(rate_constant_array,1.0);	# 167	translation_mRNA_APC_AXIN: mRNA_APC_AXIN = APC_AXIN
push!(rate_constant_array,1.0);	# 168	translation_mRNA_GSK3: mRNA_GSK3 = GSK3
push!(rate_constant_array,1.0);	# 169	translation_mRNA_PI3K: mRNA_PI3K = PI3K

# Formulate the saturation constant array -
saturation_constant_array = zeros(NREACTIONS,NSPECIES);
saturation_constant_array[1,91] = 1.0;	# TGFB12_binding: TGFB12+R =([])=> TGFB12_R	 species: TGFB12
saturation_constant_array[1,90] = 1.0;	# TGFB12_binding: TGFB12+R =([])=> TGFB12_R	 species: R
saturation_constant_array[2,94] = 1.0;	# TGFB12_binding: TGFB12_R =([])=> TGFB12+R (reverse)	 species: TGFB12_R
saturation_constant_array[3,92] = 1.0;	# TGFB3_binding: TGFB3+R =([])=> TGFB3_R	 species: TGFB3
saturation_constant_array[3,90] = 1.0;	# TGFB3_binding: TGFB3+R =([])=> TGFB3_R	 species: R
saturation_constant_array[4,95] = 1.0;	# TGFB3_binding: TGFB3_R =([])=> TGFB3+R (reverse)	 species: TGFB3_R
saturation_constant_array[5,93] = 1.0;	# VEGF_binding: VEGF+VEGFR =([])=> VEGF_VEGFR	 species: VEGF
saturation_constant_array[5,96] = 1.0;	# VEGF_binding: VEGF+VEGFR =([])=> VEGF_VEGFR	 species: VEGFR
saturation_constant_array[6,97] = 1.0;	# VEGF_binding: VEGF_VEGFR =([])=> VEGF+VEGFR (reverse)	 species: VEGF_VEGFR
saturation_constant_array[7,29] = 1.0;	# VEGF_VEGFR_PLCy: PLCy =(VEGF_VEGFR)=> Active_PLCy	 species: PLCy
saturation_constant_array[8,31] = 1.0;	# PLCy_activation_Calcineurin: Calcineurin =(Active_PLCy)=> Active_Calcineurin	 species: Calcineurin
saturation_constant_array[9,33] = 1.0;	# Calcineurin_activation_NFATc: NFATc =(Active_Calcineurin)=> Active_NFATc	 species: NFATc
saturation_constant_array[10,35] = 1.0;	# PLCy_activatation_PKC: PKC =(Active_PLCy)=> Active_PKC	 species: PKC
saturation_constant_array[11,37] = 1.0;	# PKC_activation_RAF: RAF =(Active_PKC)=> RAF_P	 species: RAF
saturation_constant_array[12,39] = 1.0;	# TGFB12_R_phosphorylation_Smad: Smad =(TGFB12_R)=> pSmad	 species: Smad
saturation_constant_array[13,39] = 1.0;	# TGFB3_R_phosphorylation_Smad: Smad =(TGFB3_R)=> pSmad	 species: Smad
saturation_constant_array[14,41] = 1.0;	# TGFB3_R_phosphorylation_PI3K: PI3K =(TGFB3_R)=> PI3K_P	 species: PI3K
saturation_constant_array[15,43] = 1.0;	# PI3K_phosphorylation_GSK3: GSK3 =(PI3K_P)=> GSK3_P	 species: GSK3
saturation_constant_array[16,45] = 1.0;	# TGFB12_R_activation_RAS: RAS =(TGFB12_R)=> RAS_GTP	 species: RAS
saturation_constant_array[17,37] = 1.0;	# RAS_GTP_phosphorylation_RAF: RAF =(RAS_GTP)=> RAF_P	 species: RAF
saturation_constant_array[18,47] = 1.0;	# RAF_P_phosphorylation_MAPK: MAPK =(RAF_P)=> pMAPK	 species: MAPK
saturation_constant_array[19,39] = 1.0;	# MAPK_phosphorylation_Smad: Smad =(pMAPK)=> pSmad	 species: Smad
saturation_constant_array[20,49] = 1.0;	# MAPK_phosphorylation_AP1_SP1: AP1_SP1 =(pMAPK)=> AP1_SP1_P	 species: AP1_SP1
saturation_constant_array[21,51] = 2.0;	# Binding_Ecadherin_Bcatenin: Ecadherin+Bcatenin =([])=> Ecadherin_Bcatenin	 species: Ecadherin
saturation_constant_array[21,52] = 0.5;	# Binding_Ecadherin_Bcatenin: Ecadherin+Bcatenin =([])=> Ecadherin_Bcatenin	 species: Bcatenin
saturation_constant_array[22,53] = 1.0;	# Binding_Ecadherin_Bcatenin: Ecadherin_Bcatenin =([])=> Ecadherin+Bcatenin (reverse)	 species: Ecadherin_Bcatenin
saturation_constant_array[23,54] = 1.0;	# Binding_APC_AXIN_forw: APC_AXIN+Bcatenin =([])=> APC_AXIN_Bcatenin	 species: APC_AXIN
saturation_constant_array[23,52] = 1.0;	# Binding_APC_AXIN_forw: APC_AXIN+Bcatenin =([])=> APC_AXIN_Bcatenin	 species: Bcatenin
saturation_constant_array[24,55] = 1.0;	# Binding_APC_AXIN_revs: APC_AXIN_Bcatenin =([])=> APC_AXIN+Bcatenin	 species: APC_AXIN_Bcatenin
saturation_constant_array[25,56] = 1.0;	# TCF4_Bcatenin_complex_formation: TCF4+Bcatenin =([])=> Active_TCF4	 species: TCF4
saturation_constant_array[25,52] = 1.0;	# TCF4_Bcatenin_complex_formation: TCF4+Bcatenin =([])=> Active_TCF4	 species: Bcatenin
saturation_constant_array[26,57] = 1.0;	# TCF4_Bcatenin_complex_formation: Active_TCF4 =([])=> TCF4+Bcatenin (reverse)	 species: Active_TCF4
saturation_constant_array[27,58] = 1.0;	# LEF1_Bcatenin_complex_formation: LEF1+Bcatenin =([])=> Active_LEF1	 species: LEF1
saturation_constant_array[27,52] = 1.0;	# LEF1_Bcatenin_complex_formation: LEF1+Bcatenin =([])=> Active_LEF1	 species: Bcatenin
saturation_constant_array[28,59] = 1.0;	# LEF1_Bcatenin_complex_formation: Active_LEF1 =([])=> LEF1+Bcatenin (reverse)	 species: Active_LEF1
saturation_constant_array[29,40] = 1.0;	# deactivation_pSmad: pSmad =(PPase)=> Smad	 species: pSmad
saturation_constant_array[30,38] = 1.0;	# deactivation_RAF_P: RAF_P =(PPase)=> RAF	 species: RAF_P
saturation_constant_array[31,46] = 1.0;	# deactivation_RAS_GTP: RAS_GTP =([])=> RAS	 species: RAS_GTP
saturation_constant_array[32,48] = 1.0;	# deactivation_pMAPK: pMAPK =(PPase)=> MAPK	 species: pMAPK
saturation_constant_array[33,50] = 1.0;	# deactivation_AP1_SP1: AP1_SP1_P =(PPase)=> AP1_SP1	 species: AP1_SP1_P
saturation_constant_array[34,30] = 1.0;	# deactivation_Active_PLCy: Active_PLCy =([])=> PLCy	 species: Active_PLCy
saturation_constant_array[35,36] = 1.0;	# deactivation_Active_PKC: Active_PKC =([])=> PKC	 species: Active_PKC
saturation_constant_array[36,32] = 1.0;	# deactivation_Active_Calcineurin: Active_Calcineurin =([])=> Calcineurin	 species: Active_Calcineurin
saturation_constant_array[37,34] = 1.0;	# deactivation_Active_NFATc: Active_NFATc =([])=> NFATc	 species: Active_NFATc
saturation_constant_array[38,44] = 1.0;	# deactivation_GSK3_P: GSK3_P =(PPase)=> GSK3	 species: GSK3_P

# Formulate control parameter array -
control_parameter_array = zeros(20,2);
# TGFB3_inhibition: TGFB3_Ab inhibition; TGFB3_binding
control_parameter_array[1,1] = 0.1;	#	1 Gain: 	TGFB3_inhibition target: TGFB3_binding actor: TGFB3_Ab type: inhibition
control_parameter_array[1,2] = 1.0;	#	1 Order: 	TGFB3_inhibition target: TGFB3_binding actor: TGFB3_Ab type: inhibition

# NFATc_inhibition: vivit inhibition; Calcineurin_activation_NFATc
control_parameter_array[2,1] = 0.1;	#	2 Gain: 	NFATc_inhibition target: Calcineurin_activation_NFATc actor: vivit type: inhibition
control_parameter_array[2,2] = 1.0;	#	2 Order: 	NFATc_inhibition target: Calcineurin_activation_NFATc actor: vivit type: inhibition

# GSK3_P_inhibition_Binding_APC_AXIN_forw: GSK3_P inhibition; Binding_APC_AXIN_forw
control_parameter_array[3,1] = 0.1;	#	3 Gain: 	GSK3_P_inhibition_Binding_APC_AXIN_forw target: Binding_APC_AXIN_forw actor: GSK3_P type: inhibition
control_parameter_array[3,2] = 1.0;	#	3 Order: 	GSK3_P_inhibition_Binding_APC_AXIN_forw target: Binding_APC_AXIN_forw actor: GSK3_P type: inhibition

# Active_TCF4_inhibition_Ecad: Ecadherin_Bcatenin inhibition; TCF4_Bcatenin_complex_formation
control_parameter_array[4,1] = 0.1;	#	4 Gain: 	Active_TCF4_inhibition_Ecad target: TCF4_Bcatenin_complex_formation actor: Ecadherin_Bcatenin type: inhibition
control_parameter_array[4,2] = 1.0;	#	4 Order: 	Active_TCF4_inhibition_Ecad target: TCF4_Bcatenin_complex_formation actor: Ecadherin_Bcatenin type: inhibition

# Active_TCF4_inhibition_APC: APC_AXIN_Bcatenin inhibition; TCF4_Bcatenin_complex_formation
control_parameter_array[5,1] = 0.1;	#	5 Gain: 	Active_TCF4_inhibition_APC target: TCF4_Bcatenin_complex_formation actor: APC_AXIN_Bcatenin type: inhibition
control_parameter_array[5,2] = 1.0;	#	5 Order: 	Active_TCF4_inhibition_APC target: TCF4_Bcatenin_complex_formation actor: APC_AXIN_Bcatenin type: inhibition

# activation_Ecadherin_by_Active_NFATc: Active_NFATc induction induction_gene_Ecadherin
control_parameter_array[6,1] = 0.1;	#	6 Gain: 	activation_Ecadherin_by_Active_NFATc target: induction_gene_Ecadherin actor: Active_NFATc type: induction
control_parameter_array[6,2] = 1.0;	#	6 Order: 	activation_Ecadherin_by_Active_NFATc target: induction_gene_Ecadherin actor: Active_NFATc type: induction

# repression_Ecadherin_by_pSmad: pSmad repression induction_gene_Ecadherin
control_parameter_array[7,1] = 0.1;	#	7 Gain: 	repression_Ecadherin_by_pSmad target: induction_gene_Ecadherin actor: pSmad type: repression
control_parameter_array[7,2] = 1.0;	#	7 Order: 	repression_Ecadherin_by_pSmad target: induction_gene_Ecadherin actor: pSmad type: repression

# repression_Ecadherin_by_SNAIL_SLUG: SNAIL_SLUG repression induction_gene_Ecadherin
control_parameter_array[8,1] = 0.1;	#	8 Gain: 	repression_Ecadherin_by_SNAIL_SLUG target: induction_gene_Ecadherin actor: SNAIL_SLUG type: repression
control_parameter_array[8,2] = 1.0;	#	8 Order: 	repression_Ecadherin_by_SNAIL_SLUG target: induction_gene_Ecadherin actor: SNAIL_SLUG type: repression

# repression_Ecadherin_by_Active_LEF1: Active_LEF1 repression induction_gene_Ecadherin
control_parameter_array[9,1] = 0.1;	#	9 Gain: 	repression_Ecadherin_by_Active_LEF1 target: induction_gene_Ecadherin actor: Active_LEF1 type: repression
control_parameter_array[9,2] = 1.0;	#	9 Order: 	repression_Ecadherin_by_Active_LEF1 target: induction_gene_Ecadherin actor: Active_LEF1 type: repression

# activation_Ecadherin_by_ecad_virus: ecad_virus induction induction_gene_Ecadherin
control_parameter_array[10,1] = 0.1;	#	10 Gain: 	activation_Ecadherin_by_ecad_virus target: induction_gene_Ecadherin actor: ecad_virus type: induction
control_parameter_array[10,2] = 1.0;	#	10 Order: 	activation_Ecadherin_by_ecad_virus target: induction_gene_Ecadherin actor: ecad_virus type: induction

# activation_TGFB3_by_TCF4: Active_TCF4 induction induction_gene_TGFB3
control_parameter_array[11,1] = 0.1;	#	11 Gain: 	activation_TGFB3_by_TCF4 target: induction_gene_TGFB3 actor: Active_TCF4 type: induction
control_parameter_array[11,2] = 1.0;	#	11 Order: 	activation_TGFB3_by_TCF4 target: induction_gene_TGFB3 actor: Active_TCF4 type: induction

# repression_TGFB3_by_YREG1: YREG1 repression induction_gene_TGFB3
control_parameter_array[12,1] = 0.1;	#	12 Gain: 	repression_TGFB3_by_YREG1 target: induction_gene_TGFB3 actor: YREG1 type: repression
control_parameter_array[12,2] = 1.0;	#	12 Order: 	repression_TGFB3_by_YREG1 target: induction_gene_TGFB3 actor: YREG1 type: repression

# repression_YREG1_by_SNAIL_SLUG: SNAIL_SLUG repression induction_gene_YREG1
control_parameter_array[13,1] = 0.1;	#	13 Gain: 	repression_YREG1_by_SNAIL_SLUG target: induction_gene_YREG1 actor: SNAIL_SLUG type: repression
control_parameter_array[13,2] = 1.0;	#	13 Order: 	repression_YREG1_by_SNAIL_SLUG target: induction_gene_YREG1 actor: SNAIL_SLUG type: repression

# activation_LEF1_by_pSmad: pSmad induction induction_gene_LEF1
control_parameter_array[14,1] = 0.1;	#	14 Gain: 	activation_LEF1_by_pSmad target: induction_gene_LEF1 actor: pSmad type: induction
control_parameter_array[14,2] = 1.0;	#	14 Order: 	activation_LEF1_by_pSmad target: induction_gene_LEF1 actor: pSmad type: induction

# repression_LEF1_by_YREG1: YREG1 repression induction_gene_LEF1
control_parameter_array[15,1] = 0.1;	#	15 Gain: 	repression_LEF1_by_YREG1 target: induction_gene_LEF1 actor: YREG1 type: repression
control_parameter_array[15,2] = 1.0;	#	15 Order: 	repression_LEF1_by_YREG1 target: induction_gene_LEF1 actor: YREG1 type: repression

# activation_vimentin_by_Active_LEF1: Active_LEF1 induction induction_gene_vimentin
control_parameter_array[16,1] = 0.1;	#	16 Gain: 	activation_vimentin_by_Active_LEF1 target: induction_gene_vimentin actor: Active_LEF1 type: induction
control_parameter_array[16,2] = 1.0;	#	16 Order: 	activation_vimentin_by_Active_LEF1 target: induction_gene_vimentin actor: Active_LEF1 type: induction

# activation_vimentin_by_AP1_SP1_P: AP1_SP1_P induction induction_gene_vimentin
control_parameter_array[17,1] = 0.1;	#	17 Gain: 	activation_vimentin_by_AP1_SP1_P target: induction_gene_vimentin actor: AP1_SP1_P type: induction
control_parameter_array[17,2] = 1.0;	#	17 Order: 	activation_vimentin_by_AP1_SP1_P target: induction_gene_vimentin actor: AP1_SP1_P type: induction

# activation_SNAIL_SLUG_by_AP1_SP1_P: AP1_SP1_P induction induction_gene_SNAIL_SLUG
control_parameter_array[18,1] = 0.1;	#	18 Gain: 	activation_SNAIL_SLUG_by_AP1_SP1_P target: induction_gene_SNAIL_SLUG actor: AP1_SP1_P type: induction
control_parameter_array[18,2] = 1.0;	#	18 Order: 	activation_SNAIL_SLUG_by_AP1_SP1_P target: induction_gene_SNAIL_SLUG actor: AP1_SP1_P type: induction

# activation_SNAIL_SLUG_by_ss_virus: ss_virus induction induction_gene_SNAIL_SLUG
control_parameter_array[19,1] = 0.1;	#	19 Gain: 	activation_SNAIL_SLUG_by_ss_virus target: induction_gene_SNAIL_SLUG actor: ss_virus type: induction
control_parameter_array[19,2] = 1.0;	#	19 Order: 	activation_SNAIL_SLUG_by_ss_virus target: induction_gene_SNAIL_SLUG actor: ss_virus type: induction

# Dilution selection array -
dilution_selection_array = ones(Int,NSPECIES)
dilution_selection_array[1] = 0.0;	#	gene_Smad
dilution_selection_array[2] = 0.0;	#	gene_Ecadherin
dilution_selection_array[3] = 0.0;	#	gene_Bcatenin
dilution_selection_array[4] = 0.0;	#	gene_PLCy
dilution_selection_array[5] = 0.0;	#	gene_PKC
dilution_selection_array[6] = 0.0;	#	gene_RAF
dilution_selection_array[7] = 0.0;	#	gene_RAS
dilution_selection_array[8] = 0.0;	#	gene_AP1_SP1
dilution_selection_array[9] = 0.0;	#	gene_GSK3
dilution_selection_array[10] = 0.0;	#	gene_tubulin
dilution_selection_array[11] = 0.0;	#	gene_SNAIL_SLUG
dilution_selection_array[12] = 0.0;	#	gene_Calcineurin
dilution_selection_array[13] = 0.0;	#	gene_APC_AXIN
dilution_selection_array[14] = 0.0;	#	gene_PI3K
dilution_selection_array[72] = 0.0;	#	gene_TCF4
dilution_selection_array[73] = 0.0;	#	gene_TGFB3
dilution_selection_array[74] = 0.0;	#	gene_YREG1
dilution_selection_array[75] = 0.0;	#	gene_LEF1
dilution_selection_array[76] = 0.0;	#	gene_vimentin
dilution_selection_array[77] = 0.0;	#	gene_R
dilution_selection_array[78] = 0.0;	#	gene_VEGFR
dilution_selection_array[79] = 0.0;	#	gene_NFATc
dilution_selection_array[80] = 0.0;	#	gene_MAPK

# Set the maximum specific growth rate -
maximum_specific_growth_rate = 0.03; # hr^-1

# Set Bionumbers scaling constants
average_transcript_length = 15000 # n
average_protein_length = 5000     # aa

max_translation_rate = 5          # aa/sec
max_transcription_rate = 6.0       # n/sec

av_number = 6.02e23               # number/mol

hl60_diameter = 12                # \mum
fraction_nucleus = 0.49           # dimensionless
VOLUME = ((1-fraction_nucleus)*(1/6)*(3.14159)*(hl60_diameter)^3)*(1e-15)

kcat_transcription = max_transcription_rate*(3600/average_transcript_length)  # hr^-1
kcat_translation = max_translation_rate*(3600/average_protein_length)  # hr^-1

saturation_transcription = (46000*10.0)*(1/av_number)*(1/VOLUME)*1e9
saturation_translation = (100000*1.0)*(1/av_number)*(1/VOLUME)*1e9

mRNA_half_life_TF = 9             # hrs
protein_half_life = 46            # hrs

const global degradation_constant_mRNA = -(1/mRNA_half_life_TF)*log(0.5)
const global degradation_constant_protein = -(1/protein_half_life)*log(0.5)

const global const_transcription = kcat_transcription/saturation_transcription
const global const_translation   = kcat_translation/saturation_translation

# ---------------------------- DO NOT EDIT BELOW THIS LINE -------------------------- #
data_dictionary = Dict();
data_dictionary["STOICHIOMETRIC_MATRIX"] = S;
data_dictionary["RATE_CONSTANT_ARRAY"] = rate_constant_array;
data_dictionary["SATURATION_CONSTANT_ARRAY"] = saturation_constant_array;
data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
data_dictionary["TIME_CONSTANT_ARRAY"] = time_constant_array;
data_dictionary["CONTROL_PARAMETER_ARRAY"] = control_parameter_array;
data_dictionary["MAXIMUM_SPECIFIC_GROWTH_RATE"] = maximum_specific_growth_rate;
data_dictionary["DILUTION_SELECTION_ARRAY"] = dilution_selection_array;
# ----------------------------------------------------------------------------------- #
return data_dictionary;
end


text:
####################################################
"""
    gmttrack(track::VisClaw.Track; kwargs...)
    gmttrack!(track::VisClaw.Track; kwargs...)
"""
gmttrack(track::VisClaw.Track, index=1:length(track.lon); kwargs...) =
GMT.plot(track.lon[index], track.lat[index]; kwargs...)
####################################################
"""
$(@doc gmttrack)
"""
gmttrack!(track::VisClaw.Track, index=1:length(track.lon); kwargs...) =
GMT.plot!(track.lon[index], track.lat[index]; kwargs...)
####################################################

ext: jl


text:
using ENN.TimedAutomata, Test

@testset "All tests" begin
@testset "Bounds tests" begin
    for T in (Rational{Int},Int,Float64)
        @testset "Time-type $(T)" begin
            b1 = TABound(true,typemin(T))
            b2 = TABound(true,-one(T))
            b3 = TABound(false,-one(T))
            b4 = TABound(true,zero(T))
            b5 = TABound(true,one(T))
            b6 = TABound(false,one(T))
            b7 = TABound(false,typemax(T))

            @test b1≤b2≤b3≤b4≤b5≤b6≤b7
            @test b1<b2<b3<b4<b5<b6<b7
            @test !(b6<b5)

            @test isinf.([b1,b2,b3,b4,b5,b6,b7]) == [true,false,false,false,false,false,true]

            b2.value = typemin(T)
            @test isinf(b2)

            @test b1==b2 && hash(b1)==hash(b2) && isequal(b1,b2)

            @test_throws ArgumentError b1+b7

            @test b1+b2 == typemin(TABound{T})
            @test b6+b7 == typemax(TABound{T})

            @test b3+b4 == TABound(true,-one(T))
            @test b3+b6 == TABound(false,zero(T))
            @test b5+b6 == TABound(true,one(T)+one(T))

            @test b3-b4 == TABound(true,-one(T))
            @test b3-b6 == TABound(false,-one(T)-one(T))
            @test b5-b6 == TABound(true,zero(T))

            @test repr.(Ref("text/plain"),[b1,b2,b3,b4,b5,b6,b7]) == ["-∞","-∞","≤ -$(one(T))","< $(zero(T))","< $(one(T))","≤ $(one(T))","∞"]
        end
    end
end

@testset "Constraint tests" begin
    c1 = TADiagonalConstraint(:x, :y, TABound(true, 5))
    c2 = TADiagonalConstraint(:y, :x, TABound(false, -5))
    c3 = TADiagonalConstraint(nothing, :x, TABound(false, -2))
    c4 = TADiagonalConstraint(:x, :y, TABound(false, -3))
    c5 = TADiagonalConstraint(:x, :y, TABound(false, 7))
    c6 = TADiagonalConstraint(:x, nothing, TABound(false, 5))

    @test repr.(Ref("text/plain"),[c1,c2,c3]) == ["x-y < 5","y-x ≤ -5","x ≥ 2"]
    @test c1 != c2
    @test (c1 == !c2) && isequal(c1,!c2) && hash(c1) == hash(!c2)

    tac = TAConstraint([c2,c3],[:x,:y])
    @test repr("text/plain", tac) == "2 ≤ x\ny-x ≤ -5"

    @test !isempty(tac)
    @test TimedAutomata.close(tac).D[1,2].value==-5
    @test isempty(tac ∩ c4)
    @test repr("text/plain", tac ∩ c4) == "∅"
    @test !isempty(tac ∩ c5 ∩ c6)
    @test repr("text/plain", tac ∩ c5 ∩ c6) == "x == 5\ny == 0\ny-x == -5"

    @test tac ∩ c5 ∩ c6 == tac ∩ TAConstraint([c5,c6],[:x,:y])
    @test tac ∩ c5 ∩ c6 ⊆ tac ∩ c5 ⊆ tac
    @test !(tac ∩ c5 ⊆ tac ∩ c6)

    tac2 = copy(tac)
    set_clock_list!(tac2, [:y,:x,:z])
    @test tac2.clocks == [:y,:x,:z] && tac2.D[[3,2],[3,2]] == tac.D[[2,3],[2,3]]

    cs = get_diagonal_constraints(tac)
    @test [c2,c3] ⊆ cs
    tac3 = TAConstraint(cs, [:x,:y])
    @test tac3 == tac && tac ⊆ tac3 && tac3 ⊆ tac

    tac4 = @constraint (y-x ≤ -5 && y ≥ 0 && 2 ≤ x) Int64
    set_clock_list!(tac4, [:x,:y])
    @test tac4 == tac3
end

@testset "Automata tests" begin
    states = ["S0", "S1", "S2", "S3"]
    initial_state = "S0"
    clocks = [:x,:y,:z]
    symbols = Set(Symbol[])
    arcs = Vector([
        @arc("S0","S1",true,TAMessage(),[:z],Int64),
        @arc("S1","S2",y>2,TAMessage(),[:y],Int64),
        @arc("S2","S3",x-z<1 && z-y<1,TAMessage(),[],Int64)
    ])
    invariants = Dict("S0"=>@constraint x-y==0 && y-z==0 && z-x==0 Int64)

    automaton1 = TA(
        states,
        initial_state,
        clocks,
        symbols,
        arcs,
        invariants
    )

    zg1 = TTS(automaton1)

    s1=TimedAutomata.TAState(3,4,5,6,7)
    s2=TimedAutomata.TAState(8)
    s3=TimedAutomata.TAState(9,10)

    s4 = s1|s2|s3
    s5 = |(s1,s2,s3)
    @test s4==s5
    @test collect(s4) == collect(s5) == collect(s4[1:end]) == [3, 4, 5, 6, 7, 8, 9, 10]

    automaton4=TA(
        ["off","dim","bright"],
        "off",
        [:x],
        Set([:press]),
        Vector([
            @arc("off", "dim", true, TAMessage(MSG_IN, :press), [:x]),
            @arc("dim", "off", x>10//1, TAMessage(MSG_IN, :press), []),
            @arc("dim", "bright", x≤10//1, TAMessage(MSG_IN, :press), []),
            @arc("bright", "off", true, TAMessage(MSG_IN, :press), [])
        ]),
        Dict{String,TAConstraint{Rational{Int}}}()
    )
    zg4 = TTS(automaton4)

    automaton5=TA(
        ["idle","t","study", "relax"],
        "idle",
        [:y],
        Set([:press]),
        Vector([
            @arc("t", "study", true, TAMessage(MSG_OUT, :press), []),
            @arc("study", "study", true, TAMessage(), []),
            @arc("study", "idle", true, TAMessage(MSG_OUT, :press), []),
            @arc("idle", "t", true, TAMessage(MSG_OUT, :press), [:y]),
            @arc("idle", "relax", true, TAMessage(MSG_OUT, :press), [:y]),
            @arc("relax", "idle", y>10//1, TAMessage(MSG_OUT, :press), [])
        ]),
        Dict(
            "t" => @constraint y<5//1
        )
    )

    automaton4p = prune(automaton4)
    automaton44 = (automaton4 | automaton4)
    automaton45 = automaton4 | automaton5
    automaton54 = automaton5 | automaton4

    automaton44p = prune(automaton44)
    automaton45p = prune(automaton45)

    @test !isnothing(find_isomorphism(automaton45,automaton54))
    @test isnothing(find_isomorphism(automaton44,automaton45))
    @test isnothing(find_isomorphism(automaton44, automaton4))

    automaton44p′ = deepcopy(automaton44p)
    TimedAutomata.set_clock_list!(automaton44p′, automaton4p.clocks; mapping=(_)->:x)
    @test !isnothing(find_isomorphism(automaton44p′, automaton4p))
    automaton44p_5 = automaton44p′ | automaton5
    @test !isnothing(find_isomorphism(automaton44p_5, automaton45))

    @test length(automaton45.states)== length(automaton4.states)*length(automaton5.states)
    @test length(automaton45p.states)== 4

    automaton445 = |(automaton4, automaton4, automaton5)
    @test length(automaton445.states) == length(automaton4.states)^2 * length(automaton5.states)

    automaton544 = |(automaton5, automaton4, automaton4)
    automaton44_5 = automaton44 | automaton5
    automaton4_45 = automaton4 | (automaton4 | automaton5)
    automaton5_44 = automaton5 | automaton44

    automata = (automaton445, automaton544, automaton44_5, automaton4_45, automaton5_44)
    for a1 in automata
        a1p=prune(a1)
        @test length(a1p.states) == length(automaton45p.states)

        for a2 in automata
            a2p=prune(a2)

            @test !isnothing(find_isomorphism(a1,a2))
            @test !isnothing(find_isomorphism(a1p,a2p))
        end
    end
end

@testset "Regex & language tests" begin
    @test_broken false # must still implement tests!
end

end


text:
"""
```
struct VarName{sym, T<:Tuple}
    indexing::T
end
```

A variable identifier. Every variable has a symbol `sym` and indices `indexing` in the format
returned by [`@vinds`](@ref).  The Julia variable in the model corresponding to `sym` can refer to a
single value or to a hierarchical array structure of univariate, multivariate or matrix
variables. `indexing` stores the indices that can access the random variable from the Julia
variable.

Examples:

- `x[1] ~ Normal()` will generate a `VarName` with `sym == :x` and `indexing == "((1,))"`.
- `x[:,1] ~ MvNormal(zeros(2))` will generate a `VarName` with `sym == :x` and
 `indexing == ((Colon(), 1))"`.
- `x[:,1][2] ~ Normal()` will generate a `VarName` with `sym == :x` and
 `indexing == ((Colon(), 1), (2,))`.
"""
struct VarName{sym, T<:Tuple}
    indexing::T
end

VarName(sym::Symbol, indexing::Tuple = ()) = VarName{sym, typeof(indexing)}(indexing)

"""
    VarName(vn::VarName, indexing)

Return a copy of `vn` with a new index `indexing`.
"""
function VarName(vn::VarName, indexing::Tuple = ())
    return VarName{getsym(vn), typeof(indexing)}(indexing)
end

"""
    getsym(vn::VarName)

Return the symbol of the Julia variable used to generate `vn`.
"""
getsym(vn::VarName{sym}) where sym = sym

"""
    getindexing(vn::VarName)

Return the indexing tuple of the Julia variable used to generate `vn`.
"""
getindexing(vn::VarName) = vn.indexing

Base.hash(vn::VarName, h::UInt) = hash((getsym(vn), getindexing(vn)), h)
Base.:(==)(x::VarName, y::VarName) = getsym(x) == getsym(y) && getindexing(x) == getindexing(y)

function Base.show(io::IO, vn::VarName)
    print(io, getsym(vn))
    for indices in getindexing(vn)
        print(io, "[")
        join(io, indices, ",")
        print(io, "]")
    end
end

"""
    Symbol(vn::VarName)

Return a `Symbol` represenation of the variable identifier `VarName`.
"""
Base.Symbol(vn::VarName) = Symbol(string(vn))  # simplified symbol

"""
    inspace(vn::Union{VarName, Symbol}, space::Tuple)

Check whether `vn`'s variable symbol is in `space`.
"""
inspace(vn, space::Tuple{}) = true # empty space is treated as universal set
inspace(vn, space::Tuple) = vn in space
inspace(vn::VarName, space::Tuple{}) = true
inspace(vn::VarName, space::Tuple) = any(_in(vn, s) for s in space)

@noinline function Base.in(vn::VarName, space::Tuple)
    Base.depwarn("`Base.in(vn::VarName, space::Tuple)` is deprecated, use `inspace(vn, space)` instead.",
                 nameof(Base.in))
    return inspace(vn, space)
end

_in(vn::VarName, s::Symbol) = getsym(vn) == s
_in(vn::VarName, s::VarName) = subsumes(s, vn)

"""
    subsumes(u::VarName, v::VarName)

Check whether the variable name `v` describes a sub-range of the variable `u`.  Supported
indexing:

- Scalar: `x` subsumes `x[1, 2]`, `x[1, 2]` subsumes `x[1, 2][3]`, etc.
- Array of scalar: `x[[1, 2], 3]` subsumes `x[1, 3]`, `x[1:3]` subsumes `x[2][1]`, etc.
  (basically everything that fulfills `issubset`).
- Slices: `x[2, :]` subsumes `x[2, 10][1]`, etc.

Currently _not_ supported are:

- Boolean indexing, literal `CartesianIndex` (these could be added, though)
- Linear indexing of multidimensional arrays: `x[4]` does not subsume `x[2, 2]` for `x` a matrix
- Trailing ones: `x[2, 1]` does not subsume `x[2]` for `x` a vector
"""
function subsumes(u::VarName, v::VarName)
    return getsym(u) == getsym(v) && subsumes(u.indexing, v.indexing)
end

subsumes(::Tuple{}, ::Tuple{}) = true  # x subsumes x
subsumes(::Tuple{}, ::Tuple) = true    # x subsumes x[1]
subsumes(::Tuple, ::Tuple{}) = false   # x[1] does not subsume x
function subsumes(t::Tuple, u::Tuple)  # does x[i]... subsume x[j]...?
    return _issubindex(first(t), first(u)) && subsumes(Base.tail(t), Base.tail(u))
end

const AnyIndex = Union{Int, AbstractVector{Int}, Colon}
_issubindex_(::Tuple{Vararg{AnyIndex}}, ::Tuple{Vararg{AnyIndex}}) = false
function _issubindex(t::NTuple{N, AnyIndex}, u::NTuple{N, AnyIndex}) where {N}
    return all(_issubrange(j, i) for (i, j) in zip(t, u))
end

const ConcreteIndex = Union{Int, AbstractVector{Int}} # this include all kinds of ranges
"""Determine whether indices `i` are contained in `j`, treating `:` as universal set."""
_issubrange(i::ConcreteIndex, j::ConcreteIndex) = issubset(i, j)
_issubrange(i::Union{ConcreteIndex, Colon}, j::Colon) = true
_issubrange(i::Colon, j::ConcreteIndex) = true

"""
    @varname(expr)

A macro that returns an instance of `VarName` given the symbol or expression of a Julia variable,
e.g. `@varname x[1,2][1+5][45][3]` returns `VarName{:x}(((1, 2), (6,), (45,), (3,)))`.

!!! compat "Julia 1.5"
    Using `begin` in an indexing expression to refer to the first index requires at least
    Julia 1.5.
"""
macro varname(expr::Union{Expr, Symbol})
    return esc(varname(expr))
end

varname(expr::Symbol) = VarName(expr)
function varname(expr::Expr)
    if Meta.isexpr(expr, :ref)
        sym, inds = vsym(expr), vinds(expr)
        return :($(DynamicPPL.VarName)($(QuoteNode(sym)), $inds))
    else
        throw("VarName: Mis-formed variable name $(expr)!")
    end
end

"""
    @vsym(expr)

A macro that returns the variable symbol given the input variable expression `expr`.
For example, `@vsym x[1]` returns `:x`.
"""
macro vsym(expr::Union{Expr, Symbol})
    return QuoteNode(vsym(expr))
end

vsym(expr::Symbol) = expr
function vsym(expr::Expr)
    if Meta.isexpr(expr, :ref)
        return vsym(expr.args[1])
    else
        throw("VarName: Mis-formed variable name $(expr)!")
    end
end

"""
    @vinds(expr)

Returns a tuple of tuples of the indices in `expr`. For example, `@vinds x[1, :][2]` returns
`((1, Colon()), (2,))`.

!!! compat "Julia 1.5"
    Using `begin` in an indexing expression to refer to the first index requires at least
    Julia 1.5.
"""
macro vinds(expr::Union{Expr, Symbol})
    return esc(vinds(expr))
end

vinds(expr::Symbol) = Expr(:tuple)
function vinds(expr::Expr)
    if Meta.isexpr(expr, :ref)
        ex = copy(expr)
        @static if VERSION < v"1.5.0-DEV.666"
            Base.replace_ref_end!(ex)
        else
            Base.replace_ref_begin_end!(ex)
        end
        last = Expr(:tuple, ex.args[2:end]...)
        init = vinds(ex.args[1]).args
        return Expr(:tuple, init..., last)
    else
        throw("VarName: Mis-formed variable name $(expr)!")
    end
end

@generated function inargnames(::VarName{s}, ::Model{_F, argnames}) where {s, argnames, _F}
    return s in argnames
end

@generated function inmissings(::VarName{s}, ::Model{_F, _a, _T, missings}) where {s, missings, _F, _a, _T}
    return s in missings
end


text:
# Log sum exp.
function lse(arr)
    max = maximum(arr)
    max == -Inf ? -Inf : max + log(sum(exp.(arr .- max)))
end

function lse(x1::Real, x2::Real)
    m = max(x1, x2)
    m == -Inf ? m : m + log(exp(x1 - m) + exp(x2 - m))
end

# Effective sample size.
function ess(lnw::Vector{Float64})
    log_ess = -lse(2. * lnw)
    return exp(log_ess)
end

# Normalize log weights.
function nw(lw::Vector{Float64})
    lt = lse(lw)
    lnw = lw .- lt
    return (lt, lnw)
end


text:
include("settings.jl")

makedocs(; sitename=sitename, settings...)

deploydocs(;
  repo="github.com/ITensor/ITensors.jl.git",
  devbranch="main",
  push_preview=true,
  deploy_config=Documenter.GitHubActions(),
)


text:
function build_gram_matrix(q::Vector,
                           monos::AbstractVector{<:MP.AbstractMonomial})
    return GramMatrix(MultivariateMoments.SymMatrix(q, length(monos)),
                         monos)
end

function build_moment_matrix(q::Vector,
                             monos::AbstractVector{<:MP.AbstractMonomial})
    return MomentMatrix(MultivariateMoments.SymMatrix(q, length(monos)),
                        monos)
end

abstract type SOSLikeCone <: PolyJuMP.PolynomialSet end
Base.broadcastable(cone::SOSLikeCone) = Ref(cone)

struct SOSPolynomialSet{DT <: AbstractSemialgebraicSet,
                        CT <: SOSLikeCone,
                        BT <: PolyJuMP.AbstractPolynomialBasis,
                        MT <: MP.AbstractMonomial,
                        MVT <: AbstractVector{MT},
                        NPT <: Tuple} <: MOI.AbstractVectorSet
    domain::DT
    cone::CT
    basis::Type{BT}
    monomials::MVT
    newton_polytope::NPT
    mindegree::Int
    maxdegree::Int
end


text:
#=##############################################################################
# DESCRIPTION
    Abstract body type definition.
# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Jun 2018
  * License   : MIT License
=###############################################################################

################################################################################
# ABSTRACT BODY TYPE
################################################################################
"""
  Implementations of AbstractBody are expected to have the following fields.
  * `grid::gt.GridTriangleSurface `     : Paneled geometry
  * `nnodes::Int64`                     : Number of nodes
  * `ncells::Int64`                     : Number of cells
  * `fields::Array{String, 1}`          : Available fields (solutions)
  * `Oaxis::Array{T,2} where {T<:Real}` : Coordinate system of original grid
  * `O::Array{T,1} where {T<:Real}`     : Position of CS of original grid

  and the following functions

  ```julia

  # Solve for distributions of potential flow
  function solve(self::BodyType, Vinfs::Array{Array{T,1},1}, args...
) where {T<:Real}
    .
    .
    .
  end

  # Returns the velocity induced by the body on the targets `targets`. It adds
  # the velocity at the i-th target to out[i].
  function _Vind(self::BodyTypes, targets::Array{Array{T,1},1},
out::Array{Array{T,1},1}, args...; optargs...) where{T<:Real}
    .
    .
    .
  end
  ```
"""
abstract type AbstractBody end

# Includes all implementations of AbstractBody
for header_name in ["nonliftingbody", "abstractliftingbody"]
  include("FLOWPanel_"*header_name*".jl")
end

# Declares implementations of AbstractGrid
const BodyTypes = Union{NonLiftingBody, NonLiftingBodyDoublet,
                                                NonLiftingBodyVRing, LBodyTypes}

##### COMMON FUNCTIONS  ########################################################

"""
  `Vind(self::BodyTypes, targets::Array{Array{T,1},1},
                  out::Array{Array{T,1},1}, args...; optargs...) where{T<:Real}`

Returns the velocity induced by the body on the targets `targets`. It adds the
velocity at the i-th target to out[i].
"""
function Vind(self::BodyTypes, targets::Array{Array{T1,1},1},
                  out::Array{Array{T2,1},1}, args...; optargs...
                                                  ) where{T1<:RType, T2<:RType}

  # ERROR CASES
  if check_solved(self)==false
    error("Body hasn't been solved yet. Please call `solve()` function first.")
  end

  _Vind(self, targets, out, args...; optargs...)
end

"""
  `save(body::BodyTypes, filename::String; opt_args...)`

Outputs a vtk file of this body. See GeometricTools.save(grid, ...) for a
description of optional arguments `opt_args...`.
"""
function save(body::BodyTypes, filename::String; out_cellindex::Bool=false,
                                                 out_cellindexdim::Array{Int64,1}=Int64[],
                                                 out_nodeindex::Bool=false,
                                                 out_controlpoints::Bool=false,
                                                 out_wake::Bool=true,
                                                 _len::RType=1.0,
                                                 _upper::Bool=true,
                                                 opt_args...)

  str = ""

  # Add special fields
  if out_cellindex
    gt.add_field(body.grid, "cellindex", "scalar",
                    [i for i in 1:body.ncells], "cell")
  end
  if out_nodeindex
    gt.add_field(body.grid, "nodeindex", "scalar",
                    [i for i in 1:body.nnodes], "node")
  end
  for dim in out_cellindexdim
    ndivs = gt.get_ndivscells(body.grid)[1:2]
    data = [ Base._ind2sub(ndivs, i)[dim] for i in 1:body.ncells]
    gt.add_field(body.grid, "cellindexdim$(dim)", "scalar", data, "cell")
  end

  # Outputs control points
  if out_controlpoints
    str *= save_controlpoints(body, filename; opt_args...)
  end

  # Outputs wake
  if out_wake
    # Case that body is not a LiftingBody
    try
       body::LBodyTypes
       if body.nnodesTE-1 != 0
           str *= _savewake(body, filename; len=_len, upper=_upper, opt_args...)
       end
     catch e
       if isa(e, TypeError)
         nothing
       else
         rethrow(e)
       end
     end
  end

  # Saves body
  return str*gt.save(body.grid, filename; opt_args...)

end

"""
  `save_controlpoints(body::BodyTypes, filename::String;
suffix::String="_cp", opt_args...)`

Outputs a vtk file with the control points of the body along with associated
normals and surface velocities.
"""
function save_controlpoints(body::BodyTypes, filename::String;
                                              suffix::String="_cp", opt_args...)
    # Control points
    CPs = [get_controlpoint(body, i) for i in 1:body.ncells]

    # Normals
    data = [
            Dict( "field_name"  => "normal",
                  "field_type"  => "vector",
                  "field_data"  => [get_normal(body, i) for i in 1:size(CPs,1)])
            ]

    # Surface velocity
    if check_solved(body)
      Vsurf = [get_fieldval(body, "Vinf", i) for i in 1:size(CPs,1)]
      Vind(body, CPs, Vsurf)

      push!(data,
              Dict( "field_name"  => "V",
                    "field_type"  => "vector",
                    "field_data"  => Vsurf)
      )
    end

    # Generates vtk
    return gt.generateVTK(filename*suffix, CPs; point_data=data, opt_args...)
end

"""
  `get_controlpoint(body::BodyTypes, i::Int64 or coor::Array{Int64,1})`

Returns the control point on the i-th panel.
"""
function get_controlpoint(body::BodyTypes, args...)
  return _get_controlpoint(body.grid, args...)
end

"""
  `get_ndivscells(body::BodyTypes)`

Returns a tuple with the number of cells in each parametric dimension
"""
get_ndivscells(body::BodyTypes) = body.grid._ndivscells

"""
  `get_ndivsnodes(body::BodyTypes)`

Returns a tuple with the number of nodes in each parametric dimension
"""
get_ndivsnodes(body::BodyTypes) = body.grid._ndivsnodes

"""
  `get_unitvectors(body::BodyTypes, i::Int64 or coor::Array{Int64,1})`

Returns the unit vectors t,n,o of the i-th panel, with t the tanget vector,
n normal, and o oblique.
"""
get_unitvectors(body::BodyTypes, args...) = gt.get_unitvectors(body.grid, args...)

"""
  `get_normal(body::BodyTypes, i::Int64 or coor::Array{Int64,1})`

Returns the normal vector the i-th panel.
"""
get_normal(body::BodyTypes, args...) = gt.get_normal(body.grid, args...)

"""
  `get_field(self::BodyTypes, field_name::String)`

Returns the requested field.
"""
function get_field(self::BodyTypes, field_name::String)
  if !(field_name in self.fields)
    error("Field $field_name not found! Available fields: $(self.fields)")
  end

  return self.grid.field[field_name]
end

"""
  `get_field(self::BodyTypes, field_name::String, i::Int64)`

Returns the requested field value. Give it `_check=false` to skip checking logic
for faster operation.
"""
function get_fieldval(self::BodyTypes, field_name::String, i::Int64;
                      _check::Bool=true)
  if _check
    if i<1
      error("Invalid index $i.")
    elseif FIELDS[field_name]["entry_type"]=="node" && i>self.nnodes
      error("Invalid index $i. Maximum is $(self.nnodes).")
    elseif FIELDS[field_name]["entry_type"]=="cell" && i>self.ncells
      error("Invalid index $i. Maximum is $(self.ncells).")
    end
  end

  return self.grid.field[field_name]["field_data"][i]
end

"""
  `add_field(self::BodyTypes, field_name::String, field_data)`

Adds a new field to the body. It overwrites the field if it already existed.
"""
function add_field(self::BodyTypes, field_name::String, field_data)

  # ERROR CASES
  if !(field_name in keys(FIELDS))
    error("Invalid field $field_name. Implemented fields are: $(keys(FIELDS))")
  end

  # Adds field to grid
  gt.add_field(self.grid, field_name, FIELDS[field_name]["field_type"],
                field_data, FIELDS[field_name]["entry_type"]; raise_warn=false)

  # Registers the field
  if !(field_name in self.fields)
    push!(self.fields, field_name)
  end

  nothing
end

"""
  `check_field(self::BodyTypes, field_name::String)`

Returns `true` of the body has the field `field_name`. Returns false otherwise.
"""
check_field(self::BodyTypes, field_name::String) = field_name in self.fields

"""
  `check_solved(self::BodyTypes)`

Returns `true` of the body has been solved. Returns false otherwise.
"""
function check_solved(self::BodyTypes)
  if check_field(self, "solved")
    return get_fieldval(self, "solved", 1)
  else
    return false
  end
end

"""
  `rotate(body::BodyTypes, roll::Real, pitch::Real, yaw::Real;
translation::Array{T, 1}=zeros(3), reset_fields::Bool=true
) where{T<:Real}`

Rotates and translates the body by the given axial angles.

NOTE: Naming follows aircraft convention, with
* roll:   rotation about x-axis.
* pitch:  rotation about y-axis.
* yaw:    rotation about z-axis.
"""
function rotate(body::BodyTypes, roll::RType, pitch::RType, yaw::RType;
                  translation::Array{T, 1}=zeros(3), reset_fields::Bool=true
                ) where{T<:RType}

  M = gt.rotation_matrix2(roll, pitch, yaw)
  gt.lintransform!(body.grid, M, translation; reset_fields=reset_fields)

  body.Oaxis[:,:] = M*body.Oaxis
  body.O[:] += translation

  nothing
end

##### COMMON INTERNAL FUNCTIONS  ###############################################
function _get_controlpoint(grid::gt.GridTriangleSurface, args...)
  return _get_controlpoint(gt.get_cellnodes(grid, args...),
                                                  gt.get_normal(grid, args...))
end

function _get_controlpoint(cellnodes::Array{Array{T1,1},1}, normal::Array{T2,1},
                                                  ) where{T1<:RType, T2<:RType}
  len = maximum([norm(cellnodes[1]-v) for v in cellnodes[2:end]])
  return mean(cellnodes) + 0.005*len*normal
end

function _solvedflag(self::BodyTypes, val::Bool)
  add_field(self, "solved", [val])
end

##### END OF ABSTRACT BODY #####################################################


text:
using Distributions
using Test

# Test for heslogpdf on univariate distributions

@test isapprox(heslogpdf(Beta(1.5, 3.0), 0.7)   , -23.242630385487523 , atol=1.0e-8)
@test isapprox(heslogpdf(Chi(5.0), 5.5)         , -1.1322314049586777 , atol=1.0e-8)
@test isapprox(heslogpdf(Chisq(7.0), 12.0)      , -0.01736111111111111, atol=1.0e-8)
@test isapprox(heslogpdf(Exponential(2.0), 7.0) ,  0.0                , atol=1.0e-8)
@test isapprox(heslogpdf(Frechet(2.0,4.0), 9.0) ,  0.02240512117055326, atol=1.0e-8)
@test isapprox(heslogpdf(Gamma(9.0, 0.5), 11.0) , -0.06611570247933884, atol=1.0e-8)
@test isapprox(heslogpdf(Gumbel(3.5, 1.0), 4.0) ,  0.6065306597126334 , atol=1.0e-8)
@test isapprox(heslogpdf(Laplace(7.0), 34.0)    ,  0.0                , atol=1.0e-8)
@test isapprox(heslogpdf(Logistic(-6.0), 1.0)   , -0.00182044236024365, atol=1.0e-8)
@test isapprox(heslogpdf(LogitNormal(0.8), 0.4) , -0.12837092681258455, atol=1.0e-8)
@test isapprox(heslogpdf(LogNormal(5.5), 2.0)   , -1.2017132048600137 , atol=1.0e-8)
@test isapprox(heslogpdf(Normal(-4.5, 2.0), 1.6), -0.25               , atol=1.0e-8)
@test isapprox(heslogpdf(TDist(8.0), 9.1)       ,  0.0816459812355031 , atol=1.0e-8)
@test isapprox(heslogpdf(Weibull(2.0), 3.5)     , -2.0816326530612246 , atol=1.0e-8)
@test isapprox(heslogpdf(Weibull(2.0), -3.5)    ,  0.0                , atol=1.0e-8)


text:
module TestCore

using Atomix: AtomicRefArray
using Test
using UnsafeAtomics: UnsafeAtomics

function test()
    A = AtomicRefArray(ones(Int, 3))
    @test UnsafeAtomics.load(A[1]) === 1
    UnsafeAtomics.store!(A[1], 123)
    @test UnsafeAtomics.load(A[1]) === 123
    @test UnsafeAtomics.modify!(A[1], +, 1) === (123, 124)
    @test UnsafeAtomics.load(A[1]) === 124
    @test UnsafeAtomics.cas!(A[1], 124, 567) === (old = 124, success = true)
    @test UnsafeAtomics.cas!(A[1], 124, 567) === (old = 567, success = false)
end

end  # module


text:
using ThreadedSparseCSR
using SparseArrays, SparseMatricesCSR
using Test

@testset "Base.Threads matrix-vec multipy" begin
    m, n = 1_000, 1_000
    d = 0.01
    num_nzs = floor(Int, m*n*d)
    rows = rand(1:m, num_nzs)
    cols = rand(1:n, num_nzs)
    vals = rand(num_nzs)

    cscA = sparse(rows, cols, vals, m, n)
    csrA = sparsecsr(rows, cols, vals, m, n)
    x = rand(n)

    @testset "5-arg, inplace mat-vec mul" begin
        y1 = rand(n)
        y2 = copy(y1)
        @test tmul!(y1, csrA, x, 0.5, 0.3) ≈ SparseArrays.mul!(y2, cscA, x, 0.5, 0.3)

        y1 = rand(n)
        y2 = copy(y1)
        @test tmul!(y1, csrA, x, 0.5, 0.3) ≈ SparseMatricesCSR.mul!(y2, csrA, x, 0.5, 0.3)
    end

    @testset "3-arg, inplace mat-vec mul" begin
        y1 = rand(n)
        y2 = copy(y1)
        @test tmul!(y1, csrA, x) ≈ SparseArrays.mul!(y2, cscA, x)

        y1 = rand(n)
        y2 = copy(y1)
        @test tmul!(y1, csrA, x) ≈ SparseMatricesCSR.mul!(y2, csrA, x)
    end

   @testset "2-arg mat-vec mul" begin
        @test tmul(csrA, x) ≈ cscA*x
        @test tmul(csrA, x) ≈ SparseMatricesCSR.:*(csrA, x)
    end

end

@testset "Polyester matrix-vec multipy" begin
    m, n = 1_000, 1_000
    d = 0.01
    num_nzs = floor(Int, m*n*d)
    rows = rand(1:m, num_nzs)
    cols = rand(1:n, num_nzs)
    vals = rand(num_nzs)

    cscA = sparse(rows, cols, vals, m, n)
    csrA = sparsecsr(rows, cols, vals, m, n)
    x = rand(n)

    @testset "5-arg, inplace mat-vec mul" begin
        y1 = rand(n)
        y2 = copy(y1)
        @test bmul!(y1, csrA, x, 0.5, 0.3) ≈ SparseArrays.mul!(y2, cscA, x, 0.5, 0.3)

        y1 = rand(n)
        y2 = copy(y1)
        @test bmul!(y1, csrA, x, 0.5, 0.3) ≈ SparseMatricesCSR.mul!(y2, csrA, x, 0.5, 0.3)
    end

    @testset "3-arg, inplace mat-vec mul" begin
        y1 = rand(n)
        y2 = copy(y1)
        @test bmul!(y1, csrA, x) ≈ SparseArrays.mul!(y2, cscA, x)

        y1 = rand(n)
        y2 = copy(y1)
        @test bmul!(y1, csrA, x) ≈ SparseMatricesCSR.mul!(y2, csrA, x)
    end

   @testset "2-arg mat-vec mul" begin
        @test bmul(csrA, x) ≈ cscA*x
        @test bmul(csrA, x) ≈ SparseMatricesCSR.:*(csrA, x)
    end

end


text:
doc"""
    Frechet(α,θ)

The *Fréchet distribution* with shape `α` and scale `θ` has probability density function

$f(x; \alpha, \theta) = \frac{\alpha}{\theta} \left( \frac{x}{\theta} \right)^{-\alpha-1}
e^{-(x/\theta)^{-\alpha}}, \quad x > 0$

```julia
Frechet()        # Fréchet distribution with unit shape and unit scale, i.e. Frechet(1.0, 1.0)
Frechet(a)       # Fréchet distribution with shape a and unit scale, i.e. Frechet(a, 1.0)
Frechet(a, b)    # Fréchet distribution with shape a and scale b

params(d)        # Get the parameters, i.e. (a, b)
shape(d)         # Get the shape parameter, i.e. a
scale(d)         # Get the scale parameter, i.e. b
```

External links

* [Fréchet_distribution on Wikipedia](http://en.wikipedia.org/wiki/Fréchet_distribution)

"""
immutable Frechet <: ContinuousUnivariateDistribution
    α::Float64
    θ::Float64

    function Frechet(α::Real, θ::Real)
    	@check_args(Frechet, α > zero(α) && θ > zero(θ))
    	new(α, θ)
    end
    Frechet(α::Real) = Frechet(α, 1.0)
    Frechet() = new(1.0, 1.0)
end

@distr_support Frechet 0.0 Inf

#### Parameters

shape(d::Frechet) = d.α
scale(d::Frechet) = d.θ
params(d::Frechet) = (d.α, d.θ)

#### Statistics

mean(d::Frechet) = (α = d.α; α > 1.0 ? d.θ * gamma(1.0 - 1.0 / α) : Inf)

median(d::Frechet) = d.θ * logtwo^(-1.0 / d.α)

mode(d::Frechet) = (iα = -1.0/d.α; d.θ * (1.0 - iα) ^ iα)

function var(d::Frechet)
    if d.α > 2.0
        iα = 1.0 / d.α
        return d.θ^2 * (gamma(1.0 - 2.0 * iα) - gamma(1.0 - iα)^2)
    else
        return Inf
    end
end

function skewness(d::Frechet)
    if d.α > 3.0
        iα = 1.0 / d.α
        g1 = gamma(1.0 - iα)
        g2 = gamma(1.0 - 2.0 * iα)
        g3 = gamma(1.0 - 3.0 * iα)
        return (g3 - 3.0 * g2 * g1 + 2 * g1^3) / ((g2 - g1^2)^1.5)
    else
        return Inf
    end
end

function kurtosis(d::Frechet)
    if d.α > 3.0
        iα = 1.0 / d.α
        g1 = gamma(1.0 - iα)
        g2 = gamma(1.0 - 2.0 * iα)
        g3 = gamma(1.0 - 3.0 * iα)
        g4 = gamma(1.0 - 4.0 * iα)
        return (g4 - 4.0 * g3 * g1 + 3 * g2^2) / ((g2 - g1^2)^2) - 6.0
    else
        return Inf
    end
end

function entropy(d::Frechet)
    const γ = 0.57721566490153286060  # γ is the Euler-Mascheroni constant
    1.0 + γ / d.α + γ + log(d.θ / d.α)
end

#### Evaluation

function logpdf(d::Frechet, x::Float64)
    (α, θ) = params(d)
    if x > 0.0
        z = θ / x
        return log(α / θ) + (1.0 + α) * log(z) - z^α
    else
        return -Inf
    end
end

pdf(d::Frechet, x::Float64) = exp(logpdf(d, x))

cdf(d::Frechet, x::Float64) = x > 0.0 ? exp(-((d.θ / x) ^ d.α)) : 0.0
ccdf(d::Frechet, x::Float64) = x > 0.0 ? -expm1(-((d.θ / x) ^ d.α)) : 1.0
logcdf(d::Frechet, x::Float64) = x > 0.0 ? -(d.θ / x) ^ d.α : -Inf
logccdf(d::Frechet, x::Float64) = x > 0.0 ? log1mexp(-((d.θ / x) ^ d.α)) : 0.0

quantile(d::Frechet, p::Float64) = d.θ * (-log(p)) ^ (-1.0 / d.α)
cquantile(d::Frechet, p::Float64) = d.θ * (-log1p(-p)) ^ (-1.0 / d.α)
invlogcdf(d::Frechet, lp::Float64) = d.θ * (-lp)^(-1.0 / d.α)
invlogccdf(d::Frechet, lp::Float64) = d.θ * (-log1mexp(lp))^(-1.0 / d.α)

function gradlogpdf(d::Frechet, x::Float64)
    (α, θ) = params(d)
    insupport(Frechet, x) ? -(α + 1.0) / x + α * (θ^α) * x^(-α-1.0)  : 0.0
end

## Sampling

rand(d::Frechet) = d.θ * randexp() ^ (-1.0 / d.α)


text:
#=
Define a transition model for the PedCarMDP problem
=#

function POMDPs.transition(mdp::PedCarMDP, s::PedCarMDPState, a::PedCarMDPAction)
    ego_ps, ego_probs = ego_transition(mdp, s, a)
    ped_ps, ped_probs = ped_transition(mdp, s, a)
    car_ps, car_routes, car_probs = car_transition(mdp, s, a)
    nps = length(car_ps)*length(ped_ps)*length(ego_ps)
    states_p = Vector{PedCarMDPState}(nps)
    states_probs = Vector{Float64}(nps)
    idx = 1
    for ie=1:length(ego_ps)
        for ip=1:length(ped_ps)
            for ic=1:length(car_ps)
                weight = ego_probs[ie]*ped_probs[ip]*car_probs[ic]
                # if !(weight ≈ 0.)
                    # collision=false
                    collision = crash(mdp, ego_ps[ie], car_ps[ic], ped_ps[ip])
                    # collision = mdp._collision_checker[(ego_ps[ie], car_ps[ic], ped_ps[ip])]
                    # collision = is_colliding(Vehicle(ego_ps[ie], mdp.ego_type, EGO_ID), Vehicle(car_ps[ic], mdp.car_type, CAR_ID)) || is_colliding(Vehicle(ego_ps[ie], mdp.ego_type, EGO_ID), Vehicle(ped_ps[ip], mdp.ped_type, PED_ID))
                    states_p[idx] = PedCarMDPState(collision, ego_ps[ie], ped_ps[ip], car_ps[ic], car_routes[ic])
                    states_probs[idx] = weight
                    idx += 1
                    # push!(states_p, PedCarMDPState(collision, ego_ps[ie], ped_ps[ip], car_ps[ic], car_routes[ic]))
                    # push!(states_probs, weight)
                # end
            end
        end
    end
    normalize!(states_probs, 1)
    return SparseCat(states_p, states_probs)
end

function ego_transition(mdp::PedCarMDP, s::PedCarMDPState, a::PedCarMDPAction)
    # ego transition
    acc = LonAccelDirection(a.acc, 4) #TODO paremeterize
    ego_p = propagate(s.ego, acc, mdp.env.roadway, mdp.ΔT)
    # interpolate ego_p in discrete space
    ego_ps, ego_probs = interpolate_state(mdp, ego_p)
    return ego_ps, ego_probs
end

function ped_transition(mdp::PedCarMDP, s::PedCarMDPState, a::PedCarMDPAction)
    # pedestrian transition
    if s.ped.posG == mdp.off_grid
        ped_ps = pedestrian_starting_states(mdp)
        push!(ped_ps, get_off_the_grid(mdp))
        ped_probs = ones(length(ped_ps))
        ped_probs[end] = 1-mdp.ped_birth
        ped_probs[1:end-1] = mdp.ped_birth/(length(ped_ps)-1)
        @assert sum(ped_probs) ≈ 1.0
    else
        ped_actions = (ConstantSpeedDawdling(0., 0.), ConstantSpeedDawdling(1., 0.), ConstantSpeedDawdling(2., 0.))
        ped_actions_probs = (1/3, 1/3, 1/3)
        ped_ps = VehicleState[]
        ped_probs = Float64[]
        for (i, ped_a) in enumerate(ped_actions)
            p_a = ped_actions_probs[i]
            ped_p = propagate(s.ped, ped_a, mdp.env.roadway, mdp.ΔT)
            if (ped_p.posF.s >= get_end(get_lane(mdp.env.roadway, ped_p)) && isapprox(ped_p.posF.ϕ, 0.)) ||
               (isapprox(ped_p.posF.s, 0., atol=0.01) && isapprox(ped_p.posF.ϕ, float(pi)))
                ped_pp = get_off_the_grid(mdp)
                index_ped_pp = findall(x -> x==ped_pp, ped_ps)
                if isempty(index_ped_pp)
                    push!(ped_ps, get_off_the_grid(mdp))
                    push!(ped_probs, p_a)
                else
                    ped_probs[index_ped_pp] += p_a
                end
                continue
            end
            # interpolate ped_p in discrete space
            itp_ped_ps, itp_ped_weights = interpolate_pedestrian(mdp, ped_p)
            for (j, ped_pss) in enumerate(itp_ped_ps)
                index_itp_state = findall(x -> x==ped_pss, ped_ps)
                if !(itp_ped_weights[j] ≈ 0.)
                    if isempty(index_itp_state)
                        push!(ped_ps, ped_pss)
                        push!(ped_probs, itp_ped_weights[j]*p_a)
                    else
                        ped_probs[index_itp_state] += itp_ped_weights[j]*p_a
                    end
                end
            end
        end
        @assert length(ped_probs) == length(ped_ps)
        normalize!(ped_probs, 1)
    end
    return ped_ps, ped_probs
end

function car_transition(mdp::PedCarMDP, s::PedCarMDPState, a::PedCarMDPAction)
    # car transition
    if s.car.posG == mdp.off_grid
        car_ps, car_routes = car_starting_states(mdp)
        push!(car_ps, get_off_the_grid(mdp))
        push!(car_routes, SVector{2, LaneTag}(LaneTag(0,0), LaneTag(0,0)))
        car_probs = ones(length(car_ps))
        car_probs[end] = 1-mdp.car_birth
        car_probs[1:end-1] = mdp.car_birth/(length(car_ps)-1)
        @assert sum(car_probs) ≈ 1.0
    else
        # update model
        set_car_model!(mdp, s, a)
        lane = get_lane(mdp.env.roadway, s.car)
        car_actions, car_actions_probs = get_distribution(mdp.car_models[s.route])
        car_ps = VehicleState[]
        car_probs = Float64[]
        car_routes = SVector{2, LaneTag}[]
        for (i, car_a) in enumerate(car_actions)
            p_a = car_actions_probs[i]
            car_p = propagate(s.car, car_a, mdp.env.roadway, mdp.ΔT)
            if car_p.posF.s >= get_end(get_lane(mdp.env.roadway, car_p)) && isempty(get_lane(mdp.env.roadway, car_p).exits)
                car_pp = get_off_the_grid(mdp)
                index_car_pp = findall(x -> x==car_pp, car_ps)
                if isempty(index_car_pp)
                    push!(car_ps, car_pp)
                    push!(car_probs, p_a)
                else
                    car_probs[index_car_pp] += p_a
                end
                car_routes = fill(s.route, length(car_ps))
                continue
            end
            # interpolate car_p in discrete space
            itp_car_ps, itp_car_weights = interpolate_state(mdp, car_p)
            for (j, car_pss) in enumerate(itp_car_ps)
                index_itp_state = findall(x -> x==car_pss, car_ps)
                if !(itp_car_weights[j] ≈ 0.)
                    if isempty(index_itp_state)
                        push!(car_ps, car_pss)
                        push!(car_probs, itp_car_weights[j]*p_a)
                    else
                        car_probs[index_itp_state] += itp_car_weights[j]*p_a
                    end
                end
            end
            car_routes = fill(s.route, length(car_ps))
        end
        @assert length(car_probs) == length(car_ps)
        @assert length(car_routes) == length(car_ps)
        normalize!(car_probs, 1)
    end
    return car_ps, car_routes, car_probs
end

function set_car_model!(mdp::PedCarMDP, s::PedCarMDPState, a::PedCarMDPAction)
    # car = s.car
    # lane = get_lane(mdp.env.roadway, car)
    # route = [mdp.env.roadway[tag] for tag in find_route(mdp.env, s.route)]
    # intersection_entrances = get_start_lanes(mdp.env.roadway)
    # if !(route[1] ∈ intersection_entrances)
    #     intersection = Lane[]
    #     intersection_exits = Lane[]
    # else
    #     intersection_exits = get_exit_lanes(mdp.env.roadway)
    #     intersection=Lane[route[1], route[2]]
    # end
    # navigator = RouteFollowingIDM(route=route, a_max=2., σ=1.)
    # intersection_driver = StopIntersectionDriver(navigator= navigator,
    #                                             intersection=intersection,
    #                                             intersection_entrances = intersection_entrances,
    #                                             intersection_exits = intersection_exits,
    #                                             stop_delta=maximum(mdp.env.params.crosswalk_width),
    #                                             accel_tol=0.,
    #                                             priorities = mdp.env.priorities)
    # # intersection_driver = TTCIntersectionDriver(navigator = navigator,
    # #                                         intersection = intersection,
    # #                                         intersection_pos = VecSE2(mdp.env.params.inter_x,
    # #                                                                     mdp.env.params.inter_y),
    # #                                         stop_delta = maximum(mdp.env.params.crosswalk_width),
    # #                                         accel_tol = 0.,
    # #                                         priorities = mdp.env.priorities,
    # #                                         ttc_threshold = (mdp.env.params.x_max - mdp.env.params.inter_x)/mdp.env.params.speed_limit
    # #                                         )
    # crosswalk_drivers = Vector{CrosswalkDriver}(length(mdp.env.crosswalks))
    # for i=1:length(mdp.env.crosswalks)
    #     cw_conflict_lanes = get_conflict_lanes(mdp.env.crosswalks[i], mdp.env.roadway)
    #     crosswalk_drivers[i] = CrosswalkDriver(navigator = navigator,
    #                             crosswalk = mdp.env.crosswalks[i],
    #                             conflict_lanes = cw_conflict_lanes,
    #                             intersection_entrances = intersection_entrances,
    #                             yield=!isempty(intersect(cw_conflict_lanes, route))
    #                             )
    #     # println(" yield to cw ", i, " ", crosswalk_drivers[i].yield)
    # end
    # mdp.car_model = UrbanDriver(navigator=navigator,
    #                                         intersection_driver=intersection_driver,
    #                                         crosswalk_drivers=crosswalk_drivers
    #                                         )
    # update model states
    # scene = Scene()
    # push!(scene, Vehicle(s.ego, mdp.ego_type, EGO_ID))
    # push!(scene, Vehicle(s.car, mdp.car_type, CAR_ID))
    # push!(scene, Vehicle(s.ped, mdp.ped_type, PED_ID))
    reset_hidden_state!(mdp.car_models[s.route])
    observe!(mdp.car_models[s.route], s.ego, s.car, s.ped, mdp.env.roadway)
    # # fill in hidden state
    # observe!(mdp.car_models[s.route], scene, mdp.env.roadway, CAR_ID)
    # # set the decision
    # observe!(mdp.car_models[s.route], scene, mdp.env.roadway, CAR_ID)
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
#---------------------------------------------------------------------#
#This function computes the high-order grid & elements in 2D.
#Written by F.X. Giraldo on 5/2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#

using Plots
include("warp_grid.jl")

function create_grid(Np,Npoin,Ne,Nboun,Nex,Ney,ξ,plot_grid,warp_grid,DFloat)

    #Initialize Global Arrays
    coord=zeros(DFloat,2,Npoin)
    intma=zeros(Int64,Np,Np,Ne)
    iboun=zeros(Int64,Nboun,1)

    #Initialize Local Arrays
    node=zeros(DFloat,Npoin,Npoin)

    #Set some constants
    xmin=-1
    xmax=+1
    ymin=-1
    ymax=+1
    dx=(xmax-xmin)/Nex
    dy=(ymax-ymin)/Ney
    N=Np-1
    nx=Nex*N + 1
    ny=Ney*N + 1

    #GENERATE COORD
    I=0
    jj=0
    for k=1:Ney
        y0=ymin + (k-1)*dy
        if (k == 1)
            l1=1
        else
            l1=2
        end
        for l=l1:Np
            y=( ξ[l]+1 )*dy/2 + y0
            jj=jj+1
            ii=0
            for i=1:Nex
                x0=xmin + (i-1)*dx
                if (i == 1)
                    j1=1
                else
                    j1=2
                end
                for j=j1:Np
                    ii=ii + 1
                    I=I + 1
                    x=( ξ[j]+1 )*dx/2 + x0
                    coord[1,I]=x
                    coord[2,I]=y
                    node[ii,jj]=I
                end #j

            end #i
        end #l
    end #k

    #GENERATE INTMA
    e=0
    for k=1:Ney
        for i=1:Nex
            e=e+1
            for l=1:Np
                jj=N*(k-1) + l
                for j=1:Np
                    ii=N*(i-1) + j
                    I=node[ii,jj]
                    intma[j,l,e]=I
                end #j
            end #l
        end #i
    end #k

    #GENERATE Boundary Conditions
    b=0
    for i=1:nx
        I=node[i,1]
        b=b + 1
        iboun[b]=I
    end
    for j=2:ny
        I=node[nx,j]
        b=b + 1
        iboun[b]=I
    end
    for i=nx-1:-1:1
        I=node[i,ny]
        b=b + 1
        iboun[b]=I
    end
    for j=ny-1:-1:2
        I=node[1,j]
        b=b + 1
        iboun[b]=I
    end

    if (b != Nboun)
        disp("Error in Create_Grid")
        return
    end

    #Warp Grid
    if (warp_grid)
        warp_grid!(coord,Npoin,DFloat)
    end

    #Plot Grid
    title_text=["Grid : Ne =  $Ne,  N = $N"]
    if (plot_grid)
        x=zeros(DFloat,5,1)
        y=zeros(DFloat,5,1)
        x[1]=-1; y[1]=-1
        x[2]=+1; y[2]=-1
        x[3]=+1; y[3]=+1
        x[4]=-1; y[4]=+1
        x[5]=-1; y[5]=-1
        p1=plot(x,y,linecolor="black",legend=false)

        for e=1:Ne
            for j=1:N
                for i=1:N
                    i1=intma[i,j,e]
                    i2=intma[i+1,j,e]
                    i3=intma[i+1,j+1,e]
                    i4=intma[i,j+1,e]
                    x[1]=coord[1,i1]; y[1]=coord[2,i1]
                    x[2]=coord[1,i2]; y[2]=coord[2,i2]
                    x[3]=coord[1,i3]; y[3]=coord[2,i3]
                    x[4]=coord[1,i4]; y[4]=coord[2,i4]
                    x[5]=coord[1,i1]; y[5]=coord[2,i1]
                    p1=plot!(p1,x,y,linecolor="red",legend=false)
                end
            end
            i1=intma[1,1,e]
            i2=intma[Np,1,e]
            i3=intma[Np,Np,e]
            i4=intma[1,Np,e]
            x[1]=coord[1,i1]; y[1]=coord[2,i1]
            x[2]=coord[1,i2]; y[2]=coord[2,i2]
            x[3]=coord[1,i3]; y[3]=coord[2,i3]
            x[4]=coord[1,i4]; y[4]=coord[2,i4]
            x[5]=coord[1,i1]; y[5]=coord[2,i1]
            p1=plot!(p1,x,y,linecolor="blue",legend=false)
        end
        display(plot(p1))
    end

    return(coord,intma,iboun)
end


text:
### A Pluto.jl notebook ###
# v0.17.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ ef41f6a0-60c0-4a86-9b6e-199b191309f1
md"""
begin

	using Pkg
	Pkg.add("Revise")
	using Revise
	Pkg.add(["HypertextLiteral","GridVisualize","PlutoUI"])
	Pkg.develop("PlutoVista")
end
""";

# ╔═╡ 60941eaa-1aea-11eb-1277-97b991548781
begin
    using PlutoUI
	using GridVisualize
	using PlutoVista
	using HypertextLiteral
	GridVisualize.default_plotter!(PlutoVista)
end;

# ╔═╡ 3ea2e28c-a323-4178-a0cf-b6b61d93672a
md"""
# Widening cells
Some js-bond-css magic from [Discourse](https://discourse.julialang.org/t/cell-width-in-pluto-notebook/49761/3)
"""

# ╔═╡ 04aa3674-60d9-40b0-9e88-5d2e14b42a0b
@bind screenWidth  ScreenWidthGrabber()

# ╔═╡ 80753ad4-3c10-41ce-ad70-6dad81266273
screenWidth

# ╔═╡ 8293eaf2-5664-4e66-bf62-f73689e3b747
PlutoCellWidener(0.9*screenWidth)

# ╔═╡ 49e13e08-3182-459c-9eb4-5d95b055a597
X0=0:0.25:10

# ╔═╡ a47d5ef4-d386-47db-8cac-965860d7ca32
scalarplot(X0,X0,X0, (x,y,z)->(sin(x)*sin(y)*sin(z)*sqrt(x*y*z)),resolution=(900,900),colormap=:rainbow)

# ╔═╡ 60a4d56f-9f53-4093-857b-62a05f2a3649
htl"""

$(scalarplot(X0,X0,X0, (x,y,z)->(sin(0.1*x*y)*sin(z)),resolution=(200,200),colormap=:hot))

$(scalarplot(X0,X0, (x,y)->sin(0.1*x*y),resolution=(200,200),colormap=:summer))

$(scalarplot(X0,X0, (x,y)->sin(0.1*x*y),resolution=(200,200),colormap=:summer,backend=:plotly))

$(scalarplot(X0, (x)->x*sin(2x),color=:red,resolution=(200,200),colormap=:summer))

"""

# ╔═╡ Cell order:
# ╠═60941eaa-1aea-11eb-1277-97b991548781
# ╟─3ea2e28c-a323-4178-a0cf-b6b61d93672a
# ╠═04aa3674-60d9-40b0-9e88-5d2e14b42a0b
# ╠═80753ad4-3c10-41ce-ad70-6dad81266273
# ╠═8293eaf2-5664-4e66-bf62-f73689e3b747
# ╟─49e13e08-3182-459c-9eb4-5d95b055a597
# ╠═a47d5ef4-d386-47db-8cac-965860d7ca32
# ╠═60a4d56f-9f53-4093-857b-62a05f2a3649
# ╠═ef41f6a0-60c0-4a86-9b6e-199b191309f1


text:
module _MSobject

export MSobject
    """
        Mobj = MSobject()

    Returns a multiscale entropy object (`Mobj`) based on that originally
    proposed by Costa et. al. (2002) using the following default parameters:
    EnType = SampEn(), embedding dimension = 2, time delay = 1,
    radius = 0.2*SD(`Sig`), logarithm = natural

        Mobj = MSobject(EnType::Function)

    Returns a multiscale entropy object by passing the entropy function
    (`EnType`) and the specifying default parameters for that entropy function.
    To see the default parameters for a particular entropy method,
    type:  `? EntropyHub.EnType`  \n
    (e.g.  `? EntropyHub.SampEn`)

        Mobj = MSobject(EnType::Function; kwargs...)

    Returns a multiscale entropy object using the specified entropy method
    (`EnType`) and the 'keyword' parameters for that particular method.
    To see the default parameters for a particular entropy method,
    type:   ? EntropyHub.EnType   (e.g.  ? EntropyHub.SampEn)

    EnType can be any of the following entropy functions:

    # Base Entropies:
        -----------------
        `ApEn`      - Approximate Entropy  \n
        `SampEn`    - Sample Entropy   \n
        `FuzzEn`    - Fuzzy Entropy    \n
        `K2En`      - Kolmogorov Entropy   \n
        `PermEn`    - Permutation Entropy	 \n
        `CondEn`    - Conditional Entropy	   \n
        `DistEn`    - Distribution Entropy	    \n
        `DispEn`    - Dispersion Entropy	    \n
        `SpecEn`    - Spectral Entropy   \n
        `SyDyEn`    - Symbolic Dynamic Entropy	  \n
        `IncrEn`    - Increment Entropy	    \n
        `CoSiEn`    - Cosine Similarity Entropy	\n
        `PhasEn`    - Phase Entropy	    \n
        `SlopEn`    - Slope Entropy       \n
        `BubbEn`    - Bubble Entropy   \n
        `GridEn`    - Gridded Distribution Entropy  \n
        `EnofEn`    - Entropy of Entropy	\n
        `AttnEn`    - Attention Entropy    \n

    # Cross Entropies:
        ------------------
        `XApEn`     - Cross-Approximate Entropy    \n
        `XSampEn`   - Cross-Sample Entropy  \n
        `XFuzzEn`   - Cross-Fuzzy Entropy   \n
        `XK2En`     - Cross-Kolmogorov Entropy  \n
        `XPermEn`   - Cross-Permutation Entropy   \n
        `XCondEn`   - Cross-Conditional Entropy    \n
        `XDistEn`   - Cross-Distribution Entropy    \n
        `XSpecEn`   - Cross-Spectral Entropy   \n

    # See also `MSEn`, `XMSEn`, `rMSEn`, `cMSEn`, `hMSEn`, `rXMSEn`, `cXMSEn`, `hXMSEn`

    """
    function MSobject(EnType::Function=SampEn; kwargs...)
    #function MSobject(EnType::Function= EntropyHub.SampEn; kwargs...)

    Chk = ["ApEn";"SampEn";"FuzzEn";"K2En";"PermEn";"CondEn";"DistEn";
        "DispEn";"SyDyEn";"IncrEn";"CoSiEn";"PhasEn";"SpecEn";"SlopEn";
        "GridEn";"BubbEn";"EnofEn";"AttnEn";"XApEn";"XSampEn";"XFuzzEn";
        "XPermEn";"XCondEn";"XDistEn";"XSpecEn";"XK2En"]
    (String(Symbol(EnType)) in Chk) ? nothing :
        error("EnType:      must be a valid entropy function name.
        For more info, type:   julia>  ? EntropyHub.MSobject")

    #Mobj = (Func=getfield(Main.EntropyHub, Symbol(EnType)), kwargs...)

    Mobj = (Func= EnType, kwargs...)

    return Mobj
    end

end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""