text:
# Join layer
"""
    Join(dim::Int64)
    Join(dim = dim::Int64)

Concatenates a tuple of arrays along a dimension `dim`. A convenient and type stable way of using `x -> cat(x..., dims = dim)`.
"""
struct Join{D}
    dim::Int64
    function Join(dim)
        if dim>3
            throw(DimensionMismatch("Dimension should be 1, 2 or 3."))
        end
        new{dim}(dim)
    end
    function Join(;dim)
        if dim>3
            throw(DimensionMismatch("Dimension should be 1, 2 or 3."))
        end
        new{dim}(dim)
    end
end
(m::Join{D})(x::NTuple{N,AbstractArray}) where {D,N} = cat(x...,dims = Val(D))

function Base.show(io::IO, l::Join)
    print(io, "Join(", "dim = ",l.dim, ")")
end

# Split layer
"""
    Split(outputs::Int64, dim::Int64)
    Split(outputs::Int64, dim = dim::Int64)

Breaks an array into a number of arrays which is equal to `outputs` along a dimension `dim`. `dim` should we divisible by `outputs` without a remainder.
"""
struct Split{O,D}
    outputs::Int64
    dim::Int64
    function Split(outputs,dim)
        if dim>3
            throw(DimensionMismatch("Dimension should be 1, 2 or 3."))
        elseif outputs<2
            throw(DomainError(outputs, "The number of outputs should be 2 or more."))
        end
        new{outputs,dim}(outputs,dim)
    end
    function Split(outputs;dim)
        if dim>3
            throw(DimensionMismatch("Dimension should be 1, 2 or 3."))
        elseif outputs<2
            throw(DomainError(outputs, "The number of outputs should be 2 or more."))
        end
        new{outputs,dim}(outputs,dim)
    end
end

function Split_func(x::T,m::Split{O,D}) where {O,D,T<:AbstractArray{<:AbstractFloat,2}}
    if D!=1
        throw(DimensionMismatch("Dimension should be 1."))
    end
    step_var = Int64(size(x, D) / O)
    f = i::Int64 -> (1+(i-1)*step_var):(i)*step_var
    vals = ntuple(i -> i, O)
    inds_vec = map(f,vals)
    x_out = map(inds -> x[inds,:],inds_vec)
    return x_out
end

function get_part(x::T,inds::NTuple{O,UnitRange{Int64}},d::Type{Val{1}}) where {O,T<:AbstractArray{<:AbstractFloat,4}}
    x_out = map(inds -> x[inds,:,:,:],inds)
end
function get_part(x::T,inds::NTuple{O,UnitRange{Int64}},d::Type{Val{2}}) where {O,T<:AbstractArray{<:AbstractFloat,4}}
    x_out = map(inds -> x[:,inds,:,:],inds)
end
function get_part(x::T,inds::NTuple{O,UnitRange{Int64}},d::Type{Val{3}}) where {O,T<:AbstractArray{<:AbstractFloat,4}}
    x_out = map(inds -> x[:,:,inds,:],inds)
end
function Split_func(x::T,m::Split{O,D}) where {O,D,T<:AbstractArray{<:AbstractFloat,4}}
    step_var = Int64(size(x, D) / O)
    f = i::Int64 -> (1+(i-1)*step_var):(i)*step_var
    vals = ntuple(i -> i, O)
    inds_tuple = map(f,vals)
    x_out = get_part(x,inds_tuple,Val{D})
    return x_out
end
(m::Split{O,D})(x::T) where {O,D,T<:AbstractArray} = Split_func(x,m)

function Base.show(io::IO, l::Split)
    print(io, "Split(", l.outputs,", dim = ",l.dim, ")")
end

# Makes Parallel layer type stable when used after Split
(m::Parallel)(xs::NTuple{N,AbstractArray}) where N = map((f,x) -> f(x), m.layers,xs)

# Addition layer
"""
    Addition()

A convenient way of using `x -> sum(x)`.
"""
struct Addition
end
(m::Addition)(x::NTuple{N,AbstractArray}) where N = sum(x)

# Activation layer
"""
    Activation(f::Function)

A convenient way of using `x -> f(x)`.
"""
struct Activation{F}
    f::F
    Activation(f::Function) = new{typeof(f)}(f)
end
(m::Activation{F})(x::AbstractArray) where F = m.f.(x)

function Base.show(io::IO, l::Activation)
    print(io, "Activation(",l.f, ")")
end

# Flatten layer
"""
    Flatten()

Flattens an array. A convenient way of using `x -> Flux.flatten(x)`.
"""
struct Flatten
end
(m::Flatten)(x::AbstractArray) = Flux.flatten(x)

# Identity layer
"""
    Identity()

Returns its input without changes. Should be used with a `Parallel` layer if one wants to have a branch that does not change its input.
"""
struct Identity
end
(m::Identity)(x::AbstractArray) = x


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
"""
    run_mpc(m::JuMP.AbstractModel, fp::String)

Solve the model predictive control problem using the `MPCScenario` defined in the JSON file stored at the file path `fp`.

Returns a Dict of results with keys matching those in the `MPCScenario`.
"""
function run_mpc(m::JuMP.AbstractModel, fp::String)
	s = MPCScenario(JSON.parsefile(fp))
	run_mpc(m, MPCInputs(s))
end

"""
    run_mpc(m::JuMP.AbstractModel,  d::Dict)

Solve the model predictive control problem using the `MPCScenario` defined in the dict `d`.

Returns a Dict of results with keys matching those in the `MPCScenario`.
"""
function run_mpc(m::JuMP.AbstractModel, d::Dict)
	s = MPCScenario(d)
	run_mpc(m, MPCInputs(s))
end

"""
    run_mpc(m::JuMP.AbstractModel, p::MPCInputs)

Solve the model predictive control problem using the `MPCInputs`.

Returns a Dict of results with keys matching those in the `MPCScenario`.
"""
function run_mpc(m::JuMP.AbstractModel, p::MPCInputs)
    build_mpc!(m, p)

    if !p.s.settings.add_soc_incentive
		@objective(m, Min, m[:Costs])
	else # Keep SOC high
		@objective(m, Min, m[:Costs] - sum(m[:dvStoredEnergy]["ElectricStorage", ts] for ts in p.time_steps) /
									   (8760. / p.hours_per_time_step)
		)
	end

	@info "Model built. Optimizing..."
	tstart = time()
	optimize!(m)
	opt_time = round(time() - tstart, digits=3)
	if termination_status(m) == MOI.TIME_LIMIT
		status = "timed-out"
    elseif termination_status(m) == MOI.OPTIMAL
        status = "optimal"
    else
		status = "not optimal"
		@warn "MPC solved with " termination_status(m), ", returning the model."
		return m
	end
	@info "MPC solved with " termination_status(m)
	@info "Solving took $(opt_time) seconds."

	tstart = time()
	results = mpc_results(m, p)
	time_elapsed = time() - tstart
	@info "Results processing took $(round(time_elapsed, digits=3)) seconds."
	results["status"] = status
	results["solver_seconds"] = opt_time
	return results
end

"""
    build_mpc!(m::JuMP.AbstractModel, p::MPCInputs)

Add variables and constraints for model predictive control model.
Similar to a REopt model but with any length of horizon (instead of one calendar year),
and the DER sizes must be provided.
"""
function build_mpc!(m::JuMP.AbstractModel, p::MPCInputs)
    add_variables!(m, p)

	for ts in p.time_steps_without_grid

		fix(m[:dvGridPurchase][ts], 0.0, force=true)

		for t in p.s.storage.types.elec
			fix(m[:dvGridToStorage][t, ts], 0.0, force=true)
		end

		for t in p.techs.elec, u in p.export_bins_by_tech[t]
			fix(m[:dvProductionToGrid][t, u, ts], 0.0, force=true)
		end
	end

	for b in p.s.storage.types.all
		if p.s.storage.attr[b].size_kw == 0 || p.s.storage.attr[b].size_kwh == 0
			@constraint(m, [ts in p.time_steps], m[:dvStoredEnergy][b, ts] == 0)
			@constraint(m, [t in p.techs.elec, ts in p.time_steps_with_grid],
						m[:dvProductionToStorage][b, t, ts] == 0)
			@constraint(m, [ts in p.time_steps], m[:dvDischargeFromStorage][b, ts] == 0)
			if b in p.s.storage.types.elec
				@constraint(m, [ts in p.time_steps], m[:dvGridToStorage][b, ts] == 0)
			end
		else
			add_general_storage_dispatch_constraints(m, p, b)
			if b in p.s.storage.types.elec
				add_elec_storage_dispatch_constraints(m, p, b)
			elseif b in p.s.storage.types.hot
				add_hot_thermal_storage_dispatch_constraints(m, p, b)
			elseif b in p.s.storage.types.cold
				add_cold_thermal_storage_dispatch_constraints(m, p, b)
			else
				@error("Invalid storage does not fall in a thermal or electrical set")
			end
		end
	end

	if any(size_kw->size_kw > 0, (p.s.storage.attr[b].size_kw for b in p.s.storage.types.all))
		add_storage_sum_constraints(m, p)
	end

	add_production_constraints(m, p)

	if !isempty(p.techs.no_turndown)
		@constraint(m, [t in p.techs.no_turndown, ts in p.time_steps],
            m[:dvRatedProduction][t,ts] == m[:dvSize][t]
        )
	end

	add_elec_load_balance_constraints(m, p)

	if !isempty(p.s.limits.grid_draw_limit_kw_by_time_step)
		add_grid_draw_limits(m, p)
	end

	if !isempty(p.s.electric_tariff.export_bins)
		add_export_constraints(m, p)
		if !isempty(p.s.limits.export_limit_kw_by_time_step)
			add_export_limits(m, p)
		end
	end

	if !isempty(p.s.electric_tariff.monthly_demand_rates)
		add_monthly_peak_constraint(m, p)
	end

	if !isempty(p.s.electric_tariff.tou_demand_ratchet_time_steps)
		add_tou_peak_constraint(m, p)
	end

	if !(p.s.electric_utility.allow_simultaneous_export_import) & !isempty(p.s.electric_tariff.export_bins)
		add_simultaneous_export_import_constraint(m, p)
	end

    m[:TotalFuelCosts] = 0.0
    m[:TotalPerUnitProdOMCosts] = 0.0

    if !isempty(p.techs.gen)
        add_gen_constraints(m, p)
		m[:TotalPerUnitProdOMCosts] += @expression(m,
			sum(p.s.generator.om_cost_per_kwh * p.hours_per_time_step *
			m[:dvRatedProduction][t, ts] for t in p.techs.gen, ts in p.time_steps)
		)
        m[:TotalGenFuelCosts] = @expression(m,
            sum(m[:dvFuelUsage][t,ts] * p.s.generator.fuel_cost_per_gallon for t in p.techs.gen, ts in p.time_steps)
        )
        m[:TotalFuelCosts] += m[:TotalGenFuelCosts]
	end

	add_elec_utility_expressions(m, p)
    add_previous_monthly_peak_constraint(m, p)
    add_previous_tou_peak_constraint(m, p)

    # TODO: random outages in MPC?
	if !isempty(p.s.electric_utility.outage_durations)
		add_dv_UnservedLoad_constraints(m,p)
		add_outage_cost_constraints(m,p)
		add_MG_production_constraints(m,p)
		if !isempty(p.s.storage.types.elec)
			add_MG_storage_dispatch_constraints(m,p)
		else
			fix_MG_storage_variables(m,p)
		end
		add_cannot_have_MG_with_only_PVwind_constraints(m,p)
		add_MG_size_constraints(m,p)

		if !isempty(p.techs.gen)
			add_MG_fuel_burn_constraints(m,p)
			add_binMGGenIsOnInTS_constraints(m,p)
		else
			m[:ExpectedMGFuelUsed] = 0
			m[:ExpectedMGFuelCost] = 0
			@constraint(m, [s in p.s.electric_utility.scenarios, tz in p.s.electric_utility.outage_start_time_steps, ts in p.s.electric_utility.outage_time_steps],
				m[:binMGGenIsOnInTS][s, tz, ts] == 0
			)
		end

		if p.s.site.min_resil_time_steps > 0
			add_min_hours_crit_ld_met_constraint(m,p)
		end
	end

	#################################  Objective Function   ########################################
	@expression(m, Costs,

		# Variable O&M
		m[:TotalPerUnitProdOMCosts] +

		# Total Generator Fuel Costs
        m[:TotalFuelCosts] +

		# Utility Bill
		m[:TotalElecBill]
	);
	if !isempty(p.s.electric_utility.outage_durations)
		add_to_expression!(Costs, m[:ExpectedOutageCost] + m[:mgTotalTechUpgradeCost] + m[:dvMGStorageUpgradeCost] + m[:ExpectedMGFuelCost])
	end
    #= Note: 0.9999*MinChargeAdder in Objective b/c when TotalMinCharge > (TotalEnergyCharges + TotalDemandCharges + TotalExportBenefit + TotalFixedCharges)
		it is arbitrary where the min charge ends up (eg. could be in TotalDemandCharges or MinChargeAdder).
		0.0001*MinChargeAdder is added back into LCC when writing to results.  =#
	nothing
