text:
#
# Contraction sequence
#

# Tools for contracting a network with a sequence
function _contract(label1::String, label2::String)
  return string("(", label1, "*", label2, ")")
end

function _contract(tensor1::ITensor, tensor2::ITensor)
  indsR = noncommoninds(tensor1, tensor2)
  return isempty(indsR) ? ITensor() : ITensor(indsR)
end

sequence_traversal(sequence) = reverse(collect(StatelessBFS(sequence)))

function contract_dict(tensors, sequence, traversal=sequence_traversal(sequence))
  net_tensors = Dict()
  traversal = reverse(collect(StatelessBFS(sequence)))
  for net in traversal
    if net isa Int
      net_tensors[net] = tensors[net]
    else # net isa Vector
      net_tensors[net] = _contract(net_tensors[net[1]], net_tensors[net[2]])
    end
  end
  return net_tensors
end

# Return all of the contractions involved in the sequence.
function contraction_sequence(
  tensors,
  sequence,
  traversal=sequence_traversal(sequence),
  contract_dict=contract_dict(tensors, sequence, traversal),
)
  all_tensors = Any[]
  tensors_1 = Vector{Union{Nothing,eltype(tensors)}}(tensors)
  net_position = Dict()
  N = length(tensors)
  n = N + 1
  for net in traversal
    if net isa Int
      net_position[net] = net
    else
      net_position[net] = n
      n += 1
    end
    if !isa(net, Int)
      for n in net
        tensors_1[net_position[n]] = nothing
      end
      push!(tensors_1, contract_dict[net])
      push!(all_tensors, copy(tensors_1))
    end
  end
  return convert.(Vector{eltype(tensors)}, filter.(!isnothing, all_tensors))
end

#
# Convert a tree to a graph
#

tree_to_graph(tr) = tree_to_graph(Tree(tr))

function tree_to_graph(tr::Tree)
  g = SimpleDiGraph()
  labels = Any[]
  walk_tree!(g, labels, tr)
  return (g, labels)
end

function walk_tree!(g, labels, tr::Tree)
  add_vertex!(g)
  top_vertex = vertices(g)[end]
  push!(labels, tr.x)
  for i in 1:length(tr.x)
    if isa(tr[i], Vector)
      child = walk_tree!(g, labels, Tree(tr[i]))
      add_edge!(g, child, top_vertex)
    else
      add_vertex!(g)
      n = vertices(g)[end]
      add_edge!(g, n, top_vertex)
      push!(labels, tr[i])
    end
  end
  return top_vertex
end

# Visualization function interface. Ultimately calls a beckend.

function visualize(g::AbstractGraph, sequence=nothing; backend=get_backend(), kwargs...)
  # TODO: do something with the sequence (show sequence, add labels indicating sequence, etc.)
  return visualize(Backend(backend), g; kwargs...)
end

function visualize(tn::Vector{ITensor}, sequence=nothing; kwargs...)
  return visualize(MetaDiGraph(tn), sequence; kwargs...)
end

function visualize(tn::Tuple{Vector{ITensor}}, args...; kwargs...)
  return visualize(only(tn), args...; kwargs...)
end
visualize(ψ::MPS, args...; kwargs...) = visualize(data(ψ), args...; kwargs...)
function visualize(tn::Tuple{Vararg{ITensor}}, args...; kwargs...)
  return visualize(collect(tn), args...; kwargs...)
end
visualize(tn::ITensor...; kwargs...) = visualize(collect(tn); kwargs...)

function visualize!(fig, g::AbstractGraph; backend=get_backend(), kwargs...)
  return visualize!(Backend(backend), fig, g; kwargs...)
end

function visualize!(fig, tn::Vector{ITensor}, sequence=nothing; kwargs...)
  return visualize!(fig, MetaDiGraph(tn); kwargs...)
end
visualize!(fig, ψ::MPS, sequence=nothing; kwargs...) = visualize!(fig, data(ψ); kwargs...)
function visualize!(fig, tn::Tuple{Vararg{ITensor}}, sequence=nothing; kwargs...)
  return visualize!(fig, collect(tn); kwargs...)
end
visualize!(fig, tn::ITensor...; kwargs...) = visualize!(fig, collect(tn); kwargs...)

function visualize!(fig, tn::Tuple{Vector{ITensor}}, sequence=nothing; kwargs...)
  return visualize!(fig, tn[1], sequence; kwargs...)
end
function visualize!(
  fig, f::Function, tn::Tuple{Vararg{ITensor}}, sequence=nothing; kwargs...
)
  return visualize!(fig, tn, sequence; kwargs...)
end

function visualize(f::Union{Function,Type}, As...; kwargs...)
  # TODO: specialize on the function type. Also accept a general collection.
  return visualize(As...; kwargs...)
end

function visualize!(fig, f::Union{Function,Type}, As...; kwargs...)
  # TODO: specialize of the function type. Also accept a general collection.
  return visualize!(fig, As...; kwargs...)
end

function _visualize_sequence!(fig, tn, sequence, n; kwargs...)
  return error("Not implemented")
end

function sequence_labels(sequence, all_sequences, vertex_labels)
  traversal = sequence_traversal(sequence)
  labels_dict = contract_dict(vertex_labels, sequence, traversal)
  all_labels = [labels_dict[s] for s in all_sequences]
  return all_labels
end

function _graphplot(backend::Backend, graph; all_labels)
  return error("Not implemented for backend $backend.")
end

function visualize_sequence(sequence, vertex_labels)
  graph, all_sequences = tree_to_graph(sequence)
  all_labels = sequence_labels(sequence, all_sequences, vertex_labels)
  fig = _graphplot(Backend"Makie"(), graph; all_labels=all_labels)
  return fig
end

function default_sequence(tn::Vector{ITensor})
  N = length(tn)
  return foldl((x, y) -> [x, y], 1:N)
end

function visualize_sequence(
  f::Union{Function,Type}, tn::Vector{ITensor}, sequence::Nothing; kwargs...
)
  return visualize_sequence(f, tn, default_sequence(tn); kwargs...)
end

function visualize_sequence(
  f::Union{Function,Type}, tn::Vector{ITensor}, sequence=default_sequence(tn); kwargs...
)
  N = length(tn)

  # TODO: clean this up a bit
  vertex_labels_prefix = get(
    kwargs,
    :vertex_labels_prefix,
    default_vertex_labels_prefix(Backend("Makie"), MetaDiGraph(tn)),
  )
  vertex_labels = get(
    kwargs,
    :vertex_labels,
    default_vertex_labels(Backend(""), MetaDiGraph(tn), vertex_labels_prefix),
  )

  fig = visualize_sequence(sequence, vertex_labels)

  visualize!(fig[1, 2], tn; vertex_labels=vertex_labels, kwargs...)

  traversal = sequence_traversal(sequence)
  labels_sequence = contraction_sequence(vertex_labels, sequence, traversal)

  tn_sequence = contraction_sequence(tn, sequence, traversal)

  for n in 1:length(tn_sequence)
    visualize!(fig[1, n + 2], tn_sequence[n]; vertex_labels=labels_sequence[n], kwargs...)
  end

  return fig
end

function visualize_sequence(
  f::Union{Function,Type}, tn::Tuple{Vector{ITensor}}, sequence; kwargs...
)
  return visualize_sequence(f, tn[1], sequence; kwargs...)
end


text:
using ArrayInterface
using ArrayInterfaceStaticArrays
using LinearAlgebra
using StaticArrays
using Static
using Test

x = @SVector [1,2,3]
@test @inferred(ArrayInterface.device(typeof(x))) === ArrayInterface.CPUTuple()

x = @MVector [1,2,3]
@test @inferred(ArrayInterface.device(typeof(x))) === ArrayInterface.CPUPointer()

@test isone(ArrayInterface.known_first(typeof(StaticArrays.SOneTo(7))))
@test ArrayInterface.known_last(typeof(StaticArrays.SOneTo(7))) == 7
@test ArrayInterface.known_length(typeof(StaticArrays.SOneTo(7))) == 7

@test ArrayInterface.parent_type(SizedVector{1, Int, Vector{Int}}) <: Vector{Int}
@test ArrayInterface.known_length(@inferred(ArrayInterface.indices(SOneTo(7)))) == 7

