text:
using .Threads

@testset "Test REPL" begin
    @testset "Use Evaluator" begin
        for (raw_input, expected) in [
            (b"1\n2\n", ["1", "2"]),
            (b"let a = 1;\na;\n", ["1", "1"]),
            (b"a\n", ["ERROR: identifier not found: a"]),
            (b"1 / 0\n", ["ERROR: divide error: division by zero"]),
            (
                b"1 ++ 2\n",
                [
                    "ERROR: parser has 1 error\nERROR: parse error: no prefix parse function for PLUS found",
                ],
            ),
            (
                b"5 + 3; 23 -- ; f((\n",
                [
                    "ERROR: parser has 3 errors\nERROR: parse error: no prefix parse function for SEMICOLON found\nERROR: parse error: no prefix parse function for EOF found\nERROR: parse error: expected next token to be RPAREN, got EOF instead",
                ],
            ),
            (b"let a=fn(x) {a(x)}; a(3)\n", ["ERROR: runtime error: stack overflow"]),
            (
                b"let a = macro(x) {x + x}; a(2)\n",
                ["ERROR: macro error: we only support returning AST-nodes from macros"],
            ),
            (
                b"let reverse = macro(a, b) { quote(unquote(b) - unquote(a)); }; reverse(2 + 2, 10 - 5);\n",
                ["1"],
            ),
            (b"puts(\"Hello, world!\")\n", ["Hello, world!\nnull"]),
        ]
            input = IOBuffer(raw_input)
            output = IOBuffer(UInt8[], read = true, write = true)
            m.start_repl(input = input, output = output)
            @test String(output.data) ==
                  m.REPL_PRELUDE *
                  "\n" *
                  join(map(x -> ">> " * x, vcat(expected, [m.REPL_FAREWELL])), "\n") *
                  "\n"
        end
    end

    @testset "Use Bytecode VM" begin
        for (raw_input, expected) in [
            (b"1\n2\n", ["1", "2"]),
            (b"let a = 1;\na;\n", ["1", "1"]),
            (b"a\n", ["ERROR: identifier not found: a"]),
            (
                b"let b = 1 / 0;\nb;\n",
                ["ERROR: divide error: division by zero", "ERROR: identifier not found: b"],
            ),
            (
                b"1 ++ 2\n",
                [
                    "ERROR: parser has 1 error\nERROR: parse error: no prefix parse function for PLUS found",
                ],
            ),
            (
                b"5 + 3; 23 -- ; f((\n",
                [
                    "ERROR: parser has 3 errors\nERROR: parse error: no prefix parse function for SEMICOLON found\nERROR: parse error: no prefix parse function for EOF found\nERROR: parse error: expected next token to be RPAREN, got EOF instead",
                ],
            ),
            # (b"let a=fn(x) {a(x)}; a(3)\n", ["ERROR: stack overflow"]),
            (
                b"let a = macro(x) {x + x}; a(2)\n",
                ["ERROR: macro error: we only support returning AST-nodes from macros"],
            ),
            (
                b"let reverse = macro(a, b) { quote(unquote(b) - unquote(a)); }; reverse(2 + 2, 10 - 5);\n",
                ["1"],
            ),
            (b"puts(\"Hello, world!\")\n", ["Hello, world!\nnull"]),
        ]
            input = IOBuffer(raw_input)
            output = IOBuffer(UInt8[], read = true, write = true)
            m.start_repl(input = input, output = output, use_vm = true)
            @test String(output.data) ==
                  m.REPL_PRELUDE *
                  "\n" *
                  join(map(x -> ">> " * x, vcat(expected, [m.REPL_FAREWELL])), "\n") *
                  "\n"
        end
    end
end


text:
using SparseArrays
using JLD2
using ShaleDrillingLikelihood

modulepath = dirname(pathof(ShaleDrillingLikelihood))
filepath = joinpath(modulepath, "..", "scratch", "simulated-data.jld2")

file = jldopen(filepath, "r")
    in_j1ptr   = file["j1ptr"]
    in_j2ptr   = file["j2ptr"]
    in_tptr    = file["tptr"]
    in_jtstart = file["jtstart"]
    in_j1chars = file["j1chars"]
    in_ichars  = file["ichars"]
    in_y       = file["y"]
    in_x       = file["x"]
    in_zchars  = file["zchars"]
    in_zrng    = file["zrng"]
    in_psis    = file["psis"]
    in_ztrans  = file["ztrans"]

    θ_drill, θ_ext, θ_α0, θ_αψ, θ_θρ = file["theta"]

close(file)

ext: jl


text:
stress = zeros(n,3*6)
for i=1:length(elem[:,1])
    l =  [1 0 0;
          0 1 0;
          0 0 1;
          0.5 0.5 0;
          0 0.5 0.5;
          0.5 0 0.5]

    i1 = elem[i,:]'
    xs = nodes[i1,1]
    ys = nodes[i1,2]
    mean = zeros(3,1)
    #u = reshape( U[[2 .* i1 (2 .* i1 .- 1)]],(2,6))
    u = U[[2 .* i1 (2 .* i1 .- 1)]]'
    stress[i,1:3] = E*b(l[1,:],xs,ys)*u
    stress[i,4:6] = E*b(l[2,:],xs,ys)*u
    stress[i,7:9] = E*b(l[3,:],xs,ys)*u
    stress[i,10:12] = E*b(l[4,:],xs,ys)*u
    stress[i,13:15] = E*b(l[5,:],xs,ys)*u
    stress[i,16:18] = E*b(l[6,:],xs,ys)*u
    for k =1:6
        B = b(l[k,:],xs,ys)

        sigma = E*B*u
        mean = mean + sigma/6;
     end
     #stress[i,:] = mean
 end

max_stars_repo_name: billyon/first

max_issues_repo_name: billyon/first

max_forks_repo_name: billyon/first


text:
using PtFEM

ProjDir = dirname(@__FILE__)

# This example is identical to p.4.1.1.jl but re-uses the Cholesky
# decomposed global stiffness matrix in the 2nd part by calling
# p41(m, data) where the loaded_nodes values are doubled.

l = 1.0       # Total length [m]
N = 5         # Number of nodes
els = N - 1   # Number of finite elements
nod = 2       # Number of nodes per finite elements
nodof = 1     # Degrees of freedom for each node
np_types = 1  # Number of proerty types
EA = 1.0e5    # Strain stiffness
nip = 1       # Number of integration points

struct_el = :Rod
fin_el = :Line

data = Dict(
  # Rod(nxe, np_types, nip, fin_el(nod, nodof))
  :struc_el => getfield(Main, Symbol(struct_el))(els, np_types, nip,
    getfield(Main, Symbol(fin_el))(nod, nodof)),
  :properties => [EA;],
  :x_coords => 0.0:l/els:l,
  :support => [(N, [0])],
  :loaded_nodes => [
      (1, [-0.625]),
      (2, [-1.25]),
      (3, [-1.25]),
      (4, [-1.25]),
      (5, [-0.625])
    ]
);

data |> display
println()

@time fem, dis_df, fm_df = p41(data)
println()

# Update (double) the loaded_nodes values

F = -10.0
data[:loaded_nodes] = [(i, [0.0]) for i in 1:N]
for node in 1:N-1
  data[:loaded_nodes][node][2][1] += F / (2els)
  data[:loaded_nodes][node+1][2][1] += F / (2els)
end

# Pass in m so m.cfgsm can be re-used with this new load vector.

@time fem2, dis_df2, fm_df2 = p41(fem, data)
println()

display(dis_df2)
println()
display(fm_df2)


text:
"""Module containing auxiliary functions."""
module Functions

using ArchGDAL
using CombinedParsers
using CombinedParsers.Regexp
using DataFrames
using GeoArrays
using Rasters

include(".\\FunctionsDB.jl")

export AbstractAnalysisObject, # Supertype of structs used for the analysis
       getindex, setindex!, convert, -, # Functions overloadings
       getOrigin, getCellDims, getSidesDistances, toCoords, toIndexes, # Utility functions for retrieving informations on rasters
       compute_result!, condition, # Empty function definitions created to allow various analysis specific methods
       compute_position, expand! # Auxiliary functions for the execution of analysis

const agd = ArchGDAL

"""
    AbstractAnalysisObject

Abstract supertype of the structs defined to reppresented some kind of analysis to be performed.
A struct defined as a subtype of `AbstractAnalysisObject` should have the following fields:
- `x::Float64`: x coordinate of the current cell to be analyzed.
- `y::Float64`: y coordinate of the current cell to be analyzed.
- `direction::Int64`: Angular direction of flow (ie wind direction or water flow direction)
It may, tho it is not strictly required, also have the following fields:
- `z::Float64`: (Optional) z coordinate of the current cell to be analyzed.

The module should also provide the following functions:
- `compute_concentration!( object::YourStruct )`: where `YourStruct` is the struct subtyping `AbstractAnalysisObject`,
    the function should return the concentration of pollutant in a single cell at the coordinates given by `x`, `y`
    and possibly `z` fields.
It may also provide the following function:
- `check_result( value::Float64, numCAS::String )`:given the concentration value for a cell and the CAS number of the pollutant, returns true if said value is abnormal
    and or potentially dangerous.
"""
abstract type AbstractAnalysisObject end

@syntax cell_dims = Sequence( "Pixel Size = (", Numeric(Float64), ",", Numeric(Float64), ")" )
@syntax raster_points = Sequence( re"[^(]+", "(", re" *", Numeric(Float64), ",", re" *", Numeric(Float64), re".+" )

Base.getindex( collection::Raster{T}, x::Float64, y::Float64 ) where {T} = collection[X(Near(x)), Y(Near(y))][1]
Base.setindex!( collection::Raster{T}, v, x::Float64, y::Float64, ) where {T} = collection[X(Near(x)), Y(Near(y))] .= v
Base.convert(::Type{Int64}, n::AbstractFloat) = round(Int64, n)

ArchGDAL.getgeotransform( raster::GeoArrays.GeoArray ) = [raster.f.translation[1], raster.f.linear[1, :]..., raster.f.translation[2], raster.f.linear[2, :]...]

"""
    writeRaster( data::Array{Float32}, driver::ArchGDAL.Driver, geotransform::Vector{Float64}, refsys::AbstractString, noDataValue::Real, output_path::AbstractString )

Given a NxMxH dimensional matrix `data`, create a raster file with H NxM bands as `output_path` file, with `refsys` and `geotransfrom` as spatial references,
using `driver` to define the format.
"""
function writeRaster( data::Array{Float32}, driver::ArchGDAL.Driver, geotransform::Vector{Float64}, refsys::AbstractString, noDataValue::Real, output_file_path::AbstractString )
    rows, cols, bands = length(size(data)) < 3 ? (size(data)..., 1) : size(data)
    agd.create(output_file_path, driver=driver, width=rows, height=cols, nbands=bands, dtype=Float32) do res_raster
        for i in 1:bands
            agd.setnodatavalue!(agd.getband(res_raster, i), noDataValue)
            agd.write!(res_raster, data[:, :, i], i)
        end
        agd.setgeotransform!(res_raster, geotransform)
        agd.setproj!(res_raster, refsys)
    end
    return nothing
end

# Additional functions

"""
    getCellDims( dtm::ArchGDAL.AbstractDataset )

Return the cells' dimentions in meters for an ArchGDAL raster dataset
"""
function getCellDims( dtm::ArchGDAL.AbstractDataset )
    info = split( agd.gdalinfo( dtm ), "\n" , keepempty=false )
    pos = findfirst(occursin.("Pixel Size", info))
    size_pars = cell_dims(info[pos])
    return size_pars[2], size_pars[4]
end

"""
    getOrigin( dtm::ArchGDAL.AbstractDataset )

Return the coordinates of the origin point of the raster
"""
function getOrigin( dtm::ArchGDAL.AbstractDataset )
    info = split( agd.gdalinfo( dtm ), "\n" , keepempty=false )
    pos = findfirst(occursin.("Origin", info))
    origin_pars = raster_points(info[pos])
    return origin_pars[4], origin_pars[7]
end

"""
    getSidesDistances( dtm::ArchGDAL.AbstractDataset )

Return distances from left, upper, right and lower sides of an ArchGDAL raster dataset
"""
function getSidesDistances( dtm::ArchGDAL.AbstractDataset )
    info = split( agd.gdalinfo( dtm ), "\n" , keepempty=false )
    pos = findfirst(occursin.("Upper Left", info))
    dists_pars = raster_points.( getindex.(Ref(info), [pos, pos+3] ) )
    return dists_pars[1][4], dists_pars[1][7], dists_pars[2][4], dists_pars[2][7]
end

"""
    toCoords( dtm::ArchGDAL.AbstractDataset, r::Int64, c::Int64 )

Convert convert the indexes `r` and `c` to the coordinates of the respective cell in `dtm` raster
"""
function toCoords( dtm::ArchGDAL.AbstractDataset, r::Int64, c::Int64 )
    gtf = agd.getgeotransform(dtm)
 #  return gtf[[1,4]] .+ ( (r + 1/2) .* gtf[[2,5]] ) .+ ( (c + 1/2) .* gtf[[3,6]] )
    return gtf[[1,4]] .+ ( r .* gtf[[2,5]] ) .+ ( c .* gtf[[3,6]] )
end

function toCoords( geotransform::Vector{Float64}, r::Int64, c::Int64 )
    return geotransform[[1,4]] .+ ( (r + 1/2) .* geotransform[[2,5]] ) .+ ( (c + 1/2) .* geotransform[[3,6]] )
end

function toCoords( dtm::Rasters.Raster{T}, r::Int64, c::Int64 ) where {T}
    return dtm.dims[1][r], dtm.dims[2][c]
end

"""
    toIndexes( dtm::ArchGDAL.AbstractDataset, x::Float64, y::Float64 )

Convert coordinates to the indexes of the respective cell in `dtm` raster
"""
function toIndexes( dtm::ArchGDAL.AbstractDataset, x::Float64, y::Float64 )
    gtf = agd.getgeotransform(dtm)
    return round.( Int64, ( ( (x, y) .- gtf[[1,4]] ) ./ gtf[[2,6]] ) )
end

