text:
using Plots
using Statistics
using StatsBase
import Base.length
# histograms

"""
    digitize(x, bins)

    Return the indices of the bins to which each value in input array belongs.
"""
function digitize(x::Vector{<:Real}, bins::LinRange{<:Real})
    return searchsortedlast.(Ref(bins), x)
end
# Notice that in this case the bins are expected to be sorted
function digitize(x::Vector{<:Real}, bins::Vector{<:Real})
    return searchsortedlast.(Ref(bins), x)
end

"""
    hist1d(x::Vector{T}, nbins::Integer, xmin::T=typemin(T), xmax::T=typemax(T))
    hist1d(x::Vector{T}, xs::String, nbins::Integer,
                    xmin::T=typemin(T), xmax::T=typemax(T); datap = true, fwl=true)
    hist1d(h::Histogram, xs::String; datap = true, markersize=3, fwl=false)

return a 1d histogram and its corresponding graphics (plots)
"""
function hist1d(x::Vector{<:Real}, bins::Vector{<:Real}, norm=false)
    h = fit(Histogram, x, bins)
    if norm
        h = StatsBase.normalize(h, mode=:density)
    end
    return h
end

function hist1d(x::Vector{T}, nbins::Integer,
                xmin::T=typemin(T), xmax::T=typemax(T), norm=false) where T
    xx   = in_range(x, xmin, xmax)
    bins = collect(LinRange(xmin, xmax, nbins + 1))
    h    = hist1d(xx, bins, norm)
    return h
end

"""
    edges(h::Histogram)

edges of the histogram
"""
function edges(h::Histogram)
    collect(h.edges[1])
end

"""
    centers(h::Histogram)

return the centres of the histogram bins.
"""
centers(h::Histogram) = centers(edges(h))

"""
    centers

Calculate the bin centers from a vector of bin edges
"""
function centers(edges::Vector{<:Real})
    edges[1:end-1] + .-(@view(edges[2:end]), @view(edges[1:end-1])) / 2
end

"""
    hist_weights

return histogram a set of data and return the bin weights.
"""
function hist_weights(edges::Vector{<:Real})
  function get_weights(y::SubArray{<:Real})
    histo = fit(Histogram, y, edges)
    return histo.weights
  end
  return get_weights
end

function hist1d(x::Vector{T}, xs::String, nbins::Integer,
                xmin::T=typemin(T), xmax::T=typemax(T);
                norm=false, datap = true, markersize=3, fwl=false,
                label="", legend=false) where T

    return hist1d(hist1d(x, nbins, xmin, xmax, norm), xs,
                         datap=datap, markersize=markersize, fwl=fwl,
                         label=label, legend=legend)
end

function hist1d(h::Histogram, xs::String; datap=true, markersize=3, fwl=false,
                label="", legend=false)

    if datap
        yg = h.weights * 1.0
        xg = centers(h)
        p = scatter(xg,yg, yerr = sqrt.(yg), markersize=markersize, label=label, legend=legend)
        if fwl
            p = plot!(p, xg,yg, yerr = sqrt.(yg), fmt = :png,
                      linewidth=1, label=label, legend=legend)
        end
    else
        p = plot(h, xlabel=xl, fmt = :png, yl="frequency")
    end
    xlabel!(xs)
    ylabel!("frequency")

    return h, p
end

function hist1d(h1::Histogram, h2::Histogram, xs::String; markersize=2, norm=false)

    if norm
        h1 = StatsBase.normalize(h1, mode=:density)
        h2 = StatsBase.normalize(h2, mode=:density)
    end

    yg1 = h1.weights * 1.0
    xg1 = centers(h1)
    yg2 = h2.weights * 1.0
    xg2 = centers(h2)

    p1 = scatter(xg1,yg1, yerr = sqrt.(yg1), fmt = :png, markersize=markersize)
    p  = scatter!(p1,xg2,yg2, yerr = sqrt.(yg2), fmt = :png, markersize=markersize)
    xlabel!(xs)
    ylabel!("frequency")

    return p
end

"""
    hist2d(x,y, nbins, xl, yl)

return a 2d histogram and its corresponding graphics (plots)
"""
function hist2d(x::Vector{T},y::Vector{T}, nbins::Integer,
                xl::String, yl::String,
                xmin::T=typemin(T), xmax::T=typemax(T),
                ymin::T=typemin(T), ymax::T=typemax(T); title="") where T
    function xy(i)
        return diff(h.edges[i])/2 .+ h.edges[i][1:end-1]
    end

    df = DataFrame(x=y,y=x)
    df1 = select_by_column_value_interval(df, "y", xmin, xmax)
    df2 = select_by_column_value_interval(df1, "x", ymin, ymax)
    data = (df2.x, df2.y)
    h = fit(Histogram, data, nbins=nbins)
    ye = xy(1)
    xe = xy(2)
    hm = heatmap(xe, ye, h.weights)
    xlabel!(xl)
    ylabel!(yl)
    if title != ""
        title!(title)
    end
    return h,hm
end

"""
    fmean_std

returns the weighted mean and std of an array given a minimum
bin weight for consideration.
"""
function fmean_std(x::Vector{<:Real}, min_prop::Float64=0.1)
  mask_func = x -> x .> min_prop * maximum(x)
  function filt_mean(w::SubArray{<:Real})
      return [mean_and_std(x[mask_func(w)],
                           FrequencyWeights(w[mask_func(w)]),
                           corrected=true)]
  end
  return filt_mean
end

"""
    p1df(x, y, nbins; ybin_width, ymin, ymax, min_proportion)

return a profile DataFrame. This is a DF in which the variable y is
histogrammed as a function of the average of variable x in each bin.
The keyword arguments allow for a filtering in the y variable by
range and min_proportion of the maximum in the bin.
TODO: Protect against fine binning that results in zeros.
"""
function p1df(x::Vector{<:Real}, y::Vector{<:Real}, nbins::Integer;
              ybin_width::Real=0.1, ymin::Real=minimum(y), ymax::Real=maximum(y),
              min_proportion::Real=0.0)
    df           = DataFrame(:y => y)
    x_upper      = nextfloat(maximum(x))
    bin_edges    = LinRange(minimum(x), x_upper, nbins + 1)
    df[!, "bin"] = digitize(x, bin_edges)
    bin_centers  = centers(collect(bin_edges))
    bin_width    = .-(@view(bin_edges[2:end]), @view(bin_edges[1:end-1]))
    ## Bin in y so filtering can be done on bin content
    y_upper      = nextfloat(ymax)
    y_bins       = collect(ymin:ybin_width:y_upper)
    y_centers    = centers(y_bins)
    y_weights    = combine(groupby(df, :bin),
                           :y => hist_weights(y_bins) => :yweights,
                           ungroup = false)
    bin_stats    = fmean_std(y_centers, min_proportion)
    ndf = combine(y_weights,
                  :yweights => bin_stats => [:y_mean, :y_std])
    ndf = ndf[:, [:y_mean, :y_std]]
    ndf[!, :x_mean] = bin_centers
    ndf[!, :x_std]  = bin_width / 2

    p1 = plot(ndf.x_mean,ndf.y_mean, yerror=ndf.y_std, fmt = :png,
              shape = :circle, color = :black, legend=false)
    return ndf, p1
end

max_stars_repo_name: jjgomezcadenas/ATools

max_issues_repo_name: jjgomezcadenas/ATools

max_forks_repo_name: jjgomezcadenas/ATools


text:
for steps in 50:500:10000
    binary = mean([binaryTokenChallenge(steps, 50, 60, 10, [Benchmarks.registryMean, Benchmarks.activeAgents, Benchmarks.accuracyBoost, Benchmarks.weightedAccuracyBoost, Benchmarks.accuracyBalanceCorrelation]) for i in 1:20])
    pro = mean([proRataTokenChallenge(steps, 50, 60, 10, [Benchmarks.registryMean, Benchmarks.activeAgents, Benchmarks.accuracyBoost, Benchmarks.weightedAccuracyBoost, Benchmarks.accuracyBalanceCorrelation]) for i in 1:20])
    println("$steps, $binary, $pro")
end

proRataTokenChallenge(100, 20, 60, 10, [Benchmarks.registryMean, Benchmarks.activeAgents, Benchmarks.accuracyBoost, Benchmarks.weightedAccuracyBoost, Benchmarks.accuracyBalanceCorrelation])

ext: jl

max_stars_repo_name: mlodykaz/tcr-simulation

max_issues_repo_name: mlodykaz/tcr-simulation

max_forks_repo_name: mlodykaz/tcr-simulation


text:
@testset "Max adj visit" begin
    gx = SimpleGraph(8)

    # Test of Min-Cut and maximum adjacency visit
    # Original example by Stoer

    wedges = [
        (1, 2, 2.),
        (1, 5, 3.),
        (2, 3, 3.),
        (2, 5, 2.),
        (2, 6, 2.),
        (3, 4, 4.),
        (3, 7, 2.),
        (4, 7, 2.),
        (4, 8, 2.),
        (5, 6, 3.),
        (6, 7, 1.),
        (7, 8, 3.)]

    m = length(wedges)
    eweights = spzeros(nv(gx), nv(gx))

    for (s, d, w) in wedges
        add_edge!(gx, s, d)
        eweights[s, d] = w
        eweights[d, s] = w
    end
    for g in testgraphs(gx)
      @test nv(g) == 8
      @test ne(g) == m

      parity, bestcut = @inferred(mincut(g, eweights))

      @test length(parity) == 8
      @test parity == [2, 2, 1, 1, 2, 2, 1, 1]
      @test bestcut == 4.0

      parity, bestcut = @inferred(mincut(g))

      @test length(parity) == 8
      @test parity == [2, 1, 1, 1, 1, 1, 1, 1]
      @test bestcut == 2.0

      v = @inferred(maximum_adjacency_visit(g))
      @test v == Vector{Int64}([1, 2, 5, 6, 3, 7, 4, 8])
    end

    g1 = SimpleGraph(1)
    for g in testgraphs(g1)
        @test @inferred(maximum_adjacency_visit(g)) == collect(vertices(g))
        @test @inferred(mincut(g)) == ([1], zero(eltype(g)))
    end
end


text:
function f(a,b)
    2a+3b
end


text:
using QuickCheck2, Base.Test

type mytype
	parameter::Float64
end

QuickCheck2.custom_generator(a::Number, b::Number, type_a::Bool, type_b::Bool, mytype) = (size) -> mytype(a+(b-a)*rand(Float64))   ## Note `QuickCheck2` prefix!

@testset "new type generator test" begin  #
	@test_formany 1<x::mytype<10, x.parameter==1
	@test_formany 1<x::mytype<10, 10>x.parameter>1
end


text:
using Clang

using CImGui.LibCImGui.CImGui_jll

const CIMPLOT_H = joinpath(dirname(CImGui_jll.libcimgui_path), "..", "include", "cimplot.h") |> normpath

ctx = DefaultContext()

# parse headers
parse_headers!(ctx, [CIMPLOT_H], args=[map(x -> "-I" * x, find_std_headers())..., "-DCIMGUI_DEFINE_ENUMS_AND_STRUCTS"], includes=[LLVM_INCLUDE])

# settings
ctx.libname = "libcimgui"
ctx.options["is_function_strictly_typed"] = false
ctx.options["is_struct_mutable"] = false

# write output
api_file = joinpath(@__DIR__, "libcimplot_api.jl")
api_stream = open(api_file, "w")

for trans_unit in ctx.trans_units
    root_cursor = getcursor(trans_unit)
    push!(ctx.cursor_stack, root_cursor)
    header = spelling(root_cursor)
    @info "wrapping header: $header ..."
    # loop over all of the child cursors and wrap them, if appropriate.
    ctx.children = children(root_cursor)
    for (i, child) in enumerate(ctx.children)
        child_name = name(child)
        child_header = filename(child)
        ctx.children_index = i
        # choose which cursor to wrap
        startswith(child_name, "__") && continue  # skip compiler definitions
        child_name in keys(ctx.common_buffer) && continue  # already wrapped
        child_header != header && continue  # skip if cursor filename is not in the headers to be wrapped

        wrap!(ctx, child)
    end

    # Remove leading 'ImPlot_' from function names
    @info "Renaming functions"
    for (i, entry) in enumerate(ctx.api_buffer)
        if entry.head == :function
            name = string(entry.args[1].args[1])
            if startswith(name, "ImPlot_")
                ctx.api_buffer[i].args[1].args[1] = Symbol(chop(name; head=7, tail=0))
            end
        end
    end

    @info "writing $(api_file)"
    println(api_stream, "# Julia wrapper for header: $(basename(header))")
    println(api_stream, "# Automatically generated using Clang.jl\n")
    print_buffer(api_stream, ctx.api_buffer)
    empty!(ctx.api_buffer)  # clean up api_buffer for the next header
end
close(api_stream)

# write "common" definitions: types, typealiases, etc.
common_file = joinpath(@__DIR__, "libcimplot_common.jl")
open(common_file, "w") do f
    println(f, "# Automatically generated using Clang.jl\n")
    print_buffer(f, dump_to_buffer(ctx.common_buffer))
end


text:
using GFX
using Test

