text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -29.23432407510253 <= q <= 25.005662503891426)

                     add_NL_constraint(m, :(gelu(-0.9584800931585375 + 0.8291526833328144*gelu(0.45833851557375827 + -0.07893784992449904*gelu(-0.22994750449881307 + -0.3590609593879259*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.05269389376863742*gelu(0.19384674006188307 + -0.9239566552941776*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.9157122134504139*gelu(0.4065836733363821 + 0.7200826994240335*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.14234137821989812*gelu(-0.18406387526761447 + 0.28290146144347617*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + -0.47308336785041494*gelu(0.8182392255084516 + 0.5838043616725095*gelu(-0.22994750449881307 + -0.3590609593879259*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.4648893036822681*gelu(0.19384674006188307 + -0.9239566552941776*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.819336300807536*gelu(0.4065836733363821 + 0.7200826994240335*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.1444146117738878*gelu(-0.18406387526761447 + 0.28290146144347617*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + 0.6108443940953072*gelu(-0.7015763069538381 + -0.6574215519391942*gelu(-0.22994750449881307 + -0.3590609593879259*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.13767395548269734*gelu(0.19384674006188307 + -0.9239566552941776*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.21605537563661414*gelu(0.4065836733363821 + 0.7200826994240335*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.3991258949611094*gelu(-0.18406387526761447 + 0.28290146144347617*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + -0.4993930565410927*gelu(-0.21931284624450997 + 0.9567749529797225*gelu(-0.22994750449881307 + -0.3590609593879259*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.8630566689068395*gelu(0.19384674006188307 + -0.9239566552941776*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.8170122079846647*gelu(0.4065836733363821 + 0.7200826994240335*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.11471605851894884*gelu(-0.18406387526761447 + 0.28290146144347617*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))) + gelu(0.475745519197857 + -0.20728788527253217*gelu(0.45833851557375827 + -0.07893784992449904*gelu(-0.22994750449881307 + -0.3590609593879259*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.05269389376863742*gelu(0.19384674006188307 + -0.9239566552941776*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.9157122134504139*gelu(0.4065836733363821 + 0.7200826994240335*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.14234137821989812*gelu(-0.18406387526761447 + 0.28290146144347617*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + 0.5536591664216424*gelu(0.8182392255084516 + 0.5838043616725095*gelu(-0.22994750449881307 + -0.3590609593879259*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.4648893036822681*gelu(0.19384674006188307 + -0.9239566552941776*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.819336300807536*gelu(0.4065836733363821 + 0.7200826994240335*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.1444146117738878*gelu(-0.18406387526761447 + 0.28290146144347617*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + 0.5772469178059145*gelu(-0.7015763069538381 + -0.6574215519391942*gelu(-0.22994750449881307 + -0.3590609593879259*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.13767395548269734*gelu(0.19384674006188307 + -0.9239566552941776*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.21605537563661414*gelu(0.4065836733363821 + 0.7200826994240335*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.3991258949611094*gelu(-0.18406387526761447 + 0.28290146144347617*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + 0.14258534816037427*gelu(-0.21931284624450997 + 0.9567749529797225*gelu(-0.22994750449881307 + -0.3590609593879259*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.8630566689068395*gelu(0.19384674006188307 + -0.9239566552941776*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.8170122079846647*gelu(0.4065836733363821 + 0.7200826994240335*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.11471605851894884*gelu(-0.18406387526761447 + 0.28290146144347617*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))) + gelu(-0.19885854025662475 + 0.9363781836864602*gelu(0.45833851557375827 + -0.07893784992449904*gelu(-0.22994750449881307 + -0.3590609593879259*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.05269389376863742*gelu(0.19384674006188307 + -0.9239566552941776*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.9157122134504139*gelu(0.4065836733363821 + 0.7200826994240335*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.14234137821989812*gelu(-0.18406387526761447 + 0.28290146144347617*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + -0.2868856241387201*gelu(0.8182392255084516 + 0.5838043616725095*gelu(-0.22994750449881307 + -0.3590609593879259*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.4648893036822681*gelu(0.19384674006188307 + -0.9239566552941776*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.819336300807536*gelu(0.4065836733363821 + 0.7200826994240335*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.1444146117738878*gelu(-0.18406387526761447 + 0.28290146144347617*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + 0.6989952318088708*gelu(-0.7015763069538381 + -0.6574215519391942*gelu(-0.22994750449881307 + -0.3590609593879259*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.13767395548269734*gelu(0.19384674006188307 + -0.9239566552941776*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.21605537563661414*gelu(0.4065836733363821 + 0.7200826994240335*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.3991258949611094*gelu(-0.18406387526761447 + 0.28290146144347617*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + 0.2430276118759056*gelu(-0.21931284624450997 + 0.9567749529797225*gelu(-0.22994750449881307 + -0.3590609593879259*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.8630566689068395*gelu(0.19384674006188307 + -0.9239566552941776*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.8170122079846647*gelu(0.4065836733363821 + 0.7200826994240335*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.11471605851894884*gelu(-0.18406387526761447 + 0.28290146144347617*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))) + gelu(-0.4748268724301301 + -0.943475605090133*gelu(0.45833851557375827 + -0.07893784992449904*gelu(-0.22994750449881307 + -0.3590609593879259*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.05269389376863742*gelu(0.19384674006188307 + -0.9239566552941776*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.9157122134504139*gelu(0.4065836733363821 + 0.7200826994240335*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.14234137821989812*gelu(-0.18406387526761447 + 0.28290146144347617*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + -0.504486324487766*gelu(0.8182392255084516 + 0.5838043616725095*gelu(-0.22994750449881307 + -0.3590609593879259*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.4648893036822681*gelu(0.19384674006188307 + -0.9239566552941776*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.819336300807536*gelu(0.4065836733363821 + 0.7200826994240335*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.1444146117738878*gelu(-0.18406387526761447 + 0.28290146144347617*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + -0.8621266685661562*gelu(-0.7015763069538381 + -0.6574215519391942*gelu(-0.22994750449881307 + -0.3590609593879259*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.13767395548269734*gelu(0.19384674006188307 + -0.9239566552941776*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.21605537563661414*gelu(0.4065836733363821 + 0.7200826994240335*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.3991258949611094*gelu(-0.18406387526761447 + 0.28290146144347617*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + -0.8772206856506037*gelu(-0.21931284624450997 + 0.9567749529797225*gelu(-0.22994750449881307 + -0.3590609593879259*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.8630566689068395*gelu(0.19384674006188307 + -0.9239566552941776*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.8170122079846647*gelu(0.4065836733363821 + 0.7200826994240335*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.11471605851894884*gelu(-0.18406387526761447 + 0.28290146144347617*gelu(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*gelu(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*gelu(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*gelu(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))) - $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:
function approxextrema(I)
    arr = sample(I)
    return (minimum(arr), maximum(arr))
end

"""
    estimatequantile(array, p)

Quickly estimates the quantile `p` of a possibly large array by using a subset of the data.
"""
function estimatequantile(array, p)
    try
        return quantile(sample(array; n=1e5), p)
    catch
        @warn "quantile could not be estimated! (maybe only NaNs)"
        return 0
    end
end

function sample(I; n=10000)
    n = min(n, length(I))
    len = ceil(Int, √n) # take len blocks of len elements
    startindices = round.(Int, range(firstindex(I) - 1, lastindex(I) - len; length=len))
    indices = vcat((i .+ (1:len) for i in startindices)...)
    ret = filter(isfinite, I[indices])
    if isempty(ret)
        ret = filter(isfinite, I)
    end
    return ret
end

function get_corner_indices(I; max_length=10)
    d = size(I)
    n = min.(max_length, ceil.(Int, d ./ 3)) # n voxels for each dim
    getrange(num, len) = [1:num, (len-num+1):len] # first and last voxels
    return collect(Iterators.product(getrange.(n, d)...))
end

# estimate noise parameters from corner without signal
"""
    estimatenoise(image::AbstractArray)

Estimates the noise from the corners of the image.
It assumes that at least one corner is without signal and only contains noise.
"""
function estimatenoise(image::AbstractArray)
    corners = get_corner_indices(image)
    (lowestmean, ind) = findmin(mean.(filter(isfinite, image[I...]) for I in corners))
    sigma = std(filter(isfinite, image[corners[ind]...]))
    if isnan(sigma) # no corner available
        # estimation that is correct if half the image is signal and half noise
        sigma = 2estimatesigma_from_quantile(image, 1/4)
        lowestmean = sigma / 2
    end
    return lowestmean, sigma
end

# sigma is only calculated for quantile (biased)
function estimatesigma_from_quantile(image, quantile)
    q = estimatequantile(image, quantile)
    samples = filter(x -> x < q, sample(image))
    return std(samples)
end

"""
    robustmask!(image)
    robustmask!(image; maskedvalue)

Creates a mask and applies it inplace.
It assumes that at least one corner is without signal and only contains noise.
"""
function robustmask!(image; maskedvalue=if eltype(image) <: AbstractFloat NaN else 0 end)
    image[.!robustmask(image)] .= maskedvalue
    image
end
"""
    robustmask(weight::AbstractArray)

Creates a mask from a weights images.
It assumes that at least one corner is without signal and only contains noise.
"""
function robustmask(weight::AbstractArray)
    μ, σ = estimatenoise(weight)
    m = mean(filter(isfinite, weight[weight .> 5σ]))
    mask = weight .> maximum((5σ, m/5))
    # remove small holes and minimally grow
    boxsizes=[[3,3] for i in 1:ndims(weight)]
    return gaussiansmooth3d(mask; nbox=2, boxsizes) .> 0.55
end

getcomplex(fnmag::AbstractString, fnphase::AbstractString) = getcomplex(niread(fnmag), niread(fnphase))

function getcomplex(mag, phase)
    minp, maxp = approxextrema(phase)
    mag .* exp.((2im * pi / (maxp - minp)) .* phase)
end

function readfromtextheader(filename, searchstring)
    for line in readlines(open(filename, "r"))
        if occursin(searchstring, line)
            # regex to search for "= " or ": " and return the following non-whitespace characters
            return match(r"(?<=(= |: ))(\S+)", line).match
        end
    end
end

# root sum of squares combination
"""
    RSS(mag; dims=ndims(mag))

Performs root-sum-of-squares combination along the last dimension of `mag`.
The dimension can be specificed via the `dims` keyword argument.
"""
RSS(mag; dims=ndims(mag)) = dropdims(.√sum(mag.^Float32(2); dims); dims)

function getscaledimage(array, div::Number, offset = 0, type::Symbol = :trans)
    array = reshape(array, size(array)[1:2]) # drops trailing singleton dimensions
    scaled = if offset != 0
        (array .- offset) .* (1 / div) .+ 0.5
    else
        array .* (1 / div)
    end
    scaled[isnan.(scaled) .| (scaled .< 0)] .= 0
    scaled[scaled .> 1] .= 1
    if type == :trans
        scaled = reverse(permutedims(scaled, [2 1]); dims = 1)
    else
    end
    scaled
end

function getscaledimage(array, type::Symbol = :trans)
    scaled = robustrescale(array, 0, 1, threshold=true)
    getscaledimage(scaled, 1, 0, type)
end

"""
    robustrescale(array, newmin, newmax; threshold=false, mask=trues(size(array)), datatype=Float64)

Rescales the image to the the new range, disregarding outliers.
Only values inside `mask` are used for estimating the rescaling option
"""
robustrescale(array, newmin, newmax; threshold=false, mask=trues(size(array)), datatype=Float64) =
    robustrescale!(datatype.(array), newmin, newmax; threshold, mask)

function robustrescale!(array, newmin, newmax; threshold=false, mask=trues(size(array)))
    mask[isnan.(array)] .= false
    q = [0.01, 0.99] # quantiles
    oldq = estimatequantile(array[mask], q)
    oldrange = (oldq[2] - oldq[1]) / (q[2] - q[1])
    oldmin = oldq[1] - q[1] * oldrange
    newrange = newmax - newmin

    array .= (array .- oldmin) .* (newrange / oldrange) .+ newmin

    if threshold
        array[array .< newmin] .= newmin
        array[array .> newmax] .= newmax
    end
    array
end

function rescale(array, newmin, newmax; datatype=eltype(array))
    rescale!(datatype.(array), newmin, newmax)
end

function rescale!(array, newmin, newmax)
    oldmin, oldmax = approxextrema(array)
    factor = (newmax - newmin) / (oldmax - oldmin)
    array .= (array .- oldmin) .* factor .+ newmin
end

"""
    to_dim(V::AbstractVector, dim::Int)

    to_dim(a::Real, dim::Int)

Converts a vector or number to a higher dimension.

# Examples
```julia-repl
julia> to_dim(5, 3)
1×1×1 Array{Int64, 3}:
[:, :, 1] =
 5
julia> to_dim([1,2], 2)
 1×2 Matrix{Int64}:
  1  2
```
"""
to_dim(a::Real, dim::Int) = to_dim([a], dim)
to_dim(V::AbstractArray, dim::Int) = reshape(V, ones(Int, dim-1)..., :)

"""
    getHIP(mag, phase; echoes=[1,2])

    getHIP(compl; echoes=[1,2])

Calculates the Hermitian Inner Product between the specified echoes.
"""
function getHIP(mag, phase; echoes=[1,2])
    e1, e2 = echoes
    compl = zeros(ComplexF64, size(mag)[1:3])
    for iCha in 1:size(mag, 5)
        compl .+= exp.(1.0im .* (phase[:,:,:,e2,iCha] .- phase[:,:,:,e1,iCha])) .* mag[:,:,:,e1,iCha] .* mag[:,:,:,e2,iCha]
    end
    compl
end

function getHIP(compl; echoes=[1,2])
    e1, e2 = echoes
    c = zeros(eltype(compl), size(compl)[1:3])
    for iCha in 1:size(compl, 5)
        c .+=  compl[:,:,:,e2,iCha] .* conj.(compl[:,:,:,e1,iCha])
    end
    return c
end

ext: jl


text:
module VermontTerrestrialPassageTool

using ArchGDAL
using DelimitedFiles
using GeoData
using Omniscape
using Shapefile
using Statistics

include("utils.jl")
include("main.jl")
include("consts.jl")

export compute_all_culverts

end


text:
# just make sure notebooks are generated properly
include("../src/notebooks.jl")

max_stars_repo_name: backwardn/diff-zoo

max_issues_repo_name: backwardn/diff-zoo

max_forks_repo_name: backwardn/diff-zoo


text:
__precompile__()
module LabConnections

    module BeagleBone
        export run_server
        import Base: read
        println("Initializing BB")
        include(joinpath("BeagleBone","BeagleBone.jl"))
        include(joinpath("BeagleBone","precompile.jl"))
        println("Precompiling BB")
        #precompile_bb()
        return
    end

    module Computer
        import Base: read, send, close, get, put!, serialize
        println("Initializing Computer")
        include(joinpath("Computer","Computer.jl"))
    end
end


text:
using LinearMaps
using Magnus

using PyCall
pygui(:qt)
using PyPlot

include("plot_convergence.jl")

type ExpEulerExponentiator <: Magnus.Exponentiator
end

function (SE::ExpEulerExponentiator)(Ω::LinearMap{T}, τ::Number,
                                     v, w) where T<:Number
    a = Ω(v)
    w[:] = v + τ*a
end

function test_stiff(A,y0,y_exact,tmax,N,name,
                    Exps = [:scalar, :exp_euler],
                    Na = 0, Nb = 3)
    t = linspace(0,tmax,N+1)
    y_exact = y_exact(t)

    figure(name)
    clf()
    subplot(311)
    plot(t, y_exact, "r", linewidth = 2.0, label="Exact")

    Exp_objs = Dict(:scalar => (ScalarExponentiator(y0),"Exact exp"),
                    :lanczos => (LanczosExponentiator(2,y0),"Lanczos"),
                    :exp_euler => (ExpEulerExponentiator(),"Explicit Euler"))
    Exps = [Exp_objs[e] for e in Exps]

    for (Exp,label) in Exps
        prop = MidpointPropagator(t -> LinearMap(A), Exp)

        y = vec(integrate(y0, tmax, N, prop, save_intermediate = true))

        subplot(311)
        plot(t,y,label=label)

        plot_convergence(y_exact[end], y0, tmax, Na, Nb,
                         prop, [312,313], label)
    end
    subplot(311)
    margins(0,0.1)
    legend(framealpha=0.75)

    subplot(312)
    xscale("log")
    yscale("log")
    xlabel(L"N")
    ylabel("Error")
    margins(0,1)

    subplot(313)
    xscale("log")
    yscale("log")
    xlabel("Error")
    ylabel("Execution time [s]")
    margins(0,1)

    tight_layout()
end

test_stiff(-15., [1.0], t -> exp(-15.0t), 1.0, 10, "stiff")

if Pkg.installed("LambertW") != nothing
    using LambertW
    δ = 0.01
    a = 1/δ - 1
    y = t -> 1./(lambertw(a*exp(a-t)) + 1)
    # Technically not a linear operator, but, hey.
    # See http://mathworks.com/company/newsletters/articles/stiff-differential-equations.html
    fireball = LinearMap{Float64}((y,x,α,β) -> y[:] = α*(x.^2 - x.^3))
    test_stiff(fireball, [δ], y, 1./δ, 200, "Fireball, half interval",
               [:exp_euler], 0, 4)
    test_stiff(fireball, [δ], y, 2./δ, 200, "Fireball, full interval",
               [:exp_euler], 0, 4)
end


text:
#####################################################
### Project of IN480 - From Python 2.7 to Julia 0.6
### --- SIMPLEXN - serial unit tests ---
### Authors: Fabio Fatelli, Emanuele Loprevite
#####################################################

include("../src/simplexn-serial.jl")

using Base.Test

@testset "Unit Tests" begin
	@testset "larExtrude1" begin
		VOID = [Int64[]],[[0]] # the empty simplicial model
		pattern1 = [1,-1]
		model1 = larExtrude1(VOID,pattern1) # 1D
		#println(model1)
		@test typeof(model1) == Tuple{Array{Array{Int64,1},1},Array{Array{Int64,1},1}}
		@test length(model1[1]) == length(pattern1)+1 # num of vertices, no rep
		@test length(model1[2]) == length(filter(x->x>0,pattern1)) # num of 1D-simplices

		pattern2 = [1,1,1,-1]
		model2 = larExtrude1(VOID,pattern2) # 1D
		#println(model2)
		@test typeof(model2) == Tuple{Array{Array{Int64,1},1},Array{Array{Int64,1},1}}
		@test length(model2[1]) == length(pattern2)+1 # num of vertices, no rep
		@test length(model2[2]) == length(filter(x->x>0,pattern2)) # num of 1D-simplices

		pattern3 = [1,1,1,1,1,-1]
		model3 = larExtrude1(VOID,pattern3) # 1D
		#println(model3)
		@test typeof(model3) == Tuple{Array{Array{Int64,1},1},Array{Array{Int64,1},1}}
		@test length(model3[1]) == length(pattern3)+1 # num of vertices, no rep
		@test length(model3[2]) == length(filter(x->x>0,pattern3)) # num of 1D-simplices

		model1 = larExtrude1(model1,pattern1) # 2D
		#println(model1)
		@test typeof(model1) == Tuple{Array{Array{Int64,1},1},Array{Array{Int64,1},1}}
		@test length(model1[1]) == (length(pattern1)+1)^2 # num of vertices, no rep
		@test length(model1[2]) == 2*length(filter(x->x>0,pattern1))^2 # num of 2D-simplices

		model2 = larExtrude1(model2,pattern2) # 2D
		#println(model2)
		@test typeof(model2) == Tuple{Array{Array{Int64,1},1},Array{Array{Int64,1},1}}
		@test length(model2[1]) == (length(pattern2)+1)^2 # num of vertices, no rep
		@test length(model2[2]) == 2*length(filter(x->x>0,pattern2))^2 # num of 2D-simplices

		model3 = larExtrude1(model3,pattern3) # 2D
		#println(model3)
		@test typeof(model3) == Tuple{Array{Array{Int64,1},1},Array{Array{Int64,1},1}}
		@test length(model3[1]) == (length(pattern3)+1)^2 # num of vertices, no rep
		@test length(model3[2]) == 2*length(filter(x->x>0,pattern3))^2 # num of 2D-simplices

		model1 = larExtrude1(model1,pattern1) # 3D
		#println(model1)
		@test typeof(model1) == Tuple{Array{Array{Int64,1},1},Array{Array{Int64,1},1}}
		@test length(model1[1]) == (length(pattern1)+1)^3 # num of vertices, no rep
		@test length(model1[2]) == 6*length(filter(x->x>0,pattern1))^3 # num of 3D-simplices

		model2 = larExtrude1(model2,pattern2) # 3D
		#println(model2)
		@test typeof(model2) == Tuple{Array{Array{Int64,1},1},Array{Array{Int64,1},1}}
		@test length(model2[1]) == (length(pattern2)+1)^3 # num of vertices, no rep
		@test length(model2[2]) == 6*length(filter(x->x>0,pattern2))^3 # num of 3D-simplices

		model3 = larExtrude1(model3,pattern3) # 3D
		#println(model3)
		@test typeof(model3) == Tuple{Array{Array{Int64,1},1},Array{Array{Int64,1},1}}
		@test length(model3[1]) == (length(pattern3)+1)^3 # num of vertices, no rep
		@test length(model3[2]) == 6*length(filter(x->x>0,pattern3))^3 # num of 3D-simplices

		model = ([[0.0,0],[0,1]],[[1,0],[1,1]]) # 2D with float
		@test typeof(model[1]) == Array{Array{Float64,1},1}
		model = larExtrude1(model,pattern1)
		#println(model)
		@test typeof(model) == Tuple{Array{Array{Float64,1},1},Array{Array{Int64,1},1}}
	end
	@testset "larSimplexGrid1" begin
		shape = [3] # 1D
		model1 = larSimplexGrid1(shape)
		#println(model1)
		@test typeof(model1) == Tuple{Array{Array{Int64,1},1},Array{Array{Int64,1},1}}
		@test length(model1[1]) == shape[1]+1 # num of vertices, no rep
		@test length(model1[2]) == shape[1] # num of 1D-simplices

		shape = [3,5] # 2D
		model2 = larSimplexGrid1(shape)
		#println(model2)
		@test typeof(model2) == Tuple{Array{Array{Int64,1},1},Array{Array{Int64,1},1}}
		# (shape[1]+1)*(shape[2]+1)
		@test length(model2[1]) == prod(shape)+sum(shape)+1 # num of vertices, no rep
		@test length(model2[2]) == 2*prod(shape) # num of 2D-simplices

		shape = [3,5,7] # 3D
		model3 = larSimplexGrid1(shape)
		#println(model3)
		@test typeof(model3) == Tuple{Array{Array{Int64,1},1},Array{Array{Int64,1},1}}
		# (shape[1]+1)*(shape[2]+1)*(shape[3]+1); num of vertices, no rep
		@test length(model3[1]) == prod(shape)+sum(prod.(collect(combinations(shape,2))))+sum(shape)+1
		@test length(model3[2]) == 6*prod(shape) # num of 3D-simplices
	end
	@testset "larSimplexFacets" begin
		s = larSimplexGrid1([3])[2] # 1D
		#println(s)
		sOut = larSimplexFacets(s)
		#println(sOut)
		@test typeof(sOut) == Array{Array{Int64,1},1}
		d1 = sum([length(s[k]) for k in 1:length(s)])/length(s)
		d2 = sum([length(sOut[k]) for k in 1:length(sOut)])/length(sOut)
		@test d1-1 == d2 # dimension
		# length(s)*binomial(length(s[1]),(length(s[1])-1))
		@test length(sOut) <= length(s)*length(s[1]) # "<=" because no rep

		s = larSimplexGrid1([3,5])[2] # 2D
		#println(s)
		sOut = larSimplexFacets(s)
		#println(sOut)
		@test typeof(sOut) == Array{Array{Int64,1},1}
		d1 = sum([length(s[k]) for k in 1:length(s)])/length(s)
		d2 = sum([length(sOut[k]) for k in 1:length(sOut)])/length(sOut)
		@test d1-1 == d2 # dimension
		# length(s)*binomial(length(s[1]),(length(s[1])-1))
		@test length(sOut) <= length(s)*length(s[1]) # "<=" because no rep

		s = larSimplexGrid1([3,5,7])[2] # 3D
		#println(s)
		sOut = larSimplexFacets(s)
		#println(sOut)
		@test typeof(sOut) == Array{Array{Int64,1},1}
		d1 = sum([length(s[k]) for k in 1:length(s)])/length(s)
		d2 = sum([length(sOut[k]) for k in 1:length(sOut)])/length(sOut)
		@test d1-1 == d2 # dimension
		# length(s)*binomial(length(s[1]),(length(s[1])-1))
		@test length(sOut) <= length(s)*length(s[1]) # "<=" because no rep
	end
	@testset "quads2tria" begin
		modIn = ([[0,0,0],[0,1,0],[1,0,0],[1,1,0]],[[0,1,2,3]]) # 2D
		modOut = quads2tria(modIn)
		#println(modOut)
		@test typeof(modIn) == Tuple{Array{Array{Int64,1},1},Array{Array{Int64,1},1}}
		@test typeof(modOut) == Tuple{Array{Array{Float64,1},1},Array{Array{Int64,1},1}}
		edges = 2*sum([length(modIn[2][k]) for k in 1:length(modIn[2])])
		EulerChar = length(modOut[1])-edges+length(modOut[2]) # V - E + F
		@test EulerChar == 1

		modIn = ([[0,0,0],[0,1,0],[1,0,0],[1,1,0],[0.5,1.5,0]],[[0,1,2,3,4]]) # 2D
		modOut = quads2tria(modIn)
		#println(modOut)
		@test typeof(modIn) == Tuple{Array{Array{Float64,1},1},Array{Array{Int64,1},1}}
		@test typeof(modOut) == Tuple{Array{Array{Float64,1},1},Array{Array{Int64,1},1}}
		edges = 2*sum([length(modIn[2][k]) for k in 1:length(modIn[2])])
		EulerChar = length(modOut[1])-edges+length(modOut[2]) # V - E + F
		@test EulerChar == 1

		modIn = ([[0,0,0],[0,1,0],[1,0,0],[1,1,0],[0,0,1],[0,1,1],[1,0,1],[1,1,1]],[[0,1,2,3,4,5,6,7]]) # 3D
		modOut = quads2tria(modIn)
		#println(modOut)
		@test typeof(modIn) == Tuple{Array{Array{Int64,1},1},Array{Array{Int64,1},1}}
		@test typeof(modOut) == Tuple{Array{Array{Float64,1},1},Array{Array{Int64,1},1}}
		edges = 2*sum([length(modIn[2][k]) for k in 1:length(modIn[2])])
		EulerChar = length(modOut[1])-edges+length(modOut[2]) # V - E + F
		@test EulerChar == 1

		modIn = ([[0,0,0],[0,1,0],[1,0,0],[1,1,0],[0,0,1],[0,1,1],[1,0,1],[1,1,1],[0.5,0.5,1.5]],[[0,1,2,3,4,5,6,7,8]]) # 3D
		modOut = quads2tria(modIn)
		#println(modOut)
		@test typeof(modIn) == Tuple{Array{Array{Float64,1},1},Array{Array{Int64,1},1}}
		@test typeof(modOut) == Tuple{Array{Array{Float64,1},1},Array{Array{Int64,1},1}}
		edges = 2*sum([length(modIn[2][k]) for k in 1:length(modIn[2])])
		EulerChar = length(modOut[1])-edges+length(modOut[2]) # V - E + F
		@test EulerChar == 1
	end
end

max_stars_repo_name: EmaLoprevite/Simplexn

max_issues_repo_name: EmaLoprevite/Simplexn

max_forks_repo_name: EmaLoprevite/Simplexn


text:
import LinearAlgebra: Adjoint, Transpose
import SparseArrays: SparseMatrixCSC
import Grisu: _show
export dirac

const SHORTEST = 1
const PRECISION = 3

"""
    set_properties(; kw...)

# Arguments
- precision::Int.
- islatex::Bool. true -> rendered by mathjax on IJulia
- displayall::Bool. true -> display all terms.
- numhead: Display the first part of terms. Default is 5 terms.
- newline::Bool.
- imag_suffix::String.
"""
function set_properties(; precision::Int=_precision,
                          islatex::Bool=_islatex,
                          displayall::Bool=_displayall,
                          numhead::Int=_numhead,
                          newline::Bool=_newline,
                          imag_suffix::AbstractString=_imag_suffix)
    global _precision = precision
    global _islatex = islatex
    global _displayall = displayall
    global _numhead = numhead
    global _newline = newline
    global _imag_suffix = imag_suffix
    nothing
end
const PureState = Union{Vector, Adjoint{T,Vector{T}}, Transpose{T,Vector{T}}} where T
const MixedState =  Union{Matrix, SparseMatrixCSC, Adjoint{T,Matrix{T}}, Transpose{T,Matrix{T}}} where T

function reset_properties()
    set_properties(precision=0, islatex=true, displayall=true, numhead=5, newline=false, imag_suffix="im")
end
reset_properties()

"""
    dirac(state::PureState, statename::String="ψ"; header::Bool=false)
    dirac(io::IO, state::PureState, statename::String="ψ"; header::Bool=false)
    dirac(state::PureState, dims::Vector{Int}, statename::String="ψ"; header::Bool=false)
    dirac(io::IO, state::PureState, dims::Vector{Int}, statename::String="ψ"; header::Bool=false)

Display a vector as Dirac notation.

# Example
```
julia> using Random; Random.seed!(0);

julia> dirac(randn(Complex{Float64}, 2))
|ψ⟩ = (0.480201+0.585777im)|0⟩+(-0.249614-0.0953561im)|1⟩

julia> qutrit = rand(Complex{Float64},3);

julia> dirac(qutrit, [3])
|ψ⟩ = (0.27888+0.203477im)|0⟩+(0.0423017+0.0682693im)|1⟩+(0.361828+0.973216im)|2⟩

dirac(qubitqutrit, [2,3])
|ψ⟩ = (-0.487131-0.539384im)|00⟩+(0.281063+0.573909im)|01⟩+(-0.24491-0.132634im)|02⟩+(-1.1365-1.75419im)|10⟩+(1.60954+0.155347im)|11⟩+(-0.0828287-0.42515im)|12⟩
```
"""
function dirac(io::IO, state::PureState, dims::Vector{Int}, statename::String="ψ"; header::Bool=false)
    if header
        println(io, summary(state))
    end
    if _islatex && isdefined(Main, :IJulia) && Main.IJulia.inited # for IJulia rendering
        if statename == "ψ"
            statename = "\\psi"
        end
        str = "\$" * sprint(io -> print_dirac(io, state, dims, statename)) * "\$"
        display("text/markdown", str)
    else
        print_dirac(io, state, dims, statename)
    end
end
dirac(state::PureState, dims::Vector{Int}, statename::String="ψ"; header::Bool=false) = dirac(stdout, state, dims, statename, header=header)
function dirac(io::IO, state::PureState, statename::String="ψ"; header::Bool=false)
    n = log2(length(state))
    @assert isinteger(n)
    dims = fill(2, (Int(n), ))
    dirac(io, state, dims, statename, header=header)
end
dirac(state::PureState, statename::String="ψ"; header::Bool=false) = dirac(stdout, state, statename, header=header)

"""
    dirac(state::MixedState, statename::String="ρ"; header::Bool=false)
    dirac(io::IO, state::MixedState, statename::String="ρ"; header::Bool=false)
    dirac(state::MixedState, dims::Vector{Int}, statename::String="ρ"; header::Bool=false)
    dirac(state::MixedState, ldims::Vector{Int}, rdims::Vector{Int}, statename::String="ρ"; header::Bool=false)
    dirac(io::IO, state::MixedState, ldims::Vector{Int}, rdims::Vector{Int}, statename::String="ρ"; header::Bool=false)

Display a matrix as Dirac notation.
# Example
```
julia> using Random; Random.seed!(0);

julia> op1 = randn(Complex{Float64}, 2,2);

julia> dirac(op1)
ρ = (0.480201+0.585777im)|0⟩⟨0|+(0.414801+0.210248im)|0⟩⟨1|+(-0.249614-0.0953561im)|1⟩⟨0|+(0.0459249-0.0770869im)|1⟩⟨1|

julia> op2 = randn(Complex{Float64}, 4,3);

julia> dirac(op2, [2,2], [3])
ρ = (-0.363602+1.11322im)|00⟩⟨0|+(-1.1365-1.75419im)|00⟩⟨1|+(0.197612+0.078787im)|00⟩⟨2|+(-0.487131-0.539384im)|01⟩⟨0|+(1.60954+0.155347im)|01⟩⟨1|+(-0.253062+0.334967im)|01⟩⟨2|+(0.281063+0.573909im)|10⟩⟨0|+(-0.0828287-0.42515im)|10⟩⟨1|+(0.212297-0.539294im)|10⟩⟨2|+(-0.24491-0.132634im)|11⟩⟨0|+(0.807711-0.0626612im)|11⟩⟨1|+(1.00625+0.288773im)|11⟩⟨2|
```
"""
function dirac(io::IO, state::MixedState, ldims::Vector{Int}, rdims::Vector{Int}, statename::String="ρ"; header::Bool=false)
    if header
        println(io, summary(state))
    end
    if _islatex && isdefined(Main, :IJulia) && Main.IJulia.inited # for IJulia rendering
        if statename == "ρ"
            statename = "\\hat{\\rho}"
        end
        str = "\$" * sprint(io -> print_dirac(io, state, ldims, rdims, statename)) * "\$"
        display("text/markdown", str)
    else
        print_dirac(io, state, ldims, rdims, statename)
    end
end
dirac(state::MixedState, ldims::Vector{Int}, rdims::Vector{Int}, statename::String="ρ"; header::Bool=false) = dirac(stdout, state, ldims, rdims, statename, header=header)
dirac(state::MixedState, dims::Vector{Int}, statename::String="ρ"; header::Bool=false) = dirac(stdout, state, dims, dims, statename, header=header)
function dirac(io::IO, state::MixedState, statename::String="ρ"; header::Bool=false)
    row, col = size(state)
    n, m = log2(row), log2(col)
    @assert isinteger(n) || isinteger(m)
    ldims, rdims = fill(2, (Int(n), )), fill(2, (Int(m), ))
    dirac(io, state, ldims, rdims, statename, header=header)
end
dirac(state::MixedState, statename::String="ρ"; header::Bool=false) = dirac(stdout, state, statename, header=header)

"""
    print_dirac(io::IO, state::PureState, dims::Vector{Int}, statename::String="ψ")
    print_dirac(io::IO, state::MixedState,
                        ldims::Vector{Int},
                        rdims::Vector{Int},
                        statename::String="ρ") where T <: Number
"""
function print_dirac(io::IO, state::PureState, dims::Vector{Int}, statename::String="ψ")
    io = IO_braket(io, state)
    braket = getindex(io, :braket)

    isfirstterm = true
    s = "$(braket[1])$(statename)$(braket[2]) = "
    print(io, s)
    num_term = 0
    for (idx, ent) in enumerate(state)
        num_term += 1
        iszero(ent) && continue
        if ent ≈ one(ent)
            ent = one(ent)
        end

        print_dirac_term(io, idx, ent, dims, isfirstterm)

        if isfirstterm
            isfirstterm = false
        end

        if !_displayall && num_term >= _numhead
            break
        end
        _newline && print(io, "\n", repeat(" ", length(s)-1))
    end
    num_term < length(state) && return println(io, " +...")
    !_newline && println(io)
    return nothing
end
function print_dirac(io::IO, state::MixedState,
                             ldims::Vector{Int},
                             rdims::Vector{Int},
                             statename::String="ρ") where T <: Number
    nrow, ncol = size(state)
    @assert nrow == prod(ldims) && ncol == prod(rdims)
    io = IO_braket(io, state)
    braket = getindex(io, :braket)

    isfirstterm = true
    print(io, statename, " = ")
    num_term = 0
    exitloop = false
    for row in 1:nrow
        for col in 1:ncol
            num_term += 1
            ent = state[row,col]
            iszero(ent) && continue
            if ent ≈ one(ent)
                ent = one(ent)
            end
            print_dirac_term(io, row, col, ent, ldims, rdims, isfirstterm)
            if isfirstterm
                isfirstterm = false
            end

            if !_displayall && num_term >= _numhead
                exitloop = true
                break
            end
        end
        exitloop && break
        _newline && print(io, "\n", repeat(" ", length(statename * " = ")-1))
    end

    num_term < length(state) && return println(io, " +...")
    !_newline && println(io)
    return nothing
end

function IO_braket(io::IO, state)
    if !haskey(io, :braket)
        braket = braket_str(state)
        io = IOContext(io, :braket => braket)
    end
    if !haskey(io, :compact)
        io = IOContext(io, :compact => true)
    end
    return io
end

function braket_str(state::PureState)
    if _islatex && isdefined(Main, :IJulia) && Main.IJulia.inited
        braket = ifelse(state isa Vector, ("|", "\\rangle"), ("\\langle", "|"))
    else
        braket = ifelse(state isa Vector, ("|", "⟩"), ("⟨", "|"))
    end
    return braket
end
function braket_str(state::MixedState)
    if _islatex && isdefined(Main, :IJulia) && Main.IJulia.inited
        braket = (("|", "\\rangle"), ("\\langle", "|"))
    else
        braket = (("|", "⟩"), ("⟨", "|"))
    end
    return braket
end

"""
    print_dirac_term(io::IO, idx::Int, value::T, dims::Vector{Int}, isfirstterm=false) where T
    print_dirac_term(io::IO, row::Int, col::Int,
                             value::T, ldims::Vector{Int}, rdims::Vector{Int}, isfirstterm=false) where T
"""
function print_dirac_term(io::IO, idx::Int, value::T, dims::Vector{Int}, isfirstterm=false) where T <: Number
    s = ""
    s *= sprint(x -> print_precision_value(io, value, isfirstterm)) # `x` makes no sense.

    if _islatex && isdefined(Main, :IJulia) && Main.IJulia.inited
        braket = get(io, :braket, ("|", "\\rangle"))
    else
        braket = get(io, :braket, ("|", "⟩"))
    end
    print(io, s, braket[1], ind2Nary(idx, dims), braket[2])
end
function print_dirac_term(io::IO, row::Int, col::Int,
                          value::T, ldims::Vector{Int}, rdims::Vector{Int}, isfirstterm=false) where T <: Number
    s = ""
    s *= sprint(x -> print_precision_value(io, value, isfirstterm)) # `x` makes no sense.

    if _islatex && isdefined(Main, :IJulia) && Main.IJulia.inited
        braket = get(io, :braket, (("|", "\\rangle"), ("\\langle", "|")) )
    else
        braket = get(io, :braket, (("|", "⟩"), ("⟨", "|")))
    end
    print(io, s, braket[1][1], ind2Nary(row, ldims), braket[1][2], braket[2][1], ind2Nary(col, rdims), braket[2][2])
end

# modify a part of Julia.
# https://github.com/JuliaLang/julia/blob/02aa9bbc258a380b5c9fbe2f2d3276a90a72abca/base/complex.jl#L182-L197.
function print_precision_value(io::IO, z::Complex, isfirstterm::Bool)
    r, i = reim(z)
    istwoterms = !iszero(r) && !iszero(i)
    istwoterms && !isfirstterm && print(io, "+") # both real and imaginary part are nonzero and not first term
    istwoterms && print(io, "(")
    compact = get(io, :compact, false)
    !iszero(r) && print_precision_value(io, r, isfirstterm, istwoterms, true)

    if istwoterms
        if signbit(i) && !isnan(i)
            i = -i
            print(io, compact ? "-" : " - ")
        else
            print(io, compact ? "+" : " + ")
        end
    end

    !iszero(i) && print_precision_value(io, i, isfirstterm, istwoterms, false)
    if !(isa(i,Integer) && !isa(i,Bool) || isa(i,AbstractFloat) && isfinite(i))
        print(io, "*")
    end
    !iszero(i) && print(io, _imag_suffix)
    !iszero(r) && !iszero(i) && print(io, ")")
    return nothing
end
function print_precision_value(io::IO, val::Union{Float64,Float32}, isfirstterm::Bool, istwoterms::Bool=false, ri::Bool=true)
    iszero(val) && return nothing
    if istwoterms
        if isone(val) || isone(-val)
            return print_one(io, val, ri)
        end
        return print_float(io, val)

    else
        if (isone(val) || isone(-val))
            isfirstterm && return print(io, val > 0 ? "" : "-")
            return print(io, val > 0 ? "+" : "-")
        end
        !isfirstterm && val > 0 && print(io, "+")
        print_float(io, val)
    end
end
function print_precision_value(io::IO, val::Integer, isfirstterm::Bool, istwoterms::Bool=false, ri::Bool=true)
    iszero(val) && return nothing
    if istwoterms
        if isone(val) || isone(-val)
            return print_one(io, val, ri)
        end
        return show(io, val)

    else
        if (isone(val) || isone(-val))
            isfirstterm && return print(io, val > 0 ? "" : "-")
            return print(io, val > 0 ? "+" : "-")
        end
        !isfirstterm && val > 0 && print(io, "+")
        show(io, val)
    end
end

# modify a part of Julia.
# https://github.com/JuliaLang/julia/blob/e2c98780e211c7b0887751b7eb2b7a141c2fd5e8/base/grisu/grisu.jl#L147-L153
function print_float(io::IO, val::Union{Float64,Float32})
    if !iszero(_precision)
        _show(io, val, PRECISION, _precision, val isa Float64, true)
    elseif get(io, :compact, false)
        _show(io, val, PRECISION, 6, val isa Float64, true)
    else
        _show(io, val, SHORTEST, 0, get(io, :typeinfo, Any) !== typeof(val), false)
    end
end
function print_one(io::IO, val::Number, ri::Bool)
    if ri # real part
        isone(val) && return print(io, val)
        isone(-val) && return print(io, val)
    end
end

"""
    ind2Nary(n::Int, dims::Vector{Int})

# Arguments
- `n::Int`: n th row(column) of qudit/operator array.
- `dims::Vector{Int}`: N-ary array
"""
function ind2Nary(m::Int, dims::Vector{Int})
    m = m - 1
    str = ""
    nq = length(dims)
    product = prod(dims[2:end])
    for ith in 1:nq-1
        d = div(m, product)
        m = m - d * product
        product = div(product, dims[ith+1])
        str *= string(d)
    end
    str *= string(m)
    return str
end


text:
function gaussradau(n::Integer)
    # RADAUPTS   Gauss-Legendre-Radau Quadrature Nodes and Weights
    if n == 1
        [-1.0], [2.0]
    elseif n == 2
        [-1.0, 1/3], [.5, 1.5]
    else
        # Compute via GaussJacobi:
        x, w = gaussjacobi(n - 1, 0.0, 1.0)
        @inbounds for i in 1:length(w)
            w[i] = w[i] / (1.0 + x[i])
        end
        pushfirst!(x, -1.0)
        pushfirst!(w, 2.0 / n^2)
        x, w
    end
end


text:
module PoolUtils

using ..CUDA

using Printf

export MEMDEBUG

const MEMDEBUG = ccall(:jl_is_memdebug, Bool, ())

## block of memory

export Block, INVALID, AVAILABLE, ALLOCATED, FREED, iswhole

@enum BlockState begin
    INVALID
    AVAILABLE
    ALLOCATED
    FREED
end

mutable struct Block
    buf::Mem.DeviceBuffer # base allocation
    sz::Int               # size into it
    off::Int              # offset into it

    state::BlockState
    prev::Union{Nothing,Block}
    next::Union{Nothing,Block}

    Block(buf, sz; off=0, state=INVALID, prev=nothing, next=nothing) =
        new(buf, sz, off, state, prev, next)
end

Base.sizeof(block::Block) = block.sz
Base.pointer(block::Block) = pointer(block.buf) + block.off

iswhole(block::Block) = block.prev === nothing && block.next === nothing

function Base.show(io::IO, block::Block)
    fields = [@sprintf("%p", Int(pointer(block)))]
    push!(fields, Base.format_bytes(sizeof(block)))
    push!(fields, "$(block.state)")
    block.off != 0 && push!(fields, "offset=$(block.off)")
    block.prev !== nothing && push!(fields, "prev=Block(offset=$(block.prev.off))")
    block.next !== nothing && push!(fields, "next=Block(offset=$(block.next.off))")

    print(io, "Block(", join(fields, ", "), ")")
end

end


text:
#=============================================================

ApproxCholPQ
It only implements pop, increment key, and decrement key.
All nodes with degrees 1 through n appear in their own doubly-linked lists.
Nodes of higher degrees are bundled together.

Also have nodes of degree 0, and lists now has a 1+ everywhere.

Will need to rename, as it is different from the version in Laplacians.
Maybe add a 0 to show that is allowed.

=============================================================#

function keyMap(x, n)
    return x <= n ? x : n + div(x,n)
end

function acPQ(a::Vector{Tind}) where {Tind}

    n = length(a)
    elems = Array{acPQ_elem{Tind}}(undef,n)
    lists = zeros(Tind, 2*n+2)
    minlist = one(n)

    for i in 1:length(a)
        key = a[i]
        head = lists[1+keyMap(key, n)]

        if head > zero(Tind)
            elems[i] = acPQ_elem{Tind}(zero(Tind), head, key)

            elems[head] = acPQ_elem{Tind}(i, elems[head].next, elems[head].key)
        else
            elems[i] = acPQ_elem{Tind}(zero(Tind), zero(Tind), key)

        end

        lists[1+keyMap(key, n)] = i
    end

    return acPQ(elems, lists, minlist, n, n)
end

function pqPop!(pq::acPQ{Tind}) where {Tind}
    if pq.nitems == 0
        error("ApproxPQ is empty")
    end
    while pq.lists[1+pq.minlist] == 0
        pq.minlist = pq.minlist + 1
    end
    i = pq.lists[1+pq.minlist]
    next = pq.elems[i].next

    pq.lists[1+pq.minlist] = next
    if next > 0
        pq.elems[next] = acPQ_elem(zero(Tind), pq.elems[next].next, pq.elems[next].key)
    end

    pq.nitems -= 1

    return i
end

function pqMove!(pq::acPQ{Tind}, i, newkey, oldlist, newlist) where {Tind}

    prev = pq.elems[i].prev
    next = pq.elems[i].next

    # remove i from its old list
    if next > zero(Tind)
        pq.elems[next] = acPQ_elem{Tind}(prev, pq.elems[next].next, pq.elems[next].key)
    end
    if prev > zero(Tind)
        pq.elems[prev] = acPQ_elem{Tind}(pq.elems[prev].prev, next, pq.elems[prev].key)

    else
        pq.lists[1+oldlist] = next
    end

    # insert i into its new list
    head = pq.lists[1+newlist]
    if head > 0
        pq.elems[head] = acPQ_elem{Tind}(i, pq.elems[head].next, pq.elems[head].key)
    end
    pq.lists[1+newlist] = i

    pq.elems[i] = acPQ_elem{Tind}(zero(Tind), head, newkey)

    #return Void
end

"""
    Decrement the key of element i
    This could crash if i exceeds the maxkey
"""
function pqDec!(pq::acPQ{Tind}, i) where {Tind}

    oldlist = keyMap(pq.elems[i].key, pq.n)
    newlist = keyMap(pq.elems[i].key - one(Tind), pq.n)

    if newlist != oldlist

        pqMove!(pq, i, pq.elems[i].key - one(Tind), oldlist, newlist)

        if newlist < pq.minlist
            pq.minlist = newlist
        end

    else
        pq.elems[i] = acPQ_elem{Tind}(pq.elems[i].prev, pq.elems[i].next, pq.elems[i].key - one(Tind))
    end

    #return Void
end

"""
    Increment the key of element i
    This could crash if i exceeds the maxkey
"""
function pqInc!(pq::acPQ{Tind}, i) where {Tind}

    oldlist = keyMap(pq.elems[i].key, pq.n)
    newlist = keyMap(pq.elems[i].key + one(Tind), pq.n)

    if newlist != oldlist

        pqMove!(pq, i, pq.elems[i].key + one(Tind), oldlist, newlist)

    else
        pq.elems[i] = acPQ_elem{Tind}(pq.elems[i].prev, pq.elems[i].next, pq.elems[i].key + one(Tind))
    end

    #return Void
end

function acPQ_rev(a::Vector{Tind}) where {Tind}

    #logwrite(fn,"pq init: $(a)")

    n = length(a)
    elems = Array{acPQ_elem{Tind}}(undef, n)
    heads = zeros(Tind, 2*n+2)
    tails = zeros(Tind, 2*n+2)
    minlist = one(n)

    for i in 1:length(a)
        key = a[i]
        tail = tails[1+key]

        if tail == zero(Tind) # if list is empty
            elems[i] = acPQ_elem{Tind}(zero(Tind), zero(Tind), key)
            tails[1+key] = i
            heads[1+key] = i
        else

            elems[i] = acPQ_elem{Tind}(tail, zero(Tind), key)

            # the following just modifies elems[tail] to point to i
            elems[tail] = acPQ_elem{Tind}(elems[tail].prev, i, elems[tail].key)

            tails[1+key] = i
        end

    end

    return acPQ_rev(elems, heads, tails, minlist, n, n)
end

function pqPop!(pq::acPQ_rev{Tind}) where {Tind}

    if pq.nitems == 0
        error("ApproxPQ is empty")
    end
    while pq.heads[1+pq.minlist] == 0
        pq.minlist = pq.minlist + 1
    end
    i = pq.heads[1+pq.minlist]
    next = pq.elems[i].next

    pq.heads[1+pq.minlist] = next
    if next > 0
        # just changing the first field
        pq.elems[next] = acPQ_elem(zero(Tind), pq.elems[next].next, pq.elems[next].key)
    else # removed the last element
        pq.tails[1+pq.minlist] = 0
    end

    pq.nitems -= 1

    #logwrite(fn,"pq pop: $(i)")

    return i
end

function pqMove!(pq::acPQ_rev{Tind}, i, newkey, oldlist, newlist) where {Tind}

    #logwrite(fn,"pq move: $(i), $(newkey), $(oldlist), $(newlist)")

    prev = pq.elems[i].prev
    next = pq.elems[i].next

    # remove i from its old list
    if next > zero(Tind)
        # really just modifying the first entry: prev
        pq.elems[next] = acPQ_elem{Tind}(prev, pq.elems[next].next, pq.elems[next].key)
    else
        # if this is the tail
        pq.tails[1+oldlist] = prev
    end

    if prev > zero(Tind)
        # really just modifying the second entry: next
        pq.elems[prev] = acPQ_elem{Tind}(pq.elems[prev].prev, next, pq.elems[prev].key)
    else
        # if this is the head
        pq.heads[1+oldlist] = next
    end

    # insert i into its new list
    tail = pq.tails[1+newlist]
    if tail > 0
        pq.elems[tail] = acPQ_elem{Tind}(pq.elems[tail].prev, i, pq.elems[tail].key)
    end
    pq.tails[1+newlist] = i

    pq.elems[i] = acPQ_elem{Tind}(tail, zero(Tind), newkey)

    if pq.heads[1+newlist] == 0
        pq.heads[1+newlist] = i
    end

    #return Void
end

"""
    Decrement the key of element i
    This could crash if i exceeds the maxkey
"""
function pqDec!(pq::acPQ_rev{Tind}, i) where {Tind}

    #logwrite(fn,"pq dec: $(i)")

    oldlist = keyMap(pq.elems[i].key, pq.n)
    newlist = keyMap(pq.elems[i].key - one(Tind), pq.n)

    if newlist != oldlist

        pqMove!(pq, i, pq.elems[i].key - one(Tind), oldlist, newlist)

        if newlist < pq.minlist
            pq.minlist = newlist
        end

    else
        pq.elems[i] = acPQ_elem{Tind}(pq.elems[i].prev, pq.elems[i].next, pq.elems[i].key - one(Tind))
    end

    #return Void
end

"""
    Increment the key of element i
    This could crash if i exceeds the maxkey
"""
function pqInc!(pq::acPQ_rev{Tind}, i) where {Tind}

    #logwrite(fn,"pq inc: $(i)")

    oldlist = keyMap(pq.elems[i].key, pq.n)
    newlist = keyMap(pq.elems[i].key + one(Tind), pq.n)

    if newlist != oldlist

        pqMove!(pq, i, pq.elems[i].key + one(Tind), oldlist, newlist)

    else
        pq.elems[i] = acPQ_elem{Tind}(pq.elems[i].prev, pq.elems[i].next, pq.elems[i].key + one(Tind))
    end

    #return Void
end


text:
abstract Tree{K,V}

type EmptyTree{K,V} <: Tree{K,V}
end

type BinaryTreeNode{K,V} <: Tree{K,V}
    key::  K
    data:: V
    left:: Tree{K,V}
    right::Tree{K,V}
end

type BinaryTree{K,V}
    root:: Tree{K,V}

    BinaryTree() = new(EmptyTree{K,V}())
end

haskey(t::EmptyTree, key) = false
haskey(t::BinaryTree, key) = haskey(t.root, key)

function haskey(t::BinaryTreeNode, key)
    if t.key == key
        true
    elseif key < t.key
        haskey(t.left, key)
    else
        haskey(t.right, key)
    end
end

getindex(t::EmptyTree, k) = throw(KeyError(k))
getindex(t::BinaryTree, k) = t.root[k]

function getindex(t::BinaryTreeNode, key)
    if t.key == key
        t.data
    elseif key < t.key
        t.left[key]
    else
        t.right[key]
    end
end

setindex!{K,V}(t::EmptyTree{K,V}, v, k) = BinaryTreeNode{K,V}(k, v, t, t)
setindex!(t::BinaryTree, v, k) = (t.root = setindex!(t.root, v, k); t)

function setindex!(t::BinaryTreeNode, v, k)
    if t.key == k
        t.data = v
    elseif k < t.key
        t.left = setindex!(t.left, v, k)
    else
        t.right = setindex!(t.right, v, k)
    end
    t
end

delete!(t::EmptyTree, k) = throw(KeyError(k))
delete!(t::BinaryTree, k) = (t.root = delete!(t.root, k); t)

function delete!(t::BinaryTreeNode, k)
    if t.key == k
        if isa(t.right,EmptyTree)
            t = t.left
        elseif isa(t.left,EmptyTree)
            t = t.right
        else
            r = t.right
            t = t.left
            treeinsert!(t, r)
        end
    elseif k < t.key
        t.left = delete!(t.left, k)
    else
        t.right = delete!(t.right, k)
    end
    t
end

treeinsert!(t::EmptyTree, r::BinaryTreeNode) = r

function treeinsert!(t::BinaryTreeNode, r::BinaryTreeNode)
    if r.key < t.key
        t.left = treeinsert!(t.left, r)
    else
        t.right = treeinsert!(t.right, r)
    end
    t
end


text:
abstract type Generator <: Injection end
const Generators = Array{<: Generator, 1}

include("generation/tech_common.jl")
include("generation/econ_common.jl")
include("generation/renewable_generation.jl")
include("generation/thermal_generation.jl")
include("generation/hydro_generation.jl")

struct GenClasses <: Component
    thermal::Union{Nothing,Array{ <: ThermalGen,1}}
    renewable::Union{Nothing,Array{ <: RenewableGen,1}}
    hydro::Union{Nothing,Array{ <: HydroGen,1}}
end
# create iterator for GenClasses
function Base.iterate(iter::GenClasses, state=1)
    if state > length(iter)
        return nothing
    else
        return (getfield(iter, state), state+1)
    end
end
Base.length(iter::GenClasses) = length(fieldnames(GenClasses))
Base.eltype(iter::GenClasses) = Union{Nothing,
                                      Array{ <: ThermalGen,1},
                                      Array{ <: RenewableGen,1},
                                      Array{ <: HydroGen,1}}

# Generator Classifier
function genclassifier(gen::Array{T}) where T <: Generator

    t = [d for d in gen if isa(d, ThermalGen)]
    r = [d for d in gen if isa(d, RenewableGen)]
    h = [d for d in gen if isa(d, HydroGen)]

    #Check for data consistency
    if isempty(t)
        t = nothing
    end

    if isempty(r)
        r = nothing
    end

    if isempty(h)
        h = nothing
    end

    generators = GenClasses(t, r, h)

    return generators
end


text:
function num_var(model)
    nq = model.nq
    nc = model.nc
    nq + nc + nc + nc + nc + nc + nc
end

function num_data(model; nf=length(friction_coefficients(model)))
    nq = model.nq
    nu = model.nu
    nw = model.nw
    2nq + nu + nw + nf + 1
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
using BenchmarkTools

f(x) = x^2
g = f
const g_const = f

F(x) = f(x)
G(x) = g(x)
G_const(x) = g_const(x)

@btime f(2.1)
@btime F(2.1)
@btime G(2.1)
@btime G_const(2.1)
println()
@btime g(2.1)
@btime $g(2.1)
println()
@btime F($(Ref(2.1))[])
@btime G($(Ref(2.1))[])
@btime G_const($(Ref(2.1))[]);

# %%
@code_llvm debuginfo=:none f(2.1)

# %%
@code_llvm debuginfo=:none F(2.1)

# %%
@code_llvm debuginfo=:none G(2.1)

# %%
@code_llvm debuginfo=:none G_const(2.1)

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
export tikz_file, tikz_print

# Code developed by Tara Abrishami

function tikz_string(G::SimpleGraph, label::Bool = false)
    #Outputs a string of tikz code to draw graph G.
    #Label = true if nodes should be labeled in the drawing, false otherwise.
    s = "\\begin{tikzpicture}\n"

    ensure_embed(G)

    xy = getxy(G)
    for v in vlist(G)
        x, y = xy[v]
        s *=
            "\\node[draw,circle,minimum size=10,inner sep=0] (" *
            string(v) *
            ") at (" *
            string(x) *
            "," *
            string(y) *
            ")"
        if label
            s *= " {" * string(v) * "}"
        else
            s *= " {}"
        end
        s *= ";\n"
    end
    for ee in elist(G)
        u, v = ee
        s *= "\\draw (" * string(u) * ") -- (" * string(v) * ");\n"
    end
    s *= "\\end{tikzpicture}\n"
    return s
end

"""
`tikz_print(G)` prints tikz code to draw `G`.
`tikz_print(G,true)` does likewise, with vertex labels drawn.
"""
function tikz_print(G::SimpleGraph, label::Bool = false)
    print(tikz_string(G, label))
end

"""
`tikz_file(G,label,filename)` writes the tikz code for drawing the graph `G`
into `filename`. If `label` is omitted (or `false`) vertex labels are not drawn.
If `filename` is omitted, it defaults to `graph.tex`.
"""
function tikz_file(G::SimpleGraph, label::Bool, filename::String = "graph.tex")
    FILE = open(filename, "w")
    print(FILE, tikz_string(G, label))
    close(FILE)
end

tikz_file(G::SimpleGraph, filename::String = "graph.tex") = tikz_file(G, false, filename)


text:
using Test
using Plots

using GCMMesh
using GCMMesh.Mesh

FT = Float64

limits = (FT(0), FT(1))

@testset "simple line mesh opposing face" begin

    @testset "assert correct element numbering" begin
        mesh = equispaced_line_mesh(limits..., 1, false)
        @test_throws AssertionError Mesh.opposing_face(mesh, 0, 1)
        @test_throws AssertionError Mesh.opposing_face(mesh, 2, 1)
    end

    @testset "1 element line mesh with all periodic boundries" begin
        mesh = equispaced_line_mesh(limits..., 1, true)
        @test Mesh.opposing_face(mesh, 1, 1) == (1, 2)
        @test Mesh.opposing_face(mesh, 1, 2) == (1, 1)
    end

    @testset "2 element line mesh with non-periodic boundaries" begin
        mesh = equispaced_line_mesh(limits..., 2, false)
        @test Mesh.opposing_face(mesh, 1, 1) == (0, 1)
        @test Mesh.opposing_face(mesh, 1, 2) == (2, 1)
        @test Mesh.opposing_face(mesh, 2, 1) == (1, 2)
        @test Mesh.opposing_face(mesh, 2, 2) == (0, 2)
    end

    @testset "2 element line mesh with periodic boundaries" begin
        mesh = equispaced_line_mesh(limits..., 2, true)
        @test Mesh.opposing_face(mesh, 1, 1) == (2, 2)
        @test Mesh.opposing_face(mesh, 1, 2) == (2, 1)
        @test Mesh.opposing_face(mesh, 2, 1) == (1, 2)
        @test Mesh.opposing_face(mesh, 2, 2) == (1, 1)
    end
end

@testset "simple line grid interior faces iterator" begin
    @testset "1×1 element line mesh with periodic boundary" begin
        mesh = equispaced_line_mesh(limits..., 1, true)
        @test length(Mesh.interior_faces(mesh)) == 1
        faces = Mesh.interior_faces(mesh)
        @test tuple(mesh.vertex_neighbors[faces[1], :]...) == (1, 2, 1, 1)
    end
    @testset "1×1 element quad mesh with non-periodic boundaries" begin
        mesh = equispaced_line_mesh(limits..., 1, false)
        @test Mesh.interior_faces(mesh) == nothing
    end
    @testset "2 element line mesh with non-periodic boundaries" begin
        mesh = equispaced_line_mesh(limits..., 2, false)
        @test length(Mesh.interior_faces(mesh)) == 1
        faces = Mesh.interior_faces(mesh)
        @test tuple(mesh.vertex_neighbors[faces[1], :]...) == (1, 2, 2, 1)
    end
end
@testset "simple line grid boundary faces iterator" begin
    @testset "1×1 element line mesh with all periodic boundries" begin
        mesh = equispaced_line_mesh(limits..., 1, true)
        @test Mesh.boundary_faces(mesh, 1) == nothing
        @test Mesh.boundary_faces(mesh, 2) == nothing
    end
    @testset "1×1 element quad mesh with non-periodic boundaries" begin
        mesh = equispaced_line_mesh(limits..., 1, false)
        faces1 = Mesh.boundary_faces(mesh, 1)
        faces2 = Mesh.boundary_faces(mesh, 2)

        @test length(faces1) == 1
        @test length(faces2) == 1
        @test [Mesh.boundary_face_info(mesh, face) for face in faces1] == [(1, 1)]
        @test [Mesh.boundary_face_info(mesh, face) for face in faces2] == [(1, 2)]
        @test Mesh.boundary_tag(mesh, :left) == 1
        @test Mesh.boundary_tag(mesh, :right) == 2
        @test Mesh.boundary_tag(mesh, :bottom) == nothing # no such boundary for line mesh
    end
    @testset "2 element line mesh with non-periodic boundaries" begin
        mesh = equispaced_line_mesh(limits..., 2, false)
        faces1 = Mesh.boundary_faces(mesh, 1)
        faces2 = Mesh.boundary_faces(mesh, 2)
        @test length(faces1) == 1
        @test length(faces2) == 1

        @test [Mesh.boundary_face_info(mesh, face) for face in faces1] == [(1, 1)]
        @test [Mesh.boundary_face_info(mesh, face) for face in faces2] == [(2, 2)]
    end
end

@testset "simple line grid vertex iterator" begin
    @testset "1×1 element line mesh with all periodic boundries" begin
        mesh = equispaced_line_mesh(limits..., 1, true)
        uverts = mesh.unique_verts
        @test length(uverts) == 1
        @test sort(Mesh.unique_vertex_connectivity(mesh, 1)) == sort([(1, 1), (1, 2)])
        @test length(Mesh.unique_vertex_connectivity(mesh, 1)) == 2
    end
    @testset "1×1 element line mesh with non-periodic boundaries" begin
        mesh = equispaced_line_mesh(limits..., 1, false)
        uverts = mesh.unique_verts
        @test length(uverts) == 2
        vc1 = Mesh.unique_vertex_connectivity(mesh, 1)
        vc2 = Mesh.unique_vertex_connectivity(mesh, 2)
        @test length(vc1) == 1
        @test collect(vc1) == [(1, 1)]
        @test length(vc2) == 1
        @test collect(vc2) == [(1, 2)]
    end
    @testset "2 element line mesh with non-periodic boundaries" begin
        mesh = equispaced_line_mesh(limits..., 2, false)
        uverts = mesh.unique_verts
        @test length(uverts) == 3
        vc1 = Mesh.unique_vertex_connectivity(mesh, 1)
        vc2 = Mesh.unique_vertex_connectivity(mesh, 2)
        vc3 = Mesh.unique_vertex_connectivity(mesh, 3)
        @test length(vc1) == 1
        @test vc1 == [(1, 1)]
        @test length(vc2) == 2
        @test vc2 == [(1, 2), (2, 1)]
        @test length(vc3) == 1
        @test vc3 == [(2, 2)]
    end
    @testset "2 element line mesh with periodic boundary" begin
        mesh = equispaced_line_mesh(limits..., 2, true)
        uverts = mesh.unique_verts
        @test length(uverts) == 2
        vc1 = Mesh.unique_vertex_connectivity(mesh, 1)
        vc2 = Mesh.unique_vertex_connectivity(mesh, 2)
        @test length(vc1) == 2
        @test sort(vc1) == [(1, 1), (2, 2)]
        @test length(vc2) == 2
        @test sort(vc2) == [(1, 2), (2, 1)]
    end
end

@testset "simple line grid coordinates" begin
    @testset "1 element line mesh with all periodic boundries" begin
        limits = (FT(0), FT(1))
        mesh = equispaced_line_mesh(limits..., 1, true)
        c1, c2 = Mesh.vertex_coordinates(mesh, 1)
        @test c1 == (0.0,)
        @test c2 == (1.0,)

        limits = (FT(-1), FT(1))
        mesh = equispaced_line_mesh(limits..., 1, false)
        c1, c2 = Mesh.vertex_coordinates(mesh, 1)
        @test c1 == (-1.0,)
        @test c2 == (1.0,)
    end

    @testset "2 element line mesh with non-periodic boundaries" begin
        limits = (FT(0), FT(1))
        mesh = equispaced_line_mesh(limits..., 2, false)
        c1, c2 = Mesh.vertex_coordinates(mesh, 1)
        @test c1 == (0.0,)
        @test c2 == (0.5,)

        c1, c2 = Mesh.vertex_coordinates(mesh, 2)
        @test c1 == (0.5,)
        @test c2 == (1.0,)
    end

    @testset "check coordinate type accuracy" begin
        limits = (big(0.0), big(1.0))
        mesh = equispaced_line_mesh(limits..., 3, false)
        c1, c2 = Mesh.vertex_coordinates(mesh, 1)
        @test eltype(c2) == BigFloat
        @test c2[1] ≈ big(1.0) / big(3.0) rtol = eps(BigFloat)
    end
end

max_stars_repo_name: CliMA/GCMMesh

max_issues_repo_name: CliMA/GCMMesh

max_forks_repo_name: CliMA/GCMMesh


text:
#= -------------------------------------------------------------------------
# @file problema2_simulacion_calentamiento.jl
#
# @date 12/07/17 15:04:34
# @author Martin Noblia
# @email mnoblia@disroot.org
#
# @brief
# simulacion de la ecuacion diferencial que modela el problema de
  calentamiento
# @detail
#
#  Licence:
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
---------------------------------------------------------------------------=#
using DifferentialEquations
using PyPlot
#=------------------------------------------------------------------------------
                              constants
------------------------------------------------------------------------------=#
const  A      =  2.0
const  R      =  2.390057361376673e-5
const  cₑ     =  4.184e6
const  V      =  0.5
const  θ_ref  =  333.0
const  θ₀     =  288.0
const  θ_a    =  293.0
G = (θ_ref - θ_a) / ( R * A)
tspan = (0.0, 1000.0)
θ_dot(t, θ) = ((G * A) / (cₑ * V)) + ((θ_a / (R * cₑ * V))) - ((θ / (R * cₑ * V)))
prob = ODEProblem(θ_dot, θ₀, tspan)
sol = solve(prob, Tsit5(), reltol=1e-12,abstol=1e-12)
#=------------------------------------------------------------------------------
                           plotting
------------------------------------------------------------------------------=#
plot(sol.t, sol[1, :] - 273.0, label=L"\theta_{inv}(t)", linewidth=3)
xlabel(L"t\,[seg]")
ylabel(L"T\,[^{\circ}C]")
grid("on")
matplotlib["rcParams"][:update](["font.size" => 14, "font.family" => "serif"])
legend()


text:
using MultivariateStats, DataFrames, CSV, HTTP, Random, LinearAlgebra, Statistics, Optim, DataFramesMeta, GLM
include("lgwt.jl")

function ps8()
	#Question 1
	url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS8-factor/nlsy.csv"
	df = CSV.read(HTTP.get(url).body)

	ols = lm(@formula(logwage ~ black + hispanic + female + schoolt + gradHS + grad4yr), df)
	print(ols)

	#print(coef(ols))

	#Question 2
	asvabMat=convert(Matrix,df[:,r"asvab"])
	correlation=cor(asvabMat)
	#for all variables:
	#cor(convert(Matrix,df))

	println(correlation)

	#Question 3
	ols2 = lm(@formula(logwage ~ black + hispanic + female + schoolt + gradHS + grad4yr+asvabAR+asvabCS+asvabMK+asvabNO+asvabPC+asvabWK), df)
	print(ols2)

	#I don't know what are these asvab variables but I assume they are some sort of proxy for log wage. They are correlated with each other.
	# if we keep them we have measurement error and if we don't we will have ommited varibale, so any way this models would be biased!

	#Question 4
	#just using transpose instead of reshape
	M = fit(PCA, asvabMat'; maxoutdim=1)

	asvabPCA = MultivariateStats.transform(M, asvabMat')

	print(asvabPCA)

	#Question 5
	M = fit(FactorAnalysis, asvabMat'; maxoutdim=1)
	asvabFA = MultivariateStats.transform(M, asvabMat')

	print(asvabFA)

	#Question 6
	#indexe of M in the formula is correct?
	#getting the random numbers:
	e=rand(Normal(0,1),size(df,1))
	Y=df[:logwage]
	X=df[:,[:black,:hispanic,:female,:schoolt,:gradHS,:grad4yr]]
	asvab=df[:,[:asvabAR, :asvabCS, :asvabMK, :asvabNO, :asvabPC, :asvabWK]]

	function loglike(sigma,m,e,y,x,asv)
		N=size(Y,1)
		J=size(asv,2)

		res_abv=zeros(N,J)
		for i=1:J
			res_abv[:,i]=residuals(lm(@formula(x1 ~ black + hispanic + female +x1_1),hcat(x,asvab[i],e,makeunique=true)))
		end

		wage=residuals(lm(@formula(x1 ~ black + hispanic + female + schoolt + gradHS + grad4yr+x1_1), hcat(x,y,e,makeunique=true)))

		#how many sigma?

		return
	end
end
ps8()


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

using Markdown
using InteractiveUtils

# ╔═╡ f96a2802-2b50-11ec-270d-355639ff9f87
begin
	using Pkg
	cd(joinpath(dirname(@__FILE__),".."))
    Pkg.activate(pwd())
end

# ╔═╡ 602d509e-1775-464f-a438-3f9abc3f0060
# erreur could not load library libGR.dll

# ╔═╡ bd156478-54b6-4085-b50e-0b72da4cab28
begin
	struct Cell
		etat::Int
		age::Int
	end

	const interpretation = Dict(0=>"bord", 1=>"agriculture",2=>"forêt", 3=>"nature", 4=>"urbain")

	Base.show(io::IO, c::Cell) = print(io, "$(interpretation[c.etat]) ($(c.age))")

	Base.zero(::Type{Cell}) = Cell(0,0)
	Base.one(::Type{Cell}) = Cell(1,1)

	"""
		endroitvoisins(i,j)

	pour le tuple d'indices (i,j) rend les indices des voisins autour
	"""
	endroitvoisins(i::Int,j::Int) = CartesianIndex.(Iterators.filter(x-> x[1]≠i || x[2] ≠ j,Iterators.product(i-1:i+1, j-1:j+1)))

	tanh_u(c::Cell) = tanh((c.age + 1/2) / 40)

	const c_status = Dict(1 => 1/4, 3 => 1/2, 2 => 1/8)

	urbanise(c::Cell, voisins::Vector{Cell}) = c_status[c.etat]/8*sum(map(tanh_u, filter(x->isequal(x.etat, 4), voisins)))

	function urbanisation(c::Cell, voisins::Vector{Cell})
		# probabilité de changer
		p_trans = urbanise(c, voisins)
		# aspect random
		if rand() < p_trans
			# on change
			return Cell(4, 0)
		else
			# on reste
			return Cell(c.etat, c.age + 1)
		end
	end

	tanh_u(c::Cell) = tanh(c.age/10)

	function ruralisation(c::Cell, voisins::Vector{Cell})
		# nature => forêt
		if isequal(c.etat,3)
			p_trans = 1/9 * (1 + length(filter(x->isequal(x.etat, 3), voisins))) * tanh( c.age/20)
			if rand() <= p_trans
				return Cell(2, 0)
			end
		# agriculture => nature
		elseif isequal(c.etat,1)
			p_trans = 1/8 * sum(map(tanh_u, filter(x->isequal(x.etat, 1), voisins)))
		else
			return Cell(c.etat, c.age + 1)
		end
	end

	# créer un monde
	function monde(m::Int,n::Int)
		M = zeros(Cell, m,n)
		M[2:end-1, 2:end-1] = reshape([Cell(rand(1:4),0) for _ in 1:(m-2)*(n-2)], m-2,:)
		return M
	end
end

# ╔═╡ 682235f8-d1d2-448a-be19-28dc6d5b13b6
let
	m = 3
	n = 6
	temp = [Cell(rand(1:4),0) for _ in 1:(m-2)*(n-2)]
	reshape(temp, m-2,:)
end

# ╔═╡ 3cb7af07-3805-40b1-a3f2-210c7e51b5e2
monde(4,6)

# ╔═╡ 2c888568-e2d1-48fb-8dfc-7c712b59a4e2
	# urbanisation
	testvec = [Cell(rand(0:4), 10) for _ in 1:8]

# ╔═╡ 6070929b-e9c7-45c7-84d8-0420e84dc280
urbanise(Cell(2,100), testvec)

# ╔═╡ 11cc6324-c6c3-4035-af7d-e14498606a24
let
	A = [1;2;3;4]
	@info filter(x-> x<=2, A), A
	@info filter!(x-> x<=2, A), A
end

# ╔═╡ e7f68f47-feee-49f8-95c9-003442038974
let
	i = 2
	j = 5

end

# ╔═╡ f294e00e-bb93-4751-bc07-ad5335c7340d
A = ones(Cell, 10,10)

# ╔═╡ 0189941a-72e1-4461-af22-c9b3eaef70af
map(x->x.age, A)

# ╔═╡ Cell order:
# ╠═602d509e-1775-464f-a438-3f9abc3f0060
# ╠═f96a2802-2b50-11ec-270d-355639ff9f87
# ╠═bd156478-54b6-4085-b50e-0b72da4cab28
# ╠═682235f8-d1d2-448a-be19-28dc6d5b13b6
# ╠═3cb7af07-3805-40b1-a3f2-210c7e51b5e2
# ╠═2c888568-e2d1-48fb-8dfc-7c712b59a4e2
# ╠═6070929b-e9c7-45c7-84d8-0420e84dc280
# ╠═11cc6324-c6c3-4035-af7d-e14498606a24
# ╠═e7f68f47-feee-49f8-95c9-003442038974
# ╠═f294e00e-bb93-4751-bc07-ad5335c7340d
# ╠═0189941a-72e1-4461-af22-c9b3eaef70af


text:
## Creation

pseudo_ideal(::Type{R},m::RingInt) where R = R(m)
pseudo_ideal(::Type{R},mm::AbstractVector) where R = Ideal(Vector{R}(mm))
pseudo_ideal(::Type{R},m::Ideal{R}) where R = m

function Ideal(m::R) where R<:Ring
    if iszero(m)
        Ideal{R}(R[])
    elseif isunit(m)
        Ideal{R}([one(R)])
    else
        Ideal{R}([m])
    end
end
Ideal(m::T) where T<:Integer = Ideal(ZZ{T}(m))
Ideal(mm::RingInt...) = Ideal(collect(Base.promote_typeof(mm...), mm))
Ideal(mm::AbstractVector{<:RingInt}) = Ideal(gcd(mm))
Ideal(mm::AbstractVector{R}) where R<:MultivariatePolynomial = Ideal{R}(groebnerbase(mm))

# Arithmetic

import Base: in, issubset, ==, intersect, sum, +, *, ^, iszero, zero, one, rem

function in(a::R, id::Ideal{R}) where R<:MultivariatePolynomial
    a, r, d = pdivrem(a, id.base)
    isone(d) && iszero(r)
end

function issubset(id1::Ideal{R}, id2::Ideal{R}) where R<:Polynomial
    for b1 in id1.base
        in(b1, id2) || return false
    end
    true
end

==(id1::Ideal{R}, id2::Ideal{S}) where {R<:Ring,S<:Ring} = id1.base == id2.base
==(id::Ideal{R}, ::Type{R}) where R<:Ring = isone(id)
==(::Type{R}, id::Ideal{R}) where R<:Ring = isone(id)

function +(id1::Ideal{R}, id2::Ideal{R}) where R<:MultivariatePolynomial
    Ideal(vcat(id1.base, id2.base))
end
+(id::Ideal{R}, a::R) where R = Ideal([id.base; a])
+(a::R, id::Ideal{R}) where R = Ideal([a; id.base])

function *(id1::Ideal{R}, id2::Ideal{R}) where R<:MultivariatePolynomial
    C = [ b1 * b2 for b1 in id1.base for b2 in id2.base]
    Ideal(C)
end

*(id::Ideal{R}, a::R) where R = Ideal(id.base .* a)
*(a::R, id::Ideal{R}) where R = Ideal(a .* id.base)

function ^(id::Ideal{R}, p::Integer) where R<:MultivariatePolynomial
    if p == 1
        id
    elseif p == 0
        Ideal([one(R)])
    elseif p == 2
        id * id
    elseif iseven(p)
        (id^(p÷2))^2
    else
        (id^((p-1)÷2))^2 * id
    end
end

# remainder when dividing by an ideal
rem(a::R, id::Ideal{R}) where R<:Ring = divrem(a, id)[2]

function intersect(id1::Ideal{R}, id2::Ideal{R}) where R<:MultivariatePolynomial
    Ideal(_intersect(id1.base, id2.base))
end

function _intersect(v1::V, v2::V) where {R<:MultivariatePolynomial,V<:AbstractVector{R}}
    t = Symbol("intersect_extension")
    Rt = lextend(R, t)
    t, = generators(Rt)
    id3 = Ideal([Rt.(v1) .* t; Rt.(v2) .* (1 - t)])
    R.([x for x in id3.base if multideg(x)[1] == 0])
end

iszero(id::Ideal) = length(id.base) == 0
isone(id::Ideal) = length(id.base) == 1 && isone(id.base[1])
zero(::Type{I}) where {R,I<:Ideal{R}} = Ideal(R[])
one(::Type{I}) where {R,I<:Ideal{R}} = Ideal(R[1])
one(::R) where R<:Ideal = one(R)
zero(::R) where R<:Ideal = zero(R)


text:
using DashBootstrapComponents

card = dbc_card(
    [
        dbc_cardimg(src = "/static/images/placeholder286x180.png", top = true),
        dbc_cardbody([
            html_h4("Card title", className = "card-title"),
            html_p(
                "Some quick example text to build on the card title and " *
                "make up the bulk of the card's content.",
                className = "card-text",
            ),
            dbc_button("Go somewhere", color = "primary"),
        ]),
    ],
    style = Dict("width" => "18rem"),
);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
module NCCL

using ..CUDA
import Libdl

if Sys.iswindows()
    const libnccl = ""
else
    const libnccl = Libdl.find_library("libnccl")
end
isempty(libnccl) && @warn "NCCL cannot be found."

function init()
    @info "NCCL API xxxx"
    # global const API_VERSION = Int(ccall((:cudnnGetVersion,libcudnn),Cint,()))
    # @info("NCCL API $API_VERSION")
end
init()

macro nccl(f, args...)
    quote
        result = ccall(($f,libnccl), Cint, $(map(esc,args)...))
        if result != 0
            p = ccall((:ncclGetErrorString,libnccl), Ptr{UInt8}, (Cint,), result)
            throw(unsafe_string(p))
        end
    end
end

include("types.jl")
include("comm.jl")
include("collective.jl")

end


text:
testnx(g, nx) = isapprox(g.nx, nx)
testny(g, ny) = isapprox(g.ny, ny)

# Physical grid tests
function testdx(g::AbstractGrid{T}) where {T, Ta}
  dxgrid = @. g.x[2:end] - g.x[1:end-1]
  dxones = g.dx*ones(T, size(dxgrid))
  isapprox(dxgrid, dxones)
end

function testdy(g::AbstractGrid{T, Ta}) where {T, Ta}
  dygrid = @. g.y[2:end] - g.y[1:end-1]
  dyones = g.dy*ones(T, size(dygrid))
  isapprox(dygrid, dyones)
end

testx(g) = isapprox(g.x[end]-g.x[1], g.Lx-g.dx)
testy(g) = isapprox(g.y[end]-g.y[1], g.Ly-g.dy)

# Test proper arrangement of fft wavenumbers. k = [0:nx/2, -nx/2+1:-1].
flippednegatives(k, mid) = -reverse(k[mid+2:end])
function testwavenumberalignment(k, nx)
  mid = Int(nx/2)
  positives = k[2:mid]
  negatives = flippednegatives(k, mid)
  isapprox(reshape(positives, mid-1), reshape(negatives, mid-1))
end

testk(g) = testwavenumberalignment(g.k, g.nx)
testl(g) = testwavenumberalignment(g.l, g.ny)
testkr(g) = isapprox(g.k[1:g.nkr], g.kr)

#testk(g) = isapprox(g.k[2:g.nkr-1], flippednegatives(g.k, g.nkr, 1))
#testk(g) = sum(g.k[2:g.nkr-1] .+ reverse(g.k[g.nkr+1:end], dims=1)) == 0.0
#testl(g) = sum(g.l[:, 2:Int(g.ny/2)] .+ reverse(g.l[:, Int(g.ny/2+2):end], dims=2)) == 0.0

function testgridpoints(g::AbstractGrid{T, Ta}) where {T, Ta}
  X, Y = gridpoints(g)
  dXgrid = @. X[2:end, :] - X[1:end-1, :]
  dYgrid = @. Y[:, 2:end] - Y[:, 1:end-1]
  dXones = g.dx*ones(T, size(dXgrid))
  dYones = g.dy*ones(T, size(dYgrid))
  isapprox(dXgrid, dXones) && isapprox(dYgrid, dYones)
end

function testdealias(g::OneDGrid)
  T = typeof(g.Lx)
  fh = ones(Complex{T}, size(g.kr))
  dealias!(fh, g)
  kmax = round(maximum(g.kr)*2/3)
  isapprox(sum(abs.(fh[ g.kr .>= kmax ])), 0)
end

function testdealias(g::TwoDGrid)
  T = typeof(g.Lx)
  fh = ones(Complex{T}, size(g.Krsq))
  dealias!(fh, g)
  kmax = round(maximum(g.kr)*2/3)
  lmax = floor(maximum(g.l)*2/3)

  temp = 0.0
  for j = 1:g.nl, i=1:g.nkr
    if ( g.kr[i] >= kmax ) & ( g.l[j] >= lmax || g.l[j] < -lmax )
      temp += abs.(fh[i, j]) #temp = sum of |fh| for aliased wavenumbers
    end
  end
  isapprox(temp, 0)
end

function testtypedonedgrid(dev::Device, nx, Lx; T=Float64)
  gr = OneDGrid(nx, Lx, T=T)
  typeof(gr.dx)==T && typeof(gr.x[1])==T && typeof(gr.Lx)==T
end

function testtypedtwodgrid(dev::Device, nx, Lx, ny=Lx, Ly=Lx; T=Float64)
  gr = TwoDGrid(nx, Lx, ny, Ly, T=T)
  typeof(gr.dx)==T && typeof(gr.dy)==T && typeof(gr.x[1])==T && typeof(gr.y[1])==T && typeof(gr.Lx)==T && typeof(gr.Ly)==T
end


text:
i = zero(BigInt)    # or i = big(0)
while true
  println(i += 1)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# NetCDFError, error check and netcdf_c.jl from NetCDF.jl (https://github.com/JuliaGeo/NetCDF.jl)
# Copyright (c) 2012-2013: Fabian Gans, Max-Planck-Institut fuer Biogeochemie, Jena, Germany
# MIT

import Base: showerror

# Exception type for error thrown by the NetCDF library
mutable struct NetCDFError <: Exception
    code::Cint
    msg::String
end

"""
    NetCDFError(code::Cint)

Construct a NetCDFError from the error code.
"""
NetCDFError(code::Cint) = NetCDFError(code, nc_strerror(code))

function Base.showerror(io::IO, err::NetCDFError)
    print(io, "NetCDF error: ")
    printstyled(io, err.msg, color=:red)
    print(io, " (NetCDF error code: $(err.code))")
end

"Check the NetCDF error code, raising an error if nonzero"
function check(code::Cint)
    # zero means success, return
    if code == Cint(0)
        return
        # otherwise throw an error message
    else
        throw(NetCDFError(code))
    end
end