function toIndexes( geotransform::Vector{Float64}, x::Float64, y::Float64 )
    return round.( Int64, ( ( (x, y) .- geotransform[[1,4]] ) ./ geotransform[[2,6]] ) )
end

function compute_concentration!( object::AbstractAnalysisObject )
    throw(DomainError(object, "No function definition for `$(typeof(object))`"))
end

"""
    expand( src_r::Int64, src_c::Int64, concentration::Float64, tollerance::Int64, dem::ArchGDAL.IDataset, object::AbstractAnalysisObject )

Return a `Vecor{Tuple{Int64, Int64, Float64}}` containing all the cells of raster `dem` that have a value of concentration considerable valid according to `tollerance`
(so that the concentration of a given cell is within `tollerance` orders of magnitude from that of other cells), the returned tuples hold the indexes of the cell and the
corresponding concentration value.
The analysis strarts from the cell (`src_r`, `src_c`) assuming that the cell holds a concentration of pollutant equal to `concentration` and using `object` to account
for the specificity of the analysis.

The algorithm creates a queue using a `Vector{Int64, Int64, Float64}`, iteratively the first element of the vector, formed by the indexes of the
cell to be examined and the result of the preceding cell, is extracted and, if its result is valid according to `tollerance`, its indexes and result are added to the
result vector, otherwise a check is operated to see if the cell's concentration value is greater than its predecessor, if so it is added along with the adjacent cell's
indexes to the queue, if not the algorithm simply procedes to the following cell, the cyclecontinues untill there are no more cells in the queue.
"""
function expand( src_r::Int64, src_c::Int64, concentration::Float64, tollerance::Int64, dem::ArchGDAL.IDataset, object::AbstractAnalysisObject )
    raster::ArchGDAL.IRasterBand{Float32} = agd.getband(dem, 1)
 # Limits for rows and columns.
    max_r, max_c = size(raster)
 # Value reppresenting a lack of data for a given cell.
    noData = something(Float32(agd.getnodatavalue(raster)), -9999.0f0)
 # Vector of the cells with a valid concentration of pollutant.
    points = Tuple{Int64, Int64, Float64}[]
 # If the cell identified by (`src_r`, `src_c`) is within the raster bounds and has an actual value.
    @inbounds if 1 <= src_r <= max_r && 1 <= src_c <= max_c && raster[src_r, src_c] != noData
        valid_magnitude = -Inf
     # Coordinates of the source.
        src_coords = toCoords(dem, src_r, src_c)
     # Direction of flow angle in radians.
        dir = deg2rad(object.direction)
     # Sine and cosine of the angle.
        sindir = sin(dir)
        cosdir = cos(dir)
     # Vector of displacements between a cell and its adjacent ones.
        disp = [(1, 0), (0, 1), (-1, 0), (0, -1)]
     # Vector of cells that have already been checked and found with an inadequate result.
        visited = Tuple{Int64, Int64}[]
     # Vector of cells that still need to be checked, initially filled with the adjacents of the source.
        new_points = Tuple{Int64, Int64, Float64}[]
        for (Δr, Δc) in disp
            1 <= src_r + Δr <= max_r &&
            1 <= src_c + Δc <= max_c &&
            raster[src_r + Δr, src_c + Δc] != noData &&
            push!(new_points, (src_r + Δr, src_c + Δc, 0.0))
        end
     # Add source point to the vector of points.
        push!(points, (src_r, src_c, concentration))
     # While there are still unchecked cells.
        while !isempty(new_points)
         # Variable used to take into account the result of a cell's predecessor (used for stop condition).
            last = new_points[1][3]
         # Update position in `object`
            Δx, Δy = src_coords .- toCoords(dem, new_points[1][1], new_points[1][2])
            object.x = (Δx * cosdir) - (Δy * sindir)
            object.y = (Δx * sindir) + (Δy * cosdir)
            if :z in fieldnames(typeof(object))
               object.z = raster[new_points[1][1], new_points[1][2]]
            end
         # Compute the resulting concentration for the current cell.
            result = compute_concentration!(object)
            magnitude = round(log10(result))
         # Check if the result's value is not valid.
            if result == 0.0 || magnitude < valid_magnitude - tollerance
             # Add the cell to the vector of already visited cells.
                push!(visited, new_points[1][1:2])
             # If the result is less than that of the preceding cell skip to the next cell to be checked.
                if result > last
                    last = result
                else
                 # Remove the first element of the array.
                    deleteat!(new_points, 1)
                    continue
                end
            else
                if magnitude > valid_magnitude
                    valid_magnitude = magnitude
                end
                push!(points, (new_points[1][1], new_points[1][2], result))
            end
         # Add the adjacent cells of the current one to the cells to be checked.
            for (Δr, Δc) in disp
             # Check the new cell is within the raster bounds
                1 <= new_points[1][1] + Δr <= max_r &&
                1 <= new_points[1][2] + Δc <= max_c &&
             # Check the cell has an actual value
                raster[new_points[1][1] + Δr, new_points[1][2] + Δc] != noData &&
             # Check that it has not been encountered yet.
                all(isnothing.(findfirst.( x -> x[1] == new_points[1][1] + Δr && x[2] == new_points[1][2] + Δc, [view(points, :), view(new_points, :)] ))) &&
                (new_points[1][1] + Δr, new_points[1][2] + Δc) ∉ visited &&
             # Add it to `new_points`.
                push!(new_points, (new_points[1][1] + Δr, new_points[1][2] + Δc, last))
            end
         # Remove the first element of the array.
            deleteat!(new_points, 1)
        end
    else
      println("\n\nThe indexes of the source must correspond to a cell inside raster `dem`\n")
   end
   return points
end

function expand( src_r::Int64, src_c::Int64, concentration::Float64, tollerance::Int64, dem::ArchGDAL.AbstractDataset, area::ArchGDAL.IGeometry{ArchGDAL.wkbPolygon},
                 object::AbstractAnalysisObject )
   raster = agd.getband(dem, 1)
 # Limits for rows and columns.
   max_r, max_c = size(raster)
 # Value reppresenting a lack of data for a given cell.
   noData = agd.getnodatavalue(raster)
 # Vector of the cells with a valid concentration of pollutant.
   points = Tuple{Int64, Int64, Float64}[]
 # If the cell identified by (`src_r`, `src_c`) is within the raster bounds and has an actual value.
   @inbounds if 1 <= src_r <= max_r && 1 <= src_c <= max_c && raster[src_r, src_c] != noData
      valid_magnitude = -Inf
    # Coordinates of the source.
      src_coords = toCoords(dem, src_r, src_c)
    # Direction of flow angle in radians.
      dir = deg2rad(object.direction)
    # Sine and cosine of the angle.
      sindir = sin(dir)
      cosdir = cos(dir)
    # Vector of displacements between a cell and its adjacent ones.
      disp = [(1, 0), (0, 1), (-1, 0), (0, -1)]
    # Vector of cells that have already been checked and found with an inadequate result.
      visited = Tuple{Int64, Int64}[]
    # Vector of cells that still need to be checked, initially filled with the adjacents of the source.
      new_points = Tuple{Int64, Int64, Float64}[]
      for (Δr, Δc) in disp
         1 <= src_r + Δr <= max_r &&
         1 <= src_c + Δc <= max_c &&
         raster[src_r + Δr, src_c + Δc] != noData &&
         push!(new_points, (src_r + Δr, src_c + Δc, 0.0))
      end
    # Add source point to the vector of points.
      push!(points, (src_r, src_c, concentration))
    # While there are still unchecked cells.
      while !isempty(new_points)
       # Coordinates of the current cell.
         coords = toCoords(dem, new_points[1][1], new_points[1][2])
       # If the polygon `area` doesn't contain the current cell skip to the next one.
         if !agd.contains(area, agd.createpoint(coords...))
          # Remove the first element of the array.
            deleteat!(new_points, 1)
            continue
         end
       # Variable used to take into account the result of a cell's predecessor (used for stop condition).
         last = new_points[1][3]
       # Update position in `object`
         Δx, Δy = src_coords .- coords
         object.x = (Δx * cosdir) - (Δy * sindir)
         object.y = (Δx * sindir) + (Δy * cosdir)
         if :z in fieldnames(typeof(object))
            object.z = raster[new_points[1][1], new_points[1][2]]
         end
       # Compute the resulting concentration for the current cell.
         result = compute_concentration!(object)
         magnitude = round(log10(result))
       # Check if the result's value is not valid.
         if result == 0 || magnitude < valid_magnitude - tollerance
          # Add the cell to the vector of already visited cells.
            push!(visited, new_points[1][1:2])
          # If the result is less than that of the preceding cell skip to the next cell to be checked.
            if result > last
               last = result
            else
             # Remove the first element of the array.
               deleteat!(new_points, 1)
               continue
            end
         else
            if magnitude > valid_magnitude
               valid_magnitude = magnitude
            end
            push!(points, (new_points[1][1], new_points[1][2], result))
         end
       # Add the adjacent cells of the current one to the cells to be checked.
         for (Δr, Δc) in disp
          # Check the new cell is within the raster bounds
            1 <= new_points[1][1] + Δr <= max_r &&
            1 <= new_points[1][2] + Δc <= max_c &&
          # Check the cell has an actual value
            raster[new_points[1][1] + Δr, new_points[1][2] + Δc] != noData &&
          # Check that it has not been encountered yet.
            all(isnothing.(findfirst.( x -> x[1] == new_points[1][1] + Δr && x[2] == new_points[1][2] + Δc, [view(points), view(new_points)] ))) &&
            (new_points[1][1] + Δr, new_points[1][2] + Δc) ∉ visited &&
          # Add it to `new_points`.
            push!(new_points, (new_points[1][1] + Δr, new_points[1][2] + Δc, last))
         end
       # Remove the first element of the array.
         deleteat!(new_points, 1)
      end
   else
      println("\n\nThe indexes of the source must correspond to a cell inside raster `dem`\n")
   end
   return points
end

end # module

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
# Helper to build function
# TODO eval -> Runtime generated
function _build_ddd_function(rhs, states, parameters, iv, eval_expression::Bool = false)

    if eval_expression
        f_oop, f_iip = eval.(build_function(rhs, value.(states), value.(parameters), [value(iv)], expression = Val{true}))
    else
        f_oop, f_iip = build_function(rhs, value.(states), value.(parameters), [value(iv)], expression = Val{false})
    end

    function f(
        u::AbstractVector{T} where T,
        p::AbstractVector{T} where T,
        t::T where T
    )
        return f_oop(u, p, t)
    end

    function f(
        du::AbstractVector{T} where T,
        u::AbstractVector{T} where T,
        p::AbstractVector{T} where T,
        t::T where T
    )
        return f_iip(du, u, p, t)
    end

    function f(
        x::AbstractMatrix{T} where T,
        p::AbstractVector{T} where T,
        t::AbstractVector{T} where T
    )
        @assert size(x, 2) == length(t) "Measurements and time points must be of equal length!"

        return hcat([f(x[:,i], p, t[i]) for i in 1:size(x, 2)]...)

    end

    function f(
        y::AbstractMatrix{T} where T,
        x::AbstractMatrix{T} where T,
        p::AbstractVector{T} where T,
        t::AbstractVector{T} where T
    )
        @assert size(x, 2) == length(t) "Measurements and time points must be of equal length!"
        @assert size(x, 2) == size(y, 2) "Measurements and preallocated output must be of equal length!"

        for i = 1:size(x, 2)
            @views f(y[:, i], x[:, i], p, t[i])
        end

        return
    end

    # Dispatch on DiffEqBase.NullParameters
    f(u,p::DiffEqBase.NullParameters, t) = f(u,[], t)
    f(du, u,p::DiffEqBase.NullParameters, t) = f(du, u,[], t)

    # And on the controls
    f(u,p,t,input) = f(u,p,t)
    f(du,u,p,t,input) = f(du,u,p,t)

    return f
end

function _build_ddd_function(
    rhs,
    states,
    parameters,
    iv,
    controls,
    eval_expression::Bool = false,
)

    length(controls) < 1 &&
        return _build_ddd_function(rhs, states, parameters, iv, eval_expression)

    # Assumes zero control is zero!

    if eval_expression
        c_oop, c_iip =
            eval.(
                build_function(
                    rhs,
                    value.(states),
                    value.(parameters),
                    [value(iv)],
                    value.(controls),
                    expression = Val{true},
                ),
            )
    else
        c_oop, c_iip = build_function(
            rhs,
            value.(states),
            value.(parameters),
            [value(iv)],
            value.(controls),
            expression = Val{false},
        )
    end

    function f(
        u::AbstractVector{T} where T,
        p::AbstractVector{T} where T,
        t::T where T,
        c::AbstractVector{T} where T = zeros(eltype(u), size(controls)...),
    )
        return c_oop(u, p, t, c)
    end

    function f(
        du::AbstractVector{T} where T,
        u::AbstractVector{T} where T,
        p::AbstractVector{T} where T,
        t::T where T,
        c::AbstractVector{T} where T= zeros(eltype(u), size(controls)...),
    )
        return c_iip(du, u, p, t, c)
    end

    function f(
        x::AbstractMatrix{T} where T,
        p::AbstractVector{T} where T,
        t::AbstractVector{T} where T
    )
        @assert size(x, 2) == length(t) "Measurements and time points must be of equal length!"

        return hcat([f(x[:,i], p, t[i]) for i in 1:size(x, 2)]...)

    end

    function f(
        y::AbstractMatrix{T} where T,
        x::AbstractMatrix{T} where T,
        p::AbstractVector{T} where T,
        t::AbstractVector{T} where T
    )
        @assert size(x, 2) == length(t) "Measurements and time points must be of equal length!"
        @assert size(x, 2) == size(y, 2) "Measurements and preallocated output must be of equal length!"

        for i = 1:size(x, 2)
            @views f(y[:, i], x[:, i], p, t[i])
        end

        return
    end

    function f(
        x::AbstractMatrix{T} where T,
        p::AbstractVector{T} where T,
        t::AbstractVector{T} where T,
        u::AbstractMatrix{T} where T
    )
        @assert size(x, 2) == length(t) "Measurements and time points must be of equal length!"
        @assert size(x, 2) == size(u, 2) "Measurements and inputs must be of equal length!"

        return hcat([f(x[:,i], p, t[i], u[:, i]) for i in 1:size(x, 2)]...)

    end

    function f(
        y::AbstractMatrix{T} where T,
        x::AbstractMatrix{T} where T,
        p::AbstractVector{T} where T,
        t::AbstractVector{T} where T,
        u::AbstractMatrix{T} where T
    )
        @assert size(x, 2) == length(t) "Measurements and time points must be of equal length!"
        @assert size(x, 2) == size(y, 2) "Measurements and preallocated output must be of equal length!"
        @assert size(x, 2) == size(u, 2) "Measurements and inputs must be of equal length!"

        for i = 1:size(x, 2)
            @views f(y[:, i], x[:, i], p, t[i], u[:, i])
        end

        return
    end

    # Dispatch on DiffEqBase.NullParameters
    f(u,p::DiffEqBase.NullParameters, t) = f(u,[], t)
    f(du, u,p::DiffEqBase.NullParameters, t) = f(du, u,[], t)

    return f
