text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -88.40110862821636 <= q <= 80.25263273980681)

                     add_NL_constraint(m, :(log(1 + exp(-0.3402051397109518 + -0.28456924845213605*log(1 + exp(-0.42610478134649465 + 0.43223652991075534*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + -0.5575457569666615*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + -0.7868073625875285*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + 0.5078849457963623*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + -0.2883931154790984*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))) + 0.6888193490874222*log(1 + exp(-0.17528834888704425 + -0.1425247358809414*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + 0.07932201643104397*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + 0.9443705506849667*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + 0.7094815218028745*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + -0.7352206438426769*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))) + -0.026837640444652866*log(1 + exp(-0.005280183583000753 + 0.6172883711874637*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + -0.07641401576689111*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + 0.42701726221355685*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + -0.5680942985776078*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + 0.9418726475599537*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))) + -0.08739699596173622*log(1 + exp(0.26121381489212014 + -0.018200625348021138*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + -0.03849358234625422*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + -0.856604860683484*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + 0.2143578715735437*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + 0.706261759527687*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))) + 0.2360509887003217*log(1 + exp(-0.09690410604182231 + 0.6473593808554208*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + -0.6752519993458765*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + 0.771479525571356*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + 0.8906247556095672*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + 0.560146720685823*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))))) + log(1 + exp(0.2351041904461102 + 0.11979579421065711*log(1 + exp(-0.42610478134649465 + 0.43223652991075534*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + -0.5575457569666615*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + -0.7868073625875285*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + 0.5078849457963623*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + -0.2883931154790984*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))) + -0.012124046222294904*log(1 + exp(-0.17528834888704425 + -0.1425247358809414*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + 0.07932201643104397*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + 0.9443705506849667*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + 0.7094815218028745*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + -0.7352206438426769*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))) + 0.704332995743346*log(1 + exp(-0.005280183583000753 + 0.6172883711874637*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + -0.07641401576689111*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + 0.42701726221355685*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + -0.5680942985776078*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + 0.9418726475599537*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))) + 0.7162932434639226*log(1 + exp(0.26121381489212014 + -0.018200625348021138*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + -0.03849358234625422*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + -0.856604860683484*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + 0.2143578715735437*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + 0.706261759527687*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))) + -0.3351546298978829*log(1 + exp(-0.09690410604182231 + 0.6473593808554208*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + -0.6752519993458765*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + 0.771479525571356*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + 0.8906247556095672*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + 0.560146720685823*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))))) + log(1 + exp(-0.4731175976624935 + -0.5144698215629462*log(1 + exp(-0.42610478134649465 + 0.43223652991075534*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + -0.5575457569666615*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + -0.7868073625875285*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + 0.5078849457963623*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + -0.2883931154790984*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))) + -0.005430014157027596*log(1 + exp(-0.17528834888704425 + -0.1425247358809414*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + 0.07932201643104397*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + 0.9443705506849667*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + 0.7094815218028745*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + -0.7352206438426769*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))) + 0.05201762086802919*log(1 + exp(-0.005280183583000753 + 0.6172883711874637*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + -0.07641401576689111*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + 0.42701726221355685*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + -0.5680942985776078*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + 0.9418726475599537*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))) + 0.06728299071980626*log(1 + exp(0.26121381489212014 + -0.018200625348021138*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + -0.03849358234625422*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + -0.856604860683484*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + 0.2143578715735437*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + 0.706261759527687*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))) + -0.8893971601390014*log(1 + exp(-0.09690410604182231 + 0.6473593808554208*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + -0.6752519993458765*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + 0.771479525571356*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + 0.8906247556095672*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + 0.560146720685823*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))))) + log(1 + exp(-0.8380258290664089 + 0.49770714260784343*log(1 + exp(-0.42610478134649465 + 0.43223652991075534*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + -0.5575457569666615*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + -0.7868073625875285*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + 0.5078849457963623*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + -0.2883931154790984*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))) + 0.732919146629087*log(1 + exp(-0.17528834888704425 + -0.1425247358809414*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + 0.07932201643104397*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + 0.9443705506849667*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + 0.7094815218028745*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + -0.7352206438426769*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))) + -0.2983691588163282*log(1 + exp(-0.005280183583000753 + 0.6172883711874637*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + -0.07641401576689111*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + 0.42701726221355685*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + -0.5680942985776078*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + 0.9418726475599537*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))) + -0.6200389337758314*log(1 + exp(0.26121381489212014 + -0.018200625348021138*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + -0.03849358234625422*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + -0.856604860683484*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + 0.2143578715735437*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + 0.706261759527687*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))) + -0.03199113104777407*log(1 + exp(-0.09690410604182231 + 0.6473593808554208*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + -0.6752519993458765*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + 0.771479525571356*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + 0.8906247556095672*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + 0.560146720685823*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))))) + log(1 + exp(0.7168755118620993 + -0.38856981611164576*log(1 + exp(-0.42610478134649465 + 0.43223652991075534*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + -0.5575457569666615*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + -0.7868073625875285*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + 0.5078849457963623*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + -0.2883931154790984*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))) + 0.9818525806719265*log(1 + exp(-0.17528834888704425 + -0.1425247358809414*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + 0.07932201643104397*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + 0.9443705506849667*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + 0.7094815218028745*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + -0.7352206438426769*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))) + 0.9400532013818674*log(1 + exp(-0.005280183583000753 + 0.6172883711874637*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + -0.07641401576689111*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + 0.42701726221355685*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + -0.5680942985776078*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + 0.9418726475599537*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))) + -0.928831604242299*log(1 + exp(0.26121381489212014 + -0.018200625348021138*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + -0.03849358234625422*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + -0.856604860683484*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + 0.2143578715735437*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + 0.706261759527687*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))) + 0.6913287958041021*log(1 + exp(-0.09690410604182231 + 0.6473593808554208*log(1 + exp(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5]))) + -0.6752519993458765*log(1 + exp(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5]))) + 0.771479525571356*log(1 + exp(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5]))) + 0.8906247556095672*log(1 + exp(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5]))) + 0.560146720685823*log(1 + exp(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
"""
    logistic_key(x_init, r, num_keys)

Generates a vector of pseudo-random keys using the Logistic Map.

The equation -

\$x_{n+1} = r * x_{n} * (1 - x_{n})\$

# Arguments
- `x_init::Float64`: Initial value of x. x ϵ (0, 1).
- `r::Float64`: A constant value. Values > 4 usually results in pseudo-random numbers.
- `num_keys::Int64`: Number of keys to be generated.

# Returns
- `keys::Vector{Int64}:`: Generated pseudo-random keys.

# Example
```jldoctest
julia> using ChaoticEncryption

julia> logistic_key(0.01, 3.97, 20)
20-element Vector{Int64}:
   0
  44
   7
  26
  14
 224
  16
 250
 162
 211
 200
 217
  97
 132
 134
 100
 135
 232
 122
 102
```
"""
function logistic_key(x_init::Float64, r::Float64, num_keys::Int64)
    keys = Vector{Int64}()
    x = x_init

    for i = 1:num_keys
        x = r * x * (1 - x)
        key = x * 10 ^ 16 % 256
        push!(keys, trunc(Int, key))
    end

    return keys
end

ext: jl


text:
using juliawer

# push!(ARGS, "arg")
juliawer.julia_main()


text:
function mass_matrix_machine_entries!(
    mass_matrix,
    machine::M,
    global_index::Base.ImmutableDict{Symbol, Int64},
) where {M <: PSY.Machine}
    @debug "Using default mass matrix entries $M"
end

"""
Model of 0-state synchronous machine in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""
function mdl_machine_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_r::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_i::AbstractArray{<:ACCEPTED_REAL_TYPES},
    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{PSY.BaseMachine, S, A, TG, P}},
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
    Sbase = get_system_base_power(dynamic_device)
    #Obtain external states inputs for component
    external_ix = get_input_port_ix(dynamic_device, PSY.BaseMachine)
    δ = device_states[external_ix[1]]

    #Obtain inner variables for component
    V_tR = inner_vars[VR_gen_var]
    V_tI = inner_vars[VI_gen_var]

    #Get parameters
    machine = PSY.get_machine(dynamic_device)
    R = PSY.get_R(machine)
    Xd_p = PSY.get_Xd_p(machine)
    eq_p = PSY.get_eq_p(machine)
    basepower = PSY.get_base_power(dynamic_device)

    #RI to dq transformation
    V_dq = ri_dq(δ) * [V_tR; V_tI]

    #Obtain electric variables
    i_d = (1.0 / (R^2 + Xd_p^2)) * (Xd_p * (eq_p - V_dq[2]) - R * V_dq[1])  #15.36
    i_q = (1.0 / (R^2 + Xd_p^2)) * (Xd_p * V_dq[1] + R * (eq_p - V_dq[2])) #15.36
    Pe = (V_dq[1] + R * i_d) * i_d + (V_dq[2] + R * i_q) * i_q      #15.35

    #Update inner_vars
    inner_vars[τe_var] = Pe #Model assume ω approx 1.0

    #Compute current from the generator to the grid
    I_RI = (basepower / Sbase) * dq_ri(δ) * [i_d; i_q]

    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]

    return
end

"""
Model of 2-state (One d- and One q-) synchronous machine in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""
function mdl_machine_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_r::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_i::AbstractArray{<:ACCEPTED_REAL_TYPES},
    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{PSY.OneDOneQMachine, S, A, TG, P}},
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
    Sbase = get_system_base_power(dynamic_device)

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, PSY.OneDOneQMachine)

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    eq_p = internal_states[1]
    ed_p = internal_states[2]

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(dynamic_device, PSY.OneDOneQMachine)
    δ = device_states[external_ix[1]]

    #Obtain inner variables for component
    V_tR = inner_vars[VR_gen_var]
    V_tI = inner_vars[VI_gen_var]
    Vf = inner_vars[Vf_var]

    #Get parameters
    machine = PSY.get_machine(dynamic_device)
    R = PSY.get_R(machine)
    Xd = PSY.get_Xd(machine)
    Xq = PSY.get_Xq(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xq_p = PSY.get_Xq_p(machine)
    Td0_p = PSY.get_Td0_p(machine)
    Tq0_p = PSY.get_Tq0_p(machine)
    basepower = PSY.get_base_power(dynamic_device)

    #RI to dq transformation
    V_dq = ri_dq(δ) * [V_tR; V_tI]

    #Obtain electric variables
    i_d = (1.0 / (R^2 + Xd_p * Xq_p)) * (Xq_p * (eq_p - V_dq[2]) + R * (ed_p - V_dq[1]))  #15.32
    i_q = (1.0 / (R^2 + Xd_p * Xq_p)) * (-Xd_p * (ed_p - V_dq[1]) + R * (eq_p - V_dq[2]))  #15.32
    Pe = (V_dq[1] + R * i_d) * i_d + (V_dq[2] + R * i_q) * i_q      #15.35

    #Compute ODEs
    output_ode[local_ix[1]] = (1.0 / Td0_p) * (-eq_p - (Xd - Xd_p) * i_d + Vf)     #15.29 eq_p
    output_ode[local_ix[2]] = (1.0 / Tq0_p) * (-ed_p + (Xq - Xq_p) * i_q)          #15.30 ed_p

    #Update inner_vars
    inner_vars[τe_var] = Pe #Model assume ω approx 1.0

    #Compute current from the generator to the grid
    I_RI = (basepower / Sbase) * dq_ri(δ) * [i_d; i_q]

    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]

    return
