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

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Revise
using ADSeismic
using ADCME
using JLD2
ADCME.options.customop.verbose = false
mpi_init()
M = N = Int(sqrt(mpi_size()))
n = 10÷M

param = MPIElasticPropagatorParams(NX=10, NY=10, n=n, NSTEP=10, DELTAT=1e-4,
    DELTAX=1.0, DELTAY=1.0,
    USE_PML_XMIN=true, USE_PML_XMAX = true, USE_PML_YMIN = true, USE_PML_YMAX = true)
compute_PML_Params!(param)
vp = 3300.
vs = 3300. / 1.732
rho = 2800.
λ, μ, ρ = compute_lame_parameters(param, vp, vs, rho)
source = Ricker(param, 50.0, 200.0)

# close("all")
# plot(run(sess, source))
# savefig("test.png")
srci = [div(param.NX,5)]
srcj = [div(param.NY,2)]
srctype = [0]
srcv = reshape(source, :, 1)
src = MPIElasticSource(param, srci, srcj, srctype, srcv)

propagator = MPIElasticPropagatorSolver(param, src, ρ, λ, μ)

sess = Session(); init(sess)
@info "starting..."
Vx, Vy, Sxx, Syy, Sxy = run(sess, [propagator.vx, propagator.vy, propagator.sigmaxx, propagator.sigmayy, propagator.sigmaxy])

# Vx = run(sess, propagator.vx)

# Vx[2,:,:]
@info "saving..."
@save "data/dat$(mpi_rank())-$(mpi_size()).jld2" Vx Vy Sxx Syy Sxy param

if mpi_size()>1
    mpi_finalize()
end
# if mpi_rank()==0 && mpi_size()==1
#     # close("all")
#     # pcolormesh(rcvv, cmap="gray")
#     # colorbar()
#     # xlabel("Location")
#     # ylabel("Time")
#     # savefig("receiver.png")

#     close("all")
#     visualize_wavefield(Vx, param)

# end

ext: jl


text:
using Combinatorics, Polyhedra
using CDDLib

v = vrep(collect(permutations([0, 1, 2, 3])))

p4 = polyhedron(v, CDDLib.Library())

# orthogonal basis for projection hyperplane
v1 = [1, -1,  0,  0]
v2 = [1,  1, -2,  0]
v3 = [1,  1,  1, -3]

p3 = project(p4, [v1 v2 v3])

# using MeshCat
# vis = Visualizer()
# open(vis)
m = Polyhedra.Mesh(p3);
setobject!(vis, m,
           MeshPhongMaterial(color=RGBA{Float32}(1, 0, 0, 0.5)))

# using Makie
# sc = Scene()
# m = Polyhedra.Mesh(p3);
# mesh!(sc, m, color=:blue)
# display(sc)
# wireframe(m)


text:
module UnivariateGaussianProcess

using DocStringExtensions
using LinearAlgebra
using Optim

include("covariance_functions.jl")
include("acquisition_functions.jl")
include("plotting.jl")

export DATA_POINTS, GAUSS_PROC_DATA
export gaussian_process, func_est_μ, func_est_σ, get_max

# DATA STRUCTURES - - - - - - - -

"""
Data about previous evaluations of the objective function.

$(TYPEDFIELDS)
"""
struct DATA_POINTS
    "Points of the objective function domain at which the function has been evaluated."
    x::AbstractArray{<:Real}
    "The values of the objective function evaluations at these points of the domain."
    y::AbstractArray{<:Real}
end

"""
All data about the performed Bayesian Optimization.

$(TYPEDFIELDS)
"""
mutable struct GAUSS_PROC_DATA
    "The objective function.\n
    x::Real -> y::Real"
    obj_func::Function
    "The deviation of the objective function value when evaluated at a single point of the domain."
    noise::Real
    "The bounds of the searched interval of the objective function domain."
    bounds::Tuple{<:Real, <:Real}
    "Data about performed objective function evaluations."
    points::DATA_POINTS
    "The covariance function used by the Gaussian Process.\n
    x1::Real, x2::Real -> cov::Real"
    cov_func::Function
    "The kernel matrix of the Gaussian Process.\n
    K = C + noise^2 * I, where C[i,j] = cov_func(xi, xj)"
    K::Matrix{<:Real}
    "The inversed kernel matrix of the Gaussian Process."
    inv_K::Matrix{<:Real}
    "The used acquisition function.\n
    x::Real, data::GAUSS_PROC_DATA -> fitness::Real"
    acq_func::Function
    "The prior belief about the objective function.\n
    x::Real -> y::Real"
    init_mean::Function
end

# return data::GAUSS_PROC_DATA
function update_data!(
    new_point::Tuple{<:Real, <:Real},
    data::GAUSS_PROC_DATA,
)
    new_points = DATA_POINTS(vcat(data.points.x, new_point[1]), vcat(data.points.y, new_point[2]))
    new_K = update_kernel(data.cov_func, data.K, data.points.x, new_point[1]; data.noise)

    data.points = new_points
    data.K = new_K
    data.inv_K = inv(new_K)
    return data
end

# GAUSSIAN PROCESS - - - - - - - -

"""
    gaussian_process(
        obj_func::Function, # x::Real -> y::Real
        bounds::Tuple{<:Real, <:Real},
        cov_func::Function, # x1::Real, x2::Real -> cov::Real
        acq_func::Function; # x::Real, data::GAUSS_PROC_DATA -> fitness::Real
        noise::Real = 1e-8, # should be > 0
        init_mean::Function = (x) -> 0., # x::Real -> y::Real
        evals::Int = 1,
        points::Union{DATA_POINTS, Nothing} = nothing,
        plot::Bool = true,
        plot_only_final::Bool = false,
    )

Perform Bayesian Optimization with Gaussian Process surrogate model.

Return a structure containing all information about the performed optimization. (::GAUSS_PROC_DATA)

See also [`GAUSS_PROC_DATA`](@ref).

# Example

```julia-repl
julia> data = gaussian_process(sin, (-π, π), matern1, ucb; iters=20, noise=1e-4, init_mean = (x) -> x);

```

# Troubleshooting

If you get error about the kernel not being positive definite or sqrt being called on a negative value,
increase the objective function noise.

The objective function noise should not be set too close to zero
to prevent these issues arising because of numerical errors.
"""
function gaussian_process(
    obj_func::Function, # x::Real -> y::Real
    bounds::Tuple{<:Real, <:Real},
    cov_func::Function, # x1::Real, x2::Real -> cov::Real
    acq_func::Function; # x::Real, data::GAUSS_PROC_DATA -> fitness::Real
    noise::Real = 1e-4, # should be > 0
    init_mean::Function = (x) -> 0., # x::Real -> y::Real
    points::Union{DATA_POINTS, Nothing} = nothing,
    evals::Int = 1, # has to be >= 1
    plot::Bool = true,
    plot_only_final::Bool = false,
)
    evals <= 0 && throw(ArgumentError("`evals` has to be >= 1"))

    if points === nothing
        evals < 2 && throw(ArgumentError("`evals` has to be >= 2 if `points === nothing`"))
        points = DATA_POINTS(collect(bounds), obj_func.(collect(bounds)))
        evals -= 2
    end

    # init data
    init_K = create_kernel(cov_func, points.x; noise)
    data = GAUSS_PROC_DATA(
        obj_func,
        noise,
        bounds,
        points,
        cov_func,
        init_K,
        inv(init_K),
        acq_func,
        init_mean,
    )

    return gaussian_process(data; evals, plot, plot_only_final)
end

"""
    gaussian_process(
        data::GAUSS_PROC_DATA;
        evals::Int = 1,
        plot::Bool = true,
        plot_only_final::Bool = false,
    )

Continue in Bayesian Optimization with Gaussian Process surrogate model given the data from the previous run.

Return a structure containing all information about the performed optimization. (::GAUSS_PROC_DATA)

See also [`GAUSS_PROC_DATA`](@ref).

# Example

```julia-repl
julia> data = gaussian_process(sin, (-π, π), matern1, ucb; iters=20, noise=1e-4, init_mean = (x) -> x);

julia> data = gaussian_process(data; iters=20);

```

# Troubleshooting

If you get error about the kernel not being positive definite or sqrt being called on a negative value,
increase the objective function noise.

The objective function noise should not be set too close to zero
to prevent these issues arising because of numerical errors.
"""
function gaussian_process(
    data::GAUSS_PROC_DATA;
    evals::Int = 1, # has to be >= 0
    plot::Bool = true,
    plot_only_final::Bool = false,
)
    evals < 0 && throw(ArgumentError("`evals` has to be >= 1"))
    evals == 0 && return data

    for _ in 1:evals
        plot && (plot_only_final || plot_state(data))

        # choose new point to evaluate
        new_x = find_acq_opt(data)[1]

        # evaluate objective function
        new_y = data.obj_func(new_x)

        # augment data by the new point
        data = update_data!((new_x, new_y), data)
    end

    plot && plot_state(data)

    return data
end

# OBJECTIVE FUNCTION APPROXIMATION - - - - - - - -

"""
    func_est_μ(x::Real, data::GAUSS_PROC_DATA)

The most probable approximation of the objective function based on the Gaussian Process with the gathered data.

Return the mean of the gaussian distribution describing the objective function value estimated by the Gaussian Process at the given point.

See also [`func_est_σ`](@ref), [`gaussian_process`](@ref), [`GAUSS_PROC_DATA`](@ref).

# Examples

```julia-repl
julia> data = gaussian_process(sin, (-π, π), matern1, ucb; iters=20, init_mean = (x) -> x);

julia> func_est_μ(1, data)
0.8406095467907111
```

```julia-repl
julia> data = gaussian_process(sin, (-π, π), matern1, ucb; iters=20, init_mean = (x) -> x);

julia> approx_f = (x) -> func_est_μ(x, data);

julia> approx_f(1)
0.8406095467907111

julia> sin(1)
0.8414709848078965
```

(Note that in this example we use the UCB acquisition function
meaning that our goal is maximizing the objective function
and not finding its best approximation.)
"""
function func_est_μ(x::Real, data::GAUSS_PROC_DATA)
    k_x = data.cov_func.(x, data.points.x)
    return data.init_mean(x) + k_x' * data.inv_K * (data.points.y - data.init_mean.(data.points.x))
end