x = view(SArray{Tuple{3,3,3}}(ones(3,3,3)), :, SOneTo(2), 2)
@test @inferred(ArrayInterface.known_length(x)) == 6
@test @inferred(ArrayInterface.known_length(x')) == 6

v = @SVector rand(8);
A = @MMatrix rand(7, 6);
T = SizedArray{Tuple{5,4,3}}(zeros(5,4,3));
@test @inferred(ArrayInterface.length(v)) === StaticInt(8)
@test @inferred(ArrayInterface.length(A)) === StaticInt(42)
@test @inferred(ArrayInterface.length(T)) === StaticInt(60)

Am = @MMatrix rand(2,10);
@test @inferred(ArrayInterface.strides(view(Am,1,:))) === (StaticInt(2),)

@test @inferred(ArrayInterface.contiguous_axis(@SArray(zeros(2,2,2)))) === ArrayInterface.StaticInt(1)
@test @inferred(ArrayInterface.contiguous_axis_indicator(@SArray(zeros(2,2,2)))) == (true,false,false)
@test @inferred(ArrayInterface.contiguous_batch_size(@SArray(zeros(2,2,2)))) === ArrayInterface.StaticInt(0)
@test @inferred(ArrayInterface.stride_rank(@SArray(zeros(2,2,2)))) == (1, 2, 3)
@test @inferred(ArrayInterface.is_column_major(@SArray(zeros(2,2,2)))) === True()
@test @inferred(ArrayInterface.dense_dims(@SArray(zeros(2,2,2)))) == (true,true,true)

ext: jl


text:
using MolSimPy
using Test
using ASE
using JuLIP.Testing

@testset "MolSimPy.jl" begin
   include("matscipy.jl")
end


text:
const ad = AutoHOOT.autodiff
const go = AutoHOOT.graphops

function get_symbol(i::Int)
  if i < 26
    return 'a' + i
  elseif i < 52
    return 'A' + i - 26
  end
  return Char(i + 140)
end

"""
Retrieve the key from the dictionary that maps AutoHOOT nodes
to tensors. Returns Nothing if key not exists.
"""
retrieve_key(dict, value) = findfirst(v -> (v === value), dict)

"""Compute the computational graph defined in AutoHOOT.
Parameters
----------
outnodes: A list of AutoHOOT einsum nodes
node_dict: A dictionary mapping AutoHOOT node to tensor
Returns
-------
A list of tensors.
"""
function compute_graph!(out_nodes, node_dict; kwargs...)
  topo_order = ad.find_topo_sort(out_nodes)
  for node in topo_order
    if haskey(node_dict, node) == false && node.name != "1.0"
      input_list = []
      for n in node.inputs
        if haskey(node_dict, n)
          push!(input_list, node_dict[n])
        else
          # the scalar that can neglect
          @assert(n.name == "1.0")
        end
      end
      input_list = Vector{typeof(input_list[1])}(input_list)
      node_dict[node] = contract(input_list; kwargs...)
    end
  end
  return [node_dict[node] for node in out_nodes]
end

function compute_graph(out_nodes, node_dict; kwargs...)
  return compute_graph!(out_nodes, copy(node_dict); kwargs...)
end

"""Extract an tensor network from an input network based on AutoHOOT einsum tree.
The input network is defined by the tensors in node_dict.
Note: this function ONLY extracts network based on the input nodes rather than einstr.
The output network can be hierarchical.
Parameters
----------
outnode: AutoHOOT einsum node
node_dict: A dictionary mapping AutoHOOT node to tensor
Returns
-------
A list representing the output network.
Example
-------
>>> extract_network(einsum("ab,bc->ac", A, einsum("bd,dc->bc", B, C)),
                    Dict(A => A_tensor, B => B_tensor, C => C_tensor))
>>> [A_tensor, [B_tensor, C_tensor]]
"""
function extract_network(out_node, node_dict)
  topo_order = ad.find_topo_sort([out_node])
  node_dict = copy(node_dict)
  for node in topo_order
    if haskey(node_dict, node) == false && node.name != "1.0"
      input_list = []
      for n in node.inputs
        if haskey(node_dict, n)
          push!(input_list, node_dict[n])
        else
          # the scalar that can neglect
          @assert(n.name == "1.0")
        end
      end
      node_dict[node] = input_list
    end
  end
  return node_dict[out_node]
end

# perform inner product between two nodes
function inner(n1, n2)
  @assert(n1.shape == n2.shape)
  str = join([get_symbol(i) for i in 1:length(n1.shape)], "")
  return ad.einsum(str * "," * str * "->", n1, n2)
end

"""Generate AutoHOOT einsum expression based on a list of input networks
Parameters
----------
network_list: An array of networks. Each network is represented by an array of tensors
Returns
-------
A list of AutoHOOT einsum node;
A dictionary mapping AutoHOOT input node to tensor
"""
function generate_einsum_expr(
  network_list::Vector{<:Vector{<:AbstractTensor}}; optimize=false
)
  node_dict = Dict()
  outnodes = [
    generate_einsum_expr!(network, node_dict; optimize=optimize) for network in network_list
  ]
  return outnodes, node_dict
end

function generate_einsum_expr(trees::Vector{SubNetwork}; optimize=false)
  node_dict = Dict()
  outnodes = [generate_einsum_expr!(tree, node_dict; optimize=optimize) for tree in trees]
  return outnodes, node_dict
end

function generate_einsum_expr!(
  network::Vector{<:AbstractTensor}, node_dict::Dict; optimize=false
)
  input_nodes = input_nodes_generation!(network, node_dict)
  einstr = einstr_generation(network)
  out = ad.einsum(einstr, input_nodes...)
  return optimize ? go.generate_optimal_tree(out) : out
end

function generate_einsum_expr!(tree::SubNetwork, node_dict::Dict; optimize=false)
  input_nodes = []
  for input in tree.inputs
    if input isa AbstractTensor
      node = input_nodes_generation!([input], node_dict)[1]
    else
      node = generate_einsum_expr!(input, node_dict; optimize=optimize)
    end
    push!(input_nodes, node)
  end
  einstr = einstr_generation(tree.inputs)
  out = ad.einsum(einstr, input_nodes...)
  return optimize ? go.generate_optimal_tree(out) : out
end

function update_dict!(node_dict::Dict, tensor)
  i = length(node_dict) + 1
  if length(inds(tensor)) != 0
    nodename = "tensor" * string(i)
    shape = [space(index) for index in inds(tensor)]
    node = ad.Variable(nodename; shape=shape)
  else
    node = ad.scalar(tensor[])
  end
  node_dict[node] = tensor
  return node
end

"""Generate AutoHOOT nodes based on tensors
Parameters
----------
network: An array of tensors
node_dict: A dictionary mapping AutoHOOT node to tensor.
node_dict will be inplace updated.
Returns
-------
node_list: An array of AutoHOOT nodes
"""
function input_nodes_generation!(network::Array, node_dict::Dict)
  node_list = []
  for (i, tensor) in enumerate(network)
    node = retrieve_key(node_dict, tensor)
    if node == nothing
      node = update_dict!(node_dict, tensor)
    end
    push!(node_list, node)
  end
  return node_list
end

function einstr_generation(network::Array)
  index_dict = Dict{Index{Int64},Char}()
  label_num = 0
  string_list = []
  # build input string
  for tensor in network
    str = ""
    for index in inds(tensor)
      if haskey(index_dict, index) == false
        index_dict[index] = get_symbol(label_num)
        label_num += 1
      end
      str = str * index_dict[index]
    end
    push!(string_list, str)
  end
  instr = join(string_list, ",")
  # build output string
  output_inds = noncommoninds(network...)
  outstr = join([index_dict[i] for i in output_inds])
  return instr * "->" * outstr
end

function generate_node_index_dict(node_dict, networks)
  node_index_dict = Dict()
  for i in 1:length(networks)
    for j in 1:length(networks[i])
      node = retrieve_key(node_dict, networks[i][j])
      node_index_dict[node] = (i => j)
    end
  end
  return node_index_dict
end


text:
using Base.Test
@everywhere using DSGE, JLD

my_tests = [
            "core",
            "parameters",
            "models/an_schorfheide/an_schorfheide",
            "models/smets_wouters/smets_wouters",
            "models/m990/m990",
            "models/m1002/m1002",
            "models/m1010/m1010",
            "data/misc",
            "data/load_data",
            "solve/gensys",
            "solve/solve",
            "estimate/filter",
            "estimate/cat",
            "estimate/posterior",
            "estimate/hessizero",
            "estimate/hessian",
            "estimate/csminwel",
            "estimate/optimize",
            "estimate/metropolis_hastings",
            "forecast/smooth",
            "forecast/forecast",
            "forecast/shock_decompositions",
            "forecast/impulse_responses",
            "forecast/io",
            "forecast/forecast_one",
            "analysis/compute_meansbands",
            "altpolicy/altpolicy",
            "scenarios/scenario",
            "scenarios/forecast",
            "scenarios/switching",
            "scenarios/drivers",
            "plot/plot"
            ]

for test in my_tests
    test_file = string("$test.jl")
    @printf " * %s\n" test_file
    include(test_file)
end


text:
"""
```
function gensys_cplus(m::AbstractDSGEModel, Γ0::Matrix{Float64}, Γ1::Matrix{Float64}, C::Vector{Float64},
                      Ψ::Matrix{Float64}, Π::Matrix{Float64}, TTT::Matrix{Float64}, RRR::Matrix{Float64},
                      CCC::Vector{Float64})
```

calculates the state space transition matrices when a temporary alternative policy applies.
"""
function gensys_cplus(m::AbstractDSGEModel, Γ0::Matrix{Float64}, Γ1::Matrix{Float64}, C::Vector{Float64},
                      Ψ::Matrix{Float64}, Π::Matrix{Float64}, TTT::Matrix{Float64}, RRR::Matrix{Float64},
                      CCC::Vector{Float64})
    T_switch = get_setting(m, :n_rule_periods) + 1 #get_setting(m, :n_regimes) - 1

    exp_eq_ind = sum(Π, dims = 2)
    Γ0_til = zeros(size(Γ0))
    Γ1_til = zeros(size(Γ1))
    Γ2_til = zeros(size(Γ0))
    C_til = C
    Ψ_til = Ψ

    for row in 1:length(exp_eq_ind)
        if exp_eq_ind[row] == 0
            # Not expectational equation
            Γ1_til[row, :] = Γ1[row, :]
            Γ0_til[row, :] = Γ0[row, :]
            Γ2_til[row, :] .= 0.
        else
            # expectational equation
            Γ0_til[row, findfirst(Γ1[row, :] .> 0)] = -1
            Γ2_til[row, findfirst(Γ0[row, :] .> 0)] = 1
        end
    end

#     Γ0_til, Γ1_til, Γ2_til, C_til, Ψ_til = gensys_to_predictable_form(Γ0, Γ1, C, Ψ, Π)

    Tcal = Vector{Matrix{Float64}}(undef, T_switch)
    Rcal = Vector{Matrix{Float64}}(undef, T_switch)
    Ccal = Vector{Vector{Float64}}(undef, T_switch)

    Tcal[end] = TTT
    Rcal[end] = RRR
    Ccal[end] = CCC

    for t = 1:(T_switch-1)
        Tcal[end-t] = (Γ2_til*TTT + Γ0_til)\Γ1_til
        Rcal[end-t] = (Γ2_til*TTT + Γ0_til)\Ψ_til
        Ccal[end-t] = (Γ2_til*TTT + Γ0_til)\(C_til - Γ2_til*CCC)

        TTT = Tcal[end-t]
        CCC = Ccal[end-t]
        RRR = Rcal[end-t]
    end
    return Tcal, Rcal, Ccal
end

function gensys_cplus(m::AbstractDSGEModel, Γ0s::Vector{Matrix{Float64}}, Γ1s::Vector{Matrix{Float64}},
                      Cs::Vector{Vector{Float64}},
                      Ψs::Vector{Matrix{Float64}}, Πs::Vector{Matrix{Float64}},
                      TTT::Matrix{Float64}, RRR::Matrix{Float64},
                      CCC::Vector{Float64})
    T_switch = get_setting(m, :n_rule_periods) + 1
    Γ0_tils = Vector{Matrix{Float64}}(undef, length(Γ0s))
    Γ1_tils = Vector{Matrix{Float64}}(undef, length(Γ0s))
    Γ2_tils = Vector{Matrix{Float64}}(undef, length(Γ0s))
    C_tils = Vector{Vector{Float64}}(undef, length(Γ0s))
    Ψ_tils = Vector{Matrix{Float64}}(undef, length(Γ0s))

    for i in 1:length(Γ0s)
        exp_eq_ind = sum(Πs[i], dims = 2)
        Γ0_tils[i] = zeros(size(Γ0s[i]))
        Γ1_tils[i] = zeros(size(Γ1s[i]))
        Γ2_tils[i] = zeros(size(Γ0s[i]))
        C_tils[i] = Cs[i]
        Ψ_tils[i] = Ψs[i]

        for row in 1:length(exp_eq_ind)
            if exp_eq_ind[row] == 0
                # Not expectational equation
                Γ1_tils[i][row, :] = Γ1s[i][row, :]
                Γ0_tils[i][row, :] = Γ0s[i][row, :]
                Γ2_tils[i][row, :] .= 0.
            else
                # expectational equation
                Γ0_tils[i][row, findfirst(Γ1s[i][row, :] .> 0)] = -1
                Γ2_tils[i][row, findfirst(Γ0s[i][row, :] .> 0)] = 1
            end
        end
    end

#=    for i in 1:length(Γ0s)
        Γ0_tils[i], Γ1_tils[i], Γ2_tils[i], C_tils[i], Ψ_tils[i] = gensys_to_predictable_form(Γ0s[i], Γ1s[i], Cs[i], Ψs[i], Πs[i])
    end=#

    Tcal = Vector{Matrix{Float64}}(undef, T_switch)
    Rcal = Vector{Matrix{Float64}}(undef, T_switch)
    Ccal = Vector{Vector{Float64}}(undef, T_switch)

    Tcal[end] = TTT
    Rcal[end] = RRR
    Ccal[end] = CCC

    for t = 1:(T_switch-1)
        Tcal[end-t] = (Γ2_tils[end-t]*TTT + Γ0_tils[end-t])\Γ1_tils[end-t]
        Rcal[end-t] = (Γ2_tils[end-t]*TTT + Γ0_tils[end-t])\Ψ_tils[end-t]
        Ccal[end-t] = (Γ2_tils[end-t]*TTT + Γ0_tils[end-t])\(C_tils[end-t] - Γ2_tils[end-t]*CCC)

        TTT = Tcal[end-t]
        CCC = Ccal[end-t]
        RRR = Rcal[end-t]
    end
    return Tcal, Rcal, Ccal
end

function gensys_to_predictable_form(Γ0::AbstractMatrix{S}, Γ1::AbstractMatrix{S}, C::AbstractVector{S},
                                    Ψ::AbstractMatrix{S}, Π::AbstractMatrix{S}) where {S <: Real}

    nonexp_eq_ind = map(i -> all(Π[i, :] .== 0.), 1:size(Π, 1))
    Γ0_til = zeros(S, size(Γ0))
    Γ1_til = zeros(S, size(Γ1))
    Γ2_til = zeros(S, size(Γ0))
    C_til  = C
    Ψ_til  = Ψ

    for row in 1:length(nonexp_eq_ind)
        if nonexp_eq_ind[row]
            # Not expectational equation
            Γ1_til[row, :]  = Γ1[row, :]
            Γ0_til[row, :]  = Γ0[row, :]
        else
            # expectational equation, assumed to take the form zₜ = Eₜ₋₁[zₜ] + ηₜ
            Γ0_til[row, findfirst(Γ1[row, :] .> 0)] = -1.
            Γ2_til[row, findfirst(Γ0[row, :] .> 0)] = 1.
        end
    end

    return Γ0_til, Γ1_til, Γ2_til, C_til, Ψ_til
end

#=

function gensys_cplus_old(m::AbstractDSGEModel, Γ0::Matrix{Float64}, Γ1::Matrix{Float64}, C::Vector{Float64},
                      Ψ::Matrix{Float64}, TTT::Matrix{Float64}, RRR::Matrix{Float64},
                      CCC::Vector{Float64})
    T_switch = 2 #get_setting(m, :T_switch)
    nstates = n_states(m)
    nant = n_mon_anticipated_shocks(m)
    n_exo = length(m.exogenous_shocks)
    n_shocks = n_exo - nant
    ind_states_exp = findall(x -> occursin("E", string(x)), collect(keys(m.endogenous_states)))
    ind_eq_exp = findall(x -> occursin("E", string(x)), collect(keys(m.equilibrium_conditions)))

    non_exp_exp_keys = map(x-> Symbol(string(x)[2:end]), collect(keys(m.endogenous_states))[ind_states_exp])
    non_exp_exp_inds = findall(x->x in non_exp_exp_keys, collect(keys(m.endogenous_states)))

    n_expected = length(ind_states_exp)

    num_z = nstates - n_expected #n_end+n_exo
    # Initializing the Γ2 matrix
    Γ2_tilde = zeros(nstates, nstates)

    # taking the Expectations
    Γ0_tilde = copy(Γ0)
    Γ0_tilde[ind_eq_exp,:] = Γ0_tilde[ind_eq_exp,:]*0.
    Γ0_tilde[:,ind_states_exp] = Γ0_tilde[:,ind_states_exp]*0.

    # moving the expectation from Γ0 to Γ2
  #  for i_ind = 1:n_expected
    Γ2_tilde[:, non_exp_exp_inds] = Γ0[:, ind_eq_exp]

    # identity map of expectational states
    Γ2_tilde[ind_eq_exp, non_exp_exp_inds] .= 1
    Γ0_tilde[ind_eq_exp, ind_states_exp] .= -1
  #  end

    #taking expectations
    #Γ2_tilde(num_z+1:end-1,:) = []
    #Γ2_tilde(:,num_z+1:end-1) = []

    #taking the Expectations
    Γ1_tilde = Γ1
    Γ1_tilde[ind_eq_exp,:] .= 0 #[]
    Γ1_tilde[:,ind_states_exp] .= 0 #[]
    # Taking expectation
    Ψ_tilde = Ψ
    Ψ_tilde[ind_eq_exp, :] =. 0 #[]
    Ψ_tilde[:,n_shocks+1:end] =. 0 #[] #taking out anticipated policy shocks
    # taking expectations out of C
    C_tilde = C #[C(1:num_z)C(end)]

    ## This solves for the TTT and RRR Matrix using the methods in "Solving
    ## Linear Rational Expectations Methods with Anticipated Policy Changes
    # taking out expectation equations
    T11 = TTT #[Not(ind_eq_exp), Not(ind_states_exp)]
    R11 = RRR #[Not(ind_eq_exp), 1:n_shocks]
    C11 = CCC #[Not(ind_eq_exp)]
#=    T11 = TTT
    T11[ind_eq_exp,:] .= NaN
    T11[:,ind_state_exp] .= NaN

    R11 = RRR
    R11[ind_eq,:] .= NaN
    R11[:,nshocks+1:end] .= NaN  #taking out anticipated policy shocks

    C11 = CCC
    C11[ind_eq_exp] .= NaN =#

    #initialize Tcal and Rcal matrices
    Tcal = zeros(size(T11, 1), size(T11, 2), T_switch+1) #zeros(num_z+1, num_z+1, T_switch+1)
    Rcal = zeros(size(R11, 1), size(R11, 2), T_switch+1) #num_z+1, nshocks, T_switch+1)
    Ccal = zeros(length(C11), T_switch+1) #num_z+1, T_switch+1)

    Tcal[:,:,end] = T11
    Rcal[:,:,end] = R11
    Ccal[:,end] = C11

    for t = 1:T_switch
        Tcal[:,:,end-t] = (Γ2_tilde*T11 + Γ0_tilde)\Γ1_tilde
        Rcal[:,:,end-t] = (Γ2_tilde*T11 + Γ0_tilde)\Ψ_tilde
        Ccal[:,end-t] = (Γ2_tilde*T11 + Γ0_tilde)\(C_tilde - Γ2_tilde*C11)

        # Tcal(:,:,end-t) = pinv(Γ2_tilde*T11 + Γ0_tilde)*Γ1_tilde
        # Rcal(:,:,end-t) = pinv(Γ2_tilde*T11 + Γ0_tilde)*PSI_tilde
        # Ccal(:,end-t) = pinv(Γ2_tilde*T11 + Γ0_tilde)*(C_tilde - Γ2_tilde*C11)

        T11 = Tcal[:,:,end-t]
    end
        C11 = Ccal[:,end-t]
        R11= Rcal[:,:,end-t]
    end
end
=#


text:
# Legendre-Gauss quadrature weights and nodes
# This code was originally written for Matlab
#  by Greg von Winckel - 02/25/2004
# posted on Mathworks website and licensed under BSD License
# http://fr.mathworks.com/matlabcentral/fileexchange/4540-legendre-gauss-quadrature-weights-and-nodes
#
#

function lgwt(N,a,b)

	N=N-1;
	N1=N+1; N2=N+2;
	xu= transpose(range(-1,stop = 1, length = N1));

	y = transpose(cos.((2*transpose(0:N).+1)*pi/(2*N+2))+(0.27/N1)*sin.(pi*xu*N/N2));

	L = zeros(N1,N2);

	Lp=zeros(N1,N2);
	y0=2;

	while maximum(abs.(y.-y0))>eps()
		L[:,1] .= 1;
		Lp[:,1] .= 0;
		L[:,2] = y;
		for k=2:N1
			L[:,k+1]=( (2*k-1)*y.*L[:,k]-(k-1)*L[:,k-1] )/k;
		end
		Lp=(N2)*( L[:,N1] - y .* L[:,N2] )./(1 .- y.^2);
		y0=y;
		y=y0-L[:,N2]./Lp;
	end

	x = (a*(1 .- y)+b*(1 .+ y))/2;

	w=(b-a)./((1 .- y.^2).*Lp.^2)*(N2/N1)^2;

	return x,w
end

max_stars_repo_name: flavioluiz/phjulia

max_issues_repo_name: flavioluiz/phjulia

max_forks_repo_name: flavioluiz/phjulia


text:
function OpenMapRealMatrix(arg0::OpenMapRealMatrix)
    return OpenMapRealMatrix((OpenMapRealMatrix,), arg0)
end

function OpenMapRealMatrix(arg0::jint, arg1::jint)
    return OpenMapRealMatrix((jint, jint), arg0, arg1)
end

function add(obj::OpenMapRealMatrix, arg0::OpenMapRealMatrix)
    return jcall(obj, "add", OpenMapRealMatrix, (OpenMapRealMatrix,), arg0)
end

function add_to_entry(obj::OpenMapRealMatrix, arg0::jint, arg1::jint, arg2::jdouble)
    return jcall(obj, "addToEntry", void, (jint, jint, jdouble), arg0, arg1, arg2)
end

function copy(obj::OpenMapRealMatrix)
    return jcall(obj, "copy", OpenMapRealMatrix, ())
end

function create_matrix(obj::OpenMapRealMatrix, arg0::jint, arg1::jint)
    return jcall(obj, "createMatrix", OpenMapRealMatrix, (jint, jint), arg0, arg1)
end

function get_column_dimension(obj::OpenMapRealMatrix)
    return jcall(obj, "getColumnDimension", jint, ())
end

function get_entry(obj::OpenMapRealMatrix, arg0::jint, arg1::jint)
    return jcall(obj, "getEntry", jdouble, (jint, jint), arg0, arg1)
end

function get_row_dimension(obj::OpenMapRealMatrix)
    return jcall(obj, "getRowDimension", jint, ())
end

function multiply(obj::OpenMapRealMatrix, arg0::OpenMapRealMatrix)
    return jcall(obj, "multiply", OpenMapRealMatrix, (OpenMapRealMatrix,), arg0)
end

function multiply(obj::OpenMapRealMatrix, arg0::RealMatrix)
    return jcall(obj, "multiply", RealMatrix, (RealMatrix,), arg0)
end

function multiply_entry(obj::OpenMapRealMatrix, arg0::jint, arg1::jint, arg2::jdouble)
    return jcall(obj, "multiplyEntry", void, (jint, jint, jdouble), arg0, arg1, arg2)
end

function multiply_transposed(obj::OpenMapRealMatrix, arg0::RealMatrix)
    return jcall(obj, "multiplyTransposed", RealMatrix, (RealMatrix,), arg0)
end

function set_entry(obj::OpenMapRealMatrix, arg0::jint, arg1::jint, arg2::jdouble)
    return jcall(obj, "setEntry", void, (jint, jint, jdouble), arg0, arg1, arg2)
end

function subtract(obj::OpenMapRealMatrix, arg0::OpenMapRealMatrix)
    return jcall(obj, "subtract", OpenMapRealMatrix, (OpenMapRealMatrix,), arg0)
end

function subtract(obj::OpenMapRealMatrix, arg0::RealMatrix)
    return jcall(obj, "subtract", OpenMapRealMatrix, (RealMatrix,), arg0)
end

function transpose_multiply(obj::OpenMapRealMatrix, arg0::RealMatrix)
    return jcall(obj, "transposeMultiply", RealMatrix, (RealMatrix,), arg0)
end


text:
module PolyBench

using BenchmarkTools

const SUITE = BenchmarkGroup()
const BENCH_DIRS = ["datamining",
                    "linear-algebra/blas",
                    "linear-algebra/kernels",
                    "linear-algebra/solvers",
                    "medley",
                    "stencils"]

function load!()
    for bench_dir in BENCH_DIRS
        bench_dir_path = joinpath(dirname(@__FILE__), bench_dir)
        for bench_file in filter(x -> endswith(x, "jl"), readdir(bench_dir_path))
            bench_file_path = joinpath(bench_dir_path, bench_file)
            bench_name = bench_file[1:end-3]
            print("loading $bench_name...")
            tic()
            eval(PolyBench, :(include($bench_file_path)))
            load_time = toq()
            println("done (took $load_time) seconds)")
        end
    end
    return SUITE
end

end # module


text:
using Random
using DABenchmarks
using LaTeXStrings

Random.seed!(0)

using Distributed
machine = [("MACHINENAME", 16)]
addprocs(machine, dir = "/home/USER", exename="/path/to/julia", exeflags="--project")
# addprocs(16, exeflags="--project") # For local workers

@everywhere using StochasticPrograms
@everywhere using LShapedSolvers
@everywhere using HydroModels
@everywhere using Gurobi

function prepare_dbenchmark(nscenarios::Integer, nsamples::Integer; timeout::Integer = 1000)
    dls = LShapedSolver(GurobiSolver(OutputFlag=0), log=false, distributed = true, subsolver = ()->GurobiSolver(OutputFlag=0), κ=1.0)
    dtr = LShapedSolver(GurobiSolver(OutputFlag=0), regularization = :tr, log=false, distributed = true, subsolver = ()->GurobiSolver(OutputFlag=0), κ=1.0)
    solvers = [dls,dtr]
    solvernames = ["Distributed L-shaped","Distributed L-shaped with trust-region"]
    # Create Day-ahead benchmark
    return Scaling(solvernames, solvers, nscenarios, nsamples; timeout = timeout)
end
@info "Running L-shaped benchmarks on 16 workers..."
dab = prepare_dbenchmark(1000, 10); warmup_benchmarks!(dab); run_benchmarks!(dab); save_results!(dab, "ls_scaling_16.json")
exit()

max_stars_repo_name: martinbiel/mpc-reproducibility

max_issues_repo_name: martinbiel/mpc-reproducibility

max_forks_repo_name: martinbiel/mpc-reproducibility


text:
using DataFrames, Distributions, Gadfly

set_default_plot_size(6inch, 3inch)

beta = Beta(2,2)
Dbeta = DataFrame(x=rand(beta, 10^4))
layer1 = layer(x->pdf(beta, x), 0, 1, Theme(default_color="black"))
p1 = plot(Dbeta, layer1,
    layer(x=:x, Geom.histogram(bincount=20, density=true, limits=(min=0,))),
)
p2 = plot(Dbeta, layer1,
    layer(x=:x, Geom.histogram(bincount=20, density=true, limits=(min=0, max=1)))
)

hstack(p1, p2)


text:
using DiffEqBase, Test

function f_lin(du,u,p,t)
  du[1] = 0.2u[1] + p[1]*u[2]
  du[2] = 0.2u[1] - p[2]*u[2]
end
p = (0.0,1.0)
prob = LinearProblem(f_lin,ones(2))
prob = LinearProblem(rand(2,2),ones(2))

function f_nonlin(du,u,p)
  du[1] = 0.2u[1] + p[1]*u[2]
  du[2] = 0.2u[1] - p[2]*u[2]
end
p = (0.0,1.0)
prob = NonlinearProblem(f_nonlin,ones(2),p)

function f_quad(du,u,p)
  du[1] = 0.2u[1] + p[1]*u[2]
  du[2] = 0.2u[1] - p[2]*u[2]
end
p = (0.0,1.0)
prob = QuadratureProblem(f_quad,zeros(2),ones(2),p)

function f(du,u,p,t)
  du[1] = 0.2u[1]
  du[2] = 0.4u[2]
end
u0 = ones(2)
tspan = (0,1.0)

prob = ODEProblem(f,u0,tspan)
@test typeof(prob.tspan) == Tuple{Float64,Float64}
prob = ODEProblem{true}(f,u0,tspan)
@test typeof(prob.tspan) == Tuple{Float64,Float64}
prob = ODEProblem(ODEFunction{true}(f),u0,tspan)
@test typeof(prob.tspan) == Tuple{Float64,Float64}
@test isinplace(prob) == true
prob = ODEProblem{false}(f,u0,tspan)
@test isinplace(prob) == false

@test_broken @inferred ODEProblem{true}(f,u0,tspan)
@test_broken @inferred ODEProblem(f,u0,tspan)

function f(dv,u,v,p,t)
  dv .= 2.0 .* v
end
u0 = ones(2)
v0 = ones(2)
tspan = (0,1.0)
prob = SecondOrderODEProblem(f,u0,v0,tspan)

prob = SDEProblem((u,p,t) -> 1.01u,(u,p,t) -> 0.87u,1/2,(0.0,1.0))

function f(du,u,p,t)
  du[1] = 0.2u[1]
  du[2] = 0.4u[2]
end
function g(du,u,p,t)
  du[1] = 0.2u[1]
  du[2] = 0.4u[2]
end
u0 = ones(2)
tspan = (0,1.0)
prob = SDEProblem(f,g,u0,tspan)
prob = SDEProblem{true}(f,g,u0,tspan)

@test_broken @inferred SDEProblem(f,g,u0,tspan)
@test_broken @inferred SDEProblem{true}(f,g,u0,tspan)

f_1delay = function (du,u,h,p,t)
    du[1] = - h(t-1)[1]
end
prob =  DDEProblem(f_1delay,ones(1),t->zeros(1),(0.0, 10.0),constant_lags = ones(1))
prob =  DDEProblem{true}(f_1delay,ones(1),t->zeros(1),(0.0, 10.0),dependent_lags = ones(1))

@test_broken @inferred DDEProblem(f_1delay,ones(1),t->zeros(1),(0.0, 10.0),constant_lags = ones(1))
@test_broken @inferred DDEProblem{true}(f_1delay,ones(1),t->zeros(1),(0.0, 10.0),dependent_lags = ones(1))

function f(r, yp, y, p,tres)
    r[1]  = -0.04*y[1] + 1.0e4*y[2]*y[3]
    r[2]  = -r[1] - 3.0e7*y[2]*y[2] - yp[2]
    r[1] -=  yp[1]
    r[3]  =  y[1] + y[2] + y[3] - 1.0
end
u0 = [1.0, 0, 0]
du0 = [-0.04, 0.04, 0.0]
prob_dae_resrob = DAEProblem(f,du0,u0,(0.0,100000.0))
prob_dae_resrob = DAEProblem{true}(f,du0,u0,(0.0,100000.0))

@test_broken @inferred DAEProblem(f,du0,u0,(0.0,100000.0))
@test_broken @inferred DAEProblem{true}(f,du0,u0,(0.0,100000.0))

f(u,t,W) = 1.01u.+0.87u.*W
u0 = 1.00
tspan = (0.0,1.0)
prob = RODEProblem(f,u0,tspan)
prob = RODEProblem{false}(f,u0,tspan)

@test_broken @inferred RODEProblem(f,u0,tspan)
@test_broken @inferred RODEProblem{false}(f,u0,tspan)

DiscreteProblem(ones(1),tspan)
f(t,u) = 0.5
DiscreteProblem{false}(f,ones(1),tspan)

@test_broken @inferred DiscreteProblem(ones(1),tspan)
@test_broken @inferred DiscreteProblem{false}(f,ones(1),tspan)

function f(du,u,t)
  du[1] = 2 - 2u[1]
  du[2] = u[1] - 4u[2]
end
u0 = zeros(2)
prob = SteadyStateProblem(f,u0)

@test_broken @inferred SteadyStateProblem(f,u0)
@test SteadyStateProblem(ODEProblem(f,u0,tspan,:param)).p == :param


text:
using StructJuMP

include("farmer_models.jl")

m = farmer_stochastic()
DSPopt.classifyConstrs(m)

@testset "Parent model" begin
    @test length(m.variables) == 3
    @test length(m.constraints) == 1
    start, index, value, rlbd, rubd, obj, clbd, cubd, ctype, cname, nqrows, linnzcnt, quadnzcnt, rhs, sense, linstart, linind, linval, quadstart, quadrow, quadcol, quadval = DSPopt.get_model_data(m)
    @test length(start) == length(m.constraints) + 1
    @test start[end] == length(index)
    @test start == [0, 3]
    @test index == [0, 1, 2]
    @test value == [1., 1., 1.]
    @test rlbd == [-Inf]
    @test rubd == [Budget]
    @test obj == Cost
    @test clbd == zeros(3)
    @test cubd == zeros(3) .+ Inf
    @test ctype == "III"
end

@testset "Child model $i" for (i,subm) in m.children
    @test i > 0
    @test m.probability[i] == probability[i]
    @test length(subm.variables) == 6
    @test length(subm.constraints) == 4
    start, index, value, rlbd, rubd, obj, clbd, cubd, ctype, cname, nqrows, linnzcnt, quadnzcnt, rhs, sense, linstart, linind, linval, quadstart, quadrow, quadcol, quadval = DSPopt.get_model_data(subm, i)
    @test length(start) == length(subm.constraints) + 1
    @test start[end] == length(index)
    @test start == [0, 3, 6, 9, 10]
    @test index == [0, 3, 5, 1, 4, 6, 2, 7, 8, 7]
    @test value == [Yield[i,1], 1., -1., Yield[i,2], 1., -1., Yield[i,3], -1., -1., 1.]
    @test rlbd == [Minreq; -Inf]
    @test rubd == [Inf, Inf, Inf, 6000]
    @test obj == [Purchase; -Sell]
    @test clbd == zeros(6)
    @test cubd == zeros(6) .+ Inf
    @test ctype == "CCCCCC"
end

DSPopt.set(WassersteinSet, 2, 1.0)

@testset "optimize!: $j" for j in [DSPopt.Dual, DSPopt.Benders]
    status = DSPopt.optimize!(m, is_stochastic = true, solve_type = j)
    @test DSPopt.termination_status(m) == MOI.OPTIMAL
    DSPopt.freeSolver(dsp)
    DSPopt.freeModel(dsp)
end


text:
## Auxiliary functions for adjusting the drift or leap step in geometric MCMC
## algorithms

# Generic function for step adjustment
function setStep(i::Int, acceptanceRatio::Float64, nMcmc::Int, nBurnin::Int,
  currentStep::Float64, steps::Vector{Float64}, lowerRatio::Vector{Float64},
  upperRatio::Vector{Float64})
  if i == 1
    steps[1]
  elseif i <= 0.1*nBurnin
    if acceptanceRatio < lowerRatio[1]
      steps[2]
    elseif upperRatio[1] <= acceptanceRatio
      steps[3]
    else
      currentStep
    end
  elseif i <= 0.15*nBurnin
    if acceptanceRatio < lowerRatio[2]
      steps[3]
    elseif upperRatio[2] <= acceptanceRatio
      steps[4]
    else
      currentStep
    end
  elseif i <= 0.2*nBurnin
    if acceptanceRatio < lowerRatio[3]
      steps[4]
    elseif upperRatio[3] <= acceptanceRatio
      steps[5]
    else
      currentStep
    end
  elseif i <= 0.25*nBurnin
    if acceptanceRatio < lowerRatio[4]
      steps[5]
    elseif upperRatio[4] <= acceptanceRatio
      steps[6]
    else
      currentStep
    end
  elseif i <= 0.3*nBurnin
    if acceptanceRatio < lowerRatio[5]
      steps[6]
    elseif upperRatio[5] <= acceptanceRatio
      steps[7]
    else
      currentStep
    end
  elseif i <= 0.35*nBurnin
    if acceptanceRatio < lowerRatio[6]
      steps[7]
    elseif upperRatio[6] <= acceptanceRatio
      steps[8]
    else
      currentStep
    end
  elseif i <= 0.4*nBurnin
    if acceptanceRatio < lowerRatio[7]
      steps[8]
    elseif upperRatio[7] <= acceptanceRatio
      steps[9]
    else
      currentStep
    end
  elseif i <= 0.45*nBurnin
    if acceptanceRatio < lowerRatio[8]
      steps[9]
    elseif upperRatio[8] <= acceptanceRatio
      steps[10]
    else
      currentStep
    end
  elseif i <= 0.5*nBurnin
    if acceptanceRatio < lowerRatio[10]
      throw("Aborted: low acceptance ratio during burn-in")
    elseif acceptanceRatio < lowerRatio[9]
      steps[10]
    elseif upperRatio[9] <= acceptanceRatio
      steps[1]
    else
      currentStep
    end
  elseif i <= 0.7*nBurnin
    currentStep
  elseif i <= 0.75*nBurnin
    if acceptanceRatio < lowerRatio[10]
      throw("Aborted: low acceptance ratio during burn-in")
    else
      currentStep
    end
  elseif i <= 0.9*nBurnin
    currentStep
  elseif i <= nBurnin
    if acceptanceRatio < lowerRatio[10]
      throw("Aborted: low acceptance ratio during burn-in")
    elseif upperRatio[10] <= acceptanceRatio
      throw("Aborted: high acceptance ratio during burn-in")
    else
      currentStep
    end
  # 0.5*(nBurnin+nMcmc) =nBurnin+0.5*nPostBurnin
  elseif i <= 0.5*(nBurnin+nMcmc)
    currentStep
  elseif i <= nMcmc
    if acceptanceRatio < lowerRatio[10]
      throw("Aborted: low acceptance ratio during post burn-in")
    elseif upperRatio[10] <= acceptanceRatio
      throw("Aborted: high acceptance ratio during post burn-in")
    else
      currentStep
    end
  end
end

# Function for adjusting MALA drift step
malaLowerRatio = [0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.1]
malaUpperRatio = [0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.85, 0.8, 0.9]

setMalaDriftStep(i::Int, acceptanceRatio::Float64, nMcmc::Int, nBurnin::Int,
  currentStep::Float64, steps::Vector{Float64}) =
  setStep(i::Int, acceptanceRatio::Float64, nMcmc::Int, nBurnin::Int,
  currentStep::Float64, steps::Vector{Float64}, malaLowerRatio, malaUpperRatio)

# Function for adjusting SMMALA drift step
smmalaLowerRatio = copy(malaLowerRatio)
smmalaUpperRatio = copy(malaUpperRatio)

setSmmalaDriftStep(i::Int, acceptanceRatio::Float64, nMcmc::Int, nBurnin::Int,
  currentStep::Float64, steps::Vector{Float64}) =
  setStep(i::Int, acceptanceRatio::Float64, nMcmc::Int, nBurnin::Int,
  currentStep::Float64, steps::Vector{Float64}, smmalaLowerRatio,
  smmalaUpperRatio)

# Function for adjusting MMALA drift step
mmalaLowerRatio = copy(malaLowerRatio)
mmalaUpperRatio = copy(malaUpperRatio)

setMmalaDriftStep(i::Int, acceptanceRatio::Float64, nMcmc::Int, nBurnin::Int,
  currentStep::Float64, steps::Vector{Float64}) =
  setStep(i::Int, acceptanceRatio::Float64, nMcmc::Int, nBurnin::Int,
  currentStep::Float64, steps::Vector{Float64}, mmalaLowerRatio,
  mmalaUpperRatio)

# Function for adjusting HMC leap step
hmcLowerRatio = [0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5, 0.5, 0.1]
hmcUpperRatio = [0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.85, 0.8, 0.925]

setHmcLeapStep(i::Int, acceptanceRatio::Float64, nMcmc::Int, nBurnin::Int,
  currentStep::Float64, steps::Vector{Float64}) =
  setStep(i::Int, acceptanceRatio::Float64, nMcmc::Int, nBurnin::Int,
  currentStep::Float64, steps::Vector{Float64}, hmcLowerRatio, hmcUpperRatio)

# Function for adjusting RMHMC leap step
rmhmcLowerRatio = [0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5, 0.6, 0.5]
rmhmcUpperRatio =
  [0.975, 0.975, 0.975, 0.975, 0.975, 0.975, 0.975, 0.95, 0.95, 0.975]

setRmhmcLeapStep(i::Int, acceptanceRatio::Float64, nMcmc::Int, nBurnin::Int,
  currentStep::Float64, steps::Vector{Float64}) =
  setStep(i::Int, acceptanceRatio::Float64, nMcmc::Int, nBurnin::Int,
  currentStep::Float64, steps::Vector{Float64}, rmhmcLowerRatio,
  rmhmcUpperRatio)


text:
module TestElectrical

using Modia
using Modia.Electric

# Desired:
#   using ModiaMath: plot
#   using Test
#
# In order that these packages need not to be defined in the user environment, they are included via Modia:
using Modia.ModiaMath: plot

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Modia.Test
end

result = simulate(Resistor, 1)
@test result["i"][end] == 0.0
@test result["v"][end] == 0.0

@model ParallelResistors begin
    R1 = Resistor(R=1, p=Pin(v=Float(start=0.0)), n=Pin(v=Float(start=0.0)))
    R2 = Resistor(R=1)
    @equations begin
        connect(R1.p, R2.p)
        connect(R1.n, R2.n)
    end
end

result = simulate(ParallelResistors, 1)
@test result["R1.i"][end] == 0.0
@test result["R1.v"][end] == 0.0

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

# result = simulate(Capacitor, 1)

@model ParallelCapacitors begin
    C1 = Capacitor(C=2, v=Float(start=1.0))
    C2 = Capacitor(C=1, v=Float(state=false))
    ground = Ground()
    @equations begin
        connect(C1.p, C2.p)
        connect(C1.n, C2.n)
        connect(C1.n, ground.p)
    end
end

result = simulate(ParallelCapacitors, 1)
@test result["C2.v"][end] == 1.0

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

# result = simulate(Inductor, 1)

@model InductorsInSeries begin
    V = ConstantVoltage(V=10)
    L1 = Inductor(L=1, i=Current(start=0.0))
    R1 = Resistor(R=100) # , n=Pin(v=Float(start=0.0)))
    R2 = Resistor(R=100) # , n=Pin(v=Float(start=0.0)))
    L2 = Inductor(L=1, i=Current(start=0.0))
    @equations begin
        connect(V.p, L1.p)
        connect(L1.n, R1.p)
        connect(L1.n, R2.p)
        connect(R1.n, L2.p)
        connect(R2.n, L2.p)
        connect(L2.n, V.n)
    end
end

# result = simulate(InductorsInSeries, 1)

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

@model ParallelCapacitorCircuit begin
    R = Resistor(R=100.0)
    C = Capacitor(C=2.5E-2, v=Float(start=0.0))
    C2 = Capacitor(C=2.5E-2, v=Float(state=false))
    V = StepVoltage(startTime=10, V=100)
    #  Vt=ConstantVoltage(V=10.0)
    ground = Ground()
    @equations begin
        connect(V.n, ground.p)
        # connect(Vt.n, ground.p)
        connect(V.p, R.p)
        connect(R.n, C.p)
        connect(C.n, V.n)
        connect(C.n, C2.n)
        connect(C.p, C2.p)
    end
end

# simulate(ParallelCapacitorCircuit, 1, useKinsol=true, removeSingularities=false, logTranslation=true)
# checkSimulation(ParallelCapacitorCircuit, 1, "C1.v", 1.0, useKinsol=true, removeSingularities=false)

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineReimageParameters <: SwaggerModel
    tempDisk::Any # spec type: Union{ Nothing, Bool } # spec name: tempDisk

    function VirtualMachineReimageParameters(;tempDisk=nothing)
        o = new()
        validate_property(VirtualMachineReimageParameters, Symbol("tempDisk"), tempDisk)
        setfield!(o, Symbol("tempDisk"), tempDisk)
        o
    end
end # type VirtualMachineReimageParameters

const _property_map_VirtualMachineReimageParameters = Dict{Symbol,Symbol}(Symbol("tempDisk")=>Symbol("tempDisk"))
const _property_types_VirtualMachineReimageParameters = Dict{Symbol,String}(Symbol("tempDisk")=>"Bool")
Base.propertynames(::Type{ VirtualMachineReimageParameters }) = collect(keys(_property_map_VirtualMachineReimageParameters))
Swagger.property_type(::Type{ VirtualMachineReimageParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineReimageParameters[name]))}
Swagger.field_name(::Type{ VirtualMachineReimageParameters }, property_name::Symbol) =  _property_map_VirtualMachineReimageParameters[property_name]

function check_required(o::VirtualMachineReimageParameters)
    true
end

function validate_property(::Type{ VirtualMachineReimageParameters }, name::Symbol, val)
end


text:
#=
-----------------------------------------------------------------------------------------------
This program  constructs points and weights using:
 1) Chebyshev,
 2) Legendre,
 3) Lobatto, and
 4) Equi-spaced points

 This constitutes Project 1 and the results are described in Secs. 3.4 and 4.4 in the textbook.
 Written by F.X. Giraldo
            Department of Applied Mathematics
            Naval Postgraduate School
            Monterey, California 93943-5216
