text:
@doc """
# Mini-SEED Support
SeisIO supports mini-SEED, the "data-only" extension of the SEED (Standard for
the Exchange of Earthquake Data) file format.

## Supported Blockette Types

| Blockette                                     | Key in `:misc`  |
|:----------------------------------------------|:----------------|
| [100] Sample Rate Blockette                   |                 |
| [201] Murdock Event Detection Blockette       | seed_event      |
| [300] Step Calibration Blockette              | seed_calib      |
| [310] Sine Calibration Blockette              | seed_calib      |
| [320] Pseudo-random Calibration Blockette     | seed_calib      |
| [390] Generic Calibration Blockette           | seed_calib      |
| [500] Timing Blockette                        | seed_timing     |
| [1000] Data Only SEED Blockette               |                 |
| [1001] Data Extension Blockette               |                 |
| [2000] Variable Length Opaque Data Blockette  | seed_opaque     |

Information unrelated to data or timing is stored in `:misc` as String arrays;
each blockette gets a single String in the named key, separated by a newline
character (\\n). `:misc["seed_opaque"]`, if present, contains raw (UInt8) byte
vectors of all data in each packet.

## Supported Data Encodings

| Format  | Data Encoding                                               |
|---------|:------------------------------------------------------------|
| 0       | ASCII text [^a]                                             |
| 1       | Int16 unencoded                                             |
| 3       | Int32 unencoded                                             |
| 4       | Float32 unencoded                                           |
| 5       | Float64 unencoded [^b]                                      |
| 10      | Steim-1                                                     |
| 11      | Steim-2                                                     |
| 12      | GEOSCOPE multiplexed, 16-bit gain ranged, 3-bit exponent    |
| 15      | GEOSCOPE multiplexed, 16-bit gain ranged, 4-bit exponent    |
| 16      | CDSN, 16-bit gain ranged                                    |
| 30      | SRO                                                         |
| 32      | DWWSSN gain ranged                                          |

[^a]: Saved to `C.misc["seed_ascii"]`; generally not used for data
[^b]: Converted to Float32
""" mseed_support
function mseed_support()
  println("\nCurrently supported mini-SEED blockette types:\n")
  printstyled(@sprintf("%10s|%40s|%18s\n", "BLOCKETTE", " ", "KEY IN :MISC"), color=:green, bold=true)
  printstyled(@sprintf("%10s|%40s|%18s\n", "[100]", "Sample Rate", " "), color=:green)
  printstyled(@sprintf("%10s|%40s|%18s\n", "[201]", "Murdock Event Detection", "seed_event"), color=:green)
  printstyled(@sprintf("%10s|%40s|%18s\n", "[300]", "Step Calibration", "seed_calib"), color=:green)
  printstyled(@sprintf("%10s|%40s|%18s\n", "[310]", "Sine Calibration", "seed_calib"), color=:green)
  printstyled(@sprintf("%10s|%40s|%18s\n", "[320]", "Pseudo-random Calibration", "seed_calib"), color=:green)
  printstyled(@sprintf("%10s|%40s|%18s\n", "[390]", "Generic Calibration", "seed_calib"), color=:green)
  printstyled(@sprintf("%10s|%40s|%18s\n", "[500]", "Timing", "seed_timing"), color=:green)
  printstyled(@sprintf("%10s|%40s|%18s\n", "[1000]", "Data Only SEED", " "), color=:green)
  printstyled(@sprintf("%10s|%40s|%18s\n", "[1001]", "Data Extension", " "), color=:green)
  printstyled(@sprintf("%10s|%40s|%18s\n", "[2000]", "Variable Length Opaque Data","seed_opaque*"), color=:green)

  return nothing
end

@doc """
# Dataless SEED Support
SeisIO supports reading of dataless SEED meta-data files.

Reading full SEED volumes is NYI.

## Supported Blockette Types

| Blockette                                             |
|:------------------------------------------------------|
| [10] Volume Identifier Blockette                      |
| [11] Volume Station Header Index Blockette            |
| [12] Volume Time Span Index Blockette                 |
| [31] Comment Description Blockette                    |
| [33] Generic Abbreviation Blockette                   |
| [34] Units Abbreviations Blockette                    |
| [41] FIR Dictionary Blockette                         |
| [43] Response (Poles & Zeros) Dictionary Blockette    |
| [44] Response (Coefficients) Dictionary Blockette     |
| [45] Response List Dictionary Blockette               |
| [47] Decimation Dictionary Blockette                  |
| [48] Channel Sensitivity/Gain Dictionary Blockette    |
| [50] Station Identifier Blockette                     |
| [52] Channel Identifier Blockette                     |
| [53] Response (Poles & Zeros) Blockette               |
| [54] Response (Coefficients) Blockette                |
| [57] Decimation Blockette                             |
| [58] Channel Sensitivity/Gain Blockette               |
| [59] Channel Comment Blockette                        |
| [60] Response Reference Blockette                     |
| [61] FIR Response Blockette                           |

## Unsupportable Blockette Types
These blockettes will probably never be supported as
their information lies outside the scope of SeisIO. At
high verbosity (v > 2), their information is dumped to
stdout.

| Blockette                                             |
|:------------------------------------------------------|
| [30] Data Format Dictionary Blockette                 |
| [32] Cited Source Dictionary Blockette                |
| [51] Station Comment Blockette                        |

"""
seed_support
function seed_support()
  println("\nCurrently supported dataless SEED blockette types:\n")
  printstyled("[10] Volume Identifier Blockette\n", color=:green)
  printstyled("[11] Volume Station Header Index Blockette\n", color=:green)
  printstyled("[12] Volume Time Span Index Blockette\n", color=:green)
  printstyled("[31] Comment Description Blockette\n", color=:green)
  printstyled("[33] Generic Abbreviation Blockette\n", color=:green)
  printstyled("[34] Units Abbreviations Blockette\n", color=:green)
  printstyled("[41] FIR Dictionary Blockette\n", color=:green)
  printstyled("[43] Response (Poles & Zeros) Dictionary Blockette\n", color=:green)
  printstyled("[44] Response (Coefficients) Dictionary Blockette\n", color=:green)
  printstyled("[45] Response List Dictionary Blockette\n", color=:green)
  printstyled("[47] Decimation Dictionary Blockette\n", color=:green)
  printstyled("[48] Channel Sensitivity/Gain Dictionary Blockette\n", color=:green)
  printstyled("[50] Station Identifier Blockette\n", color=:green)
  printstyled("[52] Channel Identifier Blockette\n", color=:green)
  printstyled("[53] Response (Poles & Zeros) Blockette\n", color=:green)
  printstyled("[54] Response (Coefficients) Blockette\n", color=:green)
  printstyled("[57] Decimation Blockette\n", color=:green)
  printstyled("[58] Channel Sensitivity/Gain Blockette\n", color=:green)
  printstyled("[59] Channel Comment Blockette\n", color=:green)
  printstyled("[60] Response Reference Blockette\n", color=:green)
  printstyled("[61] FIR Response Blockette\n", color=:green)
end

#
# ### Unsupported Data Encodings
# These have never been encountered by SeisIO and may not exist in real seismic
# data.
#
# | Format  | Data Encoding                                               |
# |---------|:------------------------------------------------------------|
# | 2       | Int24 unencoded                                             |
# | 12      | GEOSCOPE multiplexed, 24-bit integer                        |
# | 15      | US National Network                                         |
# | 17      | Graefenberg, 16-bit gain ranged                             |
# | 18      | IPG - Strasbourg, 16-bit gain ranged                        |
# | 19      | Steim-3                                                     |
# | 31      | HGLP                                                        |
# | 33      | RSTN 16-bit gain ranged                                     |


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "x264Builder"
version = v"2019.05.25"

# Collection of sources required to build x264Builder
sources = [
    "https://download.videolan.org/pub/videolan/x264/snapshots/x264-snapshot-20190525-2245-stable.tar.bz2" =>
    "638581a18bff8e9375211955422eff145011c8ccfd0994d43bd194cd82984f7a",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd x264-snapshot-20190525-2245-stable/
if [[ "${target}" == x86_64* ]] || [[ "${target}" == i686* ]]; then
    apk add nasm
    export AS=nasm
else
    export AS="${CC}"
fi
./configure --prefix=$prefix --host=$target --enable-static --enable-pic
make -j${nproc}
make install
"""

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

# The products that we will ensure are always built
products(prefix) = [
    ExecutableProduct(prefix, "x264", :x264)
]

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

]

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

ext: jl


text:
@testset "Poisson{<:ScaledLogisticLink}" begin
    test_auglik(PoissonLikelihood(ScaledLogistic(10.0)))
end


text:
"""
    StVenant(; E, ν)

Hyperelastic material
#Arguments
- `λ::Float64`: Lamé parameter
- `μ::Float64`: Lamé parameter (shear modulus)
"""

struct StVenant <: AbstractMaterial
    λ::Float64
    μ::Float64
end

struct StVenantState <: AbstractMaterialState
end

function initial_material_state(::StVenant)
    return StVenantState()
end

function StVenant(; λ::T, μ::T) where T
    return StVenant(λ, μ)
end

strainmeasure(::StVenant) = RightCauchyGreen()

function ψ(mp::StVenant, C)
    μ = mp.μ
    λ = mp.λ
    Ic = tr(C)
    IIc = tr(det(C) * inv(C)')
    return λ / 8 * (Ic - 3)^2 + μ / 4 * (Ic^2 - 2 * Ic - 2 * IIc + 3)
end

function material_response(mp::StVenant, C::SymmetricTensor{2,3}, state::StVenantState = StVenantState(),
                           Δt=nothing; cache=nothing, options=nothing)

    μ = mp.μ
    λ = mp.λ
    I = one(SymmetricTensor{2,3})

    S = λ/2 * (tr(C) - 3)*I + μ*(C-I)
    ∂S∂C = λ*(I⊗I) + μ*(otimesu(I,I) + otimesl(I,I))
    ∂S∂C = symmetric(∂S∂C)

    return S, ∂S∂C, StVenantState()
end


text:
export BrickletSoundIntensityIntensityCallbackThreshold
struct BrickletSoundIntensityIntensityCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletSoundIntensityIdentity
struct BrickletSoundIntensityIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletSoundIntensity
"""
Measures sound intensity
"""
mutable struct BrickletSoundIntensity <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletSoundIntensity(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_sound_intensity")
        deviceInternal = package.BrickletSoundIntensity(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_intensity
"""
    $(SIGNATURES)

Returns the current sound intensity.

The value corresponds to the
`upper envelop <https://en.wikipedia.org/wiki/Envelope_(waves)>`__
of the signal of the microphone capsule.

If you want to get the intensity periodically, it is recommended to use the
:cb:`Intensity` callback and set the period with
:func:`Set Intensity Callback Period`.
"""
function get_intensity(device::BrickletSoundIntensity)
    return device.deviceInternal.get_intensity()
end

export set_intensity_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Intensity` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Intensity` callback is only triggered if the intensity has changed
since the last triggering.
"""
function set_intensity_callback_period(device::BrickletSoundIntensity, period)
    device.deviceInternal.set_intensity_callback_period(period)
end

export get_intensity_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Intensity Callback Period`.
"""
function get_intensity_callback_period(device::BrickletSoundIntensity)
    return device.deviceInternal.get_intensity_callback_period()
end

export set_intensity_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Intensity Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the intensity is *outside* the min and max values"
 "'i'",    "Callback is triggered when the intensity is *inside* the min and max values"
 "'<'",    "Callback is triggered when the intensity is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the intensity is greater than the min value (max is ignored)"
"""
function set_intensity_callback_threshold(device::BrickletSoundIntensity, option, min, max)
    device.deviceInternal.set_intensity_callback_threshold(option, min, max)
end

export get_intensity_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Intensity Callback Threshold`.
"""
function get_intensity_callback_threshold(device::BrickletSoundIntensity)
    return device.deviceInternal.get_intensity_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callback

* :cb:`Intensity Reached`

is triggered, if the thresholds

* :func:`Set Intensity Callback Threshold`

keeps being reached.
"""
function set_debounce_period(device::BrickletSoundIntensity, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletSoundIntensity)
    return device.deviceInternal.get_debounce_period()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletSoundIntensity)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletSoundIntensity, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
using Documenter
using SimulationLogs

makedocs(
    modules = [SimulationLogs],
    format=Documenter.HTML(
        prettyurls = false,
        canonical = "https://jonniedie.github.io/SimulationLogs/stable",
    ),
    pages = [
        "Home" => "index.md",
        "Examples" => [
            "examples/cruise_control.md"
        ],
        "API" => "api.md",
    ],
    repo = "https://github.com/jonniedie/SimulationLogs/blob/{commit}{path}#L{line}",
    sitename = "SimulationLogs Documentation",
    authors = "Jonnie Diegelman",
    assets=String[],
)

deploydocs(
    repo = "github.com/jonniedie/SimulationLogs.git",
    devbranch = "main"
)

max_stars_repo_name: jonniedie/SimulationLogs

max_issues_repo_name: jonniedie/SimulationLogs

max_forks_repo_name: jonniedie/SimulationLogs


text:
module TestMacrosGroupBy

import Test: @testset, @test
import Volcanito: @group_by, materialize
import DataFrames: DataFrame, groupby

