text:
#=
This is an integration test using Julia's pmap parallelization.

This is an example on training the TurbulenceConvection.jl implementation
of the EDMF scheme with data generated using PyCLES (an LES solver) or
TurbulenceConvection.jl (perfect model setting).

This example parallelizes forward model evaluations using Julia's pmap() function.
This parallelization strategy is efficient for small ensembles (N_ens < 20). For
calibration processes involving a larger number of ensemble members, parallelization
using HPC resources directly is advantageous.
=#

# Import modules to all processes
using Distributed
@everywhere using Pkg
@everywhere Pkg.activate(@__DIR__)
@everywhere begin
    using ArgParse
    import Random
    Random.seed!(1234)
    using CalibrateEDMF
    using CalibrateEDMF.DistributionUtils
    using CalibrateEDMF.Pipeline
    using CalibrateEDMF.HelperFuncs
    # Import EKP modules
    using EnsembleKalmanProcesses
    using EnsembleKalmanProcesses.ParameterDistributions
end
using JLD2
using Test

s = ArgParseSettings()
@add_arg_table s begin
    "--config"
    help = "Inverse problem config file path"
    arg_type = String
end
parsed_args = parse_args(ARGS, s)
config_rel_filepath = parsed_args["config"]
config_basename = first(split(config_rel_filepath, "."))

config_filename = joinpath(@__DIR__, config_rel_filepath)
include(config_filename)
config = get_config()
algorithm = config["process"]["algorithm"]

# Initialize calibration process
outdir_path = init_calibration(config; config_path = config_filename, mode = "pmap")
priors = deserialize_prior(load(joinpath(outdir_path, "prior.jld2")))

# Dispatch SCM eval functions to workers
@everywhere begin
    scm_eval_train(version) = versioned_model_eval(version, $outdir_path, "train", $config)
    scm_eval_validation(version) = versioned_model_eval(version, $outdir_path, "validation", $config)
end
# Calibration process
N_iter = config["process"]["N_iter"]
@info "Running EK updates for $N_iter iterations"
for iteration in 1:N_iter
    @time begin
        @info "   iter = $iteration"
        versions = readlines(joinpath(outdir_path, "versions_$(iteration).txt"))
        ekp = load(ekobj_path(outdir_path, iteration))["ekp"]
        pmap(scm_eval_train, versions)
        pmap(scm_eval_validation, versions)
        ek_update(ekp, priors, iteration, config, versions, outdir_path)
    end
end

import NCDatasets
const NC = NCDatasets

nt = NC.Dataset(joinpath(outdir_path, "Diagnostics.nc"), "r") do ds
    (;
        mse_full_mean = Array(ds.group["metrics"]["mse_full_mean"]),
        mse_full_nn_mean = Array(ds.group["metrics"]["mse_full_nn_mean"]),
        mse_full_var = Array(ds.group["metrics"]["mse_full_var"]),
        val_mse_full_mean = Array(ds.group["metrics"]["val_mse_full_mean"]),
        val_mse_full_nn_mean = Array(ds.group["metrics"]["val_mse_full_nn_mean"]),
        val_mse_full_var = Array(ds.group["metrics"]["val_mse_full_var"]),
        phi_mean = Array(ds.group["ensemble_diags"]["phi_mean"]),
    )
end

mse_full_mean = nt.mse_full_mean[1:(end - 1)]
mse_full_nn_mean = nt.mse_full_nn_mean[1:(end - 1)]
mse_full_var = nt.mse_full_var[1:(end - 1)]
val_mse_full_mean = nt.val_mse_full_mean[1:(end - 1)]
val_mse_full_nn_mean = nt.val_mse_full_nn_mean[1:(end - 1)]
val_mse_full_var = nt.val_mse_full_var[1:(end - 1)]
phi_mean = nt.phi_mean

import Plots

# https://github.com/jheinen/GR.jl/issues/278#issuecomment-587090846
ENV["GKSwstype"] = "nul"

iters = 1:N_iter

@show iters
@show mse_full_mean
@show mse_full_nn_mean
@show mse_full_var
@show phi_mean

mse_full_nn_mean_std = mse_full_nn_mean .+ sqrt.(mse_full_var)
val_mse_full_nn_mean_std = val_mse_full_nn_mean .+ sqrt.(val_mse_full_var)

p1 = Plots.plot(iters, mse_full_nn_mean; label = "mean")
Plots.plot!(iters, mse_full_nn_mean_std; label = "")
Plots.plot!(iters, mse_full_nn_mean; fillrange = mse_full_nn_mean_std, fillalpha = 0.35, c = 1, label = "")
Plots.xlabel!("Iteration")
Plots.ylabel!("Train MSE (full)")
Plots.plot!(; left_margin = 40 * Plots.PlotMeasures.px)

p2 = Plots.plot(iters, val_mse_full_nn_mean; label = "mean")
Plots.plot!(iters, val_mse_full_nn_mean_std; label = "")
Plots.plot!(iters, val_mse_full_nn_mean; fillrange = val_mse_full_nn_mean_std, fillalpha = 0.35, c = 1, label = "")
Plots.xlabel!("Iteration")
Plots.ylabel!("Validation MSE (full)")
Plots.plot!(; left_margin = 40 * Plots.PlotMeasures.px)

Plots.plot(p1, p2; layout = (1, 2))

folder = joinpath(@__DIR__, "output", string(config_basename, "_julia_parallel"))
mkpath(folder)
Plots.savefig(joinpath(folder, string(config_basename, "_mse.png")))
open(string(folder, config_basename, ".txt"), "w") do io
    write(io, outdir_path)
end

using Test
@testset "Julia Parallel Calibrate" begin
    batch_size = get_entry(config["reference"], "batch_size", nothing)
    if isnothing(batch_size)
        # Test convergence
        @test mse_full_mean[end] < mse_full_mean[1]
        # Test collapse
        if algorithm != "Unscented"
            @test mse_full_var[end] < mse_full_var[1]
        end
    else
        # Test convergence
        @test val_mse_full_mean[end] < val_mse_full_mean[1]
        # Test collapse
        if algorithm != "Unscented"
            @test val_mse_full_var[end] < val_mse_full_var[1]
        end
    end
end


text:
module DoorKeyUndirectedModule

import ..GridWorlds as GW
import Random
import ReinforcementLearningBase as RLBase

#####
##### game logic
#####

const NUM_OBJECTS = 5
const AGENT = 1
const WALL = 2
const GOAL = 3
const DOOR = 4
const KEY = 5
const NUM_ACTIONS = 5

mutable struct DoorKeyUndirected{R, RNG} <: GW.AbstractGridWorld
    tile_map::BitArray{3}
    agent_position::CartesianIndex{2}
    reward::R
    rng::RNG
    done::Bool
    terminal_reward::R
    goal_position::CartesianIndex{2}
    door_position::CartesianIndex{2}
    key_position::CartesianIndex{2}
    partition_dimension::Int
    has_key::Bool
end

function DoorKeyUndirected(; R = Float32, height = 8, width = 8, rng = Random.GLOBAL_RNG)
    tile_map = falses(NUM_OBJECTS, height, width)

    tile_map[WALL, 1, :] .= true
    tile_map[WALL, height, :] .= true
    tile_map[WALL, :, 1] .= true
    tile_map[WALL, :, width] .= true

    partition_dimension = 2

    door_position = CartesianIndex(2, 3)
    tile_map[WALL, 2 : height - 1, door_position[2]] .= true
    tile_map[WALL, door_position] = false
    tile_map[DOOR, door_position] = true

    agent_position = CartesianIndex(2, 2)
    tile_map[AGENT, agent_position] = true

    key_position = CartesianIndex(3, 2)
    tile_map[KEY, key_position] = true

    goal_position = CartesianIndex(height - 1, width - 1)
    tile_map[GOAL, goal_position] = true

    reward = zero(R)
    done = false
    terminal_reward = one(R)
    has_key = false

    env = DoorKeyUndirected(tile_map, agent_position, reward, rng, done, terminal_reward, goal_position, door_position, key_position, partition_dimension, has_key)

    GW.reset!(env)

    return env
end

function GW.reset!(env::DoorKeyUndirected)
    tile_map = env.tile_map
    rng = env.rng

    _, height, width = size(tile_map)

    tile_map[AGENT, env.agent_position] = false
    tile_map[GOAL, env.goal_position] = false
    tile_map[KEY, env.key_position] = false
    tile_map[DOOR, env.door_position] = false
    if env.partition_dimension == 1
        tile_map[WALL, env.door_position[1], 2 : width - 1] .= false
    else
        tile_map[WALL, 2 : height - 1, env.door_position[2]] .= false
    end

    new_partition_dimension = rand(rng, 1:2)
    env.partition_dimension = new_partition_dimension

    if new_partition_dimension == 1
        i_new_door_position = rand(rng, 3 : height - 2)
        j_new_door_position = rand(rng, 2 : width - 1)
        tile_map[WALL, i_new_door_position, :] .= true
        regions = (CartesianIndices((2 : i_new_door_position - 1, 2 : width - 1)), CartesianIndices((i_new_door_position + 1 : height - 1, 2 : width - 1)))
    else
        i_new_door_position = rand(rng, 2 : height - 1)
        j_new_door_position = rand(rng, 3 : width - 2)
        tile_map[WALL, :, j_new_door_position] .= true
        regions = (CartesianIndices((2 : height - 1, 2 : j_new_door_position - 1)), CartesianIndices((2 : height - 1, j_new_door_position + 1 : width - 1)))
    end

    new_door_position = CartesianIndex(i_new_door_position, j_new_door_position)
    env.door_position = new_door_position
    tile_map[WALL, i_new_door_position, j_new_door_position] = false
    tile_map[DOOR, i_new_door_position, j_new_door_position] = true

    new_agent_region = rand(rng, 1:2)
    if new_agent_region == 1
        new_agent_position, new_key_position = GW.sample_two_positions_without_replacement(rng, regions[1])
        new_goal_position = rand(rng, regions[2])
    else
        new_agent_position, new_key_position = GW.sample_two_positions_without_replacement(rng, regions[2])
        new_goal_position = rand(rng, regions[1])
    end

    env.agent_position = new_agent_position
    tile_map[AGENT, new_agent_position] = true

    env.key_position = new_key_position
    tile_map[KEY, new_key_position] = true

    env.goal_position = new_goal_position
    tile_map[GOAL, new_goal_position] = true

    env.reward = zero(env.reward)
    env.done = false
    env.has_key = false

    return nothing
end

function GW.act!(env::DoorKeyUndirected, action)
    @assert action in Base.OneTo(NUM_ACTIONS) "Invalid action $(action). Action must be in Base.OneTo($(NUM_ACTIONS))"

    tile_map = env.tile_map
    agent_position = env.agent_position

    if action in Base.OneTo(4)
        if action == 1
            new_agent_position = GW.move_up(agent_position)
        elseif action == 2
            new_agent_position = GW.move_down(agent_position)
        elseif action == 3
            new_agent_position = GW.move_left(agent_position)
        else
            new_agent_position = GW.move_right(agent_position)
        end

        if !(tile_map[DOOR, new_agent_position] || tile_map[WALL, new_agent_position]) || (tile_map[DOOR, new_agent_position] && env.has_key)
            tile_map[AGENT, agent_position] = false
            env.agent_position = new_agent_position
            tile_map[AGENT, new_agent_position] = true
        end
    else
        if tile_map[KEY, agent_position]
            env.has_key = true
            tile_map[KEY, agent_position] = false
        end
    end

    if tile_map[GOAL, env.agent_position]
        env.reward = env.terminal_reward
        env.done = true
    else
        env.reward = zero(env.reward)
        env.done = false
    end

    return nothing
