text:
#module Mnli

#using GeneralFunctions
#using OffsetArrays
#using TwoFAST

#export generateMnli

function generateMnli(Pk::Function, BiasDict22::Dict{String,F}, f::Real) where {F<:Real} #; q="best") 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)

    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)))*xi00^2)/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)))*xi20^2)/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)*xi40^2)/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)*xi20^2)/972405. - (32*(3*f*bη - 5*bΠ2para)*(7*bKKpara + 4*f^2*bη + 4*f^2*bη² +
            7*bΠ2para)*xi40^2)/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)*xi00^2)/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)*xi20^2)/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)*xi40^2)/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)))*xi20^2)/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))*xi40^2)/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)*xi20^2)/16807. + (256*f^2*(bη + bη²)*(3*f*bη -
            5*bΠ2para)*xi00*xi40)/77175. - (512*f^2*(bη + bη²)*(3*f*bη - 5*bΠ2para)*xi20*xi40)/151263. + (384*f^2*(bη +
            bη²)*(3*f*bη - 5*bΠ2para)*xi40^2)/420175.

    M[0,4,0] = @. (64*(3*f*bη - 5*bΠ2para)^2*xi00^2)/25725. - (128*(3*f*bη - 5*bΠ2para)^2*xi20^2)/36015. + (64*(3*f*bη -
            5*bΠ2para)^2*xi40^2)/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)*xi20^2)/756315. - (64*(3*f*bη - 5*bΠ2para)*(7*bKKpara + 4*f^2*bη + 4*f^2*bη² +
            7*bΠ2para)*xi40^2)/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))))*xi20^2)/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)))*xi00*xi40)/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))*xi20*xi40)/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))*xi40^2)/6.0085025e7

    M[0,6,4] = @. (-320*f^2*(bη + bη²)*(3*f*bη - 5*bΠ2para)*xi20^2)/26411. + (128*f^2*(bη + bη²)*(3*f*bη -
            5*bΠ2para)*xi00*xi40)/24255. - (1280*f^2*(bη + bη²)*(3*f*bη - 5*bΠ2para)*xi20*xi40)/237699. + (192*f^2*(bη +
            bη²)*(3*f*bη - 5*bΠ2para)*xi40^2)/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))*xi20*xi40)/4851. - (32*f^2*(bη + bη²)*(15*bKKpara + 8*f^2*bη + 8*f^2*bη² +
            15*bΠ2para)*xi40^2)/24255.

    M[0,8,8] = @. (128*f^4*(bη + bη²)^2*xi40^2)/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))*xi10*xi4m1)/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))*xi10*xi4m1)/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η²)*xi30*xi4m1)/6435.

    M[1,8,7] = @. (-256*f^4*(1 + bη)*(bη + bη²)*xi30*xi4m1)/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)*xi0m2*xi00)/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)))*xi1m1^2)/28350. + (4*f^3*(1 + bη)*(-7*b1 + 3*f + 3*f*bη)*xi2m2*xi20)/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)*xi3m1^2)/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)))*xi1m1^2)/496125. + (4*f^3*(1 +
            bη)*(-7*b1 + 3*f + 3*f*bη)*xi00*xi2m2)/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)))*xi1m1*xi3m1)/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)*xi3m1^2)/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)*xi0m2*xi00)/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)))*xi1m1^2)/99225. + (8*f^3*(1 + bη)*(-7*b1 + 3*f + 3*f*bη)*xi2m2*xi20)/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)*xi3m1^2)/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)))*xi1m1^2)/99225. + (22*f^3*(1 + bη)*(-7*b1 +
            3*f + 3*f*bη)*xi00*xi2m2)/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)))*xi1m1*xi3m1)/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)*xi3m1^2)/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))*xi1m1*xi3m1)/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)*xi3m1^2)/40425. + (16*f^4*(1 + bη)^2*xi00*xi4m2)/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)*xi1m1^2)/3675. - (16*f*(f - 3*bη)*(3*f*bη -
            5*bΠ2para)*xi3m1^2)/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))*xi1m1^2)/606375. + (8*f^3*(1 + bη)*(-7*b1 +
            3*f + 3*f*bη)*xi00*xi2m2)/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)*xi1m1*xi3m1)/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)*xi3m1^2)/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)*xi1m1*xi3m1)/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)*xi3m1^2)/148225. + (52*f^4*(1 + bη)^2*xi00*xi4m2)/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η²)*xi3m1^2)/11011. + (80*f^4*(1 + bη)^2*xi20*xi4m2)/11011.

    M[2,6,2] = @. (-256*f^2*(3*f*bη - 5*bΠ2para)*xi1m1^2)/40425. - (64*f^2*(3*f*bη - 5*bΠ2para)*xi1m1*xi3m1)/13475. +
    (64*f^2*(3*f*bη - 5*bΠ2para)*xi3m1^2)/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))*xi1m1*xi3m1)/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)*xi3m1^2)/12705. + (8*f^4*(1 + bη)^2*xi00*xi4m2)/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η²)*xi3m1^2)/38115. + (664*f^4*(1 + bη)^2*xi20*xi4m2)/38115.

    M[2,8,6] = @. (64*f^4*(1 + 4*bη + bη^2 + 2*bη²)*xi3m1^2)/6435. + (64*f^4*(1 + bη)^2*xi20*xi4m2)/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)*xi0m1*xi1m2)/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)*xi1m2*xi2m1)/1125. -
    (4*f^2*(1 + bη)*(35*b1 + 9*f^2 - 21*f*bη)*xi2m1*xi3m2)/2625.

    M[3,0,3] = @. (4*f^3*(5*b1 - 3*f + (2 - 3*f)*bη + 2*bη^2)*xi1m2*xi2m1)/875. - (4*f^3*(f - bη)*(1 +
            bη)*xi0m1*xi3m2)/525. + (16*f^3*(2*f - 3*bη)*(1 + bη)*xi2m1*xi3m2)/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)*xi0m1*xi1m2)/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)*xi1m2*xi2m1)/1575. - (8*f^2*(1 + bη)*(49*b1 + 18*f^2 -
            42*f*bη)*xi2m1*xi3m2)/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)*xi1m2*xi2m1)/525. -
    (2*f^2*(1 + bη)*(9*b1 + 7*f^2 - 7*f*bη)*xi0m1*xi3m2)/315. + (8*f^2*(1 + bη)*(27*b1 + 14*f^2 -
            21*f*bη)*xi2m1*xi3m2)/4725.

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

    M[3,4,1] = @. (-8*f^2*(f - bη)*(-5*b1 + 3*f + 3*f*bη)*xi0m1*xi1m2)/525. + (16*f^2*(2*f*(-5*b1 + 3*f) + (25*b1 +
                3*f*(-3 + 2*f))*bη - 9*f*bη^2)*xi1m2*xi2m1)/2625. - (32*f^3*(3*f - 7*bη)*(1 + bη)*xi2m1*xi3m2)/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)*xi1m2*xi2m1)/9625. -
    (8*f^2*(1 + bη)*(55*b1 + 42*f^2 - 42*f*bη)*xi0m1*xi3m2)/5775. + (32*f^2*(1 + bη)*(55*b1 + 28*f^2 -
            42*f*bη)*xi2m1*xi3m2)/28875.

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

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

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

    M[3,8,5] = @. (-128*f^4*(1 + bη)*xi2m1*xi3m2)/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)*xi0m2^2)/540. + ((225*b1^2 +
            50*b1*f^2 + 6*f^4 - 30*(5*b1*f + f^3)*bη + 45*f^2*bη^2)*xi2m2^2)/1350.

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

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

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

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

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

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

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

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

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

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

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

   return r, M

end

function generateMnlismall(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)

    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)*xi20^2)/972405. - (32*(3*f*bη - 5*bΠ2para)*(7*bKKpara + 4*f^2*bη + 4*f^2*bη² +
            7*bΠ2para)*xi40^2)/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.)

   return r, M

end