@testset "@group_by" begin
    df = DataFrame(
        a = [1, 2, 3, 4],
        b = [2, 3, 4, 5],
    )

    plan = @group_by(df, a)
    res = materialize(plan)
    @test res == groupby(df, :a)
end

end


text:
export Warrior
export Archer, Knight, Pikeman

export battle!, attack!, resupply!, shoot!, mount!, dismount!

abstract type Warrior end

mutable struct Archer <: Warrior
    name::String
    health::Int
    arrows::Int
end

mutable struct Pikeman <: Warrior
    name::String
    health::Int
end

mutable struct Knight <: Warrior
    name::String
    health::Int
    mounted::Bool
end

function Knight(name::AbstractString, health::Integer)
   Knight(name, health, true)
end

function Archer(name::AbstractString, health::Integer)
    archer = Archer(name, health, 0)
    resupply!(archer)
    archer
end

"""
    shoot(a::Archer)
Shoot an arrow. Deplete nubmer of arrows in quiver.
"""
function shoot!(archer::Archer)
    if archer.arrows > 0
        archer.arrows -= 1
    end
    archer
end

"""
    resupply!(a::Archer)
Ressuply archer `a` with arrows.
"""
function resupply!(archer::Archer)
    archer.arrows = 24
    archer
end

function mount!(knight::Knight)
    knight.mounted = true
    knight
end

function dismount!(knight::Knight)
    knight.mounted = false
    knight
end

"""
    battle!(a::Warrior, b::Warrior)
Display the outcome of warrior `a` attacking warrior `b`
"""
function battle!(a::Warrior, b::Warrior)
    attack!(a, b)
    if a.health == 0 && b.health == 0
        println(a.name, " and ", b.name, " destroyed each other")
    elseif a.health == 0
        println(b.name, " defeated ", a.name)
    elseif b.health == 0
        println(a.name, " defeated ", b.name)
    else
        println(b.name, " survived attack from ", a.name)
    end
end

"""
    attack!(a::Archer, b::Archer)
Archers can do a fair amount of damage on other archers as
the archers don't have good armor to protect them.
"""
function attack!(a::Archer, b::Archer)
    if a.arrows > 0
        shoot!(a)
        damage = 6 + rand(1:6)
        b.health = max(b.health - damage, 0)
    end
    a.health, b.health
end

"""
    attack!(a::Archer, b::Knight)
Archers can do some damage to a knight when when he is mounted
because the horse is vulnurable. Hurting dismounted knights is harder
as most arrows cannot pierce plate armor.
"""
function attack!(a::Archer, b::Knight)
    if a.arrows > 0
        shoot!(a)
        damage = rand(1:6)
        if b.mounted
            damage += 3
        end
        b.health = max(b.health - damage, 0)
    end
    a.health, b.health
end

"""
    attack!(a::Knight, b::Knight)
Knights vs knights is evenly matched and both have a fairly good chance of
disengaging so casualties are not that high.
"""
function attack!(a::Knight, b::Knight)
    a.health = max(a.health - rand(1:6), 0)
    b.health = max(b.health - rand(1:6), 0)
    a.health, b.health
end

"""
    attack!(a::Archer, b::Pikeman)
Because pikemen move slowly they cannot easily get away from the rain of
arrows and thus take significant damage.
"""
function attack!(a::Archer, b::Pikeman)
    if a.arrows > 0
        shoot!(a)
        damage = 4 + rand(1:6)
        b.health = max(b.health - damage, 0)
    end
    a.health, b.health
end

"""
    attack!(a::Pikeman, b::Archer)
It is hard for pikemen to attack archers are they move slowly. Thus
only when lucky do they succeed in causing some injury.
"""
function attack!(a::Pikeman, b::Archer)
   if rand(1:6) >= 5
       b.health = max(b.health - 2, 0)
   end
   a.health, b.health
end

"""
    attack!(a::Pikeman, b::Pikeman)
Pikemen vs pikemen attacks were generally avoided as
they tended to get  brutal and with heavy causualties.
It is hard to escape such a close encounter once it begins.
"""
function attack!(a::Pikeman, b::Pikeman)
    for _ in 1:3
        a.health = max(a.health - rand(1:6), 0)
        b.health = max(b.health - rand(1:6), 0)

        if a.health > 0 && b.health > 0
            break
        end
    end
    a.health, b.health
end

"""
    attack!(a::Pikeman, b::Knight)
Due to the superior mobility of knights, a pikeman cannot easily engage a
knight which does not choose to attack, unless the knight is dismounted.
"""
function attack!(a::Pikeman, b::Knight)
    if rand(1:6) >= 6
        damage = 1
        if b.mounted
            damage += rand(1:6)
        end
        b.health = max(b.health - damage, 0)
    end
    a.health, b.health
end

"""
    attack!(a::Knight, b::Archer)
Archers attacked by knights would suffer heavy losses as they have poor
armor.
"""
function attack!(a::Knight, b::Archer)
    damage = 8 + 2*rand(1:6)
    b.health = max(b.health - damage, 0)
    a.health, b.health
end

"""
    attack!(a::Knight, b::Pikeman)
Mounted knights attacking pikemen will usually suffer heavy losses
unless attacking on a flank or the rear. Unmounted knights are more
effective as pikemen cannot attack their vulnurable horses.
"""
function attack!(a::Knight, b::Pikeman)
    damage = rand(1:6)
    b.health = max(b.health - damage, 0)

    damage = rand(1:6)
    if a.mounted
        damage += 2*rand(1:6)
    end
    a.health = max(b.health - damage, 0)

    a.health, b.health
end

max_stars_repo_name: ordovician/BattleSim

max_issues_repo_name: ordovician/BattleSim

max_forks_repo_name: ordovician/BattleSim


text:
printstyled("  nanfill, ungap\n", color=:light_green)

# Check that nanfill does not error
S = randSeisData()
Ev = SeisEvent(hdr=randSeisHdr(), data=convert(EventTraceData, deepcopy(S)))
C = deepcopy(S[1])

for i = 1:S.n
  L = length(S.x[i])
  inds = rand(1:L, div(L,2))
  S.x[i][inds] .= NaN
end
nanfill!(S)
for i = 1:S.n
  x = getindex(getfield(S, :x), i)
  @test isempty(findall(isnan.(x)))
end

for i = 1:Ev.data.n
  L = length(Ev.data.x[i])
  inds = rand(1:L, div(L,2))
  Ev.data.x[i][inds] .= NaN
end
nanfill!(Ev.data)
for i = 1:Ev.data.n
  x = getindex(getfield(getfield(Ev, :data), :x), i)
  @test isempty(findall(isnan.(x)))
end

L = length(C.x)
inds = rand(1:L, div(L,2))
C.x[inds] .= NaN
nanfill!(C)
@test isempty(findall(isnan.(C.x)))

# Test that traces of all NaNs becomes traces of all zeros
C = SeisChannel()
C.x = fill!(zeros(Float32, 128), NaN32)
nanfill!(C)
@test C.x == zeros(Float32, 128)

S = randSeisData()
U = deepcopy(S)
for i = 1:S.n
  x = getindex(getfield(S, :x), i)
  nx = lastindex(x)
  T = eltype(x)
  fill!(x, T(NaN))

  u = getindex(getfield(U, :x), i)
  fill!(u, zero(T))
end
nanfill!(S)
for i = 1:S.n
  x = getindex(getfield(S, :x), i)
  nx = lastindex(x)
  T = eltype(x)
  u = getindex(getfield(U, :x), i)
  @test T == eltype(u)
  @test S.x[i] == x == U.x[i] == u
end

# Test that ungap calls nanfill properly
Ev2 = ungap(Ev.data, tap=true)
ungap!(Ev.data, tap=true)
@test Ev.data == Ev2
ungap!(C, tap=true)
ungap!(S, tap=true)

# Ensure one segment is short enough to invoke bad behavior in ungap
Ev = randSeisEvent()
Ev.data.fs[1] = 100.0
Ev.data.x[1] = rand(1024)
Ev.data.t[1] = vcat(Ev.data.t[1][1:1,:], [5 2*ceil(S.fs[1])*sμ], [8 2*ceil(S.fs[1])*sμ], [1024 0])

redirect_stdout(out) do
  ungap!(Ev.data)
end


text:
@testset "SeisCon" begin

    global s = segy_scan(joinpath(SegyIO.myRoot,"data/"), "overthrust", ["GroupX"; "GroupY"], verbosity = 0)
    @testset "Constructor" begin
        @test typeof(s) == SeisCon
        @test length(s) == 97
        @test size(s) == (97,)
        @test s.ns == 751
        @test s.dsf == 1
    end

    @testset "Indexing" begin
        @test typeof(s[1]) <: SeisBlock
        @test typeof(s[1:10]) <: SeisBlock
        @test typeof(s[1:10:20]) <: SeisBlock
        @test typeof(s[:]) <: SeisBlock
        @test typeof(s[[1; 4; 8; 12]]) <: SeisBlock
    end

    @testset "get_header" begin
        @test size(get_header(s, "GroupX")) == (97, 2)
        @test get_header(s, "GroupX")[1] == 2400
        @test get_header(s, "GroupX")[end] == 19900
    end

    @testset "get_sources" begin
        @test size(get_sources(s)) == (97, 2)
        @test get_sources(s)[1] == 8400
        @test get_sources(s)[end] == 0
    end

    @testset "merge" begin
        b = merge(s, s)
        @test b.ns == s.ns
        @test b.dsf == s.dsf
        @test length(b) == 194
    end

    @testset "split" begin
        b = split(s, 1:10);
        c = split(b, [1; 7; 9]);
        d = split(c, 1);
        @test s.blocks[1] === d.blocks[1]
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

using FEMMaterials, Test
using DelimitedFiles

path = joinpath("test_chaboche", "unitelement_results.rpt")
data = readdlm(path, Float64; skipstart=4)
ts = data[:,1]
s11_ = data[:,2]
s12_ = data[:,3]
s13_ = data[:,4]
s22_ = data[:,5]
s23_ = data[:,6]
s33_ = data[:,7]
e11_ = data[:,8]
e12_ = data[:,9]
e13_ = data[:,10]
e22_ = data[:,11]
e23_ = data[:,12]
e33_ = data[:,13]

strains = [[e11_[i], e22_[i], e33_[i], e23_[i], e13_[i], e12_[i]] for i in 1:length(ts)]

# mat = Material(Chaboche, ())
# mat.properties.youngs_modulus = 200.0e3
# mat.properties.poissons_ratio = 0.3
# mat.properties.yield_stress = 100.0
# mat.properties.K_n = 100.0
# mat.properties.n_n = 10.0
# mat.properties.C_1 = 10000.0
# mat.properties.D_1 = 100.0
# mat.properties.C_2 = 50000.0
# mat.properties.D_2 = 1000.0
# mat.properties.Q = 50.0
# mat.properties.b = 0.1
parameters = ChabocheParameterState(E = 200.0e3,
                                    nu = 0.3,
                                    R0 = 100.0,
                                    Kn = 100.0,
                                    nn = 10.0,
                                    C1 = 10000.0,
                                    D1 = 100.0,
                                    C2 = 50000.0,
                                    D2 = 1000.0,
                                    Q = 50.0,
                                    b = 0.1)
mat = Chaboche(parameters = parameters)

sim = Simulator(mat)
FEMMaterials.initialize!(sim, strains, ts)
t0 = time()
FEMMaterials.run!(sim)
t1 = time()
@info "Test took $(t1-t0)s!"
s33s = [s[3] for s in sim.stresses]
@test isapprox(s33s, s33_; rtol=0.001)


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

using BAT
using Test
using Distributions, LinearAlgebra, DensityInterface, ValueShapes

@testset "BAT_partitioned_sampling" begin
    σ = [0.1 0 0; 0 0.1 0; 0 0 0.1]
    μ = [1 1 1; -1 1 0; -1 -1 -1; 1 -1 0]
    mixture_model = MixtureModel(MvNormal[MvNormal(μ[i,:], Matrix(Hermitian(σ)) ) for i in 1:4])
    prior = NamedTupleDist(a = [Normal(0,2), Normal(0,2), Normal(0, 2)])
    likelihood = logfuncdensity(let model = mixture_model
        params -> logpdf(model, params.a)
    end)

    posterior = PosteriorDensity(likelihood, prior)
    transformed_posterior, trafo = bat_transform(PriorToUniform(), posterior)

    #Sampling and integration algorithms
    mcmc = MCMCSampling(mcalg = MetropolisHastings(), nsteps = 10^3);
    ahmi = ahmi = AHMIntegration(whitening = BAT.NoWhitening(), max_startingIDs = 10^3)
    mcmc_exp = MCMCSampling(mcalg = MetropolisHastings(), nsteps = 400, nchains=4, strict=false)

    ps = PartitionedSampling(sampler = mcmc, npartitions=4, exploration_sampler=mcmc_exp, integrator = ahmi, nmax_resampling=5);

    #Sampling with space partition
    results = bat_sample(posterior , ps)

    #Kolmogorov-Smirnov Test
    iid_distribution = NamedTupleDist(a = mixture_model,)
    iid_samples = bat_sample(iid_distribution, IIDSampling(nsamples=10^6))#Get iid samples
    ks_test = bat_compare(iid_samples.result, results.result)#Run Kolmogorov-Smirnov test
    @test all(ks_test.result.ks_p_values .> 0.7)#Check that all the p-values are bigger than 0.7
    @testset "Array of Posteriors" begin
        posteriors_array = BAT.convert_to_posterior(transformed_posterior, results.part_tree, extend_bounds = true)

        @test posteriors_array isa AbstractVector{<:BAT.AbstractPosteriorDensity}

        #Creates a matrix whose rows represents subspaces, cols represents modes and whose entries are 1s or 0s
        #1 if mean is in the corresponding subspace, 0 otherwise
        position_matrix = [Base.in(u, post.parbounds.vol) for u in trafo.(varshape(posterior).(eachrow(μ))), post in posteriors_array]
        #Check that each mode is in a different subspace
        @test sort([index[1] for index in findall(x -> x == true, position_matrix)]) == 1:4

    end
    @testset "Combinig Samples" begin
        #Sum of the weights should be equal to sum of integrals
        @test isapprox(sum(results.info.density_integral).val, sum(results.result.weight), rtol = 1e-2)
        @test size(results.info)[1] == 4#4 partitions
    end
