text:
function generateMnliFast(Pk::Function, BiasDict22::Dict{String,F}, f::Real) where {F<:Real}
    M = OffsetArray{Array{Float64}}(undef, 0:4, 0:8, 0:8)

    b1 = BiasDict22["b1"]
    bη = BiasDict22["bη"]
    b2 = BiasDict22["b2"]
    bK² = BiasDict22["bK2"]
    bδη = BiasDict22["bδη"]
    bη² = BiasDict22["bη2"]
    bKKpara = BiasDict22["bKK∥"]
    bΠ2para = BiasDict22["bΠ2∥"]
    q = "best"
    r, xi20 = ξ(Pk,2,0,q=q)
    _, xi00 = ξ(Pk,0,0,q=q)
    _, xi1m1 = ξ(Pk,1,-1,q=q)
    _, xi3m1 = ξ(Pk,3,-1,q=q)
    _, xi2m2 = ξ(Pk,2,-2,q=q)
    _, xi4m2 = ξ(Pk,4,-2,q=q)
    _, xi10 = ξ(Pk,1,0,q=q)
    _, xi30 = ξ(Pk,3,0,q=q)
    _, xi4m1 = ξ(Pk,4,-1,q=q)
    _, xi2m1 = ξ(Pk,2,-1,q=q)
    _, xi0m1 = ξ(Pk,0,-1,q=q)
    _, xi3m2 = ξ(Pk,3,-2,q=q)
    _, xi1m2 = ξ(Pk,1,-2,q=q)
    _, xi0m2 = ξ(Pk,0,-2,q=q)
    _, xi40 = ξ(Pk,4,0,q=q)

    xi00sq = xi00.^2
    xi20sq = xi20.^2
    xi40sq = xi40.^2
    xi0040 = xi00.*xi40
    xi2040 = xi20.*xi40
    xi104m1 = xi10.*xi4m1
    xi304m1 = xi30.*xi4m1
    xi0m200 = xi0m2.*xi00
    xi2m220 = xi2m2.*xi20
    xi1m1sq = xi1m1.^2
    xi3m1sq = xi3m1.^2
    xi002m2 = xi00.*xi2m2
    xi1m13m1 = xi1m1.*xi3m1
    xi004m2 = xi00.*xi4m2
    xi204m2 = xi20.*xi4m2
    xi0m11m2 = xi0m1.*xi1m2
    xi1m22m1 = xi1m2.*xi2m1
    xi2m13m2 = xi2m1.*xi3m2
    xi0m13m2 = xi0m1.*xi3m2
    xi0m2sq = xi0m2.^2
    xi2m2sq = xi2m2.^2
    xi0m22m2 = xi0m2.*xi2m2

    M[0,0,0] = @. ((35280*bK²^2 + 6664*bKKpara^2 - 112*bKKpara*(150*b1 - 2*f*(15 + 28*f)*bη - 56*f^2*bη² - 69*bΠ2para) -
            672*bK²*(75*b1 - 35*bKKpara - f*(15 + 14*f)*bη - 14*f^2*bη² - 10*bΠ2para) + 3*(36000*b1^2 + 63000*b1*b2 +
                33075*b2^2 - 42000*b1^2*f - 44100*b1*b2*f + 20580*b1^2*f^2 + 20580*f^2*bδη^2 + 12*f^2*(216 - 196*f +
                    441*f^2)*bη^2 + 11760*b1*f^2*bη² + 14700*b2*f^2*bη² - 17640*b1*f^3*bη² + 5292*f^4*bη²^2 +
                32400*b1*bΠ2para + 35700*b2*bΠ2para - 31640*b1*f*bΠ2para + 14504*f^2*bη²*bΠ2para + 16468*bΠ2para^2 -
                140*f*bδη*(300*b1 + 315*b2 - 294*b1*f + 6*f*(-10 + 21*f)*bη + 126*f^2*bη² + 226*bΠ2para) +
                4*f*bη*(525*b2*(-6 + 7*f) - 90*b1*(40 - 56*f + 49*f^2) + 294*f^2*(-2 + 9*f)*bη² + (-2580 +
                        3626*f)*bΠ2para)))*xi00sq)/396900. + (2*(8820*bK²^2 + 833*bKKpara^2 - 7*bKKpara*(-1500*b1 -
                2205*b2 + 2058*b1*f + 2058*f*bδη + 10*(30 - 97*f)*f*bη - 970*f^2*bη² - 1620*bΠ2para) - 21*bK²*(-1500*b1 -
                2205*b2 + 1470*b1*f - 280*bKKpara + 1470*f*bδη + 300*f*bη - 602*f^2*bη - 602*f^2*bη² - 1270*bΠ2para) +
            3*(-9000*b1^2 - 11025*b1*b2 + 7350*b1^2*f + 1029*b1^2*f^2 + 1029*f^2*bδη^2 + 6*f^2*(-108 + 35*f +
                    252*f^2)*bη^2 - 1050*b1*f^2*bη² + 2058*b2*f^2*bη² - 3528*b1*f^3*bη² + 1512*f^4*bη²^2 -
                4950*b1*bΠ2para + 1470*b2*bΠ2para - 3724*b1*f*bΠ2para + 3808*f^2*bη²*bΠ2para + 1343*bΠ2para^2 -
                98*f*bδη*(-75*b1 - 21*b1*f + 3*f*(5 + 12*f)*bη + 36*f^2*bη² + 38*bΠ2para) + f*bη*(147*b2*(15 + 14*f) -
                    72*b1*(-50 + 35*f + 49*f^2) + 42*f^2*(5 + 72*f)*bη² + 2*(975 + 1904*f)*bΠ2para)))*xi20sq)/138915. +
    (8*(1125*b1^2 + 2205*bK²^2 + 294*bKKpara^2 - 450*b1*f*bη - 420*b1*f^2*bη + 81*f^2*bη^2 + 84*f^3*bη^2 + 56*f^4*bη^2 -
            420*b1*f^2*bη² + 84*f^3*bη*bη² + 112*f^4*bη*bη² + 56*f^4*bη²^2 - 14*bKKpara*(75*b1 - 3*f*(5 + 6*f)*bη -
                18*f^2*bη² - 17*bΠ2para) - 42*bK²*(75*b1 - 35*bKKpara - f*(15 + 14*f)*bη - 14*f^2*bη² - 10*bΠ2para) -
            300*b1*bΠ2para - 60*f*bη*bΠ2para + 112*f^2*bη*bΠ2para + 112*f^2*bη²*bΠ2para + 169*bΠ2para^2)*xi40sq)/77175.

    M[0,0,2] = @. (32*(3*f*bη - 5*bΠ2para)*(147*b1*f + 182*bKKpara + 147*f*bδη + 6*f^2*bη + 6*f^2*bη² +
            84*bΠ2para)*xi20sq)/972405. - (32*(3*f*bη - 5*bΠ2para)*(7*bKKpara + 4*f^2*bη + 4*f^2*bη² +
            7*bΠ2para)*xi40sq)/108045. + xi20*((32*(7*bKKpara + 3*(14*b1*f + 14*f*bδη - 8*f^2*bη - 8*f^2*bη² -
                    7*bΠ2para))*(3*f*bη - 5*bΠ2para)*xi00)/99225. + (32*(28*bKKpara - 3*(49*b1*f + 49*f*bδη - 38*f^2*bη -
                    38*f^2*bη² - 42*bΠ2para))*(3*f*bη - 5*bΠ2para)*xi40)/540225.)

    M[0,2,0] = @. (4*(2940*bK² + 980*bKKpara + 3*(-4200*b1 - 3675*b2 + 2450*b1*f + 2450*f*bδη + 2*(540 - 343*f)*f*bη -
                686*f^2*bη² - 2290*bΠ2para))*(3*f*bη - 5*bΠ2para)*xi00sq)/694575. - (4*(2940*bK² + 980*bKKpara +
            3*(-1680*b1 - 1029*b2 + 686*b1*f + 686*f*bδη + 2*(216 - 49*f)*f*bη - 98*f^2*bη² - 622*bΠ2para))*(3*f*bη -
            5*bΠ2para)*xi20sq)/194481. - (32*(525*b1 - 735*bK² - 245*bKKpara - 135*f*bη - 98*f^2*bη - 98*f^2*bη² -
            20*bΠ2para)*(3*f*bη - 5*bΠ2para)*xi40sq)/540225.

    M[0,2,2] = @. (2*(98*bKKpara^2 + 7*bKKpara*(-21840*b1 - 21609*b2 + 20874*b1*f + 20874*f*bδη + 2*(3432 - 5089*f)*f*bη -
                10178*f^2*bη² - 20526*bΠ2para) + 588*bK²*(294*b1*f + 21*bKKpara + 294*f*bδη - 184*f^2*bη - 184*f^2*bη² -
                175*bΠ2para) - 3*(41160*b1^2*f + 14406*b1^2*f^2 + 14406*f^2*bδη^2 + 48*f^3*(-11 + 441*f)*bη^2 +
                1680*b1*f^2*bη² + 28812*b2*f^2*bη² - 49392*b1*f^3*bη² + 21168*f^4*bη²^2 + 23520*b1*bΠ2para +
                50421*b2*bΠ2para - 46354*b1*f*bΠ2para + 51154*f^2*bη²*bΠ2para + 39956*bΠ2para^2 - 98*f*bδη*(-420*b1 -
                    294*b1*f + 12*f*(11 + 42*f)*bη + 504*f^2*bη² + 473*bΠ2para) + 2*f*bη*(-42*f*(134*b1 - 343*b2 +
                        588*b1*f) + 24*f^2*(-11 + 882*f)*bη² + (-3696 + 25577*f)*bΠ2para)))*xi20sq)/1.361367e6 -
    (16*(3234*bKKpara^2 + 4*f^3*(363 + 196*f)*bη^2 - 4620*b1*f^2*bη² + 6468*f^2*bK²*bη² + 784*f^4*bη²^2 -
            77*bKKpara*(105*b1 - 147*bK² - 3*f*(11 + 14*f)*bη - 42*f^2*bη² - 29*bΠ2para) - 8085*b1*bΠ2para +
            11319*bK²*bΠ2para + 814*f^2*bη²*bΠ2para - 1001*bΠ2para^2 + f*bη*(-4620*b1*f + 6468*f*bK² + 4*f^2*(363 +
                    392*f)*bη² + 2541*bΠ2para + 814*f*bΠ2para))*xi40sq)/1.663893e6 + xi20*((2*(2450*bKKpara^2 -
                28*bKKpara*(210*b1 + 294*b1*f + 294*f*bδη - 22*f*(3 + 14*f)*bη - 308*f^2*bη² - 261*bΠ2para) -
                588*bK²*(21*b1*f - 14*bKKpara + 21*f*bδη - 22*f^2*bη - 22*f^2*bη² - 28*bΠ2para) + 3*(-11760*b1^2*f -
                    15435*b1*b2*f + 11760*b1^2*f^2 + 11760*f^2*bδη^2 + 12*f^3*(-176 + 441*f)*bη^2 + 6720*b1*f^2*bη² +
                    10290*b2*f^2*bη² - 14994*b1*f^3*bη² + 5292*f^4*bη²^2 + 5880*b1*bΠ2para + 10290*b2*bΠ2para -
                    21154*b1*f*bΠ2para + 14104*f^2*bη²*bΠ2para + 8666*bΠ2para^2 - 7*f*bδη*(1680*b1 + 2205*b2 - 3360*b1*f +
                        6*f*(-88 + 357*f)*bη + 2142*f^2*bη² + 3022*bΠ2para) + 2*f*bη*(21*f*(248*b1 + 245*b2 - 357*b1*f) +
                        12*f^2*(-88 + 441*f)*bη² + (-924 + 7052*f)*bΠ2para)))*xi00)/138915. + (16*(2205*bKKpara^2 +
                147*bK²*(28*bKKpara - 3*(49*b1*f + 49*f*bδη - 38*f^2*bη - 38*f^2*bη² - 42*bΠ2para)) + bKKpara*(-2940*b1 -
                    9261*b1*f - 9261*f*bδη + 14*f*(66 + 623*f)*bη + 8722*f^2*bη² + 9044*bΠ2para) + 3*(5145*b1^2*f +
                    6*f^3*(209 + 196*f)*bη^2 - 3990*b1*f^2*bη² - 1372*b1*f^3*bη² + 1176*f^4*bη²^2 - 4410*b1*bΠ2para -
                    392*b1*f*bΠ2para + 1732*f^2*bη²*bΠ2para + 483*bΠ2para^2 - 49*f*bδη*(-105*b1 + f*(33 + 28*f)*bη +
                        28*f^2*bη² + 8*bΠ2para) + f*bη*(-7*b1*f*(801 + 196*f) + 6*f^2*(209 + 392*f)*bη² + 2*(693 +
                            866*f)*bΠ2para)))*xi40)/756315.)

    M[0,2,4] = @. (-128*f^2*(bη + bη²)*(3*f*bη - 5*bΠ2para)*xi20sq)/16807. + (256*f^2*(bη + bη²)*(3*f*bη -
            5*bΠ2para)*xi0040)/77175. - (512*f^2*(bη + bη²)*(3*f*bη - 5*bΠ2para)*xi2040)/151263. + (384*f^2*(bη +
            bη²)*(3*f*bη - 5*bΠ2para)*xi40sq)/420175.

    M[0,4,0] = @. (64*(3*f*bη - 5*bΠ2para)^2*xi00sq)/25725. - (128*(3*f*bη - 5*bΠ2para)^2*xi20sq)/36015. + (64*(3*f*bη -
            5*bΠ2para)^2*xi40sq)/60025.

    M[0,4,2] = @. (64*(3*f*bη - 5*bΠ2para)*(147*b1*f + 182*bKKpara + 147*f*bδη + 6*f^2*bη + 6*f^2*bη² +
            84*bΠ2para)*xi20sq)/756315. - (64*(3*f*bη - 5*bΠ2para)*(7*bKKpara + 4*f^2*bη + 4*f^2*bη² +
            7*bΠ2para)*xi40sq)/84035. + xi20*((64*(7*bKKpara + 3*(14*b1*f + 14*f*bδη - 8*f^2*bη - 8*f^2*bη² -
                    7*bΠ2para))*(3*f*bη - 5*bΠ2para)*xi00)/77175. + (64*(28*bKKpara - 3*(49*b1*f + 49*f*bδη - 38*f^2*bη -
                    38*f^2*bη² - 42*bΠ2para))*(3*f*bη - 5*bΠ2para)*xi40)/420175.)

    M[0,4,4] = @. (4*(18326*bKKpara^2 - 2156*bKKpara*(42*b1*f + 42*f*bδη - 40*f^2*bη - 40*f^2*bη² - 45*bΠ2para) +
            3*(79233*b1^2*f^2 + 79233*f^2*bδη^2 + 216*f^3*(-195 + 539*f)*bη^2 + 138600*b1*f^2*bη² + 158466*b2*f^2*bη² -
                271656*b1*f^3*bη² + 17248*f^2*bK²*bη² + 116424*f^4*bη²^2 - 135828*b1*f*bΠ2para + 244836*f^2*bη²*bΠ2para +
                61446*bΠ2para^2 - 22638*f*bδη*(-7*b1*f + 12*f^2*bη + 12*f^2*bη² + 6*bΠ2para) + 2*f^2*bη*(8624*bK² +
                    3*(77*(343*b2 + b1*(300 - 588*f)) + 36*f*(-195 + 1078*f)*bη² + 40806*bΠ2para))))*xi20sq)/8.319465e6 +
    (4*(2156*bKKpara^2 + 4312*bKKpara*(4*f^2*bη + 4*f^2*bη² + bΠ2para) + 3*(8085*b1^2*f^2 + 8085*f^2*bδη^2 + 12*f^3*(156 +
                    539*f)*bη^2 - 6160*b1*f^2*bη² - 10780*b1*f^3*bη² + 8624*f^2*bK²*bη² + 6468*f^4*bη²^2 -
                10780*b1*f*bΠ2para + 9816*f^2*bη²*bΠ2para + 4312*bΠ2para^2 - 5390*f*bδη*(-3*b1*f + 2*f^2*bη + 2*f^2*bη² +
                    2*bΠ2para) + 4*f^2*bη*(-1540*b1 - 2695*b1*f + 2156*bK² + 6*f*(78 + 539*f)*bη² +
                    2454*bΠ2para)))*xi0040)/848925. + (16*(1078*bKKpara^2 - 72*f^3*(39 + 98*f)*bη^2 - 7056*f^4*bη²^2 +
            6*f^2*bη²*(1540*b1 + 1372*b1*f - 2156*bK² + 1372*f*bδη - 2013*bΠ2para) + 49*bKKpara*(231*b1*f + 231*f*bδη -
                230*f^2*bη - 230*f^2*bη² - 110*bΠ2para) + 1617*(7*b1*f + 7*f*bδη - 4*bΠ2para)*bΠ2para -
            6*f^2*bη*(-1540*b1 - 1372*b1*f + 2156*bK² - 1372*f*bδη + 468*f*bη² + 2352*f^2*bη² +
                2013*bΠ2para))*xi2040)/1.663893e6 + (144*(7007*bKKpara^2 + 18*f^3*(169 + 196*f)*bη^2 -
            10010*b1*f^2*bη² + 14014*f^2*bK²*bη² + 3528*f^4*bη²^2 + 6396*f^2*bη²*bΠ2para + 7007*bΠ2para^2 +
            1274*bKKpara*(9*f^2*bη + 9*f^2*bη² + 11*bΠ2para) + 2*f^2*bη*(-5005*b1 + 7007*bK² + 9*f*(169 + 392*f)*bη² +
                3198*bΠ2para))*xi40sq)/6.0085025e7

    M[0,6,4] = @. (-320*f^2*(bη + bη²)*(3*f*bη - 5*bΠ2para)*xi20sq)/26411. + (128*f^2*(bη + bη²)*(3*f*bη -
            5*bΠ2para)*xi0040)/24255. - (1280*f^2*(bη + bη²)*(3*f*bη - 5*bΠ2para)*xi2040)/237699. + (192*f^2*(bη +
            bη²)*(3*f*bη - 5*bΠ2para)*xi40sq)/132055.

    M[0,6,6] = @. (32*f^2*(bη + bη²)*(4*bKKpara - 3*(7*b1*f + 7*f*bδη - 6*f^2*bη - 6*f^2*bη² -
                6*bΠ2para))*xi2040)/4851. - (32*f^2*(bη + bη²)*(15*bKKpara + 8*f^2*bη + 8*f^2*bη² +
            15*bΠ2para)*xi40sq)/24255.

    M[0,8,8] = @. (128*f^4*(bη + bη²)^2*xi40sq)/6435.

    M[1,0,1] = @. xi10*(-(f*(1680*f*bK²*(1 + bη) + 1120*f*bKKpara*(1 + bη) + 9*(126*f^2*(-2 + 5*f)*bη^2 -
                    350*f*bδη*(-7*b1 + 3*f + 3*f*bη) - 210*(5*b1 - 3*f)*f^2*bη² + f*bη*(1600*b1 + 875*b2 - 252*f -
                        2100*b1*f + 630*f^2 + 630*f^2*bη² + 1050*bΠ2para) + 25*(35*b2*f + 2*b1^2*(-50 + 49*f) -
                        7*b1*(15*b2 - 4*f + 6*f^2) - 6*(17*b1 - 7*f)*bΠ2para)))*xi0m1)/70875. + (4*f*(-210*bK²*(15*b1 -
                    7*f - 7*f*bη) + bKKpara*(-1785*b1 + 845*f + 845*f*bη) + 9*(18*f^2*(-1 + 10*f)*bη^2 - 5*f*bδη*(-35*b1 +
                        51*f + 51*f*bη) + f*bη*(-40*b1 + 175*b2 - 18*f - 465*b1*f + 180*f^2 + 180*f^2*bη² + 255*bΠ2para) +
                    5*(50*b1^2 + 10*b1*f + 35*b1^2*f + 35*b2*f - 51*b1*f^2 - 6*(7*b1 - 6*f)*f^2*bη² + (-33*b1 +
                            51*f)*bΠ2para)))*xi2m1)/70875.) - (4*f*xi30*(3*(875*b1^2 - 325*b1*f - 105*b1*f^2 -
                105*f^2*bδη - 640*b1*f*bη + 117*f^2*bη - 315*b1*f^2*bη + 180*f^3*bη - 105*f^2*bδη*bη + 117*f^2*bη^2 +
                180*f^3*bη^2 - 35*bK²*(35*b1 - 13*f - 13*f*bη) + bKKpara*(-490*b1 + 230*f + 230*f*bη) - 210*b1*f^2*bη² +
                180*f^3*bη² + 180*f^3*bη*bη² + 35*b1*bΠ2para + 105*f*bΠ2para + 105*f*bη*bΠ2para)*xi2m1 - 8*f*(1 +
                bη)*(-75*b1 + 105*bK² + 45*bKKpara + 27*f*bη + 20*f^2*bη + 20*f^2*bη²)*xi4m1))/55125.

    M[1,0,3] = @. xi30*((-16*f^2*(1 + bη)*(3*f*bη - 5*bΠ2para)*xi0m1)/18375. + (32*f^2*(1 + bη)*(3*f*bη -
                5*bΠ2para)*(4*xi2m1 - 3*xi4m1))/128625.) + (16*f^2*(1 + bη)*(3*f*bη - 5*bΠ2para)*xi10*(11*xi2m1 +
            4*xi4m1))/42875.

    M[1,2,1] = @. xi10*((-2*f*(2352*f*bK²*(1 + bη) + 1568*f*bKKpara*(1 + bη) + 9*(-3500*b1^2 - 3675*b1*b2 + 980*b1*f +
                    3430*b1^2*f + 1225*b2*f - 1470*b1*f^2 + 126*f^2*(-4 + 7*f)*bη^2 - 490*f*bδη*(-7*b1 + 3*f + 3*f*bη) -
                    1470*b1*f^2*bη² + 882*f^3*bη² - 4470*b1*bΠ2para + 1722*f*bΠ2para + f*bη*(2780*b1 + 1225*b2 - 504*f -
                        2940*b1*f + 882*f^2 + 882*f^2*bη² + 1722*bΠ2para)))*xi0m1)/99225. + (8*f*(-7*bKKpara*(357*b1 -
                    169*f - 169*f*bη) - 294*bK²*(15*b1 - 7*f - 7*f*bη) + 9*(350*b1^2 + 70*b1*f + 245*b1^2*f + 245*b2*f -
                    357*b1*f^2 + 36*f^2*(-1 + 7*f)*bη^2 - 7*f*bδη*(-35*b1 + 51*f + 51*f*bη) - 294*b1*f^2*bη² +
                    252*f^3*bη² - 141*b1*bΠ2para + 375*f*bΠ2para + f*bη*(-110*b1 + 245*b2 - 36*f - 651*b1*f + 252*f^2 +
                        252*f^2*bη² + 375*bΠ2para)))*xi2m1)/99225.) - (8*f*xi30*(3*(1225*b1^2 - 455*b1*f - 147*b1*f^2 -
                147*f^2*bδη - 1085*b1*f*bη + 234*f^2*bη - 441*b1*f^2*bη + 252*f^3*bη - 147*f^2*bδη*bη + 234*f^2*bη^2 +
                252*f^3*bη^2 - 49*bK²*(35*b1 - 13*f - 13*f*bη) + bKKpara*(-686*b1 + 322*f + 322*f*bη) - 294*b1*f^2*bη² +
                252*f^3*bη² + 252*f^3*bη*bη² + 364*b1*bΠ2para + 30*f*bΠ2para + 30*f*bη*bΠ2para)*xi2m1 - 8*f*(1 +
                bη)*(-105*b1 + 147*bK² + 63*bKKpara + 54*f*bη + 28*f^2*bη + 28*f^2*bη² - 27*bΠ2para)*xi4m1))/77175.

    M[1,2,3] = @. (-2*f*xi10*(3*(1980*b1*f + 2205*b1^2*f + 2205*b2*f - 3213*b1*f^2 + 2205*b1*f*bδη - 3213*f^2*bδη +
                1980*b1*f*bη + 2205*b2*f*bη - 924*f^2*bη - 5859*b1*f^2*bη + 2268*f^3*bη - 3213*f^2*bδη*bη - 924*f^2*bη^2 +
                2268*f^3*bη^2 + 168*f*bK²*(1 + bη) + bKKpara*(-294*b1 + 238*f + 238*f*bη) - 2646*b1*f^2*bη² +
                2268*f^3*bη² + 2268*f^3*bη*bη² - 1764*b1*bΠ2para + 3430*f*bΠ2para + 3430*f*bη*bΠ2para)*xi2m1 - 8*f*(1 +
                bη)*(378*bK² + 217*bKKpara - 3*(90*b1 + 70*b1*f + 70*f*bδη - 42*f*(1 + 2*f)*bη - 84*f^2*bη² -
                    49*bΠ2para))*xi4m1))/77175. + xi30*((-2*f*(-180*b1*f + 945*b1^2*f - 315*b1*f^2 + 945*b1*f*bδη -
                315*f^2*bδη - 180*b1*f*bη + 84*f^2*bη - 945*b1*f^2*bη + 378*f^3*bη - 315*f^2*bδη*bη + 84*f^2*bη^2 +
                378*f^3*bη^2 + 252*f*bK²*(1 + bη) + 168*f*bKKpara*(1 + bη) - 630*b1*f^2*bη² + 378*f^3*bη² +
                378*f^3*bη*bη² - 630*b1*bΠ2para + 238*f*bΠ2para + 238*f*bη*bΠ2para)*xi0m1)/11025. - (16*f*(11*(180*b1*f +
                    147*b1*f^2 + 147*f^2*bδη + 180*b1*f*bη - 84*f^2*bη + 441*b1*f^2*bη - 252*f^3*bη + 147*f^2*bδη*bη -
                    84*f^2*bη^2 - 252*f^3*bη^2 - 252*f*bK²*(1 + bη) + 7*bKKpara*(63*b1 - 31*f - 31*f*bη) + 294*b1*f^2*bη² -
                    252*f^3*bη² - 252*f^3*bη*bη² + 441*b1*bΠ2para - 175*f*bΠ2para - 175*f*bη*bΠ2para)*xi2m1 + 3*f*(1 +
                    bη)*(-495*b1 + 693*bK² + 462*bKKpara + 231*f*bη + 252*f^2*bη + 252*f^2*bη² +
                    77*bΠ2para)*xi4m1))/848925.)

    M[1,4,1] = @. xi10*((32*f*(-25*b1 + 7*f + 7*f*bη)*(3*f*bη - 5*bΠ2para)*xi0m1)/18375. - (64*f*(5*b1 + f +
                f*bη)*(3*f*bη - 5*bΠ2para)*xi2m1)/18375.) + (32*f*(3*f*bη - 5*bΠ2para)*xi30*((35*b1 - 13*f -
                13*f*bη)*xi2m1 + 8*f*(1 + bη)*xi4m1))/42875.

    M[1,4,3] = @. (-8*f*xi10*(33*(840*f*bK²*(1 + bη) - 70*bKKpara*(21*b1 - 17*f - 17*f*bη) + 9*(252*f^2*(-2 + 5*f)*bη^2 -
                    35*f*bδη*(-35*b1 + 51*f + 51*f*bη) + f*bη*(1100*b1 + 1225*b2 - 504*f - 3255*b1*f + 1260*f^2 +
                        1260*f^2*bη² + 1890*bΠ2para) + 5*(f*(220*b1 + 245*b1^2 + 245*b2 - 357*b1*f) - 42*(7*b1 -
                            6*f)*f^2*bη² + (-196*b1 + 378*f)*bΠ2para)))*xi2m1 - 8*f*(1 + bη)*(20790*bK² + 11935*bKKpara -
                3*(4950*b1 + 3850*b1*f + 3850*f*bδη - 2268*f*bη - 4620*f^2*bη - 4620*f^2*bη² -
                    2765*bΠ2para))*xi4m1))/1.2733875e7 + xi30*((-8*f*(4620*f*bK²*(1 + bη) + 3080*f*bKKpara*(1 + bη) +
                3*(275*b1*(-4 + 21*b1 - 7*f)*f + 42*f^2*(12 + 55*f)*bη^2 - 1925*f*bδη*(-3*b1 + f + f*bη) - 770*(5*b1 -
                        3*f)*f^2*bη² - 70*(55*b1 - 21*f)*bΠ2para + f*bη*(-1100*b1 + 504*f - 5775*b1*f + 2310*f^2 +
                        2310*f^2*bη² + 1470*bΠ2para)))*xi0m1)/606375. - (64*f*((-13860*f*bK²*(1 + bη) +
                    385*bKKpara*(63*b1 - 31*f - 31*f*bη) + 3*(55*b1*f*(60 + 49*f) - 84*f^2*(18 + 55*f)*bη^2 +
                        2695*f^2*bδη*(1 + bη) + 770*(7*b1 - 6*f)*f^2*bη² + 105*(77*b1 - 31*f)*bΠ2para - 3*f*bη*(-1100*b1 +
                            504*f - 2695*b1*f + 1540*f^2 + 1540*f^2*bη² + 1085*bΠ2para)))*xi2m1 + 9*f*(1 + bη)*(-825*b1 +
                    1155*bK² + 770*bKKpara + 378*f*bη + 420*f^2*bη + 420*f^2*bη² + 140*bΠ2para)*xi4m1))/1.2733875e7)

    M[1,4,5] = @. (-8*f^2*(1 + bη)*(2*bKKpara - 3*(5*b1*f + 5*f*bδη - 6*f^2*bη - 6*f^2*bη² -
                4*bΠ2para))*xi104m1)/2079. + (8*f^2*xi30*(65*(-4*bKKpara*(1 + bη) + 3*(7*b1*f - 12*f^2*bη^2 +
                    7*f*bδη*(1 + bη) + 14*b1*f*bη² - 12*f^2*bη² + 3*bη*(7*b1*f - 4*f^2 - 4*f^2*bη² - 2*bΠ2para) -
                    6*bΠ2para))*xi2m1 + 36*(1 + bη)*(13*bKKpara + 12*f^2*bη + 12*f^2*bη² + 13*bΠ2para)*xi4m1))/189189.

    M[1,6,3] = @. xi30*((-64*f^2*(1 + bη)*(3*f*bη - 5*bΠ2para)*xi0m1)/24255. + (128*f^2*(1 + bη)*(3*f*bη -
                5*bΠ2para)*(4*xi2m1 - 3*xi4m1))/169785.) + (64*f^2*(1 + bη)*(3*f*bη - 5*bΠ2para)*xi10*(11*xi2m1 +
            4*xi4m1))/56595.

    M[1,6,5] = @. (-16*f^2*(1 + bη)*(2*bKKpara - 3*(5*b1*f + 5*f*bδη - 6*f^2*bη - 6*f^2*bη² -
                4*bΠ2para))*xi104m1)/3465. + (16*f^2*xi30*(65*(-4*bKKpara*(1 + bη) + 3*(7*b1*f - 12*f^2*bη^2 +
                    7*f*bδη*(1 + bη) + 14*b1*f*bη² - 12*f^2*bη² + 3*bη*(7*b1*f - 4*f^2 - 4*f^2*bη² - 2*bΠ2para) -
                    6*bΠ2para))*xi2m1 + 36*(1 + bη)*(13*bKKpara + 12*f^2*bη + 12*f^2*bη² + 13*bΠ2para)*xi4m1))/315315.

    M[1,6,7] = @. (-224*f^4*(1 + bη)*(bη + bη²)*xi304m1)/6435.

    M[1,8,7] = @. (-256*f^4*(1 + bη)*(bη + bη²)*xi304m1)/6435.

    M[2,0,0] = @. (f^2*(15*b1^2 - 10*b1*f + 3*f^2 + 2*f*(-5*b1 + 3*f)*bη + 3*f^2*bη^2)*xi0m200)/270. +
    ((-28*bKKpara*(90*b1 + 17*f^2 - 30*f*bη) - 840*bK²*(9*b1 + f^2 - 3*f*bη) - 9*(900*b1^2 + 1575*b1*b2 - 1050*b1^2*f +
                100*b1*f^2 + 175*b1^2*f^2 + 175*b2*f^2 - 420*b1*f^3 + 63*f^4 + f^2*(108 - 154*f + 63*f^2)*bη^2 -
                70*f*bδη*(3*(5*b1 + f^2) - 5*f*bη) + 462*b1*f^2*bη² + 126*f^4*bη² + 930*b1*bΠ2para + 214*f^2*bΠ2para -
                f*bη*(480*b1 + 525*b2 - 812*b1*f + 36*f^2 + 210*b1*f^2 - 252*f^3 + 154*f^2*bη² +
                    390*bΠ2para)))*xi1m1sq)/28350. + (4*f^3*(1 + bη)*(-7*b1 + 3*f + 3*f*bη)*xi2m220)/945. +
    (2*(225*b1^2 + 25*b1*f^2 - 3*f^4 - 120*b1*f*bη - 42*b1*f^2*bη - 9*f^3*bη - 12*f^4*bη + 27*f^2*bη^2 + 14*f^3*bη^2 -
            3*f^4*bη^2 - 7*bKKpara*(15*b1 + 2*f^2 - 5*f*bη) - 35*bK²*(9*b1 + f^2 - 3*f*bη) - 42*b1*f^2*bη² - 6*f^4*bη² +
            14*f^3*bη*bη² - 30*b1*bΠ2para + f^2*bΠ2para - 10*f*bη*bΠ2para)*xi3m1sq)/1575.

    M[2,0,2] = @. (2*f*(588*f*bK² + 686*bKKpara*(f - 3*bη) + 9*(-490*f*bδη*(3*f - 2*bη) + 49*f^2*(-16 + 9*f)*bη^2 -
                2*bη*(f*(18*(16 - 49*f)*f + 245*b1*(-2 + 3*f)) + 392*f^2*bη² + 441*bΠ2para) + f*(7*(160*b1 + 175*b1^2 +
                        175*b2 - 420*b1*f + 63*f^2) + 882*f^2*bη² + 1744*bΠ2para)))*xi1m1sq)/496125. + (4*f^3*(1 +
            bη)*(-7*b1 + 3*f + 3*f*bη)*xi002m2)/945. - (4*f*(588*f*bK² + 14*bKKpara*(24*f - 7*bη) + 3*(-35*f*bδη*(3*f -
                    7*bη) + 7*f^2*(-26 + 9*f)*bη^2 + bη*(f*(-35*b1*(-7 + 3*f) + 36*f*(2 + 7*f)) - 182*f^2*bη² -
                    196*bΠ2para) + f*(-7*(20*b1 + 30*b1*f - 9*f^2) + 126*f^2*bη² + 62*bΠ2para)))*xi1m13m1)/55125. +
    (16*f*(-35*b1*f + 7*f^3 + 49*f*bK² + 7*bKKpara*(4*f - 7*bη) + 18*f^2*bη + 28*f^3*bη - 28*f^2*bη^2 + 7*f^3*bη^2 +
            14*f^3*bη² - 28*f^2*bη*bη² - 2*f*bΠ2para - 49*bη*bΠ2para)*xi3m1sq)/55125. + xi20*((2*f^2*(15*b1^2 - 10*b1*f +
                3*f^2 + 2*f*(-5*b1 + 3*f)*bη + 3*f^2*bη^2)*xi0m2)/675. - (16*f^3*(1 + bη)*(-7*b1 + 3*f +
                3*f*bη)*xi2m2)/6615. + (32*f^4*(1 + bη)^2*xi4m2)/11025.)

    M[2,2,0] = @. (f^2*(15*b1^2 - 10*b1*f + 3*f^2 + 2*f*(-5*b1 + 3*f)*bη + 3*f^2*bη^2)*xi0m200)/135. +
    ((-196*f*bKKpara*(17*f - 30*bη) - 5880*f*bK²*(f - 3*bη) - 9*(700*b1*f^2 + 1225*b1^2*f^2 + 1225*b2*f^2 - 2940*b1*f^3 +
                441*f^4 - 490*f^2*bδη*(3*f - 5*bη) + f^2*(1080 - 1078*f + 441*f^2)*bη^2 + 882*f^4*bη² + 2100*b1*bΠ2para +
                1678*f^2*bΠ2para - f*bη*(3360*b1 + 3675*b2 - 2450*b1*f + 360*f^2 + 1470*b1*f^2 - 1764*f^3 + 1078*f^2*bη² +
                    3270*bΠ2para)))*xi1m1sq)/99225. + (8*f^3*(1 + bη)*(-7*b1 + 3*f + 3*f*bη)*xi2m220)/945. +
    (4*(175*b1*f^2 - 21*f^4 - 49*f*bKKpara*(2*f - 5*bη) - 245*f*bK²*(f - 3*bη) - 840*b1*f*bη - 90*f^3*bη - 84*f^4*bη +
            270*f^2*bη^2 + 98*f^3*bη^2 - 21*f^4*bη^2 - 42*f^4*bη² + 98*f^3*bη*bη² + 525*b1*bΠ2para + 52*f^2*bΠ2para -
            205*f*bη*bΠ2para)*xi3m1sq)/11025.

    M[2,2,2] = @. ((924*f^2*bK² + 98*bKKpara*(63*b1 + 11*f^2 - 33*f*bη) + 9*(-2940*b1^2*f + 1760*b1*f^2 + 1925*b1^2*f^2 +
                1925*b2*f^2 - 4620*b1*f^3 + 693*f^4 + 77*f^3*(-16 + 9*f)*bη^2 - 70*f*bδη*(42*b1 + 33*f^2 - 22*f*bη) +
                2352*b1*f^2*bη² + 1386*f^4*bη² + 2646*b1*bΠ2para + 2672*f^2*bΠ2para - 2*f*bη*(f*(18*(24 - 77*f)*f +
                        7*b1*(-278 + 165*f)) + 616*f^2*bη² + 693*bΠ2para)))*xi1m1sq)/99225. + (22*f^3*(1 + bη)*(-7*b1 +
            3*f + 3*f*bη)*xi002m2)/1323. + (2*(-924*f^2*bK² - 2*bKKpara*(147*b1 + 264*f^2 - 77*f*bη) + 3*(735*b1^2*f +
                220*b1*f^2 + 330*b1*f^3 - 99*f^4 + 11*(26 - 9*f)*f^3*bη^2 + 5*f*bδη*(147*b1 + 33*f^2 - 77*f*bη) -
                546*b1*f^2*bη² - 198*f^4*bη² - 588*b1*bΠ2para - 106*f^2*bΠ2para + f*bη*(-(f*(b1*(931 - 165*f) + 36*f*(3 +
                                11*f))) + 286*f^2*bη² + 308*bΠ2para)))*xi1m13m1)/11025. + (8*(-55*b1*f^2 + 11*f^4 +
            77*f^2*bK² + 84*b1*f^2*bη + 27*f^3*bη + 44*f^4*bη - 44*f^3*bη^2 + 11*f^4*bη^2 + bKKpara*(147*b1 + 44*f^2 -
                77*f*bη) + 84*b1*f^2*bη² + 22*f^4*bη² - 44*f^3*bη*bη² + 147*b1*bΠ2para - f^2*bΠ2para -
            77*f*bη*bΠ2para)*xi3m1sq)/11025. + xi20*((11*f^2*(15*b1^2 - 10*b1*f + 3*f^2 + 2*f*(-5*b1 + 3*f)*bη +
                3*f^2*bη^2)*xi0m2)/945. - (88*f^3*(1 + bη)*(-7*b1 + 3*f + 3*f*bη)*xi2m2)/9261. + (176*f^4*(1 +
                bη)^2*xi4m2)/15435.)

    M[2,2,4] = @. (32*f^2*(2*bKKpara - 3*(10*b1*f - 3*f^2 + 5*f*bδη + (4 - 3*f)*f*bη^2 - 6*f^2*bη² + f*bη*(5*b1 - 12*f +
                    4*bη²) - 4*bΠ2para))*xi1m13m1)/11025. - (16*f^2*(9*f^2 + 22*bKKpara + f*(-22 + 9*f)*bη^2 +
            2*f*bη*(18*f - 11*bη²) + 18*f^2*bη² + 22*bΠ2para)*xi3m1sq)/40425. + (16*f^4*(1 + bη)^2*xi004m2)/2205. +
    xi20*((32*f^3*(1 + bη)*(-7*b1 + 3*f + 3*f*bη)*xi2m2)/5145. - (128*f^4*(1 + bη)^2*xi4m2)/33957.)

    M[2,4,0] = @. (16*f*(f - 3*bη)*(3*f*bη - 5*bΠ2para)*xi1m1sq)/3675. - (16*f*(f - 3*bη)*(3*f*bη -
            5*bΠ2para)*xi3m1sq)/3675.

    M[2,4,2] = @. (4*f*(924*f*bK² + 1078*bKKpara*(f - 3*bη) - 6930*f*bδη*(3*f - 2*bη) + 693*f^2*(-16 + 9*f)*bη^2 -
            18*bη*(f*(2*(272 - 693*f)*f + 385*b1*(-2 + 3*f)) + 616*f^2*bη² + 693*bΠ2para) + 3*f*(33*(160*b1 + 175*b1^2 +
                    175*b2 - 420*b1*f + 63*f^2) + 4158*f^2*bη² + 9136*bΠ2para))*xi1m1sq)/606375. + (8*f^3*(1 + bη)*(-7*b1 +
            3*f + 3*f*bη)*xi002m2)/735. - (8*f*(-660*b1*f - 990*b1*f^2 + 297*f^3 + 924*f*bK² - 495*f^2*bδη +
            22*bKKpara*(24*f - 7*bη) + 1155*b1*f*bη + 408*f^2*bη - 495*b1*f^2*bη + 1188*f^3*bη + 1155*f*bδη*bη -
            858*f^2*bη^2 + 297*f^3*bη^2 + 594*f^3*bη² - 858*f^2*bη*bη² + 178*f*bΠ2para -
            924*bη*bΠ2para)*xi1m13m1)/67375. + (32*f*(-165*b1*f + 33*f^3 + 231*f*bK² + 33*bKKpara*(4*f - 7*bη) +
            102*f^2*bη + 132*f^3*bη - 132*f^2*bη^2 + 33*f^3*bη^2 + 66*f^3*bη² - 132*f^2*bη*bη² - 38*f*bΠ2para -
            231*bη*bΠ2para)*xi3m1sq)/202125. + xi20*((4*f^2*(15*b1^2 - 10*b1*f + 3*f^2 + 2*f*(-5*b1 + 3*f)*bη +
                3*f^2*bη^2)*xi0m2)/525. - (32*f^3*(1 + bη)*(-7*b1 + 3*f + 3*f*bη)*xi2m2)/5145. + (64*f^4*(1 +
                bη)^2*xi4m2)/8575.)

    M[2,4,4] = @. (8*f^2*(-390*b1*f + 117*f^2 + 26*bKKpara - 195*f*bδη + 308*b1*bη - 195*b1*f*bη + 468*f^2*bη -
            156*f*bη^2 + 117*f^2*bη^2 + 308*b1*bη² + 234*f^2*bη² - 156*f*bη*bη² + 156*bΠ2para)*xi1m13m1)/13475. -
    (4*f^2*(351*f^2 + 858*bKKpara + 39*f*(-22 + 9*f)*bη^2 + 1694*b1*bη² + 702*f^2*bη² + 2*bη*(847*b1 + 702*f^2 -
                429*f*bη²) + 858*bΠ2para)*xi3m1sq)/148225. + (52*f^4*(1 + bη)^2*xi004m2)/2695. + xi20*((312*f^3*(1 +
                bη)*(-7*b1 + 3*f + 3*f*bη)*xi2m2)/18865. - (416*f^4*(1 + bη)^2*xi4m2)/41503.)

    M[2,4,6] = @. (80*f^4*(1 + 4*bη + bη^2 + 2*bη²)*xi3m1sq)/11011. + (80*f^4*(1 + bη)^2*xi204m2)/11011.

    M[2,6,2] = @. (-256*f^2*(3*f*bη - 5*bΠ2para)*xi1m1sq)/40425. - (64*f^2*(3*f*bη - 5*bΠ2para)*xi1m13m1)/13475. +
    (64*f^2*(3*f*bη - 5*bΠ2para)*xi3m1sq)/40425.

    M[2,6,4] = @. (16*f^2*(2*bKKpara - 3*(10*b1*f - 3*f^2 + 5*f*bδη + (4 - 3*f)*f*bη^2 - 6*f^2*bη² + f*bη*(5*b1 - 12*f +
                    4*bη²) - 4*bΠ2para))*xi1m13m1)/3465. - (8*f^2*(9*f^2 + 22*bKKpara + f*(-22 + 9*f)*bη^2 +
            2*f*bη*(18*f - 11*bη²) + 18*f^2*bη² + 22*bΠ2para)*xi3m1sq)/12705. + (8*f^4*(1 + bη)^2*xi004m2)/693. +
    xi20*((16*f^3*(1 + bη)*(-7*b1 + 3*f + 3*f*bη)*xi2m2)/1617. - (320*f^4*(1 + bη)^2*xi4m2)/53361.)

    M[2,6,6] = @. (664*f^4*(1 + 4*bη + bη^2 + 2*bη²)*xi3m1sq)/38115. + (664*f^4*(1 + bη)^2*xi204m2)/38115.

    M[2,8,6] = @. (64*f^4*(1 + 4*bη + bη^2 + 2*bη²)*xi3m1sq)/6435. + (64*f^4*(1 + bη)^2*xi204m2)/6435.

    M[3,0,1] = @. (f*(25*b1^2 + 15*b1*(-1 + f)*f - 9*f^3 - 3*f*(10*b1 + 3*(-1 + f)*f)*bη + 9*f^2*bη^2)*xi0m11m2)/225. +
    (2*f*(-125*b1^2 + 15*b1*(3 - 2*f)*f + 18*f^3 + 3*f*(40*b1 + 3*f*(-3 + 2*f))*bη - 27*f^2*bη^2)*xi1m22m1)/1125. -
    (4*f^2*(1 + bη)*(35*b1 + 9*f^2 - 21*f*bη)*xi2m13m2)/2625.

    M[3,0,3] = @. (4*f^3*(5*b1 - 3*f + (2 - 3*f)*bη + 2*bη^2)*xi1m22m1)/875. - (4*f^3*(f - bη)*(1 +
            bη)*xi0m13m2)/525. + (16*f^3*(2*f - 3*bη)*(1 + bη)*xi2m13m2)/7875.

    M[3,2,1] = @. (2*f*(35*b1^2 - 18*f^3 + 3*b1*f*(-7 + 10*f) - 3*f*(17*b1 + 6*(-1 + f)*f)*bη +
            18*f^2*bη^2)*xi0m11m2)/315. + (4*f*(-175*b1^2 + 3*b1*(21 - 20*f)*f + 36*f^3 + 3*f*(71*b1 + 6*f*(-3 +
                    2*f))*bη - 54*f^2*bη^2)*xi1m22m1)/1575. - (8*f^2*(1 + bη)*(49*b1 + 18*f^2 -
            42*f*bη)*xi2m13m2)/3675.

    M[3,2,3] = @. (2*f^2*(-21*f^2 + b1*(-18 + 35*f) + (-18*b1 + 7*(2 - 3*f)*f)*bη + 14*f*bη^2)*xi1m22m1)/525. -
    (2*f^2*(1 + bη)*(9*b1 + 7*f^2 - 7*f*bη)*xi0m13m2)/315. + (8*f^2*(1 + bη)*(27*b1 + 14*f^2 -
            21*f*bη)*xi2m13m2)/4725.

    M[3,2,5] = @. (-160*f^4*(1 + bη)*xi2m13m2)/14553.

    M[3,4,1] = @. (-8*f^2*(f - bη)*(-5*b1 + 3*f + 3*f*bη)*xi0m11m2)/525. + (16*f^2*(2*f*(-5*b1 + 3*f) + (25*b1 +
                3*f*(-3 + 2*f))*bη - 9*f*bη^2)*xi1m22m1)/2625. - (32*f^3*(3*f - 7*bη)*(1 + bη)*xi2m13m2)/6125.

    M[3,4,3] = @. (16*f^2*(-63*f^2 + 5*b1*(-11 + 21*f) + (-55*b1 + 21*(2 - 3*f)*f)*bη + 42*f*bη^2)*xi1m22m1)/9625. -
    (8*f^2*(1 + bη)*(55*b1 + 42*f^2 - 42*f*bη)*xi0m13m2)/5775. + (32*f^2*(1 + bη)*(55*b1 + 28*f^2 -
            42*f*bη)*xi2m13m2)/28875.

    M[3,4,5] = @. (-920*f^4*(1 + bη)*xi2m13m2)/21021.

    M[3,6,3] = @. (16*f^3*(5*b1 - 3*f + (2 - 3*f)*bη + 2*bη^2)*xi1m22m1)/1155. - (16*f^3*(f - bη)*(1 +
            bη)*xi0m13m2)/693. + (64*f^3*(2*f - 3*bη)*(1 + bη)*xi2m13m2)/10395.

    M[3,6,5] = @. (-544*f^4*(1 + bη)*xi2m13m2)/10395.

    M[3,8,5] = @. (-128*f^4*(1 + bη)*xi2m13m2)/6435.

    M[4,0,0] = @. ((45*b1^2 + 10*b1*f^2 + 3*f^4 - 6*(5*b1*f + f^3)*bη + 9*f^2*bη^2)*xi0m2sq)/540. + ((225*b1^2 +
            50*b1*f^2 + 6*f^4 - 30*(5*b1*f + f^3)*bη + 45*f^2*bη^2)*xi2m2sq)/1350.

    M[4,0,2] = @. (4*f^2*(7*b1 + 3*f^2 - 6*f*bη)*xi0m22m2)/945. - (2*f^2*(49*b1 + 12*f^2 - 42*f*bη)*xi2m2sq)/6615.

    M[4,0,4] = @. (16*f^4*xi2m2sq)/11025.

    M[4,2,0] = @. (f*(7*b1*f + 3*f^3 - 3*(7*b1 + 2*f^2)*bη + 9*f*bη^2)*xi0m2sq)/189. + (2*f*(35*b1*f + 6*f^3 - 15*(7*b1 +
                2*f^2)*bη + 45*f*bη^2)*xi2m2sq)/945.

    M[4,2,2] = @. (f^2*(22*b1 + 9*f^2 - 18*f*bη)*xi0m22m2)/189. - (f^2*(77*b1 + 18*f^2 - 63*f*bη)*xi2m2sq)/1323.

    M[4,2,4] = @. (272*f^4*xi2m2sq)/24255.

    M[4,4,0] = @. (2*f^2*(f^2 - 2*f*bη + 3*bη^2)*xi0m2sq)/315. + (4*f^2*(2*f^2 - 10*f*bη + 15*bη^2)*xi2m2sq)/1575.

    M[4,4,2] = @. (8*f^2*(33*b1 + 17*f^2 - 34*f*bη)*xi0m22m2)/3465. - (4*f^2*(231*b1 + 68*f^2 -
            238*f*bη)*xi2m2sq)/24255.

    M[4,4,4] = @. (3858*f^4*xi2m2sq)/175175.

    M[4,6,2] = @. (8*f^3*(f - 2*bη)*xi0m22m2)/693. - (8*f^3*(2*f - 7*bη)*xi2m2sq)/4851.

    M[4,6,4] = @. (424*f^4*xi2m2sq)/24255.

    M[4,8,4] = @. (32*f^4*xi2m2sq)/6435.

   return r, M

