text:
const  LDA_X                      =      1 # Slater exchange
const  LDA_C_WIGNER               =      2 # Wigner
const  LDA_C_RPA                  =      3 # Random Phase Approximation (RPA)
const  LDA_C_HL                   =      4 # Hedin & Lundqvist
const  LDA_C_GL                   =      5 # Gunnarson & Lundqvist
const  LDA_C_XALPHA               =      6 # Slater's Xalpha
const  LDA_C_VWN                  =      7 # Vosko, Wilk & Nusair (VWN5)
const  LDA_C_VWN_RPA              =      8 # Vosko, Wilk & Nusair (VWN5_RPA)
const  LDA_C_PZ                   =      9 # Perdew & Zunger
const  LDA_C_PZ_MOD               =     10 # Perdew & Zunger (Modified)
const  LDA_C_OB_PZ                =     11 # Ortiz & Ballone (PZ parametrization)
const  LDA_C_PW                   =     12 # Perdew & Wang
const  LDA_C_PW_MOD               =     13 # Perdew & Wang (modified)
const  LDA_C_OB_PW                =     14 # Ortiz & Ballone (PW parametrization)
const  LDA_C_2D_AMGB              =     15 # AMGB (for 2D systems)
const  LDA_C_2D_PRM               =     16 # PRM (for 2D systems)
const  LDA_C_VBH                  =     17 # von Barth & Hedin
const  LDA_C_1D_CSC               =     18 # Casula, Sorella & Senatore
const  LDA_X_2D                   =     19 # Slater exchange
const  LDA_XC_TETER93             =     20 # Teter 93
const  LDA_X_1D                   =     21 # Exchange in 1D
const  LDA_C_ML1                  =     22 # Modified LSD (version 1) of Proynov and Salahub
const  LDA_C_ML2                  =     23 # Modified LSD (version 2) of Proynov and Salahub
const  LDA_C_GOMBAS               =     24 # Gombas
const  LDA_C_PW_RPA               =     25 # Perdew & Wang (fit to the RPA energy)
const  LDA_C_1D_LOOS              =     26 # P-F Loos correlation LDA
const  LDA_C_RC04                 =     27 # Ragot-Cortona
const  LDA_C_VWN_1                =     28 # Vosko, Wilk & Nusair (VWN1)
const  LDA_C_VWN_2                =     29 # Vosko, Wilk & Nusair (VWN2)
const  LDA_C_VWN_3                =     30 # Vosko, Wilk & Nusair (VWN3)
const  LDA_C_VWN_4                =     31 # Vosko, Wilk & Nusair (VWN4)
const  GGA_X_GAM                  =     32 # Minnesota GAM exhange functional
const  GGA_C_GAM                  =     33 # GAM functional from Minnesota
const  GGA_X_HCTH_A               =     34 # HCTH-A
const  GGA_X_EV93                 =     35 # Engel and Vosko
const  HYB_MGGA_X_DLDF            =     36 # Dispersionless Density Functional
const  MGGA_C_DLDF                =     37 # Dispersionless Density Functional
const  GGA_X_BCGP                 =     38 # Burke, Cancio, Gould, and Pittalis
const  GGA_C_BCGP                 =     39 # Burke, Cancio, Gould, and Pittalis
const  GGA_X_LAMBDA_OC2_N         =     40 # lambda_OC2(N) version of PBE
const  GGA_X_B86_R                =     41 # Revised Becke 86 with modified gradient correction
const  MGGA_XC_ZLP                =     42 # Zhao, Levy & Parr, Eq. (21)
const  LDA_XC_ZLP                 =     43 # Zhao, Levy & Parr, Eq. (20)
const  GGA_X_LAMBDA_CH_N          =     44 # lambda_CH(N) version of PBE
const  GGA_X_LAMBDA_LO_N          =     45 # lambda_LO(N) version of PBE
const  GGA_X_HJS_B88_V2           =     46 # HJS screened exchange B88 corrected version
const  GGA_C_Q2D                  =     47 # Chiodo et al
const  GGA_X_Q2D                  =     48 # Chiodo et al
const  GGA_X_PBE_MOL              =     49 # Reparametrized PBE by del Campo, Gazquez, Trickey & Vela
const  LDA_K_TF                   =     50 # Thomas-Fermi kinetic energy
const  LDA_K_LP                   =     51 # Lee and Parr Gaussian ansatz for the kinetic energy
const  GGA_K_TFVW                 =     52 # Thomas-Fermi plus von Weiszaecker correction
const  GGA_K_REVAPBEINT           =     53 # interpolated version of revAPBE
const  GGA_K_APBEINT              =     54 # interpolated version of APBE
const  GGA_K_REVAPBE              =     55 # revised APBE
const  GGA_X_AK13                 =     56 # Armiento & Kuemmel 2013
const  GGA_K_MEYER                =     57 # Meyer,  Wang, and Young
const  GGA_X_LV_RPW86             =     58 # Berland and Hyldgaard
const  GGA_X_PBE_TCA              =     59 # PBE revised by Tognetti et al
const  GGA_X_PBEINT               =     60 # PBE for hybrid interfaces
const  GGA_C_ZPBEINT              =     61 # spin-dependent gradient correction to PBEint
const  GGA_C_PBEINT               =     62 # PBE for hybrid interfaces
const  GGA_C_ZPBESOL              =     63 # spin-dependent gradient correction to PBEsol
const  MGGA_XC_OTPSS_D            =     64 # oTPSS-D functional of Goerigk and Grimme
const  GGA_XC_OPBE_D              =     65 # oPBE-D functional of Goerigk and Grimme
const  GGA_XC_OPWLYP_D            =     66 # oPWLYP-D functional of Goerigk and Grimme
const  GGA_XC_OBLYP_D             =     67 # oBLYP-D functional of Goerigk and Grimme
const  GGA_X_VMT84_GE             =     68 # VMT{8,4} with constraint satisfaction with mu = mu_GE
const  GGA_X_VMT84_PBE            =     69 # VMT{8,4} with constraint satisfaction with mu = mu_PBE
const  GGA_X_VMT_GE               =     70 # Vela, Medel, and Trickey with mu = mu_GE
const  GGA_X_VMT_PBE              =     71 # Vela, Medel, and Trickey with mu = mu_PBE
const  MGGA_C_CS                  =     72 # Colle and Salvetti
const  MGGA_C_MN12_SX             =     73 # Minnesota MN12-SX correlation functional
const  MGGA_C_MN12_L              =     74 # Minnesota MN12-L correlation functional
const  MGGA_C_M11_L               =     75 # Minnesota M11-L correlation functional
const  MGGA_C_M11                 =     76 # Minnesota M11 correlation functional
const  MGGA_C_M08_SO              =     77 # Minnesota M08-SO correlation functional
const  MGGA_C_M08_HX              =     78 # Minnesota M08 correlation functional
const  GGA_C_N12_SX               =     79 # Minnesota N12-SX functional
const  GGA_C_N12                  =     80 # Minnesota N12 functional
const  HYB_GGA_X_N12_SX           =     81 # Minnesota N12-SX exchange functional
const  GGA_X_N12                  =     82 # Minnesota N12 exchange functional
const  GGA_C_REGTPSS              =     83 # regularized TPSS correlation
const  GGA_C_OP_XALPHA            =     84 # one-parameter progressive functional (Xalpha version)
const  GGA_C_OP_G96               =     85 # one-parameter progressive functional (G96 version)
const  GGA_C_OP_PBE               =     86 # one-parameter progressive functional (PBE version)
const  GGA_C_OP_B88               =     87 # one-parameter progressive functional (B88 version)
const  GGA_C_FT97                 =     88 # Filatov & Thiel correlation
const  GGA_C_SPBE                 =     89 # PBE correlation to be used with the SSB exchange
const  GGA_X_SSB_SW               =     90 # Swart, Sola and Bickelhaupt correction to PBE
const  GGA_X_SSB                  =     91 # Swart, Sola and Bickelhaupt
const  GGA_X_SSB_D                =     92 # Swart, Sola and Bickelhaupt dispersion
const  GGA_XC_HCTH_407P           =     93 # HCTH/407+
const  GGA_XC_HCTH_P76            =     94 # HCTH p=7/6
const  GGA_XC_HCTH_P14            =     95 # HCTH p=1/4
const  GGA_XC_B97_GGA1            =     96 # Becke 97 GGA-1
const  GGA_C_HCTH_A               =     97 # HCTH-A
const  GGA_X_BPCCAC               =     98 # BPCCAC (GRAC for the energy)
const  GGA_C_REVTCA               =     99 # Tognetti, Cortona, Adamo (revised)
const  GGA_C_TCA                  =    100 # Tognetti, Cortona, Adamo
const  GGA_X_PBE                  =    101 # Perdew, Burke & Ernzerhof
const  GGA_X_PBE_R                =    102 # Revised PBE from Zhang & Yang
const  GGA_X_B86                  =    103 # Becke 86
const  GGA_X_HERMAN               =    104 # Herman Xalphabeta GGA
const  GGA_X_B86_MGC              =    105 # Becke 86 with modified gradient correction
const  GGA_X_B88                  =    106 # Becke 88
const  GGA_X_G96                  =    107 # Gill 96
const  GGA_X_PW86                 =    108 # Perdew & Wang 86
const  GGA_X_PW91                 =    109 # Perdew & Wang 91
const  GGA_X_OPTX                 =    110 # Handy & Cohen OPTX 01
const  GGA_X_DK87_R1              =    111 # dePristo & Kress 87 version R1
const  GGA_X_DK87_R2              =    112 # dePristo & Kress 87 version R2
const  GGA_X_LG93                 =    113 # Lacks & Gordon 93
const  GGA_X_FT97_A               =    114 # Filatov & Thiel 97 (version A)
const  GGA_X_FT97_B               =    115 # Filatov & Thiel 97 (version B)
const  GGA_X_PBE_SOL              =    116 # Perdew, Burke & Ernzerhof SOL
const  GGA_X_RPBE                 =    117 # Hammer, Hansen, and Norskov
const  GGA_X_WC                   =    118 # Wu & Cohen
const  GGA_X_MPW91                =    119 # mPW91 of Adamo & Barone
const  GGA_X_AM05                 =    120 # Armiento & Mattsson 05
const  GGA_X_PBEA                 =    121 # Madsen 07
const  GGA_X_MPBE                 =    122 # Adamo & Barone modification to PBE
const  GGA_X_XPBE                 =    123 # Extended PBE by Xu & Goddard III
const  GGA_X_2D_B86_MGC           =    124 # Becke 86 with modified gradient correction for 2D
const  GGA_X_BAYESIAN             =    125 # Bayesian best fit for the enhancement factor
const  GGA_X_PBE_JSJR             =    126 # Reparametrized PBE by Pedroza, Silva & Capelle
const  GGA_X_2D_B88               =    127 # Becke 88 in 2D
const  GGA_X_2D_B86               =    128 # Becke 86 in 2D
const  GGA_X_2D_PBE               =    129 # Perdew, Burke & Ernzerhof in 2D
const  GGA_C_PBE                  =    130 # Perdew, Burke & Ernzerhof
const  GGA_C_LYP                  =    131 # Lee, Yang & Parr
const  GGA_C_P86                  =    132 # Perdew 86
const  GGA_C_PBE_SOL              =    133 # Perdew, Burke & Ernzerhof SOL
const  GGA_C_PW91                 =    134 # Perdew & Wang 91
const  GGA_C_AM05                 =    135 # Armiento & Mattsson 05
const  GGA_C_XPBE                 =    136 # Extended PBE by Xu & Goddard III
const  GGA_C_LM                   =    137 # Langreth & Mehl
const  GGA_C_PBE_JRGX             =    138 # Reparametrized PBE by Pedroza, Silva & Capelle
const  GGA_X_OPTB88_VDW           =    139 # opt-Becke 88 for vdW
const  GGA_X_PBEK1_VDW            =    140 # Reparametrized PBE for vdW
const  GGA_X_OPTPBE_VDW           =    141 # Reparametrized PBE for vdW
const  GGA_X_RGE2                 =    142 # Regularized PBE
const  GGA_C_RGE2                 =    143 # Regularized PBE
const  GGA_X_RPW86                =    144 # Refitted Perdew & Wang 86
const  GGA_X_KT1                  =    145 # Exchange part of Keal and Tozer version 1
const  GGA_XC_KT2                 =    146 # Keal and Tozer, version 2
const  GGA_C_WL                   =    147 # Wilson & Levy
const  GGA_C_WI                   =    148 # Wilson & Ivanov
const  GGA_X_MB88                 =    149 # Modified Becke 88 for proton transfer
const  GGA_X_SOGGA                =    150 # Second-order generalized gradient approximation
const  GGA_X_SOGGA11              =    151 # Second-order generalized gradient approximation 2011
const  GGA_C_SOGGA11              =    152 # Second-order generalized gradient approximation 2011
const  GGA_C_WI0                  =    153 # Wilson & Ivanov initial version
const  GGA_XC_TH1                 =    154 # Tozer and Handy v. 1
const  GGA_XC_TH2                 =    155 # Tozer and Handy v. 2
const  GGA_XC_TH3                 =    156 # Tozer and Handy v. 3
const  GGA_XC_TH4                 =    157 # Tozer and Handy v. 4
const  GGA_X_C09X                 =    158 # C09x to be used with the VdW of Rutgers-Chalmers
const  GGA_C_SOGGA11_X            =    159 # To be used with HYB_GGA_X_SOGGA11_X
const  GGA_X_LB                   =    160 # van Leeuwen & Baerends
const  GGA_XC_HCTH_93             =    161 # HCTH/93
const  GGA_XC_HCTH_120            =    162 # HCTH/120
const  GGA_XC_HCTH_147            =    163 # HCTH/147
const  GGA_XC_HCTH_407            =    164 # HCTH/407
const  GGA_XC_EDF1                =    165 # EDF1
const  GGA_XC_XLYP                =    166 # XLYP
const  GGA_XC_KT1                 =    167 # Keal and Tozer, version 1
const  GGA_XC_B97_D               =    170 # Becke 97-D
const  GGA_XC_PBE1W               =    173 # PBE1W
const  GGA_XC_MPWLYP1W            =    174 # mPWLYP1w
const  GGA_XC_PBELYP1W            =    175 # PBELYP1W
const  GGA_X_LBM                  =    182 # van Leeuwen & Baerends modified
const  GGA_X_OL2                  =    183 # Exchange form based on Ou-Yang and Levy v.2
const  GGA_X_APBE                 =    184 # mu fixed from the semiclassical neutral atom
const  GGA_K_APBE                 =    185 # mu fixed from the semiclassical neutral atom
const  GGA_C_APBE                 =    186 # mu fixed from the semiclassical neutral atom
const  GGA_K_TW1                  =    187 # Tran and Wesolowski set 1 (Table II)
const  GGA_K_TW2                  =    188 # Tran and Wesolowski set 2 (Table II)
const  GGA_K_TW3                  =    189 # Tran and Wesolowski set 3 (Table II)
const  GGA_K_TW4                  =    190 # Tran and Wesolowski set 4 (Table II)
const  GGA_X_HTBS                 =    191 # Haas, Tran, Blaha, and Schwarz
const  GGA_X_AIRY                 =    192 # Constantin et al based on the Airy gas
const  GGA_X_LAG                  =    193 # Local Airy Gas
const  GGA_XC_MOHLYP              =    194 # Functional for organometallic chemistry
const  GGA_XC_MOHLYP2             =    195 # Functional for barrier heights
const  GGA_XC_TH_FL               =    196 # Tozer and Handy v. FL
const  GGA_XC_TH_FC               =    197 # Tozer and Handy v. FC
const  GGA_XC_TH_FCFO             =    198 # Tozer and Handy v. FCFO
const  GGA_XC_TH_FCO              =    199 # Tozer and Handy v. FCO
const  GGA_C_OPTC                 =    200 # Optimized correlation functional of Cohen and Handy
const  MGGA_X_LTA                 =    201 # Local tau approximation
const  MGGA_X_TPSS                =    202 # Tao, Perdew, Staroverov & Scuseria
const  MGGA_X_M06_L               =    203 # Minnesota M06-L exchange functional
const  MGGA_X_GVT4                =    204 # GVT4 (X part of VSXC)
const  MGGA_X_TAU_HCTH            =    205 # tau-HCTH from Boese and Handy
const  MGGA_X_BR89                =    206 # Becke-Roussel 89
const  MGGA_X_BJ06                =    207 # Becke & Johnson 06
const  MGGA_X_TB09                =    208 # Tran & Blaha 09
const  MGGA_X_RPP09               =    209 # Rasanen, Pittalis & Proetto 09
const  MGGA_X_2D_PRHG07           =    210 # Pittalis-Rasanen-Helbig-Gross 2007
const  MGGA_X_2D_PRHG07_PRP10     =    211 # PRHG07 with Pittalis-Rasanen-Proetto 2010 correction
const  MGGA_X_REVTPSS             =    212 # revised Tao, Perdew, Staroverov & Scuseria
const  MGGA_X_PKZB                =    213 # Perdew, Kurth, Zupan, and Blaha
const  MGGA_X_MS0                 =    221 # MS exchange of Sun, Xiao, and Ruzsinszky
const  MGGA_X_MS1                 =    222 # MS1 exchange of Sun, et al
const  MGGA_X_MS2                 =    223 # MS2 exchange of Sun, et al
const  HYB_MGGA_X_MS2H            =    224 # MS2 hybrid exchange of Sun, et al
const  MGGA_X_M11_L               =    226 # Minnesota M11-L exchange functional
const  MGGA_X_MN12_L              =    227 # Minnesota MN12-L exchange functional
const  MGGA_XC_CC06               =    229 # Cancio and Chou 2006
const  MGGA_X_MK00                =    230 # Exchange for accurate virtual orbital energies
const  MGGA_C_TPSS                =    231 # Tao, Perdew, Staroverov & Scuseria
const  MGGA_C_VSXC                =    232 # VSXC (correlation part)
const  MGGA_C_M06_L               =    233 # Minnesota M06-L correlation functional
const  MGGA_C_M06_HF              =    234 # Minnesota M06-HF correlation functional
const  MGGA_C_M06                 =    235 # Minnesota M06 correlation functional
const  MGGA_C_M06_2X              =    236 # Minnesota M06-2X correlation functional
const  MGGA_C_M05                 =    237 # Minnesota M05 correlation functional
const  MGGA_C_M05_2X              =    238 # Minnesota M05-2X correlation functional
const  MGGA_C_PKZB                =    239 # Perdew, Kurth, Zupan, and Blaha
const  MGGA_C_BC95                =    240 # Becke correlation 95
const  MGGA_C_REVTPSS             =    241 # revised TPSS correlation
const  MGGA_XC_TPSSLYP1W          =    242 # TPSSLYP1W
const  MGGA_X_MK00B               =    243 # Exchange for accurate virtual orbital energies (v. B)
const  MGGA_X_BLOC                =    244 # functional with balanced localization
const  MGGA_X_MODTPSS             =    245 # Modified Tao, Perdew, Staroverov & Scuseria
const  GGA_C_PBELOC               =    246 # Semilocal dynamical correlation
const  MGGA_C_TPSSLOC             =    247 # Semilocal dynamical correlation
const  HYB_MGGA_X_MN12_SX         =    248 # Minnesota MN12-SX hybrid exchange functional
const  MGGA_X_MBEEF               =    249 # mBEEF exchange
const  MGGA_X_MBEEFVDW            =    250 # mBEEF-vdW exchange
const  MGGA_XC_B97M_V             =    254 # B97M-V exchange-correlation functional
const  GGA_XC_VV10                =    255 # Vydrov and Van Voorhis
const  MGGA_X_MVS                 =    257 # MVS exchange of Sun, Perdew, and Ruzsinszky
const  GGA_C_PBEFE                =    258 # PBE for formation energies
const  LDA_XC_KSDT                =    259 # Karasiev, Sjostrom, Dufty & Trickey
const  MGGA_X_MN15_L              =    260 # Minnesota MN15-L exchange functional
const  MGGA_C_MN15_L              =    261 # Minnesota MN15-L correlation functional
const  GGA_C_OP_PW91              =    262 # one-parameter progressive functional (PW91 version)
const  MGGA_X_SCAN                =    263 # SCAN exchange of Sun, Ruzsinszky, and Perdew
const  HYB_MGGA_X_SCAN0           =    264 # SCAN hybrid exchange (SCAN0)
const  GGA_X_PBEFE                =    265 # PBE for formation energies
const  HYB_GGA_XC_B97_1P          =    266 # version of B97 by Cohen and Handy
const  MGGA_C_SCAN                =    267 # SCAN correlation of Sun, Ruzsinszky, and Perdew
const  HYB_MGGA_X_MN15            =    268 # Minnesota MN15 hybrid exchange functional
const  MGGA_C_MN15                =    269 # Minnesota MN15 correlation functional
const  GGA_X_CAP                  =    270 # Correct Asymptotic Potential
const  GGA_X_EB88                 =    271 # Non-empirical (excogitated) B88 functional of Becke and Elliott
const  GGA_C_PBE_MOL              =    272 # Reparametrized PBE by del Campo, Gazquez, Trickey & Vela
const  HYB_GGA_XC_PBE_MOL0        =    273 # PBEmol0
const  HYB_GGA_XC_PBE_SOL0        =    274 # PBEsol0
const  HYB_GGA_XC_PBEB0           =    275 # PBEbeta0
const  HYB_GGA_XC_PBE_MOLB0       =    276 # PBEmolbeta0
const  GGA_K_ABSP3                =    277 # gamma-TFvW form by Acharya et al [g = 1 - 1.513/N^0.35]
const  GGA_K_ABSP4                =    278 # gamma-TFvW form by Acharya et al [g = l = 1/(1 + 1.332/N^(1/3))]
const  HYB_MGGA_X_BMK             =    279 # Boese-Martin for kinetics
const  GGA_C_BMK                  =    280 # Boese-Martin for kinetics
const  GGA_C_TAU_HCTH             =    281 # correlation part of tau-hcth
const  HYB_MGGA_X_TAU_HCTH        =    282 # Hybrid version of tau-HCTH
const  GGA_C_HYB_TAU_HCTH         =    283 # correlation part of hyb-tau-hcth
const  MGGA_X_B00                 =    284 # Becke 2000
const  GGA_X_BEEFVDW              =    285 # BEEF-vdW exchange
const  GGA_XC_BEEFVDW             =    286 # BEEF-vdW exchange-correlation
const  LDA_C_CHACHIYO             =    287 # Chachiyo simple 2 parameter correlation
const  MGGA_XC_HLE17              =    288 # high local exchange 2017
const  LDA_C_LP96                 =    289 # Liu-Parr correlation
const  HYB_GGA_XC_PBE50           =    290 # PBE50
const  GGA_X_PBETRANS             =    291 # Gradient-regulated connection-based correction for the PBE exchange
const  MGGA_C_SCAN_RVV10          =    292 # SCAN + rVV10 correlation
const  MGGA_X_REVM06_L            =    293 # Minnesota revM06-L exchange functional
const  MGGA_C_REVM06_L            =    294 # Minnesota revM06-L correlation functional
const  HYB_MGGA_X_M08_HX          =    295 # Minnesota M08-HX hybrid exchange functional
const  HYB_MGGA_X_M08_SO          =    296 # Minnesota M08-SO hybrid exchange functional
const  HYB_MGGA_X_M11             =    297 # Minnesota M11 hybrid exchange functional
const  GGA_X_CHACHIYO             =    298 # Chachiyo exchange
const  HYB_GGA_XC_B3PW91          =    401 # The original (ACM, B3PW91) hybrid of Becke
const  HYB_GGA_XC_B3LYP           =    402 # B3LYP
const  HYB_GGA_XC_B3P86           =    403 # B3P86
const  HYB_GGA_XC_O3LYP           =    404 # O3LYP
const  HYB_GGA_XC_MPW1K           =    405 # mPW1K
const  HYB_GGA_XC_PBEH            =    406 # PBEH (PBE0)
const  HYB_GGA_XC_B97             =    407 # Becke 97
const  HYB_GGA_XC_B97_1           =    408 # Becke 97-1
const  HYB_GGA_XC_B97_2           =    410 # Becke 97-2
const  HYB_GGA_XC_X3LYP           =    411 # X3LYP
const  HYB_GGA_XC_B1WC            =    412 # B1WC
const  HYB_GGA_XC_B97_K           =    413 # Boese-Martin for Kinetics
const  HYB_GGA_XC_B97_3           =    414 # Becke 97-3
const  HYB_GGA_XC_MPW3PW          =    415 # MPW3PW of Adamo & Barone
const  HYB_GGA_XC_B1LYP           =    416 # B1LYP
const  HYB_GGA_XC_B1PW91          =    417 # B1PW91
const  HYB_GGA_XC_MPW1PW          =    418 # mPW1PW
const  HYB_GGA_XC_MPW3LYP         =    419 # MPW3LYP
const  HYB_GGA_XC_SB98_1A         =    420 # SB98 (1a)
const  HYB_GGA_XC_SB98_1B         =    421 # SB98 (1b)
const  HYB_GGA_XC_SB98_1C         =    422 # SB98 (1c)
const  HYB_GGA_XC_SB98_2A         =    423 # SB98 (2a)
const  HYB_GGA_XC_SB98_2B         =    424 # SB98 (2b)
const  HYB_GGA_XC_SB98_2C         =    425 # SB98 (2c)
const  HYB_GGA_X_SOGGA11_X        =    426 # Hybrid based on SOGGA11 form
const  HYB_GGA_XC_HSE03           =    427 # HSE03
const  HYB_GGA_XC_HSE06           =    428 # HSE06
const  HYB_GGA_XC_HJS_PBE         =    429 # HJS hybrid screened exchange PBE version
const  HYB_GGA_XC_HJS_PBE_SOL     =    430 # HJS hybrid screened exchange PBE_SOL version
const  HYB_GGA_XC_HJS_B88         =    431 # HJS hybrid screened exchange B88 version
const  HYB_GGA_XC_HJS_B97X        =    432 # HJS hybrid screened exchange B97x version
const  HYB_GGA_XC_CAM_B3LYP       =    433 # CAM version of B3LYP
const  HYB_GGA_XC_TUNED_CAM_B3LYP =    434 # CAM version of B3LYP, tuned for excitations and properties
const  HYB_GGA_XC_BHANDH          =    435 # BHandH
const  HYB_GGA_XC_BHANDHLYP       =    436 # BHandHLYP
const  HYB_GGA_XC_MB3LYP_RC04     =    437 # B3LYP with RC04 LDA
const  HYB_MGGA_X_M05             =    438 # Minnesota M05 hybrid exchange functional
const  HYB_MGGA_X_M05_2X          =    439 # Minnesota M05-2X hybrid exchange functional
const  HYB_MGGA_XC_B88B95         =    440 # Mixture of B88 with BC95 (B1B95)
const  HYB_MGGA_XC_B86B95         =    441 # Mixture of B86 with BC95
const  HYB_MGGA_XC_PW86B95        =    442 # Mixture of PW86 with BC95
const  HYB_MGGA_XC_BB1K           =    443 # Mixture of B88 with BC95 from Zhao and Truhlar
const  HYB_MGGA_X_M06_HF          =    444 # Minnesota M06-HF hybrid exchange functional
const  HYB_MGGA_XC_MPW1B95        =    445 # Mixture of mPW91 with BC95 from Zhao and Truhlar
const  HYB_MGGA_XC_MPWB1K         =    446 # Mixture of mPW91 with BC95 for kinetics
const  HYB_MGGA_XC_X1B95          =    447 # Mixture of X with BC95
const  HYB_MGGA_XC_XB1K           =    448 # Mixture of X with BC95 for kinetics
const  HYB_MGGA_X_M06             =    449 # Minnesota M06 hybrid exchange functional
const  HYB_MGGA_X_M06_2X          =    450 # Minnesota M06-2X hybrid exchange functional
const  HYB_MGGA_XC_PW6B95         =    451 # Mixture of PW91 with BC95 from Zhao and Truhlar
const  HYB_MGGA_XC_PWB6K          =    452 # Mixture of PW91 with BC95 from Zhao and Truhlar for kinetics
const  HYB_GGA_XC_MPWLYP1M        =    453 # MPW with 1 par. for metals/LYP
const  HYB_GGA_XC_REVB3LYP        =    454 # Revised B3LYP
const  HYB_GGA_XC_CAMY_BLYP       =    455 # CAMY version of BLYP
const  HYB_GGA_XC_PBE0_13         =    456 # PBE0-1/3
const  HYB_MGGA_XC_TPSSH          =    457 # TPSSh
const  HYB_MGGA_XC_REVTPSSH       =    458 # revTPSSh
const  HYB_GGA_XC_B3LYPS          =    459 # B3LYP*
const  HYB_GGA_XC_WB97            =    463 # wB97 range-separated functional
const  HYB_GGA_XC_WB97X           =    464 # wB97X range-separated functional
const  HYB_GGA_XC_LRC_WPBEH       =    465 # Long-range corrected short-range hybrid PBE (LRC-wPBEh) by Rohrdanz, Martins and Herbert
const  HYB_GGA_XC_WB97X_V         =    466 # wB97X-V range-separated functional
const  HYB_GGA_XC_LCY_PBE         =    467 # LCY version of PBE
const  HYB_GGA_XC_LCY_BLYP        =    468 # LCY version of BLYP
const  HYB_GGA_XC_LC_VV10         =    469 # Vydrov and Van Voorhis
const  HYB_GGA_XC_CAMY_B3LYP      =    470 # CAMY version of B3LYP
const  HYB_GGA_XC_WB97X_D         =    471 # wB97D range-separated functional
const  HYB_GGA_XC_HPBEINT         =    472 # hPBEint
const  HYB_GGA_XC_LRC_WPBE        =    473 # Long-range corrected PBE (LRC-wPBE) by Rohrdanz, Martins and Herbert
const  HYB_MGGA_X_MVSH            =    474 # MVSh hybrid exchange functional
const  HYB_GGA_XC_B3LYP5          =    475 # B3LYP with VWN functional 5 instead of RPA
const  HYB_GGA_XC_EDF2            =    476 # EDF2
const  HYB_GGA_XC_CAP0            =    477 # Correct Asymptotic Potential hybrid
const  HYB_GGA_XC_LC_WPBE         =    478 # Long-range corrected PBE (LC-wPBE) by Vydrov and Scuseria
const  HYB_GGA_XC_HSE12           =    479 # HSE12
const  HYB_GGA_XC_HSE12S          =    480 # HSE12 (short-range version)
const  HYB_GGA_XC_HSE_SOL         =    481 # HSEsol
const  HYB_GGA_XC_CAM_QTP_01      =    482 # CAM-B3LYP retuned using ionization potentials of water
const  HYB_GGA_XC_MPW1LYP         =    483 # mPW1LYP
const  HYB_GGA_XC_MPW1PBE         =    484 # mPW1PBE
const  HYB_GGA_XC_KMLYP           =    485 # Kang-Musgrave hybrid
const  GGA_K_VW                   =    500 # von Weiszaecker correction to Thomas-Fermi
const  GGA_K_GE2                  =    501 # Second-order gradient expansion of the kinetic energy density
const  GGA_K_GOLDEN               =    502 # TF-lambda-vW form by Golden (l = 13/45)
const  GGA_K_YT65                 =    503 # TF-lambda-vW form by Yonei and Tomishima (l = 1/5)
const  GGA_K_BALTIN               =    504 # TF-lambda-vW form by Baltin (l = 5/9)
const  GGA_K_LIEB                 =    505 # TF-lambda-vW form by Lieb (l = 0.185909191)
const  GGA_K_ABSP1                =    506 # gamma-TFvW form by Acharya et al [g = 1 - 1.412/N^(1/3)]
const  GGA_K_ABSP2                =    507 # gamma-TFvW form by Acharya et al [g = 1 - 1.332/N^(1/3)]
const  GGA_K_GR                   =    508 # gamma-TFvW form by Gazquez and Robles
const  GGA_K_LUDENA               =    509 # gamma-TFvW form by Ludena
const  GGA_K_GP85                 =    510 # gamma-TFvW form by Ghosh and Parr
const  GGA_K_PEARSON              =    511 # Pearson 1992
const  GGA_K_OL1                  =    512 # Ou-Yang and Levy v.1
const  GGA_K_OL2                  =    513 # Ou-Yang and Levy v.2
const  GGA_K_FR_B88               =    514 # Fuentealba & Reyes (B88 version)
const  GGA_K_FR_PW86              =    515 # Fuentealba & Reyes (PW86 version)
const  GGA_K_DK                   =    516 # DePristo and Kress
const  GGA_K_PERDEW               =    517 # Perdew
const  GGA_K_VSK                  =    518 # Vitos, Skriver, and Kollar
const  GGA_K_VJKS                 =    519 # Vitos, Johansson, Kollar, and Skriver
const  GGA_K_ERNZERHOF            =    520 # Ernzerhof
const  GGA_K_LC94                 =    521 # Lembarki & Chermette
const  GGA_K_LLP                  =    522 # Lee, Lee & Parr
const  GGA_K_THAKKAR              =    523 # Thakkar 1992
const  GGA_X_WPBEH                =    524 # short-range part of the PBE (default w=0 gives PBEh)
const  GGA_X_HJS_PBE              =    525 # HJS screened exchange PBE version
const  GGA_X_HJS_PBE_SOL          =    526 # HJS screened exchange PBE_SOL version
const  GGA_X_HJS_B88              =    527 # HJS screened exchange B88 version
const  GGA_X_HJS_B97X             =    528 # HJS screened exchange B97x version
const  GGA_X_ITYH                 =    529 # Short-range recipe for exchange GGA functionals
const  GGA_X_SFAT                 =    530 # Short-range recipe for exchange GGA functionals - Yukawa
const  HYB_MGGA_XC_WB97M_V        =    531 # wB97M-V exchange-correlation functional
const  LDA_X_REL                  =    532 # Slater exchange with relativistic corrections
const  GGA_X_SG4                  =    533 # Semiclassical GGA at fourth order
const  GGA_C_SG4                  =    534 # Semiclassical GGA at fourth order
const  GGA_X_GG99                 =    535 # Gilbert and Gill 1999
const  LDA_XC_1D_EHWLRG_1         =    536 # LDA constructed from slab-like systems of 1 electron
const  LDA_XC_1D_EHWLRG_2         =    537 # LDA constructed from slab-like systems of 2 electrons
const  LDA_XC_1D_EHWLRG_3         =    538 # LDA constructed from slab-like systems of 3 electrons
const  GGA_X_PBEPOW               =    539 # PBE power
const  MGGA_X_TM                  =    540 # Tao and Mo 2016
const  MGGA_X_VT84                =    541 # meta-GGA version of VT{8,4} GGA
const  MGGA_X_SA_TPSS             =    542 # TPSS with correct surface asymptotics
const  MGGA_K_PC07                =    543 # Perdew and Constantin 2007
const  GGA_X_KGG99                =    544 # Gilbert and Gill 1999 (mixed)
const  GGA_XC_HLE16               =    545 # high local exchange 2016
const  LDA_X_ERF                  =    546 # Attenuated exchange LDA (erf)
const  LDA_XC_LP_A                =    547 # Lee-Parr reparametrization A
const  LDA_XC_LP_B                =    548 # Lee-Parr reparametrization B
const  LDA_X_RAE                  =    549 # Rae self-energy corrected exchange
const  LDA_K_ZLP                  =    550 # Wigner including kinetic energy contribution
const  LDA_C_MCWEENY              =    551 # McWeeny 76
const  LDA_C_BR78                 =    552 # Brual & Rothstein 78
const  GGA_C_SCAN_E0              =    553 # GGA component of SCAN
const  LDA_C_PK09                 =    554 # Proynov and Kong 2009
const  GGA_C_GAPC                 =    555 # GapC
const  GGA_C_GAPLOC               =    556 # Gaploc
const  GGA_C_ZVPBEINT             =    557 # another spin-dependent correction to PBEint
const  GGA_C_ZVPBESOL             =    558 # another spin-dependent correction to PBEsol
const  GGA_C_TM_LYP               =    559 # Takkar and McCarthy reparametrization
const  GGA_C_TM_PBE               =    560 # Thakkar and McCarthy reparametrization
const  GGA_C_W94                  =    561 # Wilson 94 (Eq. 25)
const  MGGA_C_KCIS                =    562 # Krieger, Chen, Iafrate, and Savin
const  HYB_MGGA_XC_B0KCIS         =    563 # Hybrid based on KCIS
const  MGGA_XC_LP90               =    564 # Lee & Parr, Eq. (56)
const  GGA_C_CS1                  =    565 # A dynamical correlation functional
const  HYB_MGGA_XC_MPW1KCIS       =    566 # MPW1KCIS for barrier heights
const  HYB_MGGA_XC_MPWKCIS1K      =    567 # MPWKCIS1K for barrier heights
const  HYB_MGGA_XC_PBE1KCIS       =    568 # PBE1KCIS for binding energies
const  HYB_MGGA_XC_TPSS1KCIS      =    569 # TPSS1KCIS for thermochemistry and kinetics
const  GGA_X_B88M                 =    570 # Becke 88 reoptimized to be used with tau1
const  MGGA_C_B88                 =    571 # Meta-GGA correlation by Becke
const  HYB_GGA_XC_B5050LYP        =    572 # B5050LYP
const  LDA_C_OW_LYP               =    573 # Wigner with corresponding LYP parameters
const  LDA_C_OW                   =    574 # Optimized Wigner
const  MGGA_X_GX                  =    575 # GX functional of Loos
const  MGGA_X_PBE_GX              =    576 # PBE-GX functional of Loos
const  LDA_XC_GDSMFB              =    577 # Groth, Dornheim, Sjostrom, Malone, Foulkes, Bonitz
const  LDA_C_GK72                 =    578 # Gordon and Kim 1972
const  LDA_C_KARASIEV             =    579 # Karasiev reparameterization of Chachiyo
const  LDA_K_LP96                 =    580 # Liu-Parr kinetic
const  MGGA_X_REVSCAN             =    581 # revised SCAN
const  MGGA_C_REVSCAN             =    582 # revised SCAN
const  HYB_MGGA_X_REVSCAN0        =    583 # revised SCAN hybrid exchange (SCAN0)
const  MGGA_C_SCAN_VV10           =    584 # SCAN + VV10 correlation
const  MGGA_C_REVSCAN_VV10        =    585 # REVSCAN + VV10 correlation
const  MGGA_X_BR89_EXPLICIT       =    586 # Becke-Roussel 89
const  GGA_K_PBE3                 =    595 # Three parameter PBE-like expansion
const  GGA_K_PBE4                 =    596 # Four parameter PBE-like expansion
const  GGA_K_EXP4                 =    597 # Intermediate form between PBE3 and PBE4
const  HYB_MGGA_XC_B98            =    598 # Becke 98