end


text:
function loadgrid(filename::AbstractString)
    len = length(readline(filename))
    A = falses(len, len)
    for (y, line) in Iterators.enumerate(eachline(filename))
        for (x, c) in Iterators.enumerate(line)
            if c == '#'
                A[y, x] = true
            end
        end
    end
    A
end

function prettyprint(A)
    for row in eachrow(A)
        println((b ? '#' : '.' for b in row)...)
    end
end

function eachneighbour(A, index)
    (
        A[i+index]
        for i in CartesianIndices(((-1:1), (-1:1)))
        if !iszero(i) && checkbounds(Bool, A, i + index)
    )
end

function step(A, len = 1)
    for _ in 1:len
        A = [
            begin
                ns = sum(eachneighbour(A, i))
                ns == 3 || ns == 2 && A[i]
            end
            for i in CartesianIndices(A)
        ]
    end
    A
end

part1() = count(step(loadgrid("input.txt"), 100))

function step2(A, len = 1)
    h, w = size(A)
    A[1, 1] = true
    A[1, w] = true
    A[h, 1] = true
    A[h, w] = true
    for _ in 1:len
        A = [
            begin
                ns = sum(eachneighbour(A, i))
                ns == 3 || ns == 2 && A[i]
            end
            for i in CartesianIndices(A)
        ]
        A[1, 1] = true
        A[1, w] = true
        A[h, 1] = true
        A[h, w] = true
    end
    A
end

part2() = count(step2(loadgrid("input.txt"), 100))

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
using HDF5
using HDF5.Filters
using Test
using H5Zblosc, H5Zlz4, H5Zbzip2, H5Zzstd

using HDF5.Filters: UnknownFilter, isavailable, isencoderenabled, isdecoderenabled

@testset "filter" begin

# Create a new file
fn = tempname()

# Create test data
data = rand(1000, 1000)

# Open temp file for writing
f = h5open(fn, "w")

# Create datasets
dsdeflate = create_dataset(f, "deflate", datatype(data), dataspace(data),
                           chunk=(100, 100), deflate=3)

dsshufdef = create_dataset(f, "shufdef", datatype(data), dataspace(data),
                           chunk=(100, 100), shuffle=true, deflate=3)

dsfiltdef = create_dataset(f, "filtdef", datatype(data), dataspace(data),
                           chunk=(100, 100), filters=Filters.Deflate(3))

dsfiltshufdef = create_dataset(f, "filtshufdef", datatype(data), dataspace(data),
                               chunk=(100, 100), filters=[Filters.Shuffle(), Filters.Deflate(3)])

# Write data
write(dsdeflate, data)
write(dsshufdef, data)
write(dsfiltdef, data)
write(dsfiltshufdef, data)

# Test compression filters

compressionFilters = Dict(
    "blosc" => BloscFilter,
    "bzip2" => Bzip2Filter,
    "lz4" => Lz4Filter,
    "zstd" => ZstdFilter
)

for (name, filter) in compressionFilters

    ds = create_dataset(
        f, name, datatype(data), dataspace(data),
        chunk=(100,100), filters=filter()
    )
    write(ds, data)

    ds = create_dataset(
        f, "shuffle+"*name, datatype(data), dataspace(data),
        chunk=(100,100), filters=[Filters.Shuffle(), filter()]
    )
    write(ds, data)

end

ds = create_dataset(
    f, "blosc_bitshuffle", datatype(data), dataspace(data),
    chunk=(100,100), filters=BloscFilter(shuffle=H5Zblosc.BITSHUFFLE)
)
write(ds, data)

# Close and re-open file for reading
close(f)
f = h5open(fn)

# Read datasets and test for equality
for name in keys(f)
    ds = f[name]
    @testset "$name" begin
        @debug "Filter Dataset" HDF5.name(ds)
        @test ds[] == data
        filters = HDF5.get_create_properties(ds).filters
        if startswith(name, "shuffle+")
            @test filters[1] isa Shuffle
            @test filters[2] isa compressionFilters[name[9:end]]
        elseif haskey(compressionFilters, name) || name == "blosc_bitshuffle"
            name = replace(name, r"_.*"=>"")
            @test filters[1] isa compressionFilters[name]
        end
    end
end

close(f)

# Issue #896 and https://github.com/JuliaIO/HDF5.jl/issues/285#issuecomment-1002243321
# Create an UnknownFilter from a Tuple
h5open(fn, "w") do f
    data = rand(UInt8, 512, 16, 512)
    # Tuple of integers should become an Unknown Filter
    ds, dt = create_dataset(f, "data", data, chunk=(256,1,256), filter=(H5Z_FILTER_BZIP2, 0))
    # Tuple of Filters should get pushed into the pipeline one by one
    dsfiltshufdef = create_dataset(f, "filtshufdef", datatype(data), dataspace(data),
                               chunk=(128, 4, 128), filters=(Filters.Shuffle(), Filters.Deflate(3)))
    write(ds, data)
    close(ds)
    write(dsfiltshufdef, data)
    close(dsfiltshufdef)
end

h5open(fn, "r") do f
    @test f["data"][] == data
    @test f["filtshufdef"][] == data
end

# Filter Pipeline test for UnknownFilter
FILTERS_backup = copy(HDF5.Filters.FILTERS)
empty!(HDF5.Filters.FILTERS)
h5open(fn, "w") do f
    data = collect(1:128)
    filter = UnknownFilter(H5Z_FILTER_LZ4, 0, Cuint[0, 2, 4, 6, 8, 10], "Unknown LZ4", 0)
    ds, dt = create_dataset(f, "data", data, chunk=(32,), filters=filter)
    dcpl = HDF5.get_create_properties(ds)
    pipeline = HDF5.Filters.FilterPipeline(dcpl)
    @test pipeline[1].data == filter.data
end
merge!(HDF5.Filters.FILTERS, FILTERS_backup)

@test HDF5.API.h5z_filter_avail(HDF5.API.H5Z_FILTER_DEFLATE)
@test HDF5.API.h5z_filter_avail(HDF5.API.H5Z_FILTER_FLETCHER32)
@test HDF5.API.h5z_filter_avail(HDF5.API.H5Z_FILTER_NBIT)
@test HDF5.API.h5z_filter_avail(HDF5.API.H5Z_FILTER_SCALEOFFSET)
@test HDF5.API.h5z_filter_avail(HDF5.API.H5Z_FILTER_SHUFFLE)
@static if Sys.iswindows() || VERSION ≤ v"1.6"
    @test HDF5.API.h5z_filter_avail(HDF5.API.H5Z_FILTER_SZIP)
else
    # These are missing in the macOS and Linux JLLs
    @test_broken HDF5.API.h5z_filter_avail(HDF5.API.H5Z_FILTER_SZIP)
end
@test HDF5.API.h5z_filter_avail(H5Z_FILTER_BZIP2)
@test HDF5.API.h5z_filter_avail(H5Z_FILTER_LZ4)
@test HDF5.API.h5z_filter_avail(H5Z_FILTER_ZSTD)
@test HDF5.API.h5z_filter_avail(H5Z_FILTER_BLOSC)
HDF5.API.h5z_unregister(H5Z_FILTER_LZ4)
HDF5.Filters.register_filter(H5Zlz4.Lz4Filter)
@test isavailable(H5Z_FILTER_LZ4)
@test isavailable(Lz4Filter)
@test isencoderenabled(H5Z_FILTER_LZ4)
@test isdecoderenabled(H5Z_FILTER_LZ4)
@test isencoderenabled(Lz4Filter)
@test isdecoderenabled(Lz4Filter)

end # @testset "filter"


text:
#SESSION FILE
#Managing Session Attribute, workspace, usefull datas about the current work.
#---------------------------------------------------------------------


text:
using WriteVTK
using Documenter

DocMeta.setdocmeta!(
    WriteVTK, :DocTestSetup,
    quote
        using WriteVTK
        using StaticArrays  # not yet used...
    end;
    recursive = true,
)

makedocs(;
    modules = [WriteVTK],
    authors = "Juan Ignacio Polanco <jipolanc@gmail.com> and contributors",
    repo = "https://github.com/jipolanco/WriteVTK.jl/blob/{commit}{path}#L{line}",
    sitename = "WriteVTK.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://jipolanco.github.io/WriteVTK.jl",
        assets = [
            "assets/tomate.js",
        ],
        mathengine = KaTeX(),
    ),
    pages=[
        "Home" => "index.md",
        "Grid formats" => [
            "grids/syntax.md",
            "grids/structured.md",
            "grids/unstructured.md",
            "grids/datasets.md",
        ],
        "Metadata formats" => [
            "metadata/multiblock.md",
            "metadata/paraview_collections.md",
            "metadata/parallel.md",
        ],
        "Additional tools" => [
            "tools/write_array.md",
            "tools/readvtk.md",
        ],
        "API Reference" => "API.md",
    ],
)

deploydocs(;
    repo = "github.com/jipolanco/WriteVTK.jl.git",
    forcepush = true,
)


text:
module TestCopulas

using Copulas
using Base.Test

epsTol = 1e-12

###############################
## CDVine consistency checks ##
###############################

## load CDVine test results
##-------------------------

dataPath = joinpath(Pkg.dir("Copulas"),
                    "test/data/r_cop_funcs_results.csv")

testResults = readcsv(dataPath, header = true)[1]
resDict = [testResults[ii, 1] => testResults[ii, 2]
           for ii=1:size(testResults, 1)]

function tName(func, cop, copInd, ptInd)
    return string(func, "_", cop, "_", copInd, "_", ptInd)
end

## single function tests
##----------------------

cop = Copulas.ClaytonPC_Cpp([4.1])

u1, u2 = (0.3, 0.4)

testName = tName("pdf", "clay", 1, 1)
val = Copulas.pdf(cop, u1, u2)
expOut = resDict[testName]
@test_approx_eq_eps val expOut epsTol

testName = tName("cdf", "clay", 1, 1)
val = Copulas.cdf(cop, u1, u2)
expOut = resDict[testName]
@test_approx_eq_eps val expOut epsTol

testName = tName("hfun", "clay", 1, 1)
val = Copulas.hfun(cop, u1, u2)
expOut = resDict[testName]
@test_approx_eq_eps val expOut epsTol

testName = tName("vfun", "clay", 1, 1)
val = Copulas.vfun(cop, u2, u1)
expOut = resDict[testName]
@test_approx_eq_eps val expOut epsTol

## tests using metaprogramming
##----------------------------

testPoints = [(0.3, 0.4),
              (0.2, 0.8)]

funcs = [:pdf,
         :cdf,
         :hfun]

cops = ["gauss" Copulas.GaussianPC_Cpp([0.5]);
        "gauss" Copulas.GaussianPC_Cpp([-0.5]);
        "stud" Copulas.tPC_Cpp([0.5, 4.0]);
        "stud" Copulas.tPC_Cpp([-0.5, 5.0]);
        "clay" Copulas.ClaytonPC_Cpp([4.1]);
        "clay" Copulas.ClaytonPC_Cpp([6.3]);
        "gumb" Copulas.GumbelPC_Cpp([4.3]);
        "gumb" Copulas.GumbelPC_Cpp([1.2]);
        "frank" Copulas.FrankPC_Cpp([4.3]);
        "frank" Copulas.FrankPC_Cpp([1.2]);
        "joe" Copulas.JoePC_Cpp([4.3]);
        "joe" Copulas.JoePC_Cpp([1.2]);
        "bb1" Copulas.BB1PC_Cpp([0.3, 2.2]);
        "bb1" Copulas.BB1PC_Cpp([1.2, 2.8]);
        "bb6" Copulas.BB6PC_Cpp([4.3, 2.2]);
        "bb6" Copulas.BB6PC_Cpp([1.2, 1.1]);
        "bb7" Copulas.BB7PC_Cpp([1.9, 0.3]);
        "bb7" Copulas.BB7PC_Cpp([1.2, 8.2]);
        "bb8" Copulas.BB8PC_Cpp([4.3, 0.2]);
        "bb8" Copulas.BB8PC_Cpp([1.2, 0.4])]

for ptInd=1:length(testPoints)
    for copInd=1:size(cops, 1)
        for f in funcs
            copName = cops[copInd, 1]
            cop = cops[copInd, 2]
            u1, u2 = testPoints[ptInd]
            mod(copInd, 2) == 0 ? ind = 2 : ind = 1
            testName = tName(string(f), copName, ind, ptInd)