end

"""
Model of 6-state (MarconatoMachine) synchronous machine in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""
function mdl_machine_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_r::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_i::AbstractArray{<:ACCEPTED_REAL_TYPES},
    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{PSY.MarconatoMachine, S, A, TG, P}},
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
    Sbase = get_system_base_power(dynamic_device)
    f0 = get_system_base_frequency(dynamic_device)

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, PSY.MarconatoMachine)

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    ψq = internal_states[1]
    ψd = internal_states[2]
    eq_p = internal_states[3]
    ed_p = internal_states[4]
    eq_pp = internal_states[5]
    ed_pp = internal_states[6]

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(dynamic_device, PSY.MarconatoMachine)
    δ = device_states[external_ix[1]]
    ω = device_states[external_ix[2]]

    #Obtain inner variables for component
    V_tR = inner_vars[VR_gen_var]
    V_tI = inner_vars[VI_gen_var]
    Vf = inner_vars[Vf_var]

    #Get parameters
    machine = PSY.get_machine(dynamic_device)
    R = PSY.get_R(machine)
    Xd = PSY.get_Xd(machine)
    Xq = PSY.get_Xq(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xq_p = PSY.get_Xq_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = PSY.get_Xq_pp(machine)
    Td0_p = PSY.get_Td0_p(machine)
    Tq0_p = PSY.get_Tq0_p(machine)
    Td0_pp = PSY.get_Td0_pp(machine)
    Tq0_pp = PSY.get_Tq0_pp(machine)
    T_AA = PSY.get_T_AA(machine)
    γd = PSY.get_γd(machine)
    γq = PSY.get_γq(machine)
    basepower = PSY.get_base_power(dynamic_device)

    #RI to dq transformation
    V_dq = ri_dq(δ) * [V_tR; V_tI]

    #Obtain electric variables
    i_d = (1.0 / Xd_pp) * (eq_pp - ψd)      #15.18
    i_q = (1.0 / Xq_pp) * (-ed_pp - ψq)     #15.18
    τ_e = ψd * i_q - ψq * i_d               #15.6

    #Compute ODEs
    output_ode[local_ix[1]] = 2 * π * f0 * (R * i_q - ω * ψd + V_dq[2])                        #15.9 ψq
    output_ode[local_ix[2]] = 2 * π * f0 * (R * i_d + ω * ψq + V_dq[1])                        #15.9 ψd
    output_ode[local_ix[3]] =
        ((1.0 / Td0_p) * (-eq_p - (Xd - Xd_p - γd) * i_d + (1 - (T_AA / Td0_p)) * Vf))                             #15.16 eq_p
    output_ode[local_ix[4]] = (1.0 / Tq0_p) * (-ed_p + (Xq - Xq_p - γq) * i_q)             #15.16 ed_p
    output_ode[local_ix[5]] =
        ((1.0 / Td0_pp) * (-eq_pp + eq_p - (Xd_p - Xd_pp + γd) * i_d + (T_AA / Td0_p) * Vf))       #15.16 eq_pp
    output_ode[local_ix[6]] = (1.0 / Tq0_pp) * (-ed_pp + ed_p + (Xq_p - Xq_pp + γq) * i_q) #15.16 ed_pp

    #Update inner_vars
    inner_vars[τe_var] = τ_e

    #Compute current from the generator to the grid
    I_RI = (basepower / Sbase) * dq_ri(δ) * [i_d; i_q]

    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]

    return
end

"""
Model of 4-state (SimpleMarconatoMachine) synchronous machine in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""
function mdl_machine_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_r::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_i::AbstractArray{<:ACCEPTED_REAL_TYPES},
    dynamic_device::DynamicWrapper{
        PSY.DynamicGenerator{PSY.SimpleMarconatoMachine, S, A, TG, P},
    },
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
    Sbase = get_system_base_power(dynamic_device)
    f0 = get_system_base_frequency(dynamic_device)
    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, PSY.SimpleMarconatoMachine)

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    eq_p = internal_states[1]
    ed_p = internal_states[2]
    eq_pp = internal_states[3]
    ed_pp = internal_states[4]

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(dynamic_device, PSY.SimpleMarconatoMachine)
    δ = device_states[external_ix[1]]

    #Obtain inner variables for component
    V_tR = inner_vars[VR_gen_var]
    V_tI = inner_vars[VI_gen_var]
    Vf = inner_vars[Vf_var]

    #Get parameters
    machine = PSY.get_machine(dynamic_device)
    R = PSY.get_R(machine)
    Xd = PSY.get_Xd(machine)
    Xq = PSY.get_Xq(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xq_p = PSY.get_Xq_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = PSY.get_Xq_pp(machine)
    Td0_p = PSY.get_Td0_p(machine)
    Tq0_p = PSY.get_Tq0_p(machine)
    Td0_pp = PSY.get_Td0_pp(machine)
    Tq0_pp = PSY.get_Tq0_pp(machine)
    T_AA = PSY.get_T_AA(machine)
    γd = PSY.get_γd(machine)
    γq = PSY.get_γq(machine)
    basepower = PSY.get_base_power(dynamic_device)

    #RI to dq transformation
    V_dq = ri_dq(δ) * [V_tR; V_tI]

    #Obtain electric variables
    #i_dq = inv([[R -Xq_pp]; [Xq_pp R]])*[ed_pp - V_dq[1]; eq_pp - V_dq[2]]
    i_d =
        (1.0 / (R^2 + Xd_pp * Xq_pp)) * (Xq_pp * (eq_pp - V_dq[2]) + R * (ed_pp - V_dq[1]))      #15.25
    i_q =
        (1.0 / (R^2 + Xd_pp * Xq_pp)) * (-Xd_pp * (ed_pp - V_dq[1]) + R * (eq_pp - V_dq[2]))      #15.25
    τ_e = (V_dq[1] + R * i_d) * i_d + (V_dq[2] + R * i_q) * i_q               #15.6

    #Compute ODEs
    output_ode[local_ix[1]] =
        ((1.0 / Td0_p) * (-eq_p - (Xd - Xd_p - γd) * i_d + (1 - (T_AA / Td0_p)) * Vf))                             #15.16 eq_p
    output_ode[local_ix[2]] = (1.0 / Tq0_p) * (-ed_p + (Xq - Xq_p - γq) * i_q)             #15.16 ed_p
    output_ode[local_ix[3]] =
        ((1.0 / Td0_pp) * (-eq_pp + eq_p - (Xd_p - Xd_pp + γd) * i_d + (T_AA / Td0_p) * Vf))       #15.16 eq_pp
    output_ode[local_ix[4]] = (1.0 / Tq0_pp) * (-ed_pp + ed_p + (Xq_p - Xq_pp + γq) * i_q) #15.16 ed_pp

    #Update inner_vars
    inner_vars[τe_var] = τ_e

    #Compute current from the generator to the grid
    I_RI = (basepower / Sbase) * dq_ri(δ) * [i_d; i_q]

    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]

    return
end

