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, :(swish(-0.3402051397109518 + -0.28456924845213605*swish(-0.42610478134649465 + 0.43223652991075534*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.5575457569666615*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + -0.7868073625875285*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.5078849457963623*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + -0.2883931154790984*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + 0.6888193490874222*swish(-0.17528834888704425 + -0.1425247358809414*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + 0.07932201643104397*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.9443705506849667*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.7094815218028745*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + -0.7352206438426769*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + -0.026837640444652866*swish(-0.005280183583000753 + 0.6172883711874637*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.07641401576689111*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.42701726221355685*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + -0.5680942985776078*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.9418726475599537*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + -0.08739699596173622*swish(0.26121381489212014 + -0.018200625348021138*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.03849358234625422*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + -0.856604860683484*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.2143578715735437*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.706261759527687*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + 0.2360509887003217*swish(-0.09690410604182231 + 0.6473593808554208*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.6752519993458765*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.771479525571356*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.8906247556095672*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.560146720685823*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5])))) + swish(0.2351041904461102 + 0.11979579421065711*swish(-0.42610478134649465 + 0.43223652991075534*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.5575457569666615*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + -0.7868073625875285*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.5078849457963623*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + -0.2883931154790984*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + -0.012124046222294904*swish(-0.17528834888704425 + -0.1425247358809414*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + 0.07932201643104397*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.9443705506849667*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.7094815218028745*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + -0.7352206438426769*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + 0.704332995743346*swish(-0.005280183583000753 + 0.6172883711874637*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.07641401576689111*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.42701726221355685*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + -0.5680942985776078*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.9418726475599537*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + 0.7162932434639226*swish(0.26121381489212014 + -0.018200625348021138*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.03849358234625422*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + -0.856604860683484*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.2143578715735437*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.706261759527687*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + -0.3351546298978829*swish(-0.09690410604182231 + 0.6473593808554208*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.6752519993458765*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.771479525571356*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.8906247556095672*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.560146720685823*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5])))) + swish(-0.4731175976624935 + -0.5144698215629462*swish(-0.42610478134649465 + 0.43223652991075534*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.5575457569666615*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + -0.7868073625875285*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.5078849457963623*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + -0.2883931154790984*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + -0.005430014157027596*swish(-0.17528834888704425 + -0.1425247358809414*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + 0.07932201643104397*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.9443705506849667*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.7094815218028745*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + -0.7352206438426769*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + 0.05201762086802919*swish(-0.005280183583000753 + 0.6172883711874637*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.07641401576689111*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.42701726221355685*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + -0.5680942985776078*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.9418726475599537*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + 0.06728299071980626*swish(0.26121381489212014 + -0.018200625348021138*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.03849358234625422*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + -0.856604860683484*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.2143578715735437*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.706261759527687*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + -0.8893971601390014*swish(-0.09690410604182231 + 0.6473593808554208*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.6752519993458765*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.771479525571356*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.8906247556095672*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.560146720685823*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5])))) + swish(-0.8380258290664089 + 0.49770714260784343*swish(-0.42610478134649465 + 0.43223652991075534*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.5575457569666615*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + -0.7868073625875285*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.5078849457963623*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + -0.2883931154790984*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + 0.732919146629087*swish(-0.17528834888704425 + -0.1425247358809414*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + 0.07932201643104397*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.9443705506849667*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.7094815218028745*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + -0.7352206438426769*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + -0.2983691588163282*swish(-0.005280183583000753 + 0.6172883711874637*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.07641401576689111*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.42701726221355685*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + -0.5680942985776078*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.9418726475599537*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + -0.6200389337758314*swish(0.26121381489212014 + -0.018200625348021138*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.03849358234625422*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + -0.856604860683484*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.2143578715735437*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.706261759527687*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + -0.03199113104777407*swish(-0.09690410604182231 + 0.6473593808554208*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.6752519993458765*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.771479525571356*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.8906247556095672*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.560146720685823*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5])))) + swish(0.7168755118620993 + -0.38856981611164576*swish(-0.42610478134649465 + 0.43223652991075534*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.5575457569666615*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + -0.7868073625875285*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.5078849457963623*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + -0.2883931154790984*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + 0.9818525806719265*swish(-0.17528834888704425 + -0.1425247358809414*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + 0.07932201643104397*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.9443705506849667*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.7094815218028745*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + -0.7352206438426769*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + 0.9400532013818674*swish(-0.005280183583000753 + 0.6172883711874637*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.07641401576689111*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.42701726221355685*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + -0.5680942985776078*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.9418726475599537*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + -0.928831604242299*swish(0.26121381489212014 + -0.018200625348021138*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.03849358234625422*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + -0.856604860683484*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.2143578715735437*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.706261759527687*swish(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + 0.6913287958041021*swish(-0.09690410604182231 + 0.6473593808554208*swish(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.6752519993458765*swish(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.771479525571356*swish(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.8906247556095672*swish(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.560146720685823*swish(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:
# by lazarusA # HIDE
# using GLMakie # HIDE
using ESDL, GLMakie
GLMakie.activate!()# HIDE
# this is a raw approach to plot data cubes
let
    cube = esdc()
    #YAXArray with the following dimensions
    #lon                 Axis with 1440 Elements from -179.875 to 179.875
    #lat                 Axis with 720 Elements from 89.875 to -89.875
    #time                Axis with 1702 Elements from 1980-01-05T00:00:00 to 2016-12-30T00:00:00
    #Variable            Axis with 75 elements: soil_moisture leaf_area_index .. snow_sublimation Rg
    #Total size: 493.03 GB
    #then from here we just call what we need.
    function getCubeFaces(datCubes, timeRange, lonRange, latRange;
        variable = "air_temperature_2m" )
        #front face
        cubef = datCubes[
            var  = variable,
            time = timeRange[1],
            lon  = lonRange,
            lat  = latRange].data[:,:]
        #back face
        cubeb = datCubes[
            var  = variable,
            time = timeRange[2],
            lon  = lonRange,
            lat  = latRange].data[:,:];
        #rigth face
        cuber = datCubes[
            var  = variable,
            time = timeRange,
            lon  = lonRange[2],
            lat  = latRange].data[:,:];
        #left face
        cubel = datCubes[
            var   = variable,
            time  = timeRange,
            lon   = lonRange[1],
            lat   = latRange].data[:,:];
        #top face
        cubet = datCubes[
            var  = variable,
            time = timeRange,
            lon  = lonRange,
            lat  = latRange[2]].data[:,:];
        #bottom face
        cubebt = datCubes[
            var  = variable,
            time = timeRange,
            lon  = lonRange,
            lat  = latRange[1]].data[:,:]

        cubef  = replace(cubef[:,end:-1:1], missing => NaN)
        cubeb  = replace(cubeb[:,end:-1:1], missing => NaN)
        cubet  = replace(cubet', missing => NaN)
        cubebt = replace(cubebt', missing => NaN)
        cubel  = replace(cubel'[:,end:-1:1], missing => NaN)
        cuber = replace(cuber'[:,end:-1:1], missing => NaN) #
        return cubef, cubeb, cuber, cubel, cubet, cubebt
    end

    function cube!(ax; a1 = 1, a2 = 1, timeLonLat = (40,50,20), cxyz = (0,0,0),
        cmap = :plasma, transparent = false, limits = limits,
        topF    = randn(timeLonLat[1], timeLonLat[2]),
        bottomF = randn(timeLonLat[1], timeLonLat[2]),
        frontF  = randn(timeLonLat[2], timeLonLat[3]),
        backF   = randn(timeLonLat[2], timeLonLat[3]),
        leftF   = randn(timeLonLat[1], timeLonLat[3]),
        rightF  =  randn(timeLonLat[1], timeLonLat[3])) # highLimit=0, lowLimit=1

        highLimit = maximum(maximum.(x->isnan(x) ? -Inf : x, (topF, bottomF, frontF, backF, leftF, rightF)))
        lowLimit =  minimum(minimum.(x->isnan(x) ?  Inf : x, (topF, bottomF, frontF, backF, leftF, rightF)))

        x = LinRange(-a1, a1, size(topF)[1])
        y = LinRange(-a1, a1, size(topF)[2])
        z = LinRange(-a1, a1, size(leftF)[2])

        bsurf = heatmap!(ax, y .+ cxyz[2], z .+ cxyz[3], backF, transformation=(:yz, -a2 + cxyz[1]),
        colormap = cmap, transparency = transparent, limits = limits, colorrange = (lowLimit, highLimit)) # back, limits = FRect((-1,-1,-1), (2,2,2))
        heatmap!(ax, x .+ cxyz[1], z .+ cxyz[3], leftF, transformation=(:xz, -a2 + cxyz[2]),
            colormap = cmap, transparency = transparent, colorrange = (lowLimit, highLimit)) # left
        heatmap!(ax, x .+ cxyz[1], y .+ cxyz[2], bottomF, transformation=(:xy, -a2 + cxyz[3]),
            colormap = cmap, transparency = transparent, colorrange = (lowLimit, highLimit)) # bottom
        heatmap!(ax, y .+ cxyz[2], z .+ cxyz[3], frontF, transformation=(:yz,  a2 + cxyz[1]),
            colormap = cmap, transparency = transparent, colorrange = (lowLimit, highLimit)) # front
        heatmap!(ax, x .+ cxyz[1], z .+ cxyz[3], rightF, transformation=(:xz,  a2 + cxyz[2]),
            colormap = cmap, transparency = transparent, colorrange = (lowLimit, highLimit)) # right
        heatmap!(ax, x .+ cxyz[1], y .+ cxyz[2], topF, transformation = (:xy,  a2 + cxyz[3]),
            colormap = cmap, transparency = transparent, colorrange = (lowLimit, highLimit)) # top
        bsurf
    end

    ## Just tell me your variable, time range, lon range and lat range.
    timeRange = (Date("2003-01-01"), Date("2003-12-31"))#(Date("2003-12-31"), Date("2010-01-23"))
    lonRange = (-20, 60)
    latRange = (-40, 40)
    allFaces = getCubeFaces(cube, timeRange, lonRange, latRange; variable = "air_temperature_2m")
    cFront, cBack, cRight, cLeft, cTop, cBottom = allFaces;
    ## Just tell me your variable, time, lon and lat.
    function cubePlot()
        limits = Node(Rect(Vec3f0(-1), Vec3f0(2)))
        fig = Figure(resolution = (850, 800))
        ax = Axis3(fig, aspect = :data, perspectiveness = 0.5, xlabel = "time",
            ylabel = "longitude", zlabel = "latitude")
        cbox = cube!(ax; a1 = 1, a2 = 1, frontF = cFront, backF = cBack,
            topF = cTop, bottomF = cBottom, leftF   = cLeft, rightF  = cRight,
            limits = limits, cmap = :CMRmap, transparent=false)
        ax.yticks = (-1:2:1,string.([lonRange[1], lonRange[2]]))
        ax.zticks = (-1:2:1,string.([latRange[1], latRange[2]]))
        ax.xticks = (-1:2:1,string.([timeRange[1],timeRange[2]]))
        cbar  = Colorbar(fig, cbox,  label = "temperature", height = Relative(0.5))
        #limits[] = Rect(Vec3f0(-2,-2,-2), Vec3f0(4.5))
        fig[1, 1] = ax
        fig[1, 2] = cbar
        fig
    end
    fig = with_theme(cubePlot, theme_black())
    save(joinpath(@__DIR__, "output", "temCubeESDL.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["GLMakie","ESDL"]) # HIDE

ext: jl

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
let
    epc = Epoch(2018,1,1,12,0,0)

    a   = R_EARTH + 500e3
    oe  = [a, 0, 0, 0, 0, sqrt(GM_EARTH/a)]
    eci = sOSCtoCART(oe, use_degrees=true)

    r_rtn   = rRTNtoECI(eci)
    r_rtn_t = rECItoRTN(eci)

    r = r_rtn * r_rtn_t

    tol = 1e-8
    @test isapprox(r[1, 1], 1.0, atol=tol)
    @test isapprox(r[1, 2], 0.0, atol=tol)
    @test isapprox(r[1, 3], 0.0, atol=tol)

    @test isapprox(r[2, 1], 0.0, atol=tol)
    @test isapprox(r[2, 2], 1.0, atol=tol)
    @test isapprox(r[2, 3], 0.0, atol=tol)

    @test isapprox(r[3, 1], 0.0, atol=tol)
    @test isapprox(r[3, 2], 0.0, atol=tol)
    @test isapprox(r[3, 3], 1.0, atol=tol)
end

let
    epc = Epoch(2018,1,1,12,0,0)

    oe  = [R_EARTH + 500e3, 0, 0, 0, 0, 0]
    eci = sOSCtoCART(oe, use_degrees=true)

    xt = deepcopy(eci) + [100, 0, 0, 0, 0, 0]

    x_rtn = sECItoRTN(eci, xt)

    tol = 1e-8
    @test isapprox(x_rtn[1], 100.0, atol=tol)
    @test isapprox(x_rtn[2], 0.0, atol=tol)
    @test isapprox(x_rtn[3], 0.0, atol=tol)
    @test isapprox(x_rtn[4], 0.0, atol=tol)
    @test isapprox(x_rtn[5], 0.0, atol=0.5)
    @test isapprox(x_rtn[6], 0.0, atol=tol)

    xt2 = sRTNtoECI(eci, x_rtn)

    @test isapprox(xt[1], xt2[1], atol=tol)
    @test isapprox(xt[2], xt2[2], atol=tol)
    @test isapprox(xt[3], xt2[3], atol=tol)
    @test isapprox(xt[4], xt2[4], atol=tol)
    @test isapprox(xt[5], xt2[5], atol=tol)
    @test isapprox(xt[6], xt2[6], atol=tol)
end

let
    epc = Epoch(2007, 4, 5, 12, 0, 0, tsys="UTC")

    set_eop(54195.5, -0.072073685, 0.0349282, 0.4833163)

    rc2i = bias_precession_nutation(epc)

    tol = 1e-8
    @test isapprox(rc2i[1, 1], +0.999999746339445, atol=tol)
    @test isapprox(rc2i[1, 2], -0.000000005138822, atol=tol)
    @test isapprox(rc2i[1, 3], -0.000712264730072, atol=tol)

    @test isapprox(rc2i[2, 1], -0.000000026475227, atol=tol)
    @test isapprox(rc2i[2, 2], +0.999999999014975, atol=tol)
    @test isapprox(rc2i[2, 3], -0.000044385242827, atol=tol)

    @test isapprox(rc2i[3, 1], +0.000712264729599, atol=tol)
    @test isapprox(rc2i[3, 2], +0.000044385250426, atol=tol)
    @test isapprox(rc2i[3, 3], +0.999999745354420, atol=tol)
end

let
    epc = Epoch(2007, 4, 5, 12, 0, 0, tsys="UTC")

    set_eop(54195.5, -0.072073685, 0.0349282, 0.4833163)

    r = earth_rotation(epc) * bias_precession_nutation(epc)

    tol = 1e-8
    @test isapprox(r[1, 1], +0.973104317573127, atol=tol)
    @test isapprox(r[1, 2], +0.230363826247709, atol=tol)
    @test isapprox(r[1, 3], -0.000703332818845, atol=tol)

    @test isapprox(r[2, 1], -0.230363798804182, atol=tol)
    @test isapprox(r[2, 2], +0.973104570735574, atol=tol)
    @test isapprox(r[2, 3], +0.000120888549586, atol=tol)

    @test isapprox(r[3, 1], +0.000712264729599, atol=tol)
    @test isapprox(r[3, 2], +0.000044385250426, atol=tol)
    @test isapprox(r[3, 3], +0.999999745354420, atol=tol)
end

let
    epc = Epoch(2007, 4, 5, 12, 0, 0, tsys="UTC")

    set_eop(54195.5, -0.072073685, 0.0349282, 0.4833163)

    r = rECItoECEF(epc)

    tol = 1e-8
    @test isapprox(r[1, 1], +0.973104317697535, atol=tol)
    @test isapprox(r[1, 2], +0.230363826239128, atol=tol)
    @test isapprox(r[1, 3], -0.000703163482198, atol=tol)

    @test isapprox(r[2, 1], -0.230363800456037, atol=tol)
    @test isapprox(r[2, 2], +0.973104570632801, atol=tol)
    @test isapprox(r[2, 3], +0.000118545366625, atol=tol)

    @test isapprox(r[3, 1], +0.000711560162668, atol=tol)
    @test isapprox(r[3, 2], +0.000046626403995, atol=tol)
    @test isapprox(r[3, 3], +0.999999745754024, atol=tol)
end

let
    epc = Epoch(2007, 4, 5, 12, 0, 0, tsys="UTC")

    set_eop(54195.5, -0.072073685, 0.0349282, 0.4833163)

    r = rECEFtoECI(epc)

    tol = 1e-8
    @test isapprox(r[1, 1], +0.973104317697535, atol=tol)
    @test isapprox(r[1, 2], -0.230363800456037, atol=tol)
    @test isapprox(r[1, 3], +0.000711560162668, atol=tol)

    @test isapprox(r[2, 1], +0.230363826239128, atol=tol)
    @test isapprox(r[2, 2], +0.973104570632801, atol=tol)
    @test isapprox(r[2, 3], +0.000046626403995, atol=tol)

    @test isapprox(r[3, 1], -0.000703163482198, atol=tol)
    @test isapprox(r[3, 2], +0.000118545366625, atol=tol)
    @test isapprox(r[3, 3], +0.999999745754024, atol=tol)
end

let
    epc = Epoch(2018,1,1,12,0,0)

    oe  = [R_EARTH + 500e3, 1e-3, 97.8, 75, 25, 45]
    eci = sOSCtoCART(oe, use_degrees=true)

    # Perform circular transformations
    ecef  = sECItoECEF(epc, eci)
    eci2  = sECEFtoECI(epc, ecef)
    ecef2 = sECItoECEF(epc, eci2)

    tol=1e-6
    # Check equivalence of ECI transforms
    @test isapprox(eci2[1], eci[1],  atol=tol)
    @test isapprox(eci2[2], eci[2],  atol=tol)
    @test isapprox(eci2[3], eci[3],  atol=tol)
    @test isapprox(eci2[4], eci[4],  atol=tol)
    @test isapprox(eci2[5], eci[5],  atol=tol)
    @test isapprox(eci2[6], eci[6],  atol=tol)

    # Check equivalence of ECEF transforms
    @test isapprox(ecef2[1], ecef[1], atol=tol)
    @test isapprox(ecef2[2], ecef[2], atol=tol)
    @test isapprox(ecef2[3], ecef[3], atol=tol)
    @test isapprox(ecef2[4], ecef[4], atol=tol)
    @test isapprox(ecef2[5], ecef[5], atol=tol)
    @test isapprox(ecef2[6], ecef[6], atol=tol)
end


text:
"""
3. Veloctiy and stress only (ρ_0 = 1.0)
f₀(x) = f_eq(ρ₀(x), u₀(x), p₀(x))
ρ₀ = 1.0
u₀ = u₀(x)
p₀ = p₀(x)
"""
struct AnalyticalVelocity <: InitializationStrategy end

# TODO compute the density
average_density(problem::FluidFlowProblem) = 1.0
average_pressure(problem::FluidFlowProblem) = 1.0
function density(
    q::Quadrature,
    problem::FluidFlowProblem,
    x,
    y,
    p̄ = average_pressure(problem),
    ρ̄ = average_density(problem),
)
    δρ = (pressure(problem, x, y) - p̄) / q.speed_of_sound_squared

    return ρ̄ + δρ
end
# TODO Move computing the averages to the initialize function so that these can be computed globaly
function initial_condition(
    ::AnalyticalVelocity,
    q::Quadrature,
    problem::FluidFlowProblem,
    x::T,
    y::T,
) where {T <: Real}
    p̄ = average_pressure(problem)
    ρ̄ = average_density(problem)

    # Compute the density from the pressure, average pressure and average density
    ρ = density(q, problem, x, y, p̄, ρ̄)

    u = lattice_velocity(q, problem, x, y), temperature = pressure(problem, x, y) / ρ

    return equilibrium(q, ρ, u, temperature)
end


text:
#Defines a realistic model of a nozzle, to be used in the RealisticThruster type

#Define AbstractNozzle type
abstract type AbstractNozzle end

#= Define Nozzle type, which includes
    - Flag that indicates whether to use tabulated data or a function to determine the
        thrust as function of fuel rate and pressure
    - Matrices that determines the exhaust velocity and pressure as function of fuel rate
    - Function that determines the exhaust velocity and pressure as function of fuel ratee
    - minimum fuel rate for choked flow
    - maximum fuel rate
    - Exhaust Area

=#
struct Nozzle <: AbstractNozzle
    ExhaustVelocityData
    ExhaustPressureData
    DataFlag::Bool
    ExhaustArea::Float64
    MinFuelRate::Float64
    MaxFuelRate::Float64

    #Constructor function
    function Nozzle(ExhaustVelocityData, ExhaustPressureData, DataFlag::Bool,
        ExhaustArea::Number, MinFuelRate::Number, MaxFuelRate::Number)
        if DataFlag == false
            if ((typeof(ExhaustVelocityData) != Matrix{Float64}) ||
                (typeof(ExhaustPressureData) != Matrix{Float64}))
                error("Expected Exhaust velocity and pressure data to be matrices.")
            end
        else
            if ((typeof(ExhaustVelocityData) != Function ||
                (typeof(ExhaustPressureData) != Function))
                error("Expected Exhaust velocity and pressure data to be functions.")
            end
        end
        new(ExhaustVelocityData, ExhaustPressureData, DataFlag, Float64(ExhaustArea),
                Float64(MinFuelRate), Float64(MaxFuelRate))
    end #Constructor

end

function Nozzle()
    return nothing
end


text:
# utility functions needed for discerning topology from geometry

using Combinatorics

using dataset_pb2

export bohr_to_angstroms, angstroms_to_bohr, distance_between_atoms, bonded, distances
export add_atom!, add_bond!
export smu_atype_to_atomic_number, number_to_smu_btype, smu_btype_to_number, smu_atom_type_to_max_con
export atomic_number_to_smu_atom_type
export add_atom_by_atomic_number!, add_atom_by_type!
export distances, bonded
export canonical_bond_topology!, same_bond_topology
export is_single_fragment

function bohr_to_angstroms(distance)
  return distance * 0.529177249
end

function angstroms_to_bohr(distance)
  return distance / 0.529177249
end

function distance_between_atoms(geom::dataset_pb2.Geometry, a1, a2)::Float32
  return bohr_to_angstroms(sqrt(
    (geom.atom_positions[a1].x - geom.atom_positions[a2].x) *
    (geom.atom_positions[a1].x - geom.atom_positions[a2].x) +
    (geom.atom_positions[a1].y - geom.atom_positions[a2].y) *
    (geom.atom_positions[a1].y - geom.atom_positions[a2].y) +
    (geom.atom_positions[a1].z - geom.atom_positions[a2].z) *
    (geom.atom_positions[a1].z - geom.atom_positions[a2].z)
  ))
end

_smu_atype_to_atomic_number = Dict{Int32, Int32}(
  BondTopology_AtomType.ATOM_H => 1,
  BondTopology_AtomType.ATOM_C => 6,
  BondTopology_AtomType.ATOM_N => 7,
  BondTopology_AtomType.ATOM_NPOS => 7,
  BondTopology_AtomType.ATOM_O => 8,
  BondTopology_AtomType.ATOM_ONEG => 8,
  BondTopology_AtomType.ATOM_F => 9
)

function smu_atype_to_atomic_number(smu_type::T)::Int32 where {T<:Integer}
  return _smu_atype_to_atomic_number[smu_type]
end

"""Given an atomic number and charge, return the BondTopology_AtomType
Args:
  atomic_number:
  charge:
Returns:
  BondTopology_AtomType
"""
function atomic_number_to_smu_atom_type(atomic_number::Integer, charge::Integer=0)::Int32
  if atomic_number == 6
    return BondTopology_AtomType.ATOM_C
  elseif atomic_number == 7
    if charge == 0
      return BondTopology_AtomType.ATOM_N
    else
      return BondTopology_AtomType.ATOM_NPOS
    end
  elseif atomic_number == 8
    if charge == 0
      return BondTopology_AtomType.ATOM_O
    else
      return BondTopology_AtomType.ATOM_ONEG
    end
  elseif atomic_number == 9
    return BondTopology_AtomType.ATOM_F
  elseif atomic_number == 1
    return BondTopology_AtomType.ATOM_H
  else
    @error("Unrecognized atomic number $(atomic_number)")
  end
end

_number_to_smu_btype = Dict{Int32, Int32}(
  0 => BondTopology_BondType.BOND_UNDEFINED,
  1 => BondTopology_BondType.BOND_SINGLE,
  2 => BondTopology_BondType.BOND_DOUBLE,
  3 => BondTopology_BondType.BOND_TRIPLE
)

function number_to_smu_btype(n::T)::Int32 where {T<:Integer}
  return _number_to_smu_btype[n]
end

_smu_btype_to_number = Dict{Int32, Int32}(
  BondTopology_BondType.BOND_UNDEFINED => 0,
  BondTopology_BondType.BOND_SINGLE => 1,
  BondTopology_BondType.BOND_DOUBLE => 2,
  BondTopology_BondType.BOND_TRIPLE => 3
)

_smu_atom_type_to_max_con = Dict{Int32, Int32}(
  BondTopology_AtomType.ATOM_H => 1,
  BondTopology_AtomType.ATOM_C => 4,
  BondTopology_AtomType.ATOM_N => 3,
  BondTopology_AtomType.ATOM_NPOS => 3,
  BondTopology_AtomType.ATOM_O => 2,
  BondTopology_AtomType.ATOM_ONEG => 1,
  BondTopology_AtomType.ATOM_F => 1
)

# Return the maximum number of bonds for `smu_atype`.
# Unfortunately the type declaration on the next line does not work.
# function smu_atom_type_to_max_con(smu_atype::BondTopology_AtomType)::Int32

function smu_atom_type_to_max_con(smu_atype::T)::Int32 where {T <:Integer}
  return _smu_atom_type_to_max_con[smu_atype]
end

function smu_btype_to_number(smu_btype::T)::Int32 where {T<:Integer}
  return _smu_btype_to_number[smu_btype]
end

function add_bond!(a1::T, a2::T, btype, bond_topology::dataset_pb2.BondTopology) where {T<:Integer}
  smu_btype = number_to_smu_btype(btype)

  new_bond = BondTopology_Bond(atom_a=a1, atom_b=a2, bond_type=smu_btype)

  if !hasproperty(bond_topology, :bonds)
    setproperty!(bond_topology, :bonds, [new_bond])
  else
    push!(bond_topology.bonds, new_bond)
  end
end

function add_atom_by_atomic_number!(bond_topology::dataset_pb2.BondTopology, atomic_number::T, charge::T=0) where {T<:Integer}
  smu_atype = atomic_number_to_smu_atom_type(atomic_number, charge)
  add_atom_by_type!(bond_topology, smu_atype)
end

function add_atom_by_type!(bond_topology::dataset_pb2.BondTopology, smu_atype::T) where {T <:Integer}
  if !hasproperty(bond_topology, :atoms)
    setproperty!(bond_topology, :atoms, [smu_atype])
  else
    push!(bond_topology.atoms, smu_atype)
  end
end

"""Return an int array of the bonded atoms in `bond_topology`.
The returned array is square[natoms,natoms]. If there is no bond
between the atoms, the value is zero. If the atoms are bonded,
the entry will be the bond type.
Args:
  bond_topology:
Returns:
  Vector of BondType's
"""
function bonded(bond_topology::dataset_pb2.BondTopology)::Array{Int32, 2}
  natoms = length(bond_topology.atoms)
  connected = zeros(Int32, (natoms, natoms))  # to be returned
  # Note need to convert to 1 indexing for Julia
  for bond in bond_topology.bonds
    a1 = bond.atom_a + 1
    a2 = bond.atom_b + 1
    connected[a1, a2] = connected[a2, a1] = bond.bond_type
  end
  return connected
end

"""Given a connection matrix `bonded`, return a vector of
vectors, one for each atom, containing the indices of the
bonded atoms.
Args:
  bonded:
Returns:
"""
function connections(bonded::Array{Int32, 2})::Vector{Vector{Int32}}
  natoms = size(bonded, 1)
  result = Vector{Vector{Int32}}(undef, (natoms))
  for i in 1:natoms
    nbrs = Vector{Int32}()
    for j in 1:natoms
      if bonded[i,j] > 0
        push!(nbrs, j)
      end
    end
    result[i] = nbrs
  end
  return result
end

"""Return a float array of the interatomic distances in `geometry`.
Args:
  geometry:
Returns:
  a numpy array of distances
"""
function distances(geometry)::Array{Float32, 2}
  natoms = length(geometry.atom_positions)
  result = zeros(Float32, (natoms, natoms))
  for atoms in combinations(1:natoms, 2)
    i = atoms[1]
    j = atoms[2]
    result[i, j] = result[j, i] = distance_between_atoms(geometry, i, j)
  end
  return result
end

"""Transform the bonds attribute of `bond_topology` to a canonical form.

Args:
  bond_topology:
Returns:
  BondTopology
"""
function canonical_bond_topology!(bond_topology)
  if length(bond_topology.bonds) < 2
    return
  end

  for bond in bond_topology.bonds
    if bond.atom_a > bond.atom_b
      bond.atom_a, bond.atom_b = bond.atom_b, bond.atom_a
    end
  end

  sort!(bond_topology.bonds, by = b -> (b.atom_a, b.atom_b))
end

Base.:(==)(b1::BondTopology_Bond, b2::BondTopology_Bond)::Bool = b1.atom_a == b2.atom_a &&
                b1.atom_b == b2.atom_b && b1.bond_type == b2.bond_type

"""Return True if bt1 == bt2.
Note that there is no attempt to canonialise the protos.
Args:
Returns:
"""
function same_bond_topology(bt1, bt2)::Bool
  length(bt2.atoms) == length(bt2.atoms) || return false
  length(bt2.bonds) == length(bt2.bonds) || return false

  # https://stackoverflow.com/questions/47564825/check-if-all-the-elements-of-a-julia-array-are-equal
  # all(y->y==x[1], x);

  atoms1 = bt1.atoms
  if ! all(atoms1->atoms1==bt2.atoms[1], bt2.atoms)
    return false
  end

  # This should work, but raises an exception.
  # bonds1 = bt1.bonds
  # bonds2 = bt2.bonds
  # return all(bonds1->bonds1==bonds2[1], bonds2)

  for i in 1:length(bt1.bonds)
    bt1.bonds[i] != bt2.bonds[i] && return false
  end
  return true
end

"""Recusrively visit nodes in the graph defined by `nbrs`.
Args:
  nbrs:
  atom:
  visited:
Returns:
  The number of nodes visited - including `atom`.
"""
function visit(nbrs::Vector, atom::T, visited::Vector{Bool})::Int32 where {T<:Integer}
  visited[atom] = true
  result = 1    # To be returned.
  for nbr in nbrs[atom]
    if visited[nbr] > 0
      continue
    end
    result += visit(nbrs, nbr, visited)
  end

  return result
end

"""Return True if `bond_topology` is a single fragment.
Args:
  bond_topology:
Returns:
  True if `bond_topology` is a single fragment.
"""
function is_single_fragment(bond_topology::BondTopology)::Bool
  natoms = length(bond_topology.atoms)
  nbonds = length(bond_topology.bonds)
  # Some special cases are easy.
  natoms == 1 && return true
  if natoms == 2 && nbonds == 1
    return true
  end
  if natoms == 3 && nbonds == 2
    return true
  end
  if natoms == nbonds && natoms <= 4
    return true
  end

  attached = connections(bonded(bond_topology))
  # Any atom with zero neighbors means a detached atom.
  any(n->length(n) == 0, attached) && return false

  visited = zeros(Bool, natoms)
  # Mark anything with a single connection as visited.
  # Record the index of an atom that has multiple connections.
  a_multiply_connected_atom = -1
  for i in 1:natoms
    if bond_topology.atoms[i] == BondTopology_AtomType.ATOM_H
      visited[i] = 1
      continue
    end

    if length(attached[i]) > 1
      a_multiply_connected_atom = i
      continue
    end

    # A singly connected heavy atom. Mark visited if not part of a two atom fragment.
    if length(attached[attached[i][1]]) > 1
      visited[i] = 1
    end
  end

  # Not sure this can happen.
  a_multiply_connected_atom < 0 && return false

  number_visited = sum(visited) + visit(attached, a_multiply_connected_atom, visited)
  return number_visited == natoms
end

max_stars_repo_name: IanAWatson/smu

max_issues_repo_name: IanAWatson/smu

max_forks_repo_name: IanAWatson/smu


text:
# -------------- #
# Symbolic Model #
# -------------- #
#
function _get_args(sm, spec)
    # get args
    args = OrderedDict{Symbol,Vector{Tuple{Symbol,Int}}}()
    for (grp, shift, nm) in spec[:eqs]
        grp == "parameters" && continue
        syms = sm.symbols[Symbol(grp)]
        args[Symbol(nm)] = Tuple{Symbol,Int}[(v, shift) for v in syms]
    end
    args
end

# ----- #
# Tools #
# ----- #

# Used for extracting complementarity conditions from arbitrage equations
function _handle_arbitrage(arb, controls)
    controls_lb = Expr[]
    controls_ub = Expr[]
    arbitrage = Expr[]
    for (i, v) in enumerate(arb)
        parts = split(v, "|")

        if length(parts) == 1
            push!(arbitrage, _to_expr(parts[1]))
            push!(controls_lb, _to_expr(-Inf))
            push!(controls_ub, _to_expr(Inf))

        # We have complementarity conditions to deal with
        elseif length(parts) == 2
            # convert arbitrage equation to an expression
            push!(arbitrage, _to_expr(parts[1]))

            # Now deal with complementarities
            c_parts = split(parts[2], "<=")
            n_c = length(c_parts)

            if n_c == 3
                push!(controls_lb, _to_expr(inf_to_Inf(Meta.parse(c_parts[1]))))
                push!(controls_ub, _to_expr(inf_to_Inf(Meta.parse(c_parts[3]))))

                # verify that the control is what we want
                mid = Meta.parse(c_parts[2])
                if mid != controls[i]
                    msg = string("Error in complementarity condition. ",
                                 "Expected $(controls[i]) found $mid")
                    error(msg)
                end
            elseif n_c == 2
                # only have one_sided condition. Need to work a bit harder
                ind = 0
                for (j, ex) in enumerate(c_parts)
                    if Meta.parse(ex) == controls[i]
                        ind = j
                        break
                    end
                end

                # we didn't find the control, throw an error
                if ind == 0
                    msg = string("Error in complementarity condition. ",
                                 "Expected $(controls[i]), but did not find.")
                    error(msg)
                end

                # otherwise continue on
                if ind == 1
                    push!(controls_lb, _to_expr(-Inf))
                    push!(controls_ub, _to_expr(inf_to_Inf(parse(c_parts[2]))))
                else
                    push!(controls_lb, _to_expr(inf_to_Inf(parse(c_parts[1]))))
                    push!(controls_ub, _to_expr(Inf))
                end
            else
                msg = string("Malformed complementarity condition. ",
                             "Expected 1 or 2 `<=`, found $(n_c-1)")
                error(msg)
            end

        else

            msg = string("Malformed arbitrage equation. ",
                         "Only have one occurance of `|` allowed.")
            error(msg)
        end

    end
    controls_lb, controls_ub, arbitrage
end


text:
#using GLMakie
#AbstractPlotting.inline!(true)
using CairoMakie

#using ColorSchemes
#using FixedPointNumbers

module DrifterViz

using OceanRobots, MeshArrays, CSV, DataFrames, Makie
import Proj4

export proj_map, subset_positions, animate_positions
export background_stuff, update_stuff!

"""
    animate_positions()

Plot positions in a loop, generate a movie

```
#tt=collect(2000:0.1:2003)
animate_positions(B)
```
"""
function animate_positions(B)
    time = Node(2005.1)
    #timestamps = 2005.1:0.01:2018.9
    timestamps = 0.0:0.01:3.0
    framerate = 20

    #dt=0.1
    #ii = @lift(findall( (B.df[1].t .> $time-dt).*(B.df[1].t .< $time+dt) ))
    #tt = @lift( string($time) )

    dt=3.0
    dd=@lift( 360.0*($time-floor.($time)) )
    ii = @lift(findall( (abs.(B.df[1].d .- $dd).<dt).|(abs.(B.df[1].d .- $dd).> 360.0-dt) ))
    tt = @lift( string(Int(floor($dd))) )

    nmax=100000
    x = @lift( [B.df[1].x[$ii];fill(NaN,nmax-length($ii))])
    y = @lift( [B.df[1].y[$ii];fill(NaN,nmax-length($ii))])
    if true
        c = @lift( [ sqrt.(B.df[1].u[$ii].^2+B.df[1].v[$ii].^2) ;fill(NaN,nmax-length($ii))])
        cr=(0.,5.0)
    else
        c = @lift( [B.df[1].u[$ii];fill(NaN,nmax-length($ii))])
        cr = (-0.4, 0.4)
    end

    F=DrifterViz.proj_map(B)
    ax = F[1, 1]
    pnts = scatter!(ax, x, y, show_axis = false, color=c, colorrange=cr,
    colormap=:turbo, markersize=2, strokewidth=0.0)

    xy=Proj4.transform(B.source, B.dest, [70.0 50.0])
    txt=text!(ax,tt,position = (xy[1],xy[2]),color=:red) #textsize = 1e6

    N=record(F, joinpath(tempdir(),"tmp.mp4"), timestamps; framerate = framerate) do t
        time[] = t
        #println(t)
    end

    return F,N
end

"""
    update_stuff!(B,t0,t1)
"""
function update_stuff!(B,t0,t1)
    tst=(maximum(B.df[1].t)>t1)||(minimum(B.df[1].t)<t0)
    if tst
        pth,lst=DrifterViz.drifters_ElipotEtAl16()
        tmp=readdir(pth)
        tmp=tmp[findall(occursin.("csv/drifters_",tmp))]
        y0=Int(floor(minimum(B.df[1].t)))

        df1=DataFrame(CSV.File(pth*"csv/drifters_$(y0).csv"))
        df2=DataFrame(CSV.File(pth*"csv/drifters_$(y0+1).csv"))
        B.df[1]=vcat(df1,df2)

        tt=sort(unique(B.df[1].t))[1:24:end]
        B.df[1]=filter(row -> sum(row.t .== tt)>0 ,B.df[1])
    end
end

"""
    subset_positions(B,t0,t1)
"""
#update the observables
function subset_positions(B,t0,t1)
    #update_stuff!(B,t0,t1)
    #
    df=B.df[1]
    df=df[df.t .> t0,:]
    df=df[df.t .<= t1,:]

    #
    nmax=100000
    xy=fill(NaN,nmax,2)
    npos=length(df.lon)
    npos>nmax ? println("nmax exceeded; skipping data") : nothing #println("$npos")
    npos=min(npos,nmax)
    xy[1:npos,:]=Proj4.transform(B.source, B.dest, [vec(df.lon[1:npos]) vec(df.lat[1:npos])])
    #
    return xy
end

"""
    proj_map(B::NamedTuple)

Geographic projection map of c using geomakie

```
B=DrifterViz.background_stuff();
F=DrifterViz.proj_map(B);
#df=DrifterViz.subset_positions(B,2005.15,2005.16)
F,N=DrifterViz.animate_positions(B)
F
```
"""
function proj_map(B::NamedTuple)
#    @unpack x, y, DL= B

    col=deepcopy(B.z)
    msk=findall(isnan.(B.z))
    col[msk].=Inf

#    if isa(B.rng,Tuple)
#        col=get(ColorSchemes.balance, col, B.rng)
#        col=to_color.(col)
#        col[msk].=RGBA{Float64}(1.,1.,1.,0.)
#    end

    f = Figure()
    ax = f[1, 1] = Axis(f)

    surf = surface!(ax,B.x,B.y,0*B.x; color=col, colorrange = B.rng, colormap=:grayC,
        shading = false, scale_plot = false, axis = (xticks = [0.0],))
    #cbar=Colorbar(f,surf, width = 20)
    #f[1, 2] = cbar

    ii=[i for i in -180:45:180, j in -78.5:1.0:78.5]';
    jj=[j for i in -180:45:180, j in -78.5:1.0:78.5]';
    xl=vcat([[ii[:,i]; NaN] for i in 1:size(ii,2)]...)
    yl=vcat([[jj[:,i]; NaN] for i in 1:size(ii,2)]...)
    tmp=Proj4.transform(B.source, B.dest,[xl[:] yl[:]])
    xl=tmp[:,1]; yl=tmp[:,2]
    lines!(xl,yl, color = :white, linewidth = 0.5)

    tmp=circshift(-179.5:1.0:179.5,(-200))
    ii=[i for i in tmp, j in -75:15:75];
    jj=[j for i in tmp, j in -75:15:75];
    xl=vcat([[ii[:,i]; NaN] for i in 1:size(ii,2)]...)
    yl=vcat([[jj[:,i]; NaN] for i in 1:size(ii,2)]...)
    tmp=Proj4.transform(B.source, B.dest,[xl[:] yl[:]])
    xl=tmp[:,1]; yl=tmp[:,2]
    lines!(xl,yl, color = :white, linewidth = 0.5)

    hidespines!(ax)
    hidedecorations!.(ax)

    return f
end

##

"""
    background_stuff()
"""
function background_stuff()
    Γ=GridLoad(GridSpec("LatLonCap",MeshArrays.GRID_LLC90))

#Global Ocean
    lo=[collect(20.5:1.0:179.5); collect(-179.5:1.0:19.5)]
    la=collect(-78.5:0.5:78.5)

#North Pacific
#    lo=[collect(110.0:0.5:180.0); collect(-180.0:0.5:-90.0)]
#    la=collect(-20.0:0.25:50.0)

    lon=[i for i in lo, j in la]
    lat=[j for i in lo, j in la]

    (f,i,j,w)=InterpolationFactors(Γ,vec(lon),vec(lat))

    DL=log10.(Interpolate(Γ.Depth,f,i,j,w))
    DL[findall((!isfinite).(DL))].=NaN
    DL=reshape(DL,size(lon))
    rng=(3.,5.)

    ##

    source=Proj4.Projection("+proj=longlat +datum=WGS84")
    dest=Proj4.Projection("+proj=eqearth +lon_0=200.0 +lat_1=0.0 +x_0=0.0 +y_0=0.0 +ellps=GRS80")

    tmp=Proj4.transform(source, dest, [lon[:] lat[:]])
    x=reshape(tmp[:,1],size(lon))
    y=reshape(tmp[:,2],size(lon))
    z=DL

    ##

    pth,lst=drifters_ElipotEtAl16();

    y0=2005
    y1=2018
    df=[DataFrame()]
    for y=y0:y1
        println("loading $(y)")
        tmp=DataFrame(CSV.File(pth*"csv/drifters_$(y).csv"))
        tt=sort(unique(tmp.t))[1:24:end]
        tmp=filter(row -> sum(row.t .== tt)>0 ,tmp)

        #geographical projection for plotting
        xy=Proj4.transform(source, dest, [vec(tmp.lon) vec(tmp.lat)])
        tmp.x=xy[:,1]
        tmp.y=xy[:,2]

        #pseudo-day of 360-day-year
        tmp.d=360.0*(tmp.t-floor.(tmp.t))

        df[1]=vcat(df[1],tmp)
    end

    return (;lon,lat,DL,x,y,z,rng,source,dest,df)
end

end


text:
#=

  Euler #9 in Julia.

  Problem 9
  """
  A Pythagorean triplet is a set of three natural numbers, a  b  c, for which,
  a^2 + b^2 = c^2

  For example, 3^2 + 4^2 = 9 + 16 = 25 = 5^2.

  There exists exactly one Pythagorean triplet for which a + b + c = 1000.
  Find the product abc.
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/Julia/

=#

include("Euler.jl")

# 0.00683682s
function euler9a()
    for c in 1:500
        for b in 1:c
            for a in 1:b
                if a+b+c === 1000 && a^2 + b^2 - c^2 == 0
                    return a*b*c;
                end
            end
        end
    end
end

run_euler(euler9a)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
function tikzlattice(model::UnitCellQModel,
                     filename::String,
                     coordchange::Bool=false)
    isfile(filename) && error("file already exists: $filename")
    D = dimension(model)
    D == 2 || error("only 2D lattices can be drawn with tikz!")
    lattice = model.lattice

    bondcolors = [
        "blue", "orange", "cyan", "olive", "red", "brown", "green", "yellow" ]
    ballcolors = [
        "green", "blue", "red", "yellow", "brown", "cyan", "orange", "olive"]

    a1 = lattice.unitc.as[1]
    a2 = lattice.unitc.as[2]

    open(filename, "w") do f
        write(f, """\\documentclass[tikz, border=1mm]{standalone}
  \\usepackage{tikz,pgf}
  \\usetikzlibrary{calc,arrows,arrows.meta}
  \\usetikzlibrary{decorations.markings}

  \\begin{document}

  \\begin{tikzpicture}
  """)
        if coordchange
            write(f, "[
        x={$a1},
        y={$a2},
        ]\n")
        end

        write(f, "  \\tikzset{\n")

        bonds  = [inter for inter in model.inters if support(inter)==2]
        n3bodys = [inter for inter in model.inters if support(inter)==3]
        n4odys = [inter for inter in model.inters if support(inter)==4]

        for i in eachindex(bonds)
            write(f, "    bond$i/.style={thin, double, $(bondcolors[i])!70!black},\n")
        end
        for i in 1:lattice.unitc.n
            write(f, "    ball$i/.style={inner color=blue, ball color=$(ballcolors[i])!20!black},\n")
        end

        write(f, """
    ball/.style={inner color=blue, ball color=green!20!black}
  }

  \\def\\l{1cm}
  \\def\\radius{\\l/10}\n
""")

        pad = 0.3
        ly, lx = model.lattice.sizes
        write(f, "  \\clip (-$pad, -$pad) -- ++(0.0, $(ly-1+2*pad)) -- ++($(lx-1+2*pad), 0.0) -- ++(0.0, -$(ly-1+2*pad)) -- cycle;\n")

        for is in Iterators.product([0:l+1 for l in model.lattice.sizes]...)
            for index in eachindex(bonds)
                interaction = bonds[index]
                #support(interaction) == 2 || continue
                n1, n2 = interaction.ucidxs
                off1, off2 = interaction.offsets
                #if lattice.bc == :OBC
                x_uc1 = [off1[i] .+ is[i] for i=1:D]
                x_uc2 = [off2[i] .+ is[i] for i=1:D]
                if all([0 <= x_uc1[i] <= lattice.sizes[i]+1 for i=1:D]) &
                    all([0 <= x_uc2[i] <= lattice.sizes[i]+1 for i=1:D])
                    one = reduce(coordsum,
                                 #[(x_uc1 .- 1)[i] .* lattice.unitc.as[i] for i=1:D],
                                 [reverse(x_uc1 .- 1)[i] .* lattice.unitc.as[i] for i=1:D],
                                 init=(0,0)).+ lattice.unitc.sites[n1]
                    two = reduce(coordsum,
                                 #[(x_uc2 .- 1)[i] .* lattice.unitc.as[i] for i=1:D],
                                 [reverse(x_uc2 .- 1)[i] .* lattice.unitc.as[i] for i=1:D],
                                 init=(0,0)).+ lattice.unitc.sites[n2]
                    write(f, "  \\draw [bond$index] ($(one[1])*\\l, $(one[2])*\\l) -- ($(two[1])*\\l, $(two[2])*\\l);\n")
                end
                #else
                #    error("boundary not supported yet!")
                #end
            end
        end

        for ls in Iterators.product([1:l for l in model.lattice.sizes]...)
            x_uc = reduce(coordsum,
                          #[(ls .- 1)[i] .* lattice.unitc.as[i] for i=1:D],
                          [reverse(ls .- 1)[i] .* lattice.unitc.as[i] for i=1:D],
                          init=(0,0))
            for n in 1:lattice.unitc.n
                x_site = x_uc .+ lattice.unitc.sites[n]
                write(f, "  \\shade [ball$n] ($(x_site[1])*\\l, $(x_site[2])*\\l) circle (\\radius);\n")
            end
        end

        for index in eachindex(bonds)
            interaction = bonds[index]
            ns = interaction.ucidxs
            offs = interaction.offsets
            for is in Iterators.product([1:l for l in model.lattice.sizes]...)
                x_ucs = [[offs[n][i] .+ is[i] for i=1:D] for n in 1:length(ns)]
                if all([
                    all([1 <= x_ucs[n][i] <= lattice.sizes[i]
                         for i=1:D])
                    for n in 1:length(ns)])
                    sites = [reduce(coordsum,
                                    #[(x_uc1 .- 1)[i] .* lattice.unitc.as[i] for i=1:D],
                                    [reverse(x_ucs[n] .- 1)[i] .* lattice.unitc.as[i] for i=1:D],
                                    init=(0,0)) .+ lattice.unitc.sites[ns[n]]
                             for n in 1:length(ns)]
                    #write(f, "  \\draw [->, shorten >=\\radius, shorten <=\\radius] ($(sites[1][1])*\\l, $(sites[1][2])*\\l) to [out=160, in=20, looseness=1] ($(sites[2][1])*\\l, $(sites[2][2])*\\l);\n")
                    mid = (sites[1] .+ sites[2]) ./ 2
                    write(f, "  \\node at ($(mid[1])*\\l, $(mid[2])*\\l) {\$t_$index\$};\n")
                    break
                end
            end
        end
        write(f, "\\end{tikzpicture}\n")
        write(f, "\\end{document}")
    end
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
using Documenter;

if get(ENV, "CIRCLECI", "") == "true"
	sha = try
		readchomp(`git rev-parse --short HEAD`)
	catch
		"(not-git-repo)"
	end

	mktempdir() do temp
		Documenter.git_push(Documenter.Utilities.currentdir(), temp, "github.com/marekdedic/ThreadedIterables.jl.git"; target = "docs/build", tag = get(ENV, "CIRCLE_TAG", ""), sha = sha, versions = ["stable" => "v^", "v#.#", "dev" => "dev"]);
	end
end


text:
# Autogenerated wrapper script for FlowCutterPACE17_jll for i686-w64-mingw32-cxx03
export flow_cutter_pace17

JLLWrappers.@generate_wrapper_header("FlowCutterPACE17")
JLLWrappers.@declare_executable_product(flow_cutter_pace17)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        flow_cutter_pace17,
        "bin\\flow_cutter_pace17.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# ROCm-specific runtime libraries

## GPU runtime library

# reset the runtime cache from global scope, so that any change triggers recompilation
GPUCompiler.reset_runtime()

function signal_exception()
    flag_ptr = get_global_pointer(Val(:__global_exception_flag), Int64)
    unsafe_store!(flag_ptr, 1)
    # TODO: Can we sync here?
    endpgm() # stop this wavefront
    return
end

function report_exception(ex)
    # Add kernel ID and exception string to exception ring buffer
    ring_ptr = get_global_pointer(Val(:__global_exception_ring), LLVMPtr{ExceptionEntry,AS.Global})
    ring_ptr = unsafe_load(ring_ptr)
    our_signal = _completion_signal()
    prev = UInt64(1)
    while prev != UInt64(0)
        # Try to write to this slot, and skip if we fail (because another wavefront wrote first)
        prev = atomic_cas!(reinterpret(LLVMPtr{UInt64,AS.Global}, ring_ptr), UInt64(0), our_signal)
        if prev == UInt64(0)
            # We get a ReadOnlyMemoryError without this copy because the data is pinned to the device
            # TODO: Don't use an expensive host malloc
            ex_len = string_length(ex)::Csize_t
            copy_buf = reinterpret(LLVMPtr{UInt8,AS.Global}, malloc(ex_len))
            memcpy!(copy_buf, reinterpret(LLVMPtr{UInt8,AS.Global}, ex), ex_len)
            Base.unsafe_store!(reinterpret(LLVMPtr{UInt64,AS.Global}, ring_ptr+sizeof(UInt64)), UInt64(copy_buf))
            break
        elseif prev == UInt64(1)
            # Tail slot, give up
            break
        end
        ring_ptr += sizeof(ExceptionEntry)
    end
    return
end

report_oom(sz) = @rocprintf("ERROR: Out of dynamic GPU memory (trying to allocate %i bytes)\n", sz)

function device_string_to_host(ex)
    ex_ptr = reinterpret(LLVMPtr{UInt8,1}, ex)
    ex_len = string_length(ex_ptr)
    ex_str = reinterpret(LLVMPtr{UInt8,1}, malloc(Csize_t(ex_len+1)))
    if reinterpret(UInt64, ex_str) == 0
        @rocprintf("Device-to-host string conversion failed\n")
        return reinterpret(Cstring, 0)
    end
    memcpy!(ex_str, ex_ptr, ex_len)
    unsafe_store!(ex_str+ex_len, UInt8(0))
    reinterpret(Cstring, ex_str)
end

function report_exception_name(ex)
    # Pass argument in host buffer
    ex_str = device_string_to_host(ex)
    @rocprintf("""
        ERROR: a %s was thrown during kernel execution.
        Stacktrace:
        """, ex_str)
    free(reinterpret(Ptr{Cvoid}, ex_str))
    return
end

function report_exception_frame(idx, func, file, line)
    # Pass arguments in host buffers
    func_str = device_string_to_host(func)
    file_str = device_string_to_host(file)
    @rocprintf(" [%i] %s at %s:%i\n", idx, func_str, file_str, line)
    free(reinterpret(Ptr{Cvoid}, func_str))
    free(reinterpret(Ptr{Cvoid}, file_str))
    return
end

## ROCm device library

function load_and_link!(mod, path)
    @debug "Loading device library" path
    ctx = LLVM.context(mod)
    lib = parse(LLVM.Module, read(path); ctx)

    for f in LLVM.functions(lib)
        # FIXME: We should be able to inline this, that we can't means
        #        we are inserting calls to it late.
        LLVM.name(f) == "__ockl_hsa_signal_store" && continue
        LLVM.name(f) == "__ockl_hsa_signal_load" && continue
        attrs = function_attributes(f)
        push!(attrs, EnumAttribute("alwaysinline"; ctx))
    end

    # override triple and datalayout to avoid warnings
    triple!(lib, triple(mod))
    datalayout!(lib, datalayout(mod))

    LLVM.link!(mod, lib)
end

function locate_lib(file)
    file_path = joinpath(device_libs_path, file*".bc")
    if !ispath(file_path)
        file_path = joinpath(device_libs_path, file*".amdgcn.bc")
        if !ispath(file_path)
            # failed to find matching bitcode file
            return nothing
        end
    end
    return file_path
end

function link_device_libs!(target, mod::LLVM.Module)
    # TODO: only link if used
    # TODO: make these globally/locally configurable

    device_libs_path === nothing && return

    # https://github.com/RadeonOpenCompute/ROCm-Device-Libs/blob/9420f6380990b09851edc2a5f9cbfaa88742b449/doc/OCML.md#controls
    # Note: It seems we need to load in reverse order, to avoid LLVM deleting the globals from the module, before we use them.

    # 1. Load other libraries
    libs = (
        "hc",
        "hip",
        "irif",
        "ockl",
        "opencl",
        "ocml",
    )

    for lib in libs
        lib_path = locate_lib(lib)
        lib_path === nothing && continue
        try
            load_and_link!(mod, lib_path)
        catch err
            @warn "Failed to load/link $lib" err=err
        end
    end

    # 2. Load OCLC library
    isa_short = replace(target.dev_isa, "gfx"=>"")
    lib = locate_lib("oclc_isa_version_$isa_short")
    @assert lib !== nothing
    try
        load_and_link!(mod, lib)
    catch err
        @warn "Failed to load/link OCLC core library for ISA $(target.dev_isa)" err=err
    end

    # 3. Load options libraries
    options = Dict(
        :finite_only => false,
        :unsafe_math => false,
        :correctly_rounded_sqrt => true,
        :daz_opt => false,
        :wavefrontsize64 => true
    )
    for (option, value) in options
        toggle = value ? "on" : "off"
        lib = locate_lib("oclc_$(option)_$(toggle)")
        if lib === nothing
            @warn "Could not find OCLC library for option $option=$value"
            continue
        end
        try
            load_and_link!(mod, lib)
        catch err
            @warn "Failed to load/link OCLC library for option $option=$value" err=err
        end
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Clang_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Clang")
JLLWrappers.@generate_main_file("Clang", UUID("0ee61d77-7f21-5576-8119-9fcc46b10100"))
end  # module Clang_jll


text:
__precompile__()

module JFVM

# global mumps_solver
# using PyPlot
try
  import MUMPS
  global mumps_solver = MUMPS
catch
  @info "MUMPS solver (optional) is not available."
end

using SparseArrays, FFTW
# using PyCall
# I prefer not to use the following command for the issues that it has on windows machines
# pygui_start(:wx)
# mayavis=0
# try
#   @pyimport mayavi.mlab as m
#   mayavis=m
# catch
#   warn("Mayavi is not installed or could not be imported.")
# end

import Base: +, -, *, /, ^, ==, >, >=, <, <=, broadcast, sin, cos, tan, cot, abs, exp, log, log10
export MeshStructure, BoundaryCondition, CellValue, FaceValue, CellVector,
       arithmeticMean, geometricMean, harmonicMean, upwindMean, linearMean,
       tvdMean, createBC, boundaryConditionTerm, cellBoundary!, solvePDE,
       divergenceTerm, gradientTerm, convectionUpwindTerm, createCellVector,
       convectionTerm, convectionTvdTerm, diffusionTerm, createCellVariable,
       createFaceVariable, copyCell, fluxLimiter, createMesh1D,
       createMesh2D, createMesh3D, createMeshRadial2D, createMeshCylindrical2D,
       createMeshCylindrical3D, createMeshCylindrical1D, solveLinearPDE,
       linearSourceTerm, constantSourceTerm, transientTerm,
       solveMUMPSLinearPDE, faceEval, cellEval, permfieldlogrndg, permfieldlogrnde,
       JFVM_test, solveExplicitPDE, reshapeCell,
       cellVolume, reshapeInternalCell, internalCells, domainInt, convectionTvdRHS,
       linearMean!, update!, solveLinearPDE!

      #  visualizeCells, visualizeCellVectors, plot, imshow, xlabel, ylabel, figure, legend, pcolor, contour, colorbar,

include("fvmToolTypes.jl")
include("meshstructure.jl")
include("boundarycondition.jl")
include("domainVariables.jl")
include("diffusionterms.jl")
include("transientTerms.jl")
include("domainOperators.jl")
include("convectionTerms.jl")
include("averagingTerms.jl")
include("calculusTerms.jl")
include("sourceTerms.jl")
include("solveVisualizePDE.jl")
include("JFVMtools.jl")
include("jfvm_test.jl")

end # module

max_stars_repo_name: simulkade/JFVM

max_forks_repo_name: simulkade/JFVM


text:
for n = 1 : 10
    print(n)
end
println()

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using POMDPModels
using POMDPs
using POMDPModelTools
using StaticArrays
using Random
using DiscreteValueIteration
using GridInterpolations
using LocalFunctionApproximation
using LocalApproximationValueIteration
using Test

@testset "all" begin

    println("Testing Requirements")
    gifa = LocalGIFunctionApproximator(RectangleGrid([0.0, 1.0], [0.0, 1.0]))
    @requirements_info LocalApproximationValueIterationSolver(gifa) SimpleGridWorld()

    @testset "integration" begin
        include("runtests_versus_discrete_vi.jl")
    end
end


text:
using Distributions, StatsBase
using MixtureModels

function test_emstep(::Type{Geometric}, N=1000)
    p = rand()
    samples = rand(Geometric(p), N)
    weights = uweights(N)

    samples_hist = unique(samples)
    weights_hist = FrequencyWeights([ count(x -> x == s, samples) for s in samples_hist ])

    pred = emstep(Geometric(), samples, weights, ones(N))
    pred_hist = emstep(Geometric(), samples_hist, weights_hist, ones(N))

    isapprox(pred.p, p, atol=0.05) && isapprox(pred_hist.p, p, atol=0.05)
end

function test_emstep(::Type{Poisson}, N=1000)
    λ = exp(randn())
    samples = rand(Poisson(λ), N)
    weights = uweights(N)

    samples_hist = unique(samples)
    weights_hist = FrequencyWeights([ count(x -> x == s, samples) for s in samples_hist ])

    pred = emstep(Poisson(), samples, weights, ones(N))
    pred_hist = emstep(Poisson(), samples_hist, weights_hist, ones(N))

    isapprox(pred.λ, λ, rtol=0.05) && isapprox(pred_hist.λ, λ, rtol=0.05)
end

function test_emstep(::Type{Bernoulli}, N=1000)
    N = 1000
    p = rand()
    samples = rand(Bernoulli(p), N)
    weights = uweights(N)

    samples_hist = unique(samples)
    weights_hist = FrequencyWeights([ count(x -> x == s, samples) for s in samples_hist ])

    pred = emstep(Bernoulli(), samples, weights, ones(N))
    pred_hist = emstep(Bernoulli(), samples_hist, weights_hist, ones(N))

    isapprox(pred.p, p, atol=0.05) && isapprox(pred.p, p, atol=0.05)
end

function test_emstep(::Type{Binomial}, n, N=1000)
    p = rand()
    samples = rand(Binomial(n, p), N)
    weights = uweights(N)

    samples_hist = unique(samples)
    weights_hist = FrequencyWeights([ count(x -> x == s, samples) for s in samples_hist ])

    pred = emstep(Binomial(n), samples, weights, ones(N))
    pred_hist = emstep(Binomial(n), samples_hist, weights_hist, ones(N))

    pred.n == n && isapprox(pred.p, p, atol=0.05) && pred_hist.n == n && isapprox(pred_hist.p, p, atol=0.05)
end

function test_emstep(::Type{NegativeBinomial}, N=10000, max_iter=100)
    r = exp(randn())
    p = rand()
    samples = rand(NegativeBinomial(r, p), N)
    weights = uweights(N)

    samples_hist = unique(samples)
    weights_hist = FrequencyWeights([ count(x -> x == s, samples) for s in samples_hist ])

    pred = pred_hist = NegativeBinomial(1, 1)
    for i in 1:max_iter
        pred = emstep(pred, samples, weights, ones(N))
        pred_hist = emstep(pred, samples_hist, weights_hist, ones(N))
    end

    isapprox(pred.r, r, atol=0.1) && isapprox(pred.p, p, atol=0.05) &&  isapprox(pred_hist.r, r, atol=0.1) && isapprox(pred_hist.p, p, atol=0.05)
end

N_rep = 100
p_success = 0.8

for T in [ Geometric, Poisson, Bernoulli ]
    @test count(test_emstep(T) for i in 1:N_rep) >= p_success * N_rep
end

@test count(test_emstep(Binomial, 5) for i in 1:N_rep) >= p_success * N_rep
@test count(test_emstep(Binomial, 20) for i in 1:N_rep) >= p_success * N_rep

@test count(test_emstep(NegativeBinomial) for i in 1:N_rep) >= 0.7 * N_rep


text:
using OptimalTransport

using Distances
using ForwardDiff
using ReverseDiff
using LogExpFunctions
using PythonOT: PythonOT

using LinearAlgebra
using Random
using Test

const POT = PythonOT

Random.seed!(100)

@testset "sinkhorn_gibbs.jl" begin
    # size of source and target
    M = 250
    N = 200

    # create two random histograms
    μ = normalize!(rand(M), 1)
    ν = normalize!(rand(N), 1)

    # create random cost matrix
    C = pairwise(SqEuclidean(), rand(1, M), rand(1, N); dims=2)

    # regularization parameter
    ε = 0.01

    @testset "example" begin
        # compute optimal transport plan and optimal transport cost
        γ = sinkhorn(μ, ν, C, ε, SinkhornGibbs(); maxiter=5_000, rtol=1e-9)
        c = sinkhorn2(μ, ν, C, ε, SinkhornGibbs(); maxiter=5_000, rtol=1e-9)

        # check that plan and cost are consistent
        @test c ≈ dot(γ, C)

        # compare with default algorithm
        γ_default = sinkhorn(μ, ν, C, ε; maxiter=5_000, rtol=1e-9)
        c_default = sinkhorn2(μ, ν, C, ε; maxiter=5_000, rtol=1e-9)
        @test γ_default == γ
        @test c_default == c

        # compare with POT
        γ_pot = POT.sinkhorn(μ, ν, C, ε; numItermax=5_000, stopThr=1e-9)
        c_pot = POT.sinkhorn2(μ, ν, C, ε; numItermax=5_000, stopThr=1e-9)[1]
        @test γ_pot ≈ γ rtol = 1e-6
        @test c_pot ≈ c rtol = 1e-7

        # compute optimal transport cost with regularization term
        c_w_regularization = sinkhorn2(
            μ, ν, C, ε, SinkhornGibbs(); maxiter=5_000, regularization=true
        )
        @test c_w_regularization ≈ c + ε * sum(x -> iszero(x) ? x : x * log(x), γ)
        @test c_w_regularization ==
              sinkhorn2(μ, ν, C, ε; maxiter=5_000, regularization=true)

        # ensure that provided plan is used and correct
        c2 = sinkhorn2(similar(μ), similar(ν), C, rand(), SinkhornGibbs(); plan=γ)
        @test c2 ≈ c
        @test c2 == sinkhorn2(similar(μ), similar(ν), C, rand(); plan=γ)
        c2_w_regularization = sinkhorn2(
            similar(μ), similar(ν), C, ε, SinkhornGibbs(); plan=γ, regularization=true
        )
        @test c2_w_regularization ≈ c_w_regularization
        @test c2_w_regularization ==
              sinkhorn2(similar(μ), similar(ν), C, ε; plan=γ, regularization=true)

        # batches of histograms
        d = 10
        for (size2_μ, size2_ν) in
            (((), (d,)), ((1,), (d,)), ((d,), ()), ((d,), (1,)), ((d,), (d,)))
            # generate batches of histograms
            μ_batch = repeat(μ, 1, size2_μ...)
            ν_batch = repeat(ν, 1, size2_ν...)

            # compute optimal transport plan and check that it is consistent with the
            # plan for individual histograms
            γ_all = sinkhorn(
                μ_batch, ν_batch, C, ε, SinkhornGibbs(); maxiter=5_000, rtol=1e-9
            )
            @test size(γ_all) == (M, N, d)
            @test all(view(γ_all, :, :, i) ≈ γ for i in axes(γ_all, 3))
            @test γ_all == sinkhorn(μ_batch, ν_batch, C, ε; maxiter=5_000, rtol=1e-9)

            # compute optimal transport cost and check that it is consistent with the
            # cost for individual histograms
            c_all = sinkhorn2(
                μ_batch, ν_batch, C, ε, SinkhornGibbs(); maxiter=5_000, rtol=1e-9
            )
            @test size(c_all) == (d,)
            @test all(x ≈ c for x in c_all)
            @test c_all == sinkhorn2(μ_batch, ν_batch, C, ε; maxiter=5_000, rtol=1e-9)
        end
    end

    # different element type
    @testset "Float32" begin
        # create histograms and cost matrix with element type `Float32`
        μ32 = map(Float32, μ)
        ν32 = map(Float32, ν)
        C32 = map(Float32, C)
        ε32 = Float32(ε)

        # compute optimal transport plan and optimal transport cost
        γ = sinkhorn(μ32, ν32, C32, ε32, SinkhornGibbs(); maxiter=5_000, rtol=1e-6)
        c = sinkhorn2(μ32, ν32, C32, ε32, SinkhornGibbs(); maxiter=5_000, rtol=1e-6)
        @test eltype(γ) === Float32
        @test typeof(c) === Float32

        # check that plan and cost are consistent
        @test c ≈ dot(γ, C32)

        # compare with default algorithm
        γ_default = sinkhorn(μ32, ν32, C32, ε32; maxiter=5_000, rtol=1e-6)
        c_default = sinkhorn2(μ32, ν32, C32, ε32; maxiter=5_000, rtol=1e-6)
        @test γ_default == γ
        @test c_default == c

        # compare with POT
        γ_pot = POT.sinkhorn(μ32, ν32, C32, ε32; numItermax=5_000, stopThr=1e-6)
        c_pot = POT.sinkhorn2(μ32, ν32, C32, ε32; numItermax=5_000, stopThr=1e-6)[1]
        @test map(Float32, γ_pot) ≈ γ rtol = 1e-3
        @test Float32(c_pot) ≈ c rtol = 1e-3

        # batches of histograms
        d = 10
        for (size2_μ, size2_ν) in
            (((), (d,)), ((1,), (d,)), ((d,), ()), ((d,), (1,)), ((d,), (d,)))
            # generate batches of histograms
            μ32_batch = repeat(μ32, 1, size2_μ...)
            ν32_batch = repeat(ν32, 1, size2_ν...)

            # compute optimal transport plan and check that it is consistent with the
            # plan for individual histograms
            γ_all = sinkhorn(
                μ32_batch, ν32_batch, C32, ε32, SinkhornGibbs(); maxiter=5_000, rtol=1e-6
            )
            @test size(γ_all) == (M, N, d)
            @test all(view(γ_all, :, :, i) ≈ γ for i in axes(γ_all, 3))
            @test γ_all ==
                  sinkhorn(μ32_batch, ν32_batch, C32, ε32; maxiter=5_000, rtol=1e-6)

            # compute optimal transport cost and check that it is consistent with the
            # cost for individual histograms
            c_all = sinkhorn2(
                μ32_batch, ν32_batch, C32, ε32, SinkhornGibbs(); maxiter=5_000, rtol=1e-6
            )
            @test size(c_all) == (d,)
            @test all(x ≈ c for x in c_all)
            @test c_all ==
                  sinkhorn2(μ32_batch, ν32_batch, C32, ε32; maxiter=5_000, rtol=1e-6)
        end
    end

    # https://github.com/JuliaOptimalTransport/OptimalTransport.jl/issues/86
    @testset "AD" begin
        # compute gradients with respect to source and target marginals separately and
        # together. test against gradient computed using analytic formula of Proposition 2.3 of
        # Cuturi, Marco, and Gabriel Peyré. "A smoothed dual approach for variational Wasserstein problems." SIAM Journal on Imaging Sciences 9.1 (2016): 320-343.
        #
        ε = 0.05 # use a larger ε to avoid having to do many iterations
        # target marginal
        for Diff in [ReverseDiff, ForwardDiff]
            ∇ = Diff.gradient(log.(ν)) do xs
                sinkhorn2(μ, softmax(xs), C, ε, SinkhornGibbs(); regularization=true)
            end
            ∇default = Diff.gradient(log.(ν)) do xs
                sinkhorn2(μ, softmax(xs), C, ε; regularization=true)
            end
            @test ∇ == ∇default

            solver = OptimalTransport.build_solver(μ, ν, C, ε, SinkhornGibbs())
            OptimalTransport.solve!(solver)
            # helper function
            function dualvar_to_grad(x, ε)
                x = -ε * log.(x)
                x .-= sum(x) / size(x, 1)
                return -x
            end
            ∇_ot = dualvar_to_grad(solver.cache.v, ε)
            # chain rule because target measure parameterised by softmax
            J_softmax = ForwardDiff.jacobian(log.(ν)) do xs
                softmax(xs)
            end
            ∇analytic_target = J_softmax * ∇_ot
            # check that gradient obtained by AD matches the analytic formula
            @test ∇ ≈ ∇analytic_target rtol = 1e-6

            # source marginal
            ∇ = Diff.gradient(log.(μ)) do xs
                sinkhorn2(softmax(xs), ν, C, ε, SinkhornGibbs(); regularization=true)
            end
            ∇default = Diff.gradient(log.(μ)) do xs
                sinkhorn2(softmax(xs), ν, C, ε; regularization=true)
            end
            @test ∇ == ∇default

            # check that gradient obtained by AD matches the analytic formula
            solver = OptimalTransport.build_solver(μ, ν, C, ε, SinkhornGibbs())
            OptimalTransport.solve!(solver)
            J_softmax = ForwardDiff.jacobian(log.(μ)) do xs
                softmax(xs)
            end
            ∇_ot = dualvar_to_grad(solver.cache.u, ε)
            ∇analytic_source = J_softmax * ∇_ot
            @test ∇ ≈ ∇analytic_source rtol = 1e-6

            # both marginals
            ∇ = Diff.gradient(log.(vcat(μ, ν))) do xs
                sinkhorn2(
                    softmax(xs[1:M]),
                    softmax(xs[(M + 1):end]),
                    C,
                    ε,
                    SinkhornGibbs();
                    regularization=true,
                )
            end
            ∇default = Diff.gradient(log.(vcat(μ, ν))) do xs
                sinkhorn2(softmax(xs[1:M]), softmax(xs[(M + 1):end]), C, ε; regularization=true)
            end
            @test ∇ == ∇default
            ∇analytic = vcat(∇analytic_source, ∇analytic_target)
            @test ∇ ≈ ∇analytic rtol = 1e-6
        end
    end

    @testset "deprecations" begin
        # check `sinkhorn`
        γ = sinkhorn(μ, ν, C, ε; atol=1e-6)
        @test (@test_deprecated sinkhorn(μ, ν, C, ε; tol=1e-6)) == γ
        γ = sinkhorn(μ, ν, C, ε; check_convergence=5)
        @test (@test_deprecated sinkhorn(μ, ν, C, ε; check_marginal_step=5)) == γ

        # check `sinkhorn2`
        c = sinkhorn2(μ, ν, C, ε; atol=1e-6)
        @test (@test_deprecated sinkhorn2(μ, ν, C, ε; tol=1e-6)) == c
        c = sinkhorn2(μ, ν, C, ε; check_convergence=5)
        @test (@test_deprecated sinkhorn2(μ, ν, C, ε; check_marginal_step=5)) == c
    end
end


text:
module ArbitrarySampler
export get_sample, @spec, generate
using MLStyle
import InteractiveUtils: subtypes
# for typed test generator

struct Iterator{T}
    f :: Function
end

function next_item!(iter::Iterator{T}):: Union{Nothing, Some{T}} where T
    e  = iter.f() :: Union{Nothing, Some{<:T}}
    e === nothing && return nothing
    e.value isa T && return e
    Some{T}(e.value)
end

function take_all!(iter::Iterator{T})::Vector{T} where T
    xs = T[]
    while (a = next_item!(iter); a !== nothing)
        push!(xs, a.value)
    end
    return xs
end

function take_n!(iter::Iterator{T}, n :: Int)::Vector{T} where T
    xs = T[]
    while (a = next_item!(iter); a !== nothing && length(xs) < n)
        push!(xs, a.value)
    end
    return xs
end

abstract type SeqElt{T}  end
abstract type Spec{T} end

→(A, B) = Function

@data SeqElt{T} begin
    One(Spec{T})
    Ellipsis(Spec{Vector{T}})
end

@data Spec{T} begin
    S_of{T}           :: (() → T)           => Spec{T}
    S_modif{T}        :: ((T → T), Spec{T}) => Spec{T}
    S_rep_of{T}       :: Spec{T}            => Spec{Iterator{T}}
    S_app{F}          :: (f, args)          => Spec{F}
    S_tp{Tp <: Tuple} :: tup                => Spec{Tp} where
      { Tp = Tuple{map(get_spec_type, tup)...} }

    S_guard{T}  :: (filter::(T → Bool), spec::Spec{T}) => Spec{T}
    S_map{A, B} :: (f::(A → B), spec::Spec{A})         => Spec{B}
    S_or{A, B}  :: (Spec{A}, Spec{B})                  => Spec{Union{A, B}}

    S_vec_of{T} :: (xs::Vector{SeqElt})                => Spec{Vector{T}} where
      { T =  isempty(xs) ? Any : reduce(typejoin, map(get_elt_type, xs)) }
end

S_tp(hd, args...) = S_tp((hd, args...))
get_spec_type(::Spec{T}) where T = T
get_elt_type(::SeqElt{T}) where T = T

function int_mod(n::Int, spec::Spec{Int})
    function modify(orig::Int)::Int
        orig % n
    end
    S_modif(modify, spec)
end

@nospecialize
function _take_method_arg_type(t)
    while t isa UnionAll
        t = t.body
    end
    t = t.parameters[2].parameters[1]
    while t isa TypeVar
        t = t.ub
    end
    return t
end
@specialize

get_subtypes(t) = get_subtypes(Val(Base.get_world_counter()), t)
@generated function get_subtypes(::Val{WorldAge}, ::Type{T}) where {WorldAge, T}
    subtypes(T)
end

get_registered_types() = get_registered_types(Val(Base.get_world_counter()))
@generated function get_registered_types(::Val{WorldAge}) where WorldAge
    meths = methods(get_sample)
    Any[_take_method_arg_type(m.sig) for m in meths.ms]
end

get_registered_datatypes() = get_registered_datatypes(Val(Base.get_world_counter()))
@generated function get_registered_datatypes(::Val{WorldAge}) where WorldAge
    regs = get_registered_types(Val(WorldAge))
    DataType[t for t in regs if t isa DataType]
end

function _get_all_concrete!(sup, cts::Set)
    if isabstracttype(sup)
        if Core.Compiler.unwrap_unionall(sup).name.module === Core
            push!(cts, sup)
        end
        for sub in get_subtypes(sup)
            _get_all_concrete!(sub, cts)
        end
    else
        if Core.Compiler.unwrap_unionall(sup).name.module === Core
            push!(cts, sup)
        end
    end
end
const Cores = Set{Any}([BigFloat, BigInt])
_get_all_concrete!(Real, Cores)

rand′(x) = rand(x)
rand′(::Type{BigInt}) = BigInt(rand(Int))
rand′(::Type{<:Irrational}) = pi
rand′(::Type{<:Irrational}) = pi
rand′(::Type{Complex}) =
    let T = rand([Int, Float64])
        rand(Complex{T})
    end
rand′(::Type{Rational}) =
    let T = rand([Integer])
        get_sample(T) // (abs(get_sample(T)) + 1)
    end

function _get_sample_generic(x)
    regs = get_registered_types()
    if isabstracttype(x)
        possible = [t for t in get_subtypes(x) if t in Cores || t in regs]
        isempty(possible) && return nothing
        return get_sample(rand(possible))
    else
        (x in Cores || x in regs) && return rand′(x)
        return nothing
    end
end

get_sample(::Type{Complex}) = rand′(Complex)
function get_sample(::Type{T}) where T <: Real
    ret = _get_sample_generic(T)
    ret === nothing && error("cannot find random generator for $T.")
    return ret
end
get_sample(::Type{String}) =
    String([rand('0':'z')  for i in 1:get_sample(Int) % 100])
get_sample(::Type{Symbol}) =
    Symbol(String([rand('0':'z')  for i in 1:get_sample(Int) % 7]))

function get_sample(::Type{DataType})
    dts = get_registered_datatypes()
    dts[rand(1:length(dts))]
end

function get_sample(::Type{Any})
    get_sample(rand(get_registered_types()))
end

function generate(spec::Spec{T})::T where T
    @switch spec begin
        @case S_of(gen)
            return gen()

        @case S_or(a, b)
            return rand(Bool) === true ?
                generate(a) :
                generate(b)

        @case S_modif(modif, spec)
            return modif(generate(spec))

        @case S_rep_of{A}(unit) where A
            # infinite generator
            return Iterator{A}() do
                return Some(generate(unit))
            end

        @case S_app(f, args)
            return f(generate(args)...)

        @case S_tp(tup)
            tp = Any[]
            for e in tup
                push!(tp, generate(e))
            end
            return Tuple(tp)

        @case S_guard(filt, spec)
            a = nothing
            while (a = generate(spec); !filt(a))
            end
            return a

        @case S_map(f, spec)
            a = generate(spec)
            return f(a)

        @case S_vec_of{G}(elts) where G
            ret = G[]
            for each in elts
                @switch each begin
                    @case One(spec)
                        push!(ret, generate(spec))
                        continue
                    @case Ellipsis(spec)
                        append!(ret, generate(spec))
                end
            end
            return ret
    end
end

const any_spec = S_of{Any}() do
    get_sample(Any)
end

const sym_spec = S_of{Symbol}() do
    get_sample(Symbol)
end

function gen_spec(s, __module__::Module)
    S_of{typeof(s)}() do
        s
    end
end

function gen_spec(s::QuoteNode, __module__::Module)
    s.value isa Symbol && return let s=s.value
            S_of{Symbol}() do
                s
            end
    end

    gen_spec(s.value)
end

function gen_spec(s::Symbol, __module__::Module)
    if s === :_
        return any_spec
    end
    error("Only '_' is a supported Symbol for spec syntax so far.")
end

function qt2ex(ex::Any)
    if ex isa Expr
        Meta.isexpr(ex, :$) && return ex.args[1]
        ret = Expr(:call)
        push!(ret.args, Expr, QuoteNode(ex.head))
        for each in ex.args
            if each isa LineNumberNode
                continue
            end
            push!(ret.args, qt2ex(each))
        end
        ret
    elseif ex isa Symbol
        QuoteNode(ex)
    else
        ex
    end
end

function gen_spec(ex::Expr, __module__::Module)
    rec(ex) = gen_spec(ex, __module__)
    @switch ex begin
        @case Expr(:quote, a)
            return rec(qt2ex(a))
        @case Expr(:(::), val, ty)
            T = __module__.eval(ty)
            val = T(val)
            return S_of{T}() do
                val
            end
        @case Expr(:(::), ty)
            T = __module__.eval(ty)
            return S_of{T}() do
                get_sample(T)
            end
        @case Expr(:$, a)
            spec = __module__.eval(a)
            spec isa Spec || error("invalid spec")
            return spec
        @case Expr(:||, a, b)
            return S_or(rec(a), rec(b))
        @case Expr(:tuple, args...)
            return S_tp(Tuple(rec(arg) for arg in args))
        @case Expr(:vect, args...)
            specs = SeqElt[]
            for arg in args
                @switch arg begin
                    @case :($a...)
                        push!(specs, Ellipsis(rec(a)))
                        continue
                    @case _
                        push!(specs, One(rec(arg)))
                end
            end
            return S_vec_of(specs)
        @case Expr(:comparison, args...)
            isodd(length(args)) || error("invalid spec: $ex")
            args[2] === :isa || error("invalid spec: $ex, require 'a isa b isa c'")
            return foldl(3:2:length(args), init=rec(args[1])) do last, i
                    arg = args[i]
                    args[i-1] === :isa || error("invalid spec: $ex, require 'a isa b isa c'")
                    @match arg begin
                        :($f.?) => S_guard(__module__.eval(f), last)
                        _ => S_modif(__module__.eval(arg), last)
                    end
            end
        @case :($a isa $f)
            return @match f begin
                :($f.?) => S_guard(__module__.eval(f), rec(a))
                _ => S_modif(__module__.eval(f), rec(a))
            end

        @case :($a.?($f))
            return S_guard(__module__.eval(f), rec(a))

        @case :($a{$min_v, $max_v}) || :($a{$min_v}) && let max_v = min_v end
            min_v = __module__.eval(min_v)::Int
            max_v = __module__.eval(max_v)::Int
            spec = rec(a)
            T = get_spec_type(spec)
            rep_spec = S_rep_of(spec)
            if min_v === max_v
                return  S_map{Iterator{T}, Vector{T}}(rep_spec) do iter
                    take_n!(iter, min_v)
                end
            else
                return  S_map{Iterator{T}, Vector{T}}(rep_spec) do iter
                    take_n!(iter, rand(min_v:max_v))
                end
            end
        @case Expr(:call, f, args...)
            f = __module__.eval(f)
            return S_app{f}(f, rec(Expr(:vect, args...)))
        @case _
            error("not supported syntax: $ex")
    end
end

macro spec(ex)
    gen_spec(ex, __module__)
end

# mod′(n) = x -> x % n
# len_lt(n) = x -> length(x) < n && length(x) !== 0
# not_num(x) = !(x isa Number)
# struct S{T}
#     a :: T
#     v :: Int
# end

# spec0 = @spec (::Int) isa iseven.?
# spec1 = @spec [_, (2||3, _, 0)]
# spec2 = @spec S(::Symbol, ::Int isa mod′(12))
# spec3 = @spec 1 || ::DataType
# spec4 = @spec [1, (2, ::String isa len_lt(6).?), 3]
# spec5 = @spec (_, S(::Bool, 2), ::Int isa mod′(3))

# for (id, spec) in [0=>spec0, 1=>spec1, 2=>spec2, 3=>spec3, 4=>spec4, 5=>spec5]
#     for nth in 1:2
#         sample = generate(spec)
#         println("spec$id, ($nth)th sample: $sample")
#     end
# end
end