end

#####
##### miscellaneous
#####

GW.get_height(env::DoorKeyUndirected) = size(env.tile_map, 2)
GW.get_width(env::DoorKeyUndirected) = size(env.tile_map, 3)

GW.get_action_names(env::DoorKeyUndirected) = (:MOVE_UP, :MOVE_DOWN, :MOVE_LEFT, :MOVE_RIGHT, :PICK_UP)
GW.get_object_names(env::DoorKeyUndirected) = (:AGENT, :WALL, :GOAL, :DOOR, :KEY)

function GW.get_pretty_tile_map(env::DoorKeyUndirected, position::CartesianIndex{2})
    characters = ('☻', '█', '♥', '▒', '⚷', '⋅')

    object = findfirst(@view env.tile_map[:, position])
    if isnothing(object)
        return characters[end]
    else
        return characters[object]
    end
end

function GW.get_pretty_sub_tile_map(env::DoorKeyUndirected, window_size, position::CartesianIndex{2})
    tile_map = env.tile_map
    agent_position = env.agent_position

    characters = ('☻', '█', '♥', '▒', '⚷', '⋅')

    sub_tile_map = GW.get_sub_tile_map(tile_map, agent_position, window_size)

    object = findfirst(@view sub_tile_map[:, position])
    if isnothing(object)
        return characters[end]
    else
        return characters[object]
    end
end

function Base.show(io::IO, ::MIME"text/plain", env::DoorKeyUndirected)
    str = "tile_map:\n"
    str = str * GW.get_pretty_tile_map(env)
    str = str * "\nsub_tile_map:\n"
    str = str * GW.get_pretty_sub_tile_map(env, GW.get_window_size(env))
    str = str * "\nreward: $(env.reward)"
    str = str * "\ndone: $(env.done)"
    str = str * "\nhas_key: $(env.has_key)"
    str = str * "\naction_names: $(GW.get_action_names(env))"
    str = str * "\nobject_names: $(GW.get_object_names(env))"
    print(io, str)
    return nothing
end

GW.get_action_keys(env::DoorKeyUndirected) = ('w', 's', 'a', 'd', 'p')

#####
##### RLBase API
#####

RLBase.StateStyle(env::GW.RLBaseEnv{E}) where {E <: DoorKeyUndirected} = RLBase.InternalState{Any}()
RLBase.state_space(env::GW.RLBaseEnv{E}, ::RLBase.InternalState) where {E <: DoorKeyUndirected} = nothing
RLBase.state(env::GW.RLBaseEnv{E}, ::RLBase.InternalState) where {E <: DoorKeyUndirected} = env.env.tile_map

RLBase.reset!(env::GW.RLBaseEnv{E}) where {E <: DoorKeyUndirected} = GW.reset!(env.env)

RLBase.action_space(env::GW.RLBaseEnv{E}) where {E <: DoorKeyUndirected} = Base.OneTo(NUM_ACTIONS)
(env::GW.RLBaseEnv{E})(action) where {E <: DoorKeyUndirected} = GW.act!(env.env, action)

RLBase.reward(env::GW.RLBaseEnv{E}) where {E <: DoorKeyUndirected} = env.env.reward
RLBase.is_terminated(env::GW.RLBaseEnv{E}) where {E <: DoorKeyUndirected} = env.env.done

end # module

ext: jl


text:
module testGenerate

# package code goes here
function test()
	println("Testing to see if a generated package can work.")
	file = joinpath("$(Pkg.dir("testGenerate"))", "deps", "open.so")
	handle = Libdl.dlopen(file)

	funcPtr = Libdl.dlsym(handle, :hello)
	ccall(funcPtr, Void, ())

end

end # module

using testGenerate
testGenerate.test()

max_stars_repo_name: camruggles/testGenerate

max_issues_repo_name: camruggles/testGenerate

max_forks_repo_name: camruggles/testGenerate


text:
global neurons = []

type Neuron
	connections::Dict{UInt64,Float16}
	potential::Float16
	error::Float16

	function Neuron(arg1,arg2,arg3)
		self = new(arg1,arg2,arg3)
		push!(neurons, self)
	end

end

function fully_connect(self)
	for neuron in neurons
		if object_id(neuron) != object_id(self)
			self.connections[object_id(neuron)] = rand(1:100)/100
			#push!(self.connections, rand(1:100)/100)
		end
	end
end

function partially_connect(self)
	if isempty(self.connections)
		neuron_count = length(neurons)
		#for neuron in neurons
		elected = rand(neurons,100)
		for neuron in elected
			if object_id(neuron) != object_id(self)
				#if rand(1:neuron_count/100) == 1
				self.connections[object_id(neuron)] = rand(1:100)/100
					#push!(self.connections, rand(1:100)/100)
				#end
			end
		end
		println("Neuron ID: ",object_id(self))
		println("    Potential: ",self.potential)
		println("    Error: ",self.error)
		println("    Connections: ",length(self.connections))
	end
end

function Build()
	for i = 1:1000000
		Neuron(Dict(),0.0,0.0)
	end
	println(length(neurons), " neurons created.")
	n = 0
	@parallel for neuron in neurons
		n += 1
		partially_connect(neuron)
		println("Counter: ",n)
	end
end

Build()

max_stars_repo_name: Maxprofs/Cerebrum


text:
# vinicius- creation of the system directory
# Legend of the systems: 020 = 0.5, 0.35 = 1.0, 0.50 = 1.5, 0.65 = 2.0, 0.80 = 2.5, 0.95 = 3.0

  include("input_gen.jl")
" Example :: EMIMDCABF4 - The system will contain the ions EMIM (two mols), DCA and BF4 (each one with one mole)"
  function folder(name::String,MM1,MM2)
    concs  = ["0.50", "1.00", "1.50", "2.00", "2.50", "3.00"]
    replic = ["00", "01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11","12", "13", "14", "15", "16", "17", "18", "19", "20"]
     gmds = ["gmd_an1.inp","gmd_an2.inp","gmd_sol.inp", "gmd_ca1.inp"]

    println("entrou")

    main_dir = pwd()
    for c in concs
      for r in replic

        work_dir = joinpath(main_dir,name,c,r)
        cil = parse(Float64,c)
        cd(work_dir)
        conc(MM1,"$(name[1:4])","$(name[5:7])",MM2,"EMIM","$(name[8:10])",cil)

        for file in gmds
          w = open(file,"r+")
          new = open("gmdx_$(file[5:7]).inp","w")

          # change names
          for line in eachline(w)
            if occursin("AN1",line)
              if "$(name[5:7])" == "DCA"
                println(new,replace(line, "AN1" => "NC", count=1))
              else
                println(new,replace(line, "AN1" => "$(name[8:10])", count=1))
              end
            elseif occursin("dist",line)
              println(new,replace(line, "dist" => "20", count=1))
            elseif occursin("AN2",line)
              if "$(name[8:10])" == "DCA"
                println(new,replace(line, "AN2" => "NC", count=1))
              else
                println(new,replace(line, "AN2" => "$(name[8:10])", count=1))
              end
            elseif  occursin("CA1",line)
              println(new,replace(line, "CA1" => "$(name[1:3])", count=1))
            else
              println(new,line)
            end

          end
          close(new)
          close(w)

          rm(file)

        end

        cd(main_dir)
      end
    end
  end

  folder("EMIMDCABF4",177.21, 197.97)


text:
using Taxonomy
using Documenter

DocMeta.setdocmeta!(Taxonomy, :DocTestSetup, :(using Taxonomy); recursive=true)

function readme2index()
    readme_path = "README.md"
    index_path = "docs/src/index.md"
    f = open(readme_path,"r")
    g = open(index_path, "w")
    try
        readme = read(f,String)
        index = replace(readme,"![](docs/src/img" => "![](img")
        write(g,readme)
    finally
        close(f)
        close(g)
    end
end

readme2index()

makedocs(;
    modules=[Taxonomy],
    authors="banhbio <ban@kuicr.kyoto-u.ac.jp>",
    repo="https://github.com/banhbio/Taxonomy.jl/blob/{commit}{path}#{line}",
    sitename="Taxonomy.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://banhbio.github.io/Taxonomy.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/banhbio/Taxonomy.jl.git",
    devbranch="main",
)


text:
include("../src/components/grosseconomy_component.jl")

@testset "grosseconomy" begin

    Precision = 1.0e-11
    T = length(MimiDICE2010.model_years)
    f = readxlsx(joinpath(@__DIR__, "..", "data", "DICE2010_082710d.xlsx"))

    m = Model()

    set_dimension!(m, :time, MimiDICE2010.model_years)

    add_comp!(m, grosseconomy, :grosseconomy)

# Set the parameters that would normally be internal connection from their Excel values
    set_param!(m, :grosseconomy, :I, read_params(f, "B101:BI101", T))

# Load the rest of the external parameters
    p = dice2010_excel_parameters(joinpath(@__DIR__, "..", "data", "DICE2010_082710d.xlsx"))
    set_param!(m, :grosseconomy, :al, p[:al])
    set_param!(m, :grosseconomy, :l, p[:l])
    set_param!(m, :grosseconomy, :gama, p[:gama])
    set_param!(m, :grosseconomy, :dk, p[:dk])
    set_param!(m, :grosseconomy, :k0, p[:k0])

# Run the one-component model
    run(m)

# Extract the generated variables
    K       = m[:grosseconomy, :K]
    YGROSS  = m[:grosseconomy, :YGROSS]

# Extract the true values
    True_K      = read_params(f, "B102:BI102", T)
    True_YGROSS = read_params(f, "B92:BI92", T)

# Test that the values are the same
    @test maximum(abs, K .- True_K) ≈ 0. atol = Precision
    @test maximum(abs, YGROSS .- True_YGROSS) ≈ 0. atol = Precision

end


text:
using .PDMats

Base.:\(A::PDMat, B::TrackedVecOrMat) = Tracker.track(\, A, B)
Base.:\(A::PDiagMat, B::TrackedVecOrMat) = Tracker.track(\, A, B)
Base.:\(A::ScalMat, B::TrackedVecOrMat) = Tracker.track(\, A, B)


text:
"""
    getkeldata(data::T) where T <: PKSubject
"""
function getkeldata(data::T) where T <: PKSubject
    data.keldata
end
"""
    getkeldata(data::T) where T <: PKSubject
"""
function getkeldata(ncar::T) where T <: NCAResult
    getkeldata(ncar.data)
end


text:
abstract ClassicalLattice{T} <: Lattice{T}

function buildClassicalLattice(latticeinfo::LatticeInfo, spinmodel::SpinModel)
	if isTwoDimension(latticeinfo) & isSquareLattice(latticeinfo)
		return Classical2dSquareLattice(spinmodel)
	elseif isTwoDimension(latticeinfo) & isFractalLattice(latticeinfo)
		return Classical2dFractalLattice(spinmodel, getLegextension(latticeinfo))
	elseif getDimension(latticeinfo,3) & isSquareLattice(latticeinfo)
		return Classical3dSquareLattice(spinmodel)
	else
		error("wrong latticeinfo")
	end
end


text:
function tick!(
	agents::AbstractArray,
	network::AbstractGraph
)
	acting_agent = StatsBase.sample(agents)
	interaction_partner = draw_interaction_partner(agents, network, acting_agent)
	if !acting_agent.stubborn
		similarity = compute_similarity(acting_agent, interaction_partner)
		if (rand() < similarity) && (similarity != 1)
			assimilate!(agents, acting_agent, interaction_partner)
		end
	end
	return (agents, network)