"""
Model of 6-state (AndersonFouadMachine) synchronous machine in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""
function mdl_machine_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_r::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_i::AbstractArray{<:ACCEPTED_REAL_TYPES},
    dynamic_device::DynamicWrapper{
        PSY.DynamicGenerator{PSY.AndersonFouadMachine, S, A, TG, P},
    },
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
    Sbase = get_system_base_power(dynamic_device)
    f0 = get_system_base_frequency(dynamic_device)

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, PSY.AndersonFouadMachine)

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    ψq = internal_states[1]
    ψd = internal_states[2]
    eq_p = internal_states[3]
    ed_p = internal_states[4]
    eq_pp = internal_states[5]
    ed_pp = internal_states[6]

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(dynamic_device, PSY.AndersonFouadMachine)
    δ = device_states[external_ix[1]]
    ω = device_states[external_ix[2]]

    #Obtain inner variables for component
    V_tR = inner_vars[VR_gen_var]
    V_tI = inner_vars[VI_gen_var]
    Vf = inner_vars[Vf_var]

    #Get parameters
    machine = PSY.get_machine(dynamic_device)
    R = PSY.get_R(machine)
    Xd = PSY.get_Xd(machine)
    Xq = PSY.get_Xq(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xq_p = PSY.get_Xq_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = PSY.get_Xq_pp(machine)
    Td0_p = PSY.get_Td0_p(machine)
    Tq0_p = PSY.get_Tq0_p(machine)
    Td0_pp = PSY.get_Td0_pp(machine)
    Tq0_pp = PSY.get_Tq0_pp(machine)
    basepower = PSY.get_base_power(dynamic_device)

    #RI to dq transformation
    V_dq = ri_dq(δ) * [V_tR; V_tI]

    #Obtain electric variables
    i_d = (1.0 / Xd_pp) * (eq_pp - ψd)      #15.18
    i_q = (1.0 / Xq_pp) * (-ed_pp - ψq)     #15.18
    τ_e = ψd * i_q - ψq * i_d               #15.6

    #Compute ODEs
    output_ode[local_ix[1]] = 2 * π * f0 * (R * i_q - ω * ψd + V_dq[2])                        #15.9 ψq
    output_ode[local_ix[2]] = 2 * π * f0 * (R * i_d + ω * ψq + V_dq[1])                        #15.9 ψd
    output_ode[local_ix[3]] = (1.0 / Td0_p) * (-eq_p + (Xd - Xd_p) * i_d + Vf)             #15.19 eq_p
    output_ode[local_ix[4]] = (1.0 / Tq0_p) * (-ed_p + (Xq - Xq_p) * i_q)                  #15.19 ed_p
    output_ode[local_ix[5]] = (1.0 / Td0_pp) * (-eq_pp + eq_p - (Xd_p - Xd_pp) * i_d)      #15.19 eq_pp
    output_ode[local_ix[6]] = (1.0 / Tq0_pp) * (-ed_pp + ed_p + (Xq_p - Xq_pp) * i_q)      #15.19 ed_pp

    #Update inner_vars
    inner_vars[τe_var] = τ_e

    #Compute current from the generator to the grid
    I_RI = (basepower / Sbase) * dq_ri(δ) * [i_d; i_q]

    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]

    return
end

"""
Model of 4-state (SimpleAFMachine) synchronous machine in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""
function mdl_machine_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_r::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_i::AbstractArray{<:ACCEPTED_REAL_TYPES},
    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{PSY.SimpleAFMachine, S, A, TG, P}},
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
    Sbase = get_system_base_power(dynamic_device)
    f0 = get_system_base_frequency(dynamic_device)

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, PSY.SimpleAFMachine)

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    eq_p = internal_states[1]
    ed_p = internal_states[2]
    eq_pp = internal_states[3]
    ed_pp = internal_states[4]

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(dynamic_device, PSY.SimpleAFMachine)
    δ = device_states[external_ix[1]]
    ω = device_states[external_ix[2]]

    #Obtain inner variables for component
    V_tR = inner_vars[VR_gen_var]
    V_tI = inner_vars[VI_gen_var]
    Vf = inner_vars[Vf_var]

    #Get parameters
    machine = PSY.get_machine(dynamic_device)
    R = PSY.get_R(machine)
    Xd = PSY.get_Xd(machine)
    Xq = PSY.get_Xq(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xq_p = PSY.get_Xq_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = PSY.get_Xq_pp(machine)
    Td0_p = PSY.get_Td0_p(machine)
    Tq0_p = PSY.get_Tq0_p(machine)
    Td0_pp = PSY.get_Td0_pp(machine)
    Tq0_pp = PSY.get_Tq0_pp(machine)
    basepower = PSY.get_base_power(dynamic_device)

    #RI to dq transformation
    V_dq = ri_dq(δ) * [V_tR; V_tI]

    #Obtain electric variables
    i_d =
        (1.0 / (R^2 + Xd_pp * Xq_pp)) * (Xq_pp * (eq_pp - V_dq[2]) + R * (ed_pp - V_dq[1]))      #15.25
    i_q =
        (1.0 / (R^2 + Xd_pp * Xq_pp)) * (-Xd_pp * (ed_pp - V_dq[1]) + R * (eq_pp - V_dq[2]))      #15.25
    τ_e = (V_dq[1] + R * i_d) * i_d + (V_dq[2] + R * i_q) * i_q               #15.6

    #Compute ODEs
    output_ode[local_ix[1]] = (1.0 / Td0_p) * (-eq_p + (Xd - Xd_p) * i_d + Vf)            #15.19 eq_p
    output_ode[local_ix[2]] = (1.0 / Tq0_p) * (-ed_p + (Xq - Xq_p) * i_q)                 #15.19 ed_p
    output_ode[local_ix[3]] = (1.0 / Td0_pp) * (-eq_pp + eq_p - (Xd_p - Xd_pp) * i_d)     #15.19 eq_pp
    output_ode[local_ix[4]] = (1.0 / Tq0_pp) * (-ed_pp + ed_p + (Xq_p - Xq_pp) * i_q)     #15.19 ed_pp

    #Update inner_vars
    inner_vars[τe_var] = τ_e

    #Compute current from the generator to the grid
    I_RI = (basepower / Sbase) * dq_ri(δ) * [i_d; i_q]

    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]

    return
end

function mdl_machine_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_r::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_i::AbstractArray{<:ACCEPTED_REAL_TYPES},
    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, A, TG, P}},
) where {
    M <: Union{PSY.RoundRotorQuadratic, PSY.RoundRotorExponential},
    S <: PSY.Shaft,
    A <: PSY.AVR,
    TG <: PSY.TurbineGov,
    P <: PSY.PSS,
}
    Sbase = get_system_base_power(dynamic_device)
    f0 = get_system_base_frequency(dynamic_device)

    #Obtain indices for component w/r to device
    machine = PSY.get_machine(dynamic_device)
    local_ix = get_local_state_ix(dynamic_device, M)

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    eq_p = internal_states[1]
    ed_p = internal_states[2]
    ψ_kd = internal_states[3]
    ψ_kq = internal_states[4]

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(dynamic_device, M)
    δ = device_states[external_ix[1]]

    #Obtain inner variables for component
    V_tR = inner_vars[VR_gen_var]
    V_tI = inner_vars[VI_gen_var]
    Vf = inner_vars[Vf_var] #E_fd: Field voltage

    #Get parameters
    R = PSY.get_R(machine)
    Td0_p = PSY.get_Td0_p(machine)
    Td0_pp = PSY.get_Td0_pp(machine)
    Tq0_p = PSY.get_Tq0_p(machine)
    Tq0_pp = PSY.get_Tq0_pp(machine)
    Xd = PSY.get_Xd(machine)
    Xq = PSY.get_Xq(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xq_p = PSY.get_Xq_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = Xd_pp
    Xl = PSY.get_Xl(machine)
    γ_d1 = PSY.get_γ_d1(machine)
    γ_q1 = PSY.get_γ_q1(machine)
    γ_d2 = PSY.get_γ_d2(machine)
    γ_q2 = PSY.get_γ_q2(machine)
    γ_qd = PSY.get_γ_qd(machine)
    basepower = PSY.get_base_power(dynamic_device)

    #RI to dq transformation
    V_dq = ri_dq(δ) * [V_tR; V_tI]

    #Additional Fluxes
    ψq_pp = γ_q1 * ed_p + ψ_kq * (1 - γ_q1)
    ψd_pp = γ_d1 * eq_p + γ_d2 * (Xd_p - Xl) * ψ_kd
    ψ_pp = sqrt(ψd_pp^2 + ψq_pp^2)
    #Currents
    I_d =
        (1.0 / (R^2 + Xq_pp * Xd_pp)) *
        (-R * (V_dq[1] - ψq_pp) + Xq_pp * (-V_dq[2] + ψd_pp))
    I_q =
        (1.0 / (R^2 + Xq_pp * Xd_pp)) * (Xd_pp * (V_dq[1] - ψq_pp) + R * (-V_dq[2] + ψd_pp))
    Se = saturation_function(machine, ψ_pp)
    Xad_Ifd = eq_p + (Xd - Xd_p) * (γ_d1 * I_d - γ_d2 * ψ_kd + γ_d2 * eq_p) + Se * ψd_pp
    Xaq_I1q =
        ed_p + (Xq - Xq_p) * (γ_q2 * ed_p - γ_q2 * ψ_kq - γ_q1 * I_q) + Se * ψq_pp * γ_qd
    τ_e = I_d * (V_dq[1] + I_d * R) + I_q * (V_dq[2] + I_q * R)

    #Compute ODEs
    output_ode[local_ix[1]] = (1.0 / Td0_p) * (Vf - Xad_Ifd)                        #2.20b eq_p
    output_ode[local_ix[2]] = (1.0 / Tq0_p) * (-Xaq_I1q)                            #2.21b ed_p
    output_ode[local_ix[3]] = (1.0 / Td0_pp) * (-ψ_kd + eq_p - (Xd_p - Xl) * I_d)   #2.20c ψ_kd
    output_ode[local_ix[4]] = (1.0 / Tq0_pp) * (-ψ_kq + ed_p + (Xq_p - Xl) * I_q)   #2.21c ψ_kq

    #Update inner_vars
    inner_vars[τe_var] = τ_e
    inner_vars[Xad_Ifd_var] = Xad_Ifd

    #Compute current from the generator to the grid
    I_RI = (basepower / Sbase) * dq_ri(δ) * [I_d; I_q]

    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]

    return
end

