text:
using Random
using Distances
using Plots

include("..\\elements\\Elements.jl")
include("..\\file_input\\FileLoaders.jl")
include("..\\configurations\\MasterConfiguration.jl")
include("..\\distances\\DistanceFunction.jl")
include("..\\instances\\MasterInstance.jl")
include("..\\slave\\ODP_SlaveAlgorithm.jl")
include("..\\net_elements\\DataInfo.jl")
include("..\\net_elements\\IterationResult.jl")
include("..\\utils\\Result.jl")

# Create algorithm master side.
mutable struct ODPMasterAlgorithm
    # MasterInstance
    instance::MasterInstance
    datasetName::String
    # Active points at the current iteration
    currentActivePoints::Int64

    # Heap representing the solution set
    solution::Heap
    # Solving Set
    solvingSet::Set{Point}

    # Set of candidate points
    C::Set{Point}

    # for future implementations
    #port::Int64
    #conntype::String

    # Default costructor
    ODPMasterAlgorithm(cfg::MasterConfiguration) = new(createMasterInstance(cfg),"",0, MinHeap(),Set{Point}(),Set{Point}())
    ODPMasterAlgorithm() = new()
end

# Factory method to create an ODPMasterAlgorithm instance from a configuration
function getMasterAlgorithm(cfg::MasterConfiguration)
    alg = ODPMasterAlgorithm()
    instance = createMasterInstance(cfg)
    alg.instance = instance
    alg.datasetName = ""
    alg.currentActivePoints = 0
    alg.solution = MinHeap()
    alg.solvingSet = Set{Point}()
    alg.C = Set{Point}()
    return alg
end

# Run the SolvingSet algorithm
function runODPAlgorithm(this::ODPMasterAlgorithm, filepath::String)

    # Initializing Result struct
    pltDataset = Array{Vector{Float64},1}()
    pltSolvingSet = Array{Vector{Float64},1}()
    pltSolution = Array{Vector{Float64},1}()
    dist_count = 0

    # Utilities to print current iteration candidate set and current iteration solving set
    iterationC = Array{Vector{Float64},1}()
    itC = []
    iterationS = Array{Vector{Float64},1}()
    itS = []
    iteration_counter = 0

    # Creating SolvingSet and Solution
    this.solvingSet = Set{Point}()
    this.solution = MinHeap(this.instance.n)
    # Init IterationResult
    iterationResult = IterationResult()

    # Lower bound of solution set
    minOut = .0
    # Active points at the current iteration
    currentActivePoints = 0

    # Create ODPSlaveAlgorithm from a path representing the dataset location
    slaveAlg = getSlaveAlgorithm(filepath)

    ### Sending configuration parameters to the Slave ###
    setParam(slaveAlg, this.instance.n, this.instance.k, this.instance.m, this.instance.dist)

    # Filling pltDataset
    for p in slaveAlg.instance.dataset
        push!(pltDataset, p.coords)
    end

    ###### Get dataset info from Slave ###
    datainfo = getDataInfo(slaveAlg)

    # Get dataset name
    datasetName = datainfo.datasetName
    this.datasetName = datainfo.datasetName
    # Get dataset point dimensions
    pointDim = datainfo.pointDim
    this.instance.pointDim = pointDim
    # Get dataset dimension
    datasetSize = datainfo.datasetSize
    currentActivePoints = datasetSize

    # Filling candidate Set with m random samples of dataset
    this.C = initCandidates(slaveAlg)
    #println("Dataset Name: ", string(datasetName), "\nDataset Dimension: ", string(datasetSize), "\nPoint Dimension: ", string(pointDim),"\n\nCurrent active Points = ", string(currentActivePoints))

    while (!isempty(this.C) && currentActivePoints>0)

        # Adding candidate points to the SolvingSet
        union!(this.solvingSet, this.C)

        # Asking the Slave to compute an iteration
        iterationResult = computeIteration(slaveAlg, this.C, minOut)

        # Retrieving data from the IterationResult
        dist_count += iterationResult.distCount
        currentActivePoints = iterationResult.leftActiveData
        #println("\nCurrent active Points = ", string(currentActivePoints))

        # Updating solution heap
        for q in this.C
            updateMax(this.solution,q)
        end

        # Updatinf solution's lower bound
        minOut = top(this.solution.valtree).nn.weight

        # Filling itC and itS for iteration plot
        if iteration_counter in [1,2,3,4]
            for x in this.C
                push!(iterationC, x.coords)
            end
            push!(itC, deepcopy(iterationC))
            setdiff!(iterationC, iterationC)

            for y in this.solvingSet
                push!(iterationS, y.coords)
            end
            push!(itS, deepcopy(iterationS))
            setdiff!(iterationS, iterationS)
        end

        # Extracting nextCand form IterationResult, and setting them as candidate points for the next itersation
        setdiff!(this.C, this.C)
        tmp = extract_all!(iterationResult.nextCand.valtree)
        union!(this.C,tmp)

        #println(iteration_counter)
        iteration_counter +=1

        if isempty(this.C)
            println("TERMINATION CONDITION SATISFIED: No more candiates..")
        end
        if currentActivePoints == 0
            println("TERMINATION CONDITION SATISFIED: No more active points..")
        end
    end

    # Filling pltSolvingSet for plotting
    for x in this.solvingSet
        push!(pltSolvingSet, x.coords)
    end

    # Filling pltSolution for plotting
    for x in extract_all!(this.solution.valtree)
        push!(pltSolution, x.coords)
    end

    res_header = "ODP <"*"D"*", "*string(this.instance.dist)*", "*string(this.instance.n)*", "*string(this.instance.k)*", "*string(this.instance.m)*">"

    res = Result(string(res_header), pltDataset, pltSolvingSet, pltSolution, itC, itS, dist_count, 0)

    return res
end


text:
# This is part of Quiz.jl. License is MIT.

# Interface defined by empty generic functions

"Returns the value of its quiz element argument."
function the end

"Returns an instantiation expression (iE) of its quiz element argument."
function iE end

export the, iE

ext: jl


text:
using MathProgBase

import MathProgBase.LinearQuadraticModel, MathProgBase.loadproblem!, MathProgBase.optimize!, MathProgBase.status, MathProgBase.getobjval, MathProgBase.getsolution, MathProgBase.getunboundedray, MathProgBase.linprog

abstract AbstractPolyhedraModel{N, T} <: MathProgBase.AbstractLinearQuadraticModel

loadproblem!(m::AbstractPolyhedraModel, hrep::HRep, c, sense) = error("loadproblem! not implemented")

include("lphrep.jl")

function loadproblem!(m::AbstractPolyhedraModel, A, l, u, c, lb, ub, sense)
  loadproblem!(m, LPHRepresentation(A, l, u, lb, ub), c, sense)
end
function loadproblem!(m::MathProgBase.AbstractLinearQuadraticModel, lp::LPHRepresentation, c, sense)
  loadproblem!(p, lp.A, lp.l, lp.u, c, lp.lb, lp.ub, sense)
end
function loadproblem!(m::MathProgBase.AbstractLinearQuadraticModel, hrep::HRep, c, sense)
  loadproblem!(p, LPHRepresentation(hrep), c, sense)
end

export LPPolyhedron, SimpleVRepPolyhedraModel

abstract LPPolyhedron{N, T}

type SimpleVRepPolyhedraModel{N, T} <: AbstractPolyhedraModel{N, T}
  vrep::Nullable{VRep{N, T}}
  obj::Nullable{Vector{T}}
  sense::Nullable{Symbol}

  objval::Nullable{T}
  solution::Nullable{Vector{T}}
  status::Symbol

  function SimpleVRepPolyhedraModel()
    new(nothing, nothing, nothing, nothing, nothing, :Undecided)
  end
end

function loadproblem!{N,T}(lpm::SimpleVRepPolyhedraModel{N,T}, vrep::VRep{N,T}, obj::Vector{T}, sense)
  if !(sense in [:Max, :Min])
    error("sense should be :Max or :Min")
  end
  if !myeqzero(obj)
    lpm.obj = copy(obj)
    lpm.sense = sense
  end
end
loadproblem!{N,T}(lpm::SimpleVRepPolyhedraModel{N,T}, vrep::VRep{N,T}, obj::Vector, sense) = loadproblem!(lpm, vrep, Vector{T}(obj), sense)
loadproblem!{N,T}(lpm::SimpleVRepPolyhedraModel{N,T}, vrep::VRep{N}, obj::Vector, sense) = loadproblem!(lpm, VRepresentation{N,T}(vrep), Vector{T}(obj), sense)

