text:
using DifferentialEquations
using SimpleDiffEq
using DiffEqSensitivity
using Random
using Distributions
using NestedSamplers
using StatsBase: sample, Weights
using MCMCChains: Chains, describe
using StatsPlots

function sir_ode!(du,u,p,t)
    (S,I,R,C) = u
    (β,c,γ) = p
    N = S+I+R
    infection = β*c*I/N*S
    recovery = γ*I
    @inbounds begin
        du[1] = -infection
        du[2] = infection - recovery
        du[3] = recovery
        du[4] = infection
    end
    nothing
end;

δt = 1.0
tmax = 40.0
tspan = (0.0,tmax)
obstimes = 1.0:1.0:tmax;

u0 = [990.0,10.0,0.0,0.0]; # S,I.R,Y

p = [0.05,10.0,0.25]; # β,c,γ

prob_ode = ODEProblem(sir_ode!,u0,tspan,p)
sol_ode = solve(prob_ode,Tsit5(),saveat=δt);

out = Array(sol_ode)
C = out[4,:];

X = C[2:end] .- C[1:(end-1)];

Random.seed!(1234);

Y = rand.(Poisson.(X));

function ll(x)
    (i0,β) = x
    I = i0*1000.0
    prob = remake(prob_ode,u0=[1000.0-I,I,0.0,0.0],p=[β,10.0,0.25])
    sol = solve(prob,Tsit5(),saveat=δt)
    out = Array(sol)
    C = out[4,:]
    X = C[2:end] .- C[1:(end-1)]
    nonpos = sum(X .<= 0)
    if nonpos > 0
        return Inf
    end
    sum(logpdf.(Poisson.(X),Y))
end;

priors = [
    Uniform(0, 0.1),
    Uniform(0, 0.1)
]

model = NestedModel(ll, priors);

spl = Nested(2, 10000, bounds=Bounds.MultiEllipsoid)

chain = sample(model, spl;
               param_names=["i0", "β"],
               chain_type=Chains)

describe(chain)

plot(chain)

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

max_stars_repo_name: snowdj/sir-julia

max_issues_repo_name: snowdj/sir-julia

max_forks_repo_name: snowdj/sir-julia


text:
export SSparseMatrixCSC

@static if VERSION < v"1.4.0"

    """
        SSparseMatrixCSC{Tv,Ti<:Integer, NNZ, NP} <: AbstractSparseMatrix{Tv,Ti}

    static version of SparseMatrixCSC
    """
    struct SSparseMatrixCSC{Tv,Ti<:Integer,NNZ,NP} <: AbstractSparseMatrix{Tv,Ti}
        m::Int                  # Number of rows
        n::Int                  # Number of columns
        colptr::SVector{NP,Ti}      # Column i is in colptr[i]:(colptr[i+1]-1)
        rowval::SVector{NNZ,Ti}      # Row values of nonzeros
        nzval::SVector{NNZ,Tv}       # Nonzero values

        function SSparseMatrixCSC{Tv,Ti,NNZ,NP}(
            m::Integer,
            n::Integer,
            colptr::SVector{NP,Ti},
            rowval::SVector{NNZ,Ti},
            nzval::SVector{NNZ,Tv},
        ) where {Tv,Ti<:Integer,NNZ,NP}
            m < 0 && throw(ArgumentError("number of rows (m) must be ≥ 0, got $m"))
            n < 0 && throw(ArgumentError("number of columns (n) must be ≥ 0, got $n"))
            new(Int(m), Int(n), colptr, rowval, nzval)
        end
    end

else
    # NOTE: from 1.4.0, by subtyping AbstractSparseMatrixCSC, things like sparse broadcast
    # should just work.

    """
        SSparseMatrixCSC{Tv,Ti<:Integer, NNZ, NP} <: AbstractSparseMatrix{Tv,Ti}

    static version of SparseMatrixCSC
    """
    struct SSparseMatrixCSC{Tv,Ti<:Integer,NNZ,NP} <:
           SparseArrays.AbstractSparseMatrixCSC{Tv,Ti}
        m::Int                  # Number of rows
        n::Int                  # Number of columns
        colptr::SVector{NP,Ti}      # Column i is in colptr[i]:(colptr[i+1]-1)
        rowval::SVector{NNZ,Ti}      # Row values of nonzeros
        nzval::SVector{NNZ,Tv}       # Nonzero values

        function SSparseMatrixCSC{Tv,Ti,NNZ,NP}(
            m::Integer,
            n::Integer,
            colptr::SVector{NP,Ti},
            rowval::SVector{NNZ,Ti},
            nzval::SVector{NNZ,Tv},
        ) where {Tv,Ti<:Integer,NNZ,NP}
            m < 0 && throw(ArgumentError("number of rows (m) must be ≥ 0, got $m"))
            n < 0 && throw(ArgumentError("number of columns (n) must be ≥ 0, got $n"))
            new(Int(m), Int(n), colptr, rowval, nzval)
        end
    end
    SparseArrays.getcolptr(M::SSparseMatrixCSC) = M.colptr
    SparseArrays.rowvals(M::SSparseMatrixCSC) = M.rowval
end # @static

function SSparseMatrixCSC(
    m::Integer,
    n::Integer,
    colptr::SVector,
    rowval::SVector,
    nzval::SVector,
)
    Tv = eltype(nzval)
    Ti = promote_type(eltype(colptr), eltype(rowval))
    SSparseMatrixCSC{Tv,Ti,length(nzval),n + 1}(m, n, colptr, rowval, nzval)
end

function size(spm::SSparseMatrixCSC{Tv,Ti,NNZ,NP}, i::Integer) where {Tv,Ti,NNZ,NP}
    i == 1 ? spm.m : (i == 2 ? NP - 1 : throw(ArgumentError("dimension out of bound!")))
end
size(spm::SSparseMatrixCSC{Tv,Ti,NNZ,NP}) where {Tv,Ti,NNZ,NP} = (spm.m, NP - 1)

function getindex(ssp::SSparseMatrixCSC{Tv}, i, j) where {Tv}
    S = ssp.colptr[j]
    E = ssp.colptr[j+1] - 1
    for ii = S:E
        if i == ssp.rowval[ii]
            return ssp.nzval[ii]
        end
    end
    return Tv(0)
end

