text:
using StaticArrays
using LinearAlgebra
include("/home/heindelj/dev/julia_development/wally/src/molecule_tools/units.jl")
include("constants.jl")

let # scope all the constants we need for this potential
    global pot_nasa, dms_nasa!

    ####### Partridge Schwenk indexing stuff #######
    idx1 = [
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2,
        2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
        2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
        3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3,
        3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
        4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4,
        4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6,
        6, 6, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5,
        6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 5, 5,
        5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,
        7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6,
        6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 9, 9,
        9, 9, 9, 9, 9]

    idx2 = [
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
        2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2,
        2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3,
        3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
        2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3,
        3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3,
        2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4,
        4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2,
        2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4,
        4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 1, 1,
        1, 1, 1, 1, 1]

    idx3 = [
        1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15, 1, 2, 3, 4, 5,
        6, 7, 8, 9,10,11,12,13,14, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,
        12,13, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13, 1, 2, 3, 4, 5,
        6, 7, 8, 9,10,11,12, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12, 1,
        2, 3, 4, 5, 6, 7, 8, 9,10,11, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,
        11, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11, 1, 2, 3, 4, 5, 6, 7, 8,
        9,10, 1, 2, 3, 4, 5, 6, 7, 8, 9,10, 1, 2, 3, 4, 5, 6, 7, 8,
        9,10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9,
        1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2,
        3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6,
        7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3,
        4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2,
        3, 4, 5, 6, 7]

    idxD0 = [
        1, 1, 1, 2, 1, 1, 1, 2, 2, 3, 1, 1, 1, 1, 2, 2, 2, 3, 3, 4,
        1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 4, 4, 5, 1, 1, 1, 1, 1,
        1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5, 6, 1, 1, 1, 1,
        1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5,
        5, 6, 6, 7]

    idxD1 = [
        1, 1, 2, 1, 1, 2, 3, 1, 2, 1, 1, 2, 3, 4, 1, 2, 3, 1, 2, 1,
        1, 2, 3, 4, 5, 1, 2, 3, 4, 1, 2, 3, 1, 2, 1, 1, 2, 3, 4, 5,
        6, 1, 2, 3, 4, 5, 1, 2, 3, 4, 1, 2, 3, 1, 2, 1, 1, 2, 3, 4,
        5, 6, 7, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 1, 2, 3, 4, 1, 2,
        3, 1, 2, 1]

    idxD2 = [
        1, 2, 1, 1, 3, 2, 1, 2, 1, 1, 4, 3, 2, 1, 3, 2, 1, 2, 1, 1,
        5, 4, 3, 2, 1, 4, 3, 2, 1, 3, 2, 1, 2, 1, 1, 6, 5, 4, 3, 2,
        1, 5, 4, 3, 2, 1, 4, 3, 2, 1, 3, 2, 1, 2, 1, 1, 7, 6, 5, 4,
        3, 2, 1, 6, 5, 4, 3, 2, 1, 5, 4, 3, 2, 1, 4, 3, 2, 1, 3, 2,
        1, 2, 1, 1]

    coefD = [
        -2.1689686086730e-03, 1.4910379754728e-02, 5.3546078430060e-02,
        -7.4055995388666e-02,-3.7764333017616e-03, 1.4089887256484e-01,
        -6.2584207687264e-02,-1.1260393113022e-01,-5.7824159269319e-02,
        1.4360743650655e-02,-1.5469680141070e-02,-1.3036350092795e-02,
        2.7515837781556e-02, 1.4098478875076e-01,-2.7663168397781e-02,
        -5.2378176254797e-03,-1.0237198381792e-02, 8.9571999265473e-02,
        7.2920263098603e-03,-2.6873260551686e-01, 2.0220870325864e-02,
        -7.0764766270927e-02, 1.2140640273760e-01, 2.0978491966341e-02,
        -1.9443840512668e-01, 4.0826835370618e-02,-4.5365190474650e-02,
        6.2779900072132e-02,-1.3194351021000e-01,-1.4673032718563e-01,
        1.1894031277247e-01,-6.4952851564679e-03, 8.8503610374493e-02,
        1.4899437409291e-01, 1.3962841511565e-01,-2.6459446720450e-02,
        -5.0128914532773e-02, 1.8329676428116e-01,-1.5559089125095e-01,
        -4.0176879767592e-02, 3.6192059996636e-01, 1.0202887240343e-01,
        1.9318668580051e-01,-4.3435977107932e-01,-4.2080828803311e-02,
        1.9144626027273e-01,-1.7851138969948e-01, 1.0524533875070e-01,
        -1.7954071602185e-02, 5.2022455612120e-02,-2.8891891146828e-01,
        -4.7452036576319e-02,-1.0939400546289e-01, 3.5916564473568e-01,
        -2.0162789820172e-01,-3.5838629543696e-01, 5.6706523551202e-03,
        1.3849337488211e-01,-4.1733982195604e-01, 4.1641570764241e-01,
        -1.2243429796296e-01, 4.7141730971228e-02,-1.8224510249551e-01,
        -1.8880981556620e-01,-3.1992359561800e-01,-1.8567550546587e-01,
        6.1850530431280e-01,-6.1142756235141e-02,-1.6996135584933e-01,
        5.4252879499871e-01, 6.6128603899427e-01, 1.2107016404639e-02,
        -1.9633639729189e-01, 2.7652059420824e-03,-2.2684111109778e-01,
        -4.7924491598635e-01, 2.4287790137314e-01,-1.4296023329441e-01,
        8.9664665907006e-02,-1.4003228575602e-01,-1.3321543452254e-01,
        -1.8340983193745e-01, 2.3426707273520e-01, 1.5141050914514e-01]

    function get_c5z()
    """
    The PS potential has a number of basis expansion coefficients which are combined
    to form a single vector of coefficients. This function just shows how the linear
    combination of coefficients is formed. This code is not actually used besides in
    initialization of the PS_Constants.
    """
    # scaling factors for contributions to emperical potential
    f5z::Float64 = 0.999677885;
    fbasis::Float64 = 0.15860145369897;
    fcore::Float64 = -1.6351695982132;
    frest::Float64 = 1.0;

    c5zA = [
            4.2278462684916e+04, 4.5859382909906e-02, 9.4804986183058e+03,
            7.5485566680955e+02, 1.9865052511496e+03, 4.3768071560862e+02,
            1.4466054104131e+03, 1.3591924557890e+02,-1.4299027252645e+03,
            6.6966329416373e+02, 3.8065088734195e+03,-5.0582552618154e+02,
        -3.2067534385604e+03, 6.9673382568135e+02, 1.6789085874578e+03,
        -3.5387509130093e+03,-1.2902326455736e+04,-6.4271125232353e+03,
        -6.9346876863641e+03,-4.9765266152649e+02,-3.4380943579627e+03,
            3.9925274973255e+03,-1.2703668547457e+04,-1.5831591056092e+04,
            2.9431777405339e+04, 2.5071411925779e+04,-4.8518811956397e+04,
        -1.4430705306580e+04, 2.5844109323395e+04,-2.3371683301770e+03,
            1.2333872678202e+04, 6.6525207018832e+03,-2.0884209672231e+03,
        -6.3008463062877e+03, 4.2548148298119e+04, 2.1561445953347e+04,
        -1.5517277060400e+05, 2.9277086555691e+04, 2.6154026873478e+05,
        -1.3093666159230e+05,-1.6260425387088e+05, 1.2311652217133e+05,
        -5.1764697159603e+04, 2.5287599662992e+03, 3.0114701659513e+04,
        -2.0580084492150e+03, 3.3617940269402e+04, 1.3503379582016e+04,
        -1.0401149481887e+05,-6.3248258344140e+04, 2.4576697811922e+05,
            8.9685253338525e+04,-2.3910076031416e+05,-6.5265145723160e+04,
            8.9184290973880e+04,-8.0850272976101e+03,-3.1054961140464e+04,
        -1.3684354599285e+04, 9.3754012976495e+03,-7.4676475789329e+04,
        -1.8122270942076e+05, 2.6987309391410e+05, 4.0582251904706e+05,
        -4.7103517814752e+05,-3.6115503974010e+05, 3.2284775325099e+05,
            1.3264691929787e+04, 1.8025253924335e+05,-1.2235925565102e+04,
        -9.1363898120735e+03,-4.1294242946858e+04,-3.4995730900098e+04,
            3.1769893347165e+05, 2.8395605362570e+05,-1.0784536354219e+06,
        -5.9451106980882e+05, 1.5215430060937e+06, 4.5943167339298e+05,
        -7.9957883936866e+05,-9.2432840622294e+04, 5.5825423140341e+03,
            3.0673594098716e+03, 8.7439532014842e+04, 1.9113438435651e+05,
        -3.4306742659939e+05,-3.0711488132651e+05, 6.2118702580693e+05,
        -1.5805976377422e+04,-4.2038045404190e+05, 3.4847108834282e+05,
        -1.3486811106770e+04, 3.1256632170871e+04, 5.3344700235019e+03,
            2.6384242145376e+04, 1.2917121516510e+05,-1.3160848301195e+05,
        -4.5853998051192e+05, 3.5760105069089e+05, 6.4570143281747e+05,
        -3.6980075904167e+05,-3.2941029518332e+05,-3.5042507366553e+05,
            2.1513919629391e+03, 6.3403845616538e+04, 6.2152822008047e+04,
        -4.8805335375295e+05,-6.3261951398766e+05, 1.8433340786742e+06,
            1.4650263449690e+06,-2.9204939728308e+06,-1.1011338105757e+06,
            1.7270664922758e+06, 3.4925947462024e+05,-1.9526251371308e+04,
        -3.2271030511683e+04,-3.7601575719875e+05, 1.8295007005531e+05,
            1.5005699079799e+06,-1.2350076538617e+06,-1.8221938812193e+06,
            1.5438780841786e+06,-3.2729150692367e+03, 1.0546285883943e+04,
        -4.7118461673723e+04,-1.1458551385925e+05, 2.7704588008958e+05,
            7.4145816862032e+05,-6.6864945408289e+05,-1.6992324545166e+06,
            6.7487333473248e+05, 1.4361670430046e+06,-2.0837555267331e+05,
            4.7678355561019e+05,-1.5194821786066e+04,-1.1987249931134e+05,
            1.3007675671713e+05, 9.6641544907323e+05,-5.3379849922258e+05,
        -2.4303858824867e+06, 1.5261649025605e+06, 2.0186755858342e+06,
        -1.6429544469130e+06,-1.7921520714752e+04, 1.4125624734639e+04,
        -2.5345006031695e+04, 1.7853375909076e+05,-5.4318156343922e+04,
        -3.6889685715963e+05, 4.2449670705837e+05, 3.5020329799394e+05,
            9.3825886484788e+03,-8.0012127425648e+05, 9.8554789856472e+04,
            4.9210554266522e+05,-6.4038493953446e+05,-2.8398085766046e+06,
            2.1390360019254e+06, 6.3452935017176e+06,-2.3677386290925e+06,
        -3.9697874352050e+06,-1.9490691547041e+04, 4.4213579019433e+04,
            1.6113884156437e+05,-7.1247665213713e+05,-1.1808376404616e+06,
            3.0815171952564e+06, 1.3519809705593e+06,-3.4457898745450e+06,
            2.0705775494050e+05,-4.3778169926622e+05, 8.7041260169714e+03,
            1.8982512628535e+05,-2.9708215504578e+05,-8.8213012222074e+05,
            8.6031109049755e+05, 1.0968800857081e+06,-1.0114716732602e+06,
            1.9367263614108e+05, 2.8678295007137e+05,-9.4347729862989e+04,
            4.4154039394108e+04, 5.3686756196439e+05, 1.7254041770855e+05,
        -2.5310674462399e+06,-2.0381171865455e+06, 3.3780796258176e+06,
            7.8836220768478e+05,-1.5307728782887e+05,-3.7573362053757e+05,
            1.0124501604626e+06, 2.0929686545723e+06,-5.7305706586465e+06,
        -2.6200352535413e+06, 7.1543745536691e+06,-1.9733601879064e+04,
            8.5273008477607e+04, 6.1062454495045e+04,-2.2642508675984e+05,
            2.4581653864150e+05,-9.0376851105383e+05,-4.4367930945690e+05,
            1.5740351463593e+06, 2.4563041445249e+05,-3.4697646046367e+03,
        -2.1391370322552e+05, 4.2358948404842e+05, 5.6270081955003e+05,
        -8.5007851251980e+05,-6.1182429537130e+05, 5.6690751824341e+05,
        -3.5617502919487e+05,-8.1875263381402e+02,-2.4506258140060e+05,
            2.5830513731509e+05, 6.0646114465433e+05,-6.9676584616955e+05,
            5.1937406389690e+05, 1.7261913546007e+05,-1.7405787307472e+04,
        -3.8301842660567e+05, 5.4227693205154e+05, 2.5442083515211e+06,
        -1.1837755702370e+06,-1.9381959088092e+06,-4.0642141553575e+05,
            1.1840693827934e+04,-1.5334500255967e+05, 4.9098619510989e+05,
            6.1688992640977e+05, 2.2351144690009e+05,-1.8550462739570e+06,
            9.6815110649918e+03,-8.1526584681055e+04,-8.0810433155289e+04,
            3.4520506615177e+05, 2.5509863381419e+05,-1.3331224992157e+05,
        -4.3119301071653e+05,-5.9818343115856e+04, 1.7863692414573e+03,
            8.9440694919836e+04,-2.5558967650731e+05,-2.2130423988459e+04,
            4.4973674518316e+05,-2.2094939343618e+05]

    # expansion coefficients for basis correction
        cbasis = [
    6.9770019624764e-04,-2.4209870001642e+01, 1.8113927151562e+01,
    3.5107416275981e+01,-5.4600021126735e+00,-4.8731149608386e+01,
    3.6007189184766e+01, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    -7.7178474355102e+01,-3.8460795013977e+01,-4.6622480912340e+01,
    5.5684951167513e+01, 1.2274939911242e+02,-1.4325154752086e+02,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00,-6.0800589055949e+00,
    8.6171499453475e+01,-8.4066835441327e+01,-5.8228085624620e+01,
    2.0237393793875e+02, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    3.3525582670313e+02, 7.0056962392208e+01,-4.5312502936708e+01,
    -3.0441141194247e+02, 2.8111438108965e+02, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00,-1.2983583774779e+02, 3.9781671212935e+01,
    -6.6793945229609e+01,-1.9259805675433e+02, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00,-8.2855757669957e+02,-5.7003072730941e+01,
    -3.5604806670066e+01, 9.6277766002709e+01, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 8.8645622149112e+02,-7.6908409772041e+01,
    6.8111763314154e+01, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    2.5090493428062e+02,-2.3622141780572e+02, 5.8155647658455e+02,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 2.8919570295095e+03,
    -1.7871014635921e+02,-1.3515667622500e+02, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00,-3.6965613754734e+03, 2.1148158286617e+02,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00,-1.4795670139431e+03,
    3.6210798138768e+02, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    -5.3552886800881e+03, 3.1006384016202e+02, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 1.6241824368764e+03, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 4.3764909606382e+03, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 1.0940849243716e+03, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 3.0743267832931e+03, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00]

    # expansion coefficients for core correction
        ccore = [
    2.4332191647159e-02,-2.9749090113656e+01, 1.8638980892831e+01,
    -6.1272361746520e+00, 2.1567487597605e+00,-1.5552044084945e+01,
    8.9752150543954e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    -3.5693557878741e+02,-3.0398393196894e+00,-6.5936553294576e+00,
    1.6056619388911e+01, 7.8061422868204e+01,-8.6270891686359e+01,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00,-3.1688002530217e+01,
    3.7586725583944e+01,-3.2725765966657e+01,-5.6458213299259e+00,
    2.1502613314595e+01, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    5.2789943583277e+02,-4.2461079404962e+00,-2.4937638543122e+01,
    -1.1963809321312e+02, 2.0240663228078e+02, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00,-6.2574211352272e+02,-6.9617539465382e+00,
    -5.9440243471241e+01, 1.4944220180218e+01, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00,-1.2851139918332e+03,-6.5043516710835e+00,
    4.0410829440249e+01,-6.7162452402027e+01, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 1.0031942127832e+03, 7.6137226541944e+01,
    -2.7279242226902e+01, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    -3.3059000871075e+01, 2.4384498749480e+01,-1.4597931874215e+02,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 1.6559579606045e+03,
    1.5038996611400e+02,-7.3865347730818e+01, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00,-1.9738401290808e+03,-1.4149993809415e+02,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00,-1.2756627454888e+02,
    4.1487702227579e+01, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    -1.7406770966429e+03,-9.3812204399266e+01, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00,-1.1890301282216e+03, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 2.3723447727360e+03, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00,-1.0279968223292e+03, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 5.7153838472603e+02, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    0.0000000000000e+00, 0.0000000000000e+00]

        crest = [
        0.0000000000000e+00,-4.7430930170000e+00,-1.4422132560000e+01,
    -1.8061146510000e+01, 7.5186735000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
    -2.7962099800000e+02, 1.7616414260000e+01,-9.9741392630000e+01,
        7.1402447000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00,-7.8571336480000e+01,
        5.2434353250000e+01, 7.7696745000000e+01, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        1.7799123760000e+02, 1.4564532380000e+02, 2.2347226000000e+02,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00,-4.3823284100000e+02,-7.2846553000000e+02,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00,-2.6752313750000e+02, 3.6170310000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
        0.0000000000000e+00, 0.0000000000000e+00]

        return f5z * c5zA + fbasis * cbasis + fcore * ccore + frest * crest
    end

    PS_Constants = (;
        reoh    = 0.958649,
        thetae  = 104.3475,
        b1      = 2.0,
        roh     = 0.9519607159623009,
        alphaoh = 2.587949757553683,
        deohA   = 42290.92019288289,
        phh1A   = 16.94879431193463,
        phh2    = 12.66426998162947,
        f5z     = 0.999677885,
        c5z     = get_c5z(),
        b1D     = 1.0,
        a       = 0.2999,
        b       = -0.6932,
        c0      = 1.0099,
        c1      = -0.1801,
        c2      = 0.0892,
        #ttm3 dms parameters
        dms_param1 = 0.5,
        dms_param2 = 0.9578,
        dms_param3 = 0.012
    )

    C = PS_Constants # just an alias for convenience

    function pot_nasa(coords::SMatrix{3, 3, Float64}, grads::Union{AbstractMatrix{Float64}, Nothing}=nothing)
        rOH1::SVector{3, Float64} = coords[:, 1] - coords[:, 2]
        rOH2::SVector{3, Float64} = coords[:, 1] - coords[:, 3]
        rHH ::SVector{3, Float64} = coords[:, 2] - coords[:, 3]

        cosθ::Float64 = rOH1 ⋅ rOH2 / (norm(rOH1) * norm(rOH2))

        deoh::Float64   = C.f5z * C.deohA
        phh1::Float64   = C.f5z * C.phh1A * exp(C.phh2)
        cosθ_e::Float64 = -.24780227221366464506
        exp1::Float64   = exp(-C.alphaoh * (norm(rOH1) - C.roh))
        exp2::Float64   = exp(-C.alphaoh * (norm(rOH2) - C.roh))
        Va::Float64     = deoh * (exp1 * (exp1 - 2.0) + exp2 * (exp2 - 2.0))
        Vb::Float64     = phh1 * exp(-C.phh2 * norm(rHH))
        dVa1::Float64   = 2.0 * C.alphaoh * deoh * exp1 * (1.0 - exp1) / norm(rOH1)
        dVa2::Float64   = 2.0 * C.alphaoh * deoh * exp2 * (1.0 - exp2) / norm(rOH2)
        dVb::Float64    = -C.phh2 * Vb / norm(rHH)
        x1::Float64     = (norm(rOH1) - C.reoh) / C.reoh
        x2::Float64     = (norm(rOH2) - C.reoh) / C.reoh
        x3::Float64     = cosθ - cosθ_e

        # SPEED: might be able to speed this up by making it 16,3?
        fmat = zeros(16, 3)
        for i in 1:3
            fmat[1, i] = 0.0
            fmat[2, i] = 1.0
        end

        for i in 3:16
            fmat[i, 1] = fmat[i-1, 1] * x1
            fmat[i, 2] = fmat[i-1, 2] * x2
            fmat[i, 3] = fmat[i-1, 3] * x3
        end

        efac::Float64 = exp(-C.b1 * ((norm(rOH1) - C.reoh)^2 + (norm(rOH2) - C.reoh)^2))

        sum0::Float64 = 0.0
        sum1::Float64 = 0.0
        sum2::Float64 = 0.0
        sum3::Float64 = 0.0

        for j in 1:244
            inI::Int = idx1[j+1]
            inJ::Int = idx2[j+1]
            inK::Int = idx3[j+1]

            sum0 += C.c5z[j+1]*(fmat[inI+1, 1] * fmat[inJ+1, 2] + fmat[inJ+1, 1]*fmat[inI+1, 2])*fmat[inK+1,3]

            sum1 += C.c5z[j+1]*((inI - 1)*fmat[inI,1]*fmat[inJ+1,2] + (inJ - 1)*fmat[inJ, 1]*fmat[inI+1, 2])*fmat[inK+1, 3]

            sum2 += C.c5z[j+1]*((inJ - 1)*fmat[inI+1, 1]*fmat[inJ, 2] + (inI - 1)*fmat[inJ+1, 1]*fmat[inI, 2])*fmat[inK+1, 3]

            sum3 += C.c5z[j+1]*(fmat[inI+1, 1]*fmat[inJ+1, 2] + fmat[inJ+1, 1]*fmat[inI+1, 2]) *(inK - 1)*fmat[inK, 3]
        end

        # Energy
        Vc::Float64 = 2 * C.c5z[1] + efac * sum0;
        energy::Float64 = Va + Vb + Vc;

        energy += 0.44739574026257 # correction
        energy *= conversion(:wavenumbers, :kcal)
        if grads !== nothing
            dVcdr1::Float64 = (-2 * C.b1 * efac * (norm(rOH1) - C.reoh) * sum0 + efac * sum1 / C.reoh) / norm(rOH1)

            dVcdr2::Float64 = (-2 * C.b1 * efac * (norm(rOH2) - C.reoh) * sum0 + efac * sum2 / C.reoh) / norm(rOH2)

            dVcdcth::Float64 = efac * sum3;

            # H1
            @views grads[:, 2] = dVa1 * rOH1 + dVb * rHH + dVcdr1 * rOH1 + dVcdcth * (rOH2 / (norm(rOH1) * norm(rOH2)) - cosθ * rOH1 / (norm(rOH1) * norm(rOH1)))
            # H2
            @views grads[:, 3] = dVa2 * rOH2 - dVb * rHH + dVcdr2 * rOH2 + dVcdcth * (rOH1 / (norm(rOH1) * norm(rOH2)) - cosθ * rOH2 / (norm(rOH2) * norm(rOH2)))
            # O
            @views grads[:, 1] = -(grads[:, 2] + grads[:, 3])
            grads .*= -conversion(:wavenumbers, :kcal) # also converts from forces to gradients
            return energy
        end

        return energy
    end

    function dms_nasa!(coords::SMatrix{3, 3, Float64}, q::MVector{3, Float64}, q_derivative::Union{MArray{Tuple{3, 3, 3}, Float64, 3, 27}, Nothing}=nothing, ttm3::Bool=false)
        rOH1::SVector{3, Float64} = coords[:, 1] - coords[:, 2]
        rOH2::SVector{3, Float64} = coords[:, 1] - coords[:, 3]

        cosθ::Float64 = rOH1 ⋅ rOH2 / (norm(rOH1) * norm(rOH2))

        cosθ_e::Float64 = -0.24780227221366464506
        aθ_0::Float64   =  1.82400520401572996557

        x1::Float64     = (norm(rOH1) - C.reoh) / C.reoh
        x2::Float64     = (norm(rOH2) - C.reoh) / C.reoh
        x3::Float64     = cosθ - cosθ_e

        efac::Float64 = exp(-C.b1D * ((norm(rOH1) - C.reoh)^2 + (norm(rOH2) - C.reoh)^2))

        fmat = zeros(16, 3)
        for i in 1:3
            fmat[1, i] = 0.0
            fmat[2, i] = 1.0
        end

        for i in 3:16
            fmat[i, 1] = fmat[i-1, 1] * x1
            fmat[i, 2] = fmat[i-1, 2] * x2
            fmat[i, 3] = fmat[i-1, 3] * x3
        end

        # Calculate the dipole moment

        p1::Float64 = 0.0
        p2::Float64 = 0.0
        pl1 = cosθ
        pl2 = 0.5*(3 * pl1 * pl1 - 1.0)

        dp1dr1::Float64 = 0.0
        dp1dr2::Float64 = 0.0
        dp1dcabc::Float64 = 0.0
        dp2dr1::Float64 = 0.0
        dp2dr2::Float64 = 0.0
        dp2dcabc::Float64 = 0.0

        for j in 1:83
            inI = idxD0[j+1]
            inJ = idxD1[j+1]
            inK = idxD2[j+1]

            p1 += coefD[j+1] * fmat[inI+1, 1] * fmat[inJ+1, 2] * fmat[inK+1, 3]
            p2 += coefD[j+1] * fmat[inJ+1, 1] * fmat[inI+1, 2] * fmat[inK+1, 3]

            if (q_derivative === nothing) # skip derivatives
                continue
            end

            dp1dr1   += coefD[j+1] * (inI-1) * fmat[inI, 1]         * fmat[inJ + 1, 2]     * fmat[inK + 1, 3]
            dp1dr2   += coefD[j+1] * (inJ-1) * fmat[inI + 1, 1]     * fmat[inJ, 2]         * fmat[inK + 1, 3]
            dp1dcabc += coefD[j+1] * (inK-1) * fmat[inI + 1, 1]     * fmat[inJ + 1, 2]     * fmat[inK, 3]
            dp2dr1   += coefD[j+1] * (inJ-1) * fmat[inJ, 1]         * fmat[inI + 1, 2]     * fmat[inK + 1, 3]
            dp2dr2   += coefD[j+1] * (inI-1) * fmat[inJ + 1, 1]     * fmat[inI, 2]         * fmat[inK + 1, 3]
            dp2dcabc += coefD[j+1] * (inK-1) * fmat[inJ + 1, 1]     * fmat[inI + 1, 2]     * fmat[inK, 3]
        end

        xx::Float64  = 5.291772109200e-01
        xx2::Float64 = xx * xx

        dp1dr1 /= (C.reoh / xx)
        dp1dr2 /= (C.reoh / xx)
        dp2dr1 /= (C.reoh / xx)
        dp2dr2 /= (C.reoh / xx)

        pc0::Float64 = C.a * (norm(rOH1)^C.b + norm(rOH2)^C.b) * (C.c0 + pl1 * C.c1 + pl2 * C.c2)

        dpc0dr1::Float64 = C.a * C.b * norm(rOH1)^(C.b - 1) * (C.c0 + pl1 * C.c1 + pl2 * C.c2) * xx2
        dpc0dr2::Float64 = C.a * C.b * norm(rOH2)^(C.b - 1) * (C.c0 + pl1 * C.c1 + pl2 * C.c2) * xx2
        dpc0dcabc::Float64 = C.a * (norm(rOH1)^C.b + norm(rOH2)^C.b) * (C.c1 + 0.5 * (6.0 * pl1) * C.c2) * xx

        defacdr1::Float64 = -2.0 * C.b1D * (norm(rOH1) - C.reoh) * efac * xx
        defacdr2::Float64 = -2.0 * C.b1D * (norm(rOH2) - C.reoh) * efac * xx

        dp1dr1 = dp1dr1*efac + p1*defacdr1 + dpc0dr1;
        dp1dr2 = dp1dr2*efac + p1*defacdr2 + dpc0dr2;
        dp1dcabc = dp1dcabc*efac + dpc0dcabc;
        dp2dr1 = dp2dr1*efac + p2*defacdr1 + dpc0dr1;
        dp2dr2 = dp2dr2*efac + p2*defacdr2 + dpc0dr2;
        dp2dcabc = dp2dcabc*efac + dpc0dcabc;

        p1 = coefD[1] + p1 * efac + pc0 * xx # q^H1 in TTM2-F
        p2 = coefD[1] + p2 * efac + pc0 * xx # q^H2 paper

        q[1] = -(p1 + p2)  # Oxygen
        q[2] = p1          # Hydrogen-1
        q[3] = p2          # Hydrogen-2

        dp1dr1 /= xx
        dp1dr2 /= xx
        dp2dr1 /= xx
        dp2dr2 /= xx

        if (ttm3)

            #---------------------------------------------------------------------
            #........ Modification of the gas-phase dipole moment surface.........
            #---------------------------------------------------------------------

            #AxB = @MVector zeros(3)

            AxB = @SVector ([rOH1[2] * rOH2[3] - rOH1[3] * rOH2[2],
                            -rOH1[1] * rOH2[3] + rOH1[3] * rOH2[1],
                             rOH1[1] * rOH2[2] - rOH1[2] * rOH2[1]])

            sum0::Float64 = sum(AxB.^2)

            sinθ::Float64 = sqrt(sum0) / (norm(rOH1) * norm(rOH2))
            angle::Float64   = atan(sinθ, cosθ)

            p1 = C.dms_param1 * (norm(rOH1) - C.dms_param2) + C.dms_param3 * (angle - aθ_0)
            p2 = C.dms_param1 * (norm(rOH2) - C.dms_param2) + C.dms_param3 * (angle - aθ_0)

            q[1] = q[1] - (p1 + p2)
            q[2] = q[2] + p1
            q[3] = q[3] + p2

            dp1dr1 += C.dms_param1
            dp2dr2 += C.dms_param1
            dp1dcabc -= C.dms_param3 / sinθ
            dp2dcabc -= C.dms_param3 / sinθ

        end # ttm3

        if (q_derivative === nothing)
            return
        end

        f1q1r13::Float64 = (dp1dr1 - (dp1dcabc * cosθ / norm(rOH1))) / norm(rOH1)
        f1q1r23::Float64 = dp1dcabc/(norm(rOH1) * norm(rOH2))
        f2q1r23::Float64 = (dp1dr2 - (dp1dcabc * cosθ / norm(rOH2))) / norm(rOH2)
        f2q1r13::Float64 = dp1dcabc/(norm(rOH1) * norm(rOH2))
        f1q2r13::Float64 = (dp2dr1 - (dp2dcabc * cosθ / norm(rOH1))) / norm(rOH1)
        f1q2r23::Float64 = dp2dcabc/(norm(rOH1) * norm(rOH2))
        f2q2r23::Float64 = (dp2dr2 - (dp2dcabc * cosθ / norm(rOH2))) / norm(rOH2)
        f2q2r13::Float64 = dp2dcabc / (norm(rOH1) * norm(rOH2))

        # first index is atom w.r.t. to which the derivative is taken
        # second index is the charge being differentiated

        # gradient of charge h1(second index) wrt displacement of h1(third index)
        q_derivative[1, 1, 1] = f1q1r13*rOH1[1] + f1q1r23*rOH2[1]
        q_derivative[2, 1, 1] = f1q1r13*rOH1[2] + f1q1r23*rOH2[2]
        q_derivative[3, 1, 1] = f1q1r13*rOH1[3] + f1q1r23*rOH2[3]

        # gradient of charge h1 wrt displacement of h2
        q_derivative[1, 1, 2] = f2q1r13*rOH1[1] + f2q1r23*rOH2[1]
        q_derivative[2, 1, 2] = f2q1r13*rOH1[2] + f2q1r23*rOH2[2]
        q_derivative[3, 1, 2] = f2q1r13*rOH1[3] + f2q1r23*rOH2[3]

        # gradient of charge h1 wrt displacement of O
        q_derivative[1, 1, 3] = -(q_derivative[1, 1, 1] + q_derivative[1, 1, 2])
        q_derivative[2, 1, 3] = -(q_derivative[2, 1, 1] + q_derivative[2, 1, 2])
        q_derivative[3, 1, 3] = -(q_derivative[3, 1, 1] + q_derivative[3, 1, 2])

        # gradient of charge h2 wrt displacement of h1
        q_derivative[1, 2, 1] = f1q2r13*rOH1[1] + f1q2r23*rOH2[1]
        q_derivative[2, 2, 1] = f1q2r13*rOH1[2] + f1q2r23*rOH2[2]
        q_derivative[3, 2, 1] = f1q2r13*rOH1[3] + f1q2r23*rOH2[3]

        # gradient of charge h2 wrt displacement of h2
        q_derivative[1, 2, 2] = f2q2r13*rOH1[1] + f2q2r23*rOH2[1]
        q_derivative[2, 2, 2] = f2q2r13*rOH1[2] + f2q2r23*rOH2[2]
        q_derivative[3, 2, 2] = f2q2r13*rOH1[3] + f2q2r23*rOH2[3]

        # gradient of charge h2 wrt displacement of O
        q_derivative[1, 2, 3] = -(q_derivative[1, 2, 1] + q_derivative[1, 2, 2])
        q_derivative[2, 2, 3] = -(q_derivative[2, 2, 1] + q_derivative[2, 2, 2])
        q_derivative[3, 2, 3] = -(q_derivative[3, 2, 1] + q_derivative[3, 2, 2])

        # gradient of charge O wrt displacement of h1
        q_derivative[1, 3, 1] = -(q_derivative[1, 1, 1] + q_derivative[1, 2, 1])
        q_derivative[2, 3, 1] = -(q_derivative[2, 1, 1] + q_derivative[2, 2, 1])
        q_derivative[3, 3, 1] = -(q_derivative[3, 1, 1] + q_derivative[3, 2, 1])

        # gradient of charge O wrt displacement of h2
        q_derivative[1, 3, 2] = -(q_derivative[1, 1, 2] + q_derivative[1, 2, 2])
        q_derivative[2, 3, 2] = -(q_derivative[2, 1, 2] + q_derivative[2, 2, 2])
        q_derivative[3, 3, 2] = -(q_derivative[3, 1, 2] + q_derivative[3, 2, 2])

        # gradient of charge O wrt displacement of O
        q_derivative[1, 3, 3] = -(q_derivative[1, 1, 3] + q_derivative[1, 2, 3])
        q_derivative[2, 3, 3] = -(q_derivative[2, 1, 3] + q_derivative[2, 2, 3])
        q_derivative[3, 3, 3] = -(q_derivative[3, 1, 3] + q_derivative[3, 2, 3])
        # for some reason I'm off by a minus sign compared to reference so multiply by -1...
        @views q_derivative[:, :, :] *= -1
        return
    end