text:
@testset "106.construct-binary-tree-from-inorder-and-postorder-traversal.jl" begin
    inorder = [9, 3, 15, 20, 7]
    postorder = [9, 15, 7, 20, 3]
    tree = build_tree_inpost(inorder, postorder)
    @test inorder_traversal(tree) == inorder && postorder_traversal(tree) == postorder
end

ext: jl


text:
function findbranchpoints(g)
	vs = findall(outdegree(g) .> 1)
	idx = map(vs) do v
		onb = outneighbors(g, v)
		any(outdegree(g, onb) .> 0) && return true
		return false
	end
	vs[idx]
end

function findmergepoints(g)
	vs = findall(indegree(g) .> 1)
	idx = map(vs) do v
		onb = inneighbors(g, v)
		any(indegree(g, onb) .> 0) && return true
		return false
	end
	vs[idx]
end

findjunctions(mp) = intersect(findmergepoints(mp), findbranchpoints(mp))

struct MainPathSegment{T}
    intermediates::Vector{T}
    startknot::T
    endknot::T
end

Base.length(seg::MainPathSegment) = length(seg.intermediates)

function _collect_segment(g::AbstractGraph{T}, v::T, nbfun, branchind, mergeind) where T <: Integer
    (branchind[v] || mergeind[v]) && return T[], v

    intermediates = T[]
    current = v

    while !(branchind[current] || mergeind[current])
        push!(intermediates, current)
        next = nbfun(g, current)

        if length(next) == 0
            return intermediates, 0
        elseif length(next) > 1
            append!(intermediates, next)
            return intermediates, 0
        elseif (branchind[next[1]] || mergeind[next[1]])
            return intermediates, next[1]
        end

        current = next[1]
    end
