text:
"""
    CoverParams(; [sel_tax=0.0], [set_shape=1.0],
                  [min_weight=1E-2], [mask_discount=0.9], [setXset_factor=1.0],
                  [uncovered_factor=0.1], [covered_factor=0.001],
                  [set_relevance_shape=0.5], [set_relevance_min=0.5]) -> CoverParams

Specify parameters for the [`AbstractCoverProblem`](@ref).
One can specify a non-default value for a particular parameter.

See [cover quality score](@ref cover_quality) for the detailed description of the parameters.

# Arguments
  * `sel_tax`: ``\\epsilon``, the constant added to the set score of each selected set to penalize
     the number of sets in the cover
  * `set_shape`: ``\\beta``, applied to set or set×set scores
  * `min_weight`: minimal non-zero set probability ???
  * `mask_discount`: ``\\alpha``, how much the overlap score of each subsequent mask (from most to less enriched)
    is discounted
  * `setXset_factor`: ``w_r``, the weight of *redundancy* [cover quality](@ref cover_quality) component
    (`setXset_score` scale), 0 = no redunancy penalty
  * `uncovered_factor`: ``w_u``, the weight of *uncovered hits* component in [cover quality](@ref cover_quality)
  * `covered_factor`: ``w_c``, the weight of *covered non-hits* component in [cover quality](@ref cover_quality)
  * `set_relevance_shape`: ``\\beta_L``, how much set relevance affects set score, 0 = no effect
  * `set_relevance_min`: ``L_{\\min``, if shaped relevance is below, it's set to `set_relevance_min`
"""
struct CoverParams
    sel_tax::Float64
    set_shape::Float64
    min_weight::Float64
    mask_discount::Float64
    setXset_factor::Float64
    uncovered_factor::Float64
    covered_factor::Float64
    set_relevance_shape::Float64
    set_relevance_min::Float64

    function CoverParams(;
                         sel_tax::Real=0.0, set_shape::Real=1.0,
                         min_weight::Real=1E-2, mask_discount::Real=0.9,
                         setXset_factor::Real=1.0,
                         uncovered_factor::Real=0.1, covered_factor::Real=0.001,
                         set_relevance_shape::Real=0.5,
                         set_relevance_min::Real=0.5)
        (0.0 < min_weight <= 1.0) || throw(ArgumentError("`min_weight` must be within (0,1] range"))
        (0.0 <= mask_discount <= 1.0) || throw(ArgumentError("`mask_discount` must be within [0,1] range"))
        (0.0 <= set_shape) || throw(ArgumentError("`set_shape` must be ≥0"))
        (0.0 <= setXset_factor) || throw(ArgumentError("`setXset_factor` must be ≥0"))
        (0.0 <= uncovered_factor) || throw(ArgumentError("`uncovered_factor` must be ≥0"))
        (0.0 <= covered_factor) || throw(ArgumentError("`covered_factor` must be ≥0"))
        (0.0 <= set_relevance_shape) || throw(ArgumentError("`set_relevance_shape` must be ≥0"))
        (0.0 <= set_relevance_min <= 1) || throw(ArgumentError("`set_relevance_min` must be within [0, 1] range"))
        new(sel_tax, set_shape, min_weight, mask_discount,
            setXset_factor, uncovered_factor, covered_factor,
            set_relevance_shape, set_relevance_min)
    end
end

"""
Linear component of an individual masked set score for the `CoverProblem`.
Doesn't take into account the overlap with the other selected sets.
"""
function overlap_score(masked::Number, set::Number, total_masked::Number, total::Number,
                       relevance::Number, params::CoverParams)
    # FIXME is it just tail=:both for one set of parameters
    #= P-value for masked-vs-set overlap enriched =# res = -(-logpvalue(masked, set, total_masked, total))^params.set_shape*max(relevance^params.set_relevance_shape, params.set_relevance_min) #-
    #= P-value for unmasked-vs-set overlap enriched =# #logpvalue(set - masked, set, total - total_masked, total)
    @assert !isnan(res) "masked=$masked set=$set total_masked=$total_masked total=$total relevance=$relevance res=NaN"
    return res
end

function _setweight(mosaic::AbstractWeightedSetMosaic,
                    setix::Integer, expix::Integer,
                    relevance::Number, params::CoverParams;
                    filter::Bool = true)
    w = _setweight(mosaic, setix, expix, filter=filter)
    ismissing(w) && return w
    res = -(-w)^params.set_shape *
            max(relevance^params.set_relevance_shape, params.set_relevance_min) #-
    #= P-value for unmasked-vs-set overlap enriched =# #logpvalue(set - masked, set, total - total_masked, total)
    @assert !isnan(res) "masked=$masked set=$set total_masked=$total_masked total=$total relevance=$relevance res=NaN"
    return res
end

var2set(mosaic::AbstractWeightedSetMosaic) = eachindex(mosaic.loc2glob_setix)

# calculates var scores
# var score is:
# the sum of overlap scores with all masked sets
# minus the log(var selection probability)
# plus nets*log(nsets), where nsets arr all masked sets overlapping with var
function var_scores(mosaic::AbstractWeightedSetMosaic, params::CoverParams,
                    var2setix::AbstractVector{Int} = var2set(mosaic))
    # find relevant tiles and
    # sum masked elements per tile in all masks and unmasked elements per tile in union mask
    # note: in masks there could be element not covered by any set/tile, so these counts refer to coverable elements
    # calculate the sum of scores of given set in each mask
    v_scores = Vector{Float64}(undef, length(var2setix))
    scores = Vector{Float64}() # temp array containing var overlap scores for each mask it overlaps

    @inbounds for (varix, setix) in enumerate(var2setix)
        empty!(scores)
        glob_setix = mosaic.loc2glob_setix[setix]
        setrelevance = originalmosaic(mosaic).set_relevances[glob_setix]
        for expix in 1:nexperiments(mosaic)
            w = _setweight(mosaic, setix, expix, setrelevance, params)
            ismissing(w) || push!(scores, w)
        end
        sort!(scores)
        scoresum = 0.0
        w = 1.0
        @inbounds for score in scores
            scoresum += score * w
            w *= params.mask_discount
        end
        v_scores[varix] = scoresum + params.sel_tax
    end
    return v_scores
end

# prepares
# 1) tile-to-var mapping (SparseMaskMatrix, only relevant tiles are considered)
# 2) var-to-(tile-in-mask) mapping (SparseMatrixCSC, only relevant tiles are considered)
# 3) the total number of masked elements in all masks for each tile
# 4) the total number of unmasked elements in all active masks (overlapping with the tile) for each tile
function varXtiles(mosaic::MaskedSetMosaic,
                   var2setix::AbstractVector{Int},
                   params::CoverParams)
    # find relevant tiles and
    # sum masked elements per tile in all masks and unmasked elements per tile in union mask
    # note: in masks there could be element not covered by any set/tile, so these counts refer to coverable elements
    tileixs = Vector{Int}()
    tile2nmasked = Vector{Int}()
    tile2nunmasked = Vector{Int}()
    vXmt_els = Vector{Tuple{Int, Int, Int, Float64}}()
    # max rank of the setXexp overlap that is still relevant
    setXexp_olaps = fill(0.0, nexperiments(mosaic))

    @inbounds for (varix, setix) in enumerate(var2setix)
        globsetix = mosaic.loc2glob_setix[setix]
        for expix in 1:nexperiments(mosaic)
            setXexp_olaps[expix] = coalesce(_setweight(mosaic, setix, expix), Inf)
        end
        # rank 1 = most significant (weights are negative)
        setXexp_olap_ranks = denserank(setXexp_olaps)

        set_tiles = view(mosaic.original.tileXset, :, globsetix)
        for tileix in set_tiles
            tile_elms = view(mosaic.original.elmXtile, :, tileix)
            tilepos = searchsortedfirst(tileixs, tileix)
            if tilepos > length(tileixs) || tileixs[tilepos] != tileix # do this once per tile
                insert!(tileixs, tilepos, tileix)
                insert!(tile2nunmasked, tilepos, length(tile_elms) - sum(view(mosaic.elunionmask, tile_elms)))
                nmasked = 0
                for (expix, olap) in enumerate(view(mosaic.setXexp_olaps, setix, :))
                    if !ismissing(olap)
                        nmasked += sum(view(mosaic.elmasks, tile_elms, expix))
                    end
                end
                insert!(tile2nmasked, tilepos, nmasked)
            end
            for (expix, olap) in enumerate(view(mosaic.setXexp_olaps, setix, :))
                ismissing(olap) && continue
                nmaskxtile = sum(view(mosaic.elmasks, tile_elms, expix))
                olap_rank = setXexp_olap_ranks[expix]
                olap_w = params.mask_discount ^ (olap_rank - 1)
                if nmaskxtile < length(tile_elms) && # there are unmasked tile elements
                   olap_w >= 0.1                     # the mask has relevant overlap with the var
                    push!(vXmt_els, (varix, expix, tileix, olap_w * (length(tile_elms) - nmaskxtile)))
                end
            end
        end
    end

    tileXvar = mosaic.original.tileXset[tileixs, mosaic.loc2glob_setix[var2setix]]
    if !isempty(tileXvar)
        # try to optimize tileXvar: find tiles that refer to the identical set of vars and merge their stats
        varXtile = transpose(tileXvar)
        vars2tiles = Vector{Tuple{Vector{Int}, Vector{Int}}}()
        for tileix in axes(varXtile, 2)
            varixs = varXtile[:, tileix]
            v2tpos = searchsortedfirst(vars2tiles, (varixs, 0), by=first)
            if v2tpos > length(vars2tiles) || first(vars2tiles[v2tpos]) != varixs
                insert!(vars2tiles, v2tpos, (varixs, [tileix]))
            else
                push!(vars2tiles[v2tpos][2], tileix)
            end
        end
        if length(vars2tiles) < size(tileXvar, 1) # there are tiles to merge
            # take the 1st (it's also minimal since tileixs are sorted) of each merged tiles
            tileXvar = tileXvar[[first(tileixs) for (_, tileixs) in vars2tiles], :]
            # aggregate n[un]masked over duplicated tiles
            tile2nmasked = [sum(i -> tile2nmasked[i], tileixs) for (_, tileixs) in vars2tiles]
            tile2nunmasked = [sum(i -> tile2nunmasked[i], tileixs) for (_, tileixs) in vars2tiles]
        end
    end
    if true # optimize varXmaskxtile: find maskXtile that refers to the identical set of vars and sum their vals
        maskxtile2vars = Dict{Tuple{Int, Int}, Vector{Int}}()
        for (varix, expix, tileix, val) in vXmt_els
            varixs = get!(() -> Vector{Int}(), maskxtile2vars, (expix, tileix))
            push!(varixs, varix)
        end
        # each varixs should be sorted, because they are traversed in asc order when constructing vXmt_els
        vars2maskxtiles = Vector{Tuple{Vector{Int}, Vector{Tuple{Int, Int}}}}()
        for (maskxtile, varixs) in pairs(maskxtile2vars)
            v2mtpos = searchsortedfirst(vars2maskxtiles, (varixs, 0), by=first)
            if v2mtpos > length(vars2maskxtiles) || first(vars2maskxtiles[v2mtpos]) != varixs
                insert!(vars2maskxtiles, v2mtpos, (varixs, [maskxtile]))
            else
                # FIXME use OrderedDict to avoid searchsorted?
                maskxtiles = last(vars2maskxtiles[v2mtpos])
                insert!(maskxtiles, searchsortedfirst(maskxtiles, maskxtile), maskxtile)
            end
        end
        maskxtile2ix = Dict{Tuple{Int, Int}, Int}()
        for (i, (_, maskxtiles)) in enumerate(vars2maskxtiles)
            for maskxtile in maskxtiles
                maskxtile2ix[maskxtile] = i
            end
        end
        # replace expix by maskxtile index in vxmxt_els
        @inbounds for (i, (varix, expix, tileix, v)) in enumerate(vXmt_els)
            vXmt_els[i] = (varix, maskxtile2ix[(expix, tileix)], 0, v)
        end
        # sort vxmxt_els by varix, the by maskxtile index to match SpraseMatrixCSC order
        sort!(vXmt_els, lt=(a,b) -> (a[1] < b[1]) || ((a[1] == b[1]) && (a[2] < b[2])))
        # construct sparse matrix
        colptrs = Vector{Int}(undef, length(var2setix) + 1)
        rowvals = Vector{Int}(undef, mapreduce(x -> length(first(x)), +, vars2maskxtiles, init=0))
        nzvals = fill!(Vector{Float64}(undef, length(rowvals)), 0)
        last_coord = (0, 0)
        i = 0
        @inbounds for (varix, maskxtileix, _, v) in vXmt_els
            if last_coord != (varix, maskxtileix)
                @assert i < length(rowvals)
                i += 1
                if last_coord[1] != varix
                    @assert last_coord[1] < varix
                    colptrs[last_coord[1]+1:varix] .= i
                else
                    @assert maskxtileix > last_coord[2]
                end
                last_coord = (varix, maskxtileix)
            end
            rowvals[i] = maskxtileix
            nzvals[i] += v
        end
        @assert i == length(rowvals)
        colptrs[last_coord[1]+1:end] .= length(rowvals) + 1
        maskxtileXvar = SparseMatrixCSC{Float64, Int}(length(vars2maskxtiles), length(var2setix),
                                                      colptrs, rowvals, nzvals)
    end
    return tileXvar, maskxtileXvar, tile2nmasked, tile2nunmasked
end

