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

# Formulate the kinetic rate vector -
rate_constant_array = data_dictionary["RATE_CONSTANT_ARRAY"];
saturation_constant_array = data_dictionary["SATURATION_CONSTANT_ARRAY"];
rate_vector = Float64[];

# 1 TGFB12_binding: TGFB12+R =([])=> TGFB12_R
tmp = rate_constant_array[1]*(TGFB12)/(saturation_constant_array[1,91] + TGFB12)*(R)/(saturation_constant_array[1,90] + R);
push!(rate_vector,tmp);

# 2 TGFB12_binding: TGFB12_R =([])=> TGFB12+R (reverse)
tmp = rate_constant_array[2]*(TGFB12_R)/(saturation_constant_array[2,94] + TGFB12_R);
push!(rate_vector,tmp);

# 3 TGFB3_binding: TGFB3+R =([])=> TGFB3_R
tmp = rate_constant_array[3]*(TGFB3)/(saturation_constant_array[3,92] + TGFB3)*(R)/(saturation_constant_array[3,90] + R);
push!(rate_vector,tmp);

# 4 TGFB3_binding: TGFB3_R =([])=> TGFB3+R (reverse)
tmp = rate_constant_array[4]*(TGFB3_R)/(saturation_constant_array[4,95] + TGFB3_R);
push!(rate_vector,tmp);

# 5 VEGF_binding: VEGF+VEGFR =([])=> VEGF_VEGFR
tmp = rate_constant_array[5]*(VEGF)/(saturation_constant_array[5,93] + VEGF)*(VEGFR)/(saturation_constant_array[5,96] + VEGFR);
push!(rate_vector,tmp);

# 6 VEGF_binding: VEGF_VEGFR =([])=> VEGF+VEGFR (reverse)
tmp = rate_constant_array[6]*(VEGF_VEGFR)/(saturation_constant_array[6,97] + VEGF_VEGFR);
push!(rate_vector,tmp);

# 7 VEGF_VEGFR_PLCy: PLCy =(VEGF_VEGFR)=> Active_PLCy
tmp = rate_constant_array[7]*(VEGF_VEGFR)*(PLCy)/(saturation_constant_array[7,29] + PLCy);
push!(rate_vector,tmp);

# 8 PLCy_activation_Calcineurin: Calcineurin =(Active_PLCy)=> Active_Calcineurin
tmp = rate_constant_array[8]*(Active_PLCy)*(Calcineurin)/(saturation_constant_array[8,31] + Calcineurin);
push!(rate_vector,tmp);

# 9 Calcineurin_activation_NFATc: NFATc =(Active_Calcineurin)=> Active_NFATc
tmp = rate_constant_array[9]*(Active_Calcineurin)*(NFATc)/(saturation_constant_array[9,33] + NFATc);
push!(rate_vector,tmp);

# 10 PLCy_activatation_PKC: PKC =(Active_PLCy)=> Active_PKC
tmp = rate_constant_array[10]*(Active_PLCy)*(PKC)/(saturation_constant_array[10,35] + PKC);
push!(rate_vector,tmp);

# 11 PKC_activation_RAF: RAF =(Active_PKC)=> RAF_P
tmp = rate_constant_array[11]*(Active_PKC)*(RAF)/(saturation_constant_array[11,37] + RAF);
push!(rate_vector,tmp);

# 12 TGFB12_R_phosphorylation_Smad: Smad =(TGFB12_R)=> pSmad
tmp = rate_constant_array[12]*(TGFB12_R)*(Smad)/(saturation_constant_array[12,39] + Smad);
push!(rate_vector,tmp);

# 13 TGFB3_R_phosphorylation_Smad: Smad =(TGFB3_R)=> pSmad
tmp = rate_constant_array[13]*(TGFB3_R)*(Smad)/(saturation_constant_array[13,39] + Smad);
push!(rate_vector,tmp);

# 14 TGFB3_R_phosphorylation_PI3K: PI3K =(TGFB3_R)=> PI3K_P
tmp = rate_constant_array[14]*(TGFB3_R)*(PI3K)/(saturation_constant_array[14,41] + PI3K);
push!(rate_vector,tmp);

# 15 PI3K_phosphorylation_GSK3: GSK3 =(PI3K_P)=> GSK3_P
tmp = rate_constant_array[15]*(PI3K_P)*(GSK3)/(saturation_constant_array[15,43] + GSK3);
push!(rate_vector,tmp);

# 16 TGFB12_R_activation_RAS: RAS =(TGFB12_R)=> RAS_GTP
tmp = rate_constant_array[16]*(TGFB12_R)*(RAS)/(saturation_constant_array[16,45] + RAS);
push!(rate_vector,tmp);

# 17 RAS_GTP_phosphorylation_RAF: RAF =(RAS_GTP)=> RAF_P
tmp = rate_constant_array[17]*(RAS_GTP)*(RAF)/(saturation_constant_array[17,37] + RAF);
push!(rate_vector,tmp);

# 18 RAF_P_phosphorylation_MAPK: MAPK =(RAF_P)=> pMAPK
tmp = rate_constant_array[18]*(RAF_P)*(MAPK)/(saturation_constant_array[18,47] + MAPK);
push!(rate_vector,tmp);

# 19 MAPK_phosphorylation_Smad: Smad =(pMAPK)=> pSmad
tmp = rate_constant_array[19]*(pMAPK)*(Smad)/(saturation_constant_array[19,39] + Smad);
push!(rate_vector,tmp);

# 20 MAPK_phosphorylation_AP1_SP1: AP1_SP1 =(pMAPK)=> AP1_SP1_P
tmp = rate_constant_array[20]*(pMAPK)*(AP1_SP1)/(saturation_constant_array[20,49] + AP1_SP1);
push!(rate_vector,tmp);

# 21 Binding_Ecadherin_Bcatenin: Ecadherin+Bcatenin =([])=> Ecadherin_Bcatenin
tmp = rate_constant_array[21]*(Ecadherin)/(saturation_constant_array[21,51] + Ecadherin)*(Bcatenin)/(saturation_constant_array[21,52] + Bcatenin);
push!(rate_vector,tmp);

# 22 Binding_Ecadherin_Bcatenin: Ecadherin_Bcatenin =([])=> Ecadherin+Bcatenin (reverse)
tmp = rate_constant_array[22]*(Ecadherin_Bcatenin)/(saturation_constant_array[22,53] + Ecadherin_Bcatenin);
push!(rate_vector,tmp);

# 23 Binding_APC_AXIN_forw: APC_AXIN+Bcatenin =([])=> APC_AXIN_Bcatenin
tmp = rate_constant_array[23]*(APC_AXIN)/(saturation_constant_array[23,54] + APC_AXIN)*(Bcatenin)/(saturation_constant_array[23,52] + Bcatenin);
push!(rate_vector,tmp);

# 24 Binding_APC_AXIN_revs: APC_AXIN_Bcatenin =([])=> APC_AXIN+Bcatenin
tmp = rate_constant_array[24]*(APC_AXIN_Bcatenin)/(saturation_constant_array[24,55] + APC_AXIN_Bcatenin);
push!(rate_vector,tmp);

# 25 TCF4_Bcatenin_complex_formation: TCF4+Bcatenin =([])=> Active_TCF4
tmp = rate_constant_array[25]*(TCF4)/(saturation_constant_array[25,56] + TCF4)*(Bcatenin)/(saturation_constant_array[25,52] + Bcatenin);
push!(rate_vector,tmp);

# 26 TCF4_Bcatenin_complex_formation: Active_TCF4 =([])=> TCF4+Bcatenin (reverse)
tmp = rate_constant_array[26]*(Active_TCF4)/(saturation_constant_array[26,57] + Active_TCF4);
push!(rate_vector,tmp);

# 27 LEF1_Bcatenin_complex_formation: LEF1+Bcatenin =([])=> Active_LEF1
tmp = rate_constant_array[27]*(LEF1)/(saturation_constant_array[27,58] + LEF1)*(Bcatenin)/(saturation_constant_array[27,52] + Bcatenin);
push!(rate_vector,tmp);

# 28 LEF1_Bcatenin_complex_formation: Active_LEF1 =([])=> LEF1+Bcatenin (reverse)
tmp = rate_constant_array[28]*(Active_LEF1)/(saturation_constant_array[28,59] + Active_LEF1);
push!(rate_vector,tmp);

# 29 deactivation_pSmad: pSmad =(PPase)=> Smad
tmp = rate_constant_array[29]*(PPase)*(pSmad)/(saturation_constant_array[29,40] + pSmad);
push!(rate_vector,tmp);

# 30 deactivation_RAF_P: RAF_P =(PPase)=> RAF
tmp = rate_constant_array[30]*(PPase)*(RAF_P)/(saturation_constant_array[30,38] + RAF_P);
push!(rate_vector,tmp);

# 31 deactivation_RAS_GTP: RAS_GTP =([])=> RAS
tmp = rate_constant_array[31]*(RAS_GTP)/(saturation_constant_array[31,46] + RAS_GTP);
push!(rate_vector,tmp);

# 32 deactivation_pMAPK: pMAPK =(PPase)=> MAPK
tmp = rate_constant_array[32]*(PPase)*(pMAPK)/(saturation_constant_array[32,48] + pMAPK);
push!(rate_vector,tmp);

# 33 deactivation_AP1_SP1: AP1_SP1_P =(PPase)=> AP1_SP1
tmp = rate_constant_array[33]*(PPase)*(AP1_SP1_P)/(saturation_constant_array[33,50] + AP1_SP1_P);
push!(rate_vector,tmp);

# 34 deactivation_Active_PLCy: Active_PLCy =([])=> PLCy
tmp = rate_constant_array[34]*(Active_PLCy)/(saturation_constant_array[34,30] + Active_PLCy);
push!(rate_vector,tmp);

# 35 deactivation_Active_PKC: Active_PKC =([])=> PKC
tmp = rate_constant_array[35]*(Active_PKC)/(saturation_constant_array[35,36] + Active_PKC);
push!(rate_vector,tmp);

# 36 deactivation_Active_Calcineurin: Active_Calcineurin =([])=> Calcineurin
tmp = rate_constant_array[36]*(Active_Calcineurin)/(saturation_constant_array[36,32] + Active_Calcineurin);
push!(rate_vector,tmp);

# 37 deactivation_Active_NFATc: Active_NFATc =([])=> NFATc
tmp = rate_constant_array[37]*(Active_NFATc)/(saturation_constant_array[37,34] + Active_NFATc);
push!(rate_vector,tmp);

# 38 deactivation_GSK3_P: GSK3_P =(PPase)=> GSK3
tmp = rate_constant_array[38]*(PPase)*(GSK3_P)/(saturation_constant_array[38,44] + GSK3_P);
push!(rate_vector,tmp);

# 39 degradation_TGFB12_R: TGFB12_R =([])=> []
tmp = rate_constant_array[39]*(TGFB12_R);
push!(rate_vector,tmp);

# 40 degradation_TGFB3_R: TGFB3_R =([])=> []
tmp = rate_constant_array[40]*(TGFB3_R);
push!(rate_vector,tmp);

# 41 degradation_TGFB3: TGFB3 =([])=> []
tmp = rate_constant_array[41]*(TGFB3);
push!(rate_vector,tmp);

# 42 VEGF_VEGFR_degradation: VEGF_VEGFR =([])=> []
tmp = rate_constant_array[42]*(VEGF_VEGFR);
push!(rate_vector,tmp);

# 43 Bcatenin_degradation: Bcatenin =([])=> []
tmp = rate_constant_array[43]*(Bcatenin);
push!(rate_vector,tmp);

# 44 Ecadherin_degrdation: Ecadherin =([])=> []
tmp = rate_constant_array[44]*(Ecadherin);
push!(rate_vector,tmp);

# 45 Smad_degradation: Smad =([])=> []
tmp = rate_constant_array[45]*(Smad);
push!(rate_vector,tmp);

# 46 TCF4_degradation: TCF4 =([])=> []
tmp = rate_constant_array[46]*(TCF4);
push!(rate_vector,tmp);

# 47 Active_TCF4_degradation: Active_TCF4 =([])=> []
tmp = rate_constant_array[47]*(Active_TCF4);
push!(rate_vector,tmp);

# 48 YREG1_degradation: YREG1 =([])=> []
tmp = rate_constant_array[48]*(YREG1);
push!(rate_vector,tmp);

# 49 SNAIL_SLUG_degradation: SNAIL_SLUG =([])=> []
tmp = rate_constant_array[49]*(SNAIL_SLUG);
push!(rate_vector,tmp);

# 50 LEF1_degradation: LEF1 =([])=> []
tmp = rate_constant_array[50]*(LEF1);
push!(rate_vector,tmp);

# 51 Active_LEF1_degradation: Active_LEF1 =([])=> []
tmp = rate_constant_array[51]*(Active_LEF1);
push!(rate_vector,tmp);

# 52 vimentin_degradation: vimentin =([])=> []
tmp = rate_constant_array[52]*(vimentin);
push!(rate_vector,tmp);

# 53 VEGFR_degradation: VEGFR =([])=> []
tmp = rate_constant_array[53]*(VEGFR);
push!(rate_vector,tmp);

