text:
# Tests for Links
l1 = Link("a to b", Intervals([(0.33,0.66)]), "a", "b", .5)
l2 = Link("a to b", Intervals([(0.33,0.66)]), "a", "b", .5)
l3 = Link("a to c", Intervals([(0.45,0.56)]), "a", "c", .5)

@test l1 == l2
@test l2 == l1
@test from_match(l1,l2)
@test from_match("a",l1)
@test from_match("a",l2)
@test from_match(l1,l3)
@test from_match("a",l3)
@test to_match(l1,l2)
@test to_match("b",l1)
@test to_match("b",l2)
@test to_match("c",l3)
@test !to_match(l1,l3)
@test !to_match("c",l1)
@test !to_match("a",l3)
@test match("a","b",l1)
@test match("a","b",l2)
@test match("a","c",l3)
@test !match("a","b",l3)
@test 0.5 ∈ l1
@test 0.67 ∉ l2
@test 0.51 ∈ l3
@test l3 ⊆ l1
@test l1 ⊈ l3
@test (l1 ∩ l3).list == [(0.45,0.56)]
@test (l3 ∪ l1).list == [(0.33,0.66)]

l1 = Link("a to b", Intervals([(0.33,0.66)]), "a", "b", .5)
l2 = Link("a to c", Intervals([(0.45,0.56)]), "a", "c", .5)
l3 = Link("b to c", Intervals([(0.5,1.0)]), "b", "c", .5)
L1 = [l3,l2]
@test from_match(l1,L1) == [2]
@test from_match("a",L1) == [2]
@test to_match(l2,L1) == [1,2]
@test to_match("c",L1) == [1,2]
@test match(l3,L1) == [1]
@test match("b","c",L1) == [1]
@test get_idx(l2,L1) == [2]
@test get_idx("b to c",L1) == [1]
@test l2 ∈ L1
@test l3 ∈ L1
@test l1 ∉ L1
L1 = [l1,l2,l3]
@test get_idx(0.4,L1) == [1]
@test get_idx(0.55,L1) == [1,2,3]
@test get_idx(1.0001,L1) == []
L2 = [l1,l3]
@test L1 ∩ L2 == [l1,l3]
@test L1 ∪ L2 == [l1,l2,l3]

l1 = Link("a to b", Intervals([(0.33,0.66),(1.33,1.66)]), "a", "b", .5)
l2 = Link("a to b", Intervals([(0.2,1.1),(1.2,1.3),(1.5,1.9)]), "a", "b", .5)
l3 = merge(l1,l2)
@test l3.name == "a to b"
@test l3.from == "a"
@test l3.to == "b"
@test l3.weight == .5
@test l3.presence == Intervals([(0.2,1.1),(1.2,1.3),(1.33,1.9)])
merge!(l1,l2)
@test l1.name == "a to b"
@test l1.from == "a"
@test l1.to == "b"
@test l1.weight == .5
@test l1.presence == Intervals([(0.2,1.1),(1.2,1.3),(1.33,1.9)])


text:
with(PNG("red-blue-square.png", 100, 100), Ink("red"), Pen(0.1)) do
    draw(Ink("blue")) do
        move(0, 0)
        line(1, 0)
        line(1, 1)
    end
    draw() do
        line(0, 1)
        line(0, 0)
    end
end

ext: jl


text:
function solve_callback!(x, opt::SS.Optimizer, callback)
    SS.initialize!(opt, x)

    while !SS.meets_stopping_criteria(opt)
        SS.next_iteration!(result(opt))
        SS.solver_step!(x, state(opt))
        SS.update!(opt, x)

        callback(opt)
    end

    SS.warn_iteration_number(status(opt), config(opt))

    return x
end

max_stars_repo_name: JuliaRCM/WassersteinDictionaries

max_issues_repo_name: JuliaRCM/WassersteinDictionaries

max_forks_repo_name: JuliaRCM/WassersteinDictionaries


text:
"""
    pdfsquaredL2norm(d::Distribution)

Return the square of the L2 norm of the probability density function ``f(x)`` of the distribution `d`:

```math
\\int_{S} f(x)^{2} \\mathrm{d} x
```

where ``S`` is the support of ``f(x)``.
"""
pdfsquaredL2norm

function pdfsquaredL2norm(d::Beta)
    α, β = params(d)
    z = beta(2 * α - 1, 2 * β - 1) / beta(α, β) ^ 2
    # L2 norm of the pdf converges only for α > 0.5 and β > 0.5
    return α > 0.5 && β > 0.5 ? z : oftype(z, Inf)
end

pdfsquaredL2norm(d::Cauchy) = inv2π / d.σ

function pdfsquaredL2norm(d::Chi)
    ν = d.ν
    z = (2 ^ (1 - ν) * gamma((2 * ν - 1) / 2)) / gamma(ν / 2) ^ 2
    # L2 norm of the pdf converges only for ν > 0.5
    return d.ν > 0.5 ? z : oftype(z, Inf)
end

function pdfsquaredL2norm(d::Chisq)
    ν = d.ν
    z = gamma(d.ν - 1) / (gamma(d.ν / 2) ^ 2 * 2 ^ d.ν)
    # L2 norm of the pdf converges only for ν > 1
    return ν > 1 ? z : oftype(z, Inf)
end

pdfsquaredL2norm(d::Exponential) = 1 / (2 * d.θ)

function pdfsquaredL2norm(d::Gamma)
    α, θ = params(d)
    z = (2^(1 - 2 * α) * gamma(2 * α - 1)) / (gamma(α) ^ 2 * θ)
    # L2 norm of the pdf converges only for α > 0.5
    return α > 0.5 ? z : oftype(z, Inf)
end

pdfsquaredL2norm(d::Logistic) = 1 / (6 * d.θ)

pdfsquaredL2norm(d::Normal) = inv(sqrt4π * d.σ)

pdfsquaredL2norm(d::Uniform) = 1 / (d.b - d.a)


text:
# set up =======================================================================

# set packages
using DelimitedFiles

include("../01_species_list_creation.jl")

# set all positional values
organism_id_pos = 1
organism_type_pos = 2
body_size_pos = 3
lower_prey_limit_pos = 4
upper_prey_limit_pos = 5
habitat_midpoint_pos = 6
habitat_lower_limit_pos = 7
habitat_upper_limit_pos = 8
nutritional_value_pos = 9
biomass_pos = 10
network_id_pos = 11
alive_pos = 12

# the different values that have different options are:
#   richness shape
#   biomass shape

i = 1
j = 1
network_id = 1

all_network_data = Array{Any}(undef, (450000, 6))
all_networks = Array{Int}(undef, (450000, 51))
all_species_lists = Array{Any}(undef, (450000, 11))

for i in 1:reps
    for richness in ["top-heavy", "bottom-heavy", "uniform"]
        for biomass in ["pyramid", "uniform", "inverted"]
            #print("starting richness = ", richness, ", biomass = ", biomass, ". ")
            temp = Array{Any}(undef, (50, 6))
            # get species list and network
            species_list, network = network_draw(50, richness)
            #print("sp list richness = ", richness, ", biomass = ", biomass, ". ")
            # save network
            all_networks[j:j+49, 1] = repeat([network_id], 50)
            all_networks[j:j+49, 2:51] = network
            # get biomass
            species_list = initialize_biomass(species_list, biomass)
            #print("sp list biomass richness = ", richness, ", biomass = ", biomass, ". ")
            # append the network id to each species list item
            species_list = hcat(species_list, fill(i, 50))
            # put the species list in the all_species_lists
            all_species_lists[j:j+49, :] = species_list
            # put each species with the network data in all_network_data
            for k in 1:size(species_list,1)
                temp[k, 1] = network_id
                temp[k, 2] = richness
                temp[k, 3] = biomass
                temp[k, 4] = species_list[k, body_size_pos]
                temp[k, 5] = species_list[k, biomass_pos]
                temp[k, 6] = species_list[k, organism_type_pos]
            end
            all_network_data[(j:j+49), :] = temp
            #print("fill network data richness = ", richness, ", biomass = ", biomass, ". ")
            # increment
            j = j + 50
            network_id = network_id + 1
        end
    end
    if (i % 5 == 0) print(i, " out of 1000 --") end
end

species_list, network = network_draw(50, richness)
for richness in ["top-heavy", "bottom-heavy", "uniform"]
    print("richness = ", richness, ". ")
    for biomass in ["pyramid", "uniform", "inverted"]
        print("biomass = ", biomass, ". ")
        species_list = initialize_biomass(species_list, biomass)
    end
end

richness = "top-heavy"
for biomass in ["pyramid", "uniform", "inverted"]
    print("biomass = ", biomass, ". ")
    species_list, network = network_draw(50, richness)
    species_list = initialize_biomass(species_list, biomass)
end
function make_species_lists(
    i::Int,
    j::Int,
    network_id::Int,
    richness::String,
    biomass::String
    )

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text: print("Hello World 👋")

max_stars_repo_name: fadhilsaheer/hello-world

max_issues_repo_name: fadhilsaheer/hello-world

max_forks_repo_name: fadhilsaheer/hello-world


text:
using Base58: base58checkdecode

@testset "Transaction" begin
    raw_tx = hex2bytes("0100000001813f79011acb80925dfe69b3def355fe914bd1d96a3f5f71bf8303c6a989c7d1000000006b483045022100ed81ff192e75a3fd2304004dcadb746fa5e24c5031ccfcf21320b0277457c98f02207a986d955c6e0cb35d446a89d3f56100f4d7f67801c31967743a9c8e10615bed01210349fc4e631e3624a545de3f89f5d8684c7b8138bd94bdd531d2e213bf016b278afeffffff02a135ef01000000001976a914bc3b654dca7e56b04dca18f2566cdaf02e8d9ada88ac99c39800000000001976a9141c4bc762dd5423e332166702cb75f40df79fea1288ac19430600")
    @testset "Parsing" begin
        @testset "Version" begin
            stream = IOBuffer(raw_tx)
            tx = txparse(stream)
            @test tx.version == 1
        end
        @testset "Inputs" begin
            stream = IOBuffer(raw_tx)
            tx = txparse(stream)
            @test length(tx.tx_ins) == 1
            want = hex2bytes("d1c789a9c60383bf715f3f6ad9d14b91fe55f3deb369fe5d9280cb1a01793f81")
            @test tx.tx_ins[1].prev_tx == want
            @test tx.tx_ins[1].prev_index == 0
            want = hex2bytes("6b483045022100ed81ff192e75a3fd2304004dcadb746fa5e24c5031ccfcf21320b0277457c98f02207a986d955c6e0cb35d446a89d3f56100f4d7f67801c31967743a9c8e10615bed01210349fc4e631e3624a545de3f89f5d8684c7b8138bd94bdd531d2e213bf016b278a")
            @test Bitcoin.scriptserialize(tx.tx_ins[1].script_sig) == want
            @test tx.tx_ins[1].sequence == 0xfffffffe
        end
        @testset "Outputs" begin
            stream = IOBuffer(raw_tx)
            tx = txparse(stream)
            @test length(tx.tx_outs) == 2
            want = 32454049
            @test tx.tx_outs[1].amount == want
            want = hex2bytes("1976a914bc3b654dca7e56b04dca18f2566cdaf02e8d9ada88ac")
            @test Bitcoin.scriptserialize(tx.tx_outs[1].script_pubkey) == want
            want = 10011545
            @test tx.tx_outs[2].amount == want
            want = hex2bytes("1976a9141c4bc762dd5423e332166702cb75f40df79fea1288ac")
            @test Bitcoin.scriptserialize(tx.tx_outs[2].script_pubkey) == want
        end
        @testset "Locktime" begin
            stream = IOBuffer(raw_tx)
            tx = txparse(stream)
            @test tx.locktime == 410393
        end
    end
    @testset "Serialize" begin
        stream = IOBuffer(raw_tx)
        tx = txparse(stream)
        @test txserialize(tx) == raw_tx
    end
    @testset  "Input Value" begin
        tx_hash = "d1c789a9c60383bf715f3f6ad9d14b91fe55f3deb369fe5d9280cb1a01793f81"
        index = 0
        want = 42505594
        tx_in = Bitcoin.TxIn(hex2bytes(tx_hash), index)
        @test Bitcoin.txinvalue(tx_in) == want
    end
    @testset "Input PubKey" begin
        tx_hash = "d1c789a9c60383bf715f3f6ad9d14b91fe55f3deb369fe5d9280cb1a01793f81"
        index = 0
        want = hex2bytes("1976a914a802fc56c704ce87c42d7c92eb75e7896bdc41ae88ac")
        tx_in = Bitcoin.TxIn(hex2bytes(tx_hash), index)
        @test Bitcoin.scriptserialize(Bitcoin.txin_scriptpubkey(tx_in)) == want
    end
    @testset "Fee" begin
        raw_tx = hex2bytes("0100000001813f79011acb80925dfe69b3def355fe914bd1d96a3f5f71bf8303c6a989c7d1000000006b483045022100ed81ff192e75a3fd2304004dcadb746fa5e24c5031ccfcf21320b0277457c98f02207a986d955c6e0cb35d446a89d3f56100f4d7f67801c31967743a9c8e10615bed01210349fc4e631e3624a545de3f89f5d8684c7b8138bd94bdd531d2e213bf016b278afeffffff02a135ef01000000001976a914bc3b654dca7e56b04dca18f2566cdaf02e8d9ada88ac99c39800000000001976a9141c4bc762dd5423e332166702cb75f40df79fea1288ac19430600")
        stream = IOBuffer(raw_tx)
        tx = txparse(stream)
        @test Bitcoin.txfee(tx) == 40000
        raw_tx = hex2bytes("010000000456919960ac691763688d3d3bcea9ad6ecaf875df5339e148a1fc61c6ed7a069e010000006a47304402204585bcdef85e6b1c6af5c2669d4830ff86e42dd205c0e089bc2a821657e951c002201024a10366077f87d6bce1f7100ad8cfa8a064b39d4e8fe4ea13a7b71aa8180f012102f0da57e85eec2934a82a585ea337ce2f4998b50ae699dd79f5880e253dafafb7feffffffeb8f51f4038dc17e6313cf831d4f02281c2a468bde0fafd37f1bf882729e7fd3000000006a47304402207899531a52d59a6de200179928ca900254a36b8dff8bb75f5f5d71b1cdc26125022008b422690b8461cb52c3cc30330b23d574351872b7c361e9aae3649071c1a7160121035d5c93d9ac96881f19ba1f686f15f009ded7c62efe85a872e6a19b43c15a2937feffffff567bf40595119d1bb8a3037c356efd56170b64cbcc160fb028fa10704b45d775000000006a47304402204c7c7818424c7f7911da6cddc59655a70af1cb5eaf17c69dadbfc74ffa0b662f02207599e08bc8023693ad4e9527dc42c34210f7a7d1d1ddfc8492b654a11e7620a0012102158b46fbdff65d0172b7989aec8850aa0dae49abfb84c81ae6e5b251a58ace5cfeffffffd63a5e6c16e620f86f375925b21cabaf736c779f88fd04dcad51d26690f7f345010000006a47304402200633ea0d3314bea0d95b3cd8dadb2ef79ea8331ffe1e61f762c0f6daea0fabde022029f23b3e9c30f080446150b23852028751635dcee2be669c2a1686a4b5edf304012103ffd6f4a67e94aba353a00882e563ff2722eb4cff0ad6006e86ee20dfe7520d55feffffff0251430f00000000001976a914ab0c0b2e98b1ab6dbf67d4750b0a56244948a87988ac005a6202000000001976a9143c82d7df364eb6c75be8c80df2b3eda8db57397088ac46430600")
        stream = IOBuffer(raw_tx)
        tx = txparse(stream)
        @test Bitcoin.txfee(tx) == 140500
    end
    @testset "Sig Hash" begin
        raw_tx = hex2bytes("0100000001813f79011acb80925dfe69b3def355fe914bd1d96a3f5f71bf8303c6a989c7d1000000006b483045022100ed81ff192e75a3fd2304004dcadb746fa5e24c5031ccfcf21320b0277457c98f02207a986d955c6e0cb35d446a89d3f56100f4d7f67801c31967743a9c8e10615bed01210349fc4e631e3624a545de3f89f5d8684c7b8138bd94bdd531d2e213bf016b278afeffffff02a135ef01000000001976a914bc3b654dca7e56b04dca18f2566cdaf02e8d9ada88ac99c39800000000001976a9141c4bc762dd5423e332166702cb75f40df79fea1288ac19430600")
        stream = IOBuffer(raw_tx)
        tx = txparse(stream)
        want = parse(BigInt, "27e0c5994dec7824e56dec6b2fcb342eb7cdb0d0957c2fce9882f715e85d81a6", base=16)
        @test Bitcoin.txsighash(tx, 0) == want
    end
    @testset "Verify P2PKH" begin
        tx = txfetch("452c629d67e41baec3ac6f04fe744b4b9617f8f859c63b3002f8684e7a4fee03")
        @test txverify(tx)
        tx = txfetch("5418099cc755cb9dd3ebc6cf1a7888ad53a1a3beb5a025bce89eb1bf7f1650a2", true)
        @test txverify(tx)
    end
    @testset "Sign Input" begin
        private_key = PrivateKey(8675309)
        tx_ins = TxIn[]
        prev_tx = hex2bytes("0025bc3c0fa8b7eb55b9437fdbd016870d18e0df0ace7bc9864efc38414147c8")
        push!(tx_ins, TxIn(prev_tx, 0))
        tx_outs = TxOut[]
        h160 = base58checkdecode(b"mzx5YhAH9kNHtcN481u6WkjeHjYtVeKVh2")[2:end]
        push!(tx_outs, TxOut(Int(0.99 * 100000000), Bitcoin.p2pkh_script(h160)))
        h160 = base58checkdecode(b"mnrVtF8DWjMu839VW3rBfgYaAfKk8983Xf")[2:end]
        push!(tx_outs, TxOut(Int(0.1 * 100000000), Bitcoin.p2pkh_script(h160)))
        tx = Tx(1, tx_ins, tx_outs, 0, true)
        @test txsigninput(tx, 0, private_key)
    end