end


text:
using    GSL, HypergeometricFunctions, Plots, Printf, SpecialFunctions, ..AngularMomentum, ..Basics, ..Continuum, ..Defaults, ..InteractionStrength, ..Radial, ..ManyElectron,
            ..Nuclear, ..Pulse, ..TableStrings

"""
`StrongField.computeSphericalAmplitudesUncoupled(comp::StrongField.Computation)`
    ... computes all necessary spherical amplitudes for the given initial and final levels, the Volkov states
        and polarization and by using the pulse envelope integrals. A newAmplitudes::Array{SphericalAmplitude,1} is returned.
"""
function  computeSphericalAmplitudesUncoupled(comp::StrongField.Computation)

    #Extract the initial orbital of the active electron from the many-electron comp.initialLevel and set quantum numbers
    initialOrbitals = comp.initialLevel.basis.orbitals

    #Find highest lying orbital (smallest ionization potential)
    defaultSubshell = [sh for (sh,or) in initialOrbitals][1] #This is not nice; must be a better way to simply get a default element from a Dict
    initialOrbital = initialOrbitals[defaultSubshell]
    minIonPotential = abs(initialOrbital.energy)
    for (subshell,orbital) in initialOrbitals
        if   abs(orbital.energy) < minIonPotential
            initialOrbital = orbital
            minIonPotential = abs(orbital.energy)
        end
    end

    initialEneV = convertUnits("energy: from atomic to eV", initialOrbital.energy)
    println("")
    println("Selected orbital $(initialOrbital.subshell) as active electron for the strong-field ionization process. Binding energy: $(initialEneV) eV.")
    println("")

    ls = LevelSymmetry(initialOrbital.subshell)
    n = initialOrbital.subshell.n;      l = (ls.J.num+1)/2;    initialEn = initialOrbital.energy
    if  (sign((-1)^l) == -1 && ls.parity == plus::Parity) || (sign((-1)^l) == 1 && ls.parity == minus::Parity)
        l = l - 1
    end
    l = floor(Int,l)

    if comp.settings.hydrogenic1s
        n = 1
        l = 0
    end

    lpMin = abs(l-1)
    lpMax = l+1

    # First determine which spherical SFA amplitudes need to be computed before the actual computation starts
    mNum = 1    #If comp.Average == false
    miAverageFactor = 1.      #If comp.Average == false
    if( comp.settings.mAverage )
        mNum = 2*l + 1
        miAverageFactor = 1.0/sqrt( mNum )
    end
    sfaAmplitudes = StrongField.determineSphericalAmplitudes(comp.observable,mNum,1)

    # Compute the requested SFA amplitudes
    for jm = 1:mNum
        m = -l + (jm-1)
        for  jAmp = 1:length(sfaAmplitudes[jm,1,:]) #amp in  sfaAmplitudes[jm,1]
            amp = sfaAmplitudes[jm,1,jAmp]
            thetap        = amp.theta;   phip = amp.phi;     energyp = amp.energy
            envIntegrals  = StrongField.computeEnvelopeIntegrals(comp.envelope, comp.beam, comp.polarization, thetap, phip, energyp, initialEn)
            fVolkovPlus, fVolkovMinus, fVolkovSquared = envIntegrals
            #
            reducedMEArray = ComplexF64[]
            for lp = lpMin:lpMax
                push!( reducedMEArray, pReducedMEHydrogenicUncoupled(energyp, lp, n, l, initialEn, comp.volkov) )
            end

            wminus = 0.0im;    wplus = 0.0im
            # Collect contributions from all l_p, q terms
            for  lp = lpMin:lpMax
                reducedME = reducedMEArray[lp-lpMin+1]

                KWignerEckart = 1.0/sqrt(2*lp+1) #Factor in the Wigner-Eckart theorem (CONVENTION that needs to fit the reduced ME!)
                for  q in [-1, 0, 1]
                    wplus = wplus + KWignerEckart * AngularMomentum.sphericalYlm(lp, m-q, amp.theta, amp.phi) * (-1)^q  *
                            Basics.determinePolarizationVector(q, comp.polarization, star=false)         *
                            AngularMomentum.ClebschGordan_old( AngularJ64(Rational(l)), AngularM64(Rational(m)), AngularJ64(1), AngularM64(Rational(-q)), AngularJ64(Rational(lp)), AngularM64(Rational(m-q)) ) * reducedME

                    wminus = wminus + KWignerEckart * AngularMomentum.sphericalYlm(lp, m+q, amp.theta, amp.phi)           *
                            Basics.determinePolarizationVector(q, comp.polarization, star=true)          *
                            AngularMomentum.ClebschGordan_old( AngularJ64(Rational(l)), AngularM64(Rational(m)), AngularJ64(1), AngularM64(Rational(q)), AngularJ64(Rational(lp)), AngularM64(Rational(m+q)) ) * reducedME
                end
            end

            scalarProd = scalarProdBoundContHydrogenicUncoupled(energyp, l, n, l, m, initialEn, comp.volkov)

            wa = (fVolkovPlus * wminus + fVolkovMinus * wplus) + fVolkovSquared * AngularMomentum.sphericalYlm(l, m, amp.theta, amp.phi) * scalarProd
            wa = - im/sqrt(2*pi) * wa
            wa = miAverageFactor * wa

            sfaAmplitudes[jm,1,jAmp] = SphericalAmplitude(amp.energy, amp.theta, amp.phi, wa)

            if( comp.settings.printAmplitudes )
                println(">> $(SphericalAmplitude(amp.energy, amp.theta, amp.phi, wa))")
            end
        end #for amp
    end #for jmj

    return( sfaAmplitudes )
end


text:
==(a::STaylor1, b::STaylor1) = (a.coeffs == b.coeffs)

iszero(a::STaylor1) = all(iszero, a.coeffs)

zero(::STaylor1{N,T}) where {N, T<:Number} = STaylor1(zero(T), Val(N-1))
one(::STaylor1{N,T}) where {N, T<:Number} = STaylor1(one(T), Val(N-1))