@testset "GFX.jl" begin
    # Write your tests here.
end

max_stars_repo_name: wmkouw/gfx-node

max_issues_repo_name: wmkouw/GFX

max_forks_repo_name: wmkouw/GFX


text:
###########################################################################################################################################################
# Declaring parameter and model descriptives
# Call functions to solve the model
#IMPORTANT: The model needs ONE steady state parameter for each variable
using DsgeSolves
# Set Seed
srand(1234)
include(joinpath(Pkg.dir("DsgeSolves"),"src","projection_method.jl"))

# Define the exogenous variables and the calibrated parameters
const Rstar = 0.01
const PFstar = 1.0
const PXstar = 1.0
const eta = 1.5
const omega = 0.25
const beta = 0.990099009901
const alpha1 = 0.15
const kappa = 0.1
const chis = 0.1
const phi0 = 0.9
const phi1 = 1.5
const Θ_s = [Rstar PFstar PXstar eta omega beta alpha1 kappa chis phi0 phi1]

# Reduced form parameters
const ρ = 0.9
const σ = 0.01
const Θ_rf = [ρ σ]

# Set out the initial steady-state values
const Blam_ss = 3.22696890071349
const C_ss = 0.45793462256871
const F_ss = 0.0
const G_ss = 0.0
const K_ss = 0.02729049017873
const L_ss = 0.74732008502373
const P_ss = 0.99999998539410
const Pk_ss = 1.0000
const R_ss = 0.0100
const S_ss = 0.99999999315490
const W_ss = 0.28812562001981
const X_ss = 0.02729049020372
const Y_ss = 0.48522511277243
const Z_ss = 1.0
const Θ_ss = [Blam_ss C_ss F_ss G_ss K_ss L_ss P_ss Pk_ss R_ss S_ss W_ss X_ss Y_ss Z_ss]

# Aproximating function
const State = 3 # number state variables
const Euler = 2  # number euler equations
const Neuron = 1 # number of neurons
const Neuronx1 = Neuron + 1
const Param = State*Euler*Neuron # number of parameters for the Neural Network function
const Start = 4
const Other = 4

# Create the shocks
const T = 1000 # length of simulated data
const K = 300  # number of simulations
const B = Int(ceil(T*0.25)) #burn-in
# const shock = initialize_shock(T,K,σ)

const Θ = Parameter(Θ_s,Θ_rf,Θ_ss)
const vDes = Descriptive(T, K, B, State, Euler, Neuron, Neuronx1, Param, Start,Other)
const d1 = Normal(0,σ)
const m = Model(Θ,vDes,d1)

gammaf = [2.8342 -1.6406 -0.0868 0.5264 -0.0000 3.0591] # starting values
model = optimize(minimize_euler_error!, vec(gammaf), BFGS())

# Evaluate the model for the bestσlected parameters
vParam_star = reshape(Optim.minimizer(model),size(gammaf))
m = evaluate_model!(vParam_star,m)

# Model assessment
irf = plot_irf!(m,0.1,50,25)
jg = judd_gaspar_stat(m) # I AM NOT SURE THIS IS CORRECT
DM = haan_marcet_statistics(m) # I AM NOT SURE THIS IS CORRECT

# @unpack endo, state_endo, euler_error, other, shock = m.Var
# (Blam,C,F,G,κ,L,P,Pk,R,S,W,X,Y,Z) = endo

# plot


text:
using EngEconomics, Roots, Plots, ColorSchemes
plotly()

# Given
constructionYears = 5 # Construction takes this many years
constructionStarts = 5 # Construction starts in this many years
constructionA = 20000 # Construction annual cost in $1000, so if this number is 20,000, the real price is 20,000,000
maintenanceStarts = 1 + constructionYears + constructionStarts # Number of years when maintenance for the project will begin. It is simply the sum of how long construction will take and when construction starts plus 1 year
maintenanceA = 2000 # Maintenance annual cost in $1000, so if this number is 2,000, then the real price is 2,000,000
inflationRate = 0.03 # Inflation rate
MARR = 0.07 # Project's MARR
projectLifeTime = 30 # I did this in my head (starts working in 2020, supposed to last till 2050)

# Find
# (a.) The present worth of the project
realMARR = (1 + MARR) / (1 + inflationRate) - 1
NPW = (-constructionA * seriesPresentAmountFactor(realMARR, constructionYears)) * presentWorthFactor(realMARR, constructionStarts) + (-maintenanceA * seriesPresentAmountFactor(realMARR, projectLifeTime)) * presentWorthFactor(realMARR, constructionYears + constructionStarts)

# (b.) Construct a sensitivity graph showing the effects of 5 percent and 10
#      percent increases and decreases in the construction costs, maintenance
#      costs and inflation rate. To which is the present worth most sensitive?
sensitivity = [0.9, 0.95, 1.0, 1.05, 1.1]

# Construction Costs
realMARR = (1 + MARR) / (1 + inflationRate) - 1
NPW_constructionCost = (-constructionA * sensitivity * seriesPresentAmountFactor(realMARR, constructionYears)) * presentWorthFactor(realMARR, constructionStarts) .+ (-maintenanceA * seriesPresentAmountFactor(realMARR, projectLifeTime)) * presentWorthFactor(realMARR, constructionYears + constructionStarts)
plot(sensitivity, NPW_constructionCost)

# Maintenance Costs
realMARR = (1 + MARR) / (1 + inflationRate) - 1
NPW_maintenanceCost = (-constructionA * seriesPresentAmountFactor(realMARR, constructionYears)) * presentWorthFactor(realMARR, constructionStarts) .+ (-maintenanceA * sensitivity * seriesPresentAmountFactor(realMARR, projectLifeTime)) * presentWorthFactor(realMARR, constructionYears + constructionStarts)
plot!(sensitivity, NPW_maintenanceCost)

# Inflation Rate
realMARR_inflationRate = (1 + MARR) ./ (1 .+ sensitivity * inflationRate) .- 1
NPW_inflationRate = (-constructionA * seriesPresentAmountFactor.(realMARR_inflationRate, constructionYears)) .* presentWorthFactor.(realMARR_inflationRate, constructionStarts) .+ (-maintenanceA * seriesPresentAmountFactor.(realMARR_inflationRate, projectLifeTime)) .* presentWorthFactor.(realMARR_inflationRate, constructionYears + constructionStarts)
plot!(sensitivity, NPW_inflationRate)

println("From the plot, it\'s clear that the NPW cost is most sensitive to the initial construction cost because the slope of that line is much bigger than that of the other lines")


text:
sbmlfile = joinpath(@__DIR__, "data", "Ec_core_flux1.xml")

if !isfile(sbmlfile)
    Downloads.download(
        "http://systemsbiology.ucsd.edu/sites/systemsbiology.ucsd.edu/files/Attachments/Images/InSilicoOrganisms/Ecoli/Ecoli_SBML/Ec_core_flux1.xml",
        sbmlfile,
    )
end

cksum = bytes2hex(sha256(open(sbmlfile)))
if cksum != "01a883b364fa60582101ca1e270515e7fcb3fb2f60084d92e5ee45f9f72bbe50"
    @warn "The downloaded E Coli core flux model seems to be different from the expected one. Tests will likely fail." cksum
end

@testset "SBML flux model loading" begin
    mdl = readSBML(sbmlfile)

    @test typeof(mdl) == Model

    @test_logs (:error, "SBML reported error: File unreadable.") @test_throws AssertionError readSBML(
        sbmlfile * ".does.not.really.exist",
    )

    @test mdl.units["mmol_per_gDW_per_hr"] == 2.7777e-7 * u"mol * g^-1 * s^-1"

    @test length(mdl.compartments) == 2

    mets, rxns, S = getS(mdl)

    @test typeof(S) <: AbstractMatrix{Float64}
    @test typeof(getS(mdl; zeros = spzeros)[3]) <: SparseMatrixCSC{Float64}
    @test typeof(getS(mdl; zeros = zeros)[3]) == Matrix{Float64}

    @test length(mets) == 77
    @test length(rxns) == 77
    @test size(S) == (length(mets), length(rxns))

    # totally arbitrary value tests
    @test isapprox(sum(S), 42.1479)
    @test mets[10:12] == ["M_akg_e", "M_fum_c", "M_pyr_c"]
    @test rxns[10:12] == ["R_H2Ot", "R_PGL", "R_EX_glc_e_"]

    lbs = getLBs(mdl)
    ubs = getUBs(mdl)
    ocs = getOCs(mdl)

    @test length(ocs) == length(mets)
    @test ocs[40] == 1.0
    deleteat!(ocs, 40)
    @test all(ocs .== 0.0)

    @test length(getLBs(mdl)) == length(rxns)
    @test length(getUBs(mdl)) == length(rxns)

    getunit((val, unit)) = unit
    @test all([broadcast(getunit, lbs) broadcast(getunit, ubs)] .== "mmol_per_gDW_per_hr")

    getval((val, unit)) = val
    lvals = broadcast(getval, lbs)
    uvals = broadcast(getval, ubs)
    @test isapprox(lvals[27], uvals[27])
    @test isapprox(lvals[27], 7.6)
    @test isapprox(lvals[12], -10)

    @test count(isapprox.(lvals, -999999)) == 40
    @test count(isapprox.(lvals, 0)) == 35
    @test count(isapprox.(uvals, 999999)) == 76
end


text:
module FunProject
using Combinatorics

"""
Each word in a sentence (string) is sorted alphabetically.

#### Examples
```julia-repl
julia> order_string("This is a sorted sentence; Who is 1 and two  ")
" This is a deorst ;ceeennst Who is 1 adn otw  "

```

"""
function order_string(s::AbstractString="Hello world")

	s0=replace(s,"+"=>" + ","-"=>" - ")
	x=split(s0," ") # Split into words delim=" "

	Y=""
	for wIdx = 1:length(x)
		# w=x[wIdx] # word
		y=string(sort([a for a in x[wIdx]])...)  # sort each word
		Y= Y * " " * y
	end
	Yout=replace(Y," " =>" ")

	return Yout
end

export order_string

# Write your package code here.
greet() = "Hello World"
export greet

"""
Adding two functions
```julia-repl
julia> add_two(3,2)

```
More bla
"""
function add_two(x,y)
    return x+y
end

export add_two

end


text:
check_real(A::AbstractArray) = nothing
function check_real(A::AbstractArray{Complex{T}}) where T
    epsT = eps(real(eltype(A)))
    rtol = 1000epsT
    atol = 100epsT

    if any(abs(imag(x)) > rtol * abs(x) && abs(imag(x)) > atol for x in A)
        relerror = imag(A) ./ abs.(A)
        relerror[map(x -> abs(imag(x)) < atol, A)] .= 0
        @warn "Large imaginary part" norm(relerror)
    end
end


text:
mutable struct SparseSDPSolution{T<:Number}
    obj::T
    dualobj::T
    primalmatrix::SparseSymmetricBlockMatrix{T}
    dualvector::Dict{Any,T}
    dualmatrix::SparseSymmetricBlockMatrix{T}
end

SparseSDPSolution(T) = SparseSDPSolution(zero(T), zero(T), Dict(), Dict(), Dict())

SparseSDPSolution() = SparseSDPSolution(Float64)

SparseSDPSolution{T<:Number}(obj::T, dualobj::T) = SparseSDPSolution{T}(obj,
                                                                        dualobj,
                                                                        SparseSymmetricBlockMatrix(T),
                                                                        Dict{Any,T}(),
                                                                        SparseSymmetricBlockMatrix(T))

setobj!{T<:Number}(sol::SparseSDPSolution{T}, v::T) = sol.obj = v

obj(sol::SparseSDPSolution) = sol.obj

setdualobj!{T<:Number}(sol::SparseSDPSolution{T}, v::T) = sol.dualobj = v

dualobj(sol::SparseSDPSolution) = sol.dualobj

primalmatrix(sol::SparseSDPSolution) = sol.primalmatrix

dualvector(sol::SparseSDPSolution) = sol.dualvector

dualmatrix(sol::SparseSDPSolution) = sol.dualmatrix

function invdict(d::Dict)
    r = Dict()
    for (k, v) in d
        r[v] = k
    end
    r
end

function readcsdpoutput!(io::IO, sol::SparseSDPSolution, cm::Dict, bm::Dict, ems::Dict)
    icm = invdict(cm)
    ibm = invdict(bm)
    iems = Dict()
    for (k, v) in ems
        iems[ibm[k]] = invdict(v)
    end

    l = readline(io)
    ls = split(l)
    for i = 1:length(ls)
        dualvector(sol)[icm[i]] = float(ls[i])
    end

    for l in readlines(io)
        ls = split(l)
        onetwo = parse(Int, ls[1])
        blockindex = ibm[parse(Int, ls[2])]
        rowindex = iems[blockindex][parse(Int, ls[3])]
        columnindex = iems[blockindex][parse(Int, ls[4])]
        value = parse(Float64, ls[5])
        if onetwo == 2
            block = get(primalmatrix(sol).blocks, blockindex, SparseSymmetricMatrix())
            block[rowindex, columnindex] = value
            primalmatrix(sol).blocks[blockindex] = block
        else
            block = get(dualmatrix(sol).blocks, blockindex, SparseSymmetricMatrix())
            block[rowindex, columnindex] = value
            dualmatrix(sol).blocks[blockindex] = block
        end
    end

    sol
