text:
using AirfoilPrep
#using Base.Test
using Test
using Xfoil

using CSV
using PyPlot
using JLD

fileLoc,_ = splitdir(@__FILE__)

# Specify plotting options
rc("figure", figsize=(4.5, 2.6))
rc("font", size=10.0)
rc("lines", linewidth=1.5)
rc("lines", markersize=3.0)
rc("legend", frameon=false)
rc("axes.spines", right=false, top=false)
rc("figure.subplot", left=0.18, bottom=0.18, top=0.97, right=0.92)
# rc("axes", color_cycle=["348ABD", "A60628", "009E73", "7A68A6", "D55E00", "CC79A7"])
color_cycle=["#348ABD", "#A60628", "#009E73", "#7A68A6", "#D55E00", "#CC79A7", "#000486","#700000","#006907","#4C0099"]

PyPlot.close("all")

modulepath,_ = splitdir(@__FILE__)
modulepath = modulepath*"/"

AP = AirfoilPrep

data_path = "data/"

#-------- TEST AIRFOILPREPY WRAPPER --------#
"Verify that 3D corrections (stall delay) to 2D airfoil
properties is done correctly. The stall delay is implemented in `airfoilprep.py`
as explained in Du, Z., & Selig, M. (1998), *A 3-D stall-delay model for
horizontal axis wind turbine performance prediction*, and is verified by
replicating the results Du reports in Fig. 5 over the FFA airfoil in 5WPX."
function verify_correction3D_2()

    Re = Int(0.5 * 10^5) # Reynolds number
    Vinf = 8.8     # (m/s) Wind speed
    RPM = 158
    D = 5.35        # (m) Rotor diameter
    TSR = (2pi * RPM / 60) * (D / 2) / (Vinf)

    fig1 = figure("correction3D")
    title("FFA 5WPX airfoil lift coefficient")
    xlabel(L"Angle of attack $\alpha (^\circ)$")
    ylabel(L"$C_l$")
    grid(true, color="0.8", linestyle="--")

    # 2D curve
    data = CSV.read(modulepath*data_path*"FFA-2Dcalc.csv";
    header=["angle", "Cl"], datarow=1)
    plot(data[:,1], data[:,2], "-k", label="2D calculated")

    polar = AP.Polar(Re, data[:,1], data[:,2],
                        zeros(length(data[:,1])),
                        zeros(length(data[:,1])),
                        Float64[], Float64[])

    # Hard code old solutions for error checking
    old_cl_cr_16 = [0.0299769, 0.0473358, 0.0641102, 0.0826366, 0.0996613, 0.116936, 0.132459, 0.149484, 0.166759, 0.183533, 0.199056, 0.216706, 0.23443, 0.2518, 0.266208, 0.283751, 0.304505, 0.328773, 0.350701, 0.369912, 0.390744, 0.411294, 0.429392, 0.447877, 0.465695, 0.485271, 0.500772, 0.517415, 0.542212, 0.565233, 0.583722, 0.599897, 0.619426, 0.634199, 0.653888, 0.673884, 0.695413, 0.716739, 0.744017, 0.772024, 0.798671, 0.823187, 0.847241, 0.873737, 0.886028, 0.884262, 0.878168, 0.873115, 0.866446, 0.860485, 0.848607, 0.838469, 0.826784, 0.815457, 0.804488, 0.794233, 0.783979, 0.772321, 0.760691, 0.750102, 0.759824, 0.773619, 0.785396]
    old_cl_cr_37 = [0.0337801, 0.0513294, 0.0677512, 0.0849978, 0.101537, 0.118195, 0.134027, 0.150567, 0.167224, 0.183646, 0.199479, 0.21546, 0.23513, 0.251588, 0.262161, 0.281975, 0.30213, 0.326092, 0.346911, 0.366866, 0.386748, 0.406807, 0.425557, 0.443637, 0.461402, 0.484257, 0.501784, 0.514734, 0.538474, 0.56623, 0.585163, 0.603007, 0.623283, 0.64132, 0.662523, 0.683871, 0.713161, 0.74411, 0.802981, 0.866457, 0.922474, 0.969816, 1.01609, 1.08993, 1.17158, 1.21422, 1.23692, 1.26182, 1.2851, 1.30616, 1.31761, 1.32563, 1.33206, 1.33866, 1.34542, 1.35253, 1.35963, 1.36693, 1.37509, 1.39523, 1.45718, 1.54739, 1.62453]
    old_cl = (old_cl_cr_16, old_cl_cr_37)

    correction3D_error_max = 0.0

    # 3D curves
    i = 1
    for (mrkr, style, c_over_r, r_over_R) in [
            ("^r","--r",16,55),("Pg","-.g",37,30)]

        data = CSV.read(modulepath*data_path*"FFA-3D0$c_over_r.csv";
        header=["angle", "Cl"], datarow=1)
        plot(data[:,1], data[:,2], mrkr,
        label="3D wind tunnel c/r=0.$c_over_r")

        newpolar = AP.correction3D(polar, r_over_R/100, c_over_r/100, TSR,
                                    alpha_linear_min=0, alpha_linear_max=7,
                                    alpha_max_corr=30)

        plot(AP.get_cl(newpolar)[1], AP.get_cl(newpolar)[2], style,
                label="3D-airfoilprep.py corrected c/r=0.$c_over_r")
        correction3D_error_max = max(correction3D_error_max, maximum(abs.(AP.get_cl(newpolar)[2] - old_cl[i])))

        i += 1
    end
    legend(loc="best")

    return correction3D_error_max
end

"Tests (but only verifies, since there is nothing to compare to) that the
extrapolation method from `airfoilprep.py` is being called correctly. This
method takes airfoil curves an extrapolates them all around 360 degrees.
NOTE: CDmax in that function is the drag coeff at alpha=90deg or
CDmax=1.11+0.18AR for a blade of aspect ratio AR<50."
function test_extrapolation(; alphas=[i for i in -10:1.0:20], iter=100,
                            verbose=true)

    Re=1 * 10^5
    CDmax = 1.3

    Vinf = 10     # (m/s) Wind speed
    RPM = 72
    D = 10        # (m) Rotor diameter
    TSR = (2pi * RPM / 60) * (D / 2) / (Vinf)
    r_over_R = 30
    c_over_r = 18

    data = CSV.read(modulepath*data_path*"FFA-2Dcalc.csv";
    header=["angle", "Cl"], datarow=1)

    polar = AP.Polar(Re, data[:,1], data[:,2],
                        zeros(length(data[:,1])),
                        zeros(length(data[:,1])),
                        Float64[], Float64[])

    newpolar = AP.correction3D(polar, r_over_R/100, c_over_r/100, TSR,
                                alpha_linear_min=0, alpha_linear_max=7,
                                alpha_max_corr=30)

    # Extrapolated polar
    # extrap_polar1 = AP.extrapolate(newpolar, CDmax)
    extrap_polar = AP.extrapolate(polar, CDmax)

    # AP.plot(extrap_polar1; cdpolar=false)
    AP.plot(extrap_polar; cdpolar=false)

    return extrap_polar
end

"""#-------- TEST NDtools Automated Table Generation from Input Function --------#

#Use S809 NREL airfoil:
#1) Verify Cl, Cd, Cm with data from airfoiltool.com
"""
function validateNDtools_from_Xfoil()

    #--- Load in Airfoiltools.com S809 Data ---#
    S809_Re2E5 = CSV.read(modulepath*data_path*"xf-s809-nr-200000.csv";
                            header=["alpha" ,"CL","CD","CDp","CM","Top_Xtr","Bot_Xtr"],delim = ",", datarow=13)

    S809_Re5E5 = CSV.read(modulepath*data_path*"xf-s809-nr-500000.csv";
                            header=["alpha" ,"CL","CD","CDp","CM","Top_Xtr","Bot_Xtr"],delim = ",", datarow=13)

    S809_Re1E6 = CSV.read(modulepath*data_path*"xf-s809-nr-1000000.csv";
                            header=["alpha" ,"CL","CD","CDp","CM","Top_Xtr","Bot_Xtr"],delim = ",", datarow=13)

    AirfoilToolsData = (S809_Re2E5,S809_Re5E5,S809_Re1E6)

    #Load in airfoil x and y points
    folder,_ = splitdir(@__FILE__)
    airfoil_file = joinpath(folder, "data", "S809.txt")
    headerlines = 2
    open(airfoil_file, "r") do f
        global  x = Array{Float64,1}(undef, 0)
        global  y = Array{Float64,1}(undef, 0)
        for (i, line) in enumerate(eachline(f))
            if i > headerlines
                x = append!(x, parse(Float64, split(chomp(line))[1]))
                y = append!(y, parse(Float64, split(chomp(line))[2]))
            else
            end
        end
    end

    # Define operating conditions
    aoas = collect(range(-15, stop=25, length=41))#linspace(-10,10,20)
    Res = [2E5, 3E5, 4E5, 5E5, 6E5, 7E5, 8E5, 9E5, 1E6]
    Ms = [0.0, 0.01]

    #Wapper function for my analysis code: Xfoil for this test
    println("Running Xfoil as Airfoil Data generator")
    function f(Re, M)
        # Note that aoas is inherited for more efficient xfoil operation
        cls, cds, cdps, cms, convs = Xfoil.xfoilsweep(x, y, aoas, Re;iter=100, npan=140, mach=M,
        percussive_maintenance=true, printdata=true, zeroinit=true, clmaxstop=true, clminstop=true)

        return cls, cds, cms, convs
    end

    var_input = (aoas, Res, Ms)
    var_names = ["aoa", "Re", "M"]
    response_names = ["cl", "cd", "cm", "convs"]

    #Since the version of Xfoil being used is more efficient if it handles the aoa sweep, we'll not generate a table with it yet.
    response_values = AirfoilPrep.genNDarray(f, response_names, var_input[2:end], var_names[2:end])
    JLD.save("$(fileLoc)/data/response_values.jld", "response_values", response_values)
    # response_values = JLD.load("$(fileLoc)/data/response_values.jld")
    # response_values = response_values["response_values"]

    # Reformat to get the ND array in the right format for my specific problem, ie
    # a table of responses lining up to aoa, re, mach
    cls = Array{Float64}(undef, length(aoas), length(Res), length(Ms))
    cds = similar(cls)
    cms = similar(cls)
    convs = similar(cls)

    for i = 1:length(aoas)
        for j = 1:length(Res)
            for k = 1:length(Ms)
                cls[i,j,k] = response_values[j,k][1][i]
                cds[i,j,k] = response_values[j,k][2][i]
                cms[i,j,k] = response_values[j,k][3][i]
                convs[i,j,k] = response_values[j,k][4][i]
            end
        end
    end

    #Put the response values into the format required by NDtools
    response_values2 = [cls, cds, cms, convs]
    NDtable = AirfoilPrep.TableND(response_values2, response_names, var_input, var_names)

    #Access the table example
    myindices = [1, 1, 1]
    cl = NDtable.response_values[1][myindices...] #Assumed cl to be first response

    # Spline the table
    #Warning, the airfoil data here has non-converged points and if they include Nans or INF, the spline will be corrupted
    splout_non_extrap = AirfoilPrep.SplineND_from_tableND(NDtable)

    #Plot the results
    Re_airfoiltools = [2E5, 5E5, 1E6]
    Re_xfoil_idx = [1, 4, 9]
    aoas21 = []
    aoas2 = []
    XfoilData_cl1 =[]
    XfoilData_cl =[]
    XfoilData_cd1 =[]
    XfoilData_cd =[]
    XfoilData_cm1 =[]
    XfoilData_cm =[]
    vars = []

    for i = 1:length(Re_airfoiltools) #length of the airfoiltools data
        for j = 1:length(aoas)
            if NDtable.response_values[4][j,i,1]==true
                vars = [aoas[j], Re_airfoiltools[i], 0.0]
                #Corrupted spline because of Inf in points, use direct response values
                push!(XfoilData_cl1, NDtable.response_values[1][j, Re_xfoil_idx[i], 1])#push!(XfoilData_cl1, AirfoilPrep.interpND(splout_non_extrap[1], vars))
                push!(XfoilData_cd1, NDtable.response_values[2][j, Re_xfoil_idx[i], 1])#push!(XfoilData_cd1, AirfoilPrep.interpND(splout_non_extrap[2], vars))
                push!(XfoilData_cm1, NDtable.response_values[3][j, Re_xfoil_idx[i], 1])#push!(XfoilData_cm1, AirfoilPrep.interpND(splout_non_extrap[3], vars))
                push!(aoas21, aoas[j])
            end
        end

        push!(XfoilData_cl, XfoilData_cl1)
        push!(XfoilData_cd, XfoilData_cd1)
        push!(XfoilData_cm, XfoilData_cm1)
        push!(aoas2, aoas21)
        XfoilData_cl1 =[]
        XfoilData_cd1 =[]
        XfoilData_cm1 =[]
        aoas21 = []
    end

    # Check the cl error from an original run
    NDSpl_maxerror_cl = 0.0
    OldXfoilData_s809_Re2E5 = [-0.473515, -0.648893, -0.528673, -0.492563, -0.357591, -0.209986, -0.110082, -0.0481071, 0.00891502, 0.113733, 0.17017, 0.235873, 0.327773, 0.41861, 0.527213, 0.633993, 0.746845, 0.87758, 0.977806, 0.954624, 0.926883]
    OldXfoilData_s809_Re5E5 = [-0.746883, -0.650123, -0.584852, -0.551076, -0.524587, -0.512261, -0.439827, -0.318644, -0.202955, -0.0871805, 0.0295255, 0.146385, 0.263179, 0.37841, 0.492835, 0.60978, 0.726073, 0.837511, 0.943375, 0.960604, 0.983048, 1.01456, 1.03991, 1.0668, 1.09374, 1.1144, 1.14403, 1.15379, 1.1609, 1.1431]
    OldXfoilData_s809_Re1E6 = [-0.809019, -0.867028, -0.75136, -0.688495, -0.653602, -0.622135, -0.57855, -0.519902, -0.448765, -0.327558, -0.208416, -0.0891631, 0.0295293, 0.148423, 0.267439, 0.386354, 0.50429, 0.621041, 0.736629, 0.851222, 0.914147, 0.97305, 1.02751, 1.07725, 1.12474, 1.17066, 1.20861, 1.23612, 1.2648, 1.29016, 1.29755, 1.31023]
    XfoilData_cl_old = (OldXfoilData_s809_Re2E5, OldXfoilData_s809_Re5E5, OldXfoilData_s809_Re1E6)

    NDSpl_maxerror_cd = 0.0
    Xfoil_cd_old_2E5 = [0.15251, 0.0622333, 0.0686069, 0.06014, 0.100844, 0.108144, 0.101768, 0.0943553, 0.0879524, 0.0834842, 0.0600548, 0.051901, 0.0521719, 0.0502692, 0.0501376, 0.0487363, 0.0441949, 0.0360625, 0.0285292, 0.0320092, 0.0485597]
    Xfoil_cd_old_5E5 = [0.173309, 0.0654769, 0.0615699, 0.0519822, 0.0404028, 0.0244429, 0.0181499, 0.0217155, 0.0206846, 0.0198172, 0.0198481, 0.0199731, 0.0201065, 0.0204247, 0.0208011, 0.0203382, 0.0186523, 0.0164196, 0.0140924, 0.0200996, 0.0269638, 0.0367127, 0.0482435, 0.0593937, 0.0709354, 0.0845088, 0.0982058, 0.116762, 0.138359, 0.168361]
    Xfoil_cd_old_1E6 = [0.150439, 0.112753, 0.0674232, 0.0618973, 0.052744, 0.0418405, 0.0318201, 0.0215296, 0.00964882, 0.0118767, 0.0118197, 0.0117866, 0.0116456, 0.011748, 0.0119901, 0.0121527, 0.0124006, 0.011345, 0.0105963, 0.00977793, 0.0139917, 0.0180614, 0.0220341, 0.0281537, 0.0365331, 0.045208, 0.054797, 0.0660075, 0.0779478, 0.0915751, 0.109754, 0.128328]
    XfoilData_cd_old = (Xfoil_cd_old_2E5, Xfoil_cd_old_5E5, Xfoil_cd_old_1E6)

    NDSpl_maxerror_cm = 0.0
    Xfoil_cm_old_2E5 = [-0.0531119, -0.0380655, -0.0304781, -0.0277283, -0.0161606, -0.0158703, -0.0182155, -0.0192237, -0.0197291, -0.0224647, -0.0268114, -0.0331048, -0.0337455, -0.0384099, -0.0386892, -0.0400551, -0.0411134, -0.0418757, -0.039125, -0.0217998, -0.0158545]
    Xfoil_cm_old_5E5 = [-0.0333839, -0.0260745, -0.0236179, -0.02154, -0.0204999, -0.0225656, -0.0292425, -0.0310045, -0.0339394, -0.0367941, -0.0393168, -0.0419707, -0.0446739, -0.0471477, -0.0492277, -0.0512713, -0.0527747, -0.0533284, -0.0526812, -0.0397054, -0.0306905, -0.0274788, -0.0246275, -0.0219109, -0.0195295, -0.0178726, -0.0178988, -0.0181141, -0.0201363, -0.0239134]
    Xfoil_cm_old_1E6 = [-0.033515, -0.037818, -0.0223279, -0.0210895, -0.0188183, -0.0178605, -0.0182186, -0.0223099, -0.0295553, -0.0321853, -0.0349487, -0.0377451, -0.0404905, -0.0431673, -0.0458499, -0.0484486, -0.0508638, -0.0527459, -0.0543582, -0.0557457, -0.0494517, -0.0430175, -0.0367772, -0.0326887, -0.0308393, -0.029174, -0.0271752, -0.0249823, -0.0235408, -0.0230162, -0.0229448, -0.024175]
    XfoilData_cm_old = (Xfoil_cm_old_2E5, Xfoil_cm_old_5E5, Xfoil_cm_old_1E6)

    #Plot the AirfoilTools Data vs the newly generated data and calculate the error
    PyPlot.figure("Verify_NDspline_cl")
    for i = 1:length(Re_airfoiltools)
        PyPlot.plot(aoas2[i], XfoilData_cl[i],"-", color = color_cycle[i], label = "NDtools Re $(round(Int, Re_airfoiltools[i]))")
        PyPlot.plot(AirfoilToolsData[i][:,1], AirfoilToolsData[i][:,2], "--", color = color_cycle[i], label = "Aifoiltools.com Re $(round(Int, Re_airfoiltools[i]))")
        NDSpl_maxerror_cl = max(NDSpl_maxerror_cl, maximum(abs.(XfoilData_cl[i] - XfoilData_cl_old[i])))
    end
    PyPlot.xlabel("AOA")
    PyPlot.ylabel("cl")
    PyPlot.legend(loc = "best")

    #Plot the AirfoilTools Data vs the newly generated data and calculate the error
    PyPlot.figure("Verify_NDspline_cd")
    for i = 1:length(Re_airfoiltools)
        PyPlot.plot(aoas2[i], XfoilData_cd[i], "-",color = color_cycle[i], label = "NDtools Re $(round(Int,Re_airfoiltools[i]))")
        PyPlot.plot(AirfoilToolsData[i][:,1], (AirfoilToolsData[i][:,3] + AirfoilToolsData[i][:,4]), "--", color = color_cycle[i], label = "Aifoiltools.com Re $(round(Int, Re_airfoiltools[i]))")
        NDSpl_maxerror_cd = max(NDSpl_maxerror_cd, maximum(abs.(XfoilData_cd[i] - XfoilData_cd_old[i])))
    end
    PyPlot.xlabel("AOA")
    PyPlot.ylabel("cd")
    PyPlot.legend(loc = "best")

    #Plot the AirfoilTools Data vs the newly generated data and calculate the error
    PyPlot.figure("Verify_NDspline_cm")
    for i = 1:length(Re_airfoiltools)
        PyPlot.plot(aoas2[i], XfoilData_cm[i], "-", color = color_cycle[i], label = "NDtools Re $(round(Int, Re_airfoiltools[i]))")
        PyPlot.plot(AirfoilToolsData[i][:,1], AirfoilToolsData[i][:,5], "--", color = color_cycle[i], label = "Aifoiltools.com Re $(round(Int, Re_airfoiltools[i]))")
        NDSpl_maxerror_cm = max(NDSpl_maxerror_cm, maximum(abs.(XfoilData_cm[i] - XfoilData_cm_old[i])))
    end
    PyPlot.xlabel("AOA")
    PyPlot.ylabel("cm")
    PyPlot.legend(loc = "best")

    return NDSpl_maxerror_cl, NDSpl_maxerror_cd, NDSpl_maxerror_cm, NDtable
