text:
export TinyVectorObservable, stddev

mutable struct TinyVectorObservable <: VectorObservable
    num :: Int64
    sum :: Vector{Float64}
    sum2 :: Vector{Float64}
end

TinyVectorObservable() = TinyVectorObservable(0, Float64[], Float64[])
zero(::Type{TinyVectorObservable}) = TinyVectorObservable()
zero(o::TinyVectorObservable) = TinyVectorObservable()
zeros(::Type{TinyVectorObservable}, dims...) = reshape([zero(TinyVectorObservable) for i in 1:prod(dims)],dims)

function reset!(obs :: TinyVectorObservable)
    obs.num = 0
    obs.sum = Float64[]
    obs.sum2 = Float64[]
    return obs
end

count(obs::TinyVectorObservable) = obs.num

function push!(obs :: TinyVectorObservable, value::Vector)
    if obs.num == 0
        obs.num = 1
        obs.sum = deepcopy(value)
        obs.sum2 = squared(value)
    else
        obs.num += 1
        for i in 1:length(obs.sum)
            obs.sum[i] += value[i]
            obs.sum2[i] += squared(value[i])
        end
    end
    return obs
end

function mean(obs::TinyVectorObservable)
    if obs.num > 0
        return obs.sum .* (1.0/obs.num)
    else
        return [NaN]
    end
end

function var(obs::TinyVectorObservable)
    if obs.num  > 1
        v = (obs.sum2 .- (obs.sum.^2)./obs.num)./(obs.num-1)
        return map!(maxzero, v)
    else
        return fill(NaN, length(obs.sum))
    end
end
stddev(obs::TinyVectorObservable) = sqrt.(var(obs))
stderror(obs::TinyVectorObservable) = sqrt.(var(obs)./count(obs))
function confidence_interval(obs::TinyVectorObservable, confidence_rate :: Real)
    if count(obs) == 0
        return [Inf]
    elseif count(obs) == 1
        return fill(Inf, length(obs.sum))
    end
    q = 0.5 + 0.5confidence_rate
    correction = quantile( TDist(obs.num - 1), q)
    serr = stderror(obs)
    return correction * serr
end

function confidence_interval(obs::TinyVectorObservable, confidence_rate_symbol::Symbol = :sigma1)
    n = parsesigma(confidence_rate_symbol)
    return confidence_interval(obs, erf(0.5n*sqrt(2.0)))
end

function merge!(obs::TinyVectorObservable, other::TinyVectorObservable)
    obs.num += other.num
    for i in 1:length(obs.sum)
        obs.sum[i] += other.sum[i]
        obs.sum2[i] += other.sum2[i]
    end
    return obs
end
merge(lhs::TinyVectorObservable, rhs::TinyVectorObservable) = merge!(deepcopy(lhs), rhs)

export TinyVectorObservableSet
const TinyVectorObservableSet = MCObservableSet{TinyVectorObservable}

function merge!(obs::TinyVectorObservableSet, other::TinyVectorObservableSet)
    obs_names = Set(keys(obs))
    union!(obs_names, Set(keys(other)))
    for name in obs_names
        if !haskey(obs, name)
            # in other only
            obs[name] = deepcopy(other[name])
        elseif haskey(other, name)
            # in both
            merge!(obs[name], other[name])

            # else
            # in obs only
            # NOTHING to do
        end
    end
    return obs
end
merge(lhs::TinyVectorObservableSet, rhs::TinyVectorObservableSet) = merge!(deepcopy(lhs), rhs)


text:
# Autogenerated wrapper script for GRASS_jll for x86_64-w64-mingw32-cxx11
export libgrass_arraystats, libgrass_bitmap, libgrass_btree, libgrass_btree2, libgrass_cairodriver, libgrass_calc, libgrass_ccmath, libgrass_cdhc, libgrass_cluster, libgrass_datetime, libgrass_dbmibase, libgrass_dbmiclient, libgrass_dbmidriver, libgrass_dbstubs, libgrass_dgl, libgrass_dig2, libgrass_display, libgrass_driver, libgrass_dspf, libgrass_g3d, libgrass_gis, libgrass_gmath, libgrass_gpde, libgrass_gproj, libgrass_htmldriver, libgrass_imagery, libgrass_interpdata, libgrass_interpfl, libgrass_iortho, libgrass_lidar, libgrass_linkm, libgrass_lrs, libgrass_manage, libgrass_neta, libgrass_pngdriver, libgrass_psdriver, libgrass_qtree, libgrass_raster, libgrass_rli, libgrass_rowio, libgrass_rtree, libgrass_segment, libgrass_shape, libgrass_sim, libgrass_sqlp, libgrass_stats, libgrass_symb, libgrass_temporal, libgrass_vector, libgrass_vedit

