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

                     @objective(m, Min, q)

                     return m


text:
include("../src/merge_segments.jl")

@testset "SegmentedImage Merge" begin
    # Set up 3 variables: test_graph, img, seg for use in tests
    g = MetaGraph(6)

    for v in 1:6
       set_prop!(g, v, :total_color, v)
       set_prop!(g, v, :pixel_count, v)
       set_prop!(g, v, :mean_color, Gray(Float64(v)))
       set_prop!(g, v, :labels, v)
    end

    for i in 1:5
       add_edge!(g, i, i+1, :weight, i + i + 1)
    end

    img = fill(1.0, (10, 10))
    img[:, 5:10] .= 2.0
    img = map(Gray, img)
    seg = fast_scanning(img, 0.5)
    test_graph = copy(g)
    # end setup

    # Test 1. Merge 2 nodes
    g = copy(test_graph)
    merge_node_props!(g, Edge(1, 2))

    # Check that labels got properly merged
    @test sort(get_prop(g, 2, :labels)) == [1, 2]

    # Check that the other node props were properly updated
    @test get_prop(g, 2, :total_color) ==  3
    @test get_prop(g, 2, :pixel_count) ==  3
    @test get_prop(g, 2, :mean_color) == 1

    # Ensure props were cleared
    @test length(props(g, 1)) == 0

    # Ensure no weights were affected.
    for v in 1:5
        @test get_prop(g, Edge(v, v+1), :weight) == (v + v + 1)
    end

    # Test 2. Merge all nodes
    g = copy(test_graph)
    for i in 1:5
        merge_node_props!(g, Edge(i, i+1))
    end

    @test sort(get_prop(g, 6, :labels)) == 1:6
    @test get_prop(g, 6, :total_color) ==  sum(1:6)
    @test get_prop(g, 6, :pixel_count) ==  sum(1:6)
    @test get_prop(g, 6, :mean_color) == 1

    for i in 1:5
        @test length(props(g, i)) == 0
    end

    # Test 3. add_neighboring_edges!
    g = copy(test_graph)
    add_edge!(g, 1, 6, :weight, 7)  # setup

    added = add_neighboring_edges!(g, Edge(1, 2))

    @test length(added) == 2
    @test has_edge(g, Edge(2, 3))
    @test has_edge(g, Edge(2, 6))

    # Test 4. seg to graph
    g = seg_to_graph(seg)
    @test length(edges(g)) == 1
    @test length(vertices(g)) == 2
    @test has_edge(g, Edge(1, 2))
    @test haskey(props(g, Edge(1, 2)), :weight)
    @test get_prop(g, Edge(1, 2), :weight) ≈ Colors.colordiff(1.0, 2.0)

    # Test 5. resegment
    g = seg_to_graph(seg)

    # identity resegment
    seg2 = resegment(seg, g)
    for f in [labels_map, segment_labels, segment_pixel_count, segment_mean]
        @test f(seg2) == f(seg)
    end

    # Test 6. merge
    seg2 = merge_segments(seg, 40)  # 40 > colordiff(1.0, 2.0)
    @test segment_labels(seg2) == [1]
end

ext: jl


text:
function free_phage_step(phage, model)
    agent = model[phage]

    nearby_cells = nearby_t(:bacterium, phage, model)
    isempty(nearby_cells) && return

    target_cell = rand(model.rng, nearby_cells)
    if rand(model.rng) < p_adsorption(target_cell, model)
        kind = agent.kind
        if kind === :temperate
            if rand(model.rng) < p_lysis(phage, model)
                agent.kind = :induced_temperate
            end
        end
        if kind === :virulent || kind === :induced_temperate
            infect(phage, target_cell, model)
        end
    end
end

function complex_step!(model)
    bacteria_death_inherent(by_single_type(:bacterium)(model), model)
    bacteria_death_lysis(by_single_type(:bacterium)(model), model)

    phages = by_single_type(:phage)(model)
    if isempty(phages)
        (model.properties.phages_count == 0) && return

        model.properties.bacteria_count = length(by_single_type(:bacterium)(model))
        model.properties.phages_count = 0
        return
    end

    filter!(id -> model[id].state === :free, phages)
    for phage ∈ phages
        free_phage_step(phage, model)
    end

    filter!(id -> model[id].state === :free, phages)
    for phage ∈ phages
        tick_phage(phage, model)
        phage_decay(phage, model)
    end

    cells = by_single_type(:bacterium)(model)
    filter!(id -> isempty(model[id].phages_inside), cells)
    p = p_grow(model)
    for cell ∈ cells
        grow(cell, p, model)
    end

    if model.properties.diffuse
        cells = by_single_type(:bacterium)(model)
        phages = by_single_type(:phage)(model)
        filter!(id -> model[id].state === :free, phages)
        if model.properties.environment === :semi_solid
            for cell ∈ cells
                diffuse(cell, model)
            end
            for phage ∈ phages
                diffuse(phage, model)
            end
        end
    end

    model.properties.bacteria_count = length(by_single_type(:bacterium)(model))
    model.properties.phages_count = length(by_single_type(:phage)(model))
end


text:
# This file allows loading (and dynamically reloading) the configs from
# the configs file: configs.jl.

config_read_timer = WallTimer()
config_print_error_timer = WallTimer()

config_reload_time_s = 0.5  # seconds

"""
    reloadConfigsFiles()

 Reload "configs.jl" and assign global variables for top-level keys.
 This should be called once for release build during compile time, and can be
 called as part of the run-loop for debug builds to allow on-the-fly editing of
 the configs file.
 """
 function reloadConfigsFiles(filenames=["configs.jl"])
     if !started(config_read_timer) || elapsed(config_read_timer) > config_reload_time_s
         for filename in filenames
             try
               include(joinpath(@__DIR__, "..", "assets/configs.jl"))
             catch e
                 print_config_error("ERROR: Failed to reload Configs File!:\n '$(e)'")
             end
         end
         start!(config_read_timer)
     end
end

function print_config_error(msg)
    if !started(config_print_error_timer) || elapsed(config_print_error_timer) > 5
        println(msg)
        start!(config_print_error_timer)
    end
end

# Run once during compile time.
reloadConfigsFiles()


text:
Base.@kwdef struct ParabolicOutsideCases <: AbstractOutsideCases
  days::Float64 = 30
  peak::Float64 = 20
  height::Float64 = 0.8
  minimum::Float64 = 0.02
  frequency::TimePoint = 0.05 |> MocosSim.TimePoint
  strain::StrainKind = ChineseStrain
  time_limit::TimePoint = typemax(TimePoint)
end

function (f::ParabolicOutsideCases)(state::AbstractSimState, ::AbstractSimParams)
  days = (f.days/2)^2
  N = length(state.individuals)
  individuals = 1:N

  for infection_time in 0.0:f.frequency:f.time_limit
    if rand(state.rng) < max((Float64(-abs(infection_time-f.peak)^2)*f.height)/days + f.height, f.minimum)
      person_id = sample(state.rng, individuals)
      event = Event(Val(OutsideInfectionEvent), infection_time, person_id, f.strain)
      push!(state.queue, event)
    end
  end

end

Base.@kwdef struct InstantOutsideCases <: AbstractOutsideCases
  num_infections::Float64 = 50
  import_time::TimePoint = 0.0 |> MocosSim.TimePoint
  strain::StrainKind = ChineseStrain
end

function (f::InstantOutsideCases)(state::AbstractSimState, ::AbstractSimParams)
  N = length(state.individuals)
  individuals = 1:N

  for _ in 1:f.num_infections
    person_id = sample(state.rng, individuals)
    event = Event(Val(OutsideInfectionEvent), f.import_time, person_id, f.strain)
    push!(state.queue, event)
  end
end

# This all to avoid using @eval and others
const imported_cases = Dict{String, Type{T} where T}(
    "ParabolicOutsideCases" => ParabolicOutsideCases,
    "InstantOutsideCases"   => InstantOutsideCases
)

make_imported_cases(name::AbstractString; args...) = imported_cases[name](;args...)


text:
#EasyPlot: Annotation objects and tools.
#-------------------------------------------------------------------------------
#=TODO
 - Use a better system for .align instead of using symbols? (ex: TOP|LEFT).
=#

#==Useful constants
===============================================================================#

#==Abstract types
===============================================================================#
#All annotation can use _apply() attributes system:
abstract type AbstractPlotAnnotation <: AbstractAttributeReceiver end

#==Main types
===============================================================================#
struct DefaultFontAnnotation <: AbstractFont; end

struct LineAttributes
	style #::Symbol #will be none if set to nothing while a glyph.shape exists
	width::PReal #[0, 10]
	color::ColorRef
end
LineAttributes(;style=nothing, width=1, color=nothing) = LineAttributes(style, width, color)

struct GlyphAttributes
#==IMPORTANT:
Edge width & color taken from LineAttributes
==#
	shape::Symbol #because "type" is reserved
	size::PReal #of glyph. Edge width taken from LineAttributes
	color::ColorRef #glyph linecolor. = nothing to match line.color
	fillcolor::ColorRef
end
GlyphAttributes(;shape=:none, size=1, color=nothing, fillcolor=nothing) =
	GlyphAttributes(shape, size, color, fillcolor)

mutable struct TextAnnotation <: AbstractPlotAnnotation
	text::String
	pos::Pos2DOffset
#	font::AbstractFont
	angle::PReal #Degrees
	align::Symbol #tl, tc, tr, cl, cc, cr, bl, bc, br
	strip::Int #0 = all
end
_TextAnnotation(text::String, x, y, angle, align, strip) =
	TextAnnotation(text, Pos2DOffset(Point2D(x,y), Vector2D(0,0), Vector2D(0,0)),
		PReal(angle), Symbol(align), Int(strip) #Enforce type to ensure WILL call main constructor.
	)
TextAnnotation(text::String, args...; x::Real=PNaN, y::Real=PNaN, angle::Real=0.0, align::Symbol=:cc, strip::Int=0, kwargs...) =
	_apply(_TextAnnotation(text, x, y, angle, align, strip), args, kwargs)