end

function segments(mp::MainPathResult)
    g = mp.mainpath
    brapoints = findbranchpoints(g)
    merpoints = findmergepoints(g)

    braind = falses(nv(g)); braind[brapoints] .= true
    merind = falses(nv(g)); merind[merpoints] .= true

    visited = falses(nv(g))
    segs = MainPathSegment[]

    for b in brapoints
        inb = inneighbors(g, b)
        onb = outneighbors(g, b)

        for i in inb
            isnothing(i) && break
            visited[i] && continue
            intermed, endpoint = _collect_segment(g, i, inneighbors, braind, merind)
            visited[intermed] .= true
            push!(segs, MainPathSegment(reverse(intermed), endpoint, b))
        end

        for o in onb
            isnothing(o) && break
            visited[o] && continue
            intermed, endpoint = _collect_segment(g, o, outneighbors, braind, merind)
            visited[intermed] .= true
            push!(segs, MainPathSegment(intermed, b, endpoint))
        end
    end

    for m in merpoints
        inb = inneighbors(g, m)
        onb = outneighbors(g, m)

        for i in inb
            isnothing(i) && break
            visited[i] && continue
            intermed, endpoint = _collect_segment(g, i, inneighbors, braind, merind)
            visited[intermed] .= true
            push!(segs, MainPathSegment(reverse(intermed), endpoint, m))
        end

        for o in onb
            isnothing(o) && break
            visited[o] && continue
            intermed, endpoint = _collect_segment(g, o, outneighbors, braind, merind)
            visited[intermed] .= true
            push!(segs, MainPathSegment(intermed, m, endpoint))
        end
    end

    segs
end

function meanweight(mp::MainPathResult, seg::MainPathSegment, w::AbstractMatrix)
    vs = seg.startknot == 0 ? seg.intermediates : vcat(seg.startknot, seg.intermediates)
    vs = seg.endknot == 0 ? vs : vcat(vs, seg.endknot)

    path, pvs = induced_subgraph(mp.mainpath, vs)
    res = 0.0
    for e in edges(path)
        s = pvs[src(e)]; d = pvs[dst(e)]
        res += w[mp.vertices[s], mp.vertices[d]]
    end

    res / ne(path)
end


text:
#=

   From BLOG example/simple-aircraft.blog

   Summary Statistics
   parameters      mean       std   naive_se      mcse        ess       rhat   ess_per_sec
       Symbol   Float64   Float64    Float64   Float64   Float64?   Float64?      Float64?

  numAircraft    1.4000    0.5831     0.0058    0.0519    57.7156     1.0166       18.9791
   numBlip[1]    2.4750    0.8212     0.0082    0.0722    61.7070     1.0233       20.2917
     numBlip1    2.4750    0.8212     0.0082    0.0722    61.7070     1.0233       20.2917
   numBlip[2]    1.3214    0.6577     0.0111    0.1033    missing    missing       missing
   numBlip[3]    1.2500    0.4334     0.0194    0.1936    missing    missing       missing

  Distributions of variable numAircraft (num:0)
  1.00000 =>    6500  (0.650000)
  2.00000 =>    3000  (0.300000)
  3.00000 =>     500  (0.050000)

  Distributions of variable numBlip1 (num:0)
  3.00000 =>    6750  (0.675000)
  1.00000 =>    1750  (0.175000)
  2.00000 =>    1375  (0.137500)
  0.00000 =>     125  (0.012500)

  Cf ~/webppl/simple_aircraft.wppl

=#
using Turing, StatsPlots, Distributions
include("jl_utils.jl")

@model function simple_aircraft(sumBlip=3)
    numAircraft ~ Poisson(5)

    # each Aircraft has a blip
    numBlip ~ filldist(Poisson(4), numAircraft)

    # We observe that there are in all 3 blips.
    sumBlip ~ Dirac(sum(numBlip))

    numBlip1 ~ Dirac(numBlip[1])

end

model = simple_aircraft()

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(5), 10)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns, :numAircraft)
show_var_dist_pct(chns, :sumBlip)
show_var_dist_pct(chns, :numBlip1)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using GaPLAC
using AbstractGPs
using KernelFunctions
using Distributions
using CSV
using DataFrames
using Turing