function generateMnliold(Pk::Function, BiasDict22::Dict{String,F}, f::Real; q="best") 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∥"]
    @show b1
    @show bη
    @show b2
    @show bK²
    @show bδη
    @show bη²
    @show bKKpara
    @show bΠ2para

    #@time begin
    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)
    #end

    M[0,0,0] = @. ((40*b1^2)/147. + (10*b1*b2)/21. + b2^2/4. - (20*b1^2*f)/63. - (b1*b2*f)/3. + (7*b1^2*f^2)/45. -
        (8*b1*bK²)/63. + (4*bK²^2)/45. - (8*b1*bKKpara)/189. + (8*bK²*bKKpara)/135. + (34*bKKpara^2)/2025. -
        (20*b1*f*bδη)/63. - (b2*f*bδη)/3. + (14*b1*f^2*bδη)/45. + (7*f^2*bδη^2)/45. - (16*b1*f*bη)/147. -
        (2*b2*f*bη)/21. + (16*b1*f^2*bη)/105. + (b2*f^2*bη)/9. - (2*b1*f^3*bη)/15. + (8*f*bK²*bη)/315. +
        (16*f^2*bK²*bη)/675. + (8*f*bKKpara*bη)/945. + (32*f^2*bKKpara*bη)/2025. + (4*f^2*bδη*bη)/63. -
        (2*f^3*bδη*bη)/15. + (24*f^2*bη^2)/1225. - (4*f^3*bη^2)/225. + (f^4*bη^2)/25. + (4*b1*f^2*bη²)/45. +
        (b2*f^2*bη²)/9. - (2*b1*f^3*bη²)/15. + (16*f^2*bK²*bη²)/675. + (32*f^2*bKKpara*bη²)/2025. -
        (2*f^3*bδη*bη²)/15. - (4*f^3*bη*bη²)/225. + (2*f^4*bη*bη²)/25. + (f^4*bη²^2)/25. -
        ((-8100*b1 - 8925*b2 + 7910*b1*f - 560*bK² - 644*bKKpara + 7910*f*bδη + 2580*f*bη - 3626*f^2*bη -
                3626*f^2*bη²)*bΠ2para)/33075. + (4117*bΠ2para^2)/33075.)*xi00^2 + ((-400*b1^2)/1029. -
        (10*b1*b2)/21. + (20*b1^2*f)/63. + (2*b1^2*f^2)/45. + (200*b1*bK²)/441. + (2*b2*bK²)/3. -
        (4*b1*f*bK²)/9. + (8*bK²^2)/63. + (200*b1*bKKpara)/1323. + (2*b2*bKKpara)/9. - (28*b1*f*bKKpara)/135. +
        (16*bK²*bKKpara)/189. + (34*bKKpara^2)/2835. + (20*b1*f*bδη)/63. + (4*b1*f^2*bδη)/45. -
        (4*f*bK²*bδη)/9. - (28*f*bKKpara*bδη)/135. + (2*f^2*bδη^2)/45. + (160*b1*f*bη)/1029. +
        (2*b2*f*bη)/21. - (16*b1*f^2*bη)/147. + (4*b2*f^2*bη)/45. - (16*b1*f^3*bη)/105. - (40*f*bK²*bη)/441. +
        (172*f^2*bK²*bη)/945. - (40*f*bKKpara*bη)/1323. + (388*f^2*bKKpara*bη)/3969. - (4*f^2*bδη*bη)/63. -
        (16*f^3*bδη*bη)/105. - (48*f^2*bη^2)/1715. + (4*f^3*bη^2)/441. + (16*f^4*bη^2)/245. -
        (20*b1*f^2*bη²)/441. + (4*b2*f^2*bη²)/45. - (16*b1*f^3*bη²)/105. + (172*f^2*bK²*bη²)/945. +
        (388*f^2*bKKpara*bη²)/3969. - (16*f^3*bδη*bη²)/105. + (4*f^3*bη*bη²)/441. + (32*f^4*bη*bη²)/245. +
        (16*f^4*bη²^2)/245. - (4*(2475*b1 - 735*b2 + 1862*b1*f - 4445*bK² - 1890*bKKpara + 1862*f*bδη -
                975*f*bη - 1904*f^2*bη - 1904*f^2*bη²)*bΠ2para)/46305. + (2686*bΠ2para^2)/46305.)*xi20^2 +
    ((8*(1125*b1^2 - 3150*b1*bK² + 2205*bK²^2 - 1050*b1*bKKpara + 1470*bK²*bKKpara + 294*bKKpara^2 - 450*b1*f*bη -
                420*b1*f^2*bη + 630*f*bK²*bη + 588*f^2*bK²*bη + 210*f*bKKpara*bη + 252*f^2*bKKpara*bη +
                81*f^2*bη^2 + 84*f^3*bη^2 + 56*f^4*bη^2 - 420*b1*f^2*bη² + 588*f^2*bK²*bη² + 252*f^2*bKKpara*bη² +
                84*f^3*bη*bη² + 112*f^4*bη*bη² + 56*f^4*bη²^2))/77175. - (16*(150*b1 - 210*bK² - 119*bKKpara +
                30*f*bη - 56*f^2*bη - 56*f^2*bη²)*bΠ2para)/77175. + (1352*bΠ2para^2)/77175.)*xi40^2

    M[0,0,2] = @. ((32*f*bη*(147*b1*f + 182*bKKpara + 147*f*bδη + 6*f^2*bη + 6*f^2*bη²))/324135. -
        (32*(735*b1*f + 910*bKKpara + 735*f*bδη - 252*f*bη + 30*f^2*bη + 30*f^2*bη²)*bΠ2para)/972405. -
        (128*bΠ2para^2)/9261.)*xi20^2 + ((-32*f*bη*(7*bKKpara + 4*f^2*bη + 4*f^2*bη²))/36015. +
        (32*(35*bKKpara - 21*f*bη + 20*f^2*bη + 20*f^2*bη²)*bΠ2para)/108045. + (32*bΠ2para^2)/3087.)*xi40^2 +
    xi20*(((32*f*bη*(42*b1*f + 7*bKKpara + 42*f*bδη - 24*f^2*bη - 24*f^2*bη²))/33075. -
            (32*(210*b1*f + 35*bKKpara + 210*f*bδη + 63*f*bη - 120*f^2*bη - 120*f^2*bη²)*bΠ2para)/99225. +
            (32*bΠ2para^2)/945.)*xi00 + ((-32*f*bη*(147*b1*f - 28*bKKpara + 147*f*bδη -
                    114*f^2*bη - 114*f^2*bη²))/180075. + (32*(735*b1*f - 140*bKKpara + 735*f*bδη + 378*f*bη -
                    570*f^2*bη - 570*f^2*bη²)*bΠ2para)/540225. - (64*bΠ2para^2)/1715.)*xi40)

    M[0,2,0] = @. ((4*f*bη*(-12600*b1 - 11025*b2 + 7350*b1*f + 2940*bK² + 980*bKKpara + 7350*f*bδη + 3240*f*bη -
                2058*f^2*bη - 2058*f^2*bη²))/231525. - (4*(-12600*b1 - 11025*b2 + 7350*b1*f + 2940*bK² +
                980*bKKpara + 7350*f*bδη + 7362*f*bη - 2058*f^2*bη - 2058*f^2*bη²)*bΠ2para)/138915. +
        (1832*bΠ2para^2)/9261.)*xi00^2 + ((-4*f*bη*(-5040*b1 - 3087*b2 + 2058*b1*f + 2940*bK² + 980*bKKpara +
                2058*f*bδη + 1296*f*bη - 294*f^2*bη - 294*f^2*bη²))/64827. + (4*(-25200*b1 - 15435*b2 +
                10290*b1*f + 14700*bK² + 4900*bKKpara + 10290*f*bδη + 12078*f*bη - 1470*f^2*bη -
                1470*f^2*bη²)*bΠ2para)/194481. - (12440*bΠ2para^2)/64827.)*xi20^2 +
    ((-32*f*bη*(525*b1 - 735*bK² - 245*bKKpara - 135*f*bη - 98*f^2*bη - 98*f^2*bη²))/180075. +
        (32*(525*b1 - 735*bK² - 245*bKKpara - 123*f*bη - 98*f^2*bη - 98*f^2*bη²)*bΠ2para)/108045. -
        (128*bΠ2para^2)/21609.)*xi40^2

    M[0,2,2] = @. ((-80*b1^2*f)/441. - (4*b1^2*f^2)/63. + (16*b1*f*bK²)/63. - (2080*b1*bKKpara)/9261. -
        (2*b2*bKKpara)/9. + (284*b1*f*bKKpara)/1323. + (8*bK²*bKKpara)/441. + (4*bKKpara^2)/27783. -
        (80*b1*f*bδη)/441. - (8*b1*f^2*bδη)/63. + (16*f*bK²*bδη)/63. + (284*f*bKKpara*bδη)/1323. -
        (4*f^2*bδη^2)/63. + (1072*b1*f^2*bη)/21609. - (8*b2*f^2*bη)/63. + (32*b1*f^3*bη)/147. -
        (1472*f^2*bK²*bη)/9261. + (4576*f*bKKpara*bη)/64827. - (2908*f^2*bKKpara*bη)/27783. +
        (176*f^2*bδη*bη)/3087. + (32*f^3*bδη*bη)/147. + (352*f^3*bη^2)/151263. - (32*f^4*bη^2)/343. -
        (160*b1*f^2*bη²)/21609. - (8*b2*f^2*bη²)/63. + (32*b1*f^3*bη²)/147. - (1472*f^2*bK²*bη²)/9261. -
        (2908*f^2*bKKpara*bη²)/27783. + (32*f^3*bδη*bη²)/147. + (352*f^3*bη*bη²)/151263. - (64*f^4*bη*bη²)/343. -
        (32*f^4*bη²^2)/343. + (2*(-23520*b1 - 50421*b2 + 46354*b1*f - 34300*bK² - 47894*bKKpara + 46354*f*bδη +
                7392*f*bη - 51154*f^2*bη - 51154*f^2*bη²)*bΠ2para)/453789. - (11416*bΠ2para^2)/64827.)*xi20^2 +
    ((16*(8085*b1*bKKpara - 11319*bK²*bKKpara - 3234*bKKpara^2 + 4620*b1*f^2*bη - 6468*f^2*bK²*bη -
                2541*f*bKKpara*bη - 3234*f^2*bKKpara*bη - 1452*f^3*bη^2 - 784*f^4*bη^2 + 4620*b1*f^2*bη² -
                6468*f^2*bK²*bη² - 3234*f^2*bKKpara*bη² - 1452*f^3*bη*bη² - 1568*f^4*bη*bη² - 784*f^4*bη²^2))/1.663893e6 +
        (16*(735*b1 - 1029*bK² - 203*bKKpara - 231*f*bη - 74*f^2*bη - 74*f^2*bη²)*bΠ2para)/151263. +
        (208*bΠ2para^2)/21609.)*xi40^2 + xi20*(((-32*b1^2*f)/63. - (2*b1*b2*f)/3. + (32*b1^2*f^2)/63. -
            (8*b1*f*bK²)/45. - (16*b1*bKKpara)/189. - (16*b1*f*bKKpara)/135. + (16*bK²*bKKpara)/135. +
            (20*bKKpara^2)/567. - (32*b1*f*bδη)/63. - (2*b2*f*bδη)/3. + (64*b1*f^2*bδη)/63. - (8*f*bK²*bδη)/45. -
            (16*f*bKKpara*bδη)/135. + (32*f^2*bδη^2)/63. + (992*b1*f^2*bη)/2205. + (4*b2*f^2*bη)/9. -
            (68*b1*f^3*bη)/105. + (176*f^2*bK²*bη)/945. + (176*f*bKKpara*bη)/6615. + (352*f^2*bKKpara*bη)/2835. +
            (352*f^2*bδη*bη)/2205. - (68*f^3*bδη*bη)/105. - (1408*f^3*bη^2)/15435. + (8*f^4*bη^2)/35. +
            (128*b1*f^2*bη²)/441. + (4*b2*f^2*bη²)/9. - (68*b1*f^3*bη²)/105. + (176*f^2*bK²*bη²)/945. +
            (352*f^2*bKKpara*bη²)/2835. - (68*f^3*bδη*bη²)/105. - (1408*f^3*bη*bη²)/15435. + (16*f^4*bη*bη²)/35. +
            (8*f^4*bη²^2)/35. - (4*(-2940*b1 - 5145*b2 + 10577*b1*f - 2744*bK² - 1218*bKKpara + 10577*f*bδη +
                    924*f*bη - 7052*f^2*bη - 7052*f^2*bη²)*bΠ2para)/46305. + (2476*bΠ2para^2)/6615.)*xi00 +
        ((16*b1^2*f)/49. - (16*b1*f*bK²)/35. - (64*b1*bKKpara)/1029. - (48*b1*f*bKKpara)/245. + (64*bK²*bKKpara)/735. +
            (16*bKKpara^2)/343. + (16*b1*f*bδη)/49. - (16*f*bK²*bδη)/35. - (48*f*bKKpara*bδη)/245. -
            (4272*b1*f^2*bη)/12005. - (64*b1*f^3*bη)/735. + (608*f^2*bK²*bη)/1715. + (704*f*bKKpara*bη)/36015. +
            (2848*f^2*bKKpara*bη)/15435. - (176*f^2*bδη*bη)/1715. - (64*f^3*bδη*bη)/735. + (6688*f^3*bη^2)/84035. +
            (128*f^4*bη^2)/1715. - (608*b1*f^2*bη²)/2401. - (64*b1*f^3*bη²)/735. + (608*f^2*bK²*bη²)/1715. +
            (2848*f^2*bKKpara*bη²)/15435. - (64*f^3*bδη*bη²)/735. + (6688*f^3*bη*bη²)/84035. + (256*f^4*bη*bη²)/1715. +
            (128*f^4*bη²^2)/1715. - (32*(6615*b1 + 588*b1*f - 9261*bK² - 4522*bKKpara + 588*f*bδη - 2079*f*bη -
                    2598*f^2*bη - 2598*f^2*bη²)*bΠ2para)/756315. + (368*bΠ2para^2)/12005.)*xi40)

    M[0,2,4] = @. ((-384*f^3*bη*(bη + bη²))/16807. + (640*f^2*(bη + bη²)*bΠ2para)/16807.)*xi20^2 +
    ((256*f^3*bη*(bη + bη²))/25725. - (256*f^2*(bη + bη²)*bΠ2para)/15435.)*xi00*xi40 +
    ((-512*f^3*bη*(bη + bη²))/50421. + (2560*f^2*(bη + bη²)*bΠ2para)/151263.)*xi20*xi40 +
    ((1152*f^3*bη*(bη + bη²))/420175. - (384*f^2*(bη + bη²)*bΠ2para)/84035.)*xi40^2

    M[0,4,0] = @. ((192*f^2*bη^2)/8575. - (128*f*bη*bΠ2para)/1715. + (64*bΠ2para^2)/1029.)*xi00^2 +
    ((-384*f^2*bη^2)/12005. + (256*f*bη*bΠ2para)/2401. - (640*bΠ2para^2)/7203.)*xi20^2 +
    ((576*f^2*bη^2)/60025. - (384*f*bη*bΠ2para)/12005. + (64*bΠ2para^2)/2401.)*xi40^2

    M[0,4,2] = @. ((64*f*bη*(147*b1*f + 182*bKKpara + 147*f*bδη + 6*f^2*bη + 6*f^2*bη²))/252105. -
        (64*(735*b1*f + 910*bKKpara + 735*f*bδη - 252*f*bη + 30*f^2*bη + 30*f^2*bη²)*bΠ2para)/756315. -
        (256*bΠ2para^2)/7203.)*xi20^2 + ((-192*f*bη*(7*bKKpara + 4*f^2*bη + 4*f^2*bη²))/84035. +
        (64*(35*bKKpara - 21*f*bη + 20*f^2*bη + 20*f^2*bη²)*bΠ2para)/84035. + (64*bΠ2para^2)/2401.)*xi40^2 +
    xi20*(((64*f*bη*(42*b1*f + 7*bKKpara + 42*f*bδη - 24*f^2*bη - 24*f^2*bη²))/25725. -
            (64*(210*b1*f + 35*bKKpara + 210*f*bδη + 63*f*bη - 120*f^2*bη - 120*f^2*bη²)*bΠ2para)/77175. +
            (64*bΠ2para^2)/735.)*xi00 + ((-192*f*bη*(147*b1*f - 28*bKKpara + 147*f*bδη - 114*f^2*bη -
                    114*f^2*bη²))/420175. + (64*(735*b1*f - 140*bKKpara + 735*f*bδη + 378*f*bη - 570*f^2*bη -
                    570*f^2*bη²)*bΠ2para)/420175. - (1152*bΠ2para^2)/12005.)*xi40)

    M[0,4,4] = @. ((4*b1^2*f^2)/35. - (32*b1*f*bKKpara)/735. + (136*bKKpara^2)/15435. + (8*b1*f^2*bδη)/35. -
        (32*f*bKKpara*bδη)/735. + (4*f^2*bδη^2)/35. + (480*b1*f^2*bη)/2401. + (8*b2*f^2*bη)/35. -
        (96*b1*f^3*bη)/245. + (128*f^2*bK²*bη)/5145. + (128*f^2*bKKpara*bη)/3087. - (96*f^3*bδη*bη)/245. -
        (11232*f^3*bη^2)/184877. + (288*f^4*bη^2)/1715. + (480*b1*f^2*bη²)/2401. + (8*b2*f^2*bη²)/35. -
        (96*b1*f^3*bη²)/245. + (128*f^2*bK²*bη²)/5145. + (128*f^2*bKKpara*bη²)/3087. - (96*f^3*bδη*bη²)/245. -
        (11232*f^3*bη*bη²)/184877. + (576*f^4*bη*bη²)/1715. + (288*f^4*bη²^2)/1715. - (16*(11319*b1*f -
                2695*bKKpara + 11319*f*bδη - 20403*f^2*bη - 20403*f^2*bη²)*bΠ2para)/924385. +
        (152*bΠ2para^2)/1715.)*xi20^2 + ((4*(24255*b1^2*f^2 + 2156*bKKpara^2 + 48510*b1*f^2*bδη +
                24255*f^2*bδη^2 - 18480*b1*f^2*bη - 32340*b1*f^3*bη + 25872*f^2*bK²*bη + 17248*f^2*bKKpara*bη -
                32340*f^3*bδη*bη + 5616*f^3*bη^2 + 19404*f^4*bη^2 - 18480*b1*f^2*bη² - 32340*b1*f^3*bη² +
                25872*f^2*bK²*bη² + 17248*f^2*bKKpara*bη² - 32340*f^3*bδη*bη² + 5616*f^3*bη*bη² + 38808*f^4*bη*bη² +
                19404*f^4*bη²^2))/848925. - (16*(8085*b1*f - 1078*bKKpara + 8085*f*bδη - 7362*f^2*bη -
                7362*f^2*bη²)*bΠ2para)/848925. + (32*bΠ2para^2)/525.)*xi00*xi40 + ((16*(11319*b1*f*bKKpara +
                1078*bKKpara^2 + 11319*f*bKKpara*bδη + 9240*b1*f^2*bη + 8232*b1*f^3*bη - 12936*f^2*bK²*bη -
                11270*f^2*bKKpara*bη + 8232*f^3*bδη*bη - 2808*f^3*bη^2 - 7056*f^4*bη^2 + 9240*b1*f^2*bη² +
                8232*b1*f^3*bη² - 12936*f^2*bK²*bη² - 11270*f^2*bKKpara*bη² + 8232*f^3*bδη*bη² - 2808*f^3*bη*bη² -
                14112*f^4*bη*bη² - 7056*f^4*bη²^2))/1.663893e6 + (16*(1029*b1*f - 490*bKKpara + 1029*f*bδη -
                1098*f^2*bη - 1098*f^2*bη²)*bΠ2para)/151263. - (64*bΠ2para^2)/1029.)*xi20*xi40 +
    ((144*(7007*bKKpara^2 - 10010*b1*f^2*bη + 14014*f^2*bK²*bη + 11466*f^2*bKKpara*bη + 3042*f^3*bη^2 +
                3528*f^4*bη^2 - 10010*b1*f^2*bη² + 14014*f^2*bK²*bη² + 11466*f^2*bKKpara*bη² + 3042*f^3*bη*bη² +
                7056*f^4*bη*bη² + 3528*f^4*bη²^2))/6.0085025e7 + (288*(539*bKKpara + 246*f^2*bη +
                246*f^2*bη²)*bΠ2para)/4.621925e6 + (144*bΠ2para^2)/8575.)*xi40^2

    M[0,6,4] = @. ((-960*f^3*bη*(bη + bη²))/26411. + (1600*f^2*(bη + bη²)*bΠ2para)/26411.)*xi20^2 +
    ((128*f^3*bη*(bη + bη²))/8085. - (128*f^2*(bη + bη²)*bΠ2para)/4851.)*xi00*xi40 +
    ((-1280*f^3*bη*(bη + bη²))/79233. + (6400*f^2*(bη + bη²)*bΠ2para)/237699.)*xi20*xi40 +
    ((576*f^3*bη*(bη + bη²))/132055. - (192*f^2*(bη + bη²)*bΠ2para)/26411.)*xi40^2

    M[0,6,6] = @. ((-32*f^2*(bη + bη²)*(21*b1*f - 4*bKKpara + 21*f*bδη - 18*f^2*bη - 18*f^2*bη²))/4851. +
        (64*f^2*(bη + bη²)*bΠ2para)/539.)*xi20*xi40 + ((-32*f^2*(bη + bη²)*(15*bKKpara + 8*f^2*bη + 8*f^2*bη²))/24255. -
        (32*f^2*(bη + bη²)*bΠ2para)/1617.)*xi40^2

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

    M[1,0,1] = @. xi10*(((20*b1^2*f)/63. + (b1*b2*f)/3. - (4*b1*f^2)/45. - (14*b1^2*f^2)/45. - (b2*f^2)/9. +
            (2*b1*f^3)/15. - (16*f^2*bK²)/675. - (32*f^2*bKKpara)/2025. - (14*b1*f^2*bδη)/45. + (2*f^3*bδη)/15. -
            (64*b1*f^2*bη)/315. - (b2*f^2*bη)/9. + (4*f^3*bη)/125. + (4*b1*f^3*bη)/15. - (2*f^4*bη)/25. -
            (16*f^2*bK²*bη)/675. - (32*f^2*bKKpara*bη)/2025. + (2*f^3*bδη*bη)/15. + (4*f^3*bη^2)/125. -
            (2*f^4*bη^2)/25. + (2*b1*f^3*bη²)/15. - (2*f^4*bη²)/25. - (2*f^4*bη*bη²)/25. +
            (2*f*(17*b1 - 7*f - 7*f*bη)*bΠ2para)/105.)*xi0m1 + ((8*b1^2*f)/63. + (8*b1*f^2)/315. +
            (4*b1^2*f^2)/45. + (4*b2*f^2)/45. - (68*b1*f^3)/525. - (8*b1*f*bK²)/45. + (56*f^2*bK²)/675. -
            (68*b1*f*bKKpara)/675. + (676*f^2*bKKpara)/14175. + (4*b1*f^2*bδη)/45. - (68*f^3*bδη)/525. -
            (32*b1*f^2*bη)/1575. + (4*b2*f^2*bη)/45. - (8*f^3*bη)/875. - (124*b1*f^3*bη)/525. + (16*f^4*bη)/175. +
            (56*f^2*bK²*bη)/675. + (676*f^2*bKKpara*bη)/14175. - (68*f^3*bδη*bη)/525. - (8*f^3*bη^2)/875. +
            (16*f^4*bη^2)/175. - (8*b1*f^3*bη²)/75. + (16*f^4*bη²)/175. + (16*f^4*bη*bη²)/175. -
            (44*b1*f*bΠ2para)/525. + (68*f^2*bΠ2para)/525. + (68*f^2*bη*bΠ2para)/525.)*xi2m1) +
    xi30*((-4*b1^2*f*xi2m1)/21. + (52*b1*f^2*xi2m1)/735. + (4*b1*f^3*xi2m1)/175. + (4*b1*f*bK²*xi2m1)/15. -
        (52*f^2*bK²*xi2m1)/525. + (8*b1*f*bKKpara*xi2m1)/75. - (184*f^2*bKKpara*xi2m1)/3675. + (4*f^3*bδη*xi2m1)/175. +
        (512*b1*f^2*bη*xi2m1)/3675. - (156*f^3*bη*xi2m1)/6125. + (12*b1*f^3*bη*xi2m1)/175. - (48*f^4*bη*xi2m1)/1225. -
        (52*f^2*bK²*bη*xi2m1)/525. - (184*f^2*bKKpara*bη*xi2m1)/3675. + (4*f^3*bδη*bη*xi2m1)/175. -
        (156*f^3*bη^2*xi2m1)/6125. - (48*f^4*bη^2*xi2m1)/1225. + (8*b1*f^3*bη²*xi2m1)/175. - (48*f^4*bη²*xi2m1)/1225. -
        (48*f^4*bη*bη²*xi2m1)/1225. - (4*b1*f*bΠ2para*xi2m1)/525. - (4*f^2*bΠ2para*xi2m1)/175. -
        (4*f^2*bη*bΠ2para*xi2m1)/175. - (32*f^2*(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^3*bη*(1 + bη))/6125. + (16*f^2*(1 + bη)*bΠ2para)/3675.)*xi0m1 +
        (128*f^2*(1 + bη)*(3*f*bη - 5*bΠ2para)*xi2m1)/128625. - (32*f^2*(1 + bη)*(3*f*bη - 5*bΠ2para)*xi4m1)/42875.) +
    xi10*((176*f^2*(1 + bη)*(3*f*bη - 5*bΠ2para)*xi2m1)/42875. + (64*f^2*(1 + bη)*(3*f*bη - 5*bΠ2para)*xi4m1)/42875.)

    M[1,2,1] = @. xi10*(((40*b1^2*f)/63. + (2*b1*b2*f)/3. - (8*b1*f^2)/45. - (28*b1^2*f^2)/45. - (2*b2*f^2)/9. +
            (4*b1*f^3)/15. - (32*f^2*bK²)/675. - (64*f^2*bKKpara)/2025. - (28*b1*f^2*bδη)/45. + (4*f^3*bδη)/15. -
            (1112*b1*f^2*bη)/2205. - (2*b2*f^2*bη)/9. + (16*f^3*bη)/175. + (8*b1*f^3*bη)/15. - (4*f^4*bη)/25. -
            (32*f^2*bK²*bη)/675. - (64*f^2*bKKpara*bη)/2025. + (4*f^3*bδη*bη)/15. + (16*f^3*bη^2)/175. -
            (4*f^4*bη^2)/25. + (4*b1*f^3*bη²)/15. - (4*f^4*bη²)/25. - (4*f^4*bη*bη²)/25. +
            (4*f*(745*b1 - 287*f - 287*f*bη)*bΠ2para)/3675.)*xi0m1 + ((16*b1^2*f)/63. + (16*b1*f^2)/315. +
            (8*b1^2*f^2)/45. + (8*b2*f^2)/45. - (136*b1*f^3)/525. - (16*b1*f*bK²)/45. + (112*f^2*bK²)/675. -
            (136*b1*f*bKKpara)/675. + (1352*f^2*bKKpara)/14175. + (8*b1*f^2*bδη)/45. - (136*f^3*bδη)/525. -
            (176*b1*f^2*bη)/2205. + (8*b2*f^2*bη)/45. - (32*f^3*bη)/1225. - (248*b1*f^3*bη)/525. + (32*f^4*bη)/175. +
            (112*f^2*bK²*bη)/675. + (1352*f^2*bKKpara*bη)/14175. - (136*f^3*bδη*bη)/525. - (32*f^3*bη^2)/1225. +
            (32*f^4*bη^2)/175. - (16*b1*f^3*bη²)/75. + (32*f^4*bη²)/175. + (32*f^4*bη*bη²)/175. -
            (376*b1*f*bΠ2para)/3675. + (40*f^2*bΠ2para)/147. + (40*f^2*bη*bΠ2para)/147.)*xi2m1) +
    xi30*((-8*b1^2*f*xi2m1)/21. + (104*b1*f^2*xi2m1)/735. + (8*b1*f^3*xi2m1)/175. + (8*b1*f*bK²*xi2m1)/15. -
        (104*f^2*bK²*xi2m1)/525. + (16*b1*f*bKKpara*xi2m1)/75. - (368*f^2*bKKpara*xi2m1)/3675. +
        (8*f^3*bδη*xi2m1)/175. + (248*b1*f^2*bη*xi2m1)/735. - (624*f^3*bη*xi2m1)/8575. + (24*b1*f^3*bη*xi2m1)/175. -
        (96*f^4*bη*xi2m1)/1225. - (104*f^2*bK²*bη*xi2m1)/525. - (368*f^2*bKKpara*bη*xi2m1)/3675. +
        (8*f^3*bδη*bη*xi2m1)/175. - (624*f^3*bη^2*xi2m1)/8575. - (96*f^4*bη^2*xi2m1)/1225. + (16*b1*f^3*bη²*xi2m1)/175. -
        (96*f^4*bη²*xi2m1)/1225. - (96*f^4*bη*bη²*xi2m1)/1225. - (416*b1*f*bΠ2para*xi2m1)/3675. -
        (16*f^2*bΠ2para*xi2m1)/1715. - (16*f^2*bη*bΠ2para*xi2m1)/1715. - (64*f^2*(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] = @. xi10*((-264*b1*f^2*xi2m1)/1715. - (6*b1^2*f^2*xi2m1)/35. - (6*b2*f^2*xi2m1)/35. +
        (306*b1*f^3*xi2m1)/1225. - (16*f^2*bK²*xi2m1)/1225. + (4*b1*f*bKKpara*xi2m1)/175. -
        (68*f^2*bKKpara*xi2m1)/3675. - (6*b1*f^2*bδη*xi2m1)/35. + (306*f^3*bδη*xi2m1)/1225. -
        (264*b1*f^2*bη*xi2m1)/1715. - (6*b2*f^2*bη*xi2m1)/35. + (88*f^3*bη*xi2m1)/1225. + (558*b1*f^3*bη*xi2m1)/1225. -
        (216*f^4*bη*xi2m1)/1225. - (16*f^2*bK²*bη*xi2m1)/1225. - (68*f^2*bKKpara*bη*xi2m1)/3675. +
        (306*f^3*bδη*bη*xi2m1)/1225. + (88*f^3*bη^2*xi2m1)/1225. - (216*f^4*bη^2*xi2m1)/1225. +
        (36*b1*f^3*bη²*xi2m1)/175. - (216*f^4*bη²*xi2m1)/1225. - (216*f^4*bη*bη²*xi2m1)/1225. +
        (24*b1*f*bΠ2para*xi2m1)/175. - (4*f^2*bΠ2para*xi2m1)/15. - (4*f^2*bη*bΠ2para*xi2m1)/15. -
        (16*f^2*(1 + bη)*(270*b1 + 210*b1*f - 378*bK² - 217*bKKpara + 210*f*bδη - 126*f*bη - 252*f^2*bη - 252*f^2*bη² -
                147*bΠ2para)*xi4m1)/77175.) + xi30*(((-2*f^2*(-60*b1 + 315*b1^2 - 105*b1*f + 84*bK² + 56*bKKpara +
                    315*b1*bδη - 105*f*bδη - 60*b1*bη + 28*f*bη - 315*b1*f*bη + 126*f^2*bη + 84*bK²*bη + 56*bKKpara*bη -
                    105*f*bδη*bη + 28*f*bη^2 + 126*f^2*bη^2 - 210*b1*f*bη² + 126*f^2*bη² + 126*f^2*bη*bη²))/3675. +
            (4*f*(45*b1 - 17*f - 17*f*bη)*bΠ2para)/1575.)*xi0m1 - (64*b1*f^2*xi2m1)/1715. - (16*b1*f^3*xi2m1)/525. +
        (64*f^2*bK²*xi2m1)/1225. - (16*b1*f*bKKpara*xi2m1)/175. + (496*f^2*bKKpara*xi2m1)/11025. -
        (16*f^3*bδη*xi2m1)/525. - (64*b1*f^2*bη*xi2m1)/1715. + (64*f^3*bη*xi2m1)/3675. - (16*b1*f^3*bη*xi2m1)/175. +
        (64*f^4*bη*xi2m1)/1225. + (64*f^2*bK²*bη*xi2m1)/1225. + (496*f^2*bKKpara*bη*xi2m1)/11025. -
        (16*f^3*bδη*bη*xi2m1)/525. + (64*f^3*bη^2*xi2m1)/3675. + (64*f^4*bη^2*xi2m1)/1225. - (32*b1*f^3*bη²*xi2m1)/525. +
        (64*f^4*bη²*xi2m1)/1225. + (64*f^4*bη*bη²*xi2m1)/1225. - (16*b1*f*bΠ2para*xi2m1)/175. +
        (16*f^2*bΠ2para*xi2m1)/441. + (16*f^2*bη*bΠ2para*xi2m1)/441. + (16*f^2*(1 + bη)*(495*b1 - 693*bK² - 462*bKKpara -
                231*f*bη - 252*f^2*bη - 252*f^2*bη² - 77*bΠ2para)*xi4m1)/282975.)

    M[1,4,1] = @. xi10*(((-32*f^2*bη*(25*b1 - 7*f - 7*f*bη))/6125. + (32*f*(25*b1 - 7*f - 7*f*bη)*bΠ2para)/3675.)*xi0m1 -
        (64*f*(5*b1 + f + f*bη)*(3*f*bη - 5*bΠ2para)*xi2m1)/18375.) + xi30*((32*f*(35*b1 - 13*f - 13*f*bη)*(3*f*bη -
                5*bΠ2para)*xi2m1)/42875. + (256*f^2*(1 + bη)*(3*f*bη - 5*bΠ2para)*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η)*(15*b1*f - 2*bKKpara + 15*f*bδη - 18*f^2*bη - 18*f^2*bη² - 12*bΠ2para)*xi10*xi4m1)/2079. +
    xi30*((40*f^2*(21*b1*f - 4*bKKpara + 21*f*bδη + 63*b1*f*bη - 36*f^2*bη - 4*bKKpara*bη + 21*f*bδη*bη - 36*f^2*bη^2 +
                42*b1*f*bη² - 36*f^2*bη² - 36*f^2*bη*bη² - 18*bΠ2para - 18*bη*bΠ2para)*xi2m1)/14553. +
        (32*f^2*(1 + bη)*(13*bKKpara + 12*f^2*bη + 12*f^2*bη² + 13*bΠ2para)*xi4m1)/21021.)

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

    M[1,6,5] = @. (16*f^2*(1 + bη)*(15*b1*f - 2*bKKpara + 15*f*bδη - 18*f^2*bη - 18*f^2*bη² -
            12*bΠ2para)*xi10*xi4m1)/3465. + xi30*((16*f^2*(21*b1*f - 4*bKKpara + 21*f*bδη + 63*b1*f*bη - 36*f^2*bη -
                4*bKKpara*bη + 21*f*bδη*bη - 36*f^2*bη^2 + 42*b1*f*bη² - 36*f^2*bη² - 36*f^2*bη*bη² - 18*bΠ2para -
                18*bη*bΠ2para)*xi2m1)/4851. + (64*f^2*(1 + bη)*(13*bKKpara + 12*f^2*bη + 12*f^2*bη² +
                13*bΠ2para)*xi4m1)/35035.)

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

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

    M[2,0,0] = @. ((f^2*(15*b1^2 - 10*b1*f + 3*f^2))/270. - ((5*b1 - 3*f)*f^3*bη)/135. + (f^4*bη^2)/90.)*xi0m2*xi00 +
    ((-2*b1^2)/7. - (b1*b2)/2. + (b1^2*f)/3. - (2*b1*f^2)/63. - (b1^2*f^2)/18. - (b2*f^2)/18. + (2*b1*f^3)/15. -
        f^4/50. - (4*b1*bK²)/15. - (4*f^2*bK²)/135. - (4*b1*bKKpara)/45. - (34*f^2*bKKpara)/2025. + (b1*f*bδη)/3. +
        (f^3*bδη)/15. + (16*b1*f*bη)/105. + (b2*f*bη)/6. - (58*b1*f^2*bη)/225. + (2*f^3*bη)/175. + (b1*f^3*bη)/15. -
        (2*f^4*bη)/25. + (4*f*bK²*bη)/45. + (4*f*bKKpara*bη)/135. - (f^2*bδη*bη)/9. - (6*f^2*bη^2)/175. +
        (11*f^3*bη^2)/225. - (f^4*bη^2)/50. - (11*b1*f^2*bη²)/75. - (f^4*bη²)/25. + (11*f^3*bη*bη²)/225. -
        ((465*b1 + 107*f^2 - 195*f*bη)*bΠ2para)/1575.)*xi1m1^2 + ((-4*(7*b1 - 3*f)*f^3)/945. -
        (4*(7*b1 - 6*f)*f^3*bη)/945. + (4*f^4*bη^2)/315.)*xi2m2*xi20 + ((2*(225*b1^2 + 25*b1*f^2 - 3*f^4 - 315*b1*bK² -
                35*f^2*bK² - 105*b1*bKKpara - 14*f^2*bKKpara - 120*b1*f*bη - 42*b1*f^2*bη - 9*f^3*bη - 12*f^4*bη +
                105*f*bK²*bη + 35*f*bKKpara*bη + 27*f^2*bη^2 + 14*f^3*bη^2 - 3*f^4*bη^2 - 42*b1*f^2*bη² - 6*f^4*bη² +
                14*f^3*bη*bη²))/1575. - (2*(30*b1 - f^2 + 10*f*bη)*bΠ2para)/1575.)*xi3m1^2

    M[2,0,2] = @. ((2*f*(10080*b1*f + 11025*b1^2*f + 11025*b2*f - 26460*b1*f^2 + 3969*f^3 + 588*f*bK² + 686*f*bKKpara -
                13230*f^2*bδη + 8820*b1*f*bη - 5184*f^2*bη - 13230*b1*f^2*bη + 15876*f^3*bη - 2058*bKKpara*bη +
                8820*f*bδη*bη - 7056*f^2*bη^2 + 3969*f^3*bη^2 + 7938*f^3*bη² - 7056*f^2*bη*bη²))/496125. +
        (4*f*(872*f - 441*bη)*bΠ2para)/55125.)*xi1m1^2 + ((-4*(7*b1 - 3*f)*f^3)/945. - (4*(7*b1 - 6*f)*f^3*bη)/945. +
        (4*f^4*bη^2)/315.)*xi00*xi2m2 + ((4*f*(420*b1*f + 630*b1*f^2 - 189*f^3 - 588*f*bK² - 336*f*bKKpara + 315*f^2*bδη -
                735*b1*f*bη - 216*f^2*bη + 315*b1*f^2*bη - 756*f^3*bη + 98*bKKpara*bη - 735*f*bδη*bη + 546*f^2*bη^2 -
                189*f^3*bη^2 - 378*f^3*bη² + 546*f^2*bη*bη²))/55125. - (8*f*(31*f - 98*bη)*bΠ2para)/18375.)*xi1m1*xi3m1 +
    ((-16*f*(35*b1*f - 7*f^3 - 49*f*bK² - 28*f*bKKpara - 18*f^2*bη - 28*f^3*bη + 49*bKKpara*bη + 28*f^2*bη^2 -
                7*f^3*bη^2 - 14*f^3*bη² + 28*f^2*bη*bη²))/55125. - (16*f*(2*f + 49*bη)*bΠ2para)/55125.)*xi3m1^2 +
    xi20*(((2*f^2*(15*b1^2 - 10*b1*f + 3*f^2))/675. - (4*(5*b1 - 3*f)*f^3*bη)/675. + (2*f^4*bη^2)/225.)*xi0m2 +
        ((16*(7*b1 - 3*f)*f^3)/6615. + (16*(7*b1 - 6*f)*f^3*bη)/6615. - (16*f^4*bη^2)/2205.)*xi2m2 +
        ((32*f^4)/11025. + (64*f^4*bη)/11025. + (32*f^4*bη^2)/11025.)*xi4m2)

    M[2,2,0] = @. ((f^2*(15*b1^2 - 10*b1*f + 3*f^2))/135. - (2*(5*b1 - 3*f)*f^3*bη)/135. + (f^4*bη^2)/45.)*xi0m2*xi00 +
    (-(f*(6300*b1*f + 11025*b1^2*f + 11025*b2*f - 26460*b1*f^2 + 3969*f^3 + 5880*f*bK² + 3332*f*bKKpara - 13230*f^2*bδη -
                30240*b1*bη - 33075*b2*bη + 22050*b1*f*bη - 3240*f^2*bη - 13230*b1*f^2*bη + 15876*f^3*bη - 17640*bK²*bη -
                5880*bKKpara*bη + 22050*f*bδη*bη + 9720*f*bη^2 - 9702*f^2*bη^2 + 3969*f^3*bη^2 + 7938*f^3*bη² -
                9702*f^2*bη*bη²))/99225. - (2*(1050*b1 + 839*f^2 - 1635*f*bη)*bΠ2para)/11025.)*xi1m1^2 +
    ((-8*(7*b1 - 3*f)*f^3)/945. - (8*(7*b1 - 6*f)*f^3*bη)/945. + (8*f^4*bη^2)/315.)*xi2m2*xi20 +
    ((4*f*(175*b1*f - 21*f^3 - 245*f*bK² - 98*f*bKKpara - 840*b1*bη - 90*f^2*bη - 84*f^3*bη + 735*bK²*bη + 245*bKKpara*bη +
                270*f*bη^2 + 98*f^2*bη^2 - 21*f^3*bη^2 - 42*f^3*bη² + 98*f^2*bη*bη²))/11025. +
        (4*(525*b1 + 52*f^2 - 205*f*bη)*bΠ2para)/11025.)*xi3m1^2

    M[2,2,2] = @. ((-26460*b1^2*f + 15840*b1*f^2 + 17325*b1^2*f^2 + 17325*b2*f^2 - 41580*b1*f^3 + 6237*f^4 + 924*f^2*bK² +
            6174*b1*bKKpara + 1078*f^2*bKKpara - 26460*b1*f*bδη - 20790*f^3*bδη + 35028*b1*f^2*bη - 7776*f^3*bη -
            20790*b1*f^3*bη + 24948*f^4*bη - 3234*f*bKKpara*bη + 13860*f^2*bδη*bη - 11088*f^3*bη^2 + 6237*f^4*bη^2 +
            21168*b1*f^2*bη² + 12474*f^4*bη² - 11088*f^3*bη*bη²)/99225. + (2*(1323*b1 + 1336*f^2 -
                693*f*bη)*bΠ2para)/11025.)*xi1m1^2 + ((-22*(7*b1 - 3*f)*f^3)/1323. - (22*(7*b1 - 6*f)*f^3*bη)/1323. +
        (22*f^4*bη^2)/441.)*xi00*xi2m2 + ((2*(2205*b1^2*f + 660*b1*f^2 + 990*b1*f^3 - 297*f^4 - 924*f^2*bK² -
                294*b1*bKKpara - 528*f^2*bKKpara + 2205*b1*f*bδη + 495*f^3*bδη - 2793*b1*f^2*bη - 324*f^3*bη +
                495*b1*f^3*bη - 1188*f^4*bη + 154*f*bKKpara*bη - 1155*f^2*bδη*bη + 858*f^3*bη^2 - 297*f^4*bη^2 -
                1638*b1*f^2*bη² - 594*f^4*bη² + 858*f^3*bη*bη²))/11025. - (4*(294*b1 + 53*f^2 -
                154*f*bη)*bΠ2para)/3675.)*xi1m1*xi3m1 + ((8*(-55*b1*f^2 + 11*f^4 + 77*f^2*bK² + 147*b1*bKKpara +
                44*f^2*bKKpara + 84*b1*f^2*bη + 27*f^3*bη + 44*f^4*bη - 77*f*bKKpara*bη - 44*f^3*bη^2 + 11*f^4*bη^2 +
                84*b1*f^2*bη² + 22*f^4*bη² - 44*f^3*bη*bη²))/11025. + (8*(147*b1 - f^2 - 77*f*bη)*bΠ2para)/11025.)*xi3m1^2 +
    xi20*(((11*f^2*(15*b1^2 - 10*b1*f + 3*f^2))/945. - (22*(5*b1 - 3*f)*f^3*bη)/945. + (11*f^4*bη^2)/315.)*xi0m2 +
        ((88*(7*b1 - 3*f)*f^3)/9261. + (88*(7*b1 - 6*f)*f^3*bη)/9261. - (88*f^4*bη^2)/3087.)*xi2m2 + ((176*f^4)/15435. +
            (352*f^4*bη)/15435. + (176*f^4*bη^2)/15435.)*xi4m2)

    M[2,2,4] = @. ((-32*f^2*(30*b1*f - 9*f^2 - 2*bKKpara + 15*f*bδη + 15*b1*f*bη - 36*f^2*bη + 12*f*bη^2 - 9*f^2*bη^2 -
                18*f^2*bη² + 12*f*bη*bη²))/11025. + (128*f^2*bΠ2para)/3675.)*xi1m1*xi3m1 + ((-16*f^2*(9*f^2 + 22*bKKpara +
                36*f^2*bη - 22*f*bη^2 + 9*f^2*bη^2 + 18*f^2*bη² - 22*f*bη*bη²))/40425. - (32*f^2*bΠ2para)/3675.)*xi3m1^2 +
    ((16*f^4)/2205. + (32*f^4*bη)/2205. + (16*f^4*bη^2)/2205.)*xi00*xi4m2 + xi20*(((-32*(7*b1 - 3*f)*f^3)/5145. -
            (32*(7*b1 - 6*f)*f^3*bη)/5145. + (32*f^4*bη^2)/1715.)*xi2m2 + ((-128*f^4)/33957. - (256*f^4*bη)/33957. -
            (128*f^4*bη^2)/33957.)*xi4m2)

    M[2,4,0] = @. ((16*f^2*(f - 3*bη)*bη)/1225. - (16*f*(f - 3*bη)*bΠ2para)/735.)*xi1m1^2 +
    ((-16*f^2*(f - 3*bη)*bη)/1225. + (16*f*(f - 3*bη)*bΠ2para)/735.)*xi3m1^2

    M[2,4,2] = @. ((4*f*(15840*b1*f + 17325*b1^2*f + 17325*b2*f - 41580*b1*f^2 + 6237*f^3 + 924*f*bK² + 1078*f*bKKpara -
                20790*f^2*bδη + 13860*b1*f*bη - 9792*f^2*bη - 20790*b1*f^2*bη + 24948*f^3*bη - 3234*bKKpara*bη +
                13860*f*bδη*bη - 11088*f^2*bη^2 + 6237*f^3*bη^2 + 12474*f^3*bη² - 11088*f^2*bη*bη²))/606375. +
        (8*f*(4568*f - 2079*bη)*bΠ2para)/202125.)*xi1m1^2 + ((-8*(7*b1 - 3*f)*f^3)/735. - (8*(7*b1 - 6*f)*f^3*bη)/735. +
        (8*f^4*bη^2)/245.)*xi00*xi2m2 + ((8*f*(660*b1*f + 990*b1*f^2 - 297*f^3 - 924*f*bK² - 528*f*bKKpara + 495*f^2*bδη -
                1155*b1*f*bη - 408*f^2*bη + 495*b1*f^2*bη - 1188*f^3*bη + 154*bKKpara*bη - 1155*f*bδη*bη + 858*f^2*bη^2 -
                297*f^3*bη^2 - 594*f^3*bη² + 858*f^2*bη*bη²))/67375. - (16*f*(89*f - 462*bη)*bΠ2para)/67375.)*xi1m1*xi3m1 +
    ((-32*f*(55*b1*f - 11*f^3 - 77*f*bK² - 44*f*bKKpara - 34*f^2*bη - 44*f^3*bη + 77*bKKpara*bη + 44*f^2*bη^2 -
                11*f^3*bη^2 - 22*f^3*bη² + 44*f^2*bη*bη²))/67375. - (32*f*(38*f + 231*bη)*bΠ2para)/202125.)*xi3m1^2 +
    xi20*(((4*f^2*(15*b1^2 - 10*b1*f + 3*f^2))/525. - (8*(5*b1 - 3*f)*f^3*bη)/525. + (4*f^4*bη^2)/175.)*xi0m2 +
        ((32*(7*b1 - 3*f)*f^3)/5145. + (32*(7*b1 - 6*f)*f^3*bη)/5145. - (32*f^4*bη^2)/1715.)*xi2m2 + ((64*f^4)/8575. +
            (128*f^4*bη)/8575. + (64*f^4*bη^2)/8575.)*xi4m2)

    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η²))/13475. +
        (1248*f^2*bΠ2para)/13475.)*xi1m1*xi3m1 + ((-4*f^2*(351*f^2 + 858*bKKpara + 1694*b1*bη + 1404*f^2*bη -
                858*f*bη^2 + 351*f^2*bη^2 + 1694*b1*bη² + 702*f^2*bη² - 858*f*bη*bη²))/148225. -
        (312*f^2*bΠ2para)/13475.)*xi3m1^2 + ((52*f^4)/2695. + (104*f^4*bη)/2695. + (52*f^4*bη^2)/2695.)*xi00*xi4m2 +
    xi20*(((-312*(7*b1 - 3*f)*f^3)/18865. - (312*(7*b1 - 6*f)*f^3*bη)/18865. + (936*f^4*bη^2)/18865.)*xi2m2 +
        ((-416*f^4)/41503. - (832*f^4*bη)/41503. - (416*f^4*bη^2)/41503.)*xi4m2)

    M[2,4,6] = @. ((80*f^4*(1 + 4*bη + bη^2))/11011. + (160*f^4*bη²)/11011.)*xi3m1^2 +
    ((80*f^4)/11011. + (160*f^4*bη)/11011. + (80*f^4*bη^2)/11011.)*xi20*xi4m2

    M[2,6,2] = @. ((-256*f^3*bη)/13475. + (256*f^2*bΠ2para)/8085.)*xi1m1^2 + ((-192*f^3*bη)/13475. +
        (64*f^2*bΠ2para)/2695.)*xi1m1*xi3m1 + ((64*f^3*bη)/13475. - (64*f^2*bΠ2para)/8085.)*xi3m1^2

    M[2,6,4] = @. ((-16*f^2*(30*b1*f - 9*f^2 - 2*bKKpara + 15*f*bδη + 15*b1*f*bη - 36*f^2*bη + 12*f*bη^2 - 9*f^2*bη^2 -
                18*f^2*bη² + 12*f*bη*bη²))/3465. + (64*f^2*bΠ2para)/1155.)*xi1m1*xi3m1 + ((-8*f^2*(9*f^2 + 22*bKKpara +
                36*f^2*bη - 22*f*bη^2 + 9*f^2*bη^2 + 18*f^2*bη² - 22*f*bη*bη²))/12705. - (16*f^2*bΠ2para)/1155.)*xi3m1^2 +
    ((8*f^4)/693. + (16*f^4*bη)/693. + (8*f^4*bη^2)/693.)*xi00*xi4m2 + xi20*(((-16*(7*b1 - 3*f)*f^3)/1617. -
            (16*(7*b1 - 6*f)*f^3*bη)/1617. + (16*f^4*bη^2)/539.)*xi2m2 +
        ((-320*f^4)/53361. - (640*f^4*bη)/53361. - (320*f^4*bη^2)/53361.)*xi4m2)

    M[2,6,6] = @. ((664*f^4*(1 + 4*bη + bη^2))/38115. + (1328*f^4*bη²)/38115.)*xi3m1^2 + ((664*f^4)/38115. +
        (1328*f^4*bη)/38115. + (664*f^4*bη^2)/38115.)*xi20*xi4m2

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

    M[3,0,1] = @. ((f*(25*b1^2 - 15*b1*f + 15*b1*f^2 - 9*f^3))/225. - (f^2*(10*b1 - 3*f + 3*f^2)*bη)/75. +
        (f^3*bη^2)/25.)*xi0m1*xi1m2 - (2*f*(125*b1^2 - 45*b1*f + 30*b1*f^2 - 18*f^3 - 120*b1*f*bη + 27*f^2*bη -
            18*f^3*bη + 27*f^2*bη^2)*xi1m2*xi2m1)/1125. - (4*f^2*(1 + bη)*(35*b1 + 9*f^2 - 21*f*bη)*xi2m1*xi3m2)/2625.

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

    M[3,2,1] = @. ((2*f*(35*b1^2 - 21*b1*f + 30*b1*f^2 - 18*f^3))/315. - (2*f^2*(17*b1 - 6*f + 6*f^2)*bη)/105. +
        (4*f^3*bη^2)/35.)*xi0m1*xi1m2 - (4*f*(175*b1^2 - 63*b1*f + 60*b1*f^2 - 36*f^3 - 213*b1*f*bη + 54*f^2*bη -
            36*f^3*bη + 54*f^2*bη^2)*xi1m2*xi2m1)/1575. - (8*f^2*(1 + bη)*(49*b1 + 18*f^2 - 42*f*bη)*xi2m1*xi3m2)/3675.

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

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

    M[3,4,1] = @. ((8*(5*b1 - 3*f)*f^3)/525. - (8*f^2*(5*b1 - 3*f + 3*f^2)*bη)/525. + (8*f^3*bη^2)/175.)*xi0m1*xi1m2 -
    (16*f^2*(10*b1*f - 6*f^2 - 25*b1*bη + 9*f*bη - 6*f^2*bη + 9*f*bη^2)*xi1m2*xi2m1)/2625. -
    (32*f^3*(3*f - 7*bη + 3*f*bη - 7*bη^2)*xi2m1*xi3m2)/6125.

    M[3,4,3] = @. (16*f^2*(-55*b1 + 105*b1*f - 63*f^2 - 55*b1*bη + 42*f*bη - 63*f^2*bη + 42*f*bη^2)*xi1m2*xi2m1)/9625. +
    ((-8*f^2*(55*b1 + 42*f^2))/5775. - (8*f^2*(55*b1 - 42*f + 42*f^2)*bη)/5775. + (16*f^3*bη^2)/275.)*xi0m1*xi3m2 +
    (32*f^2*(1 + bη)*(55*b1 + 28*f^2 - 42*f*bη)*xi2m1*xi3m2)/28875.

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

    M[3,6,3] = @. (16*f^3*(5*b1 - 3*f + 2*bη - 3*f*bη + 2*bη^2)*xi1m2*xi2m1)/1155. +
    ((-16*f^4)/693. - (16*(-1 + f)*f^3*bη)/693. + (16*f^3*bη^2)/693.)*xi0m1*xi3m2 +
    (64*f^3*(2*f - 3*bη + 2*f*bη - 3*bη^2)*xi2m1*xi3m2)/10395.

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

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

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

    M[4,0,2] = @. ((4*f^2*(7*b1 + 3*f^2))/945. - (8*f^3*bη)/315.)*xi0m2*xi2m2 +
    ((-2*f^2*(49*b1 + 12*f^2))/6615. + (4*f^3*bη)/315.)*xi2m2^2

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

    M[4,2,0] = @. ((f^2*(7*b1 + 3*f^2))/189. - (f*(7*b1 + 2*f^2)*bη)/63. + (f^2*bη^2)/21.)*xi0m2^2 +
    ((2*f^2*(35*b1 + 6*f^2))/945. - (2*f*(7*b1 + 2*f^2)*bη)/63. + (2*f^2*bη^2)/21.)*xi2m2^2

    M[4,2,2] = @. ((f^2*(22*b1 + 9*f^2))/189. - (2*f^3*bη)/21.)*xi0m2*xi2m2 +
    (-(f^2*(77*b1 + 18*f^2))/1323. + (f^3*bη)/21.)*xi2m2^2

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

    M[4,4,0] = @. ((2*f^4)/315. - (4*f^3*bη)/315. + (2*f^2*bη^2)/105.)*xi0m2^2 +
    ((8*f^4)/1575. - (8*f^3*bη)/315. + (4*f^2*bη^2)/105.)*xi2m2^2

    M[4,4,2] = @. ((8*f^2*(33*b1 + 17*f^2))/3465. - (272*f^3*bη)/3465.)*xi0m2*xi2m2 +
    ((-4*f^2*(231*b1 + 68*f^2))/24255. + (136*f^3*bη)/3465.)*xi2m2^2

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

    M[4,6,2] = @. ((8*f^4)/693. - (16*f^3*bη)/693.)*xi0m2*xi2m2 + ((-16*f^4)/4851. + (8*f^3*bη)/693.)*xi2m2^2

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

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

    return r, M
