text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -51.708688423784224 <= q <= 53.49505517188359)

                     add_NL_constraint(m, :(sigmoid(0.12651985253467002 + 0.8800990177661969*sigmoid(0.9367289958343452 + 0.010085610540752743*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.46191185446397043*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.4430556996324988*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.9244167282970754*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.15008923933940022*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.9059509477084742*sigmoid(-0.9076856947519341 + 0.6427398296015978*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.5974696839876872*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.3127996929073351*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.8921846648994118*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.8861141904898706*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.7662733112940292*sigmoid(0.014219790982957559 + -0.1830460372609095*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.18937048380167898*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.48856019606524415*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.6243584215460145*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.6378733350620909*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.3066544631255512*sigmoid(0.37873800689777815 + 0.3786798527526116*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.15214183164812667*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.4531653277321981*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + -0.2698943947473773*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.8205030207972928*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.673637701918695*sigmoid(0.41028851710489267 + 0.3789873528571355*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.13395473418456794*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.9211828416113521*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.005265087561311166*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.6679558602332873*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4])))) + sigmoid(-0.8950204772394628 + -0.19603919659654778*sigmoid(0.9367289958343452 + 0.010085610540752743*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.46191185446397043*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.4430556996324988*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.9244167282970754*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.15008923933940022*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.902347490498435*sigmoid(-0.9076856947519341 + 0.6427398296015978*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.5974696839876872*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.3127996929073351*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.8921846648994118*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.8861141904898706*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.30440629677970144*sigmoid(0.014219790982957559 + -0.1830460372609095*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.18937048380167898*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.48856019606524415*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.6243584215460145*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.6378733350620909*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.6770808003403048*sigmoid(0.37873800689777815 + 0.3786798527526116*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.15214183164812667*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.4531653277321981*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + -0.2698943947473773*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.8205030207972928*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.03232239144278681*sigmoid(0.41028851710489267 + 0.3789873528571355*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.13395473418456794*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.9211828416113521*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.005265087561311166*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.6679558602332873*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4])))) + sigmoid(-0.11076390433121786 + 0.7661103402200116*sigmoid(0.9367289958343452 + 0.010085610540752743*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.46191185446397043*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.4430556996324988*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.9244167282970754*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.15008923933940022*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.9009199905889229*sigmoid(-0.9076856947519341 + 0.6427398296015978*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.5974696839876872*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.3127996929073351*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.8921846648994118*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.8861141904898706*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.5631507797475948*sigmoid(0.014219790982957559 + -0.1830460372609095*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.18937048380167898*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.48856019606524415*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.6243584215460145*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.6378733350620909*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.4163802251888926*sigmoid(0.37873800689777815 + 0.3786798527526116*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.15214183164812667*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.4531653277321981*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + -0.2698943947473773*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.8205030207972928*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.17662226846197893*sigmoid(0.41028851710489267 + 0.3789873528571355*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.13395473418456794*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.9211828416113521*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.005265087561311166*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.6679558602332873*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4])))) + sigmoid(0.8242960442017702 + -0.7506191216499181*sigmoid(0.9367289958343452 + 0.010085610540752743*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.46191185446397043*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.4430556996324988*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.9244167282970754*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.15008923933940022*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.31521888864028513*sigmoid(-0.9076856947519341 + 0.6427398296015978*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.5974696839876872*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.3127996929073351*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.8921846648994118*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.8861141904898706*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.5308637261158768*sigmoid(0.014219790982957559 + -0.1830460372609095*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.18937048380167898*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.48856019606524415*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.6243584215460145*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.6378733350620909*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.6654841763336146*sigmoid(0.37873800689777815 + 0.3786798527526116*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.15214183164812667*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.4531653277321981*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + -0.2698943947473773*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.8205030207972928*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.2645412007156045*sigmoid(0.41028851710489267 + 0.3789873528571355*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.13395473418456794*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.9211828416113521*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.005265087561311166*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.6679558602332873*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4])))) + sigmoid(-0.6657453597374303 + 0.46323784626870257*sigmoid(0.9367289958343452 + 0.010085610540752743*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.46191185446397043*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.4430556996324988*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.9244167282970754*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.15008923933940022*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.7128966868012561*sigmoid(-0.9076856947519341 + 0.6427398296015978*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.5974696839876872*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.3127996929073351*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.8921846648994118*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.8861141904898706*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.43240821263487783*sigmoid(0.014219790982957559 + -0.1830460372609095*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.18937048380167898*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.48856019606524415*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.6243584215460145*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.6378733350620909*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.6998749584965918*sigmoid(0.37873800689777815 + 0.3786798527526116*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.15214183164812667*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.4531653277321981*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + -0.2698943947473773*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.8205030207972928*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.6177298271304656*sigmoid(0.41028851710489267 + 0.3789873528571355*sigmoid(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.13395473418456794*sigmoid(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.9211828416113521*sigmoid(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.005265087561311166*sigmoid(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.6679558602332873*sigmoid(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m


text:
module OptimNewton

export optimize_newton

include("newton.jl")

end # module

ext: jl


text:
module DSToolsDemo

# # DSTOOLSDEMO    Demonstration of the Descriptor System Tools (DSTOOLS).
# #                A. Varga.

using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

function pause()
    # Comment out the following two lines for a non-interactive execution
    println("Press any key to continue")
    readline()
end

s = rtf('s'); z = rtf('z'); # define the complex variables s and z
println("s = rtf('s'); z = rtf('z'); # define the complex variables s and z ")
pause()

println("Gc = [s^2 s/(s+1); 0 1/s]   # define the 2-by-2 improper Gc(s)")
Gc = [s^2 s/(s+1); 0 1/s]   # define the 2-by-2 improper Gc(s)
Gd = [z^2 z/(z-2); 0 1/z]   # define the 2-by-2 improper Gd(z)
println(Gc)
pause()

println("Gd = [z^2 z/(z-2); 0 1/z]   # define the 2-by-2 improper Gd(z)")
println(Gd)
pause()

println("Build LTI descriptor realizations of Gc(s) and Gd(z) ")
println("sysc = dss(Gc, minimal = true);  # build continuous-time descriptor system realization")
sysc = dss(Gc, minimal = true);  # build continuous-time descriptor system realization
println(sysc)
pause()

println("dss(Gd, minimal = true)  # build discrete-time descriptor system realization")
sysd = dss(Gd, minimal = true);   # build continuous-time descriptor system realization
println(sysd)
pause()

println("Compute poles and zeros of Gc(s)")
println("gpole(sysc)   #  computes all poles (finite and infinite)")
display(gpole(sysc))
pause()

println("gzero(sysc)   #  computes all zeros (finite and infinite)")
display(gzero(sysc))
pause()

println("Descriptor systems (A-λE,B,C,D) with singular pole pencil A-λE are not supported in DSTOOLS!")
println("Non-regular A-λE can be detected in several ways.\n ")
println("(1) non-regular A-λE can be detected when building a descriptor system ")
println("A = zeros(1,1); E = zeros(1,1); B = ones(1,1); C = ones(1,1); D = zeros(1,1);")
println("dss(A,E,B,C,D,check_reg = true)")
A = zeros(1,1); E = zeros(1,1); B = ones(1,1); C = ones(1,1); D = zeros(1,1);
try
  dss(A,E,B,C,D,check_reg = true)  # check regularity when building the descriptor system
catch err
  println(err)
end
pause()

println("(2) non-regular A-λE can be explicitly detected using isregular")
println("syst = dss(A,E,B,C,D); # no regularity check performed when building a descriptor system")
println("isregular(syst) == false")
syst = dss(A,E,B,C,D);
println(@test isregular(syst) == false)
pause()

println("(3) non-regular A-λE can be detected using isproper")
println("isproper(syst) == false")
syst = dss(A,E,B,C,D);
println(@test isproper(syst) == false)
pause()

println("(4) gpole can also be used to check regularity of the pole pencil")
println("gpole(syst)    # the system has some poles equal to NaN")
display(gpole(syst))
pause()

println("Let's try some decomposition and factorization functions\n")

println("Use gsdec to compute the separation of proper and polynomial parts of Gc(s)")
println("sysp, sysi = gsdec(sysc); # Gc(s) = Gcp(s) + Gci(s)")
sysp, sysi = gsdec(sysc);
println("@test iszero(sysc-sysp-sysi)  # checking the decomposition")
println(@test iszero(sysc-sysp-sysi))
pause()
println("dss2rm(sysp)  # display proper part Gcp(s)")
println(dss2rm(sysp))
pause()
println("dss2rm(sysi)  # display polynomial part Gci(s)")
println(dss2rm(sysi))
pause()

println("Compute the stable and proper right coprime factorization of Gc(s) as")
println("Gc(s) = N(s)*inv(M(s)), with N(s) and M(s) stable and proper,")
println(" having a stability degree -1 and poles assigned to [-2,-3,-4]")
println("sysn, sysm = grcf(sysc,evals = [-2,-3,-4],sdeg = -1,mindeg = true, mininf = true);")
sysn, sysm = grcf(sysc,evals = [-2,-3,-4], sdeg = -1, mindeg = true,  mininf = true);
println("@test iszero(sysc*sysm-sysn)  # checking the factorization")
println(@test iszero(sysc*sysm-sysn))
pause()

println("@test isproper(sysm) && isproper(sysn) # checking properness of factors")
println(@test isproper(sysm) && isproper(sysn) )
pause()

println("@test (maximum(real(gpole(sysm))) < 0 && maximum(real(gpole(sysn))) < 0) # checking stability of poles")
println(@test (maximum(real(gpole(sysm))) < 0 && maximum(real(gpole(sysn))) < 0) )
pause()

println("@test isempty(gzero(gminreal([sysn;sysm]))) # checking coprimeness:  [N(s);M(s)] has no zeros")
println(@test isempty(gzero(gminreal([sysn;sysm]))))
pause()

println("Compute the right coprime factorization with inner denominator of Gd(z) as")
println("Gd(z) = N(z)*inv(M(z)), with N(z) and M(z) stable and proper, and M(z) inner")
println("sysni, sysmi = grcfid(sysd);")
sysni, sysmi = grcfid(sysd);
println("@test iszero(sysc*sysm-sysn)  # checking the factorization")
println(@test iszero(sysc*sysm-sysn))
pause()

println("@test maximum(abs.(gpole(sysmi))) < 1 && maximum(abs.(gpole(sysni))) < 1) # checking stability of poles")
println(@test maximum(abs.(gpole(sysmi))) < 1 && maximum(abs.(gpole(sysni))) < 1  )
pause()

println("@test iszero(sysmi'*sysmi-I) # checking the innerness of M(z): conj(M(z))*M(z)-I = 0")
println(@test iszero(sysmi'*sysmi-I))
pause()

println("@test isempty(gzero(gminreal([sysn;sysm]))) # checking coprimeness:  [N(z);M(z)] has no zeros")
println(@test isempty(gzero(gminreal([sysn;sysm]))))
pause()

println("Compute the inner-outer factorization of stable 3x3 G(s) with rank 2 ")
println("s = Polynomial([0, 1],'s'); ")
s = Polynomial([0, 1],'s');
println("num = [(s-1) s 1; 0 (s-2) (s-2); (s-1) (s^2+2*s-2) (2*s-1)]; # numerators")
num = [(s-1) s 1; 0 (s-2) (s-2); (s-1) (s^2+2*s-2) (2*s-1)]; # numerators
pause()
println("den = [(s+2) (s+2) (s+2); 1 (s+1)^2 (s+1)^2; (s+2) (s+1)*(s+2) (s+1)*(s+2)]; # denominators")
den = [(s+2) (s+2) (s+2); 1 (s+1)^2 (s+1)^2; (s+2) (s+1)*(s+2) (s+1)*(s+2)]; # denominators
pause()
println("sys = dss(num,den,minimal = true, atol = 1.e-7);  # build a minimal descriptor realization")
sys = dss(num,den,minimal = true, atol = 1.e-7);
# println(sys)
pause()
println("Analysis of some properties")
println("gpole(sys);  # the system is stable")
display(gpole(sys))    # the system is stable
pause()
println("gzero(sys,atol1=1.e-7)    # the system has 2 unstable zeros and an infinite zero")
display(gzero(sys,atol1=1.e-7) )   # the system has 2 unstable zeros and an infinite zero
pause()
println("nr = gnrank(sys,atol1=1.e-7)    # the normal rank of G(s) is 2")
nr = gnrank(sys,atol1=1.e-7);  # the normal rank of G(s) is 2
println(nr)
pause()

println("Compute the inner-quasi-outer factorization of G(s) as G(s) = Gi(s)*[Go(s);0],")
println("with Gi(s) inner and square, and Go(s) quasi-outer (i.e., full row rank, ")
println("and without zeros in the open right-half plane), can be computed by using")
println("sysi, syso, info = giofac(sys)")
sysi, syso, info = giofac(sys);
pause()

println("Checking the factorization: Gi(:,1:nr)(s)*Go(s) = G(s)")
println("@test iszero(sysi[:,1:nr]*syso-sys,atol1=1.e-7)")
println(@test iszero(sysi[:,1:nr]*syso-sys,atol1=1.e-7))
pause()

println("Checking the innerness of Gi(s)")
println("@test iszero(sysi'*sysi-I,atol1=1.e-7)")
println(@test iszero(sysi'*sysi-I,atol1=1.e-7))
pause()

println("If G(s) contains a so-called free inner factor, then this factor")
println("is included in Gi(s), but the realization of Go(s) is not minimal")
println("@test order(syso)-order(gir(syso)) > 0  # a free inner factor is present in G(s)")
println(@test order(syso)-order(gir(syso)) > 0 ) # a free inner factor is present in G(s)
pause()

println("Checking that Go(s) has no zeros in the open right-half plane")
println("zer = gzero(syso,atol1 = 1.e-7);  # compute zeros of Go")
zer = gzero(syso,atol1 = 1.e-7);
println("@test all(real.(zer[.!isinf.(zer)]) .< 0 )")
println(@test all(real.(zer[.!isinf.(zer)]) .< 0 ))
pause()

##
println("Solution of linear rational equations G(s)*X(s) = F(s)")
println("Consider the Wang and Davison example (IEEE Trans. Autom. Contr.,1973)")
println("to determine a right inverse of G(s) by solving G(s)*X(s) = I.")
println("s = Polynomial([0, 1],'s'); # define polynomial s")
s = Polynomial([0, 1],'s');
println("num =  [ s+1 s+3 s^2+3*s; s+2 s^2+2*s 0 ]; # numerators")
num =  [ s+1 s+3 s^2+3*s; s+2 s^2+2*s 0 ]; # numerators
# println(num)
pause()
println("den = [s^2+3*s+2 s^2+3*s+2 s^2+3*s+2; s^2+3*s+2 s^2+3*s+2 s^2+3*s+2]; # denominators")
den = [s^2+3*s+2 s^2+3*s+2 s^2+3*s+2; s^2+3*s+2 s^2+3*s+2 s^2+3*s+2]; # denominators
# println(den)
pause()
println("sysg = dss(num,den,minimal = true, atol = 1.e-7);  # build a minimal descriptor realization of G(s)")
sysg = dss(num,den,minimal = true, atol = 1.e-7);
# println(sysg)
pause()
println("sysf = dss([1 0;0 1.]);  # build a minimal descriptor realization of F(s) = I")
sysf = dss([1 0;0 1.]);
# println(sysf)
pause()

println("zer = gzero(sysg)   # G(s) has no zeros, thus stable right inverses exist")
zer = gzero(sysg)         # G(s) has no zeros, thus stable right inverses exist
println("zer = $zer")
pause()

println("Compute a stable right inverse with poles in [-1 -2 -3] using")
println("sysx, info = grsol(sysg,sysf,poles = [-1, -2, -3]);")
sysx, info = grsol(sysg, sysf, poles = [-1, -2, -3]);
# println(sysx)
pause()

println("Checking the solution: G(s)*X(s) - I = 0")
println("@test iszero(sysg*sysx-I,atol=1.e-7)")
println(@test iszero(sysg*sysx-I,atol=1.e-7))   #  G(s)*X(s) - I = 0
pause()

println("Check assigned poles: gpole(sysx) = [-1, -2, -3]")
println(@test sort(gpole(sysx)) ≈ sort([-1, -2, -3]))
pause()

println("We can also compute a right inverse of least order = 2")
println("sysxmin, = grsol(sysg,sysf,mindeg = true); order(sysxmin) == 2")
sysxmin, = grsol(sysg,sysf,mindeg = true);
println(@test iszero(sysg*sysxmin-I,atol=1.e-7) && order(sysxmin) == 2)
pause() # Press any key to continue ...

gpole(sysxmin)
println("The least order inverse is unstable: max(real(gpole(sysxmin))) > 0")
println(@test maximum(real(gpole(sysxmin))) > 0)
pause()

##
println("\nWe can try to explicitly determine a least order right inverse")
println("using stable generators (X0(s),XN(s)), where X0(s) is a particular")
println("stable right inverse satisfying G(s)*X0(s) = I and XN(s) is a stable")
println("basis of the right nullspace of G(s), satisfying G(s)*XN(s) = 0;")
println("all solutions of G(s)*X(s) = I are given by X(s) = X0(s)+XN(s)*Y(s),")
println("where Y(s) is arbitrary\n")

println("Determine the generators X0 and XN with poles assigned to [-1 -2 -3]")
println("_,_,sysgen = grsol(sysg,sysf,poles = [-1, -2, -3],solgen = true);")
_,_,sysgen = grsol(sysg,sysf,poles = [-1, -2, -3],solgen = true);
pause()

println("Compute a least order solution X2(s) = X0(s)+XN(s)*Y2(s) using")
println("order reduction based on minimal dynamic covers of Type 2")

println("sysx2,sysy2, = grmcover2(sysgen,2); order(gminreal(sysx2)) == 2")
sysx2,sysy2, = grmcover2(sysgen,2);
println(@test order(gminreal(sysx2)) == 2)
pause()

println("Checking the solution: G(s)*X2(s) = I")
println("@test iszero(sysg*sysx2-I,atol=1.e-7)" )
println(@test iszero(sysg*sysx2-I,atol=1.e-7) )  #  G(s)*X2(s) - I = 0
pause()

println("The least order inverse is unstable: max(real(gpole(sysx2))) > 0")
println(@test maximum(real(gpole(sysx2))) > 0)
pause()

println("Checking the minimal cover reduction results: X2(s) = X0(s)+XN(s)*Y2(s)")
println("@test iszero(sysx2-sysgen[:,1:2]-sysgen[:,3]*sysy2,atol=1.e-7)")
println(@test iszero(sysx2-sysgen[:,1:2]-sysgen[:,3]*sysy2,atol=1.e-7) )
pause()

println("\nSolution of an H∞ model matching problem  min||X(s)*G(s)-F(s)||")
println("\nSetting up an example taken from Francis' book (1987)")
println("s = Polynomial([0, 1],'s');")
println("W = (s+1)/(10s+1); # weighting function")
println("G = dss([ -(s-1)/(s^2+s+1)*W; (s^2-2*s)/(s^2+s+1)*W], minimal=true);")
println("F = dss([ W; 0 ]);")

s = Polynomial([0, 1],'s');
W = (s+1)/(10s+1); # weighting function
G = dss([ -(s-1)/(s^2+s+1)*W; (s^2-2*s)/(s^2+s+1)*W], minimal=true);
F = dss([ W; 0 ]);
pause()

println("\nWe employ the γ-iteration based solution proposed by Francis (1987)")

println("Xopt, info = grasol(G, F, mindeg = true, atol = 1.e-7);")
println("This will take some time at first execution!!!")
Xopt, info = grasol(G, F, mindeg = true, atol = 1.e-7);
println("Optimum of approximation error: info.mindist = $(info.mindist)")
pause()

println("This fully agrees with info.mindist: ghinfnorm(G*Xopt-F) ≈ info.mindist")
println(@test ghinfnorm(G*Xopt-F)[1] ≈ info.mindist)
pause()

println("\nThe solution in Francis' book corresponds to a suboptimal solution for γ = 0.2729")
println("Xsub, info = grasol(G, F, 0.2729; mindeg = true, atol = 1.e-7);")
Xsub, info = grasol(G, F, 0.2729; mindeg = true, atol = 1.e-7);
println("The suboptimal approximation error: info.mindist = $(info.mindist)")
pause()

println("This fully agrees with info.mindist: ghinfnorm(G*Xsub-F) ≈ info.mindist")
println(@test ghinfnorm(G*Xsub-F)[1] ≈ info.mindist)
pause()

println("\nLet's illustrate the use of the bilinear transformation")

println("s = rtf('s');")
println("G  = [s^2 s/(s+1); 0 1/s]")
println("sys = dss(G, minimal = true);")
s = rtf('s');                 # define the complex variable s
G  = [s^2 s/(s+1); 0 1/s]     # define the 2-by-2 improper G(s)
sys = dss(G, minimal = true); # build continuous-time descriptor system realization
pause()

println("Pole-zero analysis shows that G(s) has poles and zeros in the origin and at infinity")
println("pol = gpole(sys);")
println("zer = gzero(sys);")
pol = gpole(sys);
zer = gzero(sys);
println("pol = ")
display(pol)
println("zer = ")
display(zer)
pause()

println("Define a bilinear transformation with g(s) = (s+0.01)/(1+0.01*s) to make")
println("all poles and zeros stable and finite using")
println("g = (s+0.01)/(1+0.01*s);")
g = (s+0.01)/(1+0.01*s);
println("\ng = \n")
print(g)
pause()

println("Compute the transformed system and resulting stable poles and zeros")
println("syst = gbilin(sys,g, atol = 1.e-7,minimal = true)[1];")
syst = gbilin(sys,g, atol = 1.e-7,minimal = true)[1];
pol_new = gpole(syst)
zer_new = gzero(syst)
println("pol_new = ")
display(pol_new)
println("zer_new = ")
display(zer_new)
pause()

println("Compute the ν-gap distance between the two models using Vinnicombe's formula")
println("nugap = gnugap(sys,syst)[1]")
nugap = gnugap(sys,syst)[1]
println("The resulting ν-gap distance $nugap reflects well to the s-variable perturbation")

end # module


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    g1 ∩ g2

Return the intersection of two geometries `g1` and `g2`.
"""
Base.intersect(g1::Geometry, g2::Geometry) = get(intersecttype(g1, g2))

"""
    intersecttype(g1, g2)

Return the intersection type of two geometries `g1` and `g2`.
"""
function intersecttype end

"""
    Intersection

An intersection type (e.g. crossing line segments, overlapping boxes).
"""
abstract type Intersection end

Base.get(I::Intersection) = I.value

# ------------------------
# LINE-LINE INTERSECTIONS
# ------------------------

struct CrossingLines{P<:Point} <: Intersection
  value::P
end

struct OverlappingLines{L<:Line} <: Intersection
  value::L
end

# ------------------------------
# SEGMENT-SEGMENT INTERSECTIONS
# ------------------------------

struct CrossingSegments{P<:Point} <: Intersection
  value::P
end

struct MidTouchingSegments{P<:Point} <: Intersection
  value::P
end

struct CornerTouchingSegments{P<:Point} <: Intersection
  value::P
end

struct OverlappingSegments{S<:Segment} <: Intersection
  value::S
end

# ----------------------
# BOX-BOX INTERSECTIONS
# ----------------------

struct OverlappingBoxes{B<:Box} <: Intersection
  value::B
end

struct FaceTouchingBoxes{B<:Box} <: Intersection
  value::B
end

struct CornerTouchingBoxes{P<:Point} <: Intersection
  value::P
end

# -------------------------------
# SEGMENT-TRIANGLE INTERSECTIONS
# -------------------------------

struct IntersectingSegmentTriangle{P<:Point} <: Intersection
  value::P
end

# ----------------------------
# SEGMENT-PLANE INTERSECTIONS
# ----------------------------

struct CrossingSegmentPlane{P<:Point} <: Intersection
  value::P
end

struct TouchingSegmentPlane{P<:Point} <: Intersection
  value::P
end

struct OverlappingSegmentPlane{S<:Segment} <: Intersection
  value::S
end

# ------------
# CORNER CASE
# ------------

struct NoIntersection <: Intersection end

Base.get(::NoIntersection) = nothing

# ----------------
# IMPLEMENTATIONS
# ----------------

include("intersections/lines.jl")
include("intersections/segments.jl")
include("intersections/planes.jl")
include("intersections/triangles.jl")
include("intersections/boxes.jl")

"""
    hasintersect(g1, g2)

Return `true` if geometries `g1` and `g2` intersect and `false` otherwise.

The algorithm works with any geometry that has a well-defined [`supportfun`](@ref).

## References

* Gilbert, E., Johnson, D., Keerthi, S. 1988. [A fast
  Procedure for Computing the Distance Between Complex
  Objects in Three-Dimensional Space]
  (https://ieeexplore.ieee.org/document/2083)
"""
function hasintersect(g1::Geometry{Dim,T}, g2::Geometry{Dim,T}) where {Dim,T}
  # handle non-convex geometries
  if !isconvex(g1)
    d1 = discretize(g1, FIST())
    return hasintersect(d1, g2)
  elseif !isconvex(g2)
    d2 = discretize(g2, FIST())
    return hasintersect(g1, d2)
  end

  # initial direction
  c1, c2 = centroid(g1), centroid(g2)
  d = c1 == c2 ? rand(Vec{Dim,T}) : c2 - c1

  # first point in Minkowski difference
  P = minkowskipoint(g1, g2, d)

  # origin of coordinate system
  O = minkowskiorigin(Dim, T)

  # initialize simplex vertices
  points = [P]

  # move towards the origin
  d = O - P
  while true
    P = minkowskipoint(g1, g2, d)
    if (P - O) ⋅ d < zero(T)
      return false
    end
    push!(points, P)

    # line segment case
    if length(points) == 2
      B, A = points
      AB = B - A
      AO = O - A
      d  = AB × AO × AB
    else # simplex case
      C, B, A = points
      AB  = B - A
      AC  = C - A
      AO  = O - A
      ABᵀ = AC × AB × AB
      ACᵀ = AB × AC × AC
      if ABᵀ ⋅ AO > zero(T)
        popat!(points, 1) # pop C
        d = ABᵀ
      elseif ACᵀ ⋅ AO > zero(T)
        popat!(points, 2) # pop B
        d = ACᵀ
      else
        return true
      end
    end
  end
end

hasintersect(d1::Domain, g2::Geometry) =
  any(g1 -> hasintersect(g1, g2), d1)

hasintersect(g1::Geometry, d2::Domain) =
  hasintersect(d2, g1)

function hasintersect(d1::Domain, d2::Domain)
  for g1 in d1, g2 in d2
    hasintersect(g1, g2) && return true
  end
  return false
end

# support point in Minkowski difference
function minkowskipoint(g1::Geometry{Dim,T}, g2::Geometry{Dim,T}, d) where {Dim,T}
  n = Vec{Dim}(d[1:Dim])
  v = supportfun(g1, n) - supportfun(g2, -n)
  Point(ntuple(i-> i ≤ Dim ? v[i] : zero(T), max(Dim, 3)))
end

# origin of coordinate system
minkowskiorigin(Dim, T) = Point(ntuple(i->zero(T), max(Dim, 3)))


text:
# testing ERB filterbank design
using Compat.DelimitedFiles
fb = make_erb_filterbank(16000, 23, 150)
fb_matlab = readdlm(open(joinpath(dirname(@__FILE__), "data", "ERB_filter_coeffs.csv")), ',')
@test length(fb.filters) == size(fb_matlab, 1)
for k=1:length(fb.filters)
	sos = fb.filters[k]
	biquads = sos.biquads
	g = sos.g
	for b=1:4
		@test biquads[b].b0 ≈ fb_matlab[k,1]
		@test biquads[b].b2 ≈ fb_matlab[k,6]
		@test biquads[b].a1 ≈ fb_matlab[k,8]
		@test biquads[b].a2 ≈ fb_matlab[k,9]
	end
	@test biquads[1].b1 ≈ fb_matlab[k,2]
	@test biquads[2].b1 ≈ fb_matlab[k,3]
	@test biquads[3].b1 ≈ fb_matlab[k,4]
	@test biquads[4].b1 ≈ fb_matlab[k,5]
	@test 1/g ≈ fb_matlab[k,10]
end

# testing ERB filterbank response
x = zeros(10000)
x[1] = 1.0
y = filt(fb, x)
y_matlab = readdlm(open(joinpath(dirname(@__FILE__), "data", "ERB_filter_response.csv")), ',')
@test y ≈ y_matlab'


text:
using Documenter, Includeonce

makedocs(
    modules = [Includeonce],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Jacques David",
    sitename = "Includeonce.jl",
    pages = Any["index.md"]
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

deploydocs(
    repo = "github.com/jdadavid/Includeonce.jl.git",
    push_preview = true
)


text:
module StatPlots

using Reexport
@reexport using Plots
import Plots: _cycle
using Plots.PlotMeasures
using StatsBase
using Distributions
import IterableTables
import DataValues: DataValue
import TableTraits: column_types, column_names, getiterator, isiterabletable
import TableTraitsUtils: create_columns_from_iterabletable

import KernelDensity
import Loess
@recipe f(k::KernelDensity.UnivariateKDE) = k.x, k.density
@recipe f(k::KernelDensity.BivariateKDE) = k.x, k.y, k.density.'

@shorthands cdensity

export @df

include("df.jl")
include("corrplot.jl")
include("cornerplot.jl")
include("distributions.jl")
include("boxplot.jl")
include("violin.jl")
include("hist.jl")
include("marginalhist.jl")
include("bar.jl")

end # module


text:
using SossMLJExamples
using Test

@testset "SossMLJExamples.jl" begin
    # Write your tests here.
end


text:
#=
# CG solver
=#
module CGSolver

export solve, nonlinear_solve

import ..Femshop: JULIA, CPP, MATLAB, DENDRO, HOMG, CUSTOM_GEN_TARGET,
        SQUARE, IRREGULAR, UNIFORM_GRID, TREE, UNSTRUCTURED,
        CG, DG, HDG, FV,
        NODAL, MODAL, CELL, LEGENDRE, UNIFORM, GAUSS, LOBATTO,
        NONLINEAR_NEWTON, NONLINEAR_SOMETHING,
        EULER_EXPLICIT, EULER_IMPLICIT, CRANK_NICHOLSON, RK4, LSRK4, ABM4,
        DEFAULT_SOLVER, PETSC,
        VTK, RAW_OUTPUT, CUSTOM_OUTPUT,
        DIRICHLET, NEUMANN, ROBIN, NO_BC, FLUX,
        MSH_V2, MSH_V4,
        SCALAR, VECTOR, TENSOR, SYM_TENSOR, VAR_ARRAY,
        LHS, RHS,
        LINEMESH, QUADMESH, HEXMESH

import ..Femshop: log_entry, printerr
import ..Femshop: config, prob, variables, mesh_data, grid_data, refel, time_stepper, elemental_order, genfunctions
import ..Femshop: Variable, Coefficient, GenFunction
import ..Femshop: GeometricFactors, geo_factors, geometric_factors, build_deriv_matrix

using LinearAlgebra, SparseArrays

include("fe_boundary.jl");
include("nonlinear.jl")
include("cg_matrixfree.jl");
include("level_benchmark.jl");

function linear_solve(var, bilinear, linear, stepper=nothing; assemble_func=nothing)
    if config.linalg_matrixfree
        return solve_matrix_free_sym(var, bilinear, linear, stepper);
        #return solve_matrix_free_asym(var, bilinear, linear, stepper);
    end

    # If more than one variable
    if typeof(var) <: Array
        # multiple variables being solved for simultaneously
        dofs_per_node = 0;
        dofs_per_loop = 0;
        for vi=1:length(var)
            dofs_per_loop += length(var[vi].symvar);
            dofs_per_node += var[vi].total_components;
        end
    else
        # one variable
        dofs_per_loop = length(var.symvar);
        dofs_per_node = var.total_components;
    end
    N1 = size(grid_data.allnodes,2);
    Nn = dofs_per_node * N1;
    Np = refel.Np;
    nel = mesh_data.nel;

    # Allocate arrays that will be used by assemble
    # These vectors will hold the integrated values(one per cell).
    # They will later be combined.
    rhsvec = zeros(Nn);
    lhsmatI = zeros(Int, nel*dofs_per_node*Np*dofs_per_node*Np);
    lhsmatJ = zeros(Int, nel*dofs_per_node*Np*dofs_per_node*Np);
    lhsmatV = zeros(nel*dofs_per_node*Np*dofs_per_node*Np);
    allocated_vecs = [rhsvec, lhsmatI, lhsmatJ, lhsmatV];

    if prob.time_dependent && !(stepper === nothing)
        # First assemble both lhs and rhs
        assemble_t = @elapsed((A, b) = assemble(var, bilinear, linear, allocated_vecs, dofs_per_node, dofs_per_loop, 0, stepper.dt, assemble_loops=assemble_func));

        log_entry("Initial assembly took "*string(assemble_t)*" seconds");

        log_entry("Beginning "*string(stepper.Nsteps)*" time steps.");
        t = 0;
        sol = [];
        start_t = Base.Libc.time();
        assemble_t = 0;
        linsolve_t = 0;
        last2update = 0;
        last10update = 0;

        # allocate any temporary storage needed by steppers
        if stepper.type == LSRK4
            tmppi = zeros(size(b));
            tmpki = zeros(size(b));
        elseif stepper.stages > 1
            last_result = zeros(size(b));
            tmpresult = zeros(size(b));
            tmpki = zeros(length(b), stepper.stages);
        end
        print("Time stepping progress(%): 0");
        for i=1:stepper.Nsteps
            if stepper.stages > 1
                # LSRK4 is a special case, low storage
                if stepper.type == LSRK4
                    # Low storage RK4:
                    # p0 = u
                    #   ki = ai*k(i-1) + dt*f(p(i-1), t+ci*dt)
                    #   pi = p(i-1) + bi*ki
                    # u = p5

                    tmppi = get_var_vals(var, tmppi);
                    for rki=1:stepper.stages
                        rktime = t + stepper.c[rki]*stepper.dt;
                        # p(i-1) is currently in u
                        assemble_t += @elapsed(b = assemble(var, nothing, linear, allocated_vecs, dofs_per_node, dofs_per_loop, rktime, stepper.dt; rhs_only = true, assemble_loops=assemble_func));

                        linsolve_t += @elapsed(sol = A\b);
                        if rki == 1 # because a1 == 0
                            tmpki = stepper.dt .* sol;
                        else
                            tmpki = stepper.a[rki].*tmpki + stepper.dt.*sol;
                        end
                        tmppi = tmppi + stepper.b[rki].*tmpki

                        place_sol_in_vars(var, tmppi, stepper);
                    end

                else
                    # Explicit multi-stage methods:
                    # x = x + dt*sum(bi*ki)
                    # ki = rhs(t+ci*dt, x+dt*sum(aij*kj)))   j < i

                    # will hold the final result
                    last_result = get_var_vals(var, last_result);
                    # will be placed in var.values for each stage
                    # tmpvals = tmpresult;
                    # Storage for each stage
                    # tmpki = zeros(length(b), stepper.stages);
                    for stage=1:stepper.stages
                        stime = t + stepper.c[stage]*stepper.dt;

                        assemble_t += @elapsed(b = assemble(var, nothing, linear, allocated_vecs, dofs_per_node, dofs_per_loop, stime, stepper.dt; rhs_only = true, assemble_loops=assemble_func));

                        linsolve_t += @elapsed(tmpki[:,stage] = A\b);

                        tmpresult = get_var_vals(var, tmpresult);
                        for j=1:(stage-1)
                            if stepper.a[stage, j] > 0
                                tmpresult += stepper.dt * stepper.a[stage, j] .* tmpki[:,j];
                            end
                        end
                        if stage < stepper.stages
                            place_sol_in_vars(var, tmpresult, stepper);
                        end
                    end
                    for stage=1:stepper.stages
                        last_result += stepper.dt * stepper.b[stage] .* tmpki[:, stage];
                    end
                    place_sol_in_vars(var, last_result, stepper);
                end

            else # single stage methods such as Euler
                assemble_t += @elapsed(b = assemble(var, nothing, linear, allocated_vecs, dofs_per_node, dofs_per_loop, t, stepper.dt; rhs_only = true, assemble_loops=assemble_func));

                linsolve_t += @elapsed(sol = A\b);

                place_sol_in_vars(var, sol, stepper);
            end

            t += stepper.dt;

            progressPercent = Int(floor(i*100.0/stepper.Nsteps));
            if progressPercent - last2update >= 2
                last2update = progressPercent;
                if progressPercent - last10update >= 10
                    print(progressPercent);
                    last10update = progressPercent;
                else
                    print(".");
                end
            end
        end
        println("");
        end_t = Base.Libc.time();

        log_entry("Stepping took "*string(end_t-start_t)*" seconds. ("*string(assemble_t)*" for assembly, "*string(linsolve_t)*" for linear solve)");
        #display(sol);
		# outfile = "linear_sol.txt"
		# open(outfile, "w") do f
  		# 	for ii in sol
    	# 		println(f, ii)
  		# 	end
		# end # the file f is automatically closed after this block finishes
        return sol;

    else
        assemble_t = @elapsed((A, b) = assemble(var, bilinear, linear, allocated_vecs, dofs_per_node, dofs_per_loop, assemble_loops=assemble_func));
        # uncomment to look at A
        # global Amat = A;

        sol_t = @elapsed(sol = A\b);

        log_entry("Assembly took "*string(assemble_t)*" seconds");
        log_entry("Linear solve took "*string(sol_t)*" seconds");
        # display(A);
        # display(b);
        # display(sol);
        return sol;
    end
end

function nonlinear_solve(var, nlvar, bilinear, linear, stepper=nothing; assemble_loops=nothing)
    if prob.time_dependent && !(stepper === nothing)
        #TODO time dependent coefficients

        log_entry("Beginning "*string(stepper.Nsteps)*" time steps.");
        t = 0;
        start_t = Base.Libc.time();
        nl = nonlinear(100, 1e-9, 1e-9);
        init_nonlinear(nl, var, nlvar, bilinear, linear);
        for i=1:stepper.Nsteps
			println("solve for time step ", i);
			newton(nl, assemble, assemble_rhs_only, nlvar, t, stepper.dt);
            t += stepper.dt;
			nlvar[4].values = copy(nlvar[1].values);
			nlvar[5].values = copy(nlvar[2].values);
			if (i % 10 ==0)
				outfile = string("LD_u_",i);
				open(outfile, "w") do f
					for ii=1:length(nlvar[1].values)
						println(f, nlvar[1].values[ii])
					end
					close(f)
				end
				outfile = string("LD_v_",i);
				open(outfile, "w") do f
					for ii=1:length(nlvar[2].values)
						println(f, nlvar[2].values[ii])
					end
					close(f)
				end
			end
        end
        end_t = Base.Libc.time();

        log_entry("Solve took "*string(end_t-start_t)*" seconds");
        #display(sol);
        #return nlvar.values;
        return [];

    else
        start_t = Base.Libc.time();
        nl = nonlinear(100, 1e-12, 1e-12);
        init_nonlinear(nl, var, nlvar, bilinear, linear);
        newton(nl, assemble, assemble_rhs_only, nlvar);
        end_t = Base.Libc.time();

        log_entry("Solve took "*string(end_t-start_t)*" seconds");

        #return nlvar.values;
        return [];
    end
end

# assembles the A and b in Au=b
function assemble(var, bilinear, linear, allocated_vecs, dofs_per_node=1, dofs_per_loop=1, t=0, dt=0; assemble_loops=nothing, rhs_only = false)
    # If an assembly loop function was provided, use it
    if !(assemble_loops === nothing)
        return assemble_loops.func(var, bilinear, linear, allocated_vecs, dofs_per_node, dofs_per_loop, t, dt; rhs_only=rhs_only);
    end

    # Label things that were allocated externally
    b = allocated_vecs[1];
    if !rhs_only
        AI = allocated_vecs[2];
        AJ = allocated_vecs[3];
        AV = allocated_vecs[4];
    end
    # zero b
    b .= 0;

    Np = refel.Np;
    nel = mesh_data.nel;

    # Stiffness and mass are precomputed for uniform grid meshes
    precomputed_mass_stiffness = config.mesh_type == UNIFORM_GRID && config.geometry == SQUARE
    if precomputed_mass_stiffness
        wdetj = refel.wg .* geo_factors.detJ[1];
        J = geo_factors.J[1];

        if config.dimension == 1
            (RQ1, RD1) = @level_bench Level1 build_deriv_matrix(refel, J);
            TRQ1 = RQ1';
            stiffness = [(TRQ1 * diagm(wdetj) * RQ1)];
        elseif config.dimension == 2
            (RQ1, RQ2, RD1, RD2) = build_deriv_matrix(refel, J);
            (TRQ1, TRQ2) = (RQ1', RQ2');
            stiffness = [(TRQ1 * diagm(wdetj) * RQ1) , (TRQ2 * diagm(wdetj) * RQ2)];
        else
            (RQ1, RQ2, RQ3, RD1, RD2, RD3) = build_deriv_matrix(refel, J);
            (TRQ1, TRQ2, TRQ3) = (RQ1', RQ2', RQ3');
            stiffness = [(TRQ1 * diagm(wdetj) * RQ1) , (TRQ2 * diagm(wdetj) * RQ2) , (TRQ3 * diagm(wdetj) * RQ3)];
        end
        mass = (refel.Q)' * diagm(wdetj) * refel.Q;
    else
        stiffness = 0;
        mass = 0;
    end

    loop_time = Base.Libc.time();
    # Elemental loop follows elemental ordering
    for ei=1:nel
        eid = elemental_order[ei];
        loc2glb = grid_data.loc2glb[:,eid]; # global indices of this element's nodes

        if !rhs_only
            Astart = (eid-1)*Np*dofs_per_node*Np*dofs_per_node + 1; # The segment of AI, AJ, AV for this element
        end

        volargs = (var, eid, 0, grid_data, geo_factors, refel, t, dt, stiffness, mass);

        if dofs_per_node == 1
            linchunk = @level_bench Level1 linear.func(volargs);  # get the elemental linear part
            b[loc2glb] += linchunk;

            if !rhs_only
                bilinchunk = bilinear.func(volargs); # the elemental bilinear part
                #A[glb, glb] .+= bilinchunk;         # This will be very inefficient for sparse A
                for jj=1:Np
                    offset = Astart - 1 + (jj-1)*Np;
                    for ii=1:Np
                        AI[offset + ii] = loc2glb[ii];
                        AJ[offset + ii] = loc2glb[jj];
                        AV[offset + ii] = bilinchunk[ii, jj];
                    end
                end
            end

        else # more than one DOF per node
            linchunk = linear.func(volargs);
            insert_linear!(b, linchunk, loc2glb, 1:dofs_per_node, dofs_per_node);

            if !rhs_only
                bilinchunk = bilinear.func(volargs);
                @level_bench Level1  insert_bilinear!(AI, AJ, AV, Astart, bilinchunk, loc2glb, 1:dofs_per_node, dofs_per_node);
            end
        end
    end
    loop_time = Base.Libc.time() - loop_time;

    if !rhs_only
        # Build the sparse A. Uses default + to combine overlaps
        A = @level_bench Level1 sparse(AI, AJ, AV);
    end

    # Boundary conditions
    bc_time = Base.Libc.time();
    if rhs_only
        b =@level_bench Level1  apply_boundary_conditions_rhs_only(var, b, t);
    else
        (A, b) = @level_bench Level1 apply_boundary_conditions_lhs_rhs(var, A, b, t);
    end
    bc_time = Base.Libc.time() - bc_time;

    if rhs_only
        return b;

    else
        log_entry("Elemental loop time:     "*string(loop_time));
        log_entry("Boundary condition time: "*string(bc_time));
        return (A, b);
    end
end

function assemble_rhs_only(var, linear, allocated_vecs, dofs_per_node=1, dofs_per_loop=1, t=0, dt=0; assemble_loops=nothing)
    return assemble(var, nothing, linear, allocated_vecs, dofs_per_node, dofs_per_loop, t, dt; assemble_loops=assemble_loops, rhs_only=true)
end

# ######################################################
# # To be romoved. Fix nonlinear solve first.
# # assembles the A and b in Au=b
# function assemble_rhs_only(var, linear, t=0.0, dt=0.0; keep_geometric_factors = false, saved_geometric_factors = nothing)
#     Np = refel.Np;
#     nel = mesh_data.nel;
#     N1 = size(grid_data.allnodes,2);
#     multivar = typeof(var) <: Array;
#     maxvarindex = 0;
#     if multivar
#         # multiple variables being solved for simultaneously
#         dofs_per_node = 0;
#         var_to_dofs = [];
#         for vi=1:length(var)
#             tmp = dofs_per_node;
#             dofs_per_node += length(var[vi].symvar);
#             push!(var_to_dofs, (tmp+1):dofs_per_node);
#             maxvarindex = max(maxvarindex,var[vi].index);
#         end
#     else
#         # one variable
#         dofs_per_node = length(var.symvar);
#         maxvarindex = var.index;
#     end
#     Nn = dofs_per_node * N1;

#     b = zeros(Nn);

#     # Save geometric factors if desired
#     if keep_geometric_factors
#         saved_wdetj = [];
#         saved_J = [];
#     end

#     # Stiffness and mass are precomputed for uniform grid meshes
#     precomputed_mass_stiffness = config.mesh_type == UNIFORM_GRID && config.geometry == SQUARE
#     if precomputed_mass_stiffness
#         glb = grid_data.loc2glb[:,1];
#         xe = grid_data.allnodes[:,glb[:]];
#         (detJ, J) = geometric_factors(refel, xe);
#         wdetj = refel.wg .* detJ;
#         if keep_geometric_factors
#             saved_wdetj = wdetj;
#             saved_J = J;
#         end
#         if config.dimension == 1
#             (RQ1, RD1) = build_deriv_matrix(refel, J);
#             TRQ1 = RQ1';
#             stiffness = [(TRQ1 * diagm(wdetj) * RQ1)];
#         elseif config.dimension == 2
#             (RQ1, RQ2, RD1, RD2) = build_deriv_matrix(refel, J);
#             (TRQ1, TRQ2) = (RQ1', RQ2');
#             stiffness = [(TRQ1 * diagm(wdetj) * RQ1) , (TRQ2 * diagm(wdetj) * RQ2)];
#         else
#             (RQ1, RQ2, RQ3, RD1, RD2, RD3) = build_deriv_matrix(refel, J);
#             (TRQ1, TRQ2, TRQ3) = (RQ1', RQ2', RQ3');
#             stiffness = [(TRQ1 * diagm(wdetj) * RQ1) , (TRQ2 * diagm(wdetj) * RQ2) , (TRQ3 * diagm(wdetj) * RQ3)];
#         end
#         mass = (refel.Q)' * diagm(wdetj) * refel.Q;
#     else
#         stiffness = 0;
#         mass = 0;
#     end

#     # Elemental loop follows elemental ordering
#     for e=elemental_order;
#         glb = grid_data.loc2glb[:,e];       # global indices of this element's nodes for extracting values from var arrays
#         xe = grid_data.allnodes[:,glb[:]];  # coordinates of this element's nodes for evaluating coefficient functions

#         if !precomputed_mass_stiffness
#             if saved_geometric_factors === nothing
#                 (detJ, J) = geometric_factors(refel, xe);
#                 wdetj = refel.wg .* detJ;

#                 if keep_geometric_factors
#                     push!(saved_wdetj, wdetj);
#                     push!(saved_J, J);
#                 end
#             else
#                 wdetj = saved_geometric_factors[1][e];
#                 J = saved_geometric_factors[2][e];
#             end
#         end

#         rhsargs = (var, xe, glb, refel, wdetj, J, RHS, t, dt, stiffness, mass);

#         #linchunk = linear.func(args);  # get the elemental linear part
#         if dofs_per_node == 1
#             linchunk = linear.func(rhsargs);  # get the elemental linear part
#             b[glb] .+= linchunk;

#         elseif typeof(var) == Variable
#             # only one variable, but more than one dof
#             linchunk = linear.func(rhsargs);
#             insert_linear!(b, linchunk, glb, 1:dofs_per_node, dofs_per_node);

#         else
#             linchunk = linear.func(rhsargs);
#             insert_linear!(b, linchunk, glb, 1:dofs_per_node, dofs_per_node);

#         end
#     end

#     # Boundary conditions
#     bidcount = length(grid_data.bids); # the number of BIDs
#     if dofs_per_node > 1
#         if multivar
#             dofind = 0;
#             for vi=1:length(var)
#                 for compo=1:length(var[vi].symvar)
#                     dofind = dofind + 1;
#                     for bid=1:bidcount
#                         if prob.bc_type[var[vi].index, bid] == NO_BC
#                             # do nothing
#                         else
#                             b = dirichlet_bc_rhs_only(b, prob.bc_func[var[vi].index, bid][compo], grid_data.bdry[bid], t, dofind, dofs_per_node);
#                         end
#                     end
#                 end
#             end
#         else
#             for d=1:dofs_per_node
#                 #rows = ((d-1)*length(glb)+1):(d*length(glb));
#                 dofind = d;
#                 for bid=1:bidcount
#                     if prob.bc_type[var.index, bid] == NO_BC
#                         # do nothing
#                     else
#                         b = dirichlet_bc_rhs_only(b, prob.bc_func[var.index, bid][d], grid_data.bdry[bid], t, d, dofs_per_node);
#                     end
#                 end
#             end
#         end
#     else
#         for bid=1:bidcount
#             if prob.bc_type[var.index, bid] == NO_BC
#                 # do nothing
#             else
#                 b = dirichlet_bc_rhs_only(b, prob.bc_func[var.index, bid][1], grid_data.bdry[bid], t);
#             end
#         end
#     end

#     # Reference points
#     if size(prob.ref_point,1) >= maxvarindex
#         if multivar
#             posind = zeros(Int,0);
#             vals = zeros(0);
#             for vi=1:length(var)
#                 if prob.ref_point[var[vi].index,1]
#                     eii = prob.ref_point[var[vi].index, 2];
#                     tmp = (grid_data.glbvertex[eii[1], eii[2]] - 1)*dofs_per_node + var_to_dofs[vi][1];
#                     if length(prob.ref_point[var[vi].index, 3]) > 1
#                         tmp = tmp:(tmp+length(prob.ref_point[var[vi].index, 3])-1);
#                     end
#                     posind = [posind; tmp];
#                     vals = [vals; prob.ref_point[var[vi].index, 3]];
#                 end
#             end
#             if length(vals) > 0
#                 b[posind] = vals;
#             end

#         else
#             if prob.ref_point[var.index,1]
#                 eii = prob.ref_point[var.index, 2];
#                 posind = (grid_data.glbvertex[eii[1], eii[2]] - 1)*dofs_per_node + 1;
#                 if length(prob.ref_point[var.index, 3]) > 1
#                     posind = posind:(posind+length(prob.ref_point[var[vi].index, 3])-1);
#                 end
#                 b[posind] = prob.ref_point[var.index, 3];
#             end
#         end
#     end

#     return b;
# end
# ##########################################################

# Inset the single dof into the greater construct
function insert_linear!(b, bel, glb, dof, Ndofs)
    # group nodal dofs
    for d=1:length(dof)
        ind = glb.*Ndofs .- (Ndofs-dof[d]);
        ind2 = ((d-1)*length(glb)+1):(d*length(glb));

        b[ind] = b[ind] + bel[ind2];
    end
end

function insert_bilinear!(AI, AJ, AV, Astart, ael, glb, dof, Ndofs)
    Np = length(glb);
    # group nodal dofs
    for dj=1:length(dof)
        indj = glb.*Ndofs .- (Ndofs-dof[dj]);
        indj2 = ((dj-1)*Np+1):(dj*Np);
        for di=1:length(dof)
            indi = glb.*Ndofs .- (Ndofs-dof[di]);
            indi2 = ((di-1)*Np+1):(di*Np);

            #a[indi, indj] = a[indi, indj] + ael[indi2, indj2];

            for jj=1:Np
                offset = Astart + (jj-1 + Np*(dj-1))*Np*Ndofs + Np*(di-1) - 1;
                for ii=1:Np
                    AI[offset + ii] = indi[ii];
                    AJ[offset + ii] = indj[jj];
                    AV[offset + ii] = ael[indi2[ii], indj2[jj]];
                end
            end
        end
    end

end

function place_sol_in_vars(var, sol, stepper)
    # place the values in the variable value arrays
    if typeof(var) <: Array
        tmp = 0;
        totalcomponents = 0;
        for vi=1:length(var)
            totalcomponents = totalcomponents + var[vi].total_components;
        end
        for vi=1:length(var)
            components = var[vi].total_components;
            for compi=1:components
                if stepper.type == EULER_EXPLICIT
                    var[vi].values[compi,:] = sol[(compi+tmp):totalcomponents:end];
                else
                    var[vi].values[compi,:] = sol[(compi+tmp):totalcomponents:end];
                end
                tmp = tmp + 1;
            end
        end
    else
        components = var.total_components;
        for compi=1:components
            if stepper.type == EULER_EXPLICIT
                var.values[compi,:] = sol[compi:components:end];
            else
                var.values[compi,:] = sol[compi:components:end];
            end
        end
    end
end

function get_var_vals(var, vect=nothing)
    # place the variable values in a vector
    if typeof(var) <: Array
        tmp = 0;
        totalcomponents = 0;
        for vi=1:length(var)
            totalcomponents = totalcomponents + var[vi].total_components;
        end
        if vect === nothing
            vect = zeros(totalcomponents * length(var[1].values[1,:]));
        end
        for vi=1:length(var)
            components = var[vi].total_components;
            for compi=1:components
                vect[(compi+tmp):totalcomponents:end] = var[vi].values[compi,:];
                tmp = tmp + 1;
            end
        end
    else
        components = var.total_components;
        if vect === nothing
            vect = zeros(totalcomponents * length(var.values[1,:]));
        end
        for compi=1:components
            vect[compi:components:end] = var.values[compi,:];
        end
    end

    return vect;
end

end #module


text:
function prepare_parameters_NN_training(𝒟train, f, Nz, g, α, ν₀, ν₋, Riᶜ, ΔRi, Pr, κ, conditions)
    H = abs(𝒟train.uw.z[end] - 𝒟train.uw.z[1])
    τ = abs(𝒟train.t[:,1][end] - 𝒟train.t[:,1][1])
    u_scaling = 𝒟train.scalings["u"]
    v_scaling = 𝒟train.scalings["v"]
    T_scaling = 𝒟train.scalings["T"]
    uw_scaling = 𝒟train.scalings["uw"]
    vw_scaling = 𝒟train.scalings["vw"]
    wT_scaling = 𝒟train.scalings["wT"]

    if conditions.modified_pacanowski_philander
        constants = (H=H, τ=τ, f=f, Nz=Nz, g=g, α=α, ν₀=ν₀, ν₋=ν₋, Riᶜ=Riᶜ, ΔRi=ΔRi, Pr=Pr)
    elseif conditions.convective_adjustment
        constants = (H=H, τ=τ, f=f, Nz=Nz, g=g, α=α, κ=κ)
    else
        constants = (H=H, τ=τ, f=f, Nz=Nz, g=g, α=α)
    end
    scalings = (u=u_scaling, v=v_scaling, T=T_scaling, uw=uw_scaling, vw=vw_scaling, wT=wT_scaling)
    derivatives = (cell=Float32.(Dᶜ(Nz, 1 / Nz)), face=Float32.(Dᶠ(Nz, 1 / Nz)))

    filters = (cell=WindMixing.smoothing_filter(Nz, 3), face=WindMixing.smoothing_filter(Nz+1, 3), interior=WindMixing.smoothing_filter(Nz-1, 3))
    return constants, scalings, derivatives, filters
end

function predict_uw(NN, x, BCs, conditions, scalings, constants, derivatives, filters)
    Nz, H, τ, f = constants.Nz, constants.H, constants.τ, constants.f
    uw_scaling, vw_scaling, wT_scaling = scalings.uw, scalings.vw, scalings.wT
    σ_uw, σ_vw, σ_wT = uw_scaling.σ, vw_scaling.σ, wT_scaling.σ
    μ_u, μ_v, σ_u, σ_v, σ_T = scalings.u.μ, scalings.v.μ, scalings.u.σ, scalings.v.σ, scalings.T.σ
    D_cell, D_face = derivatives.cell, derivatives.face

    u = @view x[1:Nz]
    v = @view x[Nz + 1:2Nz]
    T = @view x[2Nz + 1:3Nz]

    interior = NN(x)

    if conditions.smooth_NN
        interior = filters.interior * interior
    end

    if conditions.zero_weights
        uw = [0f0; interior; 0f0]
    else
        uw = [BCs.uw.bottom; interior; BCs.uw.top]
    end

    if conditions.modified_pacanowski_philander
        ϵ = 1f-7
        ∂u∂z = D_face * u
        ∂v∂z = D_face * v
        ∂T∂z = D_face * T
        Ri = local_richardson.(∂u∂z .+ ϵ, ∂v∂z .+ ϵ, ∂T∂z .+ ϵ, constants.H, constants.g, constants.α, scalings.u.σ, scalings.v.σ, scalings.T.σ)

        if conditions.smooth_Ri
            Ri = filters.face * Ri
        end

        ν = constants.ν₀ .+ constants.ν₋ .* tanh_step.((Ri .- constants.Riᶜ) ./ constants.ΔRi)

        if conditions.zero_weights
            ν∂u∂z = [-(BCs.uw.bottom - scalings.uw(0f0)); σ_u / σ_uw / H .* ν[2:end-1] .* ∂u∂z[2:end-1]; -(BCs.uw.top - scalings.uw(0f0))]
        else
            ν∂u∂z = σ_u / σ_uw / H .* ν .* ∂u∂z
        end

        return uw .- ν∂u∂z
    else
        return uw
    end
end

function predict_vw(NN, x, BCs, conditions, scalings, constants, derivatives, filters)
    Nz, H, τ, f = constants.Nz, constants.H, constants.τ, constants.f
    uw_scaling, vw_scaling, wT_scaling = scalings.uw, scalings.vw, scalings.wT
    σ_uw, σ_vw, σ_wT = uw_scaling.σ, vw_scaling.σ, wT_scaling.σ
    μ_u, μ_v, σ_u, σ_v, σ_T = scalings.u.μ, scalings.v.μ, scalings.u.σ, scalings.v.σ, scalings.T.σ
    D_cell, D_face = derivatives.cell, derivatives.face

    u = @view x[1:Nz]
    v = @view x[Nz + 1:2Nz]
    T = @view x[2Nz + 1:3Nz]

    interior = NN(x)

    if conditions.smooth_NN
        interior = filters.interior * interior
    end

    if conditions.zero_weights
        vw = [0f0; interior; 0f0]
    else
        vw = [BCs.vw.bottom; interior; BCs.vw.top]
    end

    if conditions.modified_pacanowski_philander
        ϵ = 1f-7
        ∂u∂z = D_face * u
        ∂v∂z = D_face * v
        ∂T∂z = D_face * T
        Ri = local_richardson.(∂u∂z .+ ϵ, ∂v∂z .+ ϵ, ∂T∂z .+ ϵ, constants.H, constants.g, constants.α, scalings.u.σ, scalings.v.σ, scalings.T.σ)

        if conditions.smooth_Ri
            Ri = filters.face * Ri
        end

        ν = constants.ν₀ .+ constants.ν₋ .* tanh_step.((Ri .- constants.Riᶜ) ./ constants.ΔRi)
        if conditions.zero_weights
            ν∂v∂z = [-(BCs.vw.bottom - scalings.vw(0f0)); σ_v / σ_vw / H .* ν[2:end-1] .* ∂v∂z[2:end-1]; -(BCs.vw.top - scalings.vw(0f0))]
        else
            ν∂v∂z = σ_v / σ_vw / H .* ν .* ∂v∂z
        end

        return vw .- ν∂v∂z
    else
        return vw
    end
end

function predict_wT(NN, x, BCs, conditions, scalings, constants, derivatives, filters)
    Nz, H, τ, f = constants.Nz, constants.H, constants.τ, constants.f
    uw_scaling, vw_scaling, wT_scaling = scalings.uw, scalings.vw, scalings.wT
    σ_uw, σ_vw, σ_wT = uw_scaling.σ, vw_scaling.σ, wT_scaling.σ
    μ_u, μ_v, σ_u, σ_v, σ_T = scalings.u.μ, scalings.v.μ, scalings.u.σ, scalings.v.σ, scalings.T.σ
    D_cell, D_face = derivatives.cell, derivatives.face

    u = @view x[1:Nz]
    v = @view x[Nz + 1:2Nz]
    T = @view x[2Nz + 1:3Nz]

    interior = NN(x)

    if conditions.smooth_NN
        interior = filters.interior * interior
    end

    if conditions.zero_weights
        wT = [0f0; interior; 0f0]
    else
        wT = [BCs.wT.bottom; interior; BCs.wT.top]
    end

    if conditions.modified_pacanowski_philander
        ϵ = 1f-7
        ∂u∂z = D_face * u
        ∂v∂z = D_face * v
        ∂T∂z = D_face * T
        Ri = local_richardson.(∂u∂z .+ ϵ, ∂v∂z .+ ϵ, ∂T∂z .+ ϵ, constants.H, constants.g, constants.α, scalings.u.σ, scalings.v.σ, scalings.T.σ)

        if conditions.smooth_Ri
            Ri = filters.face * Ri
        end

        ν = constants.ν₀ .+ constants.ν₋ .* tanh_step.((Ri .- constants.Riᶜ) ./ constants.ΔRi)
        if conditions.zero_weights
            ν∂T∂z = [-(BCs.wT.bottom - scalings.wT(0f0)); σ_T / σ_wT / H .* ν[2:end-1] ./ constants.Pr .* ∂T∂z[2:end-1]; -(BCs.wT.top - scalings.wT(0f0))]
        else
            ν∂T∂z = σ_T / σ_wT / H .* ν .* ∂T∂z ./ constants.Pr
        end

        return wT .- ν∂T∂z
    elseif conditions.convective_adjustment
        ∂T∂z = D_face * T
        κ∂T∂z = σ_T / σ_wT / H .* κ .* min.(0f0, ∂T∂z)
        return wT .- κ∂T∂z
    else
        return wT
    end
end

function predict_NN(NN, x, y)
    interior = NN(x)
    return [y[1]; interior; y[end]]
end

function save_NN_weights(weights, FILE_PATH, filename)
    NN_params = Dict(:weights => weights)
    bson(joinpath(FILE_PATH, "$filename.bson"), NN_params)
end

function loss_flux_gradient(NN, )

end

function prepare_NN_training_data(𝒟, NN_type, derivatives)
    @inline training_data_BCs(𝒟, i) = (uw=(top=𝒟.uw.scaled[end,i], bottom=𝒟.uw.scaled[1,i]),
                                       vw=(top=𝒟.vw.scaled[end,i], bottom=𝒟.vw.scaled[1,i]),
                                       wT=(top=𝒟.wT.scaled[end,i], bottom=𝒟.wT.scaled[1,i]))
    @inline training_data_uw(𝒟, i) = ((profile=𝒟.uvT_scaled[:,i], BCs=training_data_BCs(𝒟,i)), (flux=𝒟.uw.scaled[:,i], flux_gradient=calculate_flux_gradient(𝒟.uw.scaled[:,i], derivatives)))
    @inline training_data_vw(𝒟, i) = ((profile=𝒟.uvT_scaled[:,i], BCs=training_data_BCs(𝒟,i)), (flux=𝒟.vw.scaled[:,i], flux_gradient=calculate_flux_gradient(𝒟.vw.scaled[:,i], derivatives)))
    @inline training_data_wT(𝒟, i) = ((profile=𝒟.uvT_scaled[:,i], BCs=training_data_BCs(𝒟,i)), (flux=𝒟.wT.scaled[:,i], flux_gradient=calculate_flux_gradient(𝒟.wT.scaled[:,i], derivatives)))

    if NN_type == "uw"
        data = [training_data_uw(𝒟, i) for i in 1:size(𝒟.uw.scaled, 2)]
    elseif NN_type == "vw"
        data = [training_data_vw(𝒟, i) for i in 1:size(𝒟.vw.scaled, 2)]
    else
        data = [training_data_wT(𝒟, i) for i in 1:size(𝒟.wT.scaled, 2)]
    end
    return shuffle(data)
end

function calculate_flux_gradient(flux, derivatives)
    return derivatives.cell * flux
end

function train_NN(NN, 𝒟train, optimizers, train_epochs, FILE_PATH, NN_type; ν₀=1f-4, ν₋=1f-1, ΔRi=1f0, Riᶜ=0.25, Pr=1f0, κ=10f0, f=1f-4, α=1.67f-4, g=9.81f0,
                  modified_pacanowski_philander=false, convective_adjustment=false, smooth_profile=false, smooth_NN=false, smooth_Ri=false, train_gradient=false,
                  zero_weights=false, gradient_scaling=1f-4)
    Nz = length(𝒟train.u.z)

    conditions = (modified_pacanowski_philander=modified_pacanowski_philander, convective_adjustment=convective_adjustment,
                  smooth_profile=smooth_profile, smooth_NN=smooth_NN, smooth_Ri=smooth_Ri,
                  train_gradient=train_gradient, zero_weights=zero_weights)

    constants, scalings, derivatives, filters = prepare_parameters_NN_training(𝒟train, f, Nz, g, α, ν₀, ν₋, Riᶜ, ΔRi, Pr, κ, conditions)
    training_data = prepare_NN_training_data(𝒟train, NN_type, derivatives)

    function NN_loss(input, output)
        if NN_type == "uw"
            NN_flux = predict_uw(NN, input.profile, input.BCs, conditions, scalings, constants, derivatives, filters)
        elseif NN_type == "vw"
            NN_flux = predict_vw(NN, input.profile, input.BCs, conditions, scalings, constants, derivatives, filters)
        else
            NN_flux = predict_wT(NN, input.profile, input.BCs, conditions, scalings, constants, derivatives, filters)
        end
        ∂z_NN_flux = calculate_flux_gradient(NN_flux, derivatives)
        # return loss(NN_flux, output.flux)

        return loss(NN_flux, output.flux) + gradient_scaling * loss(output.flux_gradient, ∂z_NN_flux)
    end

    function total_loss(training_data)
        return mean([NN_loss(data[1], data[2]) for data in training_data])
    end

    # loss(x, y) = loss_gradient(x, y, calculate_flux_gradient(x, derivatives), calculate_flux_gradient(y, derivatives), gradient_scaling)

    for i in 1:length(optimizers), epoch in 1:train_epochs[i]
        opt = optimizers[i]
        function cb()
            @info "$NN_type NN, loss = $(total_loss(training_data)), opt $i/$(length(optimizers)), epoch $epoch/$(train_epochs[i])"
        end
        Flux.train!(NN_loss, Flux.params(NN), training_data, opt, cb=Flux.throttle(cb,10))
        write_data_NN_training(FILE_PATH, total_loss(training_data), NN)
    end

    return Flux.destructure(NN)[1]
end


text:
function layerBeamPatternRose(bear::BallTreeDensity;
                              scale::Float64=1.0,
                              c=colorant"magenta",
                              wRr=TU.R(0.0),
                              wTRr=zeros(2) )
  #
  tp = reshape([0:0.01:(2pi);], 1, :)
  belp = scale*bear(tp)
  Gadfly.plot(x=tp, y=belp, Geom.path)
  belRose = zeros(2, length(tp))
  belRose[1,:] = belp

  idx = 0
  for rRc in TU.R.(tp)
    idx += 1
    belRose[:,idx] = wRr*rRc*belRose[:,idx] + wTRr
  end

  Gadfly.layer(x=belRose[1,:], y=belRose[2,:], Geom.path, Theme(default_color=c))
end

function plotSandsharkFromDFG(dfg::G; radix::Float64=1.8, scale::Float64=0.4) where G <: AbstractDFG
  # bearing range factors
  brf = ls(dfg, Pose2Point2BearingRange)
  #variables and poses connected to br factors
  brvars = union(map(x->ls(dfg, x), brf)...)
  brposes = setdiff(brvars, ls(dfg, r"l"))

  # get body frame beam estimate
  brfacts = Dict{Symbol, BallTreeDensity}()

  PL =  drawPosesLandms(dfg, spscale=2.5)

  for i in 1:length(brposes)
    fct = getFactorType(getFactor(dfg, brf[1]))
    brfacts[brposes[i]] = fct.bearing
    rVo = getKDEMax(getKDE(dfg, brposes[i]))
    pll = AMP.plotKDECircular([fct.bearing;], rVo=rVo, radix=radix, scale=scale)
    for lyr in pll.layers
      push!(PL.layers, lyr)
    end
  end

  PL
end

function drawPosesLandmarksAndOdo(fg::G,
                                  ppbrDict::Dict,
                                  navkeys::Array,
                                  X::Array,
                                  Y::Array,
                                  lblX::Array,
                                  lblY::Array,
                                  gTitle::String="",
                                  drawhist=true  ) where G <: AbstractDFG
    #
    PLL = []
    # xx = sortVarNested(ls(fg, r"x"))
    # idx = 0
    # for idx in 1:length(xx)
    #   ep = epochs[idx]
    #   theta = getKDEMax(getVertKDE(fg, xx[idx]))
    #   wRr = TU.R(theta[3])
    #   pll = layerBeamPatternRose(ppbrDict[ep].bearing, wTRr=theta[1:2], scale=5.0 , wRr=wRr)
    #   push!(PLL, pll)
    # end
    #
    # pllandmarks = drawPosesLandms(fg, spscale=2.5, drawhist=drawhist)
    pllandmarks = plotSandsharkFromDFG(fg, scale=1.8)

    # Add odo
    navdf = DataFrame(
      ts = navkeys,
      x = X,
      y = Y
    )
    # pl = Gadfly.layer(navdf, x=:x, y=:y, Geom.path())
    push!(pllandmarks.layers, Gadfly.layer(navdf, x=:x, y=:y, Geom.path(), Theme(default_color=colorant"gray80"))[1])
    lbldf = DataFrame(
      ts = lblkeys .- lblkeys[1],
      x = lblX,
      y = lblY
    )
    push!(pllandmarks.layers, Gadfly.layer(lbldf, x=:x, y=:y, Geom.path(), Theme(default_color=colorant"magenta"))[1])
    # Gadfly.plot(pl.layers)

    # Make X/Y range same so no distorted
    # push!(PLL, Coord.Cartesian(xmin=-160.0,xmax=10.0,ymin=-135.0,ymax=35.0))
    pla = Gadfly.plot([PLL...;pllandmarks.layers; ]...,
        Guide.manual_color_key("Legend", ["Dead reckoning", "LBL Path", "Non-Gaussian SLAM"], ["gray80", "green", "light blue"]),
        Guide.title(gTitle))
    return pla
end

#


text:
Pkg.update()
Pkg.add("PyPlot")
Pkg.add("StatsFuns")
Pkg.add("SpecialFunctions")
Pkg.add("Distributions")
Pkg.add("PDMats")
Pkg.add("ProgressMeter")
Pkg.add("DataFrames")
Pkg.add("HDF5")
Pkg.add("JLD")
Pkg.add("IJulia")


text:
abstract type AbstractLoadFormulation <: AbstractDeviceFormulation end

abstract type AbstractControllablePowerLoadFormulation <: AbstractLoadFormulation end

struct StaticPowerLoad <: AbstractLoadFormulation end

struct InterruptiblePowerLoad <: AbstractControllablePowerLoadFormulation end

struct DispatchablePowerLoad <: AbstractControllablePowerLoadFormulation end

########################### dispatchable load variables ####################################
function activepower_variables!(canonical::Canonical,
                               devices::IS.FlattenIteratorWrapper{L}) where L<:PSY.ElectricLoad
    add_variable(canonical,
                 devices,
                 Symbol("P_$(L)"),
                 false,
                 :nodal_balance_active, -1.0;
                 ub_value = x -> PSY.get_maxactivepower(x),
                 lb_value = x -> 0.0)

    return

end

function reactivepower_variables!(canonical::Canonical,
                                 devices::IS.FlattenIteratorWrapper{L}) where L<:PSY.ElectricLoad
    add_variable(canonical,
                 devices,
                 Symbol("Q_$(L)"),
                 false,
                 :nodal_balance_reactive, -1.0;
                 ub_value = x -> PSY.get_maxreactivepower(x),
                 lb_value = x -> 0.0)

    return

end

function commitment_variables!(canonical::Canonical,
                              devices::IS.FlattenIteratorWrapper{L}) where L<:PSY.ElectricLoad

    add_variable(canonical,
                 devices,
                 Symbol("ON_$(L)"),
                 true)

    return

end

####################################### Reactive Power Constraints #########################
"""
Reactive Power Constraints on Loads Assume Constant PowerFactor
"""
function reactivepower_constraints!(canonical::Canonical,
                                   devices::IS.FlattenIteratorWrapper{L},
                                   device_formulation::Type{<:AbstractControllablePowerLoadFormulation},
                                   system_formulation::Type{<:PM.AbstractPowerModel}) where L<:PSY.ElectricLoad

    time_steps = model_time_steps(canonical)
    key = Symbol("reactive_$(L)")
    canonical.constraints[key] = JuMPConstraintArray(undef, (PSY.get_name(d) for d in devices), time_steps)

    for t in time_steps, d in devices
        name = PSY.get_name(d)
        pf = sin(atan((PSY.get_maxreactivepower(d)/PSY.get_maxactivepower(d))))
        canonical.constraints[key][PSY.get_name(d), t] = JuMP.@constraint(canonical.JuMPmodel,
                        canonical.variables[Symbol("Q_$(L)")][name, t] == canonical.variables[Symbol("P_$(L)")][name, t]*pf)
    end

    return

end

######################## output constraints without Time Series ############################
function _get_time_series(canonical::Canonical,
                          devices::IS.FlattenIteratorWrapper{<:PSY.ElectricLoad})

    initial_time = model_initial_time(canonical)
    use_forecast_data = model_uses_forecasts(canonical)
    time_steps = model_time_steps(canonical)
    device_total = length(devices)
    ts_data_active = Vector{Tuple{String, Int64, Float64, Vector{Float64}}}(undef, device_total)
    ts_data_reactive = Vector{Tuple{String, Int64, Float64, Vector{Float64}}}(undef, device_total)

    for (ix, device) in enumerate(devices)
        bus_number = PSY.get_number(PSY.get_bus(device))
        name = PSY.get_name(device)
        active_power = use_forecast_data ? PSY.get_maxactivepower(device) : PSY.get_activepower(device)
        reactive_power = use_forecast_data ? PSY.get_maxreactivepower(device) : PSY.get_reactivepower(device)
        if use_forecast_data
            forecast = PSY.get_forecast(PSY.Deterministic,
                                        device,
                                        initial_time,
                                        "get_maxactivepower",
                                        length(time_steps))
            ts_vector = TS.values(PSY.get_data(forecast))
        else
            ts_vector = ones(time_steps[end])
        end
        ts_data_active[ix] = (name, bus_number, active_power, ts_vector)
        ts_data_reactive[ix] = (name, bus_number, reactive_power, ts_vector)
    end

    return ts_data_active, ts_data_reactive

end

function activepower_constraints!(canonical::Canonical,
                                 devices::IS.FlattenIteratorWrapper{L},
                                 device_formulation::Type{DispatchablePowerLoad},
                                 system_formulation::Type{<:PM.AbstractPowerModel}) where L<:PSY.ElectricLoad

    parameters = model_has_parameters(canonical)
    use_forecast_data = model_uses_forecasts(canonical)

    if !parameters && !use_forecast_data
        range_data = [(PSY.get_name(d), (min = 0.0, max = PSY.get_activepower(d))) for d in devices]
        device_range(canonical,
                    range_data,
                    Symbol("activerange_$(L)"),
                    Symbol("P_$(L)"))
        return
    end

    ts_data_active, _ = _get_time_series(canonical, devices)
    if parameters
        device_timeseries_param_ub(canonical,
                                   ts_data_active,
                                   Symbol("active_$(L)"),
                                   UpdateRef{L}("get_maxactivepower"),
                                   Symbol("P_$(L)"))
    else
        device_timeseries_ub(canonical,
                            ts_data_active,
                            Symbol("active_$(L)"),
                            Symbol("P_$(L)"))
    end

    return

end

function activepower_constraints!(canonical::Canonical,
                                 devices::IS.FlattenIteratorWrapper{L},
                                 device_formulation::Type{InterruptiblePowerLoad},
                                 system_formulation::Type{<:PM.AbstractPowerModel}) where L<:PSY.ElectricLoad

    parameters = model_has_parameters(canonical)
    use_forecast_data = model_uses_forecasts(canonical)

    if !parameters && !use_forecast_data
        range_data = [(PSY.get_name(d), (min = 0.0, max = PSY.get_activepower(d))) for d in devices]
        device_range(canonical,
                    range_data,
                    Symbol("activerange_$(L)"),
                    Symbol("P_$(L)"))
        return
    end

    ts_data_active, _ = _get_time_series(canonical, devices)
    if parameters
        device_timeseries_ub_bigM(canonical,
                                 ts_data_active,
                                 Symbol("active_$(L)"),
                                 Symbol("P_$(L)"),
                                 UpdateRef{L}("get_maxactivepower"),
                                 Symbol("ON_$(L)"))
    else
        device_timeseries_ub_bin(canonical,
                                ts_data_active,
                                Symbol("active_$(L)"),
                                Symbol("P_$(L)"),
                                Symbol("ON_$(L)"))
    end

    return

end

########################## Addition of to the nodal balances ###############################
function nodal_expression!(canonical::Canonical,
                           devices::IS.FlattenIteratorWrapper{L},
                           system_formulation::Type{<:PM.AbstractPowerModel}) where L<:PSY.ElectricLoad

    ts_data_active, ts_data_reactive = _get_time_series(canonical, devices)
    parameters = model_has_parameters(canonical)

    if parameters
        include_parameters(canonical,
                        ts_data_active,
                        UpdateRef{L}("get_maxactivepower"),
                        :nodal_balance_active,
                        -1.0)
        include_parameters(canonical,
                        ts_data_reactive,
                        UpdateRef{L}("get_maxactivepower"),
                        :nodal_balance_reactive,
                        -1.0)
        return
    end

    for t in model_time_steps(canonical)
        for device_value in ts_data_active
            _add_to_expression!(canonical.expressions[:nodal_balance_active],
                            device_value[2],
                            t,
                            -device_value[3]*device_value[4][t])
        end
        for device_value in ts_data_reactive
            _add_to_expression!(canonical.expressions[:nodal_balance_reactive],
                            device_value[2],
                            t,
                            -device_value[3]*device_value[4][t])
        end
    end

    return

end

function nodal_expression!(canonical::Canonical,
                           devices::IS.FlattenIteratorWrapper{L},
                           system_formulation::Type{<:PM.AbstractActivePowerModel}) where L<:PSY.ElectricLoad

    parameters = model_has_parameters(canonical)
    ts_data_active, _ = _get_time_series(canonical, devices)

    if parameters
        include_parameters(canonical,
                        ts_data_active,
                        UpdateRef{L}("get_maxactivepower"),
                        :nodal_balance_active,
                        -1.0)
        return
    end

    for t in model_time_steps(canonical), device_value in ts_data_active
        _add_to_expression!(canonical.expressions[:nodal_balance_active],
                            device_value[2],
                            t,
                            -device_value[3]*device_value[4][t])
    end

    return
end

############################## FormulationControllable Load Cost ###########################
function cost_function(canonical::Canonical,
                       devices::IS.FlattenIteratorWrapper{L},
                       device_formulation::Type{DispatchablePowerLoad},
                       system_formulation::Type{<:PM.AbstractPowerModel}) where L<:PSY.ControllableLoad

    add_to_cost(canonical,
                devices,
                Symbol("P_$(L)"),
                :variable,
                -1.0)

    return

end

function cost_function(canonical::Canonical,
                       devices::IS.FlattenIteratorWrapper{L},
                       device_formulation::Type{InterruptiblePowerLoad},
                       system_formulation::Type{<:PM.AbstractPowerModel}) where L<:PSY.ControllableLoad

    add_to_cost(canonical,
                devices,
                Symbol("ON_$(L)"),
                :fixed,
                -1.0)

    return

end


text:
import Pkg
import JSON

@show LOAD_PATH
@show DEPOT_PATH

# JSON.parse - string or stream to Julia data structures
s = "{\"a_number\" : 5.0, \"an_array\" : [\"string\", 9]}"
j = JSON.parse(s)
#  Dict{AbstractString,Any} with 2 entries:
#    "an_array" => {"string",9}
#    "a_number" => 5.0

# JSON.json - Julia data structures to a string
@show JSON.json([2,3])
#  "[2,3]"
@show JSON.json(j)


text:
# ============================================================================ #
mutable struct Player
    h::PyCall.PyObject
    loc::Vector{Float64} # the actaul x,y location of the player
    facing::Bearing
end

function Player(x::Integer, y::Integer, facing::Bearing=East)
    pts = [
        -0.2 -0.2;
        -0.2 +0.2;
        +0.2 +0.0;
    ]
    loc = [x+.5 y+.5]
    h = plt.Polygon(pts .+ loc, facecolor="black")
    p = Player(h, vec(loc), facing)

    return set_to_face!(p, facing)
end

# get the x,y location of the player in board corrdinates (i.e. which grid square)
location(p::Player) = (Int(p.loc[1]-.5), Int(p.loc[2]-.5))
clear!(p::Player) = p.h.remove()

moveto!(p::Player, x::Real, y::Real) = moveto!(p, [x, y])
function moveto!(p::Player, xy::Vector{<:Real})
    df = xy .+ .5 .- p.loc
    p.h.set_xy(p.h.get_xy() .+ df')
    p.loc = xy .+ .5
    return p
end

function set_to_face!(p::Player, facing::Bearing)
    if facing != p.facing
        df = rem(Int(facing) - Int(p.facing), 4)
        for k = 1:abs(df)
            turn!(p, Direction(sign(df)))
        end
    end
    return p
end

forward!(p::Player) = moveto!(p, destination(p))

function destination(p::Player)
    if p.facing == North
        return p.loc .+ [-.5, .5]
    elseif p.facing == East
        return p.loc .+ [.5, -.5]
    elseif p.facing == South
        return p.loc .+ [-.5, -1.5]
    elseif p.facing == West
        return p.loc .+ [-1.5, -.5]
    end
end

function turn!(p::Player, d::Direction)
    if d == Left
        xfm = [0 -1; 1 0]
    elseif d == Right
        xfm = [0 1; -1 0]
    else
        # we only accept relative turns (i.e. LEFT and RIGHT)
        xfm = [1 0; 0 1]
    end
    pts = (permutedims(p.h.get_xy(), (2, 1)) .- p.loc)
    p.h.set_xy(permutedims(xfm * pts .+ p.loc, (2, 1)))

    p.facing = Bearing(mod(Int(p.facing) + Int(d), 4))

    return p
end


text:
const LSS{T} = Vector{<:AbstractSubsystem{T}}

Base.sort(s::LSS) = sort(s, by = x -> x.index)

function _getdimarrays(subs::LSS)
    N = length(subs)
    n = zeros(Int,N)
    m = zeros(Int,N)
    p = zeros(Int,N)
    for i in 1:N
        n[i] = subs[i].nx
        m[i] = subs[i].nu
        p[i] = subs[i].ny
    end
    (n,m,p)
end

function tomonolith(subs::LSS{T}) where T
    # sorted = sort(subs)
    n,m,p = _getdimarrays(subs)
    A = BlockArray{T}(zeros(sum(n), sum(n)), n, n)
    B = BlockArray{T}(zeros(sum(n), sum(m)), n, m)
    C = BlockArray{T}(zeros(sum(p), sum(n)), p, n)
    for s in subs
        i = s.index
        view(A, Block(i,i)) .= s.A
        view(B, Block(i,i)) .= s.B
        view(C, Block(i,i)) .= s.C
        for j in 1:length(s.neighbours)
           view(A, Block(i, s.neighbours[j].index)) .= s.Aij[j]
        end
    end
    (A=A, B=B, C=C)
end


text:
using Measurements

d(a_in, a_out) = (1 - a_in*a_in) / (1 - a_in*a_out) * a_out - a_in
d̃(a_in, a_out) = a_out - a_in

# Au(111)
# a_in, a_out...

P_in = [
    # Au(111)
     22    0  -22
    # Au(poly.)
     -4    2    6
     -4    2    6
     -4    2    6
    # Au(poly.) (not reported)
     -4    2    6
    # Al(poly.)
      0    0    0
    # Au(poly.) (not reported)
     -4    2    6
     -4    2    6
     -4    2    6
    # Cu(332)
      0    0    0
      0    0    0
    # Ag(110)
      3    0   -3
      3    0   -3
    # Au(110)
     27    0  -24
     27    0  -24
    # Au(110)
    missing  0  missing
    missing  0  missing
] ./ 100

σ_in = [
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    missing  5  missing
    missing  5  missing
] ./ 100

P_out = [
    # Au(111)
    -29 -31 -35
    # Au(poly.)
    -40 -38 -35
    -39 -36 -32
    -61 -57 -55
    # Au(poly.)
      8  14  17
    # Al(poly.)
     14.9 15 15.1
    # Au(poly.)
    -17 -11 -6
    -20 -14 -10
    -22 -18 -14
    # Cu(332)
    -6.3 -6.7 -6
     11  12  12
    # Ag(110)
     -7  -9 -12
    6.2 7.1 6.3
    # Au(111)
     22  -8 -35
     34   8 -16
    # Au(111)
    missing  -12  missing
    missing    4  missing
] ./ 100

σ_out = [
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    5  5  5
    missing  5  missing
    missing  5  missing
] ./ 100

# Polycrystalline
# a_in, a_out...

D(a_in, a_out, idx=2) = (a_out[idx] - a_in[idx]) / (1 - a_in[idx]*a_out[idx])
A(idx, a_in, a_out) = (1 - a_in[idx]*a_out[idx]) * D(a_in, a_out) - (a_out[idx] - a_in[idx])
ΔT(a_in, a_out) = (A(3, a_in, a_out) + A(1, a_in, a_out)) / (a_out[3] - a_out[1])
ΔD(a_in, a_out) = ΔT(a_in, a_out) * (a_out[3] + a_out[1])/2 - (A(3, a_in, a_out) - A(1, a_in, a_out)) / 2

D2(idx, a_in, a_out) = (a_out[idx] - a_in[idx] + (idx-2)*(ΔT(a_in, a_out)*a_out[idx] - ΔD(a_in, a_out))) / (1 - a_in[idx]*a_out[idx])

Dcons(a_in, a_out) = (D(a_in, a_out, 1) + D(a_in, a_out, 3)) / 2 - D(a_in, a_out, 2)


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id: DoublePipeIncr.mso								$
#*------------------------------------------------------------------
type Main_DoublePipe
	Main_DoublePipe()=begin
		new(
			DoublePipe_HeatTransfer(Dict{Symbol,Any}(
				:Brief=>"Double Pipe Heat Transfer"
			)),
			DoublePipe_PressureDrop(Dict{Symbol,Any}(
				:Brief=>"Double Pipe Pressure Drop"
			)),
			Physical_Properties(Dict{Symbol,Any}(
				:Brief=>"Double Pipe Properties"
			)),
			[:HeatTransfer,:PressureDrop,:Properties,]
		)
	end
	HeatTransfer::DoublePipe_HeatTransfer
	PressureDrop::DoublePipe_PressureDrop
	Properties::Physical_Properties
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Main_DoublePipe
function atributes(in::Main_DoublePipe,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="to be documented"
	fields[:Info]="to be documented"
	drive!(fields,_)
	return fields
end
Main_DoublePipe(_::Dict{Symbol,Any})=begin
	newModel=Main_DoublePipe()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
println("\n\n\nRunning tests ...")

using LaTeX
import LaTeX: getrequirements
using Base.Test

# Check that package dependency is working
@test haskey(
    getrequirements(Section("Code", [Code("""test_getrequirements()""")])),
    "texments")

# Test a minimal document
doc = document([
    DocumentClass("article", ["11pt", "letterpaper"]),
    Date(2015, 11, 23),
    Title("Sample Article"),
    Author("Andi Andromeda"),
    Section("Code", [Code("""# minimal example""")])])
@test split(doc, '\n')[1] == "\\documentclass[11pt,letterpaper]{article}"
@test contains(doc, "\\date{2015-11-23}")
@test contains(doc, "\\author{Andi Andromeda}")
@test contains(doc, "\\section")  # articles have section as top-level
@test contains(doc, "\\usestyle{default}")  # needed for code

# Test a minimal report
rep = report(Section("Test", "Lorum ipsum."))
@test split(rep, '\n')[1] == "\\documentclass[11pt,a4paper]{report}"
@test contains(rep, "\\chapter")  # reports have chapter as top-level

println("   ... done running tests!")