SparseArrays.issparse(::SSparseMatrixCSC) = true
SparseArrays.nonzeros(M::SSparseMatrixCSC) = M.nzval
SparseArrays.nnz(spm::SSparseMatrixCSC{Tv,Ti,NNZ}) where {Tv,Ti,NNZ} = NNZ
function SparseArrays.findnz(S::SSparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
    numnz = nnz(S)
    I = Vector{Ti}(undef, numnz)
    J = Vector{Ti}(undef, numnz)
    V = Vector{Tv}(undef, numnz)

    count = 1
    @inbounds for col = 1:S.n, k = S.colptr[col]:(S.colptr[col+1]-1)
        I[count] = S.rowval[k]
        J[count] = col
        V[count] = S.nzval[k]
        count += 1
    end

    return (I, J, V)
end
SparseArrays.dropzeros!(M::SSparseMatrixCSC; trim::Bool = false) = M

ext: jl


text:
using CStructures
using Test

@testset "CStructures" begin
    include("cstruct.jl")
    include("layout.jl")
end


text:
module BytePairEncoding

using WordTokenizers
using InternedStrings

export normalize, UtfNormalizer
export BPELearner, learn!, add!, emit
export set_endsym, set_tokenizer, tokenize, whitespace_tokenize
export isolate_gloss
export Bpe, process_line, segment, segment_token

include("./utfnorm.jl")
include("./stats.jl")
include("./learn.jl")
include("./glossary.jl")
include("./bpe.jl")

end # module


text:
@testset "BoolVarView.jl" begin
    @testset "BoolVarViewNot" begin
        @testset "isempty()" begin
            trailer = SeaPearl.Trailer()
            orig = SeaPearl.BoolDomain(trailer)
            dom = SeaPearl.BoolDomainViewNot(orig)

            @test !isempty(dom)

            SeaPearl.remove!(orig, false)
            SeaPearl.remove!(orig, true)
            @test isempty(dom)

            trailer = SeaPearl.Trailer()
            orig = SeaPearl.BoolDomain(trailer)
            dom = SeaPearl.BoolDomainViewNot(orig)

            @test !isempty(dom)

            SeaPearl.remove!(dom, false)
            SeaPearl.remove!(dom, true)
            @test isempty(dom)
        end

        @testset "length()" begin
            trailer = SeaPearl.Trailer()
            orig = SeaPearl.BoolDomain(trailer)
            dom = SeaPearl.BoolDomainViewNot(orig)
            @test length(dom) == 2

            SeaPearl.remove!(dom, false)
            @test length(dom) == 1

            SeaPearl.remove!(dom, true)
            @test length(dom) == 0

            trailer = SeaPearl.Trailer()
            orig = SeaPearl.BoolDomain(trailer)
            dom = SeaPearl.BoolDomainViewNot(orig)
            @test length(dom) == 2

            SeaPearl.remove!(orig, false)
            @test length(dom) == 1

            SeaPearl.remove!(orig, true)
            @test length(dom) == 0
        end

        @testset "in()" begin
            trailer = SeaPearl.Trailer()
            orig = SeaPearl.BoolDomain(trailer)
            dom = SeaPearl.BoolDomainViewNot(orig)

            @test true in dom
            @test false in dom
            SeaPearl.remove!(orig, false)
            @test !(true in dom)
            @test false in dom
        end

        @testset "remove!()" begin
            trailer = SeaPearl.Trailer()
            orig = SeaPearl.BoolDomain(trailer)
            dom = SeaPearl.BoolDomainViewNot(orig)

            SeaPearl.remove!(dom, false)

            @test !(false in dom)
            @test length(dom) == 1
        end

        @testset "assign!()" begin
            trailer = SeaPearl.Trailer()
            orig = SeaPearl.BoolDomain(trailer)
            dom = SeaPearl.BoolDomainViewNot(orig)

            removed = SeaPearl.assign!(dom, true)

            @test !(false in dom)
            @test true in dom
            @test length(dom) == 1
            @test sort(removed) == [false]
            @test false in orig
            @test !(true in orig)
        end

        @testset "iterate()" begin
            trailer = SeaPearl.Trailer()
            orig = SeaPearl.BoolDomain(trailer)
            dom = SeaPearl.BoolDomainViewNot(orig)

            j = 1
            for i in dom
                if i == 1
                    @test i
                else
                    @test !i
                end
                j += 1
            end
            @test j == 3
        end

        @testset "removeAll!()" begin
            trailer = SeaPearl.Trailer()
            orig = SeaPearl.BoolDomain(trailer)
            dom = SeaPearl.BoolDomainViewNot(orig)

            removed = SeaPearl.removeAll!(dom)

            @test isempty(dom)
            @test removed == [true, false]
        end
    end
end


text:
#=
Reconstruct "brain" from Fourier measurements.
The measurements are generated with the Matlab files from
http://bigwww.epfl.ch/algorithms/mriphantom
=#

# ------------------------------------------------------------
# Load data

using JLD
data_file = joinpath(dirname(@__FILE__), "brain512.jld")
F = load( data_file )
f = F["freqs"]

# Sampling locations
M = size( f )
J = Int( log2(M[1]) ) - 1

# ------------------------------------------------------------
# Compute reconstruction

using GeneralizedSampling

xi = GeneralizedSampling.grid(M)
T = Freq2Wave( xi, "db4", J )
w = T \ f

# ------------------------------------------------------------
# Plot reconstruction

using IntervalWavelets, Winston

y = weval( abs(w), "db4", 10 )

miny = minimum(y)
maxy = maximum(y)
imagesc( y, (miny,maxy) )


text:
# the julia iteration function

function julia(z, maxiter::Int64)
    c=-0.8+0.156im
    for n = 1:maxiter
        if abs(z) > 2
            return n-1
        end
        z = z^2 + c
    end
    return maxiter
end

max_stars_repo_name: rabarar/juliaFractals

max_issues_repo_name: rabarar/juliaFractals

max_forks_repo_name: rabarar/juliaFractals


text:
using Telegrambot
using Test
@test true


text:
# Some thresholds, implemented as types so we can save them with JLD2
# (which does not allow saving functions directly)

# Threshold from Optimal Best Arm Identification with Fixed Confidence
# (Garivier and Kaufmann 2016). Recommended in section 6.
struct GK16
    δ;
end

function (β::GK16)(t)
    log((log(t)+1)/β.δ)
end


text:
# This is a test file to make codecov happy, technically all of the
# tests here are already done / integrated within other tests.

@testset "strings" begin
    st = "blah"

    @test F.str(st) == "blah"

    sst = SubString("blahblah", 1:4)
    @test sst == "blah"
    @test F.str(sst) == "blahblah"

    sst = SubString("blah✅💕and etcσ⭒ but ∃⫙∀ done", 1:27)
    @test F.to(sst) == 27

    s = "aabccabcdefabcg"
    for m ∈ eachmatch(r"abc", s)
        @test s[F.matchrange(m)] == "abc"
    end
end

@testset "ocblock" begin
    st = "This is a block  and done"
    τ = F.find_tokens(st, F.MD_TOKENS, F.MD_1C_TOKENS)
    ocb = F.OCBlock(:COMMENT, (τ[1]=>τ[2]))
    @test F.otok(ocb) == τ[1]
    @test F.ctok(ocb) == τ[2]
end

@testset "isexactly" begin
    steps, b, λ = F.isexactly("


text:
using Test
using LinearAlgebra
using DataFlowTasks
using DataFlowTasks: data_dependency,R,W,RW

@testset "Macros" begin
    m,n  = 100,100
    A    = rand(n,m)
    b    = 2
    # verbose way of defining task
    t    = @dtask rmul!(A,b) (A,) (RW,)
    tmp  = A*b
    @test A !== tmp
    schedule(t)
    wait(t)
    @test A == tmp
end

@testset "Data dependency" begin
    m,n  = 100,100
    A    = rand(n,m)
    a    = -π
    b    = 2
    t1 = @dtask rmul!(A,b) (A,b) (RW,R)
    t2 = @dtask rmul!(A,b) (A,b) (RW,R)
    @test data_dependency(t1,t2) == true
    # rebind A to a different matrix, and make sure independency is inferred
    A = rand(100,100)
    t3 = @dtask rmul!(A,b) (A,b) (RW,R)
    @test data_dependency(t1,t3) == false
end


text:
if haskey(ENV, "CI")
    ENV["PLOTS_TEST"] = "true"
    ENV["GKSwstype"] = "100" # gr segfault workaround
end
using Plots
using RobustAndOptimalControl
using LinearAlgebra
using Test

@testset "RobustAndOptimalControl.jl" begin

    @testset "extendedstatespace" begin
        @info "Testing extendedstatespace"
        include("test_extendedstatespace.jl")
    end

    @testset "utils" begin
        @info "Testing utils"
        include("test_utils.jl")
    end

    @testset "descriptor" begin
        @info "Testing descriptor"
        include("test_descriptor.jl")
    end

    @testset "uncertainty" begin
        @info "Testing uncertainty"
        include("test_uncertainty.jl")
    end

    @testset "diskmargin" begin
        @info "Testing diskmargin"
        include("test_diskmargin.jl")
    end

    @testset "weights" begin
        @info "Testing weights"
        include("test_weights.jl")
    end

    @testset "hinfpartition" begin
        @info "Testing hinfpartition"
        include("test_hinfpartition.jl")
    end

    @testset "H∞ design" begin
        @info "Testing hinf_design"
        include("test_hinf_design.jl")
    end

    @testset "H2 design" begin
        @info "Testing H2 design"
        include("test_h2_design.jl")
    end

    @testset "LQG" begin
        @info "Testing LQG"
        include("test_lqg.jl")
    end

    @testset "Named systems" begin
        @info "Testing Named systems"
        include("test_named_systems2.jl")
    end

    @testset "find_lft" begin
        @info "Testing find_lft"
        include("test_find_lft.jl")
    end

    @testset "reduction" begin
        @info "Testing test_reduction"
        include("test_reduction.jl")
    end

    @testset "augmentation" begin
        @info "Testing augmentation"
        include("test_augmentation.jl")
    end

    @testset "glover_mcfarlane" begin
        @info "Testing glover_mcfarlane"
        include("test_glover_mcfarlane.jl")
    end

    @testset "hinfgrad" begin
        @info "Testing hinfgrad"
        include("test_hinfgrad.jl")
    end

end


text:
"""
Essential tools used throughout Iris.

Submodules:
  * [`Boundaries`](@ref)
  * [`BoundaryLayers`](@ref)
  * [`BoundaryConditions`](@ref)
  * [`Curlcurls`](@ref)
  * [`DielectricFunctions`](@ref)
  * [`Dispersions`](@ref)
  * [`Domains`](@ref)
  * [`ElectricFields`](@ref)
  * [`Lattices`](@ref)
  * [`LU_Factorizations`](@ref)
  * [`Plotting`](@ref)
  * [`Points`](@ref)
  * [`PumpFunctions`](@ref)
  * [`SelfEnergies`](@ref)
  * [`Shapes`](@ref)
  * [`Simulations`](@ref)
"""
module Common

"""
	Symmetric
"""
struct Symmetric end

"""
	Unsymmetric
"""
struct Unsymmetric end

# `Defaults`: contains all default parameters in one location for uniformity and easy alteration.
include("../Defaults.jl")

include("Points.jl") # ✅
using .Points
export Point
export Cartesian
export Polar
export Spherical

include("VectorFields.jl") # ✅
using .VectorFields
export ScalarField
export ElectricField
export VectorField
export update!
export component

include("Shapes.jl") # ✅
using .Shapes
export AbstractShape
export Interval
export Circle
export Square
export Rectangle

include("BoundaryLayers.jl") # ✅
using .BoundaryLayers
export PML
export cPML
export noBL

include("BoundaryConditions.jl") # ✅
using .BoundaryConditions
export noBC
export DirichletBC
export NeumannBC
export FloquetBC
export MatchedBC
export BCHermiticity
export HermitianBC
export NonHermitianBC

include("Boundaries.jl") # ✅
using .Boundaries
export Boundary

include("DielectricFunctions.jl") # ✅
using .DielectricFunctions
export AbstractDielectricFunction

include("PumpFunctions.jl") # ✅
using .PumpFunctions
export AbstractPumpFunction

include("Dispersions.jl") # ✅
using .Dispersions
export NoDispersion
export TwoLevelSystem
export susceptability

include("Lattices.jl") # ✅ Except for Bravais
using .Lattices
export Lattices
export Lattice
export latticeindex

include("Domains.jl") # ✅ Only 1D, 2D Cartesian
using .Domains
export Domains
export LatticeDomain
export NondispersiveDomain
export DispersiveDomain

include("Laplacians.jl") # ✅ Only 1D, 2D Cartesian
using .Laplacians
export Laplacian

include("Curlcurls.jl")
using .Curlcurls
export Curlcurl

# # include("Curls.jl")
# # using .Curls
# # export Curl

include("SelfEnergies.jl") # ✅ Only 1D, 2D Cartesian Symmetric Dirichlet & Open, functions need work
using .SelfEnergies

include("Simulations.jl") # ✅ Only 1D, 2D Cartesian Symmetric sans Curlcurls and no smoothing in 2D
using .Simulations
export Simulation
export smooth!
export update_dielectric!
export update_pump!

include("HelmholtzOperators.jl") # ✅ Only 1D, 2D Cartesian Symmetric
using .HelmholtzOperators
export Helmholtz

include("MaxwellOperators.jl")
using .MaxwellOperators
export Maxwell

include("LU_Factorizations.jl") # ✅ checked for UMFPACK, not for Pardiso or MUMPS
using .LU_Factorizations
export DEFAULT_LUPACK
export MSolver
export PSolver
export USolver
export AbstractLUPACK

include("Plotting.jl") # ✅ basic 2D and 1D simulation plotting
using .Plotting

end # module


text:
# module moore2018multipropopt

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

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

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

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

printiter = 1
#TODO verify this works with parallel

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

# # include("$path/../../BeamFEA.jl/src/BeamFEA.jl")
# rc("figure", figsize=(5.8, 2.6))
# rc("font", size=10.0)
# rc("lines", linewidth=1.5)
# rc("lines", markersize=3.0)
# rc("legend", frameon=false)
# rc("axes.spines", right=false, top=false)
# rc("figure.subplot", left=0.17, bottom=0.18, top=0.97, right=.69)
# rc("axes", color_cycle=["348ABD", "A60628", "009E73", "7A68A6", "D55E00", "CC79A7"])
# plot_cycle=["#348ABD", "#A60628", "#009E73", "#7A68A6", "#D55E00", "#CC79A7"]
# close("all")

# function optimize(numprops,to_factor,db_factor)

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

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

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

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

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

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

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

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

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

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

    plyprops = plyproperties()

    matnames = plyprops.names

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

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

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

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

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

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

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

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

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

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

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

    )

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

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

    p_lam_t = [0.00296032, 0.0015862, 0.000162432, 0.000559344, 0.000218, 0.000218, 0.000218, 0.000218]
    p_chord = [0.152775, 0.209865, 0.0599548, 0.0146912]
    p_twist_d = [43.7258, 36.3345, 24.8117, 15.5926]
    p_pitch_d_accel = 1.4774821755019851
    p_pitch_d_climb = 0.07394619362544558
    p_pitch_d_cruise = 1.2375953673499482
    velocity_accel = 19.472605774436627
    velocity_climb = 19.727538432606607
    velocity_cruise = 44.93261088033517
    p_Rtip = 0.7373054880998322
    p_rpm_accel = 1521.2024594724965
    p_rpm_climb = 1662.377310410668
    p_rpm_cruise = 2030.4950665640044
    number_blades = 2.0#2.3818948496117045
    kv = 9.01975436914855
    i0 = 5.942697466865178
    w_aoa_d_accel = 26.783254994297305
    w_aoa_d_climb = 10.682348122387333
    w_aoa_d_cruise = 1.9046025617909106
    batt_mass_accel = 0.01
    batt_mass_climb = 0.44404948094245933
    batt_mass_cruise = 119.39074611391598

    p_lam_t = [0.00336306, 0.0024627, 0.00221994, 0.000152, 0.000218, 0.000218, 0.000218, 0.000218]
    p_chord = [0.155659, 0.201413, 0.0971402, 0.0100891]
    p_twist_d = [45.0, 33.5675, 22.6151, 14.0563]
    p_pitch_d_accel = 1.4829935332414101
    p_pitch_d_climb = -7.026136800337877
    p_pitch_d_cruise = 1.2383248802411353
    velocity_accel = 19.24363792774069
    velocity_climb = 12.153346213425108
    velocity_cruise = 44.81388239980551
    p_Rtip = 0.7782175243251072
    p_rpm_accel = 1543.53387273026
    p_rpm_climb = 2019.2972133656845
    p_rpm_cruise = 2030.550389292149
    kv = 7.219833930609759
    i0 = 6.0
    w_aoa_d_accel = 26.8000091392822
    w_aoa_d_climb = 18.92926786023156
    w_aoa_d_cruise = 1.9213751825354857
    batt_mass_accel = 0.01
    batt_mass_climb = 0.2641049737151147
    batt_mass_cruise = 120.95934410306408

    p_lam_t = [0.00265251, 0.000843842, 0.000152, 0.000152, 0.000218, 0.000218, 0.000218, 0.000218]
    p_chord = [0.155659, 0.220257, 0.0117648, 0.0100891]
    p_twist_d = [45.0, 40.023, 28.6333, 1.88758]
    p_pitch_d_accel = 1.4829935332414101
    p_pitch_d_climb = -1.635119122794496
    p_pitch_d_cruise = -2.9277256334203314
    velocity_accel = 19.24363792774069
    velocity_climb = 13.643713798726896
    velocity_cruise = 44.73775558213799
    p_Rtip = 0.7118124546271428
    p_rpm_accel = 1976.2496570673752
    p_rpm_climb = 1943.6196478113399
    p_rpm_cruise = 2030.550389292149
    number_blades = 5.0#4.018472113116541
    kv = 8.155714770621731
    i0 = 5.997437563806129
    w_aoa_d_accel = 11.40529865732635
    w_aoa_d_climb = 20.23462263172518
    w_aoa_d_cruise = 1.931104840946755
    batt_mass_accel = 0.01
    batt_mass_climb = 0.16314428345511667
    batt_mass_cruise = 143.03066784675627

    p_lam_t = [0.000152, 0.000152, 0.000152, 0.000152, 0.000218, 0.000218, 0.000218, 0.000218]
    p_chord = [0.10877, 0.105743, 0.152789, 0.05]
    p_twist_d = [36.8235, 42.9044, 39.6085, 45.0]
    p_pitch_d_accel = -19.82051083412791
    p_pitch_d_climb = -20.188093719254248
    p_pitch_d_cruise = 12.113489425294997
    velocity_accel = 16.525742913704143
    velocity_climb = 20.494356096152448
    p_Rtip = 0.3352916814504692
    p_rpm_accel = 2798.6019520339005
    p_rpm_climb = 3906.0184495335675
    p_rpm_cruise = 1937.9861230380968
    kv = 7.283787837965981
    i0 = 3.246244436771027
    w_aoa_d_accel = 33.56793200298914
    w_aoa_d_climb = 17.149964223958825
    w_aoa_d_cruise = 4.6635813147838405
    batt_mass_accel = 0.001
    batt_mass_climb = 1.8299675145048462
    batt_mass_cruise = 124.65256349906039

    p_lam_t = [0.000152, 0.000152, 0.000154056, 0.000152, 0.000218, 0.000218, 0.000218, 0.000218]
    p_chord = [0.166552, 0.0975438, 0.158812, 0.05]
    p_twist_d = [36.1394, 43.9595, 39.3576, 44.9413]
    p_pitch_d_accel = -19.82051083412791
    p_pitch_d_climb = -19.673931974404663
    p_pitch_d_cruise = 12.113489425294997
    velocity_accel = 16.525742913704143
    velocity_climb = 17.92453072013977
    p_Rtip = 0.33358244398851317
    p_rpm_accel = 2798.6019520339005
    p_rpm_climb = 3982.97243760044
    p_rpm_cruise = 1937.9861230380968
    kv = 7.283787837965981
    i0 = 3.246244436771027
    w_aoa_d_accel = 33.55720131194915
    w_aoa_d_climb = 19.060521666106922
    w_aoa_d_cruise = 4.625279599698307
    batt_mass_accel = 0.001
    batt_mass_climb = 0.7513826889091813
    batt_mass_cruise = 120.56189810246892

    #Fixed aoa

    p_lam_t = [0.000152, 0.000152, 0.000153858, 0.000152, 0.000230554, 0.000218945, 0.000219017, 0.000218]
    p_chord = [0.0653822, 0.0484466, 0.0557039, 0.0841311]
    p_twist_d = [44.6798, 40.0471, 38.8732, 45.0]
    p_pitch_d_accel = -11.704791952680312
    p_pitch_d_climb = -11.40522334751054
    p_pitch_d_cruise = 14.075048962188848
    velocity_accel = 7.331401489528655
    velocity_climb = 8.03095019944203
    p_Rtip = 0.17272727272727256
    p_rpm_accel = 8003.862941356319
    p_rpm_climb = 7953.681497013595
    p_rpm_cruise = 4581.2196064094005
    kv = 15.02410595281662
    i0 = 3.2703661356139597
    w_aoa_d_accel = 30.939955827796926
    w_aoa_d_climb = 31.245140891571026
    w_aoa_d_cruise = 4.084570950676237
    batt_mass_accel = 0.001
    batt_mass_climb = 2.8935239395666583
    batt_mass_cruise = 141.3814881300007

#Fixed structure
p_lam_t = [0.000152, 0.000152, 0.000152, 0.000152, 0.000230262, 0.000218, 0.000218, 0.000218]
p_chord = [0.0520819, 0.0377408, 0.0484461, 0.036489]
p_twist_d = [45.0, 38.1963, 37.847, 35.1547]
p_pitch_d_accel = -13.405913062940945
p_pitch_d_climb = -13.527980424350305
p_pitch_d_cruise = 14.29920595639354
velocity_accel = 6.71505379850055
velocity_climb = 7.770514402557935
p_Rtip = 0.17272727272727284
p_rpm_accel = 9515.932176417387
p_rpm_climb = 9577.811551408988
p_rpm_cruise = 4840.679484070415
kv = 5.0001001868354304
i0 = 0.9035878002679382
w_aoa_d_accel = 30.905912967879324
w_aoa_d_climb = 31.289389251880156
w_aoa_d_cruise = 4.001062592416955
batt_mass_accel = 0.001
batt_mass_climb = 2.166271103964091
batt_mass_cruise = 123.44410819370431

    number_blades = 5.0

    w_halfspan = 5.7

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

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

    # sparce vars 900kg and ~50m/s vinf all
#
# xwarm = [0.01000, 0.01000, 0.01000, 0.01000, 0.01000, 0.01000, 0.01000, 0.01000, 21.62153, 7.95808, 20.34969, 5.00338, 4.33754, 4.49540, 4.49814, 3.09121, -248.53136, -23.91318, -30.31537, 132.97473, 0.24305, 33.52939, 306.22862, 45.14119, 2476.55374, 4.99968, 7.26288, 4.03376, 0.33206, 14.33827, 4.75615, 0.00100, 0.07764, 13.59082]
#     x0[1:length(xwarm)] = xwarm[1:end]

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    end

    optimize = true
    if optimize

        for i = 1:3

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

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

            x0 = xopt

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

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

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

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

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

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

    if false #sweep # sweep through the variables

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

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

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

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

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

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

    # end #optimize

    # end #module


text:
using Revise, Test
using Revise.CodeTracking
using Revise.LoweredCodeUtils

function isdefinedmod(mod::Module)
    # Not all modules---e.g., LibGit2---are reachable without loading the stdlib
    names = fullname(mod)
    pmod = Main
    for n in names
        isdefined(pmod, n) || return false
        pmod = getfield(pmod, n)
    end
    return true
end
function reljpath(path)
    for subdir in ("base/", "stdlib/", "test/")
        s = split(path, subdir)
        if length(s) == 2
            return subdir * s[end]
        end
    end
    return path
end
function filepredicate(file, reffiles)
    bfile = Base.find_source_file(file)
    bfile === nothing && return false  # when the file is "none"
    return reljpath(bfile) ∈ reffiles
end
function signature_diffs(mod::Module, signatures; filepredicate=nothing)
    extras = copy(signatures)
    modeval, modinclude = getfield(mod, :eval), getfield(mod, :include)
    failed = []
    nmethods = 0
    for fsym in names(mod; all=true)
        isdefined(mod, fsym) || continue
        f = getfield(mod, fsym)
        isa(f, Base.Callable) || continue
        (f === modeval || f === modinclude) && continue
        for m in methods(f)
            nmethods += 1
            if haskey(signatures, m.sig)
                delete!(extras, m.sig)
            else
                if filepredicate !== nothing
                    filepredicate(String(m.file)) || continue # find signatures only in selected files
                end
                push!(failed, m.sig)
            end
        end
    end
    return failed, extras, nmethods
end
function extracttype(T)
    p1 = T.parameters[1]
    isa(p1, Type) && return p1
    isa(p1, TypeVar) && return p1.ub
    error("unrecognized type ", T)
end
function in_module_or_core(T, mod::Module)
    if isa(T, TypeVar)
        return in_module_or_core(T.ub, mod)
    end
    if isa(T, UnionAll)
        T = Base.unwrap_unionall(T)
    end
    T === Union{} && return true
    if isa(T, Union)
        in_module_or_core(T.a, mod) || return false
        return in_module_or_core(T.b, mod)
    end
    Tname = T.name
    if Tname.name === :Type
        return in_module_or_core(extracttype(T), mod)
    end
    Tmod = Tname.module
    return Tmod === mod || Tmod === Core
end

module Lowering end

@testset ":lambda expressions" begin
    ex = quote
        mutable struct InnerC
            x::Int
            valid::Bool

            function InnerC(x; notvalid::Bool=false)
                return new(x, !notvalid)
            end
        end
    end
    sigs, _ = Revise.eval_with_signatures(Lowering, ex)
    @test length(sigs) >= 2
end

basefiles = Set{String}()
@time for (i, (mod, file)) in enumerate(Base._included_files)
    endswith(file, "sysimg.jl") && continue
    file = Revise.fixpath(file)
    push!(basefiles, reljpath(file))
    mexs = Revise.parse_source(file, mod)
    Revise.instantiate_sigs!(mexs; always_rethrow=true)
end
failed, extras, nmethods = signature_diffs(Base, CodeTracking.method_info; filepredicate = fn->filepredicate(fn, basefiles))
# In some cases, the above doesn't really select the file-of-origin. For example, anything
# defined with an @enum gets attributed to Enum.jl rather than the file in which @enum is used.
realfailed = similar(failed, 0)
for sig in failed
    ft = Base.unwrap_unionall(sig).parameters[1]
    match(r"^getfield\(Base, Symbol\(\"##\d", string(ft)) === nothing || continue  # exclude anonymous functions
    all(T->in_module_or_core(T, Base), Base.unwrap_unionall(sig).parameters[2:end]) || continue
    push!(realfailed, sig)
end
@test length(realfailed) < 40  # big enough for some cushion in case new "difficult" methods get added


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "blis"
version = v"0.9.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/flame/blis.git", "14c86f66b20901b60ee276da355c1b62642c18d2"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd blis/

for i in ./config/*/*.mk; do

    # Building in container forbids -march options <<< Settings overriden.
    # sed -i "s/-march[^ ]*//g" $i

    # Building in container forbids unsafe optimization.
    sed -i "s/-ffast-math//g" $i
    sed -i "s/-funsafe-math-optimizations//g" $i

done

case ${target} in

    *"x86_64"*"linux"*)
        export BLI_CONFIG=x86_64
        export BLI_THREAD=openmp
        ;;
    *"x86_64"*"w64"*)
        # MinGW doesn't support savexmm instructions
        # Build only for AMD processors.
        export BLI_CONFIG=amd64
        export BLI_THREAD=openmp
        ;;
    *"x86_64"*"apple"*)
        export BLI_CONFIG=x86_64
        export BLI_THREAD=openmp
        ;;
    *"x86_64"*"freebsd"*)
        export BLI_CONFIG=x86_64
        export BLI_THREAD=openmp
        ;;
    *"aarch64"*"apple"*)
        # Metaconfig arm64 is not needed here.
        # All Mac processors should have equal or higher specs then firestorm
        export BLI_CONFIG=firestorm
        export BLI_THREAD=openmp
        ;;
    *"aarch64"*"linux"*)
        export BLI_CONFIG=arm64
        export BLI_THREAD=openmp
        ;;
    *"arm"*"linux"*)
        export BLI_CONFIG=arm32
        export BLI_THREAD=none
        ;;
    *)
        # Default (Generic) configuration without optimized kernel.
        export BLI_CONFIG=generic
        export BLI_THREAD=none
        ;;

