text:
# VGG16 model for convolutional feature extraction
function vgg16(w, x; o=Dict())
    # get parameters
    fc6drop = get(o, :fc6drop, 0.0)
    fc7drop = get(o, :fc7drop, 0.0)
    mode = get(o, :cnnmode, 1)
    featuremaps = get(o, :featuremaps, false)

    conv1_1 = conv4(w[1], x; padding=1, mode=mode) .+ w[2]
    conv1_1 = relu.(conv1_1)
    conv1_2 = conv4(w[3], conv1_1; padding=1, mode=mode) .+ w[4]
    conv1_2 = relu.(conv1_2)
    pool1   = pool(conv1_2)

    conv2_1 = conv4(w[5], pool1; padding=1, mode=mode) .+ w[6]
    conv2_1 = relu.(conv2_1)
    conv2_2 = conv4(w[7], conv2_1; padding=1, mode=mode) .+ w[8]
    conv2_2 = relu.(conv2_2)
    pool2   = pool(conv2_2)

    conv3_1 = conv4(w[9], pool2; padding=1, mode=mode) .+ w[10]
    conv3_1 = relu.(conv3_1)
    conv3_2 = conv4(w[11], conv3_1; padding=1, mode=mode) .+ w[12]
    conv3_2 = relu.(conv3_2)
    conv3_3 = conv4(w[13], conv3_2; padding=1, mode=mode) .+ w[14]
    conv3_3 = relu.(conv3_3)
    pool3   = pool(conv3_3)

    conv4_1 = conv4(w[15], pool3; padding=1, mode=mode) .+ w[16]
    conv4_1 = relu.(conv4_1)
    conv4_2 = conv4(w[17], conv4_1; padding=1, mode=mode) .+ w[18]
    conv4_2 = relu.(conv4_2)
    conv4_3 = conv4(w[19], conv4_2; padding=1, mode=mode) .+ w[20]
    conv4_3 = relu.(conv4_3)
    pool4   = pool(conv4_3)

    conv5_1 = conv4(w[21], pool4; padding=1, mode=mode) .+ w[22]
    conv5_1 = relu.(conv5_1)
    conv5_2 = conv4(w[23], conv5_1; padding=1, mode=mode) .+ w[24]
    conv5_2 = relu.(conv5_2)
    conv5_3 = conv4(w[25], conv5_2; padding=1, mode=mode) .+ w[26]
    conv5_3 = relu.(conv5_3)

    if !featuremaps
        pool5 = pool(conv5_3)
        fc6 = w[27] * mat(pool5) .+ w[28]
        fc6 = relu.(fc6)
        fc6 = dropout(fc6, fc6drop)
        fc7 = w[29] * mat(fc6) .+ w[30]
        fc7 = relu.(fc7)
    end
end

# VGG16 model for convolutional feature extraction
function vgg19(w, x; o=Dict())
    # get parameters
    fc6drop = get(o, :fc6drop, 0.0)
    fc7drop = get(o, :fc7drop, 0.0)
    mode = get(o, :cnnmode, 1)
    featuremaps = get(o, :featuremaps, false)

    conv1_1 = conv4(w[1], x; padding=1, mode=mode) .+ w[2]
    conv1_1 = relu.(conv1_1)
    conv1_2 = conv4(w[3], conv1_1; padding=1, mode=mode) .+ w[4]
    conv1_2 = relu.(conv1_2)
    pool1   = pool(conv1_2)

    conv2_1 = conv4(w[5], pool1; padding=1, mode=mode) .+ w[6]
    conv2_1 = relu.(conv2_1)
    conv2_2 = conv4(w[7], conv2_1; padding=1, mode=mode) .+ w[8]
    conv2_2 = relu.(conv2_2)
    pool2   = pool(conv2_2)

    conv3_1 = conv4(w[9], pool2; padding=1, mode=mode) .+ w[10]
    conv3_1 = relu.(conv3_1)
    conv3_2 = conv4(w[11], conv3_1; padding=1, mode=mode) .+ w[12]
    conv3_2 = relu.(conv3_2)
    conv3_3 = conv4(w[13], conv3_2; padding=1, mode=mode) .+ w[14]
    conv3_3 = relu.(conv3_3)
    conv3_4 = conv4(w[15], conv3_3; padding=1, mode=mode) .+ w[16]
    conv3_4 = relu.(conv3_4)
    pool3   = pool(conv3_4)

    conv4_1 = conv4(w[17], pool3; padding=1, mode=mode) .+ w[18]
    conv4_1 = relu.(conv4_1)
    conv4_2 = conv4(w[19], conv4_1; padding=1, mode=mode) .+ w[20]
    conv4_2 = relu.(conv4_2)
    conv4_3 = conv4(w[21], conv4_2; padding=1, mode=mode) .+ w[22]
    conv4_3 = relu.(conv4_3)
    conv4_4 = conv4(w[23], conv4_3; padding=1, mode=mode) .+ w[24]
    conv4_4 = relu.(conv4_4)
    pool4   = pool(conv4_4)

    conv5_1 = conv4(w[25], pool4; padding=1, mode=mode) .+ w[26]
    conv5_1 = relu.(conv5_1)
    conv5_2 = conv4(w[27], conv5_1; padding=1, mode=mode) .+ w[28]
    conv5_2 = relu.(conv5_2)
    conv5_3 = conv4(w[29], conv5_2; padding=1, mode=mode) .+ w[30]
    conv5_3 = relu.(conv5_3)
    conv5_4 = conv4(w[31], conv5_3; padding=1, mode=mode) .+ w[32]
    conv5_4 = relu.(conv5_4)

    if !featuremaps
        pool5 = pool(conv5_4)
        fc6 = w[33] * mat(pool5) .+ w[34]
        fc6 = relu.(fc6)
        fc6 = dropout(fc6, fc6drop)
        fc7 = w[35] * mat(fc6) .+ w[36]
        fc7 = relu.(fc7)
        return fc7
    end

    return conv5_4
end

function get_vgg_weights(vggmat; last_layer="relu.7")
    ws = []
    for l in vggmat["layers"]
        if haskey(l, "weights") && length(l["weights"]) != 0
            w = l["weights"]
            if startswith(l["name"], "conv")
                w[2] = reshape(w[2], (1,1,length(w[2]),1))
            elseif startswith(l["name"], "fc")
                w[1] = mat(w[1])'
            end
            push!(ws, w...)
        end

        last_layer == l["name"] && break
    end
    map(KnetArray, ws)
end

# mode, 0=>train, 1=>test
function resnet50(w,x,ms; mode=1)
    # layer 1
    conv1  = conv4(w[1],x; padding=3, stride=2) .+ w[2]
    bn1    = batchnorm(w[3:4],conv1,ms; mode=mode)
    pool1  = pool(bn1; window=3, stride=2)

    # layer 2,3,4,5
    r2 = reslayerx5(w[5:34], pool1, ms; strides=[1,1,1,1], mode=mode)
    r3 = reslayerx5(w[35:73], r2, ms; mode=mode)
    r4 = reslayerx5(w[74:130], r3, ms; mode=mode) # 5
    r5 = reslayerx5(w[131:160], r4, ms; mode=mode)

    # fully connected layer
    pool5  = pool(r5; stride=1, window=7, mode=2)
    fc1000 = w[161] * mat(pool5) .+ w[162]
end

# mode, 0=>train, 1=>test
function resnet101(w,x,ms; mode=1)
    # layer 1
    conv1 = reslayerx1(w[1:3],x,ms; padding=3, stride=2, mode=mode)
    pool1 = pool(conv1; window=3, stride=2)

    # layer 2,3,4,5
    r2 = reslayerx5(w[4:33], pool1, ms; strides=[1,1,1,1], mode=mode)
    r3 = reslayerx5(w[34:72], r2, ms; mode=mode)
    r4 = reslayerx5(w[73:282], r3, ms; mode=mode)
    r5 = reslayerx5(w[283:312], r4, ms; mode=mode)

    # fully connected layer
    pool5  = pool(r5; stride=1, window=7, mode=2)
    fc1000 = w[313] * mat(pool5) .+ w[314]
end

# mode, 0=>train, 1=>test
function resnet152(w,x,ms; mode=1)
    # layer 1
    conv1 = reslayerx1(w[1:3],x,ms; padding=3, stride=2, mode=mode)
    pool1 = pool(conv1; window=3, stride=2)

    # layer 2,3,4,5
    r2 = reslayerx5(w[4:33], pool1, ms; strides=[1,1,1,1], mode=mode)
    r3 = reslayerx5(w[34:108], r2, ms; mode=mode)
    r4 = reslayerx5(w[109:435], r3, ms; mode=mode)
    r5 = reslayerx5(w[436:465], r4, ms; mode=mode)

    # fully connected layer
    pool5  = pool(r5; stride=1, window=7, mode=2)
    fc1000 = w[466] * mat(pool5) .+ w[467]
end

# Batch Normalization Layer
# works both for convolutional and fully connected layers
# mode, 0=>train, 1=>test
function batchnorm(w, x, ms; mode=1, epsilon=1e-5)
    mu, sigma = nothing, nothing
    if mode == 0
        d = ndims(x) == 4 ? (1,2,4) : (2,)
        s = reduce((a,b)->a*size(x,b), d)
        mu = sum(x,d) / s
        sigma = sqrt(epsilon + (sum(x.-mu,d).^2) / s)
    elseif mode == 1
        mu = shift!(ms)
        sigma = shift!(ms)
    end

    # we need getval in backpropagation
    push!(ms, AutoGrad.getval(mu), AutoGrad.getval(sigma))
    xhat = (x.-mu) ./ sigma
    return w[1] .* xhat .+ w[2]
end

function reslayerx0(w,x,ms; padding=0, stride=1, mode=1)
    b  = conv4(w[1],x; padding=padding, stride=stride)
    bx = batchnorm(w[2:3],b,ms; mode=mode)
end

function reslayerx1(w,x,ms; padding=0, stride=1, mode=1)
    relu.(reslayerx0(w,x,ms; padding=padding, stride=stride, mode=mode))
end

function reslayerx2(w,x,ms; pads=[0,1,0], strides=[1,1,1], mode=1)
    ba = reslayerx1(w[1:3],x,ms; padding=pads[1], stride=strides[1], mode=mode)
    bb = reslayerx1(w[4:6],ba,ms; padding=pads[2], stride=strides[2], mode=mode)
    bc = reslayerx0(w[7:9],bb,ms; padding=pads[3], stride=strides[3], mode=mode)
end

function reslayerx3(w,x,ms; pads=[0,0,1,0], strides=[2,2,1,1], mode=1) # 12
    a = reslayerx0(w[1:3],x,ms; stride=strides[1], padding=pads[1], mode=mode)
    b = reslayerx2(w[4:12],x,ms; strides=strides[2:4], pads=pads[2:4], mode=mode)
    relu.(a .+ b)
end

function reslayerx4(w,x,ms; pads=[0,1,0], strides=[1,1,1], mode=1)
    relu.(x .+ reslayerx2(w,x,ms; pads=pads, strides=strides, mode=mode))
end

function reslayerx5(w,x,ms; strides=[2,2,1,1], mode=1)
    x = reslayerx3(w[1:12],x,ms; strides=strides, mode=mode)
    for k = 13:9:length(w)
        x = reslayerx4(w[k:k+8],x,ms; mode=mode)
    end
    return x
end

function get_params(params)
    len = length(params["value"])
    ws, ms = [], []
    for k = 1:len
        name = params["name"][k]
        value = convert(Array{Float32}, params["value"][k])

        if endswith(name, "moments")
            push!(ms, reshape(value[:,1], (1,1,size(value,1),1)))
            push!(ms, reshape(value[:,2], (1,1,size(value,1),1)))
        elseif startswith(name, "bn")
            push!(ws, reshape(value, (1,1,length(value),1)))
        elseif startswith(name, "fc") && endswith(name, "filter")
            push!(ws, transpose(reshape(value,size(value,3,4))))
        elseif startswith(name, "conv") && endswith(name, "bias")
            push!(ws, reshape(value, (1,1,length(value),1)))
        else
            push!(ws, value)
        end
    end
    map(KnetArray, ws), map(KnetArray, ms)
end

max_stars_repo_name: ilkerkesen/kusrp

max_issues_repo_name: ilkerkesen/kusrp

max_forks_repo_name: ilkerkesen/kusrp


text:
using Pkg

Pkg.activate(normpath(joinpath(@__DIR__, "..")))

push!(LOAD_PATH, "$(pwd())/Maple/src") # Don't use this when this is an actuall module itself

TEST_RUNNING = true
include("$(pwd())/src/Ahorn.jl")
Ahorn.displayMainWindow()

ext: jl


text:
default(::SpecHandle) = :C_NULL
function default(spec::Union{SpecStructMember,SpecFuncParam})
    is_length_exception(spec) && return default_length_exception(spec)
    @match spec.requirement begin
        if spec.type ∈ spec_handles.name
        end => default(handle_by_name(spec.type))
        &POINTER_OPTIONAL || &POINTER_REQUIRED || if is_ptr(spec.type) || spec.type == :Cstring
        end => :C_NULL
        &OPTIONAL || &REQUIRED => 0
    end
end

function default_length_exception(spec::Spec)
    @match parent(spec) begin
        :VkWriteDescriptorSet => :(max($((:(pointer_length($(wrap_identifier(arg)))) for arg in arglen(spec))...)))
    end
end


text:
if VERSION <= v"1.2"
  # Just return the expression (i.e. eval in single thread)
  macro spawn(expr)
    expr
  end
  macro spawn(expr, rest...)
    expr
  end
else
  using Base.Threads: @spawn
end


text:
in_channel = 3
out_channel = 5
N = 6

sg = SimpleGraph(N)
add_edge!(sg, 1, 2); add_edge!(sg, 1, 3); add_edge!(sg, 2, 3)
add_edge!(sg, 3, 4); add_edge!(sg, 2, 5); add_edge!(sg, 3, 6)

ug = MetaGraph(sg)
dg = MetaDiGraph(sg)