-----------------------------------------------------------------------------------------------
=#

using Plots
include("QuadraturePoints.jl")

#Some Constants
DFloat = Float64
#----------------------------------------Only Change this---------------------------------------
N=4
#-----------------------------------------------------------------------------------------------

function main()

    #Allocate arrays
    ξarray=zeros(DFloat,N+1,4)
    ωarray=zeros(DFloat,N+1,4)

    #Some Output
    @show (N)

    #Chebyshev Points
    (ξ,ω) = QuadraturePoints.chebyshev_gauss(N+1)
    println("Chebyshev Points ") #output
    @show (ξ,ω)
    ξarray[:,1]=ξ
    ωarray[:,1]=ω

    #Legendre Points
    (ξ,ω) = QuadraturePoints.legendre_gauss(N+1)
    println("Legendre Points ") #output
    @show (ξ,ω)
    ξarray[:,2]=ξ
    ωarray[:,2]=ω

    #Lobatto Points
    (ξ,ω) = QuadraturePoints.lobatto_gauss(N+1)
    println("Lobatto Points ") #output
    @show (ξ,ω)
    ξarray[:,3]=ξ
    ωarray[:,3]=ω

    #Equispaced Points
    (ξ,ω) = QuadraturePoints.equispaced_points(N+1)
    println("Equispaced Points ") #output
    @show (ξ,ω)
    ξarray[:,4]=ξ
    ωarray[:,4]=ω

    #Plot Example
    plot_handle=plot(ξarray,ωarray,xlabel="Root",ylabel="Weight",legend=true,lw=3,label=["Chebyshev" "Legendre" "Lobatto" "Equispaced"],title="Roots and Weights")
    display(plot_handle)
    savefig(plot_handle,"QuadraturePoints.png")

    #Plot Interpolation
    println("Done") #output