end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module Knet
using AutoGrad
export grad
export KnetArray, gradcheck, gpu, relu, sigm, invx, logp, conv4, pool, mat
const libknet8 = Libdl.find_library(["libknet8"], [dirname(@__FILE__)])
const datapath = joinpath(dirname(@__FILE__),"..","data")

include("gpu.jl")               # gpu support
include("karray.jl")            # use KnetArrays
include("cuda1.jl")             # unary operators
include("cuda01.jl")            # scalar,array->array
include("cuda10.jl")            # array,scalar->array
include("cuda11.jl")            # array,array->array (elementwise, same shape)
include("cuda12.jl")            # array,array->array (broadcasting)
include("cuda20.jl")            # array->scalar (reductions)
include("cuda21.jl")            # array->vector (reductions)
include("cuda22.jl")            # array,array->array (linear algebra)
include("cuda44.jl")            # convolution and pooling
include("gradcheck.jl")         # gradient check

# See if we have a gpu at initialization:
function __init__()
    try
        r = gpu(true)
        info(r >= 0 ? "Knet using GPU $r" : "No GPU found, Knet using the CPU")
    catch e
        warn("$e: Knet using the CPU.")
        gpu(false)
    end
end

end # module


text:
@testset "Magnitude" begin
    @test magnitude(A₀) == 0.9
    @test magnitude2(A₀) == 0.81
    @test magnitude(1.2) == magnitude(-1.2) == 1.2

    if triplet(SIGNATURE) == (1, 3, 0)
        @test magnitude(NHB₁) == sqrt(0.8^2 - 0.3^2)
        @test magnitude2(NHB₁) == 0.3^2 - 0.8^2
    else
        @test magnitude(NHB₁) == sqrt(0.3^2 + 0.8^2)
        @test magnitude2(NHB₁) == 0.3^2 + 0.8^2
    end
end


text:
######################################################################## # Codes for the paper:
# An Anthropocene Map of Genetic Diversity
#
# Andreia Miraldo, Sen Li, Michael K. Borregaard, Alexander Floréz-Rodriguéz, Shyam Gopalakrishnan, Mirneza Risvanovic, Zhiheng Wang, Carsten Rahbek, Katharine A. Marske & David Nogués-Bravo
#
# Science, 2016
# Code in this file by Sen Li and Michael K. Borregaard
#########################################################################

# Load the necessary library and functions
include("Compare_Pairwise_Function.jl")
using DataFrames

"""
A function to calculate the summary statistic for all sites (e.g. grid cell or biome) for one species.

**Parameters**
* 'species':        A string with the name of the species
* 'species_seqs':   A matrix where the rows are aligned genetic sequences, and columns are loci. Basepairs must be coded as 1, 2, 3 or 4, or with a 0 signifying that the locus is absent from the alignment.
* 'sitenames':      A vector of strings with the names of the sites
"""
function calcspecies(species::String, species_seqs::Matrix{Int}, sitenames::Vector{String})
    res = DataFrame(species = String[], cell = String[], seq1 = Int[], seq2 = Int[], length_seq1 = Int[], length_seq2 = Int[], overlap = Float64[], commons = Int[], num_per_bp = Float64[])
    uniq_grid = unique(sitenames, 1)
    for iter_grid = 1:size(uniq_grid,1) # Go through each site in turn
        seqs_grid = species_seqs[findin(sitenames, uniq_grid[iter_grid, :]), :]

        tot_muts = compare_pairwise(seqs_grid)
        lns = size(tot_muts, 1)
        tmp = DataFrame(species = repeat([species], inner= lns), cell = repeat(vec(uniq_grid[iter_grid, :]), inner = lns)) #Expand species and cell names to the length of the resulting DataFrame
        tot_muts = hcat(tmp, tot_muts)

        res = vcat(res, tot_muts)
    end

    for sym in [:seq1, :seq2, :length_seq1, :length_seq2, :overlap, :commons, :num_per_bp]
        res[sym][res[sym] .< 0.] = NA   #Replace empty values with an NA identifier
    end
    res
end

"""
A function to create master data matrices that are used to compute genetic diversity, assess data quality and do sensitivity analyses.

**Parameters**
* 'foldername' : The name of of a folder containing the data files. There must be files of two types (file extension 'fasta' and file extension 'coords') with  the same filename, e.g. the species names (e.g. folder contents could be 'Bufo_bufo.fasta, Bufo_bufo.coords, Rana_arvalis.fasta, Rana_arvalis.coords', etc.). The .fasta files contain the sequences as an m x n integer matrix, where m is the number of sequences and n is the length of the alignments. The .coords files contain the geographic coordinates of the sequences, as an m x 2 floating point matrix with latitude in the first column and longitude in the second.  TODO: A folder with the anthrome and biome of each sequence
"""
function create_master_matrices(foldername::String, sitesfolder::String)
    species_list = [x[1:(end-6)] for x in filter(st->contains(st, ".fasta"), readdir(foldername))]      #identify unique file names ignoring the extension
    num_files = size(species_list,1)

    #equalarea = latbands = biomes = anthrome = gridcells = DataFrame(species = String[], cell = String[], seq1 = Int[], seq2 = Int[], length_seq1 = Int[], length_seq2 = Int[], overlap = Float64[], commons = Int[], num_per_bp = Float64[])   # Pre-initialize the DataFrame to ensure correct element types
    ret = DataFrame(species = String[], cell = String[], seq1 = Int[], seq2 = Int[], length_seq1 = Int[], length_seq2 = Int[], overlap = Float64[], commons = Int[], num_per_bp = Float64[])   # Pre-initialize the DataFrame to ensure correct element types

    for iter_file = 1:num_files
        # read in the sequence data for one species
        file_name_seq = joinpath(foldername, species_list[iter_file] * ".fasta")
        species_seqs = readdlm(file_name_seq, Int)
        length_seq = size(species_seqs,2)

        # read the coordinates for one species
        file_name_sites = joinpath(foldername, species_list[iter_file] * ".coords")
        species_sites = readdlm(file_name_sites, String)[:,1]
        #if eltype(species_sites) <: AbstractFloat
        #    species_sites = floor(Int, species_sites)
        #end

        ret = vcat(ret, calcspecies(species_list[iter_file], species_seqs, species_sites))
        println(iter_file)

    end
    writetable("pairwise_$(sitesfolder)_$(foldername).csv", ret)
end

function grid_coordinates(foldername = ".", xsize = 1, ysize = xsize, lowerleftoffset = (0,0))
    in("gridcells", readdir(foldername)) && error("There is already a folder named gridcells!")
    mkdir(joinpath(foldername, "gridcells"))
    species_list = [x[1:(end-7)] for x in filter(st->contains(st, ".coords"), readdir(foldername))]
    for species in species_list
        species_coords = readdlm(joinpath(foldername,"$(species).coords"))
        coords_grid = hcat(ceil(Int, (species_coords[:,1] + lowerleftoffset[1]) ./ xsize) .* xsize .- lowerleftoffset[1] .- 0.5xsize, ceil(Int, (species_coords[:,2] + lowerleftoffset[2]) ./ ysize) .* ysize .- lowerleftoffset[2] .- 0.5ysize)
        sitenames = mapslices(x -> "$(x[1])_$(x[2])", coords_grid, 2)
        writedlm(joinpath(foldername, "gridcells", species), sitenames)
    end
end

function latband_coordinates(foldername = ".", bandwidth = 1)
    in("latbands", readdir(foldername)) && error("There is already a folder named latbands!")
    mkdir(joinpath(foldername, "latbands"))
    species_list = [x[1:(end-7)] for x in filter(st->contains(st, ".coords"), readdir(foldername))]
    for species in species_list
        species_coords = readdlm(joinpath(foldername,"$(species).coords"))
        latband = bandwidth * floor(species_coords[:,1]/bandwidth)
        writedlm(joinpath(foldername, "latbands", species), latband)
    end
end

max_stars_repo_path: Miraldo et al. 2016/Create_Master_Matrices.jl

max_stars_repo_name: mkborregaard/codes-for-published-papers

max_issues_repo_name: DalavanCloud/codes-for-published-papers

max_forks_repo_name: DalavanCloud/codes-for-published-papers


text:
# Thanks to Jiahao Chen for this great example!

##
## NOTE: THIS NEEDS TO BE PORTED OVER TO THE NEW ARCHITECTURE
##

using AudioIO
import AudioIO.play

type note{S<:Real, T<:Real}
    pitch::S
    duration::T
    sustained::Bool
end

function play(A::note, samplingfreq::Real=44100, shape::Function=t->0.6sin(t)+0.2sin(2t)+.05*sin(8t))
    timesamples=0:1/samplingfreq:(A.duration*(A.sustained ? 0.98 : 0.9))
    v = Float64[shape(2π*A.pitch*t) for t in timesamples]
    if !A.sustained
        decay_length = int(length(timesamples) * 0.2)
        v[end-decay_length:end-1] = v[end-decay_length:end-1] .* linspace(1, 0, decay_length)
    end
    play(v)
    sleep(A.duration)
end

function parsevoice(melody::String; tempo=132, beatunit=4, lyrics=nothing)
    play([0]) #Force AudioIO to initialize
    lyrics_syllables = lyrics==nothing? nothing : split(lyrics)

    note_idx = 1
    oldduration = 4
    for line in split(melody, '\n')
        percent_idx = findfirst(line, '%') #Trim comment
        percent_idx == 0 || (line = line[1:percent_idx-1])
        for token in split(line)
            pitch, duration, dotted, sustained =parsetoken(token)
            duration==nothing && (duration = oldduration)
            oldduration = duration
            dotted && (duration *= 1.5)
            if lyrics_syllables!=nothing && 1<=note_idx<=length(lyrics_syllables) #Print the lyrics, omitting hyphens
                if lyrics_syllables[note_idx][end]=='-'
                    print(lyrics_syllables[note_idx][1:end-1])
                else
                    print(lyrics_syllables[note_idx], ' ')
                end
            end
            play(note(pitch, (beatunit/duration)*(60/tempo), sustained))
            note_idx += 1
        end
        println()
    end
end

function parsetoken(token::String, Atuning::Real=220)
    state = :findpitch
    pitch = 0.0
    sustain = dotted = false
    lengthbuf = Char[]
    for char in token
        if state == :findpitch
            scale_idx = findfirst('a':'g', char) + findfirst('A':'G', char)
            if scale_idx!=0
                const halfsteps = [12, 14, 3, 5, 7, 8, 10]
                pitch = Atuning*2^(halfsteps[scale_idx]/12)
                state = :findlength
            elseif char=='r'
                pitch, state = 0, :findlength
            else
                error("unknown pitch: $char")
            end
        elseif state == :findlength
            if     char == '#' ; pitch *= 2^(1/12) #sharp
            elseif char == 'b' ; pitch /= 2^(1/12) #flat
            elseif char == '\''; pitch *= 2        #higher octave
            elseif char == ',' ; pitch /= 2        #lower octave
            elseif char == '.' ; dotted = true     #dotted note
            elseif char == '~' ; sustain = true    #tied note
            else
                push!(lengthbuf, char)
                #Check for "is" and "es" suffixes for sharps and flats
                if length(lengthbuf) >= 2
                    if lengthbuf[end-1:end] == "is"
                        pitch *= 2^(1/12)
                        lengthbuf = lengthbuf[1:end-2]
                    elseif lengthbuf[end-1:end] == "es"
                        pitch /= 2^(1/12)
                        lengthbuf = lengthbuf[1:end-2]
                    end
                end
            end
        end
    end
    #finalize length
    lengthstr = convert(String, lengthbuf)
    duration = isempty(lengthstr) ? nothing : parseint(lengthstr)
    return (pitch, duration, sustain, dotted)
end

parsevoice("""
f# f# g a a g f# e d d e f# f#~ f#8 e e2
f#4 f# g a a g f# e d d e f# e~ e8 d d2
e4 e f# d e f#8~ g8 f#4 d e f#8~ g f#4 e d e a,
f#2 f#4 g a a g f# e d d e f# e~ e8 d8 d2""",
lyrics="""
Freu- de, schö- ner Göt- ter- fun- ken, Toch- ter aus E- li- - si- um!
Wir be- tre- ten feu- er- trun- ken, Himm- li- sche, dein Hei- - lig- thum!
Dei- ne Zau- ber bin den - wie- der, was die - Mo- de streng ge- theilt,
al- le mensch- en wer- den Brü- der wo dein sanf- ter Flü- - gel weilt.
""")

# And now with harmony!

soprano = @async parsevoice("""
f'#. f'#. g'. a'. a'. g'. f'#. e'~ e'8 d.'4 d.' e.' f#'. f#'.~ f#' e'8 e'4~ e'2
""", lyrics="Freu- de, schö- ner Göt- ter- fun- ken, Toch- ter aus E- li- - si- um!"
)
alto = @async parsevoice("""
a. a. a. a.  a.  a. a. a~ g8 f#.4 a.  a.  a. a.~ a a8 a4~ a2
""")
tenor = @async parsevoice("""
d. d. e. f#. f#. e. d. d~ e8 f#.4 f#. a,. d. d.~ d c#8 c#4 c#2
""")
bass = @async parsevoice("""
d. d. d. d. a,. a,. a,. b,~ c8 d. a., a., a., a., a, a8, a,4 a,2
""")
wait(soprano)
wait(alto)
wait(tenor)
wait(bass)