@testset "metagraphs" begin
    @testset "GCNConv" begin
        gc = GCNConv(ug, in_channel=>out_channel)
        @test size(gc.weight) == (out_channel, in_channel)
        @test size(gc.bias) == (out_channel,)
        @test graph(gc.fg) === ug.graph

        gc = GCNConv(dg, in_channel=>out_channel)
        @test size(gc.weight) == (out_channel, in_channel)
        @test size(gc.bias) == (out_channel,)
        @test graph(gc.fg) === dg.graph
    end

    @testset "ChebConv" begin
        k = 4
        cc = ChebConv(ug, in_channel=>out_channel, k)
        @test size(cc.weight) == (out_channel, in_channel, k)
        @test size(cc.bias) == (out_channel,)
        @test graph(cc.fg) == ug.graph
        @test cc.k == k
        @test cc.in_channel == in_channel
        @test cc.out_channel == out_channel

        cc = ChebConv(dg, in_channel=>out_channel, k)
        @test size(cc.weight) == (out_channel, in_channel, k)
        @test size(cc.bias) == (out_channel,)
        @test graph(cc.fg) == dg.graph
        @test cc.k == k
        @test cc.in_channel == in_channel
        @test cc.out_channel == out_channel
    end

    @testset "GraphConv" begin
        gc = GraphConv(ug, in_channel=>out_channel)
        @test graph(gc.fg) == ug.graph
        @test size(gc.weight1) == (out_channel, in_channel)
        @test size(gc.weight2) == (out_channel, in_channel)
        @test size(gc.bias) == (out_channel,)

        gc = GraphConv(dg, in_channel=>out_channel)
        @test graph(gc.fg) == dg.graph
        @test size(gc.weight1) == (out_channel, in_channel)
        @test size(gc.weight2) == (out_channel, in_channel)
        @test size(gc.bias) == (out_channel,)
    end

    @testset "GATConv" begin
        gat = GATConv(ug, in_channel=>out_channel)
        @test graph(gat.fg) == ug.graph
        @test size(gat.weight) == (out_channel, in_channel)
        @test size(gat.bias) == (out_channel,)

        gat = GATConv(dg, in_channel=>out_channel)
        @test graph(gat.fg) == dg.graph
        @test size(gat.weight) == (out_channel, in_channel)
        @test size(gat.bias) == (out_channel,)
    end

    @testset "GatedGraphConv" begin
        num_layers = 3
        ggc = GatedGraphConv(ug, out_channel, num_layers)
        @test graph(ggc.fg) == ug.graph
        @test size(ggc.weight) == (out_channel, out_channel, num_layers)

        ggc = GatedGraphConv(dg, out_channel, num_layers)
        @test graph(ggc.fg) == dg.graph
        @test size(ggc.weight) == (out_channel, out_channel, num_layers)
    end

    @testset "EdgeConv" begin
        ec = EdgeConv(ug, Dense(2*in_channel, out_channel))
        @test graph(ec.fg) == ug.graph

        ec = EdgeConv(dg, Dense(2*in_channel, out_channel))
        @test graph(ec.fg) == dg.graph
    end
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 20e7b5fd-0414-4998-935b-84708e1ea89e
using Plots, PlutoUI, SampledSignals

# ╔═╡ 42a53dad-e443-4292-8f35-18f7c85631e4
include("definitions/SoundPlots.jl")

# ╔═╡ 45c3da18-ba84-11eb-12d8-b30f345e6ddd
md"# Array synthesis
The experiment in the notebook PeriodicPitch.jl suggests
a simple way to synthesize complicated sounds at a given pitch.
First generate an array with the waveform of the required length,
then repeat it to get a note.
This seems especially attractive for waveforms that are expensive to compute,
since it is only computed once.
Also the repetitions perfectly matching the framerate give sounds without aliasing,
or at least only aliasing that overlaps with the harmonics,
as seen in the white noise examples."

# ╔═╡ ebfb2aac-3d80-497a-8b58-4e0c235439b3
fs = 48000

# ╔═╡ 13c78f70-2b94-4330-a819-f286d0f75f20
array = 2*Array(1:400)/400 .- 1

# ╔═╡ f67e7cb4-6f7a-401e-897e-d64adb029c2b
sawtooth = SampleBuf(repeat(array, fs÷200), fs)

# ╔═╡ e8a156ce-735e-4946-a243-9970f3640756
plotfreqs(sawtooth, fs; title="Sawtooth wave by repetition")

# ╔═╡ 8e9daa59-4e89-41ca-8a80-e357f0b67c5e
md"There are two major downsides to the technique:
- The waveform needs to be recalculated for each frequency
- The waveform can only be made at periods of $n/\text{fs}$ seconds,
  which means its frequencies can only be $\text{fs}/n$ Hz where $n$ is an integer"

# ╔═╡ 67bb4c8c-8ab2-4dc7-a43a-79b00dd07657
md"## Using Lookup Tables
Instead of just repeating the waveform we can use it as a lookup table.
Then we can change how much we increment the lookup index
to change the frequency of the note."

# ╔═╡ 45bbd3c6-2a2a-45d3-a4da-8ff1fd25d118
function sawlookup(t)
   return array[t % 400 + 1]
end

# ╔═╡ 42a5616a-ff89-4b55-8ec3-d025b01c30c4
@bind step Slider(1:10)

# ╔═╡ 9cfc65bb-c195-49da-8a11-176b54124224
plot(sawlookup.(step*(1:400)))

# ╔═╡ 36ae8b23-e946-469a-bba9-142e45421517
lookupsaw = SampleBuf(Array(step*sawlookup.(step*(1:fs))), fs)

# ╔═╡ 9446699f-6d9f-4d7d-9096-9f4e75d52aee
plotfreqs(lookupsaw, fs; title="Lookup table sawtooth wave")

# ╔═╡ bef54180-d713-4a32-a526-10a3aa5d8cfa
md"Sometimes there's a very regular pattern of aliasing and other times there isn't.
To understand this consider the pattern created by the increment.
If the increment is a divisor of the array's length it will loop around only part of the array.
If the increment is coprime it won't repeat until all array elements have been touched, and so the fundamental period is the period of the array and there can be aliasing at any of the original harmonics: This happens at increment 3 for an array of length 400. At any non-coprime non-factor a mix of the two happen and there is aliasing at a subset of the harmonics.

Lookup tables allow for frequencies of $\text{fs}/n$ and $\text{fs}/d$ where $d$ is a divisor of $n$ to be played without aliasing using the same array. Knowing this we can use array lengths that have lots of factors to get a selection of pitches from a single array. We are still limited by only getting fractions of the samplerate as freqencies and not being able to change pitch smoothly.
To overcome these limitations more sophisticated methods must be used."

# ╔═╡ Cell order:
# ╟─45c3da18-ba84-11eb-12d8-b30f345e6ddd
# ╠═20e7b5fd-0414-4998-935b-84708e1ea89e
# ╠═42a53dad-e443-4292-8f35-18f7c85631e4
# ╠═ebfb2aac-3d80-497a-8b58-4e0c235439b3
# ╠═13c78f70-2b94-4330-a819-f286d0f75f20
# ╠═f67e7cb4-6f7a-401e-897e-d64adb029c2b
# ╟─e8a156ce-735e-4946-a243-9970f3640756
# ╟─8e9daa59-4e89-41ca-8a80-e357f0b67c5e
# ╟─67bb4c8c-8ab2-4dc7-a43a-79b00dd07657
# ╠═45bbd3c6-2a2a-45d3-a4da-8ff1fd25d118
# ╠═42a5616a-ff89-4b55-8ec3-d025b01c30c4
# ╠═9cfc65bb-c195-49da-8a11-176b54124224
# ╠═36ae8b23-e946-469a-bba9-142e45421517
# ╠═9446699f-6d9f-4d7d-9096-9f4e75d52aee
# ╟─bef54180-d713-4a32-a526-10a3aa5d8cfa

max_stars_repo_name: jnha/music

max_issues_repo_name: jnha/music

max_forks_repo_name: jnha/music


text:
# In this file, we figure out how to link to Python (surprisingly complicated)
# and generate a deps/deps.jl file with the libpython name and other information
# needed for static compilation of PyCall.

# As a result, if you switch to a different version or path of Python, you
# will probably need to re-run Pkg.build("PyCall").

using VersionParsing
import Conda, Libdl

struct UseCondaPython <: Exception end

include("buildutils.jl")
include("depsutils.jl")

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

# we write configuration files only if they change, both
# to prevent unnecessary recompilation and to minimize
# problems in the unlikely event of read-only directories.
function writeifchanged(filename, str)
    if !isfile(filename) || read(filename, String) != str
        @info string(abspath(filename), " has been updated")
        write(filename, str)
    else
        @info string(abspath(filename), " has not changed")
    end
end

# return the first arg that exists in the PATH
function whichfirst(args...)
    for x in args
        if Sys.which(x) !== nothing
            return x
        end
    end
    return ""
end

const prefsfile = joinpath(first(DEPOT_PATH), "prefs", "PyCall")
mkpath(dirname(prefsfile))

try # make sure deps.jl file is removed on error
    python = try
        let py = get(ENV, "PYTHON", isfile(prefsfile) ? readchomp(prefsfile) :
                     (Sys.isunix() && !Sys.isapple()) || Sys.ARCH ∉ (:i686, :x86_64) ?
                     whichfirst("python3", "python") : "Conda"),
            vers = isempty(py) || py == "Conda" ? v"0.0" : vparse(pyconfigvar(py,"VERSION","0.0"))
            if vers < v"2.7"
                if isempty(py) || py == "Conda"
                    throw(UseCondaPython())
                else
                    error("Python version $vers < 2.7 is not supported")
                end
            end

            # check word size of Python via sys.maxsize, since a common error
            # on Windows is to link a 64-bit Julia to a 32-bit Python.
            pywordsize = parse(UInt64, pysys(py, "maxsize")) > (UInt64(1)<<32) ? 64 : 32
            if pywordsize != Sys.WORD_SIZE
                error("$py is $(pywordsize)-bit, but Julia is $(Sys.WORD_SIZE)-bit")
            end

            py
        end
    catch e1
        if Sys.ARCH in (:i686, :x86_64)
            if isa(e1, UseCondaPython)
                @info string("Using the Python distribution in the Conda package by default.\n",
                     "To use a different Python version, set ENV[\"PYTHON\"]=\"pythoncommand\" and re-run Pkg.build(\"PyCall\").")
            else
                @info string( "No system-wide Python was found; got the following error:\n",
                      "$e1\nusing the Python distribution in the Conda package")
            end
            abspath(Conda.PYTHONDIR, "python" * (Sys.iswindows() ? ".exe" : ""))
        else
            error("No system-wide Python was found; got the following error:\n",
                  "$e1")
        end
    end

    use_conda = dirname(python) == abspath(Conda.PYTHONDIR)
    if use_conda
        Conda.add("numpy")
    end

    (libpython, libpy_name) = find_libpython(python)
    programname = pysys(python, "executable")

    # Get PYTHONHOME, either from the environment or from Python
    # itself (if it is not in the environment or if we are using Conda)
    PYTHONHOME = if !haskey(ENV, "PYTHONHOME") || use_conda
        pythonhome_of(python)
    else
        ENV["PYTHONHOME"]
    end

    # cache the Python version as a Julia VersionNumber
    pyversion = vparse(pyvar(python, "platform", "python_version()"))

    @info "PyCall is using $python (Python $pyversion) at $programname, libpython = $libpy_name"

    if pyversion < v"2.7"
        error("Python 2.7 or later is required for PyCall")
    end

    writeifchanged("deps.jl", """
    const python = "$(escape_string(python))"
    const libpython = "$(escape_string(libpy_name))"
    const pyprogramname = "$(escape_string(programname))"
    const pyversion_build = $(repr(pyversion))
    const PYTHONHOME = ""

    "True if we are using the Python distribution in the Conda package."
    const conda = $use_conda
    """)

    # Make subsequent builds (e.g. Pkg.update) use the same Python by default:
    writeifchanged(prefsfile, use_conda ? "Conda" : isfile(programname) ? programname : python)

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

catch

    # remove deps.jl (if it exists) on an error, so that PyCall will
    # not load until it is properly configured.
    isfile("deps.jl") && rm("deps.jl")
    rethrow()

end


text:
function Primes(N::Int64)

	isprime = ones(Bool,N)

	for i = 2:int64(sqrt(N))
		if isprime[i]
			for j in (i*i):i:N
				isprime[j] = false
			end
		end
	end
	return filter(x -> isprime[x], 2:N)
end

function TrialDivision(N::Int64)

	primeFactors = [1]

	testPrimes = Primes(N);
	Nc = N;
	idx = 1;
	while (prod(primeFactors)!=N) && (idx < length(testPrimes))
		if mod(Nc, testPrimes[idx])==0
			println(testPrimes[idx])
			primeFactors=[primeFactors,testPrimes[idx]]
			Nc = Nc/testPrimes[idx]
		end
		idx = idx+1
	end

	return primeFactors[2:end]
end

function LargestPrimeFactor(N::Int64)

	first20 = Primes(1000)
	largest = 1
	for p in first20
		if mod(N, p) == 0
			largest = p
			N = fld(N,p)
		end
	end
	if N == 1
		return largest
	end
	trials = flipud(Primes(N+1))
	for p in trials
		if (mod(N,p) == 0)&&(p > largest)
			println(p)
			largest = p
			break
		end
	end
	if largest == 1
		return N
	else
		return largest
	end
end

function Euler3(N::Int64)

	big = LargestPrimeFactor(N)
	println("The largest prime factor of ",N," is ",big,".")

end

@time Euler3(600851475143)
#Euler3(100)

max_stars_repo_name: gribeill/ProjectEuler

max_issues_repo_name: gribeill/ProjectEuler

max_forks_repo_name: gribeill/ProjectEuler


