text:
module TestSKL

using TSML
using TSMLextra
using PyCall
using Test

const IRIS = getiris()
const X = IRIS[:,1:4] |> Matrix
const Y = IRIS[:,5] |> Vector

const XX = IRIS[:,1:1] |> Matrix
const YY = IRIS[:,4] |> Vector

const classifiers = [
    "LinearSVC","QDA","MLPClassifier","BernoulliNB",
    "RandomForestClassifier","LDA",
    "NearestCentroid","SVC","LinearSVC","NuSVC","MLPClassifier",
    "RidgeClassifierCV","SGDClassifier","KNeighborsClassifier",
    "GaussianProcessClassifier","DecisionTreeClassifier",
    "PassiveAggressiveClassifier","RidgeClassifier",
    "ExtraTreesClassifier","GradientBoostingClassifier",
    "BaggingClassifier","AdaBoostClassifier","GaussianNB","MultinomialNB",
    "ComplementNB","BernoulliNB"
 ]

const regressors = [
    "SVR",
    "Ridge",
    "RidgeCV",
    "Lasso",
    "ElasticNet",
    "Lars",
    "LassoLars",
    "OrthogonalMatchingPursuit",
    "BayesianRidge",
    "ARDRegression",
    "SGDRegressor",
    "PassiveAggressiveRegressor",
    "KernelRidge",
    "KNeighborsRegressor",
    "RadiusNeighborsRegressor",
    "GaussianProcessRegressor",
    "DecisionTreeRegressor",
    "RandomForestRegressor",
    "ExtraTreesRegressor",
    "GradientBoostingRegressor",
    "MLPRegressor",
    "AdaBoostRegressor"
#    "IsotonicRegression"
]

function fit_test(learner::String,in::T,out::Vector) where {T<:Union{Matrix,Vector}}
	_learner=SKLearner(Dict(:learner=>learner))
	fit!(_learner,in,out)
	@test _learner.model != nothing
	return _learner
end

function fit_transform_reg(model::TSLearner,in::T,out::Vector) where {T<:Union{Matrix,Vector}}
    @test sum((transform!(model,in) .- out).^2)/length(out) < 2.0
end

@testset "scikit classifiers" begin
    for cl in classifiers
	fit_test(cl,X,Y)
    end
end

@testset "scikit regressors" begin
    for rg in regressors
	model=fit_test(rg,XX,YY)
	fit_transform_reg(model,XX,YY)
    end
end

end


text:
# Script to run through the files and fit the zero-inflated negative binomial distrubution

using Plots, Distributions, DelimitedFiles, Base.Printf
include("Utils.jl")
import Main.Utils
gr()

# Set some parameters
spread = 0.95
folder = "Chains/"
Files = readdir(folder)

# Make a stats variable to be writtten to a file
Stats = Array{Any,2}(undef,length(Files)+1,13)
Stats[1,:] = ["Name" "p" "w" "r" "K" "r_min" "r_max" "p_min" "p_max" "w_min" "w_max" "K_min" "K_max"]

for (ii,file) in enumerate(Files)

    # Load and plot the distribution data
    chain = readdlm(folder*file)

	# Extract parameters and plot
    r = Utils.find_MAP(chain,idx=1)
    p = Utils.find_MAP(chain,idx=2)
    w = Utils.find_MAP(chain,idx=3)
    K = Utils.find_MAP(chain,idx=4)
    int_r = Utils.credibleintervals(chain,idx=1, spread=spread)
    int_p = Utils.credibleintervals(chain,idx=2, spread=spread)
    int_w = Utils.credibleintervals(chain,idx=3, spread=spread)
    int_K = Utils.credibleintervals(chain,idx=4, spread=spread)

    # Plot and save
    plt = Utils.plot_chain(chain, [r,p,w,K], [int_r,int_p,int_w,int_K])
    try
        Plots.pdf("MCMC/"*file)
    catch
        mkdir("MCMC")
        Plots.pdf("MCMC/"*file)
    end

    # Save data
    println(file)
    Stats[ii+1,:] = [file p w r K r-int_r[1] int_r[2]-r p-int_p[1] int_p[2]-p w-int_w[1] int_w[2]-w K-int_K[1] int_K[2]-K]

end

DelimitedFiles.writedlm("Statistics.txt", Stats)

ext: jl

max_stars_repo_name: rdbrackston/tx-analysis

max_issues_repo_name: rdbrackston/tx-analysis

max_forks_repo_name: rdbrackston/tx-analysis


text:
using LinearAlgebra
using Plots

x = randn(100)
A = I(100) .+ 1/norm(x,2) .* randn(100, 100)
y = A*x + 1/norm(A, 2) .* randn(100)

# Solve for x
x̂ = A \ y
# Compute the error
ϵ = x - x̂

# Make a plot and save
plot_1 = scatter(x, y,
    ylabel = "y", label = nothing
    )

plot_2 = scatter(x, ϵ,
    xlabel = "x", ylabel = "Error",
    label = nothing)

plot(plot_1, plot_2, layout = (2,1))
savefig(joinpath(pwd(), "figures", "linear_solve.pdf"))


text:
using Polyhedra, CDDLib
P = polyhedron(vrep(eye(2)), CDDLibrary())
P * P


text:
using Test
using Mill, BenchmarkTools, SparseArrays, Random, Flux
using Mill: NGramIterator, ngrams, string2ngrams, countngrams, mul, multrans, NGramMatrix, catobs

@testset "ngrams" begin
    x = [1,3,5,2,6,8,3]
    b = 8 + 1

    slicer = (x,n) -> map(i -> x[(max(i-n+1,1):min(i,length(x)))],1:length(x)+n-1)
    indexes = (x,b) -> mapreduce(i -> i[2]*b^(i[1]-1), +, enumerate(reverse(x)))
    function idx2vec(i,n)
        o = zeros(Int,n)
        for v in i
            o[mod(v,n)+1] +=1
        end
        o
    end

    @testset "testing ngrams on vector of Ints" begin
        @test all(ngrams(x,3,b) .== map(x -> indexes(x,b),slicer(x,3)))
        @test all(ngrams(x,2,b) .== map(x -> indexes(x,b),slicer(x,2)))
        @test all(ngrams(x,1,b) .== map(x -> indexes(x,b),slicer(x,1)))
    end

    @testset "testing frequency of ngrams on vector of Ints and on Strings" begin
        @test all(countngrams(x,3,b,10) .== idx2vec(map(x -> indexes(x,b), slicer(x,3)), 10))
        for s in split("Lorem ipsum dolor sit amet, consectetur adipiscing elit")
            @test all(countngrams(s,3,256,10) .== idx2vec(ngrams(s,3,256), 10))
        end

        s = split("Lorem ipsum dolor sit amet, consectetur adipiscing elit")
        @test all(countngrams(s,3,256,10) .== hcat(map(ss -> idx2vec(ngrams(ss,3,256), 10),s)...))
    end

    @testset "string2ngrams" begin
        @test size(string2ngrams(["","a"],3,2053)) == (2053,2)
    end

    @testset "NGramMatrix" begin
        @test all(collect(NGramIterator(codeunits("hello"), 3, 257)) .== ngrams("hello", 3, 257))

        A = randn(4, 10)
        n = size(A, 2)
        s = ["hello", "world", "!!!"]
        si = map(i -> Int.(codeunits(i)), s)
        B = NGramMatrix(s, 3, 256, n)
        Bi = NGramMatrix(si, 3, 256, n)
        @test all(A * B ≈ A*string2ngrams(s, 3, n))
        @test all(A * Bi ≈ A * B)
        A = randn(5,3)
        @test all(multrans(A , B) ≈ A*transpose(string2ngrams(s, 3, n)))
        @test all(multrans(A , B) ≈ multrans(A , B))
    end

    @testset "NGramMatrix to SparseMatrix" begin
        A = randn(4, 10)
        n = size(A, 2)
        s = ["hello", "world", "!!!"]
        si = map(i -> Int.(codeunits(i)), s)
        B = NGramMatrix(s, 3, 256, n)
        Bi = NGramMatrix(si, 3, 256, n)
        @test all(A * B ≈  A * SparseMatrixCSC(B))
        @test all(A * B ≈  A * SparseMatrixCSC(Bi))
    end

    @testset "integration with MILL & Flux" begin
        s = ["hello", "world", "!!!"]
        si = map(i -> Int.(codeunits(i)), s)
        for (a, s) in [(NGramMatrix(s, 3, 256, 2057), s), (NGramMatrix(si, 3, 256, 2057), si)]
            @test all(reduce(catobs, [a, a]).s .== vcat(s,s))
            @test all(hcat(a,a).s .== vcat(s,s))

            W = randn(40, 2057)
            @test gradcheck(x -> sum(x * a), W)

            a = ArrayNode(a, nothing)
            @test all(reduce(catobs, [a, a]).data.s .== vcat(s,s))
            a = BagNode(a, [1:3], nothing)
            @test all(reduce(catobs, [a, a]).data.data.s .== vcat(s,s))
        end
    end
end

begin
    println("Benchmarking multiplication")
    #begin block body
    A = randn(80,2053);
    s = [randstring(10) for i in 1:1000];
    B = NGramMatrix(s, 3, 256, 2053)
    C = sparse(string2ngrams(s, 3, size(A, 2)));
    println("A * B::NGramMatrix (This should be the fastest)");
    @btime A*B;																	# 526.456 μs (2002 allocations: 671.95 KiB)
    println("A * string2ngrams(s, 3, size(A, 2))")
    @btime A*string2ngrams(s, 3, size(A, 2)); 					# 154.646 ms (3013 allocations: 16.38 MiB)
    println("A * sparse(string2ngrams(s, 3, size(A, 2)))")
    @btime A*sparse(string2ngrams(s, 3, size(A, 2))); 	# 7.525 ms (3013 allocations: 16.57 MiB)
    print("A * C where C = sparse(string2ngrams(s, 3, size(A, 2)));");
    @btime A*C; 																				# 1.527 ms (2 allocations: 625.08 KiB)
end


text:
#=
main:
- Julia version: 1.3.1
- Author: abd_alazez_ahmed
- Date: 2020-04-02
=#
using Plots
pyplot()
thetaRange=range(-pi,stop=pi,length=100)
println("ϵₐ = ")
ep_a = 1200
println("ϵᵦ = ")
ep_b = 600
println("ϵᵪ = ")
ep_c = -700
println("θₐ = ")
theta_a = -60
println("θᵦ = ")
theta_ab = 60
println("θᵪ = ")
theta_ac = 120
ep = inv([sind(theta_ab)^2 0.5*sind(2*theta_ab); sind(theta_ac)^2 0.5*sind(2*theta_ac)])*[ep_b - ep_a*cosd(theta_ab)^2 ; ep_c - ep_a*cosd(theta_ab)^2]
final_ep = [cosd(theta_a)^2 sind(theta_a)^2 sind(2*theta_a);sind(theta_a)^2 cosd(theta_a)^2 -sind(2*theta_a);-0.5*sind(2*theta_a) 0.5*sind(2*theta_a) cosd(2*theta_a)]*[ep_a;ep[1];ep[2]/2]
println(ep)
final_ep[3]*=2
println(final_ep)
E=2e11
v=0.3
sigMat=E/(1-v^2)*[1 v 0;v 1 0;0 0 (1-v)/2]*final_ep
sigma_x = sigMat[1]
sigma_y = sigMat[2]
tau = sigMat[3]
r = sqrt((sigma_x-sigma_y)^2 /4 + tau^2)
plot((sigma_x+sigma_y)/2 .+ r*cos.(thetaRange), r.*sin.(thetaRange),aspect_ratio=1,framestyle = :origin)
plot!([sigma_x,sigma_x],[0,-tau])
hline!([-tau])
plot!([sigma_y,(sigma_x+sigma_y)/2 + 1.1*r],[-tau,0.1*r*(tau)/((sigma_x+sigma_y)/2 + r-sigma_y)])
plot!([sigma_y,(sigma_x+sigma_y)/2],[-tau,-r])
savefig("Stress.svg")
r = sqrt((final_ep[1]-final_ep[2])^2 /4 + final_ep[3]^2)
plot((final_ep[1]+final_ep[2])/2 .+ r*cos.(thetaRange), r.*sin.(thetaRange),aspect_ratio=1,framestyle = :origin)
plot!([final_ep[1],final_ep[1]],[0,-final_ep[3]])
hline!([-final_ep[3]])
plot!([final_ep[2],(final_ep[1]+final_ep[2])/2 + 1.1*r],[-final_ep[3],0.1*r*(final_ep[3])/((final_ep[1]+final_ep[2])/2 + r-final_ep[2])])
plot!([final_ep[2],(final_ep[1]+final_ep[2])/2],[-final_ep[3],-r])
savefig("Strain.svg")


text:
## File to hold all structures

mutable struct ModPar
  β :: Float64
  γ :: Float64
  η :: Float64
  ξ :: Float64
  χ :: Float64

  function ModPar(β=0.92, γ=1.5, η=0.3, ξ = 0.10, χ = .0 )
    new(β,γ,η,ξ,χ)
  end

end


text:
using Pkg
pkg"add MultivariateStats"
pkg"add NearestNeighbors"
pkg"add StatsBase"

include("simpleicp.jl")

dataset = "all"
write_results = false

if cmp(dataset, "Dragon") == 0 || cmp(dataset, "all") == 0
    println("Processing dataset \"Dragon\"")
    X_fix = readdlm("../data/dragon1.xyz")
    X_mov = readdlm("../data/dragon2.xyz")
    H, X_mov_transformed = simpleicp(X_fix, X_mov)
end

if cmp(dataset, "Airborne Lidar") == 0 || cmp(dataset, "all") == 0
    println("Processing dataset \"Airborne Lidar\"")
    X_fix = readdlm("../data/airborne_lidar1.xyz")
    X_mov = readdlm("../data/airborne_lidar2.xyz")
    H, X_mov_transformed = simpleicp(X_fix, X_mov)
end

if cmp(dataset, "Terrestrial Lidar") == 0 || cmp(dataset, "all") == 0
    println("Processing dataset \"Terrestrial Lidar\"")
    X_fix = readdlm("../data/terrestrial_lidar1.xyz")
    X_mov = readdlm("../data/terrestrial_lidar2.xyz")
    H, X_mov_transformed = simpleicp(X_fix, X_mov)