soprano = @async parsevoice("""
f'#.4 f'#. g'. a'. a'. g'. f'#. e'. d'. d'. e'. f'#. e'.~ e' d'8 d'4~ d'2
""", lyrics="Wir be- tre- ten feu- er- trun- ken, Himm- li- sche, dein Hei- - lig- thum!")
alto = @async parsevoice("""
a.4 a. b. c'. c'. b. a. g. f#. f#. g. f#. g.~ g4 f#8 f#~ f#2
""")
tenor = @async parsevoice("""
d.4 d. d. d. d. d. d. d. d. d. c#. d. c#.~ c# d8 d d2
""")
bass = @async parsevoice("""
d.4 d. d. d. a,. a,. a,. a., a., a., a., a., a.,~ a, a,8 d, d,2
""")

wait(soprano)
wait(alto)
wait(tenor)
wait(bass)


text:
# //  This is a gear train design problem taken from the GAMS test library
# //
# //  A compound gear train is to be designed to achieve a specific
# //  gear ratio between the driver and driven shafts. The objective
# //  of the gear train design is to find the number of teeth of the
# //  four gears and to obtain a required gear ratio of 1/6.931.
# //
# //  The problem originated from:
# //  Deb, K, and Goyal, M, Optimizing Engineering Designs Using a
# //  Combined Genetic Search. In Back, T, Ed, Proceedings of the
# //  Seventh International Conference on Genetic Algorithms. 1997,
# //  pp. 521-528.

function gear(oct::Bool = false)
    m = JuMP.Model()
    @variable(m, 12 <= i[1:4] <= 60, Int)
    for j = 1:4
        JuMP.set_start_value(i[j], 24)
    end
    @constraint(m, e2, - i[3] + i[4] >= 0)
    @constraint(m, e3, i[1] - i[2] >= 0)
    if !oct
        @NLobjective(m, Min, (6.931 - i[1]*i[2]/(i[3]*i[4]))^2 + 1)
        return m
    else
        @variable(m, obj)
        @objective(m, Min, obj)
        gm = GlobalModel(model = m, name = "gear")
        add_nonlinear_constraint(gm, :((i) -> (6.931 - i[1]*i[2]/(i[3]*i[4]))^2 + 1), dependent_var = obj)
        return gm
    end
end


text: revstring (str) = join(reverse(split(str, " ")), " ")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Example shows how to egnerate a whole batch of responses with the same label
# (volume fraction, radius and shape) but different realisations (actual
# positions of the particles). We will then extract stochastic information from
# them
using MultipleScattering

function moments_example()
    volfrac = 0.01
    radius = 1.0
    num_particles = 10
    k_arr = collect(LinRange(0.01,1.0,100))

    # region to place the particles
    shape = Rectangle(volfrac, radius, num_particles)
    # Holder for our simulations
    simulations = Vector{FrequencySimulation{Float64}}(undef,10)
    for i=1:10
        simulations[i] = FrequencySimulation(volfrac,radius,k_arr;seed=[0x7f5def91, 0x82255da3, 0xc5e461c7, UInt32(i)])
    end

    moments = StatisticalMoments(simulations)
    return moments
end


text:
# it's supposed that basis is a vector of its elements
# Each element of the basis can be a multidimensional array ()
# Internally it can be implemented differently: as one multidimensional array or as vector of arrays.
# Externally it should make no difference via the following organization.
using RecursiveArrayTools

import Base: length

"""
Abstract type representing any set of functions (`elements`) defined on some subset `ap` of a Cartesian domain.

    `elements(b::Basis)` gives vector of the basis functions.
    `aperture(b::Basis)` gives the  integer mask of non-zero elements of the array.
    `length(b::Basis)` gives the total number of elements in the basis.
"""
abstract type Basis end

abstract type OrthogonalBasis <: Basis end

abstract type OrthonormalBasis <: OrthogonalBasis end

elements(b::Basis) =  b.elements
elements(b::Basis, ind) =  b.elements[ind]

aperture(b::Basis)  = b.ap

# Todo: correct below, returns array instead of VectorOfArray
aperturedelements(b::Basis) =  (b.ap) .* b.elements
aperturedelements(b::Basis, ind) =  (b.ap) .* b.elements[ind]

norms(b::Basis) =  b.norms[:]

norms(b::Basis, ind::Vector) =  b.norms[ind]

length(b::Basis) = length(elements(b))

@doc raw"""
    compose(b, ind, coef) gives linear combination  math`\sum_{i\in ind \lambda_i f_i`

    compose(b, coef) requires `coef` be the complete vector of the coefficients.
"""
compose(b::Basis, ind::Vector, coef::Vector) = length(ind) != length(coef) ? error("Coefficients and indexes are of different length") : comb(coef, elements(b,ind))
compose(b::Basis,  coef::Vector) =  length(elements(b)) != length(coef) ? error("Coefficient vector does not match length of basis") : comb(coef, elements(b))

@doc raw"""

    decompose(a::Array, b::Basis)

Calculate coefficients of `a` in basis `b`.
"""
function decompose(a::Array, b::Basis)
    error("There are no decomposition rules of basis $b")
end

function decompose(a::Array, b::OrthonormalBasis)
    [inner(a,  f, aperture(b) ) for f in elements(b) ]
end

function decompose(a::Array, b::OrthogonalBasis)
    [inner(a,  f, aperture(b) ) / n^2 for (f,n) in zip(elements(b),norms(b)) ]
end

# Function below is for basis implemented as VectorOfArray
# Do we need the same for multidimensional array?
function comb(coef::Vector, a::VectorOfArray)
    sum = similar(a[1]) .* 0
    for i = 1: length(coef)
        sum += coef[i] * a[i]
    end
    return sum
end

# TODO rewrite so it works as tensor inner product if applied to two multidimensional arrays (or use Tensors.jl/ TensorOperations.jl?)
# a proper version of comb
# It can be realised as inner product in Mathematica
# Calculate inner product (tensor contraction) using the last index in `a` and first index in `b`
# Can be also related to scalar product in Banach space. Then inner product `(c, f)` of function `f` and
# constant `c` is defined as inner of `(f) ⋅ (c) = c f` .
"""
    inner(a, b)

Calculate inner product of functions or array of functions.
"""
inner(a::Array{T,N}, b::Array{T,N}) where {T <: Number, N} = integrate(a .* b)
inner(a::Array{T,N}, b::Array{T,N}, domain) where {T <: Number, N} = integrate(a .* b , domain)

inner(a::Number, b::Number) = dot(a,b)
# inner(c::Number, a::Array) = dot(c, sum(a))
# inner(a::Array, c::Number) = dot(sum(a), c)
inner(c::Number, a::Array) = conj(c) * a
inner(a::Array, c::Number) = conj(a) * c

inner(a::Union{Vector, VectorOfArray}, b::Union{Vector, VectorOfArray}) = sum(inner(va,vb) for (va,vb) in zip(a, b))

innermatrix(b::Basis) = [inner(elements(b,i), b.ap .* elements(b,j)) for i in eachindex(elements(b)), j in eachindex(elements(b))]

"""
    integrate(a)

calculate numerical integral of a
"""
# rectangle rule
integrate(a::Array{T,N}) where {T <: Number, N} = sum(a[:])
integrate(a::Array{T,N}, domain::Array{T,N}) where {T <: Number, N} = sum(a[:] .* domain[:])

# TODO simpson's rule on unit disk


text:
abstract type InitializationStrategyCombined <: InitializationStrategy end

struct SimpleCombinedStrategy <: InitializationStrategyCombined
    gamma_strategy::InitializationStrategyGamma
    C_strategy::InitializationStrategyC
end

FixedParameterInitialization(kernel, C) = SimpleCombinedStrategy(FixedGammaStrategy(kernel), FixedCStrategy(C))

function get_parameters(model, strategy::SimpleCombinedStrategy)
    C = calculate_C(model, strategy.C_strategy)
    gamma = calculate_gamma(model, strategy.gamma_strategy)
    return (C, MLKernels.GaussianKernel(gamma))
end

struct GammaFirstCombinedStrategy <: InitializationStrategyCombined
    gamma_strategy::InitializationStrategyGamma
    C_strategy::InitializationStrategyC
end

function get_parameters(model, strategy::GammaFirstCombinedStrategy)
    gamma = calculate_gamma(model, strategy.gamma_strategy)
    C = calculate_C(model, strategy.C_strategy, MLKernels.GaussianKernel(gamma))
    return (C, MLKernels.GaussianKernel(gamma))
end

struct SimpleSubspaceStrategy <: InitializationStrategyCombined
    gamma_strategy::InitializationStrategyGamma
    C_strategy::InitializationStrategyC
    gamma_scope::Scope

    function SimpleSubspaceStrategy(gamma_strategy, C_strategy; gamma_scope)
        new(gamma_strategy, C_strategy, gamma_scope)
    end
end

function get_parameters(model, strategy::SimpleSubspaceStrategy)
    C = calculate_C(model, strategy.C_strategy)

    if isa(strategy.gamma_scope, Val{:Subspace})
        gamma = @eachsubspace calculate_gamma(model, strategy.gamma_strategy)
    else
        tmp = calculate_gamma(model, strategy.gamma_strategy)
        gamma = fill(tmp, length(model.subspaces))
    end
    return (C, MLKernels.GaussianKernel.(gamma))
end

"""
    Determines C for WangGammaStrategy by using a provided C init strategy
"""
struct WangCombinedInitializationStrategy <: InitializationStrategyCombined
    solver::JuMP.OptimizerFactory
    gamma_search_range
    C_strategy::InitializationStrategyC
end

WangCombinedInitializationStrategy(solver, C_strategy) = WangCombinedInitializationStrategy(solver, 10.0.^range(-2, stop=2, length=50), C_strategy)

function get_parameters(model, strategy::WangCombinedInitializationStrategy)
    C = calculate_C(model, strategy.C_strategy)
    wang_strategy = WangGammaStrategy(strategy.solver, strategy.gamma_search_range, C)
    gamma = calculate_gamma(model, wang_strategy)
    return (C, MLKernels.GaussianKernel(gamma))
end


text:
# Runtime log file tracking.

"""
    readlog(file)

Read the run log file, check the iteration status and return the timestamps (exluding the
last) as well as the model running speed in physical seconds per simulated seconds.
"""
function readlog(file)
   isfile(file) || throw(ArgumentError("Cannot open \'$file\': not a file"))

   fid = open(file)
   lines = readlines(fid)
   close(fid)

   nline = findall(startswith("---"), lines)
   steps = lines[nline]
   stepsizes = [parse(Float32, split(step)[10]) for step in steps]
   t = lines[nline .- 1]

   local timestamps
   try
      timestamps = DateTime.(t, dateformat"e u d H:M:S y ")
   catch
      timestamps = DateTime.(t, dateformat"e u  d H:M:S y ")
   end

   timediff = diff(timestamps) / 1000 # duration for each step in [s]
   speed = [timediff[i].value / stepsizes[i] for i in eachindex(timediff)]
   return @views timestamps[1:end-1], speed
end


text:
t = [datetime(2011,09,13,13,53,36)
     datetime(2011,09,13,13,53,44)
     datetime(2011,09,13,13,53,49)
     datetime(2011,09,13,13,53,54)
     datetime(2011,09,13,14,32,53)
     datetime(2011,09,13,14,35,04)  ]

v1 = [ 5.8 5.83 5.9 6.0 5.95 5.88]
v2 = [ 1.0 3.0 1.0 20.0 12.4521 7.458]

ticks = FinancialTimeSeries(t, [v1, v2]', ["price","amount"], Stock(Ticker("ASDS")))

dt = DateTimeRange(datetime(2011,9,13,13,55,00), minutes(5), 4)


text:
"""
All the terminals are under my command. `Comonicon` is a CLI (Command Line Interface) generator
that features light-weight dependency (optional to have zero dependency),
fast start-up time and easy to use. See the [website](https://rogerluo.me/Comonicon.jl/)
for more info.
"""
module Comonicon

using Markdown
using Pkg
using Libdl
using ExprTools
using PackageCompiler

include("types.jl")
include("codegen/codegen.jl")
include("parse/parse.jl")

export @cast, @main

using .Types
using .Parse
using .CodeGen

using .Parse: main

include("tools/tools.jl")
include("tools/path.jl")
include("tools/build.jl")

using .BuildTools

# include("build.jl")
# include("precompile.jl")
# _precompile_()

end # module


text:
const JULIACOMPUTING_REPOS = MirrorUpdater.Types.SrcDestPair[
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JuliaComputing/Deprecations.jl",
        destination_repo_name = "JuliaComputing-Deprecations.jl",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JuliaComputing/FemtoCleaner.jl",
        destination_repo_name = "JuliaComputing-FemtoCleaner.jl",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JuliaComputing/NewPkgEval.jl",
        destination_repo_name =
            "NewPkgEval.jl-9f2e2246-6dce-11e8-3d98-4b291446da6e",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JuliaComputing/Registrator.jl",
        destination_repo_name =
            "Registrator.jl-4418983a-e44d-11e8-3aec-9789530b3b3e",
        ),
    ]