using GDAL_jll
using PROJ_jll
using Cairo_jll
using FreeType2_jll
using LAPACK_jll
using FFTW_jll
using Libiconv_jll
JLLWrappers.@generate_wrapper_header("GRASS")
JLLWrappers.@declare_library_product(libgrass_arraystats, "libgrass_arraystats.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_bitmap, "libgrass_bitmap.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_btree, "libgrass_btree.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_btree2, "libgrass_btree2.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_cairodriver, "libgrass_cairodriver.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_calc, "libgrass_calc.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_ccmath, "libgrass_ccmath.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_cdhc, "libgrass_cdhc.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_cluster, "libgrass_cluster.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_datetime, "libgrass_datetime.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_dbmibase, "libgrass_dbmibase.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_dbmiclient, "libgrass_dbmiclient.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_dbmidriver, "libgrass_dbmidriver.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_dbstubs, "libgrass_dbstubs.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_dgl, "libgrass_dgl.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_dig2, "libgrass_dig2.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_display, "libgrass_display.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_driver, "libgrass_driver.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_dspf, "libgrass_dspf.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_g3d, "libgrass_g3d.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_gis, "libgrass_gis.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_gmath, "libgrass_gmath.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_gpde, "libgrass_gpde.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_gproj, "libgrass_gproj.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_htmldriver, "libgrass_htmldriver.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_imagery, "libgrass_imagery.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_interpdata, "libgrass_interpdata.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_interpfl, "libgrass_interpfl.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_iortho, "libgrass_iortho.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_lidar, "libgrass_lidar.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_linkm, "libgrass_linkm.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_lrs, "libgrass_lrs.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_manage, "libgrass_manage.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_neta, "libgrass_neta.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_pngdriver, "libgrass_pngdriver.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_psdriver, "libgrass_psdriver.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_qtree, "libgrass_qtree.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_raster, "libgrass_raster.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_rli, "libgrass_rli.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_rowio, "libgrass_rowio.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_rtree, "libgrass_rtree.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_segment, "libgrass_segment.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_shape, "libgrass_shape.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_sim, "libgrass_sim.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_sqlp, "libgrass_sqlp.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_stats, "libgrass_stats.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_symb, "libgrass_symb.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_temporal, "libgrass_temporal.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_vector, "libgrass_vector.7.8.dll")
JLLWrappers.@declare_library_product(libgrass_vedit, "libgrass_vedit.7.8.dll")
function __init__()
    JLLWrappers.@generate_init_header(GDAL_jll, PROJ_jll, Cairo_jll, FreeType2_jll, LAPACK_jll, FFTW_jll, Libiconv_jll)
    JLLWrappers.@init_library_product(
        libgrass_arraystats,
        "bin\\libgrass_arraystats.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_bitmap,
        "bin\\libgrass_bitmap.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_btree,
        "bin\\libgrass_btree.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_btree2,
        "bin\\libgrass_btree2.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_cairodriver,
        "bin\\libgrass_cairodriver.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_calc,
        "bin\\libgrass_calc.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_ccmath,
        "bin\\libgrass_ccmath.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_cdhc,
        "bin\\libgrass_cdhc.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_cluster,
        "bin\\libgrass_cluster.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_datetime,
        "bin\\libgrass_datetime.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_dbmibase,
        "bin\\libgrass_dbmibase.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_dbmiclient,
        "bin\\libgrass_dbmiclient.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_dbmidriver,
        "bin\\libgrass_dbmidriver.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_dbstubs,
        "bin\\libgrass_dbstubs.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_dgl,
        "bin\\libgrass_dgl.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_dig2,
        "bin\\libgrass_dig2.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_display,
        "bin\\libgrass_display.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_driver,
        "bin\\libgrass_driver.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_dspf,
        "bin\\libgrass_dspf.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_g3d,
        "bin\\libgrass_g3d.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_gis,
        "bin\\libgrass_gis.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_gmath,
        "bin\\libgrass_gmath.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_gpde,
        "bin\\libgrass_gpde.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_gproj,
        "bin\\libgrass_gproj.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_htmldriver,
        "bin\\libgrass_htmldriver.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_imagery,
        "bin\\libgrass_imagery.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_interpdata,
        "bin\\libgrass_interpdata.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_interpfl,
        "bin\\libgrass_interpfl.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_iortho,
        "bin\\libgrass_iortho.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_lidar,
        "bin\\libgrass_lidar.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_linkm,
        "bin\\libgrass_linkm.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_lrs,
        "bin\\libgrass_lrs.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_manage,
        "bin\\libgrass_manage.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_neta,
        "bin\\libgrass_neta.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_pngdriver,
        "bin\\libgrass_pngdriver.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_psdriver,
        "bin\\libgrass_psdriver.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_qtree,
        "bin\\libgrass_qtree.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_raster,
        "bin\\libgrass_raster.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_rli,
        "bin\\libgrass_rli.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_rowio,
        "bin\\libgrass_rowio.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_rtree,
        "bin\\libgrass_rtree.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_segment,
        "bin\\libgrass_segment.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_shape,
        "bin\\libgrass_shape.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_sim,
        "bin\\libgrass_sim.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_sqlp,
        "bin\\libgrass_sqlp.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_stats,
        "bin\\libgrass_stats.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_symb,
        "bin\\libgrass_symb.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_temporal,
        "bin\\libgrass_temporal.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_vector,
        "bin\\libgrass_vector.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgrass_vedit,
        "bin\\libgrass_vedit.7.8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
using Test
using AutomotiveSimulator
using DataFrames
using Distributions

@testset "Vec" begin
    include("vec-tests/vec_runtests.jl")
end

@testset "AutomotiveSimulator" begin
    include("test_roadways.jl")
    include("test_agent_definitions.jl")
    include("test_scenes.jl")
    include("test_states.jl")
    include("test_collision_checkers.jl")
    include("test_actions.jl")
    include("test_features.jl")
    include("test_behaviors.jl")
    include("test_simulation.jl")
end

# Restore after AutomotiveVisualization v0.9 is tagged
@testset "doc tutorials" begin
    @testset "basics" begin
        include("../docs/lit/tutorials/straight_roadway.jl")
    end
    @testset "cameras" begin
        include("../docs/lit/tutorials/stadium.jl")
    end
    @testset "overlays" begin
        include("../docs/lit/tutorials/intersection.jl")
    end
    @testset "basics" begin
        include("../docs/lit/tutorials/crosswalk.jl")
    end
    @testset "basics" begin
        include("../docs/lit/tutorials/sidewalk.jl")
    end
end


text:
replacein(x, oldnode, newnode) = x
replacein(x::Tuple, oldnode, newnode) = tuple([replacein(m, oldnode, newnode) for m in x]...)

# macro replacein(T)
# 	quote
# 		function replacein(x::$(T), oldnode, newnode)
# 			x == oldnode && return(newnode)
# 			$(T)(map(f -> replacein(getproperty(x, f), oldnode, newnode), fieldnames(T))...)
# 		end
# 	end
# end

for T in [ArrayNode, ArrayModel, BagModel, BagNode, TreeNode, ProductModel]
	@eval begin
		function replacein(x::$T, oldnode, newnode)
			x == oldnode && return(newnode)
			$T(map(f -> replacein(getproperty(x, f), oldnode, newnode), fieldnames($T))...)
		end
	end
end

# replacein(x::ArrayModel, oldnode, newnode) = x == oldnode ? newnode : ArrayModel(replacein(x.m, oldnode, newnode))
# replacein(x::ArrayNode, oldnode, newnode) = x == oldnode ? newnode : ArrayNode(replacein(x.data, oldnode, newnode), x.metadata)

# replacein(x::BagModel, oldnode, newnode) = x == oldnode ? newnode : BagModel(replacein(x.im, oldnode, newnode), replacein(x.a, oldnode, newnode), replacein(x.bm, oldnode, newnode))
# replacein(x::BagNode, oldnode, newnode) = x == oldnode ? newnode : BagNode(replacein(x.data, oldnode, newnode), x.bags, x.metadata)

# replacein(x::TreeNode, oldnode, newnode) = x == oldnode ? newnode : TreeNode(tuple([replacein(m, oldnode, newnode) for m in x.data]...), x.metadata)
# replacein(x::ProductModel, oldnode, newnode) = x == oldnode ? newnode : ProductModel(tuple([replacein(m, oldnode, newnode) for m in x.ms]...), replacein(x.m, oldnode, newnode))


text:
using Base.Test
using ExclusionProcesses

@testset "epstate" begin

n_sites = 10

# Basic constructor
random_spins = bitrand(n_sites)
@test SingleSpinState(random_spins).spins == random_spins

# Randomly initialized states have correct number of particles
n_particles = rand(1:n_sites)
@test sum(SingleSpinState(n_sites,n_particles).spins) == n_particles

# Every configuration in a multi spin encoding has correct number of particles
configs = spin_configs(MultiSpinState(n_sites, n_particles))
@test all([sum(configs[i,:]) == n_particles for i=1:64])

end # testset

max_stars_repo_name: AustenLamacraft/ExclusionProcesses

max_issues_repo_name: AustenLamacraft/ExclusionProcesses

max_forks_repo_name: AustenLamacraft/ExclusionProcesses


text:
#==============================================================================#
# Macie.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/macie-2017-12-19.normal.json
#==============================================================================#

__precompile__()

module Macie

using AWSCore

"""
    using AWSSDK.Macie.associate_member_account
    associate_member_account([::AWSConfig], arguments::Dict)
    associate_member_account([::AWSConfig]; memberAccountId=)

    using AWSCore.Services.macie
    macie([::AWSConfig], "AssociateMemberAccount", arguments::Dict)
    macie([::AWSConfig], "AssociateMemberAccount", memberAccountId=)

# AssociateMemberAccount Operation

Associates a specified AWS account with Amazon Macie as a member account.

# Arguments

## `memberAccountId = ::String` -- *Required*
The ID of the AWS account that you want to associate with Amazon Macie as a member account.

# Exceptions

`InvalidInputException`, `LimitExceededException` or `InternalException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/macie-2017-12-19/AssociateMemberAccount)
"""
@inline associate_member_account(aws::AWSConfig=default_aws_config(); args...) = associate_member_account(aws, args)

@inline associate_member_account(aws::AWSConfig, args) = AWSCore.Services.macie(aws, "AssociateMemberAccount", args)

@inline associate_member_account(args) = associate_member_account(default_aws_config(), args)

"""
    using AWSSDK.Macie.associate_s3resources
    associate_s3resources([::AWSConfig], arguments::Dict)
    associate_s3resources([::AWSConfig]; s3Resources=, <keyword arguments>)

    using AWSCore.Services.macie
    macie([::AWSConfig], "AssociateS3Resources", arguments::Dict)
    macie([::AWSConfig], "AssociateS3Resources", s3Resources=, <keyword arguments>)

# AssociateS3Resources Operation

Associates specified S3 resources with Amazon Macie for monitoring and data classification. If memberAccountId isn't specified, the action associates specified S3 resources with Macie for the current master account. If memberAccountId is specified, the action associates specified S3 resources with Macie for the specified member account.

# Arguments

## `memberAccountId = ::String`
The ID of the Amazon Macie member account whose resources you want to associate with Macie.

## `s3Resources = [[ ... ], ...]` -- *Required*
The S3 resources that you want to associate with Amazon Macie for monitoring and data classification.
```
 s3Resources = [[
        "bucketName" => <required> ::String,
        "prefix" =>  ::String,
        "classificationType" => <required> [
            "oneTime" => <required> "FULL" or "NONE",
            "continuous" => <required> "FULL"
        ]
    ], ...]
```

# Returns

`AssociateS3ResourcesResult`

# Exceptions

`InvalidInputException`, `AccessDeniedException`, `LimitExceededException` or `InternalException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/macie-2017-12-19/AssociateS3Resources)
"""
@inline associate_s3resources(aws::AWSConfig=default_aws_config(); args...) = associate_s3resources(aws, args)

@inline associate_s3resources(aws::AWSConfig, args) = AWSCore.Services.macie(aws, "AssociateS3Resources", args)

@inline associate_s3resources(args) = associate_s3resources(default_aws_config(), args)

"""
    using AWSSDK.Macie.disassociate_member_account
    disassociate_member_account([::AWSConfig], arguments::Dict)
    disassociate_member_account([::AWSConfig]; memberAccountId=)

    using AWSCore.Services.macie
    macie([::AWSConfig], "DisassociateMemberAccount", arguments::Dict)
    macie([::AWSConfig], "DisassociateMemberAccount", memberAccountId=)

# DisassociateMemberAccount Operation

Removes the specified member account from Amazon Macie.

# Arguments

## `memberAccountId = ::String` -- *Required*
The ID of the member account that you want to remove from Amazon Macie.

# Exceptions

`InvalidInputException` or `InternalException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/macie-2017-12-19/DisassociateMemberAccount)
"""
@inline disassociate_member_account(aws::AWSConfig=default_aws_config(); args...) = disassociate_member_account(aws, args)

@inline disassociate_member_account(aws::AWSConfig, args) = AWSCore.Services.macie(aws, "DisassociateMemberAccount", args)

@inline disassociate_member_account(args) = disassociate_member_account(default_aws_config(), args)

"""
    using AWSSDK.Macie.disassociate_s3resources
    disassociate_s3resources([::AWSConfig], arguments::Dict)
    disassociate_s3resources([::AWSConfig]; associatedS3Resources=, <keyword arguments>)

    using AWSCore.Services.macie
    macie([::AWSConfig], "DisassociateS3Resources", arguments::Dict)
    macie([::AWSConfig], "DisassociateS3Resources", associatedS3Resources=, <keyword arguments>)

# DisassociateS3Resources Operation

Removes specified S3 resources from being monitored by Amazon Macie. If memberAccountId isn't specified, the action removes specified S3 resources from Macie for the current master account. If memberAccountId is specified, the action removes specified S3 resources from Macie for the specified member account.

# Arguments

## `memberAccountId = ::String`
The ID of the Amazon Macie member account whose resources you want to remove from being monitored by Amazon Macie.

## `associatedS3Resources = [[ ... ], ...]` -- *Required*
The S3 resources (buckets or prefixes) that you want to remove from being monitored and classified by Amazon Macie.
```
 associatedS3Resources = [[
        "bucketName" => <required> ::String,
        "prefix" =>  ::String
    ], ...]
```

# Returns

`DisassociateS3ResourcesResult`

# Exceptions

`InvalidInputException`, `AccessDeniedException` or `InternalException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/macie-2017-12-19/DisassociateS3Resources)
"""
@inline disassociate_s3resources(aws::AWSConfig=default_aws_config(); args...) = disassociate_s3resources(aws, args)

@inline disassociate_s3resources(aws::AWSConfig, args) = AWSCore.Services.macie(aws, "DisassociateS3Resources", args)

@inline disassociate_s3resources(args) = disassociate_s3resources(default_aws_config(), args)

"""
    using AWSSDK.Macie.list_member_accounts
    list_member_accounts([::AWSConfig], arguments::Dict)
    list_member_accounts([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.macie
    macie([::AWSConfig], "ListMemberAccounts", arguments::Dict)
    macie([::AWSConfig], "ListMemberAccounts", <keyword arguments>)

# ListMemberAccounts Operation

Lists all Amazon Macie member accounts for the current Amazon Macie master account.

# Arguments

## `nextToken = ::String`
Use this parameter when paginating results. Set the value of this parameter to null on your first call to the ListMemberAccounts action. Subsequent calls to the action fill nextToken in the request with the value of nextToken from the previous response to continue listing data.

## `maxResults = ::Int`
Use this parameter to indicate the maximum number of items that you want in the response. The default value is 250.

# Returns

`ListMemberAccountsResult`

# Exceptions

`InternalException` or `InvalidInputException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/macie-2017-12-19/ListMemberAccounts)
"""
@inline list_member_accounts(aws::AWSConfig=default_aws_config(); args...) = list_member_accounts(aws, args)

@inline list_member_accounts(aws::AWSConfig, args) = AWSCore.Services.macie(aws, "ListMemberAccounts", args)

@inline list_member_accounts(args) = list_member_accounts(default_aws_config(), args)

"""
    using AWSSDK.Macie.list_s3resources
    list_s3resources([::AWSConfig], arguments::Dict)
    list_s3resources([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.macie
    macie([::AWSConfig], "ListS3Resources", arguments::Dict)
    macie([::AWSConfig], "ListS3Resources", <keyword arguments>)

# ListS3Resources Operation

Lists all the S3 resources associated with Amazon Macie. If memberAccountId isn't specified, the action lists the S3 resources associated with Amazon Macie for the current master account. If memberAccountId is specified, the action lists the S3 resources associated with Amazon Macie for the specified member account.

# Arguments

## `memberAccountId = ::String`
The Amazon Macie member account ID whose associated S3 resources you want to list.

## `nextToken = ::String`
Use this parameter when paginating results. Set its value to null on your first call to the ListS3Resources action. Subsequent calls to the action fill nextToken in the request with the value of nextToken from the previous response to continue listing data.

## `maxResults = ::Int`
Use this parameter to indicate the maximum number of items that you want in the response. The default value is 250.

# Returns

`ListS3ResourcesResult`

# Exceptions

`InvalidInputException`, `AccessDeniedException` or `InternalException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/macie-2017-12-19/ListS3Resources)
"""
@inline list_s3resources(aws::AWSConfig=default_aws_config(); args...) = list_s3resources(aws, args)

@inline list_s3resources(aws::AWSConfig, args) = AWSCore.Services.macie(aws, "ListS3Resources", args)

@inline list_s3resources(args) = list_s3resources(default_aws_config(), args)

"""
    using AWSSDK.Macie.update_s3resources
    update_s3resources([::AWSConfig], arguments::Dict)
    update_s3resources([::AWSConfig]; s3ResourcesUpdate=, <keyword arguments>)

    using AWSCore.Services.macie
    macie([::AWSConfig], "UpdateS3Resources", arguments::Dict)
    macie([::AWSConfig], "UpdateS3Resources", s3ResourcesUpdate=, <keyword arguments>)

# UpdateS3Resources Operation

Updates the classification types for the specified S3 resources. If memberAccountId isn't specified, the action updates the classification types of the S3 resources associated with Amazon Macie for the current master account. If memberAccountId is specified, the action updates the classification types of the S3 resources associated with Amazon Macie for the specified member account.

# Arguments

## `memberAccountId = ::String`
The AWS ID of the Amazon Macie member account whose S3 resources' classification types you want to update.

## `s3ResourcesUpdate = [[ ... ], ...]` -- *Required*
The S3 resources whose classification types you want to update.
```
 s3ResourcesUpdate = [[
        "bucketName" => <required> ::String,
        "prefix" =>  ::String,
        "classificationTypeUpdate" => <required> [
            "oneTime" =>  "FULL" or "NONE",
            "continuous" =>  "FULL"
        ]
    ], ...]
```

# Returns

`UpdateS3ResourcesResult`

# Exceptions

`InvalidInputException`, `AccessDeniedException` or `InternalException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/macie-2017-12-19/UpdateS3Resources)
"""
@inline update_s3resources(aws::AWSConfig=default_aws_config(); args...) = update_s3resources(aws, args)

@inline update_s3resources(aws::AWSConfig, args) = AWSCore.Services.macie(aws, "UpdateS3Resources", args)

@inline update_s3resources(args) = update_s3resources(default_aws_config(), args)

end # module Macie

#==============================================================================#
# End of file
#==============================================================================#


text:
using ESDC
using Test


text:
mutable struct DenseSolver <: KKTSolver{Scaling}
	k0::Vector{Float64}
	k1::Vector{Float64}
	k2::Vector{Float64}
	m0::Vector{Float64}
	n0::Vector{Float64}
	n1::Vector{Float64}
	GWiWi::Matrix{Float64}
	GWiWiG::Matrix{Float64}
	Li::Matrix{Float64}
	AtLi::Matrix{Float64}
	AtLiA::Matrix{Float64}
	GWiWiGherm::Hermitian{Float64, Array{Float64, 2}}
	AtLiAherm::Hermitian{Float64, Array{Float64, 2}}
	GWiWiGfact::Union{Nothing, Cholesky{Float64, Array{Float64, 2}}}
	AtLiAfact::Union{Nothing, Cholesky{Float64, Array{Float64, 2}}}
	eyetgt::Matrix{Float64}
	AA::Matrix{Float64}
	function DenseSolver(pr::Problem)
		k0, k1, k2 = zeros(pr.k), zeros(pr.k), zeros(pr.k)
		m0 = zeros(pr.m)
		n0, n1 = zeros(pr.n), zeros(pr.n)
		GWiWi = zeros(pr.n, pr.k)
		GWiWiG = zeros(pr.n, pr.n)
		GWiWiGherm = Hermitian(GWiWiG)
		Li = zeros(pr.n, pr.n)
		AtLi = zeros(pr.m, pr.n)
		AtLiA = zeros(pr.m, pr.m)
		AtLiAherm = Hermitian(AtLiA)
		eyetgt = Matrix{Float64}(I,pr.n, pr.n)
		AA = zeros(pr.n, pr.n)
		mul!(AA, pr.A', pr.A)
		return new(k0, k1, k2, m0, n0, n1,
			GWiWi, GWiWiG, Li, AtLi, AtLiA,
			GWiWiGherm, AtLiAherm,
			nothing, nothing,
			eyetgt, AA)
	end
end

function setup_iter(ss::DenseSolver, pr::Problem{C,n,m,k,sing}, s::State, scaling::Scaling) where {C, n, m, k, sing}
	mul!(ss.GWiWi, pr.G', scaling.iWiW) #GWiWi = G'*W^-1 * W^-T
	mul!(ss.GWiWiG, ss.GWiWi, pr.G) #ss.GWiWiG = G'*W^-1*W^-T*G
	if sing # ss.GWiWiG will be not pos def if it wasn't before
		ss.GWiWiG .+= ss.AA # we know that aa's been populated by now
	end
	ss.GWiWiGfact = cholesky!(ss.GWiWiGherm)
	ldiv!(ss.Li, ss.GWiWiGfact, et) #iR = L^-T L^-1
	mul!(ss.ALi, pr.A, ss.Li)
	mul!(ss.AtLiA, ss.ALi, At)
	ss.AtLiAfact = cholesky!(ss.AtLiAherm)
end

function solve_kkt(ss::DenseSolver, pr::Problem{C,n,m,k,sing}, s::State, scaling::Scaling,
				   dx::Vector{Float64}, dy::Vector{Float64}, dz::Vector{Float64}, ds::Vector{Float64},
				   cx::Vector{Float64}, cy::Vector{Float64}, cz::Vector{Float64}, cs::Vector{Float64}) where {C, n, m, k, sing}
	W,iW,l = scaling.W,scaling.iW,scaling.l
	et = ss.eyetgt
	aa = ss.AA

	iprod!(pr.cones, ss.k0, l, ds)
	scale!(pr.cones, scaling, ss.k0, ss.k1)
	bx = dx
	by = dy
	ss.k2 .= dz .- ss.k1
	mul!(ss.n0, ss.GWiWi, ss.k2)
	ss.n0 .+= bx
	At = pr.A'
	if ss.issng
		ss.n0 .+= At*by
	end

	mul!(ss.m0, ss.AtLi, ss.n0)
	ss.m0 .-= by
	ldiv!(cy, ss.AtLiAfact, ss.m0)
	if ss.issng
		ss.m0 .= by .- cy
	else
		ss.m0 .= .-cy
	end
	mul!(ss.n1, At, ss.m0)
	ss.n0 .+= ss.n1
	mul!(cx, ss.Li, ss.n0)
	mul!(ss.k1, pr.G, cx)
	ss.k1 .-= ss.k2
	mul!(cz, scaling.iWiW, ss.k1)
	scale!(pr.cones, scaling, cz, ss.k1)
	ss.k0 .-= ss.k1
	scale!(pr.cones, scaling, ss.k0, cs)
end


text:
using XGPaint
using Documenter
using Healpix
using Cosmology

makedocs(;
    modules=[XGPaint],
    authors="Zack Li",
    repo="https://github.com/xzackli/XGPaint.jl/blob/{commit}{path}#L{line}",
    sitename="XGPaint.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://xzackli.github.io/XGPaint.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "CIB" => "cib.md",
        "Index" => "api.md",
    ],
)

deploydocs(;
    repo="github.com/xzackli/XGPaint.jl",
    devbranch = "main"
)


text:
using TaylorSeries: set_variables
using LazySets.Approximations: box_approximation

function post(𝒜::TMJets,
              𝑃::InitialValueProblem{<:Union{BBCS, CBBCS, CBBCCS}, <:LazySet},
              𝑂_global::Options)

    # ==================================
    # Initialization
    # ==================================

    𝑂 = merge(𝒜.options.defaults, 𝑂_global, 𝒜.options.specified)

    # system of ODEs
    f! = 𝑃.s.f
    n = statedim(𝑃)

    # initial time and time horizon
    t0 = 0.0
    T = 𝑂[:T]

    # maximum allowed number of steps
    max_steps = 𝑂[:max_steps]

    # unrap algorithm-specific options
    abs_tol, orderQ, orderT = 𝑂[:abs_tol], 𝑂[:orderQ], 𝑂[:orderT]

    # initial sets
    box_x0 = box_approximation(𝑃.x0)
    q0 = center(box_x0)
    δq0 = IntervalBox(low(box_x0)-q0, high(box_x0)-q0)

    # fix the working variables and maximum order in the global
    # parameters struct (_params_TaylorN_)
    set_variables("x", numvars=length(q0), order=2*orderQ)

    # define the property
    if 𝑂[:mode] == "check"
        property = 𝑂[:property]
    elseif 𝑂[:mode] == "reach"
        property = (t, x) -> true
    end

    # =====================
    # Flowpipe computation
    # =====================

    info("Reachable States Computation...")
    @timing begin
        tTM, xTM = validated_integ(f!, q0, δq0, t0, T, orderQ, orderT, abs_tol,
                                   maxsteps=max_steps, check_property=property)
    end

    # convert to hyperrectangle and wrap around the reach solution
    N = length(xTM)
    Rsets = Vector{ReachSet{Hyperrectangle{Float64}}}(undef, N-1)
    @inbounds for i in 1:N-1
        Hi = convert(Hyperrectangle, xTM[i])
        t0 = tTM[i]; t1 = tTM[i+1]
        Rsets[i] = ReachSet(Hi, t0, t1)
    end

    Rsol = ReachSolution(Rsets, 𝑂)

    # ===========
    # Projection
    # ===========

    if 𝑂[:project_reachset]
        info("Projection...")
        Rsol = @timing project(Rsol)
    end

    return Rsol
end


text:
__precompile__()

module PlanerTransformer
include("string_literals.jl")
include("types.jl")
include("TransformerAnalysis.jl")
include("functionalize_fieldnames.jl")
include("methods.jl")

# create ferrite and core dictionaries
include("DataEntry.jl")
include("FerriteData.jl")
include("CoreData.jl")

end # module


text:
function encode(tamanho, elemento) # Função que faz o encode dos aminoácidos e das classes
    a = zeros(Int64,tamanho)
    a[elemento] = 1
    return(a)
end

siglas = [ 'A', 'R', 'N', 'D', 'E', 'C', 'G', 'Q', 'H', 'I', 'L', 'K', 'M', 'F', 'P', 'S', 'Y', 'T', 'W', 'V', 'X' ]
codificadordeaminoacidos = Dict()

for (indice, elemento) in enumerate(siglas)
    codificadordeaminoacidos[elemento] = encode(length(siglas), indice)
end

const codificadordeSS8 = Dict(
    'H' => 'H',
    'G' => 'H',
    'I' => 'H',
    'P' => 'H',
    'E' => 'E',
    'B' => 'E',
    'C' => 'C',
    'T' => 'C',
    'S' => 'C',
    ' ' => 'C'
)

const codificadordeSS3 = Dict(
    'H' => [1,0,0],
    'E' => [0,1,0],
    'C' => [0,0,1]
)

function build_matrix_X(sequence::Vector{Char}, windowsize::Int)

    global codificadordeaminoacidos

    skip = Int64((windowsize-1)/2)
    X = Matrix{Int64}(undef, 21*windowsize, length(sequence)-2*skip) # que tenha tantas colunas quanto aminoácidos na sequência em questão, e tantas linhas quanto possíveis SS

    for i in 1:(length(sequence)-windowsize+1)

        fragment = sequence[i:i+windowsize-1]
        column = Vector{Int64}(undef, 21*windowsize)

        for aminoacid in 1:length(fragment)
            column[21*(aminoacid-1)+1:21*aminoacid] = codificadordeaminoacidos[fragment[aminoacid]]
        end

        X[:,i] = column

    end

    return(X)

end

function build_matrix_Y(structure::Vector{Char}, windowsize::Int)

    global codificadordeSS8
    global codificadordeSS3

    skip = Int64((windowsize-1)/2)
    Y = Matrix{Int64}(undef, 3, length(structure)-2*skip) # que tenha tantas colunas quanto aminoácidos na sequência em questão, e tantas linhas quanto possíveis SS

    structure_ss8 = [ codificadordeSS8[el] for el in structure ]

    for i in 1:(size(Y, 2)) # build matrix Y
        Y[:,i] = codificadordeSS3[structure_ss8[i+skip]]
    end  # build matrix Y

    return(Y)
end

using ProgressMeter

function prepare_data(listofsequences::Vector{Vector{Char}}, listofstructures::Vector{Vector{Char}}, windowsize::Int, porcentagem_split=0.6)

    skip = Int64((windowsize-1)/2)

    numero_de_proteinas = length(listofsequences)
    indice_split = round(Int,numero_de_proteinas*porcentagem_split)
    sequencias_de_treino = listofsequences[1:indice_split]
    sequencias_de_teste = listofsequences[(indice_split+1):end]
    estruturas_de_treino = listofstructures[1:indice_split]
    estruturas_de_teste = listofstructures[(indice_split+1):end]

    X_train = Matrix{Int64}(undef, 21*windowsize, 0)
    X_test = copy(X_train)
    Y_train = Matrix{Int64}(undef, 3, 0)
    Y_test = copy(Y_train)

    @showprogress for (sequence, structure) in zip(sequencias_de_treino, estruturas_de_treino) # for sequence, structure in (train)
        if length(sequence) >= windowsize
            X_train = hcat(X_train, build_matrix_X(sequence, windowsize))
            Y_train = hcat(Y_train, build_matrix_Y(structure, windowsize))
        else
            continue
        end
    end # for sequence,structure in train
    println("Done building the Train Datasets")

    @showprogress for (sequence, structure) in zip(sequencias_de_teste, estruturas_de_teste) # for sequence, structure in (train)
        if length(sequence) >= windowsize
            X_test = hcat(X_test, build_matrix_X(sequence, windowsize))
            Y_test = hcat(Y_test, build_matrix_Y(structure, windowsize))
        else
            continue
        end
    end # for sequence,structure in train
    println("Done building the Test Datasets")

    return X_train, Y_train, X_test, Y_test
end


text:
"""
    normalize(x,m,s)

Substract scalar `m` from array `x` and rescale by scalar `s`.
"""
function normalize(x,m,s)
    return (x.-m)./s
end

"""
    normalize_imgs(imgs, labels; is_train=true, m=0, s=1)

Take a onehot encoding of the labels.
If `is_train=true`, compute the mean and std of the images and normalize them.
Return the normalized images and the mean and std computed.
Otherwise, normalize the imgs with `m` and `s`.
Should be used with `is_train=true` on the train set and with `is_train=false`
otherwise with the parameters `m` and `s` computed on the training set.
"""
function normalize_imgs(imgs, labels; is_train=true, m=0, s=1)
    X = Array{Float32}(undef, size(imgs[1])..., 1, length(imgs))
    for i in 1:length(imgs)
        X[:, :, :, i] = Float32.(imgs[i])
    end
    Y = onehotbatch(labels, 0:9)
    if is_train
        m = mean(X)
        s = std(X)
        return (normalize(X,m,s), Y), m, s
    else
        return (normalize(X,m,s), Y)
    end
end

function get_processed_data_MNIST()
    # Load labels and images from Flux.Data.MNIST
    train_labels = MNIST.labels()
    train_imgs = MNIST.images()
    train_set, m, s = normalize_imgs(train_imgs,train_labels,is_train=true)

    test_imgs = MNIST.images(:test)
    test_labels = MNIST.labels(:test)
    test_set = normalize_imgs(test_imgs,test_labels,is_train=false,m=m,s=s)
    return train_set, test_set
end

abstract type AbstractDatabunch end

"""
    Databunch(train_loader::Flux.Data.DataLoader,val_loader::Flux.Data.DataLoader)

Two dataloaders, one for the train set and one for the validation set.
"""
struct Databunch <: AbstractDatabunch
    train_loader::Flux.Data.DataLoader
    val_loader::Flux.Data.DataLoader
end


text:
using Test
using BootlegCassette: BootlegCassette, @context, prehook, overdub, posthook, recurse
const Cassette = BootlegCassette

Cassette.@context TraceCtx

mutable struct Trace
    current::Vector{Any}
    stack::Vector{Any}
    Trace() = new(Any[], Any[])
end

function enter!(t::Trace, args...)
    pair = args => Any[]
    push!(t.current, pair)
    push!(t.stack, t.current)
    t.current = pair.second
    return nothing
end

function exit!(t::Trace)
    t.current = pop!(t.stack)
    return nothing
end

Cassette.prehook(ctx::TraceCtx, args...) = enter!(ctx.metadata, args...)
Cassette.posthook(ctx::TraceCtx, args...) = exit!(ctx.metadata)

trace = Trace()
x, y, z = rand(3)
f(x, y, z) = x*y + y*z
Cassette.overdub(TraceCtx(metadata = trace), () -> f(x, y, z))

@test trace.current == Any[
    (f,x,y,z) => Any[
        (*,x,y) => Any[(Base.mul_float,x,y)=>Any[]]
        (*,y,z) => Any[(Base.mul_float,y,z)=>Any[]]
        (+,x*y,y*z) => Any[(Base.add_float,x*y,y*z)=>Any[]]
    ]
]

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

Cassette.@context SinToCosCtx

# Override the default recursive `overdub` implementation for `sin(x)`.
# Note that there's no tricks here; this is just a normal Julia method
# overload using the normal multiple dispatch semantics.
Cassette.overdub(::SinToCosCtx, ::typeof(sin), x) = cos(x)

x = rand(10)
y = Cassette.overdub(SinToCosCtx(), sum, i -> cos(i) + sin(i), x)
@test y == sum(i -> 2 * cos(i), x)

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

fib(x) = x < 3 ? 1 : fib(x - 2) + fib(x - 1)
fibtest(n) = fib(2 * n) + n

@context MemoizeCtx

function Cassette.overdub(ctx::MemoizeCtx, ::typeof(fib), x)
    result = get(ctx.metadata, x, 0)
    if result === 0
        result = recurse(ctx, fib, x)
        ctx.metadata[x] = result
    end
    return result
end

ctx = MemoizeCtx(metadata=Dict{Int, Int}())

@test overdub(ctx, fibtest, 20) == fibtest(20)


text:
import YAML
import Dolang
path = Pkg.dir("Dolo","src","linter.jl")
include(path)

filename = Pkg.dir("Dolo","linter/models/rbc_iid.yaml")
filename

#DL = DoloLinter
d = yaml_node_from_file(filename)
errswarrs = lint(filename)

model_equation_types = keys(d[:equations])
model_equations = []
model_equations_vec = []
for key in keys(d[:equations])
    n = length(d[:equations][key].value)
    for i=1:n
        # would be cooler to do directly for sym in syms[vg]
        eq = d[:equations][key].value[i]
        push!(model_equations_vec, (key, eq.value, eq))
        push!(model_equations, eq.value)
    end
end

model_symbol_types = keys(d[:symbols])
model_symbols = []
model_symbols_vec = []
for key in keys(d[:symbols])
    n = length(d[:symbols][key].value)
    for i=1:n
        # would be cooler to do directly for sym in syms[vg]
        sym = d[:symbols][key].value[i]
        push!(model_symbols_vec, (key, sym.value, sym))
        push!(model_symbols, sym.value)
    end
end

for eq in model_equations_vec
  if eq[1] == "arbitrage"

    equ , comp = split(eq[2], " | ")
    for sym in keys(Dolang.IncidenceTable(parse(equ)).by_var)

    end
  end
end

import Dolo: RECIPES

errors = LinterWarning[]
warnings = LinterWarning[]

eqg = "arbitrage"
spec = RECIPES[:dtcc][:specs][:arbitrage][:eqs]
Symbol_set = Dict(eqg => Set(([el[2] for el in spec if el[1]==eqg])) for eqg in Set(el[1] for el in spec))

for eq_ in d[:equations][eqg].value
      eqs = eq_.value
      equ , comp = split(eqs, " | ")
      inc = Dolang.IncidenceTable(parse(equ)).by_var
      for sym in keys(inc)
        sym = string(sym)
        if !(sym in model_symbols)  && !(string(sym) in keys(d[:definitions]))
          errvalue = sym
          errtype = "Unkown Symbol"
          loc = get_loc(eq_)
          msg = string(sym, "  in the '", eqg , "' equation is not defined in symbols or definitions sections")
          println(sym)
          push!(errors, LinterWarning(errvalue, errtype, msg, loc, filename))

        elseif  (string(sym) in model_symbols)
          floc = findfirst(isequal(sym), model_symbols)
          declared_sym_type = model_symbols_vec[floc][1]
          println(declared_sym_type)

          if !(inc[sym] in Symbol_set[declared_sym_type])
            errvalue = sym
          ##  errtype = "Timing convention error"
          ##  loc = get_loc(eq_)
          ##  msg = string(sym, "  in the '", eqg , "' equation ")
          ##  println(sym)
          ##  push!(errors, LinterWarning(errvalue, errtype, msg, loc, filename))
          end
        end

      end

end

function format_human(exc::Union{LinterException, LinterWarning})
    # Output should be a nice one line human-readable strings
    # colors, spacing, etc...
end

function format_human(errors::Vector{LinterException}, warnings::Vector{LinterWarning})
    # ouptputs all errors then all warnings, line-by-line
    # example from python's `dolo-lint --format=human examples/models/rbc_iid.yaml`
    # error: 12, 30: Symbol 'beta' already declared as 'parameters'. (pos (12, 16))
end

function format_dict(exc::Union{LinterException, LinterWarning})
    # Output should be a dictionary like:
    # {"type": "error", "source": "examples/models/rbc_iid.yaml", "range": [[11, 29], [11, 33]], "text": "Symbol 'beta' already declared as 'parameters'."}
end

function format_json(errors::Vector{LinterException}, warnings::Vector{LinterWarning})
    # Structure compatible with atom linter
    # see python's `dolo-lint --format=json examples/models/rbc_iid.yaml` for an example
    # [{"type": "error", "source": "examples/models/rbc_iid.yaml", "range": [[11, 29], [11, 33]], "text": "Symbol 'beta' already declared as 'parameters'. (pos (12, 16))"}]
end


text:
function applyendbcs!(system::CVSystem,n::Int64,terms::Vector{Int64},hemoflag="no")
    # update periphery
    if hemoflag == "no"
        CVModule.coupledistal!(system,n,terms);
    elseif hemoflag == "yes"
        CVModule.coupledistal!(system,n,terms,hemoflag);
    end
    for i = 1:length(terms)
        CVModule.updateterms!(system,n,terms[i]);
    end
    CVModule.updatevc!(system,n,terms);
    if hemoflag == "no"
        CVModule.updateliver!(system,n);
    elseif hemoflag == "yes"
        CVModule.updateliver!(system,n,hemoflag);
    end

    # update right heart
    CVModule.updaterh!(system,n);

    # update pulmonary circulation
    CVModule.updatelungs!(system,n);

    # update left heart
    CVModule.updatela!(system,n);
    CVModule.coupleproximal!(system,n);

    # update heart chamber pressures
    CVModule.elastancemodel!(system,n+1);
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
using Plots
gr()

using LaTeXStrings

using Unitful
using Unitful: N, mm, inch, m, cm, MPa
Unitful.preferunits(u"mm, s"...)

using DataFrames

using LinearAlgebra

θ = LinRange(0, 2π, 100)
y = sin.(θ)

p1 = plot(θ, y,
	title = L"sin(\theta)",
	xlabel = L"\theta",
	ylabel = L"Amplitude",
	xlim = [0, 2π],
	ylim = [-1, 1],
	legend = false
);

display(p1)

p2 = scatter(θ, y,
	alpha = 0.5,
	color = :jet1,
	marker_z = y,
	marker = :hexagon,
	markersize = 5,
	legend = false
);
display(p2)

F = [1000, 3200, -4000, 4300, -2750, 6630]N
r = [300mm, 400mm, 10inch, 15cm, 0.188m, 27cm] .|> mm

df = DataFrame(r=r, F=F)


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

__precompile__(true)

module RadiationSpectra

using DelimitedFiles
using Statistics

using Distributions
using Optim
using RecipesBase
using StatsBase

import Base: print, println, show

export FitFunction
export get_ndims, get_nparams
export set_fitranges!
export set_parameter_names!
export set_initial_parameters!
export get_fitted_parameters
export get_initial_parameters

export llhfit!, lsqfit!

export peakfinder
export calibrate_spectrum

"""
    get_example_spectrum()::Histogram

Returns an uncalibrated radiation spectrum for testing and demonstrating purpose.
"""
function get_example_spectrum()::Histogram
    weights::Array{Int, 1} = readdlm(joinpath(@__DIR__, "..", "examples", "data", "hpge-spectrum.txt"))[:]
    h = Histogram(0:100:100 * length(weights), :left)
    h.weights = weights
    return h
end

include("PeakFinder/PeakFinder.jl")

include("Fitting/FitFunction.jl")
include("Fitting/general_model_functions.jl")

include("Fitting/llhfit.jl")

include("Fitting/lsqfit.jl")

include("AutoCalibration/AutoCalibration.jl")

end # module


text:
Rectangular() = Uniform(-0.5, +0.5)

function _string(kernel::Uniform)
    @unpack a, b = kernel
    "Rectangular kernel (U[$a,$b])"
end

const LocationScaleDists = Union{Cosine,Epanechnikov,SymTriangularDist,Triweight}

_string(kernel::SymTriangularDist) = "Triangular kernel"
_string(kernel::LocationScaleDists) = "$(Base.typename(typeof(kernel))) kernel"

const SupportedKernels = Union{LocationScaleDists,Uniform}

function setbandwidth(kernel::LocationScaleDists, h::Number)
    typeof(kernel)(zero(kernel.μ), h * kernel.σ)
end

function setbandwidth(kernel::Uniform, h::Number)
    @unpack a, b = kernel
    if a + b != 0
        error("Use symmetric Uniform/Rectangular kernel")
    end
    Uniform(a * h, b * h)
end

weights(::Uniform, ZsR::RunningVariable) = uweights(length(ZsR))
weights(D::Distribution, ZsR::RunningVariable) = pdf.(D, ZsR)

struct EquivalentKernel{K,ET,T}
    kernel::K #only order 1 for now
    E::ET
    kernel_mult::T
end

function _expectation(ub, n)
    rawNodes, rawWeights = gausslegendre(n)
    _nodes = map(x -> (0.5 * ub) * x + ub / 2, rawNodes)
    _weights = map(x -> x * 1 / 2 * ub, rawWeights) # (result of doing 1/(b-a) * (b-a)/2)
    f -> dot(f.(_nodes), _weights)
end

function EquivalentKernel(kernel::SupportedKernels)
    @unpack ub = support(kernel)
    E = _expectation(ub, 10)
    cov_x = [x->1 x->x; x->x x->x^2]
    kernel_mult = [1 0] * inv(map(f -> E(u -> f(u) * pdf(kernel, u)), cov_x))
    EquivalentKernel(kernel, E, kernel_mult)
end

Distributions.support(eq::EquivalentKernel) = eq.kernel

function Distributions.pdf(eq::EquivalentKernel, u)
    @unpack kernel_mult, kernel = eq
    dot(kernel_mult, [1; u]) * pdf(kernel, u)
end

function kernel_moment(kernel, ::Val{j}) where {j}
    @unpack ub = support(kernel)
    quadgk(x -> x^j * pdf(kernel, x), 0, ub)[1]
end

kernel_moment(::SupportedKernels, ::Val{0}) = 0.5
kernel_moment(::SupportedKernels, ::Val{1}) = 0
kernel_moment(::SupportedKernels, ::Val{3}) = 0

function kernel_moment(eq::EquivalentKernel, ::Val{j}) where {j}
    @unpack E = eq
    E(x -> x^j * pdf(eq, x))
end

function squared_kernel_moment(kernel, ::Val{j}) where {j}
    @unpack ub = support(kernel)
    quadgk(x -> x^j * pdf(kernel, x)^2, 0, ub)[1]
end

squared_kernel_moment(::SupportedKernels, ::Val{1}) = 0
squared_kernel_moment(::SupportedKernels, ::Val{3}) = 0

function squared_kernel_moment(eq::EquivalentKernel, ::Val{j}) where {j}
    @unpack E = eq
    E(x -> x^j * pdf(eq, x)^2)
end


text:
require("mnist.jl")
using CUDArt
using KUnet
using MNIST: xtrn, ytrn, xtst, ytst
accuracy(y,z)=mean(findmax(y,1)[2] .== findmax(z,1)[2])
KUnet.srandom(1)

net = [Conv(5,5,1,20), Bias(20), Relu(), Pool(2),
       Conv(5,5,20,50), Bias(50), Relu(), Pool(2),
       Mmul(500,800), Bias(500), Relu(),
       Mmul(10,500), Bias(10), XentLoss()]

@show x = CudaArray(reshape(xtrn[:,1:64], 28, 28, 1, 64))
for l in net
    @show x = KUnet.forw(l, x)
end
@show dy = CudaArray(ytrn[:,1:64])
@show loss(net[length(net)], dy)
for i=length(net):-1:1
    @show dy = KUnet.back(net[i], dy)
end

xtrn = reshape(xtrn, 28, 28, 1, size(xtrn, 2))
xtst = reshape(xtst, 28, 28, 1, size(xtst, 2))
setparam!(net; :lr=0.01)
for i=1:100
    train(net, xtrn, ytrn)
    println((i, accuracy(ytst, predict(net, xtst)),
                accuracy(ytrn, predict(net, xtrn))))
end

# lr = 0.01

# # x1:(28,28,1,64) using 1 channel and 64 minibatch size
# @show x1 = Tensor(reshape(MNIST.xtrn[:,1:64], 28, 28, 1, 64))

# dims1 = (5,5,1,20)  # using 5x5 filters, 1 input, 20 outputs
# conv1 = ConvLayer(w=Filter(float32(randn(dims1)*0.01)),
#                   b=Tensor(zeros(Float32,(1,1,dims1[4],1))),
#                   pd=PoolingDescriptor((2,2)),
#                   pw=UpdateParam(learningRate=lr),
#                   pb=UpdateParam(learningRate=2*lr),
#                   )

# # y1: (24,24,20,64)
# # z1: (12,12,20,64) => x2

# @show x2=KUnet.forw(conv1,x1)

# dims2 = (5,5,20,50)
# conv2 = ConvLayer(w=Filter(float32(randn(dims2)*0.01)),
#                   b=Tensor(zeros(Float32,(1,1,dims2[4],1))),
#                   pd=PoolingDescriptor((2,2)),
#                   pw=UpdateParam(learningRate=lr),
#                   pb=UpdateParam(learningRate=2*lr),
#                   )

# # y2: (8,8,50,64)
# # z2: (4,4,50,64)
# # x3: (800,64)

# @show x3=KUnet.forw(conv2,x2)

# dims3 = (500,800)
# ip3 = Op(w=CudaArray(float32(randn(dims3)*0.01)),
#             b=CudaArray(zeros(Float32, dims3[1], 1)),
#             f=relu,
#             pw=UpdateParam(learningRate=lr),
#             pb=UpdateParam(learningRate=2*lr),
#             )

# # y3: (500,64) => x4

# @show x4=KUnet.forw(ip3,x3)

# dims4 = (10,500)
# ip4 = Op(w=CudaArray(float32(randn(dims4)*0.01)),
#             b=CudaArray(zeros(Float32, dims4[1], 1)),
#             f=relu,
#             pw=UpdateParam(learningRate=lr),
#             pb=UpdateParam(learningRate=2*lr),
#             )

# @show x5=KUnet.forw(ip4,x4)

# net4 = [conv1, conv2, ip3, ip4]


text:
module Defs

using StaticArrays

export invariant
function invariant end

if VERSION < v"1.5"
    Base.:≈(x) = y -> y ≈ x
end

end


text:
abstract type AbstractAdvection end

struct UpWind2ndOrder <: AbstractAdvection
end

function (a::UpWind2ndOrder)(rhs, ρC, buvec, uvec, Tc, Tbf, ∇Tc, bf2c, bfn, bcv, f2c, fn, cv, fc, ccenter, floops::Face2CellLoop)
    NBF = nbfaces(floops)
    @inbounds for i=1:NBF
        j = bf2c[i][1]
        rhs[j] -= (ρC * Tbf[i] * (buvec[i] ⋅ bfn[i]))*bcv[i]
    end
    NF = nfaces(floops)
    @inbounds for i=1:NF
        j1, j2 = f2c[i]
        n = fn[i]
        v1, v2 = cv[i]
        udotn = uvec[i] ⋅ n
        if udotn >= 0
            Tf = Tc[j1] + (fc[i] - ccenter[i][1])⋅∇Tc[j1]
        else
            Tf = Tc[j2] + (fc[i] - ccenter[i][2])⋅∇Tc[j2]
        end
        qf = ρC*Tf*udotn
        rhs[j1] -= qf*v1
        rhs[j2] += qf*v2
    end

end

function (a::UpWind2ndOrder)(rhs,p)
    floops = p.mesh.f2cloops
    a(rhs, p.ρC, p.u, p.u, p.Tc, p.Tbf, p.∇Tc, floops.bf2c, floops.bfn, floops.bcv, floops.f2c, floops.fn, floops.cv, floops.fc, floops.ccenter, floops)
end

struct FullUpWind <: AbstractAdvection
end

function (a::FullUpWind)(rhs, ρC, buvec, uvec, Tc, Tbf, bf2c, bfn, bcv, f2c, fn, cv, fc, ccenter, floops::Face2CellLoop)
    NBF = nbfaces(floops)
    @inbounds for i=1:NBF
        j = bf2c[i][1]
        rhs[j] -= (ρC * Tbf[i] * (buvec[i] ⋅ bfn[i]))*bcv[i]
    end
    NF = nfaces(floops)
    @inbounds for i=1:NF
        j1, j2 = f2c[i]
        n = fn[i]
        v1, v2 = cv[i]
        udotn = uvec[i] ⋅ n
        if udotn >= 0
            Tf = Tc[j1]
        else
            Tf = Tc[j2]
        end
        qf = ρC*Tf*udotn
        rhs[j1] -= qf*v1
        rhs[j2] += qf*v2
    end

end

function (a::FullUpWind)(rhs,p)
    floops = p.mesh.f2cloops
    a(rhs, p.ρC, p.u, p.u, p.Tc, p.Tbf, floops.bf2c, floops.bfn, floops.bcv, floops.f2c, floops.fn, floops.cv, floops.fc, floops.ccenter, floops)
end

struct UpWindAdvection{FT,Fbf,VT,UBT,fl}
    f::FT
    fbf::Fbf
    an::FT
    anm1::FT
    u::VT
    ubf::UBT
    floops::fl
end

function UpWindAdvection(f,fbf,u,ubf,mesh)
    an = similar(f)
    anm1 = similar(an)
    fill!(an,zero(eltype(an)))
    fill!(anm1,zero(eltype(an)))
    fl = mesh.f2cloops
    return UpWindAdvection{typeof(f),typeof(fbf),typeof(u),typeof(ubf),typeof(fl)}(f,fbf,an,anm1,u,ubf,fl)
end

function mycopy!(rhs,inp)
    @inbounds for i in linearindices(rhs)
        rhs[i] = inp[i]
    end
end

function (Up::UpWindAdvection)()

    # Copy An to An-1
    mycopy!(Up.anm1,Up.an)

    rhs = Up.an
    fill!(rhs,zero(eltype(rhs)))

    floops = Up.floops
    bf2c = floops.bf2c
    Tbf = Up.fbf
    buvec = Up.ubf
    bfn = floops.bfn
    bcv = floops.bcv

    NBF = nbfaces(floops)
    @inbounds for i=1:NBF
        j = bf2c[i][1]
        rhs[j] += (Tbf[i] * (buvec[i] ⋅ bfn[i]))*bcv[i]
    end

    f2c = floops.f2c
    fn = floops.fn
    cv = floops.cv
    u = Up.u
    Tc = Up.f

    NF = nfaces(floops)
    @inbounds for i=1:NF
        j1, j2 = f2c[i]
        n = fn[i]
        v1, v2 = cv[i]
        uvec = 0.5*(u[j1]+u[j2])
        udotn = uvec ⋅ n
        if udotn >= 0
            Tf = Tc[j1]
        else
            Tf = Tc[j2]
        end
        qf = Tf*udotn
        rhs[j1] += qf*v1
        rhs[j2] -= qf*v2
    end

end

function (Up::UpWindAdvection)(rhs)
    an = Up.an
    anm1 = Up.anm1

    @inbounds for i in linearindices(rhs)
        rhs[i] -= 1.5*an[i] - 0.5*anm1[i]
        #rhs[i] -= an[i]
    end

end


text:
#==============================================================================#
#
#       ODE solver for H,F,c terms; relies on DifferentialEquations.jl
#
#==============================================================================#
"""
    struct HFcSolver{Tmode,Tsv,Tps,Tcb,Ts,T,Ta} <: AbstractGuidingTermSolver{Tmode}
        saved_values::Tsv
        problem_setup::Tps
        callback::Tcb
        solver::Ts
        HFc0::T
        access::Ta
    end

Struct for solving a (H,F,c)-system of ODEs.

        HFcSolver(
            ::Val{:inplace},
            tt,
            xT_plus,
            P,
            obs,
            choices,
        )
    Constructor for an ODE solver with in-place operations. Pre-allocates space
    and solves a (H,F,c)-system of ODEs once, backward in time on the interval
    `(tt[1], tt[end])`, with a terminal condition computed from `xT_plus`. `P`
    is the auxiliary diffusion law, `obs` is the observation made at time
    `tt[end]` and `choices` contains information about eltype and the algorithm
    for solving ODEs. H,F,c are saved on a grid of time-points `tt`.

        HFcSolver(
            ::Val{:outofplace},
            tt,
            xT_plus,
            P,
            obs,
            choices
        )
    Constructor for an ODE solver with out-of-place operations using
    StaticArrays. Initialises the object and solves a (H,F,c)-system of ODEs
    once, backward in time on the interval `(tt[1], tt[end])`, with a terminal
    condition computed from `xT_plus`. `P` is the auxiliary diffusion law, `obs`
    is the observation made at time `tt[end]` and `choices` contains information
    about eltype and the algorithm for solving ODEs. H,F,c are saved on a grid
    of time-points `tt`.
"""
struct HFcSolver{Tmode,Tsv,Ti,TT,T0,Ta} <: AbstractGuidingTermSolver{Tmode}
    saved_values::Tsv
    integrator::Ti
    HFcT::TT
    HFc0::T0
    access::Ta

    function HFcSolver(
            ::Val{:inplace},
            tt,
            xT_plus,
            P,
            obs,
            choices,
        )
        D = DiffusionDefinition.dimension(P).process
        el = choices.eltype

        function HFc_update!(du, u, p, t)
            # shorthand names for views, hopefully optimised
            # ---
            # current state
            H, F, c = u.H, u.F, u.c
            B, β, σ, a, tmat = p.buf.B, p.buf.β, p.buf.σ, p.buf.a, p.buf.mat
            tvec, P = p.buf.vec, p.P
            # increments (to-be-computed by this function)
            dH, dF, dc = du.H, du.F, du.c
            # ---

            # in-place computation of auxiliary process; stored in du
            DD.B!(B, t, P), DD.β!(β, t, P), DD.σ!(σ, t, P), DD.a!(a, t, P)

            # ODEs
            # ---
            # dH = - (_B' * H) - (H * _B) + outer(H * σ)
            mul!(tmat, H, a)
            mul!(dH, tmat, H')
            mul!(dH, H, B, -true, true)
            mul!(dH, B', H, -true, true)
            # dF = - (_B' * F) + (H * a * F) + (H * _β)
            mul!(dF, tmat, F)
            mul!(dF, B', F, -true, true)
            mul!(dF, H, β, true, true)
            # dc = dot(_β, F) + 0.5*outer(F' * _σ) - 0.5*sum(H .* _a)
            # apparently, according to @time, a faster way to compute a trace,
            # with no allocations
            dc[1] = 0.0
            for i in 1:D
                dc[1] -= tmat[i,i]
            end
            dc[1] *= 0.5
            #@time dc[1] = -0.5*tr(tmat)
            mul!(tvec, a', F)
            mul!(dc, tvec', F, 0.5, true)
            mul!(dc, β', F, true, true)
            # ---
        end

        HFcT = HFcContainer(el, D)
        update_HFc!(HFcT, xT_plus, obs)

        prob = ODEProblem{true}(
            HFc_update!,
            HFcT,
            (tt[end], tt[1]),
            (
                buf = HFcBuffer{el}(D),
                P = P,
            )
        )
        saved_values = SavedValues(Float64, Tuple{Matrix{el},Vector{el},el})
        callback = SavingCallback(
            # NOTE this might need changing, even though it seems to be working
            # correctly; from
            # https://docs.sciml.ai/latest/features/callback_library/ :
            # "save_func(u, t, integrator) returns the quantities which shall be
            # saved. Note that this should allocate the output (not as a view to
            # u)". Currently u.H & u.F provide views!
            (u,t,integrator)->(u.H, u.F, u.c[1]),
            saved_values;
            saveat=reverse(tt),
            tdir=-1
        )
        integrator = init(
            prob,
            choices.solver;
            callback=callback,
            save_everystep=false, # to prevent wasting memory allocations
        )
        sol = solve!(integrator)
        HFc0 = sol.u[end]
        Tsv, Ti = typeof(saved_values), typeof(integrator), typeof(HFc0)
        TT, T0 = typeof(HFcT), typeof(HFc0)
        new{:inplace,Tsv,Ti,TT,T0,Nothing}(
            saved_values,
            integrator,
            HFcT,
            HFc0,
            nothing
        )
    end

    function HFcSolver(
            ::Val{:outofplace},
            tt,
            xT_plus,
            P,
            obs,
            choices
        )
        access = Val{DiffusionDefinition.dimension(P).process}()

        function HFc_update(u, P, t)
            H, F, c = static_accessor_HFc(u, access)
            _B, _β, _σ, _a = DD.B(t, P), DD.β(t, P), DD.σ(t, P), DD.a(t, P)

            dH = - (_B' * H) - (H * _B) + outer(H * _σ)
            dF = - (_B' * F) + (H * _a * F) + (H * _β)
            dc = dot(_β, F) + 0.5*outer(F' * _σ) - 0.5*tr(H*_a)
            vcat(SVector(dH), dF, SVector(dc))
        end

        el = choices.eltype
        TH, TF, Tc = prepare_static_saving_types(Val{:hfc}(), access, el)
        prob = ODEProblem{false}(
            HFc_update,
            update_HFc(xT_plus, obs, access),
            (tt[end], tt[1]),
            P
        )
        saved_values = SavedValues(Float64, Tuple{TH,TF,Tc})
        callback = SavingCallback(
            (u,t,integrator) -> static_accessor_HFc(u, access),
            saved_values;
            saveat=reverse(tt),
            tdir=-1
        )
        integrator = init(
            prob,
            choices.solver,
            callback=callback,
            save_everystep=false, # to prevent wasting memory allocations
        )
        sol = solve!(integrator)
        HFc0 = MVector(sol.u[end])
        Tsv, Ti, T0 = typeof(saved_values), typeof(integrator), typeof(HFc0)
        Ta = typeof(access)
        new{:outofplace,Tsv,Ti,Nothing,T0,Ta}(
            saved_values,
            integrator,
            nothing,
            HFc0,
            deepcopy(access),
        )
    end
end

"""
    HFc0(s::HFcSolver)

Return the data containing H,F,c terms for the time 0+
"""
HFc0(s::HFcSolver) = s.HFc0

"""
    update_HFc!(u_T, u_Tplus, obs, access)

Update equations for H,F,c at the times of observations. Save the data into
`u_T`.
"""
function update_HFc!(u_T, u_Tplus, obs)
    L, Λ, Σ, v, μ = OBS.L(obs), OBS.Λ(obs), OBS.Σ(obs), OBS.ν(obs), OBS.μ(obs)
    m, d = size(L)
    u_T.H .= u_Tplus.H + L'*Λ*L
    u_T.F .= u_Tplus.F + L'*Λ*v
    u_T.c .= (u_Tplus.c .+ 0.5*( m*log(2π) + log(det(Σ)) + (v-μ)'*Λ*(v-μ) ))
end

"""
    update_HFc(u_Tplus, obs, access)

Update equations for H,F,c at the times of observations.
"""
function update_HFc(u_Tplus, obs, access)
    L, Λ, Σ, v, μ = OBS.L(obs), OBS.Λ(obs), OBS.Σ(obs), OBS.ν(obs), OBS.μ(obs)
    m, _ = size(L)
    H, F, c = static_accessor_HFc(u_Tplus, access)
    dH = L'*Λ*L
    dF = L'*Λ*v
    dc = 0.5*( m*log(2π) + log(det(Σ)) + (v-μ)'*Λ*(v-μ) )
    vcat(SVector(H + dH), (F + dF), SVector(c+dc))
end

"""
    prepare_static_saving_types(::Val{:hfc}, ::Val{D}, el) where D

Define data-types for H,F,c computed by out-of-place solver that are to be saved
internally.
"""
function prepare_static_saving_types(::Val{:hfc}, ::Val{D}, el) where D
    SMatrix{D,D,el,D*D}, SVector{D,el}, el
end

"""
    H(s::HFcSolver, i::Integer)

Return saved matrix H[i] (with H[1] indicating H at time 0+ and H[end]
indicating H at time T).
"""
H(s::HFcSolver, i::Integer) = s.saved_values.saveval[end-i+1][1]

"""
    F(s::HFcSolver, i::Integer)

Return saved vector F[i] (with F[1] indicating F at time 0+ and F[end]
indicating F at time T).
"""
F(s::HFcSolver, i::Integer) = s.saved_values.saveval[end-i+1][2]

"""
    c(s::HFcSolver, i::Integer)

Return saved scalar c[i] (with c[1] indicating c at time 0+ and c[end]
indicating c at time T).
"""
c(s::HFcSolver, i::Integer) = s.saved_values.saveval[end-i+1][3]

function recompute_guiding_term!(
        s::HFcSolver{:inplace},
        P,
        obs,
        xT_plus
    )
    update_HFc!(s.HFcT, xT_plus, obs)
    reinit!(s.integrator, s.HFcT)
    sol = solve!(s.integrator)
    s.HFc0 .= sol.u[end]
end

function recompute_guiding_term!(
        s::HFcSolver{:outofplace},
        P,
        obs,
        xT_plus
    )
    reinit!(s.integrator, update_HFc(xT_plus, obs, s.access))
    sol = solve!(s.integrator)
    s.HFc0 .= sol.u[end]
end

∇logρ(i, x, P::HFcSolver{:outofplace}) = F(P, i) - H(P, i)*x

function ∇logρ!(buffer, i, x, P::HFcSolver{:inplace})
    mul!(buffer.∇logρ, H(P, i), x, -true, false)
    _F = F(P, i)
    for j in 1:length(buffer.∇logρ)
        buffer.∇logρ[j] += _F[j]
    end
end

function proto_mutable_obs(guiding_term_solver::HFcSolver{:outofplace})
    _s2m_vec(F(guiding_term_solver, 1))
end

_s2m_vec(v::SVector{S,T}) where {S,T} = MVector{S,T}(v)

function proto_mutable_obs(guiding_term_solver::HFcSolver{:inplace})
    SVector{1,Float64}(1.0) # let's worry about this later, doesn't matter for now
end

mode(::HFcSolver{Tmode}) where Tmode = Tmode


text:
using Sialia
using Test

@testset "Sialia.jl" begin
    include("eventdispatcher_tests.jl")
    include("eventstore_tests.jl")
end


text:
export initNet

function initNet(layersNames,layersDims,netType=:classification)
  net=Layer[]
  if netType==:classification
    for (i,layerName) in enumerate(layersNames)
      if layerName==:d
        net=[net;dataLayer(layersDims[i,:]...)]
      elseif layerName==:h
        net=[net;hiddenLayer(layersDims[i,:]...)]
      elseif layerName==:p
        net=[net;predictLayer(layersDims[i,:]...)]
      elseif layerName==:l
        net=[net;lossLayer(layersDims[i,1])]
        net[i].Cost=crossEntropy;net[i].∂Cost=∂crossEntropy
      end
    end
  elseif netType==:regression
    for (i,layerName) in enumerate(layersNames)
      if layerName==:d
        net=[net;dataLayer(layersDims[i,:]...)]
      elseif layerName==:h
        net=[net;hiddenLayer(layersDims[i,:]...)]
      elseif layerName==:p
        net=[net;predictLayer(layersDims[i,:]...)]
      elseif layerName==:l
        net=[net;lossLayer(layersDims[i,1])]
      end
    end
  end
  return net
end


text:
import Pluto
Pluto.run()


text:
using EnhancedGJK
import GeometryTypes
const gt = GeometryTypes
using BenchmarkTools
using MeshIO, FileIO  # v0.3.2
using CoordinateTransformations, StaticArrays
using MeshCat, Colors

# Dimensions: (0.264 x 0.14 x 0.0885 m)
T = Float64
c = []
x, y, z = 0.05, 0.05, 0.01
push!(c, load(string(pwd(),"/test/meshes/r_foot_chull.obj")))
push!(c, gt.GLNormalMesh(gt.HyperSphere(gt.Point(0.,0,z), 0.005)))  # size is one-sigma sensor noise
poseA = Translation(SVector{3,T}(0,0,0))
poseB = Translation(SVector{3,T}(0,0,0))

# Calculate
result = []
push!(result, gjk(c[1], c[2], poseA, poseB))

# Additional points
push!(c, gt.GLNormalMesh(gt.HyperSphere(gt.Point(0,y,z), 0.005)))
push!(result, gjk(c[1], c[3], poseA, poseB))

push!(c, gt.GLNormalMesh(gt.HyperSphere(gt.Point(x,y,z), 0.005)))
push!(result, gjk(c[1], c[4], poseA, poseB))

# Construct surfaces
T = Float64
env_point = SVector{3,SVector{3,T}}(SVector{3,T}(0,0,z), SVector{3,T}(0,y,z), SVector{3,T}(x,y,z))  # inertial frame
vertices = poseA.(gt.vertices(c[1]))  # inertial frame
signed_distances = SVector{3,T}(result[1].signed_distance,
                                result[2].signed_distance,
                                result[3].signed_distance)
coeff, centers = implicit_surface(vertices, env_point, signed_distances)
# centers, distances = construct_centers(vertices, env_point, signed_distances)
# coeff = construct_surface(centers, distances)  # body frame

# Benchmark
@btime construct_centers(vertices, env_point, signed_distances)
@btime construct_surface(centers, distances)
@btime implicit_surface(vertices, env_point, signed_distances)

# Evaluate surface and check with actual gjk results
eval_surface(coeff, centers, env_point[1])

# Differentiation
pt = [0, 0, 0.01]
∇f = similar(pt)
surface_gradient!(∇f, coeff, centers, pt)
print(∇f)
# ∇f, H = surface_gradient(coeff, centers, SVector{3,T}(0,0,0.1))
# @btime surface_gradient(coeff, centers, SVector{3,T}(0,0,0.1))

# Test intersection set algorithm with sampled points around points in collision
zero_set, cnt = fzero(coeff, centers, [env_point[1]+[0,0,0.01]], ∇f)
cnt
eval_surface(coeff, centers, zero_set[1])
# eval_surface(coeff, centers, zero_set[2])
# eval_surface(coeff, centers, zero_set[3])

# Create implicit surface mesh
# surf_pts, surf_faces = make_geom(coeff, centers, SVector{3,T}(1,1,1), SVector{3,T}(2,2,2), 20)
# mesh = make_mesh(surf_pts, surf_faces)
mesh = make_mesh_sdf(coeff, centers, SVector{3,T}(-0.5,-0.5,-0.5), SVector{3,T}(2,2,2), 0.01)

# # get vertices, normals, and faces to create HomogenousMesh
# _vertices = GeometryBasics.decompose(GeometryBasics.Point{3,T}, mesh)
# _normals = GeometryBasics.decompose_normals(mesh)
# _faces = GeometryBasics.decompose(GeometryBasics.TriangleFace{Int}, mesh)
# _texture = [gt.Point(0.0f0,0) for i in 1:1692]
# implicit_mesh = gt.HomogenousMesh(_vertices, _normals, _faces, texture)

# Visualization
vis = Visualizer()
open(vis)

setobject!(vis["mesh/1"], c[1], MeshPhongMaterial(color=RGBA{Float32}(1.0,1.0,0.9,0.9)))
setobject!(vis["mesh/2"], c[2])
setobject!(vis["mesh/3"], c[3])
setobject!(vis["mesh/4"], c[4])
setobject!(vis["mesh/surface"], mesh, MeshPhongMaterial(color=RGBA{Float32}(0.0,0.0,1.0,0.9)))

settransform!(vis["mesh/1"], poseA)
settransform!(vis["mesh/2"], poseB)
settransform!(vis["mesh/3"], poseB)
settransform!(vis["mesh/4"], poseB)
settransform!(vis["mesh/surface"], poseA)

# Add gjk arrows
dist_vec = ArrowVisualizer(vis["dist/a"])
setobject!(dist_vec, MeshLambertMaterial(color=RGBA{Float32}(0,1,0,0.5)))
settransform!(dist_vec, gt.Point(poseA(result.closest_point_in_body.a)),
    gt.Point(poseB(result.closest_point_in_body.b)))  # poseB transformation

# delete tasks
delete!(vis["mesh/1"])
delete!(vis["mesh/2"])
delete!(vis["mesh/3"])
delete!(vis["mesh/4"])
delete!(vis["mesh/surface"])

# delete visualizer
delete!(vis)


text:
using MultiComponentFlash
using Documenter
using BenchmarkTools, Plots

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

makedocs(;
    modules=[MultiComponentFlash],
    authors="Olav Møyner <olav.moyner@gmail.com>",
    repo="https://github.com/moyner/MultiComponentFlash.jl/blob/{commit}{path}#{line}",
    sitename="MultiComponentFlash.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://moyner.github.io/MultiComponentFlash.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Examples" => Any[
            "Basic usage" => "examples/basics.md",
            "Advanced usage" => "examples/advanced.md"
        ],
        "API" => Any[
            "Mixtures" => "api/mixtures.md",
            "Equilibrium constants" => "api/kvalues.md",
            "Equations of state" => "api/eos.md",
            "Vapor-liquid equilibrium (Flash)" => "api/flash.md",
            "Utilities" => "api/utilities.md",
        ],
    ],
)

deploydocs(;
    repo="github.com/moyner/MultiComponentFlash.jl",
    devbranch="main",
)


text:
function normalization_jet_mixing(settings::PyObject, V_j, rho_j, A_j, Tt_j, c_0, rho_0, T_0)

    V_j_star = V_j ./ c_0
    rho_j_star = rho_j ./ rho_0
    A_j_star = A_j ./ settings.A_e
    Tt_j_star = Tt_j ./ T_0

    return V_j_star, rho_j_star, A_j_star, Tt_j_star
end

function normalization_jet_shock(settings::PyObject, V_j, A_j, Tt_j, c_0, T_0)

    V_j_star = V_j ./ c_0
    A_j_star = A_j ./ settings.A_e
    Tt_j_star = Tt_j ./ T_0

    return V_j_star, A_j_star, Tt_j_star
end

function normalization_jet(settings::PyObject, V_j, rho_j, A_j, Tt_j, c_0, rho_0, T_0)

    V_j_star = V_j ./ c_0
    rho_j_star = rho_j ./ rho_0
    A_j_star = A_j ./ settings.A_e
    Tt_j_star = Tt_j ./ T_0

    return V_j_star, rho_j_star, A_j_star, Tt_j_star
end

function normalization_core_ge(settings::PyObject, mdoti_c, Tti_c, Ttj_c, Pti_c, DTt_des_c, c_0, rho_0, T_0, p_0)

    mdoti_c_star = mdoti_c ./ (rho_0 .* c_0 .* settings.A_e)
    Tti_c_star = Tti_c ./ T_0
    Ttj_c_star = Ttj_c ./ T_0
    Pti_c_star = Pti_c ./ p_0
    DTt_des_c_star = DTt_des_c ./ T_0

    return mdoti_c_star, Tti_c_star, Ttj_c_star, Pti_c_star, DTt_des_c_star
end

function normalization_core_pw(settings::PyObject, mdoti_c, Tti_c, Ttj_c, Pti_c, rho_te_c, c_te_c, rho_ti_c, c_ti_c, c_0, rho_0, T_0, p_0)

    mdoti_c_star = mdoti_c / (c_0 * settings.A_e * rho_0)
    Tti_c_star = Tti_c ./ T_0
    Ttj_c_star = Ttj_c ./ T_0
    Pti_c_star = Pti_c ./ p_0
    rho_te_c_star = rho_te_c ./ rho_0
    c_te_c_star = c_te_c ./ c_0
    rho_ti_c_star = rho_ti_c ./ rho_0
    c_ti_c_star = c_ti_c ./ c_0

    return mdoti_c_star, Tti_c_star, Ttj_c_star, Pti_c_star, rho_te_c_star, c_te_c_star, rho_ti_c_star, c_ti_c_star
end

function normalization_fan(settings::PyObject, DTt_f ,mdot_f, N_f, A_f, d_f, c_0, rho_0, T_0)

    DTt_f_star = DTt_f ./ T_0
    mdot_f_star = mdot_f ./ (rho_0 .* c_0 .* settings.A_e)
    N_f_star = N_f ./ (c_0 ./ d_f * 60)
    A_f_star = A_f ./ settings.A_e
    d_f_star = d_f ./ sqrt(settings.A_e)

    return DTt_f_star, mdot_f_star, N_f_star, A_f_star, d_f_star
end

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA


text:
## ConstDict
"""
An immutable `AbstractDict` for which the value is the same, no matter what the key is.
"""
struct ConstDict{K, V} <: AbstractDict{K, V}
    val::V
    ConstDict{K}(val::V) where {K, V} = new{K, V}(val)
end
Base.getindex(d::ConstDict{K}, key) where K = d.val
Base.show(io::IO, ::MIME"text/plain", d::ConstDict{K}) where {K} = show(io, d)
Base.show(io::IO, d::ConstDict{K}) where {K} = print(io, "$(typeof(d)) with fixed value $(d.val)")


text:
using ITensors,
      Test

@testset "IndexSet" begin
  idim = 2
  jdim = 3
  kdim = 4
  ldim = 5
  i = Index(idim,"i")
  j = Index(jdim,"j")
  k = Index(kdim,"k")
  l = Index(ldim,"l")
  @testset "show" begin
      indices = (i, j, k, l)
      indset = IndexSet(indices)
      inner = sprint.(show, indices)
      @test sprint(show, indset) == string("IndexSet(", join(inner, ", "), ")")
  end
  @testset "Basic ctors" begin
    I = IndexSet(i,j,k)
    @test IndexSet(I) === I
    @test l ∈ IndexSet(I, l)
    @test l ∈ IndexSet(l, I)
    @test l ∈ IndexSet( (I, IndexSet(l)) )
  end
  @testset "Convert to Index" begin
    @test Index(IndexSet(i)) === i
    @test_throws ErrorException Index(IndexSet(i, j))
  end
  @testset "Index dimensions" begin
    I = IndexSet(i,j,k)
    @test dim(I) == idim*jdim*kdim
    @test dims(I) == (idim,jdim,kdim)
    @test dim(I,1) == idim
    @test dim(I,2) == jdim
    @test dim(I,3) == kdim
  end
  @testset "Set operations" begin
    I1 = IndexSet(i,j,k)
    I2 = IndexSet(k,l)
    I3 = IndexSet(j,l)
    @test hassameinds(I1,(k,j,i))
    @test uniqueindex(I1,(I2,I3)) == i
    @test uniqueindex(I1,IndexSet(k, j, i)) == Index()
    @test uniqueinds(I1,I2) == IndexSet(i,j)
    @test setdiff(I1,I2) == IndexSet(i,j)
    @test hassameinds(uniqueinds(I1,I2),(j,i))
    @test commoninds(I1,I2) == IndexSet(k)
    @test commonindex(I1,I2) == k
    @test commonindex(I1,IndexSet(l)) == Index()
    @test commoninds(I1,(j,l)) == IndexSet(j)
    @test commonindex(I1,(j,l)) == j
    @test commoninds(I1,(j,k)) == IndexSet(j,k)
    @test hassameinds(commoninds(I1,(j,k,l)),(j,k))
    @test findinds(I1,"i") == IndexSet(i)
    @test findindex(I1,"j") == j
    @test findindex(I1,"l") == Index()
  end
  @testset "commoninds index ordering" begin
    I = IndexSet(i,k,j)
    J = IndexSet(j,l,i)
    # Test that commoninds respects the ordering
    # of the indices in the first IndexSet
    @test commoninds(I,J) == IndexSet(i,j)
    @test commoninds(J,I) == IndexSet(j,i)
  end
  @testset "adjoint" begin
    I = IndexSet(i,k,j)
    @test adjoint(I) == IndexSet(i', k', j')
  end
  @testset "mapprime" begin
    I = IndexSet(i',k'',j)
    @test mapprime(I,1,5) == IndexSet(i^5,k'',j)
    @test mapprime(I,2,0) == IndexSet(i',k,j)

    J = IndexSet(i,j,k')
    @test mapprime(J,0,2) == IndexSet(i'',j'',k')

    mapprime!(J,1,5)
    @test J == IndexSet(i,j,k^5)
  end
  @testset "compute_strides" begin
    I = IndexSet(i, j)
    @test ITensors.compute_strides(I) == [1, idim]
  end
end


text:
module Nodes

using AxisIndices
using NamedDims
using FieldProperties
using GeometryBasics
using MetadataArrays
using NeuroCore.SpatialAPI
using NeuroCore.AnatomicalAPI

export PointNode, PointNodeList, MultiPointNodeList

"""
    CoordinateProperties
"""
struct CoordinateProperties{L,N,Ax,CS,P} <: AbstractPropertyList{P}
    sc::SpatialCoordinates{L,N,Ax,CS}
    properties::P
end

SpatialCoordinates(cp::CoordinateProperties) = getfield(cp, :sc)

FieldProperties.properties(cp::CoordinateProperties) = getfield(cp, :properties)

AxisIndices.axes_type(::Type{<:CoordinateProperties{L,N,Ax,CS,P}}) where {L,N,Ax,CS,P} = Ax

NamedDims.dimnames(::Type{<:CoordinateProperties{L,N,Ax,CS,P}}) where {L,N,Ax,CS,P} = L

Base.ndims(::Type{<:CoordinateProperties{L,N,Ax,CS,P}}) where {L,N,Ax,CS,P} = N

Base.axes(cp::CoordinateProperties) = axes(SpatialCoordinates(cp))

function SpatialAxes.coordinate_system(cp::CoordinateProperties)
    return coordinate_system(SpatialCoordinates(cp))
end

include("pointnodes.jl")

end


text:
"""
    AbstractDiscretizationAlgorithm

Abstract supertype for all discretization algorithms.

### Note

For implementing a custom discretization algorithm, a type definition
`struct NewDiscretizationAlgorithm <: AbstractDiscretizationAlgorithm end`
and a `_discretize` method

```julia
_discretize(::NewDiscretizationAlgorithm, ΔT::Real,
             A::AbstractMatrix, B::AbstractMatrix, c::AbstractVector, D::AbstractMatrix)
```
are required.
"""
abstract type AbstractDiscretizationAlgorithm end

"""
    ExactDiscretization <: AbstractDiscretizationAlgorithm

Exact discretization algorithm for affine systems.

### Algorithm

This algorithm consists of integrating the continuous differential equation over
a specified time interval to obtain an associated difference equation.
The algorithm applies to any system of the form `x' = Ax + Bu + c + Dw` where the
state matrix `A` is invertible, and other system types, e.g. linear systems
`x' = Ax` which are included in the above formulation.

Without loss of generality, consider a `NoisyAffineControlledContinuousSystem`
with system dynamics ``x' = Ax + Bu + c + Dw``.

The exact discretization is calculated by integrating on both sides of the
continuous ODE over the time span `[t, t + ΔT]`, for a fixed input `u` and fixed
noise realization `w` at time `t`.
The resulting discretization writes as
``x^+ = A^d x + B^d u + c^d  + D^d w`` where
``A^d = \\exp^{A ~ ΔT}``,
``B^d = A^{-1}(A^d - I)B``,
``c^d = A^{-1}(A^d - I)c`` and
``D^d = A^{-1}(A^d - I)D``.

The algorithm described above is a well known result from the literature [1].

[1] https://en.wikipedia.org/wiki/Discretization#Discretization_of_linear_state_space_models
"""
struct ExactDiscretization <: AbstractDiscretizationAlgorithm end

"""
    EulerDiscretization <: AbstractDiscretizationAlgorithm

Euler discretization algorithm for affine systems.

### Algorithm

This algorithm consists of a first-order approximation to obtain an associated
difference equation. The algorithm applies to any system of the form
`x' = Ax + Bu + c + Dw`, and other system types, e.g. linear systems `x' = Ax`
which are included in the above formulation.

Without loss of generality, consider a `NoisyAffineControlledContinuousSystem`
with system dynamics ``x' = Ax + Bu + c + Dw``.

The Euler discretization is calculated by taking the first-order approximation
of the exact discretization [`ExactDiscretization`](@ref).
The resulting discretization writes as
``x^+ = A^d x + B^d u + c^d + D^d w``
where  ``A^d = I + ΔT ~ A``, ``B^d = ΔT ~ B``,
``c^d = ΔT ~ c`` and ``D^d = ΔT ~ D``.

The algorithm described above is a well known result from the literature [1].

[1] https://en.wikipedia.org/wiki/Discretization#Approximations
"""
struct EulerDiscretization <: AbstractDiscretizationAlgorithm end

"""
    discretize(system::AbstractContinuousSystem, ΔT::Real,
               algorithm::AbstractDiscretizationAlgorithm=ExactDiscretization(),
               constructor=_default_complementary_constructor(system))

Discretization of a `isaffine` `AbstractContinuousSystem` to a
`AbstractDiscreteSystem` with sampling time `ΔT` using the discretization
method `algorithm`.

### Input

- `system`      -- an affine continuous system
- `ΔT`          -- sampling time
- `algorithm`   -- (optional, default: `ExactDiscretization()`) discretization algorithm
- `constructor` -- (optional, default: `_default_complementary_constructor(system)`) construction method

### Output

Returns a discretization of the input system `system` with discretization method
`algorithm` and sampling time `ΔT`.
"""
function discretize(system::AbstractContinuousSystem, ΔT::Real,
                    algorithm::AbstractDiscretizationAlgorithm=ExactDiscretization(),
                    constructor=_default_complementary_constructor(system))

    (!isaffine(system)) && throw(ArgumentError("system needs to be affine"))

    sets(x) = x ∈ [:X,:U,:W]
    matrices(x) = x ∈ [:A,:B,:b,:c,:D]

    # get all fields from system
    fields = collect(fieldnames(typeof(system)))

    # get the fields of `system` that are parameters of the affine system dynamics
    # i.e., all fields that need to be discretized
    cont_values = [getfield(system, f) for f in filter(matrices, fields)]

    # compute discretized values of dynamics_params_c
    disc_values = _discretize(algorithm, ΔT, cont_values...)
    # get the fields of `system` that are sets
    set_values = [getfield(system, f) for f in filter(sets, fields)]

    # build the new discrete type according to the constructor method with the
    # discretized and set values
    return constructor(disc_values, set_values)
end

# returns a function with signature 'Vector × Vector -> AbstractDiscreteSystem'
function _default_complementary_constructor(system)
    # get corresponding discrete type of `system`
    discrete_type = _complementary_type(typename(system))
    # build the new discrete type with the discretized and set values
    return (disc_values, set_values) -> discrete_type(disc_values..., set_values...)
end

"""
    _discretize(::AbstractDiscretizationAlgorithm, ΔT::Real
                A::AbstractMatrix, B::AbstractMatrix, c::AbstractVector, D::AbstractMatrix)

Implementation of the discretization algorithm defined by the first input argument
with sampling time `ΔT`.

### Input

- ``   -- discretization algorithm, used for dispatch
- `ΔT` -- sampling time
- `A`  -- state matrix
- `B`  -- input matrix
- `c`  -- vector
- `D`  -- noise matrix

### Output

Returns a vector containing the discretized input arguments `A`, `B`, `c` and `D`.

### Notes

See [`discretize`](@ref) for more details.
"""
function _discretize(::AbstractDiscretizationAlgorithm, ΔT::Real,
                     A::AbstractMatrix,
                     B::AbstractMatrix,
                     c::AbstractVector,
                     D::AbstractMatrix)
end

function _discretize(::ExactDiscretization, ΔT::Real,
                     A::AbstractMatrix,
                     B::AbstractMatrix,
                     c::AbstractVector,
                     D::AbstractMatrix)
    if  rank(A) == size(A, 1)
        A_d = exp(A*ΔT)
        Matr = inv(A)*(A_d - I)
        B_d = Matr*B
        c_d = Matr*c
        D_d = Matr*D
    else
        error("exact discretization for singular state matrix, i.e. A is non-invertible,"*
              " not implemented yet, please use algorithm `EulerDiscretization`")
    end
    return [A_d, B_d, c_d, D_d]
end

function _discretize(::EulerDiscretization, ΔT::Real,
                     A::AbstractMatrix,
                     B::AbstractMatrix,
                     c::AbstractVector,
                     D::AbstractMatrix)
    A_d = I + ΔT*A
    B_d = ΔT*B
    c_d = ΔT*c
    D_d = ΔT*D
    return [A_d, B_d, c_d, D_d]
end

"""
    _discretize(A::AbstractMatrix, ΔT::Real; algorithm=:exact)

Discretize the state matrix `A` with sampling time `ΔT` and discretization method
`algorithm`.

### Input

- `A`         -- state matrix
- `ΔT`        -- sampling time
- `algorithm` -- (optional, default: `:exact`) discretization algorithm

### Output

Returns a vector containing the discretized input argument `A`.

### Notes

See [`discretize`](@ref) for more details.
"""
function _discretize(algorithm::AbstractDiscretizationAlgorithm, ΔT::Real,
                     A::AbstractMatrix)
    n = size(A,1)
    mzero = spzeros(n, n)
    vzero = spzeros(n)
    A_d, _, _, _ = _discretize(algorithm, ΔT, A, mzero, vzero, mzero)
    return [A_d]
end

"""
    _discretize(algorithm::AbstractDiscretizationAlgorithm, ΔT::Real,
                A::AbstractMatrix, B::AbstractMatrix)

Discretize the state matrix `A` and input or noise matrix `B` with sampling time
`ΔT` and discretization method `algorithm`.

### Input

- `algorithm` -- discretization algorithm
- `ΔT`        -- sampling time
- `A`         -- state matrix
- `B`         -- input or noise matrix

### Output

Returns a vector containing the discretized input arguments `A` and `B`.

### Notes

This method signature with two arguments of type `AbstractMatrix` works both for
a noisy system with fields `(:A,:D)` and a controlled system with fields `(:A,:B)`.

See [`discretize`](@ref) for more details.
"""
function _discretize(algorithm::AbstractDiscretizationAlgorithm, ΔT::Real,
                     A::AbstractMatrix, B::AbstractMatrix)
    n = size(A,1)
    mzero = spzeros(n, n)
    vzero = spzeros(n)
    A_d, B_d, _, _ = _discretize(algorithm, ΔT, A, B, vzero, mzero)
    return [A_d, B_d]
end

"""
    _discretize(algorithm::AbstractDiscretizationAlgorithm, ΔT::Real,
                A::AbstractMatrix,c::AbstractVector)

Discretize the state matrix `A` and vector `c` with sampling time `ΔT` and
discretization method `algorithm`.

### Input

- `algorithm` -- discretization algorithm
- `ΔT`        -- sampling time
- `A`         -- state matrix
- `c`         -- vector

### Output

Returns a vector containing the discretized input arguments `A` and `c`.

### Notes

See [`discretize`](@ref) for more details.
"""
function _discretize(algorithm::AbstractDiscretizationAlgorithm, ΔT::Real,
                     A::AbstractMatrix,c::AbstractVector)
    n = size(A,1)
    mzero = spzeros(n, n)
    A_d, _, c_d, _ = _discretize(algorithm, ΔT, A, mzero, c, mzero)
    return [A_d, c_d]
end

"""
    _discretize(algorithm::AbstractDiscretizationAlgorithm, ΔT::Real,
                A::AbstractMatrix, B::AbstractMatrix, c::AbstractVector)

Discretize the state matrix `A`, input matrix `B` and vector `c` with sampling
time `ΔT` and discretization method `algorithm`.

### Input

- `algorithm` -- discretization algorithm
- `ΔT`        -- sampling time
- `A`         -- state matrix
- `B`         -- input matrix
- `c`         -- vector

### Output

Returns a vector containing the discretized input arguments `A`, `B` and `c`.

### Notes

See [`discretize`](@ref) for more details.
"""
function _discretize(algorithm::AbstractDiscretizationAlgorithm, ΔT::Real,
                     A::AbstractMatrix, B::AbstractMatrix, c::AbstractVector)
    n = size(A,1)
    mzero = spzeros(n, n)
    A_d, B_d, c_d, _ = _discretize(algorithm, ΔT, A, B, c, mzero)
    return [A_d, B_d, c_d]
end

"""
    _discretize(algorithm::AbstractDiscretizationAlgorithm, ΔT::Real,
                A::AbstractMatrix, B::AbstractMatrix, D::AbstractMatrix)

Discretize the state matrix `A`, input matrix `B` and noise matrix `C` with
sampling time `ΔT` and discretization method `algorithm`.

### Input

- `algorithm` -- discretization algorithm
- `ΔT`        -- sampling time
- `A`         -- state matrix
- `B`         -- input matrix
- `D`         -- noise matrix

### Output

Returns a vector containing the discretized input arguments `A`, `B` and `D`.

### Notes

See [`discretize`](@ref) for more details.
"""
function _discretize(algorithm::AbstractDiscretizationAlgorithm, ΔT::Real,
                     A::AbstractMatrix, B::AbstractMatrix, D::AbstractMatrix)
    n = size(A,1)
    vzero = spzeros(n)
    A_d, B_d, _, D_d = _discretize(algorithm, ΔT, A, B, vzero, D)
    return [A_d, B_d, D_d]
end


text:
module TestKmers

using Kmers
using BioSequences
using Test

const GROUP = get(ENV, "GROUP", "All")

include("utils.jl")

if GROUP == "BioSequences" || GROUP == "All"
    include("biosequences_interface.jl")
    include("conversion.jl")
    include("comparisons.jl")
    include("length.jl")
    include("access.jl")
    include("random.jl")
    include("find.jl")
    include("print.jl")
    include("transformations.jl")
    include("mismatches.jl")
    include("debruijn_neighbors.jl")
    #include("shuffle.jl")
end

end # module


text:
"""
    ActiveFilteredLogger(filter, logger)

Wraps `logger` in an active filter.
While loggers intrinsictally have in built filtering mechanisms.
Wrapping it in a `ActiveFilterLogger` allows for extract control,
at the cost of a bit of overhead.

The `ActiveFilteredLogger` has full control of what is logged,
as it sees the full message,
this does mean however it determines what to log at runtime,
which is the source of the overhead.
The [`EarlyFilteredLogger`](@ref) has less control,
but decides if to log before the message is computed.

The `filter` should be a function that returns a boolean.
`true` if the message should be logged and `false` if not.
As input it will be given a named tuple with the following fields:
`(level, message, _module, group, id, file, line, kwargs)`
See [`LoggingExtras.handle_message_args`](@ref) for more information on what each is.
"""
struct ActiveFilteredLogger{T <: AbstractLogger, F} <: AbstractLogger
    filter::F
    logger::T
end

function handle_message(logger::ActiveFilteredLogger, args...; kwargs...)
    log_args = handle_message_args(args...; kwargs...)
    if comp_handle_message_check(logger.logger, args...; kwargs...)
        if logger.filter(log_args)
            handle_message(logger.logger, args...; kwargs...)
        end
    end
end

function shouldlog(logger::ActiveFilteredLogger, args...)
    return comp_shouldlog(logger.logger, args...)
end

min_enabled_level(logger::ActiveFilteredLogger) = min_enabled_level(logger.logger)
catch_exceptions(logger::ActiveFilteredLogger) = catch_exceptions(logger.logger)

"""
    handle_message_args

This creates NamedTuple containing all the arguments the logger gives
to `handle_message`
It is the type pased to the active logger filter.
These argument come from the logging macro (@info`, `@warn` etc).

  * `level::LogLevel` Warn, Info, etc,
  * `message::String` the message to be logged
  * `_module::Module` can be used to specify a different originating module from
    the source location of the message.
  * `group::Symbol` can be used to override the message group (this is
    normally derived from the base name of the source file).
  * `id::Symbol` can be used to override the automatically generated unique
    message identifier.  This is useful if you need to very closely associate
    messages generated on different source lines.
  * `file::String` and `line::Int` can be used to override the apparent
    source location of a log message.
  * `kwargs...`: Any  keyword or position arguments passed to the logging macro
"""
function handle_message_args(args...; kwargs...)
    fieldnames = (:level, :message, :_module, :group, :id, :file, :line, :kwargs)
    fieldvals = (args..., kwargs)
    return NamedTuple{fieldnames, typeof(fieldvals)}(fieldvals)
end


text:
# # Implementing log/exponent transformations
# The CoordinateTransformations library provides a flexible and extensible framework
# for constructing coordinate transformations.

using CoordinateTransformations

struct LogTransformation{T <: Real} <: CoordinateTransformations.Transformation
    base
end

struct ExponentialTransformation{T <: Real} <: CoordinateTransformations.Transformation
    base::T
end

const ExpTransformation = ExponentialTransformation

function LogTransformation(; base::T = ℯ) where T
    return LogTransformation(base)
end

function ExponentialTransformation(; base::T = ℯ) where T
    return ExponentialTransformation(base)
end

(lt::LogTransformation)(x) = log.(lt.base, x)
(pt::ExponentialTransformation)(x) = lt.base .^ x

Base.inv(lt::LogTransformation) = ExponentialTransformation(lt.base)
Base.inv(pt::ExponentialTransformation) = LogTransformation(pt.base)

compose(lt::LogTransformation, pt::ExponentialTransformation) = lt.base == pt.base ? IdentityTransformation() : ComposedTransformation(lt, pt)
compose(pt::ExponentialTransformationlt::LogTransformation) = lt.base == pt.base ? IdentityTransformation() : ComposedTransformation(pt, lt)


text:
using Viruses

# a melt factor
ddf = 5.1
# calculate melt
meltrate(-1, ddf)


text:
using Test

# run plotting tests
include("plot_test.jl")

# run function tests
include("diagnostic_tests.jl")