function mdl_machine_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_r::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_i::AbstractArray{<:ACCEPTED_REAL_TYPES},
    dynamic_device::DynamicWrapper{
        PSY.DynamicGenerator{PSY.SalientPoleQuadratic, S, A, TG, P},
    },
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
    Sbase = get_system_base_power(dynamic_device)
    f0 = get_system_base_frequency(dynamic_device)

    #Obtain indices for component w/r to device
    machine = PSY.get_machine(dynamic_device)
    local_ix = get_local_state_ix(dynamic_device, typeof(machine))

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    eq_p = internal_states[1]
    ψ_kd = internal_states[2]
    ψq_pp = internal_states[3]

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(dynamic_device, typeof(machine))
    δ = device_states[external_ix[1]]

    #Obtain inner variables for component
    V_tR = inner_vars[VR_gen_var]
    V_tI = inner_vars[VI_gen_var]
    Vf = inner_vars[Vf_var] #E_fd: Field voltage

    #Get parameters
    R = PSY.get_R(machine)
    Td0_p = PSY.get_Td0_p(machine)
    Td0_pp = PSY.get_Td0_pp(machine)
    Tq0_pp = PSY.get_Tq0_pp(machine)
    Xd = PSY.get_Xd(machine)
    Xq = PSY.get_Xq(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = Xd_pp
    Xl = PSY.get_Xl(machine)
    γ_d1 = PSY.get_γ_d1(machine)
    γ_q1 = PSY.get_γ_q1(machine)
    γ_d2 = PSY.get_γ_d2(machine)
    basepower = PSY.get_base_power(dynamic_device)

    #RI to dq transformation
    V_d, V_q = ri_dq(δ) * [V_tR; V_tI]

    #Additional Fluxes
    ψd_pp = γ_d1 * eq_p + γ_q1 * ψ_kd

    #Currents
    I_d = (1.0 / (R^2 + Xd_pp^2)) * (-R * (V_d + ψq_pp) + Xd_pp * (ψd_pp - V_q))
    I_q = (1.0 / (R^2 + Xd_pp^2)) * (Xd_pp * (V_d + ψq_pp) + R * (ψd_pp - V_q))
    Se = saturation_function(machine, eq_p)
    Xad_Ifd =
        eq_p + Se * eq_p + (Xd - Xd_p) * (I_d + γ_d2 * (eq_p - ψ_kd - (Xd_p - Xl) * I_d))
    τ_e = I_d * (V_d + I_d * R) + I_q * (V_q + I_q * R)

    #Compute ODEs
    output_ode[local_ix[1]] = (1.0 / Td0_p) * (Vf - Xad_Ifd)                        #2.34b eq_p
    output_ode[local_ix[2]] = (1.0 / Td0_pp) * (-ψ_kd + eq_p - (Xd_p - Xl) * I_d)   #2.34c ψ_kd
    output_ode[local_ix[3]] = (1.0 / Tq0_pp) * (-ψq_pp - (Xq - Xq_pp) * I_q)        #2.35b ψq_pp

    #Update inner_vars
    inner_vars[τe_var] = τ_e
    inner_vars[Xad_Ifd_var] = Xad_Ifd

    #Compute current from the generator to the grid
    I_RI = (basepower / Sbase) * dq_ri(δ) * [I_d; I_q]

    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]

    return
end

function mdl_machine_ode!(
    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_r::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_i::AbstractArray{<:ACCEPTED_REAL_TYPES},
    dynamic_device::DynamicWrapper{
        PSY.DynamicGenerator{PSY.SalientPoleExponential, S, A, TG, P},
    },
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
    Sbase = get_system_base_power(dynamic_device)
    f0 = get_system_base_frequency(dynamic_device)

    #Obtain indices for component w/r to device
    machine = PSY.get_machine(dynamic_device)
    local_ix = get_local_state_ix(dynamic_device, typeof(machine))

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    eq_p = internal_states[1]
    ψ_kd = internal_states[2]
    ψq_pp = internal_states[3]

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(dynamic_device, typeof(machine))
    δ = device_states[external_ix[1]]

    #Obtain inner variables for component
    V_tR = inner_vars[VR_gen_var]
    V_tI = inner_vars[VI_gen_var]
    Vf = inner_vars[Vf_var] #E_fd: Field voltage

    #Get parameters
    R = PSY.get_R(machine)
    Td0_p = PSY.get_Td0_p(machine)
    Td0_pp = PSY.get_Td0_pp(machine)
    Tq0_pp = PSY.get_Tq0_pp(machine)
    Xd = PSY.get_Xd(machine)
    Xq = PSY.get_Xq(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = Xd_pp
    Xl = PSY.get_Xl(machine)
    γ_d1 = PSY.get_γ_d1(machine)
    γ_q1 = PSY.get_γ_q1(machine)
    γ_d2 = PSY.get_γ_d2(machine)
    γ_qd = (Xq - Xl) / (Xd - Xl)
    basepower = PSY.get_base_power(dynamic_device)

    #RI to dq transformation
    V_d, V_q = ri_dq(δ) * [V_tR; V_tI]

    #Additional Fluxes
    ψd_pp = γ_d1 * eq_p + γ_q1 * ψ_kd
    ψ_pp = sqrt(ψd_pp^2 + ψq_pp^2)

    #Currents
    I_d = (1.0 / (R^2 + Xd_pp^2)) * (-R * (V_d - ψq_pp) + Xq_pp * (-V_q + ψd_pp))
    I_q = (1.0 / (R^2 + Xd_pp^2)) * (Xd_pp * (V_d - ψq_pp) + R * (-V_q + ψd_pp))
    Se = saturation_function(machine, ψ_pp)
    Xad_Ifd =
        eq_p + Se * ψd_pp + (Xd - Xd_p) * (I_d + γ_d2 * (eq_p - ψ_kd - (Xd_p - Xl) * I_d))
    τ_e = I_d * (V_d + I_d * R) + I_q * (V_q + I_q * R)

    #Compute ODEs
    output_ode[local_ix[1]] = (1.0 / Td0_p) * (Vf - Xad_Ifd)                        #2.34b eq_p
    output_ode[local_ix[2]] = (1.0 / Td0_pp) * (-ψ_kd + eq_p - (Xd_p - Xl) * I_d)   #2.34c ψ_kd
    output_ode[local_ix[3]] =
        (1.0 / Tq0_pp) * (-ψq_pp + (Xq - Xq_pp) * I_q - Se * γ_qd * ψq_pp)        #2.35b ψq_pp

    #Update inner_vars
    inner_vars[τe_var] = τ_e
    inner_vars[Xad_Ifd_var] = Xad_Ifd

    #Compute current from the generator to the grid
    I_RI = (basepower / Sbase) * dq_ri(δ) * [I_d; I_q]

    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]

    return
end

#Not already implemented:
#=
"""
Model of 5-state (FullOrderMachine) synchronous machine in Julia.
Refer to Power System Dynamics: Stability and Control, by J. Machowski, J. Bialek and J. Bumby,
or Power System Stability and Control by P. Kundur, for the equations
"""

function mdl_machine_ode!(
    device_states,
    output_ode,
inner_vars,
current_r::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_i::AbstractArray{<:ACCEPTED_REAL_TYPES},
    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{PSY.FullMachine, S, A, TG, P}},
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, PSY.FullMachine)

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    ψd = internal_states[1]
    ψq = internal_states[2]
    ψf = internal_states[3]
    ψ1d = internal_states[4]
    ψ1q = internal_states[5]

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(dynamic_device, PSY.FullMachine)
    δ = device_states[external_ix[1]]
    ω = device_states[external_ix[2]]

    #Obtain inner variables for component
    V_tR = inner_vars[VR_gen_var]
    V_tI = inner_vars[VI_gen_var]
    Vf = inner_vars[Vf_var]
    #ψd = inner_vars[ψd_var]
    #ψq = inner_vars[ψq_var]

    #Get parameters
    machine = PSY.get_machine(dynamic_device)
    R = PSY.get_R(machine)
    R_f = PSY.get_R_f(machine)
    R_1d = PSY.get_R_1d(machine)
    R_1q = PSY.get_R_1q(machine)
    inv_d_fluxlink = PSY.get_inv_d_fluxlink(machine)
    inv_q_fluxlink = PSY.get_inv_q_fluxlink(machine)
    basepower = PSY.get_base_power(dynamic_device)

    #RI to dq transformation
    V_dq = ri_dq(δ) * [V_tR; V_tI]

    #Obtain electric currents] using Flux Linkage equations
    i_df1d = inv_d_fluxlink * [ψd; ψf; ψ1d]  #11.18 in Machowski (3.127, 3.130 and 3.131 in Kundur)
    i_q1q = inv_q_fluxlink * [ψq; ψ1q]        #11.19 in Machowski (3.128 and 3.132 in Kundur)

    #Compute ODEs
    output_ode[local_ix[1]] = 2 * pi * f0 * (ψq * ω + R * i_df1d[1] + V_dq[1])  #11.32 in Machowski or 3.120 in Kundur
    output_ode[local_ix[2]] = 2 * pi * f0 * (-ψd * ω + R * i_q1q[1] + V_dq[2])
    output_ode[local_ix[3]] = -R_f * i_df1d[2] + Vf      #11.33 in Machowski or 3.123 in Kundur
    output_ode[local_ix[4]] = -R_1d * i_df1d[3]           #11.33 in Machowski or 3.124 in Kundur
    output_ode[local_ix[5]] = -R_1q * i_q1q[2]            #11.33 in Machowski or 3.125 in Kundur

    #Update stator fluxes and torque
    #ψd = (1.0/ω)*(R*i_q1q[1] + V_dq[2])         #11.32 in Machowski or 3.121 in Kundur
    #ψq = (1.0/ω)*(-R*i_df1d[1] - V_dq[1])       #11.32 in Machowski or 3.120 in Kundur
    τ_e = ψd * i_q1q[1] - ψq * i_df1d[1]            #15.6 in Milano or 3.117 in Kundur

    #Update inner_vars
    inner_vars[τe_var] = τ_e
    inner_vars[ψd_var] = ψd
    inner_vars[ψq_var] = ψq

    #Compute current from the generator to the grid
    I_RI = basepower * dq_ri(δ) * [i_df1d[1]; i_q1q[1]]

    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]

    return
end

"""
Model of 3-state (FullOrderMachine) synchronous machine in Julia.
Refer to Power System Dynamics: Stability and Control, by J. Machowski, J. Bialek and J. Bumby,
or Power System Stability and Control by P. Kundur, for the equations
"""