end

#Call Main Program
main()


text:
using Blink, Mux

function collab(p)
  loadeditor(p, ver = "4.3.0")
  load!(p, "https://cdn.firebase.com/js/client/2.2.4/firebase.js")
  load!(p, "https://cdn.firebase.com/libs/firepad/1.1.0/firepad.min.js")
  @js_ p begin
    @var firebase = @new Firebase("crackling-heat-4207.firebaseio.com")
    @var firepad = Firepad.fromCodeMirror(firebase, cm)
  end
end

@app demo =
  (Mux.defaults,
   Blink.resroute,
   page("/juno", Blink.app(collab)),
   Mux.notfound())

@app demow =
  (Mux.wdefaults,
   Blink.ws_handler)

rundemo() = serve(demo, demow, 8080)


text:
using LazyArrays, Test
import LazyArrays: ApplyStyle, VectorSetStyle, SetStyle

@testset "Set operations" begin
    @test ApplyStyle(union,Int) == ApplyStyle(intersect,Int) == ApplyStyle(setdiff,Int) == VectorSetStyle()
    @test ApplyStyle(union,Set) == ApplyStyle(intersect,Set) == ApplyStyle(setdiff,Set) == SetStyle()
    @test 1 ∈ applied(union,1,2.0)
    @test 2 ∈ applied(union,1,2.0)
    @test 3 ∉ applied(union,1,2.0)
    @test issubset([1,2],  applied(union,1,2.0))

    @test eltype(apply(∪, 1, 2.0)) == Float64
    @test  apply(∪, 1, 2.0) == [1.0,2.0]
    @test apply(∪, [1,2], Set([2,4])) ==  Set([2,4]) ∪ [1,2] == apply(union, Set([2,4]) , [1,2]) ==
        apply(union, [1,2], Set([2,4])) == Set([1,2,4])
    @test Base.union([1,2], Set([2,4])) == Base.:∪([1,2], Set([2,4])) ==
                [1,2,4]

    @test 1 ∈ applied(intersect, [1,2], Set([1.0]))
    @test 2 ∉ applied(intersect, [1,2], Set([1.0]))
    @test 1 ∉ applied(setdiff, [1,2], Set([1.0]))
    @test 2 ∈ applied(setdiff, [1,2], Set([1.0]))

    @test 1 ∈ applied(intersect, [1,2], Set([1.0]), 1f0)
end


text:
import jInv.ForwardShare.getSensMat
function getSensMat(m::Vector,pFor::DirectParam)
	return pFor.Jacobian;
end

export getSensMatVec
function getSensMatVec(v::Vector,m::Vector,pFor::DirectParam)
	S = pFor.Jacobian;
	d = S'*v[:];

end

export getSensTMatVec
function getSensTMatVec(v::Vector,m::Vector,pFor::DirectParam)
	JtV = 0;
	println("len of v:",size(v))
	println("len of m:",size(m))

	S = pFor.Jacobian;
	println("len of J:",size(S))
	JtV = S'*v[:];
return JtV;
end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
function get24_hours_kp(obj::DTM2000InputParameters, arg0::AbsoluteDate)
    return jcall(obj, "get24HoursKp", jdouble, (AbsoluteDate,), arg0)
end

function get_instant_flux(obj::DTM2000InputParameters, arg0::AbsoluteDate)
    return jcall(obj, "getInstantFlux", jdouble, (AbsoluteDate,), arg0)
end

function get_max_date(obj::DTM2000InputParameters)
    return jcall(obj, "getMaxDate", AbsoluteDate, ())
end

function get_mean_flux(obj::DTM2000InputParameters, arg0::AbsoluteDate)
    return jcall(obj, "getMeanFlux", jdouble, (AbsoluteDate,), arg0)
end

function get_min_date(obj::DTM2000InputParameters)
    return jcall(obj, "getMinDate", AbsoluteDate, ())
end

function get_three_hourly_kp(obj::DTM2000InputParameters, arg0::AbsoluteDate)
    return jcall(obj, "getThreeHourlyKP", jdouble, (AbsoluteDate,), arg0)
end


text:
using Blink
Blink.AtomShell.isinstalled() || Blink.AtomShell.install()


text:
function advance!(pac::T, pap) where {T<:P_common{<:FdtdElastic}}
    for ipw in pac.activepw
        # store p for the last two steps
        # pppppp!(pap[ipw],pac.attrib_mod)
        advance_kernel!(pap[ipw], pac)
    end
    return nothing
end

# these relative indices of the arrays point to same location
# [ix,iy,iz]     --> tauxx, tauyy, and tauzz grid [1:nx,1:ny,1:nz]
# [ix+1/2,iy+1/2,iz+1]      --> tauxy
# [ix+1/2,iy+1,iz+1/2]      -->  tauxz
# [ix-1/2,iy,iz]       --> vx
# [ix,iy-1/2,iz]       --> vy
# [ix,iy,iz-1/2]       --> vz