esac

# For 64-bit builds, add _64 suffix to exported BLAS routines.
# This corresponds to ILP64 handling of OpenBLAS thus Julia.
if [ ${nbits} = 64 ]; then
    patch frame/include/bli_macro_defs.h < ${WORKSPACE}/srcdir/patches/bli_macro_defs.h.f77suffix64.patch
fi

# Include A64FX in Arm64 metaconfig.
if [ ${BLI_CONFIG} = arm64 ]; then
    # Add A64FX to the registry.
    patch config_registry < ${WORKSPACE}/srcdir/patches/config_registry.metaconfig+a64fx.patch

    # Unscreen Arm SVE code for metaconfig.
    patch kernels/armsve/bli_kernels_armsve.h \
        < ${WORKSPACE}/srcdir/patches/armsve_kernels_unscreen_arm_sve_h.patch
    patch kernels/armsve/1m/old/bli_dpackm_armsve512_int_12xk.c \
        < ${WORKSPACE}/srcdir/patches/armsve_kernels_unscreen_arm_sve_h.patch
    patch kernels/armsve/1m/bli_dpackm_armsve256_int_8xk.c \
        < ${WORKSPACE}/srcdir/patches/armsve_kernels_unscreen_arm_sve_h.patch

    # Screen out A64FX sector cache.
    patch config/a64fx/bli_cntx_init_a64fx.c \
        < ${WORKSPACE}/srcdir/patches/a64fx_config_screen_sector_cache.patch