function myobjval(obj, v::VRepElement)
  if islin(r)
    if lpm.sense == :Min
      objval = min(mydot(lpm.obj, coord(r)), -mydot(lpm.obj, coord(r)))
    else
      objval = max(mydot(lpm.obj, coord(r)), -mydot(lpm.obj, coord(r)))
    end
  else
    objval = mydot(lpm.obj, coord(r))
  end
end

function optimize!{N, T}(lpm::SimpleVRepPolyhedraModel{N, T})
  if !hasvreps(lpm.vrep)
    lpm.status = :Infeasible
  else
    if lpm.sense in [:Max, :Min]
      better(a, b) = (lpm.sense == :Max ? a > b : a < b)
      mybetter(a, b) = (lpm.sense == :Max ? mygt(a, b) : mylt(a, b))
      lpm.status = :Infeasible
      for r in rays(lpm.vrep)
        objval = myobjval(lpm.obj, r)
        if lpm.status != :Unbounded && mybetter(objval, zero(T))
          lpm.status = :Unbounded
          lpm.objval = lpm.sense == :Max ? typemax(T) : typemin(T)
          lpm.solution = r
        end
      end
      if status != :Unbounded
        for p in points(lpm.vrep)
          objval = myobjval(lpm.obj, p)
          if lpm.status == :Undecided || (lpm.status == :Optimal && better(objval, get(lpm.objval)))
            lpm.status = :Optimal
            lpm.objval = objval
            lpm.solution = p
          end
        end
      end
    else
      lpm.status = :Optimal
      lpm.objval = zero(T)
      lpm.solution = spzeros(T,N)
    end
  end
end

function MathProgBase.status(lpm::SimpleVRepPolyhedraModel)
  lpm.status
end
function MathProgBase.getobjval(lpm::SimpleVRepPolyhedraModel)
  get(lpm.objval)
end
function MathProgBase.getsolution(lpm::SimpleVRepPolyhedraModel)
  copy(get(lpm.solution))
end
function MathProgBase.getunboundedray(lpm::SimpleVRepPolyhedraModel)
  copy(get(lpm.solution))
end

type LinprogSolution
    status
    objval
    sol
    attrs
end

function MathProgBase.linprog{N}(c::Vector, p::Rep{N}, solver::MathProgBase.AbstractMathProgSolver = defaultLPsolverfor(p))
  if N != length(c)
    println("length of objective does not match dimension of polyhedron")
  end
  loadproblem!(m, c, :Min)
  optimize!(m)
  stat = status(m)
  if stat == :Optimal
    return LinprogSolution(stat, getobjval(m), getsolution(m), Dict())
  elseif stat == :Unbounded
    attrs = Dict()
    try
      attrs[:unboundedray] = getunboundedray(m)
    catch
      warn("Problem is unbounded, but unbounded ray is unavailable; check that the proper solver options are set.")
    end
    return LinprogSolution(stat, nothing, [], attrs)
  else
    return LinprogSolution(stat, nothing, [], Dict())
  end
end

function Base.isempty{N,T}(p::Polyhedron{N,T})
  linprog(zeros(T, N), p).status == :Infeasible
end

function ishredundantaux(p::Polyhedron, a, b, strongly, cert, solver)
  sol = linprog(-a, p, solver)
  if sol.status == :Unbounded
    cert ?  (false, sol.attrs[:unboundedray], :UnboundedRay) : false
  elseif sol.status == :Optimal
    if mygt(sol.objval, b)
      cert ? (false, sol.sol, :ExteriorPoint) : false
    elseif mygeq(sol.objval, b)
      if strongly
        cert ? (false, sol.sol, :BoundaryPoint) : false
      else
        cert ? (true, sol.sol, :BoundaryPoint) : true
      end
    else
      cert ? (true, nothing, :NotApplicable) : true
    end
  end
end
function ishredundant(p::Rep, h::HRepElement; strongly=false, cert=false, solver = defaultLPsolverfor(p))
  if islin(h)
    sol = ishredundantaux(p, h.a, h.β, strongly, cert, solver)
    if !sol[1]
      sol
    else
      ishredundantaux(p, -h.a, -h.β, strongly, cert, solver)
    end
  else
    ishredundantaux(p, h.a, h.β, strongly, cert, solver)
  end
end


text:
module MCLightTransport

    using LinearAlgebra: norm
    using StaticArrays

    export Ray, PointSource, makephoton
    export Photon
    export CatesianVolume, RadialVolume
    export simulate!

    include("absorbers.jl")
    include("emitters.jl")
    include("photon.jl")
    include("simulate.jl")

end # module


text:
#= Ornstein-Uhlenbeck process
   without parameter estimation

The Ornstein-Uhlenbeck (OU) process is a continuous-time continuous-space
stochastic process governed by the stochastic differential equation (SDE)

          dY(t) = (Θ_1 - Θ_2 Y(t)) dt + Θ_3 dW(t)

where W(t) is the Wiener process (standard Brownian motion). This OU process is
stationary, Gaussian and Markovian. It is mean-reverting in the sense that it
drifts towards its long-term mean of Θ_1 / Θ_2.

For a step size Δt, the transition distribution can be obtained analytically as
a Normal with mean μ(Δt) and variance σ^2(Δt):

      Y(t + Δt) | Y(t) ~ Normal(μ(Δt), σ^2(Δt))
                 μ(Δt) = Θ_1 / Θ_2 + (Y(t) - Θ_1 / Θ_2) exp(- Θ_2 Δt)
               σ^2(Δt) = Θ_3^2 / (2 * Θ_2) * (1 - exp(- 2 Θ_2 Δt))

Given a sequence of noisy observations of the process Y(t_n) ≈ y_n for
1 <= n <= N, the task is to infer the posterior distribution over each Y(t_n).
This choice is for the sake of simplicity, the example can be modified to
infer posterior distributions over any Y(t).

This example has been inspired by the LibBi example at
http://libbi.org/packages/OrnsteinUhlenbeckBridge.html. =#

using Turing

#= The parameter values of Θ_1, Θ_2, Θ_3 follow [Aït-Sahalia, 1999],
   [Sun et al., 2013] and [Del Moral et al., 2015]. =#
const Θ_1 = 0.0187
const Θ_2 = 0.2610
const Θ_3 = 0.0224

#= [Del Moral, 2015] observe the OU process directly with no observation noise,
and their task is to simulate the diffusion bridges between the observed values.
Our examples uses observation noise with standard deviation σ_y > 0, as this is
required for inference to work at this stage. =#
const σ_y = 0.01

#= The observation times ts and the corresponding observed values.
   It is assumed that the observation times are sorted increasingly. =#
ts      = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
data_ys = [0.0, -0.009114765844660457, 0.029830170268001097,
0.021755657377028295, 0.013934983652697013, 0.01910503829939146,
0.049404191951179434, 0.03253467537333853, 0.07701445971172535,
0.04618793002053633, 0.04001971135990065]
N = length(ts)

# The transition distribution of the OU process
function trans(Yt, Δt)
  mean = Θ_1 / Θ_2 + (Yt - Θ_1 / Θ_2) * exp(- Θ_2 * Δt)
  variance = Θ_3 * Θ_3 / (2 * Θ_2) * (1 - exp(- 2 * Θ_2 * Δt))
  # Normal is initialized with mean and standard deviation, not variance!
  return Normal(mean, sqrt(variance))
end

# The generative model of the OU process
@model ou begin
  # use TArray to support array copying when particles are duplicated!
  Y = TArray(Float64, N)
  Y[1] = 0.0
  for n = 2:N
    Δt = ts[n] - ts[n-1]
    @assume(Y[n] ~ trans(Y[n-1], Δt))
    @observe(data_ys[n] ~ Normal(Y[n], σ_y))
  end
  @predict Y
end

#= This concludes the example. The remainder of this script verifies that
inference works properly by checking that (potentially weighted) samples from
the posterior over Y[N] (the last element of the sequence Y) approximate the
exact posterior over Y[N].

The exact posterior over Y[N] is obtained by recalling that our OU process
is a Gaussian process [Rasmussen and Williams, 2005] with constant mean function
μ(t) = Θ_1 / Θ_2 and stationary covariance function

          k(r) = Θ_3^2 / (2 * Θ_2) * exp(- Θ_2 * |r|)

which is also known as the exponential, Laplace, or Laplacian kernel. (See
Appendix B.2.1 in [Rasmussen and Williams, 2005] for a derivation.) Hence exact
posteriors can be obtained by performing standard GP regression, which involves
inverting a kernel matrix A := (K + σ_y^2 I). =#

k(x, y) = Θ_3^2 / (2 * Θ_2) * exp(- Θ_2 * abs(x - y))
K = [k(x1, x2)::Float64 for x1=ts, x2=ts]
noise_matrix = σ_y^2 * eye(N)
noise_matrix[1, 1] = 0.0 # observed 0.0 at time t=0 with no noise
A = K + noise_matrix
Ainv = inv(A)