@inline +(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(a.coeffs .+ b.coeffs)
@inline -(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(a.coeffs .- b.coeffs)
@inline +(a::STaylor1) = a
@inline -(a::STaylor1) = STaylor1(minus_tuple(a.coeffs))

function +(a::STaylor1{N,T}, b::T) where {N, T<:Number}
    STaylor1{N,T}(ntuple(i -> i == 1 ? a.coeffs[1] + b : a.coeffs[i], Val(N)))
end
function +(b::T, a::STaylor1{N,T}) where {N, T<:Number}
    STaylor1{N,T}(ntuple(i -> i == 1 ? a.coeffs[1] + b : a.coeffs[i], Val(N)))
end
function -(a::STaylor1{N,T}, b::T) where {N, T<:Number}
    STaylor1{N,T}(ntuple(i -> i == 1 ? a.coeffs[1] - b : a.coeffs[i], Val(N)))
end
-(b::T, a::STaylor1{N,T}) where {N, T<:Number}  = b + (-a)

#+(a::STaylor1{N,T}, b::S) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(a,b)...)
#+(a::STaylor1{N,T}, b::STaylor1{N,S}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(a,b)...)
#+(a::STaylor1{N,T}, b::S) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(a,b)...)
#+(b::S, a::STaylor1{N,T}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(b,a)...)

#-(a::STaylor1{N,T}, b::STaylor1{N,S}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = -(promote(a,b)...)
#-(a::STaylor1{N,T}, b::S) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = -(promote(a,b)...)
#-(b::S, a::STaylor1{N,T}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = -(promote(b,a)...)

@generated function *(x::STaylor1{N,T}, y::STaylor1{N,T}) where {T<:Number,N}
     ex_calc = quote end
     append!(ex_calc.args, Any[nothing for i in 1:N])
     syms = Symbol[Symbol("c$i") for i in 1:N]
     for j = 1:N
         ex_line = :(x.coeffs[1]*y.coeffs[$j])
         for k = 2:j
             ex_line = :($ex_line + x.coeffs[$k]*y.coeffs[$(j-k+1)])
         end
         sym = syms[j]
         ex_line = :($sym = $ex_line)
         ex_calc.args[j] = ex_line
     end
     exout = :(($(syms[1]),))
     for i = 2:N
         push!(exout.args, syms[i])
     end
     return quote
                Base.@_inline_meta
                $ex_calc
                return STaylor1{N,T}($exout)
             end
end

function *(a::STaylor1{N,T}, b::T) where {N, T<:Number}
    STaylor1{N,T}(b .* a.coeffs)
end
function *(b::T, a::STaylor1{N,T}) where {N, T<:Number}
    STaylor1{N,T}(b .* a.coeffs)
end


text:
module Elpp

# Cxx stuff
using Cxx

const ElInclude = joinpath(Pkg.dir(), "Elemental/deps/usr/include/")

addHeaderDir(ElInclude, kind=C_System)
addHeaderDir("/usr/local/include/", kind=C_System)

# Libdl.dlopen(ElInclude * "../lib/libpmrrr", Libdl.RTLD_GLOBAL)
# Libdl.dlopen(ElInclude * "../lib/libElSuiteSparse", Libdl.RTLD_GLOBAL)
# Libdl.dlopen(ElInclude * "../lib/libmetis", Libdl.RTLD_GLOBAL)
# Libdl.dlopen(ElInclude * "../lib/libparmetis", Libdl.RTLD_GLOBAL)
Libdl.dlopen(ElInclude * "../lib/libEl", Libdl.RTLD_GLOBAL)

cxxinclude("El.hpp")

# Julia setup
import Base: convert, copy!, getindex, randn!, setindex!, similar, size, svdvals, svdvals!

const ElInt = icxx"sizeof(El::Int);" == 4 ? Int32 : Int64

_initialize(argc, argv) = icxx"El::Initialize($argc, $argv);"
initialize() = _initialize(Cint(0), convert(Ptr{Ptr{Cuchar}}, C_NULL))

finalize() = icxx"El::Finalize();"

### ElementalMatrix ###
# Many of the definitions in Elemental are very similar which makes it easy to define methods for a common abstract type. However, we should be careful when calling into the library with e.g. @cxx or icxx because all leaf types might not implement the actual method.

abstract type ElementalMatrix{T} <: AbstractMatrix{T} end

_getindex(A::ElementalMatrix, i::ElInt, j::ElInt) = icxx"$(A.buf).Get($i, $j);"
getindex(A::ElementalMatrix, i::Integer, j::Integer) = _getindex(A, ElInt(i - 1), ElInt(j - 1))

function size(A::ElementalMatrix, i::ElInt)
    if i < 1
        throw(error(""))
    elseif i == 1
        return Int(icxx"$(A.buf).Height();")
    elseif i == 2
        return Int(icxx"$(A.buf).Width();")
    else
        return 1
    end
end
size(A::ElementalMatrix) = (size(A, 1), size(A, 2))

similar(A::ElementalMatrix) = typeof(A)(size(A)...)
similar(A::ElementalMatrix, sz::Tuple{Int,Int}) = typeof(A)(ElInt(sz[1]), ElInt(sz[2]))

function copy!{T<:ElementalMatrix}(dest::T, src::T)
    @cxx El::Copy(src.buf, dest.buf)
    return dest
end

function randn!(A::ElementalMatrix)
    _randn!(A, size(A)...)
    return A
end

svdvals(A::ElementalMatrix) = svdvals!(copy(A))

### AbstractDistMatrix ###
abstract type AbstractDistMatrix{T} <: ElementalMatrix{T} end

_resize!(A::AbstractDistMatrix, i::ElInt, j::ElInt) = icxx"$(A.buf).Resize($i, $j);"
resize!(A::AbstractDistMatrix, i::Integer, j::Integer) = _resize!(A, ElInt(i), ElInt(j))

_reserve(A::AbstractDistMatrix, n::ElInt) = icxx"$(A.buf).Reserve($n);"
reserve(A::AbstractDistMatrix, n::Integer) = _reserve(A, ElInt(n))

_queueUpdate(A::AbstractDistMatrix{T}, i::ElInt, j::ElInt, x::T) where {T} = icxx"$(A.buf).QueueUpdate($i, $j, $x);"
queueUpdate(A::AbstractDistMatrix{T}, i::Integer, j::Integer, x) where {T} = _queueUpdate(A, ElInt(i - 1), ElInt(j - 1), T(x))

processQueues!(A::AbstractDistMatrix) = icxx"$(A.buf).ProcessQueues();"

_zeros!(A::AbstractDistMatrix, m::ElInt, n::ElInt) = icxx"El::Zeros($(A.buf), $m, $n);"
zeros!(A::AbstractDistMatrix, m::Integer = size(A,1), n::Integer = size(A,2)) = _zeros!(A, ElInt(m), ElInt(n))

### Matrix ###

struct Matrix{T} <: ElementalMatrix{T}
    # buf::Cxx.CppValue{Cxx.CxxQualType{Cxx.CppTemplate{Cxx.CppBaseType{symbol("El::Matrix")},Tuple{T}},(false,false,false)},56}
    buf::Any
end

convert(::Type{Matrix{Float32}}, m::ElInt = 0, n::ElInt = 0) = Matrix{Float32}(icxx"El::Matrix<float>($m,$n);")
convert(::Type{Matrix{Float64}}, m::ElInt = 0, n::ElInt = 0) = Matrix{Float64}(icxx"El::Matrix<double>($m,$n);")

_randn!(A::Matrix, i::ElInt, j::ElInt) = icxx"Gaussian($(A.buf), $i, $j);"

function LinearAlgebra.svdvals!(A::Matrix{T}) where {T}
    s = Matrix{T}(min(size(A)...),1)
    @cxx El::SVD(A.buf, s.buf)
    return s
end

# immutable Grid

# DistMatrix

struct DistMatrix{T,U,V} <: AbstractDistMatrix{T}
    # buf::Cxx.CppValue{Cxx.CxxQualType{Cxx.CppTemplate{Cxx.CppBaseType{symbol("El::DistMatrix")},Tuple{T,U,V,U}},(false,false,false)},144}
    buf::Any
end

const MC = Cxx.CppEnum{symbol("El::DistNS::Dist")}(0)
const MR = Cxx.CppEnum{symbol("El::DistNS::Dist")}(2)

# These converstions for C++ enums are necessary
convert(::Type{Cxx.CppEnum{symbol("El::DistNS::Dist")}}, x::Int) = Cxx.CppEnum{symbol("El::DistNS::Dist")}(x)
convert(::Type{Cxx.CppEnum{symbol("El::DistWrapNS::DistWrap")}}, x::Int) = Cxx.CppEnum{symbol("El::DistWrapNS::DistWrap")}(x)

convert(::Type{DistMatrix{Float32,MC,MR}}, m::ElInt = 0, n::ElInt = 0) = DistMatrix{Float32,MC,MR}(icxx"El::DistMatrix<float,El::MC,El::MR>($m,$n);")
convert(::Type{DistMatrix{Float64,MC,MR}}, m::ElInt = 0, n::ElInt = 0) = DistMatrix{Float64,MC,MR}(icxx"El::DistMatrix<double,El::MC,El::MR>($m,$n);")
convert(::Type{DistMatrix{T}}, m::ElInt = 0, n::ElInt = 0) where {T} = DistMatrix{T,MC,MR}(m, n)

_randn!(A::DistMatrix, i::ElInt, j::ElInt) = icxx"Gaussian($(A.buf), $i, $j);"

function LinearAlgebra.svdvals!(A::DistMatrix{T}) where {T}
    s = DistMatrix{T}(min(size(A)...),1)
    @cxx El::SVD(A.buf, s.buf)
    return s
end

# DistSparseMatrix
struct DistSparseMatrix{T} <: AbstractDistMatrix{T}
    buf::Any
end

convert(::Type{DistSparseMatrix{Float32}}, m::ElInt = 0, n::ElInt = 0) = DistSparseMatrix{Float32}(icxx"El::DistSparseMatrix<float>($m,$n);")
convert(::Type{DistSparseMatrix{Float64}}, m::ElInt = 0, n::ElInt = 0) = DistSparseMatrix{Float64}(icxx"El::DistSparseMatrix<double>($m,$n);")

processLocalQueues!(A::DistSparseMatrix) = icxx"$(A.buf).ProcessLocalQueues();"

# DistMultiVec
struct DistMultiVec{T} <: AbstractDistMatrix{T}
    buf::Any
end

convert(::Type{DistMultiVec{ElInt}}, m::ElInt = 0, n::ElInt = 0) = DistMultiVec{ElInt}(icxx"El::DistMultiVec<El::Int>($m,$n);")
convert(::Type{DistMultiVec{Float32}}, m::ElInt = 0, n::ElInt = 0) = DistMultiVec{Float32}(icxx"El::DistMultiVec<float>($m,$n);")
convert(::Type{DistMultiVec{Float64}}, m::ElInt = 0, n::ElInt = 0) = DistMultiVec{Float64}(icxx"El::DistMultiVec<double>($m,$n);")

_setindex!(A::DistMultiVec{T}, x::T, i::ElInt, j::ElInt) where {T} = icxx"$(A.buf).Set($i, $j, $x);"
setindex!(A::DistMultiVec{T}, x, i::Integer, j::Integer) where {T} = _setindex!(A, T(x), ElInt(i - 1), ElInt(j - 1))

### SOCP ###
socp(A::DistSparseMatrix, G::DistSparseMatrix, b::DistMultiVec, c::DistMultiVec, h::DistMultiVec, orders::DistMultiVec{ElInt}, firstInds::DistMultiVec{ElInt}, x::DistMultiVec, y::DistMultiVec, z::DistMultiVec, s::DistMultiVec) = icxx"El::SOCP($(A.buf), $(G.buf), $(b.buf), $(c.buf), $(h.buf), $(orders.buf), $(firstInds.buf), $(x.buf), $(y.buf), $(z.buf), $(s.buf));"

# DArray stuff

using DistributedArrays
function toback(A::DArray{Float64,2})
    rs = Array{Any}(size(A.chunks))
    for p in eachindex(A.chunks)
        ind = A.indexes[p]
        rs[p] = remotecall(A.pids[p], () -> begin
            lA = localpart(A)
            AlA = Elpp.DistMatrix{Float64}()
            Elpp.zeros!(AlA, size(A)...)
            for j = 1:size(lA, 2)
                for i = 1:size(lA, 1)
                    Elpp.queueUpdate(AlA, start(ind[1]) + i - 1, start(ind[2]) + j - 1, lA[i,j])
                end
            end
            Elpp.processQueues!(Main.AlA)
            AlA
        end)
    end
    rs
end

end


text:
module StatisticalVerificationOfIsingEnergies

using Distributions
using Cumulants
using StatsBase

export estimate_ground_state_energy, squared_error, bootstrap_get_pvalue,
       bootstrap_hists_of_mins, estiamte_temperature

"""
    estimate_ground_state_energy(energies::Vector{Float64}, α::Float64)

Returns Float64, estimated ground state energy.

```jldoctest

julia> energies = [0.1, 1., 1.2, 2., 2.2, 3., 3.1, 3.5];

julia> estimate_ground_state_energy(energies, 0.19)
9.745023110079497

```
"""
function estimate_ground_state_energy(energies::Vector{Float64}, α::Float64)
    μ = mean(energies)
    σ = std(energies)
    η = skewness(energies)

    μ-(α+2)/(α+1)*σ/η
end

"""
    squared_error(energies::Vector{Float64}, α::Float64)

Returns Float64, estimated standard deviation of the estimation of the ground state energy

```jldoctest

julia> energies = [1. ,1.5, 2., 3., 4.];

julia> squared_error(energies, 0.19)
8.633622351519547
```
"""
function squared_error(energies::Vector{Float64}, α::Float64)
    N = length(energies)
    B = (α+2)/(α+1)
    c = [Array(x)[1] for x in cumulants(reshape(energies, :, 1), 6, 1)]
    δc2 = sqrt((c[4] + 2*c[2]^2)/N)
    δc3 = sqrt((c[6] + 9*c[2]*c[4]+ 9*c[3]^2 + 6*c[2]^3)/N)
    e2 = B*2*(c[2]/c[3])*δc2
    e3 = B*(c[2]/c[3])^2*δc3
    (e2^2+e3^2)^0.5
end

"""
    bootstrap_hists_of_mins(energies::Vector{Float64}, α::Float64, s::Int, l::Int=length(x))

Returns s-long vector of Float64, estimated ground state energies from l-long samples of dootstrap resampling

julia> energies = [1. ,1.5, 2., 3., 4.];

julia> Random.seed!(1234);

```jldoctest
julia> h = bootstrap_hists_of_mins(energies, 0.19, 3)
3-element Vector{Float64}:
 10.327274688631618
 -4.407048476374387
 -4.721357273543485
```
"""
function bootstrap_hists_of_mins(energies::Vector{Float64}, α::Float64, s::Int, l::Int=length(energies))
    ret = zeros(s)
    Threads.@threads for i ∈ 1:s
        ret[i] = estimate_ground_state_energy(rand(energies, l), α)
    end
    ret
end

"""
    bootstrap_get_pvalue(energies::Vector{Float64}, α::Float64, s::Int=1000)

returns Float64, p-value that energies contains the ground state if the Ising sampler,
s is the number ob bootstrap samples.

```jldoctest
julia> Random.seed!(1234);

julia> energies = [0.1, 1., 1.2, 2., 2.2, 3., 3.1, 3.5];

julia> bootstrap_get_pvalue(energies, 0.19)
0.655
```
"""
function bootstrap_get_pvalue(energies::Vector{Float64}, α::Float64, s::Int=1000)
    y = bootstrap_hists_of_mins(energies, α, s)
    f = ecdf(y)
    Hmin = minimum(energies)
    1.0 - f(Hmin)
end

"""
    estiamte_temperature(energies::Vector{Float64}, ground_energy::Float64)

returns Float64, estimated temperature.

```jldoctest
julia> energies = [1. ,1.5, 2., 3., 4.];

julia> estiamte_temperature(energies, 0.)
0.8958957699200034
```
"""
function estiamte_temperature(energies::Vector{Float64}, ground_energy::Float64)
    EH = ground_energy - mean(energies)
    EH/(std(energies)*skewness(energies)*EH-var(energies))
end

end


text:
uniform(dim, n) = rand(dim, n)

true_uniform(n, m) = reduce(hcat, collect.(Iterators.product((1:n) / n .- 1 / 2n, (1:m) / m .- 1 / 2m)))

function net(nline, npoint)
    dx = 1 / npoint
    xs = range(dx / 2, length = npoint, step = dx)
    dy = 1 / nline
    ys = range(dy / 2, length = nline, step = dy)

    a = [[x, y] for x in xs, y in ys]
    b = [[y, x] for x in xs, y in ys]
    reduce(hcat, [a[:]; b[:]])
end

function circle_cluster(n)
    ws = range(0., length=n, step=2π / n)
    [f(w) for f in (cos, sin), w in ws]
end

normal_cluster(dim, n) = randn(dim, n)

function cross_cluster(w)
    T = [
        cos(w) sin(w)
        -sin(w) cos(w)
    ]
    x = [
        -2 -1 0 1 2 0 0 0 0
        0 0 0 0 0 -2 -1 1 2
    ]
    mapreduce(x -> T * x, hcat, eachcol(x))
end

clusters(centers, cluster) = mod.(reduce(hcat, [cluster() .+ c for c in eachcol(centers)]), 1.)

function voronoi_metric(x, support)
    a = colwise(PeriodicEuclidean(ones(2)), support, x) |> sort
    abs(a[2] - a[1])
end

function voronoi(n, support)
    ndims = size(support, 1)
    all_points = rand(ndims, 1_000_000)

    a = [voronoi_metric(x, support) for x in eachcol(all_points)]

    ix = sortperm(a)
    points = all_points[:, ix[1:n]]

    # scatter(points[1, :], points[2, :], markersize=2)
    # scatter!(support[1, :], support[2, :]) |>  display

    points
end


text:
"""
    p = dscut(p1, y1, p2, y2)

return the point `p` at which y = 0, assuming y varies linearly
from `y1` at `p1` to `y2` at `p2`.
"""
function dscut(p1, y1, p2, y2)
    k1 = y2 / (y2 - y1)
    k2 = 1 - k1
    n = size(p1, 1)
    return k1 .* p1 + k2 .* p2
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
#Functions to find an input point that predicts a desired output
import ProgressMeter
using Optim

function search{D}(target::Real, model::GaussianProcess{D})
    moveinfo = Array(Vector{Float64}, length(model))
    ProgressMeter.@showprogress for i=1:length(model)
        moveinfo[i] = calcmove(i, model, target)
    end

    bidx = indmin([m[1] for m in moveinfo])
    if isinf(moveinfo[bidx][1])
        return NaN*ones(D)
    end
    adj_idx = round(Int, moveinfo[bidx][2])
    adj_val = moveinfo[bidx][3]

    npt = adjustparam(getx(model,bidx), adj_idx, adj_val)
    npt
end

function searchgreedy{N}(target, tol, gp::GaussianProcess{N}, startx)
    Cd = 1.0
    count = 0
    function optf(pt, grad)
        count += 1
        length(grad) == 0 || error("Request for grad")
        mpt, sig2pt = query(pt, gp)
        Cprob = 1 - normprob(target-tol, target+tol, mpt, sig2pt)
        Cdist = Cd*abs(mpt - target)
        Cprob + Cdist
    end
    opt_pts = zeros(N, size(startx,2))
    opt_vals = zeros(size(startx,2))
    for i=1:size(startx,2)
        #Global search
        optprob = NLopt.Opt(:GN_CRS2_LM, N)
        #optprob = NLopt.Opt(:LN_COBYLA, N)
        NLopt.lower_bounds!(optprob, -1.0*ones(N))
        NLopt.upper_bounds!(optprob,  1.0*ones(N))
        NLopt.min_objective!(optprob, optf)
        NLopt.ftol_abs!(optprob, 1e-8)
        NLopt.maxeval!(optprob, 500)
        minx = zeros(N)
        minf = 0.0
        ret = 0
        try
            (minf, minx, ret) = NLopt.optimize(optprob, startx[:,i])
        catch ex
            println("startx: ", startx)
            rethrow()
        end

        #Local refinement
        optprob = NLopt.Opt(:LN_COBYLA, N)
        NLopt.lower_bounds!(optprob, -1.0*ones(N))
        NLopt.upper_bounds!(optprob,  1.0*ones(N))
        NLopt.min_objective!(optprob, optf)
        NLopt.ftol_abs!(optprob, 1e-8)
        NLopt.maxeval!(optprob, 500)
        try
            (minf, minx, ret) = NLopt.optimize(optprob, minx)
        catch ex
            println("minx: ", minx)
            rethrow()
        end

        opt_pts[:,i] = minx
        opt_vals[i] = minf
    end
    opt_pts, opt_vals
end

function searchgrad{D}(target, model::GaussianProcess{D}, startx)
    count = 0
    function optf(pt, grad)
        count += 1
        mpt, sig2pt = query(pt, model)
        if length(grad) > 0
            gradient!(grad, pt, model)
            grad[:] = 2.0*(mpt - target)*grad
        end
        #println("#$count")
        #println("  x: $pt")
        #println("  g: $grad")
        (mpt - target)^2
    end

    function optg!(pt::Vector, grad)
        mpt, sig2pt = query(pt, model)
        gradient!(grad, pt, model)
        grad[:] = 2.0*(mpt - target)*grad
    end

    opt_pts = zeros(D, size(startx,2))
    opt_vals = zeros(size(startx,2))
    for i=1:size(startx,2)
        optprob = NLopt.Opt(:LD_LBFGS, D)
        NLopt.lower_bounds!(optprob, -1.0*ones(D))
        NLopt.upper_bounds!(optprob,  1.0*ones(D))
        NLopt.min_objective!(optprob, optf)
        NLopt.ftol_abs!(optprob, 1e-8)
        NLopt.maxeval!(optprob, 500)
        try
            (minf, minx, ret) = NLopt.optimize(optprob, startx[:,i])
            #res = optimize(optf, optg!, startx[:,i];
                           #method = LBFGS(),
                           #iterations = 100,
                           #)
            #minx = Optim.minimizer(res)
            #minf = Optim.minimum(res)
            opt_pts[:,i] = minx
            opt_vals[i] = minf
        catch ex
            println("startx: ", startx[:,i])
            rethrow()
        end
    end
    opt_pts, opt_vals
end

function searchline{D}(
    target, tol, model::GaussianProcess{D}, func;
    maxDist=0.5, dT=0.0005
    )
    params = [-D:-1; 1:D]
    param_out = zeros(length(params), length(model))
    prerr_out = Inf*ones(length(params), length(model))
    trerr_out = Inf*ones(length(params), length(model))

    println("GPR: Searching over $(length(model)) points")
    ProgressMeter.@showprogress for midx=1:length(model)
        #println("Base point: ", midx)
        xbase = getx(model,midx)
        for (pidx,p) in enumerate(params)
            #print("\tParam ",p," : ")
            dx = zeros(D)
            dx[abs(p)] = sign(p)
            for t=0:dT:maxDist
                xtest = xbase + t*dx
                val,err = query(xtest, model)
                if abs(val - target) < tol
                    terr = abs(val - func(xtest))

                    param_out[pidx, midx] = t
                    prerr_out[pidx, midx] = err
                    trerr_out[pidx, midx] = terr
                    #println(round(Int,t/dT)+1)
                    break
                end
            end
            #Check if no solution was found for the target
            if trerr_out[pidx, midx] < 0.0
                #println(-1)
            end
        end
    end
    param_out, prerr_out, trerr_out
end

function calcmove{D}(idx, model::GaussianProcess{D}, target::Real)
    opt_repeat_param_tol = Main.options["repeat_param_tol"]

    params = [-D:-1; 1:D]
    param_adjs = zeros(2D,3)
    for pidx=1:length(params)
        oneparam!(idx, pidx, model, target, params, param_adjs)
    end
    bidx = indmin(param_adjs[:,1])

    error_val, dx, expct = param_adjs[bidx, :]
    adjinfo = [error_val, params[bidx], dx, expct]

    if isfinite(error_val)
        xval = getx(model,idx)
        nx = adjustparam(xval, params[bidx], dx)

        rpt = isrepeat(nx, model, opt_repeat_param_tol)
        if rpt > 0
            adjinfo[1] = Inf
        end
    end
    adjinfo
end

#Calculate the move for one param at one base point
function oneparam!{D}(
    idx::Integer,
    pidx::Integer,
    model::GaussianProcess{D},
    target::Real,
    params,
    param_adjs
    )
    opt_search_dx    = Main.options["gpr_search_dx"]
    opt_search_maxdx = Main.options["gpr_search_maxdx"]
    opt_search_tol   = Main.options["gpr_search_tol"]

    xval = getx(model,idx)
    fval = getf(model,idx)
    curr_err = fval - target

    param_adjs[pidx,:] = [Inf 0.0 0.0]
    for t=0:opt_search_dx:opt_search_maxdx
        xtest = adjustparam(xval, params[pidx], t)
        val,err = query(xtest, model)
        if abs(val - target) < opt_search_tol
            param_adjs[pidx,:] = [err t val]
            break
        end
    end
end


text:
using Plots

struct Source_Point
    posX             #array [xpos, ypos]
    posY
    Λ               #Magnitude. Lambda
end
struct Sink_Point
    posX
    posY
    Λ
end
struct Freestream
    Λ
    direction       #to the right
end
struct Vortex_Point
    posX             #[posx,posy]
    posY
    Λ
    direction
end
struct Point_Matrix         #This is to make the points we are testing.
    posX
    posY
    velX
    velY
end

function distance_finder(x1, y1, x2, y2)
    r = sqrt((x1-x2)^2 + (y1-y2)^2);

    return r;
end
function main()
    xVals = 0;   #X Values
    yVals = 0;   #Y Values

    # Points = [xVals, yVals]

    # scatter(xVals, yVals);

    source1 = Source_Point(1,0,2);  #located at (1,1) and has strength  of 2.
    sink1 = Sink_Point(9,0,1);       #Sink located at (9,1), stength of 1.
    vortex1 = Vortex_Point(5,0,1,1);   #Vortex located at (2,2), stength of 1, direction out of chart.

    # source_Plot = scatter([source1.posX, source1.posY], label = "Source")
    # sink_Plot = scatter!([sink1.posX, sink1.posY], label = "Sink")
    # vortex_Plot = scatter!([vortex1.posX, vortex1.posY], label = "Vortex")

    evaluation_Points  = Point_Matrix[][];

    #This is where we need to bulid a mesh around the airfoil
    for i = 0 to 20
        for j = 6 to 20
            evaluation_Points[i][j - 6] = Point_Matrix(i, j, 0, 0);  #constructor
        end
    end
    for i in evaluation_Points
        for j in evaluation_Points[i]
            evalutation_Points[i][(-j)+14] = Point_Matrix(i, -j, 0, 0);
        end
    end

    plotname = plot();
    p1 = display(scatter!([source1.posX], [source1.posY], label = "Source", markersize = 8, markercolor = [:red]))
    p2 = display(scatter!([sink1.posX], [sink1.posY], label = "Sink", markersize = 8, markercolor = [:black]))
    p3 = display(scatter!([vortex1.posX], [vortex1.posY], label = "Vortex", markersize = 8, markercolor = [:blue]))
    for i = 0 to 20
        for j = 0 to 40
            p4 = display(scatter!(evaluation_Points[i][j].posX, evaluation_Points[i][j].posY))
        end
    end
    xlims!(0, 10);
    ylims!(-10, 10);

    # savefig(plotname, "plots.pdf")

    #Velocity of each point around airfoil has an x and y component.

    #Sources - Influenced VELOCITY
    # V_vector = Λ / (2*π*r)
    # Vx = (Λ)cos(θ) / (2*π*r)        (1)
    # Vy = (Λ)sin(θ) / (2*π*r)         (2)

    #Vortex = Influenced Velocity
    # V_vector = -(Γ) / 2*π*r
    # Vx = Γ*sin(θ) / 2*π*r     #Components
    # Vy = -Γ*cos(θ) / 2*π*r

    #This is what I understand:
    #   We have a lattice of points around and airfoil representing the air passing over the airfoil.
    #   At each point, the air speed/direction is affected by the Vortex, source, sink, and Freestream.
    #The Flow Field adds a straight up velocity to the right.
    #The Source and sink act as a force that pulls/pushes on the point. Like the vortex their influence on the point, or their induced velocity, is expressed in Equations (1) and (2). r is the distance between the point and the anomoly.

    #Therefore we need to:
    #   1. Determine the location of the source, sink, and vortex.
    #   2. Determine the magnitude of the source, sink, vortex, and filed flow
    #   3. Define area of Airfoil based upon the locations of the anomolies.
    #   4. Create a lattice of points around airfoil.
    #   5. Determine resultant velocity at each point based on their relative position to each anomoly.
    #   6. Plot the results as a quiver().
    #   6b. Download the .dat file of used airfoil, and superimpose of plot.
    #   7. Determine Pressure at each point (eq 2.102) SUPER EASY. Try Making a Color Plot, if you can.

end

main();


text:
using ToyRockets

payload = Payload(22.8e3)

rocket = Rocket(payload, SmallTank(), Rutherford())

println("Mass of simple rocket: ", mass(rocket))
println("Thrust of simple rocket: ", mass(rocket))


text:
using PyPlot

function traffic(uinit="step", method="godunov", m=400; T=2, kmul=0.8)
    # Discretization
    h = 4 / m
    x = ((collect(1:m).-0.5) ./(m/4)).-2.0

    k = kmul*h
    N = ceil(Int, T/k)
    a = findmin(map(x->abs(x+(h/2)),x))[2]

    # Initial condition, double step function
    if uinit == "step"
        r = @. Float64(x<0.0)*0.8
	r = zeros(size(x)[1])
        r[1] = 0.5
        "u = zeros(size(x)[1])"
        "u = @. 2 * Float64(abs(x-0.5)<0.25) - 0.5"
    else
        error("Unknown initial condition type")
    end

    clf(); axis([-2, 2, -1, 1.5]); grid(true); ph, = plot(x,r)  # Setup plotting
    #N=5
    for it = 1:N
        rr = r
        rl =  r[[end; 1:end-1]] # Periodic boundary conditions

        F = zeros(m)
        if method == "godunov"
           for i = 1:m
                #if i<=4
                #   println("i=",i," rl[i]=",rl[i]," rr[i]=",rr[i])
                #end
		if rl[i] > rr[i]
		   F[i] = max(rl[i]*(1-rl[i]), rr[i]*(1-rr[i]))
		elseif rl[i] < rr[i]
		   F[i] = min(rl[i]*(1-rl[i]), rr[i]*(1-rr[i]))
		else
		 F[i] = 0
		end

           end
        elseif method == "upwind"
           for i = 1:m
               if (ur[i] + ul[i])/2 > 0
                    F[i] = ul[i]^2/2
               else
                    F[i] = ur[i]^2/2
               end
           end
        elseif method == "roe"
           for i = 1:m
              F[i] = 0.5*((rl[i]*(1-rl[i])) + (rr[i]*(1-rr[i]))) - 0.5*(abs(1.01-rl[i]-rr[i])*(rr[i] - rl[i]))
           end
        else
           error("Unknown method type")
        end

        if div(it,125)%2 == 0
            F[a] = 0
        end
	if ceil(Int,div(it,237.5))%2 == 0
            F[216] = 0
        end
	println(F[a])
        @. r -= k/h * (F[[2:end; 1]] - F)
        #println("it=",it," k=",k," h=",h,"  r=",r[1:4]," F=",F[1:4])

        r[1] = 0.5
        ph.set_data(x,r); pause(1e-3) # Update plot
    end
    return r
end


text:
abstract Class{Name}

function call{T <: Class}(::Type{T}, brick::Brick, path::AbstractString)
    attr = T()
    for name in fieldnames(T)
        field_type = fieldtype(T, name)
        setfield!(attr, name, field_type(brick, path, name))
    end
    attr
end

type TachoMotorClass <: Class{:tacho_motor}
    address::ReadOnly{ASCIIString}
    commands::ReadOnly{Vector{ASCIIString}}
    driver_name::ReadOnly{ASCIIString}
    command::WriteOnly{ASCIIString}
    count_per_rot::ReadOnly{Int}
    duty_cycle::ReadOnly{Int}
    duty_cycle_sp::ReadWrite{Int}
    speed_sp::ReadWrite{Int}
    speed_regulation::ReadWrite{ASCIIString}
    position_sp::ReadWrite{Int}
    position::ReadWrite{Int}
    polarity::ReadWrite{ASCIIString}
    stop_command::ReadWrite{ASCIIString}
    stop_commands::ReadOnly{Vector{ASCIIString}}

    TachoMotorClass() = new()
end

type LegoSensorClass <: Class{:lego_sensor}
    address::ReadOnly{ASCIIString}
    commands::ReadOnly{Vector{ASCIIString}}
    driver_name::ReadOnly{ASCIIString}
    command::WriteOnly{ASCIIString}
    value0::ReadOnly{Int}
    value1::ReadOnly{Int}
    value2::ReadOnly{Int}
    value3::ReadOnly{Int}
    value4::ReadOnly{Int}
    value5::ReadOnly{Int}
    value6::ReadOnly{Int}
    value7::ReadOnly{Int}
    num_values::ReadOnly{Int}
    decimals::ReadOnly{Int}
    modes::ReadOnly{Vector{ASCIIString}}
    mode::ReadWrite{ASCIIString}
    bin_data::ReadOnly{ASCIIString}
    bin_data_format::ReadOnly{ASCIIString}
    poll_ms::ReadWrite{Int}

    LegoSensorClass() = new()
end

@generated name{C <: Class}(::Type{C}) = :("$(C.super.parameters[1])")
@generated class_path{C <: Class}(::Type{C}) = :("$(joinpath("sys/class", replace(string(name(C)), '_', '-')))")


text:
function build_robot!(vis::Visualizer, model::WallPendulum12; name::Symbol=:WallPendulum12,
		r = 0.01, l = 1.0, n_leg::Int = 30, α = 1.0)

	l = model.l
	d = model.d

	r = convert(Float32, r)
	r_contact = convert(Float32, r * 2)
	# r_arm = convert(Float32, r * 0.3)
	wall_thickness = 0.15
	# l_joint = 0.16
	body_mat = MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, α))
	contact_mat = MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0.0, α))
	wall_mat = MeshPhongMaterial(color = RGBA(0.7, 0.7, 0.7, 1.0))

	default_background!(vis)

	link = Cylinder(Point3f0(0.0), Point3f0(0.0, 0.0, l), r)
	setobject!(vis[name][:robot]["link"], link, body_mat)
	setobject!(vis[name][:robot]["endeffector"], Sphere(Point3f0(0.0), r_contact), contact_mat)
	setobject!(vis[name][:robot]["linkbase"], Sphere(Point3f0(0.0), r), body_mat)

	# setobject!(vis[name][:robot]["contact1"], Sphere(Point3f0(0.0), r_contact), contact_mat)
	#
	# for i = 1:n_leg
	# 	setobject!(vis[name][:robot]["arm$i"], Sphere(Point3f0(0), r_arm), body_mat)
	# end

	setobject!(vis[name][:env]["fixed_wall1"],
		Rect(Vec(0, 0, 0),Vec(wall_thickness, wall_thickness, 1.5)), wall_mat)
	setobject!(vis[name][:env]["fixed_wall2"],
		Rect(Vec(0, 0, 0),Vec(wall_thickness, wall_thickness, 1.5)), wall_mat)

	setobject!(vis[name][:env][:spring_bundle1]["spring_wall1"],
		Rect(Vec(0, 0, 0),Vec(wall_thickness/10, wall_thickness, 0.25)), wall_mat)
	setobject!(vis[name][:env][:spring_bundle2]["spring_wall2"],
		Rect(Vec(0, 0, 0),Vec(wall_thickness/10, wall_thickness, 0.25)), wall_mat)

	spring = Cylinder(Point3f0(0.0), Point3f0(0.0, 0.0, wall_thickness), r)
	setobject!(vis[name][:env][:spring_bundle1]["spring1"], spring, body_mat)
	setobject!(vis[name][:env][:spring_bundle2]["spring2"], spring, body_mat)

	p_spring1 = [                     0*l*sin(d) + r_contact; -wall_thickness/2; l-0.125]
	p_spring2 = [-wall_thickness/10 - 0*l*sin(d) - r_contact; -wall_thickness/2; l-0.125]
	settransform!(vis[name][:env]["fixed_wall1"], Translation([                  0.2 + 1.5 * r; -wall_thickness/2; 0.0]))
	settransform!(vis[name][:env]["fixed_wall2"], Translation([-wall_thickness - 0.2 - 1.5 * r; -wall_thickness/2; 0.0]))
	settransform!(vis[name][:env][:spring_bundle1]["spring_wall1"], Translation(p_spring1))
	settransform!(vis[name][:env][:spring_bundle2]["spring_wall2"], Translation(p_spring2))
	settransform!(vis[name][:env][:spring_bundle1]["spring1"], compose(Translation(p_spring1 + [wall_thickness/10, wall_thickness/2, 0.125]), LinearMap(RotY(π/2))))
	settransform!(vis[name][:env][:spring_bundle2]["spring2"], compose(Translation(p_spring2 + [0,                 wall_thickness/2, 0.125]), LinearMap(RotY(-π/2))))

	return nothing