fi

# Import libblastrampoline-style nthreads setter.
cp ${WORKSPACE}/srcdir/nthreads64_.c frame/compat/nthreads64_.c

export BLI_F77BITS=${nbits}
./configure -p ${prefix} -t ${BLI_THREAD} -b ${BLI_F77BITS} ${BLI_CONFIG}
make -j${nproc}
make install

# Static library is not needed.
rm ${prefix}/lib/libblis.a

# Rename .dll for Windows targets.
if [[ "${target}" == *"x86_64"*"w64"* ]]; then
    mkdir -p ${libdir}
    mv ${prefix}/lib/libblis.4.dll ${libdir}/libblis.dll
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="musl"),
    Platform("armv7l", "linux"; libc="glibc"),
    Platform("x86_64", "windows"),
    Platform("x86_64", "macos"),
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("aarch64", "macos"),
    Platform("x86_64", "freebsd")
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libblis", :blis)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    # For OpenMP we use libomp from `LLVMOpenMP_jll` where we use LLVM as compiler (BSD
    # systems), and libgomp from `CompilerSupportLibraries_jll` everywhere else.
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"); platforms=filter(!Sys.isbsd, platforms)),
    Dependency(PackageSpec(name="LLVMOpenMP_jll", uuid="1d63c593-3942-5779-bab2-d838dc0a180e"); platforms=filter(Sys.isbsd, platforms)),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"11", lock_microarchitecture=false, julia_compat="1.6")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
## build_network.jl
run(pipeline(`julia ../scripts/build_network_runner.jl data/small.ChromNet.jld --quiet`, "data/small.json"))

# ## build_network.jl
# M = trues(4, ChromNet.totalBins)
# M[:,:] = rand(4, ChromNet.totalBins) .> 0.5
# header = ["A01", "B01", "C01", "D01"]
# isfile("data/small.ChromNet") && rm("data/small.ChromNet")
# save_chromnet_matrix("data/small.ChromNet", M, header)
# metadata = Dict(
#     "A01" => Dict(
#         "name" => "A",
#         "id" => "A01",
#         "description" => "A desc",
#         "cellType" => "A cellType",
#         "lab" => "A cellType",
#         "organism" => "Homo sapiens",
#         "lifeStage" => "None",
#         "treatments" => "None",
#         "antibody" => "Unknown"
#     ),
#     "B01" => Dict(
#         "name" => "B",
#         "id" => "B01",
#         "description" => "B desc",
#         "cellType" => "B cellType",
#         "lab" => "B cellType",
#         "organism" => "Homo sapiens",
#         "lifeStage" => "None",
#         "treatments" => "None",
#         "antibody" => "Unknown"
#     ),
#     "C01" => Dict(
#         "name" => "C",
#         "id" => "C01",
#         "description" => "C desc",
#         "cellType" => "C cellType",
#         "lab" => "C cellType",
#         "organism" => "Homo sapiens",
#         "lifeStage" => "None",
#         "treatments" => "None",
#         "antibody" => "Unknown"
#     ),
#     "D01" => Dict(
#         "name" => "D",
#         "id" => "D01",
#         "description" => "D desc",
#         "cellType" => "D cellType",
#         "lab" => "D cellType",
#         "organism" => "Homo sapiens",
#         "lifeStage" => "None",
#         "treatments" => "None",
#         "antibody" => "Unknown"
#     )
# )
# open(f->println(f, json(metadata)), "data/small.ChromNet/metadata", "w")
# run(pipeline(`julia ../scripts/build_network.jl data/four_bed.config --data none --quiet`, "data/small.json"))
# run(pipeline(`julia ../scripts/build_network.jl data/single_bed.config --data data/small.ChromNet --quiet`, "data/small.json"))


text:
# Per Lyndon. Needs adaptation, and/or needs redefinition of map to use functions rather
# than AbstractOp.
#function rrule(map, f, xs)
#    # Rather than 3 maps really want 1 multimap
#    ys_and_pullbacks = map(x->rrule(f, x), xs) #Take this to ys = map(f, x)
#    ys = map(first, ys_and_pullbacks)
#    pullbacks = map(last, ys_and_pullbacks)
#    function map_pullback(dys)
#        _call(f, x) = f(x)
#        dfs_and_dxs = map(_call, pullbacks, dys)
#        # but in your case you know it will be NoTangent() so can  skip
#        df = sum(first, dfs_and_dxs)
#        dxs = map(last, dfs_and_dxs)
#        return NoTangent(), df, dxs
#    end
#    return ys, map_pullback
#end


text:
immutable Erlang <: ContinuousUnivariateDistribution
    gammad::Gamma

    function Erlang(α::Real, θ::Real)
        isinteger(α) || error("Erlang shape parameter must be an integer")
        new(Gamma(α, θ))
    end

    Erlang(α::Real) = Erlang(α, 1.0)
    Erlang() = new(Gamma())
end

@distr_support Erlang 0.0 Inf

show(io::IO, d::Erlang) = ((α, β) = params(d); show_oneline(io, d, [(:α, α), (:β, β)]))

#### Parameters

shape(d::Erlang) = int(shape(d.gammad))
scale(d::Erlang) = scale(d.gammad)
rate(d::Erlang) = rate(d.gammad)
params(d::Erlang) = ((α, β) = params(d.gammad); (int(α), β))

#### Statistics

mean(d::Erlang) = mean(d.gammad)
var(d::Erlang) = var(d.gammad)
median(d::Erlang) = median(d.gammad)
skewness(d::Erlang) = skewness(d.gammad)
kurtosis(d::Erlang) = kurtosis(d.gammad)

mode(d::Erlang) = mode(d.gammad)
entropy(d::Erlang) = entropy(d.gammad)

#### Evaluation

pdf(d::Erlang, x::Float64) = pdf(d.gammad, x)
logpdf(d::Erlang, x::Float64) = logpdf(d.gammad, x)

cdf(d::Erlang, x::Float64) = cdf(d.gammad, x)
ccdf(d::Erlang, x::Float64) = ccdf(d.gammad, x)
logcdf(d::Erlang, x::Float64) = logcdf(d.gammad, x)
logccdf(d::Erlang, x::Float64) = logccdf(d.gammad, x)

quantile(d::Erlang, p::Float64) = quantile(d.gammad, p)
cquantile(d::Erlang, p::Float64) = cquantile(d.gammad, p)
invlogcdf(d::Erlang, lp::Float64) = invlogcdf(d.gammad, lp)
invlogccdf(d::Erlang, lp::Float64) = invlogccdf(d.gammad, lp)

mgf(d::Erlang, t::Real) = mgf(d.nested_gamma, t)
cf(d::Erlang, t::Real) = cf(d.nested_gamma, t)

#### Sampling

rand(d::Erlang) = rand(d.gammad)


text:
using ErlPort.Exceptions
using ErlPort.ErlTerms.Decode

include("Testing.jl")

# data for use by tests
lenheader1 = b"\0\0\0\x17"
lenheader2 = b"\0\x14"
badsizecompdata = vcat(b"\x83P\0\0\0\x16\x78\xda\xcb\x66\x10\x49\xc1\2\0\x5d",
                       b"\x60\x08\x50")