end

if cmp(dataset, "Bunny") == 0 || cmp(dataset, "all") == 0
    println("Processing dataset \"Bunny\"")
    X_fix = readdlm("../data/bunny_part1.xyz")
    X_mov = readdlm("../data/bunny_part2.xyz")
    H, X_mov_transformed = simpleicp(X_fix, X_mov, max_overlap_distance=1)
end

# Export original and adjusted point clouds to xyz files to check the result
if write_results
    target_dir = "check"
    mkpath(target_dir)
    writedlm(joinpath(target_dir, "X_fix.xyz"), X_fix)
    writedlm(joinpath(target_dir, "X_mov.xyz"), X_mov)
    writedlm(joinpath(target_dir, "X_mov_transformed.xyz"), X_mov_transformed)
end

max_stars_repo_name: mitjap/simpleICP

max_issues_repo_name: mitjap/simpleICP

max_forks_repo_name: mitjap/simpleICP


text:
### Hastings & Powell 1991
## Reproducing teacup figure (figure 2D)

#### FIGURES 2D
#### b1 = 3.0

u0 = [0.76, 0.16, 9.9]  #initial conditions for x,y,z
tspan = (0.0, 10000.0) #timespan
p_b3 = [5.0, 0.1, 3.0, 2.0, 0.4, 0.01] # parameter values (a1, a2, b1, b2, c1, c2)

# Solving the system of differential equations
prob_b3 = ODEProblem(parameterized_model!, u0, tspan, p_b3)
sol_b3 = solve(prob_b3, alg = RK4(), reltol = 1e-14) # RK4 algorithm

sol_x_b3 = sol_b3[1, :] # timeseries for x
sol_y_b3 = sol_b3[2, :] # timeseries for y
sol_z_b3 = sol_b3[3, :] # timeseries for z

# Teacup Figures 2D and 5E (Three-dimensional phase plot)
fig2d = plot(sol_x_b3, sol_y_b3, sol_z_b3, # sol_b3, vars = (1,2,3),
             linestyle = :dot,
             xlabel = "x", xflip = false, xlim = (0.0, 1.0), xticks = 0.0:0.5:1.0,
             ylabel = "y", yflip = true, ylim = (0.0, 0.5), yticks = 0.0:0.25:0.5,
             zlabel = "z", zflip = false, zlim = (7.5, 10.5), zticks = 7.5:1.5:10.5,
             tick_direction = :out, guidefontsize = 10, tickfontsize = 8,
             legend = :none, size = (600, 600), c = :black,
             right_margin = 5Plots.mm, left_margin = 2Plots.mm,
             )
# Seems like there's a bug with the axes labels when flipping y
# Only solution we could find was to substract y values from 0.5, then plot and replace labels
sol_y_b3_flip = copy(sol_y_b3)
sol_y_b3_flip = 0.5 .- sol_y_b3_flip
fig2d_flip = plot(sol_x_b3, sol_y_b3_flip, sol_z_b3, #sol_b3,vars=(1,2,3),
             linestyle = :dot,
             xlabel = "x", xflip = false, xlim = (0.0, 1.0), xticks = 0.0:0.5:1.0, xrotation = 30.0,
             ylabel = "y", yflip = false, ylim = (0.0, 0.5), yticks = (0.0:0.25:0.5, string.(0.5:-0.25:0.0)), yrotation=-10.0,
             zlabel = "z", zflip = false, zlim = (7.5, 10.5), zticks = 7.5:1.5:10.5, zrotation=0.0,
             tick_direction=:out, guidefontsize = 10, tickfontsize = 8,
             legend =:none, size = (600,600), c=:black,
             right_margin = 5Plots.mm, left_margin = 2Plots.mm,
             )
# Otherwise, it's also possible to use the pyplot backend (additionnal dependency)
#=
pyplot()
fig2d_pyplot = plot(sol_x_b3, sol_y_b3, sol_z_b3, #sol_b3,vars=(1,2,3),
             linewidth = 0.1,
             xlabel = "x", xflip = false, xlim = (0.0, 1.0), xticks = 0.0:0.5:1.0,
             ylabel = "y", yflip = true, ylim = (0.0, 0.5), yticks = 0.0:0.25:0.5,
             zlabel = "z", zflip = false, zlim = (7.5,10.5), zticks = 7.5:1.5:10.5,
             tick_direction = :out, guidefontsize = 12, tickfontsize = 12,
             legend = :none, size = (600,600), c = :black)
gr()
=#

# Export figures
# savefig(fig2d, joinpath("..", "article", "figures", "fig2d_misplaced"))
savefig(fig2d_flip, joinpath("..", "article", "figures", "fig2d.pdf"))
# savefig(fig2d_pyplot, joinpath("..", "article", "figures", "fig2d_pyplot"))


text:
module SettingsWindow

using Gtk, Gtk.ShortNames, Gtk.GConstants
using ..Ahorn, Maple

settingsWindow = nothing

function getOptions(data::Dict{Any, Any}, dropdownOptions::Dict{String, Any}, langdata::Ahorn.LangData)
    res = Ahorn.Form.Option[]

    names = get(langdata, :names)
    tooltips = get(langdata, :tooltips)

    for (attr, value) in data
        sattr = Symbol(attr)
        keyOptions = get(dropdownOptions, attr, nothing)
        displayName = haskey(names, sattr) ? names[sattr] : Ahorn.humanizeVariableName(attr)
        tooltip = Ahorn.expandTooltipText(get(tooltips, sattr, ""))

        push!(res, Ahorn.Form.suggestOption(displayName, value, tooltip=tooltip, dataName=attr, choices=keyOptions, editable=true))
    end

    return res
end

function createWindow()
    baseLangdata = get(Ahorn.langdata, :settings)

    configData = Ahorn.config.data
    configLangdata = get(baseLangdata, :config)
    configOptions = getOptions(configData, Dict{String, Any}(), configLangdata)
    configSection = Ahorn.Form.Section("General", configOptions, dataName="config")

    debugData = Ahorn.debug.config.data
    debugLangdata = get(baseLangdata, :debug)
    debugOptions = getOptions(debugData, Dict{String, Any}(), debugLangdata)
    debugSection = Ahorn.Form.Section("Debug", debugOptions, dataName="debug")

    sections = Ahorn.Form.Section[
        configSection, debugSection
    ]

    function callback(data::Dict{String, Dict{String, Any}})
        merge!(Ahorn.config.data, data["config"])
        merge!(Ahorn.debug.config.data, data["debug"])

        Ahorn.saveConfig(Ahorn.config)
        Ahorn.saveConfig(Ahorn.debug.config)

        # Todo - Specially handle some values
        #  language
        #  etc
    end

    settingsWindow = Ahorn.Form.createFormWindow("$(Ahorn.baseTitle) - Settings", sections, callback=callback)

    return settingsWindow
end

function showSettings(widget::Union{Ahorn.MenuItemsTypes, Nothing}=nothing)
    if settingsWindow !== nothing
        Gtk.destroy(settingsWindow)
    end

    @Ahorn.catchall global settingsWindow = createWindow()

    showall(settingsWindow)
end

end


text:
include(joinpath("..", "ComplexOPF.jl","src", "PowSysMod_body.jl"))
include("SDP_decomposition_functions.jl")
include("solve_SDP.jl")
include("solve_minlp.jl")
include("B&B_fixingsome1and0.jl")

struct ROPF_infos
    instance_name::String
    matpower_instance_path::String
    output_instance_path::String
    decomposition::String
    output_decomposition_path::String
    generation::String
    generation_files_path::String
end

struct BB_infos
    search_strategy::String
    branch_strategy::String
    seuil_u::Float64
    seuil_l::Float64
end

function construct_dat_file_ROPF(instance_name, matpower_instance_path, output_path)
    typeofinput = MatpowerRTEROPFSimpleInput
    OPFpbs = load_OPFproblems(typeofinput, matpower_instance_path)
    # Bulding optimization problem
    pb_global = build_globalpb!(OPFpbs)
    pb_global_real = pb_cplx2real(pb_global)
    export_to_dat(pb_global_real, output_path, filename="$(instance_name).dat")
    return
end

function generate_clique_decomposition(instance_name, matpower_instance_path, output_decomposition_path)
    sp = read_sparsity_pattern(matpower_instance_path)
    nb_nodes = size(sp,1)
    nb_edges = (nnz(sp) - size(sp,1))/2
    L, order, nb_added_edges = chordal_ext_cholesky(sp)
    H = construct_graph_from_matrix(L)
    cliques_dict = cliques_max(H,order)
    GC = weighted_graph(cliques_dict)
    A = Prim_algo(GC)
    isdir(joinpath(output_decomposition_path, "cliquetree_cholesky")) || mkpath(joinpath(output_decomposition_path, "cliquetree_cholesky"))
    f = open(joinpath(output_decomposition_path, "cliquetree_cholesky", "$(instance_name)_sdp_cliquetree.txt"), "w")
        for edge in LightGraphs.edges(A)
            clique1 = src(edge)
            clique2 = dst(edge)
            @printf(f, "%10s %20s\n", "B$clique1", "B$clique2")
        end
    close(f)
    isdir(joinpath(output_decomposition_path, "blocks_cholesky")) || mkpath(joinpath(output_decomposition_path, "blocks_cholesky"))
    f = open(joinpath(output_decomposition_path, "blocks_cholesky", "$(instance_name)_sdp_blocks.txt"), "w")
    for (clique, nodes_list) in cliques_dict
        for node in nodes_list
            @printf(f, "%10s %20s\n", "B$clique", "VOLT_$(node)_Re")
            @printf(f, "%10s %20s\n", "B$clique", "VOLT_$(node)_Im")
        end
    end
    close(f)
end

function solve1(ROPF)
    LB_plus, stat_plus = solve_SDP(ROPF, "plus")
    LB_minus, stat_minus = solve_SDP(ROPF, "minus")
    UB_minus = UB_plus = Inf
    isdir("solutions") || mkpath("solutions")
    if !(stat_plus == MOI.FEASIBLE_POINT || stat_plus == MOI.NEARLY_FEASIBLE_POINT) && ! (stat_minus == MOI.FEASIBLE_POINT || stat_minus == MOI.NEARLY_FEASIBLE_POINT)
        println("$(ROPF.instance_name) $(ROPF.generation) : SDP relaxation probably infeasible ")
    end
    if (stat_plus == MOI.FEASIBLE_POINT || stat_plus == MOI.NEARLY_FEASIBLE_POINT)
        #solve MINLP with Knitro
        UB_plus = solve_minlp(ROPF, "plus", [], Dict{String,Float64}())
        mv("knitro_solution.csv", "solutions\\solve1_solution_$(ROPF.instance_name)_plus.csv", force=true)
    end
    if (stat_minus == MOI.FEASIBLE_POINT || stat_minus == MOI.NEARLY_FEASIBLE_POINT)
        #solve MINLP with Knitro
        UB_minus = solve_minlp(ROPF, "minus", [], Dict{String,Float64}())
        mv("knitro_solution.csv", "solutions\\solve1_solution_$(ROPF.instance_name)_minus.csv", force=true)
    end
    return UB_plus, LB_plus, UB_minus, LB_minus
end

function solve2(ROPF, max_time)
    LB_plus, stat_plus = solve_SDP(ROPF, "plus")
    LB_minus, stat_minus = solve_SDP(ROPF, "minus")
    UB_minus = UB_plus = Inf
    if !(stat_plus == MOI.FEASIBLE_POINT || stat_plus == MOI.NEARLY_FEASIBLE_POINT) && ! (stat_minus == MOI.FEASIBLE_POINT || stat_minus == MOI.NEARLY_FEASIBLE_POINT)
        println("$(ROPF.instance_name) $(ROPF.generation) : SDP relaxation probably infeasible ")
    end
    if (stat_plus == MOI.FEASIBLE_POINT || stat_plus == MOI.NEARLY_FEASIBLE_POINT)
        #B&B algo
        BB_parameters = BB_infos("deepfirst", "1", 0.9, 0.0001)
        (UB_plus, nb_nodes, open_nodes) = BandB_fixingsome1and0(ROPF, "plus", BB_parameters, max_time)
    end
    if (stat_minus == MOI.FEASIBLE_POINT || stat_minus == MOI.NEARLY_FEASIBLE_POINT)
        #B&B algo
        BB_parameters = BB_infos("deepfirst", "1", 0.9, 0.0001)
        (UB_minus, nb_nodes, open_nodes) = BandB_fixingsome1and0(ROPF, "minus", BB_parameters, max_time)
    end
    return UB_plus, LB_plus, UB_minus, LB_minus
end


text:
# this script does a Monte Carlo that shows efficiency of GLS compared to OLS
# when there are AR1 errors
using Econometrics, Plots

# the function to be Monte Carlo'ed
function wrapper(args)
	n, ρ = args    # sample size and AR parameter
	x = randn(n)     # an exogenous regressor
	# make the AR1 errors
    ϵ = zeros(n)
    ϵ[1] = randn()/sqrt(1-ρ^2)
	for t = 2:n
        ϵ[t] = ρ*ϵ[t-1] + randn()
    end
	# the dep var
	y = 1. .+ x .+ ϵ
	# OLS
	x = [ones(n) x]
	β_ols = x \ y
	# estimate ρ
	e = y - x*β_ols
    # estimating ρ using correlation, to keep in stationary region
    ρhat = cor(e[1:end-1], e[2:end])
	ystar = y - ρhat.*lag(y,1)
	xstar = x - ρhat*lag(x,1)
    ystar[1] = sqrt(1-ρhat^2) .*y[1]
    xstar[1,:] = sqrt(1-ρhat^2) .* x[1,:]
	β_gls = xstar\ystar
	β_ols, β_gls
end

# do the Monte Carlo
n = 30
ρ = 0.9
args = (n, ρ)
reps = 1000
β_ols = zeros(reps,2)
β_gls = zeros(reps,2)
for i = 1:reps
    β_ols[i,:], β_gls[i,:] = wrapper(args)
end

