text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -47.02611826463841 <= q <= 46.175694332378946)

                     add_NL_constraint(m, :(sigmoid(0.14172886209520552 + -0.7326244672203441*sigmoid(-0.9862473619768446 + 0.12598046880516867*sigmoid(-0.32432918385015697 + -0.6977940842525037*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.476420397288412*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.5549201217740469*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.4315951265237059*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.9857820194451663*sigmoid(0.2633379906083797 + 0.3509284626355762*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.5977055959243147*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.642253969451728*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.10048605306401859*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.3507087497777257*sigmoid(-0.34549165187553976 + 0.5367448324029058*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.7077894233575681*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.4911236080946888*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + -0.5306931393323575*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.16930171271231975*sigmoid(0.538919292477293 + 0.8088100362525972*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.4318478690072691*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + 0.8077402924106214*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.12603114956525063*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2])))) + 0.6188053648075824*sigmoid(-0.8901994655270467 + 0.3517572387180481*sigmoid(-0.32432918385015697 + -0.6977940842525037*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.476420397288412*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.5549201217740469*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.4315951265237059*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.5057680398491988*sigmoid(0.2633379906083797 + 0.3509284626355762*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.5977055959243147*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.642253969451728*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.10048605306401859*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.15346000401268212*sigmoid(-0.34549165187553976 + 0.5367448324029058*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.7077894233575681*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.4911236080946888*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + -0.5306931393323575*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.000907155550359029*sigmoid(0.538919292477293 + 0.8088100362525972*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.4318478690072691*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + 0.8077402924106214*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.12603114956525063*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2])))) + 0.748255673773524*sigmoid(0.9418689270026719 + 0.23440750676702393*sigmoid(-0.32432918385015697 + -0.6977940842525037*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.476420397288412*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.5549201217740469*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.4315951265237059*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + -0.9578891219409944*sigmoid(0.2633379906083797 + 0.3509284626355762*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.5977055959243147*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.642253969451728*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.10048605306401859*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + -0.4332120963929875*sigmoid(-0.34549165187553976 + 0.5367448324029058*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.7077894233575681*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.4911236080946888*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + -0.5306931393323575*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.15216088947964046*sigmoid(0.538919292477293 + 0.8088100362525972*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.4318478690072691*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + 0.8077402924106214*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.12603114956525063*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2])))) + -0.19681615050825796*sigmoid(-0.3160427191298667 + -0.1453900085539872*sigmoid(-0.32432918385015697 + -0.6977940842525037*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.476420397288412*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.5549201217740469*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.4315951265237059*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.3370407689400481*sigmoid(0.2633379906083797 + 0.3509284626355762*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.5977055959243147*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.642253969451728*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.10048605306401859*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.5646332930109117*sigmoid(-0.34549165187553976 + 0.5367448324029058*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.7077894233575681*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.4911236080946888*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + -0.5306931393323575*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + -0.4312398718758548*sigmoid(0.538919292477293 + 0.8088100362525972*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.4318478690072691*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + 0.8077402924106214*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.12603114956525063*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))))) + sigmoid(0.7126847484534728 + 0.5521736820641339*sigmoid(-0.9862473619768446 + 0.12598046880516867*sigmoid(-0.32432918385015697 + -0.6977940842525037*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.476420397288412*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.5549201217740469*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.4315951265237059*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.9857820194451663*sigmoid(0.2633379906083797 + 0.3509284626355762*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.5977055959243147*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.642253969451728*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.10048605306401859*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.3507087497777257*sigmoid(-0.34549165187553976 + 0.5367448324029058*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.7077894233575681*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.4911236080946888*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + -0.5306931393323575*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.16930171271231975*sigmoid(0.538919292477293 + 0.8088100362525972*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.4318478690072691*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + 0.8077402924106214*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.12603114956525063*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2])))) + -0.5863621462894022*sigmoid(-0.8901994655270467 + 0.3517572387180481*sigmoid(-0.32432918385015697 + -0.6977940842525037*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.476420397288412*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.5549201217740469*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.4315951265237059*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.5057680398491988*sigmoid(0.2633379906083797 + 0.3509284626355762*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.5977055959243147*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.642253969451728*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.10048605306401859*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.15346000401268212*sigmoid(-0.34549165187553976 + 0.5367448324029058*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.7077894233575681*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.4911236080946888*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + -0.5306931393323575*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.000907155550359029*sigmoid(0.538919292477293 + 0.8088100362525972*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.4318478690072691*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + 0.8077402924106214*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.12603114956525063*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2])))) + -0.2365722039484086*sigmoid(0.9418689270026719 + 0.23440750676702393*sigmoid(-0.32432918385015697 + -0.6977940842525037*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.476420397288412*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.5549201217740469*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.4315951265237059*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + -0.9578891219409944*sigmoid(0.2633379906083797 + 0.3509284626355762*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.5977055959243147*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.642253969451728*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.10048605306401859*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + -0.4332120963929875*sigmoid(-0.34549165187553976 + 0.5367448324029058*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.7077894233575681*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.4911236080946888*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + -0.5306931393323575*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.15216088947964046*sigmoid(0.538919292477293 + 0.8088100362525972*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.4318478690072691*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + 0.8077402924106214*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.12603114956525063*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2])))) + -0.47134077688752996*sigmoid(-0.3160427191298667 + -0.1453900085539872*sigmoid(-0.32432918385015697 + -0.6977940842525037*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.476420397288412*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.5549201217740469*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.4315951265237059*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.3370407689400481*sigmoid(0.2633379906083797 + 0.3509284626355762*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.5977055959243147*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.642253969451728*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.10048605306401859*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.5646332930109117*sigmoid(-0.34549165187553976 + 0.5367448324029058*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.7077894233575681*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.4911236080946888*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + -0.5306931393323575*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + -0.4312398718758548*sigmoid(0.538919292477293 + 0.8088100362525972*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.4318478690072691*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + 0.8077402924106214*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.12603114956525063*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))))) + sigmoid(0.3678622626620536 + -0.645924933237755*sigmoid(-0.9862473619768446 + 0.12598046880516867*sigmoid(-0.32432918385015697 + -0.6977940842525037*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.476420397288412*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.5549201217740469*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.4315951265237059*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.9857820194451663*sigmoid(0.2633379906083797 + 0.3509284626355762*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.5977055959243147*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.642253969451728*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.10048605306401859*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.3507087497777257*sigmoid(-0.34549165187553976 + 0.5367448324029058*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.7077894233575681*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.4911236080946888*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + -0.5306931393323575*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.16930171271231975*sigmoid(0.538919292477293 + 0.8088100362525972*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.4318478690072691*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + 0.8077402924106214*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.12603114956525063*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2])))) + 0.23980134898526861*sigmoid(-0.8901994655270467 + 0.3517572387180481*sigmoid(-0.32432918385015697 + -0.6977940842525037*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.476420397288412*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.5549201217740469*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.4315951265237059*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.5057680398491988*sigmoid(0.2633379906083797 + 0.3509284626355762*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.5977055959243147*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.642253969451728*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.10048605306401859*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.15346000401268212*sigmoid(-0.34549165187553976 + 0.5367448324029058*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.7077894233575681*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.4911236080946888*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + -0.5306931393323575*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.000907155550359029*sigmoid(0.538919292477293 + 0.8088100362525972*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.4318478690072691*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + 0.8077402924106214*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.12603114956525063*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2])))) + 0.5641271464877478*sigmoid(0.9418689270026719 + 0.23440750676702393*sigmoid(-0.32432918385015697 + -0.6977940842525037*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.476420397288412*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.5549201217740469*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.4315951265237059*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + -0.9578891219409944*sigmoid(0.2633379906083797 + 0.3509284626355762*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.5977055959243147*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.642253969451728*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.10048605306401859*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + -0.4332120963929875*sigmoid(-0.34549165187553976 + 0.5367448324029058*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.7077894233575681*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.4911236080946888*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + -0.5306931393323575*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.15216088947964046*sigmoid(0.538919292477293 + 0.8088100362525972*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.4318478690072691*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + 0.8077402924106214*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.12603114956525063*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2])))) + -0.22456183587218126*sigmoid(-0.3160427191298667 + -0.1453900085539872*sigmoid(-0.32432918385015697 + -0.6977940842525037*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.476420397288412*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.5549201217740469*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.4315951265237059*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.3370407689400481*sigmoid(0.2633379906083797 + 0.3509284626355762*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.5977055959243147*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.642253969451728*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.10048605306401859*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.5646332930109117*sigmoid(-0.34549165187553976 + 0.5367448324029058*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.7077894233575681*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.4911236080946888*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + -0.5306931393323575*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + -0.4312398718758548*sigmoid(0.538919292477293 + 0.8088100362525972*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.4318478690072691*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + 0.8077402924106214*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.12603114956525063*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))))) + sigmoid(-0.9861590919257504 + -0.9031060096581198*sigmoid(-0.9862473619768446 + 0.12598046880516867*sigmoid(-0.32432918385015697 + -0.6977940842525037*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.476420397288412*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.5549201217740469*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.4315951265237059*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.9857820194451663*sigmoid(0.2633379906083797 + 0.3509284626355762*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.5977055959243147*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.642253969451728*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.10048605306401859*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.3507087497777257*sigmoid(-0.34549165187553976 + 0.5367448324029058*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.7077894233575681*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.4911236080946888*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + -0.5306931393323575*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.16930171271231975*sigmoid(0.538919292477293 + 0.8088100362525972*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.4318478690072691*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + 0.8077402924106214*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.12603114956525063*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2])))) + -0.6492392304855641*sigmoid(-0.8901994655270467 + 0.3517572387180481*sigmoid(-0.32432918385015697 + -0.6977940842525037*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.476420397288412*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.5549201217740469*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.4315951265237059*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.5057680398491988*sigmoid(0.2633379906083797 + 0.3509284626355762*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.5977055959243147*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.642253969451728*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.10048605306401859*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.15346000401268212*sigmoid(-0.34549165187553976 + 0.5367448324029058*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.7077894233575681*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.4911236080946888*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + -0.5306931393323575*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.000907155550359029*sigmoid(0.538919292477293 + 0.8088100362525972*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.4318478690072691*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + 0.8077402924106214*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.12603114956525063*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2])))) + 0.1540584284169868*sigmoid(0.9418689270026719 + 0.23440750676702393*sigmoid(-0.32432918385015697 + -0.6977940842525037*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.476420397288412*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.5549201217740469*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.4315951265237059*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + -0.9578891219409944*sigmoid(0.2633379906083797 + 0.3509284626355762*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.5977055959243147*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.642253969451728*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.10048605306401859*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + -0.4332120963929875*sigmoid(-0.34549165187553976 + 0.5367448324029058*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.7077894233575681*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.4911236080946888*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + -0.5306931393323575*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.15216088947964046*sigmoid(0.538919292477293 + 0.8088100362525972*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.4318478690072691*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + 0.8077402924106214*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.12603114956525063*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2])))) + -0.8682585596517867*sigmoid(-0.3160427191298667 + -0.1453900085539872*sigmoid(-0.32432918385015697 + -0.6977940842525037*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.476420397288412*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.5549201217740469*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.4315951265237059*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.3370407689400481*sigmoid(0.2633379906083797 + 0.3509284626355762*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.5977055959243147*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.642253969451728*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.10048605306401859*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + 0.5646332930109117*sigmoid(-0.34549165187553976 + 0.5367448324029058*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.7077894233575681*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + -0.4911236080946888*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + -0.5306931393323575*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(x[2]))) + -0.4312398718758548*sigmoid(0.538919292477293 + 0.8088100362525972*sigmoid(-0.9860262540836815 + 0.8911186935536937*$(x[1]) + 0.64267061852934*$(x[2])) + -0.4318478690072691*sigmoid(-0.12607269491856643 + -0.6901389567897818*$(x[1]) + -0.5431979251870089*$(x[2])) + 0.8077402924106214*sigmoid(0.4882924486844056 + 0.5602141072539029*$(x[1]) + 0.7222626702909976*$(x[2])) + 0.12603114956525063*sigmoid(-0.9654866490330036 + -0.27186143772128135*$(x[1]) + -0.864410033347212*$(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:
module _

import Base.map
import Base.reduce

# Collection functions
function each(a::Array, iterator::Function)
    for item in a
        iterator(item)
    end
end

# Array functions
function difference(a::Array, rest...)
    rest = vcat(rest...)
    filter((value) -> !in(value, rest), a)
end

without(a::Array, values...) = difference(a, values...)

# Function functions
partial(func::Function, args...) = (bargs...) -> func(args..., bargs...)

function memoize(func::Function, hasher)
    memo = Dict()
    function(args...)
        key = hasher(args...)
        has(memo, key) ? memo[key] : (memo[key] = func(args...))
    end
end

function memoize(func::Function)
    memo = Dict()
    function(args...)
        key = identity(args...)
        has(memo, key) ? memo[key] : (memo[key] = func(args...))
    end
end

function compose(funcs::Function...)
    return function(args...)
        for func in reverse(funcs)
            args = func(args...)
        end
        args
    end
end

end

ext: jl

max_stars_repo_name: hmbeale/dictionary

max_issues_repo_name: hmbeale/dictionary

max_forks_repo_name: hmbeale/dictionary


text:
module Psychometrics
    using Statistics
    using LinearAlgebra

    include("reliability.jl")

    export cronalpha
end


text:
# # Finding patterns in signals using Optimal Transport and Dynamic Time Warping

# In this example we will compare different distances when used as the inner metric in a dynamic time-warping of spectrograms.
# Setup
using DynamicAxisWarping, SpectralDistances, LPVSpectral, DSP, Plots, Statistics, BenchmarkTools, Distances, Random, ThreadTools, AlphaStableDistributions, SlidingDistancesBase
theme(:default)
gr();
Random.seed!(0);

# We start by creating a short patter, the *query*, we then create a long time series that contains a quite similar sound, but where the frequencies of the chirp are slightly higher. This is a realistic scenario when doing, e.g., acoustic detection. Different individuals within the same animal species might have similar calls, but slightly different pitch etc. We also add some alpha sub-gaussian noise in order to make the problem a bit harder.

N       = 48_000 ÷ 2
nfft    = 4096
f1      = range(0.01, stop = 1, length = N)
f2      = range(1, stop = 0.1, length = N) .^ 2
g(x, N) = exp(-10 * (x - N / 2)^2 / N^2)
t       = 1:N
y10     = (sin.(t .* f1) .+ sin.(t .* (f2 .+ 0.4))) .* g.(t, N)
y20     = (sin.(t .* (f1 .+ 0.5)) .+ sin.(t .* (f2 .+ 0.5))) .* g.(t, N)
q       = Float32.(y10 .+ 0.2 .* randn.())
y       = Float32.([0y20; 0y20; y20; 0y20; 0y20] .+ 1 .* randn.())
y     .+= v1(rand(AlphaSubGaussian(n=length(y))));
# y = repeat(y, 40) #src

Q = melspectrogram(q, nfft, window = hanning, fmin = 400, fs = 48_000)
Y = melspectrogram(y, nfft, window = hanning, fmin = 400, fs = 48_000)
plot(plot(Q, title = "Query"), plot(Y, title = "Data"), link = :both, layout = (2, 1))
# We now see if we can detect the pattern using DTW with the standard squared Euclidean distance
rad  = 10 # This is the maximum allowed warping radius
dist = SqEuclidean()
w    = DTWWorkspace(sqrt.(Q.power), dist, rad, Nothing)
res  = dtwnn(w, sqrt.(Y.power))
plot(Y);
vline!([Y.time[res.loc]], l = (4, :blue), primary=false, yscale=:identity)
# That did probably not go well at all! The line indicates where the smallest distance to the pattern was, i.e., where the nearest neighbor search thinks the *onset* of the pattern is.
#---
# Let's do the same with a transport-based distance
n, m = size(Q.power)
dist = DiscreteGridTransportDistance(Cityblock(), Float32, n, n)
w    = DTWWorkspace(sqrt.(Q.power), dist, rad, Nothing)
res  = dtwnn(w, sqrt.(Y.power))
plot(Y);
vline!([Y.time[res.loc]], l = (4, :blue), primary=false)
# The line should now be a better indication of the onset of the pattern
#---
# Next, we plot the cost function for each time shift to see how they behave
function naive(a, b, dist = SqEuclidean(), r = 7)
    dists = map(1:lastlength(b)-lastlength(a)) do i
        dtw_cost(a, b[!, i:i+lastlength(a)-1], dist, r)
    end
end

plot()
for dist in [DiscreteGridTransportDistance(Cityblock(), Float32, n, n), SqEuclidean()]
    @time res = naive(Q.power, Y.power, dist, rad)
    res .-= minimum(res)
    res ./= maximum(res)
    plot!(res, lab = string(typeof(dist).name))
end
fig1 = vline!([2 * size(Y.power, 2) / 5], l = (4, :blue), primary=false)
# The Euclidean distance is rather terrible for this task. The transport based distance seems to be doing the right thing.
#---
# unfortunately, calculating this using the transport-based distance takes quite a while
@btime dtwnn($w, $(sqrt.(Y.power)), prune_endpoints=false)
# while that might not seem like much, the data we're searching through only corresponds to this many seconds, if the sample rate is 48kHz
length(y) / 48000
#---
# We are thus intersted in distances that are faster to compute. We'll explore the OTRD distance, which performs optimal transport between the roots of a linear-system representation of the signal.
fm  = TimeWindow(inner = LS(na = 8, λ=1e-3), n = nfft, noverlap = nfft÷2)
Qm  = fm(q) |> change_precision(Float64)
Ym  = fm(y) |> change_precision(Float64);
# Let's see how the world looks like through the lens of autoregressive models
plot(Qm, LPVSpectral.mel_to_hz(Q.mels)./48000, rad=false)
#---
dist  = OptimalTransportRootDistance(p = 1, β = 0.5)
distf = (x, y) -> evaluate(dist, x, y, tol = 1e-3)
res   = naive(Qm.models, Ym.models, distf, rad)
res .-= minimum(res)
res ./= maximum(res)
plot!(fig1, res, lab = "OTRD")
# This should do an okay job, but is it fast?
#---
w = DTWWorkspace(Qm.models, distf, rad)
@btime dtwnn($w, $(Ym.models), prune_endpoints=false)
# About the same, but it could be much faster if the allocations were taken care of.
# Let's speed up even further, using the root distance (RD)
Qe     = Float32.(reduce(hcat, embedding.(Qm.models))) # This extracts the roots of the linear system into a real vector
Ye     = Float32.(reduce(hcat, embedding.(Ym.models)))
dist   = SqEuclidean();
# dist   = EuclideanRootDistance(p=1, weight=e->sqrt.(residueweight(e))) #src
res    = naive(Qe, Ye, dist, rad)
res  .-= minimum(res)
res  ./= maximum(res)
plot!(res, lab = "RD")
# This should do an okay job, but is it fast?
#---
w = DTWWorkspace(Qe, dist, rad)
dtwnn(w, (Ye), prune_endpoints=false)
@btime dtwnn($w, $(Ye), prune_endpoints=false)
# That should be fast very, but there are still some room for improvements through early-termination heuristics.
GC.gc(true); GC.gc(true) #src

# ## Notes on futher performance improvements
# Searching for the nearest neighbor using DTW tends to become relatively cheaper the longer the data searched through is. This is due to pruning heauristics that abort the search early based on the smallest distance found so far.
# The main thing to watch out for is the length (in time) of the query sequence. Using a larger `nfft` effectively reduces the time resolution which greatly reduces the time-complexity, 𝒪(t²) if pruning heuristics are not effective.
# Let's see how long time it takes to search through a much longer time series (note that estimating the models also takes quite a while, but this is 𝒪(t) expensive and can be easily parallelized so I don't count that for now)
yl  = Float32.(repeat(y, 400) .+ 0.2 .* randn.())
println("Length of singal: ", length(yl)/48_000/60, " minutes, at 48kHz")
Yml = fm(yl) |> change_precision(Float32);
Yel = reduce(hcat, embedding.(Yml.models))
sleep(2) #src
@btime dtwnn($w, $(Yel), prune_endpoints=true)

# ### Future work
# The mass normalization for spectrograms has to be improved. Currently, each spectrum (per time point) is normalized to sum to one. A perhaps better strategy is to normalize the entire spectrogram and then use unbalanced mass transport per time point.

GC.gc(true); GC.gc(true) #src

#---
# # Transport in time
# We can also consider using optimal transport along the time axis instead of DTW. Let's do the same stuff as above for the TimeDistance from SpectralDistances.jl
dist = TimeDistance(inner=OptimalTransportRootDistance(p=1, β=0.5),tp=1,c=0.1)
res   = SpectralDistances.distance_profile(dist, Qm, Ym)
res .-= minimum(res)
res ./= maximum(res)
plot!(fig1, res, lab = "Time transport")
#---
function naive_time(a, b, dist)
    n = length(a)
    dists = map(1:length(b.models)-length(a.models)) do i
        Y = TimeVaryingAR(b.models[i:i+length(a.models)-1])
        dist(a, Y, tol=1e-3, check_interval=2)
    end
end
@btime naive_time($Qm, $Ym, $dist);
@btime SpectralDistances.distance_profile($dist, $(change_precision(Float32,Qm)), $(change_precision(Float32,Ym)), check_interval=5)
# It appears to be surprisingly competitive, both in terms of accuracy and performance, even though it's evaluated using a rather naive method. Since the solver used to find the `TimeDistance` internally solves the *dual* problem, a lower bound of the objective function is always available for free during the optimization. This lower bound could be used to terminate the optimization early if it rises above the smallest distance found so far.

@test mean(abs2, SpectralDistances.distance_profile(dist, Qm, Ym) - naive_time(Qm, Ym, dist)) < 1e-4


text:
# ====================================================================
#                          TYPE DEFINITIONS
# ====================================================================

# ____________________________________________________________________
# Abstract types
#
abstract type AbstractDomain end
abstract type AbstractLinearDomain    <: AbstractDomain end
abstract type AbstractCartesianDomain <: AbstractDomain end

abstract type AbstractData end
abstract type AbstractMeasures <: AbstractData end
abstract type AbstractCounts   <: AbstractData end

abstract type AbstractComponent end
abstract type AbstractComponentData end

abstract type AbstractMinimizer end

# ____________________________________________________________________
# Macro to define domain, data and associated methods for any
# dimensionality
#
macro code_ndim(ndim::Int)
    @assert ndim >= 1 "Number of dimensions must be >= 1"
    out = Expr(:block)

    # Structure def.
    if ndim > 1
        name = Symbol(:CartesianDomain_, ndim, :D)
        push!(out.args, :(
            struct $name <: AbstractCartesianDomain
                axis::Vector{Vector{Float64}}
                size::NTuple{$(ndim), Int}
                index::Vector{Int}
            end;
            Base.ndims(dom::$name) = $ndim;
            Base.size(dom::$name) = dom.size;
            Base.size(dom::$name, dim::Int) = dom.size[dim];
            Base.getindex(dom::$name, dim::Int) = dom.axis[dim];
        ))

        # Constructors
        a = ["d$(i)::AbstractVector{Float64}" for i in 1:ndim]
        b = ["deepcopy(d$(i))" for i in 1:ndim]
        c = ["length(d$(i))" for i in 1:ndim]
        s = "CartesianDomain(" * join(a, ", ") * "; index=1:" * join(c, " * ") *
            ") = CartesianDomain_$(ndim)D([" * join(b, ", ") * "], (" * join(c, ", ") * ",), index)"
        push!(out.args, Meta.parse(s))

        a = ["n$(i)::Int" for i in 1:ndim]
        b = ["one(Float64):one(Float64):n$(i)" for i in 1:ndim]
        c = ["n$(i)" for i in 1:ndim]
        s = "CartesianDomain(" * join(a, ", ") * "; index=1:" * join(c, " * ") *
            ") = CartesianDomain(" * join(b, ", ") * ", index=index)"
        push!(out.args, Meta.parse(s))
    end

    # Structure def.
    name = Symbol(:Domain_, ndim, :D)
    tmp = (ndim > 1  ?  2  :  1)
    push!(out.args, :(
        struct $name <: AbstractLinearDomain
            axis::Array{Float64, $tmp}
            vmin::Vector{Float64}
            vmax::Vector{Float64}
            length::Int
        end;
        Base.ndims(dom::$name) = $ndim;
    ))

    # getindex
    if ndim == 1
        push!(out.args, :(Base.getindex(dom::$name, dim::Int) = dom.axis))
    else
        push!(out.args, :(Base.getindex(dom::$name, dim::Int) = dom.axis[dim,:]))
    end

    # Constructors
    a = ["d$(i)::AbstractVector{Float64}" for i in 1:ndim]
    b = ["deepcopy(d$(i))" for i in 1:ndim]
    c = ["d$(i)" for i in 1:ndim]
    s = "function Domain(" * join(a, ", ") * ")\n"
    if ndim > 1
        s *= "  @assert " * join("length(" .* c .* ")", " == ") * " \"Arrays must have same length\" \n"
    end
    s *= "  return Domain_$(ndim)D(" *
        (ndim == 1  ?  "d1"  :  "[" * join(b, " ") * "]'") * ", " *
        "[" * join("minimum(" .* c .* ")", ", ") * "], " *
        "[" * join("maximum(" .* c .* ")", ", ") * "], " *
        "length(d1))\n"
    s *= "end"
    push!(out.args, Meta.parse(s))

    a = ["n$(i)::Int" for i in 1:ndim]
    b = ["one(Float64):one(Float64):n$(i)" for i in 1:ndim]
    s = "Domain(" * join(a, ", ") * ") = Domain(" * join(b, ", ") * ")"
    push!(out.args, Meta.parse(s))

    # flatten
    if ndim > 1
        s = "function flatten(dom::CartesianDomain_$(ndim)D)::Domain_$(ndim)D\n"
        s *= "  out = Matrix{Float64}(undef, $ndim, length(dom))\n"
        s *= "  for i in 1:length(dom)\n"
        a = ["d$(i)" for i in 1:ndim]
        s *= "  (" * join(a, ", ") * ") = Tuple(CartesianIndices(size(dom))[dom.index[i]])\n"
        for i in 1:ndim
            s *= "  out[$i, i] = dom.axis[$i][d$(i)]\n"
        end
        s *= "  end\n"
        a = ["out[$(i), :]" for i in 1:ndim]
        s *= "  return Domain(" * join(a, ", ") * ")\n"
        s *= "end\n"
        push!(out.args, Meta.parse(s))
    end

    name = Symbol(:Measures_, ndim, :D)
    push!(out.args, :(
        struct $name <: AbstractMeasures
            val::Array{Float64, $(ndim)}
            unc::Array{Float64, $(ndim)}
        end;
        Base.ndims(dom::$name) = $ndim;
    ))

    push!(out.args, :(
        function Measures(measure::Array{Float64, $ndim}, uncert::Array{Float64, $(ndim)})
            @assert length(measure) == length(uncert) "Measure and uncertainty arrays must have same size"
            $(name)(measure, uncert)
        end;
    ))

    push!(out.args, :(
        function Measures(measure::Array{Float64, $ndim}, uncert::Float64=one(Float64))
            u = similar(measure)
            fill!(u, uncert)
            $(name)(measure, u)
        end;
    ))

    name = Symbol(:Counts_, ndim, :D)
    push!(out.args, :(
        struct $name <: AbstractCounts
            val::Array{Int, $(ndim)}
        end;
        Base.ndims(dom::$name) = $ndim;
    ))

    push!(out.args, :(
        function Counts(measure::Array{Int, $ndim})
            $(name)(measure)
        end
    ))

    s = Vector{String}()
    push!(s, "function evaluate!(cdata::FuncWrap_cdata, output::Vector{Float64}, domain::Domain_$(ndim)D, params...)")
    push!(s, "  output .= cdata.func(" * join("domain[" .* string.(collect(1:ndim)) .* "]", ", ") * ", params...)")
    push!(s, "end")
    push!(out.args, Meta.parse(join(s, "\n")))

    return esc(out)
end

# The following is here, rather than in components.jl, to avoid a
# compilation error: `FuncWrap_cdata` is used in `@code_ndim`
struct FuncWrap_cdata <: AbstractComponentData
    func::Function
end

@code_ndim 1
@code_ndim 2
@code_ndim 3

getaxismin(dom::AbstractLinearDomain, dim::Int) = dom.vmin[dim]
getaxismax(dom::AbstractLinearDomain, dim::Int) = dom.vmax[dim]
getaxisextrema(dom::AbstractLinearDomain, dim::Int) = (dom.vmin[dim], dom.vmax[dim])
length(dom::AbstractLinearDomain) = dom.length
length(dom::AbstractCartesianDomain) = length(dom.index)
length(data::AbstractData) = length(data.val)
size(data::AbstractData) = size(data.val)
size(data::AbstractData, dim::Int) = size(data.val)[dim]

# Methods to "flatten" a multidimensional object <: AbstractData into a 1D one
flatten(dom::AbstractLinearDomain) = dom
flatten(data::Measures_1D, dom::AbstractLinearDomain)::Measures_1D = data
flatten(data::Counts_1D, dom::AbstractLinearDomain)::Counts_1D = data
flatten(data::AbstractMeasures, dom::AbstractLinearDomain)::Measures_1D    = Measures_1D(data.val[:], data.unc[:])
flatten(data::AbstractMeasures, dom::AbstractCartesianDomain)::Measures_1D = Measures_1D(data.val[dom.index], data.unc[dom.index])
flatten(data::AbstractCounts  , dom::AbstractLinearDomain)::Counts_1D      = Counts_1D(data.val)
flatten(data::AbstractCounts  , dom::AbstractCartesianDomain)::Counts_1D   = Counts_1D(data.val[dom.index])

# function append!(dest::T, source::T) where T <: AbstractMeasures
#     append!(dest.val, source.val)
#     append!(dest.unc, source.unc)
#     return dest
# end

# """
# # reshape
#
# Reshape an array according to the size of a CartesianDomain object
# """
# function reshape(array::AbstractArray, dom::AbstractCartesianDomain)
#     @assert length(array) == length(dom) "Domain and array must have the same length"
#     out = Array{Float64}(undef, size(dom)...)
#     out .= NaN
#     out[dom.index] .= array
#     return out
# end

# ____________________________________________________________________
# Instrument
#
mutable struct CompEvaluation
    fixed::Bool
    npar::Int
    log::Vector{Bool}
    counter::Int
    cdata::AbstractComponentData
    lastParams::Vector{Float64}
    result::Vector{Float64}
    value::Vector{Float64}
end

mutable struct Instrument
    compile::Bool
    ldomain::AbstractLinearDomain
    domain::AbstractDomain

    code::String
    funct::Function
    counter::Int

    compnames::Vector{Symbol}  # only involved components
    compevals::Vector{CompEvaluation}

    exprnames::Vector{Symbol}
    exprs::Vector{Expr}
    exprcmp::Vector{Bool}
    exprevals::Vector{Vector{Float64}}
end

# ____________________________________________________________________
# Wrapper for components
#
mutable struct WComponent
    cname::Symbol
    comp::AbstractComponent
    fixed::Bool
end

# ____________________________________________________________________
# Main Model structure
#
mutable struct Model
    comp::OrderedDict{Symbol, WComponent}
    instruments::Vector{Instrument}
    index1d::Vector{Int}
    buffer1d::Vector{Float64}
    Model(::Nothing) = new(OrderedDict{Symbol, WComponent}(),
                           Vector{Instrument}(),
                           Vector{Int}(), Vector{Float64}())
end

# ____________________________________________________________________
# ParameterPrivate: it allows to distinguish a scalar parameter
# (`index` = 0) from a vector of parameters and to access parent Model
mutable struct ParameterPrivate
    model::Union{Nothing,Model}
    cname::Symbol
    pname::Symbol
    index::Int
    fitval::Float64
    fitunc::Float64
    ParameterPrivate() = new(nothing, :-, :-, 0, NaN, NaN)
end

# Parameter
mutable struct Parameter
    _private::ParameterPrivate
    val::Float64
    low::Float64              # lower limit value
    high::Float64             # upper limit value
    step::Float64
    fixed::Bool               # true = fixed; false = free
    cfixed::Bool
    log::Bool
    expr::String
    Parameter(value::Number) = new(ParameterPrivate(), float(value), -Inf, +Inf, NaN, false, false, false, "")
end

# ____________________________________________________________________
# Fit results
#
struct FitParam
    val::Float64
    unc::Float64
    fixed::Bool
end

struct FitComp
    params::OrderedDict{Symbol, Union{FitParam, Vector{FitParam}}}
end

struct FitResult
    bestfit::OrderedDict{Symbol, FitComp}
    ndata::Int
    dof::Int
    cost::Float64
    status::Symbol      #:Optimal, :NonOptimal, :Warn, :Error
    log10testprob::Float64
    elapsed::Float64
end


text:
import JLD2
import FileIO
import Geodesy: LLA, euclidean_distance, ENUfromLLA, wgs84, LLAfromENU, ENU
import Printf: @printf
import Dates: Dates, unix2datetime, DateTime
import TimeZones: localzone, ZonedDateTime, astimezone, @tz_str
import DataFrames: DataFrame, leftjoin
import OrderedCollections: OrderedDict

export cached, AOI, distance, minute, minutes, second, seconds, dts, @dts_str, xleftjoin, closest
export LLA, ZonedDateTime, geo2pos, pos2geo

### caching with JLD2 files

function cached(f, filename)
  if isfile(filename)
    return JLD2.jldopen(filename) do f
      JLD2.read(f, "df")
    end
  else
    df = f()
    JLD2.@save filename df
    return df
  end
end

### Geodesy utils

struct AOI{T}
  topleft::LLA{Float64}
  bottomright::LLA{Float64}
  mapimg::T
end

# set environment variable GEOAPIFY_APIKEY to auto-download maps (see https://www.geoapify.com)
# see https://apidocs.geoapify.com/docs/maps/map-tiles/ for list of available styles

function AOI(topleft::LLA{Float64}, bottomright::LLA{Float64}; style="osm-bright", width=1024, cachedir=nothing)
  mapimg = nothing
  lat1, lon1, lat2, lon2 = topleft.lat, topleft.lon, bottomright.lat, bottomright.lon
  height = round(Int, width * abs(lat2 - lat1) / abs(lon2 - lon1))
  try
    cachedir == nothing && (cachedir = joinpath(tempdir(), "arltoolkit", "cache"))
    filename = joinpath(cachedir, "map-$style-$lon1-$lat1-$lon2-$lat2-$width-$height.png")
    if isfile(filename)
      mapimg = FileIO.load(filename)
    else
      if "GEOAPIFY_APIKEY" ∈ keys(ENV)
        apikey = ENV["GEOAPIFY_APIKEY"]
        mkpath(cachedir)
        download("https://maps.geoapify.com/v1/staticmap?style=$style&format=png&" *
          "area=rect:$lon1,$lat1,$lon2,$lat2&width=$width&height=$height&apiKey=$apikey",
          filename)
        mapimg = FileIO.load(filename)
      end
    end
  catch ex
    @warn "Could not load map: $ex"
  end
  AOI(topleft, bottomright, mapimg)
end

function Base.in(x::LLA, aoi::AOI)
  (x.lat < aoi.bottomright.lat || x.lat > aoi.topleft.lat) && return false
  (x.lon < aoi.topleft.lon || x.lon > aoi.bottomright.lon) && return false
  true
end

distance(p1::LLA, p2::LLA) = euclidean_distance(p1, p2)

function prettyprint(io::IO, x::LLA)
  @printf(io, "%0.6f°%c %0.6f°%c", abs(x.lat), x.lat < 0 ? 'S' : 'N', abs(x.lon), x.lon < 0 ? 'W' : 'E')
  iszero(x.alt) || @printf(io, " %+0.3fm", x.alt)
end

"""
    geo2pos(geo::LLA, origin::LLA)

Convert latitude/longitude to local coordinates. If altitude information
is non-zero, the returned position vector is a 3-vector, otherwise a 2-vector.
"""
function geo2pos(geo::LLA, origin::LLA)
  pos = ENUfromLLA(origin, wgs84)(geo)[1:3]
  geo.alt == 0.0 && (pos = pos[1:2])
  pos
end

"""
    pos2geo(pos::AbstractArray{<:Real}, origin::LLA)

Convert local coordinates to latitude/longitude. The input may be a 2-vector
or a 3-vector. If the input is a 2-vector, the altitude is assumed to be zero.
"""
function pos2geo(pos::AbstractArray{<:Real}, origin::LLA)
  length(pos) == 2 && (pos = vcat(pos, zero(eltype(pos))))
  length(pos) == 3 || throw(ArgumentError("Bad position vector"))
  LLAfromENU(origin, wgs84)(ENU(pos))
end

### datetime conversions

const minute = Dates.Minute(1)
const second = Dates.Second(1)
const minutes = minute
const seconds = second

function dts(yyyy, mm, dd, HH, MM, SS; tz=localzone())
  t = ZonedDateTime(yyyy, mm, dd, HH, MM, SS, tz"UTC")
  astimezone(t, tz)
end

function dts(s::String; tz=localzone(), format="yyyy-mm-dd HH:MM:SS")
  ZonedDateTime(DateTime(s, format), tz)
end

function dts(t; tz=localzone())
  t = ZonedDateTime(unix2datetime(t), tz"UTC")
  astimezone(t, tz)
end

macro dts_str(s)
  dts(s)
end

Base.isapprox(t1::ZonedDateTime, t2::ZonedDateTime; atol=10seconds) = abs(t2 - t1) ≤ atol

function prettyprint(io::IO, x::ZonedDateTime)
  zone = x.zone.offset.std.value == 28800 ? "SGT" : string(x.zone)
  print(io, Dates.format(x, "yyyy-mm-dd HH:MM:SS"), " ", zone)
end

### dataframe join utils

function closest(a, b; atol=Inf64)
  Δab = abs.(a .- b)
  ndx = argmin(Δab)
  Δab[ndx] ≤ atol && return ndx
  nothing
end

closest(; atol) = (a, b) -> closest(a, b; atol)

function xleftjoin(df1, df2, on, cols; compare=isequal, find=(a, b)->findfirst((compare).(a, b)))
  df3 = DataFrame()
  for row1 ∈ eachrow(df1)
    ndx = find(row1[on], df2[!,on])
    if ndx !== nothing
      data = OrderedDict{Symbol,Any}(pairs(df2[ndx,cols]))
      data[:__key] = row1[on]
      push!(df3, data; cols=:union)
    end
  end
  leftjoin(df1, df3; on = on=>:__key)
end


text:
export arrangeRemoteCallDataIntoLocalData,getSourcesIndicesOfKthWorker

function arrangeRemoteCallDataIntoLocalData(DremoteCall::Union{Array{Array{Float64,2},1},Array{Array{Complex128,2},1}})
N = length(DremoteCall);
l = 0;
for k=1:N
	l+=size(DremoteCall[k],2);
end
Dlocal = zeros(eltype(DremoteCall[1]),size(DremoteCall[1],1),l)
for k=1:N
	Dlocal[:,k:N:end] = DremoteCall[k];
end
return Dlocal;
end

function getSourcesIndicesOfKthWorker(nworkers::Int64,k::Int64,nsrc::Int64)
return k:nworkers:nsrc;
end


text:
module Spells

using AutoHashEquals
using Dictionaries
using MLStyle
using SPIRV
using SPIRV: ID
using Graphs

const Optional{T} = Union{Nothing, T}
const magic_number = 0x12349876

include("utils.jl")
include("types.jl")
include("emit/types.jl")
include("emit/cfg_constructs.jl")
include("emit/functions.jl")

export
        Context,
        Func,
        PrimitiveType,
        CompositeType

end


text:
using Test, LinearAlgebra
using ThreeDeconv

function adjointnessCheck(A::ThreeDeconv.LinearOperator)
    x = randn(eltype(A), A.input_shape)
    forward_output = A.forward(x)

    y = randn(eltype(A), A.output_shape)
    adjoint_output = A.adjoint(y)

    return dot(forward_output, y) ≈ dot(x, adjoint_output)
end

function adjointnessCheck(A::ThreeDeconv.StackedLinearOperator)
    x = randn(A.input_shape)
    forward_output = vcat((v[:] for v in A.forward(x))...)

    y = [randn(output_shape) for output_shape in A.output_shape]
    adjoint_output = A.adjoint(y)

    return dot(forward_output, vcat(y...)) ≈ dot(x, adjoint_output)
end

@testset "linearoperator" begin
    A = randn(12, 10)
    LinOpA = ThreeDeconv.LinearOperator(
        (10, 1),
        (12, 1),
        10,
        12,
        x -> A * x,
        y -> A' * y,
        false,
        Float64,
    )
    @test adjointnessCheck(LinOpA)
    @test opnorm(A) ≈ ThreeDeconv.opnorm(LinOpA)

    B = randn(12, 10)
    LinOpB = ThreeDeconv.LinearOperator(
        (10, 1),
        (12, 1),
        10,
        12,
        x -> B * x,
        y -> B' * y,
        false,
        Float64,
    )

    LinOpK = ThreeDeconv.StackedLinearOperator([LinOpA, LinOpB])
    K = vcat(A, B)
    @test adjointnessCheck(LinOpK)
    @test opnorm(A) ≈ ThreeDeconv.opnorm(LinOpA)
    @test opnorm(B) ≈ ThreeDeconv.opnorm(LinOpB)
    @test opnorm(K) ≈ ThreeDeconv.opnorm(LinOpK)

end


text:
# =========================================================================================
"""
    struct LocalBEMResult{T, E} <: BEMResult{T, E}
        model::Model{T, E}
        u    ::Vector{T}   # [γ₀int(φ*)](ξ)    ∀ ξ ∈ Ξ; premultiplied by 4π⋅ε0
        q    ::Vector{T}   # [γ₁int(φ*)](ξ)    ∀ ξ ∈ Ξ; premultiplied by 4π⋅ε0
        umol ::Vector{T}   # [γ₀int(φ*mol)](ξ) ∀ ξ ∈ Ξ; premultiplied by 4π⋅ε0
        qmol ::Vector{T}   # [γ₁int(φ*mol)](ξ) ∀ ξ ∈ Ξ; premultiplied by 4π⋅ε0
    end

Result data of the local solving process to be used for potential computation and
post-processing, with `Ξ` being the list of observation points, that is, the set of
triangle centroids.
"""
struct LocalBEMResult{T, E} <: BEMResult{T, E}
    """Surface model"""
    model::Model{T, E}
    """[γ₀int(φ*)](ξ) for all observation points ξ"""
    u::Vector{T}
    """[γ₁int(φ*)](ξ) for all observation points ξ"""
    q::Vector{T}
    """[γ₀int(φ*mol)](ξ) for all observation points ξ"""
    umol::Vector{T}
    """[γ₁int(φ*mol)](ξ) for all observation points ξ"""
    qmol::Vector{T}
end

# =========================================================================================
"""
    solve{T, L <: LocalityType}(
             ::L,
        model::Model{T, Triangle{T}}
    )

Computes the full local or nonlocal cauchy data on the surface of the biomolecule using an *explicit* representation of the BEM system.

!!! note
    The explicit system representations require a substantial amount of memory and might easily
    exceed the available resources on your machine (the memory demands scale quadratically with
    the number of surface elements). In you encounter this problem, please consider switching to
    [`solve_implicit`](@ref) or our CUDA-accelerated solvers from the
    [`CuNESSie.jl`](https://github.com/tkemmer/CuNESSie.jl) package instead, which provide the
    same function interface and are freely interchangeable.

# Return type
`LocalBEMResult{T, Triangle{T}}` or `NonlocalBEMResult{T, Triangle{T}}`
"""
function solve(
             ::Type{LocalES},
        model::Model{T, Triangle{T}}
    ) where T
    # observation points ξ
    Ξ = [e.center for e in model.elements]

    # compute molecular potentials for the point charges;
    # molecular potentials are initially premultiplied by 4π⋅ε0⋅εΩ
    umol = model.params.εΩ \   φmol(model)
    qmol = model.params.εΩ \ ∂ₙφmol(model)

    # convenience aliases
    εΩ = model.params.εΩ
    εΣ = model.params.εΣ
    numelem = length(model.elements)

    # system matrix M
    m = zeros(T, numelem, numelem)
    k = Array{T}(undef, numelem, numelem)
    v = zeros(T, numelem, numelem)

    # right-hand sides bᵤ and b_q
    b = zeros(T, numelem)

    # Mᵤ = (1 + εΩ/εΣ) ⋅ σ;
    # since all other components of the system matrix will be premultiplied by 4π,
    # do the same for σ here
    pluseye!(m, (1 + εΩ/εΣ) * T(4π * σ))

    #=
        generate and apply V
    =#
    # M_q = V
    Rjasanow.laplacecoll!(SingleLayer, v, model.elements, Ξ)

    # bᵤ -= εΩ/εΣ ⋅ V ⋅ qmol
    gemv!(-εΩ/εΣ, v, qmol, b)

    #=
        generate and apply K
    =#
    Rjasanow.laplacecoll!(DoubleLayer, k, model.elements, Ξ)

    # Mᵤ += (εΩ/εΣ - 1) ⋅ K
    axpy!(εΩ/εΣ - 1, k, m)

    # bᵤ = (K - σ) ⋅ umol
    # again, we apply a prefactor of 4π to σ to match the other components of the vector
    pluseye!(k, -T(4π * σ))
    gemv!(one(T), k, umol, b)

    #=
        u = b / M, with
        b = (K - σ) ⋅ umol - εΩ/εΣ ⋅ V ⋅ qmol, and
        M = (1 + εΩ/εΣ) ⋅ σ + (εΩ/εΣ - 1) ⋅ K.
    =#
    u = m \ b

    # b_q = (σ + K) ⋅ u
    fill!(b, zero(T))
    pluseye!(k, T(8π * σ)) # meh...
    gemv!(one(T), k, u, b)

    #=
        q = V^{-1}⋅(σ + K)u
    =#
    q = v \ b

    LocalBEMResult(model, u, q, umol, qmol)
end

# =========================================================================================
"""
    struct LocalSystemMatrix{T} <: AbstractArray{T, 2}
        K     ::InteractionMatrix{T, Vector{T}, Triangle{T}, Kfun{T}}   # Potential matrix `K`
        params::Option{T}                                               # System constants
    end

Implicit representation of the first local BEM system matrix.
"""
struct LocalSystemMatrix{T} <: AbstractArray{T, 2}
    """Potential matrix `K`"""
    K     ::InteractionMatrix{T, Vector{T}, Triangle{T}, Kfun{T}}
    """System constants"""
    params::Option{T}
end

Base.size(A::LocalSystemMatrix{T}) where T = size(A.K)

function LinearAlgebra.diag(
    A::LocalSystemMatrix{T},
    k::Int=0
) where T
    k != 0 && error("diag not defined for k != 0 on ", typeof(A))
    (T(2π) * (1 + A.params.εΩ / A.params.εΣ)) .* ones(T, size(A, 1))
end

function LinearAlgebra.mul!(
    dst::AbstractArray{T, 1},
    A::LocalSystemMatrix{T},
    x::AbstractArray{T, 1}
) where T
    dst .= A * x
end

function Base.:*(
    A::LocalSystemMatrix{T},
    x::AbstractArray{T, 1}
) where T
    frac = A.params.εΩ / A.params.εΣ
    (T(2π) * (1 + frac)) .* x .+ ((frac - 1) .* (A.K * x))
end

# =========================================================================================
"""
    solve_implicit{T, L <: LocalityType}(
             ::L,
        model::Model{T, Triangle{T}}
    )

Computes the full local or nonlocal cauchy data on the surface of the biomolecule using an *implicit* representation of the BEM system.

!!! note
    The implicit system representations provide a small memory footprint even for large
    biomolecular systems (the memory demands scale only linearly with the number of surface
    elements). However, this comes at the cost of additional computation time, which might exceed
    several hours, depending on the system size. In you encounter this problem, please consider
    switching to [`solve`](@ref) or our CUDA-accelerated solvers from the
    [`CuNESSie.jl`](https://github.com/tkemmer/CuNESSie.jl) package instead, which provide the
    same function interface and are freely interchangeable.

# Return type
`LocalBEMResult{T, Triangle{T}}` or `NonlocalBEMResult{T, Triangle{T}}`
"""
function solve_implicit(
             ::Type{LocalES},
        model::Model{T, Triangle{T}}
    ) where T
    # observation points ξ
    Ξ = [e.center for e in model.elements]

    # compute molecular potentials for the point charges;
    # molecular potentials are initially premultiplied by 4π⋅ε0⋅εΩ
    umol = model.params.εΩ \   φmol(model, tolerance=_etol(T))
    qmol = model.params.εΩ \ ∂ₙφmol(model)

    # potential matrices
    K = InteractionMatrix(Ξ, model.elements, Kfun{T}())
    V = InteractionMatrix(Ξ, model.elements, Vfun{T}())

    # first system
    b = K * umol .- (T(2π) .* umol) .- (model.params.εΩ/model.params.εΣ .* (V * qmol))
    A = LocalSystemMatrix(K, model.params)
    u = _solve_linear_system(A, b)

    # second system
    b .= T(2π) .* u .+ (K * u)
    q = _solve_linear_system(V, b)

    LocalBEMResult(model, u, q, umol, qmol)
end


text:
# This is tools for testing ChainRulesTestUtils itself
# if they were less nasty in implementation we might consider moving them to a package
# MetaTesting.jl

# need to bring this into scope explictly so can use in @testset nonpassing_results
using Test: DefaultTestSet

"""
    nonpassing_results(f)

`f` should be a function that takes no argument, and calls some code that used `@test`.
Invoking it via `nonpassing_results(f)` will prevent those `@test` being added to the
current testset, and will return a collection of all nonpassing test results.
"""
function nonpassing_results(f)
    mute() do
        nonpasses = []
        # Specify testset type incase parent testset is some other typer
        @testset DefaultTestSet "nonpassing internal" begin
            f()
            ts = Test.get_testset()  # this is the current testset "nonpassing internal"
            nonpasses = _extract_nonpasses(ts)
            # Prevent the failure being recorded in parent testset.
            empty!(ts.results)
            ts.anynonpass = false
        end
        # Note: we allow the "nonpassing internal" testset to still be pushed as an empty
        # passing testset in its parent testset. We could remove that if we wanted
        return nonpasses
    end
end

"""
    mute(f)

Calls `f()` silencing stdout.
"""
function mute(f)
    # TODO: once we are on Julia 1.6 this can be change to just use
    # `redirect_stdout(devnull)` See: https://github.com/JuliaLang/julia/pull/36146
    mktemp() do path, tempio
        redirect_stdout(tempio) do
            f()
        end
    end
end

"extracts as flat collection of failures from a (potential nested) testset"
_extract_nonpasses(x::Test.Result) = [x,]
_extract_nonpasses(x::Test.Pass) = Test.Result[]
_extract_nonpasses(ts::Test.DefaultTestSet) = _extract_nonpasses(ts.results)
function _extract_nonpasses(xs::Vector)
    if isempty(xs)
        return Test.Result[]
    else
        return mapreduce(_extract_nonpasses, vcat, xs)
    end
end

"""
    fails(f)

`f` should be a function that takes no argument, and calls some code that used `@test`.
`fails(f)` returns true if at least 1 `@test` fails.
If a test errors then it will display that error and throw an error of its own.
"""
function fails(f)
    results = nonpassing_results(f)
    did_fail = false
    for result in results
        did_fail |= result isa Test.Fail
        if result isa Test.Error
            # Log a error message, with original backtrace
            show(result)
            # Sadly we can't throw the original exception as it is only stored as a String
            error("Error occurred during `fails`")
        end
    end
    return did_fail
end

#Meta Meta tests
@testset "meta_testing_tools.jl" begin
    @testset "Checking for non-passes" begin
        @testset "No Tests" begin
            fails = nonpassing_results(()->nothing)
            @test length(fails) === 0
        end

        @testset "No Failures" begin
            fails = nonpassing_results(()->@test true)
            @test length(fails) === 0
        end

        @testset "Single Test" begin
            fails = nonpassing_results(()->@test false)
            @test length(fails) === 1
            @test fails[1].orig_expr == false
        end

        @testset "Single Testset" begin
            fails = nonpassing_results() do
                @testset "inner" begin
                    @test false == true
                    @test true == false
                end
            end
            @test length(fails) === 2
            @test fails[1].orig_expr == :(false==true)
            @test fails[2].orig_expr == :(true==false)
        end
    end

    @testset "fails" begin
        @test !fails(()->@test true)
        @test fails(()->@test false)
        @test !fails(()->@test_broken false)

        @test fails() do
            @testset "eg" begin
                @test true
                @test false
                @test true
            end
        end

        @test_throws Exception mute() do  # mute it so we don't see the reprinted error.
            fails(()->@test error("Bad"))
        end
    end
end