compdata1 = b"\x83P\0\0\0\x17\x78\xda\xcb\x66\x10\x49\xc1\2\0\x5d\x60\x08\x50"
compdata2 = vcat(b"\x83P\0\0\0\x17\x78\xda\xcb\x66\x10\x49\xc1\2\0\x5d\x60",
                 b"\x08\x50tail")

# tests for supporting functions
testcase() do
    @test int4unpack(lenheader1) == 23
    @test int2unpack(lenheader2) == 20
    expected = [107,0,20,100,100,100,100,100,100,100,100,100,100,
                100,100,100,100,100,100,100,100,100,100]
    @test decompressterm(compdata1) == expected
end

# basic decode errors
testcase() do
    @test_throws UnsupportedType decode("")
    @test_throws IncompleteData decode(b"")
    @test_throws UnknownProtocolVersion decode(b"\0")
    @test_throws IncompleteData decode(b"\x83")
end

# decode nil

# decode atoms
testcase() do
    @test_throws IncompleteData decode(b"\x83d")
    @test_throws IncompleteData decode(b"\x83d\0")
    @test_throws IncompleteData decode(b"\x83d\0\1")
    @test decode(b"\x83d\0\0") == (symbol(""), [])
    @test decode(b"\x83d\0\0tail") == (symbol(""), b"tail")
    @test decode(b"\x83d\0\4test") == (:test, b"")
    @test decode(b"\x83d\0\4testtail") == (:test, b"tail")
    @test decode(b"\x83d\0\4true") == (true, b"")
    @test decode(b"\x83d\0\5false") == (false, b"")
    @test decode(b"\x83d\0\x09undefined") == (nothing, b"")
    @test decodeatom(b"d\0\0") == (symbol(""), [])
    @test decodeatom(b"d\0\0tail") == (symbol(""), b"tail")
    @test decodeatom(b"d\0\4test") == (:test, b"")
    @test decodeatom(b"d\0\4testtail") == (:test, b"tail")
    @test decodeatom(b"d\0\4true") == (true, b"")
    @test decodeatom(b"d\0\5false") == (false, b"")
    @test decodeatom(b"d\0\x09undefined") == (nothing, b"")
end

# decode predefined atoms

# decode empty list

# decode string list
testcase() do
    @test_throws IncompleteData decodestring(b"")
    @test_throws IncompleteData decodestring(b"\0")
    @test_throws IncompleteData decodestring(b"\0\0")
end

# decode list

# decode improper list

# decode small tuple

# decode large tuple

# decode opaque object

# decode small integer
testcase() do
    @test_throws IncompleteData decode(b"\x83a")
    @test decode(b"\x83a\0") == (0, b"")
    @test decode(b"\x83a\0tail") == (0, b"tail")
    @test decode(b"\x83a\xff") == (255, b"")
    @test decode(b"\x83a\xfftail") == (255, b"tail")
    @test_throws IncompleteData decodesmallint(b"a")
    @test decodesmallint(b"a\0") == (0, b"")
    @test decodesmallint(b"a\0tail") == (0, b"tail")
    @test decodesmallint(b"a\xff") == (255, b"")
    @test decodesmallint(b"a\xfftail") == (255, b"tail")
end

# decode integer
testcase() do
    @test_throws IncompleteData decode(b"\x83b")
    @test_throws IncompleteData decode(b"\x83b\0")
    @test_throws IncompleteData decode(b"\x83b\0\0")
    @test_throws IncompleteData decode(b"\x83b\0\0\0")
    @test decode(b"\x83b\0\0\0\0") == (0, [])
    @test decode(b"\x83b\0\0\0\0tail") == (0, b"tail")
    @test decode(b"\x83b\x7f\xff\xff\xff") == (2147483647, [])
    @test decode(b"\x83b\x7f\xff\xff\xfftail") == (2147483647, b"tail")
    @test decode(b"\x83b\xff\xff\xff\xff") == (-1, [])
    @test decode(b"\x83b\xff\xff\xff\xfftail") == (-1, b"tail")
    @test_throws IncompleteData decodeint(b"b")
    @test_throws IncompleteData decodeint(b"b\0")
    @test_throws IncompleteData decodeint(b"b\0\0")
    @test_throws IncompleteData decodeint(b"b\0\0\0")
    @test decodeint(b"b\0\0\0\0") == (0, [])
    @test decodeint(b"b\0\0\0\0tail") == (0, b"tail")
    @test decodeint(b"b\x7f\xff\xff\xff") == (2147483647, [])
    @test decodeint(b"b\x7f\xff\xff\xfftail") == (2147483647, b"tail")
    @test decodeint(b"b\xff\xff\xff\xff") == (-1, [])
    @test decodeint(b"b\xff\xff\xff\xfftail") == (-1, b"tail")
end

# decode binary
testcase() do
    @test_throws IncompleteData decode(b"\x83m")
    @test_throws IncompleteData decode(b"\x83m\0")
    @test_throws IncompleteData decode(b"\x83m\0\0")
    @test_throws IncompleteData decode(b"\x83m\0\0\0")
    @test_throws IncompleteData decode(b"\x83m\0\0\0\1")
    @test decode(b"\x83m\0\0\0\0") == (b"", b"")
    @test decode(b"\x83m\0\0\0\0tail") == (b"", b"tail")
    @test decode(b"\x83m\0\0\0\4data") == (b"data", b"")
    @test decode(b"\x83m\0\0\0\4datatail") == (b"data", b"tail")
    @test decodebin(b"m\0\0\0\0") == (b"", b"")
    @test decodebin(b"m\0\0\0\0tail") == (b"", b"tail")
    @test decodebin(b"m\0\0\0\4data") == (b"data", b"")
    @test decodebin(b"m\0\0\0\4datatail") == (b"data", b"tail")
    @test_throws IncompleteData decodebin(b"m")
    @test_throws IncompleteData decodebin(b"m\0")
    @test_throws IncompleteData decodebin(b"m\0\0")
    @test_throws IncompleteData decodebin(b"m\0\0\0")
    @test_throws IncompleteData decodebin(b"m\0\0\0\1")
end

# decode float
testcase() do
    @test_throws IncompleteData decode(b"\x83F")
    @test_throws IncompleteData decode(b"\x83F\0")
    @test_throws IncompleteData decode(b"\x83F\0\0")
    @test_throws IncompleteData decode(b"\x83F\0\0\0")
    @test_throws IncompleteData decode(b"\x83F\0\0\0\0")
    @test_throws IncompleteData decode(b"\x83F\0\0\0\0\0")
    @test_throws IncompleteData decode(b"\x83F\0\0\0\0\0\0")
    @test_throws IncompleteData decode(b"\x83F\0\0\0\0\0\0\0")
    @test decode(b"\x83F\0\0\0\0\0\0\0\0") == (0.0, b"")
    @test decode(b"\x83F\0\0\0\0\0\0\0\0tail") == (0.0, b"tail")
    @test decode(b"\x83F?\xf8\0\0\0\0\0\0") == (1.5, b"")
    @test decode(b"\x83F?\xf8\0\0\0\0\0\0tail") == (1.5, b"tail")
    @test decodefloat(b"F\0\0\0\0\0\0\0\0") == (0.0, b"")
    @test decodefloat(b"F\0\0\0\0\0\0\0\0tail") == (0.0, b"tail")
    @test decodefloat(b"F?\xf8\0\0\0\0\0\0") == (1.5, b"")
    @test decodefloat(b"F?\xf8\0\0\0\0\0\0tail") == (1.5, b"tail")
    @test_throws IncompleteData decodefloat(b"F")
    @test_throws IncompleteData decodefloat(b"F\0")
    @test_throws IncompleteData decodefloat(b"F\0\0")
    @test_throws IncompleteData decodefloat(b"F\0\0\0")
    @test_throws IncompleteData decodefloat(b"F\0\0\0\0")
    @test_throws IncompleteData decodefloat(b"F\0\0\0\0\0")
    @test_throws IncompleteData decodefloat(b"F\0\0\0\0\0\0")
    @test_throws IncompleteData decodefloat(b"F\0\0\0\0\0\0\0")
end

# decode small big integer

# decode big integer

# decode compressed term
testcase() do
    @test_throws IncompleteData decode(b"\x83P")
    @test_throws IncompleteData decode(b"\x83P\0")
    @test_throws IncompleteData decode(b"\x83P\0\0")
    @test_throws IncompleteData decode(b"\x83P\0\0\0")
    @test_throws IncompleteData decode(b"\x83P\0\0\0\0")
    @test_throws InvalidCompressedTag decode(badsizecompdata)
    @test decode(compdata1) == ([100,100,100,100,100,100,100,100,100,100,
                                 100,100,100,100, 100,100,100,100,100,100],
                                Uint8[])
    # XXX the following test fails because the Zlib library for Julia doesn't
    # provide a flush-like mechanism that the Python zlib library does
    # @test decode(compdata2) == ([100,100,100,100,100,100,100,100,100,100,
    #                              100,100,100,100, 100,100,100,100,100,100],
    #                             b"tail")
end


text:
using QuasiGeostrophy

@testset "Domain Tests" begin
    Ω = S¹(0,2π)
    @test dim(Ω)==1
    @test dim(Ω×Ω)==2
    @test dim(Ω×Ω×Ω)==3
end

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
using LTC
using Plots
gr()
using DiffEqFlux
using BenchmarkTools
using DiffEqSensitivity
using OrdinaryDiffEq
using GalacticOptim
using Juno
using Cthulhu
using Profile
using BlackBoxOptim
using ModelingToolkit
using ProfileView

function generate_data()
    in_features = 2
    out_features = 1
    N = 48
    data_x = [sin.(range(0,stop=3π,length=N)), cos.(range(0,stop=3π,length=N))]
    data_x = [reshape([Float32(data_x[1][i]),Float32(data_x[2][i])],2,1) for i in 1:N]# |> f32
    data_y = [reshape([Float32(y)],1,1) for y in sin.(range(0,stop=6π,length=N))]# |> f32
    DataLoader((data_x, data_y), batchsize=N)
end

function data(iter; data_x=nothing, data_y=nothing, short=false, noisy=false)
    ncycle(generate_data(), iter)
end

function loss(p,x::Vector{Matrix{T}},y::Vector{Matrix{T}}, m) where T

  # ŷ = bcm(m,x,p)
  # ŷ = m.(x, [p])

  ŷb = GalacticOptim.Flux.Zygote.Buffer(y[1], size(y,1), size(y[1])[1],size(y[1])[2])
  for (i, xi) in enumerate(x)
    ŷi = m(xi, p)#::Matrix{T}
    ŷb[i,:,:] = ŷi
    Inf32 ∈ ŷi && return Inf32, Flux.unstack(copy(ŷb)::Array{T,3},1)::Vector{Matrix{T}}, y # TODO: what if a layer after MTKRecur can handle Infs?
  end

  # ŷ = copy(ŷb)
  # ŷ = [ŷ[i,:,:] for i in 1:size(ŷ)[1]]
  ŷ = Flux.unstack(copy(ŷb)::Array{T,3},1)::Vector{Matrix{T}}

  LTC.reset_state!(m, p)
  # mean(bcl(GalacticOptim.Flux.Losses.mse, ŷ, y)), ŷ, y
  # GalacticOptim.Flux.Losses.mse(ŷ, y, agg=mean)[1], ŷ, y
  mean(sum.(abs2, (ŷ .- y))), ŷ, y
  # bcl(GalacticOptim.Flux.Losses.mse, ŷ, y), ŷ, y