end
#generateMnli(x->x, f=1.0, BiasDict22=Dict([("b1",0), ("bη", 0), ("b2", 0), ("bK2", 0), ("bδη", 0),
#                           ("bη2", 0), ("bKK∥", 0), ("bΠ2∥", 0)]));
#precompile(ξ, (Function, Int64, Int64, Int64))
#precompile(ξ, (Function, Int64, Int64))
#precompile(generateMnli, (Function, Dict{String,Int64}, Float64))
#end

max_stars_repo_name: JosephTomlinson/GeneralBiasPk

max_issues_repo_name: JosephTomlinson/GeneralBiasPk

max_forks_repo_name: JosephTomlinson/GeneralBiasPk


text:
abstract type AbstractEA <: AbstractAlgorithm end
using LaTeXStrings

using DataFrames

struct ea1pλwith2rates <: AbstractEA
    λ::Integer
    stop::AbstractStop
    name::LaTeXString
    function ea1pλwith2rates(;
        λ::Integer = 10,
        stop::AbstractStop = FixedBudget(1000),
        name::LaTeXString = L"(1+λ)EA with 2 rates",
    )
        new(λ, stop, name)
    end
end

function optimize(x, setting::ea1pλwith2rates)
    λ = setting.λ
    trace = Trace(setting, x)
    x = initial(x)
    n = length(x)
    # bitrand returns a random bit string.
    r = 2

    for iter ∈ 1:niterations(setting.stop)
        α = copy(x)
        rα = r
        for i = 1:λ
            if i ≤ λ / 2
                ry = r / (2n)
                y = mutation(x, UniformlyIndependentMutation(ry))
            else
                ry = 2r / n
                y = mutation(x, UniformlyIndependentMutation(ry))
            end

            if fitness(y) ≥ fitness(α)
                α = y
                rα = ry
            end
        end

        # The second condition is for implementing "breaking ties randomly".
        if fitness(α) ≥ fitness(x)
            if fitness(α) > fitness(x)
                record(trace, α, iter, isoptimum(α))
            end
            x = α

            if isoptimum(x)
                return trace
            end
        end

        if rand() < 0.5
            r = rα
        else
            r = if (rand() < 0.5)
                r = r / 2
            else
                r = 2r
            end
        end
        r = min(max(r, 2), n / 4)
    end
    trace