@parallel function compute_dtau!(
    tauxx::Data.Array{3},
    tauyy,
    tauzz,
    tauxy,
    tauxz,
    tauyz,
    dtauxxdx,
    dtauxydx,
    dtauxzdx,
    dtauyydy,
    dtauxydy,
    dtauyzdy,
    dtauzzdz,
    dtauyzdz,
    dtauxzdz,
    dxI,
    dyI,
    dzI,
)

    @all(dtauxxdx) = @d_xi(tauxx) * dxI # at [ix+1/2,iy+1,iz+1] with indices []
    @all(dtauxydx) = @d_xa(tauxy) * dxI # at [ix+1,iy+1/2,iz+1]
    @all(dtauxzdx) = @d_xa(tauxz) * dxI # at [ix+1,iy+1,iz+1/2]

    @all(dtauyydy) = @d_yi(tauyy) * dyI # at [ix+1,iy+1/2,iz+1]
    @all(dtauxydy) = @d_ya(tauxy) * dyI # at [ix+1/2,iy+1,iz+1]
    @all(dtauyzdy) = @d_ya(tauyz) * dyI # at [ix+1,iy+1,iz+1/2]

    @all(dtauzzdz) = @d_zi(tauzz) * dzI # at [ix+1,iy+1,iz+1/2]
    @all(dtauxzdz) = @d_za(tauxz) * dzI # at [ix+1/2,iy+1,iz+1]
    @all(dtauyzdz) = @d_za(tauyz) * dzI # at [ix+1,iy+1/2,iz+1]

    return
end

@parallel function compute_dtau!(
    tauxx::Data.Array{2},
    tauzz,
    tauxz,
    dtauxxdx,
    dtauxzdx,
    dtauzzdz,
    dtauxzdz,
    dxI,
    dzI,
)

    @all(dtauxxdx) = @d_xi(tauxx) * dxI # at [ix+1/2,iz+1] with indices []
    @all(dtauxzdx) = @d_xa(tauxz) * dxI # at [ix+1,iz+1/2]

    @all(dtauzzdz) = @d_zi(tauzz) * dzI # at [ix+1,iz+1/2]
    @all(dtauxzdz) = @d_za(tauxz) * dzI # at [ix+1/2,iz+1]

    return
end

@parallel function compute_v!(
    vx::Data.Array{3},
    vy,
    vz,
    dtauxxdx,
    dtauxydx,
    dtauxzdx,
    dtauyydy,
    dtauxydy,
    dtauyzdy,
    dtauzzdz,
    dtauyzdz,
    dtauxzdz,
    dt,
    rho,
)
    @inn(vx) =
        @inn(vx) - dt / @av_xi(rho) * (@all(dtauxxdx) + @all(dtauxydy) + @all(dtauxzdz))
    @inn(vy) =
        @inn(vy) - dt / @av_yi(rho) * (@all(dtauxydx) + @all(dtauyydy) + @all(dtauyzdz))
    @inn(vz) =
        @inn(vz) - dt / @av_zi(rho) * (@all(dtauxzdx) + @all(dtauyzdy) + @all(dtauzzdz))

    return
end
@parallel function compute_v!(
    vx::Data.Array{2},
    vz,
    dtauxxdx,
    dtauxzdx,
    dtauzzdz,
    dtauxzdz,
    dt,
    rho,
)
    @inn(vx) = @inn(vx) - dt / @av_xi(rho) * (@all(dtauxxdx) + @all(dtauxzdz))
    @inn(vz) = @inn(vz) - dt / @av_zi(rho) * (@all(dtauxzdx) + @all(dtauzzdz))

    return
end

@parallel function compute_dv!(
    vx::Data.Array{3},
    vy,
    vz,
    dvxdx,
    dvydy,
    dvzdz,
    dvxdy,
    dvxdz,
    dvydx,
    dvydz,
    dvzdx,
    dvzdy,
    dxI,
    dyI,
    dzI,
)
    @all(dvxdx) = @d_xa(vx) * dxI # at [ix,iy,iz]
    @all(dvydy) = @d_ya(vy) * dyI # at      "
    @all(dvzdz) = @d_za(vz) * dzI # at      "

    @all(dvxdy) = @d_yi(vx) * dyI # at [ix+1/2,iy+1/2,iz+1]
    @all(dvxdz) = @d_zi(vx) * dzI # at [ix+1/2,iy+1,iz+1/2]

    @all(dvydz) = @d_zi(vy) * dzI # at [ix+1,iy+1/2,iz+1/2]
    @all(dvydx) = @d_xi(vy) * dxI # at [ix+1/2,iy+1/2,iz+1]

    @all(dvzdx) = @d_xi(vz) * dxI # at [ix+1/2,iy+1,iz+1/2]
    @all(dvzdy) = @d_yi(vz) * dyI # at [ix+1,iy+1/2,iz+1/2]

    return
end

@parallel function compute_dv!(vx::Data.Array{2}, vz, dvxdx, dvzdz, dvxdz, dvzdx, dxI, dzI)
    @all(dvxdx) = @d_xa(vx) * dxI # at [ix,iz]
    @all(dvzdz) = @d_za(vz) * dzI # at      "

    @all(dvxdz) = @d_zi(vx) * dzI # at [ix+1/2,iz+1/2]
    @all(dvzdx) = @d_xi(vz) * dxI # at [ix+1/2,iy+1,iz+1/2]

    return
end

@parallel function compute_tauii!(
    tauxx::Data.Array{3},
    tauyy,
    tauzz,
    dvxdx,
    dvydy,
    dvzdz,
    dt,
    M,
    lambda,
)

    @all(tauxx) =
        @all(tauxx) -
        dt * ((@all(M) * @all(dvxdx)) + (@all(lambda) * (@all(dvydy) + @all(dvzdz))))
    @all(tauyy) =
        @all(tauyy) -
        dt * ((@all(M) * @all(dvydy)) + (@all(lambda) * (@all(dvxdx) + @all(dvzdz))))
    @all(tauzz) =
        @all(tauzz) -
        dt * ((@all(M) * @all(dvzdz)) + (@all(lambda) * (@all(dvydy) + @all(dvxdx))))

    return
end

@parallel function compute_tauii!(tauxx::Data.Array{2}, tauzz, dvxdx, dvzdz, dt, M, lambda)
    @all(tauxx) =
        @all(tauxx) - dt * ((@all(M) * @all(dvxdx)) + (@all(lambda) * (@all(dvzdz))))
    @all(tauzz) =
        @all(tauzz) - dt * ((@all(M) * @all(dvzdz)) + (@all(lambda) * (@all(dvxdx))))

    return
end
@parallel function compute_tauij!(
    tauxy::Data.Array{3},
    tauxz,
    tauyz,
    dvxdy,
    dvxdz,
    dvydx,
    dvydz,
    dvzdx,
    dvzdy,
    dt,
    mu,
)
    @all(tauxz) = @all(tauxz) - dt * (@av_xzi(mu) * (@all(dvxdz) + @all(dvzdx)))
    @all(tauxy) = @all(tauxy) - dt * (@av_yzi(mu) * (@all(dvxdy) + @all(dvydx)))
    @all(tauyz) = @all(tauyz) - dt * (@av_xyi(mu) * (@all(dvydz) + @all(dvzdy)))

    return
end
@parallel function compute_tauij!(tauxz::Data.Array{2}, dvxdz, dvzdx, dt, mu)
    @all(tauxz) = @all(tauxz) - dt * (@av(mu) * (@all(dvxdz) + @all(dvzdx)))
    return
end

@parallel_indices (iy, ix) function free_surface!(tau::Data.Array{3}, izfree)
    tau[izfree, iy, ix] = 0.0
    return
end
@parallel_indices (iy, ix) function free_surface_mirror!(tau::Data.Array{3}, izfree)
    tau[izfree-1, iy, ix] = -1 * tau[izfree, iy, ix]
    return
end
@parallel_indices (ix) function free_surface!(tau::Data.Array{2}, izfree)
    tau[izfree, ix] = 0.0
    return
end
@parallel_indices (ix) function free_surface_mirror!(tau::Data.Array{2}, izfree)
    tau[izfree-1, ix] = -1 * tau[izfree, ix]
    return
end

function advance_kernel!(pap, pac::T) where {T<:P_common{FdtdElastic,3}}
    w1t = pap.w1[:t]
    mw = pap.memory_pml
    pml = pac.pml

    @parallel compute_dtau!(
        w1t[:tauxx],
        w1t[:tauyy],
        w1t[:tauzz],
        w1t[:tauxy],
        w1t[:tauxz],
        w1t[:tauyz],
        w1t[:dtauxxdx],
        w1t[:dtauxydx],
        w1t[:dtauxzdx],
        w1t[:dtauyydy],
        w1t[:dtauxydy],
        w1t[:dtauyzdy],
        w1t[:dtauzzdz],
        w1t[:dtauyzdz],
        w1t[:dtauxzdz],
        pac.fc[:dxI],
        pac.fc[:dyI],
        pac.fc[:dzI],
    )

    memoryx!(
        mw[:dtauxxdx],
        w1t[:dtauxxdx],
        pml[:x][:a_half],
        pml[:x][:b_half],
        pml[:x][:k_halfI],
    )
    memory1x!(mw[:dtauxydx], w1t[:dtauxydx], pml[:x][:a], pml[:x][:b], pml[:x][:kI])
    memory1x!(mw[:dtauxzdx], w1t[:dtauxzdx], pml[:x][:a], pml[:x][:b], pml[:x][:kI])

    memoryy!(
        mw[:dtauyydy],
        w1t[:dtauyydy],
        pml[:y][:a_half],
        pml[:y][:b_half],
        pml[:y][:k_halfI],
    )
    memory1y!(mw[:dtauxydy], w1t[:dtauxydy], pml[:y][:a], pml[:y][:b], pml[:y][:kI])
    memory1y!(mw[:dtauyzdy], w1t[:dtauyzdy], pml[:y][:a], pml[:y][:b], pml[:y][:kI])

    memoryz!(
        mw[:dtauzzdz],
        w1t[:dtauzzdz],
        pml[:z][:a_half],
        pml[:z][:b_half],
        pml[:z][:k_halfI],
    )
    memory1z!(mw[:dtauyzdz], w1t[:dtauyzdz], pml[:z][:a], pml[:z][:b], pml[:z][:kI])
    memory1z!(mw[:dtauxzdz], w1t[:dtauxzdz], pml[:z][:a], pml[:z][:b], pml[:z][:kI])

    @parallel compute_v!(
        w1t[:vx],
        w1t[:vy],
        w1t[:vz],
        w1t[:dtauxxdx],
        w1t[:dtauxydx],
        w1t[:dtauxzdx],
        w1t[:dtauyydy],
        w1t[:dtauxydy],
        w1t[:dtauyzdy],
        w1t[:dtauzzdz],
        w1t[:dtauyzdz],
        w1t[:dtauxzdz],
        pac.fc[:dt],
        pac.mod[:rho],
    )

    @parallel (1:pac.ic[:nz], 1:pac.ic[:ny]) dirichletx!(
        w1t[:vx],
        w1t[:vz],
        w1t[:vy],
        pac.ic[:nx],
    )
    @parallel (1:pac.ic[:nz], 1:pac.ic[:nx]) dirichlety!(
        w1t[:vy],
        w1t[:vz],
        w1t[:vx],
        pac.ic[:ny],
    )
    @parallel (1:pac.ic[:ny], 1:pac.ic[:nx]) dirichletz!(
        w1t[:vz],
        w1t[:vy],
        w1t[:vx],
        pac.ic[:nz],
    )
    @parallel compute_dv!(
        w1t[:vx],
        w1t[:vy],
        w1t[:vz],
        w1t[:dvxdx],
        w1t[:dvydy],
        w1t[:dvzdz],
        w1t[:dvxdy],
        w1t[:dvxdz],
        w1t[:dvydx],
        w1t[:dvydz],
        w1t[:dvzdx],
        w1t[:dvzdy],
        pac.fc[:dxI],
        pac.fc[:dyI],
        pac.fc[:dzI],
    )
    memoryx!(mw[:dvxdx], w1t[:dvxdx], pml[:x][:a], pml[:x][:b], pml[:x][:kI])
    memoryy!(mw[:dvydy], w1t[:dvydy], pml[:y][:a], pml[:y][:b], pml[:y][:kI])
    memoryz!(mw[:dvzdz], w1t[:dvzdz], pml[:z][:a], pml[:z][:b], pml[:z][:kI])

    memoryy!(mw[:dvxdy], w1t[:dvxdy], pml[:y][:a_half], pml[:y][:b_half], pml[:y][:k_halfI])
    memoryz!(mw[:dvxdz], w1t[:dvxdz], pml[:z][:a_half], pml[:z][:b_half], pml[:z][:k_halfI])
    memoryx!(mw[:dvydx], w1t[:dvydx], pml[:x][:a_half], pml[:x][:b_half], pml[:x][:k_halfI])

    memoryz!(mw[:dvydz], w1t[:dvydz], pml[:z][:a_half], pml[:z][:b_half], pml[:z][:k_halfI])
    memoryx!(mw[:dvzdx], w1t[:dvzdx], pml[:x][:a_half], pml[:x][:b_half], pml[:x][:k_halfI])
    memoryy!(mw[:dvzdy], w1t[:dvzdy], pml[:y][:a_half], pml[:y][:b_half], pml[:y][:k_halfI])

    @parallel compute_tauii!(
        w1t[:tauxx],
        w1t[:tauyy],
        w1t[:tauzz],
        w1t[:dvxdx],
        w1t[:dvydy],
        w1t[:dvzdz],
        pac.fc[:dt],
        pac.mod[:M],
        pac.mod[:lambda],
    )
    @parallel compute_tauij!(
        w1t[:tauxy],
        w1t[:tauxz],
        w1t[:tauyz],
        w1t[:dvxdy],
        w1t[:dvxdz],
        w1t[:dvydx],
        w1t[:dvydz],
        w1t[:dvzdx],
        w1t[:dvzdy],
        pac.fc[:dt],
        pac.mod[:mu],
    )

    if (pac.ic[:izfree] ≠ 0)
        @parallel (1:size(w1t[:tauzz], 2), 1:size(w1t[:tauzz], 3)) free_surface_mirror!(
            w1t[:tauzz],
            pac.ic[:izfree] + 1,
        )
        @parallel (1:size(w1t[:tauxz], 2), 1:size(w1t[:tauxz], 3)) free_surface!(
            w1t[:tauxz],
            pac.ic[:izfree],
        )
        @parallel (1:size(w1t[:tauyz], 2), 1:size(w1t[:tauyz], 3)) free_surface!(
            w1t[:tauyz],
            pac.ic[:izfree],
        )
    end

end

