text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -140.37215752831116 <= q <= 140.41841564400983)

                     add_NL_constraint(m, :(1/(1 + exp(-(0.9038900968817085 + 0.1284043831226005*1/(1 + exp(-(0.06165362792774953 + 0.6866075626761576*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + -0.8737233458006566*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + -0.3942911833193632*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + 0.6872447034250162*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + 0.5711071882051262*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + 0.26042276084039084*1/(1 + exp(-(0.8058455478649482 + 0.7751639240502244*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + -0.13474005943220124*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + 0.6437669970537465*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + 0.5932192098608322*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + 0.5325023046383324*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + 0.6458808355119907*1/(1 + exp(-(-0.1540966309548315 + 0.690406062481447*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + 0.4718956794952942*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + 0.8703243890856593*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + -0.8079540211479861*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + 0.683716372184894*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + -0.6567358944011805*1/(1 + exp(-(-0.2486694730930381 + -0.9006272851005916*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + -0.9841514005441985*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + 0.7501520395123897*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + -0.43486617405885397*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + -0.4797454164387962*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + -0.41502277707459934*1/(1 + exp(-(-0.9928729682863509 + 0.9588446333785003*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + 0.7154937103302039*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + -0.17921833371514895*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + -0.7009330672719543*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + 0.15311755902026114*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5])))))))))) + 1/(1 + exp(-(-0.8820041740138964 + 0.8813585352462416*1/(1 + exp(-(0.06165362792774953 + 0.6866075626761576*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + -0.8737233458006566*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + -0.3942911833193632*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + 0.6872447034250162*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + 0.5711071882051262*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + -0.30334290832683886*1/(1 + exp(-(0.8058455478649482 + 0.7751639240502244*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + -0.13474005943220124*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + 0.6437669970537465*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + 0.5932192098608322*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + 0.5325023046383324*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + 0.6303754843616343*1/(1 + exp(-(-0.1540966309548315 + 0.690406062481447*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + 0.4718956794952942*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + 0.8703243890856593*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + -0.8079540211479861*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + 0.683716372184894*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + 0.9344814932534002*1/(1 + exp(-(-0.2486694730930381 + -0.9006272851005916*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + -0.9841514005441985*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + 0.7501520395123897*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + -0.43486617405885397*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + -0.4797454164387962*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + -0.8925753470795712*1/(1 + exp(-(-0.9928729682863509 + 0.9588446333785003*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + 0.7154937103302039*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + -0.17921833371514895*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + -0.7009330672719543*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + 0.15311755902026114*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5])))))))))) + 1/(1 + exp(-(-0.6413973695460733 + -0.5291254814289021*1/(1 + exp(-(0.06165362792774953 + 0.6866075626761576*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + -0.8737233458006566*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + -0.3942911833193632*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + 0.6872447034250162*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + 0.5711071882051262*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + -0.3315514300310589*1/(1 + exp(-(0.8058455478649482 + 0.7751639240502244*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + -0.13474005943220124*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + 0.6437669970537465*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + 0.5932192098608322*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + 0.5325023046383324*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + 0.9076257425527015*1/(1 + exp(-(-0.1540966309548315 + 0.690406062481447*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + 0.4718956794952942*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + 0.8703243890856593*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + -0.8079540211479861*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + 0.683716372184894*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + -0.9763381142518175*1/(1 + exp(-(-0.2486694730930381 + -0.9006272851005916*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + -0.9841514005441985*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + 0.7501520395123897*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + -0.43486617405885397*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + -0.4797454164387962*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + -0.06362970408750801*1/(1 + exp(-(-0.9928729682863509 + 0.9588446333785003*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + 0.7154937103302039*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + -0.17921833371514895*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + -0.7009330672719543*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + 0.15311755902026114*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5])))))))))) + 1/(1 + exp(-(0.9122346082558659 + 0.2864216745322903*1/(1 + exp(-(0.06165362792774953 + 0.6866075626761576*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + -0.8737233458006566*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + -0.3942911833193632*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + 0.6872447034250162*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + 0.5711071882051262*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + 0.7014729451110155*1/(1 + exp(-(0.8058455478649482 + 0.7751639240502244*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + -0.13474005943220124*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + 0.6437669970537465*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + 0.5932192098608322*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + 0.5325023046383324*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + -0.27913045286088156*1/(1 + exp(-(-0.1540966309548315 + 0.690406062481447*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + 0.4718956794952942*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + 0.8703243890856593*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + -0.8079540211479861*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + 0.683716372184894*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + 0.014698309183167169*1/(1 + exp(-(-0.2486694730930381 + -0.9006272851005916*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + -0.9841514005441985*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + 0.7501520395123897*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + -0.43486617405885397*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + -0.4797454164387962*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + -0.9872323021033105*1/(1 + exp(-(-0.9928729682863509 + 0.9588446333785003*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + 0.7154937103302039*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + -0.17921833371514895*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + -0.7009330672719543*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + 0.15311755902026114*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5])))))))))) + 1/(1 + exp(-(0.6835137249354726 + 0.06184589975771315*1/(1 + exp(-(0.06165362792774953 + 0.6866075626761576*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + -0.8737233458006566*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + -0.3942911833193632*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + 0.6872447034250162*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + 0.5711071882051262*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + 0.6514470819403977*1/(1 + exp(-(0.8058455478649482 + 0.7751639240502244*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + -0.13474005943220124*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + 0.6437669970537465*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + 0.5932192098608322*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + 0.5325023046383324*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + 0.1907919388573478*1/(1 + exp(-(-0.1540966309548315 + 0.690406062481447*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + 0.4718956794952942*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + 0.8703243890856593*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + -0.8079540211479861*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + 0.683716372184894*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + 0.3162639033084287*1/(1 + exp(-(-0.2486694730930381 + -0.9006272851005916*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + -0.9841514005441985*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + 0.7501520395123897*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + -0.43486617405885397*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + -0.4797454164387962*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5]))))))) + 0.2338844135686573*1/(1 + exp(-(-0.9928729682863509 + 0.9588446333785003*1/(1 + exp(-(0.28804974063451105 + 0.7911631812636903*$(x[1]) + 0.6415196610101699*$(x[2]) + -0.09319048786154349*$(x[3]) + 0.5261383650120517*$(x[4]) + -0.5505775194565721*$(x[5])))) + 0.7154937103302039*1/(1 + exp(-(0.7342872210188394 + 0.669146525148411*$(x[1]) + 0.2634123220650846*$(x[2]) + -0.9804313836692913*$(x[3]) + 0.770312239258653*$(x[4]) + 0.24263357188873957*$(x[5])))) + -0.17921833371514895*1/(1 + exp(-(-0.6886620050375161 + 0.8091467022606786*$(x[1]) + 0.910808980108254*$(x[2]) + -0.3183669784784451*$(x[3]) + 0.5442365253482557*$(x[4]) + -0.6054042374016482*$(x[5])))) + -0.7009330672719543*1/(1 + exp(-(0.5813944367886887 + 0.520001186683257*$(x[1]) + -0.3951433421155799*$(x[2]) + 0.803037085072233*$(x[3]) + 0.9067157149850686*$(x[4]) + 0.001120584088625698*$(x[5])))) + 0.15311755902026114*1/(1 + exp(-(-0.16131179027810427 + 0.7148051479719681*$(x[1]) + 0.7364049274130124*$(x[2]) + -0.647777979207961*$(x[3]) + 0.24628792618437867*$(x[4]) + 0.6449790156418587*$(x[5])))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module ITensorGaussianMPS

using ITensors
using LinearAlgebra

import LinearAlgebra: Givens

export slater_determinant_to_mps,
       slater_determinant_to_gmps,
       hopping_hamiltonian,
       slater_determinant_matrix

include("gmps.jl")

end

ext: jl


text:
# Instead of using cudaStreamSynchronize, this implements
# a version that waits on a condition. This allows you to
# use this with Julia Task interface. A practical application
# is that one Julia process can feed multiple devices, with the jobs
# running asynchronously.
import Base: AsyncCondition
@compat abstract type AbstractStream end

type Stream <: AbstractStream
    inner::CUDAdrv.CuStream
    c::AsyncCondition
end
function Stream()
    p = Ref{Ptr{Void}}()
    rt.cudaStreamCreate(p)
    ctx = contexts[device()]
    hnd = CUDAdrv.CuStream(p[], ctx)
    Stream(hnd, AsyncCondition())
end
NullStream() = Stream(CUDAdrv.CuDefaultStream(), AsyncCondition())
const null_stream = NullStream()

convert(::Type{CUDAdrv.CuStream}, s::Stream) = s.inner

unsafe_convert(::Type{Ptr{Void}}, s::Stream) = unsafe_convert(Ptr{Void}, s.inner)

synchronize(s::Stream) = synchronize(s.inner)

function wait(s::Stream)
    rt.cudaStreamAddCallback(s, c_async_send_cudastream[], Base.unsafe_convert(Ptr{Void}, s.c), 0)
    wait(s.c)
    rt.checkerror(rt.cudaGetLastError())
end

async_send_cudastream(hnd, status, data) = (ccall(:uv_async_send,Cint,(Ptr{Void},),data); nothing)

const c_async_send_cudastream = Ref{Ptr{Void}}()


text:
#!/usr/bin/env julia

using SEM

using Plots, LinearAlgebra, SparseArrays

#--------------------------------------#
# setup
#--------------------------------------#
nx1 = 10; Ex = 2;
ny1 = 10; Ey = 2;
#--------------------------------------#
# nodal operators
#--------------------------------------#
zr1,wr1 = SEM.FastGaussQuadrature.gausslobatto(nx1)
zs1,ws1 = SEM.FastGaussQuadrature.gausslobatto(ny1)

Dr1 = derivMat(zr1); Ds1 = derivMat(zs1)

#----------------------------------------------------------------------#
# local operators
#----------------------------------------------------------------------#
wx1 = kron(ones(Ex,1),wr1); wy1 = kron(ones(Ey,1),ws1);

Iex = Matrix(I,Ex,Ex);
Iey = Matrix(I,Ey,Ey);

Dx1 = kron(Iex,Dr1); Dy1 = kron(Iey,Ds1);

#----------------------------------------------------------------------#
# boundary conditions
#----------------------------------------------------------------------#
ifperiodicX = false
ifperiodicY = false

Ix1 = Matrix(I,Ex*nx1,Ex*nx1);
Iy1 = Matrix(I,Ey*ny1,Ey*ny1);

Rx1 = Ix1[2:end-1,:];
Ry1 = Iy1[2:end-1,:];

if(ifperiodicX) Rx1 = Ix1; end
if(ifperiodicY) Ry1 = Iy1; end

M1 = diag(Rx1'*Rx1) * diag(Ry1'*Ry1)';

Ix1 = nothing
Iy1 = nothing
Rx1 = nothing
Ry1 = nothing
#----------------------------------------------------------------------#
# mapping
#----------------------------------------------------------------------#

# Q: global -> local op, Q': local -> global
Qx1 = semq(Ex,nx1,ifperiodicX)
Qy1 = semq(Ey,ny1,ifperiodicY)

# gather scatter op
QQtx1 = Qx1*Qx1';
QQty1 = Qy1*Qy1';

# mult
mult = ones(nx1*Ex,ny1*Ey);
mult = gatherScatter(mult,QQtx1,QQty1);
mult = @. 1 / mult;

mult1 = reshape(mult,:)

#----------------------------------------------------------------------#
# geometry
#----------------------------------------------------------------------#
x1e,_ = semmesh(Ex,nx1); y1e,_ = semmesh(Ey,ny1);

x1,y1 = ndgrid(x1e,y1e);

# deform grid with gordonhall
x1 = @. 0.5 * (x1 + 1); y1 = @. 0.5 * (y1 + 1);

Jac1,Jaci1,rx1,ry1,sx1,sy1 = jac(x1,y1,Dx1,Dy1);

#----------------------------------------------------------------------#
# case setup
#----------------------------------------------------------------------#
# prescribe forcing, true solution, boundary data
kx=1.
ky=1.
ut = @. sin(kx*pi*x1)*sin(ky*pi*y1) # true solution
f  = @. ut*((kx^2+ky^2)*pi^2);      # forcing/RHS
ub = copy(ut);                      # boundary data

visc = @. 1+0*x1;

#f = @. 1+0*x1;
ub= @. 0+0*x1;
#----------------------------------------------------------------------#
# setup
#----------------------------------------------------------------------#

# diagonal mass matrix
B1  = Jac1 .* (wx1*wy1');
Bi1 = 1 ./ B1;

# set up matrices for Laplace op.
G11 = @. visc * B1 * (rx1 * rx1 + ry1 * ry1);
G12 = @. visc * B1 * (rx1 * sx1 + ry1 * sy1);
G22 = @. visc * B1 * (sx1 * sx1 + sy1 * sy1);
#----------------------------------------------------------------------#
# Explicit system
#----------------------------------------------------------------------#
# operators
Ixl = Matrix(I,Ex*nx1,Ex*nx1);              # Identity on local vectors
Iyl = Matrix(I,Ey*ny1,Ey*ny1);
Ixg = Matrix(I,Ex*(nx1-1)+1,Ex*(nx1-1)+1);  # Identity on global vectors
Iyg = Matrix(I,Ey*(ny1-1)+1,Ey*(ny1-1)+1);
Rxl = Ixl[2:end-1,:]; Rxg = Ixg[2:end-1,:]; # restriction on loc/gl vectors
Ryl = Iyl[2:end-1,:]; Ryg = Iyg[2:end-1,:];
Rl  = kron(Ryl,Rxl);  Rg  = kron(Ryg,Rxg);
Q   = kron(Qy1,Qx1);
Dr  = kron(Iyl,Dx1);
Ds  = kron(Dy1,Ixl);
Drs = [Dr;Ds];
B   = Diagonal(reshape(B1 ,:))
g11 = Diagonal(reshape(G11,:))
g12 = Diagonal(reshape(G12,:))
g22 = Diagonal(reshape(G22,:))
G   = [g11 g12; g12 g22]
A   = Drs' * G * Drs
fl  = reshape(f,:) # forcing

# rank deficient system acting on local, unrestricted operators
M    = Rl'*Rl
QQt  = Q*Q'
Aloc = QQt * M * A # equivale to (Q * Rg' * Rg * Q') * A
Bloc = QQt * M * B
bloc = Bloc * fl
@time uloc = pcg(bloc,Aloc;mult=mult1)
uloc = reshape(uloc,Ex*nx1,Ey*ny1)

#----------------------------------------------------------------------#
laplOp(v) = lapl(v,M1,[],[],QQtx1,QQty1,Dx1,Dy1,G11,G12,G22,mult);
#laplFdmOp(v) = lapl_fdm(v,Bi1,Sx,Sy,Sxi,Syi,Lfdmi);

b   = mass(f,[],B1,[],[],[],[],mult)
b .-= lapl(ub,[],[],[],[],[],Dx1,Dy1,G11,G12,G22,mult)
b   = mass(b,M1,[],[],[],QQtx1,QQty1,mult)

@time u = pcg(b,laplOp;mult=mult)
u = u + ub
#----------------------------------------------------------------------#
print("Rank deficient system, er: ",norm(ut-uloc,Inf),"\n")
print("Functions as operators, er: ",norm(ut-u,Inf),"\n")
p = meshplt(x1,y1,uu)
display(p)
#----------------------------------------------------------------------#
return


text:
"""
    solve(problem::Laplace)

Compute the Laplace eigenvalues, eigenfunctions and first order moments of products of pairs of eigenfunctions.
"""
function solve(laplace::Laplace)
    (; model, matrices, neig_max) = laplace
    (; mesh) = model
    (; M, S, R, Mx, Q) = matrices

    # Compute at most all eigenvalues in the given domain
    neig = Int(min(neig_max, size(M, 1)))

    @info join(
        [
            "Solving Laplace eigenvalue problem, computing $neig eigenvalues.",
            "Problem size: $(size(M, 1)) points.",
        ],
        "\n",
    )

    # Solve generalized eigenproblem, computing the smallest eigenvalues only.
    # If 2 * neig_max_domain >= nnode, a full decomposition is performed,
    # calling the eig function inside eigs
    λ, ϕ = eigs(S + Q, M, nev = neig, which = :SR)
    # λ, ϕ = eigs(Hermitian(S + Q), Hermitian(M), nev = neig, which = :SR);
    # λ, ϕ = eigen(Hermitian(Matrix(S + Q)), Hermitian(Matrix(M)))

    # TODO: Consider KrylovKit.jl

    # All Laplace eigenvalues are nonnegative
    all(0 .≤ λ) ||
        @warn "Obtained negative eigenvalues for Laplace operator." findall(λ .< 0) λ[λ.<0]

    # Normalize eigenfunctions with mass weighting
    ϕ ./= .√sum(ϕ .* (M * ϕ), dims = 1)

    # Compute first order moments of product of pairs of eigenfunctions
    moments = [ϕ' * Mx[dim] * ϕ for dim = 1:3]

    # Compute Laplace relaxation matrix
    massrelax = ϕ' * R * ϕ

    values = λ
    funcs = ϕ

    (; values, funcs, moments, massrelax)
end


text:
# CSCI5573Final/config.jl
# Licensed under the MIT License. See LICENSE.md file in the project root for
# full license information.

SIM_RUN = "FCFS"

macro pick_scheduler(NAME)
    quote
        @schedule at 0 $(esc(Symbol(NAME * "Scheduler()")))
    end
end

TRIALS                = 1
N_BUSSES              = 10
MAX_TASKS             = 250
CLOCK_CYCLE           = 1
COMM_TIMEOUT          = 100
COMM_INTERRUPT_CYCLES = 10

PROCESSOR_POWERS = [
    0.51,
    0.5,
    0.49,
    # 0.3
]

DAGS = [
    "1fat0.3100.dag",
    # "2fat0.3100.dag",
    # "3fat0.3100.dag",
    # "4fat0.3100.dag",
    # "5fat0.3100.dag",
    # "1fat0.9100.dag",
    # "2fat0.9100.dag",
    # "3fat0.9100.dag",
    # "4fat0.9100.dag",
    # "5fat0.9100.dag",
    # "1300.dag",
    # "2300.dag",
    # "3300.dag",
    # "4300.dag",
    # "5300.dag",
    # "6300.dag",
]


text:
"""
    checkMethodsDefined(GenJ)

documentation
"""
function checkParametersDefined(GenJ::GeneticJulia, index::Int = 0)

    if index == 0
        if !isdefined(GenJ._experimentInfo, :_individualType)
            error("The experiment individualType must be set")
        end

        if !isdefined(GenJ._experimentInfo, :_algorithm)
            error("The experiment algorithm must be set")
        end

        if !isdefined(GenJ, :_evaluator)
            error("The evaluator must be set")
        end

        if !isdefined(GenJ, :_generator)
            error("The generator must be set")
        end

        if !isdefined(GenJ, :_selector)
            error("The selector must be set")
        end

        if !isdefined(GenJ, :_crossoverOp)
            error("The crossover operator must be set")
        end

        if !isdefined(GenJ, :_mutationOp)
            error("The mutation operator must be set")
        end

        if !isdefined(GenJ, :_replacementOp)
            error("The replacement operator must be set")
        end
    else
        if !isdefined(GenJ._experimentInfo, :_individualType)
            error("The experiment individualType must be set " * string(index))
        end

        if !isdefined(GenJ._experimentInfo, :_algorithm)
            error("The experiment algorithm must be set in experiment " * string(index))
        end

        if !isdefined(GenJ, :_evaluator)
            error("The evaluator must be set in experiment " * string(index))
        end

        if !isdefined(GenJ, :_generator)
            error("The generator must be set in experiment " * string(index))
        end

        if !isdefined(GenJ, :_selector)
            error("The selector must be set in experiment " * string(index))
        end

        if !isdefined(GenJ, :_crossoverOp)
            error("The crossover operator must be set in experiment " * string(index))
        end

        if !isdefined(GenJ, :_mutationOp)
            error("The mutation operator must be set in experiment " * string(index))
        end

        if !isdefined(GenJ, :_replacementOp)
            error("The replacement operator must be set in experiment " * string(index))
        end
    end
end # function

function checkParametersCorrectGP(GenJ::GeneticJulia, index::Int = 0)

    rng = Random.MersenneTwister(1234)
    if index == 0
        if !isdefined(GenJ._experimentInfo, :_GPExperimentInfo)
            error("Information about Genetic Programming problem has not been given for a GP problem ($(GenJ._experimentInfo._individualType)), this information must be specified in setGPInfo")
        end
    else
        if !isdefined(GenJ._experimentInfo, :_GPExperimentInfo)
            error("Information about Genetic Programming problem has not been given for a GP problem ($(GenJ._experimentInfo._individualType)), this information must be specified in setGPInfo in experiment $index")
        end
    end

    functionSet = GenJ._experimentInfo._GPExperimentInfo._functionSet
    terminalSet = GenJ._experimentInfo._GPExperimentInfo._terminalSet

    OneByOne = true
    All = true
    ind = 0

    try
        ind = GenJ._generator._method(functionSet, terminalSet, 3, rng, GenJ._generator._varArgs...)
    catch e
        if isa(e, MethodError)
            OneByOne = false
        end
    end

    try
        ind = GenJ._generator._method(1, functionSet, terminalSet, 3, rng, GenJ._generator._varArgs...)
    catch e
        if isa(e, MethodError)
            All = false
        end
    end

    if index == 0

        if !(All || OneByOne)
            error("Generation method ($(GenJ._generator._method)) introduced is not compatible with arguments specified")
        end

        ind = All ? GenJ._generator._method(1, functionSet, terminalSet, 3, rng, GenJ._generator._varArgs...) : GenJ._generator._method(functionSet, terminalSet, 3, rng, GenJ._generator._varArgs...)

        if !OneByOne && GenJ._generator._generateOneByOne
            error("Generation method ($(GenJ._generator._method)) receives the population size, however, generateOneByOne is set to true")
        end

        if typeof(ind) <: Array && OneByOne
            error("Generation method ($(GenJ._generator._method)) returns an array of individuals while generateOneByOne is set to true, it should return a single individual")
        end

        if !All && !GenJ._generator._generateOneByOne
            error("Generation method ($(GenJ._generator._method)) doesn't receive the population size, however, generateOneByOne is set to false")
        end

        if !(typeof(ind) <: Array) && All
            error("Generation method ($(GenJ._generator._method)) returns a single individual while generateOneByOne is set to false, it should return an array of individuals")
        end
    else

        if !(All || OneByOne)
            error("Generation method ($(GenJ._generator._method)) introduced is not compatible with arguments specified in experiment $index")
        end

        ind = All ? GenJ._generator._method(1, functionSet, terminalSet, 3, rng, GenJ._generator._varArgs...) : GenJ._generator._method(functionSet, terminalSet, 3, rng, GenJ._generator._varArgs...)

        if !OneByOne && GenJ._generator._generateOneByOne
            error("Generation method ($(GenJ._generator._method)) receives the population size, however, generateOneByOne is set to true in experiment $index")
        end

        if typeof(ind) <: Array && OneByOne
            error("Generation method ($(GenJ._generator._method)) returns an array of individuals while generateOneByOne is set to true, it should return a single individual in experiment $index")
        end

        if !All && !GenJ._generator._generateOneByOne
            error("Generation method ($(GenJ._generator._method)) doesn't receive the population size, however, generateOneByOne is set to false in experiment $index")
        end

        if !(typeof(ind) <: Array) && All
            error("Generation method ($(GenJ._generator._method)) returns a single individual while generateOneByOne is set to false, it should return an array of individuals in experiment $index")
        end
    end

    margin = 1e-6
    nParents = GenJ._crossoverOp._nParents
    nChildren = GenJ._crossoverOp._nChildren
    individuals = Array{Genotype}(undef, nParents)
    popSize = GenJ._generator._popSize
    k = 0
    fitness = 0

    nSelectedParents = getNSelectedParents(GenJ._selector)
    if !(typeof(nSelectedParents) <: Int)
        nSelectedParents = convert(Int, round(nSelectedParents * popSize))
        remainder = nSelectedParents % nParents

        if remainder != 0
            nSelectedParents = nSelectedParents + nParents - remainder
        end

        if !GenJ._selector._samplingWithRep && nSelectedParents > popSize
            nSelectedParents -= nParents
        end
    end

    if GenJ._generator._generateOneByOne
        for i=1:nParents
            individuals[i] = GenJ._generator._method(GenJ._experimentInfo._GPExperimentInfo._functionSet, GenJ._experimentInfo._GPExperimentInfo._terminalSet, 3, rng, GenJ._generator._varArgs...)
        end

    else
        individuals = GenJ._generator._method(nParents, GenJ._experimentInfo._GPExperimentInfo._functionSet, GenJ._experimentInfo._GPExperimentInfo._terminalSet, 3, rng, GenJ._generator._varArgs...)
    end

    if index == 0
        if GenJ._crossoverOp._probability < margin && GenJ._mutationOp._probability < margin
            error("Either crossover or mutation or both must have a probability above 0")
        end

        if typeof(individuals[1]) != GenJ._experimentInfo._individualType
            error("Generation method ($(GenJ._generator._method)) should generate individuals (you generate $(typeof(individuals[1]))) of the same type as the one specified in individualType ($(GenJ._experimentInfo._individualType))")
        end
#fitness = fitnessEval(GenJ._evaluator._fitnessFunctions[1], GenJ._experimentInfo._GPExperimentInfo, individuals[1])
        try
            len = length(GenJ._evaluator._fitnessFunctions)
            k = 1

            for i=1:len
                fitness = fitnessEval(GenJ._evaluator._fitnessFunctions[i], GenJ._experimentInfo._GPExperimentInfo, individuals[1])

                if !(typeof(fitness) <: Number)
                    throw(TypeError(fitness, ""))
                end
                k = k + 1
            end
        catch e

            if isa(e, MethodError)
                error("Fitness function ($(GenJ._evaluator._fitnessFunctions[k]._fitnessFunction)) is not compatible with the individual type ($(GenJ._experimentInfo._individualType))")
            elseif isa(e, TypeError)
                error("The fitness ($fitness) of an individual generated by ($(GenJ._evaluator._fitnessFunctions[k]._fitnessFunction)) must be numerical")
            end
        end

        if (nSelectedParents % nParents) != 0
            error("The number of selected parents ($nSelectedParents) is not divisible by number of parents per cross ($nParents)")
        end

        if !GenJ._selector._samplingWithRep && nSelectedParents > popSize
            error("Number of selected parents ($nSelectedParents) should not be greater than population size ($popSize) when selection method is without replacement")
        end

        try
            childCross = GenJ._crossoverOp._method(individuals..., 3, rng, GenJ._crossoverOp._varArgs...)
            len = length(childCross)
            if len != GenJ._crossoverOp._nChildren
                error("Children generated by crossover ($len) must match the number of children given as argument ($(GenJ._crossoverOp._nChildren))")
            end
            typeassert(childCross[1], typeof(individuals[1]))
        catch e
            if isa(e, MethodError)
                error("The crossover method ($(GenJ._crossoverOp._method)) is not compatible with the individual type generated ($(GenJ._experimentInfo._individualType))")
            elseif isa(e, TypeError)
                error("The type of the child ($(typeof(childCross[1]))) generated by crossover does not match the type of its parents ($(GenJ._experimentInfo._individualType))")
            end
        end

        try
            childMut = GenJ._mutationOp._method(individuals[1], 3, GenJ._generator._functionSet, GenJ._generator._terminalSet, rng, GenJ._mutationOp._varArgs...)
            typeassert(childMut, typeof(individuals[1]))
        catch e
            if isa(e, MethodError)
                error("The mutation method ($(GenJ._mutationOp._method)) is not compatible with the individual type generated ($(GenJ._experimentInfo._individualType))")
            elseif isa(e, TypeError)
                error("The type of the child ($(typeof(childMut))) generated by mutation does not match the type of its parent ($(GenJ._experimentInfo._individualType))")
            end
        end

    else
        if GenJ._crossoverOp._probability < margin && GenJ._mutationOp._probability < margin
            error("Either crossover or mutation or both must have a probability above 0 in experiment $index")
        end

        if typeof(individuals[1]) != GenJ._experimentInfo._individualType
            error("Generation method ($(GenJ._generator._method)) should generate individuals (you generate $(typeof(individuals[1]))) of the same type as the one specified in individualType ($(GenJ._experimentInfo._individualType)) in experiment $index")
        end

        try
            len = length(GenJ._evaluator._fitnessFunctions)
            k = 1

            for i=1:len
                fitness = fitnessEval(GenJ._evaluator._fitnessFunctions[i], GenJ._experimentInfo._GPExperimentInfo, individuals[1])

                if !(typeof(fitness) <: Number)
                    throw(TypeError(fitness, ""))
                end
                k = k + 1
            end
        catch e
            if isa(e, MethodError)
                error("Fitness function ($(GenJ._evaluator._fitnessFunctions[k]._fitnessFunction)) is not compatible with the individual type ($(GenJ._experimentInfo._individualType)) in experiment $index")
            elseif isa(e, TypeError)
                error("The fitness ($fitness) of an individual generated by ($(GenJ._evaluator._fitnessFunctions[k]._fitnessFunction)) must be numerical in experiment $index")
            end
        end

        if (nSelectedParents % nParents) != 0
            error("The number of selected parents ($nSelectedParents) is not divisible by number of parents per cross ($nParents) in experiment $index")
        end

        if !GenJ._selector._samplingWithRep && nSelectedParents > popSize
            error("Number of selected parents ($nSelectedParents) should not be greater than population size ($popSize) when selection method is without replacement in experiment $index")
        end

        try
            childCross = GenJ._crossoverOp._method(individuals..., 3, rng, GenJ._crossoverOp._varArgs...)
            len = length(childCross)
            if len != GenJ._crossoverOp._nChildren
                error("Children generated by crossover ($len) must match the number of children given as argument ($(GenJ._crossoverOp._nChildren)) in experiment $index")
            end
            typeassert(childCross[1], typeof(individuals[1]))
        catch e
            if isa(e, MethodError)
                error("The crossover method ($(GenJ._crossoverOp._method)) is not compatible with the individual type generated ($(GenJ._experimentInfo._individualType)) in experiment $index")
            elseif isa(e, TypeError)
                error("The type of the child ($(typeof(childCross[1]))) generated by crossover does not match the type of its parents ($(GenJ._experimentInfo._individualType)) in experiment $index")
            end
        end

        try
            childMut = GenJ._mutationOp._method(individuals[1], 3, GenJ._generator._functionSet, GenJ._generator._terminalSet, rng, GenJ._mutationOp._varArgs...)
            typeassert(childMut, typeof(individuals[1]))
        catch e
            if isa(e, MethodError)
                error("The mutation method ($(GenJ._mutationOp._method)) is not compatible with the individual type generated ($(GenJ._experimentInfo._individualType)) in experiment $index")
            elseif isa(e, TypeError)
                error("The type of the child ($(typeof(childMut))) generated by mutation does not match the type of its parent ($(GenJ._experimentInfo._individualType)) in experiment $index")
            end
        end

    end

end

function checkParametersCorrectGA(GenJ::GeneticJulia, index::Int = 0)
    rng = Random.MersenneTwister(1234)
    OneByOne = true
    All = true
    ind = 0

    if index == 0
        if isdefined( GenJ._experimentInfo, :_GPExperimentInfo)
            @warn "Information about Genetic Programming problem has been given for a non-GP problem ($(GenJ._experimentInfo._individualType)), this information will be dismissed"
        end
    else
        if isdefined( GenJ._experimentInfo, :_GPExperimentInfo)
            @warn "Information about Genetic Programming problem has been given for a non-GP problem ($(GenJ._experimentInfo._individualType)), this information will be dismissed in experiment $index"
        end
    end

    try
        ind = GenJ._generator._method(rng, GenJ._generator._varArgs...)
    catch e
        if isa(e, MethodError)
            OneByOne = false
        end
    end

    try
        ind = GenJ._generator._method(1, rng, GenJ._generator._varArgs...)
    catch e
        if isa(e, MethodError)
            All = false
        end
    end

    if index == 0

        if !(All || OneByOne)
            error("Generation method ($(GenJ._generator._method)) introduced is not compatible with arguments specified")
        end

        ind = All ? GenJ._generator._method(1, rng, GenJ._generator._varArgs...) : GenJ._generator._method(rng, GenJ._generator._varArgs...)

        if !OneByOne && GenJ._generator._generateOneByOne
            error("Generation method ($(GenJ._generator._method)) receives the population size, however, generateOneByOne is set to true")
        end

        if typeof(ind) <: Array && OneByOne
            error("Generation method ($(GenJ._generator._method)) returns an array of individuals while generateOneByOne is set to true, it should return a single individual")
        end

        if !All && !GenJ._generator._generateOneByOne
            error("Generation method ($(GenJ._generator._method)) doesn't receive the population size, however, generateOneByOne is set to false")
        end

        if !(typeof(ind) <: Array) && All
            error("Generation method ($(GenJ._generator._method)) returns a single individual while generateOneByOne is set to false, it should return an array of individuals")
        end
    else

        if !(All || OneByOne)
            error("Generation method ($(GenJ._generator._method)) introduced is not compatible with arguments specified in experiment $index")
        end

        ind = All ? GenJ._generator._method(1, rng, GenJ._generator._varArgs...) : GenJ._generator._method(rng, GenJ._generator._varArgs...)

        if !OneByOne && GenJ._generator._generateOneByOne
            error("Generation method ($(GenJ._generator._method)) receives the population size, however, generateOneByOne is set to true in experiment $index")
        end

        if typeof(ind) <: Array && OneByOne
            error("Generation method ($(GenJ._generator._method)) returns an array of individuals while generateOneByOne is set to true, it should return a single individual in experiment $index")
        end

        if !All && !GenJ._generator._generateOneByOne
            error("Generation method ($(GenJ._generator._method)) doesn't receive the population size, however, generateOneByOne is set to false in experiment $index")
        end

        if !(typeof(ind) <: Array) && All
            error("Generation method ($(GenJ._generator._method)) returns a single individual while generateOneByOne is set to false, it should return an array of individuals in experiment $index")
        end
    end

    margin = 1e-6
    nParents = GenJ._crossoverOp._nParents
    nChildren = GenJ._crossoverOp._nChildren
    individuals = Array{Genotype}(undef, nParents)
    popSize = GenJ._generator._popSize
    fitness = 0
    k = 0

    nSelectedParents = getNSelectedParents(GenJ._selector)
    if !(typeof(nSelectedParents) <: Int)
        nSelectedParents = convert(Int, round(nSelectedParents * popSize))
        remainder = nSelectedParents % nParents

        if remainder != 0
            nSelectedParents = nSelectedParents + nParents - remainder
        end

        if !GenJ._selector._samplingWithRep && nSelectedParents > popSize
            nSelectedParents -= nParents
        end
    end

    if GenJ._generator._generateOneByOne
        for i=1:nParents
            individuals[i] = GenJ._generator._method(rng, GenJ._generator._varArgs...)
        end

    else
        individuals = GenJ._generator._method(nParents, rng, GenJ._generator._varArgs...)
    end

    if index == 0
        if GenJ._crossoverOp._probability < margin && GenJ._mutationOp._probability < margin
            error("Either crossover or mutation or both must have a probability above 0")
        end

        if typeof(individuals[1]) != GenJ._experimentInfo._individualType
            error("Generation method ($(GenJ._generator._method)) should generate individuals (you generate $(typeof(individuals[1]))) of the same type as the one specified in individualType ($(GenJ._experimentInfo._individualType))")
        end
fitnessEval(GenJ._evaluator._fitnessFunctions[1], individuals[1])
        try
            len = length(GenJ._evaluator._fitnessFunctions)
            k = 1

            for i=1:len
                fitness = fitnessEval(GenJ._evaluator._fitnessFunctions[i], individuals[1])

                if !(typeof(fitness) <: Number)
                    throw(TypeError(fitness, ""))
                end
                k = k + 1
            end
        catch e
            if isa(e, MethodError)
                error("Fitness function ($(GenJ._evaluator._fitnessFunctions[k]._fitnessFunction)) is not compatible with the individual type ($(GenJ._experimentInfo._individualType))")
            elseif isa(e, TypeError)
                error("The fitness ($fitness) of an individual generated by ($(GenJ._evaluator._fitnessFunctions[k]._fitnessFunction)) must be numerical")
            end
        end

        if (nSelectedParents % nParents) != 0
            error("The number of selected parents ($nSelectedParents) is not divisible by number of parents per cross ($nParents)")
        end

        if !GenJ._selector._samplingWithRep && nSelectedParents > popSize
            error("Number of selected parents ($nSelectedParents) should not be greater than population size ($popSize) when selection method is without replacement")
        end

        try
            childCross = GenJ._crossoverOp._method(individuals..., rng, GenJ._crossoverOp._varArgs...)
            len = length(childCross)
            if len != GenJ._crossoverOp._nChildren
                error("Children generated by crossover ($len) must match the number of children given as argument ($(GenJ._crossoverOp._nChildren))")
            end
            typeassert(childCross[1], typeof(individuals[1]))
        catch e
            if isa(e, MethodError)
                error("The crossover method ($(GenJ._crossoverOp._method)) is not compatible with the individual type generated ($(GenJ._experimentInfo._individualType))")
            elseif isa(e, TypeError)
                error("The type of the child ($(typeof(childCross[1]))) generated by crossover does not match the type of its parents ($(GenJ._experimentInfo._individualType))")
            end
        end

        try
            childMut = GenJ._mutationOp._method(individuals[1], rng, GenJ._mutationOp._varArgs...)
            typeassert(childMut, typeof(individuals[1]))
        catch e
            if isa(e, MethodError)
                error("The mutation method ($(GenJ._mutationOp._method)) is not compatible with the individual type generated ($(GenJ._experimentInfo._individualType))")
            elseif isa(e, TypeError)
                error("The type of the child ($(typeof(childMut))) generated by mutation does not match the type of its parent ($(GenJ._experimentInfo._individualType))")
            end
        end
    else

        if GenJ._crossoverOp._probability < margin && GenJ._mutationOp._probability < margin
            error("Either crossover or mutation or both must have a probability above 0 in experiment $index")
        end

        if typeof(individuals[1]) != GenJ._experimentInfo._individualType
            error("Generation method ($(GenJ._generator._method)) should generate individuals (you generate $(typeof(individuals[1]))) of the same type as the one specified in individualType ($(GenJ._experimentInfo._individualType)) in experiment $index")
        end

        try
            len = length(GenJ._evaluator._fitnessFunctions)
            k = 1

            for i=1:len
                fitness = fitnessEval(GenJ._evaluator._fitnessFunctions[i], individuals[1])

                if !(typeof(fitness) <: Number)
                    throw(TypeError(fitness, ""))
                end
                k = k + 1
            end
        catch e
            if isa(e, MethodError)
                error("Fitness function ($(GenJ._evaluator._fitnessFunctions[k]._fitnessFunction)) is not compatible with the individual type ($(GenJ._experimentInfo._individualType)) in experiment $index")
            elseif isa(e, TypeError)
                error("The fitness ($fitness) of an individual generated by ($(GenJ._evaluator._fitnessFunctions[k]._fitnessFunction)) must be numerical in experiment $index")
            end
        end

        if (nSelectedParents % nParents) != 0
            error("The number of selected parents ($nSelectedParents) is not divisible by number of parents per cross ($nParents) in experiment $index")
        end

        if !GenJ._selector._samplingWithRep && nSelectedParents > popSize
            error("Number of selected parents ($nSelectedParents) should not be greater than population size ($popSize) when selection method is without replacement in experiment $index")
        end

        try
            childCross = GenJ._crossoverOp._method(individuals..., rng, GenJ._crossoverOp._varArgs...)
            len = length(childCross)
            if len != GenJ._crossoverOp._nChildren
                error("Children generated by crossover ($len) must match the number of children given as argument ($(GenJ._crossoverOp._nChildren)) in experiment $index")
            end
            typeassert(childCross[1], typeof(individuals[1]))
        catch e
            if isa(e, MethodError)
                error("The crossover method ($(GenJ._crossoverOp._method)) is not compatible with the individual type generated ($(GenJ._experimentInfo._individualType)) in experiment $index")
            elseif isa(e, TypeError)
                error("The type of the child ($(typeof(childCross[1]))) generated by crossover does not match the type of its parents ($(GenJ._experimentInfo._individualType)) in experiment $index")
            end
        end

        try
            childMut = GenJ._mutationOp._method(individuals[1], rng, GenJ._mutationOp._varArgs...)
            typeassert(childMut, typeof(individuals[1]))
        catch e
            if isa(e, MethodError)
                error("The mutation method ($(GenJ._mutationOp._method)) is not compatible with the individual type generated ($(GenJ._experimentInfo._individualType)) in experiment $index")
            elseif isa(e, TypeError)
                error("The type of the child ($(typeof(childMut))) generated by mutation does not match the type of its parent ($(GenJ._experimentInfo._individualType)) in experiment $index")
            end
        end
    end
end

"""
    checkParameters(GenJ)

Checks that everything has been set correctly before running an experiment.
"""
function checkParametersCorrect(GenJ::GeneticJulia, index::Int = 0)

    if GenJ._experimentInfo._individualType <: GPGenotype
        checkParametersCorrectGP(GenJ, index)
    else
        checkParametersCorrectGA(GenJ, index)
    end

end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
struct SymmetricTensor{T}
    xx::Matrix{T}
    zz::Matrix{T}
    xz::Matrix{T}
end

struct Gradient{T}
    ∂x::Matrix{T}
    ∂z::Matrix{T}
end

function initstress(nel; nip = 6)
    F   = [@SMatrix [1.0 0.0; 0.0 1.0]  for _ in 1:nel, _ in 1:nip]
    F0  = deepcopy(F)
    τ   = SymmetricTensor(fill(0.0,nel,nip),fill(0.0,nel,nip),fill(0.0,nel,nip))
    ε   = SymmetricTensor(fill(0.0,nel,nip),fill(0.0,nel,nip),fill(0.0,nel,nip))
    τII = fill(0.0,nel,nip)
    εII = fill(0.0,nel,nip)
    ∇T  = Gradient(fill(0.0,nel,nip), fill(0.0,nel,nip))
    return F, F0, τ, ε, τII, εII, ∇T
end

function _stress!(
    # F,
    ε,
    εII,
    U,
    e2n,
    iel,
    DoF_U,
    θ,
    r,
    SF_Stress::ShapeFunctionsStress,
    # Δt
)

    # unpack shape functions
    ∇N, dN3ds, N3 = SF_Stress.∇N, SF_Stress.dN3ds, SF_Stress.N3

    # Polar coordinates of element nodes
    θ_el = @SVector [θ[i, iel] for i in 1:3]
    r_el = @SVector [r[i, iel] for i in 1:3]
    coords = SMatrix{3,2}([θ_el r_el]) # could directly allocate just this one

    # Jacobian n. 1 (p:=polar, l:=local): reference element --> current element
    J_pl = dN3ds * coords
    # invdetJ_pl = 1/mydet(J_pl)
    invdetJ_pl = 1/det(J_pl)

    # the Jacobian ∂ξ∂θ to transform local (ξ, η) into global (θ,r) derivatives
    #     ∂ξ∂θ = [ R_31    -R_21
    #             -Th_31    Th_21] / detJa_PL
    R_21 = r_el[2] - r_el[1]  # = -detJa_PL*deta_dth
    R_31 = r_el[3] - r_el[1]  # =  detJa_PL*dxi_dth
    Th_31 = θ_el[3] - θ_el[1] # = -detJa_PL*dxi_dr
    Th_21 = θ_el[2] - θ_el[1] # =  detJa_PL*deta_dr

    Udofs = DoF_U[iel]
    Uel_x = @SVector [
        U[Udofs[1]], U[Udofs[3]],  U[Udofs[5]],  U[Udofs[7]], U[Udofs[9]], U[Udofs[11]]
    ]
    Uel_z = @SVector [
        U[Udofs[2]], U[Udofs[4]],  U[Udofs[6]], U[Udofs[8]], U[Udofs[10]], U[Udofs[12]]
    ]

    # INTEGRATION LOOP
    # onesixth = 1/6
    @inbounds for ip in 1:7

        # Unpack shape functions
        N3_ip = N3[ip]
        ∇N_ip = ∇N[ip]

        # ρ at ith integration point
#         η_ip = mydot(SVector{3}([η[e2n[i, iel]] for i in 1:3]), N3_ip)

        # Polar coordinates of the integration points
        θ_ip = mydot(θ_el, N3_ip)
        r_ip = mydot(r_el, N3_ip)
        invr_ip = 1/r_ip
        sin_ip, cos_ip = sincos(θ_ip)
        cos_ip_r_ip = cos_ip * invr_ip
        sin_ip_r_ip = sin_ip * invr_ip

        # Build inverse of the 2nd Jacobian
        invJ_double = @SMatrix [
            @muladd( R_31*cos_ip_r_ip-Th_31*sin_ip)   @muladd(-R_21*cos_ip_r_ip+Th_21*sin_ip)
            @muladd(-R_31*sin_ip_r_ip-Th_31*cos_ip)   @muladd( R_21*sin_ip_r_ip+Th_21*cos_ip)
        ]

        # Partial derivatives and strain rate components
        ∂N∂x = invJ_double * ∇N_ip * invdetJ_pl
        εxx = dot(view(∂N∂x,1,:), Uel_x) # = ∂Ux∂x
        εzz = dot(view(∂N∂x,2,:), Uel_z) # = ∂Uz∂z
        ∂Ux∂z = dot(view(∂N∂x,2,:), Uel_x)
        ∂Uz∂x = dot(view(∂N∂x,1,:), Uel_z)
        εxz = 0.5*(∂Ux∂z+∂Uz∂x)

        # update stain rate tensor
        ε.xx[iel, ip] = εxx
        ε.zz[iel, ip] = εzz
        ε.xz[iel, ip] = εxz

        # strain rate second invariant
        εII[iel, ip] =  √( @muladd(0.5*(εxx*εxx + εzz*εzz) + εxz*εxz))

#         # deviatoric stress
#         τxx = η_ip * ( 4*εxx/3 + -2*εzz/3)
#         τzz = η_ip * (-2*εxx/3 +  4*εzz/3)
#         τxz = η_ip *   2*εxz

#         τ.xx[iel, ip] = τxx
#         τ.zz[iel, ip] = τzz
#         τ.xz[iel, ip] = τxz

        # transpose of the velocity gradient
        # ∇Uᵀ = @SMatrix [
        #     ∂Ux∂x ∂Ux∂z
        #     ∂Uz∂x ∂Uz∂z
        # ]
        # Δt∇Uᵀ= Δt*∇Uᵀ

        # # Integrate and update gradient of deformation tensor F
        # F0 = F[iel, ip]
        # k1 = Δt∇Uᵀ * F0
        # Fi = k1*0.5 .+ F0
        # k2 = Δt∇Uᵀ * Fi
        # Fi = k2*0.5 .+ F0
        # k3 = Δt∇Uᵀ * Fi
        # Fi = k3 .+ F0
        # k4 = Δt∇Uᵀ * Fi
        # F[iel, ip] = F0 + (k1 + 2*(k2 + k3) + k4)*onesixth

    end

end

function stress!(ε, εII, U, e2n, nel, DoF_U, coordinates, SF_Stress)
    @batch per=core for iel in 1:nel
        _stress!(ε, εII, U, e2n, iel, DoF_U, coordinates.θ, coordinates.r, SF_Stress)
        # _stress!(F, ε, εII, τ, η.node, U, e2n, iel, DoF_U, coordinates.θ, coordinates.r, SF_Stress, Δt)
    end
end

function _deformation_gradient!(
    F,
    U,
    iel,
    DoF_U,
    θ,
    r,
    SF_Stress::ShapeFunctionsStress,
    Δt
)

    # unpack shape functions
    ∇N, dN3ds, N3 = SF_Stress.∇N, SF_Stress.dN3ds, SF_Stress.N3

    # Polar coordinates of element nodes
    θ_el = @SVector [θ[i, iel] for i in 1:3]
    r_el = @SVector [r[i, iel] for i in 1:3]
    coords = SMatrix{3,2}([θ_el r_el]) # could directly allocate just this one

    # Jacobian n. 1 (p:=polar, l:=local): reference element -> current element
    J_pl = dN3ds * coords
    invdetJ_pl = 1/det(J_pl)

    # the Jacobian ∂ξ∂θ to transform local (ξ, η) into global (θ,r) derivatives
    #     ∂ξ∂θ = [ R_31    -R_21
    #             -Th_31    Th_21] / detJa_PL
    R_21 = r_el[2] - r_el[1]  # = -detJa_PL*deta_dth
    R_31 = r_el[3] - r_el[1]  # =  detJa_PL*dxi_dth
    Th_31 = θ_el[3] - θ_el[1] # = -detJa_PL*dxi_dr
    Th_21 = θ_el[2] - θ_el[1] # =  detJa_PL*deta_dr

    Udofs = DoF_U[iel]
    Uel_x = @SVector [
        U[Udofs[1]], U[Udofs[3]],  U[Udofs[5]],  U[Udofs[7]], U[Udofs[9]], U[Udofs[11]]
    ]
    Uel_z = @SVector [
        U[Udofs[2]], U[Udofs[4]],  U[Udofs[6]], U[Udofs[8]], U[Udofs[10]], U[Udofs[12]]
    ]

    # INTEGRATION LOOP
    onesixth = 1/6
    @inbounds for ip in 1:7

        # Unpack shape functions
        N3_ip = N3[ip]
        ∇N_ip = ∇N[ip]

        # Polar coordinates of the integration points
        θ_ip = mydot(θ_el, N3_ip)
        r_ip = mydot(r_el, N3_ip)
        invr_ip = 1/r_ip
        sin_ip, cos_ip = sincos(θ_ip)
        cos_ip_r_ip = cos_ip * invr_ip
        sin_ip_r_ip = sin_ip * invr_ip

        # Build inverse of the 2nd Jacobian
        invJ_double = @SMatrix [
            @muladd( R_31*cos_ip_r_ip - Th_31*sin_ip)   @muladd(-R_21*cos_ip_r_ip + Th_21*sin_ip)
            @muladd(-R_31*sin_ip_r_ip - Th_31*cos_ip)   @muladd( R_21*sin_ip_r_ip + Th_21*cos_ip)
        ]

        # Partial derivatives and strain rate components
        ∂N∂x = invJ_double * ∇N_ip * invdetJ_pl
        ∂Ux∂x = dot(view(∂N∂x,1,:), Uel_x)
        ∂Uz∂z = dot(view(∂N∂x,2,:), Uel_z)
        ∂Ux∂z = dot(view(∂N∂x,2,:), Uel_x)
        ∂Uz∂x = dot(view(∂N∂x,1,:), Uel_z)

        # transpose of the velocity gradient
        ∇Uᵀ = @SMatrix [
            ∂Ux∂x ∂Ux∂z
            ∂Uz∂x ∂Uz∂z
        ]

        # precompute velocity gradient × Δt
        Δt∇Uᵀ= Δt*∇Uᵀ

        # Integrate and update gradient of deformation tensor F
        F0 = F[iel, ip]
        k1 = Δt∇Uᵀ * F0
        Fi = k1*0.5 .+ F0
        k2 = Δt∇Uᵀ * Fi
        Fi = k2*0.5 .+ F0
        k3 = Δt∇Uᵀ * Fi
        Fi = k3 .+ F0
        k4 = Δt∇Uᵀ * Fi
        F[iel, ip] = @muladd(F0 + (k1 + 2*(k2 + k3) + k4)*onesixth)

    end

end

function deformation_gradient!(F, U, nel, DoF_U, coordinates, SF_Stress, Δt)
    @batch per=core for iel in 1:nel
        _deformation_gradient!(F, U, iel, DoF_U, coordinates.θ, coordinates.r, SF_Stress, Δt)
    end
end

secondinvariant(xx,zz,xz) = sqrt( 0.5*(xx^2 + zz^2) + xz^2 )

@inline function secondinvariant(A::SymmetricTensor)
    II  = Matrix{Float64}(undef, size(A.xx))
    xx, zz, xz = A.xx, A.zz, A.xz
    @batch for i in eachindex(xx)
         II[i] = √( @muladd(0.5*(xx[i]*xx[i] + zz[i]*zz[i]) + xz[i]*xz[i]))
    end
    return II
end

@inline function secondinvariant!(II, A::SymmetricTensor)
    xx, zz, xz = A.xx, A.zz, A.xz
    @batch for i in eachindex(xx)
        @inbounds II[i] = √( @muladd(0.5*(xx[i]*xx[i] + zz[i]*zz[i]) + xz[i]*xz[i]))
    end
end

function shearheating(τ::SymmetricTensor{T}, ε::SymmetricTensor{T}) where {T}
    shear_heating = fill(0.0, size(τ.xx))
    # trace = @. (ε.xx + ε.zz)/3
    @avx for i in CartesianIndices(shear_heating) # faster than threading for current matrices size
        shear_heating[i] = τ.xx[i]*(ε.xx[i]- (ε.xx[i] + ε.zz[i])/3) +
                           τ.zz[i]*(ε.zz[i]- (ε.xx[i] + ε.zz[i])/3) +
                           2τ.xz[i]*ε.xz[i]
    end
    shear_heating
end

function rowdot!(C,A,B)
    @inbounds for i in 1:size(A,1)
        C[i] = mydot(view(A,i,:),view(B,i,:))
    end
end

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
struct JohnsonSU{T1<:Real,T2<:Real,T3<:Real,T4<:Real} <: Distributions.Distribution{Univariate,Continuous}
    gamma::T1 #positive: left tail, negativ: right tail, 0: no tail
    delta::T2 #height and width
    xi::T3 #shift along the x-axis
    lambda::T4 #height and width
end

function JohnsonSU(; gamma::Real = 10, delta::Real = 3.5, xi::Real = 10, lambda::Real = 1)
    gamma, delta, xi, lambda = promote(gamma, delta, xi, lambda)
    JohnsonSU(gamma, delta, xi, lambda)
end

Base.minimum(d::JohnsonSU) = -Inf
Base.maximum(d::JohnsonSU) = +Inf

StatsBase.params(d::JohnsonSU) = (d.gamma, d.delta, d.xi, d.lambda)
@inline Distributions.partype(d::JohnsonSU{T}) where {T<:Real} = T

Distributions.location(d::JohnsonSU) = mean(d)
Distributions.scale(d::JohnsonSU) = var(d)

Statistics.mean(d::JohnsonSU) = d.xi - d.lambda * exp(d.delta^(-2)/2)*sinh(d.gamma/d.delta)
Statistics.median(d::JohnsonSU) = d.xi + d.lambda * sinh(-d.gamma/d.delta)
Statistics.var(d::JohnsonSU) = (d.lambda^2)/2 * (exp(d.delta^(-2)) - 1) * (exp(d.delta^(-2))*cosh(2*d.gamma/d.delta) + 1)

function johnsontrafo(x::T1, gamma::T2, delta::T3, xi::T4, lambda::T5) where {T1<:Real,T2<:Real,T3<:Real,T4<:Real, T5<:Real}
    R = float(promote_type(T1, T2, T3, T4, T5))
    convert(R, (gamma + delta*asinh((x - xi)/lambda)))
end

function johnsontrafo_inv(x::T1, gamma::T2, delta::T3, xi::T4, lambda::T5) where {T1<:Real,T2<:Real,T3<:Real,T4<:Real, T5<:Real}
    R = float(promote_type(T1, T2, T3, T4, T5))
    convert(R, lambda*sinh((x - gamma)/delta) + xi)
end

function deriv_johnsontrafo(x::T1, gamma::T2, delta::T3, xi::T4, lambda::T5) where {T1<:Real,T2<:Real,T3<:Real,T4<:Real, T5<:Real}
    R = float(promote_type(T1, T2, T3, T4, T5))
    convert(R, (delta/lambda) * (1/(sqrt(1 + ((x - xi)/lambda)^2))))
end

function deriv_johnsontrafo_inv(x::T1, gamma::T2, delta::T3, xi::T4, lambda::T5) where {T1<:Real,T2<:Real,T3<:Real,T4<:Real, T5<:Real}
    R = float(promote_type(T1, T2, T3, T4, T5))
    convert(R, lambda*cosh((x - gamma)/delta)/delta)
end

function johnsontrafo_ladj(x::T1, gamma::T2, delta::T3, xi::T4, lambda::T5) where {T1<:Real,T2<:Real,T3<:Real,T4<:Real, T5<:Real}
    R = float(promote_type(T1, T2, T3, T4, T5))
    convert(R, log(abs(deriv_johnsontrafo(x,gamma,delta,xi,lambda))))
end

function johnsontrafo_inv_ladj(x::T1, gamma::T2, delta::T3, xi::T4, lambda::T5) where {T1<:Real,T2<:Real,T3<:Real,T4<:Real, T5<:Real}
    R = float(promote_type(T1, T2, T3, T4, T5))
    convert(R, log(abs(deriv_johnsontrafo_inv(x,gamma,delta,xi,lambda))))
end

@with_kw struct JohnsonTrafo{T<:Union{Real,AbstractVector{<:Real}}} <: Function
    gamma::T = 10.0
    delta::T = 3.5
    xi::T = 10.0
    lambda::T = 1.0
end

Base.:(==)(a::JohnsonTrafo, b::JohnsonTrafo) = a.gamma == b.gamma && a.delta == b.delta && a.xi == b.xi && a.lambda == b.lambda
Base.isequal(a::JohnsonTrafo, b::JohnsonTrafo) = isequal(a.gamma, b.gamma) && isequal(a.delta, b.delta) && isequal(a.xi, b.xi) && isequal(a.lambda, b.lambda)
Base.hash(x::JohnsonTrafo, h::UInt) = hash(x.lambda, hash(x.xi, hash(x.delta, hash(x.gamma, hash(:JohnsonTrafo, hash(:EuclidianNormalizingFlows, h))))))

@functor JohnsonTrafo

(f::JohnsonTrafo)(x::Union{Real,AbstractVecOrMat{<:Real}}) = johnsontrafo.(x, f.gamma, f.delta, f.xi, f.lambda)

function with_logabsdet_jacobian(f::JohnsonTrafo, x::Union{Real,AbstractVecOrMat{<:Real}})
    y = f(x)
    ladjs = johnsontrafo_ladj.(x, f.gamma, f.delta, f.xi, f.lambda)
    (y, sum_ladjs(ladjs))
end

inverse(f::JohnsonTrafo) = JohnsonTrafoInv(f.gamma, f.delta, f.xi, f.lambda)

@with_kw struct JohnsonTrafoInv{T<:Union{Real,AbstractVector{<:Real}}} <: Function
    gamma::T = 10.0
    delta::T = 3.5
    xi::T = 10.0
    lambda::T = 1.0
end

@functor JohnsonTrafoInv

Base.:(==)(a::JohnsonTrafoInv, b::JohnsonTrafoInv) = a.gamma == b.gamma && a.delta == b.delta && a.xi == b.xi && a.lambda == b.lambda
Base.isequal(a::JohnsonTrafoInv, b::JohnsonTrafoInv) = isequal(a.gamma, b.gamma) && isequal(a.delta, b.delta) && isequal(a.xi, b.xi) && isequal(a.lambda, b.lambda)
Base.hash(x::JohnsonTrafoInv, h::UInt) = hash(x.lambda, hash(x.xi, hash(x.delta, hash(x.gamma, hash(:JohnsonTrafoInv, hash(:EuclidianNormalizingFlows, h))))))

(f::JohnsonTrafoInv)(x::Union{Real,AbstractVecOrMat{<:Real}}) = johnsontrafo_inv.(x, f.gamma, f.delta, f.xi, f.lambda)

function with_logabsdet_jacobian(f::JohnsonTrafoInv, x::Union{Real,AbstractVecOrMat{<:Real}})
    y = f(x)
    neg_ladjs = johnsontrafo_ladj.(y, f.gamma, f.delta, f.xi, f.lambda)
    (y, - sum_ladjs(neg_ladjs))
end

inverse(f::JohnsonTrafoInv) = JohnsonTrafo(f.gamma, f.delta, f.xi, f.lambda)

#johnsontrafo(d::JohnsonSU, x::Real) = (d.gamma + d.delta*asinh((x - d.xi)/d.lambda))
#johnsontrafo_inv(d::JohnsonSU, x::Real) = d.lambda*sinh((x - d.gamma)/d.delta) + d.xi
#deriv_johnsontrafo(d::JohnsonSU, x::Real) = (d.delta/d.lambda) * (1/(sqrt(1 + ((x - d.xi)/d.lambda)^2)))
#deriv_johnsontrafo_inv(d::JohnsonSU, x::Real) = d.lambda*cosh((x - d.gamma)/d.delta)/d.delta
#johnsontrafo_ladj(d::JohnsonSU, x::Real) = log(abs(deriv_johnsontrafo(d::JohnsonSU, x::Real)))
#johnsontrafo_inv_ladj(d::JohnsonSU, x::Real) = log(abs(deriv_johnsontrafo_inv(d::JohnsonSU, x::Real)))

#Distributions.pdf(d::JohnsonSU, x::Real) = deriv_johnsontrafo(d,x) * 1/sqrt(2*pi) * exp(-0.5 * johnsontrafo(d,x)^2)
Distributions.pdf(d::JohnsonSU, x::Real) = deriv_johnsontrafo(x, d.gamma, d.delta, d.xi, d.lambda) * pdf(Normal(), johnsontrafo(x, d.gamma, d.delta, d.xi, d.lambda))
Distributions.cdf(d::JohnsonSU, x::Real) = cdf(Normal(), johnsontrafo(x, d.gamma, d.delta, d.xi, d.lambda))
#Distributions.logpdf(d::JohnsonSU, x::Real) = log(deriv_johnsontrafo(d,x) * 1/sqrt(2*pi) * exp(-0.5 * johnsontrafo(d,x)^2))
Distributions.logpdf(d::JohnsonSU, x::Real) = log(deriv_johnsontrafo(x, d.gamma, d.delta, d.xi, d.lambda) * pdf(Normal(), johnsontrafo(x, d.gamma, d.delta, d.xi, d.lambda)))
Distributions.logcdf(d::JohnsonSU, x::Real) = logcdf(Normal(), johnsontrafo(x, d.gamma, d.delta, d.xi, d.lambda))
Distributions.ccdf(d::JohnsonSU, x::Real) = 1 - Distributions.cdf(d, x)
Distributions.logccdf(d::JohnsonSU, x::Real) = log(1 - Distributions.cdf(d, x))

Statistics.quantile(d::JohnsonSU, x::Real) = johnsontrafo_inv(quantile(Normal(),x), d.gamma, d.delta, d.xi, d.lambda)