text:
using HDF5, JLD
using Test

module JLDTemp
using HDF5, JLD
include("JLDTest.jl")  # Private module inside JLDTemp

function create()
    x = JLDTest.Object(convert(Int16, 5))  # int16 makes this work on 0.2
    jldopen("require.jld", "w") do file
        addrequire(file, :JLDTest)
        truncate_module_path(file, JLDTemp)
        write(file, "x", x)
    end
end
end # module

JLDTemp.create()

push!(LOAD_PATH, splitdir(@__FILE__)[1])
x = jldopen("require.jld") do file
    read(file, "x")
end
@test typeof(x) == JLDTest.Object
@test x.data == 5
pop!(LOAD_PATH)
rm("require.jld")


text:
using Test
using Petri
using LabelledArrays

@testset "Equality" begin
    sir_1 = Petri.Model([1,2,3],[([1,1], [2]),
                                 ([1],   [1])])
    sir_2 = Petri.Model(1:3,[([1,1], [2]),
                             ([1],   [1])])
    sir_3 = Petri.Model(3,[([1,1], [2]),
                           ([1],   [1])])
    @test typeof(Graph(sir_1)) == Graph
    @test sir_1 == sir_2
    @test sir_1 == sir_3

    @test EmptyPetri(5) == Petri.Model(1:5, [])

    @test EmptyPetri(5) == EmptyPetri(collect(1:5))

    sir = Petri.Model([:S,:I,:R],[(LVector(S=1,I=1), LVector(I=2)),
                                  (LVector(I=1),     LVector(R=1))])
    @test sir == sir
    seir = Petri.Model([:S,:E,:I,:R],[(LVector(S=1,I=1), LVector(I=1,E=1)),
                                      (LVector(E=1),     LVector(I=1)),
                                      (LVector(I=1),     LVector(R=1))])
    x = Petri.Model([:S,:E,:I,:R],[(LVector(S=1,I=1), LVector(I=2)),
                                   (LVector(I=1),     LVector(R=1))])
    y = Petri.Model([:S,:E,:I,:R],[(LVector(S=1,I=1), LVector(I=1,E=1)),
                                   (LVector(E=1),     LVector(I=1)),
                                   (LVector(I=1),     LVector(R=1)),
                                   (LVector(R=1),     LVector(S=1))])
    @test sir != seir
    @test seir != x
    @test seir != y

end

include("solvers.jl")

ext: jl


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

mutable struct IoK8sApiAuthenticationV1beta1UserInfo <: SwaggerModel
    extra::Any # spec type: Union{ Nothing, Dict{String, Vector{String}} } # spec name: extra
    groups::Any # spec type: Union{ Nothing, Vector{String} } # spec name: groups
    uid::Any # spec type: Union{ Nothing, String } # spec name: uid
    username::Any # spec type: Union{ Nothing, String } # spec name: username

    function IoK8sApiAuthenticationV1beta1UserInfo(;extra=nothing, groups=nothing, uid=nothing, username=nothing)
        o = new()
        validate_property(IoK8sApiAuthenticationV1beta1UserInfo, Symbol("extra"), extra)
        setfield!(o, Symbol("extra"), extra)
        validate_property(IoK8sApiAuthenticationV1beta1UserInfo, Symbol("groups"), groups)
        setfield!(o, Symbol("groups"), groups)
        validate_property(IoK8sApiAuthenticationV1beta1UserInfo, Symbol("uid"), uid)
        setfield!(o, Symbol("uid"), uid)
        validate_property(IoK8sApiAuthenticationV1beta1UserInfo, Symbol("username"), username)
        setfield!(o, Symbol("username"), username)
        o
    end
end # type IoK8sApiAuthenticationV1beta1UserInfo

const _property_map_IoK8sApiAuthenticationV1beta1UserInfo = Dict{Symbol,Symbol}(Symbol("extra")=>Symbol("extra"), Symbol("groups")=>Symbol("groups"), Symbol("uid")=>Symbol("uid"), Symbol("username")=>Symbol("username"))
const _property_types_IoK8sApiAuthenticationV1beta1UserInfo = Dict{Symbol,String}(Symbol("extra")=>"Dict{String, Vector{String}}", Symbol("groups")=>"Vector{String}", Symbol("uid")=>"String", Symbol("username")=>"String")
Base.propertynames(::Type{ IoK8sApiAuthenticationV1beta1UserInfo }) = collect(keys(_property_map_IoK8sApiAuthenticationV1beta1UserInfo))
Swagger.property_type(::Type{ IoK8sApiAuthenticationV1beta1UserInfo }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAuthenticationV1beta1UserInfo[name]))}
Swagger.field_name(::Type{ IoK8sApiAuthenticationV1beta1UserInfo }, property_name::Symbol) =  _property_map_IoK8sApiAuthenticationV1beta1UserInfo[property_name]

function check_required(o::IoK8sApiAuthenticationV1beta1UserInfo)
    true
end

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


text:
using Test
using Random
using Inflate
using CodecZlib: GzipCompressorStream, ZlibCompressorStream,
                 DeflateCompressorStream, CodecZlib

empty_string = ""
short_string = "This is a short string."
medium_string = read(pathof(Inflate), String)
long_string = join(fill(medium_string, 1000), short_string)

@testset "Text strings" begin
    for s in [empty_string, short_string, medium_string, long_string]
        @test String(inflate(read(DeflateCompressorStream(IOBuffer(s))))) == s
        @test String(inflate_zlib(read(ZlibCompressorStream(IOBuffer(s))))) == s
        @test String(inflate_gzip(read(GzipCompressorStream(IOBuffer(s))))) == s
        @test read(InflateStream(DeflateCompressorStream(IOBuffer(s))), String) == s
        @test read(InflateZlibStream(ZlibCompressorStream(IOBuffer(s))), String) == s
        @test read(InflateGzipStream(GzipCompressorStream(IOBuffer(s))), String) == s
    end
end

@testset "Incompressible data" begin
    Random.seed!(1)
    for n in [0, 1, 10, 100, 1000, 10000, 100000, 1000000]
        data = rand(UInt8, n)
        @test inflate(read(DeflateCompressorStream(IOBuffer(data)))) == data
        @test inflate_zlib(read(ZlibCompressorStream(IOBuffer(data)))) == data
        @test inflate_gzip(read(GzipCompressorStream(IOBuffer(data)))) == data
        @test read(InflateStream(DeflateCompressorStream(IOBuffer(data)))) == data
        @test read(InflateZlibStream(ZlibCompressorStream(IOBuffer(data)))) == data
        @test read(InflateGzipStream(GzipCompressorStream(IOBuffer(data)))) == data
    end
end

@testset "Huffman compressible data" begin
    Random.seed!(1)
    for n in [0, 1, 10, 100, 1000, 10000, 100000, 1000000]
        data = rand(UInt8, n) .& 0x0f
        @test inflate(read(DeflateCompressorStream(IOBuffer(data)))) == data
        @test inflate_zlib(read(ZlibCompressorStream(IOBuffer(data)))) == data
        @test inflate_gzip(read(GzipCompressorStream(IOBuffer(data)))) == data
        @test read(InflateStream(DeflateCompressorStream(IOBuffer(data)))) == data
        @test read(InflateZlibStream(ZlibCompressorStream(IOBuffer(data)))) == data
        @test read(InflateGzipStream(GzipCompressorStream(IOBuffer(data)))) == data
    end
end

# Test gzip headers, including header CRC.
include("gzip_with_header.jl")
os = 255   # Unknown
mtime = 1537006040
fextra = UInt8[0x00, 0xbe, 0xef, 0x00]
fname = "foo.txt"
fcomment = "example string"
gz = gzip_with_header("foo", mtime, os, fextra, fname, fcomment, true)

@testset "Gzip headers" begin
    headers = Dict{String, Any}()
    @test String(inflate_gzip(gz, headers = headers)) == "foo"
    @test headers["mtime"] == mtime
    @test headers["os"] == os
    @test headers["fextra"] == fextra
    @test headers["fname"] == fname
    @test headers["fcomment"] == fcomment

    headers = Dict{String, Any}()
    @test read(InflateGzipStream(IOBuffer(gz), headers = headers), String) == "foo"
    @test headers["mtime"] == mtime
    @test headers["os"] == os
    @test headers["fextra"] == fextra
    @test headers["fname"] == fname
    @test headers["fcomment"] == fcomment
end

# Test readline interface and get coverage of read(stream, UInt8).
@testset "readline" begin
    s = "one line of text\n"
    @test readline(InflateGzipStream(GzipCompressorStream(IOBuffer(s))),
                   keep = true) == s
end

# Test failure cases, mostly corrupt data.
include("provoke_errors.jl")


text:
## Load the MNIST Graph and draw it
##
using DelimitedFiles
using SparseArrays
using LinearAlgebra
ei,ej = copy.(eachcol(Int.(readdlm("mnist-train-4.edges"))))
n = maximum(ei)
A = sparse(ei,ej,1,n,n)
@assert issymmetric(A)
##
using PyCall
using Conda
pyimport_conda("igraph","python-igraph","conda-forge")
##
using PyCall
igraph = pyimport("igraph")

function igraph_layout(A::SparseMatrixCSC{T}, layoutname::AbstractString="drl") where T
    ei,ej,ew = findnz(A)
    edgelist = [(ei[i]-1,ej[i]-1) for i = 1:length(ei)]
    nverts = size(A)
    G = igraph.Graph(nverts, edges=edgelist, directed=false)
    xy = G.layout(layoutname)
    xy = [Float64(xy[i][j]) for i in 1:length(xy),  j in 1:length(xy[1])]
end

xy = igraph_layout(A, "drl")
##
writedlm("mnist-train-4.xy", xy)
##
xy = readdlm("mnist-train-4.xy")
##
using GraphRecipes
using LinearAlgebra
using Plots
graphplot(ei, ej, x =xy[:,1], y=xy[:,2],
  markercolor=:black, markerstrokecolor=:white,
  size=(1200,1200),dpi=400,
  markersize=2, linecolor=1, linealpha=0.2, linewidth=0.5,
  markeralpha=0.3,
  #colorbar=true,
  marker_z=train_y,
  axis_buffer=0.02, background=nothing)
for i=1:100
  annotate!(xy[i,1],xy[i,2],"$(train_y[i])",color=:white)
end
plot!()
##
savefig("mnist-train-4.png")


text:
export style

function style()
  """
  html, body {
    height: 100%;
  }

  html {
    overflow-y: scroll;
  }

  body {
    padding: 0;
    margin: 0;
    background: -moz-linear-gradient(top, rgba(238,238,238,1) 0%, rgba(238,238,238,0) 100%);
    background: -webkit-gradient(linear, left top, left bottom, color-stop(0%,rgba(238,238,238,1)), color-stop(100%,rgba(238,238,238,0)));
    background: -webkit-linear-gradient(top, rgba(238,238,238,1) 0%,rgba(238,238,238,0) 100%);
    background: -o-linear-gradient(top, rgba(238,238,238,1) 0%,rgba(238,238,238,0) 100%);
    background: -ms-linear-gradient(top, rgba(238,238,238,1) 0%,rgba(238,238,238,0) 100%);
    background: linear-gradient(to bottom, rgba(238,238,238,1) 0%,rgba(238,238,238,0) 100%);
    filter: progid:DXImageTransform.Microsoft.gradient( startColorstr='#eeeeee', endColorstr='#00eeeeee',GradientType=0 );
  }

  #splash_header {
    font: 300 18px 'light','Helvetica Neue',Arial,Helvetica,sans-serif;
  }
  #room_entry {
    text-align: center;
    width: 100%;
    display: block;
  }

  #username_header {
    font: 300 18px 'light','Helvetica Neue',Arial,Helvetica,sans-serif;
  }

  #username_warning_header {
    font: 300 18px 'light','Helvetica Neue',Arial,Helvetica,sans-serif;
  }

  #warning {
     margin-top: 5px;
  }

  #warning td {
     border: 1px solid black;
     vertical-align: top;
  }

  #colors {
     border-spacing: 0px;
  }

  #good {
     background-color: rgb(190, 231, 190);
  }

  #yellow {
     background-color: rgb(247, 247, 195);
  }

  #bad {
     background-color: rgb(240, 203, 203);
  }

  #chatinput {
    background: none repeat scroll 0 0 #F9F9F9;
    border: 1px solid #CCCCCC;
    border-radius: 0px 0px 5px 5px;
    font: 300 18px/40px 'light','Helvetica Neue',Arial,Helvetica,sans-serif;
    height: 40px;
    letter-spacing: 1px;
    margin-bottom: 0px;
    padding: 5px 10px;
    width: 100%;
  }

  #chatbox {
    height: 95%;
  }

  #chatbox-panel {
    height: calc(100% - 45px);
  }

  #messages {
    overflow: scroll;
    background-color: rgb(249,249,249,0.6);
    border: 1px solid #CCCCCC;
    border-radius: 5px 5px 0px 0px;
    font: 300 18px 'light','Helvetica Neue',Arial,Helvetica,sans-serif;
    height: calc(100% - 45px);
    min-height:100px;
    margin-bottom: 0px;
    padding: 5px 10px;
    width: 100%;
  }

  .small {
      font-size: 12px;
      color: rgb(131, 131, 131);
  }

  .flip {
  -webkit-transform: rotateY(180deg);
  }

  .buttonBox {
    width: 100%;
    height: 100%;
    display: block;
    position: absolute;
    top: 2px;
    left: 2px;
    z-index: 1;
  }
  .button {
    -moz-box-shadow:inset 0px 1px 0px 0px #ffffff;
    -webkit-box-shadow:inset 0px 1px 0px 0px #ffffff;
    box-shadow:inset 0px 1px 0px 0px #ffffff;
    background:-webkit-gradient( linear, left top, left bottom, color-stop(0.05, #ededed), color-stop(1, #dfdfdf) );
    background:-moz-linear-gradient( center top, #ededed 5%, #dfdfdf 100% );
    filter:progid:DXImageTransform.Microsoft.gradient(startColorstr='#ededed', endColorstr='#dfdfdf');
    background-color:#ededed;
    -moz-border-radius:6px;
    -webkit-border-radius:6px;
    border-radius:6px;
    border:1px solid #dcdcdc;
    display:inline-block;
    color:rgb(78, 78, 78);
    font-family:arial;
    font-size:15px;
    font-weight:bold;
    padding:6px 24px;
    text-decoration:none;
    text-shadow:1px 1px 0px #ffffff;
  }

  .button:hover {
    background:-webkit-gradient( linear, left top, left bottom, color-stop(0.05, #dfdfdf), color-stop(1, #ededed) );
    background:-moz-linear-gradient( center top, #dfdfdf 5%, #ededed 100% );
    filter:progid:DXImageTransform.Microsoft.gradient(startColorstr='#dfdfdf', endColorstr='#ededed');
    background-color:#dfdfdf;
  }
  .button:active {
    position:relative;
    top:1px;
  }

  /* offcanvas userlist on smaller screens, ie. mobile */

  @media screen and (max-width: 768px) {
    .row-offcanvas {
      position: relative;
      -webkit-transition: all 0.25s ease-out;
      -moz-transition: all 0.25s ease-out;
      transition: all 0.25s ease-out;
    }

    .row-offcanvas-right
    .sidebar-offcanvas {
      right: -50%; /* 6 columns */
    }

    .row-offcanvas-left
    .sidebar-offcanvas {
      left: -50%; /* 6 columns */
    }

    .row-offcanvas-right.active {
      right: 50%; /* 6 columns */
    }

    .row-offcanvas-left.active {
      left: 50%; /* 6 columns */
    }

    .sidebar-offcanvas {
      position: absolute;
      top: 0;
      width: 50%; /* 6 columns */
    }

    .btn_off {
      float: left;
      position: absolute;
      top: -20px;
    }
  }

  /* bootstrap overrides */

  body {
  	padding-top:70px; /* header shift */
  }

  .container {
  	max-width: 2170px;
  }

  .panel {
  	background-color:rgba(255,255,255,0.6);
  }

  .panel-default > .panel-heading {
  	background-color:rgba(245,245,245,0.6);
  }

  input[type="file"] {
  	display:inline;
  }

  .navbar-inverse .navbar-brand {
    color: #fff;
  }

  @media screen and (min-width: 768px) {
  	.modal-dialog {
  		width: 680px;
  	}
  }

  .frontpage_header {
  	text-align:center;
  }

  .support {
    color:rgb(179, 47, 47);
    display:none;
  }

  .connect-btn {
    margin-top:22px;
  }

  .username-panel {
    width: 360px;
    margin-left: auto;
    margin-right: auto;
  }

  .form-signin {
    max-width: 330px;
    padding: 15px;
    margin: 0 auto;
  }
  .form-signin .form-signin-heading,
  .form-signin .checkbox {
    margin-bottom: 10px;
  }
  .form-signin .checkbox {
    font-weight: normal;
  }
  .form-signin .form-control {
    position: relative;
    font-size: 16px;
    height: auto;
    padding: 10px;
    -webkit-box-sizing: border-box;
       -moz-box-sizing: border-box;
            box-sizing: border-box;
  }
  .form-signin .form-control:focus {
    z-index: 2;
  }
  .form-signin input[type="text"] {
    margin-bottom: -1px;
    border-bottom-left-radius: 0;
    border-bottom-right-radius: 0;
  }
  .form-signin input[type="password"] {
    margin-bottom: 10px;
    border-top-left-radius: 0;
    border-top-right-radius: 0;
  }
  """