function advance_kernel!(pap, pac::T) where {T<:P_common{FdtdElastic,2}}
    w1t = pap.w1[:t]
    mw = pap.memory_pml
    pml = pac.pml

    @parallel compute_dtau!(
        w1t[:tauxx],
        w1t[:tauzz],
        w1t[:tauxz],
        w1t[:dtauxxdx],
        w1t[:dtauxzdx],
        w1t[:dtauzzdz],
        w1t[:dtauxzdz],
        pac.fc[:dxI],
        pac.fc[:dzI],
    )

    memoryx!(
        mw[:dtauxxdx],
        w1t[:dtauxxdx],
        pml[:x][:a_half],
        pml[:x][:b_half],
        pml[:x][:k_halfI],
    )
    memory1x!(mw[:dtauxzdx], w1t[:dtauxzdx], pml[:x][:a], pml[:x][:b], pml[:x][:kI])

    memoryz!(
        mw[:dtauzzdz],
        w1t[:dtauzzdz],
        pml[:z][:a_half],
        pml[:z][:b_half],
        pml[:z][:k_halfI],
    )
    memory1z!(mw[:dtauxzdz], w1t[:dtauxzdz], pml[:z][:a], pml[:z][:b], pml[:z][:kI])

    @parallel compute_v!(
        w1t[:vx],
        w1t[:vz],
        w1t[:dtauxxdx],
        w1t[:dtauxzdx],
        w1t[:dtauzzdz],
        w1t[:dtauxzdz],
        pac.fc[:dt],
        pac.mod[:rho],
    )

    @parallel (1:pac.ic[:nz]) dirichletx!(w1t[:vx], w1t[:vz], pac.ic[:nx])
    @parallel (1:pac.ic[:nx]) dirichletz!(w1t[:vz], w1t[:vx], pac.ic[:nz])
    @parallel compute_dv!(
        w1t[:vx],
        w1t[:vz],
        w1t[:dvxdx],
        w1t[:dvzdz],
        w1t[:dvxdz],
        w1t[:dvzdx],
        pac.fc[:dxI],
        pac.fc[:dzI],
    )
    memoryx!(mw[:dvxdx], w1t[:dvxdx], pml[:x][:a], pml[:x][:b], pml[:x][:kI])
    memoryz!(mw[:dvzdz], w1t[:dvzdz], pml[:z][:a], pml[:z][:b], pml[:z][:kI])

    memoryz!(mw[:dvxdz], w1t[:dvxdz], pml[:z][:a_half], pml[:z][:b_half], pml[:z][:k_halfI])

    memoryx!(mw[:dvzdx], w1t[:dvzdx], pml[:x][:a_half], pml[:x][:b_half], pml[:x][:k_halfI])

    @parallel compute_tauii!(
        w1t[:tauxx],
        w1t[:tauzz],
        w1t[:dvxdx],
        w1t[:dvzdz],
        pac.fc[:dt],
        pac.mod[:M],
        pac.mod[:lambda],
    )
    @parallel compute_tauij!(
        w1t[:tauxz],
        w1t[:dvxdz],
        w1t[:dvzdx],
        pac.fc[:dt],
        pac.mod[:mu],
    )

    if (pac.ic[:izfree] ≠ 0)
        @parallel (1:size(w1t[:tauzz], 2)) free_surface_mirror!(
            w1t[:tauzz],
            pac.ic[:izfree] + 1,
        )
        @parallel (1:size(w1t[:tauxz], 2)) free_surface!(w1t[:tauxz], pac.ic[:izfree])
    end

end


text:
function DotTest(m_rand,d_rand,operators,parameters)
	# Dot product test for a vector of linear operators

	m_adj = LinearOperator(d_rand,operators,parameters,adj=true)
	d_fwd = LinearOperator(m_rand,operators,parameters,adj=false)
	inner1 = InnerProduct(d_rand,d_fwd)
	inner2 = InnerProduct(m_rand,m_adj)
	println("<d_rand,d_fwd> = ",inner1)
	println("<m_rand,m_adj> = ",inner2)
	return inner1,inner2

end

function DotTest(m_rand::AbstractString,d_rand::AbstractString,operators,parameters)
	# Dot product test for a vector of linear operators

	rand_string = string(round(Int,rand()*100000))
	m_adj = join(["tmp_DotTest_m_adj_",rand_string])
	d_fwd = join(["tmp_DotTest_d_adj_",rand_string])
	LinearOperator(m_adj,d_rand,operators,parameters,adj=true)
	LinearOperator(m_rand,d_fwd,operators,parameters,adj=false)
	inner1 = InnerProduct(d_rand,d_fwd)
	inner2 = InnerProduct(m_rand,m_adj)
	println("These values should be close to each other:")
	println("<d_rand,d_fwd> = ",inner1)
	println("<m_rand,m_adj> = ",inner2)
	SeisRemove(m_adj)
	SeisRemove(d_fwd)

end

function DotTest(m_rand::Array{AbstractString,1},d_rand::Array{AbstractString,1},operators,parameters)
	# Dot product test for a vector of linear operators
	# for 3 component data

	rand_string = string(round(Int,rand()*100000))
	m1_adj = join(["tmp_DotTest_m1_adj_",rand_string])
	m2_adj = join(["tmp_DotTest_m2_adj_",rand_string])
	d1_fwd = join(["tmp_DotTest_d1_adj_",rand_string])
	d2_fwd = join(["tmp_DotTest_d2_adj_",rand_string])
	m_adj = [m1_adj;m2_adj]
	d_fwd = [d1_fwd;d2_fwd]
	LinearOperator(m_adj,d_rand,operators,parameters,adj=true)
	LinearOperator(m_rand,d_fwd,operators,parameters,adj=false)
	inner1 = InnerProduct(d_rand,d_fwd)
	inner2 = InnerProduct(m_rand,m_adj)
	println("These values should be close to each other:")
	println("<d_rand,d_fwd> = ",inner1)
	println("<m_rand,m_adj> = ",inner2)
	SeisRemove(m_adj);
	SeisRemove(d_fwd);

end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
#defining functions
function tria6(l,xs,ys)
    (l1,l2,l3) = l
    N1 = (2*l1-1)*l1; N2 = (2*l2-1)*l2; N3 = (2*l3-1)*l3;
    N4 = 4*l1*l2; N5 = 4*l2*l3; N6 = 4*l3*l1;
    x = dot(xs,[N1 N2 N3 N4 N5 N6])
    y = dot(ys,[N1 N2 N3 N4 N5 N6])
    return x,y
end
function BC(f, xedge, yedge, map)
    #f is vertical array
    #(fx,fy)
    if map[1] == 0
        l1 = zeros(7, 1)
        l2 = [-0.949108 -0.741531 -0.405845 0.405845 0.741531 0.949108 0.0]
        w = [0.129485 0.279705 0.38183 0.38183 0.279705 0.129485 0.417959]
        #transform to 0~1
        l2 = 0.5 .* l2 .+ 0.5
        l3 = 1 .- l2
        w = 0.5 .* w
        road_coord = 3
        map = 1
    elseif map[2] == 0
        l2 = zeros(7, 1)
        l3 = [-0.949108 -0.741531 -0.405845 0.405845 0.741531 0.949108 0.0]
        w = [0.129485 0.279705 0.38183 0.38183 0.279705 0.129485 0.417959]
        #transform to 0~1
        l3 = 0.5 .* l3 .+ 0.5
        l1 = 1 .- l3
        w = 0.5 .* w
        road_coord =3
        map = 2
    elseif map[3] == 0
        l3 = zeros(7, 1)
        l1 = [-0.949108 -0.741531 -0.405845 0.405845 0.741531 0.949108 0.0]
        w = [0.129485 0.279705 0.38183 0.38183 0.279705 0.129485 0.417959]
        #transform to 0~1
        l1 = 0.5 .* l1 .+ 0.5
        l2 = 1 .- l1
        w = 0.5 .* w
        road_coord  = 2
        map = 3
    else
        println("this is not a valid edge map")
        exit()
    end

    bc = zeros(12, 1)
    for i = 1:7
        l = (l1[i], l2[i], l3[i])
        (x, y) = tria6(l, xedge, yedge)

        road = NL1L2L3(l,xedge,yedge) #gives [xl1 yl1; xl2 yl2; xl3 yl3]
        if road_coord == 1
            road = road[1,:]
        elseif road_coord ==2
            road = road[2,:]
        else
            road = road[3,:]
        end

        bc = bc + w[i]*N(l)' * f([x, y])*sqrt(road[1]^2+road[2]^2)
    end
    return bc
end

#SHITS ABOUT TO GET REAL
function N(l)
    (l1, l2, l3) = l
    N1 = (2 * l1 - 1) * l1
    N2 = (2 * l2 - 1) * l2
    N3 = (2 * l3 - 1) * l3
    N4 = 4 * l1 * l2
    N5 = 4 * l2 * l3
    N6 = 4 * l3 * l1

    N = [N1 N2 N3 N4 N5 N6 0 0 0 0 0 0;
         0 0 0 0 0 0 N1 N2 N3 N4 N5 N6]
    return N