text:
@testset "Intersections" begin
  @testset "Segments" begin
    # segments in 2D
    s1 = Segment(P2(0,0), P2(1,0))
    s2 = Segment(P2(0.5,0.0), P2(2,0))
    @test s1 ∩ s2 == Segment(P2(0.5,0.0), P2(1,0))
    @test s2 ∩ s1 == Segment(P2(0.5,0.0), P2(1,0))

    s1 = Segment(P2(0,0), P2(1,0))
    s2 = Segment(P2(0,0), P2(0,1))
    @test s1 ∩ s2 == P2(0,0)
    @test s2 ∩ s1 == P2(0,0)

    s1 = Segment(P2(0,0), P2(1,0))
    s2 = Segment(P2(0,0), P2(-1,0))
    @test s1 ∩ s2 == P2(0,0)
    @test s2 ∩ s1 == P2(0,0)

    s1 = Segment(P2(0,0), P2(0,1))
    s2 = Segment(P2(0,0), P2(0,-1))
    @test s1 ∩ s2 == P2(0,0)
    @test s2 ∩ s1 == P2(0,0)

    s1 = Segment(P2(1,1), P2(1,2))
    s2 = Segment(P2(1,1), P2(1,0))
    @test s1 ∩ s2 == P2(1,1)
    @test s2 ∩ s1 == P2(1,1)

    s1 = Segment(P2(1,1), P2(2,1))
    s2 = Segment(P2(1,0), P2(3,0))
    @test s1 ∩ s2 === nothing
    @test s2 ∩ s1 === nothing

    s1 = Segment(P2(0.181429364026879, 0.546811355144474),
                  P2(0.38282226144778, 0.107781953228536))
    s2 = Segment(P2(0.412498700935005, 0.212081819871479),
                  P2(0.395936725690311, 0.252041094122474))
    @test s1 ∩ s2 === nothing
    @test s2 ∩ s1 === nothing

    s1 = Segment(P2(1,2), P2(1,0))
    s2 = Segment(P2(1,0), P2(1,1))
    @test s1 ∩ s2 == Segment(P2(1,0), P2(1,1))

    s1 = Segment(P2(0,0), P2(2,0))
    s2 = Segment(P2(-2,0), P2(-1,0))
    s3 = Segment(P2(-1,0), P2(-2,0))
    @test s1 ∩ s2 === s2 ∩ s1 === nothing
    @test s1 ∩ s3 === s3 ∩ s1 === nothing

    s1 = Segment(P2(-1,0), P2(0,0))
    s2 = Segment(P2(0,0), P2(2,0))
    @test s1 ∩ s2 == s2 ∩ s1 == P2(0,0)

    s1 = Segment(P2(-1,0), P2(1,0))
    s2 = Segment(P2(0,0), P2(3,0))
    @test s1 ∩ s2 == s2 ∩ s1 == Segment(P2(0,0), P2(1,0))

    s1 = Segment(P2(0,0), P2(1,0))
    s2 = Segment(P2(0,0), P2(2,0))
    @test s1 ∩ s2 == s2 ∩ s1 == Segment(P2(0,0), P2(1,0))

    s1 = Segment(P2(0,0), P2(3,0))
    s2 = Segment(P2(1,0), P2(2,0))
    @test s1 ∩ s2 == s2 ∩ s1 == s2

    s1 = Segment(P2(0,0), P2(2,0))
    s2 = Segment(P2(1,0), P2(2,0))
    @test s1 ∩ s2 == s2 ∩ s1 == s2

    s1 = Segment(P2(0,0), P2(2,0))
    s2 = Segment(P2(1,0), P2(3,0))
    @test s1 ∩ s2 == s2 ∩ s1 == Segment(P2(1,0), P2(2,0))

    s1 = Segment(P2(0,0), P2(2,0))
    s2 = Segment(P2(2,0), P2(3,0))
    @test s1 ∩ s2 == s2 ∩ s1 == P2(2,0)

    s1 = Segment(P2(0,0), P2(2,0))
    s2 = Segment(P2(3,0), P2(4,0))
    @test s1 ∩ s2 === s2 ∩ s1 === nothing

    s1 = Segment(P2(2,1), P2(1,2))
    s2 = Segment(P2(1,0), P2(1,1))
    @test s1 ∩ s2 === s2 ∩ s1 === nothing

    s1 = Segment(P2(1.5,1.5), P2(3.0,1.5))
    s2 = Segment(P2(3.0,1.0), P2(2.0,2.0))
    @test s1 ∩ s2 == s2 ∩ s1 == P2(2.5,1.5)

    s1 = Segment(P2(0.94495744, 0.53224397), P2(0.94798386, 0.5344541))
    s2 = Segment(P2(0.94798386, 0.5344541), P2(0.9472896, 0.5340202))
    @test s1 ∩ s2 == s2 ∩ s1 == P2(0.94798386, 0.5344541)

    # segments in 3D
    s1 = Segment(P3(0.0, 0.0, 0.0), P3(1.0, 0.0, 0.0))
    s2 = Segment(P3(0.5, 1.0, 0.0), P3(0.5, -1.0, 0.0))
    s3 = Segment(P3(0.5, 0.0, 0.0), P3(1.5, 0.0, 0.0))
    s4 = Segment(P3(0.0, 1.0, 0.0), P3(0.0, -2.0, 0.0))
    s5 = Segment(P3(-1.0, 1.0, 0.0), P3(2.0, -2.0, 0.0))
    s6 = Segment(P3(0.0, 0.0, 0.0), P3(0.0, 1.0, 0.0))
    s7 = Segment(P3(-1.0, 1.0, 0.0), P3(-1.0, -1.0, 0.0))
    s8 = Segment(P3(-1.0, 1.0, 1.0), P3(-1.0, -1.0, 1.0))
    s9 = Segment(P3(0.5, 1.0, 1.0), P3(0.5, -1.0, 1.0))
    s10 = Segment(P3(0.0, 1.0, 0.0), P3(1.0, 1.0, 0.0))
    s11 = Segment(P3(1.5, 0.0, 0.0), P3(2.5, 0.0, 0.0))
    s12 = Segment(P3(1.0, 0.0, 0.0), P3(2.0, 0.0, 0.0))

    @test isa(intersecttype(s1, s2), CrossingSegments)        # CrossingSegments
    @test s1 ∩ s2 == P3(0.5, 0.0, 0.0)
    @test isa(intersecttype(s1, s3), OverlappingSegments)     # OverlappingSegments
    @test s1 ∩ s3 == Segment(P3(0.5, 0.0, 0.0), P3(1.0, 0.0, 0.0))
    @test isa(intersecttype(s1, s4), MidTouchingSegments)     # MidTouchingSegments (perpendicular)
    @test s1 ∩ s4 == P3(0.0, 0.0, 0.0)
    @test isa(intersecttype(s1, s5), MidTouchingSegments)     # MidTouchingSegments (obtuse)
    @test s1 ∩ s5 == P3(0.0, 0.0, 0.0)
    @test isa(intersecttype(s1, s6), CornerTouchingSegments)  # CornerTouchingSegments
    @test s1 ∩ s6 == P3(0.0, 0.0, 0.0)
    @test isa(intersecttype(s1, s7), NoIntersection)          # NoIntersection (but coplanar)
    @test isnothing(s1 ∩ s7)
    @test isa(intersecttype(s1, s8), NoIntersection)          # NoIntersection (non-coplanar)
    @test isnothing(s1 ∩ s8)
    @test isa(intersecttype(s1, s9), NoIntersection)          # NoIntersection (non-coplanar)
    @test isnothing(s1 ∩ s9)
    @test isa(intersecttype(s1, s10), NoIntersection)         # NoIntersection (parallel)
    @test isnothing(s1 ∩ s10)
    @test isa(intersecttype(s1, s11), NoIntersection)         # NoIntersection (colinear, not-overlapping)
    @test isnothing(s1 ∩ s11)
    @test isa(intersecttype(s1, s12), CornerTouchingSegments) # CornerTouchingSegments (colinear)
    @test s1 ∩ s12 == P3(1.0, 0.0, 0.0)
  end

  @testset "Lines" begin
    l1 = Line(P2(0,0), P2(1,0))
    l2 = Line(P2(-1,-1), P2(-1,1))
    @test l1 ∩ l2 == l2 ∩ l1 == P2(-1,0)

    l1 = Line(P2(0,0), P2(1,0))
    l2 = Line(P2(0,1), P2(1,1))
    @test l1 ∩ l2 === l2 ∩ l1 === nothing

    l1 = Line(P2(0,0), P2(1,0))
    l2 = Line(P2(1,0), P2(2,0))
    @test l1 == l2
    @test l1 ∩ l2 == l2 ∩ l1 == l1
  end

  @testset "Boxes" begin
    b1 = Box(P2(0,0), P2(1,1))
    b2 = Box(P2(0.5,0.5), P2(2,2))
    b3 = Box(P2(2,2), P2(3,3))
    b4 = Box(P2(1,1), P2(2,2))
    b5 = Box(P2(1.0,0.5), P2(2,2))
    @test intersecttype(b1, b2) isa OverlappingBoxes
    @test b1 ∩ b2 == Box(P2(0.5,0.5), P2(1,1))
    @test intersecttype(b1, b3) isa NoIntersection
    @test b1 ∩ b3 === nothing
    @test intersecttype(b1, b4) isa CornerTouchingBoxes
    @test b1 ∩ b4 == P2(1,1)
    @test intersecttype(b1, b5) isa FaceTouchingBoxes
    @test b1 ∩ b5 == Box(P2(1.0,0.5), P2(1,1))
  end

  @testset "Misc" begin
    t = Triangle(P2[(0,0),(1,0),(0,1)])
    q = Quadrangle(P2[(1,1),(2,1),(2,2),(1,2)])
    @test hasintersect(t, t)
    @test hasintersect(q, q)
    @test !hasintersect(t, q)
    @test !hasintersect(q, t)

    t = Triangle(P2[(1,0),(2,0),(1,1)])
    q = Quadrangle(P2[(1.3,0.5),(2.3,0.5),(2.3,1.5),(1.3,1.5)])
    @test hasintersect(t, t)
    @test hasintersect(q, q)
    @test hasintersect(t, q)
    @test hasintersect(q, t)

    t = Triangle(P2[(0,0),(1,0),(0,1)])
    b = Ball(P2(0,0), T(1))
    @test hasintersect(t, t)
    @test hasintersect(b, b)
    @test hasintersect(t, b)
    @test hasintersect(b, t)

    t = Triangle(P2[(1,0),(2,0),(1,1)])
    b = Ball(P2(0,0), T(1))
    @test hasintersect(t, t)
    @test hasintersect(b, b)
    @test hasintersect(t, b)
    @test hasintersect(b, t)

    t = Triangle(P2[(1,0),(2,0),(1,1)])
    b = Ball(P2(-0.01,0), T(1))
    @test hasintersect(t, t)
    @test hasintersect(b, b)
    @test !hasintersect(t, b)
    @test !hasintersect(b, t)

    outer = P2[(0,0),(1,0),(1,1),(0,1),(0,0)]
    hole1 = P2[(0.2,0.2),(0.4,0.2),(0.4,0.4),(0.2,0.4),(0.2,0.2)]
    hole2 = P2[(0.6,0.2),(0.8,0.2),(0.8,0.4),(0.6,0.4),(0.6,0.2)]
    poly1 = PolyArea(outer)
    poly2 = PolyArea(outer, [hole1, hole2])
    ball1 = Ball(P2(0.5,0.5), T(0.05))
    ball2 = Ball(P2(0.3,0.3), T(0.05))
    ball3 = Ball(P2(0.7,0.3), T(0.05))
    ball4 = Ball(P2(0.3,0.3), T(0.15))
    @test hasintersect(poly1, poly1)
    @test hasintersect(poly2, poly2)
    @test hasintersect(poly1, poly2)
    @test hasintersect(poly2, poly1)
    @test hasintersect(poly1, ball1)
    @test hasintersect(poly2, ball1)
    @test hasintersect(poly1, ball2)
    @test !hasintersect(poly2, ball2)
    @test hasintersect(poly1, ball3)
    @test !hasintersect(poly2, ball3)
    @test hasintersect(poly1, ball4)
    @test hasintersect(poly2, ball4)
    mesh1 = discretize(poly1, Dehn1899())
    mesh2 = discretize(poly2, Dehn1899())
    @test hasintersect(mesh1, mesh1)
    @test hasintersect(mesh2, mesh2)
    @test hasintersect(mesh1, mesh2)
    @test hasintersect(mesh2, mesh1)
  end
end


text:
IP_legends_en = Dict("title" => "Information Plane (MNEAT)", "cbartitle" => "Generation")
IP_legends_pt = Dict("title" => "Plano de informação (MNEAT)", "cbartitle" => "Geração")

fitness_legends_en = Dict("title" => "Evolution of fitness", "ylabel" => "Fitness",
                    "xlabel" => "Generation", "label1" => "Best fitness",
                    "label2" => "Population's average fitness")
fitness_legends_pt = Dict("title" => "Evolução do fitness", "ylabel" => "Fitness",
                    "xlabel" => "Geração", "label1" => "Melhor fitness",
                    "label2" => "Fitness médio da população")

output_dist_legends_en = Dict("title1" => "Best MLP distribution", "title2" => "Expected distribution",
                        "ylabel" => "Probability", "xlabel" => "Class")
output_dist_legends_pt = Dict("title1" => "Distribuição da melhor MLP", "title2" => "Distribuição esperada",
                        "ylabel" => "Probabilidade", "xlabel" => "Classe")

complexity_legends_en = Dict("title" => "Topological complexity evolution (MNEAT)",
                         "ylabel" => "Number of weights", "xlabel" => "Generation")
complexity_legends_pt = Dict("title" => "Evolução da complexidade topológica (MNEAT)",
                        "ylabel" => "Número de pesos", "xlabel" => "Geração")

accuracy_legends_en = Dict("title" => "Accuracy evolution (MNEAT)",
                            "ylabel" => "Accuracy", "xlabel" => "Generation",
                            "label1" => "Training set accuracy", "label2" => "Test set accuracy")
accuracy_legends_pt = Dict("title" => "Evolução da acurácia (MNEAT)",
                      "ylabel" => "Acurácia", "xlabel" => "Geração",
                      "label1" => "Acurácia do conjunto de treinamento",
                      "label2" => "Acurácia do conjunto de teste")

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
"""
An aggregator interface for SSA-like algorithms.

### Required Fields
- `next_jump`          # the next jump to execute
- `prev_jump`          # the previous jump that was executed
- `next_jump_time`     # the time of the next jump
- `end_time`           # the time to stop a simulation
- `cur_rates`          # vector of current propensity values
- `sum_rate`           # sum of current propensity values
- `ma_jumps`           # any MassActionJumps for the system (scalar form)
- `rates`              # vector of rate functions for ConstantRateJumps
- `affects!`           # vector of affect functions for ConstantRateJumps
- `save_positions`     # tuple for whether to save the jumps before and/or after event
- `rng`                # random number generator

### Optional fields:
- `dep_gr`             # dependency graph, dep_gr[i] = indices of reactions that should
                       # be updated when rx i occurs.
"""
abstract type AbstractSSAJumpAggregator <: AbstractJumpAggregator end

DiscreteCallback(c::AbstractSSAJumpAggregator) = DiscreteCallback(c, c, initialize = c, save_positions = c.save_positions)

########### The following routines are templates for all SSAs ###########
########### Generally they should not need to be overloaded.  ###########

## Users will normally define (see direct.jl for examples):
# aggregate
# initialize!
# execute_jumps!
# generate_jumps!

# condition for jump to occur
@inline function (p::AbstractSSAJumpAggregator)(u, t, integrator)
    p.next_jump_time == t
end

# executing jump at the next jump time
function (p::AbstractSSAJumpAggregator)(integrator)
    execute_jumps!(p, integrator, integrator.u, integrator.p, integrator.t)
    generate_jumps!(p, integrator, integrator.u, integrator.p, integrator.t)
    register_next_jump_time!(integrator, p, integrator.t)
    nothing
end

# setting up a new simulation
function (p::AbstractSSAJumpAggregator)(dj, u, t, integrator) # initialize
    initialize!(p, integrator, u, integrator.p, t)
    register_next_jump_time!(integrator, p, integrator.t)
    u_modified!(integrator,false)
    nothing
end

############################## Generic Routines ###############################

"""
    register_next_jump_time!(integrator, p::AbstractSSAJumpAggregator, t)

Adds a `tstop` to the integrator at the next jump time.
"""
@inline function register_next_jump_time!(integrator, p::AbstractSSAJumpAggregator, t)
    if p.next_jump_time < p.end_time
        add_tstop!(integrator, p.next_jump_time)
    end
    nothing
end

