text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""A StatefulSetSpec is the specification of a StatefulSet.

    IoK8sApiAppsV1beta1StatefulSetSpec(;
        podManagementPolicy=nothing,
        replicas=nothing,
        revisionHistoryLimit=nothing,
        selector=nothing,
        serviceName=nothing,
        template=nothing,
        updateStrategy=nothing,
        volumeClaimTemplates=nothing,
    )

    - podManagementPolicy::String : podManagementPolicy controls how pods are created during initial scale up, when replacing pods on nodes, or when scaling down. The default policy is `OrderedReady`, where pods are created in increasing order (pod-0, then pod-1, etc) and the controller will wait until each pod is ready before continuing. When scaling down, the pods are removed in the opposite order. The alternative policy is `Parallel` which will create pods in parallel to match the desired scale without waiting, and on scale down will delete all pods at once.
    - replicas::Int32 : replicas is the desired number of replicas of the given Template. These are replicas in the sense that they are instantiations of the same Template, but individual replicas also have a consistent identity. If unspecified, defaults to 1.
    - revisionHistoryLimit::Int32 : revisionHistoryLimit is the maximum number of revisions that will be maintained in the StatefulSet's revision history. The revision history consists of all revisions not represented by a currently applied StatefulSetSpec version. The default value is 10.
    - selector::IoK8sApimachineryPkgApisMetaV1LabelSelector : selector is a label query over pods that should match the replica count. If empty, defaulted to labels on the pod template. More info: https://kubernetes.io/docs/concepts/overview/working-with-objects/labels/#label-selectors
    - serviceName::String : serviceName is the name of the service that governs this StatefulSet. This service must exist before the StatefulSet, and is responsible for the network identity of the set. Pods get DNS/hostnames that follow the pattern: pod-specific-string.serviceName.default.svc.cluster.local where \"pod-specific-string\" is managed by the StatefulSet controller.
    - template::IoK8sApiCoreV1PodTemplateSpec : template is the object that describes the pod that will be created if insufficient replicas are detected. Each pod stamped out by the StatefulSet will fulfill this Template, but have a unique identity from the rest of the StatefulSet.
    - updateStrategy::IoK8sApiAppsV1beta1StatefulSetUpdateStrategy : updateStrategy indicates the StatefulSetUpdateStrategy that will be employed to update Pods in the StatefulSet when a revision is made to Template.
    - volumeClaimTemplates::Vector{IoK8sApiCoreV1PersistentVolumeClaim} : volumeClaimTemplates is a list of claims that pods are allowed to reference. The StatefulSet controller is responsible for mapping network identities to claims in a way that maintains the identity of a pod. Every claim in this list must have at least one matching (by name) volumeMount in one container in the template. A claim in this list takes precedence over any volumes in the template, with the same name.
"""
mutable struct IoK8sApiAppsV1beta1StatefulSetSpec <: SwaggerModel
    podManagementPolicy::Any # spec type: Union{ Nothing, String } # spec name: podManagementPolicy
    replicas::Any # spec type: Union{ Nothing, Int32 } # spec name: replicas
    revisionHistoryLimit::Any # spec type: Union{ Nothing, Int32 } # spec name: revisionHistoryLimit
    selector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: selector
    serviceName::Any # spec type: Union{ Nothing, String } # spec name: serviceName
    template::Any # spec type: Union{ Nothing, IoK8sApiCoreV1PodTemplateSpec } # spec name: template
    updateStrategy::Any # spec type: Union{ Nothing, IoK8sApiAppsV1beta1StatefulSetUpdateStrategy } # spec name: updateStrategy
    volumeClaimTemplates::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1PersistentVolumeClaim} } # spec name: volumeClaimTemplates

    function IoK8sApiAppsV1beta1StatefulSetSpec(;podManagementPolicy=nothing, replicas=nothing, revisionHistoryLimit=nothing, selector=nothing, serviceName=nothing, template=nothing, updateStrategy=nothing, volumeClaimTemplates=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1beta1StatefulSetSpec, Symbol("podManagementPolicy"), podManagementPolicy)
        setfield!(o, Symbol("podManagementPolicy"), podManagementPolicy)
        validate_property(IoK8sApiAppsV1beta1StatefulSetSpec, Symbol("replicas"), replicas)
        setfield!(o, Symbol("replicas"), replicas)
        validate_property(IoK8sApiAppsV1beta1StatefulSetSpec, Symbol("revisionHistoryLimit"), revisionHistoryLimit)
        setfield!(o, Symbol("revisionHistoryLimit"), revisionHistoryLimit)
        validate_property(IoK8sApiAppsV1beta1StatefulSetSpec, Symbol("selector"), selector)
        setfield!(o, Symbol("selector"), selector)
        validate_property(IoK8sApiAppsV1beta1StatefulSetSpec, Symbol("serviceName"), serviceName)
        setfield!(o, Symbol("serviceName"), serviceName)
        validate_property(IoK8sApiAppsV1beta1StatefulSetSpec, Symbol("template"), template)
        setfield!(o, Symbol("template"), template)
        validate_property(IoK8sApiAppsV1beta1StatefulSetSpec, Symbol("updateStrategy"), updateStrategy)
        setfield!(o, Symbol("updateStrategy"), updateStrategy)
        validate_property(IoK8sApiAppsV1beta1StatefulSetSpec, Symbol("volumeClaimTemplates"), volumeClaimTemplates)
        setfield!(o, Symbol("volumeClaimTemplates"), volumeClaimTemplates)
        o
    end
end # type IoK8sApiAppsV1beta1StatefulSetSpec

const _property_map_IoK8sApiAppsV1beta1StatefulSetSpec = Dict{Symbol,Symbol}(Symbol("podManagementPolicy")=>Symbol("podManagementPolicy"), Symbol("replicas")=>Symbol("replicas"), Symbol("revisionHistoryLimit")=>Symbol("revisionHistoryLimit"), Symbol("selector")=>Symbol("selector"), Symbol("serviceName")=>Symbol("serviceName"), Symbol("template")=>Symbol("template"), Symbol("updateStrategy")=>Symbol("updateStrategy"), Symbol("volumeClaimTemplates")=>Symbol("volumeClaimTemplates"))
const _property_types_IoK8sApiAppsV1beta1StatefulSetSpec = Dict{Symbol,String}(Symbol("podManagementPolicy")=>"String", Symbol("replicas")=>"Int32", Symbol("revisionHistoryLimit")=>"Int32", Symbol("selector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("serviceName")=>"String", Symbol("template")=>"IoK8sApiCoreV1PodTemplateSpec", Symbol("updateStrategy")=>"IoK8sApiAppsV1beta1StatefulSetUpdateStrategy", Symbol("volumeClaimTemplates")=>"Vector{IoK8sApiCoreV1PersistentVolumeClaim}")
Base.propertynames(::Type{ IoK8sApiAppsV1beta1StatefulSetSpec }) = collect(keys(_property_map_IoK8sApiAppsV1beta1StatefulSetSpec))
Swagger.property_type(::Type{ IoK8sApiAppsV1beta1StatefulSetSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAppsV1beta1StatefulSetSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1beta1StatefulSetSpec }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1beta1StatefulSetSpec[property_name]

function check_required(o::IoK8sApiAppsV1beta1StatefulSetSpec)
    (getproperty(o, Symbol("serviceName")) === nothing) && (return false)
    (getproperty(o, Symbol("template")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAppsV1beta1StatefulSetSpec }, name::Symbol, val)
end


text:
import MechanicalSketch
import MechanicalSketch: empty_figure, sethue, O, WI, HE, EM, FS, finish
import MechanicalSketch: PALETTE, settext, setline, set_figure_height
import MechanicalSketch: circle, arrow, normalize_datarange, line
import MechanicalSketch: @import_expand, @layer, Point, Table
import MechanicalSketch: ∙, °, trace_rotate_hue, noise_between_wavelengths
import MechanicalSketch: color_from_palette, draw_barplot
import Statistics: cor
let

set_figure_height(3777)
empty_figure(filename = joinpath(@__DIR__, "test_31.png"));

if !@isdefined m²
    @import_expand(m, s)
end

global curpoint = O + (-WI / 2 + EM, -HE / 2 + 4.5EM)
include("test_functions_31.jl")
n_velocity = 4
v_calib = collect(range(0.05m/s, 0.5m/s, length = n_velocity))

# Duration for a streamline
Δt = 2.0s
# Number of sample points per set (1 to 9 from tracking backward, 10 (start) to including 20 forward)
n = 20

# Sampling frequency.
f_s = n / Δt

# How far a particle travels over the duration - both directions
streamlinelengths = v_calib * Δt

wave_per_streamline = 2

# Noise spectrum desired spectrum wavelengths
λ_min = streamlinelengths[1] / wave_per_streamline
λ_max = streamlinelengths[end] / wave_per_streamline

# Frequency of interest
f_0 = wave_per_streamline / Δt

# Number of sample sets per velocity
sets_per_v = 10
settext("
For each velocity
    vcalib = $v_calib, we pick $sets_per_v sets of n = $n streamline noise samples over Δt = $(Δt).
    The noise spectrum is coloured between (λmin, λmax) = ($λ_min , $λ_max) for $wave_per_streamline periods over Δt.
    The last sample sets vary harmonically with time.
", curpoint, markup = true)

# Make a matrix with a streamline sample set per column.
# We add more columns below.
sample_matrix = zeros(n, n_velocity * sets_per_v)
for i in 1:n_velocity
    for j in 1:sets_per_v
        col = (i-1) * sets_per_v + j
        sample_matrix[:, col] = samplepoints(streamlinelengths[i], λ_min , λ_max, n)
    end
end
sample_matrix = normalize_datarange(sample_matrix)

# Add periodic (with respect to time) sample sets (add columns to sample_matrix)
n_columns = (n_velocity + 1 ) * sets_per_v
sample_matrix = hcat(sample_matrix, zeros(n, sets_per_v))
linfractions = range(0.0, step = 1 / 3,  length = sets_per_v)
f_periodic = f_0 * linfractions.^2
for (i, f) in enumerate(f_periodic)
    col = n_velocity * (sets_per_v ) + i
    ω = 2π∙f
    # Time range over the window with n samples
    trng = range(1 / (n - 1), 1, length = n) * Δt
    ϕ = i > 1 ? 2π ∙ rand() : -π / 2
    sample_matrix[:, col] = [0.5 + 0.5sin(ω ∙ tim - ϕ) for tim in trng]
end

# Show the sample sets as individual bar plots
for i in 1:(n_velocity + 1)
    global curpoint += (0.0, 2.5EM)
    # Text for this velocity sample set or harmonic set
    if i <= n_velocity
        settext("$(v_calib[i])", curpoint)
    else
        settext(" 0.5 + sin(2π∙f - ϕ)
         ϕ is random, f = ...",
            curpoint + (-EM, 0), markup = true, valign = "bottom")
    end
    # Plot  this velocity sample set or harmonic set
    for j in 1:sets_per_v
        col = (i - 1) * sets_per_v + j
        origo = curpoint + (8.5EM +EM * 3.25 * (j - 1), 0.0)
        draw_barplot(origo, sample_matrix[:, col] * 2EM, 2.25EM, first_sample_no = -9)
        if i > n_velocity
            str = round(s⁻¹, f_periodic[j], digits = 2)
            settext("$str", origo)
        end
    end
end

curpoint += (0, 9EM)

settext("
We will z-transform each sample set to a complex number: Amplitude and phase for $wave_per_streamline wave per set.
Sampling frequency is
    fs = n  / Δt  =  $f_s
We want to extract a value based on each sample set's component with frequency
    f0 = 2 / Δt = $f_0.
After a Fourier transformation, normalized angular frequencies are:
", curpoint, markup = true)

curpoint += (0, 3EM)
ω_s = 2π ∙ f_s
ω_Nyq = π
ω_0n = 2π ∙ f_0 / f_s
explanationcolumn = ["Sampling frequency",
    "    ...unitless, normalized to radians per sample",
    "Nyquist frequency (highest detectable)",
    "Frequency for which to find amplitude and phase",
    "     ...unitless, normalized",
    "     ...as a fraction of Nyquist"]
mathcolumn = ["ωs = 2πfs = $f_s",
    "ωsn = 2π ",
    "ωNyq = π",
    "ω0 = 2π∙f0 = $(round(s⁻¹, 2π∙f_0, digits = 2))",
    "ω0n = 2π∙f0 / fs = $(round(2π∙f_0 / f_s, digits = 2)) = $(round(2∙f_0 / f_s, digits = 2))π ",
    "ω0n / ωNyq = $(ω_0n / π)"]
ta = Table(6, 2, 20EM, 1.2EM)  # 5 rows, 10 columns, colwidth, colheight
for r in 1:size(ta)[1]
    for c in 1:size(ta)[2]
        settext(c ==1 ? explanationcolumn[r] : mathcolumn[r],
            curpoint + (11EM, 0) + ta[r, c], markup = true)
    end
end

curpoint += (0, 6 * 1.2EM + 1EM)
settext("
We will show the extraction of magnitude and phase of frequency ω0n = $(ω_0n / π)∙π
by use of the z-transform. The average magnitude is noted, and the result is the white point:
", curpoint, markup = true)

vmagn1 = []
for i in 1:(n_velocity + 1)
    # Plot z-transformed sample sets for this row
    global curpoint += (0.0, 2.5EM)
    @layer begin
        magn_avg = 0.0
        for j in 1:sets_per_v
            col = (i - 1) * sets_per_v + j
            origo = curpoint + (8.5EM + EM * 3.25 * (j - 0.5), -0.5 * EM)
            vs =  sample_matrix[:, col]
            vŝ = z_transform_contributions(vs, exp(im * ω_0n))
            # Plot contributions while circling and changing hue from start to finish.
            sethue(color_from_palette("red"))
            setline(1)
            trace_rotate_hue(origo, EM * real.(vŝ), EM * imag.(vŝ))
            # Plot coordinate system
            sethue(PALETTE[3])
            arrow(origo, origo + (0.0,  -EM))
            arrow(origo, origo + (EM, 0.0))
            # Find resultant
            ŝ = sum(vŝ)
            magn_avg += hypot(ŝ)  / sets_per_v
            # Plot transform result
            ŝ_pt = EM * Point(real(ŝ), imag(ŝ))
            setline(2)
            sethue("white")
            line(origo, origo + ŝ_pt, :stroke)
            circle(origo + ŝ_pt, 0.03m, :stroke)
        end
    end
    if i <= n_velocity
        push!(vmagn1, magn_avg)
    end
    # text for this row of transformed sample sets
    str_amp = round(magn_avg, digits = 2)
    if i <= n_velocity
        settext("$(v_calib[i]) magn ω0 = $str_amp",
        curpoint, markup = true)
    else
        settext(" 0.5 + sin(2π∙f - ϕ)
            magnitude ω0 = $str_amp",
               curpoint + (-EM, 0), markup = true, valign = "bottom")
    end
end

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

curpoint += (0, 6 * 1.2EM )
correl1 = round(cor(v_calib * s / m , vmagn1), digits = 3)
settext("
The correlation is $(correl1) between the four magnitudes and velocities.

Let us apply a Hanning window, which improves phase detection, and check the impact on magnitude correlation:
", curpoint, markup = true)
sample_matrix = [sample_matrix[i, j] * (sin(π * i / n)^2) for i = 1:n, j = 1:n_columns]

# Show the sample sets as individual plots
for i in 1:(n_velocity + 1)
    global curpoint += (0.0, 2.5EM)
    # Text of this row
    if i <= n_velocity
        settext("$(v_calib[i])", curpoint)
    else
        settext("0.5 + sin(2π∙f), f = ...",
            curpoint + (-EM, 0), markup = true, valign = "bottom")
    end
    # Plot sample set for this row
    for j in 1:sets_per_v
        col = (i - 1) * sets_per_v + j
        origo = curpoint + (8.5EM +EM * 3.25 * (j - 1), 0.0)
        draw_barplot(origo, sample_matrix[:, col] * 2EM, 2.25EM, first_sample_no = -9)
        if i > n_velocity
            str = round(s⁻¹, f_periodic[j], digits = 2)
            settext("$str", origo)
        end
    end
end

curpoint += (0, 0EM)

vmagn2 = []
for i in 1:(n_velocity + 1)
    # Plot z-transformed sample sets for this row
    global curpoint += (0.0, 2.5EM)
    @layer begin
        magn_avg = 0.0
        for j in 1:sets_per_v
            col = (i - 1) * sets_per_v + j
            origo = curpoint + (8.5EM + EM * 3.25 * (j - 0.5), -0.5 * EM)
            vs =  sample_matrix[:, col]
            vŝ = z_transform_contributions(vs, exp(-im * ω_0n))
            # Plot contributions while circling
            sethue(color_from_palette("red"))
            setline(1)
            trace_rotate_hue(origo, EM * real.(vŝ), EM * imag.(vŝ))
            # Plot coordinate system
            sethue(PALETTE[3])
            arrow(origo, origo + (0.0,  -EM))
            arrow(origo, origo + (EM, 0.0))
            # Plot resultant
            ŝ = sum(vŝ)
            magn_avg += hypot(ŝ)  / sets_per_v
            ŝ_pt = EM * Point(real(ŝ), imag(ŝ))
            setline(2)
            sethue("white")
            line(origo, origo + ŝ_pt, :stroke)
            circle(origo + ŝ_pt, 0.03m, :stroke)
        end
    end
    if i <= n_velocity
        push!(vmagn2, magn_avg)
    end
    # text for this row of transformed sample sets
    str_amp = round(magn_avg, digits = 2)
    if i <= n_velocity
        settext("$(v_calib[i]) magn ω0 = $str_amp",
        curpoint, markup = true)
    else
        settext(" 0.5 + sin(2π∙f - ϕ)
            magnitude ω0 = $str_amp",
               curpoint + (-EM, 0), markup = true, valign = "bottom")
    end
end

curpoint += (0, 4EM)
correl2 = round(cor(v_calib * s / m , vmagn2), digits = 3)
settext("
The correlation is $correl2 between the four magnitudes and velocities, which is $(correl1 > correl2 ? "worse" : "better") than $(correl1).
In most cases, the correlation is worse.
", curpoint, markup = true)

finish()
set_figure_height(round(Int, 3777 / 2))
correl1
end

ext: jl


text:
using TerrainMap
using Base.Test

# write your own tests here
@test 1 == 1

@test typeof(Terrain_Map()) == Terrain_Map
#Test not working@test
#@test typeof(Terrain_Map("TestJulia.map")) == Terrain_Map


text:
macro try_printerror(ex)
    return quote
        try
            $(esc(ex))
        catch e
            if isa(e, InterruptException)
                rethrow()
            end
            @error(e)
            for (exc, bt) in Base.catch_stack()
                showerror(stdout, exc, bt)
                println()
            end
        end
    end
end

"""
    tee_task(io_in, io_outs...)