end # end of local scope

max_stars_repo_name: heindelj/MolecularPotentials

max_issues_repo_name: heindelj/MolecularPotentials

max_forks_repo_name: heindelj/MolecularPotentials


text:
for p in ("Knet","ArgParse", "Compat", "GZip", "Images","ImageCore")
    Pkg.installed(p) == nothing && Pkg.add(p)
end

using Knet, Compat,GZip, ArgParse, Images, ImageCore

include(Pkg.dir("Knet","data","mnist.jl"))

function main(args=ARGS)
    s = ArgParseSettings()
    s.description="Generative Adversarial Networks Knet Implementation."
    s.exc_handler=ArgParse.debug_handler
    @add_arg_table s begin
        ("--seed"; arg_type=Int; default=-1; help="random number seed: use a nonnegative int for repeatable results")
        ("--batchsize"; arg_type=Int; default=32; help="minibatch size")
        ("--lr"; arg_type=Float64; default=0.2; help="learning rate")
        ("--k"; arg_type=Int; default=1; help="k")
        ("--geninputsize"; arg_type=Int; default=784; help="size of the generator's input")
        ("--genhidden"; arg_type=Int; default=512; help="sizes of the generator hidden layer")
        ("--dischidden"; arg_type=Int; default=256; help="sizes of the discriminator hidden layer")
        ("--fast"; action=:store_true; help="skip loss printing for faster run")
        ("--epochs"; arg_type=Int; default=3; help="number of epochs for training")
        ("--iters"; arg_type=Int; default=typemax(Int); help="maximum number of updates for training")
        ("--gcheck"; arg_type=Int; default=0; help="check N random gradients per parameter")
        ("--atype"; default=(gpu()>=0 ? "KnetArray{Float32}" : "Array{Float32}"); help="array and float type to use")
    end
    isa(args, AbstractString) && (args=split(args))
    if in("--help", args) || in("-h", args)
        ArgParse.show_help(s; exit_when_done=false)
        return
    end
    println(s.description)
    o = parse_args(args, s; as_symbols=true)
    println("o=",[(k,v) for (k,v) in o]...)
    o[:seed] > 0 && srand(o[:seed])
    global atype = eval(parse(o[:atype]))
    # if atype <: Array; warn("CPU conv4 support is experimental and very slow."); end

    xtrn,ytrn,xtst,ytst = mnist()
    global dtrn = minibatch(xtrn, ytrn, o[:batchsize]; xtype=atype)
    global dtst = minibatch(xtst, ytst, o[:batchsize]; xtype=atype)
    imsize = 28*28

    w = initweights(o[:genhidden], o[:geninputsize], imsize)
    discw = initweights(o[:dischidden], imsize, 1)
    append!(w, discw)

    println(avgloss(w,dtst,o))
    train!(w, dtrn, o)
    println(avgloss(w,dtst,o))
    generate(w,10,o)

    return w,o,dtrn,dtst
