text:
@get_units @with_kw struct HBV{T}
    Δt::T | "s"                     # Model time step [s]
    n::Int | "-"                    # Number of cells
    fc::Vector{T} | "mm"            # Field capacity [mm]
    betaseepage::Vector{T} | "-"    # Exponent in soil runoff generation equation [-]
    lp::Vector{T} | "-"             # Fraction of field capacity below which actual evaporation=potential evaporation [-]
    threshold::Vector{T} | "mm"     # Threshold soilwater storage above which AE=PE [mm]
    k4::Vector{T} | "Δt-1"          # Recession constant baseflow [Δt⁻¹]
    kquickflow::Vector{T} | "Δt-1"  # Recession constant upper reservoir [Δt⁻¹]
    suz::Vector{T} | "mm"           # Level over which k0 is used [mm]
    k0::Vector{T} | "Δt-1"          # Recession constant upper reservoir [Δt⁻¹]
    khq::Vector{T} | "Δt-1"         # Recession rate at flow hq [Δt⁻¹]
    hq::Vector{T}                   # High flow rate hq for which recession rate of upper reservoir is known [mm Δt⁻¹]
    alphanl::Vector{T}              # Measure of non-linearity of upper reservoir
    perc::Vector{T}                 # Percolation from upper to lower zone [mm Δt⁻¹]
    cfr::Vector{T} | "-"            # Refreezing efficiency constant in refreezing of freewater in snow [-]
    pcorr::Vector{T} | "-"          # Correction factor for precipitation [-]
    rfcf::Vector{T} | "-"           # Correction factor for rainfall [-]
    sfcf::Vector{T} | "-"           # Correction factor for snowfall [-]
    cflux::Vector{T}                # Maximum capillary rise from runoff response routine to soil moisture routine [mm Δt⁻¹]
    icf::Vector{T} | "mm"           # Maximum interception storage (in forested and non-forested areas) [mm]
    cevpf::Vector{T} | "-"          # Correction factor for potential evaporation [-]
    epf::Vector{T} | "mm-1"         # Exponent of correction factor for evaporation on days with precipitation
    ecorr::Vector{T} | "-"          # Evap correction [-]
    tti::Vector{T} | "ᵒC"           # Critical temperature for snowmelt and refreezing [ᵒC]
    tt::Vector{T} | "ᵒC"            # Defines interval in which precipitation falls as rainfall and snowfall [ᵒC]
    ttm::Vector{T} | "ᵒC"           # Threshold temperature for snowmelt [ᵒC]
    cfmax::Vector{T} | "mm ᵒC-1 Δt-1"   # Meltconstant in temperature-index [-]
    whc::Vector{T} | "-"                # Fraction of snow volume that can store water [-]
    g_tt::Vector{T} | "ᵒC"              # Threshold temperature for snowfall above glacier [ᵒC]
    g_cfmax::Vector{T} | "mm ᵒC-1 Δt-1" # Degree-day factor [mm ᵒC⁻¹ Δt⁻¹] for glacier
    g_sifrac::Vector{T} | "-"           # Fraction of the snowpack on top of the glacier converted into ice [-]
    glacierstore::Vector{T} | "mm"      # Water within the glacier [mm]
    glacierfrac::Vector{T} | "-"        # Fraction covered by a glacier [-]
    precipitation::Vector{T}            # Precipitation [mm Δt⁻¹]
    temperature::Vector{T} | "ᵒC"       # Temperature [ᵒC]
    potential_evaporation::Vector{T}    # Potential evapotranspiration [mm Δt⁻¹]
    potsoilevap::Vector{T}              # Potential soil evaporation [mm Δt⁻¹]
    soilevap::Vector{T}                 # Soil evaporation [mm Δt⁻¹]
    intevap::Vector{T}                  # Evaporation from interception storage [mm Δt⁻¹]
    actevap::Vector{T}                  # Total actual evapotranspiration (intevap + soilevap) [mm Δt⁻¹]
    interceptionstorage::Vector{T} | "mm"  # Actual interception storage [mm]
    snowwater::Vector{T} | "mm"            # Available free water in snow [mm]
    snow::Vector{T} | "mm"                 # Snow pack [mm]
    rainfallplusmelt::Vector{T}            # Snow melt + precipitation as rainfall [mm Δt⁻¹]
    soilmoisture::Vector{T} | "mm"         # Actual soil moisture [mm]
    directrunoff::Vector{T}             # Direct runoff to upper zone [mm Δt⁻¹]
    hbv_seepage::Vector{T}              # Recharge to upper zone [mm Δt⁻¹]
    in_upperzone::Vector{T}             # Water inflow into upper zone [mm Δt⁻¹]
    upperzonestorage::Vector{T} | "mm"  # Water content of the upper zone [mm]
    quickflow::Vector{T}                # Specific runoff (quickflow part) [mm Δt⁻¹]
    real_quickflow::Vector{T}           # Specific runoff (quickflow), if K upper zone is precalculated [mm Δt⁻¹]
    percolation::Vector{T}              # Actual percolation to the lower zone [mm Δt⁻¹]
    capflux::Vector{T}                  # Capillary rise [mm Δt⁻¹]
    lowerzonestorage::Vector{T} | "mm"  # Water content of the lower zone [mm]
    baseflow::Vector{T}                 # Specific runoff (baseflow part) per cell [mm Δt⁻¹]
    runoff::Vector{T}                   # Total specific runoff per cell [mm Δt⁻¹]
end

statevars(::HBV) = (
    :soilmoisture,
    :snow,
    :snowwater,
    :upperzonestorage,
    :lowerzonestorage,
    :interceptionstorage,
)

function update_until_snow(hbv::HBV, config)

    for i = 1:hbv.n
        precipitation = hbv.precipitation[i] * hbv.pcorr[i]
        # fraction of precipitation which falls as rain
        rainfrac = if iszero(hbv.tti[i])
            Float64(hbv.temperature[i] > hbv.tt[i])
        else
            frac = (hbv.temperature[i] - (hbv.tt[i] - hbv.tti[i] / 2.0)) / hbv.tti[i]
            min(frac, 1.0)
        end
        rainfrac = max(rainfrac, 0.0)

        # fraction of precipitation which falls as snow
        snowfrac = 1.0 - rainfrac
        # different correction for rainfall and snowfall
        precipitation =
            hbv.sfcf[i] * snowfrac * precipitation + hbv.rfcf[i] * rainfrac * precipitation

        # interception
        interception = min(precipitation, hbv.icf[i] - hbv.interceptionstorage[i])
        # current interception storage
        interceptionstorage = hbv.interceptionstorage[i] + interception
        precipitation = precipitation - interception

        # correction for potential evaporation on wet days
        potevap =
            exp(-hbv.epf[i] * precipitation) * hbv.ecorr[i] * hbv.potential_evaporation[i]
        # correct per landuse
        potevap = hbv.cevpf[i] * potevap

        # evaporation from interception storage
        intevap = min(interceptionstorage, potevap)
        interceptionstorage = interceptionstorage - intevap
        restevap = max(0.0, potevap - intevap)

        snow, snowwater, snowmelt, rainfallplusmelt, snowfall = snowpack_hbv(
            hbv.snow[i],
            hbv.snowwater[i],
            precipitation,
            hbv.temperature[i],
            hbv.tti[i],
            hbv.tt[i],
            hbv.ttm[i],
            hbv.cfmax[i],
            hbv.whc[i],
        )

        # update the outputs and states
        hbv.rainfallplusmelt[i] = rainfallplusmelt
        hbv.snowwater[i] = snowwater
        hbv.snow[i] = snow
        hbv.interceptionstorage[i] = interceptionstorage
        hbv.potsoilevap[i] = restevap
        hbv.intevap[i] = intevap
    end
end

function update_after_snow(hbv::HBV, config)

    modelglacier = get(config.model, "glacier", false)::Bool
    set_kquickflow = get(config.model, "set_kquickflow", false)::Bool
    external_qbase = get(config.model, "external_qbase", false)::Bool

    for i = 1:hbv.n
        if modelglacier
            # Run Glacier module and add the snowpack on-top of it.
            # Estimate the fraction of snow turned into ice (HBV-light).
            # Estimate glacier melt.

            hbv.snow[i], snow2glacier, hbv.glacierstore[i], glaciermelt = glacier_hbv(
                hbv.glacierfrac[i],
                hbv.glacierstore[i],
                hbv.snow[i],
                hbv.temperature[i],
                hbv.g_tt[i],
                hbv.g_cfmax[i],
                hbv.g_sifrac[i],
                Δt,
            )
            # Convert to mm per grid cell and add to snowmelt
            glaciermelt = glaciermelt * hbv.glacierfrac[i]
            rainfallplusmelt = hbv.rainfallplusmelt[i] + glaciermelt
        else
            rainfallplusmelt = hbv.rainfallplusmelt[i]
        end

        soilmoisture = hbv.soilmoisture[i] + rainfallplusmelt
        # if soil is filled to capacity: abundant water runs of directly
        directrunoff = max(soilmoisture - hbv.fc[i], 0.0)
        soilmoisture = soilmoisture - directrunoff
        # net water which infiltrates into soil
        netinsoil = rainfallplusmelt - directrunoff

        # soil evapotranspiration
        soilevap =
            soilmoisture > hbv.threshold[i] ? min(soilmoisture, hbv.potsoilevap[i]) :
            min(hbv.potsoilevap[i] * (soilmoisture / hbv.threshold[i]))
        # evaporation from soil moisture storage
        soilmoisture = soilmoisture - soilevap
        # sum of evaporation components (IntEvap+SoilEvap)
        actevap = hbv.intevap[i] + soilevap
        # runoff water from soil
        hbv_seepage =
            pow(min(soilmoisture / hbv.fc[i], 1.0), hbv.betaseepage[i]) * netinsoil
        soilmoisture = soilmoisture - hbv_seepage
        # correction for extremely wet periods: soil is filled to capacity
        back_tosoil = min(hbv.fc[i] - soilmoisture, directrunoff)
        directrunoff = directrunoff - back_tosoil
        soilmoisture = soilmoisture + back_tosoil
        # total water available for runoff
        in_upperzone = directrunoff + hbv_seepage
        ### calculations for upper zone ###
        upperzonestorage = hbv.upperzonestorage[i] + in_upperzone
        percolation = min(hbv.perc[i], upperzonestorage - in_upperzone / 2.0)
        upperzonestorage = upperzonestorage - percolation
        # capillary flux flowing back to soil
        capflux = hbv.cflux[i] * ((hbv.fc[i] - soilmoisture) / hbv.fc[i])
        capflux = min(hbv.fc[i] - soilmoisture, capflux)
        upperzonestorage = upperzonestorage - capflux
        soilmoisture = soilmoisture + capflux

        real_quickflow = 0.0
        if set_kquickflow == false

            if percolation < hbv.perc[i]
                quickflow = 0.0
            else
                quickflow = min(
                    hbv.kquickflow[i] * pow(
                        (upperzonestorage - min(in_upperzone / 2.0, upperzonestorage)),
                        1.0 + hbv.alphanl[i],
                    ),
                    upperzonestorage,
                )

                upperzonestorage =
                    percolation < hbv.perc[i] ? upperzonestorage :
                    max(upperzonestorage - quickflow, 0.0)
            end
        else
            quickflow = hbv.kquickflow[i] * upperzonestorage
            real_quickflow = max(0.0, hbv.k0[i] * (upperzonestorage - hbv.suz[i]))
            upperzonestorage = upperzonestorage - quickflow - real_quickflow

        end

        ### calculations for lower zone ###
        lowerzonestorage = hbv.lowerzonestorage[i] + percolation
        # baseflow in mm/timestep
        baseflow = min(lowerzonestorage, hbv.k4[i] * lowerzonestorage)
        lowerzonestorage = lowerzonestorage - baseflow

        if external_qbase
            directrunoffstorage = quickflow + seepage + real_quickflow
        else
            directrunoffstorage = quickflow + baseflow + real_quickflow
        end

        runoff = max(0.0, directrunoffstorage)

        # update the outputs and states
        hbv.runoff[i] = runoff
        hbv.rainfallplusmelt[i] = rainfallplusmelt
        hbv.soilmoisture[i] = soilmoisture
        hbv.soilevap[i] = soilevap
        hbv.actevap[i] = actevap
        hbv.quickflow[i] = quickflow
        hbv.real_quickflow[i] = real_quickflow
        hbv.upperzonestorage[i] = upperzonestorage
        hbv.lowerzonestorage[i] = lowerzonestorage
        hbv.baseflow[i] = baseflow
        hbv.hbv_seepage[i] = hbv_seepage
        hbv.percolation[i] = percolation
        hbv.capflux[i] = capflux
        hbv.in_upperzone[i] = in_upperzone
    end

