text:
"""
Module to store dynamic data. Dynamic data in general is comprised of a set of data points, each of which is associated to
a timestamp. Both the type of data and timestamp can be chosen freely by using the appropriate parameterization.
"""
module DynamicData
using DataStructures
import Base: <, <=, ==, length, isempty, start, next, done, size, eltype, get, get!, getindex, endof

# delegate macro by John Myles White
macro delegate(source, targets)
    typename = esc(source.args[1])
    fieldname = source.args[2].args[1]
    funcnames = targets.args
    n = length(funcnames)
    fdefs = Array(Any, n)
    for i in 1:n
        funcname = esc(funcnames[i])
        fdefs[i] = quote
                     ($funcname)(a::($typename), args...) =
                       ($funcname)(a.$fieldname, args...)
                   end
    end
    return Expr(:block, fdefs...)
end

"""
DataPoint type. A DataPoint consists of the actual data of type D and an associated timestamp of type T.
"""
type DataPoint{D, T}
	id::Int64
	timestamp::T
	value::D
	function DataPoint()
		new(-1, D(), T())
	end
	function DataPoint{D, T}(_t::T, _v::D)
		new(-1, _t, _v)
	end
end
DataPoints{D,T} = Vector{DataPoint{D,T}}

"""
The DynamicDataSet consists of a set of DataPoints of a common value and timestamp types.
Data structures and methods to access DataPoints directly or via their timestamps are provided.
Basic functionality like length and indexing are inherited from Base.
"""
type DynamicDataSet{D, T}
	time_dict::Dict{Int64, T}
	reverse_time_dict::Dict{T, Int64}
	data_points::DataPoints{D, T}
	data_indexes::Dict{T, DataPoints{D,T}}

	function DynamicDataSet()
		time_dict = Dict{Int64, T}()
		reverse_time_dict = Dict{T, Int64}()
		data_points = DataPoints{D, T}()
		data_indexes = Dict{T, DataPoints{D,T}}()
		new(time_dict, reverse_time_dict, data_points, data_indexes)
	end

	function DynamicDataSet(data::Array{D, 1}, timestamps::Array{T, 1})
		time_dict = Dict{Int64, T}()
		reverse_time_dict = Dict{T, Int64}()
		data_points = DataPoints{D, T}()
		data_indexes = Dict{T, DataPoints{D,T}}()
		dynDat = new(time_dict, reverse_time_dict, data_points, data_indexes)

		for (i, dp) in enumerate(data)
			add(dynDat, dp, timestamps[i], false)
		end
		rebuild_indexes(dynDat)
		dynDat
	end
end

@delegate DynamicDataSet.data_points [ get, get!, getindex, length, isempty, eltype, endof, size, start, next, done]

"""
Creates and adds a new DataPoint to the DynamicDataSet based on the provided data of type D, and the associated timestamp of type T.
The parameter ri indicates whether to rebuild the indexes after the insert operation. If you are filling a new dataset with the add method,
set ri to false and explicitly call rebuild_indexes after all insert operations are done for better performance.
"""
function add{D,T}(d::DynamicDataSet{D,T}, data::D, timestamp::T, ri = true)
	#check whether we already have data for the provided timestamp
	if !haskey(d.reverse_time_dict, timestamp)
		#if not, append the timestamp to the reverse_time_dict
		d.reverse_time_dict[timestamp] = length(d.reverse_time_dict) + 1
		d.data_indexes[timestamp] = DataPoints{D,T}()
	end
	#create data point
	dp = DataPoint{D,T}(timestamp, data)
	push!(d.data_points, dp)
	push!(d.data_indexes[timestamp], dp)
	if ri
		rebuild_indexes(d)
	end
	dp
end

function remove!{D,T}(d::DynamicDataSet{D,T}, dp::DataPoint{D,T})
	dp_id = findin(d.data_points, [dp])
	remove!(d, dp_id[1])
end

function remove!{D,T}(d::DynamicDataSet{D,T}, dp_id::Int64)
	d.data_indexes[d.data_points[dp_id].timestamp] = setdiff(d.data_indexes[d.data_points[dp_id].timestamp], [d.data_points[dp_id]])
	deleteat!(d.data_points, dp_id)
end

"""
Rebuilds the indexes after insertion operations.
"""
function rebuild_indexes(d::DynamicDataSet)
	#get all timestamps
	used_timestamps = get_unique_timestamps(d)
	#update the time_dict structure
	for (i, ts) in enumerate(used_timestamps)
		setindex!(d.reverse_time_dict, i, ts)
		setindex!(d.time_dict, ts, i)
	end
end

function convert_datapoint_array{T}(dps::Array{DataPoint{Array{Float64,1},T}, 1})
	ret = zeros(length(dps), length(dps[1].value))
	for (i, dp) in enumerate(dps)
		ret[i,:] = dp.value
	end
	ret
end

"""
Gets all unique timestamps (sorted).
"""
function get_unique_timestamps(d::DynamicData.DynamicDataSet)
	sort(collect(keys(d.reverse_time_dict)))
end

"""
Timestamp getters returning a set of timestamps based on the full data set, a datapoint subset defined by their
IDs in the data set or by an array of elements of type DataPoint.
"""
get_timestamps(d::DynamicDataSet) = map((x) -> x.timestamp, d.data_points)
get_timestamps_for_subset(d::DynamicDataSet, ids::Vector{Int64}) = map((x) -> x.timestamp, d.data_points[ids])
get_timestamps_for_subset{D,T}(pts::Vector{DataPoint{D,T}}) = map((x)->x.timestamp, pts)
get_unique_timestamps_for_subset(d::DynamicDataSet, ids::Vector{Int64}) = unique(get_timestamps_for_subset(d,ids))
get_unique_timestamps_for_subset{D,T}(pts::Vector{DataPoint{D,T}}) = unique(get_timestamps_for_subset(pts))

"""
Gets the timestamp value based on its index.
"""
get_timestamp(d::DynamicDataSet, ti::Int64) = d.time_dict[ti]

"""
Gets the timestamp index based on its value.
"""
get_timestamp_index{T}(d::DynamicDataSet, ts::T) = d.reverse_time_dict[ts]

"""
Gets timestamp index for a given datapoint index.
"""
get_timestamp_index_for_datapoint(d::DynamicDataSet, dp_index::Int64) = d.reverse_time_dict[d.data_points[dp_index].timestamp]

"""
Gets the timestamp indexes of a given index set of data points.
"""
function get_timestamp_indexes_for_subset(d::DynamicDataSet, ids::Vector{Int64})
	timestamp_ids = Array(Int64, length(ids))
	for (i,dp) in enumerate(d[ids])
		timestamp_ids[i] = d.reverse_time_dict[dp.timestamp]
	end
	IntSet(timestamp_ids)
end

"""
Gets the timestamp indexes for a given array of elements of type DataPoint.
"""
function get_timestamp_indexes_for_subset{D,T}(d::DynamicDataSet{D,T}, dps::Vector{DataPoint{D,T}})
	timestamp_ids = Vector{Int64}(length(dps))
	for (i, dp) in enumerate(dps)
		timestamp_ids[i] = d.reverse_time_dict[dp.timestamp]
	end
	IntSet(timestamp_ids)
end

"""
Gets the datapoints associated with a given timestamp.
"""
get_datapoints_for_timestamp{T}(d::DynamicDataSet, ts::T) = d.data_indexes[ts]

"""
Get the datapoint values for an array of elements of type DataPoint.
"""
get_datapoint_values{D,T}(pts::Array{DataPoint{D,T},1}) = map((x)->x.value, pts)

"""
Gets the datapoints associated with the given timestamp based on the provided array of datapoints.
"""
function get_datapoints_for_timestamp_in_subset{D,T}(dps::Array{DataPoint{D,T}, 1}, ts::T)
	ret = Array(DataPoint, 0)
	for dp in dps
		if dp.timestamp == ts
			ret = [ret;dp]
		end
	end
	ret
end

"""
Gets the indexes of datapoints that are associated with a given timestamp based on the provided set of datapoints.
"""
function get_datapoint_idx_for_timestamp_in_subset{D,T}(dps::Array{DataPoint{D,T}, 1}, ts::T)
	ret = Int64[]
	for (i, dp) in enumerate(dps)
		if dp.timestamp == ts
			ret = [ret;i]
		end
	end
	ret
end

"""
Returns an array mapping the timestamp indexes to the number of datapoints associated with them.
"""
function get_datacounts(d::DynamicDataSet)
	sorted_time_ids = sort(collect(keys(d.time_dict)))
	ret = Array(Int64, length(sorted_time_ids))
	for t_id in sorted_time_ids
		ret[t_id] = length(d.data_indexes[d.time_dict[t_id]])
	end
	ret
end

"""
Returns an array mapping the timestamps used in the given subset (provided as DataPoint indexes)
to the number of datapoints associated with them in the subset.
"""
function get_datacounts_for_subset(d::DynamicDataSet, ids::Array{Int64, 1})
	ret = zeros(Int64, length(d.time_dict))
	for id in ids
		t_i = get_timestamp_index_for_datapoint(d, id)
		ret[t_i] +=1
	end
	ret
end

"""
Returns an array mapping the timestamps used in the given subset (provided as an array of DataPoint elements)
to the number of datapoints associated with them in the subset.
"""
function get_datacounts_for_subset{D,T}(d::DynamicDataSet{D,T}, dps::Array{DataPoint{D,T}, 1})
	ret = zeros(Int64, length(d.time_dict))
	for dp in dps
		t_i = get_timestamp_index(d, dp.timestamp)
		ret[t_i] +=1
	end
	ret
end

export DynamicDataSet, DataPoint, add

end


text:
###############################
#
#           WENO5
#
###############################

##---------------------------
## Positive

@doc raw"""
    posWENO5!(flux, U)

Updates the vector "flux" so that flux[i] represents the positive WENO5
flux at x_{i+1/2}, i.e. a WENO5 approximation of the flux going from the
left cell to the right cell. Necessarily, `flux` is supposed to be an
OffsetArray with indices from 0 to N.

The leftmost flux is not computed, as it depends on the boundary
condition. For the stencil to be valid for the second- and
third-to-leftmost fluxes as well as the rightmost and second-to-rightmost
fluxes, Upwind and WENO3 approximations are used instead.
"""
function posWENO5!(flux, U)
    N = length(U)
    flux[1] = posUpwind(U, 1)
    flux[2] = posWENO3(U, 2)
    flux[N-1] = posWENO3(U, N-1)
    flux[N]   = posUpwind(U, N)

    local ε = 10^-8
    ## Smoothness indicators
    local S0 = take(zip(U, drop(U,1), drop(U,2)), N-4) # Stencils for β0
    local S1 = take(zip(drop(U,1), drop(U,2), drop(U,3)), N-4)
    local S2 = zip(drop(U,2), drop(U,3), drop(U,4))

    local γ = 1.08333333333333333333 # ≈ 13/12

    local β0 = [γ*(Ug - 2.0*Um + Ud)^2 + 0.25*(Ug - 4.0*Um + 3.0*Ud)^2
                for (Ug, Um, Ud) ∈ S0]

    local β1 = [γ*(Ug - 2.0*Um + Ud)^2 + 0.25*(Ug - Ud)^2
                for (Ug, Um, Ud) ∈ S1]

    local β2 = [γ*(Ug - 2.0*Um + Ud)^2 + 0.25*(3.0*Ug - 4.0*Um + Ud)^2
                for (Ug, Um, Ud) ∈ S2]

    ## Non-normalized weights
    local ω0 = 0.3 ./ (ε .+ β0).^2
    local ω1 = 0.6 ./ (ε .+ β1).^2
    local ω2 = 0.1 ./ (ε .+ β2).^2
    # normalize them
    local Σα = ω0 + ω1 + ω2
    @. ω0 /= Σα
    @. ω1 /= Σα
    @. ω2 /= Σα

    ## Pointwise values
    local c0 = SA[ 0.333333333333333333, -1.16666666666666667,  1.833333333333333333]
    local c1 = SA[-0.166666666666666667,  0.83333333333333333,  0.333333333333333333]
    local c2 = SA[ 0.333333333333333333,  0.83333333333333333, -0.166666666666666667]

    local U0 = [c0[1]*Ug + c0[2]*Um + c0[2]*Ud for (Ug, Um, Ud) in S0]
    local U1 = [c1[1]*Ug + c1[2]*Um + c1[2]*Ud for (Ug, Um, Ud) in S1]
    local U2 = [c2[1]*Ug + c2[2]*Um + c2[2]*Ud for (Ug, Um, Ud) in S2]

    @. flux[3:N-2] = ω0*U0 + ω1*U1 + ω2*U2

    nothing
end

@doc raw"""
    posWENO5(U,i::Int)

Computes the positive WENO5 flux of U at x_{i+1/2}.

Ref: https://www.cs.usask.ca/faculty/spiteri/63786-gg.pdf
"""
function posWENO5(U,i::Int)
    local ε = 10^-8

    ## Smoothness indicators
    local γ = 1.08333333333333333333 # ≈ 13/12
    local β0 = γ*(U[i-2] - 2.0*U[i-1] + U[i])^2 +
                0.25*(U[i-2] - 4.0*U[i-1] + 3.0*U[i])^2
    local β1 = γ*(U[i-1] - 2.0*U[i] + U[i+1])^2 +
                0.25*(U[i-1] - U[i+1])^2
    local β2 = γ*(U[i] - 2.0*U[i+1] + U[i+2])^2 +
                0.25*(3.0*U[i] - 4.0*U[i+1] + U[i+2])^2

    ## Non-normalized weights
    local α0 = 0.3/(ε + β0)^2
    local α1 = 0.6/(ε + β1)^2
    local α2 = 0.1/(ε + β2)^2
    local Σα = α0 + α1 + α2

    ## Weights (normalized)
    local ω0, ω1, ω2 = α0/Σα, α1/Σα, α2/Σα

    ## Reconstruction of pointwise values
    local c0 = SA[ 0.333333333333333333, -1.16666666666666667,  1.833333333333333333]
    local c1 = SA[-0.166666666666666667,  0.83333333333333333,  0.333333333333333333]
    local c2 = SA[ 0.333333333333333333,  0.83333333333333333, -0.166666666666666667]

    local U0 = c0[1]*U[i-2] + c0[2]*U[i-1] + c0[3]*U[i]
    local U1 = c1[1]*U[i-1] + c1[2]*U[i]   + c1[3]*U[i+1]
    local U2 = c2[1]*U[i]   + c2[2]*U[i+1] + c2[3]*U[i+2]

    ## Final computation
    ω0*U0 + ω1*U1 + ω2*U2