end


text:
const _attribute_defaults =  Dict(:Series => _series_defaults,
                        :Subplot => _subplot_defaults,
                        :Plot => _plot_defaults,
                        :Axis => _axis_defaults)

attrtypes() = join(keys(_attribute_defaults), ", ")
attributes(attrtype::Symbol) = sort(collect(keys(_attribute_defaults[attrtype])))

function lookup_aliases(attrtype, attribute)
    attribute = Symbol(attribute)
    attribute = in(attribute, keys(_keyAliases)) ? _keyAliases[attribute] : attribute
    in(attribute, keys(_attribute_defaults[attrtype])) && return attribute
    error("There is no attribute named $attribute in $attrtype")
end

"""
    plotattr([attr])

Look up the properties of a Plots attribute, or specify an attribute type. Call `plotattr()` for options.
The information is the same as that given on https://juliaplots.github.io/attributes/.
"""
function plotattr()
    println("Specify an attribute type to get a list of supported attributes. Options are $(attrtypes())")
end

function plotattr(attrtype::Symbol)
    in(attrtype, keys(_attribute_defaults)) || error("Viable options are $(attrtypes())")
    println("Defined $attrtype attributes are:\n$(join(attributes(attrtype), ", "))")
end

function plotattr(attribute::AbstractString)
    attribute = Symbol(attribute)
    attribute = in(attribute, keys(_keyAliases)) ? _keyAliases[attribute] : attribute
    for (k, v) in _attribute_defaults
        if in(attribute, keys(v))
            return plotattr(k, "$attribute")
        end
    end
    error("There is no attribute named $attribute")
end

printnothing(x) = x
printnothing(x::Nothing) = "nothing"

function plotattr(attrtype::Symbol, attribute::AbstractString)
    in(attrtype, keys(_attribute_defaults)) || ArgumentError("`attrtype` must match one of $(attrtypes())")

    attribute = Symbol(lookup_aliases(attrtype, attribute))

    desc = get(_arg_desc, attribute, "")
    first_period_idx = findfirst(isequal('.'), desc)
    typedesc = desc[1:first_period_idx-1]
    desc = strip(desc[first_period_idx+1:end])
    als = keys(filter(x->x[2]==attribute, _keyAliases)) |> collect |> sort
    als = join(map(string,als), ", ")
    def = _attribute_defaults[attrtype][attribute]

    # Looks up the different elements and plots them
    println("$(printnothing(attribute)) ", typedesc == "" ? "" : "{$(printnothing(typedesc))}", "\n",
        als == "" ? "" : "$(printnothing(als))\n",
        "\n$(printnothing(desc))\n",
        "$(printnothing(attrtype)) attribute, ", def == "" ? "" : " default: $(printnothing(def))")
end


text:
################################################################################
#       0. CONFIGURATION
################################################################################

# Setup initial configuration
begin
    num_workers = 4
    include("setup.jl")
end
# NOTE: first time running can take some minutes, packages installing

################################################################################
#       1. HYPERPARAMETER SEARCH
################################################################################

# Grid of values
@everywhere begin
    # Define potential hyperparameter values
    LearningRates = [[0.001,0.0025],[0.001,0.005],[0.001,0.01],[0.01,0.025],[0.01,0.05],[0.01,0.1]]
    Minibatches = [1,5,10]
    # Parameters for run
    max_epochs = 1000
    early_stoping = 100
    report = 50
    # Iterator over hyperparams, params and repetitions
    inx = collect(product(enumerate(LearningRates),enumerate(Minibatches)))[:]
    # Select "Drosophila_Gut" and "Soil_Vitro" as examples
    pars = real_data[1:2]
    # NOTE: variable `real_data` and values imported from module cNODE
    "search hyperparameters..." |> println
end

for (i,DATA) in enumerate(pars)
    # Load percentage of dataset
    path = "./data/real/$DATA/P.csv"
    Z,P = import_data(path,2/3)
    N,M = size(Z)
    # Explore hyperparameters in small dataset
    "training $DATA..." |> println
    for it in inx
        ((j,lr),(k,mb)) = it
        mb = Minibatches[k]
        "LR: $lr MB: $mb" |> println
        Qtst = SharedArray{Float64}(M,N)
        Ptst = SharedArray{Float64}(M,N)
        LossTrain = SharedArray{Float64}(M)
        LossTest = SharedArray{Float64}(M)
        # Use Leave-one-out cross validation
        LeaveOneOut = kfolds((Z,P); k = M) |> enumerate |> collect
        @sync @distributed for l in LeaveOneOut
            (l,((ztrn,ptrn),(ztst,ptst))) = l
            "training $l..."|>println
            # Get cNODE
            cnode = getModel(N)
            # Train cNODE
            W, loss_train, loss_val, loss_test = train_reptile(
                                                    cnode, max_epochs,
                                                    mb, lr,
                                                    ztrn, ptrn, ztst, ptst, ztst, ptst,
                                                    report, early_stoping
                                                )
            # Save
            Ptst[l,:] = ptst
            Qtst[l,:] = predict(cnode,ztst)
            LossTrain[l] = loss_train[end]
            LossTest[l] = loss_test[end]
            # Report
            println(l,'\t',loss_train[end],'\t',loss_test[end])
            println('#'^30)
        end
        # Save results
        results = "./results/real/$DATA/hyperparameters/"
        !ispath(results) && mkpath(results)
        writedlm(results*"real_sample_$(j)$(k).csv", Ptst, ',')
        writedlm(results*"pred_sample_$(j)$(k).csv", Qtst, ',')
        writedlm(results*"test_loss_$(j)$(k).csv",   LossTest, ',')
        writedlm(results*"train_loss_$(j)$(k).csv",  LossTrain, ',')
    end
end

################################################################################
#       2. Experimental Validation
################################################################################

for (i,DATA) in enumerate(pars)
    # Import full dataset
    path = "./data/real/$DATA/P.csv"
    Z,P = import_data(path)
    N,M = size(Z)
    # Select hyperparameters
    results = "./results/real/$DATA/hyperparameters/"
    _mean = [ mean(readdlm(results*"test_loss_$i$j.csv",',',Float64,'\n')) for i in 1:6, j in 1:3] |> argmin
    mb = Minibatches[_mean[2]]
    lr = LearningRates[_mean[1]]
    # Allocate variables
    Qtst = SharedArray{Float64}(M,N)
    Ptst = SharedArray{Float64}(M,N)
    LossTrain = SharedArray{Float64}(M)
    LossTest = SharedArray{Float64}(M)
    # Run validation
    results = "./results/real/$DATA/"
    LeaveOneOut = kfolds((Z,P); k = M) |> enumerate |> collect
    "training $DATA..." |> println
    @sync @distributed for l in LeaveOneOut
        (l,((ztrn,ptrn),(ztst,ptst))) = l
        "training $l..."|>println
        # Get cNODE
        cnode = getModel(N)
        # Train cNODE
        W, loss_train, loss_val, loss_test = train_reptile(
                                                cnode, max_epochs,
                                                mb, lr,
                                                ztrn, ptrn, ztst, ptst, ztst, ptst,
                                                report, early_stoping
                                            )
        # Save values
        Ptst[l,:] = ptst
        Qtst[l,:] = predict(cnode,ztst)
        LossTrain[l] = loss_train[end]
        LossTest[l] = loss_test[end]
        # Save realization
        !ispath(results*"loss_epochs/") && mkpath(results*"loss_epochs/")
        writedlm(results*"loss_epochs/train$l.csv",loss_train, ',')
        writedlm(results*"loss_epochs/test$l.csv",loss_test, ',')
        writedlm(results*"loss_epochs/val$l.csv",loss_val, ',')
        # Report
        println(i,'\t',loss_train[end],'\t',loss_test[end])
    end
    # Write full results
    writedlm(results*"real_sample.csv",Ptst, ',')
    writedlm(results*"pred_sample.csv",Qtst, ',')
    writedlm(results*"test_loss.csv",LossTest, ',')
    writedlm(results*"train_loss.csv",LossTrain, ',')
end


text:
using DictIO
using Base.Test

# ASCIIString,Int64 dict
d = Dict{ASCIIString,Int64}()
d["a"] = 2
d["b"] = 6
d["g"] = 1
writedict("tmp.txt", d)
d2 = readdict("tmp.txt")
@test d2["a"] == 2
@test d2["b"] == 6
@test d2["g"] == 1
@test typeof(d2["g"]) == Int64
@test typeof(first(keys(d2))) == ASCIIString

# Float64,Any dict reversed
d = Dict{Float64,Any}()
d[1.0] = 8
d[100.0] = "6"
open(f->writedict(f, d), "tmp.txt", "w")
d2 = readdict("tmp.txt", reversed=true)
@test d2["8"] == 1.0
@test d2["6"] == 100.0
@test typeof(d2["8"]) == Float64


text:
x = Variable(4, 3)
y = Variable()

x + y == 3
3 == y
rand(4, 3) == x - 4*y
sum(y, 1) == rand(1, 3)

# lists of constraints
list = []
list += x + y == ones(4, 3)

list2 = EqConstraint[]
list2 += y == 3
list2 += list

info("All equality constraints tests passed")


text:
abstract type RTMap{Name,MT,K,V,ME,F} end
abstract type AbstractHashMap{Name,MT,K,V,ME,F} <: RTMap{Name,MT,K,V,ME,F} end
abstract type AbstractArrayMap{Name,MT,K,V,ME,F} <: RTMap{Name,MT,K,V,ME,F} end

struct HashMap{Name,MT,K,V,ME,F} <: AbstractHashMap{Name,MT,K,V,ME,F} end
maptype_to_jltype(::Val{API.BPF_MAP_TYPE_HASH}) = HashMap
maptype_to_jltype(::Val{API.BPF_MAP_TYPE_STACK_TRACE}) = HashMap
struct ArrayMap{Name,MT,K,V,ME,F} <: AbstractArrayMap{Name,MT,K,V,ME,F} end
maptype_to_jltype(::Val{API.BPF_MAP_TYPE_ARRAY}) = ArrayMap

function RTMap(; name, maptype, keytype, valuetype, maxentries=1, flags=0)
    jltype = maptype_to_jltype(Val(maptype))
    jltype{Symbol(name), maptype, keytype, valuetype, maxentries, flags}()
end