mutable struct HVMarker <: AbstractPlotAnnotation
	isvert::Bool #else: horizontal
	pos::PReal
	line::LineAttributes
	strip::Int #0 = all
end
_HVMarker(isvert, pos, strip) =
	HVMarker(isvert, pos, LineAttributes(), strip)
vmarker(pos::Real, args...; strip=0, kwargs...) =
	_apply(_HVMarker(true, pos, strip), args, kwargs)
hmarker(pos::Real, args...; strip=1, kwargs...) =
	_apply(_HVMarker(false, pos, strip), args, kwargs)

#==Register constructors with cons() interface
===============================================================================#
cons(::DS{:vmarker}, args...; kwargs...) = vmarker(args...; kwargs...)
cons(::DS{:hmarker}, args...; kwargs...) = hmarker(args...; kwargs...)
cons(::DS{:atext}, args...; kwargs...) = TextAnnotation(args...; kwargs...)

#==Implement _apply() methods
===============================================================================#
function _apply(a::TextAnnotation, ::DS{:prop}; x=nooverwrite, y=nooverwrite,
		angle=nooverwrite, align=nooverwrite, strip=nooverwrite)
	_x = a.pos.v.x; _y= a.pos.v.y;
		NoOverwrite(x) || (_x = x)
		NoOverwrite(y) || (_y = y)
		a.pos.v = Point2D(_x, _y)
	NoOverwrite(angle) || (a.angle = angle)
	NoOverwrite(align) || (a.align = align)
	NoOverwrite(strip) || (a.strip = strip)
	return a
end
function _apply(a::TextAnnotation, ::DS{:offset}; x=nooverwrite, y=nooverwrite)
	_x = a.pos.offset.x; _y= a.pos.offset.y;
		NoOverwrite(x) || (_x = x)
		NoOverwrite(y) || (_y = y)
		a.pos.offset = Vector2D(_x, _y)
	return a
end
function _apply(a::TextAnnotation, ::DS{:reloffset}; x=nooverwrite, y=nooverwrite)
	_x = a.pos.reloffset.x; _y= a.pos.reloffset.y;
		NoOverwrite(x) || (_x = x)
		NoOverwrite(y) || (_y = y)
		a.pos.reloffset = Vector2D(_x, _y)
	return a
end

function _apply(m::HVMarker, ::DS{:line}; style=nooverwrite, width=nooverwrite, color=nooverwrite)
	_style = m.line.style; _width = m.line.width; _color = m.line.color
		NoOverwrite(style) || (_style = style)
		NoOverwrite(width) || (_width = width)
		NoOverwrite(color) || (_color = _resolve_ColorRef(color))
		m.line = LineAttributes(style=_style, width=_width, color=_color)
	return m
end

#Last line


