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

Evaluate the Laguerre polynomial of degree `p` at `x` using the three term recursion.
"""
function laguerre(x, p::Integer)
    T = typeof( 1-x )
    p₀ = one(T)
    p₁ = 1-x

    if p <= 0
        return p₀
    elseif p == 1
        return p₁
    end

    for n in 2:p
        p₀, p₁ = p₁, ( (2n-1-x)*p₁ - (n-1)*p₀ ) / n
    end

    p₁
end

"""
    laguerre(x, p::Integer, α)

Evaluate the generalised Laguerre polynomial with parameter `α` of degree `p`
at `x` using the three term recursion.
"""
function laguerre(x, p::Integer, α)
    T = typeof( α+1-x )
    p₀ = one(T)
    p₁ = α+1-x

    if p <= 0
        return p₀
    elseif p == 1
        return p₁
    end

    for n in 2:p
        p₀, p₁ = p₁, ( (α+2n-1-x)*p₁ - (α+n-1)*p₀ ) / n
    end

    p₁
end

ext: jl


text:
# DexPtr type holds a gpu (dev>=0) or cpu (dev=-1) allocated pointer.
# We try to minimize the number of actual allocations, which are slow,
# by reusing preallocated but garbage collected pointers.

type DexPtr
    ptr::Cptr
    len::Int
    dev::Int
    parent
end

# We use the DexPtrs type to keep track of allocated and garbage
# collected pointers: We keep one DexPtrs struct per size per device.

type DexPtrs
    used::Int                   # number of allocated pointers
    free::Array{Cptr,1}         # pointers available for reuse
    DexPtrs()=new(0,Array(Cptr,0))
end

# DexFree[dev+2] will hold a dictionary from sizes to DexPtrs for
# device dev.  DexFree[1] reserved for the cpu.  It is initialized in
# DexPtr(nbytes).

DexFree = nothing
initDexFree()=(global DexFree=[ Dict{Int,DexPtrs}() for i=1:gpuCount()+1 ])

# When Julia gc reclaims a DexPtr object, the following special
# finalizer does not actually release the memory, but inserts it in
# the DexFree[dev+2] dict keyed by length in bytes so it can be
# reused.

function freeDexPtr(p::DexPtr)
    ptrs = DexFree[p.dev+2][p.len]
    push!(ptrs.free,p.ptr)
end

# This is the low level DexPtr constructor, it adds the finalizer and
# sets parent to `nothing` which is only needed for shared pointers.

function DexPtr(ptr::Cptr,len::Int,dev::Int)
    kp = DexPtr(ptr,len,dev,nothing)
    finalizer(kp, freeDexPtr)
    return kp
end

# This constructor is used to create a shared pointer.  We need to
# keep the parent field to prevent premature gc of the parent.  The
# child does not need a special finalizer.

function DexPtr(parent::DexPtr, offs::Integer, len::Integer)
    if len < 0 || offs < 1 || offs+len-1 > parent.len
        throw(BoundsError())
    end
    DexPtr(parent.ptr+offs-1, len, parent.dev, parent)
end

# This the main DexPtr constructor.  It tries to avoid actual
# allocation which is slow.  First it tries to find a previously
# allocated and garbage collected pointer in DexFree[dev+2].  If not
# available it tries to allocate a new one (about 10 μs).  Otherwise
# it tries running gc() and see if we get a pointer back (about 75
# ms).  Finally if all else fails, it calls dexgc which cleans up all
# allocated and garbage collected DexPtrs on the current device and
# tries allocation one last time.

function DexPtr(nbytes::Integer)
    DexFree==nothing && initDexFree()
    dev = gpu()
    ptrs = get!(DexPtrs,DexFree[dev+2],nbytes)
    if !isempty(ptrs.free)
        return DexPtr(pop!(ptrs.free),nbytes,dev)
    end
    ptr = dexMalloc(nbytes)
    if ptr != nothing
        ptrs.used += 1
        return DexPtr(ptr,nbytes,dev)
    end
    gc(); print(".")
    if !isempty(ptrs.free)
        return DexPtr(pop!(ptrs.free),nbytes,dev)
    end
    dexgc(); print("+")
    ptr = dexMalloc(nbytes)
    if ptr != nothing
        ptrs.used += 1
        return DexPtr(ptr,nbytes,dev)
    end
    error("Out of gpu memory")
end

# This does the actual allocation, returns `nothing` in case of error
function dexMalloc(nbytes::Int)
    gpu() >= 0 || return(convert(Cptr, pointer(Array(UInt8,nbytes))))
    ptr = Cptr[0]
    ret = ccall(("cudaMalloc","libcudart"),UInt32,(Ptr{Cptr},Csize_t),ptr,nbytes)
    if ret == 0
        return ptr[1]
    else # error("cudaMalloc($nbytes) error $ret")
        return nothing
    end
end

# If you really want to clean up memory you need to call dexgc: Note
# that this only cleans the current device.  It frees all pointers
# previously allocated and garbage collected.
function dexgc(dev=gpu())
    gc_enable(false)
    for v in values(DexFree[dev+2])
        if dev >= 0
            for p in v.free
                @cuda(cudart,cudaFree,(Cptr,),p)
            end
        end
        v.used -= length(v.free)
        empty!(v.free)
    end
    gc_enable(true)
end

# Some utilities
meminfo(i=gpu())=[(k,v.used,length(v.free)) for (k,v) in DexFree[i+2]]
gpufree(i=gpu())=nvmlDeviceGetMemoryInfo(i)[2]

function gpuinfo(msg="",dev=gpu())
    msg != "" && print("$msg ")
    g = nvmlDeviceGetMemoryInfo(dev)
    println((:dev,dev,:total,g[1],:free,g[2],:used,g[3]))
    for (s,u,f) in sort(meminfo(dev), by=(x->x[1]), rev=true)
        println((:size,s,:alloc,u,:avail,f))
    end
end

function memdbg(msg="")
    m = nvmlDeviceGetMemoryInfo()
    c = cudaGetMemInfo()
    x = [0,0]
    for (s,u,f) in meminfo()
        x[1] += s*u
        x[2] += s*f
    end
    println("$msg: ntotal: $(m[1]) nfree: $(m[2]) nused: $(m[3]) ktotal: $(x[1]) kfree: $(x[2])\n nfree+used: $(m[2]+m[3]) ctotal: $(c[2]) cfree: $(c[1]) ktotal-kfree: $(x[1]-x[2]) nused-ktotal: $(m[3]-x[1])")
end


text:
using Pkg
Pkg.activate(".")
using LTMSim
using DataFrames
using SimpleHypergraphs
using Statistics
using Plots
using PyPlot
using Random
using Serialization

h = hg_load("data/got.hgf")

nvalues = [nhv(h)]
runs = 10
data = Dict{String, Vector{Vector{Int}}}()

push!(data, "BinarySearch(H)"=>Vector{Vector{Int}}())
push!(data, "Greedy(H)"=>Vector{Vector{Int}}())
push!(data, "Greedy([H]₂)"=>Vector{Vector{Int}}())
push!(data, "SubTSS(H)"=>Vector{Vector{Int}}())

for n=nvalues

    results1 = Vector{Int}()
    results2 = Vector{Int}()
    results3 = Vector{Int}()
	results4 = Vector{Int}()

    for run=1:runs
        metaV = randMetaV(h)
        metaE = proportionalMetaE(h, 0.5)

        r1 = greedy_tss_2section(h, metaV, metaE)
        r2 = bisect(h, metaV, metaE)
		r3 = greedy_tss(h, metaV, metaE)
        r4 = sub_tss_opt2(h,metaV,metaE)

        push!(results1, r1)
        push!(results2, r2)
        push!(results3, r3)
		push!(results4, r4[1])

    end
    push!(data["Greedy([H]₂)"], results1)
    push!(data["BinarySearch(H)"], results2)
    push!(data["Greedy(H)"], results3)
	push!(data["SubTSS(H)"], results4)
    println("end ", n)

end

#serialize("res/paper/got/got-random.data", data)
data = deserialize("res/paper/got/got-random.data")

labels = Dict{String, String}(
    "BinarySearch(H)" => "StaticGreedy",
    "Greedy(H)" => "DynamicGreedy",
    "Greedy([H]₂)" => L"DynamicGreedy_{[H]_2}",
	"SubTSS(H)" => "SubTSS"
)

ticks = [1]

function set_box_color(bp, color)
    plt.setp(bp["boxes"], color=color)
    plt.setp(bp["whiskers"], color=color)
    plt.setp(bp["caps"], color=color)
    plt.setp(bp["medians"], color=color)
end

clf()

plt.figure(figsize=(7,5))

val = -0.8
c = 1
pos = 1.5

colorz=["#2C7BB6", "#D7191C", "#FF8900", "#33CC33"]

for algo in ["BinarySearch(H)", "Greedy([H]₂)", "Greedy(H)", "SubTSS(H)"]#keys(data)
    global val, c
    b = plt.boxplot(
        data[algo],
        #positions=collect(range(0, stop=length(data[algo])-1)).*2.0.+val,
        positions = [pos+val],
        sym="",
        widths=0.2
    )
    set_box_color(b, colorz[c])
    c+=1
    val+=0.8
end

# draw temporary red and blue lines and use them to create a legend
plt.plot([], c="#2C7BB6", label=labels["BinarySearch(H)"])#label=collect(keys(data))[1])
plt.plot([], c="#D7191C", label=labels["Greedy([H]₂)"])#label=collect(keys(data))[2])
plt.plot([], c="#FF8900", label=labels["Greedy(H)"])#label=collect(keys(data))[3])
plt.plot([], c="#33CC33", label=labels["SubTSS(H)"])#label=collect(keys(data))[3])
plt.legend(fontsize="x-large") #, loc = "lower right"

plt.xticks(range(0.7, 3, step=.8), [1,2,3], fontsize="large", color="#FFFFFF")
plt.yticks(fontsize="x-large")

tick_params(length=0)

#plt.xlim(-2, length(ticks)*2)
plt.ylim(0, 600)

ylabel("Influence set size", fontstyle = "italic", fontsize="xx-large", labelpad=10) #fontweight="semibold",
#xlabel("Thresholds", color="#FFFFFF", fontsize="x-large", fontweight="semibold", labelpad=10)

plt.tight_layout()

gcf()

PyPlot.savefig("res/paper/got/got-random.png")


text:
using BenchmarkTools
using MRFingerprintingRecon
using Test

# TODO
@test 1 == 1


text:
# Run forward model

# Load it
#########
include("svalbard-setup.jl")
plotyes = true
gl,gb,pp,pm,pn,pl,dt = init_svalbard_gl(1081, update_cache=true) #tunabreen
@show gl.gid
@show tunabreen

# Run it
########
pm = BM.MPara(pm, iv_nye=true, iv_dist_exp2=0.01)
print("Time to run forward model:"); gc()
@time fwdsol = BM.fwdm(gb, pp, pm, pn, gb.bdot-gb.dhdt, gb.fsl, gb.temp)

qtot1d, qd1d = BM.postproc1d(gb, pp, pm)

if plotyes
    print("Time to plot:")
    @time display(BM.plot_run(fwdsol, reuse=false, plot2d=[:h,:iv,:ivlog][1]))
end


text:
module TestCategoricalAlgebra

using Test

@testset "FreeDiagrams" begin
  include("FreeDiagrams.jl")
end

@testset "CommutativeDiagrams" begin
  include("CommutativeDiagrams.jl")
end

@testset "Limits" begin
  include("Limits.jl")
end

@testset "Matrices" begin
  include("Matrices.jl")
end

@testset "Sets" begin
  include("Sets.jl")
  include("FinSets.jl")
end

@testset "Relations" begin
  include("FinRelations.jl")
end

@testset "Permutations" begin
  include("Permutations.jl")
end

@testset "CSets" begin
  include("CSetDataStructures.jl")
  include("CSets.jl")
end

@testset "StructuredCospans" begin
  include("StructuredCospans.jl")
end

@testset "CatElements" begin
  include("CatElements.jl")
end

@testset "DataMigration" begin
  include("DataMigration.jl")
end

@testset "DPO" begin
  include("DPO.jl")
end

end


text:
module ParameterTest

using Test
using Parametron

import Parametron: setdirty!, mock_model

@testset "Parameter" begin
    # out-of-place parameter
    val = Ref(true)
    model = mock_model()
    p1 = Parameter{Bool}(() -> val[], model)
    @test p1()

    # test caching
    val[] = false
    @test p1()
    setdirty!(p1)
    @test !p1()

    # in-place parameter
    A = zeros(3, 4)
    Aval = Ref(1.0)
    f = let Aval = Aval # https://github.com/JuliaLang/julia/issues/15276
        m -> m .= Aval[]
    end
    p2 = Parameter(f, A, model)
    @test p2() === A
    @test all(p2() .== 1)

    # caching
    Aval[] = 2.0
    @test p2() === A
    @test all(p2() .== 1.0)
    setdirty!(p2)
    @test all(p2() .== 2.0)
    setdirty!(p2)
    allocs = @allocated p2()
    @test allocs == 0
end

end


text:
module FunctionAccuracyTests

using Base.Math: significand_bits, exponent_bias
using Test, Printf

export FloatIterator, test_acc

struct FloatIterator{T}<:AbstractVector{T}
    min::T
    max::T
    FloatIterator{T}(min, max) where T = min > max ? error("max less than min") : new{T}(min,max)
end

FloatIterator(min::T,max::T) where T = FloatIterator{T}(min,max)
FloatIterator{T}() where T = FloatIterator{T}(T(-Inf),T(Inf))
Base.iterate(it::FloatIterator) = (it.min, nextfloat(it.min))
function Base.iterate(it::FloatIterator{T}, el) where T
    return el == it.max ? nothing : (el, nextfloat(el))
end

function Base.size(it::FloatIterator{T}) where T
    if isless(it.max, 0.0)
        size = reinterpret(Base.uinttype(T), it.max) - reinterpret(Base.uinttype(T), it.min)
    elseif isless(-0.0, it.min)
        size = reinterpret(Base.uinttype(T), it.max) - reinterpret(Base.uinttype(T), it.min)
    else
        size  = reinterpret(Base.uinttype(T), it.max)
        size += reinterpret(Base.uinttype(T), it.min) - reinterpret(Base.uinttype(T), T(-0.0))
    end
    return (size+1,)
end

Base.getindex(it::FloatIterator{T}, i::Int) where T = nextfloat(it.min, i-1)
Base.show(io::IO, mime::MIME"text/plain", it::FloatIterator) =  show(io, it)
Base.show(io::IO, it::FloatIterator) = print(io, "FloatIterator(", it.min, ", ", it.max, ")")

# the following compares the ulp between x and y.
# First it promotes them to the larger of the two types x,y
const infh(T) = prevfloat(T(Inf),4)
function countulp(T, x::AbstractFloat, y::AbstractFloat)
    X, Y = promote(x, y)
    x, y = T(X), T(Y) # Cast to smaller type
    (isnan(x) && isnan(y)) && return 0
    (isnan(x) || isnan(y)) && return 10000
    if isinf(x)
        (sign(x) == sign(y) && abs(y) > infh(T)) && return 0 # relaxed infinity handling
        return 10001
    end
    (x ==  Inf && y ==  Inf) && return 0
    (x == -Inf && y == -Inf) && return 0
    if y == 0
        x == 0 && return 0
        return 10002
    end
    if isfinite(x) && isfinite(y)
        return T(abs(X - Y) / ulp(y))
    end
    return 10003
end

const DENORMAL_MIN(::Type{Float64}) = 2.0^-1074
const DENORMAL_MIN(::Type{Float32}) = 2f0^-149
const DENORMAL_MIN(::Type{Float16}) = Float16(6.0e-8)

function ulp(x::T) where {T<:AbstractFloat}
    x = abs(x)
    x == T(0.0) && return DENORMAL_MIN(T)
    val, e = frexp(x)
    return max(ldexp(T(1.0), e - significand_bits(T) - 1), DENORMAL_MIN(T))
end

countulp(x::T, y::T) where {T <: AbstractFloat} = countulp(T, x, y)
strip_module_name(f::Function) = last(split(string(f), '.')) # strip module name from function f

function test_acc(fun_table::Dict, xx; tol=1.5, debug = true, tol_debug = 5)
    T = eltype(xx)
    @testset "accuracy $(strip_module_name(xfun))" for (xfun, fun) in fun_table
        rmax = 0.0
        rmean = 0.0
        xmax = map(zero, first(xx))
        tol_debug_failed = 0
        for x in xx
            q = xfun(x...)
            c = fun(map(big, x)...)
            u = countulp(T, q, c)
            rmax = max(rmax, u)
            xmax = rmax == u ? x : xmax
            rmean += u
            if debug && u > tol_debug
                tol_debug_failed += 1
                #@printf("%s = %.20g\n%s  = %.20g\nx = %.20g\nulp = %g\n", strip_module_name(xfun), q, strip_module_name(fun), T(c), x, u)
            end
        end
        if debug
            println("Tol debug failed $(100tol_debug_failed / length(xx))% of the time.")
        end
        rmean = rmean / length(xx)
        println(strip_module_name(xfun))
        println("max $rmax at x = $xmax")
        println("mean $rmean")
        t = @test rmax <= tol
    end
end

test_acc(f::Function, xx; tol = 1.5, debug = true, tol_debug = 5) = test_acc(Dict(f=>f), xx; tol = tol, debug = debug, tol_debug = tol_debug)
test_acc(f::Function, min, max; tol = 1.5, debug = true, tol_debug = 5) = test_acc(Dict(f=>f), FloatIterator(min,max); tol = tol, debug = debug, tol_debug = tol_debug)
test_acc(f::Function, T::Type; tol = 1.5, debug = true, tol_debug = 5) = test_acc(Dict(f=>f), FloatIterator(typemin(T),typemax(T); tol = tol, debug = debug, tol_debug = tol_debug))
end # module


text:
# Reluval, DLV, Sherlock, BaB
# network with one output node

using NeuralVerification
using Test

macro no_error(ex)
    quote
        try $(esc(ex))
            true
        catch
            false
        end
    end
end

at = @__DIR__

# Wide MNIST network with 1-d output

mnist_wide_1 = read_nnet("$at/../examples/networks/mnist-1-100-1.nnet")
input_center = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,121,254,136,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,13,230,253,248,99,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,118,253,253,225,42,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,61,253,253,253,74,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,32,206,253,253,186,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,211,253,253,239,69,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,254,253,253,133,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,142,255,253,186,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,149,229,254,207,21,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,54,229,253,254,105,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,152,254,254,213,26,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,112,251,253,253,26,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,29,212,253,250,149,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,36,214,253,253,137,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,75,253,253,253,59,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,93,253,253,189,17,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,224,253,253,84,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,43,235,253,126,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,99,248,253,119,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,225,235,49,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
output_center = [2713.019467032775]

in_epsilon = 1 #0-255
out_epsilon = 10 #logit domain

input_low = input_center .- in_epsilon
input_high = input_center .+ in_epsilon

output_low = output_center .- out_epsilon
output_high = output_center .+ out_epsilon

inputSet = Hyperrectangle(low=input_low, high=input_high)
outputSet = Hyperrectangle(low=output_low, high=output_high)

problem_hyperrect_wide_1 = Problem(mnist_wide_1, inputSet, outputSet)

# Reluval

solver = ReluVal(max_iter = 1)
print("\nReluval - Wide-1")
timed_result = @timed solve(solver, problem_hyperrect_wide_1)
t = timed_result[2]
out = timed_result[1]
print(" - Time: $t s - Output: $out")

# DLV
optimizer = GLPKSolverMIP()
solver = DLV(1.0)
print("\nDLV - Wide-1")
timed_result = @timed solve(solver, problem_hyperrect_wide_1)
t = timed_result[2]
out = timed_result[1]
print(" - Time: $t s - Output: $out")

# Sherlock

optimizer = GLPKSolverMIP()
solver = Sherlock(optimizer, 1.0)
print("\nSherlock - Wide-1")
timed_result = @timed solve(solver, problem_hyperrect_wide_1)
t = timed_result[2]
out = timed_result[1]
print(" - Time: $t s - Output: $out")

# BaB
optimizer = GLPKSolverMIP()
solver = BaB(optimizer, 0.1)
#print("\nBaB - Wide-1")
#timed_result = @timed solve(solver, problem_hyperrect_wide_1)
t = timed_result[2]
out = timed_result[1]
#print(" - Time: $t s - Output: $out")


text:
# -*- coding: utf-8 -*-
import Revise
import GLRenderer as GL
import ThreeDP3 as T
import Images as I
import MiniGSG as S
import NearestNeighbors as NN
import LinearAlgebra
import PoseComposition: Pose, IDENTITY_POSE, IDENTITY_ORN
using Gen
import Rotations as R
import DataStructures as DS
import MeshCatViz

import MeshCatViz
MeshCatViz.setup_visualizer()

YCB_DIR = "/home/nishadg/mcs/ThreeDVision.jl/data/ycbv2"
world_scaling_factor = 100.0
id_to_cloud, id_to_shift, id_to_box  = T.load_ycbv_models_adjusted(YCB_DIR, world_scaling_factor);
all_ids = sort(collect(keys(id_to_cloud)));

IDX =200
@show T.get_ycb_scene_frame_id_from_idx(YCB_DIR,IDX)
gt_poses, ids, rgb_image, gt_depth_image, cam_pose, original_camera = T.load_ycbv_scene_adjusted(
    YCB_DIR, IDX, world_scaling_factor, id_to_shift
);
camera = T.scale_down_camera(original_camera, 4)
@show camera
img = I.colorview(I.Gray, gt_depth_image ./ maximum(gt_depth_image))
I.colorview(I.RGB, permutedims(Float64.(rgb_image)./255.0, (3,1,2)))

# +
resolution = 0.5

renderer = GL.setup_renderer(camera, GL.DepthMode())
for id in all_ids
    cloud = id_to_cloud[id]
    v,n,f = GL.mesh_from_voxelized_cloud(GL.voxelize(cloud, resolution), resolution)
    GL.load_object!(renderer, v, f)
end
depth_image = GL.gl_render(
    renderer, ids, gt_poses, T.IDENTITY_POSE)
depth_image[depth_image .== 50000.0] .= 200.0

function get_cloud_from_ids_and_poses(ids, poses, camera_pose)
    depth_image = GL.gl_render(renderer, ids, poses, camera_pose)
    cloud = GL.flatten_point_cloud(GL.depth_image_to_point_cloud(depth_image, camera))
    cloud = cloud[:,cloud[3,:] .< (camera.far - 1.0)]
    if size(cloud)[2] == 0
        cloud = zeros(3,1)
    end
    cloud
end

img = I.colorview(I.Gray, depth_image ./ maximum(depth_image))
# -

function reflection_matrix_4x4(center, direction)
    T =  Matrix{Float64}(LinearAlgebra.I,4,4)
    T[1:3,4] = -center

    dir = direction ./ LinearAlgebra.norm(direction)
    R = [
        1-2*dir[1]^2   -2*dir[2]*dir[1] -2*dir[3]*dir[1] 0;
        -2*dir[2]*dir[1]  1-2*dir[2]^2 -2*dir[2]*dir[3] 0;
        -2*dir[3]*dir[1]  -2*dir[2]*dir[3] 1-2*dir[3]^2 0;
        0 0 0 1
    ]
    inv(T) * R * T
end

# +
IDX = 3
obj_pose = T.get_c_relative_to_a(cam_pose,gt_poses[IDX])
depth_image = GL.gl_render(renderer, [ids[IDX]], [obj_pose], cam_pose)

object_cloud = GL.flatten_point_cloud(GL.depth_image_to_point_cloud(depth_image, camera))
object_cloud = object_cloud[:,object_cloud[3,:] .< (camera.far - 1.0)]
object_cloud = T.move_points_to_frame_b(object_cloud, cam_pose)
mins,maxs = T.min_max(object_cloud)
mins = mins .- 6.5
maxs = maxs .+ 6.5

RESOLUTION = resolution
mins,maxs = floor.(Int,mins ./ RESOLUTION) * RESOLUTION , ceil.(Int, maxs ./RESOLUTION) * RESOLUTION
dimensions = [length(collect(mins[1]:RESOLUTION:maxs[1])),
              length(collect(mins[2]:RESOLUTION:maxs[2])),
              length(collect(mins[3]:RESOLUTION:maxs[3]))]
ALL_VOXELS = hcat([[a,b,c] for a in collect(mins[1]:RESOLUTION:maxs[1])
                           for b in collect(mins[2]:RESOLUTION:maxs[2])
                           for c in collect(mins[3]:RESOLUTION:maxs[3])]...)

occupied, occluded, free = T.get_occ_ocl_free(T.get_points_in_frame_b(ALL_VOXELS,cam_pose), camera, depth_image, resolution)
MeshCatViz.reset_visualizer()
MeshCatViz.viz(ALL_VOXELS[:,occupied]./200.0;channel_name=:occupied,color=I.colorant"red")
MeshCatViz.viz(ALL_VOXELS[:,occluded]./200.0;channel_name=:occluded,color=I.colorant"black")
cloud = ALL_VOXELS[:,occupied]

tree_occupied = NN.KDTree(ALL_VOXELS[:, occupied]);
tree_non_free = NN.KDTree(ALL_VOXELS[:, occupied .| occluded]);
# -

plane_points = hcat([
    [x,y,0.0]
    for x in -10.0:1.0:10.0 for y in -10.0:1.0:10.0
]...)
@show size(plane_points)

function score_reflection_plane(center, direction; verbose=false)
   Reflect = reflection_matrix_4x4(
        center, direction
    )
    new_c = (Reflect * (vcat(cloud, ones(1,size(cloud)[2]))))[1:3,:]
    new_c = T.voxelize(new_c, resolution)
    all_idxs = NN.inrange(tree_non_free, new_c, resolution ./ 2);
    matched = [length(x)>0 for x in all_idxs]
    score = sum(matched) / length(matched)
    if verbose
        println(score)
    end

    all_idxs = NN.inrange(tree_non_free, new_c, resolution ./ 2);
    matched = [length(x)==0 for x in all_idxs]
    if verbose
        println(-sum(matched) / length(matched))
    end

    score -= sum(matched) / length(matched)
    return score, new_c
end

MeshCatViz.reset_visualizer()
MeshCatViz.viz(cloud./500.0;channel_name=:occupied,color=I.colorant"red")
T.min_max(cloud)

# +
param_sweep = [
    ([x,y,10.0],[sin(θ),cos(θ),0.0])
    for x in -6.0:resolution:-2.0 for y in -1.0:resolution:5.0 for θ in 0:0.05:(2*pi)
];
@show size(param_sweep)

# for (c,d) in param_sweep[1:500]
#     MeshCatViz.viz(get_plane_points(c,d) ./ 500.0;channel_name=:plane, color=I.colorant"green")
#     sleep(0.005)
# end
# -

scores = [score_reflection_plane(c, d)[1] for (c,d) in param_sweep];

c,d = param_sweep[argmax(scores)]
score, new_c = score_reflection_plane(c,d;verbose=true)
@show score
MeshCatViz.reset_visualizer()
MeshCatViz.viz(cloud./200.0;channel_name=:occupied,color=I.colorant"red")
MeshCatViz.viz(new_c./200.0;channel_name=:occluded,color=I.colorant"black")

# +

function get_plane_points(center, direction)
    u = [0,0,1]
    v = direction
    u = u ./ LinearAlgebra.norm(u)
    v = v ./ LinearAlgebra.norm(v)
    half = (u+v) ./ LinearAlgebra.norm(u+v)
    q = R.QuatRotation(LinearAlgebra.dot(u,half),LinearAlgebra.cross(u,half)...)
    q * plane_points .+ center
end

MeshCatViz.viz(get_plane_points(c,d) ./ 200.0;channel_name=:plane, color=I.colorant"green")
# -

perm = sortperm(-1.0 .* scores)
c,d = param_sweep[perm[5]]
score, new_c = score_reflection_plane(c,d;verbose=true)
@show score
MeshCatViz.reset_visualizer()
MeshCatViz.viz(cloud./300.0;channel_name=:occupied,color=I.colorant"red")
MeshCatViz.viz(new_c./100.0;channel_name=:occluded,color=I.colorant"black")

# +
plane_cloud = [

]
# -

score_reflection_plane(best_grid_center, d)

all_idxs

obs_cloud = GL.flatten_point_cloud(GL.depth_image_to_point_cloud(gt_depth_image, original_camera))
obs_cloud = GL.move_points_to_frame_b(obs_cloud, cam_pose)
obs_cloud =obs_cloud[:, 1.0 .< obs_cloud[3, :] .< 40.0];
obs_cloud =obs_cloud[:, -25.0 .< obs_cloud[1, :] .< 25.0];
obs_cloud =obs_cloud[:, -25.0 .< obs_cloud[2, :] .< 25.0];
voxelized_obs_cloud = T.voxelize(obs_cloud, resolution)
MeshCatViz.viz(voxelized_obs_cloud ./ 100.0)

get_cloud_func(id,p) = T.move_points_to_frame_b(T.voxelize(get_cloud_from_ids_and_poses(
                                [id], [p], cam_pose),resolution), cam_pose)

blobs = T.get_entities_from_assignment(voxelized_obs_cloud,
    T.dbscan_cluster(voxelized_obs_cloud, radius=0.6))

@show ids
obs_blob = blobs[2]
MeshCatViz.reset_visualizer()
MeshCatViz.viz(obs_blob ./ 100.0)

# +
parent_box = S.Box(40.0,40.0,0.1)
parent_pose = Pose(0.0, 0.0, -1.0)
parent_face = :top

c1_tree = NN.KDTree(obs_blob)
id_to_face_cloud_and_score =  DS.DefaultDict{Any,Vector{Any}}(()->[])

for id in all_ids
    print(id," ")
    for child_face in S.BOX_SURFACE_IDS
        child_box = id_to_box[id]
        (x,y,ang) = T.icp_project_to_planar_contact(
            obs_blob, p -> get_cloud_func(id,p),
            parent_pose, parent_box, parent_face,
            child_box, child_face,
            0.0, 0.0, 0.0;
            c1_tree= c1_tree, outer_iterations=10
        )

        contact = S.ShapeContact(parent_face, Real[], child_face, Real[], S.PlanarContact(x,y,ang))
        p = parent_pose * S.getRelativePoseFromContact(parent_box, child_box, contact)
        c = get_cloud_func(id,p)
        score = Gen.logpdf(
            T.uniform_mixture_from_template,
            obs_blob, c, 0.001, 0.5*2,(-1000.0, 1000.0, -1000.0,1000.0,-1000.0,1000.0)
        )

        push!(id_to_face_cloud_and_score[id], (face=child_face, cloud=c, score=score, pose=p, contact_params=(x,y,ang)))
    end
#     println(size(id_to_face_cloud_and_score[id]))
    sort!(id_to_face_cloud_and_score[id],by=x->x.score,rev=true)
end

# +
scores = [id_to_face_cloud_and_score[id][1].score for id in all_ids]
perm = sortperm(-scores)
best_score_id = perm[1]
@show best_score_id
c = id_to_face_cloud_and_score[best_score_id][2].cloud

MeshCatViz.reset_visualizer()
MeshCatViz.viz(obs_blob ./100.0; color=I.colorant"black", channel_name=:obs_cloud)
MeshCatViz.viz(c ./100.0; color=I.colorant"red", channel_name=:gen_cloud)

# +
c = id_to_face_cloud_and_score[12][1].cloud

MeshCatViz.reset_visualizer()
MeshCatViz.viz(obs_blob ./100.0; color=I.colorant"black", channel_name=:obs_cloud)
MeshCatViz.viz(c ./100.0; color=I.colorant"red", channel_name=:gen_cloud)

# +

id_and_id_to_face_cloud_and_score =  DS.DefaultDict{Any,Vector{Any}}(()->[])

for id_1 in all_ids
    print(id_1," ")

    data = id_to_face_cloud_and_score[id_1][1]

    unexplained_idxs = T.get_unexplained_points(obs_blob, data.cloud)
    unexplained_points = obs_blob[:, unexplained_idxs]

    c1_tree = NN.KDTree(unexplained_points)

    parent_pose = data.pose
    parent_face = T.opposite_face(data.face)
    parent_box = id_to_box[id_1]

    for id in all_ids
        for child_face in S.BOX_SURFACE_IDS
            child_box = id_to_box[id]

            (x,y,ang) = T.icp_project_to_planar_contact(
                unexplained_points, p -> get_cloud_func(id,p), parent_pose, parent_box, parent_face,
                child_box, child_face,
                0.0,0.0,0.0;
                c1_tree= c1_tree
            )

            contact = S.ShapeContact(parent_face, Real[], child_face, Real[], S.PlanarContact(x,y,ang))
            p = parent_pose * S.getRelativePoseFromContact(parent_box, child_box, contact)
            c = get_cloud_func(id,p)
            score = Gen.logpdf(
                T.uniform_mixture_from_template,
                unexplained_points, c, 0.01, 0.5*2,(-1000.0, 1000.0, -1000.0,1000.0,-1000.0,1000.0)
            )

            push!(id_and_id_to_face_cloud_and_score[(id_1,id)],
                (face=child_face, cloud=c, score=score, pose=p, contact_params=(x,y,ang)))
        end
        sort!(id_and_id_to_face_cloud_and_score[(id_1,id)],by=x->x.score,rev=true)
    end
end

# -

id_and_id_to_score = Dict()
for id_1 in all_ids
    for id in all_ids
        pose_1 = id_to_face_cloud_and_score[id_1][1].pose
        pose_2 = id_and_id_to_face_cloud_and_score[(id_1,id)][1].pose

        c = T.move_points_to_frame_b(T.voxelize(get_cloud_from_ids_and_poses(
                                [id_1,id], [pose_1, pose_2], cam_pose), 0.5), cam_pose)

        score = Gen.logpdf(
            T.uniform_mixture_from_template,
            obs_blob, c, 0.01, 0.5*2,(-1000.0, 1000.0, -1000.0,1000.0,-1000.0,1000.0)
        )
        id_and_id_to_score[(id_1,id)] = score
    end
end

# +
(id_1,id) = argmax(id_and_id_to_score)
@show (id_1,id)
pose_1 = id_to_face_cloud_and_score[id_1][1].pose
pose_2 = id_and_id_to_face_cloud_and_score[(id_1,id)][1].pose

c = T.move_points_to_frame_b(T.voxelize(get_cloud_from_ids_and_poses(
                        [id_1,id], [pose_1, pose_2], cam_pose),0.5), cam_pose)
MeshCatViz.viz(obs_blob ./100.0; color=I.colorant"black", channel_name=:obs_cloud)
MeshCatViz.viz(c ./100.0; color=I.colorant"red", channel_name=:gen_cloud)
# -

c

ids

ids

# +
id_1 = 13
    data = id_to_face_cloud_and_score[id_1][1]

    unexplained_idxs = T.get_unexplained_points(obs_blob, data.cloud)
    unexplained_points = obs_blob[:, unexplained_idxs]

    c1_tree = NN.KDTree(unexplained_points)

    parent_pose = data.pose
    parent_face = T.opposite_face(data.face)
    parent_box = id_to_box[id_1]

id = 3
child_face=:left
(x,y,ang) = T.icp_project_to_planar_contact(
    unexplained_points, p -> get_cloud_func(id,p), parent_pose, parent_box, parent_face,
    child_box, child_face,
    0.0,0.0,0.0;
    c1_tree= c1_tree, outer_iterations=5, iterations=10
)

contact = S.ShapeContact(parent_face, Real[], child_face, Real[], S.PlanarContact(x,y,ang))
p = parent_pose * S.getRelativePoseFromContact(parent_box, child_box, contact)
c = get_cloud_func(id,p)

MeshCatViz.reset_visualizer()
MeshCatViz.viz(unexplained_points; color=I.colorant"black", channel_name=:obs_cloud)
MeshCatViz.viz(c; color=I.colorant"red", channel_name=:gen_cloud)

# +

c = id_and_id_to_face_cloud_and_score[(13,3)][1].cloud

MeshCatViz.reset_visualizer()
MeshCatViz.viz(obs_blob; color=I.colorant"black", channel_name=:obs_cloud)
MeshCatViz.viz(c; color=I.colorant"red", channel_name=:gen_cloud)
# -

unexplained_idxs = T.get_unexplained_points(obs_blob, data.cloud)
unexplained_points = obs_blob[:, unexplained_idxs]
MeshCatViz.viz(unexplained_points; color=I.colorant"red", channel_name=:gen_cloud)

ids

pair = (13,:right)
@show id_and_face_to_score[pair]
c = id_and_face_to_cloud[pair]
c = id_and_face_to_cloud[best_pair]

MeshCatViz.reset_visualizer()
MeshCatViz.viz(obs_blob; color=I.colorant"black", channel_name=:obs_cloud)
MeshCatViz.viz(c; color=I.colorant"red", channel_name=:gen_cloud)

ids

clouds = [[
    let
        p = T.icp_object_pose(IDENTITY_POSE, obs_blob,
                              p-> T.move_points_to_frame_b(T.voxelize(get_cloud_from_ids_and_poses(
                                [id], [p], cam_pose),0.5), cam_pose), true)

        c = T.move_points_to_frame_b(T.voxelize(
                get_cloud_from_ids_and_poses([id], [p], cam_pose), 0.5), cam_pose)
    end
    for _ in 1:10
] for id in all_ids];
scores = [
    maximum([
        Gen.logpdf(T.uniform_mixture_from_template,obs_blob, c, 0.01, 0.5*2,(-1000.0, 1000.0, -1000.0,1000.0,-1000.0,1000.0))
        for c in x
    ])
    for x in clouds
]

best_id = argmax(scores)
@show best_id
c = clouds[best_id][1]
MeshCatViz.viz(c; color=I.colorant"red", channel_name=:h1)
MeshCatViz.viz(obs_blob; color=I.colorant"black", channel_name=:h2)

ids

ids

c

MeshCatViz.reset_visualizer()
MeshCatViz.viz(T.get_obs_cloud(trace); color=I.colorant"black", channel_name=:obs_cloud)
MeshCatViz.viz(T.get_gen_cloud(trace); color=I.colorant"red", channel_name=:gen_cloud)

traces = [trace];

for i in 1:length(ids)
    if !T.isFloating(T.get_structure(trace),i)
        continue
    end
    pose_addr = T.floating_pose_addr(i)
    for _ in 1:50
        trace, acc = T.drift_move(trace, pose_addr, 0.5, 10.0)
        if acc traces = push!(traces, trace) end

        trace, acc = T.drift_move(trace, pose_addr, 1.0, 100.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.drift_move(trace, pose_addr, 1.5, 1000.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.drift_move(trace, pose_addr, 0.1, 1000.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.drift_move(trace, pose_addr, 0.1, 100.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.drift_move(trace, pose_addr, 0.5, 5000.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.pose_flip_move(trace, pose_addr, 1, 1000.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.pose_flip_move(trace, pose_addr, 2, 1000.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.pose_flip_move(trace, pose_addr, 3, 1000.0)
        if acc traces = push!(traces, trace) end
    end
end

CRAZY_POSE = Pose(-100.0 * ones(3), T.IDENTITY_ORIENTATION)
for i in 1:length(ids)
    if !T.isFloating(T.get_structure(trace),i)
        continue
    end
    addr = T.floating_pose_addr(i)

    t, = Gen.update(trace, Gen.choicemap(T.floating_pose_addr(i) => CRAZY_POSE))
    valid = T.get_unexplained_obs_cloud(t)

    p = T.icp_object_pose(trace[addr], trace[:obs][:,valid],
                          p-> get_cloud_from_ids_and_poses(
                            [ids[i]], [p], cam_pose), false)
    for _ in 1:3
        trace, acc = T.pose_mixture_move(trace, addr, [trace[addr], p], [0.5, 0.5], 0.001, 5000.0)
        if acc traces = push!(traces, trace) end
    end
end

MeshCatViz.reset_visualizer()
MeshCatViz.viz(T.get_obs_cloud(trace); color=I.colorant"black", channel_name=:obs_cloud)
MeshCatViz.viz(T.get_gen_cloud(trace); color=I.colorant"red", channel_name=:gen_cloud)

depth_image = GL.gl_render(
    renderer, ids, T.get_poses(trace)[1:end-1], cam_pose)
depth_image[depth_image .== 50000.0] .= 200.0
img = I.colorview(I.Gray, depth_image ./ maximum(depth_image))

# +
renderer_texture = GL.setup_renderer(original_camera, GL.TextureMode())
obj_paths = T.load_ycb_model_obj_file_paths(YCB_DIR)
texture_paths = T.load_ycb_model_texture_file_paths(YCB_DIR)
for id in all_ids
    v,n,f,t = renderer_texture.gl_instance.load_obj_parameters(
        obj_paths[id]
    )
    v = v * world_scaling_factor
    v .-= id_to_shift[id]'

    GL.load_object!(renderer_texture, v, n, f, t,
        texture_paths[id]
    )
end

# -

rgb_image, depth_image = GL.gl_render(
    renderer_texture, ids, T.get_poses(trace)[1:end-1], cam_pose)
I.colorview(I.RGBA, permutedims(rgb_image,(3,1,2)))


text:
"""
    is_key(candidate)