# analyze results
p1 = histogram(β_ols[:,2],bins=30, legend=false, xlabel="OLS")
p2 = histogram(β_gls[:,2],bins=30, legend=false, xlabel="GLS")
plot(p1,p2)
#savefig("AR1errors_OLSvsGLS.png")
gui()

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
"""
Bridges `CP.AbsoluteValue` to linear constraints and integer constraints.

The implemented model is the most generic one, so that the absolute value
always has a well-defined value. This requires the use of a binary variable.
In many cases, this is not necessary, and simpler models could be used, but
checking this automatically would require access to the whole model.

Based on Mosek's
[modelling cookbook](https://docs.mosek.com/modeling-cookbook/mio.html#exact-absolute-value).

YALMIP implements [a similar formulation, with big-M and small-m values
computed based on the bounds for each variable](https://github.com/yalmip/YALMIP/blob/develop/operators/alldifferent.m).
"""
struct AbsoluteValue2MILPBridge{T} <: MOIBC.AbstractBridge
    var_bin::MOI.VariableIndex
    var_bin_con::MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}
    var_pos::MOI.VariableIndex
    var_pos_con::MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{T}}
    var_neg::MOI.VariableIndex
    var_neg_con::MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{T}}

    con_original_var::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}
    con_abs_var::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}
    con_pos_var_big_m::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}
    con_neg_var_big_m::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}
end

function MOIBC.bridge_constraint(
    ::Type{AbsoluteValue2MILPBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.AbsoluteValue,
) where {T}
    return MOIBC.bridge_constraint(
        AbsoluteValue2MILPBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{AbsoluteValue2MILPBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.AbsoluteValue,
) where {T}
    f_scalars = MOIU.scalarize(f)

    # For this formulation work, both lower and upper bounds are required on
    # the argument of the absolute value.
    @assert CP.has_lower_bound(model, f_scalars[2])
    @assert CP.has_upper_bound(model, f_scalars[2])

    # Create the variables.
    var_bin, var_bin_con = MOI.add_constrained_variable(model, MOI.ZeroOne())
    var_pos, var_pos_con = MOI.add_constrained_variable(model, MOI.GreaterThan(zero(T)))
    var_neg, var_neg_con = MOI.add_constrained_variable(model, MOI.GreaterThan(zero(T)))

    # Relate the positive and negative parts to the set variables.
    con_original_var = MOI.add_constraint(
        model,
        f_scalars[2] - var_pos + var_neg,
        MOI.EqualTo(zero(T))
    )
    con_abs_var = MOI.add_constraint(
        model,
        f_scalars[1] - var_pos - var_neg,
        MOI.EqualTo(zero(T))
    )

    # Compute a bound on M.
    M = max(
        abs(CP.get_lower_bound(model, f_scalars[2])),
        abs(CP.get_upper_bound(model, f_scalars[2])),
    )

    # Big-M constraints.
    con_pos_var_big_m = MOI.add_constraint(
        model,
        var_pos - M * var_bin,
        MOI.LessThan(zero(T))
    )
    con_neg_var_big_m = MOI.add_constraint(
        model,
        var_neg - M * (1 - var_bin),
        MOI.LessThan(zero(T))
    )

    return AbsoluteValue2MILPBridge(
        var_bin, var_bin_con, var_pos, var_pos_con, var_neg, var_neg_con,
        con_original_var, con_abs_var, con_pos_var_big_m, con_neg_var_big_m,
    )
end

function MOI.supports_constraint(
    ::Type{AbsoluteValue2MILPBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.AbsoluteValue},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{AbsoluteValue2MILPBridge{T}}) where {T}
    return [(MOI.ZeroOne), (MOI.GreaterThan{T})]
end

function MOIB.added_constraint_types(::Type{AbsoluteValue2MILPBridge{T}}) where {T}
    return [
        (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
        (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}),
    ]
end

function MOI.get(::AbsoluteValue2MILPBridge, ::MOI.NumberOfVariables)
    return 3
end

function MOI.get(
    ::AbsoluteValue2MILPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return 1
end

function MOI.get(
    ::AbsoluteValue2MILPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VariableIndex, MOI.GreaterThan{T},
    },
) where {T}
    return 2
end

function MOI.get(
    ::AbsoluteValue2MILPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T},
    },
) where {T}
    return 2
end

function MOI.get(
    ::AbsoluteValue2MILPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.LessThan{T},
    },
) where {T}
    return 2
end

function MOI.get(
    b::AbsoluteValue2MILPBridge{T},
    ::MOI.ListOfVariableIndices,
) where {T}
    return [b.var_bin, b.var_pos, b.var_neg]
end

function MOI.get(
    b::AbsoluteValue2MILPBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return [b.var_bin_con]
end

function MOI.get(
    b::AbsoluteValue2MILPBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VariableIndex, MOI.GreaterThan{T},
    },
) where {T}
    return [b.var_pos_con, b.var_neg_con]
end

function MOI.get(
    b::AbsoluteValue2MILPBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T},
    },
) where {T}
    return [b.con_original_var, b.con_abs_var]
end

function MOI.get(
    b::AbsoluteValue2MILPBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.LessThan{T},
    },
) where {T}
    return [b.con_pos_var_big_m, b.con_neg_var_big_m]
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

using  Base.MultiplicativeInverses: SignedMultiplicativeInverse

struct ReshapedArray{T,N,P<:AbstractArray,MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int}}}} <: AbstractArray{T,N}
    parent::P
    dims::NTuple{N,Int}
    mi::MI
end
ReshapedArray{T,N}(parent::AbstractArray{T}, dims::NTuple{N,Int}, mi) = ReshapedArray{T,N,typeof(parent),typeof(mi)}(parent, dims, mi)

# IndexLinear ReshapedArray
ReshapedArrayLF{T,N,P<:AbstractArray} = ReshapedArray{T,N,P,Tuple{}}

# Fast iteration on ReshapedArrays: use the parent iterator
struct ReshapedArrayIterator{I,M}
    iter::I
    mi::NTuple{M,SignedMultiplicativeInverse{Int}}
end
ReshapedArrayIterator(A::ReshapedArray) = _rs_iterator(parent(A), A.mi)
function _rs_iterator{M}(P, mi::NTuple{M})
    iter = eachindex(P)
    ReshapedArrayIterator{typeof(iter),M}(iter, mi)
end

struct ReshapedIndex{T}
    parentindex::T
end

# eachindex(A::ReshapedArray) = ReshapedArrayIterator(A)  # TODO: uncomment this line
start(R::ReshapedArrayIterator) = start(R.iter)
@inline done(R::ReshapedArrayIterator, i) = done(R.iter, i)
@inline function next(R::ReshapedArrayIterator, i)
    item, inext = next(R.iter, i)
    ReshapedIndex(item), inext
end
length(R::ReshapedArrayIterator) = length(R.iter)

"""
    reshape(A, dims...)
    reshape(A, dims)

Return an array with the same data as the given array, but with different dimensions.

The new dimensions may be specified either as a list of arguments or as a shape
tuple. At most one dimension may be specified with a `:`, in which case its
length is computed such that its product with all the specified dimensions is
equal to the length of the original array A.

```jldoctest
julia> A = collect(1:16)
16-element Array{Int64,1}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16

julia> reshape(A, (4, 4))
4×4 Array{Int64,2}:
 1  5   9  13
 2  6  10  14
 3  7  11  15
 4  8  12  16

julia> reshape(A, 2, :)
2×8 Array{Int64,2}:
 1  3  5  7   9  11  13  15
 2  4  6  8  10  12  14  16
```

"""
reshape

reshape(parent::AbstractArray, dims::IntOrInd...) = reshape(parent, dims)
reshape(parent::AbstractArray, shp::NeedsShaping) = reshape(parent, to_shape(shp))
reshape(parent::AbstractArray, dims::Dims)        = _reshape(parent, dims)

# Allow missing dimensions with Colon():
reshape(parent::AbstractArray, dims::Int...) = reshape(parent, dims)
reshape(parent::AbstractArray, dims::Union{Int,Colon}...) = reshape(parent, dims)
reshape(parent::AbstractArray, dims::Tuple{Vararg{Union{Int,Colon}}}) = _reshape(parent, _reshape_uncolon(parent, dims))
# Recursively move dimensions to pre and post tuples, splitting on the Colon
@inline _reshape_uncolon(A, dims) = _reshape_uncolon(A, (), nothing, (), dims)
@inline _reshape_uncolon(A, pre, c::Void,  post, dims::Tuple{Any, Vararg{Any}}) =
    _reshape_uncolon(A, (pre..., dims[1]), c, post, tail(dims))
@inline _reshape_uncolon(A, pre, c::Void,  post, dims::Tuple{Colon, Vararg{Any}}) =
    _reshape_uncolon(A, pre, dims[1], post, tail(dims))
@inline _reshape_uncolon(A, pre, c::Colon, post, dims::Tuple{Any, Vararg{Any}}) =
    _reshape_uncolon(A, pre, c, (post..., dims[1]), tail(dims))
_reshape_uncolon(A, pre, c::Colon, post, dims::Tuple{Colon, Vararg{Any}}) =
    throw(DimensionMismatch("new dimensions $((pre..., c, post..., dims...)) may only have at most one omitted dimension specified by Colon()"))
@inline function _reshape_uncolon(A, pre, c::Colon, post, dims::Tuple{})
    sz, remainder = divrem(length(A), prod(pre)*prod(post))
    remainder == 0 || _throw_reshape_colon_dimmismatch(A, pre, post)
    (pre..., sz, post...)
end
_throw_reshape_colon_dimmismatch(A, pre, post) =
    throw(DimensionMismatch("array size $(length(A)) must be divisible by the product of the new dimensions $((pre..., :, post...))"))

reshape{T,N}(parent::AbstractArray{T,N}, ndims::Type{Val{N}}) = parent
function reshape{N}(parent::AbstractArray, ndims::Type{Val{N}})
    reshape(parent, rdims((), indices(parent), Val{N}))
end
# Move elements from inds to out until out reaches the desired
# dimensionality N, either filling with OneTo(1) or collapsing the
# product of trailing dims into the last element
@pure rdims{N}(out::NTuple{N,Any}, inds::Tuple{}, ::Type{Val{N}}) = out
@pure function rdims{N}(out::NTuple{N,Any}, inds::Tuple{Any, Vararg{Any}}, ::Type{Val{N}})
    l = length(last(out)) * prod(map(length, inds))
    (front(out)..., OneTo(l))
end
@pure rdims{N}(out::Tuple, inds::Tuple{}, ::Type{Val{N}}) = rdims((out..., OneTo(1)), (), Val{N})
@pure rdims{N}(out::Tuple, inds::Tuple{Any, Vararg{Any}}, ::Type{Val{N}}) = rdims((out..., first(inds)), tail(inds), Val{N})

# _reshape on Array returns an Array
_reshape(parent::Vector, dims::Dims{1}) = parent
_reshape(parent::Array, dims::Dims{1}) = reshape(parent, dims)
_reshape(parent::Array, dims::Dims) = reshape(parent, dims)

# When reshaping Vector->Vector, don't wrap with a ReshapedArray
function _reshape(v::AbstractVector, dims::Dims{1})
    len = dims[1]
    len == length(v) || throw(DimensionMismatch("parent has $(length(v)) elements, which is incompatible with length $len"))
    v
end
# General reshape
function _reshape(parent::AbstractArray, dims::Dims)
    n = _length(parent)
    prod(dims) == n || throw(DimensionMismatch("parent has $n elements, which is incompatible with size $dims"))
    __reshape((parent, IndexStyle(parent)), dims)
end

# Reshaping a ReshapedArray
_reshape(v::ReshapedArray{<:Any,1}, dims::Dims{1}) = _reshape(v.parent, dims)
_reshape(R::ReshapedArray, dims::Dims) = _reshape(R.parent, dims)

function __reshape(p::Tuple{AbstractArray,IndexCartesian}, dims::Dims)
    parent = p[1]
    strds = front(size_strides(parent))
    strds1 = map(s->max(1,s), strds)  # for resizing empty arrays
    mi = map(SignedMultiplicativeInverse, strds1)
    ReshapedArray(parent, dims, reverse(mi))
end

function __reshape(p::Tuple{AbstractArray,IndexLinear}, dims::Dims)
    parent = p[1]
    ReshapedArray(parent, dims, ())
end

@inline size_strides(A::AbstractArray) = tail(size_strides((1,), size(A)...))
size_strides(out::Tuple) = out
@inline size_strides(out, s, sz...) = size_strides((out..., out[end]*s), sz...)

size(A::ReshapedArray) = A.dims
similar(A::ReshapedArray, eltype::Type, dims::Dims) = similar(parent(A), eltype, dims)
IndexStyle(::Type{<:ReshapedArrayLF}) = IndexLinear()
parent(A::ReshapedArray) = A.parent
parentindexes(A::ReshapedArray) = map(s->1:s, size(parent(A)))
reinterpret{T}(::Type{T}, A::ReshapedArray, dims::Dims) = reinterpret(T, parent(A), dims)

@inline ind2sub_rs(::Tuple{}, i::Int) = i
@inline ind2sub_rs(strds, i) = ind2sub_rs((), strds, i-1)
@inline ind2sub_rs(out, ::Tuple{}, ind) = (ind+1, out...)
@inline function ind2sub_rs(out, strds, ind)
    d, r = divrem(ind, strds[1])
    ind2sub_rs((d+1, out...), tail(strds), r)
end

@inline function getindex(A::ReshapedArrayLF, index::Int)
    @boundscheck checkbounds(A, index)
    @inbounds ret = parent(A)[index]
    ret
end
@inline function getindex(A::ReshapedArray, indexes::Int...)
    @boundscheck checkbounds(A, indexes...)
    _unsafe_getindex(A, indexes...)
end
@inline function getindex(A::ReshapedArray, index::ReshapedIndex)
    @boundscheck checkbounds(parent(A), index.parentindex)
    @inbounds ret = parent(A)[index.parentindex]
    ret
end

@inline function _unsafe_getindex(A::ReshapedArray, indexes::Int...)
    @inbounds ret = parent(A)[ind2sub_rs(A.mi, sub2ind(size(A), indexes...))...]
    ret
end
@inline function _unsafe_getindex(A::ReshapedArrayLF, indexes::Int...)
    @inbounds ret = parent(A)[sub2ind(size(A), indexes...)]
    ret
end

@inline function setindex!(A::ReshapedArrayLF, val, index::Int)
    @boundscheck checkbounds(A, index)
    @inbounds parent(A)[index] = val
    val
end
@inline function setindex!(A::ReshapedArray, val, indexes::Int...)
    @boundscheck checkbounds(A, indexes...)
    _unsafe_setindex!(A, val, indexes...)