#            display(testName)
            eval(:(val = (Copulas.$f)(cop, u1, u2)))
            expOut = resDict[testName]
            @test_approx_eq_eps val expOut epsTol
        end
    end
end

## v functions
for ptInd=1:length(testPoints)
    for copInd=1:size(cops, 1)
        copName = cops[copInd, 1]
        cop = cops[copInd, 2]
        u1, u2 = testPoints[ptInd]
        mod(copInd, 2) == 0 ? ind = 2 : ind = 1
        testName = tName("vfun", copName, ind, ptInd)
        #            display(testName)
        val = Copulas.vfun(cop, u2, u1)
        expOut = resDict[testName]
        @test_approx_eq_eps val expOut epsTol
    end
end

########################
## consistency checks ##
########################

epsTol = 1e-7

uLow = 0.001
uHigh = 0.999
testGrid = [uLow, 0.1:0.1:0.9, uHigh]

cops = [
        Copulas.IndepPC_Cpp([0.3]),
        Copulas.IndepPC_Cpp([0.3]),
        Copulas.AMHPC_Cpp([0.4]),
        Copulas.AMHPC_Cpp([-0.4]),
        Copulas.AsymFGMPC_Cpp([0.2]),
        Copulas.AsymFGMPC_Cpp([0.9]),
        Copulas.BB1PC_Cpp([0.2, 5.4]),
        Copulas.BB1PC_Cpp([5.4, 1.2]),
        Copulas.BB6PC_Cpp([1.3, 4.3]),
        ## Copulas.BB6PC_Cpp([5.4, 0.2]),
        Copulas.BB7PC_Cpp([4.5, 4.4]),
        Copulas.BB7PC_Cpp([1.5, 4.2]),
        Copulas.BB8PC_Cpp([1.2, 0.2]),
        Copulas.BB8PC_Cpp([4.5, 0.8]),
        Copulas.ClaytonPC_Cpp([4.3]),
        Copulas.ClaytonPC_Cpp([2]),
        Copulas.FGMPC_Cpp([0.4]),
        Copulas.FGMPC_Cpp([-0.7]),
        Copulas.FrankPC_Cpp([-20]),
        Copulas.FrankPC_Cpp([22]),
        Copulas.GaussianPC_Cpp([0.6]),
        Copulas.GaussianPC_Cpp([-0.2]),
        Copulas.GumbelPC_Cpp([1.2]),
        Copulas.GumbelPC_Cpp([4.3]),
        Copulas.IteratedFGMPC_Cpp([0.5, 0.5]),
        Copulas.IteratedFGMPC_Cpp([-0.5, -0.4]),
        Copulas.JoePC_Cpp([4.2]),
        Copulas.JoePC_Cpp([8]),
        Copulas.PartialFrankPC_Cpp([1.2]),
        Copulas.PartialFrankPC_Cpp([12.5]),
        Copulas.PlackettPC_Cpp([1.2]),
        Copulas.PlackettPC_Cpp([18]),
        Copulas.Tawn1PC_Cpp([1.2, 0.2]),
        Copulas.Tawn1PC_Cpp([4.3, 0.4]),
        Copulas.Tawn2PC_Cpp([5.2, 0.8]),
        Copulas.Tawn2PC_Cpp([1.2, 0.4 ]),
        Copulas.TawnPC_Cpp([4.4, 0.3, 0.8]),
        Copulas.TawnPC_Cpp([12.8, 0.4, 0.5]),
        Copulas.tPC_Cpp([0.5, 4.0]),
        Copulas.tPC_Cpp([-0.3, 8.0])
        ]

## pdf, cdf, h- and v-functions
##-----------------------------

for cop in cops
    show(cop)
    println("")
    println(" *  cdf tests:")
    Copulas.testCdf(cop, epsTol, testGrid)
    println(" *  h/v tests:")
    Copulas.testHandV(cop, epsTol, testGrid)
    println(" *  pdf tests:")
    Copulas.testPdf(cop, testGrid)
end

## inverse h- and v-functions
##---------------------------

testPoints = [(0.3, 0.7),
              (0.00002, 0.00098),
              (0.1, 0.998),
              ## (0.8, 0.2132323543),
              (0.8, 0.2132),
              ## (0.2349248, 0.983475),
              (0.999, 0.999),
              (0.0001, 0.0001)]

println("-------------------------------")
println("inverse h- and v-function tests")
println("-------------------------------")

for pts in testPoints
    u1, u2 = pts
    for cop in cops
        copShown = false

        q = Copulas.hfun(cop, u1, u2)
        u1_hat = Copulas.hinv(cop, q, [u2])
        if abs(u1 - u1_hat)[1] > epsTol
            # show copula
            if !copShown
                show(cop)
                println("")
                copShown = true
            end

            # at least both values should be mapped to same hfun value
            println("u2=$u2: $u1 -> $u1_hat")
            if !(q == Copulas.hfun(cop, u1_hat, u2))
                println("different h-function values")
                h1 = Copulas.hfun(cop, u1, u2)
                h2 = Copulas.hfun(cop, u1_hat, u2)
                println("  h($u1) = $h1")
                println("  h($u1_hat) = $h2")
            end
        end

        q = Copulas.vfun(cop, u1, u2)
        u1_hat = Copulas.vinv(cop, q, [u2])
        if abs(u1 - u1_hat)[1] > epsTol
            # show copula
            if !copShown
                show(cop)
                println("")
                copShown = true
            end

            # at least both values should be mapped to same hfun value
            println("u2=$u2: $u1 -> $u1_hat")
            if !(q == Copulas.vfun(cop, u1_hat, u2))
                println("different v-function values")
                h1 = Copulas.vfun(cop, u1, u2)
                h2 = Copulas.vfun(cop, u1_hat, u2)
                println("  v($u1) = $h1")
                println("  v($u1_hat) = $h2")
            end
        end

        ## q = Copulas.hinv(cop, u1, u2)
        ## u1_hat = Copulas.hfun(cop, q, [u2])
        ## hdiff2 = u1 - u1_hat
        ## ## @test_approx_eq_eps u1 u1_hat epsTol
        ## q = Copulas.vfun(cop, u2, u1)
        ## u2_hat = Copulas.vinv(cop, q, [u1])
        ## vdiff1 = u2 - u2_hat
        ## q = Copulas.vfun(cop, u2, u1)
        ## u2_hat = Copulas.vinv(cop, q, [u1])
        ## vdiff2 = u2 - u2_hat
        ## ## @test_approx_eq_eps u2 u2_hat epsTol
        ## maxDiff = max(abs(hdiff1[1]),
        ##               abs(hdiff2[1]),
        ##               abs(vdiff1[1]),
        ##               abs(vdiff2[1]))
        ## if abs(maxDiff) > 1e-10
        ##     show(cop)
        ##     println("")
        ##     println("diffh1: $hdiff1")
        ##     println("diffh2: $hdiff2")
        ##     println("diffv1: $vdiff1")
        ##     println("diffv2: $vdiff2")
        ## end
    end
end

## ID functions
##-------------

cop = Copulas.GaussianPC_Cpp([0.5])
@test 10 == Copulas.getVineCppId(cop)
@test 2 == Copulas.getVineCppId(:AsymFGM)
@test 0 == Copulas.getVineCppId(:Indep)

@test :Gumbel == Copulas.getCopNam(11)

@test :IteratedFGMPC_Cpp == Copulas.getCopType_Cpp(12)

## display functions
##------------------

cop = Copulas.GaussianPC_Cpp([0.5])
@test "GaussianPC_Cpp: 0.5" == Copulas.displayShortString(cop)

cop = Copulas.tPC_Cpp([0.5, 12.])
@test "tPC_Cpp: 0.5, 12.0" == Copulas.displayShortString(cop)

display(cop)

## equality
##---------

cop1 = Copulas.GaussianPC_Cpp([0.5])
cop2 = Copulas.GaussianPC_Cpp([0.5])
@test cop1 == cop2

end


text:
function _promote(args...)
    _type = Base.promote_eltype(args...)
    return map(x->convert.(_type, x), args)
end

## Utilities

"""
$(SIGNATURES)

Check if the problem has control inputs.
"""
is_autonomous(::AbstractDataDrivenProblem{N, U, C}) where {N,U,C} = U

"""
$(SIGNATURES)

Check if the problem is time discrete.
"""
is_discrete(::AbstractDataDrivenProblem{N, U, C}) where {N,U,C} = C == DDProbType(2)

"""
$(SIGNATURES)

Check if the problem is direct.
"""
is_direct(::AbstractDataDrivenProblem{N, U, C}) where {N,U,C} = C == DDProbType(1)

"""
$(SIGNATURES)

Check if the problem is time continuous.
"""
is_continuous(::AbstractDataDrivenProblem{N, U, C}) where {N,U,C} = C == DDProbType(3)

"""
$(SIGNATURES)

Check if the problem is parameterized.
"""
is_parametrized(x::AbstractDataDrivenProblem{N, U, C}) where {N,U,C} = hasfield(typeof(x), :p) && !isempty(x.p)

"""
$(SIGNATURES)

Check if the problem has associated measurement times.
"""
has_timepoints(x::AbstractDataDrivenProblem{N,U,C}) where {N,U,C} = !isempty(x.t)

## Concrete Type
"""
$(TYPEDEF)

The `DataDrivenProblem` defines a general estimation problem given measurements, inputs and (in the near future) observations.
Two construction methods are available:

+ `DiscreteDataDrivenProblem` for time discrete systems
+ `ContinousDataDrivenProblem` for systems continouos in time

both are aliases for constructing a problem.

# Fields
$(FIELDS)

# Signatures
$(SIGNATURES)

# Example

```julia
X, DX, t = data...

# Define a discrete time problem
prob = DiscreteDataDrivenProblem(X)

# Define a continous time problem without explicit time points
prob = ContinuousDataDrivenProblem(X, DX)

# Define a continous time problem without explict derivatives
prob = ContinuousDataDrivenProblem(X, t)

# Define a discrete time problem with an input function as a function
input_signal(u,p,t) = t^2
prob = DiscreteDataDrivenProblem(X, t, input_signal)
```
"""
struct DataDrivenProblem{dType, cType, probType} <: AbstractDataDrivenProblem{dType, cType, probType}

    # Data
    """State measurements"""
    X::AbstractMatrix{dType}
    """Time measurements (optional)"""
    t::AbstractVector{dType}

    """Differental state measurements (optional); Used for time continuous problems"""
    DX::AbstractMatrix{dType}
    """Output measurements (optional); Used for direct problems"""
    Y::AbstractMatrix{dType}
    """Input measurements (optional); Used for non-autonoumous problems"""
    U::AbstractMatrix{dType}

    """Parameters associated with the problem (optional)"""
    p::AbstractVector{dType}
end

function DataDrivenProblem(X, t, DX, Y, U, p)
    dType = Base.promote_eltype(X, t, DX, Y, U, p)
    cType = isempty(U)
    # We assume a discrete Problem
    probType = DDProbType(2)
    if (isempty(DX) && !isempty(Y))
        probType = DDProbType(1) # Direct problem
    elseif !isempty(DX)
        probType = DDProbType(3) # Continouos
    end
    return DataDrivenProblem{dType, cType, probType}(_promote(X,t,DX,Y,U,p)...)
end

function DataDrivenProblem(X, t, DX, Y, U::F, p) where F <: Function
    # Generate the input as a Matrix

    ts = isempty(t) ? zeros(eltype(X), size(X,2)) : t

    u_ = hcat(map(i->U(X[:,i], p, ts[i]), 1:size(X,2))...)

    return DataDrivenProblem(_promote(X,t,DX,Y,u_,p)...)
end

function DataDrivenProblem(X::AbstractMatrix;
    t::AbstractVector = zeros(eltype(X), size(X,2)),
    DX::AbstractMatrix = Array{eltype(X)}(undef, 0, 0),
    Y::AbstractMatrix = Array{eltype(X)}(undef, 0,0),
    U::F = Array{eltype(X)}(undef, 0,0),
    p::AbstractVector = Array{eltype(X)}(undef, 0)
    ) where F <: Union{AbstractMatrix, Function}

    return DataDrivenProblem(X,t,DX,Y,U,p)
end

function Base.summary(io::IO, x::DataDrivenProblem{N,C,P}) where {N,C,P}
    print(io, "$P DataDrivenProblem{$N}")
    C ? nothing : print(io, " with controls")
    return
end

function Base.print(io, x::DataDrivenProblem{N,C,P}) where {N,C,P}
    summary(io, x)
end

Base.show(io::IO, x::DataDrivenProblem{N,C,P}) where {N,C,P} = summary(io, x)

## Discrete Constructors
"""
A time discrete `DataDrivenProblem` useable for problems of the form `f(x[i],p,t,u) ↦ x[i+1]`.

$(SIGNATURES)
"""
function DiscreteDataDrivenProblem(X::AbstractMatrix; kwargs...)
    DataDrivenProblem(X; kwargs...)
end

function DiscreteDataDrivenProblem(X::AbstractMatrix, t::AbstractVector; kwargs...)
    DataDrivenProblem(X, t=t ; kwargs...)
end

function DiscreteDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, U::AbstractMatrix; kwargs...)
    return DataDrivenProblem(X, t=t, U = U ; kwargs...)
end

function DiscreteDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, U::Function; kwargs...)
    return DataDrivenProblem(X, t=t, U = U ; kwargs...)
end

