text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -16.424759279035918 <= q <= 20.23452089871623)

                     add_NL_constraint(m, :(1/(1 + exp(-(0.8420310557117703 + 0.34198141401583015*1/(1 + exp(-(-0.5389217244927917 + 0.5798025520034438*1/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + 0.4428632765655749*1/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + 0.9033224385695577*1/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) + -0.1887519896617813*1/(1 + exp(-(-0.23921485861490321 + 0.5030319708107833*1/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.2081481501163327*1/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.2939216839670169*1/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) + -0.2882490310085508*1/(1 + exp(-(-0.6545847926105979 + 0.9774628615021363*1/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.21180887804098658*1/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.7013494059895882*1/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))))))) + 1/(1 + exp(-(-0.35188345522271236 + -0.11415509244934263*1/(1 + exp(-(-0.5389217244927917 + 0.5798025520034438*1/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + 0.4428632765655749*1/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + 0.9033224385695577*1/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) + -0.8088728193024375*1/(1 + exp(-(-0.23921485861490321 + 0.5030319708107833*1/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.2081481501163327*1/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.2939216839670169*1/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) + 0.8271426964720727*1/(1 + exp(-(-0.6545847926105979 + 0.9774628615021363*1/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.21180887804098658*1/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.7013494059895882*1/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))))))) + 1/(1 + exp(-(0.4798469380896484 + 0.8134725364327324*1/(1 + exp(-(-0.5389217244927917 + 0.5798025520034438*1/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + 0.4428632765655749*1/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + 0.9033224385695577*1/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) + -0.12585836968848696*1/(1 + exp(-(-0.23921485861490321 + 0.5030319708107833*1/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.2081481501163327*1/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.2939216839670169*1/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) + 0.613011911779838*1/(1 + exp(-(-0.6545847926105979 + 0.9774628615021363*1/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.21180887804098658*1/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.7013494059895882*1/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
command_pool = CommandPool(device, 0)
cbuffer = first(unwrap(allocate_command_buffers(device, CommandBufferAllocateInfo(command_pool, COMMAND_BUFFER_LEVEL_PRIMARY, 1))))
@test cbuffer isa CommandBuffer

buffer = Buffer(device, 100, BUFFER_USAGE_TRANSFER_DST_BIT, SHARING_MODE_EXCLUSIVE, [])
@test buffer isa Buffer
@test get_buffer_memory_requirements_2_khr ≠ get_buffer_memory_requirements_2

ext: jl

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
export
	# clusterpipeline.jl

	denoise,
    	denoise_pipeline,
	alignment_consensus,
	get_centroid,

	consensus_seq,
	refine_ref,
	consensus_viz,
	disagreements,
	diff_in_homopolymer_region,
	get_coarse_centroid,
	trim_to_refs,


text:
include(joinpath("..", "ComplexOPF.jl","src", "PowSysMod_body.jl"))
include("SDP_decomposition_functions.jl")
include("solve_SDP.jl")
include("solve_minlp.jl")
include("B&B_fixingsome1and0.jl")

struct ROPF_infos
    instance_name::String
    matpower_instance_path::String
    output_instance_path::String
    decomposition::String
    output_decomposition_path::String
    generation::String
    generation_files_path::String
end

struct BB_infos
    search_strategy::String
    branch_strategy::String
    seuil_u::Float64
    seuil_l::Float64
end

function construct_dat_file_ROPF(instance_name, matpower_instance_path, output_path)
    typeofinput = MatpowerRTEROPFSimpleInput
    OPFpbs = load_OPFproblems(typeofinput, matpower_instance_path)
    # Bulding optimization problem
    pb_global = build_globalpb!(OPFpbs)
    pb_global_real = pb_cplx2real(pb_global)
    export_to_dat(pb_global_real, output_path, filename="$(instance_name).dat")
    return
end

function generate_clique_decomposition(instance_name, matpower_instance_path, output_decomposition_path)
    sp = read_sparsity_pattern(matpower_instance_path)
    nb_nodes = size(sp,1)
    nb_edges = (nnz(sp) - size(sp,1))/2
    L, order, nb_added_edges = chordal_ext_cholesky(sp)
    H = construct_graph_from_matrix(L)
    cliques_dict = cliques_max(H,order)
    GC = weighted_graph(cliques_dict)
    A = Prim_algo(GC)
    isdir(joinpath(output_decomposition_path, "cliquetree_cholesky")) || mkpath(joinpath(output_decomposition_path, "cliquetree_cholesky"))
    f = open(joinpath(output_decomposition_path, "cliquetree_cholesky", "$(instance_name)_sdp_cliquetree.txt"), "w")
        for edge in LightGraphs.edges(A)
            clique1 = src(edge)
            clique2 = dst(edge)
            @printf(f, "%10s %20s\n", "B$clique1", "B$clique2")
        end
    close(f)
    isdir(joinpath(output_decomposition_path, "blocks_cholesky")) || mkpath(joinpath(output_decomposition_path, "blocks_cholesky"))
    f = open(joinpath(output_decomposition_path, "blocks_cholesky", "$(instance_name)_sdp_blocks.txt"), "w")
    for (clique, nodes_list) in cliques_dict
        for node in nodes_list
            @printf(f, "%10s %20s\n", "B$clique", "VOLT_$(node)_Re")
            @printf(f, "%10s %20s\n", "B$clique", "VOLT_$(node)_Im")
        end
    end
    close(f)
end

function generate_clique_decomposition_MAT(instance_name, matpower_instance_path, output_decomposition_path)
    sp = read_sparsity_pattern_MAT(matpower_instance_path)
    nb_nodes = size(sp,1)
    nb_edges = (nnz(sp) - size(sp,1))/2
    L, order, nb_added_edges = chordal_ext_cholesky(sp)
    H = construct_graph_from_matrix(L)
    cliques_dict = cliques_max(H,order)
    GC = weighted_graph(cliques_dict)
    A = Prim_algo(GC)
    isdir(joinpath(output_decomposition_path, "cliquetree_cholesky")) || mkpath(joinpath(output_decomposition_path, "cliquetree_cholesky"))
    f = open(joinpath(output_decomposition_path, "cliquetree_cholesky", "$(instance_name)_sdp_cliquetree.txt"), "w")
        for edge in LightGraphs.edges(A)
            clique1 = src(edge)
            clique2 = dst(edge)
            @printf(f, "%10s %20s\n", "B$clique1", "B$clique2")
        end
    close(f)
    isdir(joinpath(output_decomposition_path, "blocks_cholesky")) || mkpath(joinpath(output_decomposition_path, "blocks_cholesky"))
    f = open(joinpath(output_decomposition_path, "blocks_cholesky", "$(instance_name)_sdp_blocks.txt"), "w")
    for (clique, nodes_list) in cliques_dict
        for node in nodes_list
            @printf(f, "%10s %20s\n", "B$clique", "VOLT_$(node)_Re")
            @printf(f, "%10s %20s\n", "B$clique", "VOLT_$(node)_Im")
        end
    end
    close(f)
end

function solve1(ROPF)
    LB_plus, stat_plus = solve_SDP(ROPF, "plus")
    LB_minus, stat_minus = solve_SDP(ROPF, "minus")
    UB_minus = UB_plus = Inf
    isdir("solutions") || mkpath("solutions")
    if !(stat_plus == MOI.FEASIBLE_POINT || stat_plus == MOI.NEARLY_FEASIBLE_POINT) && ! (stat_minus == MOI.FEASIBLE_POINT || stat_minus == MOI.NEARLY_FEASIBLE_POINT)
        println("$(ROPF.instance_name) $(ROPF.generation) : SDP relaxation probably infeasible ")
    end
    if (stat_plus == MOI.FEASIBLE_POINT || stat_plus == MOI.NEARLY_FEASIBLE_POINT)
        #solve MINLP with Knitro
        UB_plus = solve_minlp(ROPF, "plus", [], Dict{String,Float64}())
        mv("knitro_solution.csv", "solutions\\solve1_solution_$(ROPF.instance_name)_plus.csv", force=true)
    end
    if (stat_minus == MOI.FEASIBLE_POINT || stat_minus == MOI.NEARLY_FEASIBLE_POINT)
        #solve MINLP with Knitro
        UB_minus = solve_minlp(ROPF, "minus", [], Dict{String,Float64}())
        mv("knitro_solution.csv", "solutions\\solve1_solution_$(ROPF.instance_name)_minus.csv", force=true)
    end
    return UB_plus, LB_plus, UB_minus, LB_minus
end

function solve1_refgen(ROPF)
    LB, stat = solve_SDP_refgen(ROPF)
    UB = Inf
    isdir("solutions") || mkpath("solutions")
    if !(stat == MOI.FEASIBLE_POINT || stat == MOI.NEARLY_FEASIBLE_POINT)
        println("$(ROPF.instance_name) $(ROPF.generation) : SDP relaxation probably infeasible ")
    end
    if (stat == MOI.FEASIBLE_POINT || stat == MOI.NEARLY_FEASIBLE_POINT)
        #solve MINLP with Knitro
        UB = solve_minlp_refgen(ROPF, [], Dict{String,Float64}())
        mv("knitro_solution.csv", "solutions\\solve1_refgen_solution_$(ROPF.instance_name).csv", force=true)
    end
    return UB, LB
end

function solve1_refgen_allshunts1(ROPF)
    output_instance_path = ROPF.output_instance_path
    output_decomposition_path = ROPF.output_decomposition_path
    generation_files_path = ROPF.generation_files_path
    FORMULATION = ROPF.decomposition
    INSTANCE_NAME = ROPF.instance_name
    generation = ROPF.generation
    originalSTDOUT = stdout
    outpath = joinpath("Mosek_runs")
    isdir(outpath) || mkpath(outpath)
    outlog = open(joinpath(outpath,"$(INSTANCE_NAME)_$(generation)_$(FORMULATION).log"), "w")
    redirect_stdout(outlog)
    instance_dat_file_path = joinpath(output_instance_path, "$(INSTANCE_NAME).dat")
    Pinput_csv_file = joinpath(generation_files_path,"$(INSTANCE_NAME)_$(generation).csv")
    solution_file = ""
    refbus = readlines(joinpath(generation_files_path, "$(INSTANCE_NAME)_refbus.csv"))[1]
    λ, Sgen_var_list, SDP_var_list, Bin_var_list, dict_quad_ctr, dict_bounds_ctr, dict_constants_ctr, dict_Bin_ctr,
     dict_MONO, dict_linear_ctr = read_dat_file(instance_dat_file_path)
    cliques_dict, CLIQUE_TREE = read_blocks(output_decomposition_path, FORMULATION, INSTANCE_NAME)
     value_bins, value_SDP_var, LB, stat = construct_SDP_refgen(cliques_dict, CLIQUE_TREE, Pinput_csv_file, refbus, Sgen_var_list, SDP_var_list, Bin_var_list,
          dict_quad_ctr, dict_linear_ctr, dict_bounds_ctr, dict_constants_ctr, dict_Bin_ctr, dict_MONO,  Dict(bin => 1 for bin in Bin_var_list), solution_file)
      index_var = Dict{String, Int64}()
      nb_bin = length(Bin_var_list)
      for i=1:nb_bin
          var = Bin_var_list[i]
          index_var[var] = i
      end
    UB = Inf
    isdir("solutions") || mkpath("solutions")
    if !(stat == MOI.FEASIBLE_POINT || stat == MOI.NEARLY_FEASIBLE_POINT)
        println("$(ROPF.instance_name) $(ROPF.generation) : SDP relaxation probably infeasible ")
    end
    if (stat == MOI.FEASIBLE_POINT || stat == MOI.NEARLY_FEASIBLE_POINT)
        #solve MINLP with Knitro
        UB = solve_minlp_refgen(ROPF, ones(nb_bin), index_var)
        mv("knitro_solution.csv", "solutions\\solve1_allshunts1_refgen_solution_$(ROPF.instance_name).csv", force=true)
    end
    return UB, LB
end

function solve2(ROPF, max_time)
    LB_plus, stat_plus = solve_SDP(ROPF, "plus")
    LB_minus, stat_minus = solve_SDP(ROPF, "minus")
    UB_minus = UB_plus = Inf
    if !(stat_plus == MOI.FEASIBLE_POINT || stat_plus == MOI.NEARLY_FEASIBLE_POINT) && ! (stat_minus == MOI.FEASIBLE_POINT || stat_minus == MOI.NEARLY_FEASIBLE_POINT)
        println("$(ROPF.instance_name) $(ROPF.generation) : SDP relaxation probably infeasible ")
    end
    if (stat_plus == MOI.FEASIBLE_POINT || stat_plus == MOI.NEARLY_FEASIBLE_POINT)
        #B&B algo
        BB_parameters = BB_infos("deepfirst", "1", 0.9, 0.0001)
        (UB_plus, nb_nodes, open_nodes) = BandB_fixingsome1and0(ROPF, "plus", BB_parameters, max_time)
    end
    if (stat_minus == MOI.FEASIBLE_POINT || stat_minus == MOI.NEARLY_FEASIBLE_POINT)
        #B&B algo
        BB_parameters = BB_infos("deepfirst", "1", 0.9, 0.0001)
        (UB_minus, nb_nodes, open_nodes) = BandB_fixingsome1and0(ROPF, "minus", BB_parameters, max_time)
    end
    return UB_plus, LB_plus, UB_minus, LB_minus
end

function solve2_refgen(ROPF, max_time)
    LB, stat = solve_SDP_refgen(ROPF)
    UB = Inf
    isdir("solutions") || mkpath("solutions")
    if !(stat == MOI.FEASIBLE_POINT || stat == MOI.NEARLY_FEASIBLE_POINT)
        println("$(ROPF.instance_name) $(ROPF.generation) : SDP relaxation probably infeasible ")
    end
    if (stat == MOI.FEASIBLE_POINT || stat == MOI.NEARLY_FEASIBLE_POINT)
        #B&B algo
        BB_parameters = BB_infos("deepfirst", "1", 0.9, 0.0001)
        (UB, nb_nodes, open_nodes) = BandB_fixingsome1and0_refgen(ROPF, BB_parameters, max_time)
    end
    return UB, LB
end


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

JLLWrappers.@generate_main_file_header("LibTorchCAPI")
JLLWrappers.@generate_main_file("LibTorchCAPI", UUID("35125b28-eb97-566c-9d83-c787ae32e58f"))
end  # module LibTorchCAPI_jll


text:
#Generate Julia composite types for ROS messages

using PyCall

export @rosimport, rostypegen, rostypereset

#Type definitions
#Composite types for internal use. Keeps track of the imported types and helps
#keep code generation orderly.
abstract type ROSModule end
mutable struct ROSPackage
    name::String
    msg::ROSModule
    srv::ROSModule
    function ROSPackage(pkgname::String)
        pkg = new(pkgname)
        pkg.msg = ROSMsgModule(pkg)
        pkg.srv = ROSSrvModule(pkg)
        pkg
    end
end
struct ROSMsgModule <: ROSModule
    pkg::ROSPackage
    members::Vector{String}
    deps::Set{String}
    ROSMsgModule(pkg) = new(pkg, String[], Set{String}())
end
struct ROSSrvModule <: ROSModule
    pkg::ROSPackage
    members::Vector{String}
    deps::Set{String}
    ROSSrvModule(pkg) = new(pkg, String[], Set{String}())
end

#These global objects maintain the hierarchy from multiple calls to
#`@rosimport` and keep the link to the Python objects whenever communication
#goes between RobotOS and rospy.
const _rospy_imports = Dict{String,ROSPackage}()
const _rospy_objects = Dict{String,PyObject}()
const _rospy_modules = Dict{String,PyObject}()

const _ros_builtin_types = Dict{String,DataType}(
    "bool"    => Bool,
    "int8"    => Int8,
    "int16"   => Int16,
    "int32"   => Int32,
    "int64"   => Int64,
    "uint8"   => UInt8,
    "uint16"  => UInt16,
    "uint32"  => UInt32,
    "uint64"  => UInt64,
    "float32" => Float32,
    "float64" => Float64,
    "string"  => String,
    "time"    => Time,
    "duration"=> Duration,
    #Deprecated by ROS but supported here
    "char"    => UInt8,
    "byte"    => Int8,
    )

#Abstract supertypes of all generated types
abstract type AbstractMsg end
abstract type AbstractSrv end
abstract type AbstractService end

_is_tuple_expr(input) = input isa Expr && input.head == :tuple
_is_colon_expr(input) = input isa Expr && input.head == :call && input.args[1] == :(:)
_is_dot_expr(input) = input isa Expr && input.head == :(.)

"""
    @rosimport

Import ROS message or service types into Julia. Call `rostypegen()` after all `@rosimport` calls.
Package or type dependencies are also imported automatically as needed.

Example usages:
```julia
  @rosimport geometry_msgs.msg.PoseStamped
  @rosimport sensor_msgs.msg: Image, Imu
  @rosimport nav_msgs.srv.GetPlan
```
"""
macro rosimport(input)
    #Rearranges the expression into a RobotOS._rosimport call. Input comes in as a single package
    #qualified expression, or as a tuple expression where the first element is the same as the
    #single expression case. Most of the code is just error checking that the input takes that form.
    @assert _is_tuple_expr(input) || _is_colon_expr(input) || _is_dot_expr(input) "Improper @rosimport input"
    if _is_tuple_expr(input)
        @assert _is_colon_expr(input.args[1]) "Improper @rosimport input, first argument needs ':' following"
        pkg, ismsg, typ = _pkgtype_import(input.args[1])
        types = String[typ]
        for t in input.args[2:end]
            @assert isa(t, Symbol) "Type name ($(string(t))) not a symbol"
            push!(types, string(t))
        end
        return :(_rosimport($pkg, $ismsg, $types...))
    else
        pkg, ismsg, typ = _pkgtype_import(input)
        return :(_rosimport($pkg, $ismsg, $typ))
    end
end
#Return the pkg and types strings for a single expression of form:
#  pkg.[msg|srv].type or pkg.[msg|srv]:type
function _pkgtype_import(input::Expr)
    @assert _is_colon_expr(input) || _is_dot_expr(input) "Improper @rosimport input"
    p_ms, t = _is_colon_expr(input) ? (input.args[2], input.args[3]) : (input.args[1], input.args[2])
    @assert _is_dot_expr(p_ms) "Improper @rosimport input"
    p = p_ms.args[1]
    @assert isa(p, Symbol) "Package name ($(string(p))) not a symbol"
    @assert isa(p_ms.args[2], QuoteNode) "Improper @rosimport input"
    m_or_s = p_ms.args[2].value
    @assert m_or_s in (:msg,:srv) "Improper @rosimport input"
    ps = string(p)
    msb = m_or_s == :msg
    ts = ""
    if isa(t, Symbol)
        ts = string(t)
    elseif isa(t, Expr)
        @assert length(t.args) == 1 "Type name ($(t)) not a symbol"
        tsym = t.args[1]
        @assert isa(tsym, Symbol) "Type name ($(string(tsym))) not a symbol"
        ts = string(tsym)
    elseif isa(t, QuoteNode)
        tsym = t.value
        @assert isa(tsym, Symbol) "Type name ($(string(tsym))) not a symbol"
        ts = string(tsym)
    end
    return ps,msb,ts
end
#Import a set of types from a single package
function _rosimport(package::String, ismsg::Bool, names::String...)
    global _rospy_imports
    if ! haskey(_rospy_imports, package)
        @debug("Importing new package: ",package,".", ismsg ? "msg" : "srv")
        _rospy_imports[package] = ROSPackage(package)
    end
    rospypkg = _rospy_imports[package]
    for n in names
        addtype!(ismsg ? rospypkg.msg : rospypkg.srv, n)
    end
end

"""
    rostypegen(rosrootmod::Module=Main)

Initiate the Julia type generation process after importing some ROS types. Creates modules in
rootrosmod (default is `Main`) with the same behavior as imported ROS modules in python.
Should only be called once, after all `@rosimport` statements are done.
"""
function rostypegen(rosrootmod::Module=Main)
    global _rospy_imports
    pkgdeps = _collectdeps(_rospy_imports)
    pkglist = _order(pkgdeps)
    for pkg in pkglist
        buildpackage(_rospy_imports[pkg], rosrootmod)
    end
end

"""
    rostypereset()

Clear out the previous `@rosimport`s, returning the type generation to its original state. Cannot do
anything about already generated modules in `Main`.
"""
function rostypereset()
    global _rospy_imports
    global _rospy_objects
    empty!(_rospy_imports)
    empty!(_rospy_objects)
    nothing
end

#Populate the module with a new message type. Import and add dependencies first
#so they will appear first in the generated code.
function addtype!(mod::ROSMsgModule, typ::String)
    global _rospy_objects
    if !(typ in mod.members)
        @debug("Message type import: ", _fullname(mod), ".", typ)
        if _nameconflicts(typ)
            @warn("Message type '$typ' conflicts with Julia builtin, " *
                  "will be imported as '$(_jl_safe_name(typ,"Msg"))'")
        end
        pymod, pyobj = _pyvars(_fullname(mod), typ)

        deptypes = pyobj[:_slot_types]
        _importdeps!(mod, deptypes)

        push!(mod.members, typ)
        _rospy_objects[_rostypestr(mod, typ)] = pyobj
    end
end

#Populate the module with a new service type. Import and add dependencies
#first.
function addtype!(mod::ROSSrvModule, typ::String)
    global _rospy_objects
    if !(typ in mod.members)
        @debug("Service type import: ", _fullname(mod), ".", typ)
        pymod, pyobj = _pyvars(_fullname(mod), typ)

        if ! haskey(pyobj, "_request_class")
            error(string("Incorrect service name: ", typ))
        end

        #Immediately import dependencies from the Request/Response classes
        #Repeats are OK
        req_obj = pymod[string(typ,"Request")]
        resp_obj = pymod[string(typ,"Response")]
        deptypes = [req_obj[:_slot_types]; resp_obj[:_slot_types]]
        _importdeps!(mod, deptypes)

        push!(mod.members, typ)
        fulltypestr = _rostypestr(mod, typ)
        _rospy_objects[fulltypestr] = pyobj
        _rospy_objects[string(fulltypestr,"Request")] = req_obj
        _rospy_objects[string(fulltypestr,"Response")] = resp_obj
    end
end

#Return the python module and python object for a particular type
function _pyvars(modname::String, typ::String)
    pymod = _import_rospy_pkg(modname)
    pyobj =
        try pymod[typ]
        catch ex
            isa(ex, KeyError) || rethrow(ex)
            error("Message type '$typ' not found in ROS package '$modname', ",
                  "check the corresponding @rosimport call")
        end
    pymod, pyobj
end

#Continue the import process on a list of dependencies. Called by `addtype!`
#and calls `addtype!` to complete the dependency recursion.
function _importdeps!(mod::ROSModule, deps::Vector)
    global _rospy_imports
    for d in deps
        #We don't care about array types when doing dependency resolution
        dclean = _check_array_type(d)[1]
        if ! haskey(_ros_builtin_types, dclean)
            @debug("Dependency: ", d)
            pkgname, typename = _splittypestr(dclean)

            @debug_addindent
            #Create a new ROSPackage if needed
            if ! haskey(_rospy_imports, pkgname)
                @debug("Creating new package: ", pkgname)
                _rospy_imports[pkgname] = ROSPackage(pkgname)
            end

            #Dependencies will always be messages only
            depmod = _rospy_imports[pkgname].msg
            #pushing to a set does not create duplicates
            push!(mod.deps, _name(depmod))

            addtype!(depmod, typename)
            @debug_subindent
        end
    end
end

#Bring in the python modules as needed
function _import_rospy_pkg(package::String)
    global _rospy_modules
    if ! haskey(_rospy_modules, package)
        @debug("Importing python package: ", package)
        try
            _rospy_modules[package] = pyimport(package)
        catch ex
            show(ex)
            error("python import error: $(ex.val[:args][1])")
        end
    end
    _rospy_modules[package]
end

#The function that creates and fills the generated top-level modules
function buildpackage(pkg::ROSPackage, rosrootmod::Module)
    @debug("Building package: ", _name(pkg))

    #Create the top-level module for the package in Main
    pkgsym = Symbol(_name(pkg))
    pkgcode = :(module ($pkgsym) end)
    pkginitcode = :(function __init__() end)

    #Add msg and srv submodules if needed
    @debug_addindent
    if length(pkg.msg.members) > 0
        msgmod = :(module msg end)
        msgcode = modulecode(pkg.msg, rosrootmod)
        for expr in msgcode
            push!(msgmod.args[3].args, expr)
        end
        push!(pkgcode.args[3].args, msgmod)
        for typ in pkg.msg.members
            push!(pkginitcode.args[2].args, :(@rosimport $(pkgsym).msg: $(Symbol(typ))))
        end
    end
    if length(pkg.srv.members) > 0
        srvmod = :(module srv end)
        srvcode = modulecode(pkg.srv, rosrootmod)
        for expr in srvcode
            push!(srvmod.args[3].args, expr)
        end
        push!(pkgcode.args[3].args, srvmod)
        for typ in pkg.srv.members
            push!(pkginitcode.args[2].args, :(@rosimport $(pkgsym).srv: $(Symbol(typ))))
        end
    end
    push!(pkgcode.args[3].args, :(import RobotOS.@rosimport))
    push!(pkgcode.args[3].args, pkginitcode)
    pkgcode = Expr(:toplevel, pkgcode)
    rosrootmod.eval(pkgcode)
    @debug_subindent
end

#Generate all code for a .msg or .srv module
function modulecode(mod::ROSModule, rosrootmod::Module)
    @debug("submodule: ", _fullname(mod))
    modcode = Expr[]

    #Common imports
    push!(modcode,
        quote
            using PyCall
            import Base: convert, getindex
            import RobotOS
            import RobotOS.Time
            import RobotOS.Duration
            import RobotOS._typedefault
            import RobotOS._typerepr
        end
    )
    #Import statement specific to the module
    append!(modcode, _importexprs(mod, rosrootmod))
    #The exported names
    push!(modcode, _exportexpr(mod))

    #The generated type codes
    @debug_addindent
    for typ in mod.members
        typecode = buildtype(mod, typ)
        append!(modcode, typecode)
    end
    @debug_subindent
    modcode
end

#The imports specific to each module, including dependant packages
function _importexprs(mod::ROSMsgModule, rosrootmod::Module)
    imports = Expr[:(import RobotOS.AbstractMsg)]
    othermods = filter(d -> d != _name(mod), mod.deps)
    append!(imports, [Expr(:using, Expr(:., fullname(rosrootmod)..., Symbol(m), :msg)) for m in othermods])
    imports
end
function _importexprs(mod::ROSSrvModule, rosrootmod::Module)
    imports = Expr[
        :(import RobotOS.AbstractSrv),
        :(import RobotOS.AbstractService),
        :(import RobotOS._srv_reqtype),
        :(import RobotOS._srv_resptype)
    ]
    append!(imports, [Expr(:using, Expr(:., fullname(rosrootmod)..., Symbol(m), :msg)) for m in mod.deps])
    imports
end

#The exported names for each module
function _exportexpr(mod::ROSMsgModule)
    exportexpr = Expr(:export)
    for m in mod.members
        push!(exportexpr.args, Symbol(_jl_safe_name(m,"Msg")))
    end
    exportexpr
end
function _exportexpr(mod::ROSSrvModule)
    exportexpr = Expr(:export)
    for typ in mod.members
        push!(exportexpr.args,
            Symbol(typ),
            Symbol(string(typ,"Request")),
            Symbol(string(typ,"Response"))
        )
    end
    exportexpr
end

#All the generated code for a generated message type
function buildtype(mod::ROSMsgModule, typename::String)
    global _rospy_objects
    fulltypestr = _rostypestr(mod, typename)
    pyobj = _rospy_objects[fulltypestr]
    memnames = pyobj[:__slots__]
    memtypes = pyobj[:_slot_types]
    members = collect(zip(memnames, memtypes))

    typecode(fulltypestr, :AbstractMsg, members)
end

#All the generated code for a generated service type
#Will create 3 different composite types.
function buildtype(mod::ROSSrvModule, typename::String)
    global _rospy_objects

    req_typestr = _rostypestr(mod, string(typename,"Request"))
    reqobj = _rospy_objects[req_typestr]
    memnames = reqobj[:__slots__]
    memtypes = reqobj[:_slot_types]
    reqmems = collect(zip(memnames, memtypes))
    pyreq  = :(RobotOS._rospy_objects[$req_typestr])
    reqexprs  = typecode(req_typestr, :AbstractSrv, reqmems)

    resp_typestr = _rostypestr(mod, string(typename,"Response"))
    respobj = _rospy_objects[resp_typestr]
    memnames = respobj[:__slots__]
    memtypes = respobj[:_slot_types]
    respmems = collect(zip(memnames, memtypes))
    pyresp = :(RobotOS._rospy_objects[$resp_typestr])
    respexprs = typecode(resp_typestr, :AbstractSrv, respmems)

    defsym = Symbol(typename)
    reqsym = Symbol(string(typename,"Request"))
    respsym = Symbol(string(typename,"Response"))
    srvexprs = Expr[
        :(struct $defsym <: AbstractService end),
        :(_typerepr(::Type{$defsym}) = $(_rostypestr(mod,typename))),
        :(_srv_reqtype(::Type{$defsym}) = $reqsym),
        :(_srv_resptype(::Type{$defsym}) = $respsym),
    ]
    [reqexprs; respexprs; srvexprs]
end

#Create the core generated expressions for a native Julia message type that has
#data fields and interchanges with a python counterpart:
# (1) the 'type ... end' block
# (2) Default outer constructer with no arguments
# (3) convert(PyObject, ...)
# (4) convert(..., o::PyObject)
# (5) getindex for accessing member constants
function typecode(rosname::String, super::Symbol, members::Vector)
    tname = _splittypestr(rosname)[2]
    @debug("Type: ", tname)

    #generated code should not conflict with julia built-ins
    #some messages need renaming
    suffix = if super == :AbstractMsg; "Msg"
         elseif super == :AbstractSrv; "Srv"
         else; "ROS" end
    jlsym = Symbol(_jl_safe_name(tname,suffix))

    exprs = Expr[]
    #First the empty expressions
    #(1) Type declaration
    push!(exprs, :(
        mutable struct $jlsym <: $super
            #Generated code here
        end
    ))
    #(2) Default constructor, but only if the type has members
    if length(members) > 0
        push!(exprs, :(
            function $jlsym()
                $jlsym() #Generated code inside parens here
            end
        ))
    else
        push!(exprs, :())
    end
    #(3) Convert to PyObject
    push!(exprs, :(
        function convert(::Type{PyObject}, o::$jlsym)
            py = pycall(RobotOS._rospy_objects[$rosname], PyObject)
            #Generated code here
            py
        end
    ))
    #(4) Convert from PyObject
    push!(exprs, :(
        function convert(jlt::Type{$jlsym}, o::PyObject)
            if convert(String, o["_type"]) != _typerepr(jlt)
                throw(InexactError(:convert, $jlsym, o))
            end
            jl = $jlsym()
            #Generated code here
            jl
        end
    ))
    #(5) Accessing member variables through getindex
    push!(exprs, :(
        getindex(::Type{$jlsym}, s::Symbol) = RobotOS._rospy_objects[$rosname][s]
    ))

    #Now add the meat to the empty expressions above
    for (namestr,typ) in members
        @debug_addindent
        _addtypemember!(exprs, namestr, typ)
        @debug_subindent
    end
    push!(exprs, :(_typerepr(::Type{$jlsym}) = $rosname))
    exprs
end

#Add the generated expression from a single member of a type, either built-in
#or ROS type. `exprs` is the Expr objects of the items created in `typecode`.
#Maybe this can be factored into something nicer.
function _addtypemember!(exprs, namestr, typestr)
    @debug("$namestr :: $typestr")
    typeargs  = exprs[1].args[3].args
    consargs  = exprs[2].args[2].args[2].args
    pyconargs = exprs[3].args[2].args
    jlconargs = exprs[4].args[2].args

    if typestr == "char" || typestr == "byte"
        @warn("Use of type '$typestr' is deprecated in message definitions, " *
              "use '$(lowercase(string(_ros_builtin_types[typestr])))' instead.")
    end

    typestr, arraylen = _check_array_type(typestr)
    if _isrostype(typestr)
        j_typ = Symbol(_splittypestr(typestr)[2])
        #Default has to be deferred until the types exist
        j_def = Expr(:call, j_typ)
    else
        if ! haskey(_ros_builtin_types, typestr)
            error("Message generation; unknown type '$typestr'")
        end
        j_typ = _ros_builtin_types[typestr]
        #Compute the default value now
        j_def = _typedefault(j_typ)
    end

    namesym = Symbol(namestr)
    if arraylen >= 0
        memexpr = :($namesym::Array{$j_typ,1})
        defexpr = :([$j_def for i = 1:$arraylen])
        jlconexpr = :(jl.$namesym = convert(Array{$j_typ,1}, o[$namestr]))

        #uint8[] is string in rospy and PyCall's conversion to bytearray is
        #rejected by ROS
        if j_typ == :UInt8
            pyconexpr = :(py[$namestr] =
                pycall(pybuiltin("str"), PyObject, PyObject(o.$namesym))
            )
        elseif _isrostype(typestr)
            pyconexpr = :(py[$namestr] =
                convert(Array{PyObject,1}, o.$namesym))
        else
            pyconexpr = :(py[$namestr] = o.$namesym)
        end
    else
        memexpr = :($namesym::$j_typ)
        defexpr = j_def
        jlconexpr = :(jl.$namesym = convert($j_typ, o[$namestr]))
        pyconexpr = :(py[$namestr] = convert(PyObject, o.$namesym))
    end
    push!(typeargs, memexpr)
    insert!(jlconargs, length(jlconargs), jlconexpr)
    insert!(pyconargs, length(pyconargs), pyconexpr)
    push!(consargs, defexpr)
end

#Build a String => Iterable{String} object from the individual package
#dependencies.
function _collectdeps(pkgs::Dict{S, ROSPackage}) where S <: AbstractString
    deps = Dict{S, Set{S}}()
    for pname in keys(pkgs)
        if ! haskey(deps, pname)
            deps[pname] = Set{S}()
        end
        union!(deps[pname], pkgs[pname].msg.deps)
        union!(deps[pname], pkgs[pname].srv.deps)
    end
    deps
end

#Produce an order of the keys of d that respect their dependencies.
#Assumed to be Dict(String => Iterable{String})
function _order(d::Dict)
    trecurse!(currlist, d, t) = begin
        if !(t in currlist)
            if haskey(d, t) #do dependencies first
                for dt in d[t]
                    if dt != t
                        trecurse!(currlist, d, dt)
                    end
                end
                #Now it's ok to add it
                push!(currlist, t)
            end
        end
    end
    tlist = String[]
    for t in keys(d)
        trecurse!(tlist, d, t)
    end
    tlist
end

_rostypestr(mod::ROSModule, name::String) = string(_name(mod),"/",name)
function _splittypestr(typestr::String)
    if ! _isrostype(typestr)
        error(string("Invalid message type '$typestr', ",
                     "use 'package_name/type_name'"))
    end
    rospkg, typ = map(ascii, split(typestr, '/'))
    rospkg, typ
end
#Valid ROS type string is all word chars split by a single forward slash, with
#optional square brackets for array types
_isrostype(s::String) = occursin(r"^\w+/\w+(?:\[\d*\])?$", s)

#Sanitize a string by checking for and removing brackets if they are present
#Return the sanitized type and the number inside the brackets if it is a fixed
#size type. Returns 0 if variable size (no number), -1 if no brackets
function _check_array_type(typ::String)
    arraylen = -1
    m = match(r"^([\w/]+)\[(\d*)\]$", typ)
    if m != nothing
        btype = m.captures[1]
        if isempty(m.captures[2])
            arraylen = 0
        else
            arraylen = parse(Int, m.captures[2])
        end
    else
        btype = typ
    end
    ascii(btype), arraylen
end

#Get the rospy PyObject corresponding to a generated type
function _get_rospy_class(typ::DataType)
    global _rospy_objects
    rospycls =
        try
            _rospy_objects[_typerepr(typ)]
        catch ex
            if isa(ex, KeyError)
                error("Type ($typ) is not generated")
            else
                error("Type ($typ) is not a valid message type")
            end
        end
    rospycls
end

#Overwrite PyCall's default constructor to call the `convert` functions generated here
PyCall.PyObject(m::AbstractMsg) = convert(PyCall.PyObject, m)
PyCall.PyObject(s::AbstractSrv) = convert(PyCall.PyObject, s)
PyCall.PyObject(s::AbstractService) = convert(PyCall.PyObject, s)

_jl_safe_name(name::AbstractString, suffix) = _nameconflicts(name) ?
    string(name,suffix) :
    name

#Check if the type name conflicts with a Julia builtin. Currently this is only
#some of the messages from the std_msgs.msg package
_nameconflicts(typename::String) = isdefined(Base, Symbol(typename))

#Get a default value for any builtin ROS type
_typedefault(::Type{T}) where {T <: Real} = zero(T)
_typedefault(::Type{String}) = ""
_typedefault(::Type{Time}) = Time(0,0)
_typedefault(::Type{Duration}) = Duration(0,0)

#Default method to get the "pkg/type" string from a generated DataType.
#Extended by the generated modules.
_typerepr(::Type{T}) where {T} = error("Not a ROS type")

#Default method to get the request/response datatypes for a generated service
_srv_reqtype( ::Type{T}) where {T} = error("Not a ROS Service type")
_srv_resptype(::Type{T}) where {T} = error("Not a ROS Service type")

#Accessors for the package name
_name(p::ROSPackage) = p.name
_name(m::ROSModule) = _name(m.pkg)

#Get the full ROS name for a module (e.g., 'std_msgs.msg' or nav_msgs.srv')
_fullname(m::ROSMsgModule) = string(_name(m), ".msg")
_fullname(m::ROSSrvModule) = string(_name(m), ".srv")


text:
################################################################################
#  Copyright 2022, Jose Antonio de la O Serna (UANL), Tom Van Acker (BASF)     #
################################################################################
# TFT.jl                                                                       #
# A Julia package for Taylor-Fourier Transform.                                #
# See http://github.com/timmyfaraday/TFT.jl                                    #
################################################################################

"""
    TFT.harmonic_estimator(prob::AbstractDTFTProblem, H::Int)

Function to obtain the up-to-Dth-degree derivative of the Hth-harmonic dynamic
phasor.

Input:
- `prob::AbstractDTFTProblem`   | DTFT problem struct
- `H::Int`                      | harmonic number [-]

Output:
- `X::Matrix{<:Complex}`        | up-to-Dth-degree derivative of the
                                | Hth-harmonic dynamic phasor
"""
function harmonic_estimator(prob::AbstractDTFTProblem, H::Int)
    # define the extended normalized time range `rU`
    Δu  = 1 / prob.N
    rU  = 0.0:Δu:(prob.K + 1 - Δu)

    # samples of up-to-Dth-degree derivative of the Kth-degree o-spline `Φ`
    Φ   = sample_ospline(prob.D, prob.K, prob.N)
    # update `Φ` to reflect the samples of up-to-Dth-degree derivative of the
    # Hth-harmonic bandpass filter
    Y   = Φ .* exp.((2 * pi * im * H) .* rU) ./ prob.N .*
            (prob.F).^collect(0:prob.D)'

    # define the appropriate range in the convolution `rC` using an offset `O`:
    # if (K+1)*N is even    → O = ((K+1)*N) / 2 + 1
    # if (K+1)*N is oneven  → O = ((K+1)*N - 1) / 2 + 1
    # NB: the Julia DSP pkg does not allow for the keyword 'same' in its conv-
    # function, as is the case in MATLAB, the range `rC` mimics that behavior
    O   = floor(Int, ((prob.K + 1) * prob.N) / 2) + 1
    rC  = O:O+length(prob.s)-1

    # return the up-to-Dth-degree derivative of the H-th harmonic dynamic phasor
    return _DSP.conv(prob.s, Y)[rC, :]
end

"""
    TFT.fast_estimator(prob::AbstractDTFTProblem)

Function to obtain the Fast Taylor Fourier transform FTFT. Uses the samples of
the up-to-Dth-degree derivatives of Kth-degree o-spline in matrix Φ.

Input:
- `prob::AbstractDTFTProblem`       | DTFT problem struct

Output:
- `X::Dict{Int,Matrix{<:Complex}}`  | dictionary contain up-to-Dth-degree
                                    | derivative of the dynamic phasors for the
                                    | full set of harmonics
"""
function fast_estimator(prob::AbstractDTFTProblem)
    # initialize a matrix of the size (|s|,|h|,D) to store the up-to-Dth-degree
    # derivative of the dynamic phasors for the full set of harmonics
    X   = zeros(Number, length(prob.s), prob.N, prob.D+1)

    # samples of up-to-Dth-degree derivative of the Kth-degree o-spline `Φ`
    Φ   = sample_ospline(prob.D, prob.K, prob.N)

    # determine the shift to center the 'current' signal point
    Δs  = floor(Int, prob.N * (prob.K+1) / 2)

    for (ni,ns) in enumerate(prob.s) if Δs < ni <= length(prob.s) - Δs
        # determine the selected range of the signal with 'center' ni
        rS  = (ni-Δs):(ni+Δs-1)
        # determine the Hadamar product of the osplines and the selected signal
        hd  = Φ .* prob.s[rS] .* (-prob.F).^collect(0:prob.D)'
        # sum the Hadamar product over o-spline degrees
        shd = [sum(hd[nh:prob.N:end, :], dims=1) for nh in 1:prob.N]
        # perform the fft and store the result at ni
        X[ni,:,:] = _FFTW.fft(reduce(vcat,shd), [1]) / prob.N
    end end

    # return a dictionary of spliced arrays over the harmonic numbers
    return Dict(nh => X[:,nh+1,:] for nh in prob.h)
end


text:
function solve_tsp(tsp_file_path::AbstractString, parameters=AlgorithmParameters(); verbose=true)
    tsp = TSPLIB.readTSP(tsp_file_path)
    return solve_tsp(tsp, parameters; verbose=verbose, use_dist_mtx=true)
end

function solve_tsp(tsp::TSP, parameters=AlgorithmParameters(); verbose=true, use_dist_mtx=false)
    n = tsp.dimension
    x = tsp.nodes[:, 1]
    y = tsp.nodes[:, 2]
    serv_time = zeros(size(x))
    dem = ones(size(x))
    vehicleCapacity = n
    maximum_number_of_vehicles = 1

    c_dist_mtx = Matrix(tsp.weights')
    for i in 1:n
        c_dist_mtx[i, i] = 0.0
    end

    if use_dist_mtx
        # need to input dist_mtx' instead of dist_mtx
        # Julia: column-first indexing
        # C: row-first indexing
        return c_api_solve_cvrp_dist_mtx(n, x, y, c_dist_mtx, serv_time, dem, vehicleCapacity, maximum_number_of_vehicles, parameters, verbose)
    else
        return c_api_solve_cvrp(n, x, y, serv_time, dem, vehicleCapacity, maximum_number_of_vehicles, parameters, verbose)
    end
end

function solve_tsp(dist_mtx::Matrix, parameters=AlgorithmParameters(); verbose=true, x_coords=zeros(size(dist_mtx,1)), y_coords=zeros(size(dist_mtx,1)))
    n = size(dist_mtx, 1)

    serv_time = zeros(n)
    dem = ones(n)
    vehicleCapacity = n
    maximum_number_of_vehicles = 1

    c_dist_mtx = Matrix(dist_mtx')
    for i in 1:n
        c_dist_mtx[i, i] = 0.0
    end

    # need to input dist_mtx' instead of dist_mtx
    # Julia: column-first indexing
    # C: row-first indexing
    return c_api_solve_cvrp_dist_mtx(n, x_coords, y_coords, c_dist_mtx, serv_time, dem, vehicleCapacity, maximum_number_of_vehicles, parameters, verbose)
end
function solve_tsp(x::Vector, y:: Vector, parameters=AlgorithmParameters(); verbose=true)
    n = length(x)
    serv_time = zeros(n)
    dem = ones(n)
    vehicleCapacity = n
    maximum_number_of_vehicles = 1

    return c_api_solve_cvrp(n, x, y, serv_time, dem, vehicleCapacity, maximum_number_of_vehicles, parameters, verbose)
end


text:
using Turing

# test case 1: stack allocated objects are deep copied.
function f()
  t = 0;
  while true
    produce(t)
    t = 1 + t
  end
end

t = Task(f)

consume(t); consume(t)
a = copy(t);
consume(a); consume(a)

# test case 2: heap allocated objects are shallowly copied.

function f2()
  t = [0 1 2];
  while true
    #println(pointer_from_objref(t));
    produce(t[1])
    t[1] = 1 + t[1]
  end
end

t = Task(f2)

consume(t); consume(t)
a = copy(t);
consume(a); consume(a)

# more: add code in copy() to handle invalid cases for cloning tasks.

function f3()
  t = [0];
  o = (x) -> x + 1;  # not heap allocated?
  while true
    produce(t[1])
    t[1] = 1 + t[1]
  end
  return o
end

t = Task(f3)

consume(t); consume(t);
a = copy(t);
consume(a); consume(a)


text:
abstract type SphericalHarmonicPlan{T} end

function *(P::SphericalHarmonicPlan, X::AbstractMatrix)
    mul!(zero(X), P, X)
end

if VERSION < v"0.7-"
    function \(P::SphericalHarmonicPlan, X::AbstractMatrix)
        At_mul_B!(zero(X), P, X)
    end
else
    function \(P::SphericalHarmonicPlan, X::AbstractMatrix)
        mul!(zero(X), transpose(P), X)
    end
end

include("sphfunctions.jl")
include("slowplan.jl")
include("Butterfly.jl")
include("fastplan.jl")
include("thinplan.jl")
include("synthesisanalysis.jl")
include("vectorfield.jl")

function plan_sph2fourier(A::AbstractMatrix; opts...)
    M, N = size(A)
    if M ≤ 1023
        SlowSphericalHarmonicPlan(A)
    else
        ThinSphericalHarmonicPlan(A; opts...)
    end
end

sph2fourier(A::AbstractMatrix; opts...) = plan_sph2fourier(A; opts...)*A
fourier2sph(A::AbstractMatrix; opts...) = plan_sph2fourier(A; opts...)\A

"""
Computes the bivariate Fourier series given by the spherical harmonic expansion:

```math
{\\rm SHT} : \\sum_{\\ell=0}^n\\sum_{m=-\\ell}^{\\ell} f_{\\ell}^m Y_{\\ell}^m(\\theta,\\varphi) \\to \\sum_{\\ell=0}^n\\sum_{m=-n}^{n} g_{\\ell}^m \\left\\{\\begin{array}{ccc}\\cos\\ell\\theta & {\\rm for} & m{\\rm~even}\\\\ \\sin(\\ell+1)\\theta & {\\rm for} & m{\\rm~odd}\\end{array}\\right\\}\\times \\sqrt{\\frac{2-\\delta_{m,0}}{2\\pi}} \\left\\{\\begin{array}{ccc} \\cos m\\varphi & {\\rm for} & m \\ge 0,\\\\ \\sin(-m\\varphi) & {\\rm for} & m < 0.\\end{array}\\right.
```

The spherical harmonic expansion coefficients are organized as follows:

```math
F = \\begin{pmatrix}
f_0^0 & f_1^{-1} & f_1^1 & f_2^{-2} & f_2^2 & \\cdots & f_n^{-n} & f_n^n\\\\
f_1^0 & f_2^{-1} & f_2^1 & f_3^{-2} & f_3^2 & \\cdots & 0 & 0\\\\
\\vdots & \\vdots & \\vdots &  \\vdots &  \\vdots & \\ddots & \\vdots & \\vdots\\\\
f_{n-2}^0 & f_{n-1}^{-1} & f_{n-1}^1 & f_n^{-2} & f_n^2 &  & \\vdots & \\vdots\\\\
f_{n-1}^0 & f_n^{-1} & f_n^1 & 0 & 0 & \\cdots & 0 & 0\\\\
f_n^0 & 0 & 0 & 0 & 0 & \\cdots & 0 & 0\\\\
\\end{pmatrix},
```

and the Fourier coefficients are organized similarly:

```math
G = \\begin{pmatrix}
g_0^0 & g_0^{-1} & g_0^1 & \\cdots & g_0^{-n} & g_0^n\\\\
g_1^0 & g_1^{-1} & g_1^1 & \\cdots & g_1^{-n} & g_1^n\\\\
\\vdots & \\vdots & \\vdots & \\ddots & \\vdots & \\vdots\\\\
g_{n-1}^0 & g_{n-1}^{-1} & g_{n-1}^1& \\cdots & g_{n-1}^{-n} & g_{n-1}^n\\\\
g_n^0 & 0 & 0 & \\cdots & g_n^{-n} & g_n^n\\\\
\\end{pmatrix}.
```
"""
sph2fourier(::AbstractMatrix; opts...)

"""
Computes the spherical harmonic expansion given by the bivariate Fourier series:

```math
{\\rm iSHT} : \\sum_{\\ell=0}^n\\sum_{m=-n}^{n} g_{\\ell}^m \\left\\{\\begin{array}{ccc}\\cos\\ell\\theta & {\\rm for} & m{\\rm~even}\\\\ \\sin(\\ell+1)\\theta & {\\rm for} & m{\\rm~odd}\\end{array}\\right\\}\\times \\sqrt{\\frac{2-\\delta_{m,0}}{2\\pi}} \\left\\{\\begin{array}{ccc} \\cos m\\varphi & {\\rm for} & m \\ge 0,\\\\ \\sin(-m\\varphi) & {\\rm for} & m < 0,\\end{array}\\right. \\to \\sum_{\\ell=0}^n\\sum_{m=-\\ell}^{\\ell} f_{\\ell}^m Y_{\\ell}^m(\\theta,\\varphi).
```

The spherical harmonic expansion coefficients are organized as follows:

```math
F = \\begin{pmatrix}
f_0^0 & f_1^{-1} & f_1^1 & f_2^{-2} & f_2^2 & \\cdots & f_n^{-n} & f_n^n\\\\
f_1^0 & f_2^{-1} & f_2^1 & f_3^{-2} & f_3^2 & \\cdots & 0 & 0\\\\
\\vdots & \\vdots & \\vdots &  \\vdots &  \\vdots & \\ddots & \\vdots & \\vdots\\\\
f_{n-2}^0 & f_{n-1}^{-1} & f_{n-1}^1 & f_n^{-2} & f_n^2 &  & \\vdots & \\vdots\\\\
f_{n-1}^0 & f_n^{-1} & f_n^1 & 0 & 0 & \\cdots & 0 & 0\\\\
f_n^0 & 0 & 0 & 0 & 0 & \\cdots & 0 & 0\\\\
\\end{pmatrix},
```

and the Fourier coefficients are organized similarly:

```math
G = \\begin{pmatrix}
g_0^0 & g_0^{-1} & g_0^1 & \\cdots & g_0^{-n} & g_0^n\\\\
g_1^0 & g_1^{-1} & g_1^1 & \\cdots & g_1^{-n} & g_1^n\\\\
\\vdots & \\vdots & \\vdots & \\ddots & \\vdots & \\vdots\\\\
g_{n-1}^0 & g_{n-1}^{-1} & g_{n-1}^1& \\cdots & g_{n-1}^{-n} & g_{n-1}^n\\\\
g_n^0 & 0 & 0 & \\cdots & g_n^{-n} & g_n^n\\\\
\\end{pmatrix}.
```
"""
fourier2sph(::AbstractMatrix; opts...)

"""
Pre-computes the spherical harmonic transform.
"""
plan_sph2fourier(::AbstractMatrix; opts...)


text:
include(pwd()*"/src/SFFM.jl")
using Plots

## define 1D wave equation model
T = [0.0]
C = [1]

rfun(x) = x.*0
Rfun(x) = r(x)

r = (
    r = function (x)
        rfun(x)
    end,
    R = function (x)
        Rfun(x)
    end
)

bounds = [0 12; -Inf Inf]
model = SFFM.Model( T, C, r, Bounds = bounds)

## approximation stuff
order = 15 # must be odd
Δ = 1
me = SFFM.MakeME(SFFM.CMEParams[order], mean = Δ)
nodes = collect(0:Δ:12)
mesh = SFFM.DGMesh(model, nodes, order)

## Initial condition
λ = -0.25# me.S[1,1]
f(x::Real) = sin(3*pi*x)#-λ*exp(λ*x)/(1-exp(λ*10))#1/12 #

## approximate integral of A(x)f(x)dx on [xₖ, xₖ₊₁]
nevals = 10000
params = SFFM.CMEParams[order]

α = zeros(order)
α[1] = params["c"]
a = params["a"]
b = params["b"]
ω =  params["omega"]
for k in 1:params["n"]
    kω = k*ω
    α[2*k] = (1/2)*( a[k]*(1+kω) - b[k]*(1-kω) )/(1+kω^2)
    α[2*k+1] = (1/2)*( a[k]*(1-kω) + b[k]*(1+kω) )/(1+kω^2)
end

A = zeros(Float64,SFFM.NIntervals(mesh),order)
for m in 1:SFFM.NIntervals(mesh)
    aInt = zeros(Float64,order)

    orbit_LHS = α
    orbit_RHS = zeros(order)
    x = mesh.Nodes[m]
    h = Δ/nevals
    for n in 2:nevals
        orbit_RHS[1] = α[1]
        x += h
        for k in 1:params["n"]
            kωx = k*ω*((n-1)*h).*-me.S[1,1]
            idx = 2*k
            idx2 = idx+1
            temp_cos = cos(kωx)
            temp_sin = sin(kωx)
            orbit_RHS[idx] = α[idx]*temp_cos + α[idx2]*temp_sin
            orbit_RHS[idx2] = -α[idx]*temp_sin + α[idx2]*temp_cos
        end
        orbit_RHS = orbit_RHS./sum(orbit_RHS)
        orbit = (orbit_LHS+orbit_RHS)./2

        aInt += orbit.*h.*(f(x)+f(x-h))/2

        orbit_LHS = copy(orbit_RHS)
    end
    A[m,:] = aInt
end

xvec = collect(mesh.Nodes[1]:0.05:mesh.Nodes[end])
p = plot(xvec,f.(xvec),label = "f")

xvec = collect(0:0.05:Δ)
for k in 1:SFFM.NIntervals(mesh)
    mek = SFFM.ME(Matrix(A[k,:]'),me.S,me.s,D = me.D)
    thePDF = SFFM.pdf(mek,xvec)
    thePDFplusΔ = SFFM.pdf(mek,xvec.+Δ)[end:-1:1]
    plot!(mesh.Nodes[k+1] .- xvec, thePDF+thePDFplusΔ, label = "A")
end
p

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
export autocor_sol, autocor_sol_corr0t, thermalization, thermalization_corr0t, get_statistics,
        distribution_avg,distribution_avg2,distribution_corr0t, get_corrections

##################################################

using StatsBase,Statistics,Measurements,Bootstrap

"""
plot AUTO-CORRELATION from solution of DifferentialEquation.jl solver
"""
function autocor_sol(sol,TT;therm=1, trajectories=[1],inxs::Array{Int64}=[1])
    fig = plot(legend=false)
    for tr in trajectories
        D = transpose(hcat(sol[tr].u[therm:TT:end]...))
        plot!(fig,autocor(D),label=string("Tr: ",tr))
    end
    display(fig)
end

function autocor_sol_corr0t(sol,TT;therm=1, trajectories=[1])
    t_steps = div(size(sol[1])[1],2)
    fig = plot(legend=false)
    for tr in trajectories
        SMu = hcat(sol[tr].u[therm:TT:end]...)
        corr0t_Re_1 = SMu[1,:] .* SMu[1,:] .- SMu[t_steps+1,:] .* SMu[t_steps+1,:]

        plot!(fig,mean(autocor(corr0t_Re_1),dims=2))
    end
    display(fig)
end

function thermalization(sol,TT::Integer,therm::Integer; trajectory=1)
    thermalization(sol,TT,therm,collect(1:size(sol[1])[1]); trajectory=trajectory)
end

function thermalization(sol,TT::Integer,therm::Integer, inxs::Array{Int64}; trajectory=1)
    D = hcat(sol[trajectory].u[therm:TT:end]...)[inxs,:]
    println("Size of solution after thermalization: ", size(D))
    plot(transpose(D),legend=true)#,xlim=[0,100])
end

function thermalization_corr0t(sol,TT::Integer,therm::Integer,inxs::Array{Int64})
    t_steps = div(size(sol[1])[1],2)
    SMu = hcat(sol[trajectory].u[therm:TT:end]...)
    corr0t_Re_1 = SMu[1,:] .* SMu[1,:] .- SMu[t_steps+1,:] .* SMu[t_steps+1,:]
    corr0t_Re_2 = SMu[1,:] .* SMu[4,:] .- SMu[t_steps + 1 ,:] .* SMu[t_steps+4,:]
    corr0t_Re_3 = SMu[1,:] .* SMu[floor(Int64,t_steps/3),:] .- SMu[t_steps + 1 ,:] .* SMu[floor(Int64,4*t_steps/3),:]
    corr0t_Re_4 = SMu[1,:] .* SMu[floor(Int64,t_steps/2),:] .- SMu[t_steps + 1 ,:] .* SMu[floor(Int64,3*t_steps/2),:]
    plot(corr0t_Re_1)#,legend=false)#,xlim=[0,10])
    plot!(corr0t_Re_2)
    plot!(corr0t_Re_3)
    plot!(corr0t_Re_4)#,xlim=[0,100])
end

function distribution_avg(sol,TT::Integer,therm::Integer, inxs::Array{Int64})
    N_Tr = size(sol)[end]

    D = hcat(sol[1].u[therm:TT:end]...)[inxs,:]
    for tr in 2:N_Tr
        D = hcat(D,hcat(sol[tr].u[therm:TT:end]...)[inxs,:])
    end
    histogram(transpose(D),bins=100,normalize=true)
end

function distribution_avg2(sol,TT::Integer,therm::Integer, inxs::Array{Int64}; Re=true, Im=false)
    t_steps = div(size(sol[1])[1],2)
    N_Tr = size(sol)[end]

    D = hcat(sol[1].u[therm:TT:end]...)[inxs,:]
    D2 =  hcat(sol[1].u[therm:TT:end]...)[t_steps .+ inxs,:]
    for tr in 2:N_Tr
        D  = hcat(D, hcat(sol[tr].u[therm:TT:end]...)[inxs,:])
        D2 = hcat(D2,hcat(sol[tr].u[therm:TT:end]...)[t_steps .+ inxs,:])
    end

    if Re
        V = transpose(D.^2 .+ D2.^2)
    elseif Im
        V = transpose(2*D .* D2)
    end

    bin_size = 0.1
    bins = floor(Int64,(maximum(V)-minimum(V))/bin_size)

    #histogram(V,bins=bins,normalize=true, xlim=[-3,3], seriestype=:stephist)
    StatsPlots.density(V,xlim=[-2,2])
end

function distribution_corr0t(sol,TT::Integer,therm::Integer, inxs::Array{Int64})#, inxs::Array{Int64})
    t_steps = div(size(sol[1])[1],2)
    N_Tr = size(sol)[end]

    D = hcat(sol[1].u[therm:TT:end]...)[inxs,:]
    D2 =  hcat(sol[1].u[therm:TT:end]...)[t_steps .+ inxs,:]
    for tr in 2:N_Tr
        D  = hcat(D, hcat(sol[tr].u[therm:TT:end]...)[inxs,:])
        D2 = hcat(D2,hcat(sol[tr].u[therm:TT:end]...)[t_steps .+ inxs,:])
    end
    histogram(D,bins=70,normalize=true,legend=false)
end

## Statistics
function get_statistics(sol;TT=1,therm=1,thermal = 0)

    N_Tr = length(sol)
    t_steps = div(size(sol[1])[1],2)
    n_boot  = 50

    avg_Re = zeros(Float64,N_Tr,t_steps)
    avg_Im = zeros(Float64,N_Tr,t_steps)
    avg2_Re = zeros(Float64,N_Tr,t_steps)
    avg2_Im = zeros(Float64,N_Tr,t_steps)
    corr0t_Re = zeros(Float64,N_Tr,t_steps)
    corr0t_Im = zeros(Float64,N_Tr,t_steps)
    Gpm_Re = zeros(Float64,N_Tr,t_steps)
    Gpm_Im = zeros(Float64,N_Tr,t_steps)
    Gmp_Re = zeros(Float64,N_Tr,t_steps)
    Gmp_Im = zeros(Float64,N_Tr,t_steps)

    Threads.@threads for tr in 1:N_Tr
    #for tr in 1:N_Tr
        SMu = hcat(sol[tr].u[therm:TT:end]...)

        for i in 1:t_steps-thermal
            avg_Re[tr,i] = mean(SMu[i,:])
            avg_Im[tr,i] = mean(SMu[t_steps+i,:])

            avg2_Re[tr,i] = mean(SMu[i,:].^2 - SMu[t_steps+i,:].^2)
            avg2_Im[tr,i] = mean(2*SMu[i,:] .* SMu[t_steps+i,:])

            corr0t_Re_tmp = SMu[1,:] .* SMu[i,:] .- SMu[t_steps+1,:] .* SMu[t_steps+i,:]
            corr0t_Re[tr,i] = mean(corr0t_Re_tmp)

            corr0t_Im_tmp =  SMu[1,:] .* SMu[t_steps+i,:] .+ SMu[i,:] .* SMu[t_steps+1,:]
            corr0t_Im[tr,i] = mean(corr0t_Im_tmp)

            Gpm_Re_tmp = SMu[1,:] .* SMu[t_steps-i-thermal+1,:] .- SMu[t_steps+1,:] .* SMu[2*t_steps-i-thermal+1,:]
            Gpm_Re[tr,i] = mean(Gpm_Re_tmp)

            Gpm_Im_tmp =  SMu[1,:] .* SMu[2t_steps-i-thermal+1,:] .+ SMu[t_steps-i-thermal+1,:] .* SMu[t_steps+1,:]
            Gpm_Im[tr,i] = mean(Gpm_Im_tmp)

            Gmp_Re_tmp = SMu[t_steps-thermal,:] .* SMu[i,:] .- SMu[2*t_steps-thermal,:] .* SMu[t_steps+i,:]
            Gmp_Re[tr,i] = mean(Gmp_Re_tmp)

            Gmp_Im_tmp =  SMu[t_steps-thermal,:] .* SMu[t_steps+i,:] .+ SMu[i,:] .* SMu[2*t_steps-thermal,:]
            Gmp_Im[tr,i] = mean(Gmp_Im_tmp)

        end
    end

    d = 1

    mean_d(x) = mean(x,dims=d)
    avg_Re_bts = bootstrap(mean_d, avg_Re, BalancedSampling(n_boot))
    avg_Im_bts = bootstrap(mean_d, avg_Im, BalancedSampling(n_boot))
    avg2_Re_bts = bootstrap(mean_d, avg2_Re, BalancedSampling(n_boot))
    avg2_Im_bts = bootstrap(mean_d, avg2_Im, BalancedSampling(n_boot))
    corr0t_Re_bts = bootstrap(mean_d, corr0t_Re, BalancedSampling(n_boot))
    corr0t_Im_bts = bootstrap(mean_d, corr0t_Im, BalancedSampling(n_boot))
    Gpm_Re_bts = bootstrap(mean_d, Gpm_Re, BalancedSampling(n_boot))
    Gpm_Im_bts = bootstrap(mean_d, Gpm_Im, BalancedSampling(n_boot))
    Gmp_Re_bts = bootstrap(mean_d, Gmp_Re, BalancedSampling(n_boot))
    Gmp_Im_bts = bootstrap(mean_d, Gmp_Im, BalancedSampling(n_boot))

    std_err(x) = std(x,dims=d,corrected=true) ./ sqrt(size(x)[2])
    avg_Re_bts_std = bootstrap(std_err, avg_Re, BalancedSampling(n_boot))
    avg_Im_bts_std = bootstrap(std_err, avg_Im, BalancedSampling(n_boot))
    avg2_Re_bts_std = bootstrap(std_err, avg2_Re, BasicSampling(n_boot))
    avg2_Im_bts_std = bootstrap(std_err, avg2_Im, BalancedSampling(n_boot))
    corr0t_Re_bts_std = bootstrap(std_err, corr0t_Re, BalancedSampling(n_boot))
    corr0t_Im_bts_std = bootstrap(std_err, corr0t_Im, BalancedSampling(n_boot))
    Gpm_Re_bts_std = bootstrap(std_err, Gpm_Re, BalancedSampling(n_boot))
    Gpm_Im_bts_std = bootstrap(std_err, Gpm_Im, BalancedSampling(n_boot))
    Gmp_Re_bts_std = bootstrap(std_err, Gmp_Re, BalancedSampling(n_boot))
    Gmp_Im_bts_std = bootstrap(std_err, Gmp_Im, BalancedSampling(n_boot))

    avg_Re = collect(original(avg_Re_bts)) .± original(avg_Re_bts_std)
    avg_Im = collect(original(avg_Im_bts)) .± original(avg_Im_bts_std)
    avg2_Re = collect(original(avg2_Re_bts)) .± (original(avg2_Re_bts_std))
    avg2_Im = collect(original(avg2_Im_bts)) .± original(avg2_Im_bts_std)
    corr0t_Re = (collect(original(corr0t_Re_bts)) .± original(corr0t_Re_bts_std)) .- avg_Re.*avg_Re[1]
    corr0t_Im = (collect(original(corr0t_Im_bts)) .± original(corr0t_Im_bts_std)) .- avg_Im.*avg_Im[1]
    Gpm_Re = (collect(original(Gpm_Re_bts)) .± original(Gpm_Re_bts_std)) .- avg_Re.*avg_Re[1]
    Gpm_Im = (collect(original(Gpm_Im_bts)) .± original(Gpm_Im_bts_std)) .- avg_Im.*avg_Im[1]
    Gmp_Re = (collect(original(Gmp_Re_bts)) .± original(Gmp_Re_bts_std)) .- avg_Re.*avg_Re[t_steps]
    Gmp_Im = (collect(original(Gmp_Im_bts)) .± original(Gmp_Im_bts_std)) .- avg_Im.*avg_Im[t_steps]

    return avg_Re, avg_Im, avg2_Re, avg2_Im, corr0t_Re, corr0t_Im, Gpm_Re, Gpm_Im, Gmp_Re, Gmp_Im  #twoPoint_Re, twoPoint_Im

end

function get_corrections(sol,params,args)

    N_Tr = length(sol)
    t_steps = div(size(sol[1])[1],2)

    du = zeros(Float64,t_steps*2)
    SS_Re_tmp = zeros(Float64,length(sol),t_steps)
    SS2_Re_tmp = zeros(Float64,length(sol),t_steps)
    SS_Im_tmp = zeros(Float64,length(sol),t_steps)
    SS2_Im_tmp = zeros(Float64,length(sol),t_steps)
    for tr in 1:N_Tr
        for i in 1:length(sol[tr][1,:])
            a_AHO_real(du,sol[tr][:,i],params,1)
            SS_Re_tmp[tr,:]  += du[1:t_steps].^2 .- du[t_steps+1:2*t_steps].^2
            SS2_Re_tmp[tr,:] +=(du[1:t_steps] .* sol[tr][1:t_steps,i]
                               .- du[t_steps+1:2*t_steps] .* sol[tr][t_steps+1:2*t_steps,i])
            SS_Im_tmp[tr,:] += 2*du[1:t_steps] .* du[t_steps+1:2*t_steps]
            SS2_Im_tmp[tr,:] += (du[1:t_steps] .* sol[tr][t_steps+1:2*t_steps,i]
                               .+ du[t_steps+1:2*t_steps] .* sol[tr][1:t_steps,i])
        end
    end

    SS_Re = transpose(mean(((args.dt/2)*(0.5-args.θ))^2 .* SS_Re_tmp ./ length(sol[1][1,:]),dims=1))
    SS2_Re = -transpose(mean((args.dt*(0.5-args.θ)) .* SS2_Re_tmp ./ length(sol[1][1,:]),dims=1))
    SS_Im = transpose(mean(((args.dt/2)*(0.5-args.θ))^2 .* SS_Im_tmp ./ length(sol[1][1,:]),dims=1))
    SS2_Im = -transpose(mean((args.dt*(0.5-args.θ)) .* SS2_Im_tmp ./ length(sol[1][1,:]),dims=1))

    return SS_Re, SS2_Re, SS_Im, SS2_Im
end


text:
using HMatrices
using LinearAlgebra
using Test
using StaticArrays
using HMatrices: RkMatrix, compression_ratio

@testset "RkMatrix" begin
    @testset "Scalar entries" begin
        m = 20
        n = 30
        r = 5
        A = rand(ComplexF64,m,r)
        B = rand(ComplexF64,n,r)
        R  = RkMatrix(A,B);
        Ra  = adjoint(R);
        M  = A*adjoint(B)
        Ma  = adjoint(M)

        ## basic tests
        @test size(R) == (m,n)
        @test rank(R) == r
        @test Matrix(R) ≈ M
        @test compression_ratio(R) ≈ m*n / (r*(m+n))
        @test R[:,5] ≈ M[:,5]
        @test Ra[:,5] ≈ Ma[:,5]
    end

    @testset "Matrix entries" begin
        m = 20
        n = 30
        r = 10
        T = SMatrix{3,3,ComplexF64,9}
        A = rand(T,m,r)
        B = rand(T,n,r)
        R  = RkMatrix(A,B)
        ## basic tests
        @test size(R) == (m,n)
        @test rank(R) == r
        @test compression_ratio(R) ≈ m*n / (r*(m+n))
    end
end

max_stars_repo_name: WaveProp/HMatrices

max_forks_repo_name: WaveProp/HMatrices


text:
####################
# value extraction #
####################

@inline extract_value!(out::DiffResult, ydual) = DiffBase.value!(value, out, ydual)
@inline extract_value!(out, ydual) = nothing

@inline function extract_value!(out, y, ydual)
    map!(value, y, ydual)
    copy_value!(out, y)
end

@inline copy_value!(out::DiffResult, y) = DiffBase.value!(out, y)
@inline copy_value!(out, y) = nothing

###################################
# vector mode function evaluation #
###################################

@generated function dualize(::F, x::SArray{S,V,D,N}) where {F,S,V,D,N}
    tag = Tag(F, x)
    dx = Expr(:tuple, [:(Dual{T}(x[$i], chunk, Val{$i}())) for i in 1:N]...)
    return quote
        chunk = Chunk{N}()
        T = typeof($tag)
        $(Expr(:meta, :inline))
        return SArray{S}($(dx))
    end
end

@inline vector_mode_dual_eval(f::F, x::SArray) where {F} = f(dualize(f, x))

function vector_mode_dual_eval(f::F, x, cfg::Union{JacobianConfig,GradientConfig}) where F
    xdual = cfg.duals
    seed!(xdual, x, cfg.seeds)
    return f(xdual)
end

function vector_mode_dual_eval(f!::F, y, x, cfg::JacobianConfig) where F
    ydual, xdual = cfg.duals
    seed!(xdual, x, cfg.seeds)
    seed!(ydual, y)
    f!(ydual, xdual)
    return ydual
end

##################################
# seed construction/manipulation #
##################################

@generated function construct_seeds(::Type{Partials{N,V}}) where {N,V}
    return Expr(:tuple, [:(single_seed(Partials{N,V}, Val{$i}())) for i in 1:N]...)
end

function seed!(duals::AbstractArray{Dual{T,V,N}}, x,
               seed::Partials{N,V} = zero(Partials{N,V})) where {T,V,N}
    for i in eachindex(duals)
        duals[i] = Dual{T,V,N}(x[i], seed)
    end
    return duals
end

function seed!(duals::AbstractArray{Dual{T,V,N}}, x,
               seeds::NTuple{N,Partials{N,V}}) where {T,V,N}
    for i in 1:N
        duals[i] = Dual{T,V,N}(x[i], seeds[i])
    end
    return duals
end

function seed!(duals::AbstractArray{Dual{T,V,N}}, x, index,
               seed::Partials{N,V} = zero(Partials{N,V})) where {T,V,N}
    offset = index - 1
    for i in 1:N
        j = i + offset
        duals[j] = Dual{T,V,N}(x[j], seed)
    end
    return duals
end

function seed!(duals::AbstractArray{Dual{T,V,N}}, x, index,
               seeds::NTuple{N,Partials{N,V}}, chunksize = N) where {T,V,N}
    offset = index - 1
    for i in 1:chunksize
        j = i + offset
        duals[j] = Dual{T,V,N}(x[j], seeds[i])
    end
    return duals
end


text:
function evaluate_definitions(model, simul::AxisArray{Tf,3}, params=model.calibration[:parameters]) where Tf

    p_ = SVector(params...)

    # @assert axisnames(simul) == (:N,:V,:T)
    T = length( simul[Axis{:T}].val )
    @assert simul[Axis{:T}].val == 1:T

    vars = cat(model.symbols[:exogenous], model.symbols[:states], model.symbols[:controls]; dims=1)

    sim = simul[Axis{:V}(vars)]

    T = length(sim[Axis{:T}].val)

    past = permutedims(sim[Axis{:T}([1;1:T-1])], [2,1,3])
    present = permutedims(sim[Axis{:T}(1:T)], [2,1,3])
    future = permutedims(sim[Axis{:T}([2:T;T])], [2,1,3])

    n_v,N,T = size(past)

    x_past = reshape(reinterpret(Point{n_v}, vec(past.data)), (T*N,))
    x_present = reshape(reinterpret(Point{n_v}, vec(present.data)), (T*N,))
    x_future = reshape(reinterpret(Point{n_v}, vec(future.data)), (T*N,))

    y_ = evaluate_definitions(model, x_past, x_present, x_future, p_)

    auxiliaries = [Dolang.arg_name(e) for e in keys(model.definitions)]
    n_y = length(auxiliaries)

    data = permutedims( reshape(reinterpret(Float64, vec(y_)), (n_y,N,T)), [2,1,3])

    array = AxisArray(copy(data), Axis{:N}(1:N), Axis{:V}(auxiliaries), Axis{:T}(1:T))

end

function evaluate_definitions(model, _simul::AxisArray{__T,2}, params=model.calibration[:parameters]) where __T

    p_ = SVector(params...)

    all_v = _simul[Axis{:V}].val
    all_t = _simul[Axis{:T}].val

    # TODO: why are we even doing that ?
    T = length(all_t)
    if all(all_t .== 0:(T-1))
        all_t = 1:T
    end
    if all(all_t .!= 1:T)
        msg = "Can only evaluate definitions if the T Axis goes from 1:T or 0:(T-1)"
        error(msg)
    end

    # expect dimensions to be (:V,:T) otherwise we transpose
    if axisnames(_simul) == (:T,:V)
        data = _simul.data'
    else
        data = _simul.data
    end

    simul = AxisArray(data, Axis{:V}(all_v), Axis{:T}(all_t))

    vars = cat(model.symbols[:exogenous], model.symbols[:states], model.symbols[:controls]; dims=1)

    sim = simul[Axis{:V}(vars)]

    T = length(sim[Axis{:T}].val)

    past = sim[Axis{:T}([1; 1:T-1])]
    present = sim[Axis{:T}(1:T)]
    future = sim[Axis{:T}([2:T; T])]

    n_v = size(past, 1)

    x_past = reshape(reinterpret(Point{n_v}, vec(past.data)), (T,))
    x_present = reshape(reinterpret(Point{n_v}, vec(present.data)), (T,))
    x_future = reshape(reinterpret(Point{n_v}, vec(future.data)), (T,))

    y_ = evaluate_definitions(model, x_past, x_present, x_future, p_)

    auxiliaries = [Dolang.arg_name(e) for e in keys(model.definitions)]
    n_y = length(auxiliaries)

    data = reshape(reinterpret(Float64, vec(y_)), (n_y, T))

    array = AxisArray(copy(data), Axis{:V}(auxiliaries), simul[Axis{:T}])

end

########################################################################
# For Float

function simulate(model::AbstractModel, dr::AbstractDecisionRule,
                  driving_process::AbstractArray{Float64,3}; s0::AbstractVector=model.calibration[:states])

    # driving_process: (ne, N, T)
    driving_process = convert(Array{Float64,3}, driving_process) # in case arg is an axisarray

    # extract data from model
    calib = model.calibration
    params = SVector(calib[:parameters]...)

    N = size(driving_process, 2)
    T = size(driving_process, 3)
    epsilons = permutedims(driving_process, [2, 1, 3]) # (N, ne, T)
    # calculate initial controls using decision rule

    x0 = dr(epsilons[1, :, 1], s0)

    # get number of states and controls
    ns = length(s0)
    nx = length(x0)
    nsx = nx+ns

    s_simul = Array{Float64}(undef, N, ns, T)
    x_simul = Array{Float64}(undef, N, nx, T)
    for i in 1:N
      s_simul[i, :, 1] = s0
      x_simul[i, :, 1] = x0
    end

    for t in 1:T
        s = copy(to_LOP(s_simul[:, :, t]))
        m = copy(to_LOP(epsilons[:, :, t]))
        x = dr(m, s)
        x_simul[:, :, t] = from_LOP(x)
        if t < T
          M = copy(to_LOP(epsilons[:, :, t+1]))
          ss = transition(model, m, s, x, M, params)
          s_simul[:,:,t+1] = from_LOP(ss)
        end
    end
    sim = cat(epsilons, s_simul, x_simul; dims=2)::Array{Float64,3}

    Ac = cat(model.symbols[:exogenous], model.symbols[:states], model.symbols[:controls]; dims=1)
    ll = [Symbol(i) for i in Ac]

    sim_aa = AxisArray(sim, Axis{:N}(1:N), Axis{:V}(ll), Axis{:T}(1:T))

    if length(model.definitions)>0
        sim_def= evaluate_definitions(model, sim_aa, model.calibration[:parameters])
        return merge(sim_aa,sim_def)
    else
        return sim_aa
    end

end

###########################################################################
# For Int (MC)
function simulate(model::AbstractModel, dr::AbstractDecisionRule,
                  driving_process::AbstractMatrix{Int}, dp_process::DiscreteMarkovProcess;
                  s0::AbstractVector=model.calibration[:states])

    # driving_process: (ne, N, T)

    # extract data from model
    # calib = model.calibration
    params = model.calibration[:parameters]

    N = size(driving_process, 2)
    T = size(driving_process, 1)
    epsilons = zeros(Int, N, 1, T)
    for _n in 1:N
        epsilons[_n, 1, :] = driving_process[:, _n]
    end

    # calculate initial controls using decision rule
    x0 = dr(epsilons[1, :, 1], s0)
    # get number of states and controls
    ns = length(s0)
    nx = length(x0)
    nsx = nx+ns
    nm = length(model.calibration[:exogenous])

    s_simul = Array{Float64}(undef, N, ns, T)
    x_simul = Array{Float64}(undef, N, nx, T)
    for i in 1:N
        s_simul[i, :, 1] = s0
        x_simul[i, :, 1] = x0
    end

    for t in 1:T
        s = s_simul[:, :, t]
        m = epsilons[:, :, t]

        m_ind=cat(m; dims=1)[:, 1]
        m_val= dp_process.values[m_ind, :]
        x = dr(m_ind, s)
        x_simul[:, :, t] = x
        if t < T
            M = epsilons[:, :, t+1]
            M_cat = cat(M, dims=1)[:, 1]
            M_val = dp_process.values[M_cat, :]
            s_simul[:,:,t+1] = transition(model, m_val, s, x, M_val, params)
        end
    end

    epsilons_values = zeros(N,nm,T)
    for n=1:N
      eps = dp_process.values[epsilons[n,:,:],:]
      epsilons_values[n,:,:] = permutedims(eps, [1, 3, 2])
    end

    sim = cat(epsilons, epsilons_values, s_simul, x_simul; dims=2)::Array{Float64,3}

    model_sym = :mc_process
    Ac = cat(model_sym, model.symbols[:exogenous], model.symbols[:states], model.symbols[:controls], dims=1)
    ll = [Symbol(i) for i in Ac]
    sim_aa = AxisArray(sim, Axis{:N}(1:N), Axis{:V}(ll), Axis{:T}(1:T))

    if length(model.definitions)>0
        sim_def=evaluate_definitions(model, sim_aa, model.calibration[:parameters])
        return merge(sim_aa,sim_def)
    else
        return sim_aa
    end

end

"""
This function simulates a model given a decision rule.

# Arguments
* `model::NumericModel`: Model object that describes the current model environment.
* `dr`: Solved decision rule.
* optional:
  * `driving process` simulated series of a model's exogenous process.
  If  `driving process` is not provided, then:
  * `dprocess`: exogenous processes; default: model.exogenous.
  * `s0::ListOfPoints`: List of initial state variable values; default: model.calibration[:states]
  * `N`: number of simulations; default: 1.
  * `T`: number of periods of simulations; default: 40.
  * `i0`: initial state of a discretized exogenous process; default: default_index(dprocess).
  or
  * `m0::ListOfPoints`: List of initial values of a continuous exogenous process; default: model.calibration[:exogenous].
# Returns
* `simulate`: simulated time series.
"""
function simulate(model::AbstractModel, dr::AbstractDecisionRule, dprocess::DiscreteMarkovProcess;
                  i0::Int=default_index(dprocess), s0::AbstractVector=model.calibration[:states], N::Int=1, T::Int=40)
    driving_process = simulate(dprocess; N=N, T=T, i0=i0)
    return simulate(model, dr, driving_process, dprocess; s0=s0)
end

function simulate(model::AbstractModel, dr::AbstractDecisionRule, dprocess::ContinuousProcess;
                  m0::AbstractVector=model.calibration[:exogenous], s0::AbstractVector=model.calibration[:states],
                  N::Int=1, T::Int=40)
    driving_process = simulate(dprocess, m0; N=N, T=T)
    return simulate(model, dr, driving_process; s0 = s0)
end

function simulate(model::AbstractModel, dr::AbstractDecisionRule; kwargs...)
    dprocess = model.exogenous
    return simulate(model, dr, dprocess; kwargs...)
end

"""
This is the one we document.
"""
simulate
################################################################################
## Impulse response functions

function response(model::AbstractModel,  dr::AbstractDecisionRule;
                  s0::AbstractVector, e1::AbstractVector, T::Int=40)
    m_sim = response(model.exogenous, e1; T=T)
    m_simul = reshape(m_sim, size(m_sim,1), 1, size(m_sim,2))
    sim = simulate(model, dr, m_simul; s0=s0)
    sim[1, :, :] # This is now an AxisArray which seems just fine !
end

function response(model::AbstractModel,  dr::AbstractDecisionRule,
    s0::AbstractVector, e1::AbstractVector; T::Int=40)

    response(model,  dr; s0=s0, e1=e1, T=T)
end

function response(model::AbstractModel,  dr::AbstractDecisionRule,
                  e1::AbstractVector; T::Int=40)
    s0 = model.calibration[:states]
    response(model, dr; s0=s0, e1=e1, T=T)
end

function response(model::AbstractModel,  dr::AbstractDecisionRule,
                  s0::AbstractVector, shock_name::Symbol; T::Int=40)
    index_s = something(findfirst(isequal(shock_name), model.symbols[:exogenous]), 0)
    # e1 = zeros(length(model.exogenous.mu))
    e1 = zeros(length(model.calibration[:exogenous]))
    Impulse = sqrt(diag(model.exogenous.Σ)[index_s])
    e1[index_s] = Impulse
    response(model, dr; s0=s0, e1=e1, T=T)
end

function response(model::AbstractModel,  dr::AbstractDecisionRule,
                  shock_name::Symbol; kwargs...)
    s0 = model.calibration[:states]
    response(model, dr, s0, shock_name;  kwargs...)
end

"""
Function "response" computes the impulse response functions with several major options:
- the user can provide a vector with the first values of the model's exogenous processes, e1.
- the user can provide a name of the shock of interest and the size of the shock_name.
- the user can provide only a name of the shock of interest. The size of the shock is assumed to be a one standard deviation given in the yaml file.

# Arguments
* `model::NumericModel`: Model object that describes the current model environment.
* `dr`: Solved decision rule.
* `e1::ListOfPoints`: List of initial model's exogenous processes values.
* If e1 is not provided, then:
  * `shock_name`: the name of the shock of interest.
* optional:
  * `Impulse`: the size of the shock; default: one standard deviation.
  * `s0::ListOfPoints`: List of initial state variable values; default: model.calibration[:states]
# Returns
* `response`: Impulse response function.
"""
function response(model::AbstractModel,  dr::AbstractDecisionRule,
                  s0::AbstractVector, shock_name::Symbol, Impulse::Float64; T::Int=40)
    index_s = something(findfirst(isequal(shock_name), model.symbols[:exogenous]), 0)
    e1 = zeros(length(model.calibration[:exogenous]))
    e1[index_s] = Impulse
    s0 = model.calibration[:states]
    response(model, dr; s0=s0, e1=e1, T=T)
end

function response(model::AbstractModel,  dr::AbstractDecisionRule,
                  shock_name::Symbol, Impulse::Float64; T::Int=40)
    s0 = model.calibration[:states]
    response(model, dr, s0, shock_name, Impulse; T=T)
end

"""
Function "tabulate" produces a 2-dimensional AxisArray{Float64,2,...}  with 2 axes : V (containing all the variables of the model) and the name of the state variable chosen in input.
"""
function tabulate(model::AbstractModel, dr::AbstractDecisionRule, state::Symbol,
                  bounds::Array{Float64,1}, s0::AbstractVector,
                  m0::Union{Int,AbstractVector}; n_steps::Int=100)

    index = findfirst(isequal(state), model.symbols[:states])
    Svalues = range(bounds[1], stop=bounds[2]; length=n_steps)
    svec = vcat([e' for e in fill(s0, n_steps)]...)
    svec[:, index] = Svalues
    s = copy(to_LOP(svec))

    if isa(dr.grid_exo, UnstructuredGrid)
                model_sym = cat([:mc_process], model.symbols[:exogenous]; dims=1)
                xvec = dr(m0, s)
                mv = node(dr.grid_exo, m0)
                l1 = [SVector( m0, mv..., e..., f...) for (e,f) in zip(s,xvec)]
    else
                xvec = dr(m0, s)
                # l1 = concat.(mm, svec, xvec)
                l1 = [SVector( m0..., e..., f...) for (e,f) in zip(s,xvec)]
                model_sym = model.symbols[:exogenous]
    end
    tb = copy(from_LOP(l1))

    l2 = cat(model_sym , model.symbols[:states], model.symbols[:controls]; dims=1)
    tab_AA = AxisArray(reshape(tb,1,size(tb)...), Axis{:T}(1:1), Axis{:N}(1:n_steps), Axis{:V}(l2))

    ## add definitions
    tab_AAA = permutedims(tab_AA, [2,3,1])

    if length(model.definitions)>0
        tab_defs = evaluate_definitions(model, tab_AAA)
        tab_ = merge(tab_AAA, tab_defs)[Axis{:T}(1)]
    else
        tab_ = tab_AAA[Axis{:T}(1)]
    end
    #change axis names
    res = AxisArray(tab_.data, Axis{state}(tb[:, index+1]), tab_[Axis{:V}])
    res' # so that we can index it directly
end

"""
Function "tabulate" produces a 2-dimensional AxisArray{Float64,2,...}  with 2 axes : V (containing all the variables of the model) and the name of the state variable chosen in input.
"""
function tabulate_old(model::AbstractModel, dr::AbstractDecisionRule, state::Symbol,
                  bounds::Array{Float64,1}, s0::AbstractVector,
                  m0::Union{Int,AbstractVector}; n_steps::Int=100)

    index = findfirst(isequal(state), model.symbols[:states])
    Svalues = range(bounds[1], stop=bounds[2], length=n_steps)
    svec = vcat([e' for e in fill(s0, n_steps)]...)
    svec[:, index] = Svalues

    if isa(dr.grid_exo, UnstructuredGrid)
        model_sym = :mc_process
    else
        xvec = dr(m0, svec)
        mm = vcat([e' for e in fill(m0, n_steps)]...)
        l1 = [mm, svec, xvec]
        tb = hcat([e' for e in l1']...)
        model_sym = model.symbols[:exogenous]
    end

    l2 = cat(model_sym , model.symbols[:states], model.symbols[:controls]; dims=1)
    tab_AA = AxisArray(reshape(tb,1,size(tb)...), Axis{:T}(1:1), Axis{:N}(1:n_steps), Axis{:V}(l2))

    ## add definitions
    tab_AAA = permutedims(tab_AA, [2,3,1])

    if length(model.definitions)>0
        tab_defs = evaluate_definitions(model, tab_AAA)
        tab_ = merge(tab_AAA, tab_defs)[Axis{:T}(1)]
    else
        tab_ = tab_AAA[Axis{:T}(1)]
    end
    #change axis names
    res = AxisArray(tab_.data, Axis{state}(tb[:, index+1]), tab_[Axis{:V}])
    res' # so that we can index it directly
end

function tabulate(model::AbstractModel, dr::AbstractDecisionRule, state::Symbol,
                  s0::AbstractVector, m0::Union{Int,AbstractVector};  n_steps=100)
    index = findfirst(isequal(state), model.symbols[:states])
    bounds = [model.domain.endo.min[index], model.domain.endo.max[index]]
    tabulate(model, dr, state, bounds, s0, m0;  n_steps=100)
end

function tabulate(model::AbstractModel, dr::AbstractDecisionRule, state::Symbol,
                  s0::AbstractVector; n_steps=100)

     if isa(dr.grid_exo, UnstructuredGrid)
        m0 = 1
    else
        m0 = SVector(model.calibration[:exogenous]...)
    end
    tabulate(model, dr, state, s0, m0;  n_steps=100)
end

function tabulate(model::AbstractModel, dr::AbstractDecisionRule, state::Symbol; n_steps=100)
    s0 = model.calibration[:states]
    tabulate(model, dr, state, s0;  n_steps=100)
end

"""
That's how we tabulate functions.
"""


text:
function ImGui_ImplOpenGL3_CreateDeviceObjects(ctx::Context)
    vertex_shader_glsl_130 = """
        #version $(ctx.GlslVersion[])
        uniform mat4 ProjMtx;
        in vec2 Position;
        in vec2 UV;
        in vec4 Color;
        out vec2 Frag_UV;
        out vec4 Frag_Color;
        void main()
        {
            Frag_UV = UV;
            Frag_Color = Color;
            gl_Position = ProjMtx * vec4(Position.xy,0,1);
        }"""

    vertex_shader_glsl_410_core = """
        #version $(ctx.GlslVersion[])
        layout (location = 0) in vec2 Position;
        layout (location = 1) in vec2 UV;
        layout (location = 2) in vec4 Color;
        uniform mat4 ProjMtx;
        out vec2 Frag_UV;
        out vec4 Frag_Color;
        void main()
        {
            Frag_UV = UV;
            Frag_Color = Color;
            gl_Position = ProjMtx * vec4(Position.xy,0,1);
        }"""

    fragment_shader_glsl_130 = """
        #version $(ctx.GlslVersion[])
        uniform sampler2D Texture;
        in vec2 Frag_UV;
        in vec4 Frag_Color;
        out vec4 Out_Color;
        void main()
        {
            Out_Color = Frag_Color * texture(Texture, Frag_UV.st);
        }"""

    fragment_shader_glsl_410_core = """
        #version $(ctx.GlslVersion[])
        in vec2 Frag_UV;
        in vec4 Frag_Color;
        uniform sampler2D Texture;
        layout (location = 0) out vec4 Out_Color;
        void main()
        {
            Out_Color = Frag_Color * texture(Texture, Frag_UV.st);
        }"""

    if ctx.GlslVersion == 410
        vertex_shader = vertex_shader_glsl_410_core
        fragment_shader = fragment_shader_glsl_410_core
    else
        vertex_shader = vertex_shader_glsl_130
        fragment_shader = fragment_shader_glsl_130
    end

    # backup GL state
    last_texture, last_array_buffer, last_vertex_array_object = GLint(0), GLint(0), GLint(0)
    @c glGetIntegerv(GL_TEXTURE_BINDING_2D, &last_texture)
    @c glGetIntegerv(GL_ARRAY_BUFFER_BINDING, &last_array_buffer)
    @c glGetIntegerv(GL_VERTEX_ARRAY_BINDING, &last_vertex_array_object)

    ctx.VertHandle = glCreateShader(GL_VERTEX_SHADER)
    glShaderSource(ctx.VertHandle, 1, Ptr{GLchar}[pointer(vertex_shader)], C_NULL)
    glCompileShader(ctx.VertHandle)
    status = GLint(-1)
    @c glGetShaderiv(ctx.VertHandle, GL_COMPILE_STATUS, &status)
    if status != GL_TRUE
        max_length = GLsizei(0)
        @c glGetShaderiv(ctx.VertHandle, GL_INFO_LOG_LENGTH, &max_length)
        actual_length = GLsizei(0)
        log = Vector{GLchar}(undef, max_length)
        @c glGetShaderInfoLog(ctx.VertHandle, max_length, &actual_length, log)
        @error "[GL]: failed to compile vertex shader: $(ctx.VertHandle): $(String(log))"
    end

    ctx.FragHandle = glCreateShader(GL_FRAGMENT_SHADER)
    glShaderSource(ctx.FragHandle, 1, Ptr{GLchar}[pointer(fragment_shader)], C_NULL)
    glCompileShader(ctx.FragHandle)
    status = GLint(-1)
    @c glGetShaderiv(ctx.FragHandle, GL_COMPILE_STATUS, &status)
    if status != GL_TRUE
        max_length = GLsizei(0)
        @c glGetShaderiv(ctx.FragHandle, GL_INFO_LOG_LENGTH, &max_length)
        actual_length = GLsizei(0)
        log = Vector{GLchar}(undef, max_length)
        @c glGetShaderInfoLog(ctx.FragHandle, max_length, &actual_length, log)
        @error "[GL]: failed to compile fragment shader: $(ctx.FragHandle): $(String(log))"
    end

    ctx.ShaderHandle = glCreateProgram()
    glAttachShader(ctx.ShaderHandle, ctx.VertHandle)
    glAttachShader(ctx.ShaderHandle, ctx.FragHandle)
    glLinkProgram(ctx.ShaderHandle)
    status = GLint(-1)
    @c glGetProgramiv(ctx.ShaderHandle, GL_LINK_STATUS, &status)
    @assert status == GL_TRUE

    ctx.AttribLocationTex = glGetUniformLocation(ctx.ShaderHandle, "Texture")
    ctx.AttribLocationProjMtx = glGetUniformLocation(ctx.ShaderHandle, "ProjMtx")
    ctx.AttribLocationVtxPos = glGetAttribLocation(ctx.ShaderHandle, "Position")
    ctx.AttribLocationVtxUV = glGetAttribLocation(ctx.ShaderHandle, "UV")
    ctx.AttribLocationVtxColor = glGetAttribLocation(ctx.ShaderHandle, "Color")

    # create buffers
    @c glGenBuffers(1, &ctx.VboHandle)
    @c glGenBuffers(1, &ctx.ElementsHandle)

    ImGui_ImplOpenGL3_CreateFontsTexture(ctx)

    # restore modified GL state
    glBindTexture(GL_TEXTURE_2D, last_texture)
    glBindBuffer(GL_ARRAY_BUFFER, last_array_buffer)
    glBindVertexArray(last_vertex_array_object)

    return true;
end

function ImGui_ImplOpenGL3_DestroyDeviceObjects(ctx::Context)
    if ctx.VboHandle != 0
        @c glDeleteBuffers(1, &ctx.VboHandle)
        ctx.VboHandle = 0
    end

    if ctx.ElementsHandle != 0
        @c glDeleteBuffers(1, &ctx.ElementsHandle)
        ctx.ElementsHandle = 0
    end

    if ctx.ShaderHandle != 0 && ctx.VertHandle != 0
        glDetachShader(ctx.ShaderHandle, ctx.VertHandle)
    end

    if ctx.ShaderHandle != 0 && ctx.FragHandle != 0
        glDetachShader(ctx.ShaderHandle, ctx.FragHandle)
    end

    if ctx.VertHandle[] != 0
        glDeleteShader(ctx.VertHandle)
        ctx.VertHandle = 0
    end

    if ctx.FragHandle[] != 0
        glDeleteShader(ctx.FragHandle)
        ctx.FragHandle = 0
    end

    if ctx.ShaderHandle != 0
        glDeleteProgram(ctx.ShaderHandle)
        ctx.ShaderHandle = 0
    end

    ImGui_ImplOpenGL3_DestroyFontsTexture(ctx)

    return true
end


text:
# Randomized SVD

function rsvd(X,    # (m,n) array
              k::Int,           # target rank
              p::Int=1,         # oversampling parameter
              q::Int=5)         # power iterations parameter

    Q, B = rqb(X, k, p, q)
    U, S, V = svd(B)
    U = Q * U

    return U[:, 1:k], S[1:k], V[:, 1:k]
end


text:
@testset "$TEST $G" begin

g3 = PathDiGraph(3, DG)
complete!(g3)
@test ne(g3) == 4
@test nv(g3) == 3
@test has_edge(g3, 3, 2)
@test has_edge(g3, 2, 1)
@test g3 == complete(PathDiGraph(3, DG))

g3 = PathGraph(5, G)
g4 = PathDiGraph(5, DG)
g5 = DG(4)
add_edge!(g5,1,2); add_edge!(g5,2,3); add_edge!(g5,1,3); add_edge!(g5,3,4)

h1 = G(5)
h2 = G(3)
h3 = G()
h4 = DG(7)
h5 = DG()

c3 = complement(g3)
c4 = complement(g4)

@test nv(c3) == 5
@test ne(c3) == 6
@test nv(c4) == 5
@test ne(c4) == 16

g = blockdiag(g3, g3)
@test nv(g) == 10
@test ne(g) == 8

h = PathGraph(2, G)
@test intersect(g3, h) == h

h = PathGraph(4, G)
z = difference(g3, h)
@test nv(z) == 5
@test ne(z) == 1
z = difference(h, g3)
@test nv(z) == 4
@test ne(z) == 0

z = symmetric_difference(h,g3)
@test z == symmetric_difference(g3,h)
@test nv(z) == 5
@test ne(z) == 1

h = G(6)
add_edge!(h, 5, 6)
e = E(5, 6)
z = union(g3, h)
@test has_edge(z, e)
@test z == PathGraph(6, G)

h = DG(6)
add_edge!(h, 5, 6)
e = E(5, 6)
z = union(g4, h)
@test has_edge(z, e)
@test z == PathDiGraph(6, DG)

g10 = CompleteGraph(2, G)
h10 = CompleteGraph(2, G)
z = blockdiag(g10, h10)
@test nv(z) == nv(g10) + nv(h10)
@test ne(z) == ne(g10) + ne(h10)
@test has_edge(z, 1, 2)
@test has_edge(z, 3, 4)
@test !has_edge(z, 1, 3)
@test !has_edge(z, 1, 4)
@test !has_edge(z, 2, 3)
@test !has_edge(z, 2, 4)

g10 = G(2)
h10 = G(2)
z = join(g10, h10)
@test nv(z) == nv(g10) + nv(h10)
@test ne(z) == 4
@test !has_edge(z, 1, 2)
@test !has_edge(z, 3, 4)
@test has_edge(z, 1, 3)
@test has_edge(z, 1, 4)
@test has_edge(z, 2, 3)
@test has_edge(z, 2, 4)

p = PathGraph(10, G)
x = p*ones(10)
@test  x[1]==1.0 && all(x[2:end-1].==2.0) && x[end]==1.0

@test size(p) == (10,10)
@test size(p, 1) == size(p, 2) == 10
@test size(p, 3) == 1

g5 = DG(4)
add_edge!(g5,1,2); add_edge!(g5,2,3); add_edge!(g5,1,3); add_edge!(g5,3,4)
@test g5 * ones(nv(g5)) == [2.0, 1.0, 1.0, 0.0]
@test sum(g5, 1) ==  [0, 1, 2, 1]
@test sum(g5, 2) ==  [2, 1, 1, 0]
@test sum(g5) == 4
@test sum(p,1) == sum(p,2)
@test_throws ErrorException sum(p,3)

@test sparse(p) == adjacency_matrix(p)
@test eltype(p) == vertextype(p)
@test length(p) == 100
@test ndims(p) == 2
@test issymmetric(p)
@test !issymmetric(g5)

g22 = CompleteGraph(2, G)
h = cartesian_product(g22, g22)
@test nv(h) == 4
@test ne(h)== 4

g22 = CompleteGraph(2, G)
h = tensor_product(g22, g22)
@test nv(h) == 4
@test ne(h) == 1

nx = 20; ny = 21
gg = PathGraph(nx, G)
hh = PathGraph(ny, G)
c = cartesian_product(gg, hh)
g = crosspath(PathGraph(nx, G), ny)
@test g == c

## test subgraphs ##

g = graph(:bull, G)
n = 3
h = g[1:n]
@test nv(h) == n
@test ne(h) == 3

h = g[[1,2,4]]
@test nv(h) == n
@test ne(h) == 2

h = g[[1,5]]
@test nv(h) == 2
@test ne(h) == 0
@test typeof(h) == typeof(g)

g = DG(100,200)
h = g[5:26]
@test nv(h) == 22
@test typeof(h) == typeof(g)
@test_throws ErrorException g[[1,1]]

r = 5:26
h2, vm = subgraph(g, r)
@test h2 == h
@test vm == r
@test h2 == g[r]

sg, vm = subgraph(CompleteGraph(10, G), 5:8)
@test nv(sg) == 4
@test ne(sg) == 6

gg = CompleteGraph(10, G)
r = 5:8
edg = (e for e in edges(gg) if (src(e) ∈ r && dst(e) ∈ r))
sg, vm = subgraph(gg, edg)
@test nv(sg) == 4
@test ne(sg) == 6

@test_throws ErrorException subgraph(CompleteGraph(10, G), [5,5,5])

sg2, vm = subgraph(CompleteGraph(10, G), [5,6,7,8])
@test sg2 == sg
@test vm[4] == 8

gg5 = CompleteGraph(10, G)
elist = [E(1,2),E(2,3),E(3,4),E(4,5),E(5,1)]
sg, vm = subgraph(gg5, elist)
@test sg == CycleGraph(5, G)
@test sort(vm) == [1:5;]

g10 = StarGraph(10, G)
@test egonet(g10, 1, 0) == G(1, 0)
@test egonet(g10, 1, 1) == g10

g1 = CompleteGraph(10, G)
g2 = copy(g1)
@test contract!(g1, 1, 2, 2) == contract!(g2, 1:2)
@test g1 == g2 == CompleteGraph(9, G)

g1 = CompleteDiGraph(10, DG)
g2 = copy(g1)
@test contract!(g1, 1, 2, 2) == contract!(g2, 1:2)
@test g1 == g2 == CompleteDiGraph(9, DG)

if G <: ANetwork
    g = readnetwork(:lesmis, G)

    h, vm = subgraph(g, 1:2)
    @test length(vprop(h)) == 0

    h, vm = subnetwork(g,1:2)
    length(vprop(h)) == 2
    @test  vprop(g, "label")[1] == "Myriel"
    @test  vprop(g, "label")[2] == "Napoleon"
    @test length(eprop(h)) == 1
    @test eprop(h, "value")[1,2] == 1

    elist = collect(edges(g))[1:10]
    h, vm = subnetwork(g, elist)
    @test ne(h) == 10
    for e in elist
        @test src(e) ∈ vm && dst(e) ∈ vm
    end
end

end # testset


text:
using StaticArrays
using OffsetArrays
import Base.Iterators: drop, take
import LinearAlgebra: dot

include("space_derivatives/upwind.jl")
include("space_derivatives/weno3.jl")
include("space_derivatives/weno5.jl")
include("space_derivatives/second_order.jl")

function compute_flux!(U :: ScalarQuantity)
    posWENO3!(U.posFlux, U.data)
    negWENO3!(U.negFlux, U.data)
    U.posFlux[0], U.negFlux[end] = U.posFlux[1], U.negFlux[end-1]
end

function compute_fluxP!(pb :: RelaxPb)
    compute_flux!(pb.U)
    @. pb.fluxP = 0.5 * (pb.p(pb.U.posFlux) + pb.p(pb.U.negFlux))
    nothing
end

function compute_fluxV!(pb :: RelaxPb)
    compute_flux!(pb.V)
    @. pb.fluxV = 0.5 * (pb.V.posFlux + pb.V.negFlux)
    nothing
end

function compute_dxP!(pb :: RelaxPb, dx⁻¹)
    compute_flux!(pb.U)
    @. pb.fluxP = 0.5 * ( pb.p(pb.U.posFlux) + pb.p(pb.U.negFlux) )

    for i in 1 : length(pb.dxP)
        pb.dxP[i] = (pb.fluxP[i] - pb.fluxP[i-1]) * dx⁻¹
    end

    nothing
end

function compute_dxV!(pb :: RelaxPb, dx⁻¹)
    compute_flux!(pb.V)
    @. pb.fluxV = 0.5 * ( pb.V.posFlux + pb.V.negFlux )

    for i in 1 : length(pb.dxV)
        pb.dxV[i] = (pb.fluxV[i] - pb.fluxV[i-1]) * dx⁻¹
    end

    nothing
end

max_stars_repo_name: tremelow/StiffKinetic

max_issues_repo_name: tremelow/StiffKinetic

max_forks_repo_name: tremelow/StiffKinetic


text:
# ---
# title: 912. Sort an Array
# id: problem912
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories:
# link: <https://leetcode.com/problems/sort-an-array/description/>
# hidden: true
# ---
#
# Given an array of integers `nums`, sort the array in ascending order.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [5,2,3,1]
#     Output: [1,2,3,5]
#
#
# **Example 2:**
#
#
#
#     Input: nums = [5,1,1,2,0,0]
#     Output: [0,0,1,1,2,5]
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 50000`
#   * `-50000 <= nums[i] <= 50000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end