text:
using LightGraphs
using Test

include("../../src/utils/filters.jl")
include("../../src/utils/maximalcliques.jl")
include("../../src/Generate/filteredges.jl")

@testset "filteredges" begin
    @testset "degree" begin
        graph = path_graph(10)
        add_edge!(graph, 1, 4)
        add_edge!(graph, 3, 8)
        add_edge!(graph, 4, 8)
        add_edge!(graph, 5, 7)
        add_edge!(graph, 6, 9)
        v = vertices(graph)
        @test filter_vertices_degree_max(vertices=v, graph=graph, degree=2) == [1, 2, 10]
        @test filter_vertices_degree_max(vertices=v, graph=graph, degree=0) == []
        @test filter_vertices_degree_min(vertices=v, graph=graph, degree=3) == [3, 4, 5, 6, 7, 8, 9]
        @test filter_vertices_degree_min(vertices=v, graph=graph, degree=5) == []
        @test filter_vertices_degree(vertices=v, graph=graph, degrees=[3]) == [3, 5, 6, 7, 9]
        @test filter_vertices_degree(vertices=v, graph=graph, degrees=[2, 3]) == [1, 2, 3, 5, 6, 7, 9]
        @test filter_vertices_degree(vertices=v, graph=graph, degrees=[0, 3]) == [3, 5, 6, 7, 9]
        @test filter_vertices_degree(vertices=v, graph=graph, degrees=[0]) == []
    end
    @testset "biggest clique" begin
        graph = path_graph(11)
        add_edge!(graph, 1, 10)
        add_edge!(graph, 1, 3)
        add_edge!(graph, 1, 9)
        add_edge!(graph, 3, 5)
        add_edge!(graph, 5, 7)
        add_edge!(graph, 7, 9)
        add_edge!(graph, 3, 9)
        add_edge!(graph, 3, 7)
        v = vertices(graph)
        cliques = get_maximal_cliques(graph)
        @test filter_vertices_size_biggest_clique_max(vertices=v, cliques=cliques, size=3) == 1:11
        @test filter_vertices_size_biggest_clique_max(vertices=v, cliques=cliques, size=2) == [11]
        @test filter_vertices_size_biggest_clique_max(vertices=v, cliques=cliques, size=1) == []
        @test filter_vertices_size_biggest_clique_min(vertices=v, cliques=cliques, size=2) == 1:11
        @test filter_vertices_size_biggest_clique_min(vertices=v, cliques=cliques, size=3) == 1:10
        @test filter_vertices_size_biggest_clique_min(vertices=v, cliques=cliques, size=4) == []
        @test filter_vertices_size_biggest_clique(vertices=v, cliques=cliques, sizes=[2]) == [11]
        @test filter_vertices_size_biggest_clique(vertices=v, cliques=cliques, sizes=[3]) == 1:10
        @test filter_vertices_size_biggest_clique(vertices=v, cliques=cliques, sizes=[4]) == []
    end
    @testset "smallest clique" begin
        graph = path_graph(11)
        add_edge!(graph, 1, 10)
        add_edge!(graph, 1, 3)
        add_edge!(graph, 1, 9)
        add_edge!(graph, 3, 5)
        add_edge!(graph, 5, 7)
        add_edge!(graph, 7, 9)
        add_edge!(graph, 3, 9)
        add_edge!(graph, 3, 7)
        v = vertices(graph)
        cliques = get_maximal_cliques(graph)
        @test filter_vertices_size_smallest_clique_max(vertices=v, cliques=cliques, size=3) == 1:11
        @test filter_vertices_size_smallest_clique_max(vertices=v, cliques=cliques, size=2) == [10, 11]
        @test filter_vertices_size_smallest_clique_max(vertices=v, cliques=cliques, size=1) == []
        @test filter_vertices_size_smallest_clique_min(vertices=v, cliques=cliques, size=2) == 1:11
        @test filter_vertices_size_smallest_clique_min(vertices=v, cliques=cliques, size=3) == 1:9
        @test filter_vertices_size_smallest_clique_min(vertices=v, cliques=cliques, size=4) == []
        @test filter_vertices_size_smallest_clique(vertices=v, cliques=cliques, sizes=[2]) == [10, 11]
        @test filter_vertices_size_smallest_clique(vertices=v, cliques=cliques, sizes=[3]) == 1:9
        @test filter_vertices_size_smallest_clique(vertices=v, cliques=cliques, sizes=[4]) == []
    end
    @testset "distance" begin
        graph = path_graph(10)
        add_edge!(graph, 1, 4)
        add_edge!(graph, 3, 8)
        add_edge!(graph, 4, 8)
        add_edge!(graph, 5, 7)
        add_edge!(graph, 6, 9)
        v = vertices(graph)
        @test filter_vertices_distance_max(vertices=v, src=1, graph=graph, dist=4) == 1:10
        @test filter_vertices_distance_max(vertices=v, src=1, graph=graph, dist=1) == [1, 2, 4]
        @test filter_vertices_distance_max(vertices=v, src=1, graph=graph, dist=-1) == []
        @test filter_vertices_distance_min(vertices=v, src=1, graph=graph, dist=1) == 2:10
        @test filter_vertices_distance_min(vertices=v, src=1, graph=graph, dist=3) == [6, 7, 9, 10]
        @test filter_vertices_distance_min(vertices=v, src=1, graph=graph, dist=6) == []
        @test filter_vertices_distance(vertices=v, src=1, graph=graph, dists=[3]) == [6, 7, 9]
        @test filter_vertices_distance(vertices=v, src=1, graph=graph, dists=[4]) == [10]
        @test filter_vertices_distance(vertices=v, src=1, graph=graph, dists=[5]) == []
    end
    @testset "same clique" begin
        graph = path_graph(11)
        add_edge!(graph, 1, 10)
        add_edge!(graph, 1, 3)
        add_edge!(graph, 1, 9)
        add_edge!(graph, 3, 5)
        add_edge!(graph, 5, 7)
        add_edge!(graph, 7, 9)
        add_edge!(graph, 3, 9)
        add_edge!(graph, 3, 7)
        v = vertices(graph)
        cliques = get_maximal_cliques(graph)
        @test filter_vertices_same_clique(vertices=v, src=1, cliques=cliques) == [1, 2, 3, 9, 10]
        @test filter_vertices_same_clique(vertices=v, src=11, cliques=cliques) == [10, 11]
        @test filter_vertices_same_clique(vertices=v, src=12, cliques=cliques) == []
    end
    @testset "has/hasn't edge" begin
        graph = path_graph(10)
        add_vertex!(graph)
        add_edge!(graph, 1, 4)
        add_edge!(graph, 3, 8)
        add_edge!(graph, 4, 8)
        add_edge!(graph, 5, 7)
        add_edge!(graph, 6, 9)
        v = vertices(graph)
        @test filter_vertices_no_edge(vertices=v, src=1, choleskygraph=graph) == [1, 3, 5, 6, 7, 8, 9, 10, 11]
        @test filter_vertices_no_edge(vertices=v, src=11, choleskygraph=graph) == 1:11
        @test filter_vertices_has_edge(vertices=v, src=1, choleskygraph=graph) == [2, 4]
        @test filter_vertices_has_edge(vertices=v, src=11, choleskygraph=graph) == []
    end
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
module TensorRefinement

using Reexport

include("Auxiliary.jl")
@reexport using .Auxiliary

include("TensorTrain.jl")
@reexport using .TensorTrain

include("Exponential.jl")
@reexport using .Exponential

include("FEM.jl")
@reexport using .FEM

end

ext: jl


text:
"""
    ABM4OSN

Provide agent-based simulation environment for Online Social Networks.

# Examples
```julia-repl
julia>using ABM4OSN
```
"""
module ABM4OSN

    # requirements
    using LightGraphs
    using DataFrames
    using Statistics
    using Random
    using StatsBase
    using JLD
    using JLD2
    using JSON
    using CSV
    using GraphIO
    import ParserCombinator

    # contain structs
    include("post.jl")
    include("agent.jl")
    include("config.jl")
    include("network.jl")
    include("simulation.jl")

    # contain functions
    include("add_input.jl")
    include("convert_results.jl")
    include("create_agents.jl")
    include("dislike.jl")
    include("drop_input.jl")
    include("evaluate_results.jl")
    include("generate_agent_props.jl")
    include("like.jl")
    include("log_network.jl")
    include("publish_post.jl")
    include("share.jl")
    include("set_inactive.jl")
    include("update_check_regularity.jl")
    include("update_feed.jl")
    include("update_input.jl")
    include("update_opinion.jl")
    include("update_perceiv_publ_opinion.jl")

    # exports for user interaction
    export Config
    export Simulation
    export cfg_ag
    export cfg_feed
    export cfg_net
    export cfg_ot
    export cfg_sim
    export cfg_mech
    export run!
    export run_batch
    export rerun_single!
    export evaluate_results
    export convert_results

end # module


text:
abstract type AbstractStorageFormulation <: AbstractDeviceFormulation end

struct BookKeeping <: AbstractStorageFormulation end

struct BookKeepingwReservation <: AbstractStorageFormulation end

#################################################Storage Variables#################################

function active_power_variables(canonical_model::CanonicalModel,
                                devices::IS.FlattenIteratorWrapper{St}) where {St<:PSY.Storage}

    add_variable(canonical_model,
                 devices,
                 Symbol("Psin_$(St)"),
                 false,
                 :nodal_balance_active,
                 -1.0;
                 lb_value = d -> 0.0,)
    add_variable(canonical_model,
                 devices,
                 Symbol("Psout_$(St)"),
                 false,
                 :nodal_balance_active;
                 lb_value = d -> 0.0,)

    return

end

function reactive_power_variables(canonical_model::CanonicalModel,
                                  devices::IS.FlattenIteratorWrapper{St}) where {St<:PSY.Storage}
    add_variable(canonical_model,
                 devices,
                 Symbol("Qst_$(St)"),
                 false,
                 :nodal_balance_reactive)

    return

end

function energy_storage_variables(canonical_model::CanonicalModel,
                                  devices::IS.FlattenIteratorWrapper{St}) where St<:PSY.Storage

    add_variable(canonical_model,
                 devices,
                 Symbol("Est_$(St)"),
                 false;
                 lb_value = d -> 0.0,)

    return

end

function storage_reservation_variables(canonical_model::CanonicalModel,
                                       devices::IS.FlattenIteratorWrapper{St}) where St<:PSY.Storage

    add_variable(canonical_model,
                 devices,
                 Symbol("Rst_$(St)"),
                 true)

    return

end

###################################################### output power constraints#################################

function active_power_constraints(canonical_model::CanonicalModel,
                                  devices::IS.FlattenIteratorWrapper{St},
                                  ::Type{BookKeeping},
                                  ::Type{S}) where {St<:PSY.Storage,
                                                                      S<:PM.AbstractPowerModel}

    range_data_in = [(PSY.get_name(s), PSY.get_inputactivepowerlimits(s)) for s in devices]
    range_data_out = [(PSY.get_name(s), PSY.get_outputactivepowerlimits(s)) for s in devices]

    device_range(canonical_model,
                 range_data_in,
                 Symbol("inputpower_range_$(St)"),
                 Symbol("Psin_$(St)"))

    device_range(canonical_model,
                range_data_out,
                Symbol("outputpower_range_$(St)"),
                Symbol("Psout_$(St)"))

    return

end

function active_power_constraints(canonical_model::CanonicalModel,
                                  devices::IS.FlattenIteratorWrapper{St},
                                  ::Type{BookKeepingwReservation},
                                  ::Type{S}) where {St<:PSY.Storage,
                                                                      S<:PM.AbstractPowerModel}

    range_data_in = [(PSY.get_name(s), PSY.get_inputactivepowerlimits(s)) for s in devices]
    range_data_out = [(PSY.get_name(s), PSY.get_outputactivepowerlimits(s)) for s in devices]

    reserve_device_semicontinuousrange(canonical_model,
                                       range_data_in,
                                       Symbol("inputpower_range_$(St)"),
                                       Symbol("Psin_$(St)"),
                                       Symbol("Rst_$(St)"))

    reserve_device_semicontinuousrange(canonical_model,
                                       range_data_out,
                                       Symbol("outputpower_range_$(St)"),
                                       Symbol("Psout_$(St)"),
                                       Symbol("Rst_$(St)"))

    return

end

"""
This function adds the reactive  power limits of generators when there are CommitmentVariables
"""
function reactive_power_constraints(canonical_model::CanonicalModel,
                                   devices::IS.FlattenIteratorWrapper{St},
                                   ::Type{D},
                                   ::Type{S}) where {St<:PSY.Storage,
                                                                       D<:AbstractStorageFormulation,
                                                                       S<:PM.AbstractPowerModel}

    range_data = [(PSY.get_name(s), PSY.get_reactivepowerlimits(s)) for s in devices]

    device_range(canonical_model,
                 range_data,
                 Symbol("reactiverange_$(St)"),
                 Symbol("Qst_$(St)"))

    return

