text:
module Dubstep

using Cassette
using LinearAlgebra
export construct, TracedRun, trace, TraceCtx, LPCtx, replacenorm,
    GraftCtx, replacefunc

function construct(T::Type, args...)
    @info "constructing a model $T"
    return T(args...)
end

Cassette.@context TraceCtx

"""    TraceCtx

builds dynamic analysis traces of a model for information extraction
"""
TraceCtx

function Cassette.overdub(ctx::TraceCtx, args...)
    subtrace = Any[]
    push!(ctx.metadata, args => subtrace)
    if Cassette.canrecurse(ctx, args...)
        newctx = Cassette.similarcontext(ctx, metadata = subtrace)
        return Cassette.recurse(newctx, args...)
    else
        return Cassette.fallback(ctx, args...)
    end
end

function Cassette.overdub(ctx::TraceCtx, f::typeof(Base.vect), args...)
    @info "constructing a vector"
    push!(ctx.metadata, (f, args))
    return Cassette.fallback(ctx, f, args...)
end

function Cassette.overdub(ctx::TraceCtx, f::typeof(Core.apply_type), args...)
    # @info "applying a type $(args)"
    push!(ctx.metadata, (f, args))
    return Cassette.fallback(ctx, f, args...)
end

# TODO: support calls like construct(T, a, f(b))
function Cassette.overdub(ctx::TraceCtx, f::typeof(construct), args...)
    @info "constructing with type $f"
    push!(ctx.metadata, (f, args))
    y = Cassette.fallback(ctx, f, args...)
    @info "constructed model: $y"
    return y
end

"""    TracedRun{T,V}

captures the dataflow of a code execution. We store the trace and the value.

see also `trace`.
"""
struct TracedRun{T,V}
    trace::T
    value::V
end

"""    trace(f)

run the function f and return a TracedRun containing the trace and the output.
"""
function trace(f::Function)
    trace = Any[]
    val = Cassette.recurse(TraceCtx(metadata=trace), f)
    return TracedRun(trace, val)
end

Cassette.@context LPCtx

"""    LPCtx

replaces all calls to `LinearAlgebra.norm` with a different `p`.

This context is useful for modifying statistical codes or machine learning regularizers.
"""
LPCtx

function Cassette.overdub(ctx::LPCtx, args...)
    if Cassette.canrecurse(ctx, args...)
        newctx = Cassette.similarcontext(ctx, metadata = ctx.metadata)
        return Cassette.recurse(newctx, args...)
    else
        return Cassette.fallback(ctx, args...)
    end
end

function Cassette.overdub(ctx::LPCtx, f::typeof(norm), arg, power)
    p = get(ctx.metadata, power, power)
    return f(arg, p)
end

"""    replacenorm(f::Function, d::AbstractDict)

run f, but replace every call to norm using the mapping in d.
"""
function replacenorm(f::Function, d::AbstractDict)
    ctx = LPCtx(metadata=d)
    return Cassette.recurse(ctx, f)
end

Cassette.@context GraftCtx

"""    GraftCtx

grafts an expression from one simulation onto another

This context is useful for modifying simulations by changing out components to add features

see also: [`Dubstep.LPCtx`](@ref)
"""
GraftCtx

function Cassette.overdub(ctx::GraftCtx, f, args...)
    if Cassette.canrecurse(ctx, f, args...)
        newctx = Cassette.similarcontext(ctx, metadata = ctx.metadata)
        return Cassette.recurse(newctx, f, args...)
    else
        return Cassette.fallback(ctx, f, args...)
    end
end

"""    replacefunc(f::Function, d::AbstractDict)

run f, but replace every call to f using the context GraftCtx.
in order to change the behavior you overload overdub based on the context.
Metadata used to influence the context is stored in d.

see also: `bin/graft.jl` for an example.
"""
function replacefunc(f::Function, d::AbstractDict)
    ctx = GraftCtx(metadata=d)
    return Cassette.recurse(ctx, f)
end

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

# Cassette.@context TypeCtx

# """   TypeCtx

# creates a MetaDiGraph tracking the types of args and ret values throughout a script

# """
# TypeCtx

# """   FCollector(depth::Int,frame::function,data::FCollector)

# struct to collect all the "frames" called throughout a script

# """
# mutable struct FCollector{I,F,C}
#     depth::I
#     frame::F
#     data::Vector{C}
# end

# """   FCollector(depth::Int,frame::Frame)

# this is an initialization funtion for the FCollector

# """
# function FCollector(d::Int, f)
#     FCollector(d, f, FCollector[])
# end

# """ Frame(func, args, ret, subtrace)

# a structure to hold metadata for recursive type information for each function call
# Every frame can be thought of as a single stack frame when a function is called

# """
# mutable struct Frame{F,T,U}
#     func::F
#     args::T
#     ret::U
# end

# function Cassette.overdub(ctx::TypeCtx, f, args...) # add boilerplate for functionality
#     c = FCollector(ctx.metadata.depth-1, Frame(f, args, Any))
#     push!(ctx.metadata.data, c)
#     if c.depth > 0 && Cassette.canrecurse(ctx, f, args...)
#         newctx = Cassette.similarcontext(ctx, metadata = c)
#         z = Cassette.recurse(newctx, f, args...)
#         c.frame.ret = typeof(z)
#         return z
#     else
#         z = Cassette.fallback(ctx, f, args...)
#         c.frame.ret = typeof(z)
#         return z
#     end
# end

# Cassette.canrecurse(ctx::TypeCtx,::typeof(Base.vect), args...) = false # limit the stacktrace in terms of which to recurse on
# Cassette.canrecurse(ctx::TypeCtx,::typeof(FCollector)) = false
# Cassette.canrecurse(ctx::TypeCtx,::typeof(Frame)) = false

# """    buildgraph

# internal function used in the typegraphfrompath
# takes the collector object and returns a metagraph

# """
# function buildgraph(g,collector)
#     try
#         add_vertex!(g,:name,collector.frame.args)
#     catch
#         nothing
#     end
#     try
#         add_vertex!(g,:name,collector.frame.ret)
#     catch
#         nothing
#     end
#     try
#         add_edge!(g,g[collector.frame.args,:name],g[collector.frame.ret,:name],:name,collector.frame.func)
#     catch
#         nothing
#     end
#     for frame in collector.data
#         buildgraph(g,frame)
#     end
#     return g
# end

# """    typegraph(path::AbstractString,maxdepth::Int)

# This is a function that takes in an array of script and produces a MetaDiGraph descibing the system.
# takes in optional parameter of recursion depth on the stacktrace defaulted to 3

# """
# function typegraph(m::Module,maxdepth::Int=3)

#     extractor = FCollector(maxdepth, Frame(nothing, (), nothing,)) # init the collector object
#     ctx = TypeCtx(metadata = extractor);     # init the context we want
#     Cassette.overdub(ctx,m.main);    # run the script internally and build the extractor data structure
#     g = MetaDiGraph()    # crete a graph where we will init our tree
#     set_indexing_prop!(g,:name)    # we want to set this metagraph to be able to index by the names
#     return buildgraph(g,extractor)    # pass the collector ds to make the acutal metagraph

# end

end #module


text:
function _shoot(prob::CauchyProblem; i, itol)
    direction = monotonicity(prob)
    limit = i + direction*itol

    odesol = _integrate(prob, limit=limit)

    residual = direction*typemax(i)

    if odesol.retcode === :Terminated &&
                    direction*odesol.u[end][1] <= direction*limit
        residual = odesol.u[end][1] - i
    end

    return odesol, residual
end

function _shoot(prob::DirichletProblem; d_dϕb, itol)
    _shoot(CauchyProblem(prob.eq, b=prob.b, d_dϕb=d_dϕb, ϕb=prob.ϕb),
           i=prob.i, itol=itol)
end