text:
module Helium_HF_FEM_SCF
    include("helium_hf_fem_eigen.jl")
    include("helium_vh_fem.jl")
    using Printf
    using .Helium_HF_FEM_Eigen
    using .Helium_Vh_FEM

    const DR = 0.01
    const MAXITER = 1000
    const CUDA_THRESHOLD = 5.0E-4
    const THRESHOLD = 1.0E-12

    function scfloop()
        # Schrödinger方程式を解くルーチンの初期処理
        hfem_param, hfem_val = Helium_HF_FEM_Eigen.construct()

        # CUDAで計算する際のフラグ
        cuda_flag = true

        # 有限要素法のデータのみ生成
        Helium_HF_FEM_Eigen.make_wavefunction(cuda_flag, 0, hfem_param, hfem_val, nothing)

        # Poisson方程式を解くルーチンの初期処理
        vh_param, vh_val = Helium_Vh_FEM.construct(hfem_param)

        # 仮の電子密度でPoisson方程式を解く
        Helium_Vh_FEM.solvepoisson!(0, hfem_param, hfem_val, vh_val)

        # 新しく計算されたエネルギー
        enew = 0.0

        for iter in 1:MAXITER
            eigenenergy = Helium_HF_FEM_Eigen.make_wavefunction(cuda_flag, iter, hfem_param, hfem_val, vh_val)

            # 前回のSCF計算のエネルギーを保管
            eold = enew

            # 今回のSCF計算のエネルギーを計算する
            enew = Helium_HF_FEM_Eigen.get_totalenergy(eigenenergy, hfem_param, hfem_val, vh_val)

            @printf "Iteration # %2d: HF eigenvalue = %.14f, energy = %.14f\n" iter eigenenergy enew

            # 今回のSCF計算のエネルギーと前回のSCF計算のエネルギーの差の絶対値
            ediff = abs(enew - eold)

            # SCF計算が収束したかどうか
            if ediff <= THRESHOLD
                # 波動関数を規格化
                hfem_val.phi ./= sqrt(4.0 * pi)

                # 収束したのでhfem_param, hfem_val, エネルギーを返す
                return hfem_param, hfem_val, enew
            elseif ediff <= CUDA_THRESHOLD
                cuda_flag = false
            end

            # Poisson方程式を解く
            Helium_Vh_FEM.solvepoisson!(scfloop, hfem_param, hfem_val, vh_val)
        end

        return nothing
    end

    save_result(hfem_param, hfem_val, filename) = let
        open(filename, "w" ) do fp
            len = length(hfem_val.node_r_glo)
            imax = floor(Int, (hfem_val.node_r_glo[len] - hfem_val.node_r_glo[1]) / DR)
            for i = 0:imax
                r = hfem_val.node_r_glo[1] + DR * float(i)
                println(fp, @sprintf "%.14f, %.14f" r Helium_Vh_FEM.phi(hfem_param, hfem_val, r))
            end
        end
    end
end


text:
# Should probably be a AbstractArray{T}??
abstract type PushStack{T} end

struct Stack{T} <: PushStack{T}
    st::Vector{T}
    Stack{T}() where {T} = new{T}(T[])
end

import Base: length, push!, pop!, last, empty!, getindex, lastindex

@inline length(s::Stack) = length(s.st)
@inline push!(s::Stack{T}, o::T) where {T} = push!(s.st, o)
@inline push!(s::Stack, o) = error("Objects of type $(typeof(o)) cannot be put on a $(typeof(s))")
@inline pop!(s::Stack) = pop!(s.st)
@inline function pop2!(s::Stack)
    a = pop!(s.st)
    b = pop!(s.st)
    return (a, b) # Note that element at the back is returned first
end
@inline last(s::Stack) = last(s.st)
@inline empty!(s::Stack) = empty!(s.st)
@inline getindex(s::Stack, args...) = getindex(s.st, args...)
@inline lastindex(s::Stack) = lastindex(s.st)


text:
function xtype(entry::BibInternal.AbstractEntry)
    str = "other"
    if typeof(entry) == Article
        str = "journal"
    elseif typeof(entry) == InProceedings
        str = "conference"
    end
    return str
end

function xtitle(entry::T) where T <: BibInternal.AbstractEntry
    return :title ∈ fieldnames(typeof(entry)) ? entry.title : get(entry.other, "title", "")
end

function xnames(entry::BibInternal.AbstractEntry, editors::Bool=false)
    if !editors && typeof(entry) ∈ [Proceedings]
        return xnames(entry, true)
    end
    names = editors ? entry.editor : entry.author
    str = ""
    start = true
    stop = 1
    for s in split(names, r"[\n\r ]and[\n\r ]")
        aux = split(s, r"[\n\r ]*,[\n\r ]*")
        str *= start ? "" : ", "
        start = false
        str *= aux[end]
        for t in aux[1:end - 1]
            str *= " " * t
        end
    end
    return str
end

# handle all entries
function xin(entry::BibInternal.AbstractEntry)
    fn = fieldnames(typeof(entry))
    str = ""
    if typeof(entry) == Article
        str *= entry.journal * ", " * entry.volume
        if entry.number != ""
            str *= "($(entry.number))"
        end
        if entry.pages != ""
            str *= ", " * entry.pages
        end
    elseif typeof(entry) == InProceedings
        str *= "Proceedings of the " * entry.booktitle
    elseif typeof(entry) == InBook
        str *= get(entry.other, "booktitle", "")
    end
    str *= str == "" ? "" : "."
    return str
end

function xyear(entry::BibInternal.AbstractEntry)
    return :year ∈ fieldnames(typeof(entry)) ? entry.year : get(entry.other, "year", "")
end

function xlink(entry::BibInternal.AbstractEntry)
    fn = fieldnames(typeof(entry))
    str = ""
    if :doi ∈ fn
        str = "https://doi.org/" * entry.doi
    elseif "doi" ∈ keys(entry.other)
        str = "https://doi.org/" * entry.other["doi"]
    elseif :url ∈ fn
        str = entry.url
    elseif "url" ∈ keys(entry.other)
        str = entry.other["url"]
    end
    return str
end

function xfile(entry::BibInternal.AbstractEntry)
    return "files/$(entry.id).pdf"
end

function xcite(entry::BibInternal.AbstractEntry)
    string(entry)
end

function xlabels(entry::BibInternal.AbstractEntry)
    str = get(entry.other, "labels", "")
    return split(str, r"[\n\r ]*,[\n\r ]*")
end

struct Publication
    id::AbstractString
    type::AbstractString
    title::AbstractString
    names::AbstractString
    in::AbstractString
    year::AbstractString
    link::AbstractString
    file::AbstractString
    cite::AbstractString
    labels::Vector{AbstractString}
end

function Publication(entry::T) where T <: BibInternal.AbstractEntry
    id = entry.id
    type = xtype(entry)
    title = xtitle(entry)
    names = xnames(entry)
    in_ = xin(entry)
    year = xyear(entry)
    link = xlink(entry)
    file = xfile(entry)
    cite = export_bibtex(entry)
    labels = xlabels(entry)
    return Publication(id, type, title, names, in_, year, link, file, cite, labels)
end

function export_web(bibliography::DataStructures.OrderedSet{BibInternal.AbstractEntry})
    entries = Vector{Publication}()
    for entry in bibliography
        p = Publication(entry)
        push!(entries, p)
    end
    return entries
end


text:
iter_or_array(x) = repeated(x)
iter_or_array(x::Repeated) = x
iter_or_array(x::AbstractArray) = x
# We treat staticarrays as scalar
iter_or_array(x::StaticArray) = repeated(x)

function metrics_bb(char::Char, font::FTFont, pixel_size)
    extent = get_extent(font, char) .* Vec2f0(pixel_size)
    mini = bearing(extent)
    return Rect2D(mini, Vec2f0(extent.scale)), extent
end

function boundingbox(char::Char, font::FTFont, pixel_size)
    bb, extent = metrics_bb(char, font, pixel_size)
    return bb
end

function glyph_ink_size(char::Char, font::FTFont, pixel_size)
    bb, extent = metrics_bb(char, font, pixel_size)
    return widths(bb)
end

"""
    iterate_extents(f, line::AbstractString, fonts, scales)
Iterates over the extends of the characters (glyphs) in line!
Newlines will be drawn like any other character.
`fonts` can be a vector of fonts, or a single font.
`scales` can be a single float or a Vec2, or a vector of any of those.

`f` will get called with `(char::Char, glyph_box::Rec2D, glyph_advance::Point2f0)`.

`char` is the currently iterated char.

`glyph_box` is the boundingbox of the glyph.
widths(box) will be the size of the bitmap, while minimum(box) is where one starts drawing the glyph.
For the minimum at y position, 0 is the where e.g. `m` starts, so `g` will start in the negative, while `^` will start positive.

`glyph_advance` The amount one advances after glyph, before drawing next glyph.
"""
function iterate_extents(f, line::AbstractString, fonts, scales)
    iterator = zip(line, iter_or_array(scales), iter_or_array(fonts))
    lastpos = 0.0
    for (char, scale, font) in iterator
        glyph_box, extent = metrics_bb(char, font, scale)
        mini = minimum(glyph_box) .+ Vec2f0(lastpos, 0.0)
        glyph_box = Rect2D(mini, widths(glyph_box))
        glyph_advance = Point2f0(extent.advance)
        lastpos += glyph_advance[1]
        f(char, glyph_box, glyph_advance)
    end
end

function glyph_rects(line::AbstractString, fonts, scales)
    rects = Rect2D[]
    iterate_extents(line, fonts, scales) do char, box, advance
        push!(rects, box)
    end
    return rects
end

function boundingbox(line::AbstractString, fonts, scales)
    return reduce(union, glyph_rects(line, fonts, scales))
end

function inkboundingbox(ext::FontExtent)
    l = leftinkbound(ext)
    r = rightinkbound(ext)
    b = bottominkbound(ext)
    t = topinkbound(ext)
    return FRect2D((l, b), (r - l, t - b))
end

function height_insensitive_boundingbox(ext::FontExtent, font::FTFont)
    l = leftinkbound(ext)
    r = rightinkbound(ext)
    b = descender(font)
    t = ascender(font)
    return FRect2D((l, b), (r - l, t - b))
end


text:
function gradient(obj::ParametricUnivariateFunction, arg0::jdouble, arg1::Vector{jdouble})
    return jcall(obj, "gradient", Vector{jdouble}, (jdouble, Vector{jdouble}), arg0, arg1)
end

function value(obj::ParametricUnivariateFunction, arg0::jdouble, arg1::Vector{jdouble})
    return jcall(obj, "value", jdouble, (jdouble, Vector{jdouble}), arg0, arg1)
end


text:
using Test
using Deductive

include("common_algorithms.jl")
include("expression.jl")
include("manipulation.jl")
include("truthtable.jl")
include("sets.jl")
include("indexes.jl")
include("proof.jl")


text:
module TestConstructors

using MLJIteration
using MLJBase
using ..DummyModel
using Test

struct Foo <: MLJBase.Unsupervised end
struct Bar <: MLJBase.Deterministic end

@testset "constructors" begin
    model = DummyIterativeModel()
    @test_throws MLJIteration.ERR_NO_MODEL IteratedModel()
    @test_throws MLJIteration.ERR_NOT_SUPERVISED IteratedModel(model=Foo())
    @test_throws MLJIteration.ERR_NOT_SUPERVISED IteratedModel(model=Int)
    @test_throws MLJIteration.ERR_NEED_MEASURE IteratedModel(model=Bar())
    @test_throws MLJIteration.ERR_NEED_PARAMETER IteratedModel(model=Bar(),
                                                                measure=rms)
    iterated_model = @test_logs((:info, r"No measure"),
               IteratedModel(model=model))
    @test iterated_model.measure == RootMeanSquaredError()
    @test_logs IteratedModel(model=model, measure=mae)

    @test_logs IteratedModel(model=model, resampling=nothing)

    @test_logs((:info, r"`resampling` must be"),
               IteratedModel(model=model,
                             resampling=CV(),
                             measure=rms))
    @test_logs((:info, r"A `resampling` vector may contain"),
               IteratedModel(model=model,
                             resampling=[([1, 2], [3, 4]),
                                         ([3, 4], [1, 2])],
                             measure=rms))
    @test_logs IteratedModel(model=model,
                             resampling=[([1, 2], [3, 4]),],
                             measure=rms)

    @test_throws(MLJIteration.ERR_MISSING_TRAINING_CONTROL,
                 IteratedModel(model=model,
                               resampling=nothing,
                               controls=[Patience(), InvalidValue()]))
end

end

true

max_stars_repo_name: ablaom/MLJIteration

max_issues_repo_name: ablaom/MLJIteration

max_forks_repo_name: ablaom/MLJIteration