end

########################## Make initial Conditions for a Model #############################
function initial_conditions!(canonical_model::CanonicalModel,
                            devices::IS.FlattenIteratorWrapper{St},
                            ::Type{D}) where {St<:PSY.Storage,
                                                                D<:AbstractStorageFormulation}

    storage_energy_init(canonical_model, devices)

return

end

###################################################### Energy Capacity constraints##########

function energy_capacity_constraints(canonical_model::CanonicalModel,
                                    devices::IS.FlattenIteratorWrapper{St},
                                    ::Type{D},
                                    ::Type{S}) where {St<:PSY.Storage,
                                                                        D<:AbstractStorageFormulation,
                                                                        S<:PM.AbstractPowerModel}

    range_data = [(PSY.get_name(s), PSY.get_capacity(s)) for s in devices]

    device_range(canonical_model,
                 range_data,
                 Symbol("energy_capacity_$(St)"),
                 Symbol("Est_$(St)"))
    return

end

###################################################### book keeping constraints ############

function make_efficiency_data(devices::IS.FlattenIteratorWrapper{St}) where {St<:PSY.Storage}

    names = Vector{String}(undef, length(devices))
    in_out = Vector{InOut}(undef, length(devices))

    for (ix, d) in enumerate(devices)
        names[ix] = PSY.get_name(d)
        in_out[ix] = PSY.get_efficiency(d)
    end

    return names, in_out

end

function energy_balance_constraint(canonical_model::CanonicalModel,
                                   devices::IS.FlattenIteratorWrapper{St},
                                   ::Type{D},
                                   ::Type{S}) where {St<:PSY.Storage,
                                                            D<:AbstractStorageFormulation,
                                                            S<:PM.AbstractPowerModel}

    key = ICKey(DeviceEnergy, St)

    if !(key in keys(canonical_model.initial_conditions))
        error("Initial Conditions for $(St) Energy Constraints not in the model")
    end

    efficiency_data = make_efficiency_data(devices)

    energy_balance(canonical_model,
                   canonical_model.initial_conditions[key],
                   efficiency_data,
                   Symbol("energy_balance_$(St)"),
                   (Symbol("Psout_$(St)"), Symbol("Psin_$(St)"), Symbol("Est_$(St)")))

    return

end


text:
module KernelFunctions

export kernelmatrix, kernelmatrix!, kerneldiagmatrix, kerneldiagmatrix!, kappa, kernelpdmat # Main matrix functions
export params, duplicate, set! # Helpers

export Kernel
export ConstantKernel, WhiteKernel, ZeroKernel
export SqExponentialKernel, ExponentialKernel, GammaExponentialKernel
export ExponentiatedKernel
export MaternKernel, Matern32Kernel, Matern52Kernel
export LinearKernel, PolynomialKernel
export RationalQuadraticKernel, GammaRationalQuadraticKernel
export KernelSum, KernelProduct

export Transform, SelectTransform, ChainTransform, ScaleTransform, LowRankTransform, IdentityTransform, FunctionTransform

using Compat
using Distances, LinearAlgebra
using SpecialFunctions: lgamma, besselk
using Zygote: @adjoint
using StatsFuns: logtwo
using PDMats: PDMat

const defaultobs = 2

"""
Abstract type defining a slice-wise transformation on an input matrix
"""
abstract type Transform end
abstract type Kernel{T,Tr<:Transform} end

include("utils.jl")
include("distances/dotproduct.jl")
include("distances/delta.jl")
include("transform/transform.jl")
kernels = ["exponential","matern","polynomial","constant","rationalquad","exponentiated"]
for k in kernels
    include(joinpath("kernels",k*".jl"))
end
include("matrix/kernelmatrix.jl")
include("matrix/kernelpdmat.jl")
include("kernels/kernelsum.jl")
include("kernels/kernelproduct.jl")

include("generic.jl")

include("zygote_adjoints.jl")

end


text:
using Dash
using PlotlyJS
using DashBootstrapComponents

#
# NOT USED! ATTEMPT TO USE ALL-PLOTLY DASHBOARD
#

#, DashHtmlComponents, DashCoreComponents

using Markdown
using DataFrames

include( "runner_libs.jl")
include( "dash_libs.jl")

const FORM_EXTRA =Dict(
	"border-bottom"=>"1px dashed #aaaaaa",
	"margin-bottom"=>"5px",
	"margin-top"=>"5px")

const PREAMBLE = """
"""

const INFO = """

* Created with [Julia](https://julialang.org/) | [Dash](https://dash-julia.plotly.com/) | [Plotly](https://plotly.com/julia/) | [Inequality an Poverty Measures](https://github.com/grahamstark/PovertyAndInequalityMeasures.jl);
* Part of the [Scottish Tax Benefit Model](https://github.com/grahamstark/ScottishTaxBenefitModel.jl);
* Open Source software released under the [MIT Licence](https://github.com/grahamstark/Visualisations.jl/blob/main/LICENSE). [Source Code](https://github.com/grahamstark/Visualisations.jl).
"""

"""
Create the block of sliders and radios on the LHS
"""
function get_input_block()
	return dbc_form([
		dbc_row([
			dbc_col(
				dbc_label("Basic Rate"; html_for="basic_rate"), width=3
			),
			dbc_col(
				dcc_slider(
					id = "basic_rate",
					min = 0,
					max = 100,
					marks = Dict([Symbol("$v") => Symbol("$v") for v in 0:10:100]),
					value = 10.0,
					step = 1
				)), # col
		]),
		dbc_row([
			dbc_col([
				dbc_button(
					id = "submit-button",
					class_name="primary",
					color = "primary",
					name = "Run",
					value = "Run",
					children = "submit"
					)
			]) # col
		]) # row
	]) # form
end

app = dash(external_stylesheets=[dbc_themes.UNITED])
# BOOTSTRAP|SIMPLEX|MINTY|COSMO|SANDSTONE|UNITED|SLATE|SOLAR|UNITED|
app.layout = dbc_container(fluid=true, className="p-5") do
	html_title( "You are The Finance Minister")
	html_h1("You are The Finance Minister"),
	dbc_row([
		dbc_col( dcc_markdown( PREAMBLE ), width=10)
	]), # row
	dbc_row([
    	dbc_col(get_input_block(), width=4),
		dbc_col([
			dcc_loading(
				id="model_running",
				type="default",
				children = [
					html_div(
						id="loading-output-1"
						#=
						children=(
							dcc_graph( id = "bc-1", figure=f1 )
						) # child graph
						=#
					) # div
				] # children
			) # loading
		]) # col
	]), # row

	dbc_row([
		dbc_col([
			dcc_graph( id = "bc-1" )])
	]),
	dbc_row([
		dbc_col( dcc_markdown( INFO ), width=10)
	]) # row
end # layout

function do_output( br )
	br /= 100.0
	sys = deepcopy( BASE_STATE.sys )
	println("sys.it.non_savings_rates[2] $(sys.it.non_savings_rates[2])")
	incr = br-sys.it.non_savings_rates[2]
	sys.it.non_savings_rates[1:3] .+= incr
	results = do_run( sys )
	gls1 = gain_lose_table(
		results.gain_lose )
	gls2 = gain_lose_table(
		results.gain_lose )
	gls3 = gain_lose_table(
		results.gain_lose )
	lorenz = draw_lorenz(
		BASE_STATE.summary.deciles[1][:,2],
		results.summary.deciles[1][:,2] )
	gbd = drawDeciles(
		results.summary.deciles[1][:,3],
		BASE_STATE.summary.deciles[1][:,3] )

	fig = make_subplots(
		rows=3,
		cols=3,
		column_widths=[0.3, 0.3, 0.4],
		row_heights=[0.33, 0.33, 0.33],
		specs=[
			Spec(kind="xy")  Spec(kind="table", rowspan=1, colspan=1)  Spec(kind="xy");
			Spec(kind="xy")  Spec(kind="table", rowspan=1, colspan=1)  Spec(kind="xy");
			Spec(kind="xy")  Spec(kind="xy")  Spec(kind="table", rowspan=1, colspan=1)]
	)
	for row in 1:3
		for col in 1:3
			if row > 1 && (row == col)
				add_trace!( fig, gls1, row=row, col=col)
			elseif row == 1 && col == 2
				add_trace!( fig, gls1, row=row, col=col)
			elseif rand(1:2) == 1
				add_trace!( fig, gbd, row=row, col=col )
			else
				for i in lorenz
					add_trace!( fig, i, row=row, col=col )
				end
			end
		end
	end
		# specs=[
		#		Spec(kind="table") Spec(kind="xy")
		#
		#		missing Spec(kind= "scene")
		#
		#	]

		#

	return fig

end

callback!(
    app,
    Output("model_running",  "children"),
	Output("bc-1", "figure"),
	Input("submit-button", "n_clicks"),
	State( "basic_rate", "value")) do n_clicks, basic_rate
	println( "n_clicks = $n_clicks")
	if ! isnothing( n_clicks )
		return [nothing,do_output( basic_rate )]
	end
	[nothing,do_output(20)]
end

run_server(app, "0.0.0.0", 8051; debug=true )


text:
@testset "phreds" begin

    @testset "phred_to_p" begin
        out = phred_to_p(Int8(20))
        @test out == 0.01
    end

    @testset "p_to_phred" begin
        out = p_to_phred(Float64(0.01))
        @test out == 20
    end

end


text:
import Base: run

using ReinforcementLearningEnvironments: reset!

export run

"""
    run(agent::AbstractAgent, env::AbstractEnv, stop_condition; hook = EmptyHook())

See also [Stop Conditions](@ref), [Hooks](@ref)
"""
function run(agent::AbstractAgent, env::AbstractEnv, stop_condition; hook = EmptyHook())

    reset!(env)
    obs = reset(observe(env); terminal = true)
    hook(PRE_EPISODE_STAGE, agent, env, obs)

    while true
        action = agent(obs)
        hook(PRE_ACT_STAGE, agent, env, obs => action)

        # async here?
        update!(agent, obs => action)
        env(action)
        obs = observe(env)
        hook(POST_ACT_STAGE, agent, env, action => obs)

        if stop_condition(agent, env, obs)
            if get_terminal(obs)
                hook(POST_EPISODE_STAGE, agent, env, obs)
            end
            update!(agent, obs => agent(obs))  # push reward terminal into buffer, state and action are not that important
            break
        end

        if get_terminal(obs)
            hook(POST_EPISODE_STAGE, agent, env, obs)
            r = get_reward(obs)  # !!! deepcopy?
            reset!(env)
            obs = reset(observe(env); reward = r, terminal = true)
            hook(PRE_EPISODE_STAGE, agent, env, obs)
        end
    end
    hook
end

function run(
    agents::Tuple{Vararg{<:AbstractAgent}},
    env::AbstractEnv,
    stop_condition;
    hook = [EmptyHook() for _ in agents],
)
    roles = [agent.role for agent in agents]
    reset!(env)
    observations = [reset(observe(env, agent.role); terminal = true) for agent in agents]
    for (h, agent, obs) in zip(hook, agents, observations)
        h(PRE_EPISODE_STAGE, agent, env, obs)
    end

    while true
        # async here?
        actions = [agent(obs) for (agent, obs) in zip(agents, observations)]
        for (h, agent, obs, action) in zip(hook, agents, observations, actions)
            h(PRE_ACT_STAGE, agent, env, obs => action)
            update!(agent, obs => action)
        end

        env(roles => actions)

        observations = [observe(env, agent.role) for agent in agents]
        for (h, agent, action, obs) in zip(hook, agents, actions, observations)
            h(POST_ACT_STAGE, agents, env, action => obs)
        end

        if stop_condition(agents, env, observations)
            if get_terminal(observations)
                for (h, agent, obs) in zip(hook, agents, observations)
                    h(POST_EPISODE_STAGE, agent, env, obs)
                end
            end
            for (agent, obs) in zip(agents, observations)
                update!(agent, obs => agent(obs))  # push reward terminal into buffer, state and action are not that important
            end
            break
        end

        if get_terminal(observations)
            for (h, agent, obs) in zip(hook, agents, observations)
                h(POST_EPISODE_STAGE, agent, env, obs)
            end

            reset!(env)

            observations = [reset(
                observe(env, agent.role);
                reward = observations[i].reward,
                terminal = observations[i].terminal,
            ) for (i, agent) in enumerate(agents)]

            for (h, agent, obs) in zip(hook, agents, observations)
                h(PRE_EPISODE_STAGE, agent, env, obs)
            end
        end
    end
    hook