end
@inline function setindex!(A::ReshapedArray, val, index::ReshapedIndex)
    @boundscheck checkbounds(parent(A), index.parentindex)
    @inbounds parent(A)[index.parentindex] = val
    val
end

@inline function _unsafe_setindex!(A::ReshapedArray, val, indexes::Int...)
    @inbounds parent(A)[ind2sub_rs(A.mi, sub2ind(size(A), indexes...))...] = val
    val
end
@inline function _unsafe_setindex!(A::ReshapedArrayLF, val, indexes::Int...)
    @inbounds parent(A)[sub2ind(size(A), indexes...)] = val
    val
end

# helpful error message for a common failure case
ReshapedRange{T,N,A<:Range} = ReshapedArray{T,N,A,Tuple{}}
setindex!(A::ReshapedRange, val, index::Int) = _rs_setindex!_err()
setindex!(A::ReshapedRange, val, indexes::Int...) = _rs_setindex!_err()
setindex!(A::ReshapedRange, val, index::ReshapedIndex) = _rs_setindex!_err()

_rs_setindex!_err() = error("indexed assignment fails for a reshaped range; consider calling collect")

unsafe_convert{T}(::Type{Ptr{T}}, a::ReshapedArray{T}) = unsafe_convert(Ptr{T}, parent(a))

max_stars_repo_name: felixrehren/julia

max_issues_repo_name: felixrehren/julia

max_forks_repo_name: felixrehren/julia


text:
using JLD, Images

function main(dim)
    points = load("./out/20k.jld", "points")
    img = zeros(RGB,dim,dim)
    for c in points
        z = c
        while abs(z) < 2
            x = trunc(Int,(real(z) + 2)*dim/3); y = trunc(Int,(imag(z) + 1.5)*dim/3)
            0<x<dim>y>0 && (img[x,y] += RGB(0.0005,0,0))
            z = z^2 + c
        end
    end
    return map(clamp01nan, img)
end

@time save("./images/t1.png", main(1500))

max_stars_repo_name: BoardTM/fractals

max_issues_repo_name: BoardTM/fractals

max_forks_repo_name: BoardTM/fractals


text:
@testset "Logging" begin
    @testset "LocalEventLog" begin
        @testset "ArrayOp" begin
            ctx = Context()
            log = Dagger.LocalEventLog()
            ctx.log_sink = log

            sz = 2^4
            bsz = 2^2
            X = rand(Float32, sz, sz)
            XD = Distribute(Blocks(bsz, bsz), X)

            dag = XD*XD
            collect(ctx, dag)
            logs = Dagger.get_logs!(log)
            @test logs isa Vector{Dagger.Timespan}
            @test length(logs) > 0
            plan = Dagger.show_plan(logs)
            @test plan isa String
            @test occursin("digraph {", plan)
            @test occursin("Move:", plan)
            @test endswith(plan, "}\n")
        end

        @testset "Argument Merging" begin
            ctx = Context()
            log = Dagger.LocalEventLog()
            ctx.log_sink = log

            X = rand(Float32, 4, 3)
            a = delayed(sum)(X)
            b = delayed(sum)(X)
            c = delayed(+)(a,b)

            X1 = rand(Float32, 5, 10)
            X2 = rand(Float32, 5, 10)
            f = delayed(sum)(X1)
            g = delayed(sum)(X2)
            h = delayed(+)(f,g)

            j = delayed(+)(c,h)
            collect(ctx, j)
            logs = Dagger.get_logs!(log)
            plan = Dagger.show_plan(logs, j)
        end

        @testset "Automatic Plan Rendering" begin
            x = compute(rand(Blocks(2,2),4,4))
            mktemp() do path, io
                ctx = Context(;log_sink=Dagger.LocalEventLog(),log_file=path)
                compute(ctx, x * x)
                plan = String(read(io))
                @test occursin("digraph {", plan)
                @test occursin("Move:", plan)
                @test endswith(plan, "}\n")
            end
        end
    end
    @testset "MultiEventLog" begin
        ctx = Context()
        ml = Dagger.MultiEventLog()
        ml[:core] = Dagger.Events.CoreMetrics()
        ml[:id] = Dagger.Events.IDMetrics()
        ml[:timeline] = Dagger.Events.TimelineMetrics()
        ml[:wsat] = Dagger.Events.WorkerSaturation()
        ml[:loadavg] = Dagger.Events.CPULoadAverages()
        ml[:bytes] = Dagger.Events.BytesAllocd()
        ml[:mem] = Dagger.Events.MemoryFree()
        ml[:esat] = Dagger.Events.EventSaturation()
        ml[:psat] = Dagger.Events.ProcessorSaturation()
        ctx.log_sink = ml

        A = rand(Blocks(4, 4), 16, 16)
        collect(ctx, A*A)

        logs = Dagger.get_logs!(ml)
        for w in keys(logs)
            len = length(logs[w][:core])
            if w == 1
                @test len > 1
            end
            for c in (:core, :id, :timeline, :wsat, :loadavg, :bytes, :mem, :esat, :psat)
                @test haskey(logs[w], c)
                @test length(logs[w][c]) == len
            end
        end
        @test length(keys(logs)) > 1

        l1 = logs[1]
        core = l1[:core]
        @test !any(isnothing, core)
        esat = l1[:esat]
        @test any(e->haskey(e, :scheduler_init), esat)
        @test any(e->haskey(e, :schedule), esat)
        @test any(e->haskey(e, :fire), esat)
        @test any(e->haskey(e, :take), esat)
        @test any(e->haskey(e, :finish), esat)
        # Note: May one day be true as scheduler evolves
        @test !any(e->haskey(e, :compute), esat)
        @test !any(e->haskey(e, :move), esat)
        psat = l1[:psat]
        # Note: May become false
        @test all(e->length(e) == 0, psat)

        had_psat_proc = 0
        for wo in filter(w->w != 1, keys(logs))
            lo = logs[wo]
            esat = lo[:esat]
            @test !any(e->haskey(e, :scheduler_init), esat)
            @test !any(e->haskey(e, :schedule), esat)
            @test !any(e->haskey(e, :fire), esat)
            @test !any(e->haskey(e, :take), esat)
            @test !any(e->haskey(e, :finish), esat)
            psat = lo[:psat]
            if any(e->length(e) > 0, psat)
                had_psat_proc += 1
                @test any(e->haskey(e, :compute), esat)
                @test any(e->haskey(e, :move), esat)
            end
        end
        @test had_psat_proc > 0

        logs = Dagger.get_logs!(ml)
        for w in keys(logs)
            for c in keys(logs[w])
                @test isempty(logs[w][c])
            end
        end
    end
end


text:
"""
Module exports functions based on [documentation structure](https://developer.spotify.com/documentation/web-api/reference/)

    using Spotify
    # See all functions :
    varinfo(Spotify; recursive=true)
    Or use tab_complation after typing

    Spotify.+tab
    Tracks.+tab

"""
module Spotify
using HTTP, Parameters, Dates, IniFile, HTTP, StructTypes, JSON3, URIs, Sockets
using HTTP: handle
import Dates
import Base64
using Base64: base64encode
import Base: show, show_vector, typeinfo_implicit
export authorize, strip_embed_code
"For the client credentials flow"
const AUTH_URL = "https://accounts.spotify.com/api/token"
"For the 'Implicit grant flow"
const OAUTH_AUTHORIZE_URL = "https://accounts.spotify.com/authorize"
const DEFAULT_REDIRECT_URI = "http://127.0.0.1:8080/api"
const NOT_ACTUAL = "Paste_32_bytes_in_inifile"
"A list of potentially available browsers, to be tried in succession if present"
const BROWSERS = ["firefox", "edge", "chrome", "safari", "phantomjs"]

@with_kw mutable struct SpotifyCredentials
    client_id::String = ""
    client_secret::String = ""
    encoded_credentials::String = ""
    access_token::String = ""
    token_type::String = ""
    expires_at::String = string(Dates.DateTime(0))
    redirect::String = ""
    ig_access_token::String = ""
end

const SPOTCRED  = Ref{SpotifyCredentials}(SpotifyCredentials())

spotcred() = SPOTCRED[]

# Helper types - colorful aliases for 'String'.
include("types.jl")

# Authorization
include("authorization/client_credentials.jl")
include("authorization/open_a_browser.jl")
include("authorization/implicit_grant_flow.jl")

# Util
include("util/utilities.jl")
include("util/request.jl")

# Browse
include("by_reference_doc/browse.jl")

# Tracks
include("by_reference_doc/tracks.jl")

# Shows
include("by_reference_doc/shows.jl")

# Albums
include("by_reference_doc/albums.jl")

# Artist
include("by_reference_doc/artist.jl")

# Episodes
include("by_reference_doc/episodes.jl")

# Library
include("by_reference_doc/library.jl")

# Personalization
include("by_reference_doc/personalization.jl")

# Follow
include("by_reference_doc/follow.jl")

# Export structure
include("export_structure.jl")

function __init__()
    authorize()
end
end # module


text:
"""
    Wing(; foils, chords, twists,
           spans, dihedrals, LE_sweeps,
           position, angle, axis)

    Wing(left :: HalfWing, right :: HalfWing)
    Wing(half :: HalfWing)

A composite type consisting of fields `left` and `right`, each a `HalfWing` for constructing a wing.

A single argument generates a symmetric `Wing`, viz. `left === right`.

# Arguments
- `foils :: Vector{Foil}`
- `chords`
- `twists`
- `spans`
- `dihedrals`
- `LE_sweeps `
- `position = zeros(3) `
- `angle    = 0 `
- `axis     = [0.,1.,0.] `

"""
struct Wing{M <: AbstractWing} <: AbstractWing
    left  :: M
    right :: M
    # function Wing(left :: T, right :: T) where T <: AbstractWing
    #     @assert first(left.chords) == first(right.chords) && first(left.twists) == first(right.twists) "Properties of first section must match."
    #     # Add more conditions...
    #     Wing{T}(left, right)
    # end
end

# Getters
left(wing :: Wing)      = wing.left
right(wing :: Wing)     = wing.right
foils(wing :: Wing)     = @views [ (foils ∘ left)(wing)[end:-1:2]     ; (foils ∘ right)(wing)     ]
chords(wing :: Wing)    = @views [ (chords ∘ left)(wing)[end:-1:2]    ; (chords ∘ right)(wing)    ]
twists(wing :: Wing)    = @views [ (twists ∘ left)(wing)[end:-1:2]    ; (twists ∘ right)(wing)    ]
spans(wing :: Wing)     = @views [ (reverse ∘ spans ∘ left)(wing)     ; (spans ∘ right)(wing)     ]
dihedrals(wing :: Wing) = @views [ (reverse ∘ dihedrals ∘ left)(wing) ; (dihedrals ∘ right)(wing) ]

sweeps(wing :: Wing, w = 0.) = [ reverse(sweeps(left(wing), w)); sweeps(right(wing), w) ]

Base.position(wing :: Wing) = right(wing).affline.translation
orientation(wing :: Wing)   = right(wing).affine.linear

affine_transformation(wing :: Wing) = affine_transformation(right(wing))

# Symmetric wings
Wing(bing :: HalfWing) = Wing(bing, bing)

Wing(;
        chords,
        twists    = zero(chords),
        spans     = ones(length(chords) - 1) / (length(chords) - 1),
        dihedrals = zero(spans),
        sweeps    = zero(spans),
        w_sweep   = 0.0,
        foils     = fill(naca4(0,0,1,2), length(chords)),
        position  = zeros(3),
        angle     = 0.,
        axis      = SVector(1., 0., 0.)
    ) = Wing(HalfWing(foils = foils, chords = chords, twists = twists, spans = spans / 2, dihedrals = dihedrals, sweeps = sweeps, w_sweep = w_sweep, position = position, angle = angle, axis = axis))

"""
    WingSection(; span, dihedral, sweep, taper, root_chord,
                  root_twist, tip_twist, root_foil, tip_foil,
                  position, angle, axis)

Define a `Wing` consisting of two trapezoidal sections with ``y``-``z`` reflection symmetry.

# Arguments
- `span       :: Real         = 1.`: Span length
- `dihedral   :: Real         = 1.`: Dihedral angle (degrees)
- `LE_sweep   :: Real         = 0.`: Leading-edge sweep angle (degrees)
- `taper      :: Real         = 1.`: Taper ratio of tip to root chord
- `root_chord :: Real         = 1.`: Root chord length
- `root_twist :: Real         = 0.`: Twist angle at root (degrees)
- `tip_twist  :: Real         = 0.`: Twist angle at tip (degrees)
- `root_foil  :: Foil         = naca4((0,0,1,2))`: Foil coordinates at root
- `tip_foil   :: Foil         = root_foil`: Foil coordinates at tip
- `position   :: Vector{Real} = zeros(3)`: Position
- `angle      :: Real         = 0.`: Angle of rotation (degrees)
- `axis       :: Vector{Real} = [0.,1.,0.]`: Axis of rotation
"""
WingSection(;
        span       = 1.,
        dihedral   = 0.,
        sweep      = 0.,
        w_sweep    = 0.,
        taper      = 1.,
        root_chord = 1.,
        root_twist = 0.,
        tip_twist  = 0.,
        root_foil  = naca4((0,0,1,2)),
        tip_foil   = root_foil,
        position   = zeros(3),
        angle      = 0.,
        axis       = SVector(1., 0., 0.)
    ) = Wing(HalfWingSection(span = span / 2, dihedral = dihedral, sweep = sweep, w_sweep = w_sweep, taper = taper, root_chord = root_chord, root_twist = root_twist, tip_twist = tip_twist, root_foil = root_foil, tip_foil = tip_foil, position = position, angle = angle, axis = axis))

span(wing :: Wing) = (span ∘ left)(wing) + (span ∘ right)(wing)

taper_ratio(wing :: Wing) = (taper_ratio ∘ left)(wing), (taper_ratio ∘ right)(wing)

projected_area(wing :: Wing) = (projected_area ∘ left)(wing) + (projected_area ∘ right)(wing)

mean_aerodynamic_chord(wing :: Wing) = ((mean_aerodynamic_chord ∘ left)(wing) + (mean_aerodynamic_chord ∘ right)(wing)) / 2