end

function add_variables!(m::JuMP.AbstractModel, p::MPCInputs)
    @variables m begin
		# dvSize[p.techs.all] >= 0  # System Size of Technology t [kW]
		# dvPurchaseSize[p.techs.all] >= 0  # system kW beyond existing_kw that must be purchased
		dvGridPurchase[p.time_steps] >= 0  # Power from grid dispatched to meet electrical load [kW]
		dvRatedProduction[p.techs.all, p.time_steps] >= 0  # Rated production of technology t [kW]
		dvCurtail[p.techs.all, p.time_steps] >= 0  # [kW]
		dvProductionToStorage[p.s.storage.types.all, p.techs.all, p.time_steps] >= 0  # Power from technology t used to charge storage system b [kW]
		dvDischargeFromStorage[p.s.storage.types.all, p.time_steps] >= 0 # Power discharged from storage system b [kW]
		dvGridToStorage[p.s.storage.types.elec, p.time_steps] >= 0 # Electrical power delivered to storage by the grid [kW]
		dvStoredEnergy[p.s.storage.types.all, 0:p.time_steps[end]] >= 0  # State of charge of storage system b
		dvStoragePower[p.s.storage.types.all] >= 0   # Power capacity of storage system b [kW]
		dvStorageEnergy[p.s.storage.types.all] >= 0   # Energy capacity of storage system b [kWh]
		# TODO rm dvStoragePower/Energy dv's
		dvPeakDemandTOU[p.ratchets, 1:1] >= 0  # Peak electrical power demand during ratchet r [kW]
		dvPeakDemandMonth[p.months] >= 0  # Peak electrical power demand during month m [kW]
		# MinChargeAdder >= 0
	end
	# TODO: tiers in MPC tariffs and variables?

	if !isempty(p.s.electric_tariff.export_bins)
		@variable(m, dvProductionToGrid[p.techs.elec, p.s.electric_tariff.export_bins, p.time_steps] >= 0)
	end

    m[:dvSize] = p.existing_sizes

	for b in p.s.storage.types.all
		fix(m[:dvStoragePower][b], p.s.storage.attr["ElectricStorage"].size_kw, force=true)
		fix(m[:dvStorageEnergy][b], p.s.storage.attr["ElectricStorage"].size_kwh, force=true)
	end

	# not modeling min charges since control does not affect them
    m[:MinChargeAdder] = 0

	if !isempty(p.techs.gen)  # Problem becomes a MILP
		@warn """Adding binary variable to model gas generator.
				 Some solvers are very slow with integer variables"""
		@variables m begin
			dvFuelUsage[p.techs.gen, p.time_steps] >= 0 # Fuel burned by technology t in each time step [kWh]
			binGenIsOnInTS[p.techs.gen, p.time_steps], Bin  # 1 If technology t is operating in time step h; 0 otherwise
		end
	end

	if !(p.s.electric_utility.allow_simultaneous_export_import)
		@warn """Adding binary variable to prevent simultaneous grid import/export.
				 Some solvers are very slow with integer variables"""
		@variable(m, binNoGridPurchases[p.time_steps], Bin)
	end

	if !isempty(p.s.electric_utility.outage_durations) # add dvUnserved Load if there is at least one outage
		@warn """Adding binary variable to model outages.
				 Some solvers are very slow with integer variables"""
		max_outage_duration = maximum(p.s.electric_utility.outage_durations)
		outage_time_steps = p.s.electric_utility.outage_time_steps
		tZeros = p.s.electric_utility.outage_start_time_steps
		S = p.s.electric_utility.scenarios
		# TODO: currently defining more decision variables than necessary b/c using rectangular arrays, could use dicts of decision variables instead
		@variables m begin # if there is more than one specified outage, there can be more othan one outage start time
			dvUnservedLoad[S, tZeros, outage_time_steps] >= 0 # unserved load not met by system
			dvMGProductionToStorage[p.techs.all, S, tZeros, outage_time_steps] >= 0 # Electricity going to the storage system during each time_step
			dvMGDischargeFromStorage[S, tZeros, outage_time_steps] >= 0 # Electricity coming from the storage system during each time_step
			dvMGRatedProduction[p.techs.all, S, tZeros, outage_time_steps]  # MG Rated Production at every time_step.  Multiply by ProdFactor to get actual energy
			dvMGStoredEnergy[S, tZeros, 0:max_outage_duration] >= 0 # State of charge of the MG storage system
			dvMaxOutageCost[S] >= 0 # maximum outage cost dependent on number of outage durations
			# dvMGTechUpgradeCost[p.techs.all] >= 0
			# dvMGStorageUpgradeCost >= 0
			# dvMGsize[p.techs.all] >= 0

			dvMGFuelUsed[p.techs.all, S, tZeros] >= 0
			dvMGMaxFuelUsage[S] >= 0
			dvMGMaxFuelCost[S] >= 0
			dvMGCurtail[p.techs.all, S, tZeros, outage_time_steps] >= 0

			# binMGStorageUsed, Bin # 1 if MG storage battery used, 0 otherwise
			# binMGTechUsed[p.techs.all], Bin # 1 if MG tech used, 0 otherwise
			binMGGenIsOnInTS[S, tZeros, outage_time_steps], Bin
		end
	end
end

ext: jl

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
#############################################################################
#############################################################################
#
# This file contains several extensions for the `MutableBinaryMaxHeap` type
#
#############################################################################
#############################################################################

"""
Check that two heaps have the same values.
"""
==(h1::MutableBinaryMaxHeap, h2::MutableBinaryMaxHeap)::Bool = extract_all!(deepcopy(h1)) == extract_all!(deepcopy(h2))

"""
Apply a mapping to the heap values.
"""
function map!(f, h::MutableBinaryMaxHeap)
    for n in h.nodes
        update!(h, n.handle, f(n.value))
    end
    return h
end

"""
Apply a mapping to the heap values and create a new heap.
"""
map(f,h::MutableBinaryMaxHeap) = map!(f,deepcopy(h))

"""
Iterate over the heap. Implements the iteratable inteface.
"""
iterate(h::MutableBinaryMaxHeap, state=1) = state > length(h) ? nothing : (h.nodes[state].value, state+1)


text:
module PrototypeExample

mutable struct Account
    id::Int
    balance::Float64
end

struct Customer
    name::String
    savingsAccount::Account
    checkingAccount::Account
end

function sample_customers()
    a1 = Account(1, 100.0)
    a2 = Account(2, 200.0)
    c1 = Customer("John Doe", a1, a2)

    a3 = Account(3, 300.0)
    a4 = Account(4, 400.0)
    c2 = Customer("Brandon King", a1, a2)

    return [c1, c2]
end

function test(copy_function::Function)
    println("--- testing ", string(copy_function), " ---")
    customers = sample_customers()
    c = copy_function(customers)
    c[1].checkingAccount.balance += 500
    println("orig: ", customers[1].checkingAccount.balance)
    println("new:  ", c[1].checkingAccount.balance)
end

end #module

using .PrototypeExample
PrototypeExample.test(copy)
PrototypeExample.test(deepcopy)


text:
# extend basic ops to SIMD wrappers

# credit: basically copied from DoubleFloats.jl and AccurateArithmetic.jl

# aka twosum
@inline function add_(x::Vec{N,T}, y::Vec{N,T}) where {N, T<:AbstractFloat}
    hi  = x + y
    v = hi - x
    lo = (x - (hi - v)) + (y - v)
    return hi, lo
end

# aka quicktwosum
@inline function add_hilo_(x::Vec{N,T}, y::Vec{N,T}) where {N, T<:AbstractFloat}
    hi = x + y
    lo = y - (hi - x)
    return hi, lo
end

# aka twodiff
@inline function sub_(x::Vec{N,T}, y::Vec{N,T}) where {N, T<:AbstractFloat}
    hi  = x - y
    v = hi - x
    lo = (x - (hi - v)) - (y + v)
    return hi, lo
end

# aka quicktwodiff
@inline function sub_hilo_(x::Vec{N,T}, y::Vec{N,T}) where {N, T<:AbstractFloat}
    hi = x - y
    lo = (x - hi) - y
    return hi, lo
end

@inline function mul_(x::Vec{N,T}, y::Vec{N,T}) where {N, T<:AbstractFloat}
    hi = x * y
    lo = fma(x, y, -hi)
    return hi, lo
end

@inline function dfvadd(xhi::Vec{N,T}, xlo::Vec{N,T}, yhi::Vec{N,T}, ylo::Vec{N,T}) where {N,T<:AbstractFloat}
    hi, lo = add_(xhi, yhi)
    thi, tlo = add_(xlo, ylo)
    c = lo + thi
    hi, lo = add_hilo_(hi, c)
    c = tlo + lo
    hi, lo = add_hilo_(hi, c)
    return hi, lo
end

@inline function dfvsub(xhi::Vec{N,T}, xlo::Vec{N,T}, yhi::Vec{N,T}, ylo::Vec{N,T}) where {N,T<:AbstractFloat}
    hi, lo = sub_(xhi, yhi)
    thi, tlo = sub_(xlo, ylo)
    c = lo + thi
    hi, lo = add_hilo_(hi, c)
    c = tlo + lo
    hi, lo = add_hilo_(hi, c)
    return hi, lo
end

@inline function dfvmul(xhi::Vec{N,T}, xlo::Vec{N,T}, yhi::Vec{N,T}, ylo::Vec{N,T}) where {N,T<:AbstractFloat}
    hi, lo = mul_(xhi, yhi)
    t = xlo * ylo
    t = fma(xhi, ylo, t)
    t = fma(xlo, yhi, t)
    t = lo + t
    hi, lo = add_hilo_(hi, t)
    return hi, lo
end

max_stars_repo_name: RalphAS/DoubleBLAS


text:
using ClearStacktrace
using Test