end


text:
"""
    chain = mcmc(θ, reps, burnin, Prior, lnL, Proposal)

    Simple MH MCMC

    You must set the needed functions, e.g.,:
    Prior = θ -> your_prior(θ, whatever other args)
    lnL = θ -> your_log_likelihood(θ, whatever other args)
    Proposal = θ -> your_proposal(θ, whatever other args)
    (optionally) mcmcProposalDensity = (θᵗ,θ) -> your_proposal_density

    then get the chain using the above syntax, or optionally,
    (non-symmetric proposal) chain = mcmc(θ, reps, burnin, Prior, lnL, Proposal, ProposalDensity, report=true)
    (example code) mcmc(): runs a simple example. edit(mcmc,()) to see the code

"""
function mcmc()
    println("mcmc(), called with no arguments, runs a simple example")
    println("execute edit(mcmc,()) to see the code")
    # sample is from exponential, prior is lognormal, proposal is random walk lognormal
    y = rand(Exponential(3.0),30)
    # set prior, likelihood and proposal
    Prior = θ -> pdf.(Ref(LogNormal(1.0,1.0)), θ)
    lnL = θ -> sum(logpdf.(Ref(Exponential(θ)), y))
    tuning = 0.5
    Proposal = θ -> rand(LogNormal(log(θ),tuning))
    # get the chain, plot posterior, and descriptive stats
    chain = mcmc(1.0, 100000, 10000, Prior, lnL, Proposal, true) # start value, chain length, and burnin
    p = npdensity(chain[:,1]) # nonparametric plot of posterior density
    plot!(p, title="posterior density, simple MCMC example: true value = 3.0", show=true) # add a title
    dstats(chain)
    return
end

# method using threads and symmetric proposal
@views function mcmc(θ, reps::Int64, burnin::Int64, Prior::Function, lnL::Function, Proposal::Function, report::Bool, nthreads::Int64)
    perthread = reps ÷ nthreads
    chain = zeros(reps, size(θ,1)+1)
    Threads.@threads for t = 1:nthreads # collect the results from the threads
        chain[t*perthread-perthread+1:t*perthread,:] = mcmc(θ, perthread, burnin, Prior, lnL, Proposal, report)
    end
    return chain
end

# method symmetric proposal
# the main loop
@views function mcmc(θ, reps::Int64, burnin::Int64, Prior::Function, lnL::Function, Proposal::Function, report::Bool=true)
    reportevery = Int((reps+burnin)/10)
    lnLθ = lnL(θ)
    chain = zeros(reps, size(θ,1)+1)
    naccept = zeros(size(θ))
    for rep = 1:reps+burnin
        θᵗ = Proposal(θ) # new trial value
        if report
            changed = Int.(.!(θᵗ .== θ)) # find which changed
        end
        # MH accept/reject: only evaluate logL if proposal is in support of prior (avoid crashes)
        pt = Prior(θᵗ)
        accept = false
        if pt > 0.0
            lnLθᵗ = lnL(θᵗ)
            accept = rand() < exp(lnLθᵗ-lnLθ) * pt/Prior(θ)
            if accept
                θ = θᵗ
                lnLθ = lnLθᵗ
            end
        end
        if report
            naccept = naccept .+ changed .* Int.(accept)
        end
        if (mod(rep,reportevery)==0 && report)
            println("current parameters: ", round.(θ,digits=3))
            println("  acceptance rates: ", round.(naccept/reportevery,digits=3))
            naccept = naccept - naccept
        end
        if rep > burnin
            chain[rep-burnin,:] = vcat(θ, accept)
        end
    end
    return chain
end

#=
# the main loop
function mcmc(θ, reps, burnin, Prior, lnL, Proposal::Function, ProposalDensity::Function, report=true::Bool)
    reportevery = Int((reps+burnin)/10)
    lnLθ = lnL(θ)
    chain = zeros(reps, size(θ,1)+1)
    naccept = zeros(size(θ))
    for rep = 1:reps+burnin
        θᵗ = Proposal(θ) # new trial value
        if report
            changed = Int.(.!(θᵗ .== θ)) # find which changed
        end
        # MH accept/reject: only evaluate logL if proposal is in support of prior (avoid crashes)
        pt = Prior(θᵗ)
        accept = false
        if pt > 0.0
            lnLθᵗ = lnL(θᵗ)
            accept = rand() <
            exp(lnLθᵗ-lnLθ)*
            pt/Prior(θ)*
            ProposalDensity(θ,θᵗ)/ProposalDensity(θᵗ,θ)
            if accept
                θ = θᵗ
                lnLθ = lnLθᵗ
            end
        end
        if report
            naccept = naccept .+ changed .* Int.(accept)
        end
        if (mod(rep,reportevery)==0 && report)
            println("current parameters: ", round.(θ,digits=3))
            println("  acceptance rates: ", round.(naccept/reportevery,digits=3))
            naccept = naccept - naccept
        end
        if rep > burnin
            chain[rep-burnin,:] = [θ; accept]
        end
    end
    return chain
end
=#

max_stars_repo_name: mcreel/Econometrics

max_issues_repo_name: mcreel/Econometrics

max_forks_repo_name: mcreel/Econometrics


text:
function SHEMS_optimizer_peer(sh, pv, hp, fh, hw, b, m, p)
    flows = [:PV_DE, :B_DE, :GR_DE, :PV_B, :PV_GR, :PV_HP, :GR_HP, :B_HP, :HP_FH, :HP_HW, :PV_PM, :B_PM, :PM_DE, :PM_HP];
    p_concr = 2400.0;   # kg/m^3
    c_concr = 1.0;      # kJ/(kg*°C)
    p_water = 997.0;    # kg/m^3
    c_water = 4.184;    # kJ/(kg*°C)
    M=1000;

    # Input data__________________________________________________________________________________________________________________________________________________________
    df = CSV.read("data/200124_datafile_all_details_right_timestamp.csv");
    h_last = sh.h_start + m.h_predict -1;                     # optimization horizon

    # all peers have the same demand
    d_e = df[sh.h_start:h_last,:electkwh];                      # electricity demand 1 year from BeOpt
    d_fh = df[sh.h_start:h_last,:heatingkwh];                   # heating demand 1 year from BeOpt
    d_hw = df[sh.h_start:h_last,:hotwaterkwh];                  # hot water demand 1 year from BeOpt
    g_e = repeat([1 1 0 0 0] .*df[sh.h_start:h_last,:PV_generation], inner=(1, sh.n_market));       # electricity generation (from Renewable ninjas), consumer has no generation
    t_outside = df[sh.h_start:h_last,:Temperature];             # °C (from Renewable ninjas)

    # Calculate coefficients of performance for every time period (1:h_predict)___________________________________________
    cop_fh = COPcalc(fh, t_outside);
    cop_hw = COPcalc(hw, t_outside);

    # Model start________________________________________________________________________________________________________________________________________________________
    #Define Model und Solver and settings
    if m.solver == "Cbc"
        model = Model(optimizer_with_attributes(() -> Cbc.Optimizer()));
        set_optimizer_attribute(model, "ratioGap", m.mip_gap);
        set_optimizer_attribute(model, "logLevel", m.output_flag);
    elseif m.solver == "Gurobi"
        model = Model(optimizer_with_attributes(() -> Gurobi.Optimizer(GUROBI_ENV)));
        set_optimizer_attribute(model, "MIPGap", m.mip_gap);
        set_optimizer_attribute(model, "Presolve", m.presolve_flag);
        set_optimizer_attribute(model, "OutputFlag", m.output_flag);
    end

    @variables(model, begin
        X[1:m.h_predict, 1:length(flows), 1:sh.n_peers*sh.n_market] >= 0;    #1:PV_DE, 2:B_DE, 3:GR_DE, 4:PV_B, 5:PV_GR, 6:PV_HP, 7:GR_HP, 8:B_HP, 9:HP_FH, 10:HP_HW, 11:PV_PM, 12:B_PM, 13:PM_DE, 14:PM_HP
        SOC_b[1:m.h_predict+1, 1:sh.n_peers*sh.n_market]  >= 0;
        B_switch[1:m.h_predict, 1:sh.n_peers*sh.n_market], Bin;
        T_fh[1:m.h_predict+1, 1:sh.n_peers*sh.n_market] >= 0; V_hw[1:m.h_predict+1, 1:sh.n_peers*sh.n_market] >= 0;
        Mod_fh[1:m.h_predict, 1:sh.n_peers*sh.n_market] >= 0; Mod_hw[1:m.h_predict, 1:sh.n_peers*sh.n_market]   >= 0;
        T_fh_plus[1:m.h_predict, 1:sh.n_peers*sh.n_market] >= 0; T_fh_minus[1:m.h_predict, 1:sh.n_peers*sh.n_market] >= 0;
        V_hw_plus[1:m.h_predict, 1:sh.n_peers*sh.n_market] >= 0; V_hw_minus[1:m.h_predict, 1:sh.n_peers*sh.n_market] >= 0;
        HP_switch[1:m.h_predict, 1:sh.n_peers*sh.n_market], Bin;
        Hot[1:m.h_predict, 1:sh.n_peers*sh.n_market], Bin;
        Peer_switch[1:m.h_predict, 1:sh.n_peers*sh.n_market], Bin;
    end)

    # Fix start SoCs
    fix.(SOC_b[1,:],  sh.soc_b[:];  force=true);
    fix.(T_fh[1,:], sh.soc_fh[:]; force=true);
    fix.(V_hw[1,:], sh.soc_hw[:]; force=true);

    # Peer 1 with full functionality PV and battery (Prosumager)
    if sh.n_peers >= 2
        # Peer 2 has no Battery (Prosumer)
        fix.(SOC_b[:, (1*sh.n_market)+1:(2*sh.n_market)], 0; force=true);
        fix.(X[:, [2,4,8,12], (1*sh.n_market)+1:(2*sh.n_market)], 0; force=true);  # 2:B_DE, 4:PV_B, 8:B_HP, 12:B_PM

        if sh.n_peers >= 3
            # Peer 3 has no PV and no Battery (Consumer)
            fix.(SOC_b[:, (2*sh.n_market)+1:(3*sh.n_market)], 0; force=true);
            fix.(X[:, [1,2,4,5,6,8,11,12], (2*sh.n_market)+1:(3*sh.n_market)], 0; force=true);  # 1:PV_DE, 2:B_DE, 4:PV_B, 5:PV_GR, 6:PV_HP, 8:B_HP, 11:PV_PM, 12:B_PM

            if sh.n_peers >= 4
                # Peer 4 has no PV but a Battery (Consumer + B)
                fix.(SOC_b[:, (3*sh.n_market)+1:(4*sh.n_market)], 0; force=true);
                fix.(X[:, [1,2,4,5,6,8,11,12], (3*sh.n_market)+1:(4*sh.n_market)], 0; force=true);  # 1:PV_DE, 2:B_DE, 4:PV_B, 5:PV_GR, 6:PV_HP, 8:B_HP, 11:PV_PM, 12:B_PM

                if sh.n_peers >= 5
                    # Peer 5 has no PV and no Battery (Consumer)
                    fix.(SOC_b[:, (4*sh.n_market)+1:(5*sh.n_market)], 0; force=true);
                    fix.(X[:, [1,2,4,5,6,8,11,12], (4*sh.n_market)+1:(5*sh.n_market)], 0; force=true);  # 1:PV_DE, 2:B_DE, 4:PV_B, 5:PV_GR, 6:PV_HP, 8:B_HP, 11:PV_PM, 12:B_PM
                end

            end
        end
    end

    if sh.market_flag < 1 # market is switched off
            fix.(X[:, [11,12,13,14], :], 0; force=true);
    end

    #1:PV_DE, 2:B_DE, 3:GR_DE, 4:PV_B, 5:PV_GR, 6:PV_HP, 7:GR_HP, 8:B_HP, 9:HP_FH, 10:HP_HW, 11:PV_PM, 12:B_PM, 13:PM_DE, 14:PM_HP
    # Objective function: maximize profit, minimize comfort violations_____________________________________________________________________
    @objective(model, Max, sum( sum( sum(
                    ([p.feedin p.grid p.peer p.network p.tax p.eeg p.others] *p.matrix)[i] *X[h,i,n] for i=1:14) -
                    (100 *(V_hw_plus[h,n] +V_hw_minus[h,n] +T_fh_plus[h,n] +T_fh_minus[h,n]))  for h=1:m.h_predict) for n=1:sh.n_peers*sh.n_market));

    # Electricity demand, generation, market clearing___________________________________
    @constraints(model, begin
        [h=1:m.h_predict, n=1:sh.n_peers*sh.n_market],      sum(X[h,i,n] for i=[1,2,3,13]) == d_e[h];                                   # fulfill energy demand (== for max. self-consmption)
        [h=1:m.h_predict, n=1:sh.n_peers*sh.n_market],      sum(X[h,i,n] for i=[1,4,5,6,11]) /pv.eta_i == g_e[h,n];                     # restricted by PV generation
        [h=1:m.h_predict],                                  sum(sum(X[h,i,n] for i=[11,12]) for n=1:sh.n_peers*sh.n_market) ==
                                                                    sum(sum(X[h,i,n] for i=[13,14]) for n=1:sh.n_peers*sh.n_market);    # market clearing peer market
    end)
    # Battery__________________________________________________________________________
    @constraints(model, begin
        [h=1:m.h_predict-1, n=1:sh.n_peers*sh.n_market],    SOC_b[h+1,n] == ((1 -b.loss) *SOC_b[h,n])+
                                                                X[h,4,n] -sum(X[h,i,n] for i=[2,8,12]) /b.eta_b;            # State of Charge, loss for unique solutions
        [h=1:m.h_predict, n=1:sh.n_peers*sh.n_market],      b.soc_min <= SOC_b[h,n] <= b.soc_max;                           # Limits Battery usable capacity
        [h=1:m.h_predict, n=1:sh.n_peers*sh.n_market],      X[h,4,n] <= B_switch[h,n]*b.rate_max;                           # limit charging to nominal power (never at same time)
        [h=1:m.h_predict, n=1:sh.n_peers*sh.n_market],      sum(X[h,i,n] for i=[2,8,12]) <= (1-B_switch[h,n])*b.rate_max;   # limit discharging to nominal power (never at same time)
    end)
    # Heat pump_______________________________________________________________________________
    @constraints(model, begin
        [h=1:m.h_predict, n=1:sh.n_peers*sh.n_market],      sum(X[h,i,n] for i=[9,10]) ==
                                                                sum(X[h,i,n] for i=[6,7,8,14])                   # level power heat pump=and out (== for max. self-consmption)
        [h=1:m.h_predict, n=1:sh.n_peers*sh.n_market],      X[h,9,n] == Mod_fh[h,n] *hp.rate_max;                # heating energy FH max cap
        [h=1:m.h_predict, n=1:sh.n_peers*sh.n_market],      X[h,10,n] == Mod_hw[h,n] *hp.rate_max;               # heating energy HW max cap
        #[h=1:m.h_predict, n=1:sh.n_peers*sh.n_market],     Mod_hw[h,n] + Mod_fh[h,n] >= 0.1;                    # minimum level on?
        [h=1:m.h_predict, n=1:sh.n_peers*sh.n_market],      Mod_fh[h,n] <= 1 -HP_switch[h,n];                    # switch FH or HW
        [h=1:m.h_predict, n=1:sh.n_peers*sh.n_market],      Mod_hw[h,n] <= HP_switch[h,n];                       # switch HW or FH
    end)
    # Floor_heating__________________________________________________________________________________________
    @constraints(model, begin
        [h=1:m.h_predict-1, n=1:sh.n_peers*sh.n_market],    T_fh[h+1,n] == T_fh[h,n]+
                                                                (60*60)/(p_concr *fh.volume *c_concr)*
                                                                ((cop_fh[h] *X[h,9,n]) -d_fh[h]-
                                                                ((1-Hot[h,n]) *fh.loss) +(Hot[h,n] *fh.loss));       # SoC floor heating (temperature)
        [h=1:m.h_predict, n=1:sh.n_peers*sh.n_market],      T_fh[h,n] -((1 -Hot[h,n])*m.big) <= t_outside[h];      # force hot binary on if hotter outside than inside
        [h=1:m.h_predict, n=1:sh.n_peers*sh.n_market],      t_outside[h] -(Hot[h,n]*m.big) <= T_fh[h,n];
        [h=1:m.h_predict, n=1:sh.n_peers*sh.n_market],      T_fh[h,n] <= fh.soc_max +T_fh_plus[h,n];                # Limits temperature FH max
        [h=1:m.h_predict, n=1:sh.n_peers*sh.n_market],      fh.soc_min -T_fh_minus[h,n] <= T_fh[h,n];               # Limits temperature FH min
    end)
    # Hot water_____________________________________________________________________________________________
    @constraints(model, begin
        [h=1:m.h_predict-1, n=1:sh.n_peers*sh.n_market],    V_hw[h+1,n] == V_hw[h,n]+
                                                                (60*60)/((p_water *hw.t_supply *c_water)/1000)*
                                                                ((cop_hw[h] *X[h,10,n]) -d_hw[h] -hw.loss);         # SoC hot water (volume)
        [h=1:m.h_predict, n=1:sh.n_peers*sh.n_market],      V_hw[h,n] <= hw.soc_max +V_hw_plus[h,n];                # Limits volume HW max
        [h=1:m.h_predict, n=1:sh.n_peers*sh.n_market],      hw.soc_min - V_hw_minus[h,n] <= V_hw[h,n];              # Limits volume HW min

    end)
    # Peer market__________________________________________________________________________
    @constraints(model, begin
        [h=1:m.h_predict, n=1:sh.n_peers*sh.n_market],      sum(X[h,i,n] for i=[11,12]) <= Peer_switch[h,n]*M;             # Feed-in peer market
        [h=1:m.h_predict, n=1:sh.n_peers*sh.n_market],      sum(X[h,i,n] for i=[13,14]) <= (1-Peer_switch[h,n])*M;      # Purchase peer market
    end)

    JuMP.optimize!(model);

    # collect returns
    results = Array{Float64}(undef,(m.h_control*sh.n_peers*sh.n_market, 11+14+3+1));    # optimization variables + flows + month + day + hour
    for n=1:sh.n_peers*sh.n_market

        profits = sum( ([p.feedin p.grid p.peer p.network p.tax p.eeg p.others] *p.matrix)[i] .*JuMP.value.(X[1:m.h_control,i,n]) for i=1:14);

        results[((n-1)*m.h_control)+1:(n*m.h_control),:] =
            hcat(JuMP.value.(T_fh[1:m.h_control,n]), JuMP.value.(V_hw[1:m.h_control,n]), JuMP.value.(SOC_b[1:m.h_control,n]),
            JuMP.value.(T_fh_plus[1:m.h_control,n]), JuMP.value.(T_fh_minus[1:m.h_control,n]),
            JuMP.value.(V_hw_plus[1:m.h_control,n]), JuMP.value.(V_hw_minus[1:m.h_control,n]),
            profits,
            JuMP.value.(Peer_switch[1:m.h_control,n]),
            cop_fh[1:m.h_control], cop_hw[1:m.h_control], (ones(m.h_control)*(mod((div(n-1, sh.n_market)+1)-1,sh.n_peers)+1)),
            JuMP.value.(X[1:m.h_control,:,n]),                              # flow variables
            Matrix(df[sh.h_start:(sh.h_start + m.h_control -1), 13:15]));   # month + day + hour
    end

    return round.(JuMP.value.(SOC_b[m.h_control+1,:]), digits=5), JuMP.value.(T_fh[m.h_control+1,:]), JuMP.value.(V_hw[m.h_control+1,:]), results