varXvar_score(setXset::Real, params::CoverParams, scale::Bool = false) =
    ifelse(scale, -(-setXset)^params.set_shape * params.setXset_factor, setXset)

function varXvar_scores(mosaic::AbstractWeightedSetMosaic, params::CoverParams,
                        var2setix::AbstractVector{Int} = var2set(mosaic),
                        scale::Bool = false)
    glob_setixs = mosaic.loc2glob_setix[var2setix]
    vXv_scores = varXvar_score.(mosaic.original.setXset_scores[glob_setixs, glob_setixs], Ref(params), scale)
    for i in eachindex(var2setix) # don't consider self-intersections
        @inbounds vXv_scores[i, i] = zero(eltype(vXv_scores))
    end
    # find minimum finite varXvar_scores element
    vXv_min = Inf
    @inbounds for vXv in vXv_scores
        if isfinite(vXv) && vXv < vXv_min
            vXv_min = vXv
        end
    end
    # replace infinite varXvar score with the minimal finite setXset score
    vXv_subst = varXvar_score(1.25*vXv_min, params, true)
    vXv_cartixs = CartesianIndices(vXv_scores)
    @inbounds for i in eachindex(vXv_scores)
        if !isfinite(vXv_scores[i])
            v1, v2 = Tuple(vXv_cartixs[i])
            @warn "var[$v1]×var[$v2] score is $(vXv_scores[i])"
            if vXv_scores[i] < 0.0
                vXv_scores[i] = vXv_subst
            end
        end
    end
    return vXv_scores
end

# detailed cover score:
# 1 => covered set scores,
# 2 => covered setXset penalties,
# 3 => uncovered masked elements penalties
# 4 => covered unmasked elements penalties
const DetailedScore = NTuple{4, Float64}

aggscore(s::DetailedScore, params::CoverParams) =
    s[1] +
    params.setXset_factor * s[2] +
    params.uncovered_factor * s[3] +
    params.covered_factor * s[4]

"""
    AbstractCoverProblem{T}

Optimal Enriched-Set Cover problem -- choose the sets from the collection `𝒞` to cover
the masked(selected) elements `M`.
The optimal sets cover `C = {c₁, c₂, ..., cₙ} ⊂ 𝒞` has to deliver 3 goals:
* be relevant (i.e. minimize the P-values of `M` and `cᵢ` sets overlap)
* be minimal (i.e. minimize the number of sets in `C`)
* be non-redundant (i.e. minimize the P-values of the pairwise non-overlap of `C` sets with each other).

Fuzzy set selection is possible -- each set is assigned a weight from `[0, 1]` range.
"""
abstract type AbstractCoverProblem{T} end

"""
Total number of problem variables.
"""
nvars(problem::AbstractCoverProblem) = length(problem.var_scores)

# clamps weights to 0..1 range and returns sum(|w - fixed_w|)
function fix_cover_weights!(weights::Matrix{Float64})
    pen = 0.0
    @inbounds for i in eachindex(weights)
        w = uncov_probs[i]
        w_new = clamp(w, 0.0, 1.0)
        pen += abs2(w_new - w)
        weights[i] = w_new
    end
    return pen
end

__check_vars(w::AbstractVector{Float64}, problem::AbstractCoverProblem) =
    length(w) == nvars(problem) ||
        throw(DimensionMismatch("Incorrect length of parameters vector: $(length(w)) ($(nvars(problem)) expected)"))

function selectvars(problem::AbstractCoverProblem,
                    _::AbstractWeightedSetMosaic,
                    weights::AbstractVector{Float64})
    @assert length(weights) == nvars(problem)
    return findall(>(problem.params.min_weight), weights) # > to make it work with min_weight=0
end

abstract type AbstractOptimizerParams{P <: AbstractCoverProblem} end;

problemtype(_::AbstractOptimizerParams{P}) where P = P


text:
#=##############################################################################################
Filename: ccblade.jl
Author: Ryan Anderson
Contact: rymanderson@gmail.com
README: this file is meant to provide convenience functions for building CCBlade rotor objects.
        This includes running all necessary airfoil analysis and store files in an accessible
        directory structure.
=###############################################################################################

"""
    CCBladeSystem{V1,V2,V3,V4,V5,V6,V7}

Defines CCBlade rotor geometry.

# Fields:

* `rotors::Vector{CCBlade.Rotor}`
* `sectionlists::Vector{Vector{CCBlade.Section}}`
* `rlists::Vector{Vector{Float64}}`
* `index::Vector{Int}`
* `positions::Vector{Vector{Float64}}` : x-y-z position in the body frame for each rotor
* `orientations::Vector{Vector{Float64}}` : unit vectors defining the axis of rotation in the direction of positive thrust for each rotor
* `spindirection::Vector{Bool}` : `true` if the rotation vector is in the same direction as positive thrust (by the right-hand rule)

"""
struct CCBladeSystem{TF,TF1,TF2,TF3,TAF}
        # T1 <: Union{Nothing, CC.Mach_correction}, T2 <: Union{Nothing, CC.ReCorrection},
        # T3 <: Union{Nothing, CC.RotationCorrection}, T4 <: Union{Nothing, CC.TipCorrection}}
    rotors::Vector{CC.Rotor}#{TF, Int64, Bool, T1, T2, T3, T4}}
    sectionlists::Vector{Vector{CCBlade.Section{TF1, TF2, TF3, TAF}}}
    rlists::Vector{Vector{TF}}
    index::Vector{Int64}
    positions::Vector{Vector{TF}}
    orientations::Vector{Vector{TF}}
    spin_directions::Vector{Bool}
end

# AircraftSystems.CCBladeSystem(
#     ::Array{CCBlade.Rotor,1},
#     ::Array{Array{CCBlade.Section{Float64,Float64,Float64,CCBlade.AlphaReAF{Float64,String}},1},1},
#     ::Array{Array{Float64,1},1},
#     ::Array{Int64,1},
#     ::Array{Array{Float64,1},1},
#     ::Array{Array{Float64,1},1},
#     ::Array{Bool,1})

"""
    CCBladeSystem(nblades_list, rhub_list, rtip_list, radii_list, chords_list, twists_list,
        airfoilcontours_list, airfoilnames_list, index, positions, orientations, spin_directions,
        Res_list = [fill([5e4, 1e5, 1e6], length(radii_list[i])) for i in 1:length(nblades_list)];
        kwargs...)

Convenience constructor for `CCBladeSystem`.

# Arguments:

* `nblades_list::Vector{Int64}`
* `rhub_list::Vector{Float64}`
* `rtip_list::Vector{Float64}`
* `radii_list::Vector{Vector{Float64}}`
* `chords_list::Vector{Vector{Float64}}`
* `twists_list::Vector{Vector{Float64}}`
* `airfoilcontours_list::Vector{Vector{String}}`
* `airfoilnames_list::Vector{Vector{String}}`
* `index::Vector{Int64}`
* `positions::Vector{Vector{Float64}}`
* `orientations::Vector{Vector{Float64}}`
* `spindirection::Vector{Bool}`

# Keyword Arguments:

* `Res_list`: vector of vectors describes Reynold's numbers at which to run Xfoil for each rotor section

# Optional Arguments:

* `mach_correction` : CCBlade Mach correction object
* `skipstart = 1`: number of header lines in airfoil contour file
* `xfoil_alpha = range(-20.0, stop=20.0, length= 161)`: angles of attack at which to run Xfoil
* `M = 0`: Mach number (WARNING: REQUIRES FURTHER DEVELOPMENT FOR NONZERO MACH NUMBER)
* `ν = 1.5e-5`: kinematic viscosity of the freestream
* `Re_digits = -4`: number of digits to round Reynolds number
* `runxfoil=true`: toggle running Xfoil
* `xfoil_iter = 300`: max iterations per angle of attack for Xfoil
* `xfoil_npan = 200`: number of panels for Xfoil to interpolate
* `xfoil_clmaxstop = true`: automatically detect cl_max and stop Xfoil sweep
* `xfoil_clminstop = true`: automatically detect cl_min and stop Xfoil sweep
* `viternaextrapolation=true`: apply the viterna extrapolation
* `rotationcorrection=true`: apply rotational correction
* `rotationcorrection_J = 2.0`: advance ratio at which to apply rotational correction
* `radians = false`: units of written files
* `savefiles = true`: save airfoil files to be read later
* `useoldfiles = true`: read old files if available
* `polar_directory = joinpath(topdirectory, "data", "airfoil", "polars", TODAY)`: save directory for airfoil files
* `plotoutput = true`: plot airfoil data; quickly show which polars may need to be modified
* `saveplots = true`: save airfoil plots
* `plot_extension = ".pdf"`: file extension for saved airfoil plots
* `verbose = true`: verbose output
* `v_lvl = 0`: number of indentations in verbose output

# Returns:

* `rotor_system::CCBladeSystem`: constructs a `CCBladeSystem` struct

"""
function CCBladeSystem(rotors, sectionlists, index, positions, orientations, spin_directions)

    # check input sizes
    for (iorientation, orientation) in enumerate(orientations)
        @assert isapprox(LA.norm(orientation), 1.0) "`orientations` must be unit vectors; check the $(iorientation)th vector"
    end

    @assert length(index) == length(positions) "length of `index` must match the length of `positions`"
    @assert length(index) == length(orientations) "length of `index` must match the length of `orientations`"
    @assert maximum(index) <= length(rotors) "maximum rotor index cannot exceed the number of rotors"

    rlists = sections2radii.(sectionlists)

    return CCBladeSystem(rotors, sectionlists, rlists, index, positions, orientations, spin_directions)
end

function CCBladeSystem(nblades_list, rhub_list, rtip_list, mach_correction,
    sectionlists, rlists, index, positions, orientations, spin_directions; kwargs...
)
    # build CCBlade.Rotor objects
    rotors = build_rotors(nblades_list, rhub_list, rtip_list, mach_correction; kwargs...)

    return CCBladeSystem(rotors, sectionlists, rlists, index, positions, orientations, spin_directions)
end

function CCBladeSystem(rotors,
    radii_list, chords_list, twists_list, ccblade_af_objects_list,
    index, positions, orientations, spin_directions
)
    # build CCBlade.Section lists
    sectionlists = build_sections(radii_list, chords_list, twists_list, ccblade_af_objects_list)

    return CCBladeSystem(rotors, sectionlists, index, positions, orientations, spin_directions)
end

function CCBladeSystem(nblades_list, rhub_list, rtip_list, mach_correction,
    radii_list, chords_list, twists_list, ccblade_af_objects_list,
    index, positions, orientations, spin_directions;
        kwargs...
)
    # build CCBlade.Rotor objects
    rotors = build_rotors(nblades_list, rhub_list, rtip_list, mach_correction; kwargs...)

    return CCBladeSystem(rotors,
        radii_list, chords_list, twists_list, ccblade_af_objects_list,
        index, positions, orientations, spin_directions
    )
end

function CCBladeSystem(nblades_list, rhub_list, rtip_list, mach_correction,
    radii_list, chords_list, twists_list,
    Res_lists, Ms_lists, contour_paths_list, uncorrected_polar_paths_lists, corrected_polar_paths_lists,
    index, positions, orientations, spin_directions;
        kwargs...
)
    # build af objects
    ccblade_af_objects_list = generate_polars.(rtip_list, radii_list, chords_list, twists_list, Res_lists, Ms_lists, contour_paths_list, uncorrected_polar_paths_lists, corrected_polar_paths_lists; kwargs...)

    return CCBladeSystem(nblades_list, rhub_list, rtip_list, mach_correction,
        radii_list, chords_list, twists_list, ccblade_af_objects_list,
        index, positions, orientations, spin_directions
    )
end

function CCBladeSystem(nblades_list, rhub_list, rtip_list, mach_correction,
    radii_list, chords_list, twists_list,
    Res_lists, Ms_lists, contour_paths_list::Vector{String},
    index, positions, orientations, spin_directions;
        polar_directory = joinpath(topdirectory, "data", "airfoil", "polars", TODAY),
        kwargs...
)

    uncorrected_polar_paths_lists = Vector{Array{String,2}}(undef, length(nblades_list))
    corrected_polar_paths_lists = Vector{Array{String,2}}(undef, length(nblades_list))
    for i_rotor in 1:length(nblades_list)
        polar_base_names = joinpath(polar_directory, splitext(splitpath(contour_paths_list[i_rotor])[end])[1])
        uncorrected_polar_paths_lists[i_rotor] = airfoil_2_polar_names(polar_base_names, Res_lists[i_rotor], Ms_lists[i_rotor];
            viterna_extrapolation=false, rotation_correction=false,
            extension=".txt"
        )
        corrected_polar_paths_lists[i_rotor] = airfoil_2_polar_names(polar_base_names, Res_lists[i_rotor], Ms_lists[i_rotor];
            viterna_extrapolation=true, rotation_correction=true,
            extension=".txt"
        )
    end

    return CCBladeSystem(nblades_list, rhub_list, rtip_list, mach_correction,
        radii_list, chords_list, twists_list,
        Res_lists, Ms_lists, contour_paths_list, uncorrected_polar_paths_lists, corrected_polar_paths_lists,
        index, positions, orientations, spin_directions;
            kwargs...
    )
end

