text:
module JointMarginalsTest

using Test
using ForneyLab
using ForneyLab: generateId, addNode!, associate!, inferMarginalRule, isApplicable, Cluster, Product, setTargets!, outboundType, messagePassingSchedule

# Integration helper
mutable struct MockNode <: FactorNode
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Int,Interface}

    function MockNode(vars::Vector{Variable}; id=generateId(MockNode))
        n_interfaces = length(vars)
        self = new(id, Vector{Interface}(undef, n_interfaces), Dict{Int,Interface}())
        addNode!(currentGraph(), self)

        for idx = 1:n_interfaces
            self.i[idx] = self.interfaces[idx] = associate!(Interface(self), vars[idx])
        end

        return self
    end
end

@marginalRule(:node_type     => MockNode,
              :inbound_types => (Message, Message, ProbabilityDistribution),
              :name          => MMockMMD)

@testset "@marginalRule" begin
    @test MMockMMD <: MarginalRule{MockNode}
    @test isApplicable(MMockMMD, [Message, Message, ProbabilityDistribution])
    @test !isApplicable(MMockMMD, [Message, Message, ProbabilityDistribution, ProbabilityDistribution])
end

@testset "inferMarginalRule" begin
    FactorGraph()
    nd = MockNode([constant(0.0), constant(0.0), constant(0.0)])
    cluster = Cluster(nd, [nd.i[1].edge, nd.i[2].edge])

    @test inferMarginalRule(cluster, [Message, Message, ProbabilityDistribution]) == MMockMMD
end

mutable struct SPMockNode <: SumProductRule{MockNode} end
ForneyLab.outboundType(SPMockNode) = Message
ForneyLab.isApplicable(SPMockNode, input_types::Vector{<:Type}) = true

mutable struct SVBMockNode <: StructuredVariationalRule{MockNode} end
ForneyLab.outboundType(SVBMockNode) = Message
ForneyLab.isApplicable(SVBMockNode, input_types::Vector{<:Type}) = true

@testset "marginalTable" begin
    fg = FactorGraph()
    v1 = Variable()
    v2 = Variable()
    v3 = constant(0.0)
    nd1 = MockNode([v1])
    nd2 = MockNode([v2])
    nd3 = MockNode([v1, v2, v3])

    pfz = PosteriorFactorization()
    pf_12 = PosteriorFactor([v1, v2])
    setTargets!(pf_12, pfz, external_targets=true)

    pf_12.schedule = messagePassingSchedule(pf_12)
    marginal_table = marginalTable(pf_12)

    @test length(marginal_table) == 1
    @test marginal_table[1].target == first(pf_12.target_clusters)
    @test marginal_table[1].interfaces[1] == nd1.i[1]
    @test marginal_table[1].interfaces[2] == nd2.i[1]
    @test marginal_table[1].marginal_update_rule == MMockMMD
end

end # module


text:
module Tests

include("../ntl.jl")

import Random
using LinearAlgebra

function changepoint_gibbs_test(;n=10, iterations=10)
    Random.seed!(1)
    data_covariance = Matrix{Float64}(0.1I, 2, 2)
    prior_covariance = Matrix{Float64}(I, 2, 2)
    prior_mean = Vector{Float64}(zeros(2))
    data_parameters = Ntl.Models.GaussianParameters(data_covariance, prior_mean, prior_covariance)

    psi_prior = Vector{Float64}([1, 1])
    phi_prior = Vector{Float64}([1, 1])
    geometric_arrival = Ntl.Models.GeometricArrivals(phi_prior)
    ntl_cluster_parameters = Ntl.Models.NtlParameters(psi_prior, geometric_arrival)

    changepoint_model = Ntl.Models.Changepoint(ntl_cluster_parameters, data_parameters)

    changepoint = Ntl.Generate.generate(changepoint_model, n=n)
    data = Matrix(transpose(changepoint[:, 2:end]))

    gibbs_sampler = Ntl.Samplers.GibbsSampler(num_iterations=iterations)
    results = Ntl.Fitter.fit(data, changepoint_model, gibbs_sampler)
end

function changepoint_smc_test(;num_particles=10, ess_threshold=0.5)
    Random.seed!(1)
    data_covariance = Matrix{Float64}(0.1I, 2, 2)
    prior_covariance = Matrix{Float64}(I, 2, 2)
    prior_mean = Vector{Float64}(zeros(2))
    data_parameters = Ntl.Models.GaussianParameters(data_covariance, prior_mean, prior_covariance)

    psi_prior = Vector{Float64}([1, 1])
    phi_prior = Vector{Float64}([1, 1])
    geometric_arrival = Ntl.Models.GeometricArrivals(phi_prior)
    ntl_cluster_parameters = Ntl.Models.NtlParameters(psi_prior, geometric_arrival)

    changepoint_model = Ntl.Models.Changepoint(ntl_cluster_parameters, data_parameters)

    changepoint = Ntl.Generate.generate(changepoint_model, n=100)
    data = Matrix(transpose(changepoint[:, 2:end]))
    smc = Ntl.Samplers.SequentialMonteCarlo(num_particles=num_particles, ess_threshold=ess_threshold)
    results = Ntl.Fitter.fit(data, changepoint_model, smc)
end

function dp_hmm_test(;n=10, iterations=10)
    Random.seed!(1)

    data_covariance = Matrix{Float64}(0.1I, 2, 2)
    prior_covariance = Matrix{Float64}(I, 2, 2)
    prior_mean = Vector{Float64}(zeros(2))
    data_parameters = Ntl.Models.GaussianParameters(data_covariance, prior_mean, prior_covariance)

    psi_prior = Vector{Float64}([1, 1])
    phi_prior = Vector{Float64}([1, 1])
    geometric_arrival = Ntl.Models.GeometricArrivals(phi_prior)
    ntl_cluster_parameters = Ntl.Models.NtlParameters(psi_prior, geometric_arrival)

    mixture_model = Ntl.Models.Mixture(ntl_cluster_parameters, data_parameters)
    mixture = Ntl.Generate.generate(mixture_model, n=n)
    data = Matrix(transpose(mixture[:, 2:end]))

    dp_arrivals = Ntl.Models.PitmanYorArrivals()
    dp_parameters = Ntl.Models.BetaNtlParameters(0., dp_arrivals)
    gibbs_sampler = Ntl.Samplers.GibbsSampler(num_iterations=iterations)
    hmm_model = Ntl.Models.HiddenMarkovModel(dp_parameters, data_parameters)
    markov_chain = Ntl.Fitter.fit(data, hmm_model, gibbs_sampler)
end

function ntl_hmm_test(;n=10, iterations=10)
    Random.seed!(1)

    data_covariance = Matrix{Float64}(0.1I, 2, 2)
    prior_covariance = Matrix{Float64}(I, 2, 2)
    prior_mean = Vector{Float64}(zeros(2))
    data_parameters = Ntl.Models.GaussianParameters(data_covariance, prior_mean, prior_covariance)

    psi_prior = Vector{Float64}([1, 1])
    phi_prior = Vector{Float64}([1, 1])
    geometric_arrival = Ntl.Models.GeometricArrivals(phi_prior)
    ntl_cluster_parameters = Ntl.Models.NtlParameters(psi_prior, geometric_arrival)

    mixture_model = Ntl.Models.Mixture(ntl_cluster_parameters, data_parameters)
    mixture = Ntl.Generate.generate(mixture_model, n=n)
    data = Matrix(transpose(mixture[:, 2:end]))

    gibbs_sampler = Ntl.Samplers.GibbsSampler(num_iterations=iterations)
    hmm_model = Ntl.Models.HiddenMarkovModel(ntl_cluster_parameters, data_parameters)
    markov_chain = Ntl.Fitter.fit(data, hmm_model, gibbs_sampler)
end

function ntl_mixture_test(;n=10, iterations=10)
    Random.seed!(1)

    scale_matrix = Matrix{Float64}(I, 2, 2)
    scale = 10000
    dof = 2
    prior_mean = Vector{Float64}(zeros(2))
    data_parameters = Ntl.Models.GaussianWishartParameters(prior_mean, scale_matrix, scale, dof)

    psi_prior = Vector{Float64}([1, 1])
    phi_prior = Vector{Float64}([1, 1])
    geometric_arrival = Ntl.Models.GeometricArrivals(phi_prior)
    ntl_cluster_parameters = Ntl.Models.NtlParameters(psi_prior, geometric_arrival)

    mixture_model = Ntl.Models.Mixture(ntl_cluster_parameters, data_parameters)
    mixture = Ntl.Generate.generate(mixture_model, n=n)
    data = Matrix(transpose(mixture[:, 2:end]))

    gibbs_sampler = Ntl.Samplers.GibbsSampler(num_iterations=iterations)
    model = Ntl.Models.Mixture(ntl_cluster_parameters, data_parameters)
    markov_chain = Ntl.Fitter.fit(data, model, gibbs_sampler)
end

end

ext: jl

max_stars_repo_name: realseanla/ntl-stick-breaking-julia

max_issues_repo_name: realseanla/ntl-stick-breaking-julia

max_forks_repo_name: realseanla/ntl-stick-breaking-julia


text:
using JSON
using GasModels
using Base.Filesystem

get_simulation_params(; T::Float64 = 288.70599999999996) =
    Dict{String,Any}(
    "Temperature (K):" => T,
    "Gas specific gravity (G):" => 0.6,
    "Specific heat capacity ratio" => 1.4,
    "units (SI = 0, standard = 1)" => 0.0
)

function generate_steady_sim_inputs_from_gaslib(file::AbstractString,
    input_folder::AbstractString, output_folder::AbstractString;
    compressor_ratio::Float64 = 1.5, control_valve_ratio::Float64 = 1.0)

    data = parse_file(input_folder * file)

    (!isdir(output_folder)) && (mkdir(output_folder))

    network_data = Dict{String,Any}(
        "nodes" => Dict{String,Any}(),
        "pipes" => Dict{String,Any}(),
        "compressors" => Dict{String,Any}(),
    )

    bc = Dict{String,Any}(
        "boundary_compressor" => Dict{String,Any}(),
        "boundary_nonslack_flow" => Dict{String,Any}(),
        "boundary_pslack" => Dict{String,Any}()
    )

    withdrawal = Dict{String,Any}()

    for (i, _) in data["junction"]
        withdrawal[i] = 0.0
    end

    for (_, receipt) in data["receipt"]
        node = string(receipt["junction_id"])
        withdrawal[node] -= receipt["injection_nominal"] * data["base_flow"]
    end

    for (_, delivery) in data["delivery"]
        node = string(delivery["junction_id"])
        withdrawal[node] += delivery["withdrawal_nominal"] * data["base_flow"]
    end

    withdrawal_values = [v for (k, v) in withdrawal]
    max_injection = minimum(withdrawal_values)
    nodes_with_max_injection = sort([k for (k,v) in withdrawal if v == max_injection])
    slack_node = nodes_with_max_injection[1]

    for (i, node) in data["junction"]
        network_data["nodes"][i] = Dict{String,Any}(
            "node_id" => parse(Int, i),
            "node_name" => "n" * i,
            "x_coord" => node["lat"],
            "y_coord" => node["lon"],
            "min_pressure" => node["p_min"] * data["base_pressure"],
            "max_pressure" => node["p_max"] * data["base_pressure"]
        )
        if i == string(slack_node)
            network_data["nodes"][i]["slack_bool"] = 1
            bc["boundary_pslack"][i] = 5000000
        else
            network_data["nodes"][i]["slack_bool"] = 0
            if !(withdrawal[i] == 0)
                bc["boundary_nonslack_flow"][i] = withdrawal[i]
            end
        end
    end

    for (i, pipe) in data["pipe"]
        network_data["pipes"][i] = Dict{String,Any}(
            "pipe_id" => parse(Int, i),
            "pipe_name" => "p" * i,
            "from_node" => pipe["fr_junction"],
            "to_node" => pipe["to_junction"],
            "diameter" => pipe["diameter"] * data["base_diameter"],
            "length" => pipe["length"] * data["base_length"],
            "friction_factor" => pipe["friction_factor"],
        )
    end

    for (i, compressor) in data["compressor"]
        network_data["compressors"][i] = Dict{String,Any}(
            "comp_id" => parse(Int, i),
            "comp_name" => "c" * i,
            "from_node" => compressor["fr_junction"],
            "to_node" => compressor["to_junction"]
        )
        bc["boundary_compressor"][i] = Dict{String,Any}(
            "control_type" => 0, "value" => compressor_ratio
        )
    end

    network_data["control_valves"] = Dict{String,Any}()
    bc["boundary_control_valve"] = Dict{String,Any}("on" => [], "off" => [])
    for (i, cv) in data["regulator"]
        network_data["control_valves"][i] = Dict{String,Any}(
        "control_valve_id" => parse(Int, i),
        "control_valve_name" => "cv" * i,
        "from_node" => cv["fr_junction"],
        "to_node" => cv["to_junction"],
        "bypass_required" => cv["bypass_required"]
        )
        push!(bc["boundary_control_valve"]["on"], parse(Int, i))
        bc["boundary_control_valve"][i] = Dict{String,Any}(
            "control_type" => 0, "value" => control_valve_ratio
        )
    end

    network_data["valves"] = Dict{String,Any}()
    bc["boundary_valve"] = Dict{String,Any}("on" => [], "off" => [])
    for (i, v) in data["valve"]
        network_data["valves"][i] = Dict{String,Any}(
            "valve_id" => parse(Int, i),
            "valve_name" => "v" * i,
            "from_node" => v["fr_junction"],
            "to_node" => v["to_junction"]
        )
        push!(bc["boundary_valve"]["on"], parse(Int, i))
    end

    network_data["resistors"] = Dict{String,Any}()
    network_data["loss_resistors"] = Dict{String,Any}()
    network_data["short_pipes"] = Dict{String,Any}()
    for (i, res) in data["resistor"]
        network_data["resistors"][i] = Dict{String,Any}(
            "resistor_id" => parse(Int, i),
            "resistor_name" => "r" * i,
            "from_node" => res["fr_junction"],
            "to_node" => res["to_junction"],
            "drag" => res["drag"],
            "diameter" => res["diameter"]
        )
    end

    for (i, res) in data["loss_resistor"]
        network_data["loss_resistors"][i] = Dict{String,Any}(
            "loss_resistor_id" => parse(Int, i),
            "loss_resistor_name" => "lr" * i,
            "from_node" => res["fr_junction"],
            "to_node" => res["to_junction"],
            "p_loss" => res["p_loss"] * data["base_pressure"]
        )
    end

    for (i, pipe) in data["short_pipe"]
        network_data["short_pipes"][i] = Dict{String,Any}(
            "short_pipe_id" => parse(Int, i),
            "short_pipe_name" => "sp" * i,
            "from_node" => pipe["fr_junction"],
            "to_node" => pipe["to_junction"]
        )
    end

    open(output_folder * "params.json", "w") do f
        JSON.print(f, Dict("simulation_params" => get_simulation_params()), 2)
    end

    open(output_folder * "network.json", "w") do f
        JSON.print(f, network_data, 2)
    end

    open(output_folder * "bc.json", "w") do f
        JSON.print(f, bc, 2)
    end