end

function set_robot!(vis::Visualizer, model::WallPendulum12, x::AbstractVector;
		name::Symbol=:WallPendulum12, r = 0.01, n_leg::Int=30)

	l = model.l
	d = model.d

	r = convert(Float32, r)
	r_contact = convert(Float32, r * 2)

	p = [kinematics(x)[1], 0.0, kinematics(x)[2]]

	# settransform!(vis[name][:robot]["joint1"], compose(cable_transform(zeros(3), pee), Translation(0,0,0)))
	# settransform!(vis[name][:robot]["joint2"], cable_transform(zeros(3), pee))
	settransform!(vis[name][:robot]["link"], cable_transform(p, zeros(3)))
	settransform!(vis[name][:robot]["endeffector"], Translation(p))
	if p[1] < -l * sin(d) # left
		settransform!(vis[name][:env][:spring_bundle2], Translation(p[1], 0, 0))
	else
		settransform!(vis[name][:env][:spring_bundle2], Translation(-l*sin(d), 0, 0))
	end
	if p[1] > l * sin(d) # right
		# settransform!(vis[name][:env][:spring_bundle1], Translation(p[1], 0, 0))
		settransform!(vis[name][:env][:spring_bundle1], Translation(l*sin(d), 0, 0))
	else
		settransform!(vis[name][:env][:spring_bundle1], Translation(l*sin(d), 0, 0))
	end

	# settransform!(vis[name][:robot]["contact1"], Translation(p1))

	# r_range = range(0, stop = q[2], length = n_leg)
	# for i = 1:n_leg
	# 	q_tmp = Array(copy(q))
	# 	q_tmp[2] = r_range[i]
	# 	p_leg = [_kinematics(model, q_tmp, mode = :d)[1], 0.0, _kinematics(model, q_tmp, mode = :d)[2]]
    #     settransform!(vis[name][:robot]["arm$i"], Translation(p_leg))
    # end

	return nothing