function CCBladeSystem(nblades_list, rhub_list, rtip_list, mach_correction,
    radii_list, chords_list, twists_list,
    Res_lists, Ms_lists, contour_paths_list,
    index, positions, orientations, spin_directions;
        polar_directory = joinpath(topdirectory, "data", "airfoil", "polars", TODAY),
        kwargs...
)

    uncorrected_polar_paths_lists = Vector{Vector{Array{String,2}}}(undef, length(nblades_list))
    corrected_polar_paths_lists = Vector{Vector{Array{String,2}}}(undef, length(nblades_list))
    for i_rotor in 1:length(nblades_list)
        polar_base_names = joinpath.(Ref(polar_directory), [splitext(splitpath(polar_i)[end])[1] for polar_i in contour_paths_list[i_rotor]])
        uncorrected_polar_paths_lists[i_rotor] = airfoil_2_polar_names(polar_base_names, Res_lists[i_rotor], Ms_lists[i_rotor];
            viterna_extrapolation=false, rotation_correction=false,
            extension=".txt"
        )
        corrected_polar_paths_lists[i_rotor] = airfoil_2_polar_names(polar_base_names, Res_lists[i_rotor], Ms_lists[i_rotor];
            viterna_extrapolation=true, rotation_correction=true,
            extension=".txt"
        )
    end

    return CCBladeSystem(nblades_list, rhub_list, rtip_list, mach_correction,
        radii_list, chords_list, twists_list,
        Res_lists, Ms_lists, contour_paths_list, uncorrected_polar_paths_lists, corrected_polar_paths_lists,
        index, positions, orientations, spin_directions;
            kwargs...
    )
end

"""
    sections2radii(sections)

# Arguments

* `sections`

# Returns

"""
@inline function sections2radii(sections)
    return [section.r for section in sections]
end

@inline function build_rotors(nblades_list, rhub_list, rtip_list, mach_correction = nothing; kwargs...)

    @assert length(nblades_list) == length(rhub_list) "list of nblades and list of hub radii must be the same length"
    @assert length(nblades_list) == length(rtip_list) "list of nblades and list of tip radii must be the same length"

    rotors = Vector{CC.Rotor}(undef,length(nblades_list))
    for i in 1:length(rotors)
        # note: Mach correction on the fly not required
        # rhub_list[i] =
        rotors[i] = CC.Rotor(rhub_list[i], rtip_list[i], nblades_list[i]; precone=0.0, turbine=false, mach=mach_correction, re=nothing, rotation=nothing, tip=CC.PrandtlTipHub()) #mach=CC.PrandtlGlauert()
    end
    return rotors
end

@inline function build_sections(radii_list, chords_list, twists_list, ccblade_af_objects_list)

    @assert length(chords_list) == length(twists_list) "chords_list and twists must be the same length"
    @assert length(chords_list) == length(ccblade_af_objects_list) "chords_list and ccblade_af_objects_list must be the same length"

    return rotor_2_sections.(radii_list, chords_list, twists_list, ccblade_af_objects_list)
end

@inline function rotor_2_sections(radii, chords, twists, ccblade_af_objects)
    sections = CC.Section.(radii, chords, twists, ccblade_af_objects)
end

"""
    generate_polars(rtip, radii, chords, twists, Res_list, Ms_list,
        contour_paths, uncorrected_polar_paths_list, corrected_polar_paths_list;
            kwargs...
    )

Builds airfoil objects for each radial section of a rotor.

# Arguments

* rtip::Float64
* radii::Vector{Float64}
* chords::Vector{Float64}
* twists::Vector{Float64}
* Res_list::Vector{Vector{Float64}}
* Ms_list::Vector{Vector{Float64}}
* contour_paths::Vector{String}
* uncorrected_polar_paths_list::Vector{Array{String,2}}
* corrected_polar_paths_list::Vector{Array{String,2}}

"""
function generate_polars(rtip, radii, chords, twists, Res_list, Ms_list, contour_paths, uncorrected_polar_paths_list, corrected_polar_paths_list; kwargs...)
    cr75 = FM.linear(radii ./ rtip, chords, 0.75) / FM.linear(radii ./ rtip, radii, 0.75)

    polars = rotor_2_polars(radii, chords, cr75, Res_list, Ms_list, contour_paths, uncorrected_polar_paths_list, corrected_polar_paths_list;
        use_old_files=true, viterna_extrapolation=true, rotation_correction=true, radians_polar=false,
        plot_polars=true, save_figure=true, extension=".png", close_figure=true,
        kwargs...
    )

    return polars
end

"""
    rotor_2_polars(radii, chords, cr75, Res_list, Ms_list, contour_paths, uncorrected_polar_paths_list, corrected_polar_paths_list;
        use_old_files=true, viterna_extrapolation=true, rotation_correction=true, radians_polar=false,
        plot_polars=true, save_figure=true, extension=".png", close_figure=true,
        kwargs...
    )

Build a vector of airfoil functions for each rotor section.

# Arguments:

* `radii::Vector{Float64}` : radii of each rotor section in absolute units
* `chords::Vector{Float64}` : chord of each rotor section in absolute units
* `cr75::Float64` : local chord to tip radius at 75% r/R location
* `Res_list::Vector{Vector{Int64}}`
* `Ms_list::Vector{Vector{Int64}}`
* `contour_paths::Vector{String}` : vector of paths to .dat airfoil contour file
* `uncorrected_polar_paths_list::Vector{Array{String,2}}` : [i][j,k] is the path to the uncorrected polar file of the ith radial station at the jth Reynold's number, kth Mach number
* `corrected_polar_paths_list::Vector{Array{String,2}}` : [i][j,k] is the path to the corrected polar file of the ith radial station at the jth Reynold's number, kth Mach number

# Keyword Arguments:

* `use_old_files=true` : re-use existing polar files
* `viterna_extrapolation=true` : perform the Viterna extrapolation on airfoil data
* `rotation_correction=true` : perform a rotation correction to airfoil data
* `radians_polar=false` : `true` indicates aoa's of polars is in radians
* `plot_polars=true` : plot polars
* `save_figure=true` : save polar plots
* `extension=".png"` : extension of polar plots
* `close_figure=true` : close polar plots after generation

# Returns:

* `airfoil_objects`: vector of `<:CCBlade.AFType` objects

"""
function rotor_2_polars(radii, chords, cr75, Res, Ms, contour_paths::String, uncorrected_polar_paths, corrected_polar_paths;
    use_old_files=true, use_old_corrected_files = true, viterna_extrapolation=true, rotation_correction=true, radians_polar=false,
    plot_polars=true, save_figure=true, extension=".png", close_figure=true,
    kwargs...
)

    @assert length(radii) == length(chords) "length of radii and chords inconsistent"
    @assert length(uncorrected_polar_paths) == length(Res) * length(Ms) "size of uncorrected_polar_paths is $(size(uncorrected_polar_paths)); expected $((length(Res), length(Ms)))"
    @assert length(corrected_polar_paths) == length(Res) * length(Ms) "size of corrected_polar_paths is $(size(corrected_polar_paths)); expected $((length(Res), length(Ms)))"

    # TODO: re-write this to use multiple dispatch rather than if statements
    if length(Res) == 1
        if length(Ms) == 1
            af_type = CC.AlphaAF
        else
            af_type = CC.AlphaMachAF
        end
    else
        if length(Ms) == 1
            af_type = CC.AlphaReAF
        else
            af_type = CC.AlphaReMachAF
        end
    end

    # check if corrected files exist
    # corrected_polar_paths = add_tags.(uncorrected_polar_paths; viterna_extrapolation, rotation_correction)
    i_existing_files_corrected = isfile.(corrected_polar_paths)
    if !(use_old_files && prod(i_existing_files_corrected) && use_old_corrected_files) # if notall files exist
        # write files
        airfoil_2_files(Res, Ms, contour_paths, uncorrected_polar_paths, corrected_polar_paths;
            radians_polar, cr75, use_old_files, use_old_corrected_files,
            plot_polars, save_figure, extension, close_figure,
            viterna_extrapolation, rotation_correction,
            kwargs...
        )
    end

    airfoil_objects = fill(Res_Ms_2_polar(Res, Ms, af_type, corrected_polar_paths, radians_polar), length(radii))

    return airfoil_objects
end

function rotor_2_polars(radii, chords, cr75, Res_list, Ms_list, contour_paths, uncorrected_polar_paths_list, corrected_polar_paths_list;
    use_old_files=true, use_old_corrected_files = true, viterna_extrapolation=true, rotation_correction=true, radians_polar=false,
    plot_polars=true, save_figure=true, extension=".png", close_figure=true,
    kwargs...
)

    @assert length(radii) == length(chords) "length of radii and chords inconsistent"
    @assert length(radii) == length(Res_list) "length of radii ($(length(radii))) and Res_list ($(length(Res_list))) inconsistent"
    @assert length(radii) == length(contour_paths) "length of radii and contour_paths inconsistent"
    @assert length(radii) == length(uncorrected_polar_paths_list) "length of radii and uncorrected_polar_files_list inconsistent"
    @assert size(uncorrected_polar_paths_list[1]) == (length(Res_list[1]), length(Ms_list[1])) "size of uncorrected_polar_files_list[1] is $(size(uncorrected_polar_files[1])); expected $((length(Res_list[1]), length(Ms_list[1])))"
    @assert size(corrected_polar_paths_list[1]) == (length(Res_list[1]), length(Ms_list[1])) "size of corrected_polar_files_list[1] is $(size(corrected_polar_files[1])); expected $((length(Res_list[1]), length(Ms_list[1])))"

    # TODO: re-write this to use multiple dispatch rather than if statements
    if length(Res_list[1]) == 1
        if length(Ms_list[1]) == 1
            af_type = CC.AlphaAF
        else
            af_type = CC.AlphaMachAF
        end
    else
        if length(Ms_list[1]) == 1
            af_type = CC.AlphaReAF
        else
            af_type = CC.AlphaReMachAF
        end
    end

    if !(typeof(contour_paths[1]) <: AbstractArray) && prod(contour_paths .== contour_paths[1]) && # single airfoil
        prod([isequal(Res_list[1], Res_list[i]) for i in 2:length(Res_list)]) && # same Res
        prod([isequal(Ms_list[1], Ms_list[i]) for i in 2:length(Ms_list)]) # same Ms

        airfoil_objects = fill(Res_Ms_2_polar(Res_list[1], Ms_list[1], af_type, corrected_polar_paths_list[1], radians_polar))
    else
        airfoil_objects = Array{af_type,1}(undef,length(radii))

        # get polars for each radial station
        for (i_radius, radius) in enumerate(radii)
            Res = Res_list[i_radius]
            Ms = Ms_list[i_radius]
            uncorrected_polar_paths = uncorrected_polar_paths_list[i_radius]
            corrected_polar_paths = corrected_polar_paths_list[i_radius]

            # check if corrected files exist
            # corrected_polar_paths = add_tags.(uncorrected_polar_paths; viterna_extrapolation, rotation_correction)
            i_existing_files_corrected = isfile.(corrected_polar_paths)
            if use_old_files && prod(i_existing_files_corrected) && use_old_corrected_files # if all files exist
                airfoil_objects[i_radius] = Res_Ms_2_polar(Res, Ms, af_type, corrected_polar_paths, radians_polar)
            else
                # write files
                airfoil_2_files(Res, Ms, contour_paths[i_radius], uncorrected_polar_paths, corrected_polar_paths;
                    radians_polar, cr75, use_old_files, use_old_corrected_files,
                    plot_polars, save_figure, extension, close_figure,
                    viterna_extrapolation, rotation_correction,
                    kwargs...
                )

                # build and store objects
                airfoil_objects[i_radius] = Res_Ms_2_polar(Res, Ms, af_type, corrected_polar_paths, radians_polar)
            end
        end
    end

    return airfoil_objects
end

function Res_Ms_2_polar(Res, Ms, af_type, corrected_polar_paths, radians_polar)
    if length(Res) == 1 && length(Ms) == 1
        return af_type(corrected_polar_paths[1]; radians = radians_polar)
    elseif length(Res) == 1 || length(Ms) == 1
        return af_type(corrected_polar_paths[:]; radians = radians_polar)
    else
        return af_type(corrected_polar_paths; radians = radians_polar)
    end
end