"""
    func_est_σ(x::Real, data::GAUSS_PROC_DATA)

The uncertainty of the objective function approximation.

Return the deviation of the gaussian distribution describing the objective function value estimated by the Gaussian Process at the given point.

See also [`func_est_μ`](@ref), [`gaussian_process`](@ref), [`GAUSS_PROC_DATA`](@ref).

# Examples

```julia-repl
julia> data = gaussian_process(sin, (-π, π), matern1, ucb; iters=20, init_mean = (x) -> x);

julia> func_est_σ(1, data)
0.20269550183507606
```

```julia-repl
julia> data = gaussian_process(sin, (-π, π), matern1, ucb; iters=20, init_mean = (x) -> x);

julia> uncertainty = (x) -> func_est_σ(x, data);

julia> uncertainty(1)
0.20269550183507606
```
"""
function func_est_σ(x::Real, data::GAUSS_PROC_DATA)
    K_x_x = data.cov_func(x, x)
    K_x = data.cov_func.(x, data.points.x)

    return sqrt(K_x_x - K_x' * data.inv_K * K_x)
end

"""
    get_max(data::GAUSS_PROC_DATA)

Return the best-so-far maximum found by the Bayesian optimization.

See also [`GAUSS_PROC_DATA`](@ref).

# Example

```julia-repl
julia> data = gaussian_process(sin, (-π, π), matern1, ucb; iters=20, init_mean = (x) -> x);

julia> get_max(data)
(1.575360217679531, 0.9999895854680737)
```
"""
function get_max(data::GAUSS_PROC_DATA)
    points = data.points
    max_i = argmax(points.y)
    return points.x[max_i], points.y[max_i]
end

# KERNEL COMPUTATION - - - - - - - -

# return K::Matrix{<:Real}
function create_kernel(
    cov_func::Function, # x1::Real, x2::Real -> cov::Real
    data_x::AbstractArray{<:Real};
    noise::Real,
)
    return [cov_func(x1, x2) for x1 in data_x, x2 in data_x] + noise^2 * I
end

# return K::Matrix{<:Real}
function update_kernel(
    cov_func::Function,
    K::Matrix{<:Real},
    old_data_x::AbstractArray{<:Real},
    new_x::Real;
    noise::Real,
)
    k_data_x = cov_func.(old_data_x, new_x)
    k_x_x = cov_func(new_x, new_x) + noise^2

    return vcat(hcat(K, k_data_x), vcat(k_data_x, k_x_x)')
end

# ACQ_FUNC OPTIMIZATION - - - - - - - -

# return opt_x::Real, opt_val::Real
function find_acq_opt(data::GAUSS_PROC_DATA)
    acq_func = x -> - data.acq_func(x[1], data)
    data_len = length(data.points.x)

    local_optima::Vector{Tuple{Float64, Float64}} = []
    sorted = sortperm(data.points.x)
    push!(local_optima, find_local_opt(acq_func, data.bounds[1], data.points.x[sorted[1]]))
    for i in 2:data_len
        push!(local_optima, find_local_opt(acq_func, data.points.x[sorted[i-1]], data.points.x[sorted[i]]))
    end
    push!(local_optima, find_local_opt(acq_func, data.points.x[sorted[data_len]], data.bounds[2]))

    optimum = reduce(local_optima) do a, b
        a[2] <= b[2] ? a : b
    end
    return optimum
end

# return opt_x::Real, opt_val::Real
function find_local_opt(f::Function, min_bound::Real, max_bound::Real)
    res = optimize(f, min_bound, max_bound)
    return Optim.minimizer(res)[1], Optim.minimum(res)
end

end


text:
immutable EnumerableJoin{T,TKey,TI,SO,SI,OKS,IKS,RS} <: Enumerable{T}
    outer::SO
    inner::SI
    outerKeySelector::OKS
    innerKeySelector::IKS
    resultSelector::RS
end

function join{TO,TI}(outer::Enumerable{TO}, inner::Enumerable{TI}, f_outerKeySelector::Function, outerKeySelector::Expr, f_innerKeySelector::Function, innerKeySelector::Expr, f_resultSelector::Function, resultSelector::Expr)
    TKeyOuter = Base.return_types(f_outerKeySelector, (TO,))[1]
    TKeyInner = Base.return_types(f_innerKeySelector, (TI,))[1]

    if TKeyOuter!=TKeyInner
        error("The keys in the join clause have different types, $TKeyOuter and $TKeyInner.")
    end

    SO = typeof(outer)
    SI = typeof(inner)

    T = Base.return_types(f_resultSelector, (TO,TI))[1]

    return EnumerableJoin{T,TKeyOuter,TI,SO,SI,FunctionWrapper{TKeyOuter,Tuple{TO}},FunctionWrapper{TKeyInner,Tuple{TI}},FunctionWrapper{T,Tuple{TO,TI}}}(outer,inner,f_outerKeySelector,f_innerKeySelector,f_resultSelector)
end

# TODO This should be changed to a lazy implementation
function start{T,TKeyOuter,TI,SO,SI,OKS,IKS,RS}(iter::EnumerableJoin{T,TKeyOuter,TI,SO,SI,OKS,IKS,RS})
    results = Array(T,0)

    inner_dict = OrderedDict{TKeyOuter,Array{TI,1}}()
    for i in iter.inner
        key = iter.innerKeySelector(i)
        if !haskey(inner_dict, key)
            inner_dict[key] = Array(TI,0)
        end
        push!(inner_dict[key], i)
    end

    for i in iter.outer
        outerKey = iter.outerKeySelector(i)
        if haskey(inner_dict,outerKey)
            for j in inner_dict[outerKey]
                push!(results, iter.resultSelector(i,j))
            end
        end
    end

    return results,1
end

function next{T,TKeyOuter,TI,SO,SI,OKS,IKS,RS}(iter::EnumerableJoin{T,TKeyOuter,TI,SO,SI,OKS,IKS,RS},state)
    results = state[1]
    curr_index = state[2]
    return results[curr_index], (results, curr_index+1)
end

function done{T,TKeyOuter,TI,SO,SI,OKS,IKS,RS}(iter::EnumerableJoin{T,TKeyOuter,TI,SO,SI,OKS,IKS,RS},state)
    results = state[1]
    curr_index = state[2]
    return curr_index > length(results)
end


text:
using Dates
using LazyArtifacts
using PkgUtility
using Test

# function and struct used to test
f(x) = x^2;

struct TestStruct
    a::Any
    b::Any
end

@testset "PkgUtility --- Artifact" begin
    # deploy an artifact per file
    mkdir("temp");
    cp("example.toml", "temp/test_1.txt"; force=true);
    cp("example.toml", "temp/test_2.txt"; force=true);
    deploy_artifact!("temp/TempArtifacts.toml", "example_artifact", "$(pwd())/temp", "$(pwd())", ["url_1", "url_2"]);
    deploy_artifact!("temp/TempArtifacts.toml", "example_artifact", "$(pwd())/temp", "$(pwd())", ["url_1", "url_2"]);

    # remove temp files
    meta = artifact_meta("example_artifact", "temp/TempArtifacts.toml");
    hash = meta["git-tree-sha1"];
    rm("$(homedir())/.julia/artifacts/$(hash)"; recursive=true);
    rm("temp"; recursive=true);
    rm("example_artifact.tar.gz");
    @test true;

    # deploy an artifact per folder
    mkdir("temp");
    cp("example.toml", "temp/test_1.txt"; force=true);
    cp("example.toml", "temp/test_2.txt"; force=true);
    deploy_artifact!("temp/TempArtifacts.toml", "example_artifact", "$(pwd())/temp", ["test_1.txt", "test_2.txt"], "$(pwd())", ["url_1", "url_2"]);

    # remove temp files
    meta = artifact_meta("example_artifact", "temp/TempArtifacts.toml");
    hash = meta["git-tree-sha1"];
    rm("$(homedir())/.julia/artifacts/$(hash)"; recursive=true);
    rm("temp"; recursive=true);
    rm("example_artifact.tar.gz");
    @test true;

    # predownload the artifact directly from the given URL
    # the manual steps are only tested on Windows
    predownload_artifact!("CI_PFT_2X_1Y_V1", "example.toml");
    predownload_artifact!("CI_PFT_2X_1Y_V1", "example.toml");
    @test true;
end

println();
@testset "PkgUtility --- Date" begin
    @test parse_timestamp("20000401"; out_format="DOY") == 92;
    @test parse_timestamp("20010401"; out_format="DOY") == 91;
    @test parse_timestamp("20010401000000"; in_format="YYYYMMDDhhmmss", out_format="FDOY") == 91;
    @test parse_timestamp("20010401"; out_format="DATE") == Date("2001-04-01");
    @test parse_timestamp("20010401"; out_format="DATETIME") == DateTime("2001-04-01T00:00:00");
    @test parse_timestamp(2000, 100) == "20000409";
    @test parse_timestamp(2001, 100) == "20010410";
    @test month_days(2019, 2) == 28;
    @test month_days(2020, 2) == 29;
    @test month_ind(2019, 60) == 3;
    @test month_ind(2020, 60) == 2;

    # some time labeled information
    err_info = terror("This is an error!");
    err_info = tinfo("This is an info!");
    err_info = twarn("This is a warning!");
    @test true;
end

println();
@testset "PkgUtility --- Display" begin
    xxx = [
        "a"    => "a",
        "asa"  => "asdasda",
        "haha" => [
                   "a"    => "a",
                   "asa"  => "asdasda",
                   "haha" => [
                              1   => "a",
                              "a" => "asdasda",
                             ],
                  ],
    ];
    @info tinfo("Display the dict in a pretty way:");
    pretty_display!(xxx);
    @test true;
end

println();
@testset "PkgUtility --- Math" begin
    # test integral function
    numerical∫(rand(5), rand(5));
    @info tinfo("Expecting dimension mismatch warning here...");
    numerical∫(rand(5), rand(6));
    numerical∫(rand(5), 0.1);
    @test true;

    # test quadratic solvers
    for FT in [Float32, Float64]
        @test lower_quadratic(FT( 1), FT(-3), FT( 2)) == FT(1);
        @test lower_quadratic(FT(-1), FT( 3), FT(-2)) == FT(1);
        @test isnan(lower_quadratic(FT( 1), FT(-3), FT(10)));
        @test upper_quadratic(FT( 1), FT(-3), FT( 2)) == FT(2);
        @test upper_quadratic(FT(-1), FT( 3), FT(-2)) == FT(2);
        @test isnan(upper_quadratic(FT( 1), FT(-3), FT(10)));
    end

    # test statistics
    xx = rand(10); xx[1]=NaN;
    yy = rand(10);
    nanmax(xx);
    nanmean(xx);
    nanmedian(xx);
    nanmin(xx);
    nanpercentile(xx, 50);
    nanstd(xx);
    mae(xx, yy);
    mape(xx, yy);
    mase(xx, yy);
    rmse(xx, yy);
    @test true;
end

println();
@testset "PkgUtility --- IO" begin
    # test the dataframe
    df = dataframe();
    df = dataframe([rand(3), rand(3)], ["A","B"]);
    @test true;

    # test the read_nc
    path = artifact_path(artifact_hash("CI_PFT_2X_1Y_V1", "example.toml"));
    file = path * "/CI_PFT_2X_1Y_V1.nc";
    data = read_nc(Float32, file, "clump");
    @test true;
end

println();
@testset "PkgUtility --- FT test" begin
    for FT in [Float32, Float64]
        sa = TestStruct(ones(FT,5), 2);
        @test FT_test(ones(FT,5), FT);
        @test FT_test(f, FT);
        @test FT_test(Test, FT);
        @test FT_test(FT, FT);
        @test FT_test(1, FT);
        @test FT_test(FT(1), FT);
        @test FT_test("a", FT);
        @test FT_test(sa, FT);
    end
    @test FT_test([1f0,2f0], Float64) == false;
    @test FT_test([[1,2],2f0], Float64) == false;
end

println();
@testset "PkgUtility --- NaN test" begin
    for FT in [Float32, Float64]
        sa = TestStruct(ones(FT,5), 2);
        sb = TestStruct(ones(FT,5), NaN);
        @test NaN_test(ones(FT,5));
        @test NaN_test(f);
        @test NaN_test(Test);
        @test NaN_test(FT);
        @test NaN_test(1);
        @test NaN_test("a");
        @test NaN_test(sa);
        @test NaN_test(NaN) == false;
        @test NaN_test([1,2,NaN]) == false;
        @test NaN_test([[1,2,NaN], 2]) == false;
        @test NaN_test(sb) == false;
    end
end

println();
@testset "PkgUtility --- Deprecation warnings" begin
    # deprecated function
    @info tinfo("Expecting deprecation warnings here");
    @test true;
end


text:
## function to sort the SKUs after the number of sales
function SORTSALES(trans::SparseMatrixCSC{Float64, Int64})
    sales = Array{Int64,2}(undef,size(trans,2),4) .= 0
    sales[:,3] = sum(trans,dims=1)
    sales[:,4] .= 0
    for i = 1:size(trans,2)
        sales[i,1] = i
    end
    sales = sortslices(sales,dims=1,by=x->x[3],rev=true)
    for i = 1:size(trans,2)
        sales[i,2] = i
    end
    # 1 column: index of SKU in allocation matrix
    # 2 column: index of SKU in sorted sales
    # 3 column: number of total sold SKU
    # 4 column: indicator whether SKU is allocated
    return sales::Array{Int64,2}
end

## function to sort the SKU pairs after the highest number of coappearances
function SORTPAIRS(Q::Array{Int64,2})
    pairs = Array{Int64,2}(undef,round(Int64,(size(Q,1)*size(Q,1)-size(Q,1))/2),3)
    iter = 1
    for i = 2:size(Q,1)
        for j = 1 : i-1
            pairs[iter,1] = i
            pairs[iter,2] = j
            pairs[iter,3] = Q[i,j]
            iter += 1
        end
    end
    #pairs = sortslices(pairs,dims=1,by=x->x[3],rev=true)
    pairs = pairs[sortperm(pairs[:,3],rev=true),:]
    return pairs::Array{Int64,2}
end

## function to check the number of coappearances in the selected warehouse
## for the selected product
function CURRENTCOAPP(X::Array{Int64,2},
                      warehouse::Int64,
                      sku::Int64,
                      Q::Array{Int64,2})
    sum(X[:,warehouse] .* Q[sku,:])
end

## function to start place the seeds in the greedy seeds heuristic
function GREEDYSEEDSTART!(sales::Array{Int64,2},
                          X::Array{Int64,2},
                          Q::Array{Int64,2},
                          capacity_left::Array{Int64,1})
    for k = 2:size(capacity_left,1)
        # list all SKUs that have not been allocated yet and
        top_sales = sales[sales[:,4] .== 0,:]
        # discard from the list those SKUs that are not in the top decile
        top_sales = top_sales[1:round(Int64, size(top_sales,1)/10),:]
        # find the SKU from the list with the least coappearances to the allocated SKUs
        coapp = Array{Int64,1}(undef,size(top_sales,1)) .= 0
        for i in 1:size(top_sales,1)
            for j in 1:size(capacity_left,1)-1
                coapp[i] = CURRENTCOAPP(X,j,i,Q)
            end
        end
        # allocate the SKU to the DC
        sales[top_sales[findmin(coapp)[2],2],4] = 1
        X[top_sales[findmin(coapp)[2],1],k] = 1
    end
    return sales::Array{Int64,2},
           X::Array{Int64,2}
end

## function to allocate the SKU with the highest potential
function GREEDYSEEDMAIN!(sales::Array{Int64,2},
                         X::Array{Int64,2},
                         Q::Array{Int64,2},
                         capacity_left::Array{Int64,1})
    for i = 1:size(Q,1)
        if sales[i,4] == 0
            best_allocation = Array{Int64,1}(undef,size(capacity_left,1)) .= 0
            for d = 1:size(capacity_left,1)
                if capacity_left[d] > 1
                    # retrieve the list of SKUs allocated to d so far
                    # calculate the average of coapperances with the SKUs in d
                    best_allocation[d] = CURRENTCOAPP(X,d,i,Q)
                end
            end
            if findmax(best_allocation)[1] > 0
                X[sales[i,1],findmax(best_allocation)[2]] = 1
                capacity_left[findmax(best_allocation)[2]] -= 1
            else
                X[sales[i,1],findmax(capacity_left)[2]] = 1
                capacity_left[findmax(capacity_left)[2]] -= 1
            end
            sales[i,4] = 1
        end
    end
    return sales::Array{Int64,2},
           X::Array{Int64,2},
           capacity_left::Array{Int64,1}
end

## function to allocate the SKUs in the greedy pairs heuristic
function GREEDYPAIRSMAIN!(pairs::Array{Int64,2},
                          X::Array{Int64,2},
                          capacity_left::Array{Int64,1})
    ## Sort the DCs by decreasing capacity
    for i = 1:size(pairs,1)
        ## For each pair of SKUs in sorted list of pairs do
        for j = 1:2
            for d in 1:size(capacity_left,1)
                if capacity_left[d] > 0
                    ## if SKU n is not allocated to DC j then allocate SKU s to DC d
                    if sum(X[pairs[i,j],:]) == 0
                        X[pairs[i,j],d] = 1
                        capacity_left[d] -= 1
                    end
                end
            end
        end
        if sum(X) == size(X,1)
            break
        end
    end
    return X::Array{Int64,2},
           capacity_left::Array{Int64,1}
end

## function to fill up the warehouses with the potentially best SKUs
## after each SKU has already been allocated once
function FILLUP(X::Array{Int64,2},
                Q::Array{Int64,2},
                capacity_left::Array{Int64,1})
    for d = 1:size(capacity_left,1)
        while capacity_left[d] > 0
            best_allocation = Array{Int64,1}(undef,size(Q,1)) .= 0
            for i = 1:size(Q,1)
                if X[i,d] == 0
                    a = @view X[:,d]
                    b = @view Q[:,i]
                    best_allocation[i] = dot(a,b)
                end
            end
            if findmax(best_allocation)[1] > 0
                X[findmax(best_allocation)[2],d] = 1
                capacity_left[d] -= 1
            else
                capacity_left[d] = 0
            end
        end
    end
    return X::Array{Int64,2}
end

## function to calculate the number of SKUs that can be allocated
## multiple times
function BESTSELLING_B(skus::Int64, capacity_left::Array{Int64,1})
    floor(Int64,(sum(capacity_left)-skus)/(size(capacity_left,1)-1))
end

## function to initialise the bestselling heuristic
function BESTSELLINGSTART!(sales::Array{Int64,2},
                           capacity_left::Array{Int64,1},
                           B::Int64,
                           X::Array{Int64,2})
    for d = 1:size(capacity_left,1)
        if capacity_left[d] < B
            for i = 1:capacity_left[d]
                X[sales[i,1],d] = 1
                sales[sales[i,2],4] = 1
                capacity_left[d] -= 1
            end
        end
        B = BESTSELLING_B(size(X,1),capacity_left)
    end
    return sales::Array{Int64,2},
           capacity_left::Array{Int64,1},
           B::Int64,
           X::Array{Int64,2}
end

## function to allocate the SKUs that have been sold most in each
## warehouse until the remaining SKUs can only be allocated once
function BESTSELLINGTOP!(sales::Array{Int64,2},
                         capacity_left::Array{Int64,1},
                         B::Int64,
                         X::Array{Int64,2})
    for d = 1:size(capacity_left,1)
        if capacity_left[d] >= B
            for i = 1:B
                X[sales[i,1],d] = 1
                sales[sales[i,2],4] = 1
                capacity_left[d] -= 1
            end
        end
    end
    return sales::Array{Int64,2},
           capacity_left::Array{Int64,1},
           X::Array{Int64,2}
end

max_stars_repo_name: beyondatlas/splitdeliveries

max_issues_repo_name: beyondatlas/splitdeliveries

max_forks_repo_name: beyondatlas/splitdeliveries


text:
module NCP

include("NCP_Exact.jl")
include("NCP_EM.jl")
include("NCP_VB.jl")
include("NCP_MC.jl")
include("BackwardKernel.jl")

using .NCP_Exact, .NCP_EM, .NCP_VB, .NCP_MC, .BackwardKernel

export allocations, full, generate
export log_marginal, alloc_dist, EP_dist
export log_posterior_T, log_PT
export standard_EM
export standard_VB
export Particle, smc_weight, filtering_proposal, resample
export EventQueue, sum, length, eltype, iterate

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
function getStripPack(num, widthMax, heightMax, wCap, seed)
    Random.seed!(seed)
    vec_w = rand(1:widthMax, num)  # [width of items]
    vec_h = rand(1:heightMax, num)  # [height of items]
    println("Whether wCap >= widthMax? $(wCap >= widthMax)")
    return (wCap, vec_w, vec_h)
end

function solveStripPack(wCap, vec_w, vec_h)
    n = length(vec_w)
    hCap = sum(vec_h)
    println("################################################################################") ########################
    println("**** hCap = $(hCap)")
    timeStart = time()
    IP = Model(solver = GurobiSolver())
    # IP = Model(solver = CplexSolver())
    @variable(IP, 0 <= x[i = 1:n] <= wCap - vec_w[i], Int )
    @variable(IP, 0 <= y[i = 1:n] <= hCap - vec_h[i], Int )
    @variable(IP, matAlpha[i=1:n, j = 1:n], Bin )
    @variable(IP, matBeta[i = 1:n, j = 1:n], Bin )
    @variable(IP, 0 <= z <= hCap, Int )
    @objective(IP, Min, z)
    @constraint(IP, [i = 1:n, j = 1:n; i < j], matAlpha[i,j] + matAlpha[j,i] + matBeta[i,j] + matBeta[j,i] >= 1)
    @constraint(IP, [i = 1:n], y[i] + vec_h[i] <= z)
    # given by variable bounds:
    # @constraint(IP, [i=1:n], x[i] + vec_w[i] <= wCap)
    @constraint(IP, [i = 1:n, j = 1:n; i != j], x[i] + vec_w[i] <= x[j] + wCap * (1 - matAlpha[i,j]))
    @constraint(IP, [i = 1:n, j = 1:n; i != j], y[i] + vec_h[i] <= y[j] + hCap * (1 - matBeta[i,j]))
    # writeLP(IP,"StripPacking.lp",genericnames=false)
    status = solve(IP, relaxation=true)
    println("**** Objective value, relaxation: ", getobjectivevalue(IP))
    #println("y = ", getvalue(y))
    status = solve(IP)
    println("**** Objective value IP: ", getobjectivevalue(IP))
    println("**** x = ", getvalue(x))
    println("**** y = ", getvalue(y))
    println("################################################################################\n",
            "**** Elapsed time is $(time() - timeStart) seconds.\n",
            "#### End #######################################################################") ########################
end

function getSeriesMatrix(i::Int64, j::Int64, num::Int64 = num)
    if i == j
        throw("!!!! Error: i = j.")
    elseif i < j
        s = (i - 1) * (num - 1) + j - 1
    elseif i > j
        s = (i - 1) * (num - 1) + j
    end
    return s
end

function getCoefMatrix(wCap, vec_w, vec_h)
    num = length(vec_w)
    hCap = sum(vec_h)
    vecNumCol = Array{Int64}(undef, 5)
    vecNumCol[1:2] .= convert(Int64, num^2 - num)
    vecNumCol[3:4] .= convert(Int64, num)
    vecNumCol[5] = convert(Int64, 1)
    vecNumRow = Array{Int64}(undef, 5)
    vecNumRow[1] = convert(Int64, sum(collect(1:(num - 1))))
    vecNumRow[2] = convert(Int64, num)
    vecNumRow[3] = convert(Int64, num)
    vecNumRow[4] = convert(Int64, num^2 - num)
    vecNumRow[5] = convert(Int64, num^2 - num)
    ##
    vec_c = zeros(convert(Int64, sum(vecNumCol)))
    vec_c[convert(Int64, sum(vecNumCol))] = -1
    ## mat_a
    mat_a11 = zeros(vecNumRow[1], vecNumCol[1])
    whiRow = 1
    for i = 1:num
        for j = 1:num
            if i < j
                s1 = getSeriesMatrix(i, j, num)
                mat_a11[whiRow, s1] = -1
                s2 = getSeriesMatrix(j, i, num)
                mat_a11[whiRow, s2] = -1
                whiRow += 1
            end
        end
    end
    mat_a1 = [mat_a11 mat_a11 zeros(vecNumRow[1], sum(vecNumCol[3:5]))]
    mat_a2 = [zeros(num, sum(vecNumCol[1:3])) Matrix{Int64}(I, num, num) -ones(num, 1)]
    mat_a3 = [zeros(num, sum(vecNumCol[1:2])) Matrix{Int64}(I, num, num) zeros(num, sum(vecNumCol[4:5]))]
    mat_a41 = zeros(vecNumRow[4], vecNumCol[3])
    whiRow = 1
    for i = 1:num
        for j = 1:num
            if i != j
                mat_a41[whiRow, i] = 1
                mat_a41[whiRow, j] = -1
                whiRow += 1
            end
        end
    end
    mat_a4 = hcat(wCap * Matrix{Int64}(I, vecNumCol[1], vecNumCol[1]), zeros(vecNumRow[4], vecNumCol[2]),
        mat_a41, zeros(vecNumRow[4], sum(vecNumCol[4:5])))
    mat_a5 = hcat(zeros(vecNumRow[5], vecNumCol[1]), hCap * Matrix{Int64}(I, vecNumCol[1], vecNumCol[1]),
        zeros(vecNumRow[5], vecNumCol[3]), mat_a41, zeros(vecNumRow[5], vecNumCol[5]))
    mat_a = vcat(mat_a1, mat_a2, mat_a3, mat_a4, mat_a5)
    ##
    vec_b = hcat(zeros(sum(vecNumRow), 1))
    vec_b[1:vecNumRow[1], 1] .= -1
    vec_b[(vecNumRow[1] + 1):sum(vecNumRow[1:2]), 1] .= -vec_h
    vec_b[(sum(vecNumRow[1:2]) + 1):sum(vecNumRow[1:3]), 1] = wCap .- vec_w
    for i = 1:num
        vec_b[(sum(vecNumRow[1:3]) + (i-1) * (num - 1) + 1):(sum(vecNumRow[1:3]) + i * (num - 1)), 1] .= wCap - vec_w[i]
    end
    for i = 1:num
        vec_b[(sum(vecNumRow[1:4]) + (i-1) * (num - 1) + 1):(sum(vecNumRow[1:4]) + i * (num - 1)), 1] .= hCap - vec_h[i]
    end
    ## Parameters for Sub-problems
    vecVecIndexSub = [vcat(collect(1:vecNumCol[1]), collect((sum(vecNumCol[1:2]) + 1):sum(vecNumCol[1:3]))),
        vcat(collect(sum(vecNumCol[1] + 1):sum(vecNumCol[1:2])),
        collect((sum(vecNumCol[1:3]) + 1):sum(vecNumCol[1:5])))]
    numSub = 2
    m_mat_h = sum(vecNumRow[1:3])
    vecRowMatD = [vecNumRow[4], vecNumRow[5]]
    vecVecIndexBinInSub = [collect(1:vecNumCol[1]), collect(1:vecNumCol[2])]
    # vecVecIndexSub = [collect(1:sum(vecNumCol[1:5]))]
    # numSub = 1
    # m_mat_h = sum(vecNumRow[1:3])
    # vecRowMatD = [sum(vecNumRow[4:5])]
    # vecVecIndexBinInSub = [collect(1:sum(vecNumCol[1:2]))]
    return (vec_c, mat_a, vec_b, vecVecIndexSub, numSub, m_mat_h, vecRowMatD, vecVecIndexBinInSub)
end

function getData()
    ## Test the algorihtm using SSCFLP
    vec_c = - [292  453  359  219  268  736  291  443  403  498  400  967]
    mat_a = [
        Matrix{Float64}(I, 5, 5) zeros(5, 1) Matrix{Float64}(I, 5, 5) zeros(5, 1);
        14 20 6 16 10 -43 zeros(1, 6);
        zeros(1, 6) 14 20 6 16 10 -43
        ]
    vec_b = hcat([ones(5,1); 0; 0])
    vecVecIndexSub = [collect(1:6), collect(7:12)]
    numSub = 2
    m_mat_h = 5
    vecRowMatD = [1, 1]
    vecVecIndexBinInSub = [collect(1:6), collect(1:6)]
    return (vec_c, mat_a, vec_b, vecVecIndexSub, numSub, m_mat_h, vecRowMatD, vecVecIndexBinInSub)
end


text:
using MLStyle
using Pipe
using Statistics

const data = readlines()

const char_pair = Dict([
    '(' => ')',
    '[' => ']',
    '{' => '}',
    '<' => '>',
])

const char_score = Dict([
    ')' => 3,
    ']' => 57,
    '}' => 1197,
    '>' => 25137,
    '(' => 1,
    '[' => 2,
    '{' => 3,
    '<' => 4,
])

function first_err_char(l)
    stack = []
    for c in l
        @match c begin
            cc && if cc ∈ keys(char_pair) end => push!(stack, cc)
            cc && if char_pair[stack[end]] == cc end => pop!(stack)
            cc => return cc, nothing
        end
    end
    return nothing, reverse(stack)
end

# Part One
@pipe data .|> first_err_char .|> get(char_score, _[1], 0) |> sum |> println

compute_score_completion(::Nothing) = 0
compute_score_completion(chars) = mapreduce((acc, x) -> 5acc + x,
                                            chars;
                                            init=0) do c
    get(char_score, c, nothing)
end

# Part Two
@pipe data .|>
      first_err_char .|>
      compute_score_completion(_[2]) |>
      filter(>(0), _) |>
      median |>
      convert(Int, _) |>
      println

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
import ThresholdStability: unstatedep, unconstrained

Σ = [randn(2, 2) for _ in 1:4]
G = LightAutomaton(4); add_transition!(G,1,1,1); add_transition!(G,1,2,1); add_transition!(G,2,3,2); add_transition!(G, 3, 4, 3); add_transition!(G, 4, 1, 4)
X = [1,2,3,4]
s = discreteswitchedsystem(Σ, G, X)

@testset "Switched system type conversion" begin
    t = unstatedep(s)
    @test typeof(t) <: ConstrainedDiscreteSwitchedLinearSystem
    @test typeof(unconstrained(s)) <: StateDepDiscreteSwitchedLinearSystem
    @test typeof(unconstrained(t)) <: DiscreteSwitchedLinearSystem
end


text:
# ---
# title: 799. Champagne Tower
# id: problem799
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/champagne-tower/description/>
# hidden: true
# ---
#
# We stack glasses in a pyramid, where the **first** row has `1` glass, the
# **second** row has `2` glasses, and so on until the 100th row.  Each glass
# holds one cup of champagne.
#
# Then, some champagne is poured into the first glass at the top.  When the
# topmost glass is full, any excess liquid poured will fall equally to the glass
# immediately to the left and right of it.  When those glasses become full, any
# excess champagne will fall equally to the left and right of those glasses, and
# so on.  (A glass at the bottom row has its excess champagne fall on the
# floor.)
#
# For example, after one cup of champagne is poured, the top most glass is full.
# After two cups of champagne are poured, the two glasses on the second row are
# half full.  After three cups of champagne are poured, those two cups become
# full - there are 3 full glasses total now.  After four cups of champagne are
# poured, the third row has the middle glass half full, and the two outside
# glasses are a quarter full, as pictured below.
#
# ![](https://s3-lc-upload.s3.amazonaws.com/uploads/2018/03/09/tower.png)
#
# Now after pouring some non-negative integer cups of champagne, return how full
# the `jth` glass in the `ith` row is (both `i` and `j` are 0-indexed.)
#
#
#
# **Example 1:**
#
#
#
#     Input: poured = 1, query_row = 1, query_glass = 1
#     Output: 0.00000
#     Explanation: We poured 1 cup of champange to the top glass of the tower (which is indexed as (0, 0)). There will be no excess liquid so all the glasses under the top glass will remain empty.
#
#
# **Example 2:**
#
#
#
#     Input: poured = 2, query_row = 1, query_glass = 1
#     Output: 0.50000
#     Explanation: We poured 2 cups of champange to the top glass of the tower (which is indexed as (0, 0)). There is one cup of excess liquid. The glass indexed as (1, 0) and the glass indexed as (1, 1) will share the excess liquid equally, and each will get half cup of champange.
#
#
# **Example 3:**
#
#
#
#     Input: poured = 100000009, query_row = 33, query_glass = 17
#     Output: 1.00000
#
#
#
#
# **Constraints:**
#
#   * `0 <= poured <= 109`
#   * `0 <= query_glass <= query_row < 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
    basicblock(inplanes, outplanes, downsample = false)

Create a basic residual block
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments:
- `inplanes`: the number of input feature maps
- `outplanes`: a list of the number of output feature maps for each convolution
               within the residual block
- `downsample`: set to `true` to downsample the input
"""
function basicblock(inplanes, outplanes, downsample = false)
  stride = downsample ? 2 : 1
  Chain(conv_bn((3, 3), inplanes, outplanes[1]; stride = stride, pad = 1, bias = false)...,
        conv_bn((3, 3), outplanes[1], outplanes[2], identity; stride = 1, pad = 1, bias = false)...)
end

"""
    bottleneck(inplanes, outplanes, downsample = false)

Create a bottleneck residual block
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments:
- `inplanes`: the number of input feature maps
- `outplanes`: a list of the number of output feature maps for each convolution
               within the residual block
- `downsample`: set to `true` to downsample the input
"""
function bottleneck(inplanes, outplanes, downsample = false)
  stride = downsample ? 2 : 1
  Chain(conv_bn((1, 1), inplanes, outplanes[1]; stride = 1, bias = false)...,
        conv_bn((3, 3), outplanes[1], outplanes[2]; stride = stride, pad = 1, bias = false)...,
        conv_bn((1, 1), outplanes[2], outplanes[3], identity; stride = 1, bias = false)...)
end

"""
    resnet(block, residuals::NTuple{2, Any}, connection = addrelu;
           channel_config, block_config, nclasses = 1000)

Create a ResNet model
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments
- `block`: a function with input `(inplanes, outplanes, downsample=false)` that returns
           a new residual block (see [`Metalhead.basicblock`](#) and [`Metalhead.bottleneck`](#))
- `residuals`: a 2-tuple of functions with input `(inplanes, outplanes, downsample=false)`,
               each of which will return a function that will be used as a new "skip" path to match a residual block.
              [`Metalhead.skip_identity`](#) and [`Metalhead.skip_projection`](#) can be used here.
- `connection`: the binary function applied to the output of residual and skip paths in a block
- `channel_config`: the growth rate of the output feature maps within a residual block
- `block_config`: a list of the number of residual blocks at each stage
- `nclasses`: the number of output classes
"""
function resnet(block, residuals::NTuple{2, Any}, connection = addrelu;
                channel_config, block_config, nclasses = 1000)
  inplanes = 64
  baseplanes = 64
  layers = []
  append!(layers, conv_bn((7, 7), 3, inplanes; stride = 2, pad = 3, bias = false))
  push!(layers, MaxPool((3, 3), stride = (2, 2), pad = (1, 1)))
  for (i, nrepeats) in enumerate(block_config)
    # output planes within a block
    outplanes = baseplanes .* channel_config
    # push first skip connection on using first residual
    # downsample the residual path if this is the first repetition of a block
    push!(layers, Parallel(connection, block(inplanes, outplanes, i != 1),
                                       residuals[1](inplanes, outplanes[end], i != 1)))
    # push remaining skip connections on using second residual
    inplanes = outplanes[end]
    for _ in 2:nrepeats
      push!(layers, Parallel(connection, block(inplanes, outplanes, false),
                                         residuals[2](inplanes, outplanes[end], false)))
      inplanes = outplanes[end]
    end
    # next set of output plane base is doubled
    baseplanes *= 2
  end

  return Chain(Chain(layers...),
               Chain(AdaptiveMeanPool((1, 1)), flatten, Dense(inplanes, nclasses)))
end

"""
    resnet(block, shortcut_config::Symbol, connection = addrelu;
           channel_config, block_config, nclasses = 1000)

Create a ResNet model
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments
- `block`: a function with input `(inplanes, outplanes, downsample=false)` that returns
           a new residual block (see [`Metalhead.basicblock`](#) and [`Metalhead.bottleneck`](#))
- `shortcut_config`: the type of shortcut style (either `:A`, `:B`, or `:C`)
    - `:A`: uses a [`Metalhead.skip_identity`](#) for all residual blocks
    - `:B`: uses a [`Metalhead.skip_projection`](#) for the first residual block
            and [`Metalhead.skip_identity`](@) for the remaining residual blocks
    - `:C`: uses a [`Metalhead.skip_projection`](#) for all residual blocks
- `connection`: the binary function applied to the output of residual and skip paths in a block
- `channel_config`: the growth rate of the output feature maps within a residual block
- `block_config`: a list of the number of residual blocks at each stage
- `nclasses`: the number of output classes
"""
function resnet(block, shortcut_config::Symbol, args...; kwargs...)
  shortcut = if shortcut_config == :A
      (skip_identity, skip_identity)
    elseif shortcut_config == :B
      (skip_projection, skip_identity)
    elseif shortcut_config == :C
      (skip_projection, skip_projection)
    else
      error("Unrecognized shortcut_config ($shortcut_config) passed to `resnet` (use :A, :B, or :C).")
  end
  resnet(block, shortcut, args...; kwargs...)
end

const resnet_config =
  Dict(18 => (([1, 1], [2, 2, 2, 2], :A), basicblock),
       34 => (([1, 1], [3, 4, 6, 3], :A), basicblock),
       50 => (([1, 1, 4], [3, 4, 6, 3], :B), bottleneck),
       101 => (([1, 1, 4], [3, 4, 23, 3], :B), bottleneck),
       152 => (([1, 1, 4], [3, 8, 36, 3], :B), bottleneck))

"""
    ResNet(channel_config, block_config, shortcut_config;
           block, connection = addrelu, nclasses = 1000)

Create a `ResNet` model
([reference](https://arxiv.org/abs/1512.03385v1)).
See also [`resnet`](#).

# Arguments
- `channel_config`: the growth rate of the output feature maps within a residual block
- `block_config`: a list of the number of residual blocks at each stage
- `shortcut_config`: the type of shortcut style (either `:A`, `:B`, or `:C`)
- `block`: a function with input `(inplanes, outplanes, downsample=false)` that returns
           a new residual block (see [`Metalhead.basicblock`](#) and [`Metalhead.bottleneck`](#))
- `connection`: the binary function applied to the output of residual and skip paths in a block
- `nclasses`: the number of output classes
"""
struct ResNet
  layers
end

function ResNet(channel_config, block_config, shortcut_config;
                block, connection = addrelu, nclasses = 1000)
  layers = resnet(block,
                  shortcut_config,
                  connection;
                  channel_config = channel_config,
                  block_config = block_config,
                  nclasses = nclasses)

  ResNet(layers)
end

@functor ResNet

(m::ResNet)(x) = m.layers(x)

backbone(m::ResNet) = m.layers[1]
classifier(m::ResNet) = m.layers[2]

"""
    ResNet(depth = 50; pretrain = false, nclasses = 1000)

Create a ResNet model with a specified depth
([reference](https://arxiv.org/abs/1512.03385v1)).
See also [`Metalhead.resnet`](#).

# Arguments
- `depth`: depth of the ResNet model. Options include (18, 34, 50, 101, 152).
- `nclasses`: the number of output classes

!!! warning
    Only `ResNet(50)` currently supports pretrained weights.
"""
function ResNet(depth::Int = 50; pretrain = false, nclasses = 1000)
    @assert depth in keys(resnet_config) "`depth` must be one of $(sort(collect(keys(resnet_config))))"

    config, block = resnet_config[depth]
    model = ResNet(config...; block = block, nclasses = nclasses)
    pretrain && loadpretrain!(model, string("ResNet", depth))
    model
end

# Compat with Methalhead 0.6; remove in 0.7
@deprecate ResNet18(; kw...) ResNet(18; kw...)
@deprecate ResNet34(; kw...) ResNet(34; kw...)
@deprecate ResNet50(; kw...) ResNet(50; kw...)
@deprecate ResNet101(; kw...) ResNet(101; kw...)
@deprecate ResNet152(; kw...) ResNet(152; kw...)


text:
using Oceananigans: instantiated_location
using Oceananigans.Architectures: arch_array, device_event, device_copy_to!
using Oceananigans.Operators: assumed_field_location
using Oceananigans.Fields: reduced_dimensions

using Oceananigans.BoundaryConditions:
            ContinuousBoundaryFunction,
            DiscreteBoundaryFunction,
            permute_boundary_conditions,
            fill_halo_event!,
            CBCT,
            CBC

import Oceananigans.Fields: tupled_fill_halo_regions!, boundary_conditions, data

import Oceananigans.BoundaryConditions:
            fill_halo_regions!,
            fill_west_and_east_halo!,
            fill_south_and_north_halo!,
            fill_bottom_and_top_halo!,
            fill_west_halo!,
            fill_east_halo!,
            fill_south_halo!,
            fill_north_halo!

@inline bc_str(::MultiRegionObject) = "MultiRegion Boundary Conditions"

@inline extract_field_buffers(field::Field)          = field.boundary_buffers
@inline boundary_conditions(field::MultiRegionField) = field.boundary_conditions

# This can be implemented once we have a buffer for field_tuples
@inline function tupled_fill_halo_regions!(full_fields, grid::MultiRegionGrid, args...; kwargs...)
    for field in full_fields
        fill_halo_regions!(field, args...; kwargs...)
    end
end

function fill_halo_regions!(field::MultiRegionField, args...; kwargs...)
    reduced_dims = reduced_dimensions(field)

    return fill_halo_regions!(field.data,
                              field.boundary_conditions,
                              instantiated_location(field),
                              field.grid,
                              field.boundary_buffers,
                              args...;
                              reduced_dimensions = reduced_dims,
                              kwargs...)
end

fill_halo_regions!(c::MultiRegionObject, ::Nothing, args...; kwargs...) = nothing

#####
##### fill_halo_regions! for a MultiRegionObject
#####

# this can be used once we don't need to fill Value, Flux and Gradient anymore!
# fill_halo_regions!(c::MultiRegionObject, bcs, loc, mrg::MultiRegionGrid, buffers, args...; kwargs...) =
#     apply_regionally!(fill_halo_regions!, c, bcs, loc, mrg, Reference(c.regions), Reference(buffers.regions), args...; kwargs...)

function fill_halo_regions!(c::MultiRegionObject, bcs, loc, mrg::MultiRegionGrid, buffers, args...; kwargs...)

    arch = architecture(mrg)

    halo_tuple = construct_regionally(permute_boundary_conditions, bcs)

    for task = 1:3
        barrier = device_event(arch)
        apply_regionally!(fill_halo_event!, task, halo_tuple,
                          c, loc, arch, barrier, mrg, Reference(c.regions), Reference(buffers.regions),
                          args...; kwargs...)
    end

    return nothing
end

#####
##### fill_halo! for Communicating boundary condition
#####

## Fill communicating boundary condition halos
for (lside, rside) in zip([:west, :south, :bottom], [:east, :north, :bottom])
    fill_both_halo! = Symbol(:fill_, lside, :_and_, rside, :_halo!)
    fill_left_halo!  = Symbol(:fill_, lside, :_halo!)
    fill_right_halo! = Symbol(:fill_, rside, :_halo!)

    @eval begin
        function $fill_both_halo!(c, left_bc::CBC, right_bc, loc, arch, dep, grid, args...; kwargs...)
            event = $fill_right_halo!(c, right_bc, arch, dep, grid, args...; kwargs...)
            $fill_left_halo!(c,  left_bc, arch, event, grid, args...; kwargs...)
            return NoneEvent()
        end
        function $fill_both_halo!(c, left_bc, right_bc::CBC, loc, arch, dep, grid, args...; kwargs...)
            event = $fill_left_halo!(c, left_bc, arch, dep, grid, args...; kwargs...)
            $fill_right_halo!(c, right_bc, arch, event, grid, args...; kwargs...)
            return NoneEvent()
        end
    end
end

function fill_west_and_east_halo!(c, westbc::CBC, eastbc::CBC, loc, arch, dep, grid, neighbors, buffers, args...; kwargs...)

    H = halo_size(grid)[1]
    N = size(grid)[1]
    w = neighbors[westbc.condition.from_rank]
    e = neighbors[eastbc.condition.from_rank]

    westdst = buffers[westbc.condition.rank].west.recv
    eastdst = buffers[eastbc.condition.rank].east.recv

    wait(device(arch), dep)

    switch_device!(getdevice(w))
    westsrc = buffers[westbc.condition.from_rank].east.send
    westsrc .= view(parent(w), N+1:N+H, :, :)

    switch_device!(getdevice(e))
    eastsrc = buffers[eastbc.condition.from_rank].west.send
    eastsrc .= view(parent(e), H+1:2H, :, :)

    switch_device!(getdevice(c))
    device_copy_to!(westdst, westsrc)
    device_copy_to!(eastdst, eastsrc)

    view(parent(c), 1:H, :, :)        .= westdst
    view(parent(c), N+H+1:N+2H, :, :) .= eastdst

    return NoneEvent()
end

function fill_south_and_north_halo!(c, southbc::CBC, northbc::CBC, loc, arch, dep, grid, neighbors, buffers, args...; kwargs...)

    H = halo_size(grid)[2]
    N = size(grid)[2]
    s = neighbors[southbc.condition.from_rank]
    n = neighbors[northbc.condition.from_rank]

    southdst = buffers[southbc.condition.rank].south.recv
    northdst = buffers[northbc.condition.rank].north.recv

    wait(device(arch), dep)

    switch_device!(getdevice(s))
    southsrc = buffers[westbc.condition.from_rank].south.send
    southsrc .= view(parent(s), :, N+1:N+H, :)

    switch_device!(getdevice(n))
    northsrc = buffers[eastbc.condition.from_rank].north.send
    northsrc .= view(parent(n), :, H+1:2H, :)

    switch_device!(getdevice(c))
    device_copy_to!(southdst, southsrc)
    device_copy_to!(northdst, northsrc)

    view(parent(c), :, 1:H, :, :)        .= southdst
    view(parent(c), :, N+H+1:N+2H, :, :) .= northdst

    return NoneEvent()
end

#####
##### Single fill_halo! for Communicating boundary condition
#####

function fill_west_halo!(c, bc::CBC, arch, dep, grid, neighbors, buffers, args...; kwargs...)

    H = halo_size(grid)[1]
    N = size(grid)[1]
    w = neighbors[bc.condition.from_rank]
    dst = buffers[bc.condition.rank].west.recv

    wait(device(arch), dep)

    switch_device!(getdevice(w))
    src = buffers[bc.condition.from_rank].east.send
    src .= view(parent(w), N+1:N+H, :, :)
    sync_device!(getdevice(w))

    switch_device!(getdevice(c))
    device_copy_to!(dst, src)

    p  = view(parent(c), 1:H, :, :)
    p .= dst

    return nothing
end

function fill_east_halo!(c, bc::CBC, arch, dep, grid, neighbors, buffers, args...; kwargs...)

    H = halo_size(grid)[1]
    N = size(grid)[1]
    e = neighbors[bc.condition.from_rank]
    dst = buffers[bc.condition.rank].east.recv

    wait(device(arch), dep)

    switch_device!(getdevice(e))
    src = buffers[bc.condition.from_rank].west.send
    src .= view(parent(e), H+1:2H, :, :)
    sync_device!(getdevice(e))

    switch_device!(getdevice(c))
    device_copy_to!(dst, src)

    p  = view(parent(c), N+H+1:N+2H, :, :)
    p .= dst

    return nothing
end

function fill_south_halo!(c, bc::CBC, arch, dep, grid, neighbors, buffers, args...; kwargs...)

    H = halo_size(grid)[2]
    N = size(grid)[2]
    s = neighbors[bc.condition.from_rank]
    dst = buffers[bc.condition.rank].south.recv

    wait(device(arch), dep)

    switch_device!(getdevice(s))
    src = buffers[bc.condition.from_rank].north.send
    src .= view(parent(s), :, N+1:N+H, :)
    sync_device!(getdevice(s))

    switch_device!(getdevice(c))
    device_copy_to!(dst, src)

    p  = view(parent(c), :, 1:H, :)
    p .= dst

    return nothing
end

function fill_north_halo!(c, bc::CBC, arch, dep, grid, neighbors, buffers, args...; kwargs...)

    H = halo_size(grid)[2]
    N = size(grid)[2]
    n = neighbors[bc.condition.from_rank]
    dst = buffers[bc.condition.rank].north.recv

    wait(device(arch), dep)

    switch_device!(getdevice(n))
    src = buffers[bc.condition.from_rank].south.send
    src .= view(parent(n), :, H+1:2H, :)
    sync_device!(getdevice(n))

    switch_device!(getdevice(c))
    device_copy_to!(dst, src)

    p  = view(parent(c), :, N+H+1:N+2H, :)
    p .= dst

    return nothing
end

#####
##### MultiRegion boundary condition utils
#####

@inline getregion(fc::FieldBoundaryConditions, i) =
        FieldBoundaryConditions(_getregion(fc.west, i),
                                _getregion(fc.east, i),
                                _getregion(fc.south, i),
                                _getregion(fc.north, i),
                                _getregion(fc.bottom, i),
                                _getregion(fc.top, i),
                                fc.immersed)

@inline getregion(bc::BoundaryCondition, i) = BoundaryCondition(bc.classification, _getregion(bc.condition, i))

@inline getregion(cf::ContinuousBoundaryFunction{X, Y, Z, I}, i) where {X, Y, Z, I} =
    ContinuousBoundaryFunction{X, Y, Z, I}(cf.func::F,
                                           _getregion(cf.parameters, i),
                                           cf.field_dependencies,
                                           cf.field_dependencies_indices,
                                           cf.field_dependencies_interp)

@inline getregion(df::DiscreteBoundaryFunction, i) =
    DiscreteBoundaryFunction(df.func, _getregion(df.parameters, i))

@inline _getregion(fc::FieldBoundaryConditions, i) =
FieldBoundaryConditions(getregion(fc.west, i),
                        getregion(fc.east, i),
                        getregion(fc.south, i),
                        getregion(fc.north, i),
                        getregion(fc.bottom, i),
                        getregion(fc.top, i),
                        fc.immersed)

@inline _getregion(bc::BoundaryCondition, i) = BoundaryCondition(bc.classification, getregion(bc.condition, i))

@inline _getregion(cf::ContinuousBoundaryFunction{X, Y, Z, I}, i) where {X, Y, Z, I} =
ContinuousBoundaryFunction{X, Y, Z, I}(cf.func::F,
                                   getregion(cf.parameters, i),
                                   cf.field_dependencies,
                                   cf.field_dependencies_indices,
                                   cf.field_dependencies_interp)

@inline _getregion(df::DiscreteBoundaryFunction, i) =
DiscreteBoundaryFunction(df.func, getregion(df.parameters, i))

# Everything goes for multi-region BC
validate_boundary_condition_location(::MultiRegionObject, ::Center, side)       = nothing
validate_boundary_condition_location(::MultiRegionObject, ::Face, side)         = nothing
validate_boundary_condition_topology(::MultiRegionObject, topo::Periodic, side) = nothing
validate_boundary_condition_topology(::MultiRegionObject, topo::Flat,     side) = nothing


text:
#Pkg.add("Example")

using Example

hello("it's me. I was wondering if after all these years you'd like to meet.")

#Pkg.add("Colors")

using Colors

bandeja = distinguishable_colors(100)

rand(bandeja, 3, 3)

max_stars_repo_name: miguelraz/JuliaTutorials

max_issues_repo_name: miguelraz/JuliaTutorials

max_forks_repo_name: miguelraz/JuliaTutorials


text:
using Test, CUDA, Random, Knet, AutoGrad
using CUDA.CUDNN:
    cudnnRNNForward,
    cudnnRNNForward!,
    cudnnRNNBackwardData_v8,
    cudnnRNNBackwardWeights_v8,
    cudnnRNNDescriptor,
    cudnnRNNDescriptor_t,
    cudnnSetRNNDescriptor_v8,
    cudnnGetRNNWeightSpaceSize,
    cudnnGetRNNTempSpaceSizes,
    cudnnRNNAlgo_t,
        CUDNN_RNN_ALGO_STANDARD,        # 0, robust performance across a wide range of network parameters
        CUDNN_RNN_ALGO_PERSIST_STATIC,  # 1, fast when the first dimension of the input tensor is small (meaning, a small minibatch), cc>=6.0
        CUDNN_RNN_ALGO_PERSIST_DYNAMIC, # 2, similar to static, optimize using the specific parameters of the network and active GPU, cc>=6.0
        CUDNN_RNN_ALGO_COUNT,           # 3
    cudnnRNNMode_t,
        CUDNN_RNN_RELU, # 0, /* basic RNN cell type with ReLu activation */
        CUDNN_RNN_TANH, # 1, /* basic RNN cell type with tanh activation */
        CUDNN_LSTM,     # 2, /* LSTM with optional recurrent projection and clipping */
        CUDNN_GRU,      # 3, /* Using h' = tanh(r * Uh(t-1) + Wx) and h = (1 - z) * h' + z * h(t-1); */
    cudnnRNNBiasMode_t,
        CUDNN_RNN_NO_BIAS,         # 0, /* rnn cell formulas do not use biases */
        CUDNN_RNN_SINGLE_INP_BIAS, # 1, /* rnn cell formulas use one input bias in input GEMM */
        CUDNN_RNN_DOUBLE_BIAS,     # 2, /* default, rnn cell formulas use two bias vectors */
        CUDNN_RNN_SINGLE_REC_BIAS, # 3  /* rnn cell formulas use one recurrent bias in recurrent GEMM */
    cudnnDirectionMode_t,
        CUDNN_UNIDIRECTIONAL, # 0, /* single direction network */
        CUDNN_BIDIRECTIONAL,  # 1, /* output concatination at each layer */
    cudnnRNNInputMode_t,
        CUDNN_LINEAR_INPUT, # 0, /* adjustable weight matrix in first layer input GEMM */
        CUDNN_SKIP_INPUT,   # 1, /* fixed identity matrix in the first layer input GEMM */
    cudnnMathType_t,
        CUDNN_DEFAULT_MATH,                    # 0,
        CUDNN_TENSOR_OP_MATH,                  # 1,
        CUDNN_TENSOR_OP_MATH_ALLOW_CONVERSION, # 2,
        CUDNN_FMA_MATH,                        # 3,
    #/* For auxFlags in cudnnSetRNNDescriptor_v8() and cudnnSetRNNPaddingMode() */
        CUDNN_RNN_PADDED_IO_DISABLED, # 0
        CUDNN_RNN_PADDED_IO_ENABLED,  # (1U << 0)
    cudnnForwardMode_t,
        CUDNN_FWD_MODE_INFERENCE, # 0
        CUDNN_FWD_MODE_TRAINING,  # 1
    cudnnRNNDataDescriptor_t,
    cudnnSetRNNDataDescriptor,
    cudnnRNNDataLayout_t,
        CUDNN_RNN_DATA_LAYOUT_SEQ_MAJOR_UNPACKED,   # 0, /* padded, outer stride from one time-step to the next */
        CUDNN_RNN_DATA_LAYOUT_SEQ_MAJOR_PACKED,     # 1, /* sequence length sorted and packed as in basic RNN api */
        CUDNN_RNN_DATA_LAYOUT_BATCH_MAJOR_UNPACKED, # 2, /* padded, outer stride from one batch to the next */
    cudnnWgradMode_t,
        CUDNN_WGRAD_MODE_ADD, # 0, /* add partial gradients to wgrad output buffers */
        CUDNN_WGRAD_MODE_SET, # 1, /* write partial gradients to wgrad output buffers */
    cudnnTensorDescriptor,
    cudnnDropoutDescriptor,
    cudnnDataType,
    math_mode,
    handle

@testset "cudnn/rnn" begin
    X,H,B,T,L,F = 4,4,4,4,1,Float64
    aw,ax,ah,ac = randn(F,10000), randn(F,X,B,T), randn(F,H,B,L), randn(F,H,B,L)
    cw,cx,ch,cc = Param.(CuArray.((aw,ax,ah,ac)))

    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H)
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,hx=ch)
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,cx=cc)
    @test @gcheck (hy=Ref{Any}(); y=cudnnRNNForward(cw,cx;hiddenSize=H,hy); sum(y)+sum(hy[]))
    @test @gcheck (cy=Ref{Any}(); y=cudnnRNNForward(cw,cx;hiddenSize=H,cy); sum(y)+sum(cy[]))
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,layout=CUDNN_RNN_DATA_LAYOUT_SEQ_MAJOR_PACKED)
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,layout=CUDNN_RNN_DATA_LAYOUT_BATCH_MAJOR_UNPACKED)
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,seqLengthArray=Cint[1,2,1,2])
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,fwdMode=CUDNN_FWD_MODE_TRAINING)
    #@test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,algo=CUDNN_RNN_ALGO_PERSIST_STATIC)  # not supported
    #@test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,algo=CUDNN_RNN_ALGO_PERSIST_DYNAMIC) # causes segfault
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,cellMode=CUDNN_RNN_RELU)
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,cellMode=CUDNN_RNN_TANH)
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,cellMode=CUDNN_GRU)
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,biasMode=CUDNN_RNN_NO_BIAS)
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,biasMode=CUDNN_RNN_SINGLE_INP_BIAS)
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,biasMode=CUDNN_RNN_SINGLE_REC_BIAS)
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,dirMode=CUDNN_BIDIRECTIONAL)
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,inputMode=CUDNN_SKIP_INPUT)
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,mathPrec=Float64)
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,mathType=CUDNN_DEFAULT_MATH)
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,mathType=CUDNN_TENSOR_OP_MATH)
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,mathType=CUDNN_TENSOR_OP_MATH_ALLOW_CONVERSION)
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,projSize=H-1)
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,numLayers=2)
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,dropout=0.5)
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,auxFlags=CUDNN_RNN_PADDED_IO_DISABLED)
    @test @gcheck cudnnRNNForward(cw,cx;hiddenSize=H,auxFlags=CUDNN_RNN_PADDED_IO_ENABLED)
end


text:
#################################################################
# Randomized estimators of elementary linear algebraic quantities
#################################################################

export rcond, reigmax, reigmin, rnorm, rnorms

@doc doc"""
Randomized Gaussian matrices normalized by column

Input:
    `el`: element type
    `m`: number of rows
    `n`: number of columns or nothing
    `normalize`: whether or not to normalize (default: `true`)

Output:
    `Ω`: a matrix of dimensions `m` x `n` containing Gaussian random numbers of
    type `el`.
""" ->
function randnn(el, m::Int, normalize::Bool=true)
    if el <: Real
        Ω = randn(m)
    elseif el <: Complex
        Ω = randn(m) + im*randn(m)
    else
        throw(ValueError("Unsupported element type: $el"))
    end
    normalize ? Ω/norm(Ω) : Ω
end
function randnn(el, m::Int, n::Int, normalize::Bool=true)
    if el <: Real
        Ω = randn(m, n)
    elseif el <: Complex
        Ω = randn(m, n) + im*randn(m, n)
    else
        throw(ValueError("Unsupported element type: $el"))
    end
    normalize || return Ω
    for i=1:n
        Ω[:, i] /= norm(sub(Ω, :, i))
    end
    Ω
end

@doc doc"""
Randomized matrix norm estimator

Computes a probabilistic upper bound on the norm of a matrix `A`.

\cite[Lemma 4.1]{Halko2011} states (with slight notational change) that

‖A‖ ≤ α √(2/π) maxᵢ ‖Aωᵢ‖

with probability $$p=α^{-r}$$.

Inputs:
    `A`: Matrix whose norm to estimate
    `r`: Number of matrix-vector products to compute
    `p`: Probability of upper bound failing (default: 0.05)

Output:
    Estimate of ‖A‖.

See also:
    `rnorms()` for a different estimator that uses
    premultiplying by both `A` and `A'`
""" ->
function rnorm(A, r::Int, p::Real=0.05)
    @assert 0<p≤1
    α = p^(-1.0/r)
    m, n = size(A)
    Ω = randnn(eltype(A), n, r, false)
    AΩ = A*Ω
    mx = maximum([norm(sub(AΩ, :, j)) for j=1:r])
    α * √(2/π) * mx
end

@doc doc"""
Randomized matrix norm estimator using `A'A`

Computes a probabilistic upper bound on the norm of a matrix `A`.

\cite[Appendix]{Liberty2007} states (with minor change in notation) that

ρ = √(‖(A'A)ʲω‖/‖(A'A)ʲ⁻¹ω‖)

which is an estimate of the spectral norm of A produced by j
steps of the power method starting with normalized ω, is a lower
bound on the true norm by a factor

ρ ≤ α ‖A‖

with probability greater than $$1 - p$$, where
$$p = 4\sqrt{n/(j-1)} α^{-2j}$$.

Inputs:
    `A`: Matrix whose norm to estimate
    `j`: Number of power iterations to perform. (Default: 1)
    `p`: Probability of upper bound failing. (Default: 0.05)
    `At` (optional keyword): Transpose of `A`. (Default: `A'`)

Output:
    Estimate of ‖A‖.

Reference:
    Appendix of \cite{Liberty2007}.

    @article{Liberty2007,
        authors = {Edo Liberty and Franco Woolfe and Per-Gunnar Martinsson
	    and Vladimir Rokhlin and Mark Tygert},
        title = {Randomized algorithms for the low-rank approximation of matrices},
        journal = {Proceedings of the National Academy of Sciences},
	volume = {104},
	issue = {51},
	year = {2007},
	pages = {20167--20172},
        doi  = {10.1073/pnas.0709640104}
    }

Comment:
    see `rnorm()` for a different estimator that does not require
    premultiplying by `A'`
""" ->
function rnorms(A, j::Int=1, p::Real=0.05; At = A')
    @assert 0<p≤1
     m, n = size(A)
    α = ((j-1)/n*(p/4)^2)^(-1/(4j))

    Ωold = Ω = randnn(eltype(A), n)
    for i=1:j #Power iterations
        Ω, Ωold = At*(A*Ω), Ω
    end
    ρ = √(norm(Ω)/norm(Ωold))
    α*ρ
end

@doc doc"""
Randomized condition number estimator

Inputs:
    `A`: Matrix whose condition number to estimate.
         Must be square and support premultiply (`A*⋅`) and solve (`A\⋅`)
    `k`: Number of power iterations to run. (Default: 1, recommended: `k ≤ 3`)
    `p`: Probability that estimate fails to hold as an upper bound
       (Default: 0.05)

Output:

    The interval `(x, y)` which contains `κ(A)` with probability $$1 - p$$.

Implementation note:

    \cite{Dixon1983} originally describes this as a computation that
    can be done by computing the necessary number of power iterations given p
    and the desired accuracy parameter θ=y/x. However, these bounds were only
    derived under the assumptions of exact arithmetic. Empirically, k≥4 has
    been seen to result in incorrect results in that the computed interval does
    not contain the true condition number. This implemention therefore makes `k`
    an explicitly user-controllable parameter from which to infer the accuracy
    parameter and hence the interval containing κ(A).

Reference:
    \cite[Theorem 2]{Dixon1983}

    @article{Dixon1983,
        author = {Dixon, John D},
        doi = {10.1137/0720053},
        journal = {SIAM Journal on Numerical Analysis},
        number = {4},
        pages = {812--814},
	title = {Estimating Extremal Eigenvalues and Condition Numbers of
		Matrices},
	volume = {20},
        year = {1983}
    }
""" ->
function rcond(A, k::Int=1, p::Real=0.05)
    @assert 0<p≤1
    m, n = size(A)
    @assert m==n
    θ = (8n/(π*p^2))^(1/k)
    x = randnn(eltype(A), n)
    for i=1:k
        x = A*x
    end
    y = randnn(eltype(A), n)
    for i=1:k
        y = A\y
    end
    φ = ((x⋅x)*(y⋅y))^(1/(2k))
    (φ, θ*φ)
end

@doc doc"""
Randomized maximal eigenvalue estimator

Inputs:

    `A`: Matrix whose maximal eigenvalue to estimate.
         Must be square and support premultiply (`A*⋅`)
    `k`: Number of power iterations to run. (Default: 1, recommended: k ≤ 3)
    `p`: Probability that estimate fails to hold as an upper bound
       (Default: 0.05)

Output:

    The interval `(x, y)` which contains the maximal eigenvalue of `A` with
    probability $$1 - p$$.

Reference:

    \cite[Corollary of Theorem 1]{Dixon1983}.
""" ->
function reigmax(A, k::Int=1, p::Real=0.05)
    @assert 0<p≤1
    m, n = size(A)
    @assert m==n
    θ = (2n/(π*p^2))^(1/k)
    y = x = randnn(eltype(A), n)
    for i=1:k
        x = A*x
    end
    φ = y⋅x
    (φ, θ*φ)
end

@doc doc"""
Randomized minimal eigenvalue estimator

Inputs:

    `A`: Matrix whose minimal eigenvalue to estimate.
         Must be square and support backslash (`A\⋅`)
    `k`: Number of power iterations to run. (Default: 1, recommended: k ≤ 3)
    `p`: Probability that estimate fails to hold as an upper bound
         (Default: 0.05)

Output:

    The interval `(x, y)` which contains the minimal eigenvalue of `A` with
    probability $$1 - p$$.

Reference:

    \cite[Corollary of Theorem 1]{Dixon1983}.
""" ->
function reigmin(A, k::Int=1, p::Real=0.05)
    @assert 0<p≤1
    m, n = size(A)
    @assert m==n
    θ = (2n/(π*p^2))^(1/k)
    y = x = randnn(eltype(A), n)
    for i=1:k
        x = A\x
    end
    φ = y⋅x
    (inv(θ*φ), inv(φ))
end

@doc doc"""
A subsampled random Fourier transform

Parameter:

    l :: Number of vectors to return
""" ->
immutable srft{T<:Integer}
    l :: T
end

(*)(A::AbstractMatrixFcn, B::srft) = error("method only defined to avoid ambiguity. If you need this method please open a pull request")

@doc doc"""
Applies a subsampled random Fourier transform to the columns of `A`

Inputs:

    `A`: A matrix to transform
    `Ω`: A `srft` type

Output:
    `B`: A matrix of dimensions size(A,1) x Ω.l

Reference:

    \[Equation 4.6]{Halko2011}
""" ->
function *(A, Ω::srft)
    m, n = size(A)
    B = A*Diagonal(exp(2π*im*rand(n))/√Ω.l)
    B = vcat([fft(A[i,:]) for i=1:m]...) #Factor of √n cancels out
    B[:, randperm(n)[1:Ω.l]]
end


text:
struct LDResampler
    lv::LowVarianceResampler
end

LDResampler(n::Int) = LDResampler(LowVarianceResampler(n))

ParticleFilters.n_init_samples(r::LDResampler) = n_init_samples(r.lv)

function ParticleFilters.resample(r::LDResampler,
                                  bp::WeightedParticleBelief,
                                  pm::LightDark1D,
                                  rm::LightDark1D,
                                  b,
                                  a,
                                  o,
                                  rng)
    if a == 0
        return ParticleCollection([LightDark1DState(-1, 0.0)])
    end
    return resample(r.lv, bp, rng)
end

@testset "domain_specific" begin

n = 100
m = LightDark1D()
up = BasicParticleFilter(m, LDResampler(n), n)
p = FunctionPolicy(b->0)

bp = first(stepthrough(m, p, up, "bp"))
@test first(particles(bp)) == LightDark1DState(-1, 0.0)

p2 = FunctionPolicy(b->2)
for bp in stepthrough(m, p2, up, "bp", max_steps=3)
    @test bp isa ParticleCollection{LightDark1DState}
    @test n_particles(bp) == n
end

end


text:
#!/usr/bin/env julia

using Distributed
using Random
using LinearAlgebra
procs_id = addprocs(4)
using DatagenCopulaBased
@everywhere using Distributions
@everywhere using Cumulants
@everywhere using SymmetricTensors
using CumulantsFeatures
using JLD2
using FileIO
using ArgParse
@everywhere import CumulantsFeatures: reduceband
@everywhere using DatagenCopulaBased
@everywhere using CumulantsFeatures
@everywhere using HypothesisTests

Random.seed!(1234)

@everywhere function gmarg2t(X::Matrix{T}, nu::Int) where T <: AbstractFloat
  # transform univariate gasussian marginals into t-Student with nu degrees of freedom
  Y = copy(X)
  for i = 1:size(X, 2)
    x = X[:,i]
    s = var(x)
    mu = mean(x)
    d = Normal(mu, s)
    u = cdf.(d, x)
    pvalue(ExactOneSampleKSTest(u,Uniform(0,1)))>0.0001 || throw(AssertionError("$i marg. not unif."))
    Y[:,i] = quantile.(TDist(nu), u)
  end
  return Y
end

function main(args)
  s = ArgParseSettings("description")
  @add_arg_table! s begin
    "--nu", "-n"
    default = 6
    help = "the number of degrees of freedom for the t-Student copula"
    arg_type = Int

    "--nuu", "-u"
    default = 6
    help = "the number of degrees of freedom for the t-Student marginal"
    arg_type = Int

    "--reals", "-N"
    default = 5
    help = "number of tests realisations"
    arg_type = Int
  end
  parsed_args = parse_args(s)
  ν = parsed_args["nu"]
  νu = parsed_args["nuu"]
  # number of generated data sets
  test_number = parsed_args["reals"]

  println("copula's degreen of freedom = ", ν)
  println("matginal's degree of freedom = ", νu)
  # parameters, data sie
  @everywhere t = 1000
  # no of marginals
  @everywhere n = 30
  # no of marginals with t-Student copula
  @everywhere malf_size = 15
  # outliers places on the beginning of data
  @everywhere a = 100
  data_dir = "."

  filename = "tstudent_$(ν)_marg$(νu)-t_size-$(n)_malfsize-$malf_size-t_$(t)_$a.jld2"

  data = Dict{String, Any}("variables_no" => n,
                         "sample_number" => t,
                         "ν" => ν,
                         "test_number" => test_number,
                         "malf_size" => malf_size,
                         "a" => a,
                         "data" => Dict{String, Dict{String,Any}}())

  known_data_size = 0
  if isfile("./$filename")
   data["data"] = load("./$filename")["data"]
   known_data_size += length(data["data"])
   println("Already have $known_data_size samples \n Will generate $(test_number-known_data_size) more")
  end

  println("Calculation started")
  for m=(known_data_size+1):test_number
    @time begin
      println(" > $m ($ν)")
      malf = randperm(n)[1:malf_size]
      Σ = cormatgen_rand(n)
      samples_orig = rand(MvNormal(Σ), t)'

      # gcop2tstudent(x[1:a, :], malf, ν) - outliers, given gaussian multivariate copula is changed for Gaussian
      # x[a+1:end, :] - ordinary data
      # gmarg2t( ..., νu) - univariate marginals are changed to t-Student
      # "malf" - dict key

      versions = [(x->gmarg2t(vcat(gcop2tstudent(x[1:a, :], malf, ν), x[a+1:end, :]), νu), "malf")]

      cur_dict = Dict{String, Any}("malf" => malf,
                                   "cor_source" => Σ)

      data_dict = @distributed (merge) for (sampler, label)=versions
        println(label)
        samples = sampler(samples_orig)
        Σ_malf = cov(samples)
        Dict("cor_$label" => Σ_malf,
             "x_$label" => samples)
      end

      data["data"]["$m"] = merge(cur_dict, data_dict)
      save("./$filename", data)
    end
  end
end

main(ARGS)