end

max_stars_repo_name: lilanger/PEERS

max_issues_repo_name: lilanger/PEERS

max_forks_repo_name: lilanger/PEERS


text:
using Test, OrdinaryDiffEq

prob = ODEProblem(function (du,u,p,t)
                    du[1] = 10.0(u[2]-u[1])
                    du[2] = u[1]*(28.0-u[3]) - u[2]
                    du[3] = u[1]*u[2] - (8/3)*u[3]
                  end,
                  [1.0;0.0;0.0], (0.0, 100.0))

function buildW(...)::NTuple{3}
    nothing #= left preconditioner =#, nothing #= right preconditioner =#
end

nlsolve = NLNewton(buildW, )
sol = solve(prob, TRBDF2(nlsolve=nlsolve))


text:
__precompile__()

"""
Module for studying stream graphs and link streams.
"""

module StreamGraphs

#######################################
# Imports

import Base.≈
import Base.==
import Base.∈
import Base.∪
import Base.∩
import Base.⊆
import Base.+
import Base.-
import Base: length, count, match, merge, merge!, string

abstract type AbstractStream end
abstract type AbstractDirectedStream <: AbstractStream end
abstract type AbstractUndirectedStream <: AbstractStream end
abstract type StreamObject end
abstract type AbstractPath end
abstract type Event end

########################################
# Exports

export

		# Stream components
		Intervals, Node, Link, LinkStream, StreamGraph,
		DirectedLinkStream, DirectedStreamGraph,
		times, nodes, links, links_from, links_to,
		neighborhood, degree, average_node_degree,
		average_time_degree, State, Transition, TimeCursor,
		degrees, average_expected_degree, expected_degree,

		# Paths
		Jump, DurationJump, Path, DurationPath,

		# Adding to streams
		add_node!, add_link!, record!, load!, parse_line,
		record_event!,

		# Durations
		duration, node_duration, link_duration, contribution,
		number_of_nodes, number_of_links, density, node_contribution,
		link_contribution,

		node_clustering, time_clustering, clustering, coverage,
		uniformity,

		Δnodes, Δlinks, apply!, next_transition, previous_transition,
		next!, previous!, goto!, start!, end!, Event, NodeEvent, LinkEvent,

		parse_line_auv, parse_line_abuv, parse_line, parse_to_events,
		string, dump,

		# functions
		count, length, merge, merge!, push, clean, from_match, to_match,
		match, get_idx, is_connected, is_valid, start, finish,

		# Operators
		==, ≈, ∈, ∪, ∩, ⊆, +, -, ×, ⊗

#######################################
# Includes

include("sets.jl")
include("intervals.jl")
include("base.jl")
include("nodes.jl")
include("links.jl")
include("times.jl")
include("size.jl")
include("durations.jl")
include("contribution.jl")
include("coverage.jl")
include("compactness.jl")
include("uniformity.jl")
include("density.jl")
include("clustering.jl")
include("degree.jl")
include("additions.jl")
include("parsers.jl")

end


text:
struct SnpLinAlg{T} <: AbstractMatrix{T}
    s::SnpArray
    model::Union{Val{1}, Val{2}, Val{3}}
    center::Bool
    scale::Bool
    impute::Bool
    μ::Vector{T}
    σinv::Vector{T}
    storagev1::Vector{T} # size(s, 1)
    storagev2::Vector{T} # size(s, 2)
    storagev3::Vector{T} # size(s, 2)
end

AbstractSnpLinAlg = Union{SnpLinAlg, SubArray{T, 1, SnpLinAlg{T}},
    SubArray{T, 2, SnpLinAlg{T}}} where T

"""
    SnpLinAlg{T}(s; model=ADDITIVE_MODEL, center=false, scale=false, impute=true)

Pad a `SnpArray` with some parameters for linear algebera operation.

# Arguments
- s: a `SnpArray`.
- model: one of `ADDITIVE_MODEL`(default), `DOMINANT_MODEL`, `RECESSIVE_MODEL`.
- center: whether to center (default: false).
- scale: whether to scale to standard deviation 1 (default: false).
- impute: whether to impute missing value with column mean (default: true).
"""
function SnpLinAlg{T}(
    s::AbstractSnpArray;
    model = ADDITIVE_MODEL,
    center::Bool = false,
    scale::Bool = false,
    impute::Bool = true) where T <: AbstractFloat
    μ = dropdims(mean(s; dims=1, model=model), dims=1)
    σinv = Vector{T}(undef, size(s, 2))
    if model == ADDITIVE_MODEL
        @inbounds @simd for j in 1:size(s, 2)
            σinv[j] = sqrt(μ[j] * (1 - μ[j] / 2))
            σinv[j] = σinv[j] > 0 ? inv(σinv[j]) : one(T)
        end
    elseif model == DOMINANT_MODEL || model == RECESSIVE_MODEL
        @inbounds @simd for j in 1:size(s, 2)
            σinv[j] = sqrt(μ[j] * (1 - μ[j]))
            σinv[j] = σinv[j] > 0 ? inv(σinv[j]) : one(T)
        end
    else
        throw(ArgumentError("unrecognized model $model"))
    end
    storagev1 = Vector{T}(undef, size(s, 1))
    storagev2 = Vector{T}(undef, size(s, 2))
    storagev3 = Vector{T}(undef, size(s, 2))
    SnpLinAlg{T}(s, model, center, scale, impute, μ, σinv, storagev1, storagev2, storagev3)
end

Base.size(sla::SnpLinAlg) = size(sla.s)
Base.size(sla::SnpLinAlg, k::Integer) = size(sla.s, k)

eltype(bm::SnpLinAlg) = eltype(bm.μ)