function _shoot(prob::FlowrateProblem; b, itol, ϕbtol)
    if isindomain(prob.eq, b)
        ϕb = !iszero(prob.ϕb) ? prob.ϕb : ϕbtol

        d_dϕb = d_dϕ(prob, :b, b=b, ϕb=ϕb)

        return _shoot(CauchyProblem(prob.eq, b=b, d_dϕb=d_dϕb, ϕb=ϕb),
                      i=prob.i, itol=itol)
    end
    return nothing, -monotonicity(prob)*typemax(prob.i)
end

function solve(prob::DirichletProblem; d_dϕb_hint=nothing,
                                       itol=1e-3,
                                       maxiter=100)

    @argcheck itol≥zero(itol)
    @argcheck maxiter≥0

    @argcheck isindomain(prob.eq, prob.b) DomainError(prob.b, "prob.b not valid for the given equation")

    residual = prob.b - prob.i

    if abs(residual) ≤ itol
        odesol, _ = _shoot(prob, d_dϕb=zero(prob.b/prob.ϕb), itol=itol)
        return Solution(odesol, prob.eq, 0)
    end

    @argcheck isindomain(prob.eq, prob.i - monotonicity(prob)*itol) DomainError(prob.i, "prob.i not valid for the given equation and itol")

    if !isnothing(d_dϕb_hint)
        @argcheck sign(d_dϕb_hint) == monotonicity(prob) "sign of d_dϕb_hint must be consistent with initial and boundary conditions"
    else
        d_dϕb_hint = d_dϕ(prob, :b_hint)
    end

    d_dϕb_trial = bracket_bisect(zero(d_dϕb_hint), d_dϕb_hint, residual)

    for iterations in 1:maxiter
        odesol, residual = _shoot(prob, d_dϕb=d_dϕb_trial(residual), itol=itol)

        if abs(residual) ≤ itol
            return Solution(odesol, prob.eq, iterations)
        end
    end

    throw(SolvingError("failed to converge within $maxiter iterations"))
end

function solve(prob::FlowrateProblem; b_hint=nothing,
                                      itol=1e-3,
                                      ϕbtol=1e-6,
                                      maxiter=100)

    @argcheck itol≥zero(itol)
    if iszero(prob.ϕb)
        @argcheck ϕbtol>zero(ϕbtol)
    else
        @argcheck ϕbtol≥zero(ϕbtol)
    end
    @argcheck maxiter≥0

    if monotonicity(prob) == 0
        odesol, residual = _shoot(prob, b=prob.i, itol=itol, ϕbtol=ϕbtol)
        @assert iszero(residual)
        return Solution(odesol, prob, 0)
    end

    if !isnothing(b_hint)
        @argcheck sign(prob.i - b_hint) == monotonicity(prob) "sign of b_hint must be consistent with initial and boundary conditions"
    else
        b_hint = prob.i - oneunit(prob.i)*monotonicity(prob)
    end

    prob.i - oneunit(prob.i)*monotonicity(prob)
    b_trial = bracket_bisect(prob.i, b_hint)
    residual = nothing

    for iterations in 1:maxiter
        odesol, residual = _shoot(prob, b=b_trial(residual), itol=itol, ϕbtol=ϕbtol)

        if abs(residual) ≤ itol
            return Solution(odesol, prob.eq, iterations)
        end
    end

    throw(SolvingError("failed to converge within $maxiter iterations"))
end

export solve

ext: jl


text:
# This file is part of the Julia package ModularForms.jl
#
# Copyright (c) 2018-2019: Lynn Engelberts and Alexandru Ghitza.

"""
    poly_to_power_series(f, K, prec=10)

Convert the polynomial `f` over `K` to a relative power series
to precision `prec` as a polynomial.

The power series has the same coefficient ring and variable
name as f. The input polynomial f must be correct up to
precision prec.

# Arguments
' `f::PolyElem`: polynomial over K, correct up to prec
- `K`: base ring of the output
- `prec::Integer=10`: precision of the output

# Examples still missing
"""
function poly_to_power_series(f::PolyElem, K, prec::Int=10)

   #initialize
   varname = string(gen(parent(f)))
   R, q = PowerSeriesRing(K, prec, varname)

   if f == 0
      return R(0)
   end

   d = min(degree(f), prec-1)
   c = [coeff(f, i) for i in 0:d]
   power_series = R(c, d+1, prec, 0)

   return power_series
end

function poly_to_power_series(f::PolyElem, prec::Int=10)
   return poly_to_power_series(f, base_ring(f), prec)
end


text:
using ODEInterfaceDiffEq, DiffEqProblemLibrary, DiffEqBase
using Test

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear

prob = prob_ode_linear
sol =solve(prob,dopri5(),dt=1//2^(4))

sol =solve(prob,dopri5())
#plot(sol,plot_analytic=true)
sol =solve(prob,dopri5(),save_everystep=false)
@test sol.t == [0.0,1.0]

sol =solve(prob,dopri5(),saveat = 0.1)
@test sol.t == collect(0:0.1:1)

sol = solve(prob,dopri5(),save_on=false,save_start=false)
@test isempty(sol.t) && isempty(sol.u)


text:
struct ZeroMean{T<:Real} <: PriorMean{T} end

"""
    ZeroMean()

Construct a mean prior set to `0` and which cannot be updated.
"""
function ZeroMean()
    return ZeroMean{Float64}()
end

Base.show(io::IO, ::MIME"text/plain", ::ZeroMean) = print(io, "Zero Mean Prior")

update!(::ZeroMean{T}, ::Any) where {T<:Real} = nothing

(μ::ZeroMean{T})(x::Real) where {T<:Real} = zero(T)
(μ::ZeroMean{T})(x::AbstractMatrix) where {T} = zeros(T, size(x, 1))
(μ::ZeroMean{T})(x::AbstractVector) where {T} = zeros(T, length(x))


text:
function download_reactions_from_kegg(path_to_ec_array,path_to_reaction_file)

    # initalize -
    record_array = String[]

    # load the ec_array file -
    ec_array = readdlm(path_to_ec_array,',')

    # how many ec's do we have?
    (number_of_ec_numbers,nc) = size(ec_array)

    for ec_number_index = 1:number_of_ec_numbers

        # get the ec_number -
        rxn_index = ec_array[ec_number_index,1]
        ec_number = ec_array[ec_number_index,2]

        # use the KEGG API to download the reaction -
        ec_rn_linkage = readstring(`curl -X GET http://rest.kegg.jp/link/rn/$(ec_number)`)

        # if the server returned something, it should be in the form ec# rn#
        # we need the rn# to pull down the reaction string -
        if (length(ec_rn_linkage)>1)

            # ec_rn_linkage can have more than one association -
            # split around the \n, and then process each item
            P = reverse(split(chomp(ec_rn_linkage),'\n'))
            while (!isempty(P))

                # pop -
                local_record = pop!(P)

                # The second fragment is the rn# -
                # split around the \t
                rn_number = string(split(local_record,'\t')[2])

                # ok, run a second query to pull the reaction string down -
                reaction_string = chomp(readstring(`curl -X GET http://rest.kegg.jp/list/$(rn_number)`))

                # if we have something in reaction string, create a string record -
                if (length(reaction_string)!=0)
                    tmp_string = "$(rxn_index),$(ec_number),$(reaction_string)"
                    push!(record_array,tmp_string)
                end
            end
        end
    end

    # dump reaction array to disk -
    open(path_to_reaction_file, "w") do f

        for line_item in record_array
            write(f,"$(line_item)\n")
        end
    end
end

function download_ec_number_from_kegg(gene_location::String,kegg_organism_code::String)

    # TODO: Check for crazy strings -

    # get the sequence -
    ecdata = read(`curl -s -X GET http://rest.kegg.jp/link/ec/$(gene_location)`, String)

    # return this record -
    return ecdata
end

function download_gene_sequence_from_kegg(gene_location::String,kegg_organism_code::String)

    # get the sequence -
    ntseq = read(`curl -s -X GET http://rest.kegg.jp/get/$(gene_location)/ntseq`,String)

    # remove the header -
    P = split(ntseq,'\n')

    # return the sequence -
    return P[3:end-1]
end

function download_protein_sequence_from_kegg(gene_location::String,kegg_organism_code::String)

    # get the sequence -
    aaseq = read(`curl -s -X GET http://rest.kegg.jp/get/$(gene_location)/aaseq`,String)

    # remove the header -
    P = split(aaseq,'\n')

    # return the sequence -
    return P[3:end-1]
end

function write_ec_numbers_to_disk(path_to_input_gene_file::String,path_to_output_file::String,kegg_organism_code::String)

    # initalize -
    ec_number_buffer = String[]

    # get -
    gene_list = readdlm(path_to_input_gene_file,',')

    # how many genes do we have?
    (number_of_genes,number_of_cols) = size(gene_list)

    # create a progress meter -
    p = Progress(number_of_genes,color=:yellow)

    # start -
    @info "Downloading gene = ecnumber pairs from KEGG ..."

    # iterate through the list of tags -
    for gene_index = 1:number_of_genes

        # get the gene tag -
        gene_tag = String(gene_list[gene_index,2])

        # get the buffer -
        line_buffer = download_ec_number_from_kegg(gene_tag,kegg_organism_code)

        # check - does this have multiple records?
        P = split(line_buffer,'\n')
        number_of_fragments = length(P)
        for fragment_index = 1:number_of_fragments

            fragment_string = string(P[fragment_index])
            if (length(fragment_string)>0)

                # add line to overall buffer -
                push!(ec_number_buffer,fragment_string)

                # message -
                msg = "$(gene_tag) => $(fragment_string) ($(gene_index) of $(number_of_genes))"

                # update the progress bar -
                ProgressMeter.next!(p; showvalues = [(:status,msg)])
            end
        end
    end

    # dump to disk -
    open("$(path_to_output_file)/ec_numbers.dat", "w") do f

        for line_item in ec_number_buffer
            write(f,"$(line_item)\n")
        end
    end

    @info "Completed ...\r"

end

function write_protein_sequences_to_disk(path_to_input_gene_file::String,path_to_output_file::String,kegg_organism_code::String)

    # get -
    gene_list = readdlm(path_to_input_gene_file,',')

    # how many genes do we have?
    (number_of_genes, number_of_cols) = size(gene_list)

    # start -
    @info "Downloading protein sequences from KEGG ..."

    # create a progress meter -
    p = Progress(number_of_genes,color=:yellow)

    # iterate through the list of tags -
    for gene_index = 1:number_of_genes

        # get the gene tag -
        gene_tag = String(gene_list[gene_index,2])

        # get the buffer -
        buffer = download_protein_sequence_from_kegg(gene_tag,kegg_organism_code)

        # dump to disk -
        file_name = "$(path_to_output_file)/P_$(gene_tag).seq"
        writedlm(file_name,buffer)

        # update the progress bar -
        ProgressMeter.next!(p; showvalues = [(:gene_index,("$(gene_index) of $(number_of_genes)")), (:gene,gene_tag)])
    end

    @info "Completed ...\r"
end

function write_gene_sequences_to_disk(path_to_input_gene_file::String,path_to_output_file::String,kegg_organism_code::String)

    # get -
    gene_list = readdlm(path_to_input_gene_file,',')

    # how many genes do we have?
    (number_of_genes, number_of_cols) = size(gene_list)

    # start -
    @info "Downloading gene sequences from KEGG ..."

    # create a progress meter -
    p = Progress(number_of_genes,color=:yellow)

    # iterate through the list of tags -
    for gene_index in 1:number_of_genes

        # get the gene tag -
        gene_tag = String(gene_list[gene_index,2])

        # get the buffer -
        buffer = download_gene_sequence_from_kegg(gene_tag,kegg_organism_code)

        # dump to disk -
        file_name = "$(path_to_output_file)/$(gene_tag).seq"
        writedlm(file_name,buffer)

        # update the progress bar -
        ProgressMeter.next!(p; showvalues = [(:gene_index,("$(gene_index) of $(number_of_genes)")), (:gene,gene_tag)])

        # message -
        # msg = "$(gene_tag) => downloaded ($(gene_index) of $(number_of_genes))"
        # println(msg)
    end

    @info "Completed ...\r"
end

function write_gene_symbols_to_disk(cobra_dictionary, path_to_output_gene_file::String, organism_id::String)

    # get list of genes from the cobra_dictionary -
    gene_list = cobra_dictionary["genes"]

    # initialize buffer -
    gene_index = 1
    gene_name_buffer = String[]
    for gene_name in gene_list
        line = "$(gene_index),$(organism_id):$(gene_name)"
        push!(gene_name_buffer,line)
        gene_index = gene_index + 1
    end

    # dump to disk -
    open("$(path_to_output_gene_file)/gene_list.dat", "w") do f

        for line_item in gene_name_buffer
            write(f,"$(line_item)\n")
        end
    end

    # start -
    @info "Wrote gene symbols to $(path_to_output_gene_file)/gene_list.dat"
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libplist_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libplist")
JLLWrappers.@generate_main_file("libplist", UUID("7bef5416-8376-5e75-ab4d-144910ca2ef4"))
end  # module libplist_jll


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

mutable struct ApplicationGatewayBackendAddressPool <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id

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

const _property_map_ApplicationGatewayBackendAddressPool = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"))
const _property_types_ApplicationGatewayBackendAddressPool = Dict{Symbol,String}(Symbol("id")=>"String")
Base.propertynames(::Type{ ApplicationGatewayBackendAddressPool }) = collect(keys(_property_map_ApplicationGatewayBackendAddressPool))
Swagger.property_type(::Type{ ApplicationGatewayBackendAddressPool }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ApplicationGatewayBackendAddressPool[name]))}
Swagger.field_name(::Type{ ApplicationGatewayBackendAddressPool }, property_name::Symbol) =  _property_map_ApplicationGatewayBackendAddressPool[property_name]

function check_required(o::ApplicationGatewayBackendAddressPool)
    true
end

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


text:
# Various specialised operations using the Cholesky factorisation.