end

function cbg(p,l,pred,y;doplot=false)
  display(l)
  if doplot
    fig = plot([ŷ[end,1] for ŷ in pred], label="ŷ")
    plot!(fig, [yi[end,1] for yi in y], label="y")
    #frame(anim)
    display(fig)
  end
  return false
end

function test()
  n = 5
  solver=Tsit5()
  sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true))
  batchsize = 1

  wiring = LTC.FWiring(2,1)
  net = LTC.Net(wiring; name=:net)
  recurmtk = LTC.RecurMTK(LTC.MTKCell(wiring.n_in, wiring.n_out, net, solver, sensealg))

  model = DiffEqFlux.FastChain(LTC.Mapper(wiring.n_in),
                               recurmtk,
                               LTC.Mapper(wiring.n_out),
                               )

  train_dl = data(n)
  opt = GalacticOptim.Flux.Optimiser(ClipValue(0.5), ADAM(0.01))
  # opt = BBO()
  AD = GalacticOptim.AutoZygote()
  # AD = GalacticOptim.AutoReverseDiff()

  x1,y1 = first(train_dl)

  p_mtk = initial_params(recurmtk)
  lb_mtk, ub_mtk = get_bounds(recurmtk)

  # @code_warntype loss(p_mtk,x1,y1,recurmtk)
  # @descend_code_warntype loss(p_mtk,x1,y1,recurmtk)
  # Juno.@profiler loss(p_mtk,x1,y1,recurmtk)
  # loss(p_mtk,x1,y1,recurmtk)

  optfun = GalacticOptim.OptimizationFunction((θ, p, x, y) -> loss(θ,x,y, recurmtk), AD)
  optfunc = GalacticOptim.instantiate_function(optfun, p_mtk, AD, nothing)
  optprob = GalacticOptim.OptimizationProblem(optfunc, p_mtk, lb=lb_mtk, ub=ub_mtk,
                                #grad = true, hess = true, sparse = true,
                                #parallel=ModelingToolkit.MultithreadedForm()
                                )

  # optfun = GalacticOptim.OptimizationFunction((θ, p, x, y) -> loss(θ,x,y, recurmtk),AutoModelingToolkit(),p_mtk,DiffEqBase.NullParameters(),
  #                    #grad = true, hess = true, sparse = true,
  #                    #checkbounds = false,
  #                    #parallel=ModelingToolkit.MultithreadedForm(),
  #                    )
  # optfunc = GalacticOptim.instantiate_function(optfun, p_mtk, AutoModelingToolkit(), DiffEqBase.NullParameters())
  # optprob = GalacticOptim.OptimizationProblem(optfunc, p_mtk, lb=lb_mtk, ub=ub_mtk,
  #                               #grad = true, hess = true, sparse = true,
  #                               #parallel=ModelingToolkit.MultithreadedForm()
  #                               )

  # GalacticOptim.solve(optprob, opt, train_dl, cb = cbg)

  res = Vector{}(undef,length(p_mtk))
  @time optfunc.grad(res,p_mtk,x1,y1)
  @time optfunc.grad(res,p_mtk,x1,y1)
  @time optfunc.grad(res,p_mtk,x1,y1)
  @time optfunc.grad(res,p_mtk,x1,y1)
end
test()

n = 5
solver=Tsit5()
sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true))
batchsize = 1

wiring = LTC.FWiring(2,1)
net = LTC.Net(wiring; name=:net)
# sys = ModelingToolkit.structural_simplify(net)
recurmtk = LTC.RecurMTK(LTC.MTKCell(wiring.n_in, wiring.n_out, net, solver, sensealg))

model = DiffEqFlux.FastChain(LTC.Mapper(wiring.n_in),
                             recurmtk,
                             LTC.Mapper(wiring.n_out),
                             )

train_dl = data(n)
opt = GalacticOptim.Flux.Optimiser(ClipValue(0.5), ADAM(0.01))
# opt = BBO()
AD = GalacticOptim.AutoZygote()
# AD = GalacticOptim.AutoReverseDiff()

x1,y1 = first(train_dl)

p_mtk = initial_params(recurmtk)
lb_mtk, ub_mtk = get_bounds(recurmtk)
h_mtk = recurmtk.state

# h_mtkcell, ŷ_mtkcell = recurmtk.cell(h_mtk, x1[1], p_mtk)
# ŷ_mtkcell = recurmtk(x1[1], p_mtk)
# LTC.reset_state!(recurmtk, p_mtk)

@code_warntype recurmtk.(x1, [p_mtk])

# brc(m, x, p) = m.(x, [p])
# @code_warntype brc(recurmtk, x1, p_mtk)

@code_warntype loss(p_mtk,x1,y1,recurmtk)
@descend_code_warntype loss(p_mtk,x1,y1,recurmtk)

ŷ_mtk = recurmtk.(x1, [p_mtk])
@code_warntype mean(GalacticOptim.Flux.Losses.mse.(ŷ_mtk, y1))
@code_warntype GalacticOptim.Flux.Losses.mse(ŷ_mtk, y1)[1]
@code_warntype LTC.reset_state!(recurmtk, p_mtk)

@code_warntype LTC.solve_ensemble(recurmtk.cell, h_mtk, x1[1], p_mtk)
@descend_code_warntype LTC.solve_ensemble(recurmtk.cell, h_mtk, x1[1], p_mtk)

@descend_code_warntype LTC.optimize(recurmtk, loss, cbg, opt, AD, train_dl)

optfun = GalacticOptim.OptimizationFunction((θ, p, x, y) -> loss(θ,x,y, recurmtk), AD)
optfunc = GalacticOptim.instantiate_function(optfun, p_mtk, AD, nothing)
optprob = GalacticOptim.OptimizationProblem(optfunc, p_mtk, lb=lb_mtk, ub=ub_mtk,
                              #grad = true, hess = true, sparse = true,
                              #parallel=ModelingToolkit.MultithreadedForm()
                              )
# GalacticOptim.solve(optprob, opt, train_dl, cb = cbg)

res = Vector{}(undef,length(p_mtk))
@time optfunc.grad(res,p_mtk,x1,y1)
@code_warntype LTC.optimize(recurmtk, loss, cbg, opt, AD, train_dl)

ps = Flux.params(pp)
gs = Flux.Zygote.gradient(ps) do
    x = prob.f(θ,prob.p, d...)
    first(x)
  end

# p = initial_params(model)
# lb, ub = get_bounds(model)
#
# ŷ = model.(x1, [p])
# loss(p, x1, y1, model)
#
# LTC.optimize(model, loss, cbg, opt, AD, train_dl)


text:
const JAVA_HOME_CANDIDATES = ["/usr/lib/jvm/default-java/",
                              "/usr/lib/jvm/default/"]

struct JavaCallError <: Exception
    msg::String
end

function javahome_winreg()
    keys = ["SOFTWARE\\JavaSoft\\Java Runtime Environment", "SOFTWARE\\JavaSoft\\Java Development Kit", "SOFTWARE\\JavaSoft\\JRE", "SOFTWARE\\JavaSoft\\JDK"]

    for key in keys
        try
            value = querykey(WinReg.HKEY_LOCAL_MACHINE, key, "CurrentVersion")
            key *= "\\" * value
            return querykey(WinReg.HKEY_LOCAL_MACHINE, key, "JavaHome")
        catch
            # Try the next value in the loop.
        end
    end

    error("Cannot find an installation of Java in the Windows Registry. Please install a JRE/JDK, or set the JAVA_HOME environment variable if one is already installed.")
end

@static Sys.isunix() ? (global const libname = "libjvm") : (global const libname = "jvm")

function findjvm()
    javahomes = Any[]
    libpaths = Any[]

    if haskey(ENV, "JAVA_HOME")
        push!(javahomes, ENV["JAVA_HOME"])
    else
        @static if Sys.iswindows()
            ENV["JAVA_HOME"] = javahome_winreg()
            push!(javahomes, ENV["JAVA_HOME"])
        end
        @static if Sys.isunix()
            # Find default javahome by checking location of the java command
            try
                javapath = chomp(read(`which java`,String))
                while(islink(javapath))
                    javapath = readlink(javapath)
                end
                javapath = dirname(javapath)
                javapath = match(r"(.*)((/jre/bin)|(/bin))+",javapath)[1]
                push!(javahomes,javapath)
            catch err
                @debug "JavaCall could not determine javapath from `which java`" err
            end
        end
    end
    isfile("/usr/libexec/java_home") && push!(javahomes, chomp(read(`/usr/libexec/java_home`, String)))

    for fname ∈ JAVA_HOME_CANDIDATES
        isdir(fname) && push!(javahomes, fname)
    end

    push!(libpaths, pwd())
    for n in javahomes
        @static if Sys.iswindows()
            push!(libpaths, joinpath(n, "bin", "server"))
            push!(libpaths, joinpath(n, "jre", "bin", "server"))
            push!(libpaths, joinpath(n, "bin", "client"))
        end
        @static if Sys.islinux()
            if Sys.WORD_SIZE==64
                push!(libpaths, joinpath(n, "jre", "lib", "amd64", "server"))
                push!(libpaths, joinpath(n, "lib", "amd64", "server"))
            elseif Sys.WORD_SIZE==32
                push!(libpaths, joinpath(n, "jre", "lib", "i386", "server"))

                push!(libpaths, joinpath(n, "lib", "i386", "server"))
            end
        end
        push!(libpaths, joinpath(n, "jre", "lib", "server"))
        push!(libpaths, joinpath(n, "lib", "server"))
    end

    ext = @static Sys.iswindows() ? "dll" : (@static Sys.isapple() ? "dylib" : "so")
    ext = "."*ext

    try
        for n in libpaths
            libpath = joinpath(n,libname*ext);
            if isfile(libpath)
                if Sys.iswindows()
                    bindir = dirname(dirname(libpath))
                    m = filter(x -> occursin(r"msvcr(?:.*).dll",x), readdir(bindir))
                    if !isempty(m)
                        return (joinpath(bindir,m[1]),libpath)
                    end
                end
                return (libpath,)
            end
        end
    catch err
        throw(err)
    end

    errorMsg =
        [
         "Cannot find java library $libname$ext\n",
         "Search Path:"
         ];
    for path in libpaths
        push!(errorMsg,"\n   $path")
    end
    throw(JavaCallError(reduce(*,errorMsg)))
end

@static Sys.isunix() ? (const sep = ":") : nothing
@static Sys.iswindows() ? (const sep = ";") : nothing
cp = OrderedSet{String}()
opts = OrderedSet{String}()