function Base.getindex(s::SnpLinAlg{T}, i::Int, j::Int) where T
    x = SnpArrays.convert(T, getindex(s.s, i, j), s.model)
    s.impute && isnan(x) && return s.μ[j]
    s.center && (x -= s.μ[j])
    s.scale && (x *= s.σinv[j])
    return x
end

# macros taken from Gaius.jl
macro _spawn(ex)
    if Threads.nthreads() > 1
        esc(Expr(:macrocall, Expr(:(.), :Threads, QuoteNode(Symbol("@spawn"))), __source__, ex))
    else
        esc(ex)
    end
end
macro _sync(ex)
    if Threads.nthreads() > 1
        esc(Expr(:macrocall, Symbol("@sync"), __source__, ex))
    else
        esc(ex)
    end
end

"""
    LinearAlgebra.mul!(out, sla::SnpLinAlg, v)

In-place matrix-vector multiplication.
"""
function mul!(
    out::AbstractVector{T},
    sla::SnpLinAlg{T},
    v::AbstractVector{T}) where T <: AbstractFloat
    @assert length(out) == size(sla, 1) && length(v) == size(sla, 2)
    if sla.scale
        sla.storagev2 .= sla.σinv .* v
        w = sla.storagev2
    else
        w = v
    end
    fill!(out, zero(eltype(out)))
    s = sla.s

    _snparray_ax_tile!(out, s.data, w, sla.model, sla.μ, sla.impute, s.m)

    if sla.center
        return out .-= dot(sla.μ, w)
    else
        return out
    end
end

"""
    LinearAlgebra.mul!(out, sla::SnpLinAlg, V)

In-place matrix-matrix multiplication.
"""
function mul!(
    out::AbstractMatrix{T},
    sla::SnpLinAlg{T},
    V::AbstractMatrix{T}) where T <: AbstractFloat
    @assert size(out, 1) == size(sla, 1) && size(out, 2) == size(V, 2) && size(sla, 2) == size(V, 1)
    sla.storagev2 .= sla.center ? sla.μ : zero(T) # for on-the-fly centering
    sla.storagev3 .= sla.scale ? sla.σinv : one(T) # for on-the-fly scaling

    fill!(out, zero(eltype(out)))
    s = sla.s

    _snparray_AX_tile!(out, s.data, V, sla.model, sla.storagev2, sla.μ, sla.impute, s.m, sla.storagev3)

    return out
end

"""
    LinearAlgebra.mul!(out, st::Union{Transpose{T, SnpLinAlg{T}}, Adjoint{T, SnpLinAlg{T}}}, v)

In-place matrix-vector multiplication, with transposed `SnpLinAlg`.
"""
function mul!(
    out::AbstractVector{T},
    st::Union{Transpose{T, SnpLinAlg{T}}, Adjoint{T, SnpLinAlg{T}}},
    v::AbstractVector{T}) where T <: AbstractFloat
    @assert length(out) == size(st, 1) && length(v) == size(st, 2)
    sla = st.parent
    s = sla.s
    fill!(out, zero(eltype(out)))

    _snparray_atx_tile!(out, s.data, v, sla.model, sla.μ, sla.impute, s.m)
    if sla.center
        out .-= sum(v) .* sla.μ
    end
    if sla.scale
        return out .*= sla.σinv
    else
        return out
    end
end

"""
    LinearAlgebra.mul!(out, st::Union{Transpose{T, SnpLinAlg{T}}, Adjoint{T, SnpLinAlg{T}}}, V)

In-place matrix-matrix multiplication, with transposed `SnpLinAlg`.
"""
function mul!(
    out::AbstractMatrix{T},
    st::Union{Transpose{T, SnpLinAlg{T}}, Adjoint{T, SnpLinAlg{T}}},
    V::AbstractMatrix{T}) where T <: AbstractFloat
    sla = st.parent
    @assert size(out, 1) == size(sla, 2) && size(out, 2) == size(V, 2) && size(sla, 1) == size(V, 1)

    sla.storagev2 .= sla.center ? sla.μ : zero(T) # for on-the-fly centering
    sla.storagev3 .= sla.scale ? sla.σinv : one(T) # for on-the-fly scaling

    s = sla.s
    fill!(out, zero(eltype(out)))

    _snparray_AtX_tile!(out, s.data, V, sla.model, sla.storagev2, sla.μ, sla.impute, s.m, sla.storagev3)

    return out
end

wait(::Nothing) = nothing

function _snparray_ax_tile!(c, A, b, model, μ, impute, rows_filled)
    vstep = 1024
    hstep = 1024
    vstep_log2 = 10
    hstep_log2 = 10

    if !impute
        if model == ADDITIVE_MODEL
            _ftn! = _snparray_ax_additive!
        elseif model == DOMINANT_MODEL
            _ftn! = _snparray_ax_dominant!
        else
            _ftn! = _snparray_ax_recessive!
        end
    else
        if model == ADDITIVE_MODEL
            _ftn! = _snparray_ax_additive_meanimpute!
        elseif model == DOMINANT_MODEL
            _ftn! = _snparray_ax_dominant_meanimpute!
        else
            _ftn! = _snparray_ax_recessive_meanimpute!
        end
    end

    M = length(c) >> 2
    N = size(A, 2)
    Miter = M >>> vstep_log2 # fast div(M, 1024)
    Mrem = rows_filled & (vstep << 2 - 1) # fast rem(rows_filled, 4vstep)
    Niter = N >>> hstep_log2
    Nrem = N & (hstep - 1)
    taskarray = Array{Any}(undef, Miter + 1)
    fill!(taskarray, nothing)
    @_sync begin
        GC.@preserve c A b for n in 0:Niter - 1
            for m in 0:Miter - 1
                wait(taskarray[m+1])
                taskarray[m+1] = @_spawn _ftn!(
                    gesp(stridedpointer(c), (4 * vstep * m,)),
                    gesp(stridedpointer(A), (vstep * m, hstep * n)),
                    gesp(stridedpointer(b), (hstep * n,)),
                    vstep << 2, hstep, @view(μ[hstep * n + 1:hstep * (n + 1)])
                )
            end
            if Mrem != 0
                wait(taskarray[Miter+1])
                taskarray[Miter+1] = @_spawn _ftn!(
                    @view(c[4 * vstep * Miter + 1:end]),
                    @view(A[vstep * Miter + 1:end, hstep * n + 1:hstep * (n + 1)]),
                    @view(b[hstep * n + 1:hstep * (n + 1)]),
                    length(c) - 4 * vstep * Miter, hstep,
                    @view(μ[hstep * n + 1:hstep * (n + 1)])
                )
            end
        end
        if Nrem != 0
            for m in 0:Miter-1
                wait(taskarray[m+1])
                taskarray[m+1] = @_spawn _ftn!(
                    @view(c[4 * vstep * m + 1:4 * vstep * (m + 1)]),
                    @view(A[vstep * m + 1:vstep * (m + 1), hstep * Niter + 1:end]),
                    @view(b[hstep * Niter + 1:end]),
                    vstep << 2,
                    Nrem, @view(μ[hstep * Niter + 1:end])
                )
            end
            if Mrem != 0
                wait(taskarray[Miter + 1])
                taskarray[Miter + 1] = @_spawn _ftn!(
                    @view(c[4 * vstep * Miter+1:end]),
                    @view(A[vstep * Miter + 1:end, hstep * Niter + 1:end]),
                    @view(b[hstep * Niter + 1:end]),
                    length(c) - 4 * vstep * Miter,
                    Nrem, @view(μ[hstep * Niter + 1:end])
                )
            end
        end
    end
end

# μ[i] is mean of SNP i, and μimpute[i] is used to impute missings for SNP i
function _snparray_AX_tile!(C, A, B, model, μ, μimpute, impute, rows_filled, σinv)
    vstep = 256
    hstep = 256
    pstep = 256
    vstep_log2 = 8
    hstep_log2 = 8
    pstep_log2 = 8

    if !impute
        if model == ADDITIVE_MODEL
            _ftn! = _snparray_AX_additive!
        elseif model == DOMINANT_MODEL
            _ftn! = _snparray_AX_dominant!
        else
            _ftn! = _snparray_AX_recessive!
        end
    else
        if model == ADDITIVE_MODEL
            _ftn! = _snparray_AX_additive_meanimpute!
        elseif model == DOMINANT_MODEL
            _ftn! = _snparray_AX_dominant_meanimpute!
        else
            _ftn! = _snparray_AX_recessive_meanimpute!
        end
    end

    M = size(C, 1) >> 2
    N = size(A, 2)
    P = size(C, 2)
    Miter = M >>> vstep_log2 # fast div(M, 1024)
    Mrem = rows_filled & (vstep << 2 - 1) # fast rem(rows_filled, 4vstep)
    Niter = N >>> hstep_log2
    Nrem = N & (hstep - 1)
    Piter = P >>> pstep_log2
    Prem = P & (pstep - 1)
    taskarray = Array{Any}(undef, Miter + 1)
    fill!(taskarray, nothing)
    @_sync begin
        GC.@preserve C A B for p in 0:Piter - 1
            for n in 0:Niter - 1
                for m in 0:Miter - 1
                    wait(taskarray[m+1])
                    taskarray[m+1] = @_spawn _ftn!(
                        @view(C[4 * vstep * m + 1:4 * vstep * (m + 1), pstep * p + 1:pstep * (p + 1)]),
                        @view(A[vstep * m + 1:vstep * (m + 1), hstep * n + 1:hstep * (n + 1)]),
                        @view(B[hstep * n + 1:hstep * (n + 1), pstep * p + 1:pstep * (p + 1)]),
                        vstep << 2, hstep, pstep,
                        @view(μ[hstep * n + 1:hstep * (n + 1)]),
                        @view(μimpute[hstep * n + 1:hstep * (n + 1)]),
                        @view(σinv[hstep * n + 1:hstep * (n + 1)])
                    )
                end
                if Mrem != 0
                    wait(taskarray[Miter+1])
                    taskarray[Miter+1] = @_spawn _ftn!(
                        @view(C[4 * vstep * Miter + 1:end, pstep * p + 1:pstep * (p + 1)]),
                        @view(A[vstep * Miter + 1:end, hstep * n + 1:hstep * (n + 1)]),
                        @view(B[hstep * n + 1:hstep * (n + 1), pstep * p + 1:pstep * (p + 1)]),
                        size(C, 1) - 4 * vstep * Miter, hstep, pstep,
                        @view(μ[hstep * n + 1:hstep * (n + 1)]),
                        @view(μimpute[hstep * n + 1:hstep * (n + 1)]),
                        @view(σinv[hstep * n + 1:hstep * (n + 1)])
                    )
                end
            end
            if Nrem != 0
                for m in 0:Miter-1
                    wait(taskarray[m+1])
                    taskarray[m+1] = @_spawn _ftn!(
                        @view(C[4 * vstep * m + 1:4 * vstep * (m + 1), pstep * p + 1:pstep * (p + 1)]),
                        @view(A[vstep * m + 1:vstep * (m + 1), hstep * Niter + 1:end]),
                        @view(B[hstep * Niter + 1:end, pstep * p + 1:pstep * (p + 1)]),
                        vstep << 2, Nrem, pstep,
                        @view(μ[hstep * Niter + 1:end]),
                        @view(μimpute[hstep * Niter + 1:end]),
                        @view(σinv[hstep * Niter + 1:end])
                    )
                end
                if Mrem != 0
                    wait(taskarray[Miter + 1])
                    taskarray[Miter + 1] = @_spawn _ftn!(
                        @view(C[4 * vstep * Miter+1:end, pstep * p + 1:pstep * (p + 1)]),
                        @view(A[vstep * Miter + 1:end, hstep * Niter + 1:end]),
                        @view(B[hstep * Niter + 1:end, pstep * p + 1:pstep * (p + 1)]),
                        size(C, 1) - 4 * vstep * Miter, Nrem, pstep,
                        @view(μ[hstep * Niter + 1:end]),
                        @view(μimpute[hstep * Niter + 1:end]),
                        @view(σinv[hstep * Niter + 1:end])
                    )
                end
            end
        end
        if Prem != 0
            for n in 0:Niter - 1
                for m in 0:Miter - 1
                    wait(taskarray[m+1])
                    taskarray[m+1] = @_spawn _ftn!(
                        @view(C[4 * vstep * m + 1:4 * vstep * (m + 1), pstep * Piter + 1:end]),
                        @view(A[vstep * m + 1:vstep * (m + 1), hstep * n + 1:hstep * (n + 1)]),
                        @view(B[hstep * n + 1:hstep * (n + 1), pstep * Piter + 1:end]),
                        vstep << 2, hstep, Prem,
                        @view(μ[hstep * n + 1:hstep * (n + 1)]),
                        @view(μimpute[hstep * n + 1:hstep * (n + 1)]),
                        @view(σinv[hstep * n + 1:hstep * (n + 1)])
                    )
                end
                if Mrem != 0
                    wait(taskarray[Miter+1])
                    taskarray[Miter+1] = @_spawn _ftn!(
                        @view(C[4 * vstep * Miter + 1:end, pstep * Piter + 1:end]),
                        @view(A[vstep * Miter + 1:end, hstep * n + 1:hstep * (n + 1)]),
                        @view(B[hstep * n + 1:hstep * (n + 1), pstep * Piter + 1:end]),
                        size(C, 1) - 4 * vstep * Miter, hstep, Prem,
                        @view(μ[hstep * n + 1:hstep * (n + 1)]),
                        @view(μimpute[hstep * n + 1:hstep * (n + 1)]),
                        @view(σinv[hstep * n + 1:hstep * (n + 1)])
                    )
                end
            end
            if Nrem != 0
                for m in 0:Miter-1
                    wait(taskarray[m+1])
                    taskarray[m+1] = @_spawn _ftn!(
                        @view(C[4 * vstep * m + 1:4 * vstep * (m + 1), pstep * Piter + 1:end]),
                        @view(A[vstep * m + 1:vstep * (m + 1), hstep * Niter + 1:end]),
                        @view(B[hstep * Niter + 1:end, pstep * Piter + 1:end]),
                        vstep << 2, Nrem, Prem,
                        @view(μ[hstep * Niter + 1:end]),
                        @view(μimpute[hstep * Niter + 1:end]),
                        @view(σinv[hstep * Niter + 1:end])
                    )
                end
                if Mrem != 0
                    wait(taskarray[Miter + 1])
                    taskarray[Miter + 1] = @_spawn _ftn!(
                        @view(C[4 * vstep * Miter+1:end, pstep * Piter + 1:end]),
                        @view(A[vstep * Miter + 1:end, hstep * Niter + 1:end]),
                        @view(B[hstep * Niter + 1:end, pstep * Piter + 1:end]),
                        size(C, 1) - 4 * vstep * Miter, Nrem, Prem,
                        @view(μ[hstep * Niter + 1:end]),
                        @view(μimpute[hstep * Niter + 1:end]),
                        @view(σinv[hstep * Niter + 1:end])
                    )
                end
            end
        end
    end
