text:
"""
`munsell_specs()`

Returns an array containing all 1269 munsell spectra, measured from 380 to 800 nm in 1 nm steps.
"""
function munsell_specs()
    filename = "munsell.dat"
    fp=@__DIR__
    cd(fp)
    cd("..\\spectra")
    sfile=readlines(filename)
    map(x->parse(Float64,x),sfile)
end

"""
`munsell_spec(idx::Int, speclist)`

Returns reflectance spectrum of the Munsell dataset at index `idx`.
Possible indices go from 1 to 1269

# Example

```jldoctest
julia> munsell_spec(437, speclist)
RSpec(Real[380, 381, 382, 383, 384, 385, 386, 387, 388, 389  …  791, 792, 793, 794, 795, 796, 797, 798, 799, 800], Real[0.1184, 0.1208, 0.124, 0.1301, 0.1349, 0.1431, 0.1498, 0.1579, 0.1631, 0.1731  …  0.3355, 0.3364, 0.3345, 0.3337, 0.3327, 0.3348, 0.3392, 0.3377, 0.335, 0.3351])
```
"""
function munsell_spec(idx::Int, speclist = munsell_specs())
    n=421
    0 < idx < 1270 ? (i=(idx-1)*n+1) : throw(DomainError(idx, "Index is outside the range of 1 ≤ idx ≤ 1269."))
    y=speclist[i:i+n-1]
    RSpec(collect(380:800),y)
end

function munsell_spec(munsellcolor::AbstractString)
    idx = get(munsell_dict(), munsellcolor, 9999)
    ifelse(idx == 9999, throw(DomainError(idx, "Munsell Color does not exist")), munsell_spec(idx, munsell_specs()))
end

function spec_color(col::AbstractString,colordict::AbstractString)
    ind = get(munsell_dict(), col, "Color key $(col) does not exist")
    if haskey(td,col)
        return spec_idx(ind, munsell_specs())
    else throw(DomainError(ind, "Color key does not exist"))
    end
end

function munsell_dict()
    d = Dict()
    n=length(munsell_keys)
    for i in 1:n
        push!(d, munsell_keys[i] => i)
    end
    d
end

function plotmunsell(index)
    λ=collect(380:800)
    n=421
    plot(λ, zeros(n))
    plot!(λ, spec_idx(index, munsell_specs())[i:i+n-1])
end