"""
    build_jump_aggregation(jump_agg_type, u, p, t, end_time, ma_jumps, rates,
                           affects!, save_positions, rng; kwargs...)

Helper routine for setting up standard fields of SSA jump aggregations.
"""
function build_jump_aggregation(jump_agg_type, u, p, t, end_time, ma_jumps, rates,
                                affects!, save_positions, rng; kwargs...)

    # mass action jumps
    majumps = ma_jumps
    if majumps === nothing
        majumps = MassActionJump(Vector{typeof(t)}(),
                             Vector{Vector{Pair{Int,eltype(u)}}}(),
                             Vector{Vector{Pair{Int,eltype(u)}}}())
    end

    # current jump rates, allows mass action rates and constant jumps
    cur_rates = Vector{typeof(t)}(undef, get_num_majumps(majumps) + length(rates))

    sum_rate = zero(typeof(t))
    next_jump = 0
    next_jump_time = typemax(typeof(t))
    jump_agg_type(next_jump, next_jump_time, end_time, cur_rates, sum_rate,
                majumps, rates, affects!, save_positions, rng; kwargs...)
end

"""
    fill_rates_and_sum!(p::AbstractSSAJumpAggregator, u, params, t)

Reevaluate all rates and their sum.
"""
function fill_rates_and_sum!(p::AbstractSSAJumpAggregator, u, params, t)
    sum_rate = zero(typeof(p.sum_rate))

    # mass action jumps
    majumps   = p.ma_jumps
    cur_rates = p.cur_rates
    @inbounds for i in 1:get_num_majumps(majumps)
        cur_rates[i] = evalrxrate(u, i, majumps)
        sum_rate    += cur_rates[i]
    end

    # constant rates
    rates = p.rates
    idx   = get_num_majumps(majumps) + 1
    @inbounds for rate in rates
        cur_rates[idx] = rate(u, params, t)
        sum_rate += cur_rates[idx]
        idx += 1
    end

    p.sum_rate = sum_rate
    nothing
end

"""
    calculate_jump_rate(ma_jumps, rates, u, params, t, rx)

Recalculate the rate for the jump with index `rx`.
"""
@inline function calculate_jump_rate(ma_jumps, num_majumps, rates, u, params, t, rx)
    if rx <= num_majumps
        return evalrxrate(u, rx, ma_jumps)
    else
        @inbounds return rates[rx - num_majumps](u, params, t)
    end
end

"""
    update_dependent_rates!(p::AbstractSSAJumpAggregator, u, params, t)

Recalculate jump rates for jumps that depend on the just executed jump.

Notes:
    - Intended for methods that have a dependency graph, i.e. define `p.dep_gr`.
"""
function update_dependent_rates!(p::AbstractSSAJumpAggregator, u, params, t)
    @inbounds dep_rxs = p.dep_gr[p.next_jump]
    cur_rates   = p.cur_rates
    sum_rate    = p.sum_rate
    num_majumps = get_num_majumps(p.ma_jumps)
    @inbounds for rx in dep_rxs
        sum_rate -= cur_rates[rx]
        @inbounds cur_rates[rx] = calculate_jump_rate(p.ma_jumps, num_majumps, p.rates, u,params,t,rx)
        sum_rate += cur_rates[rx]
    end

    p.sum_rate = sum_rate
    nothing
end

"""
    update_state!(p::AbstractSSAJumpAggregator, integrator, u)

Execute `p.next_jump`.
"""
@inline function update_state!(p::AbstractSSAJumpAggregator, integrator, u)
    @unpack ma_jumps, next_jump = p
    num_ma_rates = get_num_majumps(ma_jumps)
    if next_jump <= num_ma_rates # is next jump a mass action jump
        if u isa SVector
            integrator.u = executerx(u, next_jump, ma_jumps)
        else
            @inbounds executerx!(u, next_jump, ma_jumps)
        end
    else
        idx = next_jump - num_ma_rates
        @inbounds p.affects![idx](integrator)
    end

    # save jump that was just executed
    p.prev_jump = next_jump
    return integrator.u
end

"""
    nomorejumps!(p, sum_rate) :: Bool

Check if the total rate is zero, and if it is, make the next jump time Inf.
"""
@inline function nomorejumps!(p, sum_rate) :: Bool
    if sum_rate < eps(typeof(sum_rate))
        p.next_jump      = zero(p.next_jump)
        p.next_jump_time = convert(typeof(sum_rate), Inf)
        return true
    end
    return false
end

"""
    linear_search(iterator, r)

Perform linear search for `r` over iterator. Output index j s.t. sum(array[1:j-1])
< r <= sum(array[1:j]), where array=collect(iterator)

Notes:
- The iterator must have all positive numbers
- Returns index zero if the search is unsuccessful. Assumes this corresponds to
  the case of an infinite next reaction time and so the jump index does not
  matter.
"""
@inline function linear_search(iterator, r)
    jidx = 0
    idx = 1
    for element in iterator
        r -= element
        if r < zero(r)
            jidx = idx
            break
        end
        idx += 1
    end
    return jidx
end

"""
    rejectrx(ma_jumps, rates, cur_rate_high, cur_rate_low, rng, u, jidx, params, t)

Perform rejection sampling test (used in RSSA methods).
"""
@inline function rejectrx(ma_jumps, num_majumps, rates, cur_rate_high, cur_rate_low, rng, u, jidx, params, t)
    # rejection test
    @inbounds r2     = rand(rng) * cur_rate_high[jidx]
    @inbounds crlow  = cur_rate_low[jidx]

    if crlow > zero(crlow) && r2 <= crlow
        return false
    else
        # calculate actual propensity, split up for type stability
        crate = calculate_jump_rate(ma_jumps, num_majumps, rates, u, params, t, jidx)
        if crate > zero(crate) && r2 <= crate
            return false
        end
    end
    return true
end


text:
# Copyright © 2021, UChicago Argonne, LLC
# All Rights Reserved.

# This source code is licensed under the BSD-style license found in the
# LICENSE file in the root directory of this source tree.

using DataFrames, DataStructures, InteractiveUtils, Logging, Printf
import CSV, DelimitedFiles, SCIP, Statistics, Glob
using JuMP, MathOptInterface
const MOI = MathOptInterface

max_stars_repo_name: jongeunkim/STreCH

max_issues_repo_name: jongeunkim/STreCH

max_forks_repo_name: jongeunkim/STreCH


text:
# using Revise, Plots, Test
using PseudoArcLengthContinuation, LinearAlgebra, Setfield, SparseArrays, ForwardDiff, Parameters
const PALC = PseudoArcLengthContinuation
norminf = x -> norm(x, Inf)

function Fbp(x, p)
	return [x[1] * (3.23 .* p.μ - 0.12 * x[1] + 0.234 * x[1]^2) + x[2], x[2]]
end

par = (μ = -0.2, ν = 0)
####################################################################################################
opt_newton = PALC.NewtonPar(tol = 1e-8, verbose = false, maxIter = 20)
opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = 0.01, pMax = 0.1, pMin = -0.3, detectBifurcation = 2, nev = 2, newtonOptions = opt_newton, maxSteps = 100)

	br, _ = @time PALC.continuation(
		Fbp, [0.1, 0.1], par, (@lens _.μ),
		printSolution = (x, p) -> norminf(x),
		opts_br; plot = false, verbosity = 0, normC = norminf)

####################################################################################################
# normal form computation
D(f, x, p, dx) = ForwardDiff.derivative(t->f(x .+ t .* dx, p), 0.)

d1F(x,p,dx1)         = D((z, p0) -> Fbp(z, p0), x, p, dx1)
d2F(x,p,dx1,dx2)     = D((z, p0) -> d1F(z, p0, dx1), x, p, dx2)
d3F(x,p,dx1,dx2,dx3) = D((z, p0) -> d2F(z, p0, dx1, dx2), x, p, dx3)
jet = (Fbp,
	(x, p) -> PALC.finiteDifferences(z -> Fbp(z, p), x),
	(x, p, dx1, dx2) -> d2F(x, p, dx1, dx2),
	(x, p, dx1, dx2, dx3) -> d3F(x, p, dx1, dx2, dx3))

bp = PALC.computeNormalForm(jet..., br, 1; verbose=true)

# normal form
nf = bp.nf

@test norm(nf[1]) < 1e-10
	@test norm(nf[2] - 3.23) < 1e-10
	@test norm(nf[3]/2 - -0.12) < 1e-10
	@test norm(nf[4]/6 - 0.234) < 1e-10

####################################################################################################
# Automatic branch switching
br, _ = continuation(jet..., br, 1, opts_br; verbosity = 0)

####################################################################################################
function Fbp2d(x, p)
	return [ x[1] * (3.23 .* p.μ - 0.123 * x[1]^2 - 0.234 * x[2]^2),
			 x[2] * (3.23 .* p.μ - 0.456 * x[1]^2 - 0.123 * x[2]^2),
			 x[3]]
end

d1F2d(x,p,dx1) = D((z, p0) -> Fbp2d(z, p0), x, p, dx1)
	d2F2d(x,p,dx1,dx2)     = D((z, p0) -> d1F2d(z, p0, dx1), x, p, dx2)
	d3F2d(x,p,dx1,dx2,dx3) = D((z, p0) -> d2F2d(z, p0, dx1, dx2), x, p, dx3)

jet = (Fbp2d, (x, p) -> ForwardDiff.jacobian(z -> Fbp2d(z, p), x), d2F2d, d3F2d)

par = (μ = -0.2, ν = 0)

br, _ = @time PALC.continuation(
	Fbp2d, [0.01, 0.01, 0.01], par, (@lens _.μ),
	printSolution = (x, p) -> norminf(x),
	setproperties(opts_br; nInversion = 2); plot = false, verbosity = 0, normC = norminf)

# we have to be careful to have the same basis as for Fbp2d or the NF will not match Fbp2d
bp2d = @time PALC.computeNormalForm(jet..., br, 1; ζs = [[1, 0, 0.], [0, 1, 0.]]);

PALC.nf(bp2d)
bp2d(rand(2), 0.2)
bp2d(Val(:reducedForm), rand(2), 0.2)

@test abs(bp2d.nf.b3[1,1,1,1] / 6 - -0.123) < 1e-10
@test abs(bp2d.nf.b3[1,1,2,2] / 2 - -0.234) < 1e-10
@test abs(bp2d.nf.b3[1,1,1,2] / 2 - -0.0)   < 1e-10
@test abs(bp2d.nf.b3[2,1,1,2] / 2 - -0.456) < 1e-10
@test norm(bp2d.nf.b2, Inf) < 3e-6
@test norm(bp2d.nf.b1 - 3.23 * I, Inf) < 1e-10
@test norm(bp2d.nf.a, Inf) < 1e-6
####################################################################################################
# test of the Hopf normal form
function Fsl2!(f, u, p, t)
	@unpack r, μ, ν, c3, c5 = p
	u1 = u[1]
	u2 = u[2]
	ua = u1^2 + u2^2

	f[1] = r * u1 - ν * u2 - ua * (c3 * u1 - μ * u2) - c5 * ua^2 * u1
	f[2] = r * u2 + ν * u1 - ua * (c3 * u2 + μ * u1) - c5 * ua^2 * u2

	return f
end

Fsl2(x, p) = Fsl2!(similar(x), x, p, 0.)
par_sl = (r = 0.5, μ = 0.132, ν = 1.0, c3 = 1.123, c5 = 0.2)

d1Fsl(x,p,dx1)         = D((z, p0) -> Fsl2(z, p0), x, p, dx1)
d2Fsl(x,p,dx1,dx2)     = D((z, p0) -> d1Fsl(z, p0, dx1), x, p, dx2)
d3Fsl(x,p,dx1,dx2,dx3) = D((z, p0) -> d2Fsl(z, p0, dx1, dx2), x, p, dx3)

# detect hopf bifurcation
opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.02, ds = 0.01, pMax = 0.1, pMin = -0.3, detectBifurcation = 2, nev = 2, newtonOptions = (@set opt_newton.verbose = true), maxSteps = 100)

br, _ = @time PALC.continuation(
	Fsl2, [0.0, 0.0], (@set par_sl.r = -0.1), (@lens _.r),
	printSolution = (x, p) -> norminf(x),
	opts_br; plot = false, verbosity = 3, normC = norminf)

hp = PALC.computeNormalForm(
	(x, p) -> Fsl2(x, p),
	(x, p) -> ForwardDiff.jacobian(z -> Fsl2(z, p), x),
	(x, p, dx1, dx2) -> 	 d2Fsl(x, p, dx1, dx2),
	(x, p, dx1, dx2, dx3) -> d3Fsl(x, p, dx1, dx2, dx3),
	br, 1)

nf = hp.nf

@test abs(nf.a - 1) < 1e-10
@test abs(nf.b/2 - (-par_sl.c3 + im*par_sl.μ)) < 1e-10


text:
abstract type AbstractEigenfunction end

abstract type AbstractSphericalEigenfunction <: AbstractEigenfunction end

abstract type AbstractPlanarEigenfunction <: AbstractEigenfunction end

###
### Spherical eigenfunctions
###

mutable struct SphericalVertexEigenfunction{T<:Union{fmpq,arb}} <:
               AbstractSphericalEigenfunction
    domain::SphericalTriangle{T}
    vertex::Int
    stride::Int
    coefficients::Vector{arb}
end

mutable struct SphericalInteriorEigenfunction{T<:Union{fmpq,arb}} <:
               AbstractSphericalEigenfunction
    domain::SphericalTriangle{T}
    θ::arb # θ value for the interior point
    ϕ::arb # ϕ value for the interior point
    stride::Int
    coefficients::Vector{arb}
end

mutable struct SphericalCombinedEigenfunction{T<:Union{fmpq,arb}} <:
               AbstractSphericalEigenfunction
    domain::SphericalTriangle{T}
    us::Vector{<:AbstractSphericalEigenfunction}
    orders::Vector{Int}
end

mutable struct KrewerasEigenfunction{T<:Union{fmpq,arb}} <: AbstractSphericalEigenfunction
    domain::SphericalTriangle{T}
    coefficients::Vector{arb}
end

###
### Planar eigenfunctions
###

struct StandaloneVertexEigenfunction{T<:Union{fmpq,arb}} <: AbstractPlanarEigenfunction
    vertex::SVector{2,arb}
    orientation::T
    θ::T
    stride::Int
    offset::Int
    reversed::Bool
    coefficients::Vector{arb}
    parent::ArbField

    function StandaloneVertexEigenfunction(
        vertex::SVector{2,arb},
        orientation::T,
        θ::T;
        stride::Integer = 1,
        offset::Integer = 0,
        reversed::Bool = false,
        parent::ArbField = parent(vertex[1]),
    ) where {T<:Union{Rational,arb,fmpq}}
        if T == fmpq || T <: Rational
            orientation = fmpq(orientation)
            S = fmpq
        else
            S = arb
        end
        return new{S}(vertex, orientation, θ, stride, offset, reversed, arb[], parent)
    end
end