# Compute true posterior over Y[N], the value of Y(t) at time t=ts[N]
x_star = ts[N]
k_star = [k(x, x_star)::Float64 for x=ts]
f_star_mean = transpose(k_star) * (A \ (data_ys - Θ_1 / Θ_2)) + Θ_1 / Θ_2
f_star_variance = k(x_star, x_star) - transpose(k_star) * (A \ k_star)
f_star_std = sqrt(f_star_variance)
f_star_posterior = Normal(f_star_mean[1], f_star_std[1])
println("True posterior of last state: ", f_star_posterior)

# KL-divergence between a Normal MLE-fitted to samples and true posterior
function ou_divergence(samples, weights)
  fitted = fit_mle(Normal, samples, weights)
  return Turing.kl(fitted, f_star_posterior)
end
function ou_divergence(samples :: Vector{Float64})
  weights = fill(Float64(1), length(samples))
  return ou_divergence(samples, weights)
end

#= Compute the KL-divergence between moment-matched Normal from samples and
exact posterior, as the number of samples increases =#
num_samples = [10, 30, 100, 300, 1000, 3000, 10000]
samples = sample(ou, SMC(maximum(num_samples)))
samples_lasty = [res[N]::Float64 for res=samples[:Y]]
kl_divergences = map(n -> ou_divergence(samples_lasty[1:n]), num_samples)

# Plot the results using Gadfly
using Gadfly
plot(x=num_samples, y=kl_divergences,
     Scale.x_log10, Scale.y_log10,
     Guide.title("KL-divergence between exact and moment matched posterior"),
     Guide.xlabel("number of samples"), Guide.ylabel("KL-divergence"))

#= References:
- Y. Aït-Sahalia. Transition densities for interest rate and other nonlinear
  diffusions. The Journal of Finance, 54(4):1361–1395, 1999. ISSN 1540-6261.
  doi: 10.1111/0022-1082.00149.
- L. Sun, C. Lee, and J. A. Hoeting. Penalized importance sampling for parameter
  estimation in stochastic differential equations. 2013.
  URL http://arxiv.org/abs/1305.4390.
- Del Moral, Pierre, and Lawrence M. Murray. "Sequential Monte Carlo with highly
  informative observations." SIAM/ASA Journal on Uncertainty Quantification 3,
  no. 1 (2015): 969-997. http://arxiv.org/pdf/1405.4081v2.pdf
- C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning
  (Adaptive Computation and Machine Learning). The MIT Press, 2005.
  ISBN 978-0-262-18253-9. URL http://www.gaussianprocess.org/gpml/.
=#


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
using Printf

# =================================================================================================
"""
    struct Tsv
        level::Int
        page::Int
        block::Int
        paragraph::Int
        line::Int
        word::Int
        left::Int
        top::Int
        width::Int
        height::Int
        conf::Float32
        text::String
    end

This structure holds the details of a line in the TSV formatted text provided by the tesseract
library.

__Values:__

| Name       | Description
| :--------- | :----------
| level      | Identifies what the line describes.
| page       | This is the page number passed into the `tess_tsv()` method.
| block      | Identifies the block on the page.
| paragraph  | The paragraph number in the block.
| line       | The line in the paragraph.
| word       | The word in the line.
| left       | Left edge of the item in pixels.
| top        | Top edge of the item in pixels.
| width      | Width of the item in pixels.
| height     | Height of the item in pixels.
| conf       | How confident the OCR engine is of the word (`0` - `100`). `-1` if level is not 5.
| text       | The word that was decoded from the image.

__Details:__

Level identifies what information the line is providing:

  * `1` - Page information, added at the start of the page.
  * `2` - Block information, added at the start of a block.
  * `3` - Paragraph information, added at the start of a paragraph.
  * `4` - Line information, added at the start of a line.
  * `5` - Word information, identifies a word that was read from the page.

The left, top, width, and height values define a box in pixels that encompases the item.  So if the
level is 1, the box describes the whole image.  If the level is `1`, then the box encloses the block
that was extracted, and so on down to the word that was extracted.

See also: [`tess_parsed_tsv`](@ref)
"""
struct Tsv
    level::Int
    page::Int
    block::Int
    paragraph::Int
    line::Int
    word::Int
    left::Int
    top::Int
    width::Int
    height::Int
    conf::Float32
    text::String
end

# =================================================================================================
"""
    show(
        io::IO,
        data::Tsv
    )::Nothing

Display summary information about the object.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | io   |         | The stream to write the data to.
| R | data |         | The object to display the information for.
"""
function Base.show(
            io::IO,
            data::Tsv
        )::Nothing
    if data.level == 1
        print(io, "Page:  [$(data.page)] ($(data.width), $(data.height))")
    elseif data.level == 2
        print(io, "Block: [$(data.page)/$(data.block)] ($(data.left), $(data.top)) - ($(data.width), $(data.height))")
    elseif data.level == 3
        print(io, "Para:  [$(data.page)/$(data.block)/$(data.paragraph)] ($(data.left), $(data.top)) - ($(data.width), $(data.height))")
    elseif data.level == 4
        print(io, "Line:  [$(data.page)/$(data.block)/$(data.paragraph)/$(data.line)] ($(data.left), $(data.top)) - ($(data.width), $(data.height))")
    elseif data.level == 5
        print(io, "Word:  [$(data.page)/$(data.block)/$(data.paragraph)/$(data.line)/$(data.word)] $(data.text) ($(data.conf)%)")
    else
        print(io, "Unhandled level: $(data.level)")
    end
    nothing
end

# =================================================================================================
"""
    tsv_from(
        line::AbstractString
    )::Tsv

Create a Tsv object from a line of text.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | line |         | The line to extract the information from for the Tsv.

__Details:__

The line must have 12 values separated by tabs.
"""
function tsv_from(
            line::AbstractString
        )::Tsv
    local a = split(line, '\t')

    if length(a) != 12
        throw("Line doesn't contain 12 values $(length(a)): $line")
    end

    Tsv(
        parse(Int, a[1]),
        parse(Int, a[2]),
        parse(Int, a[3]),
        parse(Int, a[4]),
        parse(Int, a[5]),
        parse(Int, a[6]),
        parse(Int, a[7]),
        parse(Int, a[8]),
        parse(Int, a[9]),
        parse(Int, a[10]),
        parse(Float32, a[11]),
        a[12]
        )
end

# =================================================================================================
"""
    tess_parsed_tsv(
        inst::TessInst,
        page::Integer = Int32(1)
    )::Union{Vector{Tsv}, Nothing}

Retrieve the TSV results from an recognition and parse it into a list of Tsv objects.  Returns
nothing if there is an error.

__Arguments:__

| T | Name | Default    | Description
|:--| :--- | :--------- | :----------
| R | inst |            | The Tesseract instance to call.
| O | page | `Int32(1)` | The page to get the data for.

__Example:__

```jldoctest
using Tesseract

download_languages()

instance = TessInst()
pix = sample_pix()

tess_image(instance, pix)
tess_resolution(instance, 72)

tsv = tess_parsed_tsv(instance);

# output

292-element Vector{Tsv}:
 Page:  [1] (500, 600)
 Block: [1/1] (10, 9) - (479, 514)
 Para:  [1/1/1] (11, 9) - (406, 14)
 Line:  [1/1/1/1] (11, 9) - (406, 14)
 Word:  [1/1/1/1/1] No (95.79193%)
 Word:  [1/1/1/1/2] one (95.79193%)
 Word:  [1/1/1/1/3] would (92.95379%)
 Word:  [1/1/1/1/4] have (92.95379%)
 Word:  [1/1/1/1/5] believed (96.81915%)
 Word:  [1/1/1/1/6] in (96.770004%)
 ⋮
 Word:  [1/1/7/3/5] twentieth (96.62509%)
 Word:  [1/1/7/3/6] century (96.81123%)
 Word:  [1/1/7/3/7] came (96.81123%)
 Word:  [1/1/7/3/8] the (96.72949%)
 Word:  [1/1/7/3/9] great (96.72949%)
 Para:  [1/1/8] (11, 509) - (172, 14)
 Line:  [1/1/8/1] (11, 509) - (172, 14)
 Word:  [1/1/8/1/1] great (86.34806%)
 Word:  [1/1/8/1/2] disillusionment. (86.34806%)
```

See also: [`tess_tsv`](@ref), [`tess_confidences`](@ref)
"""
function tess_parsed_tsv(
            inst::TessInst,
            page::Integer = Int32(1)
        )::Union{Vector{Tsv}, Nothing}
    local text = tess_tsv(inst, page)

    if text == nothing
        return nothing
    end

    local retval = Vector{Tsv}()
    for line in split(text, '\n')
        if isempty(line) == false
            push!(retval, tsv_from(line))
        end
    end

    return retval