end

max_stars_repo_name: JosephTomlinson/GeneralBiasPk

max_issues_repo_name: JosephTomlinson/GeneralBiasPk

max_forks_repo_name: JosephTomlinson/GeneralBiasPk


text:
import RvSpectMLBase
using Test

@testset "Spectra Timeseries" begin
    using RvSpectMLBase.TheoreticalInstrument
    using DataFrames
    num_pixels = 1000
    num_orders = 10
    num_obs = 3
    lambda2d = collect(reshape(range(5000.0,stop=6000.0,length=num_pixels*num_orders),num_pixels,num_orders))
    inst = TheoreticalInstrument2D(lambda2d)
    times = sort(rand(num_obs))
    rvs = range(0.0, stop=10.0, length=length(times))
    df = DataFrame(:lambda=>[5501], :weight=>[0.5])
    @test_nowarn TheoreticalInstrument.generate_spectra_timeseries(times,df,inst, rvs)
    spectra = TheoreticalInstrument.generate_spectra_timeseries(times,df,inst, rvs)
end

ext: jl

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
using StatsFuns

@testset "Reasonable return from varprop function." begin
    sigma = Matrix{Int}(I, 3, 3)
    recAbund = ones(5)

    retval = runCkineVarProp(tps, rxntfR, sigma, recAbund)

    @test length(retval) == length(tps)
    @test all(retval .>= 0.0)