""""
    StandaloneInteriorEigenfunction(vertex, orientation, parent; stride = 1, even = false)

An eigenfunction consisting of the functions `bessel_j(ν,
r*√λ)*sin(ν*θ)` and `bessel_j(ν, ν*√λ)*cos(j*θ)` for `ν` = 0, 1, 2,
... Here `r` and `θ` are polar coordinates around `vertex` take so
that `θ = 0` at the given orientation.

The values of `ν` that are used are 0, `stride`, `2stride`,
`3stride`,...

If `even` is true then use only the function with `cos(j*θ)`.
"""
struct StandaloneInteriorEigenfunction{T<:Union{fmpq,arb}} <: AbstractPlanarEigenfunction
    vertex::SVector{2,arb}
    orientation::T
    stride::Int
    offset::Int
    even::Bool
    odd::Bool
    coefficients::Vector{arb}
    parent::ArbField

    function StandaloneInteriorEigenfunction(
        vertex::SVector{2,arb},
        orientation::T = fmpq(0),
        parent::ArbField = parent(first(vertex));
        offset::Integer = 0,
        stride::Integer = 1,
        even::Bool = false,
        odd::Bool = false,
    ) where {T<:Union{Rational,fmpq,arb}}
        even && odd && throw(ArgumentError("eigenfunction can't be both even and odd"))
        if T == fmpq || T <: Rational
            orientation = fmpq(orientation)
            S = fmpq
        else
            S = arb
        end
        return new{S}(vertex, orientation, stride, offset, even, odd, arb[], parent)
    end
end

struct StandaloneLightningEigenfunction{T<:Union{AbstractFloat,arb},S<:Union{fmpq,arb}} <:
       AbstractPlanarEigenfunction
    vertex::SVector{2,T}
    orientation::S
    θ::S
    l::T
    σ::T
    even::Bool
    odd::Bool
    reversed::Bool
    coefficients::Vector{T}
    parent::ArbField

    function StandaloneLightningEigenfunction(
        vertex::SVector{2,arb},
        orientation::S,
        θ::S,
        parent::ArbField = parent(vertex[1]);
        l::arb = parent(1),
        σ::arb = parent(4),
        even::Bool = false,
        odd::Bool = false,
        reversed::Bool = false,
    ) where {S<:Union{arb,fmpq}}
        even && odd && throw(ArgumentError("eigenfunction can't be both even and odd"))
        return new{arb,S}(vertex, orientation, θ, l, σ, even, odd, reversed, arb[], parent)
    end

    function StandaloneLightningEigenfunction{arb,S}(
        vertex::SVector{2,arb},
        orientation::S,
        θ::S,
        parent::ArbField = parent(vertex[1]);
        l = 1,
        σ = 4,
        even::Bool = false,
        odd::Bool = false,
        reversed::Bool = false,
    ) where {S<:Union{arb,fmpq}}
        even && odd && throw(ArgumentError("eigenfunction can't be both even and odd"))
        return new{arb,S}(
            vertex,
            orientation,
            θ,
            parent(l),
            parent(σ),
            even,
            odd,
            reversed,
            arb[],
            parent,
        )
    end

    function StandaloneLightningEigenfunction{T,S}(
        vertex::AbstractVector,
        orientation::S,
        θ::S;
        l = one(T),
        σ = 4one(T),
        even::Bool = false,
        odd::Bool = false,
        reversed::Bool = false,
    ) where {T<:AbstractFloat,S<:Union{arb,fmpq}}
        even && odd && throw(ArgumentError("eigenfunction can't be both even and odd"))
        return new{T,S}(
            convert(SVector{2,T}, vertex),
            orientation,
            θ,
            convert(T, l),
            convert(T, σ),
            even,
            odd,
            reversed,
            T[],
            ArbField(precision(T)),
        )
    end
end

struct LinkedEigenfunction{T<:AbstractPlanarEigenfunction} <: AbstractPlanarEigenfunction
    us::Vector{T}
    extra_coefficients::Vector{arb}
    excluded_boundaries::Vector{BitSet}

    function LinkedEigenfunction(
        us::Vector{T},
        extra_coefficients::Vector = ones(length(us));
        excluded_boundaries::Vector = fill(BitSet(), length(us)),
    ) where {T<:AbstractPlanarEigenfunction}
        isempty(us) && throw(ArgumentError("us must not be empty"))
        length(us) == length(extra_coefficients) ||
            throw(ArgumentError("us and extra_coefficients must have the same size"))

        extra_coefficients = first(us).parent.(extra_coefficients)

        return new{T}(us, extra_coefficients, excluded_boundaries)
    end
end

mutable struct CombinedEigenfunction <: AbstractPlanarEigenfunction
    domain::AbstractPlanarDomain
    us::Vector{<:AbstractPlanarEigenfunction}
    boundary_to_us::OrderedDict{Int,BitSet}
    even_boundaries::BitSet
    orders::Vector{Int}
end

mutable struct LShapeEigenfunction <: AbstractPlanarEigenfunction
    domain::LShape
    stride::Int
    coefficients::Vector{arb}
end


text:
# defines how to fold detailed score into N components
abstract type DetailedScoreFolding{N} end

# soft Fold2D score
const FoldedScore = NTuple{2, Float64}

BlackBoxOptim.fitness_scheme_type(scorefold::DetailedScoreFolding) =
    fitness_scheme_type(typeof(scorefold))

BlackBoxOptim.fitness_scheme_type(::Type{SF}) where SF <: DetailedScoreFolding{N} where N =
    ParetoFitnessScheme{N, Float64, true, MultiobjCoverProblemScoreAggregator{SF}}

struct MultiobjCoverProblemScoreAggregator{SF <: DetailedScoreFolding}
    score_folding::SF
    setXset_factor::Float64
    uncovered_factor::Float64
    covered_factor::Float64

    MultiobjCoverProblemScoreAggregator(params::CoverParams,
                                        score_folding::SF) where SF =
        new{SF}(score_folding,
                params.setXset_factor,
                params.uncovered_factor,
                params.covered_factor)
end

#=
# no folding of fitness components
struct MultiobjProblemNoFolding <: DetailedScoreFolding{4}
end
(::MultiobjProblemNoFolding)(fitness::DetailedScore) = fitness

(agg::MultiobjCoverProblemScoreAggregator{MultiobjProblemNoFolding})(score::DetailedScore) =
    score[1] + agg.setXset_factor * score[2] +
    agg.uncovered_factor * score[3] + agg.covered_factor * score[4]
=#

"""
Transforms the 4-component cover quality score into 2-component score that makes
the highly redundant solutions dominated by any less redundant ones.

See ["Cover score convolution"](@ref cover_score_convolution) section for the discussion.

# Arguments
  * `setXset_factor`: ``w_r``, same as in [`CoverParams`](@ref)
  * `uncovered_factor`: ``w_u``, same as in [`CoverParams`](@ref)
  * `covered_factor`: ``w_c``, same as in [`CoverParams`](@ref)
  * `ratio_threshold`: ``k_{\\max}``, defaults to 1
  * `shape`: ``\\alpha_k``, defaults to 0.5.
"""
struct MultiobjProblemSoftFold2d <: DetailedScoreFolding{2}
    setXset_factor::Float64
    uncovered_factor::Float64
    covered_factor::Float64
    ratio_threshold::Float64    # at what set_score/setXset_score ratio the folding kicks in
    shape::Float64              # controls folding smoothness (smaller are smoother)

    MultiobjProblemSoftFold2d(params::CoverParams,
                              ratio_threshold::Float64 = 1.0,
                              shape::Float64=0.5) =
        new(params.setXset_factor, params.uncovered_factor, params.covered_factor,
            ratio_threshold, shape)
end

function (f::MultiobjProblemSoftFold2d)(fitness::DetailedScore)
    a = fitness[1] +
        f.uncovered_factor * fitness[3] +
        f.covered_factor * fitness[4]
    b = fitness[2]
    if f.ratio_threshold !== nothing
        k = b/max(1, -a)
        (k <= f.ratio_threshold) && return (a, b)

        s = f.shape * (k - f.ratio_threshold)
        # 1/(1 + exp(1/s - s)) allows smooth transition from s = 0 (f(0)=0) to f(∞)=1
        t = inv(k * (f.setXset_factor + exp(inv(s) - s)))
        return (a + t * f.setXset_factor * b, (1.0 - t) * b)
    else
        return (a, b)
    end
end

(agg::MultiobjCoverProblemScoreAggregator{MultiobjProblemSoftFold2d})(score::FoldedScore) =
    score[1] + agg.setXset_factor * score[2]

#=
# Sums uncovered (3rd) and covered (4th) components of the score into a single score
# that to every score component mixes in the other 2 components with the score k
# The transform keeps the aggregated score intact
struct MultiobjProblemFold3d <: DetailedScoreFolding{3}
    k::Float64
    k_sXs::Float64
    k_u::Float64
    setXset_factor::Float64
    uncovered_factor::Float64
    covered_factor::Float64

    MultiobjProblemFold3d(params::CoverParams, k::Float64=0.1) =
        new(params.setXset_factor > 0.0 ? k : 0.0,
            params.setXset_factor > 0.0 ? k/params.setXset_factor : 0.0,
            params.uncovered_factor > 0.0 ? k/params.uncovered_factor : 0.0,
            params.setXset_factor,
            params.uncovered_factor,
            params.uncovered_factor > 0.0 ? params.covered_factor/params.uncovered_factor : 1.0)
end

function (f::MultiobjProblemFold3d)(fitness::DetailedScore)
    a = fitness[1]
    b = fitness[2]
    c = fitness[3] + f.covered_factor * fitness[4]
    return (0.5 * f.k * (f.setXset_factor * b + f.uncovered_factor * c) + (1-f.k) * a,
            0.5 * f.k_sXs * (a + f.uncovered_factor * c) + (1-f.k) * b,
            0.5 * f.k_u * (a + f.setXset_factor * b) + (1-f.k) * c)
end

(agg::MultiobjCoverProblemScoreAggregator{MultiobjProblemFold3d})(score::NTuple{3, Float64}) =
    score[1] + agg.setXset_factor * score[2] + agg.uncovered_factor * score[3]
=#

"""
Abstract multi-objective optimal Enriched-Set Cover problem.

See ["Method Description"](@ref method) for more details.
"""
abstract type MultiobjCoverProblem <: AbstractCoverProblem{DetailedScore} end

function score_scales(problem::MultiobjCoverProblem)
    w_min, w_max = extrema(problem.var_scores)
    sXs_min, sXs_max = length(problem.varXvar_scores) > 1 ? (Inf, -Inf) : (0.0, 0.0)
    return (max(w_max - w_min, 0.1),
            max(sXs_max - sXs_min, 0.1))
end

aggscore(w::AbstractVector{Float64}, problem::MultiobjCoverProblem) =
    aggscore(score(w, problem), problem.params)

# exclude given vars from the problem and
# adjust the remaining scores
function mask_vars(problem::MultiobjCoverProblem,
                   vars::AbstractVector{Int};
                   penalize_overlaps::Bool = true)
    varmask = fill(true, nvars(problem))
    varmask[vars] .= false
    v_scores = problem.var_scores[varmask]
    if penalize_overlaps
        # penalize overlapping sets
        pweights = fill(0.0, nvars(problem))
        pweights[vars] .= 1.0
        v_penalties = similar(pweights)
        minplus_bilinear!(take2, v_penalties,
                          problem.varXvar_scores,
                          pweights, pweights)
        v_scores .-= view(v_penalties, varmask) .* problem.params.setXset_factor
    end
    return varmask, v_scores
end

"""
Multi-objective optimal Enriched-Set Cover problem that uses MaskedSetMosaic
to define the coverage scores based on the set weights and masked elements coverage.

See ["Method Description"](@ref method) for more details.
"""
struct MaskedSetCoverProblem <: MultiobjCoverProblem
    params::CoverParams

    var2set::Vector{Int}

    var_scores::Vector{Float64}
    varXvar_scores::Matrix{Float64}

    tileXvar::SparseMaskMatrix      # tile X set(var)
    maskxtileXvar::SparseMatrixCSC{Float64} # tile-in-mask X set(var)
    nmasked_tile::Vector{Int}       # sum of masked elements for each tile in all masks (els counted for each mask separately)
    nunmasked_tile::Vector{Int}     # number of unmasked elements for each tile in the union of all masks (each el counted once)

    arrpool::ArrayPool{Float64}   # pool of vectors to reuse for evaluation FIXME move to evaluator

    function MaskedSetCoverProblem(params::CoverParams,
                        var2set::AbstractVector{Int},
                        var_scores::AbstractVector{Float64},
                        varXvar_scores::AbstractMatrix{Float64},
                        tileXvar::SparseMaskMatrix,
                        maskxtileXvar::SparseMatrixCSC,
                        nmasked_tile::AbstractVector{Int},
                        nunmasked_tile::AbstractVector{Int}
    )
        length(var2set) == length(var_scores) ==
        size(varXvar_scores, 1) == size(varXvar_scores, 2) == size(maskxtileXvar, 2) ==
        size(tileXvar, 2) ||
            throw(ArgumentError("var2set, var_scores and varXvar_scores counts do not match"))
        size(tileXvar, 1) == length(nmasked_tile) == length(nunmasked_tile) ||
            throw(ArgumentError("tileXvar and nmasked_tile rows count does not match"))
        new(params, var2set,
            var_scores, varXvar_scores,
            tileXvar, maskxtileXvar, nmasked_tile, nunmasked_tile,
            ArrayPool{Float64}())
    end
end

function MultiobjCoverProblem(mosaic::MaskedSetMosaic, params::CoverParams = CoverParams())
    v2set = var2set(mosaic)
    MaskedSetCoverProblem(params, mosaic.loc2glob_setix[v2set],
                          var_scores(mosaic, params, v2set),
                          varXvar_scores(mosaic, params, v2set, false),
                          varXtiles(mosaic, v2set, params)...)
end

ntiles(problem::MaskedSetCoverProblem) = length(problem.nmasked_tile)

function exclude_vars(problem::MaskedSetCoverProblem,
                      vars::AbstractVector{Int};
                      kwargs...)
    varmask, v_scores = mask_vars(problem, vars; kwargs...)
    return MaskedSetCoverProblem(problem.params,
                problem.var2set[varmask],
                v_scores, problem.varXvar_scores[varmask, varmask],
                problem.tileXvar[:, varmask], # FIXME can also remove unused tiles
                problem.maskxtileXvar[:, varmask], # FIXME can also remove unused tiles
                problem.nmasked_tile, problem.nunmasked_tile)
end

# the tuple of:
# - total number of uncovered masked elements (in all masks overlapping with the cover)
# - total number of covered unmasked elements (in all masks overlapping with the cover)
function miscover_score(w::AbstractVector{Float64}, problem::MaskedSetCoverProblem)
    __check_vars(w, problem)
    isempty(problem.nmasked_tile) && return (0.0, 0.0)
    # calculate the tile coverage weights
    wtile = fill!(acquire!(problem.arrpool, ntiles(problem)), 0.0)
    @inbounds for (varix, wvar) in enumerate(w)
        if wvar == 1.0
            wtile[view(problem.tileXvar, :, varix)] .= 1.0
        elseif wvar > 0.0
            for tileix in view(problem.tileXvar, :, varix)
                wtile[tileix] = max(wtile[tileix], wvar)
            end
        end
    end
    release!(problem.arrpool, wtile)
    wmaskxtile = acquire!(problem.arrpool, size(problem.maskxtileXvar, 1))
    res = sum(problem.nmasked_tile) - dot(problem.nmasked_tile, wtile),
          sum(maxplus_linear!(wmaskxtile, problem.maskxtileXvar, w))
    release!(problem.arrpool, wmaskxtile)
    return res