end


text:
#Função de substituição reversa, onde retorna um x que é o o que queremos encontrar de Ax=b

function sr00(U,b);

	(m,n) = size(U);
	x = zeros(m);

	x[m] = b[m]/U[m,m];

	for i = m-1:-1:1
		soma = 0;

		for j = i+1:n

			soma = soma + U[i,j]*x[j];
		end
		x[i] = (b[i] - soma)/U[i,i];
	end

	return x;

end


text:
using GLVisualize, GeometryTypes, GLAbstraction, ModernGL, FileIO, Reactive

if !isdefined(:runtests)
    window = glscreen()
    timesignal = loop(linspace(0f0, 1f0, 360))
end

description = """
Showing bilboard particles, which always face the camera. Also showing how
to use Images as a particle primitive.
"""

rotation_angle = const_lift(*, timesignal, 2f0*pi)
rotation = map(rotationmatrix_z, rotation_angle)

const b = Point3f0[(rand(Point3f0)*2)-1 for i=1:64]

sprites = visualize(
    (Circle(Point2f0(0), 0.25f0), b),
    billboard=true, image=loadasset("foxy.png"),
    model=rotation
)

_view(sprites, window, camera=:perspective)

if !isdefined(:runtests)
    renderloop(window)
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -6.355742389270071 <= q <= 5.291892702058364)

                     add_NL_constraint(m, :(swish(-0.046987465146302476 + 0.9722813338228051*$(x[1]) + -0.7814511744189074*$(x[2])) + swish(-0.9625414337890552 + 0.5872522747234861*$(x[1]) + 0.38591887418496684*$(x[2])) + swish(0.22754410958881088 + -0.15392955055702906*$(x[1]) + -0.6389708691593214*$(x[2])) + swish(-0.10351168983583303 + 0.8921688569394703*$(x[1]) + -0.5909582391668393*$(x[2])) + swish(0.35357163557652704 + -0.6460660745152582*$(x[1]) + 0.17482029817613265*$(x[2])) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function cuda_geometry(:: ActivationFunction, output :: Blob)
  len = length(output)

  x_block = round(Int, ceil(convert(Float64, len)/1024));
  return ((x_block,1024), (len,))
end

############################################################
# Rectified-Linear
############################################################
function forward(backend :: GPUBackend, neuron :: Neurons.ReLU, output :: Blob)
  cuda_dim, blob_dim = cuda_geometry(neuron, output)
  data_type = eltype(output)
  if data_type == Float32
    kernel = backend.mocha.relu_forward_float
  elseif data_type == Float64
    kernel = backend.mocha.relu_forward_double
  else
    error("Unsupported data type $data_type")
  end
  CUDA.launch(kernel, cuda_dim..., tuple(convert(data_type, neuron.epsilon), output.ptr.p, blob_dim...))
end

function backward(backend :: GPUBackend, neuron :: Neurons.ReLU, output :: Blob, gradient :: Blob)
  cuda_dim, blob_dim = cuda_geometry(neuron, output)
  data_type = eltype(output)
  if data_type == Float32
    kernel = backend.mocha.relu_backward_float
  elseif data_type == Float64
    kernel = backend.mocha.relu_backward_double
  else
    error("Unsupported data type $data_type")
  end
  CUDA.launch(kernel, cuda_dim..., tuple(convert(data_type, neuron.epsilon), output.ptr.p, gradient.ptr.p, blob_dim...))
end

############################################################
# Leaky Rectified-Linear
############################################################
function forward(backend :: GPUBackend, neuron :: Neurons.LReLU, output :: Blob)
  cuda_dim, blob_dim = cuda_geometry(neuron, output)
  data_type = eltype(output)
  if data_type == Float32
    kernel = backend.mocha.lrelu_forward_float
  elseif data_type == Float64
    kernel = backend.mocha.lrelu_forward_double
  else
    error("Unsupported data type $data_type")
  end
  CUDA.launch(kernel, cuda_dim..., tuple(output.ptr.p, blob_dim...))
end

function backward(backend :: GPUBackend, neuron :: Neurons.LReLU, output :: Blob, gradient :: Blob)
  cuda_dim, blob_dim = cuda_geometry(neuron, output)
  data_type = eltype(output)
  if data_type == Float32
    kernel = backend.mocha.lrelu_backward_float
  elseif data_type == Float64
    kernel = backend.mocha.lrelu_backward_double
  else
    error("Unsupported data type $data_type")
  end
  CUDA.launch(kernel, cuda_dim..., tuple(output.ptr.p, gradient.ptr.p, blob_dim...))
end

function forward(backend :: GPUBackend, neuron :: Neurons.Sigmoid, output :: Blob)
  cuda_dim, blob_dim = cuda_geometry(neuron, output)
  data_type = eltype(output)
  if data_type == Float32
    kernel = backend.mocha.sigmoid_forward_float
  elseif data_type == Float64
    kernel = backend.mocha.sigmoid_forward_double
  else
    error("Unsupported data type $data_type")
  end
  CUDA.launch(kernel, cuda_dim..., tuple(output.ptr.p, blob_dim...))
end

function backward(backend :: GPUBackend, neuron :: Neurons.Sigmoid, output :: Blob, gradient :: Blob)
  cuda_dim, blob_dim = cuda_geometry(neuron, output)
  data_type = eltype(output)
  if data_type == Float32
    kernel = backend.mocha.sigmoid_backward_float
  elseif data_type == Float64
    kernel = backend.mocha.sigmoid_backward_double
  else
    error("Unsupported data type $data_type")
  end
  CUDA.launch(kernel, cuda_dim..., tuple(output.ptr.p, gradient.ptr.p, blob_dim...))
end

function forward(backend :: GPUBackend, neuron :: Neurons.Tanh, output :: Blob)
  cuda_dim, blob_dim = cuda_geometry(neuron, output)
  data_type = eltype(output)
  if data_type == Float32
    kernel = backend.mocha.tanh_forward_float
  elseif data_type == Float64
    kernel = backend.mocha.tanh_forward_double
  else
    error("Unsupported data type $data_type")
  end
  CUDA.launch(kernel, cuda_dim..., tuple(output.ptr.p, blob_dim...))
end

function backward(backend :: GPUBackend, neuron :: Neurons.Tanh, output :: Blob, gradient :: Blob)
  cuda_dim, blob_dim = cuda_geometry(neuron, output)
  data_type = eltype(output)
  if data_type == Float32
    kernel = backend.mocha.tanh_backward_float
  elseif data_type == Float64
    kernel = backend.mocha.tanh_backward_double
  else
    error("Unsupported data type $data_type")
  end
  CUDA.launch(kernel, cuda_dim..., tuple(output.ptr.p, gradient.ptr.p, blob_dim...))
end

## Exponential

function forward(backend :: GPUBackend, neuron :: Neurons.Exponential, output :: Blob)
  CuVec.exp!(backend, output)
end

function backward(backend :: GPUBackend, neuron :: Neurons.Exponential, output :: Blob, gradient :: Blob)
  CuVec.mul!(backend, gradient, output)
end


text:
"""
    derivative_at(x::Real, values, nodes, baryweights)

Compute the derivative of the function represented by `values` on the `nodes`
at `x` using the corresponding barycentric weights `baryweights`.
[Kopriva, Implementing Spectral Methods for PDEs, Algorithm 36].
"""
function derivative_at(x::Real, values, nodes, baryweights)
    @boundscheck begin
        @assert size(values) == size(nodes) == size(baryweights)
    end
    T = promote_type(typeof(x), eltype(values), eltype(nodes), eltype(baryweights))
    at_node = false
    num = zero(T)
    den = zero(T)
    i = typemin(Int)
    @inbounds for j in eachindex(nodes)
        if x ≈ nodes[j]
            at_node = true
            p = values[j]
            den = -baryweights[j]
            i = j
        end
    end
    if at_node
        @inbounds for j in eachindex(nodes)
            if j != i
                num += baryweights[j] * (p - values[j]) / (x - nodes[j])
            end
        end
    else
        @inbounds p = interpolate(x, values, nodes, baryweights)
        @inbounds for j in eachindex(nodes)
            t = baryweights[j] / (x - nodes[j])
            num += t * (p - values[j]) / (x - nodes[j])
            den += t
        end
    end
    num / den
end

function derivative_at(x::AbstractArray, values, nodes, baryweights)
    @boundscheck begin
        @assert size(values) == size(nodes) == size(baryweights)
    end
    T = promote_type(eltype(x), eltype(values), eltype(nodes), eltype(baryweights))
    res = Array{T}(undef, size(x))
    @inbounds derivative_at!(res, x, values, nodes, baryweights)
    res
end

function derivative_at!(res::AbstractArray, x::AbstractArray, values, nodes, baryweights)
    @boundscheck begin
        @assert size(res) == size(x)
        @assert size(values) == size(nodes) == size(baryweights)
    end
    @inbounds for idx in eachindex(res)
        res[idx] = derivative_at(x[idx], values, nodes, baryweights)
    end
    nothing
end

function derivative_at(x, values, basis::NodalBasis)
    @unpack nodes, baryweights = basis
    derivative_at(x, values, nodes, baryweights)
end

function derivative_at!(res, x, values, basis::NodalBasis)
    @unpack nodes, baryweights = basis
    derivative_at!(res, x, values, nodes, baryweights)
end

"""
    derivative_matrix(nodes, baryweights)

Compute the derivative matrix corresponding to `nodes` and `baryweights`.
[Kopriva, Implementing Spectral Methods for PDEs, Algorithm 37].
"""
function derivative_matrix(nodes, baryweights)
    @boundscheck begin
        @assert length(nodes) == length(baryweights)
    end
    T = promote_type(eltype(nodes), eltype(baryweights))
    mat = Array{T}(undef, length(nodes), length(nodes))
    @inbounds derivative_matrix!(mat, nodes, baryweights)
    mat
end

derivative_matrix(nodes) = derivative_matrix(nodes, barycentric_weights(nodes))

function derivative_matrix!(D::AbstractMatrix, nodes, baryweights)
    @boundscheck begin
        @assert size(D,1) == size(D,2) == length(nodes) == length(baryweights)
    end
    @inbounds for i in eachindex(nodes)
        D[i,i] = 0
        for j in eachindex(nodes)
            if i != j
                D[i,j] = baryweights[j] / baryweights[i] / (nodes[i] - nodes[j])
                D[i,i] -= D[i,j]
            end
        end
    end
    nothing
end


text:
# can define alternate method to change which types are blocked for setproperty! etc
_restricted_fields(::AbstractSample) = (:name, :metadata)

# Samples should at minimum have name and metadata

"""
    name(t::Union{AbstractSample, AbstractFeature})

Get the `name` field from an `AbstractSample` or `AbstractFeature`.
"""
name(as::AbstractSample) = as.name

"""
    metadata(t::AbstractSample)

Get the `metadata` field from an `AbstractSample`.
Note that this is not a copy, so modifications to the returned value
will update the parent `AbstractSample` as well.
"""
metadata(as::AbstractSample) = as.metadata

name(as::AbstractFeature) = as.name

# Allows Strings to be used to index (since `Thing("name", missing)` will be == `Thing("name", something)`)
Base.String(as::AbstractSample) = name(as)
Base.String(af::AbstractFeature) = name(af)

Base.:(==)(s1::T, s2::T) where {T <: Union{AbstractSample, AbstractFeature}} = name(s1) == name(s2)

## getters and setters for metadata ##
function Base.getindex(as::AbstractSample, prop::Symbol)
    prop in _restricted_fields(as) && error("Do not use getindex to access $prop of $(typeof(as)). Use accessor function or getfield instead.")
    getindex(as.metadata, prop)
end

function Base.setindex!(as::AbstractSample, val, prop::Symbol)
    prop in _restricted_fields(as) && error("Do not use setindex! to change $prop of $(typeof(as)).")
    setindex!(as.metadata, val, prop)
    return as
end

function Base.getproperty(as::AbstractSample, prop::Symbol)
    prop in _restricted_fields(as) ? getfield(as, prop) : as.metadata[prop]
end

function Base.setproperty!(as::AbstractSample, prop::Symbol, val)
    prop in _restricted_fields(as) && error("Cannot change $prop of $(typeof(as)) using setproperty!")
    setindex!(as.metadata, val, prop)
    return as
end

function set!(as::AbstractSample, prop::Symbol, val)
    prop in _restricted_fields(as) && error("Cannot set! $prop for $(typeof(as)).")
    set!(as.metadata, prop, val)
    return as
end

function unset!(as::AbstractSample, prop::Symbol)
    prop in _restricted_fields(as) && error("Cannot unset! $prop for $(typeof(as)).")
    unset!(as.metadata, prop)
    return as
end

function insert!(as::AbstractSample, prop::Symbol, val)
    prop in _restricted_fields(as) && error("Cannot insert! $prop for $(typeof(as)).")
    insert!(as.metadata, prop, val)
    return as
end

function delete!(as::AbstractSample, prop::Symbol)
    prop in _restricted_fields(as) && error("Cannot delete! $prop for $(typeof(as)).")
    delete!(as.metadata, prop)
    return as
end

Base.keys(as::AbstractSample) = keys(metadata(as))
Base.haskey(as::AbstractSample, key::Symbol) = in(key, keys(as))

"""
    MicrobiomeSample(name::String, metadata::Dictionary{Symbol, T}) <: AbstractSample
    MicrobiomeSample(name::String)

Microbiome sample type that includes a name and a [`Dictionary`](https://github.com/andyferris/Dictionaries.jl)
of arbitrary metadata using `Symbol`s (other than `:name` or `:metadata`) as keys.

Metadata can be accessed using `getproperty` or `getindex` on the sample itself.

```jldoctest MicrobiomeSample
julia> ms = MicrobiomeSample("sample1", Dictionary([:gender, :age], ["female", 180]))
MicrobiomeSample("sample1", {:gender │ "female", :age │ 180})

julia> name(ms)
"sample1"

julia> ms.name
"sample1"

julia> ms.gender
"female"

julia> ms.age
180
```

Samples can be instantiated with only a name, leaving the `metadata` `Dictionary` blank

```jldoctest MicrobiomeSample
julia> ms2 = MicrobiomeSample("sample2")
MicrobiomeSample("sample2", {})
```

Adding or changing metadata follows the same rules as for the normal `Dictionary` type.

- to change a value use `setproperty` or `setindex`.
  Note that this will fail if the key does not already exist.
- to add or remove a value with validation that it does or does not exist already,
  use `insert!` and `delete!` respectively.
- to add or remove a value without validation ("upsert"),
  use `set!` and `unset!` respectively.

```jldoctest MicrobiomeSample
julia> ms
MicrobiomeSample("sample1", {:gender │ "female", :age │ 180})

julia> ms.age = 16 * 365 # or `ms[:age] = 16 * 365`
5840

julia> set!(ms, :gender, "nonbinary")
MicrobiomeSample("sample1", {:gender │ "nonbinary", :age │ 5840})

julia> insert!(ms, :occupation, "clerk")
MicrobiomeSample("sample1", {:gender │ "nonbinary", :age │ 5840, :occupation │ "clerk"})

julia> insert!(ms, :occupation, "bagger")
ERROR: IndexError("Dictionary already contains index: occupation")

julia> delete!(ms, :occupation)
MicrobiomeSample("sample1", {:gender │ "nonbinary", :age │ 5840})

julia> delete!(ms, :occupation)
ERROR: IndexError("Index doesn't exist: occupation")

julia> unset!(ms, :occupation)
MicrobiomeSample("sample1", {:gender │ "nonbinary", :age │ 5840})
```
"""
struct MicrobiomeSample <: AbstractSample
    name::String
    metadata::Dictionary{Symbol, T} where {T <: Any} # currently non-functional
end

MicrobiomeSample(n::AbstractString) = MicrobiomeSample(n, Dictionary{Symbol, Any}())
# Base.convert(::Type{MicrobiomeSample}, s::AbstractString) = MicrobiomeSample(s)

const _clades = (
    domain     = 0,
    kingdom    = 1,
    phylum     = 2,
    class      = 3,
    order      = 4,
    family     = 5,
    genus      = 6,
    species    = 7,
    subspecies = 8,
    strain     = 9
)

"""
    Taxon(name::String, clade::Union{Missing, Symbol, Int}) <: AbstractFeature
    Taxon(name::String)

Microbial taxon with a name and a clade that can be one of

0. `:domain`
1. `:kingom`
2. `:phylum`
3. `:class`
4. `:order`
5. `:faamily`
6. `:genus`
7. `:species`
8. `:subspecies`
9. `:strain`

or `missing`. Contructors can also use numbers 0-9, or pass a string alone
(in which case the `taxon` will be stored as `missing`)
"""
struct Taxon <: AbstractFeature
    name::String
    clade::Union{Missing, Symbol}

    Taxon(s::AbstractString, ::Missing) = new(s, missing)
    Taxon(s::AbstractString, clade::Symbol) = in(clade, keys(_clades)) ?
                                                        new(s, clade)  :
                                                        error("Invalid clade $clade, must be one of $(keys(_clades))")
end

Taxon(n::AbstractString, clade::Int) = 0 <= clade <= 9 ?
                                            Taxon(n, keys(_clades)[clade+1]) :
                                            error("Invalid clade $clade, must be one of $_clades")
Taxon(n::AbstractString) = Taxon(n, missing)

"""
    clade(t::Union{Taxon, missing})

Get the `clade` field from an `Taxon`.
Returns `missing` if the clade is not set.
"""
clade(t::Taxon) = t.clade
clade(::Missing) = missing

"""
    hasclade(t::Taxon)::Bool

Pretty self-explanatory.
"""
hasclade(t::Taxon) = !ismissing(clade(t))

"""
    GeneFunction(name::String, taxon::Union{Taxon, String, Missing}) <: AbstractFeature
    GeneFunction(name::String)

Microbial gene function object with optional stratification (`taxon`).
"""
struct GeneFunction <: AbstractFeature
    name::String
    taxon::Union{Missing, Taxon}
end

GeneFunction(n::AbstractString) = GeneFunction(n, missing)
GeneFunction(n::AbstractString, t::AbstractString) = GeneFunction(n, Taxon(t))

"""
    taxon(t::Union{GeneFunction, missing})

Get the `taxon` field from a `GeneFunction`.
Returns `missing` if the taxon is not set.
"""
taxon(gf::GeneFunction) = gf.taxon

"""
    hastaxon(t::GeneFunction)::Bool

Pretty self-explanatory.
"""
hastaxon(gf::GeneFunction) = !ismissing(taxon(gf))


text:
module F90

using Libdl
import ..ParticleInCell: TwoDGrid, CloudInCell

const piclib = joinpath(@__DIR__, "libpic")

open(joinpath(@__DIR__, "pic.f90")) do f90file
    open(
        `gfortran -fPIC -w -O3 -shared -x f95 -o $(piclib * "." * Libdl.dlext) -`,
        "w",
    ) do f
        print(f, read(f90file, String))
    end
end

function compute_current!(jx, jy, m::TwoDGrid, kernel :: CloudInCell, p )

    nx = Int32(m.nx)
    ny = Int32(m.ny)
    dx = m.dx
    dy = m.dy
    nbpart = Int32(size(p.array)[2])

    ccall(
        (:deposition, piclib),
        Cvoid,
        (
            Ref{Int32},
            Ref{Int32},
            Ref{Int32},
            Ref{Float64},
            Ref{Float64},
            Ptr{Float64},
            Ptr{Float64},
            Ptr{Float64},
        ),
        nbpart,
        nx,
        ny,
        dx,
        dy,
        p.array,
        jx,
        jy,
    )

    for i = 1:nx+1
        jx[i, 1] += jx[i, ny+1]
        jx[i, ny+1] = jx[i, 1]
    end
    for j = 1:ny+1
        jx[1, j] += jx[nx+1, j]
        jx[nx+1, j] = jx[1, j]
    end

end

function push_x!(p, mesh::TwoDGrid, dt::Float64)

    nbpart = Int32(size(p.array)[2])
    dimx = mesh.dimx
    dimy = mesh.dimy

    ccall(
        (:push_x, piclib),
        Cvoid,
        (Ref{Int32}, Ref{Float64}, Ref{Float64}, Ptr{Float64}, Ref{Float64}),
        nbpart,
        dimx,
        dimy,
        p.array,
        dt,
    )

end

export push_v!

function push_v!(p, kernel::CloudInCell, m::TwoDGrid, ex, ey, bz, dt::Float64)

    nbpart = Int32(size(p.array)[2])
    nx = Int32(m.nx)
    ny = Int32(m.ny)
    dx = m.dx
    dy = m.dy

    ccall((:push_v, piclib), Cvoid, (Ref{Int32}, Ptr{Float64}, Ref{Float64},
            Ref{Int32},
            Ref{Int32},
            Ref{Float64},
            Ref{Float64},
            Ptr{Float64},
            Ptr{Float64},
            Ptr{Float64}), nbpart, p.array, dt, nx, ny, dx, dy, ex, ey, bz)

end

end


text:
##################################################################################
# This file is part of StateSpaceEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020-2022, Bank of Canada
# All rights reserved.
##################################################################################

# export StackedTimeSolverData
"""
    StackedTimeSolverData

The data structure used in the stacked time algorithm.

**TODO** Add all the details here.
"""
struct StackedTimeSolverData # <: AbstractSolverData
    "Number of time periods"
    NT::Int
    "Number of variables"
    NV::Int
    "Number of shocks"
    NS::Int
    "Number of unknowns = variables + shocks + aux variables"
    NU::Int
    "Final condition type: `fcgiven` or `fclevel` or `fcslope`."
    FC::Vector{<:FinalCondition}
    "List of vectors for the time indexes of blocks of the matrix."
    TT::Vector{Vector{Int}}
    "List of vectors for the row indexes of each block."
    II::Vector{Vector{Int}}
    "List of vectors for the column indexes of each block."
    JJ::Vector{Vector{Int}}
    "Inverse lookup for the column indexes."
    BI::Vector{Vector{Int}}
    "Jacobian matrix cache."
    J::SparseMatrixCSC{Float64,Int}
    "Correction to J in case FC == fcslope. See description of algorithm."
    CiSc::SparseMatrixCSC{Float64,Int}   # Correction to J to account for FC in the case of fcslope
    "The steady state data, used for FC ∈ (fclevel, fcslope)."
    SS::AbstractArray{Float64,2}
    "Mask of variables with constant growth rate in the steady state"
    log_mask::AbstractVector{Bool}
    "Mask of variables with constant change in the steady state"
    lin_mask::AbstractVector{Bool}
    "The `evaldata` from the model, used to evaluate RJ and R!"
    evaldata::ModelBaseEcon.AbstractModelEvaluationData
    "Keep track of which variables are set by exogenous constraints."
    exog_mask::AbstractVector{Bool}
    "Keep track of which variables are set by final conditions."
    fc_mask::AbstractVector{Bool}
    "Keep track of which variables are \"active\" for the solver."
    solve_mask::AbstractVector{Bool}
    "Cache the LU decomposition of the active part of J"
    luJ::Ref{Any}
end

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

# var_CiSc(::StackedTimeSolverData, ::ModelVariable, ::FinalCondition) = error("Missing var_CiSc for $(typeof(fc))")
# assign_fc!(::AbstractVector{Float64}, ::AbstractVector{Float64}, ::Int, fc::FinalCondition) = error("Missing assign_fc! for $(typeof(fc))")

#######

@inline function var_CiSc(::StackedTimeSolverData, ::ModelVariable, ::FCNone)
    # No Jacobian correction
    return Dict{Int,Vector{Float64}}()
end

@inline function assign_fc!(x::AbstractVector{Float64}, exog::AbstractVector{Float64}, ::Int, sd::StackedTimeSolverData, ::FCNone)
    # nothing to see here
    return x
end

#######

@inline function var_CiSc(::StackedTimeSolverData, ::ModelVariable, ::FCGiven)
    # No Jacobian correction
    return Dict{Int,Vector{Float64}}()
end

@inline function assign_fc!(x::AbstractVector{Float64}, exog::AbstractVector{Float64}, ::Int, sd::StackedTimeSolverData, ::FCGiven)
    # values given exogenously
    if x !== exog
        x[sd.TT[end]] .= exog[sd.TT[end]]
    end
    return x
end

#######

@inline function var_CiSc(::StackedTimeSolverData, ::ModelVariable, ::FCMatchSSLevel)
    # No Jacobian correction
    return Dict{Int,Vector{Float64}}()
end

@inline function assign_fc!(x::AbstractVector{Float64}, ::AbstractVector{Float64}, vind::Int, sd::StackedTimeSolverData, ::FCMatchSSLevel)
    # values come from the steady state
    x[sd.TT[end]] .= sd.SS[sd.TT[end], vind]
    return x
end

#######

@inline function var_CiSc(sd::StackedTimeSolverData, ::ModelVariable, ::FCMatchSSRate)
    # Matrix C is composed of blocks looking like this:
    #    1  0  0  0
    #   -1  1  0  0
    #    0 -1  1  0
    #    0  0 -1  1
    # The corresponding rows in Sc look like this
    #   0 0 0 0 0 -1
    #   0 0 0 0 0  0
    #   0 0 0 0 0  0
    #   0 0 0 0 0  0
    # The C^{-1} * Sc matrix then looks like this
    #   0 0 0 0 0 -1
    #   0 0 0 0 0 -1
    #   0 0 0 0 0 -1
    #   0 0 0 0 0 -1
    # The number of rows equals NTFC (the number of final condition periods)
    # The column with the -1 is the one corresponding to the last simulation period of the variable.
    return Dict{Int,Vector{Float64}}(0 => fill(-1.0, length(sd.TT[end])))
end

@inline function assign_fc!(x::AbstractVector{Float64}, ::AbstractVector{Float64}, vind::Int, sd::StackedTimeSolverData, ::FCMatchSSRate)
    # compute values using the slope from the steady state
    for t in sd.TT[end]
        x[t] = x[t-1] + sd.SS[t, vind] - sd.SS[t-1, vind]
    end
    return x
end

#######

@inline function var_CiSc(sd::StackedTimeSolverData, ::ModelVariable, ::FCConstRate)
    # Matrix C is composed of blocks looking like this:
    #    -1  0  0  0  0
    #     2 -1  1  0  0
    #    -1  2 -1  0  0
    #     0 -1  2 -1  0
    #     0  0 -1  2 -1
    # The corresponding rows in Sc look like this
    #     0 ..  0 -1  2
    #     0 ..  0  0 -1
    #     0 ..  0  0  0
    #     0 ..  0  0  0
    #     0 ..  0  0  0
    # The C^{-1} * Sc matrix then looks like this
    #     0 ..  0  1 -2
    #     0 ..  0  2 -3
    #     0 ..  0  3 -4
    #     0 ..  0  4 -5
    #     0 ..  0  5 -6

    # The number of rows equals sd.TT[end] (the number of final condition periods)
    # The column with -2,-3,-4... has offset 0 -- the last simulation period.
    # The column with  1, 2, 3... has offset -1 -- the period before the last simulation period.
    foo = collect(Float64, 1:length(sd.TT[end]))
    return Dict{Int,Vector{Float64}}(-1 => foo, 0 => -1.0 .- foo)
end

@inline function assign_fc!(x::AbstractVector{Float64}, ::AbstractVector{Float64}, ::Int, sd::StackedTimeSolverData, ::FCConstRate)
    # compute values using the slope at the end of the simulation
    for t in sd.TT[end]
        x[t] = 2x[t-1] - x[t-2]
    end
    return x
end

"""
    update_plan!(sd::StackedTimeSolverData, model, plan; changed=false)

Implementation for the Stacked Time algorithm. Plan must have the same range as
the original plan.

By default the data structure is updated only if an actual change in the plan is
detected. Setting the `changed` flag to `true` forces the update even if the
plan seems unchanged. This is necessary only in rare circumstances.

"""
function update_plan!(sd::StackedTimeSolverData, m::Model, p::Plan; changed = false)
    if sd.NT != length(p.range)
        error("Unable to update using a simulation plan of different length.")
    end

    # LinearIndices used for indexing the columns of the global matrix
    LI = LinearIndices((p.range, 1:sd.NU))

    sim = m.maxlag+1:sd.NT-m.maxlead
    NTFC = m.maxlead

    # Assume initial conditions are set correctly to exogenous in the constructor
    # We only update the masks during the simulation periods
    foo = BitArray(undef, sd.NU)
    for t in sim
        # s = p[t]
        # si = indexin(s, unknowns)
        si = p[t, Val(:inds)]
        fill!(foo, false)
        foo[si] .= true
        if !all(foo .== sd.exog_mask[LI[t, :]])
            sd.exog_mask[LI[t, :]] .= foo
            changed = true
        end
    end

    # Update the solve_mask array
    if changed
        @. sd.solve_mask = !(sd.fc_mask | sd.exog_mask)
        @assert !any(sd.exog_mask .& sd.fc_mask)
        @assert sum(sd.solve_mask) == size(sd.J, 1)

        # Update the Jacobian correction matrix, if exogenous plan changed
        II, JJ, VV = Int[], Int[], Float64[]
        last_sim_t = last(sim)
        for (vi, (v, fc)) in enumerate(zip(m.allvars, sd.FC))
            @assert vi == v.index
            # var_CiSc returns a Dict with keys equal to the column offset relative to last_sim_t
            # and values containing the column values.
            for (offset, values) in var_CiSc(sd, v, fc)
                # col_ind is the column index in the global matrix
                col_ind = LI[last_sim_t+offset, vi]
                if !sd.solve_mask[col_ind]
                    # this value is exogenous, no column for it in Sc
                    continue
                end
                # The column index in the active matrix
                JJvar = let
                    # Full matrix has sd.NU*sd.NT columns
                    # We take a vector with this many `false`s and put a `true` only in `col_ind` position
                    # Then we select only the active positions in this vector (that's foo[sd.solve_mask])
                    # Then we find the index of the `true` value - that's the column index in the active sub-matrix.
                    foo = falses(sd.NU * sd.NT)
                    foo[col_ind] = true
                    ind = findall(foo[sd.solve_mask])
                    repeat(ind, NTFC)
                end
                # The row indices for this block
                #   the block for variable 1 gets row indexes 1 .. NTFC
                #   the block for variable 2 gets row indexes NTFC+1 .. 2*NTFC
                #   the block for variable vi gets row indexes (vi-1)*NTFC+1 .. vi*NTFC
                IIvar = ((vi-1)*NTFC+1):(vi*NTFC)
                # The values in this column came from the var_CiSc call above
                VVvar = values
                append!(II, IIvar)
                append!(JJ, JJvar)
                append!(VV, VVvar)
            end
            # Construct the sparse matrix
            sd.CiSc .= sparse(II, JJ, VV, NTFC * sd.NU, size(sd.J, 1))
        end
        # cached lu is no longer valid, since active columns have changed
        sd.luJ[] = nothing
    end

    return sd
end

import ..steadystatearray

"""
    make_BI(J, II)

Prepares the `BI` array for the solver data. Called from the constructor of
`StackedTimeSolverData`.

"""
function make_BI(JMAT::SparseMatrixCSC{Float64,Int}, II::AbstractVector{<:AbstractVector{Int}})
    # Computes the set of indexes in JMAT.nzval corresponding to blocks of equations in II

    # Bar is the inverse map of II, i.e. Bar[j] = i <=> II[i] contains j
    local Bar = zeros(Int, size(JMAT, 1))
    for i in axes(II, 1)
        Bar[II[i]] .= i
    end

    # start with empty lists in BI
    local BI = [[] for _ in II]

    # iterate the non-zero elements of JMAT and record their row numbers in BI
    local rows = rowvals(JMAT)
    local n = size(JMAT, 2)
    for i = 1:n
        for j in nzrange(JMAT, i)
            row = rows[j]
            push!(BI[Bar[row]], j)
        end
    end
    return BI
end

"""
    StackedTimeSolverData(model, plan, fctype)

"""
@inline StackedTimeSolverData(m::Model, p::Plan, fctype::FinalCondition) = StackedTimeSolverData(m, p, setfc(m, fctype))
function StackedTimeSolverData(m::Model, p::Plan, fctype::AbstractVector{FinalCondition})

    NT = length(p.range)
    init = 1:m.maxlag
    term = NT-m.maxlead+1:NT
    sim = m.maxlag+1:NT-m.maxlead
    NTFC = length(term)
    NTSIM = length(sim)
    NTIC = length(init)

    need_SS = false
    for fc in fctype
        if (fc === fclevel || fc === fcslope)
            need_SS = true
        end
        if (fc === fcnatural) && (length(sim) + length(init) < 2)
            throw(ArgumentError("Simulation must include at least 2 periods for `fcnatural`."))
        end
    end

    if need_SS && !issssolved(m)
        # NOTE: we do not verify the steady state, just make sure it's been assigned
        throw(ArgumentError("Steady state must be solved for `fclevel` or `fcslope`."))
    end

    unknowns = m.allvars
    nunknowns = length(unknowns)

    nvars = length(m.variables)
    nshks = length(m.shocks)
    nauxs = length(m.auxvars)

    equations = m.alleqns
    nequations = length(equations)

    # LinearIndices used for indexing the columns of the global matrix
    LI = LinearIndices((p.range, 1:nunknowns))

    # Initialize empty arrays
    TT = Vector{Int}[] # the time indexes of the block, used when updating values in global_RJ
    II = Vector{Int}[] # the row indexes
    JJ = Vector{Int}[] # the column indexes

    # Prep the Jacobian matrix
    neq = 0 # running counter of equations added to matrix
    # Model equations are the same for each sim period, just shifted according to t
    Jblock = [ti + NT * (vi - 1) for eqn in equations for (ti, vi) in eqn.vinds]
    Iblock = [i for (i, eqn) in enumerate(equations) for _ in eqn.vinds]
    Tblock = -m.maxlag:m.maxlead
    @timer for t in sim
        push!(TT, t .+ Tblock)
        push!(II, neq .+ Iblock)
        push!(JJ, t .+ Jblock)
        neq += nequations
    end
    # Construct the
    @timer begin
        I = vcat(II...)
        J = vcat(JJ...)
        JMAT = sparse(I, J, fill(NaN64, size(I)), NTSIM * nequations, NT * nunknowns)
    end

    # We no longer need the exact indexes of all non-zero entires in the Jacobian matrix.
    # We do however need the set of equation indexes for each sim period
    @timer foreach(unique! ∘ sort!, II)

    # BI holds the indexes in JMAT.nzval for each block of equations
    BI = make_BI(JMAT, II)  # same as the two lines above, but faster

    # We also need the times of the final conditions
    push!(TT, term)

    #
    exog_mask = falses(nunknowns * NT)
    # Initial conditions are set as exogenous values
    exog_mask[vec(LI[init, :])] .= true
    # The exogenous values during sim are set in update_plan!() call below.

    # Final conditions are complicated
    fc_mask = falses(nunknowns * NT)
    fc_mask[vec(LI[term, :])] .= true

    # The solve_mask is redundant. We pre-compute and store it for speed
    solve_mask = .!(fc_mask .| exog_mask)

    sd = StackedTimeSolverData(NT, nvars, nshks, nunknowns, fctype, TT, II, JJ, BI, JMAT,
        sparse([], [], Float64[], NTFC * nunknowns, size(JMAT, 1)),
        need_SS ? transform(steadystatearray(m, p), m) : zeros(0, 0),
        islog.(m.allvars) .| isneglog.(m.allvars), islin.(m.allvars),
        m.evaldata, exog_mask, fc_mask, solve_mask,
        Ref{Any}(nothing))

    return @timer update_plan!(sd, m, p; changed = true)
end

function assign_exog_data!(x::AbstractArray{Float64,2}, exog::AbstractArray{Float64,2}, sd::StackedTimeSolverData)
    # @assert size(x,1) == size(exog,1) == sd.NT
    x[sd.exog_mask] = exog[sd.exog_mask]
    assign_final_condition!(x, exog, sd)
    return x
end

function assign_final_condition!(x::AbstractArray{Float64,2}, exog::AbstractArray{Float64,2}, sd::StackedTimeSolverData)
    @timer for (vi, fc) in enumerate(sd.FC)
        assign_fc!(view(x, :, vi), exog[:, vi], vi, sd, fc)
    end
    return x
end

function global_R!(res::AbstractArray{Float64,1}, point::AbstractArray{Float64}, exog_data::AbstractArray{Float64}, sd::StackedTimeSolverData)
    @assert size(point) == size(exog_data) == (sd.NT, sd.NU)
    # point = reshape(point, sd.NT, sd.NU)
    # exog_data = reshape(exog_data, sd.NT, sd.NU)
    @assert(length(res) == size(sd.J, 1), "Length of residual vector doesn't match.")
    @timer "global_R!" begin
        for (ii, tt) in zip(sd.II, sd.TT)
            eval_R!(view(res, ii), point[tt, :], sd.evaldata)
        end
    end
    return res
end

#= disable - this is not necessary with log-transformed variables

@inline update_CiSc!(x::AbstractArray{Float64,2}, sd::StackedTimeSolverData) = any(sd.log_mask) ? update_CiSc!(x, sd, Val(sd.FC)) : nothing

@inline update_CiSc!(x, sd, ::Val{fcgiven}) = nothing
@inline update_CiSc!(x, sd, ::Val{fclevel}) = nothing
function update_CiSc!(x, sd::StackedTimeSolverData, ::Val{fcslope})
    # LinearIndices used for indexing the columns of the global matrix
    local LI = LinearIndices((1:sd.NT, 1:sd.NU))
    # The last simulation period
    local term_Ts = sd.TT[end]
    if isempty(term_Ts)
        return
    end
    local last_T = term_Ts[1] - 1
    local NTFC = length(term_Ts)

    # @info "last_T = $(last_T), NTFC = $(NTFC)"

    # Iterate over the non-zero values of CiSc and update
    # values as necessary
    local rows = rowvals(sd.CiSc)
    local vals = nonzeros(sd.CiSc)
    local ncols = size(sd.CiSc, 2)
    for col = 1:ncols
        nzr = nzrange(sd.CiSc, col)
        if isempty(nzr)
            continue
        end
        # Compute variable and time index from column index
        var, tm = divrem(LI[sd.solve_mask][col] - 1, sd.NT)  .+ 1

        # @info "var = $(var), tm = $(tm)"

        if !sd.log_mask[var]
            continue
        end
        if tm != last_T
            error("Last simulation times don't match")
        end
        for j = nzr
            row = rows[j]
            # update vals[j]
            local T = rem(row - 1, NTFC) + 1 + last_T

            # @info "Updating $((row, col)) with T=$(T), x[T, var] = $(x[T, var]), x[last_T, var] = $(x[last_T, var])"

            vals[j] = - x[T, var] / x[last_T, var]
        end
    end
    return nothing
end

function update_CiSc!(x, sd, ::Val{fcnatural})
    # LinearIndices used for indexing the columns of the global matrix
    local LI = LinearIndices((1:sd.NT, 1:sd.NU))
    # The last simulation period
    local term_Ts = sd.TT[end]
    if isempty(term_Ts)
        return
    end
    local last_T = term_Ts[1] - 1
    local NTFC = length(term_Ts)

    # @info "last_T = $(last_T), NTFC = $(NTFC)"

    # Iterate over the non-zero values of CiSc and update
    # values as necessary
    local rows = rowvals(sd.CiSc)
    local vals = nonzeros(sd.CiSc)
    local ncols = size(sd.CiSc, 2)
    for col = 1:ncols
        nzr = nzrange(sd.CiSc, col)
        if isempty(nzr)
            continue
        end
        # Compute variable and time index from column index
        var, tm = divrem(LI[sd.solve_mask][col] - 1, sd.NT)  .+ 1

        # @info "var = $(var), tm = $(tm)"

        if !sd.log_mask[var]
            continue
        end
        if tm == last_T
            s = -1
            c = 2
        elseif tm == last_T - 1
            s = +1
            c = 1
        else
            error("Last simulation times don't match")
        end
        for j = nzr
            row = rows[j]
            # update vals[j]
            local T = rem(row - 1, NTFC) + 1 + last_T

            # @info "Updating $((row, col)) with T=$(T), x[T, var] = $(x[T, var]), x[last_T, var] = $(x[last_T, var])"

            vals[j] = s * (c + T - last_T - 1 ) * x[T, var] / x[tm, var]
        end
    end
    return nothing
end =#

function global_RJ(point::AbstractArray{Float64}, exog_data::AbstractArray{Float64}, sd::StackedTimeSolverData;
    debugging = false)
    nunknowns = sd.NU
    @assert size(point) == size(exog_data) == (sd.NT, nunknowns)
    # point = reshape(point, sd.NT, nunknowns)
    # exog_data = reshape(exog_data, sd.NT, nunknowns)
    JAC = sd.J
    RES = Vector{Float64}(undef, size(JAC, 1))
    # Model equations @ [1] to [end-1]
    haveJ = isa(sd.evaldata, ModelBaseEcon.LinearizedModelEvaluationData) && !any(isnan.(JAC.nzval))
    haveLU = sd.luJ[] !== nothing
    if haveJ
        # update only RES
        @timer "globalRJ/evalR!" for i = 1:length(sd.BI)
            eval_R!(view(RES, sd.II[i]), point[sd.TT[i], :], sd.evaldata)
        end
    else
        # update both RES and JAC
        @timer "globalRJ/evalRJ" for i = 1:length(sd.BI)
            R, J = eval_RJ(point[sd.TT[i], :], sd.evaldata)
            RES[sd.II[i]] .= R
            JAC.nzval[sd.BI[i]] .= J.nzval
        end
        # update_CiSc!(point, sd)
        haveLU = false
    end
    if !haveLU
        if nnz(sd.CiSc) > 0
            JJ = JAC[:, sd.solve_mask] - JAC[:, sd.fc_mask] * sd.CiSc
        else
            JJ = JAC[:, sd.solve_mask]
        end
        try
            # compute lu decomposition of the active part of J and cache it.
            @timer "LU decomposition" sd.luJ[] = lu(JJ)
        catch e
            if e isa SingularException
                @error("The system is underdetermined with the given set of equations and final conditions.")
            end
            debugging && return RES, deepcopy(JJ)
            rethrow()
        end
    end
    return RES, sd.luJ[]
end

function assign_update_step!(x::AbstractArray{Float64}, λ::Float64, Δx::AbstractArray{Float64}, sd::StackedTimeSolverData)
    x[sd.solve_mask] .+= λ .* Δx
    if nnz(sd.CiSc) > 0
        x[sd.fc_mask] .-= λ .* (sd.CiSc * Δx)
        # if sd.FC == fcnatural && any(sd.log_mask) && !isempty(sd.TT[end])
        #     for t = sd.TT[end]
        #         @. x[t, sd.log_mask] = x[t - 1, sd.log_mask]^2 / x[t - 2,sd.log_mask]
        #     end
        # end
        # assign_final_condition!(x, zeros(0,0), sd, Val(fcslope))
    end
    return x
end


text:
# Dantzig-Wolfe Reformulation and Column Generation
# Data for General Assignment Problem 1
# Edward J. Xu
# 2019.5.1

function getData()
    vec_c = [#=
        =# 9 9 1 4 6 4 2 3 5 1 9 9 7 6 3  #=
        =# 9 1 8 6 7 8 2 6 6 5 3 4 7 5 3  #=
        =# 4 8 1 8 1 4 1 4 2 6 2 1 2 5 1  #=
        =# 7 9 7 8 3 6 2 3 4 7 9 9 3 9 1
        ]
    mat_a = [
        Matrix{Float64}(I, 15, 15) Matrix{Float64}(I, 15, 15) Matrix{Float64}(I, 15, 15) Matrix{Float64}(I, 15, 15);
        8 6 1 7 7 7 5 5 7 7 3 3 8 5 4 zeros(1,45) ;
        zeros(1,15) 5 5 3 8 5 6 5 9 9 6 1 9 5 6 3 zeros(1,30);
        zeros(1,30) 3 2 4 4 9 1 7 3 3 3 5 3 7 7 1 zeros(1,15);
        zeros(1,45) 6 7 9 9 3 5 2 1 5 5 4 4 6 2 1
        ]
    vec_b = hcat([ones(15,1); 22.0; 18.0; 18.0; 19.0])
    indexSub = [collect(1: 15), collect(16: 30), collect(31: 45), collect(46: 60)]
    numSub = 4
    num_x_sub = 15
    num_border = 15
    #
    return (vec_c, mat_a, vec_b, indexSub, numSub, num_x_sub, num_border, vecStrNameVar)
end

(vec_c, mat_a, vec_b, indexSub, numSub, num_x_sub, num_border, vecStrNameVar) = getData()


text:
using DataFrames
using HDF5
using PhysicalConstants
using PhysicalConstants.CODATA2018
using Unitful

import Unitful:
    nm, μm, mm, cm, ns, μs, ms, ps, s

@enum Dtsel dtfirst dtminimum dtaverage
@enum Possel postrue posreco posrecall

"""
	MlemLor

Struct representing a LOR
"""
struct MlemLor{T}
    dx::T
    x1::T
    y1::T
    z1::T
    x2::T
    y2::T
    z2::T
end

"""
	write_lors_hdf5(filename, mlor)
	Write lors in a hdf5 format required by petalorust (mlem algo)
"""
function write_lors_hdf5(filename, mlor)

    function set_datatype(::Type{MlemLor{Float32}})
        dtype = HDF5.h5t_create(HDF5.H5T_COMPOUND, sizeof(MlemLor{Float32}))
        HDF5.h5t_insert(dtype, "dx", fieldoffset(MlemLor{Float32}, 1),
                        datatype(Float32))
        HDF5.h5t_insert(dtype, "x1", fieldoffset(MlemLor{Float32}, 2),
                        datatype(Float32))
        HDF5.h5t_insert(dtype, "y1", fieldoffset(MlemLor{Float32}, 3),
                        datatype(Float32))
        HDF5.h5t_insert(dtype, "z1", fieldoffset(MlemLor{Float32}, 4),
                        datatype(Float32))
        HDF5.h5t_insert(dtype, "x2", fieldoffset(MlemLor{Float32}, 5),
                        datatype(Float32))
        HDF5.h5t_insert(dtype, "y2", fieldoffset(MlemLor{Float32}, 6),
                        datatype(Float32))
        HDF5.h5t_insert(dtype, "z2", fieldoffset(MlemLor{Float32}, 7),
                        datatype(Float32))

        HDF5.Datatype(dtype)
    end

    h5f = h5open(filename, "w")

    dtype = set_datatype(MlemLor{Float32})
    dspace = dataspace(mlor)
    grp = create_group(h5f, "true_info")
    dset = create_dataset(grp, "lors", dtype, dspace)
    write_dataset(dset, dtype, mlor)

    close(h5f)
end

"""
	radial_correction(xb::Vector{T}, yb::Vector{T},zb::Vector{T},
	                            r::Vector{T}) where T
Computes the radial correction from barycenter
"""
function radial_correction(xb::Vector{T}, yb::Vector{T},zb::Vector{T},
                            r::Vector{T}) where T
    ϕ = atan.(yb,xb)
    x = r .* cos.(ϕ)
    y = r .* sin.(ϕ)
    return x,y,zb
end

"""
    deltatime(df::DataFrame, t::Dtsel=dtfirst)
    @enum Dtsel dtfirst dtminimum dtaverage

Return t2 - t1 where t1 and t2 are:
dtfirst for nominal (true) time
dtminimum for minimum time stamp of all SiPMs
dtaverage for average time stamp of 5 minimum SiPMs
"""
function deltatime(df::DataFrame, t::Dtsel=dtfirst)
    if t     == dtminimum
        return  uconvert.(ps, df.tr2 - df.tr1)
    elseif t == dtaverage
        return uconvert.(ps,df.ta2 - df.ta1)
    else
        return uconvert.(ps, df.t2 - df.t1)

    end
end

"""
    cdoi(df::DataFrame, position::Possel=postrue, nlxe::Number=1.6)
    @enum Possel postrue posreco

Return the time correction associated to DOI (depth of interaction)
position:
  posreco: if using reconstructed values
  postrue: if using true values

"""
function cdoi(df::DataFrame, position::Possel=postrue, nlxe::Number=1.6)
    clxe = SpeedOfLightInVacuum/nlxe

    if position == posreco
        dxrb1 = [dxyz([df.x1[i], df.y1[i], df.z1[i]],
                            [df.xb1[i], df.yb1[i], df.zb1[i]]) for i in 1:nrow(df)]

        dxrb2 = [dxyz([df.x2[i], df.y2[i], df.z2[i]],
                            [df.xb2[i], df.yb2[i], df.zb2[i]]) for i in 1:nrow(df)]

    elseif position == posrecall

        xq1,yq1,zq1 =  radial_correction(df.xr1./mm, df.yr1./mm, df.zr1./mm, df.r1x./mm)
        xq2,yq2,zq2 =  radial_correction(df.xr2./mm, df.yr2./mm, df.zr2./mm, df.r2x./mm)

        dxrb1 = [dxyz([xq1[i]*mm, yq1[i]*mm, zq1[i]*mm],
                            [df.xb1[i], df.yb1[i], df.zb1[i]]) for i in 1:nrow(df)]

        dxrb2 = [dxyz([xq2[i]*mm, yq2[i]*mm, zq2[i]*mm],
                            [df.xb2[i], df.yb2[i], df.zb2[i]]) for i in 1:nrow(df)]

    else
        dxrb1 = [dxyz([df.xt1[i], df.yt1[i], df.zt1[i]],
                            [df.xb1[i], df.yb1[i], df.zb1[i]]) for i in 1:nrow(df)]

        dxrb2 = [dxyz([df.xt2[i], df.yt2[i], df.zt2[i]],
                            [df.xb2[i], df.yb2[i], df.zb2[i]]) for i in 1:nrow(df)]

    end
    return uconvert.(ps, (dxrb2 - dxrb1)/clxe)
end

function ctsr(df::DataFrame, position::Possel=postrue)
    cc = float(SpeedOfLightInVacuum)
    if position == postrue
        tsr1 = [dxyz([df.xt1[i], df.yt1[i], df.zt1[i]],
                             [df.xs[i], df.ys[i], df.zs[i]]) for i in 1:nrow(df)]/cc
        tsr2 = [dxyz([df.xt2[i], df.yt2[i], df.zt2[i]],
                             [df.xs[i], df.ys[i], df.zs[i]]) for i in 1:nrow(df)]/cc

    else
        tsr1 = [dxyz([df.x1[i], df.y1[i], df.z1[i]],
                             [df.xs[i], df.ys[i], df.zs[i]]) for i in 1:nrow(df)]/cc
        tsr2 = [dxyz([df.x2[i], df.y2[i], df.z2[i]],
                             [df.xs[i], df.ys[i], df.zs[i]]) for i in 1:nrow(df)]/cc
    end
    return uconvert.(ps, tsr2 - tsr1)
end

"""
    crt(dfu, dtsel=dtfirst, posel=postrue)

Return the CRT of the system
"""
function crt(dfu, dtsel=dtfirst, posel=postrue)

    dt12 = deltatime(dfu, dtsel)
    t12 = dt12./ps

    dtsr12 = ctsr(dfu, postrue)   # this is a nominal position for CRT
    tsr12 = dtsr12./ps

    dtrb12 = cdoi(dfu,  posel)
    trb12 = dtrb12 ./ps

    dt = t12 - tsr12 - trb12

end

function dxf(df::DataFrame, t::Dtsel=dtfirst, position::Possel=postrue)
    dt12  = deltatime(df,t)
    dtdoi = cdoi(df, position)
    # compute dx from time and speed of light, ensure that the result is in mm
    uconvert.(mm, (dt12 - dtdoi) * SpeedOfLightInVacuum)
end

"""
    dftolor(df::DataFrame, t::Dtsel=dtfirst, position::Possel=postrue, nlxe::Number=1.6)

Take dataframe df and return a vector of MlemLor.
Since MlemLor will be written to file, remove units (e.g, use implicit units, in this
case mm) and transform to Float32
"""
function dftolor(df::DataFrame, t::Dtsel=dtfirst, position::Possel=postrue; f32::Bool=true)

    function tof32(l)
        Float32.(l/mm)
    end

    function tof64(l)
        Float64.(l/mm)
    end

    @debug "f32 =", f32

    dx    = dxf(df, t, position)

    if position == postrue
        x1, x2, y1, y2, z1, z2 = df.xt1, df.xt2, df.yt1, df.yt2, df.zt1, df.zt2
        if f32
            return MlemLor{Float32}.(tof32(dx),tof32(x1),tof32(y1),tof32(z1),
                            tof32(x2), tof32(y2), tof32(z2))
        else
            return MlemLor{Float64}.(tof64(dx),tof64(x1),tof64(y1),tof64(z1),
                            tof64(x2), tof64(y2), tof64(z2))
        end
    elseif position == posreco
        x1, x2, y1, y2, z1, z2 = df.x1,  df.x2,  df.y1,  df.y2,  df.z1,  df.z2
        if f32
            return MlemLor{Float32}.(tof32(dx),tof32(x1),tof32(y1),tof32(z1),
                            tof32(x2), tof32(y2), tof32(z2))
        else
            return MlemLor{Float64}.(tof64(dx),tof64(x1),tof64(y1),tof64(z1),
                            tof64(x2), tof64(y2), tof64(z2))
        end
    else
        x1,y1,z1 =  radial_correction(df.xr1./mm, df.yr1./mm, df.zr1./mm, df.r1x./mm)
        x2,y2,z2 =  radial_correction(df.xr2./mm, df.yr2./mm, df.zr2./mm, df.r2x./mm)
        if f32
            return MlemLor{Float32}.(tof32(dx),
                            Float32.(x1),Float32.(y1),Float32.(z1),
                            Float32.(x2), Float32.(y2), Float32.(z2))
        else
            @debug "f32 =", f32
            @debug "dx1 =", tof64(dx)[1]
            @debug "x1  =", Float64.(x1)[1]

            return MlemLor{Float64}.(tof64(dx),Float64.(x1),Float64.(y1),Float64.(z1),
                            Float64.(x2), Float64.(y2), Float64.(z2))
        end
    end
end

max_stars_repo_name: jjgomezcadenas/mytest

max_issues_repo_name: jjgomezcadenas/mytest

max_forks_repo_name: jjgomezcadenas/mytest


text:
include("BUDE.jl")
using StaticArrays
using Enzyme
using FiniteDifferences
import Base.Experimental: @aliasscope

# Enzyme.API.printperf!(true)
# Enzyme.API.printall!(true)

const Device = (undef, "CPU", "Threaded")

function devices()
  return [Device]
end

const Float32Max = typemax(Float32)

function run(params::Params, deck::Deck, _::DeviceWithRepr)
  println("Using max $(Threads.nthreads()) threads")
  if params.ppwi != DefaultPPWI
    @warn "Threaded implementation only uses wgsize, the PPWI argument is ignored"
  end

  poses = size(deck.poses)[2]
  etotals = Array{Float32}(undef, poses)

  d_protein = fill!(similar(deck.protein), Atom(0f0, 0f0, 0f0, 0))
  d_etotals = fill!(similar(etotals), 0)
  # warmup
  if params.enzyme
     d_etotals[1] = 1
     @time autodiff(fasten_main, Const,
        Val(convert(Int, params.wgsize)),
        Duplicated(deck.protein, d_protein),
        Const(deck.ligand),
        Const(deck.forcefield),
        Const(deck.poses),
        Duplicated(etotals, d_etotals)
      )
      if params.verify

        etotals0 = Array{Float32}(undef, poses)
	fasten_main(
	      Val(convert(Int, params.wgsize)),
	      deck.protein,
	      deck.ligand,
	      deck.forcefield,
	      deck.poses,
	      etotals0,
       )
        function fwdwrap(x)
           tmp = copy(deck.protein)
	   tmp[13] = Atom(x, tmp[13].y, tmp[13].z, tmp[13].type)
	   fetotals = Array{Float32}(undef, poses)
           fasten_main(Val(convert(Int, params.wgsize)), tmp, deck.ligand, deck.forcefield, deck.poses, fetotals)
           fetotals[1]
        end
	fd = FiniteDifferences.central_fdm(5, 1)(fwdwrap, deck.protein[13].x)
	@show d_protein[13], fd
      end
  else
    @time fasten_main(
      Val(convert(Int, params.wgsize)),
      deck.protein,
      deck.ligand,
      deck.forcefield,
      deck.poses,
      etotals,
    )
  end

  if params.enzyme
    elapsed = @elapsed for _ = 1:params.iterations
      autodiff(fasten_main, Const,
        Val(convert(Int, params.wgsize)),
        Duplicated(deck.protein, d_protein),
        Const(deck.ligand),
        Const(deck.forcefield),
        Const(deck.poses),
        Duplicated(etotals, d_etotals)
      )
    end
  else
    elapsed = @elapsed for _ = 1:params.iterations
      fasten_main(
        Val(convert(Int, params.wgsize)),
        deck.protein,
        deck.ligand,
        deck.forcefield,
        deck.poses,
        etotals,
      )
    end
  end

  (etotals, elapsed, params.wgsize)
end

@eval @fastmath function fasten_main(
  ::Val{WGSIZE},
  protein::Array{Atom},
  ligand::Array{Atom},
  forcefield::Array{FFParams},
  poses::Array{Float32,2},
  etotals::Array{Float32}, # Only output
) where {WGSIZE}
  nposes::Int = size(poses)[2]
  numGroups::Int = nposes ÷ WGSIZE

  Enzyme.@parallel protein ligand forcefield poses etotals for group = 1:numGroups
    protein = Base.Experimental.Const(protein)
    ligand = Base.Experimental.Const(ligand)
    forcefield = Base.Experimental.Const(forcefield)
    poses = Base.Experimental.Const(poses)
    nligand::Int = length(ligand)
    nprotein::Int = length(protein)
    etot = MArray{Tuple{WGSIZE}, Float32}(undef)
    transform = MArray{Tuple{WGSIZE, 3, 4},Float32}(undef)

    @aliasscope @simd ivdep for i = 1:WGSIZE
    # for i = 1:WGSIZE
      ix = (group - 1) * (WGSIZE) + i
      @inbounds sx::Float32 = sin(poses[1, ix])
      @inbounds cx::Float32 = cos(poses[1, ix])
      @inbounds sy::Float32 = sin(poses[2, ix])
      @inbounds cy::Float32 = cos(poses[2, ix])
      @inbounds sz::Float32 = sin(poses[3, ix])
      @inbounds cz::Float32 = cos(poses[3, ix])
      @inbounds transform[i, 1, 1] = cy * cz
      @inbounds transform[i, 1, 2] = sx * sy * cz - cx * sz
      @inbounds transform[i, 1, 3] = cx * sy * cz + sx * sz
      @inbounds transform[i, 1, 4] = poses[4, ix]
      @inbounds transform[i, 2, 1] = cy * sz
      @inbounds transform[i, 2, 2] = sx * sy * sz + cx * cz
      @inbounds transform[i, 2, 3] = cx * sy * sz - sx * cz
      @inbounds transform[i, 2, 4] = poses[5, ix]
      @inbounds transform[i, 3, 1] = -sy
      @inbounds transform[i, 3, 2] = sx * cy
      @inbounds transform[i, 3, 3] = cx * cy
      @inbounds transform[i, 3, 4] = poses[6, ix]
      @inbounds etot[i] = Zero
      # $(Expr(:loopinfo, Symbol("julia.simdloop"), (Symbol("llvm.loop.unroll.full"),)))
    end

    for il = 1:nligand
      @inbounds l_atom::Atom = ligand[il]

      @inbounds l_params::FFParams = forcefield[l_atom.type+1]
      lhphb_ltz::Bool = l_params.hphb < Zero
      lhphb_gtz::Bool = l_params.hphb > Zero

      lpos = MArray{Tuple{WGSIZE, 3}, Float32}(undef)

      @simd ivdep for i = 1:WGSIZE
      # for i = 1:WGSIZE
        @inbounds lpos[i, 1] = (
          transform[i, 1, 4] +
          l_atom.x * transform[i, 1, 1] +
          l_atom.y * transform[i, 1, 2] +
          l_atom.z * transform[i, 1, 3]
        )
        @inbounds lpos[i, 2] = (
          transform[i, 2, 4] +
          l_atom.x * transform[i, 2, 1] +
          l_atom.y * transform[i, 2, 2] +
          l_atom.z * transform[i, 2, 3]
        )
        @inbounds lpos[i, 3] = (
          transform[i, 3, 4] +
          l_atom.x * transform[i, 3, 1] +
          l_atom.y * transform[i, 3, 2] +
          l_atom.z * transform[i, 3, 3]
        )
        # $(Expr(:loopinfo, Symbol("julia.simdloop"), (Symbol("llvm.loop.unroll.full"),)))
      end

      @aliasscope for ip = 1:nprotein
        @inbounds p_atom = protein[ip]
        @inbounds p_params::FFParams = forcefield[p_atom.type+1]

        radij::Float32 = p_params.radius + l_params.radius
        r_radij::Float32 = One / radij

        elcdst::Float32 = (p_params.hbtype == HbtypeF && l_params.hbtype == HbtypeF) ? Four : Two
        elcdst1::Float32 =
          (p_params.hbtype == HbtypeF && l_params.hbtype == HbtypeF) ? Quarter : Half
        type_E::Bool = ((p_params.hbtype == HbtypeE || l_params.hbtype == HbtypeE))

        phphb_ltz::Bool = p_params.hphb < Zero
        phphb_gtz::Bool = p_params.hphb > Zero
        phphb_nz::Bool = p_params.hphb != Zero
        p_hphb::Float32 = p_params.hphb * (phphb_ltz && lhphb_gtz ? -One : One)
        l_hphb::Float32 = l_params.hphb * (phphb_gtz && lhphb_ltz ? -One : One)
        distdslv::Float32 =
          (phphb_ltz ? (lhphb_ltz ? Npnpdist : Nppdist) : (lhphb_ltz ? Nppdist : -Float32Max))
        r_distdslv::Float32 = One / distdslv

        chrg_init::Float32 = l_params.elsc * p_params.elsc
        dslv_init::Float32 = p_hphb + l_hphb

        @simd ivdep for i = 1:WGSIZE
        # for i = 1:WGSIZE
          @inbounds x::Float32 = lpos[i, 1] - p_atom.x
          @inbounds y::Float32 = lpos[i, 2] - p_atom.y
          @inbounds z::Float32 = lpos[i, 3] - p_atom.z

          distij::Float32 = sqrt(x ^ 2 + y ^ 2+ z ^ 2)

          # Calculate the sum of the sphere radii
          distbb::Float32 = distij - radij
          zone1::Bool = (distbb < Zero)

          # Calculate steric energy
          @inbounds etot[i] += (One - (distij * r_radij)) * (zone1 ? Two * Hardness : Zero)

          # Calculate formal and dipole charge interactions
          chrg_e::Float32 =
            chrg_init * ((zone1 ? One : (One - distbb * elcdst1)) * (distbb < elcdst ? One : Zero))
          neg_chrg_e::Float32 = -abs(chrg_e)
          chrg_e = type_E ? neg_chrg_e : chrg_e
          @inbounds etot[i] += chrg_e * Cnstnt

          # Calculate the two cases for Nonpolar-Polar repulsive interactions
          coeff::Float32 = (One - (distbb * r_distdslv))
          dslv_e::Float32 = dslv_init * ((distbb < distdslv && phphb_nz) ? One : Zero)
          dslv_e *= (zone1 ? One : coeff)
          @inbounds etot[i] += dslv_e
          # $(Expr(:loopinfo, Symbol("julia.simdloop"), (Symbol("llvm.loop.unroll.full"),)))
        end
      end
    end
    @simd ivdep for i = 1:WGSIZE
    # for i = 1:WGSIZE
      ix = (group - 1) * WGSIZE + i
      @inbounds etotals[ix] = etot[i] * Half
      # $(Expr(:loopinfo, Symbol("julia.simdloop"), (Symbol("llvm.loop.unroll.full"),)))
    end
  end
end

main()

max_stars_repo_name: wsmoses/Enzyme-BUDE

max_issues_repo_name: wsmoses/Enzyme-BUDE

max_forks_repo_name: wsmoses/Enzyme-BUDE


text:
module SteadyState

using QuantumOptics
using IterativeSolvers, LinearMaps, LinearAlgebra, SparseArrays

const T_blas = Union{Float64,Float32,ComplexF64,ComplexF32}
const DecayRates = Union{Vector{T} where T<:Number, Matrix{T} where T<:Number, Nothing}

isblascompatible(M::Array{T,2}) where {T<:T_blas} = true
isblascompatible(M::AbstractArray{T,2}) where T = false

include("interface.jl")
include("generic_method.jl")

end # module


text:
# # Tutorial
#
# First, let us load the relevant packages and an example dataset:
#

using Statistics, StatsBase, StatsPlots, JuliaDB
using OnlineStats: Mean, Variance
using Recombinase
using Recombinase: cumulative, density, hazard, prediction

data = loadtable(joinpath(Recombinase.datafolder, "school.csv"))

#
# ### Simple scatter plots
#
# Then we can compute a simple scatter plot of one variable against an other. This is done in two steps: first the positional and named arguments of the plot call are computed, then they are passed to a plotting function:
#

args, kwargs = series2D(
    data,
    Group(:Sx),
    select = (:MAch, :SSS),
    )
scatter(args...; kwargs...)

# This creates an overcrowded plot. We could instead compute the average value of our columns of interest for each school and then plot just one point per school (with error bars representing variability within the school):
#

args, kwargs = series2D(
    data,
    Group(:Sx),
    error = :School,
    select = (:MAch, :SSS),
    )
scatter(args...; kwargs...)

# By default, this computes the mean and standard error, we can pass `stats = Mean()` to only compute the mean.
#
# ### Splitting by many variables
#
# We can use different attributes to split the data as follows:
#

args, kwargs = series2D(
    data,
    Group(color = :Sx, markershape = :Sector),
    error = :School,
    select = (:MAch, :SSS),
    stats = Mean(),
    )
scatter(args...; kwargs...)

# ### Styling the plot
#
# There are two ways in which we can style the plot: first, we can pass a custom set of colors instead of the default palette:
#

args, kwargs = series2D(
    data,
    Group(:Sx),
    error = :School,
    select = (:MAch, :SSS),
    stats = Mean(),
    color = [:red, :blue]
    )
scatter(args...; kwargs...)

# Second, we can style plat attributes as we would normally do:
#

args, kwargs = series2D(
    data,
    Group(:Sx),
    error = :School,
    select = (:MAch, :SSS),
    stats = Mean(),
    )
scatter(args...; legend = :topleft, markersize = 10, kwargs...)

# ### Computing summaries
#
# It is also possible to get average value and variability of a given analysis (density, cumulative, hazard rate and local regression are supported so far, but one can also add their own function) across groups.
#
# For example (here we use `ribbon` to signal we want a shaded ribbon to denote the error estimate):
#

args, kwargs = series2D(
    cumulative,
    data,
    Group(:Sx),
    error = :School,
    select = :MAch,
    ribbon = true
   )
plot(args...; kwargs..., legend = :topleft)

# Note that extra keyword arguments can be passed to the analysis:
#

args, kwargs = series2D(
    density(bandwidth = 1),
    data,
    Group(color=:Sx, linestyle=:Sector),
    error = :School,
    select = :MAch,
    ribbon = true
   )
plot(args...; kwargs..., legend = :bottom)

# If we do not specify `error`, it defaults to the "analyses specific error". For discrete prediction it is the standard error of the mean across observations.
#

args, kwargs = series2D(
    prediction,
    data,
    Group(color = :Minrty),
    select = (:Sx, :MAch),
)
groupedbar(args...; kwargs...)

# ### Axis style selection
#
# Analyses try to infer the axis type (continuous if the variable is numeric, categorical otherwise). If that is not appropriate for your data you can use `discrete(prediction)` or `continuous(prediction)` (works for `hazard`, `density` and `cumulative` as well).
# A special type of axis type is vectorial: the `x` and `y` axes can be contain `AbstractArray` elements, in which case
# we take views of elements of `y` corresponding to elements of `x`. This can be useful to compute averages of time varying signals.

using Recombinase: offsetrange
signal = sin.(range(0, 10π, length = 1000)) .+ 0.1 .* rand.()
events = range(50, 950, step = 100)
z = repeat([true, false], outer = 5)
x = [offsetrange(signal, ev) for ev in events]
y = fill(signal, length(x))
t = table(x, y, z, names = [:offsets, :signals, :peak])

args, kwargs = series2D(
    prediction(axis = -60:60),
    t,
    Group(:peak),
    select = (:offsets, :signals),
    ribbon = true,
)
plot(args...; kwargs...)

# ### Post processing
#
# Finally, for some analyses it can be useful to postprocess the result. For example, in the case
# of the "signal plot" above, we may wish to rescale the `x` axis. This is done by passing a named
# tuple of functions as a `postprocess` keyword argument, which will be applied element-wise to the relative column of the output. For example, if our signal was sampled at `60 Hz`, we may wish to divide the `:offsets` column by `60` to show the result in seconds:

args, kwargs = series2D(
    prediction(axis = -60:60),
    t,
    Group(:peak),
    select = (:offsets, :signals),
    ribbon = true,
    postprocess = (; offsets = t -> t/60),
)
plot(args...; kwargs...)


text:
using Plots, GLVisualize, GeometryTypes, Colors

pl_size = if !isempty(GLVisualize.get_screens())
    widths(GLVisualize.current_screen())
else
    (1000, 1000)
end
glvisualize(size = pl_size)

description = """
How to create different 3D plots with Plots.jl
"""

function lorenz(t0, a, b, c, h)
    Point3f0(
        t0[1] + h * a * (t0[2] - t0[1]),
        t0[2] + h * (t0[1] * (b - t0[3]) - t0[2]),
        t0[3] + h * (t0[1] * t0[2] - c * t0[3]),
    )
end
# step through the `time`
function lorenz(array::Vector, a=5.0,b=2.0,c=6.0,d=0.01)
    t0 = Point3f0(0.1, 0, 0)
    for i = eachindex(array)
        t0 = lorenz(t0, a,b,c,d)
        array[i] = t0
    end
    array
end

N = 1000
xyz = lorenz(zeros(Point3f0, N), 26., 28., 9.0, 0.01)
x,y,z = map(first, xyz), map(x-> x[2], xyz), map(last, xyz)
p1 = plot(x,y,z,
    line=0.5,
    marker=(:circle, 0.4), ms = linspace(2.0, 10.0, N),
    markercolor = colormap("RdBu", N),
    leg = false,
    grid = false
)

dphi, dtheta = pi/200.0f0, pi/200.0f0
function mgrid(dim1, dim2)
    X = [i for i in dim1, j in dim2]
    Y = [j for i in dim1, j in dim2]
    return X,Y
end

phi,theta = mgrid(0f0:dphi:(pi+dphi*1.5f0), 0f0:dtheta:(2f0*pi+dtheta*1.5f0));
m0 = 4f0; m1 = 3f0; m2 = 2f0; m3 = 3f0; m4 = 6f0; m5 = 2f0; m6 = 6f0; m7 = 4f0;
a = sin.(m0*phi).^m1;
b = cos.(m2*phi).^m3;
c = sin.(m4*theta).^m5;
d = cos.(m6*theta).^m7;
r = a + b + c + d;
x = r.*sin.(phi).*cos.(theta);
y = r.*cos.(phi);
z = r.*sin.(phi).*sin.(theta);
surface1 = Plots.surface(x,y,z)

function xy_data(x, y, i, N)
    x = ((x / N)-0.5) * i
    y = ((y / N)-0.5) * i
    r = sqrt(x * x + y * y)
    sin(r) / r
end

surface2 = Plots.surface(Float32[xy_data(x, y, 20f0, N) for x=1:100, y=1:100])
surface3 = Plots.wireframe(
    Float32[xy_data(x * 10, y * 10, 40f0, N) for x=1:30, y=1:30],
    line = 4.0
)

p = plot(surface1, surface2, surface3, p1, show = true)

gui()


text:
# Various constants used to determine problem type

const RASTER = ["raster", "Raster"]
const PAIRWISE = ["pairwise", "Pairwise"]
const ADVANCED = ["advanced", "Advanced"]
const ONETOALL = ["one-to-all", "one_to_all"]
const ALLTOONE = ["all-to-one", "all_to_one"]
const SINGLE = ["single", "Single"]
const DOUBLE = ["double", "Double"]

# Solver constants
const AMG = ["cg+amg", "amg+cg"]
const CHOLMOD = ["cholmod", "cholesky", "cholfact"]
const MKLPARDISO = ["mklpardiso", "MKLPardiso", "MKLPARDISO", "pardiso"]

# Constants used in IO
const AAGRID = 2
const TXTLIST = 3
const PAIRS_AAGRID = 4
const PAIRS_LIST = 5
const TRUELIST = ["True", "true", "1"]

const FILE_TYPE_NPY = 1
const FILE_TYPE_AAGRID = 2
const FILE_TYPE_TXTLIST = 3
const FILE_TYPE_INCL_PAIRS_AAGRID = 4
const FILE_TYPE_INCL_PAIRS = 5
const FILE_TYPE_GEOTIFF = 6

const FILE_HDR_GZIP = "\x1f\x8b\x08"
const FILE_HDR_NPY = "\x93NUMPY"
const FILE_HDR_AAGRID = "ncols"
const FILE_HDR_INCL_PAIRS_AAGRID = "min"
const FILE_HDR_INCL_PAIRS = "mode"

# Constants for logging
const NONE = ["NONE", "None", "none"]
const INFO = ["INFO", "info", "Info"]
const DEBUG = ["DEBUG", "debug", "Debug"]
const WARNING = ["WARNING", "WARN", "warning", "warn", "Warning", "Warn"]
const CRITICAL = ["CRITICAL", "critical", "Critical"]

# Constants for tolerance when checking matrix norms
const TOL_SINGLE = 1e-3
const TOL_DOUBLE = 1e-5


text:
s = open("input.txt") do file
    readlines(file)
end

function CountTrees(lines, deltax, deltay)
    treeCount = 0;

    repeat = length(lines[1]);
    x = 1 + deltax;
    i = 1 + deltay;
    while i <= length(lines)
        if lines[i][x] == '#'
            treeCount += 1;
        end
        x = ((x + deltax - 1) % (repeat)) + 1;
        i += deltay;
    end

    return treeCount;
end

mult = 1;
mult *= CountTrees(s, 1, 1);
mult *= CountTrees(s, 3, 1);
mult *= CountTrees(s, 5, 1);
mult *= CountTrees(s, 7, 1);
mult *= CountTrees(s, 1, 2);
mult


text:
import StatsBase: stderror, sample

"""

QuapResult

$(FIELDS)

where:

coef                            : NamedTuple with parameter MAP estimates,
vcov                            : Covariance matric,
converged                       : Simple check that multiple chains converged,
distr                           : Distributike to sample from (Normal or MvNormal),
params                          : Vector of parameter symbols.

"""
struct QuapResult{
       N <: NamedTuple
} <: StatsBase.StatisticalModel
       coef :: N
       vcov :: Array{Float64, 2}
       converged :: Bool
       distr :: Union{Normal, MvNormal}
       params :: Array{Symbol, 1}
end

"""

Compute the quadratic approximation to the posterior distribution.

$(SIGNATURES)

### Required arguments
```julia
* `name::String`                : Name for SampleModel
* `model::String`               : Stan Language model
```

### Keyword arguments
```julia
* `data`                        : Data for model (NamedTuple or Duct)
* `init`                        : Initial values for parameters (NamedTuple or Dict)
```

### Reyrns
```julia
* `res::QuapResult`              : Returned object
```

In general using `init` results in better behavior.

"""
function stan_quap(
    name::AbstractString,
    model::AbstractString;
    kwargs...)

    om = OptimizeModel(name, model)
    rc = stan_optimize(om; kwargs...)

    if success(rc)
        tmp, cnames = read_optimize(om)
        optim = Dict()
        for key in keys(tmp)
            if !(key in ["lp__", :stan_version])
                optim[Symbol(key)] = tmp[key]
            end
        end
        sm = SampleModel(name, model, om.tmpdir)
        rc2 = stan_sample(sm; kwargs...)
    else
        return ((nothing, nothing, nothing))
    end

    if success(rc2)
        qm = quap(sm, optim, cnames)
        return((qm, sm, (optim=optim, cnames=cnames)))
    else
        return((nothing, sm, nothing))
    end
end

"""

Compute the quadratic approximation to the posterior distribution.

$(SIGNATURES)

Not exported
"""
function quap(
    sm_sam::SampleModel,
    optim::Dict,
    cnames::Vector{String})

    samples = read_samples(sm_sam, :dataframe)

    n = Symbol.(names(samples))
    coefnames = tuple(n...,)
    c = [optim[Symbol(coefname)][1] for coefname in coefnames]
    cvals = reshape(c, 1, length(n))
    coefvalues = reshape(c, length(n))
    v = Statistics.covm(Array(samples[:, n]), cvals)

    distr = if length(coefnames) == 1
        Normal(coefvalues[1], √v[1])  # Normal expects stddev
    else
        MvNormal(coefvalues, v)       # MvNormal expects variance matrix
    end

    converged = true
    for coefname in coefnames
        o = optim[Symbol(coefname)]
        converged = abs(sum(o) - 4 * o[1]) < 0.001 * abs(o[1])
        !converged && break
    end

    ntcoef = namedtuple(coefnames, coefvalues)
    return(QuapResult(ntcoef, v, converged, distr, n))
end

"""

Sample from a quadratic approximation to the posterior distribution.

$(SIGNATURES)

### Required arguments
```julia
* `qm`                          : QuapResult object (see: `?QuapResult`)
```

### Keyword arguments
```julia
* `nsamples = 4000`             : Number of smaples taken from distribution
```

"""
function sample(qm::QuapResult; nsamples=4000)
  df = DataFrame()
  p = Particles(nsamples, qm.distr)
  for (indx, coef) in enumerate(qm.params)
    if length(qm.params) == 1
      df[!, coef] = p.particles
    else
      df[!, coef] = p[indx].particles
    end
  end
  df
end

# TEMPORARILY:

# Used by Max. Need to check if this works in general
function sample(qr::QuapResult, count::Int)::DataFrame
    names = qr.params                 # StatsBase.coefnames(mode_result) in Turing
    means = values(qr.coef)           # StatsBase.coef(mode_result) in Turing
    sigmas = Diagonal(qr.vcov)        # StatsBase.stderr(mode_result) in Turing

    DataFrame([
        name => rand(Normal(μ, σ), count)
            for (name, μ, σ) ∈ zip(names, means, sigmas)
    ])
end

# Will be deprecated, use QuapResult object
function sample(qm::NamedTuple; nsamples=4000)
  df = DataFrame()
  p = Particles(nsamples, qm.distr)
  for (indx, coef) in enumerate(qm.params)
    if length(qm.params) == 1
      df[!, coef] = p.particles
    else
      df[!, coef] = p[indx].particles
    end
  end
  df
end

export
    QuapResult,
    stan_quap,
    sample


text:
using TrajOptPlots
using Test
using StaticArrays
# using RecipesBase
using Plots
using RobotDynamics
using TrajectoryOptimization
using RobotZoo
# using Altro
using MeshCat
using Rotations
using LinearAlgebra
const TO = TrajectoryOptimization

include("problems.jl")

@testset "2D Plots" begin
    include("2d_plots.jl")
end

@testset "3D Visualization" begin
    include("3d_vis.jl")
    include("models.jl")
end


text:
# ------------------------------------------------------------------------------------------
# ## Getting started with Jupyter notebooks
#
# ### Running a cell
# To execute code within a cell, select that cell and either (1) hit `Shift` and `Enter` or
# (2) hit the run button (the right pointing arrow) above.
# ------------------------------------------------------------------------------------------

1 + 1
2 + 2

# ------------------------------------------------------------------------------------------
# If you're new to jupyter notebooks, note that only the last line of a cell prints by
# default when you execute that cell and that you can suppress this output with a semicolon
# ------------------------------------------------------------------------------------------

1 + 1
2 + 2;

# ------------------------------------------------------------------------------------------
# ### How to get docs for Julia functions
#
# To get docs for a function you're not familiar with, precede it with a question mark.
# (This works in the terminal too!)
# ------------------------------------------------------------------------------------------

?println

# ------------------------------------------------------------------------------------------
# ### How to use shell commands
#
# Type `;` and then you can use shell commands. For example,
# ------------------------------------------------------------------------------------------

;ls

;pwd

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
"""
```
label_components(ContourTracing(), binary_image)
```

Computes the connected components of a binary image using 8-connectivity.

# Output

Returns a two-dimensional `Int` array that represents the label of every pixel.

# Example

Computes the labels for the "blobs" image in the `TestImages` package.

```julia
using Images, TestImages, ImageBinarization, ImageComponentAnalysis

image = testimage("blobs")
binary_image = binarize(MinimumError(), image)

labels = label_components(ContourTracing(), binary_image)
```

# Reference

[1] F. Chang, C.-J. Chen, and C.-J. Lu, “A linear-time component-labeling algorithm using contour tracing technique,” Computer Vision and Image Understanding, vol. 93, no. 2, pp. 206–220, Feb. 2004. [doi:10.1016/j.cviu.2003.09.002](https://doi.org/10.1016/j.cviu.2003.09.002)
"""
function label_components(algorithm::ContourTracing, binary_image::AbstractArray)

  height, width = size(binary_image)
  add_dummy_row = false

  # Check if all of the pixels are black in the first row.
  for col = 1:width
    if (add_dummy_row = binary_image[1, col] == 1) break end
  end

  # Add a dummy row on top if there is at least one white pixel in the first row.
  if add_dummy_row
    dummy_row = zeros(eltype(binary_image), (1, width))
    binary_image = vcat(dummy_row, binary_image)
    height += 1
  end

  labels = zeros(Int, (height, width))
  labelindex = 1

  # Assign label and perform contour checks for each white pixel.
  for row = 1:height
    for col = 1:width
      if (binary_image[row, col] == 1)

        # Assign current label, perform contour tracing and increment
        # the current label by one if it is an external contour.
        if (external_contour_check(row, col, height, width, binary_image, labels))
          labels[row, col] = labelindex
          external_contour_tracing(row, col, height, width, binary_image, labels, labelindex)
          labelindex += 1
        end

        # Perform contour tracing if it is an internal contour.
        if (internal_contour_check(row, col, height, width, binary_image))
          # Assign the west neighbour's label if the current pixel is not an external contour.
          if (labels[row, col] == 0)
            assign_west_label(row, col, height, width, labels)
          end
          internal_contour_tracing(row, col, height, width, binary_image, labels)
        end

        # Assign the west neighbour's label if the current pixel is not a contour point.
        if (labels[row, col] == 0)
          assign_west_label(row, col, height, width, labels)
        end

      end
    end
  end

  # Remove the dummy row if it was added.
  if add_dummy_row
    labels = labels[2:end, :]
  end

  labels
end

# Perform contour tracing from northeast (index 7) with current label index.
function external_contour_tracing(row, col, height, width, binary_image, labels, labelindex)
  contour_tracing(7, row, col, height, width, binary_image, labels, labelindex)
end

# Check if the pixel is unlabelled and the north neighbour is black and in bounds.
function external_contour_check(row, col, height, width, binary_image, labels)
  north_neighbour = get_north_neighbour(row, col, height, width)
  if (labels[row, col] != 0) return false end
  return north_neighbour != (0, 0) && binary_image[north_neighbour...] == 0
end

# Perform contour tracing from southwest (index 3) with current pixel's label.
function internal_contour_tracing(row, col, height, width, binary_image, labels)
  contour_tracing(3, row, col, height, width, binary_image, labels, labels[row, col])
end

# Check if south neighbour is black and in bounds.
function internal_contour_check(row, col, height, width, binary_image)
  south_neighbour = get_south_neighbour(row, col, height, width)
  return south_neighbour != (0, 0) && binary_image[south_neighbour...] == 0
end

# Assign the west neighbour's label if the current pixel is unlabelled.
function assign_west_label(row, col, height, width, labels)
  if ((west_neighbour = get_west_neighbour(row, col, height, width)) != (0, 0))
    labels[row, col] = labels[west_neighbour...]
  end
end

# Get the neighbour to the north (index 6).
function get_north_neighbour(row, col, height, width)
  north_neighbour = get_neighbour(6, row, col, height, width)
end

# Get the neighbour to the south (index 2).
function get_south_neighbour(row, col, height, width)
  south_neighbour = get_neighbour(2, row, col, height, width)
end

# Get the neighbour to the west (index 4).
function get_west_neighbour(row, col, height, width)
  west_neighbour = get_neighbour(4, row, col, height, width)
end

# Return a neighbour by the index or (0, 0) if out of bounds.
function get_neighbour(index, row, col, height, width)

  if (1 <= index <= 3) row += 1 end
  if (3 <= index <= 5) col -= 1 end
  if (5 <= index <= 7) row -= 1 end
  if (index == 7|| 0 <= index <= 1) col += 1 end

  inbounds(row, col, height, width) ? (row, col) : (0, 0)
end

# Check if the current pixel is in bounds.
function inbounds(row, col, height, width)
  return 1 <= row <= height && 1 <= col <= width
end

# Trace an external or internal contour.
function contour_tracing(search_index, row, col, height, width, binary_image, labels, labelindex)

  initial_position = (row, col)
  previous_index = 0

  while (row, col) != (0, 0)
    row, col, previous_index = tracer(search_index, row, col, height, width, binary_image, labels, labelindex)
    # Stop tracing if the tracer has returned to the starting point.
    if ((row, col) == initial_position) break end

    # Increase search index by two around the pixel.
    search_index = (previous_index + 2) % 8
  end
end

# Search the pixel's neighbours and assign labels starting from search_index.
function tracer(search_index, row, col, height, width, binary_image, labels, labelindex)

  for search_iteration = 1:8

    neighbour = get_neighbour(search_index, row, col, height, width)

    # Look for the next neighbour if the current neighbour is out of bounds.
    if (neighbour == (0, 0))
      # Increase search index by one around the pixel.
      search_index = (search_index + 1) % 8
      continue
    end

    # Assign label if the neighbour is white and -1 if it is black.
    if (binary_image[neighbour...] == 1)
      labels[neighbour...] = labelindex
      # Return the neighbour's position and the index relative to the current pixel.
      return (neighbour..., (search_index + 4) % 8)
    else
      labels[neighbour...] = -1
      # Increase search index by one around the pixel.
      search_index = (search_index + 1) % 8
    end
  end
  # Return zeros if the pixel is an isolated point.
  (0, 0, 0)
end


text:
struct NNPool
    nn_pools::Dict{Symbol, Dict{Int, Int}}
    pools::Vector{Symbol} # TODO: use labelmap

    nn::Vector{Vector{Int}}
    reverse_nn::Vector{Vector{Int}}
    symmetric_nn::Vector{Vector{Int}}
end

function NNPool(data; k = 20)
    nn_pools = Dict{Symbol, Dict{Int, Int}}()
    nn_pools[:Lin] = Dict{Int, Int}()
    nn_pools[:Lout] = Dict{Int, Int}()

    kdtree = KDTree(data; leafsize = k)
    nn = knn(kdtree, data, k, true)[1]
    reverse_nn = reverseKNN(nn)
    symmetric_nn = [n ∩ r for (n,r) in zip(nn, reverse_nn)]

    NNPool(nn_pools, fill(:U, size(data,2)), nn, reverse_nn, symmetric_nn)
end

function Base.push!(nnpool::NNPool, id::Int, pool_label::Symbol)
    if haskey(nnpool.nn_pools[pool_label], id)
        nnpool.nn_pools[pool_label][id] = nnpool.nn_pools[pool_label][id] + 1
    else
        nnpool.nn_pools[pool_label][id] = 1
    end
end

function add_to_pools!(nnpool::NNPool, id, label::Symbol, k::Int)
    pool_label = ifelse(label == :inlier, :Lin, :Lout)

    nnpool.pools[id] = pool_label

    if label == :inlier
        neighbors = nnpool.nn[id][1:k]
    else
        snn_elements = min(length(nnpool.symmetric_nn[id]), k)
        neighbors = nnpool.symmetric_nn[id][1:snn_elements]
    end
    for n in neighbors
        push!(nnpool, n, pool_label)
    end
end

function reverseKNN(nn)
    map(x -> findall(in.(x, nn)), eachindex(nn))
end

max_stars_repo_name: holtri/lama

max_issues_repo_name: holtri/lama

max_forks_repo_name: holtri/lama


text:
@testset "Fitness Complexity algorithm tests" begin

	# Test fitness algorithm --------------------------
	# Initialize test objects
	INCIDENCE_MAT = [0 1 0 1 1 0 0 0 0 0;
			 0 1 1 0 0 0 1 0 0 0;
			 1 1 0 1 1 0 0 0 1 1;
			 1 0 0 1 0 1 0 0 1 1;
			 0 0 1 1 0 0 0 1 0 0]

	INCIDENCE_MAT = convert(Array{Float64, 2}, INCIDENCE_MAT)
	ITS = 10

	countries = ["c$x" for x in 1:size(INCIDENCE_MAT, 1)]
	products = ["p$x" for x in 1:size(INCIDENCE_MAT, 2)]

	# Check inputs:
	@test_throws AssertionError get_fitness(INCIDENCE_MAT, iterations=ITS, product_codes = products, country_codes = countries[1:4]) # Wrong length of country codes
	@test_throws AssertionError get_fitness(INCIDENCE_MAT, iterations=ITS, product_codes = products[1:5], country_codes = countries) # Wrong length of product codes
	@test_throws MethodError get_fitness(INCIDENCE_MAT, iterations=ITS, product_codes = [i for i in 1:10], country_codes = countries) # Non-string product codes
	@test_throws MethodError get_fitness(INCIDENCE_MAT, iterations=ITS, product_codes = products, country_codes = [i for i in 1:5]) # Non-string country codes
	@test_throws MethodError get_fitness(DataFrame(INCIDENCE_MAT), iterations=ITS, product_codes = products, country_codes = countries) # non-matrix input
	@test_throws MethodError get_fitness(INCIDENCE_MAT, iterations=4.7, product_codes = products, country_codes = countries) # non-integer iterations

	# Check output:
	# Output both Q and F
	@test all([i ∈ keys(get_fitness(INCIDENCE_MAT, iterations = ITS, product_codes = products, country_codes = countries)) for i in ["Q", "F"]])
	# Test presence of product- and country codes
	@test get_fitness(INCIDENCE_MAT, iterations = ITS, product_codes = products, country_codes = countries)["Q"].product_code == products
	@test get_fitness(INCIDENCE_MAT, iterations = ITS, product_codes = products, country_codes = countries)["F"].country_code == countries

	# Check results:
	# Create toy results -----------
	# Define initial conditions:
	F_current = ones(size(INCIDENCE_MAT, 1), 1)
	Q_current = ones(size(INCIDENCE_MAT, 2), 1)

	# Run algorithm
	for i in 1:ITS
		F = INCIDENCE_MAT * Q_current
		Q = 1 ./ (INCIDENCE_MAT' * (1 ./ F_current))
		F_current = F / (sum(F) / size(F, 1))
		Q_current = Q / (sum(Q) / size(Q, 1))
	end

	complexity_tbl = DataFrames.DataFrame(product_code = products, complexity = [Q_current...])
	fitness_tbl = DataFrames.DataFrame(country_code = countries, fitness = [F_current...])

	@test get_fitness(INCIDENCE_MAT, iterations = ITS, product_codes = products, country_codes = countries)["Q"].complexity == complexity_tbl.complexity
	@test get_fitness(INCIDENCE_MAT, iterations = ITS, product_codes = products, country_codes = countries)["F"].fitness == fitness_tbl.fitness
	@test get_fitness(INCIDENCE_MAT, iterations = ITS, product_codes = products, country_codes = countries)["Q"].product_code == complexity_tbl.product_code
	@test get_fitness(INCIDENCE_MAT, iterations = ITS, product_codes = products, country_codes = countries)["F"].country_code == fitness_tbl.country_code

end

max_stars_repo_name: sbpost/compost

max_issues_repo_name: sbpost/compost

max_forks_repo_name: sbpost/compost