end

@testset "Run fitting." begin
    p = gcSolver.getUnkVec()
    outtG = similar(p)

    outt = gcSolver.resids(p)
    ForwardDiff.gradient!(outtG, gcSolver.resids, p)

    @time ForwardDiff.gradient!(outtG, gcSolver.resids, p)

    @test isfinite(outt)
    @test all(isfinite.(outtG))
end

@testset "Test that residuals can be calculated using Fit." begin
    fitVec = gcSolver.importFit()
    fitVec = convert(Vector{Float64}, fitVec[!, :Fit])
    fitVec = softplus.(fitVec)
    resids = gcSolver.resids(fitVec)
    @test isfinite(resids)
end


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

"""
    EmpiricalDistribution(values)

An empirical distribution holding continuous values.
"""
struct EmpiricalDistribution{T<:Real} <: ContinuousUnivariateDistribution
  values::Vector{T}

  function EmpiricalDistribution{T}(values) where {T<:Real}
    @assert !isempty(values) "values must be provided"
    new(sort(values))
  end
end

EmpiricalDistribution(values) = EmpiricalDistribution{eltype(values)}(values)

quantile(d::EmpiricalDistribution, p::Real) = quantile(d.values, p, sorted=true)

function cdf(d::EmpiricalDistribution{T}, x::T) where {T<:Real}
  v = d.values
  n = length(v)

  head, mid, tail = 1, 1, n
  while tail - head > 1
    mid = (head + tail) ÷ 2
    if x < v[mid]
      tail = mid
    else
      head = mid
    end
  end

  l, u = v[head], v[tail]

  if x < l
    return zero(T)
  elseif x > u
    return one(T)
  else
    if head == tail
      return head / n
    else
      pl, pu = head / n, tail / n
      return (pu - pl) * (x - l) / (u - l) + pl
    end
  end