end

function draw_interaction_partner(
	agents::AbstractArray,
	network::AbstractGraph,
	agent::Agent
)
	agent_neighbors = LightGraphs.neighbors(network, agent.id)
	interaction_partner_id = rand(agent_neighbors)
	interaction_partner = agents[interaction_partner_id]
	return interaction_partner
end

function compute_similarity(
	acting_agent::Agent,
	interaction_partner::Agent
)
	shared_traits = sum(
		[
			i == j
			for (i, j) in zip(acting_agent.culture, interaction_partner.culture)
		]
	)
	similarity = shared_traits / length(acting_agent.culture)
	return similarity
end

function assimilate!(
	agents::AbstractArray,
	acting_agent::Agent,
	interaction_partner::Agent
)
	random_attr = rand(1:length(acting_agent.culture))
	if !(acting_agent.culture[random_attr] == interaction_partner.culture[random_attr])
		acting_agent.culture[random_attr] = interaction_partner.culture[random_attr]
	else
		assimilate!(agents, acting_agent, interaction_partner)
	end
	return acting_agent
end

max_stars_repo_name: JohannesNakayama/axelrodnet

max_issues_repo_name: JohannesNakayama/axelrodnet

max_forks_repo_name: JohannesNakayama/axelrodnet


text:
using Printf
using Plots
using ProgressMeter

function animate(input, traj, filename; size = [800,400], fps :: Int64 = 10, last=0, dpi=150)

  ENV["GKSwstype"]="nul"

  U = traj.U
  S = traj.S
  D = traj.D
  I = traj.I
  snenc = @sprintf("%5.2f",sum(traj.nenc)/traj.nsteps)

  if last < 1
    last = traj.nsteps
  end

  # for subplot 1
  lims=[-input.side/2,input.side/2]

  # for subplot 2
  time = traj.time
  xmin = 0.
  xmax = time[last]
  ymin = 0.
  ymax = 1.

  p = Progress(last,1)
  anim = @animate for i in 1:last

    x = traj.atoms[i].x[:,1];
    y = traj.atoms[i].x[:,2];
    c = Vector{String}(undef,traj.n)
    for j in 1:traj.n
      if traj.atoms[i].status[j] == 0
        c[j] = "blue"
      elseif traj.atoms[i].status[j] == 1
        c[j] = "red"
      elseif traj.atoms[i].status[j] == 2
        c[j] = "white"
      elseif traj.atoms[i].status[j] == 3
        c[j] = "green"
      end
    end

    plot(size=size,layout=(1,2),framestyle=:box,dpi=dpi)

    plot!(xlabel=@sprintf("Time: %4i",time[i]),subplot=1)
    scatter!(x,y,label="",color=c,xlim=lims,ylim=lims,subplot=1,markersize=2,xticks=:none,yticks=:none)

    #plot!(title=@sprintf("Tempo: %4i",i),subplot=2)
    plot!(time[1:i],S[1:i],subplot=2,linewidth=2,label="",color="red")
    plot!(time[1:i],U[1:i],subplot=2,linewidth=2,label="",color="blue")
    plot!(time[1:i],D[1:i],subplot=2,linewidth=2,label="",color="black")
    plot!(time[1:i],I[1:i],subplot=2,linewidth=2,label="",color="darkgreen")
    plot!(xlim=[0,xmax],ylim=[ymin,ymax],subplot=2)
    plot!(xlabel="Time",ylabel="Fraction of population",subplot=2)

    fontsize=8
    plot!(rectangle( 0, 0.37*xmax, 0.815*ymax, 1.02*ymax), opacity=0.9,label="",color=:white,subplot=2)
    x = 0.05*(xmax-xmin)+xmin
    d = 0.05*ymax ; y = [ ymax + d ]
    annotate!(x,yd!(y,d),text("Susceptible: $(@sprintf("%4.2f%%",100*U[i]))",:left,fontsize,:serif,:blue),subplot=2)
    annotate!(x,yd!(y,d),text("Sick: $(@sprintf("%4.2f%%",100*S[i]))",:left,fontsize,:serif,:red),subplot=2)
    annotate!(x,yd!(y,d),text("Dead: $(@sprintf("%4.2f%%",100*D[i]))",:left,fontsize,:serif,:black),subplot=2)
    annotate!(x,yd!(y,d),text("Immune: $(@sprintf("%4.2f%%",100*I[i]))",:left,fontsize,:serif,:darkgreen),subplot=2)

    plot!(rectangle( 0.53*xmax, 1.00*xmax, 0.810*ymax, 1.02*ymax), opacity=0.9,label="",color=:white,subplot=2)
    x = 0.98*(xmax-xmin)+xmin
    d = 0.05*ymax ; y = [ ymax + d ]
    annotate!(x,yd!(y,d),text("\"Temperature\": $(input.kavg_target)",:right,fontsize,:serif,:black),subplot=2)
    annotate!(x,yd!(y,d),text("New encounters/step: $snenc",:right,fontsize,:serif,:black),subplot=2)
    annotate!(x,yd!(y,d),text("Cross-section: $(input.xsec)",:right,fontsize,:serif,:black),subplot=2)
    annotate!(x,yd!(y,d),text("P(contamination): $(input.pcont)",:right,fontsize,:serif,:black),subplot=2)

    next!(p)

  end

  gif(anim, filename, fps = fps)

end


text:
module Training

export trainepochs!, StopException

using Flux
# import Flux.Tracker: Params, gradient, data, update!
import MLDataUtils: getobs

function _update_params_full!(pscb::AbstractArray, xs; kwargs...)
  #  @show pscb, kwargs
    for c in pscb
        apply!(c, xs; kwargs...)
    end
end

_update_params_full!(pscb, xs; kwargs...) = _update_params_full!([pscb], xs; kwargs...)

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

#from Flux/optimise/train.jl
# function update!(opt, x, x̄)
    # update!(x, -Flux.Optimise.apply!(opt, x, data(x̄)))
# end

#from Flux/optimise/train.jl
# function update!(opt, xs::Params, gs)
    # for x in xs
        # update!(opt, x, gs[x])
    # end
# end

# macro interrupts(ex)
    # :(try $(esc(ex))
    # catch e
        # e isa InterruptException || rethrow()
        # throw(e)
    # end)
# end

#from Flux/optimise/train.jl
# struct StopException <: Exception end

#from Flux/optimise/train.jl
# function stop()
    # throw(StopException())
# end

function trainepochs!(loss, ps, data, nbatches, opt;
    pscb = ParamsIdentity(),
    epochcb = nothing,
    batchcb = nothing)

    epoch = 0
    isnothing(epochcb) || epochcb(epoch = epoch, epoch_start_time = time(), ps = ps)
    epoch_start_time = time()
    for (b, s) in enumerate(data)
        batch = (b-1) % nbatches + 1

        # @show batch
        # s = getobs(s)
        # @show summary.(s)
        # @show summary(loss(s...))
        # @show summary.(ps)
        try
            gs = gradient(ps) do
                loss(s...)
            end
            isnothing(batchcb) || batchcb(ps, epoch = epoch, batch = batch, batch_size = size(samples[1])[end])
            update!(opt, ps, gs)
        catch ex
            if ex isa Flux.Optimise.StopException
                break
            else
                rethrow(ex)
            end
        end
        if  batch == nbatches
            epoch += 1
            try
                isnothing(epochcb) || epochcb(epoch = epoch, epoch_start_time = epoch_start_time, ps = ps)
            catch ex
                if ex isa Flux.Optimise.StopException
                    break
                else
                    rethrow(ex)
                end
            end
            epoch_start_time = time()
        end
    end
end

mutable struct ParamsIdentity
end

apply!(o::ParamsIdentity, xs; kwargs...) = xs

mutable struct ParamsReport
    batches
end

function apply!(o::ParamsReport, xs; epoch, batch, kwargs...)
    if batch ∉ o.batches
        return xs
    end
    param_min = flatten((abs.(x.data) for x in xs)) |> minimum
    param_max = flatten((abs.(x.data) for x in xs)) |> maximum
    param_mean = flatten((abs.(x.data) for x in xs)) |> mean

    grad_min = flatten((x.grad for x in xs)) |> minimum
    grad_max = flatten((x.grad for x in xs)) |> maximum
    grad_mean = flatten((x.grad for x in xs)) |> mean
    grad_L2 = collect(flatten((x.grad for x in xs))) |> norm
    println("$batch: |θ| [$param_min;$param_max] mean: $param_mean; Δ [$grad_min;$grad_max] mean: $grad_mean L2: $grad_L2")
    xs
end

mutable struct ClipGradient
    lb::Float32
    ub::Float32
end

ClipGradient() = ClipGradient(-1f0, 1f0)

function apply!(o::ClipGradient, xs; kwargs...)
    for x in xs
        x.grad .= min.(max.(x.grad, o.lb), o.ub)
    end
    xs
end

mutable struct ClipGradientNorm
    max_norm::Float32
    norm_type::Float32
end

ClipGradientNorm() = ClipGradientNorm(1.0f0, 2f0)

  # see https://pytorch.org/docs/stable/_modules/torch/nn/utils/clip_grad.html
function apply!(o::ClipGradientNorm, xs; kwargs...)
    total_norm = 0
    for x in xs
        total_norm += norm(x.grad)^o.norm_type
    end
    total_norm = total_norm^(1. / o.norm_type)
    clip_coef = o.max_norm / (total_norm + 1e-6)

    if clip_coef < 1
  #    @show "PRE", collect(flatten((x.grad for x in xs))) |> norm, clip_coef
        for x in xs
            x.grad .*= clip_coef
        end
  #    @show "AFT", collect(flatten((x.grad for x in xs))) |> norm, clip_coef
    end
    xs
end

end

max_stars_repo_name: dhonza/FluxGoodies

max_issues_repo_name: dhonza/FluxGoodies

max_forks_repo_name: dhonza/FluxGoodies


text:
function add!(d::Dag, parents)  # TODO rename sum!
    # TODO parents is similar to `combine!` ... need a type info for doc
    cn = combine!(d, parents)

    # combine! return a BVal{Tuple{...}}
    # the retuned node type after add! is promote_type of all types in tuple
    CT = promote_type(eltype(eltype(cn)).parameters...)  # TODO method of BVal
    @assert promote_type(map(p -> eltype(p[2]), parents)...) == CT

    return node!(d, CT, [(cn, v -> Nullable(sum(v.val)))])
end

# TODO doc; similar to combining version, but all nodes always trigger and
#      always use Latest value of previous nodes
function add!(d::Dag, nodes...)
    return apply!(d, sum, 0, nodes...)
end

function first!(d::Dag, node::Node{T}) where {T<:Tuple}
    nth!(d, node, 1)
end

function last!(d::Dag, node::Node{T}) where {T<:Tuple}
    # Cannot call nth! directly as below
    # because `end` is not recognize as last elt index
    # nth!(d, node, end)
    nth!(d, node, length(eltype(node).types))
end

function nth!(d::Dag, node::Node{T}, n) where {T<:Tuple}
    EE = eltype(node).types[n]
    node!(d, EE, [(node, v -> Nullable(v[n]))])
end

