text:
# Autogenerated wrapper script for libCZI_jll for x86_64-linux-gnu-cxx11
export CZIcmd, libCZI

using Zlib_jll
using libpng_jll
JLLWrappers.@generate_wrapper_header("libCZI")
JLLWrappers.@declare_executable_product(CZIcmd)
JLLWrappers.@declare_library_product(libCZI, "liblibCZI.so")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, libpng_jll)
    JLLWrappers.@init_executable_product(
        CZIcmd,
        "bin/CZIcmd",
    )

    JLLWrappers.@init_library_product(
        libCZI,
        "lib/liblibCZI.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Test, Agents, Random, Graphs, DataFrames
using StatsBase:mean
using StableRNGs

mutable struct Agent0 <: AbstractAgent
    id::Int
end

mutable struct Agent1 <: AbstractAgent
    id::Int
    pos::Dims{2}
end

mutable struct Agent2 <: AbstractAgent
    id::Int
    weight::Float64
end

mutable struct Agent3 <: AbstractAgent
    id::Int
    pos::Dims{2}
    weight::Float64
end

mutable struct Agent4 <: AbstractAgent
    id::Int
    pos::Dims{2}
    p::Int
end

mutable struct Agent5 <: AbstractAgent
    id::Int
    pos::Int
    weight::Float64
end

mutable struct Agent6 <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    vel::NTuple{2,Float64}
    weight::Float64
end

mutable struct Agent7 <: AbstractAgent
    id::Int
    pos::Int
    f1::Bool
    f2::Int
end

Agent7(id, pos; f1, f2) = Agent7(id, pos, f1, f2)

mutable struct Agent8 <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    f1::Bool
    f2::Int
end

Agent8(id, pos; f1, f2) = Agent8(id, pos, f1, f2)

mutable struct Agent9 <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    vel::NTuple{2,Float64}
    f1::Union{Int,Nothing}
end

mutable struct Agent10 <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int,Float64}
end

@testset "Agents.jl Tests" begin
    include("model_creation_tests.jl")
    include("api_tests.jl")
    include("randomness_tests.jl")
    include("scheduler_tests.jl")
    include("model_access.jl")
    include("space_test.jl")
    include("collect_tests.jl")
    include("continuousSpace_tests.jl")
    include("osm_tests.jl")
    include("astar_tests.jl")
    include("collisions_tests.jl")
    include("graph_tests.jl")
    include("csv_tests.jl")
    include("jld2_tests.jl")
end

ext: jl


text:
#=

    This is a port of the R2 model DigitRecognition.cs

    with the datafiles
    - digit_recognotion_input.csv  28 x 28
    - digit_recognotion_params.csv 28 x 28 x 10
    - digit_recognotion_prior.csv  for 0..9

    The image (that would be a 3):
        0000000000000000000000000000
        0000000000001111100000000000
        0000000111111111111000000000
        0000001111111111111000000000
        0000011111111111111000000000
        0000001111111111111000000000
        0000001110000001111000000000
        0000000000000001111000000000
        0000000000000011111000000000
        0000000000000111111000000000
        0000000000001111110000000000
        0000000000011111110000000000
        0000000000111111111000000000
        0000000000111111111100000000
        0000000000011111111110000000
        0000000000000000111110000000
        0000000000000000011110000000
        0000000000000000011110000000
        0000000000000000111110000000
        0000000000000011111110000000
        0011000000001111111100000000
        0111111111111111111000000000
        0111111111111111110000000000
        0111111111111111000000000000
        0011111111111100000000000000
        0000000000000000000000000000
        0000000000000000000000000000
        0000000000000000000000000000

    The R2 model outputs:
    ```
    Mean: 3.021
    Variance: 0.0405996
    Number of accepted samples = 89
    ```

    This Turing.jl model, were y0 is the digits (0..9) (y is 1..10)

    Summary Statistics
    parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
        Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

            y    3.9970    0.0949     0.0030    0.0030   1000.0080    1.0000       59.9238
            y0    2.9970    0.0949     0.0030    0.0030   1000.0080    1.0000       59.9238

    Distributions of variable y (num:0)
    4.00000 =>     999  (0.999000)
    1.00000 =>       1  (0.001000)
    Distributions of variable y0 (num:0)
    3.00000 =>     999  (0.999000)
    0.00000 =>       1  (0.001000)

=#

using Turing, StatsPlots, Distributions, StatsBase
using CSV
include("jl_utils.jl")

@model function digit_recognition(x,params,prior)
    x_len = length(x)

    y ~ Categorical(simplex(prior)) # this is 1..10
    # A uniform prior give the same result
    # y ~ Categorical(simplex([1/10 for i in 1:10]))
    y0 ~ Dirac(y-1) # adjust to digits 0..9
    for i in 1:x_len
        x[i] ~ Bernoulli(params[y,i])
    end
end

input = parse.(Int,split(readline("digit_recognition_input.csv"),","))
input_img = reshape(input,28,28)
for i in 1:28
    println(join(string.(input_img)[i,:]))
end
println()

params = Matrix(CSV.read("digit_recognition_params.csv",header=false,DataFrame))
prior = parse.(Float64,split(readline("digit_recognition_prior.csv"),","))

println("prior:$prior len:$(length(prior)) sum(prior):$(sum(prior))")

model = digit_recognition(input,params,prior)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 1_000)
chns = sample(model, PG(5), 1_000)
# chns = sample(model, SMC(), 1_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, HMC(0.1,6), 1_000)
# chns = sample(model, NUTS(), 1_000)

display(chns)

show_var_dist_pct(chns,:y)
show_var_dist_pct(chns,:y0)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Interpolations, Test, LinearAlgebra, Random