Checks if provided candidate key matches the format for AAindex accession
numbers.
"""
function is_key(candidate::String)::Bool
    match(r"\w{4}\d{6}", candidate) !== nothing
end

"""
    aaindex_by_id(id)

Load an AAindex by its id.
"""
function aaindex_by_id(id::String)
    if !is_key(id)
        throw(ArgumentError("$id is not a valid AAindex identifier"))
    end

    try
        jldopen(joinpath(datadep"AAindex", "aaindex.jld2"), "r") do file
            return file["entries"][id]
        end
    catch
        throw(ArgumentError("$id is not a valid AAindex identifier"))
    end
end

function search(term::String)
    database_path = joinpath(datadep"AAindex", "aaindex.jld2")

    ids, descriptions, types = jldopen(database_path, "r") do file
        file["ids"], file["descriptions"], file["types"]
    end

    match_indices = Set()

    for (i, id) in enumerate(ids)
        if id == term
            push!(match_indices, i)
        end
    end

    for (i, description) in enumerate(descriptions)
        if occursin(term |> lowercase, description |> lowercase)
            push!(match_indices, i)
        end
    end

    [(ids[i], descriptions[i], types[i]) for i in match_indices]
end


text:
using BijectiveHilbert
using SafeTestsets
using Test

include("suite.jl")

@testset "BijectiveHilbert.jl" begin
include("test_bitops.jl")
include("test_gray_code.jl")
include("test_simple2d.jl")
include("test_hamilton.jl")
include("test_global_gray.jl")
include("test_facecontinuous.jl")
end


text:
using SparseArrays
using LinearAlgebra
using Plots
using RecipesBase

include("C:/Users/lucas/OneDrive/Documents/GitHub/GDB.jl/experimental/ThermalLBM.jl")

function boundary_thermal_lattice_boltzmann_method(streamtensor, Dx, Dy, Dxx, Dyy, xres, yres, f, g, fb, gb, K1, K2, tau, tauc, deltax, deltat)
    c = deltax/deltat
    nu = tau*c^2/3 # kinetic viscosity

    fnew = streamtensor * f
    gnew = streamtensor * g

    rho = zeros((xres*yres))
    for i=1:xres*yres
        rho[i] = sum(fnew[i,:])
    end

    p = zeros(xres*yres)
    for i=1:xres*yres
        p[i] = c^2/3 * rho[i]
    end

    u = zeros(xres*yres, 2)
    for i=1:xres*yres
        ux = (fnew[i,2] - fnew[i,4] + 1/sqrt(2) * (fnew[i,6] + fnew[i,9] - fnew[i,7] - fnew[i,8])) / rho[i]
        uy = (fnew[i,3] - fnew[i,5] + 1/sqrt(2) * (fnew[i,6] + fnew[i,8] - fnew[i,9] - fnew[i,8])) / rho[i]
        u[i,1] = ux
        u[i,2] = uy
    end

    q = zeros(xres*yres, 9)

    grad = [Dx, Dy]
    divu = Dx*u[:,1] + Dy*u[:,2]
    gradp = zeros(xres*yres, 2)
    gradp[:,1] = Dx*p[:]
    gradp[:,2] = Dy*p[:]
    laplu = zeros(xres*yres, 2)
    laplu[:,1] = (Dxx+Dyy)*u[:,1]
    laplu[:,2] = (Dxx+Dyy)*u[:,2]

    unitvecs = zeros(9,2)
    unitvecs[:,1] = [0, c, 0, -c, 0, 1/sqrt(2)*c, -1/sqrt(2)*c, -1/sqrt(2)*c, 1/sqrt(2)*c]
    unitvecs[:,2] = [0, 0, c, 0, -c, 1/sqrt(2)*c, 1/sqrt(2)*c, -1/sqrt(2)*c, -1/sqrt(2)*c]

    vgradu = zeros(xres*yres, 9, 2)
    for k=1:9
        vgradu[:,k,1] = (unitvecs[k,1] .- u[:,1]) .* (Dx*u[:,1]) + (unitvecs[k,2] .- u[:,2]) .* (Dy*u[:,1])
        vgradu[:,k,2] = (unitvecs[k,1] .- u[:,1]) .* (Dx*u[:,2]) + (unitvecs[k,2] .- u[:,2]) .* (Dy*u[:,2])
    end

    for i=1:xres*yres
        for a=1:9
            q[i,a] = - 1/rho[i]*dot(unitvecs[a,:]-u[i,:], gradp[i,:])
            + nu*(dot(unitvecs[a,:]-u[i,:], laplu[i,:]))
            + dot(unitvecs[a,:]-u[i,:], vgradu[i,a,:])
        end
    end

    w = zeros(9)
    w[1] = 4/9
    w[2] = 1/9
    w[3] = 1/9
    w[4] = 1/9
    w[5] = 1/9
    w[6] = 1/36
    w[7] = 1/36
    w[8] = 1/36
    w[9] = 1/36

    s = zeros((xres*yres, 9))
    usquared = zeros((xres*yres))
    for i=1:xres*yres
        usquared[i] = dot(u[i,:],u[i,:])
        s[i,:] = w[:] .* (3*(unitvecs*u[i,:])/c^2 .+ 9/2*((unitvecs*u[i,:])).^2/c^3 .- 3/2*usquared[i]/c^2)
    end

    feq = zeros((xres*yres, 9))
    for i=1:xres*yres
        feq[i,:] = w[:].*rho[i] .+ rho[i].*s[i,:]
    end

    lambda = zeros(xres*yres)
    fnext = zeros((xres*yres, 9))
    for i=1:xres*yres
        lambda[i] = (1 - K1[i]*deltat/tau + K2[i]*deltat)
        fnext[i,:] = (fnew[i,:] .- K1[i]*(deltat/tau).*feq[i,:] .+ K2[i]*deltat.*fb[i,:])./lambda[i]
    end

    rhoepsilon = zeros(xres*yres)
    for i=1:xres*yres
        rhoepsilon[i] = sum(gnew[i,:] .- deltat/2 .*fnext[i,:].*q[i,:])
    end

    geq = zeros((xres*yres, 9))
    for i=1:xres*yres
        geq[i,1] = -2/3*rhoepsilon[i] * usquared[i]/c^2
        geq[i,2] = rhoepsilon[i]/9 * (1.5 + 1.5*dot(unitvecs[2,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[2,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
        geq[i,3] = rhoepsilon[i]/9 * (1.5 + 1.5*dot(unitvecs[3,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[3,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
        geq[i,4] = rhoepsilon[i]/9 * (1.5 + 1.5*dot(unitvecs[4,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[4,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
        geq[i,5] = rhoepsilon[i]/9 * (1.5 + 1.5*dot(unitvecs[5,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[5,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
        geq[i,6] = rhoepsilon[i]/36 * (3 + 6*dot(unitvecs[6,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[6,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
        geq[i,7] = rhoepsilon[i]/36 * (3 + 6*dot(unitvecs[7,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[7,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
        geq[i,8] = rhoepsilon[i]/36 * (3 + 6*dot(unitvecs[8,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[8,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
        geq[i,9] = rhoepsilon[i]/36 * (3 + 6*dot(unitvecs[9,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[9,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
    end

    gnext = zeros(xres*yres, 9)
    for i=1:xres*yres
        gnext[i,:] = (gnew[i,:] .- K1[i]*deltat/tauc.*geq[i,:] .+ K1[i]*deltat.*fnext[i,:].*q[i,:] .+ K2[i]*deltat.*gb[i,:])./lambda[i]
    end

    return fnext, gnext, q, rho, u, feq, geq
end

function stepclamp(x, lowerlim, upperlim)
    if x < lowerlim
        x = lowerlim
    elseif x > upperlim
        x = upperlim
    end
    return x
end

function stepsmooth(innerrad, outerrad, radius)
    x = stepclamp((radius-innerrad)/(outerrad-innerrad), 0.0, 1.0)
    return x
end

N = 20
tau = 1000
tauc = 1000
deltax = 1
deltat = 1
c = deltax/deltat
nu = c^2/3
xres = 100
yres = 100
periodictensor = PeriodicStreamTensor(xres, yres)
innerrad = .5
outerrad = .8
#K2 = 1 - K1 times some large (~100?) constant
#K1 = 1 in domain of interest
#step function in radius - 1 for values within domain of interest
#circle of radius res/3 centered at (xres/2,yres/2)

radius = zeros(xres,yres)
K1 = zeros(xres*yres)
for j=1:yres
    for i=1:xres
        radius[i,j] = sqrt(((i-50)/50)^2 + ((j-50)/50)^2)
        K1[(j-1)*xres+i] = 1 - stepsmooth(innerrad, outerrad, radius[i,j])
    end
end
K1trans = zeros(xres, yres)
for j=1:yres
    for i=1:xres
        K1trans[i,j] = K1[(j-1)*xres+i]
    end
end

K2 = zeros(xres*yres)
K2[:] .= (1 .- K1[:]) .* 1000

K2trans = zeros(xres, yres)
for j=1:yres
    for i=1:xres
        K2trans[i,j] = K2[(j-1)*xres+i]
    end
end

sx = 5
sy = 5
dx = 1
dy = 1
pDx = periodic_derivative(1, 0, sx, sy, xres, yres, dx, dy)
pDy = periodic_derivative(0, 1, sx, sy, xres, yres, dx, dy)
pDxx = periodic_derivative(2, 0, sx, sy, xres, yres, dx, dy)
pDyy = periodic_derivative(0, 2, sx, sy, xres, yres, dx, dy)

fsine = zeros(xres, yres, 9)
gsine = zeros(xres, yres, 9)
rhoinit = zeros(xres, yres)
uinit = zeros(xres, yres, 9)
for i=1:xres
    for j=1:yres
        for k=1:9
            fsine[i,j,k] = .5 + .4*sin(i*pi/(xres))*sin(j*pi/(yres))
        end
        rhoinit[i,j] = sum(fsine[i,j,:])
        uinit[i,j,:] = fsine[i,j,:] ./ rhoinit[i,j]
        for k=1:9
            gsine[i,j,k] = (1 - uinit[i,j,k])^2/2 * fsine[i,j,k]
        end
    end
end

fs2 = zeros(xres*yres, 9)
gs2 = zeros(xres*yres, 9)

for j=1:yres
    for i=1:xres
        fs2[(j-1)*xres+i,:] = fsine[i,j,:]
        gs2[(j-1)*xres+i,:] = gsine[i,j,:]
    end
end

fb = zeros(xres*yres, 9)
gb = zeros(xres*yres, 9)

rho = zeros((xres*yres))
for i=1:xres*yres
    rho[i] = sum(fs2[i,:])
end
p = zeros(xres*yres)
for i=1:xres*yres
    p[i] = c^2/3 * rho[i]
end
u = zeros(xres*yres, 2)
for i=1:xres*yres
    ux = (fs2[i,2] - fs2[i,4] + 1/sqrt(2) * (fs2[i,6] + fs2[i,9] - fs2[i,7] - fs2[i,8])) / rho[i]
    uy = (fs2[i,3] - fs2[i,5] + 1/sqrt(2) * (fs2[i,6] + fs2[i,8] - fs2[i,9] - fs2[i,8])) / rho[i]
    u[i,1] = ux
    u[i,2] = uy
end
unitvecs = zeros(9,2)
unitvecs[:,1] = [0, c, 0, -c, 0, 1/sqrt(2)*c, -1/sqrt(2)*c, -1/sqrt(2)*c, 1/sqrt(2)*c]
unitvecs[:,2] = [0, 0, c, 0, -c, 1/sqrt(2)*c, 1/sqrt(2)*c, -1/sqrt(2)*c, -1/sqrt(2)*c]
w = zeros(9)
w[1] = 4/9
w[2] = 1/9
w[3] = 1/9
w[4] = 1/9
w[5] = 1/9
w[6] = 1/36
w[7] = 1/36
w[8] = 1/36
w[9] = 1/36
s = zeros((xres*yres, 9))
usquared = zeros((xres*yres))
for i=1:xres*yres
    usquared[i] = dot(u[i,:],u[i,:])
    s[i,:] = w[:] .* (3*(unitvecs*u[i,:])/c^2 .+ 9/2*((unitvecs*u[i,:])).^2/c^3 .- 3/2*usquared[i]/c^2)
end
feqb = zeros((xres*yres, 9))
for i=1:xres*yres
    feqb[i,:] = w[:].*rho[i] .+ rho[i].*s[i,:]
end
q = zeros(xres*yres, 9)

grad = [pDx, pDy]
divu = pDx*u[:,1] + pDy*u[:,2]
gradp = zeros(xres*yres, 2)
gradp[:,1] = pDx*p[:]
gradp[:,2] = pDy*p[:]
laplu = zeros(xres*yres, 2)
laplu[:,1] = (pDxx+pDyy)*u[:,1]
laplu[:,2] = (pDxx+pDyy)*u[:,2]

vgradu = zeros(xres*yres, 9, 2)
for k=1:9
    vgradu[:,k,1] = (unitvecs[k,1] .- u[:,1]) .* (pDx*u[:,1]) + (unitvecs[k,2] .- u[:,2]) .* (pDy*u[:,1])
    vgradu[:,k,2] = (unitvecs[k,1] .- u[:,1]) .* (pDx*u[:,2]) + (unitvecs[k,2] .- u[:,2]) .* (pDy*u[:,2])
end

for i=1:xres*yres
    for a=1:9
        q[i,a] = - 1/rho[i]*dot(unitvecs[a,:]-u[i,:], gradp[i,:])
        + nu*(dot(unitvecs[a,:]-u[i,:], laplu[i,:]))
        + dot(unitvecs[a,:]-u[i,:], vgradu[i,a,:])
    end
end
lambda = zeros(xres*yres)
rhoepsilon = zeros(xres*yres)
for i=1:xres*yres
    rhoepsilon[i] = sum(gs2[i,:] .- deltat/2 .*fs2[i,:].*q[i,:])
end
geq = zeros((xres*yres, 9))
for i=1:xres*yres
    geq[i,1] = -2/3*rhoepsilon[i] * usquared[i]/c^2
    geq[i,2] = rhoepsilon[i]/9 * (1.5 + 1.5*dot(unitvecs[2,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[2,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
    geq[i,3] = rhoepsilon[i]/9 * (1.5 + 1.5*dot(unitvecs[3,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[3,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
    geq[i,4] = rhoepsilon[i]/9 * (1.5 + 1.5*dot(unitvecs[4,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[4,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
    geq[i,5] = rhoepsilon[i]/9 * (1.5 + 1.5*dot(unitvecs[5,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[5,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
    geq[i,6] = rhoepsilon[i]/36 * (3 + 6*dot(unitvecs[6,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[6,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
    geq[i,7] = rhoepsilon[i]/36 * (3 + 6*dot(unitvecs[7,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[7,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
    geq[i,8] = rhoepsilon[i]/36 * (3 + 6*dot(unitvecs[8,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[8,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
    geq[i,9] = rhoepsilon[i]/36 * (3 + 6*dot(unitvecs[9,:],u[i,:])/c^2 + 4.5*(dot(unitvecs[9,:],u[i,:]))^2/c^4 - 1.5*usquared[i]/c^2)
end
fb .= feqb
gb .= geq

anim = @animate for i=1:N
    fs2, gs2, ex, ex1, ex2 = boundary_thermal_lattice_boltzmann_method(periodictensor, pDx, pDy, pDxx, pDyy, xres, yres, fs2, gs2, fb, gb, K1, K2, tau, tauc, deltax, deltat)
    fplot = zeros(xres, yres)
    gplot = zeros(xres, yres)
    explot = zeros(xres, yres)
    exp1 = zeros(xres, yres)
    exp2 = zeros(xres, yres)
    ucheck = zeros(xres, yres)
    for i=1:xres
        for j=1:yres
            exp1[i,j] = ex1[(j-1)*xres+i]
            exp2[i,j] = ex2[(j-1)*xres+i,2]
            for k=1:9
                fplot[i,j] = fplot[i,j] + fs2[(j-1)*xres+i,k]
                gplot[i,j] = gplot[i,j] + gs2[(j-1)*xres+i,k]
                explot[i,j] = explot[i,j] + ex[(j-1)*xres+i,k]
                ucheck[i,j] = ucheck[i,j] + ex[(j-1)*xres+i,k]/ex1[(j-1)*xres+i]
            end
        end
    end
    println(fplot[Int32(xres/2),Int32(yres/2)])
    println(gplot[Int32(xres/2),Int32(yres/2)])
    heatmap(gplot, clims=(1,4))
end
gif(anim, fps=10)
# map = heatmap(K2trans)
# display(map)


text:
using Test

using DataStructures
using GeoInterface
using LightGraphs

using PlanarVisibility
const PV = PlanarVisibility

@testset "point sets" begin
    # empty set
    set = push!(PV.PointSet(), [])
    @test length(set) == 0

    # add individual points
    set = push!(PV.PointSet(), Point([0.0, 0.0]))
    @test length(set) == 1
    push!(set, Point([0.0, 0.0]))
    @test length(set) == 1
    push!(set, Point([0.0, 1.0]))
    @test length(set) == 2

    # add positions
    set = push!(PV.PointSet(), [0.0, 0.0])
    @test length(set) == 1

    set = push!(PV.PointSet(), [0.0, 0.0, 0.0])
    @test length(set) == 1

    # add line string coordinates
    set = push!(PV.PointSet(), [[0.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
    @test length(set) == 3

    # add polygon coordinates
    set = push!(PV.PointSet(), [[[0.0, 0.0], [0.0, 1.0], [1.0, 0.0]]])
    @test length(set) == 3

    # add polygon with hole
    set = push!(PV.PointSet(),
                Polygon([[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [0.0, 0.0]],
                         [[0.1, 0.1], [0.1, 0.9], [0.9, 0.1], [0.1, 0.1]]]))
    @test length(set) == 6

    # test getindex
    @test set[1] == Point([0.0, 0.0])
    @test set[Point([0.0, 0.0])] == 1
    @test set[[0.0, 0.0]] == 1
    @test set[4] == Point([0.1, 0.1])
    @test set[Point([0.1, 0.1])] == 4
    @test set[[0.1, 0.1]] == 4
end

@testset "extract_points" begin
    @testset "empty environment has no points" begin
        env = PV.Environment([])
        set = PV.extract_points(env)
        @test isempty(set.points)
        @test isempty(set.indices)
    end

    @testset "single polygon" begin
        env = PV.Environment([
            Polygon([[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [0.0, 0.0]]])
        ])
        set = PV.extract_points(env)
        @test length(set.points) == 3
        @test length(set.indices) == 3
        @test set.points[1] == Point([0.0, 0.0])
        @test all(set.indices[set.points[i]] == i for i in 1:3)
    end

    @testset "polygon with hole" begin
        env = PV.Environment([
            Polygon([[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [0.0, 0.0]],
                     [[0.1, 0.1], [0.1, 0.9], [0.9, 0.1], [0.1, 0.1]]])
        ])
        set = PV.extract_points(env)
        @test length(set.points) == 6
        @test length(set.indices) == 6
        @test set.points[1] == Point([0.0, 0.0])
        @test all(set.indices[set.points[i]] == i for i in 1:6)
    end

    @testset "two touching polygons" begin
        env = PV.Environment([
            Polygon([
                [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [0.0, 0.0]],
                [[1.0, 0.0], [2.0, 0.0], [2.0, 1.0], [1.0, 1.0], [1.0, 0.0]]
            ])
        ])
        set = PV.extract_points(env)
        @test length(set.points) == 6
        @test length(set.indices) == 6
        @test set.points[1] == Point([0.0, 0.0])
        @test set.points[6] == Point([2.0, 1.0])
        @test all(set.indices[set.points[i]] == i for i in 1:6)
    end
end

@testset "extract_edges" begin
    # empty environment
    env = PV.Environment([])
    set = PV.extract_points(env)
    graph = PV.extract_edges(env, set)
    @test nv(graph) == 0
    @test ne(graph) == 0

    # polygon with hole
    env = PV.Environment([
        Polygon([[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [0.0, 0.0]],
                 [[0.1, 0.1], [0.1, 0.9], [0.9, 0.1], [0.1, 0.1]]])])
    set = PV.extract_points(env)
    graph = PV.extract_edges(env, set)
    @test nv(graph) == 6
    @test ne(graph) == 6
    @test neighbors(graph, 1) == [2, 3]
    @test neighbors(graph, 2) == [1, 3]
    @test neighbors(graph, 3) == [1, 2]
    @test neighbors(graph, 4) == [5, 6]
    @test neighbors(graph, 5) == [4, 6]
    @test neighbors(graph, 6) == [4, 5]
end

@testset "extract_components" begin
    # two polygons
    env = PV.Environment([
        Polygon([[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [0.0, 0.0]]]),
        Polygon([[[2.0, 0.0], [3.0, 0.0], [2.0, 1.0], [2.0, 0.0]]])])
    set = PV.extract_points(env)
    comps = PV.extract_components(env, set)
    @test in_same_set(comps, 1, 2)
    @test in_same_set(comps, 1, 3)
    @test in_same_set(comps, 2, 3)
    @test !in_same_set(comps, 1, 4)
    @test !in_same_set(comps, 1, 5)
    @test !in_same_set(comps, 1, 6)
    @test in_same_set(comps, 4, 5)
    @test in_same_set(comps, 5, 6)
    @test in_same_set(comps, 6, 4)
end

@testset "compute angles" begin
    # 4  3
    # 1  2
    p1, p2, p3, p4 = Point.([[0, 0], [1, 0], [1, 1], [0, 1]])
    @test PV.angle_to_east(p1, p1) == 0.0
    @test PV.angle_to_east(p1, p2) ≈ 0.0
    @test PV.angle_to_east(p1, p3) ≈ 2π/8
    @test PV.angle_to_east(p1, p4) ≈ 2π/4
    @test PV.angle_to_east(p2, p1) ≈ 2π/2
    @test PV.angle_to_east(p2, p3) ≈ 2π/4
    @test PV.angle_to_east(p2, p4) ≈ 3/8 * 2π
    @test PV.angle_to_east(p3, p1) ≈ 5/8 * 2π
    @test PV.angle_to_east(p3, p2) ≈ 3/4 * 2π
    @test PV.angle_to_east(p3, p4) ≈ 2π/2
    @test PV.angle_to_east(p4, p1) ≈ 3/4 * 2π
    @test PV.angle_to_east(p4, p2) ≈ 7/8 * 2π
    @test PV.angle_to_east(p4, p3) ≈ 0.0
end

@testset "sorting points in counter clockwise order" begin
    #  4  3
    #  1  2
    p1, p2, p3, p4 = Point.([[0, 0], [1, 0], [1, 1], [0, 1]])

    @test PV.sortperm_ccw([p2, p3, p4], p1) == [1, 2, 3]
    @test PV.sortperm_ccw([p4, p2, p3], p1) == [2, 3, 1]
    @test PV.sortperm_ccw([p3, p2, p1], p4) == [1, 3, 2]
    @test PV.sortperm_ccw([p1, p2, p3, p4], p1) == [1, 2, 3, 4]

    #  1 2
    # 3   4
    #  5 6
    q1, q2, q3, q4, q5, q6 = Point.(
        [[1., 2.], [2., 2.], [0., 1.], [3., 1.], [1., 0.], [2., 0.]])

    @test PV.sortperm_ccw([q1, q2, q3, q4, q5, q6], q1) == [1, 2, 3, 5, 6, 4]
    @test PV.sortperm_ccw([q1, q2, q3, q4, q5, q6], q3) == [3, 4, 2, 1, 5, 6]
    @test PV.sortperm_ccw([q1, q2, q3, q4, q5, q6], q5) == [5, 6, 4, 2, 1, 3]
end

@testset "sorted edges" begin
    @test PV.sorted(Edge(1, 2)) == Edge(1, 2)
    @test PV.sorted(Edge(1, 1)) == Edge(1, 1)
    @test PV.sorted(Edge(2, 1)) == Edge(1, 2)
end

@testset "triangle orientation" begin
    # p -- q -- r
    p, q, r = Point.([[0.0, 0.0], [1.0, 2.0], [2.0, 4.0]])
    @test PV.orientation(p, q, r) == PV.COLLINEAR
    @test PV.orientation(r, q, p) == PV.COLLINEAR
    @test PV.orientation(q, r, p) == PV.COLLINEAR
    @test PV.orientation(r, p, q) == PV.COLLINEAR
    @test PV.orientation(q, p, r) == PV.COLLINEAR
    @test PV.orientation(p, r, q) == PV.COLLINEAR

    #   ,q.
    # p'---`r
    p, q, r = Point.([[0.0, 0.0], [1.0, 1.0], [2.0, 0.0]])
    @test PV.orientation(p, q, r) == PV.CLOCKWISE
    @test PV.orientation(r, q, p) == PV.COUNTERCLOCKWISE
    @test PV.orientation(q, r, p) == PV.CLOCKWISE
    @test PV.orientation(r, p, q) == PV.CLOCKWISE
    @test PV.orientation(q, p, r) == PV.COUNTERCLOCKWISE
    @test PV.orientation(p, r, q) == PV.COUNTERCLOCKWISE
end

@testset "collinear intersection" begin
    # p -- q -- r
    p, q, r = Point.([[0.0, 0.0], [1.0, 2.0], [2.0, 4.0]])
    @test PV.collinear_intersect(p, q, r) == true
    @test PV.collinear_intersect(r, q, p) == true
    @test PV.collinear_intersect(q, r, p) == false
    @test PV.collinear_intersect(r, p, q) == false
    @test PV.collinear_intersect(q, p, r) == false
    @test PV.collinear_intersect(p, r, q) == false
end

@testset "intersecting segments" begin
    # \/
    # /\
    points = Point.([[0., 0.], [1., 1.], [1., 0.], [0., 1.]])
    @test PV.intersect_segments(points...) == true
    # | |
    # | |
    points = Point.([[0., 0.], [0., 1.], [1., 0.], [1., 1.]])
    @test PV.intersect_segments(points...) == false
    # |
    # |__
    points = Point.([[0., 0.], [0., 1.], [0., 0.], [1., 0.]])
    @test PV.intersect_segments(points...) == true
    # |
    # | __
    points = Point.([[0., 0.], [0., 1.], [1., 0.], [2., 0.]])
    @test PV.intersect_segments(points...) == false
    # -- --
    points = Point.([[0., 0.], [0., 2.], [0., 3.], [0., 5.]])
    @test PV.intersect_segments(points...) == false
    # -=-  (parallel but intersecting)
    points = Point.([[0., 0.], [0., 2.], [0., 1.], [0., 5.]])
    @test PV.intersect_segments(points...) == true
end

@testset "visible points -- lines" begin
    #  1-2
    # 3---4
    #  5-6
    points = PV.PointSet()
    push!(points, [[1.0, 2.0], [2.0, 2.0],
                   [0.0, 1.0], [3.0, 1.0],
                   [1.0, 0.0], [2.0, 0.0]])
    @test length(points) == 6

    graph = SimpleGraph(6)
    add_edge!(graph, 1, 2)
    add_edge!(graph, 3, 4)
    add_edge!(graph, 5, 6)

    pc = IntDisjointSets(6) # no polygons

    @test PV.visible_points(points, graph, pc, 1) == [2, 3, 4]
    @test PV.visible_points(points, graph, pc, 2) == [1, 3, 4]
    @test PV.visible_points(points, graph, pc, 3) == [4, 2, 1, 5, 6]
    @test PV.visible_points(points, graph, pc, 4) == [2, 1, 3, 5, 6]
    @test PV.visible_points(points, graph, pc, 5) == [6, 4, 3]
    @test PV.visible_points(points, graph, pc, 6) == [4, 3, 5]
end

@testset "visible points -- boxes" begin
    #      8--7
    # 4--3 |  |
    # |  | 5--6
    # 1--2
    env = PV.Environment([
        Polygon([[[0., 0.], [1., 0.], [1., 2.], [0., 2.], [0., 0.]]]),
        Polygon([[[2., 1.], [3., 1.], [4., 3.], [2., 3.], [2., 1.]]])])
    points = PV.extract_points(env)
    graph = PV.extract_edges(env, points)
    comps = PV.extract_components(env, points)

    @test length(points) == 8
    @test nv(graph) == 8
    @test ne(graph) == 8

    vis = i -> PV.visible_points(points, graph, comps, i)
    @test vis(1) == [2, 4]
    @test vis(2) == [6, 5, 8, 3, 1]
    @test vis(3) == [8, 4, 2, 5]
    @test vis(4) == [3, 8, 1]
end


text:
using Statistics
using JLD2
using Printf
using Plots
using Oceananigans
using Oceananigans.Units

using Oceananigans.Fields: interpolate, Field
using Oceananigans.Architectures: arch_array
using Oceananigans.Coriolis: HydrostaticSphericalCoriolis
using Oceananigans.BoundaryConditions
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom, inactive_node, peripheral_node
using CUDA: @allowscalar, device!
using Oceananigans.Operators
using Oceananigans.Operators: Δzᵃᵃᶜ
using Oceananigans: prognostic_fields

using Oceananigans.Models.HydrostaticFreeSurfaceModels: VerticalVorticityField

@inline function visualize(field, lev, dims)
    (dims == 1) && (idx = (lev, :, :))
    (dims == 2) && (idx = (:, lev, :))
    (dims == 3) && (idx = (:, :, lev))

    r = deepcopy(Array(interior(field)))[idx...]
    r[ r.==0 ] .= NaN
    return r
end

#####
##### Grid
#####

arch = GPU()
reference_density = 1029

latitude = (-75, 75)

# 0.25 degree resolution
Nx = 1440
Ny = 600
Nz = 1

const Nyears  = 1
const Nmonths = 12
const thirty_days = 30days

output_prefix = "near_global_lat_lon_$(Nx)_$(Ny)_$(Nz)_fine"
pickup_file   = false

#####
##### Load forcing files and inital conditions from ECCO version 4
##### https://ecco.jpl.nasa.gov/drive/files
##### Bathymetry is interpolated from ETOPO1 https://www.ngdc.noaa.gov/mgg/global/
#####

using DataDeps

path = "https://github.com/CliMA/OceananigansArtifacts.jl/raw/ss/new_hydrostatic_data_after_cleared_bugs/quarter_degree_near_global_input_data/"

datanames = ["z_faces-50-levels",
             "bathymetry-1440x600",
             "temp-1440x600-latitude-75",
             "salt-1440x600-latitude-75",
             "tau_x-1440x600-latitude-75",
             "tau_y-1440x600-latitude-75",
             "initial_conditions"]

dh = DataDep("quarter_degree_near_global_lat_lon",
    "Forcing data for global latitude longitude simulation",
    [path * data * ".jld2" for data in datanames]
)

DataDeps.register(dh)

datadep"quarter_degree_near_global_lat_lon"

files = [:file_z_faces, :file_bathymetry, :file_temp, :file_salt, :file_tau_x, :file_tau_y, :file_init]
for (data, file) in zip(datanames, files)
    datadep_path = @datadep_str "quarter_degree_near_global_lat_lon/" * data * ".jld2"
    @eval $file = jldopen($datadep_path)
end

bathymetry = file_bathymetry["bathymetry"]
bathymetry[bathymetry .< 0] .= -10e3

τˣ = zeros(Nx, Ny, Nmonths)
τʸ = zeros(Nx, Ny, Nmonths)

# Files contain 1 year (1992) of 12 monthly averages
τˣ = file_tau_x["field"] ./ reference_density
τʸ = file_tau_y["field"] ./ reference_density

# Remember the convention!! On the surface a negative flux increases a positive decreases
bathymetry = arch_array(arch, bathymetry)

# A spherical domain
@show underlying_grid = LatitudeLongitudeGrid(arch,
                                              size = (Nx, Ny, Nz),
                                              longitude = (-180, 180),
                                              latitude = latitude,
                                              halo = (3, 3, 3),
                                              z = (-10e3, 0),
                                              #z = z_faces,
                                              precompute_metrics = true)

grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bathymetry))

τˣ = arch_array(arch, - τˣ)
τʸ = arch_array(arch, - τʸ)

#####
##### Boundary conditions / time-dependent fluxes
#####

@inline current_time_index(time, tot_months)     = mod(unsafe_trunc(Int32, time / thirty_days), tot_months) + 1
@inline next_time_index(time, tot_months)        = mod(unsafe_trunc(Int32, time / thirty_days) + 1, tot_months) + 1
@inline cyclic_interpolate(u₁::Number, u₂, time) = u₁ + mod(time / thirty_days, 1) * (u₂ - u₁)

Δz_top    = @allowscalar Δzᵃᵃᶜ(1, 1, grid.Nz, grid.underlying_grid)
Δz_bottom = @allowscalar Δzᵃᵃᶜ(1, 1, 1, grid.underlying_grid)

@inline function surface_wind_stress(i, j, grid, clock, fields, τ)
    time = clock.time
    n₁ = current_time_index(time, Nmonths)
    n₂ = next_time_index(time, Nmonths)

    @inbounds begin
        τ₁ = τ[i, j, n₁]
        τ₂ = τ[i, j, n₂]
    end

    return cyclic_interpolate(τ₁, τ₂, time)
end

u_wind_stress_bc = FluxBoundaryCondition(surface_wind_stress, discrete_form = true, parameters = τˣ)
v_wind_stress_bc = FluxBoundaryCondition(surface_wind_stress, discrete_form = true, parameters = τʸ)

@inline u_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.u[i, j, 1]
@inline v_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.v[i, j, 1]

# Linear bottom drag:
μ = 0.001 # ms⁻¹

@inline is_immersed_drag_u(i, j, k, grid) = Int(peripheral_node(Face(), Center(), Center(), i, j, k-1, grid) & !inactive_node(Face(), Center(), Center(), i, j, k, grid))
@inline is_immersed_drag_v(i, j, k, grid) = Int(peripheral_node(Center(), Face(), Center(), i, j, k-1, grid) & !inactive_node(Center(), Face(), Center(), i, j, k, grid))

# Keep a constant linear drag parameter independent on vertical level
@inline u_immersed_bottom_drag(i, j, k, grid, clock, fields, μ) = @inbounds - μ * is_immersed_drag_u(i, j, k, grid) * fields.u[i, j, k] / Δzᵃᵃᶜ(i, j, k, grid)
@inline v_immersed_bottom_drag(i, j, k, grid, clock, fields, μ) = @inbounds - μ * is_immersed_drag_v(i, j, k, grid) * fields.v[i, j, k] / Δzᵃᵃᶜ(i, j, k, grid)

Fu = Forcing(u_immersed_bottom_drag, discrete_form = true, parameters = μ)
Fv = Forcing(v_immersed_bottom_drag, discrete_form = true, parameters = μ)

u_bottom_drag_bc = FluxBoundaryCondition(u_bottom_drag, discrete_form = true, parameters = μ)
v_bottom_drag_bc = FluxBoundaryCondition(v_bottom_drag, discrete_form = true, parameters = μ)

u_bcs = FieldBoundaryConditions(top = u_wind_stress_bc, bottom = u_bottom_drag_bc)
v_bcs = FieldBoundaryConditions(top = v_wind_stress_bc, bottom = v_bottom_drag_bc)

free_surface = ImplicitFreeSurface(solver_method=:HeptadiagonalIterativeSolver)
# free_surface = ExplicitFreeSurface()

buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState())

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    free_surface = free_surface,
                                    momentum_advection = VectorInvariant(),
                                    coriolis = HydrostaticSphericalCoriolis(),
                                    buoyancy = nothing,
                                    boundary_conditions = (u=u_bcs, v=v_bcs))

#####
##### Initial condition:
#####

u, v, w = model.velocities
η = model.free_surface.η

@info "model initialized"

#####
##### Simulation setup
#####

ζ = VerticalVorticityField(model)
compute!(ζ)

Δt = 6minutes  # for initialization, then we can go up to 6 minutes?
simulation = Simulation(model, Δt = Δt, stop_time = Nyears*years)
start_time = [time_ns()]

using Oceananigans.Utils

function progress(sim)
    wall_time = (time_ns() - start_time[1]) * 1e-9

    u = sim.model.velocities.u
    η = sim.model.free_surface.η

    @info @sprintf("Time: % 12s, iteration: %d, max(|u|): %.2e ms⁻¹, wall time: %s",
                    prettytime(sim.model.clock.time),
                    sim.model.clock.iteration, maximum(abs, u),
                    prettytime(wall_time))

    start_time[1] = time_ns()

    return nothing
end

simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))

u, v, w = model.velocities
η = model.free_surface.η

output_fields = merge(model.velocities, (η=model.free_surface.η, ζ=ζ))
save_interval = 1days

simulation.output_writers[:surface_fields] = JLD2OutputWriter(model, output_fields, #(; u, v, T, S, η),
                                                              schedule = TimeInterval(save_interval),
                                                              filename = output_prefix * "_surface",
                                                              indices = (:, :, grid.Nz),
                                                              overwrite_existing = true)

simulation.output_writers[:checkpointer] = Checkpointer(model,
                                                        schedule = TimeInterval(1day),
                                                        prefix = output_prefix * "_checkpoint",
                                                        overwrite_existing = true)

@info "Running with Δt = $(prettytime(simulation.Δt))"

run!(simulation, pickup = pickup_file)

@info """

    Simulation took $(prettytime(simulation.run_wall_time))
    Free surface: $(typeof(model.free_surface).name.wrapper)
    Time step: $(prettytime(Δt))
"""


text:
note!(R::SeisSrc, s::String) = push!(R.notes, tnote(s))
note!(H::SeisHdr, s::String) = push!(H.notes, tnote(s))
clear_notes!(U::Union{SeisHdr,SeisSrc}) = (U.notes = Array{String,1}(undef,1); U.notes[1] = tnote("notes cleared."); return nothing)
clear_notes!(Ev::SeisEvent) = (clear_notes!(Ev.hdr); clear_notes!(Ev.source); clear_notes!(Ev.data))