end


text:
path,_ = splitdir(@__FILE__)
import Interpolations
using PyPlot
import JLD
import CSV
import Dierckx

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

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

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

include("$path/../setup.jl") # for types
# rc("figure", figsize=(4.0, 2.5))
rc("figure", figsize=(3.5, 3.2))
rc("font", size=10.0)
rc("lines", linewidth=1.5)
rc("lines", markersize=3.0)
rc("legend", frameon=false)
rc("axes.spines", right=false, top=false)
# rc("figure.subplot", left=0.17, bottom=0.21, top=0.97, right=.65)
rc("figure.subplot", left=0.22, bottom=0.32, top=0.99, right=.95)
rc("axes", color_cycle=["348ABD", "A60628", "009E73", "7A68A6", "D55E00", "CC79A7"])
plot_cycle=["#348ABD", "#A60628", "#009E73", "#7A68A6", "#D55E00", "#CC79A7"]
close("all")

numprops = [2,4,8,16,32]

cr_factor = 100#reverse([50,100, 150])
db_factor = reverse([76,65,63,60])

db_factor_plot = reverse(["<76","65","63","60"])

DETAILED = false

detailed_output_accel = []
detailed_output_climb = []
detailed_output_cruise = []
design_accel = []
design_climb = []
design_cruise = []
detailed_system = []

ranges = zeros(length(numprops),length(db_factor))
CLs_accel = zeros(length(numprops),length(db_factor))
CLs_climb = zeros(length(numprops),length(db_factor))
CLs_cruise = zeros(length(numprops),length(db_factor))

Drags_accel = zeros(length(numprops),length(db_factor))
Drags_climb = zeros(length(numprops),length(db_factor))
Drags_cruise = zeros(length(numprops),length(db_factor))

CDs_accel = zeros(length(numprops),length(db_factor))
CDs_climb = zeros(length(numprops),length(db_factor))
CDs_cruise = zeros(length(numprops),length(db_factor))

total_weight = zeros(length(numprops),length(db_factor))
climb_thrust = zeros(length(numprops),length(db_factor))
gama = zeros(length(numprops),length(db_factor))

ETAs_accel = zeros(length(numprops),length(db_factor))
ETAs_climb = zeros(length(numprops),length(db_factor))
ETAs_cruise = zeros(length(numprops),length(db_factor))

BattMasses_accel = zeros(length(numprops),length(db_factor))
BattMasses_climb = zeros(length(numprops),length(db_factor))
BattMasses_cruise = zeros(length(numprops),length(db_factor))

Power_accel = zeros(length(numprops),length(db_factor))
Power_climb = zeros(length(numprops),length(db_factor))
Power_cruise = zeros(length(numprops),length(db_factor))

DBs_accel = zeros(length(numprops),length(db_factor))
DBs_climb = zeros(length(numprops),length(db_factor))
DBs_cruise = zeros(length(numprops),length(db_factor))

Thrust_accel = zeros(length(numprops),length(db_factor))
Thrust_climb = zeros(length(numprops),length(db_factor))
Thrust_cruise = zeros(length(numprops),length(db_factor))

Velocity_accel = zeros(length(numprops),length(db_factor))
Velocity_climb = zeros(length(numprops),length(db_factor))
Velocity_cruise = zeros(length(numprops),length(db_factor))

Rtips = zeros(length(numprops),length(db_factor))

Propulsion_Mass = zeros(length(numprops),length(db_factor))

takeoff_dist = zeros(length(numprops),length(db_factor))
dist_climb = zeros(length(numprops),length(db_factor))
dist_trans = zeros(length(numprops),length(db_factor))
dist_accel = zeros(length(numprops),length(db_factor))

tip_mach_accel = zeros(length(numprops),length(db_factor))
tip_mach_climb = zeros(length(numprops),length(db_factor))
tip_mach_cruise = zeros(length(numprops),length(db_factor))

aoa_accel = zeros(length(numprops),length(db_factor))
aoa_climb = zeros(length(numprops),length(db_factor))
aoa_cruise = zeros(length(numprops),length(db_factor))

mass_prop = zeros(length(numprops),length(db_factor))
blades_prop = zeros(length(numprops),length(db_factor))

span = zeros(length(numprops),length(db_factor))

mean_chord = zeros(length(numprops),length(db_factor))

for i = 1:length(numprops)
    for j = 1:length(db_factor)
    # try
        filename = "$(path)/p$(numprops[i])_cr$(cr_factor)_db$(db_factor[j])correct_span.jld"
        summary = JLD.load(filename)

        detailed_output_accel = summary["detailed_output_accel"]
        detailed_output_climb = summary["detailed_output_climb"]
        detailed_output_cruise = summary["detailed_output_cruise"]
        design_accel = summary["design_accel"]
        design_climb = summary["design_climb"]
        design_cruise = summary["design_cruise"]
        detailed_system = summary["detailed_system"]

        ranges[i,j] = detailed_system.range_true
        CLs_accel[i,j] = detailed_output_accel.Lift/(0.5*1.225*design_accel.velocity^2*design_accel.w_halfspan*2*design_accel.w_chord_pts[1]) #TODO: rho, Sref
        CLs_climb[i,j] = detailed_output_climb.Lift/(0.5*1.225*design_climb.velocity^2*design_climb.w_halfspan*2*design_climb.w_chord_pts[1]) #TODO: rho, Sref
        CLs_cruise[i,j] = detailed_output_cruise.Lift/(0.5*1.225*design_cruise.velocity^2*design_cruise.w_halfspan*2*design_cruise.w_chord_pts[1]) #TODO: rho, Sref

        CDs_accel[i,j] = detailed_output_accel.Drag/(0.5*1.225*design_accel.velocity^2*design_accel.w_halfspan*2*design_accel.w_chord_pts[1]) #TODO: rho, Sref
        CDs_climb[i,j] = detailed_output_climb.Drag/(0.5*1.225*design_climb.velocity^2*design_climb.w_halfspan*2*design_climb.w_chord_pts[1]) #TODO: rho, Sref
        CDs_cruise[i,j] = detailed_output_cruise.Drag/(0.5*1.225*design_cruise.velocity^2*design_cruise.w_halfspan*2*design_cruise.w_chord_pts[1]) #TODO: rho, Sref

        total_weight[i,j] = detailed_system.total_weight

        climb_thrust[i,j] = detailed_output_climb.total_thrust
        gama[i,j] = detailed_system.gama

        ETAs_accel[i,j] = detailed_output_accel.ETA_total
        ETAs_climb[i,j] = detailed_output_climb.ETA_total
        ETAs_cruise[i,j] = detailed_output_cruise.ETA_total

        BattMasses_accel[i,j] = detailed_system.batt_mass_accel
        BattMasses_climb[i,j] = detailed_system.batt_mass_climb
        BattMasses_cruise[i,j] = detailed_system.batt_mass_cruise

        Rtips[i,j] = design_cruise.p_Rtip

        span[i,j] = design_cruise.w_halfspan*2

        Propulsion_Mass[i,j] = max(detailed_output_accel.propwing_mass,detailed_output_climb.propwing_mass,detailed_output_cruise.propwing_mass)

        Power_accel[i,j] = detailed_output_accel.total_thrust*design_accel.velocity
        Power_climb[i,j] = detailed_output_climb.total_thrust*design_climb.velocity
        Power_cruise[i,j] = detailed_output_cruise.total_thrust*design_cruise.velocity

        Thrust_accel[i,j] = detailed_output_accel.total_thrust
        Thrust_climb[i,j] = detailed_output_climb.total_thrust
        Thrust_cruise[i,j] = detailed_output_cruise.total_thrust

        Drags_accel[i,j] = detailed_output_accel.Drag
        Drags_climb[i,j] = detailed_output_climb.Drag
        Drags_cruise[i,j] = detailed_output_cruise.Drag

        Velocity_accel[i,j] = design_accel.velocity
        Velocity_climb[i,j] = design_climb.velocity
        Velocity_cruise[i,j] = design_cruise.velocity

        DBs_accel[i,j] = maximum(detailed_output_accel.db_test)
        DBs_climb[i,j] = maximum(detailed_output_climb.db_test)
        DBs_cruise[i,j] = maximum(detailed_output_cruise.db_test)

        takeoff_dist[i,j] = detailed_system.dist_climb+detailed_system.dist_accel
        dist_accel[i,j] = detailed_system.dist_accel
        dist_trans[i,j] = detailed_system.xtr
        dist_climb[i,j] = detailed_system.x_left

        tip_mach_accel[i,j] = design_accel.p_Rtip*design_accel.p_rpm/60*2*pi/detailed_output_accel.a
        tip_mach_climb[i,j] = design_climb.p_Rtip*design_climb.p_rpm/60*2*pi/detailed_output_climb.a
        tip_mach_cruise[i,j] = design_cruise.p_Rtip*design_cruise.p_rpm/60*2*pi/detailed_output_cruise.a

        aoa_accel[i,j] = design_accel.w_aoa
        aoa_climb[i,j] = design_climb.w_aoa
        aoa_cruise[i,j] = design_cruise.w_aoa

        mean_chord[i,j] = mean(design_cruise.p_chord_pts)

        mass_prop[i,j] = detailed_output_cruise.mass_prop*design_cruise.numprops*design_cruise.blades #same for all three
        blades_prop[i,j] = design_cruise.blades #same for all three

        if DETAILED
            figname = "cl_p_db$(numprops[i])_cr$(cr_factor)_db$(db_factor[j])"
            figure(figname)
            # maxy = zeros(length(db_factor))
            plot(detailed_output_climb.panels.y/maximum(detailed_output_climb.panels.y),detailed_output_climb.cllocal,".-",label = "local cl")
            plot(detailed_output_climb.panels.y/maximum(detailed_output_climb.panels.y),detailed_output_climb.cl,".-",label = "freestream cl")
            xlabel("wing halfspan")
            ylabel("cl")
            # ylim([0,1.05*ceil(Int,maximum(maxy))])
            # legend(loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
            # savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

        end

    # end
    end
end