# The full set of 1269 munsell color names
const munsell_keys=[
    "2.5R 9/2";
    "2.5R 8/2";
    "2.5R 7/2";
    "2.5R 6/2";
    "2.5R 5/2";
    "2.5R 4/2";
    "2.5R 3/2";
    "2.5R 2.5/2";
    "2.5R 8/4";
    "2.5R 7/4";
    "2.5R 6/4";
    "2.5R 5/4";
    "2.5R 4/4";
    "2.5R 3/4";
    "2.5R 7/6";
    "2.5R 6/6";
    "2.5R 5/6";
    "2.5R 4/6";
    "2.5R 3/6";
    "2.5R 7/8";
    "2.5R 6/8";
    "2.5R 5/8";
    "2.5R 4/8";
    "2.5R 7/10";
    "2.5R 6/10";
    "2.5R 5/10";
    "2.5R 4/10";
    "2.5R 6/12";
    "2.5R 5/12";
    "2.5R 4/12";
    "5R 9/1";
    "5R 8/1";
    "5R 7/1";
    "5R 6/1";
    "5R 5/1";
    "5R 4/1";
    "5R 3/1";
    "5R 2.5/1";
    "5R 9/2";
    "5R 8/2";
    "5R 7/2";
    "5R 6/2";
    "5R 5/2";
    "5R 4/2";
    "5R 3/2";
    "5R 2.5/2";
    "5R 8/4";
    "5R 7/4";
    "5R 6/4";
    "5R 5/4";
    "5R 4/4";
    "5R 3/4";
    "5R 7/6";
    "5R 6/6";
    "5R 5/6";
    "5R 4/6";
    "5R 3/6";
    "5R 7/8";
    "5R 6/8";
    "5R 5/8";
    "5R 4/8";
    "5R 7/10";
    "5R 6/10";
    "5R 5/10";
    "5R 4/10";
    "5R 6/12";
    "5R 5/12";
    "5R 4/12";
    "5R 6/14";
    "5R 5/14";
    "5R 4/14";
    "7.5R 9/2";
    "7.5R 8/2";
    "7.5R 7/2";
    "7.5R 6/2";
    "7.5R 5/2";
    "7.5R 4/2";
    "7.5R 3/2";
    "7.5R 2.5/2";
    "7.5R 8/4";
    "7.5R 7/4";
    "7.5R 6/4";
    "7.5R 5/4";
    "7.5R 4/4";
    "7.5R 3/4";
    "7.5R 7/6";
    "7.5R 6/6";
    "7.5R 5/6";
    "7.5R 4/6";
    "7.5R 3/6";
    "7.5R 7/8";
    "7.5R 6/8";
    "7.5R 5/8";
    "7.5R 4/8";
    "7.5R 7/10";
    "7.5R 6/10";
    "7.5R 5/10";
    "7.5R 4/10";
    "7.5R 6/12";
    "7.5R 5/12";
    "7.5R 4/12";
    "10R 9/1";
    "10R 8/1";
    "10R 7/1";
    "10R 6/1";
    "10R 5/1";
    "10R 4/1";
    "10R 3/1";
    "10R 2.5/1";
    "10R 9/2";
    "10R 8/2";
    "10R 7/2";
    "10R 6/2";
    "10R 5/2";
    "10R 4/2";
    "10R 3/2";
    "10R 2.5/2";
    "10R 8/4";
    "10R 7/4";
    "10R 6/4";
    "10R 5/4";
    "10R 4/4";
    "10R 3/4";
    "10R 8/6";
    "10R 7/6";
    "10R 6/6";
    "10R 5/6";
    "10R 4/6";
    "10R 8/8";
    "10R 7/8";
    "10R 6/8";
    "10R 5/8";
    "10R 4/8";
    "10R 7/10";
    "10R 6/10";
    "10R 5/10";
    "10R 7/12";
    "10R 6/12";
    "10R 5/12";
    "2.5YR 9/2";
    "2.5YR 8/2";
    "2.5YR 7/2";
    "2.5YR 6/2";
    "2.5YR 5/2";
    "2.5YR 4/2";
    "2.5YR 3/2";
    "2.5YR 2.5/2";
    "2.5YR 8/4";
    "2.5YR 7/4";
    "2.5YR 6/4";
    "2.5YR 5/4";
    "2.5YR 4/4";
    "2.5YR 8/6";
    "2.5YR 7/6";
    "2.5YR 6/6";
    "2.5YR 5/6";
    "2.5YR 4/6";
    "2.5YR 8/8";
    "2.5YR 7/8";
    "2.5YR 6/8";
    "2.5YR 5/8";
    "2.5YR 4/8";
    "2.5YR 7/10";
    "2.5YR 6/10";
    "2.5YR 5/10";
    "2.5YR 7/12";
    "2.5YR 6/12";
    "2.5YR 6/14";
    "5YR 9/1";
    "5YR 8/1";
    "5YR 7/1";
    "5YR 6/1";
    "5YR 5/1";
    "5YR 4/1";
    "5YR 3/1";
    "5YR 2.5/1";
    "5YR 9/2";
    "5YR 8/2";
    "5YR 7/2";
    "5YR 6/2";
    "5YR 5/2";
    "5YR 4/2";
    "5YR 3/2";
    "5YR 9/4";
    "5YR 8/4";
    "5YR 7/4";
    "5YR 6/4";
    "5YR 5/4";
    "5YR 4/4";
    "5YR 8/6";
    "5YR 7/6";
    "5YR 6/6";
    "5YR 5/6";
    "5YR 4/6";
    "5YR 8/8";
    "5YR 7/8";
    "5YR 6/8";
    "5YR 5/8";
    "5YR 7/10";
    "5YR 6/10";
    "5YR 7/12";
    "5YR 6/12";
    "7.5YR 9/2";
    "7.5YR 8/2";
    "7.5YR 7/2";
    "7.5YR 6/2";
    "7.5YR 5/2";
    "7.5YR 4/2";
    "7.5YR 3/2";
    "7.5YR 9/4";
    "7.5YR 8/4";
    "7.5YR 7/4";
    "7.5YR 6/4";
    "7.5YR 5/4";
    "7.5YR 4/4";
    "7.5YR 8/6";
    "7.5YR 7/6";
    "7.5YR 6/6";
    "7.5YR 5/6";
    "7.5YR 4/6";
    "7.5YR 8/8";
    "7.5YR 7/8";
    "7.5YR 6/8";
    "7.5YR 5/8";
    "7.5YR 8/10";
    "7.5YR 7/10";
    "7.5YR 6/10";
    "7.5YR 7/12";
    "10YR 9/1";
    "10YR 8/1";
    "10YR 7/1";
    "10YR 6/1";
    "10YR 5/1";
    "10YR 4/1";
    "10YR 3/1";
    "10YR 2.5/1";
    "10YR 9/2";
    "10YR 8/2";
    "10YR 7/2";
    "10YR 6/2";
    "10YR 5/2";
    "10YR 4/2";
    "10YR 3/2";
    "10YR 9/4";
    "10YR 8/4";
    "10YR 7/4";
    "10YR 6/4";
    "10YR 5/4";
    "10YR 4/4";
    "10YR 8/6";
    "10YR 7/6";
    "10YR 6/6";
    "10YR 5/6";
    "10YR 8/8";
    "10YR 7/8";
    "10YR 6/8";
    "10YR 5/8";
    "10YR 8/10";
    "10YR 7/10";
    "10YR 6/10";
    "10YR 7/12";
    "2.5Y 9/2";
    "2.5Y 8.5/2";
    "2.5Y 8/2";
    "2.5Y 7/2";
    "2.5Y 6/2";
    "2.5Y 5/2";
    "2.5Y 4/2";
    "2.5Y 3/2";
    "2.5Y 9/4";
    "2.5Y 8.5/4";
    "2.5Y 8/4";
    "2.5Y 7/4";
    "2.5Y 6/4";
    "2.5Y 5/4";
    "2.5Y 4/4";
    "2.5Y 9/6";
    "2.5Y 8.5/6";
    "2.5Y 8/6";
    "2.5Y 7/6";
    "2.5Y 6/6";
    "2.5Y 5/6";
    "2.5Y 8.5/8";
    "2.5Y 8/8";
    "2.5Y 7/8";
    "2.5Y 6/8";
    "2.5Y 8.5/10";
    "2.5Y 8/10";
    "2.5Y 7/10";
    "2.5Y 8/12";
    "2.5Y 7/12";
    "5Y 9/1";
    "5Y 8.5/1";
    "5Y 8/1";
    "5Y 7/1";
    "5Y 6/1";
    "5Y 5/1";
    "5Y 4/1";
    "5Y 3/1";
    "5Y 2.5/1";
    "5Y 9/2";
    "5Y 8.5/2";
    "5Y 8/2";
    "5Y 7/2";
    "5Y 6/2";
    "5Y 5/2";
    "5Y 4/2";
    "5Y 3/2";
    "5Y 9/4";
    "5Y 8.5/4";
    "5Y 8/4";
    "5Y 7/4";
    "5Y 6/4";
    "5Y 5/4";
    "5Y 4/4";
    "5Y 9/6";
    "5Y 8.5/6";
    "5Y 8/6";
    "5Y 7/6";
    "5Y 6/6";
    "5Y 5/6";
    "5Y 9/8";
    "5Y 8.5/8";
    "5Y 8/8";
    "5Y 7/8";
    "5Y 6/8";
    "5Y 8.5/10";
    "5Y 8/10";
    "5Y 7/10";
    "5Y 8.5/12";
    "5Y 8/12";
    "5Y 7/12";
    "7.5Y 9/2";
    "7.5Y 8.5/2";
    "7.5Y 8/2";
    "7.5Y 7/2";
    "7.5Y 6/2";
    "7.5Y 5/2";
    "7.5Y 4/2";
    "7.5Y 3/2";
    "7.5Y 9/4";
    "7.5Y 8.5/4";
    "7.5Y 8/4";
    "7.5Y 7/4";
    "7.5Y 6/4";
    "7.5Y 5/4";
    "7.5Y 4/4";
    "7.5Y 9/6";
    "7.5Y 8.5/6";
    "7.5Y 8/6";
    "7.5Y 7/6";
    "7.5Y 6/6";
    "7.5Y 5/6";
    "7.5Y 9/8";
    "7.5Y 8.5/8";
    "7.5Y 8/8";
    "7.5Y 7/8";
    "7.5Y 6/8";
    "7.5Y 9/10";
    "7.5Y 8.5/10";
    "7.5Y 8/10";
    "7.5Y 7/10";
    "7.5Y 8.5/12";
    "7.5Y 8/12";
    "10Y 9/1";
    "10Y 8.5/1";
    "10Y 8/1";
    "10Y 7/1";
    "10Y 6/1";
    "10Y 5/1";
    "10Y 4/1";
    "10Y 3/1";
    "10Y 2.5/1";
    "10Y 9/2";
    "10Y 8.5/2";
    "10Y 8/2";
    "10Y 7/2";
    "10Y 6/2";
    "10Y 5/2";
    "10Y 4/2";
    "10Y 3/2";
    "10Y 9/4";
    "10Y 8.5/4";
    "10Y 8/4";
    "10Y 7/4";
    "10Y 6/4";
    "10Y 5/4";
    "10Y 4/4";
    "10Y 9/6";
    "10Y 8.5/6";
    "10Y 8/6";
    "10Y 7/6";
    "10Y 6/6";
    "10Y 5/6";
    "10Y 9/8";
    "10Y 8.5/8";
    "10Y 8/8";
    "10Y 7/8";
    "10Y 6/8";
    "10Y 9/10";
    "10Y 8.5/10";
    "10Y 8/10";
    "10Y 7/10";
    "10Y 8/12";
    "2.5GY 9/2";
    "2.5GY 8.5/2";
    "2.5GY 8/2";
    "2.5GY 7/2";
    "2.5GY 6/2";
    "2.5GY 5/2";
    "2.5GY 4/2";
    "2.5GY 3/2";
    "2.5GY 9/4";
    "2.5GY 8.5/4";
    "2.5GY 8/4";
    "2.5GY 7/4";
    "2.5GY 6/4";
    "2.5GY 5/4";
    "2.5GY 4/4";
    "2.5GY 9/6";
    "2.5GY 8.5/6";
    "2.5GY 8/6";
    "2.5GY 7/6";
    "2.5GY 6/6";
    "2.5GY 5/6";
    "2.5GY 9/8";
    "2.5GY 8.5/8";
    "2.5GY 8/8";
    "2.5GY 7/8";
    "2.5GY 6/8";
    "2.5GY 8.5/10";
    "2.5GY 8/10";
    "2.5GY 7/10";
    "5GY 9/1";
    "5GY 8.5/1";
    "5GY 8/1";
    "5GY 7/1";
    "5GY 6/1";
    "5GY 5/1";
    "5GY 4/1";
    "5GY 3/1";
    "5GY 2.5/1";
    "5GY 9/2";
    "5GY 8.5/2";
    "5GY 8/2";
    "5GY 7/2";
    "5GY 6/2";
    "5GY 5/2";
    "5GY 4/2";
    "5GY 3/2";
    "5GY 9/4";
    "5GY 8.5/4";
    "5GY 8/4";
    "5GY 7/4";
    "5GY 6/4";
    "5GY 5/4";
    "5GY 4/4";
    "5GY 9/6";
    "5GY 8.5/6";
    "5GY 8/6";
    "5GY 7/6";
    "5GY 6/6";
    "5GY 5/6";
    "5GY 8.5/8";
    "5GY 8/8";
    "5GY 7/8";
    "5GY 6/8";
    "5GY 5/8";
    "5GY 8.5/10";
    "5GY 8/10";
    "7.5GY 9/2";
    "7.5GY 8.5/2";
    "7.5GY 8/2";
    "7.5GY 7/2";
    "7.5GY 6/2";
    "7.5GY 5/2";
    "7.5GY 4/2";
    "7.5GY 3/2";
    "7.5GY 9/4";
    "7.5GY 8.5/4";
    "7.5GY 8/4";
    "7.5GY 7/4";
    "7.5GY 6/4";
    "7.5GY 5/4";
    "7.5GY 4/4";
    "7.5GY 8.5/6";
    "7.5GY 8/6";
    "7.5GY 7/6";
    "7.5GY 6/6";
    "7.5GY 5/6";
    "7.5GY 8/8";
    "7.5GY 7/8";
    "7.5GY 6/8";
    "7.5GY 5/8";
    "7.5GY 7/10";
    "7.5GY 6/10";
    "10GY 9/1";
    "10GY 8.5/1";
    "10GY 8/1";
    "10GY 7/1";
    "10GY 6/1";
    "10GY 5/1";
    "10GY 4/1";
    "10GY 3/1";
    "10GY 2.5/1";
    "10GY 9/2";
    "10GY 8.5/2";
    "10GY 8/2";
    "10GY 7/2";
    "10GY 6/2";
    "10GY 5/2";
    "10GY 4/2";
    "10GY 3/2";
    "10GY 2.5/2";
    "10GY 9/4";
    "10GY 8.5/4";
    "10GY 8/4";
    "10GY 7/4";
    "10GY 6/4";
    "10GY 5/4";
    "10GY 4/4";
    "10GY 8.5/6";
    "10GY 8/6";
    "10GY 7/6";
    "10GY 6/6";
    "10GY 5/6";
    "10GY 7/8";
    "10GY 6/8";
    "10GY 5/8";
    "10GY 7/10";
    "10GY 6/10";
    "2.5G 9/2";
    "2.5G 8/2";
    "2.5G 7/2";
    "2.5G 6/2";
    "2.5G 5/2";
    "2.5G 4/2";
    "2.5G 3/2";
    "2.5G 2.5/2";
    "2.5G 8/4";
    "2.5G 7/4";
    "2.5G 6/4";
    "2.5G 5/4";
    "2.5G 4/4";
    "2.5G 3/4";
    "2.5G 8/6";
    "2.5G 7/6";
    "2.5G 6/6";
    "2.5G 5/6";
    "2.5G 4/6";
    "2.5G 7/8";
    "2.5G 6/8";
    "2.5G 5/8";
    "2.5G 7/10";
    "2.5G 6/10";
    "2.5G 5/10";
    "2.5G 6/12";
    "5G 9/1";
    "5G 8/1";
    "5G 7/1";
    "5G 6/1";
    "5G 5/1";
    "5G 4/1";
    "5G 3/1";
    "5G 2.5/1";
    "5G 9/2";
    "5G 8/2";
    "5G 7/2";
    "5G 6/2";
    "5G 5/2";
    "5G 4/2";
    "5G 3/2";
    "5G 2.5/2";
    "5G 8/4";
    "5G 7/4";
    "5G 6/4";
    "5G 5/4";
    "5G 4/4";
    "5G 3/4";
    "5G 8/6";
    "5G 7/6";
    "5G 6/6";
    "5G 5/6";
    "5G 4/6";
    "5G 7/8";
    "5G 6/8";
    "5G 5/8";
    "5G 4/8";
    "5G 7/10";
    "5G 6/10";
    "7.5G 9/2";
    "7.5G 8/2";
    "7.5G 7/2";
    "7.5G 6/2";
    "7.5G 5/2";
    "7.5G 4/2";
    "7.5G 3/2";
    "7.5G 2.5/2";
    "7.5G 8/4";
    "7.5G 7/4";
    "7.5G 6/4";
    "7.5G 5/4";
    "7.5G 4/4";
    "7.5G 3/4";
    "7.5G 8/6";
    "7.5G 7/6";
    "7.5G 6/6";
    "7.5G 5/6";
    "7.5G 4/6";
    "7.5G 7/8";
    "7.5G 6/8";
    "7.5G 5/8";
    "7.5G 4/8";
    "7.5G 7/10";
    "7.5G 6/10";
    "10G 9/1";
    "10G 8/1";
    "10G 7/1";
    "10G 6/1";
    "10G 5/1";
    "10G 4/1";
    "10G 3/1";
    "10G 2.5/1";
    "10G 9/2";
    "10G 8/2";
    "10G 7/2";
    "10G 6/2";
    "10G 5/2";
    "10G 4/2";
    "10G 3/2";
    "10G 2.5/2";
    "10G 8/4";
    "10G 7/4";
    "10G 6/4";
    "10G 5/4";
    "10G 4/4";
    "10G 3/4";
    "10G 7/6";
    "10G 6/6";
    "10G 5/6";
    "10G 4/6";
    "10G 7/8";
    "10G 6/8";
    "10G 5/8";
    "10G 4/8";
    "10G 6/10";
    "2.5BG 9/2";
    "2.5BG 8/2";
    "2.5BG 7/2";
    "2.5BG 6/2";
    "2.5BG 5/2";
    "2.5BG 4/2";
    "2.5BG 3/2";
    "2.5BG 2.5/2";
    "2.5BG 8/4";
    "2.5BG 7/4";
    "2.5BG 6/4";
    "2.5BG 5/4";
    "2.5BG 4/4";
    "2.5BG 3/4";
    "2.5BG 7/6";
    "2.5BG 6/6";
    "2.5BG 5/6";
    "2.5BG 4/6";
    "2.5BG 3/6";
    "2.5BG 7/8";
    "2.5BG 6/8";
    "2.5BG 5/8";
    "2.5BG 4/8";
    "2.5BG 6/10";
    "5BG 9/1";
    "5BG 8/1";
    "5BG 7/1";
    "5BG 6/1";
    "5BG 5/1";
    "5BG 4/1";
    "5BG 3/1";
    "5BG 2.5/1";
    "5BG 9/2";
    "5BG 8/2";
    "5BG 7/2";
    "5BG 6/2";
    "5BG 5/2";
    "5BG 4/2";
    "5BG 3/2";
    "5BG 2.5/2";
    "5BG 8/4";
    "5BG 7/4";
    "5BG 6/4";
    "5BG 5/4";
    "5BG 4/4";
    "5BG 3/4";
    "5BG 7/6";
    "5BG 6/6";
    "5BG 5/6";
    "5BG 4/6";
    "5BG 7/8";
    "5BG 6/8";
    "5BG 5/8";
    "5BG 4/8";
    "7.5BG 9/2";
    "7.5BG 8/2";
    "7.5BG 7/2";
    "7.5BG 6/2";
    "7.5BG 5/2";
    "7.5BG 4/2";
    "7.5BG 3/2";
    "7.5BG 2.5/2";
    "7.5BG 8/4";
    "7.5BG 7/4";
    "7.5BG 6/4";
    "7.5BG 5/4";
    "7.5BG 4/4";
    "7.5BG 3/4";
    "7.5BG 7/6";
    "7.5BG 6/6";
    "7.5BG 5/6";
    "7.5BG 4/6";
    "7.5BG 7/8";
    "7.5BG 6/8";
    "7.5BG 5/8";
    "7.5BG 4/8";
    "10BG 9/1";
    "10BG 8/1";
    "10BG 7/1";
    "10BG 6/1";
    "10BG 5/1";
    "10BG 4/1";
    "10BG 3/1";
    "10BG 2.5/1";
    "10BG 9/2";
    "10BG 8/2";
    "10BG 7/2";
    "10BG 6/2";
    "10BG 5/2";
    "10BG 4/2";
    "10BG 3/2";
    "10BG 2.5/2";
    "10BG 8/4";
    "10BG 7/4";
    "10BG 6/4";
    "10BG 5/4";
    "10BG 4/4";
    "10BG 3/4";
    "10BG 7/6";
    "10BG 6/6";
    "10BG 5/6";
    "10BG 4/6";
    "10BG 7/8";
    "10BG 6/8";
    "10BG 5/8";
    "10BG 4/8";
    "2.5B 9/2";
    "2.5B 8/2";
    "2.5B 7/2";
    "2.5B 6/2";
    "2.5B 5/2";
    "2.5B 4/2";
    "2.5B 3/2";
    "2.5B 2.5/2";
    "2.5B 8/4";
    "2.5B 7/4";
    "2.5B 6/4";
    "2.5B 5/4";
    "2.5B 4/4";
    "2.5B 3/4";
    "2.5B 7/6";
    "2.5B 6/6";
    "2.5B 5/6";
    "2.5B 4/6";
    "2.5B 7/8";
    "2.5B 6/8";
    "2.5B 5/8";
    "2.5B 4/8";
    "5B 9/1";
    "5B 8/1";
    "5B 7/1";
    "5B 6/1";
    "5B 5/1";
    "5B 4/1";
    "5B 3/1";
    "5B 2.5/1";
    "5B 9/2";
    "5B 8/2";
    "5B 7/2";
    "5B 6/2";
    "5B 5/2";
    "5B 4/2";
    "5B 3/2";
    "5B 2.5/2";
    "5B 8/4";
    "5B 7/4";
    "5B 6/4";
    "5B 5/4";
    "5B 4/4";
    "5B 3/4";
    "5B 7/6";
    "5B 6/6";
    "5B 5/6";
    "5B 4/6";
    "5B 3/6";
    "5B 7/8";
    "5B 6/8";
    "5B 5/8";
    "5B 4/8";
    "7.5B 9/2";
    "7.5B 8/2";
    "7.5B 7/2";
    "7.5B 6/2";
    "7.5B 5/2";
    "7.5B 4/2";
    "7.5B 3/2";
    "7.5B 2.5/2";
    "7.5B 8/4";
    "7.5B 7/4";
    "7.5B 6/4";
    "7.5B 5/4";
    "7.5B 4/4";
    "7.5B 3/4";
    "7.5B 7/6";
    "7.5B 6/6";
    "7.5B 5/6";
    "7.5B 4/6";
    "7.5B 3/6";
    "7.5B 7/8";
    "7.5B 6/8";
    "7.5B 5/8";
    "7.5B 4/8";
    "7.5B 6/10";
    "7.5B 5/10";
    "10B 9/1";
    "10B 8/1";
    "10B 7/1";
    "10B 6/1";
    "10B 5/1";
    "10B 4/1";
    "10B 3/1";
    "10B 2.5/1";
    "10B 9/2";
    "10B 8/2";
    "10B 7/2";
    "10B 6/2";
    "10B 5/2";
    "10B 4/2";
    "10B 3/2";
    "10B 2.5/2";
    "10B 8/4";
    "10B 7/4";
    "10B 6/4";
    "10B 5/4";
    "10B 4/4";
    "10B 3/4";
    "10B 2.5/4";
    "10B 7/6";
    "10B 6/6";
    "10B 5/6";
    "10B 4/6";
    "10B 3/6";
    "10B 7/8";
    "10B 6/8";
    "10B 5/8";
    "10B 4/8";
    "10B 6/10";
    "10B 5/10";
    "2.5PB 9/2";
    "2.5PB 8/2";
    "2.5PB 7/2";
    "2.5PB 6/2";
    "2.5PB 5/2";
    "2.5PB 4/2";
    "2.5PB 3/2";
    "2.5PB 2.5/2";
    "2.5PB 8/4";
    "2.5PB 7/4";
    "2.5PB 6/4";
    "2.5PB 5/4";
    "2.5PB 4/4";
    "2.5PB 3/4";
    "2.5PB 2.5/4";
    "2.5PB 8/6";
    "2.5PB 7/6";
    "2.5PB 6/6";
    "2.5PB 5/6";
    "2.5PB 4/6";
    "2.5PB 3/6";
    "2.5PB 7/8";
    "2.5PB 6/8";
    "2.5PB 5/8";
    "2.5PB 4/8";
    "2.5PB 6/10";
    "2.5PB 5/10";
    "2.5PB 4/10";
    "5PB 9/1";
    "5PB 8/1";
    "5PB 7/1";
    "5PB 6/1";
    "5PB 5/1";
    "5PB 4/1";
    "5PB 3/1";
    "5PB 2.5/1";
    "5PB 9/2";
    "5PB 8/2";
    "5PB 7/2";
    "5PB 6/2";
    "5PB 5/2";
    "5PB 4/2";
    "5PB 3/2";
    "5PB 2.5/2";
    "5PB 8/4";
    "5PB 7/4";
    "5PB 6/4";
    "5PB 5/4";
    "5PB 4/4";
    "5PB 3/4";
    "5PB 2.5/4";
    "5PB 8/6";
    "5PB 7/6";
    "5PB 6/6";
    "5PB 5/6";
    "5PB 4/6";
    "5PB 3/6";
    "5PB 7/8";
    "5PB 6/8";
    "5PB 5/8";
    "5PB 4/8";
    "5PB 3/8";
    "5PB 6/10";
    "5PB 5/10";
    "5PB 4/10";
    "5PB 5/12";
    "5PB 4/12";
    "7.5PB 9/2";
    "7.5PB 8/2";
    "7.5PB 7/2";
    "7.5PB 6/2";
    "7.5PB 5/2";
    "7.5PB 4/2";
    "7.5PB 3/2";
    "7.5PB 2.5/2";
    "7.5PB 8/4";
    "7.5PB 7/4";
    "7.5PB 6/4";
    "7.5PB 5/4";
    "7.5PB 4/4";
    "7.5PB 3/4";
    "7.5PB 2.5/4";
    "7.5PB 8/6";
    "7.5PB 7/6";
    "7.5PB 6/6";
    "7.5PB 5/6";
    "7.5PB 4/6";
    "7.5PB 3/6";
    "7.5PB 2.5/6";
    "7.5PB 7/8";
    "7.5PB 6/8";
    "7.5PB 5/8";
    "7.5PB 4/8";
    "7.5PB 3/8";
    "7.5PB 6/10";
    "7.5PB 5/10";
    "7.5PB 4/10";
    "7.5PB 3/10";
    "7.5PB 5/12";
    "7.5PB 4/12";
    "10PB 9/1";
    "10PB 8/1";
    "10PB 7/1";
    "10PB 6/1";
    "10PB 5/1";
    "10PB 4/1";
    "10PB 3/1";
    "10PB 2.5/1";
    "10PB 9/2";
    "10PB 8/2";
    "10PB 7/2";
    "10PB 6/2";
    "10PB 5/2";
    "10PB 4/2";
    "10PB 3/2";
    "10PB 2.5/2";
    "10PB 8/4";
    "10PB 7/4";
    "10PB 6/4";
    "10PB 5/4";
    "10PB 4/4";
    "10PB 3/4";
    "10PB 2.5/4";
    "10PB 7/6";
    "10PB 6/6";
    "10PB 5/6";
    "10PB 4/6";
    "10PB 3/6";
    "10PB 2.5/6";
    "10PB 7/8";
    "10PB 6/8";
    "10PB 5/8";
    "10PB 4/8";
    "10PB 3/8";
    "10PB 6/10";
    "10PB 5/10";
    "10PB 4/10";
    "2.5P 9/2";
    "2.5P 8/2";
    "2.5P 7/2";
    "2.5P 6/2";
    "2.5P 5/2";
    "2.5P 4/2";
    "2.5P 3/2";
    "2.5P 2.5/2";
    "2.5P 8/4";
    "2.5P 7/4";
    "2.5P 6/4";
    "2.5P 5/4";
    "2.5P 4/4";
    "2.5P 3/4";
    "2.5P 2.5/4";
    "2.5P 7/6";
    "2.5P 6/6";
    "2.5P 5/6";
    "2.5P 4/6";
    "2.5P 3/6";
    "2.5P 7/8";
    "2.5P 6/8";
    "2.5P 5/8";
    "2.5P 4/8";
    "2.5P 3/8";
    "2.5P 6/10";
    "2.5P 5/10";
    "2.5P 4/10";
    "5P 9/1";
    "5P 8/1";
    "5P 7/1";
    "5P 6/1";
    "5P 5/1";
    "5P 4/1";
    "5P 3/1";
    "5P 2.5/1";
    "5P 9/2";
    "5P 8/2";
    "5P 7/2";
    "5P 6/2";
    "5P 5/2";
    "5P 4/2";
    "5P 3/2";
    "5P 2.5/2";
    "5P 8/4";
    "5P 7/4";
    "5P 6/4";
    "5P 5/4";
    "5P 4/4";
    "5P 3/4";
    "5P 2.5/4";
    "5P 7/6";
    "5P 6/6";
    "5P 5/6";
    "5P 4/6";
    "5P 3/6";
    "5P 2.5/6";
    "5P 7/8";
    "5P 6/8";
    "5P 5/8";
    "5P 4/8";
    "5P 3/8";
    "5P 5/10";
    "5P 4/10";
    "7.5P 9/2";
    "7.5P 8/2";
    "7.5P 7/2";
    "7.5P 6/2";
    "7.5P 5/2";
    "7.5P 4/2";
    "7.5P 3/2";
    "7.5P 2.5/2";
    "7.5P 8/4";
    "7.5P 7/4";
    "7.5P 6/4";
    "7.5P 5/4";
    "7.5P 4/4";
    "7.5P 3/4";
    "7.5P 2.5/4";
    "7.5P 8/6";
    "7.5P 7/6";
    "7.5P 6/6";
    "7.5P 5/6";
    "7.5P 4/6";
    "7.5P 3/6";
    "7.5P 7/8";
    "7.5P 6/8";
    "7.5P 5/8";
    "7.5P 4/8";
    "7.5P 3/8";
    "7.5P 6/10";
    "7.5P 5/10";
    "7.5P 4/10";
    "10P 9/1";
    "10P 8/1";
    "10P 7/1";
    "10P 6/1";
    "10P 5/1";
    "10P 4/1";
    "10P 3/1";
    "10P 2.5/1";
    "10P 9/2";
    "10P 8/2";
    "10P 7/2";
    "10P 6/2";
    "10P 5/2";
    "10P 4/2";
    "10P 3/2";
    "10P 2.5/2";
    "10P 8/4";
    "10P 7/4";
    "10P 6/4";
    "10P 5/4";
    "10P 4/4";
    "10P 3/4";
    "10P 2.5/4";
    "10P 8/6";
    "10P 7/6";
    "10P 6/6";
    "10P 5/6";
    "10P 4/6";
    "10P 3/6";
    "10P 7/8";
    "10P 6/8";
    "10P 5/8";
    "10P 4/8";
    "10P 3/8";
    "10P 6/10";
    "10P 5/10";
    "10P 4/10";
    "10P 5/12";
    "2.5RP 9/2";
    "2.5RP 8/2";
    "2.5RP 7/2";
    "2.5RP 6/2";
    "2.5RP 5/2";
    "2.5RP 4/2";
    "2.5RP 3/2";
    "2.5RP 2.5/2";
    "2.5RP 8/4";
    "2.5RP 7/4";
    "2.5RP 6/4";
    "2.5RP 5/4";
    "2.5RP 4/4";
    "2.5RP 3/4";
    "2.5RP 2.5/4";
    "2.5RP 8/6";
    "2.5RP 7/6";
    "2.5RP 6/6";
    "2.5RP 5/6";
    "2.5RP 4/6";
    "2.5RP 3/6";
    "2.5RP 7/8";
    "2.5RP 6/8";
    "2.5RP 5/8";
    "2.5RP 4/8";
    "2.5RP 7/10";
    "2.5RP 6/10";
    "2.5RP 5/10";
    "2.5RP 4/10";
    "2.5RP 6/12";
    "2.5RP 5/12";
    "5RP 9/1";
    "5RP 8/1";
    "5RP 7/1";
    "5RP 6/1";
    "5RP 5/1";
    "5RP 4/1";
    "5RP 3/1";
    "5RP 2.5/1";
    "5RP 9/2";
    "5RP 8/2";
    "5RP 7/2";
    "5RP 6/2";
    "5RP 5/2";
    "5RP 4/2";
    "5RP 3/2";
    "5RP 2.5/2";
    "5RP 8/4";
    "5RP 7/4";
    "5RP 6/4";
    "5RP 5/4";
    "5RP 4/4";
    "5RP 3/4";
    "5RP 2.5/4";
    "5RP 8/6";
    "5RP 7/6";
    "5RP 6/6";
    "5RP 5/6";
    "5RP 4/6";
    "5RP 3/6";
    "5RP 7/8";
    "5RP 6/8";
    "5RP 5/8";
    "5RP 4/8";
    "5RP 6/10";
    "5RP 5/10";
    "5RP 4/10";
    "5RP 5/12";
    "5RP 4/12";
    "7.5RP 9/2";
    "7.5RP 8/2";
    "7.5RP 7/2";
    "7.5RP 6/2";
    "7.5RP 5/2";
    "7.5RP 4/2";
    "7.5RP 3/2";
    "7.5RP 2.5/2";
    "7.5RP 8/4";
    "7.5RP 7/4";
    "7.5RP 6/4";
    "7.5RP 5/4";
    "7.5RP 4/4";
    "7.5RP 3/4";
    "7.5RP 2.5/4";
    "7.5RP 8/6";
    "7.5RP 7/6";
    "7.5RP 6/6";
    "7.5RP 5/6";
    "7.5RP 4/6";
    "7.5RP 3/6";
    "7.5RP 7/8";
    "7.5RP 6/8";
    "7.5RP 5/8";
    "7.5RP 4/8";
    "7.5RP 6/10";
    "7.5RP 5/10";
    "7.5RP 4/10";
    "7.5RP 5/12";
    "7.5RP 4/12";
    "10RP 9/1";
    "10RP 8/1";
    "10RP 7/1";
    "10RP 6/1";
    "10RP 5/1";
    "10RP 4/1";
    "10RP 3/1";
    "10RP 2.5/1";
    "10RP 9/2";
    "10RP 8/2";
    "10RP 7/2";
    "10RP 6/2";
    "10RP 5/2";
    "10RP 4/2";
    "10RP 3/2";
    "10RP 2.5/2";
    "10RP 8/4";
    "10RP 7/4";
    "10RP 6/4";
    "10RP 5/4";
    "10RP 4/4";
    "10RP 3/4";
    "10RP 2.5/4";
    "10RP 8/6";
    "10RP 7/6";
    "10RP 6/6";
    "10RP 5/6";
    "10RP 4/6";
    "10RP 3/6";
    "10RP 7/8";
    "10RP 6/8";
    "10RP 5/8";
    "10RP 4/8";
    "10RP 6/10";
    "10RP 5/10";
    "10RP 4/10";
    "10RP 5/12";
    "10RP 4/12"]