text:
mutable struct IntervalTreeNode{R,S,T<:Real}
        min           ::R
        max           ::S
        stepSize      ::T
        subtreeMax    ::Float64
        rightNode     ::Union{Nothing,IntervalTreeNode}
        leftNode      ::Union{Nothing,IntervalTreeNode}
        parentNode    ::Union{Nothing, IntervalTreeNode}
end

mutable struct Interval{R,S,T<:Real}
        min           ::R
        max           ::S
        stepSize      ::T
end

Interval(min,max) = Interval(min,max,1)
IntervalTreeNode(min,max) = IntervalTreeNode(min, max, 1, convert(Float64, max), nothing, nothing, nothing)
IntervalTreeNode(min,max,stepSize) = IntervalTreeNode(min, max, stepSize, convert(Float64, max), nothing, nothing, nothing)

Base.show(io::IO,node::IntervalTreeNode) = Base.print(io, "{ [",node.min,",",node.max,"]"," stepSize:",node.stepSize," subtreeMax:",node.subtreeMax, " }")

Base.:(==)(x::IntervalTreeNode,y::IntervalTreeNode) = x.min==y.min && x.max==y.max && x.stepSize==y.stepSize ? true : false

Base.isless(x::IntervalTreeNode,y::IntervalTreeNode) = (x.min<y.min ? true : (x.min==y.min && x.max<y.max) ? true : false)


text:
"""
	loadeop(filein)
(Download &) read EOP data (C04)

**Input**
* filein: either URL to EOP C04 or file name with downloaded EOP C04 parameters

**Output**
* DataFrame with DateTime, and all columns of the input file

**Example**
```
# Download data automatically
eop = loadeop("http://hpiers.obspm.fr/iers/eop/eopc04/eopc04_IAU2000.62-now");
# Or read data from already downloaded file
eop = loadeop("/test/input/eop_parameters.c04");
```
"""
function loadeop(filein::String)
	# Download file if necessary
	downto = downfile(filein);
	if !isempty(downto)
		channels = eopchannels(downto);
		dataall = readdlm(downto,skipstart=14)
		eopfile = DataFrame(datetime=FileTools.mat2time(dataall[:,1:3]));
		for (i,v) in enumerate(channels)
			eopfile[v] = dataall[:,i+3];
		end
		# delete downloaded file
		delfile(filein);
		return eopfile;
	else
		return DataFrame(datetime=0);
	end
end

"""
	downfile(filein::String)
Auxiliary function for downloading of files on internet. The downloaded files
will be stored in a temporary file
"""
function downfile(filein::String)
	downto = "tempdownfile.down";
	if occursin("http:",filein) || occursin("ftp:",filein)
		try
			download(filein,downto)
		catch
			downto = "";
		end
		return downto;
	else
		return filein;
	end
end

"""
	delfile(filein)
Delete file IF downloaded
"""
function delfile(filein::String)
	if occursin("http:",filein) || occursin("ftp:",filein)
		if isfile("tempdownfile.down")
			rm("tempdownfile.down")
		end
	end
end

"""
	eopchannels(filein)
Get EOP header=channels
"""
function eopchannels(filein::String)
	channels = Array{Symbol,1}();
	open(filein,"r") do fid
		row = "";
		for i in 1:11
			row = readline(fid);
		end
		temp = split(row,"  ");
		for i in temp
			if (i != "") && (i != "Date")
				push!(channels,Symbol(replace(i," "=>"")));
			end
		end
	end
	return channels;
end


text:
using StatsBase: countmap

function compute(path::String)::Int
    is_flush(suits::Vector{Char})::Bool = length(Set(suits)) == 1
    is_royal(c_values::Vector{Int})::Bool = sum(c_values) == 55
    is_consecutive(c_values::Vector{Int}, hand_value::Int)::Bool = sum(c_values) == (5 * (hand_value - 2))
    get_rank_value(counted_values::Dict{Int, Int})::Int = maximum(value for (value, kinds) ∈ counted_values if kinds == maximum(collect(values(counted_values))))

    function get_rank(cards::Vector{SubString{String}})::Tuple{String, Int, Int}
        c_values, suits = map(x -> card_ranks[x], getindex.(cards, 1)), getindex.(cards, 2)
        counted_values, hand_value = countmap(c_values), maximum(c_values)
        unique_values = length(counted_values)
        if unique_values == 5
            if is_flush(suits)
                if is_consecutive(c_values, hand_value)
                    if is_royal(c_values)
                        rank = "Royal Flush"
                    else
                        rank = "Straight Flush"
                    end
                else
                    rank = "Flush"
                end
            elseif is_consecutive(c_values, hand_value)
                rank = "Straight"
            else
                rank = "High Card"
            end
            rank_value = hand_value
        else
            if unique_values == 4
                rank = "One Pair"
            elseif unique_values == 3
                if maximum(collect(values(counted_values))) == 3
                    rank = "Three of a Kind"
                else
                    rank = "Two Pairs"
                end
            else
                if maximum(collect(values(counted_values))) == 4
                    rank = "Four of a Kind"
                else
                    rank = "Full House"
                end
            end
            rank_value = get_rank_value(counted_values)
        end
        return rank, rank_value, hand_value
    end

    function is_winner(hands::Vector{SubString{String}})::Bool
        player_1, player_2 = map(get_rank, [hands[1:5], hands[6:end]])
        if hand_ranks[player_1[1]] > hand_ranks[player_2[1]]
            return true
        elseif hand_ranks[player_1[1]] == hand_ranks[player_2[1]]
            if player_1[2] > player_2[2]
                return true
            elseif player_1[2] == player_2[2]
                return player_1[3] > player_2[3]
            end
        end
        return false
    end

    hand_ranks = Dict(hand => rank for (rank, hand) ∈ enumerate(["High Card", "One Pair",
                                                                 "Two Pairs", "Three of a Kind",
                                                                 "Straight", "Flush", "Full House",
                                                                 "Four of a Kind", "Straight Flush",
                                                                 "Royal Flush"]))
    card_ranks = Dict(card => rank for (rank, card) ∈ enumerate(['2', '3', '4', '5', '6', '7', '8',
                                                                 '9', 'T', 'J', 'Q', 'K', 'A']))
    return sum(map(x -> is_winner(split(x, " ")), split(read(path, String), "\n")))
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
tests = ["reflecttests.jl",
         "datatests.jl",
         "modeltests.jl",
         "objectivetests.jl"
         ]

for test in tests
    include(test)
end


text:
# This file includes plot recipes for Makie

export trisurf

@recipe(Trisurf, msh, f) do scene
    Makie.Attributes(
        wireframe2 = false,
        wfcolor = :black,
        wflinewidth = 2,
        vmarkercolor = :red,
        wflinewidth3 = 3,
        vmarkercolor3 = :orange,
        vmarkersize3 = 20,
        meshcolor3 = :black,
        colormap = :viridis,
        visible = true
    )
end

function Makie.plot!(plt::Trisurf)
    msh3 = plt[1][]
    Makie.mesh!(plt, msh3, color=plt.meshcolor3, colormap=plt.colormap, visible=plt.visible)
    Makie.wireframe!(plt, msh3, linewidth=plt.wflinewidth3)
    plt
end

function Makie.convert_arguments(::Type{<:Trisurf}, msh2::GeometryBasics.Mesh, f)
    msh3 = GeometryBasics.Mesh([Point(pnt..., f(pnt...)) for pnt in msh2.position], faces(msh2))
    return (msh3, msh2, f)
end

function Makie.convert_arguments(plt::Type{<:Trisurf},
                                            pnts::AbstractVector{<:AbstractPoint{2,T}},
                                            f::Function) where T
    Makie.convert_arguments(plt, pnts, map(pnt -> Point(f(pnt...)...), pnts))
end

function Makie.convert_arguments(plt::Type{<:Trisurf},
                                            pnts2d::AbstractVector{<:AbstractPoint{2,T}},
                                            pnts1d::AbstractVector{<:Point1}) where {T}
    Makie.convert_arguments(plt, [Point(pnt2[1], pnt2[2], pnt1[1]) for (pnt2, pnt1) in zip(pnts2d, pnts1d)])
end

function Makie.convert_arguments(::Type{<:Trisurf},
                                            pnts3d::AbstractVector{<:AbstractPoint{3, T}}) where {T}
    pnts2d = project(pnts3d)
    tess = spt.Delaunay(pnts2d)
    fcs = [TriangleFace(val[1], val[2], val[3]) for val in eachrow(tess.simplices .+ 1)]
    msh3 = GeometryBasics.Mesh(pnts3d, fcs)
    return (msh3, pnts3d)
end


text:
#This script illustate how to set up constraints and project a 2D model onto an intersection
# with PARSDMM in serial, parallel or multilevel (serial or parallel)

using Distributed
@everywhere using DistributedArrays, SparseArrays, LinearAlgebra, Random
@everywhere using JOLI
@everywhere using SetIntersectionProjection
using MAT
using PyPlot

@everywhere mutable struct compgrid
  d :: Tuple
  n :: Tuple
end

#PARSDMM options:
options    = PARSDMM_options()
options.FL = Float32
#options=default_PARSDMM_options(options,options.FL)
options.adjust_gamma           = true
options.adjust_rho             = true
options.adjust_feasibility_rho = true
options.Blas_active            = true
options.maxit                  = 500

set_zero_subnormals(true)
BLAS.set_num_threads(2)

#select working precision
if options.FL==Float64
  TF = Float64
  TI = Int64
elseif options.FL==Float32
  TF = Float32
  TI = Int32
end

#load image to project
file = matopen(joinpath(dirname(pathof(SetIntersectionProjection)), "../examples/Data/compass_velocity.mat"));
m    = read(file, "Data");close(file);
m    = m[1:341,200:600];
m    = permutedims(m,[2,1]);

#set up computational grid (25 and 6 m are the original distances between grid points)
comp_grid = compgrid((TF(25.0), TF(6.0)),(size(m,1), size(m,2)))
m         = convert(Vector{TF},vec(m));

#define axis limits and colorbar limits
xmax = comp_grid.d[1]*comp_grid.n[1]
zmax = comp_grid.d[2]*comp_grid.n[2]
vmi  = 1500
vma  = 4500

#constraints
constraint = Vector{SetIntersectionProjection.set_definitions}()

#bounds:
m_min     = 1480.0
m_max     = 4500.0
set_type  = "bounds"
TD_OP     = "identity"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP));

#slope constraints (vertical)
m_min     = 0.0
m_max     = 1e6
set_type  = "bounds"
TD_OP     = "D_z"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP));

options.parallel       = false
(P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL);
(TD_OP,AtA,l,y)        = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options);

println("")
println("PARSDMM serial (bounds and bounds on D_z):")
@time (x1,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);
@time (x1,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);
@time (x1,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);