end

"""#-------- TEST NDtools Integration with AirfoilPrep.py 3D correction and extrapolation--------#

#Use S809 NREL airfoil data from previous test:
#2) Verify 3d correction
#3) Verify extrapolation
"""
function verifyNDtable_extrap(NDtable)
    # Test airfoilpreppy on the ND table
    r_over_R = 0.1
    c_over_r = 0.3
    TSR = 10.0
    grid_alphas=[i for i in -180:1.0:180]

    # Includes convergence checking for airfoil data

    NDextrap3D_3Dtable = AirfoilPrep.NDTable_correction3D_extrap(NDtable, r_over_R, c_over_r, TSR;grid_alphas=grid_alphas)

    # Spline the new table
    splout_extrap = AirfoilPrep.SplineND_from_tableND(NDextrap3D_3Dtable)

    #Plot the output and get the error
    Re_airfoiltools = [2E5, 5E5, 1E6]
    aoas = collect(range(-15, stop=25, length=41))#linspace(-10, 10, 20)
    extrap_aoas = range(-180, stop=180, length=360)
    vars = []

    #Verify cl
    cl_re2E5 = [-6.93889e-18, 0.0712869, 0.142574, 0.213861, 0.285148, 0.356435, 0.427719, 0.49902, 0.570246, 0.641786, 0.712012, 0.704806, 0.642183, 0.609179, 0.599366, 0.587846, 0.574187, 0.564356, 0.558051, 0.553071, 0.548993, 0.546057, 0.544112, 0.542846, 0.542093, 0.541765, 0.541762, 0.54198, 0.542338, 0.542767, 0.543206, 0.5436, 0.543901, 0.544066, 0.544055, 0.543837, 0.543378, 0.542654, 0.541638, 0.54031, 0.53865, 0.53664, 0.534267, 0.531516, 0.528375, 0.524836, 0.520888, 0.516527, 0.511744, 0.506537, 0.500902, 0.494837, 0.488341, 0.481414, 0.474057, 0.466272, 0.458063, 0.449433, 0.440387, 0.43093, 0.42107, 0.410814, 0.400169, 0.389144, 0.377749, 0.365994, 0.35389, 0.341448, 0.328681, 0.3156, 0.30222, 0.288554, 0.274616, 0.260421, 0.245985, 0.231322, 0.21645, 0.201384, 0.186142, 0.170741, 0.155198, 0.139531, 0.123759, 0.107899, 0.0919706, 0.0759927, 0.0599814, 0.0439581, 0.0279509, 0.0119729, -0.00399037, -0.0199576, -0.03595, -0.0519703, -0.0679905, -0.0839865, -0.0999417, -0.115838, -0.131657, -0.147379, -0.162986, -0.17846, -0.193784, -0.20894, -0.223911, -0.238681, -0.253232, -0.26755, -0.281618, -0.295422, -0.308947, -0.322179, -0.335104, -0.347711, -0.359985, -0.371916, -0.383492, -0.394703, -0.405539, -0.415991, -0.42605, -0.43571, -0.444961, -0.4538, -0.46222, -0.470218, -0.477789, -0.484931, -0.491643, -0.497924, -0.503774, -0.509194, -0.514188, -0.51876, -0.522913, -0.526656, -0.529995, -0.532939, -0.5355, -0.53769, -0.539522, -0.541014, -0.542184, -0.543051, -0.543639, -0.543974, -0.544084, -0.544003, -0.543767, -0.543413, -0.542981, -0.542545, -0.542177, -0.541848, -0.541609, -0.541918, -0.542884, -0.543142, -0.543181, -0.548885, -0.557653, -0.547701, -0.536676, -0.607032, -0.761634, -0.890489, -0.936298, -0.914037, -0.842143, -0.739127, -0.624639, -0.544183, -0.462012, -0.341172, -0.234641, -0.14943, -0.0773951, 0.0109539, 0.0968394, 0.153127, 0.219456, 0.306987, 0.397432, 0.49195, 0.590919, 0.689044, 0.794008, 0.902318, 0.964385, 1.00225, 1.00745, 0.949038, 0.893275, 0.865697, 0.847402, 0.828633, 0.812808, 0.80164, 0.793416, 0.786908, 0.781998, 0.778572, 0.776292, 0.774866, 0.774118, 0.773902, 0.77407, 0.774494, 0.775068, 0.775698, 0.776303, 0.776807, 0.777147, 0.777264, 0.777106, 0.776627, 0.775787, 0.774548, 0.772878, 0.770746, 0.768128, 0.765, 0.761342, 0.757135, 0.752365, 0.747019, 0.741085, 0.734555, 0.72742, 0.719676, 0.711319, 0.702347, 0.692759, 0.682556, 0.67174, 0.660315, 0.648286, 0.635659, 0.622442, 0.608643, 0.594273, 0.579342, 0.563862, 0.547846, 0.531309, 0.514264, 0.49673, 0.478721, 0.460255, 0.441353, 0.422031, 0.402311, 0.382214, 0.36176, 0.340973, 0.319873, 0.298486, 0.276835, 0.254945, 0.232837, 0.210536, 0.188081, 0.165503, 0.14278, 0.119909, 0.0970974, 0.0745006, 0.0515093, 0.0275967, 0.00501722, -0.0131631, -0.0282698, -0.0436438, -0.0598812, -0.0760753, -0.0920015, -0.107877, -0.123749, -0.139537, -0.155201, -0.170739, -0.186141, -0.201385, -0.21645, -0.231322, -0.245985, -0.260421, -0.274616, -0.288554, -0.30222, -0.3156, -0.328681, -0.341448, -0.35389, -0.365994, -0.377749, -0.389144, -0.400169, -0.410814, -0.42107, -0.43093, -0.440387, -0.449433, -0.458063, -0.466272, -0.474057, -0.481414, -0.488341, -0.494837, -0.500902, -0.506537, -0.511744, -0.516527, -0.520888, -0.524836, -0.528375, -0.531516, -0.534267, -0.53664, -0.53865, -0.54031, -0.541638, -0.542654, -0.543378, -0.543837, -0.544055, -0.544066, -0.543901, -0.5436, -0.543206, -0.542767, -0.542338, -0.54198, -0.541762, -0.541765, -0.542093, -0.542846, -0.544112, -0.546057, -0.548993, -0.553071, -0.558051, -0.564356, -0.574187, -0.587846, -0.599366, -0.609179, -0.642183, -0.704806, -0.712012, -0.641786, -0.570246, -0.49902, -0.427719, -0.356435, -0.285148, -0.213861, -0.142574, -0.0712869, 1.04083e-17]
    cl_re5E5 = [-1.73472e-18, 0.0603834, 0.120767, 0.18115, 0.241534, 0.301917, 0.362301, 0.422684, 0.483068, 0.543451, 0.603834, 0.664218, 0.724601, 0.784984, 0.84537, 0.905744, 0.966165, 1.02641, 1.08741, 1.14542, 1.12358, 1.0635, 1.03199, 1.00947, 0.982833, 0.958005, 0.936738, 0.916978, 0.898289, 0.880918, 0.864672, 0.849298, 0.834684, 0.820734, 0.807338, 0.794398, 0.781836, 0.769581, 0.757565, 0.745733, 0.73403, 0.722413, 0.710837, 0.699268, 0.687671, 0.676017, 0.66428, 0.652436, 0.640465, 0.628349, 0.616073, 0.603623, 0.590988, 0.578158, 0.565126, 0.551887, 0.538435, 0.524768, 0.510884, 0.496783, 0.482467, 0.467937, 0.453198, 0.438253, 0.423109, 0.407771, 0.392248, 0.376548, 0.360679, 0.344653, 0.328479, 0.31217, 0.295737, 0.279194, 0.262553, 0.24583, 0.229038, 0.212192, 0.195309, 0.178405, 0.161496, 0.144598, 0.127731, 0.110909, 0.0941532, 0.0774827, 0.0609107, 0.0444523, 0.0281522, 0.0120287, -0.00400706, -0.0200726, -0.0362784, -0.0526653, -0.069184, -0.0858062, -0.102522, -0.119313, -0.13616, -0.153044, -0.16995, -0.186859, -0.203754, -0.220621, -0.237441, -0.254201, -0.270885, -0.287478, -0.303968, -0.320341, -0.336584, -0.352685, -0.368634, -0.38442, -0.400033, -0.415464, -0.430706, -0.445751, -0.460594, -0.475228, -0.489652, -0.50386, -0.517853, -0.531628, -0.545187, -0.558533, -0.571668, -0.584598, -0.597329, -0.60987, -0.622232, -0.634426, -0.646467, -0.658372, -0.67016, -0.681853, -0.693475, -0.705054, -0.716622, -0.728214, -0.739868, -0.751629, -0.763547, -0.775674, -0.788074, -0.800815, -0.813973, -0.827634, -0.841901, -0.856881, -0.872678, -0.889455, -0.907451, -0.926703, -0.947199, -0.969969, -0.995742, -1.02129, -1.04708, -1.08992, -1.13516, -1.14114, -1.12589, -1.11455, -1.10236, -1.08919, -1.06718, -1.02653, -0.958479, -0.855646, -0.752121, -0.677341, -0.610599, -0.553211, -0.484306, -0.376134, -0.257181, -0.142558, -0.0265601, 0.0895154, 0.205717, 0.321438, 0.436201, 0.551542, 0.66826, 0.781476, 0.897821, 0.980528, 1.03529, 1.10177, 1.16687, 1.23023, 1.29534, 1.35677, 1.42126, 1.48105, 1.53539, 1.57866, 1.62264, 1.62179, 1.55562, 1.49669, 1.45901, 1.42233, 1.38567, 1.35318, 1.32385, 1.29635, 1.27065, 1.24668, 1.22411, 1.20272, 1.18233, 1.16282, 1.14402, 1.12582, 1.10811, 1.09078, 1.07376, 1.05695, 1.04031, 1.02375, 1.00722, 0.990678, 0.974075, 0.957372, 0.940532, 0.923525, 0.906323, 0.888903, 0.871243, 0.853327, 0.835139, 0.816668, 0.797904, 0.778839, 0.759469, 0.73979, 0.719801, 0.699503, 0.678898, 0.657991, 0.636787, 0.615294, 0.59352, 0.571475, 0.549171, 0.52662, 0.503836, 0.480834, 0.45763, 0.43424, 0.410684, 0.386978, 0.363144, 0.339202, 0.315172, 0.291078, 0.266941, 0.242786, 0.218635, 0.194509, 0.170448, 0.14648, 0.122563, 0.0987549, 0.075357, 0.0520911, 0.0279976, 0.0050607, -0.0130771, -0.0282399, -0.0441876, -0.060929, -0.0775417, -0.0941403, -0.110897, -0.127735, -0.1446, -0.161494, -0.178405, -0.19531, -0.212192, -0.229037, -0.24583, -0.262553, -0.279194, -0.295737, -0.31217, -0.328479, -0.344653, -0.360679, -0.376548, -0.392248, -0.407771, -0.423109, -0.438253, -0.453198, -0.467937, -0.482467, -0.496783, -0.510884, -0.524768, -0.538435, -0.551887, -0.565126, -0.578158, -0.590988, -0.603623, -0.616073, -0.628349, -0.640465, -0.652436, -0.66428, -0.676017, -0.687671, -0.699268, -0.710837, -0.722413, -0.73403, -0.745733, -0.757565, -0.769581, -0.781836, -0.794398, -0.807338, -0.820734, -0.834684, -0.849298, -0.864672, -0.880918, -0.898289, -0.916978, -0.936738, -0.958005, -0.982833, -1.00947, -1.03199, -1.0635, -1.12358, -1.14542, -1.08741, -1.02641, -0.966165, -0.905744, -0.84537, -0.784984, -0.724601, -0.664218, -0.603834, -0.543451, -0.483068, -0.422684, -0.362301, -0.301917, -0.241534, -0.18115, -0.120767, -0.0603834, 0.0]
    cl_re1E6 = [0.0, 0.0613656, 0.122731, 0.184097, 0.245462, 0.306828, 0.368194, 0.429559, 0.490925, 0.55229, 0.613656, 0.675022, 0.736387, 0.797753, 0.859118, 0.920486, 0.981841, 1.04325, 1.10444, 1.16643, 1.22523, 1.20018, 1.13713, 1.1036, 1.0781, 1.04845, 1.02105, 0.997171, 0.974721, 0.953443, 0.933561, 0.914837, 0.897034, 0.88005, 0.863783, 0.848117, 0.832961, 0.818232, 0.803859, 0.789774, 0.775919, 0.762241, 0.748694, 0.735233, 0.721823, 0.70843, 0.695022, 0.681574, 0.668062, 0.654465, 0.640765, 0.626946, 0.612995, 0.5989, 0.584651, 0.570242, 0.555666, 0.540919, 0.525998, 0.510902, 0.49563, 0.480184, 0.464567, 0.448782, 0.432834, 0.416728, 0.400472, 0.384073, 0.36754, 0.350882, 0.334109, 0.317233, 0.300266, 0.283219, 0.266105, 0.24894, 0.231736, 0.214509, 0.197275, 0.180048, 0.162846, 0.145685, 0.128582, 0.111555, 0.0946211, 0.0778019, 0.06111, 0.0445585, 0.0281951, 0.01204, -0.00401035, -0.0200964, -0.0363487, -0.0528147, -0.0694398, -0.0861961, -0.103075, -0.120058, -0.137125, -0.154259, -0.171443, -0.188659, -0.205892, -0.223125, -0.240342, -0.257528, -0.274669, -0.291751, -0.30876, -0.325683, -0.342509, -0.359226, -0.375823, -0.39229, -0.408619, -0.4248, -0.440828, -0.456695, -0.472397, -0.487929, -0.503288, -0.518472, -0.53348, -0.548314, -0.562975, -0.577467, -0.591795, -0.605966, -0.619988, -0.633871, -0.647629, -0.661275, -0.674827, -0.688305, -0.701729, -0.715126, -0.728524, -0.741955, -0.755454, -0.769061, -0.782821, -0.796784, -0.811006, -0.825548, -0.840481, -0.855881, -0.871835, -0.888446, -0.905826, -0.924074, -0.943344, -0.963902, -0.985795, -1.00888, -1.03428, -1.06312, -1.09133, -1.11851, -1.16638, -1.2172, -1.21542, -1.18445, -1.15879, -1.13218, -1.10407, -1.08385, -1.05639, -0.993347, -0.910789, -0.824305, -0.748475, -0.683404, -0.614931, -0.537761, -0.455748, -0.354322, -0.237352, -0.124335, -0.0104386, 0.103048, 0.216744, 0.330471, 0.443843, 0.556827, 0.66822, 0.782547, 0.879329, 0.958737, 1.03914, 1.11602, 1.19116, 1.2656, 1.33769, 1.40345, 1.46809, 1.53279, 1.58813, 1.64362, 1.69128, 1.73758, 1.73359, 1.66296, 1.60054, 1.55913, 1.51818, 1.47761, 1.44138, 1.40824, 1.37698, 1.34765, 1.32011, 1.29403, 1.26921, 1.24548, 1.22269, 1.20069, 1.17935, 1.15858, 1.13826, 1.11832, 1.09866, 1.07922, 1.05994, 1.04075, 1.02161, 1.00247, 0.983292, 0.964039, 0.944679, 0.925184, 0.90553, 0.885696, 0.865665, 0.845422, 0.824953, 0.80425, 0.783306, 0.762115, 0.740674, 0.718982, 0.697041, 0.674853, 0.652422, 0.629754, 0.606858, 0.583741, 0.560414, 0.53689, 0.51318, 0.489299, 0.465262, 0.441086, 0.416787, 0.392385, 0.367898, 0.343346, 0.31875, 0.294131, 0.269513, 0.244919, 0.220371, 0.195889, 0.17151, 0.147269, 0.123127, 0.0991208, 0.0755543, 0.0521978, 0.0280585, 0.00506784, -0.013071, -0.028261, -0.0443046, -0.0611374, -0.077855, -0.094606, -0.111545, -0.128587, -0.145686, -0.162844, -0.180048, -0.197275, -0.214509, -0.231736, -0.24894, -0.266105, -0.283219, -0.300266, -0.317233, -0.334109, -0.350882, -0.36754, -0.384073, -0.400472, -0.416728, -0.432834, -0.448782, -0.464567, -0.480184, -0.49563, -0.510902, -0.525998, -0.540919, -0.555666, -0.570242, -0.584651, -0.5989, -0.612995, -0.626946, -0.640765, -0.654465, -0.668062, -0.681574, -0.695022, -0.70843, -0.721823, -0.735233, -0.748694, -0.762241, -0.775919, -0.789774, -0.803859, -0.818232, -0.832961, -0.848117, -0.863783, -0.88005, -0.897034, -0.914837, -0.933561, -0.953443, -0.974721, -0.997171, -1.02105, -1.04845, -1.0781, -1.1036, -1.13713, -1.20018, -1.22523, -1.16643, -1.10444, -1.04325, -0.981841, -0.920486, -0.859118, -0.797753, -0.736387, -0.675022, -0.613656, -0.55229, -0.490925, -0.429559, -0.368194, -0.306828, -0.245462, -0.184097, -0.122731, -0.0613656, 0.0]

    extrap_cl_old = (cl_re2E5, cl_re5E5, cl_re1E6)
    ND_corr3Dextr_maxerror_cl = 0.0

    #Plot the output and get the error
    extrap_cl = Array{Float64}(undef, length(extrap_aoas), length(Re_airfoiltools))
    PyPlot.figure("Verify_cl")
    for i = 1:length(Re_airfoiltools) #length of the airfoiltools data
        for j = 1:length(extrap_aoas)
            vars = [extrap_aoas[j], Re_airfoiltools[i], 0.0]
            extrap_cl[j,i] = AirfoilPrep.interpND(splout_extrap[1], vars)
        end

        PyPlot.plot(extrap_aoas,extrap_cl[:,i], label = "Re: $(Re_airfoiltools[i])")
        # println(extrap_cl[:,i])
        ND_corr3Dextr_maxerror_cl = max(ND_corr3Dextr_maxerror_cl, maximum(abs.(extrap_cl[:,i] - extrap_cl_old[i])))
    end
    PyPlot.xlabel("aoa (deg)")
    PyPlot.ylabel("cl")

    # Verify cd

    cd_re2E5 = [0.0144709, 0.0148668, 0.0160543, 0.0180316, 0.0207965, 0.0243454, 0.0286742, 0.0337774, 0.0396487, 0.046281, 0.0536661, 0.0617949, 0.0706573, 0.0802425, 0.0905387, 0.101533, 0.113212, 0.125562, 0.138567, 0.152211, 0.166478, 0.181349, 0.196807, 0.212833, 0.229406, 0.246507, 0.264115, 0.282206, 0.300761, 0.319755, 0.339165, 0.358968, 0.379138, 0.399652, 0.420483, 0.441607, 0.462996, 0.484626, 0.506469, 0.528498, 0.550686, 0.573007, 0.595431, 0.617932, 0.640483, 0.663054, 0.685619, 0.708149, 0.730618, 0.752996, 0.775257, 0.797373, 0.819316, 0.841061, 0.862579, 0.883845, 0.904832, 0.925514, 0.945865, 0.965861, 0.985477, 1.00469, 1.02347, 1.0418, 1.05966, 1.07702, 1.09387, 1.11018, 1.12592, 1.14109, 1.15566, 1.16962, 1.18295, 1.19562, 1.20763, 1.21897, 1.2296, 1.23954, 1.24875, 1.25723, 1.26497, 1.27196, 1.27819, 1.28365, 1.28833, 1.29223, 1.29536, 1.29771, 1.2992, 1.29987, 1.29999, 1.29964, 1.29857, 1.29663, 1.29389, 1.29038, 1.28609, 1.28101, 1.27517, 1.26856, 1.26119, 1.25308, 1.24423, 1.23466, 1.22437, 1.21339, 1.20171, 1.18937, 1.17636, 1.16272, 1.14845, 1.13358, 1.11812, 1.10209, 1.08551, 1.06841, 1.05079, 1.0327, 1.01414, 0.995135, 0.975718, 0.955909, 0.935732, 0.915212, 0.894375, 0.873245, 0.85185, 0.830215, 0.808368, 0.786335, 0.764143, 0.74182, 0.719393, 0.69689, 0.674339, 0.651767, 0.629203, 0.606674, 0.584208, 0.561832, 0.539574, 0.517462, 0.495522, 0.473783, 0.45227, 0.43101, 0.410029, 0.389353, 0.36901, 0.34902, 0.329403, 0.310203, 0.291453, 0.273101, 0.255145, 0.237921, 0.221456, 0.204535, 0.187472, 0.175069, 0.165009, 0.138946, 0.111359, 0.150791, 0.254135, 0.31801, 0.297911, 0.227085, 0.142559, 0.0814296, 0.0770119, 0.0711075, 0.0831208, 0.123413, 0.122508, 0.112931, 0.100453, 0.0941793, 0.0859922, 0.0673867, 0.0522863, 0.0493302, 0.0495132, 0.0488907, 0.049421, 0.0468843, 0.0404307, 0.030929, 0.0280917, 0.0434035, 0.0611058, 0.068031, 0.0744847, 0.0847454, 0.0961685, 0.10745, 0.119251, 0.131937, 0.145323, 0.159281, 0.173836, 0.189002, 0.204751, 0.221054, 0.237892, 0.255249, 0.273101, 0.291427, 0.310204, 0.329409, 0.349019, 0.369008, 0.389353, 0.410029, 0.43101, 0.45227, 0.473783, 0.495522, 0.517462, 0.539574, 0.561832, 0.584208, 0.606674, 0.629203, 0.651767, 0.674339, 0.69689, 0.719393, 0.74182, 0.764143, 0.786335, 0.808368, 0.830215, 0.85185, 0.873245, 0.894375, 0.915212, 0.935732, 0.955909, 0.975718, 0.995135, 1.01414, 1.0327, 1.05079, 1.06841, 1.08551, 1.10209, 1.11812, 1.13358, 1.14845, 1.16272, 1.17636, 1.18937, 1.20171, 1.21339, 1.22437, 1.23466, 1.24423, 1.25308, 1.26119, 1.26856, 1.27517, 1.28101, 1.28609, 1.29038, 1.29389, 1.29664, 1.29857, 1.29963, 1.29999, 1.29988, 1.2992, 1.29771, 1.29536, 1.29223, 1.28833, 1.28365, 1.27819, 1.27196, 1.26497, 1.25723, 1.24875, 1.23954, 1.2296, 1.21897, 1.20763, 1.19562, 1.18295, 1.16962, 1.15566, 1.14109, 1.12592, 1.11018, 1.09387, 1.07702, 1.05966, 1.0418, 1.02347, 1.00469, 0.985477, 0.965861, 0.945865, 0.925514, 0.904832, 0.883845, 0.862579, 0.841061, 0.819316, 0.797373, 0.775257, 0.752996, 0.730618, 0.708149, 0.685619, 0.663054, 0.640483, 0.617932, 0.595431, 0.573007, 0.550686, 0.528498, 0.506469, 0.484626, 0.462996, 0.441607, 0.420483, 0.399652, 0.379138, 0.358968, 0.339165, 0.319755, 0.300761, 0.282206, 0.264115, 0.246507, 0.229406, 0.212833, 0.196807, 0.181349, 0.166478, 0.152211, 0.138567, 0.125562, 0.113212, 0.101533, 0.0905387, 0.0802425, 0.0706573, 0.0617949, 0.0536661, 0.046281, 0.0396487, 0.0337774, 0.0286742, 0.0243454, 0.0207965, 0.0180316, 0.0160543, 0.0148668, 0.0144709]
    cd_re5E5 = [0.186304, 0.186674, 0.187782, 0.189628, 0.192209, 0.195521, 0.199561, 0.204323, 0.2098, 0.215987, 0.222874, 0.230454, 0.238715, 0.247647, 0.257239, 0.267479, 0.278352, 0.289846, 0.301944, 0.314631, 0.327891, 0.341707, 0.35606, 0.370931, 0.386302, 0.402153, 0.418462, 0.435209, 0.452371, 0.469926, 0.487852, 0.506124, 0.524719, 0.543613, 0.56278, 0.582196, 0.601835, 0.621672, 0.64168, 0.661832, 0.682103, 0.702466, 0.722893, 0.743358, 0.763834, 0.784293, 0.804708, 0.825052, 0.845298, 0.865419, 0.885389, 0.905179, 0.924765, 0.944119, 0.963215, 0.982028, 1.00053, 1.0187, 1.03651, 1.05394, 1.07096, 1.08755, 1.10369, 1.11935, 1.13451, 1.14915, 1.16326, 1.1768, 1.18977, 1.20213, 1.21388, 1.225, 1.23547, 1.24528, 1.2544, 1.26283, 1.27056, 1.27756, 1.28384, 1.28937, 1.29415, 1.29817, 1.30142, 1.30391, 1.30561, 1.30645, 1.30661, 1.30621, 1.3044, 1.30123, 1.30016, 1.30281, 1.30557, 1.3065, 1.30659, 1.30613, 1.30487, 1.30276, 1.29989, 1.29625, 1.29185, 1.2867, 1.28079, 1.27415, 1.26679, 1.25871, 1.24993, 1.24046, 1.23032, 1.21952, 1.20809, 1.19602, 1.18336, 1.1701, 1.15627, 1.1419, 1.12699, 1.11158, 1.09568, 1.07931, 1.0625, 1.04527, 1.02765, 1.00966, 0.99132, 0.972658, 0.953701, 0.934472, 0.914999, 0.895308, 0.875425, 0.855376, 0.835189, 0.81489, 0.794507, 0.774067, 0.753596, 0.733123, 0.712673, 0.692275, 0.671955, 0.65174, 0.631656, 0.611731, 0.59199, 0.572459, 0.553164, 0.534131, 0.515382, 0.496944, 0.478846, 0.461105, 0.443728, 0.42677, 0.410309, 0.39418, 0.378302, 0.363532, 0.349881, 0.333499, 0.318073, 0.313197, 0.313937, 0.312951, 0.314932, 0.303717, 0.254403, 0.183532, 0.125173, 0.108017, 0.106282, 0.0874062, 0.0660878, 0.040154, 0.0192421, 0.0183175, 0.0210847, 0.0194241, 0.0193908, 0.0196866, 0.0199539, 0.0202434, 0.0206622, 0.0207652, 0.0196164, 0.0177589, 0.0145042, 0.0168699, 0.0255201, 0.036466, 0.0516947, 0.0687207, 0.0865783, 0.107853, 0.130738, 0.158357, 0.192239, 0.233799, 0.290907, 0.329778, 0.337402, 0.346336, 0.363171, 0.379171, 0.394141, 0.410113, 0.42681, 0.443767, 0.46109, 0.478839, 0.496949, 0.515383, 0.534129, 0.553164, 0.572459, 0.59199, 0.611731, 0.631656, 0.65174, 0.671955, 0.692275, 0.712673, 0.733123, 0.753596, 0.774067, 0.794507, 0.81489, 0.835189, 0.855376, 0.875425, 0.895308, 0.914999, 0.934472, 0.953701, 0.972658, 0.99132, 1.00966, 1.02765, 1.04527, 1.0625, 1.07931, 1.09568, 1.11158, 1.12699, 1.1419, 1.15627, 1.1701, 1.18336, 1.19602, 1.20809, 1.21952, 1.23032, 1.24046, 1.24993, 1.25871, 1.26679, 1.27415, 1.28079, 1.2867, 1.29185, 1.29625, 1.29989, 1.30276, 1.30487, 1.30613, 1.30659, 1.3065, 1.30557, 1.30281, 1.30016, 1.30123, 1.3044, 1.30621, 1.30661, 1.30645, 1.30561, 1.30391, 1.30142, 1.29817, 1.29415, 1.28937, 1.28384, 1.27756, 1.27056, 1.26283, 1.2544, 1.24528, 1.23547, 1.225, 1.21388, 1.20213, 1.18977, 1.1768, 1.16326, 1.14915, 1.13451, 1.11935, 1.10369, 1.08755, 1.07096, 1.05394, 1.03651, 1.0187, 1.00053, 0.982028, 0.963215, 0.944119, 0.924765, 0.905179, 0.885389, 0.865419, 0.845298, 0.825052, 0.804708, 0.784293, 0.763834, 0.743358, 0.722893, 0.702466, 0.682103, 0.661832, 0.64168, 0.621672, 0.601835, 0.582196, 0.56278, 0.543613, 0.524719, 0.506124, 0.487852, 0.469926, 0.452371, 0.435209, 0.418462, 0.402153, 0.386302, 0.370931, 0.35606, 0.341707, 0.327891, 0.314631, 0.301944, 0.289846, 0.278352, 0.267479, 0.257239, 0.247647, 0.238715, 0.230454, 0.222874, 0.215987, 0.2098, 0.204323, 0.199561, 0.195521, 0.192209, 0.189628, 0.187782, 0.186674, 0.186304]
    cd_re1E6 = [0.153083, 0.153458, 0.154582, 0.156453, 0.159069, 0.162427, 0.166523, 0.171351, 0.176905, 0.183177, 0.190161, 0.197846, 0.206224, 0.215283, 0.225011, 0.235396, 0.246426, 0.258084, 0.270358, 0.28323, 0.296685, 0.310704, 0.325271, 0.340366, 0.355969, 0.372062, 0.388622, 0.405628, 0.42306, 0.440893, 0.459106, 0.477674, 0.496574, 0.51578, 0.53527, 0.555016, 0.574993, 0.595177, 0.615539, 0.636055, 0.656696, 0.677437, 0.698251, 0.719109, 0.739986, 0.760853, 0.781684, 0.802451, 0.823127, 0.843684, 0.864097, 0.884337, 0.904378, 0.924194, 0.943759, 0.963046, 0.982029, 1.00068, 1.01899, 1.03691, 1.05443, 1.07153, 1.08818, 1.10435, 1.12004, 1.13521, 1.14984, 1.16392, 1.17742, 1.19033, 1.20263, 1.2143, 1.22532, 1.23568, 1.24536, 1.25435, 1.26264, 1.27021, 1.27705, 1.28315, 1.28851, 1.2931, 1.29693, 1.29999, 1.30227, 1.30371, 1.30443, 1.30455, 1.30342, 1.30098, 1.30013, 1.30222, 1.30422, 1.30458, 1.30414, 1.30309, 1.30123, 1.29855, 1.29511, 1.2909, 1.28593, 1.2802, 1.27372, 1.26652, 1.25859, 1.24994, 1.24061, 1.23058, 1.21989, 1.20854, 1.19656, 1.18395, 1.17074, 1.15695, 1.14259, 1.12769, 1.11226, 1.09633, 1.07991, 1.06304, 1.04572, 1.028, 1.00988, 0.9914, 0.972577, 0.953439, 0.93401, 0.914316, 0.894384, 0.87424, 0.85391, 0.833422, 0.812802, 0.792077, 0.771275, 0.750422, 0.729547, 0.708676, 0.687837, 0.667056, 0.646361, 0.625779, 0.605337, 0.585061, 0.564977, 0.545112, 0.525491, 0.50614, 0.487083, 0.468347, 0.449956, 0.431923, 0.414284, 0.397098, 0.380299, 0.363824, 0.348135, 0.333333, 0.317296, 0.301905, 0.292678, 0.286496, 0.280298, 0.272097, 0.271958, 0.235097, 0.160574, 0.123544, 0.10981, 0.0981588, 0.0811789, 0.0596421, 0.0398281, 0.024133, 0.00724492, 0.00188064, 0.00676104, 0.00746499, 0.00885054, 0.0098063, 0.010922, 0.0117597, 0.0125712, 0.0123365, 0.0114936, 0.010155, 0.011649, 0.0164234, 0.0210747, 0.0278185, 0.0381566, 0.0506324, 0.0645944, 0.0818512, 0.101146, 0.122921, 0.151598, 0.183467, 0.221009, 0.274224, 0.311358, 0.32013, 0.330536, 0.347971, 0.36448, 0.380227, 0.396962, 0.414325, 0.431946, 0.449942, 0.468344, 0.487087, 0.50614, 0.525491, 0.545112, 0.564977, 0.585061, 0.605337, 0.625779, 0.646361, 0.667056, 0.687837, 0.708676, 0.729547, 0.750422, 0.771275, 0.792077, 0.812802, 0.833422, 0.85391, 0.87424, 0.894384, 0.914316, 0.93401, 0.953439, 0.972577, 0.9914, 1.00988, 1.028, 1.04572, 1.06304, 1.07991, 1.09633, 1.11226, 1.12769, 1.14259, 1.15695, 1.17074, 1.18395, 1.19656, 1.20854, 1.21989, 1.23058, 1.24061, 1.24994, 1.25859, 1.26652, 1.27372, 1.2802, 1.28593, 1.2909, 1.29511, 1.29855, 1.30123, 1.30309, 1.30414, 1.30458, 1.30422, 1.30222, 1.30013, 1.30098, 1.30342, 1.30455, 1.30443, 1.30371, 1.30227, 1.29999, 1.29693, 1.2931, 1.28851, 1.28315, 1.27705, 1.27021, 1.26264, 1.25435, 1.24536, 1.23568, 1.22532, 1.2143, 1.20263, 1.19033, 1.17742, 1.16392, 1.14984, 1.13521, 1.12004, 1.10435, 1.08818, 1.07153, 1.05443, 1.03691, 1.01899, 1.00068, 0.982029, 0.963046, 0.943759, 0.924194, 0.904378, 0.884337, 0.864097, 0.843684, 0.823127, 0.802451, 0.781684, 0.760853, 0.739986, 0.719109, 0.698251, 0.677437, 0.656696, 0.636055, 0.615539, 0.595177, 0.574993, 0.555016, 0.53527, 0.51578, 0.496574, 0.477674, 0.459106, 0.440893, 0.42306, 0.405628, 0.388622, 0.372062, 0.355969, 0.340366, 0.325271, 0.310704, 0.296685, 0.28323, 0.270358, 0.258084, 0.246426, 0.235396, 0.225011, 0.215283, 0.206224, 0.197846, 0.190161, 0.183177, 0.176905, 0.171351, 0.166523, 0.162427, 0.159069, 0.156453, 0.154582, 0.153458, 0.153083]

    extrap_cd_old = (cd_re2E5, cd_re5E5, cd_re1E6)
    ND_corr3Dextr_maxerror_cd = 0.0

    #Plot the output and get the error
    extrap_cd = Array{Float64}(undef, length(extrap_aoas), length(Re_airfoiltools))
    PyPlot.figure("Verify_cd")
    for i = 1:length(Re_airfoiltools) #length of the airfoiltools data
        for j = 1:length(extrap_aoas)
            vars = [extrap_aoas[j], Re_airfoiltools[i], 0.0]
            extrap_cd[j,i] = AirfoilPrep.interpND(splout_extrap[2], vars)
        end

        # println("cd")
        # println(extrap_cd[:,i])
        PyPlot.plot(extrap_aoas, extrap_cd[:,i], label = "Re: $(Re_airfoiltools[i])")
        ND_corr3Dextr_maxerror_cd = max(ND_corr3Dextr_maxerror_cd, maximum(abs.(extrap_cd[:,i] - extrap_cd_old[i])))
    end
    PyPlot.xlabel("aoa (deg)")
    PyPlot.ylabel("cd")

    #Verify cm
    cm_re2E5 = [0.0, 0.0401114, 0.0802228, 0.120334, 0.160446, 0.200557, 0.240667, 0.280786, 0.320866, 0.361109, 0.400666, 0.399322, 0.366641, 0.345388, 0.33302, 0.31785, 0.300225, 0.283491, 0.267061, 0.249874, 0.236365, 0.231921, 0.233325, 0.234063, 0.233147, 0.232362, 0.232035, 0.231672, 0.231194, 0.23083, 0.231063, 0.232505, 0.234934, 0.237435, 0.239672, 0.241888, 0.244174, 0.246447, 0.248704, 0.251045, 0.253574, 0.256388, 0.259466, 0.262604, 0.26568, 0.268737, 0.271811, 0.274892, 0.277971, 0.281039, 0.284091, 0.287124, 0.290141, 0.293154, 0.296173, 0.299185, 0.302185, 0.305216, 0.308288, 0.311258, 0.313977, 0.316475, 0.318839, 0.321154, 0.323486, 0.325822, 0.328126, 0.330469, 0.332921, 0.33524, 0.337062, 0.338422, 0.339614, 0.340781, 0.341959, 0.343136, 0.344266, 0.345422, 0.346749, 0.348006, 0.348586, 0.348416, 0.348042, 0.347752, 0.34749, 0.347204, 0.346862, 0.34652, 0.346363, 0.346255, 0.34548, 0.343777, 0.341845, 0.340126, 0.338459, 0.336732, 0.334961, 0.333226, 0.331633, 0.329924, 0.327588, 0.324674, 0.321618, 0.318577, 0.315536, 0.312486, 0.30942, 0.306394, 0.303416, 0.30028, 0.2968, 0.293025, 0.289076, 0.285081, 0.281117, 0.277161, 0.273195, 0.269242, 0.265288, 0.261273, 0.257137, 0.252835, 0.248414, 0.243991, 0.239597, 0.235202, 0.2308, 0.226398, 0.221998, 0.217599, 0.213204, 0.20882, 0.20444, 0.200059, 0.195675, 0.191289, 0.186914, 0.182544, 0.178129, 0.173709, 0.169471, 0.165447, 0.161455, 0.157411, 0.153361, 0.149325, 0.145299, 0.141259, 0.137189, 0.133162, 0.129283, 0.125495, 0.121695, 0.11797, 0.114353, 0.110405, 0.106109, 0.10324, 0.101272, 0.0937018, 0.0763387, 0.051795, 0.0229491, -0.00917348, -0.0400863, -0.0608273, -0.0680262, -0.0651137, -0.055863, -0.0440503, -0.0335397, -0.0297563, -0.0221719, -0.0146715, -0.0171736, -0.0188859, -0.0192802, -0.0210238, -0.0232045, -0.0246555, -0.030328, -0.0335065, -0.0360356, -0.0389232, -0.0392011, -0.0406963, -0.0415387, -0.0416685, -0.0319461, -0.00858319, 0.000901334, -0.0145788, -0.0325834, -0.0440667, -0.0546156, -0.0669185, -0.0794428, -0.0914897, -0.103857, -0.115837, -0.124449, -0.129001, -0.132465, -0.136813, -0.141524, -0.145995, -0.150362, -0.154791, -0.159246, -0.163644, -0.167896, -0.17197, -0.175968, -0.179995, -0.184044, -0.188085, -0.192123, -0.196167, -0.200203, -0.204205, -0.208147, -0.212014, -0.215839, -0.219671, -0.223513, -0.22735, -0.231189, -0.235046, -0.238886, -0.242645, -0.246297, -0.249845, -0.253327, -0.256811, -0.26031, -0.263799, -0.267288, -0.270816, -0.274334, -0.277704, -0.280882, -0.283924, -0.286896, -0.289861, -0.29284, -0.295808, -0.298768, -0.301786, -0.304824, -0.307658, -0.310174, -0.312505, -0.314791, -0.317077, -0.31937, -0.321647, -0.323904, -0.326229, -0.328626, -0.330801, -0.3325, -0.333914, -0.335328, -0.336784, -0.338241, -0.339669, -0.341057, -0.342514, -0.344131, -0.34553, -0.346197, -0.346347, -0.346536, -0.346867, -0.3472, -0.347488, -0.347753, -0.348042, -0.348416, -0.348586, -0.348006, -0.346749, -0.345422, -0.344266, -0.343136, -0.341959, -0.340781, -0.339614, -0.338422, -0.337062, -0.33524, -0.332921, -0.330469, -0.328126, -0.325822, -0.323486, -0.321154, -0.318839, -0.316475, -0.313977, -0.311258, -0.308288, -0.305216, -0.302185, -0.299185, -0.296173, -0.293154, -0.290141, -0.287124, -0.284091, -0.281039, -0.277971, -0.274892, -0.271811, -0.268737, -0.26568, -0.262604, -0.259466, -0.256388, -0.253574, -0.251045, -0.248704, -0.246447, -0.244174, -0.241888, -0.239672, -0.237435, -0.234935, -0.232507, -0.23106, -0.230823, -0.231204, -0.2317, -0.232001, -0.232256, -0.233254, -0.234462, -0.232991, -0.230442, -0.237363, -0.258978, -0.286613, -0.313062, -0.339806, -0.36782, -0.392544, -0.413458, -0.446391, -0.4948, -0.500786, -0.451398, -0.40108, -0.350983, -0.300834, -0.250697, -0.200557, -0.150418, -0.100279, -0.0501393, 0.0]

    cm_re5E5 = [-2.60209e-18, 0.0401114, 0.0802228, 0.120334, 0.160446, 0.200557, 0.240668, 0.280782, 0.320878, 0.360973, 0.399207, 0.420928, 0.441516, 0.461973, 0.482448, 0.502915, 0.523389, 0.543842, 0.564381, 0.584585, 0.593738, 0.587841, 0.571248, 0.552658, 0.536268, 0.519885, 0.503103, 0.486823, 0.469807, 0.451396, 0.437477, 0.431922, 0.428842, 0.424259, 0.419466, 0.415072, 0.410716, 0.406197, 0.401537, 0.397303, 0.394142, 0.392302, 0.391575, 0.391045, 0.390273, 0.389482, 0.388761, 0.388075, 0.387291, 0.386395, 0.385936, 0.386235, 0.386735, 0.387085, 0.38742, 0.387793, 0.388161, 0.388525, 0.388897, 0.389251, 0.389556, 0.389804, 0.390007, 0.3902, 0.390405, 0.39062, 0.390819, 0.390985, 0.39123, 0.391567, 0.391509, 0.3908, 0.389939, 0.389216, 0.38851, 0.387763, 0.387001, 0.386282, 0.385615, 0.384768, 0.383499, 0.381765, 0.379683, 0.377553, 0.3755, 0.373445, 0.371374, 0.369309, 0.367239, 0.36516, 0.363122, 0.36115, 0.359191, 0.357217, 0.355247, 0.353275, 0.351284, 0.349331, 0.347428, 0.345337, 0.342826, 0.339878, 0.336619, 0.333321, 0.330094, 0.326873, 0.323623, 0.32035, 0.317155, 0.314012, 0.310497, 0.306446, 0.302298, 0.298273, 0.294253, 0.290202, 0.28615, 0.282114, 0.278089, 0.274012, 0.26983, 0.265543, 0.261187, 0.256823, 0.252475, 0.248126, 0.243774, 0.23942, 0.235072, 0.230725, 0.226358, 0.221965, 0.217567, 0.213176, 0.208787, 0.204393, 0.199993, 0.195611, 0.191243, 0.186798, 0.1822, 0.177459, 0.172627, 0.167779, 0.162961, 0.158174, 0.153317, 0.148363, 0.143762, 0.13925, 0.132826, 0.123976, 0.114784, 0.106157, 0.0975518, 0.0887678, 0.0799063, 0.0713554, 0.0629929, 0.0534064, 0.0408751, 0.0252501, 0.00855849, -0.00782357, -0.0247386, -0.039237, -0.0455294, -0.0438599, -0.0375492, -0.0294267, -0.0245347, -0.0226448, -0.0208246, -0.0208057, -0.0258966, -0.0304884, -0.0322203, -0.0354409, -0.0380452, -0.0406207, -0.0433354, -0.0459706, -0.0482285, -0.0502751, -0.0522104, -0.0529389, -0.0541183, -0.0466567, -0.0339674, -0.0286748, -0.0260163, -0.0231364, -0.0205969, -0.0184467, -0.0177887, -0.0179169, -0.018951, -0.0220533, -0.0267795, -0.0390694, -0.0577945, -0.0723581, -0.0825314, -0.0934588, -0.105365, -0.116977, -0.128212, -0.139537, -0.151786, -0.163467, -0.171232, -0.175913, -0.180738, -0.186132, -0.191372, -0.196493, -0.20166, -0.206906, -0.212115, -0.21705, -0.221553, -0.225594, -0.2294, -0.233251, -0.237156, -0.241045, -0.244908, -0.248778, -0.252718, -0.256608, -0.2602, -0.263571, -0.266962, -0.270397, -0.273819, -0.27723, -0.280643, -0.284077, -0.287501, -0.290838, -0.294042, -0.297114, -0.300121, -0.303142, -0.30618, -0.309216, -0.312228, -0.315245, -0.318348, -0.321383, -0.324042, -0.32644, -0.328875, -0.331364, -0.333837, -0.336291, -0.338749, -0.341252, -0.343729, -0.346016, -0.348031, -0.349787, -0.351427, -0.353099, -0.354792, -0.356471, -0.358162, -0.359853, -0.361483, -0.363166, -0.365115, -0.367236, -0.369321, -0.371374, -0.373443, -0.375501, -0.377554, -0.379683, -0.381765, -0.383499, -0.384768, -0.385615, -0.386282, -0.387001, -0.387763, -0.38851, -0.389216, -0.389939, -0.3908, -0.391509, -0.391567, -0.39123, -0.390985, -0.390819, -0.39062, -0.390405, -0.3902, -0.390007, -0.389804, -0.389556, -0.389251, -0.388897, -0.388525, -0.388161, -0.387793, -0.38742, -0.387085, -0.386735, -0.386235, -0.385936, -0.386395, -0.387291, -0.388075, -0.388761, -0.389482, -0.390273, -0.391045, -0.391575, -0.392302, -0.394142, -0.397303, -0.401537, -0.406197, -0.410716, -0.415072, -0.419466, -0.424259, -0.428842, -0.431923, -0.437478, -0.451393, -0.469807, -0.486834, -0.503094, -0.519842, -0.536335, -0.5528, -0.570867, -0.587485, -0.596293, -0.594112, -0.583864, -0.573372, -0.562942, -0.552497, -0.542059, -0.531607, -0.52119, -0.510509, -0.497538, -0.451192, -0.401086, -0.350979, -0.300835, -0.250697, -0.200557, -0.150418, -0.100279, -0.0501393, 1.73472e-18]

    cm_re1E6 = [0.0, 0.0401114, 0.0802228, 0.120334, 0.160446, 0.200557, 0.240668, 0.280781, 0.320881, 0.360971, 0.399788, 0.427753, 0.454915, 0.481987, 0.509069, 0.536149, 0.563223, 0.590321, 0.61732, 0.644685, 0.670585, 0.657609, 0.625808, 0.605623, 0.587618, 0.566618, 0.545771, 0.525384, 0.504346, 0.484576, 0.468724, 0.458265, 0.452444, 0.447224, 0.441001, 0.43478, 0.428844, 0.422983, 0.416696, 0.41019, 0.405626, 0.403671, 0.402364, 0.40085, 0.399355, 0.397908, 0.396453, 0.394975, 0.393493, 0.392099, 0.391039, 0.390686, 0.390844, 0.390946, 0.390943, 0.390973, 0.391021, 0.391058, 0.39109, 0.391126, 0.391183, 0.391258, 0.39134, 0.391423, 0.391505, 0.39159, 0.391674, 0.391748, 0.391832, 0.391938, 0.391818, 0.391198, 0.390353, 0.389604, 0.388895, 0.388156, 0.38739, 0.386657, 0.386006, 0.385191, 0.383901, 0.382162, 0.380139, 0.378057, 0.376033, 0.374014, 0.371979, 0.369939, 0.367917, 0.365911, 0.363826, 0.361618, 0.359386, 0.357182, 0.354983, 0.352771, 0.350549, 0.34836, 0.346197, 0.343875, 0.341199, 0.338065, 0.334601, 0.331136, 0.327747, 0.324349, 0.320921, 0.317489, 0.314135, 0.310764, 0.307043, 0.302968, 0.298841, 0.294744, 0.290639, 0.28653, 0.282422, 0.278318, 0.274212, 0.270091, 0.265911, 0.261631, 0.257297, 0.252977, 0.248668, 0.244353, 0.240037, 0.235723, 0.231407, 0.227093, 0.222785, 0.218483, 0.214182, 0.209881, 0.20558, 0.201281, 0.19698, 0.192671, 0.188378, 0.184083, 0.179573, 0.17472, 0.169763, 0.164881, 0.160034, 0.155142, 0.150161, 0.145371, 0.140813, 0.135334, 0.127956, 0.11894, 0.108996, 0.0989541, 0.0891387, 0.0793745, 0.0694069, 0.059248, 0.0500114, 0.0407189, 0.028469, 0.0146397, 0.00104805, -0.0124024, -0.0266581, -0.0369221, -0.0361635, -0.0290357, -0.0235478, -0.0217691, -0.0200313, -0.0181787, -0.0177643, -0.0195665, -0.0260045, -0.0312463, -0.0334143, -0.03636, -0.0391082, -0.04183, -0.0445148, -0.0471749, -0.0497118, -0.0519225, -0.0534362, -0.055705, -0.0530825, -0.0459715, -0.0396464, -0.0342765, -0.0315643, -0.0299929, -0.0281516, -0.0259659, -0.0240525, -0.0230925, -0.0231209, -0.0225389, -0.0288657, -0.0119173, 0.00063643, -0.0227834, -0.0448096, -0.0586258, -0.0742675, -0.0913822, -0.107784, -0.124392, -0.141024, -0.155611, -0.16649, -0.17314, -0.177751, -0.182877, -0.188426, -0.193778, -0.199022, -0.20435, -0.209877, -0.215126, -0.219467, -0.223233, -0.226996, -0.230801, -0.234584, -0.238361, -0.242145, -0.245934, -0.249714, -0.253445, -0.257019, -0.260396, -0.263716, -0.267076, -0.270445, -0.273801, -0.27715, -0.280516, -0.283901, -0.287212, -0.290386, -0.293481, -0.296554, -0.299626, -0.302702, -0.30578, -0.308853, -0.311922, -0.31501, -0.318064, -0.320892, -0.323496, -0.326089, -0.328731, -0.331371, -0.333992, -0.336609, -0.33927, -0.341931, -0.344396, -0.34659, -0.348575, -0.350466, -0.352372, -0.354297, -0.356214, -0.358132, -0.360052, -0.361954, -0.363871, -0.365867, -0.367914, -0.36995, -0.371978, -0.374012, -0.376034, -0.378058, -0.380138, -0.382161, -0.383901, -0.385191, -0.386006, -0.386657, -0.38739, -0.388156, -0.388895, -0.389604, -0.390353, -0.391198, -0.391818, -0.391938, -0.391832, -0.391748, -0.391674, -0.39159, -0.391505, -0.391423, -0.39134, -0.391258, -0.391183, -0.391126, -0.39109, -0.391058, -0.391021, -0.390973, -0.390943, -0.390946, -0.390844, -0.390686, -0.391039, -0.392099, -0.393493, -0.394975, -0.396453, -0.397908, -0.399355, -0.40085, -0.402364, -0.403671, -0.405626, -0.41019, -0.416696, -0.422983, -0.428844, -0.43478, -0.441001, -0.447225, -0.452444, -0.458264, -0.468726, -0.484578, -0.504335, -0.525389, -0.545817, -0.566559, -0.587457, -0.605996, -0.626247, -0.655667, -0.67023, -0.654109, -0.636832, -0.619844, -0.602778, -0.58573, -0.56868, -0.551621, -0.534588, -0.517325, -0.498027, -0.451172, -0.401088, -0.350978, -0.300835, -0.250697, -0.200557, -0.150418, -0.100279, -0.0501393, 5.20417e-18]

    extrap_cm_old = (cm_re2E5, cm_re5E5, cm_re1E6)
    ND_corr3Dextr_maxerror_cm = 0.0

    #Plot the output and get the error
    extrap_cm = Array{Float64}(undef, length(extrap_aoas), length(Re_airfoiltools))
    PyPlot.figure("Verify_cm")
    for i = 1:length(Re_airfoiltools) #length of the airfoiltools data
        for j = 1:length(extrap_aoas)
            vars = [extrap_aoas[j], Re_airfoiltools[i], 0.0]
            extrap_cm[j,i] = AirfoilPrep.interpND(splout_extrap[3], vars)
        end

        # println("cm")
        # println(extrap_cm[:,i])
        PyPlot.plot(extrap_aoas, extrap_cm[:,i], label = "Re: $(Re_airfoiltools[i])")
        ND_corr3Dextr_maxerror_cm = max(ND_corr3Dextr_maxerror_cm, maximum(abs.(extrap_cm[:,i] - extrap_cm_old[i])))
    end
    PyPlot.xlabel("aoa (deg)")
    PyPlot.ylabel("cm")

    return ND_corr3Dextr_maxerror_cl, ND_corr3Dextr_maxerror_cd, ND_corr3Dextr_maxerror_cm