end

function generate(w,ninstance,o;genfolder="gen/")
    z = samplez(o[:geninputsize],ninstance)
    gz = generator(w[1:4],z)
    gz = reshape(gz,(28,28,ninstance))
    gz = gz .> 0.5
    gz = permutedims(gz,(2,1,3))
    # imtype = Array{Gray{N0f8},2}
    for i=1:ninstance
        # save(genfolder*string(i)*".png",convert(imtype,gz[:,:,i]))
        # A = rand(3, 10, 10)
        img = colorview(Gray, convert(Array{Float32,2},gz[:,:,i]))
        save(genfolder*string(i)*".png", img)
    end
end

function avgloss(model,data,o)
    b = o[:batchsize]
    geninputsize = o[:geninputsize]
    total=dloss=gloss=0.0
    for (x,_) in data
        for k = 1:o[:k]
            z = samplez(geninputsize,b)
            gz = generator(model[1:4], z)
            dloss += 2*b*discloss(model[5:end],x,gz)
        end
        z = samplez(geninputsize,2*b)
        gloss += 2*b*genloss(model[1:4], model[5:end], z)
        total += 2*b
    end

    return dloss/total,gloss/total
end

function train!(model, data, o)
    b = o[:batchsize]
    geninputsize = o[:geninputsize]
    for i=1:o[:epochs]
        for (x,_) in data
            for k = 1:o[:k]
                z = samplez(geninputsize,b)
                gz = generator(model[1:4], z)
                g = discgradient(model[5:end],x,gz)
                update!(model[5:end], g, lr = o[:lr])
            end
            z = samplez(geninputsize,2*b)
            g = gengradient(model[1:4], model[5:end], z)
            update!(model[1:4], g, lr = o[:lr])
        end
        println(avgloss(model,dtrn,o))
    end
    # @save "genparams.jld2" model
    # save("ganparam.jld",model)