## Continouos Constructors
"""
A time continuous `DataDrivenProblem` useable for problems of the form `f(x,p,t,u) ↦ dx/dt`.

$(SIGNATURES)

Automatically constructs derivatives via an additional collocation method, which can be either a collocation
or an interpolation from `DataInterpolations.jl` wrapped by an `InterpolationMethod`.
"""
function ContinuousDataDrivenProblem(X::AbstractMatrix, DX::AbstractMatrix; kwargs...)
    return DataDrivenProblem(X, DX = DX; kwargs...)
end

function ContinuousDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, DX::AbstractMatrix; kwargs...)
    return DataDrivenProblem(X, t = t, DX = DX; kwargs...)
end

function ContinuousDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, DX::AbstractMatrix, U::AbstractMatrix; kwargs...)
    return DataDrivenProblem(X, t = t, DX = DX, U = U; kwargs...)
end

function ContinuousDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, DX::AbstractMatrix, U::F; kwargs...) where {F <: Function}
    return DataDrivenProblem(X, t = t, DX = DX, U = U; kwargs...)
end

function ContinuousDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, collocation = InterpolationMethod(); kwargs...)
    dx, x = collocate_data(X, t, collocation)
    return DataDrivenProblem(x, t = t, DX = dx; kwargs...)
end

function ContinuousDataDrivenProblem(X::AbstractMatrix, t::AbstractVector,  U::AbstractMatrix, collocation; kwargs...)
    dx, x = collocate_data(X, t, collocation)
    return DataDrivenProblem(x, t = t, DX = dx, U = U; kwargs...)
end

## Direct Constructors
"""
A direct `DataDrivenProblem` useable for problems of the form `f(x,p,t,u) ↦ y`.

$(SIGNATURES)
"""
function DirectDataDrivenProblem(X::AbstractMatrix, Y::AbstractMatrix; kwargs...)
    return DataDrivenProblem(X, Y = Y; kwargs...)
end

function DirectDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, Y::AbstractMatrix; kwargs...)
    return DataDrivenProblem(X, t = t, Y = Y; kwargs...)
end

function DirectDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, Y::AbstractMatrix, U; kwargs...)
    return DataDrivenProblem(X, t = t, Y = Y, U = U; kwargs...)
end

## Utils

# Check for nans, inf etc
check_domain(x) =  @assert all(.~isnan.(x)) && all(.~isinf.(x)) ("One or more measurements contain `NaN` or `Inf`.")
check_lengths(args...) = @assert all(map(x->size(x)[end], args) .== size(first(args))[end]) "One or more measurements are not sized equally."

# Return the target variables
get_target(x::AbstractDirectProb{N,C}) where {N,C} = x.Y
get_target(x::AbstractDiscreteProb{N,C}) where {N,C} = x.X[:,2:end]
get_target(x::AbstracContProb{N,C}) where {N,C} = x.DX

get_oop_args(x::AbstractDataDrivenProblem{N,C,P}) where {N <: Number, C, P} = map(f->getfield(x, f), (:X, :p, :t, :U))

function get_oop_args(x::AbstractDiscreteProb{N,C}) where {N <: Number, C}
    return  (
        x.X[:, 1:end-1],
        x.p,
        x.t[1:end-1],
        x.U[:, 1:end-1]
    )
end

function get_implicit_oop_args(x::AbstractDirectProb{N,C}) where {N <: Number, C}
    return  (
        [x.X; x.Y],
        x.p,
        x.t,
        x.U
    )
end

function get_implicit_oop_args(x::AbstracContProb{N,C}) where {N <: Number, C}
    return  (
        [x.X; x.DX],
        x.p,
        x.t,
        x.U
    )
end

function get_implicit_oop_args(x::AbstractDiscreteProb{N,C}) where {N <: Number, C}
    return  (
        [x.X[:, 1:end-1]; x.X[:, 2:end]],
        x.p,
        x.t[1:end-1],
        x.U[:, 1:end-1]
    )
end

"""
$(SIGNATURES)

Checks if a `DataDrivenProblem` is valid by checking if the data contains `NaN`, `Inf` and
if the number of measurements is consistent.

# Example

```julia
is_valid(problem)
```
"""
function is_valid(x::AbstractDirectProb{N,C}) where {N <: Number,C}
    map(check_domain, (x.X, x.Y, x.U, x.t, x.p))
    if !C && !isempty(x.t)
        check_lengths(x.X, x.Y, x.U, x.t)
    elseif !C
        check_lengths(x.X, x.Y, x.U)
    elseif C && !isempty(x.t)
        check_lengths(x.X, x.Y, x.t)
    else
        check_lengths(x.X, x.Y)
    end
    return true
end

function is_valid(x::AbstractDiscreteProb{N,C}) where {N <: Number,C}
    map(check_domain, (x.X, x.U, x.t, x.p))
    if !C && !isempty(x.t)
        check_lengths(x.X, x.U, x.t)
    elseif !C
        check_lengths(x.X, x.U)
    elseif C && !isempty(x.t)
        check_lengths(x.X, x.t)
    else
        check_lengths(x.X)
    end
    return true
end

function is_valid(x::AbstracContProb{N,C}) where {N <: Number,C}
    map(check_domain, (x.X, x.DX, x.U, x.t, x.p))
    if !C && !isempty(x.t)
        check_lengths(x.X, x.DX, x.U, x.t)
    elseif !C
        check_lengths(x.X, x.DX, x.U)
    elseif C && !isempty(x.t)
        check_lengths(x.X, x.DX, x.t)
    else
        check_lengths(x.X, x.DX)
    end
    return true
end

## DESolution dispatch

ContinuousDataDrivenProblem(sol::T; kwargs...) where T <: DiffEqBase.DESolution = DataDrivenProblem(sol; kwargs...)
DiscreteDataDrivenProblem(sol::T; kwargs...) where T <: DiffEqBase.DESolution = DataDrivenProblem(sol; kwargs...)

function DataDrivenProblem(sol::T; use_interpolation = false, kwargs...) where T <: DiffEqBase.DESolution
    if sol.retcode != :Success
        throw(AssertionError("The solution is not successful. Abort."))
        return
    end

    X = Array(sol)
    t = sol.t
    p = sol.prob.p

    p = isa(p, DiffEqBase.NullParameters) ? eltype(X)[] : p

    if isdiscrete(sol.alg)

        return DiscreteDataDrivenProblem(
            X, t, p = p; kwargs...
        )

    else
        if use_interpolation
            DX = Array(sol(sol.t, Val{1}))
        else
            DX = similar(X)
            if DiffEqBase.isinplace(sol.prob.f)
                @views map(i->sol.prob.f(DX[:, i], X[:, i], p, t[i]), 1:size(X,2))
            else
                map(i->DX[:, i] .= sol.prob.f(X[:, i], p, t[i]), 1:size(X,2))
            end
        end

        return ContinuousDataDrivenProblem(
            X, t, DX = DX, p = p; kwargs...
        )
    end

end


text:
module ResortRoofEnding

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Resort Roof Ending" => Ahorn.EntityPlacement(
        Maple.ResortRoofEnding,
        "rectangle"
    )
)

Ahorn.minimumSize(entity::Maple.ResortRoofEnding) = 8, 0
Ahorn.resizable(entity::Maple.ResortRoofEnding) = true, false

function Ahorn.selection(entity::Maple.ResortRoofEnding)
    x, y = Ahorn.position(entity)
    width = get(entity.data, "width", 8)

    return Ahorn.Rectangle(x, y, ceil(Int, width / 16) * 16 + 16, 16)
end

const centerTextures = String[
    "decals/3-resort/roofCenter.png",
    "decals/3-resort/roofCenter_b.png",
    "decals/3-resort/roofCenter_c.png",
    "decals/3-resort/roofCenter_d.png",
]

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.ResortRoofEnding, room::Maple.Room)
    rng = Ahorn.getSimpleEntityRng(entity)

    width = get(entity.data, "width", 8)
    Ahorn.drawSprite(ctx, "decals/3-resort/roofEdge_d.png", 8, 4)

    offset = 0
    while offset < width
        Ahorn.drawSprite(ctx, rand(rng, centerTextures), offset + 8, 4)

        offset += 16
    end

    Ahorn.drawSprite(ctx, "decals/3-resort/roofEdge.png", offset + 8, 4)
end

end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
#------------------------------------------------------------------------------#
#                                   Includes                                   #
#------------------------------------------------------------------------------#

using Test
using IsoOrthoTensor

#------------------------------------------------------------------------------#
#                            Isotropic Tensor Tests                            #
#------------------------------------------------------------------------------#

@testset "Isotropic Tensor Tests                            " begin
    @test Δ(1) == K(2)
    @test Δ(2) == 𝕔((K(2), K(2)), (1,))
    @test Δ(3) == 𝕔((K(2), Δ(2)), (1,))
    @test Δ(4) == 𝕔((K(2), Δ(3)), (1,))
    @test Δ(5) == 𝕔((K(2), Δ(4)), (1,))
end

#------------------------------------------------------------------------------#
#                          Orthogonality Tensor Tests                          #
#------------------------------------------------------------------------------#

@testset "Orthogonality Tensor Tests                        " begin
    @test Ο(1) == K(2)
    @test Ο(2) == 𝕡((K(2), K(2)), (1, 3))
    @test Ο(3) == 𝕡((K(2), K(2), K(2)), (1, 3, 5))
    @test Ο(4) == 𝕡((K(2), K(2), K(2), K(2)), (1, 3, 5, 7))
    @test Ο(5) == 𝕡((K(2), K(2), K(2), K(2), K(2)), (1, 3, 5, 7, 9))
end


text:
using Graphs
using Base.Test

g = simple_complete_graph(5)
@test num_vertices(g) == 5 && num_edges(g) == 20
g = simple_complete_graph(5, is_directed=false)
@test num_vertices(g) == 5 && num_edges(g) == 10

g = simple_star_graph(5)
@test num_vertices(g) == 5 && num_edges(g) == 4
g = simple_star_graph(5, is_directed=false)
@test num_vertices(g) == 5 && num_edges(g) == 4
g = simple_star_graph(1)
@test num_vertices(g) == 1 && num_edges(g) == 0

g = simple_path_graph(5)
@test num_vertices(g) == 5 && num_edges(g) == 4
g = simple_path_graph(5; is_directed=false)
@test num_vertices(g) == 5 && num_edges(g) == 4

g = simple_wheel_graph(5)
@test num_vertices(g) == 5 && num_edges(g) == 8
g = simple_wheel_graph(5, is_directed=false)
@test num_vertices(g) == 5 && num_edges(g) == 8

g = simple_diamond_graph()
@test num_vertices(g) == 4 && num_edges(g) == 5

g = simple_bull_graph()
@test num_vertices(g) == 5 && num_edges(g) == 5

g = simple_chvatal_graph()
@test num_vertices(g) == 12 && num_edges(g) == 24

g = simple_cubical_graph()
@test num_vertices(g) == 8 && num_edges(g) == 12

g = simple_desargues_graph()
@test num_vertices(g) == 20 && num_edges(g) == 30

g = simple_dodecahedral_graph()
@test num_vertices(g) == 20 && num_edges(g) == 30

g = simple_frucht_graph()
@test num_vertices(g) == 20 && num_edges(g) == 18

g = simple_heawood_graph()
@test num_vertices(g) == 14 && num_edges(g) == 21

g = simple_house_graph()
@test num_vertices(g) == 5 && num_edges(g) == 6

g = simple_house_x_graph()
@test num_vertices(g) == 5 && num_edges(g) == 8

g = simple_icosahedral_graph()
@test num_vertices(g) == 12 && num_edges(g) == 30

g = simple_krackhardt_kite_graph()
@test num_vertices(g) == 10 && num_edges(g) == 18

g = moebius_kantor_graph()
@test num_vertices(g) == 16 && num_edges(g) == 24

g = simple_octahedral_graph()
@test num_vertices(g) == 6 && num_edges(g) == 12

g = simple_pappus_graph()
@test num_vertices(g) == 18 && num_edges(g) == 27

g = simple_petersen_graph()
@test num_vertices(g) == 10 && num_edges(g) == 15

g = simple_sedgewick_maze_graph()
@test num_vertices(g) == 8 && num_edges(g) == 10

g = simple_tetrahedral_graph()
@test num_vertices(g) == 4 && num_edges(g) == 6

g = simple_truncated_cube_graph()
@test num_vertices(g) == 24 && num_edges(g) == 36

g = simple_truncated_tetrahedron_graph()
@test num_vertices(g) == 12 && num_edges(g) == 18

g = simple_tutte_graph()
@test num_vertices(g) == 46 && num_edges(g) == 69


text:
using Documenter
using Latexify
using Markdown
using Markdown: MD, Paragraph, LineBreak
using Plots

using GeometricProblems.Diagnostics
using GeometricProblems.PlotRecipes

using GeometricIntegrators
using GeometricIntegrators.Integrators.SPARK
using GeometricIntegrators.Integrators.SPARK: get_ã_vspark_primary,
                                                get_α_vspark_primary,
                                                compute_ã_vspark_primary,
                                                compute_α_vspark_primary