end

function var_scores(w::AbstractVector{Float64}, problem::MultiobjCoverProblem)
    __check_vars(w, problem)
    a = dot(problem.var_scores, w)
    if problem.params.setXset_factor == 0.0
        b = 0.0
    else
        # skip set interactions as it would be aggregated to zero
        setXset = fill!(acquire!(problem.arrpool, length(w)), 0.0)
        minplus_bilinear!(min, setXset, problem.varXvar_scores, w, w)
        b = -sum(setXset)
        release!(problem.arrpool, setXset)
    end
    return a, b
end

"""
    score(w::AbstractVector{Float64}, problem) -> NTuple{4, Float64}

Unfolded multiobjective score (fitness) of the OESC coverage.

* `w`: probabilities of the sets being covered

See ["Cover quality"](@ref cover_quality).
"""
function score(w::AbstractVector{Float64}, problem::MaskedSetCoverProblem)
    a, b = var_scores(w, problem)
    if problem.params.covered_factor == 0.0 &&
       problem.params.uncovered_factor == 0.0
        # skip miscover score calculation if it's excluded from aggregation
        c = 0.0
        d = 0.0
    else
        c, d = miscover_score(w, problem)
    end
    return (a, b, c, d)
end

struct WeightedSetCoverProblem <: MultiobjCoverProblem
    params::CoverParams

    var2set::Vector{Int}

    var_scores::Vector{Float64}
    varXvar_scores::Matrix{Float64}

    arrpool::ArrayPool{Float64}   # pool of vectors to reuse for evaluation FIXME move to evaluator

    function WeightedSetCoverProblem(params::CoverParams,
                        var2set::AbstractVector{Int},
                        var_scores::AbstractVector{Float64},
                        varXvar_scores::AbstractMatrix{Float64},
    )
        length(var2set) == length(var_scores) ==
            size(varXvar_scores, 1) == size(varXvar_scores, 2)
        new(params, var2set, var_scores, varXvar_scores, ArrayPool{Float64}())
    end
end

function MultiobjCoverProblem(mosaic::WeightedSetMosaic, params::CoverParams = CoverParams())
    v2set = var2set(mosaic)
    WeightedSetCoverProblem(params, mosaic.loc2glob_setix[v2set],
                            var_scores(mosaic, params, v2set),
                            varXvar_scores(mosaic, params, v2set, false))
end

function exclude_vars(problem::WeightedSetCoverProblem,
                      vars::AbstractVector{Int};
                      kwargs...)
    varmask, v_scores = mask_vars(problem, vars; kwargs...)
    return WeightedSetCoverProblem(problem.params,
                problem.var2set[varmask],
                v_scores, problem.varXvar_scores[varmask, varmask])
end

score(w::AbstractVector{Float64}, problem::WeightedSetCoverProblem) =
    var_scores(w, problem)..., 0.0, 0.0


text:
using Revise
using Distributed, MPIClusterManagers

manager = MPIManager(np=4)
addprocs(manager)

@everywhere ENV["OMP_NUM_THREADS"] = 4

@everywhere using Revise
@everywhere using Devito
@everywhere using MPI
@everywhere configuration!("log-level", "DEBUG")
@everywhere configuration!("language", "openmp")
@everywhere configuration!("mpi", true)

@everywhere function ricker(f, _t, t₀)
    t = reshape(_t, length(_t), 1)
    (1.0 .- 2.0 * (pi * f * (t .- t₀)).^2) .* exp.(-(pi * f * (t .- t₀)).^2)
end

@everywhere function model()
    write(stdout, "inside model()\n")
    grid = Grid(
        shape = (121,121,121),
        origin = (0.0,0.0,0.0),
        extent = (1210.0,1210.0,1210.0),
        dtype = Float32)

    b = Devito.Function(name="b", grid=grid, space_order=8)
    v = Devito.Function(name="vel", grid=grid, space_order=8)
    q = Devito.Function(name="woverq", grid=grid, space_order=8)

    b_data = data(b)
    v_data = data(v)
    q_data = data(q)

    b_data .= 1
    v_data .= 1.5
    q_data .= 1/1000

    time_range = 0.0f0:1.0f0:750.0f0

    p = TimeFunction(name="p", grid=grid, time_order=2, space_order=8)
    z,y,x,t = dimensions(p)

    src = SparseTimeFunction(name="src", grid=grid, f0=0.01f0, npoint=1, nt=length(time_range), coordinates=[605.0 605.0 10.0])
    src_data = data(src)
    w = ricker(0.01, collect(time_range), 125)
    copy!(src_data, w)
    src_term = inject(src; field=forward(p), expr=src * spacing(t)^2 * v^2 / b)

    nz,ny,nx,δz,δy,δx = size(grid)...,spacing(grid)...
    rec = SparseTimeFunction(name="rec", grid=grid, npoint=nx, nt=length(time_range))
    rec_coords = coordinates(rec)
    _rec_coords = zeros(Float32, length(dimensions(p))-1, nx)
    _rec_coords[1,:] .= δx*(0:nx-1)
    _rec_coords[2,:] .= 605
    _rec_coords[3,:] .= 20
    copy!(rec_coords, _rec_coords)
    rec_term = interpolate(rec, expr=p)

    g1(field) = dx(field,x0=x+spacing(x)/2)
    g2(field) = dy(field,x0=y+spacing(y)/2)
    g3(field) = dz(field,x0=z+spacing(z)/2)
    g1_tilde(field) = dx(field,x0=x-spacing(x)/2)
    g2_tilde(field) = dy(field,x0=y-spacing(y)/2)
    g3_tilde(field) = dz(field,x0=z-spacing(z)/2)

    # write the time update equation for u_x
    update_p = spacing(t)^2 * v^2 / b *
        (g1_tilde(b * g1(p)) + g2_tilde(b * g2(p)) + g3_tilde(b * g3(p))) +
        (2 - spacing(t) * q) * p +
        (spacing(t) * q - 1) * backward(p)

    stencil_p = Eq(forward(p), update_p)

    # update the dimension spacing_map to include the time dimension
    # these symbols will be replaced with the relevant scalars by the Operator
    dt = step(time_range)
    smap = spacing_map(grid)
    smap[spacing(t)] = dt

    op = Operator([stencil_p, src_term, rec_term], subs=smap, name="OpExampleIso")

    bx,by = 19,8
    t_apply = @elapsed begin
        summary = apply(op; x0_blk0_size=bx, y0_blk0_size=by)
    end

    write(stdout, "t_appy=$t_apply\n")

    _p = data_with_halo(p)
    _d = data(rec)

    __d = convert(Array, _d)

    if MPI.Comm_rank(MPI.COMM_WORLD) == 0
        write("d.bin", __d)
    end
    MPI.Barrier(MPI.COMM_WORLD)

    nothing
end

@mpi_do manager model()

d = read!("d.bin", Array{Float32}(undef,751,121))

using PyPlot
figure(); imshow(d); display(gcf())


text:
function _list_commits(repodir::String, branch::String; count = 10)
    !_check_gitdir(repodir) && return String[]
    out = _run("git -C $(repodir) rev-list $(branch) --max-count=$(count) 2>&1"; verbose = false, ignorestatus = true)
    return split(strip(out))
end


text:
# Copyright (c) 2019 Iain Dunning, and contributors                             #src
#                                                                               #src
# This work is licensed under a Creative Commons Attribution-ShareAlike 4.0     #src
# International License. URL: http://creativecommons.org/licenses/by-sa/4.0.    #src

# # Rocket Control

# **Originally Contributed by**: Iain Dunning

# This tutorial shows how to solve a nonlinear rocketry control problem.
# The problem was drawn from the [COPS3](http://www.mcs.anl.gov/~more/cops/cops3.pdf)
# benchmark.

# Our goal is to maximize the final altitude of a vertically launched rocket.

# We can control the thrust of the rocket, and must take account of the rocket
# mass, fuel consumption rate, gravity, and aerodynamic drag.

# Let us consider the basic description of the model (for the full description,
# including parameters for the rocket, see the COPS3 PDF)

# ### Overview

# We will use a discretized model of time, with a fixed number of time steps,
# $n$.

# We will make the time step size $\Delta t$, and thus the final time
# $t_f = n \cdot \Delta t$, a variable in the problem. To approximate the
# derivatives in the problem we will use the [trapezoidal rule](http://en.wikipedia.org/wiki/Trapezoidal_rule).

# ### State and Control

# We will have three state variables:
#
# * Velocity, $v$
# * Altitude, $h$
# * Mass of rocket and remaining fuel, $m$
#
# and a single control variable, thrust $T$.

# Our goal is thus to maximize $h(t_f)$.

# Each of these corresponds to a JuMP variable indexed by the time step.

# ### Dynamics

# We have three equations that control the dynamics of the rocket:
#
# Rate of ascent: $$h^\prime = v$$
# Acceleration: $$v^\prime = \frac{T - D(h,v)}{m} - g(h)$$
# Rate of mass loss: $$m^\prime = -\frac{T}{c}$$
#
# where drag $D(h,v)$ is a function of altitude and velocity, and gravity
# $g(h)$ is a function of altitude.

# These forces are defined as
#
# $$D(h,v) = D_c v^2 exp\left( -h_c \left( \frac{h-h(0)}{h(0)} \right) \right)$$
# and
# $$g(h) = g_0 \left( \frac{h(0)}{h} \right)^2$$
#
# The three rate equations correspond to JuMP constraints, and for convenience
# we will represent the forces with nonlinear expressions.

using JuMP
import Ipopt
import Plots

# Create JuMP model, using Ipopt as the solver

rocket = Model(Ipopt.Optimizer)
set_silent(rocket)

# ## Constants

# Note that all parameters in the model have been normalized
# to be dimensionless. See the COPS3 paper for more info.

h_0 = 1    # Initial height
v_0 = 0    # Initial velocity
m_0 = 1    # Initial mass
g_0 = 1    # Gravity at the surface

T_c = 3.5  # Used for thrust
h_c = 500  # Used for drag
v_c = 620  # Used for drag
m_c = 0.6  # Fraction of initial mass left at end

c     = 0.5 * sqrt(g_0 * h_0)  # Thrust-to-fuel mass
m_f   = m_c * m_0              # Final mass
D_c   = 0.5 * v_c * m_0 / g_0  # Drag scaling
T_max = T_c * g_0 * m_0        # Maximum thrust

n = 800    # Time steps

# ## Decision variables

@variables(rocket, begin
    Δt ≥ 0, (start = 1/n) # Time step
    ## State variables
    v[1:n] ≥ 0            # Velocity
    h[1:n] ≥ h_0          # Height
    m_f ≤ m[1:n] ≤ m_0    # Mass
    ## Control variables
    0 ≤ T[1:n] ≤ T_max    # Thrust
end)

# ## Objective

# The objective is to maximize altitude at end of time of flight.

@objective(rocket, Max, h[n])

# ## Initial conditions

fix(v[1], v_0; force = true)
fix(h[1], h_0; force = true)
fix(m[1], m_0; force = true)
fix(m[n], m_f; force = true)

# ## Forces

@NLexpressions(rocket, begin
    ## Drag(h,v) = Dc v^2 exp( -hc * (h - h0) / h0 )
    drag[j = 1:n], D_c * (v[j]^2) * exp(-h_c * (h[j] - h_0) / h_0)
    ## Grav(h)   = go * (h0 / h)^2
    grav[j = 1:n], g_0 * (h_0 / h[j])^2
    ## Time of flight
    t_f, Δt * n
end)

# ## Dynamics

for j in 2:n
    ## h' = v
    ## Rectangular integration
    ## @NLconstraint(rocket, h[j] == h[j - 1] + Δt * v[j - 1])
    ## Trapezoidal integration
    @NLconstraint(rocket, h[j] == h[j - 1] + 0.5 * Δt * (v[j] + v[j - 1]))
    ## v' = (T-D(h,v))/m - g(h)
    ## Rectangular integration
    ## @NLconstraint(
    ##     rocket,
    ##     v[j] == v[j - 1] + Δt *((T[j - 1] - drag[j - 1]) / m[j - 1] - grav[j - 1])
    ## )
    ## Trapezoidal integration
    @NLconstraint(
        rocket,
        v[j] == v[j-1] +
            0.5 * Δt * (
                (T[j] - drag[j] - m[j] * grav[j]) / m[j] +
                (T[j - 1] - drag[j - 1] - m[j - 1] * grav[j - 1]) / m[j - 1]
            )
    )
    ## m' = -T/c
    ## Rectangular integration
    ## @NLconstraint(rocket, m[j] == m[j - 1] - Δt * T[j - 1] / c)
    ## Trapezoidal integration
    @NLconstraint(rocket, m[j] == m[j - 1] - 0.5 * Δt * (T[j] + T[j-1]) / c)
end

# Solve for the control and state
println("Solving...")
status = optimize!(rocket)

# ## Display results

println("Max height: ", objective_value(rocket))

#-

function my_plot(y, ylabel)
    return Plots.plot(
        (1:n) * value.(Δt),
        value.(y)[:];
        xlabel = "Time (s)",
        ylabel = ylabel,
    )
end

Plots.plot(
    my_plot(h, "Altitude"),
    my_plot(m, "Mass"),
    my_plot(v, "Velocity"),
    my_plot(T, "Thrust");
    layout = (2, 2),
    legend = false,
    margin = 1Plots.cm,
)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Random
using Autumn
include("generativemodel.jl")

a = au"""(program
      (= GRID_SIZE 16)

      (object Particle (Cell 0 0 "blue"))

      (: particles (List Particle))
      (= particles
        (initnext (list)
                  (updateObj (prev particles) (--> obj (Particle (uniformChoice (adjPositions (.. obj origin))))))))

      (on clicked (= particles (addObj (prev particles) (Particle (Position (.. click x) (.. click y))))))
      )"""

mod = eval(compiletojulia(a))

rng = MersenneTwister(0)
state = mod.init(nothing, nothing, nothing, nothing, nothing, rng);

# add particles
for i in 1:3
  global state = mod.next(state, mod.Click(rand(0:15),rand(0:15)), nothing, nothing, nothing, nothing); mod.render(state.scene)
end

# let particles move
for i in 1:5
  global state = mod.next(state, mod.Click(rand(0:15),rand(0:15)), nothing, nothing, nothing, nothing); mod.render(state.scene)
end

# render
rendering = mod.render(state.scene)

# parse scene rendering into types and objects
types_and_objects = parsescene_autumn_singlecell(rendering)

# # generate random dynamics on object decomposition (types and objects)
# println(generateprogram_given_objects(types_and_objects, group=true))


text:
# Header size in bytes
const HEADER_SIZE = 512

type HeaderSizeError <: Exception end