# High level results
figname = "Takeoff_Speed_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],takeoff_dist[1:end-1,i],"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(takeoff_dist[:,i])
        plot(numprops,takeoff_dist[:,i],"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Takeoff Distance (m)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Range_Props_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],ranges[1:end-1,i]/1000.0,"o-",label="$(db_factor_plot[i])")#"1/$(cr_factor[i])")
    else
        maxy[i] = maximum(ranges[:,i]/1000.0)
        plot(numprops,ranges[:,i]/1000.0,"o-",label="$(db_factor_plot[i])")#"1/$(cr_factor[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Range (km)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Accel_Speed_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],dist_accel[1:end-1,i],"o-",label="$(db_factor_plot[i])")#"1/$(cr_factor[i])")
    else
        maxy[i] = maximum(dist_accel[:,i])
        plot(numprops,dist_accel[:,i],"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Ground Roll Distance (m)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "span_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],span[1:end-1,i],"o-",label="$(db_factor_plot[i])")#"1/$(cr_factor[i])")
    else
        maxy[i] = maximum(span[:,i])
        plot(numprops,span[:,i],"o-",label="$(db_factor_plot[i])")#"1/$(cr_factor[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Span (m)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "CL_climb_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],CLs_climb[1:end-1,i],"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(CLs_climb[:,i])
        plot(numprops,CLs_climb[:,i],"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Wing Lift Coefficient (CL)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "CD_cruise_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],CDs_cruise[1:end-1,i],"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(CDs_cruise[:,i])
        plot(numprops,CDs_cruise[:,i],"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Wing Drag Coefficient (CD)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Noise_climb_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],DBs_climb[1:end-1,i],"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(DBs_climb[:,i])
        plot(numprops,DBs_climb[:,i],"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Sound Pressure Level (dBa)")
# ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "FlightPathAngle_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],gama[1:end-1,i]*180/pi,"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(gama[:,i]*180/pi)
        plot(numprops,gama[:,i]*180/pi,"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Flight Path Angle (degrees)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Total_Mass_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],total_weight[1:end-1,i]/9.81,"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(total_weight[:,i]/9.81)
        plot(numprops,total_weight[:,i]/9.81,"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Total Mass (kg)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Propulsion_Mass_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],Propulsion_Mass[1:end-1,i],"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(Propulsion_Mass[:,i])
        plot(numprops,Propulsion_Mass[:,i],"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Propulsion Mass (kg)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

Propulsion_Mass

figname = "Rtip_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],Rtips[1:end-1,i],"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(Rtips[:,i])
        plot(numprops,Rtips[:,i],"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Propeller Radius (m)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Battery_Cruise_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],BattMasses_cruise[1:end-1,i],"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(BattMasses_cruise[:,i])
        plot(numprops,BattMasses_cruise[:,i],"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Battery Mass (kg)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Battery_Total_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],BattMasses_cruise[1:end-1,i]+BattMasses_accel[1:end-1,i]+BattMasses_climb[1:end-1,i],"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(BattMasses_cruise[:,i]+BattMasses_accel[:,i]+BattMasses_climb[:,i])
        plot(numprops,BattMasses_cruise[:,i]+BattMasses_accel[:,i]+BattMasses_climb[:,i],"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Battery Mass (kg)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Propulsion_ETA_Cruise_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],ETAs_cruise[1:end-1,i]*100.0,"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(ETAs_cruise[:,i]*100.0)
        plot(numprops,ETAs_cruise[:,i]*100.0,"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Propulsion Efficiency (%)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Propulsion_ETA_Accel_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],ETAs_accel[1:end-1,i]*100.0,"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(ETAs_accel[:,i]*100.0)
        plot(numprops,ETAs_accel[:,i]*100.0,"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Propulsion Efficiency (%)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Propulsion_ETA_Climb_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],ETAs_climb[1:end-1,i]*100.0,"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(ETAs_climb[:,i]*100.0)
        plot(numprops,ETAs_climb[:,i]*100.0,"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Propulsion Efficiency (%)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Climb_Power_In_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],Power_climb[1:end-1,i]./ETAs_climb[1:end-1,i],"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(Power_climb[:,i]./ETAs_climb[:,i])
        plot(numprops,Power_climb[:,i]./ETAs_climb[:,i],"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Climb Power In (watts)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Climb_Power_Out_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],Power_climb[1:end-1,i]/1000.0,"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(Power_climb[:,i]/1000.0)
        plot(numprops,Power_climb[:,i]/1000.0,"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Climb Power Out (kW)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Climb_Thrust_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],Thrust_climb[1:end-1,i],"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(Thrust_climb[:,i])
        plot(numprops,Thrust_climb[:,i],"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Total Thrust (N)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Climb_Thrust2Weight_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],Thrust_climb[1:end-1,i]./(total_weight[1:end-1,i]),"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(Thrust_climb[:,i]./(total_weight[:,i]))
        plot(numprops,Thrust_climb[:,i]./(total_weight[:,i]),"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Thrust/Weight")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Climb_NetThrust2Weight_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],(Thrust_climb[1:end-1,i]-Drags_climb[1:end-1,i])./(total_weight[1:end-1,i]),"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum((Thrust_climb[:,i]-Drags_climb[:,i])./(total_weight[:,i]))
        plot(numprops,(Thrust_climb[:,i]-Drags_climb[:,i])./(total_weight[:,i]),"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Net-Thrust/Weight")
ylim([0,0.9])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Takeoff_Energy_Percent_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    # if i = || i==2= 1 #Don't plot 16 props for tight db constraint, very bad
    #     to_percent = 100[1:end-1].0*1(end-1BattMasses_accel[:,i]+BattMasses_climb[:,i])./(BattMasses_accel[:,i]+BattMasses_climb[:,i]+BattMasses_cruise[:,i])
    # else
        to_percent = 100.0*(BattMasses_accel[:,i]+BattMasses_climb[:,i])./(BattMasses_accel[:,i]+BattMasses_climb[:,i]+BattMasses_cruise[:,i])
    # end
    semilogy(numprops,to_percent,"o-",label="$(db_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Takeoff Energy (%)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Velocity_cruise_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],Velocity_cruise[1:end-1,i],"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(Velocity_cruise[:,i])
        plot(numprops,Velocity_cruise[:,i],"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Cruise Velocity (m/s)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Velocity_cruiseMPH_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],Velocity_cruise[1:end-1,i]*2.237,"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(Velocity_cruise[:,i]*2.237)
        plot(numprops,Velocity_cruise[:,i]*2.237,"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Cruise Velocity (mph)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "tip_mach_climb_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],tip_mach_climb[1:end-1,i],"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(tip_mach_climb[:,i])
        plot(numprops,tip_mach_climb[:,i],"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Tip Speed (mach)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "aoa_cruise_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],aoa_cruise[1:end-1,i]*180/pi,"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(aoa_cruise[:,i]*180/pi)
        plot(numprops,aoa_cruise[:,i]*180/pi,"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Wing Angle of Attack (deg)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "aoa_accel_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],aoa_accel[1:end-1,i]*180/pi,"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(aoa_accel[:,i]*180/pi)
        plot(numprops,aoa_accel[:,i]*180/pi,"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Wing Angle of Attack (deg)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "aoa_climb_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],aoa_climb[1:end-1,i]*180/pi,"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(aoa_climb[:,i]*180/pi)
        plot(numprops,aoa_climb[:,i]*180/pi,"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Wing Angle of Attack (deg)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "prop_mean_chord_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],mean_chord[1:end-1,i],"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(mean_chord[:,i])
        plot(numprops,mean_chord[:,i],"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Mean Propeller Chord (m)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "mass_prop_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],mass_prop[1:end-1,i],"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(mass_prop[:,i])
        plot(numprops,mass_prop[:,i],"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Propeller Mass (kg)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "blades_prop_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],blades_prop[1:end-1,i],"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(blades_prop[:,i])
        plot(numprops,blades_prop[:,i],"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Propeller Blades (#)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Velocity_accel_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=1:length(db_factor)
    if false #i == 1 || i==2 #Don't plot 16 props for tight db constraint, very bad
        plot(numprops[1:end-1],Velocity_accel[1:end-1,i],"o-",label="$(db_factor_plot[i])")
    else
        maxy[i] = maximum(Velocity_accel[:,i])
        plot(numprops,Velocity_accel[:,i],"o-",label="$(db_factor_plot[i])")
    end
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Liftoff Speed (m/s)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Relative_Distance_cr100_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=2#1:length(db_factor)

    plot(numprops,dist_accel[:,i],"o-",label="Accel")

    plot(numprops,dist_trans[:,i],"o-",label="Transition")

    plot(numprops,dist_climb[:,i],"o-",label="Climb")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Distance (m)")
# ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Percent_Relative_Distance_db65_db"
figure(figname)
maxy = zeros(length(db_factor))
for i=2#1:length(db_factor)

    plot(numprops,dist_accel[:,i]./takeoff_dist[:,i]*100,"o-",label="Accel")

    plot(numprops,dist_trans[:,i]./takeoff_dist[:,i]*100,"o-",label="Transition")

    plot(numprops,dist_climb[:,i]./takeoff_dist[:,i]*100,"o-",label="Climb")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Relative Distance (%)")
# ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "distance_noise_8"
figure(figname)
maxy = zeros(length(db_factor))
distance_noise_8 = zeros(length(db_factor))
for i=1:length(db_factor)
    distance_noise_8[i] = takeoff_dist[3,i]
end
db_factor[end] = round(Int,DBs_climb[3,4])
plot(db_factor,distance_noise_8,"ko-")
xticks(db_factor)
yticks([0.0;round(distance_noise_8[[1,2,4]],-1)])
xlabel("Sound Pressure Level (dBa)")
ylabel("Takeoff Distance (m)")
ylim([0,1.05*ceil(Int,1.1*maximum(distance_noise_8))])
# legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "distance_noise_16"
figure(figname)
maxy = zeros(length(db_factor))
distance_noise_16 = zeros(length(db_factor))
for i=1:length(db_factor)
    distance_noise_16[i] = takeoff_dist[4,i]
end
db_factor[end] = round(Int,DBs_climb[4,4])
plot(db_factor,distance_noise_16,"ko-")
xticks(db_factor)
yticks([0.0;round(distance_noise_16[[1,2,4]],-1)])
xlabel("Sound Pressure Level (dBa)")
ylabel("Takeoff Distance (m)")
ylim([0,1.05*ceil(Int,1.1*maximum(distance_noise_16))])
# legend(title = "Noise Constraint (dBa)",loc="upper center", ncol = length(db_factor), bbox_to_anchor=(0.4, -.20))#loc="center left", bbox_to_anchor=(1, 0.5))
# savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

ext: jl


text:
function period(a, μ)
    return 2π * sqrt(a^3 / μ)
end

using Test

# Test that the period function works
const a = 450e3 + 6371e3
const μ = 3.986004418e14
@test period(a, μ) ≈ 5606.386763415304

using SPICE: recsph

func = recsph

# Test that `func` is the correct function
xyz = [1.0, 2.0, 3.0]

act = func(xyz)
exp = (3.741657386773941, 0.6405223126794246, 1.1071487177940904)
@test act[1] ≈ exp[1]
@test act[2] ≈ exp[2]
@test act[3] ≈ exp[3]

using LinearAlgebra: norm

function rectangular_to_spherical(xyz)
    x, y, z = xyz
    r = norm(xyz)
    θ = acos(z / r)
    ϕ = atan(y, x)
    return r, θ, ϕ
end

act = rectangular_to_spherical(xyz)
@test act[1] ≈ exp[1]
@test act[2] ≈ exp[2]
@test act[3] ≈ exp[3]


text:
# https://grandjam.shinyapps.io/sirj/
# https://grandjam.shinyapps.io/sirj/


text:
"""
	ALConstraintSet{T}

An [`AbstractConstraintSet`](@ref) that stores the constraint values as well as Lagrange
multiplier and penalty terms for each constraint.

The cost associated with constraint terms in the augmented Lagrangian can be evaluated for

	cost!(J::Vector, ::ALConstraintSet)

which adds the cost at each time step to the vector `J` of length `N`.

The cost expansion for these terms is evaluated along the trajectory `Z` using

	cost_expansion!(E::Objective, conSet::ALConstraintSet, Z)

which also adds the expansion terms to the terms in `E`.

The penalty and multiplier terms can be updated using

	penalty_update!(::ALConstraintSet)
	dual_update!(::ALConstraintSet)

The current set of active constraint (with tolerance `tol`) can be re-calculated using

	update_active_set!(::ALConstraintSet, ::Val{tol})

The maximum penalty can be queried using `max_penalty(::ALConstraintSet)`, and the
penalties and/or multipliers can be reset using

	reset!(::ALConstraintSet)
	reset_penalties!(::ALConstraintSet)
	reset_duals!(::ALConstraintSet)

# Constructor
	ALConstraintSet(::ConstraintList, ::AbstractModel)
	ALConstraintSet(::Problem)
"""
struct ALConstraintSet{T} <: TO.AbstractConstraintSet
    convals::Vector{ALConVal}
    errvals::Vector{ALConVal}          # TODO: is this needed?
	# ∇c_proj::Vector{<:Vector}        # Jacobians of projected constraints
    # λ::Vector{<:Vector}
    # μ::Vector{<:Vector}
	# active::Vector{<:Vector}
    c_max::Vector{T}
    μ_max::Vector{T}
    μ_maxes::Vector{Vector{T}}
	# params::Vector{TO.ConstraintParams{T}}
	p::Vector{Int}
end

function ALConstraintSet(cons::TO.ConstraintList, model::RD.AbstractModel)
    n,m = cons.n, cons.m
    n̄ = RobotDynamics.state_diff_size(model)
    ncon = length(cons)
    useG = model isa RD.LieGroupModel
	errvals = map(1:ncon) do i
        C,c = TO.gen_convals(n̄, m, cons[i], cons.inds[i])
        ALConVal(n̄, m, cons[i], cons.inds[i], C, c, useG)
    end
    convals = map(errvals) do errval
        ALConVal(n, m, errval)
    end
	errvals = convert(Vector{ALConVal}, errvals)
	convals = convert(Vector{ALConVal}, convals)
	# ∇c_proj = map(1:ncon) do i
	# 	[copy(errvals[i].jac[j,1]) for j in eachindex(cons.inds[i])]
	# end
    # λ = map(1:ncon) do i
	# 	p = length(cons[i])
	# 	convals[i].λ
    #     # [@SVector zeros(p) for j in cons.inds[i]]
    # end
    # μ = map(1:ncon) do i
	# 	p = length(cons[i])
	# 	convals[i].μ
    #     # [@SVector ones(p) for j in cons.inds[i]]
    # end
    # a = map(1:ncon) do i
    #     p = length(cons[i])
    #     [@SVector ones(Bool,p) for j in cons.inds[i]]
	# end
	# if ncon == 0
	# 	∇c_proj = Vector{Float64}[]
	# 	λ = Vector{Float64}[]
	# 	μ = Vector{Float64}[]
	# 	a = Vector{Float64}[]
	# end
    c_max = zeros(ncon)
    μ_max = zeros(ncon)
    μ_maxes = [zeros(length(ind)) for ind in cons.inds]
	# params = [TO.ConstraintParams() for con in cons.constraints]
    # ALConstraintSet(convals, errvals, ∇c_proj, λ, μ, a, c_max, μ_max, μ_maxes, params, copy(cons.p))
    ALConstraintSet(convals, errvals, c_max, μ_max, μ_maxes, copy(cons.p))
end

@inline ALConstraintSet(prob::Problem) = ALConstraintSet(prob.constraints, prob.model)

@inline TO.get_convals(conSet::ALConstraintSet) = conSet.convals
@inline TO.get_errvals(conSet::ALConstraintSet) = conSet.errvals

get_duals(conSet::ALConstraintSet) = [cval.λ for cval in conSet.convals]
function set_duals!(conSet::ALConstraintSet, λ)
	for i = 1:length(conSet)
		set_duals!(conSet.convals[i], λ[i])
	end
end
function set_duals!(cval::ALConVal, λ)
	for i = 1:length(cval.λ)
		cval.λ[i] .= λ[i]
	end
end

for method in (:violation!, :∇violation!, :dual_update!, :penalty_update!, :reset_duals!, :reset_penalties!)
	@eval function $method(conSet::ALConstraintSet)
		for i = 1:length(conSet)
			$method(conSet.convals[i])
		end
	end
end

############################################################################################
#                             Augmented Lagrangian Updates
############################################################################################

# function dual_update!(conSet::ALConstraintSet)
#     for i in eachindex(conSet.λ)
#         dual_update!(conSet.convals[i])
# 	end
# end

function dual_update!(conval::ALConVal)
	c = conval.vals
	λ = conval.λ
	μ = conval.μ
	λ_max = conval.params.λ_max
	cone = TO.sense(conval.con)
	# λ_min = TO.sense(conval.con) == Equality() ? -λ_max : zero(λ_max)
	for i in eachindex(conval.inds)
		λ[i] .= dual_update(cone, SVector(λ[i]), SVector(c[i]), SVector(μ[i]), λ_max)
	end
end

function dual_update(::Equality, λ, c, μ, λmax)
	λbar = λ + μ .* c
	return clamp.(λbar, -λmax, λmax)
end

function dual_update(::Inequality, λ, c, μ, λmax)
 	λbar = λ + μ .* c
	return clamp.(λbar, 0, λmax)  # project onto the dual cone via max(0,x)
end

function dual_update(cone::SecondOrderCone, λ, c, μ, λmax)
	 λbar = λ - μ .* c
	 return TO.projection(cone, λbar)  # project onto the dual cone
end

# function penalty_update!(conSet::ALConstraintSet)
# 	for i in eachindex(conSet.μ)
# 		penalty_update!(conSet.convals[i])
# 	end
# end

function penalty_update!(cval::ALConVal)
	μ = cval.μ
	ϕ = cval.params.ϕ
	μ_max = cval.params.μ_max
	for i = 1:length(μ)
		μ[i] .*= ϕ
		clamp!(μ[i], 0, μ_max)
	end
end

# Active Set
function update_active_set!(conSet::ALConstraintSet, val::Val{tol}=Val(0.0)) where tol
	for i in 1:length(conSet)
		# update_active_set!(conSet.active[i], conSet.λ[i], conSet.convals[i], val)
		update_active_set!(conSet.convals[i], val)
	end
end

function update_active_set!(conval::ALConVal, ::Val{tol}) where tol
	a = conval.active
	λ = conval.λ
	if TO.sense(conval.con) == TO.Inequality()
		for i in eachindex(a)
			a[i] = @. (conval.vals[i] >= -tol) | (λ[i] > zero(tol))
		end
	end
end

"""
	max_penalty(conSet::ALConstraintSet)

Calculate the maximum constrained penalty across all constraints.
"""
function max_penalty(conSet::ALConstraintSet)
	max_penalty!(conSet)
	maximum(conSet.μ_max)
end

function max_penalty!(conSet::ALConstraintSet{T}) where T
    conSet.c_max .*= 0
	for i in 1:length(conSet)
		maxes = conSet.μ_maxes[i]::Vector{T}
        max_penalty!(maxes, conSet.convals[i])
        conSet.μ_max[i] = maximum(maxes)
    end
end

function max_penalty!(μ_max::Vector{<:Real}, cval::ALConVal)
    for i in eachindex(cval.μ)
        μ_max[i] = maximum(cval.μ[i])
    end
    return nothing
end

############################################################################################
#                             Constraint Penalties
############################################################################################
# function projection_jacobians!(conSet::ALConstraintSet)
# 	for i in eachindex(conSet.convals)
# 		projection_jacobians!(conSet.convals[i])
# 	end
# end

# """
# 	constraint_penalty!(conSet::ALConstraintSet)

# Evaluate the penalty of the constraints. For equality constraints, this is the constraint
# itself. For generalized inequalities, this is the projection of the constraint onto the
# appropriate cone.
# """
# function constraint_penalty!(conSet::ALConstraintSet)
# 	for i in eachindex(conSet.convals)
# 		constraint_penalty!(conSet.convals[i], conSet.λ[i])
# 	end
# end

# function constraint_penalty!(conval::TO.ConVal, λ)
# 	for i in eachindex(conval.inds)
# 		conval.vals2[i] = penalty(TO.sense(conval.con), conval.vals[i], λ[i])
# 	end
# end

# """
# 	penalty_jacobian!(conSet::ALConstraintSet)

# Evaluate the Jacobian of the constraint penalty. For equality constraints, this is simply
# the Jacobian of the constraint itself. For generalized inequalities, this is `∇proj*∇c`
# where `∇proj` is the Jacobian of the projection of the constraint onto the corresponding
# cone.
# """
# function penalty_jacobian!(conSet::ALConstraintSet)
# 	for i in eachindex(conSet.errvals)
# 		penalty_jacobian!(conSet.∇c_proj[i], conSet.errvals[i], conSet.λ[i])
# 	end
# end

# function penalty_jacobian!(∇c_proj::Vector, conval::TO.ConVal, λ::Vector)
# 	if size(conval.jac, 2) > 1
# 		throw(ErrorException("Constraint projection not supported for CoupledConstraints"))
# 	end
# 	for i in eachindex(conval.inds)
# 		∇penalty!(TO.sense(conval.con), ∇c_proj[i], conval.jac[i], conval.vals[i], λ[i])
# 	end
# end

############################################################################################
#                                        Cost
############################################################################################

# function TO.cost!(J::Vector{<:Real}, conSet::ALConstraintSet)
# 	for i in eachindex(conSet.convals)
# 		TO.cost!(J, conSet.convals[i], conSet.λ[i], conSet.μ[i], conSet.active[i])
# 	end
# end

# function TO.cost!(J::Vector{<:Real}, conval::ALConVal, λ::Vector{<:StaticVector},
# 		μ::Vector{<:StaticVector}, a::Vector{<:StaticVector})
# 	for (i,k) in enumerate(conval.inds)
# 		c = SVector(conval.vals[i])
# 		Iμ = Diagonal(SVector(μ[i] .* a[i]))
# 		J[k] += λ[i]'c .+ 0.5*c'Iμ*c
# 	end
# end

# function TO.cost_expansion!(E::Objective, conSet::ALConstraintSet, Z::AbstractTrajectory,
# 		init::Bool=false, rezero::Bool=false)
# 	for i in eachindex(conSet.errvals)
# 		TO.cost_expansion!(E, conSet.convals[i], conSet.λ[i], conSet.μ[i], conSet.active[i])
# 	end
# end

# @generated function TO.cost_expansion!(E::QuadraticObjective{n,m}, conval::ALConVal{C}, λ, μ, a) where {n,m,C}
# 	if C <: TO.StateConstraint
# 		expansion = quote
# 			cx = ∇c
# 			E[k].Q .+= cx'Iμ*cx
# 			E[k].q .+= cx'g
# 		end
# 	elseif C <: TO.ControlConstraint
# 		expansion = quote
# 			cu = ∇c
# 			E[k].R .+= cu'Iμ*cu
# 			E[k].r .+= cu'g
# 		end
# 	elseif C<: TO.StageConstraint
# 		ix = SVector{n}(1:n)
# 		iu = SVector{m}(n .+ (1:m))
# 		expansion = quote
# 			cx = ∇c[:,$ix]
# 			cu = ∇c[:,$iu]
# 			E[k].Q .+= cx'Iμ*cx
# 			E[k].q .+= cx'g
# 			E[k].H .+= cu'Iμ*cx
# 			E[k].R .+= cu'Iμ*cu
# 			E[k].r .+= cu'g
# 		end
# 	else
# 		throw(ArgumentError("cost expansion not supported for CoupledConstraints"))
# 	end
# 	quote
# 		for (i,k) in enumerate(conval.inds)
# 			∇c = SMatrix(conval.jac[i])
# 			c = conval.vals[i]
# 			Iμ = Diagonal(a[i] .* μ[i])
# 			g = Iμ*c .+ λ[i]

# 			$expansion
# 		end
# 	end
# end

############################################################################################
#                                  RESET
############################################################################################

function reset!(conSet::ALConstraintSet)
    reset_duals!(conSet)
    reset_penalties!(conSet)
end

# function reset_duals!(conSet::ALConstraintSet)
# 	for i = 1:length(conSet)
# 		reset_duals!(conSet.convals[i])
# 	end
# end

# function reset_penalties!(conSet::ALConstraintSet)
# 	for i = 1:length(conSet)
# 		reset_penalties!(conSet.convals[i])
# 	end
# end

"""
	link_constraints!(set1, set2)

Link any common constraints between `set1` and `set2` by setting elements in `set1` to point
to elements in `set2`
"""
function link_constraints!(set1::ALConstraintSet, set2::ALConstraintSet)
	# Find common constraints
	links = Tuple{Int,Int}[]
	for (i,con1) in enumerate(set1)
		for (j,con2) in enumerate(set2)
			if con1 === con2
				push!(links, (i,j))
			end
		end
	end

	# Link values
	for (i,j) in links
		set1.convals[i] = set2.convals[j]
		set1.errvals[i] = set2.errvals[j]
		# set1.active[i] = set2.active[j]
		# set1.λ[i] = set2.λ[j]
		# set1.μ[i] = set2.μ[j]
	end
	return links
end

function shift_fill!(conSet::ALConstraintSet, n=1)
	for i = 1:length(conSet)
		shift_fill!(conSet.convals[i], n)
	end
end

############################################################################################
#                                    Solver Options
############################################################################################
function reset!(conSet::ALConstraintSet{T}, opts::SolverOptions{T}) where T
    # if !isnan(opts.dual_max)
    #     for params in conSet.params
    #         params.λ_max = opts.dual_max
    #     end
    # end
    # if !isnan(opts.penalty_max)
    #     for params in conSet.params
    #         params.μ_max = opts.penalty_max
    #     end
    # end
    # if !isnan(opts.penalty_initial)
    #     for params in conSet.params
    #         params.μ0 = opts.penalty_initial
    #     end
    # end
    # if !isnan(opts.penalty_scaling)
    #     for params in conSet.params
    #         params.ϕ = opts.penalty_scaling
    #     end
	# end
	for i = 1:length(conSet)
		set_params!(conSet.convals[i], opts)
	end
    if opts.reset_duals
        reset_duals!(conSet)
    end
    if opts.reset_penalties
        reset_penalties!(conSet)
    end
end


text:
using LinearAlgebra
using Test

include("../egrss.jl")

include("test_generators.jl")
include("test_potrf.jl")
include("test_ldl.jl")
include("test_full.jl")
include("test_full_tril.jl")
include("test_gemv.jl")
include("test_symv.jl")
include("test_trmv.jl")
include("test_trsv.jl")
include("test_trtri.jl")
include("test_trnrms.jl")
include("test_bigfloat.jl")
include("test_rational.jl")

max_stars_repo_name: martinandersen/EGRSS

max_issues_repo_name: martinandersen/EGRSS

max_forks_repo_name: martinandersen/EGRSS


text:
struct Group{NT}
    kwargs::NT
    function Group(; kwargs...)
        nt = values(kwargs)
        NT = typeof(nt)
        return new{NT}(nt)
    end
end

Group(s) = Group(color = s)
apply(f::Function, g::Group) = Group(; map(t -> apply(f, t), g.kwargs)...)

to_string(t::Tuple) = join(t, ", ")
to_string(t::Any) = string(t)
to_string(nt::NamedTuple) = join(("$a = $b" for (a, b) in pairs(nt)), ", ")

struct Observations; end
const observations = Observations()
Base.string(::Observations) = "observations"

function sortpermby(t::IndexedTable, ::Observations; return_keys = false)
    perm = Base.OneTo(length(t))
    return return_keys ? (perm, perm) : perm
end

function apply_postprocess(t::IndexedTable, res; select, postprocess)
    cols = Tuple(fieldarrays(res))
    colinds = to_tuple(lowerselection(t, select))
    N = length(colinds)
    cols_trimmed = cols[1:N]
    res = map(cols_trimmed, colinds) do col, ind
        isa(ind, Integer) && ind > 0 || return col
        name = colnames(t)[ind]
        haskey(postprocess, name) || return col
        f = postprocess[name]
        return collect_columns(apply(f, el) for el in col)
    end
    StructArray((res..., cols[N+1:end]...))
end

function series2D(s::StructVector; ribbon = false)
    kwargs = Dict{Symbol, Any}()
    xcols, ycols = map(columntuple, fieldarrays(s))
    x, y = xcols[1], ycols[1]
    yerr = ifelse(ribbon, :ribbon, :yerr)
    length(xcols) == 2 && (kwargs[:xerr] = xcols[2])
    length(ycols) == 2 && (kwargs[yerr] = ycols[2])
    return (x, y), kwargs
end

series2D(t, g = Group(); kwargs...) = series2D(nothing, t, g; kwargs...)

function series2D(f::Union{Nothing, FunctionOrAnalysis}, t, g = Group();
                  select, error = automatic, kwargs...)
    isa(g, Group) || (g = Group(g))
    by = _flatten(g.kwargs)
    err_cols = error === automatic ? () : to_tuple(error)
    sel_cols = (to_tuple(by)..., to_tuple(select)..., err_cols...)
    coltable = Tables.columntable(TableOperations.select(t, sel_cols...))
    coldict = Dict(zip(sel_cols, keys(coltable)))
    to_symbol = i -> coldict[i]
    t = table(coltable, copy=false)
    return series2D(f, t, apply(to_symbol, g);
                    select = apply(to_symbol, select),
                    error = apply(to_symbol, error),
                    kwargs...)
end

function series2D(f::Union{Nothing, FunctionOrAnalysis}, t′::IndexedTable, g = Group(); select, postprocess = NamedTuple(),
    error = automatic, ribbon=false, stats=summary, filter=isfinitevalue, transform=identity, min_nobs=2, kwargs...)

    t = dropmissing(t′, select)
    isa(g, Group) || (g = Group(g))
    group = g.kwargs
    if isempty(group)
        itr = ("" => :,)
    else
        by = _flatten(group)
        perm, group_rows = sortpermby(t, by, return_keys = true)
        itr = finduniquesorted(group_rows, perm)
    end
    keys = similar(group_rows, 0)
    vals = collect_columns_flattened(
        Base.Generator(itr) do (key, idxs)
            v = compute_summary(
                f,
                view(t, idxs),
                error;
                min_nobs=min_nobs,
                select=select,
                stats=stats,
                filter=filter,
                transform=transform,
            )
            append!(keys, fill(key, length(v)))
            return v
        end
    )
    res = apply_postprocess(t, vals; select = select, postprocess = postprocess)
    plot_args, plot_kwargs = series2D(res; ribbon = ribbon)
    plot_kwargs[:group] = columns(keys)
    grpd = collect_columns(key for (key, _) in itr)
    style_kwargs = Dict(kwargs)
    for (key, val) in pairs(group)
        col = rows(grpd, val)
        s = unique(sort(col))
        d = Dict(zip(s, 1:length(s)))
        style = get(style_kwargs, key) do
            style_dict[key]
        end
        plot_kwargs[key] = permutedims(access_style(style, getindex.(Ref(d), col)))
    end
    get!(plot_kwargs, :color, "black")
    plot_args, plot_kwargs
end

function access_style(st, n::AbstractArray)
    [access_style(st, i) for i in n]
end

function access_style(st, n::Integer)
    v = vec(st)
    m = ((n-1) % length(v))+1
    v[m]
end

_flatten(t) = IterTools.imap(to_tuple, t) |>
    Iterators.flatten |>
    IterTools.distinct |>
    Tuple


text:
using FEPreProcessing
using Test

@testset "FEPreProcessing.jl" begin
    # Write your tests here.
end


text:
# This file was generated, do not modify it. # hide
using Gnuplot, Random
Random.seed!(123)
goffset = "set offsets graph .05, graph .05, graph .05, graph .05"
@gp(rand(30), rand(30), "w p pt 6 ps 2 lc '#D55E00' lw 2 t 'marker'",
    goffset, title = "scatter plot", xlabel ="X", ylabel = "Y",
    "set key b l box opaque")
    save(term="pngcairo size 600,400", output="plt_ex11.png") # hide


text:
expr_removebackslash = GeneralTools.expr_removebackslash;

@testset "Check if regular expression pattern `$expr_removebackslash` works fine:" begin
	@test replace("\\\\using_gadfly_10_1.png", expr_removebackslash => "") == "using_gadfly_10_1.png"
	@test replace("\\blas\\(0)cu2\\using_gadfly_10_1.png", expr_removebackslash => "") == "blas\\(0)cu2\\using_gadfly_10_1.png"
	@test replace("/using_gadfly\\asdf_fd/_10_1.png", expr_removebackslash => "") == "using_gadfly\\asdf_fd/_10_1.png"
	@test replace("/using_gadfly/asdf_fd/_10_1.png", expr_removebackslash => "") == "using_gadfly/asdf_fd/_10_1.png"
	@test replace("using_gadfly/asdf_fd/_10_1.png", expr_removebackslash => "") == "using_gadfly/asdf_fd/_10_1.png"
	@test replace("using_gadfly_10_1.png", expr_removebackslash => "") == "using_gadfly_10_1.png"
end

@testset "Check getdoc:" begin
	s = [
	"% DOC BEGIN",
	"% first line",
	"% bla haha",
	"% foobar is good",
	"% DOC END",
	"function blablabla(hello, world)",
	"x = 1; ",
	"y = 1; % with comment. ",
	"% y = 1;",
	"%% only comment"
	];
	doc = GeneralTools._getrawcomment(s);

	@test doc[1] == "% DOC BEGIN"
	@test doc[end] == "% DOC END"

end

max_stars_repo_name: CGRG-lab/HypertextTools

max_issues_repo_name: CGRG-lab/HypertextTools

max_forks_repo_name: CGRG-lab/HypertextTools


text:
#
# Iterator to join left and right streams on position
#

#
# uses vector as queue for right iterator
#
# iterator state is
# (lelt_s, relt_s, rq, rqidx)
# where lelt_s is current elt from left and state
#       relt_s is element, state that just did not make it into rq
#       rq is queue of rows from right to be reported withh lelt_s[1]
#       rqidx is index into rq to be reported on next iteration
#

#
# How do we combine elts from left and right stream?
#
# If we use an inner join, can simply use mergex from utils.jl.
# If we use a left join, all columns from right stream can be missing,
# i.e. we need to declare their types as Union{Missing, coltype}
# implemented as mergeleftx(left, right) or mergeleftx(left, R) in utils.jl
#

export gorjoin

using Tables

abstract type AbstractJoin{L,R} end

struct InnerJoin{L,R} <: AbstractJoin{L,R}
    left::L
    right::R
    window::Int64
end

struct LeftJoin{L,R} <: AbstractJoin{L,R}
    left::L
    right::R
    window::Int64
end

Base.IteratorEltype(::Type{<:AbstractJoin{L,R}}) where {L,R} = Base.HasEltype()

Base.eltype(::Type{InnerJoin{L,R}}) where {L,R} =
    Base.return_types(mergex, (eltype(L), eltype(R)))[1]

Base.eltype(::Type{LeftJoin{L,R}}) where {L,R} =
    Base.return_types(mergeleftx, (eltype(L), Type{eltype(R)}))[1]

Base.IteratorSize(::Type{<:AbstractJoin{L,R}}) where {L,R} = Base.SizeUnknown()

Tables.istable(j::AbstractJoin{L,R}) where {L,R} = true
Tables.rowaccess(j::AbstractJoin{L,R}) where {L,R} = true
Tables.schema(j::AbstractJoin{L,R}) where {L,R} = Tables.Schema(eltype(j))
Tables.rows(j::AbstractJoin{L,R}) where {L,R} = j

"""
joins streams `left` and `right` on (:Chrom, :Pos) allowing
`window` difference in position

leftjoin: should left join be performed
"""
gorjoin(left, right; leftjoin = false, window = 0) =
    leftjoin ?
    LeftJoin{typeof(left), typeof(right)}(left, right, window) :
    InnerJoin{typeof(left), typeof(right)}(left, right, window)

gorjoin(right; leftjoin = false, window = 0) =
    left -> gorjoin(left, right; leftjoin = leftjoin, window = window)

# given the current lelt_s,
# update rq to contain all rows to be reported

function advance_right!(j::AbstractJoin, state)
    lelt_s, relt_s, rq, rqidx = state

    # compute new criteria based on position for rows to be included in rq
    min = (lelt_s[1].Chrom, lelt_s[1].Pos - j.window)
    max = (lelt_s[1].Chrom, lelt_s[1].Pos + j.window)

    # remove elts from front of rq that do not fulfil criteria
    while length(rq) > 0 && (rq[1].Chrom, rq[1].Pos) < min
	popfirst!(rq)
    end

    # add elts from right to rq if they fulfil criteria
    while relt_s ≠ nothing && (relt_s[1].Chrom, relt_s[1].Pos) <= max
	push!(rq, relt_s[1])
	relt_s = iterate(j.right, relt_s[2])
    end

    # return new state
    return (lelt_s, relt_s, rq, 1)
end

function Base.iterate(j::AbstractJoin)

    # prime state for first iteration
    lelt_s = iterate(j.left)
    relt_s = iterate(j.right)

    lelt_s == nothing && return nothing

    rq = Vector{eltype(j.right)}()

    state = advance_right!(j, (lelt_s, relt_s, rq, 1))

    # and iterate
    iterate(j, state)
end

function Base.iterate(j::InnerJoin{L,R}, state) where {L,R}

    lelt_s, relt_s, rq, rqidx = state

    while true

	if rqidx <= length(rq)
	    # could put further join conditions here
	    # if fulfilled, return
	    # else rqidx += 1
	    return mergex(lelt_s[1], rq[rqidx]), (lelt_s, relt_s, rq, rqidx + 1)
	end

	if rqidx > length(rq)
	    lelt_s = iterate(j.left, lelt_s[2])
	    lelt_s == nothing && return nothing
	    lelt_s, relt_s, rq, rqidx = advance_right!(j, (lelt_s, relt_s, rq, rqidx))
	end

    end
end

function Base.iterate(j::LeftJoin{L,R}, state) where {L,R}

    lelt_s, relt_s, rq, rqidx = state

    while true

	if length(rq) == 0 && rqidx == 1
	    return mergeleftx(lelt_s[1], eltype(R)), (lelt_s, relt_s, rq, rqidx + 1)
	end

	if rqidx <= length(rq)
	    # could put further join conditions here
	    # if fulfilled, return
	    # else rqidx += 1
	    return mergeleftx(lelt_s[1], rq[rqidx]), (lelt_s, relt_s, rq, rqidx + 1)
	end

	if rqidx > length(rq)
	    lelt_s = iterate(j.left, lelt_s[2])
	    lelt_s == nothing && return nothing
	    lelt_s, relt_s, rq, rqidx = advance_right!(j, (lelt_s, relt_s, rq, rqidx))
	end

    end
end


text:
function test_NfOrd_Idl()
  print("Testing NfOrd/Idl ... ")

  Qx, x = PolynomialRing(FlintQQ, "x")

  K1, a1 = NumberField(x^3 - 2, "a")
  O1 = Order(K1, Hecke.FakeFmpqMat(ZZ[1 0 0; 0 2 0; 0 0 4], one(ZZ)))

  # construction

  I = @inferred ideal(O1, -17)
  @test order(I) == O1
  @test I.princ_gen_special[1] == 1
  @test I.princ_gen_special[2] == 17
  @test I.princ_gen == O1(-17)
  @test I.basis_mat == MatrixSpace(ZZ, 3, 3)(17)

  J = @inferred ideal(O1, ZZ(-17))
  @test order(J) == O1
  @test J.princ_gen_special[1] == 2
  @test J.princ_gen_special[3] == ZZ(17)
  @test J.princ_gen == O1(-17)
  @test J.basis_mat == MatrixSpace(ZZ, 3, 3)(17)

  K = @inferred ideal(O1, O1(-17))
  @test K.princ_gen == O1(-17)
  @test K.basis_mat == MatrixSpace(ZZ, 3, 3)(17)

  M = @inferred O1(-17)*O1
  L = @inferred O1*O1(-17)

  @test I == J && J == K && K == M && M == L && L == I

  # minimum

  @test 17 == @inferred minimum(I)
  @test 17 == @inferred minimum(J)
  @test 17 == @inferred minimum(K)
  @test 17 == @inferred minimum(I)
  @test 17 == @inferred minimum(J)
  @test 17 == @inferred minimum(K)

  # test

  @test 4913 == @inferred norm(I)
  @test 4913 == @inferred norm(J)
  @test 4913 == @inferred norm(K)
  @test 4913 == @inferred norm(I)
  @test 4913 == @inferred norm(J)
  @test 4913 == @inferred norm(K)

  # deepcopy

  L = @inferred deepcopy(I)
  @test L == I
  @test order(L) == O1

  # basis

  b = @inferred basis(I)
  @test b == NfOrdElem{NfOrdGen}[ O1(17), O1(34*a1), O1(68*a1^2) ]

  # basis matrix

  M = @inferred ideal(O1, O1(4*a1^2))

  b = @inferred basis_mat(M)
  @test b == ZZ[16 0 0; 0 16 0; 0 0 1]

  b = @inferred basis_mat_inv(M)
  @test b == Hecke.FakeFmpqMat(ZZ[1 0 0; 0 1 0; 0 0 16], ZZ(16))

  b = @inferred basis_mat(M)
  @test b == ZZ[16 0 0; 0 16 0; 0 0 1]

  b = @inferred basis_mat_inv(M)
  @test b == Hecke.FakeFmpqMat(ZZ[1 0 0; 0 1 0; 0 0 16], ZZ(16))

  # inclusion

  M = @inferred ideal(O1, O1(4*a1^2))
  @test @inferred in(O1(4*a1^2), M)
  @test @inferred !in(O1(2*a1), M)

  # binary

  K = @inferred I + J # I == J
  @test K == I
  K = @inferred ideal(O1, 15) + ideal(O1, 10)
  @test K == ideal(O1, 5)

  K = @inferred intersection(ideal(O1, 15), ideal(O1, 10))
  @test K == ideal(O1, 30)
  K = @inferred lcm(ideal(O1, 15), ideal(O1, 10))
  @test K == ideal(O1, 30)

  K = @inferred ideal(O1, 15) * ideal(O1, 10)
  @test K == ideal(O1, 150)
  K = @inferred ideal(O1, O1(2*a1)) * ideal(O1, O1(4*a1^2))
  @test K == ideal(O1, 16)

  # adhoc binary

  I = ideal(O1, O1(2*a1))
  J = ideal(O1, 3)
  K = ideal(O1, fmpz(3))

  @test ideal(O1, O1(10*a1)) == @inferred 5 * I
  @test ideal(O1, O1(10*a1)) == @inferred I * 5

  @test ideal(O1, O1(10*a1)) == @inferred fmpz(5) * I
  @test ideal(O1, O1(10*a1)) == @inferred I * fmpz(5)

  @test ideal(O1, 30) == @inferred 10 * J
  @test ideal(O1, 30) == @inferred J * 10

  @test ideal(O1, 30) == @inferred fmpz(10) * J
  @test ideal(O1, 30) == @inferred J * fmpz(10)

  @test ideal(O1, 30) == @inferred 10 * K
  @test ideal(O1, 30) == @inferred K * 10

  @test ideal(O1, 30) == @inferred fmpz(10) * K
  @test ideal(O1, 30) == @inferred K * fmpz(10)

  @test ideal(O1, O1(4*a1^2)) == @inferred I * O1(2*a1)
  @test ideal(O1, O1(4*a1^2)) == @inferred O1(2*a1) * I

  # idempotents

  I = ideal(O1, 2)
  J = ideal(O1, 3)
  e, f = @inferred idempotents(I, J)
  @test in(e, I)
  @test in(f, J)
  @test e + f == 1

  I = ideal(O1, O1(2*a1))
  J = ideal(O1, O1(1 - 2*a1))
  e, f = @inferred idempotents(I, J)
  @test in(e, I)
  @test in(f, J)
  @test e + f == 1

  @test_throws ErrorException idempotents(I, I)

  # modular reduction

  I = ideal(O1, 10)
  a = O1([11, 11, 11])
  b = @inferred mod(a, I)
  @test b == O1([1, 1, 1])

  I = ideal(O1, O1(2*a1 + 4*a1^2))
  a = O1([544, 1, 1])
  b = @inferred mod(a, I)
  @test iszero(b)

  c = [ Hecke.fmpz_preinvn_struct(fmpz(544)),
        Hecke.fmpz_preinvn_struct(fmpz(1)),
        Hecke.fmpz_preinvn_struct(fmpz(1)) ]

  b = @inferred mod(a, I, c)
  @test iszero(b)

  # pradical

  I = @inferred pradical(O1, 2)
  @test I == ideal(O1, ZZ[2 0 0; 0 1 0; 0 0 1])

  println("PASS")
end


text:
include("io.jl")
include("functions.jl")
include("radiation.jl")
include("voronoi_utils.jl")
include("characteristics.jl")
include("irregular_ray_tracing.jl")

function J_λ_regular(S_λ::AbstractArray,
                     α_cont::AbstractArray,
                     atmos::Atmosphere,
                     quadrature::String)

    # Ω = (θ, φ), space angle
    weights, θ_array, ϕ_array, n_points = read_quadrature(quadrature)

    J = zero(S_λ)

    for i in 1:n_points
        θ = θ_array[i]
        ϕ = ϕ_array[i]
        k = [cos(θ*π/180), cos(ϕ*π/180)*sin(θ*π/180), sin(ϕ*π/180)*sin(θ*π/180)]
        if θ > 90
            I_0 = blackbody_λ.(500u"nm", atmos.temperature[1,:,:])
            J += weights[i]*short_characteristics_up(k, S_λ, α_cont, atmos, I_0=I_0)
        elseif θ < 90
            I_0 = zero(S_λ[1, :, :])
            J += weights[i]*short_characteristics_down(k, S_λ, α_cont, atmos, I_0=I_0)
        end
    end
    return J
end

function J_λ_voronoi(S_λ::AbstractArray,
                     α_cont::AbstractArray,
                     sites::VoronoiSites,
                     quadrature::String)

    # Ω = (θ, φ), space angle
    weights, θ_array, ϕ_array, n_points = read_quadrature(quadrature)

    J = zero(S_λ)

    n_sweeps = 3
    Nran = 1

    for i in 1:n_points
        θ = θ_array[i]
        ϕ = ϕ_array[i]
        k = [cos(θ*π/180), cos(ϕ*π/180)*sin(θ*π/180), sin(ϕ*π/180)*sin(θ*π/180)]
        if θ > 90
            bottom_layer = sites.layers_up[2] - 1
            bottom_layer_idx = sites.perm_up[1:bottom_layer]
            I_0 = blackbody_λ.(500u"nm", sites.temperature[bottom_layer_idx])
            J += weights[i]*Delaunay_upII(k, S_λ, I_0, α_cont, sites, n_sweeps)
        elseif θ < 90
            top_layer = sites.layers_down[2] - 1
            I_0 = zeros(top_layer)u"kW*nm^-1*m^-2"
            J += weights[i]*Delaunay_downII(k, S_λ, I_0, α_cont, sites, n_sweeps)
        end
    end
    return J
end

function Λ_regular(ϵ::AbstractFloat,
                   maxiter::Integer,
                   atmos::Atmosphere,
                   quadrature::String)

    # choose a wavelength
    λ = 500u"nm"  # nm

    # Lte populations
    LTE_pops = LTE_populations(atmos)

    # Find continuum extinction (only with Thomson and Rayleigh)
    α_s = α_scattering(λ, atmos.temperature, atmos.electron_density*1.0,
                       LTE_pops[:,:,:,1])

    α_a = α_absorption.(λ, atmos.temperature, atmos.electron_density*1.0,
                        LTE_pops[:,:,:,1].+LTE_pops[:,:,:,2],
                        LTE_pops[:,:,:,3])

    α_cont = α_s + α_a
    # destruction
    ε_λ = α_a./α_cont

    thick = ε_λ .> 1e-4

    # Start with the source function as the Planck function
    B_0 = blackbody_λ.(λ, atmos.temperature)
    S_new = B_0

    S_old = zero(S_new)

    i=0

    local J_new
    # check where ε < 1e-2, cut above heights
    while criterion(S_new, S_old, ϵ, i, maxiter, thick)
        S_old = copy(S_new)
        J_new = J_λ_regular(S_old, α_cont, atmos, quadrature)
        S_new = (1 .- ε_λ).*J_new .+ ε_λ.*B_0
        i+=1
    end

    if i == maxiter
        println("Did not converge inside scope")
        return J_new, S_new, α_cont
    end

    println("Converged in $i iterations")

    return J_new, S_new, α_cont
end

function Λ_voronoi(ϵ::AbstractFloat,
                   maxiter::Integer,
                   sites::VoronoiSites,
                   quadrature::String)
    println("---Iterating---")

    # choose a wavelength
    λ = 500u"nm"  # nm

    # Lte populations
    LTE_pops = LTE_populations(sites)

    # Find continuum extinction (only with Thomson and Rayleigh)
    α_s = α_scattering(λ, sites.temperature, sites.electron_density*1.0,
                       LTE_pops[:,1])

    α_a = α_absorption.(λ, sites.temperature, sites.electron_density*1.0,
                        LTE_pops[:,1].+LTE_pops[:,2],
                        LTE_pops[:,3])

    α_cont = α_s + α_a

    # destruction
    ε_λ = α_a ./ α_cont

    thick = ε_λ .> 1e-4

    # Start with the source function as the Planck function
    B_0 = blackbody_λ.(λ, sites.temperature)
    S_new = B_0

    S_old = zero(S_new)

    i=0

    local J_new
    # check where ε < 1e-2, cut above heights
    while criterion(S_new, S_old, ϵ, i, maxiter, thick)
        S_old = copy(S_new)
        J_new = J_λ_voronoi(S_old, α_cont, sites, quadrature)
        S_new = (1 .- ε_λ).*J_new .+ ε_λ.*B_0
        i+=1
    end

    if i == maxiter
        println("Did not converge inside scope")
        return J_new, S_new, α_cont
    end

    println("Converged in $i iterations")

    return J_new, S_new, α_cont
end

function criterion(S_new::Array{<:UnitsIntensity_λ, 3},
                   S_old::Array{<:UnitsIntensity_λ, 3},
                   ϵ::Float64,
                   i::Int,
                   maxiter::Int,
                   indcs)

    diff = maximum(abs.(1 .- S_old[indcs]./S_new[indcs])) |> Unitful.NoUnits
    if isnan(diff)
        println("NaN DIFF!, index $l")
    end
    if i > 0
        println("   Rel. diff.: $diff")
    end
    println("Iteration $(i+1)...")

    diff > ϵ && i < maxiter
end

function criterion(S_new::Vector{<:UnitsIntensity_λ},
                   S_old::Vector{<:UnitsIntensity_λ},
                   ϵ::Float64,
                   i::Int,
                   maxiter::Int,
                   indcs)

    diff = maximum(abs.(1 .- S_old[indcs]./S_new[indcs])) |> Unitful.NoUnits
    if isnan(diff)
        println("NaN DIFF!, index $l")
    end
    if i > 0
        println("   Rel. diff.: $diff")
    end
    println("Iteration $(i+1)...")

    diff > ϵ && i < maxiter
end

max_stars_repo_name: meudnaes/VoronoiRT

max_issues_repo_name: meudnaes/VoronoiRT

max_forks_repo_name: meudnaes/VoronoiRT


text:
# include module
include(joinpath("src", "axelrodnet.jl"))

run_experiment(
    experiment_name="low_agentcount",
    agentcount=100,
    n_iter=1_000_000,
    nettopology=1,
    networkprops=Dict("grid_height" => 10),
    socialbotfrac=0.0,
    rndseed=0,
    repcount=10,
    export_every_n=10_000,
    export_experiment=true,
    aggregated=1,
    keep_rawnet=false
)

run_experiment(
    experiment_name="high_agentcount",
    agentcount=10_000,
    n_iter=2_000_000,
    nettopology=1,
    networkprops=Dict("grid_height" => 100),
    socialbotfrac=0.0,
    rndseed=0,
    repcount=10,
    export_every_n=50_000,
    export_experiment=true,
    aggregated=1,
    keep_rawnet=false
)

max_stars_repo_name: JohannesNakayama/axelrodnet

max_issues_repo_name: JohannesNakayama/axelrodnet

max_forks_repo_name: JohannesNakayama/axelrodnet


text:
module Test

using FactCheck, PyPlot, Distributions, StatsFuns, JLD, HDF5, DataFrames, Copulas, NLopt, BlackBoxOptim

include("test.jl")

FactCheck.exitstatus()

end


text:
#
## Block methods
#

function call(b::Sidef_Types_Block_Block, args...)
    (b.value)(args...)
end

function *(b::Sidef_Types_Block_Block, n::Sidef_Types_Number_Number)
    b = b.value
    for i = 1:Int(n.value)
        b(Sidef_Types_Number_Number(i))
    end
end

max_stars_repo_name: trizen/Sidef-Deparse-Julia

max_issues_repo_name: trizen/Sidef-Deparse-Julia

max_forks_repo_name: trizen/Sidef-Deparse-Julia


text:
module SimplePolynomials
using Mods, Primes, Multisets

import Base: getindex, (==), show, zero, one, eltype, adjoint
import Base: numerator, denominator, big, inv

import Polynomials: degree, Polynomial, coeffs, roots, derivative

export SimplePolynomial, degree, coeffs, getx, Polynomial, roots
export derivative, integral, lead, make_function

# IntegerX is any sort of real or Gaussian integer
IntegerX = Union{S,Complex{S},AbstractMod} where {S<:Integer}

# RationalX is a Rational or Complex Rational based on integers
RationalX = Union{Rational{S},Complex{Rational{S}}} where {S<:Integer}

# CoefX is the type we allow for coefficients ... all exact
CoefX = Union{IntegerX,RationalX}

_enlarge(x::Integer) = big(x)
_enlarge(x::Rational) = big(x)
_enlarge(x::Complex) = big(x)
_enlarge(x::AbstractMod) = x

struct SimplePolynomial
    data::Vector
    func::Function
    function SimplePolynomial(list::Vector{T}) where {T<:CoefX}
        list = _enlarge.(_chomp(list))
        f = x -> @evalpoly(x, list...)
        new(_chomp(list), f)
    end
end

function SimplePolynomial(c...)
    n = length(c)
    T = typeof(sum(c))
    data = zeros(T, n)
    for j = 1:n
        data[j] = c[j]
    end
    SimplePolynomial(data)
end

"""
`make_function(p::SimplePolynomial)` or
`make_function(f::SimpleRationalFunction)`
returns a Julia `Function` that evaluates `p` or `f`.
Of course one can use `p(x)` or `f(x)`, but the function returned
can then be passed to (for example) `plot`.

Please note that `p(x)` [or `f(x)`] is already quite efficient.
"""
make_function(p::SimplePolynomial) = p.func

SimplePolynomial() = SimplePolynomial([0])
SimplePolynomial(p::SimplePolynomial) = SimplePolynomial(p.data)

# equality checking

(==)(p::SimplePolynomial, q::SimplePolynomial) = p.data == q.data
(==)(p::SimplePolynomial, a::T) where {T} = p == SimplePolynomial(a)
(==)(a::T, p::SimplePolynomial) where {T} = p == SimplePolynomial(a)

coeffs(p::SimplePolynomial) = copy(p.data)
eltype(p::SimplePolynomial) = eltype(p.data)

"""
`lead(p::SimplePolynomial)` returns the coefficient of the highest
power of `x` in the polynomial `p`. This is nonzero unless `p` is
the zero polynomial.
"""
lead(p::SimplePolynomial) = p.data[end]

# conversion to/from Polynomial type
SimplePolynomial(P::Polynomial) = SimplePolynomial(P.coeffs)
Polynomial(p::SimplePolynomial) = Polynomial(small.(p.data))

"""
`_chomp(list)` creates a new list by removing all trailing zeros
unless that would remove everything, in which case we leave a single zero.
"""
function _chomp(list::Vector{T})::Vector{T} where {T<:Number}
    if length(list) == 0 || all(list .== 0)
        return T[0]
    end
    k = findlast(list .!= 0)
    return list[1:k]
end

"""
`_chomp!(list)` modifies list by removing all trailing zeros, unless
the list is entirely zeros in which case we leave a single 0.
"""
function _chomp!(list::Vector{T}) where {T<:Number}
    while length(list) > 1
        n = length(list)
        if list[n] != 0
            return
        end
        deleteat!(list, n)
    end
end

"""
`degree(p::SimplePolynomial)` returns the degree of `p`; that is,
the highest exponent with a nonzero coefficient. If `p`
is the zero polynomial, then `-1` is returned.
"""
function degree(p::SimplePolynomial)
    n = length(p.data)
    if n > 1 || p.data[1] != 0
        return n - 1
    end
    return -1
end

"""
For a `SimplePolynomial`, `p`, use `p[k]` to read the coefficient
of `x^k`. In particular, `p[0]` is the constant term. If `k` is
larger than the degree of the polynomial, `0` is returned.
"""
function getindex(p::SimplePolynomial, k::Int)
    n = length(p.data)
    if k < 0
        error("index [$k] must be nonnegative")
    end

    if k >= n
        return 0
    end

    return p.data[k+1]
end

# This implements evaluation
function (p::SimplePolynomial)(x)
    try
        return p.func(x)
    catch
    end
    result = 0
    n = degree(p)
    for j = n:-1:0
        result *= x
        result += p[j]
    end
    return result
end

# special case when arg is a matrix
function (p::SimplePolynomial)(x::AbstractMatrix{T}) where {T}
    r, c = size(x)
    if r != c
        error("number of rows $r must equal the number of columns $c")
    end
    S = eltype(p)
    ST = promote_type(S, T)
    id = zeros(ST, r, r)
    for i = 1:r
        id[i, i] = 1
    end
    result = zeros(ST, r, r)
    n = degree(p)
    for j = n:-1:0
        result *= x
        result += p[j] * id
    end
    return result
end

"""
`getx()` is a convenient way to create `SimplePolynomial(0,1)`.
Typical use:
```
x = getx()
p = 2 + x - 3*x^2
```
Use `getx(T)` for coefficients of type `T`, e.g., `getx(Mod{13})`.
"""
function getx(T = Int)
    return SimplePolynomial(T(0), T(1))
end

zero(::Type{SimplePolynomial}) = SimplePolynomial(zero(Int))
one(::Type{SimplePolynomial}) = SimplePolynomial(one(Int))

roots(p::SimplePolynomial) = roots(Polynomial(p))

"""
`derivative(P::SimplePolynomial)` returns the derivative of `P`.
May also be found as `P'`.
"""
function derivative(P::SimplePolynomial)
    if degree(P) < 1
        return SimplePolynomial(0)
    end
    data = [k * P[k] for k = 1:degree(P)]
    return SimplePolynomial(data)
end
adjoint(P::SimplePolynomial) = derivative(P)

"""
`integral(P::SimplePolynomial)` returns the integral of `P`
with constant term equal to zero.
"""
function integral(P::SimplePolynomial)
    if P == 0
        return P
    end
    data = [P[k-1] // k for k = 1:degree(P)+1]
    prepend!(data, 0)
    return SimplePolynomial(data)
end

include("SimpleRationalFunctions.jl")
include("show.jl")
include("arithmetic.jl")
include("rat_arith.jl")
include("gcd.jl")
include("small.jl")
include("rational_roots.jl")
include("newton.jl")

end


text:
using CartesianGP
using Base.Test

# Half adder

numinputs0 = 2
numoutputs0 = 2
numperlevel0 = 2
numlevels0 = 1
numlevelsback0 = 1

f0 = default_funcs()
p0 = Parameters(numinputs0, numoutputs0, numperlevel0, numlevels0, numlevelsback0)
c0 = Chromosome(p0)

c0.inputs = [InputNode(1), InputNode(2)]
c0.interiors = [InteriorNode(XOR, [(0, 1), (0, 2)]) InteriorNode(AND, [(0, 1), (0, 2)])]
c0.outputs = [OutputNode((1, 1)), OutputNode((1, 2))]

@test execute_chromosome(c0) == BitString[0x6,0x8]
@test c0.number_active_nodes == 6

numinputs1 = 2
numoutputs1 = 2
numperlevel1 = 1
numlevels1 = 5
numlevelsback1 = 5

p1 = Parameters(numinputs1, numoutputs1, numperlevel1, numlevels1, numlevelsback1)
c1 = Chromosome(p1)

c1.inputs = [InputNode(1), InputNode(2)]
c1.interiors = Array(InteriorNode, p1.numlevels, p1.numperlevel)
c1.interiors[1,1] = InteriorNode(OR, [(0, 2), (0, 1)])
c1.interiors[2,1] = InteriorNode(XOR, [(0, 2), (0, 1)])
c1.interiors[3,1] = InteriorNode(OR, [(2, 1), (0, 2)])
c1.interiors[4,1] = InteriorNode(AND, [(1, 1), (0, 2)])
c1.interiors[5,1] = InteriorNode(AND, [(2, 1), (0, 2)])
c1.outputs = [OutputNode((5, 1)),OutputNode((4,1))]
@test execute_chromosome(c1) == BitString[0x2,0xa]
@test c1.number_active_nodes == 8
@test c1.interiors[1,1].active == true
@test c1.interiors[2,1].active == true
@test c1.interiors[3,1].active == false
@test c1.interiors[4,1].active == true
@test c1.interiors[5,1].active == true


text:
import Base: *

struct OneHotVector <: AbstractVector{Bool}
  ix::UInt32
  of::UInt32
end

Base.size(xs::OneHotVector) = (Int64(xs.of),)

Base.getindex(xs::OneHotVector, i::Integer) = i == xs.ix

Base.getindex(xs::OneHotVector, ::Colon) = OneHotVector(xs.ix, xs.of)

A::AbstractMatrix * b::OneHotVector = A[:, b.ix]

struct OneHotMatrix{A<:AbstractVector{OneHotVector}} <: AbstractMatrix{Bool}
  height::Int
  data::A
end

Base.size(xs::OneHotMatrix) = (Int64(xs.height),length(xs.data))

Base.getindex(xs::OneHotMatrix, i::Union{Integer, AbstractVector}, j::Integer) = xs.data[j][i]
Base.getindex(xs::OneHotMatrix, ::Colon, i::Integer) = xs.data[i]
Base.getindex(xs::OneHotMatrix, ::Colon, i::AbstractArray) = OneHotMatrix(xs.height, xs.data[i])
Base.getindex(xs::OneHotMatrix, ::Colon, ::Colon) = OneHotMatrix(xs.height, copy(xs.data))

Base.getindex(xs::OneHotMatrix, i::Integer, ::Colon) = map(x -> x[i], xs.data)

# remove workaround when https://github.com/JuliaGPU/CuArrays.jl/issues/676 is fixed
A::AbstractMatrix * B::OneHotMatrix = A[:, cpu(map(x->x.ix, B.data))]

Base.hcat(x::OneHotVector, xs::OneHotVector...) = OneHotMatrix(length(x), [x, xs...])

batch(xs::AbstractArray{<:OneHotVector}) = OneHotMatrix(length(first(xs)), xs)

import Adapt: adapt, adapt_structure

adapt_structure(T, xs::OneHotMatrix) = OneHotMatrix(xs.height, adapt(T, xs.data))

import .CUDA: CuArray, CuArrayStyle, cudaconvert
import Base.Broadcast: BroadcastStyle, ArrayStyle
BroadcastStyle(::Type{<:OneHotMatrix{<:CuArray}}) = CuArrayStyle{2}()
cudaconvert(x::OneHotMatrix{<:CuArray}) = OneHotMatrix(x.height, cudaconvert(x.data))

"""
    onehot(l, labels[, unk])

Create a `OneHotVector` with its `l`-th element `true` based on the
possible set of `labels`.
If `unk` is given, return `onehot(unk, labels)` if the input label `l` is not found
in `labels`; otherwise, it will raise an error.

# Examples
```jldoctest
julia> Flux.onehot(:b, [:a, :b, :c])
3-element Flux.OneHotVector:
 0
 1
 0

julia> Flux.onehot(:c, [:a, :b, :c])
3-element Flux.OneHotVector:
 0
 0
 1
```
"""
function onehot(l, labels)
  i = something(findfirst(isequal(l), labels), 0)
  i > 0 || error("Value $l is not in labels")
  OneHotVector(i, length(labels))
end

function onehot(l, labels, unk)
  i = something(findfirst(isequal(l), labels), 0)
  i > 0 || return onehot(unk, labels)
  OneHotVector(i, length(labels))
end

"""
    onehotbatch(ls, labels[, unk...])

Create a `OneHotMatrix` with a batch of labels based on the
possible set of `labels`.
If `unk` is given, return [`onehot(unk, labels)`](@ref) if one of the input
labels `ls` is not found in `labels`; otherwise it will error.

# Examples
```jldoctest
julia> Flux.onehotbatch([:b, :a, :b], [:a, :b, :c])
3×3 Flux.OneHotMatrix{Array{Flux.OneHotVector,1}}:
 0  1  0
 1  0  1
 0  0  0
```
"""
onehotbatch(ls, labels, unk...) =
  OneHotMatrix(length(labels), [onehot(l, labels, unk...) for l in ls])

Base.argmax(xs::OneHotVector) = xs.ix

"""
    onecold(y[, labels = 1:length(y)])

Inverse operations of [`onehot`](@ref).

# Examples
```jldoctest
julia> Flux.onecold([true, false, false], [:a, :b, :c])
:a

julia> Flux.onecold([0.3, 0.2, 0.5], [:a, :b, :c])
:c
```
"""
onecold(y::AbstractVector, labels = 1:length(y)) = labels[Base.argmax(y)]

onecold(y::AbstractMatrix, labels...) =
  dropdims(mapslices(y -> onecold(y, labels...), y, dims=1), dims=1)

onecold(y::OneHotMatrix, labels...) = map(x -> Flux.onecold(x, labels...), y.data)

@nograd onecold, onehot, onehotbatch


text:
using Distributions
using StatsBase

flip(bit::Bool) = ~bit
function flip!(individual::BitArray, index::Integer)
    individual[index] = ~individual[index]
end

abstract type Mutation end

struct UniformlyIndependentMutation <: Mutation
    probability::Real
end

function uiimutationpositions(x::Instance, probability::Real)
    if probability == 0
        return []
    elseif probability == 1
        return collect(1:length(x))
    end

    positions::Array{Integer,1} = []
    geom = Geometric(probability)
    bit = rand(geom) + 1

    while bit ≤ length(x)
        push!(positions, bit)
        bit += rand(geom) + 1
    end
    positions
end

function mutation(x::Instance, mut::UniformlyIndependentMutation)
    positions = uiimutationpositions(x, mut.probability)
    length(positions) == 0 && return x

    y = copy(x.individual)
    for bit ∈ positions
        flip!(y, bit)
    end
    Instance(y, x.problem, fitnessvalue = fitness(x, positions))
end

struct KBitFlip <: Mutation
    K::Integer
end

function mutation(x::Instance, mut::KBitFlip)
    positions = sample(1:length(x), mut.K, replace = false)
    y = copy(x.individual)
    for bit ∈ positions
        flip!(y, bit)
    end
    Instance(y, x.problem, fitnessvalue = fitness(x, positions))
end

function discretepowerlaw(β::Real, n)
    if β ≤ 1
        error("β must be greater than 1")
    end
    n2::Integer = floor(n)

    pdf = collect(Float64, 1:n2) .^ (-β)

    pdf = pdf ./ sum(pdf)

    Categorical(pdf)
end

struct HeavyTailedMutation <: Mutation
    β::Float64
    n::Integer
    dpl::DiscreteNonParametric{Integer,Float64,Base.OneTo{Integer},Array{Float64,1}}
    function HeavyTailedMutation(β::Float64, n::Integer; dpl = discretepowerlaw(β, n))
        new(β, n, dpl)
    end
end

function mutation(x::Instance, mut::HeavyTailedMutation)
    positions = uiimutationpositions(x, rand(mut.dpl) // mut.n)
    length(positions) == 0 && return x

    y = copy(x.individual)
    for bit ∈ positions
        flip!(y, bit)
    end
    Instance(y, x.problem, fitnessvalue = fitness(x, positions))
end

max_stars_repo_name: DTUComputeTONIA/SDRLS

max_issues_repo_name: DTUComputeTONIA/SDRLS

max_forks_repo_name: DTUComputeTONIA/SDRLS


text:
module ZeroInflatedDistributions

using Random
using Statistics

using Distributions
using StatsFuns

export AbstractZeroInflatedLink,
    LogitLogLink,
    PoissonLink,
    IdentityLink,
    encprob,
    posrate,
    ZeroInflatedDistribution

include("zeroinflated-link-funs.jl")
include("zeroinflated-distribution.jl")
end # module


text:
import ImageView
using Colors
using TestImages

c = ImageView.canvasgrid(2,2)
ImageView.view(c[1,1], testimage("lighthouse"), pixelspacing=[1,1])
ImageView.view(c[1,2], testimage("mountainstream"), pixelspacing=[1,1])
ImageView.view(c[2,1], testimage("moonsurface"), pixelspacing=[1,1])
ImageView.view(c[2,2], testimage("mandrill"), pixelspacing=[1,1])