end

function _snparray_atx_tile!(c, A, b, model, μ, impute, rows_filled)
    vstep = 2048
    hstep = 2048
    vstep_log2 = 11
    hstep_log2 = 11

    if !impute
        if model == ADDITIVE_MODEL
            _ftn! = _snparray_atx_additive!
        elseif model == DOMINANT_MODEL
            _ftn! = _snparray_atx_dominant!
        else
            _ftn! = _snparray_atx_recessive!
        end
    else
        if model == ADDITIVE_MODEL
            _ftn! = _snparray_atx_additive_meanimpute!
        elseif model == DOMINANT_MODEL
            _ftn! = _snparray_atx_dominant_meanimpute!
        else
            _ftn! = _snparray_atx_recessive_meanimpute!
        end
    end

    M = length(b) >> 2
    N = size(A, 2)
    Miter = M >>> vstep_log2 # fast div(M, 1024)
    Mrem = rows_filled & (vstep << 2 - 1) # fast rem(rows_filled, 4vstep)
    Niter = N >>> hstep_log2
    Nrem = N & (hstep - 1)
    taskarray = Array{Any}(undef, Niter+1)
    fill!(taskarray, nothing)
    @_sync begin
        GC.@preserve c A b for m in 0:Miter - 1
            for n in 0:Niter - 1
                wait(taskarray[n + 1])
                taskarray[n + 1] = @_spawn _ftn!(
                    gesp(stridedpointer(c), (hstep * n,)),
                    gesp(stridedpointer(A), (vstep * m, hstep * n)),
                    gesp(stridedpointer(b), (4 * vstep * m,)),
                    vstep << 2, hstep, @view(μ[hstep * n + 1:hstep * (n + 1)])
                )
            end
            if Nrem != 0
                wait(taskarray[Niter + 1])
                taskarray[Niter + 1] = @_spawn _ftn!(
                    @view(c[hstep * Niter + 1:end]),
                    @view(A[vstep * m + 1:vstep * (m + 1), hstep * Niter + 1:end]),
                    @view(b[4 * vstep * m + 1:4 * vstep * (m + 1)]),
                    vstep << 2, Nrem, @view(μ[hstep * Niter + 1:end])
                )
            end

        end
        if Mrem != 0
            for n in 0:Niter - 1
                wait(taskarray[n + 1])
                taskarray[n + 1] = @_spawn _ftn!(
                    @view(c[hstep * n + 1:hstep * (n + 1)]),
                    @view(A[vstep * Miter + 1:end, hstep * n + 1:hstep * (n + 1)]),
                    @view(b[4 * vstep * Miter + 1:end]),
                    length(b) - 4 * vstep * Miter,
                    hstep, @view(μ[hstep * n + 1:hstep * (n + 1)])
                )
            end
            if Nrem != 0
                wait(taskarray[Niter + 1])
                taskarray[Niter + 1] = @_spawn _ftn!(
                    @view(c[hstep * Niter + 1:end]),
                    @view(A[vstep * Miter + 1:end, hstep * Niter + 1:end]),
                    @view(b[4 * vstep * Miter + 1:end]),
                    length(b) - 4 * vstep * Miter,
                    Nrem, @view(μ[hstep * Niter + 1:end])
                )
            end
        end
    end
end

# μ[i] is mean of SNP i, and μimpute[i] is used to impute missings for SNP i
function _snparray_AtX_tile!(C, A, B, model, μ, μimpute, impute, rows_filled, σinv)
    vstep = 2048
    hstep = 2048
    pstep = 2048
    vstep_log2 = 11
    hstep_log2 = 11
    pstep_log2 = 11

    if !impute
        if model == ADDITIVE_MODEL
            _ftn! = _snparray_AtX_additive!
        elseif model == DOMINANT_MODEL
            _ftn! = _snparray_AtX_dominant!
        else
            _ftn! = _snparray_AtX_recessive!
        end
    else
        if model == ADDITIVE_MODEL
            _ftn! = _snparray_AtX_additive_meanimpute!
        elseif model == DOMINANT_MODEL
            _ftn! = _snparray_AtX_dominant_meanimpute!
        else
            _ftn! = _snparray_AtX_recessive_meanimpute!
        end
    end

    M = size(B, 1) >> 2
    N = size(A, 2)
    P = size(C, 2)
    Miter = M >>> vstep_log2 # fast div(M, 1024)
    Mrem = rows_filled & (vstep << 2 - 1) # fast rem(rows_filled, 4vstep)
    Niter = N >>> hstep_log2
    Nrem = N & (hstep - 1)
    Piter = P >>> pstep_log2
    Prem = P & (pstep - 1)
    taskarray = Array{Any}(undef, Niter + 1)
    fill!(taskarray, nothing)
    @_sync begin
        GC.@preserve C A B for p in 0:Piter - 1
            for m in 0:Miter - 1
                for n in 0:Niter - 1
                    wait(taskarray[n + 1])
                    taskarray[n + 1] = @_spawn _ftn!(
                        @view(C[hstep * n + 1:hstep * (n + 1), pstep * p + 1:pstep * (p + 1)]),
                        @view(A[vstep * m + 1:vstep * (m + 1), hstep * n + 1:hstep * (n + 1)]),
                        @view(B[4 * vstep * m + 1:4 * vstep * (m + 1), pstep * p + 1:pstep * (p + 1)]),
                        vstep << 2, hstep, pstep,
                        @view(μ[hstep * n + 1:hstep * (n + 1)]),
                        @view(μimpute[hstep * n + 1:hstep * (n + 1)]),
                        @view(σinv[hstep * n + 1:hstep * (n + 1)])
                    )
                end
                if Nrem != 0
                    wait(taskarray[Niter + 1])
                    taskarray[Niter + 1] = @_spawn _ftn!(
                        @view(C[hstep * Niter + 1:end, pstep * p + 1:pstep * (p + 1)]),
                        @view(A[vstep * m + 1:vstep * (m + 1), hstep * Niter + 1:end]),
                        @view(B[4 * vstep * m + 1:4 * vstep * (m + 1), pstep * p + 1:pstep * (p + 1)]),
                        vstep << 2, Nrem, pstep,
                        @view(μ[hstep * Niter + 1:end]),
                        @view(μimpute[hstep * Niter + 1:end]),
                        @view(σinv[hstep * Niter + 1:end])
                    )
                end
            end
            if Mrem != 0
                for n in 0:Niter - 1
                    wait(taskarray[n + 1])
                    taskarray[n + 1] = @_spawn _ftn!(
                        @view(C[hstep * n + 1:hstep * (n + 1), pstep * p + 1:pstep * (p + 1)]),
                        @view(A[vstep * Miter + 1:end, hstep * n + 1:hstep * (n + 1)]),
                        @view(B[4 * vstep * Miter + 1:end, pstep * p + 1:pstep * (p + 1)]),
                        size(B, 1) - 4 * vstep * Miter, hstep, pstep,
                        @view(μ[hstep * n + 1:hstep * (n + 1)]),
                        @view(μimpute[hstep * n + 1:hstep * (n + 1)]),
                        @view(σinv[hstep * n + 1:hstep * (n + 1)])
                    )
                end
                if Nrem != 0
                    wait(taskarray[Niter + 1])
                    taskarray[Niter + 1] = @_spawn _ftn!(
                        @view(C[hstep * Niter + 1:end, pstep * p + 1:pstep * (p + 1)]),
                        @view(A[vstep * Miter + 1:end, hstep * Niter + 1:end]),
                        @view(B[4 * vstep * Miter + 1:end, pstep * p + 1:pstep * (p + 1)]),
                        size(B, 1) - 4 * vstep * Miter, Nrem, pstep,
                        @view(μ[hstep * Niter + 1:end]),
                        @view(μimpute[hstep * Niter + 1:end]),
                        @view(σinv[hstep * Niter + 1:end])
                    )
                end
            end
        end
        if Prem != 0
            for m in 0:Miter - 1
                for n in 0:Niter - 1
                    wait(taskarray[n + 1])
                    taskarray[n + 1] = @_spawn _ftn!(
                        @view(C[hstep * n + 1:hstep * (n + 1), pstep * Piter + 1:end]),
                        @view(A[vstep * m + 1:vstep * (m + 1), hstep * n + 1:hstep * (n + 1)]),
                        @view(B[4 * vstep * m + 1:4 * vstep * (m + 1), pstep * Piter + 1:end]),
                        vstep << 2, hstep, Prem,
                        @view(μ[hstep * n + 1:hstep * (n + 1)]),
                        @view(μimpute[hstep * n + 1:hstep * (n + 1)]),
                        @view(σinv[hstep * n + 1:hstep * (n + 1)])
                    )
                end
                if Nrem != 0
                    wait(taskarray[Niter + 1])
                    taskarray[Niter + 1] = @_spawn _ftn!(
                        @view(C[hstep * Niter + 1:end, pstep * Piter + 1:end]),
                        @view(A[vstep * m + 1:vstep * (m + 1), hstep * Niter + 1:end]),
                        @view(B[4 * vstep * m + 1:4 * vstep * (m + 1), pstep * Piter + 1:end]),
                        vstep << 2, Nrem, Prem,
                        @view(μ[hstep * Niter + 1:end]),
                        @view(μimpute[hstep * Niter + 1:end]),
                        @view(σinv[hstep * Niter + 1:end])
                    )
                end
            end
            if Mrem != 0
                for n in 0:Niter - 1
                    wait(taskarray[n + 1])
                    taskarray[n + 1] = @_spawn _ftn!(
                        @view(C[hstep * n + 1:hstep * (n + 1), pstep * Piter + 1:end]),
                        @view(A[vstep * Miter + 1:end, hstep * n + 1:hstep * (n + 1)]),
                        @view(B[4 * vstep * Miter + 1:end, pstep * Piter + 1:end]),
                        size(B, 1) - 4 * vstep * Miter, hstep, Prem,
                        @view(μ[hstep * n + 1:hstep * (n + 1)]),
                        @view(μimpute[hstep * n + 1:hstep * (n + 1)]),
                        @view(σinv[hstep * n + 1:hstep * (n + 1)])
                    )
                end
                if Nrem != 0
                    wait(taskarray[Niter + 1])
                    taskarray[Niter + 1] = @_spawn _ftn!(
                        @view(C[hstep * Niter + 1:end, pstep * Piter + 1:end]),
                        @view(A[vstep * Miter + 1:end, hstep * Niter + 1:end]),
                        @view(B[4 * vstep * Miter + 1:end, pstep * Piter + 1:end]),
                        size(B, 1) - 4 * vstep * Miter, Nrem, Prem,
                        @view(μ[hstep * Niter + 1:end]),
                        @view(μimpute[hstep * Niter + 1:end]),
                        @view(σinv[hstep * Niter + 1:end])
                    )
                end
            end
        end
    end