Creates an asynchronous task that reads in from `io_in` in buffer chunks, writing
available bytes out to all elements of `io_outs` as quickly as possible until `io_in` is
closed.  Closes all elements of `io_outs` once that is finished.  Returns the `Task`.
"""
function tee_task(io_in, io_outs...)
    return @async begin
        @try_printerror begin
            total_size = 0
            while !eof(io_in)
                chunk = readavailable(io_in)
                total_size += length(chunk)
                for io_out in io_outs
                    write(io_out, chunk)
                end
            end
            for io_out in io_outs
                close(io_out)
            end
            return total_size
        end
    end
end

"""
    wait_first(args...)

Return the first waitable (Task, Condition, etc...) that becomes ready.
"""
function wait_first(args...)
    if isempty(args)
        return nothing
    end

    c = Channel()
    for arg in args
        @async begin
            wait(arg)
            put!(c, arg)
        end
    end
    return take!(c)
end

function try_open(args...; kwargs...)
    try
        return open(args...; kwargs...)
    catch e
        if isa(e, InterruptException)
            rethrow(e)
        end
    end
    return nothing
end

max_stars_repo_name: StefanKarpinski/PkgServer

max_issues_repo_name: StefanKarpinski/PkgServer

max_forks_repo_name: StefanKarpinski/PkgServer


text:
@testset "319.bulb-switcher.jl" begin
    @test bulb_switch(3) == 1
    @test bulb_switch(4) == 2
    @test bulb_switch(5) == 2
    @test bulb_switch(31096836) == 5576
end


text:
# We try to test approximation for all function sets, except those that
# are currently known to fail for lack of an implementation.
supports_approximation(s::Dictionary) = true
supports_interpolation(s::Dictionary) = isbasis(s)
# Pick a simple function to approximate
suitable_function(s::Dictionary1d) = x->exp(x/supremum(support(s)))
suitable_function(s::Dictionary) = x -> prod(x)

# Make a tensor product of suitable functions
function suitable_function(dict::TensorProductDict)
    if dimension(dict) == 2
        f1 = suitable_function(component(dict,1))
        f2 = suitable_function(component(dict,2))
        (x,y) -> f1(x)*f2(y)
    elseif dimension(dict) == 3
        f1 = suitable_function(component(dict,1))
        f2 = suitable_function(component(dict,2))
        f3 = suitable_function(component(dict,3))
        (x,y,z) -> f1(x)*f2(y)*f3(z)
    end
    # We should never get here
end
# Make a suitable function by undoing the map
function suitable_function(dict::MappedDict)
    f = suitable_function(superdict(dict))
    m = inverse_map(dict)
    x -> f(m(x))
end
function suitable_function(dict::BasisFunctions.MappedDict2d)
    f = suitable_function(superdict(dict))
    m = inverse_map(dict)
    (x,y) -> f(m(x,y)...)
end
function suitable_function(dict::WeightedDict1d)
    f = suitable_function(superdict(dict))
    g = weightfunction(dict)
    x -> g(x) * f(x)
end
function suitable_function(dict::WeightedDict2d)
    f = suitable_function(superdict(dict))
    g = weightfunction(dict)
    (x,y) -> g(x, y) * f(x, y)
end
suitable_function(dict::OperatedDict) = suitable_function(src(dict))

function suitable_interpolation_grid(basis::Dictionary)
    if BF.hasinterpolationgrid(basis)
        interpolation_grid(basis)
    else
        T = domaintype(basis)
        # A midpoint grid avoids duplication of the endpoints for a periodic basis
        MidpointEquispacedGrid(length(basis), affine_point_in_domain(basis, T(0)), affine_point_in_domain(basis, T(1)))
    end
end

function test_generic_dict_boundschecking(basis)
    ELT = coefficienttype(basis)
    # Bounds checking
    # disable periodic splines for now, since sometimes left(basis,idx) is not
    # in_support currently...
    if (dimension(basis) == 1)
        s = support(basis)
        l = infimum(s)
        r = supremum(s)
        if ~isinf(l)
            @test in_support(basis, l)
            # This assumes test_tolerance(ELT) is equal to tolerance(dict), should probably use the latter here
            @test in_support(basis, l.-1/10*test_tolerance(ELT))
            @test ~in_support(basis, l.-1)
        end
        if ~isinf(r)
            @test in_support(basis, r)
            @test in_support(basis, r.+1/10*test_tolerance(ELT))
            @test ~in_support(basis, r.+1)
        end
        if ~isinf(l) && ~isinf(r) && ~isa(basis,BasisFunctions.PiecewiseDict)
            @test in_support(basis, 1, 1/2*(l + r))
        end
    end
end

function test_generic_dict_native_indexing(basis)
    ## Test iteration over the set
    n = length(basis)
    equality = true
    l = 0
    for i in eachindex(basis)
        l += 1
        # Is conversion from logical to native index and back a bijection?
        idxn = native_index(basis, l)
        idx = linear_index(basis, idxn)
        equality = equality & (idx == l)
    end
    @test l == n
    @test equality
end

function test_generic_dict_indexing(basis)
    ## Does indexing work as intended?
    idx = random_index(basis)
    bf = basis[idx]
    # Test below disabled because the assumption that the set of a basis function
    # is always the set that was indexed is false, e.g., for multidicts.
    #@test dictionary(bf) == basis

    @test linear_index(basis, lastindex(basis)) == length(basis)

    # Is a boundserror thrown when the index is too large?
    @test try
        basis[length(basis)+1]
        false
    catch
        true
    end

    # Error thrown for a range that is out of range?
    @test try
        basis[1:length(basis)+1]
        false
    catch
        true
    end
end

function test_generic_dict_evaluation(basis)
    ELT = coefficienttype(basis)
    idx = random_index(basis)
    bf = basis[idx]
    x = fixed_point_in_domain(basis)
    @test bf(x) ≈ eval_element(basis, idx, x)

    if ! (support(basis) isa DomainSets.FullSpace)
        x_outside = point_outside_domain(basis)
        @test bf(x_outside) == 0
    end

    # Create a random expansion in the basis to test expansion interface
    e = random_expansion(basis)
    coef = coefficients(e)

    ## Does evaluating an expansion equal the sum of coefficients times basis function calls?
    x = fixed_point_in_domain(basis)
    @test e(x) ≈ sum([coef[i] * basis[i](x) for i in eachindex(coef)])

    ## Test evaluation on an array
    x_array = [random_point_in_domain(basis) for i in 1:10]
    z = map(e, x_array)
    @test  z ≈ ELT[ e(x_array[i]) for i in eachindex(x_array) ]

    # Test dictionary evaluation
    x = random_point_in_domain(basis)
    @test norm(basis(x) - [eval_element(basis, i, x) for i in eachindex(basis)]) < test_tolerance(ELT)
end

function test_generic_dict_coefficient_linearization(basis)
    e = random_expansion(basis)
    coef = coefficients(e)
    # Test linearization of coefficients
    linear_coefs = zeros(coefficienttype(basis), length(basis))
    BasisFunctions.linearize_coefficients!(basis, linear_coefs, coef)
    coef2 = BasisFunctions.delinearize_coefficients(basis, linear_coefs)
    @test coef ≈ coef2
end

function test_generic_dict_grid(basis)
    grid1 = interpolation_grid(basis)
    @test length(grid1) == length(basis)
    e = random_expansion(basis)
    z1 = e(grid1)
    z2 = [ e(grid1[i]) for i in eachindex(grid1) ]
    @test z1 ≈ z2
    E = evaluation(basis, GridBasis(basis))
    z3 = E * coefficients(e)
    @test z1 ≈ z3
end

function test_generic_dict_codomaintype(basis)
    T = domaintype(basis)
    FT = prectype(T)
    n = length(basis)
    types_correct = true
    # The comma in the line below is important, otherwise the two static vectors
    # are combined into a statix matrix.
    for x in [ fixed_point_in_domain(basis), rationalize(affine_point_in_domain(basis, FT(0.5))) ]
        if length(basis) > 1
            indices = [1 2 n>>1 n-1 n]
        else
            # For a singleton subset we can only use index 1
            indices = 1
        end
        for idx in indices
            z = eval_element(basis, idx, x)
            types_correct = types_correct & (typeof(z) == codomaintype(basis))
        end
    end
    @test types_correct
end

function test_generic_dict_approximation(basis)
    A = approximation(basis)
    f = suitable_function(basis)
    e = Expansion(basis, A*f)

    # We choose a fairly large error, because the ndof's can be very small.
    # We don't want to test convergence, only that something terrible did
    # not happen, so an error of 2e-3 will do.
    x = random_point_in_domain(basis)
    @test abs(e(x)-f(x...)) < 2e-3
end

truncate(d::Domain) = d
truncate(d::FullSpace{T}) where {T} = -T(10)..T(10)
truncate(d::ProductDomain) = ProductDomain(map(truncate, components(d))...)

function test_gram_projection(basis)
    if hasmeasure(basis)
        G = gram(basis)
        @test src(G) == basis
        @test dest(G) == basis
        n = length(basis)
        @test size(G) == (n,n)
        z = zero(prectype(basis))
        for i in 1:n
            for j in 1:n
                DD = abs(G[i,j] - innerproduct(basis, i, basis, j))
                z = max(z,abs(G[i,j] - innerproduct(basis, i, basis, j)))
            end
        end
        @test z < 1000test_tolerance(prectype(basis))

        # No efficient implementation for BigFloat to construct full gram matrix.
        if !(prectype(basis)==BigFloat)
            f = suitable_function(basis)
            μ = measure(basis)
            # Do we compute the projection integrals accurately?
            Z = integral(t->f(t...)*BasisFunctions.unsafe_eval_element(basis, 1, t)*DomainIntegrals.unsafe_weightfun(μ,t), truncate(support(basis, 1)))
            @test abs(innerproduct(t->f(t...), basis[1]) - Z)/abs(Z) < 1e-1
            @test abs(innerproduct(basis[1], t->f(t...)) - Z)/abs(Z) < 1e-1
            e = approximate(basis, t->f(t...); discrete=false, rtol=1e-6, atol=1e-6)
            x = random_point_in_domain(basis)
            @test abs(e(x)-f(x...)) < 2e-3
        end
    end
end

function test_generic_dict_interpolation(basis)
    ELT = coefficienttype(basis)
    g = suitable_interpolation_grid(basis)
    I = interpolation(basis, g)
    x = rand(GridBasis{coefficienttype(basis)}(g))
    e = Expansion(basis, I*x)
    @test maximum(abs.(e(g)-x)) < 100test_tolerance(ELT)
end

function test_generic_dict_domaintype(basis)
    T = domaintype(basis)
    # Test type promotion
    @test domaintype(convert(Dictionary{T}, basis)) == T
    try
        # We attempt to widen the type. This will throw an exception if
        # the dictionary does not implement `similar` for domain types different
        # from its own. We just ignore that. Otherwise, if no exception is thrown,
        # we test whether the promotion succeeded.
        T2 = widen_type(T)
        basis2 = convert(Dictionary{T2}, basis)
        @test domaintype(basis2) == T2
    catch
    end
end

function test_generic_dict_evaluation_different_grid(basis)
    ELT = coefficienttype(basis)
    T = domaintype(basis)
    n = length(basis)
    e = random_expansion(basis)

    # Test evaluation on a different grid on the support of the basis
    a = infimum(support(basis))
    b = supremum(support(basis))

    if isinf(a)
        a = -T(1)
    end
    if isinf(b)
        b = T(1)
    end

    grid2 = EquispacedGrid(n+3, T(a)+test_tolerance(ELT), T(b)-test_tolerance(ELT))
    z = e(grid2; warnslow = false)
    @test z ≈ ELT[ e(grid2[i]) for i in eachindex(grid2) ]

    # Test evaluation on a grid of a single basis function
    idx = random_index(basis)
    z = eval_element(basis, idx, grid2)
    @test  z ≈ ELT[ basis[idx](grid2[i]) for i in eachindex(grid2) ]
end

function test_generic_dict_transform(basis, grid = interpolation_grid(basis))
    T = domaintype(basis)
    EPS = test_tolerance(T)

    tbasis = GridBasis{coefficienttype(basis)}(grid)

    @test hastransform(basis, grid)
    @test hastransform(basis, tbasis)

    t = transform(tbasis, basis)
    it = transform(basis, tbasis)

    # - transform from grid
    x = random_expansion(tbasis)
    e = t*x
    @test abs(e(grid[1]) - x.coefficients[1]) < EPS
    @test maximum(abs.(e(grid)-x.coefficients)) < EPS
    # - transform to grid
    e = random_expansion(basis)
    x = it*e
    @test abs(e(grid[1]) - x.coefficients[1]) < EPS
    @test abs(e(grid[end]) - x.coefficients[end]) < EPS

    @test maximum(abs.(matrix(t)'-matrix(t'))) < EPS
    @test maximum(abs.(matrix(it)'-matrix(it'))) < EPS
end

function test_generic_dict_evaluation_operator(basis)
    ELT = coefficienttype(basis)
    ## Test evaluation operator
    g = suitable_interpolation_grid(basis)
    E = evaluation(basis, g)
    e = random_expansion(basis)
    y = E*e
    @test maximum([abs.(e(g[i])-y[i]) for i in eachindex(g)]) < test_tolerance(ELT)
end

function test_generic_dict_antiderivative(basis)
    ELT = coefficienttype(basis)
    T = domaintype(basis)
    FT = prectype(T)
    coef = coefficients(random_expansion(basis))

    for dim in 1:dimension(basis)
        if dimension(basis) == 1
            D = antidifferentiation(basis)
        else
            N = dimension(basis)
            order = dimension_tuple(Val{N}(), dim)
            D = antidifferentiation(basis, order)
        end
        @test basis == src(D)
        antidiff_dest = dest(D)

        coef1 = random_expansion(basis)
        coef2 = D*coef
        e1 = Expansion(basis, coef)
        e2 = Expansion(antidiff_dest, coef2)

        x = fixed_point_in_domain(basis)
        delta = test_tolerance(ELT)/10
        N = dimension(basis)
        if N > 1
            unit_vector = zeros(FT, dimension(basis))
            unit_vector[dim] = 1
            x2 = x + SVector{N}(delta*unit_vector)
        else
            x2 = x+delta
        end
        @test abs( (e2(x2)-e2(x))/delta - e1(x) ) / abs(e1(x)) < 2000*test_tolerance(ELT)
    end
end

function test_generic_dict_derivative(basis)
    ELT = coefficienttype(basis)
    T = domaintype(basis)
    FT = prectype(T)
    for dim in 1:dimension(basis)
        # TODO: Sort out problem with dim and multidict
        if dimension(basis) == 1
            D = differentiation(basis)
        else
            N = dimension(basis)
            order = dimension_tuple(Val{N}(), dim)
            D = differentiation(basis, order)
        end
        if dimension(basis)>1
            D = differentiation(basis; dim=dim)
        else
            D = differentiation(basis)
        end
#        @test basis == src(D)
        diff_dest = dest(D)
        coef = coefficients(random_expansion(basis))

        coef1 = random_expansion(basis)
        coef2 = D*coef
        e1 = Expansion(basis, coef)
        e2 = Expansion(diff_dest, coef2)

        x = fixed_point_in_domain(basis)
        delta = test_tolerance(ELT)/10
        N = dimension(basis)
        if N > 1
            unit_vector = zeros(FT, dimension(basis))
            unit_vector[dim] = 1
            x2 = x + SVector{N}(delta*unit_vector)
        else
            x2 = x+delta
        end
        @test abs( (e1(x2)-e1(x))/delta - e2(x) ) / abs(e2(x)) < 2000*test_tolerance(ELT)
    end

    if dimension(basis) == 1
        @test derivative_dict(basis, 0) == basis

        x = fixed_point_in_domain(basis)
        D = differentiation(basis)
        # Verify derivatives in three basis functions: the first, the last,
        # and the middle one
        i1 = 1
        i2 = length(basis)
        i3 = (i1+i2) >> 1

        c1 = Expansion(basis,zeros(basis))
        c1[i1] = 1
        u1 = D*c1
        @test abs(u1(x) - eval_element_derivative(basis, i1, x)) < test_tolerance(ELT)

        c2 = Expansion(basis,zeros(basis))
        c2[i2] = 1
        u2 = D*c2
        @test abs(u2(x) - eval_element_derivative(basis, i2, x)) < test_tolerance(ELT)

        c3 = Expansion(basis,zeros(basis))
        c3[i3] = 1
        u3 = D*c3
        @test abs(u3(x) - eval_element_derivative(basis, i3, x)) < test_tolerance(ELT)
    end

    # TODO: experiment with this test and enable
    # if dimension(basis) == 1 && hasderivative(basis,2)
    #     D2 = differentiation(basis, 2)
    #     e = random_expansion(basis)
    #     f = D2*e
    #     x = fixed_point_in_domain(basis)
    #     delta = sqrt(test_tolerance(ELT))
    #     @test abs(f(x+delta) -2f(x)+f(x-delta))/delta^2 < 100sqrt(test_tolerance(ELT))
    # end
end

function test_generic_dict_interface(basis)
    ELT = coefficienttype(basis)
    T = domaintype(basis)
    FT = prectype(T)
    RT = codomaintype(basis)

    n = length(basis)
    if hasmeasure(basis)
        @test isorthogonal(basis) == isorthogonal(basis, measure(basis))
        @test isorthonormal(basis) == isorthonormal(basis, measure(basis))
        @test isbiorthogonal(basis) == isbiorthogonal(basis, measure(basis))
    end

    if isorthonormal(basis)
        @test isorthogonal(basis)
    end

    if isorthogonal(basis)
        @test isbiorthogonal(basis)
    end

    test_generic_dict_domaintype(basis)

    ## Test dimensions
    s = size(basis)
    for i = 1:length(s)
        @test size(basis, i) == s[i]
    end

    test_generic_dict_boundschecking(basis)

    test_generic_dict_native_indexing(basis)

    test_generic_dict_indexing(basis)

    test_generic_dict_evaluation(basis)

    test_generic_dict_coefficient_linearization(basis)

    ## Verify evaluation on the associated grid
    if BF.hasinterpolationgrid(basis)
        test_generic_dict_grid(basis)
    end

    ## Test output type of calling function
    test_generic_dict_codomaintype(basis)

    if dimension(basis) == 1
        test_generic_dict_evaluation_different_grid(basis)
    end

    ## Test extensions
    if hasextension(basis)
        n2 = extensionsize(basis)
        basis2 = resize(basis, n2)
        E = extension(basis, basis2)
        e1 = random_expansion(basis)
        e2 = E * e1
        x1 = affine_point_in_domain(basis, 1/2)
        @test e1(x1) ≈ e2(x1)
        x2 = affine_point_in_domain(basis, 0.3)
        @test e1(x2) ≈ e2(x2)

        R = restriction(basis2, basis)
        e3 = R * e2
        @test e1(x1) ≈ e3(x1)
        @test e1(x2) ≈ e3(x2)
    end

    # Verify whether evaluation in a larger grid works
    if BF.hasextension(basis) && BF.hasinterpolationgrid(basis)
        basisext = extend(basis)
        grid_ext = interpolation_grid(basisext)
        L = evaluation(basis, grid_ext)
        e = random_expansion(basis)
        z = L*e
        L2 = evaluation(basisext, grid_ext) * extension(basis, basisext)
        z2 = L2*e
        @test maximum(abs.(z-z2)) < 20test_tolerance(ELT)
        # In the future, when we can test for 'fastness' of operators
        # @test isfast(L2) == isfast(L)
    end

    ## Test derivatives
    if BF.hasderivative(basis)
        test_generic_dict_derivative(basis)
    end

    ## Test antiderivatives
    if BF.hasantiderivative(basis)
        test_generic_dict_antiderivative(basis)
    end

    ## Test associated transform
    if BF.hastransform(basis)
        test_generic_dict_transform(basis)
    end

    ## Test interpolation operator on a suitable interpolation grid
    if supports_interpolation(basis)
        test_generic_dict_interpolation(basis)
    end

    test_generic_dict_evaluation_operator(basis)

    ## Test approximation operator
    if supports_approximation(basis)
        test_generic_dict_approximation(basis)
        test_gram_projection(basis)
    end

    if hasmeasure(basis)
        test_orthogonality_orthonormality(basis)
    end
end


text:
"""
    HiddenStateModel{S, T<:TimeType} <: AbstractModel{S}