function _genmap!(mod::LLVM.Module, ::Type{<:RTMap{Name,MT,K,V,ME,F}}, ctx) where {Name,MT,K,V,ME,F}
    T_i32 = LLVM.Int32Type(ctx)
    T_map = LLVM.StructType([T_i32, T_i32, T_i32, T_i32, T_i32]; ctx)
    name = string(Name)
    gv = GlobalVariable(mod, T_map, name)
    section!(gv, "maps")
    alignment!(gv, 4)
    vec = Any[Int32(MT),Int32(sizeof(K)),Int32(sizeof(V)),Int32(ME),Int32(F)]
    init = ConstantStruct([ConstantInt(v; ctx) for v in vec]; ctx)
    initializer!(gv, init)
    linkage!(gv, LLVM.API.LLVMLinkOnceODRLinkage)
    return gv
end

# From AMDGPU.jl/src/device/gcn/memory_static.jl
@inline function _memset!(builder, ctx, mod, dest, value, len, volatile)
    T_nothing = LLVM.VoidType(ctx)
    T_dest = llvmtype(dest)
    T_int8 = convert(LLVMType, UInt8; ctx)
    T_int64 = convert(LLVMType, UInt64; ctx)
    T_int1 = LLVM.Int1Type(ctx)

    T_intr = LLVM.FunctionType(T_nothing, [T_dest, T_int8, T_int64, T_int1])
    intr = LLVM.Function(mod, "llvm.memset.p$(Int(addrspace(T_dest)))i8.i64", T_intr)
    call!(builder, intr, [dest, value, len, volatile])
end
@inline @generated function memset!(dest_ptr::LLVMPtr{UInt8,DestAS}, value::UInt8, len::LT) where {DestAS,LT<:Union{Int64,UInt64}}
    Context() do ctx
        T_nothing = LLVM.VoidType(ctx)
        T_pint8_dest = convert(LLVMType, dest_ptr; ctx)
        T_int8 = convert(LLVMType, UInt8; ctx)
        T_int64 = convert(LLVMType, UInt64; ctx)
        T_int1 = LLVM.Int1Type(ctx)

        llvm_f, _ = create_function(T_nothing, [T_pint8_dest, T_int8, T_int64])
        mod = LLVM.parent(llvm_f)
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            _memset!(builder, ctx, mod, parameters(llvm_f)[1], parameters(llvm_f)[2], parameters(llvm_f)[3], ConstantInt(T_int1, 0))
            ret!(builder)
        end
        call_function(llvm_f, Nothing, Tuple{LLVMPtr{UInt8,DestAS},UInt8,LT}, :dest_ptr, :value, :len)
    end
end
@inline memset!(dest_ptr::LLVMPtr{T,DestAS}, value::UInt8, len::Integer) where {T,DestAS} =
    memset!(reinterpret(LLVMPtr{UInt8,DestAS}, dest_ptr), value, UInt64(len))
@inline function ZeroInitRef(T, val)
    ref = Ref{T}()
    ref_llptr = reinterpret(LLVMPtr{T,0}, Base.unsafe_convert(Ptr{T}, ref))
    memset!(ref_llptr, UInt8(0), sizeof(T))
    ref[] = val
    ref
end


text:
### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# ╔═╡ 938f32d0-4dc2-11eb-29e1-19b4e28bf3c8
using BenchmarkTools, Random, DataFrames, Plots

# ╔═╡ 12341f32-4dbf-11eb-1753-6fac96fd666a
md"""
# Ancient Algorithms

The concept of algorithm has existed since ancient ages. Babylonian (c. 2500 BC) and Egyptian (c. 1550 BC) used arithmetic algorithms for example division. Greek mathematicians used the sieve of Eratosthenes for finding prime numbers, and the Euclidean algorithm for finding the greatest common divisor of two numbers.

The word algorithm itself was derived from the name of the 9th-century Persian mathematician Muḥammad ibn Mūsā al-Khwārizmī (c. 780–850) , whose name (Khwarazm) was Latinized as Algoritmi.
"""

# ╔═╡ 10ef1d92-4dc0-11eb-3524-bd5cb089a31b
md"""
## The Peasant Multiplication

Origin: Russian and Egyptian

Idea: To multiply two positive integers a and b,
- we can instead compute (a/2)*(2b) if a is even and
- ((a − 1)/2) * (2b) + b if a is odd.

**Algorithm Analysis:**
- a is reduced by at least a factor of 2, at the cost of one halving, one doubling, and possibly one addition.
- If we assume that halving, doubling, and addition are all constant-time operations, then the total computational complexity of the algorithm is proportional to the number of binary digits of a.
- Converges to correct answer
"""

# ╔═╡ 09744cea-4dc0-11eb-11b0-c7c6b42683bc
html"""
<iframe width="560" height="315" src="https://www.youtube.com/embed/HJ_PP5rqLg0" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>
"""

# ╔═╡ aae2f0aa-4dbf-11eb-3f50-39e262cc1845
function peasantprod(a::T, b::T) where T <: Integer
    c = zero(T)
    while a > one(T)
        if isodd(a)
            c = c + b
        end
        a = a >> 1 # divide a by 2
        b = b << 1 # multiply b by 2
    end
    return c + b
end

# ╔═╡ 81d2d826-4dc0-11eb-2efb-9f5ef5161edd
peasantprod(10, 33)

# ╔═╡ 87f00406-4dc0-11eb-29b9-6b5c19d03b53
md"""
## Babylonian Method for finding square root

Still used by many calculators.

Idea:

$$x \cdot \frac{a}{x} = a = \sqrt{a} \sqrt{a}$$

- If $x < \sqrt{a}$, then $\frac{a}{x} > \sqrt{a}$
- If $x > \sqrt{a}$, then $\frac{a}{x} < \sqrt{a}$

Take averge of the two

$X_{n+1} = \frac{X_n + \frac{a}{X_n}}{2} = \frac{X_n^2 + a}{2 X_n}$

"""

# ╔═╡ a78d5570-4dc0-11eb-37a8-adc45a4964ee
function babylonian(a::T, tol::T) where T <: Real
    x = one(T)  # start x at 1
    while abs(x^2 - a) > tol  # convergence test
        x = (x + c / x) / 2
    end
    return x
end

# ╔═╡ c0c78d52-4dc1-11eb-325d-81452aa6d218

# ╔═╡ 0dd9ea7c-4dc2-11eb-1b99-f7b120c21a70
md"""
## [Sieve of Eratosthenes](https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes) for finding primes

Iteratively mark of all composites of the identified prime.

![Sieve of Eratosthenes](https://upload.wikimedia.org/wikipedia/commons/b/b9/Sieve_of_Eratosthenes_animation.gif)

"""

# ╔═╡ 371f8004-4dc2-11eb-340e-c367fc0d718d
begin
	function es(n::Int)
	    isprime = trues(n)  # n-element vector of true-s
	    isprime[1] = false  # 1 is not a prime
	    for i in 2:isqrt(n)  # loop integers less or equal than sqrt(n)
	        if isprime[i]  # conditional evaluation
	            for j in i^2:i:n  # sequence with step i
	                isprime[j] = false
	                end
	            end
	        end
	    return filter(x -> isprime[x], 1:n)# filter using an anonymous function
	    end

	println(es(100))# print all primes less or equal than 100
	@time length(es(10^6))# check function execution time and memory usage
end

# ╔═╡ 5a9e6414-4dc2-11eb-2ed3-8f6abe2b7847
md"""
## Greatest Common Divisor (GCD)

For integers a and b (likely to be very large positive integers), GCD(a,b) is the largest integer d that divides both a and b. GCD is very important in number theory such as factorization, prime numbers and hence cryptography.
"""

# ╔═╡ 61a2ce56-4dc2-11eb-3b58-df1a7f2bf4d6
begin
	"""GCD using the definition of GCD"""
	function GCDnaive(a,b)
	    largest = 0
	    for d in 1:min(a,b)   # runtime ~min(a,b)
	        if a%d==0 && b%d==0
	            largest = d
	        end
	    end
	    return largest
	end

	@assert GCDnaive(8,12)==4
	@benchmark GCDnaive(3918848,1653264)
end

# ╔═╡ c3fdeae2-4dc2-11eb-0d57-598c117205c7
md"""
### Euclidean GCD

**Lemma:**
> Let $a'$ be the remainder when a is divided by b, then $$gcd(a,b) = gcd(a',b)= gcd(b,a')$$

**Proof**
- when a is divided by b, $a = b*q + a'$
- now d divides a and b if and only if it divides $a'$ and $b$

Each step reduces the size of the numbers by about a factor of 2. Takes about log(ab) steps
"""

# ╔═╡ e0401ae0-4dc2-11eb-2ba6-9d18e1915fa1
begin
	function GCDEuclid(a,b)
	    if b == 0
	        return a
	    else
	        return GCDEuclid(b,a%b)
	    end
	end

	# onliner
	GCD(a,b) = b==0 ? a : GCD(b , a%b)

	@benchmark GCDEuclid(3918848,1653264)
end

# ╔═╡ f40e1e6e-4dc2-11eb-2c04-d7020e727f7f
md"""
## Fibonacci Numbers

Developed by to study Rabbit Population after multiple genertions

$$F_n = \begin{cases} 0, \qquad n=0, \\ 1, \qquad n=1, \\ F_{n+1} + F_{n+2}, \quad n>1 \end{cases}$$

This has a rapid growth. $F_n \geq 2^{n/2} \text{ for } n \geq 6$

"""

# ╔═╡ 02317fca-4dc3-11eb-34c4-a5e009b09775
begin
	"""Fibbonaci Implementation directly from definition"""
	function fibbo_naive(n::Integer)
	    n >= 0 || throw(DomainError(n, "Domain: n ≥ 0"))
	    if n <= 1
	        return n
	    else
	        return fibbo_naive(n-1) + fibbo_naive(n-2)
	    end

	end

	# Do not try huge numbers here
	@assert fibbo_naive(1) == 1
	@assert fibbo_naive(7) == 13
end

# ╔═╡ 1b865a74-4dc3-11eb-25ac-d56e25e7b4b6
md"""
If T(n) is the Lines of Code executed by the function on input n. If you step through the debugger you can calculate the number of lines  for the implementation.

$$T(n) = \begin{cases} 2, \qquad n \leq 1, \\ F_{n+1} + F_{n+2} + 3, \quad n >1 \end{cases}$$
$$\therefore T(n) \geq F_n$$

$T(100) \approx 1.77 * 10^{21} \approx 56000$ years at 1GHz

Why is this so slow?
![Fibbonacci call stack](https://idawud.tech/assets/img/support/fibonacci_1.png)
For evaluating any number same number has to be evaluated multiple times eg f(3), f(2) has been repeatedly evaluated so many times, this increases exponentially when n is very large.

When we calculate manually we don't redo such calculation but rather use only the last two numbers.
"""

# ╔═╡ 0eea1918-4dc3-11eb-1484-b7f0db8d3a9a
begin
	function fibbo_efficient(n::Integer)
	    n >= 0 || throw(DomainError(n, "Domain: n ≥ 0"))
	    n <= 1 && return n

	    Fn = zeros(Integer, n + 1)
	    Fn[1:2] = [0, 1]  # you can only store two numbers

	    for i = 3:n+1
	        Fn[i] = Fn[i-1] + Fn[i-2]
	    end

	    return Fn[end]
	end

	@assert fibbo_efficient(7) == 13

	@benchmark fibbo_efficient(50)
end

# ╔═╡ b98372a2-4dc3-11eb-2956-69554b9ff7b8
md"""
$$T(n) = \begin{cases} 2, \qquad n \leq 1, \\ 2n + 3, \quad n >1 \end{cases}$$

$T(100) \approx 202$ Incredibly fast compared to above
"""

# ╔═╡ Cell order:
# ╟─12341f32-4dbf-11eb-1753-6fac96fd666a
# ╠═938f32d0-4dc2-11eb-29e1-19b4e28bf3c8
# ╟─10ef1d92-4dc0-11eb-3524-bd5cb089a31b
# ╟─09744cea-4dc0-11eb-11b0-c7c6b42683bc
# ╠═aae2f0aa-4dbf-11eb-3f50-39e262cc1845
# ╠═81d2d826-4dc0-11eb-2efb-9f5ef5161edd
# ╟─87f00406-4dc0-11eb-29b9-6b5c19d03b53
# ╠═a78d5570-4dc0-11eb-37a8-adc45a4964ee
# ╠═c0c78d52-4dc1-11eb-325d-81452aa6d218
# ╠═0dd9ea7c-4dc2-11eb-1b99-f7b120c21a70
# ╠═371f8004-4dc2-11eb-340e-c367fc0d718d
# ╟─5a9e6414-4dc2-11eb-2ed3-8f6abe2b7847
# ╠═61a2ce56-4dc2-11eb-3b58-df1a7f2bf4d6
# ╟─c3fdeae2-4dc2-11eb-0d57-598c117205c7
# ╠═e0401ae0-4dc2-11eb-2ba6-9d18e1915fa1
# ╟─f40e1e6e-4dc2-11eb-2c04-d7020e727f7f
# ╠═02317fca-4dc3-11eb-34c4-a5e009b09775
# ╟─1b865a74-4dc3-11eb-25ac-d56e25e7b4b6
# ╠═0eea1918-4dc3-11eb-1484-b7f0db8d3a9a
# ╟─b98372a2-4dc3-11eb-2956-69554b9ff7b8


text:
using BinaryBuilder

name="Cubature"
version=v"1.0.3"

# Collection of sources required to build Cubature
sources = [
    GitSource("https://github.com/stevengj/cubature.git",
              "a61dd32094bcf925fbc4a93d02e0f0f147f8b734"), # v1.0.3
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/cubature*
mkdir -p ${libdir}
${CC} ${LDFLAGS} -shared -fPIC -O3 hcubature.c pcubature.c -o ${libdir}/libcubature.${dlext}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libcubature", :libcubature),
]

# Dependencies that must be installed before this package can be built
dependencies = []

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
struct allFilms
    query::String