"""
    JavaCall.getClassPath()

    Obtains the Java classpath.

    Before JavaCall.init(), this reports the classpath constructed using
    JavaCall.addClasspath.

    After JavaCall.init(), this reports System.getProperty("java.class.path")
"""
function getClassPath()
    if isloaded()
        # jls = @jimport java.lang.System
        jls = JavaObject{Symbol("java.lang.System")}
        return jcall(jls, "getProperty", JString, (JString,), "java.class.path")::String
    else
        ccp = collect(cp)
        return join(ccp,sep)
    end
end

"""
    JavaCall.addClassPath(s::String)

    Add a string to the classpath. Must be called before JavaCall.init

    foo/*.jar will add all the jars in the directory foo to the classpath
    foo/*     will add all the jars and directories recursively including foo

    See also addClassPathRecursive and addJarsToClassPath
"""
function addClassPath(s::String)
    if isloaded()
        @warn("JVM already initialised. This call has no effect")
        return
    end
    if s==""; return; end
    dirname, pattern = splitdir(s)
    if pattern == "*.jar" && isdir(dirname)
        _addJarsToClassPath(dirname)
    elseif pattern == "*" && isdir(dirname)
        _addClassPathRecursive(dirname)
    else
        push!(cp, s)
    end
    return
end

function _addClassPathRecursive(dirname::String)
    push!(cp,dirname)
    for (root,dirs,filenames) in walkdir(dirname)
        addJarsToClassPath(filenames,root)
        union!(cp, [joinpath(root,dir) for dir in dirs])
    end
end
"""
    JavaCall.addClassPathRecursive(dirname::String)

    Adds dirname and all jars and directories in dirname recursively to the classpath
"""
addClassPathRecursive(dirname) = isloaded() ?
    @warn("JVM already initialized. This call has no effect") :
    _addClassPathRecursive(dirname)

function _addJarsToClassPath(files::Array{String,1},path::String="")
    for filename in files
        if endswith(filename,".jar")
            push!(cp, joinpath(path,filename) )
        end
    end
end
_addJarsToClassPath(dirname::String) = addJarsToClassPath(readdir(dirname),dirname)
"""
    JavaCall.addJarsToClassPath(dirname::String)

    Add jars in dirname to the classpath.
    Equivalent to JavaCall.addJarsToClasspath(readdir(dirname),dirname)

    JavaCall.addJarsToClassPath(files::Array{String,1}, [path::String])

    Add files which end in ".jar" to the classpath prefixed by an optional path
"""
addJarsToClassPath(args...) = isloaded() ?
    @warn("JVM already initialized. This call has no effect") :
    _addJarsToClassPath(args...)

function addOpts(s::String)
    if isloaded()
        @warn("JVM already initialised. This call has no effect")
    else
        m = match(r"^-Djava.class.path=(.*)",s)
        if m != nothing
            addClassPath(String(m.captures[1]))
        else
            push!(opts, s)
        end
    end
end

const ROOT_TASK_ERROR = JavaCallError(
    "Either the environmental variable JULIA_COPY_STACKS must be 1 " *
    "OR JavaCall must be used on the root Task.")

const JULIA_COPY_STACKS_ON_WINDOWS_ERROR = JavaCallError(
    "JULIA_COPY_STACKS should not be set on Windows.")

const THREADID_NOT_ONE_WINDOWS_ERROR = JavaCallError(
    "JavaCall must be used on Thread 1 only in Windows. Multithreading JavaCall is not supported on Windows."
)

# JavaCall must run on the root Task or JULIA_COPY_STACKS is enabled
isroottask() = Base.roottask === Base.current_task()
@static if Sys.iswindows()
    isgoodenv() = ( ! JULIA_COPY_STACKS ) && Base.Threads.threadid() == 1
    assertroottask_or_goodenv() = isgoodenv() ? true : Base.Threads.threadid() == 1 ?
        throw(JULIA_COPY_STACKS_ON_WINDOWS_ERROR) : throw(THREADID_NOT_ONE_WINDOWS_ERROR)
else
    isgoodenv() = JULIA_COPY_STACKS || isroottask()
    assertroottask_or_goodenv() = isgoodenv() ? true : throw(ROOT_TASK_ERROR)
end

isloaded() = JNI.is_jni_loaded() && JNI.is_env_loaded()

assertloaded() = isloaded() ? true : throw(JavaCallError("JVM not initialised. Please run init()"))
assertnotloaded() = isloaded() ? throw(JavaCallError("JVM already initialised")) : true

"""
    JavaCall.init(opts::Array{String,1})
    JavaCall.init(opt1::String,opt2::String, ...)
    JavaCall.init()

    Initialize JavaCall with JVM options.

    As of JavaCall v0.7.4 new options passed to init will be appended
    to previous options added with addClasspath and addOpts.

    Once init() is called, addClasspath and addOpts no longer have any effect.

    See http://juliainterop.github.io/JavaCall.jl/methods.html

    Example
    JavaCall.init(["-Xmx512M", "-Djava.class.path=$(@__DIR__)", "-verbose:jni", "-verbose:gc"])
"""
function init(opts::Array{String,1})
    addOpts.(opts)
    init()
end

# Accept options strings as a set of arguments
init(opts::Array{AbstractString,1}) = init(String.(opts))
init(opts::Vararg{AbstractString,N}) where N = init(String.([opts...]))

function init()
    if isempty(cp)
        _init(opts)
    else
        ccp = collect(cp)
        options = collect(opts)
        classpath = "-Djava.class.path="*join(ccp,sep)
        _init(vcat(options, classpath))
    end
end

# Below is the original main initialization option
# Pointer to pointer to pointer to pointer alert! Hurrah for unsafe load
function _init(opts)
    assertnotloaded()
    assertroottask_or_goodenv()
    JNI.init_new_vm(findjvm(),opts);
end

"""
    init_current_vm()

Allow initialization from running VM. Uses the first VM it finds.

# Example using JNA

```java
package zot.julia;

import com.sun.jna.Native;

public class Julia {
    static {
        Native.register("julia");
        jl_init__threading();
    }

    public static double bubba = Math.random();

    public static native void jl_init__threading();
    public static native void jl_eval_string(String code);
    public static native void jl_atexit_hook(int status);

    public static void main(String args[]) {
        System.out.println("test");
        jl_eval_string("println(\"test from Julia\")");
        jl_eval_string("using JavaCall");
        jl_eval_string("JavaCall.init_current_vm()");
        jl_eval_string("println(\"initialized VM\")");
        jl_eval_string("jlm = @jimport java.lang.Math");
        jl_eval_string("println(jcall(jlm, \"sin\", jdouble, (jdouble,), pi/2))");
        jl_eval_string("jl = @jimport zot.julia.Julia");
        System.out.println("Bubba should be " + bubba);
        jl_eval_string("println(\"bubba: \", jfield(jl, \"bubba\", jdouble))");
        jl_eval_string("println(\"Done with tests\")");
        jl_atexit_hook(0);
    }
}
```
"""
function init_current_vm()
    JNI.init_current_vm(findjvm())
end

function destroy()
    assertroottask_or_goodenv()
    JNI.destroy()
end


text:
module GRAYTR

include("geometry.jl")
using .Geometry

using LinearAlgebra
using Random
using Distributed
using FileIO
using Images
using StaticArrays

include("abstracts.jl")
include("utility.jl")
include("rays.jl")
include("diffgeom.jl")
include("shape.jl")
include("samplers.jl")
include("bounding.jl")

include("spectrum.jl")
include("brdf.jl")

include("lights.jl")
include("primitives.jl")
include("dumb_aggregate.jl")
include("accelerators.jl")
include("scene.jl")
include("integrator.jl")
include("renderer.jl")
include("filters.jl")

include("cameras/projective.jl")

include("films/imagefilm.jl")
include("films/delayfilm.jl")

export Vector3, Point3, Normal3
export X_AXIS, Y_AXIS, Z_AXIS, NEG_X_AXIS, NEG_Y_AXIS, NEG_Z_AXIS, IDENTITY_TRANSFORM
export Transformation, rotation, translation, scaling

export ImageFilm, DelayFilm, PhotometricFilm
export TriangleFilter, BoxFilter
export OrthographicCamera
export StratifiedSampler, WhittedIntegrator, PressureIntegrator
export render
export write_image, write_txt, reset!

end #module

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
# This is Subset Tournament CMSA-ES as proposed by Robert Feldt in the paper:
#  R. Feldt, "Covariate Subset Tournaments for High-Dimensional Blackbox Optimization with Covariance Matrix Adapting Evolutionary Strategies", 2014

function normalize_utilities(utilities)
  utilities / sum(utilities)
end

function linear_utilities(mu, lambda)
  normalize_utilities(hcat(ones(1, mu), zeros(1, lambda-mu)))
end

