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, :(1/(1 + exp(-(-0.3402051397109518 + -0.28456924845213605*1/(1 + exp(-(-0.42610478134649465 + 0.43223652991075534*1/(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*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(-0.17528834888704425 + -0.1425247358809414*1/(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*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(-0.005280183583000753 + 0.6172883711874637*1/(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*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(0.26121381489212014 + -0.018200625348021138*1/(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*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(-0.09690410604182231 + 0.6473593808554208*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5])))))))))) + 1/(1 + exp(-(0.2351041904461102 + 0.11979579421065711*1/(1 + exp(-(-0.42610478134649465 + 0.43223652991075534*1/(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*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(-0.17528834888704425 + -0.1425247358809414*1/(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*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(-0.005280183583000753 + 0.6172883711874637*1/(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*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(0.26121381489212014 + -0.018200625348021138*1/(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*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(-0.09690410604182231 + 0.6473593808554208*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5])))))))))) + 1/(1 + exp(-(-0.4731175976624935 + -0.5144698215629462*1/(1 + exp(-(-0.42610478134649465 + 0.43223652991075534*1/(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*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(-0.17528834888704425 + -0.1425247358809414*1/(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*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(-0.005280183583000753 + 0.6172883711874637*1/(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*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(0.26121381489212014 + -0.018200625348021138*1/(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*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(-0.09690410604182231 + 0.6473593808554208*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5])))))))))) + 1/(1 + exp(-(-0.8380258290664089 + 0.49770714260784343*1/(1 + exp(-(-0.42610478134649465 + 0.43223652991075534*1/(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*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(-0.17528834888704425 + -0.1425247358809414*1/(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*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(-0.005280183583000753 + 0.6172883711874637*1/(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*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(0.26121381489212014 + -0.018200625348021138*1/(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*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(-0.09690410604182231 + 0.6473593808554208*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5])))))))))) + 1/(1 + exp(-(0.7168755118620993 + -0.38856981611164576*1/(1 + exp(-(-0.42610478134649465 + 0.43223652991075534*1/(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*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(-0.17528834888704425 + -0.1425247358809414*1/(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*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(-0.005280183583000753 + 0.6172883711874637*1/(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*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(0.26121381489212014 + -0.018200625348021138*1/(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*1/(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*1/(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*1/(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*1/(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*1/(1 + exp(-(-0.09690410604182231 + 0.6473593808554208*1/(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*1/(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*1/(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*1/(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*1/(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:
## Si-Midorikawa 1999 GMPE tests
@testset "Si-Midorikawa 1999 GMPE PGA" begin
  # init model parameters
  include("../examples/si-midorikawa-1999.conf")
  ## test at epicenter on grid M7.0
  @test gmpe_simidorikawa1999(eq_7,config_simidorikawa1999_crustal_pga,grid_epicenter)[1].pga == 59.04
  ## run PGA modeling on grid withoit minpga Depth <= 30 M6.0
  S_c = gmpe_simidorikawa1999(eq_6,config_simidorikawa1999_crustal_pga,grid)
  @test length(S_c) == TEST_GRID_SIZE
  @test round(sum([S_c[i].pga for i=1:length(S_c)]),digits=2) == 6.98
  S_intp = gmpe_simidorikawa1999(eq_6,config_simidorikawa1999_interplate_pga,grid)
  @test length(S_intp) == TEST_GRID_SIZE
  @test round(sum([S_intp[i].pga for i=1:length(S_intp)]),digits=2) == 8.38
  S_intra = gmpe_simidorikawa1999(eq_6,config_simidorikawa1999_intraplate_pga,grid)
  @test length(S_intra) == TEST_GRID_SIZE
  @test round(sum([S_intra[i].pga for i=1:length(S_intra)]),digits=2) == 13.9
  ## run PGA modeling on grid with minpga Depth <= 30 M6.0
  S_c = gmpe_simidorikawa1999(eq_6,config_simidorikawa1999_crustal_pga,grid,min_val=0.34)
  @test length(S_c) == WITH_MINPGA
  @test round(sum([S_c[i].pga for i=1:length(S_c)]),digits=2) == 4.61
  ## run PGA modeling on grid with minpga Depth > 30 M6.0
  eq_30 = Earthquake(143.04,51.92,35,6.0)
  S_c = gmpe_simidorikawa1999(eq_30,config_simidorikawa1999_crustal_pga,grid,min_val=0.15)
  @test length(S_c) == WITH_MINPGA
  @test round(sum([S_c[i].pga for i=1:length(S_c)]),digits=2) == 2.04
  ## run PGA modeling for plotting Depth <= 30 M6.0
  S_c = gmpe_simidorikawa1999(eq_6,config_simidorikawa1999_crustal_pga)
  @test length(S_c) == SIMULATION_ARRAY_SIZE
  @test round(sum(S_c),digits=2) == 1096.79
  S_intp = gmpe_simidorikawa1999(eq_6,config_simidorikawa1999_interplate_pga)
  @test round(sum(S_intp),digits=2) == 1318.74
  S_intra = gmpe_simidorikawa1999(eq_6,config_simidorikawa1999_intraplate_pga)
  @test round(sum(S_intra),digits=2) == 2188.89
  ## run PGA modeling for plotting Depth > 30 M6.0
  S_c = gmpe_simidorikawa1999(eq_30,config_simidorikawa1999_crustal_pga)
  @test round(sum(S_c),digits=2) == 722.93
end

ext: jl


text:
# ---
# title: 329. Longest Increasing Path in a Matrix
# id: problem329
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Depth-first Search, Topological Sort, Memoization
# link: <https://leetcode.com/problems/longest-increasing-path-in-a-matrix/description/>
# hidden: true
# ---
#
# Given an integer matrix, find the length of the longest increasing path.
#
# From each cell, you can either move to four directions: left, right, up or
# down. You may NOT move diagonally or move outside of the boundary (i.e. wrap-
# around is not allowed).
#
# **Example 1:**
#
#
#
#     Input: nums =
#     [
#       [9,9,4],
#       [6,6,8],
#       [2,1,1]
#     ]
#     Output: 4
#     Explanation: The longest increasing path is [1, 2, 6, 9].
#
#
# **Example 2:**
#
#
#
#     Input: nums =
#     [
#       [3,4,5],
#       [3,2,6],
#       [2,2,1]
#     ]
#     Output: 4
#     Explanation: The longest increasing path is [3, 4, 5, 6]. Moving diagonally is not allowed.
#
#
#
## @lc code=start
using LeetCode

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


text:
using ApproxFun, SpectralMeasures

# randn and 1/n^2 decay gives good error bounds
n = 50
a = randn(n)./(1:n).^2
b = max(.05,.5+randn(n-1)./(1:n-1).^2)
spectrum, maxerr = validatedspectrum(a,b)
maxerr

# randn and 1/n decay gives okay error bounds
n = 50
a = randn(n)./(1:n)
b = max(.05,.5+randn(n-1)./(1:n-1))
spectrum, maxerr = validatedspectrum(a,b)
maxerr

# randn with no decay gives terrible error bounds
n = 20
a = randn(n)
b = max(.05,.5+randn(n-1))
spectrum, maxerr = validatedspectrum(a,b)
maxerr


text:
module TropicalCycloneRecords

using HTTP: get
using StringEncodings: decode
using Dates: DateTime

# Write your package code here.
include("download_data.jl")

end


text:
export fmtstar!

# status = 'U', 'O', or 'C' (unvisited, open, closed)
const FMTNodeInfo{X,D} = NamedTuple{(:is_free, :status, :parent, :cost_to_come),Tuple{Bool,Char,X,D}}
function FMTNodeInfo{X,D}(; is_free=false, status='U', parent=zero(X), cost_to_come=zero(D)) where {X,D}
    (is_free=is_free, status=status, parent=parent, cost_to_come=cost_to_come)
end

function fmtstar!(P::MPProblem; r, ensure_goal_count=0, compute_full_metadata=true)
    P.status, P.solution = fmtstar!(P.state_space, P.bvp, P.init, P.goal, P.collision_checker,
                                    isdefined(P, :graph) ? P.graph : error("TODO"), r, ensure_goal_count)
    if P.status === :solved && compute_full_metadata
        metadata = P.solution.metadata

        X = indextype(P.graph)
        node_info = metadata[:node_info]
        metadata[:tree] = Dict{X,X}(k => node_info[k].parent for k in keys(node_info) if node_info[k].status != 'U')
    end
    P.status, P.solution
end

function fmtstar!(state_space::StateSpace,
                  bvp::SteeringBVP,
                  init::State,
                  goal::Goal,
                  collision_checker::CollisionChecker,
                  graph::NearNeighborGraph{NeighborInfo{X,D,U}},
                  r, ensure_goal_count=0) where {X,D,U}
    tic = time()
    collision_checker.motion_count[] = collision_checker.edge_count[] = 0    # reset!(collision_checker)

    if !is_free_state(collision_checker, init)
        warn("Initial state is infeasible!")
        status = :infeasible
        solution = MPSolution(status, D(Inf), time() - tic, Dict{Symbol,Any}())
        return status, solution
    end

    setinit!(graph, init, r)
    setgoal!(graph, () -> rand_free_state(collision_checker, goal), r, count=ensure_goal_count)
    node_info = node_info_datastructure(graph.nodes, FMTNodeInfo{X,D},
        x -> FMTNodeInfo{X,D}(is_free=is_free_state(collision_checker, x) && x in state_space))
    is_unvisited = let node_info=node_info; n -> (ni = node_info[n.index]; ni.is_free && ni.status == 'U'); end
    is_open      = let node_info=node_info; n -> node_info[n.index].status == 'O';  end    # JuliaLang/julia#15276

    open_queue = PriorityQueue{X,D}()
    z = X(0)
    node_info[z] = (node_info[z]..., status='O')

    while !(graph[z] in goal)
        for (x, _, _) in neighbors(is_unvisited, graph, z, r, dir=Val(:F))
            isempty(neighbors(is_open, graph, x, r, dir=Val(:B))) && @show (z, x, isempty(neighbors(is_open, graph, x, r, dir=Val(:B))))
            y_min, c_min, u_min = mapreduce(let node_info=node_info    # JuliaLang/julia#15276
                                                n -> (index=n.index,
                                                      cost=node_info[n.index].cost_to_come + n.cost,
                                                      controls=n.controls)
                                            end,
                                            (a, b) -> ifelse(a.cost < b.cost, a, b),
                                            neighbors(is_open, graph, x, r, dir=Val(:B)))
            if is_free_edge(collision_checker, bvp, graph[y_min], graph[x], u_min)
                node_info[x] = (is_free=true, status='O', parent=y_min, cost_to_come=c_min)
                open_queue[x] = c_min
            end
        end
        node_info[z] = (node_info[z]..., status='C')
        isempty(open_queue) ? break : z = dequeue!(open_queue)
    end

    solution_nodes = [z]
    costs_to_come = [node_info[z].cost_to_come]
    while linear_index(solution_nodes[1]) != 0
        pushfirst!(solution_nodes, node_info[solution_nodes[1]].parent)
        pushfirst!(costs_to_come, node_info[solution_nodes[1]].cost_to_come)
    end

    status = graph[z] in goal ? :solved : :failed
    solution_metadata = Dict{Symbol,Any}(
        :r => r,
        :motion_checks => collision_checker.motion_count[],
        :edge_checks => collision_checker.edge_count[],
        :cost => costs_to_come[end],
        :costs_to_come => costs_to_come,
        :planner => :FMTstar,
        :solved => status === :solved,
        :solution_nodes => solution_nodes,
        :node_info => node_info
    )
    solution = MPSolution(status, costs_to_come[end], time() - tic, solution_metadata)
    status, solution
end

max_stars_repo_name: pbpf/clean

max_issues_repo_name: pbpf/clean

max_forks_repo_name: pbpf/clean


text:
module TestSorting

using Test
using RLEVectors

@testset begin

# issorted
@test issorted( RLEVector([3,4,2],[3,6,9])) == false
@test issorted( RLEVector([2,3,4],[3,6,9])) == true
@test issorted( RLEVector(["bob","fred","joe"],[3,6,9])) == true

# sort, sort!
x = RLEVector([2,4,3,5],[4,7,9,10])
@test sort(x) == RLEVector([2,3,4,5],[4,6,9,10])
@test sort!(x) == RLEVector([2,3,4,5],[4,6,9,10])

# sortperm
x = RLEVector([2,4,3,5],[4,7,9,10])
@test sortperm(x) == RLEVector([1,3,2,4],[4,6,9,10])

# reverse
x = RLEVector([2,4,3,5],[4,7,9,10])
@test reverse(x) == RLEVector([5,3,4,2],[1,3,6,10])
@test reverse!(x) == RLEVector([5,3,4,2],[1,3,6,10])

# permute runs
x = RLEVector([2,4,3,5],[4,7,9,10])
y = RLEVector([4,5,3,2],[3,4,6,10])
@test permute_runs(x, [2, 4, 3, 1]) == y

end # testset

end # module


text:
#####
##### Function to parse hitran par-file data
#####

using ..Architectures: GPU

"""
    read_hitran(filepath::String, mol::Int=-1, iso::Int=-1, ν_min::Real=0, ν_max::Real=Inf)

Read/parse a HITRAN data file and return the data in [`HitranTable`](@ref) format

"""
function read_hitran(filepath::String; mol::Int=-1, iso::Int=-1,
                     ν_min::Real=0, ν_max::Real=Inf,
                     min_strength::Real=0)

    # Infer type from input ν_min
    FT = eltype(AbstractFloat(ν_min))

    # Declare some constant properties of the hitran line parameters, such as
    # names, lengths, and types
    varNames = ["molec_id", "local_iso_id", "nu", "sw", "a", "gamma_air",
                "gamma_self", "elower", "n_air", "delta_air", "global_upper_quanta",
                "global_lower_quanta", "local_upper_quanta", "local_lower_quanta",
                "ierr", "iref", "line_mixing_flag", "gp", "gpp"]
    varLengths = [2, 1, 12, 10, 10, 5, 5, 10, 4, 8, 15, 15, 15, 15, 6, 12, 1, 7, 7]
    varTypes = [Int64, Int64, FT, FT, FT, FT, FT, FT, FT, FT, String, String, String, String, String, String, String, FT, FT]

    # Take the line parameters' lengths and assemble an index-ranges list
    # Each value is the starting index, minus 1, of var i and ending index of var i-1
    # i.e.
    #   [0, 2, 3, 15, 25, 35, 40, 45, 55, 59, 67, 82, 97, 112, 127, 133, 145, 146, 153, 160]
    # Starting at 0 correctly produces (1,2), (3,3), (4,15), ..., (154, 160)
    idxRanges = append!([0], [sum(varLengths[1:i]) for i in 1:length(varLengths)])

    # Open the file specified
    open(filepath, "r") do file

        # Where to hold the matching rows
        rows = []

        # Loop through every line
        for ln in eachline(file)

            # Go from a line to a list of values
            values = [varTypes[i] == String ? ln[(idxRanges[i]+1):idxRanges[i+1]] : something(tryparse(varTypes[i], ln[(idxRanges[i]+1):idxRanges[i+1]]), varTypes[i](0)) for i in 1:(length(varLengths))]

            # Check that the search criteria are met (molecule, isotope, wavenumber range, and min. line-strength)
            if((values[1] == mol || mol == -1) && (values[2] == iso || iso == -1)
                && (ν_min <= values[3] <= ν_max) && values[4] >= min_strength)

                # Add this row to the list of rows
                rows = append!(rows, [values])
            end
        end

        # Convert the row-major list of lists to a column-major list of lists
        cols = [[rows[i][j] for i in 1:length(rows)] for j in 1:length(varNames)]

        # Check that there is at least 1 matching record in the HITRAN file
        length(cols[1]) == 0 ? throw(HitranEmptyError()) : nothing

        # Return a HitranTable struct populated with the obtained columns
        return HitranTable(mol=cols[1], iso=cols[2], νᵢ=cols[3], Sᵢ=cols[4], Aᵢ=cols[5], γ_air=cols[6], γ_self=cols[7], E″=cols[8], n_air=cols[9], δ_air=cols[10], global_upper_quanta=cols[11], global_lower_quanta=cols[12], local_upper_quanta=cols[13], local_lower_quanta=cols[14], ierr=cols[15], iref=cols[16], line_mixing_flag=cols[17], g′=cols[18], g″=cols[19])

    end

end

"""
    make_hitran_model(hitran::HitranTable,
                      broadening::AbstractBroadeningFunction;
                      wing_cutoff::Real=40,
                      vmr::Real=0,
                      CEF::AbstractComplexErrorFunction=HumlicekWeidemann32SDErrorFunction(),
                      architecture = default_architecture)

Convenience function to make a HitranModel out of the parameters (Matches make_interpolation_model)

"""
function make_hitran_model(hitran::HitranTable,
                           broadening::AbstractBroadeningFunction;
                           wing_cutoff::Integer=40,
                           vmr::Real=0,
                           CEF::AbstractComplexErrorFunction=HumlicekWeidemann32SDErrorFunction(),
                           architecture = default_architecture)

    if architecture isa GPU && !(CEF isa HumlicekWeidemann32SDErrorFunction)
        @warn "Cross-section calculations on GPU may or may not work with this CEF (use HumlicekWeidemann32SDErrorFunction if you encounter issues)"
    end

    return HitranModel(hitran=hitran, broadening=broadening , wing_cutoff=wing_cutoff , vmr=vmr, CEF=CEF, architecture=architecture)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct VirtualMachineScaleSetVMExtensionsApi <: SwaggerApi
    client::Swagger.Client
end

"""

The operation to create or update the VMSS VM extension.
Param: resourceGroupName::String (required)
Param: vmScaleSetName::String (required)
Param: instanceId::String (required)
Param: vmExtensionName::String (required)
Param: extensionParameters::VirtualMachineExtension (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: VirtualMachineExtension
"""
function _swaggerinternal_virtualMachineScaleSetVMExtensionsCreateOrUpdate(_api::VirtualMachineScaleSetVMExtensionsApi, resourceGroupName::String, vmScaleSetName::String, instanceId::String, vmExtensionName::String, extensionParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", VirtualMachineExtension, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/virtualMachineScaleSets/{vmScaleSetName}/virtualMachines/{instanceId}/extensions/{vmExtensionName}", ["azure_auth"], extensionParameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "vmScaleSetName", vmScaleSetName)  # type String
    Swagger.set_param(_ctx.path, "instanceId", instanceId)  # type String
    Swagger.set_param(_ctx.path, "vmExtensionName", vmExtensionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualMachineScaleSetVMExtensionsCreateOrUpdate(_api::VirtualMachineScaleSetVMExtensionsApi, resourceGroupName::String, vmScaleSetName::String, instanceId::String, vmExtensionName::String, extensionParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualMachineScaleSetVMExtensionsCreateOrUpdate(_api, resourceGroupName, vmScaleSetName, instanceId, vmExtensionName, extensionParameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualMachineScaleSetVMExtensionsCreateOrUpdate(_api::VirtualMachineScaleSetVMExtensionsApi, response_stream::Channel, resourceGroupName::String, vmScaleSetName::String, instanceId::String, vmExtensionName::String, extensionParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualMachineScaleSetVMExtensionsCreateOrUpdate(_api, resourceGroupName, vmScaleSetName, instanceId, vmExtensionName, extensionParameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

The operation to delete the VMSS VM extension.
Param: resourceGroupName::String (required)
Param: vmScaleSetName::String (required)
Param: instanceId::String (required)
Param: vmExtensionName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_virtualMachineScaleSetVMExtensionsDelete(_api::VirtualMachineScaleSetVMExtensionsApi, resourceGroupName::String, vmScaleSetName::String, instanceId::String, vmExtensionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/virtualMachineScaleSets/{vmScaleSetName}/virtualMachines/{instanceId}/extensions/{vmExtensionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "vmScaleSetName", vmScaleSetName)  # type String
    Swagger.set_param(_ctx.path, "instanceId", instanceId)  # type String
    Swagger.set_param(_ctx.path, "vmExtensionName", vmExtensionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualMachineScaleSetVMExtensionsDelete(_api::VirtualMachineScaleSetVMExtensionsApi, resourceGroupName::String, vmScaleSetName::String, instanceId::String, vmExtensionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualMachineScaleSetVMExtensionsDelete(_api, resourceGroupName, vmScaleSetName, instanceId, vmExtensionName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualMachineScaleSetVMExtensionsDelete(_api::VirtualMachineScaleSetVMExtensionsApi, response_stream::Channel, resourceGroupName::String, vmScaleSetName::String, instanceId::String, vmExtensionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualMachineScaleSetVMExtensionsDelete(_api, resourceGroupName, vmScaleSetName, instanceId, vmExtensionName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

The operation to get the VMSS VM extension.
Param: resourceGroupName::String (required)
Param: vmScaleSetName::String (required)
Param: instanceId::String (required)
Param: vmExtensionName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: expand::String
Return: VirtualMachineExtension
"""
function _swaggerinternal_virtualMachineScaleSetVMExtensionsGet(_api::VirtualMachineScaleSetVMExtensionsApi, resourceGroupName::String, vmScaleSetName::String, instanceId::String, vmExtensionName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", VirtualMachineExtension, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/virtualMachineScaleSets/{vmScaleSetName}/virtualMachines/{instanceId}/extensions/{vmExtensionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "vmScaleSetName", vmScaleSetName)  # type String
    Swagger.set_param(_ctx.path, "instanceId", instanceId)  # type String
    Swagger.set_param(_ctx.path, "vmExtensionName", vmExtensionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualMachineScaleSetVMExtensionsGet(_api::VirtualMachineScaleSetVMExtensionsApi, resourceGroupName::String, vmScaleSetName::String, instanceId::String, vmExtensionName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_virtualMachineScaleSetVMExtensionsGet(_api, resourceGroupName, vmScaleSetName, instanceId, vmExtensionName, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualMachineScaleSetVMExtensionsGet(_api::VirtualMachineScaleSetVMExtensionsApi, response_stream::Channel, resourceGroupName::String, vmScaleSetName::String, instanceId::String, vmExtensionName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_virtualMachineScaleSetVMExtensionsGet(_api, resourceGroupName, vmScaleSetName, instanceId, vmExtensionName, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

The operation to get all extensions of an instance in Virtual Machine Scaleset.
Param: resourceGroupName::String (required)
Param: vmScaleSetName::String (required)
Param: instanceId::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: expand::String
Return: VirtualMachineExtensionsListResult
"""
function _swaggerinternal_virtualMachineScaleSetVMExtensionsList(_api::VirtualMachineScaleSetVMExtensionsApi, resourceGroupName::String, vmScaleSetName::String, instanceId::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", VirtualMachineExtensionsListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/virtualMachineScaleSets/{vmScaleSetName}/virtualMachines/{instanceId}/extensions", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "vmScaleSetName", vmScaleSetName)  # type String
    Swagger.set_param(_ctx.path, "instanceId", instanceId)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualMachineScaleSetVMExtensionsList(_api::VirtualMachineScaleSetVMExtensionsApi, resourceGroupName::String, vmScaleSetName::String, instanceId::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_virtualMachineScaleSetVMExtensionsList(_api, resourceGroupName, vmScaleSetName, instanceId, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualMachineScaleSetVMExtensionsList(_api::VirtualMachineScaleSetVMExtensionsApi, response_stream::Channel, resourceGroupName::String, vmScaleSetName::String, instanceId::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_virtualMachineScaleSetVMExtensionsList(_api, resourceGroupName, vmScaleSetName, instanceId, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

The operation to update the VMSS VM extension.
Param: resourceGroupName::String (required)
Param: vmScaleSetName::String (required)
Param: instanceId::String (required)
Param: vmExtensionName::String (required)
Param: extensionParameters::VirtualMachineExtensionUpdate (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: VirtualMachineExtension
"""
function _swaggerinternal_virtualMachineScaleSetVMExtensionsUpdate(_api::VirtualMachineScaleSetVMExtensionsApi, resourceGroupName::String, vmScaleSetName::String, instanceId::String, vmExtensionName::String, extensionParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PATCH", VirtualMachineExtension, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/virtualMachineScaleSets/{vmScaleSetName}/virtualMachines/{instanceId}/extensions/{vmExtensionName}", ["azure_auth"], extensionParameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "vmScaleSetName", vmScaleSetName)  # type String
    Swagger.set_param(_ctx.path, "instanceId", instanceId)  # type String
    Swagger.set_param(_ctx.path, "vmExtensionName", vmExtensionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualMachineScaleSetVMExtensionsUpdate(_api::VirtualMachineScaleSetVMExtensionsApi, resourceGroupName::String, vmScaleSetName::String, instanceId::String, vmExtensionName::String, extensionParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualMachineScaleSetVMExtensionsUpdate(_api, resourceGroupName, vmScaleSetName, instanceId, vmExtensionName, extensionParameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualMachineScaleSetVMExtensionsUpdate(_api::VirtualMachineScaleSetVMExtensionsApi, response_stream::Channel, resourceGroupName::String, vmScaleSetName::String, instanceId::String, vmExtensionName::String, extensionParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualMachineScaleSetVMExtensionsUpdate(_api, resourceGroupName, vmScaleSetName, instanceId, vmExtensionName, extensionParameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export virtualMachineScaleSetVMExtensionsCreateOrUpdate, virtualMachineScaleSetVMExtensionsDelete, virtualMachineScaleSetVMExtensionsGet, virtualMachineScaleSetVMExtensionsList, virtualMachineScaleSetVMExtensionsUpdate


text:
#=
Power-law slope plot for the paper.
=#
using DrWatson
@quickactivate "VelocityMTD"
using MusicManipulations, PyPlot, Statistics
include(srcdir("style.jl"))

datasets = ["tapping", "uwe", "pgmusic", "playalongs"]
datanames = ["tapping",  "Uwe Meile (piano)", "PG music (piano)", "drums playalongs"]

fig = figure(figsize = (2figx, 3figx/4))
axβ = subplot(311)
axδ = subplot(312)
axb = subplot(313)
axs = (axβ, axδ, axb)

alldata = load(datadir("psd", "alldatasets.bson"))

# %% plot
pcount = 0
for ax in axs; ax.clear(); end

for i in 1:length(datasets)
    dataset = datasets[i]
    data = alldata[dataset]

    βs = data[:βs]
    δs = data[:δs]
    bs = data[:bs]

    L = length(bs)
    l = pcount + 1
    r = l + L - 1

    for (j, (ax, q)) in enumerate(zip(axs, (βs, δs, bs)))
        ax.plot(l:r, q, color = "C$(i)", ls = "None", marker = "o")
        j != 3 && ax.plot([l, r], fill(mean(q), 2), color = "C$(i)", lw = 2, alpha = 0.75)
    end

    axb.text((r+l)/2, 3.5, datanames[i], ha="center", va="center", color = "C$i")

    global pcount += L+1
    # break
end

for ax in axs; ax.set_xlim(0, 54); end
for ax in axs[1:2]; ax.set_xticklabels([]); end

axβ.set_ylabel("\$\\beta\$")
axβ.set_yticks(-1.2:0.4:0.2)
axβ.set_ylim(-1.2, 0.2)
axβ.axhline(0, zorder = - 99)
axβ.invert_yaxis()
axβ.yaxis.grid(true)
# axβ.axhline(0, color = "black", ls = "dashed", lw = .5)
# axβ.axhline(-1, color = "black", ls = "dashed", lw = .5)

axδ.set_ylabel("\$\\delta\$")
axδ.set_yticks(0:0.4:1.2)
axδ.axhline(0, zorder = - 99)

axδ.set_ylim(-0.2, 1.2)
axδ.yaxis.grid(true)
# axδ.axhline(0, color = "black", ls = "dashed", lw = .5)
# axδ.axhline(1, color = "black", ls = "dashed", lw = .5)

axb.set_ylabel("\$b_c\$")
axb.set_yticks(1:4)
axb.set_xlabel("recording #")

add_identifiers!(fig; loc = (0.995, 0.975))
fig.align_ylabels()
fig.tight_layout()
fig.subplots_adjust(hspace = 0.2)
# fig.savefig(plotsdir("psd", "allpsd"))

max_stars_repo_name: Datseris/VelocityMTD

max_issues_repo_name: Datseris/VelocityMTD

max_forks_repo_name: Datseris/VelocityMTD


text:
# using Optim
# using Elliptic
# using Struve
# using SpecialFunctions

function CalcFlux(Ae,Le,μr,N,V,ω)
    μ0 = 4*pi*1e-7
    MagnetizingInductance = N^2*μr*μ0*Ae /Le
    Impedance = ω*MagnetizingInductance
    I = V/Impedance
    Flux = N*I*μr*μ0*Ae /Le

    B = Flux/Ae
    return B
end

function CalcPermeability(AL,Ae,Le)
    #This function takes in the nominal inductance (inductance at one turn),effective length and area
    #And calculates the permeabily
    # L = N^2 μ_r*μ_o A_e /L_e ⟶ N=1, A_L = μA_e/L_e
    μ0 = 4*pi*1e-7
    return AL*Le/Ae/μ0
end

function BMax(V,f,N,A;k=(2*π/√2))
    return V ./ (k.*N.*A.*f)
end

function MaxTurns(CoreID,WireD)
    InnerWindingCircum = pi*(CoreID-WireD)
    return InnerWindingCircum/WireD
end

function RunMagnetizingFlux(WireD,CoreID,Ae,Le,AL,NRatio,VPk,ω)
    μ0 = 4*pi*1e-7
    μr = CalcPermeability(AL,Ae,Le)
    NTot  = MaxTurns(CoreID,WireD)
    NPrimary = NTot* (1/(1+NRatio))
    BCore = CalcFlux(Ae,Le,μr,NPrimary,VPk,ω)
    println("
    μr              = $(round(μr;sigdigits=3))
    Primary Turns   = $(round(NPrimary;sigdigits=3))
    Total Turns     = $(round(NTot;sigdigits=3))
    BCore(mT)       = $(round(BCore*1e3;sigdigits=3))
    ")
end

"""
Calculate Magnetizing Inductance, units in meters
"""
function MagInductance(Aₑ,Lₑ,N,μᵣ)
    μ₀ = 4*pi*1e-7
    L = N^2* μᵣ*μ₀ * Aₑ /Lₑ
    return L
end

"""
Calculate Magnetizing Inductance, units in meters
"""
function MagInductance_Gap(Aₑ,Lₑ,N,μᵣ,LGap)
    μ₀ = 4*pi*1e-7
    ℛgap = LGap/(μ₀*Aₑ)#\scrR
    ℛcore = Lₑ/(μ₀*μᵣ*Aₑ)#\scrR

    L = N^2 / (ℛcore+ℛgap)

    return L
end

function MagnetizingCurrent(L,f)
    return 2*π*f*L
end
""""
determine the resistance of a winding on a rectangular core with N layers of copper film winding

N is the number of layers
t the the thickness of copper
L is the length of the core
w is the width of the core (Inner dimension)
H is the height of the core
ρ is the resistivity of copper
"""
function FilmCoreResist_Rect(N,t,L,w,H,ρ)
    Length = N* 2*( (w+N/2*t) + (L+N/2*t))
    println(Length)
    R = ρ*Length / (H*t)
    return R
end

""""
determine the resistance of a winding on a rectangular core with Nₗ layers and Nₜ turns of copper circular wire

Nₗ is the number of layers
Nₜ is the number of turns
D the the diamter of copper wire
L is the length of the core
w is the width of the core (Inner dimension)
H is the height of the core
ρ is the resistivity of copper
"""
function WireWoundCoreResist_Rect(Nₗ,Nₜ,D,L,w,H,ρ)
    Length = Nₜ* 2*( (w+Nₗ/2*D) + (L+Nₗ/2*D))
    println(Length)
    R = ρ*Length / (π/4*D^2)
    return R
end

"""
The minimum allowable window area given some current density J₀, Maximum current, and turns (N)
"""
function MinWindowArea(Imax, J₀,N)
    return Imax*N/J₀
end

function MinWindowArea(Vmax,Imax, BMax, Acore,f, J₀)
    return Imax*MinTurns(Vmax, BMax, Acore,f) /J₀
end

function MinTurns(Vmax, BMax, Acore,f)
    return Vmax*√2 /(2*π*f* BMax * Acore)
end

function MeasuredInductance2Params(ZₚSecₒ,ZₚSecₛ,ZₛPriₒ,ZₛPriₛ)

    ZₚEqiv_SecLoad(Zₗ1,Zₘ1,Zₗ2,N,RTerm) = Zₗ1+ Par(Zₘ1,(1 / N^2 * (Zₗ2+RTerm))) #Primary referred impedance with some secondary termination (RTerm)
    ZₛEqiv_PriLoad(Zₗ1,Zₘ1,Zₗ2,N,RTerm) = Zₗ2+ N^2 * Par(Zₘ1,(Zₗ1+RTerm)) #Primary referred impedance with some secondary termination (RTerm)

    Cost(Zₗ1,Zₘ1,Zₗ2,N) = abs((ZₚSecₒ - ZₚEqiv_SecLoad(Zₗ1,Zₘ1,Zₗ2,N,1e9)) + (ZₚSecₛ - ZₚEqiv_SecLoad(Zₗ1,Zₘ1,Zₗ2,N,0)) +
                         (ZₛPriₒ - ZₛEqiv_PriLoad(Zₗ1,Zₘ1,Zₗ2,N,1e9)) + (ZₛPriₛ - ZₛEqiv_PriLoad(Zₗ1,Zₘ1,Zₗ2,N,0)))

    CostVecIn(Vec) = Cost(Vec[1],Vec[2],Vec[3],Vec[4])
    optimize(CostVecIn, [10.0 , 40.0, 10.0, 2.0])

end

# """
# Function to determine the leakage inductnace of a toroidal transformer
# Ref. `Calculation of self and mutual impedances between sections of tranformer windings`, Wilcox et al. IEE Proceedings, 1989

# """
# function CalcLeakageInductance_Toroid()

# end
# """
# Function to determine the leakage inductnace of a toroidal transformer
# Ref. `Calculation of self and mutual impedances between sections of tranformer windings`, Wilcox et al. IEE Proceedings, 1989

# """
# function CalcLₖₘ(Nₖ,Nₘ,r,a,z)
#     μ₀ = 4*π*1e-7
#     κ =κFunk(a,r,z)
#     Lₖₘ = μ₀*Nₖ*Nₘ*√(r*a)*2/κ*(  (1 - κ^2/2) * Elliptic.K(κ) - Elliptic.E(κ) )

# end

# function Z₁₍ₖₘ₎(Nₖ,Nₘ,b,μᵣ,ω,λ,ρ)

#     μ₁ = 4*π*1e-7 #for consistency with eqn. using 1 subscript
#     μₓ =4*π*1e-7*μᵣ
#     m = √(im*ω*μₓ/ρ)
#     Z = 1*im*ω*Nₖ*Nₘ *π*b^2/λ * (2*μₓ*besseli(1,m*b) / (m*b*besseli(0,m*b)) - μ₁)

# end

# """
# Equation 12,Ref. `Calculation of self and mutual impedances between sections of tranformer windings`, Wilcox et al. IEE Proceedings, 1989
# """
# function BigP₁(x,y,β)
#     (p₁(x) - p₁(y)) ./ (β.^2)

# end

# function p₁(α)
#     π*α/2 * (bessely(0,α)*Struve.L(0, α)+bessely(1,α)*Struve.L(1, α))
# end

# """
# Equation 10,Ref. `Calculation of self and mutual impedances between sections of tranformer windings`, Wilcox et al. IEE Proceedings, 1989
# """
# function κFunk(a,r,z)
#     √(4*a*r ./ (z^2 + (a+r)^2))
# end


text:
function read_yahoo_for_sixty_months(stock::String, fm::Int, fd::Int, fy::Int, tm::Int, td::Int, ty::Int, period::String)

# take care of yahoo's 0 indexing for month
  fm-=1
  tm-=1

  ydata     = readlines(`curl -s "http://ichart.finance.yahoo.com/table.csv?s=$stock&a=$fm&b=$fd&c=$fy&d=$tm&e=$td&f=$ty&g=$period"`)
  numstring = ydata[2:end]
  sa  = split(numstring[1], ',')'

  for i in 2:length(numstring)
    sa  = [sa ; split(numstring[i], ',')']
  end

  time_conversion = map(x -> parse("yyyy-MM-dd", x), convert(Array{UTF16String}, sa[:,1]))

  df = DataFrame(quote
     Date  = $time_conversion
     Open  = float($sa[:,2])
     High  = float($sa[:,3])
     Low   = float($sa[:,4])
     Close = float($sa[:,5])
     Vol   =   int($sa[:,6])
     Adj   = float($sa[:,7])
     end)

  flipud(df)

end

# default to last three years, monthly data

read_yahoo_for_sixty_months(stock::String) = read_yahoo_for_sixty_months(stock::String, month(now()), day(now()), year(now())-5, month(now()),  day(now()), year(now()), "m")

function stock_returns{T <: Union(ASCIIString, UTF8String)}(list::Array{T})

  str1 = list[1]
  str2 = list[2]
  str3 = list[3]
  str4 = list[4]
  str5 = list[5]

  s1 = read_yahoo_for_sixty_months(list[1])
  s2 = read_yahoo_for_sixty_months(list[2])
  s3 = read_yahoo_for_sixty_months(list[3])
  s4 = read_yahoo_for_sixty_months(list[4])
  s5 = read_yahoo_for_sixty_months(list[5])

  r1 = log_return(s1["Adj"])
  r2 = log_return(s2["Adj"])
  r3 = log_return(s3["Adj"])
  r4 = log_return(s4["Adj"])
  r5 = log_return(s5["Adj"])

  df = DataFrame(quote
                 $str1 = $r1
                 $str2 = $r2
                 $str3 = $r3
                 $str4 = $r4
                 $str5 = $r5
                 end)
  df[2:end, :]
end

function read_monthlies(filename::String)

df  = read_table(filename)
time_conversion =  Calendar.parse("MMM-yy", df[:,1])

within!(df, quote
        Date = $(time_conversion)
        end)

#flipud(df)
end

function read_yahoo(stock::String, fm::Int, fd::Int, fy::Int, tm::Int, td::Int, ty::Int, period::String)

# take care of yahoo's 0 indexing for month
  fm-=1
  tm-=1

  ydata = readlines(`curl -s "http://ichart.finance.yahoo.com/table.csv?s=$stock&a=$fm&b=$fd&c=$fy&d=$tm&e=$td&f=$ty&g=$period"`)
  numstring = ydata[2:end]
  sa  = split(numstring[1], ',')'

  for i in 2:length(numstring)
    sa  = [sa ; split(numstring[i], ',')']
  end

  time_conversion = map(x -> parse("yyyy-MM-dd", x), convert(Array{UTF16String}, sa[:,1]))

  df = DataFrame(quote
     Date  = $time_conversion
     Open  = float($sa[:,2])
     High  = float($sa[:,3])
     Low   = float($sa[:,4])
     Close = float($sa[:,5])
     Vol   =   int($sa[:,6])
     Adj   = float($sa[:,7])
     end)

  flipud(df)
end

# default to last three years, daily data
read_yahoo(stock::String) = read_yahoo(stock::String, month(now()), day(now()), year(now())-3, month(now()),  day(now()), year(now()), "d")


text:
struct CKSAFTParam <: EoSParam
    segment::SingleParam{Float64}
    sigma::PairParam{Float64}
    epsilon::PairParam{Float64}
    c::SingleParam{Float64}
    epsilon_assoc::AssocParam{Float64}
    bondvol::AssocParam{Float64}
end

abstract type CKSAFTModel <: SAFTModel end
@newmodel CKSAFT CKSAFTModel CKSAFTParam

export CKSAFT
function CKSAFT(components; idealmodel=BasicIdeal, userlocations=String[], ideal_userlocations=String[], verbose=false)
    params,sites = getparams(components, ["SAFT/CKSAFT","properties/molarmass.csv"]; userlocations=userlocations, verbose=verbose)
    segment = params["m"]
    c = params["c"]
    k = params["k"]
    sigma = params["vol"]
    sigma.values .*= 6*0.74048/N_A/1e6/π
    sigma.values .^= 1/3
    sigma = sigma_LorentzBerthelot(sigma)
    epsilon = epsilon_LorentzBerthelot(params["epsilon"], k)
    epsilon_assoc = params["epsilon_assoc"]
    bondvol = params["bond_vol"]
    packagedparams = CKSAFTParam(segment, sigma, epsilon,c, epsilon_assoc, bondvol)
    references = ["TODO CKSAFT", "TODO CKSAFT"]

    model = CKSAFT(packagedparams, sites, idealmodel; ideal_userlocations=ideal_userlocations, references=references, verbose=verbose)
    return model
end

function a_res(model::CKSAFTModel, V, T, z)
    return @f(a_seg) + @f(a_chain) + @f(a_assoc)
end

function a_seg(model::CKSAFTModel, V, T, z)
    Σz = ∑(z)
    comps = @comps
    m = model.params.segment.values
    m̄ = sum(0.5*(m[i]+m[j])*z[i]*z[j] for i ∈ comps for j ∈ comps)/(Σz*Σz)
    return m̄*(@f(a_hs) + @f(a_disp))
end

function a_hs(model::CKSAFTModel, V, T, z)
    ζ0 = @f(ζ,0)
    ζ1 = @f(ζ,1)
    ζ2 = @f(ζ,2)
    ζ3 = @f(ζ,3)
    return 1/ζ0 * (3ζ1*ζ2/(1-ζ3) + ζ2^3/(ζ3*(1-ζ3)^2) + (ζ2^3/ζ3^2-ζ0)*log(1-ζ3))
end

function a_disp(model::CKSAFTModel, V, T, z)
    ϵ̄ = @f(ū)
    η = @f(ζ,3)
    τ = 0.74048
    D1 = CKSAFT_consts.D1
    D2 = CKSAFT_consts.D2
    D3 = CKSAFT_consts.D3
    D4 = CKSAFT_consts.D4
    eT = (ϵ̄/T)
    A1 = ∑(D1[j]*eT*(η/τ)^j for j ∈ 1:6)
    A2 = ∑(D2[j]*eT^2*(η/τ)^j for j ∈ 1:9)
    A3 = ∑(D3[j]*eT^3*(η/τ)^j for j ∈ 1:5)
    A4 = ∑(D4[j]*eT^4*(η/τ)^j for j ∈ 1:4)
    return A1+A2+A3+A4
end

function d(model::CKSAFTModel, V, T,z, i)
    ϵ = model.params.epsilon.diagvalues[i]
    σ = model.params.sigma.diagvalues[i]
    return σ * (1 - 0.12exp(-3ϵ/T))
end

function d(model::CKSAFTModel, V, T, z, i,j)
    #return (d(model::CKSAFTModel, z, v, T, i)+d(model::CKSAFTModel, z, v, T, j))/2
    return (@f(d,i) + @f(d,j))*0.5
end

function u(model::CKSAFTModel, V, T, z, i,j)
    ϵ0 = model.params.epsilon.values[i,j]
    c1 = model.params.c.values[i]
    c2 = model.params.c.values[j]
    return ϵ0*sqrt((1+c1/T)*(1+c2/T))
end

function ū(model::CKSAFTModel, V, T, z)
    Σz = ∑(z)
    x = z.*(1/Σz)
    it = @comps #only to stop vscode from recognizing this as an error
    m = model.params.segment.values
    num = ∑(
            x[i]*x[j]*m[i]*m[j]*@f(u,i,j)*@f(d,i,i)^3*@f(d,j,j)^3
            for i ∈ it for j ∈ it)

    denom = ∑(
            x[i]*m[i]*@f(d,i,i)^3
            for i ∈ it)
    return num/denom^2
end

function ζ(model::CKSAFTModel, V, T, z, n)
    ∑z = ∑(z)
    m = model.params.segment.values
    return N_A*π/6/V * ∑(z[i]*m[i]*@f(d,i)^n for i ∈ @comps)
end

function a_chain(model::CKSAFTModel, V, T, z)
    m = model.params.segment.values
    return ∑(z[i]*(1-m[i])*log(@f(g_hsij,i,i)) for i ∈ @comps)/∑(z)
end

function g_hsij(model::CKSAFTModel, V, T, z, i, j)
    di = @f(d,i)
    dj = @f(d,j)
    ζ2 = @f(ζ,2)
    ζ3 = @f(ζ,3)
    return 1/(1-ζ3) + di*dj/(di+dj)*3ζ2/(1-ζ3)^2 + (di*dj/(di+dj))^2*2ζ2^2/(1-ζ3)^3
end

function Δ(model::CKSAFTModel, V, T, z, i, j, a, b)
    ϵ_associjab = model.params.epsilon_assoc.values[i,j][a,b]
    κijab = model.params.bondvol.values[i,j][a,b]
    σij = model.params.sigma.values[i,j]
    gij = @f(g_hsij,i,j)
    return gij*σij^3*(exp(ϵ_associjab/T)-1)*κijab
end

const CKSAFT_consts =(
    D1 = [-8.8043,4.164627,-48.203555,140.4362,-195.23339,113.515],
    D2 = [2.9396,-6.0865383,40.137956,-76.230797,-133.70055,860.25349,-1535.3224,1221.4261,-409.10539],
    D3 = [-2.8225,4.7600148,11.257177,-66.382743,69.248785],
    D4 = [0.34,-3.1875014,12.231796,-12.110681],
    D =
    [0.9105631445 -0.3084016918 -0.0906148351;
    0.6361281449 0.1860531159 0.4527842806;
    2.6861347891 -2.5030047259 0.5962700728;
    -26.547362491 21.419793629 -1.7241829131;
    97.759208784 -65.255885330 -4.1302112531;
    -159.59154087 83.318680481 13.776631870;
    91.297774084 -33.746922930 -8.6728470368],
)


text:
using Tutorials
using Test

function get_title(files,filename)
  for (title,fn) in files
    if fn == filename
      return title
    end
  end
  error("File $filename not found!")
end

if (length(ARGS) != 0)
  files = [get_title(Tutorials.files,filename)=>filename for filename in ARGS]
else
  files = Tutorials.files
end

for (title,filename) in files
    # Create temporal modules to isolate and protect test scopes
    tmpdir=mktempdir(;cleanup=true)
    filename_wo_extension=split(filename,".")[1]
    tmpmod = filename_wo_extension
    tmpfile = joinpath(tmpdir,tmpmod)
    isfile(tmpfile) && error("File $tmpfile already exists!")
    testpath = joinpath(@__DIR__,"../src", filename)
    open(tmpfile,"w") do f
      println(f, "# This file is automatically generated")
      println(f, "# Do not edit")
      println(f)
      println(f, "module $tmpmod include(\"$testpath\") end")
    end
    @time @testset "$title" begin include(tmpfile) end
end

# module fsi_tutorial
# using Test
# @time @testset "fsi_tutorial" begin include("../src/fsi_tutorial.jl") end
# end # module


text:
#
# Correctness Tests
#

using Calculus
using Test
using LinearAlgebra

tests = ["finite_difference",
         "derivative",
         "check_derivative",
         "symbolic",
         "deparse"]

for t in tests
    @testset "$t" begin
        include("$(t).jl")
    end
end


text:
using MPI
using OrderedCollections
using StaticArrays
using Statistics
using Dierckx
using Test
using Pkg.Artifacts
using DelimitedFiles

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using ClimateMachine
using ClimateMachine.Land
using ClimateMachine.Land.SoilWaterParameterizations
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.DGMethods: BalanceLaw, LocalGeometry
using ClimateMachine.MPIStateArrays
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates
using ClimateMachine.SingleStackUtils
using ClimateMachine.BalanceLaws:
    BalanceLaw, Prognostic, Auxiliary, Gradient, GradientFlux, vars_state

@testset "hydrostatic test 1" begin
    ClimateMachine.init()
    FT = Float64

    function init_soil_water!(land, state, aux, localgeo, time)
        myfloat = eltype(aux)
        state.soil.water.ϑ_l = myfloat(land.soil.water.initialϑ_l(aux))
        state.soil.water.θ_i = myfloat(land.soil.water.initialθ_i(aux))
    end

    soil_heat_model = PrescribedTemperatureModel()

    Ksat = (aux) -> eltype(aux)(0.0443 / (3600 * 100))
    S_s = (aux) -> eltype(aux)((1e-3) * exp(-0.2 * aux.z))
    vgn = FT(2)
    wpf = WaterParamFunctions(FT; Ksat = Ksat, S_s = S_s)

    soil_param_functions = SoilParamFunctions(FT; porosity = 0.495, water = wpf)
    bottom_flux = (aux, t) -> eltype(aux)(0.0)
    surface_flux = bottom_flux
    ϑ_l0 = (aux) -> eltype(aux)(0.494)
    bc = LandDomainBC(
        bottom_bc = LandComponentBC(soil_water = Neumann(bottom_flux)),
        surface_bc = LandComponentBC(soil_water = Neumann(surface_flux)),
    )
    soil_water_model = SoilWaterModel(
        FT;
        moisture_factor = MoistureDependent{FT}(),
        hydraulics = vanGenuchten(FT; n = vgn),
        initialϑ_l = ϑ_l0,
    )

    m_soil = SoilModel(soil_param_functions, soil_water_model, soil_heat_model)
    sources = ()
    m = LandModel(
        param_set,
        m_soil;
        boundary_conditions = bc,
        source = sources,
        init_state_prognostic = init_soil_water!,
    )

    N_poly = 2
    nelem_vert = 20

    # Specify the domain boundaries
    zmax = FT(0)
    zmin = FT(-10)

    driver_config = ClimateMachine.SingleStackConfiguration(
        "LandModel",
        N_poly,
        nelem_vert,
        zmax,
        param_set,
        m;
        zmin = zmin,
        numerical_flux_first_order = CentralNumericalFluxFirstOrder(),
    )

    t0 = FT(0)
    timeend = FT(60 * 60 * 24 * 200)

    dt = FT(500)
    n_outputs = 3
    every_x_simulation_time = ceil(Int, timeend / n_outputs)
    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        ode_dt = dt,
    )
    aux = solver_config.dg.state_auxiliary
    state_types = (Prognostic(), Auxiliary())
    dons_arr =
        Dict[dict_of_nodal_states(solver_config, state_types; interp = true)]
    time_data = FT[0]
    callback = GenericCallbacks.EveryXSimulationTime(every_x_simulation_time) do
        dons = dict_of_nodal_states(solver_config, state_types; interp = true)
        push!(dons_arr, dons)
        push!(time_data, gettime(solver_config.solver))
        nothing
    end
    ClimateMachine.invoke!(solver_config; user_callbacks = (callback,))
    z = dons_arr[1]["z"]
    interface_z = -1.0395
    function hydrostatic_profile(z, zm, porosity, n, α)
        myf = eltype(z)
        m = FT(1 - 1 / n)
        S = FT((FT(1) + (α * (z - zm))^n)^(-m))
        return FT(S * porosity)
    end
    function soln(z, interface, porosity, n, α, δ, S_s)
        if z < interface
            return porosity + S_s * (interface - z) * exp(-δ * z)
        else
            return hydrostatic_profile(z, interface, porosity, n, α)
        end
    end

    MSE = mean(
        (
            soln.(z, interface_z, 0.495, vgn, 2.6, 0.2, 1e-3) .-
            dons_arr[4]["soil.water.ϑ_l"]
        ) .^ 2.0,
    )
    @test MSE < 1e-4
end

@testset "hydrostatic test 2" begin
    ClimateMachine.init()
    FT = Float64

    function init_soil_water!(land, state, aux, localgeo, time)
        myfloat = eltype(aux)
        state.soil.water.ϑ_l = myfloat(land.soil.water.initialϑ_l(aux))
        state.soil.water.θ_i = myfloat(land.soil.water.initialθ_i(aux))
    end

    soil_heat_model = PrescribedTemperatureModel()

    Ksat = (4.42 / 3600 / 100)
    S_s = 1e-3
    wpf = WaterParamFunctions(FT; Ksat = Ksat, S_s = S_s)

    soil_param_functions = SoilParamFunctions(FT, porosity = 0.6, water = wpf)
    bottom_flux = (aux, t) -> eltype(aux)(0.0)
    surface_flux = bottom_flux
    bc = LandDomainBC(
        bottom_bc = LandComponentBC(soil_water = Neumann(bottom_flux)),
        surface_bc = LandComponentBC(soil_water = Neumann(surface_flux)),
    )
    sigmoid(x, offset, width) =
        typeof(x)(exp((x - offset) / width) / (1 + exp((x - offset) / width)))
    function soln(z::f, interface::f, porosity::f) where {f}
        function hydrostatic_profile(
            z::f,
            interface::f,
            porosity::f,
            n::f,
            α::f,
            m::f,
        )
            ψ_interface = f(-1)
            ψ = -(z - interface) + ψ_interface
            S = (f(1) + (-α * ψ)^n)^(-m)
            return S * porosity
        end
        if z < interface
            return hydrostatic_profile(
                z,
                interface,
                porosity,
                f(1.31),
                f(1.9),
                f(1) - f(1) / f(1.31),
            )
        else
            return hydrostatic_profile(
                z,
                interface,
                porosity,
                f(1.89),
                f(7.5),
                f(1) - f(1) / f(1.89),
            )
        end
    end
    ϑ_l0 = (aux) -> soln(aux.z, -1.0, 0.6)

    vgα(aux) = aux.z < -1.0 ? 1.9 : 7.5
    vgn(aux) = aux.z < -1.0 ? 1.31 : 1.89

    soil_water_model = SoilWaterModel(
        FT;
        moisture_factor = MoistureDependent{FT}(),
        hydraulics = vanGenuchten(FT; n = vgn, α = vgα),
        initialϑ_l = ϑ_l0,
    )

    m_soil = SoilModel(soil_param_functions, soil_water_model, soil_heat_model)
    sources = ()
    m = LandModel(
        param_set,
        m_soil;
        boundary_conditions = bc,
        source = sources,
        init_state_prognostic = init_soil_water!,
    )

    N_poly = 1
    nelem_vert = 80

    # Specify the domain boundaries
    zmax = FT(0)
    zmin = FT(-2)

    driver_config = ClimateMachine.SingleStackConfiguration(
        "LandModel",
        N_poly,
        nelem_vert,
        zmax,
        param_set,
        m;
        zmin = zmin,
        numerical_flux_first_order = CentralNumericalFluxFirstOrder(),
    )

    t0 = FT(0)
    timeend = FT(60 * 60 * 12)

    dt = FT(5)
    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        ode_dt = dt,
    )

    ClimateMachine.invoke!(solver_config)

    state_types = (Prognostic(), Auxiliary())
    dons = dict_of_nodal_states(solver_config, state_types; interp = true)
    z = dons["z"]

    RMSE =
        mean(
            (soln.(z, Ref(-1.0), Ref(0.6)) .- dons["soil.water.ϑ_l"]) .^ 2.0,
        )^0.5
    @test RMSE < 2.0 * eps(FT)
end


text:
synth_imgs = [
    "Float64_0" => rand(127, 257),
    "Float64_1" => rand(127, 257, 1),
    "Float64_2" => rand(127, 257, 2),
    "Float64_3" => rand(127, 257, 3),
    "Float64_4" => rand(127, 257, 4),
    "Bool_0" => rand(Bool, 127, 257),
    "Bool_1" => rand(Bool, 127, 257, 1),
    "Bool_2" => rand(Bool, 127, 257, 2),
    "Bool_3" => rand(Bool, 127, 257, 3),
    "Bool_4" => rand(Bool, 127, 257, 4),
    "UInt8_0" => rand(UInt8, 127, 257),
    "UInt8_1" => rand(UInt8, 127, 257, 1),
    "UInt8_2" => rand(UInt8, 127, 257, 2),
    "UInt8_3" => rand(UInt8, 127, 257, 3),
    "UInt8_4" => rand(UInt8, 127, 257, 4),
    "UInt16_0" => rand(UInt16, 127, 257),
    "UInt16_1" => rand(UInt16, 127, 257, 1),
    "UInt16_2" => rand(UInt16, 127, 257, 2),
    "UInt16_3" => rand(UInt16, 127, 257, 3),
    "UInt16_4" => rand(UInt16, 127, 257, 4),
    "N0f8_0" => rand(N0f8, 127, 257),
    "N0f8_1" => rand(N0f8, 127, 257, 1),
    "N0f8_2" => rand(N0f8, 127, 257, 2),
    "N0f8_3" => rand(N0f8, 127, 257, 3),
    "N0f8_4" => rand(N0f8, 127, 257, 4),
    "N0f16_0" => rand(N0f16, 127, 257),
    "N0f16_1" => rand(N0f16, 127, 257, 1),
    "N0f16_2" => rand(N0f16, 127, 257, 2),
    "N0f16_3" => rand(N0f16, 127, 257, 3),
    "N0f16_4" => rand(N0f16, 127, 257, 4),
    "Gray" => rand(Gray, 127, 257),
    "GrayA" => rand(GrayA, 127, 257),
    "RGB" => rand(RGB, 127, 257),
    "RGBA" => rand(RGBA, 127, 257),
    "ARGB32" => reinterpret(ARGB32, rand(UInt32, 127, 257)),
    "Gray-Bool" => rand(Gray{Bool}, 127, 257),
    "Gray-N0f8" => rand(Gray{N0f8}, 127, 257),
    "GrayA-N0f8" => rand(GrayA{N0f8}, 127, 257),
    "RGB-N0f8" => rand(RGB{N0f8}, 127, 257),
    "RGBA-N0f8" => rand(RGBA{N0f8}, 127, 257),
    "Gray-N0f16" => rand(Gray{N0f16}, 127, 257),
    "GrayA-N0f16" => rand(GrayA{N0f16}, 127, 257),
    "RGB-N0f16" => rand(RGB{N0f16}, 127, 257),
    "RGBA-N0f16" => rand(RGBA{N0f16}, 127, 257),
    "Gray-N0f8" => rand(Gray{N0f8}, 127, 257),
    "GrayA-N0f8" => rand(GrayA{N0f8}, 127, 257),
    "Gray-Float64" => rand(Gray{Float64}, 127, 257),
    "GrayA-Float64" => rand(GrayA{Float64}, 127, 257),
    "RGB-Float64" => rand(RGB{Float64}, 127, 257),
    "RGBA-Float64" => rand(RGBA{Float64}, 127, 257),
]

edge_case_imgs = [
    ("BitArray_0", x -> _convert(Gray, N7f1, x), randn(127, 257) .> 0),
    ("BitArray_1", x -> _convert(Gray, N7f1, x), randn(127, 257, 1) .> 0),
    ("BitArray_2", x -> _convert(GrayA, N7f1, x), randn(127, 257, 2) .> 0),
    ("BitArray_3", x -> _convert(RGB, N7f1, x), randn(127, 257, 3) .> 0),
    ("BitArray_4", x -> _convert(RGBA, N7f1, x), randn(127, 257, 4) .> 0),
    ("N4f12_0", x -> _convert(Gray, N0f16, x), rand(N4f12, 127, 257)),
    ("N4f12_1", x -> _convert(Gray, N0f16, x), rand(N4f12, 127, 257, 1)),
    ("N4f12_2", x -> _convert(GrayA, N0f16, x), rand(N4f12, 127, 257, 2)),
    ("N4f12_3", x -> _convert(RGB, N0f16, x), rand(N4f12, 127, 257, 3)),
    ("N4f12_4", x -> _convert(RGBA, N0f16, x), rand(N4f12, 127, 257, 4)),
    ("BGR-N0f8", identity, rand(BGR{N0f8}, 127, 257)),
    ("BGR-N0f16", identity, rand(BGR{N0f16}, 127, 257)),
    ("BGR-Float64", identity, BGR{Float64}.(rand(BGR{N0f8}, 127, 257))),
    ("BGRA-N0f8", identity, rand(BGRA{N0f8}, 127, 257)),
    ("BGRA-N0f16", identity, rand(BGRA{N0f16}, 127, 257)),
    ("BGRA-Float64", identity, BGRA{Float64}.(rand(BGRA{N0f8}, 127, 257))),
    ("ABGR-N0f8", identity, rand(ABGR{N0f8}, 127, 257)),
    ("ABGR-N0f16", identity, rand(ABGR{N0f16}, 127, 257)),
    ("ABGR-Float64", identity, ABGR{Float64}.(rand(ABGR{N0f8}, 127, 257))),
    ("ARGB-N0f8", identity, rand(ARGB{N0f8}, 127, 257)),
    ("ARGB-N0f16", identity, rand(ARGB{N0f16}, 127, 257)),
    ("ARGB-Float64", identity, ARGB{Float64}.(rand(ARGB{N0f8}, 127, 257))),
]


text:
#
# Concatenation.
#

function MergeQuery(qs::Vector{Query})
    dom = obound(Domain[domain(output(q)) for q in qs])
    optional = all(isoptional(q) for q in qs)
    q = RecordQuery(qs)
    return q >> Query(MergeSig(), Input(domain(output(q))), Output(dom) |> setoptional(optional) |> setplural())
end

MergeQuery(q1::Query, qrest...) =
    MergeQuery([q1], qrest...)

MergeQuery(qs1::Vector{Query}, q2::Query, qrest...) =
    MergeQuery([qs1..., q2], qrest...)

MergeQuery(qs1::Vector{Query}, qs2::Vector{Query}, qrest...) =
    MergeQuery([qs1..., qs2...], qrest...)

immutable MergeSig <: AbstractPrimitive
end

ev(::MergeSig, ::Vector{Query}, ::Input, oty::Output, iflow::InputFlow) =
    OutputFlow(oty, merge_impl(datatype(domain(oty)), length(iflow), columns(values(iflow))))

function merge_impl(T::DataType, len::Int, cols::Vector{Column})
    optional = true
    size = 0
    for col in cols
        size += length(col.vals)
        if !isoptional(col)
            optional = false
        end
    end
    offs = Vector{Int}(len+1)
    offs[1] = 1
    vals = Vector{T}(size)
    n = 1
    for k = 1:len
        for col in cols
            cr = cursor(col, k)
            copy!(vals, n, cr)
            n += length(cr)
        end
        offs[k+1] = n
    end
    return Column{optional,true}(offs, vals)
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
# using Revise, Plots
using PseudoArcLengthContinuation, LinearAlgebra, Setfield, SparseArrays
const PALC = PseudoArcLengthContinuation

k = 2
N = 10
F = (x, p) -> p[1] .* x .+ x.^(k+1)/(k+1) .+ 0.01
Jac_m = (x, p) -> diagm(0 => p[1]  .+ x.^k)

normInf = x -> norm(x, Inf)

opts = PALC.ContinuationPar(dsmax = 0.051, dsmin = 1e-3, ds=0.001, maxSteps = 140, pMin = -3., saveSolEveryNsteps = 0, newtonOptions = NewtonPar(tol = 1e-8, verbose = false), computeEigenValues = false, saveEigenvectors = false)
x0 = 0.01 * ones(N)

opts = @set opts.doArcLengthScaling = true
br0, sol, _ = @time PALC.continuation(F,Jac_m,x0, -1.5, (@lens _),opts,verbosity=0) #(15.28 k allocations: 1.055 MiB)

# test with callbacks
br0, sol, _ = @time PALC.continuation(F,Jac_m,x0, -1.5, (@lens _), (@set opts.maxSteps = 3), verbosity=2, callbackN = (x, f, J, res, iteration, itlinear, optionsN; kwargs...)->(@show x;true))

###### Used to check type stability of the methods
# using RecursiveArrayTools
# iter = PALC.PALCIterable(F,Jac_m,x0,-1.5, (@lens _), opts,verbosity=0)
# state = iterate(iter)[1]
# contRes = PALC.initContRes(iter, state)
# @time continuation!(iter, state, contRes)
#
# typeof(contRes)
#
# state = iterate(iter)[1]
#     contRes = PALC.initContRes(iter, state)
#     @code_warntype continuation!(iter, state, contRes)
#####

opts = @set opts.doArcLengthScaling = false
br1, sol, _ = @time PALC.continuation(F,Jac_m,x0,-1.5, (@lens _),opts,verbosity=0) #(14.28 k allocations: 1001.500 KiB)
show(br1)

br2, sol, _ = @time PALC.continuation(F,Jac_m,x0,-1.5, (@lens _),opts,verbosity=0, printSolution = (x,p) -> norm(x,2))

# test for different norms
br3, sol, _ = @time PALC.continuation(F,Jac_m,x0,-1.5, (@lens _),opts,verbosity=0, normC = normInf)

# test for linesearch in Newton method
opts = @set opts.newtonOptions.linesearch = true
br4, sol, _ = @time PALC.continuation(F,Jac_m,x0,-1.5, (@lens _),opts,verbosity=0, normC = normInf) # (15.61 k allocations: 1.020 MiB)

# test for different ways to solve the bordered linear system arising during the continuation step
opts = @set opts.newtonOptions.linesearch = false
br5, sol, _ = @time PALC.continuation(F,Jac_m,x0,-1.5, (@lens _),opts,verbosity=0, normC = normInf, linearAlgo = PALC.BorderingBLS())

br5, sol, _ = @time PALC.continuation(F,Jac_m,x0,-1.5, (@lens _),opts,verbosity=0, normC = normInf, linearAlgo = PALC.MatrixBLS())

# test for stopping continuation based on user defined function
finaliseSolution = (z, tau, step, contResult) -> (step < 20)
br5a, sol, _ = @time PALC.continuation(F,Jac_m,x0,-1.5, (@lens _),opts,verbosity=2, finaliseSolution = finaliseSolution)
@test length(br5a.branch) == 21

# test for different predictors
br6, sol, _ = @time PALC.continuation(F,Jac_m,x0,-1.5, (@lens _),opts,verbosity=0, tangentAlgo = PALC.SecantPred())

optsnat = setproperties(opts; ds = 0.001, dsmax = 0.02, dsmin = 0.0001)
br7, sol, _ = @time PALC.continuation(F,Jac_m,x0,-1.5, (@lens _),optsnat,verbosity=0, tangentAlgo = PALC.NaturalPred())

# tangent prediction with Bordered predictor
br8, sol, _ = @time PALC.continuation(F,Jac_m,x0,-1.5, (@lens _),opts,verbosity=0, tangentAlgo = PALC.BorderedPred())

# further testing with sparse Jacobian operator
Jac_sp_simple = (x, p) -> SparseArrays.spdiagm(0 => p  .+ x.^k)
brsp, sol, _ = @time PALC.continuation(F,Jac_sp_simple,x0,-1.5, (@lens _),opts,verbosity=0)
brsp, sol, _ = @time PALC.continuation(F,Jac_sp_simple,x0,-1.5, (@lens _),opts,verbosity=0, printSolution = (x,p) -> norm(x,2))
brsp, sol, _ = @time PALC.continuation(F,Jac_sp_simple,x0,-1.5, (@lens _),opts,verbosity=0,linearAlgo = PALC.BorderingBLS())
brsp, sol, _ = @time PALC.continuation(F,Jac_sp_simple,x0,-1.5, (@lens _),opts,verbosity=0,linearAlgo = PALC.MatrixBLS())
# plotBranch(br1,marker=:d);title!("")
# plotBranch!(br8,marker=:d);title!("")


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module DelehamDelta

using Nemo, NumberTheory, Triangles

export ModuleDelehamDelta
export DeléhamΔ, T084938, T060693, T106566, T094665, T090238, T225478
export T055883, T184962, T088969, T090981, T011117

"""

Philippe Deléham’s Δ-operation maps, similar to the Riordan product, two integer sequences on a lower triangular matrix.
It effectively computes a continued fraction depending on the two input sequences!

Applying Deléham's Δ-operation often gives an additional first column or an additional main diagonal in the resulting triangle compared to what is listed in the OEIS.

[Introduction to the Riordan Square](http://luschny.de/math/seq/RiordanSquare.html)

* DeléhamΔ, T084938, T060693, T106566, T094665, T090238, T225478, T055883, T184962, T088969, T090981, T011117
"""
const ModuleDelehamDelta = ""

"""

Return the product of two integer sequences introduced by Philippe Deléham in A084938.
"""
function DeléhamΔ(n::Int, S::Function, T::Function)
    n ≤ 0 && return fmpz[]
    R, x = PolynomialRing(ZZ, "x")
    A = [R(S(k) + x * T(k)) for k ∈ 0:n-2]
    C = [R(1) for i ∈ 0:n]
    C[1] = R(0)
    M = ZTriangle(n)
    m = 1

    for k ∈ 0:n-1
        for j ∈ k + 1:-1:2
            C[j] = C[j-1] + C[j+1] * A[j-1]
        end
        for j ∈ 0:k
            M[m] = coeff(C[2], j)
            m += 1
        end
    end
    M
end

"""

Return the number of permutations of ``{1,2,...,n}`` having ``k`` cycles such that the elements of each cycle of the permutation form an interval. (Ran Pan)
"""
T084938(n::Int) = DeléhamΔ(n, i -> div(i + 1, 2), i -> 0^i)

"""

Return the number of lattice paths from ``(0,0)`` to ``(x,y)`` that never pass below ``y = x`` and use step set ``{(0,1), (1,0), (2,0), (3,0), ...}``.
"""
T011117(n::Int) = DeléhamΔ(n, i -> 0^i, i -> isodd(i) ? 1 : (i > 0 ? 2 : 0))

"""

Return the number of Schroeder paths (i.e., consisting of steps ``U=(1,1), D=(1,-1), H=(2,0)`` and never going below the x-axis) from ``(0,0)`` to ``(2n,0)``, having ``k`` peaks.
"""
T060693(n::Int) = DeléhamΔ(n, i -> 1, i -> isodd(i) ? 0 : 1)

"""

Return the the Catalan convolution triangle.
"""
T106566(n::Int) = DeléhamΔ(n, i -> i == 0 ? 0 : 1, i -> i == 0 ? 1 : 0)

"""

Return the number of increasing 0-2 trees (A002105) on ``2n`` edges in which the minimal path from the root has length ``k``.
"""
T094665(n::Int) = DeléhamΔ(n, i -> div(i * (i + 1), 2), i -> i + 1)

"""

Return the number of lists of ``k`` unlabeled permutations whose total length is ``n``.
"""
T090238(n::Int) =
    DeléhamΔ(n, i -> div(i, 2) + (isodd(i) ? 2 : 0), i -> i == 0 ? 1 : 0)

"""

Return the triangle ``4^k S_4(n, k)`` where ``S_m(n, k)`` are the Stirling-Frobenius cycle numbers of order ``m``.
"""
T225478(n::Int) =
    DeléhamΔ(n, i -> 2(i + 1) + (i + 1) % 2, i -> isodd(i) ? 0 : 4)

"""

Return the exponential transform of Pascal's triangle.
"""
T055883(n::Int) =
    DeléhamΔ(
        n,
        i -> isodd(i) ? div(i + 1, 2) : 1,
        i -> isodd(i) ? div(i + 1, 2) : 1
    )

"""

Return the number of Schroeder paths of length ``2n`` and having ``k`` ascents.
"""
T090981(n::Int) =
    DeléhamΔ(n, i -> i == 0 ? 1 : (isodd(i) ? 0 : 2), i -> isodd(i) ? 1 : 0)

"""

Return a triangle related to the median Euler numbers.
"""
T088969(n::Int) =
    DeléhamΔ(n, i -> i^2, i -> isodd(i) ? 3div(i, 2) + 2 : 5div(i, 2) + 1)

"""

Return the Bell transform of the Fubini numbers.
"""
T184962(n::Int) =
    DeléhamΔ(
        n,
        i -> div((i + 1) - (i + 1) % 2, 2 - (i + 1) % 2),
        i -> isodd(i) ? 0 : 1
    )

#START-TEST-########################################################
# References to the OEIS A-numbers are always approximately only!

using Test

function test()

    Data = Dict{Int, Array{fmpz}}(
        084938 => [1, 0, 1, 0, 1, 1, 0, 2, 2, 1],
        060693 => [1, 1, 1, 2, 3, 1, 5, 10, 6, 1],
        106566 => [1, 0, 1, 0, 1, 1, 0, 2, 2, 1],
        094665 => [1, 0, 1, 0, 1, 3, 0, 4, 15, 15],
        090238 => [1, 0, 1, 0, 2, 1, 0, 6, 4, 1],
        225478 => [1, 3, 4, 21, 40, 16, 231, 524, 336, 64],
        055883 => [1, 1, 1, 2, 4, 2, 5, 15, 15, 5],
        184962 => [1, 0, 1, 0, 1, 1, 0, 3, 3, 1],
        088969 => [1, 0, 1, 0, 1, 3, 0, 5, 20, 21],
        090981 => [1, 1, 0, 1, 1, 0, 1, 4, 1, 0],
        011117 => [1, 1, 0, 1, 1, 0, 1, 2, 3, 0]
    )

    # [0] 1
    # [1] 1    1
    # [2] 2    4     2
    # [3] 5    15    15    5
    # [4] 15   60    90    60    15
    # [5] 52   260   520   520   260   52
    # [6] 203  1218  3045  4060  3045  1218  203

    @testset "Deléham" begin

        n = 7
        B = [bell(n) * binomial(n, j) for j ∈ 0:n]
        R = Row(T055883(n + 1), n)
        @test all(B .== R)

        a = fmpz[0, 5040, 2208, 828, 272, 70, 12, 1]
        b = Row(T090238(8), 7)
        @test all(a .== b)

        SeqNum(seq) = parse(Int, string(seq)[end-5:end])
        Seq = [
            T084938,
            T060693,
            T106566,
            T094665,
            T090238,
            T225478,
            T055883,
            T184962,
            T088969
        ]

        for seq ∈ Seq
            S = seq(20)
            anum = SeqNum(seq)
            data = Data[anum]
            @test all(S[1:9] .== data[1:9])
        end
    end
end

function demo()
    dim = 5
    for n ∈ 0:dim
        println(T055883(n))
    end

    ShowAsΔ(T055883(dim))
end

"""

T225478(100) :: 0.071950 seconds (15.27 k allocations: 638.438 KiB)
T055883(100) :: 0.076918 seconds (15.27 k allocations: 638.438 KiB)
"""
function perf()
    GC.gc()
    @time T225478(100)
    @time T055883(100)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
"""
    function kinship(ped, i, j)
---
This function is handy if just to calculate relationship of a few (pairs of) ID.
It can also speed up by adding `Thread.@threads` before your pair loop.
"""
function kinship(ped, i, j)
    (i == 0 || j == 0) && return 0
    ipa, ima = ped[i, :]          # used both for below and the last
    i == j && (return 1 + .5kinship(ped, ipa, ima))
    if i < j
        jpa, jma = ped[j, :]
        return .5(kinship(ped, i, jpa) + kinship(ped, i, jma))
    end
    return .5(kinship(ped, j, ipa) + kinship(ped, j, ima))
end


text:
@testset "complex" begin
  (1 + im) / ArbNumerics.ArbComplex(1, 1) == 1
end


text:
# Copyright 2020 Google LLC

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

#     https://www.apache.org/licenses/LICENSE-2.0

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

function parseline(line::String)
    words = split(line)
    string_range = words[1]
    nums = map(s->parse(Int, s), split(string_range, '-'))

    return (nums[1]:nums[2], words[2][1], words[3])
end

parsed = map(parseline, input)

function matches1(tuple)
    req_num, character, password = tuple
    return count(c->c==character, password) ∈ req_num
end

part1(input) = count(matches1, input)

function matches2(tuple)
    req_num, character, password = tuple
    return xor(password[req_num[1]] == character, password[req_num[end]] == character)
end

part2(input) = count(matches2, input)

print(part1(parsed))
print(part2(parsed))


text:
using Pkg
Pkg.activate(".")

using Serialization

filenames = (train_image="train-images-idx3-ubyte",
             train_label="train-labels-idx1-ubyte",
             test_image="t10k-images-idx3-ubyte",
             test_label="t10k-labels-idx1-ubyte")

function read_labels(fname)
    f = open(fname)
    magic = bswap(read(f, UInt32))
    @assert magic == 2049
    nobs = bswap(read(f, UInt32))
    @assert filesize(f) == 8 + nobs
    labels = read(f)
    @assert length(labels) == nobs
    close(f)
    return Int.(labels)
end

function read_images(fname)
    f = open(fname)
    magic = bswap(read(f, UInt32))
    @assert magic == 2051
    nobs = bswap(read(f, UInt32))
    @assert bswap(read(f, UInt32)) == 28
    @assert bswap(read(f, UInt32)) == 28
    @assert filesize(f) == 16 + nobs*28*28
    images = Matrix{Float64}[]
    for i in 1:nobs
        data = read(f, 28*28)
        mx = permutedims(reshape(data, 28, 28)) / 255
        push!(images, mx)
    end
    close(f)
    return images
end

train_images = read_images(filenames.train_image)
train_labels = read_labels(filenames.train_label)
test_images = read_images(filenames.test_image)
test_labels = read_labels(filenames.test_label)

serialize("mnist.data", (train_images, train_labels,
                         test_images, test_labels))


text:
"""

"""
find_circuits(g::DiGraph, ::Type{Val{:johnson}}) = johnson(g)

"""

"""
function johnson(g::DiGraph) end

"""

"""
function findcycles(g::DiGraph, v, s, blocked, stack, B, cycles) end

"""

"""
function unblock!(v, blocked, B)
end


text:
export Qubit, Q0, Q1, RandomQubit

"""
A `Qubit` represents a single qubit.

* `Qubit()` creates the `|0>` value; equivalent to `Qubit(1,0)`.
* `Qubit(0)` and `Qubit(1)` create the `|0>` and `|1>` values, respecively. They are also
  available as `Q0` and `Q1`.
* `Qubit(w::Complex,z::Complex)` creates a qubit whose state is `[w;z]` (normalized)
"""
mutable struct Qubit <: QuantumState
    vec::Vector{CC}
    function Qubit(w, z)
        v = [w; z]
        v /= _norm(v)
        new(v)
    end
end

Qubit(b::Integer = 0) = b == 0 ? Qubit(1, 0) : Qubit(0, 1)

"""
`Q0` is the value `Qubit(1,0)` which is a quantum zero.
"""
const Q0 = Qubit(0)

"""
`Q1` is the value `Qubit(0,1)` which is a quantum one.
"""
const Q1 = Qubit(1)

function show(io::IO, q::Qubit)
    x = q.vec[1]
    y = q.vec[2]
    print(io, "Qubit($x,$y)")
end

"""
`measure!(q::Qubit)` performs a quantum measurement on the qubit `q`
returning `0` or `1` and overwriting `q` to be either `Q0` or `Q1`
respectively.
"""
function measure!(q::Qubit)::Int
    p = abs(q.vec[2])^2
    val::Int = 0
    if rand() < p
        val = 1
        q.vec = Q1.vec
    else
        q.vec = Q0.vec
    end
    return val
end

"""
`RandomQubit()` returns a random qubit.
"""
function RandomQubit()
    a = randn() + randn() * im
    b = randn() + randn() * im
    return Qubit(a, b)
end

length(Q::Qubit) = 1