gpspec = GaPLAC.gp_formula("y ~| SExp(:x)")

df = CSV.read("data.tsv", DataFrame)
eq, vars = GaPLAC._apply_vars(gp_form)

kernels = GaPLAC._walk_kernel(eq)
inferable = [:x]

y = df[!, resp]
x = Matrix(df[!, vars])

@model function inference_engine(Y, X, eq, inferable)
    ℓ ~ TruncatedNormal(1, 5, 0, 30)
    gp, = GaPLAC._apply_vars(eq; hyperparams=Dict(v=> ℓ for v in inferable))

    fx ~ AbstractGPs.FiniteGP(GP(gp), RowVecs(X), 0.1)
    Y .~ Normal.(fx, 1)
end

m = inference_engine(y, x, gp_form, inferable)

chain = sample(m, NUTS(0.65), 200, progress=true)

max_stars_repo_name: biobakery/GaPLAC

max_issues_repo_name: biobakery/GaPLAC

max_forks_repo_name: biobakery/GaPLAC


text:
using Test: @testset, @test, @test_throws, @test_broken
using AlignedArrays

@testset "AlignedArrays" begin
	@testset "AlignedArrays" begin
		a = AlignedVector{Int, 256}(undef, 3)
		@test eltype(a) === Int
		@test length(a) === 3
		@test reinterpret(Int, pointer(a)) % 256 == 0

		a[1] = 1234
		@test a[1] == 1234
		a .= zeros(Int, 3)
		@test a[1] == a[2] == a[3] == 0
	end

	@testset "PageAlignedArrays" begin
		a = PageAlignedVector{Int}(undef, 3)
		@test eltype(a) === Int
		@test length(a) === 3
		@test reinterpret(Int, pointer(a)) % AlignedArrays.PAGESIZE == 0

		a[1] = 1234
		@test a[1] == 1234
		a .= zeros(Int, 3)
		@test a[1] == a[2] == a[3] == 0
	end
end


text:
module EMLratio
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...types.EMLtypes
	import EMLtypes.length
	include("ratio/Ratio.jl")
end


text:
module Viz
include("./observable.jl")

using PyPlot

const cm = PyPlot.cm.get_cmap("tab20")
options = [
    Dict(
        "divided_by" => 1.0,
        "xlim" => (),
        "xticks" => nothing,
        "xlabel" => "Time",
        "ylim" => (),
        "yticks" => nothing,
        "ylabel" => replace(replace(observables[i], "__" => "\n"), "_" => " "),
        "exp_data" => true,
        "legend_loc" => nothing,
        "cmap" => [cm.colors[j] for j in 1:20],
        "shape" => ["o", "v", "^", "<", ">", "8", "s",
                    "p", "*", "h", "H", "D", "d", "P", "X"],
        "dont_show" => [],
    ) for i in 1:length(observables)]
# ---
options[observables_index("Nuclear_NFkB")]["xlim"] = (0, 360)
options[observables_index("Nuclear_NFkB")]["xticks"] = [0, 120, 240,  360]
options[observables_index("Nuclear_NFkB")]["xlabel"] = "Time (min)"
options[observables_index("Nuclear_NFkB")]["ylim"] = (0.0, 1.1)
options[observables_index("Nuclear_NFkB")]["yticks"] = [0.0, 0.5, 1.0]
options[observables_index("Nuclear_NFkB")]["ylabel"] = "NFκB activity"

function set_rcParams()
    rc("figure",figsize = (4,3))
    rc("font",family = "Arial")
    rc("font",size = 20)
    rc("axes",linewidth = 1.5)
    rc("xtick.major",width = 1.5)
    rc("ytick.major",width = 1.5)
    rc("lines",linewidth = 1.8)
    rc("lines",markersize = 8)
end

# ---
end # module

max_stars_repo_name: himoto/paramEstim

max_issues_repo_name: himoto/paramEstim

max_forks_repo_name: himoto/paramEstim


text:
# Geometric.
geo(p::Float64) = rand(:flip, Bernoulli(p)) ? 0 : 1 + rand(:geo, geo, p)

@testset "Trace" begin
    ret, cl = simulate(geo, 0.5)
    @test haskey(cl.trace, :flip)
end

@testset "Generate" begin
    sel = selection([(:flip, ) => true])
    ret, cl, w = generate(sel, geo, 0.5)
    @test haskey(cl.trace, :flip)
    @test get_ret(cl[:flip]) == true
end

@testset "Update" begin
    ret, cl = simulate(geo, 0.5)
    sel = selection([(:flip, ) => true])
    ret, cl, w, retdiff, d = update(sel, cl)
    @test get_ret(cl[:flip]) == true
    ret, cl, w, retdiff, d = update(sel, cl, UndefinedChange(), 0.1)
    @test get_ret(cl[:flip]) == true
end


text:
push!(LOAD_PATH, "../src/serial")
using Pairwise, Multiagent

# visualize pairwise encounter
d = get_pomdp()
gen_pairwise_policy(d, 1.0, "example")
viz_pairwise_policy(d)

# load Q(s, a) table and grid
alpha = read_alpha()
g = get_grid(d.pomdp.states)

# visualize policy for simple 3-uavs scenario
viz_policy(alpha, g)

# simple stress test
uavs = randuavs(3)
trajs, nlms, nlmsb, time = stress!(utm, maxmin, true, NT, uavs, alpha, g)
trajs, nlms, nlmsb, time = stress!(uncrd, maxmin, true, NT, uavs, alpha, g)
trajs, nlms, nlmsb, time = stress!(coord, maxsum, true, NT, uavs, alpha, g)
trajs, nlms, nlmsb, time = stress!(naive, maxmin, true, NT, uavs, alpha, g)
trajs, nlms, nlmsb, time = stress!(distr, maxsum, true, NT, uavs, alpha, g)
plot_trajs(trajs)

# bulk stress test
nuavs = 2:5
nsim = 100
bulk_test(nuavs, nsim, alpha, g)


text:
using GAP_pkg_unipot
using Test

import GAP_pkg_unipot
println("TODO")


text:
#!/usr/bin/env julia

using SpectralElemenets, OrdinaryDiffEq, Flux
using LinearAlgebra,Plots
#----------------------------------#
Ex = 2; nr1 = 8;
Ey = 2; ns1 = 8;

m1 = Mesh(nr1,ns1,Ex,Ey,[true,true])
bc = ['D','D','D','D']
M  = generateMask(bc,m1)
opM(u) = m1.Bi .* u
opB(u) = gatherScatter(mass(u,m1),m1)
#----------------------------------#
kx=1.0
ky=1.0
kt=1.0
ν =1.0

utrue(x,y,t) = @. sin(kx*pi*x)*sin.(ky*pi*y)*cos(kt*pi*t)
setIC(x,y)   = utrue(x,y,0.0)

function setForcing(x,y,t)
    ut = utrue(x,y,t)
    f  = ut*((kx^2+ky^2)*pi^2*ν)
    f -= sin(kx*pi*x)*sin(ky*pi*y)*sin(kt*pi*t)*(kt*pi)
#   f  = zero(x)
end

#----------------------------------#
it = 0
function cond(u,t,integrator)
    global it
    cond = (it % 100) == 0
    it += 1
    return true # cond
end

function affect!(integrator)
    global it
    u = integrator.u
    p = integrator.p
    t = integrator.t
    dt = integrator.dt

    ut = utrue(m1.x,m1.y,t)
    er = norm(ut-u,Inf)
    println("Step=$it, Time=$t, dt=$dt, er=$er")
    plt = meshplt(u,m1)
    plt = plot!(zlims=(-1,1))
    display(plt)
    return
end

cb = DiscreteCallback(cond,affect!,save_positions=(false,false))

#----------------------------------#
function f1(u,p,t) # implicit
    f = setForcing.(m1.x,m1.y,t)

    rhs = -lapl(u,m1)
    rhs = gatherScatter(rhs,m1)

    dudt = pcg(rhs,opB;opM=opM,mult=m1.mult,ifv=false)

    # do inversion here

    return dudt
end

function f2(u,p,t) # explicit
    f = setForcing.(m1.x,m1.y,t)

    rhs = mass(f,m1)
    rhs = gatherScatter(rhs,m1)

    dudt = pcg(rhs,opB;opM=opM,mult=m1.mult,ifv=false)

    return dudt
end

u0  = setIC(m1.x,m1.y)
dt  = 0.01
tsp = (0.0,1.0)

prob = SplitODEProblem(f1,f2,u0,tsp)
sol  = solve(prob,SBDF3();dt=dt,saveat=0.1,callback=cb)

err = sol.u - [utrue.(m1.x,m1.y,sol.t[i]) for i=axes(sol.t,1)]
ee = maximum.(abs.(err[i]) for i=axes(sol.t,1))
#----------------------------------#
return ee


text:
"""
    make_api_url(root::String)
    make_api_url(root::String, apiversion::Int)

Produce a GETable api url. Root is the base url, apiversion is the version of
the api to use. If apiversion is not specified, the default behavior when a
querry is performed is to use the latest version.
"""
function make_api_url(root::String)
    apiversion = ""
    make_api_url(root, apiversion)
end

function make_api_url(root::String, apiversion::Int)
    @assert apiversion > 0 "API version must be larger than 0"
    apiversion = string(apiversion)
    make_api_url(root, apiversion)
end

function make_api_url(root::String, apiversion::String)
    if apiversion != ""
        apiversion = "/" * apiversion
    end
    root * "/api" * apiversion * "/"
end

"""
    shorten_url(url::String, squeeze::Int=20)

Shorten the url for printing.
"""
function shorten_url(url::String, squeeze::Int=20)
    url[1:squeeze] * "..." * url[(end-squeeze):end]
end

"""
    test_url(url::String)

Test url validity.
"""
function test_url(url::String)
    @info "Testing url $url"
    status = http_get(url).status
    if status == 200
        @info "URL is valid [status: $status]"
        return true
    else
        @info "URL is invalid [status: $status]"
        return false
    end
end


text:
module TestSimulation
using GrandGameOfLife
using Test

@testset "flashlight" begin
    state = to_cell.([(10, 10), (9, 10), (11, 10)])
    state = LifeCells(state)

    newstate = makestep(state)
    @test !has(newstate, to_cell(9, 10))
    @test !has(newstate, to_cell(11, 10))
    @test has(newstate, to_cell(10, 10))
    @test has(newstate, to_cell(10, 11))
    @test has(newstate, to_cell(10, 9))
end

end # module


text:
"
  Name     : c4_04_PyPlot_julia.jl
  Book     : Hands-on Data Science with Anaconda )
  Publisher: Packt Publishing Ltd.
  Author   : Yuxing Yan and James Yan
  Date     : 1/25/2018
  email    : yany@canisius.edu
             paulyxy@hotmail.com
"

# Pkg.add("Plots")

using PyPlot
x = linspace(0,2*pi,1000)
y = sin(2*x + 3*cos.(1.5*x));
plot(x, y, color="green", linewidth=2.0, linestyle="--")


text:
using SymPy

function binary_code( f::Symbol )
    x = Sym( "x" )
    y = Sym( "y" )
    ux = Sym( "u(x)" )
    uy = Sym( "u(y)" )
    dx = Sym( "d(x)" )
    dy = Sym( "d(y)" )
    ex = eval( quote( $f( x, y ) end ))
    code1 = string( f, "{Ux,Dx,Uy,Dy}( x::Dual{Ux,Dx}, y::Dual{Uy,Dy} ) = Dual",
                   (replace( replace( ex, x, ux ), y, uy ),
                    replace( replace( dx*diff( ex, x )+dy*diff( ex, y ), x, ux), y, uy)))
    code2 = string( f, "{Ux,Dx,Y}( x::Dual{Ux,Dx}, y::Y ) = Dual",
                   (replace( replace( ex, x, ux ), y, uy ),
                    replace( replace( dx*diff( ex, x )+dy*diff( ex, y ), x, ux), y, uy)))
    code3 = string( f, "{X,Uy,Dy}( x::X, y::Dual{Uy,Dy} ) = Dual",
                   (replace( ex, y, uy ),
                    replace( replace( dx*diff( ex, x )+dy*diff( ex, y ), x, ux), y, uy)))

end

function binary_code( f::Symbol )
    x = Sym( "x" )
    y = Sym( "y" )
    ux = Sym( "u(x)" )
    dx = Sym( "d(x)" )
    ex = eval( quote( $f( x, y ) end ))
    string( f, "{Ux,Dx,Y}( x::Dual{Ux,Dx}, y::Y ) = Dual",
           (replace( replace( ex, x, ux ), y, uy ),
            replace( replace( dx*diff( ex, x )+dy*diff( ex, y ), x, ux), y, uy)))
end

for op in (:+, :-, :*, :/)
    binary_code( op )
end