@testset "ScalingNoInterpTests" begin
    xs = -pi:2pi/10:pi
    f1(x) = sin(x)
    f2(x) = cos(x)
    f3(x) = sin(x) .* cos(x)
    f(x,y) = y == 1 ? f1(x) : (y == 2 ? f2(x) : (y == 3 ? f3(x) : error("invalid value for y (must be 1, 2 or 3, you used $y)")))
    ys = 1:3

    A = hcat(map(f1, xs), map(f2, xs), map(f3, xs))

    itp = interpolate(A, (BSpline(Quadratic(Periodic(OnGrid()))), NoInterp()))
    sitp = scale(itp, xs, ys)

    for (ix,x0) in enumerate(xs[1:end-1]), y0 in ys
        x,y = x0, y0
        @test ≈(sitp(x,y),f(x,y),atol=0.05)
    end

    @test length(Interpolations.gradient(sitp, pi/3, 2)) == 1

    # check for case where initial/middle indices are NoInterp but later ones are <:BSpline
    isdefined(Random, :seed!) ? Random.seed!(1234) : srand(1234) # `srand` was renamed to `seed!`
    z0 = rand(10,10)
    za = copy(z0)
    zb = copy(z0')

    itpa = interpolate(za, (BSpline(Linear()), NoInterp()))
    itpb = interpolate(zb, (NoInterp(), BSpline(Linear())))

    rng = range(1.0, stop=19.0, length=10)
    sitpa = scale(itpa, rng, 1:10)
    sitpb = scale(itpb, 1:10, rng)
    @test Interpolations.gradient(sitpa, 3.0, 3) ==  Interpolations.gradient(sitpb, 3, 3.0)

end


text:
include("TableProcessing.jl")
include("FindKm.jl")
include("Simulate.jl")
include("Plots.jl")
include("pValues.jl")
import Cairo
import Fontconfig
import DataFrames
using Distributions
using DelimitedFiles
using HDF5, JLD

#Set probability(FP|N) and probability(FN|P)
pFN = 0.04
pFP = 0.01

#Function that performs global level predictions
function pDropoutInfo(countsMatrixPath::String, GeneIsoTablePath::String,
      output::String, filteringThreshold::Int64)

    #Read in counts matrix
    println("Data wrangling commencing!")
    println("Reading in counts matrix.") #Quick
    countsMatrix = ReadCountsMatrix(countsMatrixPath)

    #Read in gene-isoform table
    println("Reading in gene-isoform relationships")
    GeneIsoformRelationships = ReadGeneIsoformRelationships(GeneIsoTablePath)

    #Make sure gene-isoform table matches counts matrix
    if CheckGenesMatchIsoforms(countsMatrix, GeneIsoformRelationships) == false
        throw(ArgumentError("Your counts matrix doesn't match your list of genes
        and isoforms - there is something wrong with your input files."))
    end

    #Add gene information to isoforms table
    println("Adding Gene information to countsMatrix")
    countsMatrixWithGenes = AddGeneInfo(countsMatrix, GeneIsoformRelationships)

    #Clean data for finding Km
    println("Cleaning data for finding Km")
    MMdf = CleanDataForMM(countsMatrix)

    #Find Km
    println("Calculating Km....")
    MMParams = FindKm(MMdf)
    Km = MMParams[1]
    sigma = MMParams[2]

    println("Converting counts to dict")
    GenesWithNIsoforms = KeepGenesWithNIsoforms(countsMatrixWithGenes,
    filteringThreshold)
    workingDict = GenesWithNIsoforms[1]
    unexpr_dict = GenesWithNIsoforms[2]
    MaxNumIsos = GenesWithNIsoforms[3]

    #Get ranking and dropout info for each gene
    pDropoutArr::Array{Float64} = []
    pDropoutDict = Dict()
    ExprArr::Array{Float64} =[]
    ExprDict = Dict()

    for (key, value) in workingDict
        ExprArr = []
        for (k,v) in value
            S = mean(v)
            push!(ExprArr, S)
        end

        append!(pDropoutArr, FindPDropout(value, Km))
        pDropoutDict[key] = FindPDropout(value, Km)
        ExprDict[key] = ExprArr
    end

    save(output * "_dict.jld", "data", pDropoutDict)
    save(output * "expr_dict.jld", "data", ExprDict)

    pDropoutDistributionPlot = pDropoutDistribution(pDropoutArr, output)
    params = fit(Beta, pDropoutArr)
    return (params.α, params.β, pDropoutDistributionPlot)
end

H9_96 = pDropoutInfo("data/H9_96_scnorm_counts.txt", "data/gencode_human_gene_isoform_relationships.txt", "figures/pdrop_scnorm_H9_96.png", 4)
H9_24 = pDropoutInfo("data/H9_24_scnorm_counts.txt", "data/gencode_human_gene_isoform_relationships.txt", "figures/pdrop_scnorm_H9_24.png", 4)
H1_96 = pDropoutInfo("data/H1_96_scnorm_counts.txt", "data/gencode_human_gene_isoform_relationships.txt", "figures/pdrop_scnorm_H1_96.png", 4)
H1_24 = pDropoutInfo("data/H1_24_scnorm_counts.txt", "data/gencode_human_gene_isoform_relationships.txt", "figures/pdrop_scnorm_H1_24.png", 4)

p1 = plot([H1_24[1], H1_96[1], H9_96[1], H9_24[1]], [H1_24[2], H1_96[2], H9_96[2], H9_24[2]], seriestype=:scatter)
savefig(p1, "figures/scnorm_beta_params.png")

v0 = pdf(Beta(0.45, 1.45), 0:0.01:1)
v1 = pdf(Beta(0.7, 1.2), 0:0.01:1)
v2 = pdf(Beta(0.95, 0.95), 0:0.01:1)
v3 = pdf(Beta(1.2, 0.7), 0:0.01:1)
v4 = pdf(Beta(1.45, 0.45), 0:0.01:1)
p_all = plot(0:0.01:1, [v0,v1,v2,v3,v4], label = ["alpha=0.45, beta=1.45",
"alpha=0.7, beta=1.2",
"alpha=0.95, beta=0.95",
"alpha=1.2, beta=0.7",
"alpha=1.45, beta=0.45"])
savefig(p_all, "figures/proposed_beta_distributions.png")

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
module Pandas

using PyCall
using PyPlot

import Base.getindex, Base.setindex!, Base.length, Base.size, Base.mean, Base.std, Base.show, Base.merge, Base.convert, Base.hist, Base.join, Base.replace, Base.endof, Base.start, Base.next, Base.done, Base.sum, Base.var
import PyPlot.plot

export DataFrame, Iloc, Loc, Ix, Series, MultiIndex, Index, GroupBy
export mean, std, agg, aggregate, median, var, ohlc, transform, groups, indices, get_group
export iloc,loc,reset_index,index,head,xs,plot,hist,join,align,drop,drop_duplicates,duplicated,filter,first,idxmax,idxmin,last,reindex,reindex_axis,reindex_like,rename,tail,set_index,select,take,truncate,abs,any,clip,clip_lower,clip_upper,corr,corrwith,count,cov,cummax,cummin,cumprod,cumsum,describe,diff,mean,median,min,mode,pct_change,rank,quantile,sum,skew,var,std,dropna,fillna,replace,delevel,pivot,reodrer_levels,sort,sort_index,sortlevel,swaplevel,stack,unstack,T,boxplot
export to_clipboard,to_csv,to_dense,to_dict,to_excel,to_gbq,to_hdf,to_html,to_json,to_latex,to_msgpack,to_panel,to_pickle,to_records,to_sparse,to_sql,to_string,query, groupby, columns, app, values, from_arrays, from_tuples
export read_csv, read_html, read_json, read_excel, read_table, save, stats,  melt, rolling_count, rolling_sum, rolling_window, rolling_quantile, ewma, setcolumns!, concat, read_pickle
export pivot_table, crosstab, cut, qcut, get_dummies, deletecolumn!, siz, name, setname!
export argsort,order,asfreq,asof,shift,first_valid_index,last_valid_index,weekday,resample,tz_conert,tz_localize
export resample,date_range,to_datetime,to_timedelta,bdate_range,period_range,ewma,ewmstd,ewmvar,ewmcorr,ewmcov
export rolling_count, expanding_count, rolling_sum, expanding_sum, rolling_mean, expanding_mean, rolling_median, expanding_median, rolling_var, expanding_var, rolling_std, expanding_std, rolling_min, expanding_min, rolling_max, expanding_max, rolling_corr, expanding_corr, rolling_corr_pairwise, expanding_corr_pairwise, rolling_cov, expanding_cov, rolling_skew, expanding_skew, rolling_kurt, expanding_kurt, rolling_apply, expanding_apply, rolling_quantile, expanding_quantile

np = pyimport("numpy")
pandas_raw = pyimport("pandas")
pandas_mod = pywrap(pandas_raw)
type_map = []

abstract PandasWrapped

convert(::Type{PyObject}, x::PandasWrapped) = x.pyo
PyCall.PyObject(x::PandasWrapped) = x.pyo

macro pytype(name, class)
    quote
        immutable $(name) <: PandasWrapped
            pyo::PyObject
            $(esc(name))(pyo::PyObject) = new(pyo)
            function $(esc(name))(args...; kwargs...)
                pandas_method = pandas_mod.$name
                new(pandas_method(args...; kwargs...))
            end
        end
        $(esc(:start))(x::$name) = start(x.pyo)
        function $(esc(:next))(x::$name, state)
            new_val, new_state = next(x.pyo, state)
            return pandas_wrap(new_val), new_state
        end
        $(esc(:done))(x::$name, state) = done(x.pyo, state)
        push!(type_map, ($class, $name))
    end
end

quot(x) = Expr(:quote, x)

function Base.Array(x::PandasWrapped)
    c = np[:asarray](x)
    if typeof(c).parameters[1] == PyObject
        out = cell(size(x))
        for row=1:size(x,1)
            for col=1:size(x,2)
                out[row, col] = convert(PyAny, c[row, col])
            end
        end
        out
    else
        c
    end
end

Base.values(x::PandasWrapped) = convert(PyArray, x.pyo["values"])

function pandas_wrap(pyo::PyObject)
    for (pyt, pyv) in type_map
        if pyisinstance(pyo, pyt)
            return pyv(pyo)
        end
    end
    return convert(PyAny, pyo)
end

pandas_wrap(x::Union{AbstractArray, Tuple}) = [pandas_wrap(_) for _ in x]

pandas_wrap(pyo) = pyo

#fix_arg(x::Range1) = pyeval(@sprintf "slice(%d, %d)" x.start x.start+x.len)
fix_arg(x::Range) = pyeval(@sprintf "slice(%d, %d, %d)" x.start x.start+x.len*x.step x.step)
fix_arg(x) = x

function pyattr(class, method, orig_method)
    if orig_method == :nothing
        m_quote = quot(method)
    else
        m_quote = quot(orig_method)
    end
    quote
        function $(esc(method))(pyt::$class, args...; kwargs...)
            pyo = pyt.pyo[string($m_quote)]
            if pytype_query(pyo) == Function
                new_args = [fix_arg(arg) for arg in args]
                pyo = pyt.pyo[$m_quote](new_args...; kwargs...)
            else
                pyo = pyt.pyo[string($m_quote)]
            end

            wrapped = pandas_wrap(pyo)
        end
    end
end

function delegate(new_func, orig_func, escape=false)
    @eval begin
        function $(escape ? esc(new_func) : new_func)(args...; kwargs...)
            f = $(orig_func)
            pyo = f(args...; kwargs...)
            pandas_wrap(pyo)
        end
    end
end

macro delegate(new_func, orig_func)
    delegate(new_func, orig_func, true)
end

macro pyattr(class, method, orig_method)
    pyattr(class, method, orig_method)
end

macro df_pyattrs(methods...)
    classes = [:DataFrame, :Series]
    m_exprs = Expr[]
    for method in methods
        exprs = Array(Expr, length(classes))
        for (i, class) in enumerate(classes)
            exprs[i] = pyattr(class, method, :nothing)
        end
        push!(m_exprs, Expr(:block, exprs...))
    end
    Expr(:block, m_exprs...)
end

macro df_pyattrs_eval(methods...)
    m_exprs = Expr[]
    for method in methods
        push!(m_exprs, quote
            function $(esc(method))(df::PandasWrapped, arg)
                res = pyeval(@sprintf("df.%s(\"%s\")", $method, arg), df=df)
                pandas_wrap(res)
            end
        end)
    end
    Expr(:block, m_exprs...)
end

macro gb_pyattrs(methods...)
    classes = [:GroupBy, :SeriesGroupBy]
    m_exprs = Expr[]
    for method in methods
        exprs = Array(Expr, length(classes))
        for (i, class) in enumerate(classes)
            exprs[i] = pyattr(class, method, :nothing)
        end
        push!(m_exprs, Expr(:block, exprs...))
    end
    Expr(:block, m_exprs...)
end

macro pyasvec(class, offset)
    offset = eval(offset)
    if offset
        index_expr = quote
            function $(esc(:getindex))(pyt::$class, args...)
                new_args = tuple([fix_arg(arg-1) for arg in args]...)
                pyo = pyt.pyo[:__getitem__](length(new_args)==1 ? new_args[1] : new_args)
                pandas_wrap(pyo)
            end

            function $(esc(:setindex!))(pyt::$class, value, idxs...)
                new_idx = [fix_arg(idx-1) for idx in idxs]
                pyt.pyo[:__setitem__](tuple(new_idx...), value)
            end
        end
    else
        index_expr = quote
            function $(esc(:getindex))(pyt::$class, args...)
                new_args = tuple([fix_arg(arg) for arg in args]...)
                pyo = pyt.pyo[:__getitem__](length(new_args)==1 ? new_args[1] : new_args)
                pandas_wrap(pyo)
            end

            function $(esc(:setindex!))(pyt::$class, value, idxs...)
                new_idx = [fix_arg(idx) for idx in idxs]
                if length(new_idx) > 1
                    pyt.pyo[:__setitem__](tuple(new_idx...), value)
                else
                    pyt.pyo[:__setitem__](new_idx[1], value)
                end
            end
        end
    end
    if class in [:Iloc, :Loc, :Ix]
        length_expr = quote
            function $(esc(:length))(x::$class)
                x.pyo[:obj][:__len__]()
            end
        end
    else
        length_expr = quote
            @pyattr $class length __len__
        end
    end
    quote
        $index_expr
        $length_expr
        function $(esc(:endof))(x::$class)
            length(x)
        end
    end
end

@pytype DataFrame pandas_mod.core[:frame]["DataFrame"]
@pytype Iloc pandas_mod.core[:indexing]["_iLocIndexer"]
@pytype Loc pandas_mod.core[:indexing]["_LocIndexer"]
@pytype Ix pandas_mod.core[:indexing]["_IXIndexer"]
@pytype Series pandas_mod.core[:series]["Series"]
@pytype MultiIndex pandas_mod.core[:index]["MultiIndex"]
@pytype Index pandas_mod.core[:index]["Index"]
@pytype GroupBy pandas_mod.core[:groupby]["DataFrameGroupBy"]
@pytype SeriesGroupBy pandas_mod.core[:groupby]["SeriesGroupBy"]

@pyattr DataFrame groupby nothing
@pyattr DataFrame columns nothing
@pyattr GroupBy app apply

@gb_pyattrs mean std agg aggregate median var ohlc transform groups indices get_group hist plot count

siz(gb::GroupBy) = gb.pyo[:size]()

@df_pyattrs iloc loc reset_index index head xs plot hist join align drop drop_duplicates duplicated filter first idxmax idxmin last reindex reindex_axis reindex_like rename tail set_index select take truncate abs any clip clip_lower clip_upper corr corrwith count cov cummax cummin cumprod cumsum describe diff mean median min mode pct_change rank quantile sum skew var std dropna fillna replace delevel pivot reodrer_levels sort sort_index sortlevel swaplevel stack unstack T boxplot argsort order asfreq asof shift first_valid_index last_valid_index weekday resample tz_conert tz_localize

@df_pyattrs_eval to_clipboard to_csv to_dense to_dict to_excel to_gbq to_hdf to_html to_json to_latex to_msgpack to_panel to_pickle to_records to_sparse to_sql to_string query

Base.size(x::Union{Loc, Iloc, Ix}) = x.pyo[:obj][:shape]
Base.size(df::PandasWrapped, i::Integer) = size(df)[i]
Base.size(df::PandasWrapped) = df.pyo[:shape]

@pyasvec Series false
@pyasvec Loc false
@pyasvec Ix false
@pyasvec Iloc true
@pyasvec DataFrame false
@pyasvec Index true
@pyasvec GroupBy false

Base.ndims(df::Union{DataFrame, Series}) = length(size(df))

for m in [:read_pickle, :read_csv, :read_html, :read_json, :read_excel, :read_table, :save, :stats,  :melt, :ewma, :concat, :merge, :pivot_table, :crosstab, :cut, :qcut, :get_dummies, :resample, :date_range, :to_datetime, :to_timedelta, :bdate_range, :period_range, :ewma, :ewmstd, :ewmvar, :ewmcorr, :ewmcov, :rolling_count, :expanding_count, :rolling_sum, :expanding_sum, :rolling_mean, :expanding_mean, :rolling_median, :expanding_median, :rolling_var, :expanding_var, :rolling_std, :expanding_std, :rolling_min, :expanding_min, :rolling_max, :expanding_max, :rolling_corr, :expanding_corr, :rolling_corr_pairwise, :expanding_corr_pairwise, :rolling_cov, :expanding_cov, :rolling_skew, :expanding_skew, :rolling_kurt, :expanding_kurt, :rolling_apply, :expanding_apply, :rolling_quantile, :expanding_quantile, :rolling_window]
    delegate(m, quote pandas_mod.$m end)
end

# for m in ["count", "sum", "mean", "median", "var", "std", "min", "max", "corr", "corr_pairwise", "cov", "skew", "kurt", "apply", "quantile"]
#     for f in ["rolling", "expanding"]
#         s = symbol(string(f, "_", m))
#         delegate(m, quote pandas_mod.$s end)
#     end
# end

function show(io::IO, df::PandasWrapped)
    s = df.pyo[:__str__]()
    println(io, s)
end

function query(df::DataFrame, e::Expr) # This whole method is a terrible hack
    s = string(e)
    s = s[3:end-1]
    s = replace(s, "&&", "&")
    s = replace(s, "||", "|")
    query(df, s)
end

for m in [:from_arrays, :from_tuples]
    @eval function $m(args...; kwargs...)
        f = pandas_raw["MultiIndex"][string($(quot(m)))]
        res = pycall(f, PyObject, args...; kwargs...)
        pandas_wrap(res)
    end
end

for op in [(:+, :__add__), (:*, :__mul__), (:/, :__div__), (:-, :__sub__)]
    @eval begin
        function $(op[1])(x::PandasWrapped, y::PandasWrapped)
            f = $(quot(op[2]))
            py_f = pyeval(string("x.", f), x=x.pyo)
            res = py_f(y)
            return pandas_wrap(res)
        end

        function $(op[1])(x::PandasWrapped, y)
            f = $(quot(op[2]))
            py_f = pyeval(string("x.", f), x=x.pyo)
            res = py_f(y)
            return pandas_wrap(res)
        end

        $(op[1])(y, x::PandasWrapped) = $(op[1])(x, y)
    end
end

for op in [(:-, :__neg__)]
    @eval begin
        $(op[1])(x::PandasWrapped) = pandas_wrap(x.pyo[$(quot(op[2]))]())
    end
end

function setcolumns!(df::PandasWrapped, new_columns)
    df.pyo[:__setattr__]("columns", new_columns)
end

function deletecolumn!(df::DataFrame, column)
    df.pyo[:__delitem__](column)
end

name(s::Series) = s.pyo[:name]
setname!(s::Series, name) = s.pyo[:name] = name

end


text:
module Retirement

using XLSX, DataFrames, ArgParse, Pkg, Dates, Logging

export run

"Global Constants"
const version=Pkg.project().version
const today=Dates.today()
const started=Dates.now()
const startDir=pwd()
const program=@__MODULE__

const monthPeriod=12
const weekPeriod=52
const period=monthPeriod
@enum accountStatus active inactive none

"Default Arguments"
const dataDir=pwd()*"/data"
const outDir=pwd()*"/output"
const dataFile="Retirement.xlsx"
const outFile="RetirementRun"*string(started)*".xlsx"

about() =  println("$program Version: $version $(string(started))")

"Command line arguments"
s=ArgParseSettings()
@add_arg_table! s begin
"-D"
  help = "Set debugging"
  action = :store_true
"--save"
  help = "path to save output"
  default = joinpath(outDir, outFile)
  required=false
"path"
  help = "path for input files"
  default = joinpath(dataDir, dataFile)
  required=false
end

logmsg(msg)=string(now(), " ", msg)
rows(df::DataFrame)=size(df, 1)
readxl(file, sheet)=DataFrame(XLSX.readtable(file, sheet)...)

"Growth Type "
abstract type Growth end
abstract type ConstGrowth <: Growth end
abstract type ZeroGrowth <: Growth end
abstract type FixedLoan <: Growth end

"Generic growth rule.  Empty"
grow(g::Growth, base, factors::NamedTuple)=error("Illegal call to abstract type Growth $g")

"Const Growth Rule"
grow(g::Type{ConstGrowth}, base, factors::NamedTuple,  period=monthPeriod)=base*(factors.rate/period + 1.0)

"Zero Growth Rule"
grow(g::Type{ZeroGrowth}, base::Float64, factors::NamedTuple)=  base

"""
Loan factor[1] = rate, factor[2] = payment loan = interest - payment
Returns state of the loan (active, closed)
"""
function grow(g::Type{FixedLoan}, base, factors::NamedTuple, period=monthPeriod)
  if (base <= factors.payment)
    closeLoan(factors.name)
    return 0.0
  end
  base*(1+factors.rate/period)-factors.payment
end

function calcInitial(sheet, calcType)
  @debug logmsg("$calcType rows: $(rows(sheet))")
  initialSum=sum(sheet.Initial)
  @debug logmsg("Initial $calcType: $initialSum")
  initialSum
end

"Save row with name / value from dataframe df into dictionary d"
pushRow!(d, df, row)=push!(d, df[row,1] => df[row,2])

"Create a dictionary from name value pairs in a dataframe"
function df2Dict!(df)
  d = Dict{String, Any}()
  map(r -> pushRow!(d, df, r), 1:nrow(df))
  d
end

"Set initial values from input config file"
function getInitialValues(dataFileName)
  income=readxl(dataFileName, "Income")
  @debug logmsg("Income rows: $(rows(income))")
  expenses=readxl(dataFileName, "Expenses")
  @debug logmsg("Expense rows: $(rows(expenses))")
  assets=readxl(dataFileName, "Assets")
  @debug logmsg("Asset rows: $(rows(assets))")
  liabilities=readxl(dataFileName, "Liabilities")
  @debug logmsg("Liability rows: $(rows(liabilities))")
  constantSheet=readxl(dataFileName, "Constants")
  @debug logmsg("Constant rows: $(rows(constantSheet))")
  @info logmsg("Completed reading input file")

  constants = df2Dict!(constantSheet)
  @debug logmsg("Constants: $constants")
  income, expenses, assets, liabilities, constants
end

"Convert Decimal Age to Year month"
function ageToYearMonth(age::Number)
  y=trunc(Int, age)
  m=trunc(Int, (age-y)*12)
  m,y
end

"Calculate the retirement date from birth date and age at retirement"
function getRetDate(retAgeYear, birthDate)
  (retMonth, retYear) = ageToYearMonth(retAgeYear)
  retDate=birthDate+Year(retYear)+Month(retMonth)
  retDate
end

"Calculate the number of periods in the model"
function calcPeriods(type, modelLength)
  periods  = if type == "Month"; 12 elseif type == "Week"; 52 end
  periods*modelLength
end

"""
Create a DataFrame for saving outputs from the model run
The dataframe is structured with a row for each instance in the model.
"""
function initializeDF(initialValues)
  output= DataFrame(Date=Date[], NetWorth=Float64[], Assets=Float64[], Liabilties=Float64[], NetIncome=Float64[], Income=Float64[], Expenses=Float64[] )
  push!(output, initialValues)
end

"Reshape a ledger table so that the Description Column becomes the column headings and the Initial column becomes the first row in the table"
function reshapeLedger(ledgerDF::DataFrame, desc::Symbol, initial::Symbol)
  df = select(ledgerDF, [desc, initial])
  transform!(df, desc => ByRow(Symbol) => desc)
  convertFloat(v)=convert(Float64,v)
  transform!(df, initial => ByRow(convertFloat) => initial)
  newDF=permutedims(df, 1)
  select!(newDF, Not(1)) # Drop the description column
  newDF
end

"Add a row to the ledger for a df to model growth during the next period"
function growNextLedger!(ledgerDF, factorMap)
  lastRow=ledgerDF[end,:]
  newRow=copy(lastRow)
  cols=names(ledgerDF)
  g=map(name -> factorMap[name],  cols)
  updated=map(n -> grow(g[n].rule, newRow[n], g[n]), 1:length(newRow))
  push!(ledgerDF, updated)
end

"If there are any closed loans, set the expense to zero "
function zeroClosedLoans(liability, expense)
  lastLia=liability[end,:]
  lastexp=expense[end,:]
  for n in names(expense)
    status = checkLoanStatus(n)
    if status==inactive
      lastexp[Symbol(n)] = 0.0
    end
  end
end

"Grow each entry in each ledger by the growth rule"
function growLedgers!(income, incMap, expense, expMap, asset, asstMap, liabilty, liaMap, periods)
  for p in 1:periods-1
    growNextLedger!(income, incMap)
    growNextLedger!(expense, expMap)
    growNextLedger!(asset, asstMap)
    growNextLedger!(liabilty, liaMap)
    zeroClosedLoans(liabilty, expense)
  end
end

"Update the ruleMap from the row of initial conditions sheet"
function mapGrowth!(factorMap, row, descr=:Description, rule=:Rule, factor1=:Factor1, factor2=:Factor2)
  f=eval(Symbol(row[rule])) # Use the description to find the function
  if row[rule] == "ConstGrowth"
    factors= (rule=f, rate=convert(Float64, row[factor1]))
  elseif row[rule] == "ZeroGrowth"
    factors=(rule=f, name=row[descr], rate=missing)
  elseif row[rule] == "FixedLoan"
    rate=convert(Float64, row[factor1])
    payment=convert(Float64, row[factor2])
    factors=(rule=f, rate=rate, payment=payment, name=row[descr])
  else
    error("Illegal growth rule: $(row[rule])")
  end
  push!(factorMap, row[descr] => factors)
end

"Manage loans. Loans are used to connect liabilties with payments"
const loanState=Dict{String, accountStatus}()
openLoan(name)=  push!(loanState, name => active)
closeLoan(name)= (loanState[name] = inactive)
function checkLoanStatus(name)
  if name in keys(loanState)
    return loanState[name]
  else
    return none
  end
end

"Open a loan for every liability"
openLoans(lMap)=map(l -> openLoan(l.name), values(lMap))

"Create a list of factors for each column in a ledger"
function mapGrowthFactors(df::DataFrame)
  factorMap=Dict{String, NamedTuple}()
  map(r -> mapGrowth!(factorMap, r), eachrow(df))
  factorMap
end

"Return a DF with the sums across all values for a ledger by period"
sumLedger(df, colName)=select(df, AsTable(:) =>  ByRow(sum) =>  colName)

"Return a DF with the net income by period"
function calcNet(incomeDF, expenseDF, net=:NetIncome)
  sumDF=hcat(incomeDF, expenseDF)
  transform!(sumDF, [1,2] => ByRow(-) => net)
  select!(sumDF, net)
  sumDF
end

function addNetToAssets(netIncome, assetsDF, cash, net)
  assetsPlusNet=hcat(netIncome, assetsDF)
  assetsNet=select(assetsPlusNet,  [cash, net] => ByRow(+) => :NetCash)
  temp=hcat(assetsDF, assetsNet)
  select(temp, Not(cash))
end

"Calculate loan payment amount"
payment(principal, rate, periods)=principal * (rate * (1+rate)^periods )/ ((1 + rate)^periods -1)
fv(pmt, rate, periods)=pmt*(1+rate)^periods
pv(pmt, rate, periods)=pmt/(1+rate)^periods

"Save the results to an excel file"
function saveRunExcel(incomeDF, expenseDF, assetsDF, liabilityDF, netIncomeDF, outputPath)
  XLSX.writetable(outputPath,
    Income=( collect(DataFrames.eachcol(incomeDF)), DataFrames.names(incomeDF) ),
    Expense=( collect(DataFrames.eachcol(expenseDF)), DataFrames.names(expenseDF)),
    Assets=( collect(DataFrames.eachcol(assetsDF)), DataFrames.names(assetsDF)),
    Liabilties=( collect(DataFrames.eachcol(liabilityDF)), DataFrames.names(liabilityDF)),
    Net=( collect(DataFrames.eachcol(netIncomeDF)), DataFrames.names(netIncomeDF)),
    )
end

"Main workflow method"
function run(ARGS)
  #Get the configuration
  args=parse_args(ARGS, s)
  about()
  @info logmsg("Started")
  debug=args["D"]
  outputPath=args["save"]

  debug ? ENV["JULIA_DEBUG"] = Main  : nothing
  @debug logmsg(args)
  @debug logmsg("Current directory: $startDir")
  dataFileName=args["path"]
  @info logmsg("Reading configuration from file $dataFileName")

  #Set initial values
  initialIncome, initialExpenses, initialAssets, initialLiabilities, constants =  getInitialValues(dataFileName)

  initIncomeSum=calcInitial(initialIncome, string(:Income))
  initExpensesSum=calcInitial(initialExpenses, string(:Expenses))
  initalNetIncome=initIncomeSum - initExpensesSum

  initAssetSum=calcInitial(initialAssets, string(:Assets))
  initLiabiltySum=calcInitial(initialLiabilities, string(:Liabilties))
  initialNetWoth=initAssetSum-initLiabiltySum

  retAgeYear=constants["Retirement Age"]
  birthDate=constants["Birth Date"]
  modelStart=constants["Model Start"]
  modelLength=constants["Model Length"]
  periodType=constants["Period Type"]
  retDate=getRetDate(retAgeYear, birthDate)
  @info logmsg("retAgeYear: $retAgeYear  birthDate $birthDate")
  @info logmsg("Planned retirement Date: $retDate")
  periods=calcPeriods(periodType, modelLength)
  @info logmsg("There are $periods periods in the model")
  initialValues= [modelStart initialNetWoth initAssetSum initLiabiltySum initalNetIncome initIncomeSum initExpensesSum]
  @debug "Initial Values $initialValues"
  output=initializeDF(initialValues)
  @debug "Output: $output"

  #Reshape tables with categories in columns and one row per period (Month)
  incomeDF=reshapeLedger(initialIncome, :Description, :Initial)
  expenseDF=reshapeLedger(initialExpenses, :Description, :Initial)
  assetsDF=reshapeLedger(initialAssets, :Description, :Initial)
  liabilityDF=reshapeLedger(initialLiabilities, :Description, :Initial)

  #Grow each category for the number of periods by the rule and factors associated with it
  incomeMap=mapGrowthFactors(initialIncome)
  @debug "incomemap: $incomeMap"
  expenseMap=mapGrowthFactors(initialExpenses)
  @debug "expenseMap: $expenseMap"
  assetsMap=mapGrowthFactors(initialAssets)
  @debug "assetMap: $assetMap"
  liabilityMap=mapGrowthFactors(initialLiabilities)
  @debug "liabilityMap: $liabilityMap"
  openLoans(liabilityMap)
  growLedgers!(incomeDF, incomeMap, expenseDF, expenseMap, assetsDF, assetsMap, liabilityDF, liabilityMap, periods)

  #Calculate totals
  totalLiabilty=sumLedger(liabilityDF, :TotalLiability)
  totalIncome=sumLedger(incomeDF, :TotalIncome)
  totalExpense=sumLedger(expenseDF, :TotalExpense)

  #Calculate Nets after expenses and income
  netIncome=calcNet(totalIncome, totalExpense)
  assetsNet=addNetToAssets(netIncome, assetsDF, :Cash, :NetIncome)
  totalAssets=sumLedger(assetsNet, :TotalAssets)
  netWorth=calcNet(totalAssets, totalLiabilty, :NetWorth)
  net=hcat(netIncome, totalIncome, totalExpense, netWorth)
  @info logmsg("Completed model run")

  #Save model run
  if !isdir(outDir)
    mkdir(outDir)
  end
  @info logmsg("Saving output to $outputPath")
  saveRunExcel(incomeDF, expenseDF, assetsDF, liabilityDF, net, outputPath)
end

#ENV["JULIA_DEBUG"]=all
run(ARGS)
end


text:
"""
    @jupyter expr default=nothing

Execute `expr` if it's currently in Jupyter.
"""
macro jupyter(expr, default=nothing)
    return esc(
        quote
            if isdefined(Main, :IJulia) && Main.IJulia.inited
                return $expr
            else
                return $default
            end
        end
    )
end

"""
    @script expr default=nothing

Execute `expr` if it's currently not in Jupyter.
"""
macro script(expr, default=nothing)
    return esc(
        quote
            if !(isdefined(Main, :IJulia) && Main.IJulia.inited)
                return $expr
            else
                return $default
            end
        end
    )
end

"""
    @tb expr default=nothing

Execute `expr` if the current logger is a `TBLogger`.
"""
macro tb(expr, default=nothing)
    return esc(
        quote
            !isdefined(Main, :current_logger) && import Logging
            !isdefined(Main, :TBLogger) && import TensorBoardLogger
            if Logging.current_logger() isa TensorBoardLogger.TBLogger
                return $expr
            else
                return $default
            end
        end
    )
end

"""
    @wb expr default=nothing

Execute `expr` if the current logger is a `WBLogger`.
"""
macro wb(expr, default=nothing)
    return esc(
        quote
            !isdefined(Main, :current_logger) && import Logging
            !isdefined(Main, :WBLogger) && import WeightsAndBiasLogger
            if Logging.current_logger() isa WeightsAndBiasLogger.WBLogger
                return $expr
            else
                return $default
            end
        end
    )
end

function checknumerics_strs(vcheck, vmonitor...; vcheckname=nothing, vmonitornames=["monitor$i" for i in 1:length(vmonitor)])
    for i in eachindex(vmonitor)
        @assert size(vcheck) == size(vmonitor[i]) "All variables in `vmonitor` should have the same size as `vcheck`; vmonitor[$i] doesn't."
    end
    # Get indics for NaN and Inf
    idcs = filter(eachindex(vcheck)) do i
        isnan(vcheck[i]) || isinf(vcheck[i])
    end
    # Convert to Cartesian
    c = CartesianIndices(size(vcheck))
    warnstr_list = []
    for i in idcs
        Istr = string(c[i].I)
        Istr_flat = replace(Istr, r"(\s)|(\()|(\))" => "")  # "(1, 2)" => "1,2"
        # Make checkstr
        if isnothing(vcheckname)
            warnstr = "Index $Istr is $(vcheck[i])"
        else
            warnstr = "$vcheckname[$Istr_flat] = $(vcheck[i])"
        end
        if length(vmonitor) > 0
            # Make monitorstr
            monitorstr = map(v -> v[i], vmonitor)
            monitorstr = collect(zip(map(mn -> "$mn[$Istr_flat]", vmonitornames), monitorstr))
            monitorstr = map(t -> "$(t[1]) = $(t[2])", monitorstr)
            monitorstr = join(monitorstr, "\n  ")
            warnstr = "$warnstr\n  $monitorstr"
        end
        push!(warnstr_list, warnstr)
    end
    return warnstr_list
end

"""
    @checknumerics scheck smonitor...

Check if each entry in `scheck` is `NaN` or `Inf`.
If so, report the index and the value of that index in all variables.

Usage:
```julia
vc = [1 Inf 3; NaN 5 6]
vm1 = [1 2 3; 3 4 5]
vm2 = [11 22 33; 44 55 66]
@checknumerics vc vm1 vm2
```
"""
macro checknumerics(scheck, smonitor...)
    local vcheck = esc(scheck)
    local vmonitor = esc.(smonitor)
    vcheckname = string(scheck)
    vmonitornames = string.(smonitor)
    return quote
        for warnstr in checknumerics_strs($vcheck, $(vmonitor...); vcheckname=$vcheckname, vmonitornames=$vmonitornames)
            @warn warnstr
        end
    end
end


text:
using Test
using NeXLParticle
using StaticArrays
using CoordinateTransformations
using Random

@testset "Translate" begin
    r = MersenneTwister(0xBADF00D)
    oldPts = [ SA[ 10.0 * rand(r, Float64,2)...] for _ in 1:20 ]
    θ, sc, off = π/3, SA[ 1.1 0.0; 0.0 0.9], SA[10.0,-10.0]
    tr = AffineMap( SA[ cos(θ) -sin(θ) ; sin(θ) cos(θ) ] * sc, off)
    newPts = (x->SA[ (tr(x) .+ 0.1*rand(r, 2))...]).(oldPts)
    # Estimate the transform that takes `oldPts` to `newPts`
    # tr and tr′ should be very similar
    tr′ = translate(oldPts, newPts)
    diff = tr′.(oldPts) .- newPts
    @test all(d->abs(d[1])<0.1, diff)
    @test all(d->abs(d[2])<0.1, diff)
    # Check that inverting the transform brings `newPts` back to `oldPts`
    itr′ = inv(tr′)
    diff = itr′.(newPts) .- oldPts
    @test all(d->abs(d[1])<0.1, diff)
    @test all(d->abs(d[2])<0.1, diff)

    zep = Zeppelin(joinpath(datadep"ZepTestArtifact","test.hdz"))
    zeptr = translate(zep, tr)
    # Check that `tr` and `tr′` both produce a similar forward transformation
    zeptr′ = translate(zep, tr′)
    @test all(rr->abs(rr[1].XABS-rr[2].XABS)<0.1, zip(eachrow(zeptr.data), eachrow(zeptr′.data)))
    @test all(rr->abs(rr[1].YABS-rr[2].YABS)<0.1, zip(eachrow(zeptr.data), eachrow(zeptr′.data)))

    # Check that `inv(tr′)` brings us back to close to `zep`
    izeptr′ = translate(zeptr, inv(tr′))
    @test all(rr->abs(rr[1].XABS-rr[2].XABS)<0.1, zip(eachrow(zep.data), eachrow(izeptr′.data)))
    @test all(rr->abs(rr[1].YABS-rr[2].YABS)<0.1, zip(eachrow(zep.data), eachrow(izeptr′.data)))
end

max_stars_repo_name: NicholasWMRitchie/NeXLParticle

max_issues_repo_name: NicholasWMRitchie/NeXLParticle

max_forks_repo_name: NicholasWMRitchie/NeXLParticle


text:
mutable struct Context
  cache::Union{IdDict{Any,Any},Nothing}
end

Context() = Context(nothing)

cache(cx::Context) = cx.cache == nothing ? (cx.cache = IdDict()) : cx.cache

struct J{S,T}
  t::T
end

J{S}(x) where S = J{S,typeof(x)}(x)

Base.show(io::IO, j::J{S}) where S = print(io, "J($(j.t[1]))")

# interface2.jl

# Wrappers

_forward(f, args...) = _forward(Context(), f, args...)

tailmemaybe(::Nothing) = nothing
tailmemaybe(x::Tuple) = Base.tail(x)

function forward(f, args...)
  y, back = _forward(f, args...)
  y, Δ -> tailmemaybe(back(Δ))
end

function gradient(f, args...)
  y, J = forward(f, args...)
  y isa Real || error("Function output is not scalar")
  return J(1)
end

derivative(f, x) = gradient(f, x)[1]

Base.adjoint(f::Function) = x -> derivative(f, x)

# Param-style wrappers

# TODO store ids only
struct Params
  params::IdSet{Any}
  Params(xs) = new(IdSet(xs))
end

@forward Params.params Base.iterate

struct Grads
  grads::IdDict{Any,Any}
end

Base.show(io::IO, ps::Grads) = print(io, "Grads(...)")

@forward Grads.grads Base.setindex!, Base.getindex, Base.haskey

function forward(f, ps::Params)
  cx = Context()
  y, back = _forward(cx, f)
  y, function (Δ)
    for p in ps
      cache(cx)[p] = nothing
    end
    back(Δ)
    Grads(cx.cache) # TODO make a copy
  end
end

# Reflection

# function code_grad(f, T)
#   forw = code_typed(_forward, Tuple{Context,Typeof(f),T.parameters...})[1]
#   Y, J = forw[2].parameters
#   back = typed_meta(Tuple{J,Y}, optimize=true)
#   back = back.code=>back.ret
#   (forw, back)
# end

# macro code_grad(ex)
#   isexpr(ex, :call) || error("@code_grad f(args...)")
#   f, args = ex.args[1], ex.args[2:end]
#   :(code_grad($(esc(f)), typesof($(esc.(args)...))))
# end


text:
###############
# H-OBCA: Hierarchical Optimization-based Collision Avoidance - a path planner for autonomous parking
# Copyright (C) 2018
# Atsushi SAKAI [atsushisakai@global.komatsu; Komatsu Ltd / MPC Lab]
# Alexander LINIGER [liniger@control.ee.ethz.ch; Automatic Control Lab, ETH Zurich]
# Xiaojing ZHANG [xiaojing.zhang@berkeley.edu; MPC Lab, UC Berkeley]
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
###############
# The paper describing the theory can be found here:
# 	X. Zhang, A. Liniger and F. Borrelli; "Optimization-Based Collision Avoidance"; Technical Report, 2017, [https://arxiv.org/abs/1711.03449]
#   X. Zhang, A. Liniger, A. Sakai and F. Borrelli; "Autonomous  Parking  using  Optimization-Based  Collision  Avoidance"; Technical Report, 2018 [add URL]
###############

module collision_check

using NearestNeighbors
using PyPlot

const B = 1.0 #[m] distance from rear to vehicle back end
const C = 3.7 #[m] distance from rear to vehicle front end
const I = 2.0 #[m] width of vehicle
const WBUBBLE_DIST = (B+C)/2.0-B #[m] distance from rear and the center of whole bubble
const WBUBBLE_R = (B+C)/2.0 #[m] whole bubble radius

const vrx = [C, C, -B, -B, C ]
const vry = [-I/2.0, I/2.0, I/2.0, -I/2.0, -I/2.0]

function check_collision(x, y, yaw, kdtree, ox, oy)

    for (ix, iy, iyaw) in zip(x, y, yaw)
        cx = ix + WBUBBLE_DIST*cos(iyaw)
        cy = iy + WBUBBLE_DIST*sin(iyaw)

        # Whole bubble check
        ids = inrange(kdtree, [cx, cy], WBUBBLE_R, true)
        if length(ids) == 0 continue end

        if !rect_check(ix, iy, iyaw, ox[ids], oy[ids])
            return false #collision
        end
    end

    return true #OK

end

function rect_check(ix, iy, iyaw, ox, oy)

    c = cos(-iyaw)
    s = sin(-iyaw)

    for (iox, ioy) in zip(ox, oy)
        tx = iox - ix
        ty = ioy - iy
        lx = (c*tx - s*ty)
        ly = (s*tx + c*ty)

        sumangle = 0.0
        for i in 1:length(vrx)-1
            x1 = vrx[i] - lx
            y1 = vry[i] - ly
            x2 = vrx[i+1] - lx
            y2 = vry[i+1] - ly
            d1 = hypot(x1,y1)
            d2 = hypot(x2,y2)
            theta1 = atan2(y1,x1)
            tty = (-sin(theta1)*x2 + cos(theta1)*y2)

            tmp = (x1*x2+y1*y2)/(d1*d2)
            if tmp >= 1.0 tmp = 1.0 end

            if tty >= 0.0
                sumangle += acos(tmp)
            else
                sumangle -= acos(tmp)
            end
        end

        if sumangle >= pi
            return false
        end
    end

    return true #OK
end

function main()

    ox = rand(3)*30.0 - 30.0/2.0
    oy = rand(3)*30.0 - 30.0/2.0

    kdtree = KDTree(hcat(ox, oy)')

    x = [10.0, 5.0]
    y = [10.0, 5.0]
    yaw = [deg2rad(10.0), deg2rad(0.0)]

    flag = check_collision(x, y, yaw, kdtree, ox, oy)
    if flag
        println("OK")
    else
        println("Collision")
    end

    plot(ox, oy, ".r")
    grid(true)
    axis("equal")
    show()

end

if length(PROGRAM_FILE)!=0 &&
    contains(@__FILE__, PROGRAM_FILE)

    @time main()
end

end #module

max_stars_repo_name: treyfortmuller/barc

max_issues_repo_name: treyfortmuller/barc


text:
# Exercise 1-1
# 1. Leaving out () shows function information. Missing parenthesis will result in incomplete expression or error.
# 2. Similar to 1.
# 3. Putting + before numbers will just print the value of the number. 2 ++ 2 is not valid, but 2 + + 2 is valid.
# 4. Leading zeroes will be removed when printed as numbered values.
# 5. Invalid expression, eg. 2 2 is not valid.

# Exercise 1-2
println("Number of seconds in 42min 42s = ", 42*60 + 42, 's')
println("Number of miles in 10km = ", 10 * 1.41, " miles")
println("For a 10km race in 37min and 48s,")
println("Average pace = ", (37 * 60 + 48) / 10, " s/mile")
println("or ", (37 + 48/60) / 10, " min/mile")
println("Average speed = ", 10 / (37 + 48/60) / 60, " miles/hr")


text:
# decompose formula into normal vs iv part
function is_iv(ex::Expr)
	if (ex.head == :macrocall && ex.args[1] == Symbol("@~")) || (ex.head == :call && ex.args[1] == :(~))
		# expression with ~
		length(ex.args) == 3 || error("malformed expression in formula")
		has_iv = true
		endos = ex.args[2]
		ivs = ex.args[3]
	else
		has_iv = false
		endos = nothing
		ivs = nothing
	end
	return has_iv, endos, ivs
end
is_iv(ex) = false, nothing, nothing

function decompose_iv!(rf::Formula)
	iv_formula = nothing
	iv_terms = nothing
	endo_formula = nothing
	endo_terms = nothing
	if !isa(rf.rhs, Expr)
		# symbol or Int 0
		# case with exactly 1 exogenous, 0 endogeneous
		has_iv = false
	else
		has_iv, endos, ivs = is_iv(rf.rhs)
		if has_iv
			# case without exogeneous variables
			has_iv = true
			endo_formula = @eval(@formula($nothing ~ $(endos)))
			iv_formula = @eval(@formula($nothing ~ $(ivs)))
			rf.rhs = :1
		elseif (rf.rhs.head == :call) && (rf.rhs.args[1] == :(+))
			# case with exogeneous variable(s)
			i = 2
			while !has_iv && i <= length(rf.rhs.args)
				has_iv, endos, ivs = is_iv(rf.rhs.args[i])
				if has_iv
					endo_formula = @eval(@formula($nothing ~ $(endos)))
					iv_formula = @eval(@formula($nothing ~  $(ivs)))
					splice!(rf.rhs.args, i)
				else
					i += 1
				end
			end
		end
	end
	if has_iv
		iv_terms = Terms(iv_formula)
		iv_terms.intercept = false
		endo_terms = Terms(endo_formula)
		endo_terms.intercept = false
	end
	return(has_iv, iv_formula, iv_terms, endo_formula, endo_terms)
end

function secondstage!(rf::Formula)
	(has_iv, iv_formula, iv_terms, endo_formula, endo_terms) = decompose_iv!(rf)
	if has_iv
	    if typeof(rf.rhs) == Symbol
	        rf.rhs = endo_formula.rhs
	    else
	        push!(rf.rhs.args, endo_formula.rhs)
	    end
	end
end

##############################################################################
##
## build model
##
##############################################################################
function nonmissing(mf::ModelFrame)
	if  :msng ∈ fieldnames(typeof(mf))
		mf.msng
	else
		mf.nonmissing
	end
end

# copy ModelFrame but modify to have esample as argument
function ModelFrame2(trms::Terms, d::AbstractDataFrame, esample; contrasts::Dict = Dict())
	df = DataFrame(map(x -> d[x], trms.eterms), Symbol.(string.(trms.eterms)))
	if !all(esample)
		# doing view takes more time and memory for some reason (both in in ModelFrame2 and in evalConstrasts)
		# df = view(df, esample, names(df))
		df = df[esample, :]
	end
	evaledContrasts = evalcontrasts(df, contrasts)
	## Check for non-redundant terms, modifying terms in place
	check_non_redundancy!(trms, df)
	ModelFrame(df, trms, esample, evaledContrasts)
end

# Directly from DataFrames.jl
function dropresponse(trms::Terms)
    if trms.response
        ckeep = 2:size(trms.factors, 2)
        rkeep = vec(any(trms.factors[:, ckeep], dims = 2))
        Terms(trms.terms, trms.eterms[rkeep], trms.factors[rkeep, ckeep],
              trms.is_non_redundant[rkeep, ckeep], trms.order[ckeep], false, trms.intercept)
    else
        trms
    end
end

function allvars(ex::Expr)
    if ex.head != :call error("Non-call expression encountered") end
    vcat([allvars(a) for a in ex.args[2:end]]...)
end
allvars(f::Formula) = unique(vcat(allvars(f.rhs), allvars(f.lhs)))
allvars(sym::Symbol) = [sym]
allvars(::Any) = Array{Symbol}(undef, 0)


text:
using JLD2, Distributions, Plots, StatPlots, AverageShiftedHistograms, DynNets
 plotly()

 N_est = 5
 NGWpoints = 3
 Tvals =  [  600 400 200 ]#[50 250 1250]#
 N_Tind = length(Tvals)[1]
 scoreRescType = ""#"FISHER-EWMA"#
 useStartVal = false# true# false
 Nmin = 10
 Nspacing =  10
 Nmax = 70
Nvals = Vector(Nmax:-Nspacing:Nmin)
 N_Nind = length(Nvals)

 save_fold = "./data/estimatesTest/compComplexity/"
 file_name = "EstTest_varN_Nest_$(N_est)_$(Nmin)_$(Nspacing)_$(Nmax)_varGW_$(NGWpoints)_T_$(Tvals)_corectSpec"* scoreRescType *".jld"
 save_path = save_fold*file_name#
 @load(save_path,est_conv_flag,est_times,est_flag,Nvals,GWvalStore ,useStartVal)

## histograms of static parameters

function Bias(A::Array{<:Real,2},B::Array{<:Real,2})
    return  Bias =    (A.-B)./std(A .-B)
end
Bias(A::SharedArray{<:Real,2},B::SharedArray{<:Real,2}) = BiasRel( convert(Array{Real},A),convert(Array{Real},B))

pW1 = Plots.Plot{Plots.PlotlyBackend};
pB1 = Plots.Plot{Plots.PlotlyBackend};
pA1 = Plots.Plot{Plots.PlotlyBackend};
ind = 0;
for ind = 1:length(Tvals)
#ind=1
errA =  Bias(simPar_3[ind,:,:],estPar_3[ind,:,:])
errB =   Bias(simPar_2[ind,:,:],estPar_2[ind,:,:])
errW =  (simPar_1[ind,:,:].- estPar_1[ind,:,:])./(repmat(std(estPar_1[ind,:,:],1),N_est,1))

rng_space = -5:.1:5

err= errW
tmpAsh = ash(err;rng =rng_space )
pW1 = plot(tmpAsh,
    title = scoreRescType *"  T= $(Tvals[ind]),
    W     ")
plot!(Normal(0,1),linewidth = 5)

err= errB
tmpAsh = ash(err;rng = rng_space)
pB1 = plot(tmpAsh,
title =    " $(N_est) Est   B   ")
plot!(Normal(0,1),linewidth = 5)

err= errA
rng_space = -3:.1:3
tmpAsh = ash(err;rng = rng_space)
pA1 = plot(tmpAsh,
title = "StartVal:" * startValFlag * "  GW: $(GWest), GBA: $(GBAest), N: $(N),   A   ")
plot!(Normal(0,1),linewidth = 5)
exp = "pl_T_$(ind)=  plot(pW1,pB1,pA1,layout=(3,1) )"
eval(parse(exp))
end

exp2 = " " ; for i = 1:ind exp2 *=  "pl_T_$(i) , " end

exp3 = "plot( " * exp2 * " layout = (1,ind),size=(1350,600),legend =:none)"
eval(parse(exp3))
#plot(pl_T_1 , pl_T_2,pl_T_3 , pl_T_3, layout = (1,4),size=(1350,600),legend =:none)

##

RavgMSE(A,B) = sqrt.(mean((A .- B).^2,1))
AvgBias(A,B) = mean((A .- B),1)
tmpModelEst = DynNets.SdErgm1(zeros(10,10),[ones(GWest),ones(GBAest),ones(GBAest)],groupsIndsEst,scoreRescType)
tmpModelDgp = DynNets.SdErgm1(zeros(10,10),[ones(GWdgp),ones(GBAdgp),ones(GBAdgp)],groupsIndsDgp,scoreRescType)

rmseTotGas = zeros(N_est,N,N_ind)
biasTotGas = zeros(N_est,N,N_ind)
rmseTotSnap = zeros(N_est,N,N_ind)
biasTotSnap = zeros(N_est,N,N_ind)
dgpTvParT = Array{Array{Float64,3},1}(N_ind)
filTvParT = Array{Array{Float64,3},1}(N_ind)
ParSnapT =  Array{Array{Float64,3},1}(N_ind)
# Compute MSE Bias and Snap estimates
@time for indT = 1:N_ind
    #indT=1
            T =Tvals[indT]
            tmpParDgp = zeros(T,N,N_est)
            tmpParFil = zeros(T,N,N_est)
            tmpParSnap = zeros(T,N,N_est)
        for n=1:N_est
    #n=100
    obsT = allObs[:,1:T,indT,n]'
    vecParEst =  [estPar_1[indT,n,:]; estPar_2[indT,n];  estPar_3[indT,n]]
    vecParDgp =  [simPar_1[indT,n,:]; simPar_2[indT,n];  simPar_3[indT,n]]
    GroupIndsEst = Array{Array{<:Int,1},1}(2)
    GroupIndsEst[1] = groupsIndsEst[1]; GroupIndsEst[2] = groupsIndsEst[2]
    tmpParFil[:,:,n] =   DynNets.score_driven_filter_or_dgpAndLikeliood(tmpModelEst,vecParEst,
                                                degsT = obsT,
                                                groupsInds = GroupIndsEst )[1]

    fitSnap = DynNets.estSnapSeq(DynNets.SdErgm1(obsT)) # sequence of snapshots estimate
    # remove the infinites
    for i=1:5
        infInd = find(.!isfinite.(fitSnap))

        replInd = infInd[infInd .> 2]
        fitSnap[replInd] = fitSnap[replInd.-1]
        infInd = find(.!isfinite.(fitSnap))
        replInd = infInd[infInd .< (Int(N*T) -2) ]
        fitSnap[replInd] = fitSnap[replInd.+1]
    end
    tmpParSnap[:,:,n] = fitSnap

    GroupIndsDgp = Array{Array{<:Int,1},1}(2)
    GroupIndsDgp[1] = groupsIndsDgp[1]; GroupIndsDgp[2] = groupsIndsDgp[2]
    tmpParDgp[:,:,n] = DynNets.score_driven_filter_or_dgpAndLikeliood(tmpModelDgp,vecParDgp,
                                                degsT = obsT,
                                                groupsInds = GroupIndsDgp )[1]

            rmseTotGas[n,:,indT] = RavgMSE(tmpParFil[:,:,n],tmpParDgp[:,:,n] )
            biasTotGas[n,:,indT] = AvgBias(tmpParFil[:,:,n],tmpParDgp[:,:,n] )
            rmseTotSnap[n,:,indT] = RavgMSE(fitSnap,tmpParDgp[:,:,n] )
            biasTotSnap[n,:,indT] = AvgBias(fitSnap,tmpParDgp[:,:,n] )
        end
    dgpTvParT[indT]  = tmpParDgp
    filTvParT[indT]  = tmpParFil
    ParSnapT[indT] = tmpParSnap
end

## Example single realization plots
Nlines=10
indPl = round(Int,linspace(1,N,Nlines))
C(g::ColorGradient) = RGB[g[z] for z=linspace(0,1,Nlines)]
lineColors = C(cgrad(:rainbow))

indT = 2
n=4
T = Tvals[indT]

uncMeansFit = simPar_1[indT,n,:]./(1-simPar_2[indT,n,:])
fit_T = [uncMeansFit'; dgpTvParT[indT][:,:,n] ]
#Compute the skewness of each degree's distribution
probMat =  1./(1+ Base.exp(uncMeansFit .+ uncMeansFit' ));DynNets.putZeroDiag!(probMat );
    degVar = sum(probMat.*(1-probMat),2);
    degSkew = round(degVar.^(-3/2).*sum( probMat.*(1-probMat).*(1-2.*probMat),2),3)
uncMeansFilFit = estPar_1[indT,n,groupsIndsEst[1]]./(1-estPar_2[indT,n,:])
fil_fit_T = [uncMeansFilFit'; filTvParT[indT][:,:,n] ]
fitSnap = [zeros(1,N); ParSnapT[indT][:,:,n]]
 p=plot()
for i = 1:Nlines p = plot!(fit_T[:,indPl[i]],  color = lineColors[i] , lab = "$(degSkew[i])" ) end #lineColors
for i = 1:Nlines p = plot!(fil_fit_T[:,indPl[i]],linestyle=:dash,  color = lineColors[i]  ) end #lineColors
for i = 1:Nlines p= plot!(fitSnap[:,indPl[i]],linestyle=:dot,markershape=:circle,
    markersize=1,linewidth = 0, markerstrokecolor=  lineColors[i], color = lineColors[i]  ) end

p =plot(p,titlefont = font(16),xlims = (0,50) ,titleloc=:center,
    title = "Gas DGP,N $(N) T = $(T) B = $(simPar_2[indT,n,1]), A = $(simPar_3[indT,n,1])",size=(1350,600))
##
p = plot()
 for i = 1:Nlines p = plot!(ObsAllT[indT][:,indPl[i]],  color = lineColors[i] ) end #lineColors
 plot(p)

## Box Plots
yLim = (0,1)
    pRT1 = boxplot( rmseTotGas[:,:,1]  , ylims = yLim  , legend = false, title = "T = $(Tvals[1]), N = $(N), GWdgp =  $(GWdgp), GWest = $(GWest)")
    pRT2 = boxplot( rmseTotGas[:,:,2]  , ylims = yLim  , legend = false, title = "T = $(Tvals[2])")
      yLim = (-0.8,1.1)
    pBT1 = boxplot( biasTotGas[:,:,1] , ylims = yLim , legend = false, title = "T = $(Tvals[1])")
    pBT2 = boxplot( biasTotGas[:,:,2] , ylims = yLim , legend = false, title = "T = $(Tvals[2])")
     ptot = plot(pRT1,pRT2, pBT1,pBT2, titlefont = font(12) ,titleloc=:center,size=(1350,600) )
##
yLim = (0,1.5)
    pRT1 = boxplot( rmseTotSnap[:,:,1]  , ylims = yLim  , legend = false, title = "T = $(Tvals[1]), N = $(N), GWdgp =  $(GWdgp), Single Snapshot ")
    pRT2 = boxplot( rmseTotSnap[:,:,2]  , ylims = yLim  , legend = false, title = "T = $(Tvals[2])")
          yLim = (-1.4,0.8)
    pBT1 = boxplot( biasTotSnap[:,:,1] , ylims = yLim , legend = false, title = "T = $(Tvals[1])")
    pBT2 = boxplot( biasTotSnap[:,:,2] , ylims = yLim , legend = false, title = "T = $(Tvals[2])")
     ptot = plot(pRT1,pRT2 ,pBT1,pBT2, titlefont = font(12) ,titleloc=:center,size=(1350,600) )

##

##


text:
using SolveDSGE

filename = "model3b.txt"
path3 = joinpath(@__DIR__,filename)
process_model(path3)

processed_filename = "model3b_processed.txt"
processed_path3 =  joinpath(@__DIR__,processed_filename)

dsge3 = retrieve_processed_model(processed_path3)

x3 = [0.05, 3.05, 0.7, 0.7, 2.8]

tol = 1e-8
maxiters = 1000

ss3_obj = compute_steady_state(dsge3,x3,tol,maxiters)
ss3 = ss3_obj.zero

N   = PerturbationScheme(ss3,1.0,"first")
NN  = PerturbationScheme(ss3,1.0,"second")
NNN = PerturbationScheme(ss3,1.0,"third")

soln_fo3 = solve_model(dsge3,N)
soln_so3 = solve_model(dsge3,NN)
soln_to3 = solve_model(dsge3,NNN)

P = ChebyshevSchemeDet(ss3,chebyshev_nodes,[11,11,11],4,[0.07 3.2 0.8; -0.07 2.9 0.6],tol,1e-6,maxiters)
PP = ChebyshevSchemeDet(ss3,chebyshev_nodes,[21,21,21],4,[0.07 3.2 0.8; -0.07 2.9 0.6],tol,1e-6,maxiters)
PPP = ChebyshevSchemeDet(ss3,chebyshev_nodes,[31,31,31],5,[0.07 3.2 0.8; -0.07 2.9 0.6],tol,1e-6,maxiters)

soln_nl3a = solve_model(dsge3,P)
soln_nl3b = solve_model(dsge3,soln_so3,PP)
soln_nl3c = solve_model(dsge3,soln_nl3b,PPP)

L = SmolyakSchemeDet(ss3,chebyshev_gauss_lobatto,3,[0.07 3.2 0.8; -0.07 2.9 0.6],tol,1e-6,maxiters)
LL = SmolyakSchemeDet(ss3,chebyshev_gauss_lobatto,3,[0.07 3.2 0.8; -0.07 2.9 0.6],tol,1e-6,maxiters)
LLL = SmolyakSchemeDet(ss3,clenshaw_curtis_equidistant,3,[0.07 3.2 0.8; -0.07 2.9 0.6],tol,1e-6,maxiters)

soln_nl3d = solve_model(dsge3,L)
soln_nl3e = solve_model(dsge3,soln_to3,LL)
soln_nl3f = solve_model(dsge3,soln_nl3a,LLL)

M = PiecewiseLinearSchemeDet(ss3,[11,11,11],[0.07 3.2 0.8; -0.07 2.9 0.6],tol,1e-6,maxiters)

soln_nl3h = solve_model(dsge3,M)
soln_nl3i = solve_model(dsge3,soln_to3,M)
soln_nl3j = solve_model(dsge3,soln_nl3c,M)

simulated_data3 = simulate(soln_nl3f,ss3[1:3],100)


text:
export @spike_str
macro spike_str(spike)
  spike=split(spike,"\n")
  if length(spike)==2
    threshold=spike[1]
    refractory=spike[2]
  else
    threshold=spike[1]
    refractory=0
  end
  ###############################################
  #initialize tSpike && construct function spike
  #which include spike,reset,refractory
  str="net.tSpike=Array(Vector{Float64},net.N)\n"
  str*="fill!(net.tSpike,[-$(refractory)])\n"
  print(str)
  map(x->eval(parse(x)),split(str,"\n"))

  threshold=replace(threshold,r"(.+)(?=[\>\<≂])", s"net.\1")

  str="function spike!(net::Network$timestr)\n"
  str*="isExceedThreshold=$(threshold)\n"
  if refractory!=0
    str*="isNotRefractory= (net.t).>=(map(x->x[end],net.tSpike)+$(refractory))\n"
    str*="isSpike=isExceedThreshold & isNotRefractory\n"
  else
    str*="isSpike=isExceedThreshold\n"
  end
  str*="net.i=find(isSpike)\n"
  str*="for spikeNeuron in net.i\n"
  str*="net.tSpike[spikeNeuron]=[net.tSpike[spikeNeuron];net.t]\n"
  str*="end\n"
  #   str*="map(x->push!(x,net.t),net.tSpike[net.i])\n"
  str*="end\n"
  print(str)
  eval(parse(str))
end


text:
@testset "stdlib" begin
    orig_logger = Base.CoreLogging.global_logger()

    try
        logger = getlogger("Memento")
        @test_log(
            logger,
            "notice",
            "Substituting global logging",
            Memento.config!("info"; substitute=true)
        )
    finally
        Base.CoreLogging.global_logger(orig_logger)
    end
end


text:
# Copyright (c) 2019 Yury Dvorkin, Miles Lubin, and contributors                #src
#                                                                               #src
# This work is licensed under a Creative Commons Attribution-ShareAlike 4.0     #src
# International License. URL: http://creativecommons.org/licenses/by-sa/4.0.    #src

# # Power Systems

# This tutorial demonstrates how to formulate basic power systems engineering
# models in JuMP.

# We will consider basic "economic dispatch" and "unit commitment" models
# without taking into account transmission constraints.

# ## Illustrative example

# For this example, we set the following characteristics of generators,
# transmission lines, wind farms and demands:

# | Quantity         | Generator 1 | Generator 2 |
# |:-----------------|:-----------:|------------:|
# | $g_{min}$, MW    | 0           | 300         |
# | $g_{max}$, MW    | 1000        | 1000        |
# | $c^g$, \$/MWh    | 50          | 100         |
# | $c^{g0}$, \$/MWh | 1000        | 0           |

# | Quantity      | Line 1 | Line 2 |
# |:--------------|:------:|-------:|
# | $f^{max}$, MW | 100    | 1000   |
# | $x$, p.u.     | 0.001  | 0.001  |

# | Quantity        | Wind farm 1 | Wind farm 2 |
# |:----------------|:-----------:|------------:|
# | $w^{f}$, MW     | 150         | 50          |
# | $c^{w}$, \$/MWh | 50          | 50          |

# | Quantity | Bus 1 | Bus 2 | Bus 3 |
# |:---------|:-----:|:-----:|------:|
# | $d$, MW  | 0     | 0     | 15000 |

# ## Economic dispatch

# Economic dispatch (ED) is an optimization problem that minimizes the cost of
# supplying energy demand subject to operational constraints on power system
# assets. In its simplest modification, ED is an LP problem solved for an
# aggregated load and wind forecast and for a single infinitesimal moment.

# Mathematically, the ED problem can be written as follows:

# ```math
# \min \sum_{i \in I} c^g_{i} \cdot g_{i} + c^w \cdot w,
# ```

# where $c_{i}$ and $g_{i}$ are the incremental cost (\\$/MWh) and power
# output (MW) of the $i^{th}$ generator, respectively, and $c^w$ and $w$ are the
# incremental cost (\\$/MWh) and wind power injection (MW), respectively.

# Subject to the constraints:

# * Minimum ($g^{\min}$) and maximum ($g^{\max}$) limits on power outputs of
#   generators:
#   $g^{\min}_{i} \leq g_{i} \leq g^{\max}_{i}.$
# * Constraint on the wind power injection:
#   $0 \leq w \leq w^f,$
#   where $w$ and $w^f$ are the wind power injection and wind power forecast,
#   respectively.
# * Power balance constraint:
#   $\sum_{i \in I} g_{i} + w = d^f,$
#   where $d^f$ is the demand forecast.

# Further reading on ED models can be found in A. J. Wood, B. F. Wollenberg, and
# G. B. Sheblé, "Power Generation, Operation and Control", Wiley, 2013.

# ## JuMP Implementation of Economic Dispatch

using JuMP
import DataFrames
import GLPK
using Plots
using StatsPlots

# Define some input data about the test system.

# We define some thermal generators:

function ThermalGenerator(
    min::Float64,
    max::Float64,
    fixed_cost::Float64,
    variable_cost::Float64,
)
    return (
        min = min,
        max = max,
        fixed_cost = fixed_cost,
        variable_cost = variable_cost,
    )
end

generators = [
    ThermalGenerator(0.0, 1000.0, 1000.0, 50.0),
    ThermalGenerator(300.0, 1000.0, 0.0, 100.0),
]

# A wind generator

WindGenerator(variable_cost::Float64) = (variable_cost = variable_cost,)

wind_generator = WindGenerator(50.0)

# And a scenario

function Scenario(
    demand::Float64,
    wind::Float64,
)
    return (demand = demand, wind = wind)
end

scenario = Scenario(1500.0, 200.0)

# Create a function solve_ed, which solves the economic dispatch problem for a
# given set of input parameters.

function solve_ed(generators::Vector, wind, scenario)
    ## Define the economic dispatch (ED) model
    ed = Model(GLPK.Optimizer)
    ## Define decision variables
    ## power output of generators
    N = length(generators)
    @variable(ed, generators[i].min <= g[i = 1:N] <= generators[i].max)
    ## wind power injection
    @variable(ed, 0 <= w <= scenario.wind)
    ## Define the objective function
    @objective(
        ed,
        Min,
        sum(generators[i].variable_cost * g[i] for i in 1:N) +
        wind.variable_cost * w,
    )
    ## Define the power balance constraint
    @constraint(ed, sum(g[i] for i in 1:N) + w == scenario.demand)
    ## Solve statement
    optimize!(ed)
    ## return the optimal value of the objective function and its minimizers
    return (
        g = value.(g),
        w = value(w),
        wind_spill = scenario.wind - value(w),
        total_cost = objective_value(ed),
    )
end

# Solve the economic dispatch problem

solution = solve_ed(generators, wind_generator, scenario);

println("Dispatch of Generators: ", solution.g, " MW")
println("Dispatch of Wind: ", solution.w, " MW")
println("Wind spillage: ", solution.wind_spill, " MW")
println("\n")
println("Total cost: \$", solution.total_cost)

# ### Economic dispatch with adjustable incremental costs

# In the following exercise we adjust the incremental cost of generator G1 and
# observe its impact on the total cost.

function scale_generator_cost(g, scale)
    return ThermalGenerator(
        g.min,
        g.max,
        g.fixed_cost,
        scale * g.variable_cost,
    )
end

start = time()
c_g_scale_df = DataFrames.DataFrame(
    ## Scale factor
    scale = Float64[],
    ## Dispatch of Generator 1 [MW]
    dispatch_G1 = Float64[],
    ## Dispatch of Generator 2 [MW]
    dispatch_G2 = Float64[],
    ## Dispatch of Wind [MW]
    dispatch_wind = Float64[],
    ## Spillage of Wind [MW]
    spillage_wind = Float64[],
    ## Total cost [$]
    total_cost = Float64[],
)
for c_g1_scale in 0.5:0.1:3.0
    ## Update the incremental cost of the first generator at every iteration.
    new_generators = scale_generator_cost.(generators, [c_g1_scale, 1.0])
    ## Solve the ed problem with the updated incremental cost
    sol = solve_ed(new_generators, wind_generator, scenario)
    push!(
        c_g_scale_df,
        (c_g1_scale, sol.g[1], sol.g[2], sol.w, sol.wind_spill, sol.total_cost),
    )
end
print(string("elapsed time: ", time() - start, " seconds"))

#-

c_g_scale_df

# ## Modifying the JuMP model in place

# Note that in the previous exercise we entirely rebuilt the optimization model
# at every iteration of the internal loop, which incurs an additional
# computational burden. This burden can be alleviated if instead of re-building
# the entire model, we modify a specific constraint(s) or the objective
# function, as it shown in the example below.

# Compare the computing time in case of the above and below models.

function solve_ed_inplace(
    generators::Vector,
    wind,
    scenario,
    scale::AbstractVector{Float64},
)
    start = time()
    obj_out = Float64[]
    w_out = Float64[]
    g1_out = Float64[]
    g2_out = Float64[]
    ## This function only works for two generators
    @assert length(generators) == 2
    ed = Model(GLPK.Optimizer)
    N = length(generators)
    @variable(ed, generators[i].min <= g[i = 1:N] <= generators[i].max)
    @variable(ed, 0 <= w <= scenario.wind)
    @objective(
        ed,
        Min,
        sum(generators[i].variable_cost * g[i] for i in 1:N) +
        wind.variable_cost * w,
    )
    @constraint(ed, sum(g[i] for i in 1:N) + w == scenario.demand)
    for c_g1_scale in scale
        @objective(
            ed,
            Min,
            c_g1_scale * generators[1].variable_cost * g[1] +
            generators[2].variable_cost * g[2] +
            wind.variable_cost * w,
        )
        optimize!(ed)
        push!(obj_out, objective_value(ed))
        push!(w_out, value(w))
        push!(g1_out, value(g[1]))
        push!(g2_out, value(g[2]))
    end
    df = DataFrames.DataFrame(
        scale = scale,
        dispatch_G1 = g1_out,
        dispatch_G2 = g2_out,
        dispatch_wind = w_out,
        spillage_wind = scenario.wind .- w_out,
        total_cost = obj_out,
    )
    elapsed = time() - start
    print(string("elapsed time: ", elapsed, " seconds"))
    return df
end

solve_ed_inplace(generators, wind_generator, scenario, 0.5:0.1:3.0)

# Adjusting specific constraints and/or the objective function is faster than
# re-building the entire model.

# ## A few practical limitations of the economic dispatch model

# ### Inefficient usage of wind generators

# The economic dispatch problem does not perform commitment decisions and, thus,
# assumes that all generators must be dispatched at least at their minimum power
# output limit. This approach is not cost efficient and may lead to absurd
# decisions. For example, if $ d = \sum_{i \in I} g^{\min}_{i}$, the wind power
# injection must be zero, i.e. all available wind generation is spilled, to meet
# the minimum power output constraints on generators.

# In the following example, we adjust the total demand and observed how it
# affects wind spillage.

demand_scale_df = DataFrames.DataFrame(
    demand = Float64[],
    dispatch_G1 = Float64[],
    dispatch_G2 = Float64[],
    dispatch_wind = Float64[],
    spillage_wind = Float64[],
    total_cost = Float64[],
)

function scale_demand(scenario, scale)
    return Scenario(scale * scenario.demand, scenario.wind)
end

for demand_scale in 0.2:0.1:1.5
    new_scenario = scale_demand(scenario, demand_scale)
    sol = solve_ed(generators, wind_generator, new_scenario)
    push!(
        demand_scale_df,
        (
            new_scenario.demand,
            sol.g[1],
            sol.g[2],
            sol.w,
            sol.wind_spill,
            sol.total_cost,
        ),
    )
end

demand_scale_df

#-

dispatch_plot = @df(
    demand_scale_df,
    Plots.plot(
        :demand,
        [:dispatch_G1, :dispatch_G2],
        labels = ["G1" "G2"],
        title = "Thermal Dispatch",
        legend = :bottomright,
        linewidth = 3,
        xlabel = "Demand",
        ylabel = "Dispatch [MW]"
    ),
)

wind_plot = @df(
    demand_scale_df,
    Plots.plot(
        :demand,
        [:dispatch_wind, :spillage_wind],
        labels = ["Dispatch" "Spillage"],
        title = "Wind",
        legend = :bottomright,
        linewidth = 3,
        xlabel = "Demand [MW]",
        ylabel = "Energy [MW]"
    ),
)

Plots.plot(dispatch_plot, wind_plot)

# This particular drawback can be overcome by introducing binary decisions on
# the "on/off" status of generators. This model is called unit commitment and
# considered later in these notes.

# For further reading on the interplay between wind generation and the minimum
# power output constraints of generators, we refer interested readers to R.
# Baldick, "Wind and Energy Markets: A Case Study of Texas," IEEE Systems
# Journal, vol. 6, pp. 27-34, 2012.

# ### Transmission-infeasible solution

# The ED solution is entirely market-based and disrespects limitations of the
# transmission network. Indeed, the flows in transmission lines would attain the
# following values:

# ```math
# f_{1-2} = 150 MW \leq f_{1-2}^{\max} = 100 MW
# ```

# ```math
# f_{2-3} = 1200 MW \leq f_{2-3}^{\max} = 1000 MW
# ```

# Thus, if this ED solution was enforced in practice, the power flow limits on
# both lines would be violated. Therefore, in the following section we consider
# the optimal power flow model, which amends the ED model with network
# constraints.

# The importance of the transmission-aware decisions is emphasized in E.
# Lannoye, D. Flynn, and M. O'Malley, "Transmission, Variable Generation, and
# Power System Flexibility," IEEE Transactions on Power Systems, vol. 30, pp.
# 57-66, 2015.

# ## Unit Commitment model

# The Unit Commitment (UC) model can be obtained from ED model by introducing
# binary variable associated with each generator. This binary variable can
# attain two values: if it is "1", the generator is synchronized and, thus, can
# be dispatched, otherwise, i.e. if the binary variable is "0", that generator
# is not synchronized and its power output is set to 0.

# To obtain the mathematical formulation of the UC model, we will modify the
# constraints of the ED model as follows:

# ```math
# g^{\min}_{i} \cdot u_{t,i} \leq g_{i} \leq g^{\max}_{i} \cdot u_{t,i},
# ```
# where $u_{i} \in \{0,1\}.$ In this constraint, if $u_{i} = 0$, then
# $g_{i}  = 0$. On the other hand, if $u_{i} = 1$, then
# $g^{min}_{i} \leq g_{i} \leq g^{max}_{i}$.

# For further reading on the UC problem we refer interested readers to G.
# Morales-Espana, J. M. Latorre, and A. Ramos, "Tight and Compact MILP
# Formulation for the Thermal Unit Commitment Problem," IEEE Transactions on
# Power Systems, vol. 28, pp. 4897-4908, 2013.

# In the following example we convert the ED model explained above to the UC
# model.

function solve_uc(generators::Vector, wind, scenario)
    uc = Model(GLPK.Optimizer)
    N = length(generators)
    @variable(uc, generators[i].min <= g[i = 1:N] <= generators[i].max)
    @variable(uc, 0 <= w <= scenario.wind)
    @constraint(uc, sum(g[i] for i in 1:N) + w == scenario.demand)
    ## !!! New: add binary on-off variables for each generator
    @variable(uc, u[i = 1:N], Bin)
    @constraint(uc, [i = 1:N], g[i] <= generators[i].max * u[i])
    @constraint(uc, [i = 1:N], g[i] >= generators[i].min * u[i])
    @objective(
        uc,
        Min,
        sum(generators[i].variable_cost * g[i] for i in 1:N) +
        wind.variable_cost * w +
        ## !!! new
        sum(generators[i].fixed_cost * u[i] for i in 1:N)
    )
    optimize!(uc)
    status = termination_status(uc)
    if status != MOI.OPTIMAL
        return (status = status,)
    end
    return (
        status = status,
        g = value.(g),
        w = value(w),
        wind_spill = scenario.wind - value(w),
        u = value.(u),
        total_cost = objective_value(uc),
    )
end

# Solve the economic dispatch problem
solution = solve_uc(generators, wind_generator, scenario)

println("Dispatch of Generators: ", solution.g, " MW")
println("Commitments of Generators: ", solution.u)
println("Dispatch of Wind: ", solution.w, " MW")
println("Wind spillage: ", solution.wind_spill, " MW")
println("\n")
println("Total cost: \$", solution.total_cost)

# ### Unit Commitment as a function of demand

# After implementing the UC model, we can now assess the interplay between the
# minimum power output constraints on generators and wind generation.

uc_df = DataFrames.DataFrame(
    demand = Float64[],
    commitment_G1 = Float64[],
    commitment_G2 = Float64[],
    dispatch_G1 = Float64[],
    dispatch_G2 = Float64[],
    dispatch_wind = Float64[],
    spillage_wind = Float64[],
    total_cost = Float64[],
)

for demand_scale in 0.2:0.1:1.5
    new_scenario = scale_demand(scenario, demand_scale)
    sol = solve_uc(generators, wind_generator, new_scenario)
    if sol.status == MOI.OPTIMAL
        push!(
            uc_df,
            (
                new_scenario.demand,
                sol.u[1],
                sol.u[2],
                sol.g[1],
                sol.g[2],
                sol.w,
                sol.wind_spill,
                sol.total_cost,
            ),
        )
    end
    println("Status: $(sol.status) for demand_scale = $(demand_scale)")
end

#-

uc_df

#-

commitment_plot = @df(
    uc_df,
    Plots.plot(
        :demand,
        [:commitment_G1, :commitment_G2],
        labels = ["G1" "G2"],
        title = "Committment",
        legend = :bottomright,
        linewidth = 3,
        xlabel = "Demand [MW]",
        ylabel = "Commitment decision {0, 1}"
    ),
)

dispatch_plot = @df(
    uc_df,
    Plots.plot(
        :demand,
        [:dispatch_G1, :dispatch_G2, :dispatch_wind],
        labels = ["G1" "G2" "Wind"],
        title = "Dispatch [MW]",
        legend = :bottomright,
        linewidth = 3,
        xlabel = "Demand",
        ylabel = "Dispatch [MW]"
    ),
)

Plots.plot(commitment_plot, dispatch_plot)

# ## Nonlinear economic dispatch

# As a final example, we modify our economic dispatch problem in two ways:
#
#  * The thermal cost function is user-defined
#  * The output of the wind is only the square-root of the dispatch

import Ipopt

"""
    thermal_cost_function(g)

A user-defined thermal cost function in pure-Julia! You can include
nonlinearities, and even things like control flow.

!!! warning
    It's still up to you to make sure that the function has a meaningful
    derivative.
"""
function thermal_cost_function(g)
    if g <= 500
        return g
    else
        return g + 1e-2 * (g - 500)^2
    end
end

function solve_nonlinear_ed(
    generators::Vector,
    wind,
    scenario;
    silent::Bool = false,
)
    model = Model(Ipopt.Optimizer)
    if silent
        set_silent(model)
    end
    register(model, :tcf, 1, thermal_cost_function; autodiff = true)
    N = length(generators)
    @variable(model, generators[i].min <= g[i = 1:N] <= generators[i].max)
    @variable(model, 0 <= w <= scenario.wind)
    @NLobjective(
        model,
        Min,
        sum(generators[i].variable_cost * tcf(g[i]) for i in 1:N) +
        wind.variable_cost * w,
    )
    @NLconstraint(model, sum(g[i] for i in 1:N) + sqrt(w) == scenario.demand)
    optimize!(model)
    return (
        g = value.(g),
        w = value(w),
        wind_spill = scenario.wind - value(w),
        total_cost = objective_value(model),
    )
end

solution = solve_nonlinear_ed(generators, wind_generator, scenario)

# Now let's see how the wind is dispatched as a function of the cost:

wind_cost = 0.0:1:100
wind_dispatch = Float64[]
for c in wind_cost
    sol = solve_nonlinear_ed(
        generators,
        WindGenerator(c),
        scenario;
        silent = true,
    )
    push!(wind_dispatch, sol.w)
end

Plots.plot(
    wind_cost,
    wind_dispatch,
    xlabel = "Cost",
    ylabel = "Dispatch [MW]",
    label = false,
 )

max_stars_repo_name: blegat/jump-training-materials

max_issues_repo_name: blegat/jump-training-materials

max_forks_repo_name: blegat/jump-training-materials


text:
using Printf

# @todo does Node need to be mutable?
mutable struct Node
    visit_count::Int32
    prior::Float32
    value_sum::Float32
    children::Dict{Int32,Node}
    hidden_state::Union{Nothing,AbstractArray{Float32, 3}}
    reward::Float32
    info::String
    state::AbstractEnvState # @todo - remove state altogether
    prior_policy::Dict
    action_scores::Dict
    function Node(prior, state; visit_count=0, value_sum=0,
        children=Dict(), hidden_state=nothing, reward=0)
        new(visit_count, prior, value_sum, children, hidden_state, reward, "", state, Dict(), Dict())
    end
end

gethiddenstate(n::Node) = n.hidden_state

gettoplay(n::Node) = n.state.to_play

function json(x::Node, visit_pr=0.0)
    d = Dict()
    noisy_prior_policy = Dict()
    for (action, node) in x.children
        noisy_prior_policy[action] = node.prior
    end
    d["visit_pr"] = visit_pr
    d["visit_count"] = x.visit_count
    d["to_play"] = get_player_label(x)
    d["prior"] = x.prior
    d["value_sum"] = x.value_sum
    d["prior_policy"] = x.prior_policy
    d["noisy_prior_policy"] = noisy_prior_policy
    d["reward"] = x.reward
    d["children"] = Dict()
    d["action_scores"] = x.action_scores
    total_child_visit_count = sum([child.visit_count for child in values(x.children)])
    for (action, child) in x.children
        visit_pr = child.visit_count / total_child_visit_count
        d["children"][action] = json(child, visit_pr)
    end
    d
end

get_player_label(x::Node) = x.state.to_play.idx == PLAYER_ID ? "A" : "B"

function Base.show(io::IO, x::Node)
    show_indented(io, x)
end

function show_indented(io::IO, x::Node, level=0, action=nothing)
    #show(io, x.channels)
    print(" "^(level*5))
    action_label = isnothing(action) ? "" : "A($action)"
    hex_state = "hex_str" #hex_str(x.state)
    @printf "%s %s pr: %.2f {c: %i, r: %i, vs: %.3f}|%s %s" action_label hex_state x.prior x.visit_count x.reward x.value_sum get_player_label(x) x.info
    println(io)
    actions = [action for action in keys(x.children)]
    for action in sort(actions)
        child = x.children[action]
        show_indented(io, child, level+1, action)
    end
end

expanded(x::Node) = length(x.children) > 0

function calc_value(x::Node)
    Float32(x.visit_count == 0 ? 0.0 : x.value_sum / x.visit_count)
end

function priors(x::Node)
    sort([(action, child.prior) for (action, child) in x.children])
end

#=
We expand a node using the value, reward and policy prediction obtained from
the neural network.
=#

function expand_node!(node::Node, actions::Vector, output::SingleNetworkOutput, reward::Number)
    p_logits = policy_logits(output)
    hidden = hidden_state(output)
    node.hidden_state = hidden
    node.reward = Float32(reward)

    for (action, pr) in enumerate(softmax(p_logits))
        node.prior_policy[action] = pr
    end

    if !isempty(actions)
        policy = Dict([a=>exp(p_logits[a.idx]) for a in actions])
        policy_sum = sum(values(policy))
        for (action, p) in policy
            node.children[action.idx] = Node(p / policy_sum, getnextstate(node.state, action)) #reward=get_reward(node.state, action)
        end
    end

end

function expand_node!(
    nodes::Vector{Node},
    actionSets::Vector,
    output::NetworkOutput,
    reward::Vector{Float32}
    )
    setfield!.(
        nodes,
        :hidden_state,
        eachslice(hidden_state(output), dims=4)
    )
    setfield!.(
        nodes,
        :reward,
        Float32.(reward)
    )
    #### @todo - efficiency could be much improved
    p_logits_tensor = policy_logits(output)
    for (idx, actions) in enumerate(actionSets)
        node = nodes[idx]
        p_logits = view(p_logits_tensor, :, idx)
        for (action, pr) in enumerate(softmax(p_logits))
            node.prior_policy[action] = pr
        end
        if !isempty(actions)
            policy = Dict([a=>exp(p_logits[a.idx]) for a in actions])
            policy_sum = sum(values(policy))
            for (action, p) in policy
                node.children[action.idx] = Node(p / policy_sum, getnextstate(node.state, action))
            end
        end
    end
    ####
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""SupplementalGroupsStrategyOptions defines the strategy type and options used to create the strategy. Deprecated: use SupplementalGroupsStrategyOptions from policy API Group instead.

    IoK8sApiExtensionsV1beta1SupplementalGroupsStrategyOptions(;
        ranges=nothing,
        rule=nothing,
    )

    - ranges::Vector{IoK8sApiExtensionsV1beta1IDRange} : ranges are the allowed ranges of supplemental groups.  If you would like to force a single supplemental group then supply a single range with the same start and end. Required for MustRunAs.
    - rule::String : rule is the strategy that will dictate what supplemental groups is used in the SecurityContext.
"""
mutable struct IoK8sApiExtensionsV1beta1SupplementalGroupsStrategyOptions <: SwaggerModel
    ranges::Any # spec type: Union{ Nothing, Vector{IoK8sApiExtensionsV1beta1IDRange} } # spec name: ranges
    rule::Any # spec type: Union{ Nothing, String } # spec name: rule

    function IoK8sApiExtensionsV1beta1SupplementalGroupsStrategyOptions(;ranges=nothing, rule=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1SupplementalGroupsStrategyOptions, Symbol("ranges"), ranges)
        setfield!(o, Symbol("ranges"), ranges)
        validate_property(IoK8sApiExtensionsV1beta1SupplementalGroupsStrategyOptions, Symbol("rule"), rule)
        setfield!(o, Symbol("rule"), rule)
        o
    end
end # type IoK8sApiExtensionsV1beta1SupplementalGroupsStrategyOptions

const _property_map_IoK8sApiExtensionsV1beta1SupplementalGroupsStrategyOptions = Dict{Symbol,Symbol}(Symbol("ranges")=>Symbol("ranges"), Symbol("rule")=>Symbol("rule"))
const _property_types_IoK8sApiExtensionsV1beta1SupplementalGroupsStrategyOptions = Dict{Symbol,String}(Symbol("ranges")=>"Vector{IoK8sApiExtensionsV1beta1IDRange}", Symbol("rule")=>"String")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1SupplementalGroupsStrategyOptions }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1SupplementalGroupsStrategyOptions))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1SupplementalGroupsStrategyOptions }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiExtensionsV1beta1SupplementalGroupsStrategyOptions[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1SupplementalGroupsStrategyOptions }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1SupplementalGroupsStrategyOptions[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1SupplementalGroupsStrategyOptions)
    true
end

function validate_property(::Type{ IoK8sApiExtensionsV1beta1SupplementalGroupsStrategyOptions }, name::Symbol, val)
end


text:
"""
    Cylinder{N, T}

A `Cylinder` is a 2D rectangle or a 3D cylinder defined by its origin point,
its extremity and a radius. `origin`, `extremity` and `r`, must be specified.
"""
struct Cylinder{N,T} <: GeometryPrimitive{N,T}
    origin::Point{N,T}
    extremity::Point{N,T}
    r::T
end

"""
    Cylinder2{T}
    Cylinder3{T}

A `Cylinder2` or `Cylinder3` is a 2D/3D cylinder defined by its origin point,
its extremity and a radius. `origin`, `extremity` and `r`, must be specified.
"""
const Cylinder2{T} = Cylinder{2,T}
const Cylinder3{T} = Cylinder{3,T}

origin(c::Cylinder{N,T}) where {N,T} = c.origin
extremity(c::Cylinder{N,T}) where {N,T} = c.extremity
radius(c::Cylinder{N,T}) where {N,T} = c.r
height(c::Cylinder{N,T}) where {N,T} = norm(c.extremity - c.origin)
direction(c::Cylinder{N,T}) where {N,T} = (c.extremity .- c.origin) ./ height(c)

function rotation(c::Cylinder{2,T}) where {T}
    d2 = direction(c)
    u = @SVector [d2[1], d2[2], T(0)]
    v = @MVector [u[2], -u[1], T(0)]
    normalize!(v)
    return hcat(v, u, @SVector T[0, 0, 1])
end

function rotation(c::Cylinder{3,T}) where {T}
    d3 = direction(c)
    u = @SVector [d3[1], d3[2], d3[3]]
    if abs(u[1]) > 0 || abs(u[2]) > 0
        v = @MVector [u[2], -u[1], T(0)]
    else
        v = @MVector [T(0), -u[3], u[2]]
    end
    normalize!(v)
    w = @SVector [u[2] * v[3] - u[3] * v[2], -u[1] * v[3] + u[3] * v[1],
                  u[1] * v[2] - u[2] * v[1]]
    return hcat(v, w, u)
end

function coordinates(c::Cylinder{2,T}, nvertices=(2, 2)) where {T}
    r = Rect(c.origin[1] - c.r / 2, c.origin[2], c.r, height(c))
    M = rotation(c)
    points = coordinates(r, nvertices)
    vo = to_pointn(Point3{T}, origin(c))
    return (M * (to_pointn(Point3{T}, point) .- vo) .+ vo for point in points)
end

function faces(sphere::Cylinder{2}, nvertices=(2, 2))
    return faces(Rect(0, 0, 1, 1), nvertices)
end

function coordinates(c::Cylinder{3,T}, nvertices=30) where {T}
    if isodd(nvertices)
        nvertices = 2 * (nvertices ÷ 2)
    end
    nvertices = max(8, nvertices)
    nbv = nvertices ÷ 2

    M = rotation(c)
    h = height(c)
    range = 1:(2 * nbv + 2)
    function inner(i)
        return if i == length(range)
            c.extremity
        elseif i == length(range) - 1
            origin(c)
        else
            phi = T((2π * (((i + 1) ÷ 2) - 1)) / nbv)
            up = ifelse(isodd(i), 0, h)
            (M * Point(c.r * cos(phi), c.r * sin(phi), up)) .+ c.origin
        end
    end

    return (inner(i) for i in range)
end

function faces(c::Cylinder{3}, facets=30)
    isodd(facets) ? facets = 2 * div(facets, 2) : nothing
    facets < 8 ? facets = 8 : nothing
    nbv = Int(facets / 2)
    indexes = Vector{TriangleFace{Int}}(undef, facets)
    index = 1
    for j in 1:(nbv - 1)
        indexes[index] = (index + 2, index + 1, index)
        indexes[index + 1] = (index + 3, index + 1, index + 2)
        index += 2
    end
    indexes[index] = (1, index + 1, index)
    indexes[index + 1] = (2, index + 1, 1)

    for i in 1:length(indexes)
        i % 2 == 1 ? push!(indexes, (indexes[i][1], indexes[i][3], 2 * nbv + 1)) :
        push!(indexes, (indexes[i][2], indexes[i][1], 2 * nbv + 2))
    end
    return indexes
end


text:
"""
    get_prior_neighbors(g, ordering, i, exclude=[])
Return the neighbors of `ordering[i]` in `g` with index inferior to i.
"""
function get_prior_neighbors(g, ordering, i, exclude=[])
    prior_neighbors = []
    for j in 1:i-1
        if  !(j in exclude) && (ordering[j] in neighbors(g, ordering[i]))
            push!(prior_neighbors, ordering[j])
        end
    end
    return prior_neighbors
end

"""
    get_later_neighbors(g, ordering, i, exclude=[])
Return the neighbors of `ordering[i]` in `g` with index superior to i.
"""
function get_later_neighbors(g, ordering, i, exclude=[])
    later_neighbors = []
    for j in i+1:length(ordering)
        if !(j in exclude) && (ordering[j] in neighbors(g, ordering[i]))
            push!(later_neighbors, ordering[j])
        end
    end
    return later_neighbors
end

"""
    get_closest_neighbor(g::AbstractGraph, ordering, i)
Return the index of the closest neighbor of the value in `ordering` corresponding to the index `i`.
The return index is inferior to `j`
"""
function get_closest_neighbor(g, ordering, i)
    for j in i-1:-1:1
        if ordering[j] in neighbors(g, ordering[i])
            return j
        end
    end
    return nothing
end

"""
    refine(S, X)
Refine the set `S` into tow sets using `X`.
"""
function refine(S, X)
    intersection = intersect(S, X)
    difference = setdiff(S, X)
    if isempty(intersection)
        return [difference]
    elseif isempty(difference)
        return [intersection]
    else
        return [intersection, difference]
    end
end

"""
    refineset(S, X)
Refine all the sets inside `S` with the set `X`.
"""
function refineset(S, X)
    new_S = []
    for set in S
        append!(new_S, refine(set, X))
    end
    return new_S
end

"""
    lbfs(g)
Return a lexicogic ordering of `G`.
"""
function lbfs(g)
    S = [[vertices(g)...]]
    visited = []
    ordering = []
    for v in vertices(g)
        current_vertex = S[1][1]
        push!(ordering, current_vertex)
        setdiff!(S[1], current_vertex)
        if isempty(S[1])
            popfirst!(S)
        end
        nghbrs = setdiff(neighbors(g, current_vertex), ordering)
        S = refineset(S, nghbrs)
    end
    return ordering
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
using .IntervalArithmetic

function evaluate(a::STaylor1{N,T}, dx::Interval) where {N, T <: Number}
    order = N - 1
    uno = one(dx)
    dx2 = dx^2
    if iseven(order)
        kend = order-2
        @inbounds sum_even = a[end]*uno
        @inbounds sum_odd = a[end-1]*zero(dx)
    else
        kend = order-3
        @inbounds sum_odd = a[end]*uno
        @inbounds sum_even = a[end-1]*uno
    end
    @inbounds for k = kend:-2:0
        sum_odd = sum_odd*dx2 + a[k + 1]
        sum_even = sum_even*dx2 + a[k]
    end
    return sum_even + sum_odd*dx
end

normalize_taylor(a::STaylor1{N,T}, I::Interval{T}, symI::Bool=true) where {N, T <: Number} =
    _normalize(a, I, Val(symI))

function _normalize(a::STaylor1{N,T}, I::Interval{T}, ::Val{true}) where {N, T <: Number}
    t = STaylor1(zero(T), Val{N}())
    tnew = mid(I) + t*radius(I)
    return a(tnew)
end

function _normalize(a::STaylor1{N,T}, I::Interval{T}, ::Val{false}) where {N, T <: Number}
    t = STaylor1(zero(promote_type(S,T)), Val{N}())
    tnew = inf(I) + t*diam(I)
    return a(tnew)
end


text:
using DataInterpolations, Test

@testset "Interface" begin include("interface.jl") end
@testset "Interpolation Tests" begin include("interpolation_tests.jl") end


text:
using ApproxManifoldProducts
# using Distributions

##

# these do work
ptsA = [-pi/2; pi/2]

# these dont work with difftheta, but fine with `-`
ptsB = [-pi/2; pi/2+1e-10]

pA = kde!(ptsA, [0.1], (addtheta,), (difftheta,))
pB = kde!(ptsB, [0.1], (addtheta,), (difftheta,))

##

ept = ones(1,1)
ept[1,1] = pi/2

resA = evaluateDualTree(pA, ept, false, 1e-3, (addtheta,), (difftheta,))[1]

println("")
println("")

resB = evaluateDualTree(pB, ept, false, 1e-3, (addtheta,), (difftheta,))[1]

##

pA.bt.centers

pB.bt.centers

#


text:
module PhasePlots

using RecipesBase, ColorTypes

@userplot PhasePlot

@recipe function f(p::PhasePlot; clim=nothing, γ=2)
    x, y, z = p.args
    isnothing(clim) && (clim = maximum(abs, z))
    extend(x, n) = (h = (x[end]-x[begin])/2(n-1); [x[begin]-h, x[end]+h])
    x = extend(x, size(z,1))
    xlims --> x
    y = extend(y, size(z,2))
    ylims --> y
    # The docs say :yflip defaults to false, Plots.default agrees, but the plot with yflip omitted is the same as when it is set to true.
    xflip --> !issorted(x)
    yflip --> !issorted(y)
    aspect_ratio --> 1
    grid --> false
    tick_direction --> :out

    N = transpose(z)
    x, y, @. HSVA(360/2π*angle(N), 1, !iszero(N), (abs(N)/clim)^γ)
end

@doc """
    phaseplot(x, y, z; clim=maximum(abs.(z)))

Plot a complex function with phase as hue

The fourth roots of unity are red, green, cyan and purple.  Modulus is mapped to the alpha channel, so the classic complex phase plot is `phaseplot(..., bg=:black)`.  Alpha becomes opaque at the modulus clim; this is useful when multiple plots should have the same brightness scale.

Plots.jl is dreadfully inconsistent at mapping points in the plane to raster indices.  For consistency, phaseplot has a different default to every existing function.  The first dimension is the real/horizontal axis, the second the imaginary/vertical one, and by default x increases rightward and y upward.

When Plots.jl draws a raster, the extrema of x and y are the coordinates of points half a pixel outside the edge of the grid.  This function moves them to the edges, which is almost always what you will want.
""" phaseplot

end # module


text:
# EnsembleMethods
using Distributions
using DataFrames
using Base.Threads

import Base.convert

const NO_BEST = (0, 0)

immutable Leaf
    majority::Any
    values::Vector
end

immutable Node
    col_idx::Int
    split_value::Any
    surrogates::Array{Tuple{Int, Real}, 1}                # This holds our surrogate vars

    # pointers to daughter nodes
    left::Union{Leaf, Node}
    right::Union{Leaf, Node}
end

typealias LeafOrNode Union{Leaf, Node}

immutable Ensemble
    trees::Vector{Node}
end

function convert(::Type{Node}, x::Leaf)
    warn("converting leaf to node")
    return Node(0, nothing, Array{Tuple{Int, Real}, 1}(5), x, Leaf(nothing,[nothing]))
end

promote_rule(::Type{Node}, ::Type{Leaf}) = Node
promote_rule(::Type{Leaf}, ::Type{Node}) = Node

include("utils.jl")
include("measures.jl")
include("classification.jl")
include("surrogates.jl")
include("apply_models.jl")
include("random_forest_regression.jl")


text:
using DiffRules, SpecialFunctions, NaNMath

for (M, f, arity) in DiffRules.diffrules()
  arity == 1 || continue
  @eval begin
    @adjoint $M.$f(x::Real) = $M.$f(x),
      Δ -> (Δ * $(DiffRules.diffrule(M, f, :x)),)
  end
end

for (M, f, arity) in DiffRules.diffrules()
  arity == 2 || continue
  da, db = DiffRules.diffrule(M, f, :a, :b)
  @eval begin
    @adjoint $M.$f(a::Real, b::Real) = $M.$f(a, b),
      Δ -> (Δ * $da, Δ * $db)
  end
end

@adjoint Base.convert(T::Type{<:Real}, x::Real) = convert(T, x), Δ -> (nothing, Δ)

@adjoint Base.:+(xs...) = +(xs...), Δ -> map(_ -> Δ, xs)

@adjoint function sincos(x)
  s, c = sincos(x)
  (s, c), ((s̄, c̄),) -> (s̄*c - c̄*s,)
end


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

using DynamicHMCModels

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

ProjDir = @__DIR__
cd(ProjDir)

# Read the milk data

df = CSV.read(joinpath("..", "..", "data", "milk.csv"), DataFrame)
df = filter(row -> !(row[:neocortex_perc] == "NA"), df)
#df[:, :kcal_per_g] = convert(Vector{Float64}, df[:, :kcal_per_g])
df[!, :log_mass] = log.(convert(Vector{Float64}, df[:, :mass]))

# Define the model struct

Base.@kwdef mutable struct MilkModel{Ty <: AbstractVector, Tx <: AbstractMatrix}
    "Observations."
    y::Ty
    "Covariates"
    x::Tx
end

# Write a function to return properly dimensioned transformation.

function make_transformation(model::MilkModel)
  as((β = as(Array, size(model.x, 2)), σ = asℝ₊))
end

# Instantiate the model with data and inits.

x = hcat(ones(size(df, 1)), df[:, :log_mass]);
model = MilkModel(;y=df[:, :kcal_per_g], x=x)

# Make the type callable with the parameters *as a single argument*.

function (model::MilkModel)(θ)
    @unpack y, x, = model   # extract the data
    @unpack β, σ = θ            # works on the named tuple too
    ll = 0.0
    ll += logpdf(Normal(0, 100), x[1])
    ll += logpdf(Normal(0, 1), x[2])
    ll += logpdf(TDist(1.0), σ)
    ll += loglikelihood(Normal(0, σ), y .- x*β)
    ll
end

println()
model((β = [1.0, 2.0], σ = 1.0))
println()

# Wrap the problem with a transformation, then use Flux for the gradient.

P = TransformedLogDensity(make_transformation(model), model)
∇P = ADgradient(:ForwardDiff, P);

# Tune and sample.

results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000)
posterior = P.transformation.(results.chain)

p = as_particles(posterior)
display(p)

stan_result = "
Iterations = 1:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 1000

Empirical Posterior Estimates:
          Mean         SD        Naive SE       MCSE      ESS
    a  0.70472876 0.057040655 0.00090189195 0.0011398893 1000
   bm -0.03150330 0.023642759 0.00037382484 0.0004712342 1000
sigma  0.18378372 0.039212805 0.00062000888 0.0011395979 1000

Quantiles:
          2.5%       25.0%       50.0%        75.0%       97.5%
    a  0.59112968  0.66848775  0.70444950  0.741410500 0.81915225
   bm -0.07729257 -0.04708425 -0.03104865 -0.015942925 0.01424901
sigma  0.12638780  0.15605950  0.17800600  0.204319250 0.27993590
";

# End of `05/5.6d.jl`


text:
using LightXML, Compat

xdoc = XMLDocument()

xroot = create_root(xdoc, "States")

xs1 = new_child(xroot, "State")
add_text(xs1, "Massachusetts")
set_attribute(xs1, "tag", "MA")

xs2 = new_child(xroot, "State")
add_text(xs2, "Illinois")
set_attributes(xs2, @compat Dict{Any,Any}("tag"=>"IL", "cap"=>"Springfield"))

xs3 = new_child(xroot, "State")
add_text(xs3, "California")
set_attributes(xs3; tag="CA", cap="Sacramento")

rtxt = """
<?xml version="1.0" encoding="utf-8"?>
<States>
  <State tag="MA">Massachusetts</State>
  <State tag="IL" cap="Springfield">Illinois</State>
  <State tag="CA" cap="Sacramento">California</State>
</States>
"""

@assert strip(string(xdoc)) == strip(rtxt)

free(xdoc)