end
function NL1L2L3(l,xs,ys)
    J = JacobianL1L2L3(l,xs,ys)
    J = [1 1 1; J']
    P = J\[0 0;1 0; 0 1]
    nxy = Nxy(l,xs,ys)'

    nl1l2l3 = P'\nxy'
    xyl1l2l3 = nl1l2l3*[xs' ys']
    return xyl1l2l3
end
function Nl(l)
#    to be removed

    (l1, l2, l3) = l
    Nl1 = [(4 * l1 - 1) 0 0 (4 * l2) 0 (4 * l3)]
    Nl2 = [0 (4 * l2 - 1) 0 (4 * l1) (4 * l3) 0]
    Nl3 = [0 0 (4 * l3 - 1) 0 (4 * l2) (4 * l1)]
    Nl = [Nl1; Nl2; Nl3]
    #returns partials of Ni at l
    return Nl
end
function JacobianL1L2L3(l,xs,ys)

    J = Nl(l)*[xs' ys']
    return J
end
#=
function Nxyl1l2l3(l,xs,ys)
    #row 1 has x derivatives
    #row 2 has y derivatives
    nl = Nl(l)
    J = JacobianL1L2L3(l,xs,ys)
    nxy = J\nl
    return nxy
end
=#
function Nxy(l,xs,ys)
    #row 1 has x derivatives
    #row 2 has y derivatives
    Bo = Nl1l2(l)
    J = JacobianL1L2(l,xs,ys)
    nxy = J\Bo
    return nxy
end
function LN(l,xs,ys)
    nxy = Nxy(l,xs,ys)
    nx = nxy[1,:]'
    ny = nxy[2,:]'
    nx = [nx zeros(1,6);
          zeros(1,6) nx]
    ny = [ny zeros(1,6);
          zeros(1,6) ny]

    LN = [1 0; 0 0; 0 1]*nx + [0 0; 0 1; 1 0]*ny

    return LN
end

function Ematrix(e, v, cond)
    if cond == "stress"
        E = (e / (1 - v * v)) * [
            1 v 0
            v 1 0
            0 0 (1 - v) / 2]
    elseif cond == "strain"
        E =
            (e * (1 - v) / ((1 + v) * (1 - 2 * v))) * [
                1 v / (v - 1) 0
                v / (v - 1) 1 0
                0 0 (1 - 2 * v) / (2 - 2 * v)]
    else
        println("no planestress or planestrain condition")
        exit()
    end
    return E
end
function Nl1l2(l)
    (l1,l2,l3) = l
    nl1 = [(4*l1-1) 0 (-4*(1-l1-l2)+1) 4*l2 -4*l2 4*(1-l2-2*l1)]
    nl2 = [0 (4*l2-1) (-4*(1-l1-l2)+1) 4*l1 4*(1-l1-2*l2) -4*l1]
    return [nl1 ; nl2]
end
function JacobianL1L2(l,xs,ys)
    (l1,l2,l3) = l
    nl1 = [(4*l1-1) 0 (-4*(1-l1-l2)+1) 4*l2 -4*l2 4*(1-l2-2*l1)]
    nl2 = [0 (4*l2-1) (-4*(1-l1-l2)+1) 4*l1 4*(1-l1-2*l2) -4*l1]
    J = [nl1 ; nl2]*[xs' ys']
    return J
end
function ke(E,xs,ys)
    a1 = 0.0597158717
    b1 = 0.4701420641
    a2 = 0.7974269853
    b2 = 0.1012865073
    w1 = 0.2250000000
    w2 = 0.1323941527
    w3 = 0.1323941527
    w4 = 0.1323941527
    w5 = 0.1259391805
    w6 = 0.1259391805
    w7 = 0.1259391805
    g1 = [1 1 1] ./ 3
    g2 = [a1 b1 b1]
    g3 = [b1 a1 b1]
    g4 = [b1 b1 a1]
    g5 = [a2 b2 b2]
    g6 = [b2 a2 b2]
    g7 = [b2 b2 a2]
    w = [w1 w2 w3 w4 w5 w6 w7]
    g = [[g1] [g2] [g3] [g4] [g5] [g6] [g7]]

    ke = zeros(12, 12)
    for i = 1:7
        li = g[i]
        wi = w[i]
        jl1l2 = JacobianL1L2(li,xs,ys)
        ln = LN(li,xs,ys)
        ke = ke + 0.5*wi*ln'*E*ln*abs(det(jl1l2))
    end
    return ke
end

function me(xs,ys)
    a1 = 0.0597158717
    b1 = 0.4701420641
    a2 = 0.7974269853
    b2 = 0.1012865073
    w1 = 0.2250000000
    w2 = 0.1323941527
    w3 = 0.1323941527
    w4 = 0.1323941527
    w5 = 0.1259391805
    w6 = 0.1259391805
    w7 = 0.1259391805
    g1 = [1 1 1] ./ 3
    g2 = [a1 b1 b1]
    g3 = [b1 a1 b1]
    g4 = [b1 b1 a1]
    g5 = [a2 b2 b2]
    g6 = [b2 a2 b2]
    g7 = [b2 b2 a2]
    w = [w1 w2 w3 w4 w5 w6 w7]
    g = [[g1] [g2] [g3] [g4] [g5] [g6] [g7]]

    me = zeros(12, 12)
    for i = 1:7
        li = g[i]
        wi = w[i]
        jl1l2 = JacobianL1L2(li,xs,ys)
        ni = N(li)
        me = me + 0.5*wi*ni'*ni*abs(det(jl1l2))
    end
    return me
end
function IDs(bc)
    #returns
    #nbc: number of boundary conditions
    #ix: ids of edges for each condition
    #ni: contains number of edges for each condition
    nbc = count(x->x==-1,bc)+1
    ix = findall(x->x==-1,bc)
    ix = [ix[i][1] for i=1:nbc-1]
    prepend!(ix,0)
    ix = [collect(ix[i-1]+1:ix[i]-1) for i=2:nbc]

    ni = [length(i) for i in ix]
    ix = [ [i for i in bc[ix[j]]] for j=1:nbc-1]
    ix = [edge[ix[i],:] for i=1:nbc-1]
    return Int64(nbc-1),ni,ix
end
#=
xs = [1 0 0 0.5 0 0.5]
ys = [0 -.5 .5 -.25 0 .25]
E = Ematrix(220e9,0.5,"stress")
figure(1)
for l1=0:0.1:1
    for l2=0:0.1:1
        l3 = 1 -l1-l2
        if l3<0 break end
        l = [l1 l2 l3]
        (x,y) = tria6(l,xs,ys)
        f  = sum(N(l),dims=2)
        scatter3D(x,y,f[1])
    end
end
=#
xs = [1 0 0 0.5 0 0.5]
ys = [0 -.5 .5 -.25 0 .25]

max_stars_repo_name: billyon/first

max_issues_repo_name: billyon/first

max_forks_repo_name: billyon/first


text:
@testset "get_quantified_variables_form" begin
    @customtest String(get_quantified_variables_form(Formula("*{foo}P(foo(foo), foo, bar) & E(baz, bar)"))) ==
        "?{baz}?{bar}(*{foo}P(foo(foo), foo, bar) & E(baz, bar))"
    @customtest String(get_quantified_variables_form(Formula("*{x}?{y}*{z}?{w}(~P(a, w) | Q(f(x), y))"))) ==
        "?{a}*{x}?{y}*{z}?{w}(~P(a, w) | Q(f(x), y))"
end


text:
module MultistartOptimization

export MinimizationProblem, LocationValue, NLoptLocalMethod, local_minimization, TikTak,
    multistart_minimization

using ArgCheck: @argcheck
using Base.Threads: @spawn, fetch
using DocStringExtensions: FIELDS, SIGNATURES, TYPEDEF
using NLopt: NLopt
using Parameters: @unpack
using Sobol: SobolSeq, Sobol
using ProgressMeter

####
#### structures for the problem and results
####

"""
$(TYPEDEF)

Wrapper for a minimization problem.

# Fields

$(FIELDS)
"""
struct MinimizationProblem{F,T}
    "The function to be minimized."
    objective::F
    "Lower bounds (a vector of real numbers)."
    lower_bounds::T
    "Upper bounds (a vector of real numbers)."
    upper_bounds::T
    # FIXME constructor checks
end

"""
$(TYPEDEF)

A location-value pair.

# Fields

$(FIELDS)
"""
struct LocationValue{T <: AbstractVector{<:Real}, S <: Real}
    "Location (a vector of real numbers)."
    location::T
    "The value of the objective at `location`."
    value::S
end

####
#### internal utilities
####

"""
$(SIGNATURES)

Evaluate and return points of an `N`-element Sobol sequence.

An effort is made to parallelize the code using `Threads` when available.
"""
function sobol_starting_points(minimization_problem::MinimizationProblem, N::Integer; parallel_map = true)
    @unpack objective, lower_bounds, upper_bounds = minimization_problem
    s = SobolSeq(lower_bounds, upper_bounds)
    skip(s, N)                  # better uniformity
    if parallel_map
        map(fetch, map(x -> @spawn(LocationValue(x, objective(x))), Iterators.take(s, N)))
    else
        map(x -> LocationValue(x, objective(x)), Iterators.take(s, N))
    end
end

"""
$(SIGNATURES)

Helper function to keep the `N` points with the lowest `value`.
"""
function _keep_lowest(xs, N)
    @argcheck 1 ≤ N ≤ length(xs)
    partialsort(xs, 1:N, by = p -> p.value)
end

####
#### local minimization
####

Base.@kwdef struct NLoptLocalMethod
    algorithm::NLopt.Algorithm
    xtol_abs::Float64 = 1e-8
    xtol_rel::Float64 = 1e-8
    maxeval::Int = 0
    maxtime::Float64 = 0.0
end

"""
$(SIGNATURES)

A wrapper for algorithms supported by `NLopt`. Used to construct the corresponding
optimization problem.

See the NLopt documentation for the options. Defaults are changed slightly.
"""
function NLoptLocalMethod(algorithm::NLopt.Algorithm; options...)
    NLoptLocalMethod(; algorithm = algorithm, options...)
end

"""
$(SIGNATURES)

Solve `minimization_problem` using `local_method`, starting from `x`. Return a
`LocationValue`.
"""
function local_minimization(local_method::NLoptLocalMethod,
                            minimization_problem::MinimizationProblem, x)
    @unpack algorithm, xtol_abs, xtol_rel, maxeval, maxtime = local_method
    @unpack objective, lower_bounds, upper_bounds = minimization_problem
    opt = NLopt.Opt(algorithm, length(x))
    opt.lower_bounds = lower_bounds
    opt.upper_bounds = upper_bounds
    function f̃(x, grad)         # wrapper for NLopt
        @argcheck isempty(grad) # ensure no derivatives are asked for
        objective(x)
    end
    opt.min_objective = f̃
    opt.xtol_abs = xtol_abs
    opt.xtol_rel = xtol_rel
    opt.maxeval = maxeval
    opt.maxtime = maxtime
    optf, optx, ret = NLopt.optimize(opt, x)
    # FIXME ret is ignored, save it? is it useful?
    LocationValue(optx, optf)
end

####
#### multistart minimization
####

struct TikTak
    quasirandom_N::Int
    initial_N::Int
    θ_min::Float64
    θ_max::Float64
    θ_pow::Float64
end

"""
$(SIGNATURES)

The “TikTak” multistart method, as described in *Arnoud, Guvenen, and Kleineberg (2019)*.

This implements the *multistart* part, can be called with arbitrary local methods, see
[`multistart_minimization`](@ref).

# Arguments

- `quasirandom_N`: the number of quasirandom points for the first pass (using a Sobol
  sequence).

# Keyword arguments

- `keep_ratio`: the fraction of best quasirandom points which are kept

- `θ_min` and `θ_max` clamp the weight parameter, `θ_pow` determines the power it is raised
  to.

The defaults are from the paper cited above.
"""
function TikTak(quasirandom_N; keep_ratio = 0.1, θ_min = 0.1, θ_max = 0.995, θ_pow = 0.5)
    @argcheck 0 < keep_ratio ≤ 1
    TikTak(quasirandom_N, ceil(keep_ratio * quasirandom_N), θ_min, θ_max, θ_pow)
end

function _weight_parameter(t::TikTak, i)
    @unpack initial_N, θ_min, θ_max, θ_pow = t
    clamp((i / initial_N)^θ_pow, θ_min, θ_max)
end

"""
$(SIGNATURES)

Solve `minimization_problem` by using `local_method` within `multistart_method`.
"""
function multistart_minimization(multistart_method::TikTak, local_method,
                                 minimization_problem; progress=false, parallel_sobol_points=true)
    @unpack quasirandom_N, initial_N, θ_min, θ_max, θ_pow = multistart_method
    quasirandom_points = sobol_starting_points(minimization_problem, quasirandom_N; parallel_map=parallel_sobol_points)
    initial_points = _keep_lowest(quasirandom_points, initial_N)
    if progress
        prog = Progress(length(initial_points))
    end

    function _step(visited_minimum, (i, initial_point))
        θ = _weight_parameter(multistart_method, i)
        x = @. (1 - θ) * initial_point.location + θ * visited_minimum.location
        local_minimum = local_minimization(local_method, minimization_problem, x)
        if @isdefined(prog)
            ProgressMeter.next!(prog)
        end
        local_minimum.value < visited_minimum.value ? local_minimum : visited_minimum
    end
    foldl(_step, enumerate(initial_points); init = first(initial_points))
end

end # module


text:
module Sylvia

import AbstractInstances
import Cassette
import Combinatorics: permutations
import DataStructures: OrderedDict
import LinearAlgebra
import MacroTools: striplines, unblock

export @S_str, @λ, @scope, @!, @sym,
    Mock, Sym, Wild,
    gather, substitute, tagof

##################################################
# Misc
##################################################

@inline istrue(x) = x == true
@inline isfalse(x) = x == false
@inline commuteswith(op, x, y) = op(x, y) == op(y, x)

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

include("cassettectx.jl")

include("sym.jl")
include("wild.jl")
include("match.jl")
include("promotion.jl")
include("expr.jl")
include("show.jl")
include("compile.jl")
include("context.jl")
include("register.jl")
include("apply.jl")
include("mock.jl")
include("protect.jl")

include("substitute.jl")
include("simplify.jl")
include("array.jl")
include("ops.jl")
include("rules.jl")

end # module


text:
using PoincareInvariants

export guiding_center_4d_loop_ode,
       guiding_center_4d_loop_iode,
       guiding_center_4d_loop_vode

export guiding_center_4d_ode_poincare_invariant_1st,
       guiding_center_4d_iode_poincare_invariant_1st,
       guiding_center_4d_vode_poincare_invariant_1st

function f_loop(i, n)
   f_loop(i/n)
end

function initial_conditions_loop(n)
   q₀ = zeros(4, n)

   for i in axes(q₀,2)
       q₀[:,i] .= f_loop(i, n)
   end

   return q₀
end

guiding_center_4d_loop_ode_init(q₀) = guiding_center_4d_ode(q₀, (μ = μ_loop(),); periodic=false)
guiding_center_4d_loop_iode_init(q₀) = guiding_center_4d_iode(q₀, (μ = μ_loop(),); periodic=false)
guiding_center_4d_loop_vode_init(q₀) = guiding_center_4d_vode_formal_lagrangian(q₀, (μ = μ_loop(),); periodic=false)

function guiding_center_4d_loop_ode(n)
   guiding_center_4d_loop_ode_init(initial_conditions_loop(n))
end

function guiding_center_4d_loop_iode(n)
   guiding_center_4d_loop_iode_init(initial_conditions_loop(n))
end

function guiding_center_4d_loop_vode(n)
   guiding_center_4d_loop_iode_init(initial_conditions_loop(n))
end

function guiding_center_4d_ode_poincare_invariant_1st(Δt, nloop, ntime, nsave, DT=Float64)
   PoincareInvariant1st(guiding_center_4d_loop_ode_init, f_loop, ϑ, Δt, 4, nloop, ntime, nsave, DT)
end

function guiding_center_4d_iode_poincare_invariant_1st(Δt, nloop, ntime, nsave, DT=Float64)
   PoincareInvariant1st(guiding_center_4d_loop_iode_init, f_loop, ϑ, Δt, 4, nloop, ntime, nsave, DT)
end

function guiding_center_4d_vode_poincare_invariant_1st(Δt, nloop, ntime, nsave, DT=Float64)
   PoincareInvariant1st(guiding_center_4d_loop_vode_init, f_loop, ϑ, Δt, 4, nloop, ntime, nsave, DT)
end


text:
mutable struct Parameters{T}
    diam::T
    rmax::T
    ρ::T
    mr::Integer
    nρ::Integer
    δρ::T
    ktemp::T
end

function Parameters(diam, rmax, ρ, mr::Integer, nr, ktemp)
    (d, rm, r, n) = promote(diam, rmax, ρ, nr)

    return Parameters{typeof(d)}(d, rm, r, mr, n, zero(typeof(d)), ktemp)
end

struct Structure{T}
    r::T
    q::T
    pot::Potential
end

function Structure(p::Parameters, pot::Potential)
    dr = p.rmax / p.mr
    dq = π / p.rmax
    r = range(0.0, (p.mr - 1) * dr; length=p.mr)
    q = range(0.0, (p.mr - 1) * dq; length=p.mr)

    return Structure{typeof(r)}(r, q, pot)
end

struct Result{T<:AbstractArray}
    r::T
    q::T
    gr::T
    cr::T
    sq::T
end

function Result()
    vecs = [zeros(2) for _ in 1:5]
    return Result{Vector{Float64}}(vecs...)
end

function Result(cr, gam, st::Structure, params::Parameters)
    # First, we deal with the RDF
    gr = @. gam[2:end] + cr[2:end] + 1.0
    result_cr = copy(cr)

    # Next, we deal with the structure factor
    fits = continuous!(cr, st.r, params.diam)

    r, p = make_fft_plan(cr)
    ck = fft_oz(cr, r, params.rmax, params.mr, p)
    remove_continuous!(ck, st.q, params.diam, fits)

    sq_inv = @. 1.0 - params.ρ * ck[2:end]

    # Also save the new structure information
    new_r = collect(st.r[2:end])
    new_q = collect(st.q[2:end])

    return Result{typeof(gr)}(new_r, new_q, gr, result_cr, 1.0 ./ sq_inv)
end

struct Interaction
    params::Parameters
    structure::Structure
    bridge::Closure
end


text:
function sqrt_only_positive(;x)
    return (;sqrt=sqrt(x))
end

fns = (;sqrt_only_positive, ) # implicit field name: see https://docs.julialang.org/en/v1/base/base/#Core.NamedTuple


text:
module Sand_03
using LinearAlgebra, Random

#export ring, tstep!, step!

#using StaticArrays
"""
    Sand 03: Hyperbolic diffusion, with no advection (U=0)
"""

JIGGLE = 0.9 # how much to randomly tweak the node positions, to test algo; 0 ≤ JIGGLE < 1.0
q = 1.0 # hyperbolic wave speed
s = 3.0 # relaxation time
r₊ = [1.0, √q]
r₋ = [1.0,-√q]

mutable struct ring
    N::Integer              # number of cells
    Δt::Float64             # time step
    X_::Array{Float64,1}    # positions (nodes)
    ΔX_::Array{Float64,1}   # cell sizes (cell)
#    U_::Array{Float64,1}    # velocity (on nodes) XXX IGNORED
#    ∇U_::Array{Float64,1}   # divergence of velocity (on cells) XXX IGNORED
    Σ_::Array{Float64,1}    # surface density
    S_::Array{Float64,1}    # integral of the above
    J_::Array{Float64,1}    # Flux
    F_::Array{Float64,1}    # integral of flux
    α⁺::Array{Float64,1}    #
    α⁻::Array{Float64,1}    #
    t::Float64              # time
end
function ring(N=100; U = 1.0, xi=0.0, xf=10.0)
    xtemp_ = collect(range(xi,xf,length=N+1))
    xtemp_[2:end-1] += (Random.rand(N-1).-0.5) * (xf-xi)/N * JIGGLE
    @assert prod(diff(xtemp_).>0)  # checks to make sure xtemp_ monotonic
    X_ = xtemp_#[2:end]
    ΔX_ = diff(xtemp_)

    # The next several lines are superfluous since
#    U_ = U * ones(N+1)
#    """ XXX putting step-down in speed XXX """
#    U_[3*N÷4:end] .*= 0.5
#    U_[1:N÷4]     .*= 0.5
#    ∇U_ = diff(U_) ./ ΔX_
    # init w/top-hat density:
    Σ_ = ones(N); Σ_[1] = 0.0;
    Σ_[1:(2*N÷5)] .= 0.0; Σ_[(3*N÷5):end] .= 0.0
    S_ = Σ_ .* ΔX_
    J_ = zeros(N); F_ = J_ * 1.0
    #
    #Δt = 0.5 / max(maximum(U_[1:end-1]./ΔX_),-minimum(U_[2:end]./ΔX_))
    Δt = 0.5 / max( maximum(q./ΔX_), -minimum(-q./ΔX_) , s)
    #
    return ring(N,Δt,X_,ΔX_,Σ_,S_,J_,F_,zeros(N+1),zeros(N+1),0.0)
end

function getAlphas!(r::ring)
    N = r.N; Σ_ = r.Σ_; J_ = r.J_; α⁺ = r.α⁺; α⁻ = r.α⁻
    ΔΣ_ = cat(Σ_[1]-Σ_[end], diff(Σ_), Σ_[1]-Σ_[end], dims=1)
    ΔJ_ = cat(J_[1]-J_[end], diff(J_), J_[1]-J_[end], dims=1)
    for i in 1:(N+1) # Tried to write this as two one-liners instead of loop, but probs w/converting Vec{Array} to Array...
        α⁺[i] = 0.5 * [1  1/√q] ⋅ [ΔΣ_[i], ΔJ_[i]]
        α⁻[i] = 0.5 * [1 -1/√q] ⋅ [ΔΣ_[i], ΔJ_[i]]
    end
    r.α⁺ = α⁺
    r.α⁻ = α⁻
end#function

function step!(r::ring; Δt=r.Δt)
    Σ_=r.Σ_; N = r.N; ΔX_ = r.ΔX_; S_ = r.S_; J_=r.J_; F_ =r.F_
    #
    getAlphas!(r)
    #
    # get integral quantities
    S_ = Σ_ .* ΔX_
    F_ = J_ .* ΔX_
#    U⁺_ = U_[1:end-1] .* (U_[1:end-1].>0) # right-traveling
#    U⁻_ = U_[2:end  ] .* (U_[2:end  ].<0) # left-traveling
    # right-traveling. nb: √q is wave speed
    S_ -= r.α⁺[1:end-1] * Δt * r₊[1] * √q
    F_ -= r.α⁺[1:end-1] * Δt * r₊[2] * √q
    # left-traveling
    S_ += r.α⁻[2:end] * Δt * r₋[1] * √q
    F_ += r.α⁻[2:end] * Δt * r₋[2] * √q

    # This part kinda sucks b/c it doesn't conserve very well:
    Σ_ = S_ ./ ΔX_;  J_ = F_ ./ ΔX_
    #Σ_ += -Δt * r.∇U_ .* Σ_ # <- forward Euler
    #  Σ_ ./= (1 .+ Δt*r.∇U_)  # <- backward Euler  - conserves a bit better
    J_ ./= (1 + Δt / s)
    #
    S_ = Σ_ .* ΔX_; F_ = J_ .* ΔX_

    #
    r.S_ = S_
    r.Σ_ = S_ ./ ΔX_
    r.F_ = F_
    r.J_ = F_ ./ ΔX_

    r.t += Δt
    s1 = sum(S_); s2 = sum(F_)
    #print("Sums: $s1,  $s2 \n")
    return (sum(S_),sum(F_))
end

function tstep!(r::ring,Δt)
    M = floor(Δt / r.Δt)
    Δt_last = Δt - M*r.Δt
    @assert Δt_last ≥ 0
    for i in 1:M
        step!(r)
    end
    step!(r, Δt = Δt_last)
end

end#module

max_stars_repo_name: diakosmos/razor

max_issues_repo_name: diakosmos/razor

max_forks_repo_name: diakosmos/razor


text:
"""
Internal implementation detail.
"""
module Common

using Compat

include("bytes.jl")
include("errors.jl")

end


text:
using HTTP
using SHA

function lookup_pwned_api(pwd::AbstractString)
    # Compute the SHA1 hash of pwd and split it at the fifth position
    sha1pwd = bytes2hex(sha1(pwd))
    head,tail = sha1pwd[1:5], sha1pwd[6:end]

    # Query the pwnedpasswords API
    url = "https://api.pwnedpasswords.com/range/$head"
    res = HTTP.request("GET",url)
    if (res.status != 200)
        @error "The API lookup failed!"
    end

    # The hashes are stored in the body of the response as hex characters (notice: NOT integers)
    # To split them, look for line-feeds, which in this case are given by CRLF (boo)
    hashes = split(String(res.body),"\r\n")
    hashes = [split(st,':') for st in hashes]

    # Find whether or not pwd (or rather its hash) appears in the list of all hashes!
    # Because the characters in the response are uppercase (01233456789ABCDEF)
    # and bytes2hex produces lowercase (0123456789abcdef), comparing strings is not good enough.
    # Instead do a RegEx match with the i flag set!
    found = 0
    re = Regex(tail,"i")
    for h in hashes
        # h[1] is the tail of the hash, h[2] is the number of occurances
        if match(re,h[1]) != nothing
            found = h[2]
            break
        end
    end
    return sha1pwd,found
end

function check_pwd(pwd::AbstractString)
    # Check whether (and if so, how often) the string pwd has been pwned.
    sha1pwd,cnt = "",0
    try
        sha1pwd,cnt = lookup_pwned_api(pwd)
    catch e
        # Error Message is printed elsewhere. Just return silently
        return
    end
    # cnt is the number of occurances. If pwd has not been pwned yet, cnt is zero.
    print("Password $pwd with the hash $sha1pwd was ")
    if (cnt != 0)
        print("found $cnt times!\n")
    else
        print("not found!\n")
    end
end

# ARGS stores the command-line arguments passed to pwned.jl
# If none were given, ARGS is empty and pwned.jl will read from stdin
if ARGS == []
    print("Reading passwords from stdin. Type exit or ^C to quit.\n")
    while true
        pwd = readline()
        if pwd == "exit"
            break
        end
        check_pwd(pwd)
    end
else
    @warn "Entering passwords in plain text on a terminal is a bad idea!i\nInstead run this program without any arguments,\nin which case it reads from stdin and no history will be created.\nIt is highly recommended to clear your terminal's history after using this program!"
    for pwd in ARGS
        check_pwd(pwd)
    end
end

max_stars_repo_name: ManAtTheMoon/pwned-search

max_issues_repo_name: ManAtTheMoon/pwned-search

max_forks_repo_name: ManAtTheMoon/pwned-search


text:
abstract type ShaderResource end

struct SampledImage <: ShaderResource
    image::Allocated{Created{Image,ImageCreateInfo},DeviceMemory}
    view::Allocated{Created{ImageView,ImageViewCreateInfo},DeviceMemory}
    sampler::Sampler
end

Vulkan.DescriptorType(resource::ShaderResource) = DescriptorType(typeof(resource))
Vulkan.DescriptorType(::Type{SampledImage}) = DESCRIPTOR_TYPE_COMBINED_IMAGE_SAMPLER

function Vulkan.WriteDescriptorSet(descriptor::Descriptor, resource::SampledImage)
    WriteDescriptorSet(
        descriptor.set,
        descriptor.binding,
        descriptor.index - 1, # 1 to 0-based indexing
        DescriptorType(resource),
        [DescriptorImageInfo(resource.sampler, resource.view, IMAGE_LAYOUT_SHADER_READ_ONLY_OPTIMAL)],
        [],
        [],
    )
end

struct StorageBuffer <: ShaderResource
    buffer::Allocated{Buffer,DeviceMemory}
end

Vulkan.DescriptorType(::Type{StorageBuffer}) = DESCRIPTOR_TYPE_STORAGE_BUFFER


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

include("../../../src/Sparso.jl")
using Sparso

function pagerank(A, p, r, d_inv, maxiter) # p: initial rank, r: damping factor
#  set_matrix_property(Dict(
#    :A => SA_SYMM_STRUCTURED | SA_SYMM_VALUED | SA_STRUCTURE_ONLY))

  bytes = maxiter*(nnz(A)*4 + size(A, 1)*4*8)

  p = p.*d_inv
  d_inv = copy(d_inv)
  Ap = zeros(size(A, 1))

  t = time()

  for i = 1:maxiter
    Ap = (1-r)*A*p + r
    Ap = Ap.*d_inv

    if i == maxiter
      err = norm(Ap - p)/norm(p)
      println("error = $err")
    end

    temp = Ap
    Ap = p
    p = temp
  end

  t = time() - t
  println("pagerank takes $t sec ($(bytes/t/1e9) gbps)")

  p
end

function pagerank_with_init(A0, r, d_inv, maxiter) # p: initial rank, r: damping factor
  set_matrix_property(Dict(
    :A0 => SA_SYMM_VALUED)
  )

  A0 = spones(A0)

  m = size(A0, 1)
  p = repmat([1/m], m)

  d = max(convert(Array{eltype(A0),1}, vec(sum(A0, 2))), 1) # num of neighbors
  A = scale(A0,1./d)

  bytes = maxiter*(nnz(A)*4 + m*4*8)

  p = p.*d_inv
  d_inv = copy(d_inv)
  Ap = zeros(size(A, 1))

  t = time()

  for i = 1:maxiter
    Ap = (1-r)*A*p + r
    Ap = Ap.*d_inv

    if i == maxiter
      err = norm(Ap - p)/norm(p)
      println("error = $err")
    end

    temp = Ap
    Ap = p
    p = temp
  end

  t = time() - t
  println("pagerank takes $t sec ($(bytes/t/1e9) gbps)")

  p
end

A0 = matrix_market_read(ARGS[1], true, true)
A0 = spones(A0)

m = size(A0, 1)
p = repmat([1/m], m)
r = 0.15

d = max(convert(Array{eltype(A0),1}, vec(sum(A0, 2))), 1) # num of neighbors
A = scale(A0,1./d)

maxiter = 100
bytes = maxiter*(nnz(A)*4 + m*4*8)

if length(ARGS) == 2
  test = ARGS[2]
else
  test = "julia"
end

d_inv = 1./d

if test == "call-repl"
  set_options(SA_ENABLE, SA_USE_SPMP, SA_REPLACE_CALLS)
elseif test == "context"
  set_options(SA_ENABLE, SA_USE_SPMP, SA_CONTEXT, SA_REPLACE_CALLS)
elseif test == "reorder"
  set_options(SA_ENABLE, SA_USE_SPMP, SA_CONTEXT, SA_REORDER, SA_REPLACE_CALLS)
elseif test == "verbose"
  set_options(SA_ENABLE, SA_USE_SPMP, SA_VERBOSE, SA_CONTEXT, SA_REORDER, SA_REPLACE_CALLS)
end

println("compiler warmup (ignored): ")
if test == "julia"
  x = pagerank(A0, p, r, d_inv, maxiter)
else
  @acc x= pagerank_with_init(A0, r, d_inv, maxiter)
end

println("\nRUN: ")
A2 = copy(A0)
p = repmat([1/m], m)
if test == "julia"
  x = pagerank(A2, p, r, d_inv, maxiter)
else
  #Sparso.reset_spmp_spmv_time()
  #Sparso.reset_knob_spmv_time()
  @acc x = pagerank_with_init(A2, r, d_inv, maxiter)
end

  #t = Sparso.get_spmp_spmv_time()
  #println("time spent on spmp spmv $t sec ($(bytes/t/1e9) gbps)")
  #t = Sparso.get_knob_spmv_time()
  #println("time spent on knob spmv $t sec ($(bytes/t/1e9) gbps)")

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
const NOTSOLVED = :NotSolved
const SOVLERERROR = :SolverError
const OPTIMAL = :Optimal
const UNBOUNDED = :Unbounded
const INFEASIBLE = :Infeasible
const UNBNDORINF = :UnboundedOrInfeasible
const NORESULTS = :ResultsNotRetrieved
const SOLVED = :Solved

struct NEOSSolverError <: Exception
    msg::String
end

abstract type NEOSModel <: MPB.AbstractMathProgModel end

addemail!(m::NEOSModel, email::String) = addemail!(m.solver, email)

function MPB.optimize!(m::NEOSModel)
    neos_writexmlmodel!(m)
    add_solver_xml!(m.solver, m)
    if m.solver.requires_email && m.solver.server.email==""
	throw(NEOSSolverError("$(typeof(m.solver)) requires that NEOS users supply a valid email"))
    end
    if m.solver.server.email != ""
	if match(r"<email>.*</email>", m.xmlmodel) == nothing
	    m.xmlmodel = replace(m.xmlmodel, r"<document>" => "<document>\n<emaixl>$(m.solver.server.email)</email>")
	else
	    m.xmlmodel = replace(m.xmlmodel, r"<email>.*</email>" => "<email>$(m.solver.server.email)</email>")
	end
    end
    # println(m.xmlmodel)
    starttime = time()
    job = submitJob(m.solver.server, m.xmlmodel, m.solver.print_level)
    #println("job $(job.number) $(job.password)")
    iterations = 0
    pollingperiod = 5.0
    offset = 0
    results_changed = true
    while true
        stat = getJobStatus(m.solver.server, job)
	if results_changed && m.solver.print_level >= 1
	    (intermediate_results, new_offset) = getIntermediateResultsNonBlocking(m.solver.server, job; offset=offset)
	    results_changed = new_offset != offset
	    offset = new_offset
	    println(intermediate_results)
	end
        # println("Waiting for results. Status: $(stat)")
        if stat == "Done"
	    if time() - starttime > 28_800 # NEOS 8hr limit
		@error("""You've reached the NEOS 8hr timelimit. No meaningful
		          results can be returned. You should set a timelimit that is less
			  than 8 hours using the solver parameters.""")
	    end
            break
        else
            sleep(pollingperiod)
        end
	iterations += 1
	if mod(iterations, 10) == 0
	    pollingperiod += 10.0
	end

    end
    m.last_results = getFinalResults(m.solver.server, job)
    if m.solver.print_level >= 2 || m.solver.print_results
	println(m.last_results)
    end
    if m.solver.result_file != ""
	open(m.solver.result_file, "w") do f
	    write(f, m.last_results)
	end
    end
    parseresults!(m, job)
    if m.status == NORESULTS
        @info("Try retrieving problem results from NEOS: https://neos-server.org/neos/jobs/$(10_000 * floor(Int, job.number / 10_000))/$(job.number).html")
    end
    return m.status
end


text:
using Test, AiEmpires

t1 = 1+2
@test t1 == 3


text:
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

const VOV = MOI.VectorOfVariables
const VI = MOI.VariableIndex

const OPTIMIZER_CONSTRUCTOR = MOI.OptimizerWithAttributes(
    CBLS.Optimizer, MOI.Silent() => true
)
const OPTIMIZER = MOI.instantiate(OPTIMIZER_CONSTRUCTOR)

@testset "LocalSearchSolvers" begin
    @test MOI.get(OPTIMIZER, MOI.SolverName()) == "LocalSearchSolvers"
end

# @testset "supports_default_copy_to" begin
#     @test MOIU.supports_default_copy_to(OPTIMIZER, false)
#     # Use `@test !...` if names are not supported
#     @test !MOIU.supports_default_copy_to(OPTIMIZER, true)
# end

const BRIDGED = MOI.instantiate(
    OPTIMIZER_CONSTRUCTOR, with_bridge_type = Float64
)
const CONFIG = MOIT.TestConfig(atol=1e-6, rtol=1e-6)

# @testset "Unit" begin
#     # Test all the functions included in dictionary `MOI.Test.unittests`,
#     # except functions "number_threads" and "solve_qcp_edge_cases."
#     MOIT.unittest(
#         BRIDGED,
#         CONFIG,
#         ["number_threads", "solve_qcp_edge_cases"]
#     )
# end

# @testset "Modification" begin
#     MOIT.modificationtest(BRIDGED, CONFIG)
# end

# @testset "Continuous Linear" begin
#     MOIT.contlineartest(BRIDGED, CONFIG)
# end

# @testset "Continuous Conic" begin
#     MOIT.contlineartest(BRIDGED, CONFIG)
# end

# @testset "Integer Conic" begin
#     MOIT.intconictest(BRIDGED, CONFIG)
# end
@testset "MOI: examples" begin
    # m = LocalSearchSolvers.Optimizer()
    # MOI.add_variables(m, 3)
    # MOI.add_constraint(m, VI(1), LS.DiscreteSet([1,2,3]))
    # MOI.add_constraint(m, VI(2), LS.DiscreteSet([1,2,3]))
    # MOI.add_constraint(m, VI(3), LS.DiscreteSet([1,2,3]))

    # MOI.add_constraint(m, VOV([VI(1),VI(2)]), LS.MOIPredicate(allunique))
    # MOI.add_constraint(m, VOV([VI(2),VI(3)]), LS.MOIAllDifferent(2))

    # MOI.set(m, MOI.ObjectiveFunction{LS.ScalarFunction}(), LS.ScalarFunction(sum, VI(1)))

    # MOI.optimize!(m)

    m1 = CBLS.Optimizer()
    MOI.add_variable(m1)
    MOI.add_constraint(m1, VI(1), CBLS.DiscreteSet([1,2,3]))

    m2 = CBLS.Optimizer()
    MOI.add_constrained_variable(m2, CBLS.DiscreteSet([1,2,3]))

    # opt = CBLS.sudoku(3, modeler = :MOI)
    # MOI.optimize!(opt)
    # @info solution(opt)
end