text:
# # Lightning tour of MLJ

# *For a more elementary introduction to MLJ, see [Getting
# Started](https://alan-turing-institute.github.io/MLJ.jl/dev/getting_started/).*

# **Note.** Be sure this file has not been separated from the
# accompanying Project.toml and Manifest.toml files, which should not
# should be altered unless you know what you are doing. Using them,
# the following code block instantiates a julia environment with a tested
# bundle of packages known to work with the rest of the script:

using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

# Assuming Julia 1.6

# In MLJ a *model* is just a container for hyper-parameters, and that's
# all. Here we will apply several kinds of model composition before
# binding the resulting "meta-model" to data in a *machine* for
# evaluation, using cross-validation.

# Loading and instantiating a gradient tree-boosting model:

using MLJ
MLJ.color_off()

Booster = @load EvoTreeRegressor # loads code defining a model type
booster = Booster(max_depth=2)   # specify hyper-parameter at construction

#-

booster.nrounds=50               # or mutate post facto
booster

# This model is an example of an iterative model. As is stands, the
# number of iterations `nrounds` is fixed.

# ### Composition 1: Wrapping the model to make it "self-iterating"

# Let's create a new model that automatically learns the number of iterations,
# using the `NumberSinceBest(3)` criterion, as applied to an
# out-of-sample `l1` loss:

using MLJIteration
iterated_booster = IteratedModel(model=booster,
                                 resampling=Holdout(fraction_train=0.8),
                                 controls=[Step(2), NumberSinceBest(3), NumberLimit(300)],
                                 measure=l1,
                                 retrain=true)

# ### Composition 2: Preprocess the input features

# Combining the model with categorical feature encoding:

pipe = @pipeline ContinuousEncoder iterated_booster

# ### Composition 3: Wrapping the model to make it "self-tuning"

# First, we define a hyper-parameter range for optimization of a
# (nested) hyper-parameter:

max_depth_range = range(pipe,
                        :(deterministic_iterated_model.model.max_depth),
                        lower = 1,
                        upper = 10)

# Now we can wrap the pipeline model in an optimization strategy to make
# it "self-tuning":

self_tuning_pipe = TunedModel(model=pipe,
                              tuning=RandomSearch(),
                              ranges = max_depth_range,
                              resampling=CV(nfolds=3, rng=456),
                              measure=l1,
                              acceleration=CPUThreads(),
                              n=50)

# ### Binding to data and evaluating performance

# Loading a selection of features and labels from the Ames
# House Price dataset:

X, y = @load_reduced_ames;

# Binding the "self-tuning" pipeline model to data in a *machine* (which
# will additionally store *learned* parameters):

mach = machine(self_tuning_pipe, X, y)

# Evaluating the "self-tuning" pipeline model's performance using 5-fold
# cross-validation (implies multiple layers of nested resampling):

evaluate!(mach,
          measures=[l1, l2],
          resampling=CV(nfolds=5, rng=123),
          acceleration=CPUThreads())

using Literate #src
Literate.markdown(@__FILE__, @__DIR__, execute=false) #src
Literate.notebook(@__FILE__, @__DIR__, execute=true) #src

ext: jl


text:
function setref!(x::EXPR, binding::Binding)
    if !hasmeta(x)
        x.meta = Meta()
    end
    x.meta.ref = binding
    binding.refs !== nothing && push!(binding.refs, x)
end

function setref!(x::EXPR, binding)
    if !hasmeta(x)
        x.meta = Meta()
    end
    x.meta.ref = binding
end

# Main function to be called. Given the `state` tries to determine what `x`
# refers to. If it remains unresolved and is in a delayed evaluation scope
# (i.e. a function) it gets pushed to list (.urefs) to be resolved after we've
# run over the entire top-level scope.
function _resolve_ref(x, state)
    if !(parentof(x) isa EXPR && typof(parentof(x)) == CSTParser.Quotenode)
        resolved = resolve_ref(x, state.scope, state, [])
        if !resolved && (state.delayed || isglobal(valof(x), state.scope))
            push!(state.urefs, x)
        end
    end
end

# The first method that is tried. Searches the current scope for local bindings
# that match `x`. Steps:
# 1. Check whether we've already checked this scope (inifinite loops are
# possible when traversing nested modules.)
# 2. Check what sort of EXPR we're dealing with, separate name from EXPR that
# binds.
# 3. Look in the scope's variable list for a binding matching the name.
# 4. If 3. is unsuccessful, check whether the scope imports any modules then check them.
# 5. If no match is found within this scope check the parent scope.
# The return value is a boolean that is false if x should point to something but
# can't be resolved.

function resolve_ref(x::EXPR, scope::Scope, state::State, visited_scopes)::Bool
    hasref(x) && return true

    resolved = false
    if (typof(scope.expr) === CSTParser.ModuleH || typof(scope.expr) === CSTParser.BareModule) && CSTParser.length(scope.expr) > 1 && CSTParser.typof(scope.expr[2]) === IDENTIFIER
        s_m_name = scope.expr[2].val isa String ? scope.expr[2].val : ""
        if s_m_name in visited_scopes
            return resolved
        else
            push!(visited_scopes, s_m_name)
        end
    end

    if is_getfield(x)
        return resolve_getfield(x, scope, state)
    elseif isidentifier(x)
        if typof(x) === IDENTIFIER
            mn = valof(x)
            x1 = x
        else
            # NONSTDIDENTIFIER, e.g. var"name"
            mn = valof(x[2])
            x1 = x
        end
        if (mn == "__source__" || mn == "__module__") && _in_macro_def(x)
            setref!(x, Binding(noname, nothing, nothing, [], nothing, nothing))
            return true
        end
    elseif resolvable_macroname(x)
        x1 = x[2]
        mn = string("@", valof(x1))
    elseif typof(x) === x_Str
        if typof(x[1]) === IDENTIFIER
            x1 = x[1]
            mn = string("@", valof(x1), "_str")
        else
            return false
        end
    elseif typof(x) === CSTParser.Kw
        # Note to self: this seems wronge - Binding should be attached to entire Kw EXPR.
        if typof(x[1]) === IDENTIFIER
            setref!(x[1], Binding(noname, nothing, nothing, [], nothing, nothing))
        elseif typof(x[1]) === BinaryOpCall && kindof(x[1][2]) === CSTParser.Tokens.DECLARATION && typof(x[1][1]) === IDENTIFIER
            setref!(x[1][1], Binding(noname, nothing, nothing, [], nothing, nothing))
        end
        return true
    else
        return true
    end
    if scopehasbinding(scope, mn)
        setref!(x1, scope.names[mn])
        resolved = true
    elseif scope.modules isa Dict && length(scope.modules) > 0
        for m in scope.modules
            resolved = resolve_ref(x, m[2], state, visited_scopes)
            resolved && return true
        end
    end
    if !resolved && !scope.ismodule && parentof(scope) isa Scope
        return resolve_ref(x, parentof(scope), state, visited_scopes)
    end
    return resolved
end