end

struct ea1p1 <: AbstractEA
    mutation::Mutation
    stop::AbstractStop
    name::LaTeXString
    function ea1p1(;
        stop::AbstractStop = FixedBudget(1000),
        mutation::Mutation = UniformlyIndependentMutation(0.5),
        name::LaTeXString = L"(1+1)EA",
    )
        new(mutation, stop, name)
    end
end

function optimize(x, setting::ea1p1)
    trace = Trace(setting, x)
    x = initial(x)
    n = length(x)
    # bitrand returns a random bit string.

    for iter ∈ 1:niterations(setting.stop)
        # println("Enters")

        y = mutation(x, setting.mutation)

        # The second condition is for implementing "breaking ties randomly".
        if fitness(y) ≥ fitness(x)
            if fitness(y) > fitness(x)
                # println(fitness(y))
                record(trace, y, iter, isoptimum(y))
            end
            x = y

            if isoptimum(x)
                # println("Optimum found")
                return trace
            end
        end
    end
    # println("Finished")
    trace
end

struct ea1p1SD <: AbstractEA
    R::Real
    stop::AbstractStop
    thresholds::Array
    name::LaTeXString
    function ea1p1SD(;
        R::Real = 1,
        stop::AbstractStop = FixedBudget(1000),
        thresholds = [typemax(Int) for _ = 1:10],
        name::LaTeXString = name=L"SD-(1+1)EA"
    )
        new(R, stop, thresholds, name)
    end