"""
    airfoil_2_files(Res, Ms, contour_path, uncorrected_polar_paths, corrected_polar_paths;
        viterna_extrapolation=true, rotation_correction=true,
        skipstart=1, # 1 line header in .dat contour file
        iter=300, npan=200, clmaxstop=true, clminstop=true,
        xfoil_alpha = range(-20.0, stop=20.0, length= 161), radians_xfoil = false, radians_polar = true,
        cr75 = 0.0,
        use_old_files=true,
        plot_polars=true, save_figure=true, extension=".pdf", close_figure=true,
        verbose=true, v_lvl=0,
        kwargs...
    )

Writes uncorrected and corrected polar files from an airfoil contour file at designated Reynolds and Mach numbers.

# Arguments:

* `Res::Vector` : vector of Reynolds numbers at which to run Xfoil
* `Ms::Vector` : vector of Mach numbers at which to run Xfoil
* `contour_path::String` : path to airfoil contour file
* `uncorrected_polar_paths::Array{String,2}` : [i,j]th element is the path to the uncorrected polar of the ith Reynold's number, jth Mach number
* `corrected_polar_paths::Array{String,2}` : [i,j]th element is the path to the corrected polar of the ith Reynold's number, jth Mach number

# Keyword Arguments:

* `skipstart = 1` : number of header lines in contour file
* `xfoil_alpha = range(-20.0, stop=20.0, length= 161)` : angles of attack at which to run Xfoil
* `M = 0` : Mach number
* `ν = NU` : kinematic viscocity
* `Re_digits = -4` : number of digits to round the Reynolds number
* `xfoil_iter = 300` : max iterations for Xfoil
* `xfoil_npan = 200` : number of panels for Xfoil
* `xfoil_clmaxstop = true` : stop alpha sweep when cl reaches a maximum
* `xfoil_clminstop = true` : stop alpha sweep when cl reaches a minimum
* `radians_xfoil = false` : units of aoa's input to xfoil (`xfoil_alpha`)
* `radians_polar = false` : desired units for polar files
* `cr75 = 0.0` : r/R=0.75 local chord over local radius (used for rotation correction)
* `use_old_files = true` : use old polar files (if they exist)
* `plot_polars = true` : toggle ploting polars
* `save_figure = true` : toggle saving plots
* `extension = ".pdf"` : extension of saved plots
* `close_figure = true` : close the figure after generation
* `verbose = true` : toggle verbose output
* `v_lvl = 0` : set verbosity level

# Returns:

* nothing

"""
function airfoil_2_files(Res, Ms, contour_path, uncorrected_polar_paths, corrected_polar_paths;
    viterna_extrapolation=true, rotation_correction=true,
    skipstart=1, # 1 line header in .dat contour file
    iter=300, npan=200, clmaxstop=true, clminstop=true,
    xfoil_alpha = range(-20.0, stop=20.0, length= 161), radians_xfoil = false, radians_polar = true,
    cr75 = 0.0, rotation_correction_J = 2.0,
    use_old_files=true, use_old_corrected_files = true,
    plot_polars=true, save_figure=true, extension=".pdf", close_figure=true,
    verbose=true, v_lvl=0, kwargs...
)

    @assert isfile(contour_path) "`contour_path` $(contour_path) does not exist"
    @assert length(uncorrected_polar_paths) == length(Res) * length(Ms) "size of uncorrected_polar_paths is $(size(uncorrected_polar_paths)); expected ($(length(Res)), $(length(Ms)))"
    @assert length(corrected_polar_paths) == length(Res) * length(Ms) "size of corrected_polar_paths is $(size(corrected_polar_paths)); expected ($(length(Res)), $(length(Ms)))"

    if verbose; println("\t"^v_lvl, "Reading $contour_path\n"); end
    xy = DF.readdlm(contour_path, ',', Float64, '\n'; skipstart = skipstart)
    percussive_maintenance = true
    α_ref = nothing
    for (i_M, M) in enumerate(Ms)
        for (i_Re, Re) in enumerate(Res)
            if !isfile(uncorrected_polar_paths[i_Re, i_M]) || !use_old_files # if the file already exists
                α, cl, cd = xfoil_2_file(uncorrected_polar_paths[i_Re, i_M], xy[:,1], xy[:,2], Re, M;
                    xfoil_alpha, radians_xfoil, radians_polar,
                    iter, npan, percussive_maintenance, verbose, clmaxstop, clminstop,
                )

                if plot_polars
                    plot_airfoil(α, cl, cd, uncorrected_polar_paths[i_Re, i_M];
                        radians=radians_polar,
                        viterna_extrapolation=false, rotation_correction=false,
                        save_figure, extension, clear_figure = true, close_figure
                    )
                    plot_airfoil(α, cl, cd, corrected_polar_paths[i_Re, i_M];
                        radians=radians_polar,
                        viterna_extrapolation=false, rotation_correction=false,
                        save_figure = save_figure && !viterna_extrapolation && !rotation_correction,
                        extension, clear_figure = true, close_figure = false
                    )
                end
            else # read the file
                polar = CC.AlphaAF(uncorrected_polar_paths[i_Re, i_M]; radians = radians_polar)
                α = polar.alpha
                cl = polar.cl
                cd = polar.cd
            end
            if !isfile(corrected_polar_paths[i_Re, i_M]) || !use_old_corrected_files # if the file already exists, skip this
                convert_to_radians = radians_polar ? 1.0 : pi / 180.0
                if viterna_extrapolation
                    α, cl, cd = CC.viterna(α * convert_to_radians, cl, cd, cr75) # only use converged values
                    α .*= 180/pi
                    if plot_polars
                        plot_airfoil(α, cl, cd, corrected_polar_paths[i_Re, i_M];
                        radians = radians_polar,
                        viterna_extrapolation=true, rotation_correction=false,
                        save_figure = save_figure && !rotation_correction,
                        extension, clear_figure = false, close_figure = false
                        )
                    end
                end

                if rotation_correction
                    for (i_α, this_α) in enumerate(α)
                        # rotation_correction(rc::RotationCorrection, cl, cd, cr, rR, tsr, alpha, phi=alpha, alpha_max_corr=30*pi/180)
                        cl[i_α], cd[i_α] = CC.rotation_correction(CC.DuSeligEggers(), cl[i_α], cd[i_α], cr75, 0.75, pi/rotation_correction_J, this_α * convert_to_radians)
                    end
                    if plot_polars
                        plot_airfoil(α, cl, cd, corrected_polar_paths[i_Re, i_M];
                        radians = radians_polar,
                        viterna_extrapolation, rotation_correction=true,
                        save_figure, extension, clear_figure = false, close_figure = false
                        )
                    end
                end

                # interpolate to consistent aoas
                if i_M == i_Re && i_M == 1; global α_ref = deepcopy(α); end
                cl_aligned = FM.linear(α, cl, α_ref)
                cd_aligned = FM.linear(α, cd, α_ref)
                if plot_polars
                    plot_airfoil(α_ref, cl_aligned, cd_aligned, corrected_polar_paths[i_Re, i_M];
                    radians = radians_polar,
                    viterna_extrapolation, rotation_correction, aligned = true,
                    save_figure, extension, clear_figure = false, close_figure
                    )
                end

                corrected_polar = CC.AlphaAF(α_ref, cl_aligned, cd_aligned, splitpath(corrected_polar_paths[i_Re, i_M])[end], Re, M) # TODO: maybe Re * 1.0, mach * 1.0
                CC.write_af(corrected_polar_paths[i_Re, i_M], corrected_polar)
            end
        end
    end

    return nothing
end

"""

# Arguments

* `radians_xfoil::Bool` : true if the aoas provided are in radians
* `radians_polar::Bool` : true if aoas of the polar object should be in radians

"""
function xfoil_2_file(polar_path, x, y, Re, mach;
    xfoil_alpha = range(-20.0, stop=20.0, length= 161), radians_xfoil = false, radians_polar = true,
    iter=300, npan=200, percussive_maintenance=true, verbose=true, clmaxstop=true, clminstop=true,
)
    if verbose; println("\n", "==================== Re = $Re, M = $mach ===================="); end

    # write uncorrected polar
    angle_conversion_xfoil = radians_xfoil ? 180.0 / pi : 1
    cl_raw, cd_raw, cdp, cm, conv = XF.alpha_sweep(x, y, xfoil_alpha .* angle_conversion_xfoil, Re;
                                iter,
                                npan,
                                mach,
                                percussive_maintenance,
                                printdata=verbose,
                                zeroinit=true,
                                clmaxstop,
                                clminstop)

    # create object and save files
    angle_conversion_polar = radians_polar ? pi/180 : 1.0
    α = xfoil_alpha[conv] .* angle_conversion_polar
    cl = cl_raw[conv]
    cd = cd_raw[conv]
    polar = CC.AlphaAF(α, cl, cd, splitpath(polar_path)[end], Re, mach) # TODO: maybe Re * 1.0, mach * 1.0
    if verbose; println("Writing to $polar_path"); end
    println("pathcheck = $(splitpath(polar_path)[1])")
    if !isdir(splitpath(polar_path)[1]); mkpath(splitpath(polar_path)[1]); end
    CC.write_af(polar_path, polar)

    return α, cl, cd
end

"""
    airfoil_2_polar_names(airfoilname, Res, Ms;
        M=0, viternaextrapolation=false, rotationcorrection=false,
        extension=".txt")

Convenience function to create polar file names.

# Arguments:

* `airfoilname::String`
* `Res`: Reynolds numbers
* `Ms`: Mach numbers

# Keyword Arguments:

* `viternaextrapolation::Bool`
* `rotationcorrection::Bool`
* `aoaset::Bool`
* `extension::String`

"""
function airfoil_2_polar_names(airfoil_name::String, Res, Ms;
            viterna_extrapolation=false, rotation_correction=false,
            extension=".txt")

    tags = ""
    if viterna_extrapolation; tags *= "_ext"; end
    if rotation_correction; tags *= "_rot"; end
    [airfoil_name * "_Re$(Int(round(Re; digits=0)))" * "_M$(M)" * tags * extension for Re in Res, M in Ms]
end

"Multiple dispatch for a vector of airfoil_names."
function airfoil_2_polar_names(airfoil_names, Res, Ms;
    viterna_extrapolation=false, rotation_correction=false,
    extension=".txt")
    [airfoil_2_polar_names(airfoil_name, Res[i], Ms[i];
        viterna_extrapolation, rotation_correction, extension) for (i, airfoil_name) in enumerate(airfoil_names)]
end

# """
# Multiple dispatch for M=0.
# """
# function airfoil_2_polar_names(airfoilname, Res;
#             viterna_extrapolation=false, rotation_correction=false,
#             extension=".txt")

#     tags = ""
#     if viterna_extrapolation; tags *= "_ext"; end
#     if rotation_correction; tags *= "_rot"; end

#     [airfoilname * "_Re$(Int(round(Re; digits=0)))" * tags * extension for Re in Res, M in Ms]
# end

function add_tags(path, tags=""; viterna_extrapolation = false, rotation_correction = false)
    file, ext = splitext(path)
    if viterna_extrapolation; tags *= "_ext"; end
    if rotation_correction; tags *= "_rot"; end

    return file * tags * ext
end

"""
    plot_airfoil(α, cl, cd, filename, airfoilname, Re, M;
        radians=false,
        viternaextrapolation=false, rotationcorrection=false,
        savefigure=true, savepath=joinpath(topdirectory, "data", "airfoil", "polars", TODAY),
        extension=".pdf", tag="", clearfigure=true, closefigure=true)

Plots airfoil polars.

# Arguments:

* `α::Vector`: vector of angles of attack
* `cl::Vector`: vector of lift coefficents
* `cd::Vector`: vector of drag coefficents
* `filename::String`: name of the associated CCBlade airfoil polar file
* `airfoilname::String`: name of the airfoil
* `Re`: Reynolds number
* `M`: Mach number

# Keyword Arguments:

* `radians::Bool`
* `viternaextrapolation::Bool`: toggle viterna extrapolation
* `rotationcorrection::Bool`: toggle rotation correction
* `savefigure::Bool`: save the figure
* `savepath::String`: directory to save the figure
* `extension::String`: figure extension
* `tag::String`: string to append to the series legend
* `clearfigure::Bool`: clear the figure before plotting (set false if you want to compare series)
* `closefigure::Bool`:

"""
function plot_airfoil(α, cl, cd, file_name;
    radians=false,
    viterna_extrapolation=false, rotation_correction=false, aligned = false,
    save_figure=true,
    extension=".pdf", tag="", clear_figure=true, close_figure=true
)

    # set up units
    aoaunits = L" [^\circ]"
    conversionfactor = radians ? 180.0 / pi : 1

    # set up figure
    fig = plt.figure(file_name)
    if !clear_figure && length(fig.get_axes()) == 2
        axes = fig.get_axes()
    else
        fig.clear()
        fig.suptitle(splitpath(file_name)[end])
        fig.add_subplot(211)
        fig.add_subplot(212)
        axes = fig.get_axes()
        axes[1].set_ylabel(L"c_l")
        axes[2].set_ylabel(L"c_d")
        axes[2].set_xlabel(L"\alpha" * aoaunits)
    end

    # set labels
    labeltag = ""
    viterna_extrapolation ? labeltag *= "Ext." : labeltag *= "Not ext."
    rotation_correction ? labeltag *= " w/ rot. corr." : labeltag *= " w/o rot. corr."
    if aligned; labeltag *= ", aoa ali."; end
    labeltag *= tag

    # plot data
    axes[1].plot(α .* conversionfactor, cl, label = labeltag)
    axes[2].plot(α .* conversionfactor, cd, label = labeltag)

    # show figure
    axes[1].legend(loc="upper left", bbox_to_anchor=(1.03,0.9), prop=Dict("size" => 9))
    fig.subplots_adjust(left=0.1, right=0.7, top=0.9, bottom=0.1)
    # fig.show()

    # save plot
    if save_figure
        save_name = splitext(file_name)[1] * extension
        save_path = joinpath(splitpath(file_name)[1:end-1]...)
        if !isdir(save_path); mkpath(save_path); end
        fig.savefig(joinpath(save_path, save_name), bbox_inches="tight")
    end

    if close_figure; plt.close(file_name); end

    return nothing
end

"Overload `plotairfoil` to accept a CCBlade.AlphaAF object."
function plot_airfoil(airfoilobject::CCBlade.AlphaAF, file_name; kwargs...)
    plot_airfoil(airfoilobject.alpha, airfoilobject.cl, airfoilobject.cd, file_name; kwargs...)
end