end

for (_ftn!, _ftn_rem!, expr) in [
        (:_snparray_ax_additive!, :_snparray_ax_additive_rem!,
            :(((Aij >= 2) + (Aij == 3)) * v[j])),
        (:_snparray_ax_dominant!, :_snparray_ax_dominant_rem!,
            :((Aij >= 2)  * v[j])),
        (:_snparray_ax_recessive!, :_snparray_ax_recessive_rem!,
            :((Aij == 3) * v[j])),
        (:_snparray_ax_additive_meanimpute!, :_snparray_ax_additive_meanimpute_rem!,
            :(((Aij >= 2) * 1.0 + (Aij == 3) * 1.0 + (Aij == 1) * μ[j]) * v[j])),
        (:_snparray_ax_dominant_meanimpute!, :_snparray_ax_dominant_meanimpute_rem!,
            :((Aij >= 2) * v[j] + (Aij == 1) * μ[j] * v[j])),
        (:_snparray_ax_recessive_meanimpute!, :_snparray_ax_recessive_meanimpute_rem!,
            :((Aij == 3) * v[j] + (Aij == 1) * μ[j] * v[j]))
    ]
    @eval begin
        function ($_ftn_rem!)(out, s, v, μ)
            maxp = length(out)
            @avx for j in eachindex(v)
                block = s[1, j]
                for p in 1:maxp
                    Aij = (block >> (2 * (p - 1))) & 3
                    out[p] += $expr
                end
            end
        end

        function ($_ftn!)(out, s, v, rows, cols, μ)
            k = rows >> 2
            rem = rows & 3

            if k ≥ 1 # avoid avxing over empty collection
                @avx for j ∈ 1:cols
                    for l in 1:k
                        block = s[l, j]

                        for p in 1:4
                            Aij = (block >> (2 * (p - 1))) & 3
                            out[4 * (l - 1) + p] += $expr
                        end

                    end
                end
            end
            if rem != 0
                ($_ftn_rem!)(@view(out[4k + 1:end]), @view(s[k + 1:k + 1, :]), v, μ)
            end
            nothing
        end
    end
end

for (_ftn!, _ftn_rem!, expr) in [
        (:_snparray_atx_additive!, :_snparray_atx_additive_rem!,
            :(((Aij >= 2) + (Aij == 3)) * v[4 * (l - 1) + p])),
        (:_snparray_atx_dominant!, :_snparray_atx_dominant_rem!,
            :((Aij >= 2)  * v[4 * (l - 1) + p])),
        (:_snparray_atx_recessive!, :_snparray_atx_recessive_rem!,
            :((Aij == 3) * v[4 * (l - 1) + p])),
        (:_snparray_atx_additive_meanimpute!, :_snparray_atx_additive_meanimpute_rem!,
            :(((Aij >= 2) * 1.0 + (Aij == 3) * 1.0 + (Aij == 1) * μ[i]) *  v[4 * (l - 1) + p])),
        (:_snparray_atx_dominant_meanimpute!, :_snparray_atx_dominant_meanimpute_rem!,
            :((Aij >= 2) * v[4 * (l - 1) + p] + (Aij == 1) * μ[i] * v[4 * (l - 1) + p])),
        (:_snparray_atx_recessive_meanimpute!, :_snparray_atx_recessive_meanimpute_rem!,
            :((Aij == 3) * v[4 * (l - 1) + p] + (Aij == 1) * μ[i] * v[4 * (l - 1) + p]))
    ]
    @eval begin
        function $(_ftn_rem!)(out, s, v, μ)
            maxp = length(v)
            l = 1
            @avx for i in eachindex(out)
                block = s[1, i]
                for p in 1:maxp
                    Aij = (block >> (2 * (p - 1))) & 3
                    out[i] += $expr
                end
            end
        end

        function $(_ftn!)(out, s, v, rows, cols, μ)
            k = rows >> 2
            rem = rows & 3

            if k ≥ 1 # avoid avxing over empty collection
                @avx for i ∈ 1:cols
                    for l in 1:k
                        block = s[l, i]

                        for p in 1:4
                            Aij = (block >> (2 * (p - 1))) & 3
                            out[i] += $expr
                        end
                    end
                end
            end
            if rem != 0
                $(_ftn_rem!)(out, @view(s[k + 1:k + 1, :]), @view(v[4k + 1:end]), μ)
            end
            nothing
        end
    end
end

for (_ftn!, _ftn_rem!, expr) in [
        (:_snparray_AX_additive!, :_snparray_AX_additive_rem!,
            :(((Aij >= 2) + (Aij == 3) - μ[j]) * σinv[j] * V[j, c])),
        (:_snparray_AX_dominant!, :_snparray_AX_dominant_rem!,
            :(((Aij >= 2) - μ[j]) * σinv[j] * V[j, c])),
        (:_snparray_AX_recessive!, :_snparray_AX_recessive_rem!,
            :(((Aij == 3) - μ[j]) * σinv[j] * V[j, c])),
        (:_snparray_AX_additive_meanimpute!, :_snparray_AX_additive_meanimpute_rem!,
            :(((((Aij >= 2) + (Aij == 3) - μ[j]) + (Aij == 1) * μimpute[j]) * σinv[j] * V[j, c]))),
        (:_snparray_AX_dominant_meanimpute!, :_snparray_AX_dominant_meanimpute_rem!,
            :(((Aij >= 2) - μ[j]) * σinv[j] * V[j, c] + (Aij == 1) * μimpute[j] * σinv[j] * V[j, c])),
            (:_snparray_AX_recessive_meanimpute!, :_snparray_AX_recessive_meanimpute_rem!,
            :(((Aij == 3) - μ[j]) * V[j, c] * σinv[j] + (Aij == 1) * μimpute[j] * σinv[j] * V[j, c]))
    ]
    @eval begin
        function ($_ftn_rem!)(out, s, V, μ, μimpute, σinv, c)
            maxp = size(out, 1)
            @avx for j in 1:size(V, 1)
                block = s[1, j]
                for p in 1:maxp
                    Aij = (block >> (2 * (p - 1))) & 3
                    out[p, c] += $expr
                end
            end
        end

        function ($_ftn!)(out, s, V, srows, scols, Vcols, μ, μimpute, σinv)
            k = srows >> 2 # fast div(srows, 4)
            rem = srows & 3 # fast rem(srows, 4)

            # compute out[i, c] = s[i, j] * V[j, c] for j in 1:n
            if k ≥ 1 # avoid avxing over empty collection
                @avx for c in 1:Vcols
                    for j in 1:scols
                        for l in 1:k
                            block = s[l, j]
                            for p in 1:4
                                Aij = (block >> (2 * (p - 1))) & 3
                                out[4*(l - 1) + p, c] += $expr
                            end
                        end
                    end
                end
            end
            if rem != 0
                for c in 1:Vcols
                    ($_ftn_rem!)(@view(out[4k + 1:end, :]), @view(s[k + 1:k + 1, :]),
                        V, μ, μimpute, σinv, c)
                end
            end
            nothing
        end
    end
end

for (_ftn!, _ftn_rem!, expr) in [
    (:_snparray_AtX_additive!, :_snparray_AtX_additive_rem!,
        :(((Aij >= 2) + (Aij == 3) - μ[i]) * σinv[i] * V[4 * (l - 1) + p, c])),
    (:_snparray_AtX_dominant!, :_snparray_AtX_dominant_rem!,
        :(((Aij >= 2) - μ[i]) * σinv[i] * V[4 * (l - 1) + p, c])),
    (:_snparray_AtX_recessive!, :_snparray_AtX_recessive_rem!,
        :(((Aij == 3) - μ[i]) * σinv[i] * V[4 * (l - 1) + p, c])),
    (:_snparray_AtX_additive_meanimpute!, :_snparray_AtX_additive_meanimpute_rem!,
        :(((((Aij >= 2) + (Aij == 3) - μ[i]) + (Aij == 1) * μimpute[i]) * σinv[i] * V[4 * (l - 1) + p, c]))),
    (:_snparray_AtX_dominant_meanimpute!, :_snparray_AtX_dominant_meanimpute_rem!,
        :(((Aij >= 2) - μ[i]) * σinv[i] * V[4 * (l - 1) + p, c] + (Aij == 1) * μimpute[i] * σinv[i] * V[4 * (l - 1) + p, c])),
    (:_snparray_AtX_recessive_meanimpute!, :_snparray_AtX_recessive_meanimpute_rem!,
        :(((Aij == 3) - μ[i]) * σinv[i] * V[4 * (l - 1) + p, c] + (Aij == 1) * μimpute[i] * σinv[i] * V[4 * (l - 1) + p, c]))
]
    @eval begin
        function ($_ftn_rem!)(out, s, V, μ, μimpute, σinv, c)
            maxp = size(V, 1)
            l = 1
            @avx for i in 1:size(out, 1)
                block = s[1, i]
                for p in 1:maxp
                    Aij = (block >> (2 * (p - 1))) & 3
                    out[i, c] += $expr
                end
            end
        end

        function ($_ftn!)(out, s, V, srows, scols, Vcols, μ, μimpute, σinv)
            k = srows >> 2 # fast div(srows, 4)
            rem = srows & 3 # fast rem(srows, 4)

            if k ≥ 1 # avoid avxing over empty collection
                @avx for c in 1:Vcols
                    for i in 1:scols
                        for l in 1:k
                            block = s[l, i]
                            for p in 1:4
                                Aij = (block >> (2 * (p - 1))) & 3
                                out[i, c] += $expr
                            end
                        end
                    end
                end
            end
            if rem != 0
                for c in 1:Vcols
                    ($_ftn_rem!)(out, @view(s[k + 1:k + 1, :]),
                        @view(V[4k + 1:end, :]), μ, μimpute, σinv, c)
                end
            end
            nothing
        end
    end
end

"""
    Base.copyto!(v, s)

Copy SnpLinAlg `s` to numeric vector or matrix `v`. If `s` is centered/scaled,
`v` will be centered/scaled using precomputed column mean `s.μ` and inverse std
`s.σinv`.
"""
function Base.copyto!(
    v::AbstractVecOrMat{T},
    s::AbstractSnpLinAlg
    ) where T <: AbstractFloat
    m, n = size(s, 1), size(s, 2)
    @inbounds for j in 1:n
        @simd for i in 1:m
            v[i, j] = s[i, j]
        end
    end
    return v
end

"""
    Base.convert(t, s)

Convert a AbstractSnpLinAlg `s` to a numeric vector or matrix of same shape as `s`.
If `s` is centered/scaled, `v` will be centered/scaled using precomputed column
mean `s.μ` and inverse std `s.σinv`.

# Arguments
- `t::Type{AbstractVecOrMat{T}}`: Vector or matrix type.
"""
Base.convert(::Type{T}, s::AbstractSnpLinAlg) where T <: Array = T(s)
Array{T,N}(s::AbstractSnpLinAlg) where {T,N} =
    copyto!(Array{T,N}(undef, size(s)), s)


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

if length(ARGS)!=1
    throw(ValueError(length(ARGS), "wrong number of arguments"))
end

c = Coverage.process_file(ARGS[1])
lines = split(c.source, '\n')
for (c,l) in zip(c.coverage, lines)
    if c==nothing
        print("-\t")
    elseif c==0
        print("0 ***\t")
    else
        print("$c\t")
    end
    println(l)
end


text:
"""
The base type of a fractional integral sense, in FractionalCalculus.jl, all of the fractional integral senses belong to ```FracIntAlg```
"""
abstract type FracIntAlg end

include("RL.jl")
include("Hadamard.jl")

include("SymbolicInt.jl")

"""
    fracint(f, α, point, h, FracIntAlg())

The general function used for fractional integral computing.

### Example

```julia-repl
julia> fracint(x->x^2, 0.5, 1, 0.001, RLInt_Approx())
0.6017877646029814
```

There are many algorithms can be used to computing fractional integral, please refer to our manual for more details: [manual](http://scifracx.org/FractionalCalculus.jl/dev/Integral/integralapi/).
"""
function fracint(f, α, point, h, ::FracIntAlg)
end

## Macros for convenient computing.
"""
    @fracint(f, α, point)

Return the α-order integral of ``f`` at specific point.

```julia-repl
julia> @fracint(x->x, 0.5, 1)
0.7522525439593486
```
"""
macro fracint(f, α, point)
    return :(fracint($f, $α, $point, 0.0001, RLIntApprox()))
end