function mdl_machine_ode!(
    device_states,
    output_ode,
inner_vars,
current_r::AbstractArray{<:ACCEPTED_REAL_TYPES},
    current_i::AbstractArray{<:ACCEPTED_REAL_TYPES},
    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{PSY.SimpleFullMachine, S, A, TG, P}},
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, PSY.SimpleFullMachine)

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    ψf = internal_states[1]
    ψ1d = internal_states[2]
    ψ1q = internal_states[3]

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(dynamic_device, PSY.SimpleFullMachine)
    δ = device_states[external_ix[1]]
    ω = device_states[external_ix[2]]

    #Obtain inner variables for component
    V_tR = inner_vars[VR_gen_var]
    V_tI = inner_vars[VI_gen_var]
    Vf = inner_vars[Vf_var]
    ψd = inner_vars[ψd_var]
    ψq = inner_vars[ψq_var]

    #Get parameters
    machine = PSY.get_machine(dynamic_device)
    R = PSY.get_R(machine)
    R_f = PSY.get_R_f(machine)
    R_1d = PSY.get_R_1d(machine)
    R_1q = PSY.get_R_1q(machine)
    inv_d_fluxlink = PSY.get_inv_d_fluxlink(machine)
    inv_q_fluxlink = PSY.get_inv_q_fluxlink(machine)
    BaseMVA = PSY.get_MVABase(machine)

    #RI to dq transformation
    V_dq = ri_dq(δ) * [V_tR; V_tI]

    #Obtain electric currents] using Flux Linkage equations
    i_df1d = inv_d_fluxlink * [ψd; ψf; ψ1d]  #11.18 in Machowski (3.127, 3.130 and 3.131 in Kundur)
    i_q1q = inv_q_fluxlink * [ψq; ψ1q]        #11.19 in Machowski (3.128 and 3.132 in Kundur)

    #Compute ODEs
    output_ode[local_ix[1]] = -R_f * i_df1d[2] + Vf      #11.33 in Machowski or 3.123 in Kundur
    output_ode[local_ix[2]] = -R_1d * i_df1d[3]           #11.33 in Machowski or 3.124 in Kundur
    output_ode[local_ix[3]] = -R_1q * i_q1q[2]            #11.33 in Machowski or 3.125 in Kundur

    #Update stator fluxes and torque
    ψd = (1.0 / ω) * (R * i_q1q[1] + V_dq[2])         #11.32 in Machowski or 3.121 in Kundur
    ψq = (1.0 / ω) * (-R * i_df1d[1] - V_dq[1])       #11.32 in Machowski or 3.120 in Kundur
    τ_e = ψd * i_q1q[1] - ψq * i_df1d[1]            #15.6 in Milano or 3.117 in Kundur

    #Update inner_vars
    inner_vars[τe_var] = τ_e
    inner_vars[ψd_var] = ψd
    inner_vars[ψq_var] = ψq

    #Compute current from the generator to the grid
    I_RI = (BaseMVA / Sbase) * dq_ri(δ) * [i_df1d[1]; i_q1q[1]]

    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]

    return
end
=#


text:
__precompile__(true)
module GaussianFilters

using LinearAlgebra
using ForwardDiff
using Random
import Random: rand

# Kalman, Extended Kalman, Unscented Kalman Filters

export
	DynamicsModel,
	LinearDynamicsModel,
	NonlinearDynamicsModel,
	ObservationModel,
	LinearObservationModel,
	NonlinearObservationModel,
	predict,
	measure
include("models.jl")

export
	AbstractFilter,
	KalmanFilter,
	ExtendedKalmanFilter,
	UnscentedKalmanFilter,
	GaussianBelief

include("kf_classes.jl")

export
	simulate_step,
	simulation,
	run_filter,
	unpack

include("simulate.jl")

export	update
include("kf.jl")
include("ekf.jl")

export
	unscented_transform,
	unscented_transform_inverse
include("ukf.jl")

# Gaussian Mixture PHD Filter

export
	Measurement,
	Dynamics,
	GaussianMixture,
	Spawn,
	PHDFilter
include("gmphd_classes.jl")

export
	update,
	predict,
	measure,
	prune
include("gmphd.jl")

export
	multiple_target_state_extraction
include("gmphd_extraction.jl")

# TODO: Add Labeled Multi-Bernoulli Filter

# Utilities
export
	belief_ellipse
include("utils.jl")

end # module


text:
"""
Abstract type for time series storage implementations.

All subtypes must implement:
- serialize_time_series!
- add_time_series_reference!
- remove_time_series!
- deserialize_time_series
- clear_time_series!
- get_num_time_series
- check_read_only
- is_read_only
"""
abstract type TimeSeriesStorage end

function make_time_series_storage(;
    in_memory = false,
    filename = nothing,
    directory = nothing,
)
    if in_memory
        storage = InMemoryTimeSeriesStorage()
    elseif !isnothing(filename)
        storage = Hdf5TimeSeriesStorage(; filename = filename)
    else
        storage = Hdf5TimeSeriesStorage(true; directory = directory)
    end

    return storage
end

const COMPONENT_name_DELIMITER = "__"

function make_component_name(component_uuid::UUIDs.UUID, name::AbstractString)
    return string(component_uuid) * COMPONENT_name_DELIMITER * name
end

function deserialize_component_name(component_name::AbstractString)
    data = split(component_name, COMPONENT_name_DELIMITER)
    component = UUIDs.UUID(data[1])
    name = data[2]
    return component, name
end

function serialize(storage::TimeSeriesStorage, file_path::AbstractString)
    if storage isa Hdf5TimeSeriesStorage
        if abspath(get_file_path(storage)) == abspath(file_path)
            if !is_read_only(storage)
                error("Attempting to overwrite identical time series file")
            end

            @debug "Skip time series serialization because the paths are identical"
            return
        end

        # The data is currently in a temp file, so we can just make a copy.
        copy_file(get_file_path(storage), file_path)
    elseif storage isa InMemoryTimeSeriesStorage
        convert_to_hdf5(storage, file_path)
    else
        error("unsupported type $(typeof(storage))")
    end

    @info "Serialized time series data to $file_path."
end


text:
using Test, Lattices

lattice = Chain(4)
nameof(lattice)
@test nameof(lattice) == "Chain Lattice"
@test collect(sites(lattice)) == [1, 2, 3, 4]
@test collect(edges(lattice)) == [(1, 2), (2, 3), (3, 4), (4, 1)]

@test lattice[1] == 1
@test lattice[2] == 2
@test lattice[3] == 3
@test lattice[4] == 4
@test lattice[5] == 1

lattice = Chain(4, boundary=Fixed)
@test_throws BoundsError lattice[5]


text:
"""
**Bayesian SVM**

The [Bayesian SVM](https://arxiv.org/abs/1707.05532) is a Bayesian interpretation of the classical SVM.
``p(y|f) \\propto \\exp\\left(2\\max(1-yf,0)\\right)``

```julia
BayesianSVM()
```
---
For the analytic version of the likelihood, it is augmented via:
```math
p(y|f,\\omega) = \\frac{1}{\\sqrt{2\\pi\\omega}}\\exp\\left(-\\frac{1}{2}\\frac{(1+\\omega-yf)^2}{\\omega}\\right)
```
where ``\\omega\\sim 1_{[0,\\infty]}`` has an improper prior (his posterior is however has a valid distribution (Generalized Inverse Gaussian)). For reference [see this paper](http://ecmlpkdd2017.ijs.si/papers/paperID502.pdf)
"""
struct BayesianSVM{T<:Real} <: ClassificationLikelihood{T}
    ω::AbstractVector{AbstractVector{T}}
    θ::AbstractVector{AbstractVector{T}}
    function BayesianSVM{T}() where {T<:Real}
        new{T}()
    end
    function BayesianSVM{T}(ω::AbstractVector{<:AbstractVector{<:Real}},θ::AbstractVector{<:AbstractVector{<:Real}}) where {T<:Real}
        new{T}(ω,θ)
    end
end

function BayesianSVM()
    BayesianSVM{Float64}()
end

function init_likelihood(likelihood::BayesianSVM{T},inference::Inference{T},nLatent::Integer,nSamplesUsed) where T
    BayesianSVM{T}([abs.(rand(T,nSamplesUsed)) for _ in 1:nLatent],[zeros(T,nSamplesUsed) for _ in 1:nLatent])
end
function pdf(l::BayesianSVM,y::Real,f::Real)
    svmlikelihood(y*f)
end

function Base.show(io::IO,model::BayesianSVM{T}) where T
    print(io,"Bayesian SVM")
end

"""Return likelihood equivalent to SVM hinge loss"""
function svmlikelihood(f::Real)
    pos = svmpseudolikelihood(f)
    return pos./(pos.+svmpseudolikelihood(-f))
end

"""Return the pseudo likelihood of the SVM hinge loss"""
function svmpseudolikelihood(f::Real)
    return exp(-2.0*max.(1.0-f,0))
end

function compute_proba(l::BayesianSVM{T},μ::Vector{T},σ²::Vector{T}) where {T<:Real}
    N = length(μ)
    pred = zeros(T,N)
    for i in 1:N
        if σ²[i] <= 0.0
            pred[i] = svmlikelihood(μ[i])
        else
            pred[i] =  expectation(svmlikelihood,Normal(μ[i],sqrt(σ²[i])))
        end
    end
    return pred
end

###############################################################################

function local_updates!(model::VGP{BayesianSVM{T},<:AnalyticVI}) where {T<:Real}
    model.likelihood.ω .= broadcast((μ,Σ,y)->abs2.(one(T) .- y.*μ) + Σ ,model.μ,diag.(model.Σ),model.y)
    model.likelihood.θ .= broadcast(b->one(T)./sqrt.(b),model.likelihood.ω)
end