end
function (f::allFilms)(;
    after::Union{String,Missing,Nothing} = nothing,
    first::Union{Int32,Missing,Nothing} = nothing,
    before::Union{String,Missing,Nothing} = nothing,
    last::Union{Int32,Missing,Nothing} = nothing,
)
    q = (
        inp -> begin
            s = "query AllFilms(after: String, first: Int, before: String, last: Int) {\n    allFilms(after: \$after, first: \$first, before: \$before, last: \$last) {\n"
            s *= inp
            s *= "    }\n}\n"
        end
    )
    query = q(f.query)
    variables = Dict("after" => after, "first" => first, "before" => before, "last" => last)
    filter!((v -> !(isnothing(v[2]))), variables)
    return (; query, variables)
end

struct film
    query::String
end
function (f::film)(;
    id::Union{String,Missing,Nothing} = nothing,
    filmID::Union{String,Missing,Nothing} = nothing,
)
    q = (
        inp -> begin
            s = "query Film(id: ID, filmID: ID) {\n    film(id: \$id, filmID: \$filmID) {\n"
            s *= inp
            s *= "    }\n}\n"
        end
    )
    query = q(f.query)
    variables = Dict("id" => id, "filmID" => filmID)
    filter!((v -> !(isnothing(v[2]))), variables)
    return (; query, variables)
end

struct allPeople
    query::String
end
function (f::allPeople)(;
    after::Union{String,Missing,Nothing} = nothing,
    first::Union{Int32,Missing,Nothing} = nothing,
    before::Union{String,Missing,Nothing} = nothing,
    last::Union{Int32,Missing,Nothing} = nothing,
)
    q = (
        inp -> begin
            s = "query AllPeople(after: String, first: Int, before: String, last: Int) {\n    allPeople(after: \$after, first: \$first, before: \$before, last: \$last) {\n"
            s *= inp
            s *= "    }\n}\n"
        end
    )
    query = q(f.query)
    variables = Dict("after" => after, "first" => first, "before" => before, "last" => last)
    filter!((v -> !(isnothing(v[2]))), variables)
    return (; query, variables)
end

struct person
    query::String
end
function (f::person)(;
    id::Union{String,Missing,Nothing} = nothing,
    personID::Union{String,Missing,Nothing} = nothing,
)
    q = (
        inp -> begin
            s = "query Person(id: ID, personID: ID) {\n    person(id: \$id, personID: \$personID) {\n"
            s *= inp
            s *= "    }\n}\n"
        end
    )
    query = q(f.query)
    variables = Dict("id" => id, "personID" => personID)
    filter!((v -> !(isnothing(v[2]))), variables)
    return (; query, variables)
end

struct allPlanets
    query::String
end
function (f::allPlanets)(;
    after::Union{String,Missing,Nothing} = nothing,
    first::Union{Int32,Missing,Nothing} = nothing,
    before::Union{String,Missing,Nothing} = nothing,
    last::Union{Int32,Missing,Nothing} = nothing,
)
    q = (
        inp -> begin
            s = "query AllPlanets(after: String, first: Int, before: String, last: Int) {\n    allPlanets(after: \$after, first: \$first, before: \$before, last: \$last) {\n"
            s *= inp
            s *= "    }\n}\n"
        end
    )
    query = q(f.query)
    variables = Dict("after" => after, "first" => first, "before" => before, "last" => last)
    filter!((v -> !(isnothing(v[2]))), variables)
    return (; query, variables)
end

struct planet
    query::String
end
function (f::planet)(;
    id::Union{String,Missing,Nothing} = nothing,
    planetID::Union{String,Missing,Nothing} = nothing,
)
    q = (
        inp -> begin
            s = "query Planet(id: ID, planetID: ID) {\n    planet(id: \$id, planetID: \$planetID) {\n"
            s *= inp
            s *= "    }\n}\n"
        end
    )
    query = q(f.query)
    variables = Dict("id" => id, "planetID" => planetID)
    filter!((v -> !(isnothing(v[2]))), variables)
    return (; query, variables)
end

struct allSpecies
    query::String
end
function (f::allSpecies)(;
    after::Union{String,Missing,Nothing} = nothing,
    first::Union{Int32,Missing,Nothing} = nothing,
    before::Union{String,Missing,Nothing} = nothing,
    last::Union{Int32,Missing,Nothing} = nothing,
)
    q = (
        inp -> begin
            s = "query AllSpecies(after: String, first: Int, before: String, last: Int) {\n    allSpecies(after: \$after, first: \$first, before: \$before, last: \$last) {\n"
            s *= inp
            s *= "    }\n}\n"
        end
    )
    query = q(f.query)
    variables = Dict("after" => after, "first" => first, "before" => before, "last" => last)
    filter!((v -> !(isnothing(v[2]))), variables)
    return (; query, variables)
end

struct species
    query::String
end
function (f::species)(;
    id::Union{String,Missing,Nothing} = nothing,
    speciesID::Union{String,Missing,Nothing} = nothing,
)
    q = (
        inp -> begin
            s = "query Species(id: ID, speciesID: ID) {\n    species(id: \$id, speciesID: \$speciesID) {\n"
            s *= inp
            s *= "    }\n}\n"
        end
    )
    query = q(f.query)
    variables = Dict("id" => id, "speciesID" => speciesID)
    filter!((v -> !(isnothing(v[2]))), variables)
    return (; query, variables)
end

struct allStarships
    query::String
end
function (f::allStarships)(;
    after::Union{String,Missing,Nothing} = nothing,
    first::Union{Int32,Missing,Nothing} = nothing,
    before::Union{String,Missing,Nothing} = nothing,
    last::Union{Int32,Missing,Nothing} = nothing,
)
    q = (
        inp -> begin
            s = "query AllStarships(after: String, first: Int, before: String, last: Int) {\n    allStarships(after: \$after, first: \$first, before: \$before, last: \$last) {\n"
            s *= inp
            s *= "    }\n}\n"
        end
    )
    query = q(f.query)
    variables = Dict("after" => after, "first" => first, "before" => before, "last" => last)
    filter!((v -> !(isnothing(v[2]))), variables)
    return (; query, variables)
end

struct starship
    query::String
end
function (f::starship)(;
    id::Union{String,Missing,Nothing} = nothing,
    starshipID::Union{String,Missing,Nothing} = nothing,
)
    q = (
        inp -> begin
            s = "query Starship(id: ID, starshipID: ID) {\n    starship(id: \$id, starshipID: \$starshipID) {\n"
            s *= inp
            s *= "    }\n}\n"
        end
    )
    query = q(f.query)
    variables = Dict("id" => id, "starshipID" => starshipID)
    filter!((v -> !(isnothing(v[2]))), variables)
    return (; query, variables)
end

struct allVehicles
    query::String
end
function (f::allVehicles)(;
    after::Union{String,Missing,Nothing} = nothing,
    first::Union{Int32,Missing,Nothing} = nothing,
    before::Union{String,Missing,Nothing} = nothing,
    last::Union{Int32,Missing,Nothing} = nothing,
)
    q = (
        inp -> begin
            s = "query AllVehicles(after: String, first: Int, before: String, last: Int) {\n    allVehicles(after: \$after, first: \$first, before: \$before, last: \$last) {\n"
            s *= inp
            s *= "    }\n}\n"
        end
    )
    query = q(f.query)
    variables = Dict("after" => after, "first" => first, "before" => before, "last" => last)
    filter!((v -> !(isnothing(v[2]))), variables)
    return (; query, variables)
end

struct vehicle
    query::String
end
function (f::vehicle)(;
    id::Union{String,Missing,Nothing} = nothing,
    vehicleID::Union{String,Missing,Nothing} = nothing,
)
    q = (
        inp -> begin
            s = "query Vehicle(id: ID, vehicleID: ID) {\n    vehicle(id: \$id, vehicleID: \$vehicleID) {\n"
            s *= inp
            s *= "    }\n}\n"
        end
    )
    query = q(f.query)
    variables = Dict("id" => id, "vehicleID" => vehicleID)
    filter!((v -> !(isnothing(v[2]))), variables)
    return (; query, variables)
end

struct node
    query::String
end
"""
Fetches an object given its ID
"""
function (f::node)(id::String;)
    q = (inp -> begin
        s = "query Node(id: ID!) {\n    node(id: \$id) {\n"
        s *= inp
        s *= "    }\n}\n"
    end)
    query = q(f.query)
    variables = Dict("id" => id)
    filter!((v -> !(isnothing(v[2]))), variables)
    return (; query, variables)
end


text:
include("common.jl")

@testset "Basics" begin include("basics.jl") end
@testset "Exploration" begin include("ifsandbuts.jl") end
@testset "Hessian sparsity" begin include("hessian.jl") end
@testset "Paraboloid example" begin include("paraboloid.jl") end


text:
using Printf
using FiniteDiff
using ForwardDiff
using KernelAbstractions
using LinearAlgebra
using Random
using SparseArrays
using Test
using TimerOutputs
using ExaPF
import ExaPF: PowerSystem, AutoDiff

const PS = PowerSystem
const KA = KernelAbstractions