end

#------- RUN TESTS ------#
ERROR_TOL = 1E-1
NDSpl_maxerror_cl, NDSpl_maxerror_cd, NDSpl_maxerror_cm, NDtable = validateNDtools_from_Xfoil()
ND_corr3Dextr_maxerror_cl, ND_corr3Dextr_maxerror_cd, ND_corr3Dextr_maxerror_cm = verifyNDtable_extrap(NDtable)
correction3D_error_max = verify_correction3D_2()

@test correction3D_error_max <= ERROR_TOL
@test NDSpl_maxerror_cl <= ERROR_TOL
@test NDSpl_maxerror_cd <= ERROR_TOL
@test NDSpl_maxerror_cm <= ERROR_TOL
@test ND_corr3Dextr_maxerror_cl <= ERROR_TOL
@test ND_corr3Dextr_maxerror_cd <= ERROR_TOL
@test ND_corr3Dextr_maxerror_cm <= ERROR_TOL
println("Tests Passed")


text:
module FHIRClientJSON

import Base64
import HTTP
import JSON3

include("types.jl")

include("credentials.jl")
include("headers.jl")
include("requests.jl")

end # end module FHIRClientJSON

ext: jl


text:
using Mill, Flux, Duff

"""
	The explainer works by estimating a Shappley values on the levels of instances.
	The goal is to identify a "skeleton" of the tree that is responsible for a sample
	having a high score. This, the core of the explainer is explanation of BagNode.

	In Shappley values, we need to have an ability to turn on and off instances. Since the
	rest of the sample is not interesting, we want to process only a subset.

	We use `replacein` ability to create only and upper part of the model and the datamodel,
	allows efficient sampling. Also, we replace BagNode with a WeightedBagNode, such that
	we can do masking (WeightedBagNode would need a little more work).

	Once we have prepared upper model with upper data, we can run and estimator of Shappley
	values. We keep only those with positive values

	The main function is `explain(root_data, root_model, ϕ,  explaining_fun, recurse_depth)`
	where

	`root_data` is a root of a hierarchical explained sample
	`root_model` is a corresponding root of a hierarchical model
	`ϕ(mask, bagnode, x, mm)` is a function evaluating an effect when instances `x` set to `falses` in `mask`
			are removed from a `bagnode` and evaluated by a model `mm`. See function `onoff` for an example
			of such function.

	`explaining_fun(d, f)` estimate importance of `d` components of function `f`. Note that explaining function
	does not have any clue what `d` and `f` are. They just measure importance

"""
function explain(ds, model, i)
	sub_ds = deepcopy(ds)
	for index in setdiff(list(model),"")
		try
			sub_ds = explain(sub_ds[index], model[index], sub_ds, model, (x...) -> onoff(x..., i), explaining_fun)
		catch
		end
	end
	sub_ds