# Searches a module store for a binding/variable that matches the reference `x1`.
function resolve_ref(x1::EXPR, m::SymbolServer.ModuleStore, state::State, visited_scopes)::Bool
    hasref(x1) && return true
    if isidentifier(x1)
        x = x1
        if Symbol(valof(x)) == m.name.name
            setref!(x, m)
            return true
        elseif isexportedby(x, m)
            val = m[Symbol(valof(x))]
            if val isa SymbolServer.VarRef
                val1 = SymbolServer._lookup(val, getsymbolserver(state.server))
                if val1 !== nothing
                    setref!(x, val1)
                    return true
                else
                    return false
                end
            else
                setref!(x, val)
                return true
            end
        end
    elseif typof(x1) === MacroName
        x = x1[2]
        mn = Symbol("@", valof(x))
        if isexportedby(mn, m)
            setref!(x, m[mn])
            return true
        end
    elseif typof(x1) === x_Str
        mac = x1
        mn = Symbol("@", valof(mac[1]), "_str")
        if isexportedby(mn, m)
            setref!(mac[1], m[mn])
            return true
        end
    end
    return false
end

# Fallback method
function resolve_ref(x::EXPR, m, state::State, visited_scopes)::Bool
    return hasref(x)::Bool
end

"""
    resolve_getfield(x::EXPR, parent::Union{EXPR,Scope,ModuleStore,Binding}, state::State)::Bool

Given an expression of the form `parent.x` try to resolve `x`. The method
called with `parent::EXPR` resolves the reference for `parent`, other methods
then check whether the Binding/Scope/ModuleStore to which `parent` points has
a field matching `x`.
"""
function resolve_getfield(x::EXPR, scope::Scope, state::State)::Bool
    hasref(x) && return true
    resolved = false
    if isidentifier(x[1])
        resolved = resolve_ref(x[1], scope, state, [])
        if resolved && typof(x[3]) === Quotenode && isidentifier(x[3][1])
            resolved = resolve_getfield(x[3][1], refof(x[1]), state)
        end
    elseif is_getfield_w_quotenode(x[1])
        resolved = resolve_ref(x[1], scope, state, [])
        if resolved && typof(x[3]) === Quotenode && length(x[3]) > 0 && isidentifier(x[3][1])
            resolved = resolve_getfield(x[3][1], refof(x[1][3][1]), state)
        end
    end
    return resolved
end

function resolve_getfield(x::EXPR, parent_type::EXPR, state::State)::Bool
    hasref(x) && return true
    resolved = false
    if CSTParser.isidentifier(x)
        if CSTParser.defines_module(parent_type) && scopeof(parent_type) isa Scope
            resolved = resolve_ref(x, scopeof(parent_type), state, [])
        elseif CSTParser.defines_struct(parent_type)
            if scopehasbinding(scopeof(parent_type), valof(x))
                setref!(x, scopeof(parent_type).names[valof(x)])
                resolved = true
            end
        end
    end
    return resolved
end

function resolve_getfield(x::EXPR, b::Binding, state::State)::Bool
    hasref(x) && return true
    resolved = false
    if b.type isa Binding
        resolved = resolve_getfield(x, b.type.val, state)
    elseif b.val isa SymbolServer.ModuleStore
        resolved = resolve_getfield(x, b.val, state)
    elseif b.type isa SymbolServer.DataTypeStore
        resolved = resolve_getfield(x, b.type, state)
    elseif b.val isa EXPR && CSTParser.defines_module(b.val)
        resolved = resolve_getfield(x, b.val, state)
    elseif b.val isa Binding && b.val.val isa EXPR && CSTParser.defines_module(b.val.val)
        resolved = resolve_getfield(x, b.val.val, state)
    end
    return resolved
end

function resolve_getfield(x::EXPR, parent_type, state::State)::Bool
    hasref(x) && return true
    return false
end

function resolve_getfield(x::EXPR, m::SymbolServer.ModuleStore, state::State)::Bool
    hasref(x) && return true
    resolved = false
    if CSTParser.isidentifier(x) && haskey(m, Symbol(valof(x)))
        val = m[Symbol(valof(x))]
        if val isa SymbolServer.VarRef
            val = SymbolServer._lookup(val, getsymbolserver(state.server))
            !(val isa SymbolServer.SymStore) && return false
        end
        setref!(x, val)
        resolved = true
    end
    return resolved
end

function resolve_getfield(x::EXPR, parent::SymbolServer.DataTypeStore, state::State)::Bool
    hasref(x) && return true
    resolved = false
    if CSTParser.isidentifier(x) && Symbol(valof(x)) in parent.fieldnames
        fi = findfirst(f->Symbol(valof(x)) == f, parent.fieldnames)
        ft = parent.types[fi]
        val = SymbolServer._lookup(ft, getsymbolserver(state.server))
        if val !== nothing
            setref!(x, Binding(noname, nothing, val, [], nothing, nothing))
            resolved = true
        end
    end
    return resolved
end

resolvable_macroname(x::EXPR) = typof(x) === MacroName && length(x) == 2 && isidentifier(x[2]) && refof(x[2]) === nothing

function _in_macro_def(x::EXPR)
    if typof(x) === CSTParser.Macro
        return true
    elseif parentof(x) isa EXPR
        return _in_macro_def(parentof(x))
    else
        return false
    end
end


text:
##############################################################################
##
## Combination behaves like [A B C ...] without forming it
##
##############################################################################

type Combination{N}
    A::NTuple{N, Matrix{Float64}}
    cumlength::Vector{Int}
end

function Combination(A::Matrix{Float64}...)
    cumlength = cumsum([size(x, 2) for x in A])
    Combination(A, cumlength)
end

function size(c::Combination, i)
    if i == 1
        size(c.A[1], 1)
    elseif i == 2
        c.cumlength[end]
    end
end

function view(c::Combination, ::Colon, j)
    index = searchsortedfirst(c.cumlength, j)
    newj = j
    if index > 1
        newj = j - c.cumlength[index-1]
    end
    view(c.A[index], :, newj)
end

##############################################################################
##
## Crossprod computes [A B C ...]' [A B C ...] without forming it
##
##############################################################################

# Construct [A B C]'[A B C] without generating [A B C]
function crossprod{N}(c::Combination{N})
    out = Array(Float64,  size(c, 2), size(c, 2))
    idx = 0
    for j in 1:size(c, 2)
        viewj = view(c, :, j)
        @inbounds for i in j:size(c, 2)
            idx += 1
            out[i, j] = dot(viewj, view(c, :, i))
        end
    end
    # make symmetric
    @inbounds for j in 1:size(c, 2), i in 1:(j-1)
        out[i, j] = out[j, i]
    end
    return out
end
crossprod(A::Matrix{Float64}) = A'A
crossprod(A::Matrix{Float64}...) = crossprod(Combination(A...))

##############################################################################
##
## Returns base of [A B C ...]
##
##############################################################################

# rank(A) == rank(A'A)
function basecol(X::Matrix{Float64}...)
    chol = cholfact!(crossprod(X...), :U, Val{true})
    ipermute!(diag(chol.factors) .> 0, chol.piv)
end

function getcols(X::Matrix{Float64},  basecolX::BitArray{1})
    if sum(basecolX) == size(X, 2)
        return X
    else
        return X[:, basecolX]
    end
end


text:
function extend_to_face(x,faceid,cell::IntervalBox{2})
    if faceid == 1
        return extend(x,2,cell[2].lo)
    elseif faceid == 2
        return extend(x,1,cell[1].hi)
    elseif faceid == 3
        return extend(x,2,cell[2].hi)
    elseif faceid == 4
        return extend(x,1,cell[1].lo)
    else
        throw(ArgumentError("Expected faceid ∈ {1,2,3,4}, got faceid = $faceid"))
    end
end

function extend_face_roots(roots,faceids,cell)
    extended_roots = hcat([extend_to_face(r,f,cell) for (r,f) in zip(roots,faceids)]...)
    return extended_roots
end

function roots_without_end(f,xL,xR,atol)
    r = sort!(find_zeros(f,xL,xR))
    if length(r) > 0
        isapprox(r[end],xR,atol=atol) ? pop!(r) : nothing
    end
    return r
end

function roots_without_start(f,xL,xR,atol)
    r = sort!(find_zeros(f,xL,xR))
    if length(r) > 0
        isapprox(r[1],xL,atol=atol) ? popfirst!(r) : nothing
    end
    return r
end

function roots_on_edges(funcs,xL,xR,yL,yR;atol=0.0)
    r1 = roots_without_end(funcs[1],xL,xR,atol)
    f1 = repeat([1],length(r1))

    r2 = roots_without_end(funcs[2],yL,yR,atol)
    f2 = repeat([2],length(r2))

    r3 = roots_without_start(funcs[3],xL,xR,atol)
    f3 = repeat([3],length(r3))

    r4 = roots_without_start(funcs[4],yL,yR,atol)
    f4 = repeat([4],length(r4))

    r = vcat(r1,r2,r3,r4)
    f = vcat(f1,f2,f3,f4)

    return r,f
end

function element_face_intersections(poly,cell::IntervalBox{2})

    xL,xR = (cell[1].lo,cell[1].hi)
    yL,yR = (cell[2].lo,cell[2].hi)

    funcs = restrict_on_faces(poly,cell)
    roots,faceids = roots_on_edges(funcs,xL,xR,yL,yR)

    @assert length(roots) == 2

    x1 = extend_to_face(roots[1],faceids[1],cell)
    x2 = extend_to_face(roots[2],faceids[2],cell)

    return x1,x2
end

function gradient_descent_to_zero!(x,P,∇P,atol,maxiter)
    @assert atol > 0.0
    @assert maxiter > 0

    val = P(x)
    iter = 1
    while abs(val) > atol && iter < maxiter
        dval = vec(∇P(x))
        alpha = val/dot(dval,dval)
        x .-= alpha*dval
        val = P(x)
        iter += 1
    end
    if iter == maxiter
        throw(ErrorException("Newton iteration failed to converge in $maxiter iterations"))
    end
end

function gradient_descent_to_zero!(x,P::InterpolatingPolynomial,atol,maxiter)
    gradient_descent_to_zero!(x,P,x->gradient(P,x),atol,maxiter)
end

function resolve_zero_levelset(poly,refpoints,xL,xR;atol=1e-12,maxiter=50)
    dim = dimension(poly)
    @assert dim == 2
    np = length(refpoints)

    points = zeros(2,np)
    dx = xR - xL

    for (idx,p) in enumerate(refpoints)
        x = xL + 0.5*(1.0+p)*dx
        gradient_descent_to_zero!(x,poly,atol,maxiter)
        points[:,idx] = x
    end
    return points
end

function fit_zero_levelset(poly,basis1d,quad1d,mass,cell;atol=1e-12,maxiter=50)
    xL,xR = element_face_intersections(poly,cell)
    intp = resolve_zero_levelset(poly,quad1d.points,xL,xR,atol=atol,maxiter=maxiter)
    rhs = linear_form(intp,basis1d,quad1d)
    return mass\rhs
end


text:
function find_library(ruby::AbstractString)
  v = rbexpr(ruby, "RUBY_VERSION", "")
  libs = [ dlprefix*"ruby"*v, dlprefix*"ruby" ]
  lib = rbconfig_expand(ruby, "\$(LIBRUBY_SO)")
  lib != nothing && unshift!(libs, splitext(lib)[1])
  libs = unique(libs)

  libpaths = [ rbconfig_expand(ruby, "\$(libdir)") ]
  ruby_path = rbexpr(ruby, "RbConfig.ruby")
  if ruby_path != nothing
    if is_windows()
      push!(libpaths, dirname(ruby_path))
    else
      push!(libpaths, joinpath(dirname(dirname(ruby_path)), "lib"))
    end
  end
  if is_apple()
    # TODO: The framework directory of the system ruby should be added in libpaths
  end

  exec_prefix = rbconfig(ruby, "exec_prefix")
  push!(libpaths, exec_prefix)
  push!(libpaths, joinpath(exec_prefix, "lib"))

  error_strings = Compat.String[]

  # find libruby
  for lib in libs
    for libpath in libpaths
      libpath_lib = joinpath(libpath, lib)
      if isfile(libpath_lib*"."*Libdl.dlext)
        try
          return (Libdl.dlopen(libpath_lib,
                               Libdl.RTLD_LAZY | Libdl.RTLD_DEEPBIND | Libdl.RTLD_GLOBAL),
                  libpath_lib)
        catch e
          push!(error_strings, string("dlopen($libpath_lib) ==> ", e))
        end
      end
    end
  end

  # find libruby from the system library path
  for lib in libs
    lib = splitext(lib)[1]
    try
      return (Libdl.dlopen(lib, Libdl.RTLD_LAZY | Libdl.RTLD_DEEPBIND | Libdl.RTLD_GLOBAL),
              lib)
    catch e
      push!(error_strings, string("dlopen($lib) ==> ", e))
    end
  end
end

const ruby = "ruby"

const (libruby_handle, libruby) = find_library(ruby)
const programname = rbexpr(ruby, "RbConfig.ruby")
const ruby_version = convert(VersionNumber, rbexpr(ruby, "RUBY_VERSION"))


text:
# Quate
function BioAFM_Q(afmData, pdbList, quateList, sigma)::BioAFMResult
    imgHeight, imgWidth = size(afmData)
    pdbNum = size(pdbList, 1)
    quateNum = size(quateList, 1)
    posteriors = Float64[]
    mostPropability = Float64(0)
    mostPropabilityPDB = Int(0)
    mostPropabilityQuate = Int(0)

    for pdbId in 1:pdbNum
        sumPropability = Float64(0)
        for quateId in 1:quateNum

            pdbRotated = MDToolbox.rotate(pdbList[pdbId], quateList[quateId, :])
            nowAfmData = MDToolbox.afmize(pdbRotated, (150.0, 150.0), (imgHeight, imgWidth))
            nowPropability = gaussian(nowAfmData, afmData, sigma)
            sumPropability += nowPropability
            if mostPropability < nowPropability
                mostPropability = nowPropability
                mostPropabilityPDB = pdbId
                mostPropabilityQuate = quateId
            end
        end
        push!(posteriors, sumPropability)
    end
    posteriors ./= sum(posteriors)
    return BioAFMResult(posteriors, mostPropabilityPDB, mostPropabilityQuate)
end

max_stars_repo_name: matsunagalab/BioAFM

max_issues_repo_name: matsunagalab/BioAFM

max_forks_repo_name: matsunagalab/BioAFM


text:
# type in Julia: Pkg.test("EmABCDQR")

using EmABCDQR
using Base.Test

include("../src/ABCDQR_simula_test.jl")
include("../src/ABCDQR_kfilter_test.jl")
include("../src/ABCDQR_em_test.jl")
include("../src/ABCDQR_em1_test.jl")


text:
# module moore2018multipropopt

path,_ = splitdir(@__FILE__); path = "$path/../"
# using PyPlot
#Packages I haven't developed
import Interpolations
# using PyPlot
import JLD
import CSV
import Dierckx

# Packages I played a very minor role
import PreComp
import Composites
import VTKtools

# Packages I played a minor role
import OptParams
import Atmosphere
import CCBlade #new CCBlade package w/o ND splined af data
import BeamFEA
import Snopt
import Gradients

# Packages I played a major role
import BrentMin
import MotorPower
import Akima
# import VLM
include("$path/oldvlm.jl")
import BPM
import AirfoilPrep

printiter = 37171
#TODO verify this works with parallel

include("$path/objective.jl")
include("$path/propstructmisc.jl")
include("$path/propwingmisc.jl")
include("$path/setup.jl")
include("$path/save_vtk.jl")