end

function pusher_arc(center::AbstractVector, radius::Real, α::Real; sense::Symbol=:trig)
	p = center + radius*[cos(α), 0, sin(α)]
	if sense == :trig
		θ = α + π/2
	elseif sense == :antitrig
		θ = α - π/2
	else
		error("Invalid sense.")
	end
	pθ = [p; θ]
	return pθ
end

function pusher_arc_traj(center::AbstractVector, radius::Real, α_start::Real, α_end::Real,
	 	N::Int; sense::Symbol=:trig)
	pθ = []
	for α in range(α_start, stop=α_end, length=N)
		pθ_ = pusher_arc(center, radius, α, sense=sense)
		push!(pθ, pθ_)
	end
	return pθ
end

function generate_pusher_traj(ind::Vector{Int}, w::Vector, x::Vector; side::Symbol=:right,
		center=zeros(3), radius::Real=0.8, N_arc::Int=17)
	H = length(x)
	if side == :right
		sense = :trig
		θ_rest = 0.0
	elseif side == :left
		sense = :antitrig
		θ_rest = π
	else
		error("Invalid side.")
	end

	pθ = [pusher_arc(center, radius, θ_rest; sense=sense) for t=1:H]

	for (i, idx) in enumerate(ind)
		wi = w[i][1]
		θ_impact = x[idx+1][1] + π/2
		θ_next = x[idx+2][1] + π/2

		if ((wi > 0.0) && (side == :right)) || ((wi < 0.0) && (side == :left))
			forward_arc = pusher_arc_traj(center, radius, θ_rest, θ_impact, N_arc; sense=sense)
			backward_arc = pusher_arc_traj(center, radius, θ_next, θ_rest, N_arc; sense=sense)
			arc = [forward_arc; backward_arc] # 2N
			pθ[idx-N_arc-0:idx+N_arc-1] .= arc
		end
	end
	return pθ
end

function build_disturbance!(vis::Visualizer, model::WallPendulum12; name::Symbol=:Pusher, r=0.01, α=1.0)
	r1 = convert(Float32, r)
	r2 = convert(Float32, 6r)
	pusher_mat = MeshPhongMaterial(color = RGBA(1.0, 0.0, 0.0, α))

	# build pusher
	p1 = Cylinder(Point3f0(0.0), Point3f0(0.0, 0.0, 0.150), r1)
	p2 = Cylinder(Point3f0(0.0), Point3f0(0.0, 0.0, 0.03), r2)
	setobject!(vis[name, "pusher1"], p1, pusher_mat)
	setobject!(vis[name, "pusher2"], p2, pusher_mat)
	return anim
end

function set_disturbance!(vis::Visualizer, model::WallPendulum12, pθ::AbstractVector; name::Symbol=:Pusher, offset::Real=0.0)
	p = pθ[1:3]
	θ = pθ[4]
	settransform!(vis[name],
		compose(Translation(p...),
		compose(LinearMap(RotY(-θ-pi/2)),
				Translation(0,0,offset))))
	return nothing
end

function animate_disturbance!(vis::Visualizer, anim::MeshCat.Animation, model::WallPendulum12,
		pθ::AbstractVector; name::Symbol=:Pusher, offset::Real=0.0)
	H = length(pθ)
	for t in 1:H
		MeshCat.atframe(anim, t) do
			set_disturbance!(vis, model, pθ[t], name=name, offset=offset)
		end
	end
	setanimation!(vis, anim)
	return anim
end

function visualize_disturbance!(vis::Visualizer, model::ContactModel, pθ::AbstractVector;
		h=0.01, α=1.0,
		sample=max(1, Int(floor(length(pθ) / 100))),
		anim::MeshCat.Animation=MeshCat.Animation(Int(floor(1/h))),
		name::Symbol=:Pusher, offset::Real=0.0)

	build_disturbance!(vis, model, name=name, α=α)
	animate_disturbance!(vis, anim, model, pθ[1:sample:end], name=name, offset=offset)
	return anim
end


text:
## Wing analysis case
using AeroMDAO

## Wing section setup
wing_right = HalfWing(
                      foils     = fill(naca4((0,0,1,2)), 3),
                      chords    = [2.0, 1.6, 0.2],
                      twists    = [0.0, 0.0, 0.0],
                      spans     = [2.5, 0.5],
                      dihedrals = [5., 5.],
                      sweeps    = [10., 10.],
                      w_sweep   = 0.25 # Quarter-chord
                     );
wing = Wing(wing_right, wing_right)
x_w, y_w, z_w = wing_mac = mean_aerodynamic_center(wing)
print_info(wing, "Wing")

## Meshing and assembly
wing_mesh = WingMesh(wing, [12,6], 6,
                     span_spacing = Cosine()
                    );
aircraft = ComponentVector(wing = make_horseshoes(wing_mesh))

# Freestream conditions
fs  = Freestream(
                 alpha = 1.0,
                 beta  = 0.0,
                 omega = [0.,0.,0.]
                )

# Reference values
refs = References(
                  speed     = 150,
                  density   = 1.225,
                  viscosity = 1.5e-5,
                  area      = projected_area(wing),
                  span      = span(wing),
                  chord     = mean_aerodynamic_chord(wing),
                  location  = mean_aerodynamic_center(wing)
                 )

## Solve system
system  = solve_case(aircraft, fs, refs;
                     print = true
                    );

## Compute dynamics
ax       = Wind()
CFs, CMs = surface_coefficients(system; axes = ax)
# Fs       = surface_forces(system)
# Ms       = surface_moments(system)
# Fs, Ms   = surface_dynamics(system; axes = ax)

## Viscous drag prediction using empirical models

CFs, CMs = surface_coefficients(system; axes = ax)
FFs = farfield_coefficients(system)

# Create array of nearfield and farfield coefficients for each component as a row vector.
comp_coeffs = mapreduce(name -> [ sum(CFs[name]); sum(CMs[name]); FFs[name] ], hcat, keys(system.vortices))

## Equivalent flat-plate skin-friction estimation
x_tr        = fill(0.98, 4)              # Transition locations over sections
CDv_plate   = profile_drag_coefficient(wing_mesh, x_tr, refs)

## Local-dissipation drag estimation (WRONG???)
cam_panels  = camber_panels(wing_mesh)
edge_speeds = surface_velocities(system).wing
CDv_diss    = profile_drag_coefficient(wing_mesh, x_tr, edge_speeds, refs)

## Viscous drag coefficient
CDv = CDv_plate

## Total force coefficients with viscous drag prediction
CDi_nf, CY_nf, CL_nf, Cl, Cm, Cn = nf = nearfield(system)
CDi_ff, CY_ff, CL_ff = ff = farfield(system)

nf_v = (CD = CDi_nf + CDv, CDv = CDv, nf...)
ff_v = (CD = CDi_ff + CDv, CDv = CDv, ff...)

## Evaluate case with stability derivatives
dvs = solve_case_derivatives(aircraft, fs, refs;
                             print_components = true);

## Plotting
using Plots
gr()

## Streamlines

# Spanwise distribution
span_points = 20
init        = chop_trailing_edge(wing, span_points)
dx, dy, dz  = 0, 0, 1e-3
seed        = [ init .+ Ref([dx, dy,  dz])  ;
                init .+ Ref([dx, dy, -dz]) ];

distance = 5
num_stream_points = 100
streams = plot_streamlines(system, seed, distance, num_stream_points);

## Display
horseshoe_panels = chord_panels(wing_mesh)
horseshoe_coords = plot_panels(horseshoe_panels)
wing_coords      = plot_planform(wing);
horseshoe_points = Tuple.(collocation_point.(system.vortices))
ys               = getindex.(horseshoe_points, 2)

## Coordinates
z_limit = 5
plot(xaxis = "x", yaxis = "y", zaxis = "z",
     aspect_ratio = 1,
     camera = (30, 60),
     zlim = (-0.1, z_limit),
     size = (800, 600))
# plot!(wing_coords[:,1], wing_coords[:,2], wing_coords[:,3])
plot!.(horseshoe_coords, color = :black, label = :none)
scatter!(vec(horseshoe_points), marker = 1, color = :black, label = :none)
[ plot!(stream, color = :green, label = :none) for stream in eachcol(Tuple.(streams)) ]
plot!()

## Compute spanwise loads
@code_warntype spanwise_loading(horseshoe_panels, CFs.wing, projected_area(wing))

##
CL_loads   = vec(sum(system.circulations.wing, dims = 1)) / (0.5 * refs.speed * refs.chord)

## Plot spanwise loadings
plot_CD = plot(span_loads[:,1], span_loads[:,2], label = :none, ylabel = "CDi")
plot_CY = plot(span_loads[:,1], span_loads[:,3], label = :none, ylabel = "CY")
plot_CL = begin
            plot(span_loads[:,1], span_loads[:,4], label = :none, xlabel = "y", ylabel = "CL")
            plot!(span_loads[:,1], CL_loads, label = "Normalized", xlabel = "y")
          end
plot(plot_CD, plot_CY, plot_CL, size = (800, 700), layout = (3,1))

## Lift distribution

# Exaggerated CF distribution for plot
hs_pts = vec(Tuple.(bound_leg_center.(system.vortices)))

plot(xaxis = "x", yaxis = "y", zaxis = "z",
     aspect_ratio = 1,
     camera = (60, 60),
     zlim = (-0.1, z_limit),
     title = "Forces (Exaggerated)"
    )
plot!.(horseshoe_coords, color = :gray, label = :none)
# scatter!(cz_pts, zcolor = vec(CLs), marker = 2, label = "CL (Exaggerated)")
quiver!(hs_pts, quiver=(span_loads[:,2], span_loads[:,3], span_loads[:,4]) .* 10)
plot!(size = (800, 600))

## VARIABLE ANALYSES
#=========================================================#

using Setfield
using Base.Iterators: product

## Speed sweep
Vs = 1.0:10:300
res_Vs = permutedims(combinedimsview(
    map(Vs) do V
        ref = @set refs.speed = V
        sys = solve_case(aircraft, fs, ref)
        [ V; farfield(sys)[:]; nearfield(sys) ]
    end
))

plot(
    res_Vs[:,1], res_Vs[:,2:end],
    layout = (3,3), size = (900, 800),
    xlabel = "V",
    labels = ["CD_ff" "CY_ff" "CL_ff" "CD" "CY" "CL" "Cl" "Cm" "Cn"]
)