function local_updates!(model::SVGP{BayesianSVM{T},<:AnalyticVI}) where {T<:Real}
    model.likelihood.ω .= broadcast((κ,μ,Σ,y,K̃)->abs2.(one(T) .- y[model.inference.MBIndices].*(κ*μ)) + opt_diag(κ*Σ,κ) + K̃,model.κ,model.μ,model.Σ,model.y,model.K̃)
    model.likelihood.θ .= broadcast(b->one(T)./sqrt.(b),model.likelihood.ω)
end

""" Return the gradient of the expectation for latent GP `index` """
function cond_mean(model::VGP{BayesianSVM{T}},index::Integer) where {T<:Real}
    return model.y[index].*(model.likelihood.θ[index] .+ one(T))
end

function ∇μ(model::VGP{BayesianSVM{T}}) where {T<:Real}
    return broadcast((y,θ)->y.*(θ.+one(T)),model.y,model.likelihood.θ)
end

""" Return the gradient of the expectation for latent GP `index` """
function cond_mean(model::SVGP{BayesianSVM{T}},index::Integer) where {T<:Real}
    return model.y[index][model.inference.MBIndices].*(model.likelihood.θ[index].+one(T))
end

function ∇μ(model::SVGP{BayesianSVM{T}}) where {T<:Real}
    return broadcast((y,θ)->y[model.inference.MBIndices].*(θ.+one(T)),model.y,model.likelihood.θ)
end

function ∇Σ(model::AbstractGP{BayesianSVM{T}}) where {T<:Real}
    return model.likelihood.θ
end

function ELBO(model::AbstractGP{<:BayesianSVM})
    return expecLogLikelihood(model) - GaussianKL(model) - GIGEntropy(model)
end

function expecLogLikelihood(model::VGP{BayesianSVM{T},AnalyticVI{T}}) where {T<:Real}
    tot = -model.nLatent*(0.5*model.nSample*logtwo)
    tot += sum(broadcast((μ,y,θ,Σ)->(sum(μ.*y)-0.5*dot(θ,Σ+abs2.(one(T).-y.*μ))),
                        model.μ,model.y,model.likelihood.θ,diag.(model.Σ)))
    return tot
end

function expecLogLikelihood(model::SVGP{BayesianSVM{T},AnalyticVI{T}}) where {T<:Real}
    tot = -model.nLatent*(0.5*model.nSample*logtwo)
    tot += sum(broadcast((κμ,y,θ,κΣκ,K̃)->(sum(κμ.*y[model.inference.MBIndices])-0.5*dot(θ,K̃+κΣκ+abs2.(one(T).-y[model.inference.MBIndices].*κμ))),
                        model.κ.*model.μ,model.y,model.likelihood.θ,opt_diag.(model.κ.*model.Σ,model.κ),model.K̃))
    return model.inference.ρ*tot
end


text:
@testset "PARAFAC2" begin

Random.seed!(1)
r = 2
A = randn(r, 10)
T = [(randn(4 * (i + 1), r) * A) for i in 1:5]
@time factors = parafac2(T, 2, compute_error=true)
@test size.(factors.factors) == ntuple(i -> (size(T[i], 1), r), length(T))
@test size(factors.B) == (r, length(T))
@test size(factors.A) == (size(T[1], 2), r)
@test rel_residue(factors) < 0.05