# # include("$path/../../BeamFEA.jl/src/BeamFEA.jl")
# rc("figure", figsize=(5.8, 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.17, bottom=0.18, top=0.97, right=.69)
# rc("axes", color_cycle=["348ABD", "A60628", "009E73", "7A68A6", "D55E00", "CC79A7"])
# plot_cycle=["#348ABD", "#A60628", "#009E73", "#7A68A6", "#D55E00", "#CC79A7"]
# close("all")

# function optimize(numprops,to_factor,db_factor)

    # define run controls
    ctlparams = controlparams(
    objective= :takeoff, #:range, #:mass,#:range #:ETA_total_cruise,
    printfreq = 70,
    savevtk = 70,
    printobjective=true,
    printdetailedoutput=true,
    printconstraintviolations=true,
    )

    designvars =  [
    # :p_orientation,
    :p_lam_t, #1-12
    :p_chord, #13-18
    :p_twist_d, #19-24
    :p_pitch_d_accel, #25
    :p_pitch_d_climb, #25
    :p_pitch_d_cruise, #25

    :velocity_accel,
    :velocity_climb,
    # :velocity_cruise,
    :p_Rtip,
    :p_rpm_accel,
    :p_rpm_climb,
    :p_rpm_cruise,
    # :number_blades,
    :kv, #28
    :i0, #29

    # :w_orientation,
    # :w_chord, #30-32
    :w_aoa_d_accel,   #multi
    :w_aoa_d_climb,   #multi
    :w_aoa_d_cruise,   #multi
    # :w_twist_d, #33-35
    # :w_halfspan,
    # :w_sweep_d,
    # :w_dihedral_d,
    # :w_lam_t, #36-end
    :batt_mass_accel,
    :batt_mass_climb,
    :batt_mass_cruise,
    ]

    # Define general constraints
    constraints = [
    :p_materialfailure,
    :p_localbuckling,
    :p_chordlamthick,
    # :w_materialfailure,
    # :w_localbuckling,
    # :w_chordlamthick,

    :totalmass,
    :machtip,
    :Lift,
    :Drag,
    :noise,

    :range,
    # :takeoff_dist,
    :cruisespeed,
    :radius,
    :flightpathangle,
    :stall,
    :capacity,
    :p_alphas,
    :p_separation,
    # :span,
    ]

    n_CP = 120.0
    payload = 684.0 + 245.0 #airframe + payload
    max_mass = 1230.0
    gravity = 9.81
    h_set = 15.24  #meters or 50 feet
    specif_energy = 300.0  #wh/kg
    required_range = 50.0*1000  #km
    altitude_accel = 0.0
    altitude_climb = 0.0
    altitude_cruise = 155.0
    cruise_factor = 100.0
    required_cruise = cruise_factor * 0.447 #MPH TO m/s

    # designvars, constraints, ctlparams = getoptparams()
    numprops = 8
    to_factor = 1
    db_factor = 1

    noise_factor = 76.0
    required_noise = noise_factor#68.0/db_factor
    required_takeoff_dist = 395.0/to_factor  #meters

    plyprops = plyproperties()

    matnames = plyprops.names

    name = "highmodulus_uni"
    idx = find(matnames -> matnames == name,matnames)
    uni_t = plyprops.plies[idx[1]].t

    name = "highmodulus_weave"
    idx = find(matnames -> matnames == name,matnames)
    weave_t = plyprops.plies[idx[1]].t

    name = "taylor_foam"
    idx = find(matnames -> matnames == name,matnames)
    foam_t = plyprops.plies[idx[1]].t

    # Prop inputs
    p_n_akima = 8
    Rtip = 0.5
    p_radii=[0.1,0.4,0.7,1.0]#[0.1,0.25,0.4,0.6,0.8,1.0]
    full_twist = ones(p_radii)*30.0
    full_chord = ones(p_radii)*0.07
    full_p_x_offset = ones(p_radii)*0.25
    full_p_aerocenter = ones(p_radii)*0.275

    p_lam_t = [
    ones(p_radii)*uni_t*10;#t1 highmodulus_uni
    ones(p_radii)*weave_t*85;#t2 highmodulus_weave
    ]
    p_usedmaterials = ["highmodulus_uni","highmodulus_weave"]
    p_orientation = [0.0,45.0]

    #Wing inputs
    # w_n_linear = 5
    # w_chord = [0.805,0.705,0.68]
    # w_halfspan = 10.5
    # w_nondim_halfspany = [0.0,.2,1.0]
    # full_w_x_offset = ones(w_nondim_halfspany)*0.25
    # full_w_aerocenter = ones(w_nondim_halfspany)*0.275
    # w_dihedral_d = zeros(w_chord)
    # w_sweep_d = ones(w_chord)*0.0

    #TODO: Technam inputs
    w_n_linear = 5
    w_halfspan = 11.4
    w_chord = [1.8/7,1.8/7,1.0/7]*w_halfspan
    w_nondim_halfspany = [0.0,4.0/7,1.0]
    full_w_x_offset = ones(w_nondim_halfspany)*0.25
    full_w_x_offset[end] = 1.0
    full_w_aerocenter = ones(w_nondim_halfspany)*0.275
    w_dihedral_d = zeros(w_chord)
    w_sweep_d = ones(w_chord)*0.0

    w_lam_t = [
    ones(w_nondim_halfspany)*uni_t*20;#t1 highmodulus_uni
    ones(w_nondim_halfspany)*weave_t*20;#t2 highmodulus_weave
    ones(w_nondim_halfspany)*foam_t*2;#t2 highmodulus_weave
    ones(w_nondim_halfspany)*weave_t*20;#t2 highmodulus_weave
    ones(w_nondim_halfspany)*foam_t*2;#t2 highmodulus_weave
    ]

    w_usedmaterials = ["highmodulus_uni","highmodulus_weave","taylor_foam","highmodulus_weave","taylor_foam"]
    w_orientation = [0.0,45.0,45.0] #uni, weave, webweave

    designparams = designparameters(;
    p_n_akima = p_n_akima,
    p_radii=p_radii,
    p_webloc=[0.999],
    altitude=0.0,
    plyprops=plyprops,
    p_usedmaterials = p_usedmaterials,
    p_orientation=p_orientation, #an p_orientation for each material at each section
    # structural design variables
    p_lam_t = p_lam_t,
    # geometric design variables
    p_chord = full_chord,
    p_sweep_d = fill(0.0,length(p_radii)-1),
    p_twist_d = full_twist,
    p_pitch_d = 1.0,
    p_dihedral_d = fill(0.0,length(p_radii)-1),
    p_airfoilthickness = fill(1.0,length(p_radii)), #not used
    p_aerocenter = full_p_aerocenter,
    p_x_offset = full_p_x_offset,
    # operating point design variables
    velocity = 25.0,
    # propulsion design variables
    p_Rtip = Rtip,
    blades = 2,
    numprops = numprops,
    machtip = 0.8,
    p_rpm = 1000.0,
    kv = 50.271513933940295,
    i0 = 2.0,
    batt_mass = 100.0,
    # Wing Design variables
    w_n_linear = w_n_linear, w_webloc=[0.5],
    w_lam_t = w_lam_t,
    w_chord = w_chord,
    w_aoa_d = 12.0,
    w_twist_d = ones(w_chord)*0.0,
    w_halfspan = w_halfspan,
    w_nondim_halfspany = w_nondim_halfspany,
    w_sweep_d = w_sweep_d,
    w_dihedral_d = w_dihedral_d,
    w_usedmaterials = w_usedmaterials,
    w_orientation = w_orientation,
    w_aerocenter = full_w_aerocenter,
    w_x_offset = full_w_x_offset,

    n_CP = n_CP, # 60.0,
    payload = payload, # 684.0,
    max_mass = max_mass, # 1500.0,
    gravity = gravity, # 9.81,
    h_set = h_set, # 20.0, #meters
    specif_energy = specif_energy, # 300.0, #wh/kg
    required_range = required_range, # 150000.0, #km
    required_takeoff_dist = required_takeoff_dist, # 395.0, #meters
    altitude_accel = altitude_accel, # 0.0,
    altitude_climb = altitude_climb, # 0.0,
    altitude_cruise = altitude_cruise, # 3000.0,
    required_noise = required_noise,
    required_cruise = required_cruise,

    )

    optxbounds = OptParams.initdict()
    getoptxbounds!(optxbounds,designparams)

    # get inputs to objective function
    x0,lb,ub = OptParams.assembleinput(designvars,optxbounds)
    rangevar = OptParams.getrangedict(designvars,optxbounds)