function unary_code( f::Symbol )
    x = Sym( "x" )
    u = Sym( "u(x)" )
    d = Sym( "d(x)" )
    ex = eval( quote $f( x ) end )
    string( f, "{U,D}(x::Dual{U,D}) = Dual",

end

for specfun in (:sin, :cos, :tan, :cot, :sec, :csc)
    @eval unary_code( specfun )
end

for specfun in (:asin :acos, :atan, :acot, :asec, :acsc)
    @eval unary_code( specfun )
end

max_stars_repo_name: sabjohnso/AutoDiff

max_issues_repo_name: sabjohnso/AutoDiff

max_forks_repo_name: sabjohnso/AutoDiff


text:
# Bryson Denham

## Packages that will be used
using NLOptControl, PrettyPlots

n=define(numStates=2,numControls=1,X0=[0.,1],XF=[0.,-1.],XL=[0.,NaN],XU=[1/9,NaN]);

dx=[:(x2[j]),:(u1[j])]
dynamics!(n,dx)

configure!(n;(:integrationScheme=>:bkwEuler),(:finalTimeDV=>true));

obj=integrate!(n,:(0.5*u1[j]^2));
@NLobjective(n.ocp.mdl,Min,obj);

optimize!(n);


text:
##########
# TArray #
##########

"""
    TArray{T}(dims, ...)

Implementation of data structures that automatically perform copy-on-write after task copying.

If current_task is an existing key in `s`, then return `s[current_task]`. Otherwise, return `s[current_task] = s[last_task]`.

Usage:

```julia
TArray(dim)
```

Example:

```julia
ta = TArray(4)              # init
for i in 1:4 ta[i] = i end  # assign
Array(ta)                   # convert to 4-element Array{Int64,1}: [1, 2, 3, 4]
```
"""
struct TArray{T,N} <: DenseArray{T,N}
  ref :: Symbol  # object_id
  TArray{T,N}() where {T,N} = new(gensym())
end

TArray{T}() where T = TArray(T,  d)
TArray{T,1}(d::Integer) where T = TArray(T,  d)
TArray{T}(d::Integer...) where T = TArray(T, convert(Tuple{Vararg{Int}}, d))
TArray{T,N}(d::Integer...) where {T,N} = length(d)==N ? TArray(T, convert(Tuple{Vararg{Int}}, d)) : error("malformed dims")
TArray{T,N}(dim::NTuple{N,Int}) where {T,N} = TArray(T, dim)

function TArray(T::Type, dim)
  res = TArray{T,length(dim)}();
  n = n_copies()
  d = Array{T}(undef, dim)
  task_local_storage(res.ref, (n,d))
  res
end

# pass through getindex and setindex!
# duplicate TArray if task id does not match current_task
function Base.getindex(S::TArray, i::Real)
  t, d = task_local_storage(S.ref)
  newd = d
#   ct = current_task()
#   if t != ct
#     # println("[getindex]: $(S.ref ) copying data")
#     newd = deepcopy(d)
#     task_local_storage(S.ref, (ct, newd))
#   end
  getindex(newd, i)
end

function Base.setindex!(S::TArray, x, i::Real)
  n, d = task_local_storage(S.ref)
  cn   = n_copies()
  newd = d
  if cn > n
    # println("[setindex!]: $(S.ref) copying data")
    newd = deepcopy(d)
    task_local_storage(S.ref, (cn, newd))
  end
  setindex!(newd, x, i)
end

function Base.push!(S::TArray, x)
  n, d = task_local_storage(S.ref)
  cn   = n_copies()
  newd = d
  if cn > n
    newd = deepcopy(d)
    task_local_storage(S.ref, (cn, newd))
  end
  push!(newd, x)
end

function Base.pop!(S::TArray)
  n, d = task_local_storage(S.ref)
  cn   = n_copies()
  newd = d
  if cn > n
    newd = deepcopy(d)
    task_local_storage(S.ref, (cn, newd))
  end
  pop!(d)
end

function Base.convert(::Type{TArray}, x::Array)
  res = TArray{typeof(x[1]),ndims(x)}();
  n   = n_copies()
  task_local_storage(res.ref, (n,x))
  res
end

function Base.convert(::Array, x::Type{TArray})
  n,d = task_local_storage(S.ref)
  c = deepcopy(d)
  return c
end

Base.show(io::IO, S::TArray) = Base.show(io::IO, task_local_storage(S.ref)[2])
Base.size(S::TArray) = Base.size(task_local_storage(S.ref)[2])
Base.ndims(S::TArray) = Base.ndims(task_local_storage(S.ref)[2])

# Base.get(t::Task, S) = S
# Base.get(t::Task, S::TArray) = (t.storage[S.ref][2])
Base.get(S::TArray) = (current_task().storage[S.ref][2])

##########
# tzeros #
##########

"""
     tzeros(dims, ...)

Construct a distributed array of zeros.
Trailing arguments are the same as those accepted by `TArray`.

```julia
tzeros(dim)
```

Example:

```julia
tz = tzeros(4)              # construct
Array(tz)                   # convert to 4-element Array{Int64,1}: [0, 0, 0, 0]
```
"""
function tzeros(T::Type, dim)
  res = TArray{T,length(dim)}();
  n = n_copies()
  d = zeros(T,dim)
  task_local_storage(res.ref, (n,d))
  res
end

tzeros(::Type{T}, d1::Integer, drest::Integer...) where T = tzeros(T, convert(Dims, tuple(d1, drest...)))
tzeros(d1::Integer, drest::Integer...) = tzeros(Float64, convert(Dims, tuple(d1, drest...)))
tzeros(d::Dims) = tzeros(Float64, d)


text:
using BenchmarkTools

lines = replace.(split(open(f -> read(f, String), "4/input.txt"), "\n\n"), "\n" => " ");

required_fields = ["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"];

function convert_line_to_record(line)
    parse = reduce(hcat, split.(string.(split(rstrip(line), " ")), ":"));
    return [parse[1,:] parse[2,:]]
end

function validate_required_fields(field_list, required_fields)
    return isempty(setdiff(required_fields, field_list)) ? 1 : 0
end

function validate_byr(value)
    return length(value) == 4 && (1920 <= parse(Int, value) <= 2002);
end

function validate_iyr(value)
    return length(value) == 4 && (2010 <= parse(Int, value) <= 2020);
end

function validate_eyr(value)
    return length(value) == 4 && (2020 <= parse(Int, value) <= 2030);
end

function validate_hgt(value)
    if value[end - 1:end] == "cm"
        return 150 <= parse(Int, value[1:end - 2]) <= 193
    elseif value[end - 1:end] == "in"
        return 59 <= parse(Int, value[1:end - 2]) <= 76
    else
        return false;
    end
end

function validate_hcl(value)
    return match(r"#[a-fA-F\d]{6}", value) !== nothing;
end

function validate_ecl(value)
    return value in ["amb" "blu" "brn" "gry" "grn" "hzl" "oth"];
end

function validate_pid(value)
    return length(value) == 9 && tryparse(Int, value) !== nothing;
end

function validate_passport(record)
    if validate_required_fields(record[:,1], required_fields) == 1
        for (item, value) in zip(record[:,1], record[:,2])
            if item == "byr"
                if !validate_byr(value)
                    return 0
                end
            elseif item == "iyr"
                if !validate_iyr(value)
                    return 0
                end
            elseif item == "eyr"
                if !validate_eyr(value)
                    return 0
                end
            elseif item == "hgt"
                if !validate_hgt(value)
                    return 0
                end
            elseif item == "hcl"
                if !validate_hcl(value)
                    return 0
                end
            elseif item == "ecl"
                if !validate_ecl(value)
                    return 0
                end
            elseif item == "pid"
                if !validate_pid(value)
                    return 0
                end
            end
        end
        return 1
    else
        return 0
    end
end

function part1(input)
    count = 0;
    for line in input
        field_list = convert_line_to_record(line)[:,1];
        count += validate_required_fields(field_list, required_fields);
    end
    return count;
end

function part2(input)
    count = 0;
    for line in input
        record = convert_line_to_record(line);
        count += validate_passport(record);
    end
    return count;
end

println("Part 1: $(part1(lines))");
println("Part 2: $(part2(lines))");

@benchmark part2(lines)


text:
module SlimPlotting

using PyPlot, Statistics, ColorSchemes

const cc = PyPlot.PyNULL()

function __init__()
    try
        copy!(cc, PyPlot.pyimport_conda("colorcet", "colorcet"))
    catch e
        # Not using julia's conda and not installed. Installing and loading
        run(Cmd([PyPlot.PyCall.pyprogramname, "-m", "pip", "install", "-U", "--user", "colorcet"]))
        copy!(cc, PyPlot.pyimport("colorcet"))
    end
end

export plot_fslice, plot_velocity, plot_simage, plot_sdata, wiggle_plot
export colorschemes

"""
    _plot_with_units(image, spacing; perc=95, cmap=:cet_CET_L1,
                     o=(0, 0), interp="hanning", aspect=nothing, d_scale=0,
                     positive=false, labels=(:X, :Depth), cbar=false,
                     units=(:m, :m), name="RTM", new_fig=true, save=nothing)

Plot a 2D grided image with physical units defined by the grid spacing `spacing`.

# Arguments
  - `image::Array{T, 2}`: image to be plotted
  - `spacing::Tuple`: grid spacing in physical units
  - `perc::Int`: (Optional) Clipping percentile, default=95
  - `cmap::Symbol`: (Optional) Color map, default=:linear_grey_10_95_c0_n256
  - `o::Tuple`: (Optional) Origin of the image, default=(0, 0)
  - `interp::String`: (Optional) Interpolation method, default="hanning"
  - `aspect::Symbol`: (Optional) Aspect ratio, default=:auto
  - `d_scale::Float`: (Optional) Depth scaling, default=1.5. Applied scaling is `(1:max_depth).^d_scale`.
  - `positive::Bool`: (Optional) Plot positive only image (clip `[0:max(image)]`), default=false
  - `labels::Tuple`: (Optional) Labels for the axes, default=(:X, :Depth)
  - `name::String`: (Optional) Figure title, default="RTM"
  - `units::Tuple(String)`: (Optional) Physical units of each axis, default=(:m, :m).
  - `new_fig::Bool`: (Optional) Create a new figure, default=true
  - `save::String`: (Optional) Save figure to file, default=nothing doesn't save the figure
  - `cbar::Bool`: (Optional) Show colorbar, default=false
"""
function _plot_with_units(image, spacing; perc=95, cmap=:cet_CET_L1, vmin=nothing, vmax=nothing,
                          o=(0, 0), interp="hanning", aspect=nothing, d_scale=0,
                          positive=false, labels=(:X, :Depth), cbar=false,
                          units=(:m, :m), name="RTM", new_fig=true, save=nothing, ax=nothing)
    nz, nx = size(image)
    dz, dx = spacing
    oz, ox = o
    depth = range(oz, oz + (nz - 1)*spacing[2], length=nz).^d_scale
    scaled = image .* depth

    a = positive ? maximum(scaled) : quantile(abs.(vec(scaled)), perc/100)
    isnothing(vmax) || (a = vmax)
    ma = positive ? minimum(scaled) : -a
    isnothing(vmin) || (ma = vmin)
    extent = [ox, ox+ (nx-1)*dx, oz+(nz-1)*dz, oz]
    isnothing(aspect) && (aspect = :auto)

    # color map
    cmap = try ColorMap(cmap); catch; ColorMap(colorschemes[cmap].colors); end
    new_fig && figure()
    # Plot
    if isnothing(ax)
        imshow(scaled, vmin=ma, vmax=a, cmap=cmap, aspect=aspect, interpolation=interp, extent=extent)
    else
        ax.imshow(scaled, vmin=ma, vmax=a, cmap=cmap, aspect=aspect, interpolation=interp, extent=extent)
    end
    xlabel("$(labels[1]) [$(units[1])]")
    ylabel("$(labels[2]) [$(units[2])]")
    title("$name")
    cbar && colorbar(fraction=0.046, pad=0.04)

    if ~isnothing(save)
    save == true ? filename=name : filename=save
    savefig(filename, bbox_inches="tight", dpi=150)
    end
end

"""
    plot_simage(image, spacing; perc=98, cmap=:linear_grey_10_95_c0_n256,
                o=(0, 0), interp="hanning", aspect=nothing, d_scale=1.5,
                labels=(:X, :Depth), name="RTM", units=(:m, :m), new_fig=true,
                save=nothing, cbar=false)

Plot a 2D seismic image with a grid spacing `spacing`. Calls [`_plot_with_units`](@ref).

# Arguments
  - `image::Array{T, 2}`: image to be plotted
  - `spacing::Tuple`: grid spacing in physical units
  - `perc::Int`: (Optional) Clipping percentile, default=95
  - `cmap::Symbol`: (Optional) Color map, default=:linear_grey_10_95_c0_n256
  - `o::Tuple`: (Optional) Origin of the image, default=(0, 0)
  - `interp::String`: (Optional) Interpolation method, default="hanning"
  - `aspect::Symbol`: (Optional) Aspect ratio, default=:auto
  - `d_scale::Float`: (Optional) Depth scaling, default=1.5. Applied scaling is `(1:max_depth).^d_scale`.
  - `labels::Tuple`: (Optional) Labels for the axes, default=(:X, :Depth)
  - `name::String`: (Optional) Figure title, default="RTM"
  - `units::Tuple(String)`: (Optional) Physical units of each axis, default=(:m, :m).
  - `new_fig::Bool`: (Optional) Create a new figure, default=true
  - `save::String`: (Optional) Save figure to file, default=nothing doesn't save the figure
  - `cbar::Bool`: (Optional) Show colorbar, default=false

"""
function plot_simage(image; kw...)
    d = try image.d; catch nothing; end
    if isnothing(d)
        @warn "No grid spacing specified, plotting with a 1m grid spacing"
        d = (1, 1)
    end
    plot_simage(image.data, d; kw...)
end

"""
    plot_fslice(image, spacing; perc=98, cmap=:diverging_bwr_20_95_c54_n256,
                o=(0, 0), interp="hanning", aspect=nothing, d_scale=1.5,
                name="RTM", units="m", new_fig=true, save=nothing)

Plot a 2D frequency slice of seismic data. Calls [`_plot_with_units`](@ref).

# Arguments
  - `image::Array{T, 2}`: image to be plotted
  - `spacing::Tuple`: grid spacing in physical units
  - `perc::Int`: (Optional) Clipping percentile, default=95
  - `cmap::Symbol`: (Optional) Color map, default=:diverging_bwr_20_95_c54_n256
  - `o::Tuple`: (Optional) Origin of the image, default=(0, 0)
  - `interp::String`: (Optional) Interpolation method, default="hanning"
  - `aspect::Symbol`: (Optional) Aspect ratio, default=:auto
  - `d_scale::Float`: (Optional) Depth scaling, default=1.5. Applied scaling is `(1:max_depth).^d_scale`.
  - `labels::Tuple`: (Optional) Labels for the axes, default=(:X, :Depth)
  - `name::String`: (Optional) Figure title, default="RTM"
  - `units::Tuple(String)`: (Optional) Physical units of each axis, default=(:m, :m).
  - `new_fig::Bool`: (Optional) Create a new figure, default=true
  - `save::String`: (Optional) Save figure to file, default=nothing doesn't save the figure
  - `cbar::Bool`: (Optional) Show colorbar, default=false

"""
function plot_fslice(image; kw...)
    @warn "No grid spacing specified, plotting with a 1m grid spacing"
    d = (1, 1)
    plot_fslice(image.data, d; kw...)
end

plot_fslice(image::AbstractArray{T}, args...; kw...) where {T<:Complex} = plot_fslice(real.(image), args...; kw...)

"""
    plot_velocity(image, spacing; perc=98, cmap=:cet_rainbow,
                o=(0, 0), interp="hanning", aspect=nothing, d_scale=1.5,
                name="RTM", units="m", new_fig=true, save=nothing)

Plot a velocity model. Calls [`_plot_with_units`](@ref).

# Arguments
  - `image::Array{T, 2}`: image to be plotted
  - `spacing::Tuple`: grid spacing in physical units
  - `perc::Int`: (Optional) Clipping percentile, default=95
  - `cmap::Symbol`: (Optional) Color map, default=:cet_rainbow
  - `o::Tuple`: (Optional) Origin of the image, default=(0, 0)
  - `interp::String`: (Optional) Interpolation method, default="hanning"
  - `aspect::Symbol`: (Optional) Aspect ratio, default=:auto
  - `d_scale::Float`: (Optional) Depth scaling, default=1.5. Applied scaling is `(1:max_depth).^d_scale`.
  - `labels::Tuple`: (Optional) Labels for the axes, default=(:X, :Depth)
  - `name::String`: (Optional) Figure title, default="RTM"
  - `units::Tuple(String)`: (Optional) Physical units of each axis, default=(:m, :m).
  - `new_fig::Bool`: (Optional) Create a new figure, default=true
  - `save::String`: (Optional) Save figure to file, default=nothing doesn't save the figure
  - `cbar::Bool`: (Optional) Show colorbar, default=false

"""
function plot_velocity(image; kw...)
    d = try image.d; catch nothing; end
    if isnothing(d)
        @warn "No grid spacing specified, plotting with a 1m grid spacing"
        d = (1, 1)
    end
    plot_velocity(image.data, d; kw...)
end

"""
    plot_sdata(image, spacing; perc=98, cmap=:linear_grey_10_95_c0_n256,
                o=(0, 0), interp="hanning", aspect=nothing, d_scale=1.5,
                name="RTM", units="m", new_fig=true, save=nothing)

Plot seismic data gather (i.e shot record). Calls [`_plot_with_units`](@ref).

# Arguments
  - `image::Array{T, 2}`: image to be plotted
  - `spacing::Tuple`: grid spacing in physical units
  - `perc::Int`: (Optional) Clipping percentile, default=95
  - `cmap::Symbol`: (Optional) Color map, default=:linear_grey_10_95_c0_n256
  - `o::Tuple`: (Optional) Origin of the image, default=(0, 0)
  - `interp::String`: (Optional) Interpolation method, default="hanning"
  - `aspect::Symbol`: (Optional) Aspect ratio, default=:auto
  - `d_scale::Float`: (Optional) Depth scaling, default=1.5. Applied scaling is `(1:max_depth).^d_scale`.
  - `labels::Tuple`: (Optional) Labels for the axes, default=(:X, :Depth)
  - `name::String`: (Optional) Figure title, default="RTM"
  - `units::Tuple(String)`: (Optional) Physical units of each axis, default=(:m, :m).
  - `new_fig::Bool`: (Optional) Create a new figure, default=true
  - `save::String`: (Optional) Save figure to file, default=nothing doesn't save the figure
  - `cbar::Bool`: (Optional) Show colorbar, default=false

"""
function plot_sdata(image; kw...)
    dt = try image.dt; catch nothing; end
    dx = try
        geom = try image.geometry; catch nothing; end
        geom = hasproperty(geom, :xloc) ? geom : Geometry(geom)
        diff(geom.xloc[1])[1]
    catch
        nothing
    end
    if isnothing(dt)
        @warn "No grid spacing specified, plotting with a 1m grid spacing"
        d = (1, 1)
    else
        d = (dt, dx)
    end
    plot_sdata(image.data, d; kw...)
end

# Defaults
_funcs = [:plot_simage, :plot_fslice, :plot_velocity, :plot_sdata]
_default_colors = [:cet_CET_L1, :cet_CET_D1A, :cet_rainbow4, :cet_CET_L1]
_names = ["RTM", "Frequency slice", "Velocity", "Shot record"]
_units = [(:m, :m), (:m, :m), (:m, :m), (:m, :s)]
_labels = [(:X, :Depth), (:Xsrc, :Xrec), (:X, :Depth), (:Xrec, :T)]

for (func, cmap, pname, u ,l) ∈ zip(_funcs, _default_colors, _names, _units, _labels)
    @eval begin
        function $func(image, spacing; kw...)
            cmap = haskey(kw, :cmap) ? kw[:cmap] : $(Meta.quot(cmap))
            pname = haskey(kw, :name) ? kw[:name] : $(Meta.quot(pname))
            u = haskey(kw, :units) ? [:units] : $(Meta.quot(u))
            l = haskey(kw, :labels) ? [:labels] : $(Meta.quot(l))
            positive = $func == plot_velocity
            _plot_with_units(image, spacing; kw...,positive=positive, cmap=cmap, name=pname, units=u, labels=l)
        end
    end
end

"""
    wiggle_plot(image, xrec, time_axis; t_scale=1.5,
                new_fig=true)

wiggle_plot of a seismic traces.

# Arguments
  - `image::Array{T, 2}`: Shot record to be plotted
  - `xrec::Array{T, 1}`: Receiver coordinates
  - `time_axis::Array{T, 1}`: Time axis
  - `t_scale::Float`: (Optional) Time scaling, default=1.5. Applied scaling is `(1:max_time).^t_scale`.
  - `new_fig::Bool`: (Optional) Create a new figure, default=true
"""
function wiggle_plot(data::Array{Td, 2}, xrec=nothing, time_axis=nothing;
                     t_scale=1.5, new_fig=true) where Td
    # X axis
    if isnothing(xrec)
        @info "No X coordinates prvided, using 1:ntrace"
        xrec = range(0, size(data, 2), length=size(data, 2))
    end
    length(xrec) == size(data, 2) || error("xrec must be the same length as the number of columns in data");
    dx = diff(xrec); dx = 2 .* vcat(dx[1], dx)
    # time axis
    if isnothing(time_axis)
        @info "No time axis provided, using 1:ntime"
        time_axis = range(0, size(data, 1), length=size(data, 1))
    end
    length(time_axis) == size(data, 1) || error("time_axis must be the same length as the number of rows in data");
    # Time gain
    tg = time_axis .^ t_scale;
    ax = new_fig ? subplots()[2] : gca()

    ax.set_ylim(maximum(time_axis), minimum(time_axis))
    ax.set_xlim(minimum(xrec), maximum(xrec))
    for (i, xr) ∈ enumerate(xrec)
        x = tg.* data[:, i]
        x = dx[i] * x ./ maximum(x) .+ xr
        # rescale to avoid large spikes
        ax.plot(x, time_axis, "k-")
        ax.fill_betweenx(time_axis, xr, x, where=(x.>xr), color="k")
    end
    ax.set_xlabel("X")
    ax.set_ylabel("Time")
end

end # module


text:
#Photosynthesis model for given total conductance (stomatal+mesophyll), gₜ, irradiance incident on the leaf (Iᵢ),
#and Jₘₐₓ
function Farquhar(gₜ::T, Iᵢ::T, Jₘₐₓ::T,photo::PhotoPar,env::EnvironmentStruct) where {T<:Float64}
    Pₜₒₜ,cₐ = env.P,env.Cₐ

    gₜₚ = gₜ/Pₜₒₜ

    J = (photo.α*Iᵢ+Jₘₐₓ-sqrt(photo.α^2*Iᵢ^2+2*photo.α*Iᵢ*Jₘₐₓ*(1-2*photo.θ)+Jₘₐₓ^2))/(2*photo.θ)
    Aⱼ = J/4

    #Intercellular carbon dioxide concentration
    cᵢ = (gₜₚ*(cₐ-2*photo.Γ)-Aⱼ*(1-photo.b_r)+sqrt((gₜₚ*(cₐ-2*photo.Γ)-Aⱼ*(1-photo.b_r))^2+4*Aⱼ*gₜₚ*(photo.Γ+photo.b_r*photo.K)+8*gₜₚ^2*cₐ*photo.Γ))/(2*gₜₚ)

    #Carbon assimilation rate
    A = gₜₚ*(cₐ-cᵢ)

    return A, cᵢ
end

#Calculate the C assimilation (kg C year⁻¹ m⁻² leaf area)
function C_assimilation(gₜ::T,Iᵢ::T,Jmax::T,growthlength::T,model::CCPHStruct) where {T<:Float64}
    A = model.treepar.Xₜ*model.cons.M_C*Farquhar(gₜ,Iᵢ,Jmax,model.photopar,model.env)[1]*growthlength

    return A
end

#Calculate per tree canopy gross primary production
GPP(gₜ::T,Iᵢ::T,Jmax::T,LAI::T,growthlength::T,model::CCPHStruct) where {T<:Float64} =
C_assimilation(gₜ,Iᵢ,Jmax,growthlength,model)*(1-exp(-model.treepar.k*LAI))/(model.treesize.N*model.treepar.k)

#Calcualte leaf performance at crown base
function Calc_Δ_leaf(gₜ::T,Iᵢ::T,LAI::T,growthlength::T,Nₘ_f::T,Jmax::T,model::CCPHStruct) where {T<:Float64}
    Iᵢ_b = Iᵢ*exp(-model.treepar.k*LAI)
    Jmax_b = Jmax*exp(-model.treepar.k*LAI)
    A_b = C_assimilation(gₜ,Iᵢ_b,Jmax_b,growthlength,model)
    Δ_leaf = model.treepar.y*(A_b-model.treepar.rₘ*Nₘ_f*model.treepar.LMA)-model.treepar.LMA/model.treepar.Tf #Bottom leaf performance
    return Δ_leaf
end


text:
module TestSuggestItem

include("preamble.jl")

using SymSpellChecker: SuggestItem

@testset "relations" begin
    si1 = SuggestItem("foo", 1, 5)
    si2 = SuggestItem("bar", 2, 20)
    @test si1 < si2

    si3 = SuggestItem("bzr", 2, 15)
    @test si2 < si3

    si4 = SuggestItem("bbr", 2, 20)
    @test si2 < si4
end

end


text:
using Revise # reduce recompilation time
using Plots
# using Documenter
using LinearAlgebra
using SparseArrays
using BenchmarkTools
using UnPack
using StaticArrays
using DelimitedFiles
using Polyester
using MuladdMacro

push!(LOAD_PATH, "./src")
using Basis1D
using CommonUtils
# using Basis2DTri
# using UniformTriMesh
# using NodesAndModes
# using NodesAndModes.Tri
# using CommonUtils
# using Basis1D
using Basis2DQuad
using UniformQuadMesh

using SetupDG

push!(LOAD_PATH, "./examples/EntropyStableEuler.jl/src")
include("../../EntropyStableEuler.jl/src/logmean.jl")
include("../SBP_quad_data.jl")
using EntropyStableEuler
using EntropyStableEuler.Fluxes2D

@inline function pfun(rho,rhou,E)
    return (γ-1)*(E-.5*rhou^2/rho)
end

@inline function pfun(rho,rhou,rhov,E)
    return (γ-1)*(E-.5*(rhou^2+rhov^2)/rho)
end

@inline function Efun(rho,u,v,p)
    return p/(γ-1) + .5*rho*(u^2+v^2)
end

@inline function wavespeed_1D(rho,rhou,E)
    p = pfun(rho,rhou,E)
    return abs(rhou/rho) + sqrt(γ*p/rho)
end

@inline function logmean(aL,aR,logL,logR)

    # "from: Entropy stable num. approx. for the isothermal and polytropic Euler"

    da = aR-aL;
    aavg = .5*(aR+aL);
    f = da/aavg;
    v = f^2;
    if abs(f)<1e-4
        # numbers assume the specific value γ = 1.4
        return aavg*(1 + v*(-.2-v*(.0512 - v*0.026038857142857)))
    else
        return -da/(logL-logR)
    end
end

@inline function euler_fluxes_2D(rhoL,uL,vL,betaL,rhologL,betalogL,
                                 rhoR,uR,vR,betaR,rhologR,betalogR)

    rholog  = logmean(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = .5*(rhoL+rhoR)
    uavg   = .5*(uL+uR)
    vavg   = .5*(vL+vR)

    unorm = uL*uR + vL*vR
    pa    = rhoavg/(betaL+betaR)
    f4aux = rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm

    FxS1 = rholog*uavg
    FxS2 = FxS1*uavg + pa
    FxS3 = FxS1*vavg
    FxS4 = f4aux*uavg

    FyS1 = rholog*vavg
    FyS2 = FxS3
    FyS3 = FyS1*vavg + pa
    FyS4 = f4aux*vavg

    return FxS1,FxS2,FxS3,FxS4,FyS1,FyS2,FyS3,FyS4
end

@inline function euler_fluxes_2D_x(rhoL,uL,vL,betaL,rhologL,betalogL,
                                   rhoR,uR,vR,betaR,rhologR,betalogR)

    rholog  = logmean(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = .5*(rhoL+rhoR)
    uavg   = .5*(uL+uR)
    vavg   = .5*(vL+vR)

    unorm = uL*uR + vL*vR
    pa    = rhoavg/(betaL+betaR)
    f4aux = rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm

    FxS1 = rholog*uavg
    FxS2 = FxS1*uavg + pa
    FxS3 = FxS1*vavg
    FxS4 = f4aux*uavg

    return FxS1,FxS2,FxS3,FxS4
end

@inline function euler_fluxes_2D_y(rhoL,uL,vL,betaL,rhologL,betalogL,
                                   rhoR,uR,vR,betaR,rhologR,betalogR)

    rholog  = logmean(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = .5*(rhoL+rhoR)
    uavg   = .5*(uL+uR)
    vavg   = .5*(vL+vR)

    unorm = uL*uR + vL*vR
    pa    = rhoavg/(betaL+betaR)
    f4aux = rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm

    FyS1 = rholog*vavg
    FyS2 = FyS1*uavg
    FyS3 = FyS1*vavg + pa
    FyS4 = f4aux*vavg

    return FyS1,FyS2,FyS3,FyS4
end

@inline function inviscid_flux_prim(rho,u,v,p)
    E = Efun(rho,u,v,p)

    rhou  = rho*u
    rhov  = rho*v
    rhouv = rho*u*v
    Ep    = E+p

    fx1 = rhou
    fx2 = rhou*u+p
    fx3 = rhouv
    fx4 = u*Ep

    fy1 = rhov
    fy2 = rhouv
    fy3 = rhov*v+p
    fy4 = v*Ep

    return fx1,fx2,fx3,fx4,fy1,fy2,fy3,fy4
end

@inline function limiting_param(rhoL,rhouL,rhovL,EL,rhoP,rhouP,rhovP,EP)
    # L - low order, P - P_ij
    l = 1.0
    # Limit density
    if rhoL + rhoP < -TOL
        l = max((-rhoL+POSTOL)/rhoP, 0.0)
    end

    # limiting internal energy (via quadratic function)
    a = rhoP*EP-(rhouP^2+rhovP^2)/2.0
    b = rhoP*EL+rhoL*EP-rhouL*rhouP-rhovL*rhovP
    c = rhoL*EL-(rhouL^2+rhovL^2)/2.0-POSTOL

    d = 1.0/(2.0*a)
    e = b^2-4.0*a*c
    g = -b*d

    l_eps_ij = 1.0
    if e >= 0
        f = sqrt(e)
        h = f*d
        r1 = g+h
        r2 = g-h
        if r1 > TOL && r2 > TOL
            l_eps_ij = min(r1,r2)
        elseif r1 > TOL && r2 < -TOL
            l_eps_ij = r1
        elseif r2 > TOL && r1 < -TOL
            l_eps_ij = r2
        end
    end

    l = min(l,l_eps_ij)
    return l
end

function build_meshfree_sbp(rq,sq,wq,rf,sf,wf,nrJ,nsJ,α)
    # [-1,1,0], [-1,-1,sqrt(4/3)]
    equilateral_map(r,s) = (@. .5*(2*r+1*s+1), @. sqrt(3)*(1+s)/2 - 1/sqrt(3) )
    req,seq = equilateral_map(rq,sq)
    ref,sef = equilateral_map(rf,sf)
    barycentric_coords(r,s) = ((@. (1+r)/2), (@. (1+s)/2), (@. -(r+s)/2))
    λ1,λ2,λ3 = barycentric_coords(rq,sq)
    λ1f,λ2f,λ3f = barycentric_coords(rf,sf)

    Br = diagm(nrJ.*wf)
    Bs = diagm(nsJ.*wf)

    # build extrapolation matrix
    E = zeros(length(rf),length(rq))
    for i = 1:length(rf)
        # d = @. (λ1 - λ1f[i])^2 + (λ2 - λ2f[i])^2 + (λ3 - λ3f[i])^2
        d2 = @. (req-ref[i])^2 + (seq-sef[i])^2
        p = sortperm(d2)
        h2 = (wf[i]/sum(wf))*2/pi # set so that h = radius of circle with area w_i = face weight
        nnbrs = min(4,max(3,count(d2[p] .< h2))) # find 3 closest points
        p = p[1:nnbrs]
        Ei = vandermonde_2D(1,[rf[i]],[sf[i]])/vandermonde_2D(1,rq[p],sq[p])
        E[i,p] = Ei
    end
    E = Matrix(droptol!(sparse(E),1e-13))

    # build stencil
    A = spzeros(length(req),length(req))
    for i = 1:length(req)
        d2 = @. (req-req[i])^2 + (seq-seq[i])^2
        p = sortperm(d2)

        # h^2 = wq[i]/pi = radius of circle with area wq[i]
        # h2 =     (sqrt(3)/sum(wq))*wq[i]/pi
        h2 = α^2*(sqrt(3)/sum(wq))*wq[i]/pi

        nnbrs = count(d2[p] .< h2)
        nbrs = p[1:nnbrs]
        A[i,nbrs] .= one(eltype(A))
    end
    A = (A+A')
    A.nzval .= one(eltype(A)) # bool-ish

    # build graph Laplacian
    L1 = (A-diagm(diag(A))) # ignore
    L1 -= diagm(vec(sum(L1,dims=2)))

    b1r = -sum(.5*E'*Br*E,dims=2)
    b1s = -sum(.5*E'*Bs*E,dims=2)
    ψ1r = pinv(L1)*b1r
    ψ1s = pinv(L1)*b1s

    function fillQ(adj,ψ)
        Np = length(ψ)
        S = zeros(Np,Np)
        for i = 1:Np
            for j = 1:Np
                if adj[i,j] != 0
                        S[i,j] += (ψ[j]-ψ[i])
                end
            end
        end
        return S
    end

    S1r,S1s = fillQ.((A,A),(ψ1r,ψ1s))
    Qr = Matrix(droptol!(sparse(S1r + .5*E'*Br*E),1e-14))
    Qs = Matrix(droptol!(sparse(S1s + .5*E'*Bs*E),1e-14))

    return Qr,Qs,E,Br,Bs,A
end

function init_reference_tri_sbp_GQ(N, qnode_choice)
    include("SBP_quad_data.jl")
    # initialize a new reference element data struct
    rd = RefElemData()

    fv = tri_face_vertices() # set faces for triangle
    Nfaces = length(fv)
    @pack! rd = fv, Nfaces

    # Construct matrices on reference elements
    r, s = nodes(Tri(),N)
    VDM = vandermonde(Tri(),N, r, s)
    Vr, Vs = grad_vandermonde(Tri(),N, r, s)
    Dr = Vr/VDM
    Ds = Vs/VDM
    @pack! rd = r,s,VDM,Dr,Ds

    # low order interpolation nodes
    r1,s1 = nodes(Tri(),1)
    V1 = vandermonde(Tri(),1,r,s)/vandermonde(Tri(),1,r1,s1)
    @pack! rd = V1

    #Nodes on faces, and face node coordinate
    if qnode_choice == "GQ"
        r1D, w1D = gauss_quad(0,0,N)
    elseif qnode_choice == "GL" || qnode_choice == "tri_diage"
        r1D, w1D = gauss_lobatto_quad(0,0,N+1)
    end
    Nfp = length(r1D) # number of points per face
    e = ones(Nfp) # vector of all ones
    z = zeros(Nfp) # vector of all zeros
    rf = [r1D; -r1D; -e];
    sf = [-e; r1D; -r1D];
    wf = vec(repeat(w1D,3,1));
    nrJ = [z; e; -e]
    nsJ = [-e; e; z]
    @pack! rd = rf,sf,wf,nrJ,nsJ

    if qnode_choice == "GQ"
        rq,sq,wq = GQ_SBP[N];
    elseif qnode_choice == "GL"
        rq,sq,wq = GL_SBP[N];
    elseif qnode_choice == "tri_diage"
        rq,sq,wq = Tri_diage[N];
    end
    # rq,sq,wq = GQ_SBP[N]
    # rq,sq,wq = GL_SBP[N]
    # rq,sq,wq = Tri_diage[N]
    Vq = vandermonde(Tri(),N,rq,sq)/VDM
    M = Vq'*diagm(wq)*Vq
    Pq = M\(Vq'*diagm(wq))
    @pack! rd = rq,sq,wq,Vq,M,Pq

    Vf = vandermonde(Tri(),N,rf,sf)/VDM # interpolates from nodes to face nodes
    LIFT = M\(Vf'*diagm(wf)) # lift matrix used in rhs evaluation
    @pack! rd = Vf,LIFT

    # plotting nodes
    rp, sp = equi_nodes(Tri(),10)
    Vp = vandermonde(Tri(),N,rp,sp)/VDM
    @pack! rd = rp,sp,Vp

    return rd
end

const TOL = 5e-16
const POSTOL = 1e-14
const WALLPT = 1.0/6.0
const Nc = 4 # number of components
"Approximation parameters"
const N = 2
const K1D = 10
const T = 0.5#1e-4
const dt0 = 1e-3
const XLENGTH = 1.5#2.0
const CFL = 1.0
const NUM_THREADS = Threads.nthreads()

# Initial condition 2D shocktube
const γ = 1.4

"Mesh related variables"
# VX, VY, EToV = uniform_quad_mesh(Int(round(XLENGTH*K1D)),K1D)
# @. VX = (VX+1)/2*XLENGTH*5
# @. VY = (VY+1)/2*5
Kx = Int(round(XLENGTH*K1D))
Ky = 2*K1D
VX, VY, EToV = uniform_quad_mesh(Kx,Ky)
# @. VX = (VX+1)/2*XLENGTH*5
# @. VY = (VY+1)/2*5
@. VX = 15*(1+VX)/2
@. VY = 5*VY

rd = init_reference_quad(N,gauss_lobatto_quad(0,0,N))
"Initialize reference element"
r,_ = gauss_lobatto_quad(0,0,N)   # Reference nodes
VDM = vandermonde_1D(N,r)         # modal to nodal
Dr = grad_vandermonde_1D(N,r)/VDM # nodal differentiation
V1 = vandermonde_1D(1,r)/vandermonde_1D(1,[-1;1]) # nodal linear interpolation

Nq = N
rq,wq = gauss_lobatto_quad(0,0,Nq)
M = diagm(wq)
Mlump = zeros(size(M))
Mlump_inv = zeros(size(M))
for i = 1:Nq+1
    Mlump[i,i] = sum(M[i,:])
    Mlump_inv[i,i] = 1.0/Mlump[i,i]
end

# 1D operators
M1D = Mlump
M1D_inv = Mlump_inv
Q1D = M*Dr
B1D = zeros(N+1,N+1)
B1D[1,1] = -1
B1D[end,end] = 1
L = Array(spdiagm(0=>-2*ones(N+1), 1=>ones(N), -1=>ones(N)))
L[1,1] = -1
L[end,end] = -1
psi = pinv(L)*-1/2*B1D*ones(N+1)
S0 = zeros(N+1,N+1)
for i = 1:N+1
    for j = 1:N+1
        if L[i,j] != 0
            S0[i,j] = psi[j] - psi[i]
        end
    end
end
Q01D = S0+1/2*B1D

# Drop zeros
Q1D = Matrix(droptol!(sparse(Q1D),TOL))
Q01D = Matrix(droptol!(sparse(Q01D),TOL))
B1D = Matrix(droptol!(sparse(B1D),TOL))
S01D = Matrix(droptol!(sparse(S0),TOL))
S1D = Matrix(droptol!(sparse((Q1D-Q1D')/2),TOL))

# Tensor product operators
M = droptol!(sparse(kron(M1D,M1D)),TOL)
Minv = spdiagm(0 => 1 ./ diag(M))
Qr = droptol!(sparse(kron(M1D,Q1D)),TOL)
Qs = droptol!(sparse(kron(Q1D,M1D)),TOL)
Qr0 = droptol!(sparse(kron(M1D,Q01D)),TOL)
Qs0 = droptol!(sparse(kron(Q01D,M1D)),TOL)
Sr = droptol!(sparse((Qr-Qr')/2),TOL)
Ss = droptol!(sparse((Qs-Qs')/2),TOL)
S0r = droptol!(sparse(kron(M1D,S01D)),TOL)
S0s = droptol!(sparse(kron(S01D,M1D)),TOL)
Br = droptol!(sparse(Qr+Qr'),TOL)
Bs = droptol!(sparse(Qs+Qs'),TOL)
Br_halved = -sum(S0r,dims=2)
Bs_halved = -sum(S0s,dims=2)

@unpack Vf,Dr,Ds,LIFT,nrJ,nsJ,wf = rd
md = init_mesh((VX,VY),EToV,rd)
@unpack xf,yf,mapM,mapP,mapB,nxJ,nyJ,x,y = md
xb,yb = (x->x[mapB]).((xf,yf))

const K  = size(x,2)
const Nfaces = 4

# Make domain periodic
@unpack Vf = rd
@unpack xf,yf,mapM,mapP,mapB,rxJ,ryJ,sxJ,syJ,J,sJ,nxJ,nyJ = md
LX,LY = (x->maximum(x)-minimum(x)).((VX,VY)) # find lengths of domain
mapPB = build_periodic_boundary_maps(xf,yf,LX,LY,Nfaces*K,mapM,mapP,mapB)
mapP[mapB] = mapPB
@pack! md = mapP

xq = x
yq = y

nx = nxJ ./ sJ
ny = nyJ ./ sJ

E = Vf
Fmask = zeros(Int64,size(Vf,1))
for i = 1:size(Vf,1)
    tmparr = findall(abs.(E[i,:] .- 1) .< 1e-12)
    Fmask[i] = tmparr[1]
end

@inline function get_infoP(mapP,Fmask,i,k)
    gP = mapP[i,k]           # exterior global face node number
    kP = fld1(gP,Nfp)        # exterior element number
    iP = Fmask[mod1(gP,Nfp)] # exterior node number
    return iP,kP
end

function rhs_ESDG!(U,rhsU,ops,geom)
    Sr,Ss,S0r,S0s,Minv,Br_halved,Bs_halved,E = ops
    mapP,Fmask,xq,yq,rxJ,ryJ,sxJ,syJ,J,sJ,nx,ny,nrJ,nsJ,wf = geom

    FHx = zeros(Nc,Np,Np)
    FHy = zeros(Nc,Np,Np)
    FHPx = zeros(Nc,Nfp)
    FHPy = zeros(Nc,Nfp)
    Sx  = zeros(Np,Np)
    Sy  = zeros(Np,Np)
    EBx  = zeros(Np,Np)
    EBy  = zeros(Np,Np)
    rhsvec = zeros(Nc,Np)

    Br = diagm(wf.*nrJ)
    Bs = diagm(wf.*nsJ)
    for k = 1:K
        for j = 1:Np
            for i = 1:Np
                rhoM     = U[1,i,k]
                rhouM    = U[2,i,k]
                rhovM    = U[3,i,k]
                EM       = U[4,i,k]
                uM       = rhouM/rhoM
                vM       = rhovM/rhoM
                pM       = pfun(rhoM,rhouM,rhovM,EM)
                betaM    = rhoM/(2*pM)
                rhologM  = log(rhoM)
                betalogM = log(betaM)
                rhoP     = U[1,j,k]
                rhouP    = U[2,j,k]
                rhovP    = U[3,j,k]
                EP       = U[4,j,k]
                uP       = rhouP/rhoP
                vP       = rhovP/rhoP
                pP       = pfun(rhoP,rhouP,rhovP,EP)
                betaP    = rhoP/(2*pP)
                rhologP  = log(rhoP)
                betalogP = log(betaP)

                Fx1,Fx2,Fx3,Fx4,Fy1,Fy2,Fy3,Fy4 = euler_fluxes_2D(rhoM,uM,vM,betaM,rhologM,betalogM,
                                                                  rhoP,uP,vP,betaP,rhologP,betalogP)

                FHx[1,i,j] = Fx1
                FHx[2,i,j] = Fx2
                FHx[3,i,j] = Fx3
                FHx[4,i,j] = Fx4
                FHy[1,i,j] = Fy1
                FHy[2,i,j] = Fy2
                FHy[3,i,j] = Fy3
                FHy[4,i,j] = Fy4
            end
        end

        for i = 1:Nfp
            iM       = Fmask[i]
            rhoM     = U[1,iM,k]
            rhouM    = U[2,iM,k]
            rhovM    = U[3,iM,k]
            EM       = U[4,iM,k]
            uM       = rhouM/rhoM
            vM       = rhovM/rhoM
            pM       = pfun(rhoM,rhouM,rhovM,EM)
            betaM    = rhoM/(2*pM)
            rhologM  = log(rhoM)
            betalogM = log(betaM)

            iP,kP = get_infoP(mapP,Fmask,i,k)
            rhoP     = U[1,iP,kP]
            rhouP    = U[2,iP,kP]
            rhovP    = U[3,iP,kP]
            EP       = U[4,iP,kP]
            uP       = rhouP/rhoP
            vP       = rhovP/rhoP
            pP       = pfun(rhoP,rhouP,rhovP,EP)
            betaP    = rhoP/(2*pP)
            rhologP  = log(rhoP)
            betalogP = log(betaP)

            Fx1,Fx2,Fx3,Fx4,Fy1,Fy2,Fy3,Fy4 = euler_fluxes_2D(rhoM,uM,vM,betaM,rhologM,betalogM,
                                                              rhoP,uP,vP,betaP,rhologP,betalogP)

            FHPx[1,i] = Fx1
            FHPx[2,i] = Fx2
            FHPx[3,i] = Fx3
            FHPx[4,i] = Fx4
            FHPy[1,i] = Fy1
            FHPy[2,i] = Fy2
            FHPy[3,i] = Fy3
            FHPy[4,i] = Fy4
        end

        Sx  =     rxJ[1,k]*Sr + sxJ[1,k]*Ss
        Sy  =     ryJ[1,k]*Sr + syJ[1,k]*Ss
        EBx = E'*(rxJ[1,k]*Br + sxJ[1,k]*Bs)
        EBy = E'*(ryJ[1,k]*Br + syJ[1,k]*Bs)

        for c = 1:Nc
            rhsvec[c,:] = 2*(Sx .* FHx[c,:,:])*ones(Np) + 2*(Sy .* FHy[c,:,:])*ones(Np) + EBx*FHPx[c,:] + EBy*FHPy[c,:]
            rhsvec[c,:] = -(1/J[1,k]*Minv)*rhsvec[c,:]
        end

        for c = 1:Nc
            rhsU[c,:,k] = rhsvec[c,:]
        end
    end
end

# Initial condition 2D shocktube
# at_left(x,y) = y-sqrt(3)*x+sqrt(3)/6 > 0.0
function vortex_sol(x,y,t)
    x0 = 5
    y0 = 0
    beta = 5
    r2 = @. (x-x0-t)^2 + (y-y0)^2

    u = @. 1 - beta*exp(1-r2)*(y-y0)/(2*pi)
    v = @. beta*exp(1-r2)*(x-x0-t)/(2*pi)
    rho = @. 1 - (1/(8*γ*pi^2))*(γ-1)/2*(beta*exp(1-r2))^2
    rho = @. rho^(1/(γ-1))
    p = @. rho^γ

    return (rho, u, v, p)
end

const Np  = size(xq,1)
const Nfp = size(Vf,1)
U = zeros(Nc,Np,K)
for k = 1:K
    for i = 1:Np
        rho,u,v,p = vortex_sol(xq[i,k],yq[i,k],0.0)
        U[1,i,k] = rho
        U[2,i,k] = rho*u
        U[3,i,k] = rho*v
        U[4,i,k] = Efun(rho,u,v,p)
    end
end

ops      = (Sr,Ss,S0r,S0s,Minv,Br_halved,Bs_halved,E)
geom     = (mapP,Fmask,xq,yq,rxJ,ryJ,sxJ,syJ,J,sJ,nx,ny,nrJ,nsJ,wf)

# Time stepping
"Time integration"
t = 0.0
U = collect(U)
rhsU = zeros(size(U))
resW = zeros(size(U))

dt_hist = []
i = 1

@time while t < T
#while i < 2
    @show t
    # SSPRK(3,3)
    dt = dt0
    rhs_ESDG!(U,rhsU,ops,geom)
    dt = min(CFL*dt,T-t)
    @. resW = U + dt*rhsU
    rhs_ESDG!(resW,rhsU,ops,geom)
    @. resW = resW+dt*rhsU
    @. resW = 3/4*U+1/4*resW
    rhs_ESDG!(resW,rhsU,ops,geom)
    @. resW = resW+dt*rhsU
    @. U = 1/3*U+2/3*resW

    push!(dt_hist,dt)
    global t = t + dt
    println("Current time $t with time step size $dt, and final time $T, at step $i")
    flush(stdout)
    global i = i + 1
end

exact_U = @. vortex_sol.(xq,yq,T)
exact_rho = [x[1] for x in exact_U]
exact_u = [x[2] for x in exact_U]
exact_v = [x[3] for x in exact_U]
exact_rhou = exact_rho .* exact_u
exact_rhov = exact_rho .* exact_v
exact_p = [x[4] for x in exact_U]
exact_E = Efun.(exact_rho,exact_u,exact_v,exact_p)

rho = U[1,:,:]
u = U[2,:,:]./U[1,:,:]
v = U[3,:,:]./U[1,:,:]
rhou = U[2,:,:]
rhov = U[3,:,:]
E = U[4,:,:]

#plotting nodes
# @unpack Vp = rd
gr(aspect_ratio=:equal,legend=false,
   markerstrokewidth=0,markersize=2)

@unpack r,s,rf,sf,wf,rq,sq,wq,nrJ,nsJ = rd
@unpack VDM,V1,Vq,Vf,Dr,Ds,M,Pq,LIFT,Vp = rd
vv = Vp*U[1,:,:]
scatter(Vp*xq,Vp*yq,vv,zcolor=vv,camera=(0,90))
savefig("~/Desktop/tmpquad.png")

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
# TODO: make it easier to combine schedules and update
# schedule mid-training

# Schedule data structures

@with_kw struct Schedules
    data::Dict = Dict()
end

@with_kw struct Schedule{T}
    nsteps::Int
    from::T
    to::T
    annealfn = anneal_linear
end

duration(schedule::Schedule) = schedule.nsteps
duration(schedules::Vector{Schedule}) = sum(duration.(schedules))
duration(schedules::Schedules) = maximum(duration.(values(schedules.data)))

function schedulevalue(schedule::Schedule, step)
    pctg = min(1, step / schedule.nsteps)
    return schedule.annealfn(pctg, schedule.from, schedule.to)
end

function schedulevalue(schedules::Vector{<:Schedule}, step)
    length(schedules) > 0 || error("Missing schedule")
    i = 1
    # FIXME: BoundsError when given too short schedule
    while duration(schedules[i]) < step
        step -= duration(schedules[i])
        i += 1
    end
    return schedulevalue(schedules[i], step)
end

# Scheduler Callback

struct ParamScheduler <: SafeCallback end

canwrite(::ParamScheduler) = (:opt,)
order(::Type{ParamScheduler}) = -80

"""
    on(::BatchBegin, phase::AbstractTrainingPhase, cb::ParamScheduler, learner)

Every batch, set all hyperparameters according to `cb.scheduledict`.
"""
function on(::BatchBegin, phase::AbstractTrainingPhase, cb::ParamScheduler, learner)
    for (P, schedule) in learner.state.schedule.data
        setoptimparam!(
            learner.opt,
            P,
            schedulevalue(schedule, learner.state.history.nsteps + 1))
    end
end

# Sample schedules

function onecycleschedule(
        nsteps,
        maxlr;
        startfactor = 1/5,
        endfactor = 1/10,
        pctstart = 0.1)

    upschedule = Schedule(
        round(Int, nsteps * pctstart),
        maxlr * startfactor,
        maxlr,
        anneal_cosine)
    downschedule = Schedule(
        round(Int, nsteps * (1-pctstart)),
        maxlr,
        maxlr * endfactor,
        anneal_cosine)

    return Dict(LR => [upschedule, downschedule])
end

# Utilities
# TODO: update

#=
function sampleschedule(
        schedule::Union{ParamSchedule, AbstractVector{ParamSchedule}};
        epochs = duration(schedule))
    xs = collect(0:0.01:epochs)
    ys = map(xs) do x
        schedulevalue(schedule, x)
    end
    xs, ys
end

function delayschedule(schedule::Dict, nepochs::Real)
    newschedule = Dict()
    for (P, schedules) in schedule
        v = schedules[1].startvalue
        delay = ParamSchedule(nepochs, v, v, anneal_const)
        newschedule[P] = vcat(delay, schedules)
    end
    return newschedule
end
=#


text:
function inv_log10(logvals)
    log10(1.0 - exp10(logvals))
end

function check_all_cols(A, B, codon_moves)
    @test A[end, end] ≈ B[1, 1]
    expected = A[end, end]
    ncols = size(A)[2]
    # if codon_moves is true, we cannot expect every column to contain
    # the correct score. every three columns should, though.
    if codon_moves
        for j in 1:ncols-2
            score = maximum(A[:, j:j+2] + B[:, j:j+2])
            @test score ≈ expected
        end
    else
        for j in 1:ncols
            score = maximum(A[:, j] + B[:, j])
            @test score ≈ expected
        end
    end
end


text:
# =============================
# Vector Field for continuous system
# =============================

"""
    vector_field(system::ContinuousIdentitySystem, x::AbstractVector)

Return the vector field state of a `ContinuousIdentitySystem`.

### Input

- `system` -- `ContinuousIdentitySystem`
- `x`      -- state (it should be any vector type)

### Output

A zeros vector of dimension `statedim`.
"""
function vector_field(system::ContinuousIdentitySystem, x::AbstractVector)
    !_is_conformable_state(system, x) && _argument_error(:x)
    return zeros(statedim(system))
end

"""
    vector_field(system::ConstrainedContinuousIdentitySystem, x::AbstractVector;
              [check_constraints]=true)

Return the vector field state of a `ConstrainedContinuousIdentitySystem`.

### Input

- `system`            -- `ConstrainedContinuousIdentitySystem`
- `x`                 -- state (it should be any vector type)
- `check_constraints` -- (optional, default: `true`) check if the state belongs to
                         the state set

### Output

A zeros vector of dimension `statedim`.
"""
function vector_field(system::ConstrainedContinuousIdentitySystem, x::AbstractVector;
                   check_constraints::Bool=true)
    !_is_conformable_state(system, x) && _argument_error(:x)
    if check_constraints
        !_in_stateset(system, x) && _argument_error(:x,:X)
    end
    return zeros(statedim(system))
end

"""
    vector_field(system::AbstractContinuousSystem, x::AbstractVector;
              [check_constraints]=true)

Return the vector field state of an `AbstractContinuousSystem`.

### Input

- `system`            -- `AbstractContinuousSystem`
- `x`                 -- state (it should be any vector type)
- `check_constraints` -- (optional, default: `true`) check if the state belongs to
                         the state set

### Output

The vector field of the system at state `x`.
"""
vector_field(system::AbstractContinuousSystem, x::AbstractVector; kwargs...) =
    _instantiate(system, x; kwargs...)

"""
    vector_field(system::AbstractContinuousSystem, x::AbstractVector, u::AbstractVector;
              [check_constraints]=true)

Return the vector field state of an `AbstractContinuousSystem`.

### Input

- `system`            -- `AbstractContinuousSystem`
- `x`                 -- state (it should be any vector type)
- `u`                 -- input (it should be any vector type) or noise, if `system`
                         is not controlled
- `check_constraints` -- (optional, default: `true`) check if the state belongs to
                         the state set

### Output

The vector field of the system at state `x` and applying input `u`.

### Notes

If the system is not controlled but noisy, the input `u` is interpreted as noise.
    """
vector_field(system::AbstractContinuousSystem, x::AbstractVector, u::AbstractVector; kwargs...) =
    _instantiate(system, x, u; kwargs...)

"""
    vector_field(system::AbstractContinuousSystem,
              x::AbstractVector, u::AbstractVector, w::AbstractVector; [check_constraints]=true)

Return the vector field state of an `AbstractContinuousSystem`.

### Input

- `system`            -- `AbstractContinuousSystem`
- `x`                 -- state (it should be any vector type)
- `u`                 -- input (it should be any vector type)
- `w`                 -- noise (it should be any vector type)
- `check_constraints` -- (optional, default: `true`) check if the state belongs to
                         the state set

### Output

The vector field of the system at state `x` and applying input `u` and noise `w`.
"""
vector_field(system::AbstractContinuousSystem, x::AbstractVector, u::AbstractVector, w::AbstractVector; kwargs...) =
    _instantiate(system, x, u, w; kwargs...)

# =============================
# VectorField type for continuous system
# =============================
"""
    VectorField{T<:Function}

Type that computes the vector field of an `AbstractContinuousSystem`.

### Fields

- `field`  -- function for calculating the vector field
"""
struct VectorField{T<:Function}
    field::T
end

# function-like evaluation
@inline function (V::VectorField)(args...)
    evaluate(V, args...)
end

function evaluate(V::VectorField, args...)
    return V.field(args...)
end

"""
    VectorField(sys::AbstractContinuousSystem)

Constructor for creating a `VectorField` from an `AbstractContinuousSystem`.

### Input

- `sys` -- `AbstractContinuousSystem`

### Ouptut

The `VectorField` for the continuous system `sys`.
"""
function VectorField(sys::AbstractContinuousSystem)
    if inputdim(sys) == 0 && noisedim(sys) == 0
        field = (x) -> vector_field(sys, x)
    elseif inputdim(sys) == 0 || noisedim(sys) == 0
        field = (x, u) -> vector_field(sys, x, u)
    else
        field = (x, u, w) -> vector_field(sys, x, u, w)
    end
    return VectorField(field)
end

function _position_value_list(V::VectorField, grid, dimx, dimy)
    X, Y = grid
    x1 = X[1]
    y1 = Y[1]
    N_POS = typeof(x1)
    N_VAL = typeof(V([x1, y1])[dimx])
    m = length(X) * length(Y)

    x = Vector{N_POS}(undef, m)  # x coordinates of each
    y = similar(x)
    vx = Vector{N_VAL}(undef, m)
    vy = similar(vx)
    k = 1
    max_entry = N_VAL(-Inf)
    for xi in X
        for yj in Y
            x[k] = xi
            y[k] = yj
            val = V([xi, yj])
            vx_k = val[dimx]
            vy_k = val[dimy]
            max_entry = max(max_entry, abs(vx_k), abs(vy_k))
            vx[k] = vx_k
            vy[k] = vy_k
            k += 1
        end
    end

    # normalize
    if max_entry > zero(N_VAL)
        step_x = (maximum(X) - minimum(X)) / length(X)
        step_y = (maximum(Y) - minimum(Y)) / length(Y)
        max_entry_x = max_entry / step_x
        max_entry_y = max_entry / step_y
        max_entry = max(max_entry_x, max_entry_y)
        for k in 1:m
            vx[k] /= max_entry
            vy[k] /= max_entry
        end
    end

    # filter out (0, 0) entries
    for k in m:-1:1
        if vx[k] == vy[k] == zero(N_VAL)
            deleteat!(x, k)
            deleteat!(y, k)
            deleteat!(vx, k)
            deleteat!(vy, k)
        end
    end

    # center arrows in the grid points
    for k in 1:length(x)
        x[k] -= vx[k] / 2
        y[k] -= vy[k] / 2
    end

    return x, y, vx, vy
end

### plot recipe
# arguments:
# - grid_points: pair containing the x and y coordinates of the grid points
#   example: the pair `([-1, 0, 1], [2, 3])` represents the grid of points
#   `[-1, 2], [0, 2], [1, 2], [-1, 3], [0, 3], [1, 3]`
# - dims: the two dimensions to plot
@recipe function plot(V::VectorField;
                      grid_points=[range(-3, stop=3, length=21),
                                   range(-3, stop=3, length=21)],
                      dims=[1, 2])
    seriestype := :quiver

    x, y, vx, vy = _position_value_list(V, grid_points, dims[1], dims[2])
    quiver := (vx, vy)

    (x, y)
end


text:
# ---
# title: 139. Word Break
# id: problem139
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/word-break/description/>
# hidden: true
# ---
#
# Given a **non-empty** string _s_ and a dictionary _wordDict_ containing a list
# of **non-empty** words, determine if _s_ can be segmented into a space-
# separated sequence of one or more dictionary words.
#
# **Note:**
#
#   * The same word in the dictionary may be reused multiple times in the segmentation.
#   * You may assume the dictionary does not contain duplicate words.
#
# **Example 1:**
#
#
#
#     Input: s = "leetcode", wordDict = ["leet", "code"]
#     Output: true
#     Explanation: Return true because "leetcode" can be segmented as "leet code".
#
#
# **Example 2:**
#
#
#
#     Input: s = "applepenapple", wordDict = ["apple", "pen"]
#     Output: true
#     Explanation: Return true because "applepenapple" can be segmented as "apple pen apple".
#                  Note that you are allowed to reuse a dictionary word.
#
#
# **Example 3:**
#
#
#
#     Input: s = "catsandog", wordDict = ["cats", "dog", "sand", "and", "cat"]
#     Output: false
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end