# 54 NFATc_degradation: NFATc =([])=> []
tmp = rate_constant_array[54]*(NFATc);
push!(rate_vector,tmp);

# 55 Active_NFATc_degradation: Active_NFATc =([])=> []
tmp = rate_constant_array[55]*(Active_NFATc);
push!(rate_vector,tmp);

# 56 degradation_MAPK: MAPK =([])=> []
tmp = rate_constant_array[56]*(MAPK);
push!(rate_vector,tmp);

# 57 degradation_R: R =([])=> []
tmp = rate_constant_array[57]*(R);
push!(rate_vector,tmp);

# 58 Ecadherin_Bcatenin_degradation: Ecadherin_Bcatenin =([])=> []
tmp = rate_constant_array[58]*(Ecadherin_Bcatenin);
push!(rate_vector,tmp);

# 59 AP1_SP1_degradation: AP1_SP1 =([])=> []
tmp = rate_constant_array[59]*(AP1_SP1);
push!(rate_vector,tmp);

# 60 AP1_SP1_P_degradation: AP1_SP1_P =([])=> []
tmp = rate_constant_array[60]*(AP1_SP1_P);
push!(rate_vector,tmp);

# 61 RAS_degradation: RAS =([])=> []
tmp = rate_constant_array[61]*(RAS);
push!(rate_vector,tmp);

# 62 RAF_degradation: RAF =([])=> []
tmp = rate_constant_array[62]*(RAF);
push!(rate_vector,tmp);

# 63 RAS_GTP_degradation: RAS_GTP =([])=> []
tmp = rate_constant_array[63]*(RAS_GTP);
push!(rate_vector,tmp);

# 64 RAF_P_degradation: RAF_P =([])=> []
tmp = rate_constant_array[64]*(RAF_P);
push!(rate_vector,tmp);

# 65 pMAPK_degradation: pMAPK =([])=> []
tmp = rate_constant_array[65]*(pMAPK);
push!(rate_vector,tmp);

# 66 pSmad_degradation: pSmad =([])=> []
tmp = rate_constant_array[66]*(pSmad);
push!(rate_vector,tmp);

# 67 PLCy_degradation: PLCy =([])=> []
tmp = rate_constant_array[67]*(PLCy);
push!(rate_vector,tmp);

# 68 Active_PLCy_degradation: Active_PLCy =([])=> []
tmp = rate_constant_array[68]*(Active_PLCy);
push!(rate_vector,tmp);

# 69 PKC_degradation: PKC =([])=> []
tmp = rate_constant_array[69]*(PKC);
push!(rate_vector,tmp);

# 70 Active_PKC_degradation: Active_PKC =([])=> []
tmp = rate_constant_array[70]*(Active_PKC);
push!(rate_vector,tmp);

# 71 Calcineurin_degradation: Calcineurin =([])=> []
tmp = rate_constant_array[71]*(Calcineurin);
push!(rate_vector,tmp);

# 72 Active_Calcineurin_degradation: Active_Calcineurin =([])=> []
tmp = rate_constant_array[72]*(Active_Calcineurin);
push!(rate_vector,tmp);

# 73 degradation_APC_AXIN: APC_AXIN =([])=> []
tmp = rate_constant_array[73]*(APC_AXIN);
push!(rate_vector,tmp);

# 74 degradation_GSK3_P: GSK3_P =([])=> []
tmp = rate_constant_array[74]*(GSK3_P);
push!(rate_vector,tmp);

# 75 degradation_GSK3: GSK3 =([])=> []
tmp = rate_constant_array[75]*(GSK3);
push!(rate_vector,tmp);

# 76 degradation_PI3K: PI3K =([])=> []
tmp = rate_constant_array[76]*(PI3K);
push!(rate_vector,tmp);

# 77 degradation_APC_AXIN_Bcatenin: APC_AXIN_Bcatenin =([])=> []
tmp = rate_constant_array[77]*(APC_AXIN_Bcatenin);
push!(rate_vector,tmp);

# 78 tubulin_degradation: tubulin =([])=> []
tmp = rate_constant_array[78]*(tubulin);
push!(rate_vector,tmp);

# 79 generation_RNAP: RNAP =([])=> []
tmp = rate_constant_array[79]*(RNAP);
push!(rate_vector,tmp);

# 80 generation_RNAP: [] =([])=> RNAP (reverse)
tmp = rate_constant_array[80];
push!(rate_vector,tmp);

# 81 generation_RIBOSOME: RIBOSOME =([])=> []
tmp = rate_constant_array[81]*(RIBOSOME);
push!(rate_vector,tmp);

# 82 generation_RIBOSOME: [] =([])=> RIBOSOME (reverse)
tmp = rate_constant_array[82];
push!(rate_vector,tmp);

# 83 generation_PPase: PPase =([])=> []
tmp = rate_constant_array[83]*(PPase);
push!(rate_vector,tmp);

# 84 generation_PPase: [] =([])=> PPase (reverse)
tmp = rate_constant_array[84];
push!(rate_vector,tmp);

# 85 addition_TGFB3_Ab: TGFB3_Ab =([])=> []
tmp = rate_constant_array[85]*(TGFB3_Ab);
push!(rate_vector,tmp);

# 86 addition_TGFB3_Ab: [] =([])=> TGFB3_Ab (reverse)
tmp = rate_constant_array[86];
push!(rate_vector,tmp);

# 87 addition_ss_virus: ss_virus =([])=> []
tmp = rate_constant_array[87]*(ss_virus);
push!(rate_vector,tmp);

# 88 addition_ss_virus: [] =([])=> ss_virus (reverse)
tmp = rate_constant_array[88];
push!(rate_vector,tmp);

# 89 addition_ecad_virus: ecad_virus =([])=> []
tmp = rate_constant_array[89]*(ecad_virus);
push!(rate_vector,tmp);

# 90 addition_ecad_virus: [] =([])=> ecad_virus (reverse)
tmp = rate_constant_array[90];
push!(rate_vector,tmp);

# 91 addition_TGFb12: [] =([])=> TGFB12
tmp = rate_constant_array[91];
push!(rate_vector,tmp);

# 92 addition_TGFb12: TGFB12 =([])=> [] (reverse)
tmp = rate_constant_array[92]*(TGFB12);
push!(rate_vector,tmp);

# 93 addition_VEGF: [] =([])=> VEGF
tmp = rate_constant_array[93];
push!(rate_vector,tmp);

# 94 addition_VEGF: VEGF =([])=> [] (reverse)
tmp = rate_constant_array[94]*(VEGF);
push!(rate_vector,tmp);

# 95 addition_TGFB3: [] =([])=> TGFB3
tmp = rate_constant_array[95];
push!(rate_vector,tmp);

# 96 addition_TGFB3: TGFB3 =([])=> [] (reverse)
tmp = rate_constant_array[96]*(TGFB3);
push!(rate_vector,tmp);

# 97 addition_U0126: U0126 =([])=> []
tmp = rate_constant_array[97]*(U0126);
push!(rate_vector,tmp);

# 98 addition_U0126: [] =([])=> U0126 (reverse)
tmp = rate_constant_array[98];
push!(rate_vector,tmp);

# 99 addition_vivit: vivit =([])=> []
tmp = rate_constant_array[99]*(vivit);
push!(rate_vector,tmp);

# 100 addition_vivit: [] =([])=> vivit (reverse)
tmp = rate_constant_array[100];
push!(rate_vector,tmp);

# 101 degradation_mRNA_Smad: mRNA_Smad = []
tmp = rate_constant_array[101]*(mRNA_Smad);
push!(rate_vector,tmp);

# 102 degradation_mRNA_Ecadherin: mRNA_Ecadherin = []
tmp = rate_constant_array[102]*(mRNA_Ecadherin);
push!(rate_vector,tmp);

# 103 degradation_mRNA_Bcatenin: mRNA_Bcatenin = []
tmp = rate_constant_array[103]*(mRNA_Bcatenin);
push!(rate_vector,tmp);

# 104 degradation_mRNA_TCF4: mRNA_TCF4 = []
tmp = rate_constant_array[104]*(mRNA_TCF4);
push!(rate_vector,tmp);

# 105 degradation_mRNA_TGFB3: mRNA_TGFB3 = []
tmp = rate_constant_array[105]*(mRNA_TGFB3);
push!(rate_vector,tmp);

# 106 degradation_mRNA_YREG1: mRNA_YREG1 = []
tmp = rate_constant_array[106]*(mRNA_YREG1);
push!(rate_vector,tmp);

# 107 degradation_mRNA_LEF1: mRNA_LEF1 = []
tmp = rate_constant_array[107]*(mRNA_LEF1);
push!(rate_vector,tmp);

# 108 degradation_mRNA_vimentin: mRNA_vimentin = []
tmp = rate_constant_array[108]*(mRNA_vimentin);
push!(rate_vector,tmp);

# 109 degradation_mRNA_VEGFR: mRNA_VEGFR = []
tmp = rate_constant_array[109]*(mRNA_VEGFR);
push!(rate_vector,tmp);

# 110 degradation_mRNA_NFATc: mRNA_NFATc = []
tmp = rate_constant_array[110]*(mRNA_NFATc);
push!(rate_vector,tmp);

# 111 degradation_mRNA_MAPK: mRNA_MAPK = []
tmp = rate_constant_array[111]*(mRNA_MAPK);
push!(rate_vector,tmp);

# 112 degradation_mRNA_R: mRNA_R = []
tmp = rate_constant_array[112]*(mRNA_R);
push!(rate_vector,tmp);

# 113 degradation_mRNA_PLCy: mRNA_PLCy = []
tmp = rate_constant_array[113]*(mRNA_PLCy);
push!(rate_vector,tmp);

# 114 degradation_mRNA_PKC: mRNA_PKC = []
tmp = rate_constant_array[114]*(mRNA_PKC);
push!(rate_vector,tmp);

# 115 degradation_mRNA_RAF: mRNA_RAF = []
tmp = rate_constant_array[115]*(mRNA_RAF);
push!(rate_vector,tmp);

# 116 degradation_mRNA_RAS: mRNA_RAS = []
tmp = rate_constant_array[116]*(mRNA_RAS);
push!(rate_vector,tmp);

# 117 degradation_mRNA_AP1_SP1: mRNA_AP1_SP1 = []
tmp = rate_constant_array[117]*(mRNA_AP1_SP1);
push!(rate_vector,tmp);

# 118 degradation_mRNA_Calcineurin: mRNA_Calcineurin = []
tmp = rate_constant_array[118]*(mRNA_Calcineurin);
push!(rate_vector,tmp);

# 119 degradation_mRNA_SNAIL_SLUG: mRNA_SNAIL_SLUG = []
tmp = rate_constant_array[119]*(mRNA_SNAIL_SLUG);
push!(rate_vector,tmp);

# 120 degradation_mRNA_tubulin: mRNA_tubulin = []
tmp = rate_constant_array[120]*(mRNA_tubulin);
push!(rate_vector,tmp);

# 121 degradation_mRNA_APC_AXIN: mRNA_APC_AXIN = []
tmp = rate_constant_array[121]*(mRNA_APC_AXIN);
push!(rate_vector,tmp);

# 122 degradation_mRNA_GSK3: mRNA_GSK3 = []
tmp = rate_constant_array[122]*(mRNA_GSK3);
push!(rate_vector,tmp);

# 123 degradation_mRNA_PI3K: mRNA_PI3K = []
tmp = rate_constant_array[123]*(mRNA_PI3K);
push!(rate_vector,tmp);

# 124 induction_gene_Smad: gene_Smad = mRNA_Smad
tmp = rate_constant_array[124]*gene_Smad*RNAP;
push!(rate_vector,tmp);

# 125 induction_gene_Ecadherin: gene_Ecadherin = mRNA_Ecadherin
tmp = rate_constant_array[125]*gene_Ecadherin*RNAP;
push!(rate_vector,tmp);

# 126 induction_gene_Bcatenin: gene_Bcatenin = mRNA_Bcatenin
tmp = rate_constant_array[126]*gene_Bcatenin*RNAP;
push!(rate_vector,tmp);

# 127 induction_gene_TCF4: gene_TCF4 = mRNA_TCF4
tmp = rate_constant_array[127]*gene_TCF4*RNAP;
push!(rate_vector,tmp);

# 128 induction_gene_TGFB3: gene_TGFB3 = mRNA_TGFB3
tmp = rate_constant_array[128]*gene_TGFB3*RNAP;
push!(rate_vector,tmp);

# 129 induction_gene_YREG1: gene_YREG1 = mRNA_YREG1
tmp = rate_constant_array[129]*gene_YREG1*RNAP;
push!(rate_vector,tmp);

# 130 induction_gene_LEF1: gene_LEF1 = mRNA_LEF1
tmp = rate_constant_array[130]*gene_LEF1*RNAP;
push!(rate_vector,tmp);

# 131 induction_gene_R: gene_R = mRNA_R
tmp = rate_constant_array[131]*gene_R*RNAP;
push!(rate_vector,tmp);

# 132 induction_gene_vimentin: gene_vimentin = mRNA_vimentin
tmp = rate_constant_array[132]*gene_vimentin*RNAP;
push!(rate_vector,tmp);