"""
    @semifracint(f, point)

Return the semi-integral of ``f`` at specific point.

```julia-repl
julia> @semifracint(x->x, 1)
0.7522525439593486
```
"""
macro semifracint(f, point)
    return :(fracint($f, 0.5, $point, 0.0001, RLIntApprox()))
end


text:
export View

using Glob

struct Criteria
    instrument_type::Union{DataType,UnionAll,Nothing}
    instrument_name::Union{Glob.FilenameMatch{String},Nothing}
    meter_name::Union{String,Nothing}
    meter_version::Union{VersionNumber,Nothing}
    meter_schema_url::Union{String,Nothing}
end

function Base.occursin(ins::AbstractInstrument, c::Criteria)
    if !isnothing(c.instrument_type)
        if !(ins isa c.instrument_type)
            return false
        end
    end

    if !isnothing(c.instrument_name)
        if !occursin(c.instrument_name, ins.name)
            return false
        end
    end

    if !isnothing(c.meter_name)
        if ins.meter.name != c.meter_name
            return false
        end
    end

    if !isnothing(c.meter_version)
        if ins.meter.version != c.meter_version
            return false
        end
    end

    if !isnothing(c.meter_schema_url)
        if ins.meter.schema_url != c.meter_schema_url
            return false
        end
    end

    return true
end

struct View{A}
    name::Union{String,Nothing}
    description::Union{String,Nothing}
    criteria::Criteria
    attribute_keys::Union{Tuple{Vararg{String}},Nothing}
    extra_dimensions::StaticAttrs
    aggregation::A
end

function View(
    name = nothing;
    description = nothing,
    attribute_keys = nothing,
    extra_dimensions = StaticAttrs(),
    aggregation = nothing,
    # criteria
    instrument_name = nothing,
    instrument_type = nothing,
    meter_name = nothing,
    meter_version = nothing,
    meter_schema_url = nothing,
)
    if !isnothing(name)
        if isnothing(instrument_name)
            instrument_name = name
            if occursin(r"[*?\[\]]", name)
                # name contains wildcard
                name = nothing  # so that the name will be instrument name automatically
            end
        end
    end

    something(instrument_name, instrument_type, meter_name, meter_version, meter_schema_url)
    criteria = Criteria(
        instrument_type,
        isnothing(instrument_name) ? nothing : Glob.FilenameMatch(instrument_name),
        meter_name,
        meter_version,
        meter_schema_url,
    )
    View(name, description, criteria, attribute_keys, extra_dimensions, aggregation)
end


text:
include("test_params.jl")
using SymbolicRegression, SymbolicUtils, Test, Random, ForwardDiff
using SymbolicRegression:
    Options, string_tree, eval_tree_array, Dataset, differentiable_eval_tree_array
using SymbolicRegression: print_tree, pow, eval_loss, score_func, Node
using SymbolicRegression:
    plus,
    sub,
    mult,
    square,
    cube,
    div,
    log_abs,
    log2_abs,
    log10_abs,
    sqrt_abs,
    acosh_abs,
    neg,
    greater,
    greater,
    relu,
    logical_or,
    logical_and,
    gamma
using SymbolicRegression: node_to_symbolic, symbolic_to_node
using SymbolicRegression: check_constraints, loss

x1 = 2.0
# Initialize functions in Base....
for unaop in [cos, exp, log_abs, log2_abs, log10_abs, relu, gamma, acosh_abs]
    for binop in [sub]
        function make_options(; kw...)
            return Options(;
                default_params...,
                binary_operators=(+, *, ^, /, binop),
                unary_operators=(unaop,),
                npopulations=4,
                verbosity=(unaop == gamma) ? 0 : Int(1e9),
                kw...,
            )
        end
        make_options()

        # for unaop in
        f_true = (x,) -> binop((3.0 * unaop(x))^2.0, -1.2)

        # binop at outside:
        tree = Node(5, (Node(3.0) * Node(1, Node("x1")))^2.0, -1.2)
        tree_bad = Node(5, (Node(3.0) * Node(1, Node("x1")))^2.1, -1.3)
        n = count_nodes(tree)

        true_result = f_true(x1)

        result = eval(Meta.parse(string_tree(tree, make_options())))

        # Test Basics
        @test n == 8
        @test result == true_result

        types_to_test = [Float32, Float64, BigFloat]
        if unaop == cos
            # Other unary operators produce numbers too large
            # to do meaningful tests
            types_to_test = [Float16, types_to_test...]
        end
        for T in types_to_test
            if T == Float16 || unaop == gamma
                zero_tolerance = 3e-2
            else
                zero_tolerance = 1e-6
            end

            Random.seed!(0)
            N = 100
            X = T.(randn(MersenneTwister(0), Float64, 5, N) / 3)
            X = X + sign.(X) * T(0.1)
            y = T.(f_true.(X[1, :]))
            dataset = Dataset(X, y)
            test_y, complete = eval_tree_array(tree, X, make_options())
            test_y2, complete2 = differentiable_eval_tree_array(tree, X, make_options())

            # Test Evaluation
            @test complete == true
            @test all(abs.(test_y .- y) / N .< zero_tolerance)
            @test complete2 == true
            @test all(abs.(test_y2 .- y) / N .< zero_tolerance)

            # Test loss:
            @test abs(eval_loss(tree, dataset, make_options())) < zero_tolerance
            @test eval_loss(tree, dataset, make_options()) ==
                score_func(dataset, one(T), tree, make_options())[2]

            #Test Scoring
            @test abs(score_func(dataset, one(T), tree, make_options(; parsimony=0.0))[1]) <
                zero_tolerance
            @test score_func(dataset, one(T), tree, make_options(; parsimony=1.0))[1] > 1.0
            @test score_func(dataset, one(T), tree, make_options())[1] <
                score_func(dataset, one(T), tree_bad, make_options())[1]
            @test score_func(dataset, one(T) * 10, tree_bad, make_options())[1] <
                score_func(dataset, one(T), tree_bad, make_options())[1]

            # Test gradients:
            df_true = x -> ForwardDiff.derivative(f_true, x)
            dy = T.(df_true.(X[1, :]))
            test_dy = ForwardDiff.gradient(
                _x -> sum(differentiable_eval_tree_array(tree, _x, make_options())[1]), X
            )
            test_dy = test_dy[1, 1:end]
            @test all(abs.(test_dy .- dy) / N .< zero_tolerance)
        end
    end
end

# Generic operator tests
types_to_test = [Float16, Float32, Float64, BigFloat]
for T in types_to_test
    val = T(0.5)
    val2 = T(3.2)
    @test sqrt_abs(val) == sqrt_abs(-val)
    @test abs(log_abs(-val) - log(val)) < 1e-6
    @test abs(log2_abs(-val) - log2(val)) < 1e-6
    @test abs(log10_abs(-val) - log10(val)) < 1e-6
    @test neg(-val) == val
    @test sqrt_abs(val) == sqrt(val)
    @test mult(val, val2) == val * val2
    @test plus(val, val2) == val + val2
    @test sub(val, val2) == val - val2
    @test square(val) == val * val
    @test cube(val) == val * val * val
    @test div(val, val2) == val / val2
    @test greater(val, val2) == T(0.0)
    @test greater(val2, val) == T(1.0)
    @test relu(-val) == T(0.0)
    @test relu(val) == val
    @test logical_or(val, val2) == T(1.0)
    @test logical_or(T(0.0), val2) == T(1.0)
    @test logical_and(T(0.0), val2) == T(0.0)
end

# Test SymbolicUtils interface
_inv(x) = 1 / x
options = Options(;
    default_params...,
    binary_operators=(+, *, ^, /, greater),
    unary_operators=(_inv,),
    constraints=(_inv => 4,),
    npopulations=4,
)
tree = Node(5, (Node(3.0) * Node(1, Node("x1")))^2.0, -1.2)

eqn = node_to_symbolic(tree, options; varMap=["energy"], index_functions=true)
tree2 = symbolic_to_node(eqn, options; varMap=["energy"])

@test string_tree(tree, options) == string_tree(tree2, options)

# Test constraint-checking interface
tree = Node(5, (Node(3.0) * Node(1, Node("x1")))^2.0, -1.2)
violating_tree = Node(1, tree)

@test check_constraints(tree, options) == true
@test check_constraints(violating_tree, options) == false

# Test different loss functions
customloss(x, y) = abs(x - y)^2.5
customloss(x, y, w) = w * (abs(x - y)^2.5)
testl1(x, y) = abs(x - y)
testl1(x, y, w) = abs(x - y) * w

for (loss_fnc, evaluator) in [(L1DistLoss(), testl1), (customloss, customloss)]
    local options = Options(;
        default_params...,
        binary_operators=(+, *, -, /),
        unary_operators=(cos, exp),
        npopulations=4,
        loss=loss_fnc,
    )
    x = randn(MersenneTwister(0), Float32, 100)
    y = randn(MersenneTwister(1), Float32, 100)
    w = abs.(randn(MersenneTwister(2), Float32, 100))
    @test abs(loss(x, y, options) - sum(evaluator.(x, y)) / length(x)) < 1e-6
    @test abs(loss(x, y, w, options) - sum(evaluator.(x, y, w)) / sum(w)) < 1e-6
end

# Test derivatives
include("test_derivatives.jl")

# Test simplification:
include("test_simplification.jl")

# Test `print`:
include("test_print.jl")

include("test_evaluation.jl")

include("test_prob_pick_first.jl")

include("test_crossover.jl")

include("test_nan_detection.jl")

include("test_constraints.jl")

include("test_complexity.jl")


text:
using StenoGraphs
@testset "Conversion" begin
    @test convert(Node,:a) == Node(:a)
end

struct Start{S <: Number} <: EdgeModifier
    s::S
end

struct Label <: NodeModifier
    l
end

@testset "Promotion" begin
    @test eltype([Node(:a), :b]) <: Node
    @test eltype([Node(:a), :b, Start(1) * Node(:c)]) <: AbstractNode
    @test eltype([:b, Start(1) * Node(:c)]) <: AbstractNode
    @test eltype([Start(1) * Node(:b), Start(1) * Node(:c)]) <: MetaNode

    @test eltype(promote(Node(:a), ModifiedNode(Node(:b), Label("hi")))) <: ModifiedNode
    ns = Node.([:a, :b]) .^ Ref(Label("test"))
    ns[1] = Node(:a)
    @test eltype(ns) <: ModifiedNode

    @test eltype([ModifiedEdge(Edge(Node(:a), Node(:b)), Start(3.0)), Edge(Node(:a), Node(:c))]) <: ModifiedEdge
end


text:
@testset "GraphDataFrameBridge" begin

    # Test data from constructor
    df = DataFrame(Dict("start" => ["a", "b", "a", "d"],
                        "finish" => ["b", "c", "e", "e"],
                        "weights" => 1:4,
                        "extras" => 5:8))

    @inferred(MetaGraph(df, :start, :finish))
    mg = MetaGraph(df, :start, :finish)
    @test get_prop(mg, 1, :name) == "a"
    @test get_prop(mg, 2, :name) == "b"

    @test neighbors(mg, 1)[1] == 2
    @test length(neighbors(mg, 1)) == 2

    @test neighbors(mg, 2)[1] == 1
    @test neighbors(mg, 2)[2] == 3
    @test get_prop(mg, 3, :name) == "c"
    @test length(neighbors(mg, 2)) == 2

    @inferred(MetaDiGraph(df, :start, :finish))
    mg = MetaDiGraph(df, :start, :finish)
    @test neighbors(mg, 2)[1] == 3
    @test length(neighbors(mg, 2)) == 1

    @inferred (MetaGraph(df, :start, :finish,
                         weight=:weights,
                         edge_attributes=:extras))
    mg = MetaGraph(df, :start, :finish,
                   weight=:weights,
                   edge_attributes=:extras)
    @test length(neighbors(mg, 2)) == 2
    @test get_prop(mg, Edge(1, 2), :weight) == 1
    @test get_prop(mg, Edge(4, 5), :weight) == 4
    @test get_prop(mg, Edge(4, 5), :extras) == 8

    # Test with different column names to ensure nothing name sensitive
    df2 = DataFrame(Dict("alice" => ["a", "b", "a", "d"],
                        "bob" => ["b", "c", "e", "e"],
                        "weightz" => 1:4,
                        "extraz" => 5:8))

    mg = MetaGraph(df2, :alice, :bob,
             weight=:weightz,
             edge_attributes=:extraz)
    @test length(neighbors(mg, 2)) == 2
    @test get_prop(mg, Edge(1, 2), :weight) == 1
    @test get_prop(mg, Edge(4, 5), :weight) == 4
    @test get_prop(mg, Edge(4, 5), :extraz) == 8

    # Test name column is indexed (Issue #9)
    @test mg["a", :name] == 1
    @test mg["b", :name] == 2

end