text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -58.37772561872491 <= q <= 51.748312049639644)

                     add_NL_constraint(m, :(1/(1 + exp(-(-0.3976822916974827 + -0.11888104611683081*1/(1 + exp(-(-0.09559544092481564 + 0.23201220435158865*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + -0.9126138494455787*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + 0.24324641666113456*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + 0.1690500900766101*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + 0.20833731206579253*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))))))) + -0.8277675607937618*1/(1 + exp(-(-0.41653739657549194 + 0.7547021537813423*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + 0.5787262183714494*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + -0.11492508246438904*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + 0.30062983013792843*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + -0.1276250114011983*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))))))) + 0.37430250255255215*1/(1 + exp(-(0.9557400948867283 + -0.3128610884258398*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + 0.29518845150992634*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + -0.18188148649835556*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + 0.8330182900675651*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + -0.8804319372999063*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))))))) + 0.9159784725268092*1/(1 + exp(-(-0.32868009819829513 + 0.7041500634234112*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + -0.44147758504740464*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + -0.5937281428689025*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + -0.3198554528784734*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + -0.047330375221362964*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))))))) + -0.21239133839693736*1/(1 + exp(-(-0.6560670298284208 + 0.11579976545709414*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + 0.7227941593613254*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + 0.5041216754022764*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + -0.9501104307252297*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + 0.7033943599032555*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5])))))))))) + 1/(1 + exp(-(-0.962366690533889 + 0.21288820713294498*1/(1 + exp(-(-0.09559544092481564 + 0.23201220435158865*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + -0.9126138494455787*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + 0.24324641666113456*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + 0.1690500900766101*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + 0.20833731206579253*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))))))) + 0.04621026197884115*1/(1 + exp(-(-0.41653739657549194 + 0.7547021537813423*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + 0.5787262183714494*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + -0.11492508246438904*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + 0.30062983013792843*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + -0.1276250114011983*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))))))) + 0.3579305330052933*1/(1 + exp(-(0.9557400948867283 + -0.3128610884258398*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + 0.29518845150992634*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + -0.18188148649835556*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + 0.8330182900675651*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + -0.8804319372999063*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))))))) + -0.40247732796428703*1/(1 + exp(-(-0.32868009819829513 + 0.7041500634234112*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + -0.44147758504740464*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + -0.5937281428689025*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + -0.3198554528784734*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + -0.047330375221362964*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))))))) + -0.3937461182921762*1/(1 + exp(-(-0.6560670298284208 + 0.11579976545709414*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + 0.7227941593613254*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + 0.5041216754022764*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + -0.9501104307252297*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + 0.7033943599032555*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5])))))))))) + 1/(1 + exp(-(-0.9189915189234159 + -0.5125130010222394*1/(1 + exp(-(-0.09559544092481564 + 0.23201220435158865*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + -0.9126138494455787*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + 0.24324641666113456*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + 0.1690500900766101*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + 0.20833731206579253*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))))))) + -0.3502499398907637*1/(1 + exp(-(-0.41653739657549194 + 0.7547021537813423*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + 0.5787262183714494*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + -0.11492508246438904*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + 0.30062983013792843*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + -0.1276250114011983*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))))))) + -0.8274375080497247*1/(1 + exp(-(0.9557400948867283 + -0.3128610884258398*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + 0.29518845150992634*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + -0.18188148649835556*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + 0.8330182900675651*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + -0.8804319372999063*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))))))) + -0.7757421130805153*1/(1 + exp(-(-0.32868009819829513 + 0.7041500634234112*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + -0.44147758504740464*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + -0.5937281428689025*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + -0.3198554528784734*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + -0.047330375221362964*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))))))) + 0.6602598595684599*1/(1 + exp(-(-0.6560670298284208 + 0.11579976545709414*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + 0.7227941593613254*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + 0.5041216754022764*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + -0.9501104307252297*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + 0.7033943599032555*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5])))))))))) + 1/(1 + exp(-(-0.760456357412739 + 0.8974231726544617*1/(1 + exp(-(-0.09559544092481564 + 0.23201220435158865*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + -0.9126138494455787*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + 0.24324641666113456*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + 0.1690500900766101*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + 0.20833731206579253*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))))))) + -0.7589867374120325*1/(1 + exp(-(-0.41653739657549194 + 0.7547021537813423*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + 0.5787262183714494*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + -0.11492508246438904*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + 0.30062983013792843*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + -0.1276250114011983*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))))))) + -0.3014899851721653*1/(1 + exp(-(0.9557400948867283 + -0.3128610884258398*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + 0.29518845150992634*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + -0.18188148649835556*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + 0.8330182900675651*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + -0.8804319372999063*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))))))) + 0.03263132145111802*1/(1 + exp(-(-0.32868009819829513 + 0.7041500634234112*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + -0.44147758504740464*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + -0.5937281428689025*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + -0.3198554528784734*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + -0.047330375221362964*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))))))) + 0.8460059521417822*1/(1 + exp(-(-0.6560670298284208 + 0.11579976545709414*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + 0.7227941593613254*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + 0.5041216754022764*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + -0.9501104307252297*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + 0.7033943599032555*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5])))))))))) + 1/(1 + exp(-(-0.1364886291423204 + -0.7689732133011344*1/(1 + exp(-(-0.09559544092481564 + 0.23201220435158865*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + -0.9126138494455787*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + 0.24324641666113456*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + 0.1690500900766101*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + 0.20833731206579253*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))))))) + 0.9275671307049969*1/(1 + exp(-(-0.41653739657549194 + 0.7547021537813423*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + 0.5787262183714494*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + -0.11492508246438904*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + 0.30062983013792843*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + -0.1276250114011983*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))))))) + -0.311285022430249*1/(1 + exp(-(0.9557400948867283 + -0.3128610884258398*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + 0.29518845150992634*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + -0.18188148649835556*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + 0.8330182900675651*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + -0.8804319372999063*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))))))) + -0.2373305013184268*1/(1 + exp(-(-0.32868009819829513 + 0.7041500634234112*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + -0.44147758504740464*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + -0.5937281428689025*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + -0.3198554528784734*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + -0.047330375221362964*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(x[5]))))))) + 0.7676258102518863*1/(1 + exp(-(-0.6560670298284208 + 0.11579976545709414*1/(1 + exp(-(0.22080335056838862 + -0.11207096478827117*$(x[1]) + 0.9802919969431398*$(x[2]) + -0.17783835331857878*$(x[3]) + 0.5553411614422452*$(x[4]) + 0.8384080073044262*$(x[5])))) + 0.7227941593613254*1/(1 + exp(-(-0.03619259662306584 + -0.3774088093292831*$(x[1]) + 0.859114754276578*$(x[2]) + -0.4844223734792257*$(x[3]) + 0.08654340008535533*$(x[4]) + -0.46594229223618644*$(x[5])))) + 0.5041216754022764*1/(1 + exp(-(-0.05960623749041094 + -0.4311407494411239*$(x[1]) + -0.025139997052890273*$(x[2]) + 0.060401892548068936*$(x[3]) + 0.5875386419954256*$(x[4]) + 0.5553096462934759*$(x[5])))) + -0.9501104307252297*1/(1 + exp(-(0.12047524138880839 + -0.4988789730921872*$(x[1]) + -0.41177595995379823*$(x[2]) + -0.6913890776083331*$(x[3]) + 0.032955544793602076*$(x[4]) + -0.5657303188428733*$(x[5])))) + 0.7033943599032555*1/(1 + exp(-(-0.14422850857319736 + 0.5362813944125526*$(x[1]) + 0.29620144485119626*$(x[2]) + -0.18264867223264059*$(x[3]) + 0.5265937438977235*$(x[4]) + 0.6008528580664168*$(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:
using NetCDF, GeoMakie
const Node = Observable

# Open the NetCDF file
f = NetCDF.open("/Users/anshul/Downloads/CERES_EBAF-TOA_Ed4.1_Subset_200003-201910.nc"; readdimvar = true)

lat = f.dim["lat"].vals
lon = f.dim["lon"].vals

t = 1 # time
f.vars["toa_net_all_mon"][:, :, t][:, :, 1] # this hack seems necessary for NetCDF, but I'm not sure why

field = Node{Matrix{Float64}}(f.vars["toa_net_all_mon"][:, :, t][:, :, 1])

cf = Node(copy(field[]))

source = Projection("+proj=lonlat +lon_0=180 +pm=180")
dest   = Projection("+proj=moll +lon_0=0")
projtup = Node([source, dest])

pm = Node(0)
oxs, oys = xygrid(lon, lat)

xs, ys = Node(copy(oxs)), Node(copy(oys))

on(projtup) do projtup
    # restore the arrays to their original values
    xs.val .= oxs
    ys.val .= oys
    source, dest = projtup
    # Transform the underlying array
    Proj4.transform!(source, dest, vec(xs.val), vec(ys.val))
    # Notify the arrays that they've been updated
    xs[] = xs[]
    ys[] = ys[]
end

projtup[] = [source, dest]

scene = surface(xs, ys; color = cf, shading = false, show_axis = false)

ph = on(pm) do pm
    projtup[] = [Projection("+proj=lonlat +lon_0=180 +pm=$pm"), Projection("+proj=moll +lon_0=-pm +pm=$pm")]
    circshift!(cf.val, field[], pm)
    cf[] = cf[]
end

record(scene, "test.mp4", 0:10:180; framerate = 2) do i
    pm[] = i
end

ext: jl


text:
using Conda
Conda.add("vsearch", :vsearch, channel="bioconda")
Conda.add("TaxonKit", :vsearch, channel="bioconda")


text:
@defcomp climatedynamics begin
    # Total radiative forcing
    radforc = Parameter(index=[time])

    # Average global temperature
    temp = Variable(index=[time])

    # lifetempconst
    lifetempconst = Parameter()

    # lifetemplin
    lifetemplin = Parameter()

    # lifetempqd
    lifetempqd = Parameter()

    # Climate sensitivity
    climatesensitivity = Parameter(default = 2.999999803762826)

    function run_timestep(p, v, d, t)

        if is_first(t)
            v.temp[t] = 0.20
        else
            LifeTemp = max(p.lifetempconst + p.lifetemplin * p.climatesensitivity + p.lifetempqd * p.climatesensitivity^2.0, 1.0)

            delaytemp = 1.0 / LifeTemp

            temps = p.climatesensitivity / 5.35 / log(2.0)

            # Calculate temperature
            dtemp = delaytemp * temps * p.radforc[t] - delaytemp * v.temp[t - 1]

            v.temp[t] = v.temp[t - 1] + dtemp
        end
    end
end


text:
using AlgebraOfGraphics, Makie, Test
using AlgebraOfGraphics: Sorted
using AlgebraOfGraphics: separate
using AlgebraOfGraphics: Arguments, NamedArguments

@testset "utils" begin
    v1 = [1, 2, 7, 11]
    v2 = [1, 3, 4, 5.1]
    @test AlgebraOfGraphics.mergesorted(v1, v2) == [1, 2, 3, 4, 5.1, 7, 11]
    @test AlgebraOfGraphics.mergesorted(v2, v1) == [1, 2, 3, 4, 5.1, 7, 11]
    @test_throws ArgumentError AlgebraOfGraphics.mergesorted([2, 1], [1, 3])
    @test_throws ArgumentError AlgebraOfGraphics.mergesorted([1, 2], [10, 3])

    e1 = (-3, 11)
    e2 = (-5, 10)
    @test AlgebraOfGraphics.extend_extrema(e1, e2) == (-5, 11)

    @test AlgebraOfGraphics.midpoints(1:10) == 1.5:9.5

    s = Any[1, 2]
    t = NamedArguments((a=3, b=4))
    u = NamedArguments((b=5, c=6))
    v = Any[]
    w = AlgebraOfGraphics.concatenate_values(s, t, u)
    @test w == Any[1, 2, 3, 4, 5, 6]
end

@testset "arguments" begin
    s = Arguments([10, 20, 30])
    @test s == Any[10, 20, 30]
    t = map(keys(s), s) do k, v
        return "key $k and value $v"
    end
    @test t == ["key 1 and value 10", "key 2 and value 20", "key 3 and value 30"]

    s = NamedArguments([:a, :b, :c], [10, 20, 30])
    @test s == NamedArguments((a=10, b=20, c=30))
    t = map(keys(s), s) do k, v
        return "key $k and value $v"
    end
    @test t == NamedArguments(
        [:a, :b, :c],
        ["key a and value 10", "key b and value 20", "key c and value 30"]
    )

    s = NamedArguments([:a, :b, :c], [1, 2, 3])
    odd, even = separate(isodd, s)
    @test odd == NamedArguments((a=1, c=3))
    @test even == NamedArguments((b=2,))

    t = AlgebraOfGraphics.set(s, :d => 5, :a => 3)
    @test t == NamedArguments([:a, :b, :c, :d], [3, 2, 3, 5])

    s = NamedArguments([:a, :b, :c], [1, 2, 3])
    t = AlgebraOfGraphics.set(s)
    @test s == t
    @test keys(s) !== keys(t)

    s = NamedArguments([:a, :b, :c], [1, 2, 3])
    t = AlgebraOfGraphics.filterkeys(!in((:a, :b)), s)
    @test t == NamedArguments([:c], [3])
end

@testset "layers" begin
    df = (x = rand(1000), y = rand(1000), c = rand(["a", "b", "c"], 1000))
    d = mapping(:x, :y, color=:c)
    s = visual(color=:red) + mapping(markersize=:c)
    layers = data(df) * d * s
    @test layers[1].transformation isa AlgebraOfGraphics.Visual
    @test layers[1].transformation.attributes[:color] == :red
    @test layers[1].positional == Any[:x, :y]
    @test layers[1].named == NamedArguments((color=:c,))
    @test layers[1].data == df
end

@testset "process_mappings" begin
    df = (x=rand(1000), y=rand(1000), z=rand(1000), c=rand(["a", "b", "c"], 1000))
    d = mapping(:x => exp, [:y, :z], color=:c, marker = dims(1) => renamer(["a", "b"]))
    layer = data(df) * d
    processedlayer = AlgebraOfGraphics.process_mappings(layer)
    @test processedlayer.positional[1] == fill(map(exp, df.x))
    @test processedlayer.positional[2] == [df.y, df.z]
    @test processedlayer.primary[:color] == fill(df.c)
    @test processedlayer.primary[:marker] == [fill(Sorted(1, "a")), fill(Sorted(2, "b"))]
    @test processedlayer.named == NamedArguments((;))
    @test processedlayer.labels[1] == fill("x")
    @test processedlayer.labels[2] ==  ["y", "z"]
    @test processedlayer.labels[:color] == fill("c")
    @test processedlayer.labels[:marker] == ""
end

@testset "grouping" begin
    df = (x=rand(1000), y=rand(1000), z=rand(1000), w=rand(1000), c=rand(["a", "b", "c"], 1000))
    df.c[1:3] .= ["a", "b", "c"] # ensure all three values exist
    d = mapping(:x => exp, [:y, :z], color=:c, marker=dims(1) => t -> ["1", "2"][t], markersize=:w)
    layer = data(df) * d
    processedlayer = AlgebraOfGraphics.ProcessedLayer(layer)
    processedlayers = map(CartesianIndices(AlgebraOfGraphics.shape(processedlayer))) do c
        primary, positional, named = map((processedlayer.primary, processedlayer.positional, processedlayer.named)) do tup
            return map(v -> v[c], tup)
        end
        labels = map(l -> AlgebraOfGraphics.getnewindex(l, c), processedlayer.labels)
        return ProcessedLayer(processedlayer; primary, positional, named, labels)
    end
    @test length(processedlayers) == 6
    for i in 1: 6
        @test processedlayers[i].plottype === Any
        @test isempty(processedlayers[i].attributes)
        @test processedlayers[i].labels[1] == "x"
        @test processedlayers[i].labels[2] == (i ≤ 3 ? "y" : "z")
        @test processedlayers[i].labels[:color] == "c"
        @test processedlayers[i].labels[:marker] == ""
        @test processedlayers[i].labels[:markersize] == "w"
    end

    @test processedlayers[1].primary[:color] == "a"
    @test processedlayers[1].primary[:marker] == "1"
    @test processedlayers[1].positional[1] == exp.(df.x[df.c .== "a"])
    @test processedlayers[1].positional[2] == df.y[df.c .== "a"]
    @test processedlayers[1].named[:markersize] == df.w[df.c .== "a"]

    @test processedlayers[2].primary[:color] == "b"
    @test processedlayers[2].primary[:marker] == "1"
    @test processedlayers[2].positional[1] == exp.(df.x[df.c .== "b"])
    @test processedlayers[2].positional[2] == df.y[df.c .== "b"]
    @test processedlayers[2].named[:markersize] == df.w[df.c .== "b"]

    @test processedlayers[3].primary[:color] == "c"
    @test processedlayers[3].primary[:marker] == "1"
    @test processedlayers[3].positional[1] == exp.(df.x[df.c .== "c"])
    @test processedlayers[3].positional[2] == df.y[df.c .== "c"]
    @test processedlayers[3].named[:markersize] == df.w[df.c .== "c"]

    @test processedlayers[4].primary[:color] == "a"
    @test processedlayers[4].primary[:marker] == "2"
    @test processedlayers[4].positional[1] == exp.(df.x[df.c .== "a"])
    @test processedlayers[4].positional[2] == df.z[df.c .== "a"]
    @test processedlayers[4].named[:markersize] == df.w[df.c .== "a"]

    @test processedlayers[5].primary[:color] == "b"
    @test processedlayers[5].primary[:marker] == "2"
    @test processedlayers[5].positional[1] == exp.(df.x[df.c .== "b"])
    @test processedlayers[5].positional[2] == df.z[df.c .== "b"]
    @test processedlayers[5].named[:markersize] == df.w[df.c .== "b"]

    @test processedlayers[6].primary[:color] == "c"
    @test processedlayers[6].primary[:marker] == "2"
    @test processedlayers[6].positional[1] == exp.(df.x[df.c .== "c"])
    @test processedlayers[6].positional[2] == df.z[df.c .== "c"]
    @test processedlayers[6].named[:markersize] == df.w[df.c .== "c"]
end

@testset "helpers" begin
    r = renamer("a" => "A", "b" => "B", "c" => "C")
    @test r("a") == Sorted(1, "A")
    @test r("b") == Sorted(2, "B")
    @test r("c") == Sorted(3, "C")
    @test_throws KeyError r("d")
    @test string(r("a")) == "A"
    @test string(r("b")) == "B"
    @test string(r("c")) == "C"
    @test r("a") < r("b") < r("c")
    @test r("a") == r("a")
    @test r("a") != r("b")
    r̂ = renamer(["a" => "A", "b" => "B", "c" => "C"])
    @test r̂("a") == Sorted(1, "A")
    @test r̂("b") == Sorted(2, "B")
    @test r̂("c") == Sorted(3, "C")
    @test_throws KeyError r̂("d")
    @test string(r̂("a")) == "A"
    @test string(r̂("b")) == "B"
    @test string(r̂("c")) == "C"
    @test r̂("a") < r̂("b") < r̂("c")
    @test r̂("a") == r̂("a")
    @test r̂("a") != r̂("b")

    s = sorter("b", "c", "a")
    @test s("a") == Sorted(3, "a")
    @test s("b") == Sorted(1, "b")
    @test s("c") == Sorted(2, "c")
    @test_throws KeyError s("d")
    @test string(s("a")) == "a"
    @test string(s("b")) == "b"
    @test string(s("c")) == "c"
    @test s("b") < s("c") < s("a")
    @test s("a") == s("a")
    @test s("a") != s("b")
    ŝ = sorter(["b", "c", "a"])
    @test ŝ("a") == Sorted(3, "a")
    @test ŝ("b") == Sorted(1, "b")
    @test ŝ("c") == Sorted(2, "c")
    @test_throws KeyError ŝ("d")
    @test string(ŝ("a")) == "a"
    @test string(ŝ("b")) == "b"
    @test string(ŝ("c")) == "c"
    @test ŝ("b") < ŝ("c") < ŝ("a")
    @test ŝ("a") == ŝ("a")
    @test ŝ("a") != ŝ("b")

    a = Sorted(1, [1, 2])
    b = Sorted(1, [1, 2])
    c = Sorted(1, [1, 3])
    @test a == b
    @test hash(a) == hash(b)
    @test a != c
    @test hash(a) != hash(c)

    @test string(nonnumeric(1)) == "1"
    @test isless(nonnumeric(1), nonnumeric(2))
end

@testset "legend_merging" begin
    mwe_data = (; t = 1:20, y1 = ones(20), y2 = rand(20))
    mwe_names = ["f_n", "f_d"]
    plt = data(mwe_data) *
        mapping("t", ["y1", "y2"] .=> "y"; color = dims(1) => (i -> mwe_names[i]), marker=dims(1)) *
        (visual(Lines) + visual(Scatter))
    @test_throws ArgumentError draw(plt)
end


text:
type MVector{S, T} <: StaticVector{T}
    data::NTuple{S, T}

    function MVector(in::NTuple{S, T})
        new(in)
    end

    function MVector(in::NTuple{S})
        new(convert_ntuple(T,in))
    end

    function MVector(in::T)
        new((in,))
    end

    function MVector()
        new()
    end
end

@inline (::Type{MVector}){S}(x::NTuple{S}) = MVector{S}(x)
@inline (::Type{MVector{S}}){S, T}(x::NTuple{S,T}) = MVector{S,T}(x)
@inline (::Type{MVector{S}}){S, T <: Tuple}(x::T) = MVector{S,promote_tuple_eltype(T)}(x)

# Some more advanced constructor-like functions
@inline zeros{N}(::Type{MVector{N}}) = zeros(MVector{N,Float64})
@inline ones{N}(::Type{MVector{N}}) = ones(MVector{N,Float64})

#####################
## MVector methods ##
#####################

@pure size{S}(::Union{MVector{S},Type{MVector{S}}}) = (S, )
@pure size{S,T}(::Type{MVector{S,T}}) = (S,)

@propagate_inbounds function getindex(v::MVector, i::Integer)
    v.data[i]
end

# Mutating setindex!
@propagate_inbounds setindex!{S,T}(v::MVector{S,T}, val, i::Integer) = setindex!(v, convert(T, val), i)
@inline function setindex!{S,T}(v::MVector{S,T}, val::T, i::Integer)
    @boundscheck if i < 1 || i > length(v)
        throw(BoundsError())
    end

    if isbits(T)
        unsafe_store!(Base.unsafe_convert(Ptr{T}, Base.data_pointer_from_objref(v)), val, i)
    else
        # This one is thought to be unsafe (#27)
        #unsafe_store!(Base.unsafe_convert(Ptr{Ptr{Void}}, Base.data_pointer_from_objref(v.data)), Base.data_pointer_from_objref(val), i)
        error("setindex!() with non-isbits eltype is not supported by StaticArrays")
    end

    return val
end

@inline Tuple(v::MVector) = v.data

@inline function Base.unsafe_convert{N,T}(::Type{Ptr{T}}, v::MVector{N,T})
    Base.unsafe_convert(Ptr{T}, Base.data_pointer_from_objref(v))
end

macro MVector(ex)
    if isa(ex, Expr) && ex.head == :vect
        return esc(Expr(:call, MVector{length(ex.args)}, Expr(:tuple, ex.args...)))
    elseif isa(ex, Expr) && ex.head == :ref
        return esc(Expr(:call, Expr(:curly, :MVector, length(ex.args[2:end]), ex.args[1]), Expr(:tuple, ex.args[2:end]...)))
    elseif isa(ex, Expr) && ex.head == :comprehension
        if length(ex.args) != 1 || !isa(ex.args[1], Expr) || ex.args[1].head != :generator
            error("Expected generator in comprehension, e.g. [f(i) for i = 1:3]")
        end
        ex = ex.args[1]
        if length(ex.args) != 2
            error("Use a one-dimensional comprehension for @MVector")
        end

        rng = eval(current_module(), ex.args[2].args[2])
        f = gensym()
        f_expr = :($f = ($(ex.args[2].args[1]) -> $(ex.args[1])))
        exprs = [:($f($j)) for j in rng]

        return quote
            $(Expr(:meta, :inline))
            $(esc(f_expr))
            $(esc(Expr(:call, Expr(:curly, :MVector, length(rng)), Expr(:tuple, exprs...))))
        end
    elseif isa(ex, Expr) && ex.head == :typed_comprehension
        if length(ex.args) != 2 || !isa(ex.args[2], Expr) !! ex.args[2].head != :generator
            error("Expected generator in typed comprehension, e.g. Float64[f(i) for i = 1:3]")
        end
        T = ex.args[1]
        ex = ex.args[2]
        if length(ex.args) != 2
            error("Use a one-dimensional comprehension for @MVector")
        end

        rng = eval(current_module(), ex.args[2].args[2])
        f = gensym()
        f_expr = :($f = ($(ex.args[2].args[1]) -> $(ex.args[1])))
        exprs = [:($f($j)) for j in rng]

        return quote
            $(Expr(:meta, :inline))
            $(esc(f_expr))
            $(esc(Expr(:call, Expr(:curly, :MVector, length(rng), T), Expr(:tuple, exprs...))))
        end
    elseif isa(ex, Expr) && ex.head == :call
        if ex.args[1] == :zeros || ex.args[1] == :ones || ex.args[1] == :rand ||ex.args[1] == :randn
            if length(ex.args) == 2
                return quote
                    $(Expr(:meta, :inline))
                    $(esc(ex.args[1]))(MVector{$(esc(ex.args[2]))})
                end
            elseif length(ex.args) == 3
                return quote
                    $(Expr(:meta, :inline))
                    $(esc(ex.args[1]))(MVector{$(esc(ex.args[3])), $(esc(ex.args[2]))})
                end
            else
                error("@MVector expected a 1-dimensional array expression")
            end
        elseif ex.args[1] == :fill
            if length(ex.args) == 3
                return quote
                    $(Expr(:meta, :inline))
                    $(esc(ex.args[1]))($(esc(ex.args[2])), MVector{$(esc(ex.args[3]))})
                end
            else
                error("@MVector expected a 1-dimensional array expression")
            end
        else
            error("@MVector only supports the zeros(), ones(), fill(), rand() and randn() functions.")
        end
    else
        error("Use @MVector [a,b,c] or @MVector([a,b,c])")
    end
end


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ 2fdecb39-bff4-4f3c-b3f5-35438541a29d
# import necessary libraries
using Random

# ╔═╡ 29802e3a-8dfe-41c7-a25e-2f0961949854
using Statistics

# ╔═╡ cba2ccb0-c5a2-11eb-11be-41eb0278cb72
md"""
### Logistic Regression

[Reference 1](https://towardsdatascience.com/introduction-to-logistic-regression-66248243c148)\
[Reference 2](https://towardsdatascience.com/building-a-logistic-regression-in-python-301d27367c24)\
[ML Glossary](https://ml-cheatsheet.readthedocs.io/en/latest/logistic_regression.html)

#### Sigmoid Function

$$f(x)=\frac{1}{1+e^{-(x)}}$$

which maps predicted values to probabilities

#### Hypothesis Representation

$$Z=\beta_0+\beta_1X$$

$$h\Theta(x)=\text{sigmoid}(Z)=\frac{1}{1+e^{-(\beta_0+\beta_1X)}}$$

#### Cost Function

$$\text{Cost}(h_\theta(x),y)=
    \begin{cases}
    -\log(h_\theta(x)) & \quad \text{if } y = 1\\
    -\log(1-h_\theta(x)) & \quad \text{if } y = 0
    \end{cases}$$

which is also

$$J(\theta)=-\frac{1}{m}\sum \Big[ y^{(i)}\log(h\theta(x^{(i)})) + (1-y^{(i)})\log(1-h\theta(x^{(i)})) \Big]$$

which should be minimized

#### Gradient Descent

$$\theta_j := \theta_j - \alpha\frac{\partial}{\partial\theta_j}J(\theta)$$

In this case:

$$\theta_j := \theta_j - \alpha\sum^m_{i=1}(h_\theta(x^{(i)})-y^{(i)})x_j^{(i)}$$
"""

# ╔═╡ 3761501d-f02e-41cb-9003-74cb68283d74
begin
X_data = rand(Float32, (100, 5))
Y_data = rand(0:1, 100)
size(X_data), size(Y_data)
end

# ╔═╡ ad8d7761-f843-46f4-a7c6-4a5442d1e17a
function sigmoid(Z::Number)::Float32
    k = 1 + MathConstants.e ^ (-Float32(Z))
    return 1 / k
end

# ╔═╡ d457cb21-4e23-4527-a2e4-9654579f714a
typeof(X_data)

# ╔═╡ fcf0ea6c-72c4-40b9-bca9-d1ad40b67dd2
function scale(X::Array)::Array
    @assert ndims(X) == 2
    u = mean(X, dims=1) # compute mean
    s = std(X, dims=1)  # compute standard deviation
    res = (X .- u) ./ s
    return res
end

# ╔═╡ c7a7f2af-fa56-4bbf-afa1-a64b1af5b78d
X_data_scaled = scale(X_data);

# ╔═╡ 1a17435a-484c-4ea8-809b-1d3b0a548288
mean(X_data_scaled, dims=1)

# ╔═╡ 5a87d92e-643f-4032-a6fb-226ce7cab2af
std(X_data_scaled, dims=1)

# ╔═╡ 0890963a-806a-456a-89cf-ad68b11b41c3
# update sigmoid function
function sigmoid(Z::Array)::Array
    denom = 1 .+ (MathConstants.e .^ (-Z))
    return 1 ./ denom
end

# ╔═╡ bf7c00f8-1a1b-4010-b2a9-c02f33351bd5
# test sigmoid function
sigmoid(1234)

# ╔═╡ 4c1f5e23-a985-4f8f-81f5-6b780412a4ad
sigmoid(0)

# ╔═╡ 7fa3c5e2-8431-49d7-94c2-d3885a6a79f6
sigmoid(0.123)

# ╔═╡ ed47d544-9ccd-47ff-bf2a-01a9b2b365b0
# test it
sigmoid(X_data_scaled)

# ╔═╡ 8eaa7843-aa2b-4af3-9e32-629629d0fc36
begin
beta = randn(size(X_data)[2]+1)
size(beta)
end

# ╔═╡ e4d93713-c75d-4cbf-a31c-a9194f25280d
function cost(X::Array, y::Array, beta::Array)::AbstractFloat
    @assert ndims(X) == 2
    @assert ndims(y) == 1
    @assert ndims(beta) == 1
    @assert size(X) == (size(y)[1], size(beta)[1]-1)
    m = size(X)[1]
    X_extended = hcat(X, ones(size(X)[1]))
    X_combined = X_extended * reshape(beta, (size(beta)[1], 1))
    prob = sigmoid(X_combined)
    y_prep = reshape(y, (size(y)[1], 1))
    vec = y_prep .* (log.(prob)) .+ (1 .- y_prep) .* (log.(1 .- prob))
    cost = -1 / float(m) * sum(vec)
    return cost
end

# ╔═╡ 4ae75dea-6fc1-4d10-8c63-176a01193caf
# try it
cost(X_data_scaled, Y_data, beta)

# ╔═╡ 83a969a4-2704-421e-ba02-70875e5b22fc
# prediction function
function predict_proba(X::Array, beta::Array)::Array
    @assert ndims(X) == 2
    @assert ndims(beta) == 1
    @assert size(X)[2] == size(beta)[1]-1
    X_extended = hcat(X, ones(size(X)[1]))
    X_combined = X_extended * reshape(beta, (size(beta)[1], 1))
    prob = sigmoid(X_combined)
    prob = [(prob...)...] # flatten to 1d array
    return prob
end

# ╔═╡ a04e19f3-8c20-47b4-b30a-ccc0696bf54b
predict_proba(X_data_scaled, beta)

# ╔═╡ 7d94d627-941c-41fb-8fce-cc0da6972f74
function predict(X::Array, beta::Array)::Array
    @assert ndims(X) == 2
    @assert ndims(beta) == 1
    @assert size(X)[2] == size(beta)[1]-1
    X_extended = hcat(X, ones(size(X)[1]))
    X_combined = X_extended * reshape(beta, (size(beta)[1], 1))
    prob = sigmoid(X_combined)
    prob = [(prob...)...] # flatten to 1d array
    prob = map(m -> m >= 0.5 ? 1 : 0, prob)
    return prob
end

# ╔═╡ 6418946f-eff9-43fd-8654-415fc41579b5
predict(X_data_scaled, beta)

# ╔═╡ 4c6ee096-98f5-48b8-9977-f2485365759c
# inplace learning function (Gradient Descent)
function learn!(X::Array, y::Array, beta::Array, alpha::AbstractFloat)
    @assert ndims(X) == 2
    @assert ndims(y) == 1
    @assert ndims(beta) == 1
    @assert size(X) == (size(y)[1], size(beta)[1]-1)
    predictions = predict_proba(X, beta)
    offset = predictions .- y
    offset = reshape(offset, (size(offset)[1], 1))
    X_extended = hcat(X, ones(size(X)[1]))
    gradients = X_extended' * offset
    gradients = gradients ./ size(X)[1]
    gradients = gradients .* alpha
    beta .= beta .- [(gradients...)...]
    return nothing
end

# ╔═╡ addad88e-317d-4086-a5e4-a9682f80f54a
beta

# ╔═╡ 2c3b2fed-e9aa-44b5-a12d-40281ae487d9
learn!(X_data_scaled, Y_data, beta, 0.01)

# ╔═╡ 63ab48c3-8c92-4935-a4ce-36a7049e5c86
beta

# ╔═╡ c0f88234-87ec-41ca-8a90-df7d89997090
# define the logistic regression function
function train(X::Array, y::Array; learning_rate::AbstractFloat=0.01, max_iter::Integer=10, return_all::Bool=false)::Array
    @assert ndims(X) == 2
    @assert ndims(y) == 1
    @assert size(X)[1] == size(y)[1]
    @assert max_iter >= 0
    beta = Random.randn(size(X)[2]+1)
    res = nothing
    if return_all
        res = reshape(beta, (1, size(beta)[1]))
    else
        res = beta
    end
    for i = 1:max_iter
        if return_all
            beta = learn(X, y, beta, learning_rate)
            res = cat(res, reshape(beta, (1, size(beta)[1])), dims=1)
        else
            learn!(X, y, beta, learning_rate)
            res .= beta
        end
    end
    return res
end

# ╔═╡ addd77a3-6d60-4387-8bad-65c55fc0171e
# accuracy function
function accuracy(y_pred::Array, y_real::Array)::AbstractFloat
    @assert ndims(y_pred) == ndims(y_real) == 1
    @assert size(y_pred) == size(y_real)
    sum = 0
    for (m, n) in zip(y_pred, y_real)
        if m == n
            sum += 1
        end
    end
    return sum / size(y_pred)[1]
end

# ╔═╡ eb2362a6-8d78-4f38-8b53-a1b9a1f54820
accuracy(predict(X_data_scaled, beta), Y_data)

# ╔═╡ acfefff5-53c3-421c-84c8-441007cbba2e
accuracy(predict(X_data_scaled, train(X_data_scaled, Y_data, max_iter=100, learning_rate=0.5)), Y_data)

# ╔═╡ 8d99210c-eabe-4fe7-a529-4e83be5e6cf8
md"""
Since it is linear model, we should not expect it to have a high accuracy on a completely randomly generated dataset\
But we can see the improvements
"""

# ╔═╡ Cell order:
# ╟─cba2ccb0-c5a2-11eb-11be-41eb0278cb72
# ╠═2fdecb39-bff4-4f3c-b3f5-35438541a29d
# ╠═3761501d-f02e-41cb-9003-74cb68283d74
# ╠═ad8d7761-f843-46f4-a7c6-4a5442d1e17a
# ╠═bf7c00f8-1a1b-4010-b2a9-c02f33351bd5
# ╠═4c1f5e23-a985-4f8f-81f5-6b780412a4ad
# ╠═7fa3c5e2-8431-49d7-94c2-d3885a6a79f6
# ╠═d457cb21-4e23-4527-a2e4-9654579f714a
# ╠═29802e3a-8dfe-41c7-a25e-2f0961949854
# ╠═fcf0ea6c-72c4-40b9-bca9-d1ad40b67dd2
# ╠═c7a7f2af-fa56-4bbf-afa1-a64b1af5b78d
# ╠═1a17435a-484c-4ea8-809b-1d3b0a548288
# ╠═5a87d92e-643f-4032-a6fb-226ce7cab2af
# ╠═0890963a-806a-456a-89cf-ad68b11b41c3
# ╠═ed47d544-9ccd-47ff-bf2a-01a9b2b365b0
# ╠═8eaa7843-aa2b-4af3-9e32-629629d0fc36
# ╠═e4d93713-c75d-4cbf-a31c-a9194f25280d
# ╠═4ae75dea-6fc1-4d10-8c63-176a01193caf
# ╠═83a969a4-2704-421e-ba02-70875e5b22fc
# ╠═a04e19f3-8c20-47b4-b30a-ccc0696bf54b
# ╠═7d94d627-941c-41fb-8fce-cc0da6972f74
# ╠═6418946f-eff9-43fd-8654-415fc41579b5
# ╠═4c6ee096-98f5-48b8-9977-f2485365759c
# ╠═addad88e-317d-4086-a5e4-a9682f80f54a
# ╠═2c3b2fed-e9aa-44b5-a12d-40281ae487d9
# ╠═63ab48c3-8c92-4935-a4ce-36a7049e5c86
# ╠═c0f88234-87ec-41ca-8a90-df7d89997090
# ╠═addd77a3-6d60-4387-8bad-65c55fc0171e
# ╠═eb2362a6-8d78-4f38-8b53-a1b9a1f54820
# ╠═acfefff5-53c3-421c-84c8-441007cbba2e
# ╟─8d99210c-eabe-4fe7-a529-4e83be5e6cf8

max_stars_repo_name: teamclouday/JuliaLearn

max_issues_repo_name: teamclouday/JuliaLearn

max_forks_repo_name: teamclouday/JuliaLearn


text:
"""
    to_vec(x)

Transform `x` into a `Vector`, and return the vector, and a closure which inverts the
transformation.
"""
function to_vec(x::Real)
    function Real_from_vec(x_vec)
        return first(x_vec)
    end
    return [x], Real_from_vec
end

function to_vec(z::Complex)
    function Complex_from_vec(z_vec)
        return Complex(z_vec[1], z_vec[2])
    end
    return [real(z), imag(z)], Complex_from_vec
end

# Base case -- if x is already a Vector{<:Real} there's no conversion necessary.
to_vec(x::Vector{<:Real}) = (x, identity)

# get around the constructors and make the type directly
# Note this is moderately evil accessing julia's internals
if VERSION >= v"1.3"
    @generated function _force_construct(T, args...)
        return Expr(:splatnew, :T, :args)
    end
else
    @generated function _force_construct(T, args...)
        return Expr(:new, :T, Any[:(args[$i]) for i in 1:length(args)]...)
    end
end

# Fallback method for `to_vec`. Won't always do what you wanted, but should be fine a decent
# chunk of the time.
function to_vec(x::T) where {T}
    Base.isstructtype(T) || throw(error("Expected a struct type"))
    isempty(fieldnames(T)) && return (Bool[], _ -> x) # Singleton types

    val_vecs_and_backs = map(name -> to_vec(getfield(x, name)), fieldnames(T))
    vals = first.(val_vecs_and_backs)
    backs = last.(val_vecs_and_backs)

    v, vals_from_vec = to_vec(vals)
    function structtype_from_vec(v::Vector{<:Real})
        val_vecs = vals_from_vec(v)
        values = map((b, v) -> b(v), backs, val_vecs)
        try
            T(values...)
        catch MethodError
            return _force_construct(T, values...)
        end
    end
    return v, structtype_from_vec
end

function to_vec(x::DenseVector)
    x_vecs_and_backs = map(to_vec, x)
    x_vecs, backs = first.(x_vecs_and_backs), last.(x_vecs_and_backs)
    function Vector_from_vec(x_vec)
        sz = cumsum(map(length, x_vecs))
        x_Vec = [backs[n](x_vec[sz[n] - length(x_vecs[n]) + 1:sz[n]]) for n in eachindex(x)]
        return oftype(x, x_Vec)
    end
    # handle empty x
    x_vec = isempty(x_vecs) ? eltype(eltype(x_vecs))[] : reduce(vcat, x_vecs)
    return x_vec, Vector_from_vec
end

function to_vec(x::DenseArray)
    x_vec, from_vec = to_vec(vec(x))

    function Array_from_vec(x_vec)
        return oftype(x, reshape(from_vec(x_vec), size(x)))
    end

    return x_vec, Array_from_vec
end

# Some specific subtypes of AbstractArray.
function to_vec(x::Base.ReshapedArray{<:Any, 1})
    x_vec, from_vec = to_vec(parent(x))
    function ReshapedArray_from_vec(x_vec)
        p = from_vec(x_vec)
        return Base.ReshapedArray(p, x.dims, x.mi)
    end

    return x_vec, ReshapedArray_from_vec
end

# To return a SubArray we would endup needing to copy the `parent` of `x` in `from_vec`
# which doesn't seem particularly useful. So we just convert the view into a copy.
# we might be able to do something more performant but this seems good for now.
to_vec(x::Base.SubArray) = to_vec(copy(x))

function to_vec(x::T) where {T<:LinearAlgebra.AbstractTriangular}
    x_vec, back = to_vec(Matrix(x))
    function AbstractTriangular_from_vec(x_vec)
        return T(reshape(back(x_vec), size(x)))
    end
    return x_vec, AbstractTriangular_from_vec
end

function to_vec(x::T) where {T<:LinearAlgebra.HermOrSym}
    x_vec, back = to_vec(Matrix(x))
    function HermOrSym_from_vec(x_vec)
        return T(back(x_vec), x.uplo)
    end
    return x_vec, HermOrSym_from_vec
end

function to_vec(X::Diagonal)
    x_vec, back = to_vec(Matrix(X))
    function Diagonal_from_vec(x_vec)
        return Diagonal(back(x_vec))
    end
    return x_vec, Diagonal_from_vec
end

function to_vec(X::Transpose)
    x_vec, back = to_vec(Matrix(X))
    function Transpose_from_vec(x_vec)
        return Transpose(permutedims(back(x_vec)))
    end
    return x_vec, Transpose_from_vec
end

function to_vec(x::Transpose{<:Any, <:AbstractVector})
    x_vec, back = to_vec(Matrix(x))
    Transpose_from_vec(x_vec) = Transpose(vec(back(x_vec)))
    return x_vec, Transpose_from_vec
end

function to_vec(X::Adjoint)
    x_vec, back = to_vec(Matrix(X))
    function Adjoint_from_vec(x_vec)
        return Adjoint(conj!(permutedims(back(x_vec))))
    end
    return x_vec, Adjoint_from_vec
end

function to_vec(x::Adjoint{<:Any, <:AbstractVector})
    x_vec, back = to_vec(Matrix(x))
    Adjoint_from_vec(x_vec) = Adjoint(conj!(vec(back(x_vec))))
    return x_vec, Adjoint_from_vec
end

function to_vec(X::T) where {T<:PermutedDimsArray}
    x_vec, back = to_vec(parent(X))
    function PermutedDimsArray_from_vec(x_vec)
        X_parent = back(x_vec)
        return T(X_parent)
    end
    return x_vec, PermutedDimsArray_from_vec
end

# Factorizations

function to_vec(x::F) where {F <: SVD}
    # Convert the vector S to a matrix so we can work with a vector of matrices
    # only and inferrence work
    v = [x.U, reshape(x.S, length(x.S), 1), x.Vt]
    x_vec, back = to_vec(v)
    function SVD_from_vec(v)
        U, Smat, Vt = back(v)
        return F(U, vec(Smat), Vt)
    end
    return x_vec, SVD_from_vec
end

function to_vec(x::Cholesky)
    x_vec, back = to_vec(x.factors)
    function Cholesky_from_vec(v)
        return Cholesky(back(v), x.uplo, x.info)
    end
    return x_vec, Cholesky_from_vec
end

function to_vec(x::S) where {U, S <: Union{LinearAlgebra.QRCompactWYQ{U}, LinearAlgebra.QRCompactWY{U}}}
    # x.T is composed of upper triangular blocks. The subdiagonals elements
    # of the blocks are abitrary. We make sure to set all of them to zero
    # to avoid NaN.
    blocksize, cols = size(x.T)
    T = zeros(U, blocksize, cols)

    for i in 0:div(cols - 1, blocksize)
        used_cols = i * blocksize
        n = min(blocksize, cols - used_cols)
        T[1:n, (1:n) .+ used_cols] = UpperTriangular(view(x.T, 1:n, (1:n) .+ used_cols))
    end

    x_vec, back = to_vec([x.factors, T])

    function QRCompact_from_vec(v)
        factors, Tback = back(v)
        return S(factors, Tback)
    end

    return x_vec, QRCompact_from_vec
end

# Non-array data structures

function to_vec(x::Tuple)
    x_vecs_and_backs = map(to_vec, x)
    x_vecs, x_backs = first.(x_vecs_and_backs), last.(x_vecs_and_backs)
    lengths = map(length, x_vecs)
    sz = typeof(lengths)(cumsum(collect(lengths)))
    function Tuple_from_vec(v)
        map(x_backs, lengths, sz) do x_back, l, s
            return x_back(v[s - l + 1:s])
        end
    end
    return reduce(vcat, x_vecs), Tuple_from_vec
end

function to_vec(x::NamedTuple)
    x_vec, back = to_vec(values(x))
    function NamedTuple_from_vec(v)
        v_vec_vec = back(v)
        return typeof(x)(v_vec_vec)
    end
    return x_vec, NamedTuple_from_vec
end

# Convert to a vector-of-vectors to make use of existing functionality.
function to_vec(d::Dict)
    d_vec, back = to_vec(collect(values(d)))
    function Dict_from_vec(v)
        v_vec_vec = back(v)
        return Dict(key => v_vec_vec[n] for (n, key) in enumerate(keys(d)))
    end
    return d_vec, Dict_from_vec
end

# ChainRulesCore Differentials
function FiniteDifferences.to_vec(x::Tangent{P}) where{P}
    x_canon = canonicalize(x)  # to be safe, fill in every field and put in primal order.
    x_inner = ChainRulesCore.backing(x_canon)
    x_vec, back_inner = FiniteDifferences.to_vec(x_inner)
    function Tangent_from_vec(y_vec)
        y_back = back_inner(y_vec)
        return Tangent{P, typeof(y_back)}(y_back)
    end
    return x_vec, Tangent_from_vec
end

function FiniteDifferences.to_vec(x::AbstractZero)
    function AbstractZero_from_vec(x_vec::Vector)
        return x
    end
    return Bool[], AbstractZero_from_vec
end

function FiniteDifferences.to_vec(t::Thunk)
    v, back = to_vec(unthunk(t))
    Thunk_from_vec = v -> @thunk(back(v))
    return v, Thunk_from_vec
end


text:
function quad_transform(q::QDHT{p,n}, f; rmin = 0, rmax = q.R) where {p,n}
    return hquadrature(rmin, rmax) do r
        return @. r^n * f(r) * Hankel.sphbesselj(p, n, q.k * q.r) /
                  Hankel.sphbesselj_scale(n)
    end[1]
end

function quad_integrate(q::QDHT{p,n}, f; rmin = 0, rmax = q.R) where {p,n}
    return hquadrature(rmin, rmax) do r
        return @. r^n * f(r)
    end[1]
end

function test_transform(
    q::QDHT,
    f;
    fk = nothing,
    quad = true,
    nroundtrip = 1000,
    inplace = true,
    kwargs...,
)
    v = f.(q.r)
    vk = q * v
    @test eltype(vk) === eltype(v)

    inplace && @testset "inplace" begin
        vk2 = similar(vk)
        mul!(vk2, q, v)
        @test all(isapprox.(vk2, vk; kwargs...))

        v2 = similar(v)
        ldiv!(v2, q, vk)
        @test all(isapprox.(v2, v; kwargs...))
    end

    @testset "roundtrip" begin
        vv = copy(v)
        tmp = q * vv
        for _ in 1:nroundtrip
            vv = q \ tmp
            tmp = q * vv
        end
        @test all(isapprox.(v, vv; kwargs...))
    end

    fk !== nothing && @testset "analytical" begin
        vka = fk.(q.k)
        @test all(isapprox.(vk, vka; kwargs...))
    end

    quad && @testset "quadrature" begin
        vki = quad_transform(q, f)
        @test all(isapprox.(vk, vki; kwargs...))
    end
    return nothing
end

function test_l2norm(q::QDHT, f; Efr = nothing, quad = true, kwargs...)
    v = f.(q.r)
    vk = q * v
    Er = integrateR(v .^ 2, q)
    Ek = integrateK(vk .^ 2, q)
    @test all(isapprox.(Er, Ek; kwargs...))

    Efr !== nothing && @testset "analytical" begin
        Era = Efr.(q.r)
        @test all(isapprox.(Er, Era; kwargs...))
    end

    quad && @testset "quadrature" begin
        Eri = quad_integrate(q, r -> f(r)^2)
        @test all(isapprox.(Er, Eri; kwargs...))
    end
    return nothing
end

function test_integrate(q::QDHT, f; Ifr = nothing, quad = true, kwargs...)
    v = f.(q.r)
    Ir = integrateR(v, q)

    Ifr !== nothing && @testset "analytical" begin
        Ira = Ifr.(q.r)
        @test all(isapprox.(Ir, Ira; kwargs...))
    end

    quad && @testset "quadrature" begin
        Iri = quad_integrate(q, f)
        @test all(isapprox.(Ir, Iri; kwargs...))
    end
    return nothing
end

# "dynamical error" from Guizar-Sicairos and Gutierrez-Vega
dynε(ext, est) = 20 * log10.(abs.(ext .- est) ./ maximum(abs.(est)))

logbesseli(ν, x) = x > one(x) ? log(besselix(ν, x)) + x : log(besseli(ν, x))

# log(I_ν(x) / x^ν), approaches -log(2^ν * γ(ν + 1)) as x → 0
logbesselirat(ν, x) = logbesseli(ν, x) - ν * log(x)

# pdf of non-central chi distribution reparameterized with location μ and scale σ
# density wrt rⁿ dr
function pdf_ncchi(r, μ, σ, n)
    ν = n + 1
    x = r / σ
    λ = μ / σ
    logJ = -ν * log(σ) # (rⁿ dr) / (xⁿ dx) = 1 / σⁿ⁺¹
    logpdf = -(x^2 + λ^2) / 2 + logbesselirat(ν / 2 - 1, λ * x) + logJ
    return exp(logpdf)
end

# transform of pdf_ncchi. Adapted from Gradshteyn and Ryzhik 6.633.2.
function trans_pdf_ncchi(k, μ, σ, n, p)
    cninv = 1 / Hankel.sphbesselj_scale(n)
    return cninv * exp(-σ^2 / 2 * k^2) * Hankel.sphbesselj(p, n, μ * k)
end

@testset "QDHT" begin
    @testset "constructors" begin
        q = QDHT(1.0, 10)
        @test q isa QDHT{0,1,Float64}
        @test Hankel.order(q) == 0
        @test Hankel.sphericaldim(q) == 1
        q = QDHT{1}(1.0, 10)
        @test q isa QDHT{1,1,Float64}
        @test Hankel.order(q) == 1
        @test Hankel.sphericaldim(q) == 1
        q = QDHT(2, 1.0, 10)
        @test q isa QDHT{2,1,Float64}
        @test Hankel.order(q) == 2
        @test Hankel.sphericaldim(q) == 1
        q = QDHT{1, 2}(1.0, 10)
        @test q isa QDHT{1,2,Float64}
        @test Hankel.order(q) == 1
        @test Hankel.sphericaldim(q) == 2
        q = QDHT(1, 2, 1.0, 10)
        @test q isa QDHT{1,2,Float64}
        @test Hankel.order(q) == 1
        @test Hankel.sphericaldim(q) == 2
    end

    @testset "Base.inv(::QDHT)" begin
        q = QDHT{1, 2}(1.0, 10)
        qinv = inv(q)
        q2 = QDHT{1, 2}(q.K, q.N)
        @test typeof(q2) === typeof(qinv)
        @test all(map(s -> isapprox(getproperty(q2, s), getproperty(qinv, s); atol=1e-12), propertynames(q2)))
        q2inv = inv(q2)
        @test all(map(s -> isapprox(getproperty(q, s), getproperty(q2inv, s); atol=1e-12), propertynames(q)))
    end

    @testset "Base.reshape(::QDHT, ::Int)" begin
        q = QDHT{1, 2}(1.0, 10)
        q2 = reshape(q, 20)
        @test typeof(q) === typeof(q2)
        @test q2.R == q.R
        @test q2.N == 20
    end

    @testset "cylindrical" begin
        @testset "transform" begin
            R = 4e-2
            N = 256
            w0 = 1e-3
            a = 2/w0
            q = Hankel.QDHT(R, N)
            @test all(isreal.(q.T))
            f(r) = exp(-1//2 * a^2 * r^2)
            fk(k) = 1/a^2 * exp(-k^2/(2*a^2))
            v = f.(q.r)
            global vv = copy(v)
            for _ = 1:1000
                tmp = q * vv
                global vv = q \ tmp
            end
            @test all(isapprox.(v, vv, atol=8e-12))
            vk = q * v
            vka = fk.(q.k)
            fki(k) = hquadrature(r -> r.*f(r).*besselj(0, k.*r), 0, R)[1]
            @test all(isapprox.(vka, vk, atol=9e-22))
            @test fki(q.k[1]) ≈ vk[1] # doing all of them takes too long
            @test fki(q.k[128]) ≈ vk[128]
            Er = Hankel.integrateR(v.^2, q)
            Ek = Hankel.integrateK(vk.^2, q)
            @test Er ≈ Ek
            Er_c = hquadrature(r -> r.*f(r).^2, 0, 1)
            @test Er_c[1] ≈ Er
            # Test integral of function without squaring
            ir = Hankel.integrateR(v, q)
            ir_c = hquadrature(r -> r*f(r), 0, 1)
            @test ir_c[1] ≈ ir
            # Test that in-place transform works
            vk2 = similar(vk)
            vk3 = copy(v)
            mul!(vk2, q, v)
            @test all(vk2 .≈ vk)

            v2d = repeat(v, outer=(1, 16))'
            q2d = Hankel.QDHT(R, N, dim=2)
            v2dk = q2d * v2d
            @test all([all(v2dk[ii, :] ≈ vk) for ii = 1:size(v2dk, 1)])

            f0 = f(0)
            f0q = Hankel.onaxis(vk, q)
            @test f0 ≈ f0q

            f2(r) = sinc(100*r)^2
            v = f2.(q.r);
            vk = q * v
            f0 = f2(0)
            f0q = Hankel.onaxis(vk, q)
            @test f0 ≈ f0q

            # test also with a different function which isn't positive definite
            w0 = 5e-3
            a = 2/w0
            f(r) = exp(-1//2 * a^2 * r^2) * cos(16*a*r)
            v = f.(q.r);
            vk = q * v;
            Ir = Hankel.integrateR(v.^2, q)
            Ik = Hankel.integrateK(vk.^2, q)
            Ir_c = hquadrature(r -> r*f(r)^2, 0, 1)
            @test Ir_c[1] ≈ Ir ≈ Ik
            # Test integral of function without squaring
            ir = Hankel.integrateR(v, q)
            ir_c = hquadrature(r -> r*f(r), 0, 1)
            @test ir_c[1] ≈ ir

            q = QDHT(10, 128); A = exp.(-q.r.^2)
            As = symmetric(A, q)
            @test length(As) == 2*length(A)+1
            @test As[1:128] == A[128:-1:1]
            @test As[129] ≈ 1
            @test As[130:end] == A

            AA = hcat(A, A)
            AAs = symmetric(AA, q; dim=1)
            @test size(AAs, 1) == 2*length(A)+1
            for i=1:2
                @test AAs[1:128, i] == A[128:-1:1]
                @test AAs[129, i] ≈ 1
                @test AAs[130:end, i] == A
            end
        end

        @testset "non-zero orders" begin
            #= Test case from Guizar-Sicairos and Gutierrez-Vega,
              "Computation of quasi-discrete Hankel transforms of integer order for propagating
              optical wave fields" =#
            γ = 5
            f(r) = sinc(2γ*r)
            function fk(ν, p)
                if ν < γ
                    return ν^p * cos(p*π/2) / (2π*γ * sqrt(γ^2 - ν^2)*(γ + sqrt(γ^2 - ν^2))^p)
                else
                    return sin(p*asin(γ/ν))/(2π*γ*sqrt(ν^2 - γ^2))
                end
            end

            R = 3
            N = 256
            q1 = QDHT(1, R, N)
            v = f.(q1.r);
            vk = q1*v;
            vka = fk.(q1.k/2π, 1)/2π
            @test maximum(dynε(vka, vk)) < -10
            @test_throws MethodError onaxis(vk, q1)
            @test_throws DomainError symmetric(v, q1)
            @test integrateR(abs2.(v), q1) ≈ integrateK(abs2.(vk), q1)
            @test isapprox(integrateR(abs2.(v), q1), hquadrature(r -> r*f(r)^2, 0, 3)[1], rtol=1e-2)

            q4 = QDHT(4, R, N)
            v = f.(q4.r);
            vk = q4*v;
            vka = fk.(q4.k/2π, 4)/2π
            @test maximum(dynε(vka, vk)) < -10
            @test integrateR(abs2.(v), q4) ≈ integrateK(abs2.(vk), q4)

            R = 4e-2
            N = 256
            w0 = 1e-2
            a = 2/w0
            q4 = QDHT(4, R, N)
            f(r) = exp(-1//2 * a^2 * r^2)
            fk(k) = 1/a^2 * exp(-k^2/(2*a^2))
            v = f.(q4.r)
            vk = q4 * v
            @test integrateR(abs2.(v), q4) ≈ integrateK(abs2.(vk), q4)
            # note the high tolerance here is due to the lack of samples close to the axis
            @test isapprox(integrateR(abs2.(v), q4), hquadrature(r -> r*f(r)^2, 0, R)[1], rtol=1e-2)

            R = 4e-2
            N = 256
            w0 = 1e-2
            a = 2/w0
            q12 = QDHT(1/2, R, N)
            f(r) = exp(-1//2 * a^2 * r^2)
            fk(k) = 1/a^2 * exp(-k^2/(2*a^2))
            v = f.(q12.r)
            vk = q12 * v
            @test integrateR(abs2.(v), q12) ≈ integrateK(abs2.(vk), q12)
            @test isapprox(integrateR(abs2.(v), q12), hquadrature(r -> r*f(r)^2, 0, R)[1], rtol=1e-3)
        end

        @testset "big floats" begin
            R = 4e-2
            N = 256
            w0 = 4e-3
            a = 2/w0
            q = QDHT(BigFloat(R), 128)
            f(r) = exp(-1//2 * a^2 * r^2)
            fk(k) = 1/a^2 * exp(-k^2/(2*a^2))
            v = f.(q.r)
            global vv = copy(v)
            for _ = 1:100
                tmp = q * vv
                global vv = q \ tmp
            end
            @test all(isapprox.(v, vv, atol=2e-12))
            vk = q * v
            vka = fk.(q.k)
            @test all(isapprox.(vka, vk, atol=5e-22))
        end

        @testset "Gaussian divergence" begin
            q = Hankel.QDHT(12.7e-3, 512)
            λ = 800e-9
            k = 2π/λ
            kz = @. sqrt(k^2 - q.k^2)
            z = 2 # propagation distance
            prop = @. exp(1im*kz*z)
            w0 = 200e-6 # start at focus
            w1 = w0*sqrt(1+(z*λ/(π*w0^2))^2)
            Ir0 = exp.(-2*(q.r/w0).^2)
            Ir1 = exp.(-2*(q.r/w1).^2)*(w0/w1)^2 # analytical solution (in paraxial approx)
            Er0 = sqrt.(Ir0)
            Ek0 = q * Er0
            Ek1 = prop .* Ek0
            Er1 = q \ Ek1
            @test isapprox(abs2.(Er1), Ir1, rtol=1e-6)
            energy = π/2*w0^2 # analytical energy for gaussian beam
            @test energy ≈ 2π*Hankel.integrateR(Ir0, q)
            @test energy ≈ 2π*Hankel.integrateR(abs2.(Er1), q)
            @test energy ≈ 2π*Hankel.integrateK(abs2.(Ek1), q)
        end

        @testset "Oversampling" begin
            R = 4e-2
            N = 256
            w0 = 1e-2
            a = 2/w0
            q = Hankel.QDHT(R, N)
            f(r) = exp(-1//2 * a^2 * r^2)
            v = f.(q.r)
            vo, qo = Hankel.oversample(v, q, factor=4)
            @test qo.N == 4*q.N
            @test all(isapprox(vo, f.(qo.r), rtol=1e-13))
            @test integrateR(abs2.(v), q) ≈ integrateR(abs2.(vo), qo)
            @test f(0) ≈ onaxis(qo*vo, qo)

            qo2 = reshape(q, N * 4)
            vo2 = Hankel.oversample(v, q, qo2)
            @test typeof(qo2) == typeof(qo)
            @test qo2.R == qo.R
            @test qo2.N == qo.N
            @test vo2 == vo

            v2 = hcat(v, v)
            vo, qo = Hankel.oversample(v2, q, factor=4)
            @test qo.N == 4*q.N
            @test all(isapprox(vo[:, 1], f.(qo.r), rtol=1e-13))
            @test all(isapprox(vo[:, 2], f.(qo.r), rtol=1e-13))
            @test integrateR(abs2.(v2), q) ≈ integrateR(abs2.(vo), qo)

            qo2 = reshape(q, N * 4)
            vo2 = Hankel.oversample(v2, q, qo2)
            @test typeof(qo2) == typeof(qo)
            @test qo2.R == qo.R
            @test qo2.N == qo.N
            @test vo2 == vo
        end
    end

    @testset "spherical" begin
        @testset "transform" begin
            @testset "f(r) = exp(-a²r²/2)" begin
                R = 4e-2
                N = 256
                q = Hankel.QDHT{0,2}(R, N)
                w0 = 1e-3
                a = 2 / w0
                @test all(isreal.(q.T))
                f(r) = exp(-1 // 2 * a^2 * r^2)
                fk(k) = 1 / a^3 * exp(-k^2 / (2 * a^2))
                test_transform(q, f; fk = fk, atol = 1e-9)
                test_l2norm(q, f; atol = 1e-20)
                test_integrate(q, f; atol = 1e-20)
                f0 = f(0)
                f0q = Hankel.onaxis(q * f.(q.r), q)
                @test f0 ≈ f0q

                v = f.(q.r)
                vk = q * v
                @testset "2d" begin
                    v2d = repeat(v, outer = (1, 16))'
                    q2d = Hankel.QDHT{0,2}(R, N, dim = 2)
                    v2dk = q2d * v2d
                    @test all([all(v2dk[ii, :] ≈ vk) for ii in 1:size(v2dk, 1)])
                end
            end

            @testset "f(r) = sinc(100r)²" begin
                R = 4e-2
                N = 256
                q = Hankel.QDHT{0,2}(R, N)
                f(r) = sinc(100 * r)^2
                v = f.(q.r)
                vk = q * v
                f0 = f(0)
                f0q = Hankel.onaxis(vk, q)
                @test isapprox(f0, f0q; atol = 1e-6)
            end

            @testset "f(r) = exp(-a²r²/2)cos(16ar)" begin
                R = 1e-2
                N = 256
                q = Hankel.QDHT{0,2}(R, N)
                w0 = 5e-3
                a = 2 / w0
                f(r) = exp(-1 // 2 * a^2 * r^2) * cos(16 * a * r)
                test_transform(q, f; atol = 1e-10, quad = false)
                test_l2norm(q, f; atol = 1e-15)
                test_integrate(q, f; atol = 1e-10)
            end

            @testset "pdf non-central chi" begin
                R = 30
                N = 256
                q = Hankel.QDHT{0,2}(R, N)
                f(r) = pdf_ncchi(r, 10, 2, 2)
                fk(k) = trans_pdf_ncchi(k, 10, 2, 2, 0)
                test_transform(q, f; fk = fk, atol = 2e-12, quad = false)
                test_l2norm(q, f; quad = true, atol = 1e-15)
                test_integrate(q, f; Ifr = _ -> 1, quad = false, atol = 1e-15)
            end
        end

        @testset "symmetric" begin
            q = QDHT{0,2}(10, 128)
            A = exp.(-q.r .^ 2)
            As = symmetric(A, q)
            @test length(As) == 2 * length(A) + 1
            @test As[1:128] == A[128:-1:1]
            @test As[129] ≈ 1
            @test As[130:end] == A

            AA = hcat(A, A)
            AAs = symmetric(AA, q; dim = 1)
            @test size(AAs, 1) == 2 * length(A) + 1
            for i in 1:2
                @test AAs[1:128, i] == A[128:-1:1]
                @test AAs[129, i] ≈ 1
                @test AAs[130:end, i] == A
            end
        end

        @testset "multiple orders/spherical dimensions" begin
            @testset "p=$p, n=$n" for p in (0, 1 / 2, 1, 2, 3), n in (1, 2, 3)
                #= Test case from Guizar-Sicairos and Gutierrez-Vega,
                  "Computation of quasi-discrete Hankel transforms of integer order for propagating
                  optical wave fields" =#
                # Adapted from Gradshteyn and Ryzhik 6.671.1
                γ = 5
                f(r) = r^(-(n + 1) / 2) * sin(2π * γ * r)
                function fk(k)
                    s = p + (n - 1) / 2
                    return if k < 2π * γ
                        return k^s * cos(s * π / 2) / sqrt(4 * π^2 * γ^2 - k^2) /
                               (2π * γ + sqrt(4 * π^2 * γ^2 - k^2))^s / k^((n - 1) / 2)
                    else
                        return sin(s * asin(2π * γ / k)) / sqrt(k^2 - 4 * π^2 * γ^2) /
                               k^((n - 1) / 2)
                    end
                end

                R = 1
                N = 256
                q1 = QDHT{p,n}(R, N)
                v = f.(q1.r)
                vk = q1 * v
                vka = fk.(q1.k)
                err = dynε(vka, vk)
                i = findlast(ki -> ki < 2π * γ, q1.k) # find point flanking discontinuity
                @test maximum(dynε(vka, vk)[[1:(i - 1); (i + 2):end]]) < -10
                if p != 0
                    @test_throws MethodError onaxis(vk, q1)
                end
                @test integrateR(abs2.(v), q1) ≈ integrateK(abs2.(vk), q1)
                # curve is too steep for quadrature, so compare with integral computed with
                # Mathematica
                @test isapprox(integrateR(abs2.(v), q1), 2.358965372, rtol = 1e-2)

                R = 4e-2
                N = 256
                w0 = 1e-2
                a = 2 / w0
                q4 = QDHT{p,n}(R, N)
                g(r) = exp(-1 // 2 * a^2 * r^2)
                v = g.(q4.r)
                vk = q4 * v
                @test integrateR(abs2.(v), q4) ≈ integrateK(abs2.(vk), q4)
                # note the high tolerance here is due to the lack of samples close to the axis
                @test isapprox(
                    integrateR(abs2.(v), q4),
                    quad_integrate(q4, abs2 ∘ g),
                    rtol = 1e-2,
                )

                @testset "pdf non-central chi" begin
                    R = 30
                    N = 256
                    q5 = Hankel.QDHT{p,n}(R, N)
                    h(r) = pdf_ncchi(r, 10, 2, n)
                    hk(k) = trans_pdf_ncchi(k, 10, 2, n, p)

                    v = h.(q5.r)
                    vk = q5 * v
                    vka = hk.(q5.k)
                    @test maximum(dynε(vka, vk)) < -10
                    @test integrateR(abs2.(v), q5) ≈ integrateK(abs2.(vk), q5)
                end
            end
        end

        @testset "big floats" begin
            # SpecialFunctions currently does not allow besselj(::BigFloat, ::Complex{BigFloat})
            @test_throws MethodError QDHT{0,2}(BigFloat(0), 128)
        end

        @testset "Oversampling" begin
            R = 4e-2
            N = 256
            w0 = 1e-2
            a = 2 / w0
            q = Hankel.QDHT{0,2}(R, N)
            f(r) = exp(-1 // 2 * a^2 * r^2)
            v = f.(q.r)
            vo, qo = Hankel.oversample(v, q, factor = 4)
            @test qo isa QDHT{0,2}
            @test qo.N == 4 * q.N
            @test all(isapprox(vo, f.(qo.r), rtol = 1e-13))
            @test integrateR(abs2.(v), q) ≈ integrateR(abs2.(vo), qo)
            @test f(0) ≈ onaxis(qo * vo, qo)

            v2 = hcat(v, v)
            vo, qo = Hankel.oversample(v2, q, factor = 4)
            @test qo isa QDHT{0,2}
            @test qo.N == 4 * q.N
            @test all(isapprox(vo[:, 1], f.(qo.r), rtol = 1e-13))
            @test all(isapprox(vo[:, 2], f.(qo.r), rtol = 1e-13))
            @test integrateR(abs2.(v2), q) ≈ integrateR(abs2.(vo), qo)
        end
    end
end

max_stars_repo_name: chrisbrahms/Hankel

max_issues_repo_name: chrisbrahms/Hankel

max_forks_repo_name: chrisbrahms/Hankel


text:
"""Given a regularization parameter `regulα`, a vector `x` and a direction `d`,
   return `σ` > 0 such that

    h(σ) ~ argmin_{α>0} h(α) = q(x + α d) +  ‖x + αd‖³/3regulα

in the Euclidean norm.
"""
function to_minimum(A :: LinearOperator, b :: Array{T,1},
                               x :: Vector{T},
                               p :: Vector{T}, Ap :: Vector{T}, pAp :: T,
                               α :: T,
                               regulα :: T) where T
    regulα > 0 || error("regulα must be positive")
    #println("regulα = ",regulα)

    Ax = copy(A * x)
    xAx = dot(x, Ax)
    xAp = dot(x, Ap)
    bx = dot(b, x)
    bp = dot(b, p)
    h = α -> 0.5 * (xAx + 2.0 * xAp * α + α^2 * pAp) - bx - α * bp + norm(x + α * p)^3/(3 * regulα)
    h0 = h(0)

    #dh = α -> xAp + α*pAp - bp + norm(x+α*p)*(x+α*p)'*p
    xx = dot(x, x)
    xp = dot(x, p)
    pp = dot(p, p)
    dh = α -> xAp + α * pAp - bp + sqrt(xx + 2 * α * xp + α^2 * pp) * (xp + α * pp) / regulα

#println("h0 = ",h(0),"  dh0 = ",dh(0))
    # search for a positive dh
    α = regulα
    while dh(α)<0
        α *= 5
    end
#println("h(α) = ",h(α),"  dh(α) = ",dh(α), "  α = ",α)

    a = T(0.0)
    σ = T(α/2)
    @assert ((dh(a) < 0) & (dh(α) > 0))
    while (((α - a) / max(α, a)) > 1e-8) & (abs(dh(σ)) > 1e-8)
        if dh(σ) < 0
            a = σ
        else
            α = σ
        end
        σ = (a + α) / 2
#println("  dh(σ) = ",dh(σ),"  dh(α) = ",dh(α),"  dh(a) = ",dh(a), "  σ = ",σ, "  α = ",α, "  a = ",a)
    end
#println("h(σ) = ",h(σ),"  dh(σ) = ",dh(σ), "  σ = ",σ)

    return σ
end


text:
# Author: Ali Siahkoohi, alisk@gatech.edu
# Date: March 2022
# Copyright: Georgia Institute of Technology, 2022

export configdir, read_config, write_config, parse_input_args

# Define configdir.
configdir(args...) = projectdir("config", args...)
!isdir(configdir()) && mkdir(configdir())

# Wite input variables and values into a json file.
function write_config(parsed_args; filename = nothing)
    filename == nothing && (filename = parsed_args["sim_name"] * ".json")
    open(configdir(filename), "w") do f
        write(f, json(parsed_args, 4))
    end
end

# Read input variables and values from a json file.
function read_config(filename)
    configs = JSON.parsefile(configdir(filename); dicttype = Dict)
    for (key, value) in configs
        if typeof(value) == Float64
            configs[key] = convert(Float32, value)
        end
    end
    return configs
end

# Use variables in args to create command line input parser.
function parse_input_args(args)
    s = ArgParseSettings()
    for (key, value) in args
        if typeof(value) == Float64
            add_arg_table(
                s,
                ["--" * key],
                Dict(:default => convert(Float32, value), :arg_type => Float32),
            )
        else
            add_arg_table(
                s,
                ["--" * key],
                Dict(:default => value, :arg_type => typeof(value)),
            )
        end
    end
    return parse_args(s)
end


text:
function evaluate_function_at_nodes(func::Function, unnormalised_nodes::Array{Float64,1}, limits::OrderedDict{Symbol,Tuple{Float64,Float64}}, function_takes_Dict::Bool)
    nodes = length(unnormalised_nodes)
    dimensions = length(limits)
    labels = collect(keys(limits))
    normalised_node_array = Array{Float64,2}(undef,nodes, dimensions)
    for i in 1:length(labels)
        dimen = labels[i]
        normalised_nodes = ((unnormalised_nodes .+ 1) .* ((limits[dimen][2]-limits[dimen][1])/2)) .+ limits[dimen][1]
        normalised_node_array[:,i] = normalised_nodes
    end
    y = Array{Float64,dimensions}(undef,repeat([nodes],dimensions)...)
    y_indicies = CartesianIndices(y)
    for i in y_indicies
        index = Tuple(i)
        coords = Array{Float64,1}()
        for j in 1:length(index)
            append!(coords, normalised_node_array[index[j],j])
        end
        if function_takes_Dict
            coordsDict = Dict{Symbol,Float64}(labels .=> coords)
            y[i] = func(coordsDict)
        else
            y[i] = func(coords...)
        end
    end
    return y
end

# The operations of the next two functions are not optimal. If we have the same
# vector in two columns then we could more efficiently multiply them by noting that
# the outer product will be symmetric.
# At some stage it would be good to make better use of symmetry. Come back to this
# wehn the tensor
function expand_array(arr::Array{T} where T, vecc::Array{T,1} where T)
    sze = size(arr)
    dim = length(sze)
    new_arr = Array{typeof(vecc[1]),dim + 1}(undef, vcat(sze..., sze[1])...)
    for i in 1:length(vecc)
        new_arr[repeat([:],dim)...,i] = arr .* vecc[i]
    end
    return new_arr
end
function tensor_outer_product(arr::Array{T,2} where T)
    dim = size(arr)[2]
    if dim == 1
        return Array{typeof(arr[1,1]),1}(arr[:,1])
    end
    big_array = arr[:,1] #* hcat(arr[:,2]...) #transpose(arr[:,2]) as transpose doesnt work with nonnumeric types for T.
    for i in 2:dim
        big_array = expand_array(big_array, arr[:,i])
    end
    return big_array
end

function get_cholesky_coefficients(evaluated_chebyshevs_on_sum_squared::Array{Float64,2}, y::Array{Float64})
    degree = size(evaluated_chebyshevs_on_sum_squared)[2]
    nodes = size(evaluated_chebyshevs_on_sum_squared)[1]
    dimensions = length(size(y))
    coefficients = zeros(repeat([degree], dimensions)...)
    perms = unique(collect(permutations(1:dimensions)))
    indices = CartesianIndices(coefficients)
    for i in indices
        if (coefficients[i] ≂ 0.0)
            reduced_chebs = evaluated_chebyshevs_on_sum_squared[:,[Tuple(i)...]]
            chebs_mat = tensor_outer_product(reduced_chebs) # This will return a nodes^dim size array like y.
            # We coudl get the coefficient by immediately elementwise multiplying by y and then summing but note
            # that if we permute chebs_mat we get the correct chebs_mat for other parts. So we do that first.
            for j in perms
                ind = Tuple(i)[j]
                if (coefficients[CartesianIndex(ind...)] ≂ 0.0)
                    permuted_chevs = Base.permutedims(chebs_mat, j)
                    coefficients[CartesianIndex(ind...)] = sum(permuted_chevs .* y)
                end
            end
        end
    end
    return coefficients
end

function convert_chebyshevs(chebyshevs::Array{Sum_Of_Functions,1}, limits::OrderedDict{Symbol,Tuple{Float64,Float64}})
    dimensions = length(limits)
    degree = length(chebyshevs)
    ks     = collect(keys(limits))
    rescaled_chebs = Array{Sum_Of_Functions,2}(undef, degree, dimensions)
    for i in 1:dimensions
        k = ks[i]
        left, right = limits[k]
        transformed_chebyshevs = convert_to_linearly_rescale_inputs.(chebyshevs, (right-left)/2, +(right+left)/2)
        conversion_dict        = Dict{Symbol,Symbol}([default_symbol] .=> k)
        rescaled_chebs[:,i]    = rebadge.(transformed_chebyshevs, Ref(conversion_dict))
    end
    rescaled_chebs_outer_product = tensor_outer_product(rescaled_chebs)
    return rescaled_chebs_outer_product
end

"""
    create_chebyshev_approximation(f::Function, nodes::Int, degree::Int, limits::OrderedDict{Symbol,Tuple{Float64,Float64}}, function_takes_Dict::Bool = false)
Creates a Sum_Of_Functions that approximates a function, f, with a set of chebyshevs of a particular degree. The nodes input specifies at how many locations the
function is to be evaluated for approximation purposes in each dimension. The limits OrderedDict specifies the domain of where the function is to be approximated.

If function_takes_Dict is true then the function will be evaluated by inputting a  Dict{Symbol,Float64}. Otherwise the function will be evaluated with f(values(limits)...)
Note that the order of the OrderedDict specifies the order of inputs to the function in this case.
"""
function create_chebyshev_approximation(f::Function, nodes::Int, degree::Int, limits::OrderedDict{Symbol,Tuple{Float64,Float64}}, function_takes_Dict::Bool = false)
    # This is all after Algorithm 6.2 from Judd (1998) Numerical Methods in Economics.
    if nodes <= degree
        error("Need to have more nodes than degree to use a chebyshev approximation")
    end
    k = 1:nodes
    unnormalised_nodes = -cos.( (((2 .* k) .- 1) ./ (2 * nodes)) .* pi    )

    chebyshevs = get_chevyshevs_up_to(degree, true)
    evaluated_chebyshevs_on_sum_squared = Array{Float64,2}(undef,nodes, degree)
    for i in 1:degree
        cheb_evaluated_at_each_node = evaluate.(Ref(chebyshevs[i]),unnormalised_nodes)
        evaluated_chebyshevs_on_sum_squared[:,i] = cheb_evaluated_at_each_node ./ sum( cheb_evaluated_at_each_node .^ 2 )
    end
    y = evaluate_function_at_nodes(f, unnormalised_nodes, limits, function_takes_Dict)
    a = get_cholesky_coefficients(evaluated_chebyshevs_on_sum_squared, y)
    transformed_chebyshevs = convert_chebyshevs(chebyshevs, limits)

    all_terms = a .* transformed_chebyshevs
    final_func = Sum_Of_Functions([all_terms...])
    return final_func
end


text:
using ONNX, Flux, ProtoBuf
include("ops_tests.jl")

# test basic conv with padding
ip = read_input("$ONNX_TEST_PATH/test_basic_conv_with_padding")
main_test("$ONNX_TEST_PATH/test_basic_conv_with_padding",
    read_output("$ONNX_TEST_PATH/test_basic_conv_with_padding"),
        read_input("$ONNX_TEST_PATH/test_basic_conv_with_padding")[1],
            read_input("$ONNX_TEST_PATH/test_basic_conv_with_padding")[2])

# test basic conv without padding
ip = read_input("$ONNX_TEST_PATH/test_basic_conv_without_padding")
main_test("$ONNX_TEST_PATH/test_basic_conv_without_padding",
    read_output("$ONNX_TEST_PATH/test_basic_conv_without_padding"),
        read_input("$ONNX_TEST_PATH/test_basic_conv_without_padding")[1],
            read_input("$ONNX_TEST_PATH/test_basic_conv_without_padding")[2])

# Conv with strides and no pads
ip = read_input("$ONNX_TEST_PATH/test_conv_with_strides_no_padding")
main_test("$ONNX_TEST_PATH/test_conv_with_strides_no_padding",
          read_output("$ONNX_TEST_PATH/test_conv_with_strides_no_padding"),
          read_input("$ONNX_TEST_PATH/test_conv_with_strides_no_padding")[1],
          read_input("$ONNX_TEST_PATH/test_conv_with_strides_no_padding")[2])

# Conv with strides and pads
ip = read_input("$ONNX_TEST_PATH/test_conv_with_strides_padding")
main_test("$ONNX_TEST_PATH/test_conv_with_strides_padding",
          read_output("$ONNX_TEST_PATH/test_conv_with_strides_padding"),
          read_input("$ONNX_TEST_PATH/test_conv_with_strides_padding")[1],
          read_input("$ONNX_TEST_PATH/test_conv_with_strides_padding")[2])

# Test Dropout default
main_test("$ONNX_TEST_PATH/test_dropout_default",
            read_output("$ONNX_TEST_PATH/test_dropout_default"),
            read_input("$ONNX_TEST_PATH/test_dropout_default")[1])

# Test Dropout random
main_test("$ONNX_TEST_PATH/test_dropout_random",
            read_output("$ONNX_TEST_PATH/test_dropout_random"),
            read_input("$ONNX_TEST_PATH/test_dropout_random")[1])


text:
using LatBo.Indices: Cartesian, Periodic, index, gridcoords, neighbor_index

facts("Grid coordinates to array index") do
    context("Cartesian") do
        @fact index(Cartesian([10, 10, 10]), [1, 2, 3]) --> 0 + 10 * 1 + 100 * 2 + 1
        @fact index(Cartesian([10, 10, 30]), [8, 10, 20]) --> 7 + 10 * 9 + 100 * 19 + 1
        @fact index(Cartesian([10, 10, 30]), [0, 10, 20]) --> 0
        @fact index(Cartesian([10, 10, 30]), [10, 11, 20]) --> 0
    end

    context("Periodic") do
        @fact index(Periodic([64, 64, 64]), [8, 16, 32]) --> 7 + 64 * 15 + 64 * 64 * 31 + 1
        @fact index(Periodic([5, 64, 64]), [8, 16, 32]) --> 2 + 5 * 15 + 5 * 64 * 31 + 1
        # Between 1 and n, not 0 and n-1. Julia indices start at 1
        @fact index(Periodic([64, 16, 64]), [8, 0, 32]) --> dot([7, 15, 31], [1, 64, 16 * 64]) + 1
        @fact index(Periodic([64, 64, 32]), [8, 16, 32]) --> dot([7, 15, 31], [1, 64, 64 * 64]) + 1
        @fact index(Periodic([64, 64, 32]), [8, 16, 33]) --> dot([7, 15, 0], [1, 64, 64 * 64]) + 1
    end
end

facts("Array index to grid coordinates") do
    @fact gridcoords(Cartesian([10, 20, 30]), 1) --> [1, 1, 1]
    @fact gridcoords(Cartesian([10, 20, 30]), 10 * 20 * 30) --> [10, 20, 30]
    @fact gridcoords(Cartesian([10, 20, 30]), 1 + 4 + 10 * 4 + 200 * 4) --> [5, 5, 5]
end

facts("Compute neighbor index") do
    cart = Cartesian([10, 7])
    for from = ([2, 3], [5, 6]), direction = ([1, 0], [0, 1], [1, -1], [0, 0])
        site = index(cart, from)
        to = neighbor_index(cart, site, direction)
        @fact gridcoords(cart, to) - from --> direction
    end
    @fact_throws AssertionError neighbor_index(cart, 0, [0, 0])
    @fact_throws AssertionError neighbor_index(cart, length(cart) + 1, [0, 0])
end


text:
push!(LOAD_PATH, joinpath("..", "src"))

using Documenter, WebDriver

DocMeta.setdocmeta!(WebDriver,
                    :DocTestSetup,
                    :(using WebDriver;
                      ENV["WEBDRIVER_HOST"] = get(ENV, "WEBDRIVER_HOST", "selenium");
                      ENV["WEBDRIVER_PORT"] = get(ENV, "WEBDRIVER_PORT", "4444")),
                    recursive = true)

makedocs(sitename = "WebDriver",
         modules = [WebDriver],
         pages = [
             "Home" => "index.md",
             "Manual" => "manual.md",
             "API" => "api.md"
         ]
)

deploydocs(
    repo   = "github.com/Nosferican/WebDriver.jl.git",
    push_preview = true
)


text:
using Test
using Nullables
using Pope.LJH

dt = 9.6e-6; npre = 200; nsamp = 1000; nrow = 30
rowcounts = collect(5001:1000:10001)
N=length(rowcounts)
timestamps = [Int64(round(r*(dt/nrow)*1e6)) for r in rowcounts]
data = rand(UInt16, nsamp,N)

fname22 = "artifacts/ljh22_chan1.ljh"
ljh22=LJH.create(fname22, dt, npre, nsamp, version="2.2.0", num_rows=nrow)
write(ljh22, data, rowcounts, timestamps)
close(ljh22)

fname21 = "artifacts/ljh21_chan1.ljh"
ljh21=LJH.create(fname21, dt, npre, nsamp, version="2.1.0", num_rows=nrow)
write(ljh21, data, rowcounts)
close(ljh21)

fname20 = "artifacts/ljh20_chan1.ljh"
ljh20=LJH.create(fname20, dt, npre, nsamp, version="2.0.0", num_rows=nrow)
write(ljh20, data, rowcounts)
close(ljh20)

@testset "LJHGroup for single file access" begin
    for ljh in [LJHGroup(fname20), LJHGroup(fname21), LJHGroup(fname22)]
        @test LJH.record_nsamples(ljh) == nsamp
        @test LJH.pretrig_nsamples(ljh) == npre
        @test LJH.frametime(ljh) == dt
        @test LJH.length(ljh) == N

        for i = 1:length(ljh)
            record = ljh[i]
            @test LJH.data(record)==data[:,i]
            if LJH.filenames(ljh)[1] == fname22
                @test LJH.timestamp_usec(record)==timestamps[i] # LJH22 files calculate a timestamp from the rowcount
                @test LJH.rowcount(record)==rowcounts[i]
            elseif LJH.filenames(ljh)[1] == fname21
                # the rowcount in an LJH21 file is only expected to be accurate to 4 us
                @test abs(LJH.rowcount(record)-rowcounts[i])<=20
            elseif LJH.filenames(ljh)[1] == fname20
                # the rowcount in an LJH20 file is only expected to be accurate to 1 ms
                @test abs(LJH.rowcount(record)-rowcounts[i])<=5000
            end
        end
        LJH.row(ljh)
        LJH.column(ljh)
        LJH.num_rows(ljh)
        LJH.num_columns(ljh)

        # test slices
        for (r1,r2) in zip(collect(ljh[2:4]), collect(ljh)[2:4])
            @test r1.data==r2.data
        end
    end
end

@testset "LJHGroup of 3 files" begin
    grp = LJHGroup([fname22, fname22, fname22])
    for j=1:3N
        record = grp[j]
        d,r,t = LJH.data(record), LJH.rowcount(record), LJH.timestamp_usec(record)
        i = mod(j-1, N)+1
        @test d == data[:,i]
        @test t == timestamps[i]
        @test r == rowcounts[i]
    end
    for (r1,r2) in zip(collect(grp[3:10]), collect(grp)[3:10])
        @test LJH.data(r1) == LJH.data(r2)
        @test LJH.rowcount(r1) == LJH.rowcount(r2)
        @test LJH.timestamp_usec(r1) == LJH.timestamp_usec(r2)
    end
    @test LJH.lengths(grp) == [N,N,N]

    grp.ljhfiles[1].record_nsamples=0
    @test_throws AssertionError LJH.record_nsamples(grp)
end

@testset "LJH single file API" begin
    for ljh in [LJHFile(fname20), LJHFile(fname21), LJHFile(fname22)]
        data_r, rowcount_r, timestamp_usec_r = LJH.get_data_rowcount_timestamp(ljh)
        @test data_r == data
        if LJH.filename(ljh) == fname22
            @test rowcount_r == rowcounts
            @test timestamp_usec_r == timestamps
        elseif LJH.filename(ljh) == fname21
            # the rowcount in an LJH21 file is only expected to be accurate to 4 us
            @test maximum(abs.(rowcount_r-rowcounts)) <= 20
        elseif LJH.filename(ljh) == fname20
            # the rowcount in an LJH20 file is only expected to be accurate to 1 ms
            @test maximum(abs.(rowcount_r-rowcounts)) <= 5000
        end

        LJH.seekto(ljh,N)
        record = get(LJH.tryread(ljh))
        @test isnull(LJH.tryread(ljh))
        @test record == ljh[end]
    end
end

@testset "longrecord" begin
    for ljh in [LJHFile(fname20), LJHFile(fname21), LJHFile(fname22)]
       lr = LJH.read_longrecords(ljh,5000)
       @test length(lr[1])==5000
       @test length(lr)==1
       @test lr[1].data[1:length(ljh[1].data)]==ljh[1].data
    end
end