# 133 induction_gene_VEGFR: gene_VEGFR = mRNA_VEGFR
tmp = rate_constant_array[133]*gene_VEGFR*RNAP;
push!(rate_vector,tmp);

# 134 induction_gene_NFATc: gene_NFATc = mRNA_NFATc
tmp = rate_constant_array[134]*gene_NFATc*RNAP;
push!(rate_vector,tmp);

# 135 induction_gene_PLCy: gene_PLCy = mRNA_PLCy
tmp = rate_constant_array[135]*gene_PLCy*RNAP;
push!(rate_vector,tmp);

# 136 induction_gene_PKC: gene_PKC = mRNA_PKC
tmp = rate_constant_array[136]*gene_PKC*RNAP;
push!(rate_vector,tmp);

# 137 induction_gene_RAF: gene_RAF = mRNA_RAF
tmp = rate_constant_array[137]*gene_RAF*RNAP;
push!(rate_vector,tmp);

# 138 induction_gene_RAS: gene_RAS = mRNA_RAS
tmp = rate_constant_array[138]*gene_RAS*RNAP;
push!(rate_vector,tmp);

# 139 induction_gene_MAPK: gene_MAPK = mRNA_MAPK
tmp = rate_constant_array[139]*gene_MAPK*RNAP;
push!(rate_vector,tmp);

# 140 induction_gene_AP1_SP1: gene_AP1_SP1 = mRNA_AP1_SP1
tmp = rate_constant_array[140]*gene_AP1_SP1*RNAP;
push!(rate_vector,tmp);

# 141 induction_gene_GSK3: gene_GSK3 = mRNA_GSK3
tmp = rate_constant_array[141]*gene_GSK3*RNAP;
push!(rate_vector,tmp);

# 142 induction_gene_tubulin: gene_tubulin = mRNA_tubulin
tmp = rate_constant_array[142]*gene_tubulin*RNAP;
push!(rate_vector,tmp);

# 143 induction_gene_SNAIL_SLUG: gene_SNAIL_SLUG = mRNA_SNAIL_SLUG
tmp = rate_constant_array[143]*gene_SNAIL_SLUG*RNAP;
push!(rate_vector,tmp);

# 144 induction_gene_Calcineurin: gene_Calcineurin = mRNA_Calcineurin
tmp = rate_constant_array[144]*gene_Calcineurin*RNAP;
push!(rate_vector,tmp);

# 145 induction_gene_APC_AXIN: gene_APC_AXIN = mRNA_APC_AXIN
tmp = rate_constant_array[145]*gene_APC_AXIN*RNAP;
push!(rate_vector,tmp);

# 146 induction_gene_PI3K: gene_PI3K = mRNA_PI3K
tmp = rate_constant_array[146]*gene_PI3K*RNAP;
push!(rate_vector,tmp);

# 147 translation_mRNA_Smad: mRNA_Smad = Smad
tmp = rate_constant_array[147]*mRNA_Smad*RIBOSOME;
push!(rate_vector,tmp);

# 148 translation_mRNA_Ecadherin: mRNA_Ecadherin = Ecadherin
tmp = rate_constant_array[148]*mRNA_Ecadherin*RIBOSOME;
push!(rate_vector,tmp);

# 149 translation_mRNA_Bcatenin: mRNA_Bcatenin = Bcatenin
tmp = rate_constant_array[149]*mRNA_Bcatenin*RIBOSOME;
push!(rate_vector,tmp);

# 150 translation_mRNA_TCF4: mRNA_TCF4 = TCF4
tmp = rate_constant_array[150]*mRNA_TCF4*RIBOSOME;
push!(rate_vector,tmp);

# 151 translation_mRNA_TGFB3: mRNA_TGFB3 = TGFB3
tmp = rate_constant_array[151]*mRNA_TGFB3*RIBOSOME;
push!(rate_vector,tmp);

# 152 translation_mRNA_YREG1: mRNA_YREG1 = YREG1
tmp = rate_constant_array[152]*mRNA_YREG1*RIBOSOME;
push!(rate_vector,tmp);

# 153 translation_mRNA_LEF1: mRNA_LEF1 = LEF1
tmp = rate_constant_array[153]*mRNA_LEF1*RIBOSOME;
push!(rate_vector,tmp);

# 154 translation_mRNA_vimentin: mRNA_vimentin = vimentin
tmp = rate_constant_array[154]*mRNA_vimentin*RIBOSOME;
push!(rate_vector,tmp);

# 155 translation_mRNA_VEGFR: mRNA_VEGFR = VEGFR
tmp = rate_constant_array[155]*mRNA_VEGFR*RIBOSOME;
push!(rate_vector,tmp);

# 156 translation_mRNA_NFATc: mRNA_NFATc = NFATc
tmp = rate_constant_array[156]*mRNA_NFATc*RIBOSOME;
push!(rate_vector,tmp);

# 157 translation_mRNA_MAPK: mRNA_MAPK = MAPK
tmp = rate_constant_array[157]*mRNA_MAPK*RIBOSOME;
push!(rate_vector,tmp);

# 158 translation_mRNA_R: mRNA_R = R
tmp = rate_constant_array[158]*mRNA_R*RIBOSOME;
push!(rate_vector,tmp);

# 159 translation_mRNA_PLCy: mRNA_PLCy = PLCy
tmp = rate_constant_array[159]*mRNA_PLCy*RIBOSOME;
push!(rate_vector,tmp);

# 160 translation_mRNA_PKC: mRNA_PKC = PKC
tmp = rate_constant_array[160]*mRNA_PKC*RIBOSOME;
push!(rate_vector,tmp);

# 161 translation_mRNA_RAF: mRNA_RAF = RAF
tmp = rate_constant_array[161]*mRNA_RAF*RIBOSOME;
push!(rate_vector,tmp);

# 162 translation_mRNA_RAS: mRNA_RAS = RAS
tmp = rate_constant_array[162]*mRNA_RAS*RIBOSOME;
push!(rate_vector,tmp);

# 163 translation_mRNA_AP1_SP1: mRNA_AP1_SP1 = AP1_SP1
tmp = rate_constant_array[163]*mRNA_AP1_SP1*RIBOSOME;
push!(rate_vector,tmp);

# 164 translation_mRNA_Calcineurin: mRNA_Calcineurin = Calcineurin
tmp = rate_constant_array[164]*mRNA_Calcineurin*RIBOSOME;
push!(rate_vector,tmp);

# 165 translation_mRNA_SNAIL_SLUG: mRNA_SNAIL_SLUG = SNAIL_SLUG
tmp = rate_constant_array[165]*mRNA_SNAIL_SLUG*RIBOSOME;
push!(rate_vector,tmp);

# 166 translation_mRNA_tubulin: mRNA_tubulin = tubulin
tmp = rate_constant_array[166]*mRNA_tubulin*RIBOSOME;
push!(rate_vector,tmp);

# 167 translation_mRNA_APC_AXIN: mRNA_APC_AXIN = APC_AXIN
tmp = rate_constant_array[167]*mRNA_APC_AXIN*RIBOSOME;
push!(rate_vector,tmp);

# 168 translation_mRNA_GSK3: mRNA_GSK3 = GSK3
tmp = rate_constant_array[168]*mRNA_GSK3*RIBOSOME;
push!(rate_vector,tmp);

# 169 translation_mRNA_PI3K: mRNA_PI3K = PI3K
tmp = rate_constant_array[169]*mRNA_PI3K*RIBOSOME;
push!(rate_vector,tmp);

# return the kinetics vector -
return rate_vector;
end


text:
using ChainRulesCore
using Zygote: pullback

mutable struct SupervisedObjective{L, D, I, S}
    loss::L
    data_iterable::D
    next_instance::I
    next_state::S
    function SupervisedObjective(loss::L, data_iterable::D) where {L, D}
        next = iterate(data_iterable)
        next !== nothing || error("data iterable is empty")
        instance, state = next
        new{L, D, typeof(instance), typeof(state)}(loss, data_iterable, instance, state)
    end
end

Base.show(io::IO, mo::T) where T <: SupervisedObjective = print(io, "$T($(mo.loss), $(mo.data_iterable))")

function get_instance_and_update!(f::SupervisedObjective)
    instance, state = f.next_instance, f.next_state
    next = iterate(f.data_iterable, f.next_state)
    next !== nothing || error("data iterable is exhausted")
    f.next_instance, f.next_state = next
    return instance
end

function (f::SupervisedObjective)(m)
    x, y = get_instance_and_update!(f)
    return f.loss(m(x), y)
end

# NOTE we need to explicitly implement this because Zygote doesn't like mutation
function ChainRulesCore.rrule(f::SupervisedObjective, m)
    x, y = get_instance_and_update!(f)
    out, pb = pullback(m -> f.loss(m(x), y), m)
    return out, c -> (ChainRulesCore.NoTangent(), pb(c)[1])
end

ext: jl


text:
module DigestAuthRequest

using .. Base64

import ..Layer, ..request
using URIs
using ..Pairs: getkv, setkv
import ..@debug, ..DEBUG_LEVEL

"""
    request(DigestAuthLayer, method, ::URI, headers, body) -> HTTP.Response
Add `Authorization: Digest` header using credentials from url userinfo.
"""
abstract type DigestAuthLayer{Next <: Layer} <: Layer{Next} end
export DigestAuthLayer

abstract type DigestAuthLayer{Next <: Layer} <: Layer{Next} end
export DigestAuthLayer

function request(::Type{DigestAuthLayer{Next}},
                 method::String, url::URI, headers, body; kw...) where Next

    userinfo = unescapeuri(url.userinfo)

    if !isempty(userinfo) && getkv(headers, "Authorization", "") == ""
        @debug 1 "Adding Authorization: Digest header."
        setkv(headers, "Authorization", "Digest $(base64encode(userinfo))")
    end

    return request(Next, method, url, headers, body; kw...)
end

function sign_digest!(algorithm::Any,
                    method::String,
                    url::URI,
                    headers::Headers,
                    body::Vector{UInt8};
    )
end

end # module DigestAuthRequest


text:
# updateGammaz does not recognize this as bad diamond II
# Claudia April 2015

PhyloNetworks.CHECKNET || error("need CHECKNET==true in PhyloNetworks to test snaq in test_correctLik.jl")

@testset "test: bad diamond, max pseudo lik" begin
global tree, net, df, d

tree = "(6,(5,#H7:0.0):9.970714072991349,(3,(((2,1):0.2950382234364404,4):0.036924483697671304)#H7:0.00926495670648208):1.1071489442240392);"
net = readTopologyLevel1(tree);
checkNet(net)
#printNodes(net)
#printEdges(net)
@test net.node[10].number == 3 # or: wrong hybrid
@test net.node[10].hybrid # or: does not know it is hybrid
@test net.node[10].isBadDiamondII # or: does not know it is bad diamond II
##plot(net,showEdgeNumber=true)
@test [e.inCycle for e in net.edge] == [-1,-1,3,3,-1,-1,-1,-1,-1,-1,3,3] # or: error in incycle
@test [e.containRoot for e in net.edge] == [true,true,false,true,true,false,false,false,false,false,false,true] # or: error in contain root
@test [e.istIdentifiable for e in net.edge] == [false,false,true,true,false,false,false,true,false,false,true,true] # or: istIdentifiable not correct
@test (net.edge[3].hybrid && net.edge[11].hybrid) # or: hybrid edges wrong")

df=DataFrame(t1=["1","1","2","2","1","2","2","2","2","2","1","2","2","3","2"],
             t2=["3","3","3","3","3","1","5","1","1","1","5","1","1","5","3"],
             t3=["5","5","5","5","6","5","6","3","6","5","6","3","3","6","6"],
             t4=["6","4","4","6","4","6","4","5","4","4","4","6","4","4","4"],
             CF1234=[0.565,0.0005,0.0005,0.565,0.00003,0.99986,0.0410167,1,0.99987,1,0.040167,0.998667,1,0.073167,0.00003],
             CF1324=[0.0903,0.8599,0.8599,0.0903,0.8885,0.00006,0.263,0,0.00006,0,0.2630,0.00006,0,0.0424667,0.8885])
df[!,:CF1423] = 1.0 .- df[!,:CF1234] .- df[!,:CF1324]
d = readTableCF(df)

Random.seed!(345);
net2 = topologyMaxQPseudolik!(net,d, ftolRel=1e-5,ftolAbs=1e-6,xtolRel=1e-3,xtolAbs=1e-4)
@test net2.loglik < 340.5 # loglik ~ 339.95
@test net2.edge[3].istIdentifiable # or: wrong hybrid is t identifiable
# @test 9.95 < net2.edge[3].length < 9.99
# why it's not a good idea to test the branch length estimate:
# BL~4.92 when loglik~339.95
# BL>9.95 when loglik>340
# BL~0.0  when loglik~339.88 (using tolerances of 1e-8)
@test net2.edge[11].istIdentifiable # or: wrong hybrid is t identifiable
@test net2.edge[11].length < 0.01 # or: wrong bl estimated
@test net2.edge[10].length == 0.0 # or: tree edge in bad diamond II not 0
#printEdges(net2)