end


text:
#=
   The Piranha Puzzle in Turing.jl
   (From Tijms 2004)

   http://cs.ioc.ee/ewscs/2020/katoen/katoen-slides-lecture1.pdf
   """
   One fish is contained withing the confines of an opaque fishbowl.
   The fish is equally likely to be a piranha or a goldfish. A sushi
   lover throws a piranha into the fish bowl alongside the other fish.
   Then, immediately, before either fish can devour the other, one
   of the fish is blindly removed from the fishbowl. The fish that
   has been removed from the bowl turns out to be a piranha. What
   is the probability that the fish that was originaly in the bowl
   by itself was a piranha?

   The Piranha Puzzle Program

   f1 := gf [0.5] f1 := pir;
   f2 := pir;
   s := f1 [0.5] s := f2;
   observe(s = pir)

   E(f1 = pir | P terminates) = 1/2 / 3/4 = 2/3
   """

   Distributions of variable f1
   piranha    =>    6596  (0.659600)
   goldfish   =>    3404  (0.340400)

   The exact probability of first fish is a piranha is 0.6666

   Cf ~/webppl/piranha_puzzle.wppl

=#

using Turing
include("jl_utils.jl")

@model function piranha_puzzle()

    goldfish = 1
    piranha = 2

    f1 ~ Categorical([0.5,0.5])
    f2 = piranha
    fs = [f1,f2]
    s ~ Categorical([0.5,0.5])
    true ~ Dirac(fs[s] == piranha)

    f1_piranha ~ Dirac(f1 == piranha)