end

##-----------------------
## Negative

@doc raw"""
    negWENO5!(flux, U)

Updates the vector "flux" so that flux[i] represents the negative WENO5
flux at x_{i+1/2}, i.e. a WENO5 approximation of the flux going from the
left cell to the right cell. Necessarily, `flux` is supposed to be an
OffsetArray with indices from 0 to N.

The rightmost flux is not computed, as it depends on the boundary
condition. For the stencil to be valid for the second- and
third-to-rightmost fluxes as well as the leftmost and second-to-leftmost
fluxes, Upwind and WENO3 approximations are used instead.
"""
function negWENO5!(flux, U)
    N = length(U)
    flux[0] = negUpwind(U, 0)
    flux[1] = negWENO3(U, 1)
    flux[N-2] = negWENO3(U, N-2)
    flux[N-1] = negUpwind(U, N-1)

    @inbounds for i in 2 : N-3
        flux[i] = negWENO5(U, i)
    end

    nothing
end

@doc raw"""
    negWENO5(U,i::Int)

Computes the negative WENO5 flux of U at x_{i+1/2}. This could also be
obtained by symmetry from posWENO5.

Ref: https://www.cs.usask.ca/faculty/spiteri/63786-gg.pdf
"""
function negWENO5(U, i::Int)
    local ε = 10^-8

    ## Smoothness indicators
    local γ = 1.08333333333333333333 # = 13/12
    local β0 = γ*(U[i+3] - 2.0*U[i+2] + U[i+1])^2 +
                0.25*(U[i+3] - 4.0*U[i+2] + 3.0*U[i+1])^2
    local β1 = γ*(U[i+2] - 2.0*U[i+1] + U[i])^2 +
                0.25*(U[i+2] - U[i])^2
    local β2 = γ*(U[i+1] - 2.0*U[i] + U[i-1])^2 +
                0.25*(3.0*U[i+1] - 4.0*U[i] + U[i-1])^2

    ## Non-normalized weights
    local α0 = 0.3/(ε + β0)^2
    local α1 = 0.6/(ε + β1)^2
    local α2 = 0.1/(ε + β2)^2
    local Σα = α0 + α1 + α2

    ## Weights (normalized coefficients)
    local ω0, ω1, ω2 = α0/Σα, α1/Σα, α2/Σα

    ## Reconstructions of pointwise values
    # local c0 = SA[ 1.833333333333333333, -1.16666666666666667,  0.333333333333333333]
    # local c1  = SA[ 0.333333333333333333,  0.83333333333333333, -0.166666666666666667]
    # local c2  = SA[-0.166666666666666667,  0.83333333333333333,  0.333333333333333333]

    # local U0 = c0[1]*U[i+1] + c0[2]*U[i+2] + c0[3]*U[i+3]
    # local U1 = c1[1]*U[i]   + c1[2]*U[i+1] + c1[3]*U[i+2]
    # local U2 = c2[1]*U[i-1] + c2[2]*U[i]   + c2[3]*U[i+1]

    local c0 = SA[ 0.333333333333333333, -1.16666666666666667,  1.833333333333333333]
    local c1 = SA[-0.166666666666666667,  0.83333333333333333,  0.333333333333333333]
    local c2 = SA[ 0.333333333333333333,  0.83333333333333333, -0.166666666666666667]

    local U0 = c0[1]*U[i+3] + c0[2]*U[i+2] + c0[3]*U[i+1]
    local U1 = c1[1]*U[i+2] + c1[2]*U[i+1] + c1[3]*U[i]
    local U2 = c2[1]*U[i+1] + c2[2]*U[i]   + c2[3]*U[i-1]

    ## Final computation
    ω0*U0 + ω1*U1 + ω2*U2
end

ext: jl

max_stars_repo_name: tremelow/StiffKinetic

max_issues_repo_name: tremelow/StiffKinetic

max_forks_repo_name: tremelow/StiffKinetic


text:
module TestGrouping

using Test
using DataFrames
using DataFramesMeta
using Statistics
using Tables

const ≅ = isequal

d = DataFrame(n = 1:20, x = [3, 3, 3, 3, 1, 1, 1, 2, 1, 1, 2, 1, 1, 2, 2, 2, 3, 1, 1, 2])
g = groupby(d, :x, sort=true)

@test  @where(d, :x .== 3) == DataFramesMeta.where(d, x -> x.x .== 3)
@test  DataFrame(@where(g, length(:x) > 5)) == DataFrame(DataFramesMeta.where(g, x -> length(x.x) > 5))
@test  DataFrame(@where(g, length(:x) > 5))[!, :n][1:3] == [5, 6, 7]

@test  DataFrame(DataFramesMeta.orderby(g, x -> mean(x.n))) == DataFrame(@orderby(g, mean(:n)))

@test @based_on(g, nsum = sum(:n)).nsum == [99, 84, 27]

@testset "@based_on" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing]
        )

    m = [100, 200, 300, 400, 500]

    gq = :g
    iq = :i
    tq = :t
    yq = :y
    cq = :c

    gr = "g"
    ir = "i"
    tr = "t"
    yr = "y"
    cr = "c"

    gd = groupby(df, :g)

    n_str = "new_column"
    n_sym = :new_column
    n_space = "new column"

    @test @based_on(gd, n = mean(:i)).n == [2.0, 4.5]
    @test @based_on(gd, n = mean(:i) + mean(:g)).n == [3.0, 6.5]
    @test @based_on(gd, n = first(:t .* string.(:y))).n == ["av", "cy"]
    @test @based_on(gd, n = first(Symbol.(:y, ^(:t)))).n == [:vt, :yt]
    @test @based_on(gd, n = first(Symbol.(:y, ^(:body)))).n == [:vbody, :ybody]
    @test @based_on(gd, body = :i).body == df.i
    @test @based_on(gd, transform = :i).transform == df.i
    @test @based_on(gd, (n1 = [first(:i)], n2 = [first(:y)])).n1 == [1, 4]

    @test @based_on(gd, n = mean(cols(iq))).n == [2.0, 4.5]
    @test @based_on(gd, n = mean(cols(iq)) + mean(cols(gq))).n == [3.0, 6.5]
    @test @based_on(gd, n = first(cols(tq) .* string.(cols(yq)))).n == ["av", "cy"]
    @test @based_on(gd, n = first(Symbol.(cols(yq), ^(:t)))).n == [:vt, :yt]
    @test @based_on(gd, n = first(Symbol.(cols(yq), ^(:body)))).n == [:vbody, :ybody]
    @test @based_on(gd, body = cols(iq)).body == df.i
    @test @based_on(gd, transform = cols(iq)).transform == df.i
    @test @based_on(gd, (n1 = [first(cols(iq))], n2 = [first(cols(yq))])).n1 == [1, 4]

    @test @based_on(gd, n = mean(cols(ir))).n == [2.0, 4.5]
    @test @based_on(gd, n = mean(cols(ir)) + mean(cols(gr))).n == [3.0, 6.5]
    @test @based_on(gd, n = first(cols(tr) .* string.(cols(yr)))).n == ["av", "cy"]
    @test @based_on(gd, n = first(Symbol.(cols(yr), ^(:t)))).n == [:vt, :yt]
    @test @based_on(gd, n = first(Symbol.(cols(yr), ^(:body)))).n == [:vbody, :ybody]
    @test @based_on(gd, body = cols(ir)).body == df.i
    @test @based_on(gd, transform = cols(ir)).transform == df.i
    @test @based_on(gd, (n1 = [first(cols(ir))], n2 = [first(cols(yr))])).n1 == [1, 4]

    @test @based_on(gd, :i) == select(df, :g, :i)
    @test @based_on(gd, :i, :g) ≅ select(df, :g, :i)

    @test @based_on(gd, :i, n = 1).n == fill(1, nrow(df))

    @test @based_on(gd, cols("new_column") = 2).new_column == [2, 2]
    @test @based_on(gd, cols(n_str) = 2).new_column == [2, 2]
    @test @based_on(gd, cols(n_sym) = 2).new_column == [2, 2]
    @test @based_on(gd, cols(n_space) = 2)."new column" == [2, 2]
    @test @based_on(gd, cols("new" * "_" * "column") = 2)."new_column" == [2, 2]
end

# Defined outside of `@testset` due to use of `@eval`
df = DataFrame(
    g = [1, 1, 1, 2, 2],
    i = 1:5,
    t = ["a", "b", "c", "c", "e"],
    y = [:v, :w, :x, :y, :z],
    c = [:g, :quote, :body, :transform, missing]
    )

m = [100, 200, 300, 400, 500]

gq = :g
iq = :i
tq = :t
yq = :y
cq = :c

gr = "g"
ir = "i"
tr = "t"
yr = "y"
cr = "c"

gd = groupby(df, :g)

newvar = :n

@testset "Limits of @based_on" begin
    t = @based_on(gd, [:i, :g]).i_g_function
    @test t == [[1, 2, 3], [1, 1, 1], [4, 5], [2, 2]]
    @test t isa Vector{SubArray{Int64,1,Array{Int64,1},Tuple{Array{Int64,1}},false}}
    @test @based_on(gd, All()).function isa Vector{<:All}
    @test @based_on(gd, Not(:i)).i_function isa Vector{<:InvertedIndex}
    @test @based_on(gd, Not([:i, :g])).g == [1, 2]
    @test_throws MethodError @eval @based_on(gd, n = sum(Between(:i, :t)))
    @test_throws LoadError @eval @based_on(gd; n = mean(:i))
end

@testset "@by" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing]
        )

    m = [100, 200, 300, 400, 500]

    gq = :g
    iq = :i
    tq = :t
    yq = :y
    cq = :c

    gr = "g"
    ir = "i"
    tr = "t"
    yr = "y"
    cr = "c"

    gd = groupby(df, :g)

    n_str = "new_column"
    n_sym = :new_column
    n_space = "new column"

    @test @by(df, :g, n = mean(:i)).n == [2.0, 4.5]
    @test @by(df, :g, n = mean(:i) + mean(:g)).n == [3.0, 6.5]
    @test @by(df, :g, n = first(:t .* string.(:y))).n == ["av", "cy"]
    @test @by(df, :g, n = first(Symbol.(:y, ^(:t)))).n == [:vt, :yt]
    @test @by(df, :g, n = first(Symbol.(:y, ^(:body)))).n == [:vbody, :ybody]
    @test @by(df, :g, body = :i).body == df.i
    @test @by(df, :g, transform = :i).transform == df.i
    @test @by(df, :g, (n1 = [first(:i)], n2 = [first(:y)])).n1 == [1, 4]

    @test @by(df, :g, n = mean(cols(iq))).n == [2.0, 4.5]
    @test @by(df, :g, n = mean(cols(iq)) + mean(cols(gq))).n == [3.0, 6.5]
    @test @by(df, :g, n = first(cols(tq) .* string.(cols(yq)))).n == ["av", "cy"]
    @test @by(df, :g, n = first(Symbol.(cols(yq), ^(:t)))).n == [:vt, :yt]
    @test @by(df, :g, n = first(Symbol.(cols(yq), ^(:body)))).n == [:vbody, :ybody]
    @test @by(df, :g, body = cols(iq)).body == df.i
    @test @by(df, :g, transform = cols(iq)).transform == df.i
    @test @by(df, :g, (n1 = [first(cols(iq))], n2 = [first(cols(yq))])).n1 == [1, 4]

    @test @by(df, "g", n = mean(cols(ir))).n == [2.0, 4.5]
    @test @by(df, "g", n = mean(cols(ir)) + mean(cols(gr))).n == [3.0, 6.5]
    @test @by(df, "g", n = first(cols(tr) .* string.(cols(yr)))).n == ["av", "cy"]
    @test @by(df, "g", n = first(Symbol.(cols(yr), ^(:t)))).n == [:vt, :yt]
    @test @by(df, "g", n = first(Symbol.(cols(yr), ^(:body)))).n == [:vbody, :ybody]
    @test @by(df, "g", body = cols(ir)).body == df.i
    @test @by(df, "g", transform = cols(ir)).transform == df.i
    @test @by(df, "g", (n1 = [first(cols(ir))], n2 = [first(cols(yr))])).n1 == [1, 4]

    @test @by(df, :g, :i) == select(df, :g, :i)
    @test @by(df, :g, :i, :g) ≅ select(df, :g, :i)

    @test @by(df, :g, :i, n = 1).n == fill(1, nrow(df))

    @test @by(df, :g, cols("new_column") = 2).new_column == [2, 2]
    @test @by(df, :g, cols(n_str) = 2).new_column == [2, 2]
    @test @by(df, :g, cols(n_sym) = 2).new_column == [2, 2]
    @test @by(df, :g, cols(n_space) = 2)."new column" == [2, 2]
    @test @by(df, :g, cols("new" * "_" * "column") = 2)."new_column" == [2, 2]
end

# Defined outside of `@testset` due to use of `@eval`
df = DataFrame(
    g = [1, 1, 1, 2, 2],
    i = 1:5,
    t = ["a", "b", "c", "c", "e"],
    y = [:v, :w, :x, :y, :z],
    c = [:g, :quote, :body, :transform, missing]
    )

m = [100, 200, 300, 400, 500]

gq = :g
iq = :i
tq = :t
yq = :y
cq = :c

gr = "g"
ir = "i"
tr = "t"
yr = "y"
cr = "c"

gd = groupby(df, :g)

newvar = :n

@testset "limits of @by" begin
    t = @by(df, :g, [:i, :g]).i_g_function
    @test t == [[1, 2, 3], [1, 1, 1], [4, 5], [2, 2]]
    @test t isa Vector{SubArray{Int64,1,Array{Int64,1},Tuple{Array{Int64,1}},false}}
    @test @by(df, :g, All()).function isa Vector{<:All}
    @test @by(df, :g, Not(:i)).i_function isa Vector{<:InvertedIndex}
    @test @by(df, :g, Not([:i, :g])).g == [1, 2]
    @test_throws MethodError @eval @by(df, :g, n = sum(Between(:i, :t)))
    @test_throws MethodError @eval @by(df, :g; n = mean(:i))