@testset "Compute reduced gradient on CPU" begin
    @testset "Case $case" for case in ["case9.m", "case30.m"]
        datafile = joinpath(dirname(@__FILE__), "..", "..", "data", case)
        tolerance = 1e-8
        pf = PS.PowerNetwork(datafile)

        polar = PolarForm(pf, CPU())
        cache = ExaPF.get(polar, ExaPF.PhysicalState())
        ExaPF.init_buffer!(polar, cache)

        xk = ExaPF.initial(polar, State())
        u = ExaPF.initial(polar, Control())

        jx = AutoDiff.Jacobian(polar, ExaPF.power_balance, State())
        ju = AutoDiff.Jacobian(polar, ExaPF.power_balance, Control())
        ∂obj = ExaPF.AdjointStackObjective(polar)
        pbm = AutoDiff.TapeMemory(ExaPF.active_power_generation, ∂obj, nothing)

        # solve power flow
        conv = powerflow(polar, jx, cache, NewtonRaphson(tol=1e-12))
        # No need to recompute ∇gₓ
        ∇gₓ = jx.J
        ∇gᵤ = AutoDiff.jacobian!(polar, ju, cache)
        # test jacobian wrt x
        ∇gᵥ = AutoDiff.jacobian!(polar, jx, cache)
        @test isequal(∇gₓ, ∇gᵥ)

        # Test with Matpower's Jacobian
        V = cache.vmag .* exp.(im * cache.vang)
        J = ExaPF.matpower_jacobian(polar, State(), ExaPF.power_balance, V)
        @test isapprox(∇gₓ, J)

        # Test gradients
        @testset "Reduced gradient" begin
            # Refresh cache with new values of vmag and vang
            ExaPF.update!(polar, PS.Generators(), PS.ActivePower(), cache)
            # We need uk here for the closure
            uk = copy(u)
            ExaPF.gradient_objective!(polar, pbm, cache)
            ∇fₓ = ∂obj.∇fₓ
            ∇fᵤ = ∂obj.∇fᵤ

            ## ADJOINT
            # lamba calculation
            λk  = -(∇gₓ') \ ∇fₓ
            grad_adjoint = ∇fᵤ + ∇gᵤ' * λk
            # ## DIRECT
            S = - inv(Array(∇gₓ)) * ∇gᵤ
            grad_direct = ∇fᵤ + S' * ∇fₓ
            @test isapprox(grad_adjoint, grad_direct)

            # Compare with finite difference
            function reduced_cost(u_)
                # Ensure we remain in the manifold
                ExaPF.transfer!(polar, cache, u_)
                convergence = powerflow(polar, jx, cache, NewtonRaphson(tol=1e-14))
                ExaPF.update!(polar, PS.Generators(), PS.ActivePower(), cache)
                return ExaPF.cost_production(polar, cache)
            end

            grad_fd = FiniteDiff.finite_difference_gradient(reduced_cost, uk)
            @test isapprox(grad_fd, grad_adjoint, rtol=1e-4)
        end
        @testset "Gradient of line-flow constraints" begin
            # Adjoint of flow_constraints()
            nbus = length(cache.vmag)
            M = typeof(u)
            m = ExaPF.size_constraint(polar, ExaPF.flow_constraints)
            x = ExaPF.xzeros(M, 2 * nbus)
            x[1:nbus] .= cache.vmag
            x[1+nbus:2*nbus] .= cache.vang
            bus_gen = polar.indexing.index_generators
            VI = typeof(bus_gen)

            ## Example with using sum as a sort of lumping of all constraints
            function sum_constraints(x)
                VT = typeof(x)
                # Needed for ForwardDiff to have a cache with the right active type VT
                adcache = ExaPF.PolarNetworkState{VI, VT}(
                    cache.vmag, cache.vang, cache.pinj, cache.qinj,
                    cache.pg, cache.qg, cache.balance, cache.dx, bus_gen,
                )
                adcache.vmag .= x[1:nbus]
                adcache.vang .= x[1+nbus:2*nbus]
                g = VT(undef, m) ; fill!(g, 0)
                ExaPF.flow_constraints(polar, g, adcache)
                return sum(g)
            end
            adgradg = ForwardDiff.gradient(sum_constraints,x)
            fdgradg = FiniteDiff.finite_difference_gradient(sum_constraints,x)
            ## We pick sum() as the reduction function. This could be a mask function for active set or some log(x) for lumping.
            m_flows = ExaPF.size_constraint(polar, ExaPF.flow_constraints)
            weights = ones(m_flows)
            gradg = ExaPF.xzeros(M, 2 * nbus)
            ExaPF.flow_constraints_grad!(polar, gradg, cache, weights)
            @test isapprox(adgradg, fdgradg)
            # TODO: The gradient is slightly off with the handwritten adjoint
            @test_broken isapprox(gradg, fdgradg)
        end
    end
end


text:
using CPUTime
using Test

function add_and_sleep()
    x = 0
    for i in 1:10_000_000
        x += i
    end
    sleep(1)
end

function time_difference()
    dt = @elapsed begin
        CPUdt = @CPUelapsed add_and_sleep()
    end
    return dt - CPUdt
end

time_difference()  # compilation overhead should not affect results, but just in case
time_diff = time_difference()
eval(Meta.parse("@test abs($time_diff - 1.0) <= .1"))


text:
module Test_Prismatic_xAxis

using  Modia3D
import Modia3D.ModiaMath

@assembly FallingBall(;h=1.0) begin
   world  = Object3D()
   sphere = Object3D( Modia3D.Solid(Modia3D.SolidSphere(0.1), "Aluminium") )

   # Constrain sphere movement (initial placement at position [0,h,0])
   prismatic = Modia3D.Prismatic(world, sphere, axis=1, s_start=h)
end

gravField = Modia3D.UniformGravityField(n=[-1,0,0])
fallBall = FallingBall(h=3.5, sceneOptions=Modia3D.SceneOptions(gravityField=gravField,visualizeFrames=true, defaultFrameLength=0.3))
simulationModel = Modia3D.SimulationModel( fallBall, stopTime=1.0 )
result          = ModiaMath.simulate!(simulationModel)
ModiaMath.plot(result, ("prismatic.s", "prismatic.v"))

println("... success of Test_Prismatic_xAxis.jl!")

end


text:
"""
    marchonintime(W0,Z,B,I)

Solve by marching-on-in-time the causal convolution problem defined by `(W0,Z,B)`
up to timestep `I`. Here, `Z` is an array of order 3 that contains a discretisation
of a time translation invariant retarded potential operator. `W0` is the inverse of
the slice `Z[:,:,1]`.
"""
function marchonintime(W0,Z,B,I)
    T = eltype(W0)
    M,N = size(Z)
    @assert M == size(B,1)
    x = zeros(T,N,I)
    for i in 1:I
        b = B[:,i] - convolve(Z,x,i,2)
        x[:,i] += W0 * b
        (i % 10 == 0) && print(i, "[", I, "] - ")
    end
    return x
end

function marchonintime(iZ0, Z::ConvOp, B, Nt)

    T = eltype(iZ0)
    Ns = size(Z,1)
    x = zeros(T,Ns,Nt)
    csx = zeros(T,Ns,Nt)
    y = zeros(T,Ns)

    todo, done, pct = Nt, 0, 0
    for i in 1:Nt
        fill!(y,0)
        convolve!(y, Z, x, csx, i, 2, Nt)
        y .*= -1
        y .+= B[:,i]
        # @show norm(B[:,i])

        x[:,i] .+= iZ0 * y
        if i > 1
            csx[:,i] .= csx[:,i-1] .+ x[:,i]
        else
            csx[:,i] .= x[:,i]
        end

        done += 1
        new_pct = round(Int, done / todo * 100)
        new_pct > pct+9 && (println("[$new_pct]"); pct=new_pct)
    end
    x
end

using BlockArrays

function convolve(Z::BlockArray, x, i, j_start)
    # ax1 = axes(Z,1)
    ax2 = axes(Z,2)
    T = eltype(eltype(Z))
    y = PseudoBlockVector{T}(undef,blocklengths(axes(Z,1)))
    fill!(y,0)
    for I in blockaxes(Z,1)
        for J in blockaxes(Z,2)
            xJ = view(x, ax2[J], :)
            try
                ZIJ = Z[I,J].banded
                y[I] .+= convolve(ZIJ, xJ, i, j_start)
            catch
                @info "Skipping unassigned block."
                continue
            end
        end
    end
    return y
end

function convolve!(y,Z::BlockArray, x, csx, i, j_start, j_stop)
    ax1 = axes(Z,1)
    ax2 = axes(Z,2)
    T = eltype(eltype(Z))
    # y = PseudoBlockVector{T}(undef,blocklengths(axes(Z,1)))
    fill!(y,0)
    for I in blockaxes(Z,1)
        for J in blockaxes(Z,2)
            xJ = view(x, ax2[J], :)
            csxJ = view(csx, ax2[J], :)
            try
                ZIJ = Z[I,J].banded
                # y[I] .+= convolve(ZIJ, xJ, i, j_start)
                yI = view(y, ax1[I])
                convolve!(yI, ZIJ, xJ, csxJ, i, j_start, j_stop)
            catch
                @info "Skipping unassigned block."
                continue
            end
        end
    end
    return y
end

function marchonintime(W0,Z::BlockArray,B,I)

    T = eltype(W0)
    M,N = size(W0)
    @assert M == size(B,1)

    x = zeros(T,N,I)
    y = zeros(T,N)
    csx = zeros(T,N,I)

    for i in 1:I
        R = [ B[j][i] for j in 1:N ]
        # @show norm(R)
        k_start = 2
        k_stop = I

        fill!(y,0)
        convolve!(y,Z,x,csx,i,k_start,k_stop)
        b = R - y
        x[:,i] .+= W0 * b
        if i > 1
            csx[:,i] .= csx[:,i-1] .+ x[:,i]
        else
            csx[:,i] .= x[:,i]
        end

        (i % 10 == 0) && print(i, "[", I, "] - ")
    end
    return x
end


text:
data = split(strip(read(joinpath(@__DIR__, "./2.txt"), String)), '\n')

#part 1
position = Dict(:horizontal => 0, :depth => 0)

command = Dict("forward" => (:horizontal,+), "down" => (:depth,+), "up" => (:depth,-))

function executeCommands(data)
  global position,command
  direction,distance = split(data,' ')
  position[command[direction][1]] = command[direction][2](position[command[direction][1]],parse(Int,distance))
end

executeCommands.(data)
position[:horizontal]*position[:depth]

# part 2
position = Dict(:horizontal => 0, :depth => 0, :aim => 0)

command = Dict("forward" => (:horizontal,+), "down" => (:aim,+), "up" => (:aim,-))

function executeCommands2(data)
  global position
  direction,distance = split(data,' ')
  d = parse(Int,distance)
  position[command[direction][1]] = command[direction][2](position[command[direction][1]],d)
  if direction == "forward"
    position[:depth] = position[:depth] + d * position[:aim]
  end
end

executeCommands2.(data)
position[:horizontal]*position[:depth]


text:
module PlutoCanvasPlot

using UUIDs

include("api.jl")
include("internals.jl")

export command!
export CanvasPlot,  polygon!,linecolor!, fillcolor!
export textcolor!,textsize!,text!
export polyline!,linecolor!
export polygon!,fillcolor!
export axis!

end # module


text:
using GuillotineModels.Data

const SLOPP_HEADER = """
***2D Rectangular Problem***
***Instances for the Single Large Object Placement Problem (SLOPP)***
Input parameter file: SLOPP_parameters.txt
***************************************************************************************************************
Total number of instances
LargeObject.Length      LargeObject.Width
Number of different item types (i)
Item[i].Length  Item[i].Width   Item[i].LowerBoundDemand        Item[i].UpperBoundDemand        Item[i].Value
***************************************************************************************************************"""

const SLOPP_BODY_1 = """
50      60
28
21      13      4       4       264
6       45      1       2       191
5       60      2       4       264
30      39      1       1       795
14      42      5       5       311
29      17      2       3       314
10      9       4       5       73
23      31      2       3       567
16      22      4       5       181
14      26      1       1       182
24      24      2       5       476
12      56      2       2       504
12      49      3       5       337
12      47      1       3       413
18      32      4       4       520
26      42      1       2       1064
14      16      2       3       163
8       38      4       5       219
6       58      3       4       293
21      40      2       3       791
13      52      2       5       479
8       17      1       1       73
29      56      1       5       1449
18      9       1       5       101
27      27      1       1       596
6       50      1       2       255
38      21      3       4       571
27      35      2       3       697"""

SLOPP_BODY_2 = """
83      67
10
48      28      4       4       1017
68      44      1       1       2458
14      58      4       5       519
5       43      1       1       205
63      8       1       2       349
49      56      1       5       1447
6       22      3       3       95
12      30      4       4       239
50      42      4       5       1763
17      34      4       4       325 """

const SLOPP_INST_1 = SLOPP{Int,Int,Int}(
	50, 60,
	[21, 6, 5, 30, 14, 29, 10, 23, 16, 14, 24, 12, 12, 12, 18, 26, 14, 8, 6, 21, 13, 8, 29, 18, 27, 6, 38, 27],
	[13, 45, 60, 39, 42, 17, 9, 31, 22, 26, 24, 56, 49, 47, 32, 42, 16, 38, 58, 40, 52, 17, 56, 9, 27, 50, 21, 35],
	[4, 1, 2, 1, 5, 2, 4, 2, 4, 1, 2, 2, 3, 1, 4, 1, 2, 4, 3, 2, 2, 1, 1, 1, 1, 1, 3, 2],
	[4, 2, 4, 1, 5, 3, 5, 3, 5, 1, 5, 2, 5, 3, 4, 2, 3, 5, 4, 3, 5, 1, 5, 5, 1, 2, 4, 3],
	[264, 191, 264, 795, 311, 314, 73, 567, 181, 182, 476, 504, 337, 413, 520, 1064, 163, 219, 293, 791, 479, 73, 1449, 101, 596, 255, 571, 697]
)

const SLOPP_INST_2 = SLOPP{Int64,Int64,Int64}(
	83, 67,
	[48, 68, 14, 5, 63, 49, 6, 12, 50, 17],
	[28, 44, 58, 43, 8, 56, 22, 30, 42, 34],
	[4, 1, 4, 1, 1, 1, 3, 4, 4, 4],
	[4, 1, 5, 1, 2, 5, 3, 4, 5, 4],
	[1017, 2458, 519, 205, 349, 1447, 95, 239, 1763, 325]
)

const SLOPP_WITHOUT_ERRORS = [
(
"File has a single instance.",
SLOPP{Int,Int,Int}[SLOPP_INST_1],
"""
$(SLOPP_HEADER)
1
$(SLOPP_BODY_1)
"""
),
(# Two instances.
"File has two instances.",
SLOPP{Int,Int,Int}[SLOPP_INST_1, SLOPP_INST_2],
"""
$(SLOPP_HEADER)
2
$(SLOPP_BODY_1)
$(SLOPP_BODY_2)
"""
)
]

const SSSCSP_HEADER = """
***2D Rectangular Problem***
***Instances for the Single Stock Size Cutting Stock Problem (SSSCSP)***
Input parameter file: SSSCSP_parameters.txt
****************************************************************************************************
Total number of instances
LargeObject.Length      LargeObject.Width
Number of different item types (i)
Item[i].Length  Item[i].Width   Item[i].Demand
*****************************************************************************************************"""

const SLOPP_WITH_WARNINGS = [
(
"Warning because of wrong header.",
(:warn, r"Could not match the header pattern for format:"),
"""
$(SSSCSP_HEADER)
1
$(SLOPP_BODY_1)
"""
),
(
"Warning because of extra lines after parsing all instances.",
(:warn, """The file/string indicated 1 instance(s), after consuming those some lines remained. Ignoring any further lines. The last parsed line was \"$(last(split(SLOPP_BODY_1, isequal('\n'); keepempty = false)))\" (line $(count("\n", SLOPP_HEADER) + 3 + count("\n", SLOPP_BODY_1)) considering only non-empty lines)."""),
"""
$(SLOPP_HEADER)
1
$(SLOPP_BODY_1)
Additional non-empty line that triggers warning.
"""
)
]

const SLOPP_WITH_ERRORS = [
(
"Incomplete header (not two lines of asterisks).",
GenericParseError(CPG_SLOPP{Int64,Int64,Int64}(), "Parser expects to find two lines only with asterisks; found 1 line(s) like this. As the data should be after the second of these lines, the parser did not start gathering data."),
"""
$(join(split(SLOPP_HEADER, isequal('\n'); keepempty = false)[1:end-1], "\n"))
2
$SLOPP_BODY_1
$SLOPP_BODY_2
"""
),
(
"Only the header is present (no data).",
GenericParseError(CPG_SLOPP{Int64,Int64,Int64}(), "Expected non-empty line number 10 to have 'the number of problem instances' but the file/string finished before this line."),
"""
$(SLOPP_HEADER)
"""
),
(
"Number of instances line was skipped (has wrong number of tokens).",
GenericParseError(CPG_SLOPP{Int64,Int64,Int64}(), "Expected 1 token(s) (the number of problem instances) from \"50      60\" (non-empty line number 10) got 2 token(s) instead."),
"""
$(SLOPP_HEADER)
$(SLOPP_BODY_1)
"""
),
(
"Number of instances line has an invalid token.",
GenericParseError(CPG_SLOPP{Int64,Int64,Int64}(), "Expected 'Int64' (the number of problem instances) from \"some_garbage\" (non-empty line number 10) but could not parse some of the tokens to their respective types."),
"""
$(SLOPP_HEADER)
some_garbage
$(SLOPP_BODY_1)
"""
),
(
"Provided number of instances is larger than the real number.",
GenericParseError(CPG_SLOPP{Int64,Int64,Int64}(), "Expected non-empty line number 41 to have 'L and W' but the file/string finished before this line."),
"""
$(SLOPP_HEADER)
2
$(SLOPP_BODY_1)
"""
),
(
"The original plate dimensions are not two tokens.",
GenericParseError(CPG_SLOPP{Int64,Int64,Int64}(), "Expected 2 token(s) (L and W) from \"28\" (non-empty line number 11) got 1 token(s) instead."),
"""
$(SLOPP_HEADER)
1
$(join(split(SLOPP_BODY_1, isequal('\n'); keepempty = false)[2:end], "\n"))
"""
),
(
"The original plate dimensions are two invalid tokens.",
GenericParseError(CPG_SLOPP{Int64,Int64,Int64}(), "Expected 'Int64 Int64' (L and W) from \"some garbage\" (non-empty line number 11) but could not parse some of the tokens to their respective types."),
"""
$(SLOPP_HEADER)
1
some garbage
$(join(split(SLOPP_BODY_1, isequal('\n'); keepempty = false)[2:end], "\n"))
"""
),
(
"The file ends before the number of small objects inside an instance.",
GenericParseError(CPG_SLOPP{Int64,Int64,Int64}(), "Expected non-empty line number 12 to have 'the number of small objects' but the file/string finished before this line."),
"""
$(SLOPP_HEADER)
1
100 100
"""
),
(
"Number of pieces is skipped (i.e., next line is a piece).",
GenericParseError(CPG_SLOPP{Int64,Int64,Int64}(), "Expected 1 token(s) (the number of small objects) from \"100 100\" (non-empty line number 12) got 2 token(s) instead."),
"""
$(SLOPP_HEADER)
1
100 100
100 100
"""
),
#=
( # Third line not a number.
GenericParseError(Val(:CPG_SLOPP), "dummy"),
"""
$(SLOPP_HEADER)
"""
),
( # Not enough lines for the small objects.
GenericParseError(Val(:CPG_SLOPP), "dummy"),
"""
$(SLOPP_HEADER)
"""
),
( # Not five tokens in a small object line.
GenericParseError(Val(:CPG_SLOPP), "dummy"),
"""
$(SLOPP_HEADER)
"""
),
=#
]

#= used for debug
function shallow_obj_diff(x :: T, y :: T) where {T}
	println("My test 2.")
	for field in fieldnames(T)
		@show field
		xf, yf = getfield(x, field), getfield(y, field)
		@show typeof.((xf, yf))
		if xf != yf
			for (vxf, vyf) in zip(xf, yf)
				@test vxf == vyf
			end
			println("The values in the fields $field are distinct:")
			#show(IOContext(stdout, :limit => false), xf)
			dump(IOContext(stdout, :limit => false), xf)
			println()
			#show(IOContext(stdout, :limit => false), yf)
			dump(IOContext(stdout, :limit => false), yf)
			println()
		end
	end

	return nothing
end
=#

function test_SLOPP()
	@testset "Valid SLOPP dataset, no warnings." begin
		for (case, expected, s) in SLOPP_WITHOUT_ERRORS
			# If the instance is perfectly normal then no message should be printed.
			@testset "$case" begin
				obtained = @test_logs read_from_string(Val(:CPG_SLOPP), s)
				#show(IOContext(stdout, :limit => false), obtained)
				@test expected == obtained
			end
		end
	end
	@testset "SLOPP datasets that trigger warnings." begin
		for (case, expected, s) in SLOPP_WITH_WARNINGS
			@testset "$case" begin
				@test_logs expected read_from_string(Val(:CPG_SLOPP), s)
			end
		end
	end
	@testset "SLOPP datasets that trigger exceptions." begin
		for (case, exception, s) in SLOPP_WITH_ERRORS
			@testset "$case" begin
				@test_throws exception read_from_string(Val(:CPG_SLOPP), s)
			end
		end
	end
	# TODO: tests with different number types in each field
end

const ODPW_HEADER = """
***2D Rectangular Problem***
***Problem tests for the Open Dimension Problem (ODP/W)***
Input parameter file: ODPW_parameters.txt
***********************************************************************
Total number of instances
Number of different large objects (j)
LargeObject[j].Width    LargeObject[j].Available        LargeObject[j].Value
Number of different item types (i)
Item[i].Length  Item[i].Width   Item[i].Demand
***********************************************************************"""

const MHLOPPW_HEADER = """
***2D Rectangular Problem***
***Instances for the Multiple Heterogeneous Large Object Placement Problem (MHLOPP/W)***
Input parameter file: MHLOPPW_parameters.txt
***************************************************************************************************************
Total number of instances
Number of different large objects (j)
LargeObject[j].Length   LargeObject[j].Width    LargeObject[j].Available
Number of different item types (i)
Item[i].Length  Item[i].Width   Item[i].LowerBoundDemand        Item[i].UpperBoundDemand        Item[i].Value
***************************************************************************************************************"""

const SSSCSP_BODY_1 = """
82      138
20
81      7       5
6       128     2
7       127     5
39      58      1
9       130     2
28      126     3
77      87      4
6       80      1
14      80      2
35      109     2
49      62      1
5       96      1
47      6       2
32      7       2
40      109     5
33      80      5
55      114     2
10      61      5
34      138     1
61      99      5"""

const SSSCSP_BODY_2 = """
144     94
23
38      88      3
58      27      1
88      32      5
122     19      5
7       16      2
132     44      1
104     44      2
15      52      1
60      86      2
33      91      1
62      13      1
122     11      4
88      86      2
37      27      5
140     13      2
60      35      3
38      37      1
58      8       1
20      63      1
109     19      5
6       21      2
139     7       1
6       47      1"""

const SSSCSP_INST_1 = SSSCSP{Int64,Int64,Int64}(
	82, 138,
	[81, 6, 7, 39, 9, 28, 77, 6, 14, 35, 49, 5, 47, 32, 40, 33, 55,
		10, 34, 61],
	[7, 128, 127, 58, 130, 126, 87, 80, 80, 109, 62, 96, 6, 7, 109,
		80, 114, 61, 138, 99],
	[5, 2, 5, 1, 2, 3, 4, 1, 2, 2, 1, 1, 2, 2, 5, 5, 2, 5, 1, 5]
)

const SSSCSP_INST_2 = SSSCSP{Int64,Int64,Int64}(
	144, 94,
	[38, 58, 88, 122, 7, 132, 104, 15, 60, 33, 62, 122, 88, 37, 140,
		60, 38, 58, 20, 109, 6, 139, 6],
	[88, 27, 32, 19, 16, 44, 44, 52, 86, 91, 13, 11, 86, 27, 13, 35,
		37, 8, 63, 19, 21, 7, 47],
	[3, 1, 5, 5, 2, 1, 2, 1, 2, 1, 1, 4, 2, 5, 2, 3, 1, 1, 1, 5, 2, 1, 1]
)

const ODPW_BODY_1 = """
2
54      1       33
93      1       48
26
116     86      2
82      70      1
6       45      4
5       36      1
59      56      4
149     71      3
8       54      2
121     81      5
28      13      5
40      52      5
150     35      1
150     59      5
111     16      4
146     29      3
132     42      1
30      35      1
74      79      5
6       80      3
81      32      5
76      29      1
150     16      3
52      30      5
84      41      5
5       62      5
140     20      1
118     74      4"""

const ODPW_BODY_2 = """
9
133     1       104
147     1       133
71      1       36
102     1       78
145     1       108
144     1       141
87      1       76
95      1       74
57      1       54
24
142     22      1
123     20      1
30      7       2
5       58      5
24      24      2
146     53      1
117     36      2
77      21      4
142     49      5
90      22      5
57      21      4
26      42      3
150     20      1
38      52      1
111     51      5
26      42      2
19      11      5
61      21      5
57      20      5
8       22      1
48      28      1
134     40      2
67      67      5
111     26      4"""

const ODPW_INST_1 = ODPW{Int64,Int64,Int64}(
	[54, 93], [1, 1], [33, 48],
	[116, 82, 6, 5, 59, 149, 8, 121, 28, 40, 150, 150, 111, 146, 132, 30,
		74, 6, 81, 76, 150, 52, 84, 5, 140, 118],
	[86, 70, 45, 36, 56, 71, 54, 81, 13, 52, 35, 59, 16, 29, 42, 35, 79,
		80, 32, 29, 16, 30, 41, 62, 20, 74],
	[2, 1, 4, 1, 4, 3, 2, 5, 5, 5, 1, 5, 4, 3, 1, 1, 5, 3, 5, 1, 3, 5, 5,
		5, 1, 4]
)

const ODPW_INST_2 = ODPW{Int64,Int64,Int64}(
	[133, 147, 71, 102, 145, 144, 87, 95, 57],
	[1, 1, 1, 1, 1, 1, 1, 1, 1], [104, 133, 36, 78, 108, 141, 76, 74, 54],
	[142, 123, 30, 5, 24, 146, 117, 77, 142, 90, 57, 26, 150, 38, 111, 26,
		19, 61, 57, 8, 48, 134, 67, 111],
	[22, 20, 7, 58, 24, 53, 36, 21, 49, 22, 21, 42, 20, 52, 51, 42, 11,
		21, 20, 22, 28, 40, 67, 26],
	[1, 1, 2, 5, 2, 1, 2, 4, 5, 5, 4, 3, 1, 1, 5, 2, 5, 5, 5, 1, 1, 2, 5, 4]
)

const MHLOPPW_BODY_1 = """
2
138     96      39
133     125     21
5
5       83      0       3       415
58      8       0       5       464
127     9       0       3       1143
47      92      0       5       4324
77      85      0       5       6545"""

const MHLOPPW_BODY_2 = """
10
104     147     21
124     126     32
96      124     25
105     124     28
114     130     28
135     74      26
140     118     29
91      143     27
147     137     33
130     133     24
1
96      5       0       1       480"""

MHLOPPW_INST_1 = MHLOPPW{Int64,Int64,Int64}(
	[138, 133], [96, 125], [39, 21], [5, 58, 127, 47, 77],
	[83, 8, 9, 92, 85], [0, 0, 0, 0, 0], [3, 5, 3, 5, 5],
	[415, 464, 1143, 4324, 6545]
)

MHLOPPW_INST_2 = MHLOPPW{Int64,Int64,Int64}(
	[104, 124, 96, 105, 114, 135, 140, 91, 147, 130],
	[147, 126, 124, 124, 130, 74, 118, 143, 137, 133],
	[21, 32, 25, 28, 28, 26, 29, 27, 33, 24],
	[96], [5], [0], [1], [480]
)

const NON_SLOPP_NO_ERROR = Dict{Symbol, Vector{Tuple{String, String, Any}}}(
:ODPW => Tuple{String, String, Any}[
(
"File has a single instance.",
"""
$ODPW_HEADER
1
$ODPW_BODY_1
""",
ODPW{Int64,Int64,Int64}[ODPW_INST_1]
),
(
"File has two instances.",
"""
$ODPW_HEADER
2
$ODPW_BODY_1
$ODPW_BODY_2
""",
ODPW{Int64,Int64,Int64}[ODPW_INST_1, ODPW_INST_2]
),
], # END ODPW
:MHLOPPW => Tuple{String, String, Any}[
(
"File has a single instance.",
"""
$MHLOPPW_HEADER
1
$MHLOPPW_BODY_1
""",
MHLOPPW{Int64,Int64,Int64}[MHLOPPW_INST_1]
),
(
"File has two instances.",
"""
$MHLOPPW_HEADER
2
$MHLOPPW_BODY_1
$MHLOPPW_BODY_2
""",
MHLOPPW{Int64,Int64,Int64}[MHLOPPW_INST_1, MHLOPPW_INST_2]
),
], # END MHLOPPW
:SSSCSP => Tuple{String, String, Any}[
(
"File has a single instance.",
"""
$SSSCSP_HEADER
1
$SSSCSP_BODY_1
""",
SSSCSP{Int64,Int64,Int64}[SSSCSP_INST_1]
),
(
"File has two instances.",
"""
$SSSCSP_HEADER
2
$SSSCSP_BODY_1
$SSSCSP_BODY_2
""",
SSSCSP{Int64,Int64,Int64}[SSSCSP_INST_1, SSSCSP_INST_2]
),
] # END SSSCSP
) # END DICT

# Checking every possible exception for every possible type is very
# time consuming. The tests over SLOPP cover most exceptional cases.
# Here we test if the reading of valid instances of other types is
# working.
function test_non_SLOPP()
	for problem_type in (:ODPW, :MHLOPPW, :SSSCSP)
		@testset "Valid $problem_type dataset, no warnings." begin
			for (case, s, expected) in NON_SLOPP_NO_ERROR[problem_type]
				# If the instance is perfectly normal then no message should be printed.
				@testset "$case" begin
					format_symbol = Symbol("CPG_$problem_type")
					obtained = @test_logs read_from_string(Val(format_symbol), s)
					show(IOContext(stdout, :limit => false), obtained)
					@test expected == obtained
				end
			end
		end
	end
end

test_SLOPP()
test_non_SLOPP()


text:
module LinAlg
include("./graphmatrices.jl")
end

include("./spectral.jl")
include("./nonbacktracking.jl")

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
function render_duel(img1, img2)
    gui = ImageView.imshow_gui((300, 300), (1, 2))
    canvases = gui["canvas"]
    ImageView.imshow(canvases[1,1], img1)
    ImageView.imshow(canvases[1,2], img2)
    Gtk.showall(gui["window"])

    options = ["right image", "left image"]
    menu    = TerminalMenus.RadioMenu(options, pagesize=2)
    choice  = TerminalMenus.request("Choose image with better quality:", menu)
    if(choice == -1)
        ErrorException("Cancelled choice")
    end
    ImageView.closeall()
    return choice
end

function pyramid_diffusion(img, x::Vector)
    factor = 2
    pyr    = Images.gaussian_pyramid(img, 4, factor, x[1])
    pyr    = laplacian_pyramid(pyr, factor)
    pyr[1] = diffusion(pyr[1], x[2:3]..., Int64(x[4]))
    pyr[2] = diffusion(pyr[2], x[5:6]..., Int64(x[7]))
    pyr[3] = diffusion(pyr[3], x[8:9]..., Int64(x[10]))
    pyr[4] = diffusion(pyr[4], x[11:12]..., Int64(x[13]))
    img    = synthesize_pyramid(pyr, factor)
end

function lerp(x::Real, low::Real, high::Real)
    return x * (high - low) + low;
end

function pyramid_diffusion_duel(x1::Vector, x2::Vector)
    function transform_domain(x::Vector)
        res    = clamp.(x, 0, 1)
        res[1] = 2^lerp(x[1], 0, 5)

        res[2] = lerp(x[2], 0.0001, 0.25)
        res[3] = 2^lerp(x[3], -8, 0)
        res[4] = floor(2^lerp(x[4], 0, 8))

        res[5] = lerp(x[5], 0.0001, 0.25)
        res[6] = 2^lerp(x[6], -8, 0)
        res[7] = floor(2^lerp(x[7], 0, 8))

        res[8] = lerp(x[2], 0.0001, 0.25)
        res[9] = 2^lerp(x[3], -8, 0)
        res[10] = floor(2^lerp(x[4], 0, 8))

        res[11] = lerp(x[11], 0.0001, 0.25)
        res[12] = 2^lerp(x[12], -8, 0)
        res[13] = floor(2^lerp(x[13], 0, 8))
        res
    end

    img = TestImages.testimage("lena_gray_256.tif")
    img = clamp.(img + rand(eltype(img), size(img))*0.2, 0, 1)

    x1 = transform_domain(x1)
    x2 = transform_domain(x2)
    img1 = pyramid_diffusion(img, x1)
    img2 = pyramid_diffusion(img, x2)
    render_duel(img1, img2)
end

function diffusion_duel(x1::Vector, x2::Vector)
    function transform_domain(x::Vector)
        res    = deepcopy(x)
        res[1] = 2^lerp(x[1], 0, 5)
        res[2] = lerp(x[2], 0.0001, 0.249)
        res[3] = 2^lerp(x[3], -8, 0)
        res[4] = floor(2^lerp(x[4], 0, 8))
        res
    end

    img = TestImages.testimage("lena_gray_256.tif")
    img = clamp.(img + rand(eltype(img), size(img))*0.2, 0, 1)

    x1 = transform_domain(x1)
    x2 = transform_domain(x2)
    img1 = diffusion(img, x1[2:3]..., Int64(x1[4]))
    img2 = diffusion(img, x2[2:3]..., Int64(x2[4]))
    render_duel(img1, img2)
end

#pairwise_prefbo(f, dims, warmup_steps)


text:
using Random: Random, AbstractRNG, default_rng, randstring

abstract type AbstractGenerator{T} end

struct Generator{T,F} <: AbstractGenerator{T}
    gen::F
end
Generator{T}(g) where T = Generator{T,typeof(g)}(g)
Generator(el::T) where T = Generator{T}((rng)->generate(rng, el))
Generator(::Type{T}) where T = Generator{T}((rng)->generate(rng, T))

# unsure if this is really a good idea
function Base.iterate(g::Generator, state=nothing)
    element = generate(default_rng(), g)
    element === nothing && return nothing
    (element, nothing)
end
Base.IteratorEltype(::Type{<:Generator}) = Base.HasEltype()
Base.IteratorSize(::Type{<:Generator}) = Base.IsInfinite()
Base.eltype(::Type{Generator{T,F}}) where {T,F} = T

generate(rng, g::Generator) = g.gen(rng)

# fallback generator if no RNG is passed in
generate(x) = generate(Random.default_rng(), x)

struct Word
    hi::UInt
end
generate(rng, w::Word) = Word(rand(rng, 0:w.hi))
generate(rng, ::Type{Word}) = Word(rand(rng, UInt))
const mWord(w) = Generator(Word(w))
shrink(w::Word) = map(Word, shrink(w.hi))

#######################
# type based generation
#######################

generate(rng, ::Type{T}) where {T <: Number} = rand(rng, T) # numbers
generate(rng, ::Type{Bool}) = rand(rng, Bool)
generate(rng, ::Type{NTuple{N,T}}) where {N,T} = ntuple(_ -> generate(rng, T), N)

########################
# value based generation
########################

generate(rng, t::T) where {T <: Unsigned} = rand(zero(T):t) # numbers
generate(rng, t::T) where {T <: Signed} = rand(-t:t) # numbers
generate(rng, t::NTuple{N,T}) where {N,T} = ntuple(i -> generate(rng, t[i]), N)


text:
function init_pot_harm_3d( LF, ω, center )
    Npoints = LF.Nx * LF.Ny * LF.Nz
    Vpot = zeros( Float64,Npoints )
    for ip = 1:Npoints
        Vpot[ip] = 0.5 * ω^2 * norm(LF.lingrid[:,ip] - center[:])^2
    end
    return Vpot
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

# What is included in the file `FMI2_fmu_to_md.jl` (FMU to model description)?
# - wrappers to call the model description functions from a FMU-instance [exported]

"""
Returns the number of states of the FMU.
"""
function fmi2GetNumberOfStates(fmu::FMU2)
    fmi2GetNumberOfStates(fmu.modelDescription)
end

function fmi2GetModelName(fmu::FMU2)
    fmi2GetModelName(fmu.modelDescription)
end

function fmi2GetGUID(fmu::FMU2)
    fmi2GetGUID(fmu.modelDescription)
end

function fmi2GetGenerationTool(fmu::FMU2)
    fmi2GetGenerationTool(fmu.modelDescription)
end

function fmi2GetGenerationDateAndTime(fmu::FMU2)
    fmi2GetGenerationDateAndTime(fmu.modelDescription)
end

function fmi2GetVariableNamingConvention(fmu::FMU2)
    fmi2GetVariableNamingConvention(fmu.modelDescription)
end

function fmi2GetNumberOfEventIndicators(fmu::FMU2)
    fmi2GetNumberOfEventIndicators(fmu.modelDescription)
end

function fmi2CanGetSetState(fmu::FMU2)
    fmi2CanGetSetState(fmu.modelDescription)
end

function fmi2CanSerializeFMUstate(fmu::FMU2)
    fmi2CanSerializeFMUstate(fmu.modelDescription)
end

function fmi2ProvidesDirectionalDerivative(fmu::FMU2)
    fmi2ProvidesDirectionalDerivative(fmu.modelDescription)
end

function fmi2IsCoSimulation(fmu::FMU2)
    fmi2IsCoSimulation(fmu.modelDescription)
end

function fmi2IsModelExchange(fmu::FMU2)
    fmi2IsModelExchange(fmu.modelDescription)
end


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

@doc raw"""CustomResourceConversion describes how to convert different versions of a CR.

    IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion(;
        conversionReviewVersions=nothing,
        strategy=nothing,
        webhookClientConfig=nothing,
    )

    - conversionReviewVersions::Vector{String} : conversionReviewVersions is an ordered list of preferred `ConversionReview` versions the Webhook expects. The API server will use the first version in the list which it supports. If none of the versions specified in this list are supported by API server, conversion will fail for the custom resource. If a persisted Webhook configuration specifies allowed versions and does not include any versions known to the API Server, calls to the webhook will fail. Defaults to `[\"v1beta1\"]`.
    - strategy::String : strategy specifies how custom resources are converted between versions. Allowed values are: - `None`: The converter only change the apiVersion and would not touch any other field in the custom resource. - `Webhook`: API Server will call to an external webhook to do the conversion. Additional information   is needed for this option. This requires spec.preserveUnknownFields to be false, and spec.conversion.webhookClientConfig to be set.
    - webhookClientConfig::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1WebhookClientConfig : webhookClientConfig is the instructions for how to call the webhook if strategy is `Webhook`. Required when `strategy` is set to `Webhook`.
"""
mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion <: SwaggerModel
    conversionReviewVersions::Any # spec type: Union{ Nothing, Vector{String} } # spec name: conversionReviewVersions
    strategy::Any # spec type: Union{ Nothing, String } # spec name: strategy
    webhookClientConfig::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1WebhookClientConfig } # spec name: webhookClientConfig

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion(;conversionReviewVersions=nothing, strategy=nothing, webhookClientConfig=nothing)
        o = new()
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion, Symbol("conversionReviewVersions"), conversionReviewVersions)
        setfield!(o, Symbol("conversionReviewVersions"), conversionReviewVersions)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion, Symbol("strategy"), strategy)
        setfield!(o, Symbol("strategy"), strategy)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion, Symbol("webhookClientConfig"), webhookClientConfig)
        setfield!(o, Symbol("webhookClientConfig"), webhookClientConfig)
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion = Dict{Symbol,Symbol}(Symbol("conversionReviewVersions")=>Symbol("conversionReviewVersions"), Symbol("strategy")=>Symbol("strategy"), Symbol("webhookClientConfig")=>Symbol("webhookClientConfig"))
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion = Dict{Symbol,String}(Symbol("conversionReviewVersions")=>"Vector{String}", Symbol("strategy")=>"String", Symbol("webhookClientConfig")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1WebhookClientConfig")
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion)
    (getproperty(o, Symbol("strategy")) === nothing) && (return false)
    true