"""
    plotrotor(file_name, rotorname, rs, chords, twists;
        cratio=1.0, savefigure=true, savepath=joinpath(topdirectory, "data", "plots", TODAY),
        extension=".pdf", clearfigure=true, closefigure=true)

Plots rotor geometry.

# Arguments:

* `file_name::String`: used to save the file
* `rotorname::String`: label for the current plot series
* `rs::Vector{Float64}`: radial stations of the rotor
* `chords::Vector{Float64}`: chord lengths at each radial station
* `twists::Vector{Float64}`: twists at each radial station

# Keyword Arguments:

* `cratio`: value between 0 and 1 used to set the plot color
* `savefigure`: save the figure
* `savepath`: directory to save the figure
* `extension`: file extension to the saved file
* `clearfigure`: clear the figure before plotting
* `closefigure`: close the figure upon completion

"""
function plotrotor(file_name, rotorname, rs, chords, twists;
            cratio=1.0, savefigure=true, savepath=joinpath(topdirectory, "data", "plots", TODAY),
            extension=".pdf", clearfigure=true, closefigure=true)

    # set up units
    aoaunits = L" [^\circ]"
    conversionfactor = radians ? 180.0 / pi : 1

    # set up figure
    fig = plt.figure(file_name)
    if clearfigure
        fig.clear()
        fig.suptitle("Rotor Geometry")
        axes = []
        fig.add_subplot(211)
        fig.add_subplot(212)
        axes = fig.get_axes()
        axes[1].set_ylabel(L"c")
        axes[2].set_ylabel(L"\theta")
        axes[2].set_xlabel(L"\alpha" * aoaunits)
    else
        axes = fig.get_axes()
    end

    # plot data
    axes[1].plot(rs, chords, label = rotorname, color=(0.05, 0.85-cratio*0.7, 0.15 + 0.75 * cratio))
    axes[2].plot(rs, twists, label = rotorname, color=(0.05, 0.85-cratio*0.7, 0.15 + 0.75 * cratio))

    # save plot
    if savefigure
        if !isdir(savepath); mkpath(savepath); end
        savename = splitext(file_name)[1] * extension
        fig.savefig(joinpath(savepath, savename), bbox_inches="tight")
    end

    if closefigure; plt.close(file_name); end

    return nothing
end

"""
    rotorinflow(freestream, orientation)

# Arguments:

* `freestream::Freestream`: freestream object
* `orientation::Vector{Float64}`: unit vector aligned with the rotor axis of rotation in the direction of positive thrust

# Returns:

* `vinflow::Float64`: the magnitude of the freestream in the negative thrust direction.
"""
function rotorinflow(freestream, orientation)#, interstream)

    Vinf = freestream2vector(freestream)
    vinflow = -LA.dot(Vinf, orientation)

    return vinflow
end

"""
    rotor2operatingpoints(vinflow, r, omega, environment)

# Arguments:

* `vinflow`
* `r::Vector{Float64}`: vector of a rotor's radial sections
* `omega::Float64`: rotational velocity in rad/s
* `environment::Environment`: `Environment` struct

# Returns:

* `operatingpoints::Vector{CCBlade.OperatingPoint}`: vector of operating points for a single rotor

"""
function rotor2operatingpoints(vinflow, r, omega, environment)

    operatingpoints = CC.simple_op.(vinflow, omega, r, environment.ρ; pitch = 0.0, mu = environment.μ, asound = environment.a, precone = 0.0)

    return operatingpoints
end

"""
    solverotorsnondimensional(rotor_system, omegas, freestream, environment)

# Arguments:

* `rotor_system`
* `omegas`
* `freestream`
* `environment`

# Returns:

* `Js`
* `CTs`
* `CQs`
* `ηs`

"""
function solverotorsnondimensional(rotor_system, omegas, freestream, environment)#, interstream::Interstream)

    # extract rotor info
    orientations = rotor_system.orientations

    # operating conditions
    vinflows = rotorinflow.(Ref(freestream), orientations)

    # pre-allocate solution
    Js, Ts, Qs, Ps, us, vs = solverotors(rotor_system, omegas, freestream, environment)

    CTs = similar(Qs)
    CQs = similar(Qs)
    ηs = similar(Qs)

    # iterate over rotors
    for (i, rotorindex) in enumerate(rotor_system.index)

        # isolate rotor
        rotor = rotor_system.rotors[rotorindex]

        η, CT, CQ = CC.nondim(Ts[i], Qs[i], vinflows[i], omegas[i], environment.ρ, rotor, "propeller")
        CTs[i] = CT
        CQs[i] = CQ
        ηs[i] = η
    end

    return Js, CTs, CQs, ηs
end

"""
In-place version of solverotorsnondimensional.
"""
function solverotorsnondimensional!(Js, Ts, Qs, Ps, CTs, CQs, ηs, rotor_system, omegas, freestream, environment)#, interstream::Interstream)

    # extract rotor info
    orientations = rotor_system.orientations

    # operating conditions
    vinflows = rotorinflow.(Ref(freestream), orientations)

    # pre-allocate solution
    solve_rotors!(Js, Ts, Qs, Ps, rotor_system, omegas, freestream, environment)

    # iterate over rotors
    for (i, rotorindex) in enumerate(rotor_system.index)
        # isolate rotor
        rotor = rotor_system.rotors[rotorindex]
        η, CT, CQ = CC.nondim(Ts[i], Qs[i], vinflows[i], omegas[i], environment.ρ, rotor, "propeller")
        CTs[i] = CT
        CQs[i] = CQ
        ηs[i] = η
    end

    return nothing
end

"""
In-place version of solverotorsnondimensional including `us` and `vs`.
"""
function solverotorsnondimensional!(Js, Ts, Qs, Ps, CTs, CQs, ηs, us, vs, rotor_system, omegas, freestream, environment)#, interstream::Interstream)

    # extract rotor info
    orientations = rotor_system.orientations

    # operating conditions
    vinflows = rotorinflow.(Ref(freestream), orientations)

    # pre-allocate solution
    solve_rotors!(Js, Ts, Qs, Ps, us, vs, rotor_system, omegas, freestream, environment)

    # iterate over rotors
    for (i, rotorindex) in enumerate(rotor_system.index)
        # isolate rotor
        rotor = rotor_system.rotors[rotorindex]
        η, CT, CQ = CC.nondim(Ts[i], Qs[i], vinflows[i], omegas[i], environment.ρ, rotor, "propeller")
        CTs[i] = CT
        CQs[i] = CQ
        ηs[i] = η
    end

    return nothing
end

"""
    solverotors(rotor_system, omegas, freestream, environment)

Solves the rotor system at the specified RPM, freestream, and environment.

# Arguments:

* `rotor_system::CCBladeSystem`: a rotor system
* `omegas::Vector{Float64}: vector of rotational velocities in rad/s of each rotor
* `freestream::Freestream`: Freestream object
* `environment::Environment`: Environment object

# Returns:

* `Js::Vector{Float64}`: vector of advance ratios of each rotor
* `Ts::Vector{Float64}`: vector of thrust values of each rotor
* `Qs::Vector{Float64}`: vector of torque values of each rotor
* `Ps::Vector{Float64}`: vector of power supplied at each rotor
* `us::Vector{Vector{Float64}}`: [i][j] refers to the axial induced velocity at the rotor disk at the jth radial station of the ith rotor
* `vs::Vector{Vector{Float64}}`: [i][j] refers to the swirl induced velocity at the rotor disk at the jth radial station of the ith rotor

"""
function solverotors(rotor_system, omegas, freestream, environment)#, interstream::Interstream)

    # pre-allocate solution
    Js = zeros(length(rotor_system.index))
    Ts = similar(Js)
    Qs = similar(Js)
    Ps = similar(Js)
    us = [zeros(typeof(freestream.vinf), length(radii)) for radii in rotor_system.rlists[rotor_system.index]]
    vs = deepcopy(us)

    solve_rotors!(Js, Ts, Qs, Ps, us, vs, rotor_system, omegas, freestream, environment)

    return Js, Ts, Qs, Ps, us, vs
end

"""
    solve_rotors!(Js, Ts, Qs, us, vs, rotor_system, omegas, freestream, environment)

In-place instance of `solverotors`.

# Arguments:

* `Js::Vector{Float64}` : vector of advance ratios of each rotor
* `Ts::Vector{Float64}` : vector of thrust values of each rotor
* `Qs::Vector{Float64}` : vector of torque values of each rotor
* `us::Vector{Vector{Float64}}` : [i][j] refers to the axial induced velocity at the rotor disk at the jth radial station of the ith rotor
* `vs::Vector{Vector{Float64}}` : [i][j] refers to the swirl induced velocity at the rotor disk at the jth radial station of the ith rotor
* `rotor_system::CCBladeSystem` : a rotor system
* `omegas::Vector{Float64} : vector of rotational velocities in rad/s of each rotor
* `freestream::Freestream` : Freestream object
* `environment::Environment` : Environment object

# Modifies:

* `Js::Vector{Float64}`
* `Ts::Vector{Float64}`
* `Qs::Vector{Float64}`
* `us::Vector{Vector{Float64}}`
* `vs::Vector{Vector{Float64}}`

"""
function solve_rotors!(Js, Ts, Qs, Ps, us, vs, rotor_system, omegas, freestream, environment)#, interstream::Interstream)

    # benchmarks = Vector{Float64}(undef,5)

    # benchmarks[1] = @elapsed begin
    # extract rotor info
    rlists = rotor_system.rlists[rotor_system.index]
    orientations = rotor_system.orientations
    # end

    # benchmarks[2] = @elapsed begin
    # operating conditions
    vinflows = rotorinflow.(Ref(freestream), orientations)
    operatingpoints_list = rotor2operatingpoints.(vinflows, rlists, omegas, Ref(environment))
    # end

    # iterate over rotors
    for (i, rotorindex) in enumerate(rotor_system.index)
        # benchmarks[3] = @elapsed begin
        # isolate rotor
        rotor = rotor_system.rotors[rotorindex]
        sections = rotor_system.sectionlists[rotorindex]

        # isolate operating conditions
        operatingpoints = operatingpoints_list[i]
        # end

        # benchmarks[4] = @elapsed begin
        # solve CCBlade
        outputs = CC.solve.(Ref(rotor), sections, operatingpoints)
        # end

        # benchmarks[5] = @elapsed begin
        T, Q = CC.thrusttorque(rotor, sections, outputs)
        Ts[i] = T
        Qs[i] = Q
        Ps[i] = Q * omegas[i]
        us[i][:] .= outputs.u
        vs[i][:] .= outputs.v
        # end
    end

    # calculate Js
    Js .= vinflows ./ omegas * 2 * pi ./ [2 * rotor.Rtip for rotor in rotor_system.rotors[rotor_system.index]]

    # println("CHRONOS: $benchmarks")
    return nothing
end

"""
    solve_rotors!(Js, Ts, Qs, rotor_system, omegas, freestream, environment)

In-place instance of `solverotors`.

# Arguments:

* `Js::Vector{Float64}` : vector of advance ratios of each rotor
* `Ts::Vector{Float64}` : vector of thrust values of each rotor
* `Qs::Vector{Float64}` : vector of torque values of each rotor
* `us::Vector{Vector{Float64}}` : [i][j] refers to the axial induced velocity at the rotor disk at the jth radial station of the ith rotor
* `vs::Vector{Vector{Float64}}` : [i][j] refers to the swirl induced velocity at the rotor disk at the jth radial station of the ith rotor
* `rotor_system::CCBladeSystem` : a rotor system
* `omegas::Vector{Float64} : vector of rotational velocities in rad/s of each rotor
* `freestream::Freestream` : Freestream object
* `environment::Environment` : Environment object

# Modifies:

* `Js::Vector{Float64}`
* `Ts::Vector{Float64}`
* `Qs::Vector{Float64}`
* `us::Vector{Vector{Float64}}`
* `vs::Vector{Vector{Float64}}`

"""
function solve_rotors!(Js, Ts, Qs, Ps, rotor_system, omegas, freestream, environment)#, interstream::Interstream)

    # extract rotor info
    rlists = rotor_system.rlists[rotor_system.index]
    orientations = rotor_system.orientations

    # operating conditions
    vinflows = rotorinflow.(Ref(freestream), orientations)
    operatingpoints_list = rotor2operatingpoints.(vinflows, rlists, omegas, Ref(environment))

    # iterate over rotors
    for (i, rotorindex) in enumerate(rotor_system.index)
        # isolate rotor
        rotor = rotor_system.rotors[rotorindex]
        sections = rotor_system.sectionlists[rotorindex]

        # isolate operating conditions
        operatingpoints = operatingpoints_list[i]

        # solve CCBlade
        outputs = CC.solve.(Ref(rotor), sections, operatingpoints)
        T, Q = CC.thrusttorque(rotor, sections, outputs)
        Ts[i] = T
        Qs[i] = Q
        Ps[i] = Q * omegas[i]
    end

    # calculate Js
    Js .= vinflows ./ omegas * 2 * pi ./ [2 * rotor.Rtip for rotor in rotor_system.rotors[rotor_system.index]]

    return Js, Ts, Qs, Ps
end

solverotor(rotor, sections, operatingpoints) = CC.solve.(Ref(rotor), sections, operatingpoints)

"""
    solverotors_outputs(rotor_system, omegas, freestream, environment)

# Arguments
* `rotor_system`
* `omegas`
* `freestream`
* `environment`

# Returns

* `Js`
* `outputs`

"""
function solverotors_outputs(rotor_system, omegas, freestream, environment)

    # extract rotor info
    rlists = rotor_system.rlists[rotor_system.index]
    orientations = rotor_system.orientations

    # operating conditions
    vinflows = rotorinflow.(Ref(freestream), orientations)
    operatingpoints_list = rotor2operatingpoints.(vinflows, rlists, omegas, Ref(environment))

    # solve rotors
    outputs_list = solverotor.(rotor_system.rotors[rotor_system.index], rotor_system.sectionlists[rotor_system.index], operatingpoints_list)

    # calculate Js
    Js = vinflows ./ omegas * 2 * pi ./ [2 * rotor.Rtip for rotor in rotor_system.rotors[rotor_system.index]]

    return Js, outputs
end