# create a moving window node, i.e. transform input in vector of last N values
# TODO maybe have a diff output: instead of vector could also output what's
#      added and what's removed (for online stream calculation) ... maybe a
#      diff win! interface can do that
function win!(d::Dag, node::Node, n)
    EE = Vector{eltype(node)}
    # TODO Datastructures.jl CircularDeque
    vec = EE()

    function win_push!(v)
        while length(vec) >= n
            shift!(vec)
        end
        push!(vec, v)
    end

    # node!(d, EE, [(node, v -> Nullable(win_push!(v)))])
    tr!(d, node, win_push!)
end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the Lattice-Boltzmann equations for the D3Q27 scheme

L = 1 # reference length
equations = LatticeBoltzmannEquations3D(Ma=0.1, Re=1600; L=L)

initial_condition = initial_condition_taylor_green_vortex

solver = DGSEM(polydeg=3, surface_flux=flux_godunov)

coordinates_min = (-pi*L, -pi*L, -pi*L)
coordinates_max = ( pi*L,  pi*L,  pi*L)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=5,
                n_cells_max=300_000,)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 20*equations.L/equations.u0) # Final time is `20` in non-dimensional time
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 20
analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
                                     save_analysis=true,
                                     extra_analysis_integrals=(Trixi.energy_kinetic_nondimensional,))

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2macroscopic)

stepsize_callback = StepsizeCallback(cfl=0.3)

collision_callback = LBMCollisionCallback()

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        stepsize_callback,
                        collision_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks, maxiters=1e5,
            save_start=false, alias_u0=true);
summary_callback() # print the timer summary


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

# Harris Corner Detection
# Original port by Todd A. Anderson from PolyMage version.

using DocOpt
using Images

# Compute and return the application of a stencil to an input 2D array
# shape specifies the offsets from output pixel of the desired points of source
# kf takes these values in same order and computes value of output pixel
# boundary is just zeros (where any point doesn't exist)
function run(src, shape, kf)
  (s1, s2) = size(src)
  min1 = 1
  max1 = s1
  min2 = 1
  max2 = s1
  npts = length(shape)
  for p in 1:npts
    (o1, o2) = shape[p]
    min1 = max(min1, 1-o1)
    max1 = min(max1, s1-o1)
    min2 = max(min2, 1-o2)
    max2 = min(max2, s2-o2)
  end
  dst = Array(Float32, s1, s2)
  args = Array(Float32, npts)
  for i1 in 1:s1
    for i2 in 1:s2
      if i1>=min1 && i1<=max1 && i2>=min2 && i2<=max2
        for p in 1:npts
          (o1, o2) = shape[p]
          args[p] = src[i1+o1, i2+o2]
        end
        dst[i1, i2] = kf(args)
      else
        dst[i1, i2] = 0.0f0
      end
    end
  end
  return dst
end

function harrisCornerDetect(Iin)
  (w, h) = size(Iin)

  Ix_args(args) = ((args[1] * -1.0f0) + (args[2] * -2.0f0) + (args[3] * -1.0f0) + args[4] + (args[5] * 2.0f0) + args[6]) / 12.0f0
  Ix = run(Iin, [(-1,-1), (-1,0), (-1,1), (1,-1), (1,0), (1,1)], Ix_args)
  Iy_args(args) = ((args[1] * -1.0f0) + (args[2] * -2.0f0) + (args[3] * -1.0f0) + args[4] + (args[5] * 2.0f0) + args[6]) / 12.0f0
  Iy = run(Iin, [(-1,-1), (0,-1), (1,-1), (-1,1), (0,1), (1,1)], Iy_args)

  Ixx = Ix .* Ix
  Iyy = Iy .* Iy
  Ixy = Ix .* Iy

  S_args(args) = sum(args)
  Sxx = run(Ixx, [(-1,-1),(-1,0),(-1,1),(0,-1),(0,0),(0,1),(1,-1),(1,0),(1,1)], S_args)
  Syy = run(Iyy, [(-1,-1),(-1,0),(-1,1),(0,-1),(0,0),(0,1),(1,-1),(1,0),(1,1)], S_args)
  Sxy = run(Ixy, [(-1,-1),(-1,0),(-1,1),(0,-1),(0,0),(0,1),(1,-1),(1,0),(1,1)], S_args)

  det    = (Sxx .* Syy) .- (Sxy .* Sxy)
  trace  = Sxx .+ Syy
  harris = det .- (0.04f0 .* trace .* trace)

  return harris
end

function main()
    doc = """harris.jl

Harris corner detection.

Usage:
  harris.jl -h | --help
  harris.jl [--img-file=<img-file>]

Options:
  -h --help              Show this screen.
  --img-file=<img-file>  Specify a path to an input (grayscale) image file [default: ../../example.jpg].
"""
    arguments = docopt(doc)

    dir = dirname(@__FILE__)
    img_file = joinpath(dir, arguments["--img-file"])

    (fname, ext) = splitext(img_file)
    out_file = string(fname, "-corners", ".jpg")

    println("input file = ", img_file)
    println("output file = ", out_file)

    function harris(input_fname, output_fname)
        local img = convert(Matrix{Float32}, load(input_fname))
        tic()
        res = harrisCornerDetect(Matrix{Float32}())
        println("SELFPRIMED ", toq())

        tic()
        res = harrisCornerDetect(img)
        selftimed = toq()
        res = map((x,y) -> x > 0.001f0 ? 1.0f0 : (y / 2.0f0), res, img)
        save(output_fname, ufixed8sc(convert(Image, res)))
        println("checksum: ", sum(res))
        println("SELFTIMED ", selftimed)
    end

    harris(img_file, out_file)

end

main()


text:
"""
    SRational

# Examples
```jldoctest
julia> SInt(3) // SInt(5)
3//5

julia> (SInt(3) // SInt(5)) // (SInt(2) // SInt(1))
3//10
```
"""
struct SRational{N<:SInteger,D<:SInteger} <: Real end

SRational(num::SIntegerZeroType, den::SIntegerZeroType) =
    throw("invalid rational: zero($(eltype(num))//zero($(eltype(den))")

for (S,B) in SI2BI
    @eval begin
        function SRational(num::$S, den::$S)
            num2, den2 = (sign(den) < 0) ? Base.divgcd(-num, -den) : Base.divgcd(num, den)
            SRational{typeof(num),typeof(den)}()
        end
    end
end

SRational(n::SInteger, d::SInteger) = SRational(promote(n,d)...)
SRational(n::SInteger) = SRational(n,one(n))

Base.numerator(::SRational{N,D}) where {N<:SInteger,D<:SInteger} = N()
Base.denominator(::SRational{N,D}) where {N<:SInteger,D<:SInteger} = D()

Base.decompose(x::SRational) = numerator(x), SZero, denominator(x)