end

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


text:
module SimulationBasedStability

using Reexport
@reexport using LinearAlgebra
#@reexport using SparseArrays
#@reexport using StaticArrays
#@reexport using Arpack
@reexport using DifferentialEquations
@reexport using Random
#using QuadGK
#using Lazy: iterated, take

include("structures.jl")
#include("structures_input.jl")
#include("structures_result.jl")

include("functions.jl")
#include("functions_discretization.jl")
#include("functions_method.jl")

export
dynamic_problem,
histopryremapp,
CompRand
#SemiDiscretization, NumericSD,
#ProportionalMX,
#Delay,DelayMX,
#Additive,
#LDDEProblem,
#DiscreteMapping, DiscreteMapping_1step,
#fixPointOfMapping, spectralRadiusOfMapping,
#DiscreteMapping_LR,spectralRadiusOfMappingLR

end # module


text:
"""
A CPD for which the distribution never changes.
    target: name of the CPD's variable
    parents: list of parent variables.
    d: a Distributions.jl distribution

While a StaticCPD can have parents, their assignments will not affect the distribution.
"""
mutable struct StaticCPD{D} <: CPD{D}
    target::NodeName
    parents::NodeNames
    d::D
end
StaticCPD(target::NodeName, d::Distribution) = StaticCPD(target, NodeName[], d)