"""
    wing_bounds(wing :: Wing)

Return the leading and trailing edge coordinates of a `Wing`.
"""
function wing_bounds(wing :: Wing)
    left_lead, left_trail   = wing_bounds(left(wing), true)
    right_lead, right_trail = wing_bounds(right(wing))

    leading  = @views [ left_lead[1:end-1,:] ; right_lead  ]
    trailing = @views [ left_trail[1:end-1,:]; right_trail ]

    leading, trailing
end

leading_edge(wing :: Wing)  = @views [ leading_edge(left(wing), true)[1:end-1,:] ; leading_edge(right(wing))  ]
trailing_edge(wing :: Wing) = @views [ trailing_edge(left(wing), true)[1:end-1,:]; trailing_edge(right(wing)) ]

reflect_xz(vector) = SVector(vector[1], -vector[2], vector[3])
mean_aerodynamic_center(wing :: Wing, factor = 0.25) = (mean_aerodynamic_center(right(wing), factor) .+ reflect_xz(mean_aerodynamic_center(left(wing), factor))) ./ 2


text:
module MeshingBenchmarks

using FileIO
using NRRD
using Meshing
using MeshIO
using GeometryTypes
using BenchmarkTools

function benchmark()
    here = dirname(@__FILE__)

    println(pwd())
    println("CTA-cardio.nrrd loading...")
    ctacardio = @btime load($here*"/../data/CTA-cardio.nrrd")
    q = 100

    samples = 1

    println("CTA-cardio.nrrd loaded")

    println("CTA-cardio.nrrd MarchingCubes Float32 runtime")
    mc = @btime HomogenousMesh{Point{3,Float32}, Face{3,Int}}($ctacardio, MarchingCubes(iso=Float32($q), insidepositive=true))
    for i in 1:samples-1
        @btime HomogenousMesh{Point{3,Float32}, Face{3,Int}}($ctacardio, MarchingCubes(iso=Float32($q), insidepositive=true))
    end

    println("CTA-cardio.nrrd MarchingCubes Float64 runtime")
    @btime HomogenousMesh{Point{3,Float64}, Face{3,Int}}($ctacardio, MarchingCubes(iso=$q, insidepositive=true))
    for i in 1:samples-1
        @btime HomogenousMesh{Point{3,Float64}, Face{3,Int}}($ctacardio, MarchingCubes(iso=$q, insidepositive=true))
    end

    # println("CTA-cardio.nrrd MarchingTetrahedra Float32 runtime")
    # mt = @btime HomogenousMesh{Point{3,Float32},Face{3,Int}}($ctasdf, MarchingTetrahedra($q))
    # for i in 1:samples-1
    #     @btime HomogenousMesh{Point{3,Float32},Face{3,Int}}($ctasdf, MarchingTetrahedra($q))
    # end

    println("Saving files")
    save("ctacardio_mc.ply", mc)
    #save("ctacardio_mt.ply", mt)
end

function brain(samples=1, q=100)

    here = dirname(@__FILE__)

    brain = load(here*"/../Seg3DData/Brain_DataSet/MRI-brain50.nrrd")

    brainsdf = SignedDistanceField(HyperRectangle(Vec(0,0,0), Vec(10,10,10)),brain.data)

    println("Brain.nrrd MarchingCubes Float32 runtime")
    mc_brain = @btime HomogenousMesh{Point{3,Float32},Face{3,Int}}($brainsdf, MarchingCubes(iso=$q, insidepositive=false))
    for i in 1:samples-1
        @btime HomogenousMesh{Point{3,Float32},Face{3,Int}}($brainsdf, MarchingCubes(iso=$q, insidepositive=false))
    end

    save("brian_mc.ply", mc_brain)
end

end # module


text:
using Test

@testset "iseven isodd" begin
    # check sign bits
    @test iseven(1f0,-8)
    @test isodd(-1f0,-8)
    @test iseven(1.0,-11)
    @test isodd(-1.0,-11)
    @test iseven(Float16(1),-5)
    @test isodd(Float16(-1),-5)

    @test isodd(1.5f0,1)
    @test isodd(1.5,1)
    @test isodd(Float16(1.5),1)

    @test iseven(1.25f0,1)
    @test iseven(1.25,1)
    @test iseven(Float16(1.25),1)

    @test isodd(1.25f0,2)
    @test isodd(1.25,2)
    @test isodd(Float16(1.25),2)
end

@testset "Zero rounds to zero" begin
    for T in [Float16,Float32,Float64]
        for k in -5:50
            A = zeros(T,2,3)
            Ar = round(A,k)
            @test A == Ar
            @test zero(T) == round(zero(T),k)
        end
    end
end

@testset "one rounds to one" begin
    for T in [Float16,Float32,Float64]
        for k in 0:50
            A = ones(T,2,3)
            Ar = round(A,k)
            @test A == Ar
            @test one(T) == round(one(T),k)
        end
    end
end

@testset "minus one rounds to minus one" begin
    for T in [Float16,Float32,Float64]
        for k in 0:50
            A = -ones(T,2,3)
            Ar = round(A,k)
            @test A == Ar
            @test -one(T) == round(-one(T),k)
        end
    end
end

@testset "No rounding for keepbits=10,23,52" begin
    for (T,k) in zip([Float16,Float32,Float64],
                        [11,24,53])
        A = rand(T,200,300)
        Ar = round(A,k)
        @test A == Ar

        # and a single one
        r = rand(T)
        @test r == round(r,k)
    end
end

@testset "Approx equal for keepbits=5,10,25" begin
    for (T,k) in zip([Float16,Float32,Float64],
                        [5,10,25])
        A = rand(T,200,300)
        Ar = round(A,k)
        @test A ≈ Ar

        # and a single one
        r = rand(T)
        @test r ≈ round(r,k)
    end
end

@testset "Idempotence" begin
    for T in [Float16,Float32,Float64]
        for k in 0:20
            A = rand(T,200,300)
            Ar = round(A,k)
            Ar2 = round(Ar,k)
            @test Ar == Ar2
        end
    end
end

@testset "Tie to even" begin

    for T in [Float16,Float32,Float64]

        @test round(1.5f0,0) == T(2)

        @test round(1.25f0,1) == T(1)
        @test round(1.5f0,1) == T(1.5)
        @test round(1.75f0,1) == T(2)

        @test round(1.125f0,2) == T(1)
        @test round(1.375f0,2) == T(1.5)
        @test round(1.625f0,2) == T(1.5)
        @test round(1.875f0,2) == T(2)
    end

    for k in 1:10
        m = 0x8000_0000 >> (9+k)
        x = reinterpret(UInt32,one(Float32)) + m
        x = reinterpret(Float32,x)
        @test 1f0 == round(x,k)
    end

    # for T in [Float16,Float32,Float64]
    #     for k in 1:20
    #         x = randn(T)
    #         xr1 = round(x,k+1)
    #         xr = round(xr1,k)

    #         if iseven(xr,k)

    #     end
    # end
end

@testset "Round to nearest?" begin
    N = 1000
    for _ in 1:N
        for (T,UIntT) in zip([Float16,Float32,Float64],
                                [UInt16,UInt32,UInt64])
            for k in 1:9
                x = randn(T)
                xr = round(x,k)

                ulp = Base.sign_mask(T) >> (Base.exponent_bits(T)+k)
                next_xr = reinterpret(T,reinterpret(UIntT,xr) + ulp)
                prev_xr = reinterpret(T,reinterpret(UIntT,xr) - ulp)

                @test abs(next_xr - x) >= abs(xr - x)
                @test abs(prev_xr - x) >= abs(xr - x)
            end
        end
    end
end


text:
module DetrendUtilities

	function winsorize(data, trim=1)
		M, N, J, T = size(data)

		nu_guess = copy(data)
		for m = 1:M, j = 1:J, t = 1:T
			vector = nu_guess[m,:,j,t]
			smallest = sort(vector[vector .> 0])[1+trim]
			largest = sort(vector[vector .< 1])[end-trim]
			nu_guess[m,vector .< smallest,j,t] = smallest
			nu_guess[m,vector .> largest,j,t] = largest
		end
	    return nu_guess
	end

	function detrend(X, weights::Array{Float64,1})
		(ndims(X) == 2 || ndims(X) == 4) || error("The function 'detrend' can only be applied to 2 or 4-dimensional arrays")

		dims = ndims(X)

		if dims == 2
			I, S = size(X)
			J, K = 1, 1
		elseif dims == 4
			I, J, K, S = size(X)
			X = reshape(X, (I*J*K, S))
		end

		X_c = zeros(I*J*K,S)
		X_t = zeros(I*J*K,S)

		# Maximum lag
		P = size(weights,1) - 1
		w = [flipdim(weights,1); weights[2:end]]

		# Extended series with P 0's at both ends
		X_ext = zeros(I*J*K, S + 2P)
		X_ext[:,(P + 1):(P + S)] = X

		# Reflect head/tail
		X_ext[:,1:P] = repeat(X[:,1], outer = [1,P]) + flipdim(repeat(X[:,1], outer = [1,P]) - X[:,2:(P + 1)],2)
		X_ext[:,(S + P + 1):(S + 2P)] = repeat(X[:,S], outer = [1,P]) + flipdim(repeat(X[:,S],outer = [1,P]) - X[:,(S - P):(S - 1)],2)

		for s in 1:S
			X_c[:,s] = X_ext[:,s:(s + 2P)] * w
			X_t[:,s] = X_ext[:,s + P] - X_c[:,s]
		end

		if dims == 4
			X_c = reshape(X_c, (I, J, K, S))
			X_t = reshape(X_t, (I, J, K, S))
		end

		return X_c, X_t
	end
end


text:
function u_test(x, y, σ)
    exp( - ( (x-π)^2 + (y-π)^2 ) / (2σ^2) )
end

function testFourierQuadrature()
    M = 64
    N = 64

    σ = 0.5

    grid = Grid2d(M,N)
    ft   = FourierTransform(grid)
    u    = get_field(grid)
    û    = get_trans(ft)

    for i in 1:M
        for j in 1:N
            u[i,j]  = u_test(grid.x[i], grid.y[j], σ)
        end
    end

    prfft!(ft, u, û)
    ûint = fourier_quadrature(û, û, ft.μ, grid)

    @test real(ûint) - π*σ^2 ≈ zero(eltype(u)) atol=1E-15
    @test imag(ûint) ≈ zero(eltype(u)) atol=eps()
end

testFourierQuadrature()


text:
mutable struct MetaDigraph{T} <: AbstractMetaNexus{T}
    core::Digraph{T}
    nprops::Dict{T,Dict{Symbol,Any}}
    aprops::Dict{Arc{T},Dict{Symbol,Any}}
end

@forward MetaDigraph.core badj
MetaDigraph{T}() where {T} = MetaDigraph(Digraph{T}())
function MetaDigraph(core::Digraph{T}) where {T}
    MetaDigraph{T}(core, Dict{T,Dict{Symbol,Any}}(), Dict{Arc{T},Dict{Symbol,Any}}())
end
is_directed(::MetaDigraph) = true
Digraph(g::MetaDigraph{T}) where {T} = g.core

function set_props!(g::MetaDigraph{T}, a::Arc{T}, d::Dict) where {T}
    has_arc(g, a) || return false
    !_hasdict(g, a) ? g.aprops[a] = d : merge!(g.aprops[a], d)
    return true
end

rem_prop!(g::MetaDigraph{T}, a::Arc{T}, prop::Symbol) where {T} =
    delete!(g.aprops[a], prop)

clear_props!(g::MetaDigraph{T}, a::Arc{T}) where {T} =
    _hasdict(g, a) && delete!(g.aprops, a)


text:
#=
""" `construct_gp_prior()`
Returns a GP prior using a Matern 5/2 kernel and specified parameters
Optional Inputs:
- smooth_factor: Multiplies variance and length scale for GP kernel (1)
- σ²: Variance for GP kernel (0.5)
- l: length scale for GP kernel (5.8e-5)
"""
function construct_gp_prior(; smooth_factor::Real = 1, σ²::Real = 0.5, l::Real = 5.8e-5)
	σ² *= smooth_factor
	l *= smooth_factor
	k = σ² * transform(Matern52Kernel(), 1 / l)
	f_naive = GP(k)
	#f = to_sde(f_naive)   # if develop issues with StaticArrays could revert to this
	f = to_sde(f_naive, SArrayStorage(Float64))
end
=#

""" `construct_gp_prior()`
Returns a GP prior using a Matern 5/2 kernel and specified parameters
Optional Inputs:
- smooth_factor_σ²: Multiplies variance for GP kernel (1)
- smooth_factor_l: Multiplies length scale for GP kernel (1)
- σ²: Variance for GP kernel (0.5)
- l: length scale for GP kernel (5.8e-5)
"""
function construct_gp_prior(; smooth_factor_σ²::Real = 1, smooth_factor_l::Real = 1, smooth_factor::Real = 1, σ²::Real = 0.5, l::Real = 5.8e-5)
	σ² *= smooth_factor_σ² * smooth_factor
	l *= smooth_factor_l * smooth_factor
	k = σ² * transform(Matern52Kernel(), 1 / l)
	f_naive = GP(k)
	#f = to_sde(f_naive)   # if develop issues with StaticArrays could revert to this
	f = to_sde(f_naive, SArrayStorage(Float64))
end

#=
function construct_gp_prior_kernel_param_as_vector(gp_param = [0.5, 5.8e-5 ])
	σ², l = gp_param
	k = σ² * transform(Matern52Kernel(), 1 / l)
	f_naive = GP(k)
	#f = to_sde(f_naive)   # if develop issues with StaticArrays could revert to this
	f = to_sde(f_naive, SArrayStorage(Float64))
end
=#