end

function discloss(w,x,gz)
    dx = discriminator(w,x)
    dz = discriminator(w,gz)
    loss = log.(1e-9 .+ dx) + log.(1e-9 + 1 .- dz)
    return -mean(loss)/2.0
end
discgradient = grad(discloss)

function genloss(wgen, wdisc, z)
    gz = generator(wgen,z)
    dz = discriminator(wdisc,gz)
    loss = -log.(1e-9 .+ dz)
    return mean(loss)
end
gengradient = grad(genloss)

function discriminator(w,x)
    for i=1:2:length(w)
       x = w[i]*mat(x) .+ w[i+1]
       if i<length(w)-1
           x = leakyrelu(x) # max(0,x)
       end
   end
   return sigm.(x)
end

function generator(w,z)
    for i=1:2:length(w)
       z = w[i]*mat(z) .+ w[i+1]
       if i<length(w)-1
           z = leakyrelu(z) # max(0,x)
       end
   end
   return sigm.(z)
end

function leakyrelu(x;alpha = 0.2)
    return max.(0,x) + min.(0,x) * alpha
end

function initweights(hidden, ninputs, noutputs; winit=1)
    w = Any[]
    x = ninputs
    for y in [hidden...,noutputs]
        push!(w, convert(atype, winit*xavier(y,x)))
        push!(w, convert(atype, zeros(y, 1)))
        x = y
    end
    return w
end

function samplez(geninputsize, batch)
    convert(atype, randn(Float32, geninputsize, batch))
end

main(ARGS)

ext: jl

max_stars_repo_name: feyzatoksal/generative-adversarial-nets

max_issues_repo_name: feyzatoksal/generative-adversarial-nets

max_forks_repo_name: feyzatoksal/generative-adversarial-nets


text:
import Base.Test: @test

using SVMLightLoader

println("Testing isnumeric")

@test SVMLightLoader.isnumeric("4.0")
@test SVMLightLoader.isnumeric("0")
@test SVMLightLoader.isnumeric(".3")
@test !SVMLightLoader.isnumeric("#")
@test !SVMLightLoader.isnumeric("A")

println("Testing line_to_data")

line_to_data = SVMLightLoader.line_to_data

# when the format is invalid
try line_to_data(" #comment") catch err @test isa(err, NoDataException) end
try line_to_data("# comment") catch err @test isa(err, NoDataException) end
try line_to_data("\n") catch err @test isa(err, NoDataException) end
try line_to_data("-1 2:1.0 5:") catch err @test isa(err, InvalidFormatError) end
try line_to_data("-1 :3") catch err @test isa(err, InvalidFormatError) end
try line_to_data("A") catch err @test isa(err, InvalidFormatError) end