## Alpha sweep
αs = -5:0.5:5
res_αs = permutedims(combinedimsview(
    map(αs) do α
        fst = @set fs.alpha = deg2rad(α)
        sys = solve_case(aircraft, fst, refs)
        [ α; farfield(sys); nearfield(sys) ]
    end
))

plot(
    res_αs[:,1], res_αs[:,2:end],
    layout = (3,3), size = (900, 800),
    xlabel = "α",
    labels = ["CD_ff" "CY_ff" "CL_ff" "CD" "CY" "CL" "Cl" "Cm" "Cn"]
)

## (Speed, alpha) sweep
res = combinedimsview(
    map(product(Vs, αs)) do (V, α)
        ref = @set refs.speed = V
        fst = @set fs.alpha = deg2rad(α)
        sys = solve_case(aircraft, fst, ref)
        [ V; α; farfield(sys); nearfield(sys) ]
    end
)

##
res_p = permutedims(res, (3,1,2))

# CDi
plt_CDi_ff = plot(camera = (60,45))
[ plot!(
    res_p[:,1,n], res_p[:,2,n], res_p[:,3,n],
    ylabel = "α", xlabel = "V", zlabel = "CDi_ff",
    label = "", c = :black,
) for n in axes(res_p,3) ]

# CL
plt_CL_ff = plot(camera = (45,45))
[ plot!(
    res_p[:,1,n], res_p[:,2,n], res_p[:,8,n],
    ylabel = "α", xlabel = "V", zlabel = "CL_ff",
    label = "", c = :black,
) for n in axes(res_p,3) ]

plt_Cm_ff = plot(camera = (45,45))
[ plot!(
    res_p[:,1,n], res_p[:,2,n], res_p[:,10,n],
    ylabel = "α", xlabel = "V", zlabel = "Cm",
    label = "", c = :black,
) for n in axes(res_p,3) ]

##
plot(plt_CDi_ff, plt_CL_ff, plt_Cm_ff, layout = (1,3), size = (1300, 400))


text:
export Jacobi, Legendre, WeightedJacobi

immutable Jacobi{T,D<:Domain} <: PolynomialSpace{D}
    a::T
    b::T
    domain::D
end
Legendre(domain)=Jacobi(0.,0.,domain)
Legendre()=Legendre(Interval())
Jacobi(a,b,d::Domain)=Jacobi(promote(a,b)...,d)
Jacobi(a,b,d)=Jacobi(a,b,Domain(d))
Jacobi(a,b)=Jacobi(a,b,Interval())
Jacobi(A::Ultraspherical) = Jacobi(order(A)-0.5,order(A)-0.5,domain(A))
Jacobi(A::Chebyshev)=Jacobi(-0.5,-0.5,domain(A))

function Ultraspherical(J::Jacobi)
    if J.a == J.b
        Ultraspherical(J.a+0.5,domain(J))
    else
        error("Cannot construct Ultraspherical with a=$(J.a) and b=$(J.b)")
    end
end

Base.promote_rule{T,V,D}(::Type{Jacobi{T,D}},::Type{Jacobi{V,D}})=Jacobi{promote_type(T,V),D}
Base.convert{T,V,D}(::Type{Jacobi{T,D}},J::Jacobi{V,D})=Jacobi{T,D}(J.a,J.b,J.domain)

typealias WeightedJacobi{T,D} JacobiWeight{Jacobi{T,D},D}

@compat (::Type{WeightedJacobi})(α,β,d::Domain) = JacobiWeight(α,β,Jacobi(β,α,d))
@compat (::Type{WeightedJacobi})(α,β) = JacobiWeight(α,β,Jacobi(β,α))

spacescompatible(a::Jacobi,b::Jacobi) = a.a==b.a && a.b==b.b

function canonicalspace(S::Jacobi)
    if isapproxinteger(S.a+0.5) && isapproxinteger(S.b+0.5)
        Chebyshev(domain(S))
    else
        # return space with parameters in (-1,0.]
        Jacobi(mod(S.a,-1),mod(S.b,-1),domain(S))
    end
end

#####
# jacobirecA/B/C is from dlmf:
# p_{n+1} = (A_n x + B_n)p_n - C_n p_{n-1}
#####
jacobirecA{T}(::Type{T},α,β,k)=k==0&&((α+β==0)||(α+β==-1))?.5*(α+β)+one(T):(2k+α+β+one(T))*(2k+α+β+2one(T))/(2*(k+one(T))*(k+α+β+one(T)))
jacobirecB{T}(::Type{T},α,β,k)=k==0&&((α+β==0)||(α+β==-1))?.5*(α-β)*one(T):(α-β)*(α+β)*(2k+α+β+one(T))/(2*(k+one(T))*(k+α+β+one(T))*(2one(T)*k+α+β))
jacobirecC{T}(::Type{T},α,β,k)=(one(T)*k+α)*(one(T)*k+β)*(2k+α+β+2one(T))/((k+one(T))*(k+α+β+one(T))*(2one(T)*k+α+β))

#####
# jacobirecA/B/C is from dlmf:
# x p_{n-1} =γ_n p_{n-2} + α_n p_{n-1} +  p_n β_n
#####

jacobirecγ{T}(::Type{T},α,β,k)=jacobirecC(T,α,β,k-1)/jacobirecA(T,α,β,k-1)
jacobirecα{T}(::Type{T},α,β,k)=-jacobirecB(T,α,β,k-1)/jacobirecA(T,α,β,k-1)
jacobirecβ{T}(::Type{T},α,β,k)=1/jacobirecA(T,α,β,k-1)

for (REC,JREC) in ((:recα,:jacobirecα),(:recβ,:jacobirecβ),(:recγ,:jacobirecγ),
                   (:recA,:jacobirecA),(:recB,:jacobirecB),(:recC,:jacobirecC))
    @eval $REC{T}(::Type{T},sp::Jacobi,k)=$JREC(T,sp.a,sp.b,k)
end

function jacobip(r::Range,α,β,x)
    n=r[end]+1
    if n<=2
        v=[1.,.5*(α-β+(2+α+β)*x)]
    else
        T=promote_type(Float64,typeof(x))
        v=Vector{T}(n)  # x may be complex
        v[1]=1.
        v[2]=.5*(α-β+(2+α+β)*x)

        for k=2:n-1
            v[k+1]=((x-jacobirecα(T,α,β,k))*v[k] - jacobirecγ(T,α,β,k)*v[k-1])/jacobirecβ(T,α,β,k)
        end
    end
    v[r+1]
end

function jacobip(r::Range,α,β,x::Number)
    if x==1. && α==0.
        ones(length(r))
    elseif x==-1. && β==0.
        (-1.).^r
    else
        n=r[end]+1
        if n<=2
            v=[1.,.5*(α-β+(2+α+β)*x)]
        else
            T=promote_type(Float64,typeof(x))
            v=Vector{T}(n)  # x may be complex
            v[1]=1.
            v[2]=.5*(α-β+(2+α+β)*x)

            for k=2:n-1
                v[k+1]=((x-jacobirecα(T,α,β,k))*v[k] - jacobirecγ(T,α,β,k)*v[k-1])/jacobirecβ(T,α,β,k)
            end
        end
        v[r+1]
    end
end
jacobip(n::Integer,α,β,v)=jacobip(n:n,α,β,v)[1]
jacobip(n::Range,α,β,v::Vector)=transpose(hcat(map(x->jacobip(n,α,β,x),v)...))
jacobip(n::Integer,α,β,v::Vector)=map(x->jacobip(n,α,β,x),v)
jacobip(n,S::Jacobi,v)=jacobip(n,S.a,S.b,v)

include("jacobitransform.jl")
include("JacobiOperators.jl")

for op in (:(Base.ones),:(Base.zeros))
    @eval ($op){T<:Number}(::Type{T},S::Jacobi)=Fun(($op)(T,1),S)
    @eval ($op)(S::Jacobi)=Fun(($op)(1),S)
end

function identity_fun(J::Jacobi)
    if domain(J)==Interval()
        Fun([(J.b-J.a)/(2+J.a+J.b),2.0/(2+J.a+J.b)],J)
    else
        d=domain(J)
        complexlength(d)/2*(Fun([(J.b-J.a)/(2+J.a+J.b),2.0/(2+J.a+J.b)],J)+1.)+first(d)
    end
end

setdomain(S::Jacobi,d::Domain)=Jacobi(S.a,S.b,d)

# O(min(m,n)) Jacobi conjugated inner product

function conjugatedinnerproduct{S,V}(sp::Jacobi,u::Vector{S},v::Vector{V})
    T,mn = promote_type(S,V),min(length(u),length(v))
    α,β = sp.a,sp.b
    if mn > 1
        wi = 2^(α+β+1)*gamma(α+1)*gamma(β+1)/gamma(α+β+2)
        ret = u[1]*wi*v[1]
        for i=2:mn
            wi *= (α+i-1)*(β+i-1)/(i-1)/(i-1+α+β)*(2i+α+β-3)/(2i+α+β-1)
            ret += u[i]*wi*v[i]
        end
        return ret
    elseif mn > 0
        wi = 2^(α+β+1)*gamma(α+1)*gamma(β+1)/gamma(α+β+2)
        return u[1]*wi*v[1]
    else
        return zero(promote_type(eltype(u),eltype(v)))
    end
end

function bilinearform{J<:Jacobi}(f::Fun{J},g::Fun{J})
    @assert domain(f) == domain(g)
    if f.space.a == g.space.a == 0. && f.space.b == g.space.b == 0.
        return complexlength(domain(f))/2*conjugatedinnerproduct(g.space,f.coefficients,g.coefficients)
    else
        return defaultbilinearform(f,g)
    end
end

function bilinearform{J<:Jacobi,DD<:Interval}(f::Fun{JacobiWeight{J,DD}},g::Fun{J})
    @assert domain(f) == domain(g)
    if f.space.β == f.space.space.a == g.space.a && f.space.α == f.space.space.b == g.space.b
        return complexlength(domain(f))/2*conjugatedinnerproduct(g.space,f.coefficients,g.coefficients)
    else
        return defaultbilinearform(f,g)
    end
end

function bilinearform{J<:Jacobi,DD<:Interval}(f::Fun{J},g::Fun{JacobiWeight{J,DD}})
    @assert domain(f) == domain(g)
    if g.space.β == g.space.space.a == f.space.a && g.space.α == g.space.space.b == f.space.b
        return complexlength(domain(f))/2*conjugatedinnerproduct(f.space,f.coefficients,g.coefficients)
    else
        return defaultbilinearform(f,g)
    end
end

function bilinearform{J<:Jacobi,DD<:Interval}(f::Fun{JacobiWeight{J,DD}},g::Fun{JacobiWeight{J,DD}})
    @assert domain(f) == domain(g)
    if f.space.β + g.space.β == f.space.space.a == g.space.space.a && f.space.α + g.space.α == f.space.space.b == g.space.space.b
        return complexlength(domain(f))/2*conjugatedinnerproduct(f.space.space,f.coefficients,g.coefficients)
    else
        return defaultbilinearform(f,g)
    end
end

function linebilinearform{J<:Jacobi}(f::Fun{J},g::Fun{J})
    @assert domain(f) == domain(g)
    if f.space.a == g.space.a == 0. && f.space.b == g.space.b == 0.
        return arclength(domain(f))/2*conjugatedinnerproduct(g.space,f.coefficients,g.coefficients)
    else
        return defaultlinebilinearform(f,g)
    end
end

function linebilinearform{J<:Jacobi,DD<:Interval}(f::Fun{JacobiWeight{J,DD}},g::Fun{J})
    @assert domain(f) == domain(g)
    if f.space.β == f.space.space.a == g.space.a && f.space.α == f.space.space.b == g.space.b
        return arclength(domain(f))/2*conjugatedinnerproduct(g.space,f.coefficients,g.coefficients)
    else
        return defaultlinebilinearform(f,g)
    end
end

function linebilinearform{J<:Jacobi,DD<:Interval}(f::Fun{J},g::Fun{JacobiWeight{J,DD}})
    @assert domain(f) == domain(g)
    if g.space.β == g.space.space.a == f.space.a && g.space.α == g.space.space.b == f.space.b
        return arclength(domain(f))/2*conjugatedinnerproduct(f.space,f.coefficients,g.coefficients)
    else
        return defaultlinebilinearform(f,g)
    end
end

function linebilinearform{J<:Jacobi,DD<:Interval}(f::Fun{JacobiWeight{J,DD}},g::Fun{JacobiWeight{J,DD}})
    @assert domain(f) == domain(g)
    if f.space.β + g.space.β == f.space.space.a == g.space.space.a && f.space.α + g.space.α == f.space.space.b == g.space.space.b
        return arclength(domain(f))/2*conjugatedinnerproduct(f.space.space,f.coefficients,g.coefficients)
    else
        return defaultlinebilinearform(f,g)
    end
end


text:
struct OperationModelResults
    variables::Dict{Symbol, DataFrames.DataFrame}
    total_cost::Dict{Symbol, Any}
    optimizer_log::Dict{Symbol, Any}
    time_stamp::DataFrames.DataFrame

end

function get_variable(res_model::OperationModelResults, key::Symbol)
        try
            !isnothing(res_model.variables)
        catch
            error("No variable with key $(key) has been found.")
        end
    return get(res_model.variables, key, nothing)
end

function get_optimizer_log(res_model::OperationModelResults)
    return res_model.optimizer_log
end

function get_time_stamp(res_model::OperationModelResults, key::Symbol)
    return res_model.time_stamp
end

"""
    results = load_operation_results(path, folder_name)

This function can be used to load results from a folder
of results from a single-step problem, or for a single foulder
within a simulation.

# Arguments
-`path::AbstractString = folder path`
-`directory::AbstractString = "2019-10-03T09-18-00"`: the foulder name that contains
feather files of the results.

# Example
```julia
results = load_operation_results("/Users/test/", "2019-10-03T09-18-00")
```
"""