name(cpd::StaticCPD) = cpd.target
parents(cpd::StaticCPD) = cpd.parents
(cpd::StaticCPD)(a::Assignment) = cpd.d # no update
(cpd::StaticCPD)() = (cpd)(Assignment()) # cpd()
(cpd::StaticCPD)(pair::Pair{NodeName}...) = (cpd)(Assignment(pair)) # cpd(:A=>1)

function Distributions.fit(::Type{StaticCPD{D}},
                           data::DataFrame,
                           target::NodeName,
                           parents::NodeNames=NodeName[],
                           ) where {D<:Distribution}

    d = fit(D, data[target])
    StaticCPD(target, parents, d)
end

nparams(cpd::StaticCPD) = paramcount(params(cpd.d))


text:
@enum Version begin

  NO_DEFAULT_OPEN_CLOSE    = 155
  PRICE_BASED_VOLATILITY   = 156
  REPLACE_FA_END           = 157
  DURATION                 = 158
  MARKET_DATA_IN_SHARES    = 159
  POST_TO_ATS              = 160
  WSHE_CALENDAR            = 161
  AUTO_CANCEL_PARENT       = 162
  FRACTIONAL_SIZE_SUPPORT  = 163

end


text:
module Canopy

include("constants.jl")
include("water.jl")
include("air.jl")
include("physchem/physchem.jl")
include("radtrans/radtrans.jl")
include("transfer/transfer.jl")
include("leaf/leaf.jl")
include("watercycle.jl")
include("units.jl")

end  # module


text:
using CUDArt, CUSPARSE, CUDNN, CUBLAS
#using Knet,
#import Knet: forw, back, loss, ninputs, overwrites, back_reads_x, back_reads_y, gpu, axpb!, @gpu, issimilar, mul2!
#using Base.LinAlg: axpy!, scale!

include("util/gpu.jl")
include("util/cudart.jl")
include("util/curand.jl")
include("util/cusparse.jl")
include("util/linalg.jl")
include("model.jl")
include("op.jl")
include("op/actf.jl");
include("op/add.jl");
include("op/dot.jl");
include("op/input.jl");
include("op/loss.jl");
include("op/mul.jl");
include("op/par.jl");
include("net.jl")
include("net/initforw.jl")
include("net/initback.jl")
include("net/forw.jl")
include("net/back.jl")
include("net/util.jl")
include("update.jl")

using Base.Test
include("data.jl")
# include("util/dense.jl")
include("util/colops.jl")
isdefined(:MNIST) || include("../examples/mnist.jl")
setseed(42)
nbatch=100
f64(x)=map(Float64,x)

dtrn = ItemTensor(f64(MNIST.xtrn), f64(MNIST.ytrn); batch=nbatch)
dtst = ItemTensor(f64(MNIST.xtst), f64(MNIST.ytst); batch=nbatch)

x0 = copy(dtrn.data[1])
y0 = copy(dtrn.data[2])

info("Testing simple mlp")

softmax() = quote
    x1 = input()
    x2 = soft(x1)
    x3 = softloss(x2)
end

layer(;n=1,f=nothing) = quote
    x1 = input()
    w1 = par($n,0)
    x2 = dot(w1,x1)
    b2 = par(0)
    x3 = add(b2,x2)
    y3 = $(symbol(f))(x3)
end

function mlp(loss, actf, hidden...)
    net = quote
        x0 = input()
    end
    N = length(hidden)
    for n=1:N
        x1 = symbol("x$(n-1)")
        x2 = symbol("x$n")
        push!(net.args, :($x2 = layer($x1; n=$(hidden[n]), f=$(n<N ? actf : loss))))
    end
    return net
end

net = Net(mlp(softmax, relu, 64, 10))
DBG = false

setopt!(net, lr=0.5)
@time for i=1:3
    @show (l,w,g) = train(net, dtrn; gclip=0, gcheck=1000, getloss=true, getnorm=true, atol=0.0001, rtol=0.0001)
    @show (test(net, dtrn), accuracy(net, dtrn))
    @show (test(net, dtst), accuracy(net, dtst))
end

@test isequal(x0,dtrn.data[1])
@test isequal(y0,dtrn.data[2])

# DBG = true
# nout = 2
# winit = rand(4,3)
# prog = quote
#     x = input()
#     w = par($nout, 0; init=Gaussian(1,3))
#     # w = par($winit)
#     y = dot(w,x)
#     r = relu(y)
#     z = qloss(r)
# end
# net1 = Net(prog)
# x = rand(3,5)
# forw(net1, x)