Abstract type for any model of the hidden state of type `S`.
"""
abstract type HiddenStateModel{S, T} <: AbstractModel{S} end

"""
    state_dim(model::HiddenStateModel)

Returns the dimensionality of the hidden state in `model`.
"""
function state_dim(model::HiddenStateModel) end

"""
    initial_condition(model::HiddenStateModel)

Returns the specification of the initial condition in `model`.
This is either a fixed value or a samplable distribution.
"""
function initial_condition(model::HiddenStateModel) end

"""
    state_type(model::HiddenStateModel)

Returns the data type of the hidden state in `model`.
"""
state_type(model::HiddenStateModel{S,T}) where {S,T} = S

"""
    time_type(model::HiddenStateModel)

Returns the time type of `model`, e.g. `DiscreteTime' or `ContinuousTime'.
"""
time_type(model::HiddenStateModel{S,T}) where {S,T} = T

"""
    initialize(model::HiddenStateModel)

Returns a sample from the initial distribution of `model`.
"""
function initialize(state_model::HiddenStateModel{S, T})::S where {S, T}
    init = initial_condition(state_model)
    if init isa Distributions.Sampleable
        return rand(init)
    else
        return init
    end
end

"""
    propagate!(state[s], model::HiddenStateModel[, dt])

Propagates the state(s) according to the model.
For `ContinuousTime' models, a time step `dt' has to be provided.
Multiple states are given as a matrix with columns corresponding to states, and are processed i.i.d.
"""
function propagate!(state_model::HiddenStateModel{S, T})::S where {S, T} end

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
@testset "plot.jl" begin

model_test(x, a, b) = a .* exp.(-x .* b)

xdata_test = linspace(0, 10, 100)
ydata_test = model_test(xdata_test, 1, 2) + 0.01*randn(length(xdata_test))
eydata_test = 0.01;

fitter_a_test = Fitter(model_test; autoplot=false)
fitter_b_test = Fitter(model_test; autoplot=false)

fig_a_test = set_data!(fitter_a_test, xdata_test, ydata_test, eydata_test) |> plot!
fig_b_test = set_data!(fitter_b_test, xdata_test, ydata_test, eydata_test) |> plot!

@test fig_a_test ≠ fig_b_test

fit!(fitter_a_test)

@test fig_a_test == ignore_outliers!(fitter_a_test, 1) |> fit! |> plot!

@test fig_a_test ≠ plot(fitter_a_test)

end


text:
module Server

# --------------------------------
# Include
# --------------------------------

include("startegy/federated_averaging.jl")

include("client_manager.jl")
include("central_node.jl")
include("service.jl")
include("app.jl")

# --------------------------------
# Export
# --------------------------------

export federated_averaging, ClientManager, length, num_available_clients,
    register, sample_clients, CentralNode, register_client, start, build_router

end # module


text:
using JultraDark
using Test
using MPI

@testset "Actual time step" begin
    @test JultraDark.actual_time_step(0.5, 1, 2)[1] ≈ 0.5
    @test JultraDark.actual_time_step(0.5, 1, 2)[2] == 2

    @test JultraDark.actual_time_step(0.5, 1, 3)[1] ≈ 0.5
    @test JultraDark.actual_time_step(0.5, 1, 3)[2] == 2

    @test JultraDark.actual_time_step(0.01, 1, 10)[1] ≈ 0.01
    @test JultraDark.actual_time_step(0.01, 1, 10)[2] == 10

    @test JultraDark.actual_time_step(0.7, 1, 2)[1] ≈ 0.5
    @test JultraDark.actual_time_step(0.7, 1, 2)[2] == 2

    @test JultraDark.actual_time_step(0.7, 2, 3)[1] ≈ 2/3
    @test JultraDark.actual_time_step(0.7, 2, 3)[2] == 3
end

@testset "Take n steps" begin
    grids = PencilGrids(1.0, 16)
    output_config = OutputConfig(mktempdir(), [])
    @test JultraDark.take_steps!(grids, 0, 1, 10, output_config, t -> 1) == 10.0
end

# Test that evolve to arrives at right time
@testset "Evolve from $t_begin to $t_end" for t_begin in [0, 0.123, 1], t_end in [1.23, 2]
    grids = PencilGrids(1.0, 16)

    a = 1
    grids.ψx .= grids.dist ./ 1e9 # Set ψx to something non-zero
    grids.ψk .= (grids.fft_plan * grids.ψx)
    grids.ρx .= abs2.(grids.ψx)
    grids.Φk .= -4 * π * (grids.rfft_plan * grids.ρx) ./ (a * grids.rk.^2)
    grids.Φk[1, 1, 1] = 0
    grids.Φx .= grids.rfft_plan \ grids.Φk

    output_config = OutputConfig(mktempdir(), [])
    options = JultraDark.Config.SimulationConfig(10, t->1)
    @test JultraDark.evolve_to!(t_begin, t_end, grids, output_config, options) ≈ t_end
end


text:
Base.promote_rule(
    ::Type{Character{T,Cl}},
    ::Type{VirtualCharacter{S,Cl}},
) where {T,S,Cl} = VirtualCharacter{promote_type(T, S),Cl}

for f in (:+, :-)
    @eval begin
        function Base.$f(χ::Character{T,Cl}, ψ::Character{S,Cl}) where {T,S,Cl}
            @assert χ.inv_of == ψ.inv_of
            return VirtualCharacter(
                $f.(values(χ), values(ψ)),
                χ.inv_of,
                conjugacy_classes(χ),
            )
        end
    end
end

for C in (:Character, :VirtualCharacter)
    @eval begin
        Base.:*(χ::$C, c::Number) =
            VirtualCharacter(c .* values(χ), χ.inv_of, conjugacy_classes(χ))

        Base.:*(c::Number, χ::$C) =
            VirtualCharacter(c .* values(χ), χ.inv_of, conjugacy_classes(χ))

        Base.:/(χ::$C, c::Number) =
            VirtualCharacter(values(χ) ./ c, χ.inv_of, conjugacy_classes(χ))
    end
end


text:
export DNN_CTC_With_Softmax
export DNN_Batch_CTC_With_Softmax
export RNN_Batch_CTC_With_Softmax
export CRNN_Batch_CTC_With_Softmax

"""
    DNN_CTC_With_Softmax(x::Variable{Array{T}}, seq; blank=1, weight=1.0)