function Integrators.integrate!(int::Union{IntegratorSPARK{DT,TT,D,S,R},IntegratorVSPARKprimary{DT,TT,D,S,R}}, sol::Solution{DT,TT}) where {DT,TT,D,S,R}
    asol = AtomicSolution(int)

    Qi = [SDataSeries(DT, D, ntime(sol), nsamples(sol)) for i in 1:S]
    Pi = [SDataSeries(DT, D, ntime(sol), nsamples(sol)) for i in 1:S]
    Vi = [SDataSeries(DT, D, ntime(sol), nsamples(sol)) for i in 1:S]
    Φi = [SDataSeries(DT, D, ntime(sol), nsamples(sol)) for i in 1:S]

    Qp = [SDataSeries(DT, D, ntime(sol), nsamples(sol)) for i in 1:R]
    Pp = [SDataSeries(DT, D, ntime(sol), nsamples(sol)) for i in 1:R]
    Λp = [SDataSeries(DT, D, ntime(sol), nsamples(sol)) for i in 1:R]
    Φp = [SDataSeries(DT, D, ntime(sol), nsamples(sol)) for i in 1:R]

    # loop over initial conditions showing progress bar
    for m in eachsample(sol)
        # get cache from solution
        get_initial_conditions!(sol, asol, m, 1)
        initialize!(int, asol)

        # loop over time steps
        n = 0
        try
            for outer n in eachtimestep(sol)
                integrate!(int, sol, asol, m, n)

                for i in eachindex(Qi,Pi,Vi,Φi)
                    set_data!(Qi[i], asol.internal.Qi[i], n, m)
                    set_data!(Pi[i], asol.internal.Pi[i], n, m)
                    set_data!(Vi[i], asol.internal.Vi[i], n, m)
                    set_data!(Φi[i], asol.internal.Φi[i], n, m)
                end

                for i in eachindex(Qp,Pp,Λp,Φp)
                    set_data!(Qp[i], asol.internal.Qp[i], n, m)
                    set_data!(Pp[i], asol.internal.Pp[i], n, m)
                    set_data!(Λp[i], asol.internal.Λp[i], n, m)
                    set_data!(Φp[i], asol.internal.Φp[i], n, m)
                end
            end
        catch ex
            tstr = " in time step " * string(n)

            if nsamples(sol) > 1
                tstr *= " for initial condition " * string(m)
            end

            tstr *= "."

            if isa(ex, DomainError)
                show(stdout, "text/markdown", Markdown.parse("DOMAIN ERROR: Simulation crashed" * tstr))
                @warn("DOMAIN ERROR: Simulation crashed" * tstr)
            elseif isa(ex, NonlinearSolverException)
                show(stdout, "text/markdown", Markdown.parse("SOLVER ERROR: Simulation crashed" * tstr))
                show(stdout, "text/markdown", Markdown.parse(ex.msg))
                @warn("SOLVER ERROR: Simulation crashed" * tstr)
                @warn(ex.msg)
            else
                show(stdout, "text/markdown", Markdown.parse("ERROR: Simulation crashed with $(typeof(ex)) " * tstr))
                show(stdout, "text/markdown", Markdown.parse(ex.msg))
                @warn("ERROR: Simulation crashed with $(typeof(ex)) " * tstr)
                @warn(ex.msg)
                # showerror(stdout, ex, catch_backtrace())
            end
        end
    end

    (Qi=Qi, Pi=Pi, Vi=Vi, Φi=Φi, Qp=Qp, Pp=Pp, Λp=Λp, Φp=Φp)
end

# function _arr_str(a)
#     # io = IOBuffer()
#     # show(io, "text/plain", a)
#     # String(take!(io))

#     # astr = latexify(a; env=:array)
#     # astr = replace(astr, "\\begin{equation}" => "\$")
#     # astr = replace(astr, "\\end{equation}" => "\$")
#     # Markdown.parse(astr)
# end

_arr_str(a) = latexify(a; env=:mdtable, latex=false, side=collect(axes(a,1)), head=collect(axes(a,2)))
_linebreak(io) = show(io, "text/markdown", MD(Paragraph([LineBreak()])))

function write_symplecticity(tab, dir, file, name)
    if !isdir(dir)
        mkdir(dir)
    end

    symp_file = dir * "/" * file * ".md"

    symp_cond = symplecticity_conditions(tab)
    symp_arrs = check_symplecticity(tab)

    open(symp_file, "w") do f
        show(f, "text/markdown", Markdown.parse("# $name"))
        _linebreak(f)
        show(f, "text/markdown", Markdown.parse("## Symplecticity Conditions"))
        _linebreak(f)

        for i in eachindex(symp_arrs, symp_cond)
            show(f, "text/markdown", Markdown.parse(symp_cond[i]))
            _linebreak(f)
            show(f, "text/markdown", _arr_str(symp_arrs[i]))
            _linebreak(f)
        end
    end

end

function _plot_figure_md(file, name, filename)
    # if isfile(filename)
        show(file, "text/markdown", Markdown.parse("![$name]($filename)"))
        _linebreak(file)
    # else
    #     show(stdout, "text/markdown", Markdown.parse("ERROR: Plot output $filename does not exist!"))
    #     @warn("Plot output $filename does not exist!")
    # end
end

function write_plots(tab, dir, file, name, fig_suff)

    plot_file = file * ".md"

    open(plot_file, "w") do f
        show(f, "text/markdown", Markdown.parse("# $name"))
        _linebreak(f)

        _plot_figure_md(f, name, "$(dir)/$(file)_solution$(fig_suff)")
        _plot_figure_md(f, name, "$(dir)/$(file)_traces$(fig_suff)")

        show(f, "text/markdown", Markdown.parse("## Energy Error"))
        _linebreak(f)

        _plot_figure_md(f, name, "$(dir)/$(file)_energy_error$(fig_suff)")
        _plot_figure_md(f, name, "$(dir)/$(file)_energy_drift$(fig_suff)")

        show(f, "text/markdown", Markdown.parse("## Constraint"))
        _linebreak(f)

        _plot_figure_md(f, name, "$(dir)/$(file)_constraint_error$(fig_suff)")

        for i in 1:tab.s
            _plot_figure_md(f, name, "$(dir)/$(file)_constraint_error_phi_i$(i)$(fig_suff)")
        end

        for i in 1:tab.r
            _plot_figure_md(f, name, "$(dir)/$(file)_constraint_error_phi_p$(i)$(fig_suff)")
        end

        _plot_figure_md(f, name, "$(dir)/$(file)_lambda$(fig_suff)")

        for i in 1:tab.r
            _plot_figure_md(f, name, "$(dir)/$(file)_lambda_p$(i)$(fig_suff)")
        end
    end

end

function _plot(sol, stages, equ, dir, file, fig_suff)
    H, ΔH = compute_energy_error(sol.t, sol.q, equ.parameters)
    plotenergyerror(sol.t, ΔH; nt=lastentry(sol), fmt=:png)
    savefig(dir * "/" * file * "_energy_error" * fig_suff)

    plotenergydrift(compute_error_drift(sol.t, ΔH, div(ntime(sol), 10))...; nt=lastentry(sol))
    savefig(dir * "/" * file * "_energy_drift" * fig_suff)

    plotconstrainterror(sol.t, compute_momentum_error(sol.t, sol.q, sol.p, (t,q,k)->ϑ(t,q,equ.parameters,k)); nt=lastentry(sol), fmt=:png)
    savefig(dir * "/" * file * "_constraint_error" * fig_suff)

    plotlagrangemultiplier(sol.t, sol.λ; nt=lastentry(sol), fmt=:png)
    savefig(dir * "/" * file * "_lambda" * fig_suff)

    if stages !== nothing
        for i in eachindex(stages.Φi)
            plotconstrainterror(sol.t, stages.Φi[i]; nt=lastentry(sol), fmt=:png, plot_title="Φi,$(i)")
            savefig(dir * "/" * file * "_constraint_error_phi_i$(i)" * fig_suff)
        end

        for i in eachindex(stages.Φp)
            plotconstrainterror(sol.t, stages.Φp[i]; nt=lastentry(sol), fmt=:png, plot_title="Φp,$(i)")
            savefig(dir * "/" * file * "_constraint_error_phi_p$(i)" * fig_suff)
        end

        for i in eachindex(stages.Λp)
            plotlagrangemultiplier(sol.t, stages.Λp[i]; nt=lastentry(sol), fmt=:png, plot_title="Λp,$(i)")
            savefig(dir * "/" * file * "_lambda_p$(i)" * fig_suff)
        end
    end
end

function run_spark(idae, tab, nt, dir, file, fig_suff, phi_average)
    sol = Solution(idae, Δt, nt)
    int = Integrator(idae, tab, Δt)

    stages = integrate!(int, sol)

    if phi_average !== nothing
        push!(stages.Φp, SDataSeries(phi_average([parent(stages.Φp[i]) for i in eachindex(stages.Φp)])))
    end

    # if length(stages.Φp) == 2
    #     push!(stages.Φp, SDataSeries(parent(stages.Φp[1]) .+ parent(stages.Φp[2])))
    #     push!(stages.Φp, SDataSeries(parent(stages.Φp[1]) .- parent(stages.Φp[2])))
    # end

    try
        plot(sol, stages, idae, dir, file, fig_suff)
    catch ex
        if isa(ex, DomainError)
            show(stdout, "text/markdown", Markdown.parse("DOMAIN ERROR: Diagnostics crashed."))
            @warn("DOMAIN ERROR: Diagnostics crashed.")
        else
            show(stdout, "text/markdown", Markdown.parse("ERROR: Plot crashed with $(typeof(ex))"))
            show(stdout, "text/markdown", Markdown.parse(ex.msg))
            @warn("ERROR: Plot crashed with $(typeof(ex))")
            @warn(ex.msg)
            showerror(stdout, ex, catch_backtrace())
        end
    end
end

function run_list(idae, name, list, plot_dir = PLOT_DIR, symp_dir = SYMP_DIR;
                    fig_suff = ".png", phi_average = nothing)

    for run in list
        tab, file = run

        write_plots(tab, plot_dir, file, name, fig_suff)
        write_symplecticity(tab, symp_dir, file, name)

        show(stdout, "text/markdown", Markdown.parse("### $(tab.name)"))
        _linebreak(stdout)

        show(stdout, "text/markdown", Markdown.parse("[Plots]($file.md)"))
        show(stdout, "text/markdown", Markdown.parse(" • "))
        show(stdout, "text/markdown", Markdown.parse("[Symplecticity]($symp_dir/$file.md)"))
        show(stdout, "text/markdown", Markdown.parse(" • Tableau: "))
        show(stdout, "text/markdown", Markdown.parse("[`$name`](@ref)"))

        _linebreak(stdout)

        run_spark(idae, tab, length(run) ≥ 3 ? run[3] : nt, plot_dir, file, fig_suff, phi_average)
        show(stdout, "text/markdown", Markdown.parse("![$name]($plot_dir/$file$fig_suff)"))
    end

    nothing
end


text:
# High level model construction

function xpress_model(;    # solver environment
    name::String="",               # model name
    sense::Symbol=:minimize,       # :minimize or :maximize
    H::CoeffMat=emptyfmat,         # quadratic coefficient matrix
    f::FVec=emptyfvec,             # linear coefficient vector
    A::CoeffMat=emptyfmat,         # LHS of inequality constraints
    b::FVec=emptyfvec,             # RHS of inequality constraints
    Aeq::CoeffMat=emptyfmat,       # LHS of equality constraints
    beq::FVec=emptyfvec,           # RHS of equality constraints
    lb::Union{Real,Vector}=-Inf,               # upper bounds
    ub::Union{Real,Vector}=Inf)                # lower bounds

    # check f
    if isempty(f)
        error("f must be specified.")
    end

    # create model
    model = Model(name)

    # set sense
    if sense != :minimize
        set_sense!(model, sense)
    end

    # add variables
    add_cvars!(model, f, lb, ub)
    #update_model!(model)

    # add qpterms
    if !isempty(H)
        add_qpterms!(model, H)
    end

    # add constraints
    if !isempty(A) && !isempty(b)
        add_constrs!(model, A, XPRS_LEQ, b)
    end

    if !isempty(Aeq) && !isempty(beq)
        add_constrs!(model, Aeq, XPRS_EQ, beq)
    end
    #update_model!(model)

    return model
end


text:
using Documenter
using MSM
using DataStructures
using OrderedCollections
using Random
using Distributions
using Statistics
using LinearAlgebra
using Plots
using LaTeXStrings
Random.seed!(1234)  #for replicability reasons

makedocs(
         sitename = "MSM.jl",
         modules  = [MSM],
         authors = "Julien Pascal",
    pages =["Home" => "index.md",
        "Installation" => "installation.md",
        "Getting started" => "gettingstarted.md",
        "Functions and Types" => "functions.md",
        "References" => "references.md",
        "Contributing" => "contributing.md",
        "License" => "license.md"
        ],

)

deploydocs(
    repo = "github.com/JulienPascal/MSM.jl.git",
    devbranch = "main",
)


text:
using Base: @propagate_inbounds

mutable struct OffsetVec{T,VV<:AbstractVector{T}} <: AbstractVector{T}
    parent::VV
    offset::Int
end
offset(V::OffsetVec) = V.offset
Base.parent(V::OffsetVec) = V.parent
function Base.show(io::IO, ov::OffsetVec)
    print(io, "Offset ($(ov.offset)) ")
    Base.show(io, ov.parent)
end
function Base.show(io::IO, mime::MIME"text/plain", ov::OffsetVec)
    print(io, "Offset ($(ov.offset)) ")
    Base.show(io, mime, ov.parent)
