text:
import Random.shuffle
using Logging

export random_contraction_plan, inorder_contraction!
export contract_pair!, contract_network!, full_wavefunction_contraction!
export compress_tensor_chain!, decompose_tensor!, compress_bond!
export contract_mps_tensor_network_circuit!, contract_tensor_network_circuit_with_compression!
export calculate_mps_amplitudes!
export cost

include("layer2/slicing.jl")
include("layer2/netcon_contraction.jl")
include("layer2/bgreedy_contraction.jl")
include("layer2/quickbb_contraction.jl")
# include("layer2/KaHyPar_contraction.jl")

using HDF5

"""
    merge_common_bonds!(network::TensorNetworkCircuit, a_label::Symbol, b_label::Symbol)

Merge the common bonds between two nodes
"""
function merge_common_bonds!(network::TensorNetworkCircuit, a_label::Symbol, b_label::Symbol)
    a = network.nodes[a_label]
    b = network.nodes[b_label]
    common_edges = intersect(a.indices, b.indices)

    if length(common_edges) > 1
        index_map_a = Dict([v => k for (k, v) in enumerate(a.indices)])
        index_map_b = Dict([v => k for (k, v) in enumerate(b.indices)])

        a_common_indices = findall(x -> x in b.indices, a.indices)
        a_remaining_indices = findall(x -> !(x in b.indices), a.indices)
        permute_tensor(network, a_label, vcat(a_remaining_indices, a_common_indices))

        b_common_indices = findall(x -> x in a.indices, b.indices)
        b_remaining_indices = findall(x -> !(x in a.indices), b.indices)
        permute_tensor(network, b_label, vcat(b_remaining_indices, b_common_indices))

        new_edge = new_label!(network, "index")
        indices_a = vcat(a.indices[a_remaining_indices], [new_edge])
        indices_b = vcat(b.indices[b_remaining_indices], [new_edge])

        dims_map_a = Dict{Symbol, Int64}(a.indices .=> a.dims)
        dims_map_b = Dict{Symbol, Int64}(b.indices .=> b.dims)

        dims_map_a[new_edge] = dims_map_b[new_edge] = prod([dims_map_a[a.indices[ind]] for ind in a_common_indices])
        dims_a = [dims_map_a[index] for index in indices_a]
        dims_b = [dims_map_b[index] for index in indices_b]

        network.nodes[a_label] = Node(indices_a, dims_a, a_label)
        network.nodes[b_label] = Node(indices_b, dims_b, b_label)

        l, m = length(a_remaining_indices), length(a_common_indices)
        groups = [map(x -> [x], 1:l)..., collect((l+1):l+m)]
        reshape_tensor(network, a_label, groups)

        l, m = length(b_remaining_indices), length(b_common_indices)
        groups = [map(x -> [x], 1:l)..., collect((l+1):l+m)]
        reshape_tensor(network, b_label, groups)

        network.edges[new_edge] = Edge(a_label, b_label, nothing, true)
        for common_edge in common_edges
            delete!(network.edges, common_edge)
        end
    end
end

"""
    inorder_contraction!(network::TensorNetworkCircuit)

Function to contract the network in order starting from input nodes
"""
function inorder_contraction!(network::TensorNetworkCircuit)
    layer_nodes = Dict{Int, Array{Symbol,1}}()
    for (k, v) in pairs(network.node_layers)
        if haskey(layer_nodes, v)
            push!(layer_nodes[v], k)
        else
            layer_nodes[v] = [k]
        end
    end

    # create array of layers to contract in
    gate_layers = [k for k in sort(collect(keys(layer_nodes))) if k > 0]
    # the output layer has layer index -1 to disginguish it. If it exists, add it to list
    if haskey(layer_nodes, -1)
        gate_layers = vcat(gate_layers, [-1])
    end

    for layer in gate_layers
        nodes = layer_nodes[layer]

        new_nodes = Array{Symbol, 1}(undef, length(nodes))
        for (i, node) in enumerate(nodes)

            in_nodes = unique(inneighbours(network, node))
            for in_node in in_nodes
                node = contract_pair!(network, in_node, node)
            end
            new_nodes[i] = node
        end
        # merge multiple virtual bonds if present
        for (x, y) in Iterators.product(new_nodes, new_nodes)
            if y > x
                merge_common_bonds!(network, x, y)
            end
        end
    end
end

"""
    random_contraction_plan(network::TensorNetworkCircuit)

Function to create a random contraction plan
"""
function random_contraction_plan(network::TensorNetworkCircuit)
    closed_edges = [k for (k, v) in pairs(network.edges) if v.src !== nothing && v.dst !== nothing]
    shuffle(closed_edges)
end