#             number_blades = 5.0
#
#
#             #Slow
#                 kv = 5.0
#                 i0 = 2.2037300540130143
#
#                 p_lam_t = [0.0007782, 0.000674115, 0.00050298, 0.000238061, 0.000218, 0.000218, 0.000218, 0.000218]
#                 p_chord = [0.217695, 0.179515, 0.0962492, 0.0129004]
#                 p_twist_d = [39.721, 18.0451, 10.7105, 4.20071]
#                 p_pitch_d_accel = 9.021458633187988
#                 p_pitch_d_climb = 2.6779561382353223
#                 p_pitch_d_cruise = 9.955675862211843
#                 velocity_accel = 21.907192616714354
#                 velocity_climb = 15.022172822158165
#                 velocity_cruise = 26.141124740636286
#                 p_Rtip = 0.1470747752685889
#                 p_rpm_accel = 4550.787583268113
#                 p_rpm_climb = 3617.25555613192
#                 p_rpm_cruise = 4206.390824864598
#                 w_aoa_d_accel = 11.805258303183681
#                 w_aoa_d_climb = 12.338223392213303
#                 w_aoa_d_cruise = 7.444036776377783
#                 batt_mass_accel = 0.01
#                 batt_mass_climb = 1.78807293443778
#                 batt_mass_cruise = 48.4944436895349
#
#                 number_blades = 2.0
#
#                 #SLOW
#                     kv = 5.0
#                     i0 = 5.999996165805586
#
#                 p_lam_t = [0.0106565, 0.00233204, 0.000521854, 0.000643721, 0.000235849, 0.000218, 0.000221644, 0.00023394]
#     p_chord = [0.152893, 0.189842, 0.0122381, 0.0126047]
#     p_twist_d = [43.6566, 25.3336, 14.8314, 2.15019]
#     p_pitch_d_accel = -1.8829344927066909
#     p_pitch_d_climb = -6.623043953201483
#     p_pitch_d_cruise = 3.1389505191969
#     velocity_accel = 10.516263781338967
#     velocity_climb = 22.110235464233497
#     velocity_cruise = 22.34991575727592
#     p_Rtip = 1.4607155806562915
#     p_rpm_accel = 1363.801616202054
#     p_rpm_climb = 1762.6245635737948
#     p_rpm_cruise = 785.088616955683
#     number_blades = 2.018354828584329
#     kv = 5.0000000099351745
#     i0 = 5.999999998436644
#     w_aoa_d_accel = 25.7671847059628
#     w_aoa_d_climb = 24.085518334858087
#     w_aoa_d_cruise = 7.6874985270624485
#     batt_mass_accel = 0.01
#     batt_mass_climb = 1.2565648978309747
#     batt_mass_cruise = 43.031131277642885
#
#     p_lam_t = [0.003252, 0.00142364, 0.000362052, 0.00099775, 0.000235849, 0.000218, 0.000221119, 0.000230127]
#     p_chord = [0.175241, 0.162084, 0.0107099, 0.0100064]
#     p_twist_d = [31.4471, 26.7358, 19.7654, 9.63553]
#     p_pitch_d_accel = -7.719866630174941
#     p_pitch_d_climb = -13.960686428934565
#     p_pitch_d_cruise = 1.173855411701484
#     velocity_accel = 10.178355866604997
#     velocity_climb = 10.271952373597921
#     p_Rtip = 1.3668277834692588
#     p_rpm_accel = 1145.630443055776
#     p_rpm_climb = 1665.6227960199028
#     p_rpm_cruise = 788.352329791602
#     number_blades = 2.445602851739818
#     kv = 6.105269032410999
#     i0 = 4.968452555691652
#     w_aoa_d_accel = 26.193250259194144
#     w_aoa_d_climb = 23.474670351125177
#     w_aoa_d_cruise = 7.79790809369224
#     batt_mass_accel = 0.01
#     batt_mass_climb = 0.6598943736785684
#     batt_mass_cruise = 99.85217871333208
#
#     p_lam_t = [0.00349585, 0.00154944, 0.000152, 0.00099775, 0.000235849, 0.000218, 0.000221119, 0.000230127]
# p_chord = [0.175239, 0.161664, 0.0110125, 0.0100079]
# p_twist_d = [32.2629, 27.3062, 21.6001, 10.2174]
# p_pitch_d_accel = -7.719866630174941
# p_pitch_d_climb = -13.924415137156194
# p_pitch_d_cruise = 1.173855411701484
# velocity_accel = 10.178355866604997
# velocity_climb = 10.197874857065022
# p_Rtip = 1.3668274103862978
# p_rpm_accel = 1145.630443055776
# p_rpm_climb = 1674.2859783769227
# p_rpm_cruise = 788.352329791602
# number_blades = 2.1438870567106965
# kv = 6.105269032410999
# i0 = 4.968452555691652
# w_aoa_d_accel = 26.145740407471234
# w_aoa_d_climb = 23.522365831369434
# w_aoa_d_cruise = 7.797908093691838
# batt_mass_accel = 0.01
# batt_mass_climb = 0.43827470796137974
# batt_mass_cruise = 98.31020185664771
#
# number_blades = 2.0
#
# p_lam_t = [0.00296032, 0.0015862, 0.000162432, 0.000559344, 0.000218, 0.000218, 0.000218, 0.000218]
# p_chord = [0.152775, 0.209865, 0.0599548, 0.0146912]
# p_twist_d = [43.7258, 36.3345, 24.8117, 15.5926]
# p_pitch_d_accel = 1.4774821755019851
# p_pitch_d_climb = 0.07394619362544558
# p_pitch_d_cruise = 1.2375953673499482
# velocity_accel = 19.472605774436627
# velocity_climb = 21.727538432606607
# velocity_cruise = 44.93261088033517
# p_Rtip = 0.7373054880998322
# p_rpm_accel = 1521.2024594724965
# p_rpm_climb = 2662.377310410668
# p_rpm_cruise = 2030.4950665640044
# number_blades = 2.0#2.3818948496117045
# kv = 9.01975436914855
# i0 = 5.942697466865178
# w_aoa_d_accel = 26.783254994297305
# w_aoa_d_climb = 15.682348122387333
# w_aoa_d_cruise = 1.9046025617909106
# batt_mass_accel = 0.01
# batt_mass_climb = 0.44404948094245933
# batt_mass_cruise = 119.39074611391598
#
# p_lam_t = [0.00491215, 0.00207895, 0.00109735, 0.000951505, 0.000218, 0.000218, 0.000218, 0.000218]
# p_chord = [0.151192, 0.325948, 0.192366, 0.0520607]
# p_twist_d = [45.0, 42.1689, 28.5539, 24.6288]
# p_pitch_d_accel = 0.5160959129483517
# p_pitch_d_climb = -6.379369415208218
# p_pitch_d_cruise = 1.2554747468011076
# velocity_accel = 17.504799071392068
# velocity_climb = 16.011054271640106
# velocity_cruise = 45.556861878613994
# p_Rtip = 0.6333333333333337
# p_rpm_accel = 1700.7399651777887
# p_rpm_climb = 2214.9896316661266
# p_rpm_cruise = 2033.7975016015014
# kv = 9.010532347513744
# i0 = 5.983991910449766
# w_aoa_d_accel = 26.4773214375485
# w_aoa_d_climb = 24.947976808858428
# w_aoa_d_cruise = 6.09880815849964
# batt_mass_accel = 0.0012680910284242287
# batt_mass_climb = 0.2174067264581435
# batt_mass_cruise = 122.52852597671871
#
# #With own af data
# p_lam_t = [0.00893571, 0.00234997, 0.00182375, 0.000152, 0.00083697, 0.000953682, 0.000876195, 0.000540092]
# p_chord = [0.106509, 0.150221, 0.0692545, 0.05]
# p_twist_d = [45.0, 39.3542, 30.3407, 28.5409]
# p_pitch_d_accel = -14.801534303269628
# p_pitch_d_climb = -16.663949276805756
# p_pitch_d_cruise = 1.116079551638493
# velocity_accel = 9.944097279374834
# velocity_climb = 7.73005569315936
# p_rpm_accel = 3785.4038793504687
# p_rpm_climb = 4097.5298786080875
# p_rpm_cruise = 1942.562723001659
# kv = 5.0
# i0 = 3.9951737452421248
# w_aoa_d_accel = 32.04188683229715
# w_aoa_d_climb = 31.337251317158273
# w_aoa_d_cruise = 4.691196096841078
# batt_mass_accel = 0.001
# batt_mass_climb = 1.1260003355278403
# batt_mass_cruise = 74.94512715737756
#
# p_lam_t = [0.00888127, 0.00239683, 0.00121549, 0.000152, 0.00083697, 0.000953682, 0.000876195, 0.000540092]
# p_chord = [0.116602, 0.143469, 0.0767974, 0.0262864]
# p_twist_d = [45.0, 38.8659, 30.5658, 29.6107]
# p_pitch_d_accel = -15.245276832060593
# p_pitch_d_climb = -16.990543523778836
# p_pitch_d_cruise = 1.1216558284591813
# velocity_accel = 9.951767417608675
# velocity_climb = 7.696453282205797
# p_rpm_accel = 3771.664591019246
# p_rpm_climb = 4097.529878608089
# p_rpm_cruise = 1943.4767749266518
# kv = 5.000001441672897
# i0 = 3.995010265667068
# w_aoa_d_accel = 32.07085586412711
# w_aoa_d_climb = 31.354963988641288
# w_aoa_d_cruise = 4.698885062498646
# batt_mass_accel = 0.001
# batt_mass_climb = 1.009896688608764
# batt_mass_cruise = 75.72755761609312
#
# p_lam_t = [0.00888127, 0.00239683, 0.00121549, 0.000152, 0.00083697, 0.000953682, 0.000876195, 0.000540092]
# p_chord = [0.116602, 0.143469, 0.0767974, 0.0262864]
# p_twist_d = [45.0, 38.8659, 30.5658, 29.6107]
# p_pitch_d_accel = -15.245276832060593
# p_pitch_d_climb = -16.990543523778836
# p_pitch_d_cruise = 1.1216558284591813
# velocity_accel = 9.951767417608675
# velocity_climb = 7.696453282205797
# p_rpm_accel = 3771.664591019246
# p_rpm_climb = 4097.529878608089
# p_rpm_cruise = 1943.4767749266518
# kv = 5.000001441672897
# i0 = 3.995010265667068
# w_aoa_d_accel = 32.07085586412711
# w_aoa_d_climb = 31.354963988641288
# w_aoa_d_cruise = 4.698885062498646
# batt_mass_accel = 0.001
# batt_mass_climb = 1.009896688608764
# batt_mass_cruise = 75.72755761609312
#
# #fixed aoa
# p_lam_t = [0.0090107, 0.00229688, 0.00114354, 0.000152023, 0.00083697, 0.000953682, 0.000218, 0.000218]
# p_chord = [0.102455, 0.150463, 0.0727664, 0.0242614]
# p_twist_d = [45.0, 39.5646, 30.6596, 28.9851]
# p_pitch_d_accel = -15.291890229352216
# p_pitch_d_climb = -16.82857992276188
# p_pitch_d_cruise = 1.1252752346571662
# velocity_accel = 7.719281971361207
# velocity_climb = 11.864801002989012
# p_rpm_accel = 3804.787843784844
# p_rpm_climb = 4097.529878608089
# p_rpm_cruise = 1943.5718130812047
# kv = 5.0
# i0 = 3.9791093851566774
# w_aoa_d_accel = 31.23708813422462
# w_aoa_d_climb = 32.507899772141734
# w_aoa_d_cruise = 4.649161517426887
# batt_mass_accel = 0.001
# batt_mass_climb = 0.7805441172811979
# batt_mass_cruise = 74.80296024119511
#
# #Fixed Structures
# p_lam_t = [0.007538, 0.00173269, 0.00088642, 0.000152, 0.00083697, 0.000953682, 0.000218, 0.000218]
# p_chord = [0.151763, 0.151763, 0.151763, 0.151763]
# p_twist_d = [45.0, 39.1948, 30.4775, 29.6842]
# p_pitch_d_accel = -15.296567262681904
# p_pitch_d_climb = -16.56620424726087
# p_pitch_d_cruise = 1.1254471622631674
# velocity_accel = 8.248197902866446
# velocity_climb = 12.52998048101345
# p_rpm_accel = 3811.247887669918
# p_rpm_climb = 4097.529878608089
# p_rpm_cruise = 1943.6393645230926
# kv = 5.008349580516639
# i0 = 3.966580647765332
# w_aoa_d_accel = 30.770162029306185
# w_aoa_d_climb = 32.18321792176881
# w_aoa_d_cruise = 4.655881487309745
# batt_mass_accel = 0.001
# batt_mass_climb = 0.8400014116756471
# batt_mass_cruise = 74.58121110817127
#
# w_halfspan = 5.7
#
# #Example GIF
# p_lam_t = [0.007538, 0.00173269, 0.00088642, 0.000152, 0.00083697, 0.000953682, 0.000218, 0.000218]
# p_chord = [0.151763, 0.151763, 0.151763, 0.151763]
# p_twist_d = [29.0, 29.0, 29.0, 29.0]
# p_pitch_d_accel = -15.296567262681904
# p_pitch_d_climb = -16.56620424726087
# p_pitch_d_cruise = 1.1254471622631674
# velocity_accel = 8.248197902866446
# velocity_climb = 12.52998048101345
# p_rpm_accel = 3811.247887669918
# p_rpm_climb = 4097.529878608089
# p_rpm_cruise = 1943.6393645230926
# kv = 5.008349580516639
# i0 = 3.966580647765332
# w_aoa_d_accel = 30.770162029306185
# w_aoa_d_climb = 32.18321792176881
# w_aoa_d_cruise = 4.655881487309745
# batt_mass_accel = 0.001
# batt_mass_climb = 0.8400014116756471
# batt_mass_cruise = 74.58121110817127
#
# w_halfspan = 5.7
#
# #Fixed Structures
# p_lam_t = [0.007538, 0.00173269, 0.00088642, 0.000152, 0.00083697, 0.000953682, 0.000218, 0.000218]*1.5
# p_chord = [0.101763, 0.153117, 0.0724272, 0.0237822]*1.5
# p_twist_d = [45.0, 39.1948, 30.4775, 29.6842]*0.9
# p_pitch_d_accel = -15.296567262681904*1.1
# p_pitch_d_climb = -16.56620424726087*1.1
# p_pitch_d_cruise = 1.1254471622631674*1.1
# velocity_accel = 8.248197902866446*1.1
# velocity_climb = 12.52998048101345*1.1
# p_rpm_accel = 3811.247887669918*1.1
# p_rpm_climb = 4097.529878608089*1.1
# p_rpm_cruise = 1943.6393645230926*1.1
# kv = 5.008349580516639
# i0 = 3.966580647765332
# w_aoa_d_accel = 30.770162029306185*1.1
# w_aoa_d_climb = 32.18321792176881*1.1
# w_aoa_d_cruise = 4.655881487309745*1.1
# batt_mass_accel = 0.001*1.1
# batt_mass_climb = 0.8400014116756471*1.1
# batt_mass_cruise = 74.58121110817127*1.1

w_halfspan = 5.7
p_chord = ones(4)*.2
p_twist_d = ones(4)*30.0
p_pitch_d_accel = -5.0
p_pitch_d_climb = -3.0
p_pitch_d_cruise = 0.0
velocity_accel = 13.0
velocity_climb = 21.0
velocity_cruise = 45.556861878613994
p_rpm_accel = 10000.0
p_rpm_climb = 10000.0
p_rpm_cruise = 5000.0
kv = 10.0
i0 = 5.0
w_aoa_d_accel = 15.0
w_aoa_d_climb = 30.0
w_aoa_d_cruise = 5.0
batt_mass_accel = 1.0
batt_mass_climb = 1.0
batt_mass_cruise = 20.0
p_Rtip = 0.6
number_blades = 2.0

p_lam_t = [0.0103505, 0.000285035, 0.000734952, 0.000348269, 0.000388793, 0.0017972, 0.000717596, 0.000407097]
p_chord = [0.131096, 0.151933, 0.0718815, 0.0189429]
p_twist_d = [44.8067, 38.1564, 30.961, 29.4853]
p_pitch_d_accel = -15.281648370110014
p_pitch_d_climb = -16.16372411483507
p_pitch_d_cruise = 1.3095227963163993
velocity_accel = 8.148039950692041
velocity_climb = 14.409762122371184
p_Rtip = 0.5814151924837737
p_rpm_accel = 4383.71038341194
p_rpm_climb = 4256.8755486191085
p_rpm_cruise = 2225.53942029128
kv = 5.000183261626615
i0 = 3.3656590406139397
w_aoa_d_accel = 31.226885842207395
w_aoa_d_climb = 28.68626884140169
w_aoa_d_cruise = 4.829112321652306
batt_mass_accel = 0.001
batt_mass_climb = 2.1438224245890556
batt_mass_cruise = 88.06496124285786

p_lam_t2 = [0.00317757, 0.000785375, 0.000970208, 0.000323186, 0.00101117, 0.00187457, 0.000479939, 0.000345597]
p_chord2 = [0.127766, 0.142388, 0.0683798, 0.0267748]
p_twist_d2 = [44.966, 39.5187, 30.8776, 28.7751]
p_pitch_d_accel2 = -16.298004244335992
p_pitch_d_climb2 = -16.26869739354685
p_pitch_d_cruise2 = 1.2379918784894843
p_Rtip2 = 0.6333333333333337
velocity_accel2 = 7.85487134395491
velocity_climb2 = 12.545434643740208
p_rpm_accel2 = 4095.136682148108
p_rpm_climb2 = 4097.501219539165
p_rpm_cruise2 = 1935.8454972489915
kv2 = 5.000390427715033
i02 = 3.9838408230865894
w_aoa_d_accel2 = 30.519865300480163
w_aoa_d_climb2 = 32.156697969032514
w_aoa_d_cruise2 = 4.655850131710671
batt_mass_accel2 = 0.0010000086985410138
batt_mass_climb2 = 0.8418729122777846
batt_mass_cruise2 = 74.2385048375922

blend = 0.2