end

model = piranha_puzzle()

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(15), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
println("goldfish: 1 piranha: 2")
show_var_dist_pct(chns,:f1,["goldfish","piranha"])

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
"""
A basis of Legendre polynomials on the interval `[-1,1]`. These classical
polynomials are orthogonal with respect to the weight function `w(x) = 1`.
"""
struct Legendre{T} <: OPS{T}
    n   ::  Int
end

Legendre(n::Int) = Legendre{Float64}(n)

similar(b::Legendre, ::Type{T}, n::Int) where {T} = Legendre{T}(n)

support(b::Legendre{T}) where {T} = ChebyshevInterval{T}()

first_moment(b::Legendre{T}) where {T} = T(2)

jacobi_α(b::Legendre{T}) where {T} = T(0)
jacobi_β(b::Legendre{T}) where {T} = T(0)

measure(b::Legendre{T}) where {T} = LegendreWeight{T}()
isorthogonal(::Legendre, ::LegendreWeight) = true
issymmetric(::Legendre) = true
interpolation_grid(dict::Legendre{T}) where T = LegendreNodes{T}(length(dict))
iscompatible(dict::Legendre, grid::LegendreNodes) = length(dict) == length(grid)
isorthogonal(dict::Legendre, measure::GaussLegendre) = opsorthogonal(dict, measure)