end

@testset "@transform with grouped data frame" begin
	d = DataFrame(n = 1:20, x = [3, 3, 3, 3, 1, 1, 1, 2, 1, 1, 2, 1, 1, 2, 2, 2, 3, 1, 1, 2])
	g = groupby(d, :x)

	@test (@transform(g, y = :n .- median(:n)))[1,:y] == -2.0

	d = DataFrame(a = [1,1,1,2,2,3,3,1],
	              b = Any[1,2,3,missing,missing,6.0,5.0,4],
	              c = CategoricalArray([1,2,3,1,2,3,1,2]))
	g = groupby(d, :a)

      ## Scalar output
	# Type promotion Int -> Float
	t = @transform(g, t = :b[1]).t
	s = @select(g, t = :b[1]).t
	@test t ≅ s ≅ [1.0, 1.0, 1.0, missing, missing, 6.0, 6.0, 1.0] &&
	      t isa Vector{Union{Float64, Missing}}

	# Type promotion Number -> Any
	t = @transform(g, t = isequal(:b[1], 1) ? :b[1] : "a").t
	s = @select(g, t = isequal(:b[1], 1) ? :b[1] : "a").t
	@test t ≅ s ≅ [1, 1, 1, "a", "a", "a", "a", 1] &&
	      t isa Vector{Any}
	## Vector output
	# Normal use
	t = @transform(g, t = :b .- mean(:b)).t
	s = @select(g, t = :b .- mean(:b)).t
	@test t ≅ s ≅ [-1.5, -0.5, 0.5, missing, missing, 0.5, -0.5, 1.5] &&
	      t isa Vector{Union{Float64, Missing}}
	# Type promotion
	t = @transform(g, t = isequal(:b[1], 1) ? fill(1, length(:b)) : fill(2.0, length(:b))).t
	s = @transform(g, t = isequal(:b[1], 1) ? fill(1, length(:b)) : fill(2.0, length(:b))).t
	@test t ≅ s ≅ [1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 2.0, 1.0] &&
	      t isa Vector{Float64}
	# Vectors whose eltypes promote to any
	t = @transform(g, t = isequal(:b[1], 1) ? :b : fill("a", length(:b))).t
	s = @transform(g, t = isequal(:b[1], 1) ? :b : fill("a", length(:b))).t
	@test s ≅ t ≅ [1, 2, 3, "a", "a", "a", "a", 4] &&
	      t isa Vector{Any}
	# Categorical Array
	# Scalar
	t = @transform(g, t = :c[1]).t
	s = @transform(g, t = :c[1]).t
	@test t ≅ s ≅  [1, 1, 1, 1, 1, 3, 3, 1] &&
	      t isa CategoricalVector{Int}
	# Vector
	t = @transform(g, t = :c).t
	s = @transform(g, t = :c).t
	@test t ≅ s ≅ [1, 2, 3, 1, 2, 3, 1, 2] &&
	      t isa CategoricalVector{Int}

	@test @transform(g, t = :c).a ≅ d.a
	@test @select(g, :a, t = :c).a ≅ d.a
end

end # module


text:
# This does MLE and then MCMC, either using raw statistic, or using NN transform,
# depending on the argument usenn
using Flux, Econometrics, LinearAlgebra, Statistics, DelimitedFiles
using BSON:@load
include("transform.jl")
include("lnL.jl")

# uniform random walk in one dimension
function proposal1(current, tuning)
    trial = copy(current)
    i = rand(1:size(trial,1))
    trial[i] += tuning[i].*randn()
    return trial
end

# MVN random walk, or occasional draw from prior
function proposal2(current, cholV)
    current + cholV'*randn(size(current))
end