"""
    full_wavefunction_contraction!(network::TensorNetworkCircuit,
                                            output_shape::Union{String, Array{Int, 1}}="")

Function to contract a network by first contracting input nodes together, to
get the wavefunction representing the initial state, and then contracting
gates into the wavefunction in the order they appear in the circuit.
"""
function full_wavefunction_contraction!(network::TensorNetworkCircuit,
                                        output_shape::Union{String, Array{Int, 1}}="")

    # Get the input nodes of the circuit
    input_nodes = [network.edges[edge].src for edge in network.input_qubits]
    if nothing in input_nodes
        error("Please create nodes for each input before using wavefunction
               contraction")
    end

    # Contract all of the input nodes to get a tensor representing the full
    # initial wavefunction.
    wf = input_nodes[1]
    for wfi in input_nodes[2:end]
        wf = contract_pair!(network, wf, wfi)
    end

    layer_nodes = Dict{Int, Array{Symbol,1}}()
    for (k, v) in pairs(network.node_layers)
        if haskey(layer_nodes, v)
            push!(layer_nodes[v], k)
        else
            layer_nodes[v] = [k]
        end
    end

    gate_layers = [k for k in sort(collect(keys(layer_nodes))) if k > 0]
    if haskey(layer_nodes, -1)
        gate_layers = vcat(gate_layers, [-1])
    end
    for layer in gate_layers
        nodes = layer_nodes[layer]

        # contract nodes in the layer and then contract with the wf vector
        n = nodes[1]
        if length(nodes) > 1
            for n2 in nodes[2:end]
                n = contract_pair!(network, n, n2)
            end
        end
        wf = contract_pair!(network, wf, n)
    end

    # Permute the indices of the final tensor to have the correct order.
    # output_tensor = new_label!(tn, "node")
    node = network.nodes[wf]
    if length(node.indices) != 0 # if not a scalar
        # use the qubit_ordering to change order of output qubits as required
        output_qubits = network.output_qubits[network.qubit_ordering]
        order = [findfirst(x->x==ind, node.indices) for ind in output_qubits]
        permute_tensor(network, wf, order)

        # Reshape the final tensor if a shape is specified by the user.
        if output_shape == "vector"
            reshape_tensor(network, wf, [collect(1:length(node.indices))])
        elseif output_shape != ""
            reshape_tensor(network, wf, output_shape)
        end
    end

    # save the final tensor under the name "result".
    save_output(network, wf)
    wf
end

# *************************************************************************** #
#                  Functions to contract a tensor network
# *************************************************************************** #

"""
    sort_indices(A::Node, B::Node)

This function divides all of the indices of two nodes A and B into two arrays,
one array for all the shared indices between A and B and a second array for all
the indices that are not shared.
"""
function sort_indices(A::Node, B::Node)
    common_indices = intersect(A.indices, B.indices)
    remaining_indices_A = setdiff(A.indices, common_indices)
    remaining_indices_B = setdiff(B.indices, common_indices)
    uncommon_indices = union(remaining_indices_A, remaining_indices_B)

    common_indices, uncommon_indices
end

"""
    create_ncon_indices(A::Node, B::Node,
                                 common_indices::Array{Symbol, 1},
                                 uncommon_indices::Array{Symbol, 1})

This function returns ncon indices for the nodes A and B. These can then to be
passed to ncon in order to contract the tensor data associated with A and B.
"""
function create_ncon_indices(A::Node, B::Node,
                             common_indices::Array{Symbol, 1},
                             uncommon_indices::Array{Symbol, 1})
    # Contracted indices are assigned positive integers, uncontracted indices
    # are assigned negative integers.
    common_index_map = Dict([(x[2], x[1]) for x in enumerate(common_indices)])
    uncommon_indices_map = Dict([(x[2], -x[1]) for x
                                    in enumerate(uncommon_indices)])

    index_map = merge(common_index_map, uncommon_indices_map)
    A_ncon_indices = [index_map[x] for x in A.indices]
    B_ncon_indices = [index_map[x] for x in B.indices]

    return A_ncon_indices, B_ncon_indices
end

"""
    contract_network!(network::TensorNetworkCircuit,
                               plan::Array{Symbol, 1},
                               output_shape::Union{String, Array{Int, 1}})

Function to contract the given network according to the given contraction plan.
The resulting tensor will be given the shape described by 'output_shape'.
"""
function contract_network!(network::TensorNetworkCircuit,
                           plan::Array{Symbol, 1},
                           output_shape::Union{String, Array{Int, 1}}="")

    # Loop through the plan and contract each edge in sequence
    for edge in plan
        contract_pair!(network, edge)
    end

    # Contract disjoint pieces of the network, if any
    while length(network.nodes) > 1
        n1, state = iterate(network.nodes)
        n2, _ = iterate(network.nodes, state)
        contract_pair!(network, n1.first, n2.first)
    end

    # Permute the indices of the final tensor to have the correct order.
    output_tensor = Symbol("node_$(network.counters["node"])")
    node = network.nodes[output_tensor]
    if length(node.indices) != 0
        order = [findfirst(x->x==ind, node.indices) for ind in network.output_qubits]
        node.indices[:] = network.output_qubits[:]
        permute_tensor(network, output_tensor, order)

        # Reshape the final tensor if a shape is specified by the user.
        if output_shape == "vector"
            reshape_tensor(network, output_tensor, [collect(1:length(node.indices))])
        elseif output_shape != ""
            reshape_tensor(network, output_tensor, output_shape)
        end
    end

    # save the final tensor under the name "result".
    save_output(network, output_tensor)
end

"""
    contract_network!(network::TensorNetworkCircuit,
                               plan::Array{Array{Symbol, 1}, 1},
                               output_shape::Union{String, Array{Int, 1}, 1})

Function to contract the given network according to the given contraction plan.
The resulting tensor will be given the shape described by 'output_shape'.
"""
function contract_network!(network::TensorNetworkCircuit,
                           plan::Array{Array{Symbol, 1}, 1},
                           output_shape::Union{String, Array{Int, 1}}="")

    # Loop through the plan and contract each tensor pair in sequence
    for (A, B) in plan
        contract_pair!(network, A, B)
    end

    # Permute the indices of the final tensor to have the correct order.
    output_tensor = Symbol("node_$(network.counters["node"])")
    node = network.nodes[output_tensor]
    if length(node.indices) != 0
        order = [findfirst(x->x==ind, node.indices) for ind in network.output_qubits]
        node.indices[:] = network.output_qubits[:]
        permute_tensor(network, output_tensor, order)

        # Reshape the final tensor if a shape is specified by the user.
        if output_shape == "vector"
            reshape_tensor(network, output_tensor, [collect(1:length(node.indices))])
        elseif output_shape != ""
            reshape_tensor(network, output_tensor, output_shape)
        end
    end

    # save the final tensor under the name "result".
    save_output(network, output_tensor)
end

"""
    contract_pair!(network::TensorNetworkCircuit,
                            edge::Symbol)

Contract a pair of nodes connected by the given edge.
"""
function contract_pair!(network::TensorNetworkCircuit,
                        edge::Symbol)
    # sometimes edges get contracted ahead of time if connecting two
    # tensors being contracted
    if edge in keys(network.edges)
        e = network.edges[edge]
        return contract_pair!(network, e.src, e.dst)
    end
end

"""
    contract_pair!(network::TensorNetworkCircuit,
                            A_label::Symbol,
                            B_label::Symbol)

Function to contract the nodes A and B of the network.
"""
function contract_pair!(network::TensorNetworkCircuit,
                        A_label::Symbol,
                        B_label::Symbol)

    # Should only contract different tensors, so skip contraction if A = B.
    if A_label == B_label
        return nothing
    end

    # Get the tensors being contracted and denote them A and B
    A = network.nodes[A_label]
    B = network.nodes[B_label]

    # Sort the indices of A and B into shared and unshared.
    common_indices, remaining_indices = sort_indices(A, B)

    # Create array of ncon indices for A and B.
    A_ncon_indices, B_ncon_indices = create_ncon_indices(A, B,
                                                         common_indices,
                                                         remaining_indices)

    # Get the dimensions of the contracted and open indices and use them
    # to record the costs of contraction
    dims_map = Dict{Symbol, Int}([A.indices; B.indices] .=> [A.dims; B.dims])
    contracted_dims = [dims_map[ind] for ind in common_indices]
    C_dims = [dims_map[ind] for ind in remaining_indices]
    record_compute_costs!(network.backend, C_dims, contracted_dims)

    # Create and add the new node to the tensor network.
    C_label = new_label!(network, "node")
    C = Node(remaining_indices, C_dims, C_label)
    network.nodes[C_label] = C

    # remove contracted edges
    for index in common_indices
        delete!(network.edges, index)
    end

    # replumb existing edges to point to new node
    # TODO: might become time consuming for large networks
    for index in remaining_indices
        e = network.edges[index]
        if e.src in [A_label, B_label]
            e.src = C_label
        elseif e.dst in [A_label, B_label]
            e.dst = C_label
        end
    end

    # delete the nodes that were contracted
    delete!(network.nodes, A_label)
    delete!(network.nodes, B_label)

    # Call lower level function to perform contraction on tensors
    contract_tensors(network, A_label, A_ncon_indices,
                     B_label, B_ncon_indices, C_label)

    C_label
end

# *************************************************************************** #
#                       Compression of a tensor network
# *************************************************************************** #

"""
    compress_tensor_chain!(network::TensorNetworkCircuit,
                                    nodes::Array{Symbol, 1};
                                    threshold::AbstractFloat=1e-13,
                                    max_rank::Int=0)

Compress a chain of tensors by given by the array of symbols. This is achieved
by peforming forward and backward sweeps where of compression operations on each
bond. Compression of each bond proceeds by discarding singular values and
corresponding axes with singular values below the given threshold or values
beyond the max_rank (max_rank zero corresponds to infinite rank)
"""
function compress_tensor_chain!(network::TensorNetworkCircuit,
                                nodes::Array{Symbol, 1};
                                threshold::AbstractFloat=1e-13,
                                max_rank::Int=0)

    # TODO: required that tensors in the chain only have virtual bonds between
    # consecutive tensors. Add a check that enforces this

    # forward pass
    for i in 1:(length(nodes) - 1)
        compress_bond!(network, nodes[i], nodes[i+1], threshold=threshold, max_rank=max_rank)
    end

    # now do a backward pass
    for i in (length(nodes) - 1):-1:1
        compress_bond!(network, nodes[i], nodes[i+1], threshold=threshold, max_rank=max_rank)
    end

end

"""
    compress_bond!(network::TensorNetworkCircuit,
                            node_1::Symbol,
                            node_2::Symbol;
                            threshold::AbstractFloat=1e-13,
                            max_rank::Int=0)

Compress the virtual bond between the given nodes
"""
function compress_bond!(network::TensorNetworkCircuit,
                        node_1::Symbol,
                        node_2::Symbol;
                        threshold::AbstractFloat=1e-13,
                        max_rank::Int=0)

    left_node = network.nodes[node_1]
    right_node = network.nodes[node_2]

    left_indices = setdiff(left_node.indices, right_node.indices)
    right_indices = setdiff(right_node.indices, left_node.indices)

    combined_node = contract_pair!(network, node_1, node_2)

    decompose_tensor!(network, combined_node, left_indices, right_indices,
                        left_label=node_1,
                        right_label=node_2,
                        threshold=threshold,
                        max_rank=max_rank
                        )
end

"""
    decompose_tensor!(tng::TensorNetworkCircuit,
                               node::Symbol
                               left_indices::Array{Symbol, 1},
                               right_indices::Array{Symbol, 1};
                               threshold::AbstractFloat=1e-15,
                               max_rank::Int=0,
                               left_label::Union{Nothing, Symbol}=nothing,
                               right_label::Union{Nothing, Symbol}=nothing)

Decompose a tensor into two smaller tensors
"""
function decompose_tensor!(network::TensorNetworkCircuit,
                           node_label::Symbol,
                           left_indices::Array{Symbol, 1},
                           right_indices::Array{Symbol, 1};
                           threshold::AbstractFloat=1e-13,
                           max_rank::Int=0,
                           left_label::Union{Nothing, Symbol}=nothing,
                           right_label::Union{Nothing, Symbol}=nothing)

    node = network.nodes[node_label]
    index_map = Dict([v => k for (k, v) in enumerate(node.indices)])
    left_positions = [index_map[x] for x in left_indices]
    right_positions = [index_map[x] for x in right_indices]

    # plumb these nodes back into the graph and delete the original
    B_label = (left_label === nothing) ? new_label!(network, "node") : left_label
    C_label = (right_label === nothing) ? new_label!(network, "node") : right_label

    chi = decompose_tensor!(network,
                            node_label,
                            left_positions,
                            right_positions;
                            threshold=threshold,
                            max_rank=max_rank,
                            left_label=B_label,
                            right_label=C_label)

    # Work out the dimensions of the indices.
    B_dims = [node.dims[i] for i in left_positions]
    C_dims = [node.dims[i] for i in right_positions]
    if chi > 0
        virtual_dim = chi
    else
        left_dim = prod(B_dims)
        right_dim = prod(C_dims)
        virtual_dim = min(left_dim, right_dim)
        if max_rank > 0
            virtual_dim = min(virtual_dim, max_rank)
        end
    end
    append!(B_dims, virtual_dim)
    prepend!(C_dims, virtual_dim)

    index_label = new_label!(network, "index")
    B_node = Node(vcat(left_indices, [index_label,]), B_dims, B_label)
    C_node = Node(vcat([index_label,], right_indices), C_dims, C_label)

    # add the nodes
    network.nodes[B_label] = B_node
    network.nodes[C_label] = C_node

    # remap edge endpoints
    for index in left_indices
        if network.edges[index].src == node_label
            network.edges[index].src = B_label
        elseif network.edges[index].dst == node_label
            network.edges[index].dst = B_label
        end
    end
    for index in right_indices
        if network.edges[index].src == node_label
            network.edges[index].src = C_label
        elseif network.edges[index].dst == node_label
            network.edges[index].dst = C_label
        end
    end

    # add new edge
    network.edges[index_label] = Edge(B_label, C_label, nothing, true)

    delete_tensor!(network, node_label)
    delete!(network.nodes, node_label)

    (B_label, C_label)
end

"""
    contract_mps_tensor_network_circuit(network::TensorNetworkCircuit;
                                                 max_bond::Int=2,
                                                 threshold::AbstractFloat=1e-13,
                                                 max_rank::Int=0)

Contract a tensor network representing a quantum circuit using MPS techniques
"""
function contract_mps_tensor_network_circuit!(network::TensorNetworkCircuit;
                                              max_bond::Int=2,
                                              threshold::AbstractFloat=1e-13,
                                              max_rank::Integer=0,
                                              pbc::Bool=false)
    # identify the MPS nodes as the input nodes
    qubits = length(network.input_qubits)
    mps_nodes = [network.edges[x].src for x in network.input_qubits]
    @assert all([x !== nothing for x in mps_nodes]) "Input qubit values must be set"

    layer_nodes = Dict{Int, Array{Symbol,1}}()
    for (k, v) in pairs(network.node_layers)
        if haskey(layer_nodes, v)
            push!(layer_nodes[v], k)
        else
            layer_nodes[v] = [k]
        end
    end

    gate_layers = [k for k in sort(collect(keys(layer_nodes))) if k > 0]
    if haskey(layer_nodes, -1)
        gate_layers = vcat(gate_layers, [-1])
    end

    for gate_layer in gate_layers
        nodes = layer_nodes[gate_layer]

        updated_indices = Array{Int, 1}(undef, length(nodes))
        for (i, node) in enumerate(nodes)
            # find the input node to connect to it
            input_node = inneighbours(network, node)[1]
            @assert input_node in mps_nodes "$input_node not in mps nodes"

            idx = findfirst(x -> x == input_node, mps_nodes)
            mps_nodes[idx] = contract_pair!(network, input_node, node)
            updated_indices[i] = idx
         end

        sort!(updated_indices)
        for i in 1:(length(updated_indices)-1)
            distance = updated_indices[i+1] - updated_indices[i]
            if pbc
                distance = min(distance, updated_indices[i] + length(mps_nodes)  - updated_indices[i+1])
            end
            @assert distance == 1 "Gates between non-neighboring qubits"
            compress_bond!(network, mps_nodes[updated_indices[i]], mps_nodes[updated_indices[i+1]], threshold=threshold, max_rank=max_rank)
        end
    end

    # save the final mps tensors
    for (i, node) in enumerate(mps_nodes)
        save_output(network, node, String(node))
    end
    mps_nodes
end

"""
    calculate_mps_amplitudes!(network::TensorNetworkCircuit,
                                      mps_nodes::Array{Symbol, 1},
                                      result::String="result")

Calculate amplitudes from an MPS state
"""
function calculate_mps_amplitudes!(network::TensorNetworkCircuit,
                                   mps_nodes::Array{Symbol, 1},
                                   result::String="result")
    output_node = mps_nodes[1]
    for node in mps_nodes[2:end]
        output_node = contract_pair!(network, output_node, node)
    end
    permute_tensor(network, output_node, network.qubit_ordering)
    reshape_tensor(network, output_node, [collect(1:length(mps_nodes))])
    save_output(network, output_node, result)
end

"""
    function contract_tensor_network_circuit_with_compression(network::TensorNetworkCircuit;
                                                              max_bond::Integer=2,
                                                              threshold::AbstractFloat=1e-13,
                                                              max_rank::Integer=0)

Contract a tensor network representing a quantum circuit tensor network compresson on bonds
"""
function contract_tensor_network_circuit_with_compression!(network::TensorNetworkCircuit;
                                                           max_bond::Integer=2,
                                                           threshold::AbstractFloat=1e-13,
                                                           max_rank::Integer=0)
    # identify the network nodes as the input nodes
    qubits = length(network.input_qubits)
    network_nodes = [network.edges[x].src for x in network.input_qubits]
    @assert all([x !== nothing for x in network_nodes]) "Input qubit values must be set"

    layer_nodes = Dict{Int64, Array{Symbol,1}}()
    for (k, v) in pairs(network.node_layers)
        if haskey(layer_nodes, v)
            push!(layer_nodes[v], k)
        else
            layer_nodes[v] = [k]
        end
    end

    gate_layers = [k for k in sort(collect(keys(layer_nodes))) if k > 0]
    if haskey(layer_nodes, -1)
        gate_layers = vcat(gate_layers, [-1])
    end

    for gate_layer in gate_layers
        nodes = layer_nodes[gate_layer]

        updated_indices = Array{Int64, 1}(undef, length(nodes))
        for (i, node) in enumerate(nodes)
            # find the input node to connect to it
            input_node = inneighbours(network, node)[1]
            @assert input_node in network_nodes "$input_node not in network nodes"

            idx = findfirst(x -> x == input_node, network_nodes)
            network_nodes[idx] = contract_pair!(network, input_node, node)
            updated_indices[i] = idx
         end

        if length(updated_indices) > 1
            for i in 1:(length(updated_indices)-1)
                compress_bond!(network, network_nodes[updated_indices[i]], network_nodes[updated_indices[i+1]], threshold=threshold, max_rank=max_rank)
            end
        end
    end

    # save the final mps tensors
    for (i, node) in enumerate(network_nodes)
        save_output(network, node, String(node))
    end
    network_nodes
end


text:
function LineLine(l1::Geometry, l2::Geometry)
    #Ips = Array{Tuple{Point2D{Float64}, Float64, Float64}, 1}()
    dx_cx = l2.Ep.x - l2.Sp.x
    dy_cy = l2.Ep.y - l2.Sp.y
    bx_ax = l1.Ep.x - l1.Sp.x
    by_ay = l1.Ep.y - l1.Sp.y
    Δ = bx_ax * dy_cy - by_ay * dx_cx
    if abs(Δ) < 1e-4 # colinear
        if l1.Ep.x == l2.Sp.x # l2 follows l1
            return (intersect=true, type="connect", Ip=Point4To2D(l1.Ep), Dist1=l1.Length, Dist2=0.0)
        elseif l2.Ep.x == l1.Sp.x # l1 follows l2
            return (intersect=true, type="connect", Ip=Point4To2D(l1.Sp), Dist1=0.0, Dist2=l2.Length)
        elseif onSegment(l1.Sp, l1.Ep, l2.Sp)
            return (intersect=true, type="overlap", Ip=Point4To2D(l2.Sp), Dist1=Distance2D(l1.Sp, l2.Sp), Dist2=0.0)
        elseif onSegment(l2.Sp, l2.Ep, l1.Sp)
            return (intersect=true, type="overlap", Ip=Point4To2D(l1.Sp), Dist1=0.0, Dist2=Distance2D(l2.Sp, l1.Sp))
        else
            return (intersect=false, type="none")
        end
    else # not colinear
        ax_cx = l1.Sp.x - l2.Sp.x
        ay_cy = l1.Sp.y - l2.Sp.y
        r = (ay_cy * dx_cx - ax_cx * dy_cy) / Δ
        s = (ay_cy * bx_ax - ax_cx * by_ay) / Δ
        if (0.0 <= r <= 1.0 &&  0.0 <= s <= 1.0)  # intersection within 2 line-segments
            Ix = l1.Sp.x + r * bx_ax
            Iy = l1.Sp.y + r * by_ay
            Ipc = Point2D(Ix, Iy)
            dist1 = Distance2D(Point4To2D(l1.Sp), Ipc)
            dist2 = Distance2D(Point4To2D(l2.Sp), Ipc)
            return (intersect=true, type="cross", Ip=Ipc, Dist1=dist1, Dist2=dist2)
        else
            return (intersect=false, type="none")
        end
    end
end

function onSegment(Sp::Point4D{Float64}, Ep::Point4D{Float64}, Pi::Point4D{Float64})
    if (Pi.x <= max(Sp.x, Ep.x) && Pi.x >= min(Sp.x, Ep.x) &&
        Pi.y <= max(Sp.y, Ep.y) && Pi.y >= min(Sp.y, Ep.y))
        return true
    else
        return false
    end
end

ext: jl


text:
# included in runtests.jl
using Base.Test
using HTTP
using HttpServer
import WebSockets:  is_upgrade,
                    upgrade
function echows(req, ws)
    @test origin(req) == ""
    @test target(req) == "/"
    @test subprotocol(req) == ""
    while true
        data, success = readguarded(ws)
        !success && break
        !writeguarded(ws, data) && break
    end
end

const port_HTTP = 8000
const port_HTTP_ServeWS = 8001
const port_HttpServer = 8081
const TCPREF = Ref{Base.TCPServer}()

# Start HTTP listen server on port $port_HTTP"
tas = @schedule HTTP.listen("127.0.0.1", port_HTTP, tcpref = TCPREF) do s
    if WebSockets.is_upgrade(s.message)
        WebSockets.upgrade(echows, s)
    end
end
while !istaskstarted(tas);yield();end

# Start HttpServer on port $port_HttpServer
const server = Server(HttpHandler((req, res)->Response(200)),
                      WebSocketHandler(echows));
tas = @schedule run(server, port_HttpServer)
while !istaskstarted(tas);yield();end

# Start HTTP ServerWS on port $port_HTTP_ServeWS
server_WS = WebSockets.ServerWS(
    HTTP.HandlerFunction(req-> HTTP.Response(200)),
    WebSockets.WebsocketHandler(echows))

tas = @schedule WebSockets.serve(server_WS, "127.0.0.1", port_HTTP_ServeWS)
while !istaskstarted(tas);yield();end

const servers = [
    ("HTTP",        "ws://127.0.0.1:$(port_HTTP)"),
    ("HttpServer",  "ws://127.0.0.1:$(port_HttpServer)"),
    ("HTTTP ServerWS",  "ws://127.0.0.1:$(port_HTTP_ServeWS)"),
    ("ws",          "ws://echo.websocket.org"),
    ("wss",         "wss://echo.websocket.org")]

const lengths = [0, 125, 126, 127, 2000]

for (s, url) in servers, len in lengths, closestatus in [false, true]
    len == 0 && contains(url, "echo.websocket.org") && continue
    info("Testing client -> server at $(url), message length $len")
    test_str = randstring(len)
    forcecopy_str = test_str |> collect |> copy |> join
    WebSockets.open(url) do ws
        print(" -Foo-")
        write(ws, "Foo")
        @test String(read(ws)) == "Foo"
        print(" -Ping-")
        send_ping(ws)
        print(" -String length $len-\n")
        write(ws, test_str)
        @test String(read(ws)) == forcecopy_str
        closestatus && close(ws, statusnumber = 1000)
    end
    sleep(0.2)
end

# make a very simple http request for the servers with defined http handlers
resp = HTTP.request("GET", "http://127.0.0.1:$(port_HTTP_ServeWS)")
@test resp.status == 200
resp = HTTP.request("GET", "http://127.0.0.1:$(port_HttpServer)")
@test resp.status == 200

# Close the servers
close(TCPREF[])
close(server.http.sock)
put!(server_WS.in, HTTP.Servers.KILL)


text:
module CUDD

using Compat

# cudd_path = joinpath(dirname(@__FILE__), "..", "deps", "cudd-3.0.0", "cudd", ".libs")
const depfile = joinpath(@__DIR__, "..", "deps", "deps.jl")

if isfile(depfile)
    include(depfile)
else
    error("CUDD not properly installed. Please run Pkg.build(\"CUDD\")")
end

include("ADD.jl")
include("ADD_apply.jl")
include("ADD_find.jl")
include("ADD_abs.jl")
include("ADD_inverse.jl")
include("ADD_ITE.jl")
include("ADD_negate.jl")
include("ADD_walsh.jl")

end # module


text:
export produce_or_load, tagsave, @tagsave, safesave

"""
    produce_or_load([prefix="",] c, f; kwargs...) -> file, path
Let `s = savename(prefix, c, suffix)`.
If a file named `s` exists then load it and return it, along
with the path that it is saved at, `s` (default behavior).

If the file does not exist then call `file = f(c)`, save `file` as
`s` and then return `file, s`.
The function `f` must return a dictionary.
The macros [`@dict`](@ref) and [`@strdict`](@ref) can help with that.

## Keywords
* `tag = true` : Add the Git commit of the project in the saved file.
* `gitpath = projectdir()` : Path to search for a Git repo.
* `suffix = "bson"` : Used in `savename`.
* `force = false` : If `true` then don't check if file `s` exists and produce
  it and save it anyway.
* `loadfile = true` : If `false`, this function does not actually load the
  file, but only checks if it exists. The return value in this case is always
  `nothing, s`, regardless of whether the file must be produced or not.
* `verbose = true` : print info about the process.
* `kwargs...` : All other keywords are propagated to `savename`.

See also [`savename`](@ref) and [`tag!`](@ref).
"""
produce_or_load(c, f; kwargs...) = produce_or_load("", c, f; kwargs...)
function produce_or_load(prefix::String, c, f;
    tag::Bool = true, gitpath = projectdir(), loadfile = true,
    suffix = "bson", force = false, verbose = true, kwargs...)

    s = savename(prefix, c, suffix; kwargs...)

    if !force && isfile(s)
        if loadfile
            file = wload(s)
            return file, s
        else
            return nothing, s
        end
    else
        if force
            verbose && @info "Producing file $s now..."
        else
            verbose && @info "File $s does not exist. Producing it now..."
        end
        file = f(c)
        try
            mkpath(dirname(s))
            if tag
                tagsave(s, file, false, gitpath)
            else
            wsave(s, copy(file))
        end
        catch er
            @warn "Could not save file, got error $er. "*
            "\nReturning the file if `loadfile=true`."
        end
        if loadfile
            return file, s
        else
            return nothing, s
        end
    end
end

"""
    tagsave(file::String, d::Dict [, safe = false, gitpath = projectdir()])
First [`tag!`](@ref) dictionary `d` and then save `d` in `file`.
If `safe = true` save the file using [`safesave`](@ref).
"""
tagsave(file, d, p::String) = tagsave(file, d, false, p)
function tagsave(file, d, safe = false, gitpath = projectdir(), s = nothing)
    d2 = tag!(d, gitpath, s)
    mkpath(dirname(file))
    if safe
        safesave(file, copy(d2))
    else
        wsave(file, copy(d2))
    end
    return d2
end

"""
    @tagsave(file::String, d::Dict [, safe = false, gitpath = projectdir()])
Same as [`tagsave`](@ref) but also add a field `script` that records
the local path of the script that called `@tagsave`, see [`@tag!`](@ref).
"""
macro tagsave(file, d, safe::Bool = false, gitpath = projectdir())
    s = QuoteNode(__source__)
    :(tagsave($(esc(file)), $(esc(d)), $(esc(safe)), $(esc(gitpath)), $s))
end

################################################################################
#                          Backup files before saving                          #
################################################################################

# Implementation inspired by behavior of GROMACS
"""
    safesave(filename, data)

Safely save `data` in `filename` by ensuring that no existing files
are overwritten. Do this by renaming already existing data with a backup-number
ending like `#1, #2, ...`. For example if `filename = test.bson`, the first
time you `safesave` it, the file is saved normally. The second time
the existing save is renamed to `test_#1.bson` and a new file `test.bson`
is then saved.

If a backup file already exists then its backup-number is incremented
(e.g. going from `#2` to `#3`). For example safesaving `test.bson` a third time
will rename the old `test_#1.bson` to `test_#2.bson`, rename the old
`test.bson` to `test_#1.bson` and then save a new `test.bson` with the latest
`data`.

See also [`tagsave`](@ref).
"""
function safesave(f, data)
    recursively_clear_path(f)
    mkpath(dirname(f))
    wsave(f, data)
end

#take a path of a results file and increment its prefix backup number
function increment_backup_num(filepath)
    path, filename = splitdir(filepath)
    fname, suffix = splitext(filename)
    m = match(r"^(.*)_#([0-9]+)$", fname)
    if m == nothing
        return joinpath(path, "$(fname)_#1$(suffix)")
    end
    newnum = string(parse(Int, m.captures[2]) +1)
    return joinpath(path, "$(m.captures[1])_#$newnum$(suffix)")
end

#recursively move files to increased backup number
function recursively_clear_path(cur_path)
    isfile(cur_path) || return
    new_path=increment_backup_num(cur_path)
    if isfile(new_path)
        recursively_clear_path(new_path)
    end
    mv(cur_path, new_path)
end

################################################################################
#                    Compliment to dict_list: tmpsave                          #
################################################################################
export tmpsave
using Random
"""
    tmpsave(dicts::Vector{Dict} [, tmp]; kwargs...) -> r
Save each entry in `dicts` into a unique temporary file in the directory `tmp`.
Then return the list of file names (relative to `tmp`) that were used
for saving each dictionary.

`tmp` defaults to `projectdir("_research/tmp")`.

See also [`dict_list`](@ref).

## Keywords
* `l = 8` : number of characters in the random string.
* `prefix = ""` : prefix each temporary name with this.
* `suffix = "bson"` : ending of the temporary names (no need for the dot).
"""
function tmpsave(dicts, tmp = projectdir("_research/tmp");
    l = 8, suffix = "bson", prefix = "")

    mkpath(tmp)
    n = length(dicts)
    existing = readdir(tmp)
    r = String[]
    i = 0
    while i < n
        x = prefix*randstring(l)*"."*suffix
        while x ∈ r || x ∈ existing
            x = prefix*randstring(l)*"."*suffix
        end
        i += 1
        push!(r, x)
        wsave(joinpath(tmp, x), copy(dicts[i]))
    end
    r
end


text:
using Symbolics

include("./kinematics.jl")
using .Kinematics

N = 2
@variables Q[1:3N]
@variables Ξ[1:N]

println(Q[1])

Ps, Rs = calc_local(N, Q, Ξ);

Φs, Θs = calc_global(N, Ps, Rs);


text:
using ExamplePkg
using Documenter

makedocs(;
    modules=[ExamplePkg],
    authors="Al Wonder",
    repo="https://github.com/alyiwonder/ExamplePkg.jl/blob/{commit}{path}#L{line}",
    sitename="ExamplePkg.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://alyiwonder.github.io/ExamplePkg.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/alyiwonder/ExamplePkg.jl",
)


text:
using Flux, Flux.Data.MNIST, Statistics
using Flux: onehotbatch, onecold, crossentropy, throttle
using Base.Iterators: repeated

using Jax

#
N = 5000

# Classify MNIST digits with a simple multi-layer-perceptron
imgs = MNIST.images()[1:N]
# Stack images into one large batch
X = hcat(float.(reshape.(imgs, :))...) |> tojax

labels = MNIST.labels()[1:N]
# One-hot-encode the labels
Y = onehotbatch(labels, 0:9) |>collect |> tojax

m = Chain(
  Dense(28^2, 32, relu),
  Dense(32, 10),
  softmax) |> tojax

loss(m, x, y) = crossentropy(m(x), y)

dataset = repeated((X, Y), 200)
opt = ADAM()

#

vals, funs = Jax.flatten(m)
loss(w::Tuple, x, y) = loss(Jax.unflatten(funs, w), x, y)
grad(vals, X, Y) = jax.grad(loss)(vals, X, Y)

function Flux.update!(opt, xs::Tuple, gs::Tuple)
  for (x,g) in zip(xs,gs)
    Flux.update!(opt, x, g)
  end
end

function do_update(vals, X, Y)
  g = grad(vals, X, Y)
  Flux.update!(opt, vals, g)
  return vals
end

do_update(vals, X, Y)

j_du = jax.jit(do_update)

for data=dataset
  global vals = j_du(vals, data...)
  show(loss(vals, data...))
end


text:
include("../envs/RLEnvironments.jl")
using .RLEnvironments
using Plots
using Flux
using Statistics: mean

function episode(model, env, maxt, γ)
    R = zero(γ)
    reset!(env)
    for t in 1:maxt
        a = model(env.current)
        _, r, done = step!(env, a)
        R += γ^t * r
        done && break
    end
    return R
end

function run!(model, env; epochs=100, maxt=100, popsize=50, elitesize=10, σ=5f-1, γ=99f-2)
    rewards = zeros(Float32, epochs)
    p, re = Flux.destructure(model)
    W = σ .* randn(eltype(p), length(p))
    for i in 1:epochs
        Ws = [W + σ .* randn(eltype(p), length(p)) for _ in 1:popsize]
        Rs = map(Wi -> episode(re(Wi), env, maxt, γ), Ws)
        elids = sortperm(Rs, rev=true)[1:elitesize]
        W = mean(Ws[elids])

        rewards[i] = mean(Rs[elids])
    end
    return rewards
end

env = Pendulum{Float32}()

model = Chain(
    Dense(length(env.observationspace), 32, relu),
    Dense(32, length(env.actionspace), tanh)
)

r = run!(model, env)
plt = plot(r, labels="Reward")
display(plt)
#savefig(plt, "cem")

#=
function setweights!(model, W)
    curr = 1
    for layer in model.layers
        Wsize = size(layer.W)
        bsize = size(layer.b)
        Wlen = *(Wsize...)
        blen = *(bsize...)
        layer.W .= reshape(W[curr:curr+Wlen-1], Wsize)
        layer.b .= reshape(W[(curr+Wlen):(curr+Wlen+blen-1)], bsize)
        curr += Wlen + blen
    end
end

numparams = (length(env.observationspace) + 1) * 32 + (32 + 1) * length(env.actionspace)

Rs = map(Wi -> (setweights!(model, Wi); episode(model, env, maxt, γ)), Ws)
=#


text:
struct GpuIndirectSolver <: LinearSolver end

if haskey(ENV, "JULIA_SCS_LIBRARY_PATH")
    @isdefined(libscsgpuindir) && push!(available_solvers, GpuIndirectSolver)
else
    import SCS_GPU_jll
    const gpuindirect = SCS_GPU_jll.libscsgpuindir
    push!(available_solvers, GpuIndirectSolver)
end

scsint_t(::Type{GpuIndirectSolver}) = Int32

function scs_set_default_settings(
    ::Type{GpuIndirectSolver},
    data::ScsData{Int32},
)
    return ccall(
        (:scs_set_default_settings, gpuindirect),
        Cvoid,
        (Ptr{Cvoid},),
        data,
    )
end

function scs_init(
    ::Type{GpuIndirectSolver},
    data::ScsData{Int32},
    cone::ScsCone{Int32},
    info::ScsInfo{Int32},
)
    return ccall(
        (:scs_init, gpuindirect),
        Ptr{Cvoid},
        (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
        data,
        cone,
        info,
    )
end

function scs_solve(
    ::Type{GpuIndirectSolver},
    p_work::Ptr{Cvoid},
    data::ScsData{Int32},
    cone::ScsCone{Int32},
    solution::ScsSolution,
    info::ScsInfo{Int32},
)
    return ccall(
        (:scs_solve, gpuindirect),
        Int32,
        (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
        p_work,
        data,
        cone,
        solution,
        info,
    )
end

function scs_finish(::Type{GpuIndirectSolver}, p_work::Ptr{Cvoid})
    return ccall((:scs_finish, gpuindirect), Cvoid, (Ptr{Cvoid},), p_work)
end


text:
struct AutoReverseDiff <: AbstractADType end

function instantiate_function(f, x, ::AutoReverseDiff, p=SciMLBase.NullParameters(), num_cons = 0)
    num_cons != 0 && error("AutoReverseDiff does not currently support constraints")

    _f = (θ, args...) -> first(f.f(θ,p, args...))

    if f.grad === nothing
        grad = (res, θ, args...) -> ReverseDiff.gradient!(res, x -> _f(x, args...), θ, ReverseDiff.GradientConfig(θ))
    else
        grad = f.grad
    end

    if f.hess === nothing
        hess = function (res, θ, args...)
            if res isa DiffResults.DiffResult
                DiffResults.hessian!(res, ForwardDiff.jacobian(θ) do θ
                                                ReverseDiff.gradient(x -> _f(x, args...), θ)[1]
                                            end)
            else
                res .=  ForwardDiff.jacobian(θ) do θ
                    ReverseDiff.gradient(x ->_f(x, args...), θ)
                  end
            end
        end
    else
        hess = f.hess
    end

    if f.hv === nothing
        hv = function (H,θ,v, args...)
            _θ = ForwardDiff.Dual.(θ,v)
            res = DiffResults.GradientResult(_θ)
            grad(res, _θ, args...)
            H .= getindex.(ForwardDiff.partials.(DiffResults.gradient(res)),1)
        end
    else
        hv = f.hv
    end

    return OptimizationFunction{false,AutoReverseDiff,typeof(f),typeof(grad),typeof(hess),typeof(hv),Nothing,Nothing,Nothing}(f,AutoReverseDiff(),grad,hess,hv,nothing,nothing,nothing)
end


text:
"""
  propose_new_component(old_mixture::Distribution, ll::Function; algorithm="plain", method="MLE", n=10, n_est_c=10000)

Propose new component to `old_mixture` for approximating `true_density`.

* `ll`: log of true_density
* `algorithm`: `plain`, `ratio`, `rejection` or `demean`
* `method`: method for fitting the new component with weighted samples. Default to `MLE`.
* `n`: number of particles
* `n_est_c`: number of particles for estimating the normalization factor for `demean` algorithm
"""
function propose_new_component(old_mixture::Distribution, ll::Function;
                               algorithm="plain", method="MLE", n=10, n_est_c=10000)
    if algorithm == "plain"
        return propose_new_component_plain(old_mixture, ll, n=n, method=method)
    elseif algorithm == "ratio"
        return propose_new_component_ratio(old_mixture, ll, n=n, method=method)
    elseif algorithm == "rejection"
        return propose_new_component_rejection(old_mixture, ll, n=n, method=method)
    elseif algorithm == "demean"
        return propose_new_component_demean(old_mixture, ll, n=n, n_est_c=n_est_c, method=method)
    else
        error("wrong method")
    end
end

"fit new component of distribution D with weighted samples"
function fit_new_component{D <: Distribution, T <: AbstractFloat}(d::Type{D}, x::Array{T}, w::Vector{Float64}; method="MLE")
    w = w / sum(w)
    if method == "MLE"
        return fit_mle(d, x, w)
    else
        error("wrong method")
    end
end

"equal weight variant"
function fit_new_component{D <: Distribution, T <: AbstractFloat}(d::Type{D}, x::Array{T}; method="MLE")
    n = ndims(x)==1 ? length(x) : size(x, 2)
    return fit_new_component(d, x, ones(n) / n, method=method)
end

"proposal: particles ~ true distribution"
function propose_new_component_plain(old_mixture::Distribution, ll::Function; n=10, method="MLE")
  pt = sample_importance(old_mixture, ll, n=n)
  return fit_new_component(get_component_distribution(old_mixture), pt.x, pt.w, method=method)
end

"proposal: particles ~ exp(log true - log approx)"
function propose_new_component_ratio(old_mixture::Distribution, ll::Function; n=10, method="MLE")
  pt = sample_importance(old_mixture, ll, n=n, normalized=false)
  w = pt.w ./ pdf(old_mixture, pt.x)
  return fit_new_component(get_component_distribution(old_mixture), pt.x, w, method=method)
end

"max density for univariate normal"
function max_density(d::Normal)
  return pdf(d, d.μ)
end

"max density for multivariate normal"
function max_density(d::MultivariateNormal)
    return pdf(d, d.μ)
end

"max density for a mixture model"
function max_density(d::AbstractMixtureModel)
  K = length(probs(d))
  return maximum([pdf(d, d.components[i].μ) for i in 1:K])
end

"proposal: particles ~ true distribution, but rejecting those around the max mode of the approx"
function propose_new_component_rejection(old_mixture::Distribution, ll::Function; n=10, method="MLE")
  γ = max_density(old_mixture)
  pt = sample_importance(old_mixture, ll, n=n)
  rej = 1.0 - pdf(old_mixture, pt.x) ./ γ
  rej[rej .< 0.0] = 0.0 # small numerical error may make it negative
  w = pt.w .* rej
  return fit_new_component(get_component_distribution(old_mixture), pt.x, w, method=method)
end

"proposal: `demean estimator` particle ~ (true density - approx) ∨ 0"
function propose_new_component_demean(old_mixture::Distribution, ll::Function; n=10, n_est_c=1000, method="MLE")
    mean_w_est = estimate_normalizing_factor(old_mixture, ll, n=n_est_c)
    return propose_new_component_demean(old_mixture, ll, mean_w_est; n=n, method=method)
end

"proposal: `demean estimator` particle ~ (true density - approx) ∨ 0, with mean_w as an input"
function propose_new_component_demean(old_mixture::Distribution, ll::Function, mean_w_est::AbstractFloat; n=10, method="MLE")
    w = zeros(n)
    p = length(old_mixture)
    if p==1
        x = zeros(n)
    else
        x = zeros((p, n))
    end
    i = 1
    while i <= n
        if p==1
            x[i] = rand(old_mixture)
            w[i] = exp(ll(x[i]) - logpdf(old_mixture, x[i])) - mean_w_est # demean
        else
            x[:, i] = rand(old_mixture)
            w[i] = exp(ll(x[:, i]) - logpdf(old_mixture, x[:, i])) - mean_w_est # demean
        end

        if w[i] > 0
            i += 1
        end
    end
    return fit_new_component(get_component_distribution(old_mixture), x, w, method=method)
end

max_stars_repo_name: richardkwo/boosting-VB

max_issues_repo_name: richardkwo/boosting-VB

max_forks_repo_name: richardkwo/boosting-VB


text:
# ----------------------
# -- read             --
# ----------------------
function read_nodenum(skipnum)
    """
    xmax : 仮想セルも含めたnodeのxの数
    ymax : 仮想セルも含めたnodeのyの数
    """
    fff=[]
    open("grid/nodesnum", "r") do f
        fff=read(f,String)
    end
    fff=split(fff,"\n",keepempty=false)
    num_nodes=length(fff)-skipnum

    for i in 1+skipnum:length(fff)
        fff[i]=replace(fff[i]," \r" => "")
    end

    temp = split(fff[2]," ")
    xmax = parse(Int64,temp[1])
    ymax = parse(Int64,temp[2])

    return xmax, ymax
end

function read_nodes(skipnum,xmax,ymax)
    """
    nodes[i][j][k]
    i : x点の番号
    j : y点の番号
    k=1 : 点のx座標
    k=2 : 点のy座標
    """

    fff=[]
    open("grid/nodes", "r") do f
        fff=read(f,String)
    end
    fff=split(fff,"\n",keepempty=false)
    num_nodes=length(fff)-skipnum

    for i in 1+skipnum:length(fff)
        fff[i]=replace(fff[i]," \r" => "")
    end

    nodes = zeros(xmax,ymax,2)
    for i in 1:num_nodes
        temp=split(fff[i+skipnum]," ")

        xnum = parse(Int64,temp[1])
        ynum = parse(Int64,temp[2])
        nodes[xnum,ynum,1]=parse(Float64,temp[3])
        nodes[xnum,ynum,2]=parse(Float64,temp[4])
    end
    return nodes
end

function read_nodes_vtk(skipnum)

    fff=[]
    open("grid/nodes_forvtk", "r") do f
        fff=read(f,String)
    end
    fff=split(fff,"\n",keepempty=false)
    num_nodes=length(fff)-skipnum

    for i in 1+skipnum:length(fff)
        fff[i]=replace(fff[i]," \r" => "")
    end

    nodes=zeros(num_nodes,3)
    for i in 1:num_nodes
        temp=split(fff[i+skipnum]," ")

        # x = parse(Float64,temp[1])
        # y = parse(Float64,temp[2])
        # z = parse(Float64,temp[3])
        x = parse(Float64,temp[2])
        y = parse(Float64,temp[3])
        z = 0.0
        nodes[i,1] = x
        nodes[i,2] = y
        nodes[i,3] = z
    end
    return nodes
end

function read_elements_vtk(skipnum)
    fff=[]
    open("grid/element_forvtk", "r") do f
        fff=read(f,String)
    end
    fff=split(fff,"\n",keepempty=false)
    num_elements=length(fff)-skipnum

    for i in 1+skipnum:length(fff)
        fff[i]=replace(fff[i]," \r" => "")
    end

    elements = zeros(Int64,num_elements,4)
    for i in 1:num_elements
        temp=split(fff[i+skipnum]," ")

        elements[i,1] = parse(Int64,temp[2])
        elements[i,2] = parse(Int64,temp[3])
        elements[i,3] = parse(Int64,temp[4])
        elements[i,4] = parse(Int64,temp[5])

    end
    return elements
end

function read_result(skipnum)
    fff=[]
    open("test333.dat", "r") do f
        fff=read(f,String)
    end
    fff=split(fff,"\n",keepempty=false)
    num_point=length(fff)-skipnum

    for i in 1+skipnum:length(fff)
        fff[i]=replace(fff[i]," \r" => "")
    end

    readQ = zeros(num_point,5)
    for i in 1:num_point
        temp=split(fff[i+skipnum]," ")

        k = 1
        for j in 1:length(temp)
            if temp[j] != ""
                readQ[i,k] = parse(Float64,temp[j])
                k += 1
            end
        end

    end
    return readQ
end

function read_allgrid()
    skip=1
    xmax,ymax = read_nodenum(skip)
    nodes     = read_nodes(skip,xmax,ymax)
    nodes_vtk = read_nodes_vtk(skip)
    elements  = read_elements_vtk(skip)
    readQ     = read_result(skip)
    println("fin read grid")
    return xmax,ymax,nodes,nodes_vtk,elements,readQ
end

max_stars_repo_name: hide-dog/Qtree

max_issues_repo_name: hide-dog/Qtree

max_forks_repo_name: hide-dog/Qtree


text:
function ∇1e(BS::BasisSet, compute::String, iA, T::DataType = Float64)
    # Pre allocate output
    out = zeros(T, BS.nbas, BS.nbas, 3)
    return ∇1e!(out, BS, compute, iA, T)
end

function ∇1e!(out, BS::BasisSet, compute::String, iA, T::DataType = Float64)

    if size(out) != (BS.nbas, BS.nbas, 3)
        throw(DimensionMismatch("Size of the output array needs to be (nbas, nbas, 3)"))
    end

    if compute == "overlap"
        libcint_1e! =  cint1e_ipovlp_sph!
    elseif compute == "kinetic"
        libcint_1e! =  cint1e_ipkin_sph!
    elseif compute == "nuclear"
        return ∇nuclear(BS, iA)
    end

    A = BS.atoms[iA]

    # Get shell index (s0) where basis of the desired atom start
    s0 = 0
    for a in eachindex(BS.atoms)
        if BS.atoms[a] != A
            s0 += length(BS.basis[a])
        else
            break
        end
    end

    # Shell indexes for basis in the atom A (C notation: Starts from 0)
    Ashells = [(s0 + i -1) for i = 1:length(BS.basis[iA])]
    notAshells = Int[]
    for i = 1:BS.nshells
        if !((i-1) in Ashells)
            push!(notAshells, i-1)
        end
    end

    lvals = [Libcint.CINTcgtos_spheric(i-1, BS.lc_bas) for i = 1:BS.nshells]
    ao_offset = [sum(lvals[1:(i-1)]) for i = 1:BS.nshells]
    Lmax = maximum(lvals)
    buf_arrays = [zeros(Cdouble, 3*Lmax^2) for _ = 1:Threads.nthreads()]
    @sync for i in Ashells
        Threads.@spawn begin
        @inbounds begin
            Li = lvals[i+1]
            ioff = ao_offset[i+1]
            buf = buf_arrays[Threads.threadid()]
            for j in notAshells
                Lj = lvals[j+1]
                joff = ao_offset[j+1]
                Lij = Li*Lj
                # Call libcint
                libcint_1e!(buf, Cint.([i,j]), BS.lc_atoms, BS.natoms, BS.lc_bas, BS.nbas, BS.lc_env)
                I = (ioff+1):(ioff+Li)
                J = (joff+1):(joff+Lj)

                # Get strides for each cartesian
                for k in 1:3
                    r = (1+Lij*(k-1)):(k*Lij)
                    @views ∇k = buf[r]
                    out[I,J,k] .-= reshape(∇k, Int(Li), Int(Lj))
                end

                # Copy over the transpose
                out[J,I,:] .+= permutedims(out[I,J,:], (2,1,3))
            end
        end #inbounds
        end #spawn
    end #sync
    return out
end

function ∇overlap(BS::BasisSet, iA, T::DataType = Float64)
    return ∇1e(BS, "overlap", iA, T)
end

function ∇overlap!(out, BS::BasisSet, iA, T::DataType = Float64)
    return ∇1e!(out, BS, "overlap", iA, T)
end

function ∇kinetic(BS::BasisSet, iA, T::DataType = Float64)
    return ∇1e(BS, "kinetic", iA, T)
end

function ∇kinetic!(out, BS::BasisSet, iA, T::DataType = Float64)
    return ∇1e!(out, BS, "kinetic", iA, T)
end

function ∇nuclear(BS::BasisSet, iA, T::DataType = Float64)
    # Pre allocate output
    out = zeros(T, BS.nbas, BS.nbas, 3)
    return ∇nuclear!(out, BS, iA, T)
end

function ∇nuclear!(out, BS::BasisSet, iA, T::DataType = Float64)

    if size(out) != (BS.nbas, BS.nbas, 3)
        throw(DimensionMismatch("Size of the output array needs to be (nbas, nbas, 3)"))
    end

    A = BS.atoms[iA]

    # Fudge lc_atoms
    lc_atoms_A = deepcopy(BS.lc_atoms)
    lc_atoms_woA = deepcopy(BS.lc_atoms)
    for i = eachindex(BS.atoms)
        if i == iA
            lc_atoms_woA[1 + 6*(i-1)] = 0
            continue
        end
        lc_atoms_A[1 + 6*(i-1)] = 0
    end

    # Get shell index (s0) where basis of the desired atom start
    s0 = 0
    for a in eachindex(BS.atoms)
        if BS.atoms[a] != A
            s0 += length(BS.basis[a])
        else
            break
        end
    end

    # Shell indexes for basis in the atom A
    Ashells = [(s0 + i -1) for i = 1:length(BS.basis[iA])]
    notAshells = Int[]
    for i = 1:BS.nshells
        # i-1 is used to convert to C notation (starting with 0)
        if !((i-1) in Ashells)
            push!(notAshells, i-1)
        end
    end

    lvals = [Libcint.CINTcgtos_spheric(i-1, BS.lc_bas) for i = 1:BS.nshells]
    ao_offset = [sum(lvals[1:(i-1)]) for i = 1:BS.nshells]
    Lmax = maximum(lvals)
    buf_arrays = [zeros(Cdouble, 3*Lmax^2) for _ = 1:Threads.nthreads()]

    # i ∉ A & j ∉ A
    @sync for i in notAshells
        Threads.@spawn begin
        @inbounds begin
            Li = lvals[i+1]
            ioff = ao_offset[i+1]
            I = (ioff+1):(ioff+Li)
            buf = buf_arrays[Threads.threadid()]
            for j in notAshells
                Lj = lvals[j+1]
                Lij = Li*Lj
                joff = ao_offset[j+1]
                J = (joff+1):(joff+Lj)

                # + ⟨i'|Va|j⟩ + ⟨i|Va|j'⟩   (Note that Va is the potential of the nuclei A alone!!)
                cint1e_ipnuc_sph!(buf, Cint.([i,j]), lc_atoms_A, BS.natoms, BS.lc_bas, BS.nbas, BS.lc_env)

                # Get strides for each cartesian
                for k in 1:3
                    r = (1+Lij*(k-1)):(k*Lij)
                    ∇k = reshape(buf[r], Int(Li), Int(Lj))
                    out[I,J,k] .+= ∇k  # ⟨i'|Va|j ⟩
                end
            end
        end # inbounds
        end # #spawn
    end # sync

    # i ∈ A & j ∈ A
    @sync for i in Ashells
        Threads.@spawn begin
        @inbounds begin
            Li = lvals[i+1]
            ioff = ao_offset[i+1]
            I = (ioff+1):(ioff+Li)
            buf = buf_arrays[Threads.threadid()]
            for j in Ashells
                Lj = lvals[j+1]
                Lij = Li*Lj
                joff = ao_offset[j+1]
                J = (joff+1):(joff+Lj)

                # - ⟨i'|∑Vc|j⟩ - ⟨i|∑Vc|j'⟩ c != a
                cint1e_ipnuc_sph!(buf, Cint.([i,j]), lc_atoms_woA, BS.natoms, BS.lc_bas, BS.nbas, BS.lc_env)

                for k in 1:3
                    r = (1+Lij*(k-1)):(k*Lij)
                    ∇k = reshape(buf[r], Int(Li), Int(Lj))
                    out[I,J,k] .-= ∇k  # ⟨i'|∑Vc|j ⟩ c != a
                end
            end
        end #inbounds
        end #spawn
    end #sync

    # i ∈ A & j ∉ A
    @sync for i in Ashells
        Threads.@spawn begin
        @inbounds begin
            Li = lvals[i+1]
            ioff = ao_offset[i+1]
            I = (ioff+1):(ioff+Li)
            buf = buf_arrays[Threads.threadid()]
            for j in notAshells
                Lj = lvals[j+1]
                Lij = Li*Lj
                joff = ao_offset[j+1]
                J = (joff+1):(joff+Lj)

                # - ⟨i'|∑Vc|j⟩ + ⟨i|Va|j'⟩ c != a
                cint1e_ipnuc_sph!(buf, Cint.([i,j]), lc_atoms_woA, BS.natoms, BS.lc_bas, BS.nbas, BS.lc_env)
                for k in 1:3
                    r = (1+Lij*(k-1)):(k*Lij)
                    ∇k = buf[r]
                    out[I,J,k] .-= reshape(∇k, Int(Li), Int(Lj))
                end

                cint1e_ipnuc_sph!(buf, Cint.([j,i]), lc_atoms_A, BS.natoms, BS.lc_bas, BS.nbas, BS.lc_env)
                for k in 1:3
                    r = (1+Lij*(k-1)):(k*Lij)
                    ∇k = buf[r]
                    out[I,J,k] .+= transpose(reshape(∇k, Int(Lj), Int(Li)))
                end
            end
        end #inbounds
        end #spawn
    end #sync

    # Add transpose values
    # This must be done outside the threaded loops
    # to avoid race conditions.
    for k in 1:3
        out[:,:,k] .+= out[:,:,k]'
    end

    return out
end


text:
for (ityp,jtyp) ∈ [("i8", UInt8), ("i16", UInt16), ("i32", UInt32), ("i64", UInt64), ("i128", UInt128)]
    @eval begin
        @inline function _atomic_load(ptr::Ptr{$jtyp})
            Base.llvmcall($("""
              %p = inttoptr i$(8sizeof(Int)) %0 to $(ityp)*
              %v = load atomic $(ityp), $(ityp)* %p acquire, align $(Base.gc_alignment(jtyp))
              ret $(ityp) %v
            """), $jtyp, Tuple{Ptr{$jtyp}}, ptr)
        end
        @inline function _atomic_store!(ptr::Ptr{$jtyp}, x::$jtyp)
            Base.llvmcall($("""
              %p = inttoptr i$(8sizeof(Int)) %0 to $(ityp)*
              store atomic $(ityp) %1, $(ityp)* %p release, align $(Base.gc_alignment(jtyp))
              ret void
            """), Cvoid, Tuple{Ptr{$jtyp}, $jtyp}, ptr, x)
        end
        @inline function _atomic_cas_cmp!(ptr::Ptr{$jtyp}, cmp::$jtyp, newval::$jtyp)
            Base.llvmcall($("""
              %p = inttoptr i$(8sizeof(Int)) %0 to $(ityp)*
              %c = cmpxchg $(ityp)* %p, $(ityp) %1, $(ityp) %2 acq_rel acquire
              %bit = extractvalue { $ityp, i1 } %c, 1
              %bool = zext i1 %bit to i8
              ret i8 %bool
            """), Bool, Tuple{Ptr{$jtyp}, $jtyp, $jtyp}, ptr, cmp, newval)
        end
    end
end
for op ∈ ["xchg", "add", "sub", "and", "nand", "or", "xor", "max", "min", "umax", "umin"] # "fadd", "fsub"
    f = Symbol("_atomic_", op, '!')
    for (ityp,jtyp) ∈ [("i8", UInt8), ("i16", UInt16), ("i32", UInt32), ("i64", UInt64), ("i128", UInt128)]
        @eval begin
            @inline function $f(ptr::Ptr{$jtyp}, x::$jtyp)
                Base.llvmcall($("""
                  %p = inttoptr i$(8sizeof(Int)) %0 to $(ityp)*
                  %v = atomicrmw $op $(ityp)* %p, $(ityp) %1 acq_rel
                  ret $(ityp) %v
                """), $jtyp, Tuple{Ptr{$jtyp}, $jtyp}, ptr, x)
            end
        end
    end
    @eval begin
        @inline function $f(ptr::Ptr{UInt}, x::ThreadState)
            reinterpret(ThreadState, $f(reinterpret(Ptr{UInt32}, ptr), reinterpret(UInt32, x)))
        end
    end
end
@inline _atomic_state(ptr::Ptr{UInt}) = reinterpret(ThreadState, _atomic_load(reinterpret(Ptr{UInt32}, ptr)))
@inline _atomic_store!(ptr::Ptr{UInt}, x::ThreadState) = _atomic_store!(reinterpret(Ptr{UInt32}, ptr), reinterpret(UInt32, x))
@inline function _atomic_cas_cmp!(ptr::Ptr{UInt}, cmp::ThreadState, newval::ThreadState)
    _atomic_cas_cmp!(reinterpret(Ptr{UInt32}, ptr), reinterpret(UInt32, cmp), reinterpret(UInt32, newval))
end


text:
using Test

using HomogeneousCoordinate

@testset "Smoke" begin
    p1 = Point(1.0, 0.0)
    p2 = Point(1.0, 1.0)
    @test isnothing(line_through_two_points(p1, p2))
    @test is_vertical(line_through_two_points(cart2homo(p1), cart2homo(p2)))

    p1 = Point(0.0, 1.0)
    p2 = Point(1.0, 3.0)
    p3 = Point(2.0, 5.0)
    @test are_collinear(p1, p2, p3)
    @test are_collinear(cart2homo(p1), cart2homo(p2), cart2homo(p3))
    @test is_point_on_line(p1, line_through_two_points(p2, p3))
    @test is_point_on_line(cart2homo(p1), line_through_two_points(cart2homo(p2), cart2homo(p3)))

    @test isinf(PointH([1.0, 1.0, 0.0]))
    @test isinf(LineH([0.0, 0.0, 1.0]))

    @test PointH([1.0, 1.0, 1.0]) == PointH([2.0, 2.0, 2.0])
    @test LineH([1.0, 1.0, 1.0]) == LineH([2.0, 2.0, 2.0])
end

@testset "Randomized" begin
    for _ in 1:50
        coords = Float64.(rand(1:5, (3, 2)))
        points = mapslices(v -> Point(v...), coords, dims=[2])
        h_points = cart2homo.(points)

        collinear = are_collinear(points[1], points[2], points[3])
        h_collinear = are_collinear(h_points[1], h_points[2], h_points[3])

        line_p1_p2 = line_through_two_points(points[1], points[2])
        h_line_p1_p2 = line_through_two_points(h_points[1], h_points[2])

        line_p2_p3 = line_through_two_points(points[2], points[3])
        h_line_p2_p3 = line_through_two_points(h_points[2], h_points[3])

        h_p1_clone = deepcopy(h_points[1])
        normalize!(h_p1_clone)
        @test h_p1_clone ≈ h_points[1]

        if isnothing(line_p1_p2)
            @test is_vertical(h_line_p1_p2)
        else
            @test cart2homo(line_p1_p2) ≈ h_line_p1_p2

            h_line_p1_p2_clone = deepcopy(h_line_p1_p2)
            normalize!(h_line_p1_p2_clone)
            @test h_line_p1_p2_clone ≈ h_line_p1_p2
        end

        if !isnothing(line_p1_p2) && !isnothing(line_p2_p3) && !collinear
            @test intersection_of_two_lines(line_p1_p2, line_p2_p3) ≈ points[2] &&
                  intersection_of_two_lines(h_line_p1_p2, h_line_p2_p3) ≈ h_points[2]
        end
    end
end

@testset "Dual" begin
    for _ in 1:5
        p = [PointH(rand(3)) for _ in 1:2]
        l = dual.(p)

        l_of_p = line_through_two_points(p...)
        p_of_l = intersection_of_two_lines(l...)

        @test l_of_p.v == p_of_l.v
    end
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
include("TwistedHessianEC.jl")
using .TwistedHessianEC

function partI(m::W, n:: W, P::S, Q::S, Pother::S, Qother::S, l::U, e::V, E::T) where {T <: EllipticCurve, S <: EllipticCurveElem, U <: Integer, V <: Integer, W <: Integer}
    K = P*m + Q*n
    Ei = E
    Pi = Pother
    Qi = Qother
    for i in V(1):e
        Ki = (l^(e-i))*K
        phi = Isogeny(Ei, Ki)
        K = phi(K)
        Pi = phi(Pi)
        Qi = phi(Qi)
        Ei = phi.codomain
    end
    return Ei, Pi, Qi
end

function partII(m::W, n:: W, P::S, Q::S, l::U, e::V, E::T) where {T <: EllipticCurve, S <: EllipticCurveElem, U <: Integer, V <: Integer, W <: Integer}
    K = P*m + Q*n
    Ei = E
    for i in V(1):e
        Ki = (l^(e-i))*K
        phi = Isogeny(Ei, Ki)
        K = phi(K)
        Ei = phi.codomain
    end
    return jInvariant(Ei)
end

# Example:
#Public params
la = big(2)
lb = big(3)
ea = 63
eb = 41
OGp = la^ea*lb^eb*11-1
F = Fpsquare(OGp)
E = TwiHesEC(F, one(F), F([3, 2458499195747164226075539892235780232099]))

XP = F([big(264972986297324027107883370923360079519), big(483679027370140600263719808627289108176)])
YP = F([big(2550992266319579452440794756588926041584), big(2214394007774501716861032109362868970653)])
Z = one(F)
P = E(XP, YP, Z)

XQ = F([big(1287714346072412371952478415842888698876), big(427368407186067938537508979187813002488)])
YQ = F([big(528075014001520129686916245492679497447), big(857779714716146821676003610752395199871)])
Q = E(XQ, YQ, Z)

Pa = normalized(P*(big(3)^41*11))
Qa = normalized(Q*(big(3)^41*11))

Pb = normalized(P*(big(2)^63*11))
Qb = normalized(Q*(big(2)^63*11))

# Alice and Bob chooses secret params
# They cannot both be divisible prime by "their" prime...
function genSecret(la, lb)
    ma = rand(big(1):OGp-1)
    na = rand(big(1):OGp-1)
    while ma%la == 0 && na%la == 0
        ma = rand(big(1):OGp-1)
        na = rand(big(1):OGp-1)
    end
    mb = rand(big(1):OGp-1)
    nb = rand(big(1):OGp-1)
    while mb%lb == 0 && nb%lb == 0
        mb = rand(big(1):OGp-1)
        nb = rand(big(1):OGp-1)
    end
    return ma, na, mb, nb
end

ma, na, mb, nb = genSecret(la, lb)

# Alice do her part
E_A, P_B, Q_B = partI(ma, na, Pa, Qa, Pb, Qb, la, ea, E)
# Bob does his part
E_B, P_A, Q_A = partI(mb, nb, Pb, Qb, Pa, Qa, lb, eb, E)

#Both derive the shared key, hopefully it is equal
key_a = partII(ma, na, P_A, Q_A, la, ea, E_B)
key_b = partII(mb, nb, P_B, Q_B, lb, eb, E_A)
println("Alice key: $key_a")
println("Bob's key: $key_b")

# Del


text:
#Written By Johnathan Bizzano
include("../src/Computation/Computation.jl")
using Main.Computation

f(a, b) = 3 * (a + b)

function test_kernels(device)
    a = alloc(device, collect(1:10))
    b = alloc(device, collect(2:11))
    broadcast!(f, a, a, b)
    return Array(a)
end

function test_device_kernels()
    println("Starting Device Kernel Operation Test....")
    true_array = [f(i, i + 1) for i in 1:10]
    for device in get_devices()
        println("\tChecking $device")
        @assert test_kernels(device) == true_array "Failed $device Kernel Operation Test"
        println("\tSuccesfully Completed $device")
    end
    println("Successfully Compeleted Device Kernel Operation Tes!")
end

test_device_kernels()


text:
using GPUArrays
using Minibatch

function posenc(timestep::Integer, channel::Integer, nchannels::Integer)
    if iseven(channel)
        return sin(timestep/(10000^(channel/nchannels)))
    else
        return cos(timestep/(10000^((channel-1)/nchannels)))
    end
end

function posenc(idx::CartesianIndex, nchannels::Integer)
    return posenc(idx[2], idx[1], nchannels)
end

function posenc!(A::GPUArray, state, nchannels::Integer)
    idx = @cartesianidx A state
    @inbounds A[idx] += posenc(idx, nchannels)
    return A
end

function posenc!(A::GPUArray)
    nchannels = size(A, 1)
    gpucall(posenc, A, (nchannels,))
    return A
end

function posenc!(A::AbstractArray)
    nchannels = size(A, 1)
    for idx in CartesianRange(size(A))
        # @inbounds
        A[idx] += posenc(idx, nchannels)
    end
    return A
end

function posenc!(B::MaskedBatch)
    posenc!(B.data)
    B.data = B.data .* B.mask # TODO in-place
    return B
end

function posenc!(A::Flux.TrackedArray)
    return posenc!(Flux.Tracker.data(A))
end


text:
"""
    struct MatpowerShunt <: AbstractNodeLabel

Structure descendant of AbstractNodeLabel

# Fields
- `Yshunt::Complex128`: Yshunt = Gs + im Bs
"""
struct MatpowerShunt <: AbstractNodeLabel
  Yshunt::Complex128 # coeffS = Gs + im*Bs
end

## 1. Variables creation
# function create_vars!(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}) where T <: MyType
#   return
# end

"""
  [cstrname, polynom, lb, ub] = Snodal(element::T, busid::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}) where T <: MatpowerShunt

Return the polynomial power contribution for shunt `elemid` at `busid`. Will be used to construct power balance constraints in polynomial problem.\n
Return expression of Sshunt depending on the voltage `["SHUNT", conj(shunt_value)|V|^2, 0, 0]`\n

"""
function Snodal(element::T, busid::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}) where T <: MatpowerShunt
  p = conj(element.Yshunt) * abs2(bus_vars[volt_name()])
  return ["SHUNT", p, 0, 0]
end

## 3. Constraints creation
# function constraint(element::T, busid::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}, scenario::String, OPFpbs::OPFProblems) where T <: MyType
  # return Dict{String, Constraint}()
# end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
include("./cmd.jl")
include("./Potentials.jl")

using Random
using Plots

# Bilinear problem:
f(x, y) = 2 * (x .- 0.00)' * (y .- 0.00) - sum((y .- 1.0).^2)
g(x, y) = -f(x, y)

m = 1
n = 1

α = 4.0

ψ_x = squared_distance(α)
ψ_y = squared_distance(α)
x = ones(m) * α
y = ones(n) * α

# ψ_x = shannon_entropy(α)
# ψ_y = shannon_entropy(α)
# x = zeros(m) * α
# y = zeros(n) * α

cmd = CMD(x, ψ_x, f, y, ψ_y, g)

max_iter = 100
xHist = zeros(m, max_iter + 1)
yHist = zeros(n, max_iter + 1)
for k = 1 : 100
  xHist[:, k], yHist[:, k] = iterate_primal(cmd)
  step!(cmd)

  cmd.x̄ .= max.(cmd.x̄, 0.0)
  cmd.ȳ .= max.(cmd.ȳ, 0.0)
end
xHist[:, end], yHist[:, end]= iterate_primal(cmd)

max_stars_repo_name: f-t-s/CMD

max_issues_repo_name: f-t-s/CMD

max_forks_repo_name: f-t-s/CMD


text:
# # Workspace
#
#md # [`notebook`](@__NBVIEWER_ROOT_URL__notebooks/two_problems.ipynb)
#
#using WaterWaves1D,ProgressMeter,Plots,FFTW,Statistics;gr()
#include("../src/dependencies.jl")
#
#include("../src/models/PseudoSpectral.jl")
#include("../src/models/Spectral.jl")
#include("../src/models/WaterWaves.jl")
#
#include("../src/solvers/Euler.jl")
#include("../src/solvers/EulerSymp.jl")

# if there is an error at this step, try commenting the lines above and commenting out the line below
#include("../src/dependencies.jl")

#---- parameters
param = ( μ  = .25,
			ϵ  = .1,
        	N  = 2^8, 	# number of collocation points
            L  = π,	# size of the mesh (-L,L)
            T  = 50,		# final time of computation
            dt = 0.001, # timestep
			ns = 50,   	# data stored every ns step
				);
#---- initial data
g(x) = exp.(-abs.(x).^4);
z(x) = 0 .*(x .+1).*exp.(-x.^2);
init = Init(g,z);

#---- building problems
model=WaterWaves(param)
solver1=Euler(model)
solver2=EulerSymp(model,Niter=5,implicit=1)
solver3=EulerSymp(model,Niter=5,implicit=2)
solver=RK4(model)

problems = Problem[]
push!(problems, Problem(model, init, param, solver=solver ))
push!(problems, Problem(model, init, param, solver=solver1 ))
push!(problems, Problem(model, init, param, solver=solver2 ))
push!(problems, Problem(model, init, param, solver=solver3 ))

#---- computation
for problem in problems
	@time solve!( problem )
end
#solve!(problems[3])
#---- visualization
#include("../src/Figures.jl")
#using Plots
#gr()
(η,v,x,t)=solution(problems[1])
(η2,v2,x2,t)=solution(problems[2])
(η3,v3,x3,t)=solution(problems[3])
(η4,v4,x4,t)=solution(problems[4])

using Plots;gr()
plot(x1,[η-η2])
plot(x1,[η-η2 η-η3])


text:
#
# "unstructured" (and therefore serial) PIPS-NLP interface
#
include("pips_serial_cfunc.jl")

module PipsNlpInterfaceSerial

using PipsNlpSolverSerial

using StructJuMP, JuMP
using StructJuMPSolverInterface

import MathProgBase

type NonStructJuMPModel <: ModelInterface
    model::JuMP.Model
    jac_I::Vector{Int}
    jac_J::Vector{Int}
    hess_I::Vector{Int}
    hess_J::Vector{Int}
    nz_jac::Vector{Int}
    nz_hess::Vector{Int}
    nzj::Int
    nzh::Int
    g_iter::Int

    write_solution::Function
    get_x0::Function
    numvars::Function
    numcons::Function
    nele_jac::Function
    nele_hess::Function
    get_bounds::Function
    eval_f::Function
    eval_g::Function
    eval_grad_f::Function
    eval_jac_g::Function
    eval_h::Function

    function NonStructJuMPModel(model)
        instance = new(model,
            Vector{Int}(), Vector{Int}(), Vector{Int}(), Vector{Int}(),
            Vector{Int}(), Vector{Int}(), 0 , 0, 0
            )

        instance.write_solution = function(x)
            @assert length(x) == getTotalNumVars(m)
            m = instance.model
            idx = 1
            for i = 0:num_scenarios(m)
                mm = getModel(m,i)
                for j = 1:getNumVars(m,i)
                    setvalue(Variable(mm,j),x[idx])
                    idx += 1
                end
            end
        end

        instance.get_x0 = function(x)
            @assert length(x) == getTotalNumVars(m)
            m = instance.model
            idx = 1
            for i = 0:num_scenarios(m)
                mm = getModel(m,i)
                for j = 1:getNumVars(m,i)
                    v_j = getvalue(Variable(mm,j))
                    x[idx] = isnan(v_j)? 1.0:v_j
                    idx += 1
                end
            end
            # @show x
            return x
        end
        instance.numvars = function()
            return getTotalNumVars(instance.model)
        end

        instance.numcons = function()
            return getTotalNumCons(instance.model)
        end

        instance.nele_jac = function()
            @assert length(instance.jac_I) == length(instance.jac_J)
            mat = sparse(instance.jac_I,instance.jac_J,ones(Float64,length(instance.jac_I)))
            instance.nzj = length(mat.nzval)
            return instance.nzj
        end

        instance.nele_hess = function()
            @assert length(instance.hess_I) == length(instance.hess_J)
            mat = sparse(instance.hess_J,instance.hess_I,ones(Float64,length(instance.hess_I)))
            instance.nzh = length(mat.nzval)
            return instance.nzh
        end

        instance.get_bounds = function()
            m = instance.model
            nvar = getTotalNumVars(m)
            ncon = getTotalNumCons(m)
            x_L = Vector{Float64}(nvar)
            x_U = Vector{Float64}(nvar)
            g_L = Vector{Float64}(ncon)
            g_U = Vector{Float64}(ncon)

            row_start = 1
            col_start = 1
            for i = 0:num_scenarios(m)
                mm = getModel(m,i)
                nx = getNumVars(m,i)
                array_copy(mm.colUpper, 1, x_U, col_start, nx)
                array_copy(mm.colLower, 1, x_L, col_start, nx)

                lb,ub = JuMP.constraintbounds(mm)
                ncons = getNumCons(m,i)
                array_copy(lb,1,g_L,row_start,ncons)
                array_copy(ub,1,g_U,row_start,ncons)

                row_start += getNumCons(m,i)
                col_start += getNumVars(m,i)
            end
            # @show g_L, g_U
            return x_L, x_U, g_L, g_U
        end

        instance.eval_f = function(x)
            m = instance.model
            obj = 0.0
            start_idx = 1
            for i=0:num_scenarios(m)
                x_new = strip_x(m,i,x,start_idx)
                obj += MathProgBase.eval_f(get_nlp_evaluator(m,i),x_new)
                start_idx += getNumVars(m,i)
            end
            # @printf("++++++++++++++++ eval_f  \n")
            # @show obj
            # @show x
            return obj;
        end

        instance.eval_g = function(x,g)
            m = instance.model
            @assert length(g) == getTotalNumCons(m)
            start_idx = 1
            g_start_idx = 1
            for i=0:num_scenarios(m)
                x_new = strip_x(instance.model,i,x,start_idx)
                ncon = getNumCons(m,i)
                g_new = Vector{Float64}(ncon)
                e = get_nlp_evaluator(m,i)
                MathProgBase.eval_g(e,g_new,strip_x(instance.model,i,x,start_idx))
                array_copy(g_new,1,g,g_start_idx,ncon)
                g_start_idx += ncon
                start_idx += getNumVars(m,i)
            end
            # @printf("+++++++++++++ eval_g \n")
            # @show x , g
        end

        instance.eval_grad_f = function(x,grad_f)
            fill!(grad_f,0.0)
            m = instance.model
            start_idx = 1
            for i=0:num_scenarios(m)
                x_new = strip_x(instance.model,i,x,start_idx)
                e = get_nlp_evaluator(m,i)

                g_f = Vector{Float64}(length(x_new))
                MathProgBase.eval_grad_f(e,g_f,x_new)
                nx = getNumVars(m,i)

                array_copy(g_f,1,grad_f,start_idx,nx)

                othermap = getStructure(getModel(m,i)).othermap
                for i in othermap
                    pid = i[1].col
                    cid = i[2].col
                    grad_f[pid] += g_f[cid]
                    @assert pid <= getNumVars(m,0)
                end
                start_idx += nx
            end
            # @printf("+++++++++++ eval_grad_f \n")
            # @show  x, grad_f
        end

        instance.eval_jac_g = function(x,mode,rows,cols,nzvals) #x, mode, irows, kcols, values)
            # @printf("+++++++++++ eval_jac_g %s \n", mode)
            m = instance.model
            if mode==:Structure
                mat = sparse(instance.jac_I,instance.jac_J,ones(Float64,length(instance.jac_I)))
                @assert length(mat.nzval) == instance.nzj
                array_copy(mat.rowval,1,rows,1,length(mat.rowval))
                array_copy(mat.colptr,1,cols,1,length(mat.colptr))
            else
                start_idx = 1
                value_start = 1
                value = Vector{Float64}(length(instance.jac_I))
                for i = 0:num_scenarios(m)
                    e = get_nlp_evaluator(m,i)
                    x_new = strip_x(instance.model,i,x,start_idx)
                    i_nz_jac = instance.nz_jac[i+1]
                    jac_g = Vector{Float64}(i_nz_jac)
                    MathProgBase.eval_jac_g(e,jac_g,x_new)
                    array_copy(jac_g,1,value,value_start,i_nz_jac)
                    nx = getNumVars(m,i)
                    start_idx += nx
                    value_start += i_nz_jac
                end
                @assert length(instance.jac_I) == length(instance.jac_J) == length(value)
                mat = sparse(instance.jac_I,instance.jac_J,value, getTotalNumCons(instance.model), getTotalNumVars(instance.model), keepzeros=true)

                jac_I = instance.jac_I
                jac_J = instance.jac_J
                # @printf( "m=%d; n=%d; \n", getTotalNumCons(instance.model), getTotalNumVars(instance.model))
                # @show jac_I, jac_J, value
                # @printf( "sjac=sparse(jac_I,jac_J,value,m,n); \n")

                @assert length(mat.nzval) == instance.nzj
                array_copy(mat.nzval,1,nzvals,1,instance.nzj)
            end
        end

        instance.eval_h = function(x, mode, rows, cols, obj_factor, lambda, nzvals) #x, mode, irows, kcols, obj_factor, lambda, values)
            # @printf("+++++++++++ eval_h - %s, %f \n", mode, obj_factor)
            # @show x
            # @show lambda
            m = instance.model
            if mode == :Structure
                mat = sparse(instance.hess_J,instance.hess_I,ones(Float64,length(instance.hess_I)))
                # @show mat
                @assert length(mat.nzval) == instance.nzh
                array_copy(mat.rowval,1,rows,1,length(mat.rowval))
                array_copy(mat.colptr,1,cols,1,length(mat.colptr))
            else
                start_idx = 1
                value_start = 1
                lambda_start = 1
                value = Vector{Float64}(length(instance.hess_I))
                for i = 0:num_scenarios(m)
                    e = get_nlp_evaluator(m,i)
                    x_new = strip_x(instance.model,i,x,start_idx)
                    nc = getNumCons(m,i)
                    lambda_new = Vector{Float64}(nc)
                    array_copy(lambda,lambda_start, lambda_new, 1, nc)
                    i_nz_hess = instance.nz_hess[i+1]
                    h = Vector{Float64}(i_nz_hess)
                    # @show x
                    # @show x_new
                    # @show obj_factor
                    # @show lambda_new
                    MathProgBase.eval_hesslag(e,h,x_new,obj_factor,lambda_new)
                    # @show h
                    array_copy(h,1,value,value_start,i_nz_hess)
                    nx = getNumVars(m,i)
                    start_idx += nx
                    lambda_start += nc
                    value_start += i_nz_hess
                end
                # @show value
                @assert length(instance.hess_I) == length(instance.hess_J) == length(value)
                mat = sparse(instance.hess_J,instance.hess_I,value,  getTotalNumVars(instance.model), getTotalNumVars(instance.model), keepzeros=true)
                @assert length(mat.nzval) == instance.nzh
                array_copy(mat.nzval,1,nzvals,1,instance.nzh)

                hess_I = instance.hess_J
                hess_J = instance.hess_I
                # @printf("m=%d; n=%d; \n", getTotalNumVars(instance.model), getTotalNumVars(instance.model))
                # @show hess_I, hess_J, value
                # @printf("shess=sparse(hess_I,hess_J,value,m,n); \n")

                # write_x("pips",instance.g_iter,x)
                instance.g_iter += 1
                # @show nzvals
            end
        end

        # initialization  jac
        col_offset = 0
        row_offset = 0
        m = instance.model
        for i = 0:num_scenarios(m)
            reverse_map = Dict{Int,Int}()
            mm = getModel(m,i)
            for ety in getStructure(mm).othermap
                reverse_map[ety[2].col] = ety[1].col #child->parent
            end
            # @show reverse_map
            e = get_nlp_evaluator(m,i)
            # @show "after e"
            i_jac_I, i_jac_J =  MathProgBase.jac_structure(e)
            # @show "after strct jac"
            for idx = 1:length(i_jac_J)
                jj = i_jac_J[idx]
                if haskey(reverse_map,jj)
                    push!(instance.jac_J, reverse_map[jj])
                else
                    push!(instance.jac_J, jj + col_offset)
                end
                push!(instance.jac_I, i_jac_I[idx] + row_offset)
            end
            push!(instance.nz_jac, length(i_jac_J)) #offset by 1

            col_offset += getNumVars(m,i)
            row_offset += getNumCons(m,i)
        end

        #initialization hess
        offset = 0
        for i = 0:num_scenarios(m)
            reverse_map = Dict{Int,Int}()
            mm = getModel(m,i)
            for ety in getStructure(mm).othermap
                reverse_map[ety[2].col] = ety[1].col #child->parent
            end

            e = get_nlp_evaluator(m,i)
            i_hess_I, i_hess_J =  MathProgBase.hesslag_structure(e)
            for idx = 1:length(i_hess_I)
                ii = i_hess_I[idx]
                jj = i_hess_J[idx]

                if haskey(reverse_map,ii)
                    new_ii = reverse_map[ii]
                else
                    new_ii = ii + offset
                end

                if haskey(reverse_map,jj)
                    new_jj = reverse_map[jj]
                else
                    new_jj = jj + offset
                end

                if(new_ii>new_jj)
                    push!(instance.hess_I,new_ii)
                    push!(instance.hess_J,new_jj)
                else
                    push!(instance.hess_I,new_jj)
                    push!(instance.hess_J,new_ii)
                end
            end
            push!(instance.nz_hess, length(i_hess_I)) #offset by 1

            offset += getNumVars(m,i)
        end
        return instance
    end
end

function structJuMPSolve(model; suppress_warmings=false,kwargs...)
    # @show typeof(model)
    nm = NonStructJuMPModel(model)
    x_L, x_U, g_L, g_U = nm.get_bounds()
    n = getTotalNumVars(model)
    m = getTotalNumCons(model)
    nele_jac = nm.nele_jac()
    nele_hess = nm.nele_hess()

    # @show x_L, x_U
    # @show g_L, g_U
    # @show n,m
    # @show nele_jac,nele_hess

    prob = createProblem(n, m, x_L, x_U, g_L, g_U, nele_jac, nele_hess,
                         nm.eval_f, nm.eval_g, nm.eval_grad_f, nm.eval_jac_g, nm.eval_h)
    # setProblemScaling(prob,1.0)
    nm.get_x0(prob.x)
    status = solveProblem(prob)
    nm.write_solution(prob.x)

    return PIPSRetCodeToSolverInterfaceCode[status]
end

KnownSolvers["PipsNlpSerial"] = PipsNlpInterfaceSerial.structJuMPSolve

end


text:
using BinDeps

@BinDeps.setup

mpath = get(ENV, "MAGICK_HOME", "") # If MAGICK_HOME is defined, add to library search path
if !isempty(mpath)
    push!(DL_LOAD_PATH, mpath)
    push!(DL_LOAD_PATH, joinpath(mpath,"lib"))
end
libnames = ["libMagickWand", "CORE_RL_wand_"]
suffixes = ["", "-Q16", "-6.Q16", "-Q8"]
options = ["", "HDRI"]
extensions = ["", ".so.4", ".so.5"]
aliases = vec(libnames.*transpose(suffixes).*reshape(options,(1,1,length(options))).*reshape(extensions,(1,1,1,length(extensions))))
libwand = library_dependency("libwand", aliases = aliases)

@linux_only begin
    provides(AptGet, "libmagickwand4", libwand)
    provides(AptGet, "libmagickwand5", libwand)
    provides(Pacman, "imagemagick", libwand)
    provides(Yum, "ImageMagick", libwand)
end

# TODO: remove me when upstream is fixed
@windows_only push!(BinDeps.defaults, BuildProcess)

@windows_only begin
    const OS_ARCH = (WORD_SIZE == 64) ? "x64" : "x86"

    # TODO: checksums: we have gpg
    # Extract the appropriate filename to download
    magick_base = "http://www.imagemagick.org/download/binaries"
    binariesfn = download(magick_base)
    str = readall(binariesfn)
    pattern = "ImageMagick-6.9.*?-Q16-$(OS_ARCH)-dll.exe"
    m = match(Regex(pattern), str)
    magick_exe = convert(ASCIIString, m.match)

    magick_tmpdir = BinDeps.downloadsdir(libwand)
    magick_url = "$(magick_base)/$(magick_exe)"
    magick_libdir = joinpath(BinDeps.libdir(libwand), OS_ARCH)

    innounp_url = "https://bintray.com/artifact/download/julialang/generic/innounp.exe"

    provides(BuildProcess,
        (@build_steps begin
            CreateDirectory(magick_tmpdir)
            CreateDirectory(magick_libdir)
            FileDownloader(magick_url, joinpath(magick_tmpdir, magick_exe))
            FileDownloader(innounp_url, joinpath(magick_tmpdir, "innounp.exe"))
            @build_steps begin
                ChangeDirectory(magick_tmpdir)
                info("Installing ImageMagick library")
                `innounp.exe -q -y -b -e -x -d$(magick_libdir) $(magick_exe)`
            end
        end),
        libwand,
        os = :Windows,
        unpacked_dir = magick_libdir,
        preload =
            """
            ENV["MAGICK_CONFIGURE_PATH"] = \"$(escape_string(magick_libdir))\"
            ENV["MAGICK_CODER_MODULE_PATH"] = \"$(escape_string(magick_libdir))\"
            """)
end

@osx_only begin
    if Pkg.installed("Homebrew") === nothing
            error("Homebrew package not installed, please run Pkg.add(\"Homebrew\")")
    end
    using Homebrew
    provides( Homebrew.HB, "imagemagick", libwand, os = :Darwin, onload =
    """
    function __init__()
        ENV["MAGICK_CONFIGURE_PATH"] = joinpath("$(Homebrew.prefix("imagemagick"))","lib","ImageMagick","config-Q16")
        ENV["MAGICK_CODER_MODULE_PATH"] = joinpath("$(Homebrew.prefix("imagemagick"))", "lib","ImageMagick","modules-Q16","coders")
        ENV["PATH"] = joinpath("$(Homebrew.prefix("imagemagick"))", "bin") * ":" * ENV["PATH"]
    end
    """ )
end

@BinDeps.install Dict([(:libwand, :libwand)])

# Save the library version; by checking this now, we avoid a runtime dependency on libwand
# See https://github.com/timholy/Images.jl/issues/184#issuecomment-55643225
module CheckVersion
include("deps.jl")
if isdefined(:__init__)
    __init__()
end
p = ccall((:MagickQueryConfigureOption, libwand), Ptr{Uint8}, (Ptr{Uint8},), "LIB_VERSION_NUMBER")
vstr = string("v\"", join(split(bytestring(p), ',')[1:3], '.'), "\"")
open(joinpath(dirname(@__FILE__),"versioninfo.jl"), "w") do file
    write(file, "const libversion = $vstr\n")
end
end


text:
using jInv.Mesh
using Multigrid
using SparseArrays
using LinearAlgebra
using Test
println("************************************************* GMG RAP for Poisson 2D ******************************************************");

domain = [0.0, 1.0, 0.0, 1.0];
n = [128,128];
Mr = getRegularMesh(domain,n)
G = getNodalGradientMatrix(Mr);
Ar = G'*G;
Ar = Ar + 1e-4*opnorm(Ar,1)*sparse(1.0I,size(Ar,2),size(Ar,2));

levels      = 4;
numCores 	= 8;
maxIter     = 5;
relativeTol = 1e-10;
relaxType   = "Jac-GMRES";
relaxParam  = 0.75;
relaxPre 	= 1;
relaxPost   = 1;
cycleType   ='V';
coarseSolveType = "NoMUMPS";

MG = getMGparam(Float64,Int64,levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType,0.5,0.0);

N = size(Ar,2);

b = Ar*rand(N,2);
b = b/norm(b);
x = zeros(N,2);

MGsetup(Ar,Mr,MG,size(b,2),true);

println("****************************** Stand-alone GMG RAP for Poisson: ******************************")
x = solveMG(MG,b,x,true)[1];
println("*******************Outside: ",norm(Ar*x - b)/norm(b))
@test norm(Ar*x - b) < 0.005;
println("****************************** Stand-alone GMG RAP : iterative coarsest ***********************")
coarseSolveType = "GMRES"
MG = getMGparam(Float64,Int64,levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType,0.5,0.0);
MGsetup(Ar,Mr,MG,size(b,2),true);
solveMG(MG,b,x,true);
@test norm(Ar*x - b) < 0.001;
println("****************************** GMRES preconditioned with GMG, Jac-GMRES relaxation: ******************************")
x[:] .= 0.0
relaxType   = "Jac-GMRES";
coarseSolveType = "NoMUMPS"
MG = getMGparam(Float64,Int64,levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType,0.5,0.0);
MGsetup(Ar,Mr,MG,size(b,2),true);
x = solveGMRES_MG(Ar,MG,b,x,true,10,true)[1];
println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Outside: ",norm(Ar*x - b)/norm(b))
@test norm(Ar*x - b) < 0.001;

println("****************************** 3D Poisson ******************************")

domain = [0.0, 1.0, 0.0, 1.0,0.0,1.0];
n = [32,32,16];
Mr = getRegularMesh(domain,n)
G = getNodalGradientMatrix(Mr);
Ar = G'*G;
Ar = Ar + 1e-4*opnorm(Ar,1)*sparse(1.0I,size(Ar,2),size(Ar,2));

MG = getMGparam(Float64,Int64,levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType,0.5,0.0);

N = size(Ar,2);

b = Ar*rand(N,2);
b = b/norm(b);
x = zeros(N,2);

MGsetup(Ar,Mr,MG,size(b,2),true);

println("****************************** Stand-alone 3D GMG RAP Poisson: ******************************")
solveMG(MG,b,x,true);
@test norm(Ar*x - b) < 0.01;

Ar = 0;
b = 0;
x = 0;
Mr = 0;
copySolver(MG);
destroyCoarsestLU(MG);
Multigrid.clear!(MG)
MG = 0;