"""
    induced2wakefunction(rotor_system, us, vs;
        wake_shape_functions = Function[(Rtip, x) -> Rtip for i in 1:length(rotor_system.index)],
        axial_interpolations = Function[(rs, us, r, Rtip) -> FM.linear(rs, us, r) for i in 1:length(rotor_system.index)],
        swirl_interpolations = Function[(rs, vs, r, Rtip) -> FM.linear(rs, vs, r) for i in 1:length(rotor_system.index)],
        axial_multipliers = Function[(distance2plane, Rtip) -> 2 for i in 1:length(aircraft.rotor_system.index)],
        swirl_multipliers = Function[(distance2plane, Rtip) -> 1 for i in 1:length(aircraft.rotor_system.index)]
    )

Interpolates rotor-induced velocities into a wake function of X.

# Arguments:

* `rotor_system::CCBladeSystem`: system of rotors
* `us::Vector{Vector{Float64}}`: [i][j]th element is the axial induced velocity at the jth radial station of the ith rotor
* `vs::Vector{Vector{Float64}}`: [i][j]th element is the swirl induced velocity at the jth radial station of the ith rotor

# Keyword Arguments:

* `wakeshapefunction = (Rtip, x) -> Rtip,` : function accepts a rotor radius and axial coordinate and returns the local wake radius
* `axial_interpolations = (rs, us, r) -> FM.linear(rs, us, r)` : function accepts radial stations, axial velocities, and the radial coordinate and returns the unmodified axial induced velocity
* `swirl_interpolations = (rs, vs, r) -> FM.linear(rs, vs, r)` : function accepts radial stations, swirl velocities, and the radial coordinate and returns the unmodified swirl induced velocity
* `axialmultiplier = (u, distance2plane, Rtip) -> u * 2` : function multiplied by the axial interpolation function
* `swirlmultiplier = (v, distance2plane, Rtip) -> v` : function multiplied by the swirlinterpolation function

"""
function induced2wakefunction(rotor_system, us, vs;
    wake_shape_functions = Function[(Rtip, x) -> Rtip for i in 1:length(rotor_system.index)],
    axial_interpolations = Function[(rs, us, r, Rtip) -> FM.linear(rs, us, r) for i in 1:length(rotor_system.index)],
    swirl_interpolations = Function[(rs, vs, r, Rtip) -> FM.linear(rs, vs, r) for i in 1:length(rotor_system.index)],
    axial_multipliers = Function[(distance2plane, Rtip) -> 2 for i in 1:length(aircraft.rotor_system.index)],
    swirl_multipliers = Function[(distance2plane, Rtip) -> 1 for i in 1:length(aircraft.rotor_system.index)]
)
    function wakefunction(X)
        Vwake = zeros(3)
        for (irotor, rotorindex) in enumerate(rotor_system.index)

            # translate origin to rotor center (R is the vector from the rotor center to X)
            R = X - rotor_system.positions[irotor]
            r = LA.norm(R)

            # calculate distance from rotor plane (positive is in the negative thrust direction)
            distance2plane = -LA.dot(R, rotor_system.orientations[irotor])

            # angle between normal and R (\in [0,\pi])
            θ = acos(distance2plane/r)

            # calculate distance to centerline
            distance2centerline = r * sin(θ) # always positive

            # check if X is within the wake
            Rtip = rotor_system.rotors[rotorindex].Rtip
            Rhub = rotor_system.rotors[rotorindex].Rhub
            if distance2centerline < wake_shape_functions[irotor](Rtip, distance2plane) && distance2centerline > Rhub #! Quick fix, not as comprehensive. Do we need another function? Or perhaps the wakeshapefunction outputs a lower and upper bound on the radius
                # interpolate wake velocity
                u = axial_interpolations[irotor](rotor_system.rlists[rotorindex], us[irotor], distance2centerline, Rtip) * axial_multipliers[irotor](distance2plane, Rtip)
                v = swirl_interpolations[irotor](rotor_system.rlists[rotorindex], vs[irotor], distance2centerline, Rtip) * swirl_multipliers[irotor](distance2plane, Rtip)
                # get axial unit vector
                axialhat = -rotor_system.orientations[irotor]

                # get tangential unit vector
                tangentialhat = LA.cross(rotor_system.orientations[irotor], R)
                tangentialhat ./= LA.norm(tangentialhat)

                # project R onto the plane parallel to normal
                Vinduced = u * axialhat + v * tangentialhat * (rotor_system.spin_directions[irotor] ? 1 : -1)
                Vwake .+= Vinduced
            end
        end

        return Vwake
    end
end

"""
    `get_omega(;Vinf=nothing, D=nothing, J=nothing)`

Returns rotational velocity in [rad/s] given velocity and advance ratio.
"""
function get_omega(; Vinf=nothing, D=nothing, J=nothing)

    @assert Vinf != nothing "User must supply a velocity with the `Vinf` keyword argument."
    @assert D != nothing "User must supply a diameter with the `D` keyword argument."
    @assert J != nothing "User must supply an advance ratio with the `J` keyword argument."

    return 2pi * Vinf / (D * J)
end

"""
    `get_J(;Vinf=nothing, D=nothing, omega=nothing, n=nothing)`

Returns advance ratio given velocity and rotational velocity. Must supply either `omega` or `n`.
"""
function get_J(;Vinf=nothing, D=nothing, omega=nothing, n=nothing)

    @assert Vinf != nothing "User must supply a velocity with the `Vinf` keyword argument."
    @assert D != nothing "User must supply a diameter with the `D` keyword argument."

    if n == nothing
        @assert omega != nothing "User must supply either `omega` or `n` as a keyword argument."
        n = omega / (2pi)
    end

    return Vinf / (n * D)
end

"""
    `get_Vinf(;D=nothing, J=nothing, omega=nothing, n=nothing)`

Returns velocity in [m/s] given advance ratio and rotational velocity. Must supply either `omega` or `n`.
"""
function get_Vinf(;D=nothing, J=nothing, omega=nothing, n=nothing)

    @assert D != nothing "User must supply a diameter with the `D` keyword argument."
    @assert J != nothing "User must supply an advance ratio with the `J` keyword argument."

    if n == nothing
        @assert omega != nothing "User must supply either `omega` or `n` as a keyword argument."
        n = omega / (2pi)
    end

    return n * D / J
end

"""
    omega_to_rpm(omega)

Convenience function to translate rotational velocity to rotations per minute.

# Arguments

* `omega`: rotational velocty (rad/sec)

# Returns

* `rpm`: rotations per minute

"""
function omega_to_rpm(omega)
    return omega * 30/pi
end

"""
    rpm_to_omega(rpm)

Convenience function to translate rotations per minute to rotational velocity.

# Arguments

* `rpm`: rotations per minute

# Returns

* `omega`: rotational velocty (rad/sec)

"""
function rpm_to_omega(rpm)
    return rpm * pi/30
end

ext: jl

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
using Petri
using LabelledArrays
using OrdinaryDiffEq
using StochasticDiffEq
using DiffEqJump
using Random
using Plots

sir = Petri.Model([:S,:I,:R],LVector(
                                inf=(LVector(S=1,I=1), LVector(I=2)),
                                rec=(LVector(I=1),     LVector(R=1))))

Graph(sir)

tmax = 40.0
tspan = (0.0,tmax);

u0 = LVector(S=990.0, I=10.0, R=0.0)

p = LVector(inf=0.5/sum(u0), rec=0.25);

Random.seed!(1234);

prob_ode = ODEProblem(sir,u0,tspan,p)
sol_ode = solve(prob_ode, Tsit5());
plot(sol_ode)

prob_sde,cb = SDEProblem(sir,u0,tspan,p)
sol_sde = solve(prob_sde,LambaEM(),callback=cb);
plot(sol_sde)

prob_jump = JumpProblem(sir, u0, tspan, p)
sol_jump = solve(prob_jump,SSAStepper());
plot(sol_jump)

include(joinpath(@__DIR__,"tutorials","appendix.jl"))
appendix()

max_stars_repo_name: snowdj/sir-julia

max_issues_repo_name: snowdj/sir-julia

max_forks_repo_name: snowdj/sir-julia


text:
using WaterLily
using LinearAlgebra: norm2
include("ThreeD_Plots.jl")

function donut(p=6,Re=1e3)
    # Define simulation size, geometry dimensions, viscosity
    n = 2^p
    center,R,r = [n/2,n/2,n/2], n/4, n/16
    ν = R/Re
    @show R,ν

    # Apply signed distance function for a torus
    body = AutoBody() do xyz,t
        x,y,z = xyz - center
        norm2([x,norm2([y,z])-R])-r
    end

    # Initialize simulation
    Simulation((2n+2,n+2,n+2),[1.,0.,0.],R;ν,body),center
end

function ω_θ_data(sim)
    @inside sim.flow.σ[I] = WaterLily.ω_θ(I,[1,0,0],center,sim.flow.u)*sim.L/sim.U
    @view sim.flow.σ[2:end-1,2:end-1,2:end-1]
end
function body_data(sim)
    @inside sim.flow.σ[I] = sum(sim.flow.μ₀[I,i]+sim.flow.μ₀[I+δ(i,I),i] for i=1:3)
    @view sim.flow.σ[2:end-1,2:end-1,2:end-1]
end

sim,center = donut()
sim,fig = contour_video!(sim,ω_θ_data,body_data,name="donut.mp4",duration=10)


text:
module mrotationmatrixs2
using FinEtools
using LinearAlgebra: norm, I
# using BenchmarkTools
using Test
function rotmat3original!(Rmout, a)
    na = norm(a);
    thetatilde = zeros(3,3);
    skewmat!(thetatilde, a);
    a = a/na;
    ca = cos(na);
    sa = sin(na);
    aa = a*a';
    copyto!(Rmout, ca * (1.0*I-aa) + sa/na*thetatilde + aa);
    return Rmout
end
function test()
    for i in 1:10
        Rmout = rand(3, 3)
        a = rand(3)
        # Rmout = [0.056575544213388396 0.8235145288079604 0.4566335160952417; 0.6335215149040396 0.8852762577618685 0.803417762628938; 0.5940995203626245 0.35601280936101065 0.6825575362008556]
        # a = [0.11011521831193871, 0.5097478695998647, 0.8139760429477749]
        rotmat3!(Rmout, a)
        Rmout2 = rand(3, 3)
        rotmat3original!(Rmout2, a)
        # Rmout = [0.5736166081143759 -0.6670428981551146 0.47541325067375245; 0.7189364978350996 0.6881251218379588 0.0980516638109532; -0.3925484670406451
        # 0.2855478746485778 0.8742814834523946]
        @test norm(Rmout-Rmout2) <= 1.0e-6
        # @btime rotmat3!($Rmout, $a)
    end
    true
end
end
using .mrotationmatrixs2
mrotationmatrixs2.test()


text:
function SeisSmoothGathers(m::String,d::String,adj;Nsmooth=3,Nrepeat=1)

	if (adj==false)
		m1,h,ext = SeisRead(m)
		ngather = size(m1,3)
		for igather = 1 : ngather
			A = m1[:,:,igather];
			A = smooth_angles(A;Nsmooth=Nsmooth,Nrepeat=Nrepeat)
			m1[:,:,igather] = A
		end
		SeisWrite(d,m1,h,ext)
	else
		d1,h,ext = SeisRead(d)
		ngather = size(d1,3)
		for igather = 1 : ngather
			A = d1[:,:,igather];
			A = smooth_angles(A;Nsmooth=Nsmooth,Nrepeat=Nrepeat)
			d1[:,:,igather] = A
		end
		SeisWrite(m,d1,h,ext)
	end

end

function smooth_angles(A;Nsmooth=3,Nrepeat=1)

	for iz = 1 : size(A,1)
		a = A[iz,:]
		for irepeat = 1 : Nrepeat
			a = mean_filter(a,Nsmooth)
		end
		A[iz,:] = a
	end

	return A;
end

function mean_filter(a,nw)

	b = vec(zeros(length(a),1))
	for ix = 1 : length(a)
		sum  = 0.
		nsum = 0.
		for iw = 1 : nw
			index1 = ix - Int(floor(nw/2)) - 1 + iw
			if (index1>0 && index1<length(a))
				sum += a[index1]
				nsum += 1.
			end
		end
		b[ix] = sum/nsum
	end
	return b
end


text:
using PkgServerCacheWarmer
using Documenter

DocMeta.setdocmeta!(PkgServerCacheWarmer, :DocTestSetup, :(using PkgServerCacheWarmer); recursive=true)

makedocs(;
    modules=[PkgServerCacheWarmer],
    authors="Dilum Aluthge, Brown Center for Biomedical Informatics, and contributors",
    repo="https://github.com/bcbi/PkgServerCacheWarmer.jl/blob/{commit}{path}#{line}",
    sitename="PkgServerCacheWarmer.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://bcbi.github.io/PkgServerCacheWarmer.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
    strict=true,
)

deploydocs(;
    repo="github.com/bcbi/PkgServerCacheWarmer.jl",
    devbranch="main",
)


text:
# TO DO:
# test termination criterion

using Base.Test
using FactCheck

function qp_test()
    Q = speye(3);
    A = sparse(ones(1,3));
    b = 3*ones(1);
    c = ones(3);
    return class_quadratic_program(A, b, c, Q), A, b, c, Q;
end

settings = class_settings();