""" `construct_gp_posterior(xobs, yobs; sigmasq_obs, use_logx, use_logy, smooth_factor, boost_factor )`
Inputs:
- xobs: x locations where data is provided
- yobs: y values of data to condition on
Optional Inputs:
- sigmasq_obs: variances for y values being conditioned on
- use_logx: If true, take log's of x values before fitting GP
- use_logy: If true, perform log transform on y's
- smooth_factor: scales GP hyperparameters so as to result in smoother GP posterior (1)
- boost_factor: scales xobs by 1/boost_factor (1)
Returns:
- Posterior GP at locations xpred given training data
"""
function construct_gp_posterior(xobs::AA1, yobs::AA2; sigmasq_obs::AA4 = 1e-16*ones(length(xobs)),
			 	use_logx::Bool = true, #= use_logy::Bool = true, =# smooth_factor::Real = 1, boost_factor::Real = 1 )   where {
				T1<:Real, AA1<:AbstractArray{T1,1}, T2<:Real, AA2<:AbstractArray{T2,1}, T3<:Real, AA3<:AbstractArray{T3,1}, T4<:Real, AA4<:AbstractArray{T4,1} }
	f = construct_gp_prior( smooth_factor=smooth_factor )
	if boost_factor != one(boost_factor)
		xobs_trans = use_logx ? log.(xobs./boost_factor) : xobs./boost_factor
	else
		xobs_trans = use_logx ? log.(xobs) : xobs
	end
	fx = f(xobs_trans, sigmasq_obs)
	f_posterior = posterior(fx, yobs )
    #yobs_trans = use_logy ? log.(yobs) : yobs
    #sigmasq_obs_trans = use_logy ? sigmasq_obs./yobs.^2 : sigmasq_obs
	#fx = f(xobs_trans, sigmasq_obs_trans)
	#f_posterior = posterior(fx, yobs_trans )
end

""" predict_gp_mean_var(gp, xpred ; use_logx, use_logy)
Inputs:
- gp:
- xpred: Locations to predict GP at
Optional inputs:
- use_logx: If true, apply log transform to xpred before evaluating GP
- use_logy: If true, apply exp transform after evaluating GP
Returns tuple with vector of means and variance of GP posterior at locations in xpred.
"""
function predict_gp_mean_var(gp::AGP,  xpred::AA3, use_logx::Bool = true, use_logy::Bool = false )   where {
		AGP<:AbstractGP, T3<:Real, AA3<:AbstractArray{T3,1} }
		gp_post = use_log_x ? gp(log.(xpred)) : gp(xpred)
		m = mean.(gp_post)
		v = var.(gp_post)
		if use_logy
			m = exp.(m)
			v .*= m.^2
		end
		return m, v
end

""" predict_gp_mean(gp, xpred ; use_logx, use_logy)
Inputs:
- gp:
- xpred: Locations to predict GP at
Optional inputs:
- use_logx: If true, apply log transform to xpred before evaluating GP
Returns vector of means of GP posterior at locations in xpred.
"""
function predict_mean(gp::AGP )   where { AGP<:Fx_PosteriorType, T3<:Real, AA3<:AbstractArray{T3,1} }
	gp_post = marginals(gp)
	m = mean.(gp_post)
	return m
end

function predict_deriv(gp::AGP, xpred::AA3; use_logx::Bool = true )   where { AGP<:Fx_PosteriorType,
				T3<:Real, AA3<:AbstractArray{T3,1}  }
  #kobs = make_kernel_data(xobs, kernel=kernel, sigmasq_obs=sigmasq_obs, sigmasq_cor=sigmasq_cor, rho=rho)
  #kobs_pred_deriv = make_kernel_obs_pred(xobs,xpred, kernel=dkerneldx, sigmasq_cor=sigmasq_cor, rho=rho)
  #alpha = kobs \ yobs
  #pred_deriv = kobs_pred_deriv' * alpha
  xpred_trans = use_logx ? log.(xpred) : xpred
  m = predict_mean(gp) #,xpred_trans)
  #if use_logy   m .= exp(m)   end
  return numerical_deriv(xpred_trans,m)
  #=
  dfluxdlnλ = zeros(size(m))
  dfluxdlnλ[1] = (m[2]-m[1])/(xpred_trans[2]-xpred_trans[1])
  dfluxdlnλ[2:end-1] .= (m[3:end].-m[1:end-2])./(xpred_trans[3:end].-xpred_trans[1:end-2]) # exp.(m[2:end-1]).*
  dfluxdlnλ[end] = (m[end]-m[end-1])/(xpred_trans[end]-xpred_trans[end-1])
  return dfluxdlnλ
	=#
end

function predict_deriv2(gp::AGP, xpred::AA3; use_logx::Bool = true )   where { AGP<:Fx_PosteriorType,
				T3<:Real, AA3<:AbstractArray{T3,1}  }
	#kobs = make_kernel_data(xobs, kernel=kernel, sigmasq_obs=sigmasq_obs, sigmasq_cor=sigmasq_cor, rho=rho)
	#kobs_pred_deriv2 = make_kernel_obs_pred(xobs,xpred, kernel=d2kerneldx2, sigmasq_cor=sigmasq_cor, rho=rho)
	#alpha = kobs \ yobs
	#pred_deriv = kobs_pred_deriv2' * alpha

	xpred_trans = use_logx ? log.(xpred) : xpred
	m = predict_mean(gp)
	#if use_logy   m .= exp(m)   end
	d2fluxdlnλ2 = zeros(size(m))
	d2fluxdlnλ2[2:end-1] .= (m[3:end].+m[1:end-2].-2.0.*m[2:end-1])./(xpred_trans[3:end].-xpred_trans[1:end-2]).^2
	d2fluxdlnλ2[1] = d2fluxdlnλ2[2]
	d2fluxdlnλ2[end] = d2fluxdlnλ2[end-1]

	# TODO add some sort of check based on second derivative and dx that the finite difference is good enough?
	return d2fluxdlnλ2
end

function predict_mean(xobs::AA1, yobs::AA2, xpred::AA3;	sigmasq_obs::AA4 = 1e-16*ones(length(xobs)),
			use_logx::Bool = true, use_logy::Bool = false, smooth_factor::Real = 1, boost_factor::Real = 1, verbose::Bool = false  )   where {
				T1<:Real, AA1<:AbstractArray{T1,1}, T2<:Real, AA2<:AbstractArray{T2,1}, T3<:Real, AA3<:AbstractArray{T3,1}, T4<:Real, AA4<:AbstractArray{T4,1} }
	# global ncalls += 1
	@assert size(xobs) == size(yobs) == size(sigmasq_obs)
  	#println("# predict_mean (TemporalGPs): size(xobs) = ",size(xobs), "  size(xpred) = ", size(xpred))
	if verbose    println("# Entering predict_mean for TemporalGPs")    end
	tstart = now()
	if verbose    println("# extrema(yobs) = ", extrema(yobs) )   end
	mean_y = mean(yobs)
	if verbose    println("# mean_y = ", mean_y )    end
	y_trans = yobs
	sigmasq_obs_trans = sigmasq_obs
	if use_logy
		#println("# mean_y = ", mean_y, " min_y = ", minimum(yobs), " max_y = ", maximum(yobs))
		y_trans = log.(yobs./mean_y)
		sigmasq_obs_trans = sigmasq_obs./yobs.^2
	else
		y_trans = yobs./mean_y .- 1
		sigmasq_obs_trans = sigmasq_obs./mean_y.^2
	end
	if verbose    println("# extrema(xobs) = ", extrema(xobs) )    end
	if verbose    println("# extrema(y_trans) = ", extrema(y_trans) )    end
	if verbose    println("# extrema(sigmasq_obs_trans) = ", extrema(sigmasq_obs_trans) )   end
	if verbose    println("# smooth_factor = ", smooth_factor, " boost_factor = ",boost_factor )    end
	f_posterior = construct_gp_posterior(xobs,y_trans,sigmasq_obs=sigmasq_obs_trans, use_logx=use_logx, smooth_factor=smooth_factor, boost_factor=boost_factor )
	#println("typeof(f_posterior) = ",typeof(f_posterior))
	#println("f_posterior <: AbstractGP = ",typeof(f_posterior) <: AbstractGP )
	#println("f_posterior <: AbstractMvNormal = ",typeof(f_posterior) <: AbstractMvNormal )
	xpred_trans = use_logx ? log.(xpred) : xpred
	if verbose    println("# extrema(xpred_trans) = ", extrema(xpred_trans) )    end
	fx_posterior = f_posterior(xpred_trans)
	#println("typeof(fx_posterior) = ",typeof(fx_posterior))
	#println("f_posterior <: AbstractGP = ",typeof(fx_posterior) <: AbstractGP )
	#println("f_posterior <: AbstractMvNormal = ",typeof(fx_posterior) <: AbstractMvNormal )
	#println("f_posterior <: Fx_PosteriorType = ",typeof(fx_posterior) <: Fx_PosteriorType )
	#output = predict_mean(f_posterior(xpred_trans), xpred_trans ) #, use_logx=use_logx,use_logy=use_logy)
	pred_mean = predict_mean(fx_posterior)
	if verbose    println("# extrema(pred_mean) = ", extrema(pred_mean) )    end
	if use_logy
		pred_mean .= exp.(pred_mean).*mean_y
	else
		pred_mean .= (pred_mean.+1).*mean_y
	end
	if verbose    println("# extrema(pred_mean) = ", extrema(pred_mean) )    end
	#println("# predict_mean (TemporalGPs) runtime: ", now()-tstart)
	return pred_mean
end

# NEED TO TEST
function predict_deriv(xobs::AA1, yobs::AA2, xpred::AA3; sigmasq_obs::AA4 = 1e-16*ones(length(xobs))
						, use_logx::Bool = true, use_logy::Bool = false, smooth_factor::Real = 1, verbose::Bool = false  )    where { T1<:Real, AA1<:AbstractArray{T1,1}, T2<:Real, AA2<:AbstractArray{T2,1}, T3<:Real, AA3<:AbstractArray{T3,1}, T4<:Real, AA4<:AbstractArray{T4,1}  }
  #kobs = make_kernel_data(xobs, kernel=kernel, sigmasq_obs=sigmasq_obs, sigmasq_cor=sigmasq_cor, rho=rho)
  #kobs_pred_deriv = make_kernel_obs_pred(xobs,xpred, kernel=dkerneldx, sigmasq_cor=sigmasq_cor, rho=rho)
  #alpha = kobs \ yobs
  #pred_deriv = kobs_pred_deriv' * alpha
  if verbose   println("# predict_deriv (TemporalGPs): size(xobs) = ",size(xobs), "  size(xpred) = ", size(xpred))   end
  pred_mean = predict_mean(xobs,yobs,xpred, sigmasq_obs=sigmasq_obs, use_logx=use_logx, use_logy=use_logy, smooth_factor=smooth_factor)
  pred_deriv = calc_dfluxdlnlambda(pred_mean,xpred)
  return pred_deriv
  end

# NEED TO TEST
function predict_deriv2(xobs::AA1, yobs::AA2, xpred::AA3;sigmasq_obs::AA4 = 1e-16*ones(length(xobs))
						, use_logx::Bool = true, use_logy::Bool = false, smooth_factor::Real = 1  )   where { T1<:Real, AA1<:AbstractArray{T1,1}, T2<:Real, AA2<:AbstractArray{T2,1}, T3<:Real, AA3<:AbstractArray{T3,1}, T4<:Real, AA4<:AbstractArray{T4,1}  }
  #kobs = make_kernel_data(xobs, kernel=kernel, sigmasq_obs=sigmasq_obs, sigmasq_cor=sigmasq_cor, rho=rho)
  #kobs_pred_deriv2 = make_kernel_obs_pred(xobs,xpred, kernel=d2kerneldx2, sigmasq_cor=sigmasq_cor, rho=rho)
  #alpha = kobs \ yobs
  #pred_deriv = kobs_pred_deriv2' * alpha
println("# predict_deriv2 (TemporalGPs) runtime: ", now()-tstart)
  pred_mean = predict_mean(xobs,yobs,xpred, sigmasq_obs=sigmasq_obs, use_logx=use_logx, use_logy=use_logy, smooth_factor=smooth_factor)
	pred_deriv2 = calc_d2fluxdlnlambda2(pred_mean,xpred)
	return pred_deriv2
end

function predict_mean_and_deriv(xobs::AA1, yobs::AA2, xpred::AA3;sigmasq_obs::AA4 = 1e-16*ones(length(xobs))
	 							, use_logx::Bool = true, use_logy::Bool = false, smooth_factor::Real = 1, verbose::Bool = false ) where { T1<:Real, AA1<:AbstractArray{T1,1}, T2<:Real, AA2<:AbstractArray{T2,1}, T3<:Real, AA3<:AbstractArray{T3,1}, T4<:Real, AA4<:AbstractArray{T4,1}  }
	#		kernel::Function = matern52_sparse_kernel, dkerneldx::Function = dkerneldx_matern52_sparse,
	#		sigmasq_cor::Real=1.0, rho::Real=1.0)
  #kobs = make_kernel_data(xobs, kernel=kernel, sigmasq_obs=sigmasq_obs, sigmasq_cor=sigmasq_cor, rho=rho)
  #alpha = kobs \ yobs
  #kobs_pred = make_kernel_obs_pred(xobs,xpred, kernel=kernel, sigmasq_cor=sigmasq_cor, rho=rho)
  #pred_mean = kobs_pred' * alpha
  #kobs_pred_deriv = make_kernel_obs_pred(xobs,xpred, kernel=dkerneldx, sigmasq_cor=sigmasq_cor, rho=rho)
  #pred_deriv = kobs_pred_deriv' * alpha
  if verbose   println("# predict_mean_and_deriv (TemporalGPs): size(xobs) = ",size(xobs), "  size(xpred) = ", size(xpred))   end
  xpred_trans = use_logx ? log.(xpred) : xpred
  pred_mean = predict_mean(xobs,yobs,xpred, sigmasq_obs=sigmasq_obs, use_logx=use_logx, use_logy=use_logy, smooth_factor=smooth_factor)
  pred_deriv = numerical_deriv(xpred_trans,pred_mean)
  return (mean=pred_mean, deriv=pred_deriv)
end