function load_operation_results(path::AbstractString, directory::AbstractString)

    if isfile(path)
        path = dirname(path)
    end

    folder_path = joinpath(path, directory)
    files_in_folder = collect(readdir(folder_path))

        variables = setdiff(files_in_folder, ["time_stamp.feather", "optimizer_log.feather"])
        variable_dict = Dict{Symbol, DataFrames.DataFrame}()

        for i in 1:length(variables)

            variable = variables[i]
            variable_name = split("$variable", ".feather")[1]
            file_path = joinpath(folder_path,"$variable_name.feather")
            variable_dict[Symbol(variable_name)] = Feather.read("$file_path") #change key to variable

        end

        file_path = joinpath(folder_path,"optimizer_log.feather")
        optimizer = Dict{Symbol, Any}(eachcol(Feather.read("$file_path"),true))

        file_path = joinpath(folder_path,"time_stamp.feather")
        temp_time_stamp = Feather.read("$file_path")
        time_stamp = temp_time_stamp[1:(size(temp_time_stamp,1)-1),:]

        obj_value = Dict{Symbol, Any}(:OBJECTIVE_FUNCTION => optimizer[:obj_value])
        results = OperationModelResults(variable_dict, obj_value, optimizer, time_stamp)

    return results

end
# This function returns the length of time_stamp for each step
# that is unique to that step and not overlapping with the next.

function _count_time_overlap(stage::String,
                            variable::Array,
                            references::Dict{Any,Any})

    date_df = references[stage][variable[1]]
    ref = DataFrames.DataFrame()

    for (ix,time) in enumerate(date_df.Date)

        try
            file_path =  date_df[ix, :File_Path]
            time_file_path = joinpath(dirname(file_path), "time_stamp.feather")
            temp_time_stamp = DataFrames.DataFrame(Feather.read("$time_file_path"))
            t = size(temp_time_stamp, 1)
            append!(ref,temp_time_stamp[(1:t-1),:])

        catch
            @warn "The given date_range is outside the results time stamp."
        end
    end
    if size(unique(ref),2) == size(ref,2)
        return 0.0
    else
        extra_time_length = size(unique(ref),1)./(length(step)+1)
        return extra_time_length
    end
    end

"""
    results = load_simulation_results(stage, step, date_range,variable,references)
This function goes through the reference table of file paths and
aggregates the results over time into a struct of type OperationModelResults
**Note:** the array of steps should match the date range provided.
# Arguments
-`stage::String = "stage-1"``: The stage of the results getting parsed, stage-1 or stage-2
-`step::Array{String} = ["step-1", "step-2", "step-3"]`: the steps of the results getting parsed
-`date_range::StepRange = 2020/01/01T:00:00:00 : 2020/01/03:00:00:00`: the date range to be parsed
-`variable::Array{Symbol} = [:P_ThermalStandard, :P_RenewableDispatch]`: the variables to be parsed
# Example
```julia
date_range = (Dates.DateTime(2020, April, 4):Dates.Hour(24):Dates.DateTime(2020, April, 6))
stage = "stage-1"
step = ["step-1","step-2", "step-3"] # has to match the date range
variable = [:P_ThermalStandard, :P_RenewableDispatch]
results = load_simulation_results(stage,step, date_range, variable, references)
```
"""
function load_simulation_results(stage::String,
                                 references::Dict{Any,Any};
                                 kwargs...)

    variable_dict = Dict()
    variable = collect(keys(references[stage]))
    time_stamp = DataFrames.DataFrame(Range = Dates.DateTime[])
    extra_time_length = _count_time_overlap(stage,
                                            variable,
                                            references)
    for l in 1:length(variable)
        date_df = references[stage][variable[l]]
        variable_dict[(variable[l])] = DataFrames.DataFrame()

        for (ix,time) in enumerate(date_df.Date)

            file_path = date_df[ix, :File_Path]
            var = Feather.read("$file_path")
            correct_var_length = size(1:(size(var,1) - extra_time_length),1)
            variable_dict[(variable[l])] = vcat(variable_dict[(variable[l])],var[1:correct_var_length,:])
            if l == 1
                time_file_path = joinpath(dirname(file_path), "time_stamp.feather")
                temp_time_stamp = DataFrames.DataFrame(Feather.read("$time_file_path"))
                non_overlap = size((1:size(temp_time_stamp, 1) - extra_time_length - 1),1)
                time_stamp = vcat(time_stamp,temp_time_stamp[(1:non_overlap),:])

            end
        end

    end

    # first_file = references[stage][variable[1]]
    # file_path = first_file[first_file.Date .== date_range[1], :File_Path][1]
    # opt_file_path = joinpath(dirname(file_path),"optimizer_log.feather")
    optimizer = Dict{Symbol, Any}()
    obj_value = Dict{Symbol, Any}(:OBJECTIVE_FUNCTION => 0.0)
    results = OperationModelResults(variable_dict, obj_value, optimizer, time_stamp)
    file_path = references[stage][variable[1]][1,:File_Path]
    if (:write in keys(kwargs)) == true
        write_model_results(results, dirname(dirname(dirname(dirname(file_path)))),"results")
    end

    return results

end


text:
#!/usr/local/bin/julia

if length(ARGS) != 1
    @warn "usage: jl_to_md.jl jl_filename.jl"
    exit(1)
end

# Default methods.

enter(a...) = nothing
format(::Any, n, line) = line
leave(a...) = nothing

# Content Types.

struct Code end
struct Heading end
struct Doc end
struct FixMe end

format(::FixMe, n, line) = "FIXME[^FIXME$n]\n\n[^FIXME$n]: ⚠ $line"

#code_start(n) = """```julia"""
code_start(n) = """```{.julia .numberLines .lineAnchors startFrom="$n"}"""
code_end = "```"

enter(::Code, v, n) = push!(v, code_start(n))

function leave(::Code, v)
    while last(v) == ""
        pop!(v)
    end
    push!(v, code_end, "")
end

"""
Read C source code from `c_filename`.
Convert C comments to Markdown.
Convert C code to Markdown code blocks.
Write output to `\$c_filename.md`.
"""
function c_to_md(c_filename)

    md_filename = "$c_filename.md"

    v = String[]
    state = Code

    for (n, line) in enumerate(eachline(c_filename))

        for ((pattern_state, pattern), next_state) in [
            (Doc,   r"^ \"\"\"                         $"x) => nothing,
            (Doc,   r"^          (.*?)                 $"x) => Doc(),
            (Any,   r"^ [#] [ ]* FIXME (.*?) [ ]*      $"x) => FixMe(),
            (Any,   r"^ ([#]+ [ ] .+?) [ ]*            $"x) => Heading(),
            (Any,   r"^ [@]doc [ ] README\"\"\"        $"x) => Doc(),
            (Any,   r"^ [a-z]* \"\"\"                  $"x) => Doc(),
            (Code,  r"^          (.*?)                 $"x) => Code(),
            (Any,   r"^          (.+?)                 $"x) => Code()
        ]
            if state isa pattern_state
                m = match(pattern, line)
                m == nothing && continue
                state == next_state || leave(state, v)
                if state != next_state
                    state = next_state
                    push!(v, "")
                    enter(state, v, n)
                end
                if !isempty(m.captures)
                    push!(v, format(state, n, m[1]))
                end
                line = nothing
                break
            end
        end
        if line != nothing
            push!(v, line)
        end
    end
    leave(state, v)

    i = 1
    if v[i] == "---"
        i = 1 + findnext(isequal("---"), v, 2)
    end
#    f = basename(c_filename)
#    insert!(v, i, "")
#    insert!(v, i, "*Raw source: [`$f`]($(f)!)*")

    write(md_filename, join(v, "\n"))
end

c_to_md(ARGS[1])

max_stars_repo_name: notinaboat/ReadmeDocs

max_issues_repo_name: notinaboat/ReadmeDocs

max_forks_repo_name: notinaboat/ReadmeDocs


text:
# Tests for Revenue DEA Models
@testset "RevenueDEAModel" begin

    ## Test Revenue DEA Model with Cooper et al. (2007)
    # Test agains book results
    X = [3 2; 1 3; 4 6]
    Y = [3; 5; 6]
    P = [6; 6; 6]

    dearevenuecooper = dearevenue(X, Y, P, rts = :CRS)

    @test typeof(dearevenuecooper) == RevenueDEAModel

    @test nobs(dearevenuecooper) == 3
    @test ninputs(dearevenuecooper) == 2
    @test noutputs(dearevenuecooper) == 1
    @test ismonetary(dearevenuecooper) == false

    @test efficiency(dearevenuecooper, :Economic)   ≈ [0.9; 1; 0.6] atol = 1e-3
    @test efficiency(dearevenuecooper, :Technical)  ≈ [0.9; 1; 0.6  ] atol = 1e-3
    @test efficiency(dearevenuecooper, :Allocative) ≈ [1; 1; 1] atol = 1e-3

    ## Test Revenue DEA Model with Zofío and Prieto (2006) data.
    # Test agains results in R
    X = [5 3; 2 4; 4 2; 4 8; 7 9]
    Y = [7 4; 10 8; 8 10; 5 4; 3 6]
    P = [3 2; 3 2; 3 2; 3 2; 3 2.0]

    # Revnue CRS
    dearevenuecrs = dearevenue(X, Y, P, rts = :CRS)
    @test efficiency(dearevenuecrs) ≈ [0.4915254237;
                                    1.000;
                                    1.000;
                                    0.250;
                                    0.1735537190]
    @test efficiency(dearevenuecrs, :Technical) ≈ [0.6364;
                                1.000;
                                1.000;
                                0.250;
                                0.2609]  atol = 1e-3
    @test efficiency(dearevenuecrs, :Allocative) ≈ [0.7723970944;
                                1.000;
                                1.000;
                                1.000;
                                0.6652892562]

    @test efficiency(dearevenue(targets(dearevenuecrs, :X), targets(dearevenuecrs, :Y), P, rts = :CRS)) ≈ ones(5)

    # Revenue VRS
    dearevenuevrs = dearevenue(X, Y, P, rts = :VRS)
    @test efficiency(dearevenuevrs) ≈ [0.6444444444 ;
                                    1.000;
                                    1.000;
                                    0.500;
                                    0.4565217391]
    @test efficiency(dearevenuevrs, :Technical) ≈ [0.7777777778 ;
                                1.000;
                                1.000;
                                0.500 ;
                                0.6000000000]
    @test efficiency(dearevenuevrs, :Allocative) ≈ [0.8285714286;
                                1.000;
                                1.000;
                                1.000;
                                0.7608695652]

    @test efficiency(dearevenue(targets(dearevenuevrs, :X), targets(dearevenuevrs, :Y), P, rts = :VRS)) ≈ ones(5)

    # Check defaults
    @test efficiency(dearevenue(X, Y, P)) == efficiency(dearevenuevrs)
    @test efficiency(dearevenuevrs, :Economic) == efficiency(dearevenuevrs)

    # Print
    show(IOBuffer(), dearevenuecooper)

    # Test errors
    @test_throws DimensionMismatch dearevenue([1; 2 ; 3], [4 ; 5], [1; 1; 1]) #  Different number of observations
    @test_throws DimensionMismatch dearevenue([1; 2; 3], [4; 5; 6], [1; 2; 3; 4]) # Different number of observation in prices
    @test_throws DimensionMismatch dearevenue([1; 2; 3], [4 4; 5 5; 6 6], [4 4 4; 5 5 5; 6 6 6]) # Different number of output prices and outputs
    @test_throws ArgumentError dearevenue([1; 2; 3], [4; 5; 6], [1; 2; 3], rts = :Error) # Invalid returns to scale
    @test_throws ArgumentError dearevenue([1; 2; 3], [4; 5; 6], [1; 2; 3], dispos = :Error) # Invalid disposability
    @test_throws ArgumentError normfactor(dearevenue(X, Y, P)) # ERROR: RevenueDEAModel has no normalization factor

    # ------------------
    # Weak Disposability Tests
    # ------------------

    X = [1; 2; 3; 2; 4]
    Y = [2; 3; 4; 1; 3]
    P = [1; 1; 1; 1; 1]

    dearevenueStrong = dearevenue(X, Y, P, dispos = :Strong)
    @test efficiency(dearevenueStrong, :Economic) ≈ [1.0; 1.0; 1.0; 0.3333333333333333; 0.75]
    @test efficiency(dearevenueStrong, :Technical) ≈ [1.0; 1.0; 1.0; 0.3333333333333333; 0.75]
    @test efficiency(dearevenueStrong, :Allocative) ≈ [1.0; 1.0; 1.0; 1.0; 1.0]

    dearevenueWeak = dearevenue(X, Y, P, dispos = :Weak)
    @test efficiency(dearevenueWeak, :Economic) ≈ [1.0; 1.0; 1.0; 0.3333333333333333; 1.0]
    @test efficiency(dearevenueWeak, :Technical) ≈ [1.0; 1.0; 1.0; 0.3333333333333333; 1.0]
    @test efficiency(dearevenueWeak, :Allocative) ≈ [1.0; 1.0; 1.0; 1.0; 1.0]

    # ------------------
    # Test Vector and Matrix inputs and outputs
    # ------------------
    # Tests against results in R

    # Inputs is Matrix, Outputs is Vector
    X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6	8]
    Y = [1; 1; 1; 1; 1; 1; 1; 1]
    P = [1; 1; 1; 1; 1; 1; 1; 1]

    @test efficiency(dearevenue(X, Y, P)) ≈ [1; 1; 1; 1; 1; 1; 1; 1]

    # Inputs is Vector, Output is Matrix
    X = [1; 1; 1; 1; 1; 1; 1; 1]
    Y = [7 7; 4 8; 8 4; 3 5; 3 3; 8 2; 6 4; 1.5 5]
    P = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1]

    @test efficiency(dearevenue(X, Y, P)) ≈ [1; 0.8571428571; 0.8571428571; 0.5714285714; 0.4285714286; 0.7142857143; 0.7142857143; 0.4642857143]

    # Inputs is Vector, Output is Vector
    X = [2; 4; 8; 12; 6; 14; 14; 9.412]
    Y = [1; 5; 8; 9; 3; 7; 9; 2.353]
    P = [1; 1; 1; 1; 1; 1; 1; 1]

    @test efficiency(dearevenue(X, Y, P)) ≈ [1; 1; 1; 1; 0.4615384615; 0.7777777778; 1; 0.2816951993]

end