facts("core") do
  ###############################
  # test class_variables
  ###############################
  context("variables") do
    var = class_variables(3,3);

    x_values = [-1, 2, 2];
    var._v[var._x_ind] = x_values;

    @fact x(var) --> x_values
    @fact_throws ErrorException check_positive(var)

    x_values = [1, NaN, 2];
    var._v[var._x_ind] = x_values;
    @fact_throws ErrorException check_for_wrong_vals(var)

    x_values = [1, Inf, 2];
    var._v[var._x_ind] = x_values;
    @fact_throws ErrorException check_for_wrong_vals(var)

    x_values = [1, -Inf, 2];
    var._v[var._x_ind] = x_values;
    @fact_throws ErrorException check_for_wrong_vals(var)

    var = class_variables(3,3);
    @fact_throws ErrorException move!(var,-2.0,var)
  end

  qp, A, b, c, Q = qp_test();

  # test quadratic program and validate dimensions
  context("qp and validate dimensions") do
      x_var = [1.0,1.0,1.0]
      y_var = [1.0];

      @fact 4.5 --> internal_eval_c(qp, x_var)
      @fact [0.0] --> internal_eval_a(qp, x_var)
      @fact A --> internal_eval_jac_a(qp, x_var)
      @fact Q --> internal_eval_hesslag_prod(qp, x_var, y_var)
      @fact c + Q * x_var - A' * y_var --> internal_eval_gradlag(qp, x_var, y_var)

      vars2 = class_variables(3,3);
      @fact_throws ErrorException validate_dimensions(qp, vars2)
  end
end

begin
    start_advanced_timer()
    qp, A, b, c, Q = qp_test();
    nlp_vals = class_nlp_cache();
    vars = class_variables(3,1);
    update_nlp_cache!(nlp_vals, qp, vars)
    validate_dimensions(qp, vars)

    newt = class_homogeneous_newton();
    # test newton direction
    facts("newton direction") do

      begin
          newt = class_homogeneous_newton();
          initialize_newton!(newt, qp, vars, settings)
          res = class_homogeneous_residuals();
      end

      begin
          @fact 1.0 --> mu(newt, vars);

          validate_dimensions(qp, vars)
          initialize_newton!(newt, qp, vars, settings)
          update_newton!(newt, vars, settings)
          update_newton_diag!(newt, vars, settings)
          validate_dimensions( newt.nlp_vals )
          form_woodbury!(newt, vars)
          compute_newton_direction!(newt, vars, class_theta(0.5,0.5,0.5))
      end

  end

  facts("strategies") do
      # test stratagies

      # test delta heuristics
      begin
          #println("TEST IPOPT INERTIA CORRECTION")
          validate_dimensions(qp, vars)
          delta, number_of_factors = ipopt_style_inertia_correction!(newt, vars, settings)
          @fact delta --> 0.0
          @fact number_of_factors --> 1
      end
  end

  pause_advanced_timer()

    begin
        vars = class_variables(3,1);
        A = sparse(ones(1,3));
        b = 3*ones(1);
        c = [2.0, 0.0, 1.0];
        qp = class_quadratic_program(A, b, c);
        status, vars = homogeneous_algorithm(qp, vars, settings)
        @fact norm(x_scaled(vars) - [0.0,3.0,0.0]) --> less_than(1e-6)
        @fact tau(vars) --> greater_than(0.1)
        @fact 0.0 --> less_than(kappa(vars))
        @fact kappa(vars) --> less_than(1e-6)
    end

end

include("sample_nlps.jl")

max_forks_repo_name: ohinder/Everest


text:
# TODO require that T <: U?
struct ShapedView{T,N,SH,D<:AbstractVector{T}} <: DenseArray{T,N}
    data::D
    offset::Int
    Base.@propagate_inbounds function ShapedView{T,N,SH,D}(
        data::D,
        offset::Int
    ) where {T,N,SH,D<:AbstractVector{T}}
        @boundscheck check_offset_shape_inbounds(data, offset, SH)
        check_has_unit_axes(data)
        SH isa AbstractShape || throw(ArgumentError("Type parameter SH must be <: AbstractShape"))
        new{T,N,SH,D}(data, offset)
    end
end

@propagate_inbounds function ShapedView{T,N,SH}(data, offset) where {T,N,SH}
    ShapedView{T,N,SH,typeof(data)}(data, offset)
end

@propagate_inbounds function ShapedView(data::AbstractVector, offset::Int, shape::AbstractShape)
    ShapedView{eltype(data), ndims(shape), shape, typeof(data)}(data, offset)
end

@propagate_inbounds function ShapedView(data::AbstractVector, shape::AbstractShape)
    ShapedView{eltype(data), ndims(shape), shape, typeof(data)}(data, 0)
end

@propagate_inbounds ShapedView(data, shape) = ShapedView(data, 0, shape)

@propagate_inbounds (shape::AbstractShape)(data) = ShapedView(data, 0, shape)

@pure shapeof(::Type{SV}) where {T,N,SH,SV <: ShapedView{T,N,SH}} = SH
@inline shapeof(A::ShapedView) = shapeof(typeof(A))

@pure Size(::Type{SV}) where {SV <: ShapedView} = Size(shapeof(SV))
@inline Size(A::ShapedView) = Size(typeof(A))

@pure Length(::Type{SV}) where {SV <: ShapedView} = Length(shapeof(SV))
@inline Length(A::ShapedView) = Length(typeof(A))

@pure Base.size(::Type{SV}) where {SV <: ShapedView} = size(shapeof(SV))
@inline Base.size(A::ShapedView) = size(typeof(A))

@pure Base.length(::Type{SV}) where {SV <: ShapedView} = length(shapeof(SV))
@inline Base.length(A::ShapedView) = length(typeof(A))

@pure Base.ndims(::Type{SV}) where {SV <: ShapedView} = ndims(shapeof(SV))
@inline Base.ndims(A::ShapedView) = ndims(typeof(A))

@pure Base.axes(::Type{SV}) where {SV <: ShapedView} = axes(shapeof(SV))
@inline Base.axes(A::ShapedView) = axes(typeof(A))

@pure Base.IndexStyle(::Type{SV}) where {SV <: ShapedView} = IndexLinear()
@inline Base.IndexStyle(A::ShapedView) = IndexStyle(typeof(A))

@pure Base.propertynames(::Type{SV}) where {SV <: ShapedView} = propertynames(shapeof(SV))
@inline Base.propertynames(A::ShapedView) = propertynames(typeof(A))

@inline function Base.unsafe_convert(::Type{Ptr{T}}, A::ShapedView{T}) where {T}
    Base.unsafe_convert(Ptr{T}, _data(A))
end

@inline Base.dataids(A::ShapedView) = Base.dataids(_data(A))

@inline function Base.copy(A::ShapedView{T,N,SH}) where {T,N,SH}
    @inbounds ShapedView{T,N,SH}(copy(_data(A)), _offset(A))
end

@propagate_inbounds function Base.getproperty(A::ShapedView, name::Symbol)
    shape = shapeof(A)
    offset = _offset(A) + getoffset(shape, name)
    innershape = getproperty(shape, name)
    _maybe_shapedview(_data(A), offset, innershape)
end

@propagate_inbounds function _maybe_shapedview(data, offset, ::AbstractScalarShape)
    i = offset + firstindex(data)
    @boundscheck checkbounds(data, i)
    @inbounds getindex(data, i)
end

@propagate_inbounds function _maybe_shapedview(data, offset, shape)
    @boundscheck check_offset_shape_inbounds(data, offset, shape)
    @inbounds ShapedView{eltype(data), ndims(shape), shape, typeof(data)}(data, offset)
end

@propagate_inbounds function Base.setproperty!(A::ShapedView, name::Symbol, x)
    innershape = getproperty(shapeof(A), name)
    _maybe_setproperty!(innershape, A, name, x)
end

@propagate_inbounds function _maybe_setproperty!(shape::AbstractScalarShape, A, name, x)
    offset = _offset(A) + getoffset(shapeof(A), name)
    data = _data(A)
    i = firstindex(data) + offset
    @boundscheck checkbounds(data, i)
    @inbounds setindex!(data, x, i)
end

@inline function _maybe_setproperty!(shape::AbstractShape, A, name, x)
    error("Cannot call `setproperty!` for shape $name of type $(typeof(shape))")
end

@propagate_inbounds function Base.getindex(A::ShapedView, i::Int)
    data = _data(A)
    i += _offset(A)
    @boundscheck checkbounds(data, i)
    @inbounds getindex(data, i)
end

@propagate_inbounds function Base.setindex!(A::ShapedView, val, i::Int)
    data = _data(A)
    i += _offset(A)
    @boundscheck checkbounds(data, i)
    @inbounds setindex!(data, val, i)
end

function UnsafeArrays.unsafe_uview(A::ShapedView{T,N,SH}) where {T,N,SH}
    @inbounds ShapedView{T,N,SH}(UnsafeArrays.unsafe_uview(_data(A)), _offset(A))
end

function check_offset_shape_inbounds(data::AbstractVector{T}, offset::Int, shape::AbstractShape{S,U,N}) where {T,S,U,N}
    if length(shape) > 0 && !(0 <= offset < length(data))
        throw(ArgumentError("offset must be in range [0, length(data))"))
    end
    if length(data) < offset + length(shape)
        throw(ArgumentError("offset + length(shape) cannot be greater than length(data)"))
    end
end

@inline _data(A::ShapedView) = getfield(A, :data)
@inline _offset(A::ShapedView) = getfield(A, :offset)


text:
# Automatically generated using Clang.jl

const CONEMODE = 0
const INSIDE_CONE = 0
const OUTSIDE_CONE = 1
const idxint = @static Sys.iswindows() ? Int64 : Clong
const pfloat = Cdouble

struct lpcone
    p::idxint
    w::Ptr{pfloat}
    v::Ptr{pfloat}
    kkt_idx::Ptr{idxint}
end

struct socone
    p::idxint
    skbar::Ptr{pfloat}
    zkbar::Ptr{pfloat}
    a::pfloat
    d1::pfloat
    w::pfloat
    eta::pfloat
    eta_square::pfloat
    q::Ptr{pfloat}
    Didx::Ptr{idxint}
    u0::pfloat
    u1::pfloat
    v1::pfloat
end

struct expcone
    colstart::NTuple{3,idxint}
    v::NTuple{6,pfloat}
    g::NTuple{3,pfloat}
end

struct cone
    lpc::Ptr{lpcone}
    soc::Ptr{socone}
    nsoc::idxint
    expc::Ptr{expcone}
    nexc::idxint
    fexv::idxint
end

# Skipping MacroDefinition: init_ctrlc ( )
# Skipping MacroDefinition: remove_ctrlc ( )
# Skipping MacroDefinition: check_ctrlc ( ) 0

const ECOS_VERSION = "2.0.8"
const MAXIT = 100
const FEASTOL = 1.0e-8
const ABSTOL = 1.0e-8
const RELTOL = 1.0e-8
const FTOL_INACC = 0.0001
const ATOL_INACC = 5.0e-5
const RTOL_INACC = 5.0e-5
const GAMMA = 0.99
const STATICREG = 1
const DELTASTAT = 7.0e-8
const DELTA = 2.0e-7
const EPS = 1.0e-13
const VERBOSE = 1
const NITREF = 9
const IRERRFACT = 6
const LINSYSACC = 1.0e-14
const SIGMAMIN = 0.0001
const SIGMAMAX = 1.0
const STEPMIN = 1.0e-6
const STEPMAX = 0.999
const SAFEGUARD = 500
const MAX_BK = 90
const BK_SCALE = 0.8
const MIN_DISTANCE = 0.1
const CENTRALITY = 1
const EQUILIBRATE = 1
const EQUIL_ITERS = 3
const ECOS_OPTIMAL = 0
const ECOS_PINF = 1
const ECOS_DINF = 2
const ECOS_INACC_OFFSET = 10
const ECOS_MAXIT = -1
const ECOS_NUMERICS = -2
const ECOS_OUTCONE = -3
const ECOS_SIGINT = -4
const ECOS_FATAL = -7

# Skipping MacroDefinition: MAX ( X , Y ) ( ( X ) < ( Y ) ? ( Y ) : ( X ) )
# Skipping MacroDefinition: SAFEDIV_POS ( X , Y ) ( ( Y ) < EPS ? ( ( X ) / EPS ) : ( X ) / ( Y ) )

struct settings
    gamma::pfloat
    delta::pfloat
    eps::pfloat
    feastol::pfloat
    abstol::pfloat
    reltol::pfloat
    feastol_inacc::pfloat
    abstol_inacc::pfloat
    reltol_inacc::pfloat
    nitref::idxint
    maxit::idxint
    verbose::idxint
    max_bk_iter::idxint
    bk_scale::pfloat
    centrality::pfloat
end

struct stats
    pcost::pfloat
    dcost::pfloat
    pres::pfloat
    dres::pfloat
    pinf::pfloat
    dinf::pfloat
    pinfres::pfloat
    dinfres::pfloat
    gap::pfloat
    relgap::pfloat
    sigma::pfloat
    mu::pfloat
    step::pfloat
    step_aff::pfloat
    kapovert::pfloat
    iter::idxint
    nitref1::idxint
    nitref2::idxint
    nitref3::idxint
    tsetup::pfloat
    tsolve::pfloat
    pob::idxint
    cb::idxint
    cob::idxint
    pb::idxint
    db::idxint
    affBack::idxint
    cmbBack::idxint
    centrality::pfloat
end

struct spmat
    jc::Ptr{idxint}
    ir::Ptr{idxint}
    pr::Ptr{pfloat}
    n::idxint
    m::idxint
    nnz::idxint
end