end
@propagate_inbounds Base.getindex(V::OffsetVec, i::Int) = parent(V)[i-offset(V)]
@propagate_inbounds function Base.setindex!(V::OffsetVec, val, i::Int)
    parent(V)[i-offset(V)] = val
    V
end

mutable struct OffsetMat{T,MM<:AbstractMatrix{T}} <: AbstractMatrix{T}
    parent::MM
    offset1::Int
    offset2::Int
end
offset1(M::OffsetMat) = M.offset1
offset2(M::OffsetMat) = M.offset2
Base.parent(M::OffsetMat) = M.parent
function Base.show(io::IO, om::OffsetMat)
    print(io, "Offset ($(om.offset1),$(om.offset2)) ")
    Base.show(io, om.parent)
end
function Base.show(io::IO, mime::MIME"text/plain", om::OffsetMat)
    print(io, "Offset ($(om.offset1),$(om.offset2)) ")
    Base.show(io, mime, om.parent)
end
@propagate_inbounds Base.getindex(M::OffsetMat, i::Int, j::Int) = parent(M)[i-offset1(M), j-offset2(M)]
@propagate_inbounds Base.getindex(M::OffsetMat, ::Colon, j::Int) = parent(M)[:, j-offset2(M)]
@propagate_inbounds function Base.setindex!(M::OffsetMat, val, i::Int, j::Int)
    parent(M)[i-offset1(M), j-offset2(M)] = val
    M
end


text:
module TheFix
    using REPL

    function cleanse(expr; safe::Bool=true)
        return expr
    end

    function cleanse(expr::Symbol; safe::Bool = true)
        try
            Main.eval(expr)
        catch ex
            expr = replace(expr, ex; safe = safe)
        end
        return expr
    end

    function cleanse(expr::Expr; safe::Bool=true)
        start = expr.head == :(=) && expr.args[1] isa Symbol ? 2 : 1
        for i in start:length(expr.args)
            expr.args[i] = cleanse(expr.args[i]; safe = safe)
        end
        try
            Main.eval(expr)
        catch ex
            expr = replace(expr, ex; safe = safe)
        end
        return expr
    end

    function replace(expr::Symbol, ex::UndefVarError; safe::Bool=true)
        if !(first(string(expr)) == '.' && isdefined(Main, Symbol(string(expr)[2:end])))
            possibilities = REPL.levsort(String(ex.var), REPL.accessible(Main))
            if safe
                i = 1
                while i ≤ min(length(possibilities), 3)
                    if i == 1
                        @info "Couldn't find $(ex.var). Did you mean $(possibilities[i])?"
                    else
                        @info "Did you mean $(possibilities[i])?"
                    end
                    answer = lowercase(strip(readline(stdin)))
                    if isempty(answer) || answer == "y" || answer == "yes"
                        expr = Meta.parse(possibilities[i])
                        break
                    elseif answer == "n" || answer == "no"
                        i += 1
                    else
                        println(stdout, "Unrecognized answer. Answer `y` or `n`.")
                    end
                end
                if i > min(length(possibilities), 3)
                    @info "Couldn't find a fix. What did you mean?"
                    expr = Meta.parse(strip(readline(stdin)))
                end
            else
                expr = Meta.parse(first(possibilities))
            end
            @info "Fixing $ex with $expr."
        end
        return expr
    end

    function replace(expr::Expr, ex::DomainError; safe::Bool=true)
        if occursin("complex", ex.msg) || occursin("NaN result for non-NaN input.", ex.msg)
            for (i, arg) in enumerate(expr.args)
                if Main.eval(arg) isa Number && Main.eval(arg) == ex.val
                    expr.args[i] = :(Complex($arg))
                elseif arg isa Expr
                    if arg.head == :tuple
                        for (j, argt) in enumerate(arg.args)
                            if any(x -> x == ex.val, Main.eval(argt))
                                arg.args[j] = :(Complex.($argt))
                            end
                        end
                    end
                end
            end
            @info "Fixing $ex with $expr."
        elseif occursin("negative power", ex.msg)
            for (i, arg) in enumerate(expr.args)
                if arg == :(.^) || i == 1
                    continue
                elseif Main.eval(arg) isa Number && Main.eval(arg) == ex.val
                    expr.args[i] = :(float($arg))
                elseif Main.eval(arg) isa AbstractArray && any(x -> x == ex.val, Main.eval(arg))
                    expr.args[i] = :(float.($arg))
                elseif Main.eval(arg) isa Tuple && any(x -> x == ex.val, Main.eval(arg))
                    expr.args[i] = :(float.($arg))
                end
            end
            @info "Fixing $ex with $expr."
        else
            @warn "Exception $ex not implemented."
        end
        return expr
    end

    function replace(expr::Expr, ex::OverflowError; safe::Bool=true)
        if occursin("factorial", ex.msg)
            for (i, arg) in enumerate(expr.args)
                if Main.eval(arg) isa Integer
                    expr.args[i] = :(big($arg))
                elseif arg isa Expr
                    if arg.head == :tuple
                        for (j, argt) in enumerate(arg.args)
                            arg.args[j] = :(big.($argt))
                        end
                    end
                end
            end
            @info "Fixing $ex with $expr."
        elseif occursin("gcd", ex.msg)
            for (i, arg) in enumerate(expr.args)
                if Main.eval(arg) isa Integer
                    expr.args[i] = :(widen($arg))
                elseif arg isa Expr
                    if arg.head == :tuple
                        for (j, argt) in enumerate(arg.args)
                            arg.args[j] = :(widen.($argt))
                        end
                    end
                end
            end
            @info "Fixing $ex with $expr."
        elseif occursin("checked", ex.msg) || occursin("overflow", ex.msg)
            for (i, arg) in enumerate(expr.args)
                if Main.eval(arg) isa Integer || Main.eval(arg) isa Rational
                    expr.args[i] = :(widen($arg))
                elseif Main.eval(arg) isa AbstractArray
                    expr.args[i] = :(widen.($arg))
                elseif arg isa Expr
                    if arg.head == :tuple
                        for (j, argt) in enumerate(arg.args)
                            arg.args[j] = :(widen.($argt))
                        end
                    end
                end
            end
            @info "Fixing $ex with $expr."
        else
            @warn "Exception $ex not implemented."
        end
        return expr
    end

    function replace(expr::Expr, ex::DivideError; safe::Bool=true)
        for (i, arg) in enumerate(expr.args)
            if arg == :div
                while length(expr.args) > 3
                    pop!(expr.args)
                end
                expr.args[i] = :(/)
            elseif arg == :(÷) || arg == :cld || arg == :fld
                expr.args[i] = :(/)
            elseif arg == :rem || arg == :mod
                if length(expr.args) > 3
                    x = expr.args[2]
                    y = expr.args[3]
                    r = expr.args[4]
                    expr = :($x - $y*round($x/$y, $r))
                    break
                else
                    x = expr.args[2]
                    y = expr.args[3]
                    expr = :($x - $y*($x/$y))
                    break
                end
            elseif arg == :(%)
                x = expr.args[2]
                y = expr.args[3]
                expr = :($x - $y*($x/$y))
                break
            elseif arg == :divrem
                if length(expr.args) > 3
                    x = expr.args[2]
                    y = expr.args[3]
                    r = expr.args[4]
                    expr = :(($x/$y, $x - $y*round($x/$y, $r)))
                    break
                else
                    x = expr.args[2]
                    y = expr.args[3]
                    expr = :(($x/$y, $x - $y*($x/$y)))
                    break
                end
            end
        end
        @info "Fixing $ex with $expr."
        return expr
    end

    function replace(expr::Expr, ex::BoundsError; safe::Bool=true)
        if length(expr.args) == 2
            arg1 = expr.args[1]
            arg2 = expr.args[2]
            if Main.eval(arg2) isa Integer
                expr.args[2] = :(clamp($arg2, extrema(eachindex($arg1))...))
            elseif Main.eval(arg2) isa AbstractUnitRange
                expr.args[2] = :(intersect($arg2, eachindex($arg1)))
            end
        else
            arg1 = expr.args[1]
            for i in 2:length(expr.args)
                argi = expr.args[i]
                l, r = Main.eval(:(extrema(axes($arg1, $(i-1)))))
                idx = Main.eval(argi)
                if idx isa Integer
                    if first(idx) < l || last(idx) > r
                        expr.args[i] = :(clamp($argi, extrema(axes($arg1, $(i-1)))...))
                    end
                elseif idx isa AbstractUnitRange
                    if first(idx) < l || last(idx) > r
                        expr.args[i] = :(intersect($argi, axes($arg1, $(i-1))))
                    end
                end
            end
        end
        @info "Fixing $ex with $expr."
        return expr
    end

    function replace(expr, ex::Exception; safe::Bool=true)
        @warn "Exception $ex not implemented."
        return expr
    end

    function display_error(io::IO, stack::Vector)
        printstyled(io, "ERROR: "; bold=true, color=Base.error_color())
        bt = Any[ (x[1], Base.scrub_repl_backtrace(x[2])) for x in stack ]
        Base.show_exception_stack(IOContext(io, :limit => true), bt)
    end

    """
        TheFix.@safeword(fix, safe)

    Create a safe word `fix` that may be used to correct common errors in the REPL.
    This creates a struct `fix` with no fields and the fix occurs by showing it.
    The Boolean variable `safe` indicates whether or not the user requests confirmation.

    # Examples
    ```jldoctest
    julia> TheFix.@safeword fix true

    julia> sine(4)
    ERROR: UndefVarError: sine not defined
    Stacktrace:
     [1] top-level scope at REPL[3]:1

    julia> fix
    [ Info: Couldn't find sine. Did you mean sin?
    y
    [ Info: Fixing UndefVarError(:sine) with sin.

    julia> sin(4)
    -0.7568024953079282

    julia> TheFix.@safeword FIX false

    julia> logarithm(3)
    ERROR: UndefVarError: logarithm not defined
    Stacktrace:
     [1] top-level scope at REPL[7]:1

    julia> FIX
    [ Info: Fixing UndefVarError(:logarithm) with log.

    julia> log(3)
    1.0986122886681098
    ```
    """
    macro safeword(fix, safe)
        return esc(quote
            export $fix

            struct $fix end

            function Base.show(io::IO, f::Type{$fix})
                file = open(TheFix.REPL.find_hist_file(), read = true, write = true)
                str = read(file, String)
                itr = findlast("# mode: julia", str)
                itr = findprev("# mode: julia", str, first(itr))
                nxt = itr[end] + 3
                lst = first(findnext("# time:", str, nxt)) - 2
                expr = TheFix.cleanse(Meta.parse(str[nxt:lst]); safe = $safe)
                try
                    printstyled(io, "\n"*TheFix.REPL.JULIA_PROMPT, bold=true; color=202)
                    println(io, string(expr))
                    show(io, "text/plain", Main.eval(expr))
                    hist = Base.active_repl.mistate.current_mode.hist
                    !isempty(hist.history) && isequal(:julia, last(hist.modes)) && string(expr) == last(hist.history) && return
                    push!(hist.modes, :julia)
                    push!(hist.history, string(expr))
                    TheFix.REPL.history_reset_state(hist)
                    entry = """
                    # time: $(Libc.strftime("%Y-%m-%d %H:%M:%S %Z", time()))
                    # mode: julia
                    $(replace(string(expr), r"^"ms => "\t"))
                    """
                    seekend(file)
                    print(file, entry)
                    flush(file)
                catch ex
                    Base.invokelatest(TheFix.display_error, io, Base.catch_stack())
                end
            end
        end)
    end
end # module


text:
#
# Basic two-body orbit simularion written in julia
#
using PyPlot

# Universal gravitation constant
const G = 6.67e-11

infile = open("earth.dat")

#
# Read speed, position and mass of first body
#
linia = readline(infile)
pos = split(strip(linia), " ")
linia = readline(infile)
spd = split(strip(linia), " ")
mass = strip(readline(infile))

# Initialize and load position and velocity vectors
r1 = zeros(3,1)
v1 = zeros(3,1)
for i in 1:3
    r1[i] = parse(Float64,pos[i])
    v1[i] = parse(Float64,spd[i])
end
# Initialize mass
m1 = parse(Float64,mass)

#
# Read speed, position and mass of second body
#
linia = readline(infile)
pos = split(strip(linia), " ")
linia = readline(infile)
spd = split(strip(linia), " ")
mass = strip(readline(infile))

# Initialize and load position and velocity vectors
r2 = zeros(3,1)
v2 = zeros(3,1)
for i in 1:3
    r2[i] = parse(Float64,pos[i])
    v2[i] = parse(Float64,spd[i])
end
# Initialize mass
m2 = parse(Float64,mass)

# Time increment (in seconds) and number of steps
# will be a day.
linia = readline(infile)
times = split(strip(linia), " ")
nintv = parse(Int32,times[2])
dt    = parse(Float64,times[1])

close(infile)

# Newton's formula for gravitational force
function grav(r1,r2,m1,m2)
  r = norm(r1-r2)
  f1 = G * m1 * m2 * (r2 - r1) / r^3
  f2 = G * m1 * m2 * (r1 - r2) / r^3
  return f1,f2
end

# Newton's formula for acceleration
function accel(f,m)
  return f/m
end

# Matrices to contain the (x,y) coordinates computed for each iteration
# for bodies 1 and 2
C1 = zeros(nintv ,2)
C2 = zeros(nintv ,2)