end

function explain(a::BagNode, m::BagModel, root_data, root_model, ϕ,  explaining_fun)
	# println("explaining BagNode")
	nobs(a) == 0 && return(root_data)
	xim = m.im(a.data)
	bagnode = WeightedBagNode(xim, a.bags, fill(true,size(xim.data, 2)))
	x = replacein(root_data, a, bagnode)
	mm = replacein(root_model, m.im, ArrayModel(identity))

	mask = explaining_fun(length(bagnode.weights), mask -> ϕ(mask, bagnode, x, mm))
	sa = removeinstances(a, mask)
	# @show sa
	replacein(root_data, a, sa)
end

function explain(a::TreeNode, m::ProductModel, root_data, root_model, ϕ,  explaining_fun)
	root_data
end

function explain(a::ArrayNode{T,A}, m::ArrayModel, root_data, root_model, ϕ,  explaining_fun) where {T<:SparseMatrixCSC, A}
	x = a.data
	nnz(x) == 0 && return(root_data)
	original_x = deepcopy(x.nzval)

	mask = explaining_fun(nnz(x),
		mask -> ϕ(mask, x.nzval, original_x, root_data, root_model)
		)
	x.nzval[ .! mask] .= 0
	root_data
end

function explain(a::ArrayNode{T,A}, m::ArrayModel, root_data, root_model, ϕ,  explaining_fun) where {T<:Matrix, A}
	x = a.data
	all(x .== 0) && return(root_data)
	original_x = deepcopy(x)

	mask = explaining_fun(nnz(x),
		mask -> ϕ(mask, x.nzval, original_x, root_data, root_model)
		)
	x.nzval[ .! mask] .= 0
	root_data