gauss_rule(dict::Legendre{T}) where T = GaussLegendre{T}(length(dict))

function innerproduct_native(d1::Legendre, i::PolynomialDegree, d2::Legendre, j::PolynomialDegree, m::LegendreWeight; options...)
	T = coefficienttype(d1)
	if i == j
		2 / convert(T, 2*value(i)+1)
	else
		zero(T)
	end
end

# See DLMF, Table 18.9.1
# http://dlmf.nist.gov/18.9#i
rec_An(b::Legendre{T}, n::Int) where {T} = T(2*n+1)/T(n+1)
rec_Bn(b::Legendre{T}, n::Int) where {T} = zero(T)
rec_Cn(b::Legendre{T}, n::Int) where {T} = T(n)/T(n+1)

show(io::IO, b::Legendre{Float64}) = print(io, "Legendre($(length(b)))")
show(io::IO, b::Legendre{T}) where T = print(io, "Legendre{$(T)}($(length(b)))")


text:
using CairoMakie, Random, LaTeXStrings
Random.seed!(123)
m = 200
function spiral(; a = 1, n = 100, h = 0, k = 0)
    φ = LinRange(rand() + 1, 6π, n)
    h .+  rand(-1:2:1)*a*cos.(φ)./φ, k .+  rand(-1:2:1)*a*sin.(φ)./φ
end
x = raw"x=a\,\cos(\varphi)/\varphi,\quad "
y = raw"y=a\,\sin(\varphi)/\varphi,\quad "
curves = [spiral(; a = rand(), h = rand(-1:1)) for i in 1:m]

with_theme(theme_dark()) do
    fig = Figure(resolution = (600,400))
    ax = Axis(fig[1,1])
    series!(curves; color = to_colormap(:bone_1, m), linewidth = 1.5)
    text!(latexstring(x * y * "\\varphi>0"), position = (-1,-1))
    hidedecorations!(ax; grid = false)
    save("series.png", fig, px_per_unit = 2)
    display(fig)
end


text:
module DeltaRobot

include("animation.jl")
include("delta_robot_params.jl")
include("kinematics.jl")
include("simulation.jl")
include("stats.jl")

export transform_to_cartesian, animate_arm, animate_robot
export DeltaRobotParams
export delta_arm_ics, delta_robot_ics, delta_arm!, L_delta_robot!, L_delta_robot_base!
export TABLE_HEADER, summary_data

end

max_stars_repo_name: zborffs/Delta

max_issues_repo_name: zborffs/Delta

max_forks_repo_name: zborffs/Delta


text:
module PaperSimulation
    using OceanWaveSpectralFitting, Optim, ProgressMeter
    import OceanWaveSpectralFitting: WhittleLikelihoodInference.GaussianProcess
    const WindSea = JS_BWG_HNE{0} # simulation with no aliasing
    include("../freqdirspectra/FrequencyDirectionSpectra.jl")
    include("leastsquares.jl")
    include("methodofmoments.jl")
    include("whittle.jl")

    ##
    "Helper function to initialise peak"
    function initialisepeak(ts::Matrix{Float64}, Δ, ωₗ, ωᵤ)
        P = OceanWaveSpectralFitting.Periodogram(ts,Δ)
        S = P.ordinates[end÷2+1:end]
        Ω = P.Ω[end÷2+1:end]
        peakind = argmax([(ωₗ < ω < ωᵤ) ? real(s[1,1]) : -Inf for (s,ω) ∈ zip(S, Ω)])
        meandir = angle(imag(S[peakind][2,1])+1im*imag(S[peakind][3,1]))
        return Ω[peakind], mod2pi.(meandir.+pi) .- pi
    end
    ##
    "Simulation study struct"
    struct SimulationStudy
        LSmlm::Vector{Vector{Float64}}
        LSmem::Vector{Vector{Float64}}
        M::Vector{Vector{Float64}}
        DW::Vector{Vector{Float64}}
        DWuni::Vector{Vector{Float64}}
        LSmlmflag::Vector{Bool}
        LSmemflag::Vector{Bool}
        Mflag::Vector{Bool}
        DWflag::Vector{Bool}
        DWuniflag::Vector{Bool}
        θ_true::Vector{Float64}
        cutoff::Float64
        function SimulationStudy(
            nreplications::Int,
            θ::Vector{Float64};
            ωₗ::Float64 = 0.0,
            name = "Simulation ",
        )
            @assert length(θ) == 9
            Δ = 1 / 1.28
            X = GaussianProcess(WindSea(θ), 2304, Δ)
            LSmlm = Vector{Vector{Float64}}(undef, nreplications)
            LSmlmflag = Vector{Bool}(undef, nreplications)
            LSmem = similar(LSmlm)
            M = similar(LSmlm)
            DW = similar(LSmlm)
            DWuni = similar(LSmlm)
            LSmemflag = similar(LSmlmflag)
            Mflag = similar(LSmlmflag)
            DWflag = similar(LSmlmflag)
            DWuniflag = similar(LSmlmflag)
            @showprogress name for rep ∈ 1:nreplications
                ts = rand(X).ts
                peak, dir = initialisepeak(ts, Δ, ωₗ, π / Δ)
                θ_ini = [0.7, peak, 2, 4.5, dir, 4.5, 3, 0.6, 0.3]
                LSmlmflag[rep], LSmlm[rep] =
                    LSfit(FrequencyDirectionSpectra.FDJONSWAP, ts, Δ, θ_ini, fdest = FrequencyDirectionSpectra.MLM; ωₗ = ωₗ)
                LSmemflag[rep], LSmem[rep] =
                    LSfit(FrequencyDirectionSpectra.FDJONSWAP, ts, Δ, θ_ini, fdest = FrequencyDirectionSpectra.MEM; ωₗ = ωₗ)
                Mflag[rep], M[rep] = momentsfit(FrequencyDirectionSpectra.FDJONSWAP, ts, Δ, θ_ini; ωₗ = ωₗ)
                DWflag[rep], DW[rep] = DWfit(WindSea, ts, Δ, θ_ini; ωₗ = ωₗ)
                DWuniflag[rep], DWuni[rep] = DWfit(JONSWAP{0}, ts[:,1], Δ, θ_ini[1:4]; ωₗ = ωₗ)
            end
            new(
                LSmlm,
                LSmem,
                M,
                DW,
                DWuni,
                LSmlmflag,
                LSmemflag,
                Mflag,
                DWflag,
                DWuniflag,
                θ,
                ωₗ,
            )
        end
    end
    Base.show(io::IO, sim::SimulationStudy) = print(io,"Simulation study results.")

end


text:
module C

const NAMES = [
    "V1"
    "Km1"
    "V2"
    "Km2"
    "V3"
    "Km3"
    "V4"
    "Km4"
    "V5"
    "Km5"
    "V6"
    "Km6"
    "KimERK"
    "KexERK"
    "KimpERK"
    "KexpERK"
    "KimppERK"
    "KexppERK"
    "V10"
    "Km10"
    "n10"
    "p11"
    "p12"
    "p13"
    "V14"
    "Km14"
    "V15"
    "Km15"
    "p16"
    "p17"
    "KimDUSP"
    "KexDUSP"
    "KimpDUSP"
    "KexpDUSP"
    "V20"
    "Km20"
    "V21"
    "Km21"
    "p22"
    "p23"
    "V24"
    "Km24"
    "V25"
    "Km25"
    "KimRSK"
    "KexRSK"
    "V27"
    "Km27"
    "V28"
    "Km28"
    "V29"
    "Km29"
    "V30"
    "Km30"
    "V31"
    "Km31"
    "n31"
    "p32"
    "p33"
    "p34"
    "V35"
    "Km35"
    "V36"
    "Km36"
    "V37"
    "Km37"
    "p38"
    "p39"
    "KimFOS"
    "KexFOS"
    "KimpcFOS"
    "KexpcFOS"
    "V42"
    "Km42"
    "V43"
    "Km43"
    "V44"
    "Km44"
    "p45"
    "p46"
    "p47"
    "m47"
    "p48"
    "p49"
    "m49"
    "p50"
    "p51"
    "m51"
    "p52"
    "m52"
    "p53"
    "p54"
    "m54"
    "p55"
    "p56"
    "m56"
    "V57"
    "Km57"
    "n57"
    "p58"
    "p59"
    "p60"
    "p61"
    "KimF"
    "KexF"
    "p63"
    "KF31"
    "nF31"
    #
    "a"
    "Vn"
    "Vc"
    "Ligand"
    "EGF"
    "HRG"
    "no_ligand"
]