end

SDCounter(n::Integer, r::Real, R::Real) =
    (n / r)^r * (n / (n - r))^(n - r) * log(Base.MathConstants.e * n * R)
SDCounterEstimation(n::Integer, r::Real, R::Real) =
    (Base.MathConstants.e * n / r)^r * log(Base.MathConstants.e * n * R)

function threshold_gen(generator::Function, n::Integer, R::Real)
    thresh = []
    for r = 1:ceil(Integer, n)
        val = generator(n, r, R)
        push!(thresh, val)
        if val >= typemax(Int) / n
            break
        end
    end
    thresh
end

function optimize(x, setting::ea1p1SD)
    trace = Trace(setting, x)
    x = initial(x)
    n = length(x)
    thresholds = setting.thresholds
    r = 1
    u = 0
    for iter ∈ 1:niterations(setting.stop)
        y = mutation(x, UniformlyIndependentMutation(r // n))
        u += 1
        # The second condition is for implementing "breaking ties randomly".
        if fitness(y) > fitness(x)
            x = y
            u = 0
            r = 1
            # println("New rate", r)
            # println(fitness(x), " in iteration= ", iter)
            record(trace, y, iter, isoptimum(y))
            if isoptimum(x)
                # println("The Optimum is found.")
                return trace
            end
        elseif fitness(y) == fitness(x) && r == 1
            x = y
        end

        if u > thresholds[r]
            r = min(r + 1, ceil(Integer, n / 2))
            # println("New rate", r)
            u = 0
        end
    end
    trace
end

struct ea1pλSASD <: AbstractEA
    R::Real
    λ::Integer
    stop::AbstractStop
    thresholds::Array
    function ea1pλSASD(;
        R::Real = 1,
        λ::Integer = 10,
        stop::AbstractStop = FixedBudget(1000),
        thresholds = [typemax(Int) for _ = 1:10],
    )
        new(R, λ, stop, thresholds)
    end
end

function optimize(x, setting::ea1pλSASD)
    trace = Trace(setting, x)
    x = initial(x)
    n = length(x)
    thresholds = setting.thresholds
    λ = setting.λ
    r_init = 2
    r = r_init
    u = 0
    g = false

    for iter ∈ 1:niterations(setting.stop)
        u += 1
        if g == false
            y = copy(x)
            ry = r
            for i = 1:λ
                if i ≤ λ / 2
                    rα = r / (2n)
                    α = mutation(x, UniformlyIndependentMutation(rα))
                else
                    rα = 2r / n
                    α = mutation(x, UniformlyIndependentMutation(rα))
                end

                if fitness(α) ≥ fitness(y)
                    y = α
                    ry = rα
                end
            end
            if fitness(y) ≥ fitness(x)

                if fitness(y) > fitness(x)
                    # println("New rate $r in $g")
                    record(trace, y, iter, isoptimum(y))
                end
                x = y

                if isoptimum(x)
                    return trace
                end
            end

            if rand() < 0.5
                r = ry
            else
                r = if (rand() < 0.5)
                    r = r / 2
                else
                    r = 2r
                end
            end
            r = floor(Integer, min(max(r, 2), n / 4))

            if u > thresholds[1] / λ
                r = 2
                # println("New rate $r in $g")
                # println("STAG detection.")
                g = true
                u = 0
            end

        else
            y = copy(x)
            for i = 1:λ
                α = mutation(x, UniformlyIndependentMutation(r / n))
                if fitness(α) ≥ fitness(y)
                    y = α
                end
            end
            # The second condition is for implementing "breaking ties randomly".
            if fitness(y) > fitness(x)
                x = y
                r = r_init
                # println("New rate $r in $g")
                g = false
                u = 0
                record(trace, y, iter, isoptimum(y))
                if isoptimum(x)
                    # println("The Optimum is found.")
                    return trace
                end
            end

            if u > thresholds[r] / λ
                r = min(r + 1, ceil(Integer, n / 2))
                # println("New rate $r in $g")
                u = 0
            end
        end
    end
    trace
end

struct RLSSDstar <: AbstractEA
    R::Real
    stop::AbstractStop
    name::LaTeXString
    function RLSSDstar(;
        R::Real = 1,
        stop::AbstractStop = FixedBudget(1000),
        name::LaTeXString = L"SD-RLS^r",
    )
        new(R, stop, name)
    end
end

RLSSDCounter(n::Integer, r::Real, R::Real) = binomial(n, r) * log(R)

function optimize(x, setting::RLSSDstar)
    trace = Trace(setting, x)
    x = initial(x)
    n = length(x)
    thresholds = threshold_gen(RLSSDCounter, n, setting.R)
    r = 1
    s = 1
    u = 0
    for iter ∈ 1:niterations(setting.stop)
        y = mutation(x, KBitFlip(s))
        u += 1
        # The second condition is for implementing "breaking ties randomly".
        if fitness(y) > fitness(x)
            x = y
            u = 0
            r = 1
            s = 1
            record(trace, y, iter, isoptimum(y))
            if isoptimum(x)
                return trace
            end
        elseif fitness(y) == fitness(x) && s == 1
            x = y
        end

        if u > thresholds[s]
            if s == 1
                if r < n / 2
                    r = r + 1
                else
                    r = n
                end
                s = r
            else
                s = s - 1
            end
            u = 0
            # println("RLSSDstar ", r, " ", s, " ", fitness(x), isoptimum(x))
        end
    end
    trace
end

struct RLS12 <: AbstractEA
    stop::AbstractStop
    name::LaTeXString
    function RLS12(;
        stop::AbstractStop = FixedBudget(1000),
        name::LaTeXString = L"RLS^{1,2}",
    )
        new(stop, name)
    end
end

function optimize(x, setting::RLS12)
    trace = Trace(setting, x)
    x = initial(x)
    n = length(x)
    # bitrand returns a random bit string.

    for iter ∈ 1:niterations(setting.stop)
        if rand() < 0.5
            y = mutation(x, KBitFlip(1))
        else
            y = mutation(x, KBitFlip(2))
        end

        # The second condition is for implementing "breaking ties randomly".
        if fitness(y) ≥ fitness(x)
            if fitness(y) > fitness(x)
                record(trace, y, iter, isoptimum(y))
            end
            x = y

            if isoptimum(x)
                return trace
            end
        end
    end
    trace
end

struct RLSSDstarstar <: AbstractEA
    R::Real
    stop::AbstractStop
    name::LaTeXString
    function RLSSDstarstar(;
        R::Real = 1,
        stop::AbstractStop = FixedBudget(1000),
        name::LaTeXString = L"SD-RLS^m",
    )
        new(R, stop, name)
    end
end

function optimize(x, setting::RLSSDstarstar)
    trace = Trace(setting, x)
    x = initial(x)
    n = length(x)
    thresholds = threshold_gen(RLSSDCounter, n, setting.R)
    r = 1
    s = 1
    u = 0
    B = typemax(Int)
    for iter ∈ 1:niterations(setting.stop)
        y = mutation(x, KBitFlip(s))
        u += 1
        # The second condition is for implementing "breaking ties randomly".
        if fitness(y) > fitness(x)
            x = y
            r = s
            s = 1
            if r > 1
                B = u / (r - 1) * 1 / log(n)
            else
                B = typemax(Int)
            end
            u = 0
            record(trace, y, iter, isoptimum(y))
            if isoptimum(x)
                return trace
            end
        elseif fitness(y) == fitness(x) && s == 1
            x = y
        end

        if u > min(B, thresholds[s])
            if s == r
                if r < n / 2 - 1
                    r = r + 1
                else
                    r = n
                end
                s = 1
            else
                s = s + 1
                if s == r
                    B = typemax(Int)
                end
            end
            u = 0
            # println("RLSSDstarstar ", r, " ", s, " ", fitness(x))
        end
    end
    trace
end

struct fastSD <: AbstractEA
    β::Real
    γ::Real
    R::Real
    stop::AbstractStop
    name::LaTeXString
    function fastSD(;
        β::Real = 1.5,
        γ::Real = 2/3,
        R::Real = 1,
        stop::AbstractStop = FixedBudget(1000),
        name::LaTeXString = L"fast-SD",
    )
        new(β, γ, R, stop, name)
    end
end

function threshold_gen(generator::Function, n::Integer, γ::Real, R::Real)
    thresh = []
    for r = 1:ceil(Integer, n)
        val = generator(n, r, γ, R)
        push!(thresh, val)
        if val >= typemax(Int) / n
            break
        end
    end
    thresh
end

fastSDCounter(n::Integer, s::Real, γ::Real, R::Real) = binomial(n, s) * log(R) / (1-γ)

function optimize(x, setting::fastSD)
    trace = Trace(setting, x)
    x = initial(x)
    n = length(x)
    γ = setting.γ
    β = setting.β
    thresholds = threshold_gen(fastSDCounter, n, γ, setting.R)
    r = 1
    u = 0
    for iter ∈ 1:niterations(setting.stop)
        s = r
        if rand() < γ
            if rand() ≤ 0.5
                α = discretepowerlaw(β, n-r)
                s = r + rand(α)
            else
                α = discretepowerlaw(β, max(1,r-1))
                s = r - rand(α)
            end
        end
        y = mutation(x, KBitFlip(s))
        u += 1
        if fitness(y) > fitness(x)
            x = y
            u = 0
            r = 1
            record(trace, y, iter, isoptimum(y))
            if isoptimum(x)
                return trace
            end
        elseif fitness(y) == fitness(x) && r == 1
            x = y
        end

        if u > thresholds[r]
            r = min(r+1, floor(Integer, n/2))
            u = 0
        end
    end
    trace
end

ext: jl

max_stars_repo_name: DTUComputeTONIA/SDRLS

max_issues_repo_name: DTUComputeTONIA/SDRLS

max_forks_repo_name: DTUComputeTONIA/SDRLS


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((-1+2x)^-7//2*(-5+3x)^2, x) == :(-49//20*(-1+2x)^-5//2+-9//4*(-1+2x)^-1//2+7//2*(-1+2x)^-3//2)
@test integrate((-5*exp(-1*m*x)+2*exp(m*x))^-1, x) == :(-1//10*10^(1/2)*m^-1*arctanh(1//5*10^(1/2)*exp(m*x)))
@test integrate((a+b*cos(x))^-1, x) == :(2*(a+b)^-1//2*(a+-1b)^-1//2*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tan((1/2)*x)))
@test integrate((3+3*cos(x)+4*sin(x))^-1, x) == :(1//4*log(3+4*tan((1/2)*x)))
@test integrate((4+3*cos(x)+4*sin(x))^-1, x) == :(-1//3*log(4+3*cot((1/2)*x+1//4*pi)))
@test integrate((6+3*cos(x)+4*sin(x))^-1, x) == :(1//11*x*11^(1/2)+2//11*11^(1/2)*arctan((6+11^(1/2)+3*cos(x)+4*sin(x))^-1*(-3*sin(x)+4*cos(x))))
@test integrate((1/2)*log((x^2+-1*a^2)^2), x) == :(-2x+(1/2)*x*log((x^2+-1*a^2)^2)+2*a*arctanh(x*a^-1))

max_stars_repo_path: test/0 Independent test suites/Wester Problems.jl


text:
using BinaryProvider # requires BinaryProvider 0.4.0 or later

const forcecompile = get(ENV, "JULIA_SPECIALFUNCTIONS_BUILD_SOURCE", "false") == "true"

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get(filter(!isequal("verbose"), ARGS), 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, String["libopenspecfun"], :openspecfun),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaPackaging/Yggdrasil/releases/download/OpenSpecFun-v0.5.3+0"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.aarch64-linux-gnu-gcc4.tar.gz",
         "5d464227f3d45721747b21f71287614480d8b8048ae05feb44636f15b1fa0bfa"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.aarch64-linux-gnu-gcc7.tar.gz",
         "54e2c7f571ef991534aef07c0c5da1a52b29d775d4d608d2fe0b2936337a8ccf"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.aarch64-linux-gnu-gcc8.tar.gz",
         "260fa67e1903706f904445092280a984f9a2e0b97a3c5c9c289a961ddde3c9a6"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc4)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.aarch64-linux-musl-gcc4.tar.gz",
         "c18eb4a5140780ea78ac335c5c103617415fb5b009ab145c0d23831a8d52ed32"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc7)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.aarch64-linux-musl-gcc7.tar.gz",
         "ee6b805784bff8e0f3fab4fec8fee7e42fc5301e11459c21ea78680ad3b14a6c"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc8)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.aarch64-linux-musl-gcc8.tar.gz",
         "d362508cf6ac755719d2ce92a514f3465f48e75f1009eb3a7a6e3ff22c246177"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.arm-linux-gnueabihf-gcc4.tar.gz",
         "c5a5942a3bad5e65ec46533d87df747cf25d8e5e8b756d4bec592cdef6e0b2ce"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.arm-linux-gnueabihf-gcc7.tar.gz",
         "198663ea65dfb9563fb94641dc0974996f96818f35a878e14728186fb3d5d6bc"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.arm-linux-gnueabihf-gcc8.tar.gz",
         "b293d95c6ed2d0b60efbd45a59499f13df2bf5424c2660ad32c529eae4e10d2a"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.arm-linux-musleabihf-gcc4.tar.gz",
         "e8c56849d2814b493dc688e39e91ed402ffab96bd21791e74befaf4449587561"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.arm-linux-musleabihf-gcc7.tar.gz",
         "3578d09f65c52e78cc721043774988616e59b70f68bdab7233171c95b07f86f2"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.arm-linux-musleabihf-gcc8.tar.gz",
         "55dc2d0317da8c72942a050fca22637b3f715989ed9babdffd5407e3ce85957b"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.i686-linux-gnu-gcc4.tar.gz",
         "b8a0b32c3f28eb663f7cf44fa047211781fa770bb6e63d9931c083293ab99528"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.i686-linux-gnu-gcc7.tar.gz",
         "e15e699e16a5371122f55db2689dcea7066816b63c0b18fd3f66e04a612a3e70"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.i686-linux-gnu-gcc8.tar.gz",
         "98fb54908383119682f9a1544fb3ebc4b9af212ec4cfb4f10fb4611256a35128"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc4)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.i686-linux-musl-gcc4.tar.gz",
         "4eed4d21100b2849991c8d93cbed89979c007ed99a7dffb2e19616e9970b9715"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc7)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.i686-linux-musl-gcc7.tar.gz",
         "f52b4b867e3b1d3f8f968e5cf24268d154f9017955bd5880f49bf389a251161e"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc8)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.i686-linux-musl-gcc8.tar.gz",
         "169d56320b7c5f0df8f74d6dd4303195eeb2700e8a3032379cbec4f58d755b2b"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc4)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.i686-w64-mingw32-gcc4.tar.gz",
         "ceb767e8f71d1094b16d9a866f83ffd2769a57532a4bd79fc1b7586062263f16"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc7)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.i686-w64-mingw32-gcc7.tar.gz",
         "4779208e2340000d0ca7039b07fcef3a5ae58326e98b60beb79c7532bdbc05d7"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc8)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.i686-w64-mingw32-gcc8.tar.gz",
         "965ed911f2479917ff0bffc43137854e370ee961c5c3d03dd63d38ea23b855b6"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.powerpc64le-linux-gnu-gcc4.tar.gz",
         "966687bd725dfb954f8f8b6fdaf597d57a21718eacecb9249ad17fed0575b424"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.powerpc64le-linux-gnu-gcc7.tar.gz",
         "b016d073e293ee11182760489856e40970c718cb422185e5f71e31a027aeb2d8"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.powerpc64le-linux-gnu-gcc8.tar.gz",
         "beabf985b77943f77ab4a3eff413e5e689945c49a1d3a5c0c8e1f11188f9ea7a"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc4)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.x86_64-apple-darwin14-gcc4.tar.gz",
         "2e3706c0ba8c5820d327859f5cf8ece7085e8cf7e6f744e2a0a25d913086aa79"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc7)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.x86_64-apple-darwin14-gcc7.tar.gz",
         "be099a699c445f3bb94ccd421e3674e04adb8c9a85cdc48ea2b3e29b85c6a7bb"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc8)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.x86_64-apple-darwin14-gcc8.tar.gz",
         "d7e546d83e3e501f52587ed27796bae6c53d483d947cd07af97e00f000634a92"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.x86_64-linux-gnu-gcc4.tar.gz",
         "80a34541ab787872f8c9685db9be8eeca30b7959f29223555d4cbfb5d6a4893b"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.x86_64-linux-gnu-gcc7.tar.gz",
         "870bb918dc7a39550829daefd507b53d9748c9039a52a072d7d535143cd6731c"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.x86_64-linux-gnu-gcc8.tar.gz",
         "c8838499cc2f16845f10697da72e38eb8ea222cea00fbe6d0f2c72b573c4a997"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc4)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.x86_64-linux-musl-gcc4.tar.gz",
         "2775ed11e0a79e3e6941de1d6aa7f8095090bc8ffcf986b8d588fe04d3151185"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc7)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.x86_64-linux-musl-gcc7.tar.gz",
         "bb803c5ae7c7fd55cba257d6d49ad21638a421ac92a08da7be7bff86eda9593e"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc8)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.x86_64-linux-musl-gcc8.tar.gz",
         "8525a0f6fb54c1903f8e127275772ea846c6a81baedfd618e6bb60eef2d14f9b"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc4)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.x86_64-unknown-freebsd11.1-gcc4.tar.gz",
         "145a4846031f52f0777678eacbda0727649fe96f52695f101d67802539bb206a"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc7)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.x86_64-unknown-freebsd11.1-gcc7.tar.gz",
         "fb5b344a33dd7dfac5ff9a0b8d27caf5d5e9613f63e07bfed96a8acf6d311ba3"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc8)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.x86_64-unknown-freebsd11.1-gcc8.tar.gz",
         "dd944fa05a3ffde46ba88ffefbf9de62a97dea3ce43aac44fd4f8e1a5df9261c"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc4)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.x86_64-w64-mingw32-gcc4.tar.gz",
         "5b5ec3851b4d0ef8c802b6333259a4939c045ab8478f6f8b671b2129f5a02a86"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc7)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.x86_64-w64-mingw32-gcc7.tar.gz",
         "248a39a24cb36c7be5f6a2d8470f7b43ccddb026d8813e48139ea0b3eb978c3b"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc8)) =>
        ("$bin_prefix/OpenSpecFun.v0.5.3.x86_64-w64-mingw32-gcc8.tar.gz",
         "316b7beec8b7962d9114b309133863a95f80c3c04e368ebb4db76240a9594c93"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(p->!satisfied(p; verbose=verbose), products)
to_download = choose_download(download_info, platform_key_abi())
if to_download !== nothing && !forcecompile
    if !isinstalled(to_download...; prefix=prefix)
        # Download and install binaries
        install(to_download...; prefix=prefix, force=true, verbose=verbose)
        unsatisfied = any(p->!satisfied(p; verbose=verbose), products)
    end
    if unsatisfied
        rm(joinpath(@__DIR__, "usr", "lib"); force=true, recursive=true)
    end
end
if unsatisfied || forcecompile
    include("scratch.jl")
else
    # Write out a deps.jl file that will contain mappings for our products
    write_deps_file(joinpath(@__DIR__, "deps.jl"), products; verbose=verbose)
end


text:
using CmdStan, Turing, ContinuousBenchmarks
using Mamba

include(joinpath(ContinuousBenchmarks.STAN_DATA_DIR, "lda-stan.data.jl"))
include(joinpath(ContinuousBenchmarks.STAN_MODELS_DIR, "lda-stan.model.jl"))

stan_model_name = "LDA"
# ldastan = Stanmodel(Sample(save_warmup=true), name=stan_model_name, model=ldastanmodel, nchains=1);
# To understand parameters, use: ?CmdStan.Static, ?CmdStan.Hmc
ldastan = Stanmodel(Sample(algorithm=CmdStan.Hmc(CmdStan.Static(0.05),CmdStan.diag_e(),0.005,0.0),
  save_warmup=true,adapt=CmdStan.Adapt(engaged=false)),
  num_samples=3000, num_warmup=0, thin=1,
  name=stan_model_name, model=ldastanmodel, nchains=1);

rc, lda_stan_sim = stan(ldastan, ldastandata, CmdStanDir=ContinuousBenchmarks.CMDSTAN_HOME, summary=false);
# lda_stan_sim.names
V = ldastandata[1]["V"]
K = ldastandata[1]["K"]
lda_stan_d_raw = Dict()
for i = 1:K, j = 1:V
  lda_stan_d_raw["phi[$i][$j]"] = lda_stan_sim[1001:end, ["phi.$i.$j"], :].value[:]
end

lda_stan_d = Dict()
for i = 1:K
  lda_stan_d["phi[$i]"] = mean([[lda_stan_d_raw["phi[$i][$k]"][j] for k = 1:V] for j = 1:(ldastan.num_samples-1000)])
end

lda_time = get_stan_time(stan_model_name)
println("CmdStan time: ", lda_time)


text:
using Tmatrix
import Zygote

# inputs
rx, rz = 1e-6, 1.3e-6
n_max = 1
k1_r = 1e7; k1_i = 0.0
k2_r = 1.5e7; k2_i = 1e3

calculate_Tmatrix_for_spheroid_SeparateRealImag(rx, rz, n_max, k1_r, k1_i, k2_r, k2_i; rotationally_symmetric=false, symmetric_about_plane_perpendicular_z=false);

function Tmatrix_spheroid_simple(rx, rz)
    return calculate_Tmatrix_for_spheroid_SeparateRealImag(rx, rz, n_max, k1_r, k1_i, k2_r, k2_i; rotationally_symmetric=false, symmetric_about_plane_perpendicular_z=false);
end

Zygote.jacobian(Tmatrix_spheroid_simple, rx, rz)


text:
"""
    a

Scale length of the Kuzmin-Toomre disk.
Set to `1.0` by default.
"""
const a = 1.0

"""
    b

Scale length of the Toomre spherical dark halo.
Set to `2.8` by default, see `Reddish et al. (2021)`.
"""
const b = 2.8

"""
    c

Scale length of the Toomre spherical bulge.
Set to `1/20` by default, see `Reddish et al. (2021)`.
"""
const c = 1/20

"""
    error_dis

Error threshold on the convergence distance when determining
physical eigenvalues of the system.
"""
const error_dis = 0.005

"""
    error_var

Error threshold on the convergence variance when determining
physical eigenvalues of the system.
"""
const error_var = 0.005

"""
    threshold_nb_egv

Minimal number of cutoff matrix eigenvalues that must be found
to apply the selection criterion.
"""
const threshold_nb_egv = 10

"""
    threshold_GR

Lower threshold under which we don't compute the growth rate of the system.
Set to `0.1` by default, see `Reddish et al. (2021)`.
"""
const threshold_GR = 0.1

max_stars_repo_name: KerwannTEP/LiRA

max_issues_repo_name: KerwannTEP/LiRA

max_forks_repo_name: KerwannTEP/LiRA


text:
module PhysicsEngine

export fullLibraryPath

macro fullLibraryPath()
  return "./src/renderer/libphysicsenginerenderer"
end

include("vector.jl")
include("particle.jl")
include("motion.jl")
include("menu.jl")
include("renderer.jl")
include("simulation.jl")
include("grid.jl")

end # module

max_stars_repo_name: jonathanbuchanan/PhysicsEngine

max_issues_repo_name: jonathanbuchanan/PhysicsEngine

max_forks_repo_name: jonathanbuchanan/PhysicsEngine