"""
    readnimrod_hdr(fname::AbstractString)
    readnimrod_hdr(f::IOStream)::Array{Any,1}

Read the header of the Nimrod file called `fname` or from the stream
`f`. Returns an array of `Any` containing the fields of the header in the order
they are defined in according to the Nimrod file spec.
"""
readnimrod_hdr(fname::AbstractString) = open(readnimrod_hdr, fname)
function readnimrod_hdr(f::IOStream)::Array{Any,1}
    seekstart(f)

    # First 4 bytes are the length of the header in bytes. Should be 512
    if ntoh(read(f, Int32)) != HEADER_SIZE
        throw(HeaderSizeError())
    end

    # Read first 31 2 byte integer fields
    general_int_elements = map(ntoh, read(f, Int16, 31))

    # 28 4 byte Real fields
    general_real_elements = map(ntoh, read(f, Float32, 28))

    # 45 4 byte real fields
    data_real_elements = map(ntoh, read(f, Float32, 45))

    # 56 bytes of char data forming three fields
    data_char_elements = String(read(f, UInt8, 56))
    data_char_elements = [data_char_elements[1:8];
                          data_char_elements[9:32];
                          data_char_elements[33:56]];
    # Tidy null characters
    map!((s) -> strip(s, '\0'), data_char_elements)

    # Remaining 51 data specific 2 byte integer fields
    data_extra_elements = map(ntoh, read(f, Int16, 51))

    # Last 4 bytes should be size of header again.
    @assert ntoh(read(f, Int32)) == HEADER_SIZE "Last 4 bytes of header weren't its size."

    return vcat(general_int_elements,
                general_real_elements,
                data_real_elements,
                data_char_elements,
                data_extra_elements)
end