case batchsize==1 for test case. `x` is the output of a whole complete input sequence

# Inputs
`x`      : 2-D Variable, input sequence\n
`seq`    : 1-D Array, input sequence's label\n
`weight` : weight for CTC loss

# Structure
    ┌───┐
    │ │ │
    │ W ├──►─┐
    │ │ │    │
    └───┘    │
    ┌───┐    │    ┌───┐          ┌───┐
    │ │ │  ┌─┴─┐  │ │ │ softmax  │ │ │   ┌───────┐
    │ Z ├─►│ × ├─►│ X ├─────────►│ P ├──►│CTCLOSS│◄── (seqLabel)
    │ │ │  └───┘  │ │ │          │ │ │   └───┬───┘
    └───┘         └─┬─┘          └─┬─┘       │
                    │              │+        ▼
                  ┌─┴─┐            ▼       ┌─┴─┐
                  │ │ │          ┌─┴─┐ -   │ │ │
                  │ δ │◄─────────┤ - │──◄──┤ r │
                  │ │ │          └───┘     │ │ │
                  └───┘                    └───┘
"""
function DNN_CTC_With_Softmax(x::Variable{Array{T}}, seq; blank=1, weight=1.0) where T
    p = softmax(ᵛ(x); dims=1)
    L = length(seq) * 2 + 1
    r, loglikely = CTC(p, seq, blank=blank)
    if x.backprop
        function DNN_CTC_With_Softmax_Backward()
            if need2computeδ!(x)
                if weight==1.0
                    δ(x) .+=  p - r
                else
                    δ(x) .+= (p - r) .* weight
                end
            end
        end
        push!(graph.backward, DNN_CTC_With_Softmax_Backward)
    end
    return loglikely / L
end

"""
    DNN_Batch_CTC_With_Softmax(x::Variable{Array{T}}, seqlabels::Vector, inputlens; blank=1, weight=1.0) where T

a batch of concatenated input sequence is processed by neural networks into `x`

# Inputs
`x`         : 2-D Variable, inputs of softmax\n
`seqlabels` : a batch of sequential labels, like [[i,j,k],[x,y],...]\n
`inputlens` : records each input sequence's length, like [20,17,...]\n
`weight`    : weight for CTC loss

# Structure
    ┌───┐
    │ │ │
    │ W ├──►─┐
    │ │ │    │
    └───┘    │
    ┌───┐    │    ┌───┐          ┌───┐
    │ │ │  ┌─┴─┐  │ │ │ softmax  │ │ │   ┌───────┐
    │ Z ├─►│ × ├─►│ X ├─────────►│ P ├──►│CTCLOSS│◄── (seqLabel)
    │ │ │  └───┘  │ │ │          │ │ │   └───┬───┘
    └───┘         └─┬─┘          └─┬─┘       │
                    │              │+        ▼
                  ┌─┴─┐            ▼       ┌─┴─┐
                  │ │ │          ┌─┴─┐ -   │ │ │
                  │ δ │◄─────────┤ - │──◄──┤ r │
                  │ │ │          └───┘     │ │ │
                  └───┘                    └───┘
"""
function DNN_Batch_CTC_With_Softmax(x::Variable{Array{T}}, seqlabels::Vector, inputlens; blank=1, weight=1.0) where T
    batchsize = length(inputLengths)
    loglikely = zeros(T, batchsize)
    I, F = indexbounds(inputlens)
    p = softmax(ᵛ(x); dims=1)
    r = zero(p)

    Threads.@threads for b = 1:batchsize
        span = I[b]:F[b]
        r[:,span], loglikely[b] = CTC(p[:,span], seqlabels[b], blank=blank)
        loglikely[b] /= length(seqlabels[b]) * 2 + 1
    end

    if x.backprop
        function DNN_Batch_CTC_With_Softmax_Backward()
            if need2computeδ!(x)
                if weight==1.0
                    δ(x) .+=  p - r
                else
                    δ(x) .+= (p - r) .* weight
                end
            end
        end
        push!(graph.backward, DNN_Batch_CTC_With_Softmax_Backward)
    end
    return sum(loglikely)/batchsize
end

"""
    RNN_Batch_CTC_With_Softmax(x::Variable{Array{T}}, seqlabels::Vector, inputlens; blank=1, weight=1.0) where T

a batch of padded input sequence is processed by neural networks into `x`

# Inputs
`x`         : 3-D Variable with shape (featdims,timesteps,batchsize), inputs of softmax\n
`seqlabels` : a batch of sequential labels, like [[i,j,k],[x,y],...]\n
`inputlens` : each input's length, like [19,97,...]\n
`weight`    : weight for CTC loss

# Structure
    ┌───┐
    │ │ │
    │ W ├──►─┐
    │ │ │    │
    └───┘    │
    ┌───┐    │    ┌───┐          ┌───┐
    │ │ │  ┌─┴─┐  │ │ │ softmax  │ │ │   ┌───────┐
    │ Z ├─►│ × ├─►│ X ├─────────►│ P ├──►│CTCLOSS│◄── (seqLabel)
    │ │ │  └───┘  │ │ │          │ │ │   └───┬───┘
    └───┘         └─┬─┘          └─┬─┘       │
                    │              │+        ▼
                  ┌─┴─┐            ▼       ┌─┴─┐
                  │ │ │          ┌─┴─┐ -   │ │ │
                  │ δ │◄─────────┤ - │──◄──┤ r │
                  │ │ │          └───┘     │ │ │
                  └───┘                    └───┘
"""
function RNN_Batch_CTC_With_Softmax(x::Variable{Array{T}}, seqlabels::Vector, inputlens; blank=1, weight=1.0) where T
    batchsize = length(inputlens)
    loglikely = zeros(T, batchsize)
    p = zero(ᵛ(x))
    r = zero(ᵛ(x))

    Threads.@threads for b = 1:batchsize
        Tᵇ = inputlens[b]
        Lᵇ = length(seqlabels[b])
        p[:,1:Tᵇ,b] = softmax(x.value[:,1:Tᵇ,b]; dims=1)
        r[:,1:Tᵇ,b], loglikely[b] = CTC(p[:,1:Tᵇ,b], seqlabels[b], blank=blank)
        loglikely[b] /= Lᵇ * 2 + 1
    end

    if x.backprop
        function RNN_Batch_CTC_With_Softmax_Backward()
            if need2computeδ!(x)
                if weight==1.0
                    δ(x) .+=  p - r
                else
                    δ(x) .+= (p - r) .* weight
                end
            end
        end
        push!(graph.backward, RNN_Batch_CTC_With_Softmax_Backward)
    end
    return sum(loglikely)/batchsize
end

"""
    CRNN_Batch_CTC_With_Softmax(x::Variable{Array{T}}, seqlabels::Vector; blank=1, weight=1.0) where T