p_lam_t = p_lam_t*blend + p_lam_t2*(1-blend)
p_chord = p_chord*blend + p_chord2*(1-blend)
p_twist_d = p_twist_d*blend + p_twist_d2*(1-blend)
p_pitch_d_accel = p_pitch_d_accel*blend + p_pitch_d_accel2*(1-blend)
p_pitch_d_climb = p_pitch_d_climb*blend + p_pitch_d_climb2*(1-blend)
p_pitch_d_cruise = p_pitch_d_cruise*blend + p_pitch_d_cruise2*(1-blend)
p_Rtip = p_Rtip*blend + p_Rtip2*(1-blend)
velocity_accel = velocity_accel*blend + velocity_accel2*(1-blend)
velocity_climb = velocity_climb*blend + velocity_climb2*(1-blend)
p_rpm_accel = p_rpm_accel*blend + p_rpm_accel2*(1-blend)
p_rpm_climb = p_rpm_climb*blend + p_rpm_climb2*(1-blend)
p_rpm_cruise = p_rpm_cruise*blend + p_rpm_cruise2*(1-blend)
kv = kv*blend + kv2*(1-blend)
i0 = i0*blend + i02*(1-blend)
w_aoa_d_accel = w_aoa_d_accel*blend + w_aoa_d_accel2*(1-blend)
w_aoa_d_climb = w_aoa_d_climb*blend + w_aoa_d_climb2*(1-blend)
w_aoa_d_cruise = w_aoa_d_cruise*blend + w_aoa_d_cruise2*(1-blend)
batt_mass_accel = batt_mass_accel*blend + batt_mass_accel2*(1-blend)
batt_mass_climb = batt_mass_climb*blend + batt_mass_climb2*(1-blend)
batt_mass_cruise = batt_mass_cruise*blend + batt_mass_cruise2*(1-blend)

    OptParams.setvar!(:p_lam_t,x0,optxbounds,rangevar,p_lam_t)
    OptParams.setvar!(:p_chord,x0,optxbounds,rangevar,p_chord)
    OptParams.setvar!(:p_twist_d,x0,optxbounds,rangevar,p_twist_d)
    OptParams.setvar!(:p_pitch_d_accel,x0,optxbounds,rangevar,p_pitch_d_accel)
    OptParams.setvar!(:p_pitch_d_climb,x0,optxbounds,rangevar,p_pitch_d_climb)
    OptParams.setvar!(:p_pitch_d_cruise,x0,optxbounds,rangevar,p_pitch_d_cruise)
    OptParams.setvar!(:velocity_accel,x0,optxbounds,rangevar,velocity_accel)
    OptParams.setvar!(:velocity_climb,x0,optxbounds,rangevar,velocity_climb)
    OptParams.setvar!(:velocity_cruise,x0,optxbounds,rangevar,velocity_cruise)
    OptParams.setvar!(:p_Rtip,x0,optxbounds,rangevar,p_Rtip)
    OptParams.setvar!(:p_rpm_accel,x0,optxbounds,rangevar,p_rpm_accel)
    OptParams.setvar!(:p_rpm_climb,x0,optxbounds,rangevar,p_rpm_climb)
    OptParams.setvar!(:p_rpm_cruise,x0,optxbounds,rangevar,p_rpm_cruise)
    OptParams.setvar!(:kv,x0,optxbounds,rangevar,kv)
    OptParams.setvar!(:i0,x0,optxbounds,rangevar,i0)
    OptParams.setvar!(:w_aoa_d_accel,x0,optxbounds,rangevar,w_aoa_d_accel)
    OptParams.setvar!(:w_aoa_d_climb,x0,optxbounds,rangevar,w_aoa_d_climb)
    OptParams.setvar!(:w_aoa_d_cruise,x0,optxbounds,rangevar,w_aoa_d_cruise)
    OptParams.setvar!(:batt_mass_accel,x0,optxbounds,rangevar,batt_mass_accel)
    OptParams.setvar!(:batt_mass_climb,x0,optxbounds,rangevar,batt_mass_climb)
    OptParams.setvar!(:batt_mass_cruise,x0,optxbounds,rangevar,batt_mass_cruise)
    OptParams.setvar!(:number_blades,x0,optxbounds,rangevar,number_blades)

    OptParams.setvar!(:w_halfspan,x0,optxbounds,rangevar,w_halfspan)

    # sparce vars 900kg and ~50m/s vinf all
    # xwarm = [0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 3881.27, 2112.42, 434.783, 3999.49, 4499.16, 2783.33, 152.197, 15.5839, -1814.98, 526.502, 11.2334, 5151.27, 2163.88, 1201.38, 1201.38, 4685.41, 3446.64, 20000.0, 1.48825, 1.85602, 1.0, 2.97122, 3.99704, 14.1136, 168.285, 169.996, 160.905, 114.286, 0.0440764, 0.0809586, 0.0139054, 0.01, 0.01, 0.01, 0.312975, 1.00904, 0.207683, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 63.476]

    # x0[1:length(xwarm)] = xwarm[1:end]

    table_af_wing_blown = JLD.load("$path/airfoils/af_prop_E212.jld")
    table_3D_wing_blown = table_af_wing_blown["NDtable"]
    spline_3D_blown = AirfoilPrep.SplineND_from_tableND(table_3D_wing_blown)

    table_af_wing_free = JLD.load("$path/airfoils/af_prop_E212.jld")
    table_3D_wing_free = table_af_wing_free["NDtable"]
    spline_3D_freestream = AirfoilPrep.SplineND_from_tableND(table_3D_wing_free)

    # # Load the prop airfoil ND table
    # ClarkY_non_extrap = JLD.load("$(path)/airfoils/af_prop_E212.jld")
    # ClarkY_non_extrap = ClarkY_non_extrap["NDtable"]
    # AR = 1.0
    # omega = p_rpm_climb*pi/30
    # TSR = omega * p_Rtip / velocity_cruise
    #
    # grid_alphas=[i for i in -180:1.0:180]
    # Rtip = OptParams.getvar(:p_Rtip,x0,optxbounds,rangevar)
    # chord_prop_pt = OptParams.getvar(:p_chord,x0,optxbounds,rangevar)
    # rpt = linspace(.2,.911,length(chord_prop_pt)) * Rtip
    # r = linspace(.2,.911,p_n_akima) * Rtip
    # chord_prop = Akima.interp_noderiv(rpt, chord_prop_pt, r)
    # chord_prop[chord_prop.<1E-3] = 1E-3
    # af = Array{Any}(p_n_akima)
    # for i = 1:p_n_akima
    #
    #    r_over_R = r[i]
    #    c_over_r = chord_prop[i]/r[i]
    #    NDextrap3D_3Dtable = AirfoilPrep.NDTable_correction3D_extrap(ClarkY_non_extrap,r_over_R,c_over_r,TSR;grid_alphas=grid_alphas)
    #    afspl = AirfoilPrep.SplineND_from_tableND(NDextrap3D_3Dtable)
    #    af[i] = afspl
    # end
    # JLD.save("$(path)/airfoils/p$(numprops)_cr$(round(Int,required_cruise/0.447))_db$(round(Int,required_noise))af_rerun.jld", "af", af)
    af = JLD.load("$(path)/airfoils/p$(numprops)_cr$(round(Int,required_cruise/0.447))_db$(round(Int,required_noise))af_rerun.jld")
    af = af["af"]
    af = fill(af[1],p_n_akima)

    af_xy = CSV.read("$(path)/airfoils/e212-il.csv",header = true, delim = ',',datarow = 10,allowmissing=:none)
    xaf1 = reverse(af_xy[:,1])/100.0
    yaf1 = reverse(af_xy[:,2])/100.0

    # nameout = "$(path)/AirfoilData/e212.dat"
    # write(nameout,"x\ty\n")
    # open(nameout,"a") do x
    #     for i = 1:length(xaf1)
    #         write(x,"$(xaf1[i])\t$(yaf1[i])\n")
    #     end
    # end

    #Filter airfoil to get points where strain is evaluated
    n_strainpts = 10
    xafstrain = zeros(n_strainpts)
    yafstrain = zeros(n_strainpts)
    step_simplify = round(Int,length(xaf1)/n_strainpts)
    idx = 1
    for i = 1:n_strainpts
        xafstrain[i] = xaf1[idx[1]]
        yafstrain[i] = yaf1[idx[1]]
        idx += step_simplify
    end

    # Fill airfoil data for each section
    p_xaf = zeros(p_n_akima,length(xaf1))
    p_yaf = zeros(p_n_akima,length(yaf1))
    p_xafstrain = zeros(p_n_akima,length(xafstrain))
    p_yafstrain = zeros(p_n_akima,length(yafstrain))
    for i = 1:p_n_akima
        p_xaf[i,:] = copy(xaf1) #break links since using in multiple areas
        p_yaf[i,:] = copy(yaf1)
        p_xafstrain[i,:] = copy(xafstrain)
        p_yafstrain[i,:] = copy(yafstrain)
    end

    #TODO: different airfoil for wing
    w_xaf = zeros(w_n_linear,length(xaf1))
    w_yaf = zeros(w_n_linear,length(yaf1))
    w_xafstrain = zeros(w_n_linear,length(xafstrain))
    w_yafstrain = zeros(w_n_linear,length(yafstrain))
    for i = 1:w_n_linear
        w_xaf[i,:] = copy(xaf1)
        w_yaf[i,:] = copy(yaf1)
        w_xafstrain[i,:] = copy(xafstrain)
        w_yafstrain[i,:] = copy(yafstrain)
    end

    func(x) = objective(x,optxbounds,rangevar,designparams,designvars,constraints,
    ctlparams,af,p_xaf,p_yaf,p_xafstrain,p_yafstrain,w_xaf,w_yaf,w_xafstrain,
    w_yafstrain,spline_3D_freestream,spline_3D_blown;detailed=detailed,
    savefinalvtk=savefinalvtk, verification=verification, plots=plots)

    function fungrad(x)
        function func1(x)
            j,C,_ = func(x)
            return [j;C]
        end
        f, dfdx = Gradients.forwarddiff(func1,x)
        lenx = length(x)

        return f[1],f[2:end],dfdx[1,:],dfdx[2:end,:],false
    end

    function pfungrad(x)
        function func1(x)
            j,C,_ = func(x)
            return [j;C]
        end
        f, dfdx = Gradients.centraldiff(func1,x;h=1E-5,parallel = true)

        return f[1],f[2:end],dfdx[1,:],dfdx[2:end,:],false
    end

    # ----- Define Optimizer Options ----- #
    options = Dict{String, Any}()
    options["Derivative level"] = 0
    options["Verify level"] = 0
    # options["Function precision"] = 1.00E-4
    # options["Difference interval"] = 1e-4
    # options["Central difference interval"] = 1e-4
    options["Iterations limit"] = 1e8
    options["Major iterations limit"] = 1000
    # options["Minor iterations limit"]= 1e8
    options["Major optimality tolerance"] = 1e-8 #Should be scaled so it is optimal with a solid 2 digits
    # options["Minor optimality  tolerance"] = 1e-6
    options["Major feasibility tolerance"] = 5e-6
    # options["Minor feasibility tolerance"] = 1e-6
    # options["Minor print level"] = 1E8
    # options["Print frequency"] = 1E8
    options["Scale option"] = 1
    # options["Scale tolerance"] = .95
    # options["Verify level"] = 3 #only if specifying gradients

    printfile = "./SNOPT_out/SNOPT_printfile_p$(numprops)_cr$(round(Int,required_cruise/0.447))_db$(round(Int,required_noise)).out"
    sumfile = "./SNOPT_out/SNOPT_sumfile_p$(numprops)_cr$(round(Int,required_cruise/0.447))_db$(round(Int,required_noise)).out"
    time0 = 0.0
    #if true
        verification=false
        savefinalvtk=false
        detailed = true
        plots = false

        detailed_output_accel,detailed_output_climb,detailed_output_cruise,
        design_accel,design_climb,design_cruise, detailed_system, f, c = objective(x0,optxbounds,rangevar,
        designparams,designvars,constraints,ctlparams,af,p_xaf,p_yaf,p_xafstrain,
        p_yafstrain,w_xaf,w_yaf,w_xafstrain,w_yafstrain,spline_3D_freestream,spline_3D_blown;verification=verification,
        savefinalvtk=savefinalvtk,
        detailed =detailed,
        plots = plots)

        # tic()
        # finaloutput,p_geom,p_struct, f, c = objective(x0,optxbounds,rangevar,
        # designparams,designvars,constraints,ctlparams,af,p_xaf,p_yaf,p_xafstrain,
        # p_yafstrain,w_xaf,w_yaf,w_xafstrain,w_yafstrain,spline_3D_freestream,spline_3D_blown;verification=verification,
        # savefinalvtk=savefinalvtk,
        # detailed =detailed,
        # plots = plots)
        # time0 = toc()

        JLD.save("$(path)/OptSummaries/p$(numprops)_cr$(round(Int,required_cruise/0.447))_db$(round(Int,required_noise))correct_span.jld",
        "detailed_output_accel", detailed_output_accel,
        "detailed_output_climb", detailed_output_climb,
        "detailed_output_cruise", detailed_output_cruise,
        "design_accel", design_accel,
        "design_climb", design_climb,
        "design_cruise", design_cruise,
        "detailed_system", detailed_system)
    #end

    dfdx1 = 0.0
    time1 = 0.0
    nameout = "./SNOPT_out/printout.txt"
    if false#true
        verification=false
        savefinalvtk=false
        detailed = false
        plots = false

        f1,c1,dfdx1,dcdx1,state = fungrad(x0)
        tic()
        f1,c1,dfdx1,dcdx1,state = fungrad(x0)
        time1 = toc()
    end

    if false#true
        verification=false
        savefinalvtk=false
        detailed = false
        plots = false

        # f2,c2,dfdx2,dcdx2,state = pfungrad(x0)
        tic()
        f2,c2,dfdx2,dcdx2,state = pfungrad(x0)
        time2 = toc()
        println("x0: $x0")
        println("f2: $(f2)")
        println("c2: $(c2)")
        println("dfdx2: $(dfdx2)")
        for i = 1:length(c2)
            if maximum(abs.(dcdx2[i,:]))>30.0 # || mean(abs.(dcdx2[i,:]))<1E-4
                println("dc$i:
                $(dcdx2[i,:])")
            end
        end
        println(dfdx1-dfdx2)

        open(nameout,"a")do x
            write(x,"x0: $x0
            f2: $f2
            dfdx2: $(dfdx2)
            $(dfdx1-dfdx2)
            nproc: $(nprocs())
            time0: $time0
            time1: $time1
            time2: $time2
            ")
        end

    end

# x0 = [0.207583, 0.16042, 0.0980674, 0.0341654, 0.0347087, 0.0203864, 0.0538978, 0.0121262, 3.15971, 2.49569, 1.83168, 0.66743, 3.73101, 3.05527, 2.26264, 1.56159, 406.152, -7.77481, -113.17, 50.0, 0.13628, 29.4291, 20.0, 88.1265, 8814.21, 2.40003, 26.5281, 2.87881, 0.4, 28.236, 3.47529, 0.001, 0.001, 0.5]

    optimize = false
    if optimize

        verification=false
        savefinalvtk=false
        detailed = false
        plots = false
        if nprocs()==1
            optfunc = func
        else
            optfunc = pfungrad
        end

        xopt, fopt, exitmsg = Snopt.snopt(optfunc, x0, lb, ub, options,
        printfile = printfile, sumfile = sumfile)

        println(exitmsg)
        println("fopt: $fopt")
        println("xopt: $xopt")

        verification=false
        savefinalvtk=true
        detailed = true
        plots = false

        detailed_output_accel,detailed_output_climb,detailed_output_cruise,
        design_accel,design_climb,design_cruise, detailed_system, f, c = objective(xopt,optxbounds,rangevar,
        designparams,designvars,constraints,ctlparams,af,p_xaf,p_yaf,p_xafstrain,
        p_yafstrain,w_xaf,w_yaf,w_xafstrain,w_yafstrain,spline_3D_freestream,spline_3D_blown;detailed=detailed, savefinalvtk=savefinalvtk, verification=verification, plots=plots)

        JLD.save("$(path)/OptSummaries/p$(numprops)_cr$(round(Int,required_cruise/0.447))_db$(round(Int,required_noise))correct_span.jld",
        "detailed_output_accel", detailed_output_accel,
        "detailed_output_climb", detailed_output_climb,
        "detailed_output_cruise", detailed_output_cruise,
        "design_accel", design_accel,
        "design_climb", design_climb,
        "design_cruise", design_cruise,
        "detailed_system", detailed_system)

        println("----------Detailed Output----------")
        for field in fieldnames(detailed_system)
            # if field in keys(rangevar)
            if length(string(field)) < 7
                println(field,":\t\t\t",getfield(detailed_system,field))
            elseif length(string(field)) < 15
                println(field,":\t\t",getfield(detailed_system,field))
            else
                println(field,":\t",getfield(detailed_system,field))
            end
            # end
        end
    end

    C = []
    J = []
    Jin = []
    Cin = []
    variable = []
    # lb = x0*0.97
    # ub = x0*1.03

    if false #sweep # sweep through the variables

        verification=false
        savefinalvtk=false
        detailed = false
        plots = false

        N_points = 20 #discretization for sweeping
        for i = 1:length(x0)#[37] #[1,13,19,25,26,28,29,30,33,37,38] #[1,13,19,25,26,27,28,29,30,31,34,37,38]

            xsweep = copy(x0)
            J1,C1,fail1 = func(xsweep)
            variable = linspace(x0[i]*0.99,x0[i]*1.01,N_points)
            J = zeros(N_points,length(J1))
            C = zeros(N_points,length(C1))
            for j = 1:N_points
                xsweep[i] = variable[j]
                Jin,Cin,_ = func(xsweep)
                J[j,:] = Jin
                C[j,:] = Cin
                xsweep = copy(x0)
            end

            # figure("objective")
            # clf()
            # figure("contraints")
            # clf()

            figure("objective")
            PyPlot.plot(variable/maximum(abs.(variable)),J,"-",label = "J variable $i") #/maximum(abs(J))
            legend()

            # PyPlot.figure("contraints")
            # for k = 1:length(C[1,:])
            #
            #     # if !isapprox(norm(C[:,k]), norm(C[:,k+1]); atol = norm(C[:,k])*1E-1)
            #     # if (maximum(C[:,k]/maximum(abs.(C[:,k])))<=1.0 && minimum(C[:,k]/maximum(abs.(C[:,k])))>-.9) #&& minimum((C[:,k]./maximum(abs.(C[:,k]))))<-0.75
            #         PyPlot.plot(variable,C[:,k]/maximum(abs.(C[:,k])),label = "C$k var $i")
            #     # end
            #     # end
            # end
            # PyPlot.legend()
            pause(.100001)
        end
    end

    # end #optimize

    # end #module


text:
# Julia wrapper for header: sunnonlinsol_fixedpoint.h
# Automatically generated using Clang.jl

function SUNNonlinSol_FixedPoint(y::N_Vector, m::Cint)
    ccall((:SUNNonlinSol_FixedPoint, libsundials_sunnonlinsolfixedpoint), SUNNonlinearSolver, (N_Vector, Cint), y, m)
end

function SUNNonlinSol_FixedPoint(y, m)
    __y = convert(NVector, y)
    SUNNonlinSol_FixedPoint(convert(N_Vector, __y), convert(Cint, m))
end

function SUNNonlinSol_FixedPointSens(count::Cint, y::N_Vector, m::Cint)
    ccall((:SUNNonlinSol_FixedPointSens, libsundials_sunnonlinsolfixedpoint), SUNNonlinearSolver, (Cint, N_Vector, Cint), count, y, m)
end

function SUNNonlinSol_FixedPointSens(count, y, m)
    __y = convert(NVector, y)
    SUNNonlinSol_FixedPointSens(convert(Cint, count), convert(N_Vector, __y), convert(Cint, m))
end

function SUNNonlinSolGetType_FixedPoint(NLS::SUNNonlinearSolver)
    ccall((:SUNNonlinSolGetType_FixedPoint, libsundials_sunnonlinsolfixedpoint), SUNNonlinearSolver_Type, (SUNNonlinearSolver,), NLS)
end

function SUNNonlinSolInitialize_FixedPoint(NLS::SUNNonlinearSolver)
    ccall((:SUNNonlinSolInitialize_FixedPoint, libsundials_sunnonlinsolfixedpoint), Cint, (SUNNonlinearSolver,), NLS)
end

function SUNNonlinSolSolve_FixedPoint(NLS::SUNNonlinearSolver, y0::N_Vector, y::N_Vector, w::N_Vector, tol::realtype, callSetup::Cint, mem)
    ccall((:SUNNonlinSolSolve_FixedPoint, libsundials_sunnonlinsolfixedpoint), Cint, (SUNNonlinearSolver, N_Vector, N_Vector, N_Vector, realtype, Cint, Ptr{Cvoid}), NLS, y0, y, w, tol, callSetup, mem)
end

function SUNNonlinSolSolve_FixedPoint(NLS::SUNNonlinearSolver, y0, y, w, tol, callSetup, mem)
    __y0 = convert(NVector, y0)
    __y = convert(NVector, y)
    __w = convert(NVector, w)
    SUNNonlinSolSolve_FixedPoint(NLS, convert(N_Vector, __y0), convert(N_Vector, __y), convert(N_Vector, __w), tol, convert(Cint, callSetup), mem)
end

function SUNNonlinSolFree_FixedPoint(NLS::SUNNonlinearSolver)
    ccall((:SUNNonlinSolFree_FixedPoint, libsundials_sunnonlinsolfixedpoint), Cint, (SUNNonlinearSolver,), NLS)
end

function SUNNonlinSolSetSysFn_FixedPoint(NLS::SUNNonlinearSolver, SysFn::SUNNonlinSolSysFn)
    ccall((:SUNNonlinSolSetSysFn_FixedPoint, libsundials_sunnonlinsolfixedpoint), Cint, (SUNNonlinearSolver, SUNNonlinSolSysFn), NLS, SysFn)
end

function SUNNonlinSolSetConvTestFn_FixedPoint(NLS::SUNNonlinearSolver, CTestFn::SUNNonlinSolConvTestFn, ctest_data)
    ccall((:SUNNonlinSolSetConvTestFn_FixedPoint, libsundials_sunnonlinsolfixedpoint), Cint, (SUNNonlinearSolver, SUNNonlinSolConvTestFn, Ptr{Cvoid}), NLS, CTestFn, ctest_data)
end

function SUNNonlinSolSetMaxIters_FixedPoint(NLS::SUNNonlinearSolver, maxiters::Cint)
    ccall((:SUNNonlinSolSetMaxIters_FixedPoint, libsundials_sunnonlinsolfixedpoint), Cint, (SUNNonlinearSolver, Cint), NLS, maxiters)
end

function SUNNonlinSolSetMaxIters_FixedPoint(NLS, maxiters)
    SUNNonlinSolSetMaxIters_FixedPoint(NLS, convert(Cint, maxiters))
end

function SUNNonlinSolSetDamping_FixedPoint(NLS::SUNNonlinearSolver, beta::realtype)
    ccall((:SUNNonlinSolSetDamping_FixedPoint, libsundials_sunnonlinsolfixedpoint), Cint, (SUNNonlinearSolver, realtype), NLS, beta)
end

function SUNNonlinSolGetNumIters_FixedPoint(NLS::SUNNonlinearSolver, niters)
    ccall((:SUNNonlinSolGetNumIters_FixedPoint, libsundials_sunnonlinsolfixedpoint), Cint, (SUNNonlinearSolver, Ptr{Clong}), NLS, niters)
end

function SUNNonlinSolGetCurIter_FixedPoint(NLS::SUNNonlinearSolver, iter)
    ccall((:SUNNonlinSolGetCurIter_FixedPoint, libsundials_sunnonlinsolfixedpoint), Cint, (SUNNonlinearSolver, Ptr{Cint}), NLS, iter)
end

function SUNNonlinSolGetNumConvFails_FixedPoint(NLS::SUNNonlinearSolver, nconvfails)
    ccall((:SUNNonlinSolGetNumConvFails_FixedPoint, libsundials_sunnonlinsolfixedpoint), Cint, (SUNNonlinearSolver, Ptr{Clong}), NLS, nconvfails)
end

function SUNNonlinSolGetSysFn_FixedPoint(NLS::SUNNonlinearSolver, SysFn)
    ccall((:SUNNonlinSolGetSysFn_FixedPoint, libsundials_sunnonlinsolfixedpoint), Cint, (SUNNonlinearSolver, Ptr{SUNNonlinSolSysFn}), NLS, SysFn)
end
# Julia wrapper for header: sunnonlinsol_newton.h
# Automatically generated using Clang.jl

function SUNNonlinSol_Newton(y::N_Vector)
    ccall((:SUNNonlinSol_Newton, libsundials_sunnonlinsolnewton), SUNNonlinearSolver, (N_Vector,), y)
end

function SUNNonlinSol_Newton(y)
    __y = convert(NVector, y)
    SUNNonlinSol_Newton(convert(N_Vector, __y))
end

function SUNNonlinSol_NewtonSens(count::Cint, y::N_Vector)
    ccall((:SUNNonlinSol_NewtonSens, libsundials_sunnonlinsolnewton), SUNNonlinearSolver, (Cint, N_Vector), count, y)
end

function SUNNonlinSol_NewtonSens(count, y)
    __y = convert(NVector, y)
    SUNNonlinSol_NewtonSens(convert(Cint, count), convert(N_Vector, __y))
end

function SUNNonlinSolGetType_Newton(NLS::SUNNonlinearSolver)
    ccall((:SUNNonlinSolGetType_Newton, libsundials_sunnonlinsolnewton), SUNNonlinearSolver_Type, (SUNNonlinearSolver,), NLS)
end

function SUNNonlinSolInitialize_Newton(NLS::SUNNonlinearSolver)
    ccall((:SUNNonlinSolInitialize_Newton, libsundials_sunnonlinsolnewton), Cint, (SUNNonlinearSolver,), NLS)
end

function SUNNonlinSolSolve_Newton(NLS::SUNNonlinearSolver, y0::N_Vector, y::N_Vector, w::N_Vector, tol::realtype, callLSetup::Cint, mem)
    ccall((:SUNNonlinSolSolve_Newton, libsundials_sunnonlinsolnewton), Cint, (SUNNonlinearSolver, N_Vector, N_Vector, N_Vector, realtype, Cint, Ptr{Cvoid}), NLS, y0, y, w, tol, callLSetup, mem)
end

function SUNNonlinSolSolve_Newton(NLS, y0, y, w, tol, callLSetup, mem)
    __y0 = convert(NVector, y0)
    __y = convert(NVector, y)
    __w = convert(NVector, w)
    SUNNonlinSolSolve_Newton(NLS, convert(N_Vector, __y0), convert(N_Vector, __y), convert(N_Vector, __w), tol, convert(Cint, callLSetup), mem)
end

function SUNNonlinSolFree_Newton(NLS::SUNNonlinearSolver)
    ccall((:SUNNonlinSolFree_Newton, libsundials_sunnonlinsolnewton), Cint, (SUNNonlinearSolver,), NLS)
end

function SUNNonlinSolSetSysFn_Newton(NLS::SUNNonlinearSolver, SysFn::SUNNonlinSolSysFn)
    ccall((:SUNNonlinSolSetSysFn_Newton, libsundials_sunnonlinsolnewton), Cint, (SUNNonlinearSolver, SUNNonlinSolSysFn), NLS, SysFn)
end

function SUNNonlinSolSetLSetupFn_Newton(NLS::SUNNonlinearSolver, LSetupFn::SUNNonlinSolLSetupFn)
    ccall((:SUNNonlinSolSetLSetupFn_Newton, libsundials_sunnonlinsolnewton), Cint, (SUNNonlinearSolver, SUNNonlinSolLSetupFn), NLS, LSetupFn)
end

function SUNNonlinSolSetLSolveFn_Newton(NLS::SUNNonlinearSolver, LSolveFn::SUNNonlinSolLSolveFn)
    ccall((:SUNNonlinSolSetLSolveFn_Newton, libsundials_sunnonlinsolnewton), Cint, (SUNNonlinearSolver, SUNNonlinSolLSolveFn), NLS, LSolveFn)
end

function SUNNonlinSolSetConvTestFn_Newton(NLS::SUNNonlinearSolver, CTestFn::SUNNonlinSolConvTestFn, ctest_data)
    ccall((:SUNNonlinSolSetConvTestFn_Newton, libsundials_sunnonlinsolnewton), Cint, (SUNNonlinearSolver, SUNNonlinSolConvTestFn, Ptr{Cvoid}), NLS, CTestFn, ctest_data)
end

function SUNNonlinSolSetMaxIters_Newton(NLS::SUNNonlinearSolver, maxiters::Cint)
    ccall((:SUNNonlinSolSetMaxIters_Newton, libsundials_sunnonlinsolnewton), Cint, (SUNNonlinearSolver, Cint), NLS, maxiters)
end

function SUNNonlinSolSetMaxIters_Newton(NLS, maxiters)
    SUNNonlinSolSetMaxIters_Newton(NLS, convert(Cint, maxiters))
end

function SUNNonlinSolGetNumIters_Newton(NLS::SUNNonlinearSolver, niters)
    ccall((:SUNNonlinSolGetNumIters_Newton, libsundials_sunnonlinsolnewton), Cint, (SUNNonlinearSolver, Ptr{Clong}), NLS, niters)
end

function SUNNonlinSolGetCurIter_Newton(NLS::SUNNonlinearSolver, iter)
    ccall((:SUNNonlinSolGetCurIter_Newton, libsundials_sunnonlinsolnewton), Cint, (SUNNonlinearSolver, Ptr{Cint}), NLS, iter)
end

function SUNNonlinSolGetNumConvFails_Newton(NLS::SUNNonlinearSolver, nconvfails)
    ccall((:SUNNonlinSolGetNumConvFails_Newton, libsundials_sunnonlinsolnewton), Cint, (SUNNonlinearSolver, Ptr{Clong}), NLS, nconvfails)
end

function SUNNonlinSolGetSysFn_Newton(NLS::SUNNonlinearSolver, SysFn)
    ccall((:SUNNonlinSolGetSysFn_Newton, libsundials_sunnonlinsolnewton), Cint, (SUNNonlinearSolver, Ptr{SUNNonlinSolSysFn}), NLS, SysFn)
end


text:
export Polynomial

"""
    Polynomial{T<:Number}(coeffs::AbstractVector{T}, var=:x)

Construct a polynomial from its coefficients `a`, lowest order first, optionally in
terms of the given variable `x`. `x` can be a character, symbol, or string.

If ``p = a_n x^n + \\ldots + a_2 x^2 + a_1 x + a_0``, we construct this through
`Polynomial([a_0, a_1, ..., a_n])`.

The usual arithmetic operators are overloaded to work with polynomials as well as
with combinations of polynomials and scalars. However, operations involving two
polynomials of different variables causes an error except those involving a constant polynomial.

# Examples
```@meta
DocTestSetup = quote
    using Polynomials
end
```

```jldoctest
julia> Polynomial([1, 0, 3, 4])
Polynomial(1 + 3*x^2 + 4*x^3)

julia> Polynomial([1, 2, 3], :s)
Polynomial(1 + 2*s + 3*s^2)

julia> one(Polynomial)
Polynomial(1.0)
```
"""
struct Polynomial{T <: Number} <: StandardBasisPolynomial{T}
    coeffs::Vector{T}
    var::Symbol
    function Polynomial{T}(coeffs::Vector{T}, var::Symbol) where {T <: Number}
        length(coeffs) == 0 && return new{T}(zeros(T, 1), var)
        last_nz = findlast(!iszero, coeffs)
        last = max(1, last_nz === nothing ? 0 : last_nz)
        return new{T}(coeffs[1:last], var)
    end
end

@register Polynomial

"""
    (p::Polynomial)(x)

Evaluate the polynomial using [Horner's Method](https://en.wikipedia.org/wiki/Horner%27s_method), also known as synthetic division, as implemented in `evalpoly` of base `Julia`.

```@meta
DocTestSetup = quote
    using Polynomials
end
```

```@meta
DocTestSetup = quote
    using Polynomials
end
```

# Examples
```jldoctest
julia> p = Polynomial([1, 0, 3])
Polynomial(1 + 3*x^2)

julia> p(0)
1

julia> p.(0:3)
4-element Array{Int64,1}:
  1
  4
 13
 28
```
"""
(p::Polynomial{T})(x::S) where {T,S} = evalpoly(x, coeffs(p))

function Base.:+(p1::Polynomial{T}, p2::Polynomial{S}) where {T, S}
    n1, n2 = length(p1), length(p2)
    R = promote_type(T,S)
    if n1 > 1 && n2 > 1
       p1.var != p2.var && error("Polynomials must have same variable")
       if n1 >= n2
          c = R.(copy(p1.coeffs))
          for i = 1:n2
            c[i] += p2.coeffs[i]
          end
        else
            c = R.(copy(p2.coeffs))
            for i = 1:n1
              c[i] += p1.coeffs[i]
            end
        end
        return Polynomial(c, p1.var)
    elseif n1 <= 1
       c = R.(copy(p2.coeffs))
       c[1] += p1[0]
       return Polynomial(c, p2.var)
    else
       c = R.(copy(p1.coeffs))
       c[1] += p2[0]
       return Polynomial(c, p1.var)
    end
end

function Base.:*(p1::Polynomial{T}, p2::Polynomial{S}) where {T,S}

    n, m = length(p1)-1, length(p2)-1 # not degree, so pNULL works
    if n > 0 && m > 0
        p1.var != p2.var && error("Polynomials must have same variable")
        R = promote_type(T, S)
        c = zeros(R, m + n + 1)
        for i in 0:n, j in 0:m
            @inbounds c[i + j + 1] += p1[i] * p2[j]
        end
        return Polynomial(c, p1.var)
    elseif n <= 0
        return Polynomial(copy(p2.coeffs) * p1[0], p2.var)
    else
        return Polynomial(copy(p1.coeffs) * p2[0], p1.var)
    end

end