function predict_mean_and_derivs(xobs::AA1, yobs::AA2, xpred::AA3; sigmasq_obs::AA4 = 1e-16*ones(length(xobs))
								, use_logx::Bool = true, use_logy::Bool = false, smooth_factor::Real = 1  ) where { T1<:Real, AA1<:AbstractArray{T1,1}, T2<:Real, AA2<:AbstractArray{T2,1}, T3<:Real, AA3<:AbstractArray{T3,1}, T4<:Real, AA4<:AbstractArray{T4,1}  }
	#		kernel::Function = matern52_sparse_kernel, dkerneldx::Function = dkerneldx_matern52_sparse, d2kerneldx2::Function = d2kerneldx2_matern52_sparse,
	#		,	sigmasq_cor::Real=1.0, rho::Real=1.0)
  #=
  kobs = make_kernel_data(xobs, kernel=kernel, sigmasq_obs=sigmasq_obs, sigmasq_cor=sigmasq_cor, rho=rho)
  alpha = kobs \ yobs
  kobs_pred = make_kernel_obs_pred(xobs,xpred, kernel=kernel, sigmasq_cor=sigmasq_cor, rho=rho)
  pred_mean = kobs_pred' * alpha
  kobs_pred_deriv = make_kernel_obs_pred(xobs,xpred, kernel=dkerneldx, sigmasq_cor=sigmasq_cor, rho=rho)
  pred_deriv = kobs_pred_deriv' * alpha
  kobs_pred_deriv2 = make_kernel_obs_pred(xobs,xpred, kernel=d2kerneldx2, sigmasq_cor=sigmasq_cor, rho=rho)
  pred_deriv2 = kobs_pred_deriv2' * alpha
  =#
  #println("# predict_mean_and_derivs (TemporalGPs): size(xobs) = ",size(xobs), "  size(xpred) = ", size(xpred))
  pred_mean = predict_mean(xobs,yobs,xpred, sigmasq_obs=sigmasq_obs, use_logx=use_logx, use_logy=use_logy, smooth_factor=smooth_factor)
  pred_deriv = calc_dfluxdlnlambda(pred_mean,xpred)
  pred_deriv2 = calc_d2fluxdlnlambda2(pred_mean,xpred)
  return (mean=pred_mean, deriv=pred_deriv, deriv2=pred_deriv2)
end

function log_pdf_gp_posterior(xobs::AA, yobs::AA #, kernel::Function;
			; sigmasq_obs::AA = 1e-16*ones(length(xobs)),
			use_logx::Bool = true, use_logy::Bool = false, smooth_factor::Real = 1  )  where { T<:Real, AA<:AbstractArray{T,1} }
  	#kobs = make_kernel_data(xobs, kernel=kernel, sigmasq_obs=sigmasq_obs, sigmasq_cor=sigmasq_cor, rho=rho)
  	# -0.5*( invquad(kobs, yobs) + logdet(kobs) + length(xobs)*log(2pi) )
	mean_y = mean(yobs)
	xobs_trans = use_logx ? log.(xobs) : xobs
	if use_logy
		yobs_trans = log.(yobs./mean_y)
		sigmasq_obs_trans = sigmasq_obs./yobs.^2
	else
		yobs_trans = yobs./mean_y .- 1
		sigmasq_obs_trans = sigmasq_obs./mean_y.^2
	end
	f_posterior = construct_gp_posterior(xobs_trans,yobs_trans,xobs_trans,sigmasq_obs=sigmasq_obs_trans, use_logx=false, smooth_factor=smooth_factor )
	return -logpdf(f_posterior(xobs_trans, sigmasq_obs_trans, use_logx=false), yobs_trans)
end