module ModuleC
    func_c() = error("Error")
    end

    module ModuleB
    import ..ModuleC
    func_b(::String, ::Rational) = ModuleC.func_c()
    end

    module ModuleA
    import ..ModuleB
    func_a(::Int) = ModuleB.func_b("a", 1//3)
end

@testset "ClearStacktrace.jl" begin
    try
        ModuleA.func_a(1)
    catch e
        st = stacktrace(catch_backtrace())
        Base.show_backtrace(stdout, st)
    end
end


text:
module TTTexample
import Flux

using AlphaZero
using ..MuZero

using ParameterSchedulers: Scheduler, Cos

gspec = Examples.games["tictactoe"]

n=16
self_play = (;
  sim=SimParams(
    num_games=100,
    num_workers=n,
    batch_size=8,
    use_gpu=false, # not used, change device to Flux.gpu
    # todevice=Flux.cpu, #?
    reset_every=4, #not used, mcts resets everytime
    flip_probability=0.,
    alternate_colors=false),
  mcts = MctsParams(
    num_iters_per_turn=64, #1000 benchmark
    cpuct=1.25,
    temperature=ConstSchedule(1.0),
    dirichlet_noise_ϵ=0.25,
    dirichlet_noise_α=0.1),
    device=Flux.cpu) # Flux.cpu || Flux.gpu

arena = (;
  sim=SimParams(
    num_games=50,
    num_workers=n,
    batch_size=n÷2,
    use_gpu=false, # not used, change device to Flux.gpu
    reset_every=1,
    flip_probability=0., #0.5
    alternate_colors=true),
  mcts = MctsParams(
    self_play.mcts,
    temperature=ConstSchedule(0.3),
    dirichlet_noise_ϵ=0.1),
  update_threshold=0.0,
  device=Flux.cpu)# Flux.cpu || Flux.gpu

learning_params = (;
  num_unroll_steps=5, #if =0, g is not learning, muzero-general=20
  td_steps=20, # with max length=9, always go till the end of the game, rootvalues don't count
  discount=0.997,
  #// value_loss_weight = 0.25, #TODO
  l2_regularization=1f-4, #Float32
  #// l2_regularization=0f0, #Float32
  loss_computation_batch_size=512,
  batches_per_checkpoint=10,
  num_checkpoints=1,
  opt=Scheduler(
    Cos(λ0=3e-3, λ1=1e-5, period=10000), # cosine annealing, google period 2e4, generat doesn't use any
    Flux.ADAM()
  ),
  model_type = :mlp, # :mlp || :resnet
  device=Flux.cpu # Flux.cpu || Flux.gpu
)

benchmark_sim = SimParams(
    num_games=400,
    num_workers=n,
    batch_size=n÷2,
    use_gpu=false,
    reset_every=1,
    flip_probability=0.0, #0.5
    alternate_colors=true)

bench_mcts = MctsParams(
  num_iters_per_turn=400, #1000 benchmark
  cpuct=1.25,
  temperature=ConstSchedule(0.3),
  dirichlet_noise_ϵ=0.25,
  dirichlet_noise_α=0.1)

benchmark = (;
  vanilla_mcts=Benchmark.Duel(
    Mu(arena.mcts),
    Benchmark.MctsRollouts(bench_mcts),
    benchmark_sim),
  minmax_d5=Benchmark.Duel(
    Mu(arena.mcts),
    Benchmark.MinMaxTS(depth=5, amplify_rewards=true, τ=1.),
    benchmark_sim)
)

  μparams = MuParams(self_play, learning_params, arena, 5_000, 3000)

# muzero-general params
# μNetworkHP = MuNetworkHP(
#   gspec,
#   PredictionHP(hiddenstate_shape=32, width=64, depth_common=-1,
#     depth_vectorhead=0, depth_scalarhead=0, use_batch_norm=false, batch_norm_momentum=1.),
#   DynamicsHP(hiddenstate_shape=32, width=64, depth_common=-1,
#    depth_vectorhead=0, depth_scalarhead=0s,  use_batch_norm=false, batch_norm_momentum=1.),
#   RepresentationHP(hiddenstate_shape=32, width=0, depth=-1))

  hs_shape = 32
  simμNetworkHP = MuNetworkHP(
    gspec,
    PredictionSimpleHP(hiddenstate_shape=hs_shape, width=64, depth_common=-1,
      depth_vectorhead=0, depth_scalarhead=0, use_batch_norm=false, batch_norm_momentum=0.8f0),
    DynamicsSimpleHP(hiddenstate_shape=hs_shape, width=64, depth_common=-1,
     depth_vectorhead=0, depth_scalarhead=0,  use_batch_norm=false, batch_norm_momentum=0.8f0),
    RepresentationSimpleHP(hiddenstate_shape=hs_shape, width=0, depth=-1))

  env = MuEnv(gspec, μparams, MuNetwork(simμNetworkHP))

end # module


text:
#=
Popis
  Funkcia vypocita pre kazdy vrchol hodnotu jeho eigenvector centrality.

Parametre
  g : datova struktura grafu :: GenericGraph
  dist_key : nazov atributu, v ktorom je ulozena vlastnost vrchola :: String

Navratova hodnota
  eigenvector : slovnik obsahujuci mapovanie vrcholov a k nim prisluchajucich hodnot
           eigenvector centrality (index_vrchola => hodnota_eigenvector_centrality) :: Dict{Int64, Float64}()

Poznamky
  Eigen vector a eigen hodnoty su pocitane pomocou externeho procesu. Vzhladom na problemy implementacie
  je vysledok vypoctu nedeterministicky. Graf nesmie obsahovat paralelne hrany.
=#

function centrality_eigenvector(g::GenericGraph; dist_key::AbstractString="weight")

  if GraphAnalysis.contains_parallel_edges(g)
      println("Graf obsahuje paralelne hrany!")
      return false
  end

  matrix = sparse(weight_matrix(g,map(x->x.attributes[dist_key],edges(g))))
  matrix = matrix+matrix'

  (d, v) = eigs(matrix, nev=1, which=:LR, maxiter=50)

  largest = vec(v)
  x = vecnorm(largest)
  map!(v -> real(v), largest)
  norm = sign(sum(largest))*x

  eigenvector = DataStructures.OrderedDict(zip([v.index for v=vertices(g)],map(float,largest/norm)))

  return eigenvector

end

max_stars_repo_name: mpuk/GraphAnalysis

max_issues_repo_name: mpuk/GraphAnalysis

max_forks_repo_name: mpuk/GraphAnalysis


text:
abstract type AbstactFunJac{J2} end
mutable struct FunJac{F, F2, J, P, J2, uType, uType2} <: AbstactFunJac{J2}
    fun::F
    fun2::F2
    jac::J
    p::P
    jac_prototype::J2
    u::uType
    du::uType
    resid::uType2
end
FunJac(fun,jac,p,jac_prototype,u,du) = FunJac(fun,nothing,jac,p,jac_prototype,u,du,nothing)
FunJac(fun,jac,p,jac_prototype,u,du,resid) = FunJac(fun,nothing,jac,p,jac_prototype,u,du,resid)

function cvodefunjac(t::Float64,
                     u::N_Vector,
                     du::N_Vector,
                     funjac::FunJac)
    if !(pointer(funjac.u) === __N_VGetArrayPointer_Serial(u))
      #@warn "Pointer is broken to FunJac.u"
      funjac.u = convert(Vector, u)
    end
    _u = funjac.u

    if !(pointer(funjac.du) === __N_VGetArrayPointer_Serial(du))
      #@warn "Pointer is broken to FunJac.du"
      funjac.du = convert(Vector, du)
    end
    _du = funjac.du

    funjac.fun(_du, _u, funjac.p, t)
    return CV_SUCCESS
end

function cvodefunjac2(t::Float64,
                     u::N_Vector,
                     du::N_Vector,
                     funjac::FunJac)
    if !(pointer(funjac.u) === __N_VGetArrayPointer_Serial(u))
      #@warn "Pointer is broken to FunJac.u"
      funjac.u = convert(Vector, u)
    end
    _u = funjac.u

    if !(pointer(funjac.du) === __N_VGetArrayPointer_Serial(du))
      #@warn "Pointer is broken to FunJac.du"
      funjac.du = convert(Vector, du)
    end
    _du = funjac.du

    funjac.fun2(_du, _u, funjac.p, t)
    return CV_SUCCESS
end

function cvodejac(t::realtype,
                  u::N_Vector,
                  du::N_Vector,
                  J::SUNMatrix,
                  funjac::AbstactFunJac{Nothing},
                  tmp1::N_Vector,
                  tmp2::N_Vector,
                  tmp3::N_Vector)

    if !(pointer(funjac.u) === __N_VGetArrayPointer_Serial(u))
      #@warn "Pointer is broken to FunJac.u"
      funjac.u = convert(Vector, u)
    end
    _u = funjac.u

    #=
    if !(pointer(funjac.du) === __N_VGetArrayPointer_Serial(du))
      #@warn "Pointer is broken to FunJac.du"
      _du = convert(Vector, du)
    end
    =#

    funjac.jac(convert(Matrix, J), _u, funjac.p, t)
    return CV_SUCCESS
end

function cvodejac(t::realtype,
                  u::N_Vector,
                  du::N_Vector,
                  _J::SUNMatrix,
                  funjac::AbstactFunJac{<:SparseMatrixCSC},
                  tmp1::N_Vector,
                  tmp2::N_Vector,
                  tmp3::N_Vector)
    jac_prototype = funjac.jac_prototype
    J = convert(SparseMatrixCSC,_J)

    if !(pointer(funjac.u) === __N_VGetArrayPointer_Serial(u))
      #@warn "Pointer is broken to FunJac.u"
      funjac.u = convert(Vector, u)
    end
    _u = funjac.u

    funjac.jac(jac_prototype, _u, funjac.p, t)
    J.nzval .= jac_prototype.nzval
    # Sundials resets the value pointers each time, so reset it too
    @. J.rowval = jac_prototype.rowval - 1
    @. J.colptr = jac_prototype.colptr - 1
    return CV_SUCCESS
end

function idasolfun(t::Float64, u::N_Vector, du::N_Vector, resid::N_Vector, funjac::FunJac)
    if !(pointer(funjac.u) === __N_VGetArrayPointer_Serial(u))
      #@warn "Pointer is broken to FunJac.u"
      _u = convert(Vector, u)
    else
      _u = funjac.u
    end

    if !(pointer(funjac.du) === __N_VGetArrayPointer_Serial(du))
      #@warn "Pointer is broken to FunJac.du"
      _du = convert(Vector, du)
    else
      _du = funjac.du
    end

    funjac.fun(funjac.resid, _du, _u, funjac.p, t)

    unsafe_copyto!(Sundials.__N_VGetArrayPointer_Serial(resid),
                   pointer(funjac.resid),
                   length(funjac.resid))

    return IDA_SUCCESS
end

function idajac(t::realtype,
                cj::realtype,
                u::N_Vector,
                du::N_Vector,
                res::N_Vector,
                J::SUNMatrix,
                funjac::AbstactFunJac{Nothing},
                tmp1::N_Vector,
                tmp2::N_Vector,
                tmp3::N_Vector)

    if !(pointer(funjac.u) === __N_VGetArrayPointer_Serial(u))
      #@warn "Pointer is broken to FunJac.u"
      _u = convert(Vector, u)
    else
      _u = funjac.u
    end

    if !(pointer(funjac.du) === __N_VGetArrayPointer_Serial(du))
      #@warn "Pointer is broken to FunJac.du"
      _du = convert(Vector, du)
    else
      _du = funjac.du
    end

    funjac.jac(convert(Matrix, J), _du,
               _u, funjac.p, cj, t)
    return IDA_SUCCESS
end

function idajac(t::realtype,
                cj::realtype,
                u::N_Vector,
                du::N_Vector,
                res::N_Vector,
                _J::SUNMatrix,
                funjac::AbstactFunJac{<:SparseMatrixCSC},
                tmp1::N_Vector,
                tmp2::N_Vector,
                tmp3::N_Vector)

  jac_prototype = funjac.jac_prototype
  J = convert(SparseMatrixCSC,_J)

  if !(pointer(funjac.u) === __N_VGetArrayPointer_Serial(u))
    #@warn "Pointer is broken to FunJac.u"
    _u = convert(Vector, u)
  else
    _u = funjac.u
  end

  if !(pointer(funjac.du) === __N_VGetArrayPointer_Serial(du))
    #@warn "Pointer is broken to FunJac.du"
    _du = convert(Vector, du)
  else
    _du = funjac.du
  end

  funjac.jac(jac_prototype, _du, convert(Vector, _u), funjac.p, cj, t)
  J.nzval .= jac_prototype.nzval
  # Sundials resets the value pointers each time, so reset it too
  @. J.rowval = jac_prototype.rowval - 1
  @. J.colptr = jac_prototype.colptr - 1

  return IDA_SUCCESS
end


text:
# IupGetParam: Example
#  Shows a dialog with all the possible fields.

module getparam_

using IUP

export
	getparam

include("../src/libiup_h.jl")
include("../src/iupcontrols.jl")

#------------------------example html/examples/C/getparam.c ----------------------------------
function getparam()

	pboolean::Cint = 1
	pinteger::Cint = 3456
	preal::Float32 = 3.543
	pinteger2::Cint = 192
	preal2::Float32 = 0.5
	pangle::Float32 = 90
	pstring = "string text"
	pcolor = "255 0 128"
	plist::Cint = 2
	pstring2 = "second text\nsecond line"
	pfile_name = "test.jpg"

	IupOpen()                       #Initializes IUP
	IupControlsOpen()   			# so we can show the IupDial, but can be tested without it
	IupSetLanguage("ENGLISH")

	if (IupGetParam("Title", param_action, 0,
	               "Boolean: %b[No,Yes]{Boolean Tip}\n
	               Integer: %i{Integer Tip}\n
	               Real 1: %r{Real Tip}\n
	               Sep1 %t\n
	               Integer: %i[0,255]{Integer Tip 2}\n
	               Real 2: %r[-1.5,1.5]{Real Tip 2}\n
	               Sep2 %t\n
	               Angle: %a[0,360]{Angle Tip}\n
	               String: %s{String Tip}\n
	               List: %l|item1|item2|item3|{List Tip}\n
	               File: %f[OPEN|*.bmp;*.jpg|CURRENT|NO|NO]{File Tip}\n
	               Color: %c{Color Tip}\n
	               Sep3 %t\n
	               Multiline: %m{Multiline Tip}\n",
	               pointer(pboolean), pointer(pinteger), pointer(preal), pointer(pinteger2), pointer(preal2),
	               pointer(pangle), pstring, pointer(plist), pfile_name, pcolor, pstring2, C_NULL) != 0)
		return IUP_DEFAULT
	end

	IupMessagef("IupGetParam",
	          "Boolean Value: %d\n
	          Integer: %d\n
	          Real 1: %g\n
	          Integer: %d\n
	          Real 2: %g\n
	          Angle: %g\n
	          String: %s\n
	          List Index: %d\n
	          FileName: %s\n
	          Color: %s\n
	          Multiline: %s",
	          pboolean, pinteger, float64(preal), pinteger2, float64(preal2), float64(pangle),
	          pstring, plist, pfile_name, pcolor, pstring2)

	IupClose()                                  # And close it when ready
end

function param_action(dialog::Ptr{Ihandle}, param_index::Int, user_data::Ptr{Void})
	if (param_index == IUP_GETPARAM_OK)
		@printf("OK\n")
	elseif (param_index == IUP_GETPARAM_INIT)
		@printf("Map\n")
	elseif (param_index == IUP_GETPARAM_CANCEL)
		@printf("Cancel\n")
	elseif (param_index == IUP_GETPARAM_HELP)
		@printf("Help\n")
	else
		param_str = sprintf("PARAM%d", param_index)
		param = convert(Ptr{Ihandle}, IupGetAttribute(dialog, param_str))
		@printf("%s = %s\n", param_str, IupGetAttribute(param, "VALUE"))
	end
	return 1
end

end 	# module


text:
@testset "entropy_thresholding" begin
    img = testimage("cameraman")
    edges, histogram = build_histogram(img,256)
    @test_throws ErrorException find_threshold(EntropyThresholding(),histogram,edges)
    @test find_threshold(EntropyThresholding(),histogram[1:256],edges) == 0.76171875
end


text:
function test_mapreduce(AT)
    @testset "mapreduce" begin
        for ET in supported_eltypes()
            T = AT{ET}
            @testset "$ET" begin
                range = ET <: Integer ? (ET(-2):ET(2)) : ET
                @testset "mapreducedim" begin
                    for N in (2, 10)
                        @test compare(x -> sum(x, dims=2),      AT, rand(range, N, N))
                        @test compare(x -> sum(x, dims=1),      AT, rand(range, N, N))
                        @test compare(x -> sum(x, dims=(1, 2)), AT, rand(range, N, N))

                        @test compare(x -> sum(x, dims=2),      AT, rand(range, N, 10))
                        @test compare(x -> sum(x, dims=1),      AT, rand(range, N, 10))

                        @test compare(x -> sum(x, dims=2),      AT, rand(range, 10, N))
                        @test compare(x -> sum(x, dims=1),      AT, rand(range, 10, N))

                        _zero = zero(ET)
                        _addone(z) = z + one(z)
                        @test compare(x->mapreduce(_addone, +, x; dims = 2),
                                      AT, rand(range, N, N))
                        @test compare(x->mapreduce(_addone, +, x; dims = 2, init = _zero),
                                      AT, rand(range, N, N))
                    end
                end
                @testset "sum maximum minimum prod" begin
                    for dims in ((4048,), (1024,1024), (77,), (1923,209))
                        @test compare(sum,  AT, rand(range, dims))
                        @test compare(prod, AT, rand(range, dims))
                        @test compare(x -> sum(abs, x),  AT, rand(range, dims))
                        @test compare(x -> prod(abs, x), AT, rand(range, dims))
                        ET <: Complex || @test compare(maximum, AT,rand(range, dims))
                        ET <: Complex || @test compare(minimum, AT,rand(range, dims))
                    end
                end

                @testset "broadcasted arrays" begin
                    for dims in ((4048,), (1024,1024), (77,), (1923,209))
                        @test compare(x->mapreduce(z -> z + one(z), +,
                                                   Broadcast.Broadcasted(+, (x, x));
                                                   init = zero(ET)), AT, rand(range, dims))
                    end
                end
            end
        end
        @testset "any all ==" begin
            for Ac in ([false, false], [false, true], [true, true])
                A = AT(Ac)
                @test A isa AT{Bool,1}
                @test any(A) == any(Ac)
                @test all(A) == all(Ac)
                @test A == copy(A)
                @test A !== copy(A)
                @test A == deepcopy(A)
                @test A !== deepcopy(A)
            end
        end

        @testset "isapprox" begin
            for ET in supported_eltypes()
                ET <: Complex && continue
                A = fill(AT{ET}, ET(0), (100,))
                B = ones(AT{ET}, 100)
                @test !(A ≈ B)
                @test !(A ≈ Array(B))
                @test !(Array(A) ≈ B)

                ca = AT(randn(ComplexF64,3,3))
                cb = copy(ca)
                cb[1:1, 1:1] .+= 1e-7im
                @test isapprox(ca, cb, atol=1e-5)
                @test !isapprox(ca, cb, atol=1e-9)
            end
        end
    end
end


text:
using Documenter, Seismic

makedocs(
    modules = [Seismic],
    doctest = true)

deploydocs(
    deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo   = "github.com/msacchi/Seismic.jl.git",
    julia  = "0.5",
    osname = "osx")


text:
# Metida
# Copyright © 2019-2020 Vladimir Arnautov aka PharmCat <mail@pharmcat.net>

__precompile__()
module Metida

using Distributions, LinearAlgebra, StatsBase, ForwardDiff, CategoricalArrays, Polyester#, LoopVectorization
using Optim, LineSearches, MetidaBase
using StatsModels
import MetidaBase: Tables, MetidaModel, AbstractCovarianceStructure, AbstractCovmatMethod, AbstractCovarianceType, AbstractLMMDataBlocks, MetidaTable, metida_table, PrettyTables

import LinearAlgebra:checksquare
import StatsModels: @formula, termvars
import StatsBase: fit, fit!, coef, coefnames, confint, nobs, dof_residual, dof, loglikelihood, aic, bic, aicc, isfitted, vcov, stderror, modelmatrix, response
import Base:show

export @formula, @covstr,
SI, ScaledIdentity,
DIAG, Diag,
AR, Autoregressive,
ARH, HeterogeneousAutoregressive,
CS, CompoundSymmetry,
CSH, HeterogeneousCompoundSymmetry,
ARMA, AutoregressiveMovingAverage,
TOEP, Toeplitz,
TOEPP, ToeplitzParameterized,
TOEPH, HeterogeneousToeplitz,
TOEPHP, HeterogeneousToeplitzParameterized,
CovarianceType, CovmatMethod,
fit!, LMM, VarEffect, theta, logreml, m2logreml, thetalength, dof_satter, dof_contain, rankx, caic, lcontrast, typeiii, estimate, contrast,
gmatrix, rmatrix, vmatrix!,
AbstractCovarianceType, AbstractCovmatMethod, MetidaModel,
getlog

export coef, coefnames, confint, nobs, dof_residual, dof, loglikelihood, aic, bic, aicc, isfitted, vcov, stderror, modelmatrix, response

const LDCORR = sqrt(eps())
const LOGLDCORR = log(sqrt(eps()))
const NEWTON_OM = Optim.Newton(;alphaguess = LineSearches.InitialHagerZhang(), linesearch = LineSearches.HagerZhang())
const LBFGS_OM  = Optim.LBFGS(;alphaguess = LineSearches.InitialStatic(), linesearch = LineSearches.Static())

include("sweep.jl")
include("varstruct.jl")
include("gmat.jl")
include("rmat.jl")
include("linearalgebra.jl")
include("options.jl")
include("modelresult.jl")
include("lmmdata.jl")
include("lmm.jl")
include("reml.jl")
include("ml.jl")
include("fit.jl")
include("showutils.jl")
include("statsbase.jl")
include("utils.jl")
include("dof_satter.jl")
include("dof_contain.jl")
include("fvalue.jl")
include("typeiii.jl")
include("estimate.jl")

end # module


text:
#river.msh geometry, created by JIGSAW
#adapted by @miguelraz 28/03/2020 from
#https://github.com/dengwirda/jigsaw-matlab
#mshid=1
#ndims=2
#point=342
const river = [(40388958.94,3568824.25)
(40389408.97,3568742.25)
(40389692.69,3568631.75)
(40390153.88,3568354.5)
(40390541.16,3568103.25)
(40390754.84,3567948)
(40390972,3567770.25)
(40391300,3567432.75)
(40391416.09,3567291.25)
(40391583.16,3567037.5)
(40391692.88,3566873.25)
(40391845.41,3566605.75)
(40391936.72,3566417.25)
(40392189.81,3565992.25)
(40392461.75,3565558.25)
(40392562.94,3565290.75)
(40392767.97,3564852.25)
(40392867.91,3564665)
(40393000.97,3564406.25)
(40393086.63,3564221.75)
(40393205.91,3563958.5)
(40393275.56,3563767.75)
(40393382.03,3563490.25)
(40393416.16,3563308.5)
(40393526.88,3563056)
(40393831.52,3563001.625)
(40393933.06,3562983.5)
(40394162.26,3563000.909)
(40394216.13,3563005)
(40394424,3563015.75)
(40394703.38,3563061.75)
(40394931.63,3563097)
(40395216.03,3563147.5)
(40395463.63,3563146.25)
(40395738.75,3563163)
(40395979.5,3563153.75)
(40396262.63,3563124.5)
(40396321.91,3563120.653)
(40396321.91,3561801.377)
(40396062.16,3561859.25)
(40395587.53,3561941.25)
(40395134.16,3562011.5)
(40394698.41,3562037)
(40394206.91,3562022.75)
(40394017.13,3561996.75)
(40393621.03,3561937.75)
(40393231.03,3561821.25)
(40393084.55,3561788.734)
(40393029.44,3561776.5)
(40392787.41,3561633.093)
(40392680.09,3561569.5)
(40392558.53,3561492.25)
(40392335.91,3561357.75)
(40392002.25,3561114.25)
(40391830.69,3561008)
(40391592.94,3560874.5)
(40391193.22,3560592.5)
(40390765.13,3560300.5)
(40390564.59,3560166)
(40390098.81,3559965.75)
(40389615.03,3559809.5)
(40389331.25,3559691.5)
(40388908.94,3559563.5)
(40388451,3559441.5)
(40388062.34,3559243.25)
(40387893.63,3559170.5)
(40387427.78,3558946.5)
(40387293.56,3558830.5)
(40387100.84,3558659.5)
(40386681.81,3558361.25)
(40386502.34,3558231.25)
(40386131.13,3557972.75)
(40385747.5,3557659.25)
(40385358.03,3557315)
(40385173.66,3557160.75)
(40384996.41,3557018.5)
(40384780.28,3556850.5)
(40384336.47,3556621)
(40384277.52,3556597.71)
(40383870.91,3556437)
(40383535.28,3556243.75)
(40383237.95,3556115.306)
(40383135.97,3556071.25)
(40383005.67,3556011.056)
(40382922.48,3555972.628)
(40382693.04,3555866.635)
(40382677.59,3555859.5)
(40382657.27,3555850.988)
(40382616.89,3555834.07)
(40382556.19,3555808.639)
(40382465.09,3555770.469)
(40382263.47,3555686)
(40381962.15,3555490.101)
(40381676.45,3555343.167)
(40381581.21,3555294.189)
(40381379.83,3555123.419)
(40381285,3555043)
(40381029.81,3554731)
(40380851.38,3554491.5)
(40380588.59,3554122.75)
(40380327.91,3553680)
(40380059.97,3553293.5)
(40379945.75,3553070.047)
(40379842.59,3552868.25)
(40379796.84,3552596)
(40379744.78,3552367.25)
(40379722.02,3552040.953)
(40379707.22,3551828.75)
(40379690.56,3551757.25)
(40379660.04,3551708.138)
(40379523.93,3551776.191)
(40379325.59,3551875.365)
(40379133.03,3551971.645)
(40378890.19,3552093.065)
(40378743.19,3552166.563)
(40378650.45,3552212.933)
(40378618.3,3552229.006)
(40378495.39,3552290.462)
(40378315.88,3552380.216)
(40378171.71,3552452.303)
(40378123.65,3552476.332)
(40378247.97,3552624.25)
(40378587.44,3553038.25)
(40378739.44,3553236.25)
(40378974.63,3553554.437)
(40379053.03,3553660.5)
(40379186.34,3553821.5)
(40379332.17,3554068.313)
(40379449.56,3554267)
(40379548.84,3554456)
(40379783.66,3554866.5)
(40379903.56,3555045.5)
(40380221.5,3555441)
(40380464.72,3555659.491)
(40380580.78,3555763.75)
(40380683.84,3555895.25)
(40380828.33,3556103.846)
(40380964.72,3556300.75)
(40381271.02,3556632.923)
(40381546.97,3556926.5)
(40381902.53,3557309.75)
(40381926.39,3557332.637)
(40381944.21,3557349.735)
(40381985.9,3557389.72)
(40382035.57,3557437.367)
(40382047.99,3557449.279)
(40382051.1,3557452.257)
(40382051.87,3557453.002)
(40382052.13,3557453.25)
(40382161.38,3557541.714)
(40382281.84,3557639.25)
(40382486.31,3557770.75)
(40382749.21,3557976.775)
(40383040.34,3558139.904)
(40383353.56,3558315.419)
(40383520.6,3558409.018)
(40383635.37,3558473.328)
(40383719.45,3558617.007)
(40383786.25,3558731.154)
(40383856.91,3558939.325)
(40383902.75,3559160.864)
(40383897.57,3559460.749)
(40383895.11,3559603.942)
(40383788.13,3560076.101)
(40383673.57,3560581.771)
(40383581.9,3561276.945)
(40383558.98,3561544.32)
(40383536.06,3561765.859)
(40383552.19,3562278.25)
(40383367.28,3562640.062)
(40383300.75,3562770.25)
(40383214.78,3562988.5)
(40383218.19,3563455.75)
(40383148.53,3563701)
(40383023.46,3564009.283)
(40382984.63,3564105)
(40382917.75,3564300.75)
(40382866.75,3564510.75)
(40382814.88,3564752.75)
(40382757.53,3565028.5)
(40382787.69,3565249.5)
(40382826.53,3565503)
(40382958.63,3565680)
(40383207.44,3566020.5)
(40383479.36,3566261.48)
(40383548.78,3566323)
(40383672.36,3566464.433)
(40383723.53,3566523)
(40383881.38,3566665)
(40384100.53,3566841)
(40384487.81,3567137.5)
(40384855.19,3567455)
(40385238.38,3567736.25)
(40385653.19,3568020.25)
(40385813.25,3568126.75)
(40385909.59,3568184.75)
(40386335.16,3568402.25)
(40386773.63,3568696.25)
(40386977.93,3568643.188)
(40387046.03,3568625.5)
(40387120.08,3568645.14)
(40387332.56,3568701.5)
(40387822.91,3568754.75)
(40388029.03,3568791.5)
(40388270.28,3568824)
(40388520.28,3568818.75)
(40385136.56,3558533.466)
(40385531.06,3558738.939)
(40385864.18,3558962.219)
(40386169.11,3559285.007)
(40386297.96,3559403.596)
(40386379.54,3559479.002)
(40386489.09,3559580.25)
(40386635.59,3559714)
(40386823.78,3559883.75)
(40386963.87,3559993.25)
(40387010.56,3560029.75)
(40387223.5,3560199.75)
(40387411.41,3560310.5)
(40387671.25,3560461.75)
(40387857.97,3560575.25)
(40388173.81,3560681.75)
(40388294.34,3560705.25)
(40388580.28,3560742)
(40388675.59,3560754.25)
(40388877.5,3560820.75)
(40389365.09,3560882)
(40389606.72,3560897.25)
(40389746.78,3560908.5)
(40390052.29,3560996.494)
(40390152.13,3561025.25)
(40390371.69,3561155)
(40390534.63,3561247)
(40390909.72,3561450.25)
(40391169.88,3561662.406)
(40391279.13,3561751.5)
(40391650.13,3562080.25)
(40391901.63,3562225)
(40392081.97,3562335)
(40392305.34,3562457.75)
(40392515.63,3562558)
(40392740.66,3562711.75)
(40392836.07,3562819.756)
(40392798.88,3563035.75)
(40392707.58,3563211.061)
(40392683,3563258.25)
(40392572.44,3563446)
(40392471.87,3563656.563)
(40392438.34,3563726.75)
(40392357.59,3563886.75)
(40392293.34,3564162.5)
(40392233.34,3564382.75)
(40392153.69,3564629.75)
(40392003.22,3564964.063)
(40391953.06,3565075.5)
(40391893,3565295)
(40391823.88,3565544.5)
(40391791.03,3565749.25)
(40391688.72,3566038.474)
(40391632.56,3566197.25)
(40391543.35,3566400.107)
(40391448.41,3566616)
(40391217.11,3566855.377)
(40391115.78,3566960.25)
(40390957.91,3567086.75)
(40390619.09,3567337)
(40390466.28,3567425.75)
(40390250.22,3567535.5)
(40390112.19,3567642)
(40390007.04,3567735.648)
(40389926.09,3567807.75)
(40389615.75,3567844.903)
(40389498,3567859)
(40389326.07,3567813.001)
(40389280.28,3567800.75)
(40389100.05,3567843.125)
(40389039.97,3567857.25)
(40388990.26,3567860.778)
(40388789.3,3567875.042)
(40388659.56,3567884.25)
(40388461.63,3567917)
(40388244.05,3567900.424)
(40388084.25,3567888.25)
(40387839.53,3567858.25)
(40387542.31,3567826.188)
(40387410.78,3567812)
(40387268.48,3567775.532)
(40386970.84,3567699.25)
(40386790.06,3567661.75)
(40386365.84,3567522.75)
(40386137,3567437.75)
(40385734.34,3567234.25)
(40385541.56,3567107)
(40385355.44,3567006.75)
(40385231.18,3566931.056)
(40385144.5,3566878.25)
(40385053.85,3566829.221)
(40384715.56,3566646.25)
(40384512.13,3566553.25)
(40384280.88,3566437.5)
(40384110.03,3566326)
(40383922.74,3566240.5)
(40383860.31,3566212)
(40383731.5,3566125.75)
(40383553.56,3565973.5)
(40383425.21,3565709.543)
(40383371.94,3565600)
(40383371.84,3565318.25)
(40383274.28,3564910.026)
(40383245.78,3564790.75)
(40383273.06,3564571.75)
(40383335.34,3564413.428)
(40383428.94,3564175.5)
(40383531.16,3563776.25)
(40383642.44,3563321.75)
(40383810.19,3562821.75)
(40383928.86,3562532.133)
(40383993.25,3562375)
(40384073.5,3562158.25)
(40384227.84,3561709.25)
(40384321.22,3561416.25)
(40384372.13,3561191.5)
(40384462.22,3560927)
(40384494.16,3560711)
(40384531.3,3560487.548)
(40384541.53,3560426)
(40384596.53,3560224.5)
(40384634.59,3559944.5)
(40384634.2,3559831.205)
(40384634.1,3559802.336)
(40384633.85,3559730.603)
(40384633.85,3559703.137)
(40384633.85,3559667.332)
(40384633.85,3559448.881)
(40384633.85,3559389.062)
(40384618.86,3559231.416)
(40384557.12,3558952.631)
(40384561.23,3558714.283)
(40384585.89,3558586.889)
(40384709.17,3558525.248)
(40384873.55,3558471.825)
(40385009.16,3558492.372)]


text:
function simulate_logit(J,T, beta, v)
# --------------------------------------------------%
# Simulate logit with outside option
# --------------------------------------------------%
# J = number of products
# T = number of markets
# beta = coefficient on delta
# v = standard deviation of market-product demand shock

zt = 0.9 .* rand(T,J) .+ 0.05;
xit = randn(T,J).*v;
pt = 2 .*(zt .+ rand(T,J)*0.1).+xit;

denominator = 1;
s = zeros(T,J)
for j = 1:1:J
    denominator = denominator .+ exp.(beta*pt[:,j] .+ xit[:,j]);
end
for j = 1:1:J
    s[:,j] = exp.(beta*pt[:,j] .+ xit[:,j])./denominator;
end
s, pt, zt, xit
end


text:
using ExponentialUtilities, YaoArrayRegister

export TimeEvolution, time_evolve

"""
    TimeEvolution{D, TT, GT} <: PrimitiveBlock{D}

TimeEvolution, where GT is block type. input matrix should be hermitian.

!!!note:
    `TimeEvolution` contructor check hermicity of the input block by default, but sometimes it can be slow. Turn off the check manually by specifying optional parameter `check_hermicity = false`.
"""
mutable struct TimeEvolution{D,Tt,HT<:AbstractBlock{D}} <: PrimitiveBlock{D}
    H::HT
    dt::Tt
    tol::Float64

    function TimeEvolution(
        H::TH,
        dt::Tt;
        tol::Real = 1e-7,
        check_hermicity::Bool = true,
    ) where {D,Tt,TH<:AbstractBlock{D}}
        (check_hermicity && !ishermitian(H)) &&
            error("Time evolution Hamiltonian has to be a Hermitian")
        return new{D,Tt,TH}(H, dt, Float64(tol))
    end
end
nqudits(te::TimeEvolution) = nqudits(te.H)

"""
    TimeEvolution(H, dt[; tol::Real=1e-7])

Create a [`TimeEvolution`](@ref) block with Hamiltonian `H` and time step `dt`. The
`TimeEvolution` block will use Krylove based `expv` to calculate time propagation.

Optional keywords are tolerance `tol` (default is `1e-7`)
`TimeEvolution` block can also be used for
[imaginary time evolution](http://large.stanford.edu/courses/2008/ph372/behroozi2/) if dt is complex.
"""
time_evolve(M::AbstractBlock, dt; kwargs...) = TimeEvolution(M, dt; kwargs...)
time_evolve(dt; kwargs...) = @λ(M -> time_evolve(M, dt; kwargs...))

content(te::TimeEvolution) = te.H
chcontent(te::TimeEvolution, blk::AbstractBlock) = time_evolve(blk, te.dt; tol = te.tol)

function mat(::Type{T}, te::TimeEvolution) where {T}
    return exp(-im * T(te.dt) * Matrix(mat(T, te.H)))
end

struct BlockMap{T,GT<:AbstractBlock} <: AbstractArray{T,2}
    block::GT
    BlockMap(::Type{T}, block::GT) where {T,GT<:AbstractBlock} = new{T,GT}(block)
end

Base.size(bm::BlockMap{T,GT}, i::Int) where {T,D,GT<:AbstractBlock{D}} =
    0 < i <= 2 ? D^nqudits(bm.block) : DimensionMismatch("")
Base.size(bm::BlockMap{T,GT}) where {T,D,GT<:AbstractBlock{D}} = (L = D^nqudits(bm.block); (L, L))
LinearAlgebra.ishermitian(bm::BlockMap) = ishermitian(bm.block)

function LinearAlgebra.mul!(y::AbstractVector, A::BlockMap{T,GT}, x::AbstractVector) where {T,D,GT<:AbstractBlock{D}}
    copyto!(y, x)
    apply!(ArrayReg{D}(y), A.block)
    return y
end

function _apply!(reg::AbstractArrayReg{D,T}, te::TimeEvolution) where {D,T}
    st = state(reg)
    dt = real(te.dt) == 0 ? imag(te.dt) : -im * te.dt
    A = BlockMap(T, te.H)
    @inbounds for j = 1:size(st, 2)
        v = view(st, :, j)
        Ks = arnoldi(A, v; tol = te.tol, ishermitian = true, opnorm = 1.0)
        expv!(v, dt, Ks)
    end
    return reg
end

cache_key(te::TimeEvolution) = (te.dt, cache_key(te.H))

# parametric interface
niparams(::Type{<:TimeEvolution}) = 1
getiparams(x::TimeEvolution) = x.dt
setiparams!(r::TimeEvolution, param::Number) = (r.dt = param; r)
setiparams(r::TimeEvolution, param::Number) =
    TimeEvolution(r.H, param; tol = r.tol, check_hermicity = false)

function Base.:(==)(lhs::TimeEvolution, rhs::TimeEvolution)
    return lhs.H == rhs.H && lhs.dt == rhs.dt
end

function Base.adjoint(te::TimeEvolution)
    return TimeEvolution(te.H, -adjoint(te.dt); tol = te.tol)
end
Base.copy(te::TimeEvolution) = TimeEvolution(te.H, te.dt, tol = te.tol)

function YaoBase.isunitary(te::TimeEvolution)
    iszero(imag(te.dt)) || return false
    return true
end

iparams_range(::TimeEvolution{D,T}) where {D,T} = ((typemin(T), typemax(T)),)


text:
# ### m8.1stan

# m8.1stan is the first model in the Statistical Rethinking book (pp. 249) using Stan.

# Here we will use Turing's NUTS support, which is currently (2018) the originalNUTS by [Hoffman & Gelman]( http://www.stat.columbia.edu/~gelman/research/published/nuts.pdf) and not the one that's in Stan 2.18.2, i.e., Appendix A.5 in: https://arxiv.org/abs/1701.02434

# The TuringModels pkg imports modules such as CSV and DataFrames

using TuringModels, DataFrames

# use reverse mode automatic differentiation
#Turing.setadbackend(:tracker);

# Read in the `rugged` data as a DataFrame

df = CSV.read(joinpath(@__DIR__, "..", "..", "data", "rugged.csv"), DataFrame);

# Show size of the DataFrame (should be 234x51)

size(df)

# Apply log() to each element in rgdppc_2000 column and add it as a new column

df = hcat(df, map(log, df[!, Symbol("rgdppc_2000")]));

# Rename our col x1 => log_gdp

DataFrames.rename!(df, :x1 => :log_gdp);

# Now we need to drop every row where rgdppc_2000 == missing

# When this (https://github.com/JuliaData/DataFrames.jl/pull/1546) hits DataFrame it'll be conceptually easier: i.e., completecases!(d, :rgdppc_2000)

notisnan(e) = !ismissing(e)
dd = df[map(notisnan, df[:, :rgdppc_2000]), :];

# Define the Turing model

@model m8_1stan(y, x₁, x₂) = begin
    σ ~ truncated(Cauchy(0, 2), 0, Inf)
    βR ~ Normal(0, 10)
    βA ~ Normal(0, 10)
    βAR ~ Normal(0, 10)
    α ~ Normal(0, 100)

    for i ∈ 1:length(y)
        y[i] ~ Normal(α + βR * x₁[i] + βA * x₂[i] + βAR * x₁[i] * x₂[i], σ)
    end

    # alternative formulation using broadcasting
    # y .~ Normal.(α .+ βR * x₁ .+ βA * x₂ .+ βAR * x₁ .* x₂, σ)
end;

# Use Turing mcmc

chns = sample(m8_1stan(dd[:, :log_gdp], dd[:, :rugged], dd[:, :cont_africa]),
  NUTS(0.65), 1000)

# Here's the ulam() output from rethinking

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

Empirical Posterior Estimates:
          Mean         SD       Naive SE       MCSE      ESS
    a  9.22360053 0.139119116 0.0021996664 0.0034632816 1000
   bR -0.20196346 0.076106388 0.0012033477 0.0018370185 1000
   bA -1.94430980 0.227080488 0.0035904578 0.0057840746 1000
  bAR  0.39071684 0.131889143 0.0020853505 0.0032749642 1000
sigma  0.95036370 0.052161768 0.0008247500 0.0009204073 1000

Quantiles:
          2.5%       25.0%       50.0%      75.0%        97.5%
    a  8.95307475  9.12719750  9.2237750  9.31974000  9.490234250
   bR -0.35217930 -0.25334425 -0.2012855 -0.15124725 -0.054216855
   bA -2.39010825 -2.09894500 -1.9432550 -1.78643000 -1.513974250
  bAR  0.13496995  0.30095575  0.3916590  0.47887625  0.650244475
sigma  0.85376115  0.91363250  0.9484920  0.98405750  1.058573750
";

# Describe the posterior samples

chns |> display

# End of `08/m8.1t.jl`


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: lakeformation
using AWS.Compat
using AWS.UUIDs

"""
    batch_grant_permissions(entries)
    batch_grant_permissions(entries, params::Dict{String,<:Any})

Batch operation to grant permissions to the principal.

# Arguments
- `entries`: A list of up to 20 entries for resource permissions to be granted by batch
  operation to the principal.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"CatalogId"`: The identifier for the Data Catalog. By default, the account ID. The Data
  Catalog is the persistent metadata store. It contains database definitions, table
  definitions, and other control information to manage your AWS Lake Formation environment.
"""
batch_grant_permissions(Entries; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("BatchGrantPermissions", Dict{String, Any}("Entries"=>Entries); aws_config=aws_config)
batch_grant_permissions(Entries, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("BatchGrantPermissions", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Entries"=>Entries), params)); aws_config=aws_config)

"""
    batch_revoke_permissions(entries)
    batch_revoke_permissions(entries, params::Dict{String,<:Any})

Batch operation to revoke permissions from the principal.

# Arguments
- `entries`: A list of up to 20 entries for resource permissions to be revoked by batch
  operation to the principal.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"CatalogId"`: The identifier for the Data Catalog. By default, the account ID. The Data
  Catalog is the persistent metadata store. It contains database definitions, table
  definitions, and other control information to manage your AWS Lake Formation environment.
"""
batch_revoke_permissions(Entries; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("BatchRevokePermissions", Dict{String, Any}("Entries"=>Entries); aws_config=aws_config)
batch_revoke_permissions(Entries, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("BatchRevokePermissions", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Entries"=>Entries), params)); aws_config=aws_config)

"""
    deregister_resource(resource_arn)
    deregister_resource(resource_arn, params::Dict{String,<:Any})

Deregisters the resource as managed by the Data Catalog. When you deregister a path, Lake
Formation removes the path from the inline policy attached to your service-linked role.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the resource that you want to
  deregister.

"""
deregister_resource(ResourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("DeregisterResource", Dict{String, Any}("ResourceArn"=>ResourceArn); aws_config=aws_config)
deregister_resource(ResourceArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("DeregisterResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceArn"=>ResourceArn), params)); aws_config=aws_config)

"""
    describe_resource(resource_arn)
    describe_resource(resource_arn, params::Dict{String,<:Any})

Retrieves the current data access role for the given resource registered in AWS Lake
Formation.

# Arguments
- `resource_arn`: The resource ARN.

"""
describe_resource(ResourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("DescribeResource", Dict{String, Any}("ResourceArn"=>ResourceArn); aws_config=aws_config)
describe_resource(ResourceArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("DescribeResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceArn"=>ResourceArn), params)); aws_config=aws_config)

"""
    get_data_lake_settings()
    get_data_lake_settings(params::Dict{String,<:Any})

Retrieves the list of the data lake administrators of a Lake Formation-managed data lake.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"CatalogId"`: The identifier for the Data Catalog. By default, the account ID. The Data
  Catalog is the persistent metadata store. It contains database definitions, table
  definitions, and other control information to manage your AWS Lake Formation environment.
"""
get_data_lake_settings(; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("GetDataLakeSettings"; aws_config=aws_config)
get_data_lake_settings(params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("GetDataLakeSettings", params; aws_config=aws_config)

"""
    get_effective_permissions_for_path(resource_arn)
    get_effective_permissions_for_path(resource_arn, params::Dict{String,<:Any})

Returns the Lake Formation permissions for a specified table or database resource located
at a path in Amazon S3. GetEffectivePermissionsForPath will not return databases and tables
if the catalog is encrypted.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the resource for which you want to get
  permissions.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"CatalogId"`: The identifier for the Data Catalog. By default, the account ID. The Data
  Catalog is the persistent metadata store. It contains database definitions, table
  definitions, and other control information to manage your AWS Lake Formation environment.
- `"MaxResults"`: The maximum number of results to return.
- `"NextToken"`: A continuation token, if this is not the first call to retrieve this list.
"""
get_effective_permissions_for_path(ResourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("GetEffectivePermissionsForPath", Dict{String, Any}("ResourceArn"=>ResourceArn); aws_config=aws_config)
get_effective_permissions_for_path(ResourceArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("GetEffectivePermissionsForPath", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceArn"=>ResourceArn), params)); aws_config=aws_config)

"""
    grant_permissions(permissions, principal, resource)
    grant_permissions(permissions, principal, resource, params::Dict{String,<:Any})

Grants permissions to the principal to access metadata in the Data Catalog and data
organized in underlying data storage such as Amazon S3. For information about permissions,
see Security and Access Control to Metadata and Data.

# Arguments
- `permissions`: The permissions granted to the principal on the resource. AWS Lake
  Formation defines privileges to grant and revoke access to metadata in the Data Catalog and
  data organized in underlying data storage such as Amazon S3. AWS Lake Formation requires
  that each principal be authorized to perform a specific task on AWS Lake Formation
  resources.
- `principal`: The principal to be granted the permissions on the resource. Supported
  principals are IAM users or IAM roles, and they are defined by their principal type and
  their ARN. Note that if you define a resource with a particular ARN, then later delete, and
  recreate a resource with that same ARN, the resource maintains the permissions already
  granted.
- `resource`: The resource to which permissions are to be granted. Resources in AWS Lake
  Formation are the Data Catalog, databases, and tables.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"CatalogId"`: The identifier for the Data Catalog. By default, the account ID. The Data
  Catalog is the persistent metadata store. It contains database definitions, table
  definitions, and other control information to manage your AWS Lake Formation environment.
- `"PermissionsWithGrantOption"`: Indicates a list of the granted permissions that the
  principal may pass to other users. These permissions may only be a subset of the
  permissions granted in the Privileges.
"""
grant_permissions(Permissions, Principal, Resource; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("GrantPermissions", Dict{String, Any}("Permissions"=>Permissions, "Principal"=>Principal, "Resource"=>Resource); aws_config=aws_config)
grant_permissions(Permissions, Principal, Resource, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("GrantPermissions", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Permissions"=>Permissions, "Principal"=>Principal, "Resource"=>Resource), params)); aws_config=aws_config)

"""
    list_permissions()
    list_permissions(params::Dict{String,<:Any})

Returns a list of the principal permissions on the resource, filtered by the permissions of
the caller. For example, if you are granted an ALTER permission, you are able to see only
the principal permissions for ALTER. This operation returns only those permissions that
have been explicitly granted. For information about permissions, see Security and Access
Control to Metadata and Data.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"CatalogId"`: The identifier for the Data Catalog. By default, the account ID. The Data
  Catalog is the persistent metadata store. It contains database definitions, table
  definitions, and other control information to manage your AWS Lake Formation environment.
- `"MaxResults"`: The maximum number of results to return.
- `"NextToken"`: A continuation token, if this is not the first call to retrieve this list.
- `"Principal"`: Specifies a principal to filter the permissions returned.
- `"Resource"`: A resource where you will get a list of the principal permissions. This
  operation does not support getting privileges on a table with columns. Instead, call this
  operation on the table, and the operation returns the table and the table w columns.
- `"ResourceType"`: Specifies a resource type to filter the permissions returned.
"""
list_permissions(; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("ListPermissions"; aws_config=aws_config)
list_permissions(params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("ListPermissions", params; aws_config=aws_config)

"""
    list_resources()
    list_resources(params::Dict{String,<:Any})

Lists the resources registered to be managed by the Data Catalog.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"FilterConditionList"`: Any applicable row-level and/or column-level filtering
  conditions for the resources.
- `"MaxResults"`: The maximum number of resource results.
- `"NextToken"`: A continuation token, if this is not the first call to retrieve these
  resources.
"""
list_resources(; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("ListResources"; aws_config=aws_config)
list_resources(params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("ListResources", params; aws_config=aws_config)

"""
    put_data_lake_settings(data_lake_settings)
    put_data_lake_settings(data_lake_settings, params::Dict{String,<:Any})

Sets the list of data lake administrators who have admin privileges on all resources
managed by Lake Formation. For more information on admin privileges, see Granting Lake
Formation Permissions. This API replaces the current list of data lake admins with the new
list being passed. To add an admin, fetch the current list and add the new admin to that
list and pass that list in this API.

# Arguments
- `data_lake_settings`: A structure representing a list of AWS Lake Formation principals
  designated as data lake administrators.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"CatalogId"`: The identifier for the Data Catalog. By default, the account ID. The Data
  Catalog is the persistent metadata store. It contains database definitions, table
  definitions, and other control information to manage your AWS Lake Formation environment.
"""
put_data_lake_settings(DataLakeSettings; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("PutDataLakeSettings", Dict{String, Any}("DataLakeSettings"=>DataLakeSettings); aws_config=aws_config)
put_data_lake_settings(DataLakeSettings, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("PutDataLakeSettings", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("DataLakeSettings"=>DataLakeSettings), params)); aws_config=aws_config)

"""
    register_resource(resource_arn)
    register_resource(resource_arn, params::Dict{String,<:Any})

Registers the resource as managed by the Data Catalog. To add or update data, Lake
Formation needs read/write access to the chosen Amazon S3 path. Choose a role that you know
has permission to do this, or choose the AWSServiceRoleForLakeFormationDataAccess
service-linked role. When you register the first Amazon S3 path, the service-linked role
and a new inline policy are created on your behalf. Lake Formation adds the first path to
the inline policy and attaches it to the service-linked role. When you register subsequent
paths, Lake Formation adds the path to the existing policy. The following request registers
a new location and gives AWS Lake Formation permission to use the service-linked role to
access that location.  ResourceArn = arn:aws:s3:::my-bucket UseServiceLinkedRole = true  If
UseServiceLinkedRole is not set to true, you must provide or set the RoleArn:
arn:aws:iam::12345:role/my-data-access-role

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the resource that you want to register.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"RoleArn"`: The identifier for the role that registers the resource.
- `"UseServiceLinkedRole"`: Designates an AWS Identity and Access Management (IAM)
  service-linked role by registering this role with the Data Catalog. A service-linked role
  is a unique type of IAM role that is linked directly to Lake Formation. For more
  information, see Using Service-Linked Roles for Lake Formation.
"""
register_resource(ResourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("RegisterResource", Dict{String, Any}("ResourceArn"=>ResourceArn); aws_config=aws_config)
register_resource(ResourceArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("RegisterResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceArn"=>ResourceArn), params)); aws_config=aws_config)

"""
    revoke_permissions(permissions, principal, resource)
    revoke_permissions(permissions, principal, resource, params::Dict{String,<:Any})

Revokes permissions to the principal to access metadata in the Data Catalog and data
organized in underlying data storage such as Amazon S3.

# Arguments
- `permissions`: The permissions revoked to the principal on the resource. For information
  about permissions, see Security and Access Control to Metadata and Data.
- `principal`: The principal to be revoked permissions on the resource.
- `resource`: The resource to which permissions are to be revoked.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"CatalogId"`: The identifier for the Data Catalog. By default, the account ID. The Data
  Catalog is the persistent metadata store. It contains database definitions, table
  definitions, and other control information to manage your AWS Lake Formation environment.
- `"PermissionsWithGrantOption"`: Indicates a list of permissions for which to revoke the
  grant option allowing the principal to pass permissions to other principals.
"""
revoke_permissions(Permissions, Principal, Resource; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("RevokePermissions", Dict{String, Any}("Permissions"=>Permissions, "Principal"=>Principal, "Resource"=>Resource); aws_config=aws_config)
revoke_permissions(Permissions, Principal, Resource, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("RevokePermissions", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Permissions"=>Permissions, "Principal"=>Principal, "Resource"=>Resource), params)); aws_config=aws_config)

"""
    update_resource(resource_arn, role_arn)
    update_resource(resource_arn, role_arn, params::Dict{String,<:Any})

Updates the data access role used for vending access to the given (registered) resource in
AWS Lake Formation.

# Arguments
- `resource_arn`: The resource ARN.
- `role_arn`: The new role to use for the given resource registered in AWS Lake Formation.

"""
update_resource(ResourceArn, RoleArn; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("UpdateResource", Dict{String, Any}("ResourceArn"=>ResourceArn, "RoleArn"=>RoleArn); aws_config=aws_config)
update_resource(ResourceArn, RoleArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lakeformation("UpdateResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceArn"=>ResourceArn, "RoleArn"=>RoleArn), params)); aws_config=aws_config)


text:
module DoloBenchmarks
using Dolo
using Dolang
using BenchmarkTools
using JLD

suite = BenchmarkGroup()

include("ea_quest.jl")

if "tune" in ARGS || !isfile("params.jld")
    println("Tuning benchmarks")
    tune!(suite)
    JLD.save("params.jld", "suite", params(suite))
else
    println("loading saved tuning params")
    loadparams!(suite, JLD.load("params.jld", "suite"), :evals, :samples)
end

results = run(suite, verbose=true)

const git_dir = joinpath(dirname(dirname(@__FILE__)), ".git")
const sha = readchomp(`git --git-dir $(git_dir) rev-parse HEAD`)

if !isfile("results.jld")
    close(jldopen("results.jld", "w"))
end

jldopen("results.jld", "r+") do f
    nm = "results_$(sha)"
    JLD.exists(f, nm) && JLD.delete!(f, nm)
    JLD.write(f, nm, results)
end

end  # module


text:
# module for handling the sniper map
module UrbanMaps

export
    UrbanMap,
    locals,
    inbounds

using Colladas # probably don't need this
using Maps

type UrbanMap <: Map
    buildings::Vector{Building}
    nBuildings::Int
    xSize::Int
    ySize::Int
    xl::Float64
    yl::Float64
    gridSize::Int
    # for rectangular maps only
    xCenters::Vector{Float64}
    yCenters::Vector{Float64}
    xLengths::Vector{Float64}
    yLengths::Vector{Float64}

    # TODO: need to deal with normalizing the building vertices for diff x,y sizes
    function UrbanMap(coll::Collada, xSize::Int64, ySize::Int64)
        self = new()

        nBuildings = coll.nObjects

        self.xSize    = xSize
        self.ySize    = ySize
        self.gridSize = xSize * ySize

        self.xl = float(xSize)
        self.yl = float(ySize)

        self.nBuildings = nBuildings

        # assumes a square grid
        self.buildings = [Polygon(xSize * c.points, c.vertices) for c in coll.normalizedObjects]

        return self
    end

    # Generates a map with n rectangular buildings
    # all input coordinates are normalized to unit square
    function UrbanMap(n::Int64, xCenters::Vector{Float64}, yCenters::Vector{Float64},
                      xLengths::Vector{Float64}, yLengths::Vector{Float64},
                      xSize::Int64, ySize::Int64)
        self = new()

        self.xSize      = xSize
        self.ySize      = ySize
        self.gridSize   = xSize * ySize
        self.nBuildings = n

        self.xl = float(xSize)
        self.yl = float(ySize)

        buildings = Array(Polygon, 0)

        # convert centers and lengths to vertices
        xm = [-0.5, 0.5, 0.5, -0.5]
        ym = [0.5, 0.5, -0.5, -0.5]
        verts = [1,2,3,4]
        for i = 1:n
            # make vertices
            xc = xCenters[i]
            yc = yCenters[i]
            xl = xLengths[i]
            yl = yLengths[i]
            points = zeros(4,3)
            for j = 1:4
                points[j,1] = xc + xm[j] * xl
                points[j,2] = yc + ym[j] * yl
            end
            p = Polygon(xSize * points, verts)
            push!(buildings, p)
        end
        self.buildings = buildings

        self.xCenters = xCenters
        self.yCenters = yCenters
        self.xLengths = xLengths
        self.yLengths = yLengths
        return self
    end

    function UrbanMap(t::Symbol, xSize::Int64, ySize::Int64)
        self = new()

        self.xSize = xSize
        self.ySize = ySize

        if t == :cirlce
            self.nBuildings = 1
            self.buildings = [Circular((0.5, 0.5), 0.1)]
        end

        return self
    end

end

function locals(map::UrbanMap)
    return (map.xCenters, map.yCenters, map.xLengths, map.yLengths)
end

# generates a map with rectangular buildings
function generateMaps()

end

function inbounds(map::UrbanMap, x::Vector{Float64})
    if 0.0 <= x[1] <= map.xl && 0.0 <= x[2] <= map.yl
        return true
    end
    return false
end

end # module

max_stars_repo_name: etotheipluspi/ParticleFilters

max_issues_repo_name: etotheipluspi/ParticleFilters

max_forks_repo_name: etotheipluspi/ParticleFilters


text:
using Convertible
using Base.Test
using Compat

abstract type AbstractFrame end

struct Frame1 <: AbstractFrame end
struct Frame2 <: AbstractFrame end
struct Frame3 <: AbstractFrame end
struct Frame4 <: AbstractFrame end

struct State{Frame} end

struct Rotation{From,To} <: Converter{From,To}
    a::Int
    b::Int
    c::Int
end

(::Rotation{From,To})(s::State{From}) where {From,To} = State{To}()

@convertible struct AType end
@convertible struct BType end
@convertible struct CType end
@convertible struct DType end

@convert function BType(a::AType)
    BType()
end
@convert function CType(a::BType)
    CType()
end
@convert function DType(a::BType)
    DType()
end

@convert function Rotation(origin::Frame1, target::Frame2, a, b, c)
    Rotation{origin,target}(a, b, c)
end

@convert function Rotation(origin::Frame2, target::Frame3, a, b, c)
    Rotation{origin,target}(a, b, c)
end

@convert function Rotation(origin::Frame2, target::Frame4, a, b, c)
    Rotation{origin,target}(a, b, c)
end

@testset "Convertible" begin
    @test DType(AType()) == DType()
    @test CType(AType()) == CType()

    s = State{Frame1()}()
    r = Rotation(Frame1(), Frame2(), 1, 2, 3)
    @test r(s) == State{Frame2()}()

    r = Rotation(Frame1(), Frame3(), 1, 2, 3)
    @test r(s) == State{Frame3()}()

    r = Rotation(Frame1(), Frame4(), 1, 2, 3)
    @test r(s) == State{Frame4()}()
end


text:
# Convert cell-centered values to node-centered values by averaging over all
# four neighbors and making use of the periodicity of the solution
function cell2node(cell_centered_data::AbstractArray{Float64})
  # Create temporary data structure to make the averaging algorithm as simple
  # as possible (by using a ghost layer)
  tmp = similar(cell_centered_data, size(cell_centered_data) .+ (2, 2, 0))

  # Fill center with original data
  tmp[2:end-1, 2:end-1, :] .= cell_centered_data

  # Fill sides with opposite data (periodic domain)
  # x-direction
  tmp[1,   2:end-1, :] .= cell_centered_data[end, :, :]
  tmp[end, 2:end-1, :] .= cell_centered_data[1,   :, :]
  # y-direction
  tmp[2:end-1, 1,   :] .= cell_centered_data[:, end, :]
  tmp[2:end-1, end, :] .= cell_centered_data[:, 1,   :]
  # Corners
  tmp[1,   1,   :] = cell_centered_data[end, end, :]
  tmp[end, 1,   :] = cell_centered_data[1,   end, :]
  tmp[1,   end, :] = cell_centered_data[end, 1,   :]
  tmp[end, end, :] = cell_centered_data[1,   1,   :]

  # Create output data structure
  resolution_in, _, n_variables = size(cell_centered_data)
  resolution_out = resolution_in + 1
  node_centered_data = Array{Float64}(undef, resolution_out, resolution_out, n_variables)

  # Obtain node-centered value by averaging over neighboring cell-centered values
  for j in 1:resolution_out
    for i in 1:resolution_out
      node_centered_data[i, j, :] = (tmp[i,   j,   :] +
                                     tmp[i+1, j,   :] +
                                     tmp[i,   j+1, :] +
                                     tmp[i+1, j+1, :]) / 4
    end
  end

  return node_centered_data
end

# Convert 3d unstructured data to 2d slice.
# Additional to the new unstructured data updated coordinates, levels and
# center coordinates are returned.
function unstructured_2d_to_3d(unstructured_data, coordinates, levels,
                               length_level_0, center_level_0, slice_axis,
                               slice_axis_intercept)

  if slice_axis === :x
    slice_axis_dimension = 1
    other_dimensions = [2, 3]
  elseif slice_axis === :y
    slice_axis_dimension = 2
    other_dimensions = [1, 3]
  elseif slice_axis === :z
    slice_axis_dimension = 3
    other_dimensions = [1, 2]
  else
    error("illegal dimension '$slice_axis', supported dimensions are :x, :y, and :z")
  end

  # Limits of domain in slice_axis dimension
  lower_limit = center_level_0[slice_axis_dimension] - length_level_0 / 2
  upper_limit = center_level_0[slice_axis_dimension] + length_level_0 / 2

  if slice_axis_intercept < lower_limit || slice_axis_intercept > upper_limit
    error(string("Slice plane $slice_axis = $slice_axis_intercept outside of domain. ",
        "$slice_axis must be between $lower_limit and $upper_limit"))
  end

  # Extract data shape information
  n_nodes_in, _, _, n_elements, n_variables = size(unstructured_data)

  # Get node coordinates for DG locations on reference element
  nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in)

  # New unstructured data has one dimension less.
  # The redundant element ids are removed later.
  @views new_unstructured_data = similar(unstructured_data[1, ..])

  # Declare new empty arrays to fill in new coordinates and levels
  new_coordinates = Array{Float64}(undef, 2, n_elements)
  new_levels = Array{eltype(levels)}(undef, n_elements)

  # Counter for new element ids
  new_id = 0

  # Save vandermonde matrices in a Dict to prevent redundant generation
  vandermonde_to_2d = Dict()

  # Permute dimensions such that the slice axis dimension is always the
  # third dimension of the array. Below we can always interpolate in the
  # third dimension.
  if slice_axis === :x
    unstructured_data = permutedims(unstructured_data, [2, 3, 1, 4, 5])
  elseif slice_axis === :y
    unstructured_data = permutedims(unstructured_data, [1, 3, 2, 4, 5])
  end

  for element_id in 1:n_elements
    # Distance from center to border of this element (half the length)
    element_length = length_level_0 / 2^levels[element_id]
    min_coordinate = coordinates[:, element_id] .- element_length / 2
    max_coordinate = coordinates[:, element_id] .+ element_length / 2

    # Check if slice plane and current element intersect.
    # The first check uses a "greater but not equal" to only match one cell if the
    # slice plane lies between two cells.
    # The second check is needed if the slice plane is at the upper border of
    # the domain due to this.
    if !((min_coordinate[slice_axis_dimension] <= slice_axis_intercept &&
          max_coordinate[slice_axis_dimension] > slice_axis_intercept) ||
        (slice_axis_intercept == upper_limit &&
          max_coordinate[slice_axis_dimension] == upper_limit))
      # Continue for loop if they don't intersect
      continue
    end

    # This element is of interest
    new_id += 1

    # Add element to new coordinates and levels
    new_coordinates[:, new_id] = coordinates[other_dimensions, element_id]
    new_levels[new_id] = levels[element_id]

    # Construct vandermonde matrix (or load from Dict if possible)
    normalized_intercept =
        (slice_axis_intercept - min_coordinate[slice_axis_dimension]) /
        element_length * 2 - 1

    if haskey(vandermonde_to_2d, normalized_intercept)
      vandermonde = vandermonde_to_2d[normalized_intercept]
    else
      # Generate vandermonde matrix to interpolate values at nodes_in to one value
      vandermonde = polynomial_interpolation_matrix(nodes_in, [normalized_intercept])
      vandermonde_to_2d[normalized_intercept] = vandermonde
    end

    # 1D interpolation to specified slice plane
    # We permuted the dimensions above such that now the dimension in which
    # we will interpolate is always the third one.
    for i in 1:n_nodes_in
      for ii in 1:n_nodes_in
        # Interpolate in the third dimension
        data = unstructured_data[i, ii, :, element_id, :]

        value = interpolate_nodes(permutedims(data), vandermonde, n_variables)
        new_unstructured_data[i, ii, new_id, :] = value[:, 1]
      end
    end
  end

  # Remove redundant element ids
  unstructured_data = new_unstructured_data[:, :, 1:new_id, :]
  new_coordinates = new_coordinates[:, 1:new_id]
  new_levels = new_levels[1:new_id]

  center_level_0 = center_level_0[other_dimensions]

  return unstructured_data, new_coordinates, new_levels, center_level_0
end

# Interpolate unstructured DG data to structured data (cell-centered)
function unstructured2structured(unstructured_data, normalized_coordinates,
                                 levels, resolution, nvisnodes_per_level)
  # Extract data shape information
  n_nodes_in, _, n_elements, n_variables = size(unstructured_data)

  # Get node coordinates for DG locations on reference element
  nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in)

  # Calculate interpolation vandermonde matrices for each level
  max_level = length(nvisnodes_per_level) - 1
  vandermonde_per_level = []
  for l in 0:max_level
    n_nodes_out = nvisnodes_per_level[l + 1]
    dx = 2 / n_nodes_out
    nodes_out = collect(range(-1 + dx/2, 1 - dx/2, length=n_nodes_out))
    push!(vandermonde_per_level, polynomial_interpolation_matrix(nodes_in, nodes_out))
  end

  # For each element, calculate index position at which to insert data in global data structure
  lower_left_index = element2index(normalized_coordinates, levels, resolution, nvisnodes_per_level)

  # Create output data structure
  structured = Array{Float64}(undef, resolution, resolution, n_variables)

  # For each variable, interpolate element data and store to global data structure
  for v in 1:n_variables
    # Reshape data array for use in interpolate_nodes function
    reshaped_data = reshape(unstructured_data[:, :, :, v], 1, n_nodes_in, n_nodes_in, n_elements)

    for element_id in 1:n_elements
      # Extract level for convenience
      level = levels[element_id]

      # Determine target indices
      n_nodes_out = nvisnodes_per_level[level + 1]
      first = lower_left_index[:, element_id]
      last = first .+ (n_nodes_out - 1)

      # Interpolate data
      vandermonde = vandermonde_per_level[level + 1]
      structured[first[1]:last[1], first[2]:last[2], v] .= (
          reshape(interpolate_nodes(reshaped_data[:, :, :, element_id], vandermonde, 1),
                  n_nodes_out, n_nodes_out))
    end
  end

  return structured
end

# For a given normalized element coordinate, return the index of its lower left
# contribution to the global data structure
function element2index(normalized_coordinates::AbstractArray{Float64}, levels::AbstractArray{Int},
                       resolution::Int, nvisnodes_per_level::AbstractArray{Int})
  n_elements = length(levels)

  # First, determine lower left coordinate for all cells
  dx = 2 / resolution
  lower_left_coordinate = Array{Float64}(undef, ndim, n_elements)
  for element_id in 1:n_elements
    nvisnodes = nvisnodes_per_level[levels[element_id] + 1]
    lower_left_coordinate[1, element_id] = (
        normalized_coordinates[1, element_id] - (nvisnodes - 1)/2 * dx)
    lower_left_coordinate[2, element_id] = (
        normalized_coordinates[2, element_id] - (nvisnodes - 1)/2 * dx)
  end

  # Then, convert coordinate to global index
  indices = coordinate2index(lower_left_coordinate, resolution)

  return indices
end

# Find 2D array index for a 2-tuple of normalized, cell-centered coordinates (i.e., in [-1,1])
function coordinate2index(coordinate, resolution::Integer)
  # Calculate 1D normalized coordinates
  dx = 2/resolution
  mesh_coordinates = collect(range(-1 + dx/2, 1 - dx/2, length=resolution))

  # Find index
  id_x = searchsortedfirst.(Ref(mesh_coordinates), coordinate[1, :], lt=(x,y)->x .< y .- dx/2)
  id_y = searchsortedfirst.(Ref(mesh_coordinates), coordinate[2, :], lt=(x,y)->x .< y .- dx/2)
  return transpose(hcat(id_x, id_y))
end

#=function coordinate2index(coordinate, resolution::Integer)=#
#=  # Calculate 1D normalized coordinates=#
#=  dx = 2/resolution=#
#=  mesh_coordinates = collect(range(-1 + dx/2, 1 - dx/2, length=resolution))=#
#==#
#=  # Prepare output storage=#
#=  n_elements = size(coordinate)[2]=#
#=  indices = Array{Int}(undef, ndim, n_elements)=#
#==#
#=  # Find indicex=#
#=  for element_id in 1:n_elements=#
#=    indices[1] = searchsortedfirst(mesh_coordinates, coordinate[1, element_id], lt=(x,y)->x .< y .- dx/2)=#
#=    indices[2] = searchsortedfirst(mesh_coordinates, coordinate[2, element_id], lt=(x,y)->x .< y .- dx/2)=#
#=  end=#
#==#
#=  return hcat(id_x, id_y)=#
#=end=#

# Find 2D array index for a 2-tuple of normalized, cell-centered coordinates (i.e., in [-1,1])
#=function coordinate2index(coordinate, resolution::Integer)=#
#=  # Calculate 1D normalized coordinates=#
#=  dx = 2/resolution=#
#=  mesh_coordinates = collect(range(-1 + dx/2, 1 - dx/2, length=resolution))=#
#==#
#=  # Build mesh for nearest neighbor search=#
#=  mesh = Array{Float64}(undef, ndim, resolution, resolution)=#
#=  for j in 1:resolution=#
#=    for i in 1:resolution=#
#=      mesh[1, i, j] = mesh_coordinates[i]=#
#=      mesh[2, i, j] = mesh_coordinates[j]=#
#=    end=#
#=  end=#
#==#
#=  # Create tree=#
#=  tree = KDTree(mesh, leafsize=10)=#
#==#
#=  # Find index in a nearest-neighbor search=#
#=  index, _ = knn(tree, coordinate, 1)=#
#=end=#

function interpolate_data(data_in::AbstractArray, n_nodes_in::Integer, n_nodes_out::Integer)
  # Get node coordinates for input and output locations on reference element
  nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in)
  dx = 2/n_nodes_out
  #=nodes_out = collect(range(-1 + dx/2, 1 - dx/2, length=n_nodes_out))=#
  nodes_out = collect(range(-1, 1, length=n_nodes_out))

  # Get interpolation matrix
  vandermonde = polynomial_interpolation_matrix(nodes_in, nodes_out)

  # Create output data structure
  n_elements = div(size(data_in, 1), n_nodes_in^ndim)
  n_variables = size(data_in, 2)
  data_out = Array{eltype(data_in)}(undef, n_nodes_out, n_nodes_out, n_elements, n_variables)

  for n in 1:1
  # Interpolate each variable separately
  for v = 1:n_variables
    # Reshape data to fit expected format for interpolation function
    # FIXME: this "reshape here, reshape later" funny business should be implemented properly
    reshaped = reshape(data_in[:, v], 1, n_nodes_in, n_nodes_in, n_elements)

    # Interpolate data for each cell
    for element_id = 1:1#n_elements
      data_out[:, :, element_id, v] = interpolate_nodes(reshaped[:, :, :, element_id],
                                                        vandermonde, 1)
    end
  end
  end

  return reshape(data_out, n_nodes_out^ndim * n_elements, n_variables)
end

function calc_vertices(coordinates::AbstractArray{Float64, 2},
                       levels::AbstractArray{Int}, length_level_0::Float64)
  @assert ndim == 2 "Algorithm currently only works in 2D"

  # Initialize output arrays
  n_elements = length(levels)
  x = Array{Float64, 2}(undef, 2^ndim+1, n_elements)
  y = Array{Float64, 2}(undef, 2^ndim+1, n_elements)

  # Calculate vertices for all coordinates at once
  for element_id in 1:n_elements
    length = length_level_0 / 2^levels[element_id]
    x[1, element_id] = coordinates[1, element_id] - 1/2 * length
    x[2, element_id] = coordinates[1, element_id] + 1/2 * length
    x[3, element_id] = coordinates[1, element_id] + 1/2 * length
    x[4, element_id] = coordinates[1, element_id] - 1/2 * length
    x[5, element_id] = coordinates[1, element_id] - 1/2 * length

    y[1, element_id] = coordinates[2, element_id] - 1/2 * length
    y[2, element_id] = coordinates[2, element_id] - 1/2 * length
    y[3, element_id] = coordinates[2, element_id] + 1/2 * length
    y[4, element_id] = coordinates[2, element_id] + 1/2 * length
    y[5, element_id] = coordinates[2, element_id] - 1/2 * length
  end

  return x, y
end


text:
using Bijectors
using Random, ProgressBars
using Distributions, DistributionsAD
using NNlib, SpecialFunctions

"""
    invlink_parameters(θ::AbstractVector, K::Int, D::Int)

Map a vector of reals (θ) to the constrained space used for the log density problem.
The function splits the θ vector into the individual paramters and links them back to their constrained space if necessary.
"""
function invlink_parameters(θ::AbstractVector, K::Int, D::Int; pν = InverseGamma(2,3), pb_h = Uniform(-π, π))
    s = 1
    e = (s-1) + 1
    ν = invlink(pν, θ[s])

    s = e+1
    e = (s-1) + (K*D)
    w_h = reshape(θ[s:e], D, K)

    s = e+1
    e = (s-1) + (K)
    b_h = reshape(invlink.(pb_h, θ[s:e]), 1, K)

    s = e+1
    e = (s-1) + K
    w_o = θ[s:e]

    s = e+1
    b_o = θ[end]

    return (ν, w_h, b_h, w_o, b_o)
end

function invlink_parameters(θ::AbstractVector, K::Int; pν = InverseGamma(2,3), pb = Uniform(-π, π))
    s = 1
    e = (s-1) + 1
    ν = invlink(pν, θ[s])

    s = e+1
    e = (s-1) + K
    w_h = θ[s:e]

    s = e+1
    e = (s-1) + K
    b_h = reshape(invlink.(pb, θ[s:e]), 1, K)

    s = e+1
    e = (s-1) + K
    w_o = θ[s:e]

    return (ν, w_h, b_h, w_o)
end

"""
    init(K::Int, D::Int; seed = 2021)

Return initial parameters by randomly sampling them.
"""
function init(kernel::Kernel, K::Int, D::Int; seed = 2021)

    Random.seed!(seed)

    νℝ₊ = rand()
    w_h = init(kernel, D*K)
    bℝ₊_h = rand(K)
    w_o = randn(K)
    b_o = randn(K)

    θ = vcat(νℝ₊, w_h, bℝ₊_h, w_o, b_o)

    return θ
end

function init(K::Int; seed = 2021)

    Random.seed!(seed)

    νℝ₊ = rand()
    w_h = vec(randn(K))
    bℝ₊_h = vec(rand(K))
    w_o = randn(K)

    θ = vcat(νℝ₊, w_h, bℝ₊_h, w_o)

    return θ
end

# helper to compute pmf of Bernoulli
bernoulli(pj,yj) = yj == zero(yj) ? one(pj) - pj : pj

"""
    predict(Xnew, samples)

Return the mean function value for each Xnew.
"""
function predict_banana(Xnew, samples, kernel, K)

    Nnew, D = size(Xnew)
    pw_h, pw_o, pb_h, pb_o = getpriors(kernel, K, D)
    pν = InverseGamma(2,3) # prior on the lengthscale

    issincos = hasSinCosActivation(kernel)

    p = zeros(Nnew, length(samples))
    for (i, θ) in ProgressBar(enumerate(samples))
        ν, w_h, b_h, w_o, b_o = invlink_parameters(θ, K, D; pν = pν, pb_h = pb_h)

        z = issincos ? σ.(Ref(kernel), Xnew*(w_h*sqrt(ν))) : σ.(Ref(kernel), Xnew*(w_h*sqrt(ν)) .+ b_h)
        p[:,i] = logistic.(z * w_o .+ b_o)
    end

    return mean(p, dims=2), var(p, dims=2)
end

function predict_regression(xnew, samples, kernel, K)
    Nnew = length(xnew)
    pw_h, pw_o, pb = getpriors(kernel, K)
    pν = InverseGamma(2,3) # prior on the lengthscale

    p = zeros(Nnew, length(samples))
    for (i, θ) in enumerate(samples)
        ν, w_h, b_h, w_o = invlink_parameters(θ, K; pν = pν, pb = pb)
        z = σ.(Ref(kernel), xnew*(w_h') .+ b_h)
        p[:,i] = z * w_o
    end

    return mean(p, dims=2), var(p, dims=2)
end

max_stars_repo_name: AaltoML/PeriodicBNN

max_issues_repo_name: AaltoML/PeriodicBNN

max_forks_repo_name: AaltoML/PeriodicBNN


text:
# Test that basic function works and that it doesn change the mean
@test all(findcenteredmean(m, only_positive = true) .>= 0)
@test all(getmean(m) .== means)

# Try different centers and changing the mean
@test isa(findcenteredmean(m, center = :analytic), Vector)
@test all(findcenteredmean(m, center = :analytic, only_positive = true, update_mean = true) .>= 0)
@test !all(getmean(m) .== means)

# Test invalid input
@test_throws ErrorException findcenteredmean(m, center = :Analytic)

# Reset the mean back to correct value
setmean(m, means)
@test all(getmean(m) .== means)