end

generate_steady_sim_inputs_from_gaslib("GasLib-11.zip",
    "./data/GasLib-benchmarks/", "./data/GasLib-11/")
generate_steady_sim_inputs_from_gaslib("GasLib-24.zip",
    "./data/GasLib-benchmarks/", "./data/GasLib-24/")
generate_steady_sim_inputs_from_gaslib("GasLib-40.zip",
    "./data/GasLib-benchmarks/", "./data/GasLib-40/")
generate_steady_sim_inputs_from_gaslib("GasLib-134.zip",
    "./data/GasLib-benchmarks/", "./data/GasLib-134/")
generate_steady_sim_inputs_from_gaslib("GasLib-135.zip",
    "./data/GasLib-benchmarks/", "./data/GasLib-135/")

max_stars_repo_name: kaarthiksundar/GasFlowRuns

max_issues_repo_name: kaarthiksundar/GasFlowRuns

max_forks_repo_name: kaarthiksundar/GasFlowRuns


text:
@testset "Fitting" begin
    game = Kuhn()

    sol = DeepCFRSolver(game)
    train!(sol, 100)

    p = 1
    net = sol.V[p]

    DeepCFR.initialize!(net, Flux.glorot_normal)

    dv_0 = DeepCFR.optimality_distance(net, sol.Mv[p])

    hist_glorot = DeepCFR.train_net_tracked!(
        net,
        sol.Mv[p],
        256,
        1000,
        ADAM(1e-3);
        show_progress = false
    )
    dv_f = DeepCFR.optimality_distance(net, sol.Mv[p])

    @test dv_0 > dv_f

    DeepCFR.initialize!(sol.Π, Flux.glorot_normal)
    dπ_0 = DeepCFR.optimality_distance(sol.Π, sol.Mπ)
    strat_hist = DeepCFR.train_net_tracked!(
        sol.Π,
        sol.Mπ,
        256,
        500,
        ADAM(1e-3);
        show_progress = false
    )

    dπ_f = DeepCFR.optimality_distance(sol.Π, sol.Mπ)

    @test dπ_0 > dπ_f
    @test last(hist_glorot) < first(hist_glorot)
end


text:
function solis_wets(problem;
  max_evals_per_dim = 1e5, ftol = 1e-7,
  max_success_steps = 5, max_fail_steps = 3,
  rnggen = (sigma) -> Normal(0, sigma))

  start_time = time()

  n = numdims(problem)
  ss = search_space(problem)
  num_successes = num_failures = 0
  sigma = sqrt(0.30 * minimum(diameters(ss)))
  diff = bias = zeros(n)
  rng = rnggen(sigma)
  max_evals = max_evals_per_dim * n

  xbest = rand_individual(ss)
  termination_reason = "unknown"

  a = BlackBoxOptim.TopListArchive(n, 10)
  fbest = eval1(xbest, problem)
  num_fevals = 1
  add_candidate!(a, fbest, xbest[:], num_fevals)

  while(true)
    # Terminate if fitness within ftol of known fmin
    if fitness_is_within_ftol(problem, ftol, fbest)
      termination_reason = "Within ftol"
      break
    end

    if num_fevals > max_evals
      termination_reason = "Max evals budget exceeded"
      break
    end

    diff = rand(rng, n)
    xplus = xbest + bias + diff
    fxplus = eval1(xplus, problem)
    num_fevals += 1

    if fxplus < fbest

      fbest = fxplus
      add_candidate!(a, fxplus, xplus[:], num_fevals)
      num_successes += 1
      num_failures = 0
      bias = 0.2 * bias + 0.4 * (diff + bias)
      xbest = xplus
      println( "$(num_fevals): New best, fitness = $(fbest), width = $(width_of_confidence_interval(a, 0.01))")

    else

      xminus = xbest - bias - diff
      fxminus = eval1(xminus, problem)
      num_fevals += 1

      if fxminus < fbest

        fbest = fxminus
        add_candidate!(a, fxminus, xminus[:], num_fevals)
        num_successes += 1
        num_failures = 0
        bias = bias - 0.4 * (diff + bias)
        xbest = xminus
        println( "$(num_fevals): New best, fitness = $(fbest), width = $(width_of_confidence_interval(a, 0.01))")

      else
        bias = bias / 2
        num_failures += 1
        num_successes = 0
      end

    end

    if num_successes > max_success_steps
      sigma = sigma * 2
      rng = rnggen(sigma)
      num_successes = 0
    elseif num_failures > max_fail_steps
      sigma = max(1e-20, sigma / 2)
      rng = rnggen(sigma)
      num_failures = 0
    end

  end

  return xbest, best_fitness(a), num_fevals, termination_reason, a
end

#type SolisWetsOptimizer <: Optimizer
#  # Current population: 1st column is current solution, 2nd columnn is current
#  # candidate to consider.
#  pop::Array{Float64, 2}
#  N::Int
#  bias::Array{Float64, 1}   # Bias
#  diff::Array{Float64, 1}   # Latest Diff vector.
#  sigma::Float64            # Measure of spread in random sampling. Note! Sigma used in stead of rho=sigma^2
#  rng::ContinuousMultivariateDistribution
#  ss::SearchSpace
#  num_successes::Int
#  num_failures::Int
#  state_is_plus::Bool       # Indicates which state we are in, i.e. if next candidate to generate is the plus or the minus candidate.
#
#  SolisWetsOptimizer(ss::SearchSpace, x = false; bias = false, sigma = false, rng = false) = begin
#    if rng == false
#      sigma = sigma ? sigma : (0.5*mean(diameters(ss)))
#      rng = Normal(0, sigma)
#    end
#
#    N = numdims(ss)
#
#    pop = zeros(N, 2)
#    if x != false
#      pop[:, 1] = x
#    else
#      pop[:, 1] = rand(rng, N)
#    end
#
#    bias = zeros(N)
#    diff = zeros(N)
#
#    new(pop, N, bias, bdiff, sigma, rng, ss, 0, 0, true)
#  end
#end
#
## Get a new set of candidate solutions.
#function ask(sw::SolisWetsOptimizer)
#  if sw.state_is_plus
#    sw.diff = rand(sw.rng, sw.N)
#    xdiff = sw.bias + sw.diff
#  else
#    xdiff = -sw.bias - sw.diff
#  end
#
#  sw.pop[:, 2] = sw.pop[:, 1] + xdiff
#
#  return sw.pop, [2], [1]
#end
#
#function increase_sigma(sw::SolisWetsOptimizer)
#  sw.sigma = sw.sigma * 2
#  sw.rng = Normal(0, sw.sigma)
#end
#
#function decrease_sigma(sw::SolisWetsOptimizer)
#  sw.sigma = sw.sigma / 2
#  sw.rng = Normal(0, sw.sigma)
#end
#
#function tell(sw::SolisWetsOptimizer, pop, ranks)
#  success = ranks[1] == 2  # Success if the 2nd individual is ranked first.
#
#  if success
#    sw.pop[:, 1] = sw.pop[:, 2]
#    sw.num_successes += 1
#    sw.num_failures = 0
#    sign = sw.state_is_plus ? (+1) : (-1)
#    sw.bias = 0.2 * sw.bias + sign * 0.4 * (sw.diff + sw.bias)
#    # With a success we always go next in the plus direction.
#    sw.state_is_plus = true
#  else
#    if sw.state_is_plus == false
#      sw.num_successes = 0
#      sw.num_failures += 1
#      sw.state_is_plus = true
#    end
#  end
#
#  if sw.num_successes > 5
#    increase_sigma(sw)
#    sw.num_successes = 0
#  elseif sw.num_failures > 3
#    decrease_sigma(sw)
#    sw.num_failures = 0
#  end
#end


text:
using Documenter
using StatsBase
using MixedModels

makedocs(
#    root = joinpath(dirname(pathof(MixedModels)), "..", "docs"),
    sitename = "MixedModels",
    doctest = true,
    pages = [
        "index.md",
        "constructors.md",
        "optimization.md",
        "GaussHermite.md",
        "prediction.md",
        "bootstrap.md",
        "rankdeficiency.md",
        "mime.md",
        "api.md",
    ],
)

deploydocs(;repo = "github.com/JuliaStats/MixedModels.jl.git", push_preview = true, devbranch = "main")


text:
# ------------------------------------------------------------------------------------------
# # Introduction to DataFrames
# **[Bogumił Kamiński](http://bogumilkaminski.pl/about/), Apr 21, 2018**
# ------------------------------------------------------------------------------------------

using DataFrames # load package
srand(1);

# ------------------------------------------------------------------------------------------
# ## Manipulating rows of DataFrame
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Reordering rows
# ------------------------------------------------------------------------------------------

x = DataFrame(id=1:10, x = rand(10), y = [zeros(5); ones(5)]) # and we hope that x[:x] is not sorted :)

issorted(x), issorted(x, :x) # check if a DataFrame or a subset of its columns is sorted

sort!(x, :x) # sort x in place

y = sort(x, :id) # new DataFrame

sort(x, (:y, :x), rev=(true, false)) # sort by two columns, first is decreasing, second is increasing

sort(x, (order(:y, rev=true), :x)) # the same as above

sort(x, (order(:y, rev=true), order(:x, by=v->-v))) # some more fancy sorting stuff

x[shuffle(1:10), :] # reorder rows (here randomly)

sort!(x, :id)
x[[1,10],:] = x[[10,1],:] # swap rows
x

x[1,:], x[10,:] = x[10,:], x[1,:] # and swap again
x

# ------------------------------------------------------------------------------------------
# ### Merging/adding rows
# ------------------------------------------------------------------------------------------

x = DataFrame(rand(3, 5))

[x; x] # merge by rows - data frames must have the same column names; the same is vcat

y = x[reverse(names(x))] # get y with other order of names

vcat(x, y) # we get what we want as vcat does column name matching

vcat(x, y[1:3]) # but column names must still match

append!(x, x) # the same but modifies x

append!(x, y) # here column names must match exactly

push!(x, 1:5) # add one row to x at the end; must give correct number of values and correct types
x

push!(x, Dict(:x1=> 11, :x2=> 12, :x3=> 13, :x4=> 14, :x5=> 15)) # also works with dictionaries
x

# ------------------------------------------------------------------------------------------
# ### Subsetting/removing rows
# ------------------------------------------------------------------------------------------

x = DataFrame(id=1:10, val='a':'j')

x[1:2, :] # by index

view(x, 1:2) # the same but a view

x[repmat([true, false], 5), :] # by Bool, exact length required

view(x, repmat([true, false], 5), :) # view again

deleterows!(x, 7) # delete one row

deleterows!(x, 6:7) # delete a collection of rows

x = DataFrame([1:4, 2:5, 3:6])

filter(r -> r[:x1] > 2.5, x) # create a new DataFrame where filtering function operates on DataFrameRow

# in place modification of x, an example with do-block syntax
filter!(x) do r
    if r[:x1] > 2.5
        return r[:x2] < 4.5
    end
    r[:x3] < 3.5
end

# ------------------------------------------------------------------------------------------
# ### Deduplicating
# ------------------------------------------------------------------------------------------

x = DataFrame(A=[1,2], B=["x","y"])
append!(x, x)
x[:C] = 1:4
x

unique(x, [1,2]) # get first unique rows for given index

unique(x) # now we look at whole rows

nonunique(x, :A) # get indicators of non-unique rows

unique!(x, :B) # modify x in place

# ------------------------------------------------------------------------------------------
# ### Extracting one row from `DataFrame` into a vector
# ------------------------------------------------------------------------------------------

x = DataFrame(x=[1,missing,2], y=["a", "b", missing], z=[true,false,true])

cols = [:x, :y]
[x[1, col] for col in cols] # subset of columns

[[x[i, col] for col in names(x)] for i in 1:nrow(x)] # vector of vectors, each entry contains one full row of x

Tuple(x[1, col] for col in cols) # similar construct for Tuples, when ported to Julia 0.7 NamedTuples will be added

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using ObjectDetectors
using GeometryBasics