function log_utilities(mu, lambda)
  normalize_utilities(hcat((log(mu+1) - log(1:mu))', zeros(1, lambda-mu)))
end

# Optimize the n-dimensional objective func with a (mu,lambda) CMSA-ES using
# subset tournaments to optimize subsets of variables in rounds.
function st_cmsa_es(p;
  trace = true,

  # Stopping criteria related
  max_seconds = 4*numdims(p), max_evals_per_dim = 1e7,
  ftol = 1e-7, xtol = 1e-10, stol = 1e-10,
  max_rounds_without_improvement = 200,

  # Starting points, will be random unless specified
  xmean = false,

  # Algorithm specific params:
  covarMatrixSampler = SubsetCholeskyCovarSampler,
  utilitiesFunc = log_utilities,
  lambda = 4*numdims(p),
  mu = int(max(ceil(lambda/rand(4:20)), 1.0)),
  tau = 1 / sqrt(2*numdims(p)), # Equation (1) on page 5 in Beyer2008
  tau_c = 1 + numdims(p) * (numdims(p) + 1) / (2 * mu),  # Equation (2) on page 5 in Beyer2008
  sigma = 0.05*rand(1:8)*minimum(diameters(search_space(p))),
  decompose_covar_prob = 0.4,

  # Subsetting specific parameters:
  subset_size = int(floor(0.20*numdims(p))),
  subsets_per_tournament = 2,
  num_rounds_per_tournament = 1,
  num_rounds_of_optimization_between_tournaments = 30,
  subset_selection_mechanism = :random

  )

  N = numdims(p)
  ss = search_space(p)
  max_evals = max_evals_per_dim * N
  a = 1 - 1 / tau_c

  C = covarMatrixSampler(N)
  utilities = utilitiesFunc(mu, lambda)
  xbest = xmean = rand_individual(ss)   # Current best mean value.
  fbest = eval1(xbest, p)
  num_fevals = 1
  archive = BlackBoxOptim.TopListArchive(N, 10)
  add_candidate!(archive, fbest, xbest[:], num_fevals)
  fevals_last_best = num_fevals
  next_print_covar = 100
  termination_reason = "?"

  # Init for subset tournaments
  st_state = :optimization   # Can be either :tournament or :optimization and starts with :optimization to ensure we set up a new tournament
  num_optimization_rounds_for_this_subset = num_rounds_of_optimization_between_tournaments # Also ensures we are end of opt round => new tournament will be set up
  st_current_subset = 1  # When in :tournament mode this is the index to the subset currently being evaluated, can be in range 1-subsets_per_tournament
  st_subsets = Array{Int, 1}[] # Subsets that are currently in a tournament
  num_tournament_rounds_for_this_subset = 0   # Number of eval rounds we have ran with current subset in tournament mode
  fitness_per_tournament_round = zeros(num_rounds_per_tournament, subsets_per_tournament)

  # Keep stats per covariate being optimized for how effective the optimization is when they are included.
  # We save the expected relative change per round of optimization and use a history parameter in [0.0, 1.0]
  # to decide how much history should be saved when updating it, i.e. newval = oldvalue * h + delta * (1-h).
  # These stats values are used in multiple tournaments like so:
  #   1. Repeat until the right number of vars has been added to selected set
  #   2.   Randomly sample two unselected vars
  #   3.   Select the one with the best stats value (lowest value if minimizing)
  improvement_per_var = 0.01 * ones(N) # Start with some value since it will be adapted as we go...

  start_time = time()

  # Now lets optimize! Ensure we run at least one iteration.
  while(true)

    if (time() - start_time) > max_seconds
      termination_reason = "Exceeded time budget"
      break
    end

    if num_fevals > max_evals
      termination_reason = "Exceeded function eval budget"
      break
    end

    if st_state == :tournament
      if num_tournament_rounds_for_this_subset < num_rounds_per_tournament
        num_tournament_rounds_for_this_subset += 1
      else
        st_current_subset += 1
        if st_current_subset <= subsets_per_tournament
          set_subset!(C, st_subsets[st_current_subset])
          num_tournament_rounds_for_this_subset = 1
        else # No more tournaments needed, select best subset and start optimizing with it
          winning_subset = select_winning_subset(st_subsets, fitness_per_tournament_round)
          set_subset!(C, st_subsets[winning_subset])
          st_current_subset = winning_subset
          st_state = :optimization
          num_optimization_rounds_for_this_subset = 1
        end
      end
    else # In optimization mode
      if num_optimization_rounds_for_this_subset < num_rounds_of_optimization_between_tournaments
        num_optimization_rounds_for_this_subset += 1
      else
        # We have now used the current subset for a full set of optimization rounds so we set up for a new tournament
        st_subsets = select_new_subsets(st_subsets, N, subsets_per_tournament, st_current_subset, subset_selection_mechanism, subset_size)
        st_state = :tournament
        st_current_subset = 1
        set_subset!(C, st_subsets[st_current_subset])
        num_tournament_rounds_for_this_subset = 1
      end
    end

    # Decompose only with a given probability => saves time
    if rand() <= decompose_covar_prob
      decompose!(C)
    end

    # Generate new population
    sigmas = sigma * exp( tau * randn(1, lambda) )  # 1*lambda
    s = multivariate_normal_sample(C, N, lambda)
    z = broadcast(*, sigmas, s)                     # n*lambda
    xs = repmat(xmean, 1, lambda) + z               # n*lambda

    # Evaluate fitness
    fitnesses = eval_fitnesses(p, xs, lambda)
    num_fevals += lambda

    # Check if best new fitness is best ever and print some info if tracing.
    indbest = argmin(fitnesses)
    fbest_new = fitnesses[indbest]

    # Save info about the fitnesses if we are in tournament mode
    if st_state == :tournament
      fitness_per_tournament_round[num_tournament_rounds_for_this_subset, st_current_subset] = fbest_new
    end

    if fbest_new < fbest
      xbest = xs[:, indbest]
      fbest = fbest_new
      add_candidate!(archive, fbest, xbest, num_fevals)
      fevals_last_best = num_fevals

      if fitness_is_within_ftol(p, ftol, fbest)
        termination_reason = "Within ftol"
        break
      end

      if trace
        println("$(num_fevals): Best fitness = $(fbest)")
        if num_fevals > next_print_covar
          next_print_covar = num_fevals + 100
          #println("covar summary: ", sumstats(C.C, (x) -> @sprintf("%.2e", x)))
          println("sigma: ", sigma)
        end
      end
    else
      if (num_fevals - fevals_last_best) > max_rounds_without_improvement * lambda
        termination_reason = "Max rounds without improvement reached"
        break
      end
    end

    # Assign weights to the best individuals according to the utilities vector.
    weights = assign_weights(lambda, fitnesses, utilities)

    # Calculate new mean value based on weighted version of steps in population.
    xmean += (z * weights)

    # Update the covariance matrix
    uc = zeros(Float64, N, N)
    for i in 1:lambda
      if weights[i] > 0.0
        se = s[:,i]
        uc += weights[i] * (se * se')
      end
    end
    update_covariance_matrix!(C, uc, a)
    #ws = broadcast(*, weights', s)
    #update_covariance_matrix!(C, (ws * ws'), covar_learning_rate)

    # Adapt sigma for next round
    sigma = sigmas * weights

    # Terminate if sigma very small
    #println("sigma = $(sigma[1,1])")
    if sigma[1,1] < stol
      termination_reason = "Sigma too small"
      break
    end
  end

  return xbest, fbest, num_fevals, termination_reason, archive
end

function select_winning_subset(subsets, fitness_per_tournament_round)
  if length(size(fitness_per_tournament_round)) > 1
    fitness_summary = mean(fitness_per_tournament_round, 1)
    argmin(fitness_summary)
  else
    argmin(fitness_per_tournament_round)
  end
end

function select_new_subsets(subsets, n, num_subsets, winning_subset, subset_selection_mechanism, subset_size)
  if length(subsets) == 0 || subset_selection_mechanism == :random

    generate_random_subsets(num_subsets, subset_size, n)

  else # We should keep the winner

    # Generate new subsets randomly
    new_subsets = generate_random_subsets(num_subsets-1, subset_size, n)

    # Add the previous winning subset
    push!(new_subsets, subsets[winning_subset])
    new_subsets

  end
end

# Select subsets based on running binary tournaments between stats values
function select_new_subsets_based_on_stats(n, num_subsets, subset_size, stats, smaller_stats_is_better)
  [generate_subset_based_on_stats(subset_size, n, stats, smaller_stats_is_better) for i in 1:num_subsets]
end

function generate_subset_based_on_stats(subset_size, n, stats, smaller_stats_is_better = false)
  subset = Int[]
  candidates = shuffle(collect(1:n))
  len = n-1
  op = smaller_stats_is_better ? < : >

  for i in 1:subset_size
    index = rand(1:len)
    candidate1 = candidates[index]
    candidate2 = candidates[index+1]
    println("Comparing $(candidate1) and $(candidate2)")
    if op( stats[candidate1], stats[candidate2] )
      splice!(candidates, index)
      push!(subset, candidate1)
    else
      splice!(candidates, index+1)
      push!(subset, candidate2)
    end
    len -= 1
  end

  subset
end

function generate_random_subsets(num_subsets, subset_size, n)
  new_subsets = Array{Int, 1}[]
  for i in 1:num_subsets
    push!(new_subsets, sort(shuffle(collect(1:n))[1:subset_size]))
  end
  new_subsets
end

function eval_fitnesses(problem, xs, lambda = size(xs, 2))
  fitnesses = zeros(lambda)
  for i in 1:lambda
    fitnesses[i] = eval1(xs[:,i], problem)
  end
  fitnesses
end

function assign_weights(n, fitnesses, utilities; minimize = true)
  us_ordered = zeros(n, 1)
  perms = sortperm(fitnesses, rev = !minimize)
  for i in 1:n
    us_ordered[perms[i]] = utilities[i]
  end
  us_ordered
end

# We create different types of Covariance matrix samplers based on different
# decompositions.
abstract type CovarianceMatrixSampler end

function update_covariance_matrix!(cms::CovarianceMatrixSampler, delta, a)
  C = a * cms.C + (1 - a) * delta
  cms.C = triu(C) + triu(C,1)' # Ensure C is symmetric. Should not be needed, investigate...
end

mutable struct EigenCovarSampler <: CovarianceMatrixSampler
  C::Array{Float64,2}
  B::Array{Float64,2}
  diagD::Array{Float64,1}

  EigenCovarSampler(n) = begin
    new(Matrix{Float64}(I, n,n), Matrix{Float64}(I, n,n), ones(n))
  end
end

function decompose!(cms::EigenCovarSampler)
  try
    EV, B = eig(cms.C)
    cms.B = B
    cms.diagD = sqrt(EV)
  catch
    # We don't update if there is some problem
  end
end

function multivariate_normal_sample(cms::CovarianceMatrixSampler, n, m)
  cms.B * (cms.diagD .* randn(n, m))
end

mutable struct CholeskyCovarSampler <: CovarianceMatrixSampler
  C::Array{Float64,2}
  sqrtC::Array{Float64,2}

  CholeskyCovarSampler(n) = begin
    new(Matrix{Float64}(I, n,n), Matrix{Float64}(I, n,n))
  end
end

function decompose!(cms::CholeskyCovarSampler)
  try
    cms.sqrtC = chol(cms.C)'
  catch error
    # We don't update if there is some problem
  end
end

function multivariate_normal_sample(cms::CholeskyCovarSampler, n, m)
  cms.sqrtC * randn(n, m)
end

# The subset CholeskyCovarSampler only uses a subset of the variables in the
# (expensive) cholesky decomposition, the other variables are kept constant.
# However, the covariance matrix itself is always updated and saved in full
# so that the overall learning of the shape of the fitness landscape is not lost.
mutable struct SubsetCholeskyCovarSampler <: CovarianceMatrixSampler
  C::Array{Float64,2}
  sqrtC::Array{Float64,2}
  subset::Array{Int, 1} # Indices of the currently active subset of variables

  SubsetCholeskyCovarSampler(n) = begin
    new(Matrix{Float64}(I, n,n), Matrix{Float64}(I, n,n), collect(1:n))
  end
end

# Set new subset.
function set_subset!(cms::SubsetCholeskyCovarSampler, subset)
  cms.subset = subset
  decompose!(cms)
end

using JSON

function decompose!(cms::SubsetCholeskyCovarSampler)
  try
    subset_C = cms.C[cms.subset, cms.subset]
    cms.sqrtC = chol(subset_C)'
  catch error
    # We don't update if there is some problem
    println("ERROR: Could not do cholesky decomposition!")
    show(error)
  end
end

function multivariate_normal_sample(cms::SubsetCholeskyCovarSampler, n, m)
  subset_size = length(cms.subset)

  subset_inv_c = 1

  # Calc the inverted covar matrix only for the subset
  try
    subset_inv_c = cms.sqrtC * randn(subset_size, m)
  catch error
    println("Error in subset_inv_c multiplication, size(cms.sqrtC) = $(size(cms.sqrtC)), subset_size = $(subset_size), m = $(m), size(cms.C) = $(size(cms.C))")
    show(error)
  end

  # The rest should be zero. Maybe faster if we do this with a sparse matrix
  # since most of them will be zero?
  inv_c = zeros(n, m)
  for i in 1:subset_size
    si = cms.subset[i]
    for j in 1:m
      inv_c[si, j] = subset_inv_c[i, j]
    end
  end

  inv_c
end