a batch of padded input sequence is processed by neural networks into `x`

# Inputs
`x`         : 3-D Variable (featdims,timesteps,batchsize), inputs of softmax\n
`seqlabels` : a batch of sequential labels, like [[i,j,k],[x,y],...]\n
`weight`    : weight for CTC loss

# Structure
    ┌───┐
    │ │ │
    │ W ├──►─┐
    │ │ │    │
    └───┘    │
    ┌───┐    │    ┌───┐          ┌───┐
    │ │ │  ┌─┴─┐  │ │ │ softmax  │ │ │   ┌───────┐
    │ Z ├─►│ × ├─►│ X ├─────────►│ P ├──►│CTCLOSS│◄── (seqLabel)
    │ │ │  └───┘  │ │ │          │ │ │   └───┬───┘
    └───┘         └─┬─┘          └─┬─┘       │
                    │              │+        ▼
                  ┌─┴─┐            ▼       ┌─┴─┐
                  │ │ │          ┌─┴─┐ -   │ │ │
                  │ δ │◄─────────┤ - │──◄──┤ r │
                  │ │ │          └───┘     │ │ │
                  └───┘                    └───┘
"""
function CRNN_Batch_CTC_With_Softmax(x::Variable{Array{T}}, seqlabels::Vector; blank=1, weight=1.0) where T
    featdims, timesteps, batchsize = size(x)
    loglikely = zeros(T, batchsize)
    p = softmax(ᵛ(x); dims=1)
    r = zero(ᵛ(x))

    Threads.@threads for b = 1:batchsize
        r[:,:,b], loglikely[b] = CTC(p[:,:,b], seqlabels[b], blank=blank)
        loglikely[b] /= length(seqlabels[b]) * 2 + 1
    end

    if x.backprop
        function CRNN_Batch_CTC_With_Softmax_Backward()
            if need2computeδ!(x)
                if weight==1.0
                    δ(x) .+=  p - r
                else
                    δ(x) .+= (p - r) .* weight
                end
            end
        end
        push!(graph.backward, CRNN_Batch_CTC_With_Softmax_Backward)
    end
    return sum(loglikely)/batchsize
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
using Test, OpenStreetMapX
using Random
import LightGraphs

# include("test_core.jl")
include("test_walkpath.jl")


text:
using Dates

fname, _ = mktemp()

println("The modification time of $fname is ", Dates.unix2datetime(mtime(fname)))
println("\nTouch this file.")
touch(fname)
println("The modification time of $fname is now ", Dates.unix2datetime(mtime(fname)))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
###############################################################################
# This file defines commonly used constraints for power flow models
# These constraints generally assume that the model contains p and q values
# for branches flows and bus flow conservation
###############################################################################

# Generic thermal limit constraint
"`p[f_idx]^2 + q[f_idx]^2 <= rate_a^2`"
function constraint_thermal_limit_from(pm::GenericPowerModel, n::Int, f_idx, rate_a)
    p_fr = pm.var[:nw][n][:p][f_idx]
    q_fr = pm.var[:nw][n][:q][f_idx]
    @constraint(pm.model, p_fr^2 + q_fr^2 <= rate_a^2)
end

"`p[t_idx]^2 + q[t_idx]^2 <= rate_a^2`"
function constraint_thermal_limit_to(pm::GenericPowerModel, n::Int, t_idx, rate_a)
    p_to = pm.var[:nw][n][:p][t_idx]
    q_to = pm.var[:nw][n][:q][t_idx]
    @constraint(pm.model, p_to^2 + q_to^2 <= rate_a^2)
end

"`norm([p[f_idx]; q[f_idx]]) <= rate_a`"
function constraint_thermal_limit_from(pm::GenericPowerModel{T}, n::Int, f_idx, rate_a) where T<:AbstractConicPowerFormulation
    p_fr = pm.var[:nw][n][:p][f_idx]
    q_fr = pm.var[:nw][n][:q][f_idx]
    @constraint(pm.model, norm([p_fr; q_fr]) <= rate_a)
end

"`norm([p[t_idx]; q[t_idx]]) <= rate_a`"
function constraint_thermal_limit_to(pm::GenericPowerModel{T}, n::Int, t_idx, rate_a) where T<:AbstractConicPowerFormulation
    p_to = pm.var[:nw][n][:p][t_idx]
    q_to = pm.var[:nw][n][:q][t_idx]
    @constraint(pm.model, norm([p_to; q_to]) <= rate_a)
end

# Generic on/off thermal limit constraint

"`p[f_idx]^2 + q[f_idx]^2 <= (rate_a * branch_z[i])^2`"
function constraint_thermal_limit_from_on_off(pm::GenericPowerModel, n::Int, i, f_idx, rate_a)
    p_fr = pm.var[:nw][n][:p][f_idx]
    q_fr = pm.var[:nw][n][:q][f_idx]
    z = pm.var[:nw][n][:branch_z][i]
    @constraint(pm.model, p_fr^2 + q_fr^2 <= rate_a^2*z^2)
end

"`p[t_idx]^2 + q[t_idx]^2 <= (rate_a * branch_z[i])^2`"
function constraint_thermal_limit_to_on_off(pm::GenericPowerModel, n::Int, i, t_idx, rate_a)
    p_to = pm.var[:nw][n][:p][t_idx]
    q_to = pm.var[:nw][n][:q][t_idx]
    z = pm.var[:nw][n][:branch_z][i]
    @constraint(pm.model, p_to^2 + q_to^2 <= rate_a^2*z^2)
end

"`p_ne[f_idx]^2 + q_ne[f_idx]^2 <= (rate_a * branch_ne[i])^2`"
function constraint_thermal_limit_from_ne(pm::GenericPowerModel, n::Int, i, f_idx, rate_a)
    p_fr = pm.var[:nw][n][:p_ne][f_idx]
    q_fr = pm.var[:nw][n][:q_ne][f_idx]
    z = pm.var[:nw][n][:branch_ne][i]
    @constraint(pm.model, p_fr^2 + q_fr^2 <= rate_a^2*z^2)
end

"`p_ne[t_idx]^2 + q_ne[t_idx]^2 <= (rate_a * branch_ne[i])^2`"
function constraint_thermal_limit_to_ne(pm::GenericPowerModel, n::Int, i, t_idx, rate_a)
    p_to = pm.var[:nw][n][:p_ne][t_idx]
    q_to = pm.var[:nw][n][:q_ne][t_idx]
    z = pm.var[:nw][n][:branch_ne][i]
    @constraint(pm.model, p_to^2 + q_to^2 <= rate_a^2*z^2)
end

"`pg[i] == pg`"
function constraint_active_gen_setpoint(pm::GenericPowerModel, n::Int, i, pg)
    pg_var = pm.var[:nw][n][:pg][i]
    @constraint(pm.model, pg_var == pg)
end

"`qq[i] == qq`"
function constraint_reactive_gen_setpoint(pm::GenericPowerModel, n::Int, i, qg)
    qg_var = pm.var[:nw][n][:qg][i]
    @constraint(pm.model, qg_var == qg)
end

"""
Creates Line Flow constraint for DC Lines (Matpower Formulation)