# function lstm2(; n=1)
#     quote
#         x  = input()
#         i  = lin2(x,h;n=$n,f=sigm)
#         f  = lin2(x,h;n=$n,f=sigm)
#         o  = lin2(x,h;n=$n,f=sigm)
#         g  = lin2(x,h;n=$n,f=tanh)
#         ig = mul(i,g)
#         fc = mul(f,c)
#         c  = add(ig,fc)
#         tc = tanh(c)
#         h  = mul(tc,o)
#     end
# end

# function lin2(; n=1,f=sigm)
#     quote
#         w1 = par($n,0)
#         x1 = input()
#         y1 = dot(w1,x1)
#         w2 = par($n,0)
#         x2 = input()
#         y2 = dot(w2,x2)
#         w3 = par(0)
#         x3 = add(y1,y2)
#         y3 = add(w3,x3)
#         # ou = symbol($f)(y3)
#         ou = relu(y3)
#     end
# end

# net = Net(lstm2(n=10))
# forw(net, x)
# forw(net, x)
# dy = convert(Array, net.out0[end])
# back(net, dy)

# s = sprand(3,5,.5)
# snet = Net(lstm2(n=10))
# forw(snet, s)
# forw(snet, s)
# sy = convert(SparseMatrixCSC, net.out0[end])
# # TODO: this is failing
# # back(snet, dy)

# # Op			out0	dif	size	flags

# # (1,Input)		7e00 	0	(3,5)	tosave,toincr,!toback
# # (2,Par)               8000 	2400	(10,3)
# # (3,Dot,2,1)           8200+ 	1c00+	(10,5)
# # (4,Par)               8400 	2000	(10,10)
# # (5,Dot,4,38)          8800*	1c00+	(10,5)
# # (6,Par)               8a00    1e00	(10,)
# # (7,Add,3,5)           8800*	1c00+	(10,5)
# # (8,Add,6,7)           8800*	1c00+	(10,5)
# # (9,Knet.Relu,8)	8c00    1c00+	(10,5)	tosave
# # (10,Par)              8e00    1a00	(10,3)
# # (11,Dot,10,1)         8800*	1200*	(10,5)
# # (12,Par)              9000    1600	(10,10)
# # (13,Dot,12,38)        8200+	1200*	(10,5)
# # (14,Par)              9400    1400	(10,)
# # (15,Add,11,13)        8200+	1200*	(10,5)
# # (16,Add,14,15)        8200+	1200*	(10,5)
# # (17,Knet.Relu,16)	9600    1200*	(10,5)	tosave
# # (18,Par)              9800    1000	(10,3)
# # (19,Dot,18,1)         8200+	0800!	(10,5)
# # (20,Par)              9a00    0c00	(10,10)
# # (21,Dot,20,38)        8800*	0800!	(10,5)
# # (22,Par)              9e00    0a00	(10,)
# # (23,Add,19,21)        8800*	0800!	(10,5)
# # (24,Add,22,23)        8800*	0800!	(10,5)
# # (25,Knet.Relu,24)	a000    0800!	(10,5)	tosave
# # (26,Par)	        a200    0600	(10,3)
# # (27,Dot,26,1)         8800*	fe00@	(10,5)
# # (28,Par)              a400    0200	(10,10)
# # (29,Dot,28,38)        8200+	fe00@	(10,5)
# # (30,Par)              a800    0000	(10,)
# # (31,Add,27,29)        8200+	fe00@	(10,5)
# # (32,Add,30,31)        8200+	fe00@	(10,5)
# # (33,Knet.Relu,32)	aa00    fe00@	(10,5)	tosave
# # (34,Mul,9,33)         8200+	fc00	(10,5)
# # (35,Mul,17,36)        8800*	fa00	(10,5)
# # (36,Add,34,35)	ac00    f800	(10,5)	tosave,toincr,tozero,tmp=f400
# # (37,Knet.Tanh,36)	ae00    f600	(10,5)	tosave
# # (38,Mul,37,25)	b000	f200	(10,5)	tosave,toincr,tozero,tmp=f400

# # ops: 38

# # out0: 22 unique registers, 12 par, 8 tosave, 2 temp

# # dif: 20 unique registers, 12 par, 2 toincr

# # tmp: where is tmp for par? only needed if sequence.


text:
INITIALIZE_DEFAULT(cb,t,u,integrator) = nothing

immutable ContinuousCallback{F1,F2,F3,F4,T,T2,I} <: DECallback
  condition::F1
  affect!::F2
  affect_neg!::F3
  initialize::F4
  idxs::I
  rootfind::Bool
  interp_points::Int
  save_positions::Tuple{Bool,Bool}
  abstol::T
  reltol::T2
end

ContinuousCallback(condition,affect!,affect_neg!;
                   initialize = INITIALIZE_DEFAULT,
                   idxs = nothing,
                   rootfind=true,
                   save_positions=(true,true),
                   interp_points=10,
                   abstol=1e-12,reltol=0) = ContinuousCallback(
                              condition,affect!,affect_neg!,initialize,
                              idxs,
                              rootfind,interp_points,
                              save_positions,abstol,reltol)

function ContinuousCallback(condition,affect!;
                   initialize = INITIALIZE_DEFAULT,
                   idxs = nothing,
                   rootfind=true,
                   save_positions=(true,true),
                   affect_neg! = affect!,
                   interp_points=10,
                   abstol=1e-12,reltol=0)

 ContinuousCallback(
            condition,affect!,affect_neg!,initialize,idxs,
            rootfind,interp_points,
            save_positions,abstol,reltol)

end

immutable DiscreteCallback{F1,F2,F3} <: DECallback
  condition::F1
  affect!::F2
  initialize::F3
  save_positions::Tuple{Bool,Bool}
end
DiscreteCallback(condition,affect!;
        initialize = INITIALIZE_DEFAULT,save_positions=(true,true)) = DiscreteCallback(condition,affect!,initialize,save_positions)

# DiscreteCallback(condition,affect!,save_positions) = DiscreteCallback(condition,affect!,save_positions)

immutable CallbackSet{T1<:Tuple,T2<:Tuple} <: DECallback
  continuous_callbacks::T1
  discrete_callbacks::T2
end

CallbackSet(callback::DiscreteCallback) = CallbackSet((),(callback,))
CallbackSet(callback::ContinuousCallback) = CallbackSet((callback,),())
CallbackSet() = CallbackSet((),())
CallbackSet(cb::Void) = CallbackSet()

# For Varargs, use recursion to make it type-stable
CallbackSet(callbacks::Union{DECallback,Void}...) = CallbackSet(split_callbacks((), (), callbacks...)...)

@inline split_callbacks(cs, ds) = cs, ds
@inline split_callbacks(cs, ds, c::Void, args...) = split_callbacks(cs, ds, args...)
@inline split_callbacks(cs, ds, c::ContinuousCallback, args...) = split_callbacks((cs..., c), ds, args...)
@inline split_callbacks(cs, ds, d::DiscreteCallback, args...) = split_callbacks(cs, (ds..., d), args...)
@inline function split_callbacks(cs, ds, d::CallbackSet, args...)
  split_callbacks((cs...,d.continuous_callbacks...), (ds..., d.discrete_callbacks...), args...)
end

initialize!(cb::CallbackSet,t,u,integrator::DEIntegrator) = initialize!(t,u,integrator,cb.continuous_callbacks...,cb.discrete_callbacks...)
initialize!(cb::CallbackSet{Tuple{},Tuple{}},t,u,integrator::DEIntegrator) = nothing
initialize!(t,u,integrator::DEIntegrator,c::DECallback,cs::DECallback...) = (c.initialize(c,t,u,integrator); initialize!(t,u,integrator,cs...))
initialize!(t,u,integrator::DEIntegrator,c::DECallback) = c.initialize(c,t,u,integrator)


text:
using CoDa
using Test

# list of maintainers
maintainers = ["juliohm"]

# environment settings
istravis = "TRAVIS" ∈ keys(ENV)
ismaintainer = "USER" ∈ keys(ENV) && ENV["USER"] ∈ maintainers
datadir = joinpath(@__DIR__,"data")

# list of tests
testfiles = [
  "compositions.jl",
  "transforms.jl"
]

@testset "CoDa.jl" begin
  for testfile in testfiles
    include(testfile)
  end
end


text:
# Tests for conditional entropy function

entropy1 = get_entropy(arr1)
entropy2 = get_entropy(arr2)
entropy12 = get_entropy(arr1, arr2)

# Test with 1 values array conditioned on itself
@test get_conditional_entropy(arr1, arr1) ≈ 0 atol = 0.05
println("Conditional entropy with 1 array passed.")

conditional_entropy = entropy12 - entropy2

# Test conditional entropy
@test get_conditional_entropy(arr1, arr2) ≈ conditional_entropy atol = 0.05
println("Conditional entropy passed.")

# Test conditional entropy formula
@test apply_conditional_entropy_formula(entropy12, entropy2) ≈ conditional_entropy atol = 0.05
println("Conditional entropy formula passed.")


text:
# Test for weights with a small set of trials

using SubmodularMaximization
using PyPlot
using Base.Iterators
using Base.Threads
using JLD2
using Statistics
using HistogramFilters
using RosDataProcess

close("all")

experiment_name = "large_scale_weights"
data_folder = "./data"

reprocess = false

num_robots = 8:8:100

trials = 1:5
steps = 100
trial_steps = 20:steps

# Range limits
communication_range = 20
robot_target_range_limit = 12
weights_range_limit = 20 # Use a larger range for weights for accurate results

solver = x->solve_communication_range_limit(x, num_partitions=4,
                                            communication_range=
                                            communication_range)

# Plotting configs

line_styles = ["-", "--"]
experiment_names = ["large-scale", "solver-trials"]

# Use colors from the default cycle
colors = get_cmap("tab10")
weights_color = colors(0)
objective_color = colors(1)

#
# Code to run experiments
#

# We will also load a separate set of results from the experiments
results_array = []

all_tests = product(num_robots, trials)
all_configurations = num_robots

function trial_fun(x)
  num_robots, trial = x

  num_targets = default_num_targets(num_robots)

  configs = MultiRobotTargetTrackingConfigs(num_robots,
                                            robot_target_range_limit=
                                              robot_target_range_limit
                                           )

  full_trial_data = target_tracking_experiment(steps=steps,
                                               num_robots=num_robots,
                                               configs=configs,
                                               solver=solver)

  weights = map(full_trial_data[trial_steps]) do step_data
        robot_states = step_data.robot_states
        target_filters = step_data.target_filters

        problem = MultiRobotTargetTrackingProblem(robot_states, target_filters,
                                                  configs)

        compute_weight_matrix(problem,
                              robot_target_range_limit=weights_range_limit
                             )
  end

  (trial_weights=weights,
   trial_data=full_trial_data[trial_steps],
   full_trial_data=full_trial_data,
   configs=configs)
end
function print_summary(x)
  (num_robots, trial) = x
  println("Num. Robots: ", num_robots,
          " Trial: ", trial)
end

let
  @time results = run_experiments(all_tests,
                                  trial_fun=trial_fun,
                                  print_summary=print_summary,
                                  experiment_name=experiment_name,
                                  data_folder=data_folder,
                                  reprocess=reprocess,
                                  threaded=true
                                 )
  push!(results_array, (results, num_robots, trials))
end

#
# *Optionally* load the weights dataset from the full set of trials if the data
# exists
#

experiment_weights_file = string(data_folder,
                                 "/weights_by_number_of_robots.jld2")

let
  if isfile(experiment_weights_file)
    @load experiment_weights_file results

    num_robots = sort(union([first(x) for x in keys(results)]))
    trials = sort(union([last(x) for x in keys(results)]))

    push!(results_array, (results, num_robots, trials))
  end
end

#
# Weights per num robots
#

figure()

for ii in 1:length(results_array)
  (results, num_robots, trials) = results_array[ii]
  style = line_styles[ii]
  experiment_name = experiment_names[ii]

  weight_per_robot = map(num_robots) do num_robots
    trial_weights = map(trials) do trial
      trial = results[num_robots, trial]

      trial_weights = trial.trial_weights

      total_weights = map(total_weight, trial_weights)

      total_weights / num_robots
    end

    vcat(trial_weights...)
  end

  weights_series = TimeSeries(num_robots, hcat(weight_per_robot...)')
  plot_trials(weights_series, mean=true, marker=".", markersize=25, linewidth=3,
              color=weights_color, linestyle=style,
              label=string("redundancy, ", experiment_name))
end

tick_params(axis=:y, labelcolor=weights_color)
ylabel("Redundancy per robot", color=weights_color)

weights_axis = gca()

#
# Plot objective values
#

twinx()

for ii in 1:length(results_array)
  (results, num_robots, trials) = results_array[ii]
  style = line_styles[ii]
  experiment_name = experiment_names[ii]

  objective_values = map(num_robots) do num_robots
    trial_objectives = map(trials) do trial
      trial = results[num_robots, trial]

      objectives = map(x->x.objective, trial.trial_data)

      objectives / num_robots
    end

    vcat(trial_objectives...)
  end

  objective_series = TimeSeries(num_robots, hcat(objective_values...)')
  plot_trials(objective_series, mean=true, marker=".", markersize=25,
              linewidth=3, color=objective_color, linestyle=style,
              label=string("objective, ", experiment_name))
end

tick_params(axis=:y, labelcolor=objective_color)
ylabel("Objective value per robot", color=objective_color)

objective_axis = gca()

lines1, labels1 = weights_axis.get_legend_handles_labels()
lines2, labels2 = objective_axis.get_legend_handles_labels()

legend(vcat(lines1, lines2), vcat(labels1, labels2), loc="lower right")

xlabel("Number of robots")
title("Weights and objective values per num. robots")

save_fig("fig", "large_scale_weights_per_num_robots")

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing