text:
using CombinatorialSpaces
using CombinatorialSpaces.ExteriorCalculus
using AlgebraicPetri

using Catlab
using Catlab.Graphics
using Catlab.Graphics.Graphviz
using Catlab.Programs

using Decapodes.Simulations
using Decapodes.Examples
using Decapodes.Diagrams
using Decapodes.Schedules
using Decapodes.OpenDiagrams
using Decapodes.PetriNets
include("deca_schema.jl")

draw_uwd(d) = to_graphviz(d, box_labels=:name, junction_labels=:variable,
            graph_attrs=Dict(:start => "2"))

draw_equations(eq) = begin
  to_graphviz(eq,
  node_labels=true, prog="neato",
  node_attrs=Dict(:shape=>"oval"),
  graph_attrs=Dict(:nodesep=>"4.0"))
end

@present Flow2DQuantities <: ExtendedFreeExtCalc2D begin
  X::Space

  k::Hom(DualForm1(X), DualForm1(X))    # diffusivity (usually scalar multiplication)
  k⁻¹::Hom(DualForm1(X), DualForm1(X))  # inverse of diffusivity
  neg₀::Hom(Form0(X), Form0(X))
  L₀::Hom(Form1(X)⊗DualForm2(X), DualForm2(X))
end

#####################
# Example DECAPODES #
#####################

Diffusion = @decapode Flow2DQuantities begin
  T::Form0{X}
  ϕ::DualForm1{X}

  # Fick's first law
  ϕ ==  k(⋆₁{X}(d₀{X}(T)))
  # Diffusion equation
  ∂ₜ{Form0{X}}(T) == ⋆₀⁻¹{X}(dual_d₁{X}(ϕ))
end

open("diffusion.svg", "w") do io
  run_graphviz(io, draw_equations(Diffusion), format="svg")
end

UncompilableDiffusion = @decapode Flow2DQuantities begin
  T::Form0{X}
  ϕ::DualForm1{X}

  # Fick's first law
  k⁻¹(ϕ) ==  ⋆₁{X}(d₀{X}(T))
  # Diffusion equation
  ∂ₜ{Form0{X}}(T) == ⋆₀⁻¹{X}(dual_d₁{X}(ϕ))
end

open("uncomp_diffusion.svg", "w") do io
  run_graphviz(io, draw_equations(UncompilableDiffusion), format="svg")
end

##############################
# Physics Theory Composition #
##############################

Diffusion = @decapode Flow2DQuantities begin
  (T, Ṫ)::Form0{X}
  ϕ::DualForm1{X}

  # Fick's first law
  ϕ ==  k(⋆₁{X}(d₀{X}(T)))
  Ṫ == ⋆₀⁻¹{X}(dual_d₁{X}(ϕ))
end

Advection = @decapode Flow2DQuantities begin
  (T, Ṫ)::Form0{X}
  V::Form1{X}
  Ṫ == neg₀(⋆₀⁻¹{X}(L₀(V, ⋆₀{X}(T))))
end

Dynamics = @decapode Flow2DQuantities begin
  (X, Ẋ)::Form0{X}
  ∂ₜ{Form0{X}}(X) == Ẋ
end

dynam_physics = @relation (X, Ẋ) begin
  physics(X, Ẋ)
  dynamics(X, Ẋ)
end

DiffusionDynam = oapply(dynam_physics,
                  [OpenDiagram(Diffusion, [:T, :Ṫ]),
                   OpenDiagram(Dynamics, [:X, :Ẋ])])

AdvectionDynam = oapply(dynam_physics,
                   [OpenDiagram(Advection, [:T, :Ṫ]),
                    OpenDiagram(Dynamics, [:X, :Ẋ])])

draw_equations(Diffusion)

mkpath("simple_composition")

open("simple_composition/diffusion.svg", "w") do io
  run_graphviz(io, draw_equations(Diffusion), format="svg")
end

open("simple_composition/advection.svg", "w") do io
  run_graphviz(io, draw_equations(Advection), format="svg")
end

open("simple_composition/dynamics.svg", "w") do io
  run_graphviz(io, draw_equations(Dynamics), format="svg")
end

open("simple_composition/comp_pattern.svg", "w") do io
  run_graphviz(io, draw_uwd(dynam_physics), format="svg")
end

open("simple_composition/diff_dynam.svg", "w") do io
  run_graphviz(io, draw_equations(DiffusionDynam), format="svg")
end

open("simple_composition/adv_dynam.svg", "w") do io
  run_graphviz(io, draw_equations(AdvectionDynam), format="svg")
end

mkpath("multi_composition")

Superposition = @decapode Flow2DQuantities begin
  (Ẋ, Ẋ₁, Ẋ₂)::Form0{X}
  Ẋ == Ẋ₁ + Ẋ₂
end

two_phys = @relation (X, Ẋ₁, Ẋ₂) begin
  physics1(X, Ẋ₁)
  physics2(X, Ẋ₂)
  superposition(Ẋ, Ẋ₁, Ẋ₂)
end

AdvectionDiffusion = oapply(two_phys,
                   [OpenDiagram(Advection, [:T, :Ṫ]),
                    OpenDiagram(Diffusion, [:T, :Ṫ]),
                    OpenDiagram(Superposition, [:Ẋ, :Ẋ₁, :Ẋ₂])])

AdvectionDiffusionDynam = oapply(dynam_physics,
  [AdvectionDiffusion,
  OpenDiagram(Dynamics, [:X, :Ẋ])])

open("multi_composition/adv_diff_pattern.svg", "w") do io
  run_graphviz(io, draw_uwd(two_phys), format="svg")
end

open("multi_composition/advection_diffusion.svg", "w") do io
  run_graphviz(io, draw_equations(AdvectionDiffusion), format="svg")
end

open("multi_composition/superposition.svg", "w") do io
  run_graphviz(io, draw_equations(Superposition), format="svg")
end

open("multi_composition/comp_pattern.svg", "w") do io
  run_graphviz(io, draw_uwd(dynam_physics), format="svg")
end

open("multi_composition/diff_adv_dynam.svg", "w") do io
  run_graphviz(io, draw_equations(AdvectionDiffusionDynam), format="svg")
end
#=
enz_deg = LabelledReactionNet{Float64, Float64}([:AB=>0.0, :E=>0.0, :A=>0.0, :B=>0.0], (:deg=>0.05)=>((:AB,:E)=>(:A, :B, :E)))
Graph(enz_deg)

expand_pres!(Flow2DQuantities, enz_deg)
EnzymesDyn = pn2dec(Flow2DQuantities, enz_deg)
draw_equations(EnzymesDyn)

enz_phys = @relation (A, B, E, AB) begin
  physicsA(A, Ȧₚ)
  physicsB(B, Ḃₚ)
  physicsAB(AB, ∂ₜABₚ)
  enzyme(A, B, E, AB, Ȧₑ, Ḃₑ, ∂ₜABₑ)
  superposition(Ȧ, Ȧₑ, Ȧₚ)
  superposition(Ḃ, Ḃₑ, Ḃₚ)
  superposition(∂ₜAB, ∂ₜABₑ, ∂ₜABₚ)
end

draw_uwd(enz_phys)
AdvectionDiffusion = oapply(two_phys,
                   [AdvectionDiffusion,
                    AdvectionDiffusion,
                    AdvectionDiffusion,
                    OpenDiagram(EnzymesDyn, [:A, :B, :E, :AB, :Ḃ, :Ḃ, :∂ₜAB]),
                    OpenDiagram(Superposition, [:Ẋ, :Ẋ₁, :Ẋ₂]),
                    OpenDiagram(Superposition, [:Ẋ, :Ẋ₁, :Ẋ₂]),
                    OpenDiagram(Superposition, [:Ẋ, :Ẋ₁, :Ẋ₂])])

mkpath("chem_composition")=#

max_stars_repo_name: AlgebraicJulia/DECAPODES-Benchmarks

max_issues_repo_name: AlgebraicJulia/DECAPODES-Benchmarks

max_forks_repo_name: AlgebraicJulia/DECAPODES-Benchmarks


text:
export tap_on_subscribe

import Base: show

"""
    tap_on_subscribe(tapFn::F) where { F <: Function }

Creates a tap operator, which performs a side effect on the subscription on the source Observable, but return an Observable that is identical to the source.

# Arguments
- `tapFn::Function`: side-effect tap function with `() -> Nothing` signature

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

# Examples
```jldoctest
using Rocket

source = from([ 1, 2, 3 ])
subscribe!(source |> tap_on_subscribe(() -> println("Someone subscribed")), logger())
;

# output

Someone subscribed
[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Data: 3
[LogActor] Completed

```

See also: [`tap`](@ref), [`tap_on_complete`](@ref), [`logger`](@ref)
"""
tap_on_subscribe(tapFn::F) where { F <: Function } = TapOnSubscribeOperator{F}(tapFn)

struct TapOnSubscribeOperator{F} <: InferableOperator
    tapFn :: F
end

operator_right(operator::TapOnSubscribeOperator, ::Type{L}) where L = L

function on_call!(::Type{L}, ::Type{L}, operator::TapOnSubscribeOperator{F}, source) where { L, F }
    return proxy(L, source, TapOnSubscribeProxy{F}(operator.tapFn))
end

struct TapOnSubscribeProxy{F} <: SourceProxy
    tapFn :: F
end

source_proxy!(::Type{L}, proxy::TapOnSubscribeProxy{F}, source::S) where { L, S, F } = TapOnSubscribeSource{L, S, F}(proxy.tapFn, source)

struct TapOnSubscribeSource{L, S, F} <: Subscribable{L}
    tapFn  :: F
    source :: S
end

function on_subscribe!(source::TapOnSubscribeSource, actor)
    source.tapFn()
    return subscribe!(source.source, actor)
end

Base.show(io::IO, ::TapOnSubscribeOperator)          = print(io, "TapOnSubscribeOperator()")
Base.show(io::IO, ::TapOnSubscribeProxy)             = print(io, "TapOnSubscribeProxy()")
Base.show(io::IO, ::TapOnSubscribeSource{L}) where L = print(io, "TapOnSubscribeSource($L)")

ext: jl


text:
runner_path = joinpath(dirname(@__FILE__), "../src/runner.jl")
julia_exe = Base.julia_cmd()

# FIXME: ugly way of making the test pass - come up with a way to really test test selection
"Returns an expected output from Base.Test to check if tests are selected as expected"
tests_count_indicator(passed, total) =
  v"0.5" <= VERSION < v"0.6.0-dev.2370" ? "test set    |    $passed      $total" :
                                           "test set      |    $passed      $total"

cd(joinpath(dirname(@__FILE__), "test_package")) do
  @test contains(readstring(`$julia_exe $runner_path`),
                 tests_count_indicator(3,3))

  @test contains(readstring(`$julia_exe $runner_path runtests.jl/one`),
                 tests_count_indicator(1,1))

  @test contains(readstring(`$julia_exe $runner_path runtests.jl`),
                 tests_count_indicator(3,3))

  @test contains(readstring(`$julia_exe $runner_path runtests.jl/one runtests.jl/two`),
                 tests_count_indicator(2,2))
end


text:
#########################################################
# Config{N}(args...) --> Config(f, args..., Chunk{N}()) #
#########################################################

function (::Type{GradientConfig{N}})(x) where N
    msg = "GradientConfig{N}(x) is deprecated; use GradientConfig(nothing, x, Chunk{N}()) instead."
    Base.depwarn(msg, :GradientConfig)
    return GradientConfig(nothing, x, Chunk{N}())
end

function (::Type{JacobianConfig{N}})(x) where N
    msg = "JacobianConfig{N}(x) is deprecated; use JacobianConfig(nothing, x, Chunk{N}()) instead."
    Base.depwarn(msg, :JacobianConfig)
    return JacobianConfig(nothing, x, Chunk{N}())
end

function (::Type{JacobianConfig{N}})(y, x) where N
    msg = "JacobianConfig{N}(y, x) is deprecated; use JacobianConfig(nothing, y, x, Chunk{N}()) instead."
    Base.depwarn(msg, :JacobianConfig)
    return JacobianConfig(nothing, y, x, Chunk{N}())
end

function (::Type{HessianConfig{N}})(x) where N
    msg = "HessianConfig{N}(x) is deprecated; use HessianConfig(nothing, x, Chunk{N}()) instead."
    Base.depwarn(msg, :HessianConfig)
    return HessianConfig(nothing, x, Chunk{N}())
end

function (::Type{HessianConfig{N}})(out, x) where N
    msg = "HessianConfig{N}(out, x) is deprecated; use HessianConfig(nothing, out, x, Chunk{N}()) instead."
    Base.depwarn(msg, :HessianConfig)
    return HessianConfig(nothing, out, x, Chunk{N}())
end

function MultithreadConfig(cfg::AbstractConfig)
    msg = "MultithreadConfig(cfg) is deprecated; use cfg instead (ForwardDiff no longer implements experimental multithreading)."
    Base.depwarn(msg, :MultithreadConfig)
    return cfg
end


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

#
# This implements Agha's example 3.2.1
#

using Actors
import Actors: spawn, newLink

mutable struct StackNode{T,L}
    content::T
    link::L
end

struct Pop{L}
    customer::L
end

struct Push{T}
    content::T
end

forwarder = send

function (sn::StackNode)(msg::Pop)
    isnothing(sn.content) || become(forwarder, sn.link)
    send(msg.customer, sn.content)
end
(sn::StackNode)(msg::Push) = become(StackNode(msg.content, spawn(sn)))

mystack = spawn(StackNode(nothing, newLink()))
response = newLink()

for i ∈ 1:5
    send(mystack, Push(i))
end

for i ∈ 1:5
    send(mystack, Pop(response))
    println(receive(response))
end


text:
function OceanLoadingCoefficients(arg0::JString, arg1::GeodeticPoint, arg2::Vector{Vector{Tide}}, arg3::Vector{Vector{jdouble}}, arg4::Vector{Vector{jdouble}}, arg5::Vector{Vector{jdouble}}, arg6::Vector{Vector{jdouble}}, arg7::Vector{Vector{jdouble}}, arg8::Vector{Vector{jdouble}})
    return OceanLoadingCoefficients((JString, GeodeticPoint, Vector{Vector{Tide}}, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8)
end

function get_nb_species(obj::OceanLoadingCoefficients)
    return jcall(obj, "getNbSpecies", jint, ())
end

function get_nb_tides(obj::OceanLoadingCoefficients, arg0::jint)
    return jcall(obj, "getNbTides", jint, (jint,), arg0)
end

function get_site_location(obj::OceanLoadingCoefficients)
    return jcall(obj, "getSiteLocation", GeodeticPoint, ())
end

function get_site_name(obj::OceanLoadingCoefficients)
    return jcall(obj, "getSiteName", JString, ())
end

function get_south_amplitude(obj::OceanLoadingCoefficients, arg0::jint, arg1::jint)
    return jcall(obj, "getSouthAmplitude", jdouble, (jint, jint), arg0, arg1)
end

function get_south_phase(obj::OceanLoadingCoefficients, arg0::jint, arg1::jint)
    return jcall(obj, "getSouthPhase", jdouble, (jint, jint), arg0, arg1)
end

function get_tide(obj::OceanLoadingCoefficients, arg0::jint, arg1::jint)
    return jcall(obj, "getTide", Tide, (jint, jint), arg0, arg1)
end

function get_west_amplitude(obj::OceanLoadingCoefficients, arg0::jint, arg1::jint)
    return jcall(obj, "getWestAmplitude", jdouble, (jint, jint), arg0, arg1)
end

function get_west_phase(obj::OceanLoadingCoefficients, arg0::jint, arg1::jint)
    return jcall(obj, "getWestPhase", jdouble, (jint, jint), arg0, arg1)
end

function get_zenith_amplitude(obj::OceanLoadingCoefficients, arg0::jint, arg1::jint)
    return jcall(obj, "getZenithAmplitude", jdouble, (jint, jint), arg0, arg1)
end

function get_zenith_phase(obj::OceanLoadingCoefficients, arg0::jint, arg1::jint)
    return jcall(obj, "getZenithPhase", jdouble, (jint, jint), arg0, arg1)
end


text:
using ImageMagick
using RegisterQD
using Test

include("util.jl")
include("qd_random.jl")
include("qd_standard.jl")
include("gridsearch.jl")
include("initial_tfm.jl")


text:
include("network.jl")

spike_train = rand(Bool, 3,5)

nn = Network(;input_size=3, output_size=2, hidden_neurons=3, runtime=5)
nn.construct(nn; n₁=1, n₂=3)
nn.construct(nn; n₁=1, n₂=2)

nn.activate(nn, spike_train; θ=0.1, r=1, η=0.1, γᵥ=0.3, γₚ=0.3, γₘ=0.3, A₊=2.0, A₋=1.0, w₊=1.5, w₋=-0.8)

max_stars_repo_name: Louisenlund-plusMINT/NeuronSim

max_issues_repo_name: Louisenlund-plusMINT/NeuronSim

max_forks_repo_name: Louisenlund-plusMINT/NeuronSim


text:
using BenchmarkTools
using LinearAlgebra
using Random
using Suppressor

include("utils/utils.jl")

# Command line arguments
nARGS = length(ARGS)
file = ARGS[1]
solver = ARGS[2]
bench = ARGS[3]
if bench == "nthreads"
	n = parse(Int, ARGS[4])
elseif bench == "size"
	N0 = parse(Int, ARGS[4])
	ppomag = parse(Int, ARGS[5])
	nrun = parse(Int, ARGS[6])
else
	throw(ArgumentError(bench))
end
nthreads = parse(Int, ARGS[nARGS])

# Benchmark parameters
BenchmarkTools.DEFAULT_PARAMETERS.samples = parse(Int, ARGS[length(ARGS)]) # 10000 by default
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 60*60 # 5 by default

# Run benchmarks for increasing number of threads
function run_nthreads(file::String, solver::String, n::Int, nthreads::Int)

	# Initialize arrays
	Random.seed!(1)
	a = rand()
	x = rand(n)
	y = rand(n)

	# Run benchmarks and get runtimes
	if solver == "dot"
		benchoutput = @capture_out @mybtimes dot($x, $y)
	elseif solver == "axpy"
		benchoutput = @capture_out @mybtimes axpy!($a, $x, $y)
	end

	# Write benchmark runtimes to output file
	open(file, "a") do io
		write(io, string(nthreads, " ", n, benchoutput))
	end

end

# Run benchmarks for increasing vector sizes
function run_size(file::String, solver::String, N0::Int, ppomag::Int,
	nrun::Int, nthreads::Int)

	# Get vector sizes
	vecsizes = get_vsizes(N0, ppomag, nrun)

	# Write header to output file
	write(file, "nthreads size min med mean max\n")

	# Terminal output
	print("size = ")

	# Run benchmarks for all vector sizes
	for n in vecsizes

		# Terminal output
		print(n, "..")

		# Initialize arrays
		Random.seed!(1)
		a = rand()
		x = rand(n)
		y = rand(n)

		# Run benchmarks and get runtimes
		if solver == "dot"
			benchoutput = @capture_out @mybtimes dot($x, $y)
		elseif solver == "axpy"
			benchoutput = @capture_out @mybtimes axpy!($a, $x, $y)
		end

		# Write vectorsize and benchmark runtimes to output file
		open(file, "a") do io
			write(io, string(nthreads, " ", n, benchoutput))
		end

	end

	# Terminal output
	println()

end

# Run from command line
function main()
	if bench == "nthreads"
		run_nthreads(file, solver, n, nthreads)
	elseif bench == "size"
		run_size(file, solver, N0, ppomag, nrun, nthreads)
	else
		throw(ArgumentError(bench))
	end
end

main()

max_stars_repo_name: jkrch/JuliaParGMRES-BenchCLI

max_issues_repo_name: jkrch/JuliaParGMRES-BenchCLI

max_forks_repo_name: jkrch/JuliaParGMRES-BenchCLI


text:
using CompileBot

botconfig = BotConfig(
    "Plots",
    version = [v"1.6", v"1.7"],  # <<< keep versions in sync with .github/workflows/SnoopCompile.yml
    # else_version = v"1.8",
)


text:
using Pkg; Pkg.activate("."); Pkg.instantiate()
top = pwd()
cd("./examples")
include(joinpath(pwd(), "./flint_rock3_examples.jl"));
d2 =  flint_rock3_examples.free_vibration_solver_w_remeshing_ren(6, 2)
@show d2["eigenvalues"]
cd(top)

max_stars_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids

max_issues_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids

max_forks_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids


text:
include("../common.jl")

build_libcurl(ARGS, "CURL")


text:
# Copyright 2016 Martin Holters
# See accompanying license file.

using ACME

# Model of "Der Birdie", see
# http://diy.musikding.de/wp-content/uploads/2013/06/birdieschalt.pdf
# for schematics.
# Note: Contrary to the schematics, in the real circuit as well as in this
# model, C5 and D1 are connected directly to the power supply; R6 is in series
# with LED1 only, which has been omitted from the model, though.

function birdie(::Type{Circuit}; vol=nothing)
    c1 = capacitor(2.2e-9)
    c3 = capacitor(2.2e-9)
    c5 = capacitor(100e-6)

    r1 = resistor(1e6)
    r2 = resistor(43e3)
    r3 = resistor(430e3)
    r4 = resistor(390)
    r5 = resistor(10e3)
    if vol == nothing
        p1 = potentiometer(100e3)
    else
        p1 = potentiometer(100e3, vol)
    end

    d1 = diode(is=350e-12, η=1.6);
    t1 = bjt(:npn, isc=154.1e-15, ise=64.53e-15, ηc=1.10, ηe=1.06, βf=500, βr=12)

    j1 = voltagesource() # input
    j2 = voltageprobe() # output
    j3 = voltagesource(9) # 9V power supply

    circ = Circuit()
    connect!(circ, j3[:+], c5[2], d1[:-], :vcc)
    connect!(circ, j3[:-], c5[1], d1[:+], :gnd)
    connect!(circ, j1[:-], r1[2], :gnd)
    connect!(circ, j1[:+], r1[1], c1[1])
    connect!(circ, c1[2], r2[1], r3[1], t1[:base])
    connect!(circ, r2[2], :gnd)
    connect!(circ, r3[2], :vcc)
    connect!(circ, r4[1], t1[:emitter])
    connect!(circ, r4[2], :gnd)
    connect!(circ, r5[1], c3[1], t1[:collector])
    connect!(circ, r5[2], :vcc)
    connect!(circ, c3[2], p1[3])
    connect!(circ, p1[2], j2[:+])
    connect!(circ, p1[1], j2[:-], :gnd)

    return circ
end

birdie{T<:DiscreteModel}(::Type{T}=DiscreteModel; vol=nothing, fs=44100) =
    T(birdie(Circuit, vol=vol), 1/fs)


text:
# acc plotting
using PGFPlots
using QHull
using JLD2
using LazySets
using LinearAlgebra

include("../../../../models/problems.jl")
include("../../../../OverApprox/src/overapprox_nd_relational.jl")
include("../../../../OverApprox/src/overt_parser.jl")
include("../../../../MIP/src/overt_to_mip.jl")
include("../../../../MIP/src/mip_utils.jl")
include("../../../../models/tora/tora.jl")

controller = ARGS[1] #"smallest"

if controller == "smallest"
    JLD2.@load "examples/jmlr/data/tora/tora_reachability_smallest_controller_data_1step.jld2"
    one_step_state_sets = symbolic_state_sets
    one_step_meas_sets = symbolic_meas_sets

    JLD2.@load "examples/jmlr/data/tora/tora_reachability_smallest_controller_data.jld2"

elseif controller == "smaller"
    JLD2.@load "examples/jmlr/data/tora/tora_reachability_smaller_controller_data_1step.jld2"
    one_step_state_sets = symbolic_state_sets
    one_step_meas_sets = symbolic_meas_sets

    JLD2.@load "examples/jmlr/data/tora/tora_reachability_smaller_controller_data.jld2"
elseif controller == "big"
    JLD2.@load "examples/jmlr/data/tora/tora_reachability_big_controller_data_1step.jld2"
    one_step_state_sets = symbolic_state_sets
    one_step_meas_sets = symbolic_meas_sets

    JLD2.@load "examples/jmlr/data/tora/tora_reachability_big_controller_data.jld2"
end

if (!@isdefined concretization_intervals) || length(concretization_intervals) == query.ntime
    concretization_intervals = [length(s) for s in concrete_state_sets] .- 1
end

concrete_state_sets = vcat(concrete_state_sets...)
symbolic_state_sets = vcat(symbolic_state_sets...)

# monte carlo simulate
mc_state_sets, xvec, x0, mc_meas_sets, yvec, y0 = monte_carlo_simulate(query, input_set, n_sim=1000000)

# clean up sets
init, reachable_state_sets = clean_up_sets(concrete_state_sets, symbolic_state_sets, concretization_intervals)

####################################
##### plot 1: intervals of x1 ######
####################################
dim=[1]
sym_sets = get_interval_subsets(reachable_state_sets, dim)
conc_sets = get_interval_subsets(one_step_state_sets, dim)
mc_sets = get_interval_subsets(mc_state_sets, dim)
deleteat!(mc_sets, 1) # pop init set off start

# plot styles
# define_color("concrete_color", 0xc9a8ff)
# define_color("symbolic_color", 0xff5c57)
# define_color("mc_color", 0x96eaff)
define_color("concrete_color", 0x139EAB)
define_color("symbolic_color", 0x9BFF85)
define_color("mc_color", 0x38caff) #96eaff)
# define_color("concrete_color", 0x648FFF)
# define_color("symbolic_color", 0xDC267F)
# define_color("mc_color", 0xFE6100)
conc_style = "solid, concrete_color, line width=1pt, mark=none"
sym_style = "solid, symbolic_color, line width = 2 pt, mark=none"
mc_style = "solid, mc_color, line width=3pt, mark=none"
input_style = "dashed, blue, mark=none"
############################################################
# Plot: x1 and property
############################################################
fig = PGFPlots.Axis(style="width=10cm, height=10cm", ylabel="\$x_1\$", xlabel="timesteps", title="TORA, State Sets")

# AVOID SETS
push!(fig, PGFPlots.Plots.Linear([-1., 16], [-2., -2.], style="red, dashed", mark="none", legendentry="Unsafe Set"))

# plot init set
dims=[1]
input = [low(input_set)[dims]..., high(input_set)[dims]...]
push!(fig, PGFPlots.Plots.Linear([0.,0.], input, style=conc_style))
push!(fig, PGFPlots.Plots.Linear([0.,0.], input, style=sym_style))
push!(fig, PGFPlots.Plots.Linear([0.,0.], input, style=mc_style))

for t in 1:query.ntime
    timestep = Float64[t,t]
    if t == query.ntime
        push!(fig, PGFPlots.Plots.Linear( timestep, conc_sets[t], style=conc_style, legendentry="Concrete Sets"))
        push!(fig, PGFPlots.Plots.Linear( timestep, sym_sets[t], style=sym_style, legendentry="OVERT Hybrid Symbolic Sets"))
        push!(fig, PGFPlots.Plots.Linear( timestep, mc_sets[t], style=mc_style, legendentry="Monte Carlo Simulations"))
    else
        push!(fig, PGFPlots.Plots.Linear( timestep, conc_sets[t], style=conc_style))
        push!(fig, PGFPlots.Plots.Linear( timestep, sym_sets[t], style=sym_style))
        push!(fig, PGFPlots.Plots.Linear( timestep, mc_sets[t], style=mc_style))
    end
end

# second unsafe set
push!(fig, PGFPlots.Plots.Linear([-1., 16], [2., 2.], style="red, dashed", mark="none"))

fig.legendStyle =  "at={(1.05,1.0)}, anchor=north west"

save("examples/jmlr/plots/tora/tora_$(controller)_x1.tex", fig)
save("examples/jmlr/plots/tora/tora_$(controller)_x1.pdf", fig)

####################################
##### plot 2: boxes ######
####################################
conc_style_solid = "solid, concrete_color, thick, mark=none, fill=concrete_color"
conc_style_transparent = conc_style_solid*", fill opacity=0.5"
sym_style_solid = "solid, symbolic_color, thick, mark=none, fill=symbolic_color"
sym_style_transparent =  sym_style_solid*", fill opacity=0.5"
mc_style_solid = "solid, mc_color, thick, mark=none, fill=mc_color"
mc_style_transparent = mc_style_solid*", fill opacity=0.5"

fig = PGFPlots.Axis(style="width=10cm, height=10cm", ylabel="\$x_2\$", xlabel="\$x_3\$", title="TORA, State Sets")

# hack for legend
push!(fig, PGFPlots.Plots.Linear([-0.7, -0.7], [0.2, 0.2], style=mc_style_solid, legendentry="Monte Carlo Simulations"))

# plot init set in both concrete and hybrid colors
dims=[2,3]
inputx, inputy = get_rectangle(input_set, dims)
push!(fig, PGFPlots.Plots.Linear(inputx, inputy, style=conc_style_transparent, legendentry="Concrete Sets"))
push!(fig, PGFPlots.Plots.Linear(inputx, inputy, style=sym_style_transparent, legendentry="OVERT Hybrid Symbolic Sets"))

for t in 1:query.ntime
    if t == query.ntime
        push!(fig, PGFPlots.Plots.Linear( get_rectangle(one_step_state_sets[t], dims)..., style=conc_style_solid))
        push!(fig, PGFPlots.Plots.Linear( get_rectangle(reachable_state_sets[t], dims)..., style=sym_style_solid))
        # push!(fig, PGFPlots.Plots.Linear( get_rectangle(mc_state_sets[t+1], dims)..., style=mc_style, legendentry="Monte Carlo Simulations"))
        # plot mc_sim points
        points = xvec[:, t, dims]
        #plot hull
        border_idx = chull(points).vertices
        x = points[border_idx, 1]
        y = points[border_idx, 2]
        push!(fig, PGFPlots.Plots.Linear([x..., x[1]], [y..., y[1]], style=mc_style_solid))
    else
        push!(fig, PGFPlots.Plots.Linear( get_rectangle(one_step_state_sets[t], dims)..., style=conc_style_transparent))
        push!(fig, PGFPlots.Plots.Linear( get_rectangle(reachable_state_sets[t], dims)..., style=sym_style_transparent))
        points = xvec[:, t, dims]
        #plot hull
        border_idx = chull(points).vertices
        x = points[border_idx, 1]
        y = points[border_idx, 2]
        push!(fig, PGFPlots.Plots.Linear([x..., x[1]], [y..., y[1]], style=mc_style_solid))
    end
end

fig.legendStyle =  "at={(1.05,1.0)}, anchor=north west"

PGFPlots.save("examples/jmlr/plots/tora/tora_$(controller)_x23.tex", fig)
PGFPlots.save("examples/jmlr/plots/tora/tora_$(controller)_x23.pdf", fig)


text:
using ProjectiveGeometricAlgebra2d
using Test
using ProjectiveGeometricAlgebra2d: PGA2d
using ProjectiveGeometricAlgebra2d: wedge, anti_wedge, geomul
using ProjectiveGeometricAlgebra2d: dual, reverse
using ProjectiveGeometricAlgebra2d: euclidean_point, euclidean_line_abc
using ProjectiveGeometricAlgebra2d: pga_point, pga_line
using ProjectiveGeometricAlgebra2d: meet_lines, join_points
using ProjectiveGeometricAlgebra2d: randn_pga

function test_associative(op, x, y, z)
    lhs = op(op(x,y), z)
    rhs = op(x, op(y, z))
    @test typeof(lhs)  == typeof(rhs)
    if !(lhs ≈ rhs)
        @show op
        for pname in propertynames(lhs)
            val_lhs = getproperty(lhs, pname)
            val_rhs = getproperty(rhs, pname)
            if !(val_lhs ≈ val_rhs)
                msg = """
                lhs.$(pname) = $(val_lhs)
                rhs.$(pname) = $(val_rhs)
                """
                println(msg)
            end
        end
    end
    @test lhs ≈ rhs
end

function test_associative(op)

    x = randn_pga()
    y = randn_pga()
    z = randn_pga()
    test_associative(op,x,y,z)
end

@testset "conversion" begin
    @testset "point" begin
        pt_euclid = randn(2)
        pt_pga = pga_point(pt_euclid)
        @test euclidean_point(pt_pga) == pt_euclid
    end
    @testset "line" begin
        l_euclid = NamedTuple{(:a,:b,:c)}(randn(3))
        l_pga = pga_line(l_euclid)
        @test euclidean_line_abc(l_pga) == l_euclid
    end
end

@testset "join and meet" begin

    @test join_points([0,1], [1,0]) === (a=1,b=1,c=-1)
    @test join_points([0,0], [1,0]) === (a=0,b=1,c=0) # TODO correct sign?
    @test join_points([0,1], [0,0]) === (a=1,b=0,c=0) # TODO correct sign?

    l1 = (a=1, b=0, c=-2)
    l2 = (a=0, b=1, c=-1)
    @test meet_lines(l1, l2) == [2,1]
end

@testset "mwe geomul" begin
    @test geomul(PGA2d(e1=1), PGA2d(e1=1)) == PGA2d(e=1)

end

@testset "mwe associativity geomul" begin
    x = PGA2d(e2=1)
    y = PGA2d(e12=1, e012=1)
    z = PGA2d( e02=1)
    xy = geomul(x,y)
    @test xy == PGA2d(e1=-1, e01=1)
    xy_z = geomul(xy, z)
    @test xy_z == PGA2d(e012=1)

    yz = geomul(y,z)
    @test yz == PGA2d(e01=1)
    x_yz = geomul(x, yz)
    @test x_yz == PGA2d(e012=1)
    @test xy_z == x_yz
end

@testset "mwe associtivity geomul" begin
    x = PGA2d(e01=1)
    y = PGA2d(e1=1)
    z = PGA2d(e1=1)

    xy = geomul(x,y)
    @test xy == PGA2d(e0=1)
    xy_z = geomul(xy, z)
    @test xy_z == PGA2d(e01=1)
    yz = geomul(y,z)
    @test yz == PGA2d(e=1)
    x_yz = geomul(x, yz)
    @test x_yz == PGA2d(e01=1)
    @test x_yz == xy_z
end

@testset "mwe associativity wedge" begin
    x = PGA2d(e=1, e1=1, e12=1)
    y = PGA2d(e=1, e0=1, e1=1, e2=1, e01=1)
    z = PGA2d(e0=1, e1=1, e02=1, e01=1)
    @test wedge(wedge(x,y), z) == wedge(x, wedge(y,z))
end

@testset "associativity" begin

    for _ in 1:100
        test_associative(wedge)
        test_associative(anti_wedge)
        test_associative(geomul)
    end
end


text:
using Plots

f(x) = sin(x)

h = 1e-4
df(x) = (f(x+h) - f(x)) / h

x = -pi:0.01:pi
plot(x, f.(x), label="f(x)", xlabel="x", ylabel="y")
plot!(x, df.(x), label="df(x)", xlabel="x", ylabel="y")


text:
const BlockId = Int
const BranchId = Int

# "Mapping from variable to set of variables it corresponds to in inverse"
# d[(v1, id) = Set(v2)] means that v1 in ir corresponds to v2 on the invblock id
const VarMap = Dict{Tuple{Variable, BlockId}, Set{Variable}}

# A Branch on a Block -- need this because branches on different blocks share same id
const BranchBlock = NamedTuple{(:branch, :block), Tuple{BranchId, BlockId}}

"add `v` to `vm[k]`"
function add!(vm::VarMap, v::Variable, bid::BlockId, v2::Variable)
  k = (v, bid)
  if k in keys(vm)
    push!(vm[k], v2)
  else
    vm[k] = Set([v2])
  end
end
@post [@cap(vm)[k] ; v] == vm[k]

unwrap(x::Type) = x
unwrap(x::Mjolnir.Const{T}) where {T} = T

"Mapping from variables (including arguments) to inferred types"
function vartypes(ir::IR)
  vtypes = Dict{Variable, Type}()
  for b in blocks(ir)
    for (v, t) in zip(arguments(b), unwrap.(argtypes(b)))
      vtypes[v] = t
    end
  end
  for v in keys(ir)
    vtypes[v] = ir[v].type
  end
  vtypes
end

"Argument types of statement `stmt` according to vartypes `vtypes`"
function stmtargtypes(stmt::Statement, vtypes)
  t = map(stmt.expr.args[2:end]) do arg
    if isvar(arg)
      vtypes[arg]
    else
      typeof(arg)
    end
  end
  Tuple{t...}
end

#zt - move to helpers
"What branchblocks are incoming into block b"
function incomingbranches(b)
  bs = BranchBlock[]
  for pa in IRTools.predecessors(b)
    for (brid, br) in enumerate(branches(pa))
      if br.block == b.id
        push!(bs, (branch = brid, block = pa.id))
      end
    end
  end
  bs
end

"Arguments of a statement"
fargs(b::Block) =
  length(arguments(b)) > 1 ? arguments(b)[2:end] : []

"Head of expression defined by statement"
head(stmt::Statement) =  stmt.expr.args[1]

"Each statement of block `b` and associated Variable"
varstatements(b::Block) = collect((var = k, stmt = b[k]) for k in keys(b))


text:
include("./exportData.jl")
using .exportData, DifferentialEquations, DiffEqParamEstim, PyPlot, Optim, DelimitedFiles

# define fenkar system
function H(x;k=100.0)
	return 0.5*(1.0 + tanh(k*x))
end

function fenkar!(du, u, p, t)

	# parameters
	@views tsi,tv1m,tv2m,tvp,twm,twp,td,to,tr,xk,uc,uv,ucsi,t0,TI,IA = p[1:16]

	# fenkar dynamics
	du[1] = H(t-t0;k=1.0)*IA*sin(pi*(t-t0)/TI)^500 - (u[1]*H(uc-u[1])/to + H(u[1]-uc)/tr - u[2]*H(u[1]-uc)*(1.0-u[1])*(u[1]-uc)/td - u[3]*H(u[1]-ucsi;k=xk)/tsi)
	du[2] = H(uc-u[1])*(1.0-u[2])/(tv1m*H(uv-u[1]) + tv2m*H(u[1]-uv)) - H(u[1]-uc)*u[2]/tvp
	du[3] = H(uc-u[1])*(1.0-u[3])/twm - H(u[1]-uc)*u[3]/twp

	return nothing
end

# fenkar system parameters in p[1:13] (BR; from https://doi.org/10.1063/1.166311)
# and then p[14:16] = [offset_time, stimulus_period, stimulus_amplitude]
p = [29.0,19.6,1250.0,3.33,41.0,870.0,0.25,12.5,33.3,10.0,0.13,0.04,0.85,0.0,125.0,0.05]
tspan = (0.0, 3000.0)
u0 = zeros(Float64,3)

#solve and take data
prob = ODEProblem(fenkar!, u0, tspan, p)
sol = solve(prob, Tsit5(); saveat=tspan[begin]:2.0:tspan[end])

# getExpData
ind = 114
t, data, p[15] = getExpData(ind, tInds=1:1000)
p[14] = t[findfirst(data[1,:] .> 0.5)]-p[15]/2

# remake problem
prob = remake(prob, tspan=(0.0, t[end]), p=p)

# use the save_idxs to select the mid-point of the u-variable
sol = solve(prob, Tsit5(); saveat=t, save_idxs=1:1)

fig = plt.figure(figsize=(4,2))
plt.plot(t, data[:], "-", label="Exp. Sample")
plt.plot(sol.t, sol[1,:], "--", label="Initial Parameters")
plt.xlabel("\$ t \$")

# define the cost function
cost_function = build_loss_objective(prob, Tsit5(), L2Loss(t,data), maxiters=10000, verbose_opt=true, verbose_steps=50; save_idxs=1:1);

# define limits for parameters
#          tsi,  tv1m,   tv2m,   tvp,   twm,    twp,  td,    to,    tr,   xk,  uc,    uv, ucsi,        t0,     TI,   IA
#   p = [ 29.0,  19.6, 1250.0,   3.33, 41.0,  870.0, 0.25, 12.5,  33.3, 10.0, 0.13, 0.04,  0.85,      1.0 ,   BCL, 0.05]
lower = [  1.0,  10.0,  400.0,   2.33, 10.0,  100.0, 0.10, 10.0,  10.0,  9.0, 0.10, 0.01,  0.10,      0.0,   40.0, 0.00]
upper = [ 50.0,  50.0, 2000.0,   4.33, 50.0, 1000.0, 0.50, 50.0,  50.0, 11.0, 0.90, 0.90,  0.90,     t[1], 1000.0, 1.00]

@assert all(p .>= lower) && all(p .<= upper)

# fit
result = optimize(cost_function, lower, upper, p, Fminbox(NelderMead()))
@info "Initial p = $(p) \n"
@info "Minimizer = $(result.minimizer) \n"
prob = remake(prob, p=result.minimizer)
sol = solve(prob, Tsit5(); saveat=t, save_idxs=1:1)
plt.plot(sol.t, sol[1,:], "--", label="Nelder-Mead Parameters")

result = optimize(cost_function, lower, upper, result.minimizer, Fminbox(BFGS()))
@info "Initial p = $(p) \n"
@info "Minimizer = $(result.minimizer) \n"
prob = remake(prob, p=result.minimizer)
sol = solve(prob, Tsit5(); saveat=t, save_idxs=1:1)
plt.plot(sol.t, sol[1,:], "--", label="BFGS Parameters")

plt.savefig("./ode/optim/$(ind)_fittings.pdf")
plt.close()

open("./ode/optim/$(ind)_params.txt", "w") do io
	writedlm(io, ["tsi" "tv1m" "tv2m" "tvp" "twm" "twp" "td" "to" "tr" "xk" "uc" "uv" "ucsi" "t0" "TI" "IA"])
	writedlm(io, "\n")
	writedlm(io, transpose(p))
	writedlm(io, "\n")
	writedlm(io, transpose(result.minimizer))
end


text:
# # Section Allocation
# Suppose you have $n$ students in a class who need to be assigned to $m$ discussion sections. Each student needs to be assigned to exactly one section. Each discussion section should have between 6 and 10 students. Suppose an $n \times m$ preference matrix $P$ is given, where $P_{ij}$ gives student $i$'s ranking for section $j$ (1 would mean it is the student's top choice, 10,000 or a large number would mean the student can not attend that section).
#
# The goal will be to get an allocation matrix $X$, where $X_{ij} = 1$ if student $i$ is assigned to section $j$ and $0$ otherwise.

using Convex, GLPKMathProgInterface
aux(str) = joinpath(@__DIR__, "aux", str) # path to auxiliary files

# Load our preference matrix, P
include(aux("data.jl"))
summary(P)

X = Variable(size(P), :Bin)

# We want every student to be assigned to exactly one section. So, every row must have exactly one non-zero entry
# In other words, the sum of all the columns for every row is 1
# We also want each section to have between 6 and 10 students, so the sum of all the rows for every column should
# be between these
constraints = [sum(X, dims=2) == 1, sum(X, dims=1) <= 10, sum(X, dims=1) >= 6]

# Our objective is simple sum(X .* P), which can be more efficiently represented as vec(X)' * vec(P)
# Since each entry of X is either 0 or 1, this is basically summing up the rankings of students that were assigned to them.
# If all students got their first choice, this value will be the number of students since the ranking of the first choice is 1.
p = minimize(vec(X)' * vec(P), constraints)

solve!(p, GLPKSolverMIP())
p.optval


text:
function write_netcdf2d_reg(fname, var, lon, lat, varname, filenc)
    (nsx,) = size(lon)
    (nsy,) = size(lat)
    (ns, ns, nt) = size(var)

    mode = NC_NETCDF4

    ncin = NetCDF.open(filenc)
    lonatt = ncin.vars["lon"].atts
    latatt = ncin.vars["lat"].atts
    timeatt = ncin.vars["time"].atts
    varatt = ncin.vars[varname].atts
    tim = NetCDF.readvar(ncin,"time")

    #ncgetatt(filenc, "global", gatts)

    varattr = Dict()
    for k in collect(keys(varatt))
        if((k!="_FillValue")&&(k!="add_offset")&&(k!="scale_factor"))
            varattr[k] = varatt[k]
        end
    end

    isfile(fname) ? rm(fname) : nothing
    nccreate(fname, varname , "lon" , nsx, "lat", nsy, "time", tim[1:nt],
             timeatt, atts=varattr,mode=mode,t=NC_FLOAT)
    nccreate(fname, "lon" , "lon" , nsx,  atts=lonatt, mode=mode, t=NC_FLOAT)
    nccreate(fname, "lat" , "lat" , nsy,  atts=latatt, mode=mode, t=NC_FLOAT)

    ncwrite(var, fname, varname)
    ncwrite(lon, fname, "lon")
    ncwrite(lat, fname, "lat")
end


text:
abstract type Particle end;
abstract type Grid end;

mutable struct Particle2D <: Particle
    x::Vector{Float32};
    v::Vector{Float32};
    C::Matrix{Float32};
    V::Float32;
    J::Float32;
    property::Property2D;
end

mutable struct Grid2D <: Grid
    size_x::Int32;
    size_y::Int32;
    dh::Float32;
    field_v::Matrix{Vector{Float32}};
    field_m::Matrix{Float32};
end

function Particle2D(x0::Vector{Float32}, v0::Vector{Float32}, V0::Float32, p::Property2D)
    C = zeros(2, 2);
    return Particle2D(x0, v0, C, V0, 1.0, p);
end

function getMass(self::Particle2D)
    return self.property.M.rho * self.V;
end

function Grid2D(size_x::Int64, size_y::Int64, dh::Float32)
    filed_v = Matrix{Vector{Float32}}(undef, size_x, size_y);
    filed_m = Matrix{Float32}(undef, size_x, size_y);
    return Grid2D(size_x, size_y, dh, filed_v, filed_m);
end

function initGrid(self::Grid2D)
    for i = 1 : self.size_x
        for j = 1 : self.size_y
            self.field_m[i, j] = 0.0;
            self.field_v[i, j] = [0.0, 0.0];
        end
    end
end


text:
# A key principle in building up the problem structure is that there should be no mutable
# state that gets modified during the continuation run. This means that things like the size
# of (variable-size) variables is not included in the problem structure; fixed-size
# variables can be determined using initial_u.

# The key benefit of this principle is that it becomes possible to reuse a problem structure
# multiple times and, possibly more importantly, it becomes impossible to corrupt the
# problem structure during construction.

export Var, Data, Func, Problem
export glue, get_problem, get_flatproblem

"""
    $(TYPEDEF)

An abstract representation of a continuation variable, that is, state that is continually
updated during continuation.
"""
struct Var
    name::String
    initial_dim::Int64
    initial_u::Any
    initial_t::Any
    top_level::Bool
end
function Var(
    name;
    initial_u=nothing,
    initial_t=nothing,
    initial_dim=length(initial_u),
    top_level=false,
)
    return Var(name, initial_dim, initial_u, initial_t, top_level)
end

function Base.show(io::IO, mime::MIME"text/plain", var::Var)
    println(io, "$Var(\"$(var.name)\")")
    println(io, "    initial_dim → $(var.initial_dim)")
    println(io, "    initial_u → $(var.initial_u)")
    return print(io, "    initial_t → $(var.initial_t)")
end

const ALLVARS = Var("ALLVARS"; initial_dim=0, top_level=true)

"""
    $(TYPEDEF)

An abstract representation of continuation data, that is, information that may change
between two continuation steps but not during a single continuation step.
"""
struct Data
    name::String
    initial_data::Any
end

function Base.show(io::IO, mime::MIME"text/plain", data::Data)
    println(io, "$Data(\"$(data.name)\")")
    return print(io, "    initial_data → $(data.initial_data)")
end

"""
    $(TYPEDEF)

An abstract representation of a function of the form `f!(output, var, data)`.
"""
struct Func
    name::String
    func::Any
    initial_dim::Int64
    initial_f::Any
    group::Vector{Symbol}
    pass_problem::Bool
    var::Vector{Var}
    var_names::Dict{String,Int64}
    data::Vector{Data}
    data_names::Dict{String,Int64}
end

function Func(
    name,
    func,
    var=(),
    data=();
    initial_f=nothing,
    initial_dim=length(initial_f),
    group=[:embedded],
    pass_problem=false,
)
    f = Func(
        name,
        func,
        initial_dim,
        initial_f,
        group,
        pass_problem,
        Var[],
        Dict{String,Int64}(),
        Data[],
        Dict{String,Int64}(),
    )
    for v in var
        push!(f, v)
    end
    for d in data
        push!(f, d)
    end
    return f
end

function Base.show(io::IO, mime::MIME"text/plain", func::Func)
    println(io, "$Func(\"$(func.name)\")")
    println(io, "    initial_dim → $(func.initial_dim)")
    println(io, "    initial_f → $(func.initial_f)")
    println(io, "    group → $(func.group)")
    println(io, "    pass_problem → $(func.pass_problem)")
    println(io, "    var → $([nameof(v) for v in func.var])")
    return print(io, "    data → $([nameof(d) for d in func.data])")
end

function glue(var1::Var, var2::Var)
    if var1.initial_dim != var2.initial_dim
        throw(ArgumentError("Size mismatch between var1 and var2"))
    end
    return Func(
        var1.name * "=" * var2.name,
        (res, u) -> res .= u[1] .- u[2],
        (var1, var2);
        initial_dim=var1.initial_dim,
    )
end

"""
Return the underlying problem tree structure.
"""
function get_problem end

"""
    $(TYPEDEF)

An abstract type for `Problem` owners to inherit from. (Inheriting is not required but
reduces boilerplate.)
"""
abstract type ProblemOwner end

(owner::ProblemOwner)(::Signal, args...) = nothing  # fallback for signal handling

"""
    $(TYPEDEF)

An abstract representation of a continuation problem.
"""
struct Problem
    name::String
    func::Vector{Func}
    func_names::Dict{String,Int64}
    problem::Vector{Problem}
    problem_names::Dict{String,Int64}
    owner::Any
    owner_pass_indices::Vector{Union{Var,Data,Func,Problem}}
end

const ProblemTypes = Union{Var,Data,Func,Problem}  # Cannot put this earlier because of type recursion issue

function Problem(name, owner=nothing, func=(), problem=())
    prob = Problem(
        name,
        Vector{Func}[],
        Dict{String,Int64}(),
        Problem[],
        Dict{String,Int64}(),
        owner,
        ProblemTypes[],
    )
    for f in func
        push!(prob, f)
    end
    for p in problem
        push!(prob, p)
    end
    return prob
end

function Base.show(io::IO, mime::MIME"text/plain", problem::Problem)
    println(io, "$Problem(\"$(problem.name)\")")
    println(io, "    func → $([nameof(f) for f in problem.func])")
    return print(io, "    problem → $([nameof(p) for p in problem.problem])")
end

get_problem(problem::Problem) = problem

"""
    $(SIGNATURES)

Request that the indices of a particular item (Var, Data, Func, or Problem) are passed to
the problem owner when signalled.
"""
function pass_indices!(problem::Problem, item::ProblemTypes)
    push!(problem.owner_pass_indices, item)
    return problem
end

# Generic functions for all collections
for (Collection, Item, name) in (
    (Func, Var, :var),
    (Func, Data, :data),
    (Problem, Func, :func),
    (Problem, Problem, :problem),
)
    @eval function Base.push!(collection::$Collection, item::$Item)
        if item in collection.$name
            throw(ArgumentError(string($Item) * " already added: " * item.name))
        else
            push!(collection.$name, item)
            if !isempty(item.name)
                idx = lastindex(collection.$name)
                if haskey(collection.$(Symbol(name, "_names")), item.name)
                    @warn "Duplicate name" item.name item
                else
                    collection.$(Symbol(name, "_names"))[item.name] = idx
                end
            end
        end
        return collection
    end

    @eval function Base.append!(
        collection::$Collection, items::Union{NTuple{<:Any,$Item},AbstractVector{$Item}}
    )
        for item in items
            push!(collection, item)
        end
        return collection
    end

    @eval function Base.getindex(
        collection::$Collection, ::Type{$Item}, item::AbstractString
    )
        return collection.$name[collection.$(Symbol(name, "_names"))[item]]
    end
end

for (Collection, names) in ((Func, (:var, :data)), (Problem, (:func, :problem)))
    @eval function Base.getindex(collection::$Collection, item::AbstractString)
        itempath = split(item, "."; limit=2)
        for subcollection in $names
            names_dict = getfield(collection, Symbol(subcollection, "_names"))
            if haskey(names_dict, itempath[1])
                idx = names_dict[itempath[1]]
                if length(itempath) == 1
                    return getfield(collection, subcollection)[idx]
                else
                    return getfield(collection, subcollection)[idx][itempath[2]]
                end
            end
        end
        return throw(KeyError(item))
    end
end

for Item in (Var, Data, Func, Problem)
    @eval Base.nameof(item::$Item) = item.name
    @eval ownerof(item::$Item) = item.owner
end

"""
Return the underlying flattened problem structure (includes generated code).
"""
function get_flatproblem end

"""
    $(TYPEDEF)

A flattened representation of a continuation problem.
"""
struct FlatProblem{G,O}
    var::Vector{Var}
    var_names::Dict{String,Int64}
    data::Vector{Data}
    data_names::Dict{String,Int64}
    func::Vector{Func}
    func_names::Dict{String,Int64}
    group::Vector{Vector{Func}}
    group_names::Dict{Symbol,Int64}
    problem::Vector{Problem}
    problem_names::Dict{String,Int64}
    call_group::G
    call_owner::O
end

function FlatProblem()
    return FlatProblem(
        Var[],
        Dict{String,Int64}(),
        Data[],
        Dict{String,Int64}(),
        Func[],
        Dict{String,Int64}(),
        Vector{Func}[],
        Dict{Symbol,Int64}(),
        Problem[],
        Dict{String,Int64}(),
        nothing,
        nothing,
    )
end

function Base.show(io::IO, mime::MIME"text/plain", flat::FlatProblem)
    println(io, "$FlatProblem()")
    println(io, "    var → $(collect(keys(flat.var_names)))")
    println(io, "    data → $(collect(keys(flat.data_names)))")
    println(io, "    func → $(collect(keys(flat.func_names)))")
    println(io, "    problem → $(collect(keys(flat.problem_names)))")
    return print(io, "    group → $(collect(keys(flat.group_names)))")
end

function evaluate!(res, flat::FlatProblem, group::Symbol, u, data, problem)
    return flat.call_group[group](res, u, data, problem)
end

function signal!(flat::FlatProblem, signal::Signal, args...)
    return flat.call_owner(signal, args...)
end

function signal!(signal::Signal, flat::FlatProblem, args...)
    return flat.call_owner(signal, flat, args...)
end

get_problem(flat::FlatProblem) = flat.problem[1]
get_flatproblem(flat::FlatProblem) = flat

"""
    $(SIGNATURES)

Take a hierarchical `Problem` structure and return a flattened version that contains
generated code for each function group.
"""
function flatten(problem::Problem)
    flat = FlatProblem()
    _flatten!(flat, problem, "")
    call_group = NamedTuple{(keys(flat.group_names)...,)}((
        eval(_gen_call_group(flat, i)) for i in eachindex(flat.group)
    ))
    call_owner = eval(_gen_call_owner(flat))
    return FlatProblem(
        flat.var,
        flat.var_names,
        flat.data,
        flat.data_names,
        flat.func,
        flat.func_names,
        flat.group,
        flat.group_names,
        flat.problem,
        flat.problem_names,
        call_group,
        call_owner,
    )
end

const NAME_SEP = "."

function _flatten!(flat::FlatProblem, problem::Problem, basename)
    if !(problem in flat.problem)
        push!(flat.problem, problem)
        pidx = lastindex(flat.problem)
        basename *= problem.name
        haskey(flat.problem_names, basename) &&
            @warn "Duplicate Problem name" basename problem
        flat.problem_names[basename] = pidx
        # Iterate over the sub-problems (depth first)
        for subproblem in problem.problem
            _flatten!(flat, subproblem, basename * NAME_SEP)
        end
        # Iterate over functions in the problem
        for func in problem.func
            _flatten!(flat, func, basename * NAME_SEP)
        end
    else
        throw(ErrorException("Duplicate Problem within the problem structure"))
    end
    return flat
end

function _flatten!(flat::FlatProblem, func::Func, basename)
    if !(func in flat.func)
        # Func has not been previously added, so add it
        push!(flat.func, func)
        fidx = lastindex(flat.func)
        basename *= func.name
        haskey(flat.func_names, basename) && @warn "Duplicate Func name" basename func
        flat.func_names[basename] = fidx
        # Iterate over the groups the Func belongs to
        for group in func.group
            if !haskey(flat.group_names, group)
                push!(flat.group, Func[])
                gidx = lastindex(flat.group)
                flat.group_names[group] = gidx
            else
                gidx = flat.group_names[group]
            end
            push!(flat.group[gidx], func)
        end
        # Iterate over the Vars belonging to the Func
        for var in func.var
            _flatten!(flat, var, basename * NAME_SEP)
        end
        # Iterate over the Data belonging to the Func
        for data in func.data
            _flatten!(flat, data, basename * NAME_SEP)
        end
    else
        throw(ErrorException("Duplicate Func within the problem structure"))
    end
    return flat
end

function _flatten!(flat::FlatProblem, var::Var, basename)
    if !(var in flat.var)
        # Not been previously added, so add it
        push!(flat.var, var)
        idx = lastindex(flat.var)
        basename = var.top_level ? var.name : basename * var.name
        haskey(flat.var_names, basename) && @warn "Duplicate Var name" basename var
        flat.var_names[basename] = idx
    end
    return flat
end

function _flatten!(flat::FlatProblem, data::Data, basename)
    if !(data in flat.data)
        # Not been previously added, so add it
        push!(flat.data, data)
        idx = lastindex(flat.data)
        basename *= data.name
        haskey(flat.data_names, basename) && @warn "Duplicate Data name" basename data
        flat.data_names[basename] = idx
    end
    return flat
end

function _gen_call_group(flat::FlatProblem, group::Integer)
    func = :(function (res, u, data, problem) end)
    for (i, f) in enumerate(flat.group[group])
        func_f = :($(f.func)(res[$i]))
        u = :(())
        for var in f.var
            if var === ALLVARS
                push!(u.args, :(u[:]))
            else
                idx = findfirst(==(var), flat.var)
                push!(u.args, :(u[$idx]))
            end
        end
        if length(u.args) == 0
            push!(func_f.args, :(u[:]))
        elseif length(u.args) == 1
            push!(func_f.args, u.args[1])
        else
            push!(func_f.args, u)
        end
        d = :(())
        for data in f.data
            idx = findfirst(==(data), flat.data)
            push!(d.args, :(data[$idx]))
        end
        if length(d.args) == 1
            push!(func_f.args, d.args[1])
        elseif length(d.args) > 1
            push!(func_f.args, d)
        end
        if f.pass_problem
            push!(func_f.args, :problem)
        end
        push!(func.args[2].args, func_f)
    end
    push!(func.args[2].args, :nothing)
    return func
end
function _gen_call_group(flat::FlatProblem, group::Symbol)
    return _gen_call_group(flat, flat.group_names[group])
end

function _gen_call_owner(flat::FlatProblem)
    func = :(function (signal::Signal, args...) end)
    for problem in flat.problem
        if problem.owner !== nothing
            indices = Int64[]
            for item in problem.owner_pass_indices
                if item isa Var
                    idx = findfirst(==(item), flat.var)
                elseif item isa Data
                    idx = findfirst(==(item), flat.data)
                elseif item isa Func
                    idx = findfirst(==(item), flat.func)
                elseif item isa Problem
                    idx = findfirst(==(item), flat.problem)
                else
                    idx = nothing
                end
                if idx === nothing
                    throw(ErrorException("Requested item does not exist in the problem structure: $item"))
                else
                    push!(indices, idx)
                end
            end
            # Using the call signature owner(signal, problem) means that we don't have to
            # worry about whether we are adding methods to the correct function. For
            # example, we could have two covering algorithms that provide an adapt signal -
            # which function do you overload?
            if isempty(indices)
                push!(func.args[2].args, :($(problem.owner)(signal, args...)))
            else
                push!(
                    func.args[2].args,
                    :($(problem.owner)(signal, $((indices...,)), args...)),
                )
            end
        end
    end
    push!(func.args[2].args, :nothing)
    return func
end


text:
abstract type AbstractFVSolution{T,N} end

struct FVSolution{T,N,mType,M} <: AbstractFVSolution{T,N}
  ode_sol::AbstractODESolution{T,N}
  mesh::mType
end

function FVSolution(ode_sol::AbstractODESolution{T,N}, mesh::M, vars::Int) where {T,N,M}
  FVSolution{T,N,M,vars}(ode_sol, mesh)
end

function Base.show(io::IO, A::FVSolution)
  show(io,A.ode_sol)
  show(io, A.mesh)
end

getvalues(sol::FVSolution) = sol.ode_sol.u
gettimes(sol::FVSolution) = sol.ode_sol.t
getmesh(sol::FVSolution) = sol.mesh
getnvars(sol::FVSolution{T,N,m,M}) where {T,N,m,M} = M

function save_csv(sol::FVSolution, file_name::String; idx = -1)
  if !endswith(file_name,".csv")
    file_name = "$file_name.csv"
  end
  if idx == -1
    writedlm(file_name, hcat(cell_centers(getmesh(sol)),getvalues(sol)[end]), ',')
  else
    writedlm(file_name, hcat(cell_centers(getmesh(sol)),getvalues(sol)[idx]), ',')
  end
end

function get_total_u(sol::FVSolution{T}) where {T}
    masa = fill(zero(T),size(getvalues(sol),1))
    for (i,u) in enumerate(getvalues(sol))
        masa[i] = sum(u*cell_volume(getmesh(sol),i))
    end
    return masa
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    dfstree, dfstree_edges, cotree_edges, circuitrank

function dfstree(adjfunc, graph, root)
    stack = [root]
    pred = Dict{Int,Int}(root => 0)
    while !isempty(stack)
        i = pop!(stack)
        for adj in adjfunc(graph, i)
            if !haskey(pred, adj)
                pred[adj] = i
                push!(stack, adj)
            end
        end
    end
    delete!(pred, root)
    return pred
end

function dfstree_edges(nbrfunc, graph, root)
    stack = [root]
    edges = Dict{Int,Int}(root => 0)
    while !isempty(stack)
        i = pop!(stack)
        for (inc, adj) in nbrfunc(graph, i)
            if !haskey(edges, adj)
                edges[adj] = inc
                push!(stack, adj)
            end
        end
    end
    delete!(edges, root)
    return edges
end

"""
    cotree_edges(graph::UndirectedGraph, root) -> Set{Int}

Return a set of co-tree edges of a spanning tree.
"""
function cotree_edges(graph::UndirectedGraph, root)
    edges = dfstree_edges(neighbors, graph, root)
    return setdiff(edgeset(graph), values(edges))
end

"""
    circuitrank(graph::UndirectedGraph) -> Int

Return circuit rank (the number of cycles) of the graph.
"""
function circuitrank(graph::UndirectedGraph)
    nodecount(graph) == 0 && return 0
    if isdefined(graph, :cache) && haskey(graph.cache, :mincycles)
        return length(graph.cache[:mincycles])
    end
    treesize = 0
    for conn in connected_components(graph)
        root = iterate(conn)[1]
        edges = dfstree_edges(neighbors, graph, root)
        treesize += length(edges)
    end
    return edgecount(graph) - treesize
end


text:
__precompile__()
module WRS

using Distributions
using Distributed

export hd, pb2gen

# This implements (3.16) from (2017) p. 71-72
function hd(X, q=0.5, issorted=false)
    n = length(X)

    a = (n+1)*q
    b = (n+1)*(1-q)

    β = Beta(a, b)

    sorted_X = issorted ? X : sort(X)

    θ = 0.
    for i=1:n
        w = cdf(β, i/n) - cdf(β, (i-1)/n)
        θ += w*sorted_X[i]
    end

    return θ
end

function pb2gen(x, y, est, q=0.5; alpha=0.05, nboot=2000)
    bootstrapped_diff_est = pmap(1:nboot) do i
        sampled_values_x = sample(x, length(x), replace=true)
        sampled_values_y = sample(y, length(y), replace=true)

        sort!(sampled_values_x)
        sort!(sampled_values_y)

        return est(sampled_values_x, q, true) - est(sampled_values_y, q, true)
    end
    sort!(bootstrapped_diff_est)

    low = round(Int, (alpha/2)*nboot)+1
    up = nboot-low

    est_x = est(x, q)
    est_y = est(y, q)
    est_diff = est_x - est_y
    ci = (bootstrapped_diff_est[low], bootstrapped_diff_est[up])

    A = count(i->i<0,bootstrapped_diff_est)
    C = count(i->i==0,bootstrapped_diff_est)
    p_hat_star = A/nboot+0.5*C/nboot
    pvalue = 2*(min(p_hat_star,1-p_hat_star))

    se = var(bootstrapped_diff_est)

    return est_x, est_y, est_diff, ci, pvalue, se
end

end # module


text:
using ControlBenchmarks
using ControlSystems
using SparseArrays

benchmarkProb = controlbenchmark( B767FlutterCondition() )

@test ControlSystems.nstates( benchmarkProb.sys ) == 55
@test ControlSystems.ninputs( benchmarkProb.sys ) == 2

# Ensure it is a sparse matrix
@test issparse( benchmarkProb.sys.A )
@test issparse( benchmarkProb.sys.B )
@test issparse( benchmarkProb.sys.C )
@test issparse( benchmarkProb.sys.D )


text:
function convex_hull(df::AbstractArray)
    T, N = size(df, 1), size(df, 3)
    result = Matrix{Union{Float64,Nothing}}(
        nothing, T, N
    )
    @inbounds for n = 1:N
        @views m = convert.(Float64, remove_nothing(df[:, :, n]))
        # m = extract(df, n, id, [x, y])
        add_noise!(m)
        T = size(m, 1)
        # if T <= 90
        #     continue
        # end
        for t = 3:T-2
            result[t-2, n] = chull(m[t-2:t+2, :]).area
        end
    end
    result
end

function convex_hull(
    df::DataFrame;
    id = :TrackID,
    x = :POSITION_X,
    y = :POSITION_Y,
    frame = :FRAME2
)
    convex_hull(xy2matrix(df, id, x, y, frame))
end

max_stars_repo_name: SosUts/spt

max_issues_repo_name: SosUts/spt

max_forks_repo_name: SosUts/spt


text:
using LazySets.Arrays: SingleEntryVector

@testset "LGG09 algorithm: template construction" begin
    # case where the template directions is as single vector
    dirs = fill(0.0, 5); dirs[1] = 1.0 # e1
    alg = LGG09(δ=0.01, template=dirs)
    alg.template == CustomDirections([dirs])

    # case where the template directions is a vector of vectors
    dirs = [SingleEntryVector(1, 5, 1.0), SingleEntryVector(1, 5, -1.0)]
    alg = LGG09(δ=0.01, template=dirs)
    alg.template == CustomDirections(dirs)

    # alias dirs
    dirs = [SingleEntryVector(1, 5, 1.0), SingleEntryVector(1, 5, -1.0)]
    alg = LGG09(δ=0.01, dirs=dirs)
    alg.template == CustomDirections(dirs)

    # specification of vars (requires n)
    @test_throws ArgumentError LGG09(δ=0.01, vars=(1))
    alg = LGG09(δ=0.01, vars=(1), n=5)
    @test collect(alg.template) == dirs

    alg = LGG09(δ=0.01, vars=(1, 3), n=5)
    dirs_1_3 = [SingleEntryVector(1, 5, 1.0), SingleEntryVector(3, 5, 1.0),
                SingleEntryVector(1, 5, -1.0), SingleEntryVector(3, 5, -1.0)]
    @test collect(alg.template) == dirs_1_3

    # alias using dim
    alg2 = LGG09(δ=0.01, vars=(1, 3), dim=5)
    @test collect(alg2.template) == dirs_1_3
end

@testset "LGG09 algorithm: 1d" begin
    # one-dimensional
    prob, dt = exponential_1d()
    box1d = CustomDirections([[1.0], [-1.0]])
    sol = solve(prob, tspan=dt, LGG09(δ=0.01, template=box1d))
    @test isa(sol.alg, LGG09)
    @test setrep(sol) <: HPolyhedron
    @test setrep(sol) == HPolyhedron{Float64, Vector{Float64}}
    @test dim(sol) == 1
end

@testset "LGG09 algorithm: 5d" begin
    # higher-dimensional
    prob, dt = linear5D()
    box5d = BoxDirections{Float64, Vector{Float64}}(5)
    sol = solve(prob, tspan=dt, LGG09(δ=0.01, template=box5d))
    @test setrep(sol) == HPolyhedron{Float64, Vector{Float64}}
    @test dim(sol) == 5

    # equivalent algorithm definitions
    alg0 = LGG09(δ=0.01, template=box5d)
    alg1 = LGG09(δ=0.01, template=BoxDirections(5))
    alg2 = LGG09(δ=0.01, template=:box, n=5)
    alg3 = LGG09(δ=0.01, template=:box, dim=5)
    @test alg1 == alg2 == alg3

    @test !(alg0 == alg1) # alg1-3 use single-entry vectors, but alg0 uses Vector
    @test collect(alg0.template) == Vector.(collect(alg1.template))

    # other directions
    alg4 = LGG09(δ=0.01, template=:oct, dim=5)
    @test alg4.template isa OctDirections && dim(alg4.template) == 5
end

@testset "LGG09 algorithm: underapproximation" begin
    X0 = BallInf([1.0, 1.0], 0.1)
    function foo()
        A = [0.0 1.0; -1.0 0.0]
        prob = @ivp(x' = Ax, x(0) ∈ X0)
        tspan = (0.0, 20.0)
        return prob, tspan
    end
    prob, dt = foo()
    δ = 2e-1

    approx_model = ReachabilityAnalysis.ForwardDdt()
    sol = solve(prob, tspan=(0.0, 20.0),
                alg=LGG09(δ=δ, template=PolarDirections(40), approx_model=approx_model))
    approx_model = ReachabilityAnalysis.ForwardDdt(oa=false)
    sol_ua = solve(prob, tspan=(0.0, 20.0),
                   alg=LGG09(δ=δ, template=PolarDirections(40), approx_model=approx_model))
end


text:
# This file is a part of BATBase.jl, licensed under the MIT License (MIT).

__precompile__(true)

"""
    BATBase

Template for Julia packages.
"""
module BATBase

using Random

using ArraysOfArrays
using Distributions
using DocStringExtensions
using StatsBase
using ValueShapes

include("basetypes/basetypes.jl")
include("algotypes/algotypes.jl")

end # module


text:
# ---
# title: 1005. Maximize Sum Of Array After K Negations
# id: problem1005
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Greedy
# link: <https://leetcode.com/problems/maximize-sum-of-array-after-k-negations/description/>
# hidden: true
# ---
#
# Given an array `A` of integers, we **must**  modify the array in the following
# way: we choose an `i` and replace `A[i]` with `-A[i]`, and we repeat this
# process `K` times in total.  (We may choose the same index `i` multiple
# times.)
#
# Return the largest possible sum of the array after modifying it in this way.
#
#
#
# **Example 1:**
#
#
#
#     Input: A = [4,2,3], K = 1
#     Output: 5
#     Explanation: Choose indices (1,) and A becomes [4,-2,3].
#
#
# **Example 2:**
#
#
#
#     Input: A = [3,-1,0,2], K = 3
#     Output: 6
#     Explanation: Choose indices (1, 2, 2) and A becomes [3,1,0,2].
#
#
# **Example 3:**
#
#
#
#     Input: A = [2,-3,-1,5,-4], K = 2
#     Output: 13
#     Explanation: Choose indices (1, 4) and A becomes [2,3,-1,5,4].
#
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 10000`
#   2. `1 <= K <= 10000`
#   3. `-100 <= A[i] <= 100`
#
#
## @lc code=start
using LeetCode

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


text:
using Test
using BackedUpImmutable

d = StaticDict(["a" => 7, "d" => 2])
@test d["d"] == 2

d2 = StaticDict("a" => 7, "d" => 2, "e" => 4)
@test d2["e"] != 3

function testBackedUpImmutableDict()
    fibr = BackedUpImmutableDict{String,Int64}(["a" => 0, "b" => 1, "c" => 1, "d" => 2,
        "e" => 3, "f" => 5, "g" => 8, "h" => 13, "i" => 21, "j" => 34, "extra" => -1])

    x = fibr["extra"]
    @test x == -1
    fibr["extra"] = 0
    y = fibr["extra"]
    @test y == 0
    restore!(fibr, "extra")
    z = fibr["extra"]
    @test z == -1

    @test_throws String begin fibr["k"] = 55 end

    # test alternative constructor
    fibr = BackedUpImmutableDict{String, Int64}("a" => 0, "b" => 1, "c" => 1, "d" => 2,
        "e" => 3, "f" => 5, "g" => 8, "h" => 13, "i" => 21, "j" => 34, "extra" => -1)

    fibr["a"] = 9
    fibr["b"] = 7

    # test restore all to default
    restoreall!(fibr)

    @test fibr["a"] == 0
end

testBackedUpImmutableDict()


text:
# Julia wrapper for header: /usr/include/scip/pub_benders.h
# Automatically generated using Clang.jl wrap_c

function SCIPbendersComp(elem1, elem2)
    ccall((:SCIPbendersComp, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPbendersCompName(elem1, elem2)
    ccall((:SCIPbendersCompName, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPbendersGetData(benders)
    ccall((:SCIPbendersGetData, libscip), Ptr{SCIP_BENDERSDATA}, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersSetData(benders, bendersdata)
    ccall((:SCIPbendersSetData, libscip), Cvoid, (Ptr{SCIP_BENDERS}, Ptr{SCIP_BENDERSDATA}), benders, bendersdata)
end

function SCIPbendersGetName(benders)
    ccall((:SCIPbendersGetName, libscip), Cstring, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersGetDesc(benders)
    ccall((:SCIPbendersGetDesc, libscip), Cstring, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersGetPriority(benders)
    ccall((:SCIPbendersGetPriority, libscip), Cint, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersGetNSubproblems(benders)
    ccall((:SCIPbendersGetNSubproblems, libscip), Cint, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersSubproblem(benders, probnumber)
    ccall((:SCIPbendersSubproblem, libscip), Ptr{SCIP_}, (Ptr{SCIP_BENDERS}, Cint), benders, probnumber)
end

function SCIPbendersGetNCalls(benders)
    ccall((:SCIPbendersGetNCalls, libscip), Cint, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersGetNCutsFound(benders)
    ccall((:SCIPbendersGetNCutsFound, libscip), Cint, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersGetSetupTime(benders)
    ccall((:SCIPbendersGetSetupTime, libscip), Cdouble, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersGetTime(benders)
    ccall((:SCIPbendersGetTime, libscip), Cdouble, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersIsInitialized(benders)
    ccall((:SCIPbendersIsInitialized, libscip), UInt32, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersIsActive(benders)
    ccall((:SCIPbendersIsActive, libscip), UInt32, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersOnlyCheckConvexRelax(benders)
    ccall((:SCIPbendersOnlyCheckConvexRelax, libscip), UInt32, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersCutLP(benders)
    ccall((:SCIPbendersCutLP, libscip), UInt32, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersCutPseudo(benders)
    ccall((:SCIPbendersCutPseudo, libscip), UInt32, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersCutRelaxation(benders)
    ccall((:SCIPbendersCutRelaxation, libscip), UInt32, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersShareAuxVars(benders)
    ccall((:SCIPbendersShareAuxVars, libscip), UInt32, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersGetAuxiliaryVar(benders, probnumber)
    ccall((:SCIPbendersGetAuxiliaryVar, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_BENDERS}, Cint), benders, probnumber)
end

function SCIPbendersGetAuxiliaryVars(benders)
    ccall((:SCIPbendersGetAuxiliaryVars, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersSetSubproblemObjval(benders, probnumber, objval)
    ccall((:SCIPbendersSetSubproblemObjval, libscip), Cvoid, (Ptr{SCIP_BENDERS}, Cint, Cdouble), benders, probnumber, objval)
end

function SCIPbendersGetSubproblemObjval(benders, probnumber)
    ccall((:SCIPbendersGetSubproblemObjval, libscip), Cdouble, (Ptr{SCIP_BENDERS}, Cint), benders, probnumber)
end

function SCIPfindBenderscut(benders, name)
    ccall((:SCIPfindBenderscut, libscip), Ptr{SCIP_BENDERSCUT}, (Ptr{SCIP_BENDERS}, Cstring), benders, name)
end

function SCIPbendersGetBenderscuts(benders)
    ccall((:SCIPbendersGetBenderscuts, libscip), Ptr{Ptr{SCIP_BENDERSCUT}}, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersGetNBenderscuts(benders)
    ccall((:SCIPbendersGetNBenderscuts, libscip), Cint, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersSetBenderscutPriority(benders, benderscut, priority)
    ccall((:SCIPbendersSetBenderscutPriority, libscip), SCIP_RETCODE, (Ptr{SCIP_BENDERS}, Ptr{SCIP_BENDERSCUT}, Cint), benders, benderscut, priority)
end

function SCIPbendersSetSubproblemIsConvex(benders, probnumber, isconvex)
    ccall((:SCIPbendersSetSubproblemIsConvex, libscip), Cvoid, (Ptr{SCIP_BENDERS}, Cint, UInt32), benders, probnumber, isconvex)
end

function SCIPbendersSubproblemIsConvex(benders, probnumber)
    ccall((:SCIPbendersSubproblemIsConvex, libscip), UInt32, (Ptr{SCIP_BENDERS}, Cint), benders, probnumber)
end

function SCIPbendersGetNConvexSubproblems(benders)
    ccall((:SCIPbendersGetNConvexSubproblems, libscip), Cint, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersSolveSubproblemLP(scip, benders, probnumber, infeasible)
    ccall((:SCIPbendersSolveSubproblemLP, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_BENDERS}, Cint, Ptr{UInt32}), scip, benders, probnumber, infeasible)
end

function SCIPbendersSolveSubproblemCIP(scip, benders, probnumber, infeasible, type, solvecip)
    ccall((:SCIPbendersSolveSubproblemCIP, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_BENDERS}, Cint, Ptr{UInt32}, SCIP_BENDERSENFOTYPE, UInt32), scip, benders, probnumber, infeasible, type, solvecip)
end

function SCIPbendersGetNTransferredCuts(benders)
    ccall((:SCIPbendersGetNTransferredCuts, libscip), Cint, (Ptr{SCIP_BENDERS},), benders)
end

function SCIPbendersUpdateSubproblemLowerbound(benders, probnumber, lowerbound)
    ccall((:SCIPbendersUpdateSubproblemLowerbound, libscip), Cvoid, (Ptr{SCIP_BENDERS}, Cint, Cdouble), benders, probnumber, lowerbound)
end

function SCIPbendersGetSubproblemLowerbound(benders, probnumber)
    ccall((:SCIPbendersGetSubproblemLowerbound, libscip), Cdouble, (Ptr{SCIP_BENDERS}, Cint), benders, probnumber)
end

function SCIPbendersSetSubproblemIsIndependent(benders, probnumber, isindep)
    ccall((:SCIPbendersSetSubproblemIsIndependent, libscip), Cvoid, (Ptr{SCIP_BENDERS}, Cint, UInt32), benders, probnumber, isindep)
end

function SCIPbendersSubproblemIsIndependent(benders, probnumber)
    ccall((:SCIPbendersSubproblemIsIndependent, libscip), UInt32, (Ptr{SCIP_BENDERS}, Cint), benders, probnumber)
end

function SCIPbendersSubproblemIsEnabled(benders, probnumber)
    ccall((:SCIPbendersSubproblemIsEnabled, libscip), UInt32, (Ptr{SCIP_BENDERS}, Cint), benders, probnumber)
end


text:
using SymPy
function Base.show(io::IO, ::MIME"text/latex", x::SymPy.SymbolicObject)
    print(io, SymPy.as_markdown(sympy.latex(x)))
end

# Show SymPy
## Type piracy
function Base.show(io::IO, ::MIME"text/html", x::T) where {T <: SymPy.SymbolicObject}
    #write(io, "
")
    write(io, "
")
    show(io, "text/latex", x)
    write(io, "
")
end

function Base.show(io::IO, ::MIME"text/html", x::Array{T}) where {T <: SymPy.SymbolicObject}
    #write(io, "
")
    write(io, "
")
    show(io, "text/latex", x)
    write(io, "
")
end

## --------

gradient(ex::SymPy.Sym, vars::AbstractArray=free_symbols(ex)) = diff.(ex, [vars...])
divergence(F::Vector{Sym}, vars=free_symbols(F)) = sum(diff.(F, vars))
curl(F::Vector{Sym}, vars=free_symbols(F)) = curl(F.jacobian(vars))


text:
### conversions.jl
#
# Copyright (C) 2016, 2017 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: uncertainty, error propagation, physics
#
# This file is a part of Measurements.jl.
#
# License is MIT "Expat".
#
### Commentary:
#
# This file defines conversion and promotion rules involving Measurement
# objects.
#
### Code:

Base.convert(::Type{Measurement{T}}, a::Irrational) where {T<:AbstractFloat} =
    measurement(T(a))::Measurement{T}
Base.convert(::Type{Measurement{T}}, a::Rational{<:Integer}) where {T<:AbstractFloat} =
    measurement(T(a))::Measurement{T}
Base.convert(::Type{Measurement{T}}, a::Real) where {T<:AbstractFloat} =
    measurement(T(a))::Measurement{T}

Base.convert(::Type{Measurement{T}}, a::Measurement{T}) where {T<:AbstractFloat} = a
function Base.convert(::Type{Measurement{T}},
                      a::Measurement{<:AbstractFloat}) where {T<:AbstractFloat}
    newder = empty_der2(zero(T))
    for tag in keys(a.der)
        newder = Derivatives(newder, (T(tag[1]), T(tag[2]), tag[3])=>T(a.der[tag]))
    end
    return Measurement(T(a.val), T(a.err), a.tag, newder)::Measurement{T}
end

Base.convert(::Type{Measurement}, a::Measurement) = a::Measurement
Base.convert(::Type{Measurement}, a::Rational{<:Integer}) = measurement(a)::Measurement
Base.convert(::Type{Measurement}, a::Real) = measurement(a)::Measurement
# This is used in power of complex numbers.  Allow the conversion only if the uncertainty is
# zero, so there is no loss of information.
function Base.convert(::Type{Int}, a::Measurement)
    @assert isinteger(a)
    return convert(Int, a.val)::Int
end

Base.promote_rule(::Type{Measurement{T}}, ::Type{S}) where {T<:AbstractFloat, S<:Real} =
    Measurement{promote_type(T, S)}
Base.promote_rule(::Type{Measurement{T}},
                  ::Type{Measurement{S}}) where {T<:AbstractFloat, S<:AbstractFloat} =
    Measurement{promote_type(T, S)}


text:
# This file is a part of InverseFunctions.jl, licensed under the MIT License (MIT).
"""
    InverseFunctions

Lightweight package that defines an interface to invert functions.
"""
module InverseFunctions

using Test

include("functions.jl")
include("inverse.jl")
include("test.jl")

end # module


text:
push!(LOAD_PATH, joinpath(@__DIR__, ".."))

using Documenter, DocumenterCitations
using TimeMachine

bib = CitationBibliography(joinpath(@__DIR__, "refs.bib"))

makedocs(
    bib,
    sitename = "TimeMachine",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    modules = [TimeMachine],
    pages = [
        "index.md",
        "algorithms.md",
        "callbacks.md",
        "Background" => [
            "background/lsrk.md",
            "background/ssprk.md",
            "background/ark.md",
            "background/mrrk.md",
        ],
        "references.md",
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/CliMA/TimeMachine.jl.git",
    push_preview = true,
    forcepush = true,
)


text:
using Clang.Generators
using Clang.LibClang.Clang_jll
using Test

include("rewriter.jl")

@testset "Generators" begin
    INCLUDE_DIR = joinpath(Clang_jll.artifact_dir, "include") |> normpath
    CLANG_C_DIR = joinpath(INCLUDE_DIR, "clang-c")

    options = load_options(joinpath(@__DIR__, "test.toml"))

    # add compiler flags
    args = get_default_args()
    push!(args, "-I$INCLUDE_DIR")

    # search top-level headers
    headers = detect_headers(CLANG_C_DIR, args)

    # add extra definition
    @add_def time_t AbstractJuliaSIT JuliaCtime_t Ctime_t

    # create context
    ctx = create_context(headers, args, options)

    # build without printing so we can do custom rewriting
    build!(ctx, BUILDSTAGE_NO_PRINTING)

    rewrite!(ctx.dag)

    # print
    @test_logs (:info, "Done!") match_mode=:any build!(ctx, BUILDSTAGE_PRINTING_ONLY)
end


text:
const preferredTypes = Dict{Expr, Dict{String, Union{Symbol, Expr}}}()

function getPreferredType(T::Type, attr::String)::Union{Type, Nothing}
    MT = nameof(T)
    ST = Meta.quot(T.parameters[1])

    expr = :($MT{$ST})

    if haskey(preferredTypes, expr) && haskey(preferredTypes[expr], attr)
        return eval(preferredTypes[expr][attr])
    end

    return nothing
end

getPreferredType(obj, attr::String)::Union{Type, Nothing} = getPreferredType(typeof(obj), attr)

max_stars_repo_name: CelestialCartographers/Maple

max_issues_repo_name: CelestialCartographers/Maple

max_forks_repo_name: CelestialCartographers/Maple


text:
function deserialize_histogram_summary(summary)
    # custom deserialization
    if hasproperty(summary, :metadata)
        if summary.metadata.plugin_data.plugin_name == TB_PLUGIN_JLARRAY_NAME
            val = reshape(summary.histo.bucket,
                          reinterpret(Int,
                                      summary.metadata.plugin_data.content)...)

            return val
        end
    end

    # deserialize histogramproto
    hist_proto = summary.histo
    bin_edges = similar(hist_proto.bucket_limit, length(hist_proto.bucket_limit)+1)
    bin_edges[1] = hist_proto.min
    bin_edges[2:end] .= hist_proto.bucket_limit

    return Histogram(bin_edges, hist_proto.bucket)
end

# Lookahead for histograms .
# Checks the next entry if it is /re and /im, and combines them
function lookahead_deserialize_histogram_summary(old_tag, old_val, evs::Summary,
                                                 state)
    # prepare the default output (when no action is taken)
    result = old_tag, old_val, state
    # iterate to the next element
    res = iterate(evs, state + 1)
    res == nothing && return result

    # if the next event is identified, check its type
    (new_tag, summary), i_state = res
    typ = summary_type(summary)

    # if types match, check tags. If they match, return modified structure
    if typ === :histo && tags_match_re_im(old_tag, new_tag)
        val_im = deserialize_histogram_summary(summary)
        result = new_tag[1:end-3], old_val + im*val_im, i_state - 1
    else
        return result
    end
end


text:
using Test, JLD2, LinearAlgebra
#
# gradient computation by parial automatic differentiation and adjoint merthods
# of
#           E(v) = ∫ |U(v)|² dx
# where
# A(v) U(v) = F(v) the discrete version of - div(M ∇ Ua) = Fa
# and v are the vertices of the voronoi cells
#

# load polygonal mesh
filename, nc = "Lpolmesh_coarse", 100
mesh_filename = "$(@__DIR__)/data/$(filename)_$(nc).jld2"
println("read file $(mesh_filename)")
JLD2.@load(mesh_filename, pv, cellsb, cellsbt, t, pb, tb)
# boundary points / rhs
rhs, boundary_condition, meshboundary = if occursin("square", mesh_filename)
    PolygonalFem.rhs_sqr, PolygonalFem.boundary_condition_sqr, []
else
    EPS = 1e-4
    PolygonalFem.rhs_L_elas, PolygonalFem.boundary_condition_L_elas,
    findall(pv[:, 2] .> (1 - EPS))
end
np, ε = size(pv, 1), 1e-4
x, dp = pv, 2 * rand(size(pv)...) .- 1
# IMPORTANT Do not move boundary nodes
dp[meshboundary, :] .= 0.
xp, xm = x .+ ε * dp, x .- ε * dp
# finite differences
f = x -> PolygonalFem.costnormU_elas2D(x, cellsb, t, rhs, boundary_condition,
                                       meshboundary)
fp, fm = f(xp), f(xm)
@time g = PolygonalFem.∇costnormU_elas2D(x, cellsb, t, rhs, boundary_condition,
                                         meshboundary)
@time g = PolygonalFem.∇costnormU_elas2D(x, cellsb, t, rhs, boundary_condition,
                                         meshboundary)
println("")
println(fp)
println(fm)
println(" ")
println((fp - fm) / (2 * ε))
println(dot(g, dp))
println("")
@testset "test eval gradient dot function for polygonal costnormU_elas2D" begin
    @test abs(dot(g, dp) - (fp - fm) / (2 * ε)) < 1e-4
end


text:
#! /usr/bin/julia

# Rosetta Code, Day of the week

isdefined(:Date) || using Dates

lo = 2008
hi = 2121
xmas = Date(lo, 12, 25):Year(1):Date(hi, 12, 25)

smas = recur(xmas) do y
    Dates.dayofweek(y) == Dates.Sun
end

println("Years (from ", lo, " to ", hi, ") having Christmas on a Sunday:")

for y in smas
    println("    ", Dates.year(y))
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
struct ElectricPotential{T, N, S} <: AbstractArray{T, N}
    data::Array{T, N}
    grid::Grid{T, N, S}
end

@inline size(ep::ElectricPotential{T, N, S}) where {T, N, S} = size(ep.data)
@inline length(ep::ElectricPotential{T, N, S}) where {T, N, S} = length(ep.data)
@inline getindex(ep::ElectricPotential{T, N, S}, I::Vararg{Int, N}) where {T, N, S} = getindex(ep.data, I...)
@inline getindex(ep::ElectricPotential{T, N, S}, i::Int) where {T, N, S} = getindex(ep.data, i)
@inline getindex(ep::ElectricPotential{T, N, S}, s::Symbol) where {T, N, S} = getindex(ep.grid, s)

"""
    ElectricPotential(setup::PotentialSimulationSetup{T, 3, :cylindrical} ; kwargs...)::ElectricPotential{T, 3, :cylindrical}

Extracts the electric potential from `setup` and extrapolate it to an 2π grid.

For 2D grids (r and z) the user has to set the keyword `n_points_in_φ::Int`, e.g.: `n_points_in_φ = 36`.
"""
function ElectricPotential(setup::PotentialSimulationSetup{T, 3, :cylindrical} ; kwargs...)::ElectricPotential{T, 3, :cylindrical} where {T}
    return get_2π_potential(ElectricPotential{T, 3, :cylindrical}(setup.potential, setup.grid); kwargs...)
end

"""
    ElectricPotential(setup::PotentialSimulationSetup{T, 3, :cartesian} ; kwargs...)::ElectricPotential{T, 3, :cartesian}

Extracts the electric potential from `setup`.
"""
function ElectricPotential(setup::PotentialSimulationSetup{T, 3, :cartesian} )::ElectricPotential{T, 3, :cartesian} where {T}
    return ElectricPotential{T, 3, :cartesian}(setup.potential, setup.grid)
end

function NamedTuple(ep::ElectricPotential{T, 3}) where {T}
    return (
        grid = NamedTuple(ep.grid),
        values = ep.data * u"V",
    )
end
Base.convert(T::Type{NamedTuple}, x::ElectricPotential) = T(x)

function ElectricPotential(nt::NamedTuple)
    grid = Grid(nt.grid)
    T = typeof(ustrip(nt.values[1]))
    S = get_coordinate_type(grid)
    N = get_number_of_dimensions(grid)
    ElectricPotential{T, N, S}( ustrip.(uconvert.(u"V", nt.values)), grid)
end
Base.convert(T::Type{ElectricPotential}, x::NamedTuple) = T(x)

@recipe function f( ep::ElectricPotential{T, 3, :cartesian};
                    # dim = missing, dimvalue = missing,
                    x = missing,
                    y = missing,
                    z = missing,
                    contours_equal_potential=false ) where {T}
    g::Grid{T, 3, :cartesian} = ep.grid

    seriescolor --> :viridis
    st --> :heatmap
    aspect_ratio --> 1
    foreground_color_border --> nothing
    tick_direction --> :out

    cross_section::Symbol, idx::Int = if ismissing(x) && ismissing(y) && ismissing(z)
        :x, 1
    elseif !ismissing(x) && ismissing(y) && ismissing(z)
        :x, searchsortednearest(g[:x], T(x))
    elseif ismissing(x) && !ismissing(y) && ismissing(z)
        :y, searchsortednearest(g[:y], T(y))
    elseif ismissing(x) && ismissing(y) && !ismissing(z)
        :z, searchsortednearest(g[:z], T(z))
    else
        error(ArgumentError, ": Only one of the keywords `r, y, z` is allowed.")
    end
    value::T = if cross_section == :x
        g[:x][idx]
    elseif cross_section == :y
        g[:y][idx]
    elseif cross_section == :z
        g[:z][idx]
    end

    @series begin
        title --> "Electric Potential @$(cross_section) = $(round(value, sigdigits = 2))"
        if cross_section == :x
            xlabel --> "y / m"
            ylabel --> "z / m"
            g[:y], g[:z], ep.data[idx, :, :]'
        elseif cross_section == :y
            xlabel --> "x / m"
            ylabel --> "z / m"
            g[:x], g[:z], ep.data[:, idx, :]'
        elseif cross_section == :z
            xlabel --> "x / m"
            ylabel --> "y / m"
            g[:x], g[:y], ep.data[:,:,idx]
        end
    end
    if contours_equal_potential
        @series begin
            seriescolor := :thermal
            st := :contours
            if cross_section == :x
                g[:y], g[:z], ep.data[idx, :, :]'
            elseif cross_section == :y
                g[:x], g[:z], ep.data[:, idx, :]'
            elseif cross_section == :z
                g[:x], g[:y], ep.data[:,:,idx]
            end
        end
    end
end

@recipe function f( ep::ElectricPotential{T, 3, :cylindrical};
                    r = missing,
                    φ = missing,
                    z = missing,
                    contours_equal_potential=false ) where {T}
    g::Grid{T, 3, :cylindrical} = ep.grid

    seriescolor --> :viridis
    st --> :heatmap
    aspect_ratio --> 1
    foreground_color_border --> nothing
    tick_direction --> :out

    cross_section::Symbol, idx::Int = if ismissing(φ) && ismissing(r) && ismissing(z)
        :φ, 1
    elseif !ismissing(φ) && ismissing(r) && ismissing(z)
        φ_rad::T = T(deg2rad(φ))
        while !(g[:φ].interval.left <= φ_rad <= g[:φ].interval.right) && g[:φ].interval.right != g[:φ].interval.left
            if φ_rad > g[:φ].interval.right
                φ_rad -= g[:φ].interval.right - g[:φ].interval.left
            elseif φ_rad < g[:φ].interval.left
                φ_rad += g[:φ].interval.right - g[:φ].interval.left
            end
        end
        :φ, searchsortednearest(g[:φ], φ_rad)
    elseif ismissing(φ) && !ismissing(r) && ismissing(z)
        :r, searchsortednearest(g[:r], T(r))
    elseif ismissing(φ) && ismissing(r) && !ismissing(z)
        :z, searchsortednearest(g[:z], T(z))
    else
        error(ArgumentError, ": Only one of the keywords `r, φ, z` is allowed.")
    end
    value::T = if cross_section == :φ
        g[:φ][idx]
    elseif cross_section == :r
        g[:r][idx]
    elseif cross_section == :z
        g[:z][idx]
    end

    @series begin
        if cross_section == :φ
            title --> "Electric Potential @$(cross_section) = $(round(rad2deg(value), sigdigits = 2))"
            xlabel --> "r / m"
            ylabel --> "z / m"
            size --> ( 400, 350 / (g[:r][end] - g[:r][1]) * (g[:z][end] - g[:z][1]) )
            g[:r], g[:z], ep.data[:, idx,:]'
        elseif cross_section == :r
            title --> "Electric Potential @$(cross_section) = $(round(value, sigdigits = 2))"
            g[:φ], g[:z], ep.data[idx,:,:]'
        elseif cross_section == :z
            title --> "Electric Potential @$(cross_section) = $(round(value, sigdigits = 2))"
            proj --> :polar
            g[:φ], g[:r], ep.data[:,:,idx]
        end
    end
    if contours_equal_potential
        @series begin
            seriescolor := :thermal
            st := :contours
            if cross_section == :φ
                g[:r], g[:z], ep.data[:, idx,:]'
            elseif cross_section == :r
                g[:φ], g[:z], ep.data[idx,:,:]'
            elseif cross_section == :z
                proj --> :polar
                g[:φ], g[:r], ep.data[:,:,idx]
            end
        end
    end
end


text:
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("colorscheme.jl"))
using InteractiveDynamics, DynamicalSystems
import GLMakie

xmax = 1.0
xmin = -1.0
ymin = -1.0
ymax = 1.0
d = 50
const a = 1.4
const b = 0.3

ics = [
    SVector(x, y) for x in range(xmin, xmax; length = d) for y in range(ymin, ymax; length = d)
]

f1(u) = SVector(u[1], 1 - a*u[1]^2 + u[2])
f2(u) = SVector(b*u[1], u[2])
f3(u) = SVector(u[2], u[1])
function three_phases(ics)
    p1 = f1.(ics)
    p2 = f2.(p1)
    p3 = f3.(p2)
    return p1, p2, p3
end

function linear_transition!(o, p1, p2, steps, io)
    for i in 1:steps # assumes we start with o[] = p1
        o[] = p1 .+ (p2 .- p1) .* i/steps
        GLMakie.recordframe!(io)
        # sleep(0.01) # change to record frame
    end
end

iterations = 10
steps = 20

fig = GLMakie.Figure(); display(fig)
ax = GLMakie.Axis(fig[1,1])
o = GLMakie.Observable(ics)
GLMakie.scatter!(ax, o;
    color = COLORS[1], strokewidth = 0.5, strokecolor = :black, markersize = 5
)
ax.limits = ((-2.5,2.5),(-1.5,1.5))

GLMakie.record(fig, string(@__FILE__)[1:end-2]*"mp4"; framerate = 30) do io
    for j in 1:iterations
        p1, p2, p3 = three_phases(ics)
        linear_transition!(o, ics, p1, steps, io)
        GLMakie.recordframe!(io)
        linear_transition!(o, p1, p2, steps, io)
        GLMakie.recordframe!(io)
        linear_transition!(o, p2, p3, steps, io)
        GLMakie.recordframe!(io)
        global ics = p3
        GLMakie.recordframe!(io)
    end
end

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
import LinearAlgebra: I, Symmetric, cond, diagm, dot
using RunLengthArrays
import Healpix

################################################################################
# Pointing generation (used for testing)

function build_F(ntod, baseline_lengths)
    F = zeros(ntod, length(baseline_lengths))

    startidx = 1
    for (i, cur_len) in enumerate(baseline_lengths)
        nelems = min(ntod - startidx + 1, cur_len)
        F[startidx:(startidx + nelems - 1), i] .= 1
        startidx += nelems
    end

    F
end

function build_P_and_psi(pixidx, npix)
    ntod = length(pixidx)
    P = zeros(ntod, 3npix)
    psi = Array{Float64}(undef, ntod)

    curangle = 0.0
    for (idx, curpix) in enumerate(pixidx)
        psi[idx] = curangle

        P[idx, 1 + 3(curpix - 1)] = 1
        P[idx, 2 + 3(curpix - 1)] = cos(2curangle)
        P[idx, 3 + 3(curpix - 1)] = sin(2curangle)

        curangle += 2π / ntod
    end

    P, psi
end

generate_pixidx(ntod, npix) = Iterators.take(Iterators.cycle(1:npix), ntod) |> collect

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

# Here we build a number of vectors and matrices that solve the full map-making
# equation by explicitly building all the matrices in full form. The values here
# will be compared with the ones produced by the map-maker.
#
# Unlike the destriping code, we do *not* split baselines/TODs/etc. among MPI
# processes: each MPI process carries the full computation.

NTOD = 512
BASELINES = RunLengthArray{Int, Float64}(
    Int[NTOD/8, NTOD/4, NTOD/4, 3NTOD/8],
    Float64[2, -3, 1, 0],
)

# Map with NSIDE = 1
skymap = Healpix.PolarizedMap{Float64, Healpix.RingOrder}(
    Float64[1,   3,   -4,    6,    5,   -1,   -4,    2,   8,   0,   -1,   -2],
    Float64[1.0, 1.8,  0.7,  6.8, -2.2,  6.0,  0.3, -0.3, 4.7, 7.2,  5.1,  7.2],
    Float64[1.2, 0.4,  0.7, -1.4,  5.4,  4.9,  4.1, -1.0, 4.2, 5.9,  2.4,  3.8],
)

# The order is I1, Q1, U1, I2, Q2, U2…
m = reshape([skymap.i skymap.q skymap.u]', 3length(skymap.i))

NPIX = length(skymap.i)
PIXIDX = generate_pixidx(NTOD, NPIX)

diagCw = [1.1 for _ in 1:NTOD]
Cw = diagm(0 => diagCw)
invCw = inv(Cw)

P, PSI = build_P_and_psi(PIXIDX, NPIX)
F = build_F(NTOD, runs(BASELINES))
M = P' * invCw * P
invM = inv(M)
Z = I - P * invM * P' * invCw
D = F' * invCw * Z * F

y = P * m + F * values(BASELINES)

# These are the terms in the destriping equation "Ax = b"
A = F' * invCw * Z * F
b = F' * invCw * Z * y

baseline_guess = inv(F' * invCw * Z * F) * F' * invCw * Z * y
# Remove the average
baseline_guess .-= sum(baseline_guess) / length(baseline_guess)
# This is the maximum-likelihood map
map_guess = inv(P' * invCw * P) * P' * invCw * (y - F * baseline_guess)

################################################################################
# True map-maker

# The code that follows structures the data depending whether we are
# using MPI or not:
#
# 1. When we are *not* using MPI, there are two observations
# 2. When we are using MPI, we force the usage of two processes, and each
#    process holds only one observation.

obs_range1 = 1:sum(runs(BASELINES)[1:2])
obs_range2 = (sum(runs(BASELINES)[1:2]) + 1):NTOD

mpicomm = if Harlequin.use_mpi()
    # If use_mpi returns true, this means that this file
    # ("test_mapmaking.jl") is being run within "test_mpi.jl", which
    # has already imported the MPI library. So we can freely use
    # functions in the MPI package

    MPI.Init()
    MPI.COMM_WORLD
else
    nothing
end

# Range of validity of the array of baselines
baseline_range = if Harlequin.use_mpi()
    (MPI.Comm_rank(mpicomm) == 0) ? (1:2) : (3:4)
else
    1:4
end

rank_baselines = if Harlequin.use_mpi()
    @assert MPI.Comm_size(mpicomm) == 2

    RunLengthArray{Int, Float64}(
        runs(BASELINES)[baseline_range],
        values(BASELINES)[baseline_range],
    )
else
    BASELINES
end

rank_runs = if Harlequin.use_mpi()
    [runs(rank_baselines)]
else
    [runs(rank_baselines)[1:2], runs(rank_baselines)[3:4]]
end

rank_baselines = if Harlequin.use_mpi()
    [collect(rank_baselines)]
else
    [collect(BASELINES)[obs_range1], collect(BASELINES)[obs_range2]]
end

time_samp = 0.01 # Integration time for one sample; it is not really used

# This holds the same information as matrix M = P' ⋅ Cw^-1 ⋅ P
nobs_matrix = [NobsMatrixElement{Float64}() for i in 1:NPIX]

obs_list = [
    Observation{Float64}(
        pixidx = PIXIDX[obs_range1],
        psi_angle = PSI[obs_range1],
        tod = y[obs_range1],
        sigma_squared = diagCw[obs_range1],
        name = "Det1",
    )
    Observation{Float64}(
        pixidx = PIXIDX[obs_range2],
        psi_angle = PSI[obs_range2],
        tod = y[obs_range2],
        sigma_squared = diagCw[obs_range2],
        name = "Det2",
    )
]

# When using MPI, we only use two processes and assign the two
# observations to each of them
obs_list = Harlequin.use_mpi() ? [obs_list[MPI.Comm_rank(mpicomm) + 1]] : obs_list

let io = IOBuffer()
    show(io, obs_list[1])
    show(io, "text/plain", obs_list[1])
end

compute_nobs_matrix!(nobs_matrix, obs_list, comm = mpicomm)

# Test that NobsMatrixElement objects can be showed
let io = IOBuffer()
    show(io, nobs_matrix[1])
    show(io, "text/plain", nobs_matrix[1])
end

# Verify that nobs_matrix is really a representation of matrix M
for i in 1:NPIX
    startidx = 1 + 3(i - 1)
    subm = @view M[startidx:(startidx + 2), startidx:(startidx + 2)]
    @test nobs_matrix[i].invm ≈ inv(subm)
    @test Symmetric(nobs_matrix[i].m) * nobs_matrix[i].invm ≈ I
end

# Test compute_z_and_group!

@test length(runs(BASELINES)) == 4

dest_baselines = if use_mpi()
    [
        Float64[0.0, 0.0],
    ]
else
    [
        Float64[0.0, 0.0],
        Float64[0.0, 0.0],
    ]
end

destriping_data = DestripingData{Float64, Healpix.RingOrder}(
    1,
    obs_list,
    rank_runs,
    max_iterations = 15,
    threshold = 1e-14,
    use_preconditioner = true,
    comm = mpicomm,
)

# Check that the "show" method works in all its forms

let io = IOBuffer()
    show(io, destriping_data)
    show(io, "text/plain", destriping_data)
end

compute_z_and_group!(
    dest_baselines,
    rank_runs,
    destriping_data,
    obs_list,
    comm = mpicomm,
)

@test collect(Iterators.flatten(dest_baselines)) ≈ b[baseline_range]

left_side = if use_mpi()
    [
        Float64[0.0, 0.0],
    ]
else
    [
        Float64[0.0, 0.0],
        Float64[0.0, 0.0],
    ]
end

reset_maps!(destriping_data)

compute_z_and_group!(
    left_side,
    rank_baselines,
    rank_runs,
    destriping_data,
    obs_list,
    comm = mpicomm,
)

@test collect(Iterators.flatten(left_side)) ≈ (A * values(BASELINES))[baseline_range]

reset_maps!(destriping_data)

destripe!(
    obs_list,
    destriping_data,
    comm = mpicomm,
)

estimated_baselines = collect(Iterators.flatten([x.values for x in destriping_data.baselines]))
@test estimated_baselines ≈ values(BASELINES)[baseline_range]

@test (abs.(destriping_data.skymap.i .- skymap.i) |> maximum) < 0.3
@test (abs.(destriping_data.skymap.q .- skymap.q) |> maximum) < 0.3
@test (abs.(destriping_data.skymap.u .- skymap.u) |> maximum) < 0.3


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

JLLWrappers.@generate_main_file_header("ABC")
JLLWrappers.@generate_main_file("ABC", UUID("73717f16-8662-5d85-9e79-a3c1fc1c8861"))
end  # module ABC_jll


text:
module CoupledOscillators
using OrdinaryDiffEq

include("kuramoto.jl")
include("strogatz.jl")

end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -7.9284408379108005 <= q <= 14.377266175441108)

                     add_NL_constraint(m, :(log(1 + exp(0.4343485145630286 + -0.5735104656090533*log(1 + exp(-0.47832195115882303 + 0.6358617391786257*$(x[1]) + -0.4037373919488174*$(x[2]))) + -0.40833038450807946*log(1 + exp(0.4113265195541662 + -0.7994001067331897*$(x[1]) + -0.3401124060181533*$(x[2]))) + 0.09714914815406805*log(1 + exp(-0.3249145933093094 + 0.9313960721440679*$(x[1]) + 0.2643917347767659*$(x[2]))) + -0.19623354791767778*log(1 + exp(0.34085895078435247 + -0.33650259644885194*$(x[1]) + -0.9436140496760688*$(x[2]))) + -0.3237276355524821*log(1 + exp(0.8342164127508305 + -0.6721074785776127*$(x[1]) + 0.009028597341863787*$(x[2]))))) + log(1 + exp(-0.11905941659917563 + 0.17224590981590104*log(1 + exp(-0.47832195115882303 + 0.6358617391786257*$(x[1]) + -0.4037373919488174*$(x[2]))) + 0.8406566605626429*log(1 + exp(0.4113265195541662 + -0.7994001067331897*$(x[1]) + -0.3401124060181533*$(x[2]))) + 0.24719777641241691*log(1 + exp(-0.3249145933093094 + 0.9313960721440679*$(x[1]) + 0.2643917347767659*$(x[2]))) + -0.24308200319021767*log(1 + exp(0.34085895078435247 + -0.33650259644885194*$(x[1]) + -0.9436140496760688*$(x[2]))) + -0.39557038452022475*log(1 + exp(0.8342164127508305 + -0.6721074785776127*$(x[1]) + 0.009028597341863787*$(x[2]))))) + log(1 + exp(0.967203829156619 + 0.3785230261408685*log(1 + exp(-0.47832195115882303 + 0.6358617391786257*$(x[1]) + -0.4037373919488174*$(x[2]))) + -0.23118999327456224*log(1 + exp(0.4113265195541662 + -0.7994001067331897*$(x[1]) + -0.3401124060181533*$(x[2]))) + -0.5197073892915105*log(1 + exp(-0.3249145933093094 + 0.9313960721440679*$(x[1]) + 0.2643917347767659*$(x[2]))) + 0.6005003791404913*log(1 + exp(0.34085895078435247 + -0.33650259644885194*$(x[1]) + -0.9436140496760688*$(x[2]))) + 0.15463985115916623*log(1 + exp(0.8342164127508305 + -0.6721074785776127*$(x[1]) + 0.009028597341863787*$(x[2]))))) + log(1 + exp(-0.18340968019959902 + 0.01980109053282808*log(1 + exp(-0.47832195115882303 + 0.6358617391786257*$(x[1]) + -0.4037373919488174*$(x[2]))) + 0.8767310690037822*log(1 + exp(0.4113265195541662 + -0.7994001067331897*$(x[1]) + -0.3401124060181533*$(x[2]))) + 0.42243485881005904*log(1 + exp(-0.3249145933093094 + 0.9313960721440679*$(x[1]) + 0.2643917347767659*$(x[2]))) + -0.1694649191201152*log(1 + exp(0.34085895078435247 + -0.33650259644885194*$(x[1]) + -0.9436140496760688*$(x[2]))) + 0.5278959282338329*log(1 + exp(0.8342164127508305 + -0.6721074785776127*$(x[1]) + 0.009028597341863787*$(x[2]))))) + log(1 + exp(0.7043455967116397 + -0.48930630508951234*log(1 + exp(-0.47832195115882303 + 0.6358617391786257*$(x[1]) + -0.4037373919488174*$(x[2]))) + 0.7375930151357304*log(1 + exp(0.4113265195541662 + -0.7994001067331897*$(x[1]) + -0.3401124060181533*$(x[2]))) + -0.9188673755139103*log(1 + exp(-0.3249145933093094 + 0.9313960721440679*$(x[1]) + 0.2643917347767659*$(x[2]))) + 0.5801272291202126*log(1 + exp(0.34085895078435247 + -0.33650259644885194*$(x[1]) + -0.9436140496760688*$(x[2]))) + 0.06743755305487698*log(1 + exp(0.8342164127508305 + -0.6721074785776127*$(x[1]) + 0.009028597341863787*$(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