```
p_fr + p_to == loss0 + p_fr * loss1
```
"""
function constraint_dcline(pm::GenericPowerModel{T}, n::Int, f_bus, t_bus, f_idx, t_idx, loss0, loss1) where T
    p_fr = pm.var[:nw][n][:p_dc][f_idx]
    p_to = pm.var[:nw][n][:p_dc][t_idx]

    @constraint(pm.model, (1-loss1) * p_fr + (p_to - loss0) == 0)
end

"`pf[i] == pf, pt[i] == pt`"
function constraint_active_dcline_setpoint(pm::GenericPowerModel, n::Int, f_idx, t_idx, pf, pt)
    p_fr = pm.var[:nw][n][:p_dc][f_idx]
    p_to = pm.var[:nw][n][:p_dc][t_idx]

    @constraint(pm.model, p_fr == pf)
    @constraint(pm.model, p_to == pt)
end


text:
charlist = 'a':'z'

function randomstring(n::Integer)
    str = ""
    for i = 1:n
        str *= rand(charlist)
    end
    return str
end


text:
################################################################################
#
#	CONCRETE TYPE
#
#   ReciprocalPath{P} <: AbstractReciprocalPath{P}
#   --> P is the reciprocal point type (<: AbstractReciprocalPoint{D})
#       --> D is the dimension of embedding space
#
#   FILE CONTAINS
#       - concrete struct definition
#       - interface implementation
#
################################################################################

################################################################################
#
#   STRUCT DEFINITION
#
################################################################################
mutable struct ReciprocalPath{P} <: AbstractReciprocalPath{P}

    # the only field is the array of points
    points :: Vector{P}

end

# export the type
export ReciprocalPath

################################################################################
#
#	INTERFACING / ACCESSING RECIPROCAL PATHS
#
################################################################################

# default constructor interface
# used for creation of new reciprocal paths
function newReciprocalPath(
            ::Type{ReciprocalPath{RPO}},
            points :: Vector{<:RPO}
        ) :: ReciprocalPath{RPO} where {D,RPO<:AbstractReciprocalPoint{D}}

    # return the newly created object
    return ReciprocalPath{RPO}(points)
end

# export the newUnitcell interface
export newReciprocalPath

# access the list of points
function points(
            p :: ReciprocalPath{RPO}
        ) :: Vector{RPO} where {D,RPO<:AbstractReciprocalPoint{D}}

    # return the respective field
    return p.points
end
# set the list of points
function points!(
            p :: ReciprocalPath{RPO},
            points :: Vector{<:RPO}
        ) where {D,RPO<:AbstractReciprocalPoint{D}}

    # set the respective field
    p.points = points
end

# export access to points
export points, points!

function saveReciprocalPath(
        rp :: RP,
        fn :: AbstractString,
        group :: AbstractString = "reciprocal_path"
        ;
        append :: Bool = false
    ) where {D,P<:ReciprocalPoint{D},RP<:ReciprocalPath{P}}

    # determine the mode based on if one wants to append stuff
    if append
        mode = "r+"
    else
        mode = "w"
    end

    # open the file in mode
    h5open(fn, mode) do file
        # create the group in which the bonds are saved
        group_path = g_create(file, group)
        # save the parameters
        attrs(group_path)["D"] = Int64(D)
        # save all Positions (D dimensions)
        if 0 < Int64(D)
            for n in 1:Int64(D)
                group_path["point_$(n)"] = Float64[point(p)[n] for p in rp]
            end
        end
        # save all labels
        group_path["label"]       = String[string(label(p)) for p in rp]
        group_path["label_LaTeX"] = String[string(labelLaTeX(p)) for p in rp]
    end

    # return nothing
    return nothing
end

function loadReciprocalPath(
        ::Type{RPI},
        fn :: AbstractString,
        group :: AbstractString = "reciprocal_path"
    ) where {DI,PI<:ReciprocalPoint{DI},RPI<:Union{ReciprocalPath{PI},ReciprocalPath}}

    # read attribute data
    attr_data = h5readattr(fn, group)
    # determine D based on this
    D  = attr_data["D"]

    # load all remaining data
    rp_label       = h5read(fn, group*"/label")
    rp_label_LaTeX = h5read(fn, group*"/label_LaTeX")

    if D == 0
        rp_point = [Float64[] for i in 1:length(rp_label)]
    else
        rp_point_parts = [
            h5read(fn, group*"/point_"*string(j)) for j in 1:D
        ]
        rp_point = Vector{Float64}[Float64[rp_point_parts[j][i] for j in 1:D] for i in 1:length(rp_label)]
    end

    # create list of points in the path
    rpp = ReciprocalPoint{D}[
        newReciprocalPoint(ReciprocalPoint{D}, rp_point[i], rp_label[i], rp_label_LaTeX[i])
        for i in 1:length(rp_label)
    ]

    # return the path
    return newReciprocalPath(ReciprocalPath{ReciprocalPoint{D}}, rpp)
end

# convinience function for standard type
function loadReciprocalPath(
        fn :: AbstractString,
        group :: AbstractString = "reciprocal_path"
    )

    return loadReciprocalPath(ReciprocalPath, fn, group)
end


text:
mutable struct StructuredOutputFile{RT}
    B::Block{RT}
    CODE::Dict{String,TCode}
    CODEINV::Dict{TCode,String}
    LINES::Vector{String}
    TOKLINES::Vector{String}
end

# returns the
# B, # tree-representation of the data
# code, codeinv, # word-token and token-word maps
# lines,  # original file
# lines1  # tokenized file

function structurize_file(fn::String)::StructuredOutputFile{RTYPE}
    @assert isfile(fn)
    lines  = String[]
    lines1 = String[]
    lines, lines1, patts, code = load_file(fn)
    codeinv = revert(code) ;
    if length(lines)==0
        B = Multiline()
    else
        B = (   build_block_init(patts)
                    |> find_block_MostFreqSimilarSubsq
                    |> find_block_MostFreqSimilarSubsq
                    |> find_block_MostFreqSimilarSubsq
                    |> find_block_MostFreqSimilarSubsq
                    |> find_block_MostFreqSimilarSubsq
                    |> find_block_MostFreqSimilarSubsq
                    |> find_block_MostFreqSimilarSubsq
                    |> find_block_MostFreqSimilarSubsq   );
        parse_file!(0, B, lines, lines1, codeinv) ;
        @assert is_valid_block(B)
    end
    return StructuredOutputFile{RTYPE}(B, code, codeinv, lines, lines1)
end

max_stars_repo_name: algorithmx/SmartParser

max_issues_repo_name: algorithmx/SmartParser

max_forks_repo_name: algorithmx/SmartParser


text:
# Interface -------

function collaborators(owner, repo; auth = AnonymousAuth(), options...)
    collaborators(auth, owner, repo; options...)
end

function collaborators(auth::Authorization, owner, repo; headers = Dict(), options...)
    authenticate_headers(headers, auth)
    r = get(URI(API_ENDPOINT; path = "/repos/$owner/$repo/collaborators");
            headers = headers,
            options...)

    handle_error(r)

    data = JSON.parse(r.data)
    map!( u -> User(u), data)
end

function iscollaborator(owner, repo, user; auth = AnonymousAuth(), options...)
    iscollaborator(auth, owner, repo, user; options...)
end

function iscollaborator(auth::Authorization, owner, repo, user; headers = Dict(), options...)
    authenticate_headers(headers, auth)
    r = get(URI(API_ENDPOINT; path = "/repos/$owner/$repo/collaborators/$user");
            headers = headers,
            options...)

    r.status == 204 && return true
    r.status == 404 && return false

    handle_error(r)  # 404 is not an error in this case

    return false  # at this point, assume no
end

function add_collaborator(owner, repo, user; auth = AnonymousAuth(), options...)
    add_collaborator(auth, owner, repo, user; options...)
end

function add_collaborator(auth::Authorization, owner, repo, user; headers = Dict(), options...)
    authenticate_headers(headers, auth)
    r = put(URI(API_ENDPOINT; path = "/repos/$owner/$repo/collaborators/$user");
            headers = headers,
            options...)

    handle_error(r)
end

function remove_collaborator(owner, repo, user; auth = AnonymousAuth(), options...)
    remove_collaborator(auth, owner, repo, user; options...)
end

function remove_collaborator(auth::Authorization, owner, repo, user; headers = Dict(), options...)
    authenticate_headers(headers, auth)
    r = delete(URI(API_ENDPOINT; path = "/repos/$owner/$repo/collaborators/$user");
            headers = headers,
            options...)

    handle_error(r)
end


text:
## An example of how to add points, lines to a Graphics Scene and have them displayed in a GraphicsView.
## This follows some of the code of the qtpaint package from the R world by Michael Lawrence

using PySide

## copy qtpaint style
## We return graphics items generated
graphic_item_flag(item, flags) = PyCall.pyeval(join(["item.$u" for u in flags], " | "), item=item)
set_flags(item, flags) = item[:setFlags](graphic_item_flag(item, flags))

## Two constructors, on vectorized
function addPoint(sc::GraphicsScene, xs::Vector, ys::Vector, r;
                  lwd::Integer=1, pen_style::String="SolidLine", color::String="Black")
    function add_point(x, y, r, pen, brush)
        item = sc[:addEllipse](0, 0, r, r, pen, brush)
        item[:setPos](x, y)
        set_flags(item, ["ItemIsSelectable", "ItemIgnoresTransformations"])
        item
    end
    pen = Qt[:QPen](Qt[:QColor](color))
    pen[:setStyle](qt_enum(pen_style))
    pen[:setWidth](lwd)
    brush = Qt[:QBrush](qt_enum("SolidPattern"))
    [add_point(xs[i], ys[i], r, pen, brush) for i in 1:length(xs)]
end

function addLine(sc::GraphicsScene, xs, ys;
                  lwd::Integer=1, pen_style::String="SolidLine", color::String="Black")
    pen = Qt[:QPen](Qt[:QColor](color))
    pen[:setStyle](qt_enum(pen_style))
    pen[:setWidth](lwd)
    brush = Qt[:QBrush](qt_enum("SolidPattern"))
    [sc[:addLine](xs[i-1], ys[i-1], xs[i], ys[i], pen, brush) for i in 2:length(x)]
end

function addSegment(sc::GraphicsScene,  xs, ys, x1s, y1s;
                     lwd::Integer=1, pen_style::String="SolidLine", color::String="Black")
    pen = Qt[:QPen](Qt[:QColor](color))
    pen[:setStyle](qt_enum(pen_style))
    pen[:setWidth](lwd)
    brush = Qt[:QBrush](qt_enum("SolidPattern"))
    [sc[:addLine](xs[i], ys[i], x1s[i], y1s[i], pen, brush) for i in 1:length(xs)]
end

function addPolygon(sc::GraphicsScene, xs, ys,
                     lwd::Integer=1, pen_style::String="SolidLine", color::String="Black")
    pen = Qt[:QPen](Qt[:QColor](color))
    pen[:setStyle](qt_enum(pen_style))
    pen[:setWidth](lwd)
    brush = Qt[:QBrush]()

    pts = [QtCore.QPoint(xs[i], ys[i]) for i in 1:length(xs)]
    item = sc[:addPolygon](pts, pen, brush)
   item
end

function addRect(sc::GraphicsScene,  xs, ys, ws, hs,
                 lwd::Integer=1, pen_style::String="SolidLine", color::String="Black")
    pen = Qt[:QPen](Qt[:QColor](color))
    pen[:setStyle](qt_enum(pen_style))
    pen[:setWidth](lwd)
    brush = Qt[:QBrush](qt_enum("SolidPattern"))

    [sc[:addRect](x[i], y[i], w[i], h1[i], pen, brush) for i in 1:length(xs)]
end

function addText(sc::GraphicsScene, xs, ys, labels, html::Bool; font::String="Arial")
    function add_text(sc, x, y, label, html, font)
        item = sc[:addText]("")
        item[:setFont](Qt.QFont(font))
        if html
            item[:setHtml](label)
            item[:setOpenExternalLinks](true)
        else
            item[:setPlainText](label)
        end
        item[:setPos](x, y)
        set_flags(item, ["ItemIsMovable", "ItemIsSelectable", "QGraphicsItem::ItemIsFocusable", "ItemIgnoresTransformations"])
        item
    end
    [add_text(sc, xs[i], ys[i], labels[i], html, font) for i in 1:length(xs)]
end

## Return (view, scene)
function PlotView(parent)
    qnew_class("PlotViewWidget", "QtGui.QGraphicsView")

    o = qnew_class_instance("PlotViewWidget")
    qinvoke(o, :setParent, parent)
    scene = GraphicsScene(parent)

    qinvoke(scene, :setSceneRect, 0,0, 1000, 1000)
    qinvoke(o, :setScene, scene)

    ## qtpaint has other options for resizing that may be better suited.
    function resize_event(event)
        osize = qinvoke(event, :oldSize)
        size =  qinvoke(o, [:viewport, :size])

        owidth = qinvoke(osize, :width); oheight = qinvoke(osize, :height)
        width = qinvoke(size, :width);   height = qinvoke(size, :height)

        if owidth < 0
            qinvoke(o, :fitInView, qinvoke(o, :sceneRect))
        else
            sx = width / owidth
            sy = height / oheight

            anchor = qinvoke(o, :transformationAnchor)
            qinvoke(o, :setTransformationAnchor, o[:NoAnchor])
            qinvoke(o, :scale, sx, sy)
            qinvoke(o, :setTransformationAnchor, anchor)

        end
    end
    o[:resizeEvent] = resize_event

    (o, scene)
end

## Some constants, best to be elsewhere
    xpx = 1000
    ypx = 1000
    xlim = [0, 10]
    ylim = [0, 10]

function scale_xy_px(xs, ys)
    width = xpx; height = ypx
    newxs = width*(xs - xlim[1])/(xlim[2] - xlim[1])
    newys = height - height * (ys - ylim[1]) / (ylim[2] - ylim[1])
    ([newxs], [newys])                  #  make vectors
end

    ## R-like interface to adding points and lines
function points(scene, xs, ys)
    newxs, newys = scale_xy_px(xs, ys)
    addPoint(scene, newxs, newys, 5)
end

function lines(scene, xs, ys)
    newxs, newys = scale_xy_px(xs, ys)
    addLine(scene, newxs, newys, 4)
end

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

## https://groups.google.com/forum/?fromgroups=#!topic/julia-users/g9dp5fypJ88
## Example of stepping through graphics
t = linspace(0, 4pi, 512);
x = sin(t);
y = 0.1t.*cos(t);

xlim = [minimum(x), maximum(x)] * 1.05
ylim = [minimum(y), maximum(y)] * 1.05

w = Widget()
view, scene = PlotView(w)
push!(w, view)
raise(w)

i = 1
timer = QtCore[:QTimer]()
qconnect(timer, :timeout) do
global i
if i == 512
    qinvoke(timer, :stop)
else
    i = i + 1
    points(scene, x[i], y[i])
end
end
qinvoke(timer, :start, 20)


text:
using IntroJuliaDemo
using Test

@testset "IntroJuliaDemo.jl" begin
    # Write your tests here.

    # Pruebas de valores
    @test f(2,3) == 13
    @test f(2, 3.0) ≈ 13.0

    # Pruebas de tipos
    @test typeof(f(2, 3.0)) == Float64
end


text:
using Distributions
using PyPlot
using LsqFit
using CSV, DataFrames, DataFramesMeta

""" Functions for results of Fig. 2D """
# Run:
# graph_1_visits, graph_3_visits, graph_1_subjIDs, graph_3_subjIDs = SinergiafMRI_datafit.get_state_visits()
# and then (to process and save results in csv):
# SinergiafMRI_datafit.csv_state_decision_stats(1, graph_1_visits)
# SinergiafMRI_datafit.csv_state_decision_stats(3, graph_3_visits)
#
# and (for quick plotting) :
# SinergiafMRI_datafit.plot_state_decision_stats(1, graph_1_visits)
# SinergiafMRI_datafit.plot_state_decision_stats(3, graph_3_visits)
function get_state_visits()
    # Real:
    data_path = "./mcmc_rl_fit/projects/fmri/data/SARSPEICZVG_all_fMRI.csv"
    # Simulated:
    # data_path = "./mcmc_rl_fit/projects/fmri/data_sim/SARSPEICZVG_all_fMRI_Sim_exporder_SurpAC_run2.csv"

    max_episodes_per_subj = 40
    (all_data, all_stats) = RL_Fit.load_SARSPEI(data_path;max_episodes_per_subj=max_episodes_per_subj)
    all_subject_ids = all_stats.all_subject_ids
    println("all_subject_ids: $all_subject_ids")

    graph_1_visits = zeros(Int64, max_episodes_per_subj,7, 2)
    graph_3_visits = zeros(Int64, max_episodes_per_subj,7, 2)
    graph_1_subjIDs = []
    graph_3_subjIDs = []
    for subj_id in all_subject_ids
        @show subj_id
        subj_data = RL_Fit.filter_by_person(all_data, subj_id)
        graph_id = subj_data[1,Int(RL_Fit.c_graph_version)]
        all_episodes = unique(subj_data[:,Int(RL_Fit.c_episode)])
        trajectory_count_per_state = zeros(Int64, 7) #count for each state individually how many times it participated in a rewarded episode
        # per subject, do not count (sum), just set to 0 or 1.
        state_action_indicator = zeros(Int64, max_episodes_per_subj, 7, 2)

        for episode in all_episodes
#             print("\nS:$subj_id/E:$episode ")
            for s in 2:all_stats.max_state_id
                # check if state s has "contributed" to the current episode. if yes, count this trajectory.
                if any( (subj_data[:,Int(RL_Fit.c_episode)] .== episode) .& (subj_data[:,Int(RL_Fit.c_state)] .== s) )
                    trajectory_count_per_state[s] += 1
                    #             end
                    # now count the value of actions. distinguish two cases:
                    # for some states, the two actions are equally valid. In this case count the action
                    # for other states, there is a correct (1) and a wrong(-1) action
                    # action correct:
                    if any( (subj_data[:,Int(RL_Fit.c_episode)] .== episode) .& (subj_data[:,Int(RL_Fit.c_state)] .== s) .& (subj_data[:,Int(RL_Fit.c_action_value)] .== 0))
                        # action A or action B?
                        if any( (subj_data[:,Int(RL_Fit.c_episode)] .== episode) .& (subj_data[:,Int(RL_Fit.c_state)] .== s) .& (subj_data[:,Int(RL_Fit.c_action)] .== 1))
                            # A
                            state_action_indicator[trajectory_count_per_state[s], s, 1] = 1 # do not +=1
                            #                         print("A")
                        end
                        # note: it's not an ELSE here. we check for both and count both options at most once
                        if any( (subj_data[:,Int(RL_Fit.c_episode)] .== episode) .& (subj_data[:,Int(RL_Fit.c_state)] .== s) .& (subj_data[:,Int(RL_Fit.c_action)] .== 2))
                            # B
                            state_action_indicator[trajectory_count_per_state[s], s, 2] = 1 # do not +=1
                            #                         print("B")
                        end
                    else
                        # action correct or wrong?
                        if any( (subj_data[:,Int(RL_Fit.c_episode)] .== episode) .& (subj_data[:,Int(RL_Fit.c_state)] .== s) .& (subj_data[:,Int(RL_Fit.c_action_value)] .== +1))
                            # correct
                            state_action_indicator[trajectory_count_per_state[s], s, 2] = 1 # do not +=1
                            #                         print("c")
                        end
                        # note: it's not an ELSE here. we check for both and count both options at most once
                        if any( (subj_data[:,Int(RL_Fit.c_episode)] .== episode) .& (subj_data[:,Int(RL_Fit.c_state)] .== s) .& (subj_data[:,Int(RL_Fit.c_action_value)] .== -1))
                            # wrong
                            state_action_indicator[trajectory_count_per_state[s], s, 1] = 1 # do not +=1
                            #                         print("w")
                        end
                    end
                else
#                     print("N")
                end
            end # end all states
        end # end all episodes
        if graph_id == 1
            graph_1_visits += state_action_indicator
            push!(graph_1_subjIDs, subj_id)
        elseif graph_id == 3
            graph_3_visits += state_action_indicator
            push!(graph_3_subjIDs, subj_id)
        else
            error("unknown graph id: $graph_id")
        end
    end # end all subject
    graph_1_visits, graph_3_visits, graph_1_subjIDs, graph_3_subjIDs
end

model_exp(episode, param) = 1. .- param[1] * exp.(-Float64.(episode).*param[2])
# model_exp(episode, param) = 1- param[1]*exp.(-episode.*param[2])
model_lin(episode, param) = param[1] .+ episode .* param[2]
inv_model_exp(val, param) = (-1/param[2])*log((1-val)/param[1])

function get_this_graph_info(graph_info, graph_id)
    if graph_id == 1
        states = graph_info["state_order_g1"]
        minimal_dist = graph_info["minimal_dist_g1"]
        state_txt = graph_info["state_txt_g1"]
        graph = graph_info["transition_matrix_g1"]
    else
        states = graph_info["state_order_g3"]
        minimal_dist = graph_info["minimal_dist_g3"]
        state_txt = graph_info["state_txt_g3"]
        graph = graph_info["transition_matrix_g3"]
    end
    states, minimal_dist, state_txt, graph
end
function plot_state_decision_stats(graph_id::Int,
                            graph_state_counts;
                            learning_threshold=0.8,
                            min_visits_requirement = 2)

    save_path = makehomesavepath("learningcurves")
    mkdir(save_path)

    graph_info = get_graph_info()
    states, minimal_dist, state_txt, graph = get_this_graph_info(graph_info, graph_id)

    y_min_lim = 0.25
    init_vals_lin = [0.5, 0.01]
    init_vals_exp = [0.5, 0.]
    fig, axes = subplots(nrows=6, ncols=1, figsize=(8.5, 11.5))
    ax_idx = 0
    for s = states
        println(s)
        ax_idx += 1

        counts_per_episode_actionwrong = vec(graph_state_counts[:,s,1])
        counts_per_episode_actioncorrect = vec(graph_state_counts[:,s,2])
        summed_counts_per_episode = counts_per_episode_actionwrong.+counts_per_episode_actioncorrect
        episodes_min_visits = findall(summed_counts_per_episode .>= min_visits_requirement)
        ratio = counts_per_episode_actioncorrect[episodes_min_visits] ./ summed_counts_per_episode[episodes_min_visits]
        weight = (summed_counts_per_episode[episodes_min_visits])
        nr_x = length(ratio)

        x_data = 1:nr_x
        y_data = ratio
        marker_symbol = (minimal_dist[1,s] == minimal_dist[2,s]) ? "oy" : "og"
        for i in 1:nr_x #plot does not accept a vector for size.
            axes[ax_idx].plot(x_data[i], ratio[i], marker_symbol,
                                            markersize = 1+2*sqrt(weight[i]),
                                            markeredgewidth=0.0)
        end
        weight = 1.0 .*weight  # just one example why I sometimes do not like JuliaLang
        # weight = Float64.(weight)
        if minimal_dist[1,s] == minimal_dist[2,s]
            fit = curve_fit(model_lin, x_data, y_data, weight, init_vals_lin)
            y_hat = model_lin(x_data, fit.param)
            axes[ax_idx].plot(x_data, y_hat, "-r", lineWidth=2)
            axes[ax_idx].set_title("G: $graph_id | S: $(state_txt[s]) | dist:($(minimal_dist[1,s])/$(minimal_dist[2,s]))| LIN fit:$(round.(fit.param, digits=2))", fontsize=10)
            axes[ax_idx].set_yticks((0,1))
            axes[ax_idx].set_yticklabels(("A", "B"))
            axes[ax_idx].set_ylabel("Decision", fontsize=10)
            axes[ax_idx].set_ylim([-0.05, 1.05])

        else
            fit = curve_fit(model_exp, x_data, y_data, weight, init_vals_exp) # Weighted cost function
            y_hat = model_exp(x_data, fit.param)
            axes[ax_idx].plot(x_data, y_hat, "-r", lineWidth=2)
            threshold_delay = inv_model_exp(learning_threshold, fit.param)
            @show s
            @show fit.param
            @show threshold_delay
            if -10< threshold_delay && threshold_delay <= 50
                line_s = (fit.param[2] >0) ? "-" : "--"
                line_c = (fit.param[2] >0) ? "b" : "c"
                axes[ax_idx].axvline(threshold_delay, ls=line_s, c=line_c, lineWidth=2)
            end
            axes[ax_idx].set_title("G: $graph_id | S: $(state_txt[s]) | dist:($(minimal_dist[1,s])/$(minimal_dist[2,s]))| EXP fit:$(round.(fit.param, digits=2)) | threshold $learning_threshold at $(round(threshold_delay, digits=1))", fontsize=10)
            axes[ax_idx].set_yticks((0.5, 1.0))
            axes[ax_idx].set_ylabel("correct", fontsize=10)
            axes[ax_idx].set_ylim([y_min_lim, 1.05])
        end
        axes[ax_idx].set_xlabel("n-th trajectory", fontsize=10)
        axes[ax_idx].set_xlim([0,35])
#         break
    end
        fig.subplots_adjust(hspace=1.1)
        fig.savefig(joinpath(save_path, "learning_curves_G$graph_id.png"), dpi = 600)
        fig.savefig(joinpath(save_path, "learning_curves_G$graph_id.pdf"), dpi = 600)
    return axes
end

function csv_state_decision_stats(graph_id::Int,
                            graph_state_counts;
                            learning_threshold=0.8,
                            min_visits_requirement = 2,
                            # filestr = "_Sim",
                            filestr = "",
                            # savepath = "./",
                            )

    savepath = makehomesavepath("learningcurves")
    mkdir(savepath)

    dfresults = []
    graph_info = get_graph_info()
    states, minimal_dist, state_txt, graph = get_this_graph_info(graph_info, graph_id)

    init_vals_lin = [0.5, 0.01]
    init_vals_exp = [0.5, 0.]
    dfresults_threshold_all = DataFrame()
    for s = states
        dfresults = DataFrame()
        println(s)
        counts_per_episode_actionwrong = vec(graph_state_counts[:,s,1])
        counts_per_episode_actioncorrect = vec(graph_state_counts[:,s,2])
        summed_counts_per_episode = counts_per_episode_actionwrong.+counts_per_episode_actioncorrect
        episodes_min_visits = findall(summed_counts_per_episode .>= min_visits_requirement)
        ratio = counts_per_episode_actioncorrect[episodes_min_visits] ./ summed_counts_per_episode[episodes_min_visits]
        weight = (summed_counts_per_episode[episodes_min_visits])
        nr_x = length(ratio)

        x_data = 1:nr_x
        y_data = ratio
        # markersize = 1+2*sqrt(weight[i]),
        threshold_delay = 0.
        weight = Float64.(weight)
        if minimal_dist[1,s] == minimal_dist[2,s]
            fit = curve_fit(model_lin, x_data, y_data, weight, init_vals_lin)
            y_hat = model_lin(x_data, fit.param)
            # axes[ax_idx].set_title("G: $graph_id | S: $(state_txt[s]) | dist:($(minimal_dist[1,s])/$(minimal_dist[2,s]))| LIN fit:$(round.(fit.param, digits=2))", fontsize=10)
        else
            fit = curve_fit(model_exp, x_data, y_data, weight, init_vals_exp) # Weighted cost function
            y_hat = model_exp(x_data, fit.param)
            threshold_delay = inv_model_exp(learning_threshold, fit.param)
            if threshold_delay > episodes_min_visits[end]
                threshold_delay = episodes_min_visits[end]
            end
            # axes[ax_idx].set_title("G: $graph_id | S: $(state_txt[s]) | dist:($(minimal_dist[1,s])/$(minimal_dist[2,s]))| EXP fit:$(round.(fit.param, digits=2)) | threshold $learning_threshold at $(round(threshold_delay, digits=1))", fontsize=10)
        end
        @show s
        @show threshold_delay
        dfresults[!, Symbol(string(s)*"_"*state_txt[s]*"_episodenumber")] = x_data
        dfresults[!, Symbol(string(s)*"_"*state_txt[s]*"_ratiocorrect")] = y_data
        dfresults[!, Symbol(string(s)*"_"*state_txt[s]*"_nrofsubjects")] = weight
        dfresults[!, Symbol(string(s)*"_"*state_txt[s]*"_fitteddata")] = y_hat
        dfresults[!, Symbol(string(s)*"_"*state_txt[s]*"_param1")] .= [fit.param[1]]
        dfresults[!, Symbol(string(s)*"_"*state_txt[s]*"_param2")] .= [fit.param[2]]
        dfresults[!, Symbol(string(s)*"_"*state_txt[s]*"_minimaldist")] .= ["$(minimal_dist[1,s])/$(minimal_dist[2,s])"]
        dfresults[!, Symbol(string(s)*"_"*state_txt[s]*"_thresholddelay")] .= [threshold_delay]
        dfresults[!, Symbol(string(s)*"_"*state_txt[s]*"_graphid")] .= [graph_id]

        CSV.write(joinpath(savepath, "learningcurve_G" * string(graph_id) * "_" * string(s) * "_" * state_txt[s] * filestr * ".csv"), dfresults,
            delim = " ");

        min_dist = sort(minimal_dist[:,s] .- 1.)
        dfresults_threshold_all[!, Symbol(string(s)*"_"*state_txt[s]*"_thresholddelay")] = [threshold_delay]
        dfresults_threshold_all[!, Symbol(string(s)*"_"*state_txt[s]*"_dist1vsdist2")] = [string(Int(min_dist[1]))*"-or-"*string(Int(min_dist[2]))]
        dfresults_threshold_all[!, Symbol(string(s)*"_"*state_txt[s]*"_graphid")] = [graph_id]

        CSV.write(joinpath(savepath, "learningcurvebar_true_G" * string(graph_id) * filestr * ".csv"), dfresults_threshold_all,
            delim = " ");
    end
end

# axes_1, deltaV_vs_Lerning1 = plot_state_decision_stats(1, graph_1_visits, state_txt_g1, graph_1_minimal_dist, transition_matrix_g1,Vvalues_g1)
# axes_3, deltaV_vs_Lerning3 = plot_state_decision_stats(3, graph_3_visits, state_txt_g3, graph_3_minimal_dist, transition_matrix_g3,Vvalues_g3)

# NOTE: NOT USED!!!
function csv_state_decision_stats_nofitting(graph_id::Int,
                            graph_state_counts;
                            learning_threshold=0.8,
                            min_visits_requirement = 2,
                            savepath = "/Users/vasia/Desktop/")
    dfresults = []
    graph_info = get_graph_info()
    states, minimal_dist, state_txt, graph = get_this_graph_info(graph_info, graph_id)

    dfresults_threshold_all = DataFrame()
    for s = states
        dfresults = DataFrame()
        println(s)
        counts_per_episode_actionwrong = vec(graph_state_counts[:,s,1])
        counts_per_episode_actioncorrect = vec(graph_state_counts[:,s,2])
        summed_counts_per_episode = counts_per_episode_actionwrong.+counts_per_episode_actioncorrect
        episodes_min_visits = findall(summed_counts_per_episode .>= min_visits_requirement)
        ratio = counts_per_episode_actioncorrect[episodes_min_visits] ./ summed_counts_per_episode[episodes_min_visits]
        weight = (summed_counts_per_episode[episodes_min_visits])
        nr_x = length(ratio)

        x_data = 1:nr_x
        y_data = ratio
        # markersize = 1+2*sqrt(weight[i]),
        threshold_delay = 0.

        if minimal_dist[1,s] != minimal_dist[2,s]
            threshold_delay = findfirst(x-> x.>=0.8, ratio)
        end
        @show threshold_delay
        if isnothing(threshold_delay)
            threshold_delay = Float64(episodes_min_visits[end])
        end
        @show s
        @show threshold_delay
        dfresults[!, Symbol(string(s)*"_"*state_txt[s]*"_episodenumber")] = x_data
        dfresults[!, Symbol(string(s)*"_"*state_txt[s]*"_ratiocorrect")] = y_data
        dfresults[!, Symbol(string(s)*"_"*state_txt[s]*"_nrofsubjects")] = weight
        dfresults[!, Symbol(string(s)*"_"*state_txt[s]*"_minimaldist")] .= ["$(minimal_dist[1,s])/$(minimal_dist[2,s])"]
        dfresults[!, Symbol(string(s)*"_"*state_txt[s]*"_thresholddelay")] .= [threshold_delay]
        dfresults[!, Symbol(string(s)*"_"*state_txt[s]*"_graphid")] .= [graph_id]

        CSV.write(joinpath(savepath, "learningcurve_true_nofitting_G" * string(graph_id) * "_" * string(s) * "_" * state_txt[s] *".csv"), dfresults,
            delim = " ");

        min_dist = sort(minimal_dist[:,s] .- 1.)
        dfresults_threshold_all[!, Symbol(string(s)*"_"*state_txt[s]*"_thresholddelay")] = [threshold_delay]
        dfresults_threshold_all[!, Symbol(string(s)*"_"*state_txt[s]*"_dist1vsdist2")] = [string(Int(min_dist[1]))*"-or-"*string(Int(min_dist[2]))]
        dfresults_threshold_all[!, Symbol(string(s)*"_"*state_txt[s]*"_graphid")] = [graph_id]

        CSV.write(joinpath(savepath, "learningcurvebar_true_nofitting_G" * string(graph_id) *".csv"), dfresults_threshold_all,
            delim = " ");
    end
end


text:
using Distributions
using Plots
pyplot()

include("gen_data.jl")
include("particle_filters.jl")

"""Experimental parameters"""

# Visualization options
viz_particles = true

# Number of particles
N = 20

# Time horizon
T = 100

# Experimental parameters
transition = 1.0
emission = 1.0
process_noise = 0.9
measurement_noise = 0.3

# Prior state
prior = [0, 1]

# Generate signal
observations, latent_states = LGDS(transition,
                                   emission,
                                   process_noise,
                                   measurement_noise,
                                   prior;
                                   time_horizon=T)

# Check signal visually
scatter(1:T, observations, color="blue", label="observations")
plot!(1:T, latent_states[2:end], color="red", label="latent states")

"""SIS filter"""

# Nominal distributions
p_transition(a) = Normal(a, sqrt(process_noise))
p_likelihood(b) = Normal(b, sqrt(measurement_noise))

# Call SIS filter
samples, weights = sequential_importance_sampling(observations,
                                                  p_transition,
                                                  p_likelihood,
                                                  num_particles=N,
                                                  imp_processp=inv(process_noise),
                                                  imp_measurementp=inv(measurement_noise))

# Generate posterior distribution
posterior_mean = sum(weights .* samples, dims=2)
posterior_var = sum(weights .* (broadcast(-, samples, posterior_mean)).^2, dims=2)

# Visualize particles
scatter(1:T, observations, color="blue", label="observations")
plot!(1:T, latent_states[2:end], color="red", label="latent states")
plot!(1:T, posterior_mean[2:end], color="purple", label="inferred")
plot!(1:T, posterior_mean[2:end],
      ribbon=[sqrt.(posterior_var[2:end]), sqrt.(posterior_var[2:end])],
      color="purple", alpha=0.1, label="")
xlabel!("time (t)")
ylabel!("signal")
savefig(pwd()*"/viz/sis.png")

if viz_particles

    # Show sample spreads
    scatter(1:T, observations, color="blue", label="observations")
    for t = 2:T
        scatter!((t-1).*ones(1,N), samples[t,:], markersize=weights[t,:]*50, markercolor="purple", label="")
    end
    savefig(pwd()*"/viz/sis_particles.png")
end

"""SIR filter"""

# Nominal distributions
p_transition(a) = Normal(a, sqrt(process_noise))
p_likelihood(b) = Normal(b, sqrt(measurement_noise))

# Call filter
samples, weights = sequential_importance_resampling(observations,
                                                    p_transition,
                                                    p_likelihood,
                                                    num_particles=N,
                                                    imp_processp=inv(process_noise),
                                                    imp_measurementp=inv(measurement_noise))

# Generate posterior distribution
posterior_mean = sum(weights .* samples, dims=2)
posterior_var = sum(weights .* (broadcast(-, samples, posterior_mean)).^2, dims=2)

# Visualize particles
scatter(1:T, observations, color="blue", label="observations")
plot!(1:T, latent_states[2:end], color="red", label="latent states")
plot!(1:T, posterior_mean[2:end], color="purple", label="inferred")
plot!(1:T, posterior_mean[2:end],
    ribbon=[sqrt.(posterior_var[2:end]), sqrt.(posterior_var[2:end])],
    color="purple", alpha=0.1, label="")
xlabel!("time (t)")
ylabel!("signal")
savefig(pwd()*"/viz/sir.png")

if viz_particles

    # Show sample spreads
    scatter(1:T, observations, color="blue", label="observations")
    for t = 2:T
        scatter!((t-1).*ones(1,N), samples[t,:], markersize=weights[t,:]*50, markercolor="purple", label="")
    end
    savefig(pwd()*"/viz/sir_particles.png")
end

"""Bootstrap filter"""

# Nominal distributions
p_transition(a) = Normal(a, sqrt(process_noise))
p_likelihood(b) = Normal(b, sqrt(measurement_noise))

# Call filter
samples, weights = bootstrap_filter(observations,
                                    p_transition,
                                    p_likelihood,
                                    num_particles=N)

# Generate posterior distribution
posterior_mean = sum(weights .* samples, dims=2)
posterior_var = sum(weights .* (broadcast(-, samples, posterior_mean)).^2, dims=2)

# Visualize particles
scatter(1:T, observations, color="blue", label="observations")
plot!(1:T, latent_states[2:end], color="red", label="latent states")
plot!(1:T, posterior_mean[2:end], color="purple", label="inferred")
plot!(1:T, posterior_mean[2:end],
    ribbon=[sqrt.(posterior_var[2:end]), sqrt.(posterior_var[2:end])],
    color="purple", alpha=0.1, label="")
xlabel!("time (t)")
ylabel!("signal")
savefig(pwd()*"/viz/bootstrap.png")

if viz_particles

    # Show sample spreads
    scatter(1:T, observations, color="blue", label="observations")
    for t = 2:T
        scatter!((t-1).*ones(1,N), samples[t,:], markersize=weights[t,:]*50, markercolor="purple", label="")
    end
    savefig(pwd()*"/viz/bootstrap_particles.png")
end

"""RB-SIR filter"""

# Nominal distributions
α = 0.5
imp_pp = inv(process_noise)
imp_mp = inv(measurement_noise)
nominal_dynamics(u) = Normal(u, sqrt(α))
importance_dist(u,y) = Normal(inv(imp_pp + imp_mp)*(imp_pp*u + imp_mp*y), sqrt(inv(imp_pp + imp_mp)))

# Functions
A(x) = x
Q(x) = x^2
H(x) = x
R(x) = x^2

# Call filter
samples, weights, means, cvars = rao_blackwellized_particle_filter(observations,
                                                                   nominal_dynamics,
                                                                   importance_dist,
                                                                   A, Q, H, R,
                                                                   num_particles=N,
                                                                   resampling_threshold=0.5)

# Generate posterior distribution
posterior_mean = sum(weights .* samples .* means, dims=2)
posterior_cvar = sum(weights .* samples.^2 .* cvars, dims=2)

# Visualize particles
scatter(1:T, observations, color="blue", label="observations")
plot!(1:T, latent_states[2:end], color="red", label="latent states")
plot!(1:T, posterior_mean[2:end], color="purple", label="inferred")
plot!(1:T, posterior_mean[2:end],
      ribbon=[sqrt.(posterior_cvar[2:end]), sqrt.(posterior_cvar[2:end])],
      color="purple", alpha=0.1, label="")
xlabel!("time (t)")
ylabel!("signal")
savefig(pwd()*"/viz/rb-sir.png")

if viz_particles

    # Show sample spreads
    scatter(1:T, observations, color="blue", label="observations")
    for t = 2:T
        scatter!((t-1).*ones(1,N), samples[t,:], markersize=weights[t,:]*50, markercolor="purple", label="")
    end
    savefig(pwd()*"/viz/rb-sir_particles.png")
end

max_stars_repo_name: wmkouw/kalman-filters

max_issues_repo_name: wmkouw/kalman-filters

max_forks_repo_name: wmkouw/kalman-filters


text:
@symbol_func function B_rho(cur_reactor::AbstractReactor, cur_rho::AbstractSymbol, cur_theta::AbstractSymbol)
  cur_angle = cur_theta * ( cur_reactor.pi / 180.0 )

  cur_B = cur_reactor.B_0

  cur_B /= 1 + cur_reactor.epsilon * cur_rho * cos(cur_angle)

  cur_B
end


text:
"""
    ThinPlateSpline <: Kernel