end

function onoff(mask, bagnode::WeightedBagNode, ds, model, i)
	fill!(bagnode.weights, false);
	bagnode.weights[mask] .= true;
	sum(m(x).data[i,:])
end

explaining_fun(d, f) = explaining_fun(d, f, 0.5, 100)

"""
		function explaining_fun(d, f, p, nprobes)

		d
"""
function explaining_fun(d, f, p, nprobes)
	daf = Duff.onlinedaf(d, f, p, nprobes)
	scores = Duff.meanscore(daf)
	I = sortperm(scores)
	# @show scores
	mask = trues(d)
	if f(mask) <= 0.5
		# @info "nothing to explain"
		return(mask)
	end
	# @info "before removing feature / intances $(f(mask))"
	for i in 1:d
		mask[I[i]] = false
		if f(mask) <=0.5
			mask[I[i]] = true
			break
		end
	end
	# @show mask
	# @info "after removing feature / intances $(f(mask))"
	isnan(f(mask)) && @error "f(mask) is nan"
	mask
end


text:
# Autogenerated wrapper script for libigc_jll for i686-linux-gnu-cxx03
export GenX_IR, iga, libiga, libigc, libigdfcl, libopencl_clang

JLLWrappers.@generate_wrapper_header("libigc")
JLLWrappers.@declare_executable_product(GenX_IR)
JLLWrappers.@declare_executable_product(iga)
JLLWrappers.@declare_library_product(libiga, "libiga32.so.1")
JLLWrappers.@declare_library_product(libigc, "libigc.so.1")
JLLWrappers.@declare_library_product(libigdfcl, "libigdfcl.so.1")
JLLWrappers.@declare_library_product(libopencl_clang, "libopencl-clang.so.10")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        GenX_IR,
        "bin/GenX_IR",
    )

    JLLWrappers.@init_executable_product(
        iga,
        "bin/iga32",
    )

    JLLWrappers.@init_library_product(
        libiga,
        "lib/libiga32.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libigc,
        "lib/libigc.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libigdfcl,
        "lib/libigdfcl.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencl_clang,
        "lib/libopencl-clang.so.10",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function state_at(state_var, uns, pos)
    if pos < 1 || pos > length(state_var)
        return uns
    end
    return state_var[pos]
end


text:
module NumericalMath

# package code goes here

  export trisolve
  export pval, horner, pzero
  export interp1d, pchip
  export ridders, brent_dekker
  export fd_gradient, fd_jacobian, fd_hessian, fd_laplacian,
         numderiv, complex_step
  export trapz, romberg
  export agm, arc_length
  export lambertW

  include("la.jl")
  include("polynom.jl")
  include("interpolate.jl")
  include("fzero.jl")
  include("derivative.jl")
  include("integrate.jl")
  include("misc.jl")
  include("specfun.jl")

end # module


text:
using Test
using JSON
using MianalyzerBackend
const GMeansController = MianalyzerBackend.GMeansController

include("../TestSupport.jl")

@testset "GMeanController" begin
  @testset ".call" begin
    request = Dict([("records", [[1.0, 1.0], [2.0, 2.0], [3.0, 3.0], [4.0, 4.0], [5.0, 5.0]])])

    actual = JSON.parse(IOBuffer(GMeansController.call(request)))

    @test actual["clusters"] == [[[1.0, 1.0], [2.0, 2.0], [3.0, 3.0], [4.0, 4.0], [5.0, 5.0]]]
    @test actual["centers"] == [[3.0, 3.0]]
  end
end


text:
using ApproxFunBase, ApproxFunOrthogonalPolynomials, LinearAlgebra, SpecialFunctions, BlockBandedMatrices, Test
import ApproxFunBase: testbandedblockbandedoperator, testraggedbelowoperator, factor, Block, cfstype,
                    blocklengths, block, tensorizer, Vec, ArraySpace, ∞,
                    testblockbandedoperator
import ApproxFunOrthogonalPolynomials: chebyshevtransform

@testset "Multivariate" begin
    @testset "Square" begin
        S = Space(ChebyshevInterval()^2)
        @test @inferred(blocklengths(S)) ≡ Base.OneTo(∞)

        @test block(tensorizer(S), 1) == Block(1)

        @time for k=0:5,j=0:5
            ff=(x,y)->cos(k*acos(x))*cos(j*acos(y))
            f=Fun(ff,ChebyshevInterval()^2)
            @test f(0.1,0.2) ≈ ff(0.1,0.2)
        end

        @time for k=0:5,j=0:5
            ff = (x,y)->cos(k*acos(x/2))*cos(j*acos(y/2))
            f=Fun(ff,Interval(-2,2)^2)
            @test f(0.1,0.2) ≈ ff(0.1,0.2)
        end

        @time for k=0:5,j=0:5
            ff=(x,y)->cos(k*acos(x-1))*cos(j*acos(y-1))
            f=Fun(ff,Interval(0,2)^2)
            @test f(0.1,0.2) ≈ ff(0.1,0.2)
        end

        ## Try constructor variants

        ff=(x,y)->exp(-10(x+.2)^2-20(y-.1)^2)*cos(x*y)
        gg=x->exp(-10(x[1]+.2)^2-20(x[1]-.1)^2)*cos(x[1]*x[2])
        f=Fun(ff,ChebyshevInterval()^2,10000)
        @test f(0.,0.) ≈ ff(0.,0.)

        f=Fun(gg,ChebyshevInterval()^2,10000)
        @test f(0.,0.) ≈ ff(0.,0.)

        f=Fun(ff,ChebyshevInterval()^2)
        @test f(0.,0.) ≈ ff(0.,0.)
        f=Fun(gg,ChebyshevInterval()^2)
        @test f(0.,0.) ≈ ff(0.,0.)

        f=Fun(ff)
        @test f(0.,0.) ≈ ff(0.,0.)
        f=Fun(gg)
        @test f(0.,0.) ≈ ff(0.,0.)

        # Fun +-* constant
        f=Fun((x,y)->exp(x)*cos(y))

        @test f(0.1,0.2)+2 ≈ (f+2)(0.1,0.2)
        @test f(0.1,0.2)-2 ≈ (f-2)(0.1,0.2)
        @test f(0.1,0.2)*2 ≈ (f*2)(0.1,0.2)
    end

    @testset "LowRankFun" begin
        @time F = LowRankFun((x,y)->besselj0(10(y-x)),Chebyshev(),Chebyshev())

        @test F(.123,.456) ≈ besselj0(10(.456-.123))

        @time G = LowRankFun((x,y)->besselj0(10(y-x));method=:Cholesky)

        @test G(.357,.246) ≈ besselj0(10(.246-.357))

        # test "fast" grid evaluation of LowRankFun
        f = LowRankFun((x,y) -> exp(x) * cos(y)); n = 1000
        x = range(-1, stop=1, length=n); y = range(-1, stop=1, length=n)
        X = x * fill(1.0,1,n); Y = fill(1.0, n) * y'
        @time v1 = f.(X, Y);
        @time v2 = f.(x, y');
        @test v1 ≈ v2
    end

    @testset  "Vec segment" begin
        d = Segment(Vec(0.,0.) , Vec(1.,1.))
        x = Fun()
        @test (ApproxFunBase.complexlength(d)*x/2)(0.1)  ≈ (d.b - d.a)*0.1/2
        @test ApproxFunBase.fromcanonical(d,x)(0.1) ≈ (d.b+d.a)/2 + (d.b - d.a)*0.1/2

        x,y = Fun(Segment(Vec(0.,0.) , Vec(2.,1.)))
        @test x(0.2,0.1) ≈ 0.2
        @test y(0.2,0.1) ≈ 0.1

        d=Segment((0.,0.),(1.,1.))
        f=Fun(xy->exp(-xy[1]-2cos(xy[2])),d)
        @test f(0.5,0.5) ≈ exp(-0.5-2cos(0.5))
        @test f(Vec(0.5,0.5)) ≈ exp(-0.5-2cos(0.5))

        f=Fun(xy->exp(-xy[1]-2cos(xy[2])),d,20)
        @test f(0.5,0.5) ≈ exp(-0.5-2cos(0.5))

        f=Fun((x,y)->exp(-x-2cos(y)),d)
        @test f(0.5,0.5) ≈ exp(-0.5-2cos(0.5))

        f=Fun((x,y)->exp(-x-2cos(y)),d,20)
        @test f(0.5,0.5) ≈ exp(-0.5-2cos(0.5))
    end

    @testset "Multivariate calculus" begin
        ## Sum
        ff = (x,y) -> (x-y)^2*exp(-x^2/2-y^2/2)
        f=Fun(ff, (-4..4)^2)
        @test f(0.1,0.2) ≈ ff(0.1,0.2)

        @test sum(f,1)(0.1) ≈ 2.5162377980828357
        f=LowRankFun(f)
        @test evaluate(f.A,0.1) ≈ map(f->f(0.1),f.A)
    end

    @testset "KroneckerOperator" begin
        Mx = Multiplication(Fun(cos),Chebyshev())
        My = Multiplication(Fun(sin),Chebyshev())
        K = Mx⊗My

        @test BandedBlockBandedMatrix(view(K,1:10,1:10)) ≈ [K[k,j] for k=1:10,j=1:10]
        C = Conversion(Chebyshev()⊗Chebyshev(),Ultraspherical(1)⊗Ultraspherical(1))
        @test C[1:100,1:100] ≈ Float64[C[k,j] for k=1:100,j=1:100]
    end

    @time @testset "Partial derivative operators" begin
        d = Space(0..1) * Space(0..2)
        Dx = Derivative(d, [1,0])
        testbandedblockbandedoperator(Dx)
        f = Fun((x,y) -> sin(x) * cos(y), d)
        fx = Fun((x,y) -> cos(x) * cos(y), d)
        @test (Dx*f)(0.2,0.3) ≈ fx(0.2,0.3)
        Dy = Derivative(d, [0,1])
        testbandedblockbandedoperator(Dy)
        fy = Fun((x,y) -> -sin(x) * sin(y), d)
        @test (Dy*f)(0.2,0.3) ≈ fy(0.2,0.3)
        L = Dx + Dy
        testbandedblockbandedoperator(L)
        @test (L*f)(0.2,0.3) ≈ (fx(0.2,0.3)+fy(0.2,0.3))

        B=ldirichlet(factor(d,1))⊗ldirichlet(factor(d,2))
        @test abs(Number(B*f)-f(0.,0.)) ≤ 10eps()

        B=Evaluation(factor(d,1),0.1)⊗ldirichlet(factor(d,2))
        @test Number(B*f) ≈ f(0.1,0.)

        B=Evaluation(factor(d,1),0.1)⊗Evaluation(factor(d,2),0.3)
        @test Number(B*f) ≈ f(0.1,0.3)

        B=Evaluation(d,(0.1,0.3))
        @test Number(B*f) ≈ f(0.1,0.3)
    end

    @time @testset "x,y constructors" begin
        d=ChebyshevInterval()^2

        sp = ArraySpace(d,2)
        @test blocklengths(sp) == 2:2:∞
        @test block(ArraySpace(d,2),1) == Block(1)

        x,y=Fun(d)
        @test x(0.1,0.2) ≈ 0.1
        @test y(0.1,0.2) ≈ 0.2

        x,y=Fun(identity, d, 20)
        @test x(0.1,0.2) ≈ 0.1
        @test y(0.1,0.2) ≈ 0.2

        # Boundary

        x,y=Fun(identity, ∂(d), 20)
        @test x(0.1,1.0) ≈ 0.1
        @test y(1.0,0.2) ≈ 0.2

        x,y=Fun(identity, ∂(d))
        @test x(0.1,1.0) ≈ 0.1
        @test y(1.0,0.2) ≈ 0.2

        x,y=Fun(∂(d))
        @test x(0.1,1.0) ≈ 0.1
        @test y(1.0,0.2) ≈ 0.2
    end

    @testset "conversion between" begin
        dx = dy = ChebyshevInterval()
        d = dx × dy
        x,y=Fun(d)
        @test x(0.1,0.2) ≈ 0.1
        @test y(0.1,0.2) ≈ 0.2

        x,y = Fun(∂(d))
        x,y = components(x),components(y)

        g = [real(exp(x[1]-1im));0.0y[2];real(exp(x[3]+1im));real(exp(-1+1im*y[4]))]
        B = [ Operator(I,dx)⊗ldirichlet(dy);
             ldirichlet(dx)⊗Operator(I,dy);
             Operator(I,dx)⊗rdirichlet(dy);
             rneumann(dx)⊗Operator(I,dy)    ]

        @test Fun(g[1],rangespace(B)[1])(-0.1,-1.0) ≈ g[1](-0.1,-1.0)
        @test Fun(g[3],rangespace(B)[3])(-0.1,1.0)  ≈ g[3](-0.1,1.0)

        A = [B; Laplacian()]

        @test cfstype([g;0.0]) == Float64
        g2 = Fun([g;0.0],rangespace(A))
        @test cfstype(g2) == Float64

        @test g2[1](-0.1,-1.0) ≈ g[1](-0.1,-1.0)
        @test g2[3](-0.1,1.0)  ≈ g[3](-0.1,1.0)
    end

    @testset "Cheby * Interval" begin
        d = ChebyshevInterval()^2
        x,y = Fun(∂(d))

        @test ApproxFunBase.rangetype(Space(∂(d))) == Float64
        @test ApproxFunBase.rangetype(space(y)) == Float64

        @test (im*y)(1.0,0.1) ≈ 0.1im
        @test (x+im*y)(1.0,0.1) ≈ 1+0.1im

        @test exp(x+im*y)(1.0,0.1) ≈ exp(1.0+0.1im)
    end

    @testset "DefiniteIntegral" begin
        f = Fun((x,y) -> exp(-x*cos(y)))
        @test Number(DefiniteIntegral()*f) ≈ sum(f)
    end

    @testset "Piecewise Tensor" begin
        a = Fun(0..1) + Fun(2..3)
        f = a ⊗ a
        @test f(0.1,0.2) ≈ 0.1*0.2
        @test f(1.1,0.2) ≈ 0
        @test f(2.1,0.2) ≈ 2.1*0.2

        @test component(space(f),1,1) == Chebyshev(0..1)^2
        @test component(space(f),1,2) == Chebyshev(0..1)*Chebyshev(2..3)
        @test component(space(f),2,1) == Chebyshev(2..3)*Chebyshev(0..1)
        @test component(space(f),2,2) == Chebyshev(2..3)^2
    end

    @testset "Bug in chop of ProductFun" begin
        u = Fun(Chebyshev()^2,[0.0,0.0])
        @test coefficients(chop(ProductFun(u),10eps())) == zeros(0,1)

        d= (-1..1)^2
        B=[Dirichlet(factor(d,1))⊗I;I⊗ldirichlet(factor(d,2));I⊗rneumann(factor(d,2))]
        Δ=Laplacian(d)

        rs = rangespace([B;Δ])
        f = Fun((x,y)->exp(-x^2-y^2),d)
        @test_throws DimensionMismatch coefficients([0.0;0.0;0.0;0.0;f],rs)
    end

    @testset "off domain evaluate" begin
        g = Fun(1, Segment(Vec(0,-1) , Vec(π,-1)))
        @test g(0.1,-1) ≈ 1
        @test g(0.1,1) ≈ 0
    end

    @testset "Dirichlet" begin
        testblockbandedoperator(Dirichlet((0..1)^2))
        testblockbandedoperator(Dirichlet((0..1) × (0.0 .. 1)))
        testraggedbelowoperator(Dirichlet(Chebyshev()^2))
        testraggedbelowoperator(Dirichlet(Chebyshev(0..1) * Chebyshev(0.0..1)))
    end

    @testset "2d derivative (issue #346)" begin
        d = Chebyshev()^2
        f = Fun((x,y) -> sin(x) * cos(y), d)
        C=Conversion(Chebyshev()⊗Chebyshev(),Ultraspherical(1)⊗Ultraspherical(1))
        @test (C*f)(0.1,0.2) ≈ f(0.1,0.2)
        Dx = Derivative(d, [1,0])
        f = Fun((x,y) -> sin(x) * cos(y), d)
        fx = Fun((x,y) -> cos(x) * cos(y), d)
        @test (Dx*f)(0.2,0.3) ≈ fx(0.2,0.3)
        Dy = Derivative(d, [0,1])
        fy = Fun((x,y) -> -sin(x) * sin(y), d)
        @test (Dy*f)(0.2,0.3) ≈ fy(0.2,0.3)
        L=Dx+Dy
        testbandedblockbandedoperator(L)

        @test (L*f)(0.2,0.3) ≈ (fx(0.2,0.3)+fy(0.2,0.3))

        B=ldirichlet(factor(d,1))⊗ldirichlet(factor(d,2))
        @test Number(B*f) ≈ f(-1.,-1.)

        B=Evaluation(factor(d,1),0.1)⊗ldirichlet(factor(d,2))
        @test Number(B*f) ≈ f(0.1,-1.)

        B=Evaluation(factor(d,1),0.1)⊗Evaluation(factor(d,2),0.3)
        @test Number(B*f) ≈ f(0.1,0.3)

        B=Evaluation(d,(0.1,0.3))
        @test Number(B*f) ≈ f(0.1,0.3)
    end

    @testset "ProductFun" begin
        u0   = ProductFun((x,y)->cos(x)+sin(y) +exp(-50x.^2-40(y-0.1)^2)+.5exp(-30(x+0.5)^2-40(y+0.2)^2))

        @test values(u0)-values(u0|>LowRankFun)|>norm < 1000eps()
        @test chebyshevtransform(values(u0))-coefficients(u0)|>norm < 100eps()

        ##TODO: need to do adaptive to get better accuracy
        @test sin(u0)(.1,.2)-sin(u0(.1,.2))|>abs < 10e-4

        F = LowRankFun((x,y)->hankelh1(0,10abs(y-x)),Chebyshev(1.0..2.0),Chebyshev(Segment(1.0im,2.0im)))

        @test F(1.5,1.5im) ≈ hankelh1(0,10abs(1.5im-1.5))
    end

    @testset "Functional*Fun" begin
        d=ChebyshevInterval()
        B=ldirichlet(d)
        f=ProductFun((x,y)->cos(cos(x)*sin(y)),d^2)

        @test norm(B*f-Fun(y->cos(cos(-1)*sin(y)),d))<20000eps()
        @test norm(f*B-Fun(x->cos(cos(x)*sin(-1)),d))<20000eps()
    end

    @testset "matrix" begin
        f=Fun((x,y)->[exp(x*cos(y));cos(x*sin(y));2],ChebyshevInterval()^2)
        @test f(0.1,0.2) ≈ [exp(0.1*cos(0.2));cos(0.1*sin(0.2));2]

        f=Fun((x,y)->[exp(x*cos(y)) cos(x*sin(y)); 2 1],ChebyshevInterval()^2)
        @test f(0.1,0.2) ≈ [exp(0.1*cos(0.2)) cos(0.1*sin(0.2));2 1]
    end
end


text:
using Fermi.HartreeFock

import Base: show

export RCCSD

"""
    Fermi.CoupledCluster.RCCSDAlgorithm

Abstract type for RCCSD implementations.
"""
abstract type RCCSDAlgorithm end

"""
    Fermi.CoupledCluster.RCCSD

Wave function object for Restricted Coupled Cluster Singles and Doubles.

# High Level Interface
Run a RCCSD computation and return the RCCSD object:
```
julia> @energy rccsd
```
Equivalent to
```
julia> Fermi.CoupledCluster.RCCSD()
```
This function calls a constructor that runs a RCCSD computation based on the options found in `Fermi.Options`.

# Fields

| Name   |   Description     |
|--------|---------------------|
| `guessenergy` |   Energy recovered in the first iteration, normally RMP2 |
| `correlation` |   Computed RCCSD correlation energy |
| `energy`   |   Total wave function energy (Reference energy + Correlation energy)      |
| `e_conv`   | ΔE from the last iteration  |
| `t_conv`   |  Amplitudes RMS change from the last iteration|

# Relevant options

These options can be set with `@set <option> <value>`

| Option         | What it does                      | Type      | choices [default]     |
|----------------|-----------------------------------|-----------|-----------------------|
| `cc_alg`      | Picks RCCSD algorithm              | `Int`     | [1]                   |
| `cc_e_conv`   | Energy convergence criterion           | `Float64` | [10^-10]              |
| `cc_max_rms`    | Amplitudes RMS convergence criterion   | `Float64` | [10^-10]              |
| `cc_max_iter`   | Max number of CC iterations   | `Int` | [50]              |
| `cc_damp_ratio` | Fraction of old amplitudes to be kept   | `Float64` | 0.0--1.0 [0.0]              |
| `cc_diis` | Whether to use DIIS   | `Bool` | `false` [`true`]              |
| `diis_start` | Iteration number where DIIS starts   | `Int` | [3]              |
| `cc_diis_relax` | Interval between DIIS extrapolations   | `Int` | [3]              |
| `cc_ndiis` | Maximum number of stored vectors for DIIS   | `Int` | [3]              |
| `basis`       | What basis set to use             | `String`  | ["sto-3g"]            |
| `df`          | Whether to use density fitting    | `Bool`    | `true` [`false`]      |
| `rifit`       | What aux. basis set to use for RI | `String`  | ["auto"]              |
| `drop_occ`    | Number of occupied electrons to be dropped | `Int`  | [0]              |
| `drop_vir`    | Number of virtual electrons to be dropped | `Int`  | [0]              |
"""
struct RCCSD{T} <: AbstractCCWavefunction
    guessenergy::Float64
    correlation::T
    energy::Float64
    T1::AbstractArray{T,2}
    T2::AbstractArray{T,4}
    e_conv::T
    t_conv::T
end

"""
    Fermi.CoupledCluster.get_rccsd_alg

Returns a singleton type corresponding to a RCCSD implementation based on the options.
"""
function get_rccsd_alg(N::Int = Options.get("cc_alg"))
    try
        return get_rccsd_alg(Val(N))
    catch MethodError
        throw(FermiException("implementation number $N not available for RCCSD."))
    end
end

# For each implementation a singleton type must be create
struct RCCSDa <: RCCSDAlgorithm end
include("RCCSDa.jl")
# And a number is assigned to the implementation
get_rccsd_alg(x::Val{1}) = RCCSDa()

function RCCSD(x...)
    if !any(i-> i isa RCCSDAlgorithm, x)
        RCCSD(x..., get_rccsd_alg())
    else
        # Print the type of arguments given for a better feedback
        args = "("
        for a in x[1:end-1]
            args *= "$(typeof(a)), "
        end
        args = args[1:end-2]*")"
        throw(FermiException("invalid arguments for RCCSD method: $args"))
    end
end

## MISCELLANEOUS
# Pretty printing
function string_repr(X::RCCSD)
    out = ""
    out = out*" ⇒ Fermi Restricted CCSD Wave function\n"
    out = out*" ⋅ Correlation Energy:     $(X.correlation)\n"
    out = out*" ⋅ Total Energy:           $(X.energy)"
    return out
end

function show(io::IO, ::MIME"text/plain", X::RCCSD)
    print(io, string_repr(X))
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
using RowTables
using DataStructures
using DataFrames
using Test

@testset "empty RowTable" begin
    @test size(RowTable()) == (0, 0)
end

@testset "some" begin
    a = [OrderedDict(:b => 3, :a => 1), OrderedDict(:b => 3, :a => 4.0), OrderedDict(:b => 3, :a => 1)]
    b = [Dict(:a => 1, :b => 3), Dict(:a => 4.0, :b => 3), Dict(:a => 1, :b => 3)]
    @test RowTable(a) == RowTable(b, [:b, :a])
end

let
    ds = [
        Dict(:a=>2, :b=>1, :c=>7)
        Dict(:a=>1, :b=>4, :c=>8)
        Dict(:a=>9, :b=>6, :c=>4)
        Dict(:a=>10, :b=>8, :c=>7)
        Dict(:a=>10, :b=>4, :c=>3)]

    rt = RowTable(ds, [:a, :b, :c])

    @testset "construction" begin
        @test RowTable(ds) == rt
    end

    @testset "conversion" begin
        # Construct from array of dictionaries
        df = DataFrame(rt)
        @test rt == RowTable(df)
    end

    @testset "renaming" begin
        rt = RowTable(ds, [:a, :b, :c])
        rt1 = deepcopy(rt)
        @test rename!(rename!(rt1, [:a => :b, :b => :a]), [:a => :b, :b => :a]) == rt
        cyc = [:b => :a, :c => :b, :a => :c]
        @test rename!(rename!(rename!(rt1, cyc), cyc), cyc) == rt
    end

    @testset "indexing" begin
        @test rt == rt[:]
        @test rt == rt[:,:]
        @test rt[:,1] == rt[1]
        @test rt[:,[1]] == rt[[1]]
        @test rt[1,1] == 2
        @test rt[1,:] == [2, 1, 7]
        @test rt[:,:b] == [1, 4, 6, 8, 4]
        @test rt[:,:b] == rt[:,2]
        @test rt[:,:a] == rt[:,1]
        @test rt[1, 2:3] == [1, 7]
        # Construct from array of arrays.
        @test RowTable(rows(rt), [:a, :b, :c]) == rt

        rs = Any[[2, 7], [1, 8], [9, 4], [10, 7], [10, 3]]
        @test RowTable(rs, [:a, :c]) == rt[:,[:a, :c]]
        @test RowTable([[4, 8], [6, 4], [8, 7]], [:b, :c]) == rt[2:4, [:b,:c]]
    end

    @testset "named tuples" begin
        @test RowTable(cols=[[1,2], [3,4]], names = [:a, :b], tuples=true)[1,:] == (a = 1, b = 3)
        @test RowTable(cols=[[1,2], [3,4]], names = [:a, :b], tuples=false)[1,:] == [1, 3]
        @test RowTable(ds, [:a, :b, :c], tuples=true).rows[1] == (a = 2, b = 1, c = 7)
    end

    @testset "iteration" begin
        s = 0
        for r in rt
            s += sum(r)
        end
        @test s == 84
    end

    @testset "pop! popfirst!" begin
        @test pop!(copy(rt)) == [10, 4, 3]
        @test popfirst!(copy(rt)) == [2, 1, 7]
    end

    @testset "sort" begin
        df = DataFrame(rt)
        df[:d] = ["b", "v", "d", "x", "x"]
        rt = RowTable(df)
        @test sort(rt) == RowTable(sort(df))
        @test sort(rt, [2]) == RowTable(sort(df, [2]))
        @test sort(rt, [2], rev=true) == RowTable(sort(df, [2], rev=true))
        @test sort(rt, [2], rev=true) != RowTable(sort(df, [2]))
        @test sort(rt, [:d, :c], rev=true) == RowTable(sort(df, [:d, :c], rev=true))
    end
end


text:
"""
    romberg(f, a, b; maxi = 10, delta = 10^-9)
berisi program untuk mencari integral numerik dari fungsi `f(x)` pada
interval `[a,b]`. Program ini secara default berisi 3 masukan, yaitu fungsi `f(x)`, dan titik
ujung interval `[a,b]`. Secara default memiliki nilai maksimum iterasi `maxi=10` dan
toleransi error sebesar `delta = 1e-9`
# Examples
```jl
julia> f(x) = (x^2+x+1)*cos(x)
f (generic function with 1 method)

julia> sol, flag, err, R = romberg(f, 0, pi/2);

julia> sol
2.0381974270672245

julia> flag
0

julia> err
1.213065203842234e-10

julia> R
6×6 Array{Float64,2}:
 0.785398  0.0      0.0     0.0     0.0     0.0
 1.72681   2.04062  0.0     0.0     0.0     0.0
 1.96053   2.03844  2.0383  0.0     0.0     0.0
 2.01879   2.03821  2.0382  2.0382  0.0     0.0
 2.03335   2.0382   2.0382  2.0382  2.0382  0.0
 2.03698   2.0382   2.0382  2.0382  2.0382  2.0382

 julia> sol, flag, err, R = romberg(f, 0, pi/2, maxi=20, delta=1e-14);

 julia> sol
 2.038197427067236

 julia> flag
 0

 julia> err
 0.0

 julia> R
 8×8 Array{Float64,2}:
  0.785398  0.0      0.0     0.0     0.0     0.0     0.0     0.0
  1.72681   2.04062  0.0     0.0     0.0     0.0     0.0     0.0
  1.96053   2.03844  2.0383  0.0     0.0     0.0     0.0     0.0
  2.01879   2.03821  2.0382  2.0382  0.0     0.0     0.0     0.0
  2.03335   2.0382   2.0382  2.0382  2.0382  0.0     0.0     0.0
  2.03698   2.0382   2.0382  2.0382  2.0382  2.0382  0.0     0.0
  2.03789   2.0382   2.0382  2.0382  2.0382  2.0382  2.0382  0.0
  2.03812   2.0382   2.0382  2.0382  2.0382  2.0382  2.0382  2.0382
```
return nilai solusi integral numerik `sol`, status solusi
`flag`, hampiran galat `err`, dan matriks `R` yaitu matriks iterasi Romberg.
"""
function romberg(f, a, b; maxi = 10, delta = 10^-9)
  flag = 1;
  M = 1;
  h = b-a;
  err = 1;
  R = h/2*(f(a)+f(b));
  for J = 1:maxi
    # Rekursif Trapesium
    h = h/2;
    M = 2*M;
    s = 0;
    for p = 1:M/2
      x = a+h*(2*p-1);
      s = s+f(x);
    end
    R = [R zeros(size(R,1),1) ; R[J,1]/2 + h*s zeros(1,size(R,1))]
    # Aturan Romberg
    for K=1:J
      R[J+1,K+1]=R[J+1,K]+(R[J+1,K]-R[J,K])/(4^K-1);
    end
    err = abs.(R[J,J]-R[J+1,J+1]);
    if err<delta; flag=0; break; end
  end
  sol = R[end,end];
  return sol, flag, err, R
end


text:
#=
testODE:
- Julia version:
- Author: ymocquar
- Date: 2019-11-13
=#
include("../src/twoscales_pure_ab.jl")
using DifferentialEquations
using LinearAlgebra
using Plots
using Random
function twoscales_solve(par_u0::PrepareU0, order, t, nb) end
function fctMain(n_tau)
    seed = 123321
    Random.seed!(seed)
    u0 = rand(BigFloat, 4)
    B = 2rand(BigFloat, 4, 4) - ones(BigFloat, 4, 4)
    println("seed = $seed B=$B")
    nbeps = 48
    nbmaxtest = 3
    paseps = big"10.0"^0.125
    ordmin = 3
    ordmax = 14
    ordpas = 1
    t_max = big"1.0"
    A = [0 0 1 0; 0 0 0 0; -1 0 0 0; 0 0 0 0]
    for order = ordmin:ordpas:ordmax
        y = ones(Float64, nbeps, nbmaxtest)
        x = zeros(Float64, nbeps)
        for i_eps = 1:nbeps
            epsilon = big"1.0" / paseps^i_eps
            x[i_eps] = epsilon
            fct = u -> B * u
            parphi = PreparePhi(epsilon, n_tau, A, fct, B)
            println("prepareU0 eps=$epsilon n_tau=$n_tau")
            @time par_u0 = PrepareU0(parphi, order + 2, u0)
            solref = getexactsol(parphi, u0, t_max)
            eps_v = convert(Float32, epsilon)
            println("epsilon = $eps_v solref=$solref")
            for indc = 1:nbmaxtest
                nb = 10 * 10^indc
                @time parGen = PrepareTwoScalePureAB(nb, t_max, order, par_u0)
                @time sol = twoscales_pure_ab(parGen, only_end = true)
                println("solref=$solref")
                println("nb=$nb sol=$sol")
                diff = solref - sol
                println("nb=$nb dt=$(1.0/nb) normInf=$(norm(diff,Inf)) norm2=$(norm(diff))")
                y[i_eps, indc] = norm(diff, Inf)
                println("epsilon=$epsilon result=$y")
            end
        end
        labels = Array{String,2}(undef, 1, nbmaxtest)
        for i = 1:nbmaxtest
            nb = 10 * 10^i
            labels[1, i] = " delta t,order=$(1/nb), $order "
            mindiff = minimum(y[:, i])
            bornediff = min(mindiff * 1e20, 1)
            for j = 1:nbeps
                if y[j, i] >= bornediff
                    y[j, i] = NaN
                end
            end
        end
        gr()
        p = Plots.plot(
            x,
            y,
            xlabel = "epsilon",
            xaxis = :log,
            ylabel = "error",
            yaxis = :log,
            legend = :topleft,
            label = labels,
            marker = 2,
        )
        prec_v = precision(BigFloat)
        Plots.savefig(p, "out/r_$(prec_v)_$(order)_$(n_tau)_epsilon_lgn_2.pdf")
    end
end
fctMain(32)


text:
module TestLP

import MathOptInterface
using Test

const MOI = MathOptInterface
const MOIU = MOI.Utilities
const LP = MOI.FileFormats.LP
const LP_TEST_FILE = "test.lp"

function test_show()
    @test sprint(show, LP.Model()) == "A .LP-file model"
end

function test_comprehensive_write()
    model = LP.Model()
    MOIU.loadfromstring!(
        model,
        """
variables: a, x, y, z
minobjective: x
x >= -1.0
x <= 2.0
y == 3.0
z in Interval(4.0, 5.0)
c5: 1.1x + 0.0 <= 5.1
c6: 1.3x + -1.4 >= -0.1
c7: 1.5a + 1.6 == 0.2
c8: 1.7a + 1.8 in Interval(0.3, 0.4)
x in ZeroOne()
y in Integer()
""",
    )
    MOI.write_to_file(model, LP_TEST_FILE)
    @test read(LP_TEST_FILE, String) ==
          "minimize\n" *
          "obj: x\n" *
          "subject to\n" *
          "c5: 1.1 x <= 5.1\n" *
          "c6: -1.4 + 1.3 x >= -0.1\n" *
          "c7: 1.6 + 1.5 a = 0.2\n" *
          "c8: 0.3 <= 1.8 + 1.7 a <= 0.4\n" *
          "Bounds\n" *
          "x <= 2\n" *
          "x >= -1\n" *
          "y = 3\n" *
          "4 <= z <= 5\n" *
          "a free\n" *
          "General\n" *
          "y\n" *
          "Binary\n" *
          "x\n" *
          "End\n"

    @test !MOI.is_empty(model)
    MOI.empty!(model)
    @test MOI.is_empty(model)
end

function test_name_sanitization_start()
    for starting_letter in [".", "0", "E", "e"]
        model = LP.Model()
        MOI.Utilities.loadfromstring!(
            model,
            """
variables: x
minobjective: x
c1: 1.0 * x >= -1.0
""",
        )
        x = MOI.get(model, MOI.VariableIndex, "x")
        MOI.set(model, MOI.VariableName(), x, starting_letter * "x")
        c = MOI.get(model, MOI.ConstraintIndex, "c1")
        MOI.set(model, MOI.ConstraintName(), c, starting_letter * "c1")
        MOI.write_to_file(model, LP_TEST_FILE)
        @test read(LP_TEST_FILE, String) ==
              "minimize\n" *
              "obj: _$(starting_letter)x\n" *
              "subject to\n" *
              "_$(starting_letter)c1: 1 _$(starting_letter)x >= -1\n" *
              "Bounds\n" *
              "_$(starting_letter)x free\n" *
              "End\n"
    end
end

function test_name_sanitization_illegal()
    for illegal_letter in ["[", "]", "*", "^"]
        model = LP.Model()
        MOI.Utilities.loadfromstring!(
            model,
            """
variables: x
minobjective: x
c1: 1.0 * x >= -1.0
""",
        )
        x = MOI.get(model, MOI.VariableIndex, "x")
        MOI.set(model, MOI.VariableName(), x, "x$(illegal_letter)y")
        c = MOI.get(model, MOI.ConstraintIndex, "c1")
        MOI.set(model, MOI.ConstraintName(), c, "c$(illegal_letter)d")
        MOI.write_to_file(model, LP_TEST_FILE)
        @test read(LP_TEST_FILE, String) ==
              "minimize\n" *
              "obj: x_y\n" *
              "subject to\n" *
              "c_d: 1 x_y >= -1\n" *
              "Bounds\n" *
              "x_y free\n" *
              "End\n"
    end
end

function test_name_sanitization_duplicate()
    model = LP.Model()
    MOI.Utilities.loadfromstring!(
        model,
        """
variables: a, b, c, d
minobjective: a + b + c + d
c1: a + b + c + d >= -1.0
""",
    )
    variables = MOI.get(model, MOI.ListOfVariableIndices())
    MOI.set.(model, MOI.VariableName(), variables, ["a[", "a]", "a*", "a^"])
    MOI.write_to_file(model, LP_TEST_FILE)
    @test read(LP_TEST_FILE, String) ==
          "minimize\n" *
          "obj: 1 a_ + 1 a__1 + 1 a__2 + 1 a__3\n" *
          "subject to\n" *
          "c1: 1 a_ + 1 a__1 + 1 a__2 + 1 a__3 >= -1\n" *
          "Bounds\n" *
          "a_ free\n" *
          "a__1 free\n" *
          "a__2 free\n" *
          "a__3 free\n" *
          "End\n"
end

function test_name_sanitization_length()
    model = LP.Model(maximum_length = 3)
    x = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "abcdefg")
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    MOI.write_to_file(model, LP_TEST_FILE)
    @test read(LP_TEST_FILE, String) ==
          "minimize\n" *
          "obj: abc\n" *
          "subject to\n" *
          "Bounds\n" *
          "abc free\n" *
          "End\n"
end

function test_name_sanitization_too_long()
    model = LP.Model(maximum_length = 3)
    MOI.Utilities.loadfromstring!(
        model,
        """
variables: abcd, abce
minobjective: abcd + abce
c1: abcd + abce >= -1.0
""",
    )
    MOI.write_to_file(model, LP_TEST_FILE)
    @test read(LP_TEST_FILE, String) ==
          "minimize\n" *
          "obj: 1 abc + 1 abc_1\n" *
          "subject to\n" *
          "c1: 1 abc + 1 abc_1 >= -1\n" *
          "Bounds\n" *
          "abc free\n" *
          "abc_1 free\n" *
          "End\n"
end

function test_name_sanitization_other()
    model = LP.Model()
    MOIU.loadfromstring!(
        model,
        """
variables: x
maxobjective: 2.0 * x + -1.0
""",
    )
    MOI.write_to_file(model, LP_TEST_FILE)
    @test read(LP_TEST_FILE, String) ==
          "maximize\n" *
          "obj: -1 + 2 x\n" *
          "subject to\n" *
          "Bounds\n" *
          "x free\n" *
          "End\n"
end

function test_free_variables()
    model = LP.Model()
    MOIU.loadfromstring!(
        model,
        """
variables: x, y, z
maxobjective: x
x in ZeroOne()
y in Integer()
""",
    )
    MOI.write_to_file(model, LP_TEST_FILE)
    @test read(LP_TEST_FILE, String) ==
          "maximize\n" *
          "obj: x\n" *
          "subject to\n" *
          "Bounds\n" *
          "y free\n" *
          "z free\n" *
          "General\n" *
          "y\n" *
          "Binary\n" *
          "x\n" *
          "End\n"
end

function test_quadratic_objective()
    model = LP.Model()
    @test_throws(
        MOI.UnsupportedAttribute,
        MOIU.loadfromstring!(
            model,
            """
variables: x
minobjective: 1.0*x*x
""",
        )
    )
end

function test_read_example_lo1()
    model = LP.Model()
    MOI.read_from_file(model, joinpath(@__DIR__, "models", "example_lo1.lp"))
    @test MOI.get(model, MOI.NumberOfVariables()) == 4
    constraints = MOI.get(model, MOI.ListOfConstraintTypesPresent())
    @test (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) in
          constraints
    @test (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) in
          constraints
    @test (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) in
          constraints
    @test (MOI.VariableIndex, MOI.GreaterThan{Float64}) in constraints
    @test (MOI.VariableIndex, MOI.Interval{Float64}) in constraints
    return nothing
end

function test_read_model2()
    model = LP.Model()
    MOI.read_from_file(model, joinpath(@__DIR__, "models", "model2.lp"))
    @test MOI.get(model, MOI.NumberOfVariables()) == 8
    constraints = MOI.get(model, MOI.ListOfConstraintTypesPresent())
    @test (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) in
          constraints
    @test (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) in
          constraints
    @test (MOI.VariableIndex, MOI.GreaterThan{Float64}) in constraints
    @test (MOI.VariableIndex, MOI.Interval{Float64}) in constraints
    @test (MathOptInterface.VariableIndex, MathOptInterface.Integer) in
          constraints
    @test (MathOptInterface.VariableIndex, MathOptInterface.ZeroOne) in
          constraints
    # Adicionar testes dos bounds de V8
    @test MOI.get(model, MOI.VariableName(), MOI.VariableIndex(8)) == "V8"
    @test model.variables.lower[8] == -Inf
    @test model.variables.upper[8] == -3
    obj_type = MOI.get(model, MOI.ObjectiveFunctionType())
    obj_func = MOI.get(model, MOI.ObjectiveFunction{obj_type}())
    @test obj_func.constant == 2.5
    return nothing
end

function test_read_model1_tricky()
    model = LP.Model()
    MOI.read_from_file(model, joinpath(@__DIR__, "models", "model1_tricky.lp"))
    @test MOI.get(model, MOI.NumberOfVariables()) == 8
    var_names = MOI.get.(model, MOI.VariableName(), MOI.VariableIndex.(1:8))
    @test Set(var_names) ==
          Set(["Var4", "V5", "V1", "V2", "V3", "V6", "V7", "V8"])
    return nothing
end

function test_read_corrupt()
    model = LP.Model()
    @test_throws ErrorException MOI.read_from_file(
        model,
        joinpath(@__DIR__, "models", "corrupt.lp"),
    )
    return nothing
end

function test_read_invalid_variable_name()
    model = LP.Model()
    @test_throws ErrorException MOI.read_from_file(
        model,
        joinpath(@__DIR__, "models", "invalid_variable_name.lp"),
    )
    return nothing
end

function test_read_invalid_affine_term_objective()
    model = LP.Model()
    @test_throws ErrorException MOI.read_from_file(
        model,
        joinpath(@__DIR__, "models", "invalid_affine_term_objective.lp"),
    )
    return nothing
end

function test_read_invalid_affine_term_constraint()
    model = LP.Model()
    @test_throws ErrorException MOI.read_from_file(
        model,
        joinpath(@__DIR__, "models", "invalid_affine_term_constraint.lp"),
    )
    return nothing
end

function test_read_invalid_sos_set()
    model = LP.Model()
    @test_throws ErrorException MOI.read_from_file(
        model,
        joinpath(@__DIR__, "models", "invalid_sos_set.lp"),
    )
    return nothing
end

function test_read_invalid_sos_constraint()
    model = LP.Model()
    @test_throws ErrorException MOI.read_from_file(
        model,
        joinpath(@__DIR__, "models", "invalid_sos_constraint.lp"),
    )
    return nothing
end

function test_read_invalid_bound()
    model = LP.Model()
    @test_throws ErrorException MOI.read_from_file(
        model,
        joinpath(@__DIR__, "models", "invalid_bound.lp"),
    )
    return nothing
end

function test_read_invalid_constraint()
    model = LP.Model()
    @test_throws ErrorException MOI.read_from_file(
        model,
        joinpath(@__DIR__, "models", "invalid_constraint.lp"),
    )
    return nothing
end

function test_read_model1()
    model = LP.Model()
    MOI.read_from_file(model, joinpath(@__DIR__, "models", "model1.lp"))
    constraints = MOI.get(model, MOI.ListOfConstraintTypesPresent())
    @test (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) in
          constraints
    @test (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) in
          constraints
    @test (MOI.VariableIndex, MOI.GreaterThan{Float64}) in constraints
    @test (MOI.VariableIndex, MOI.Interval{Float64}) in constraints
    @test (MathOptInterface.VariableIndex, MathOptInterface.Integer) in
          constraints
    @test (MathOptInterface.VariableIndex, MathOptInterface.ZeroOne) in
          constraints
    @test (
        MathOptInterface.VectorOfVariables,
        MathOptInterface.SOS1{Float64},
    ) in constraints
    @test (
        MathOptInterface.VectorOfVariables,
        MathOptInterface.SOS2{Float64},
    ) in constraints
    return nothing
end

function runtests()
    for name in names(@__MODULE__, all = true)
        if startswith("$(name)", "test_")
            @testset "$name" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    sleep(1.0)  # Allow time for unlink to happen.
    rm(LP_TEST_FILE, force = true)
    return
end

end

TestLP.runtests()


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

"""
    GeoTable{𝒯}

A geospatial table where the underlying table of type `𝒯`
is converted lazily row-by-row to contain geometries from
Meshes.jl for computational geometry in pure Julia.

This table type implements the `Meshes.Data` trait and is
therefore compatible with the GeoStats.jl ecosystem.
"""
struct GeoTable{𝒯} <: Meshes.Data
  table::𝒯
end

function Meshes.domain(t::GeoTable)
  table = getfield(t, :table)
  gcol  = geomcolumn(table)
  geoms = Tables.getcolumn(table, gcol)
  items = geom2meshes.(geoms)
  Meshes.Collection(items)
end

function Meshes.values(t::GeoTable, rank=nothing)
  # find ranks of all geometries
  table = getfield(t, :table)
  gcol  = geomcolumn(table)
  geoms = Tables.getcolumn(table, gcol)
  items = geom2meshes.(geoms)
  ranks = paramdim.(items)

  # select geometries with given rank
  rmax  = maximum(ranks)
  rsel  = isnothing(rank) ? rmax : rank
  rind  = findall(==(rsel), ranks)

  # check if rank exists in data
  if isempty(rind)
    nothing
  else
    # if rank exists, load other columns
    rows = Tables.rows(table)
    sche = Tables.schema(rows)
    vars = setdiff(sche.names, [gcol])
    cols = map(vars) do var
      col = Tables.getcolumn(table, var)
      var => col[rind]
    end
    (; cols...)
  end
end

# helper function to find the geometry column of a table
function geomcolumn(table)
  rows = Tables.rows(table)
  sche = Tables.schema(rows)
  if :geometry ∈ sche.names
    :geometry
  elseif :geom ∈ sche.names
    :geom
  else
    throw(ErrorException("geometry column not found"))
  end
end


text:
#=
The diffusion matrix insert noise to deterministic model.
=#

"""
Generate Diffusion matrix.
For independent noise, the diffusion matrix is an identity matrix

Arguement
---------
- `len_of_u`{Integer}: the dimension of ODE model
"""
function DiffusionMatrix(len_of_u)
    return Matrix(I, len_of_u, len_of_u)
end


text:
#=
        project1.jl -- This is where the magic happens!

    All of your code must either live in this file, or be `include`d here.
=#

#=
    If you want to use packages, please do so up here.
    Note that you may use any packages in the julia standard library
    (i.e. ones that ship with the julia language) as well as Statistics
    (since we use it in the backend already anyway)
=#

# Example:
# using LinearAlgebra

#=
    If you're going to include files, please do so up here. Note that they
    must be saved in project1_jl and you must use the relative path
    (not the absolute path) of the file in the include statement.

    [Good]  include("somefile.jl")
    [Bad]   include("/pathto/project1_jl/somefile.jl")
=#

# Example
# include("myfile.jl")

"""
    optimize(f, g, x0, n, prob)

Arguments:
    - `f`: Function to be optimized
    - `g`: Gradient function for `f`
    - `x0`: (Vector) Initial position to start from
    - `n`: (Int) Number of evaluations allowed. Remember `g` costs twice of `f`
    - `prob`: (String) Name of the problem. So you can use a different strategy for each problem. E.g. "simple1", "secret2", etc.

Returns:
    - The location of the minimum
"""
function optimize(f, g, x0, n, prob)
    x_best = x0
    return x_best
end