# Compute the orbit position for each step
@sync for i in 1:nintv
    C1[i,1] = r1[1]
    C1[i,2] = r1[2]
    C2[i,1] = r2[1]
    C2[i,2] = r2[2]

    f1,f2 = grav(r1,r2,m1,m2)
    a1 = accel(f1,m1)
    a2 = accel(f2,m2)
    v1 = v1 + a1 * dt
    v2 = v2 + a2 * dt
    r1 = r1 + v1 * dt
    r2 = r2 + v2 * dt
    @printf("x=%+10.5e, y=%+10.5e, z=%+10.5e\n", r1[1],r1[2],r1[3])
    flush(STDOUT)
end

# Plot both orbits
scatter(C1[1:nintv,1],C1[1:nintv,2],marker=".")
scatter(C2[1:nintv,1],C2[1:nintv,2],marker=",")

savefig("orbit.png")

max_stars_repo_name: jguillaumes/orbit-examples

max_issues_repo_name: jguillaumes/orbit-examples

max_forks_repo_name: jguillaumes/orbit-examples


text: include("ModelParametersTest.jl")


text:
import GEMPIC: set_x!, set_v!, set_weights!
import GEMPIC: get_charge, add_charge!
import GEMPIC: push_x_accumulate_j!

@testset "Hamiltonian splitting Boris" begin
    n_particles = 2
    eta_min = 0.0
    eta_max = 4.0π
    num_cells = 16
    delta_t = 0.1
    degree_smoother = 3

    mesh = OneDGrid(eta_min, eta_max, num_cells)

    pg = ParticleGroup{1,2}(n_particles; common_weight=1.0)

    particle_info_ref = reshape(
        [
            11.780972450961723,
            -1.5341205443525459,
            0.15731068461017067,
            6.8636759376074723,
            5.4977871437821380,
            0.15731068461017067,
            -1.5341205443525459,
            5.7026946767517002,
        ],
        4,
        n_particles,
    )

    # Initialize particles from particle_info_ref

    for i_part in 1:n_particles
        xi = particle_info_ref[1, i_part]
        set_x!(pg, i_part, xi)
        xi = particle_info_ref[2:3, i_part]
        set_v!(pg, i_part, xi)
        xi = particle_info_ref[4, i_part]
        set_weights!(pg, i_part, xi[1])
    end

    # Initialize kernel smoother

    kernel_smoother_1 = ParticleMeshCoupling1D(
        mesh, n_particles, degree_smoother - 1, :galerkin
    )
    kernel_smoother_0 = ParticleMeshCoupling1D(
        mesh, n_particles, degree_smoother, :galerkin
    )

    # Initialize Maxwell solver
    maxwell_solver = Maxwell1DFEM(mesh, degree_smoother)

    efield = [
        ones(Float64, (kernel_smoother_0.n_dofs)), ones(Float64, (kernel_smoother_0.n_dofs))
    ]

    bfield = ones(Float64, (kernel_smoother_0.n_dofs))

    rho = zeros(Float64, (kernel_smoother_0.n_dofs))

    for i_part in 1:n_particles
        xi = get_x(pg, i_part)[1]

        wi = get_charge(pg, i_part)

        add_charge!(rho, kernel_smoother_0, xi, wi[1])
    end

    compute_e_from_rho!(efield[1], maxwell_solver, rho)

    propagator = HamiltonianSplittingBoris(
        maxwell_solver, kernel_smoother_0, kernel_smoother_1, pg, efield, bfield
    )

    staggering!(propagator, 0.5 * delta_t)

    strang_splitting!(propagator, delta_t, 1)

    # Compare to reference
    # Particle information after operatorV application
    particle_info_ref = reshape(
        [
            11.590250917552364,
            -1.5236851980054489,
            0.40632305631363136,
            6.8636759376074723,
            5.5028810915315152,
            1.1611806341228641e-002,
            -1.4090781780562809,
            5.7026946767517002,
        ],
        4,
        n_particles,
    )

    for i_part in 1:n_particles
        xi = get_x(pg, i_part)

        @test xi[1] ≈ particle_info_ref[1, i_part]

        xi = get_v(pg, i_part)

        @test xi[1] ≈ particle_info_ref[2, i_part]
        @test xi[2] ≈ particle_info_ref[3, i_part]

        xi = get_charge(pg, i_part)

        @test xi[1] ≈ particle_info_ref[4, i_part]
    end

    bfield_ref = [
        1.0007152502482244,
        0.99445857477120148,
        1.0126078645236389,
        0.97492340832308344,
        1.0502493982082513,
        0.88560153553873056,
        1.1141170903091009,
        0.95020604146500376,
        1.0248316646998006,
        0.98745124326634226,
        1.0056520876592949,
        0.99896954592982978,
        0.99645963114519265,
        1.0122347396384064,
        0.98742382899460557,
        1.0040980952792931,
    ]

    @test maximum(abs.(bfield .- bfield_ref)) ≈ 0.0

    efield_ref = reshape(
        [
            1.7801632633330595,
            0.51717913695931361,
            4.2412959490364575e-002,
            -1.1212332464258261,
            -1.1684637217416221,
            -3.7806966401897091,
            3.7363585338991756,
            1.1932927327940017,
            1.1033167281181586,
            -0.018180082909540971,
            -0.56595990082616909,
            -1.6702327921823792,
            -1.8045561507212267,
            -4.3141455300350842,
            4.8236042130256775,
            1.5537560432220632,
            0.98879315272141322,
            1.0323154046944396,
            0.93329346828214133,
            1.1302445587553216,
            0.73558670810944515,
            1.6340701469431378,
            0.73779661553285170,
            1.1288774513495987,
            0.93385001285635261,
            1.0324077177717137,
            0.98801632510199811,
            0.99610949244376890,
            1.0239154844069209,
            0.92782406399041817,
            1.0265970800925086,
            0.99441071503466349,
        ],
        num_cells,
        2,
    )

    @test maximum(abs.(efield[1] .- efield_ref[:, 1])) ≈ 0.0 atol = 1e-15
    @test maximum(abs.(efield[2] .- efield_ref[:, 2])) ≈ 0.0 atol = 1e-15
end


text:
#https://github.com/travigd/ComparativeAutograder.jl/blob/f9e99eaec800ecb8ca5f919e820b967dd27fe063/src/stdio.jl#L3
#Travis DePrato
"""
    @capture_stdstreams(expr)
Evaluate the specified expression and return a `(stdout, stderr)` tuple.
This macro will attempt to capture every form of output, including via the
`Logging` module (_e.g._ the `@logmsg` and friend macros), but it's unable to
do this with perfect accuracy. Any code that stores a reference to `stdout` or
`stderr` before this is executed will retain those references to the original
streams.
For example, the original log instance may be obtained using
`Logging.global_logger()` because no attempt is made to set the global logger
instance.
# Examples
```julia-repl
julia> out, err = ComparativeAutograder.@capture_stdstreams begin
           println("I'm going to out!")
           println(stderr, "I'm going to err!")
       end
("I'm going to out!\\n", "I'm going to err!\\n")
```
"""
macro capture_stdstreams(expr)
    # We use a whole bunch of `gensym`'s to avoid polluting the namespace
    out, outread, outwrite, outtask = gensym.([:stdout, :outread, :outwrite, :outtask])
    err, errread, errwrite, errtask = gensym.([:stderr, :errread, :errwrite, :errtask])
    return quote
        $out, $err = Base.stdout, Base.stderr
        $outread, $outwrite = redirect_stdout()
        $errread, $errwrite = redirect_stderr()
        $outtask = @async read($outread, String)
        $errtask = @async read($errread, String)
        try
            with_logger(ConsoleLogger(Base.stderr, Logging.Debug)) do
                $(esc(expr))
            end
        catch (e)
            rethrow(e)
        finally
            redirect_stdout($out)
            redirect_stderr($err)
            close.([$outwrite, $errwrite])
        end
        fetch($outtask), fetch($errtask)
    end
end

"""
    evaluate(ui, truth)::Tuple{ Bool, Bool }

Evaluates the output of a user's input against ground truth, and the variable declaration.
"""
function evaluate(ui, var, ans)::Tuple{ Bool, Bool }
    #try to run the users code
    syntax_error, incorrect_answer = false, false;
    ui_ex = Meta.parse( ui )
    try; Meta.eval( ui_ex ); catch; syntax_error = true; end
    # check whether the variable is defined
    try; Meta.eval( var ); catch; incorrect_answer = true; end
    # if defined check whether it is set to the correct answer
    !incorrect_answer && (incorrect_answer = !(Meta.eval(var) === ans))
    return (syntax_error, incorrect_answer)
end

"""
    evaluate(ui, truth)::Tuple{ Bool, Bool }

Evaluates the output of a user's input against ground truth.
"""
function evaluate(ui, truth)::Tuple{ Bool, Bool }
    syntax_error, incorrect_answer = false, false;
    ui_ex = Meta.parse( ui )
    result = nothing
    #maybe use Main.eval?
    try; result = @capture_stdstreams( Meta.eval( ui_ex ) ); catch; syntax_error = true; end
    # Check whether it is set to the correct answer
    if !isnothing( result)
        incorrect_answer = !( first( result ) == truth )
    else
        incorrect_answer = false
    end
    return ( syntax_error, incorrect_answer )
end


text:
module TestExamplesTransducers
include("preamble.jl")
include("../examples/transducers.jl")

@test eltype(addone_out1) === Int
@test length(recall_out1) == 5
@test length(recall_out2) == 5 + RandomRecall().history

end  # module


text:
include("Agent.jl")
include("linspecer.jl")
using Distributed

dt = .1;
A1 = [1. 0 0 0 0 0 dt 0;
      0 1 0 0 0 0 0 dt;
      0 0 1 0 dt 0 -dt 0;
      0 0 0 1 0 dt 0 -dt;
      0 0 0 0 1 0 0 0;
      0 0 0 0 0 1 0 0;
      0 0 0 0 0 0 1 0;
      0 0 0 0 0 0 0 1];

 A2 = [1. 0 0 0 0 0 dt 0;
       0 1 0 0 0 0 0 dt;
       0 0 1 0 0 0 -dt 0;
       0 0 0 1 0 0 0 -dt;
       0 0 0 0 1 0 0 0;
       0 0 0 0 0 1 0 0;
       0 0 0 0 0 0 1 0;
       0 0 0 0 0 0 0 1];

As = Array{Array{Float64,2},1}(undef,2);
As[1] = A1;
As[2] = A2;

C = [0. 0 1 0 0 0 0 0;
     0. 0 0 1 0 0 0 0;
     0. 0 0 0 0 0 1 0;
     0. 0 0 0 0 0 0 1];
muPrior = [0., 0, 10., 0, 0., 1, 1., 0];
endT = 100.;
initVar = [1., 1, 100., 100, 1., 1, 1., 1];
emitVar = [1., 1, 100., 100, 1., 1, 1., 1];
a = 1.;

SWs = empty!(Array{SimWorld,1}(undef,1));
for A in As
      push!(SWs, SimWorld(A, C, muPrior, emitVar, endT, dt))
end

KMs = (Array{KalmanModel,1}(undef,length(As)));
for i in 1:length(As)
      KMs[i] = KalmanModel(As[i], C, muPrior, initVar, a, Integer(ceil(endT/dt))+1);
end

N=100
MusLL = Array{Array{Float64,2},3}(undef, length(SWs), length(KMs), N)
for s in 1:length(SWs), i in 1:N
      Zs, Ys = getStates(SWs[s])
      for k in 1:length(KMs)
            # MusLL[s,k,i] = KMs[k].getMusLL(Zs, Ys)
            MusLL[s,k,i] = runSim(Zs, Ys, KMs[k])
            if i==N
                  println(s,":",k,":",i,"->",MusLL[s,k,i][end-1])
            end
      end
end

findmax!(maxval, ind, A)

___________________________
static = StaticWorld(A, C, muPrior, endT, dt);

flexSigValues = [10, 1, .1]
flexSigInd = 4
flexVarName = "sigmaVest"

v=1

simNames = Array{String, 1}(undef,2)
simNames[1] = "Moving"
simNames[2] = "Still"

set_zero_subnormals(true)
variationNum = length(simNames)
colors = linspecer(variationNum)

simOpts = SimOpts(sigmas, 1000)
plotOpts = PlotOpts(simNames[v], colors[v,:])
kworld = FullWorld(static, simOpts);

max_stars_repo_name: rhoadesScholar/Drugowitsch-Rotation

max_issues_repo_name: rhoadesScholar/Drugowitsch-Rotation

max_forks_repo_name: rhoadesScholar/Drugowitsch-Rotation


text:
using Gadfly, DataFrames

set_default_plot_size(6inch, 3inch)

xs = 0:0.1:20

df_cos = DataFrame(
    x=xs,
    y=cos.(xs),
    ymin=cos.(xs) .- 0.5,
    ymax=cos.(xs) .+ 0.5,
    f="cos"
)

df_sin = DataFrame(
    x=xs,
    y=sin.(xs),
    ymin=sin.(xs) .- 0.5,
    ymax=sin.(xs) .+ 0.5,
    f="sin"
)

df = vcat(df_cos, df_sin)

plot(df, x=:x, y=:y, ymin=:ymin, ymax=:ymax, color=:f, Geom.line, Geom.ribbon)