TODO
"""
struct ThinPlateSpline <: Kernel
    # TODO
end


text:
#RSDeltaSigmaPort: Generate plots of the complex Z-plane
#-------------------------------------------------------------------------------

#==Waveform builders
===============================================================================#
function plot_unitcircle!(plot; color=:black, width=2, nseg::Int=360)
	Θ = range(0, stop=2π, length=nseg+1)
	cplxcircle = exp.((2π*j)*Θ)
	circle = DataF1(real(cplxcircle), imag(cplxcircle))
	push!(plot,
		cons(:wfrm, circle, line=set(style=:solid, color=color, width=width))
	)
	return plot
end

#==Pole-zero plots
===============================================================================#
#Create empty pole-zero plot
function plotPZ()
	plot = cons(:plot, title = "Pole-Zero Plot", legend=false,
		xyaxes=set(xscale=:lin, yscale=:lin, xmin=-1.1, xmax=1.1, ymin=-1.1, ymax=1.1),
		labels=set(xaxis="Real", yaxis="Imag"),
	)
	plot_unitcircle!(plot)
	return plot
end

"""`function plotPZ!(H,color=:black,markersize=5,list=false)`
Plot the poles and zeros of a transfer function.
If list=true, a list of the poles and zeros is superimposed on the plot.
"""
function plotPZ!(plot, H::ZPKData; color=:blue, markersize=2, list::Bool=false)
	isVector = isa(color, Vector)
	color_pole = isVector ? color[1] : color
	color_zero = isVector ? color[end] : color
	attrib_poles = cons(:a,
		glyph=set(shape=:x, color=color_pole, size=markersize), line=set(style=:none, width=2)
	)
	attrib_zeros = cons(:a,
		glyph=set(shape=:o, color=color_zero, size=markersize), line=set(style=:none, width=2)
	)

	z, p, k = _zpkdata(H)

	# Plot x and o for poles and zeros, respectively:
	push!(plot, cons(:wfrm, DataF1(real(p), imag(p)), attrib_poles))
	if !isempty(z) #Add zeros
		push!(plot, cons(:wfrm, DataF1(real(z), imag(z)), attrib_zeros))
	end

#=
	if list
		# List the poles and zeros
		p = cplxpair(p);
		y = 0.05*(ceil(length(p)/2)+1);
		str = 'Poles:               ';
		text( 0, y, str, 'Hor', 'Right', 'Ver', 'Mid'); y = y - 0.1;
		for i = 1:2:length(p);
			if abs(imag(p(i))) < 1e-6
				str = sprintf('%+.4f      ', real(p(i)) );
			else
				str = sprintf( '%+.4f+/-j%.4f  ', real(p(i)), abs(imag(p(i))) );
			end
			text( 0, y, str, 'Hor', 'Right', 'Ver', 'Mid'); y = y - 0.1;
		end
		if !isempty(z)
			z = z( !isnan(z) && !isinf(z) );
			z = cplxpair(z);
			y = 0.05*(ceil(length(z)/2)+1);
			str = '        Zeros:';
			text( 0, y, str, 'Hor', 'Left', 'Ver', 'Mid'); y = y - 0.1;
			for i = 1:2:length(z);
				if abs(imag(z(i))) < 1e-6
					str = sprintf('%+.4f      ', real(z(i)) );
				else
					str = sprintf( '  %+.4f+/-j%.4f', real(z(i)), abs(imag(z(i))) );
				end
				text( 0, y, str, 'Hor', 'Left', 'Ver', 'Mid');
				y = y - 0.1;
			end
		end
	end
=#
	return plot
end

plotPZ(H::ZPKData, args...; kwargs...) = plotPZ!(plotPZ(), H::ZPKData, args...; kwargs...)

#Last line