for (idx,name) in enumerate(NAMES)
    eval(Meta.parse("const $name = $idx"))
end

const NUM = length(NAMES)

end  # module

max_stars_repo_name: himoto/paramEstim

max_issues_repo_name: himoto/paramEstim


text:
# Implementation of Armijo backtracking linesearch

"""Armijo backtracking linesearch
Armijo(model, x, s, fx, g, t, α) computes a step t satisfying Armijo condition
The function returns (xtrial, fxtrial) : the new iterate for the linesearch and
the value of the objective in xtrial
"""
function Armijo(model, x, s, fx=obj(model, x), g=grad(model, x), t=1., α=1e-4)
    @info(loggerArm, @sprintf("Armijo backtracking linesearch on problem %s", model.meta.name))
    (t ≤ 0) && (t = 1.)
    @debug(loggerArm, @sprintf("initial step t = %7.1e", t))
    (α ≤ 0 || α ≥ 1) && (α = 1e-4)
    @debug(loggerArm, @sprintf("α = %8.1e", α))
    xtrial = x + t * s
    fxtrial = obj(model, xtrial) # f(x + t * s)
    @debug(loggerArm, @sprintf("f(x) = %8.1e", fxtrial))
    αgs = α * dot(g, s)
    β = fx + t * αgs # f(x) + t * α * gᵀs
    @debug(loggerArm, @sprintf("β = %8.1e", β))

    while fxtrial > β && t > eps()
        @debug(loggerArm, @sprintf("f(x) = %8.1e > β = %8.1e", fxtrial, β))
        t = 0.5 * t
        @debug(loggerArm, @sprintf("new t = %7.1e", t))
        xtrial = x + t * s
        fxtrial = obj(model, xtrial)
        β = fx + t * αgs
    end
    @debug(loggerArm, @sprintf("final values : f(x+t*s) = %8.1e ≤ β = %8.1e", fxtrial, β))
    return xtrial, fxtrial
end


text:
struct Point{I <: Integer} <: AbstractShape
    i::I
    j::I
end

get_i_min(shape::Point) = shape.i
get_i_max(shape::Point) = shape.i

get_j_min(shape::Point) = shape.j
get_j_max(shape::Point) = shape.j

get_drawing_optimization_style(::Point) = PUT_PIXEL

function draw!(f::F, image::AbstractMatrix, shape::Point, color) where {F <: Function}
    f(image, shape.i, shape.j, color)
    return nothing
end


text:
# ---
# title: 1190. Reverse Substrings Between Each Pair of Parentheses
# id: problem1190
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Stack
# link: <https://leetcode.com/problems/reverse-substrings-between-each-pair-of-parentheses/description/>
# hidden: true
# ---
#
# You are given a string `s` that consists of lower case English letters and
# brackets.
#
# Reverse the strings in each pair of matching parentheses, starting from the
# innermost one.
#
# Your result should **not** contain any brackets.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "(abcd)"
#     Output: "dcba"
#
#
# **Example 2:**
#
#
#
#     Input: s = "(u(love)i)"
#     Output: "iloveu"
#     Explanation:  The substring "love" is reversed first, then the whole string is reversed.
#
#
# **Example 3:**
#
#
#
#     Input: s = "(ed(et(oc))el)"
#     Output: "leetcode"
#     Explanation:  First, we reverse the substring "oc", then "etco", and finally, the whole string.
#
#
# **Example 4:**
#
#
#
#     Input: s = "a(bcdefghijkl(mno)p)q"
#     Output: "apmnolkjihgfedcbq"
#
#
#
#
# **Constraints:**
#
#   * `0 <= s.length <= 2000`
#   * `s` only contains lower case English characters and parentheses.
#   * It's guaranteed that all parentheses are balanced.
#
#
## @lc code=start
using LeetCode

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


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

using Markdown
using InteractiveUtils

# ╔═╡ d54abd5a-5f44-11eb-34fe-c96ecbd1327c
begin
	# Load required packages
	using Pkg, HTTP, CSV, DataFrames, Random, Distances,
        UMAP, Clustering, Gadfly, DelimitedFiles,
        MultivariateStats, Cairo, Fontconfig, StatsBase

	# Set default plot size
    set_default_plot_size(50cm, 20cm)

	# Set working directory
	cd("/Users/martintreppner/Desktop")
end;

# ╔═╡ 59a1670e-5fb6-11eb-0cf1-6fefc45fb083
begin
	# Install and Load single-cell deep Boltzmann machine (scDBM)
	# Pkg.add(PackageSpec(url="https://github.com/MTreppner/scDBM.jl", rev="master"))
	using scDBM
end;

# ╔═╡ 09c53f44-5fb6-11eb-121b-757954f80b1d
begin
	# Load Python packages using PyCall
	using PyCall

	random = pyimport("random")
	os = pyimport("os")
	np = pyimport("numpy")
	pd = pyimport("pandas")
	scvi = pyimport("scvi")
	scvi_dataset = pyimport("scvi.dataset")
	scvi_models = pyimport("scvi.models")
	scvi_inference = pyimport("scvi.inference")
end;

# ╔═╡ d0caecf8-5f44-11eb-26fc-e749e38caaca
md"
##### Synthetic Single-Cell RNA-Sequencing Data from Small Pilot Studies using Deep Generative Models
"

# ╔═╡ db8ff29a-5f44-11eb-1178-076b084421f1
md"
###### Single-Cell Deep Boltzmann Machine (scDBM)
"

# ╔═╡ 0dc3d50e-5f45-11eb-1d35-957a2ad4984a
begin
	# Set hyperparameters
	epochsfirstlayer = 60
    learningratedbm = 1.0e-11
    dispinit = 0.5
    epochsdbm = 10
    nhiddennegbin = 4
    nhiddenbernoulli = 2
    lambdareg = 1.0
    repititions = collect(4:10:300)
end;

# ╔═╡ 641ff96c-5f4c-11eb-0a04-07bc290451a7
# Create mapping function
maplab(x1,x2,y) = map(arg -> y[findmin(sum((x1 .- x2[arg:arg,:]).^2,dims=2))[2][1]],1:size(x2,1));

# ╔═╡ 9c4dee34-5fb5-11eb-1711-3f606be0a97b
md"
Since the running time of the algorithms across all plates and iterations take a considerable amount of time, we only consider a subset of the overall analyses from the manuscript.
"