struct kkt
    PKPt::Ptr{spmat}
    L::Ptr{spmat}
    D::Ptr{pfloat}
    work1::Ptr{pfloat}
    work2::Ptr{pfloat}
    work3::Ptr{pfloat}
    work4::Ptr{pfloat}
    work5::Ptr{pfloat}
    work6::Ptr{pfloat}
    RHS1::Ptr{pfloat}
    RHS2::Ptr{pfloat}
    dx1::Ptr{pfloat}
    dx2::Ptr{pfloat}
    dy1::Ptr{pfloat}
    dy2::Ptr{pfloat}
    dz1::Ptr{pfloat}
    dz2::Ptr{pfloat}
    P::Ptr{idxint}
    Pinv::Ptr{idxint}
    PK::Ptr{idxint}
    Parent::Ptr{idxint}
    Sign::Ptr{idxint}
    Pattern::Ptr{idxint}
    Flag::Ptr{idxint}
    Lnz::Ptr{idxint}
    delta::pfloat
end

struct pwork
    n::idxint
    m::idxint
    p::idxint
    D::idxint
    x::Ptr{pfloat}
    y::Ptr{pfloat}
    z::Ptr{pfloat}
    s::Ptr{pfloat}
    lambda::Ptr{pfloat}
    kap::pfloat
    tau::pfloat
    best_x::Ptr{pfloat}
    best_y::Ptr{pfloat}
    best_z::Ptr{pfloat}
    best_s::Ptr{pfloat}
    best_kap::pfloat
    best_tau::pfloat
    best_cx::pfloat
    best_by::pfloat
    best_hz::pfloat
    best_info::Ptr{stats}
    dsaff::Ptr{pfloat}
    dzaff::Ptr{pfloat}
    W_times_dzaff::Ptr{pfloat}
    dsaff_by_W::Ptr{pfloat}
    saff::Ptr{pfloat}
    zaff::Ptr{pfloat}
    C::Ptr{cone}
    A::Ptr{spmat}
    G::Ptr{spmat}
    c::Ptr{pfloat}
    b::Ptr{pfloat}
    h::Ptr{pfloat}
    AtoK::Ptr{idxint}
    GtoK::Ptr{idxint}
    xequil::Ptr{pfloat}
    Aequil::Ptr{pfloat}
    Gequil::Ptr{pfloat}
    resx0::pfloat
    resy0::pfloat
    resz0::pfloat
    rx::Ptr{pfloat}
    ry::Ptr{pfloat}
    rz::Ptr{pfloat}
    rt::pfloat
    hresx::pfloat
    hresy::pfloat
    hresz::pfloat
    nx::pfloat
    ny::pfloat
    nz::pfloat
    ns::pfloat
    cx::pfloat
    by::pfloat
    hz::pfloat
    sz::pfloat
    KKT::Ptr{kkt}
    info::Ptr{stats}
    stgs::Ptr{settings}
end

const MI_PRINTLEVEL = 1
const MI_ABS_EPS = 1.0e-6
const MI_REL_EPS = 0.001
const MI_MAXITER = 1000
const MI_INT_TOL = FTOL_INACC
const MI_SOLVED_NON_BRANCHABLE = 3
const MI_SOLVED_BRANCHABLE = 2
const MI_NOT_SOLVED = 1
const MI_FREE = 0
const MI_ONE = 1
const MI_ZERO = 0
const MI_STAR = -1
const MI_OPTIMAL_SOLN = ECOS_OPTIMAL
const MI_INFEASIBLE = ECOS_PINF
const MI_UNBOUNDED = ECOS_DINF
const MI_MAXITER_FEASIBLE_SOLN = ECOS_OPTIMAL + ECOS_INACC_OFFSET
const MI_MAXITER_NO_SOLN = ECOS_PINF + ECOS_INACC_OFFSET
const MI_MAXITER_UNBOUNDED = ECOS_DINF + ECOS_INACC_OFFSET
const MAX_FLOAT_INT = 8388608

@cenum BRANCHING_STRATEGY::UInt32 begin
    BRANCHING_STRATEGY_MOST_INFEASIBLE = 0
    BRANCHING_STRATEGY_STRONG_BRANCHING = 1
    BRANCHING_STRATEGY_PSEUDOCOST_BRANCHING = 2
    BRANCHING_STRATEGY_RELIABILITY = 3
    BRANCHING_STRATEGY_RANDOM = 4
end

@cenum NODE_SELECTION_METHOD::UInt32 begin
    BREADTH_FIRST = 0
    DIVE_LOWER_NODE = 1
    DIVE_UPPER_NODE = 2
end

struct settings_bb
    maxit::idxint
    verbose::idxint
    abs_tol_gap::pfloat
    rel_tol_gap::pfloat
    integer_tol::pfloat
    branching_strategy::BRANCHING_STRATEGY
    reliable_eta::idxint
    node_selection_method::NODE_SELECTION_METHOD
end

struct node
    status::UInt8
    L::pfloat
    U::pfloat
    relaxation::pfloat
    split_idx::idxint
    split_val::pfloat
    prev_split_idx::idxint
    prev_split_val::pfloat
    prev_relaxation::pfloat
    up_branch_node::Cint
end

struct ecos_bb_pwork
    num_bool_vars::idxint
    num_int_vars::idxint
    nodes::Ptr{node}
    bool_node_ids::Cstring
    int_node_ids::Ptr{pfloat}
    bool_vars_idx::Ptr{idxint}
    int_vars_idx::Ptr{idxint}
    ecos_prob::Ptr{pwork}
    A::Ptr{spmat}
    G::Ptr{spmat}
    c::Ptr{pfloat}
    b::Ptr{pfloat}
    h::Ptr{pfloat}
    x::Ptr{pfloat}
    y::Ptr{pfloat}
    z::Ptr{pfloat}
    s::Ptr{pfloat}
    kap::pfloat
    tau::pfloat
    info::Ptr{stats}
    global_U::pfloat
    global_L::pfloat
    tmp_bool_node_id::Cstring
    tmp_int_node_id::Ptr{pfloat}
    iter::idxint
    dive_node_id::idxint
    tmp_branching_bool_node_id::Cstring
    tmp_branching_int_node_id::Ptr{pfloat}
    pseudocost_bin_sum::Ptr{pfloat}
    pseudocost_int_sum::Ptr{pfloat}
    pseudocost_bin_cnt::Ptr{idxint}
    pseudocost_int_cnt::Ptr{idxint}
    Gpr_new::Ptr{pfloat}
    Gjc_new::Ptr{idxint}
    Gir_new::Ptr{idxint}
    h_new::Ptr{pfloat}
    ecos_stgs::Ptr{settings}
    stgs::Ptr{settings_bb}
    default_settings::idxint
end

const PRINTLEVEL = 2
const PROFILING = 1
const DEBUG = 0
const ECOS_INFINITY = Inf
const ECOS_NAN = NaN
# const PRINTTEXT = printf
# const MALLOC = malloc
# const FREE = free
const KKT_PROBLEM = 0
const KKT_OK = 1


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.3
#     language: julia
#     name: julia-1.6
# ---

# %%
using SparseArrays

nzptr(A::SparseArrays.AbstractSparseMatrixCSC, col::Integer) = view(rowvals(A), nzrange(A, col))

using Random, LinearAlgebra
Random.seed!(12345678)

n = 20000
p = 2/n
@time A = Symmetric(sprand(Bool, n, n, p)) # test adjacency matrix
@time B = sparse(A)
@time C = nzptr.(Ref(B), axes(B, 2)) # adjacency list

# %%
B[:, 4]

# %%
Random.seed!(12345678)

n = 20000
p = 0.1
@time A = Symmetric(sprand(Bool, n, n, p)) # test adjacency matrix
@time B = sparse(A)
@time C = nzptr.(Ref(B), axes(B, 2)) # adjacency list

# %%
B[:, 1]

# %%

max_stars_repo_path: 0022/adj matrix to adj list.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0022/adj matrix to adj list.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0022/adj matrix to adj list.jl

max_forks_repo_name: genkuroki/public


text:
function GPSAlmanac(arg0::JString, arg1::jint, arg10::jdouble, arg11::jdouble, arg12::jdouble, arg13::jdouble, arg14::jint, arg15::jint, arg16::jint, arg17::AbsoluteDate, arg2::jint, arg3::jint, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return GPSAlmanac((JString, jint, jdouble, jdouble, jdouble, jdouble, jint, jint, jint, AbsoluteDate, jint, jint, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg11, arg12, arg13, arg14, arg15, arg16, arg17, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function GPSAlmanac(arg0::JString, arg1::jint, arg10::jdouble, arg11::jdouble, arg12::jdouble, arg13::jdouble, arg14::jint, arg15::jint, arg16::jint, arg2::jint, arg3::jint, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return GPSAlmanac((JString, jint, jdouble, jdouble, jdouble, jdouble, jint, jint, jint, jint, jint, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg11, arg12, arg13, arg14, arg15, arg16, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function get_af0(obj::GPSAlmanac)
    return jcall(obj, "getAf0", jdouble, ())
end

function get_af1(obj::GPSAlmanac)
    return jcall(obj, "getAf1", jdouble, ())
end

function get_cic(obj::GPSAlmanac)
    return jcall(obj, "getCic", jdouble, ())
end

function get_cis(obj::GPSAlmanac)
    return jcall(obj, "getCis", jdouble, ())
end

function get_crc(obj::GPSAlmanac)
    return jcall(obj, "getCrc", jdouble, ())
end

function get_crs(obj::GPSAlmanac)
    return jcall(obj, "getCrs", jdouble, ())
end

function get_cuc(obj::GPSAlmanac)
    return jcall(obj, "getCuc", jdouble, ())
end

function get_cus(obj::GPSAlmanac)
    return jcall(obj, "getCus", jdouble, ())
end

function get_date(obj::GPSAlmanac)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_e(obj::GPSAlmanac)
    return jcall(obj, "getE", jdouble, ())
end

function get_health(obj::GPSAlmanac)
    return jcall(obj, "getHealth", jint, ())
end

function get_i0(obj::GPSAlmanac)
    return jcall(obj, "getI0", jdouble, ())
end

function get_i_dot(obj::GPSAlmanac)
    return jcall(obj, "getIDot", jdouble, ())
end

function get_iodc(obj::GPSOrbitalElements)
    return jcall(obj, "getIODC", jint, ())
end

function get_iode(obj::GPSOrbitalElements)
    return jcall(obj, "getIODE", jint, ())
end

function get_m0(obj::GPSAlmanac)
    return jcall(obj, "getM0", jdouble, ())
end

function get_mean_motion(obj::GPSAlmanac)
    return jcall(obj, "getMeanMotion", jdouble, ())
end

function get_omega0(obj::GPSAlmanac)
    return jcall(obj, "getOmega0", jdouble, ())
end

function get_omega_dot(obj::GPSAlmanac)
    return jcall(obj, "getOmegaDot", jdouble, ())
end

function get_pa(obj::GPSAlmanac)
    return jcall(obj, "getPa", jdouble, ())
end

function get_prn(obj::GPSAlmanac)
    return jcall(obj, "getPRN", jint, ())
end

function get_sat_configuration(obj::GPSAlmanac)
    return jcall(obj, "getSatConfiguration", jint, ())
end

function get_sma(obj::GPSAlmanac)
    return jcall(obj, "getSma", jdouble, ())
end

function get_source(obj::GPSAlmanac)
    return jcall(obj, "getSource", JString, ())
end

function get_svn(obj::GPSAlmanac)
    return jcall(obj, "getSVN", jint, ())
end

function get_tgd(obj::GPSOrbitalElements)
    return jcall(obj, "getTGD", jdouble, ())
end

function get_time(obj::GPSAlmanac)
    return jcall(obj, "getTime", jdouble, ())
end

function get_ura(obj::GPSAlmanac)
    return jcall(obj, "getURA", jint, ())
end

function get_week(obj::GPSAlmanac)
    return jcall(obj, "getWeek", jint, ())
end


text:
using Unitful#, UnitfulPlots
using OrdinaryDiffEq, DiffEqBase
NON_IMPLICIT_ALGS = filter((x)->isleaftype(x) && !OrdinaryDiffEq.isimplicit(x()),union(subtypes(OrdinaryDiffEqAlgorithm),subtypes(OrdinaryDiffEqAdaptiveAlgorithm)))

f = (t,y) -> 0.5*y / 3.0u"s"
u0 = 1.0u"N"
prob = ODEProblem(f,u0,(0.0u"s",1.0u"s"))

sol =solve(prob,Euler(),dt=1u"s"/10)

sol =solve(prob,Midpoint(),dt=1u"s"/10)

sol =solve(prob,ExplicitRK())

for alg in NON_IMPLICIT_ALGS
  if !(alg <: DP5Threaded) && !(alg <: Rosenbrock23) && !(alg <: Rosenbrock32)
    @show alg
    sol = solve(prob,alg(),dt=1u"s"/10)
  end
end

println("Units for Number pass")

u0 = [1.0u"N" 2.0u"N"
      3.0u"N" 1.0u"N"]
f = (t,y,dy) -> (dy .= 0.5.*y ./ 3.0u"s")
prob = ODEProblem(f,u0,(0.0u"s",1.0u"s"))

sol =solve(prob,RK4(),dt=1u"s"/10)

sol =solve(prob,ExplicitRK())
sol =solve(prob,DP5())

sol =solve(prob,DP5Threaded())

for alg in NON_IMPLICIT_ALGS
  if !(alg <: Rosenbrock23) && !(alg <: Rosenbrock32)
    @show alg
    sol = solve(prob,alg(),dt=1u"s"/10)
  end
end

println("Units for 2D pass")