#plot
figure();imshow(permutedims(reshape(m,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]); title("model to project")
savefig("original_model.png",bbox_inches="tight")
figure();imshow(permutedims(reshape(x1,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]); title("Projection (bounds and bounds on D_z)")
savefig("projected_model.png",bbox_inches="tight")

#plot PARSDMM logs
figure();
subplot(3, 3, 3);semilogy(log_PARSDMM.r_pri)          ;title("r primal")
subplot(3, 3, 4);semilogy(log_PARSDMM.r_dual)         ;title("r dual")
subplot(3, 3, 1);semilogy(log_PARSDMM.obj)            ;title(L"$ \frac{1}{2} || \mathbf{m}-\mathbf{x} ||_2^2 $")
subplot(3, 3, 2);semilogy(log_PARSDMM.set_feasibility);title("TD feasibility violation")
subplot(3, 3, 5);plot(log_PARSDMM.cg_it)              ;title("nr. of CG iterations")
subplot(3, 3, 6);semilogy(log_PARSDMM.cg_relres)      ;title("CG rel. res.")
subplot(3, 3, 7);semilogy(log_PARSDMM.rho)            ;title("rho")
subplot(3, 3, 8);plot(log_PARSDMM.gamma)              ;title("gamma")
subplot(3, 3, 9);semilogy(log_PARSDMM.evol_x)         ;title("x evolution")
tight_layout()
#tight_layout(pad=0.0, w_pad=0.0, h_pad=1.0)
savefig("PARSDMM_logs.png",bbox_inches="tight")

#print timings in terminal
log_PARSDMM.timing

println("")
println("PARSDMM parallel (bounds and bounds on D_z):")
options.parallel       = true
(P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL);
(TD_OP,AtA,l,y)        = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options);

@time (x2,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);
@time (x2,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);
@time (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);

#plot
figure();
subplot(2,1,1);imshow(permutedims(reshape(m,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]); title("model to project")
subplot(2,1,2);imshow(permutedims(reshape(x2,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]); title("Projection (bounds and bounds on D_z)")
savefig("projected_model_ParallelPARSDMM.png",bbox_inches="tight")

#print timings in terminal
log_PARSDMM.timing

#use multilevel-serial (2-levels)
options.parallel = false

#2 levels, the gird point spacing at level 2 is 3X that of the original (level 1) grid
n_levels          = 2
coarsening_factor = 3

#set up all required quantities for each level
#(m_levels,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels)=setup_multi_level_PARSDMM(m,n_levels,coarsening_factor,comp_grid,constraint,options)
(TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels)=setup_multi_level_PARSDMM(m,n_levels,coarsening_factor,comp_grid,constraint,options)

println("")
println("PARSDMM multilevel-serial (bounds and bounds on D_z):")
@time (x,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
@time (x,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
@time (x,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);

figure();
subplot(2,1,1);imshow(permutedims(reshape(m,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]); title("model to project")
subplot(2,1,2);imshow(permutedims(reshape(x,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]); title("Projection (bounds and bounds on D_z)")
savefig("projected_model_MultilevelSerialPARSDMM.png",bbox_inches="tight")

#print timings in terminal, note that for the multi-level version, the timing is only for the final level on the original grid
log_PARSDMM.timing

# #now use multi-level with parallel PARSDMM
# options.parallel=true

# #2 levels, the gird point spacing at level 2 is 3X that of the original (level 1) grid
# n_levels=2
# coarsening_factor=3

# #set up all required quantities for each level
# (TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels)=setup_multi_level_PARSDMM(m,n_levels,coarsening_factor,comp_grid,constraint,options)

# println("")
# println("PARSDMM multilevel-parallel (bounds and bounds on D_z):")
# BLAS.set_num_threads(2)
# @time (x,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
# @time (x,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
# @time (x,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);

# figure();
# subplot(2,1,1);imshow(permutedims(reshape(m,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]); title("model to project")
# subplot(2,1,2);imshow(permutedims(reshape(x,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]); title("Projection (bounds and bounds on D_z)")
# savefig("projected_model_MultilevelParallelPARSDMM.png",bbox_inches="tight")

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
module Geom

using Colors
using Compat
using Compose
using DataArrays
using DataStructures
using Gadfly
using Measures

import Compose.combine # Prevent DataFrame.combine from taking over.
import Gadfly: render, layers, element_aesthetics, inherit, escape_id,
               default_statistic, default_scales, element_coordinate_type,
               ScaleElement, svg_color_class_from_label, isconcrete,
               concretize
import Iterators
import Iterators: cycle, product, distinct, takestrict, chain, repeated

const empty_tag = Symbol("")

function subtags(parent_tag, suffixes...)
    if parent_tag == empty_tag
        return map(s->empty_tag, suffixes)
    end
    return map(s->Symbol(parent_tag, "#", s), suffixes)
end

# Geometry that renders nothing.
immutable Nil <: Gadfly.GeometryElement
end

const nil = Nil

function render(geom::Nil, theme::Gadfly.Theme, aes::Gadfly.Aesthetics,
                data::Gadfly.Data, scales::Dict{Symbol, ScaleElement},
                subplot_layer_aess::Vector{Gadfly.Aesthetics})
end

# Subplot geometries require some more arguments to render. A simpler render
# function is defined and passed through to here for non-subplot geometries.
function render(geom::Gadfly.GeometryElement, theme::Gadfly.Theme, aes::Gadfly.Aesthetics,
                subplot_layer_aess::@compat(Union{(@compat Void), Vector{Gadfly.Aesthetics}}),
                subplot_layer_datas::@compat(Union{(@compat Void), Vector{Gadfly.Data}}),
                scales::Dict{Symbol, ScaleElement})
    render(geom, theme, aes)
end

# Catchall
function default_statistic(::Gadfly.GeometryElement)
    return Gadfly.Stat.identity()
end

include("geom/bar.jl")
include("geom/boxplot.jl")
include("geom/errorbar.jl")
include("geom/hexbin.jl")
include("geom/hline.jl")
include("geom/label.jl")
include("geom/line.jl")
include("geom/point.jl")
include("geom/rectbin.jl")
include("geom/subplot.jl")
include("geom/vline.jl")
include("geom/ribbon.jl")
include("geom/violin.jl")
include("geom/polygon.jl")
include("geom/beeswarm.jl")
include("geom/segment.jl")

end # module Geom


text:
include("approximators.jl")
include("preprocessors.jl")
include("explorers.jl")
include("policies.jl")
include("trajectories.jl")
include("agents.jl")


text:
"""
    ErrorWeights <: AbstractTimeParallelParameters

A composite type for the weights of [`ErrorControl`](@ref).

# Constructors
```julia
ErrorWeights(; w=1.0, updatew=false)
```

# Arguments
- `w::Union{Real, AbstractVector} = ` : weighting of ψ.
- `updatew :: Bool` : to update `w` using local information.

# Functions
- [`show`](@ref) : shows name and contents.
- [`summary`](@ref) : shows name.
"""
mutable struct ErrorWeights{w_T, updatew_T} <: AbstractTimeParallelParameters
    w::w_T
    updatew::updatew_T
end

ErrorWeights(; w=1.0, updatew=false) = ErrorWeights(w, updatew)

#####
##### Functions
#####

function update!(weights::ErrorWeights, U, F)
    @↓ w, updatew = weights
    if updatew
        N = length(U)
        for i = 2:N-1
            # w += norm(F[i+1] - F[i]) / norm(U[i] - U[i-1])
            w = max(w, norm(F[i+1] - F[i]) / norm(U[i] - U[i-1]))
        end
    end
    # w = max(1.0, w / (N-1))
    w = max(1.0, w)
    @↑ weights = w
    return weights
end


text:
using Documenter, Microbiome

makedocs(
    sitename = "Microbiome.jl",
    pages = [
        "Home" => "index.md",
        "Microbial Abundances" => "abundances.md",
        "Distances & Dissimilarity" => "distances.md",
        "Contributing" => "contributing.md"
    ],
    authors = "Kevin Bonham, PhD",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true")
)

deploydocs(
    repo = "github.com/BioJulia/Microbiome.jl.git",
    osname = "linux",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
#
# StabilityCounter
#

counter = StabilityCounter()
@test counter.update(1e-10) == 0
@test counter.update(1e-13) == 1
@test counter.update(1e-13) == 2
@test counter.update(1e-10) == 0

for i=1:149
    counter.update(1e-13)
end

@test counter.ishalt() == false
counter.update(1e-13)
@test counter.ishalt() == true


text:
function findfirstfromlast(arr_Bool)
	for i in Iterators.reverse(eachindex(arr_Bool))
		if arr_Bool[i]
			return i
		end
	end

# 0.5 to 2 times slower:
# function findfirstfromlast1b(arr_Bool);
# revind = Iterators.reverse(eachindex(arr_Bool));
# v = findfirst(arr_Bool[revind]);
# return revind[v]
# end
end


text:
function extract_material(file, ::Type{TF}=Float64) where {TF}
    elastic_heading_pattern = r"\*ELASTIC"
    Emu_pattern = r"(-?\d+\.?\d*)\s*,\s*(-?\d+\.?\d*)"
    line = readline(file)
    m = match(elastic_heading_pattern, line)
    if m != nothing
        line = readline(file)
        m = match(Emu_pattern, line)
        if m != nothing
            E = parse(TF, m[1])
            mu = parse(TF, m[2])
        end
    else
        throw("Material not supported.")
    end
    line = readline(file)
    return E, mu, line
end


text:
module GithubMarkdown

export rendergfm

const libpath = normpath(joinpath(@__DIR__, "..", "deps", "usr", "lib"))
const binary_path = normpath(joinpath(@__DIR__, "..", "deps", "usr", "bin", "cmark-gfm"))

@static if Sys.iswindows()
    const execenv = ("PATH" => string(libpath, ";", Sys.BINDIR))
elseif Sys.isapple()
    const execenv = ("DYLD_LIBRARY_PATH" => libpath)
else
    const execenv = ("LD_LIBRARY_PATH" => libpath)
end

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("LibFoo not installed properly, run Pkg.build(\"LibFoo\"), restart Julia and try again")
end
include(depsjl_path)

# Module initialization function
function __init__()
    # Always check your dependencies from `deps.jl`
    check_deps()
end

"""
    rendergfm(fileout::AbstractString, file::AbstractString; documenter = false, format="html")
    rendergfm(io, file::AbstractString; documenter = false, format="html")

Render the markdown document `file` to `fileout` or `io`, following the cmark-gfm spec.

- `documenter`: Wraps the output in a Documenter `@raw`-block of the specified format.
- `format`: Can be one of `html`, `xml`, `man`, `commonmark`, `plaintext`, `latex`.

Spec: https://github.github.com/gfm/
"""
function rendergfm(fileout::AbstractString, file::AbstractString; kwargs...)
    io = IOBuffer()
    rendergfm(io, file; kwargs...)

    open(fileout, "w") do file
        write(file, seekstart(io))
    end

    return fileout
end

function rendergfm(io, file::AbstractString; documenter = false, format="html")
    isfile(file) || throw(ErrorException("File not found."))

    if !(format in ["html", "xml", "man", "commonmark", "plaintext", "latex"])
        throw(ArgumentError("""
            Invalid format `$(format)`.
            Only `html`, `xml`, `man`, `commonmark`, `plaintext`, `latex` are supported.
        """))
    end

    withenv(execenv) do
        documenter && println(io, "````````````@raw ", format)
        print(io, read(`$(binary_path) $(file) -t $(format) --width 100 -e footnotes -e table -e strikethrough -e tasklist`, String))
        documenter && println(io, "\n````````````")
    end

    return nothing
end

end # module


text:
using ClearStacktrace
ClearStacktrace.CONTRACT_USER_DIR[] = true

function findpair(numbers)
    for i in numbers
        for j in numbers
            if i + j == 2020
                return(i, j)
            end
        end
    end
    error("Couldn't find pair")
end

function findtriplet(numbers)
    for i in numbers
        for j in numbers
            for k in numbers
                if i + j + k == 2020
                    return(i, j, k)
                end
            end
        end
    end
    error("Couldn't find pair")
end

function main()
    numbers = []
    open("inputs/day1.inp") do input
        numbers = sort(parse.(Int, readlines(input)))
    end
    pair = findpair(numbers)
    i, j = pair[1], pair[2]
    println(i * j)

    triplet = findtriplet(numbers)
    i, j, k = triplet[1], triplet[2], triplet[3]
    println(i * j * k)
end

main()


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    BisectFractionPartitioner(normal, fraction=0.5, maxiter=10)

A method for partitioning spatial data into two half spaces
defined by a `normal` direction and a `fraction` of points.
The partition is returned within `maxiter` bisection iterations.
"""
struct BisectFractionPartitioner{T,N} <: AbstractPartitioner
  normal::SVector{N,T}
  fraction::Float64
  maxiter::Int
end

BisectFractionPartitioner(normal::NTuple{N,T},
                          fraction::Real=0.5,
                          maxiter::Integer=10) where {T,N} =
  BisectFractionPartitioner{T,N}(normalize(SVector(normal)), fraction, maxiter)

function partition(object::AbstractSpatialObject{T,N},
                   partitioner::BisectFractionPartitioner{T,N}) where {T,N}
  bbox = boundbox(object)
  n = partitioner.normal
  f = partitioner.fraction
  c = center(bbox)
  d = diagonal(bbox)

  # maximum number of bisections
  maxiter = partitioner.maxiter

  iter = 0; p = 0
  a = c - d/2 * n
  b = c + d/2 * n
  while iter < maxiter
    m = (a + b) / 2

    p = partition(object, BisectPointPartitioner(n, m))
    g = npoints(p[1]) / npoints(object)

    g ≈ f && break
    g > f && (b = m)
    g < f && (a = m)

    iter += 1
  end

  p
end


text:
# Mixed-integer linear programming (MILP) as the expert in IL

"""
    run_MILP!(s::DataFrame, h::Dict) -> JuMP.Model

Assuming all future information is known in the scenario `s`, we run MILP to solve the optimization
problem and get the ideally optimal solutions. The solutions are written into `s` in place, while the
modified `s` is also returned for convenience. The load specifications of the home are presented by `h`.

The JuMP model after optimization is returned.
"""
function run_MILP!(s::DataFrame, h::Dict)::JuMP.Model
    ρ = s.ρ
    PPV = s.PPV

    model = build_model(h, PPV, ρ)
    optimize!(model)
    if termination_status(model) == OPTIMAL
        # @show objective_value(model)
    else
        @warn "Not OPTIMAL"
    end
    write_solution!(s, model, h)
    return model
end

function add_CR_load!(model::JuMP.Model, load::Dict, ρ::AbstractVector)
    # a CR load does not affect optimization
    # we simply store some constants here for convenience
    id = load["id"]
    P = load["P"]::Real
    model[Symbol("P_", id)] = fill(float(P), model[:T])
    model[Symbol("C_", id)] = ρ .* float(P)
    return model
end

"""
    add_AD_load!(model::JuMP.Model, load::Dict, ρ::AbstractVector; K::Integer=10) -> JuMP.Model

Add an adjustable `load` to `model`, whose quadratic cost is approximated by `K` piece linearization.
"""
function add_AD_load!(model::JuMP.Model, load::Dict, ρ::AbstractVector; K::Integer = 10)
    @assert K >= 3
    T = model[:T]
    @assert length(ρ) == T
    Pmin, Pmax = load["Pmin"], load["Pmax"]
    id = load["id"]
    model[Symbol("P_", id)] = P = @variable(model, [1:T], lower_bound = Pmin, upper_bound = Pmax)
    model[Symbol("C_", id)] = C = @variable(model, [1:T])
    # PLF approximation of the quadratic cost
    # note that we cannot separate the electricity cost and the discomfort cost here
    for t = 1:T
        f(P) = ρ[t] * P + load["α"] * (P - Pmax)^2
        xs, ks, bs = linearize(f, Pmin, Pmax, K)
        @constraint(model, C[t] .>= ks .* P[t] .+ bs)
    end
    # store discomfort cost for more info
    model[Symbol("C_", id, "_dc")] = @expression(model, load["α"] .* (P .- Pmax) .^ 2)
    return model
end

"""
    add_SU_load!(model::JuMP.Model, load::Dict, ρ::AbstractVector)

Add an SU `load` to `model` with the price profile `ρ`.
"""
function add_SU_load!(model::JuMP.Model, load::Dict, ρ::AbstractVector)
    T = model[:T]
    L = load["L"]
    ts, tf = load["ts"], load["tf"]
    id = load["id"]
    model[Symbol("z_", id)] = z = @variable(model, [2-L:T], Bin)    # DenseAxisArray
    @constraint(model, z[2-L:ts-1] .== 0)
    @constraint(model, z[tf-L+2:T] .== 0)
    # @constraint(model, sum(z) == 1)  # equiv. to the above two lines
    @constraint(model, sum(z[ts:tf-L+1]) == 1)
    # define P as affine expressions since z is the control variable
    model[Symbol("P_", id)] = P = @expression(model, zeros(AffExpr, T))
    Prate = load["P"]
    @assert Prate > 0
    for t = 1:T
        P[t] = @expression(model, Prate * sum(z[t-L+1:t]))
    end
    model[Symbol("C_", id)] = @expression(model, ρ .* P)
    return model
end

"""
    add_SI_load!(model::JuMP.Model, load::Dict, ρ::AbstractVector)

Add an SI `load` to `model` with the price profile `ρ`.
"""
function add_SI_load!(model::JuMP.Model, load::Dict, ρ::AbstractVector)
    T = model[:T]
    L = load["L"]
    ts, tf = load["ts"], load["tf"]
    id = load["id"]
    Prate = load["P"]
    # only need to define binary variables between ts and tf
    model[Symbol("z_", id)] = z = @variable(model, [1:T], Bin)
    @constraint(model, z[1:ts-1] .== 0)
    @constraint(model, z[tf+1:end] .== 0)
    @constraint(model, sum(z[ts:tf]) == L)
    model[Symbol("P_", id)] = P = @expression(model, Prate .* z)
    model[Symbol("C_", id)] = @expression(model, ρ .* P)
    return model
end

"""
    setup_power_exchange!(model::JuMP.Model, g::Dict, ρ::AbstractVector)

Setup the power exchange with the utility grid configured by `g`. The electricity purchase price
profile is given by `ρ`.
"""
function setup_power_exchange!(model::JuMP.Model, g::Dict, ρ::AbstractVector)
    Pmax = g["Pmax"] # any large value is OK
    β = g["β"]
    @assert 0 <= β <= 1
    T = model[:T]
    @variable(model, zb[1:T], Bin)
    @variable(model, zs[1:T], Bin)
    @variable(model, Pb[1:T] >= 0)
    @variable(model, Ps[1:T] >= 0)
    @constraint(model, Pb .<= zb .* Pmax)
    @constraint(model, Ps .<= zs .* Pmax)
    @constraint(model, zb .+ zs .<= 1)
    ρs = β .* ρ
    @expression(model, Cele, ρ .* Pb - ρs .* Ps)
    return model
end

"""
    build_model(h::Dict, PPV::AbstractVector, ρ::AbstractVector) -> JuMP.Model

Given the PV generation `PPV` and price `ρ` of `T` steps, build an MILP model for HEMS with home
configuration `h`.
"""
function build_model(h::Dict, PPV::AbstractVector, ρ::AbstractVector)::JuMP.Model
    T = length(PPV)
    @assert T == length(ρ)
    model = JuMP.Model(Cbc.Optimizer)
    set_optimizer_attribute(model, "logLevel", 0)
    model[:T] = T
    # add each load (no need to involve the critical load)
    for l in h["CR"]
        add_CR_load!(model, l, ρ)
    end
    for l in h["AD"]
        add_AD_load!(model, l, ρ; K = 20)
    end
    for l in h["SU"]
        add_SU_load!(model, l, ρ)
    end
    for l in h["SI"]
        add_SI_load!(model, l, ρ)
    end
    # add the utility grid
    g = h["UtilityGrid"]
    setup_power_exchange!(model, g, ρ)
    # total load consumption
    @expression(model, P_load, sum(model[Symbol("P_", l["id"])]
                                   for l in Iterators.flatten((h["CR"], h["AD"], h["SU"], h["SI"]))))
    @expression(model, C_load, sum(model[Symbol("C_", l["id"])]
                                   for l in Iterators.flatten((h["CR"], h["AD"], h["SU"], h["SI"]))))
    # set the power balance
    @constraint(model, cst_power_balance, P_load .- PPV .== model[:Pb] .- model[:Ps])
    # set up objective
    ρs = g["β"] .* ρ
    @expression(model, C, (C_load .- ρ .* PPV .+ (ρ .- ρs) .* model[:Ps]))
    @objective(model, Min, sum(C))
    return model
end

"""
    write_solution!(sc::DataFrame, model::JuMP.Model, h::Dict; CR::Bool=false) -> DataFrame

Write the optimal solution of related variables in the optimized `model` into the scenario dataframe
`sc` in place. The updated dataframe is also returned. `h` refers to the home configuration.

The keyword argument `CR` controls whether the constant power and corresponding cost of a critical
load is written into `sc` (default = `false`).
"""
function write_solution!(sc::DataFrame, model::JuMP.Model, h::Dict; CR::Bool = false)
    if termination_status(model) != OPTIMAL
        @warn "The termination status of the model is '$(termination_status(model))'"
    end
    T = model[:T]
    if CR
        for l in h["CR"]
            id = l["id"]
            for pre in ["P_", "C_"]
                key = Symbol(pre, id)
                sc[!, key] = value.(model[key])
            end
        end
    end
    for l in h["AD"]
        id = l["id"]
        for pre in ["P_", "C_"]
            key = Symbol(pre, id)
            sc[!, key] = value.(model[key])
        end
        dc_key = Symbol("C_", id, "_dc")
        sc[!, dc_key] = value.(model[dc_key])
    end
    for l in h["SU"]
        id = l["id"]
        # handle z_ specially because its container is a DenseAxisArray (instead of a Vector)
        zkey = Symbol("z_", id)
        sc[!, zkey] = [value(model[zkey][i]) for i = 1:T]
        for pre in ["P_", "C_"]
            key = Symbol(pre, id)
            sc[!, key] = value.(model[key])[1:T]    # z_ starts from a negative index
        end
    end
    for l in h["SI"]
        id = l["id"]
        for pre in ["z_", "P_", "C_"]
            key = Symbol(pre, id)
            sc[!, key] = value.(model[key])    # z_ starts from a negative index
        end
    end
    for key in (:zb, :zs, :Pb, :Ps, :Cele)
        sc[!, key] = value.(model[key])
    end
    for key in (:P_load, :C_load, :C)
        sc[!, key] = value.(model[key])
    end
    return sc
end

max_stars_repo_name: ShuhuaGao/HybridHEMS

max_issues_repo_name: ShuhuaGao/HybridHEMS

max_forks_repo_name: ShuhuaGao/HybridHEMS


text:
function part1(input)
    terrain = parse_input(input)
    temp = similar(terrain)
    i = 1
    while step!(terrain, temp)
        i += 1
    end
    return i
end

part2(input) = nothing

function parse_input(input)
    return permutedims(reduce(hcat, map_line.(eachline(input))))
end

function map_line(line)
    d = Dict('v' => 1, '>' => 2, '.' => 0)
    return map(x -> d[x], collect(line))
end

function step!(from, temp)
    h, w = size(from)
    change = false
    for j = 1:w
        for i = 1:h
            if from[i, j] == 1
                temp[i, j] = 1
            elseif from[i, j] == 0
                if from[i, mod1(j - 1, w)] == 2
                    change = true
                    temp[i, j] = 2
                else
                    temp[i, j] = 0
                end
            else
                if from[i, mod1(j + 1, w)] == 0
                    temp[i, j] = 0
                else
                    temp[i, j] = 2
                end
            end
        end
    end
    for j = 1:w
        for i = 1:h
            if temp[i, j] == 2
                from[i, j] = 2
            elseif temp[i, j] == 0
                if temp[mod1(i - 1, h), j] == 1
                    change = true
                    from[i, j] = 1
                else
                    from[i, j] = 0
                end
            else
                if temp[mod1(i + 1, h), j] == 0
                    from[i, j] = 0
                else
                    from[i, j] = 1
                end
            end
        end
    end
    return change
end


text:
using MLUtils
using MLUtils.Datasets
using MLUtils: RingBuffer, eachobsparallel
using SparseArrays
using Random, Statistics
using Test
using Transducers
using FoldsThreads: TaskPoolEx
using ChainRulesTestUtils: test_rrule
using Zygote: ZygoteRuleConfig
using ChainRulesCore: rrule_via_ad

showcompact(io, x) = show(IOContext(io, :compact => true), x)

# --------------------------------------------------------------------
# create some test data
Random.seed!(1335)
const X = rand(4, 15)
const y = repeat(["setosa","versicolor","virginica"], inner = 5)
const Y = permutedims(hcat(y,y), [2,1])
const Yt = hcat(y,y)
const yt = Y[1:1,:]
const Xv = view(X,:,:)
const yv = view(y,:)
const XX = rand(20,30,15)
const XXX = rand(3,20,30,15)
const vars = (X, Xv, yv, XX, XXX, y)
const tuples = ((X,y), (X,Y), (XX,X,y), (XXX,XX,X,y))
const Xs = sprand(10, 15, 0.5)
const ys = sprand(15, 0.5)
const X1 = hcat((1:15 for i = 1:10)...)'
const Y1 = collect(1:15)

struct EmptyType end

struct FallbackType end
Base.getindex(::FallbackType, i) = 1234
Base.length(::FallbackType) = 5678

struct CustomType end
MLUtils.numobs(::CustomType) = 15
MLUtils.getobs(::CustomType, i::Int) = i
MLUtils.getobs(::CustomType, i::AbstractVector) = collect(i)

struct CustomArray{T,N} <: AbstractArray{T,N}
    length::Int
    size::NTuple{N, Int}

    function CustomArray{T}(n...) where T
        N = length(n)
        new{T,N}(prod(n), Tuple(n))
    end
end
CustomArray(n...) = CustomArray{Float32}(n...)
Base.length(x::CustomArray) = x.length
Base.size(x::CustomArray) = x.size
Base.getindex(x::CustomArray{T}, i::Int) where T = T(1)
Base.getindex(x::CustomArray{T, 1}, i::AbstractVector{Int}) where T = CustomArray{T}(length(i))

struct CustomArrayNoView{T,N} <: AbstractArray{T,N}
    length::Int
    size::NTuple{N, Int}

    function CustomArrayNoView{T}(n...) where T
        N = length(n)
        new{T,N}(prod(n), Tuple(n))
    end
end
CustomArrayNoView(n...) = CustomArrayNoView{Float32}(n...)
Base.length(x::CustomArrayNoView) = x.length
Base.size(x::CustomArrayNoView) = x.size
Base.getindex(x::CustomArrayNoView{T}, i::Int) where T = T(1)
Base.getindex(x::CustomArrayNoView{T, 1}, i::AbstractVector{Int}) where T = CustomArrayNoView{T}(length(i))
Base.view(x::CustomArrayNoView, i...) = error("view not allowed")

struct CustomRangeIndex
    n::Int
end
Base.length(r::CustomRangeIndex) = r.n
Base.getindex(r::CustomRangeIndex, idx::Int) = idx
Base.getindex(r::CustomRangeIndex, idxs::UnitRange) = idxs

# --------------------------------------------------------------------

include("test_utils.jl")

# @testset "MLUtils.jl" begin

@testset "batchview" begin; include("batchview.jl"); end
@testset "eachobs" begin; include("eachobs.jl"); end
@testset "dataloader" begin; include("dataloader.jl"); end
@testset "folds" begin; include("folds.jl"); end
@testset "observation" begin; include("observation.jl"); end
@testset "obsview" begin; include("obsview.jl"); end
@testset "obstransform" begin; include("obstransform.jl"); end
@testset "randobs" begin; include("randobs.jl"); end
@testset "resample" begin; include("resample.jl"); end
@testset "splitobs" begin; include("splitobs.jl"); end
@testset "utils" begin; include("utils.jl"); end
@testset "eachobsparallel" begin; include("parallel.jl"); end

@testset "Datasets/datasets" begin; include("Datasets/datasets.jl"); end
@testset "Datasets/generators" begin; include("Datasets/generators.jl"); end

# end

max_stars_repo_name: JuliaML/MLBase

max_issues_repo_name: JuliaML/MLBase

max_forks_repo_name: JuliaML/MLBase