Base.:(//)(n::SInteger,  d::SInteger) = SRational(n,d)

function Base.:(//)(x::SRational, y::SInteger)
    xn, yn = Base.divgcd(numerator(x), y)
    xn//(denominator(x) * yn)
end

function Base.:(//)(x::SInteger,  y::SRational)
    xn, yn = Base.divgcd(x, numerator(y))
    (xn * denominator(y))//yn
end
function Base.:(//)(x::SRational, y::SRational)
    xn, yn = Base.divgcd(numerator(x), numerator(y))
    xd, yd = Base.divgcd(denominator(x), denominator(y))
    (xn * yd)//(xd * yn)
end

#//(x::SComplex,  y::SReal) = complex(real(x)//y,imag(x)//y)
#//(x::SNumber, y::SComplex) = x*conj(y)//abs2(y)

#//(X::AbstractArray, y::SNumber) = X .// y

#=
Base.one(::SRational{<:Any,T}) where {V,T<:BaseInteger} = SRational{Tuple{T(1)::T,T(1)::T},T}()
Base.zero(::SRational{<:Any,T}) where {V,T<:BaseInteger} = SRational{Tuple{T(0)::T,T(1)::T},T}()
Base.one(::Type{<:SRational{<:Any,T}}) where {V,T<:BaseInteger} = SRational{Tuple{T(1)::T,T(1)::T},T}()
Base.zero(::Type{<:SRational{<:Any,T}}) where {V,T<:BaseInteger} = SRational{Tuple{T(0)::T,T(1)::T},T}()
=#

# TODO
promote_rule(::Type{SRational{<:Any,T}}, ::Type{S}) where {T<:BaseInteger,S<:SInteger} = SRational{Any,promote_type(T,S)}
promote_rule(::Type{SRational{<:Any,T}}, ::Type{SRational{S}}) where {T<:BaseInteger,S<:BaseInteger} = SRational{Any,promote_type(T,S)}
promote_rule(::Type{SRational{<:Any,T}}, ::Type{S}) where {T<:BaseInteger,S<:SFloat} = promote_type(T, S)


text:
# CR3BP Parameters
struct CR3BPParams

    # Primary and Secondary Body Mass
    m1::Float64
    m2::Float64

    # Radius of Primary and Secondary Bodies
    R1::Float64
    R2::Float64

    # Total Mass
    mtot::Float64

    # Gravitational Parameters
    μ::Float64

    # Length Unit
    LU::Float64

    # Time Unit
    TU::Float64

    # Speed Unit
    VU::Float64

    # Mass Unit
    MU::Float64

    # Gravitational Constant
    G::Float64

    function CR3BPParams(m1,m2,R1,R2,MU,r12)
        G       = 6.673e-20
        mtot    = m1 + m2
        μ       = m2 / (m1 + m2)
        LU      = r12
        TU      = 1.0 / sqrt(G * (m1 + m2) / (LU^3))
        VU      = LU / TU
        MU      = MU
        new(m1,m2,R1,R2,mtot,μ,LU,TU,VU,MU,G)
    end
end

# CR3BP Eq. of Motion w/o control
function cr3bpEomNoControl(u::AbstractArray, p::CR3BPParams, t)
    # Precompute required
    xpμ = u[1] + p.μ
    r1  = sqrt(xpμ^2 + u[2]^2 + u[3]^2)
    r2  = sqrt((xpμ - 1)^2 + u[2]^2 + u[3]^2)

    # Compute g(r)
    invR13 = 1.0 / (r1^3)
    invR23 = 1.0 / (r2^3)
    gx = u[1] - (1.0 - p.μ)*xpμ*invR13 - p.μ*(xpμ - 1.0)*invR23
    gy = u[2] - (1.0 - p.μ)*u[2]*invR13 - p.μ*u[2]*invR23
    gz = -(1.0 - p.μ)*u[3]*invR13 - p.μ*u[3]*invR23

    # Compute h(v)
    hx = 2.0*u[5]
    hy = -2.0*u[4]

    return @SVector [u[4], u[5], u[6],
                     gx + hx,
                     gy + hy,
                     gz]
end
function cr3bpEomNoControl!(du::AbstractArray,u::AbstractArray,p::CR3BPParams,t)
    # Precompute required
    xpμ = u[1] + p.μ
    r1  = sqrt(xpμ^2 + u[2]^2 + u[3]^2)
    r2  = sqrt((xpμ - 1)^2 + u[2]^2 + u[3]^2)

    # Compute g(r)
    invR13 = 1.0 / (r1^3)
    invR23 = 1.0 / (r2^3)
    gx = u[1] - (1.0 - p.μ)*xpμ*invR13 - p.μ*(xpμ - 1.0)*invR23
    gy = u[2] - (1.0 - p.μ)*u[2]*invR13 - p.μ*u[2]*invR23
    gz = -(1.0 - p.μ)*u[3]*invR13 - p.μ*u[3]*invR23

    # Compute h(v)
    hx = 2.0*u[5]
    hy = -2.0*u[4]

    # Set du
    @. du[1:3] = view(u,4:6)
    du[4] = gx + hx
    du[5] = gy + hy
    du[6] = gz
end

# CR3BP Eq. of Motion w/ control
function cr3bpEomControl(u::AbstractArray,p::CR3BPParams,t,
                         a::AbstractArray)
    # Precompute required
    xpμ = u[1] + p.μ
    r1  = sqrt(xpμ^2 + u[2]^2 + u[3]^2)
    r2  = sqrt((xpμ - 1)^2 + u[2]^2 + u[3]^2)

    # Compute g(r)
    invR13 = inv(r1^3)
    invR23 = inv(r2^3)
    gx = u[1] - (1.0 - p.μ)*xpμ*invR13 - p.μ*(xpμ - 1.0)*invR23
    gy = u[2] - (1.0 - p.μ)*u[2]*invR13 - p.μ*u[2]*invR23
    gz = -(1.0 - p.μ)*u[3]*invR13 - p.μ*u[3]*invR23

    # Compute h(v)
    hx = 2.0*u[5]
    hy = -2.0*u[4]

    # Compute and return dynamics
    dx  = u[4]
    dy  = u[5]
    dz  = u[6]
    ddx = gx + hx + a[1]
    ddy = gy + hy + a[2]
    ddz = gz + a[3]

    return @SVector [dx,dy,dz,ddx,ddy,ddz]
end

function cr3bpEomControl!(du::AbstractArray,u::AbstractArray,p::CR3BPParams,t,
                          a::AbstractArray)
    # Precompute required
    xpμ = u[1] + p.μ
    r1  = sqrt(xpμ^2 + u[2]^2 + u[3]^2)
    r2  = sqrt((xpμ - 1)^2 + u[2]^2 + u[3]^2)

    # Compute g(r)
    invR13 = 1.0 / (r1^3)
    invR23 = 1.0 / (r2^3)
    gx = u[1] - (1.0 - p.μ)*xpμ*invR13 - p.μ*(xpμ - 1.0)*invR23
    gy = u[2] - (1.0 - p.μ)*u[2]*invR13 - p.μ*u[2]*invR23
    gz = -(1.0 - p.μ)*u[3]*invR13 - p.μ*u[3]*invR23

    # Compute h(v)
    hx = 2.0*u[5]
    hy = -2.0*u[4]

    # Set du
    @inbounds for i in 1:3
        du[i] = u[i + 3]
    end
    du[4] = gx + hx + a[1]
    du[5] = gy + hy + a[2]
    du[6] = gz + a[3]
end


text:
type SymFloat
sym::Array{Sym}
val::Array{Float64}
end
SymFloat(name::String) = SymFloat([Sym(name)])
SymFloat(var::Sym) = SymFloat([var])
SymFloat(sym::Array{Sym}) = SymFloat(sym,zeros(Float64,size(sym)))
SymFloat{T<:Real}(num::T) = SymFloat(Sym(num))
# SymFloat(num::Int64) = SymFloat(Sym(num))
SymFloat{T<:Real}(M::Array{T}) = SymFloat(Sym(M))
getindex(A::SymFloat,inds...) = SymFloat(A.sym[inds...],A.val[inds...])

function setindex!(a::SymFloat,arg::SymFloat,d1,d2)
  if isa(d1,Int64) && isa(d2,Int64) && isa(arg.sym,Array)
    a.sym[d1,d2] = arg.sym[1]
    a.val[d1,d2] = arg.val[1]
  else
    a.sym[d1,d2] = arg.sym
    a.val[d1,d2] = arg.val
  end
end
function setindex!(a::SymFloat,arg::SymFloat,d1)
    if isa(d1,Int64) && isa(arg.sym,Array)
      a.sym[d1] = arg.sym[1]
      a.val[d1] = arg.val[1]
    else
      a.sym[d1] = arg.sym
      a.val[d1] = arg.val
 end
end

convert(::Type{SymFloat},x::String) = SymFloat(x)
convert(::Type{SymFloat},x::Sym) = SymFloat(x)
convert(::Type{SymFloat},x::Array{Sym}) = SymFloat(x)

zeros(::Type{SymFloat},d...) = SymFloat(zeros(Sym,d...))

size(x::SymFloat) = size(x.val)
length(x::SymFloat) = length(x.val)

convert{T<:Real}(::Type{Array{T,1}}, x::T) = [x]
convert(::Type{Array{Sym,1}},x::Sym) = [x]

function ==(x1::SymFloat,x2::SymFloat)
  if x1.sym == x2.sym && x1.val == x2.val
    return true
  end
  return false
end

# function display(d::Display,M::MIME"text/plain",s::SymFloat)
#   display(d,M,s.sym)
#   display(d,M,"=")
#   display(d,M,s.val)
# end
# function show(io::IO, s::SymFloat)
#   print(io,s.sym)
# end


text:
using Statistics

logistic(x) = 1 ./ (1 .+ exp.(-x))

function logistic(t::Tensor)
    values = logistic(t.values)
    if t.requires_grad
        requires_grad = true
        dependencies = [
            Dependency(t, x -> x .* logistic(t.values) .* (1 .- logistic(t.values)))
        ]
    else
        requires_grad = false
        dependencies = nothing
    end
    Tensor(
        values;
        dependencies=dependencies,
        device=t.device
    )
end

function logloss(y::Tensor, ŷ::Tensor)
    values = [mean(-y.values .* log.(ŷ.values) - (1 .- y.values) .* log.(1 .- ŷ.values))]
    if y.requires_grad || ŷ.requires_grad
        requires_grad = true
        dependencies = [
            Dependency(y, x -> x .* (log.(1 .- ŷ.values) - log.(ŷ.values)) / length(y)),
            Dependency(ŷ, x -> x .* (ŷ.values - y.values) ./ (ŷ.values .* (1 .- ŷ.values) * length(y)))
        ]
    else
        requires_grad = false
        dependencies = nothing
    end
    Tensor(
        values;
        dependencies=dependencies,
        device=y.device
    )
end

function mse(y::Tensor, ŷ::Tensor)
    values = [mean((y.values - ŷ.values).^2 / 2)]
    if y.requires_grad || ŷ.requires_grad
        requires_grad = true
        dependencies = [
            Dependency(y, x -> x .* (y.values - ŷ.values) / length(y)),
            Dependency(ŷ, x -> x .* (ŷ.values - y.values) / length(y))
        ]
    else
        requires_grad = false
        dependencies = nothing
    end
    Tensor(
        values;
        dependencies=dependencies,
        device=y.device
    )
end

max_stars_repo_name: taraspiotr/DeepJulia

max_issues_repo_name: taraspiotr/DeepJulia

max_forks_repo_name: taraspiotr/DeepJulia


text:
@testset "830.positions-of-large-groups.jl" begin
    @test large_group_positions("abcdddeeeeaabbbcd") == [[3,5],[6,9],[12,14]]
    @test large_group_positions("abbxxxxzzy") == [[3,6]]
    @test large_group_positions("abc") == []
end


text:
using DataFrames, TableOperations, Tables, Random, Logging
using TreeParzen
using Recommenders:
    Movielens100k,
    load_dataset,
    ratio_split,
    SLIM,
    evaluate_u2i,
    MeanPrecision,
    MeanRecall,
    MeanNDCG

function main()
    ml100k = Movielens100k()
    download(ml100k)
    rating, _, _ = load_dataset(ml100k)

    rating = rating |> TableOperations.transform(Dict(:rating => x -> 1.0))
    for (i, row) in enumerate(Tables.rows(rating))
        println(row)
        if i >= 4
            break
        end
    end

    Random.seed!(1234)
    train_valid_table, test_table = ratio_split(rating, 0.8)
    train_table, valid_table = ratio_split(train_valid_table, 0.8)

    prec10 = MeanPrecision(10)
    recall10 = MeanRecall(10)
    ndcg10 = MeanNDCG(10)
    metrics = [prec10, recall10, ndcg10]

    space = Dict(
        :λminratio => HP.LogUniform(:λminratio, log(1e-6), log(0.5)),
        :l1_ratio => HP.LogUniform(:l1_ratio, log(1e-6), log(0.9999)),
        :k => HP.QuantUniform(:k, 1.0, 1000.0, 1.0),
    )

    function invert_output(params)
        @info params
        λminratio = params[:λminratio]
        l1_ratio = params[:l1_ratio]
        k = convert(Int, params[:k])

        model = SLIM(l1_ratio, λminratio, k)

        result = evaluate_u2i(
            model,
            train_table,
            valid_table,
            metrics,
            10,
            col_item = :movieid,
            drop_history = true,
        )
        @info result
        return -result[:ndcg10]
    end

    @info "Tuning start."
    best = fmin(invert_output, space, 100, logging_interval = -1)
    @info best

    @info "Evaluate best model."

    best_model = SLIM(best[:l1_ratio], best[:λminratio], convert(Int, best[:k]))
    result = evaluate_u2i(
        best_model,
        train_valid_table,
        test_table,
        metrics,
        10,
        col_item = :movieid,
        drop_history = true,
    )
    @info result

end

io = open("log_param_search_slim_ml00k.txt", "w+")

logger = SimpleLogger(io)
global_logger(logger)
main()

close(io)


text:
using AuditoryFilters, DSP, Compat, Compat.Test

@testset "ERB Filterbank" begin
    include("test_erb_fb.jl")
end

@testset "Modulation Filterbank" begin
    include("test_modulation_fb.jl")
end

@testset "Gammatonegram" begin
    include("test_gammatonegram.jl")
end


text:
include("kernel/kernel.jl")
include("tasks/tasks.jl")
include("generic.jl")

function subspace_mul!(st::AbstractVector{T}, comspace, U, subspace, offset=0) where T
    if size(U, 1) == 4
        return subspace_mul4x4!(st, comspace, U, subspace, offset)
    elseif size(U, 1) == 8
        return subspace_mul8x8!(st, comspace, U, subspace, offset)
    else
        # force compiler to specialize on comspace size
        indices = StrideArray{Int}(undef, (length(comspace), ))
        @simd ivdep for i in eachindex(indices)
            indices[i] = comspace[i] + 1
        end
        return subspace_mul_generic!(st, indices, U, subspace, offset)
    end
end

function subspace_mul!(st::AbstractMatrix, comspace, U::AbstractMatrix{T}, subspace, offset=0) where T
    # NOTE: avx version seems to be faster at 8 already
    if size(U, 1) == 4
        return subspace_mul4x4!(st, comspace, U, subspace, offset)
    else
        # force compiler to specialize on comspace size
        indices = StrideArray{Int}(undef, (length(comspace), ))
        @simd ivdep for i in eachindex(indices)
            indices[i] = comspace[i] + 1
        end
        return subspace_mul_generic!(st, indices, U, subspace, offset)
    end
end

# NOTE:
# expanded kernel is slower than avx kernel in multithreading
# thus we don't need to do any manual loop unroll.
function threaded_subspace_mul!(S::AbstractVecOrMat, comspace, U, subspace, offset=0)
    indices = StrideArray{Int}(undef, (length(comspace), ))
    @simd ivdep for i in eachindex(indices)
        indices[i] = comspace[i] + 1
    end
    return subspace_mul_generic!(S, indices, U, subspace, offset)
end

function subspace_mul4x4!(st::AbstractVector{T}, comspace, U, subspace, offset=0) where T
    Base.Cartesian.@nextract 4 indices i -> comspace[i] + 1

    Base.Cartesian.@nextract 4 U i->begin
        Base.Cartesian.@nextract 4 U_i j->U[i, j]
    end

    @inbounds for k in subspace
        Base.Cartesian.@nextract 4 idx i-> k + indices_i + offset

        Base.Cartesian.@nexprs 4 i -> begin
            y_i = zero(T)
            Base.Cartesian.@nexprs 4 j -> begin
                y_i += U_i_j * st[idx_j]
            end
        end

        Base.Cartesian.@nexprs 4 i -> begin
            st[idx_i] = y_i
        end
    end
    return st
end

function subspace_mul8x8!(st::AbstractVector{T}, comspace, U, subspace, offset=0) where T
    Base.Cartesian.@nextract 8 indices i -> comspace[i] + 1

    Base.Cartesian.@nextract 8 U i->begin
        Base.Cartesian.@nextract 8 U_i j->U[i, j]
    end

    @inbounds for k in subspace
        Base.Cartesian.@nextract 8 idx i-> k + indices_i + offset
        Base.Cartesian.@nexprs 8 i -> begin
            y_i = zero(T)
            Base.Cartesian.@nexprs 8 j -> begin
                y_i += U_i_j * st[idx_j]
            end
        end

        Base.Cartesian.@nexprs 8 i -> begin
            st[idx_i] = y_i
        end
    end
    return st
end

function subspace_mul4x4!(st::AbstractMatrix{T}, comspace, U::AbstractMatrix, subspace, offset=0) where T
    Base.Cartesian.@nextract 4 indices i -> comspace[i] + 1

    Base.Cartesian.@nextract 4 U i->begin
        Base.Cartesian.@nextract 4 U_i j->U[i, j]
    end

    @inbounds for k in subspace
        Base.Cartesian.@nextract 4 idx i-> k + indices_i + offset

        for b in 1:size(st, 1)
            Base.Cartesian.@nexprs 4 i -> begin
                y_i = zero(T)
                Base.Cartesian.@nexprs 4 j -> begin
                    y_i += U_i_j * st[b, idx_j]
                end
            end

            Base.Cartesian.@nexprs 4 i -> begin
                st[b, idx_i] = y_i
            end
        end
    end
    return st
end


text:
export MtlFence

const MTLFence = Ptr{MtFence}

mutable struct MtlFence
	handle::MTLFence
	device::MtlDevice
end

Base.unsafe_convert(::Type{MTLFence}, fen::MtlFence) = fen.handle

Base.:(==)(a::MtlFence, b::MtlFence) = a.handle == b.handle
Base.hash(ev::MtlFence, h::UInt) = hash(ev.handle, h)

function MtlFence(dev::MtlDevice)
	handle = mtDeviceNewFence(dev)
	obj = MtlFence(handle, dev)
	finalizer(unsafe_destroy!, obj)
	return obj
end

function unsafe_destroy!(fen::MtlFence)
	mtRelease(fen.handle)
end


text:
using Tables
using TableTraits
using IteratorInterfaceExtensions
using DataFrames

@testset "Tables interface tests" begin

    @testset "Chains" begin
        val = rand(1000, 8, 4)
        colnames = ["a", "b", "c", "d", "e", "f", "g", "h"]
        internal_colnames = ["c", "d", "e", "f", "g", "h"]
        chn = Chains(val, colnames, Dict(:internals => internal_colnames))

        @testset "Tables interface" begin
            @test Tables.istable(typeof(chn))
            @test Tables.columnaccess(typeof(chn))
            @test Tables.columns(chn) === chn
            @test Tables.columnnames(chn) ==
                  (:iteration, :chain, :a, :b, :c, :d, :e, :f, :g, :h)
            @test Tables.getcolumn(chn, :iteration) == [1:1000; 1:1000; 1:1000; 1:1000]
            @test Tables.getcolumn(chn, :chain) ==
                  [fill(1, 1000); fill(2, 1000); fill(3, 1000); fill(4, 1000)]
            @test Tables.getcolumn(chn, :a) == [
                vec(chn[:, :a, 1])
                vec(chn[:, :a, 2])
                vec(chn[:, :a, 3])
                vec(chn[:, :a, 4])
            ]
            @test_throws Exception Tables.getcolumn(chn, :j)
            @test Tables.getcolumn(chn, 1) == Tables.getcolumn(chn, :iteration)
            @test Tables.getcolumn(chn, 2) == Tables.getcolumn(chn, :chain)
            @test Tables.getcolumn(chn, 3) == Tables.getcolumn(chn, :a)
            @test_throws Exception Tables.getcolumn(chn, :i)
            @test_throws Exception Tables.getcolumn(chn, 11)
            @test Tables.rowaccess(typeof(chn))
            @test Tables.rows(chn) === chn
            @test length(Tables.rowtable(chn)) == 4000
            nt = Tables.rowtable(chn)[1]
            @test nt ==
                  (; (k => Tables.getcolumn(chn, k)[1] for k in Tables.columnnames(chn))...)
            @test nt == collect(Iterators.take(Tables.namedtupleiterator(chn), 1))[1]
            nt = Tables.rowtable(chn)[2]
            @test nt ==
                  (; (k => Tables.getcolumn(chn, k)[2] for k in Tables.columnnames(chn))...)
            @test nt == collect(Iterators.take(Tables.namedtupleiterator(chn), 2))[2]
            @test Tables.schema(chn) isa Tables.Schema
            @test Tables.schema(chn).names ===
                  (:iteration, :chain, :a, :b, :c, :d, :e, :f, :g, :h)
            @test Tables.schema(chn).types === (
                Int,
                Int,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
            )
            @test Tables.matrix(chn[:, :, 1])[:, 3:end] ≈ chn[:, :, 1].value
            @test Tables.matrix(chn[:, :, 2])[:, 3:end] ≈ chn[:, :, 2].value

            val = rand(1000, 2, 4)
            chn2 = Chains(val, ["iteration", "a"])
            @test_throws Exception Tables.columns(chn2)
            @test_throws Exception Tables.rows(chn2)
            @test_throws Exception Tables.schema(chn2)
            chn3 = Chains(val, ["chain", "a"])
            @test_throws Exception Tables.columns(chn3)
            @test_throws Exception Tables.rows(chn3)
            @test_throws Exception Tables.schema(chn3)
        end

        @testset "TableTraits interface" begin
            @test IteratorInterfaceExtensions.isiterable(chn)
            @test TableTraits.isiterabletable(chn)
            nt = collect(Iterators.take(IteratorInterfaceExtensions.getiterator(chn), 1))[1]
            @test nt ==
                  (; (k => Tables.getcolumn(chn, k)[1] for k in Tables.columnnames(chn))...)
            nt = collect(Iterators.take(IteratorInterfaceExtensions.getiterator(chn), 2))[2]
            @test nt ==
                  (; (k => Tables.getcolumn(chn, k)[2] for k in Tables.columnnames(chn))...)

            val = rand(1000, 2, 4)
            chn2 = Chains(val, ["iteration", "a"])
            @test_throws Exception IteratorInterfaceExtensions.getiterator(chn2)
            chn3 = Chains(val, ["chain", "a"])
            @test_throws Exception IteratorInterfaceExtensions.getiterator(chn3)
        end

        @testset "DataFrames.DataFrame constructor" begin
            @inferred DataFrame(chn)
            @test DataFrame(chn) isa DataFrame
            df = DataFrame(chn)
            @test Tables.columntable(df) == Tables.columntable(chn)
        end
    end

    @testset "ChainDataFrames" begin
        val = rand(1000, 8, 4)
        colnames = ["a", "b", "c", "d", "e", "f", "g", "h"]
        internal_colnames = ["c", "d", "e", "f", "g", "h"]
        chn = Chains(val, colnames, Dict(:internals => internal_colnames))
        cdf = describe(chn)[1]

        @testset "Tables interface" begin
            @test Tables.istable(typeof(cdf))
            @test Tables.columnaccess(typeof(cdf))
            @test Tables.columns(cdf) === cdf
            @test Tables.columnnames(cdf) == keys(cdf.nt)
            for (k, v) in pairs(cdf.nt)
                @test Tables.getcolumn(cdf, k) == v
            end
            @test Tables.getcolumn(cdf, 1) == Tables.getcolumn(cdf, keys(cdf.nt)[1])
            @test Tables.getcolumn(cdf, 2) == Tables.getcolumn(cdf, keys(cdf.nt)[2])
            @test_throws Exception Tables.getcolumn(cdf, :blah)
            @test_throws Exception Tables.getcolumn(cdf, length(cdf.nt) + 1)
            @test Tables.rowaccess(typeof(cdf))
            @test Tables.rows(cdf) === cdf
            @test length(Tables.rowtable(cdf)) == length(cdf.nt[1])
            @test Tables.columntable(cdf) == cdf.nt
            nt = Tables.rowtable(cdf)[1]
            @test nt == (; (k => v[1] for (k, v) in pairs(cdf.nt))...)
            @test nt == collect(Iterators.take(Tables.namedtupleiterator(cdf), 1))[1]
            nt = Tables.rowtable(cdf)[2]
            @test nt == (; (k => v[2] for (k, v) in pairs(cdf.nt))...)
            @test nt == collect(Iterators.take(Tables.namedtupleiterator(cdf), 2))[2]
            @test Tables.schema(cdf) isa Tables.Schema
            @test Tables.schema(cdf).names === keys(cdf.nt)
            @test Tables.schema(cdf).types === eltype.(values(cdf.nt))
        end

        @testset "TableTraits interface" begin
            @test IteratorInterfaceExtensions.isiterable(cdf)
            @test TableTraits.isiterabletable(cdf)
            nt = collect(Iterators.take(IteratorInterfaceExtensions.getiterator(cdf), 1))[1]
            @test nt == (; (k => v[1] for (k, v) in pairs(cdf.nt))...)
            nt = collect(Iterators.take(IteratorInterfaceExtensions.getiterator(cdf), 2))[2]
            @test nt == (; (k => v[2] for (k, v) in pairs(cdf.nt))...)
        end

        @testset "DataFrames.DataFrame constructor" begin
            @inferred DataFrame(cdf)
            df = DataFrame(cdf)
            @test df isa DataFrame
            @test Tables.columntable(df) == cdf.nt
        end
    end
end


text:
module FoosValidator

using SearchLight, SearchLight.Validation

function not_empty(field::Symbol, m::T, args::Vararg{Any})::ValidationResult where {T<:AbstractModel}
  isempty(getfield(m, field)) && return ValidationResult(invalid, :not_empty, "should not be empty")

  ValidationResult(valid)
end

function is_int(field::Symbol, m::T, args::Vararg{Any})::ValidationResult where {T<:AbstractModel}
  isa(getfield(m, field), Int) || return ValidationResult(invalid, :is_int, "should be an int")

  ValidationResult(valid)
end

function is_unique(field::Symbol, m::T, args::Vararg{Any})::ValidationResult where {T<:AbstractModel}
  findone(typeof(m); NamedTuple(field => getfield(m, field))... ) === nothing ||
    return ValidationResult(invalid, :is_unique, "already exists")

  ValidationResult(valid)
end

end


text:
function mandel(z)
    c = z
    maxiter = 80
    for n = 1:maxiter
        if abs(z) > 2
            return n-1
        end
        z = z^2 + c
    end
    return maxiter
end

function randmatstat(t)
    n = 5
    v = zeros(t)
    w = zeros(t)
    for i = 1:t
        a = randn(n,n)
        b = randn(n,n)
        c = randn(n,n)
        d = randn(n,n)
        P = [a b c d]
        Q = [a b; c d]
        v[i] = trace((P.'*P)^4)
        w[i] = trace((Q.'*Q)^4)
    end
    std(v)/mean(v), std(w)/mean(w)
end

max_stars_repo_name: MiincGu/DarkBint-syntax

max_issues_repo_name: MiincGu/DarkBint-syntax

max_forks_repo_name: MiincGu/DarkBint-syntax


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license

export read_step, read_step!, set_topology!, set_cell!, path
export MemoryTrajectory

__ptr(trajectory::Trajectory) = __ptr(trajectory.__handle)
__const_ptr(trajectory::Trajectory) = __const_ptr(trajectory.__handle)

"""
Opens the trajectory file at the given `path`.

The opening `mode` can be `'r'` for read, `'w'` for write or `'a'` for append,
and defaults to `'r'`. The optional `format` parameter give a specific file
format to use when opening the file.
"""
function Trajectory(path::AbstractString, mode::Char='r', format::AbstractString="")
    ptr = @__check_ptr(lib.chfl_trajectory_with_format(
        pointer(path), Int8(mode), pointer(format),
    ))
    return Trajectory(CxxPointer(ptr, is_const=false), nothing)
end

"""
Open a trajectory that read in-memory data as if it were a formatted file.

The `format` of the data is mandatory, and must match one of the known formats.
If `data` is `nothing`, this function creates a memory writer, the resulting
data can be retrieved with `buffer`. Else, this function creates a memory reader
using the given data.
"""
function MemoryTrajectory(format::AbstractString, data::DataBuffer=nothing)
    if data === nothing
        # memory writer
        ptr = @__check_ptr(lib.chfl_trajectory_memory_writer(pointer(format)))
    else
        # memory reader
        ptr = @__check_ptr(lib.chfl_trajectory_memory_reader(
            pointer(data), UInt64(sizeof(data)), pointer(format),
        ))
    end
    return Trajectory(CxxPointer(ptr, is_const=false), data)
end

"""
Apply the function `f` to the result of `Trajectory(args...)` and close the
resulting trajectory upon completion, similar to `open(f, args...)`.
"""
function Trajectory(f::Function, args...)
    tr = Trajectory(args...)
    try
        f(tr)
    finally
        close(tr)
    end
end

"""
Read the next step of the `trajectory`, and return the corresponding `Frame`.
"""
function Base.read(trajectory::Trajectory)
    @assert isopen(trajectory)
    frame = Frame()
    __check(lib.chfl_trajectory_read(
        __ptr(trajectory), __ptr(frame)
))
    return frame
end

"""
Read the next step of the `trajectory` data into the preexisting `Frame` structure.
"""
function Base.read!(trajectory::Trajectory, frame::Frame)
    @assert isopen(trajectory)
    __check(lib.chfl_trajectory_read(
        __ptr(trajectory), __ptr(frame)
    ))
end

"""
Read the given `step` of the `trajectory`, and return the corresponding
`Frame`.
"""
function read_step(trajectory::Trajectory, step::Integer)
    @assert isopen(trajectory)
    frame = Frame()
    __check(lib.chfl_trajectory_read_step(
        __ptr(trajectory), UInt64(step), __ptr(frame)
))
    return frame
end

"""
Read the given `step` of the `trajectory` into an preexisting `Frame` structure.
"""
function read_step!(trajectory::Trajectory, step::Integer, frame::Frame)
    @assert isopen(trajectory)
    __check(lib.chfl_trajectory_read_step(
        __ptr(trajectory), UInt64(step), __ptr(frame)
    ))
end

"""
Write the given `frame` to the `trajectory`.
"""
function Base.write(trajectory::Trajectory, frame::Frame)
    @assert isopen(trajectory)
    __check(lib.chfl_trajectory_write(
        __ptr(trajectory), __const_ptr(frame)
    ))
return nothing
end

"""
Obtain the memory buffer of a in-memory trajectory writer as a `Vector{UInt8}`,
without copying.

If more data is written to the trajectory, the buffer is invalidated.
"""
function Base.take!(trajectory::Trajectory)
    @assert isopen(trajectory)
    size = Ref{UInt64}(0)
    ptr = Ref{Ptr{UInt8}}()
    __check(lib.chfl_trajectory_memory_buffer(
        __const_ptr(trajectory), ptr, size
    ))
    return unsafe_wrap(Array{UInt8,1}, ptr[], (size[],); own=false)
end

"""
Set the `Topology` associated with a `trajectory`. This topology will be
used when reading and writing the file, replacing any topology in the file.
"""
function set_topology!(trajectory::Trajectory, topology::Topology)
    @assert isopen(trajectory)
    __check(lib.chfl_trajectory_set_topology(
        __ptr(trajectory), __const_ptr(topology)
    ))
return nothing
end

"""
Set the `Topology` associated with a `trajectory` by reading the first frame
of the file at `path`; and extracting the topology of this frame. The optional
`format` parameter can be used to specify the file format.
"""
function set_topology!(trajectory::Trajectory, path::AbstractString, format::AbstractString="")
    @assert isopen(trajectory)
    __check(lib.chfl_trajectory_topology_file(
        __ptr(trajectory), pointer(path), pointer(format)
    ))
return nothing
end

"""
Set the `cell` associated with a `trajectory`. This cell will be used when
reading and writing the file, replacing any unit cell in the file.
"""
function set_cell!(trajectory::Trajectory, cell::UnitCell)
    @assert isopen(trajectory)
    __check(lib.chfl_trajectory_set_cell(__ptr(trajectory), __const_ptr(cell)))
return nothing
end

"""
Get the number of steps (the number of frames) in a `trajectory`.
"""
function Base.size(trajectory::Trajectory)
    @assert isopen(trajectory)
    result = Ref{UInt64}(0)
    __check(lib.chfl_trajectory_nsteps(
        __const_ptr(trajectory), result
    ))
    return result[]
end

"""
Get the number of steps (the number of frames) in a `trajectory`.
"""
function Base.length(trajectory::Trajectory)
    size(trajectory)
end

"""
Get the path used to open a `trajectory`.
"""
function path(trajectory::Trajectory)
    @assert isopen(trajectory)
    return __call_with_growing_buffer(
        (buffer, size) -> __check(lib.chfl_trajectory_path(__const_ptr(trajectory), buffer, size))
    )
end

"""
Close a `trajectory`. This function flushes any buffer content to the hard
drive, and frees the associated memory. Necessary when running on the REPL to
finish writing.
"""
function Base.close(trajectory::Trajectory)
    # Manually free and set the pointer to NULL
    lib.chfl_trajectory_close(trajectory.__handle.__ptr)
    trajectory.__handle.__ptr = 0
return nothing
end

"""
Check if the `trajectory` is open.
"""
function Base.isopen(trajectory::Trajectory)
    return Int(trajectory.__handle.__ptr) != 0
end

# Iteration support
function Base.iterate(trajectory::Trajectory, step=0)
    if step >= length(trajectory)
        return nothing
    else
        return (read_step(trajectory, step), step + 1)
    end
end
Base.eltype(::Trajectory) = Frame


text:
using Parameters: @unpack, @with_kw

"""
    as(self, ::Type{T}) :: T

Manual Go-style type [embedding] as a replacement of "inheritance from
concrete type".

[embedding]: https://golang.org/doc/effective_go.html?#embedding
"""
as(self, ::Type{T}) where {T} = as(self.super, T) :: T
as(self::T, ::Type{T}) where {T} = self

abstract type AbstractContinuationSolver end

mutable struct ContinuationSolver{P <: AbstractContinuationProblem,
                                  C <: ContinuationCache,
                                  S <: ContinuationSolution,
                                  } <: AbstractContinuationSolver
    prob::P
    opts::ContinuationOptions
    cache::C
    sol::S
    i::Int
end

function ContinuationSolver(prob, opts)
    cache = ContinuationCache(prob, opts.h0, opts.direction)
    sol = ContinuationSolution(typeof(cache.u), prob)
    return ContinuationSolver(prob, opts, cache, sol, 0)
end

init(prob::AbstractContinuationProblem; kwargs...) =
    ContinuationSolver(prob, ContinuationOptions(; kwargs...))

solve(prob::AbstractContinuationProblem; kwargs...) =
    solve!(init(prob; kwargs...)).sol

LastPoint(solver::ContinuationSolver) =
    ContinuationLastPoint(solver.sol.sweeps[end], solver.cache.simple_bifurcation)
LastPoint(solver::AbstractContinuationSolver) =
    LastPoint(as(solver, ContinuationSolver))  # TODO: don't

function step!(solver::ContinuationSolver)
    predictor_corrector_step!(solver.cache, solver.opts)

    # Errors are now thrown in predictor_corrector_step!.  I need to
    # reconsider the error handling...
    if ! solver.cache.adaptation_success
        error("Failed to adapt steplength h.")
    end
    if ! solver.cache.corrector_success
        error("Failed in corrector loop.")
    end
    @debug "$(solver.i)-th step succeeded"

    record!(solver.sol, solver.cache)
    solver.i += 1
end

function record!(sol, cache)
    push_point!(sol, cache)
end

function step!(f, wrapper::AbstractContinuationSolver, max_steps)
    solver = as(wrapper, ContinuationSolver)
    cache = solver.cache
    cache.h = solver.opts.h0
    @progress_if solver.opts.verbose for _ in 1:max_steps
        step!(wrapper)
        f(LastPoint(wrapper))
        if ! isindomain(cache.u, cache.prob_cache)
            return true
        end
    end
    return false
end

struct SweepSetup{uType}
    direction::Int
    u0::uType
    past_points::Array{uType}
    max_steps::Int
end

SweepSetup(solver::ContinuationSolver;
           direction = solver.opts.direction,
           u0 = get_u0(solver.cache.prob_cache.prob),
           past_points = [],
           max_steps = solver.opts.max_samples,
           h = nothing,  # currently ignored  # TODO: use it?
           ) =
    SweepSetup(direction, u0, past_points, max_steps)

SweepSetup(solver::AbstractContinuationSolver; kwargs...) =
    SweepSetup(as(solver, ContinuationSolver); kwargs...)

function new_sweep!(solver::ContinuationSolver, setup::SweepSetup)
    @unpack direction, u0, past_points = setup
    cache = solver.cache
    cache.direction = direction
    cache.u = u0

    new_sweep!(solver.sol, direction)
    for u in past_points
        push_point!(solver.sol, u)
    end
    push_point!(solver.sol, u0)
end

function sweep!(f, solver::AbstractContinuationSolver; kwargs...)
    setup = SweepSetup(solver; kwargs...)
    new_sweep!(solver, setup)
    f(LastPoint(solver))
    step!(f, solver, setup.max_steps)
end

@with_kw struct NonRootException <: Exception
    msg::String = "Not a root"
    H::AbstractVector
    u::AbstractVector
end

function Base.showerror(io::IO, e::NonRootException)
    println(io, e.msg, ": norm(H) = ", norm(e.H))
    println(io, "H = ", e.H)
    println(io, "u = ", e.u)
end

function pre_solve!(wrapper::AbstractContinuationSolver)
    solver = as(wrapper, ContinuationSolver)
    opts = solver.opts
    cache = solver.cache
    H = residual!(cache.H, cache.u, cache.prob_cache)
    if ! isalmostzero(H, opts.atol)
        if opts.start_from_nearest_root
            @debug "Finding the nearest root..."
            cache.u = nearest_root!(cache, opts)
        else
            throw(NonRootException("Initial point is not a root", H, cache.u))
        end
    end
    return wrapper
end

solve!(solver::AbstractContinuationSolver) = solving!(donothing, solver)
donothing(_) = nothing

"""
    solving!(f, solver::AbstractContinuationSolver)

Same as `solve!(solve)` but call `f` after each continuation step with
the [`LastPoint`](@ref).

# Examples
```julia
solving!(solver) do point
    @show point.i_sweep point.i_point point.u[1] point.u[end]
end
```
"""
function solving!(f, wrapper::AbstractContinuationSolver)
    pre_solve!(wrapper)
    solver = as(wrapper, ContinuationSolver)
    opts = solver.opts
    cache = solver.cache

    @debug "Starting the first sweep..."
    u0 = copy(cache.u)
    sweep!(f, wrapper; u0=u0)
    if opts.bidirectional_first_sweep
        @debug "Starting the second sweep..."
        sweep!(f, wrapper; u0=u0, direction = solver.opts.direction * -1)
    end

    # TODO: Detect the case that the solution is isomorphic to the
    # circle.

    bifurcations = vcat(
        (s.simple_bifurcation for s in reverse(solver.sol.sweeps))...
    )
    for _ in 1:opts.max_branches
        if isempty(bifurcations)
            break
        end
        sbint = popfirst!(bifurcations)
        for (u0, u1, direction, h) in new_branches!(cache, opts, sbint)
            if ! (isindomain(u0, cache.prob_cache) &&
                  isindomain(u1, cache.prob_cache))
                # Stepped outside the domain.  Skip it.
                continue
            end
            @debug "Starting $(length(solver.sol.sweeps))-th sweep..."
            sweep!(f, wrapper;
                   u0 = u1,
                   past_points = [u0],
                   direction = direction,
                   h = h)
            append!(bifurcations, solver.sol.sweeps[end].simple_bifurcation)
        end
    end

    return wrapper
end

function residual(u::AbstractArray, cache::AbstractProblemCache)
    H = similar(u[1:end-1])
    return residual!(H, u, cache)
end

function residual(u::SVector, cache::AbstractProblemCache)
    return residual!(nothing, u, cache)
end

function residual_jacobian(u::AbstractArray, cache::AbstractProblemCache)
    H = similar(u[1:end-1])
    J = similar(H, (length(H), length(u)))
    return residual_jacobian!(H, J, u, cache)
end

function residual_jacobian(u::SVector, cache::AbstractProblemCache)
    return residual_jacobian!(nothing, nothing, u, cache)
end

function residual_jacobian!(cache::ContinuationCache, u::AbstractArray)
    @unpack H, J, prob_cache = cache
    H, J = residual_jacobian!(H, J, u, prob_cache)
    cache.u = u
    cache.H = H
    cache.J = J
    return H, J
end


text:
function CharactheristicExponent(mcProcess::FinancialMonteCarlo.BlackScholesProcess)
    σ = mcProcess.σ

    val_(v) = 0.5 * v * σ^2 * (1im - v)

    return val_
end

function CharactheristicExponent(mcProcess::FinancialMonteCarlo.MertonProcess)
    σ = mcProcess.σ
    sigma1 = mcProcess.σⱼᵤₘₚ
    mu1 = mcProcess.μⱼᵤₘₚ
    lam = mcProcess.λ

    val_(v) = -σ^2 * v^2 * 0.5 + lam * (exp(-sigma1^2 * v^2 / 2 + 1im * mu1 * v) - 1)

    return val_
end

function CharactheristicExponent(mcProcess::FinancialMonteCarlo.KouProcess)
    σ = mcProcess.σ
    lamp = mcProcess.λ₊
    lamm = mcProcess.λ₋
    p = mcProcess.p
    λ = mcProcess.λ

    val_(v) = -σ^2 * v^2 / 2 + 1im * v * λ * (p / (lamp - 1im * v) - (1 - p) / (lamm + 1im * v))

    return val_
end

function CharactheristicExponent(mcProcess::FinancialMonteCarlo.VarianceGammaProcess)
    σ = mcProcess.σ
    θ = mcProcess.θ
    κ = mcProcess.κ

    val_(v) = -1 / κ * log(1 + v^2 * σ^2 * κ / 2 - 1im * θ * κ * v)

    return val_
end

function CharactheristicExponent(mcProcess::FinancialMonteCarlo.NormalInverseGaussianProcess)
    σ = mcProcess.σ
    θ = mcProcess.θ
    κ = mcProcess.κ

    val_(v) = (1 - sqrt(1 + (v^2 * σ^2 - 2 * 1im * θ * v) * κ)) / κ

    return val_
end

function CharactheristicExponent(v::num_, mcProcess::proc, T::num2) where {num_ <: Number, proc <: FinancialMonteCarlo.AbstractMonteCarloProcess, num2 <: Number}
    ce = CharactheristicExponent(mcProcess)
    cf = ce(v) * T
    return cf
end

function CharactheristicFunction(mcProcess::proc, T::num2) where {proc <: FinancialMonteCarlo.AbstractMonteCarloProcess, num2 <: Number}
    # ! a function is returned
    ce = CharactheristicExponent(mcProcess)
    cf(v) = exp(ce(v) * T)
    return cf
end


text:
# dynamics functions, including the restricted three-body problem

export r3bp_hci_eom, twobody_eom, propagate_twobody, propagate_r3bp_hci

"""
Derivatives function for the restricted three-body problem in the
Heliocentric Inertial (HCI) reference frame xHCI is the stacked state
vector for the secondary planet + satellite: [xPlanet, xSat]
"""
function r3bp_hci_eom(t,
                      xHCI::Vector{T},
                      xdot,
                      mu2 = muEarth) where {T<:AbstractFloat}

    derivatives = zeros(12)
    derivatives[1:3] = xHCI[4:6]   # Planet's velocity
    derivatives[7:9] = xHCI[10:12] # Satellite velocity

    rPlanet = xHCI[1:3]
    rSat    = xHCI[7:9]
    derivatives[4:6]   = -muSun/(norm(rPlanet)^3) .* rPlanet
    derivatives[10:12] = -muSun/(norm(rSat)^3) .* rSat -
                            muEarth/(norm(rSat - rPlanet)^3) .* (rSat - rPlanet)
    xdot[:]=derivatives
end

"Derivatives function for the two-body problem"
function twobody_eom(t,
                     x::Vector{T},
                     xdot,
                     mu = muEarth) where {T<:AbstractFloat}

    xdot[1:3] = x[4:6]
    xdot[4:6] = -mu/(norm(x[1:3])^3) * x[1:3]
end

"""
Propagate the R3BP in the HCI frame vectors are expected in the HCI
frame

t should be an array of time values
"""
function propagate_r3bp_hci(x0Planet::Vector{T},
                            x0Sat::Vector{T},
                            t) where {T<:AbstractFloat}

    return Sundials.cvode(r3bp_hci_eom,
                   vcat(x0Planet, x0Sat),
                   t,
                   abstol=1e-12,reltol=1e-12)'
end

"Propagate the two-body problem"
function propagate_twobody(x0::Vector{T},
                           t) where {T<:AbstractFloat}

    return Sundials.cvode(twobody_eom,
                          x0,
                          t,
                          abstol=1e-12,reltol=1e-12)'
end


text:
function exportToPly(msh::Mesh, fn::String)
    vts = msh.vertices
    fcs = msh.faces
    nV = size(vts,1)
    nF = size(fcs,1)

    str = open(fn,"w")

    # write the header
    write(str,"ply\n")
    write(str,"format binary_little_endian 1.0\n")
    write(str,"element vertex $nV\n")
    write(str,"property float x\nproperty float y\nproperty float z\n")
    write(str,"element face $nF\n")
    write(str,"property list uchar int vertex_index\n")
    write(str,"end_header\n")

    # write the data
    for i = 1:nV
        v = vts[i]
        write(str,float32(v.e1))
        write(str,float32(v.e2))
        write(str,float32(v.e3))
    end

   for i = 1:nF
        f = fcs[i]
        write(str,uint8(3))
        write(str,int32(f.v1-1))
        write(str,int32(f.v2-1))
        write(str,int32(f.v3-1))
    end
    close(str)
end
export exportToPly

function exportToStl(msh::Mesh, fn::String)
    vts = msh.vertices
    fcs = msh.faces
    nV = size(vts,1)
    nF = size(fcs,1)

    str = open(fn,"w")

    # write the header
    write(str,"solid vcg\n")

    # write the data
    for i = 1:nF
        f = fcs[i]
        n = [0,0,0] # TODO: properly compute normal(f)
        txt = @sprintf "  facet normal %e %e %e\n" n[1] n[2] n[3]
        write(str,txt)
        write(str,"    outer loop\n")
        v = vts[f.v1]
        txt = @sprintf "      vertex  %e %e %e\n" v[1] v[2] v[3]
        write(str,txt)

        v = vts[f.v2]
        txt = @sprintf "      vertex  %e %e %e\n" v[1] v[2] v[3]
        write(str,txt)

        v = vts[f.v3]
        txt = @sprintf "      vertex  %e %e %e\n" v[1] v[2] v[3]
        write(str,txt)

        write(str,"    endloop\n")
        write(str,"  endfacet\n")
    end

    write(str,"endsolid vcg\n")
    close(str)
end
export exportToStl

# | Read a .2dm (SMS Aquaveo) mesh-file and construct a @Mesh@
function import2dm(file::String)
    parseNode(w::Array{String}) = Vertex(float64(w[3]), float64(w[4]), float64(w[5]))
    parseTriangle(w::Array{String}) = Face(int64(w[3]), int64(w[4]), int64(w[5]))
    # Qudrilateral faces are split up into triangles
    function parseQuad(w::Array{String})
        w[7] = w[3]                     # making a circle
        Face[Face(int64(w[i]), int64(w[i+1]), int64(w[i+2])) for i = [3,5]]
    end
    con = open(file, "r")
    nd =  Array(Vertex, 0)
    ele = Array(Face,0)
    for line = readlines(con)
        line = chomp(line)
        w = split(line)
        if w[1] == "ND"
            push!(nd, parseNode(w))
        elseif w[1] == "E3T"
            push!(ele,parseTriangle(w))
        elseif w[1] == "E4Q"
            append!(ele, parseQuad(w))
        else
            continue
        end
    end
    close(con)
    Mesh(nd,ele)
end
export import2dm

# | Write @Mesh@ to an IOStream
function exportTo2dm(f::IO,m::Mesh)
    function renderVertex(i::Int,v::Vertex)
        "ND $i $(v.x) $(v.y) $(v.z)\n"
    end
    function renderFace(i::Int, f::Face)
        "E3T $i $(f.v1) $(f.v2) $(f.v3) 0\n"
    end
    for i = 1:length(m.faces)
        write(f, renderFace(i, m.faces[i]))
    end
    for i = 1:length(m.vertices)
        write(f, renderVertex(i, m.vertices[i]))
    end
    nothing
end

# | Write a @Mesh@ to file in SMS-.2dm-file-format
function exportTo2dm(f::String,m::Mesh)
    con = open(f, "w")
    exportTo2dm(con, m)
    close(con)
    nothing
end
export exportTo2dm


text:
@testset "Label" begin
    @test l"abc" === Label{:abc}()
    @test l"abc" == l"abc"
    @test l"abc" != l"def"

    @test MVTs.label_to_index((:abc, :def), l"abc") == 1
    @test MVTs.label_to_index((:abc, :def), l"def") == 2
    @test_throws ErrorException MVTs.label_to_index((:abc, :def), l"ghi")
end