end


text:
function method_IFP(kx, ky, n, percentile_cutoff, rho; verbose=true)

    verbose && @info "generating data"
    data = gendata(kx, ky, n, percentile_cutoff, rho)
    verbose && @info "data has size $(size(data.gen)) with $(count(ismissing, data.obs)) missings"

    A = data.gen
    A_obs = @. ifelse(ismissing(data.obs), 0, data.obs)

    initial_matrix = rand(1:1, size(data.gen)[1], size(data.gen)[2])

    real_row_sums = [sum(data.gen[i,:]) for i in 1:size(data.gen)[1]]
    real_col_sums = [sum(data.gen[:,j]) for j in 1:size(data.gen)[2]]

    fixed_row_sums = [sum((A_obs)[i,:]) for i in 1:size(data.gen)[1]]
    fixed_col_sums = [sum((A_obs)[:,j]) for j in 1:size(data.gen)[2]]

    A_over_threshold = sign.(A_obs)
    A_under_threshold = 1 .- A_over_threshold

    test_max = 100

    while test_max >= 1.0001

        current_row_sums = [sum((initial_matrix .* A_under_threshold)[i,:]) for i in 1:size(data.gen)[1]]
        row_multiplier = (real_row_sums - fixed_row_sums) ./ current_row_sums
        replace!(row_multiplier, NaN => 1)
        initial_matrix = initial_matrix .* row_multiplier

        current_col_sums = [sum((initial_matrix .* A_under_threshold)[:,j]) for j in 1:size(data.gen)[2]]
        col_multiplier = (real_col_sums - fixed_col_sums) ./ current_col_sums
        replace!(col_multiplier, NaN => 1)
        initial_matrix = transpose(transpose(initial_matrix) .* col_multiplier)

        current_row_sums = [sum((initial_matrix .* A_under_threshold .+ data.gen .* A_over_threshold)[i,:]) for i in 1:size(data.gen)[1]]
        row_multiplier = real_row_sums ./ current_row_sums
        current_col_sums = [sum((initial_matrix .* A_under_threshold .+ data.gen .* A_over_threshold)[:,j]) for j in 1:size(data.gen)[2]]
        col_multiplier = real_col_sums ./ current_col_sums

        test_max = maximum(vcat(row_multiplier, col_multiplier))
    end

    A_obs = initial_matrix .* A_under_threshold .+ A_obs

    mae = mean_absolute_error(data, A_obs)

    verbose && @info "mae: $mae"

    return(data=data, blanks=count(ismissing, data.obs), mae=mae, vc = cramer_v(data.gen), fit = A_obs)

end

max_stars_repo_name: Michal-Kot/FillMissingEntries

max_issues_repo_name: Michal-Kot/FillMissingEntries

max_forks_repo_name: Michal-Kot/FillMissingEntries


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

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

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

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

include("$path/../setup.jl") # for types
rc("figure", figsize=(3.5, 3.20))
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.27, bottom=0.2, top=0.8, right=.95)
rc("figure.subplot", left=0.22, bottom=0.32, top=0.99, right=.95)
rc("axes", color_cycle=["348ABD", "A60628", "009E73", "7A68A6", "D55E00", "CC79A7"])
plot_cycle=["#348ABD", "#A60628", "#009E73", "#7A68A6", "#D55E00", "#CC79A7"]
close("all")

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

cr_factor = ([50,100, 150])

cr_factor_plot = ([">22","45", "67"])

db_factor = 76 #[1,2,8,32,128]

DETAILED = false

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

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

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

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

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

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

ETAs_prop_accel = zeros(length(numprops),length(cr_factor))
ETAs_prop_climb = zeros(length(numprops),length(cr_factor))
ETAs_prop_cruise = zeros(length(numprops),length(cr_factor))

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

MotorMasses_accel = zeros(length(numprops),length(cr_factor))
MotorMasses_climb = zeros(length(numprops),length(cr_factor))
MotorMasses_cruise = zeros(length(numprops),length(cr_factor))

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

Power_IN_accel = zeros(length(numprops),length(cr_factor))
Power_IN_climb = zeros(length(numprops),length(cr_factor))
Power_IN_cruise = zeros(length(numprops),length(cr_factor))

Power_Shaft_accel = zeros(length(numprops),length(cr_factor))
Power_Shaft_climb = zeros(length(numprops),length(cr_factor))
Power_Shaft_cruise = zeros(length(numprops),length(cr_factor))

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

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

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

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

R_transition = zeros(length(numprops),length(cr_factor))

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

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

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

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

t_accel = zeros(length(numprops),length(cr_factor))
t_transition = zeros(length(numprops),length(cr_factor))
t_climb = zeros(length(numprops),length(cr_factor))
t_cruise = zeros(length(numprops),length(cr_factor))

mass_prop = zeros(length(numprops),length(cr_factor))

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

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

blades_prop = zeros(length(numprops),length(cr_factor))

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

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

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

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

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

        total_weight[i,j] = detailed_system.total_weight

        R_transition[i,j] = detailed_system.radius

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

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

        ETAs_prop_accel[i,j] = detailed_output_accel.ETA_p
        ETAs_prop_climb[i,j] = detailed_output_climb.ETA_p
        ETAs_prop_cruise[i,j] = detailed_output_cruise.ETA_p

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

        MotorMasses_accel[i,j] = detailed_output_accel.mass_motor
        MotorMasses_climb[i,j] = detailed_output_climb.mass_motor
        MotorMasses_cruise[i,j] = detailed_output_cruise.mass_motor

        Rtips[i,j] = design_cruise.p_Rtip

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

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

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

        Power_IN_accel[i,j] = detailed_output_accel.V_m*detailed_output_accel.I_m*numprops[i]
        Power_IN_climb[i,j] = detailed_output_climb.V_m*detailed_output_climb.I_m*numprops[i]
        Power_IN_cruise[i,j] = detailed_output_cruise.V_m*detailed_output_cruise.I_m*numprops[i]

        Power_Shaft_accel[i,j] = Power_IN_accel[i,j] * detailed_output_accel.ETA_m
        Power_Shaft_climb[i,j] = Power_IN_climb[i,j] * detailed_output_climb.ETA_m
        Power_Shaft_cruise[i,j] = Power_IN_cruise[i,j] * detailed_output_cruise.ETA_m

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

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

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

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

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

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

        t_accel[i,j] = detailed_system.t_accel
        t_transition[i,j] = detailed_system.t_transition
        t_climb[i,j] = detailed_system.t_climb
        t_cruise[i,j] = detailed_system.t_cruise

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

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

        blades_prop[i,j] = design_cruise.blades #same for all three

        if DETAILED
            splchord = Dierckx.Spline1D(design_climb.w_nondim_halfspany,design_climb.w_chord_pts,k = 1)
            splchordpts = splchord(detailed_output_climb.panels.y/maximum(detailed_output_climb.panels.y))
            figname = "cl_p$(numprops[i])_cr$(cr_factor[j])"
            figure(figname)
            plot(detailed_output_climb.panels.y/maximum(detailed_output_climb.panels.y),detailed_output_climb.cllocal,".-",label = "local cl")
            plot(detailed_output_climb.panels.y/maximum(detailed_output_climb.panels.y),detailed_output_climb.cl,".-",label = "freestream cl")
            xlabel("wing halfspan")
            ylabel("cl")
            # legend(loc="center left", bbox_to_anchor=(1, 0.5))
            savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

            figname = "lift_p$(numprops[i])_cr$(cr_factor[j])"
            figure(figname)
            # plot(detailed_output_climb.panels.y/maximum(detailed_output_climb.panels.y),detailed_output_climb.cllocal,".-",label = "local cl")
            plot(detailed_output_climb.panels.y/maximum(detailed_output_climb.panels.y),detailed_output_climb.cl*0.5*1.225*design_climb.velocity^2.0.*splchordpts,".-",label = "freestream cl")
            xlabel("wing halfspan")
            ylabel("local lift")
            # legend(loc="center left", bbox_to_anchor=(1, 0.5))
            savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

            figname = "localq_p$(numprops[i])_cr$(cr_factor[j])"
            figure(figname)
            plot(detailed_output_climb.panels.y/maximum(detailed_output_climb.panels.y),0.5*1.225*detailed_output_climb.Vinfeff.^2,".-",label = "local q")
            xlabel("wing halfspan")
            ylabel("q")
            # legend(loc="center left", bbox_to_anchor=(1, 0.5))
            savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

            figname = "chord_p$(numprops[i])_cr$(cr_factor[j])"
            figure(figname)
            plot(design_climb.w_nondim_halfspany*design_climb.w_halfspan,design_climb.w_chord_pts,".-",label = "local cl")
            xlabel("wing halfspan (meters)")
            ylabel("chord (meters)")

            # legend(loc="center left", bbox_to_anchor=(1, 0.5))
            # ylim([0.1.2*0,maximum(design_climb.maxy1])
            savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

            figname = "geometricaoa_p$(numprops[i])_cr$(cr_factor[j])"
            figure(figname)
            plot(linspace(0,1,length(detailed_output_climb.alphaeff)),design_climb.w_aoa*ones(detailed_output_climb.alphaeff)*180/pi,".-",label = "geometric aoa")
            plot(linspace(0,1,length(detailed_output_climb.alphaeff)),detailed_output_climb.alphaeff*180/pi,".-",label = "effective aoa")
            plot(linspace(0,1,length(detailed_output_climb.alphaeff)),(detailed_output_climb.alphaeff-design_climb.w_aoa)*180/pi,".-",label = "effective - geometric aoa")
            xlabel("wing halfspan")
            ylabel("aoa (deg)")
            legend(loc="center left", bbox_to_anchor=(1, 0.5))
            savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

        end

    # end
    end