# ╔═╡ e0f4cab2-5f44-11eb-0471-8f2825d901cc
begin
	# Initialize arrays
	global dbi = DataFrame(Array{Float64,2}(undef, 1, 2));
	global dbi[:cells] = "Init"
	global ari = DataFrame(Array{Float64,2}(undef, 1, 1));
	global ari[:cells] = "Init"
	global cluster_prop = DataFrame(Array{Int64,2}(undef, 1, 4));
	global cluster_prop[:cells] = "Init"

	global Iter = ["384", "1152", "1920"]
	for t in Iter
		for i in 1:5
			# Countmatrix with 2000 highly variable genes
			Random.seed!(repititions[i])
			countmatrix = CSV.read("notebook_test/segerstolpe_hvg.csv")
			genenames = countmatrix[:,1]
			countmatrix = Array{Float64,2}(Array{Float64,2}(countmatrix[:,2:end])');

			# Sample t random cells
			cells = parse(Int64, t)
			global seedset = repititions[i]
			random_cells = rand(1:size(countmatrix,1),cells)

			# Save sub-sampled data
			countmatrix_tmp =
			DataFrame(hcat(collect(1:size(countmatrix[random_cells,:]',1)),
					countmatrix[random_cells,:]'))
			CSV.write(string("notebook_test/segerstolpe_cells",
					t,
					"_seed",
					seedset,
					".csv"),
				countmatrix_tmp,
				writeheader = true)

			# Train- test split
			data, datatest = splitdata(countmatrix[random_cells,:], 0.30);
			datadict = DataDict("Training data" => data, "Test data" => datatest);
			lr_rbm = [1.0e-11*ones(3500);0.00005*ones(1000);0.00001*ones(10000)];

			# Train scDBM
			global dbm = fitdbm(data,
				epochs = epochsdbm,
				learningrate = learningratedbm,
				batchsizepretraining = 64,
				pretraining = [
				# Negative-Binomial RBM
				TrainLayer(nhidden = nhiddennegbin,
						learningrates = lr_rbm,
						pcd = true,
						zeroinflation = false,
						epochs = epochsfirstlayer,
						rbmtype = NegativeBinomialBernoulliRBM,
						inversedispersion = (ones(size(data,2)) .* dispinit),
						fisherscoring = 1,
						estimatedispersion = "gene",
						lambda = lambdareg)
				# Bernoulli RBM
				TrainLayer(nhidden = nhiddenbernoulli,
						learningrate = 1.0e-11,
						epochs = 40)]
				)

			# Generate synthetic samples
			ngensamples = size(countmatrix,1);
			particles = initparticles(dbm, ngensamples);
			gibbssamplenegativebinomial!(countmatrix,
				particles,
				dbm, 50;
				zeroinflation=false);
			gen_data_dbm = particles[1];

			# Save generated data for one iteration
			if i == 3
				gen_data_dbm_out = DataFrame(gen_data_dbm)
				CSV.write(string("notebook_test/gen_data_scdbm",
					t,
					"_seed",
					seedset,
					".csv"),
				gen_data_dbm_out,
				writeheader = true)
			end

			# UMAP on synthetic data
			Random.seed!(123);
			umap_dbm = umap(gen_data_dbm',2;
				n_neighbors=30,
				metric=CosineDist(),
				min_dist=0.3)';

			# Load original cluster labels
			orig_labels = CSV.read("notebook_test/segerstolpe_hvg_clustering.csv");
			orig_labels = Array{Int64,1}(orig_labels[:x]);

			# UMAP on original data
			countmatrix = CSV.read("notebook_test/segerstolpe_hvg.csv");
			countmatrix = Array{Float64,2}(Array{Float64,2}(countmatrix[:,2:end])');
			Random.seed!(321);
			umap_original = umap(countmatrix',2;
				n_neighbors=30,
				metric=CosineDist(),
				min_dist=0.3)';

			# Mapping in UMAP space
			genlab_dbm = maplab(umap_original,umap_dbm,orig_labels);

			# DBI
			orig_dbindex = DBindex(countmatrix, orig_labels);
			dbm_dbindex = DBindex(gen_data_dbm, genlab_dbm);
			combined_dbi = DataFrame(vec([orig_dbindex, dbm_dbindex])')
        	combined_dbi[:cells] = t

			# ARI
			ari_index = DataFrame(
				hcat(Clustering.randindex(orig_labels,genlab_dbm)[2], t))
	        rename!(ari_index, :x2 => :cells);

			global dbi = vcat(dbi, combined_dbi)
			global ari = vcat(ari, ari_index)

			# Cluster Proportions
			countmap_dbm = countmap(genlab_dbm)
			countmap_orig = countmap(orig_labels)

			for l in 1:size(unique(orig_labels),1)
				if in(l, keys(countmap_dbm)) == true
					nothing
				else
					countmap_dbm[l] = 0
				end
			end

			tmp =	DataFrame(hcat(collect(values(countmap_orig)),
						collect(values(countmap_dbm)),
						collect(keys(countmap_orig)),
						parse.(Int64, repeat(["$i"], size(unique(orig_labels),1)))))
	        tmp[:cells] = t
    	    global cluster_prop = vcat(cluster_prop, tmp)
		end
	end
end;

# ╔═╡ 2c93f912-5f4c-11eb-3873-237144ded866
begin
	# DBI scDBM
	dbi_all = dbi[2:end,:]
	rename!(dbi_all, :x1 => :Original, :x2 => :scDBM);
	dbi_all = DataFrames.stack(dbi_all, 1:2)

	# ARI scDBM
	ari_all = ari[2:end,:]
	rename!(ari_all, :x1 => :ARI);

	# Cluster Proportions scDBM
	cluster_prop_all = cluster_prop[2:end,:]
	rename!(cluster_prop_all,
		:x1 => :Original,
		:x2 => :proportion,
		:x3 => :Cluster,
		:x4 => :Rep);
	cluster_prop_all[:model] = "scDBM"
end;

# ╔═╡ 71528b5c-5f4e-11eb-1c3e-7f62baec183a
md"
###### Single-Cell Variational Inference (scVI)
"

# ╔═╡ 87a19476-5fb6-11eb-36bf-7d0244187077
begin
	# Set fixed number of genes to avoid automatic sampling in scVI
	n_genes = 2000

	cd("notebook_test/");

	# Initialize arrays
	global dbi_scvi = DataFrame(Array{Float64,2}(undef, 1, 2));
	global dbi_scvi[:cells] = "Init";
	global ari_scvi = DataFrame(Array{Float64,2}(undef, 1, 1));
	global ari_scvi[:cells] = "Init"
	global cluster_prop_scvi = DataFrame(Array{Int64,2}(undef, 1, 4));
	global cluster_prop_scvi[:cells] = "Init"

	for t in Iter
    	for i = 1:5

			# Read sub-sampled data
			filenames = readdir()[occursin.(string("cells",t,"_seed"), readdir())];
			tmp = filenames[i];

			countmatrix_sample = scvi_dataset.CsvDataset(tmp,
					save_path = "",
					new_n_genes = n_genes);

			# Set hyperparameters
			n_latent=14
			n_hidden=64
			n_layers=2
			n_epochs=200
			lr_scvi=0.01
			use_batches = false
			use_cuda = true

			# Train the model and output model likelihood every epoch
			vae = scvi_models.VAE(countmatrix_sample.nb_genes,
				n_batch=countmatrix_sample.n_batches * use_batches,
				n_latent=n_latent,
				n_hidden=n_hidden,
				n_layers=n_layers,
				reconstruction_loss="nb",
				dispersion = "gene"
			);

			# Set up unsupervised trainer
			trainer = scvi_inference.UnsupervisedTrainer(vae,
				countmatrix_sample,
				train_size=0.7,
				use_cuda=use_cuda,
				frequency=1,
				n_epochs_kl_warmup=75,
				batch_size=36
        	);

			# Train scVI
			trainer.train(n_epochs=n_epochs,lr=lr_scvi)

			# Create posterior model
			full = trainer.create_posterior(trainer.model, countmatrix_sample);

			batch_size = 32
        	n_cells = 16

			# Generate Data
			gen_data_scvi = full.generate(batch_size=batch_size,
				n_samples=10,
				n_cells=n_cells)

			# Paste together generated cells
			n_subset = parse(Int64, t)

			n_paste = collect(1:trunc(Int,2090/n_subset))
			global gen_data_scvi_posterior = gen_data_scvi[1][:,:,1]
			for i in 2:size(n_paste,1)

				global gen_data_scvi_posterior = vcat(gen_data_scvi_posterior,
					gen_data_scvi[1][:,:,i])
			end
			rest = 2090 % n_subset

			gen_data_scvi = Array{Float64,2}(vcat(gen_data_scvi_posterior,
					gen_data_scvi[1][1:rest,:,size(n_paste,1) + 1]))

			# Save generated data for one iteration
			if i == 3
				gen_data_scvi_out = DataFrame(gen_data_scvi)
				CSV.write(string("gen_data_scvi",
					t,
					"_seed",
					seedset,
					".csv"),
				gen_data_scvi_out,
				writeheader = true)
			end

			# UMAP on synthetic data
			Random.seed!(123);
			umap_scvi = umap(gen_data_scvi',2;
				n_neighbors=30,
				metric=CosineDist(),
				min_dist=0.3)';

			# Load original cluster labels
			orig_labels = CSV.read("segerstolpe_hvg_clustering.csv")
			orig_labels = Array{Int64,1}(orig_labels[:x])

			# UMAP on original data
			countmatrix = CSV.read("segerstolpe_hvg.csv")
			countmatrix = Array{Float64,2}(Array{Float64,2}(countmatrix[:,2:end])')

			Random.seed!(321);
			umap_original = umap(countmatrix',2;
				n_neighbors=30,
				metric=CosineDist(),
				min_dist=0.3)';

			# Mapping in UMAP space
	        genlab_scvi = maplab(umap_original,umap_scvi,orig_labels);

			# DBI
			orig_dbindex_scvi = DBindex(countmatrix, orig_labels);
			scvi_dbindex = DBindex(gen_data_scvi, genlab_scvi);
			combined_dbi_scvi = DataFrame(vec([orig_dbindex_scvi, scvi_dbindex])')
        	combined_dbi_scvi[:cells] = t

			# ARI
			ari_index_scvi = DataFrame(
				hcat(Clustering.randindex(orig_labels,genlab_scvi)[2], t))
	        rename!(ari_index_scvi, :x2 => :cells);

			global dbi_scvi = vcat(dbi_scvi, combined_dbi_scvi)
			global ari_scvi = vcat(ari_scvi, ari_index_scvi)

			countmap_scvi = countmap(genlab_scvi)
			countmap_orig = countmap(orig_labels)
			for l in 1:size(unique(orig_labels),1)
				if in(l, keys(countmap_scvi)) == true
					nothing
				else
					countmap_scvi[l] = 0
				end
			end

			# Cluster Proportions
			global tmp1 = DataFrame(hcat(collect(values(countmap_orig)),
					collect(values(countmap_scvi)),
					collect(keys(countmap_orig)), parse.(Int64, repeat(["$i"],
							size(unique(orig_labels),1)))))
			global tmp1[:cells] = t
    	    global cluster_prop_scvi = vcat(cluster_prop_scvi, tmp1)
		end
	end
end;

# ╔═╡ bd161ef8-5fb6-11eb-2fb7-335fb40af173
begin
	# DBI scVI
	dbi_scvi_all = dbi_scvi[2:end,:];
	rename!(dbi_scvi_all, :x1 => :Original, :x2 => :scVI);
	dbi_scvi_all = DataFrames.stack(dbi_scvi_all, 1:2);

    # ARI scVI
    ari_scvi_all = ari_scvi[2:end,:];
    rename!(ari_scvi_all, :x1 => :ARI);

    # Cluster Proportions scVI
    cluster_prop_scvi_all = cluster_prop_scvi[2:end,:]
    rename!(cluster_prop_scvi_all,
		:x1 => :Original,
		:x2 => :proportion,
		:x3 => :Cluster,
		:x4 => :Rep);
	cluster_prop_scvi_all[:model] = "scVI"
end;

# ╔═╡ 945354fe-5fe6-11eb-3f0e-2df1f5961fac
md"
###### Combine results for both models
"

# ╔═╡ a242d864-5fe6-11eb-21f6-4b7640b2e816
begin
	dbi_plotting = vcat(dbi_all, dbi_scvi_all)
	ari_all[:model] = "scDBM"
	ari_scvi_all[:model] = "scVI"
	ari_plotting = vcat(ari_all, ari_scvi_all)
	cluster_prop_plotting = vcat(cluster_prop_all, cluster_prop_scvi_all)
end;

# ╔═╡ 961d9ff0-5fe7-11eb-14a4-5bb9ee53054a
begin
	# Save DBI data
	CSV.write(string("dbi_plotting.csv"),
				dbi_plotting,
				writeheader = true)

	# Save ARI data
	CSV.write(string("ari_plotting.csv"),
				ari_plotting,
				writeheader = true)

	# Save cluster proportion data
	CSV.write(string("cluster_prop_plotting.csv"),
				cluster_prop_plotting,
				writeheader = true)
end;

# ╔═╡ Cell order:
# ╟─d0caecf8-5f44-11eb-26fc-e749e38caaca
# ╠═d54abd5a-5f44-11eb-34fe-c96ecbd1327c
# ╠═db8ff29a-5f44-11eb-1178-076b084421f1
# ╠═59a1670e-5fb6-11eb-0cf1-6fefc45fb083
# ╠═0dc3d50e-5f45-11eb-1d35-957a2ad4984a
# ╠═641ff96c-5f4c-11eb-0a04-07bc290451a7
# ╟─9c4dee34-5fb5-11eb-1711-3f606be0a97b
# ╠═e0f4cab2-5f44-11eb-0471-8f2825d901cc
# ╠═2c93f912-5f4c-11eb-3873-237144ded866
# ╟─71528b5c-5f4e-11eb-1c3e-7f62baec183a
# ╠═09c53f44-5fb6-11eb-121b-757954f80b1d
# ╠═87a19476-5fb6-11eb-36bf-7d0244187077
# ╠═bd161ef8-5fb6-11eb-2fb7-335fb40af173
# ╟─945354fe-5fe6-11eb-3f0e-2df1f5961fac
# ╠═a242d864-5fe6-11eb-21f6-4b7640b2e816
# ╠═961d9ff0-5fe7-11eb-14a4-5bb9ee53054a

max_stars_repo_name: MTreppner/scDBM-paper

max_issues_repo_name: MTreppner/scDBM-paper

max_forks_repo_name: MTreppner/scDBM-paper


text:
function _parse_hp(block)
    dofs = OrderedDict{Symbol, PointStatus}()
    intros = :([])

    if block.head == :block
        for ex in block.args
            if ex isa Expr
                _parse_statement!(dofs, intros, ex)
            end
        end
    else
        _parse_statement!(dofs, intros, block)
    end
    return :(Hypothesis($dofs, $intros))
end

function _parse_th(block)
    th = :([])

    if block.head == :block
        for ex in block.args
            if ex isa Expr
                _parse_constraint!(th, ex)
            end
        end
    else
        if block isa Expr
            _parse_constraint!(th, block)
        end
    end

    return :(Thesis($th))
end