function MCMC(m, usenn, info, useJacobian=true)
    S = 100 # number of simulations
    lb, ub = PriorSupport()
    nParams = size(lb,1)
    # get the trained net
    if usenn
        @load "best.bson" model
        m = transform(m', info)
        m = Float64.(model(m'))
        θinit = PriorMean() # prior mean as initial θ
        lnL = θ -> LL(θ, m, S, model, info, useJacobian)
    else
        θinit = PriorMean() # prior mean as initial θ
        lnL = θ -> LL(θ, m, S, useJacobian)
    end
    # use a rapid SAMIN to get good initialization values for chain
    obj = θ -> -1.0*lnL(θ)
    θmile, junk, junk, junk = samin(obj, θinit, lb, ub; coverage_ok=0, maxevals=100000, verbosity = 0, rt = 0.5)
    prior = θ -> Prior(θ) # uniform, doesn't matter
    # define things for MCMC
    verbosity = false
    ChainLength = 1000
    MCMCburnin = 100
    tuning = 0.2/sqrt(12.0)*(ub-lb) # two tenths of a standard. dev. of prior
    Proposal = θ -> proposal1(θ, tuning)
    chain = mcmc(θmile, ChainLength, MCMCburnin, prior, lnL, Proposal, verbosity)
    # now use a MVN random walk proposal with updates of covariance and longer chain
    # on final loop
    Σ = NeweyWest(chain[:,1:nParams])
    tuning = 1.0
    MC_loops = 5
    @inbounds for j = 1:MC_loops
        P = try
            P = (cholesky(Σ)).U
        catch
            P = diagm(diag(Σ))
        end
        Proposal = θ -> proposal2(θ,tuning*P)
        if j == MC_loops
            ChainLength = 10000
        end
        θinit = mean(chain[:,1:nParams],dims=1)[:]
        chain = mcmc(θinit, ChainLength, 0, prior, lnL, Proposal, verbosity)
        if j < MC_loops
            accept = mean(chain[:,end])
            if accept > 0.35
                tuning *= 1.5
            elseif accept < 0.25
                tuning *= 0.25
            end
            Σ = 0.5*Σ + 0.5*NeweyWest(chain[:,1:nParams])
        end
    end
    chain = chain[:,1:nParams]
    return chain, θmile
end

max_stars_repo_name: wenddymacro/Econometrics

max_issues_repo_name: wenddymacro/Econometrics

max_forks_repo_name: wenddymacro/Econometrics


text:
if basename(pwd()) == "aoc"
    cd("2015/15")
end

function loadingredients(filename::AbstractString)
    map(eachline(filename)) do line
        parse.(Int, (m.match for m in eachmatch(r"-?\d+", line)))
    end
end

function part1(filename::AbstractString)
    ingredients = loadingredients(filename)
    foreach(pop!, ingredients)
    maximum(score(ingredients, p) for p in eachproportion(length(ingredients)))
end

function score(ingredients, proportions)
    prod(max.(0, sum(ingredients .* proportions)))
end

function eachproportion(len, max = 100)
    if len == 1
        max
    else
        ([i, p...] for i in 0:max for p in eachproportion(len - 1, max - i))
    end
end

function part2(filename::AbstractString)
    ingredients = loadingredients(filename)
    calories = pop!.(ingredients)
    maximum(
        score(ingredients, p)
        for p in eachproportion(length(ingredients))
        if score(calories, p) == 500
    )
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
# # Modal analysis of Argyris frame: effect of prestress

# Source code: [`argyris_frame_modal_tut.jl`](argyris_frame_modal_tut.jl)

# ## Description

# Vibration analysis of a L-shaped frame under a loading.
# The fundamental vibration frequency depends on the prestress force.

# ## Goals

# - Construct an L-shaped frame by merging individual members.
# - Compute the geometric stiffness.
# - Evaluate the effect of prestress on the fundamental frequency of vibration.

# The finite element code relies on the basic functionality implemented in this
# package.
using FinEtools
# The linear deformation code will be needed to evaluate the loading.
using FinEtoolsDeforLinear
# The functionality concerning beam models comes from these modules.
using FinEtoolsFlexBeams.CrossSectionModule: CrossSectionRectangle
using FinEtoolsFlexBeams.RotUtilModule:  update_rotation_field!
using FinEtoolsFlexBeams.MeshFrameMemberModule: frame_member, merge_members
using FinEtoolsFlexBeams.RotUtilModule: initial_Rfield
using FinEtoolsFlexBeams.FEMMCorotBeamModule: FEMMCorotBeam
using FinEtoolsFlexBeams.FEMMCorotBeamModule
CB = FEMMCorotBeamModule

# Parameters:
E = 71240.0 * phun("MPa")
nu = 0.31; # Poisson ratio
rho = 5000 * phun("kg/m^3");
# cross-sectional dimensions and length of each leg in millimeters
b = 0.6 * phun("mm"); h = 30.0 * phun("mm"); L = 240.0 * phun("mm");
# Magnitude of the total applied force, Newton
magn = 1e-5 * phun("N");

# Cross-sectional properties
cs = CrossSectionRectangle(s -> b, s -> h, s -> [0.0, 1.0, 0.0])

##
# ## Generate the discrete model

# Select the number of elements per leg.
n = 8;
members = Tuple{FENodeSet, AbstractFESet}[]
push!(members, frame_member([0 0 L; L 0 L], n, cs))
push!(members, frame_member([L 0 L; L 0 0], n, cs))
fens, fes = merge_members(members; tolerance = L / 10000)

# Construct the requisite fields, geometry and displacement
# Initialize configuration variables
geom0 = NodalField(fens.xyz)
u0 = NodalField(zeros(size(fens.xyz,1), 3))
Rfield0 = initial_Rfield(fens)
dchi = NodalField(zeros(size(fens.xyz,1), 6))

# Apply EBC's
l1 = selectnode(fens; box = [0 0 0 0 L L], tolerance = L / 10000)
for i in [1, 2, 3, 4, 5, 6]
    setebc!(dchi, l1, true, i)
end
applyebc!(dchi)
numberdofs!(dchi);

# Assemble the global discrete system

# Material properties
material = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

femm = FEMMCorotBeam(IntegDomain(fes, GaussRule(1, 2)), material)
K = CB.stiffness(femm, geom0, u0, Rfield0, dchi);
M = CB.mass(femm, geom0, u0, Rfield0, dchi);

# Construct force intensity,  loaded boundary, and assemble the load.
tipn = selectnode(fens; box=[L L 0 0  0 0], tolerance=L/n/1000)[1]
loadbdry = FESetP1(reshape([tipn], 1, 1))
lfemm = FEMMBase(IntegDomain(loadbdry, PointRule()))
fi = ForceIntensity(FFlt[-magn, 0, 0, 0, 0, 0]);
F = CB.distribloads(lfemm, geom0, dchi, fi, 3);

# Solve for the displacement under the static load.
scattersysvec!(dchi, K\F);

# Update deflections and rotations so that the initial stress can be computed.
# First the displacements:
u1 = deepcopy(u0)
u1.values .= dchi.values[:, 1:3]
# Then the rotations:
Rfield1 = deepcopy(Rfield0)
update_rotation_field!(Rfield1, dchi)

# The static deflection is now used to compute the internal forces
# which in turn lead to the geometric stiffness.
Kg = CB.geostiffness(femm, geom0, u1, Rfield1, dchi);

##
# ## Solution of the eigenvalue free-vibration problem

using Arpack

# We will solve for this many natural frequencies. Then, since they are ordered
# by magnitude, we will pick the fundamental by taking the first from the
# list.
neigvs = 4

# First we will  sweep through the loading factors that are positive, meaning
# the force points in the direction in which it was defined (towards the
# clamped end of the frame).

lfp = linearspace(0.0, 68000.0, 400)
fsp = let
    fsp = Float64[]
    for load_factor in lfp
        evals, evecs, nconv = eigs(K + load_factor .* Kg, M; nev=neigvs, which=:SM);
        f = evals[1] > 0 ? sqrt(evals[1]) / (2 * pi) : 0;
        push!(fsp, f)
    end
    fsp
end

# Next, we will sweep through arrange of negative load factors: this simply
# turns the force around so that it points away from the clamped end. This can
# also buckle the frame, but the magnitude is higher.
lfm = linearspace(-109000.0, 0.0, 400)
fsm = let
    fsm = Float64[]
    for load_factor in lfm
        evals, evecs, nconv = eigs(K + load_factor .* Kg, M; nev=neigvs, which=:SM);
        f = evals[1] > 0 ? sqrt(evals[1]) / (2 * pi) : 0;
        push!(fsm, f)
    end
    fsm
end

##
# ## Plot of the fundamental frequency is it depends on the loading factor

using Gnuplot

# We concatenate the ranges for the load factors and the calculated fundamental
# frequencies and present them in a single plot.

Gnuplot.gpexec("reset session")
@gp  "set terminal windows 0 "  :-

@gp  :- cat(collect(lfp), collect(lfm); dims=1) cat(fsp, fsm; dims=1) " lw 2 lc rgb 'red' with p title 'Fundamental frequency' "  :-

@gp  :- "set xlabel 'Loading factor P'" :-
@gp  :- "set ylabel 'Frequency(P) [Hz]'" :-
@gp  :- "set title 'Frame fundamental frequency'"

# Clearly, the curve giving the dependence of the fundamental frequency on the
# loading factor consists of two branches. These two branches correspond to two
# different buckling modes: one for the positive orientation of the force and
# one for the negative orientation.

##
# ## Visualize some fundamental mode shapes

# Here we visualize the fundamental vibration modes for different values of the
# loading factor.

using PlotlyJS
using FinEtoolsFlexBeams.VisUtilModule: plot_space_box, plot_solid, render, react!, default_layout_3d, save_to_json
scale = 0.005

vis(loading_factor, evec) = let
    tbox = plot_space_box(reshape(inflatebox!(boundingbox(fens.xyz), 0.3 * L), 2, 3))
    tenv0 = plot_solid(fens, fes; x=geom0.values, u=0.0 .* dchi.values[:, 1:3], R=Rfield0.values, facecolor="rgb(125, 155, 125)", opacity=0.3);
    plots = cat(tbox, tenv0; dims=1)
    layout = default_layout_3d(;width=600, height=600)
    layout[:scene][:aspectmode] = "data"
    pl = render(plots; layout=layout, title = "Loading factor $(loading_factor)")
    for xscale in scale .* sin.(collect(0:1:89) .* (2 * pi / 21))
        scattersysvec!(dchi, xscale .* evec)
        u1 = deepcopy(u0)
        u1.values .= dchi.values[:, 1:3]
        Rfield1 = deepcopy(Rfield0)
        update_rotation_field!(Rfield1, dchi)
        tenv1 = plot_solid(fens, fes; x=geom0.values, u=dchi.values[:, 1:3], R=Rfield1.values, facecolor="rgb(50, 55, 125)");
        plots = cat(tbox, tenv0, tenv1; dims=1)
        react!(pl, plots, pl.plot.layout)
        sleep(0.115)
    end
end

# This is the vibration mode in the lead up to the buckling mode for the
# positive orientation of the force.

loading_factor = 60000
evals, evecs, nconv = eigs(K + loading_factor .* Kg, M; nev=neigvs, which=:SM);
vis(loading_factor, evecs[:, 1])

# This is the same vibration mode for the negative orientation of the force, but
# note that the associated fundamental frequency increased due to the effect of
# the force upon the stiffening of the clamped leg of the frame that is now in
# tension, and therefore stiffer.

loading_factor = -50000
evals, evecs, nconv = eigs(K + loading_factor .* Kg, M; nev=neigvs, which=:SM);
vis(loading_factor, evecs[:, 1])

# Increasing the load factor in the negative orientation further, the
# fundamental frequency will switch: it will be a different mode shape, the one
# that is close to the buckling mode shape for this orientation of the force.

loading_factor = -100000
evals, evecs, nconv = eigs(K + loading_factor .* Kg, M; nev=neigvs, which=:SM);
vis(loading_factor, evecs[:, 1])

true


text:
Base.lock(h::Hist2D) = lock(h.hlock)
Base.unlock(h::Hist2D) = unlock(h.hlock)

"""
    bincounts(h::Hist2D)

Get the bin counts of a histogram.
"""
@inline bincounts(h::Hist2D) = h.hist.weights

"""
    binedges(h::Hist2D)

Get a 2-tuple of the bin edges of a histogram.
"""
@inline binedges(h::Hist2D) = h.hist.edges

"""
    bincenters(h::Hist2D)

Get a 2-tuple of the bin centers of a histogram.
"""
function bincenters(h::Hist2D)
    StatsBase.midpoints.(binedges(h))
end

"""
    binerrors(f::T, h::Hist2D) where T<:Function = f.(h.sumw2)
    binerrors(h::Hist2D) = binerrors(sqrt, h)

Calculate the bin errors from `sumw2` with a Gaussian default.
"""
binerrors(f::T, h::Hist2D) where T<:Function = f.(h.sumw2)
binerrors(h::Hist2D) = binerrors(sqrt, h)

"""
    nbins(h::Hist2D)

Get a 2-tuple of the number of x and y bins of a histogram.
"""
function nbins(h::Hist2D)
    size(bincounts(h))
end

"""
    integral(h::Hist2D)

Get the integral a histogram.
"""
function integral(h::Hist2D)
    sum(bincounts(h))
end

"""
    empty!(h::Hist2D)

Resets a histogram's bin counts and `sumw2`.
"""
function Base.empty!(h::Hist2D{T,E}) where {T,E}
    h.hist.weights .= zero(T)
    h.sumw2 .= 0.0
    return h
end

"""
    push!(h::Hist2D, valx::Real, valy::Real, wgt::Real=1)
    atomic_push!(h::Hist2D, valx::Real, valy::Real, wgt::Real=1)

Adding one value at a time into histogram.
`sumw2` (sum of weights^2) accumulates `wgt^2` with a default weight of 1.
`atomic_push!` is a slower version of `push!` that is thread-safe.

"""
@inline function atomic_push!(h::Hist2D{T,E}, valx::Real, valy::Real, wgt::Real=1) where {T,E}
    lock(h)
    push!(h, valx, valy, wgt)
    unlock(h)
    return nothing
end

@inline function Base.push!(h::Hist2D{T,E}, valx::Real, valy::Real, wgt::Real=1) where {T,E}
    rx, ry = binedges(h)
    Lx, Ly = nbins(h)
    binidxx = _edge_binindex(rx, valx)
    binidxy = _edge_binindex(ry, valy)
    if h.overflow
        binidxx = clamp(binidxx, 1, Lx)
        binidxy = clamp(binidxy, 1, Ly)
        @inbounds h.hist.weights[binidxx,binidxy] += wgt
        @inbounds h.sumw2[binidxx,binidxy] += wgt^2
    else
        if (unsigned(binidxx - 1) < Lx) && (unsigned(binidxy - 1) < Ly)
            @inbounds h.hist.weights[binidxx,binidxy] += wgt
            @inbounds h.sumw2[binidxx,binidxy] += wgt^2
        end
    end
    return nothing
end

Base.broadcastable(h::Hist2D) = Ref(h)

"""
    Hist2D(elT::Type{T}=Float64; binedges, overflow) where {T}

Initialize an empty histogram with bin content typed as `T` and bin edges.
To be used with [`push!`](@ref). Default overflow behavior (`false`)
will exclude values that are outside of `binedges`.
"""
function Hist2D(elT::Type{T}=Float64; bins, overflow=_default_overflow) where {T}
    counts = zeros(elT, length.(bins) .- 1)
    return Hist2D(Histogram(bins, counts); overflow=overflow)
end

"""
    Hist2D(tuple, edges::NTuple{2,AbstractRange}; overflow)
    Hist2D(tuple, edges::NTuple{2,AbstractVector}; overflow)

Create a `Hist2D` with given bin `edges` and values from
a 2-tuple of arrays of x, y values. Weight for each value is assumed to be 1.
"""
function Hist2D(A::NTuple{2,AbstractVector}, r::NTuple{2,AbstractRange}; overflow=_default_overflow)
    h = Hist2D(Int; bins=r, overflow=overflow)
    push!.(h, A[1], A[2])
    return h
end
function Hist2D(A::NTuple{2,AbstractVector}, edges::NTuple{2,AbstractVector}; overflow=_default_overflow)
    if all(_is_uniform_bins.(edges))
        r = (range(first(edges[1]), last(edges[1]), length=length(edges[1])),
             range(first(edges[2]), last(edges[2]), length=length(edges[2])))
        return Hist2D(A, r; overflow=overflow)
    else
        h = Hist2D(Int; bins=edges, overflow=overflow)
        push!.(h, A[1], A[2])
        return h
    end
end

"""
    Hist2D(tuple, wgts::AbstractWeights, edges::NTuple{2,AbstractRange}; overflow)
    Hist2D(tuple, wgts::AbstractWeights, edges::NTuple{2,AbstractVector}; overflow)

Create a `Hist2D` with given bin `edges` and values from
a 2-tuple of arrays of x, y values.
`wgts` should have the same `size` as elements of `tuple`.
"""
function Hist2D(A::NTuple{2,AbstractVector}, wgts::AbstractWeights, r::NTuple{2,AbstractRange}; overflow=_default_overflow)
    @boundscheck @assert size(A[1]) == size(A[2]) == size(wgts)
    h = Hist2D(eltype(wgts); bins=r, overflow=overflow)
    push!.(h, A[1], A[2], wgts)
    return h
end
function Hist2D(A::NTuple{2,AbstractVector}, wgts::AbstractWeights, edges::NTuple{2,AbstractVector}; overflow=_default_overflow)
    if all(_is_uniform_bins.(edges))
        r = (range(first(edges[1]), last(edges[1]), length=length(edges[1])),
             range(first(edges[2]), last(edges[2]), length=length(edges[2])))
        return Hist2D(A, wgts, r; overflow=overflow)
    else
        h = Hist2D(Int; bins=edges, overflow=overflow)
        push!.(h, A[1], A[2], wgts)
        return h
    end
end

"""
    Hist2D(A::AbstractVector{T}; nbins::NTuple{2,Integer}, overflow) where T
    Hist2D(A::AbstractVector{T}, wgts::AbstractWeights; nbins::NTuple{2,Integer}, overflow) where T

Automatically determine number of bins based on `Sturges` algo.
"""
function Hist2D(A::NTuple{2,AbstractVector{T}};
        nbins::NTuple{2,Integer}=_sturges.(A),
        overflow=_default_overflow,
    ) where {T}
    F = float(T)
    nbinsx, nbinsy = nbins
    lox, hix = minimum(A[1]), maximum(A[1])
    loy, hiy = minimum(A[2]), maximum(A[2])
    rx = StatsBase.histrange(F(lox), F(hix), nbinsx)
    ry = StatsBase.histrange(F(loy), F(hiy), nbinsy)
    r = (rx, ry)
    return Hist2D(A, r; overflow=overflow)
end

function Hist2D(A::NTuple{2,AbstractVector{T}}, wgts::AbstractWeights;
        nbins::NTuple{2,Integer}=_sturges.(A),
        overflow=_default_overflow,
    ) where {T}
    F = float(T)
    nbinsx, nbinsy = nbins
    lox, hix = minimum(A[1]), maximum(A[1])
    loy, hiy = minimum(A[2]), maximum(A[2])
    rx = StatsBase.histrange(F(lox), F(hix), nbinsx)
    ry = StatsBase.histrange(F(loy), F(hiy), nbinsy)
    r = (rx, ry)
    return Hist2D(A, wgts, r; overflow=overflow)
end

"""
    function lookup(h::Hist2D, x, y)

For given x-axis and y-axis value `x`, `y`, find the corresponding bin and return the bin content.
If a value is out of the histogram range, return `missing`.
"""
function lookup(h::Hist2D, x, y)
    rx, ry = binedges(h)
    !(first(rx) <= x <= last(rx)) && return missing
    !(first(ry) <= y <= last(ry)) && return missing
    return bincounts(h)[_edge_binindex(rx, x), _edge_binindex(ry, y)]
end

"""
    normalize(h::Hist2D)

Create a normalized histogram via division by `integral(h)`.
"""
function normalize(h::Hist2D)
    return h*(1/integral(h))
end

"""
    rebin(h::Hist2D, nx::Int=1, ny::Int=nx)
    rebin(nx::Int, ny::Int) = h::Hist2D -> rebin(h, nx, ny)

Merges `nx` (`ny`) consecutive bins into one along the x (y) axis by summing.
"""
function rebin(h::Hist2D, nx::Int=1, ny::Int=nx)
    sx, sy = nbins(h)
    @assert sx % nx == sy % ny == 0
    p1d = (x,n)->Iterators.partition(x, n)
    p2d = x->(x[i:i+(nx-1),j:j+(ny-1)] for i=1:nx:sx, j=1:ny:sy)
    counts = sum.(p2d(bincounts(h)))
    sumw2 = sum.(p2d(h.sumw2))
    ex = first.(p1d(binedges(h)[1], nx))
    ey = first.(p1d(binedges(h)[2], ny))
    _is_uniform_bins(ex) && (ex = range(first(ex), last(ex), length=length(ex)))
    _is_uniform_bins(ey) && (ey = range(first(ey), last(ey), length=length(ey)))
    return Hist2D(Histogram((ex,ey), counts), sumw2; overflow=h.overflow)
end
rebin(nx::Int, ny::Int) = h::Hist2D -> rebin(h, nx, ny)

"""
    project(h::Hist2D, axis::Symbol=:x)
    project(axis::Symbol=:x) = h::Hist2D -> project(h, axis)

Computes the `:x` (`:y`) axis projection of the 2D histogram by
summing over the y (x) axis. Returns a `Hist1D`.
"""
function project(h::Hist2D, axis::Symbol=:x)
    @assert axis ∈ (:x, :y)
    dim = axis == :x ? 2 : 1
    ex, ey = binedges(h)
    counts = [sum(bincounts(h), dims=dim)...]
    sumw2 = [sum(h.sumw2, dims=dim)...]
    edges = axis == :x ? ex : ey
    return Hist1D(Histogram(edges, counts), sumw2; overflow=h.overflow)
end

"""
    transpose(h::Hist2D)

Reverses the x and y axes.
"""
function transpose(h::Hist2D)
    edges = reverse(binedges(h))
    counts = collect(bincounts(h)')
    sumw2 = collect(h.sumw2')
    return Hist2D(Histogram(edges, counts), sumw2; overflow=h.overflow)
end

"""
    profile(h::Hist2D, axis::Symbol=:x)
    profile(axis::Symbol=:x) = h::Hist2D -> profile(h, axis)

Returns the `axis`-profile of the 2D histogram by
calculating the weighted mean over the other axis.
`profile(h, :x)` will return a `Hist1D` with the y-axis edges of `h`.
"""
function profile(h::Hist2D, axis::Symbol=:x)
    @assert axis ∈ (:x, :y)
    if axis == :y
        h = transpose(h)
    end

    edges = binedges(h)[1]
    centers = bincenters(h)[2]
    counts = bincounts(h)
    sumw2 = h.sumw2

    num = counts*centers
    den = sum(counts, dims=2)
    numerr2 = sumw2 * centers.^2
    denerr2 = sum(sumw2, dims=2)
    val = vec(num ./ den)
    sw2 = vec(@. numerr2/den^2 - denerr2*(num/den^2)^2)

    # ROOT sets the NaN entries and their error to 0
    val[isnan.(val)] .= zero(eltype(val))
    sw2[isnan.(sw2)] .= zero(eltype(sw2))

    return Hist1D(Histogram(edges, val), sw2; overflow=h.overflow)
end
profile(axis::Symbol=:x) = h::Hist2D -> profile(h, axis)


text:
using ExtremeLearning

const nhidden = 100
const act_tanh = tanh
const act_sig = ExtremeLearning.sigmoid
const beta = 0.1
const train_length = 200
const x_train_dims = 4
const y_train_dims = 2
const predict_length = 2

inputs = ones(train_length, x_train_dims)
outputs = ones(train_length, y_train_dims)

elm = ELM(nhidden,inputs,outputs)
#testing constructor
@test isequal(nhidden, size(elm.weight_matrix)[2])
@test isequal(x_train_dims, size(elm.weight_matrix)[1])
@test isequal(nhidden, elm.n_hidden_neurons)
@test isequal(act_sig, elm.activation)
@test isequal(nhidden, size(elm.bias_vector)[2])
@test isequal(y_train_dims, size(elm.output_weights)[2])
@test isequal(nhidden, size(elm.output_weights)[1])
@test isequal(zero(inputs[1]), elm.β)

fit!(elm, inputs, outputs)
#testing training
@test isequal(y_train_dims, size(elm.output_weights)[2])
@test isequal(nhidden, size(elm.output_weights)[1])

new_inputs = ones(predict_length, x_train_dims)
prediction = predict(elm, new_inputs)
#testing a prediction
@test size(prediction)[2] == y_train_dims
@test size(prediction)[1] == predict_length

elm2 = ELM(nhidden,inputs,outputs, activation = act_tanh)

@test isequal(act_tanh, elm2.activation)

elm3 = ELM(nhidden,inputs, outputs, regularization = beta)

@test isequal(beta, elm3.β)


text:
"""
    runme(code::String = ",.", tapelen = 10000, a)

The interpreter of the HackMoji language.

# Arguments:
- `code`: Your HackMoji code given as a string
- `tapelen`: Your HackMoji memory given as number of bytes
"""
function runme(code::String=",.", tapelen=10000)
    data = fill(UInt8(0), tapelen)
    index = 1
    ptr = 1
    jumps = mapjumps(code)
    jumpstarts = collect(values(jumps))
    jumpends = collect(keys(jumps))
    while index <= length(code)
        op = code[index]
        if op=='.' print(Char(data[ptr]))
        elseif op==','
            inp=readline()
            if(length(inp)>0) data[ptr] = UInt8(inp[1])
            else data[ptr] = UInt8(0)
            end
        elseif op=='+' data[ptr] += UInt8(1)
        elseif op=='-' data[ptr] -= UInt8(1)
        elseif op=='>' ptr += 1
        elseif op=='<' ptr -= 1
        elseif op=='['
            if data[ptr] == 0x00
                index = jumpends[findfirst(isequal(index), jumpstarts)]
            end
        elseif op==']'
            if data[ptr] != 0x00
                index = jumps[index]
            end
        end
        index += 1
    end
end

"""
    mapjumps(code=",.[,.,.+++[,..],--]..,")

Register the beginning and end of all loops in the HackMoji Code.
"""
function mapjumps(code::String=",.[,.,.+++[,..],--]..,")
    jumps = Dict{Int, Int}()
    stack = []
    for i in 1:length(code)
        if code[i] == '['  push!(stack, i)
        elseif code[i] == ']' jumps[i] = pop!(stack)
        end
    end
    return jumps
end


text:
"""
	MM1(μ, λ)

Tạo mô hình M/M/1
"""
struct MM1{T<:Any} <: AbstractMMCK
	λ::Union{T, Real}
	μ::Union{T, Real}
	ρ::Union{T, Real}
	function MM1(λ, μ)
		T = Union{typeof(μ), typeof(λ)}
		new{T}(λ, μ, λ/μ)
	end
	function MM1(λ::Real, μ::Real)
		if λ/μ >= 1
			@error "ρ phải < 1"
		end
		T = Union{typeof(μ), typeof(λ)}
		new{T}(λ, μ, λ/μ)
	end
end

pn(m::MM1, n::Int) = (1 - m.ρ) * m.ρ^n
L(m::MM1) = m.ρ / (1 - m.ρ)
Lq(m::MM1) = m.ρ ^ 2 / (1 - m.ρ)
W(m::MM1) = L(m) / m.λ
Wq(m::MM1) = Lq(m) / m.λ


text:
# function (load::CSVImporter{CSVSchema{Empatica, E4, SkinConductance}})(path::Path)
#     @show "Inside Empatica Skin Conductance importer..."
#     disable!(load)
# end

function create_import_context(vendor::AbstractVendor, product::AbstractProduct, data::AbstractData, keystr::String)
    # The controller is responsible for instantiating the dialog GUI,
    # handling user input and tiggering the file importer once the user
    # confirms the selection of a particular file.
    control = FileDialogControl(true)
    # The model store the path that the user has confirmed, as well as the
    # path that the user has selected but not necessarily confirmed.
    model = FileDialogModel(Path(pwd(),""), Path(pwd(),""))
    # The display properties include the title of the file dialog and what
    # action string is to be displayed on the button.
    properties = FileDialogDisplayProperties(caption = "Open File###"*keystr, action = "Open###"*keystr, width = Cfloat(640), height = Cfloat(395))
    # The importer is triggered once the user has confirmed the selection of a
    # particular path.
    load_csv = CSVImporter(false, CSVSchema(vendor, product, data))
    # Construct a context and assign roles to complete the requisite import interactions.
    ImportContext(control, model, properties, load_csv)
end

function create_import_context(vendor::Empatica, product::E4, data::SkinConductance, keystr::String)
    # The controller is responsible for instantiating the dialog GUI,
    # handling user input and tiggering the file importer once the user
    # confirms the selection of a particular file.
    control = FileDialogControl(true)
    # The model store the path that the user has confirmed, as well as the
    # path that the user has selected but not necessarily confirmed.
    model = FileDialogModel(Path(pwd(),""), Path(pwd(),""))
    # The display properties include the title of the file dialog and what
    # action string is to be displayed on the button.
    properties = FileDialogDisplayProperties(caption = "Open File###"*keystr, action = "Open###"*keystr, width = Cfloat(640), height = Cfloat(395))
    # The importer is triggered once the user has confirmed the selection of a
    # particular path.
    load_csv = CSVImporter(false, CSVSchema(vendor, product, data))
    # Construct a context and assign roles to complete the requisite import interactions.
    ImportContext(control, model, properties, load_csv)
end

function (load::CSVImporter{<:CSVSchema{<: Empatica, <: E4, <: SkinConductance}})(path::Path)
    disable!(load)
    df = load_dataframe(path, ["EDA"])
    if isnothing(df)
        CImGui.OpenPopup("Have you opened the appropriate file?")
        return nothing
    else
        data = df[:EDA]
        timestamp = data[1]
        sampling_frequency = convert(Int,data[2])
        return ElectrodermalData(timestamp , sampling_frequency, data[3:end], data[3:end])
    end
 end

 function (load::CSVImporter{<:CSVSchema{<: Empatica, <: E4, <: Tags}})(path::Path)
     disable!(load)
     df = load_dataframe(path, ["Tags"])
     if isnothing(df)
         CImGui.OpenPopup("Have you opened the appropriate file?")
         return nothing
     else
         data = df[:,1]
         return data
     end
  end


text:
# ---
# title: 1598. Crawler Log Folder
# id: problem1598
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Stack
# link: <https://leetcode.com/problems/crawler-log-folder/description/>
# hidden: true
# ---
#
# The Leetcode file system keeps a log each time some user performs a _change
# folder_ operation.
#
# The operations are described below:
#
#   * `"../"` : Move to the parent folder of the current folder. (If you are already in the main folder, **remain in the same folder** ).
#   * `"./"` : Remain in the same folder.
#   * `"x/"` : Move to the child folder named `x` (This folder is **guaranteed to always exist** ).
#
# You are given a list of strings `logs` where `logs[i]` is the operation
# performed by the user at the `ith` step.
#
# The file system starts in the main folder, then the operations in `logs` are
# performed.
#
# Return _the minimum number of operations needed to go back to the main folder
# after the change folder operations._
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/09/sample_11_1957.png)
#
#
#
#     Input: logs = ["d1/","d2/","../","d21/","./"]
#     Output: 2
#     Explanation: Use this change folder operation "../" 2 times and go back to the main folder.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/09/sample_22_1957.png)
#
#
#
#     Input: logs = ["d1/","d2/","./","d3/","../","d31/"]
#     Output: 3
#
#
# **Example 3:**
#
#
#
#     Input: logs = ["d1/","../","../","../"]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= logs.length <= 103`
#   * `2 <= logs[i].length <= 10`
#   * `logs[i]` contains lowercase English letters, digits, `'.'`, and `'/'`.
#   * `logs[i]` follows the format described in the statement.
#   * Folder names consist of lowercase English letters and digits.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#function [G1,CC,impact,fmat,fwt,ywt,gev,eu]=

    function gensysToAMA(g0, g1, cc, psi, pi, div, varargin = "" )
        # function [G1,CC,impact,fmat,fwt,ywt,gev,eu]=gensysToAMA(g0,g1,cc,psi,pi,div,
        # optionalArg)
        # gensys interface to both gensys and the Anderson-Moore algorithm.
        # Just as with gensys, system given as
        #        g0*y(t)=g1*y(t-1)+c+psi*z(t)+pi*eta(t),
        # with z an exogenous variable process and eta being endogenously determined
        # one-step-ahead expectational errors.  Returned system is
        #       y(t)=G1*y(t-1)+C+impact*z(t)+ywt*inv(I-fmat*inv(L))*fwt*z(t+1) .
        # If z(t) is i.i.d., the last term drops out.
        # If div is omitted from argument list, a div>1 is calculated.
        # eu(1)=1 for existence, eu(2)=1 for uniqueness.  eu(1)=-1 for
        # existence only with not-s.c. z; eu=[-2,-2] for coincident zeros.
        # By Christopher A. Sims
        # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
        # when called with no optional args,
        #        program first tries to find gensys on the matlab path
        #        if that fails, the program runs gensys2007 in the gensysToAMA directory
        #
        # an optional string argument may follow the original gensys arguments
        # 'gensys'   run gensys
        #        program first tries to find gensys on the matlab path
        #        if that fails, the program runs gensys2007 in the gensysToAMA directory
        #
        # 'gensys2007'   run gensys
        #        the program runs gensys2007 in the gensysToAMA directory
        #
        # 'ama'    run the anderson-moore algorithm with gensys inputs and outputs
        #
        # 'both'    run the anderson-moore algorithm and the gensys program
        #        verify that solutions are equivalent print out execution times
        #

        reps = 5 #for timing tests average

        G1     = nothing #'not set';
        CC     = nothing #'not set';
        impact = nothing #'not set';
        fmat   = nothing #'not set';
        fwt    = nothing #'not set';
        ywt    = nothing #'not set';
        gev    = nothing #'not set';
        eu     = nothing #'not set';
        div    = 1.0

        if varargin == ""  #assume calling as though original gensys

            try
                println("gensysToAMA: trying gensys on your matlab path")
                (G1, CC, impact, fmat, fwt, ywt, gev, eu) = gensys(g0, g1, cc, psi, pi, div)
                sss = run(`which gensys`) # which('gensys')
                println("used gensys located at:  $sss")
            catch
                println("gensysToAMA:  that failed, using gensys2007 in gensysToAMA dir")
                (G1, CC, impact, fmat, fwt, ywt, gev, eu) = gensys2007(g0, g1, cc, psi, pi, div)
                sss = run(`which gensys`) # which('gensys2007')
                println("used gensys2007 located at:  $sss")
            end

	elseif varargin == "ama"

            println("gensysToAMA:running ama")
            println("gensysToAMA:converting gensys inputs to ama format")
            (theHM, theH0, theHP) = convertFromGensysIn(g0, g1, pi)
            condn  = 1.e-10     #AndersonMooreAlg uses this in zero tests
            uprbnd = 1 + 1.e-6  #allow unit roots
            # theH = [theHM, theH0, theHP]
            theH = hcat(theHM, theH0, theHP)
            neq = size(theHM, 1)
            println("gensysToAMA:running ama")
            (bb, rts, ia, nexact, nnumeric, lgroots, aimcode) =
                AndersonMooreAlg(theH, neq, 1, 1, condn, uprbnd)
            eu = setEu(aimcode)

            if aimcode == 1

                #    scof = SPObstruct(theH,bb,neq,1,1);not needed for generating sims output
                phi  = inv(theH0 + theHP * bb) #inv
                theF = -phi * theHP
                ncpi = size(pi, 2)
                println("gensysToAMA:converting ama output to gensys format")
                (CC, G1, impact, ywt, fmat, fwt) =
                    convertToGensysOut(bb, phi, theF, cc, g0, g1, psi, ncpi)
            else
                println("no unique solution: not bothering to try and convert to gensys output")
            end

            println("gensysToAMA:done")

#####
#####  Implement below
#####
#============================================================
        elseif varargin == "gensys" #assume calling as though original gensys but with a trailing optional arg
            #arg1=varargin{1}
            #switch arg1
            #case 'gensys'

            println("gensysToAMA:running gensys")

            try
                println("gensysToAMA: trying gensys on your matlab path")
                (G1, CC, impact, fmat, fwt, ywt, gev, eu) = gensys(g0, g1, cc, psi, pi, div)
                sss = run(`which gensys`)  # which('gensys')
                println("used gensys located at:  $sss")
            catch
                println("gensysToAMA:  that failed, using gensys2007 in gensysToAMA dir")
                (G1, CC, impact, fmat, fwt, ywt, gev, eu) = gensys2007(g0, g1, cc, psi, pi, div)
                sss = run(`which gensys`) # which('gensys2007')
                println("used gensys2007 located at:  $sss")
                println("gensysToAMA:done")
            end

        elseif varargin == "gensys2007" # case 'gensys2007'

            println("gensysToAMA:running gensys2007")
            (G1, CC, impact, fmat, fwt, ywt, gev, eu) = gensys2007(g0, g1, cc, psi, pi, div)
            sss = run(`which gensys2007`) # which('gensys2007')
            println("used gensys2007 located at:  $sss")
            println("gensysToAMA:done")

        elseif varargin == "both"
            println("gensysToAMA:running both ama and gensys for comparison")
            parg1 = size(g0,1)
            parg2 = size(g0,2)
            parg3 = size(psi,1)
            parg4 = size(psi,2)
            parg5 = size(pi,1)
            parg6 = size(pi,2)
            println("problem dimensions: g0:$parg1 x $parg2, psi:$parg3 x $parg4, pi:$parg5 x $parg6")
            println("gensysToAMA: doing  $reps reps for timing averages")
            println("gensysToAMA:running gensys")

            try
                println("gensysToAMA: trying gensys on your matlab path")
                (sG1, sCC, simpact, sfmat, sfwt, sywt, sgev, eu) = gensys(g0, g1, cc, psi, pi, div)
                sss = run(`which gensys`) # which('gensys')
                println("used gensys located at:  $sss")
            catch
                println("gensysToAMA: that failed, using gensys2007 in gensysToAMA dir")
                (sG1, sCC, simpact, sfmat, sfwt, sywt, sgev, eu) = gensys2007(g0, g1, cc, psi, pi, div)
                sss = run(`which gensys2007`) # which('gensys2007')
                println("used gensys2007 located at:  $sss")
            end

            gensysDone=0
            for ii = 1:reps

                gensysStart=cputime
                try
                    (sG1, sCC, simpact, sfmat, sfwt, sywt, sgev, eu) = gensys(g0, g1, cc, psi, pi, div)
                catch
                    (sG1, sCC, simpact, sfmat, sfwt, sywt, sgev, eu) = gensys2007(g0, g1, cc, psi, pi, div)
                end

                gensysDone = gensysDone + (cputime - gensysStart)
            end

            gensysDone = gensysDone / reps
            println("gensysToAMA:running ama")
            convertToStart = cputime
            (theHM, theH0, theHP) = convertFromGensysIn(g0, g1, pi)
            theH = (theHM, theH0, theHP)
            condn  = 1.e-10 #AndersonMooreAlg uses this in zero tests
            uprbnd = 1 + 1.e-6 #allow unit roots
            neq = size(theHM, 1)
            convertToDone =cputime - convertToStart
            (bb, rts, ia, nexact, nnumeric, lgroots, aimcode) =
                SPAmalg(theH,neq,1,1,condn,uprbnd);
            AMADone=0

            for ii = 1:reps
                AMAStart=cputime
                (bb, rts, ia, nexact, nnumeric, lgroots, aimcode) =
                    SPAmalg(theH,neq,1,1,condn,uprbnd)
                #    scof = SPObstruct(theH,bb,neq,1,1);#not needed for generating sims output
                #phi=inv(theH0+theHP*sparse(bb));
                AMADone = AMADone + (cputime - AMAStart)
            end

            AMADone = AMADone / reps

            if aimcode == 1
                #NAStart=cputime;
                #[qq,bb,info,sinv]=numericAim(1,theH);
                #NADone=cputime-NAStart;

                AMAFStart = cputime
                theF = -((theH0 + theHP * sparse(bb)) \ theHP)
                AMAFDone = cputime - AMAFStart
                println("gensysToAMA:converting ama output to gensys format");
                AMAConvertStart = cputime
                x=1
                eu = setEu(x)
                phi = inv(theH0 + theHP * sparse(bb))
                #note only need column dimension of pi for conversion
                ncpi =size(pi, 2)
                (CC, G1, impact, ywt, fmat, fwt) =
                    convertToGensysOut(bb,phi,theF,cc,g0,g1,psi,ncpi)
                AMAConvertDone = cputime - AMAConvertStart
            else
                println("no unique solution: not bothering to try and convert to gensys output")
            end
            println("gensysToAMA:runs complete")
            #[newvv,newfmat,newNz]=smallF(sfmat);

            if aimcode == 1
                try
                    theDiffs=
                        [norm(real(sG1) -G1), norm(real(sCC) - CC), norm(real(simpact) - impact),
                         simMats(sfmat,fmat),
                         norm(real(sywt * sfwt) -ywt * fwt), norm(real(sywt * sfmat * sfwt) -ywt * fmat * fwt)]
                catch
                    println("detected a difference in sims and AMA results")
                    #println(sprintf('dims of matrices, G1''s , C''s. impact''s, fmat''s, ywt''s, fwt''s'))

                    [size(sG1) size(G1);size(sCC) size(CC);
                     size(simpact) size(impact);size(sfmat) size(fmat);
                     size(sywt) size(ywt);size(sfwt) size(fwt)]
                end

                if(any(theDiffs > 1e-7))
                    println("detected a difference in sims and AMA results")
                    theDiffs
                else
                    println("no difference in sims and AMA results")
                end

                println("gensysToAMA: after one primer run and   $reps reps")
                print("AMATime=$AMADone  AMAFTime=$AMAFDone  \n convertToTime=$convertToDone convertFromTime=$AMAConvertDone  \n genSysTime=$gensysDone \n")

            else
                println("did not compute conversion, no comparisons to gensys output")
                println("AMATime=$AMADone   genSysTime=$gensysDone")
            end
#####
##### Implement above
#####
==============================================================#

	else
	    println("gensysToAMA:unknown optional string")
	end

        return G1,CC,impact,fmat,fwt,ywt,gev,eu
    end


text:
module Phases

"""
    abstract type Phase

Abstract supertype for all phases. See `subtypes(FluxTraining.Phase)`.
A `Phase` is used in dispatch for training loop functions [`step!`](#)
and [`epoch!`](#) as well as in [`Callback`](#) handler methods [`on`](#).
"""
abstract type Phase end

"""
    abstract type AbstractTrainingPhase <: Phase

An abstract type for phases where parameter updates
are being made. This exists so callbacks can dispatch on it and work
with custom training phases.

The default implementation for supervised tasks is [`TrainingPhase`](#).
"""
abstract type AbstractTrainingPhase <: Phase end

"""
    TrainingPhase() <: AbstractTrainingPhase

A regular training phase for supervised learning. It iterates over batches
in `learner.data.training` and updates the model parameters
using `learner.optim` after calculating the gradients.

Throws the following events in this order:

- [`EpochBegin`](#) when an epoch starts,
- [`StepBegin`](#) when a step starts,
- [`LossBegin`](#) after the forward pass but before loss calculation,
- [`BackwardBegin`](#) after loss calculation but before backward pass,
- [`BackwardEnd`](#) after the bacward pass but before the optimization
    step,
- [`StepEnd`](#) when a step ends; and
- [`EpochEnd`](#) when an epoch ends

It writes the following step state to `learner.state`, grouped by the
event from which on it is available.

- `StepBegin`:
    - `xs` and `ys`: encoded input and target (batch)
- `LossBegin`:
    - `ŷs`: model output
- `BackwardBegin`:
    - `loss`: loss
- `BackwardEnd`:
    - `grads`: calculated gradients
"""
struct TrainingPhase <: AbstractTrainingPhase end

"""
    abstract type AbstractValidationPhase <: Phase

An abstract type for phases where no parameter updates
are being made. This exists so callbacks can dispatch on it and work
with custom validation phases.

The default implementation for supervised tasks is [`ValidationPhase`](#).
"""
abstract type AbstractValidationPhase <: Phase end

"""
    ValidationPhase()

A regular validation phase. It iterates over batches
in `learner.data.validation` and performs a forward pass.

Throws the following events: [`EpochBegin`](#), [`StepBegin`](#),
[`LossBegin`](#), [`StepEnd`](#), [`EpochEnd`](#).

Throws the following events in this order:

- [`EpochBegin`](#) when an epoch starts,
- [`StepBegin`](#) when a step starts,
- [`LossBegin`](#) after the forward pass but before loss calculation,
- [`StepEnd`](#) when a step ends; and
- [`EpochEnd`](#) when an epoch ends

It writes the following step state to `learner.state`, grouped by the
event from which on it is available.

- `StepBegin`:
    - `xs` and `ys`: encoded input and target (batch)
- `LossBegin`:
    - `ŷs`: model output
- `StepEnd`:
    - `loss`: loss
"""
struct ValidationPhase <: AbstractValidationPhase end

export
    Phase,
    AbstractTrainingPhase,
    TrainingPhase,
    AbstractValidationPhase,
    ValidationPhase
end # module

using .Phases


text:
using Skans
using Skans:
    TRIPLEQUOTE
using Gumbo: parsehtml
using Test
using TOML:
    parse as tomlparse

notify = false

const PAGES = [
    Skans.MockPage("url1", "a"),
    Skans.MockPage("url2", "b")
]

function pages2state(pages::Vector{<:Skans.Page})
    scans = Skans.scan.(pages)
    return Skans.State(scans)
end

include("select.jl")

@testset "skan! updating" begin
    state = pages2state(PAGES)

    repo = Skans.MockRepo(state)

    changed = skan!(repo, PAGES; notify)
    @test isempty(changed)

    page = Skans.MockPage("url1", "c")
    pagescan = Skans.scan(page)
    changed = skan!(repo, [page, PAGES[2]]; notify)
    @test length(changed) == 1
    @test first(changed).content == "c"
    state = Skans.retrieve(repo)
    @test state.scans["url1"].content == "c"
end

@testset "HTTP get" begin
    url = "http://example.com"
    page = WebPage(url)

    repo = Skans.MockRepo()

    changed = skan!(repo, [page]; notify)
    @test !isempty(changed)
end

@testset "Clone and pull GitHub repo" begin
    repo = GitHubRepo(;
        user="",
        token="",
        repo="githubtraining/hellogitworld",
        branch="master"
    )
    Skans.pull_or_clone!(repo)
    @test !isempty(readdir(repo.dir))
    # Second time goes through pull logic.
    Skans.pull_or_clone!(repo)
end

@testset "diff" begin
    @test Skans.startswith_one("+", '+') == true
    @test Skans.startswith_one("+-", '+') == true
    @test Skans.startswith_one("++", '+') == false
end

@testset "Store and read MockFileRepo" begin
    repo = Skans.MockFileRepo()
    changed = skan!(repo, PAGES; notify)
    @test length(changed) == 2
    @test first(changed).content == "a"

    changed = skan!(repo, PAGES; notify)
    @test isempty(changed)
end

@testset "Multiline TOML" begin
    state = pages2state(PAGES)
    expected = """
        "url1" = $TRIPLEQUOTE
        a$TRIPLEQUOTE
        "url2" = $TRIPLEQUOTE
        b$TRIPLEQUOTE"""
    actual = Skans.toml(state.scans)
    @test expected == actual

    trouble = raw"userAgent.match(/Firefox[\/\s](\d+\.\d+)/)"
    page = Skans.MockPage("u", trouble)
    state = pages2state([page])
    # Smoke test.
    tomlparse(Skans.toml(state.scans))

    trouble = raw"(/[^a-z0-9\-#$%&'*+.\^_`|~]"
    page = Skans.MockPage("u", trouble)
    state = pages2state([page])
    # Smoke test.
    tomlparse(Skans.toml(state.scans))
end


text:
####################
# Continuous Lagrange
####################

struct ContinuousLagrange{dim,shape,qorder,gorder} <: FiniteElement{dim,shape,qorder,gorder}
    nodal_base_coefs::Matrix{Float64}
    geom_base_coefs::Matrix{Float64}
    topology::Dict{Int,Int}
    geom_topology::Dict{Int,Int}
    M::Matrix{Float64}
end

function ContinuousLagrange{dim,shape,qorder}() where {dim,shape<:AbstractRefShape, qorder}
    return ContinuousLagrange{dim, shape, qorder,1}()
end

function ContinuousLagrange{dim,shape,qorder,gorder}() where {dim,shape<:AbstractRefShape, qorder, gorder}
    nodal_base_coefs, topology = _nodal_data(dim, shape, qorder)
    if qorder == gorder
        geom_base_coefs = nodal_base_coefs
        geom_topology = topology
    else
        geom_base_coefs, geom_topology = _nodal_data(dim, shape, gorder)
    end
    if gorder == 1
        M = one(geom_base_coefs)
    else
        M = _nodal_geom_data(dim, shape, gorder, size(geom_base_coefs,1), geom_base_coefs)
    end
    ContinuousLagrange{dim, shape, qorder, gorder}(nodal_base_coefs, geom_base_coefs, topology, geom_topology, M)
end

function _nodal_data(dim::Int,shape::Type{Shape}, order::Int) where {Shape <: AbstractRefShape}
    nodal_points, topology = get_nodal_points(shape, Val{dim}, order)
    ip_prime = getdefaultdualbasis(dim,shape,order)
    nbasefuncs = getnbasefunctions(ip_prime)
    V = [value(ip_prime, j, nodal_points[i]) for i=1:nbasefuncs,j = 1:nbasefuncs]
    return inv(V), topology
end

function _nodal_geom_data(dim::Int,shape::Type{Shape}, order::Int, n_geom_basefuncs::Int, geom_base_coefs::Matrix) where {Shape <: AbstractRefShape}
    # Matrix to get spacial coordinates
    nodal_points, topology = get_nodal_points(shape, Val{dim}, order)
    T = eltype(nodal_points[1])
    qrs = QuadratureRule{dim,shape,T}(fill(T(NaN), length(nodal_points)), nodal_points) # weights will not be used
    n_qpoints = length(getweights(qrs))
    M =    fill(zero(T)           * T(NaN), n_geom_basefuncs, n_qpoints)
    for (qp, ξ) in enumerate(qrs.points)
        for i in 1:n_geom_basefuncs
            M[i, qp] = dot(geom_base_coefs[:,i], value(defaultDualBasis(dim,shape,order), ξ))
        end
    end
    M
end

@inline getnbasefunctions(::ContinuousLagrange{1,RefTetrahedron,order}) where {order} = order + 1
@inline getnbasefunctions(::ContinuousLagrange{2,RefTetrahedron,order}) where {order} = Int((order+1)*(order+2)/2)
@inline getngeombasefunctions(::ContinuousLagrange{1,RefTetrahedron,order,gorder}) where {order,gorder} = gorder + 1
@inline getngeombasefunctions(::ContinuousLagrange{2,RefTetrahedron,order,gorder}) where {order,gorder} = Int((gorder+1)*(gorder+2)/2)
@inline gettopology(ip::ContinuousLagrange) = ip.topology
@inline getgeomtopology(ip::ContinuousLagrange) = ip.geom_topology

"""
value(ip::ContinuousLagrange{dim,shape<:AbstractRefShape,order}, k::Int, ξ::Vec{dim,T}) where {dim,shape,order, T}
Compute value of Continuous Lagrange Finite Element basis `j` at point ξ on shape `shape`
"""
function value(ip::ContinuousLagrange{dim,shape,order}, k::Int, ξ::Vec{dim,T}) where {dim,shape<:AbstractRefShape,order, T}
    if k > getnbasefunctions(ip);throw(ArgumentError("no shape function $k for interpolation $ip"));end
    n = getnbasefunctions(ip)
    dot(ip.nodal_base_coefs[:,k], value(defaultDualBasis(dim,shape,order), ξ))
end

"""
geom_value(ip::ContinuousLagrange{dim,shape<:AbstractRefShape,order}, k::Int, ξ::Vec{dim,T}) where {dim,shape,order, T}
Compute value of Continuous Lagrange Finite Element basis `j` at point ξ on shape `shape`
"""
function geom_value(ip::ContinuousLagrange{dim,shape,order,gorder}, k::Int, ξ::Vec{dim,T}) where {dim,shape<:AbstractRefShape,order,gorder, T}
    if k > getnbasefunctions(ip);throw(ArgumentError("no shape function $k for interpolation $ip"));end
    n = getnbasefunctions(ip)
    dot(ip.geom_base_coefs[:,k], value(defaultDualBasis(dim,shape,gorder), ξ))
end


text:
#ngb treewalk: recursively find gather ngbs for particle p within the searching radius
function get_ptl_in_box_gather(radius::T, p::SVector{N,T}, node::Node{N,T,D}, boxsizes::SVector{N,T}, idx_ngbs::Vector{Int64}) where {N,T,D<:AbstractData{N,T}}
    if isLeaf(node)
        #println("in a leaf node")
        if node.p != nothing
            #println("nonempty leaf")
	    	dist2 = get_distance2(node.p.pos, p, boxsizes, true)
            if dist2 < radius^2
                #println("push ", node.part)
				push!(idx_ngbs, node.p.idx)
            end
        end
    else
        #println("This is a node... check if node and radius overlap")
		@inbounds for j in 1:N
			if abs(nearest(node.center[j] - p[j], boxsizes[j])) > 0.5*node.length[j] + radius
				return  #no overlap, so node can be skipped!
			end
		end
		#=
		if get_distance2(node.center, p, boxsizes, true) + 0.75 * node.length[1]^2 < radius^2
			append!(idx_ngbs, node.n.idxs)
			return
		end
		=#
		#println("node and radius overlap, so we need to open the node and check its children")
        @inbounds for i in 1:2^N
			get_ptl_in_box_gather(radius, p, node.child[i], boxsizes, idx_ngbs)
        end
    end
end

#ngb treewalk: recursively find scatter ngbs for particle p
function get_ptl_in_box_scatter(p::SVector{N,T}, node::Node{N,T,D}, boxsizes::SVector{N,T}, idx_ngbs::Vector{Int64}) where {N,T,D<:AbstractData{N,T}}
    if isLeaf(node)
        #println("in a leaf node")
		if node.p != nothing
            #println("nonempty leaf")
	    	dist2 = get_distance2(node.p.pos, p, boxsizes, true)
            if dist2 < node.p.hsml^2
                #println("push ", node.part)
				push!(idx_ngbs, node.p.idx)
            end
        end
    else
		#println("This is a node... check if node and radius overlap")
		@inbounds for j in 1:N
			if abs(nearest(node.center[j] - p[j], boxsizes[j])) > 0.5*node.length[j] + node.n.hsml
				return  #no overlap, so node can be skipped!
			end
		end
		#println("node and radius overlap, so we need to open the node and check its children")
        @inbounds for i in 1:2^N
			get_ptl_in_box_scatter(p, node.child[i], boxsizes, idx_ngbs)
        end
    end
end

#driver function to get gather ngbs for particle x within the searching radius h
function get_gather_ngb_tree(x::SVector{N,T}, h::T, node::Node{N,T,D}, boxsizes::SVector{N,T}) where {N,T,D<:AbstractData{N,T}}
    idx_ngbs = Int64[]
    get_ptl_in_box_gather(h, x, node, boxsizes, idx_ngbs)
    return idx_ngbs
end

#driver function to get scatter ngbs for particle x; no need to specify a searching radius
function get_scatter_ngb_tree(x::SVector{N,T}, node::Node{N,T,D}, boxsizes::SVector{N,T}) where {N,T,D<:AbstractData{N,T}}
    idx_ngbs = Int64[]
    get_ptl_in_box_scatter(x, node, boxsizes, idx_ngbs)
    return idx_ngbs
end


text:
conf_file = joinpath(@__DIR__,"..")
readdir(conf_file)


text:
workspace()
using JukaiNLP
using JukaiNLP.Tagging
using Merlin

train(10, traindata, testdata)

modelpath = "C:/Users/hshindo/Desktop/postagger.h5"
h5save(modelpath, Dict("wordfun"=>t.model.wordfun,"charfun"=>t.model.charfun,"sentfun"=>t.model.sentfun))
t.model = nothing
path = "C:/Users/hshindo/Desktop/postagger.jld"
JLD.save(path, "postagger", t)
throw("finish")

modelpath = "C:/Users/hshindo/Desktop/postagger.h5"
m = h5load(modelpath)
model = Tagging.POSModel(m["wordfun"],m["charfun"],m["sentfun"])
path = "C:/Users/hshindo/Desktop/postagger.jld"
t = JLD.load(path, "postagger")
t.model = model
t(["I","have","a","pen","."]) |> println

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
"""
A node that represents a piece of a JavaScript abstract syntax tree.

Every node is either a [`JSAST`](@ref) (_i.e._, non-terminal node) or a
[`JSTerminal`](@ref).
"""
abstract type JSNode end

"""
    JSTerminal(s::JSString)

A JavaScript terminal node.
Any JavaScript literal is represented as a `JSTerminal` node.
"""
struct JSTerminal <: JSNode
    s::JSString
end

JSTerminal(x::Union{Symbol, Number, Bool}) = JSTerminal(JSString(string(x)))
JSTerminal(x::String) = JSTerminal(js"$x")
JSTerminal(::Nothing) = JSTerminal(js"null")

"""
    JSAST(head::Symbol, args::JSNode...)

A struct that represents an abstract syntax tree for a JavaScript expression.
This is meant to be analogous to Julia's `Expr`s.
"""
struct JSAST <: JSNode
    head::Symbol
    args::Vector{JSNode}
end

JSAST(h) = JSAST(h, [])
JSAST(h, body::Union{JSNode, Nothing}...) = JSAST(h, [b for b in body if nothing !== b])

# Custom show methods to make printing look less noisy.
Base.show(io:: IO, m::MIME"text/plain", t::JSTerminal) = print(io, "JSTerminal($(t.s.s))")
function Base.show(io:: IO, m::MIME"text/plain", t::JSAST)
    print(io, "JSAST(", QuoteNode(t.head))
    foreach(t.args) do x
        print(io, ", ")
        Base.show(io, m, x)
    end
    print(io, ")")
end

_unquote(x) = x
_unquote(x::QuoteNode) = _unquote(x.value)


text:
export oneplus

using Distributed

function oneplus(nin::Int64, nout::Int64, fitness::Function;
                 ctype::DataType=CGPChromo, seed::Int64=0, expert::Any=nothing,
                 id::String="")
    population = Array{ctype}(undef, Config.lambda)
    for i in eachindex(population)
        population[i] = ctype(nin, nout)
    end
    if expert != nothing
        population[1] = expert
    end
    best = population[1]
    max_fit = -Inf
    eval_count = 0
    fits = -Inf*ones(Float64, Config.lambda)

    while eval_count < Config.total_evals
        # evaluation
        log_gen = false

        new_fits = pmap(fitness, population)

        for p in eachindex(population)
            # fit = fitness(population[p])
            fit = new_fits[p]
            eval_count += 1
            if fit >= max_fit
                best = clone(population[p])
                if fit > max_fit
                    max_fit = fit
                    log_gen = true
                end
            end
            fits[p] = fit
            if eval_count == Config.total_evals
                log_gen = true
                break
            end
        end

        eval(Config.log_function)(id, seed, eval_count, max_fit, best, GA,
                                  ctype, log_gen)

        if eval_count == Config.total_evals
            break
        end

        # selection
        fits .= -Inf
        for p in eachindex(population)
            population[p] = mutate(best)
        end

        # size limit
        for i in eachindex(population)
            if length(population[i].nodes) > Config.node_size_cap
                @info("Discarding individual because above node size cap")
                population[i] = ctype(nin, nout)
                fits[i] = -Inf
            end
        end
    end

    max_fit, best.genes
end


text:
#!/usr/bin/env julia
#
using Printf

println(" Running with Julia ", VERSION)

macro VarName(arg)
    string(arg)
end

function printsum(varname, a)
    println(" ", varname, " => ", summary(a), ": ", repr(a))
end

function printline(varname, a)
    println(" ", varname, " => ", a)
end

# dictionaries can be initialised directly:
a1 = Dict(1=>"one", 2=>"two")
printsum(@VarName(a1), a1) #> Dict{Int64,String}: {2=>"two",1=>"one"}

# then added to:
a1[3]="three"
printsum(@VarName(a1), a1) #> Dict{Int64,String}: {2=>"two",3=>"three",1=>"one"}
# (note dictionaries cannot be assumed to keep their original order)

# dictionaries may also be created with the type explicitly set
a2 = Dict{Int64, AbstractString}()
a2[0]="zero"
printsum(@VarName(a2), a2)
#> Dict{Int64,AbstractString} with 1 entry: Dict{Int64,AbstractString}(Pair{Int64,AbstractString}(0,"zero"))

# dictionaries like arrays, may also be created from comprehensions
a3 = Dict([i => @sprintf("%d", i) for i = 1:10])
printsum(@VarName(a3), a3)
#> Dict{Any,Any}: {5=>"5",4=>"4",6=>"6",7=>"7",2=>"2",10=>"10",9=>"9",8=>"8",3=>"3",1=>"1"}

# as you would expect, Julia comes with all the normal helper functions
# for dictionaries
printline(@VarName(a1), haskey(a1,1)) #> true

# which is equivalent to
printline(@VarName(a1), 1 in keys(a1)) #> true
# where keys creates an iterator over the keys of the dictionary

# similar to keys, values get iterators over the dict's values:
printsum(@VarName(a1), values(a1))
#> Base.ValueIterator for a Dict{Int64,String} with 3 entries: String["two","three","one"]

# use collect to get an array:
printsum(@VarName(a1), collect(values(a1)))
#> 3-element Array{String,1}: String["two","three","one"]

max_stars_repo_name: wbeebe/julia

max_issues_repo_name: wbeebe/julia

max_forks_repo_name: wbeebe/julia


text:
############################################################################################
#                                  EXPLICIT METHODS 								       #
############################################################################################

function integrate(::Type{Euler}, model::AbstractModel, x::StaticVector, u::StaticVector, t, dt)
	xdot = dynamics(model, x, u, t)
	return x + xdot * dt
end

function integrate(::Type{RK2}, model::AbstractModel, x::StaticVector, u::StaticVector, t, dt)
	k1 = dynamics(model, x,        u, t       )*dt
	k2 = dynamics(model, x + k1/2, u, t + dt/2)*dt
	x + k2
end

function integrate(::Type{RK3}, model::AbstractModel, x::StaticVector, u::StaticVector,
		t, dt)
    k1 = dynamics(model, x,             u, t       )*dt;
    k2 = dynamics(model, x + k1/2,      u, t + dt/2)*dt;
    k3 = dynamics(model, x - k1 + 2*k2, u, t + dt  )*dt;
    x + (k1 + 4*k2 + k3)/6
end

function integrate(::Type{RK4}, model::AbstractModel, x::StaticVector, u::StaticVector, t, dt)
	k1 = dynamics(model, x,        u, t       )*dt
	k2 = dynamics(model, x + k1/2, u, t + dt/2)*dt
	k3 = dynamics(model, x + k2/2, u, t + dt/2)*dt
	k4 = dynamics(model, x + k3,   u, t + dt  )*dt
	x + (k1 + 2k2 + 2k3 + k4)/6
end

function jacobian!(::Type{RK3}, ∇f::DynamicsJacobian{n,nm}, model, z::AbstractKnotPoint, tmp) where {n,nm}
	m = nm-n
	# n,m = size(model)
	x = state(z)
	u = control(z)
	dt = z.dt
	t = z.t

	z1 = StaticKnotPoint([x; u], z._x, z._u, dt, t)
	k1 = dynamics(model, z1)*dt
	z2 = StaticKnotPoint([x + 0.5*k1; u], z._x, z._u, dt, t + 0.5*dt)
	k2 = dynamics(model, z2)*dt
	z3 = StaticKnotPoint([x - k1 + 2*k2; u], z._x, z._u, dt, t + dt)
	k3 = dynamics(model, z3)*dt
	jacobian!(tmp[1], model, z1)
	jacobian!(tmp[2], model, z2)
	jacobian!(tmp[3], model, z3)
	A1 = SMatrix{n,n}(tmp[1].A)
	A2 = SMatrix{n,n}(tmp[2].A)
	A3 = SMatrix{n,n}(tmp[3].A)
	B1 = SMatrix{n,m}(tmp[1].B)
	B2 = SMatrix{n,m}(tmp[2].B)
	B3 = SMatrix{n,m}(tmp[3].B)

	I3 = Diagonal(@SVector ones(n))
	∇f.A .= I3 + (A1*dt + 4*A2*(I3 + 0.5*A1*dt)*dt + A3*(I3 - A1*dt + 2*A2*(I3 + 0.5*A1*dt)*dt)*dt)/6
	∇f.B .= (B1*dt +
             4*(0.5*A2*B1*dt + B2)*dt +
             (A3*(-B1*dt + A2*B1*dt*dt + 2*B2*dt) + B3)*dt)/6
end


text:
const TEST_NAME = "Threefry2x64"

include("common.jl")

using RandomNumbers.Random123

r = Threefry2x(UInt64, (123, 321))

test_all(r, 100_000_000)


text:
using ParametricFunctions
using Base.Test
using ContinuousTransformations
import ForwardDiff: derivative

include("utilities.jl")

include("test-accuracy.jl")
include("test-families.jl")
include("test-collocation.jl")
include("test-plots.jl")


text:
using ReachabilityAnalysis, Symbolics, Plots

@variables x y z
const positive_orthant = HPolyhedron([x >= 0, y >= 0, z >= 0], [x, y, z])

function prod_dest_verif(sol; T=100.0, target=10.0)

    # convert to a zonotopic flowpipe
    solz = overapproximate(sol, Zonotope)

    # project the final reach-set onto the space variables x, y, z
    X = project(solz(T), vars=(1, 2, 3))

    # check that all variables are nonnegative
    nonnegative = X ⊆ positive_orthant

    # compute the volume of the last reach-set
    H = overapproximate(X, Hyperrectangle)
    vol = volume(H)

    # check that that target belongs to the minkowski sum of the reach-sets projected in each coordinate
    B = convert(IntervalBox, H) # get the product-of-intervals representation
    contains_target = target ∈ sum(B)

    return nonnegative && contains_target, vol
end

@taylorize function prod_dest_I!(du, u, params, t)
    local a = 0.3
    x, y, z = u[1], u[2], u[3]

    du[1] = - (x * y) / (1 + x)
    du[2] = (x * y) / (1 + x) - a * y
    du[3] = a * y
    return du
end

X0 = (9.5 .. 10.0) × (0.01 .. 0.01) × (0.01 .. 0.01)
prob = @ivp(x'= prod_dest_I!(x), dim:3, x(0) ∈ X0)

solI = solve(prob, T=100.0, alg=TMJets(abstol=1e-11, orderT=7, orderQ=1));

property, vol = prod_dest_verif(solI)

plot(solI, vars=(0, 3), linecolor=:orange, color=:orange, alpha=0.3, lab="I")

@taylorize function prod_dest_IP!(du, u, params, t)
    x, y, z, a = u[1], u[2], u[3], u[4]

    du[1] = - (x * y) / (1 + x)
    du[2] = (x * y) / (1 + x) - a * y
    du[3] = a * y
    du[4] = zero(x)
    return du
end

X0 = (9.98 .. 9.98) × (0.01 .. 0.01) × (0.01 .. 0.01) × (0.296 .. 0.304)
prob = @ivp(x'= prod_dest_IP!(x), dim:4, x(0) ∈ X0)

solP = solve(prob, T=100.0, alg=TMJets(abstol=1e-12, orderT=7, orderQ=1));

property, vol = prod_dest_verif(solP)

plot(solP, vars=(0, 3), linecolor=:blue, color=:blue, alpha=0.3, lab="P")

X0 = (9.5 .. 10.0) × (0.01 .. 0.01) × (0.01 .. 0.01) × (0.296 .. 0.304)
prob = @ivp(x'= prod_dest_IP!(x), dim:4, x(0) ∈ X0)

solIP = solve(prob, T=100.0, alg=TMJets(abstol=1e-11, orderT=7, orderQ=1));

property, vol = prod_dest_verif(solIP)

plot(solIP, vars=(0, 3), linecolor=:red, color=:red, alpha=0.3, lab="I & P")

@taylorize function prod_dest_I_optimized!(du, u, params, t)
    local a = 0.3
    x, y, z = u[1], u[2], u[3]

    num = x * y
    den = 1 + x
    aux = num/den
    aux2 = a * y
    du[1] = -aux
    du[2] = aux - aux2
    du[3] = aux2
    return du
end

@taylorize function prod_dest_IP_optimized!(du, u, params, t)
    x, y, z, a = u[1], u[2], u[3], u[4]

    num = x * y
    den = 1 + x
    aux = num/den
    aux2 = a * y
    du[1] = -aux
    du[2] = aux - aux2
    du[3] = aux2
    du[4] = zero(x)
    return du
end


text:
export ADFW, ADFWPath

function ADFWPath(prob::Problem{L, B, CS},
                Cs::Vector{F},
                subsolver::AbstractNextBasisSolver=defaultsolver(prob);
                warmstart=true,
                kwargs...
                ) where {F, L, B<:AbBasis{F}, CS}
    enses = Vector{FWEnsemble{F,B}}()
    for (i, c) in enumerate(Cs)
        println("C = $c")
        ens0 = warmstart && (i > 1) ? enses[end] : FWEnsemble(prob)
        push!(enses, ADFW(prob, c, subsolver; ens=ens0, kwargs...))
    end
    enses
end

function ADFW(prob::Problem{L, B, CS},
              C::F,
              subsolver::AbstractNextBasisSolver=defaultsolver(prob);
              fc=true,                              # Fully Corrective
              ad=false,                             # Alternate Descent
              ens::FWEnsemble=FWEnsemble(prob),     # Warm Start
              maxiter::Int = 300,
              maxjter::Int = 10) where {L, B, CS, F<:AbstractFloat}
  ens = deepcopy(ens)
#  sanitycheck(ens)
  if length(ens) > 0            # warmstart
      @assert fc
      lassoonA!(ens, prob, C)
  end

  iter, jter = 0, 0
  currloss = trainerr(ens)
  # llog(txt="") = printfmtln("{:2d}:{:2d}) C={:5.0f} L={:8.4f} |β|={:.2f}({}) {}", iter, jter, C, currloss, sum(abs, ens.β), countnz(ens.β), txt)
  llog(t="") = nothing
  llog("NULL")

  for iter in 1:maxiter
 #   sanitycheck(ens)
    prevloss = currloss
    objtv, grad = lossngrad(prob, ens.fit)
    ωnew, ρ, score = getnextω(prob, subsolver, grad) # TODO : Get z
    push!(ens, ωnew, 0.)
    # sanitycheck(ens)

    if !fc
        fwvanillaupdate!(ens, iter, C)
  #      sanitycheck(ens)
    else
        lassoonA!(ens, prob, C)
#        sanitycheck(ens)
    end

    currloss = trainerr(ens)
    if ad
        jter = 0
        llog("new (ω, β)")
        if !fc                              # Vanilla FW w/ AD ⇒ One update only
            ens, currloss = backpropωs(prob, ens)
#            sanitycheck(ens)
        else                                # Corrective FW w/ AD ⇒ sub-iterations
            adconverged = false
            while jter < maxjter && !adconverged
                jter += 1
                beforeadloss = currloss
                newens, currloss = backpropωs(prob, ens)
                llog("ad {ω} ∇ωs = $(ssedir(ens, newens))")
                ens = newens
 #               sanitycheck(ens)
                lassoonA!(ens, prob, C)
 #               sanitycheck(ens)
                currloss = trainerr(ens)
                adconverged = beforeadloss-currloss<1e-3
                llog("ad β (converged=$adconverged)")
            end
        end
    end

    llog("END $iter")
    abs(prevloss-currloss) < 1e-5 && break
  end

  llog("\n<--END-->")
  #sanitycheck(ens)
  ens
end


text:
using AxisKeys
using MCMCChains
using ParetoSmooth
using Statistics
using Test
using Turing

import RData

let og_array = RData.load("Example_Log_Likelihood_Array.RData")["x"]
    global log_lik_arr = copy(permutedims(og_array, [3, 1, 2]))
end
let og_weights = RData.load("Weight_Matrix.RData")["weightMatrix"]
    global r_weights = exp.(permutedims(reshape(og_weights, 500, 2, 32), [3, 1, 2]))
end
r_eff = RData.load("Rel_Eff.RData")["rel_eff"]
r_psis = RData.load("Psis_Object.RData")["x"]
r_tail_len = Int.(RData.load("Tail_Vector.RData")["tail"])
pareto_k = RData.load("Pareto_K.RData")["pareto_k"]
r_loo = RData.load("Example_Loo.RData")["example_loo"]

# Add labels, reformat
r_pointwise = KeyedArray(
    r_loo["pointwise"][:, Not(4)];
    data = 1:size(r_loo["pointwise"], 1),
    statistic=[:loo_est, :mcse, :overfit, :pareto_k],
)

r_loo["estimates"] = hcat(r_loo["estimates"], r_loo["estimates"] / size(r_pointwise, 1))
r_ests = KeyedArray(
    r_loo["estimates"][Not(3), :];
    criterion=[:loo_est, :overfit],
    statistic=[:total, :se_total, :mean, :se_mean],
)

@testset "ParetoSmooth.jl" begin
    @testset "Basic Arrays" begin

        # All of these should run
        with_r_eff = psis(log_lik_arr, r_eff)
        jul_psis = psis(log_lik_arr)
        log_lik_mat = reshape(log_lik_arr, 32, 1000)
        chain_index = vcat(fill(1, 500), fill(2, 500))
        matrix_psis = psis(log_lik_mat; chain_index=chain_index)
        log_psis = psis(log_lik_arr; log_weights=true)

        jul_loo = loo(log_lik_arr)
        r_eff_loo = psis_loo(log_lik_arr, r_eff)

        @test display(jul_psis) === nothing
        @test display(jul_loo) === nothing

        # max 20% difference in tail length calc between Julia and R
        @test maximum(abs.(log.(jul_psis.tail_len ./ r_tail_len))) ≤ .2
        @test maximum(abs.(jul_psis.tail_len .- r_tail_len)) ≤ 10
        @test maximum(abs.(with_r_eff.tail_len .- r_tail_len)) ≤ 2

        # RMSE from R version is less than .1%
        @test sqrt(mean((with_r_eff.weights ./ r_weights .- 1).^2)) ≤ .001
        # RMSE less than .2% when using InferenceDiagnostics' ESS
        @test sqrt(mean((jul_psis.weights ./ r_weights .- 1).^2)) ≤ .002
        # Max difference is 1%
        @test maximum(log_psis.weights .- log.(r_weights)) ≤ .01

        ## Test difference in loo pointwise results

        # Different r_eff
        jul_pointwise = jul_loo.pointwise([:loo_est, :mcse, :overfit, :pareto_k])
        errs = (r_pointwise - jul_pointwise).^2
        @test sqrt(mean(errs(:loo_est))) ≤ .01
        @test sqrt(mean(errs(:overfit))) ≤ .01
        @test sqrt(mean(errs(:pareto_k))) ≤ .025
        errs_mcse = log.(r_pointwise(:mcse) ./ jul_loo.pointwise(:mcse)).^2
        @test_broken sqrt(mean(errs_mcse)) ≤ .1

        # Same r_eff
        r_eff_pointwise = r_eff_loo.pointwise([:loo_est, :mcse, :overfit, :pareto_k])
        errs = (r_pointwise - r_eff_pointwise).^2
        @test sqrt(mean(errs(:loo_est))) ≤ .01
        @test sqrt(mean(errs(:overfit))) ≤ .01
        @test sqrt(mean(errs(:pareto_k))) ≤ .025
        errs_mcse = log.(r_pointwise(:mcse) ./ r_eff_loo.pointwise(:mcse)).^2
        @test_broken sqrt(mean(errs_mcse)) ≤ .1

        # Test estimates
        errs = r_ests - jul_loo.estimates(criterion=[:loo_est, :overfit])
        display(r_ests)
        display(jul_loo.estimates(criterion=[:loo_est, :overfit]))
        display(errs)
        @test maximum(abs.(errs)) ≤ .01

        errs = r_ests - r_eff_loo.estimates(criterion=[:loo_est, :overfit])
        display(errs)
        @test maximum(abs.(errs)) ≤ .01

        # Test for calling correct method
        @test jul_loo.psis_object.weights ≈ psis(-log_lik_arr).weights
        @test r_eff_loo.psis_object.weights ≈ psis(-log_lik_arr, r_eff).weights
    end

    @testset "MCMCChains and Turing utilities" begin
        using Distributions, Random
        Random.seed!(112)
        # simulated samples for μ
        samples = randn(100, 1, 1)
        data = randn(50)
        chain = Chains(samples)

        compute_loglike(μ, data) = logpdf(Normal(μ, 1), data)
        compute_loglike(μ, σ, data) = logpdf(Normal(μ, σ), data)

        pll1 = pointwise_log_likelihoods(compute_loglike, chain, data)
        pll2 = pointwise_log_likelihoods(compute_loglike, samples, data)
        # the pointwise log likehoods should be the same for both methods
        @test pll1 ≈ pll2 atol = 1e-6
        # test the dimensions: data points, samples, chains
        @test size(pll1) == (50, 100, 1)
        # test that sum of pointwise log likelihoods equals sum of log likelihoods
        @test sum(sum(map(s->logpdf.(Normal(s, 1), data), samples))) ≈ sum(pll1) atol = 1e-6
        # test that psis_loo works with MCMCChains and yields correct type
        psis_loo_output = psis_loo(compute_loglike, chain, data)
        @test isa(psis_loo_output, PsisLoo)
        # test that loo works with MCMCChains and yields correct type
        psis_output = loo(compute_loglike, chain, data)
        @test isa(psis_output, PsisLoo)
        # test that psis works with MCMCChains and yields correct type
        psis_output = psis(compute_loglike, chain, data)
        @test isa(psis_output, Psis)

        # ensure that methods work with r_eff argument
        r_eff = similar(pll2, 0)
        # test that psis_loo works with MCMCChains and yields correct type
        psis_loo_output = psis_loo(compute_loglike, chain, data, r_eff)
        @test isa(psis_loo_output, PsisLoo)
        # test that loo works with MCMCChains and yields correct type
        psis_output = loo(compute_loglike, chain, data, r_eff)
        @test isa(psis_output, PsisLoo)
        # test that psis works with MCMCChains and yields correct type
        psis_output = psis(compute_loglike, chain, data, r_eff)
        @test isa(psis_output, Psis)

        @model function model(data)
            μ ~ Normal(0, 1)
            σ ~ truncated(Cauchy(0, 1), 0, Inf)
            for i in eachindex(data)
                data[i] ~ Normal(μ, σ)
            end
        end

        chain = sample(model(data), NUTS(1000, .9), MCMCThreads(), 1000, 4)
        pw_lls_turing = pointwise_log_likelihoods(model(data), chain)
        pw_lls_loglike = pointwise_log_likelihoods(compute_loglike, chain, data)

        # test the dimensions: data points, samples, chains
        @test size(pw_lls_turing) == (50, 1000, 4)
        # test that sum of pointwise log likelihoods equals sum of log likelihoods
        turing_samples = Array(Chains(chain, :parameters).value)
        # make this more terse with @tullio or other method later
        LL = 0.0
        n_samples,n_parms,n_chains = size(turing_samples)
        for s in 1:n_samples
            for c in 1:n_chains
                LL += sum(logpdf.(Normal(turing_samples[s,:,c]...), data))
            end
        end
        @test LL ≈ sum(pw_lls_turing) atol = 1e-6
        # Turing should work the same as compute_loglike
        @test pw_lls_loglike ≈ pw_lls_turing atol = 1e-6

        # test that psis_loo works with Turing model and MCMCChains and yields correct type
        psis_loo_output = psis_loo(model(data), chain)
        @test isa(psis_loo_output, PsisLoo)
        # test that loo works with Turing model and MCMCChains and yields correct type
        psis_output = loo(model(data), chain)
        @test isa(psis_output, PsisLoo)
        # test that psis works with Turing model and MCMCChains and yields correct type
        psis_output = psis(model(data), chain)
        @test isa(psis_output, Psis)

        # ensure that methods work with r_eff argument
        r_eff = similar(pw_lls_turing, 0)
        # test that psis_loo works with Turin and gives correct type
        psis_loo_output = psis_loo(model(data), chain, r_eff)
        @test isa(psis_loo_output, PsisLoo)
        # test that loo works with Turing model and MCMCChains and gives correct type
        psis_output = loo(model(data), chain, r_eff)
        @test isa(psis_output, PsisLoo)
        # test that psis works with Turing model and MCMCChains and gives correct type
        psis_output = psis(model(data), chain, r_eff)
        @test isa(psis_output, Psis)

    end
end