"""
    readnimrod(fname::AbstractString)
    readnimrod(f::IOStream)::Nimrod

Read the Nimrod file at path `fname` or stream `f`. Returns an instance of
`Nimrod{T}` where `T` is determined from the header of the file.
"""
readnimrod(fname::AbstractString) = open(readnimrod, fname)
function readnimrod(f::IOStream)::Nimrod
    seekstart(f)

    hdr = readnimrod_hdr(f)
    nrows = Int(hdr[16])
    ncols = Int(hdr[17])

    # Data can be integers or floats of varying sizes. Figure it out from the
    # header.
    data_format_code = hdr[12]
    data_size_code = hdr[13]
    dataType = Any

    if data_format_code == 1 && data_size_code == 1
        dataType = Int8
    elseif data_format_code == 1 && data_size_code == 2
        dataType = Int16
    elseif data_format_code == 1 && data_size_code == 4
        dataType = Int32
    elseif data_format_code == 0 && data_size_code == 2
        dataType = Float16
    elseif data_format_code == 0 && data_size_code == 4
        dataType = Float32
    else
        error("Unable to process data type specified in header. Type Code: $(data_format_code), Size: $(data_size_code)")
    end

    read(f, Int32) # Skip the data size specification
    data = map(ntoh, read(f, dataType, (ncols, nrows))')

    # Calculate the coordinates of the centre of the lower left and top left
    # resolution cells
    llcenter = (hdr[36], hdr[34] - hdr[35] * (nrows - 1))
    trcenter = (hdr[36] + (ncols - 1) * hdr[37], hdr[34])

    Nimrod(hdr, data, NimrodGridType(hdr[15]), hdr[37], hdr[35], (hdr[67], hdr[66]), llcenter, trcenter)
end


text:
mutable struct Node
	coord::Tuple{Int, Int}
	parent::Int
	id::Int
	status::UInt8
end

function find_node(nodes::Array{Node, 1}, coord::Tuple{Int, Int})::Int
	for i in 1 : length(nodes)
		if nodes[i].coord == coord
			return i
		end
	end
	return 0
end

function nodify(maze::Maze)::Tuple{Array{Node, 1}, Array{Tuple{Int, Int, Int}, 1}}
	edges::Array{Tuple{Int, Int, Int}, 1} = []
	nodes::Array{Node, 1} = []
	dirs::UInt8 = 0
	if maze.entr[1] == 1
		dirs = 1
	elseif maze.entr[2] == 1
		dirs = 2
	elseif maze.entr[1] == maze.w
		dirs = 4
	else
		dirs = 8
	end
	push!(nodes, Node(maze.entr, 0, 1, 0))
	push!(nodes, Node(maze.exit, 0, 2, 0))
	set_mazecolor(maze, maze.entr, 0xff7fff00)
	set_mazecolor(maze, maze.exit, 0xff7fff00)
	nodify(maze, 1, dirs, nodes, edges)
	return nodes, edges
end

function nodify(maze::Maze, node::Int, dirs::UInt8, nodes::Array{Node, 1}, edges::Array{Tuple{Int, Int, Int}, 1})
	x, y = nodes[node].coord
	if dirs & 1 != 0
		l = 1
		x += 1
		while get_mazecolor(maze, (x + 1, y)) == 0xffffffff && get_mazecolor(maze, (x, y + 1)) == 0xff000000 && get_mazecolor(maze, (x, y - 1)) == 0xff000000
			x += 1
			l += 1
		end
		if (get_mazecolor(maze, (x + 1, y)) == 0xff000000 || get_mazecolor(maze, (x + 1, y)) == 0xffffffff) && get_mazecolor(maze, (x, y)) == 0xffffffff
			push!(nodes, Node((x, y), node, length(nodes) + 1, 0))
			push!(edges, (node, length(nodes), l))
			set_mazecolor(maze, (x, y), 0xff7fff00)
			nodify(maze, length(nodes), get_dirs(maze, (x, y)), nodes, edges)
		else
			push!(edges, (node, find_node(nodes, (x + 1, y)), l + 1))
		end
	end
	x, y = nodes[node].coord
	if dirs & 2 != 0
		l = 1
		y += 1
		while get_mazecolor(maze, (x, y + 1)) == 0xffffffff && get_mazecolor(maze, (x + 1, y)) == 0xff000000 && get_mazecolor(maze, (x - 1, y)) == 0xff000000
			y += 1
			l += 1
		end
		if (get_mazecolor(maze, (x, y + 1)) == 0xff000000 || get_mazecolor(maze, (x, y + 1)) == 0xffffffff) && get_mazecolor(maze, (x, y)) == 0xffffffff
			push!(nodes, Node((x, y), node, length(nodes) + 1, 0))
			push!(edges, (node, length(nodes), l))
			set_mazecolor(maze, (x, y), 0xff7fff00)
			nodify(maze, length(nodes), get_dirs(maze, (x, y)), nodes, edges)
		else
			push!(edges, (node, find_node(nodes, (x, y + 1)), l + 1))
		end
	end
	x, y = nodes[node].coord
	if dirs & 4 != 0
		l = 1
		x -= 1
		while get_mazecolor(maze, (x - 1, y)) == 0xffffffff && get_mazecolor(maze, (x, y + 1)) == 0xff000000 && get_mazecolor(maze, (x, y - 1)) == 0xff000000
			x -= 1
			l += 1
		end
		if (get_mazecolor(maze, (x - 1, y)) == 0xff000000 || get_mazecolor(maze, (x - 1, y)) == 0xffffffff) && get_mazecolor(maze, (x, y)) == 0xffffffff
			push!(nodes, Node((x, y), node, length(nodes) + 1, 0))
			push!(edges, (node, length(nodes), l))
			set_mazecolor(maze, (x, y), 0xff7fff00)
			nodify(maze, length(nodes), get_dirs(maze, (x, y)), nodes, edges)
		else
			push!(edges, (node, find_node(nodes, (x - 1, y)), l + 1))
		end
	end
	x, y = nodes[node].coord
	if dirs & 8 != 0
		l = 1
		y -= 1
		while get_mazecolor(maze, (x, y - 1)) == 0xffffffff && get_mazecolor(maze, (x + 1, y)) == 0xff000000 && get_mazecolor(maze, (x - 1, y)) == 0xff000000
			y -= 1
			l += 1
		end
		if (get_mazecolor(maze, (x, y - 1)) == 0xff000000 || get_mazecolor(maze, (x, y - 1)) == 0xffffffff) && get_mazecolor(maze, (x, y)) == 0xffffffff
			push!(nodes, Node((x, y), node, length(nodes) + 1, 0))
			push!(edges, (node, length(nodes), l))
			set_mazecolor(maze, (x, y), 0xff7fff00)
			nodify(maze, length(nodes), get_dirs(maze, (x, y)), nodes, edges)
		else
			push!(edges, (node, find_node(nodes, (x, y - 1)), l + 1))
		end
	end
end

function index_of_id(nodes::Array{Node, 1}, id::Int, s::Int=0, r::Int=0)::Int
	if s == 0
		s = 1
		r = length(r)
	end
	if s == r
		return s
	end
	index = fld(s + r, 2)
	if nodes[index].id == ind
		return index
	elseif nodes[index].id > id
		return index_of_id(nodes, id, s, index - 1)
	else
		return index_of_id(nodes, id, indes + 1, r)
	end
end

function remove_dead_ends(nodes::Array{Node, 1}, edges::Array{Tuple{Int, Int, Int}, 1})
	no_dead_ends = false

	while !no_dead_ends
		loners = []
		for i in 3 : length(nodes)
			neighbours::Int = 0
			for edge in edges
				if i in edge
					neighbours += 1
				end
			end
			if neighbours <= 1
				push!(loners, i)
			end
		end
		if length(loners) == 0
			no_dead_ends = true
		else
			for i in length(loners) : -1 : 1
				deleteat!(nodes, loners[i])
			end
			i = length(edges)
			while i >= 1
				if edges[i][1] in loners || edges[i][2] in loners
					deleteat!(edges, i)
					i -= 1
				end
				i -= 1
			end

			for i in 1 : length(edges)
				edges[i] = (index_of_id(nodes, edges[i][1]), index_of_id(nodes, edges[i][2]), edges[i][3])
			end

			for i in 1 : length(nodes)
				nodes[i].id = i
			end
		end
	end
end

function mark_nodes(maze::Maze, nodes::Array{Node, 1})
	for i in nodes
		set_mazecolor(maze, i.coord, 0xff7fff00)
	end
end


text:
@testset "custom deterministic generative function" begin

    struct MyDeterministicGFState
        sum::Float64
        prev::AbstractArray
        my_param::Float64
    end

    # sum(arr) * my_param
    # arr is the argument
    # my_param is a trainable parameter

    mutable struct MyDeterministicGF <: CustomDetermGF{Float64,MyDeterministicGFState}
        my_param::Float64
        my_grad::Float64
    end

    MyDeterministicGF() = MyDeterministicGF(1., 0.)

    function Gen.execute_determ(gen_fn::MyDeterministicGF, args::Tuple)
        arr = args[1]
        retval = sum(arr) * gen_fn.my_param
        state = MyDeterministicGFState(retval, arr, gen_fn.my_param)
        retval, state
    end

    function Gen.update_determ(gen_fn::MyDeterministicGF, state, args, argdiffs)
        arr = args[1]
        retval = sum(arr) * state.my_param
        state = MyDeterministicGFState(retval, arr, state.my_param)
        state, retval, UnknownChange()
    end

    function Gen.update_determ(::MyDeterministicGF, state, args, argdiffs::Tuple{NoChange})
        state, state.sum * state.my_param, NoChange()
    end

    function Gen.update_determ(::MyDeterministicGF, state, args, argdiffs::Tuple{VectorDiff})
        arr = args[1]
        retval = state.sum * state.my_param
        for i in keys(argdiffs[1].updated)
            retval += (arr[i] - state.prev[i])
        end
        prev_length = length(state.prev)
        new_length = length(arr)
        for i=prev_length+1:new_length
            retval += arr[i]
        end
        for i=new_length+1:prev_length
            retval -= arr[i]
        end
        state = MyDeterministicGFState(retval, arr, state.my_param)
        state, retval, UnknownChange()
    end

    Gen.has_argument_grads(::MyDeterministicGF) = (true,)

    Gen.accepts_output_grad(::MyDeterministicGF) = true

    function Gen.gradient_determ(::MyDeterministicGF, state, retgrad)
        arr_gradient = fill(retgrad * state.my_param, length(state.prev))
        (arr_gradient,)
    end

    function Gen.accumulate_param_gradients_determ!(gen_fn::MyDeterministicGF, state, retgrad, scaler)
        gen_fn.my_grad += (retgrad * state.sum) * scaler
        arr_gradient = fill(retgrad * state.my_param, length(state.prev))
        (arr_gradient,)
    end

    # simulate
    gen_fn = MyDeterministicGF()
    trace = simulate(gen_fn, ([1, 2, 3],))
    @test get_retval(trace) == 1 + 2 + 3
    @test get_score(trace) == 0.
    @test project(trace, EmptySelection()) == 0.
    @test isempty(get_choices(trace))
    @test get_args(trace) == ([1, 2, 3],)
    @test get_gen_fn(trace) == gen_fn

    # generate
    trace, w = generate(MyDeterministicGF(), ([1, 2, 3],))
    @test w == 0.
    @test get_retval(trace) == 1 + 2 + 3

    # update (UnknownChange)
    trace = simulate(MyDeterministicGF(), ([1, 2, 3],))
    new_trace, w, retdiff = update(trace, ([1, 2, 4],), (UnknownChange(),), EmptyChoiceMap())
    @test w == 0.
    @test get_retval(new_trace) == 1 + 2 + 4
    @test get_args(new_trace) == ([1, 2, 4],)
    @test retdiff == UnknownChange()

    # update (VectorDiff)
    trace = simulate(MyDeterministicGF(), ([1, 2, 3],))
    diff = VectorDiff(4, 3, Dict(3 => UnknownChange()))
    new_trace, w, retdiff = update(trace, ([1, 2, 4, 5],), (diff,), EmptyChoiceMap())
    @test w == 0.
    @test get_retval(new_trace) == 1 + 2 + 4 + 5
    @test get_args(new_trace) == ([1, 2, 4, 5],)
    @test retdiff == UnknownChange()

    # regenerate
    trace = simulate(MyDeterministicGF(), ([1, 2, 3],))
    new_trace, w, retdiff = regenerate(trace, ([1, 2, 4],), (UnknownChange(),), EmptySelection())
    @test w == 0.
    @test get_retval(new_trace) == 1 + 2 + 4
    @test get_args(new_trace) == ([1, 2, 4],)
    @test retdiff == UnknownChange()

    # choice gradients
    trace = simulate(MyDeterministicGF(), ([1, 2, 3],))
    arg_grads, _, _ = choice_gradients(trace, EmptySelection(), 2.)
    @test arg_grads[1] == [2., 2., 2.]

    # accumulate parameter gradients
    gen_fn = MyDeterministicGF()
    trace = simulate(gen_fn, ([1, 2, 3],))
    arg_grads = accumulate_param_gradients!(trace, 2., 1.)
    @test arg_grads[1] == [2., 2., 2.]
    @test gen_fn.my_grad == 2. * (1 + 2 + 3)

end

max_stars_repo_name: cameronraysmith/Gen

max_issues_repo_name: cameronraysmith/Gen

max_forks_repo_name: cameronraysmith/Gen


text:
function fib50()
#
#  test polynomial suggested by Goedecker
#
    n = 50;
    p = fib(n)
    z = [-.9784087187483260;
         -.9708270453807137-.1218355210630181*im; -.9708270453807137+.1218355210630181*im;
         -.9481969430968524-.2418020920979234*im; -.9481969430968524+.2418020920979234*im;
         -.9108614275411458-.3580588620158148*im; -.9108614275411458+.3580588620158148*im;
         -.8593864256832723-.4688207429141271*im; -.8593864256832723+.4688207429141271*im;
         -.7945522211328693-.5723852109383657*im; -.7945522211328693+.5723852109383657*im;
         -.7173416644328307-.6671578269195402*im; -.7173416644328307+.6671578269195402*im;
         -.6289253327402240-.7516760773713197*im; -.6289253327402240+.7516760773713197*im;
         -.5306438727753174-.8246311589916703*im; -.5306438727753174+.8246311589916703*im;
         -.4239878085206529-.8848873569443235*im; -.4239878085206529+.8848873569443235*im;
         -.3105751412168900-.9314986982161571*im; -.3105751412168900+.9314986982161571*im;
         -.1921271147696794-.9637225955421171*im; -.1921271147696794+.9637225955421171*im;
         -.07044256696033078-.9810302342100871*im; -.07044256696033078+.9810302342100871*im;
         .5262866000569139e-1-.9831134936221445*im; .5262866000569139e-1+.9831134936221445*im;
         .1752127097133740-.9698882397909325*im; .1752127097133740+.9698882397909325*im;
         .2954384705313263-.9414938807153827*im; .2954384705313263+.9414938807153827*im;
         .4114633335560060-.8982891618708299*im; .4114633335560060+.8982891618708299*im;
         .5214969674899803-.8408443381786473*im; .5214969674899803+.8408443381786473*im;
         .6238218132212410-.7699301944269942*im; .6238218132212410+.7699301944269942*im;
         .7168085997048416-.6865051209865677*im; .7168085997048416+.6865051209865677*im;
         .7989249239308063-.5917030159705396*im; .7989249239308063+.5917030159705396*im;
         .8687357356619744-.4868277941747678*im; .8687357356619744+.4868277941747678*im;
         .9248989976696934-.3733647458567231*im; .9248989976696934+.3733647458567231*im;
         .9661725191057798-.2530209270696782*im; .9661725191057798+.2530209270696782*im;
         .9914691930342274-.1277916638512058*im; .9914691930342274+.1277916638512058*im;
         1.999999999999999];

    z = [z ones(n)]
    p, PolyZeros(z)
end


text:
using DataFrames
using CSV
using SpatialEcology
using Test
using Plots

@testset "Analysis" begin
    amphdat = CSV.read(joinpath(dirname(pathof(SpatialEcology)), "..", "data", "amph_Europe.csv"), DataFrame)
    amph = Assemblage(amphdat[!,4:end], amphdat[!,1:3], sitecolumns = false)
    addtraits!(amph, asquantiles(occupancy(amph), 4), :quantile)
    gb = groupspecies(amph, :quantile) # error ?
    ps = [plot(g) for g in gb]
    p = ps[4];

    @test p[1][1][:seriestype] == :heatmap
    @test p[1][1][:z].surf[30,30] == 5.0
end


text:
using Documenter, NameResolution

makedocs(;
    modules=[NameResolution],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/thautwarm/NameResolution.jl/blob/{commit}{path}#L{line}",
    sitename="NameResolution.jl",
    authors="thautwarm",
    assets=String[],
)

deploydocs(;
    repo="github.com/thautwarm/NameResolution.jl",
)


text:
# derive epsilon given K=minpts using the median
function eps_Kmedian{T <: Real}(A::Matrix{T},k::Int)
  return median(Knearest(A,k))
end

# derive epsilon given K=minpts using the mean
function eps_Kmean{T <: Real}(A::Matrix{T},k::Int)
  return mean(Knearest(A,k))
end

#derive epsilon from unit grid of region
function eps_grid{T <: Real}(A::Matrix{T})
  gd = maximum(A)/sqrt(2)
  n = size(A,1)
  return gd/sqrt(n)
end

#derive epsilon from random [median, .85%]
function eps_randm85{T <: Real}(A::Matrix{T},k::Int)
  i = rand(0.5:0.01:0.85)
  return quantile(Knearest(A,k),i)
end

#derivce epsilon from quanitle at K level
function eps_Q{T <: Real}(A::Matrix{T},Q::Float64,k::Int)
  @assert 0.0 <= Q <= 1.0 "Quantile Q must be in range [0.0, 1.0]."
  return quantile(Knearest(A,k),Q)
end

function eps_Q{T <: Real}(A::Matrix{T},Q::Vector{Float64},k::Int)
  return quantile(Knearest(A,k),Q)
end

#= Strip all zero values and return a vector
    of all remaining entities =#
function nonZeros{T <: Real}(A::Matrix{T})
  nnzA = countnz(A)
  NZs = Vector{T}(nnzA)
  count = 1
  if nnzA > 0
    for j=indices(A,2), i=indices(A,1)
      Aij = A[i,j]
      if Aij != zero(T)
        NZs[count] = Aij
        count += 1
      end
    end
  end
  return NZs
end

function Knearest{T <: Real}(A::Matrix{T},k::Int)
  n = size(A,2)
  mins = zeros(T,n)
  for i=1:n
    a = A[:,i]
    mins[i] = a[sortperm(a)[k+1]]
  end
  return sort(mins)
end


text:
module testSep_Constraints
# sep_rank_proj_path = dirname(dirname(@__FILE__))
# Pkg.activate(sep_rank_proj_path)
using Test
using JuMP

srcDir = pwd()*"\\src\\"
include(srcDir*"Utils_states.jl")
include(srcDir*"Examples.jl")
include(srcDir*"Moments.jl")
include(srcDir*"R_constraints.jl")

using .Utils_states
using .Examples
using .Moments
using .R_constraints

examples = Examples.get_examples()

@testset "idx2var" begin
    t = (2,2)
    d = (2,2)
    n = sum(d)
    model = Model()
    mom_list = R_constraints.make_mon_expo_keys(n,t[1]) # Define variables in the moment matrix.
    @variable(model, Lx[mom_list] ) #
    @test length(Lx) == 70

    M_full = Moments.make_mon_expo(4,t;isle = true)
    M      = M_full[1:4,1:4]
    A      = reshape(1:16,4,4)
    HC =    [Lx[[0, 0, 0, 0]]  Lx[[1, 0, 0, 0]]  Lx[[0, 1, 0, 0]]  Lx[[0, 0, 1, 0]]
             Lx[[1, 0, 0, 0]]  Lx[[2, 0, 0, 0]]  Lx[[1, 1, 0, 0]]  Lx[[1, 0, 1, 0]]
             Lx[[0, 1, 0, 0]]  Lx[[1, 1, 0, 0]]  Lx[[0, 2, 0, 0]]  Lx[[0, 1, 1, 0]]
             Lx[[0, 0, 1, 0]]  Lx[[1, 0, 1, 0]]  Lx[[0, 1, 1, 0]]  Lx[[0, 0, 2, 0]]]
    @test Utils_cons.idx2var(Lx, M) ==  HC
    @test Utils_cons.idx2var(Lx, M, A) == A .* HC
    M_dict = Dict(1 => M,
                  2 => M_full[1:4,5:8],
                  3 => M_full[5:8,1:4])
    A_dict = Dict(1 => A,
                  2 => 1 ./ A,
                  3 => -A )
    @test Utils_cons.idx2var(Lx, M_dict, A_dict) == Utils_cons.idx2var(Lx,M,A) +
                                              Utils_cons.idx2var(Lx,M_full[1:4,5:8],1 ./ A) +
                                              Utils_cons.idx2var(Lx,M_full[5:8,1:4],-A)
end

## Separable states
@testset "Randd₁4d₂4" begin
    ex = "Randd₁4d₂4"
    d = Utils_states.extr_d(ex)
    @test d == (4,4)
    ρ = examples["sep"][ex]

    n = sum(d)
    model = Model()
    t = (3,3)
    list_of_keys = Moments.make_mom_expo_keys(n,t) # Define variables in the moment matrix.
    @variable(model, Lx[list_of_keys] ) #
    @test length(sep_Constraints.R_cons.make_loc_cons_S₁(ρ,d,t,Lx)) == 32
    @test length(sep_Constraints.R_cons.make_loc_cons_S₂(ρ,d,t,Lx)) == 8
    @test length(sep_Constraints.R_cons.make_loc_cons_S₃(ρ,d,t,Lx)) == 2
    @test length(sep_Constraints.R_cons.make_weakG_con(ρ,d,t,Lx)) == 4
    @test length(sep_Constraints.R_cons.make_G_con(ρ,d,t,Lx)) == 4
end

@testset "HHH1d₁2d₂2" begin
    ex = "HHH1d₁2d₂2"
    d = Utils_states.extr_d(ex)
    @test d == (2,2)
    ρ = examples["sep"][ex]

    n = sum(d)
    model = Model()
    t = (2,2)
    list_of_keys = Moments.make_mom_expo_keys(n,t) # Define variables in the moment matrix.
    @variable(model, Lx[list_of_keys] ) #
    @test length(sep_Constraints.R_cons.make_loc_cons_S₁(ρ,d,t,Lx)) == 16
    @test length(sep_Constraints.R_cons.make_loc_cons_S₂(ρ,d,t,Lx)) == 8
    @test length(sep_Constraints.R_cons.make_loc_cons_S₃(ρ,d,t,Lx)) == 2
    @test length(sep_Constraints.R_cons.make_weakG_con(ρ,d,t,Lx)) == 0
    @test length(sep_Constraints.R_cons.make_G_con(ρ,d,t,Lx)) == 4
end

@testset "DNY1d₁3d₂3" begin
    ex = "DNY1d₁3d₂3"
    d = Utils_states.extr_d(ex)
    @test d == (3,3)
    ρ = examples["sep"][ex]

    n = sum(d)
    model = Model()
    t = (4,4)
    list_of_keys = Moments.make_mom_expo_keys(n,t) # Define variables in the moment matrix.
    @variable(model, Lx[list_of_keys] ) #
    @test length(sep_Constraints.R_cons.make_loc_cons_S₁(ρ,d,t,Lx)) == 24
    @test length(sep_Constraints.R_cons.make_loc_cons_S₂(ρ,d,t,Lx)) == 8
    @test length(sep_Constraints.R_cons.make_loc_cons_S₃(ρ,d,t,Lx)) == 2
    @test length(sep_Constraints.R_cons.make_weakG_con(ρ,d,t,Lx)) == 8
    @test length(sep_Constraints.R_cons.make_G_con(ρ,d,t,Lx)) == 4
end

end

# ρ_dict = Examples.get_examples()
# ρ_sep  = ρ_dict["sep"]
# ρ = ρ_sep["sep3d3r3"]
#
# d =  Int(sqrt(size(ρ)[1]))
# n = 2*d
# t = 3
#
# @testset "Model setup" begin
#     model = Model()
# #end
#
# #@testset "Making variables" begin
#     list_of_keys = Moments.make_mom_expo_keys(n,t) # Define variables in the moment matrix.
#     @variable(model, Lx[list_of_keys] ) # Create the variables, i.e., entries of moment matrix.
#
#     @test size(Lx) == size(list_of_keys)
#     @test supertype(typeof(Lx)) == AbstractVector{VariableRef}
# #end
#
# #@testset "Build the moment matrix and constrain it to be PSD  : L([x,y]≦ₜᵀ[x,y]≦ₜ) ⪰ 0" begin
#     mom_matₜ_expo = Moments.make_mon_expo_mat(n,t,ρ,true)
#     mom_matₜ      = Utils.index_to_var(Lx, mom_matₜ_expo)
#
#     mom_matₜ₋₁_expo = Moments.make_mon_expo_mat(n,t-1,true)
#
#     @test size(mom_matₜ) == size(mom_matₜ_expo)
#     @test supertype(typeof(mom_matₜ)) ==  DenseArray{VariableRef,2}
# #end
#
# #@testset "Fourth order Moment constraints: L(xxᵀ ⊗ yyᵀ) = ρ," begin
#     xxᵀ_tens_yyᵀ    = Moments.make_xxᵀ_tens_yyᵀ(d,ρ)
#     L_xxᵀ_tens_yyᵀ  = Utils.index_to_var(Lx,xxᵀ_tens_yyᵀ)
#
#     @test typeof(L_xxᵀ_tens_yyᵀ) == Matrix{VariableRef}
#     @test size(L_xxᵀ_tens_yyᵀ) == size(Moments.prop_zero_diags(ρ))
# #end
#
# #@testset "Localizing g constraints: S₁" begin
#     loc_con_S₁ = make_loc_cons_S₁(ρ,t,d,Lx)
#     @test length(loc_con_S₁) == n
#     rand_key = [key for key in keys(loc_con_S₁)][rand(1:n)]
#     @test typeof(loc_con_S₁[rand_key]) == Matrix{GenericAffExpr{Float64,VariableRef}}
#     @test size(loc_con_S₁[rand_key]) == size(mom_matₜ₋₁_expo)
# #end
#
# #@testset "Localizing g constraints: S₂" begin
#     loc_con_S₂ = make_loc_cons_S₂(ρ,t,d,Lx)
#     @test length(loc_con_S₂) == 2
#     rand_key = [key for key in keys(loc_con_S₂)][rand(1:2)]
#     @test typeof(loc_con_S₂[rand_key]) == Matrix{GenericAffExpr{Float64,VariableRef}}
#     @test size(loc_con_S₂[rand_key]) == size(mom_matₜ₋₁_expo)
# #end
#
# #@testset "Localizing g constraints: S₃" begin
#     loc_con_S₃, g₂ = make_loc_cons_S₃(ρ,t,d,Lx)
#     @test length(loc_con_S₃) == 1
#     @test typeof(loc_con_S₃[1]) == Matrix{GenericAffExpr{Float64,VariableRef}}
#     @test size(loc_con_S₃[1]) == size(mom_matₜ₋₁_expo)
# #end
#
# #@testset "weak G Constraints" begin
#     weakG_con = make_weakG_con(ρ,t,d,Lx)
#     for key in keys(weakG_con)
#       @test size(weakG_con[key]) == size(make_mon_expo_mat(n,key,false)).*size(ρ)[1]
#       @test typeof(weakG_con[1]) == Matrix{GenericAffExpr{Float64,VariableRef}}
#     end
# #end
#
# #@testset "G Constraints" begin
#     G_con = make_G_con(ρ,t,d,Lx)
#     @test typeof(G_con) == Matrix{GenericAffExpr{Float64,VariableRef}}
#     @test size(G_con ) == size(make_mon_expo_mat(n,t-2,true)).*size(ρ)[1]
# end
#
#
#
#
# end

max_stars_repo_name: JAndriesJ/sep-rank-bounding

max_issues_repo_name: JAndriesJ/sep-rank-bounding

max_forks_repo_name: JAndriesJ/sep-rank-bounding