@test_logs show(devnull, net2)
@test_logs [show(devnull, net2.node[i]) for i in [1,3,10]];
@test_logs [show(devnull, net2.edge[i]) for i in [1,3,11]];
@test_logs show(devnull, d)
@test_logs show(devnull, d.quartet[1])
@test_logs show(devnull, d.quartet[1].qnet)
@test tipLabels(d.quartet) == ["1","2","3","4","5","6"]
a = (@test_logs fittedQuartetCF(d, :wide));
@test size(a) == (15,10)
a = (@test_logs fittedQuartetCF(d, :long));
@test size(a) == (45,7)

end

## testing readTopology----------------------------------------------------------------
## will remove this from the test because readTopology should not have to worry about istIdentifiable
## we move into updateGammaz
## tree = "(6,(5,#H7:0.0):9.970714072991349,(3,(((2,1):0.2950382234364404,4):0.036924483697671304)#H7:0.00926495670648208):1.1071489442240392);"
## net2 = readTopology(tree)
## printEdges(net2)

## for e in net2.edge
##     if(e.node[1].leaf || e.node[2].leaf)
##         !e.istIdentifiable || error("ext edge should not identifiable")
##     else
##         e.istIdentifiable || error("int edge should not identifiable")
##     end
## end

## ## this case works fine
## tree = "((((8,10))#H1,7),6,(4,#H1));" # Case I Bad diamond I
## net = readTopologyLevel1(tree)
## checkNet(net)
## net2 = readTopology(tree)
## printEdges(net2)

## for e in net2.edge
##     if(e.node[1].leaf || e.node[2].leaf)
##         !e.istIdentifiable || error("ext edge should not identifiable")
##     else
##         e.istIdentifiable || error("int edge should not identifiable")
##     end
## end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

###############
# Prepare FMU #
###############

using FMIImport.FMICore: fmi2FMUstate

pathToFMU = "https://github.com/ThummeTo/FMI.jl/raw/main/model/" * ENV["EXPORTINGTOOL"] * "/SpringPendulum1D.fmu"

myFMU = fmi2Load(pathToFMU)

comp = fmi2Instantiate!(myFMU; loggingOn=true)
@test comp != 0

@test fmi2EnterInitializationMode(comp) == 0
@test fmi2ExitInitializationMode(comp) == 0

@test fmi2SetupExperiment(comp, 0.0) == 0

###########################
# Testing state functions #
###########################

if fmi2CanGetSetState(myFMU) && fmi2CanSerializeFMUstate(myFMU)
    @test fmi2GetReal(comp, "mass.s") == 0.5
    FMUstate = fmi2GetFMUstate(comp)
    @test typeof(FMUstate) == fmi2FMUstate
    len = fmi2SerializedFMUstateSize(comp, FMUstate)
    @test len > 0
    serial = fmi2SerializeFMUstate(comp, FMUstate)
    @test length(serial) == len
    @test typeof(serial) == Array{Char,1}

    fmi2SetReal(comp, "mass.s", 10.0)
    FMUstate = fmi2GetFMUstate(comp)
    @test fmi2GetReal(comp, "mass.s") == 10.0

    FMUstate2 = fmi2DeSerializeFMUstate(comp, serial)
    @test typeof(FMUstate2) == fmi2FMUstate
    fmi2SetFMUstate(comp, FMUstate2)
    @test fmi2GetReal(comp, "mass.s") == 0.5
    fmi2SetFMUstate(comp, FMUstate)
    @test fmi2GetReal(comp, "mass.s") == 10.0
    fmi2FreeFMUstate!(comp, FMUstate)
    fmi2FreeFMUstate!(comp, FMUstate2)
else
    @info "The FMU provided from the tool `$(ENV["EXPORTINGTOOL"])` does not support state get, set, serialization and deserialization. Skipping related tests."
end

############
# Clean up #
############

@test fmi2Terminate(comp) == 0
fmi2Unload(myFMU)


text:
using MarkovTransitionMatrices
using Test
using BenchmarkTools
using SharedArrays
using Distributed
using LinearAlgebra

using Distributions

isapprox_oneval(x::AbstractArray) = all( x .≈ x[1] )

include("FarmerTodaBmark.jl")

# setup for test
testvec = -3.0:0.5:0.0
ss = collect(Base.product(testvec, testvec))
n = length(testvec)
nn = length(ss)
mvd = MvNormal(zeros(2), Matrix{Float64}(I,2,2))

# test the 1-d version
fullp1 = Matrix(markov_transition((s) -> 0., (s) -> 1., 1e-8, testvec))
ratio1 = fullp1 ./ pdf.(Normal(0,1), testvec')
@test isapprox_oneval(ratio1)

# test the 2-d version
fullp2 = Matrix(markov_transition((s) -> zeros(2), (s) -> Matrix{Float64}(I,2,2), 1e-8, testvec, testvec))
ratio2 = reshape(fullp2[1,:], n, n) ./  [pdf(mvd, [s...]) for s in ss]
@test isapprox_oneval(ratio2)

# test the markov-switching version
μswitch(r::Real,s::NTuple{N,T}) where {N,T<:Real} = r==1 ? zeros(2) : [s...]
Σswitch(r::Real,s::NTuple{N,T}) where {N,T<:Real} = r==1 ? Matrix{Float64}(I,2,2)   : eps(1.0)*Matrix(I,2,2)
πswitch = Matrix{Float64}(I, 2, 2)
fullp3 = Matrix(markovswitching_transition(μswitch, Σswitch, πswitch, 0.0, testvec, testvec))

@test all(fullp3[1:nn, 1:nn] .== fullp2)
@test all(fullp3[nn+1:2*nn, 1:nn] .== 0.0)
@test all(fullp3[1:nn, nn+1:2*nn] .== 0.0)
@test all(fullp3[nn+1:2*nn, nn+1:2*nn] .== Matrix(I,nn,nn))

# ---------------------------------------------------------------

println("Testing moment-matching")

s = -3.0:0.25:3.0
ss = -3.0:0.125:3.0, -3.0:0.125:3.0

pids = addprocs()
@everywhere using MarkovTransitionMatrices

# simulate random walks w/ moment matching
p_P_match_1, p_mom, p_er = markov_transition_moment_matching_parallel((s) -> 0.0   , (s) -> 1.0   , 1e-8, s)
s_P_match_1, s_mom, s_er = markov_transition_moment_matching_serial(  (s) -> 0.0   , (s) -> 1.0   , 1e-8, s)
@test true == all(s_P_match_1 .== p_P_match_1) == all(s_mom .== p_mom) == all(s_er[isfinite.(s_er)] .== p_er[isfinite.(p_er)])
println("Single-var serial/parallel OK")

# p_P_match_2, p_mom, p_er = markov_transition_moment_matching_parallel((s) -> [s...], (s) -> [1.0 0.0; 0.0 1.0], 1e-8, ss...)
p_P_match_2, p_mom, p_er = markov_transition_moment_matching_parallel((s) -> [s...], (s) -> Matrix{Float64}(LinearAlgebra.I,2,2), 1e-8, ss...)
s_P_match_2, s_mom, s_er = markov_transition_moment_matching_serial(  (s) -> [s...], (s) -> Matrix{Float64}(LinearAlgebra.I,2,2), 1e-8, ss...)
@test true == all(s_P_match_2 .≈ p_P_match_2) == all(s_mom .≈ p_mom) == all(s_er[isfinite.(s_er)] .== p_er[isfinite.(p_er)])
println("Multi-var serial/parallel OK")

# share of states that get moment-matching
share = sum(s_mom .> 0) / length(s_mom) * 100
@show "$share percent of states get moment-matching"
@test share .> 0.5

# without matching
P_nomatch_2 = markov_transition((s) -> [s...], (s) -> Matrix{Float64}(I,2,2), 1e-8, ss...)
number_diff = sum(s_P_match_2 .!= P_nomatch_2)
sd_diff = sum((s_P_match_2 .- P_nomatch_2).^2)/prod(size(s_P_match_2))
abs_diff = sum(abs.(s_P_match_2 .- P_nomatch_2))/prod(size(s_P_match_2))
sharediff = number_diff / prod(size(s_P_match_2))
nm = norm(vec(s_P_match_2 .- P_nomatch_2), Inf)
@show "max difference between matching & no-matching is $nm"
@show "$sharediff percent of probabilities change"

# benchmarking
println("Benchmarking parallel...")
@show @benchmark markov_transition_moment_matching_parallel(  (s) -> [s...], (s) -> Matrix{Float64}(LinearAlgebra.I,2,2), 1e-8, ss...)
println("Benchmarking serial...")
@show @benchmark markov_transition_moment_matching_serial(    (s) -> [s...], (s) -> Matrix{Float64}(I,2,2), 1e-8, ss...)

rmprocs(pids)

# using Plots
# gr()
#
# heatmap(s, s, reshape(P_match_2[10,:]  , length.(ss)...))
# heatmap(s, s, reshape(P_nomatch_2[10,:], length.(ss)...))
# heatmap(s, s, reshape(P_match_2[20,:] .- P_nomatch_2[20,:], length.(ss)...))
# heatmap(s, s, reshape(mom.>0, length.(ss)...))
# histogram(vec(clamp.(P_match_2 .- P_nomatch_2, -0.05, 0.05)) )


text:
using Test
using NQCDynamics
using OrdinaryDiffEq: Tsit5
using Random: seed!

@testset "Ehrenfest" begin
    sim = RingPolymerSimulation{Ehrenfest}(Atoms(2000), TullyModelOne(), 10; temperature=1e-3)
    u = DynamicsVariables(sim, fill(10/2000, size(sim)), fill(-5, size(sim)) .+ randn(size(sim)), SingleState(1))
    dt = 0.1

    sol = run_trajectory(u, (0, 2000.0), sim; algorithm=DynamicsMethods.IntegrationAlgorithms.BCBwithTsit5(), saveat=0:10:2000, dt=dt)
    sol1 = run_trajectory(u, (0, 2000.0), sim; algorithm=Tsit5(), saveat=0:10:2000, abstol=1e-6, reltol=1e-6)

    @test sol.u ≈ sol1.u rtol=1e-3
end

@testset "FSSH" begin
    sim = RingPolymerSimulation{FSSH}(Atoms(2000), TullyModelOne(), 10; temperature=1e-3)
    u = DynamicsVariables(sim, fill(20/2000, size(sim)), fill(-5, size(sim)) .+ randn(size(sim)), SingleState(1))
    dt = 0.1

    seed!(1)
    sol = run_trajectory(u, (0, 2000.0), sim; algorithm=DynamicsMethods.IntegrationAlgorithms.BCBwithTsit5(), saveat=0:10:2000, dt=dt)
    seed!(1)
    sol1 = run_trajectory(u, (0, 2000.0), sim; algorithm=Tsit5(), saveat=0:10:2000, dt=dt, adaptive=false)

    @test sol.u ≈ sol1.u rtol=1e-3
end


text:
using Test
using MPI

MPI.Init()

comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)
root = 0

# Defining this to make ones work for Char
Base.one(::Type{Char}) = '\01'

for typ in Base.uniontypes(MPI.MPIDatatype)

    A = ones(typ, mod(rank,2) + 1)
    counts = Cint[mod(i,2) + 1 for i in 0:(size-1)]
    B = MPI.Gatherv(A, counts, root, comm)
    if rank == root
        @test B == ones(typ, 3 * div(size,2) + mod(size,2))
    end

    # Test passing the output buffer
    A = ones(typ, mod(rank,2) + 1)
    counts = Cint[mod(i,2) + 1 for i in 0:(size-1)]
    B = ones(typ, sum(counts)-1)
    if rank == root
        @test_throws AssertionError MPI.Gatherv!(A, B, counts, root, comm)
    end

    # Test assertion when output size is too small
    A = ones(typ, mod(rank,2) + 1)
    counts = Cint[mod(i,2) + 1 for i in 0:(size-1)]
    B = ones(typ, sum(counts))
    MPI.Gatherv!(A, B, counts, root, comm)
    if rank == root
        @test B == ones(typ, 3 * div(size,2) + mod(size,2))
    end

    # Test explicit MPI_IN_PLACE
    A = ones(typ, mod(rank,2) + 1)
    counts = Cint[mod(i,2) + 1 for i in 0:(size-1)]
    B = ones(typ, sum(counts))
    start = (cumsum(counts)-counts.+1)[MPI.Comm_rank(comm)+1]
    len   = counts[MPI.Comm_rank(comm)+1]
    B[start:(start+len-1)] .= A
    MPI.Gatherv_in_place!(B, counts, root, comm)
    if rank == root
        @test B == ones(typ, 3 * div(size,2) + mod(size,2))
    end
end

MPI.Finalize()


text:
##
# layer 3 access
# read data as records which are named tuple representations of the schema

##
# Row cursor iterates through row numbers of a column
mutable struct RowCursor
    par::ParFile

    rows::UnitRange{Int64}                      # rows to scan over
    row::Int64                                  # current row

    rowgroups::Vector{RowGroup}                 # row groups in range
    rg::Union{Int,Nothing}                      # current row group
    rgrange::Union{UnitRange{Int64},Nothing}    # current rowrange

    function RowCursor(par::ParFile, rows::UnitRange{Int64}, col::Vector{String}, row::Int64=first(rows))
        rgs = rowgroups(par, col, rows)
        cursor = new(par, rows, row, rgs, nothing, nothing)
        setrow(cursor, row)
        cursor
    end
end

function setrow(cursor::RowCursor, row::Int64)
    cursor.row = row
    (cursor.rgrange !== nothing) && (cursor.row in cursor.rgrange) && return
    startrow = Int64(1)
    rowgroups = cursor.rowgroups
    for rowgroup_idx in 1:length(rowgroups)
        rowgroup_row_range = startrow:(startrow + rowgroups[rowgroup_idx].num_rows)
        if row in rowgroup_row_range
            cursor.row = row
            cursor.rg = rowgroup_idx
            cursor.rgrange = rowgroup_row_range
            return
        else
            startrow = last(rowgroup_row_range) + 1
        end
    end
    throw(BoundsError(par.path, row))
end

rowgroup_offset(cursor::RowCursor) = cursor.row - first(cursor.rgrange)

function _start(cursor::RowCursor)
    row = first(cursor.rows)
    (cursor.row === row) || setrow(cursor, row)
    row
end
_done(cursor::RowCursor, row::Int64) = (row > last(cursor.rows))
function _next(cursor::RowCursor, row::Int64)
    setrow(cursor, row)
    row, (row+1)
end

function Base.iterate(cursor::RowCursor, state)
    _done(cursor, state) && return nothing
    return _next(cursor, state)
end

function Base.iterate(cursor::RowCursor)
    r = iterate(x, _start(x))
    return r
end

##
# Column cursor iterates through all values of the column, including null values.
# Each iteration returns the value (as a Union{T,Nothing}), definition level, and repetition level for each value.
# Row can be deduced from repetition level.
mutable struct ColCursor{T}
    row::RowCursor
    colname::Vector{String}
    required_at::Vector{Bool}
    repeated_at::Vector{Bool}
    logical_converter_fn::Function
    maxdefn::Int32

    colchunks::Union{Vector{ColumnChunk},Nothing}
    cc::Union{Int,Nothing}
    ccrange::Union{UnitRange{Int64},Nothing}

    vals::Vector{T}
    valpos::Int64

    defn_levels::Vector{Int32}
    repn_levels::Vector{Int32}
    levelpos::Int64
    levelrange::UnitRange{Int}

    function ColCursor{T}(row::RowCursor, colname::Vector{String}) where T
        sch = schema(row.par)

        len_colname_parts = length(colname)
        required_at = Array{Bool}(undef, len_colname_parts)
        repeated_at = Array{Bool}(undef, len_colname_parts)
        for idx in 1:len_colname_parts
            partname = colname[1:idx]
            required_at[idx] = isrequired(sch, partname)         # determine whether field is optional and repeated
            repeated_at[idx] = isrepeated(sch, partname)
        end

        logical_converter_fn = logical_converter(sch, colname)
        maxdefn = max_definition_level(sch, colname)
        new{T}(row, colname, required_at, repeated_at, logical_converter_fn, maxdefn, nothing, nothing, nothing)
    end
end

function ColCursor(par::ParFile, colname::Vector{String}; rows::UnitRange=1:nrows(par), row::Signed=first(rows))
    rowcursor = RowCursor(par, UnitRange{Int64}(rows), colname, Int64(row))

    rowgroup = rowcursor.rowgroups[rowcursor.rg]
    colchunks = columns(par, rowgroup, colname)
    parquet_coltype = coltype(colchunks[1])
    T = PLAIN_JTYPES[parquet_coltype+1]
    if (parquet_coltype == _Type.BYTE_ARRAY) || (parquet_coltype == _Type.FIXED_LEN_BYTE_ARRAY)
        T = Vector{T}
    end

    cursor = ColCursor{T}(rowcursor, colname)
    setrow(cursor, Int64(row))
    cursor
end

eltype(cursor::ColCursor{T}) where {T} = NamedTuple{(:value, :defn_level, :repn_level),Tuple{Union{Nothing,T},Int64,Int64}}
length(cursor::ColCursor) = length(cursor.row.rows)

function setrow(cursor::ColCursor{T}, row::Int64) where {T}
    par = cursor.row.par
    rg = cursor.row.rowgroups[cursor.row.rg]
    ccincr = (row - cursor.row.row) == 1 # whether this is just an increment within the column chunk
    setrow(cursor.row, row) # set the row cursor
    if cursor.colchunks === nothing
        cursor.colchunks = columns(par, rg, cursor.colname)
    end

    # check if cursor is done
    if _done(cursor.row, row)
        cursor.cc = length(cursor.colchunks) + 1
        #cursor.ccrange = row:(row-1)
        #cursor.vals = Vector{T}()
        #cursor.repn_levels = cursor.defn_levels = Int32[]
        cursor.valpos = cursor.levelpos = 0
        cursor.levelrange = 0:-1
        return
    end

    # find the column chunk with the row
    if (cursor.ccrange === nothing) || !(row in cursor.ccrange)
        offset = rowgroup_offset(cursor.row) # the offset of row from beginning of current rowgroup
        colchunks = (cursor.colchunks)::Vector{ColumnChunk}

        startrow = row - offset
        for cc in 1:length(colchunks)
            vals, defn_levels, repn_levels = values(par, colchunks[cc], T)
            if isempty(repn_levels)
                nrowscc = length(vals) # number of values is number of rows
            else
                nrowscc = 0
                for i in 1:length(repn_levels)
                    (repn_levels[i] !== Int32(0)) && (nrowscc += 1)
                end
            end
            ccrange = startrow:(startrow + nrowscc)

            if row in ccrange
                cursor.cc = cc
                cursor.ccrange = ccrange
                cursor.vals = vals
                cursor.defn_levels = defn_levels
                cursor.repn_levels = repn_levels
                cursor.valpos = cursor.levelpos = 0
                ccincr = false
                break
            else
                startrow = last(ccrange) + 1
            end
        end
    end

    if cursor.ccrange === nothing
        # we did not find the row in this column
        cursor.valpos = cursor.levelpos = 0
        cursor.levelrange = 0:-1
    else
        # find the starting positions for values and levels
        ccrange = (cursor.ccrange)::UnitRange{Int64}
        defn_levels = cursor.defn_levels
        repn_levels = cursor.repn_levels
        levelpos = valpos = Int64(0)

        # compute the level and value pos for row
        if isempty(repn_levels)
            # no repetitions, so each entry corresponds to one full row
            levelpos = row - first(ccrange) + 1
            levelrange = levelpos:levelpos
        else
            # multiple entries may constitute one row
            idx = first(ccrange)
            levelpos = Int64(findfirst(x->x===Int32(0), repn_levels)) # NOTE: can start from cursor.levelpos to optimize, but that will prevent using setrow to go backwards
            while idx < row
                levelpos = Int64(findnext(x->x===Int32(0), repn_levels, levelpos+1))
                idx += 1
            end
            levelend = Int64(max(findnext(x->x===Int32(0), repn_levels, levelpos+1)-1, length(repn_levels)))
            levelrange = levelpos:levelend
        end

        # compute the val pos for row
        if isempty(defn_levels)
            # all entries are required, so there must be a corresponding value
            valpos = levelpos
        else
            maxdefn = cursor.maxdefn
            if ccincr
                valpos = cursor.valpos
            else
                valpos = 1
                for idx in 1:(levelpos-1)
                    (defn_levels[idx] === maxdefn) && (valpos += 1)
                end
            end
        end

        cursor.levelpos = levelpos
        cursor.levelrange = levelrange
        cursor.valpos = valpos
    end
    nothing
end

function _start(cursor::ColCursor)
    row = _start(cursor.row)
    setrow(cursor, row)
    row, cursor.levelpos
end
function _done(cursor::ColCursor, rowandlevel::Tuple{Int64,Int64})
    row, levelpos = rowandlevel
    (levelpos > last(cursor.levelrange)) || _done(cursor.row, row)
end
function _next(cursor::ColCursor{T}, rowandlevel::Tuple{Int64,Int64}) where {T}
    # find values for current row and level in row
    row, levelpos = rowandlevel
    (levelpos == cursor.levelpos) || throw(InvalidStateException("Invalid column cursor state", :levelpos))

    maxdefn = cursor.maxdefn
    defn_level = isempty(cursor.defn_levels) ? maxdefn : cursor.defn_levels[cursor.valpos]
    repn_level = isempty(cursor.repn_levels) ? 0 : cursor.repn_levels[cursor.valpos]
    cursor.levelpos += 1
    if defn_level == maxdefn
        val = (cursor.vals[cursor.valpos])::T
        cursor.valpos += 1
    else
        val = nothing
    end

    # advance row
    if cursor.levelpos > last(cursor.levelrange)
        row += 1
        setrow(cursor, row)
    end

    NamedTuple{(:value, :defn_level, :repn_level),Tuple{Union{Nothing,T},Int64,Int64}}((val, defn_level, repn_level)), (row, cursor.levelpos)
end

function Base.iterate(cursor::ColCursor{T}, state) where {T}
    _done(cursor, state) && return nothing
    return _next(cursor, state)
end

function Base.iterate(cursor::ColCursor)
    r = iterate(cursor, _start(cursor))
    return r
end

##

mutable struct BatchedColumnsCursor{T}
    par::ParFile
    colnames::Vector{Vector{String}}
    colcursors::Vector{ColCursor}
    colstates::Vector{Tuple{Int64,Int64}}
    rowgroupid::Int
    row::Int64
end

function BatchedColumnsCursor(par::ParFile)
    sch = schema(par)

    # supports only non nested columns as of now
    if !all(num_children(schemaelem) == 0 for schemaelem in sch.schema[2:end])
        error("nested schemas are not supported with BatchedColumnsCursor yet")
    end

    colcursors = [ColCursor(par, colname) for colname in colnames(par)]
    rectype = ntcolstype(sch, sch.schema[1])
    BatchedColumnsCursor{rectype}(par, colnames(par), colcursors, Array{Tuple{Int64,Int64}}(undef, length(colcursors)), 1, 1)
end

eltype(cursor::BatchedColumnsCursor{T}) where {T} = T
length(cursor::BatchedColumnsCursor) = length(rowgroups(cursor.par))

function colcursor_values(colcursor::ColCursor{T}, ::Type{Vector{Union{Missing,V}}}) where {T,V}
    defn_levels = colcursor.defn_levels
    vals = colcursor.vals
    val_idx = 0
    logical_converter_fn = colcursor.logical_converter_fn
    Union{Missing,V}[defn_levels[idx] === Int32(1) ? logical_converter_fn(vals[val_idx+=1]) : missing for idx in 1:length(defn_levels)]
end

function colcursor_values(colcursor::ColCursor{T}, ::Type{Vector{Union{Missing,T}}}) where {T}
    defn_levels = colcursor.defn_levels
    vals = colcursor.vals
    val_idx = 0
    Union{Missing,T}[defn_levels[idx] === Int32(1) ? vals[val_idx+=1] : missing for idx in 1:length(defn_levels)]
end

function colcursor_values(colcursor::ColCursor{T}, ::Type{Vector{V}}) where {T,V}
    defn_levels = colcursor.defn_levels
    vals = colcursor.vals
    logical_converter_fn = colcursor.logical_converter_fn
    V[logical_converter_fn(v) for v in vals]
end

colcursor_values(colcursor::ColCursor{T}, ::Type{Vector{T}}) where {T} = colcursor.vals

function Base.iterate(cursor::BatchedColumnsCursor{T}, rowgroupid) where {T}
    (rowgroupid > length(cursor)) && (return nothing)

    colcursors = cursor.colcursors
    for colcursor in colcursors
        setrow(colcursor, cursor.row)
    end
    coltypes = T.types
    colvals = [colcursor_values(colcursor,coltype) for (colcursor,coltype) in zip(colcursors,coltypes)]

    cursor.row += (rowgroups(cursor.par)[cursor.rowgroupid]).num_rows
    cursor.rowgroupid += 1
    T(colvals), cursor.rowgroupid
end

function Base.iterate(cursor::BatchedColumnsCursor{T}) where {T}
    cursor.row = 1
    cursor.colstates = [_start(colcursor) for colcursor in cursor.colcursors]
    iterate(cursor, cursor.rowgroupid)
end

##

mutable struct RecordCursor{T}
    par::ParFile
    colnames::Vector{Vector{String}}
    colcursors::Vector{ColCursor}
    colstates::Vector{Tuple{Int64,Int64}}
    rows::UnitRange{Int64}                      # rows to scan over
    row::Int64                                  # current row
end

function RecordCursor(par::ParFile; rows::UnitRange=1:nrows(par), colnames::Vector{Vector{String}}=colnames(par), row::Signed=first(rows))
    colcursors = [ColCursor(par, colname; rows=rows, row=row) for colname in colnames]
    sch = schema(par)
    rectype = ntelemtype(sch, sch.schema[1])
    RecordCursor{rectype}(par, colnames, colcursors, Array{Tuple{Int64,Int64}}(undef, length(colcursors)), rows, row)
end

eltype(cursor::RecordCursor{T}) where {T} = T
length(cursor::RecordCursor) = length(cursor.rows)

function Base.iterate(cursor::RecordCursor{T}, row) where {T}
    (row > last(cursor.rows)) && (return nothing)

    states = cursor.colstates
    cursors = cursor.colcursors

    colvals = Dict{Symbol,Any}()
    col_repeat_state = Dict{Tuple{Int,Int},Int}()
    for colid in 1:length(states)  # for each column
        colcursor = cursors[colid]
        colstate = states[colid]
        states[colid] = update_record(cursor.par, colvals, colid, colcursor, colstate, col_repeat_state)
    end
    cursor.row += 1
    _nt(colvals, T), cursor.row
end

function Base.iterate(cursor::RecordCursor{T}) where {T}
    cursor.row = first(cursor.rows)
    cursor.colstates = [_start(colcursor) for colcursor in cursor.colcursors]
    iterate(cursor, cursor.row)
end

function _val_or_missing(dict::Dict{Symbol,Any}, k::Symbol, ::Type{T}) where {T}
    v = get(dict, k, missing)
    if isa(v, Vector)
        elt = eltype(v)
        if Dict{Symbol,Any} <: elt
            nonmissing_elt = nonmissingtype(eltype(T))
            v = [el === missing ? el : _nt(el,nonmissing_elt) for el in v]
        end
    end
    (isa(v, Dict{Symbol,Any}) ? _nt(v, nonmissingtype(T)) : v)::T
end

@generated function _nt(dict::Dict{Symbol,Any}, ::Type{T}) where {T}
    names = fieldnames(T)
    strnames = ["$n" for n in names]
    quote
        return T(($([:(_val_or_missing(dict,Symbol($(strnames[i])),$(fieldtype(T,i)))) for i in 1:length(names)]...),))
    end
end

default_init(::Type{Vector{T}}) where {T} = Vector{T}()
default_init(::Type{Dict{Symbol,Any}}) = Dict{Symbol,Any}()
default_init(::Type{T}) where {T} = ccall(:jl_new_struct_uninit, Any, (Any,), T)::T

function update_record(par::ParFile, row::Dict{Symbol,Any}, colid::Int, colcursor::ColCursor, colcursor_state::Tuple{Int64,Int64}, col_repeat_state::Dict{Tuple{Int,Int},Int})
    if !_done(colcursor, colcursor_state)
        colval, colcursor_state = _next(colcursor, colcursor_state)                                                         # for each value, defn level, repn level in column
        update_record(par, row, colid, colcursor, colval.value, colval.defn_level, colval.repn_level, col_repeat_state)    # update record
    end
    colcursor_state # return new colcursor state
end

function update_record(par::ParFile, row::Dict{Symbol,Any}, colid::Int, colcursor::ColCursor, val, defn_level::Int64, repn_level::Int64, col_repeat_state::Dict{Tuple{Int,Int},Int})
    nameparts = colcursor.colname
    required_at = colcursor.required_at
    repeated_at = colcursor.repeated_at
    logical_converter_fn = colcursor.logical_converter_fn

    lparts = length(nameparts)
    sch = par.schema
    F = row  # the current field corresponding to the level in nameparts
    Fdefn = 0
    Frepn = 0

    # for each name part of colname (a field)
    for idx in 1:lparts
        colname = view(nameparts, 1:idx)
        #@debug("updating part $colname of $nameparts isnull:$(val === nothing), def:$(defn_level), rep:$(repn_level)")
        leaf = nameparts[idx]
        symleaf = Symbol(leaf)

        required = required_at[idx]                 # determine whether field is optional and repeated
        repeated = repeated_at[idx]
        required || (Fdefn += 1)                    # if field is optional, increment defn level
        repeated && (Frepn += 1)                    # if field can repeat, increment repn level

        defined = ((val === nothing) || (idx < lparts)) ? haskey(F, symleaf) : false
        mustdefine = defn_level >= Fdefn
        mustrepeat = repeated && (repn_level <= Frepn)
        repidx = 0
        if mustrepeat
            repkey = (colid, idx)
            repidx = get(col_repeat_state, repkey, 0)
            repidx += 1
            col_repeat_state[repkey] = repidx
        end
        nreps = (defined && isa(F[symleaf], Vector)) ? length(F[symleaf]) : 0

        #@debug("repeat:$mustrepeat, nreps:$nreps, repidx:$repidx, defined:$defined, mustdefine:$mustdefine")
        if mustrepeat && (nreps < repidx)
            if !defined && mustdefine
                Vtyp = elemtype(sch, colname)
                Vrep = F[symleaf] = default_init(Vtyp)
            else
                Vrep = F[symleaf]
            end
            if length(Vrep) < repidx
                resize!(Vrep, repidx)
                if !isbits(eltype(Vrep))
                    Vrep[repidx] = default_init(eltype(Vrep))
                end
            end
            F = Vrep[repidx]
        elseif !defined && mustdefine
            if idx == length(nameparts)
                V = logical_converter_fn(val)
            else
                Vtyp = elemtype(sch, colname)
                V = default_init(Vtyp)
            end
            F[symleaf] = V
            F = V
        else
            F = get(F, symleaf, nothing)
        end
    end
    nothing
end


text:
# This file is part of the TaylorSeries.jl Julia package, MIT license
#
# Luis Benet & David P. Sanders
# UNAM
#
# Handles Taylor series of arbitrary but finite order

"""
    TaylorSeries

A Julia package for Taylor expansions in one or more independent variables.

The basic constructors are [`Taylor1`](@ref) and [`TaylorN`](@ref);
see also [`HomogeneousPolynomial`](@ref).

"""
module TaylorSeries

using InteractiveUtils: subtypes
using SparseArrays: SparseMatrixCSC
using Markdown
import LinearAlgebra: norm, mul!
if VERSION ≥ v"1.0.0"
    # export gradient
else
    import LinearAlgebra: gradient
end

import Base: ==, +, -, *, /, ^

import Base: iterate, size, eachindex, firstindex, lastindex,
    eltype, length, getindex, setindex!

import Base: zero, one, zeros, ones, isinf, isnan, iszero,
    convert, promote_rule, promote, show,
    real, imag, conj, adjoint,
    rem, mod, mod2pi, abs, abs2,
    sqrt, exp, log, sin, cos, tan,
    asin, acos, atan, sinh, cosh, tanh,
    power_by_squaring,
    rtoldefault, isfinite, isapprox, rad2deg, deg2rad

export Taylor1, TaylorN, HomogeneousPolynomial, AbstractSeries

export getcoeff, derivative, integrate, differentiate,
    evaluate, evaluate!, inverse,
    show_params_TaylorN, show_monomials, displayBigO, use_show_default,
    get_order, get_numvars,
    set_variables, get_variables,
    get_variable_names, get_variable_symbols,
    # jacobian, hessian, jacobian!, hessian!,
    ∇, taylor_expand, update!, constant_term

include("parameters.jl")
include("hash_tables.jl")
include("constructors.jl")
include("conversion.jl")
include("auxiliary.jl")
include("arithmetic.jl")
include("power.jl")
include("functions.jl")
include("other_functions.jl")
include("evaluate.jl")
include("calculus.jl")
include("dictmutfunct.jl")
include("printing.jl")

end # module


text:
using Base.Iterators: product
using Dates
using Documenter: doctest
using Infinity: Infinite, InfExtendedReal, InfExtendedTime, InfMinusInfError, ∞
using Intervals
using Intervals: isfinite
using Serialization: deserialize
using Test
using TimeZones

const BOUND_PERMUTATIONS = product((Closed, Open), (Closed, Open))

include("test_utils.jl")

@testset "Intervals" begin
    include("inclusivity.jl")
    include("endpoint.jl")
    include("interval.jl")
    include("anchoredinterval.jl")
    include("comparisons.jl")
    include("plotting.jl")

    # Note: The output of the doctests currently requires a newer version of Julia
    # https://github.com/JuliaLang/julia/pull/34387
    # The doctests fail on x86, so only run them on 64-bit hardware
    if v"1.6" <= VERSION < v"1.7" && Sys.WORD_SIZE == 64
        doctest(Intervals)
    else
        @warn "Skipping doctests"
    end
end


text:
function get_pendulum(;
    timestep=0.01,
    gravity=-9.81,
    mass=1.0,
    len=1.0,
    spring=0.0,
    damper=0.0,
    spring_offset=szeros(1),
    orientation_offset=one(Quaternion),
    T=Float64)

    # Parameters
    joint_axis = [1.0; 0.0; 0.0]
    width, depth = 0.1, 0.1
    parent_vertex = [0.0; 0.0; 1.1 * len] # joint connection point
    child_vertex = [0.0; 0.0; 0.5 * len] # joint connection point

    # Links
    origin = Origin{T}()
    body = Box(width, depth, len, mass, name=:pendulum)

    # Constraints
    joint = JointConstraint(Revolute(origin, body, joint_axis;
        parent_vertex,
        child_vertex,
        spring,
        damper,
        rot_spring_offset=spring_offset,
        orientation_offset), name=:joint)

    bodies = [body]
    joints = [joint]

    mech = Mechanism(origin, bodies, joints;
        gravity,
        timestep)

    return mech
end

function initialize_pendulum!(mechanism::Mechanism;
    angle=0.7,
    angular_velocity=0.0) where T
    joint = mechanism.joints[1]
    set_minimal_coordinates_velocities!(mechanism, joint;
        xmin=[angle, angular_velocity])
end

function get_npendulum(;
    timestep=0.01,
    gravity=-9.81,
    mass=1.0,
    len=1.0,
    spring=0.0,
    damper=0.0,
    num_bodies=5,
    basetype=:Revolute,
    joint_type=:Revolute,
    T=Float64)

    # Parameters
    ex = [1.0; 0.0; 0.0]
    r = 0.05
    vert11 = [0.0; 0.0; len / 2.0]
    vert12 = -vert11

    # Links
    origin = Origin{T}()
    bodies = [Box(r, r, len, mass, color=RGBA(1.0, 0.0, 0.0)) for i = 1:num_bodies]

    # Constraints
    jointb1 = JointConstraint(Prototype(basetype, origin, bodies[1], ex;
        child_vertex=vert11,
        spring,
        damper))
    if num_bodies > 1
        joints = [JointConstraint(Prototype(joint_type, bodies[i - 1], bodies[i], ex;
            parent_vertex=vert12,
            child_vertex=vert11,
            spring,
            damper)) for i = 2:num_bodies]
        joints = [jointb1; joints]
    else
        joints = [jointb1]
    end

    mech = Mechanism(origin, bodies, joints;
        gravity,
        timestep)
    return mech
end

function initialize_npendulum!(mechanism::Mechanism{T};
    base_angle=π / 4.0,
    base_angular_velocity=[0.0, 0.0, 0.0],
    relative_linear_velocity=[0.0, 0.0, 0.0],
    relative_angular_velocity=[0.0, 0.0, 0.0]) where T

    pbody = mechanism.bodies[1]
    joint = mechanism.joints[1]
    vert11 = joint.translational.vertices[2]
    vert12 = -vert11

    # set position and velocities
    set_maximal_configurations!(mechanism.origin, pbody,
        child_vertex=vert11,
        Δq=RotX(base_angle))
    set_maximal_velocities!(pbody,
        ω=base_angular_velocity)

    previd = pbody.id
    for (i, body) in enumerate(Iterators.drop(mechanism.bodies, 1))
        set_maximal_configurations!(get_body(mechanism, previd), body,
            parent_vertex=vert12,
            child_vertex=vert11)
        set_maximal_velocities!(get_body(mechanism, previd), body,
            parent_vertex=vert12,
            child_vertex=vert11,
            Δv=relative_linear_velocity, Δω=1 / i * relative_angular_velocity)
        previd = body.id
    end
end


text:
function numerical_diff(f::Union{Function,Core.Function}, x::Variable, eps::Float64=1e-4)
    x0 = Variable(x.data .- eps)
    x1 = Variable(x.data .+ eps)
    y0::Variable = f(x0)
    y1::Variable = f(x1)
    return (y1.data - y0.data) / (2 * eps)
end


text:
module TestLossFunctions

using FactCheck
importall GradientBoost.LossFunctions

type DummyLossFunction <: LossFunction; end

y_examples = {
  [1,1],
  [1,0],
  [0,0],
  [-1,1],
  [-1,-1]
}
y_pred_examples = {
  [1,1],
  [1,1],
  [1,1],
  [1,1],
  [1,1]
}
bern_y_examples = {
  [1,1],
  [1,0],
  [0,1],
  [0,0],
  [0,0]
}

function test_loss(lf::LossFunction,
  y_examples, y_pred_examples, expected)

  for i = 1:size(y_examples, 1)
    @fact loss(lf, y_examples[i], y_pred_examples[i]) => roughly(expected[i])
  end
end
function test_negative_gradient(lf::LossFunction,
  y_examples, y_pred_examples, expected)

  for i = 1:size(y_examples, 1)
    @fact negative_gradient(
      lf, y_examples[i], y_pred_examples[i]
    ) => roughly(expected[i])
  end
end
function test_minimizing_scalar(lf::LossFunction,
  y_examples, expected)

  for i = 1:size(y_examples, 1)
    @fact minimizing_scalar(lf, y_examples[i]) => roughly(expected[i])
  end
end

facts("Loss functions") do
  context("not implemented functions throw an error") do
    emp_vec = Array[]
    dlf = DummyLossFunction()

    @fact_throws loss(dlf, emp_vec, emp_vec)
    @fact_throws negative_gradient(dlf, emp_vec, emp_vec)
    @fact_throws minimizing_scalar(dlf, emp_vec)
  end

  context("LeastSquares loss works") do
    lf = LeastSquares()
    expected = { 0.0, 0.5, 1.0, 2.0, 4.0 }
    test_loss(lf, y_examples, y_pred_examples, expected)
  end
  context("LeastSquares negative_gradient works") do
    lf = LeastSquares()
    expected = { [0,0], [0,-1], [-1,-1], [-2,0], [-2,-2] }
    test_negative_gradient(lf, y_examples, y_pred_examples, expected)
  end
  context("LeastSquares minimizing_scalar works") do
    lf = LeastSquares()
    expected = { 1.0, 0.5, 0.0, 0.0, -1.0 }
    test_minimizing_scalar(lf, y_examples, expected)
  end

  context("LeastAbsoluteDeviation loss works") do
    lf = LeastAbsoluteDeviation()
    expected = { 0.0, 0.5, 1.0, 1.0, 2.0 }
    test_loss(lf, y_examples, y_pred_examples, expected)
  end
  context("LeastAbsoluteDeviation negative_gradient works") do
    lf = LeastAbsoluteDeviation()
    expected = { [0,0], [0,-1], [-1,-1], [-1,0], [-1,-1] }
    test_negative_gradient(lf, y_examples, y_pred_examples, expected)
  end
  context("LeastAbsoluteDeviation minimizing_scalar works") do
    lf = LeastAbsoluteDeviation()
    expected = { 1.0, 0.5, 0.0, 0.0, -1.0 }
    test_minimizing_scalar(lf, y_examples, expected)
  end

  context("BinomialDeviance loss works") do
    lf = BinomialDeviance()
    expected = {
      0.626523, 1.626523, 1.626523, 2.626523, 2.626523,
    }
    test_loss(lf, bern_y_examples, y_pred_examples, expected)
  end
  context("BinomialDeviance negative_gradient works") do
    lf = BinomialDeviance()
    expected = {
      [0.268941, 0.268941],
      [0.268941, -0.731059],
      [-0.731059, 0.268941],
      [-0.731059, -0.731059],
      [-0.731059, -0.731059]
    }
    test_negative_gradient(lf, bern_y_examples, y_pred_examples, expected)
  end
  context("BinomialDeviance minimizing_scalar works") do
    lf = BinomialDeviance()
    expected = { Inf, 0.0, 0.0, -Inf, -Inf}
    test_minimizing_scalar(lf, bern_y_examples, expected)
  end
end

end # module


text:
export Gamma

"""
Description:

    A gamma node with shape-rate parameterization:

    f(out,a,b) = Gam(out|a,b) = 1/Γ(a) b^a out^{a - 1} exp(-b out)

Interfaces:

    1. out
    2. a (shape)
    3. b (rate)

Construction:

    Gamma(out, a, b, id=:some_id)
"""
mutable struct Gamma <: SoftFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function Gamma(out, a, b; id=generateId(Gamma))
        @ensureVariables(out, a, b)
        self = new(id, Array{Interface}(undef, 3), Dict{Symbol,Interface}())
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)
        self.i[:a] = self.interfaces[2] = associate!(Interface(self), a)
        self.i[:b] = self.interfaces[3] = associate!(Interface(self), b)

        return self
    end
end

slug(::Type{Gamma}) = "Gam"

format(dist::ProbabilityDistribution{Univariate, Gamma}) = "$(slug(Gamma))(a=$(format(dist.params[:a])), b=$(format(dist.params[:b])))"

ProbabilityDistribution(::Type{Univariate}, ::Type{Gamma}; a=1.0, b=1.0) = ProbabilityDistribution{Univariate, Gamma}(Dict(:a=>a, :b=>b))
ProbabilityDistribution(::Type{Gamma}; a=1.0, b=1.0) = ProbabilityDistribution{Univariate, Gamma}(Dict(:a=>a, :b=>b))

dims(dist::ProbabilityDistribution{Univariate, Gamma}) = ()

vague(::Type{Gamma}) = ProbabilityDistribution(Univariate, Gamma, a=1.0, b=tiny) # Flat prior leads to more stable behaviour than Jeffrey's prior

unsafeMean(dist::ProbabilityDistribution{Univariate, Gamma}) = dist.params[:a]/dist.params[:b] # unsafe mean

unsafeLogMean(dist::ProbabilityDistribution{Univariate, Gamma}) = digamma(dist.params[:a]) - log(dist.params[:b])

# https://stats.stackexchange.com/questions/457357/what-is-the-expected-value-of-x-logx-of-the-gamma-distribution
unsafeMeanLogMean(dist::ProbabilityDistribution{Univariate, Gamma}) =  (gamma(dist.params[:a]+1)/(gamma(dist.params[:a])*dist.params[:b])) * (digamma(dist.params[:a]+1) - log(dist.params[:b]))

unsafeVar(dist::ProbabilityDistribution{Univariate, Gamma}) = dist.params[:a]/dist.params[:b]^2 # unsafe variance

logPdf(dist::ProbabilityDistribution{Univariate, Gamma}, x) = dist.params[:a]*log(dist.params[:b]) - labsgamma(dist.params[:a]) + (dist.params[:a]-1)*log(x) - dist.params[:b]*x

isProper(dist::ProbabilityDistribution{Univariate, Gamma}) = (dist.params[:a] >= tiny) && (dist.params[:b] >= tiny)

function prod!( x::ProbabilityDistribution{Univariate, Gamma},
                y::ProbabilityDistribution{Univariate, Gamma},
                z::ProbabilityDistribution{Univariate, Gamma}=ProbabilityDistribution(Univariate, Gamma, a=0.0, b=0.0))

    z.params[:a] = x.params[:a] + y.params[:a] - 1.0
    z.params[:b] = x.params[:b] + y.params[:b]

    return z
end

@symmetrical function prod!(x::ProbabilityDistribution{Univariate, Gamma},
                            y::ProbabilityDistribution{Univariate, PointMass},
                            z::ProbabilityDistribution{Univariate, PointMass}=ProbabilityDistribution(Univariate, PointMass, m=0.0))

    (y.params[:m] > 0.0) || error("PointMass location $(y.params[:m]) should be positive")
    z.params[:m] = y.params[:m]

    return z
end

sample(dist::ProbabilityDistribution{Univariate, Gamma}) = gammainvcdf(dist.params[:a], 1/dist.params[:b], rand())

naturalParams(dist::ProbabilityDistribution{Univariate, Gamma}) = [dist.params[:a]-1.0, -dist.params[:b]]

standardDistribution(V::Type{Univariate}, F::Type{Gamma}; η::Vector) = ProbabilityDistribution(V, F, a=η[1]+1.0, b=-η[2])

logNormalizer(::Type{Univariate}, ::Type{Gamma}; η::Vector) = loggamma(η[1]+1.0) - (η[1]+1.0)*log(-η[2])

logPdf(V::Type{Univariate}, F::Type{Gamma}, x::Number; η::Vector) =[log(x), x]'*η - logNormalizer(V, F, η=η)

# Entropy functional
function differentialEntropy(dist::ProbabilityDistribution{Univariate, Gamma})
    labsgamma(dist.params[:a]) -
    (dist.params[:a] - 1.0)*digamma(dist.params[:a]) -
    log(dist.params[:b]) +
    dist.params[:a]
end

# Average energy functional
function averageEnergy(::Type{Gamma}, marg_out::ProbabilityDistribution{Univariate}, marg_a::ProbabilityDistribution{Univariate, PointMass}, marg_b::ProbabilityDistribution{Univariate})
    labsgamma(marg_a.params[:m]) -
    marg_a.params[:m]*unsafeLogMean(marg_b) -
    (marg_a.params[:m] - 1.0)*unsafeLogMean(marg_out) +
    unsafeMean(marg_b)*unsafeMean(marg_out)
end

# By Stirling's approximation
function averageEnergy(::Type{Gamma}, marg_out::ProbabilityDistribution{Univariate}, marg_a::ProbabilityDistribution{Univariate}, marg_b::ProbabilityDistribution{Univariate})
    unsafeMeanLogMean(marg_a) - unsafeMean(marg_a) + 0.5*(log(2*pi)-unsafeLogMean(marg_a)) -
    unsafeMean(marg_a)*unsafeLogMean(marg_b) -
    (unsafeMean(marg_a) - 1.0)*unsafeLogMean(marg_out) +
    unsafeMean(marg_b)*unsafeMean(marg_out)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 310ce685-2661-4f32-bf14-91a4f4e569ce
# Explicit use of own environment instead of a local one for each notebook
begin
	using Pkg
	cd(joinpath(dirname(@__FILE__),".."))
    Pkg.activate(pwd())
    using NativeSVG
end

# ╔═╡ bbf4ab90-eb95-11ea-174e-c3624a0d92cd
using Random

# ╔═╡ 43627847-ce0d-4969-b6e1-2d59a54d1fd2
using PlutoUI

# ╔═╡ 56f41ca0-eb93-11ea-1ea6-11b0e8bb9a7d
md"""# Cellular Automaton

Port of [Think Complexity chapter 5](http://greenteapress.com/complexity2/html/index.html) by Allen Downey."""

# ╔═╡ dab709d0-eb93-11ea-050e-3bfa6a5e1836
md"Cellular Automaton = discrete space (cells) as input for a calculation in discrete time"

# ╔═╡ f4b601b0-eb93-11ea-0969-7967a1e85c8b
md"""## A Trivial Example

0 dimensional CA, inverting its cell at each timestep (2 state values only)

rule:"""

# ╔═╡ 14a02c30-eb94-11ea-2114-5102409c8ae5
function rule0dim(x::Bool)
    !x
end

# ╔═╡ 1a915240-eb94-11ea-087f-231ade62690d
md"time evolution:"

# ╔═╡ 2393fe0e-eb94-11ea-2858-b50503395d4a
function step0dim(x₀::Bool, steps::Int64)
    xs = [x₀]
    for i in 1:steps
        push!(xs, rule0dim(xs[end]))
    end
    xs
end

# ╔═╡ 2dde4510-eb94-11ea-212a-5da1e7733bf6
md"visualisation:"

# ╔═╡ 35cbd760-eb94-11ea-201a-6b58802154a1
let
	res = step0dim(false, 10)
	Drawing(width = 50, height = 300) do
		for (i, val) in enumerate(res)
			fill = if val "grey" else "lightgrey" end
			rect(x = 20, y = 5+i*20, width = 20, height = 20, fill = fill)
		end
	end
end

# ╔═╡ 758a0c50-eb94-11ea-11a2-e3e007b089a9
md"""## Wolfram's Experiment

1 dimensional CA with 2 state values, new value of a cell depends only on state of neighbouring cells.

rule can be expressed as an integer

|prev  |111|110|101|100|011|010|001|000|
|------|---|---|---|---|---|---|---|---|
|next  |b7 |b6 |b5 |b4 |b3 |b2 |b1 |b0 |
|rule50|0  |0  |1  |1  |0  |0  |1  |0  |

b can be converted to an integer

get a rule from an integer:"""

# ╔═╡ d8481030-eb94-11ea-1af4-db838adc37ed
function inttorule1dim(val::UInt8)
    digs = BitArray(digits(val, base=2))
    for i in length(digs):7
        push!(digs, false)
    end
    digs
end

# ╔═╡ e393e8b2-eb94-11ea-3d70-d31f7ee89420
md"Apply rule to a cell knowing its own previous state and the previous state of his left and right neighbour"

# ╔═╡ eed0a600-eb94-11ea-1863-0f33980508ba
function applyrule1dim(rule::BitArray{1}, bits::BitArray{1})
    val = 1 + bits[3] + 2*bits[2] + 4*bits[1]
    rule[val]
end

# ╔═╡ 0002a810-eb95-11ea-2ba8-bb2849bdec17
md"time evolution:"

# ╔═╡ 04139c70-eb95-11ea-1759-298193ce97b0
function step1dim(x₀::BitArray{1}, rule::BitArray{1}, steps::Int64)
    xs = [x₀]
    len = length(x₀)
    for i in 1:steps
        x = copy(x₀)
        for j in 2:len-1
            x[j] = applyrule1dim(rule, xs[end][j-1:j+1])
        end
        push!(xs, x)
    end
    xs
end

# ╔═╡ 10100040-eb95-11ea-3a6d-271e63301b17
md"initialisation:"

# ╔═╡ 1a508340-eb95-11ea-08f9-6f085748d7ff
res = let
	x₀ = falses(21)
	x₀[11] = true
	step1dim(x₀, inttorule1dim(UInt8(50)), 9)
end;

# ╔═╡ 5d2a01f0-eb95-11ea-3dcd-b794e1b0d566
md"visualisation:"

# ╔═╡ 60e8deb0-eb95-11ea-2bde-d9c259432318
function visualize1dim(res, dim)
    width = dim * (length(res[1]) + 1)
    height = dim * (length(res) + 1)
    Drawing(width = width, height = height) do
        for (i, arr) in enumerate(res)
            for (j, val) in enumerate(arr)
                fill = if val "grey" else "lightgrey" end
                rect(x = j*dim, y = i*dim, width = dim, height = dim, fill = fill)
            end
        end
    end
end

# ╔═╡ 6a08d680-eb95-11ea-17b5-d199cc14b12a
visualize1dim(res, 30)

# ╔═╡ 9db5c010-eb95-11ea-3257-652ddbd3939d
md"""## Classifying CAs

### Class 1

Evolution from any starting condition to the same uniform pattern, eg. rule0"""

# ╔═╡ adfe9aa0-eb95-11ea-1ce9-e5370aeec375
let
	x₀ = bitrand(21)
	x₀[1] = false; x₀[end] = false;
	res = step1dim(x₀, inttorule1dim(UInt8(0)), 1)
	visualize1dim(res, 30)
end

# ╔═╡ 1c6b5780-eb96-11ea-2bbd-755d35c27938
md"""### Class 2

Generation of a simple pattern with nested structure, i.e. a pattern that contains many smaller versions of itself, eg. rule50.

Example that looks like a Sierpinsi triangle (fractal): rule18."""

# ╔═╡ 2bcfc580-eb96-11ea-1837-c703fac6cd69
let
	x₀ = falses(129)
	x₀[65] = true
	res = step1dim(x₀, inttorule1dim(UInt8(18)), 63);
	visualize1dim(res, 6)
end

# ╔═╡ 4d3dd0e0-eb96-11ea-0f6a-8f6fb24e8f9d
md"""### Class 3

CAs that generate randomness, eg. rule30."""

# ╔═╡ 4d2f51f0-eb96-11ea-1c16-b7b7b81a3ca3
let
	x₀ = falses(201)
	x₀[101] = true
	res = step1dim(x₀, inttorule1dim(UInt8(30)), 99);
	visualize1dim(res, 4)
end

# ╔═╡ 8eb8d100-eb96-11ea-12f5-751c018540fc
md"""Center column as a sequence of bits, is hard to distinguish from a truly random sequence: pseudo-random number generators (PRNGs).

- regularities can be detected statistically
- a PRNG with finite amount of state will eventually repeat itself (period)
- underlying process is fundamentally deterministic (unlike some physical processes: thermodynamics or quantum mechanics)

This complex behavior is surprising (chaos is often associated with non-linear behavior of continuous time and space processes)."""

# ╔═╡ a075adee-eb96-11ea-3ddc-572f2795d6b5
md"""### Class 4

CAs that are Turing complete or universal, which means that they can compute any computable function, eg. rule110."""

# ╔═╡ aed6b5b0-eb96-11ea-2a3d-7f82f61ea518
let
	x₀ = bitrand(600)
	res = step1dim(x₀, inttorule1dim(UInt8(110)), 599);
	visualize1dim(res, 1)
end

# ╔═╡ c56df180-eb96-11ea-1a57-2d0571da6ac3
md"""- After about 100 steps, simple repeating patterns emerge, but there are a number of persistent structures that appear as disturbances. Some are vertical, other are diagonal and are called spaceships.

- Collisions between spaceships yields different results depending on their type and their phase. Some collisions annihilate both ships; other leaves one ship unchanged; still other yield one or more ships of different types.

- The collisions are the basis of computation in a rule110 CA. You can think of spaceships as signales that propagate through space, and collisions as gate that compute logical operations like AND and OR."""

# ╔═╡ f2a45860-eb96-11ea-2e5f-b1850621020b
md"""## Turing State-Machines

Based on [wikipedia: Turing Machine](https://en.wikipedia.org/wiki/Turing_machine).

A Turing machine is a mathematical model of computation that defines an abstract machine, which manipulates symbols on a tape according to a table of rules. Despite the model's simplicity, given any computer algorithm, a Turing machine capable of simulating that algorithm's logic can be constructed.

- A tape divided into cells, one next to the other. Each cell contains a symbol from some finite alphabet. The alphabet contains a special blank symbol (here written as '0') and one or more other symbols. The tape is assumed to be arbitrarily extendable to the left and to the right, i.e., the Turing machine is always supplied with as much tape as it needs for its computation. Cells that have not been written before are assumed to be filled with the blank symbol. In some models the tape has a left end marked with a special symbol; the tape extends or is indefinitely extensible to the right.
- A head that can read and write symbols on the tape and move the tape left and right one (and only one) cell at a time. In some models the head moves and the tape is stationary.
- A state register that stores the state of the Turing machine, one of finitely many. Among these is the special start state with which the state register is initialized. These states, writes Turing, replace the "state of mind" a person performing computations would ordinarily be in.
- A finite table of instructions that, given the state the machine is currently in and the symbol it is reading on the tape (symbol currently under the head), tells the machine to do the following in sequence:
    - Erase or write a symbol.
    - Move the head ( 'L' for one step left or 'R' for one step right or 'N' for staying in the same place).
    - Assume the same or a new state as prescribed.

Table of rules:

| Tape Symbol | State A   | State B   | State C   |
|:-----------:|-----------|-----------|-----------|
| 0           | 1 - R - B | 1 - L - A | 1 - L - B |
| 1           | 1 - L - C | 1 - R - B | 1 - R - H |"""

# ╔═╡ fa9e94e0-eb96-11ea-061c-a36a28f5e7b1
function applyrulebusybeaver(state, read)
    if state == 'A' && read == 0
        return 1, 'R', 'B'
    elseif state == 'A' && read == 1
        return 1, 'L', 'C'
    elseif state == 'B' && read == 0
        return 1, 'L', 'A'
    elseif state == 'B' && read == 1
        return 1, 'R', 'B'
    elseif state == 'C' && read == 0
        return 1, 'L', 'B'
    elseif state == 'C' && read == 1
        return 1, 'R', 'H'
    end
end

# ╔═╡ 11fa5930-eb97-11ea-274e-3f2c45958666
md"struct to represent the Turing State-Machine:"

# ╔═╡ 1a07f6a0-eb97-11ea-0f51-a3d742fcf260
mutable struct Turing
    tape :: Array{Int64}
    position :: Int64
    state :: Char
end

# ╔═╡ 7bfdef90-eb97-11ea-03ce-4fb081079815
function Base.show(io::IO, turing::Turing)
    print(io, turing.position, " - ", turing.state, ": ", turing.tape)
end

# ╔═╡ 9b61fa20-eb97-11ea-0f1c-c922c04a796f
md"implementation of a step:"

# ╔═╡ a68bf2c0-eb97-11ea-30ea-87d20ea07175
function stepturing!(turing, applyrule)
    if turing.state == 'H'
        error("Machine has stopped!")
    end
    read = turing.tape[turing.position]
    (write, dir, turing.state) = applyrule(turing.state, read)
    turing.tape[turing.position] = write
    if dir == 'L'
        if turing.position == length(turing.tape)
            push!(turing.tape, false)
        end
        turing.position += 1
    else
        if turing.position == 1
            pushfirst!(turing.tape, false)
        else
            turing.position -= 1
        end
    end
    nothing
end

# ╔═╡ a3263280-eb97-11ea-2e91-532bfa433f6c
md"output:"

# ╔═╡ be984390-ebb1-11ea-0729-39795172949b
turing = Turing(zeros(Int64, 11), 6, 'A')

# ╔═╡ ccda0380-ebb1-11ea-1677-11d9d171897c
begin
	stepturing!(turing, applyrulebusybeaver)
	turing
end

# ╔═╡ 74832baf-5779-47f6-907e-5288420bb1b9
with_terminal() do
	turing = Turing(zeros(Int64, 11), 6, 'A')
	println("- ", turing)
	try
		while true
			stepturing!(turing, applyrulebusybeaver)
			println("- ", turing)
		end
	catch
	end
end

# ╔═╡ Cell order:
# ╟─56f41ca0-eb93-11ea-1ea6-11b0e8bb9a7d
# ╠═310ce685-2661-4f32-bf14-91a4f4e569ce
# ╟─dab709d0-eb93-11ea-050e-3bfa6a5e1836
# ╟─f4b601b0-eb93-11ea-0969-7967a1e85c8b
# ╠═14a02c30-eb94-11ea-2114-5102409c8ae5
# ╟─1a915240-eb94-11ea-087f-231ade62690d
# ╠═2393fe0e-eb94-11ea-2858-b50503395d4a
# ╟─2dde4510-eb94-11ea-212a-5da1e7733bf6
# ╠═35cbd760-eb94-11ea-201a-6b58802154a1
# ╟─758a0c50-eb94-11ea-11a2-e3e007b089a9
# ╠═d8481030-eb94-11ea-1af4-db838adc37ed
# ╟─e393e8b2-eb94-11ea-3d70-d31f7ee89420
# ╠═eed0a600-eb94-11ea-1863-0f33980508ba
# ╟─0002a810-eb95-11ea-2ba8-bb2849bdec17
# ╠═04139c70-eb95-11ea-1759-298193ce97b0
# ╟─10100040-eb95-11ea-3a6d-271e63301b17
# ╠═1a508340-eb95-11ea-08f9-6f085748d7ff
# ╟─5d2a01f0-eb95-11ea-3dcd-b794e1b0d566
# ╠═60e8deb0-eb95-11ea-2bde-d9c259432318
# ╠═6a08d680-eb95-11ea-17b5-d199cc14b12a
# ╟─9db5c010-eb95-11ea-3257-652ddbd3939d
# ╠═bbf4ab90-eb95-11ea-174e-c3624a0d92cd
# ╠═adfe9aa0-eb95-11ea-1ce9-e5370aeec375
# ╟─1c6b5780-eb96-11ea-2bbd-755d35c27938
# ╠═2bcfc580-eb96-11ea-1837-c703fac6cd69
# ╟─4d3dd0e0-eb96-11ea-0f6a-8f6fb24e8f9d
# ╠═4d2f51f0-eb96-11ea-1c16-b7b7b81a3ca3
# ╟─8eb8d100-eb96-11ea-12f5-751c018540fc
# ╟─a075adee-eb96-11ea-3ddc-572f2795d6b5
# ╠═aed6b5b0-eb96-11ea-2a3d-7f82f61ea518
# ╟─c56df180-eb96-11ea-1a57-2d0571da6ac3
# ╟─f2a45860-eb96-11ea-2e5f-b1850621020b
# ╠═fa9e94e0-eb96-11ea-061c-a36a28f5e7b1
# ╟─11fa5930-eb97-11ea-274e-3f2c45958666
# ╠═1a07f6a0-eb97-11ea-0f51-a3d742fcf260
# ╠═7bfdef90-eb97-11ea-03ce-4fb081079815
# ╟─9b61fa20-eb97-11ea-0f1c-c922c04a796f
# ╠═a68bf2c0-eb97-11ea-30ea-87d20ea07175
# ╟─a3263280-eb97-11ea-2e91-532bfa433f6c
# ╠═be984390-ebb1-11ea-0729-39795172949b
# ╠═ccda0380-ebb1-11ea-1677-11d9d171897c
# ╠═43627847-ce0d-4969-b6e1-2d59a54d1fd2
# ╠═74832baf-5779-47f6-907e-5288420bb1b9


text:
# This file was generated, do not modify it. # hide
foo(v) = v[1:2]
combine(df, :a => maximum, :b => foo)


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))

using Pluto, MLCourse, MLJ, MLJLinearModels, StatsPlots, DataFrames,
      Distributions

X, y = make_regression()

mach = fit!(machine(LinearRegressor(), X, y), verbosity = 0)

predict(mach)

plot(rand(10), rand(10))

df = DataFrame(a = rand(10), b = rand(10))
@df df scatter(:a, :b)

rand(Bernoulli(.4))

redirect_stdout(Pipe()) do

session = Pluto.ServerSession()
session.options.server.port = 40404
session.options.server.launch_browser = false
session.options.security.require_secret_for_access = false

path = tempname()
original = joinpath(@__DIR__, "..", "index.jl")
# so that we don't overwrite the file:
Pluto.readwrite(original, path)

# @info "Loading notebook"
nb = Pluto.load_notebook(Pluto.tamepath(path));
session.notebooks[nb.notebook_id] = nb;

# @info "Running notebook"
Pluto.update_save_run!(session, nb, nb.cells; run_async=false, prerender_text=true)

# nice! we ran the notebook, so we already precompiled a lot

# @info "Starting HTTP server"
# next, we'll run the HTTP server which needs a bit of nasty code
t = @async Pluto.run(session)

sleep(5)
# download("http://localhost:40404/")

# this is async because it blocks for some reason
@async Base.throwto(t, InterruptException())
sleep(2) # i am pulling these numbers out of thin air

end
@info "Warmup done"

max_stars_repo_name: baltisberger/MLCourse

max_issues_repo_name: baltisberger/MLCourse

max_forks_repo_name: baltisberger/MLCourse