text:
function solve_timeless(model::State_Space_Form{T}, obj::State_Space_Objective{T}, tol::T, maxiters::S) where {T <: AbstractFloat, S <: Int}

  nx    = copy(model.nx)
  ny    = copy(model.ny)
  beta  = copy(obj.beta)
  sigma = copy(model.sigma)

  a = copy(model.a)
  b = copy(model.b)
  c = copy(model.c)

  q = copy(obj.q)
  u = copy(obj.u)
  r = copy(obj.r)

  n = nx+ny
  np = size(b,2)
  ns = size(c,2)

  v = q+(a+a')/n
  f = zeros(np,n)

  len = Inf
  iters = 0
  while len > tol && iters < maxiters

    vn = q-u*f-f'u'+f'r*f+beta*(a-b*f)'v*(a-b*f)
    fn = (r+beta*b'vn*b)\(u'+beta*b'vn*a)

    len = maximum(maximum(abs,vn-v), maximum(abs,fn-f))
    f = fn
    v = vn

    iters += 1

  end

  if iters <= maxiters
    retcode = true
  else
    retcode = false
  end

  t = eye(n)
  t[(nx+1):n,1:nx]     = v[(nx+1):n,1:nx]
  t[(nx+1):n,(nx+1):n] = v[(nx+1):n,(nx+1):n]

  p = t*(a-b*f)/t
  k = t*[c[1:nx,:]; -(v[(nx+1):n,(nx+1):n])\v[(nx+1):n,1:nx]*c[1:nx,:]]
  h = inv(t)
  v = h'v*h
  f = -f*h

  # Now re-express the solution by eliminating the Lagrange multipliers

  p11 = p[1:nx,1:nx]
  p12 = p[1:nx,(nx+1):n]
  p21 = p[(nx+1):n,1:nx]
  p22 = p[(nx+1):n,(nx+1):n]

  k11 = k[1:nx,:]

  h21 = h[(nx+1):n,1:nx]
  h22 = h[(nx+1):n,(nx+1):n]

  f1 = f[1:np,1:nx]
  f2 = f[1:np,(nx+1):n]

  l = [h22; f2]
  g = [h21; f1]

  invl = pinv(l)

  p0                             = eye(2*nx+ny+np)
  p0[1:nx,(2*nx+1):(2*nx+ny+np)] = -p12*invl

  p1                                              = zeros(2*nx+ny+np,2*nx+ny+np)
  p1[1:nx,1:nx]                                   = p11-p12*invl*g
  p1[(nx+1):(2*nx),1:nx]                          = eye(nx)
  p1[(2*nx+1):(2*nx+ny+np),1:nx]                  = g
  p1[(2*nx+1):(2*nx+ny+np),(nx+1):(2*nx)]         = l*(p21-p22*invl*g)
  p1[(2*nx+1):(2*nx+ny+np),(2*nx+1):(2*nx+ny+np)] = l*p22*invl

  p = p0\p1

  h = p[(nx+1):(2*nx+ny),:]
  f = p[(2*nx+ny+1):(2*nx+ny+np),:]
  k = [k11; zeros((nx+ny+np),ns)]

  soln = State_Space_Soln(p,k,h,f,v,sigma,retcode)

  return soln

end

function solve_timeless(model::Generalized_State_Space_Form{T}, obj::State_Space_Objective{T}, tol::T, maxiters::S) where {T <: AbstractFloat, S <: Int}

  nx    = copy(model.nx)
  ny    = copy(model.ny)
  beta  = copy(obj.beta)
  sigma = copy(model.sigma)

  a0 = copy(model.a0)
  a  = copy(model.a)
  b  = copy(model.b)
  c  = copy(model.c)

  q = copy(obj.q)
  u = copy(obj.u)
  r = copy(obj.r)

  n  = nx+ny
  np = size(b,2)
  ns = size(c,2)

  q11 = q[1:nx,1:nx]
  q12 = q[1:nx,(nx+1):n]
  q21 = q[(nx+1):n,1:nx]
  q22 = q[(nx+1):n,(nx+1):n]

  u1 = u[1:nx,:]
  u2 = u[(nx+1):n,:]

  a11 = a[1:nx,1:nx]
  a12 = a[1:nx,(nx+1):n]
  a21 = a[(nx+1):n,1:nx]
  a22 = a[(nx+1):n,(nx+1):n]

  b1 = b[1:nx,:]
  b2 = b[(nx+1):n,:]

  a_tilda = [[a11 zeros(nx,ny)]; zeros(ny,n)]
  b_tilda = [[a12 b1 zeros(nx,ny)]; [zeros(ny,ny+np) eye(ny)]]
  c_tilda = c

  q_tilda = [[q11 zeros(nx,ny)]; zeros(ny,n)]
  u_tilda = [[q12 u1 a21']; [((-1/beta)*a0) zeros(ny,ny+np)]]
  r_tilda = [[q22 u2 a22']; [u2' r b2']; [a22 b2 zeros(ny,ny)]]

  v = q+(a+a')/n

  len = Inf
  iters = 0
  while len > tol && iters < maxiters

    f  = -(r_tilda+beta*b_tilda'v*b_tilda)\(u_tilda'+beta*b_tilda'v*a_tilda)
    vn = q_tilda+u_tilda*f+f'u_tilda'+f'r_tilda*f+beta*(a_tilda+b_tilda*f)'v*(a_tilda+b_tilda*f)

    len   = maximum(abs,vn-v)
    v     = vn

    iters += 1

  end

  if iters <= maxiters
    retcode = true
  else
    retcode = false
  end

  q_orig = q_tilda
  u_orig = [[q12 u1]; zeros(ny,ny+np)]
  r_orig = [[q22 u2]; [u2' r]]

  p = a_tilda+b_tilda*f
  v = dlyap(sqrt(beta)*p',q_orig+u_orig[:,1:(ny+np)]*f[1:(ny+np),:]+f[1:(ny+np),:]'u_orig[:,1:(ny+np)]'+f[1:(ny+np),:]'r_orig[1:(ny+np),1:(ny+np)]*f[1:(ny+np),:])
  k = c_tilda
  h = [[eye(nx) zeros(nx,ny)]; f[1:ny,:]]
  f = f[(ny+1):(ny+np),:]

  # Now re-express the solution by eliminating the multipliers

  p11 = p[1:nx,1:nx]
  p12 = p[1:nx,(nx+1):n]
  p21 = p[(nx+1):n,1:nx]
  p22 = p[(nx+1):n,(nx+1):n]

  k11 = k[1:nx,:]

  h21 = h[(nx+1):n,1:nx]
  h22 = h[(nx+1):n,(nx+1):n]

  f1 = f[1:np,1:nx]
  f2 = f[1:np,(nx+1):n]

  l = [h22; f2]
  g = [h21; f1]

  invl = pinv(l)

  p0                             = eye(2*nx+ny+np)
  p0[1:nx,(2*nx+1):(2*nx+ny+np)] = -p12*invl

  p1                                              = zeros(2*nx+ny+np,2*nx+ny+np)
  p1[1:nx,1:nx]                                   = p11-p12*invl*g
  p1[(nx+1):(2*nx),1:nx]                          = eye(nx)
  p1[(2*nx+1):(2*nx+ny+np),1:nx]                  = g
  p1[(2*nx+1):(2*nx+ny+np),(nx+1):(2*nx)]         = l*(p21-p22*invl*g)
  p1[(2*nx+1):(2*nx+ny+np),(2*nx+1):(2*nx+ny+np)] = l*p22*invl

  p = p0\p1

  h = p[(nx+1):(2*nx+ny),:]
  f = p[(2*nx+ny+1):(2*nx+ny+np),:]
  k = [k11; zeros((nx+ny+np),ns)]

  soln = State_Space_Soln(p,k,h,f,v,sigma,retcode)

  return soln

end


text:
# Various functions that work with learners.
module Util

importall Orchestra.Types
import MLBase: Kfold

export holdout,
       kfold,
       score,
       infer_eltype,
       nested_dict_to_tuples,
       nested_dict_set!,
       nested_dict_merge,
       create_transformer

# Holdout method that partitions a collection
# into two partitions.
#
# @param n Size of collection to partition.
# @param right_prop Percentage of collection placed in right partition.
# @return Two partitions of indices, left and right.
function holdout(n, right_prop)
  shuffled_indices = randperm(n)
  partition_pivot = int(right_prop * n)
  right = shuffled_indices[1:partition_pivot]
  left = shuffled_indices[partition_pivot+1:end]
  return (left, right)
end

# Returns k-fold partitions.
#
# @param num_instances Total number of instances.
# @param num_partitions Number of partitions required.
# @return Returns training set partition.
function kfold(num_instances, num_partitions)
  return collect(Kfold(num_instances, num_partitions))
end

# Score learner predictions against ground truth values.
#
# Available metrics:
# - :accuracy
#
# @param metric Metric to assess with.
# @param actual Ground truth values.
# @param predicted Predicted values.
# @return Score of learner.
function score(metric::Symbol, actual, predicted)
  if metric == :accuracy
    mean(actual .== predicted) * 100.0
  else
    error("Metric $metric not implemented for score.")
  end
end

# Returns element type of vector unless it is Any.
# If Any, returns the most specific type that can be
# inferred from the vector elements.
#
# @param vector Vector to infer element type on.
# @return Inferred element type.
function infer_eltype(vector::Vector)
  # Obtain element type of vector
  vec_eltype = eltype(vector)

  # If element type of Vector is Any and not empty,
  # and all elements are of the same type,
  # then return the vector elements' type.
  if vec_eltype == Any && !isempty(vector)
    all_elements_same_type = all(x -> typeof(x) == typeof(first(vector)), vector)
    if all_elements_same_type
      vec_eltype = typeof(first(vector))
    end
  end

  # Return inferred element type
  return vec_eltype
end

# Converts nested dictionary to list of tuples
#
# @param dict Dictionary that can have other dictionaries as values.
# @return List where elements are ([outer-key, inner-key, ...], value).
function nested_dict_to_tuples(dict::Dict)
  set = Set()
  for (entry_id, entry_val) in dict
    if typeof(entry_val) <: Dict
      inner_set = nested_dict_to_tuples(entry_val)
      for (inner_entry_id, inner_entry_val) in inner_set
        new_entry = (vcat([entry_id], inner_entry_id), inner_entry_val)
        push!(set, new_entry)
      end
    else
      new_entry = ([entry_id], entry_val)
      push!(set, new_entry)
    end
  end
  return set
end

# Set value in a nested dictionary.
#
# @param dict Nested dictionary to assign value.
# @param keys Keys to access nested dictionaries in sequence.
# @param value Value to assign.
function nested_dict_set!{T}(dict::Dict, keys::Array{T, 1}, value)
  inner_dict = dict
  for key in keys[1:end-1]
    inner_dict = inner_dict[key]
  end
  inner_dict[keys[end]] = value
end

# Second nested dictionary is merged into first.
#
# If a second dictionary's value as well as the first
# are both dictionaries, then a merge is conducted between
# the two inner dictionaries.
# Otherwise the second's value overrides the first.
#
# @param first First nested dictionary.
# @param second Second nested dictionary.
# @return Merged nested dictionary.
function nested_dict_merge(first::Dict, second::Dict)
  target = copy(first)
  for (second_key, second_value) in second
    values_both_dict =
      typeof(second_value) <: Dict &&
      typeof(get(target, second_key, nothing)) <: Dict
    if values_both_dict
      target[second_key] = nested_dict_merge(target[second_key], second_value)
    else
      target[second_key] = second_value
    end
  end
  return target
end

# Create transformer
#
# @param prototype Prototype transformer to base new transformer on.
# @param options Additional options to override prototype's options.
# @return New transformer.
function create_transformer(prototype::Transformer, options=nothing)
  new_options = copy(prototype.options)
  if options != nothing
    new_options = nested_dict_merge(new_options, options)
  end

  prototype_type = typeof(prototype)
  return prototype_type(new_options)
end

end # module


text:
const _pkg_root = dirname(dirname(@__FILE__))
const _pkg_deps = joinpath(_pkg_root, "deps")
const _pkg_assets = joinpath(_pkg_root, "assets")

!isdir(_pkg_assets) && mkdir(_pkg_assets)

include("make_schema_docs.jl")


text:
abstract type InitializationStrategyCombined <: InitializationStrategy end

struct SimpleCombinedStrategy <: InitializationStrategyCombined
    gamma_strategy::InitializationStrategyGamma
    C_strategy::InitializationStrategyC
end

FixedParameterInitialization(kernel, C) = SimpleCombinedStrategy(FixedGammaStrategy(kernel), FixedCStrategy(C))

function get_parameters(model, strategy::SimpleCombinedStrategy)
    C = calculate_C(model, strategy.C_strategy)
    gamma = calculate_gamma(model, strategy.gamma_strategy)
    return (C, MLKernels.GaussianKernel(gamma))
end

struct GammaFirstCombinedStrategy <: InitializationStrategyCombined
    gamma_strategy::InitializationStrategyGamma
    C_strategy::InitializationStrategyC
end

function get_parameters(model, strategy::GammaFirstCombinedStrategy)
    gamma = calculate_gamma(model, strategy.gamma_strategy)
    C = calculate_C(model, strategy.C_strategy, MLKernels.GaussianKernel(gamma))
    return (C, MLKernels.GaussianKernel(gamma))
end


text:
# # Polarizability using automatic differentiation
#
# Simple example for computing properties using (forward-mode)
# automatic differentation.
# For a more classical approach and more details about computing polarizabilities,
# see [Polarizability by linear response](@ref).

using DFTK
using LinearAlgebra
using ForwardDiff

## Construct PlaneWaveBasis given a particular electric field strength
## Again we take the example of a Helium atom.
function make_basis(ε::T; a=10., Ecut=30) where T
    lattice=T(a) * I(3)  # lattice is a cube of ``a`` Bohrs
    He = ElementPsp(:He, psp=load_psp("hgh/lda/He-q2"))
    atoms = [He => [[1/2; 1/2; 1/2]]]  # Helium at the center of the box

    model = model_DFT(lattice, atoms, [:lda_x, :lda_c_vwn];
                      extra_terms=[ExternalFromReal(r -> -ε * (r[1] - a/2))],
                      symmetries=false)
    PlaneWaveBasis(model; Ecut, kgrid=[1, 1, 1])  # No k-point sampling on isolated system
end

## dipole moment of a given density (assuming the current geometry)
function dipole(basis, ρ)
    @assert isdiag(basis.model.lattice)
    a  = basis.model.lattice[1, 1]
    rr = [a * (r[1] - 1/2) for r in r_vectors(basis)]
    sum(rr .* ρ) * basis.dvol
end

## Function to compute the dipole for a given field strength
function compute_dipole(ε; tol=1e-8, kwargs...)
    scfres = self_consistent_field(make_basis(ε; kwargs...), tol=tol)
    dipole(scfres.basis, scfres.ρ)
end;

# With this in place we can compute the polarizability from finite differences
# (just like in the previous example):
polarizability_fd = let
    ε = 0.01
    (compute_dipole(ε) - compute_dipole(0.0)) / ε
end

# ## Forward-mode implicit differentiation
#
# Right now DFTK has no out-of-the-box support for implicit differentiation through the SCF.
# However one can easily work around this as follows. We keep both a non-dual basis
# and a basis including duals for easy bookkeeping (but redundant computation ...).

function self_consistent_field_dual(basis::PlaneWaveBasis, basis_dual::PlaneWaveBasis{T};
                                    kwargs...) where T <: ForwardDiff.Dual
    scfres = self_consistent_field(basis; kwargs...)
    ψ, occupation = DFTK.select_occupied_orbitals(basis, scfres.ψ, scfres.occupation)

    ## promote everything eagerly to Dual numbers
    occupation_dual = [T.(occupation[1])]
    ψ_dual = [Complex.(T.(real(ψ[1])), T.(imag(ψ[1])))]
    ρ_dual = compute_density(basis_dual, ψ_dual, occupation_dual)

    _, δH = energy_hamiltonian(basis_dual, ψ_dual, occupation_dual; ρ=ρ_dual)
    δHψ = δH * ψ_dual
    δHψ = [ForwardDiff.partials.(δHψ[1], 1)]
    δψ = DFTK.solve_ΩplusK(basis, ψ, -δHψ, occupation)
    δρ = DFTK.compute_δρ(basis, ψ, δψ, occupation)
    ρ = ForwardDiff.value.(ρ_dual)
    ψ, ρ, δψ, δρ
end;

# This function is now used in the following to provide a dual version
# for the compute_dipole function:
function compute_dipole(ε::ForwardDiff.Dual; tol=1e-8, kwargs...)
    T = ForwardDiff.tagtype(ε)
    basis = make_basis(ForwardDiff.value(ε); kwargs...)
    basis_dual = make_basis(ε; kwargs...)
    ψ, ρ, δψ, δρ = self_consistent_field_dual(basis, basis_dual; tol)
    ρ_dual = ForwardDiff.Dual{T}.(ρ, δρ)
    dipole(basis_dual, ρ_dual)
end;

# This setup allows to compute the polarizability via automatic differentiation:
polarizability = ForwardDiff.derivative(compute_dipole, 0.0)
println()
println("Polarizability via ForwardDiff:       $polarizability")
println("Polarizability via finite difference: $polarizability_fd")


text:
using Documenter
#include("../src/QuantumWalk.jl")
using QuantumWalk, LightGraphs
using SparseArrays,Arpack,LinearAlgebra

# same for contributing and license
cp(normpath(@__FILE__, "../../LICENSE.md"), normpath(@__FILE__, "../src/license.md"), force=true)

makedocs(
    modules     = [QuantumWalk],
    format      = Documenter.HTML(assets = ["assets/logo.ico"]),
    sitename    = "QuantumWalk",
    clean       = true,
    doctest     = true,
    checkdocs   = :exports,
    pages       = Any[
		"Home"             => "index.md",
		"Type hierarchies" => "type_hierarchy.md",
    		"Dynamics"	   => Any[
                "Quantum walk evolution" => "quantum_walk.md",
				"Quantum search"         => "quantum_search.md"
                ],
			"Models"       => Any[
				"CTQW and CTQWDense" 	 => "ctqw.md",
				"Szegedy" 		         => "szegedy.md"
                ],
		"Contributing"	   => "contributing.md",
        	"Citing"	       => "citing.md",
		"Licence"	       => "license.md",
    ]
)

deploydocs(
    deps        = nothing,
    make        = nothing,
    repo        = "github.com/iitis/QuantumWalk.jl",
    target      = "build",
)

#rm(normpath(@__FILE__, "src/license.md"))


text:
#jl #! format: off
# # Add time series to `System`

# **Originally Contributed by**: Clayton Barrows

# ## Introduction

# An example of how to parse add time series data to a `System` using [PowerSystems.jl](github.com/NREL-SIIP/PowerSystems.jl)
#
# For example, a `System` created by [parsing a MATPOWER file](https://nbviewer.jupyter.org/github/NREL-SIIP/SIIPExamples.jl/blob/master/notebook/2_PowerSystems_examples/parse_matpower.ipynb)
# doesn't contain any time series data. So a user may want to add time series to the `System`
# ### Dependencies
# Let's use the 5-bus dataset we parsed in the MATPOWER example
using PowerSystems
using TimeSeries
using JSON3

base_dir = PowerSystems.download(PowerSystems.TestData; branch = "master");

sys = System(joinpath(base_dir, "matpower", "case5_re.m"))
sys

# ### Define pointers to time series files
# For example, if we want to add a bunch of time series files, say one for each load and
# one for each renewable generator, we need to define pointers to each .csv file containing
# the time series in the following format (PowerSystems.jl also supports a CSV format for this file)

base_dir = PowerSystems.download(PowerSystems.TestData; branch = "master");
FORECASTS_DIR = joinpath(base_dir, "forecasts", "5bus_ts")
fname = joinpath(FORECASTS_DIR, "timeseries_pointers_da.json")
open(fname, "r") do f
    JSON3.@pretty JSON3.read(f)
end

# ### Read and assign time series to `System` using these parameters.
add_time_series!(sys, fname)
sys

max_stars_repo_name: NREL-SIIP/SIIP-Tutorial

max_issues_repo_name: NREL-SIIP/SIIP-Tutorial

max_forks_repo_name: NREL-SIIP/SIIP-Tutorial


text:
@inline function check_simulation_stop!(sol, t::Float64, Y, YP, run::AbstractRun, p, bounds, opts::options_simulation_immutable{T};
    ϵ::Float64 = t < 1.0 ? opts.reltol : 0.0,
    ) where T<:Function

    if t ≥ run.tf
        run.info.flag = 0
        run.info.exit_reason = "Final time reached"
        return nothing
    end

    # continue checking the bounds or return after only evaluating the time
    !opts.check_bounds && (return nothing)

    I = calc_I(Y,p)

    check_stop_I(         p, run, sol, Y, YP, bounds, ϵ, I)
    check_stop_V(         p, run, sol, Y, YP, bounds, ϵ, I)
    check_stop_SOC(       p, run, sol, Y, YP, bounds, ϵ, I)
    check_stop_T(         p, run, sol, Y, YP, bounds, ϵ, I)
    check_stop_c_s_surf(  p, run, sol, Y, YP, bounds, ϵ, I)
    check_stop_c_e(       p, run, sol, Y, YP, bounds, ϵ, I)
    check_stop_η_plating( p, run, sol, Y, YP, bounds, ϵ, I)
    check_stop_dfilm(     p, run, sol, Y, YP, bounds, ϵ, I)
    opts.stop_function(   p, run, sol, Y, YP, bounds, ϵ, I)

    return nothing
end

@inline check_stop_I(p::R4, run::R3, sol, Y::R2, YP::R2, bounds::R5, ϵ::Float64, I::Float64
    ) where {R2<:Vector{Float64}, R3<:AbstractRun{method_I},R4<:model,R5<:boundary_stop_conditions_immutable} = nothing
@inline function check_stop_I(p::R4, run::R3, sol, Y::R2, YP::R2, bounds::R5, ϵ::Float64, I::Float64
    ) where {R2<:Vector{Float64}, method, R3<:AbstractRun{method},R4<:model,R5<:boundary_stop_conditions_immutable}

    if (I - bounds.I_max > ϵ)
        t_frac = (bounds.prev.I - bounds.I_max)/(bounds.prev.I - I)
        if t_frac < bounds.prev.t_final_interp_frac
            bounds.prev.t_final_interp_frac = t_frac
            run.info.flag = 7
            run.info.exit_reason = "Above maximum permitted C-rate"
        end
    elseif (bounds.I_min - I > ϵ)
        t_frac = (bounds.prev.I - bounds.I_min)/(bounds.prev.I - I)
        if t_frac < bounds.prev.t_final_interp_frac
            bounds.prev.t_final_interp_frac = t_frac
            run.info.flag = 8
            run.info.exit_reason = "Below minimum permitted C-rate"
        end
    end
    bounds.prev.I = I

    return nothing
end

@inline check_stop_V(p::R4, run::R3, sol, Y::R2, YP::R2, bounds::R5, ϵ::Float64, I::Float64
) where {R2<:Vector{Float64}, R3<:AbstractRun{method_V},R4<:model,R5<:boundary_stop_conditions_immutable} = nothing
@inline function check_stop_V(p::R4, run::R3, sol, Y::R2, YP::R2, bounds::R5, ϵ::Float64, I::Float64
    ) where {R2<:Vector{Float64}, method, R3<:AbstractRun{method},R4<:model,R5<:boundary_stop_conditions_immutable}

    V = calc_V(Y,p)
    if (bounds.V_min - V > ϵ) && I < 0
        t_frac = (bounds.prev.V - bounds.V_min)/(bounds.prev.V - V)
        if t_frac < bounds.prev.t_final_interp_frac
            bounds.prev.t_final_interp_frac = t_frac
            run.info.flag = 1
            run.info.exit_reason = "Below minimum voltage limit"
        end
    elseif (V - bounds.V_max > ϵ) && I > 0
        t_frac = (bounds.prev.V - bounds.V_max)/(bounds.prev.V - V)
        if t_frac < bounds.prev.t_final_interp_frac
            bounds.prev.t_final_interp_frac = t_frac
            run.info.flag = 2
            run.info.exit_reason = "Above maximum voltage limit"
        end
    end
    bounds.prev.V = V

    return nothing
end

@inline function check_stop_SOC(p::R4, run::R3, sol::solution, Y::R2, YP::R2, bounds::R5, ϵ::Float64, I::Float64
    ) where {R2<:Vector{Float64}, R3<:AbstractRun,R4<:model,R5<:boundary_stop_conditions_immutable}

    SOC = @inbounds sol.SOC[end]
    if (bounds.SOC_min - SOC > ϵ) && I < 0
        t_frac = (bounds.prev.SOC - bounds.SOC_min)/(bounds.prev.SOC - SOC)
        if t_frac < bounds.prev.t_final_interp_frac
            bounds.prev.t_final_interp_frac = t_frac
            run.info.flag = 3
            run.info.exit_reason = "Below minimum SOC limit"
        end
    elseif (SOC - bounds.SOC_max > ϵ) && I > 0
        t_frac = (bounds.prev.SOC - bounds.SOC_max)/(bounds.prev.SOC - SOC)
        if t_frac < bounds.prev.t_final_interp_frac
            bounds.prev.t_final_interp_frac = t_frac
            run.info.flag = 4
            run.info.exit_reason = "Above maximum SOC limit"
        end
    end
    bounds.prev.SOC = SOC

    return SOC
end

@inline check_stop_T(p::model_temp{false}, run, sol, Y, YP, bounds, ϵ, I) = nothing
@inline function check_stop_T(p::R4, run::R3, sol, Y::R2, YP::R2, bounds::R5, ϵ::Float64, I::Float64
    ) where {R2<:Vector{Float64}, R3<:AbstractRun,R4<:model_temp{true},R5<:boundary_stop_conditions_immutable}

    if !isnan(bounds.T_max)
        T = temperature_weighting(calc_T(Y,p),p)
        if T - bounds.T_max > ϵ
            t_frac = (bounds.prev.T - bounds.T_max)/(bounds.prev.T - T)
            if t_frac < bounds.prev.t_final_interp_frac
                bounds.prev.t_final_interp_frac = t_frac
                run.info.flag = 5
                run.info.exit_reason = "Above maximum permitted temperature"
            end
        end
        bounds.prev.T = T
    end

    return nothing
end

@inline function c_s_n_maximum(Y::Vector{Float64},p::model_solid_diff{:Fickian})
    c_s_n_max = -Inf
    @inbounds for i in 1:p.N.n
        @inbounds c_s_n_max = max(c_s_n_max, Y[p.ind.c_s_avg.n[p.N.r_n*i]])
    end
    return c_s_n_max
end
@inline function c_s_n_maximum(Y::Vector{Float64},p::Union{model_solid_diff{:polynomial},model_solid_diff{:quadratic}})
    c_s_n_max = -Inf
    @inbounds for ind in p.ind.c_s_avg.n
        @inbounds c_s_n_max = max(c_s_n_max, Y[ind])
    end
    return c_s_n_max
end

@inline function check_stop_c_s_surf(p::R4, run::R3, sol, Y::R2, YP::R2, bounds::R5, ϵ::Float64, I::Float64
    ) where {R2<:Vector{Float64}, R3<:AbstractRun,R4<:model,R5<:boundary_stop_conditions_immutable}

    if !isnan(bounds.c_s_n_max)
        c_s_n_max = c_s_n_maximum(Y,p)

        if I > 0
            if c_s_n_max - bounds.c_s_n_max*p.θ[:c_max_n] > ϵ
                t_frac = (bounds.prev.c_s_n - bounds.c_s_n_max*(p.θ[:c_max_n]))/(bounds.prev.c_s_n - c_s_n_max)
                if t_frac < bounds.prev.t_final_interp_frac
                    bounds.prev.t_final_interp_frac = t_frac
                    run.info.flag = 6
                    run.info.exit_reason = "Above c_s_n saturation threshold"
                end
            end
        end
        bounds.prev.c_s_n = c_s_n_max
    end

    return nothing
end

@inline function check_stop_c_e(p::R4, run::R3, sol, Y::R2, YP::R2, bounds::R5, ϵ::Float64, I::Float64
    ) where {R2<:Vector{Float64}, R3<:AbstractRun, R4<:model, R5<:boundary_stop_conditions_immutable}

    if !isnan(bounds.c_e_min)
        c_e_min = +Inf
        @inbounds for ind in p.ind.c_e
            @inbounds c_e_min = min(c_e_min, Y[ind])
        end
        if bounds.c_e_min - c_e_min > ϵ
            t_frac = (bounds.prev.c_e_min - bounds.c_e_min)/(bounds.prev.c_e_min - c_e_min)
            if t_frac < bounds.prev.t_final_interp_frac
                bounds.prev.t_final_interp_frac = t_frac
                run.info.flag = 9
                run.info.exit_reason = "Below minimum permitted c_e"
            end
        end
        bounds.prev.c_e_min = c_e_min
    end

    return nothing
end

@inline function check_stop_η_plating(p::R4, run::R3, sol, Y::R2, YP::R2, bounds::R5, ϵ::Float64, I::Float64
    ) where {R2<:Vector{Float64}, R3<:AbstractRun, R4<:model, R5<:boundary_stop_conditions_immutable}

    η_plating = calc_η_plating(Y,p)
    if !isnan(bounds.η_plating_min) && bounds.η_plating_min - η_plating > ϵ
        t_frac = (bounds.prev.η_plating - bounds.η_plating_min)/(bounds.prev.η_plating - η_plating)
        if t_frac < bounds.prev.t_final_interp_frac
            bounds.prev.t_final_interp_frac = t_frac
            run.info.flag = 9
            run.info.exit_reason = "Below minimum permitted η_plating"
        end
    end
    bounds.prev.η_plating = η_plating

    return nothing
end

@inline check_stop_dfilm(::model_age{false}, run, sol, Y, YP, bounds, ϵ, I) = nothing
@inline function check_stop_dfilm(p::R4, run::R3, sol, Y::R2, YP::R2, bounds::R5, ϵ::Float64, I::Float64
    ) where {R2<:Vector{Float64}, R3<:AbstractRun, R4<:model_age{:SEI}, R5<:boundary_stop_conditions_immutable}

    dfilm_max = -Inf
    @inbounds for i in 1:p.N.n
        @inbounds dfilm_max = max(dfilm_max, YP[p.ind.film[i]])
    end

    if !isnan(bounds.dfilm_max) && dfilm_max - bounds.dfilm_max > ϵ
        t_frac = (bounds.prev.dfilm - bounds.dfilm_max)/(bounds.prev.dfilm - dfilm_max)
        if t_frac < bounds.prev.t_final_interp_frac
            bounds.prev.t_final_interp_frac = t_frac
            run.info.flag = 10
            run.info.exit_reason = "Above maximum film growth rate"
        end
    end
    bounds.prev.dfilm = dfilm_max

    return nothing
end

@inline function check_solve(run::run_constant, sol::R1, int::R2, p, bounds, opts::R5, funcs, keep_Y::Bool, iter::Int64, Y::Vector{Float64}, t::Float64) where {R1<:solution,R2<:Sundials.IDAIntegrator,R5<:AbstractOptionsModel}
    if t == int.tprev
        # Sometimes the initial step at t = 0 can be too large. This reduces the step size
        if t == 0.0
            if iter == 2
                Sundials.IDASetInitStep(int.mem,opts.reltol)
            else
                error("Model failed to converge at t = $(t). Try tightening the absolute and relative tolerances.")
            end
        else
            error("Model failed to converge at t = $(t)")
        end
    elseif iter == opts.maxiters
        error("Reached max iterations of $(opts.maxiters) at t = $(t)")
    elseif within_bounds(run)
        # update Y only after checking the stop conditions. this is done to store a copy of the
        # previous sol run in case any back-interpolation is needed
        set_var!(sol.Y, keep_Y ? copy(Y) : Y, keep_Y)
    else # no errors and run.info.flag ≠ -1
        return false
    end

    return true
end

@inline function check_solve(run::run_function, sol::R1, int::R2, p::model, bounds::boundary_stop_conditions_immutable, opts::R5, funcs, keep_Y::Bool, iter::Int64, Y::Vector{Float64}, t::Float64) where {R1<:solution,R2<:Sundials.IDAIntegrator,R5<:AbstractOptionsModel}
    if iter == opts.maxiters
        error("Reached max iterations of $(opts.maxiters) at t = $(int.t)")
    elseif within_bounds(run)
        # update Y only after checking the stop conditions. this is done to store a copy of the
        # previous sol run in case any back-interpolation is needed
        set_var!(sol.Y, keep_Y ? copy(Y) : Y, keep_Y)

        # check to see if the run needs to be reinitialized
        if t - int.tprev < 1e-3opts.reltol
            check_reinitialization!(sol, int, run, p, bounds, opts, funcs)
        end

        return true
    else # no errors and run.info.flag ≠ -1
        return false
    end
end

@inline function check_reinitialization!(sol::R1, int::R2, run::R3, p::R4, bounds::R5, opts::R6, funcs) where {R1<:solution, R2<:Sundials.IDAIntegrator, R3<:AbstractRun,R4<:model,R5<:boundary_stop_conditions_immutable,R6<:AbstractOptionsModel}
    """
    Checking the current function for discontinuities.
    If there is a significant change in current after a step size of dt = reltol,
    then rerun Newton's method and reinitialize at this new time step.
    """

    Y = int.u.v
    YP = int.du.v
    # take a step of Δt = the relative tolerance
    t_new = int.t + opts.reltol

    value_old = value(run)
    value_new = run.func(t_new,Y,YP,p)

    # if the function values at t vs. t + Δt are very different (i.e., there is a discontinuity)
    # then reinitialize the DAE at t + Δt
    if !≈(value_old, value_new, atol=opts.abstol, rtol=opts.reltol)
        initialize_states!(p,Y,YP,run,opts,funcs,(@inbounds sol.SOC[end]); t=t_new)

        Sundials.IDAReInit(int.mem, t_new, Y, YP)
    end
    return nothing
end

@inline function check_errors_parameters_runtime(p::R1,opts::R2) where {R1<:model,R2<:AbstractOptionsModel}
    ϵ_sp, ϵ_sn = active_material(p)

    if ( ϵ_sp > 1 ) error("ϵ_p + ϵ_fp must be ∈ [0, 1)") end
    if ( ϵ_sn > 1 ) error("ϵ_n + ϵ_fn must be ∈ [0, 1)") end
    if ( p.θ[:θ_max_p] > p.θ[:θ_min_p] ) error("θ_max_p must be < θ_min_p") end
    if ( p.θ[:θ_min_n] > p.θ[:θ_max_n] ) error("θ_min_n must be < θ_max_n") end

    return nothing
end

function check_errors_initial(θ, numerics, N)
    if !(numerics.jacobian ∈ (:symbolic, :AD))
        error("`jacobian` can either be :symbolic or :AD")
    end

    return nothing
end

check_is_hold(x::Symbol,sol::solution) = (x == :hold) && (!isempty(sol) ? true : error("Cannot use `:hold` without a previous sol."))
check_is_hold(::Any,::solution) = false


text:
"""
    PrescribedConditions{T}

Describes the forces, moments, displacements, and/or rotations prescribed at a
point for each time step
"""
struct PrescribedConditions{T}
    force::Vector{Bool}
    value::Vector{SVector{6, T}}
    follower::Vector{SVector{6, T}}
end
Base.eltype(::PrescribedConditions{T}) where T = T

"""
    PrescribedConditions(dt=0.0; kwargs...)

Construct an object of type PrescribedConditions which stores the prescribed
conditions for a point at each time step.

Prescribed conditions may be assigned as either a scalar parameter or as a
function of time.

Prescribed Wiener-Milenkovic parameters must satisfy the following inequality:
sqrt(theta_x^2 + theta_y^2 + theta_z^2) <= 4.  Note that this restriction still
allows all possible rotations to be represented.

# Arguments
 - `dt`: Time step size.
 - `nstep`: The total length of the time vector
 - `ux`: Prescribed x-direction displacement of the point
 - `uy`: Prescribed y-direction displacement of the point
 - `uz`: Prescribed z-direction displacement of the point
 - `theta_x`: Prescribed first Wiener-Milenkovic parameter of the point
 - `theta_y`: Prescribed second Wiener-Milenkovic parameter of the point
 - `theta_z`: Prescribed third Wiener-Milenkovic parameter of the point
 - `Fx`: Prescribed force in x-direction applied on the point
 - `Fy`: Prescribed force in y-direction applied on the point
 - `Fz`: Prescribed force in z-direction applied on the point
 - `Mx`: Prescribed moment about x-axis applied on the point
 - `My`: Prescribed moment about y-axis applied on the point
 - `Mz`: Prescribed moment about z-axis applied on the point
 - `Fx_follower`: Prescribed follower force in x-direction applied on the point
 - `Fy_follower`: Prescribed follower force in y-direction applied on the point
 - `Fz_follower`: Prescribed follower force in z-direction applied on the point
 - `Mx_follower`: Prescribed follower moment about x-axis applied on the point
 - `My_follower`: Prescribed follower moment about y-axis applied on the point
 - `Mz_follower`: Prescribed follower moment about z-axis applied on the point
"""
function PrescribedConditions(dt=0.0;
    nstep = 1,
    ux = nothing,
    uy = nothing,
    uz = nothing,
    theta_x = nothing,
    theta_y = nothing,
    theta_z = nothing,
    Fx = nothing,
    Fy = nothing,
    Fz = nothing,
    Mx = nothing,
    My = nothing,
    Mz = nothing,
    Fx_follower = nothing,
    Fy_follower = nothing,
    Fz_follower = nothing,
    Mx_follower = nothing,
    My_follower = nothing,
    Mz_follower = nothing)

    force = Vector{Bool}(undef, 6)

    t = range(0, dt*(nstep-1), length=nstep)

    # Set function for first slot
    if isnothing(ux)
        force[1] = true
        v1 = isnothing(Fx) ? (t) -> 0.0 : typeof(Fx) <: Number ? (t) -> Fx : Fx
        v1_follower = isnothing(Fx_follower) ? (t) -> 0.0 : typeof(Fx_follower) <: Number ? (t) -> Fx_follower : Fx
    elseif !isnothing(ux) && (isnothing(Fx) && isnothing(Fx_follower))
        force[1] = false
        v1 = isnothing(ux) ? (t) -> 0.0 : typeof(ux) <: Number ? (t) -> ux : ux
        v1_follower = (t) -> 0.0
    else
        error("Both `ux` and `Fx` or `Fx_follower` cannot be specified at the same time")
    end

    # Set function for second slot
    if isnothing(uy)
        force[2] = true
        v2 = isnothing(Fy) ? (t) -> 0.0 : typeof(Fy) <: Number ? (t) -> Fy : Fy
        v2_follower = isnothing(Fy_follower) ? (t) -> 0.0 : typeof(Fy_follower) <: Number ? (t) -> Fy_follower : Fy
    elseif !isnothing(uy) && (isnothing(Fy) && isnothing(Fy_follower))
        force[2] = false
        v2 = isnothing(uy) ? (t) -> 0.0 : typeof(uy) <: Number ? (t) -> uy : uy
        v2_follower = (t) -> 0.0
    else
        error("Both `uy` and `Fy` or `Fy_follower` cannot be specified at the same time")
    end

    # Set function for third slot
    if isnothing(uz)
        force[3] = true
        v3 = isnothing(Fz) ? (t) -> 0.0 : typeof(Fz) <: Number ? (t) -> Fz : Fz
        v3_follower = isnothing(Fz_follower) ? (t) -> 0.0 : typeof(Fz_follower) <: Number ? (t) -> Fz_follower : Fz
    elseif !isnothing(uz) && (isnothing(Fz) && isnothing(Fz_follower))
        force[3] = false
        v3 = isnothing(uz) ? (t) -> 0.0 : typeof(uz) <: Number ? (t) -> uz : uz
        v3_follower = (t) -> 0.0
    else
        error("Both `uz` and `Fz` or `Fz_follower` cannot be specified at the same time")
    end

    # Set function for fourth slot
    if isnothing(theta_x)
        force[4] = true
        v4 = isnothing(Mx) ? (t) -> 0.0 : typeof(Mx) <: Number ? (t) -> Mx : Mx
        v4_follower = isnothing(Mx_follower) ? (t) -> 0.0 : typeof(Mx_follower) <: Number ? (t) -> Mx_follower : Mx
    elseif !isnothing(theta_x) && (isnothing(Mx) && isnothing(Mx_follower))
        force[4] = false
        v4 = isnothing(theta_x) ? (t) -> 0.0 : typeof(theta_x) <: Number ? (t) -> theta_x : theta_x
        v4_follower = (t) -> 0.0
    else
        error("Both `theta_x` and `Mx` or `Mx_follower` cannot be specified at the same time")
    end

    # Set function for fifth slot
    if isnothing(theta_y)
        force[5] = true
        v5 = isnothing(My) ? (t) -> 0.0 : typeof(My) <: Number ? (t) -> My : My
        v5_follower = isnothing(My_follower) ? (t) -> 0.0 : typeof(My_follower) <: Number ? (t) -> My_follower : My
    elseif !isnothing(theta_y) && (isnothing(My) && isnothing(My_follower))
        force[5] = false
        v5 = isnothing(theta_y) ? (t) -> 0.0 : typeof(theta_y) <: Number ? (t) -> theta_y : theta_y
        v5_follower = (t) -> 0.0
    else
        error("Both `theta_y` and `My` or `My_follower` cannot be specified at the same time")
    end

    # Set function for sixth slot
    if isnothing(theta_z)
        force[6] = true
        v6 = isnothing(Mz) ? (t) -> 0.0 : typeof(Mz) <: Number ? (t) -> Mz : Mz
        v6_follower = isnothing(Mz_follower) ? (t) -> 0.0 : typeof(Mz_follower) <: Number ? (t) -> Mz_follower : Mz
    elseif !isnothing(theta_z) && (isnothing(Mz) && isnothing(Mz_follower))
        force[6] = false
        v6 = isnothing(theta_z) ? (t) -> 0.0 : typeof(theta_z) <: Number ? (t) -> theta_z : theta_z
        v6_follower = (t) -> 0.0
    else
        error("Both `theta_z` and `Mz` or `Mz_follower` cannot be specified at the same time")
    end

    # now populate values for each time step
    value = [
        SVector{6}(
            v1(t[i]), v2(t[i]), v3(t[i]),
            v4(t[i]), v5(t[i]), v6(t[i])
            ) for i in eachindex(t)
    ]

    follower = [
        SVector{6}(
            v1_follower(t[i]), v2_follower(t[i]), v3_follower(t[i]),
            v4_follower(t[i]), v5_follower(t[i]), v6_follower(t[i])
        ) for i in eachindex(t)
    ]

    return PrescribedConditions(force, promote(value, follower)...)
end

"""
    DistributedLoads{T}

Contains the integrated distributed forces and moments for each beam element for each time step.

# Fields
 - f1: Integrated non-follower distributed force for the beam element's left endpoint for each time step
 - f2: Integrated non-follower distributed force for the beam element's right endpoint for each time step
 - m1: Integrated non-follower distributed moment for the beam element's left endpoint for each time step
 - m2: Integrated non-follower distributed moment for the beam element's right endpoint for each time step
 - f1_follower: Integrated follower distributed force for the beam element's left endpoint for each time step
 - f2_follower: Integrated follower distributed force for the beam element's right endpoint for each time step
 - m1_follower: Integrated follower distributed moment for the beam element's left endpoint for each time step
 - m2_follower: Integrated follower distributed moment for the beam element's right endpoint for each time step
"""
struct DistributedLoads{T}
    f1::Vector{SVector{3, T}}
    f2::Vector{SVector{3, T}}
    m1::Vector{SVector{3, T}}
    m2::Vector{SVector{3, T}}
    f1_follower::Vector{SVector{3, T}}
    f2_follower::Vector{SVector{3, T}}
    m1_follower::Vector{SVector{3, T}}
    m2_follower::Vector{SVector{3, T}}
end
Base.eltype(::DistributedLoads{T}) where T = T

"""
    DistributedLoads(assembly, ibeam[, dt]; kwargs...)

Integrates the specified distributed loads over the element for each time step.

# Arguments
 - `assembly`: The beam element assembly
 - `ibeam`: The index of the beam element which the distributed load is assigned to
 - `dt`: Time step size.  If omitted a single time step is assumed and specified
     functions become a function of `s` only.
 - `s1 = 0.0`: Start of beam element (used for integrating the distributed loads)
 - `s2 = 1.0`: End of beam element (used for integrating the distributed loads)
 - `nstep`: The total length of the time vector
 - `method = (f, a, b) -> gauss_quadrature(f, a, b)`: Method which integrates function
    `f` from `a` to `b`. Defaults to the Gauss-Legendre quadrature with 4 points on each element.
 - `fx = (s, t) -> 0.0`: Distributed non-follower force on beam element in x-direction
 - `fy = (s, t) -> 0.0`: Distributed non-follower force on beam element in y-direction
 - `fz = (s, t) -> 0.0`: Distributed non-follower force on beam element in z-direction
 - `mx = (s, t) -> 0.0`: Distributed non-follower moment on beam element in x-direction
 - `my = (s, t) -> 0.0`: Distributed non-follower moment on beam element in y-direction
 - `mz = (s, t) -> 0.0`: Distributed non-follower moment on beam element in z-direction
 - `fx_follower = (s, t) -> 0.0`: Distributed follower force on beam element in x-direction
 - `fy_follower = (s, t) -> 0.0`: Distributed follower force on beam element in y-direction
 - `fz_follower = (s, t) -> 0.0`: Distributed follower force on beam element in z-direction
 - `mx_follower = (s, t) -> 0.0`: Distributed follower moment on beam element in x-direction
 - `my_follower = (s, t) -> 0.0`: Distributed follower moment on beam element in y-direction
 - `mz_follower = (s, t) -> 0.0`: Distributed follower moment on beam element in z-direction
"""
function DistributedLoads(assembly, ibeam;
    s1 = 0.0,
    s2 = 1.0,
    method = gauss_quadrature,
    fx = (s) -> 0.0,
    fy = (s) -> 0.0,
    fz = (s) -> 0.0,
    mx = (s) -> 0.0,
    my = (s) -> 0.0,
    mz = (s) -> 0.0,
    fx_follower = (s) -> 0.0,
    fy_follower = (s) -> 0.0,
    fz_follower = (s) -> 0.0,
    mx_follower = (s) -> 0.0,
    my_follower = (s) -> 0.0,
    mz_follower = (s) -> 0.0,
    )

    return DistributedLoads(assembly, ibeam, 0.0;
        s1 = s1,
        s2 = s2,
        nstep = 1,
        method = method,
        fx = (s, t) -> fx(s),
        fy = (s, t) -> fy(s),
        fz = (s, t) -> fz(s),
        mx = (s, t) -> mx(s),
        my = (s, t) -> my(s),
        mz = (s, t) -> mz(s),
        fx_follower = (s, t) -> fx_follower(s),
        fy_follower = (s, t) -> fy_follower(s),
        fz_follower = (s, t) -> fz_follower(s),
        mx_follower = (s, t) -> mx_follower(s),
        my_follower = (s, t) -> my_follower(s),
        mz_follower = (s, t) -> mz_follower(s),
        )

end

function DistributedLoads(assembly, ibeam, dt;
    s1 = 0.0,
    s2 = 1.0,
    nstep = 1,
    method = gauss_quadrature,
    fx = (s, t) -> 0.0,
    fy = (s, t) -> 0.0,
    fz = (s, t) -> 0.0,
    mx = (s, t) -> 0.0,
    my = (s, t) -> 0.0,
    mz = (s, t) -> 0.0,
    fx_follower = (s, t) -> 0.0,
    fy_follower = (s, t) -> 0.0,
    fz_follower = (s, t) -> 0.0,
    mx_follower = (s, t) -> 0.0,
    my_follower = (s, t) -> 0.0,
    mz_follower = (s, t) -> 0.0,
    )

    # get element length
    ΔL = assembly.elements[ibeam].L

    # create function for general coordinate
    ξ = (s) -> (s-s1)/(s2-s1)

    # create time range
    t = range(0.0, step=dt, length=nstep+1)

    # integrate to get f1 for each time step
    f1 = [
        SVector(
            method((s)->(1-ξ(s))*fx(s, t[i])*ΔL/(s2-s1), s1, s2),
            method((s)->(1-ξ(s))*fy(s, t[i])*ΔL/(s2-s1), s1, s2),
            method((s)->(1-ξ(s))*fz(s, t[i])*ΔL/(s2-s1), s1, s2)
        ) for i in eachindex(t)
    ]

    # integrate to get f2 for each time step
    f2 = [
        SVector(
            method((s)->ξ(s)*fx(s, t[i])*ΔL/(s2-s1), s1, s2),
            method((s)->ξ(s)*fy(s, t[i])*ΔL/(s2-s1), s1, s2),
            method((s)->ξ(s)*fz(s, t[i])*ΔL/(s2-s1), s1, s2)
        ) for i in eachindex(t)
    ]

    # integrate to get m1 for each time step
    m1 = [
        SVector(
            method((s)->(1-ξ(s))*mx(s, t[i])*ΔL/(s2-s1), s1, s2),
            method((s)->(1-ξ(s))*my(s, t[i])*ΔL/(s2-s1), s1, s2),
            method((s)->(1-ξ(s))*mz(s, t[i])*ΔL/(s2-s1), s1, s2)
        ) for i in eachindex(t)
    ]

    # integrate to get m2 for each time step
    m2 = [
        SVector(
            method((s)->ξ(s)*mx(s, t[i])*ΔL/(s2-s1), s1, s2),
            method((s)->ξ(s)*my(s, t[i])*ΔL/(s2-s1), s1, s2),
            method((s)->ξ(s)*mz(s, t[i])*ΔL/(s2-s1), s1, s2)
        ) for i in eachindex(t)
    ]

    # integrate to get f1_follower for each time step
    f1_follower = [
        SVector(
            method((s)->(1-ξ(s))*fx_follower(s, t[i])*ΔL/(s2-s1), s1, s2),
            method((s)->(1-ξ(s))*fy_follower(s, t[i])*ΔL/(s2-s1), s1, s2),
            method((s)->(1-ξ(s))*fz_follower(s, t[i])*ΔL/(s2-s1), s1, s2)
        ) for i in eachindex(t)
    ]

    # integrate to get f2_follower for each time step
    f2_follower = [
        SVector(
            method((s)->ξ(s)*fx_follower(s, t[i])*ΔL/(s2-s1), s1, s2),
            method((s)->ξ(s)*fy_follower(s, t[i])*ΔL/(s2-s1), s1, s2),
            method((s)->ξ(s)*fz_follower(s, t[i])*ΔL/(s2-s1), s1, s2)
        ) for i in eachindex(t)
    ]

    # integrate to get m1_follower for each time step
    m1_follower = [
        SVector(
            method((s)->(1-ξ(s))*mx_follower(s, t[i])*ΔL/(s2-s1), s1, s2),
            method((s)->(1-ξ(s))*my_follower(s, t[i])*ΔL/(s2-s1), s1, s2),
            method((s)->(1-ξ(s))*mz_follower(s, t[i])*ΔL/(s2-s1), s1, s2)
        ) for i in eachindex(t)
    ]

    # integrate to get m2_follower for each time step
    m2_follower = [
        SVector(
            method((s)->ξ(s)*mx_follower(s, t[i])*ΔL/(s2-s1), s1, s2),
            method((s)->ξ(s)*my_follower(s, t[i])*ΔL/(s2-s1), s1, s2),
            method((s)->ξ(s)*mz_follower(s, t[i])*ΔL/(s2-s1), s1, s2)
        ) for i in eachindex(t)
    ]

    return DistributedLoads(f1, f2, m1, m2, f1_follower, f2_follower, m1_follower, m2_follower)
end


text:
using WebIO, JSExpr, GeometryTypes

struct ThreeScreen <: Base.AbstractDisplay
    context::Scope
end

function ThreeScreen()
    threejs = Scope(imports = [
        "//cdnjs.cloudflare.com/ajax/libs/three.js/84/three.min.js",
    ]);

    onimport(threejs, @js function (THREE, OrbitControlsModule)
        @var container = this.dom.querySelector("#container");
        @var scene; @var renderer;
        @var geometry; @var material; @var mesh;
        @var globalscope = this;

        function onWindowResize(event)
            globalscope.camera.aspect = window.innerWidth / window.innerHeight
            globalscope.camera.updateProjectionMatrix()
        end
        function animate()
            requestAnimationFrame(animate)
            render()
        end
        function render()
            renderer.render(scene, globalscope.camera)
        end

        function init()
            renderer = @new THREE.WebGLRenderer()
            window.globalscope = globalscope
            window.THREE = THREE
            globalscope.camera = @new THREE.PerspectiveCamera(50, window.innerWidth / window.innerHeight, 1, 5000);
            globalscope.camera.position.z = 1400;
            scene = @new THREE.Scene()
            window.scene = scene
            renderer.setPixelRatio(window.devicePixelRatio)
            renderer.setSize(window.innerWidth, window.innerHeight)
            container.appendChild(renderer.domElement)
            scene.add(@new THREE.AmbientLight("#444444"));
            @var light1 = @new THREE.DirectionalLight(0xffffff, 0.5);
            light1.position.set(1, 1, 1);
            scene.add(light1);
            window.addEventListener("resize", onWindowResize, false)
        end
        init()
        animate()
    end)
    ThreeScreen(threejs)
end
geom = GLNormalMesh(Sphere(Point3f0(0), 100f0))
gvertices = reinterpret(Float32, decompose(Point3f0, geom)) |> collect
gindices = reinterpret(UInt32, decompose(GLTriangle, geom)) |> collect
gnormals = reinterpret(Float32, decompose(Normal{3, Float32}, geom)) |> collect

using Blink
screen = ThreeScreen()
w = Window()
body!(w, screen.context(dom"div#container"()))
evaljs(screen.context, @js begin
    @var THREE = window.THREE

    @var geometry = @new THREE.BufferGeometry();
    geometry.setIndex($(gindices));
    @var vertices = @new Float32Array($(gvertices));
    @var normals = @new Float32Array($(gnormals));
    # itemSize = 3 because there are 3 values (components) per vertex
    geometry.addAttribute("position", @new THREE.BufferAttribute(vertices, 3))
    geometry.addAttribute("normal", @new THREE.BufferAttribute(normals, 3))

    @var material = @new THREE.MeshBasicMaterial(d(color = "#00ff00"))
    @var mesh = @new THREE.Mesh(geometry, material)
    window.scene.add(mesh)
end)
tools(w)


text:
#= See dataStructuresDefinitions.md for explanation
=#

mutable struct SolutionState
    cellState::Matrix{Float64}
    cellFluxes::Matrix{Float64}
    cellPrimitives::Matrix{Float64}
    fluxResiduals::Matrix{Float64}
    faceFluxes::Matrix{Float64}
end

struct Mesh
    cells::Vector{Vector{Int64}}
    cVols::Vector{Float64}
    cCenters::Vector{Vector{Float64}}
    cellSizes::Matrix{Float64}

    faces::Vector{Vector{Int64}}
    fAVecs::Vector{Vector{Float64}}
    fCenters::Vector{Vector{Float64}}
    boundaryFaces::Vector{Vector{Int64}}
end

struct Fluid
    Cp::Float64
    R::Float64
    gamma::Float64
end

max_stars_repo_name: henrystoldt/fvCFD


text:
#Pkg.checkout("PlotRecipes")
using StatsBase
using Plots
using PlotRecipes
using LightGraphs
using StatPlots
#using GraphPlot
gr(size=(1000,1000))

#pyplot(size=(1000,1000))

function mainAlex(NN=20,TT=24000,pb=0.01,cb=1,cr=4)

    MessagesFired = Array{Float64}(24,NN)
    ResponsesTriggered = zeros(NN,NN,2)

    AA = dynamicCommunicators(NN,TT,pb,cb,cr,MessagesFired,ResponsesTriggered)
    degOut = totalOut(AA,NN)
    Cbroad = dynamicCentrality(AA,NN,TT)
#    scatterPlot(AA,size(AA)[1],degOut,Cbroad)
#    rankTimeSeries(MessagesFired,NN);
#    TimeSeries(MessagesFired,NN);

   # AA2 = ringNet()
   # degOut = totalOut(AA2,21)
   # Cbroad = dynamicCentrality(AA2,21,3)
   # scatterPlot(AA2,size(AA2)[1],degOut,Cbroad)
end

function scatterPlot(AA,NN,degOut,Cbroad)
    labels = [string(x) for x in 1:NN]
    scatter(degOut,Cbroad,markersize=1,series_annotations =labels,grid=false)
    xlabel!("total deg out")
    ylabel!("broadcast centrality")
    title!(string("node num=",NN))
end

function multiHistogram(data,file,title)
    labels = ["Pre-Polarization", "Post-Polarization"]
    groupedbar(data, labels=labels, nbins=20)
    xlabel!("Nodes")
    ylabel!("Broadcast Power Relative to Weakest Node")
    title!(title)
    savefig(file)
end

function graphVisualization(G,NN,file)
    names = [string(x) for x in 1:NN]
  #  g = DiGraph(G)
    graphplot(G,nodesize=22,names=names)
    savefig(file)
    println("Graph Visualization")

end

function heatMap(AA,file,flag)
        xs = [string(x) for x in 1:20]
	ys = [string(x) for x in 1:20]
	heatmap(xs,ys,AA);
	title!(string("Send vs Recieve ", (flag==1) ? "Non Polar" : "Polarized"))
	savefig(file)
	println("Heat Map")
end

function rankTimeSeries(MessagesFired,NN)
    Rank = Array{Float64}(16,NN)

    for ii in 1:16
    # Rank[ii,jj] = Num nodes lower than node ii,jj at timepoint (ii*1000)
        for jj in 1:NN
	    Rank[ii,jj] = sum([MessagesFired[ii,jj] < x for x in MessagesFired[ii,:]])
	end
    end
    plot(Rank,label=[string("Node ", x) for x in 1:20],grid=false,background_color=RGB(.2,.2,.2))
    xlabel!("Time")
    ylabel!("Rank")
    title!(string("Number of Nodes = 20 - Rank vs Time - Polar Importance"))
    savefig("Plots/Rank_Importance_Evolving_Polar");
end

function TimeSeries(Data,NN,file,title)

    plot(Data,label=["Power Top/Bot"],grid=false,background_color=RGB(.2,.2,.2))
    xlabel!("Time")
    ylabel!("Average Power Top Nodes over Bottom Nodes")
    title!(title)
    savefig(file);
end

function ringNet()
    AA2 = zeros(21,21,3)
    AA2[1,7,1] =1;AA2[7,1,1] =1;AA2[2,7,1] =1;AA2[7,2,1] =1;
    AA2[21,20,1] =1;AA2[20,21,1] =1;AA2[21,19,1] =1;AA2[19,21,1] =1;
    AA2[2,8,2]=1;AA2[8,2,2]=1;AA2[3,8,2]=1;AA2[8,3,2]=1;
    AA2[20,18,2]=1;AA2[18,20,2]=1;AA2[20,17,2]=1;AA2[17,20,2]=1;
    AA2[19,16,2]=1;AA2[16,19,2]=1;
    AA2[4,10,3]=1;AA2[10,4,3]=1;AA2[5,10,3]=1;AA2[10,5,3]=1;
    AA2[18,15,3]=1;AA2[15,18,3]=1;AA2[18,14,3]=1;AA2[14,18,3]=1;
    AA2[17,14,3]=1;AA2[14,17,3]=1;AA2[17,13,3]=1;AA2[13,17,3]=1;
    AA2[16,13,3]=1;AA2[13,16,3]=1;AA2[16,12,3]=1;AA2[12,16,3]=1;
    return AA2
end

function dynamicCommunicators(NN,TT,pb,cb,cr,MessagesFired,ResponsesTriggered)

    imp = [x^3 for x in 1:NN]
    Mask1 = [(-1)^x for x in 1:NN]
    Mask2 = [(-1)^(x+1) for x in 1:NN]
    Messages = Array{Float64}(20)
    data = Array{Float64}(20,2)
    Avg = Array{Float64}(24)
    d = 1
    f = 1
    AA = zeros(NN,NN,TT)
    i=1;

    #INIT DATA>>
    for ii in 1:NN
        if(rand(1)[1] <= pb)
            tmp = deleteat!(collect(1:NN),ii)
            destinationN = tmp[rand(1:end)]
            AA[ii,destinationN,1] = 1
        end

    end

    #DATA LOOP: 1st put basal edges into <present> / 2nd put response from past to now
    for tt in 2:(TT)
        #println(string("tt=",tt))

        # Collect Messages Fired to Calculate Rank
        if(tt%1000 == 0)
            MessagesFired[i,:] = Messages;
	    i+=1;
        end

        #BASAL LOOP
        for ii in 1:NN
            if(rand(1)[1] <= pb)
                tmp = deleteat!(collect(1:NN),ii)
                destinationN = tmp[rand(1:end)]
                AA[ii,destinationN,tt] = 1#FIRE NOW
            end

        end

        #RESPONSE LOOP
        for ii in 1:NN
            msgsTo_ii = AA[:,ii,tt-1]#RESPOND TO PREVIOUS STATE (tt-1)
            if(sum(AA[:,ii,tt-1]) >= 1)

	        if(tt <= 8000)
                    totalImportance = sum([msgsTo_ii[x]*imp[x] for x in 1:NN])
		else
		    if(ii%2==0)
                       totalImportance = sum([msgsTo_ii[x]*imp[x]*Mask1[x] for x in 1:NN])
                    else
                       totalImportance = sum([msgsTo_ii[x]*imp[x]*Mask2[x] for x in 1:NN])
	            end
		end

                r_nNumerator = totalImportance
                r_nDenominator = 1 + (findmax(imp)[1] * sum(AA[:,ii,tt-1]))
                r_next = r_nNumerator / r_nDenominator
              #  println(r_next)
                if(r_next >=  rand(1)[1])#if so generate cr links

                    tmp = deleteat!(collect(1:NN),ii)
                    destinationNodes = tmp[randperm(length(tmp))[1:cr]]
                    AA[ii,destinationNodes,tt] = 1#FIRE NOW
		    imp =5*(msgsTo_ii) + imp;
		    Messages = Messages + msgsTo_ii;
		    if(tt < 8000)
		        ResponsesTriggered[:,ii,1] += msgsTo_ii;
		    else
		        ResponsesTriggered[:,ii,2] += msgsTo_ii;
		    end
                end
            end
        end
    end

    for ll in 1:24
    	Avg[ll] = (sum(MessagesFired[ll,11:20])/10)/(sum(MessagesFired[ll,1:10])/10)
    end

    TimeSeries(Avg,NN,"Plots/TopAvg-BotAvg-PolarvsNonPolar","Number of Nodes = 20 - Power Top/Bottom vs Time")

  #  temp1 = sum(ResponsesTriggered[:,:,1],2)
  #  temp2 = sum(ResponsesTriggered[:,:,2],2)
  #  data[:,1] = temp1./temp1[1]
  #  data[:,2] = temp2./temp2[1]

  #  heatMap(ResponsesTriggered[:,:,1], "Plots/ResponsesTriggered_1-8000",1);
  #  heatMap(ResponsesTriggered[:,:,2], "Plots/ResponsesTriggered_8001-16000",2);
  #  multiHistogram(data[:,:,1], "Plots/ResponsesTriggered_GroupedBar_Lowest", "Responses Relative to Lowest - Polarized to Non-Polarized - Nodes = 20")

    return AA
end

function totalOut(AA,NN)
    coltmp = zeros(NN,1)
    outDegVec = sum(sum(AA,3),2)
    coltmp[:] = outDegVec[:]
    return coltmp
end

function dynamicCentrality(AA,NN,TT)
    println(":")
    daynorm = zeros(TT,1);
    for ii in 1:TT
        Atemp = AA[:,:,ii]
        daynorm[ii] = findmax(abs(eigvals(Atemp)))[1]
    end
    alpha = 0.9*(1/findmax(daynorm)[1]);#90percent of the max
    println(alpha)
    #alpha = 0.15;
    BB = eye(NN,NN);
    for kk = TT:-1:1
        Aday = AA[:,:,kk]
        #println(Aday)
        BBtmp = inv(eye(NN,NN) - alpha*Aday)#\BB
        #BBtmp = (eye(NN,NN) - alpha*Aday) \ BB
        BBtmp = abs.(BBtmp)
        #BB = BBtmp / norm(BBtmp)
        BB = (BB*BBtmp) / norm(BB)
    end
    #sum across columns
    Cbroad = sum(BB',2)
    return Cbroad
end

function dynamicCentrality2(AA,NN,TT)
    println(":")
    daynorm = zeros(TT,1);
    for ii in 1:TT
        Atemp = AA[:,:,ii]
        daynorm[ii] = findmax(abs(eigvals(Atemp)))[1]
    end
    alpha = 0.9*(1/findmax(daynorm)[1]);#90percent of the max
    println(alpha)
    #alpha = 0.15;
    BB = eye(NN,NN);
    for kk = TT:-1:1
        Aday = AA[:,:,kk]
        #println(Aday)
        BBtmp = inv(eye(NN,NN) - alpha*Aday)#\BB
        #BBtmp = (eye(NN,NN) - alpha*Aday) \ BB
        BBtmp = abs.(BBtmp)
        #BB = BBtmp / norm(BBtmp)
        BB = (BB*BBtmp) / norm(BB)
    end
   return BB
end

max_stars_repo_name: mantzaris/competitiveCommunicators

max_issues_repo_name: mantzaris/competitiveCommunicators

max_forks_repo_name: mantzaris/competitiveCommunicators