Xt_A_X(A::Cholesky, x::AbstractVector) = sum(abs2, A.U * x)
function Xt_A_X(A::Cholesky, X::AbstractMatrix)
    V = A.U * X
    return Symmetric(V'V)
end

Xt_A_Y(X::AVM, A::Cholesky, Y::AVM) = (A.U * X)' * (A.U * Y)

Xt_invA_X(A::Cholesky, x::AbstractVector) = sum(abs2, A.U' \ x)
function Xt_invA_X(A::Cholesky, X::AVM)
    V = A.U' \ X
    return Symmetric(V'V)
end

Xt_invA_Y(X::AVM, A::Cholesky, Y::AVM) = (A.U' \ X)' * (A.U' \ Y)

diag_At_A(A::AVM) = vec(sum(abs2.(A); dims=1))

tr_At_A(A::AVM) = sum(abs2, A)

function diag_At_B(A::AVM, B::AVM)
    @assert size(A) == size(B)
    return vec(sum(A .* B; dims=1))
end

diag_Xt_A_X(A::Cholesky, X::AVM) = diag_At_A(A.U * X)

function diag_Xt_A_Y(X::AVM, A::Cholesky, Y::AVM)
    @assert size(X) == size(Y)
    return diag_At_B(A.U * X, A.U * Y)
end

diag_Xt_invA_X(A::Cholesky, X::AVM) = diag_At_A(A.U' \ X)

function diag_Xt_invA_Y(X::AM, A::Cholesky, Y::AM)
    @assert size(X) == size(Y)
    return diag_At_B(A.U' \ X, A.U' \ Y)
end

function Xtinv_A_Xinv(A::Cholesky, X::Cholesky)
    @assert size(A) == size(X)
    C = A.U \ (X.U' \ A.U')
    return Symmetric(C * C')
end


text:
""" Parse relation expressions in Julia syntax into undirected wiring diagrams.
"""
module RelationalPrograms
export RelationDiagram, UntypedRelationDiagram, TypedRelationDiagram,
  @relation, parse_relation_diagram

using Compat
using MLStyle: @match

using ...CategoricalAlgebra.CSets, ...Present
using ...WiringDiagrams.UndirectedWiringDiagrams
using ...WiringDiagrams.MonoidalUndirectedWiringDiagrams:
  TheoryUntypedHypergraphDiagram, TheoryHypergraphDiagram

# Data structures
#################

@present TheoryRelationDiagram <: TheoryUntypedHypergraphDiagram begin
  VarName::Data
  variable::Attr(Junction, VarName)
end

@present TheoryTypedRelationDiagram <: TheoryHypergraphDiagram begin
  VarName::Data
  variable::Attr(Junction, VarName)
end

const RelationDiagram = AbstractACSetType(TheoryRelationDiagram)
const UntypedRelationDiagram = ACSetType(TheoryRelationDiagram,
  index=[:box, :junction, :outer_junction], unique_index=[:variable])
const TypedRelationDiagram = ACSetType(TheoryTypedRelationDiagram,
  index=[:box, :junction, :outer_junction], unique_index=[:variable])

@present TheoryNamedRelationDiagram <: TheoryRelationDiagram begin
  port_name::Attr(Port, VarName)
  outer_port_name::Attr(OuterPort, VarName)
end

@present TheoryTypedNamedRelationDiagram <: TheoryTypedRelationDiagram begin
  port_name::Attr(Port, VarName)
  outer_port_name::Attr(OuterPort, VarName)
end

const UntypedNamedRelationDiagram = ACSetType(TheoryNamedRelationDiagram,
  index=[:box, :junction, :outer_junction], unique_index=[:variable])
const TypedNamedRelationDiagram = ACSetType(TheoryTypedNamedRelationDiagram,
  index=[:box, :junction, :outer_junction], unique_index=[:variable])

function RelationDiagram{Name}(ports::Int; port_names=nothing) where {Name}
  if isnothing(port_names)
    return UntypedRelationDiagram{Name,Symbol}(ports)
  end
  d = UntypedNamedRelationDiagram{Name,Symbol}(ports)
  set_subpart!(d, :outer_port_name, port_names)
  return d
end

function RelationDiagram{Name}(ports::AbstractVector{T};
                               port_names=nothing) where {T,Name}
  if isnothing(port_names)
    return TypedRelationDiagram{T,Name,Symbol}(ports)
  end
  d = TypedNamedRelationDiagram{T,Name,Symbol}(ports)
  set_subpart!(d, :outer_port_name, port_names)
  return d
end

# Relation macro
################

""" Construct an undirected wiring diagram using relation notation.

Unlike the `@program` macro for directed wiring diagrams, this macro departs
from the usual semantics of the Julia programming language. Function calls with
n arguments are now interpreted as assertions that an n-ary relation holds at a
particular point. For example, the composition of binary relations R ⊆ X × Y and
S ⊆ Y × Z can be represented as an undirected wiring diagram by the macro call

```julia
@relation (x,z) where (x::X, y::Y, z::Z) begin
  R(x,y)
  S(y,z)
end
```

In general, the context in the `where` clause defines the set of junctions in
the diagram and variable sharing defines the wiring of ports to junctions.

The ports and junctions of the diagram may be typed or untyped, and the ports
may be named or unnamed. Thus four possible types of diagrams may be returned,
with the type determined by the form of relation header:

1. Untyped, unnamed: `@relation (x,z) where (x,y,z) ...`
2. Typed, unnamed: `@relation (x,z) where (x::X, y::Y, z::Z) ...`
3. Untyped, named: `@relation (out1=x, out2=z) where (x,y,z) ...`
4. Typed, named: `@relation (out=1, out2=z) where (x::X, y::Y, z::Z) ...`
"""
macro relation(exprs...)
  :(parse_relation_diagram($((QuoteNode(expr) for expr in exprs)...)))
end

""" Parse an undirected wiring diagram from a relation expression.

For more information, see the corresponding macro [`@relation`](@ref).
"""
function parse_relation_diagram(expr::Expr)
  @match expr begin
    Expr(:function, head, body) => parse_relation_diagram(head, body)
    Expr(:->, head, body) => parse_relation_diagram(head, body)
    _ => error("Not a function or lambda expression")
  end
end

function parse_relation_diagram(head::Expr, body::Expr)
  # Parse variables and their types from context.
  outer_expr, all_vars, all_types = @match head begin
    Expr(:where, expr, context) => (expr, parse_relation_context(context)...)
    _ => (head, nothing, nothing)
  end
  var_types = if isnothing(all_types) # Untyped case.
    vars -> length(vars)
  else # Typed case.
    var_type_map = Dict{Symbol,Symbol}(zip(all_vars, all_types))
    vars -> getindex.(Ref(var_type_map), vars)
  end

  # Create wiring diagram and add outer ports and junctions.
  _, outer_port_names, outer_vars = parse_relation_call(outer_expr)
  isnothing(all_vars) || outer_vars ⊆ all_vars ||
    error("One of variables $outer_vars is not declared in context $all_vars")
  d = RelationDiagram{Symbol}(var_types(outer_vars),
                              port_names=outer_port_names)
  if isnothing(all_vars)
    new_vars = unique(outer_vars)
    add_junctions!(d, var_types(new_vars), variable=new_vars)
  else
    add_junctions!(d, var_types(all_vars), variable=all_vars)
  end
  set_junction!(d, ports(d, outer=true),
                incident(d, outer_vars, :variable), outer=true)

  # Add box to diagram for each relation call.
  body = Base.remove_linenums!(body)
  for expr in body.args
    name, port_names, vars = parse_relation_call(expr)
    isnothing(all_vars) || vars ⊆ all_vars ||
      error("One of variables $vars is not declared in context $all_vars")
    box = add_box!(d, var_types(vars), name=name)
    if !isnothing(port_names)
      set_subpart!(d, ports(d, box), :port_name, port_names)
    end
    if isnothing(all_vars)
      new_vars = setdiff(unique(vars), d[:variable])
      add_junctions!(d, var_types(new_vars), variable=new_vars)
    end
    set_junction!(d, ports(d, box), incident(d, vars, :variable))
  end
  return d
end

function parse_relation_context(context)
  terms = @match context begin
    Expr(:tuple) => return (Symbol[], nothing)
    Expr(:tuple, terms...) => terms
    _ => error("Invalid syntax in relation context $context")
  end
  vars = map(terms) do term
    @match term begin
      Expr(:(::), var::Symbol, type::Symbol) => (var => type)
      var::Symbol => var
      _ => error("Invalid syntax in term $expr of context")
    end
  end
  if vars isa AbstractVector{Symbol}
    (vars, nothing)
  elseif vars isa AbstractVector{Pair{Symbol,Symbol}}
    (first.(vars), last.(vars))
  else
    error("Context $context mixes typed and untyped variables")
  end
end

function parse_relation_call(call)
  @match call begin
    Expr(:call, name::Symbol, Expr(:parameters, args)) =>
      (name, parse_relation_kw_args(args)...)
    Expr(:call, name::Symbol) => (name, nothing, Symbol[])
    Expr(:call, name::Symbol, args...) =>
      (name, parse_relation_inferred_args(args)...)

    Expr(:tuple, Expr(:parameters, args...)) =>
      (nothing, parse_relation_kw_args(args)...)
    Expr(:tuple) => (nothing, nothing, Symbol[])
    Expr(:tuple, args...) => (nothing, parse_relation_inferred_args(args)...)
    Expr(:(=), args...) => (nothing, parse_relation_inferred_args([call])...)

    _ => error("Invalid syntax in relation $call")
  end
end

function parse_relation_kw_args(args)
  args = map(args) do arg
    @match arg begin
      Expr(:kw, name::Symbol, var::Symbol) => (name => var)
      _ => error("Expected name as keyword argument")
    end
  end
  (first.(args), last.(args))
end

function parse_relation_inferred_args(args)
  @assert !isempty(args) # Need at least one argument to infer named/unnamed.
  args = map(args) do arg
    @match arg begin
      Expr(:kw, name::Symbol, var::Symbol) => (name => var)
      Expr(:(=), name::Symbol, var::Symbol) => (name => var)
      var::Symbol => var
      _ => error("Expected name as positional or keyword argument")
    end
  end
  if args isa AbstractVector{Symbol}
    (nothing, args)
  elseif args isa AbstractVector{Pair{Symbol,Symbol}}
    (first.(args), last.(args))
  else
    error("Relation mixes named and unnamed arguments $args")
  end
end

end


text:
## 07.jl
# テンプレートによる生成文

function weather_template(x, y, z)
  return string(x)*"時の"*y*"は"*string(z)
end

println(weather_template(12,"気温",22.4))


text:
@info("-------[Testing Complementarity/MPEC]------------------------------------------")

include("../examples/mpec/bard1.jl")
include("../examples/mpec/dempe.jl")
include("../examples/mpec/gnash1m.jl")
include("bard1_test.jl")


text:
module FFDpacking

include("objects.jl")

export load_data,ffdpacking,config_container,draw_container,run_pack

#reading input data
"""
    load_data(filecont::String,fileitem::String)

Load data for use first fit decrease function (ffd).
This function return two vector. The first one is
a vector of item type (ItemT) and the second is a
vector of container type (ContainerT).

# Examples
```julia-repl
julia> output=load_data("container.txt","item.txt")
       list_cont=output[2]
       list_item=output[1]
```
"""
function load_data(filecont::String,fileitem::String)
    item_txt=readdlm(fileitem)
    n=length(item_txt[:,1])
    v_item=Array{ItemT}(n)
    for i=1:n
        v_item[i]=ItemT(item_txt[i,1],item_txt[i,2],0,0)
    end
    container_txt=readdlm(filecont)
    m=length(container_txt[:,2])
    v_container=Array{ContainerT}(m)
    for i=1:m
        v_container[i]=ContainerT(container_txt[i,1],container_txt[i,2],[])
    end
    return v_item,v_container
end

#main function for packing
"""
    ffd(list_container,list_item)

Allocate the list_item in containers type find in list_container.

# Examples

```julia-repl
julia> run_pack("container.txt","item.txt")
```
"""
function ffdpacking(list_container,list_item)
    #creating an array of container
    v_container=Array{ContainerT}(50)
    #selecting the type of container
    container=copy(list_container[2])
    #setting the first element of v_container
    for i=1:10
    v_container[i]=copy(container)
    end
    ncontainer=1
    #setting auxiliary variables
    s_height=0.0
    s_length=0.0
    size_list=length(list_item)
    updated_list=copy(list_item)
    strip=updated_list[1].height
    s_height=updated_list[1].height
    #main loop for all elements
    while !isempty(updated_list)
        i=1
        removed_item=[]
        if s_height<=v_container[ncontainer].height
            while s_length<=v_container[ncontainer].length && i<=size_list
                if s_length + updated_list[i].length<=v_container[ncontainer].length
                    s_length+=updated_list[i].length
                    removed_item=push!(removed_item,i)
                    v_container[ncontainer].content=[v_container[ncontainer].content;updated_list[i]]
                end
                i+=1
            end
            deleteat!(updated_list,removed_item)
            size_list=length(updated_list)
            #display(updated_list)
            if isempty(updated_list)
                break
            end
            s_height+=updated_list[1].height
            s_length=0.0
        else
            println("Container $ncontainer :")
            display(v_container[ncontainer])
            println("New container was created ...")
            s_height=updated_list[1].height
            s_length=0.0
            ncontainer+=1
            size_list=length(updated_list)

        end

    end
    println("Container $ncontainer :")
    display(v_container[ncontainer])
    println("Configuring containers ... ")
    for i=1:ncontainer
        v_container[i]=config_container(v_container[i])
        draw_container(v_container[i],i);
    end
    display(v_container[1:ncontainer])
end

#configuring container assuming that the allocation was done (auxiliary function)
"""
    config_container(container)

Given a contaneir with a list of items in its content property, this
function generates the position of items in that container. More specifically,
this function changes posX and posY in each item of the list in content.

"""

function config_container(container)
    s_length=0.0
    s_height=0.0
    strip=container.content[1].height
    for item in container.content
        if (s_length+item.length)>container.length
            s_length=0.0
            s_height+=strip
            strip=item.height
        end
        item.posX=s_length
        item.posY=s_height
        s_length=s_length+item.length
    end
    return container
end

#sort the list of items (auxiliary function)
"""
    sortitems(list_item)

Basic function for ordering list item.

"""

function sortitems(list_item)
    n=length(list_item)
    for i=1:n
        for j=i+1:n
            if list_item[i].height<list_item[j].height
                swap=list_item[j]
                list_item[j]=list_item[i]
                list_item[i]=swap
            end
        end
    end
    return list_item
end

#packages for drawing solutions
using PyPlot
using PyCall
@pyimport matplotlib.patches as patch
#draw the configuration of one container (auxiliary function)
function draw_container(container,nc)
    cfig = figure();
    ax = cfig[:add_subplot](1,1,1);
    ax[:set_aspect]("equal");
    ax[:set_xlim]([0,container.length]);
    ax[:set_ylim]([0,container.height]);
    for item in container.content
        c = patch.Rectangle([item.posX,item.posY],item.length,item.height,facecolor="red");
        ax[:add_patch](c);
    end
    #ax[:add_artist](c)
    cfig[:savefig]("container$nc.png");
end

#function for testing the main program run_pack
#example run_pack("container.txt","item.txt")
function run_pack(SContainer::String,SItem::String)
    output=load_data(SContainer,SItem)
    list_cont=output[2] #The list of containers
    list_item=output[1] #The list of items
    list_item=sortitems(list_item)
    ffdpacking(list_cont,list_item)
end

end

max_stars_repo_name: evcastelani/FFDpacking

max_issues_repo_name: evcastelani/FFDpacking

max_forks_repo_name: evcastelani/FFDpacking


text:
using Flux
using Flux: @epochs
using Flux: params
using Flux: throttle
using Plots
using Distributions
using DataFrames
using CSV

@inline phi(x) = @. x/(1 + exp(-x))
@inline f(x) = 2*π^2 .*sum(cos.(π*x[i,:]') for i=1:2, dims=1)

#build a block
block_layer = SkipConnection(Chain(Dense(10,10, phi),Dense(10,10, phi)), +)

#construct network skeleton
M = Chain(Dense(2,10),
block_layer,
block_layer,
block_layer,
block_layer,
Dense(10,1))

#check number of params
params(M)

#inner loss function
function inner_loss(N, space)
    a = rand(Uniform(0,1),2,N)
    loss = 0
    l1 = [1, 0]
    l2 = [0, 1]

    @inbounds loss += (sum(0.5*((M(a .+ space.*l1) .- M(a))./space).^2)/N)[1]
    @inbounds loss += (sum(0.5*((M(a .+ space.*l2) .- M(a))./space).^2)/N)[1]
    @inbounds loss += (sum(0.5*(π^2 * (M(a)).^2))/N)[1]
    @inbounds loss += (sum(-f(a) .*M(a))/N)[1]

    return loss
end

function bc_loss(β, N = 1000)

    d = (rand(1)/100)
    point = zeros(2,1000)
    p = rand(1,100)
    point = hcat(hcat(zeros(250), rand(250))', hcat(ones(250), rand(250))', hcat(rand(250), zeros(250))', hcat(rand(250), ones(250))')
    loss = 0
    for i=1:2
        if i == 1
            l1 = [1, 0]
        else
            l1 = [0, 1]
        end
        @inbounds loss += β*(sum(((1/2*M(point[:,500*(i-1)+1:500*i] .+ d.*l1) .- 1/2*M(point[:,500*(i-1)+1:500*i] .- d.*l1))./d).^2)/100)[1]
    end
    return loss
end

#add losses
function total_loss(N, space = 0.01, β = 50)
    return inner_loss(N, space) + bc_loss(β)
end

function loss_true(xx=1,yy=1,xx1=1,xx2=2)
  points=rand(2,2000)
  F_true(x)= sum(cos.(pi*x[i,:]') for i=1:2, dims=1)
  errors=sqrt(sum((M(points)-F_true(points)).^2)/2000)
  return errors
end

#train model
function train_model(N = 100, train_time = 20)
    df = DataFrame(time = Int[], in_loss = Float64[],bound_loss = Float64[],true_loss = Float64[])
    c = 1
    evalcb = function()
        @show(total_loss(N), loss_true())
        in_loss = inner_loss(N, 0.01)
        bound_loss = bc_loss(50)
        true_loss = loss_true()
        push!(df, [c, in_loss, bound_loss, true_loss])
        c += 1
        if c == train_time
            CSV.write("$(c)_2d_neumann.csv",df)
        end
    end
    opt = ADAM()
    dataset = [(N) for i=1:train_time]
    Flux.train!(total_loss, params(M), zip(dataset), opt, cb=throttle(evalcb, 50))
end

@epochs 1000 train_model(4000)

#plot the solution over a mesh in the domain
F_true(x)= sum(cos.(pi*x[i,:]') for i=1:2, dims=1)
function plot_surface()
    data = zeros(100,100)
    grid_space = [i for i in 1:100]
    for i in grid_space
        for j in grid_space
            point = [i/100 j/100]
            data[i,j] = M(point')[1]
        end
    end

    gr()
    heatmap(1:size(data,1),
        1:size(data,2), data',
        c = cgrad([:blue, :white,:red, :yellow]),
        xlabel = "X", ylabel = "Y",
        title = "2D Neumann Fitted Solution")

    savefig("2d_neumann_fitted.png")

    for i in grid_space
        for j in grid_space
            point = [i/100 j/100]
            data[i,j] = F_true(point')[1]
        end
    end

    gr()
    heatmap(1:size(data,1),
        1:size(data,2), data',
        c = cgrad([:blue, :white,:red, :yellow]),
        xlabel = "X", ylabel = "Y",
        title = "2D Neumann True Solution")

    savefig("2d_neumann_true.png")
end

plot_surface()


text:
include("../Tests/Tests.jl")
include("utilities.jl")

using Test, JuMP

@testset "Term" begin
    include("term.jl")
end
@testset "Term fixed" begin
    include("term_fixed.jl")
end
@testset "Quadratic" begin
    include("quadratic.jl")
end
@testset "Choi" begin
    include("choi.jl")
end
@testset "Concave then convex cubic" begin
    include("concave_then_convex_cubic.jl")
end
@testset "Horn" begin
    include("horn.jl")
end
@testset "Lyapunov Switched System" begin
    include("lyapunov_switched_system.jl")
end
@testset "Motzkin" begin
    include("motzkin.jl")
end
@testset "BPT12e399" begin
    include("BPT12e399.jl")
end
@testset "Max Cut" begin
    include("maxcut.jl")
end
@testset "Chebyshev" begin
    include("chebyshev.jl")
end
@testset "Quartic comparison" begin
    include("quartic_comparison.jl")
end
@testset "SOSDEMO9" begin
    include("sosdemo9.jl")
end
@testset "SOSDEMO10" begin
    include("sosdemo10.jl")
end


text:
using DataStructures

import Base.Ordering
import Base.lt
import DataStructures.eq

include("file_blacklist.jl")
include("auxiliary_functions.jl")

#=
 # Lambda function used to add custom order for Dictionaires
=#
nav_order = Lt((a, b) ->
    begin
        custom_order = ["members", "papers", "philosophy", "teaching-resources", "thesis-art", "open-projects-and-positions", "contact-us"]
        a_in = a in custom_order
        b_in = b in custom_order
        if(a_in && b_in)
            return findfirst(isequal(a), custom_order) < findfirst(isequal(b), custom_order) ? true : false
        elseif(a_in && !b_in)
             return true
        elseif(!a_in && b_in)
            return false
        else
            return isless(lowercase(a), lowercase(b))
        end
    end)

#=
 # Appends all files and folders, subfolders and their contents of a given path into the provided
 # dict, but ignoring all black listed files and folders.
 # Every folder is pointed at another dict, containing all subfiles and -folders.
=#
function append_files!(dict, folder)
    try
        for file in readdir(string("./", folder))

            # Check for .md file and add to dict, pointing to empty dict
            if contains(file, ".md") && !file_blacklisted(file)

                push!(dict, replace(file, ".md" => "") => SortedDict{String, SortedDict}(nav_order))

            # If not .md check for (black-listed) folders and add them to the dict. This dict is then
            # filled recursivly.
            elseif isdir(file) && !(folder_blacklisted(file))

                inner_dict = SortedDict{String, SortedDict}(nav_order)
                push!(dict, file => inner_dict)
                append_files!(inner_dict, file)

            end
        end
    catch e
        @info e
    end
end

# Recursivly generates the HTML code for the side nav, given a dict and a path.
function write_html_side_nav(dict, path, level)
    level_indent = "\t"^(3+level)

    html_string = ""
    for key in keys(dict)

        title = apply_formatting(key)
        href_link = string(path, key)
        #@info "SIDE_NAV_GEN: \n" title

        # Checks if the dict, the key element is pointing to is empty, implying it being a page,
        # not a folder
        if length(keys(get(dict, key, SortedDict{String, SortedDict}(nav_order)))) == 0
            level_indent = "\t"^(3+level)
            list_element_string = "$(level_indent)
$title
\n"
            html_string = string(html_string, list_element_string)

        else

            list_element_string = "$(level_indent)
$title\n $(level_indent)\t
 \n"

            inner_dynamic_string = write_html_side_nav(get(dict, key, SortedDict(nav_order)), string(path, key, "/"), level + 2)
            html_string = string(html_string, list_element_string, inner_dynamic_string, "$(level_indent)\t
\n $(level_indent)
\n")

        end
    end

    return html_string
end

function write_html_top_nav(dict, path, level)

    level_indent = "\t"^(5+level)

    html_string = ""
    for key in keys(dict)

        title = apply_formatting(key)
        href_link = string(path, key)

        # Checks if the dict, the key element is pointing to is empty, implying it being a page,
        # not a folder
        if length(keys(get(dict, key, SortedDict{String, SortedDict}(nav_order)))) == 0
            list_element_string = "$(level_indent)
$title
\n"
            html_string = string(html_string, list_element_string)

        else

            list_element_string = "$(level_indent)
$title\n$(level_indent)\t
\n"

            inner_dynamic_string = write_html_top_nav(get(dict, key, SortedDict(nav_order)), string(path, key, "/"), level + 2)
            html_string = string(html_string, list_element_string, inner_dynamic_string, "$(level_indent)\t
\n$(level_indent)
\n")

        end
    end

    return html_string
end

function generate_side_nav()
    # Static first part of HTML file
    file_start =
    "

    \t

    \t\t<form id=\"lunrSearchForm\" name=\"lunrSearchForm\">
    \t\t\t<input class=\"search-input\" name=\"q\" placeholder=\"Search...\" type=\"text\">
    \t\t\t<input type=\"submit\" value=\"Search\" formaction=\"/search/index.html\">
    \t\t</form>
    \t

    \t

    \t\tNavigation
    \t
\n"

    # Static end of file
    file_end =
    "\t

\n"

    # Index files
    @info "Indexing files..."
    file_list = SortedDict{String, SortedDict}(nav_order)
    append_files!(file_list, "")

    # Custom add Home page
    file_string_gen = ""
    file_string_gen = string(file_string_gen, "\t\t
\n")
    file_string_gen = string(file_string_gen, "\t\t\t
Home
\n")

    # Generate dynamic nav HTML
    @info "Generating side-nav..."
    dynamic_string = write_html_side_nav(file_list, "", 0)
    file_string_gen = string(file_string_gen, dynamic_string)

    # Custom add Impressum page
    file_string_gen = string(file_string_gen, "\t\t
\n")

    page = string(file_start, file_string_gen, file_end)

    # write HTML
    open("./_layout/side_nav.html", "w") do io
        write(io, page)
    end
    @info "Successfully written to HTML file!"
end

function generate_top_nav()
    # Static first part of HTML file
    file_start =
    "

    \t

    \t\t
    \t\t\t
    \t\t\t

    \t\t\t\t
Open Positions

    \t\t\t\t
Members

    \t\t\t\t
Teaching Resources

    \t\t\t

    \t\t\t
    \t\t\t\t

    \t\t\t\t\t

    \t\t\t\t\t

    \t\t\t\t\t

    \t\t\t\t

    \t\t\t
    \t\t
    \t\t

    \t\t\t
\n"

     # Static end of file
    file_end =
"\n\t\t\t\t

\t\t\t

\t\t

\t

"

    # Custom add Home page
    file_string_gen = ""
    file_string_gen = string(file_string_gen, "\t\t\t\t\t
Home
\n")

    # Index files
    @info "Indexing files..."
    file_list = SortedDict()
    try
        file_list = SortedDict{String, SortedDict}(nav_order)
        append_files!(file_list, "")
    catch e
        @info e
    end

    # Generate dynamic nav HTML
    @info "Generating top-nav..."
    dynamic_string = write_html_top_nav(file_list, "", 0)
    file_string_gen = string(file_string_gen, dynamic_string)

    # Custom add Impressum and static end of page
    file_string_gen = string(file_string_gen, "\t\t\t\t\t
Impressum
")

    page = string(file_start, file_string_gen, file_end)

    # write HTML
    open("./_layout/nav.html", "w") do io
        write(io, page)
    end
    @info "Successfully written to HTML file!"

end

generate_side_nav()
generate_top_nav()


text:
module Plot_cor

using  PyCall
using  PyPlot
import Modia3D

pyplot_rc = PyCall.PyDict(PyPlot.matplotlib."rcParams")
pyplot_rc["font.size"] = 10.0

vsmall = 0.01
cor0 = 1.0
cor1 = 0.7
cor2 = 0.3
cor3 = 0.1
cor4 = 0.01

vrela = collect(range(0.5*vsmall, 4*vsmall, length=100))
vrelb = collect(range(0         ,   vsmall, length=500))
vrelc = collect(range(0         , 4*vsmall, length=100))
vreld = collect(range(0         , 2*vsmall, length=100))

d_res0a = zeros( length(vrela) )
d_res1a = zeros( length(vrela) )
d_res2a = zeros( length(vrela) )
d_res3a = zeros( length(vrela) )
d_res4a = zeros( length(vrela) )

d_res0b = zeros( length(vrelb) )
d_res1b = zeros( length(vrelb) )
d_res2b = zeros( length(vrelb) )
d_res3b = zeros( length(vrelb) )
d_res4b = zeros( length(vrelb) )

cor_res0 = zeros( length(vrelc) )
cor_res1 = zeros( length(vrelc) )
cor_res2 = zeros( length(vrelc) )
cor_res3 = zeros( length(vrelc) )
cor_res4 = zeros( length(vrelc) )

reg = zeros( length(vreld) )
w   = 0.0
for i in 1:length(vrela)
    d_res0a[i] = Modia3D.resultantDampingCoefficient(cor0,vrela[i],vsmall)
    d_res1a[i] = Modia3D.resultantDampingCoefficient(cor1,vrela[i],vsmall)
    d_res2a[i] = Modia3D.resultantDampingCoefficient(cor2,vrela[i],vsmall)
    d_res3a[i] = Modia3D.resultantDampingCoefficient(cor3,vrela[i],vsmall)
    d_res4a[i] = Modia3D.resultantDampingCoefficient(cor4,vrela[i],vsmall)
end

for i in 1:length(vrelb)
    d_res0b[i] = Modia3D.resultantDampingCoefficient(cor0,vrelb[i],vsmall)
    d_res1b[i] = Modia3D.resultantDampingCoefficient(cor1,vrelb[i],vsmall)
    d_res2b[i] = Modia3D.resultantDampingCoefficient(cor2,vrelb[i],vsmall)
    d_res3b[i] = Modia3D.resultantDampingCoefficient(cor3,vrelb[i],vsmall)
    d_res4b[i] = Modia3D.resultantDampingCoefficient(cor4,vrelb[i],vsmall)
end

for i in 1:length(vrelc)
    cor_res0[i] = Modia3D.resultantCoefficientOfRestitution(cor0,vrelc[i],vsmall)
    cor_res1[i] = Modia3D.resultantCoefficientOfRestitution(cor1,vrelc[i],vsmall)
    cor_res2[i] = Modia3D.resultantCoefficientOfRestitution(cor2,vrelc[i],vsmall)
    cor_res3[i] = Modia3D.resultantCoefficientOfRestitution(cor3,vrelc[i],vsmall)
    cor_res4[i] = Modia3D.resultantCoefficientOfRestitution(cor4,vrelc[i],vsmall)
end

for i in 1:length(vreld)
    reg[i] = Modia3D.regularize(vreld[i],vsmall)
end

figure(1)
clf()
plot(vrelc, cor_res0, vrelc, cor_res1, vrelc, cor_res2, vrelc, cor_res3, vrelc, cor_res4)
grid(true)
xlabel("\$\\dot{\\delta}^- \\; [m/s]\$")
ylabel("\$cor_{reg}\$")
legend(["\$cor = 1.0, v_{small}=0.01 \\; m/s\$",
        "\$cor = 0.7, v_{small}=0.01 \\; m/s\$",
        "\$cor = 0.3, v_{small}=0.01 \\; m/s\$",
        "\$cor = 0.1, v_{small}=0.01 \\; m/s\$",
        "\$cor = 0.0, v_{small}=0.01 \\; m/s\$"],loc="upper right")

figure(2)
clf()
plot(vrela, d_res0a, vrela, d_res1a, vrela, d_res2a, vrela, d_res3a, vrela, d_res4a)
grid(true)
xlabel("\$\\dot{\\delta}^- \\; [m/s]\$")
ylabel("\$d \\; [Ns/m]\$")
legend(["\$cor = 1.0, v_{small}=0.01 \\; m/s\$",
        "\$cor = 0.7, v_{small}=0.01 \\; m/s\$",
        "\$cor = 0.3, v_{small}=0.01 \\; m/s\$",
        "\$cor = 0.1, v_{small}=0.01 \\; m/s\$",
        "\$cor = 0.0, v_{small}=0.01 \\; m/s\$"],loc="upper right")

figure(3)
clf()
plot(vrelb, d_res0b, vrelb, d_res1b, vrelb, d_res2b, vrelb, d_res3b, vrelb, d_res4b)
grid(true)
xlabel("\$\\dot{\\delta}^- \\; [m/s]\$")
ylabel("\$d \\; [Ns/m]\$")
legend(["\$cor = 1.0, v_{small}=0.01 \\; m/s\$",
        "\$cor = 0.7, v_{small}=0.01 \\; m/s\$",
        "\$cor = 0.3, v_{small}=0.01 \\; m/s\$",
        "\$cor = 0.1, v_{small}=0.01 \\; m/s\$",
        "\$cor = 0.0, v_{small}=0.01 \\; m/s\$"],loc="upper right")

figure(4)
clf()
plot(vreld, reg)
grid(true)
xlabel("\$v_{abs}\$")
legend(["\$reg(v_{abs},0.1)\$"])

println("... test/collision/Plot_cor.jl completed.")

end


text:
"""experimental scripts for hyper-parameter optimization"""

using MRMP
import Random: seed!
using Hyperopt
import YAML
import Printf: @printf, @sprintf
import Base.Threads
import Dates
import Logging

function main(config::Dict; pre_compile::Bool = false)

    # create result directory
    date_str = string(Dates.now())
    root_dir = joinpath(
        get(config, "root", joinpath(@__DIR__, "..", "..", "data", "exp")),
        date_str,
    )
    if !pre_compile
        @info @sprintf("result will be saved in %s", root_dir)
        !isdir(root_dir) && mkpath(root_dir)
        additional_info = Dict(
            "git_hash" => read(`git log -1 --pretty=format:"%H"`, String),
            "date" => date_str,
        )
        YAML.write_file(joinpath(root_dir, "config.yaml"), merge(config, additional_info))
    end

    # generate instance
    seed_offset = get(config, "seed_offset", 0)
    num_instances = get(config, "num_instances", 3)
    @info @sprintf("generating %d instances", num_instances)
    instances = begin
        params = Dict([(Symbol(key), val) for (key, val) in config["instance"]])
        delete!(params, Symbol("_target_"))
        target = Meta.parse(config["instance"]["_target_"])
        map(e -> begin
            seed!(e + seed_offset)
            eval(target)(; params...)
        end, 1:num_instances)
    end

    num_search_times = get(config, "num_search_times", 100)
    time_limit = get(config, "time_limit", 30)
    results = Dict()
    for solver_info in config["solvers"]
        solver_name = solver_info["_target_"]
        @info @sprintf(
            "hyper parameter search for %s with %d samples with %d threads, timeout: %f sec",
            solver_name,
            num_search_times,
            Threads.nthreads(),
            time_limit
        )
        params_key = collect(filter(key -> key != "_target_", keys(solver_info)))
        params_cands_str = Dict(map(key -> (Symbol(key), solver_info[key]), params_key))
        ho = Hyperoptimizer(num_search_times; params_cands_str...)
        for (i, params...) in ho
            score = Threads.Atomic{Float64}(0)
            Threads.@threads for k = 1:num_instances
                seed!(k)
                config_init, config_goal, obstacles, ins_params... = instances[k]
                connect = gen_connect(config_init[1], obstacles, ins_params...)
                collide = gen_collide(config_init[1], ins_params...)
                check_goal = MRMP.gen_check_goal(config_goal)
                t = @elapsed begin
                    solution, roadmaps = eval(Meta.parse(solver_name))(
                        config_init,
                        config_goal,
                        connect,
                        collide,
                        check_goal;
                        TIME_LIMIT = time_limit,
                        Dict(pairs(params))...,
                    )
                end
                isnothing(solution) && Threads.atomic_add!(score, 1.0 + t * 0.0001)
            end
            @info @sprintf(
                "fin: %04d/%04d, failure: %04d/%04d, params: %s\n",
                i,
                num_search_times,
                score[],
                num_instances,
                params
            )
            push!(ho.results, score[])
        end
        results[solver_name] = ho

        if !pre_compile
            @info @sprintf("\nsolver:%s\n%s\n", solver_name, ho)
            YAML.write_file(
                joinpath(root_dir, @sprintf("best_params_%s.yaml", solver_name)),
                Dict(solver_name => Dict(zip(ho.params, ho.minimizer))),
            )
        end
    end

    if !pre_compile
        @info "summary:"
        foreach(e -> @info(@sprintf("\nsolver:%s\n%s\n", e...)), results)
        YAML.write_file(
            joinpath(root_dir, "best_params.yaml"),
            Dict([
                (Symbol(solver_name), Dict(zip(ho.params, ho.minimizer))) for
                (solver_name, ho) in results
            ]),
        )
    end
end

function main(args...)
    if length(args) < 2
        @warn @sprintf("two arguments [param_file, eval_file] are required")
        return
    end
    param_file, eval_file = args

    # check file existence
    for file in [param_file, eval_file]
        if !isfile(file)
            @warn @sprintf("%s does not exists", file)
            return
        end
    end

    # create config
    config = merge(
        YAML.load_file(param_file),
        Dict("instance" => YAML.load_file(eval_file)["instance"]),
    )

    # run once to force compilation
    @info "pre-compilation"
    Logging.with_logger(Logging.SimpleLogger(stdout, Logging.Error)) do
        main(
            merge(config, Dict("num_instances" => 1, "num_search_times" => 1));
            pre_compile = true,
        )
    end

    # start experiment
    @info "done, start hypra search"
    main(config)
end

main() = main(ARGS)


text:
using Test
using Suppressor

#assuming being run from test/ folder
include("../src/LambdaTAME.jl")

#= Potential Kwargs:
   tol - (Float)
   iterations - (Int)
   profile- (bool)
   alphas - (Array{Float64,1})
   betas - (Array{Float64,1})
   rank - Int
=#

A = load_ThirdOrderSymTensor("test_tensors/test_tensorA.ssten")
B  = load_ThirdOrderSymTensor("test_tensors/test_tensorA.ssten")

#= TODO:

    align_tensors(graph_A_file::String,graph_B_file::String;
                           ThirdOrderSparse=true,kwargs...)

    get_TAME_ranks(graph_A_file::String,graph_B_file::String)

    distributed_pairwise_alignment(dir::String;kwargs...)

    distributed_pairwise_alignment(files::Array{String,1},dirpath::String;kwargs...)

    distributed_TAME_rank_experiment(files::Array{String,1},dirpath::String;kwargs...)

    self_alignment(dir::String;kwargs...)
=#

@testset "Type Stability" begin

    @suppress_out begin
        @inferred ΛTAME_param_search_profiled(A,B)
        @inferred LowRankTAME_param_search_profiled(A,B)
        @inferred TAME_param_search_profiled(A,B)

        @inferred ΛTAME_param_search(A,B)
        @inferred LowRankTAME_param_search(A,B)
        @inferred TAME_param_search(A,B)
    end
end


text:
module NewtonsMethod

greet() = print("Hello! This function computes the Newton's Method; starting value is 0.3.")

using LinearAlgebra, ForwardDiff

# Function using f'
function newtonroot(f, g; iv = 0.3, tol = 1E-13, maxiter = 1000)

    x_old = iv
    normdiff = Inf
    iter = 1

    while normdiff > tol && iter <= maxiter

        # handle negative first derivative
        if !(g(x_old) == 0)
            x_new = x_old - f(x_old) / g(x_old)
            normdiff = norm(x_new - x_old)
            x_old = x_new
            iter = iter + 1
        else
            return (value = nothing, normdiff = nothing, iter = nothing)
        end

    end

    # handle non convergence
    if normdiff <= tol && iter <= maxiter
        return (value = x_old, normdiff = normdiff, iter = iter)
    else
        return (value = nothing, normdiff = nothing, iter = nothing)
    end

end

# Function using ForwardDiff
function newtonroot(f; iv = 0.3, tol = 1E-13, maxiter = 1000)

    x_old = iv
    normdiff = Inf
    iter = 1

    D(f) = x -> ForwardDiff.derivative(f,x)
    g = D(f)

    sol = newtonroot(f, g; iv = iv, tol = tol, maxiter = maxiter)

    return sol
end

export newtonroot

end

max_stars_repo_name: loforteg/NewtonsMethod

max_issues_repo_name: loforteg/NewtonsMethod