(indices, values), label = line_to_data("-1")
@test (indices, values) == (Int64[], Float64[])
@test label == -1

(indices, values), label = line_to_data("-1 #comment")
@test (indices, values) == (Int64[], Float64[])
@test label == -1

# when the format is correct
(indices, values), label = line_to_data("-1 2:1.0 5:3.0 #comment")
@test (indices, values) == ([2, 5], [1.0, 3.0])
@test label == -1

(indices, values), label = line_to_data("2 2:1 5:3 #comment")
@test (indices, values) == ([2, 5], [1.0, 3.0])
@test label == 2

I = [2, 10, 15, 5, 12, 20]
J = [1, 1, 1, 2, 2, 3]
V = [2.5, -5.2, 1.5, 1.0, -3.0, 27.0]
X = sparse(I, J, V)

y = [1.0, 2.0, 3.0]

ndim = 30
Xndim = sparse(I, J, V, ndim, maximum(J))

println("Testing load_svmlight_file")

vectors, labels = load_svmlight_file("test.txt")
@test vectors == X
@test labels == y

vectors, labels = load_svmlight_file("test.txt", ndim)
@test vectors == Xndim
@test labels == y

error_thrown = false
try
    load_svmlight_file("invalid.txt")
catch error
    error_thrown = true
    @test isa(error, InvalidFormatError)
end
@test error_thrown

println("Testing iteration of SVMLightFile")

i = 0
for (vector, label) in SVMLightFile("test.txt")
    i += 1
    @test findnz(vector) == findnz(X[:, i])
    @test label == y[i]
end
@test i == size(X, 2)

i = 0
for (vector, label) in SVMLightFile("test.txt", ndim)
    i += 1
    @test findnz(vector) == findnz(X[:, i])
    @test label == y[i]
end
@test i == size(X, 2)

# empty.txt contains only newlines and comments
i = 1
for (vector, label) in SVMLightFile("empty.txt")
    i += 1
end

@test i == 1  # nothing should be loaded

error_thrown = false
try
    for (vector, label) in SVMLightFile("invalid.txt") end
catch error
    error_thrown = true
    @test isa(error, InvalidFormatError)
end
@test error_thrown

println("Testing length(s::SVMLightFile)")

@test length(SVMLightFile("test.txt")) == size(X, 2)

println("Testing List Comprehensions")

list = [label for (vector, label) in SVMLightFile("empty.txt")]
@test list == []

list = [label for (vector, label) in SVMLightFile("test.txt")]
@test list == [1.0, 2.0, 3.0]


text:
# Functions for Bayesian probit model with a Normal prior N(0, priorVar*I)

function logPrior(pars::Vector{Float64}, nPars::Int, data::Dict{Any, Any})
  return (-dot(pars,pars)/data["priorVar"]
    -nPars*log(2*pi*data["priorVar"]))/2
end

function logLikelihood(pars::Vector{Float64}, nPars::Int, data::Dict{Any, Any})
  XPars = data["X"]*pars
  return (dot(logOfNormalCdf(XPars), data["y"])
    +dot(logOfNormalCdf(-XPars), 1-data["y"]))
end

function gradLogPosterior(pars::Vector{Float64}, nPars::Int,
  data::Dict{Any, Any})
  XPars = data["X"]*pars
  return (data["X"]'
    *(data["y"].*exp(-(XPars.^2+log(2*pi))/2-logOfNormalCdf(XPars))
    -(1-data["y"]).*exp(-(XPars.^2+log(2*pi))/2-logOfNormalCdf(-XPars)))
    -pars/data["priorVar"])
end

function tensor(pars::Vector{Float64}, nPars::Int, data::Dict{Any, Any})
  XPars = data["X"]*pars
  vector =
    exp(-XPars.^2-logOfNormalCdf(XPars)-logOfNormalCdf(-XPars)-log(2*pi))
  return ((data["X"]'.*repmat(vector', nPars, 1))*data["X"]
    +(eye(nPars)/data["priorVar"]))
end

function derivTensor(pars::Vector{Float64}, nPars::Int, data::Dict{Any, Any})
  output = Array(Float64, nPars, nPars, nPars)

  XPars = data["X"]*pars
  vector01 =
    exp(-XPars.^2-2*logOfNormalCdf(XPars)-logOfNormalCdf(-XPars)-log(2*pi))

  for i = 1:nPars
    vector02 = (vector01.*(exp(-(XPars.^2+log(2*pi))/2-logOfNormalCdf(-XPars))
      -2*(pdf(Normal(0., 1.), XPars)
      +XPars.*cdf(Normal(0., 1.), XPars))).*data["X"][:, i])

    output[:, :, i] = (data["X"]'.*repmat(vector02', nPars, 1))*data["X"]
  end

  return output
end

function randPrior(nPars::Int, data::Dict{Any, Any})
  return rand(Normal(0.0, sqrt(data["priorVar"])), nPars)
end


text:
module WeakGradTestsSeq
   using Gridap
   using GridapGeosciences
   import Gridap.Fields: ∇, divergence
   using Plots
   using Test

   include("../ConvergenceAnalysisTools.jl")
   include("../WeakGradTests.jl")

   n_values=generate_n_values(2)
   hs=[2.0/n for n in n_values]
   model_args_series=[(n,) for n in n_values]

   @time ahs0,ak0errors,as0=convergence_study(compute_error_weak_grad,hs,model_args_series,0,4)
   @test as0 ≈ 0.8346320885900106

   n_values=generate_n_values(2,n_max=50)
   hs=[2.0/n for n in n_values]
   model_args_series=[(n,) for n in n_values]

   @time ahs1,ak1errors,as1=convergence_study(compute_error_weak_grad,hs,model_args_series,1,8)
   @test as1 ≈ 1.1034326200306834

end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
using .Gym
env = CartPoleEnv()
ctx = Ctx(env)

display(ctx.s)

#using Blink # when not on Juno
#body!(Blink.Window(), ctx.s)

actions = [sample(env.action_space) for i=1:1000]
i = 1
done = false
reset!(env)
while i <= length(actions) && !done
    global i, done
    a, b, done, d = step!(env, actions[i])

    render(env, ctx)
    i += 1
    sleep(0.1) # to see an animation
end


text:
using Polynomials, LinearMaps
using LinearAlgebra, SparseArrays
mutable struct CustomBOSolver{T} <: COSMO.AbstractKKTSolver
    k::Int
    σ::T
    ρ1::T
    ρ2::T
    n::Int
    m::Int
    sqrt_dim::Int
    A1::SparseMatrixCSC{T, Int}
    A::SparseMatrixCSC{T, Int}
    P::SparseMatrixCSC{T, Int}
    tmp_n::Vector{T}

    function CustomBOSolver(P::AbstractMatrix{T}, A::AbstractMatrix{T}, σ::T, ρ) where T
        n = size(P, 1)
        m = size(A, 1) - n
        ρ1 = ρ[1];
        ρ2 = ρ[end];
        k = Int(round(real(roots(Poly([-n, 0, 1/2, 1/2]))[end])))
        @assert k*(k + 1)*k/2 == n
        sqrt_dim = Int(k*(k + 1)/2)
        A1 = kron(ones(1, k), SparseMatrixCSC{T, Int}(I, sqrt_dim, sqrt_dim))
        new{T}(k, σ, ρ1, ρ2, n, m, Int(k*(k + 1)/2), A1, A, P, zeros(T, n))
    end
end

function mul_A1!(y, x, S)
    y .= 0
    idx = 1
    @inbounds for i = 1:S.k, j = 1:S.sqrt_dim
        y[j] += x[idx]
        idx += 1
    end
    #@assert norm(S.A1*x - y) <= 1e-8
end

function mul_A1t!(y, x, S)
    idx = 1
    @inbounds for i = 1:S.k, j = 1:S.sqrt_dim
        y[idx] = x[j]
        idx += 1
    end
    #@assert norm(S.A1'*x - y) <= 1e-8
end