end


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb04zd(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    NMAX = 20
    LDA = NMAX
    LDQG = NMAX
    LDU = NMAX
    LDWORK = ( NMAX+NMAX )*( NMAX+NMAX+1 )
    ZERO = 0.0e0
    ONE = 1.0e0
    A = Array{Float64,2}(undef, LDA,NMAX)
    QG = Array{Float64,2}(undef, LDQG,NMAX+1)
    U = Array{Float64,2}(undef, LDU,NMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    COMPU =  vs[2][1]
    if ( N<0 || N>NMAX )
        @error "Illegal N=$N"
        return
    end

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*(_jsz+1)÷2
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    _i0 = 0
    for j in 1:_jsz
       QG[j,j+1:_jsz+1] .= parsex.(Float64, vs[_i0+1:_i0+_jsz-j+1])
       _i0 += _jsz-j+1
    end
    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*(_jsz+1)÷2
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    _i0 = 0
    for j in 1:_jsz
       QG[j:_jsz,j] .= parsex.(Float64, vs[_i0+1:_i0+_jsz-j+1])
       _i0 += _jsz-j+1
    end
    close(f)
    # interp call 1

    INFO = SLICOT.mb04zd!(COMPU, N, A, QG, U)

    @test INFO == 0
    INFO != 0 && return

# interp output 1
    println(io, "A:")
    _nc = N
    _nr = N
    show(io, "text/plain", A[1:_nr,1:_nc])
    println(io)

# interp output 2
    println(io, "QG:")
    _nc = N
    _nr = N+1
    show(io, "text/plain", QG[1:_nr,1:_nc])
    println(io)

    Ht = zeros(2N,2N)
    for i in 1:N
        Ht[:,i] .= vcat(A[i,1:N],QG[1:i-1,i+1],QG[i,i+1:N+1])
        Ht[:,N+i] .= vcat(QG[i,1:i-1],QG[i:N,i],-A[1:N,i])
    end
    H = Matrix(Ht')
    println(io, "sq. reduced Hamiltonian:")
    show(io, "text/plain", H)
    println(io)
    # the original does this w/ BLAS in bits, but why bother?
    Hsq = H*H
    println(io, "squared sq. reduced Hamiltonian:")
    show(io, "text/plain", Hsq)
    println(io)
    llquadres = norm(Hsq[N+1:N,1:N])
    @test llquadres < 1e-12
    A2_hess_err =  norm(tril(Hsq[1:N,1:N],-2))
    @test A2_hess_err < 1e-12

end # run_mb04zd()


text:
struct Figure4
    θ::Vector{Float64}
    function Figure4(θ=[0.7,0.7,3.3,5,pi,4,2.7,0.55,0.26])
        new(θ)
    end
end

@recipe function f(f::Figure4)
    Ω = range(0,pi,length=200)
    Φ = range(0,2pi,length=200)
    model = FrequencyDirectionSpectra.FDJONSWAP(f.θ)
    layout := (1,3)
    @series begin
        subplot := 1
        xguide := "angular frequency"
        yguide := "direction"
        seriestype := :contourf
        seriescolor := :thermal
        levels := 30
        Ω, Φ, (ω,ϕ) -> FrequencyDirectionSpectra.freqdirspectra(model,ω,ϕ)
    end
    @series begin
        subplot := 2
        label := false
        xguide := "angular frequency"
        yguide := "spectral density"
        range(0,pi,length=1000), ω -> FrequencyDirectionSpectra.sdf(model, ω)
    end
    @series begin
        subplot := 3
        xguide := "angular frequency"
        yguide := "direction"
        seriestype := :contourf
        seriescolor := :acton
        levels := 30
        Ω, Φ, (ω,ϕ) -> FrequencyDirectionSpectra.spreading(model,ω,ϕ)
    end
end


text:
using Documenter
using GradientStarvation.Data
using GradientStarvation.Figures
using GradientStarvation.Experiments
using GradientStarvation.NeuralNetwork
using GradientStarvation.PlotLib
using GradientStarvation.PlotUtils
using GradientStarvation.Regularization

makedocs(
    sitename = "Gradient Starvation",
    format = Documenter.HTML(),
    modules = [
        GradientStarvation.Data,
        GradientStarvation.Figures,
        GradientStarvation.Experiments,
        GradientStarvation.NeuralNetwork,
        GradientStarvation.PlotLib,
        GradientStarvation.PlotUtils,
        GradientStarvation.Regularization
    ],
    pages = Any[
        "Data" => "man/data.md",
        "Figures" => "man/figures.md",
        "Experiments" => "man/experiments.md",
        "Neural Network" => "man/neuralnetwork.md",
        "Plots" => "man/plotlib.md",
        "Plot Utils" => "man/plotutils.md",
        "Regularization" => "man/regularization.md"
    ]
)

#=
deploydocs(
            repo = "github.com/justinsane1337/GradientStarvation.git",
        )
=#


text:
# Overall structure
#----------------------------

import ArgParse: parse_item
export parse_commandline, print_parsed_args

function arg_settings()
    s = ArgParseSettings()

    add_arg_group!(s, "Computation", "comp")
    add_arg_group!(s, "Payoffs", "payoff")
    add_arg_group!(s, "Approximation", "approx")
    add_arg_group!(s, "Parameters", "param")

    @add_arg_table! s begin

        # Dataset
        "--dataset", "-d"
            help = "Dataset to use"
            arg_type = String
            default = "data_all_leases.RData"

        # Computation control
        "--Mcnstr"
            arg_type = Int
            default = 500
            group = "comp"
        "--Mfull"
            arg_type = Int
            default  = 2000
            group = "comp"

        "--doCnstr"
            action = :store_true
            group = "comp"
        "--noFull"
            action = :store_true
            group = "comp"

        "--maxtimeCnstr"
            arg_type = Int
            default = 0 * 60^2
            group = "comp"
        "--maxtimeFull"
            arg_type = Int
            default  = 0 * 60^2
            group = "comp"
        "--noPar"
            action = :store_true
            group = "comp"

        # Parameters
        "--computeStarting"
            action = :store_true
            group = "param"
        "--theta"
            arg_type = Vector{Float64}
            default = zeros(0)
            group = "param"
        "--convertKappa"
            action = :store_true
            group = "param"

        # Payoffs
        "--cost"
            arg_type = AbstractDrillingCost
            default = DrillingCost_TimeFE(2008.5,2012.5)
            group = "payoff"
        "--extension"
            arg_type = AbstractExtensionCost
            default = ExtensionCost_Constant()
            group = "payoff"
        "--revenue"
            arg_type = DrillingRevenue
            default = DrillingRevenue(Unconstrained(), TimeTrend(), GathProcess() )
            # default = DrillingRevenue(Unconstrained(), TimeFE(2008.5,2016.5), GathProcess() )
            group = "payoff"
        "--scrap"
            arg_type = AbstractScrapValue
            default = ScrapValue_Zero()
            group = "payoff"

        "--anticipateT1EV"
            help = "T1ev shocks anticipated"
            action = :store_true
            group = "payoff"
        "--discount"
            arg_type = Float64
            default = RealDiscountRate()
            group = "payoff"

        # Approximation
        "--numP"
            arg_type=Int
            default=13
            group="approx"
        "--numPsi"
            arg_type=Int
            default=13
            group="approx"
        "--extendPriceGrid"
            arg_type=Float64
            default=log(2)
            group="approx"
        "--minTransProb"
            arg_type=Float64
            default=minp_default()
            group="approx"
    end
    return s
end

parse_commandline() = parse_args(arg_settings())

ArgParse.parse_item(::Type{AbstractDrillingCost}, x::AbstractString) =
    eval(parse(x))::AbstractDrillingCost

ArgParse.parse_item(::Type{AbstractExtensionCost}, x::AbstractString) =
    eval(parse(x))::AbstractExtensionCost

ArgParse.parse_item(::Type{AbstractScrapValue}, x::AbstractString) =
    eval(parse(x))::AbstractScrapValue

ArgParse.parse_item(::Type{DrillingRevenue}, x::AbstractString) =
    eval(parse(x))::DrillingRevenue

ArgParse.parse_item(::Type{Vector{T}}, x::AbstractString) where {T} =
    eval(parse(x))::Vector{T}

ArgParse.parse_item(::Type{T}, x::AbstractString) where {T<:Real} =
    eval(parse(x))::T

function print_parsed_args(x::Dict)
    dset = ["dataset",]
    ctrl = [
        "Mcnstr",
        "Mfull",
        "doCnstr",
        "noFull",
        "noPar",
        "maxtimeCnstr",
        "maxtimeFull",
        "computeStarting",
    ]
    parm = [
        "computeStarting",
        "theta",
        "convertKappa"
    ]
    payoff = [
        "revenue",
        "cost",
        "extension",
        "scrap",
        "anticipateT1EV",
        "discount",
    ]
    approx = [
        "numP",
        "numPsi",
        "extendPriceGrid",
        "minTransProb",
    ]

    kys = vcat(dset, ctrl, parm, payoff, approx)
    padlen = maximum(length.(kys))
    padk(k) = rpad(k,padlen)

    println("--------------------------------------")
    println("ARGUMENTS PASSED")
    println("--------------------------------------")
    println("DATASET")
    map(k -> println("  $(padk(k))  =>  $(x[k])"), dset)
    println("ESTIMATION CONTROL")
    map(k -> println("  $(padk(k))  =>  $(x[k])"), ctrl)
    println("PAYOFF")
    map(k -> println("  $(padk(k))  =>  $(x[k])"), payoff)
    println("APPROXIMATION")
    map(k -> println("  $(padk(k))  =>  $(x[k])"), approx)
    println("--------------------------------------")
    print("\n\n")
    return nothing
end


text:
#=
Utility functions for Str strings

Copyright 2018 Gandalf Software, Inc., Scott P. Jones, and other contributors to the Julia language
Licensed under MIT License, see LICENSE.md
Based initially on julia/test/strings/util.jl
=#

# starts with and ends with predicates

starts_with(a::MaybeSub{<:Str{C}}, b::MaybeSub{<:Str{C}}) where {C<:CSE} =
    (len = ncodeunits(b)) <= ncodeunits(a) &&
    (@preserve a b _memcmp(pointer(a), pointer(b), len)) == 0

ends_with(a::MaybeSub{<:Str{C}}, b::MaybeSub{<:Str{C}}) where {C<:CSE} =
    (lenb = ncodeunits(b)) <= (lena = ncodeunits(a)) &&
    (@preserve a b _memcmp(pointer(a) + lena - lenb, pointer(b), lenb)) == 0

@static if false
function chop(s::AbstractString; head::Integer = 0, tail::Integer = 1)
    SubString(s, nextind(s, firstindex(s), head), prevind(s, lastindex(s), tail))
end
end # if false

function chomp(str::MaybeSub{<:Str})
    (len = ncodeunits(str)) != 0 && @preserve str begin
        pnt = pointer(str)
        get_codeunit(pnt, len) == 0xa &&
            (len -= (len > 1 && get_codeunit(pnt, len - 1) == 0x0d) + 1)
    end
    SubString(str, 1, thisind(str, len))
end

## string padding functions ##

# Todo: make these not build two strings, and return the correct type (i.e. of the first argument)

function _lpad(cnt, pad, str)
    cnt, rem = divrem(cnt, length(pad))
    rem == 0 ? string(pad^cnt, str) : string(pad^cnt, first(pad, rem), str)
end
lpad(str::MaybeSub{<:Str}, cnt::Integer, pad::AbstractString) =
    (cnt -= length(str)) <= 0 ? str : _lpad(cnt, pad, str)
lpad(ch::Chr, cnt::Integer, pad::AbstractString) =
    (cnt -= 1) <= 0 ? string(ch) : _lpad(cnt, pad, ch)
lpad(str::MaybeSub{<:Str}, cnt::Integer, pad::AbsChar=' ') =
    (cnt -= length(str)) <= 0 ? str : string(pad^cnt, str)
lpad(ch::Chr, cnt::Integer, pad::AbstractChar=' ') =
    (cnt -= 1) <= 0 ? string(ch) : string(pad^cnt, ch)

function _rpad(cnt, pad, str)
    cnt, rem = divrem(cnt, length(pad))
    rem == 0 ? string(str, pad^cnt) : string(str, pad^cnt, first(pad, rem))
end
rpad(str::MaybeSub{<:Str}, cnt::Integer, pad::AbstractString) =
    (cnt -= length(str)) <= 0 ? str : _rpad(cnt, pad, str)
rpad(ch::Chr, cnt::Integer, pad::AbstractString) =
    (cnt -= 1) <= 0 ? string(ch) : _rpad(cnt, pad, ch)
rpad(str::MaybeSub{<:Str}, cnt::Integer, pad::AbsChar=' ') =
    (cnt -= length(str)) <= 0 ? str : string(str, pad^cnt)
rpad(ch::Chr, cnt::Integer, pad::AbsChar=' ') =
    (cnt -= 1) <= 0 ? string(ch) : string(ch, pad^cnt)

const SetOfChars = Union{Tuple{Vararg{<:AbsChar}},AbstractVector{<:AbsChar},Set{<:AbsChar}}

function __split(str, splitter, limit::Integer, keep_empty::Bool, strs::Vector)
    pos = 1
    lst = lastindex(str)
    res = find(First, splitter, str)
    if res != 0:-1
        j, k = first(res), nextind(str, last(res))
        while 0 < j <= lst && length(strs) != limit - 1
            if pos < k
                (keep_empty || pos < j) &&
                    push!(strs, SubString(str, pos, prevind(str, j)))
                pos = k
            end
            (k <= j) && (k = nextind(str, j))
            res = find(Fwd, splitter, str, k)
            res == 0:-1 && break
            j, k = first(res), nextind(str, last(res))
        end
    end
    (keep_empty || pos <= lst) ? push!(strs, SubString(str, pos)) : strs
end

function __rsplit(str, splitter, limit::Integer, keep_empty::Bool, strs::Vector)
    res = find(Last, splitter, str)
    j, k = first(res), last(res)
    pos = lastindex(str)
    while j > 0 && k > 0 && length(strs) != limit-1
        (keep_empty || k < pos) && pushfirst!(strs, SubString(str, nextind(str, k), pos))
        (pos = prevind(str, j)) > 0 || break
        res = find(Rev, splitter, str, pos)
        j, k = first(res), last(res)
    end
    (keep_empty || pos > 0) && pushfirst!(strs, SubString(str, 1, pos))
    strs
end

function checkkeep(keepempty, keep, fun)
    keep === nothing && return keepempty
    Base.depwarn("The `keep` keyword argument is deprecated; use `keepempty` instead", fun)
end

splitarr(::Type{C}) where {C<:CSE} = SubString{Str{basecse(C),Nothing,Nothing,Nothing}}[]

splitarr(::MaybeSub{String}) = SubString{String}[]
splitarr(::MaybeSub{T}) where {T<:Str} =
    SubString{Str{basecse(T),Nothing,Nothing,Nothing}}[]

Base._split(str::MaybeSub{<:Str}, splitter, limit, keepempty, vec) =
    __split(str, splitter, limit, keepempty, vec)
Base._rsplit(str::MaybeSub{<:Str}, splitter, limit, keepempty, vec) =
    __rsplit(str, splitter, limit, keepempty, vec)

split(str::MaybeSub{<:Str}, splitter;
      limit::Integer=0, keepempty::Bool=true, keep::Union{Nothing,Bool}=nothing) =
    __split(str, splitter, limit, checkkeep(keepempty, keep, :split), splitarr(str))

split(str::MaybeSub{<:Str}, splitter::AbstractChar;
      limit::Integer=0, keepempty::Bool=true, keep::Union{Nothing,Bool}=nothing) =
    __split(str, isequal(splitter), limit, checkkeep(keepempty, keep, :split), splitarr(str))

split(str::MaybeSub{<:Str}, splitter::SetOfChars;
      limit::Integer=0, keepempty::Bool=true, keep::Union{Nothing,Bool}=nothing) =
    __split(str, in(splitter), limit, checkkeep(keepempty, keep, :split), splitarr(str))

rsplit(str::MaybeSub{<:Str}, splitter;
       limit::Integer=0, keepempty::Bool=true, keep::Union{Nothing,Bool}=nothing) =
    __rsplit(str, splitter, limit, checkkeep(keepempty, keep, :rsplit), splitarr(str))
rsplit(str::MaybeSub{<:Str}, splitter::AbstractChar;
       limit::Integer=0, keepempty::Bool=true, keep::Union{Nothing,Bool}=nothing) =
    __rsplit(str, isequal(splitter), limit, checkkeep(keepempty, keep, :rsplit), splitarr(str))
rsplit(str::MaybeSub{<:Str}, splitter::SetOfChars;
       limit::Integer=0, keepempty::Bool=true, keep::Union{Nothing,Bool}=nothing) =
    __rsplit(str, in(splitter), limit, checkkeep(keepempty, keep, :rsplit), splitarr(str))

# Todo: this is using print, but it should be changed to make sure that everything is done via
# writes (i.e. no translation to UTF-8)
function __replace(str, pat_repl::Pair; count::Integer=typemax(Int))
    pattern, repl = pat_repl
    count == 0 && return str
    count < 0 && throw(DomainError(count, "`count` must be non-negative."))
    pos = 1
    lst = lastindex(str)
    res = find(First, pattern, str)
    # Just return the string if not found
    (j = first(res)) == 0 && return str

    out = get_iobuffer(floor(Int, 1.2 * sizeof(str)))

    while true
        k = last(res)
        if pos == 1 || pos <= k
            print(out, SubString(str, pos, thisind(str, j - 1)))
            #unsafe_write(out, pointer(str, pos), UInt(j - pos))
            Base._replace(out, repl, str, res, pattern)
        end
        if k < j
            pos = j
            j > lst && break
            k = nextind(str, j)
        else
            pos = k = nextind(str, k)
        end
        (count -= 1) > 0 || break
        res = find(Fwd, pattern, str, k)
        (j = first(res)) == 0 && break
    end
    print(out, SubString(str, pos))
    convert(Str{basecse(str),Nothing,Nothing,Nothing}, String(take!(out)))
end

replace(str::MaybeSub{<:Str}, pat_repl::Pair; count::Integer=typemax(Int)) =
    __replace(str, pat_repl; count=count)
replace(str::MaybeSub{<:Str}, pat_repl::Pair{<:AbstractChar}; count::Integer=typemax(Int)) =
    __replace(str, ==(first(pat_repl)) => last(pat_repl); count=count)
replace(str::MaybeSub{<:Str}, pat_repl::Pair{<:SetOfChars}; count::Integer=typemax(Int)) =
    __replace(str, in(first(pat_repl)) => last(pat_repl); count=count)


text:
using RingArrays
using FactCheck

if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

include("RingArrays.jl")

FactCheck.exitstatus()


text:
using FunctionalCollections: PersistentVector, assoc, push, pop

###########################################
# trace used by vector-shaped combinators #
###########################################

"""

U is the type of the subtrace, R is the return value type for the kernel
"""
struct VectorTrace{GenFnType,T,U} <: Trace
    gen_fn::GenerativeFunction
    subtraces::PersistentVector{U}
    retval::PersistentVector{T}
    args::Tuple

    # number of active subtraces (may be less than length of subtraces)
    len::Int

    # number of active subtraces that are nonempty (used for has_choices)
    num_nonempty::Int

    score::Float64
    noise::Float64
end

function VectorTrace{GenFnType,T,U}(gen_fn::GenerativeFunction,
                                    subtraces::PersistentVector{U},
                                    retval::PersistentVector{T},
                                    args::Tuple, score::Float64, noise::Float64,
                                    len::Int, num_nonempty::Int) where {GenFnType,T,U}
    @assert length(subtraces) == length(retval)
    @assert length(subtraces) == len
    @assert num_nonempty >= 0
    VectorTrace{GenFnType,T,U}(gen_fn, subtraces, retval, args, len,
        num_nonempty, score, noise)
end

function VectorTrace{GenFnType,T,U}(gen_fn::GenerativeFunction, args::Tuple) where {GenFnType,T,U}
    subtraces = PersistentVector{U}()
    retvals = PersistentVector{T}()
    VectorTrace{GenFnType,T,U}(gen_fn, subtraces, retvals, args, 0, 0, 0., 0.)
end

# trace API

@inline get_choices(trace::VectorTrace) = VectorTraceChoiceMap(trace)
@inline get_retval(trace::VectorTrace) = trace.retval
@inline get_args(trace::VectorTrace) = trace.args
@inline get_score(trace::VectorTrace) = trace.score
@inline get_gen_fn(trace::VectorTrace) = trace.gen_fn

@inline function Base.getindex(trace::VectorTrace{GenFnType, T, U}, addr::Pair) where {GenFnType, T, U}
    (first, rest) = addr
    subtrace = trace.subtraces[first]
    subtrace[rest]
end

@inline function Base.getindex(trace::VectorTrace, addr)
    # we expose the return values in the auxiliary state
    trace.retval[addr]
end

function project(trace::VectorTrace, selection::Selection)
    weight = 0.
    for key=1:trace.len
        subselection = selection[key]
        weight += project(trace.subtraces[key], subselection)
    end
    weight
end
project(trace::VectorTrace, ::EmptySelection) = trace.noise

struct VectorTraceChoiceMap{GenFnType, T, U} <: ChoiceMap
    trace::VectorTrace{GenFnType, T, U}
end

@inline Base.isempty(assignment::VectorTraceChoiceMap) = assignment.trace.num_nonempty == 0
@inline get_address_schema(::Type{VectorTraceChoiceMap}) = VectorAddressSchema()

@inline function get_submap(choices::VectorTraceChoiceMap, addr::Int)
    if addr <= choices.trace.len
        get_choices(choices.trace.subtraces[addr])
    else
        EmptyChoiceMap()
    end
end

@inline function get_submaps_shallow(choices::VectorTraceChoiceMap)
    ((i, get_choices(choices.trace.subtraces[i])) for i=1:choices.trace.len)
end

@inline get_submap(choices::VectorTraceChoiceMap, addr::Pair) = _get_submap(choices, addr)
@inline get_value(choices::VectorTraceChoiceMap, addr::Pair) = _get_value(choices, addr)
@inline has_value(choices::VectorTraceChoiceMap, addr::Pair) = _has_value(choices, addr)
@inline get_values_shallow(::VectorTraceChoiceMap) = ()

############################################
# code shared by vector-shaped combinators #
############################################

function get_retained_and_constrained(constraints::ChoiceMap, prev_length::Int, new_length::Int)
    keys = Set{Int}()
    for (key::Int, _) in get_submaps_shallow(constraints)
        if key > 0 && key <= new_length
            push!(keys, key)
        else
            error("Constrained address does not exist: $key")
        end
    end
    keys
end

function get_retained_and_selected(selection::EmptySelection, prev_length::Int, new_length::Int)
    Set{Int}()
end

function get_retained_and_selected(selection::HierarchicalSelection, prev_length::Int, new_length::Int)
    keys = Set{Int}()
    for (key::Int, _) in get_subselections(selection)
        if key > 0 && key <= new_length
            push!(keys, key)
        else
            error("Selected address does not exist: $key")
        end
    end
    keys
end

function get_retained_and_selected(selection::AllSelection, prev_length::Int, new_length::Int)
    Set{Int}(1:min(prev_length, new_length))
end

function vector_compute_retdiff(updated_retdiffs::Dict{Int,Diff}, new_length::Int, prev_length::Int)
    if new_length == prev_length && length(updated_retdiffs) == 0
        NoChange()
    else
        VectorDiff(new_length, prev_length, updated_retdiffs)
    end
end

function vector_update_delete(new_length::Int, prev_length::Int,
                                 prev_trace::VectorTrace)
    num_nonempty = prev_trace.num_nonempty
    discard = choicemap()
    score_decrement = 0.
    noise_decrement = 0.
    for key=new_length+1:prev_length
        subtrace = prev_trace.subtraces[key]
        score_decrement += get_score(subtrace)
        noise_decrement += project(subtrace, EmptySelection())
        if !isempty(get_choices(subtrace))
            num_nonempty -= 1
        end
        @assert num_nonempty >= 0
        set_submap!(discard, key, get_choices(subtrace))
    end
    return (discard, num_nonempty, score_decrement, noise_decrement)
end

function vector_regenerate_delete(new_length::Int, prev_length::Int,
                                  prev_trace::VectorTrace)
    num_nonempty = prev_trace.num_nonempty
    score_decrement = 0.
    noise_decrement = 0.
    for key=new_length+1:prev_length
        subtrace = prev_trace.subtraces[key]
        score_decrement += get_score(subtrace)
        noise_decrement += project(subtrace, EmptySelection())
        if !isempty(get_choices(subtrace))
            num_nonempty -= 1
        end
        @assert num_nonempty >= 0
    end
    return (num_nonempty, score_decrement, noise_decrement)
end

function vector_remove_deleted_applications(subtraces, retval, prev_length, new_length)
    for i=new_length+1:prev_length
        subtraces = pop(subtraces)
        retval = pop(retval)
    end
    (subtraces, retval)
end

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
# path grath
g = DAG(:A => :B, :B => :C)
@test is_d_separated(g, [1], [3], [2])
@test !is_d_separated(g, [1], [3], Int[])

# fork graph
g = DAG(:A => :B, :A => :C)
@test is_d_separated(g, [:B], [:C], [:A])
@test !is_d_separated(g, [:B], [:C], Symbol[])

# collider graph
g = DAG(:A => :B, :C => :B)
@test is_d_separated(g, [:A], [:C], Symbol[])
@test !is_d_separated(g, [:A], [:C], [:B])

# asia graph
g = DAG(
    :asia => :tuberculosis,
    :smoking => :cancer,
    :smoking => :bronchitis,
    :tuberculosis => :either,
    :cancer => :either,
    :either => :xray,
    :either => :dyspnea,
    :bronchitis => :dyspnea
)
# this test is wrong - dyspnea is in isolated component, but xray is on the main
@test is_d_separated(g,
    [:asia, :smoking],
    [:dyspnea, :xray],
    [:bronchitis, :either]
)
@test is_d_separated(g,
    [:tuberculosis, :cancer],
    [:bronchitis],
    [:smoking, :xray]
)

# test ignore_edges argument
g = DAG(:A => :B, :A => :C, :C => :B)
@test !is_d_separated(g, [:A], [:B], [:C])
@test is_d_separated(g, [:A], [:B], [:C], ignore_edges=[:A => :B])