end

# High level results
figname = "Takeoff_Speed"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(takeoff_dist[:,i])
    plot(numprops,takeoff_dist[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Takeoff Distance (m)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Accel_Speed"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(dist_accel[:,i])
    plot(numprops,dist_accel[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Ground Roll Distance (m)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Range_Props"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(ranges[:,i]/1000.0)
    plot(numprops,ranges[:,i]/1000.0,"o-",label="$(cr_factor_plot[i])")#"1/$(cr_factor[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Range (km)")
# ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "span"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(span[:,i])
    plot(numprops,span[:,i],"o-",label="$(cr_factor_plot[i])")#"1/$(cr_factor[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Span (m)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "CL_climb"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(CLs_climb[:,i])
    plot(numprops,CLs_climb[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Wing Lift Coefficient (CL)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "CL_cruise"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(CLs_cruise[:,i])
    plot(numprops,CLs_cruise[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Wing Lift Coefficient (CL)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "CD_cruise"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(CDs_cruise[:,i])
    plot(numprops,CDs_cruise[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Wing Drag Coefficient (CD)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Noise_climb"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(DBs_climb[:,i])
    plot(numprops,DBs_climb[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Sound Pressure Level (dBa)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "FlightPathAngle"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(gama[:,i]*180/pi)
    plot(numprops,gama[:,i]*180/pi,"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Flight Path Angle (degrees)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Total_Mass"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(total_weight[:,i])
    plot(numprops,total_weight[:,i]/9.81,"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Total Mass (kg)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Propulsion_Mass"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(Propulsion_Mass[:,i])
    plot(numprops,Propulsion_Mass[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Propulsion Mass (kg)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Motor_Mass"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)

    maxy[i] = maximum(MotorMasses_climb[:,i].*numprops)
    plot(numprops,MotorMasses_climb[:,i].*numprops,"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Motor Mass (kg)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Rtip"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(Rtips[:,i])
    plot(numprops,Rtips[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Propeller Radius (m)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Battery_Cruise"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(BattMasses_cruise[:,i])
    plot(numprops,BattMasses_cruise[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Battery Mass (kg)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Battery_Total"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(BattMasses_cruise[:,i])
    plot(numprops,BattMasses_cruise[:,i]+BattMasses_accel[:,i]+BattMasses_climb[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Battery Mass (kg)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Battery_Mass_Fraction"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum((BattMasses_cruise[:,i]+BattMasses_accel[:,i]+BattMasses_climb[:,i])./(total_weight[:,i]/9.81)*100)
    plot(numprops,(BattMasses_cruise[:,i]+BattMasses_accel[:,i]+BattMasses_climb[:,i])./(total_weight[:,i]/9.81)*100,"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Battery Mass Fraction (%)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Propulsion_ETA_Cruise"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(ETAs_cruise[:,i]*100.0)
    plot(numprops,ETAs_cruise[:,i]*100.0,"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Propulsion Efficiency (%)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Propulsion_ETA_Accel"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(ETAs_accel[:,i]*100.0)
    plot(numprops,ETAs_accel[:,i]*100.0,"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Propulsion Efficiency (%)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Propulsion_ETA_Climb"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(ETAs_climb[:,i]*100.0)
    plot(numprops,ETAs_climb[:,i]*100.0,"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Propulsion Efficiency (%)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Propeller_ETA_Climb"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(ETAs_prop_climb[:,i]*100.0)
    plot(numprops,ETAs_prop_climb[:,i]*100.0,"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Propeller Efficiency (%)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Propeller_ETA_Cruise"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(ETAs_prop_cruise[:,i]*100.0)
    plot(numprops,ETAs_prop_cruise[:,i]*100.0,"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Propeller Efficiency (%)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Climb_Thrust"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(Thrust_climb[:,i])
    plot(numprops,Thrust_climb[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Total Thrust (N)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Climb_Thrust2Weight"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(Thrust_climb[:,i])
    plot(numprops,Thrust_climb[:,i]./(total_weight[:,i]),"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Thrust/Weight")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Climb_NetThrust2Weight"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum((Thrust_climb[:,i]-Drags_climb[:,i])./(total_weight[:,i]))
    plot(numprops,(Thrust_climb[:,i]-Drags_climb[:,i])./(total_weight[:,i]),"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Net-Thrust/Weight")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Transition_Radius"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(R_transition[:,i])
    plot(numprops,R_transition[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Transition Radius (m)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Takeoff_Energy_Percent"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    to_percent = 100.0*(BattMasses_accel[:,i]+BattMasses_climb[:,i])./(BattMasses_accel[:,i]+BattMasses_climb[:,i]+BattMasses_cruise[:,i])
    semilogy(numprops,to_percent,"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Takeoff Energy (%)")
# ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Velocity_cruise"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(Velocity_cruise[:,i])
    plot(numprops,Velocity_cruise[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Cruise Velocity (m/s)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Velocity_cruiseMPH"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(Velocity_cruise[:,i])
    plot(numprops,Velocity_cruise[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Cruise Velocity (mph)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "time_accel"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(t_accel[:,i])
    plot(numprops,t_accel[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Acceleration Time (s)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "time_transition"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(t_transition[:,i])
    plot(numprops,t_transition[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Transition Time (s)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "time_climb"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(t_climb[:,i])
    plot(numprops,t_climb[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Climb Time (s)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "time_tclimb"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(t_climb[:,i]+t_transition[:,i])
    plot(numprops,t_climb[:,i]+t_transition[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Transition \& Climb Time (s)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "time_cruise"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(t_cruise[:,i])
    plot(numprops,t_cruise[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Cruise Time (s)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "tip_mach_climb"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(tip_mach_climb[:,i])
    plot(numprops,tip_mach_climb[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Tip Speed (mach)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "aoa_cruise"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(aoa_cruise[:,i])
    plot(numprops,aoa_cruise[:,i]*180/pi,"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Wing Angle of Attack (deg)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "aoa_accel"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(aoa_accel[:,i])
    plot(numprops,aoa_accel[:,i]*180/pi,"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Wing Angle of Attack (deg)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "aoa_climb"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(aoa_climb[:,i])
    plot(numprops,aoa_climb[:,i]*180/pi,"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Wing Angle of Attack (deg)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "prop_mean_chord"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(mean_chord[:,i])
    plot(numprops,mean_chord[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Mean Propeller Chord (m)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "mass_prop"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(mass_prop[:,i])
    plot(numprops,mass_prop[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Propeller Mass (kg)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "blades_prop"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(blades_prop[:,i])
    plot(numprops,blades_prop[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Propeller Blades (#)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Power_climb"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(Power_climb[:,i])
    plot(numprops,Power_climb[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Climb Power Out (watts)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Power_cruise"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(Power_cruise[:,i])
    plot(numprops,Power_cruise[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Cruise Power Out (watts)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Power_accel"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(Power_accel[:,i])
    plot(numprops,Power_accel[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Acceleration Power Out (watts)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Power_IN_climb"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(Power_IN_climb[:,i])
    plot(numprops,Power_IN_climb[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Climb Power IN (watts)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Power_IN_cruise"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(Power_IN_cruise[:,i])
    plot(numprops,Power_IN_cruise[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Cruise Power IN (watts)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Power_IN_accel"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(Power_IN_accel[:,i])
    plot(numprops,Power_IN_accel[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Acceleration Power IN (watts)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Power_Shaft_climb"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(Power_Shaft_climb[:,i])
    plot(numprops,Power_Shaft_climb[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Climb Power Shaft (watts)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Power_Shaft_cruise"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(Power_Shaft_cruise[:,i])
    plot(numprops,Power_Shaft_cruise[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Cruise Power Shaft (watts)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Power_Shaft_accel"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(Power_Shaft_accel[:,i])
    plot(numprops,Power_Shaft_accel[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Acceleration Power Shaft (watts)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Velocity_accel"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(Velocity_accel[:,i])
    plot(numprops,Velocity_accel[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Liftoff Speed (m/s)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Velocity_climb"
figure(figname)
maxy = zeros(length(cr_factor))
for i=1:length(cr_factor)
    maxy[i] = maximum(Velocity_climb[:,i])
    plot(numprops,Velocity_climb[:,i],"o-",label="$(cr_factor_plot[i])")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Liftoff Speed (m/s)")
ylim([0,1.05*ceil(Int,maximum(maxy))])
# legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, 1.30))#loc="center left", bbox_to_anchor=(1, 0.5))
legend(title = "Cruise Constraint (m/s)",loc="upper center", ncol = length(cr_factor), bbox_to_anchor=(0.5, -.20))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Relative_Distance_cr100"
figure(figname)
for i=2#1:length(cr_factor)
    plot(numprops,dist_accel[:,i],"o-",label="Accel")
    plot(numprops,dist_trans[:,i],"o-",label="Transition")
    plot(numprops,dist_climb[:,i],"o-",label="Climb")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Distance (m)")
# ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(loc="center left", bbox_to_anchor=(1, 0.5))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)

figname = "Percent_Relative_Distance_cr100"
figure(figname)
for i=2#1:length(cr_factor)
    plot(numprops,dist_accel[:,i]./takeoff_dist[:,i]*100,"o-",label="Accel")
    plot(numprops,dist_trans[:,i]./takeoff_dist[:,i]*100,"o-",label="Transition")
    plot(numprops,dist_climb[:,i]./takeoff_dist[:,i]*100,"o-",label="Climb")
end
xticks(numprops)
xlabel("Propellers (#)")
ylabel("Relative Distance (%)")
# ylim([0,1.05*ceil(Int,maximum(maxy))])
legend(loc="center left", bbox_to_anchor=(1, 0.5))
savefig("$path/../../MDAO_Paper/figures/$figname.pdf",transparent = true)


text:
# Load Julia packages (libraries) needed  for the snippets in chapter 0

using StatisticalRethinking, CmdStan
#gr(size=(600,600));

# CmdStan uses a tmp directory to store the output of cmdstan

ProjDir = rel_path("..", "scripts", "04")
cd(ProjDir)

# ### snippet 4.7

howell1 = CSV.read(rel_path("..", "data", "Howell1.csv"), delim=';')
df = convert(DataFrame, howell1);

# Use only adults

df2 = filter(row -> row[:age] >= 18, df);
mean_weight = mean(df2[:weight]);
df2[:weight_c] = df2[:weight] .- mean_weight;
first(df2, 5)

# Define the Stan language model

weightsmodel = "
data {
 int < lower = 1 > N; // Sample size
 vector[N] height; // Predictor
 vector[N] weight; // Outcome
}

parameters {
 real alpha; // Intercept
 real beta; // Slope (regression coefficients)
 real < lower = 0 > sigma; // Error SD
}

model {
 height ~ normal(alpha + weight * beta , sigma);
}

generated quantities {
}
";

# Define the Stanmodel and set the output format to :mcmcchains.

stanmodel = Stanmodel(name="weights", model=weightsmodel,
  output_format=:mcmcchains);

# Input data for cmdstan

heightsdata = Dict("N" => length(df2[:height]), "height" => df2[:height], "weight" => df2[:weight_c]);

# Sample using cmdstan

rc, chn, cnames = stan(stanmodel, heightsdata, ProjDir, diagnostics=false,
  summary=false, CmdStanDir=CMDSTAN_HOME);

# Describe the draws

MCMCChains.describe(chn)

# Plot the density of posterior draws

plot(chn)

# Plot regression line using means and observations

scatter(df2[:weight_c], df2[:height], lab="Observations",
  ylab="height [cm]", xlab="weight[kg]")
xi = -16.0:0.1:18.0
rws, vars, chns = size(chn)
alpha_vals = convert(Vector{Float64}, reshape(chn.value[:, 1, :], (rws*chns)));
beta_vals = convert(Vector{Float64}, reshape(chn.value[:, 2, :], (rws*chns)));
yi = mean(alpha_vals) .+ mean(beta_vals)*xi;
plot!(xi, yi, lab="Regression line")

# End of `clip-43s.jl`


text:
#==========================================================================================#

# PARSE WEIGHTS

for W in (AnalyticWeights, FrequencyWeights, ProbabilityWeights, Weights)
    @eval begin
        function update_weights(w::$W{T, Bool}, nonmissing::BitVector) where T
            UnitWeights{Float64}(sum(nonmissing))
        end
    end
end

function update_weights(w::FrequencyWeights, nonmissing::BitVector)
    v = w[nonmissing]
    s = Int(sum(v))
    FrequencyWeights(float(v), s)
end

function update_weights(w::UnitWeights, nonmissing::BitVector)
    UnitWeights{Float64}(sum(nonmissing))
end

function update_weights(w::AbstractWeights, nonmissing::BitVector)
    n  = sum(nonmissing)
    v  = float(w[nonmissing])
    s  = n / sum(v)
    v .= v .* s
    Weights(v, n)
end

#==========================================================================================#

# COPY

function Base.copy(corr::C) where {C <: CorrStructure}
    C((copy(getfield(corr, f)) for f in fieldnames(C))...)
end

# EQUALITY

function Base.isequal(corr₁::C, corr₂::C) where {C <: CorrStructure}
    mapreduce(i -> isequal(getfield.([corr₁, corr₂], i)...), *, fieldnames(C), init = true)
end

Base.isequal(corr₁::CorrStructure, corr₂::CorrStructure) = false

# INTERSECTION BETWEEN A FRAME AND CORRELATION STRUCTURE

update_corr(corr::Homoscedastic, nonmissing::BitVector)   = copy(corr)
update_corr(corr::Heteroscedastic, nonmissing::BitVector) = copy(corr)

function update_corr(corr::Clustered, nonmissing::BitVector)

    isequal(nonmissing, corr.nonmissing) && (return copy(corr))

    nonmissing    .*= corr.nonmissing
    touse           = nonmissing[corr.nonmissing]
    new_id_clusters = corr.id_clusters[touse]
    iter            = sort(unique(corr.id_clusters))
    new_iter        = sort(unique(new_id_clusters))
    new_n_clusters  = length(new_iter)
    new_matrix      = corr.matrix[findall((in)(new_iter), iter), touse]

    Clustered(corr.corrected, nonmissing, new_matrix, new_id_clusters, new_n_clusters)
end

function update_corr(corr::CrossCorrelated, nonmissing::BitVector)

    isequal(nonmissing, corr.nonmissing) && (return copy(corr))

    nonmissing .*= corr.nonmissing
    touse        = nonmissing[corr.nonmissing]
    new_matrix   = SuiteSparse.CHOLMOD.Sparse(Symmetric(corr.matrix[touse, touse]))

    CrossCorrelated(corr.corrected, nonmissing, new_matrix)
end

#==========================================================================================#

# RETRIEVAL

function getvector(obj::Microdata, x::Symbol)
    (length(obj.mapping[x]) == 1) || throw(ArgumentError(string(x) * " is not a vector"))
    view(obj.matrix, :, obj.mapping[x]...)
end

function getmatrix(obj::Microdata, args...)
    x = mapreduce(i -> obj.mapping[i], vcat, args, init = Vector{Int64}())
    view(obj.matrix, :, x)
end

function getnames(obj::Microdata, args...)
    n = coefnames(eterms(obj.model))
    x = mapreduce(i -> obj.mapping[i], vcat, args, init = Vector{Int64}())
    n[x]
end

getcorr(obj::Microdata)       = obj.corr
getnonmissing(obj::Microdata) = obj.nonmissing
getweights(obj::Microdata)    = obj.weights

max_stars_repo_name: lbittarello/Microeconometrics

max_issues_repo_name: lbittarello/Microeconometrics

max_forks_repo_name: lbittarello/Microeconometrics


text:
# Copyright 2021 Gravifer
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
#=
  Author: Gravifer, X. H. Zhu, H. J. Yan, H. Y. Sun
  Date: 2021-04-10
  Description: TBD
=#

module transfinite_interpolate

using StaticArrays
using OffsetArrays
include("../../src/misc-util.jl")
using .__CFD2021__misc_util__: tuplejoin

"""
    transfinite_interpolate_2d(f_lo, f_hi[, v_cr])

Return a interpolating function defined on ``[0, 1]²``.

**Notice**: does not check or warn when extrapolating.

See also: https://en.wikipedia.org/wiki/Transfinite_interpolation

TODO: generalize to higher dimensions.

# Arguments
- `f_lo::Tuple{Function, Function}`: a pair of functions.
- `f_hi::Tuple{Function, Function}`: a pair of functions.
- `v_cr::SMatrix{2, 2} = nothing`: (optional) values to take at the corners;
    use average of edge functions if not provided.
    Still, in practical use you should make sure the ends of the edges meet.

# Output
- `interpolating_function(u::NTuple{2, T} where T <: Number):: Number`:
    Take in a length-2 static vector of numbers and return a number.

# Examples
```jldoctest
julia> interp_func = transfinite_interpolate_2d( (x->2x, y->3y), (x->4x+3, y->2+5y) )
(::var"#interpolating_function#1"{Tuple{var"#2#6", var"#3#7"}, Tuple{var"#4#8", var"#5#9"}}) (generic function with 1 method)

julia> interp_func((0.5,0.5))
3.0
```
"""
function transfinite_interpolate_2d(f_lo::Tuple{Function,Function},
                                    f_hi::Tuple{Function,Function};
                                    v_cr::Union{Nothing,SMatrix{2,2}}=nothing)::Function
    #= handling of optionals suggested by
       https://stackoverflow.com/questions/42499528/julia-convention-for-optional-arguments;
       though Base.Void is deprecated and now called Nothing instead, as pointed out by
       https://discourse.julialang.org/t/void-not-defined-version-problem/14183 =#
    if typeof(v_cr) <: Nothing
        v_cr = SA[(f_lo[1](0) + f_lo[2](0))/2 (f_lo[1](1) + f_hi[2](0))/2
                  (f_hi[1](0) + f_lo[2](1))/2 (f_hi[1](1) + f_hi[2](1))/2]
    end

    function interpolating_function(u::NTuple{2,T} where {T<:Number})::Number
        #= linear-part from edges   bilinear-part from vertices =#
        (1-u[2])*f_lo[1](u[1]) - (1-u[1])*(1-u[2])*v_cr[1,1] +
        (  u[2])*f_hi[1](u[1]) - (  u[1])*(1-u[2])*v_cr[1,2] +
        (1-u[1])*f_lo[2](u[2]) - (1-u[1])*(  u[2])*v_cr[2,1] +
        (  u[1])*f_hi[2](u[2]) - (  u[1])*(  u[2])*v_cr[2,2]
    end
    return interpolating_function
end

"""
    transfinite_interpolate_2d(v_lo, v_hi[, v_cr])
    transfinite_interpolate_2d!(interpolated_array, v_lo, v_hi[, v_cr])

Return a interpolated array defined on a discretization of ``[0, 1]²``.

**Notice**: the in-place version does not check dimension-compatibility of its inputs.

TODO: generalize to higher dimensions.

# Arguments
- `v_lo::Tuple{Vector{<: Real}, Vector{<: Real}}`: a pair of vectors.
- `v_hi::Tuple{Vector{<: Real}, Vector{<: Real}}`: a pair of vectors.
- `v_cr::SMatrix{2, 2, <: Real} = nothing`: (optional) values to take at the corners;
    use average of edge values if not provided.
    Still, in practical use you should make sure the ends of the edges meet.

# Output
- `interpolated_array::Array{Real, 2}`

# Examples
```jldoctest
julia> transfinite_interpolate_2d(([0,0.5,1], [0,0.5,1]), ([1,1.5,2], [1,1.5,2]))
3×3 Matrix{Real}:
 0.0  0.5  1.0
 0.5  1.0  1.5
 1.0  1.5  2.0
```
"""
function transfinite_interpolate_2d( # for regular Arrays
    v_lo::Tuple{Vector{T}, Vector{S}} where {T <: Real, S <: Real},
    v_hi::Tuple{Vector{T}, Vector{S}} where {T <: Real, S <: Real},
    v_cr::Union{Nothing, SMatrix{2, 2, T} where T <: Real} = nothing)::Array{Real, 2}

    #= check dimension-compatibility before allocating anything =#
    length(v_lo) == length(v_hi) || throw(
        DimensionMismatch(" incompatible grid specification ($(length(v_lo))-d at one end but $(length(v_hi))-d at the other)."))
    N = length(v_lo); errmsg = "" #= In this function, dimensions are always indexed from 1, but not so for each axis of them =#
    for dim in eachindex(v_lo)  if ! (axes(v_lo[dim]) == axes(v_hi[dim]))
        errmsg *= "grid mismatch in dimension $(dim) ($(axes(v_lo[dim])) at one end but $(axes(v_hi[dim])) at the other)."
    end end; errmsg == "" || throw(DimensionMismatch(errmsg))

    interpolated_array = Array{Float64}(undef, (length(v) for v in v_lo)...)
    transfinite_interpolate_2d!(interpolated_array, v_lo, v_hi, v_cr)
    return interpolated_array
end
function transfinite_interpolate_2d( # for OffsetArrays
    v_lo::Tuple{OffsetArray, OffsetArray},
    v_hi::Tuple{OffsetArray, OffsetArray},
    v_cr::Union{Nothing, SMatrix{2, 2, T} where T <: Real} = nothing)::OffsetArray

    #= check dimension-compatibility before allocating anything =#
    length(v_lo) == length(v_hi) || throw(
        DimensionMismatch(" incompatible grid specification ($(length(v_lo))-d at one end but $(length(v_hi))-d at the other)."))
    N = length(v_lo); errmsg = "" #= In this function, dimensions are always indexed from 1, but not so for each axis of them =#
    for dim in eachindex(v_lo)  if ! (axes(v_lo[dim]) == axes(v_hi[dim]))
        errmsg *= "grid mismatch in dimension $(dim) ($(axes(v_lo[dim])) at one end but $(axes(v_hi[dim])) at the other)."
    end end; errmsg == "" || throw(DimensionMismatch(errmsg))

    interpolated_array = Array{Float64}(undef, (length(v) for v in v_lo)...)
    interpolated_array = OffsetArray(interpolated_array, tuplejoin((axes(v) for v in v_lo)...))
    transfinite_interpolate_2d!(interpolated_array, v_lo, v_hi, v_cr)
    return interpolated_array
end
function transfinite_interpolate_2d!( # for regular Arrays
    interpolated_array::Array{T, 2}   where  T <: Real,
    v_lo::Tuple{Vector{T}, Vector{S}} where {T <: Real, S <: Real},
    v_hi::Tuple{Vector{T}, Vector{S}} where {T <: Real, S <: Real},
    v_cr::Union{Nothing, SMatrix{2, 2, T} where T <: Real} = nothing)

    if typeof(v_cr) <: Nothing # e.g, the following lines still assumes range ``1:2``.
        v_cr = SA[(v_lo[1][begin]+v_lo[2][begin])/2 (v_lo[1][ end ]+v_hi[2][begin])/2
                  (v_hi[1][begin]+v_lo[2][ end ])/2 (v_hi[1][ end ]+v_hi[2][ end ])/2]
    end; M = [length(v)-1 for v in v_lo] # this is only for scaling, and never for iterating

    # loop Version
    for u in CartesianIndices(interpolated_array) #= #= Deprecated =# Iterators.product(IndexCartesian(), [eachindex(v) for v in v_lo]...) =#
        interpolated_array[u] =
        #= linear-part   from   edges      bilinear-part   from   vertices    =#
        (1-(u[2]-1)/M[2])*v_lo[1][u[1]] - (1-(u[1]-1)/M[1])*(1-(u[2]-1)/M[2])*v_cr[1,1] +
        (  (u[2]-1)/M[2])*v_hi[1][u[1]] - (  (u[1]-1)/M[1])*(1-(u[2]-1)/M[2])*v_cr[1,2] +
        (1-(u[1]-1)/M[1])*v_lo[2][u[2]] - (1-(u[1]-1)/M[1])*(  (u[2]-1)/M[2])*v_cr[2,1] +
        (  (u[1]-1)/M[1])*v_hi[2][u[2]] - (  (u[1]-1)/M[1])*(  (u[2]-1)/M[2])*v_cr[2,2]
    end

    return interpolated_array
end
function transfinite_interpolate_2d!(interpolated_array::OffsetArray,
                                     v_lo::Tuple{OffsetArray,OffsetArray},
                                     v_hi::Tuple{OffsetArray,OffsetArray},
                                     v_cr::Union{Nothing,SMatrix{2,2,T} where T<:Real}=nothing)
    if typeof(v_cr) <: Nothing # e.g, the following lines still assumes range ``1:2``.
        v_cr = SA[(v_lo[1][begin]+v_lo[2][begin])/2 (v_lo[1][ end ]+v_hi[2][begin])/2
                  (v_hi[1][begin]+v_lo[2][ end ])/2 (v_hi[1][ end ]+v_hi[2][ end ])/2]
    end; M = [length(v)-1 for v in v_lo] # this is only for scaling, and never for iterating

    # loop Version
    for u in CartesianIndices(interpolated_array) #= #= Deprecated =# Iterators.product(IndexCartesian(), [eachindex(v) for v in v_lo]...) =#
        interpolated_array[u] =
        #= linear-part   from   edges      bilinear-part   from   vertices    =#
        (1-u[2]/M[2])*v_lo[1][u[1]] - (1-u[1]/M[1])*(1-u[2]/M[2])*v_cr[1,1] +
        (  u[2]/M[2])*v_hi[1][u[1]] - (  u[1]/M[1])*(1-u[2]/M[2])*v_cr[1,2] +
        (1-u[1]/M[1])*v_lo[2][u[2]] - (1-u[1]/M[1])*(  u[2]/M[2])*v_cr[2,1] +
        (  u[1]/M[1])*v_hi[2][u[2]] - (  u[1]/M[1])*(  u[2]/M[2])*v_cr[2,2]
    end
    return interpolated_array
end

end # module transfinite_interpolate


text:
export set!
function set! end
function (sim::Simulation)(::typeof(set!),m::AbstractArray,valueF::Function,maskF =  EVERYWHERE,gridtype::GridType = p̂)
	mask = maskF.(Coordinates(sim.grid,gridtype))
	vals = valueF.(Coordinates(sim.grid,gridtype));
	m[mask] = vals[mask];
end

function (sim::Simulation)(::typeof(set!),m::AbstractArray,vals::AbstractArray,maskF =  EVERYWHERE,gridtype::GridType = p̂)
	mask = maskF.(Coordinates(sim.grid,gridtype))
	m[mask] = vals[mask];
end


text:
module PMI

using Libdl
include("util.jl")

const __libpmi = Ref{String}("libpmi.so")
function __init__()
    if haskey(ENV, "PMI_LIBRARY")
        __libpmi[] = ENV["PMI_LIBRARY"]
    end
    if haskey(ENV, "FLUX_PMI_LIBRARY_PATH")
        __libpmi[] = ENV["FLUX_PMI_LIBRARY_PATH"]
    end
end

const PMI_SUCCESS = Cint(0) # operation completed successfully
#    PMI_FAIL (-1): operation failed
#    PMI_ERR_INIT (1): PMI not initialized
#    PMI_ERR_NOMEM (2): input buffer not large enough
#    PMI_ERR_INVALID_ARG (3): invalid argument
#    PMI_ERR_INVALID_KEY (4): invalid key argument
#    PMI_ERR_INVALID_KEY_LENGTH (5): invalid key length argument
#    PMI_ERR_INVALID_VAL (6): invalid val argument
#    PMI_ERR_INVALID_VAL_LENGTH (7): invalid val length argument
#    PMI_ERR_INVALID_LENGTH (8): invalid length argument
#    PMI_ERR_INVALID_NUM_ARGS (9): invalid number of arguments
#    PMI_ERR_INVALID_ARGS (10): invalid args argument
#    PMI_ERR_INVALID_NUM_PARSED (11): invalid num_parsed length argument
#    PMI_ERR_INVALID_KEYVALP (12): invalid keyvalp argument
#    PMI_ERR_INVALID_SIZE (13): invalid size argument

struct PMIError
    code::Cint
end
macro check(ex)
    quote
        res = $(esc(ex))
        if res != PMI_SUCCESS
            throw(PMIError(res))
        end
    end
end

"""
    init()

Initialize the PMI library for this process. Upon success, the return value
will be `true` if this process was created by `PMI_Spawn_multiple()``,
or `false` if not.
"""
function init()
    res = Ref{Cint}()
    @check @runtime_ccall((:PMI_Init, __libpmi[]), Cint, (Ptr{Cint},), res)
    res[] == 1
end

"""
    initialized()

Check if the PMI library has been initialized for this process.
"""
function initialized()
    res = Ref{Cint}()
    @check @runtime_ccall((:PMI_Initialized, __libpmi[]), Cint, (Ptr{Cint},), res)
    res[] == 1
end

"""
    finalize()

Finalize the PMI library for this process.
"""
function finalize()
    @check @runtime_ccall((:PMI_Finalize, __libpmi[]), Cint, ())
end

"""
    abort(exit_code, error_msg)

Abort the process group associated with this process.
"""
function abort(exit_code, error_msg)
    @check @runtime_ccall((:PMI_Abort, __libpmi[]), Cint, (Cint, Cstring), exit_code, error_msg)
end

"""
    get_size()

Obtain the size of the process group to which the local process belongs.
"""
function get_size()
    res = Ref{Cint}()
    @check @runtime_ccall((:PMI_Get_size, __libpmi[]), Cint, (Ptr{Cint},), res)
    res[]
end

"""
    get_rank()

Obtain the rank (0…size-1) of the local process in the process group.
"""
function get_rank()
    res = Ref{Cint}()
    @check @runtime_ccall((:PMI_Get_rank, __libpmi[]), Cint, (Ptr{Cint},), res)
    res[]
end

"""
    get_clique_size()

Obtain the number of processes on the local node.
"""
function get_clique_size()
    res = Ref{Cint}()
    @check @runtime_ccall((:PMI_Get_clique_size, __libpmi[]), Cint, (Ptr{Cint},), res)
    res[]
end

"""
    get_clique_ranks()

Get the ranks of the local processes in the process group. This is a simple
topology function to distinguish between processes that can communicate through
IPC mechanisms (e.g., shared memory) and other network mechanisms.
"""
function get_clique_ranks()
    size = get_clique_size()
    ranks = Array{Cint,1}(undef, size)
    @check @runtime_ccall((:PMI_Get_clique_ranks, __libpmi[]), Cint, (Ptr{Cint}, Cint), ranks, size)
    return ranks
end

"""
    barrier()

This function is a collective call across all processes in the process group
the local process belongs to.
"""
function barrier()
    @check @runtime_ccall((:PMI_Barrier, __libpmi[]), Cint, ())
end

# int PMI_Get_id_length_max (int *length);
function get_name_length_max()
    res = Ref{Cint}()
    @check @runtime_ccall((:PMI_KVS_Get_name_length_max, __libpmi[]), Cint, (Ptr{Cint},), res)
    res[]
end

function get_key_length_max()
    res = Ref{Cint}()
    @check @runtime_ccall((:PMI_KVS_Get_key_length_max, __libpmi[]), Cint, (Ptr{Cint},), res)
    res[]
end

function get_value_length_max()
    res = Ref{Cint}()
    @check @runtime_ccall((:PMI_KVS_Get_value_length_max, __libpmi[]), Cint, (Ptr{Cint},), res)
    res[]
end

function get_my_name()
    length = get_name_length_max()
    data = Array{UInt8, 1}(undef, length)
    @check @runtime_ccall((:PMI_KVS_Get_my_name, __libpmi[]), Cint, (Ptr{Cchar}, Cint), data, length)
    return unsafe_string(pointer(data))
end
struct KVS
    name::String
    function KVS()
        new(get_my_name())
    end
end

function put!(kvs::KVS, key, value)
    @assert length(key) <= get_key_length_max()
    @assert length(value) <= get_value_length_max()

    @check @runtime_ccall((:PMI_KVS_Put, __libpmi[]), Cint, (Cstring, Cstring, Cstring), kvs.name, key, value)
end

function get(kvs::KVS, key)
    @assert length(key) <= get_key_length_max()
    len = get_value_length_max()
    data = Array{UInt8, 1}(undef, len)
    @check @runtime_ccall((:PMI_KVS_Get, __libpmi[]), Cint, (Cstring, Cstring, Ptr{Cchar}, Cint), kvs.name, key, data, len)
    return unsafe_string(pointer(data))
end

function commit!(kvs::KVS)
    @check @runtime_ccall((:PMI_KVS_Commit, __libpmi[]), Cint, (Cstring,), kvs.name)
end

# int PMI_KVS_Create( char kvsname[], int length );
# int PMI_KVS_Destroy( const char kvsname[] );
# int PMI_KVS_Iter_first(const char kvsname[], char key[], int key_len, char val[], int val_len);
# int PMI_KVS_Iter_next(const char kvsname[], char key[], int key_len, char val[], int val_len);
# int PMI_Get_universe_size (int *size);
# int PMI_Get_appnum (int *appnum);

end # module