function solve!(S::CustomBOSolver, y::AbstractVector{T}, x::AbstractVector{T}) where {T}
    # Solves the (KKT) linear system
    # | P + σI     A'  |  |y1|  =  |x1|
    # | A        -I/ρ  |  |y23|  =  |x23|
    # x1,y1: R^n, x2/y2: R^m
    # where [y1; y2] := y, [x1; x2] := x
    y1 = view(y, 1:S.n); y2 = view(y, S.n+1:S.n+S.m); y3 = view(y, S.n+S.m+1:2*S.n+S.m);
    x1 = view(x, 1:S.n); x2 = view(x, S.n+1:S.n+S.m); x3 = view(x, S.n+S.m+1:2*S.n+S.m);

    #y2 .= (S.A1*(x1 + x3*S.ρ2)/(S.ρ2 + S.σ) - x2)./(1/S.ρ1 + S.k/(S.ρ2 + S.σ))
    #y1 .= (-S.A1'*y2 + x1 + x3*S.ρ2)/(S.ρ2 + S.σ)
    #y3 .= S.ρ2* (y1 - x3)

    @. S.tmp_n = (x1 + x3*S.ρ2)/(S.ρ2 + S.σ);
    mul_A1!(y2, S.tmp_n, S)
    @. y2 = (y2 - x2)/(1/S.ρ1 + S.k/(S.ρ2 + S.σ))
    mul_A1t!(y1, y2, S)
    @. y1 = (-y1 + x1 + x3*S.ρ2)/(S.ρ2 + S.σ)
    @. y3 = S.ρ2* (y1 - x3)

    #=
    F = [S.P + S.σ*I SparseMatrixCSC(S.A'); S.A I]
    @inbounds @simd for i = (S.n + 1):size(F, 2)
        if i <= S.n + S.m
            ρ = S.ρ1
        else
            ρ = S.ρ2
        end
        F[i, i] = -1.0 / ρ
    end
    # show(stdout, "text/plain", F); println()
    @assert norm(y - F\x) <= 1e-7
    =#
end

function update_rho!(S::CustomBOSolver, ρ)
    S.ρ1 = ρ[1]
    S.ρ2 = ρ[end]
end


text:
using Test

# DOCUMENTER_TEST_EXAMPLES can be used to control which builds are performed in
# make.jl. But for the tests we need to make sure that all the relevant builds
# ran.
haskey(ENV, "DOCUMENTER_TEST_EXAMPLES") && error("DOCUMENTER_TEST_EXAMPLES env. variable is set")

# When the file is run separately we need to include make.jl which actually builds
# the docs and defines a few modules that are referred to in the docs. The make.jl
# has to be expected in the context of the Main module.
if (@__MODULE__) === Main && !@isdefined examples_root
    include("make.jl")
elseif (@__MODULE__) !== Main && isdefined(Main, :examples_root)
    using Documenter
    const examples_root = Main.examples_root
elseif (@__MODULE__) !== Main && !isdefined(Main, :examples_root)
    error("examples/make.jl has not been loaded into Main.")
end

@testset "Examples" begin
    @testset "HTML: deploy" begin
        doc = Main.examples_html_doc

        @test isa(doc, Documenter.Documents.Document)

        let build_dir = joinpath(examples_root, "builds", "html")
            @test joinpath(build_dir, "index.html") |> isfile
            @test joinpath(build_dir, "omitted", "index.html") |> isfile
            @test joinpath(build_dir, "hidden", "index.html") |> isfile
            @test joinpath(build_dir, "lib", "autodocs", "index.html") |> isfile
            @test joinpath(build_dir, "man", "style", "index.html") |> isfile

            # Test existence of some HTML elements
            man_style_html = String(read(joinpath(build_dir, "man", "style", "index.html")))
            @test occursin("is-category-myadmonition", man_style_html)

            # Assets
            @test joinpath(build_dir, "assets", "documenter.js") |> isfile

            # This build includes erlang and erlang-repl highlighting
            documenterjs = String(read(joinpath(build_dir, "assets", "documenter.js")))
            @test occursin("languages/julia.min", documenterjs)
            @test occursin("languages/julia-repl.min", documenterjs)
            @test occursin("languages/erlang-repl.min", documenterjs)
            @test occursin("languages/erlang.min", documenterjs)
        end
    end

    @testset "HTML: local" begin
        doc = Main.examples_html_local_doc

        @test isa(doc, Documenter.Documents.Document)

        let build_dir = joinpath(examples_root, "builds", "html-local")

            index_html = read(joinpath(build_dir, "index.html"), String)
            @test occursin("bold output from MarkdownOnly", index_html)

            @test isfile(joinpath(build_dir, "index.html"))
            @test isfile(joinpath(build_dir, "omitted.html"))
            @test isfile(joinpath(build_dir, "hidden.html"))
            @test isfile(joinpath(build_dir, "lib", "autodocs.html"))

            # Assets
            @test joinpath(build_dir, "assets", "documenter.js") |> isfile
            documenterjs = String(read(joinpath(build_dir, "assets", "documenter.js")))
            @test occursin("languages/julia.min", documenterjs)
            @test occursin("languages/julia-repl.min", documenterjs)
        end
    end

    @testset "Markdown" begin
        doc = Main.examples_markdown_doc

        @test isa(doc, Documenter.Documents.Document)

        let build_dir  = joinpath(examples_root, "builds", "markdown"),
            source_dir = joinpath(examples_root, "src")

            @test isdir(build_dir)
            @test isdir(joinpath(build_dir, "assets"))
            @test isdir(joinpath(build_dir, "lib"))
            @test isdir(joinpath(build_dir, "man"))

            @test isfile(joinpath(build_dir, "index.md"))
            @test isfile(joinpath(build_dir, "assets", "mathjaxhelper.js"))
            @test isfile(joinpath(build_dir, "assets", "Documenter.css"))
            @test isfile(joinpath(build_dir, "assets", "custom.css"))
            @test isfile(joinpath(build_dir, "assets", "custom.js"))
            @test isfile(joinpath(build_dir, "lib", "functions.md"))
            @test isfile(joinpath(build_dir, "man", "tutorial.md"))
            @test isfile(joinpath(build_dir, "man", "data.csv"))
            @test isfile(joinpath(build_dir, "man", "julia.svg"))

            @test (==)(
                read(joinpath(source_dir, "man", "data.csv"), String),
                read(joinpath(build_dir,  "man", "data.csv"), String),
            )
        end

        @test doc.user.root   == examples_root
        @test doc.user.source == "src"
        @test doc.user.build  == "builds/markdown"
        @test doc.user.clean  == true
        @test doc.user.format == [Documenter.Writers.MarkdownWriter.Markdown()]

        @test realpath(doc.internal.assets) == realpath(joinpath(dirname(@__FILE__), "..", "..", "assets"))

        @test length(doc.blueprint.pages) == 18

        let headers = doc.internal.headers
            @test Documenter.Anchors.exists(headers, "Documentation")
            @test Documenter.Anchors.exists(headers, "Documentation")
            @test Documenter.Anchors.exists(headers, "Index-Page")
            @test Documenter.Anchors.exists(headers, "Functions-Contents")
            @test Documenter.Anchors.exists(headers, "Tutorial-Contents")
            @test Documenter.Anchors.exists(headers, "Index")
            @test Documenter.Anchors.exists(headers, "Tutorial")
            @test Documenter.Anchors.exists(headers, "Function-Index")
            @test Documenter.Anchors.exists(headers, "Functions")
            @test Documenter.Anchors.isunique(headers, "Functions")
            @test Documenter.Anchors.isunique(headers, "Functions", joinpath("builds", "markdown", "lib", "functions.md"))
            let name = "Foo", path = joinpath("builds", "markdown", "lib", "functions.md")
                @test Documenter.Anchors.exists(headers, name, path)
                @test !Documenter.Anchors.isunique(headers, name)
                @test !Documenter.Anchors.isunique(headers, name, path)
                @test length(headers.map[name][path]) == 4
            end
        end

        @test length(doc.internal.objects) == 41
    end
end


text:
using ROOTDataFrames, DataFrames
using Base.Test

df = DataFrame(a=Float64[1.0, 2.0, 3.0], b=Int64[0, 1, 2], e=Int32[1, 1, 1])

writetree("test.root", df)
@test isfile("test.root")
print_with_color(:green, "writetree OK\n")

tdf = TreeDataFrame("test.root")
print_with_color(:green, "TreeDataFrame from file OK\n")
@test size(tdf) == size(df)
print_with_color(:green, "size OK\n")

#@test all(names(tdf) .== [:a, :b, :e])

@test all(tdf[:a] .== df[:a])
print_with_color(:green, "tdf[:a] .== df[:a] OK\n")

@test all(tdf[[:a, :b]][:a] .== df[[:a, :b]][:a])
print_with_color(:green, "tdf[[:a, :b]][:a] .== df[[:a, :b]][:a] OK\n")
@test all(tdf[[:a, :b]][:b] .== df[[:a, :b]][:b])

@test all(tdf[:b] .== df[:b])
@test all(tdf[:e][1:2] .== df[:e][1:2])

@test all(tdf[[true, true, false], [:a]][:a] .== [1.0, 2.0])

print_with_color(:green, "masked index OK\n")

tdf = TreeDataFrame(["test.root", "test.root"])
print_with_color(:green, "two-file chain opened OK\n")
@test nrow(tdf) == 6
print_with_color(:green, "two-file chain size OK\n")
@test all(tdf[:a] .== vcat(df[:a], df[:a]))
print_with_color(:green, "two-file :a OK\n")

#rm("test.root")

println("bigdf")
df = DataFrame(x = Int64[i for i=1:10^7])
writetree("big.root", df)
print_with_color(:green, "bigdf written OK\n")

tdf = TreeDataFrame("big.root")
@test nrow(tdf) == nrow(df)
print_with_color(:green, "bigdf opened OK\n")

tdf = TreeDataFrame(["big.root", "big.root"])
@test nrow(tdf) == 2 * nrow(df)
print_with_color(:green, "2x bigdf opened OK\n")

#rm("big.root")


text:
# Not tests as such;
# just some code to track down weird shit ftom the 1st 2
# simulations front ends
#
using Test

using UUIDs
using Observables
using CSV

using ScottishTaxBenefitModel
using .BCCalcs
using .Definitions
using .ExampleHelpers
using .FRSHouseholdGetter
using .GeneralTaxComponents
using .ModelHousehold
using .Monitor
using .Results
using .Runner
using .RunSettings
using .SimplePovertyCounts: GroupPoverty
using .SingleHouseholdCalculations
using .STBIncomes
using .STBOutput
using .STBParameters
using .Utils

const BASE_UUID = UUID("985c312f-129b-4acd-9e40-cb629d184183")

function load_system()::TaxBenefitSystem
	sys = load_file( joinpath( Definitions.MODEL_PARAMS_DIR, "sys_2021_22.jl" ))
	#
	# Note that as of Budget21 removing these doesn't actually happen till May 2022.
	#
	load_file!( sys, joinpath( Definitions.MODEL_PARAMS_DIR, "sys_2021-uplift-removed.jl"))
	# uc taper to 55
	load_file!( sys, joinpath( Definitions.MODEL_PARAMS_DIR, "budget_2021_uc_changes.jl"))
	weeklyise!( sys )
	return sys
end

f = open("weird_results.md", "w");

function initialise_settings()::Settings
    settings = Settings()
	settings.uuid = BASE_UUID
	settings.means_tested_routing = modelled_phase_in
    settings.run_name="run-$(date_string())"
	settings.income_data_source = ds_frs
	settings.dump_frames = true
	settings.do_marginal_rates = true
	settings.requested_threads = 4
	settings.dump_frames = true
	return settings
end

const BASE_SETTINGS = initialise_settings()

function get_and_print(
	f :: IOStream,
	hid :: BigInt,
	datayear :: Int,
	settings :: Settings,
	params   :: Vector{TaxBenefitSystem{T}} ) where T
	num_systems = size(params)[1]
	hh = get_household( hid, datayear )
	println( f, ModelHousehold.to_string(hh) )
	for sysno in 1:num_systems
		res = do_one_calc( hh, params[sysno], settings )
		println( f, "### RESULTS FOR SYSTEM $sysno")
		println( f, Results.to_string( res ))
	end
end

#
# a) increasing scottish child payment - some METRs go 20-30 -> 10-20
#
@testset "scottish child payment metrs" begin
	settings = initialise_settings()
	settings.run_name = "scottish_child_payment"
	sys = load_system()
    chsys = deepcopy( sys )
    obs = Observable(
		Progress(settings.uuid, "",0,0,0,0))
	chsys.scottish_child_payment.amount = 20.0
	params = [sys, chsys]
    results = do_one_run( settings, params, obs )
	p = collect(keys(skipmissing( results.indiv[1].metr )))
	ip1 = results.indiv[1][p,:]
	ip2 = results.indiv[2][p,:]
	d = irdiff( ip1, ip2 )
	targets = d[(abs.(d.metr) .> 0.00001),[:hid,:pid,:data_year,:metr,:means_tested_benefits]]
	println( f, "# scottish child payment metrs\n")
	CSV.write( "scottish_child_payment_metrs.csv", targets )
	for t in eachrow( targets )
		get_and_print(
			f,
			t.hid,
			t.data_year,
			settings,
			params )

	end
	n = size( results.indiv[1] )[1]
	for i in 1:n
		r1 = results.indiv[1][i,:]
		r2 = results.indiv[2][i,:]
		println( r1.hid )
		if r1.metr === missing
			@test r2.metr === missing
		else
			@test r1.metr ≈ r2.metr "r1.metr ≈ r2.metr r1=$(r1.metr) ≈ r2=$(r2.metr) "
		end
	end
	m1 = STBOutput.metrs_to_hist( results.indiv[1] )
	m2 = STBOutput.metrs_to_hist( results.indiv[2] )
	println(m1.hist.weights)
	println(m2.hist.weights)
	@test m1.mean ≈ m2.mean
	@test m1.hist.weights ≈ m2.hist.weights

end

#
# b) 1p on income tax - total raised  £412m. > total income tax +405
#
@testset "1p income tax" begin
    sys = load_system()
	settings = initialise_settings()
	settings.run_name = "plus_1p_income_tax"
    chsys = deepcopy( sys )
    obs = Observable(
		Progress(settings.uuid, "",0,0,0,0))
	chsys.it.non_savings_rates[1:3] .+= 0.01
	params = [sys, chsys]
    results = do_one_run( settings, params, obs )
	d = idiff( results.income[1], results.income[2] )
	targets = d[(d.income_tax.>0),[:hid,:pid,:data_year,:income_tax]]
	println( f, "# +1p income tax\n")
	CSV.write( "1p_income_tax.csv", targets )
	i = 0
	for t in eachrow( targets )
		i += 1
		if i < 10
			get_and_print(
				f,
				t.hid,
				t.data_year,
				settings,
				params )
		end
	end
end

#
# c) cutting basic 20-19 : a few losers
#
# Gainers	4,158	0.08
# Losers	718	0.01 <- pension contributions, no tax payable
# Unchanged	5,461,124	99.91
@testset "-1p income tax" begin
	settings = initialise_settings()
	settings.run_name = "minus_1p_income_tax"
	sys = load_system()
    chsys = deepcopy( sys )
    obs = Observable(
		Progress(settings.uuid, "",0,0,0,0))
		chsys.it.non_savings_rates[1:3] .-= 0.01
	params = [sys, chsys]
    results = do_one_run( settings, params, obs )
	d = hdiff( results.hh[1], results.hh[2] )
	targets = d[(d.bhc_net_income.<0),[:hid,:data_year,:bhc_net_income]]
	println( f, "# -1p income tax\n")
	CSV.write( "m1p_income_tax.csv", targets )
	for t in eachrow( targets )
		get_and_print(
			f,
			t.hid,
			t.data_year,
			settings,
			params )

	end
end

close( f )

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
################################################################################
###  Introduction to Scientific Programming and Machine Learning with Julia  ###
###                                                                          ###
### Run each script on a new clean Julia session                             ###
### GitHub: https://github.com/sylvaticus/IntroSPMLJuliaCourse               ###
### Licence (apply to all material of the course: scripts, videos, quizes,..)###
### Creative Commons By Attribution (CC BY 4.0), Antonello Lobianco          ###
################################################################################

# # 0302 - The Perceptron algorithm for linear classification

# ## Some stuff to set-up the environment..

cd(@__DIR__)
using Pkg
Pkg.activate(".")
## If using a Julia version different than 1.7 please uncomment and run the following line (the guarantee of reproducibility will however be lost)
## Pkg.resolve()
Pkg.instantiate()
using Random
Random.seed!(123)

# ## Perceptron elementary operations

using StatsPlots
function plot2DClassifierWithData(X,y,θ;d1=1,d2=2,origin=false,pid=1)
    colors = [y == -1 ? "red" : "green" for y in y]
    labels = [y == -1 ? "-1" : "+1" for y in y]
    minD1,maxD1 = extrema(X[:,d1])
    minD2,maxD2 = extrema(X[:,d2])
    myplot = scatter(X[:,d1],X[:,d2], colour=colors, title="Linear classifier in 2D",xlabel="Dimx: $d1", ylabel="Dimy: $d2", group=labels)
    constTerm = 0.0
    if !origin
        d1 += 1
        d2 += 1
        constTerm = -θ[1]/θ[d2]
    end
    d2Class(x) = constTerm -x * θ[d1]/θ[d2]
    if θ[d2] == 0
        vline!([0], color= "blue",label="",linewidth=5)
    else
        plot!(d2Class,min(θ[d1],minD1),max(maxD1,θ[d1]), color= "blue",label="",linewidth=5)
    end
    plot!([0,θ[d1]],[0,θ[d2]],arrow=true,color=:black,linewidth=2,label="")
    display(myplot)
    savefig("currentPlot$(pid).svg");
end
isClassificationError(θ,y,x) =  y * (θ' * x) <= eps()
perceptronUpdate(θ,y,x)      = return θ .+ y .* x

X = [ 2 4
     -6 1]
y = [-1,-1]
θ₀ = [0,0]
θ = θ₀

ϵ = isClassificationError(θ,y[1],X[1,:])
θ = perceptronUpdate(θ,y[1],X[1,:])
plot2DClassifierWithData(X,y,θ,origin=true,pid=1)
# ![](currentPlot1.svg)
ϵ = isClassificationError(θ,y[1],X[1,:])
ϵ = isClassificationError(θ,y[2],X[2,:])
θ = perceptronUpdate(θ,y[2],X[2,:])
plot2DClassifierWithData(X,y,θ,origin=true,pid=2)
# ![](currentPlot2.svg)
ϵ = isClassificationError(θ,y[2],X[2,:])
ϵ = isClassificationError(θ,y[1],X[1,:])

X = [ 2 4
     1 -2]
θ = θ₀

ϵ = isClassificationError(θ,y[1],X[1,:])
θ = perceptronUpdate(θ,y[1],X[1,:])
plot2DClassifierWithData(X,y,θ, origin=true,pid=3)
# ![](currentPlot3.svg)
ϵ = isClassificationError(θ,y[1],X[1,:])
ϵ = isClassificationError(θ,y[2],X[2,:])
θ = perceptronUpdate(θ,y[2],X[2,:])
plot2DClassifierWithData(X,y,θ, origin=true,pid=4)
# ![](currentPlot4.svg)
ϵ = isClassificationError(θ,y[1],X[1,:])
ϵ = isClassificationError(θ,y[2],X[2,:])
θ = perceptronUpdate(θ,y[2],X[2,:])
plot2DClassifierWithData(X,y,θ,origin=true,pid=5)
# ![](currentPlot5.svg)
ϵ = isClassificationError(θ,y[1],X[1,:])
ϵ = isClassificationError(θ,y[2],X[2,:])
θ

X = [ 2 4
     -2 2]
y = [-1,1]
θ = θ₀
ϵ = isClassificationError(θ,y[1],X[1,:])
θ = perceptronUpdate(θ,y[1],X[1,:])
plot2DClassifierWithData(X,y,θ, origin=true,pid=6)
# ![](currentPlot6.svg)
ϵ = isClassificationError(θ,y[2],X[2,:])
θ = perceptronUpdate(θ,y[2],X[2,:])
plot2DClassifierWithData(X,y,θ,origin=true,pid=7)
# ![](currentPlot7.svg)

# ## The complete algorithm

function perceptronOrigin(X,y,epochs=1;verbose=false)
    (nR,nD) = size(X)
    local θ = zeros(nD)
    for t in 1:epochs
        for n in 1:nR
            if verbose
                println("$n: X[n,:] \t θ: $θ")
            end
            if isClassificationError(θ,y[n],X[n,:])
                θ = perceptronUpdate(θ,y[n],X[n,:])
                if verbose
                    println("**update! New theta: $θ")
                end
            end
        end
    end
    return θ
end
θopt =  perceptronOrigin(X,y,verbose=true)
plot2DClassifierWithData(X,y,θopt, origin=true,pid=8)
# ![](currentPlot8.svg)

using BetaML, DelimitedFiles
baseDir          = joinpath(dirname(pathof(BetaML)),"..","test","data")
perceptronData   = readdlm(joinpath(dirname(pathof(BetaML)),"..","test","data","binary2DData.csv"),'\t')

nR = size(perceptronData,1)

idx = shuffle(1:nR)
perceptronData = perceptronData[idx,:]
X                = copy(perceptronData[:,[2,3]])
y                = convert(Array{Int64,1},copy(perceptronData[:,1]))
θopt             = perceptronOrigin(X,y,verbose=true)

plot2DClassifierWithData(X,y,θopt, origin=true,pid=20)
# ![](currentPlot20.svg)

# ## A better organisation

# Now we rewrite the perceptron algorithm setting all the parameters in a structure and using what could be a generic interface for any supervised model. This is the approach used by most ML libraries.
# We will see how to measure the classification error and as we are here we add the constant term with the constant addition to the data trick (note that editing every time the feature matrix is _NOT_ efficient and we do it here only for simplicity. A better way is to explicitly model the perceptron model with a constant parameter.)

abstract type SupervisedModel end
abstract type TrainingOptions end

mutable struct Perceptron <: SupervisedModel
    θ::Vector{Float64}
end

mutable struct PerceptronTrainingOptions <: TrainingOptions
    epochs::Int64
    verbose::Bool
    shuffle::Bool
    function PerceptronTrainingOptions(;epochs=1,verbose=false,shuffle=false)
        return new(epochs,verbose,shuffle)
    end
end

function predict(model::Perceptron,x::AbstractVector)
    x = vcat(1.0,x)
    x' * model.θ > eps() ? (return 1) : (return -1)
end

function predict(model::Perceptron,X::AbstractMatrix)
    return [predict(model,r) for r in eachrow(X)]
end

function update!(model::Perceptron,X::Vector,y)
    X       = vcat(1.0,X)
    model.θ = model.θ .+ y .* X
    return model.θ
end

function train!(model::Perceptron,X,y,ops=PerceptronTrainingOptions()::TrainingOptions)
    epochs  = ops.epochs
    verbose = ops.verbose
    (nR,nD) = size(X)
    nD += 1
    for t in 1:epochs
        errors = 0
        if ops.shuffle   # more efficient !
          idx = shuffle(1:nR)
          X = X[idx,:]
          y = y[idx]
        end
        for n in 1:nR
            if verbose
                println("$n: X[n,:] \t θ: $(model.θ)")
            end
            if  predict(model,X[n,:]) != y[n]
                errors += 1
                θ = update!(model,X[n,:],y[n])
                if verbose
                    println("**update! New theta: $(model.θ)")
                end
            end
        end
        if verbose
            println("Epoch $t errors: $errors")
        end
    end
    return model.θ
end

# ## Testing the Perceptron algorithm

m   = Perceptron(zeros(size(X,2)+1))
ops = PerceptronTrainingOptions()
train!(m,X,y,ops)
plot2DClassifierWithData(X,y,m.θ,pid=9)
# ![](currentPlot9.svg)

ŷ = predict(m,X)
inSampleAccuracy = sum(y .== ŷ)/length(y)

# Let's see if shuffling and increasing epochs we improve the accuracy....
ops = PerceptronTrainingOptions(verbose=false,epochs=5,shuffle=true)
m   = Perceptron(zeros(size(X,2)+1))
train!(m,X,y,ops)
plot2DClassifierWithData(X,y,m.θ,pid=10)
# ![](currentPlot10.svg)
ŷ = predict(m,X)
inSampleAccuracy = sum(y .== ŷ)/length(y)

# ## Cross-validation and hyperparameters optimisation

# Let's see now using separate training/validation
# We use the BetaML `partition()` function
((xtrain,xtest),(ytrain,ytest)) = partition([X,y],[0.6,0.4])
m             = Perceptron(zeros(size(X,2)+1))
ops           = PerceptronTrainingOptions(epochs=5,shuffle=true)
train!(m,xtrain,ytrain,ops)
plot2DClassifierWithData(xtrain,ytrain,m.θ,pid=11)
# ![](currentPlot11.svg)
ŷtrain = predict(m,xtrain)
trainAccuracy = accuracy(ŷtrain,ytrain)
sum(ytrain  .== ŷtrain)/length(ytrain)
## @edit accuracy(ŷtrain,ytrain)
ŷtest         = predict(m,xtest)
testAccuracy  = accuracy(ŷtest,ytest)
plot2DClassifierWithData(xtest,ytest,m.θ,pid=12)
# ![](currentPlot12.svg)
cfOut = ConfusionMatrix(ŷtest,ytest)
print(cfOut)

# Lets use CrossValidation
((xtrain,xvalidation,xtest),(ytrain,yvalidation,ytest)) = partition([X,y],[0.6,0.2,0.2])
## Very few records..... let's go back to using only two subsets but with CrossValidation
((xtrain,xtest),(ytrain,ytest)) = partition([X,y],[0.6,0.4])

sampler    = KFold(nSplits=10)

ops     = PerceptronTrainingOptions(epochs=10,shuffle=true)
(acc,σ) = crossValidation([xtrain,ytrain],sampler) do trainData,valData,rng
                (xtrain,ytrain) = trainData; (xval,yval) = valData
                m               = Perceptron(zeros(size(xtrain,2)+1))
                train!(m,xtrain,ytrain,ops)
                ŷval         = predict(m,xval)
                valAccuracy  = accuracy(ŷval,yval)
                return valAccuracy
            end

epochsSet  = 1:10:301
shuffleSet = [false,true]

bestE       = 0
bestShuffle = false
bestAcc     = 0.0

for e in epochsSet, s in shuffleSet
    global bestE, bestShuffle, bestAcc
    local acc
    local ops  = PerceptronTrainingOptions(epochs=e,shuffle=s)
    (acc,_)    = crossValidation([xtrain,ytrain],sampler) do trainData,valData,rng
                    (xtrain,ytrain) = trainData; (xval,yval) = valData
                    m               = Perceptron(zeros(size(xtrain,2)+1))
                    train!(m,xtrain,ytrain,ops)
                    ŷval            = predict(m,xval)
                    valAccuracy     = accuracy(ŷval,yval)
                    return valAccuracy
                end
    if acc > bestAcc
        bestAcc     = acc
        bestE       = e
        bestShuffle = s
    end
end

bestAcc
bestE
bestShuffle

ops = PerceptronTrainingOptions(epochs=bestE,shuffle=bestShuffle)
m   = Perceptron(zeros(size(xtest,2)+1))
train!(m,xtrain,ytrain,ops)
ŷtest           = predict(m,xtest)
testAccuracy    = accuracy(ŷtest,ytest)

plot2DClassifierWithData(xtest,ytest,m.θ,pid=13)
# ![](currentPlot13.svg)

max_stars_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_issues_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_forks_repo_name: sylvaticus/IntroSPMLJuliaCourse


text:
################################################################################
# Runs all .jl files in the directory of this script containing a @test or
# @testset macro, except those specified to be excluded in run_all_tests.jl
################################################################################

# This section activates the test project. It's needed for "Pkg.test()" to work for the root project.
push!(LOAD_PATH, "@stdlib") # needed for "using Pkg"
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate();

include("run_all_tests.jl")

test_name_regex = length(ARGS) < 1 ? "" : ARGS[1]
run_all_tests(test_name_regex)

max_issues_repo_name: maxbennedich/julia-regression-analysis

max_forks_repo_name: maxbennedich/julia-regression-analysis


text:
## Packages
using BenchmarkTools
using LoopVectorization

## Reference kernel
function divergence_ref(
    u, v, w,
    is, ie, js, je, ks, ke)

    divmax = 0
    @inbounds for k in ks:ke
        for j in js:je
            for i in is:ie
                div = u[i, j, k] + v[i, j, k] + w[i, j, k]
                divmax = max(abs(div), divmax)
            end
        end
    end

    return divmax
end

## Reference kernel
function divergence(
    u, v, w,
    is, ie, js, je, ks, ke)

    divmax = 0
    @tturbo for k in ks:ke
        for j in js:je
            for i in is:ie
                div = u[i, j, k] + v[i, j, k] + w[i, j, k]
                divmax = max(abs(div), divmax)
            end
        end
    end

    return divmax
end

## User input
itot = 384; jtot = 384; ktot = 384;
igc = 1; jgc = 1; kgc = 1;

icells = itot + 2igc; jcells = jtot + 2jgc; kcells = ktot + 2kgc
is = igc+1; ie = igc+itot; js = jgc+1; je = jgc+jtot; ks = kgc+1; ke = kgc+ktot

u = rand(icells, jcells, kcells)
v = rand(icells, jcells, kcells)
w = rand(icells, jcells, kcells)

## Run reference kernel
div_ref = divergence_ref(
    u, v, w,
    is, ie, js, je, ks, ke)

@btime dummy = divergence_ref(
    $u, $v, $w,
    $is, $ie, $js, $je, $ks, $ke)

## Run fast kernel
div = divergence(
    u, v, w,
    is, ie, js, je, ks, ke)

@btime divergence(
    $u, $v, $w,
    $is, $ie, $js, $je, $ks, $ke)

## Compare results
println(div_ref, " ", div, " ", isapprox(div, div_ref))


text:
function Legendre_polynomials!(x::Float64, lmax::Int64; Pₗ = nothing)
    #### returns an array of size lmax_half. Each cell contains one (even) Legendre polynomials evaluated at point x ####
    # getting all the constant we will need
    lmax_half = div(lmax, 2)

    # Preallocating Pₗ if it hasn't been given
    if isnothing(Pₗ)
        Pₗ = zeros(lmax_half) # will store all the Legendre polynomials
    end

    # Let's start the recurrence relation -- We use Bonnet's recurrence formula
    P_l = nothing # these 3 lines just to help the interpreter with scopes
    P_lmin = nothing
    P_lminmin = nothing

    for l = 0:lmax
        if (l == 0)
            P_l = 1
        elseif (l == 1)
            P_l = x
        else
            P_l = x * P_lmin * (2 * l - 1) / l - P_lminmin * (l - 1) / l # l Pₗ = x * (2 l - 1) Pₗ₋₁ + (l - 1) Pₗ₋₂
        end

        if ((l % 2 == 0) && (l > 0)) # we only store the even harmonics, and l = 0 is useless because it brings a constant to the energy
            l_half = div(l, 2)
            Pₗ[l_half] = P_l
        end

        P_lminmin = P_lmin
        P_lmin = P_l
    end
end

function Normalized_associated_Legendre_polynomials!(x::Float64, lmax::Int64; P = nothing)
    #### Returns an array of size lmax_half * (lmax + 1). Each cell contains one normalized associated Legendre polynomails Pᴺₗₘ(x)        ####
    #### The definition is Pᴺₗₘ = √((2l + 1) / (4 π)) × √(l - m)! / (l + m)! × Pₗₘ where Pₗₘ is the standard associated Legendre polynomial. ####
    # getting all the constants we will need
    lmax_half = div(lmax, 2)

    # Preallocating P if it hasn't been given
    if isnothing(P)
        P = zeros(lmax + 1, lmax_half) # Will store the normalized associated Legendre polynomials
    end

    # Let's start the recurrence -- We use the method from Numerical Recipes (eq. 6.7.9)
    Pₘₘ = nothing # this is just to help the compilator with scopes

    for m = 0:lmax
        # First let's get Pₘₘ by a first recurrence relation
        if (m == 0)
            Pₘₘ = sqrt(1 / (4 * π)) # P₀₀ = sqrt(1 / 4π)
        else
            Pₘₘ = -1.0 * sqrt((2 * m + 1) / (2 * m)) * sqrt(1 - x^(2)) * Pₘₘ # Pₘ,ₘ = - sqrt((2m + 1) / 2m) * sqrt(1 - x²) * Pₘ₋₁,ₘ₋₁
        end

        # The we go from m,m to m,l step by step
        Pₗₘ = nothing # this is just to help the compilator with scopes
        Pₗₘ_min = nothing
        Pₗₘ_minmin = nothing

        for l = m:lmax
            if (l == m)
                Pₗₘ = Pₘₘ
            elseif (l == (m + 1))
                Pₗₘ = x * sqrt(2 * m + 3) * Pₘₘ # Pₘ,ₘ₊₁ = x sqrt(2m + 3) Pₘₘ
            else
                Pₗₘ = sqrt((4 * l^2 - 1) / (l^2 - m^2)) * ( x * Pₗₘ_min - sqrt( ((l - 1)^2 - m^2) / (4 * (l - 1)^2 -1)) * Pₗₘ_minmin )
                # Pₘₗ = sqrt((4l² - 1) / (l² - m²)) * [x Pₘ,ₗ₋₁ - sqrt(((l - 1)² - m²) / (4 (l - 1)² - 1)) Pₘₗ₋₂]
            end

            if (((l % 2) == 0) && (l > 0)) # We store the coefficients only if l is even
                l_half = div(l, 2)
                P[m + 1, l_half] = Pₗₘ
            end

            Pₗₘ_minmin = Pₗₘ_min
            Pₗₘ_min = Pₗₘ
        end
    end

    return(P)
end

max_stars_repo_name: NathanMagnan/VRROOMe

max_issues_repo_name: NathanMagnan/VRROOMe

max_forks_repo_name: NathanMagnan/VRROOMe