function fit_blaze_order(blaze::AbstractArray{T1,1}, order::Integer; mask::AbstractArray{Bool,1} = trues(length(blaze)),
						smooth_factor_σ²::Real = 1, smooth_factor_l::Real = 1 ) where { T1<:Real }
   gp_pr = construct_gp_prior(smooth_factor_σ²=smooth_factor_σ², smooth_factor_l=smooth_factor_l)
   x_obs = (1:length(blaze))[mask]
   num_pix_in_mask = length(x_obs)
   sigmasq_obs = var(blaze[mask][floor(Int,num_pix_in_mask//4):ceil(Int,num_pix_in_mask*3//4)])
   println("σ²_obs = ", sigmasq_obs)
   fx = gp_pr(x_obs, sigmasq_obs)
   f_post = posterior(fx, blaze[x_obs])
   cost = logpdf(f_post, blaze[x_obs])
   println("cost = ", cost)
   gp_mean = marginals(f_post)
   mean.(gp_post)
end


text:
# output = wx + b
immutable Affine{T,P<:Params} <: Learnable
    nin::Int
    nout::Int
    input::SumNode{T,1}
    output::OutputNode{T,1}
    params::P
end

function Affine{T}(::Type{T}, nin::Int, nout::Int,
                    θ::AbstractVector = zeros(T, nout*(nin+1)),
                    ∇::AbstractVector = zeros(T, nout*(nin+1)))
    input = InputNode(T, nin)
    output = OutputNode(T, nout)
    params = Params(θ, ∇, ((nout,nin), (nout,)))
    w, b = params.views
    initialize_weights!(w)
    initialize_bias!(b)
    Affine(nin, nout, input, output, params)
end
Affine(nin::Int, nout::Int, args...) = Affine(Float64, nin, nout, args...)

function Base.show(io::IO, t::Affine)
    print(io, "Affine{$(t.nin)-->$(t.nout)}")
end

# compute output = wx + b
function transform!(t::Affine)
    w, b = t.params.views
    x = t.input.val
    y = t.output.val
    copy!(y, b)
    for o=1:t.nout
        y[o] += sum(w[o,i] * x[i] for i=1:t.nin)
    end
    y
end

# update the partial derivatives:
#   ∇x = ∂L/∂x
#   ∇w = ∂L/∂w
#   ∇b = ∂L/∂b
# use the chain rule, assuming that we've already updated ∇out = ∂L/∂y
function grad!(t::Affine)
    w, b = t.params.views
    ∇w, ∇b = t.params.∇_views
    x = t.input.val
    ∇x = t.input.∇
    ∇y = t.output.∇

    # ∇x, ∇w
    for i=1:t.nin
        ∇xᵢ = zero(eltype(∇x))
        for o=1:t.nout
            ∇xᵢ += w[o,i] * ∇y[o]
            ∇w[o,i] = x[i] * ∇y[o]
        end
        ∇x[i] = ∇xᵢ
    end

    # ∇b
    copy!(∇b, ∇y)
    return
end


text:
"""
    dc_EGM!(m::FModel,p::Param)

DCEGM algorithm as in Ishkakov et al.
"""
function dc_EGM!(m::FModel,p::Param; demo = false)

    for it in p.nT:-1:1
        # println(it)
        # @info("period = $it")

        if it==p.nT
            for id in 1:p.nD   # work of dont work
                # final period: consume everyting.
                m.c[id,it] = Envelope(MLine(vcat(p.a_lowT,p.a_high),vcat(0.0,p.a_high)) )
                # initialize value function with bequest function
                m.v[id,it] = p.ν > 0 ?  Envelope(MLine(m.avec, bequest(m.avec,id == 1,p) )) : Envelope(MLine(vcat(p.a_lowT,p.a_high),vcat(0.0,NaN)))
            end
        else
            for id in 1:p.nD   # current period dchoice
                working = id==1  # working today is id=1
                # clamp!(mm1, p.cfloor, Inf)  # dont actullay want to do that: want to allow for neg assets
                mm1 = m.m1[it][id]
                # next period consumption and values y-coords
                # for each d-choice
                cmat = fill(-Inf,p.nD,p.na*p.ny)
                vmat = fill(-Inf,p.nD,p.na*p.ny)

                # get future cons and values
                # for future nD compute consumption and value function
                for iid in 1:p.nD
                    # what's next period's cash on hand in this discrete choice?
                    # mm1 = m.m1[it+1][iid]

                    c1 = interp(m.c[iid,it+1].env, mm1[:])
                    floory!(c1,p.cfloor)   # floor negative consumption
                    cmat[iid,:] = gety(c1)  # get y-values: next period consumption
                    vmat[iid,:] = vfun(iid,it+1,cmat[iid,:],mm1[:],m.v[iid,it+1],p)
                end

                # get ccp to be a worker
                pwork = working ? ccp(vmat,p) : p.delta*ones(size(vmat)[2])

                if it == p.nT-1 && p.ν > 0
                    # expected marginal utility of bequest on RHS
                    cmat[:] = cmat .+ p.bbar  # add luxury level of bequest before we evaluate in marginal utility
                    mu1 = reshape(pwork .* up(cmat[1,:],p) .+ (1.0 .- pwork) .* up(cmat[2,:],p),p.ny,p.na)
                    RHS = p.beta * p.R * p.ν * m.ywgt' * mu1
                else
                    # get expected marginal utility of consumption t+1
                    mu1 = reshape(pwork .* up(cmat[1,:],p) .+ (1.0 .- pwork) .* up(cmat[2,:],p),p.ny,p.na)

                    # get expected marginal value of saving: RHS of euler equation
                    # beta * R * E[ u'(c_{t+1}) ]
                    # ywgt: weight on the income states
                    RHS = p.beta * p.R * m.ywgt' * mu1
                end

                # optimal consumption today: invert the RHS of euler equation
                c0 = iup(Array(RHS)[:],p)

                # set optimal consumption function today. endo grid m and cons c0
                cline = MLine(m.avec .+ c0, c0)
                # store
                m.c[id,it] = Envelope(cline)
                # consumption function done.

                # compute value function
                # ----------------------
                if working
                    ev =  m.ywgt' * reshape(logsum(vmat,p),p.ny,p.na)
                else
                    ev =  (1-p.delta) * m.ywgt' * reshape(vmat[2,:],p.ny,p.na) + p.delta * m.ywgt' * reshape(logsum(vmat,p),p.ny,p.na)
                end
                vline = MLine(m.avec .+ c0, u(c0,id==1,p) .+ p.beta * ev[:])
                if demo
                    m.vdirty[id,it] = Envelope(vline)
                    m.cdirty[id,it] = Envelope(cline)
                end

                # SECONDARY ENVELOPE COMPUTATION
                # ==============================

                if id==1   # only for workers
                    m.v[id,it], m.c[id,it] = do_secondary(vline,cline,working,ev[1],p)
                else   # if id==1
                    m.v[id,it] = Envelope(vline)
                    m.c[id,it] = Envelope(cline)
                end

                # store the expected value at the lower boundary
                # in a separate object
                m.v[id,it].vbound = ev[1]

                # this creates the credit constrained region
                prepend!(m.c[id,it].env,[Point(m.avec[1],0.0)])
                prepend!(m.v[id,it].env,[Point(m.avec[1],ev[1])])
                sortx!(m.c[id,it].env)  # sort cons by default
            end # current discrete choice
        end    # if final period
    end     # loop over time
end

"""
    dc_EGM!(m::GModel,p::Param)

DCEGM algorithm for a model with state dependence.
"""
function dc_EGM!(m::GModel,p::Param)

    cmat = fill(-Inf,p.nD,p.na)
    vmat = fill(-Inf,p.nD,p.na)
    ctmp = fill(-Inf,p.nD,p.ny,p.na)
    # vtmp = fill(-Inf,p.nD,p.ny,p.na)

    for it in p.nT:-1:1
        is_bequest = (it == p.nT-1) && (p.ν > 0)

        if it==p.nT
            for id in 1:p.nD
                for iy in 1:p.ny
                    # final period: consume everyting.
                    m.c[id,iy,it] = Envelope(MLine(vcat(p.a_lowT,p.a_high),vcat(0.0,p.a_high)) )
                    # initialize value function with vf(1) = 0
                    # the ((id > 1)*10.0 is here to make sure that working is never optimal when dead.
                    m.v[id,iy,it] = p.ν > 0 ?  Envelope(MLine(m.avec, bequest(m.avec,id == 1,p) .- ((id > 1)*10.0) )) : Envelope(MLine(vcat(p.a_lowT,p.a_high),vcat(0.0,NaN)))
                end
            end
        else
            for iy in 1:p.ny  # current state
                # next period consumption and values y-coords
                # for each d-choice
                # reset all value matrices
                fill!(vmat,0.0)
                fill!(ctmp,-Inf)
                fill!(cmat,0.0)

                # All future uncertainty to be accounted for here
                for jy in 1:p.ny # future state: owner, renter, income, etc
                    pr = m.ywgt[iy,jy]  # transprob
                    for iid in 1:p.nD  # future dchoice
                        m1 = m.m1[it+1][iid][jy,:]  # state specific mvec
                        c1 = interp(m.c[iid,jy,it+1].env, m1) # C(d',y',m')
                        floory!(c1,p.cfloor)   # floor negative consumption
                        ctmp[iid,jy,:] = gety(c1)
                        vmat[iid,:] += pr * vfun(iid,it+1,ctmp[iid,jy,:],m1,m.v[iid,jy,it+1],p)
                        if is_bequest
                            cmat[iid,:] += pr * up(ctmp[iid,jy,:] .+ p.bbar,p)
                        else
                            cmat[iid,:] += pr * up(ctmp[iid,jy,:],p)
                        end
                    end
                end # end future state

                for id in 1:p.nD  # current dchoice
                    working = id==1  # working today is id=1

                    # get ccp of choices: P(d'|iy), pwork
                    pwork = working ? ccp(vmat,p) : zeros(size(vmat)[2])

                    # get expected marginal utility of consumption t+1
                    mu1 = pwork .* cmat[1,:] .+ (1.0 .- pwork) .* cmat[2,:] # 1,na
                    #RHS
                    RHS = is_bequest ? p.beta * p.R * p.ν * mu1 : p.beta * p.R * mu1

                    #optimal cons
                    c0 = iup(RHS,p)
                    # set optimal consumption function today. endo grid m and cons c0
                    cline = MLine(m.avec .+ c0, c0)
                    # consumption function done.

                    # compute value function
                    # ----------------------
                    if working
                        ev =  logsum(vmat,p)
                    else
                        ev =  vmat[2,:]
                    end
                    vline = MLine(m.avec .+ c0, u(c0,id==1,p) .+ p.beta * ev[:])

                    # SECONDARY ENVELOPE COMPUTATION
                    # ==============================

                    if id==1   # only for workers
                        m.v[id,iy,it], m.c[id,iy,it] = do_secondary(vline,cline,working,ev[1],p)
                    else   # if id==1
                        m.v[id,iy,it] = Envelope(vline)
                        m.c[id,iy,it] = Envelope(cline)
                    end

                    # store the expected value at the lower boundary
                    # in a separate object
                    m.v[id,iy,it].vbound = ev[1]

                    # this creates the credit constrained region
                    prepend!(m.c[id,iy,it].env,[Point(m.avec[1],0.0)])
                    prepend!(m.v[id,iy,it].env,[Point(m.avec[1],ev[1])])
                    sortx!(m.c[id,iy,it].env)  # sort cons by default
                end  # current id
            end  # iy
        end # if last period
    end# it
end

"""
    vfun(id::Int,it::Int,c1::Vector{Float64},m1::Vector{Float64},en::Array{Envelope},p::Param)

Calculate the period `it`, discrete choice `id`-specific value function. Avoids interpolation in credit constrained region by using the analytic form of the value function (no need to interpolate expected value function when on the lower bound of assets.)
"""
function vfun(id::Int,it::Int,c1::Vector{Float64},m1::Vector{Float64},v::Envelope,p::Param)

    # L = en.L[id]
    # v = en[id,iy,it]

    # computes v_{it}(m) = u(c) + beta v_{it+1}(m1)

    if length(getx(v.env)) < 2
        error("need more than 2 points in envelope object")
    end

    r = fill(NaN,size(m1))

    if (it == p.nT) && (p.ν > 0)
        r[:] = bequest(c1,false, p)  # working == true
    else
        mask = m1.<getx(v.env)[2]
        mask = it==p.nT ? trues(size(mask)) : mask

        if all(mask)
            # in the credit constrained region:
            r[:] = u(c1,id==1,p) .+ p.beta * bound(v)
        elseif any(mask)
            r[mask] = u(c1[mask],id==1,p) .+ p.beta * bound(v)
            # elsewhere
            r[.!mask] = gety(interp(v.env,m1[.!mask]))
        else
            r[:] = gety(interp(v.env,m1))
        end
    end
    return r
end

function vfun(ufun::Function,it::Int,c1::Vector{Float64},m1::Vector{Float64},v::Envelope,p::Param)

    # L = en.L[id]
    # v = en[id,iy,it]

    # computes v_{it}(m) = u(c) + beta v_{it+1}(m1)

    if length(getx(v.env)) < 2
        error("need more than 2 points in envelope object")
    end

    r = fill(NaN,size(m1))
    mask = m1.<getx(v.env)[2]
    mask = it==p.nT ? trues(size(mask)) : mask

    if all(mask)
        # in the credit constrained region:
        r[:] = ufun(c1) .+ p.beta * bound(v)
    elseif any(mask)
        r[mask] = ufun(c1[mask]) .+ p.beta * bound(v)
        # elsewhere
        r[.!mask] = gety(interp(v.env,m1[.!mask]))
    else
        r[:] = gety(interp(v.env,m1))
    end

    return r
end

"""
vfun used in simulation of state dependent model
"""
function vfun(ufun::Function,it::Int,c1::Vector{Float64},m1::Vector{Float64},vv::Vector{Envelope},p::Param)

    # each entry of c1,m1 and v represents an individual
    N = length(c1)
    @assert N == length(m1) == length(vv)
    # L = en.L[id]
    # v = en[id,iy,it]

    # computes v_{it}(m) = u(c) + beta v_{it+1}(m1)

    r = fill(NaN,N)
    critx = map(x -> getx(x.env)[2],vv)
    for (i,v) in enumerate(vv)
        if length(getx(v.env)) < 2
            error("need more than 2 points in envelope object")
        end
        if m1[i] < critx[i]  # credit constrained
        # if m1[i] < getx(v.env)[2]  # credit constrained
            r[i] = ufun(c1[i]) + p.beta * bound(v)
        else
            r[i] = gety(interp(v.env,[m1[i]]))[1]
        end
    end
    r
end

"""
    ccp(x::Matrix,p::Param)

Conditional Choice probability of working
"""
function ccp(x::Matrix,p::Param)
    #choice probability of the first row in 2-row matrix
    mx = maximum(x,dims = 1)
    mxx = x.-repeat(mx,size(x)[1],1)   # center values at max for numerical stability
    vec(exp.(mxx[1,:]./p.lambda)' ./ sum(exp.(mxx./p.lambda),dims = 1))
end

"""
    logsum(x::Matrix,p::Param)

Logsum of conditional values used in Expected value function.
"""
function logsum(x::Matrix,p::Param)
    mx = maximum(x,dims = 1)
    mxx = x .- repeat(mx,size(x)[1],1)
    mx .+ p.lambda * log.( sum(exp.(mxx./p.lambda), dims = 1) )
end

function runf(;par=Dict())
    p = Param(par=par)
    # p.beta = 1/p.R
    m = FModel(p)
    dc_EGM!(m,p)
    (m,p)
end
function runfdemo(;par=Dict())
    p = Param(par=par)
    # p.beta = 1/p.R
    m = FModel(p)
    dc_EGM!(m,p,demo = true)
    (m,p)
end
function runfp(;par = Dict())
    m,p = runf(par = par)
    plot(m,p)
end

function rung(;par=Dict())
    p = Param(par=par)
    # p.beta = 1/p.R
    m = GModel(p)
    dc_EGM!(m,p)
    (m,p)
end

function run2()
    p = Param()
    m = GModel(p)
    dc_EGM!(m,p)
    (m,p)
    # plot(m,id=2,xlim=(0,4),ylim=(-100,-20))
    # plot!(m,id=1,linestyle=:dash,label="")
end
function runit()
    p = Param()
    m = GModel(p)
    dc_EGM!(m,p)
    return (m,p)
end
function pp()
    m,p = runit()
    plot(m.v[1,1].env)
end

"""
    do_secondary(vline::MLine, cline::MLine, working::Bool, ev0::Float64, p::Param)

Performs secondary envelope cleaning on the DC algo example with working true/false.

1. Checks whether there is any backward bending of calculated EGM value function beyond first grid point in x0
2. Inserts a new point if so (from analytic expression of utility)
3. Calls [`secondary_envelope`](@ref)
4. Prunes invalid points from the consumption policy function and adds new intersections ([`secondary_envelope`](@ref) is only pruning the value function)
"""
function do_secondary(vline::MLine, cline::MLine, working::Bool, ev0::Float64, p::Param)
    minx = min_x(vline)
    if vline.v[1].x <= minx
        # normal case - no bending back behind first point
        # call secondary_envelope on vline
        vout = secondary_envelope(vline)
        cout = Envelope(cline)

    else
        # non-convex region lies inside credit constraint.
        # endogenous x grid bends back before the first x grid point.
        # println("minx = $minx")
        # println("first point = $(vline.v[1].x)")
        x0 = collect(range(minx,stop = vline.v[1].x,length = floor(Integer,p.na/10))) # some points to the left of first x point
        x0 = x0[1:end-1]
        c0 = copy(x0)
        if minx < 0
            c0[:] .= c0 .+ abs(minx) .+ p.cfloor
        end
        y0 = u(c0,working,p) .+ p.beta .* ev0   # use c0: positive consumption even with neg assets!

        prepend!(vline,convert(Point,x0,y0))
        prepend!(cline,convert(Point,x0,c0))  # cons policy in credit constrained is 45 degree line
        cout = Envelope(cline)
        vout = secondary_envelope(vline)
    end

    # now we have a cleaned value function in vout

    # if any points were removed from vline, now we need to clean the consumption function as well
    if length(getr(vout)) > 0

        # analyse policy function
        # =======================

        rmidx = getr(vout)  # indices removed from value function
        insert_left = Point[]
        insert_right = Point[]

        # remove illegal points from c
        # deleteat!(m.c[id,iy,it].env.v, rmidx)

        # insert new intersections into consumption function
        isecs = gets(vout)
        consx = getx(cout.env)

        # compute any new cons-points (intersections in v)
        # ---------------------------------------

        for isec in 1:length(isecs)
            I = isecs[isec]

            # interpolate from left
            jl = findall(consx .< I.x)
            jl = jl[ jl .∉ Ref(rmidx) ]  # keep those who are not to be deleted
            if length(jl) > 0
                jl = maximum(jl)  # biggest of those
                if jl < length(consx)
                    newleft = MLine(cout.env.v[jl:jl+1])
                    # @fediskhakov: who guarantees that jl+1 is not to be deleted?
                    # more generally: why do we not delete points before we
                    # insert the new intersections?

                    sortx!(newleft)
                    tmp = getv(interp(newleft, [ I.x ] ))[1]
                    # take this new point at a minimal left shift in x
                    push!(insert_left, Point(tmp.x - 1e3 * eps(), tmp.y) )
                # else
                #     println("trying to insert after last point")
                end
            # else
            #     push!(insert_left,I)
            end

            # interpolate from right
            jr = findall(consx .> I.x)
            jr = jr[ jr .∉ Ref(rmidx) ]  # keep those who are not to be deleted
            if length(jr) > 0
                jr = minimum(jr)   # smallest of those
                # push!(insert_right, interp(m.c[id,iy,it].env[jr-1:jr], [ I.x ] ) )
                newright = MLine(cout.env.v[jr-1:jr])
                sortx!(newright)
                push!(insert_right, getv(interp(newright, [ I.x ] ))[1] )
            # else
            #     push!(insert_right,I)
            end
        end # all intersections

        # remove illegal points from c
        deleteat!(cout.env.v, rmidx)

        # add new points in twice with a slight offset from left
        # to preserve the ordering in x.
        for ileft in 1:length(insert_left)

            consx = getx(cout.env)
            j = findfirst(consx .> insert_left[ileft].x)  # first point past new intersection
            insert!(cout.env.v,j,insert_left[ileft])  # item is j-th index
            insert!(cout.env.v,j+1,insert_right[ileft])
        end
    end
    (vout, cout)
end


text:
"""
    ARTMAP.jl

Description:
    Includes all of the ARTMAP (i.e., explicitly supervised) ART modules definitions.
"""

include("common.jl")    # train!
include("DAM.jl")       # Default ARTMAP
include("FAM.jl")       # Fuzzy ARTMAP
include("SFAM.jl")      # Simplified Fuzzy ARTMAP
include("ARTSCENE.jl")  # ARTSCENE filters


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 0b19c578-762b-11eb-34b6-01e80cef1406
using Pkg, DrWatson

# ╔═╡ 3f814e9e-762b-11eb-1340-91617ca7b58a
begin
    using Distributions
    using StatsPlots
    using StatsBase
    using LaTeXStrings
    using CSV
    using DataFrames
    using LinearAlgebra
    using Random
    using StanSample
    using StatisticalRethinking
end

# ╔═╡ 3faf5b5e-762b-11eb-3a66-41ec8a3030ae
stan9_0 ="
data {
    int N;
    vector[N] x;
    vector[N] y;
}
parameters {
    real mux;
    real muy;
}
model {
    mux ~ normal(0, 0.5);
    muy ~ normal(0, 0.5);
    x ~ normal(mux, 1);
    y ~ normal(muy, 1);
}
";

# ╔═╡ 3fafc7fe-762b-11eb-37da-254b6d75f5ca
begin
	m9_0s = SampleModel("m9.0s", stan9_0)
	data = (N = 100, x = rand(Normal(), 100), y = rand(Normal(), 100))
	rc9_0s = stan_sample(m9_0s; data)

	if success(rc9_0s)
		m9_0s_df = read_samples(m9_0s, :dataframe)
		PRECIS(m9_0s_df)
	end
end

# ╔═╡ 3fccf39e-762b-11eb-0f4d-c120f0c58ffc
begin
	# test data
	#Random.seed!(123)

	y = zscore(rand(Normal(), 50))
	x = zscore(rand(Normal(), 50))
end;

# ╔═╡ 697559da-7ac6-11eb-03b2-373578580a2a
begin
    current_q = [-0.1, 0.2]
    pr = 0.5
    eps = 0.03
    L = 11
    n_samples = 4
end

# ╔═╡ 69b1fcaa-7ac6-11eb-1c91-591ff5af66e6
begin
	Random.seed!(159)
    res = Vector{NamedTuple}(undef, n_samples)
    local q = current_q
    for i in 1:n_samples
        res[i] = hmc(x, y, u, ugrad, eps, L, q)
        q = res[i].q
    end
end;

# ╔═╡ 69bde2a6-7ac6-11eb-314f-f9f9ba80de24
begin
    fig = plot(xlims=[-pr, pr], ylims=[-pr, pr], leg=false)
	k0 = zeros(L, n_samples)
	for i in 1:L
		for j in 1:n_samples
			k0[i, j] = sum(res[j].ptraj[i, :].^2)/2		# kinetic energy
		end
	end

    for i in 1:n_samples
        xpos = res[i].qtraj[:, 1]
        ypos = res[i].qtraj[:, 2]
		for j in 1:L
        	plot!([xpos[j], xpos[j+1]], [ypos[j], ypos[j+1]], color=:lightgrey,
				linewidth=k0[j, i])
		end
		scatter!(xpos, ypos, color=:grey)
		i == 1 && plot!([xpos[1]], [ypos[1]], marker = ([:x :d], 10, 0.9, :black))
		annotate!([(xpos[end] - 0.005, ypos[end] - 0.03,
			Plots.text(i, 10, :red, :right))])
    end
	fig
end

# ╔═╡ 8aef9a52-7aea-11eb-0435-e37fa4d80369
begin
    res50 = Vector{NamedTuple}(undef, 50)
    local q = current_q
    for i in 1:50
        res50[i] = hmc(x, y, u, ugrad, eps, L, q)
        q = res50[i].q
    end
end;

# ╔═╡ b669674e-7aea-11eb-1549-b969f18740e1
begin
   fig50 = plot(xlims=[-pr, pr], ylims=[-pr, pr], leg=false)
   for i in 1:50
		xpos = [res50[i].qtraj[end, 1]]
		ypos = [res50[i].qtraj[end, 2]]
		scatter!(xpos, ypos, color=res50[i].accept==1 ? :grey : :red)
		i == 1 && plot!([xpos[1]], [ypos[1]], marker = ([:x :d], 10, 0.9, :blue))
    end
	fig50
end

# ╔═╡ Cell order:
# ╠═0b19c578-762b-11eb-34b6-01e80cef1406
# ╠═3f814e9e-762b-11eb-1340-91617ca7b58a
# ╠═3faf5b5e-762b-11eb-3a66-41ec8a3030ae
# ╠═3fafc7fe-762b-11eb-37da-254b6d75f5ca
# ╠═3fccf39e-762b-11eb-0f4d-c120f0c58ffc
# ╠═697559da-7ac6-11eb-03b2-373578580a2a
# ╠═69b1fcaa-7ac6-11eb-1c91-591ff5af66e6
# ╠═69bde2a6-7ac6-11eb-314f-f9f9ba80de24
# ╠═8aef9a52-7aea-11eb-0435-e37fa4d80369
# ╠═b669674e-7aea-11eb-1549-b969f18740e1