@testset "BoundingBoxDetection" begin
    core_args = ("class", 0.7, 256, 256, MACHINE, "detector")

    box1 = BoundingBoxDetection(Point2(0.0, 1.0), 3.0, 6.0, core_args...)
    @test centroid(box1) == [1.5, 4.0]

    box2 = BoundingBoxDetection(Point2(1.0, 2.0), Point2(3.0, 4.0), core_args...)
    @test centroid(box2) == [2.0, 3.0]

    box3 = BoundingBoxDetection([Point2(0.0, 0.0), Point2(2.0, 0.0), Point2(2.0, 2.0), Point2(0.0, 2.0)], core_args...)
    @test centroid(box3) == [1.0, 1.0]

    box4 = BoundingBoxDetection(Rect2(0.0, 0.0, 1.0, 1.0), core_args...)
    @test centroid(box4) == [0.5, 0.5]

    box5 = create_with_center(Point2(3.0, 3.0), 5.0, 3.0, core_args...)
    @test centroid(box5) == [3.0, 3.0]
end


text:
using CamiXon
using IntervalSets
using LinearAlgebra
using Test

@testset "CamiXon.jl" begin
    @test sup(-5//2) == "⁻⁵ᐟ²"
    @test sub(-5//2) == "₋₅⸝₂"
    @test frac(-5//2) == "-⁵/₂"
    @test mendeleev(11) == ("Sodium", "Na")
    @test Atom("Helium ion", "⁴Heᐩ", 2, 1, 4.0026, 1//2, 0.0) == Atom("Helium ion", "⁴Heᐩ", 2, 1, 4.0026, 1//2, 0.0)
    @test createAtom(1,0,1.008,1//2,5.59) == Atom("Hydrogen", "¹H", 1, 0, 1.008, 1//2, 5.59)
    @test Term("1s ²S₁⸝₂", 1, 0, 1//2, 0, 1//2) == Term("1s ²S₁⸝₂", 1, 0, 1//2, 0, 1//2)
    @test createTerm(1,0,1//2,0,1//2) == Term("1s ²S₁⸝₂", 1, 0, 1//2, 0, 1//2)
    @test Grid(2, 0.01, 0.001, [0.0, 1.005], [1.0e-5, 1.01005]).r == [0.0, 1.005]
    @test gridfunction(10,0.01,0.001) == 9.41742836956125e-5
    @test createGrid(3; h=0.01, r0=0.001).r == [0.0, 1.0050167084168057e-5, 2.020134002675555e-5]
    @test find_all([:📑, :📌,:📢,:📌,:📞]) == [[1], [2, 4], [3], [5]]
    @test find_all([:📑, :📌,:📢,:📌,:📞]; count=true) == [1, 2, 1, 1]
    @test find_all([:📑, :📌,:📢,:📌,:📞], :📌) == [[2, 4]]
    @test find_all([:📑, :📌,:📢,:📌,:📞], :📌; count=true) == [2]
    @test find_all("aβcβd") == [[1], [2, 4], [3], [5]]
    @test find_all("aβcβd"; count=true) == [1, 2, 1, 1]
    @test find_all("aβcβd", 'β') == [[2, 4]]
    @test find_all("aβcβd", 'β'; count=true) == [2]
    @test find_first([:📑,:📌,:📢,:📌,:📞]) == [1,2,3,5]
    @test find_first([:📑,:📌,:📢,:📌,:📞]; dict=true) == [:📑 => 1,:📌 => 2,:📢 => 3,:📞 => 5]
    @test find_first([:📑,:📌,:📢,:📌,:📞], :📌) == [2]
    @test find_first([:📑,:📌,:📢,:📌,:📞], :📌; dict=true) == [:📌 => 2]
    @test find_first([:📑,:📌,:📢,:📌,:📞]) == find_first([1,2,3,2,5]) == find_first("aβcβd")
    @test find_last([:📑,:📌,:📢,:📌,:📞]) == [1,4,3,5]
    @test find_last([:📑,:📌,:📢,:📌,:📞]; dict=true) == [:📑 => 1,:📌 => 4,:📢 => 3,:📞 => 5]
    @test find_last([:📑,:📌,:📢,:📌,:📞], :📌) == [4]
    @test find_last([:📑,:📌,:📢,:📌,:📞], :📌; dict=true) == [:📌 => 4]
    @test find_last([:📑,:📌,:📢,:📌,:📞]) == find_last([1,2,3,2,5]) == find_last("aβcβd")
    @test f_diff_weight(5, 3) == -10
    @test f_diff_weights(3) == [-1, 3, -3, 1]
    @test f_diff_weights_array(3) ==  [[1], [-1, 1], [1, -2, 1], [-1, 3, -3, 1]]
    @test f_diff_expansion_coeffs_lagrange(3,-1) == [1, -1, 0, 0]
    @test f_diff_expansion_coeffs_lagrange(5,2) == [1, 2, 3, 4, 5, 6]
    @test f_diff_expansion_weights(UnitRange(0,5), f_diff_weights_array(5)) == [-5, 29, -69, 85, -55, 15]
    @test [summation_range(7,i,2,1) for i=0:6] == UnitRange{Int64}[1:3, 2:4, 3:5, 4:6, 5:7, 5:7, 5:7]
    @test f_diff_function_sequences([0,1,2,3,4,5,6],2) == [[0, 1, 2], [1, 2, 3], [2, 3, 4], [3, 4, 5], [4, 5, 6], [4, 5, 6], [4, 5, 6]]
    @test lagrange_interpolation([0.0,1,2,3,4,5,6], 0.0..1.0; k=2, m=2) == (0.0:0.08333333333333333:1.0, [0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0])
    @test lagrange_extrapolation([0.0,1,2,3,4,5,6,7], 0.0..1.0; k=2, e=1, m=2) == (1.0:0.07142857142857142:1.1428571428571428, [7.0, 7.5, 8.0])
    @test f_diff_expansion_coeffs_differentiation(2,0) == [0.0, 1.0, 0.5]
    @test lagrange_differentiation([0.0,1,2,3,4,5], 0.0..5.0; k=2, m=1) == (0.0:1.0:5.0, [1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
    @test create_lagrange_differentiation_weights(8,0) == [-761//280, 8//1, -14//1, 56//3, -35//2, 56//5, -14//3, 8//7, -1//8]
    @test create_lagrange_differentiation_matrix(3) == [-11//6 3//1 -3//2 1//3; -1//3 -1//2 1//1 -1//6; 1//6 -1//1 1//2 1//3; -1//3 3//2 -3//1 11//6]
    @test trapezoidal_weights(5; rationalize=true) == [95//288, 317//240, 23//30, 793//720, 157//160]
    @test trapezoidal_integration([1.0, 4.0, 15.0, 40.0, 85.0, 156.0], 0.0..5.0, [3//8, 7//6, 23//24]) ≈ 215.4166666
    @test create_adams_moulton_weights(3;rationalize=true) == [1//24, -5//24, 19//24, 3//8]
    @test f_diff_expansion_coeffs_adams_moulton(5) == [1//1, -1//2, -1//12, -1//24, -19//720, -3//160]
    @test f_diff_expansion_coeffs_adams_bashford(5) == [1//1, 1//2, 5//12, 3//8, 251//720, 95//288]
    @test canonical_partitions(6; header=true) == [[1, 1, 1, 1, 1, 1], [2, 2, 2], [3, 3], [4, 2], [5, 1], [6]]
    @test canonical_partitions(6) == [[1, 1, 1, 1, 1, 1], [2, 2, 2], [3, 3], [4, 2], [5, 1], [6]]
    @test canonical_partitions(6; header=true, reverse=false) == [[6], [5, 1], [4, 2], [3, 3], [2, 2, 2], [1, 1, 1, 1, 1, 1]]
    @test integer_partitions(7) == [[1, 1, 1, 1, 1, 1, 1], [2, 2, 2, 1], [2, 2, 1, 1, 1], [2, 1, 1, 1, 1, 1], [3, 3, 1], [3, 2, 2], [3, 2, 1, 1], [3, 1, 1, 1, 1], [4, 3], [4, 2, 1], [4, 1, 1, 1], [5, 2], [5, 1, 1], [6, 1], [7]]
    @test integer_partitions(9;count=true) == 30
    @test integer_partitions(9) == unique(integer_partitions(9))
    @test integer_partitions(7,4) == [[4, 3], [4, 2, 1], [4, 1, 1, 1]]
    @test integer_partitions(7,4;count=true) == 3
    @test integer_partitions(7; transpose=true) == [[7], [4, 3], [5, 2], [6, 1], [3, 2, 2], [3, 3, 1], [4, 2, 1], [5, 1, 1], [2, 2, 2, 1], [3, 2, 1, 1], [4, 1, 1, 1], [2, 2, 1, 1, 1], [3, 1, 1, 1, 1], [2, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1]]
    @test log10_mantissa.([3,30,300]) == [0.47712125471966244, 0.4771212547196624, 0.4771212547196626]
    @test log10_characteristic_power.([3,30,300]) == [0, 1, 2]
    @test bernoulli_numbers(10) == [1//1, -1//2, 1//6, 0//1, -1//30, 0//1, 1//42, 0//1, -1//30, 0//1, 5//66]
    @test faulhaber_polynom(6) == [0//1, 0//1, -1//12, 0//1, 5//12, 1//2, 1//6]
    @test faulhaber_summation(3,5) == 276
    @test harmonic_number(12) == 86021//27720
    @test harmonic_number(12, 3) == 25535765062457//21300003648000
    @test pascal_triangle(5) == [[1], [1, 1], [1, 2, 1], [1, 3, 3, 1], [1, 4, 6, 4, 1], [1, 5, 10, 10, 5, 1]]
    @test pascal_next([1, 4, 6, 4, 1]) == [1, 5, 10, 10, 5, 1]
    @test polynom(ones(6),2.0) == 63.0
    @test polynom_derivative([1,1,1,1,1]) == [1,2,3,4]
    @test polynom_derivatives([1,1,1,1,1]) == [[1, 2, 3, 4], [2, 6, 12], [6, 24], [24], [0]]
    @test polynom_power([1,1,1],2) == [1, 2, 3, 2, 1]
    @test polynom_powers([1,1,1],3) == [[1, 1, 1], [1, 2, 3, 2, 1], [1, 3, 6, 7, 6, 3, 1]]
    @test polynom_primitive([1,1,1,1,1]) == [0//1, 1//1, 1//2, 1//3, 1//4, 1//5]
    @test polynom_product([1,1],[1,- 1]) == [1, 0, -1]
    @test polynom_product([1,1],[1,- 1,2]) == [1, 0, 1, 2]
    @test polynom_product([1,- 1,2],[1,1]) == [1, 0, 1, 2]
    @test polynom_product_expansion([1,-1,1],[1,1,-1,1,1,1], 5) == [1, 0, -1, 3, -1, 1]
    @test permutations_unique_count([[1,2,3],[2,3,1,4,3]],2) == 60
    @test normalize_VectorRational([2//3,4//5]).num == [10,12]
    @test edges(1:5,2.5,2.5) == [-1.25, 1.25, 3.75, 6.25, 8.75]
    @test steps([4,2,6]) == [0, 4, 6, 12]
    @test stepcenters([4,2,6]) == [2.0, 5.0, 9.0]
    @test stepedges([4,2,6]) == [0, 4, 6, 12]
    @test select125([1,2,4,6,8,10,13,16,18,20,40,60,80,100]) == [2, 6, 10, 16, 20, 60, 100]
    @test step125.([5,10,21.3,50,100.1]) == [1, 2, 5, 10, 20]
    #@test fits_info("T01.fits") == "T01.fits: file was found (for more information set info=true)"
    #@test fits_copy("T01.fits") == "T01.fits was saved as T01 - Copy.fits"
    #@test fits_copy("T01.fits","T01a.fits";protect=false) == "T01.fits was saved as T01a.fits"
    @test fits_create()
    @test fits_extend()
    @test fits_read()
    @test fits_add_key()
    @test fits_edit_key()
    @test fits_delete_key()
    @test fits_rename_key()
end


text:
# Tables.jl interface

Tables.istable(::Type{<:TableRowIterator}) = true
Tables.rowaccess(::Type{<:TableRowIterator}) = true
Tables.rows(itr::TableRowIterator) = itr
Tables.schema(itr::TableRowIterator) = Tables.Schema(itr.index.column_labels, fill(Any, length(itr.index.column_labels)))
Tables.columnnames(tr::TableRow) = tr.index.column_labels
Tables.getcolumn(tr::TableRow, nm::Symbol) = getdata(tr, nm)
Tables.getcolumn(tr::TableRow, i::Integer) = getdata(tr, i)

function writetable(filename::AbstractString, x; kw...)

	_as_vector(y::AbstractVector) = y
	_as_vector(y) = collect(y)

    writetable(filename, Any[_as_vector(c) for c in Tables.Columns(x)], collect(Symbol, Tables.columnnames(x)); kw...)
end


text:
#==============================================================================#
# Cloud9.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/cloud9-2017-09-23.normal.json
#==============================================================================#

__precompile__()

module Cloud9

using AWSCore

"""
    using AWSSDK.Cloud9.create_environment_ec2
    create_environment_ec2([::AWSConfig], arguments::Dict)
    create_environment_ec2([::AWSConfig]; name=, instanceType=, <keyword arguments>)

    using AWSCore.Services.cloud9
    cloud9([::AWSConfig], "CreateEnvironmentEC2", arguments::Dict)
    cloud9([::AWSConfig], "CreateEnvironmentEC2", name=, instanceType=, <keyword arguments>)

# CreateEnvironmentEC2 Operation

Creates an AWS Cloud9 development environment, launches an Amazon Elastic Compute Cloud (Amazon EC2) instance, and then connects from the instance to the environment.

# Arguments

## `name = ::String` -- *Required*
The name of the environment to create.

This name is visible to other AWS IAM users in the same AWS account.

## `description = ::String`
The description of the environment to create.

## `clientRequestToken = ::String`
A unique, case-sensitive string that helps AWS Cloud9 to ensure this operation completes no more than one time.

For more information, see [Client Tokens](http://docs.aws.amazon.com/AWSEC2/latest/APIReference/Run_Instance_Idempotency.html) in the *Amazon EC2 API Reference*.

## `instanceType = ::String` -- *Required*
The type of instance to connect to the environment (for example, `t2.micro`).

## `subnetId = ::String`
The ID of the subnet in Amazon VPC that AWS Cloud9 will use to communicate with the Amazon EC2 instance.

## `automaticStopTimeMinutes = ::Int`
The number of minutes until the running instance is shut down after the environment has last been used.

## `ownerArn = ::String`
The Amazon Resource Name (ARN) of the environment owner. This ARN can be the ARN of any AWS IAM principal. If this value is not specified, the ARN defaults to this environment's creator.

# Returns

`CreateEnvironmentEC2Result`

# Exceptions

`BadRequestException`, `ConflictException`, `NotFoundException`, `ForbiddenException`, `TooManyRequestsException`, `LimitExceededException` or `InternalServerErrorException`.

# Example: CreateEnvironmentEC2

Input:
```
[
    "name" => "my-demo-environment",
    "automaticStopTimeMinutes" => 60,
    "description" => "This is my demonstration environment.",
    "instanceType" => "t2.micro",
    "ownerArn" => "arn:aws:iam::123456789012:user/MyDemoUser",
    "subnetId" => "subnet-1fab8aEX"
]
```

Output:
```
Dict(
    "environmentId" => "8d9967e2f0624182b74e7690ad69ebEX"
)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/cloud9-2017-09-23/CreateEnvironmentEC2)
"""
@inline create_environment_ec2(aws::AWSConfig=default_aws_config(); args...) = create_environment_ec2(aws, args)

@inline create_environment_ec2(aws::AWSConfig, args) = AWSCore.Services.cloud9(aws, "CreateEnvironmentEC2", args)

@inline create_environment_ec2(args) = create_environment_ec2(default_aws_config(), args)

"""
    using AWSSDK.Cloud9.create_environment_membership
    create_environment_membership([::AWSConfig], arguments::Dict)
    create_environment_membership([::AWSConfig]; environmentId=, userArn=, permissions=)

    using AWSCore.Services.cloud9
    cloud9([::AWSConfig], "CreateEnvironmentMembership", arguments::Dict)
    cloud9([::AWSConfig], "CreateEnvironmentMembership", environmentId=, userArn=, permissions=)

# CreateEnvironmentMembership Operation

Adds an environment member to an AWS Cloud9 development environment.

# Arguments

## `environmentId = ::String` -- *Required*
The ID of the environment that contains the environment member you want to add.

## `userArn = ::String` -- *Required*
The Amazon Resource Name (ARN) of the environment member you want to add.

## `permissions = "read-write" or "read-only"` -- *Required*
The type of environment member permissions you want to associate with this environment member. Available values include:

*   `read-only`: Has read-only access to the environment.

*   `read-write`: Has read-write access to the environment.

# Returns

`CreateEnvironmentMembershipResult`

# Exceptions

`BadRequestException`, `ConflictException`, `NotFoundException`, `ForbiddenException`, `TooManyRequestsException`, `LimitExceededException` or `InternalServerErrorException`.

# Example: CreateEnvironmentMembership

Input:
```
[
    "environmentId" => "8d9967e2f0624182b74e7690ad69ebEX",
    "permissions" => "read-write",
    "userArn" => "arn:aws:iam::123456789012:user/AnotherDemoUser"
]
```

Output:
```
Dict(
    "membership" => Dict(
        "environmentId" => "8d9967e2f0624182b74e7690ad69ebEX",
        "permissions" => "read-write",
        "userArn" => "arn:aws:iam::123456789012:user/AnotherDemoUser",
        "userId" => "AIDAJ3BA6O2FMJWCWXHEX"
    )
)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/cloud9-2017-09-23/CreateEnvironmentMembership)
"""
@inline create_environment_membership(aws::AWSConfig=default_aws_config(); args...) = create_environment_membership(aws, args)

@inline create_environment_membership(aws::AWSConfig, args) = AWSCore.Services.cloud9(aws, "CreateEnvironmentMembership", args)

@inline create_environment_membership(args) = create_environment_membership(default_aws_config(), args)

"""
    using AWSSDK.Cloud9.delete_environment
    delete_environment([::AWSConfig], arguments::Dict)
    delete_environment([::AWSConfig]; environmentId=)

    using AWSCore.Services.cloud9
    cloud9([::AWSConfig], "DeleteEnvironment", arguments::Dict)
    cloud9([::AWSConfig], "DeleteEnvironment", environmentId=)

# DeleteEnvironment Operation

Deletes an AWS Cloud9 development environment. If an Amazon EC2 instance is connected to the environment, also terminates the instance.

# Arguments

## `environmentId = ::String` -- *Required*
The ID of the environment to delete.

# Returns

`DeleteEnvironmentResult`

# Exceptions

`BadRequestException`, `ConflictException`, `NotFoundException`, `ForbiddenException`, `TooManyRequestsException`, `LimitExceededException` or `InternalServerErrorException`.

# Example: DeleteEnvironment

Input:
```
[
    "environmentId" => "8d9967e2f0624182b74e7690ad69ebEX"
]
```

Output:
```
Dict(

)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/cloud9-2017-09-23/DeleteEnvironment)
"""
@inline delete_environment(aws::AWSConfig=default_aws_config(); args...) = delete_environment(aws, args)

@inline delete_environment(aws::AWSConfig, args) = AWSCore.Services.cloud9(aws, "DeleteEnvironment", args)

@inline delete_environment(args) = delete_environment(default_aws_config(), args)

"""
    using AWSSDK.Cloud9.delete_environment_membership
    delete_environment_membership([::AWSConfig], arguments::Dict)
    delete_environment_membership([::AWSConfig]; environmentId=, userArn=)

    using AWSCore.Services.cloud9
    cloud9([::AWSConfig], "DeleteEnvironmentMembership", arguments::Dict)
    cloud9([::AWSConfig], "DeleteEnvironmentMembership", environmentId=, userArn=)

# DeleteEnvironmentMembership Operation

Deletes an environment member from an AWS Cloud9 development environment.

# Arguments

## `environmentId = ::String` -- *Required*
The ID of the environment to delete the environment member from.

## `userArn = ::String` -- *Required*
The Amazon Resource Name (ARN) of the environment member to delete from the environment.

# Returns

`DeleteEnvironmentMembershipResult`

# Exceptions

`BadRequestException`, `ConflictException`, `NotFoundException`, `ForbiddenException`, `TooManyRequestsException`, `LimitExceededException` or `InternalServerErrorException`.

# Example: DeleteEnvironmentMembership

Input:
```
[
    "environmentId" => "8d9967e2f0624182b74e7690ad69ebEX",
    "userArn" => "arn:aws:iam::123456789012:user/AnotherDemoUser"
]
```

Output:
```
Dict(

)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/cloud9-2017-09-23/DeleteEnvironmentMembership)
"""
@inline delete_environment_membership(aws::AWSConfig=default_aws_config(); args...) = delete_environment_membership(aws, args)

@inline delete_environment_membership(aws::AWSConfig, args) = AWSCore.Services.cloud9(aws, "DeleteEnvironmentMembership", args)

@inline delete_environment_membership(args) = delete_environment_membership(default_aws_config(), args)

"""
    using AWSSDK.Cloud9.describe_environment_memberships
    describe_environment_memberships([::AWSConfig], arguments::Dict)
    describe_environment_memberships([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.cloud9
    cloud9([::AWSConfig], "DescribeEnvironmentMemberships", arguments::Dict)
    cloud9([::AWSConfig], "DescribeEnvironmentMemberships", <keyword arguments>)

# DescribeEnvironmentMemberships Operation

Gets information about environment members for an AWS Cloud9 development environment.

# Arguments

## `userArn = ::String`
The Amazon Resource Name (ARN) of an individual environment member to get information about. If no value is specified, information about all environment members are returned.

## `environmentId = ::String`
The ID of the environment to get environment member information about.

## `permissions = ["owner", "read-write" or "read-only", ...]`
The type of environment member permissions to get information about. Available values include:

*   `owner`: Owns the environment.

*   `read-only`: Has read-only access to the environment.

*   `read-write`: Has read-write access to the environment.

If no value is specified, information about all environment members are returned.

## `nextToken = ::String`
During a previous call, if there are more than 25 items in the list, only the first 25 items are returned, along with a unique string called a *next token*. To get the next batch of items in the list, call this operation again, adding the next token to the call. To get all of the items in the list, keep calling this operation with each subsequent next token that is returned, until no more next tokens are returned.

## `maxResults = ::Int`
The maximum number of environment members to get information about.

# Returns

`DescribeEnvironmentMembershipsResult`

# Exceptions

`BadRequestException`, `ConflictException`, `NotFoundException`, `ForbiddenException`, `TooManyRequestsException`, `LimitExceededException` or `InternalServerErrorException`.

# Example: DescribeEnvironmentMemberships1

The following example gets information about all of the environment members for the specified AWS Cloud9 development environment.

Input:
```
[
    "environmentId" => "8d9967e2f0624182b74e7690ad69ebEX"
]
```

Output:
```
Dict(
    "memberships" => [
        Dict(
            "environmentId" => "8d9967e2f0624182b74e7690ad69ebEX",
            "permissions" => "read-write",
            "userArn" => "arn:aws:iam::123456789012:user/AnotherDemoUser",
            "userId" => "AIDAJ3BA6O2FMJWCWXHEX"
        ),
        Dict(
            "environmentId" => "8d9967e2f0624182b74e7690ad69ebEX",
            "permissions" => "owner",
            "userArn" => "arn:aws:iam::123456789012:user/MyDemoUser",
            "userId" => "AIDAJNUEDQAQWFELJDLEX"
        )
    ]
)
```

# Example: DescribeEnvironmentMemberships2

The following example gets information about the owner of the specified AWS Cloud9 development environment.

Input:
```
[
    "environmentId" => "8d9967e2f0624182b74e7690ad69ebEX",
    "permissions" => [
        "owner"
    ]
]
```

Output:
```
Dict(
    "memberships" => [
        Dict(
            "environmentId" => "8d9967e2f0624182b74e7690ad69ebEX",
            "permissions" => "owner",
            "userArn" => "arn:aws:iam::123456789012:user/MyDemoUser",
            "userId" => "AIDAJNUEDQAQWFELJDLEX"
        )
    ]
)
```

# Example: DescribeEnvironmentMemberships3

The following example gets AWS Cloud9 development environment membership information for the specified user.

Input:
```
[
    "userArn" => "arn:aws:iam::123456789012:user/MyDemoUser"
]
```

Output:
```
Dict(
    "memberships" => [
        Dict(
            "environmentId" => "10a75714bd494714929e7f5ec4125aEX",
            "lastAccess" => "2018-01-19T11:06:13Z",
            "permissions" => "owner",
            "userArn" => "arn:aws:iam::123456789012:user/MyDemoUser",
            "userId" => "AIDAJNUEDQAQWFELJDLEX"
        ),
        Dict(
            "environmentId" => "12bfc3cd537f41cb9776f8af5525c9EX",
            "lastAccess" => "2018-01-19T11:39:19Z",
            "permissions" => "owner",
            "userArn" => "arn:aws:iam::123456789012:user/MyDemoUser",
            "userId" => "AIDAJNUEDQAQWFELJDLEX"
        )
    ]
)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/cloud9-2017-09-23/DescribeEnvironmentMemberships)
"""
@inline describe_environment_memberships(aws::AWSConfig=default_aws_config(); args...) = describe_environment_memberships(aws, args)

@inline describe_environment_memberships(aws::AWSConfig, args) = AWSCore.Services.cloud9(aws, "DescribeEnvironmentMemberships", args)

@inline describe_environment_memberships(args) = describe_environment_memberships(default_aws_config(), args)

"""
    using AWSSDK.Cloud9.describe_environment_status
    describe_environment_status([::AWSConfig], arguments::Dict)
    describe_environment_status([::AWSConfig]; environmentId=)

    using AWSCore.Services.cloud9
    cloud9([::AWSConfig], "DescribeEnvironmentStatus", arguments::Dict)
    cloud9([::AWSConfig], "DescribeEnvironmentStatus", environmentId=)

# DescribeEnvironmentStatus Operation

Gets status information for an AWS Cloud9 development environment.

# Arguments

## `environmentId = ::String` -- *Required*
The ID of the environment to get status information about.

# Returns

`DescribeEnvironmentStatusResult`

# Exceptions

`BadRequestException`, `ConflictException`, `NotFoundException`, `ForbiddenException`, `TooManyRequestsException`, `LimitExceededException` or `InternalServerErrorException`.

# Example: DescribeEnvironmentStatus

Input:
```
[
    "environmentId" => "8d9967e2f0624182b74e7690ad69ebEX"
]
```

Output:
```
Dict(
    "message" => "Environment is ready to use",
    "status" => "ready"
)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/cloud9-2017-09-23/DescribeEnvironmentStatus)
"""
@inline describe_environment_status(aws::AWSConfig=default_aws_config(); args...) = describe_environment_status(aws, args)

@inline describe_environment_status(aws::AWSConfig, args) = AWSCore.Services.cloud9(aws, "DescribeEnvironmentStatus", args)

@inline describe_environment_status(args) = describe_environment_status(default_aws_config(), args)

"""
    using AWSSDK.Cloud9.describe_environments
    describe_environments([::AWSConfig], arguments::Dict)
    describe_environments([::AWSConfig]; environmentIds=)

    using AWSCore.Services.cloud9
    cloud9([::AWSConfig], "DescribeEnvironments", arguments::Dict)
    cloud9([::AWSConfig], "DescribeEnvironments", environmentIds=)

# DescribeEnvironments Operation

Gets information about AWS Cloud9 development environments.

# Arguments

## `environmentIds = [::String, ...]` -- *Required*
The IDs of individual environments to get information about.

# Returns

`DescribeEnvironmentsResult`

# Exceptions

`BadRequestException`, `ConflictException`, `NotFoundException`, `ForbiddenException`, `TooManyRequestsException`, `LimitExceededException` or `InternalServerErrorException`.

# Example: DescribeEnvironments

Input:
```
[
    "environmentIds" => [
        "8d9967e2f0624182b74e7690ad69ebEX",
        "349c86d4579e4e7298d500ff57a6b2EX"
    ]
]
```

Output:
```
Dict(
    "environments" => [
        Dict(
            "name" => "my-demo-environment",
            "type" => "ec2",
            "arn" => "arn:aws:cloud9:us-east-2:123456789012:environment:8d9967e2f0624182b74e7690ad69ebEX",
            "description" => "This is my demonstration environment.",
            "id" => "8d9967e2f0624182b74e7690ad69ebEX",
            "ownerArn" => "arn:aws:iam::123456789012:user/MyDemoUser"
        ),
        Dict(
            "name" => "another-demo-environment",
            "type" => "ssh",
            "arn" => "arn:aws:cloud9:us-east-2:123456789012:environment:349c86d4579e4e7298d500ff57a6b2EX",
            "id" => "349c86d4579e4e7298d500ff57a6b2EX",
            "ownerArn" => "arn:aws:sts::123456789012:assumed-role/AnotherDemoUser/AnotherDemoUser"
        )
    ]
)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/cloud9-2017-09-23/DescribeEnvironments)
"""
@inline describe_environments(aws::AWSConfig=default_aws_config(); args...) = describe_environments(aws, args)

@inline describe_environments(aws::AWSConfig, args) = AWSCore.Services.cloud9(aws, "DescribeEnvironments", args)

@inline describe_environments(args) = describe_environments(default_aws_config(), args)

"""
    using AWSSDK.Cloud9.list_environments
    list_environments([::AWSConfig], arguments::Dict)
    list_environments([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.cloud9
    cloud9([::AWSConfig], "ListEnvironments", arguments::Dict)
    cloud9([::AWSConfig], "ListEnvironments", <keyword arguments>)

# ListEnvironments Operation

Gets a list of AWS Cloud9 development environment identifiers.

# Arguments

## `nextToken = ::String`
During a previous call, if there are more than 25 items in the list, only the first 25 items are returned, along with a unique string called a *next token*. To get the next batch of items in the list, call this operation again, adding the next token to the call. To get all of the items in the list, keep calling this operation with each subsequent next token that is returned, until no more next tokens are returned.

## `maxResults = ::Int`
The maximum number of environments to get identifiers for.

# Returns

`ListEnvironmentsResult`

# Exceptions

`BadRequestException`, `ConflictException`, `NotFoundException`, `ForbiddenException`, `TooManyRequestsException`, `LimitExceededException` or `InternalServerErrorException`.

# Example: ListEnvironments

Input:
```
[

]
```

Output:
```
Dict(
    "environmentIds" => [
        "349c86d4579e4e7298d500ff57a6b2EX",
        "45a3da47af0840f2b0c0824f5ee232EX"
    ]
)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/cloud9-2017-09-23/ListEnvironments)
"""
@inline list_environments(aws::AWSConfig=default_aws_config(); args...) = list_environments(aws, args)

@inline list_environments(aws::AWSConfig, args) = AWSCore.Services.cloud9(aws, "ListEnvironments", args)

@inline list_environments(args) = list_environments(default_aws_config(), args)

"""
    using AWSSDK.Cloud9.update_environment
    update_environment([::AWSConfig], arguments::Dict)
    update_environment([::AWSConfig]; environmentId=, <keyword arguments>)

    using AWSCore.Services.cloud9
    cloud9([::AWSConfig], "UpdateEnvironment", arguments::Dict)
    cloud9([::AWSConfig], "UpdateEnvironment", environmentId=, <keyword arguments>)

# UpdateEnvironment Operation

Changes the settings of an existing AWS Cloud9 development environment.

# Arguments

## `environmentId = ::String` -- *Required*
The ID of the environment to change settings.

## `name = ::String`
A replacement name for the environment.

## `description = ::String`
Any new or replacement description for the environment.

# Returns

`UpdateEnvironmentResult`

# Exceptions

`BadRequestException`, `ConflictException`, `NotFoundException`, `ForbiddenException`, `TooManyRequestsException`, `LimitExceededException` or `InternalServerErrorException`.

# Example: UpdateEnvironment

Input:
```
[
    "name" => "my-changed-demo-environment",
    "description" => "This is my changed demonstration environment.",
    "environmentId" => "8d9967e2f0624182b74e7690ad69ebEX"
]
```

Output:
```
Dict(

)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/cloud9-2017-09-23/UpdateEnvironment)
"""
@inline update_environment(aws::AWSConfig=default_aws_config(); args...) = update_environment(aws, args)

@inline update_environment(aws::AWSConfig, args) = AWSCore.Services.cloud9(aws, "UpdateEnvironment", args)

@inline update_environment(args) = update_environment(default_aws_config(), args)

"""
    using AWSSDK.Cloud9.update_environment_membership
    update_environment_membership([::AWSConfig], arguments::Dict)
    update_environment_membership([::AWSConfig]; environmentId=, userArn=, permissions=)

    using AWSCore.Services.cloud9
    cloud9([::AWSConfig], "UpdateEnvironmentMembership", arguments::Dict)
    cloud9([::AWSConfig], "UpdateEnvironmentMembership", environmentId=, userArn=, permissions=)

# UpdateEnvironmentMembership Operation

Changes the settings of an existing environment member for an AWS Cloud9 development environment.

# Arguments

## `environmentId = ::String` -- *Required*
The ID of the environment for the environment member whose settings you want to change.

## `userArn = ::String` -- *Required*
The Amazon Resource Name (ARN) of the environment member whose settings you want to change.

## `permissions = "read-write" or "read-only"` -- *Required*
The replacement type of environment member permissions you want to associate with this environment member. Available values include:

*   `read-only`: Has read-only access to the environment.

*   `read-write`: Has read-write access to the environment.

# Returns

`UpdateEnvironmentMembershipResult`

# Exceptions

`BadRequestException`, `ConflictException`, `NotFoundException`, `ForbiddenException`, `TooManyRequestsException`, `LimitExceededException` or `InternalServerErrorException`.

# Example: UpdateEnvironmentMembership

Input:
```
[
    "environmentId" => "8d9967e2f0624182b74e7690ad69ebEX",
    "permissions" => "read-only",
    "userArn" => "arn:aws:iam::123456789012:user/AnotherDemoUser"
]
```

Output:
```
Dict(
    "membership" => Dict(
        "environmentId" => "8d9967e2f0624182b74e7690ad69eb31",
        "permissions" => "read-only",
        "userArn" => "arn:aws:iam::123456789012:user/AnotherDemoUser",
        "userId" => "AIDAJ3BA6O2FMJWCWXHEX"
    )
)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/cloud9-2017-09-23/UpdateEnvironmentMembership)
"""
@inline update_environment_membership(aws::AWSConfig=default_aws_config(); args...) = update_environment_membership(aws, args)

@inline update_environment_membership(aws::AWSConfig, args) = AWSCore.Services.cloud9(aws, "UpdateEnvironmentMembership", args)

@inline update_environment_membership(args) = update_environment_membership(default_aws_config(), args)

end # module Cloud9

#==============================================================================#
# End of file
#==============================================================================#


text:
# index.jl - tim.sterne.weiler@utoronto.ca, 1/28/16

abstract type SeqLibrary end

struct GraphLib <: SeqLibrary
   offset::Vector{CoordInt}
   names::Vector{GeneName}
   info::Vector{GeneInfo}
   lengths::Vector{Float64}
   graphs::Vector{SpliceGraph}
   edges::Edges
   index::FMIndex
   sorted::Bool
   kmer::Int64
end

function build_chrom_dict( ref::RefSet )
   ret = Dict{SeqName,Vector{GeneName}}() # refgenomeseq->geneid[]
   for g in keys(ref)
      chrom = ref[g].info.name
      if !haskey(ret, chrom)
         ret[chrom] = Vector{GeneName}()
      end
      push!(ret[chrom], g)
   end
   ret
end

function twobit_enc(seq)
   len = length(seq)
   ret = IntVector{2,UInt8}(len)
   for i in 1:len
      val = isambiguous(seq[i]) ? 0x00 : convert(UInt8, trailing_zeros(seq[i]))
      ret[i] = val
   end
   ret
end

function load_fasta( fhIter; verbose=false )
   seq = SGSequence(mutable=false)
   offset = Int[]
   names = SeqName[]
   for r in fhIter
      immutable!(r.seq)
      push!(names, r.name)
      push!(offset, length(seq)+1)
      println( STDERR, r )
      @time seq *= r.seq
   end
   seq, offset, names
end

function single_genome_index!( fhIter; verbose=false )
   seq,offset,names = load_fasta( fhIter )
   @time fm = FMIndex(twobit_enc(seq), 4, r=4, program=:SuffixArrays, mmap=true)
   println( STDERR, "Finished building Index..." )

   SeqLib(seq, offset, names, fm, false)
end

function trans_index!( fhIter, ref::RefSet; kmer=9 )
   seqdic  = build_chrom_dict( ref )
   xcript  = dna""
   xoffset = Vector{UInt64}()
   xgenes  = Vector{GeneName}()
   xinfo   = Vector{GeneInfo}()
   xlength = Vector{Float64}()
   xgraph  = Vector{SpliceGraph}()
   # set up splice graphs
   runoffset = 0
   num = 0
   for r in fhIter
      name = String(r.data[r.identifier])
      haskey(seqdic, name) || continue
      #BioSequences.immutable!(r.seq)
      sg = SGSequence( r.data[r.sequence] )

      println(STDERR, "Building Splice Graphs for $( name ).." )
      @timer for g in seqdic[name]
         #println( STDERR, "Building $g Splice Graph..." )
         curgraph = SpliceGraph( ref[g], sg, kmer )
         xcript  *= curgraph.seq
         push!(xgraph, curgraph)
         push!(xgenes, g)
         push!(xinfo, ref[g].info )
         push!(xlength, ref[g].length )
         push!(xoffset, runoffset)
         runoffset += length(curgraph.seq)
         num += 1
      end
   end
   if num > 1
      println( STDERR, "Building full sg-index from $num genes..." )
      @time fm = FMIndex(twobit_enc(xcript), 4, r=1, program=:SuffixArrays, mmap=true)
   else
      error("ERROR: No genes in the GTF file matched chromosome names in the FASTA file!")
   end

   # clean up
   xcript = DNASequence()
   gc()

   println( STDERR, "Building edges.." )
   @time edges = build_edges( xgraph, kmer ) # TODO make variable kmer

   GraphLib( xoffset, xgenes, xinfo, xlength, xgraph, edges, fm, true, kmer )
end

function fasta_to_index( filename::String, ref::RefSet; kmer=9 )
   if isgzipped( filename )
      println(STDERR, "Decompressing and Indexing $filename...")
      to_open = open( filename ) |> x->ZlibInflateInputStream(x, reset_on_end=true)
   else
      println(STDERR, "Indexing $filename...")
      to_open = open( filename ) |> BufferedInputStream
   end
   # iterate through fasta entries
   index = @time trans_index!(FASTA.Reader( to_open ), ref, kmer=kmer)
   index
end


text:
using Documenter
makedocs(sitename="TaxAnalytics.jl",
	  pages = [
		   "Introduction" => "index.md"
		   "About Taxes" => "tax.md"
		   "ELT vs ETL" => "elt.md"
		   "Self-Service Data Analysis" => "service.md"
		   "Data Visualization" => "visuals.md"
		   "Data Science Operations" => "ds.md"
		   "p-value Null hypothesis tests" => "thesis.md"
		   "Feature Engineering" => "features.md"
		   "Deploy ML Models" => "ml.md"
		   ])


text:
@testset "DataHR tests" begin show_testset_section() #Scope for test data

using MDDatasets
using Test

#==Input data
===============================================================================#
#Sample DataF1 leaf dataset:
d1 = DataF1([1,2,3], [4,5,6])
sweeplist1 = PSweep[
	PSweep("v1", [1,2])
	PSweep("v2", [4,5])
]
sweeplist2 = PSweep[ #2nd sweep is different variable than sweeplist1
	PSweep("v1", [1,2])
	PSweep("v4", [4,5])
]
sweeplist3 = PSweep[ #Points of 2nd sweep are different than sweeplist1
	PSweep("v1", [1,2])
	PSweep("v2", [3,4])
]

#Construct using low-level functions:
dhr1_ll = DataHR(sweeplist1,DataF1[d1 d1; d1 d1])

dhr1 = fill(d1, DataHR, sweeplist1)
dhr2 = fill(d1, DataHR, sweeplist2)
dhr3 = fill(d1, DataHR, sweeplist3)

#==Tests
===============================================================================#
@testset "Validate _datamatch()" begin show_testset_description()
	other = dhr1*1.0
	@test !_datamatch(dhr1, other+1) #Make sure _datamatch works
	@test _datamatch(dhr1, other)
end

@testset "Construction of DataHR objects" begin show_testset_description()
	dhr1_0 = zeros(DataHR, sweeplist1)
	dhr1_1 = ones(DataHR, sweeplist1)
	dhr1_2 = fill(2, DataHR, sweeplist1)

	@test _datamatch(dhr1_2-dhr1_1, dhr1_1)
	@test _datamatch(dhr1_2-dhr1_1, dhr1_1)
	@test _datamatch(dhr1_2-2*dhr1_1, dhr1_0)

	@test _datamatch(dhr1_0.elem, zeros(Float64, size(DataHR, sweeplist1)))
	@test _datamatch(dhr1_1.elem, ones(Float64, size(DataHR, sweeplist1)))
end

@testset "Broadcast operations on DataHR" begin show_testset_description()
	s1 = dhr1 + dhr1
	@testset "Validate s1 = dhr1+dhr2" begin
		for (es1, e1) in zip(s1.elem, dhr1.elem)
			@test es1.x == e1.x
			@test es1.y == (e1.y .+ e1.y)
		end
	end
	@test_throws ArgumentError dhr1+dhr2 #Mismatched sweeps
	@test_throws ArgumentError dhr1+dhr3 #Mismatched sweeps
end

@testset "Conversion (DataRS <=> DataHR)" begin show_testset_description()
	datahr_1 = dhr1
	datars = convert(DataRS, datahr_1)
	datahr_2 = convert(DataHR, datars)
	@test _datamatch(datahr_1, datahr_2)

	#Verify that the sweeps match:
	@test datars.sweep == dhr1.sweeps[1]
	@test datars.elem[1].sweep == datahr_1.sweeps[2]
	@test datars.elem[2].sweep == datahr_1.sweeps[2]

end

end


text:
using StaticArrays, BenchmarkTools, Unrolled

const i1 = @SVector [1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,8,9,8,9,10,10,]
const i2 = @SVector [3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,4,3,4,3,2,2,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,7,6,7,6,5,5,6,7,5,5,7,6,]
const i3 = @SVector [5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,8,9,8,9,10,10,9,8,10,10,8,9,]
x = @SVector rand(10)

tdot_A(x) = (dot(x[i1] .* x[i2], x[i3]))
tdot_B(x) = sum(x[i1] .* x[i2] .* x[i3])
# tdot_C(x) = sum(x[i1[n]] .* x[i2[n]] .* x[i3[n]] for n = 1:60)
function tdot_D(x)
    y = zero(eltype(x))
    for n = 1:length(i1)
        @inbounds a = x[i1[n]] * x[i2[n]]
        @inbounds y = muladd(a, x[i3[n]], y)
    end
    return y
end

@btime tdot_A($x)
@btime tdot_B($x)
# @btime tdot_C($x)
@btime tdot_D($x)

# @code_llvm(tdot_A(x))
@code_llvm(tdot_D(x))

@btime ForwardDiff.gradient($tdot_A, $x)
@btime ForwardDiff.gradient($tdot_D, $x)

const j1 = @SVector [1,1,1,2,2,2,3,4,3,4,4,3,4,5,6,7,5,6,7,7,8,9,8,9,10,10,9,5,5,6,7,6,7,8,9,8,9,10,10,]
const j2 = @SVector [3,4,2,1,3,4,2,2,4,3,1,1,1,2,3,4,1,1,1,1,3,4,2,2,4 ,3 ,2,7,6,7,6,5,5,6,7,5,5,7 ,6 ,]
A12_ = rand([0,1], (10,10))
const A12 = SMatrix{10,10}(A12_...)

ddot_A(x) = (dot(x[j1], x[j2]))
ddot_B(x) = (dot(x, A12 * x) )
@inline @unroll function ddot_C(x, j1,j2)
   y = zero(eltype(x))
   @unroll for n = 1:length(x)
     @inbounds y = muladd(x[j1[n]], x[j2[n]], y)
   end
   return y
end

x = @SVector rand(10)

@btime ddot_A($x)
@btime ddot_B($x)
@btime ddot_C($x, $j1, $j2)

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
#!/usr/bin/env julia

const OUTPUT = haskey(ENV, "OUTPUT")

# maximum power density possible (say 300W for a 10mm x 10mm chip)
MAX_PD = 3.0e6
# required precision in degrees
PRECISION = 0.001
SPEC_HEAT_SI = 1.75e6
# capacitance fitting factor
K_SI = 100
FACTOR_CHIP = 0.5

# chip parameters
t_chip = 0.0005
chip_height = 0.016
chip_width = 0.016
# ambient temperature, assuming no package at all
amb_temp = 80.0

#= Single iteration of the transient solver in the grid model.
 = advances the solution of the discretized difference equations
 = by one time step
 =#

function single_iteration(result::SharedArray{Float64,2}, temp::Matrix{Float64}, power::Matrix{Float64}, cap::Float64, rx::Float64, ry::Float64, rz::Float64, step::Float64)
    row = size(temp,1)
    col = size(temp,2)
    @sync @parallel for r in 1:row
        for c in 1:col
            # Corner 1
            if ((r == 1) & (c == 1))
                delta = (step / cap) * (power[1,1] + (temp[1,2] - temp[1,1]) / rx +
                                        (temp[2,1] - temp[1,1]) / ry +
                                        (amb_temp - temp[1,1]) / rz)
            # Corner 2
            elseif ((r == 1) & (c == col))
                delta =
                    (step / cap) * (power[1,c] + (temp[1,c-1] - temp[1,c]) / rx +
                         (temp[2,c] - temp[1,c]) / ry +
                         (amb_temp - temp[1,c]) / rz)
            # Corner 3
            elseif ((r == row) & (c == col))
                delta = (step / cap) *
                    (power[r,c] +
                     (temp[r,c-1] - temp[r,c]) / rx +
                     (temp[r-1,c] - temp[r,c]) / ry +
                     (amb_temp - temp[r,c]) / rz)
            # Corner 4
            elseif ((r == row) & (c == 1))
                delta =
                    (step / cap) *
                    (power[r,1] + (temp[r,2] - temp[r,1]) / rx +
                     (temp[r-1,1] - temp[r,1]) / ry +
                     (amb_temp - temp[r,1]) / rz)
            # Edge 1
            elseif (r == 1)
                delta = (step / cap) *
                        (power[c] +
                         (temp[1,c+1] + temp[1,c-1] - 2.0*temp[1,c]) / rx +
                         (temp[2,c] - temp[1,c]) / ry +
                         (amb_temp - temp[1,c]) / rz)
            # Edge 2
            elseif (c == col)
                delta = (step / cap) *
                        (power[r,c] +
                         (temp[r+1,c] + temp[r-1,c] -
                          2.0*temp[r,c]) /
                             ry +
                         (temp[r,c-1] - temp[r,c]) / rx +
                         (amb_temp - temp[r,c]) / rz)
            # Edge 3
            elseif (r == row)
                delta = (step / cap) *
                        (power[r,c] +
                         (temp[r,c+1] + temp[r,c-1] -
                          2.0*temp[r,c]) /
                             rx +
                         (temp[r-1,c] - temp[r,c]) / ry +
                         (amb_temp - temp[r,c]) / rz)
            # Edge 4
            elseif (c == 1)
                delta =
                    (step / cap) * (power[r,1] +
                                    (temp[r+1,1] + temp[r-1,1] -
                                     2.0*temp[r,1]) /
                                        ry +
                                    (temp[r,1] - temp[r,1]) / rx +
                                    (amb_temp - temp[r,1]) / rz);
            # Inside the chip
            else
                delta = (step / cap) *
                        (power[r,c] +
                         (temp[r+1,c] + temp[r-1,c] -
                          2.0*temp[r,c]) /
                             ry +
                         (temp[r,c+1] + temp[r,c-1] -
                          2.0*temp[r,c]) /
                             rx +
                         (amb_temp - temp[r,c]) / rz);
            end
            # Update Temperatures
            result[r,c] = temp[r,c]+delta
        end
    end
    temp[:] = result[:]
end

#= Transient solver driver routine: simply converts the heat
 = transfer differential equations to difference equations
 = and solves the difference equations by iterating
 =#
function compute_tran_temp(result::SharedArray{Float64,2}, num_iterations::Int32, temp::Matrix{Float64},
                       power::Matrix{Float64})
    row = size(temp,1)
    col = size(temp,2)
    grid_height = chip_height / row
    grid_width = chip_width / col

    cap = FACTOR_CHIP * SPEC_HEAT_SI * t_chip * grid_width * grid_height
    rx = grid_width / (2.0 * K_SI * t_chip * grid_height)
    ry = grid_height / (2.0 * K_SI * t_chip * grid_width)
    rz = t_chip / (K_SI * grid_height * grid_width)

    max_slope = MAX_PD / (FACTOR_CHIP * t_chip * SPEC_HEAT_SI)
    step = PRECISION / max_slope

    for i in 1:num_iterations
        single_iteration(result,temp,power,cap,rx,ry,rz,step)
    end
end

function read_input(vect::Matrix{Float64}, file)
    fp = open(file,"r")
    for r in 1:size(vect,1)
        for c in 1:size(vect,2)
            vect[r,c] = parse(Float64,readline(fp))
        end
    end
    close(fp)
end

function usage(args)
    prog = basename(Base.source_path())
# no MT
#    println(STDERR, "Usage: ",prog," <grid_rows> <grid_cols> <sim_time> <no. of threads> <temp_file> <power_file>")
    println(STDERR, "Usage: ",prog," <grid_rows> <grid_cols> <sim_time> <temp_file> <power_file>")
    println(STDERR,
            "\t<grid_rows>  - number of rows in the grid (positive integer)")
    println(
        STDERR,
        "\t<grid_cols>  - number of columns in the grid (positive integer)")
    println(STDERR, "\t<sim_time>   - number of iterations")
#    println(STDERR, "\t<no. of threads>   - number of threads")
    println(STDERR, "\t<temp_file>  - name of the file containing the initial temperature values of each cell")
    println(STDERR, "\t<power_file> - name of the file containing the dissipated power values of each cell")
    exit(1)
end

function main(args)
    # check validity of inputs
    if length(args) != 5
        usage(args)
    end
    grid_rows = get(tryparse(Int64,args[1]),0)
    grid_cols = get(tryparse(Int64,args[2]),0)
    sim_time = get(tryparse(Int32,args[3]),0)
    if (grid_rows <= 0) |
       (grid_cols <= 0) |
       (sim_time <= 0)
        usage(args)
    end

    # allocate memory for the temperature and power arrays
    temp = Matrix{Float64}(grid_rows,grid_cols)
    power = Matrix{Float64}(grid_rows,grid_cols)
    result::SharedArray{Float64,2} = SharedArray(Float64,(grid_rows,grid_cols), init = result -> result[Base.localindexes(result)] = 0)

    # read initial temperatures and input power
    tfile = args[4]
    pfile = args[5]
    read_input(temp,tfile)
    read_input(power,pfile)

    tic()
    println("Start computing the transient temperature")
    compute_tran_temp(result,sim_time,temp,power)
    println("Ending simulation")
    toc()

    # output results
    output = ""
    if OUTPUT
        if ENV["OUTPUT"] != ""
            println("Writing results output.txt")
            f = open("output.txt","w")
            for r in 1:size(temp,1)
                for c in 1:size(temp,2)
                    println(f,@sprintf("%d\t%g",(r-1)*size(temp,2)+c-1,temp[r,c]))
                end
            end
            close(f)
        else
            println("OUTPUT environment variable not set or empty, not writing results")
        end
    end
end


text:
import Base: ==, +, -, *, /, Broadcast
import Adapt
using Base.Cartesian

"""
    Operator{BL,BR,T} <: DataOperator{BL,BR}

Operator type that stores the representation of an operator on the Hilbert spaces
given by `BL` and `BR` (e.g. a Matrix).
"""
mutable struct Operator{BL,BR,T} <: DataOperator{BL,BR}
    basis_l::BL
    basis_r::BR
    data::T
    function Operator{BL,BR,T}(basis_l::BL,basis_r::BR,data::T) where {BL,BR,T}
        (length.((basis_l,basis_r))==size(data)) || throw(DimensionMismatch("Tried to assign data of size $(size(data)) to bases of length $(length(basis_l)) and $(length(basis_r))!"))
        new(basis_l,basis_r,data)
    end
end
Operator{BL,BR}(basis_l::BL,basis_r::BR,data::T) where {BL,BR,T} = Operator{BL,BR,T}(basis_l,basis_r,data)
Operator(basis_l::BL,basis_r::BR,data::T) where {BL,BR,T} = Operator{BL,BR,T}(basis_l,basis_r,data)
Operator(b::Basis,data) = Operator(b,b,data)

Base.zero(op::Operator) = Operator(op.basis_l,op.basis_r,zero(op.data))
Base.eltype(op::Operator) = eltype(op.data)
Base.eltype(::Type{T}) where {BL,BR,D,T<:Operator{BL,BR,D}} = eltype(D)
Base.size(op::Operator) = size(op.data)
Base.size(op::Operator, d::Int) = size(op.data, d)

# Convert data to CuArray with cu(::Operator)
Adapt.adapt_structure(to, x::Operator) = Operator(x.basis_l, x.basis_r, Adapt.adapt(to, x.data))

const DenseOpPureType{BL,BR} = Operator{BL,BR,<:Matrix}
const DenseOpAdjType{BL,BR} = Operator{BL,BR,<:Adjoint{<:Number,<:Matrix}}
const DenseOpType{BL,BR} = Union{DenseOpPureType{BL,BR},DenseOpAdjType{BL,BR}}
const AdjointOperator{BL,BR} = Operator{BL,BR,<:Adjoint}

"""
    DenseOperator(b1[, b2, data])

Dense array implementation of Operator. Converts any given data to a dense `Matrix`.
"""
DenseOperator(basis_l::Basis,basis_r::Basis,data::T) where T = Operator(basis_l,basis_r,Matrix(data))
DenseOperator(basis_l::Basis,basis_r::Basis,data::Matrix) = Operator(basis_l,basis_r,data)
DenseOperator(b::Basis, data) = DenseOperator(b, b, data)
DenseOperator(::Type{T},b1::Basis,b2::Basis) where T = Operator(b1,b2,zeros(T,length(b1),length(b2)))
DenseOperator(::Type{T},b::Basis) where T = Operator(b,b,zeros(T,length(b),length(b)))
DenseOperator(b1::Basis, b2::Basis) = DenseOperator(ComplexF64, b1, b2)
DenseOperator(b::Basis) = DenseOperator(ComplexF64, b)
DenseOperator(op::DataOperator) = DenseOperator(op.basis_l,op.basis_r,Matrix(op.data))

Base.copy(x::Operator) = Operator(x.basis_l, x.basis_r, copy(x.data))

"""
    dense(op::AbstractOperator)

Convert an arbitrary Operator into a [`DenseOperator`](@ref).
"""
dense(x::AbstractOperator) = DenseOperator(x)

==(x::DataOperator{BL,BR}, y::DataOperator{BL,BR}) where {BL,BR} = (samebases(x,y) && x.data==y.data)
==(x::DataOperator, y::DataOperator) = false
Base.isapprox(x::DataOperator{BL,BR}, y::DataOperator{BL,BR}; kwargs...) where {BL,BR} = (samebases(x,y) && isapprox(x.data, y.data; kwargs...))
Base.isapprox(x::DataOperator, y::DataOperator; kwargs...) = false

# Arithmetic operations
+(a::Operator{BL,BR}, b::Operator{BL,BR}) where {BL,BR} = Operator(a.basis_l, a.basis_r, a.data+b.data)
+(a::Operator, b::Operator) = throw(IncompatibleBases())

-(a::Operator) = Operator(a.basis_l, a.basis_r, -a.data)
-(a::Operator{BL,BR}, b::Operator{BL,BR}) where {BL,BR} = Operator(a.basis_l, a.basis_r, a.data-b.data)
-(a::Operator, b::Operator) = throw(IncompatibleBases())

*(a::Operator{BL,BR}, b::Ket{BR}) where {BL,BR} = Ket{BL}(a.basis_l, a.data*b.data)
*(a::DataOperator, b::Ket) = throw(IncompatibleBases())
*(a::Bra{BL}, b::Operator{BL,BR}) where {BL,BR} = Bra{BR}(b.basis_r, transpose(b.data)*a.data)
*(a::Bra, b::DataOperator) = throw(IncompatibleBases())
*(a::Operator{B1,B2}, b::Operator{B2,B3}) where {B1,B2,B3} = Operator(a.basis_l, b.basis_r, a.data*b.data)
*(a::DataOperator, b::DataOperator) = throw(IncompatibleBases())
*(a::Operator, b::Number) = Operator(a.basis_l, a.basis_r, b*a.data)
*(a::Number, b::Operator) = Operator(b.basis_l, b.basis_r, a*b.data)
function *(op1::AbstractOperator{B1,B2}, op2::Operator{B2,B3,T}) where {B1,B2,B3,T}
    result = Operator{B1,B3,T}(op1.basis_l, op2.basis_r, similar(op2.data,length(op1.basis_l),length(op2.basis_r)))
    mul!(result,op1,op2)
    return result
end
function *(op1::Operator{B1,B2,T}, op2::AbstractOperator{B2,B3}) where {B1,B2,B3,T}
    result = Operator{B1,B3,T}(op1.basis_l, op2.basis_r, similar(op1.data,length(op1.basis_l),length(op2.basis_r)))
    mul!(result,op1,op2)
    return result
end
function *(op::AbstractOperator{BL,BR}, psi::Ket{BR,T}) where {BL,BR,T}
    result = Ket{BL,T}(op.basis_l,similar(psi.data,length(op.basis_l)))
    mul!(result,op,psi)
    return result
end
function *(psi::Bra{BL,T}, op::AbstractOperator{BL,BR}) where {BL,BR,T}
    result = Bra{BR,T}(op.basis_r, similar(psi.data,length(op.basis_r)))
    mul!(result,psi,op)
    return result
end

/(a::Operator, b::Number) = Operator(a.basis_l, a.basis_r, a.data ./ b)

dagger(x::Operator) = Operator(x.basis_r,x.basis_l,adjoint(x.data))
transpose(x::Operator) = Operator(x.basis_r,x.basis_l,transpose(x.data))
ishermitian(A::DataOperator) = false
ishermitian(A::DataOperator{B,B}) where B = ishermitian(A.data)

issymmetric(A::DataOperator) = false
issymmetric(A::DataOperator{B,B}) where B = issymmetric(A.data)

Base.collect(A::Operator) = Operator(A.basis_l, A.basis_r, collect(A.data))

tensor(a::Operator, b::Operator) = Operator(tensor(a.basis_l, b.basis_l), tensor(a.basis_r, b.basis_r), kron(b.data, a.data))

conj(a::Operator) = Operator(a.basis_l, a.basis_r, conj(a.data))
conj!(a::Operator) = (conj!(a.data); a)

"""
    tensor(x::Ket, y::Bra)

Outer product ``|x⟩⟨y|`` of the given states.
"""
tensor(a::Ket, b::Bra) = Operator(a.basis, b.basis, reshape(kron(b.data, a.data), length(a.basis), length(b.basis)))

tr(op::Operator{B,B}) where B = tr(op.data)

function ptrace(a::DataOperator, indices)
    check_ptrace_arguments(a, indices)
    rank = length(a.basis_l.shape)
    result = _ptrace(Val{rank}, a.data, a.basis_l.shape, a.basis_r.shape, indices)
    return Operator(ptrace(a.basis_l, indices), ptrace(a.basis_r, indices), result)
end
ptrace(op::AdjointOperator, indices) = dagger(ptrace(op, indices))

function ptrace(psi::Ket, indices)
    check_ptrace_arguments(psi, indices)
    b = basis(psi)
    b_ = ptrace(b, indices)
    rank = length(b.shape)
    result = _ptrace_ket(Val{rank}, psi.data, b.shape, indices)::Matrix{eltype(psi)}
    return Operator(b_, b_, result)
end

function ptrace(psi::Bra, indices)
    check_ptrace_arguments(psi, indices)
    b = basis(psi)
    b_ = ptrace(b, indices)
    rank = length(b.shape)
    result = _ptrace_bra(Val{rank}, psi.data, b.shape, indices)::Matrix{eltype(psi)}
    return Operator(b_, b_, result)
end

normalize!(op::Operator) = (rmul!(op.data, 1.0/tr(op)); op)

function expect(op::DataOperator{B,B}, state::Ket{B}) where B
    dot(state.data, op.data, state.data)
end

function expect(op::DataOperator{B1,B2}, state::DataOperator{B2,B2}) where {B1,B2}
    check_samebases(op, state)
    result = zero(promote_type(eltype(op),eltype(state)))
    @inbounds for i=1:size(op.data, 1), j=1:size(op.data,2)
        result += op.data[i,j]*state.data[j,i]
    end
    result
end

function exp(op::T) where {B,T<:DenseOpType{B,B}}
    return DenseOperator(op.basis_l, op.basis_r, exp(op.data))
end

function permutesystems(a::Operator{B1,B2}, perm) where {B1<:CompositeBasis,B2<:CompositeBasis}
    @assert length(a.basis_l.bases) == length(a.basis_r.bases) == length(perm)
    @assert isperm(perm)
    data = reshape(a.data, [a.basis_l.shape; a.basis_r.shape]...)
    data = permutedims(data, [perm; perm .+ length(perm)])
    data = reshape(data, length(a.basis_l), length(a.basis_r))
    return Operator(permutesystems(a.basis_l, perm), permutesystems(a.basis_r, perm), data)
end
permutesystems(a::AdjointOperator{B1,B2}, perm) where {B1<:CompositeBasis,B2<:CompositeBasis} = dagger(permutesystems(dagger(a),perm))

identityoperator(::Type{S}, ::Type{T}, b1::Basis, b2::Basis) where {S<:DenseOpType,T<:Number} =
    Operator(b1, b2, Matrix{T}(I, length(b1), length(b2)))

"""
    projector(a::Ket, b::Bra)

Projection operator ``|a⟩⟨b|``.
"""
projector(a::Ket, b::Bra) = tensor(a, b)
"""
    projector(a::Ket)

Projection operator ``|a⟩⟨a|``.
"""
projector(a::Ket) = Operator(a.basis, a.data*a.data')
"""
    projector(a::Bra)

Projection operator ``|a⟩⟨a|``.
"""
projector(a::Bra) = projector(a')

"""
    dm(a::StateVector)

Create density matrix ``|a⟩⟨a|``. Same as `projector(a)`.
"""
dm(x::Ket) = tensor(x, dagger(x))
dm(x::Bra) = tensor(dagger(x), x)

# Partial trace implementation for dense operators.
function _strides(shape)
    N = length(shape)
    S = zeros(eltype(shape), N)
    S[1] = 1
    for m=2:N
        S[m] = S[m-1]*shape[m-1]
    end
    return S
end

# Dense operator version
@generated function _ptrace(::Type{Val{RANK}}, a,
                            shape_l, shape_r,
                            indices) where RANK
    return quote
        a_strides_l = _strides(shape_l)
        result_shape_l = copy(shape_l)
        @inbounds for idx ∈ indices
            result_shape_l[idx] = 1
        end
        result_strides_l = _strides(result_shape_l)
        a_strides_r = _strides(shape_r)
        result_shape_r = copy(shape_r)
        @inbounds for idx ∈ indices
            result_shape_r[idx] = 1
        end
        result_strides_r = _strides(result_shape_r)
        N_result_l = prod(result_shape_l)
        N_result_r = prod(result_shape_r)
        result = zeros(eltype(a), N_result_l, N_result_r)
        @nexprs 1 (d->(Jr_{$RANK}=1;Ir_{$RANK}=1))
        @nloops $RANK ir (d->1:shape_r[d]) (d->(Ir_{d-1}=Ir_d; Jr_{d-1}=Jr_d)) (d->(Ir_d+=a_strides_r[d]; if !(d in indices) Jr_d+=result_strides_r[d] end)) begin
            @nexprs 1 (d->(Jl_{$RANK}=1;Il_{$RANK}=1))
            @nloops $RANK il (k->1:shape_l[k]) (k->(Il_{k-1}=Il_k; Jl_{k-1}=Jl_k; if (k in indices && il_k!=ir_k) Il_k+=a_strides_l[k]; continue end)) (k->(Il_k+=a_strides_l[k]; if !(k in indices) Jl_k+=result_strides_l[k] end)) begin
                result[Jl_0, Jr_0] += a[Il_0, Ir_0]
            end
        end
        return result
    end
end

@generated function _ptrace_ket(::Type{Val{RANK}}, a,
                            shape, indices) where RANK
    return quote
        a_strides = _strides(shape)
        result_shape = copy(shape)
        @inbounds for idx ∈ indices
            result_shape[idx] = 1
        end
        result_strides = _strides(result_shape)
        N_result = prod(result_shape)
        result = zeros(eltype(a), N_result, N_result)
        @nexprs 1 (d->(Jr_{$RANK}=1;Ir_{$RANK}=1))
        @nloops $RANK ir (d->1:shape[d]) (d->(Ir_{d-1}=Ir_d; Jr_{d-1}=Jr_d)) (d->(Ir_d+=a_strides[d]; if !(d in indices) Jr_d+=result_strides[d] end)) begin
            @nexprs 1 (d->(Jl_{$RANK}=1;Il_{$RANK}=1))
            @nloops $RANK il (k->1:shape[k]) (k->(Il_{k-1}=Il_k; Jl_{k-1}=Jl_k; if (k in indices && il_k!=ir_k) Il_k+=a_strides[k]; continue end)) (k->(Il_k+=a_strides[k]; if !(k in indices) Jl_k+=result_strides[k] end)) begin
                result[Jl_0, Jr_0] += a[Il_0]*conj(a[Ir_0])
            end
        end
        return result
    end
end

@generated function _ptrace_bra(::Type{Val{RANK}}, a,
                            shape, indices) where RANK
    return quote
        a_strides = _strides(shape)
        result_shape = copy(shape)
        @inbounds for idx ∈ indices
            result_shape[idx] = 1
        end
        result_strides = _strides(result_shape)
        N_result = prod(result_shape)
        result = zeros(eltype(a), N_result, N_result)
        @nexprs 1 (d->(Jr_{$RANK}=1;Ir_{$RANK}=1))
        @nloops $RANK ir (d->1:shape[d]) (d->(Ir_{d-1}=Ir_d; Jr_{d-1}=Jr_d)) (d->(Ir_d+=a_strides[d]; if !(d in indices) Jr_d+=result_strides[d] end)) begin
            @nexprs 1 (d->(Jl_{$RANK}=1;Il_{$RANK}=1))
            @nloops $RANK il (k->1:shape[k]) (k->(Il_{k-1}=Il_k; Jl_{k-1}=Jl_k; if (k in indices && il_k!=ir_k) Il_k+=a_strides[k]; continue end)) (k->(Il_k+=a_strides[k]; if !(k in indices) Jl_k+=result_strides[k] end)) begin
                result[Jl_0, Jr_0] += conj(a[Il_0])*a[Ir_0]
            end
        end
        return result
    end
end

"""
    mul!(Y::DataOperator,A::AbstractOperator,B::DataOperator,alpha,beta) -> Y
    mul!(Y::StateVector,A::AbstractOperator,B::StateVector,alpha,beta) -> Y

Fast in-place multiplication of operators/state vectors. Updates `Y` as
`Y = alpha*A*B + beta*Y`. In most cases, the call gets forwarded to
Julia's 5-arg mul! implementation on the underlying data.
See also [`LinearAlgebra.mul!`](@ref).
"""
mul!(result::Operator{B1,B3},a::Operator{B1,B2},b::Operator{B2,B3},alpha,beta) where {B1,B2,B3} = (LinearAlgebra.mul!(result.data,a.data,b.data,alpha,beta); result)
mul!(result::Ket{B1},a::Operator{B1,B2},b::Ket{B2},alpha,beta) where {B1,B2} = (LinearAlgebra.mul!(result.data,a.data,b.data,alpha,beta); result)
mul!(result::Bra{B2},a::Bra{B1},b::Operator{B1,B2},alpha,beta) where {B1,B2} = (LinearAlgebra.mul!(result.data,transpose(b.data),a.data,alpha,beta); result)
rmul!(op::Operator, x) = (rmul!(op.data, x); op)

# Multiplication for Operators in terms of their gemv! implementation
function mul!(result::Operator{B1,B3},M::AbstractOperator{B1,B2},b::Operator{B2,B3},alpha,beta) where {B1,B2,B3}
    for i=1:size(b.data, 2)
        bket = Ket(b.basis_l, b.data[:,i])
        resultket = Ket(M.basis_l, result.data[:,i])
        mul!(resultket,M,bket,alpha,beta)
        result.data[:,i] = resultket.data
    end
    return result
end

function mul!(result::Operator{B1,B3},b::Operator{B1,B2},M::AbstractOperator{B2,B3},alpha,beta) where {B1,B2,B3}
    for i=1:size(b.data, 1)
        bbra = Bra(b.basis_r, vec(b.data[i,:]))
        resultbra = Bra(M.basis_r, vec(result.data[i,:]))
        mul!(resultbra,bbra,M,alpha,beta)
        result.data[i,:] = resultbra.data
    end
    return result
end

# Broadcasting
Base.size(A::DataOperator) = size(A.data)
Base.size(A::DataOperator, d) = size(A.data, d)
@inline Base.axes(A::DataOperator) = axes(A.data)
Base.broadcastable(A::DataOperator) = A

# Custom broadcasting styles
abstract type DataOperatorStyle{BL,BR} <: Broadcast.BroadcastStyle end
struct OperatorStyle{BL,BR} <: DataOperatorStyle{BL,BR} end

# Style precedence rules
Broadcast.BroadcastStyle(::Type{<:Operator{BL,BR}}) where {BL,BR} = OperatorStyle{BL,BR}()
Broadcast.BroadcastStyle(::OperatorStyle{B1,B2}, ::OperatorStyle{B3,B4}) where {B1,B2,B3,B4} = throw(IncompatibleBases())

# Out-of-place broadcasting
@inline function Base.copy(bc::Broadcast.Broadcasted{Style,Axes,F,Args}) where {BL,BR,Style<:OperatorStyle{BL,BR},Axes,F,Args<:Tuple}
    bcf = Broadcast.flatten(bc)
    bl,br = find_basis(bcf.args)
    bc_ = Broadcasted_restrict_f(bcf.f, bcf.args, axes(bcf))
    return Operator{BL,BR}(bl, br, copy(bc_))
end
find_basis(a::DataOperator, rest) = (a.basis_l, a.basis_r)

const BasicMathFunc = Union{typeof(+),typeof(-),typeof(*)}
function Broadcasted_restrict_f(f::BasicMathFunc, args::Tuple{Vararg{<:DataOperator}}, axes)
    args_ = Tuple(a.data for a=args)
    return Broadcast.Broadcasted(f, args_, axes)
end
function Broadcasted_restrict_f(f, args::Tuple{Vararg{<:DataOperator}}, axes)
    throw(error("Cannot broadcast function `$f` on type `$(eltype(args))`"))
end

# In-place broadcasting
@inline function Base.copyto!(dest::DataOperator{BL,BR}, bc::Broadcast.Broadcasted{Style,Axes,F,Args}) where {BL,BR,Style<:DataOperatorStyle{BL,BR},Axes,F,Args}
    axes(dest) == axes(bc) || Base.Broadcast.throwdm(axes(dest), axes(bc))
    # Performance optimization: broadcast!(identity, dest, A) is equivalent to copyto!(dest, A) if indices match
    if bc.f === identity && isa(bc.args, Tuple{<:DataOperator{BL,BR}}) # only a single input argument to broadcast!
        A = bc.args[1]
        if axes(dest) == axes(A)
            return copyto!(dest, A)
        end
    end
    # Get the underlying data fields of operators and broadcast them as arrays
    bcf = Broadcast.flatten(bc)
    bc_ = Broadcasted_restrict_f(bcf.f, bcf.args, axes(bcf))
    copyto!(dest.data, bc_)
    return dest
end
@inline Base.copyto!(A::DataOperator{BL,BR},B::DataOperator{BL,BR}) where {BL,BR} = (copyto!(A.data,B.data); A)
@inline Base.copyto!(dest::DataOperator{BL,BR}, bc::Broadcast.Broadcasted{Style,Axes,F,Args}) where {BL,BR,Style<:DataOperatorStyle,Axes,F,Args} =
    throw(IncompatibleBases())


text:
using AlphaZero

include("main.jl")
using .ConnectFour

computer = ConnectFour.Solver.Player()
human = Human{Game}()

interactive!(Game(), computer, human)


text:
# # This file includes the test set for ports

@testset "PortTestSet" begin
    @info "Running PortTestSet ..."

    # Constructiokn of Outport and Inport
    op = Outport(5)
    @test length(op.pins) == 5
    @test all(.!isbound.(op))
    @test eltype(op) == Outpin{Float64}
    op = Outport{Int}(5)
    @test eltype(op) == Outpin{Int}
    op = Outport()
    @test length(op.pins) == 1

    # Construction of Inport
    ip = Inport(3)
    @test length(ip.pins) == 3
    @test all(.!isbound.(ip))
    @test eltype(ip) == Inpin{Float64}
    ip = Inport{Bool}(4)
    @test eltype(ip) == Inpin{Bool}
    ip = Inport()
    @test length(ip.pins) == 1

    # Connection of Outport and Inport
    op, ip = Outport(3), Inport(3)
    ls = connect!(op, ip)
    @test typeof(ls) <: Vector{<:Link}
    @test length(ls) == 3
    for (l, _op, _ip) in zip(ls, op, ip)
        @test l.masterid == _op.id
        @test l.slaveid == _ip.id
    end
    @test isconnected(op, ip)

    # Partial connection of Outport and Inport
    op = Outport(5)
    ip1, ip2 = Inport(3), Inport(2)
    @test_throws DimensionMismatch connect!(op, ip1) # Length of op and ip1 are not same
    ls1 = connect!(op[1:3], ip1)
    ls2 = connect!(op[4:5], ip2)
    @test isconnected(op[1], ip1[1])
    @test isconnected(op[4], ip2[1])
    @test !isconnected(op[3], ip2[1])

    # Data transfer through ports.
    op, ip = Outport(2), Inport(2)
    @test_throws MethodError take!(op)
    @test_throws MethodError put!(ip, zeros(2))
    ls = connect!(op, ip)
    t = @async while true
        all(take!(ip) .=== NaN) && break
    end
    for val in 1 : 5
        put!(op, ones(2) * val)
        @test !istaskfailed(t)
    end
    put!(op, [NaN, NaN])
    @test istaskdone(t)

    # Disconnection of ports
    op, ip = Outport(), Inport()
    ls = connect!(op, ip)
    @test isconnected(op, ip)
    disconnect!(op, ip)
    @test !isconnected(op, ip)

    @info "Donke PortTestSet."
end  # testset


text:
#!/usr/bin/env julia
using DelimitedFiles: readdlm

function solve(case, R, C, grid)
    seen = BitMatrix(undef, size(grid)...)
    added = 0

    function propagate(I)
        checkbounds(Bool, grid, I) || return
        seen[I] && return

        maxc = zero(I)
        maxneighbor = 0
        neighbors = Ref(I) .+ CartesianIndex.([(0, 1), (1, 0), (0, -1), (-1, 0)])
        for _I in neighbors
            checkbounds(Bool, grid, _I) || continue
            candidate = grid[I] - grid[_I]
            if candidate > maxneighbor
                maxneighbor = candidate
                maxc = _I
            end
        end
        iszero(maxc) && return
        grid[maxc] = grid[I] - 1
        seen[I] = true
        added += abs(maxneighbor) - 1

        for neighbor in neighbors
            propagate(neighbor)
        end
    end

    maxval, _ = findmax(grid)
    maxinds = findall(==(maxval), grid)

    for maxind in maxinds
        propagate(maxind)
    end

    println("Case #$case: $added")
end

function main()
    T = parse(Int, readline(stdin))

    for case in 1:T
        R, C = map(x -> parse(Int, x), split(readline(stdin)))
        grid = Matrix{Int}(undef, R, C)
        for r in 1:R
            grid[r,:] .= map(x -> parse(Int, x), split(readline(stdin)))
        end
        solve(case, R, C, grid)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
using DelimitedFiles
using ZipFile

function conn76()
    fname = "conn76.zip"
    if ~isfile(fname)
        url = "https://github.com/the-virtual-brain/tvb-data/raw/master/tvb_data/connectivity/connectivity_76.zip"
        download(url, fname)
    end
    zf = ZipFile.Reader(fname)
    # @show zf.files
    l = readdlm(zf.files[6])
    w = readdlm(zf.files[7])
    close(zf)
    w, l
end


text:
################################################################################
## Defines linear system solvers.
## Note that linear system solvers need to be coupled with a KKT system solver
## these solvers can be found in the `kkt_system_solver` folder and
## reduces the linear system that needs to be solved, e.g., does a "primal schur complement"
################################################################################

println("Loading linear_system_solvers ... ")
@compat abstract type abstract_linear_system_solver end

include("julia.jl")
USE_HSL = true
if USE_HSL
	try
		include("hsl.jl")
	catch (e)
		println("Loading HSL failed:")
		warn(e)
		println("Continuing although you will not be able to choose HSL as a linear solver ...")
	end
end
#include("matlab.jl")
if USE_MUMPS
	include("mumps_wrapper.jl")
end

function initialize!(solver::abstract_linear_system_solver)

end

function finalize!(solver::abstract_linear_system_solver)

end

function inertia_status(pos_eigs::Int64, neg_eigs::Int64, zero_eigs::Int64, num_vars::Int64, num_constraints::Int64)
	###########################################################################################
	# INPUT:
	# pos_eigs, neg_eigs, zero_eigs = number of positive, negative and zero eigenvalues in the matrix:
	# [[ H A'];
	# [ A D ]];
	# num_vars = number of variables in problem
	# num_constraints = number of constraints in problem
	# OUTPUT:
	# is the inertia good, i.e., is H + A' D^{-1} A positive semidefinite?
	# this ensures that our direction will be a descent direction on the shifted log barrier.
	###########################################################################################

	try
		#println("inertia_status called") # ????
		if(pos_eigs + neg_eigs + zero_eigs != num_vars + num_constraints)
			println("pos_eigs = $pos_eigs")
			println("neg_eigs = $neg_eigs")
			println("zero_eigs = $pos_eigs")
			println("num_vars = $num_vars")
			println("num_constraints = $num_constraints")
			error("pos_eigs + neg_eigs + zero_eigs != num_vars + num_constraints")
		end
		# return number instead of true or false???

		if pos_eigs == num_vars && neg_eigs == num_constraints
			return true
		elseif pos_eigs > num_vars || neg_eigs + pos_eigs != num_vars + num_constraints
			# Warning: numerical instability in LDL factorization
			if pos_eigs > num_vars
				# more positive eigenvalues than variables
				return false
			else
				# zero eigenvalues
				return false
			end
		else
			return false
		end
	catch e
		println("ERROR in compute_inertia")
		throw(e)
	end
end


text:
using Test
using FFTResampling
using GeneralizedGenerated
using Base.Cartesian
using FFTW
using Random
Random.seed!(42)
include("utils.jl")
include("sinc_interpolation.jl")


text:
using ConcurrentCollections: DualLinkedConcurrentRingQueue
using ConcurrentCollectionsBenchmarks.BenchQueueHotPotato: hotpotato!, fai_stats
using JSON

function sweep(; repeat = 10, duration = 1, maxntasks = Threads.nthreads())
    # cooldown() = sleep(0.1)
    # cooldown() = GC.gc()
    cooldown() = nothing
    potatos = []
    fais = []
    for trialid in 1:repeat
        for ntasks in 1:maxntasks
            println(stderr, "Trial $trialid/$repeat #Tasks $ntasks/$maxntasks")

            cooldown()
            result = hotpotato!(DualLinkedConcurrentRingQueue{Bool}(); ntasks, duration)
            push!(potatos, (; trialid, result, ntasks, impl = :dlcrq))
            cooldown()
            result = hotpotato!(Channel{Bool}(Inf); ntasks, duration)
            push!(potatos, (; trialid, result, ntasks, impl = :base))

            cooldown()
            result = fai_stats(; ntasks, duration)
            push!(fais, (; trialid, result, ntasks))
        end
    end
    return (; potatos, fais, repeat, duration)
end

function git_info(dir = @__DIR__)
    git(cmd) = strip(read(setenv(`git $cmd`; dir), String))
    return (;
        revision = git(`rev-parse HEAD`),
        status = git(`status --short --untracked-files=no --porcelain`),
    )
end

function main(args = ARGS)
    output = get(args, 1, joinpath(@__DIR__, "build", "results.json"))
    mkpath(dirname(output))
    git = git_info()
    @info "Warmup..."
    sweep(; repeat = 1, duration = 0.1, maxntasks = 1)
    @info "Benchmarking..."
    results = sweep()
    results = (; results..., git)
    open(output, write = true) do io
        JSON.print(io, results)
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end