text:
# ==========================================
#          REGULAR 2D SHAPES
# ==========================================

function make_rpolygon(
    xyz_centre::NTuple{3,Float64},
    theta::Float64,
    side_len::Float64,
    rf::Float64,
    num_sides::Real;
    normal_vec::NTuple{3,Float64} = (0.0, 0.0, 1.0)
)::Int
    num_sides = convert(Int64, num_sides)
    xc, yc, zc = xyz_centre
    ax, ay, az = normal_vec
    α::Float64 = π / num_sides
    circum_radius::Float64 = side_len / (2.0 * sin(α))
    h::Float64 = rf / cos(α)
    #
    pointsTags = []
    for asn in 1:num_sides
        θ_i = 2.0 * α * (asn - 1.0)
        o_dimTag = gmsh.model.occ.add_point(xc + circum_radius - h, yc, zc)
        p_dimTag = gmsh.model.occ.add_point(xc + circum_radius - h + rf * cos(α), yc + (rf * sin(α)), zc)
        q_dimTag = gmsh.model.occ.add_point(xc + circum_radius - h + rf * cos(α), yc - (rf * sin(α)), zc)
        gmsh.model.occ.rotate([(0, o_dimTag)], xc, yc, zc, ax, ay, az, θ_i)
        gmsh.model.occ.rotate([(0, p_dimTag)], xc, yc, zc, ax, ay, az, θ_i)
        gmsh.model.occ.rotate([(0, q_dimTag)], xc, yc, zc, ax, ay, az, θ_i)
        push!(pointsTags, (p_dimTag, o_dimTag, q_dimTag))
    end
    curvesTags = []
    for asn in 1:num_sides
        arcTag = gmsh.model.occ.add_circle_arc(pointsTags[asn]...)
        if asn <= num_sides - 1
            lineTag = gmsh.model.occ.add_line(pointsTags[asn][1], pointsTags[asn+1][3])
        elseif asn == num_sides
            lineTag = gmsh.model.occ.add_line(pointsTags[asn][1], pointsTags[1][3])
        end
        curvesTags = [curvesTags; [arcTag, lineTag]]
    end
    rpolygon_wireTag = gmsh.model.occ.add_curve_loop(curvesTags)
    rpolygon_discTag = gmsh.model.occ.add_plane_surface([rpolygon_wireTag])
    gmsh.model.occ.rotate((2, rpolygon_discTag), xc, yc, zc, ax, ay, az, theta)
    # gmsh.model.occ.synchronize()
    #
    return rpolygon_discTag
end

function make_capsule(
    xyz::NTuple{3,Float64},
    theta::Float64,
    smj_len::Float64,
    smn_len::Float64,
)::Int
    xc, yc, zc = xyz
    p1 = gmsh.model.occ.add_point(xc + smj_len - smn_len, yc + smn_len, zc)
    p2 = gmsh.model.occ.add_point(xc - smj_len + smn_len, yc + smn_len, zc)
    p3 = gmsh.model.occ.add_point(xc - smj_len + smn_len, yc - smn_len, zc)
    p4 = gmsh.model.occ.add_point(xc + smj_len - smn_len, yc - smn_len, zc)
    c1 = gmsh.model.occ.add_point(xc + smj_len - smn_len, yc, zc)
    c2 = gmsh.model.occ.add_point(xc - smj_len + smn_len, yc, zc)
    #
    line_1 = gmsh.model.occ.add_line(p1, p2)
    arc_1 = gmsh.model.occ.add_circle_arc(p3, c2, p2)
    line_2 = gmsh.model.occ.add_line(p3, p4)
    arc_2 = gmsh.model.occ.add_circle_arc(p1, c1, p4)
    #
    caps_wire_tag = gmsh.model.occ.add_curve_loop([line_1, arc_1, line_2, arc_2])
    caps_disc_tag = gmsh.model.occ.add_plane_surface([caps_wire_tag])
    gmsh.model.occ.rotate((2, caps_disc_tag), xc, yc, zc, 0.0, 0.0, 1.0, theta)
    return caps_disc_tag
end

function make_ellipse(
    xyz::NTuple{3,Float64},
    theta::Float64,
    smj_len::Float64,
    smn_len::Float64,
)::Int
    xc, yc, zc = xyz
    ell_curve_tag = gmsh.model.occ.add_ellipse(xc, yc, zc, smj_len, smn_len)
    ell_wire_tag = gmsh.model.occ.add_wire([ell_curve_tag])
    ell_disc_tag = gmsh.model.occ.add_plane_surface([ell_wire_tag])
    gmsh.model.occ.rotate((2, ell_disc_tag), xc, yc, zc, 0.0, 0.0, 1.0, theta)
    return ell_disc_tag
end

function make_rectangle(
    xyz::NTuple{3,Float64},
    theta::Float64,
    smj_len::Float64,
    smn_len::Float64,
    corner_radius::Float64
)::Int
    xc, yc, zc = xyz
    rect_tag = gmsh.model.occ.add_rectangle(xc - smj_len, yc - smn_len, zc, 2.0 * smj_len, 2.0 * smn_len, -1, corner_radius)
    gmsh.model.occ.rotate((2, rect_tag), xc, yc, zc, 0.0, 0.0, 1.0, theta)
    return rect_tag
end

# ==========================================
#          CLUSTERS OF REGULAR 2D SHAPES
# ==========================================

function add_circular_discs(
    xyr::Matrix{Float64},
    z_min::Float64,
)::Vector{Int}
    return [gmsh.model.occ.add_disk(ax, ay, z_min, ar, ar) for (ax, ay, ar) in eachrow(xyr)]
end

function add_elliptical_discs(
    xyt_ab::Matrix{Float64},
    z_min::Float64,
)::Vector{Int}
    @assert size(xyt_ab)[2] == 5 "fibre data matrix should contain five columns with xy coordinates of ellipse centre, theta, semi major and minor legths respectively"
    return [make_ellipse((axc, ayc, z_min), ath, aa, ab) for (axc, ayc, ath, aa, ab) in eachrow(xyt_ab)]
end

function add_capsular_discs(
    xyt_ab::Matrix{Float64},
    z_min::Float64,
)::Vector{Int}
    @assert size(xyt_ab)[2] == 5 "fibre data matrix should contain five columns with xy coordinates of capsule centre, theta, semi major and minor legths respectively"
    return [make_capsule((axc, ayc, z_min), ath, aa, ab) for (axc, ayc, ath, aa, ab) in eachrow(xyt_ab)]
end

function add_rectangular_discs(
    xyt_abr::Matrix{Float64},
    z_min::Float64,
)::Vector{Int}
    @assert size(xyt_abr)[2] == 6 "fibre data matrix should contain six columns with xy coordinates of rectangle centre, theta, semi major and minor legths, corner radius respectively"
    return [make_rectangle((ax, ay, z_min), ath, aa, ab, ar) for (ax, ay, ath, aa, ab, ar) in eachrow(xyt_abr)]
end

function add_rpolygons(
    xyt_arn::Matrix{Float64},
    z_min::Float64,
)::Vector{Int}
    @assert size(xyt_arn)[2] == 6 "fibre data matrix should contain six columns with xy coordinates of rpolygon centre, theta, side legnth, corner radius and number of sides respectively"
    return [make_rpolygon((ax, ay, z_min), ath, aa, ar, an) for (ax, ay, ath, aa, ar, an) in eachrow(xyt_arn)]
end

# ==========================================
#          REGULAR 3D SHAPES
# ==========================================

function make_ellipsoid(
    x::Float64, y::Float64, z::Float64,
    azm::Float64, plr::Float64,
    a::Float64, b::Float64, c::Float64,
)::Int
    dim_tag = (3, gmsh.model.occ.add_sphere(x, y, z, 1.0))
    gmsh.model.occ.dilate(dim_tag, x, y, z, a, b, c)  #CHECK_IT major axis or semi-major axis lengths
    gmsh.model.occ.rotate(dim_tag, x, y, z, 0.0, 0.0, 1.0, azm)
    gmsh.model.occ.rotate(dim_tag, x, y, z, sin(azm), -cos(azm), 0.0, (π*0.5) - plr)  # CHECK_IT if the orientations are as expected?
    return dim_tag[2]
end

# ==========================================
#          CLUSTERS OF REGULAR 3D SHAPES
# ==========================================

function add_spheres(
    xyzr::Matrix{Float64},
)::Vector{Int}
    return [
        gmsh.model.occ.add_sphere(ax, ay, az, ar) for (
            ax, ay, az, ar) in eachrow(xyzr)
    ]
end

"""
    add_ellipsoids() -> Vector{Int}

ellipsoids data must be `Matrix{Float64}` with
x, y, z, azimuthal angle, polar angle, smjrx, intrx, smnrx axis lengths,
as columns, in order.

"""
function add_ellipsoids(
    xyz_ap_abc::Matrix{Float64}
)
    @assert size(xyt_ab, 2) == 8 "Ellipsoids data matrix should contain 8 columns with,
     xyz coordinates of ellipsoid centre,
     azimuthal angle, polar angle, semi major, intermediate and semi minor legths respectively
     BUT $(size(xyt_ab, 2)) is found."
    return [make_ellipsoid(axc, ayc, azc, aazm, aplr, aa, ab, ac) for (
        axc, ayc, azc, aazm, aplr, aa, ab, ac
    ) in eachrow(xyz_ap_abc)]
end


text:
module CombinedParserTools
using CombinedParsers
using CombinedParsers.Regexp

import CombinedParsers.Regexp: word

footnote    = re"[⁰¹²³⁴⁵⁶⁷⁸⁹]+"
quotes      = re"[\"'`]"
enum_label = re"(?:[0-9]{1,3}|[ivx]{1,6}|[[:alpha:]])[\.\)]"
emptyline          = re"[ \t]*\r?\n"
whitespace_newline = re"[ \t]*\r?\n"
email_regexp = re"[-+_.~a-zA-Z][-+_.~:a-zA-Z0-9]*@[-.a-zA-Z0-9]+"
delimiter = re"[-, _/\.;:*\|!?&]"
capitalized = Sequence(CharIn(isuppercase), Repeat1(CharIn(islowercase)))

## is this official??
author_email = Sequence(
    :name => !Repeat(CharNotIn('<')),
    " <", :email => Repeat_until(email_regexp, ">"))

## export tokenizer_regex
##function tokenizer_regex()
##    re = (
newline     = re"\r?\n"
inline     = re"[^\r\n]"
lf          = re"\n"
whitespace  = re"[ \t]+"
# indentation = re"[ \t]*"
# content_characters = re"[^\t\r\n]+"
number      = re"[0-9]+"  ## TODO alt(...) csv
letters     = re"[A-Za-z*-]*"

parenthesisP(open,close) =
    Sequence(open, Repeat_until(AnyChar(),close,wrap=JoinSubstring))
pad(x) = Sequence(2,whitespace_maybe, x, whitespace_maybe)

extension   = re"\.[[:alnum:]~#]+"

include("tokens.jl")

############################################################
## legacy constructors

@deprecate opt(x...;kw...) Optional(x...;kw...)
@deprecate alt(a...) Either(a...)
# @deprecate seq(a...; kw...) Sequence(a...; kw...)
export seq
function seq(tokens::Vararg;
             transform=nothing, kw...)
    if transform isa Integer
        Sequence(transform,tokens...)
    elseif transform===nothing
        Sequence(tokens...)
    elseif transform isa Function
        map(transform, Sequence(tokens...))
    end
end

include("deprecated.jl")

include("greedy.jl")

import CombinedParsers: Numeric
export url
"""
Parse generic URL syntax, see [Wikipedia](https://en.wikipedia.org/wiki/Uniform_Resource_Identifier#Generic_syntax).
"""
url(;
    blacklist=re"[^\v\h]",
    scheme = !!re"[a-z]+",
    match_parse = !!Repeat(AnyChar()),
    host_parse = join(match_parse,'.'),
    path_parse = match_parse, #path_file(),
    query_parse = match_parse
    # ["http","https","ftp","ftps",
    #  "mailto","irc",
    #  "file","data"]
    ) =
        Sequence(:scheme => Atomic(scheme),
                 ':', Atomic(Either(
                     "//",
                     PositiveLookbehind(caseless("mailto:")))),
                 :authority => Sequence(
                     :userinfo => Optional(Sequence(1,!!re"\w+", '@')),
                     :host => map(host_parse, !(Repeat1(blacklist .& CharNotIn("/:")))),
                     ## TODO: parse IP here?
                     :port => (':',Numeric(Int64))[2]  | 80
                 ),
                 :path => Optional(Sequence(
                     2, '/', map(path_parse, !Repeat(blacklist .& CharNotIn("?#"))))),
                 :query => Optional(
                     Sequence(2,'?',   !!(Repeat(blacklist .& CharNotIn("#"))))), #   ∘ query_parse )),
                 :fragment => Optional(Sequence(2, '#', !!Repeat(blacklist))),
                 NegativeLookahead(CharIn(",!?.)]:")) ## not the best solution
                 )

_prefix(p,x...) = ( isempty(x) || x == tuple("") ) ? tuple() : (p,x...)
_suffix(x,p) = x == "" ? "" : x*p
Base.show(io::IO, x::NamedTuple{fieldnames(result_type(url()))}) =
    print(io, x.scheme,"://",
          _suffix(x.authority.userinfo,"@")...,
          x.authority.host...,
          x.authority.port == 80 ? "" : x.authority.port,
          "/", x.path...,
          _prefix("?", x.query, _prefix("#", x.fragment...)...)...
          )

end # module

ext: jl


text:
# %% imports

using PyCall
gym = pyimport("gym")

include("../../Modules/Gridcells/StateSpikeGaussGrid.jl")
include("../../Modules/Networks/net_arch.jl")
include("simviz.jl")

# %% Build the grids

include("../../Modules/Gridcells/GridGenerator.jl")

grid_cells = grid_gen()

# %% import the trained cartpole

torch = pyimport("torch")
nn    = pyimport("torch.nn")
f     = pyimport("torch.nn.functional")

@pydef mutable struct DQLnet <: nn.Module
    function __init__(self)
        pybuiltin(:super)(DQLnet, self).__init__()
        self.fc1 = nn.Linear(4,16)
        self.fc2 = nn.Linear(16,16)
        self.fc3 = nn.Linear(16,2)
    end

    function forward(self, s)
        s = f.relu(self.fc1(s))
        s = f.relu(self.fc2(s))
        s = self.fc3(s)
        return s
    end
end

Qnet = DQLnet()
Qnet.load_state_dict(torch.load("/Users/titou/Documents/PFE/RL-STDP/Julia/CartPole/saved_cartpole.net"))

# %% Function

function count_spikes(arch::Array{Int64}, spiked::Array{Int64}, layer::Int64)
    count_layer = zeros(arch[layer])
    interval = [sum(arch[1:layer-1])+1, sum(arch[1:layer])]
    spikes_in = filter(x-> interval[1] <= x <= interval[2], spiked)
    spikes_in .-= sum(arch[1:layer-1])
    count_layer[spikes_in] .+= 1
    return count_layer
end

function left_right(count_spikes::Array{Float64}, n_left::Int64, n_right::Int64)
    sl = sum(count_spikes[1:n_left])
    sr = sum(count_spikes[n_left+1:end])
    (sl > sr) ? action = 0 : action = 1
    if sl == sr
        action = rand(0:1)
    end
    return action, sl, sr
end

# %% Simulation

global s_hist = [[],[]]

function run_sim(net::Network, n_steps::Int64, render::Bool=false;
                 method::String = "snn")
    env = gym.make("CartPole-v1")
    env.seed(0)
    obs = env.reset()
    spike_history = Array{Int64}[]
    @inbounds for j in 1:n_steps
        count_tot = zeros(net.arch[end])
        spikes = input_spikes(obs, grid_cells)
        influence = 2
        @inbounds for msec in 1:10
            (msec==10) ? train=true : train=false
            spiked = step!(net, spikes, train)
            if msec >= length(net.arch)-1
                count_tot .+= influence .* count_spikes(net.arch, spiked, length(net.arch))
                influence *= 0.95
                push!(spike_history,spiked)
            end
        end
        if method == "snn"
            n_left = div(net.arch[end],2)
            n_right = net.arch[end] - n_left
            action, sl, sr = left_right(count_tot, n_left, n_right)
            push!(s_hist[1], sl)
            push!(s_hist[2], sr)
        elseif method == "trained"
            n_left = div(net.arch[end],2)
            n_right = net.arch[end] - n_left
            action = torch.argmax(Qnet(torch.tensor(obs).float())).item()
            _, sl, sr = left_right(count_tot, n_left, n_right)
            push!(s_hist[1], sl)
            push!(s_hist[2], sr)
        elseif method == "random"
            n_left = div(net.arch[end],2)
            n_right = net.arch[end] - n_left
            action = rand(0:1)
            _, sl, sr = left_right(count_tot, n_left, n_right)
            push!(s_hist[1], sl)
            push!(s_hist[2], sr)
        end
        obs_new, reward, done, _ = env.step(action)
        if render==true
            env.render()
        end
        obs = obs_new
    end
    return spike_history
end

parameters = YAML.load(open("/Users/titou/Documents/PFE/RL-STDP/Julia/Modules/cfg.yml"))

n = Network([120,60,80], parameters)
#randomize!(n)

# %% run the simulation

@time spike_history = run_sim(n, 200, method = "snn")

for _ in 1:9
    @time spike_history .= run_sim(n, 200, method = "snn")
end

# %% Plots difference
x = collect(1:length(s_hist[1]))
y = s_hist[2] .- s_hist[1]
plot(x,y,
     xlabel = "Time step cartpole",
     ylabel = "Sum spikes left minus right",
     title = "Difference between left or right",
     legend = false)

# %% Plot s_hist
x = collect(1:length(s_hist[1]))
plot(x,s_hist,
     xlabel = "Time step cartpole",
     ylabel = "Sum spikes left and right",
     lab = ["left" "right"],
     title = "Left or Right plot")

# %% spike history plot per layer
layer_spike_hist(n.arch, spike_history, 1)

# %% spike distribution for all layers
spike_histogram(n.arch, spike_history)

# %% weights visualization heatmap
weight_viz(n, "heatmap")

# %% weights visualization histogram
weight_viz(n, "histogram")

# %% visualize network connection
heatmap(n.connection)

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
using PdbTool

function print_help()
	@printf("\nScript to write out mappings/distances between a HMM and\na PDB file. Note that this is inefficient for\nbatch jobs.\n\n")
	@printf("%sUSAGE%s:  julia backmap.jl [pdb_file] [chain] [hmm_file] [output_file] [options...]\n\n",PdbTool.KRED,PdbTool.KRES)
	@printf("\tThis will parse the PDB file [pdb_file], then map [chain] to\n\tthe Hidden-Markov model in [hmm_file] and write the mapping\n\tto [output_file].\n")
	@printf("\tFormat:\tColumn 1 = HMM position\n\t\tColumn 2 = PDB Residue Name\n\t\tColumn 3 = Amino Acid\n\n")
	@printf("%sOPTIONS%s:\n\n",PdbTool.KGRN,PdbTool.KRES)
	@printf("\t-d [distance_output_file]\n")
	@printf("\tCalculate distances between mapped residues and write\n\tto [distance_output_file].\n\tFormat:\tColumn 1 = HMM Pos 1\n\t\tColumn 2 = HMM Pos 2\n\t\tColumn 3 = PDB Residue Name 1\n\t\tColumn 4 = PDB Residue Name 2\n\t\tColumn 5 = Distance in Angstrom\n\n")
@printf("\t-dt [distance_type]\n")
	@printf("\tDistance type used with -d option. Should be 'heavyMin' (distance\n\tbetween closest heavy atoms in residues, default) or\n\t'ca' (distance between CA atoms).\n")
end

# Parse Input

if length(ARGS)==1 && ARGS[1]=="-h"
	print_help()
	exit(0)
end

if length(ARGS)<4
	@printf("%sERROR:%s expected at least 4 arguments, %d given\n",PdbTool.KRED,PdbTool.KRES,length(ARGS))
	@printf("Usage: julia backmap.jl [pdb_file] [chain] [hmm_file] [output_file] [options...]\n")
	@printf("For a list of possible options: julia backmap.jl -h\n")
	exit(1)
end

pdb_file=ARGS[1]
chain=ARGS[2]
hmm_file=ARGS[3]
output_file=ARGS[4]

MAKE_DIST=false
dist_type="heavyMin"
dist_output_file=nothing

options=("-d","-dt")
if length(ARGS)>4
	l=5
	while l<=length(ARGS)
		if !(ARGS[l] in options)
			@printf("%sERROR:%s argument \"%s\" not recognized\n",PdbTool.KRED,PdbTool.KRES,ARGS[l])
			exit(1)
		end
		if ARGS[l]=="-d"
			MAKE_DIST=true
			if length(ARGS)==l
				@printf("%sERROR:%s please provide a output file for distances after -d option\n",PdbTool.KRED,PdbTool.KRES)
				exit(1)
			end
			if (ARGS[l+1] in options)
				@printf("%sERROR:%s please provide a output file for distances after -d option\n",PdbTool.KRED,PdbTool.KRES)
				exit(1)
			end
			dist_output_file=ARGS[l+1]
			l+=2
			continue
		end
		if ARGS[l]=="-dt"
			if length(ARGS)==l
				@printf("%sERROR:%s please provide a valid distance measure (ca or heavyMin) after -dt option\n",PdbTool.KRED,PdbTool.KRES)
				exit(1)
			end
			if (ARGS[l+1] in options)
				@printf("%sERROR:%s please provide a valid distance measure (ca or heavyMin) after -dt option\n",PdbTool.KRED,PdbTool.KRES)
				exit(1)
			end
			if !(ARGS[l+1] in ("ca","heaveMin"))
				@printf("%sERROR:%s please provide a valid distance measure (ca or heavyMin) after -dt option\n",PdbTool.KRED,PdbTool.KRES)
				exit(1)
			end
			dist_type=ARGS[l+1]
			l+=2
			continue
		end
		error("Parsing arguments failed")
	end
end

if !isfile(pdb_file)
	@printf("%sERROR: %sPDB file not found\n",PdbTool.KRED,PdbTool.KRES)

	exit(1)
end

if !isfile(hmm_file)
	@printf("%sERROR: %sHMM file not found\n",PdbTool.KRED,PdbTool.KRES)
	exit(1)
end

@printf("Parsing %s ...",basename(pdb_file))
pdb=PdbTool.parsePdb(pdb_file);
@printf("\%sSUCCESS\%s\n",PdbTool.KGRN,PdbTool.KRES)

if !haskey(pdb.chain,chain)
	@printf("%sERROR: %sPDB object has no chain %s\n",PdbTool.KRED,PdbTool.KRES,chain)
	exit(1)
end

@printf("Mapping %s to chain %s ...",basename(hmm_file),chain)
PdbTool.mapChainToHmm(pdb.chain[chain],hmm_file)
if length(pdb.chain[chain].align)==0
	@printf("\%sERROR\%s: no residue could be mapped\n",PdbTool.KRED,PdbTool.KRES)
	exit(1)
end
@printf("%sSUCCESS%s (mapped %d residues)\n",PdbTool.KGRN,PdbTool.KRES,length(pdb.chain[chain].align))

@printf("Writing mapping to \"%s\" ... ",output_file)
oFid=open(output_file,"w")
for l in sort(collect(keys(pdb.chain[chain].align)))
	id=pdb.chain[chain].align[l].identifier
	aa=pdb.chain[chain].align[l].aminoAcid
	@printf(oFid,"%d\t%s\t%s\n",l,id,aa)
end
close(oFid)
@printf("done.\n")

@printf("Writing distances to \"%s\" ... ",dist_output_file)
oFid=open(dist_output_file,"w")
if MAKE_DIST
	ind=sort(collect(keys(pdb.chain[chain].align)))
	for i=1:length(ind)-1
		for j=(i+1):length(ind)
			ri=pdb.chain[chain].align[ind[i]]
			rj=pdb.chain[chain].align[ind[j]]
			d=PdbTool.residueDist(ri,rj;distType=dist_type)
			@printf(oFid,"%d\t%d\t%s\t%s\t%f\n",ind[i],ind[j],ri.identifier,rj.identifier,d)
		end
	end
end
close(oFid)
@printf("done.\n")

max_stars_repo_name: christophfeinauer/PdbTool

max_issues_repo_name: christophfeinauer/PdbTool

max_forks_repo_name: christophfeinauer/PdbTool


text:
module AstroDates

using ..AstroTime: ASTRO_ISO_FORMAT

using Dates: Dates
using Dates: year, month, day, hour, minute, second
using Dates: millisecond, microsecond, nanosecond
using Dates: dayofyear

const J2000 = 2.4515445e6

function findyear(calendar, j2000day)
    j2kday = ifelse(j2000day isa Int32, widen(j2000day), j2000day)
    if calendar == :proleptic_julian
        return -((-4 * j2kday - 2920488) ÷ 1461)
    elseif calendar == :julian
        return -((-4 * j2kday - 2921948) ÷ 1461)
    end

    year = (400 * j2kday + 292194288) ÷ 146097

    # The previous estimate is one unit too high in some rare cases
    # (240 days in the 400 years gregorian cycle, about 0.16%)
    if j2kday <= last_j2000_dayofyear(:gregorian, year - 1)
        year -= 1
    end

    return year
end

function last_j2000_dayofyear(calendar, year)
    if calendar == :proleptic_julian
        return 365 * year + (year + 1) ÷ 4 - 730123
    elseif calendar == :julian
        return 365 * year + year ÷ 4 - 730122
    end

    return 365 * year + year ÷ 4 - year ÷ 100 + year ÷ 400 - 730120
end

function isleap(calendar, year)
    if calendar in (:proleptic_julian, :julian)
        return year % 4 == 0
    end

    return year % 4 == 0 && (year % 400 == 0 || year % 100 != 0)
end

const PREVIOUS_MONTH_END_DAY_LEAP = (0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335)

const PREVIOUS_MONTH_END_DAY = (0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334)

function findmonth(dayinyear, isleap)
    offset = ifelse(isleap, 313, 323)
    return ifelse(dayinyear < 32, 1, (10 * dayinyear + offset) ÷ 306)
end

function findday(dayinyear, month, isleap)
    (!isleap && dayinyear > 365) &&
        throw(ArgumentError("Day of year cannot be 366 for a non-leap year."))
    previous_days = ifelse(isleap, PREVIOUS_MONTH_END_DAY_LEAP, PREVIOUS_MONTH_END_DAY)
    return dayinyear - previous_days[month]
end

function finddayinyear(month, day, isleap)
    previous_days = ifelse(isleap, PREVIOUS_MONTH_END_DAY_LEAP, PREVIOUS_MONTH_END_DAY)
    return day + previous_days[month]
end

function getcalendar(year, month, day)
    if year < 1583
        if year < 1
            return :proleptic_julian
        elseif year < 1582 || month < 10 || (month < 11 && day < 5)
            return :julian
        end
    end
    return :gregorian
end

struct Date
    calendar::Symbol
    year::Int
    month::Int
    day::Int
end

function Date(offset::Integer)
    calendar = :gregorian
    if offset < -152384
        if offset > -730122
            calendar = :julian
        else
            calendar = :proleptic_julian
        end
    end

    year = findyear(calendar, offset)
    dayinyear = offset - last_j2000_dayofyear(calendar, year - 1)

    month = findmonth(dayinyear, isleap(calendar, year))
    day = findday(dayinyear, month, isleap(calendar, year))

    return Date(calendar, year, month, day)
end

Date(epoch::Date, offset::Int64) = Date(j2000(epoch) + offset)

function Date(year, month, day)
    if month < 1 || month > 12
        throw(ArgumentError("Invalid month number: $month"))
    end

    check = Date(j2000(year, month, day))
    if check.year != year || check.month != month || check.day != day
        throw(ArgumentError("Invalid date."))
    end
    return Date(calendar(check), year, month, day)
end

function Date(year, dayinyear)
    calendar = :gregorian
    if year < 1583
        if year < 1
            calendar = :proleptic_julian
        else
            calendar = :julian
        end
    end
    leap = isleap(calendar, year)
    month = findmonth(dayinyear, leap)
    day = findday(dayinyear, month, leap)

    return Date(calendar, year, month, day)
end

Dates.year(s::Date) = s.year
Dates.month(s::Date) = s.month
Dates.day(s::Date) = s.day
calendar(s::Date) = s.calendar

function Base.show(io::IO, d::Date)
    return print(io, year(d), "-", lpad(month(d), 2, '0'), "-", lpad(day(d), 2, '0'))
end

Date(d::Dates.Date) = Date(Dates.year(d), Dates.month(d), Dates.day(d))
Dates.Date(d::Date) = Dates.Date(year(d), month(d), day(d))

function j2000(calendar, year, month, day)
    d1 = last_j2000_dayofyear(calendar, year - 1)
    d2 = finddayinyear(month, day, isleap(calendar, year))
    return d1 + d2
end

function j2000(year, month, day)
    calendar = getcalendar(year, month, day)
    return j2000(calendar, year, month, day)
end

j2000(d::Date) = j2000(calendar(d), year(d), month(d), day(d))

julian(d::Date) = j2000(d) + J2000

const JULIAN_EPOCH = Date(-4712, 1, 1)
const MODIFIED_JULIAN_EPOCH = Date(1858, 11, 17)
const FIFTIES_EPOCH = Date(1950, 1, 1)
const CCSDS_EPOCH = Date(1958, 1, 1)
const GALILEO_EPOCH = Date(1999, 8, 22)
const GPS_EPOCH = Date(1980, 1, 6)
const J2000_EPOCH = Date(2000, 1, 1)
const MIN_EPOCH = Date(typemin(Int32))
const MAX_EPOCH = Date(typemax(Int32))
const UNIX_EPOCH = Date(1970, 1, 1)

struct Time{T}
    hour::Int
    minute::Int
    second::Int
    fraction::T

    function Time(hour, minute, second, fraction::T) where {T}
        if hour < 0 || hour > 23
            throw(ArgumentError("`hour` must be an integer between 0 and 23."))
        elseif minute < 0 || minute > 59
            throw(ArgumentError("`minute` must be an integer between 0 and 59."))
        elseif second < 0 || second >= 61
            throw(ArgumentError("`second` must be an integer between 0 and 61."))
        elseif fraction < 0 || fraction > 1
            throw(ArgumentError("`fraction` must be a number between 0 and 1."))
        end

        return new{T}(hour, minute, second, fraction)
    end
end

function Time(hour, minute, second)
    sec, frac = divrem(second, 1)
    return Time(hour, minute, sec, frac)
end

function Time(secondinday, fraction)
    if secondinday < 0 || secondinday > 86400
        throw(ArgumentError("Seconds are out of range"))
    end

    hour = secondinday ÷ 3600
    secondinday -= 3600 * hour
    minute = secondinday ÷ 60
    secondinday -= 60 * minute

    return Time(hour, minute, secondinday, fraction)
end

function Base.isapprox(a::Time, b::Time; kwargs...)
    return a.hour == b.hour &&
           a.minute == b.minute &&
           isapprox(a.fraction + a.second, b.fraction + b.second; kwargs...)
end

function subsecond(fraction, n, r)
    n % 3 == 0 || throw(ArgumentError("`n` must be divisible by 3."))
    factor = ifelse(Int === Int32, widen(10), 10)^n
    rounded = round(fraction, r; digits=n)
    return round(Int64, rounded * factor, r) % 1000
end

function subsecond(fraction, n)
    r = ifelse(subsecond(fraction, n+3, RoundNearest) == 0, RoundNearest, RoundToZero)
    return subsecond(fraction, n, r)
end

subsecond(t::Time, n) = subsecond(t.fraction, n)

Dates.hour(t::Time) = t.hour
Dates.minute(t::Time) = t.minute
Dates.second(::Type{Float64}, t::Time) = t.fraction + t.second
Dates.second(::Type{Int}, t::Time) = t.second
Dates.second(t::Time) = second(Int, t)
Dates.millisecond(t::Time) = subsecond(t.fraction, 3)
Dates.microsecond(t::Time) = subsecond(t.fraction, 6)
Dates.nanosecond(t::Time) = subsecond(t.fraction, 9)

fractionofday(t::Time) = (t.fraction + t.second) / 86400 + t.minute / 1440 + t.hour / 24
fractionofsecond(t::Time) = t.fraction

secondinday(t::Time) = t.fraction + t.second + 60 * t.minute + 3600 * t.hour

function Time(t::Dates.Time)
    return Time(
        Dates.hour(t),
        Dates.minute(t),
        Dates.second(t),
        1e-9 * Dates.nanosecond(t) +
        1e-6 * Dates.microsecond(t) +
        1e-3 * Dates.millisecond(t),
    )
end
function Dates.Time(t::Time)
    return Dates.Time(
        hour(t),
        minute(t),
        second(t),
        millisecond(t),
        microsecond(t),
        nanosecond(t),
    )
end

const H00 = Time(0, 0, 0, 0.0)
const H12 = Time(12, 0, 0, 0.0)

function Base.show(io::IO, t::Time)
    h = lpad(hour(t), 2, '0')
    m = lpad(minute(t), 2, '0')
    s = lpad(second(t), 2, '0')
    f = lpad(millisecond(t), 3, '0')
    return print(io, h, ":", m, ":", s, ".", f)
end

struct DateTime{T} <: Dates.AbstractDateTime
    date::Date
    time::Time{T}
end

Dates.default_format(::Type{DateTime}) = ASTRO_ISO_FORMAT[]

function Base.:(==)(a::DateTime, b::DateTime)
    return a.date == b.date && a.time == b.time
end

function Base.isapprox(a::DateTime, b::DateTime; kwargs...)
    return a.date == b.date && isapprox(a.time, b.time; kwargs...)
end

Date(dt::DateTime) = dt.date
Time(dt::DateTime) = dt.time

Base.show(io::IO, dt::DateTime) = print(io, Date(dt), "T", Time(dt))

function DateTime(str::AbstractString, df::Dates.DateFormat=Dates.default_format(DateTime))
    return Dates.parse(DateTime{Float64}, str, df)
end

function DateTime(
    year::Integer,
    month::Integer,
    day::Integer,
    hour::Integer=0,
    minute::Integer=0,
    second::Integer=0,
    fraction=0.0,
)
    return DateTime(Date(year, month, day), Time(hour, minute, second, fraction))
end

function DateTime(
    year::Integer, month::Integer, day::Integer, hour::Integer, minute::Integer, second
)
    return DateTime(Date(year, month, day), Time(hour, minute, second))
end

function DateTime(
    year::Int64,
    month::Int64,
    day::Int64,
    dayofyear::Int64,
    hour::Int64,
    minute::Int64,
    second::Int64,
    milliseconds::Int64,
    fractionofsecond::T,
) where {T}
    return DateTime{T}(
        year, month, day, dayofyear, hour, minute, second, milliseconds, fractionofsecond
    )
end

function DateTime{T}(
    year::Int64,
    month::Int64,
    day::Int64,
    dayofyear::Int64,
    hour::Int64,
    minute::Int64,
    second::Int64,
    milliseconds::Int64,
    fractionofsecond::T,
) where {T}
    if dayofyear != 0
        date = Date(year, dayofyear)
    else
        date = Date(year, month, day)
    end

    if !iszero(fractionofsecond)
        time = Time(hour, minute, second, fractionofsecond)
    else
        time = Time(hour, minute, second, 1e-3milliseconds)
    end

    return DateTime{T}(date, time)
end

Dates.year(dt::DateTime) = year(Date(dt))
Dates.month(dt::DateTime) = month(Date(dt))
Dates.day(dt::DateTime) = day(Date(dt))
Dates.yearmonthday(dt::DateTime) = year(Date(dt)), month(Date(dt)), day(Date(dt))
Dates.hour(dt::DateTime) = hour(Time(dt))
Dates.minute(dt::DateTime) = minute(Time(dt))
Dates.second(typ, dt::DateTime) = second(typ, Time(dt))
Dates.second(dt::DateTime) = second(Int, Time(dt))
Dates.millisecond(dt::DateTime) = millisecond(Time(dt))
Dates.microsecond(dt::DateTime) = microsecond(Time(dt))
Dates.nanosecond(dt::DateTime) = nanosecond(Time(dt))
fractionofsecond(dt::DateTime) = fractionofsecond(Time(dt))
subsecond(dt::DateTime, n) = subsecond(Time(dt), n)
calendar(dt::DateTime) = calendar(dt.date)

function Dates.dayofyear(dt::DateTime)
    leap = isleap(calendar(dt), year(dt))
    return finddayinyear(month(dt), day(dt), leap)
end

julian(dt::DateTime) = fractionofday(Time(dt)) + julian(Date(dt))
j2000(dt::DateTime) = fractionofday(Time(dt)) + j2000(Date(dt))
julian_twopart(dt::DateTime) = julian(Date(dt)), fractionofday(Time(dt))

function DateTime(dt::Dates.DateTime)
    return DateTime(
        Dates.year(dt),
        Dates.month(dt),
        Dates.day(dt),
        Dates.hour(dt),
        Dates.minute(dt),
        1e-3Dates.millisecond(dt) + Dates.second(dt),
    )
end
function Dates.DateTime(dt::DateTime)
    y = year(dt)
    m = month(dt)
    d = day(dt)
    h = hour(dt)
    mi = minute(dt)
    s = second(Int, dt)
    ms = millisecond(dt)
    return Dates.DateTime(y, m, d, h, mi, s, ms)
end

end

max_stars_repo_name: JuliaAstrodynamics/AstronomicalTime

max_issues_repo_name: JuliaAstrodynamics/AstronomicalTime

max_forks_repo_name: JuliaAstrodynamics/AstronomicalTime


text:
# Minimal encoder

encode(codec::MinimalCodec, x) = x

# BSON encoder

encode(codec::BSONCodec, x::BSONType) = x
encode(codec::BSONCodec, x::Symbol) = codec.as_string ? string(x) : BSONSymbol(x)
encode(codec::BSONCodec, x::AbstractVector{UInt8}) = Binary(0x00, x)
encode(codec::BSONCodec, x::Union{Int8,UInt8,Int16,UInt16,UInt32}) = Int32(x)
encode(codec::BSONCodec, x::AbstractFloat) = Float64(x)
encode(codec::BSONCodec, x::Integer) = Int64(x)

function encode(codec::BSONCodec, x::AbstractDict)
    Document([Element(string(k), encode(codec, v)) for (k, v) in x])
end

function encode(codec::BSONCodec, x::AbstractVector)
    BSONArray([Element(string(k - 1), encode(codec, v)) for (k, v) in enumerate(x)])
end


text:
using Documenter, ModernRoboticsBook

makedocs(;
    modules=[ModernRoboticsBook],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Manual" => Any[
            "man/examples.md",
        ],
        "Library" => Any[
            "Public" => "lib/public.md",
        ],
    ],
    repo="https://github.com/ferrolho/ModernRoboticsBook.jl/blob/{commit}{path}#L{line}",
    sitename="ModernRoboticsBook.jl",
    authors="Henrique Ferrolho",
    analytics = "UA-72743607-4",
    assets=[],
)

deploydocs(;
    repo="github.com/ferrolho/ModernRoboticsBook.jl",
)


text:
# The structure follows https://developer.spotify.com/documentation/web-api/reference/

export Albums, Artists, Browse, Episodes, Follow, Library
export Markets, Personalization, Player, Playlists, Search
export Shows, Tracks, UsersProfile, Objects
@info "Exported Library & etc."
"import export shorthand"
macro _ie(foo...)
    expr = :()
    for f in foo
        push!(expr.args, esc(
            quote
               import ..($f)
               export $f
            end))
    end
    expr
end
macro _ie()
    nothing
end

baremodule Albums
    import ..@_ie
    @_ie album_get_single album_get_multiple album_get_tracks album_get_saved album_get_contains
end
baremodule Artists
    import ..@_ie
    @_ie artist_get artist_get_albums artist_top_tracks artist_get_related_artists
end
baremodule Browse
    import ..@_ie
    @_ie category_get_single category_get_multiple
    @_ie category_get_new_releases
    @_ie recommendations_get  recommendations_dict_parser
end
baremodule Episodes
    import ..@_ie
    @_ie episodes_get_single episodes_get_multiple episodes_get_saved episodes_get_contains
end
baremodule Follow
    import ..@_ie
    @_ie follow_check follow_check_playlist follow_artists follow_artists_users follow_playlist unfollow_artists_users unfollow_playlist
end
baremodule Library
    import ..@_ie
    @_ie library_remove_albums library_remove_shows
    @_ie library_remove_tracks library_save_album library_save_show library_save_track
end
baremodule Markets
    import ..@_ie
    @_ie
end

baremodule Personalization
    import ..@_ie
    @_ie top_tracks top_artists
end
baremodule Player
    import ..@_ie
    @_ie player_get_state player_get_devices player_get_current_track player_get_recent_tracks
end
baremodule Playlists
    import ..@_ie
    @_ie playlist_get playlist_get_tracks playlist_get_current_user playlist_get_user playlist_get_featured
    @_ie playlist_get_category playlist_get_cover_image
end
baremodule Search
    import ..@_ie
    @_ie search_get
end
baremodule Shows
    import ..@_ie
    @_ie show_get_single show_get_multiple show_get_episodes show_get_saved show_get_contains
end
baremodule Tracks
    import ..@_ie
    @_ie tracks_get_audio_analysis tracks_get_audio_features tracks_get_single tracks_get_multiple
    @_ie tracks_get_saved tracks_get_contains
end
baremodule UsersProfile
    import ..@_ie
    @_ie user_get_profile user_get_current_profile
end
baremodule Objects
    import ..@_ie
    @_ie user_get_current_profile user_get_profile
end


text:
#=
Run-length encoding is a fast and simple method of encoding strings. The basic idea is to represent repeated successive characters as a single count and character. For example, the string "AAAABBBCCDAA" would be encoded as "4A3B2C1D2A".

Implement run-length encoding and decoding. You can assume the string to be encoded have no digits and consists solely of alphabetic characters. You can assume the string to be decoded is valid.
=#

function encode(decoded_string::String)
    encoded_string = ""
    while decoded_string != "1"
        character = decoded_string[1:1]
        number = 0
        while decoded_string[1:1] == character
            number = number + 1
            if length(decoded_string) >=2
                decoded_string = decoded_string[2:end]
            else
                decoded_string = "1"
            end
        end
        encoded_string = encoded_string * string(number) * character
    end
    return encoded_string
end

function decode(encoded_string::String)
    decoded_string = ""
    while encoded_string != "finished"
        decoded_string = decoded_string * repeat(encoded_string[2:2], parse(Int8, encoded_string[1:1]))
        if length(encoded_string) >2
            encoded_string = encoded_string[3:end]
        else
            encoded_string = "finished"
        end
    end
    return decoded_string
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
    function mime_type()
        executor(data["en"]["faker"]["file"]["mime_type"])
    end

        #param category: application|audio|image|message|model|multipart|text|video

    function file_extension()
        executor(data["en"]["faker"]["file"]["extension"])
    end

  function file_name( extencion="None")
        #:param extension: file extension

    if extencion == "None"
     extencion = file_extension()
    end

    filename = lexify("????")

    "$(filename).$(extencion)"

    end


text:
module LSR

import ..ED
using ..Levels
using ..SimLib
using ..SimLib: FArray
using SharedArrays: sdata

import Statistics

export levelspacingratio, LSRData, LSRDataDescriptor, center_region, load_lsr, LevelSpacingRatio

## Data structure

"""
    struct LSRDataDescriptor

Carries the information to construct a [`EDDataDescriptor`](!ref) object.
 - geometry
 - dimension
 - system_size
 - α
 - shots
 - densities ρs
 - field strengths
 - scaling of field strengths
 - symmetries to respect
For `load`ing data only the first 4 fields are required.
Can also be constructed from a `PositionDataDescriptor` by supplying a the missing bits (α, fields).
"""
struct LSRDataDescriptor <: ED.EDDerivedDataDescriptor
    derivedfrom::ED.EDDataDescriptor
end

LSRDataDescriptor(args...; kwargs...) = LSRDataDescriptor(EDDataDescriptor(args...; kwargs...))

"""
    struct LSRData

Stores the actual level-spacing ratios (LSR) specified by the descriptor [`LSRDataDescriptor`](@ref).
Level spacing ratio is defined as lsr_i = min(r_i, 1/r_i) where r_i = (E_(i-2) - E_(i-1))/(E_(i-1) - E_i)

# Index order
The indices mean:
 - i (as in lsr_i)
 - shot index
 - field index
 - density ρ

The default save directory is "lsr".

`Statistics.mean` and `Statistics.std` are overloaded to act on the first dimension to conveniently compute
mean LSR and its variance.
"""
struct LSRData{N} <: SimLib.AbstractSimpleData
    descriptor::LSRDataDescriptor
    # [dummy, shot, h, rho]
    data::FArray{N}
end

LSRData(lsrdd::LSRDataDescriptor) = LSRData(lsrdd, FArray{4}(undef, ED.ed_size(lsrdd), lsrdd.shots, length(lsrdd.fields), length(lsrdd.ρs)))

## Saving/Loading

ED._default_folder(::LSRDataDescriptor) = "lsr"

"""
    load_lsr(edd)
    load_lsr(model[, diagtype][, location])
"""
load_lsr(args...; kwargs...) = load(LSRDataDescriptor(args...; kwargs...))

function center_indices(L, center_region)
    cutoff = floor(Int, (L*(1-center_region)/2))
    (1+cutoff):L-cutoff
end

center_region(lsr::LSRData, center) = @view lsr.data[center_indices(size(lsr.data, 1), center), :, :, :]

function levelspacingratio(levels; center=1.0)
    sizes = size(levels)
    range = center_indices(sizes[1]-2, center) .+ 2
    res = Array{Float64, length(sizes)}(undef, length(range), sizes[2:end]...)
    for I in CartesianIndices(axes(levels)[2:end])
        for (i,j) in enumerate(range)
            a,b,c = levels[j-2,I], levels[j-1,I], levels[j,I]
            if a ≈ b || b ≈ c || a ≈ c
                # prevent NaNs or -Infs
                # a=b=c -> NaN
                # a>b=c -> -Inf # this can happen if all of them are very close
                res[i, I] = 0
            else
                ratio = abs((b-a)/(c-b)) # use abs to be sure
                res[i, I] = min(ratio, 1/ratio)
            end
        end
    end
    res
end

function SimLib.create(lsrdd::LSRDataDescriptor)
    leveldata = load_or_create(LevelDataDescriptor(lsrdd.derivedfrom))
    LSRData(lsrdd, levelspacingratio(leveldata.data))
end

mutable struct LSRTask <: ED.EDTask
    data
end

LevelSpacingRatio() = LSRTask(nothing)

function ED.initialize!(task::LSRTask, arrayconstructor, spectral_size)
    task.data = arrayconstructor(Float64, spectral_size-2)
end

function ED.compute_task!(task::LSRTask, evals, evecs, inds...)
    n = min(size(task.data,1), length(evals)-2)
    task.data[1:n, inds...] .= levelspacingratio(view(evals, 1:n+2))
end

function ED.failed_task!(task::LSRTask, inds...)
    task.data[:, inds...] .= NaN64
end

function ED.assemble(task::LSRTask, edd)
    LSRData(LSRDataDescriptor(edd), sdata(task.data))
end

###
## Mean LSR
###

mutable struct MeanLSRTask <: ED.EDTask
    data
end

MeanLevelSpacingRatio() = MeanLSRTask(nothing)

function ED.initialize!(task::MeanLSRTask, arrayconstructor, spectral_size)
    task.data = arrayconstructor(Float64, 1)
end

function ED.compute_task!(task::MeanLSRTask, evals, evecs, inds...)
    task.data[1, inds...] = mean(levelspacingratio(evals))
end

function ED.failed_task!(task::MeanLSRTask, inds...)
    task.data[1, inds...] .= NaN64
end

function ED.assemble(task::MeanLSRTask, edd)
    LSRData(LSRDataDescriptor(edd), sdata(task.data))
end

end #module


text:
using DifferentialEquations, Plots

function fhn!(du, u, p, t)
	du[1] = (u[1]-u[1]*u[1])*(u[1]-p[2]) - u[2]
	du[2] = p[3]*(p[1]*u[1] - u[2])
	return nothing
end

function lor!(du, u, p, t)
	du[1] = p[1]*(u[2]-u[1])
	du[2] = u[1]*(p[2]-u[3]) - u[2]
	du[3] = u[1]*u[2] - p[3]*u[3]
	du .= du .* 0.005
	return nothing
end

function f!(du, u, p, t)
	fhn!(view(du,1:2), view(u,1:2), view(p,1:3), t)
	lor!(view(du,3:5), view(u,3:5), view(p,4:6), t)
	#du[1] = du[1] + 0.05*(u[5]-27.0)
	return nothing
end

u0 = randn(Float64,2)
u1 = randn(Float64,3)

p0 = [0.37; 0.1; 0.01]
p1 = [10.0;28.0; 8/3]

tspan = (0.0, 5000.0)

prob0 = ODEProblem(fhn!, u0, tspan, p0)
prob1 = ODEProblem(lor!, u1, tspan, p1)
prob  = ODEProblem(f!, [u0;u1], tspan, [p0;p1])

sol0  = solve(prob0, Tsit5())
sol1  = solve(prob1, Tsit5())
sol   = solve(prob,  Tsit5())

plot(sol0)
plot!(sol1)
plot!(sol, linestyle=:dash, linewidth=2)

function condition(u,t,integrator)
	u[3]
end
function affect!(integrator)
  integrator.u[1] = integrator.u[1] + rand(Float64)
end
cb = ContinuousCallback(condition,affect!)
prob = remake(prob, tspan=(0.0,10000.0))
sol   = solve(prob, Tsit5(); callback=cb)
plot(sol, vars=(0,1))

max_stars_repo_name: cmarcotte/CTESN

max_issues_repo_name: cmarcotte/CTESN

max_forks_repo_name: cmarcotte/CTESN


text:
using MyPack
using Test

@testset "MyPack.jl" begin
    # Write your own tests here.
end


text:
function createClass(dict)
   f = open("BilSoz.cpp","w")
   print(f,"//BilSoz.cpp\n#include \"bilsoz/BilSoz.hpp\"\nunordered_map<string,string> BilSoz::dict = {\n")
   i=1
   for (i,(k,w)) in enumerate(dict)
      print(f,"{\"",k,"\",\"",w, "\"}")
      i != length(dict) ? println(f,",") : println(f)
      i+=1
   end
   println(f,"};")
   close(f)
end

dict = Dict(split.(readlines(open("../bilsoz.tsv")),"\t"))
createClass(dict)

max_stars_repo_name: ekurtulus/bilsoz

max_issues_repo_name: ekurtulus/bilsoz

max_forks_repo_name: ekurtulus/bilsoz


text:
using Distributions
using StatsBase
include("./baum_welch_lib.jl")
using Hmm

# 状態数
state_num = 2
# 出力シンボル(この場合は 0 と 1 )
symbol = [0 1]
# 出力シンボルの数
symbol_num = length(symbol)

# これが本当の値
# 遷移確率行列
A = [0.9 0.1; 0.15 0.85]
# 出力確率行列
B = [0.8 0.2; 0.2 0.8]
# 初期確率
ρ = [0.01 0.99]

# これはパラメータ学習の初期値
# 遷移確率行列
eA = [0.4 0.6; 0.3 0.7]
# 出力確率行列
eB = [0.8 0.2; 0.2 0.8]
# 初期確率
eρ = [1/2 1/2]

srand(1234)

function simulate(nSteps)

  observations = zeros(nSteps)
  states = zeros(nSteps)
  states[1] = sample(symbol, Weights(vec(ρ)))
  observations[1] = sample(symbol, Weights(B[Int(states[1])+1, :]))
  for t in 2:nSteps
    states[t] = sample(symbol, Weights(A[Int(states[t-1])+1, :]))
    observations[t] = sample(symbol, Weights(B[Int(states[t])+1, :]))
  end
  return observations,states
end

o1, s = simulate(500)
o2, s = simulate(20)
o3, s = simulate(9)
o4, s = simulate(3)
obs = [o1, o2, o3, o4]

hmm = Hmm.hmm_initialization(eA, eB, eρ)
Hmm.train(hmm, obs, 1e-4, 100000)

println("Actual parameter ρ")
println(ρ)
println("Estimated parameter ρ")
println(hmm.ρ)
println("Actual parameter A")
println(A)
println("Estimated parameter A")
println(hmm.A)
println("Actual parameter B")
println(B)
println("Estimated parameter B")
println(hmm.B)


text:
using Test
using TestSetExtensions
using MetaBalls:MetaBalls

include("generators.jl")

@testset "generate_voxels" begin
    x = generate_voxels(4)

        # printstyled(x[1], "\n"; color=:blue)
        # printstyled(size(x[1], 1), "\n"; color=:blue)

    @test size(x[1], 1) == 8
    @test size(x, 1) == 27
end


text:
export Rate, reset, sleep, cycleTime, expectedCycleTime
export Duration, Time, init, isSimTime, isSystemTime, isValid, toSec
export now, setNow, shutdown, sleepUntil, useSystemTime, waitForValid

# --- Rate ---

"""
    Rate(frequency::Number)

Constructor, creates a Rate.
"""
Rate(frequency::Number) = @cxxnew ros::Rate(frequency)

"""
    Rate(duration)

Constructor, creates a Rate from a Duration.
"""
Rate(duration::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},(false, false, false)}}) = @cxxnew ros::Rate(*(duration))

"""
    reset(rate)

Sets the start time for the rate to now.
"""
reset(rate::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Rate")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Rate")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Rate")},(false, false, false)},(false, false, false)}}) = @cxx rate->reset()

"""
    sleep(rate)

Sleeps for any leftover time in a cycle. Calculated from the last time sleep, reset, or the constructor was called.
"""
sleep(rate::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Rate")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Rate")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Rate")},(false, false, false)},(false, false, false)}}) = @cxx rate->sleep()

"""
    cycleTime(rate)

Get the actual run time of a cycle from start to sleep.
"""
cycleTime(rate::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Rate")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Rate")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Rate")},(false, false, false)},(false, false, false)}}) = @cxx rate->cycleTime()

"""
    expectedCycleTime(rate)

Get the expected cycle time – one over the frequency passed in to the constructor.
"""
expectedCycleTime(rate::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Rate")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Rate")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Rate")},(false, false, false)},(false, false, false)}}) = @cxx rate->expectedCycleTime()

# --- Duration ---
"""
    Duration()

Default constructor, creates a Duration.
"""
Duration() = @cxxnew ros::Duration()

"""
    Duration(t::Number)

Constructor, creates a Duration.
"""
Duration(t::Number) = @cxxnew ros::Duration(t)

"""
    Duration(rate)

Constructor, creates a Duration from a Rate.
"""
Duration(rate::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Rate")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Rate")},(false, false, false)},(false, false, false)}}) = @cxxnew ros::Duration(*(rate))

"""
    Duration(sec::Int, nsec::Int)

Constructor, creates a Duration.
"""
Duration(sec::Int, nsec::Int) = @cxxnew ros::Duration(sec, nsec)

"""
    sleep(duration)

Sleep for the amount of time specified by this Duration. If a signal interrupts the sleep, resleeps for the time remaining.
"""
sleep(duration::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},(false, false, false)}}) = @cxx duration->sleep()

"""
    toSec(duration)

Turn Duration into floating point seconds.
"""
toSec(duration::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},(false, false, false)}}) = @cxx duration->toSec()

# --- Time ---
"""
    Time()

Default constructor, creates Time.
"""
Time() = @cxxnew ros::Time()

"""
    Time(t::Number)

Constructor, creates Time.
"""
Time(t::Number) = @cxxnew ros::Time(t)

"""
    Time(sec::Int, ncsec::Int)

Constructor, creates Time.
"""
Time(sec::Int, nsec::Int) = @cxxnew ros::Time(sec, nsec)

"""
    init(time)

Initialise Time.
"""
init(t::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},(false, false, false)}}) = @cxx t->init()

"""
    isSimTime(time)

Returns whether or not sim time is used.
"""
isSimTime(t::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},(false, false, false)}}) = @cxx t->isSimTime()

"""
    isSystemTime(time)

Returns whether or not system time is used.
"""
isSystemTime(t::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},(false, false, false)}}) = @cxx t->isSystemTime()

"""
    isValid(time)

Returns whether or not the current time source is valid. Simulation time is valid if it is non-zero.
"""
isValid(t::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},(false, false, false)}}) = @cxx t->isValid()

"""
    now(time)

Retrieve the current time. If ROS clock time is in use, this returns the time according to the ROS clock. Otherwise returns the current wall clock time.
"""
now(t::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},(false, false, false)}}) = @cxx t->now()

"""
    setNow(time, now) = @cxx t->setNow(now)

Set now to `now` time.
"""
setNow(t::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},(false, false, false)}}, now::Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N) = @cxx t->setNow(now)

setNow(t::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},(false, false, false)}}, now::Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}) = @cxx t->setNow(*(now))

setNow(t::Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N, now::Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N) = @cxx t->setNow(now)

setNow(t::Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N, now::Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}) = @cxx t->setNow(*(now))

"""
    shutdown(time)

Shuts down the Time.
"""
shutdown(t::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},(false, false, false)}}) = @cxx t->shutdown()

"""
    sleepUntil(time, _end)

Sleep until a specific time has been reached.
"""
sleepUntil(t::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},(false, false, false)}}, _end::Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N) = @cxx t->sleepUntil(_end)

sleepUntil(t::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},(false, false, false)}}, _end::Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}) = @cxx t->sleepUntil(*(_end))

sleepUntil(t::Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N, _end::Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N) = @cxx t->sleepUntil(_end)

sleepUntil(t::Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N, _end::Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}) = @cxx t->sleepUntil(*(_end))

"""
    useSystemTime(time)

Switch to system time.
"""
useSystemTime(t::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},(false, false, false)}}) = @cxx t->useSystemTime()

"""
    waitForValid(time)

Wait for time source to become valid.
"""
waitForValid(t::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},(false, false, false)}}) = @cxx t->waitForValid()

"""
    toSec(time)

Turn Time into floating point seconds.
"""
toSec(t::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},(false, false, false)}}) = @cxx t->toSec()

function Base.getproperty(c::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},(false, false, false)}, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},(false, false, false)}}, f::Symbol)
    if f == :sec
        @cxx c->sec
    elseif f == :nsec
        @cxx c->nsec
    else
        throw(ErrorException("Time/Duration does not have a field of $(f)!"))
    end
end

function Base.setproperty!(c::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},(false, false, false)}, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},(false, false, false)}}, f::Symbol, value::Number)
    if f == :sec
        if typeof(c) <: Cxx.CxxCore.CppValue
            icxx"$c.sec = $value;"
        else
            icxx"$c->sec = $value;"
        end
    elseif f == :nsec
        if typeof(c) <: Cxx.CxxCore.CppValue
            icxx"$c.nsec = $value;"
        else
            icxx"$c->nsec = $value;"
        end
    else
        throw(ErrorException("Time/Duration does not have a field of $(f)!"))
    end
end

function Base.:+(dt1::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},(false, false, false)}, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},(false, false, false)}}, dt2::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},(false, false, false)}, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},(false, false, false)}})
    if typeof(dt1) <: Cxx.CxxCore.CppPtr
        if typeof(dt2) <: Cxx.CxxCore.CppPtr
            icxx"return *$dt1 + *$dt2;"
        else
            icxx"return *$dt1 + $dt2;"
        end
    else
        if typeof(dt2) <: Cxx.CxxCore.CppPtr
            icxx"return $dt1 + *$dt2;"
        else
            icxx"return $dt1 + $dt2;"
        end
    end
end

function Base.:-(dt1::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},(false, false, false)}, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},(false, false, false)}}, dt2::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Time")},(false, false, false)},(false, false, false)}, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},N} where N, Cxx.CxxCore.CppPtr{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Duration")},(false, false, false)},(false, false, false)}})
    if typeof(dt1) <: Cxx.CxxCore.CppPtr
        if typeof(dt2) <: Cxx.CxxCore.CppPtr
            icxx"return *$dt1 - *$dt2;"
        else
            icxx"return *$dt1 - $dt2;"
        end
    else
        if typeof(dt2) <: Cxx.CxxCore.CppPtr
            icxx"return $dt1 - *$dt2;"
        else
            icxx"return $dt1 - $dt2;"
        end
    end
end


text:
function test_mobil_creation()
	println("\t\tTesting MOBIL Creation")
	p1 = MOBILParam()
	p2 = MOBILParam(p=0.5)
	p3 = MOBILParam(b_safe=5.,a_thr=0.3)
	ps = [MOBILParam(s) for s in ["cautious";"normal";"aggressive"]]
end

function test_mobil_equality()
	#println("\t\tTesting MOBIL Equality")
	p1 = MOBILParam()
	p2 = MOBILParam(p=0.5)
	p3 = MOBILParam(b_safe=5.,a_thr=0.3)
	ps = [MOBILParam(s) for s in ["cautious";"normal";"aggressive"]]
	test_equality("MOBIL",ps)
	assert(p1 == p1)
	assert(p1 != p2)
	assert(p2 != p1)
	assert(p1 == MOBILParam(0.25,4.,0.2))
	assert(MOBILParam(0.25,4.,0.2) == p1)
end

function test_mobil_hashing()
	println("\t\tTesting MOBIL Hashing")
	p1 = MOBILParam()
	p2 = MOBILParam(p=0.5)
	p3 = MOBILParam(b_safe=5.,a_thr=0.3)
	ps = [MOBILParam(s) for s in ["cautious";"normal";"aggressive"]]
	test_hashing("MOBIL",ps)
end

function test_get_neighborhood()
	println("\t\tTesting get_neighborhood")
	##TODO: assertions
	nb_lanes = 3
	pp = PhysicalParam(nb_lanes,lane_length=12.)
	bs = IDMMOBILBehavior[IDMMOBILBehavior(x[1],x[2],x[3],idx) for (idx,x) in enumerate(product(["cautious","normal","aggressive"],[pp.v_slow;pp.v_med;pp.v_fast],[pp.l_car]))]
	#CASE: just agent car
	cs = CarState[CarState(6.,2,31.,0,bs[1],1)]
	dp1 = 6.-pp.l_car
	dp2 = 6.-pp.l_car
	_a = MLAction(0,0)

	#just make sure it doesn't explode--this should be handled by get_mobil_lane_change
	dmodel = IDMMOBILModel(length(cs), pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,1)
	assert(length(nbhd) == 6)
	for jdx in nbhd
		assert(jdx == 0)
	end
	#CASE: nobody to the left
	cs = CarState[CarState(12.,2,31.,0,bs[1],2),CarState(6.,2,31.,0,bs[1],3),CarState(0.,2,31.,0,bs[1],4),CarState(12.,1,31.,0,bs[1],5),CarState(0.,1,31.,0,bs[1],6)]

	dmodel = IDMMOBILModel(length(cs), pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,2)
	assert(length(nbhd) == 6)
	for jdx in nbhd[[3;6]]
		assert(jdx == 0)
	end
	assert(nbhd[1] , 4)
	assert(nbhd[2] , 1)
	assert(nbhd[4] , 5)
	assert(nbhd[5] , 3)
	#CASE: nobody to the right
	cs = CarState[CarState(12.,2,31.,0,bs[1],7),CarState(0.,2,31.,0,bs[1],8),CarState(6.,2,31.,0,bs[1],9),CarState(12.,3,31.,0,bs[1],10),CarState(0.,3,31.,0,bs[1],11)]

	dmodel = IDMMOBILModel(length(cs), pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,3)
	assert(length(nbhd) == 6)
	for jdx in nbhd[[1;4]]
		assert(jdx , 0)
	end
	assert(nbhd[3] , 4)
	assert(nbhd[2] , 1)
	assert(nbhd[6] , 5)
	assert(nbhd[5] , 2)
	#CASE: no one ahead
	cs = CarState[CarState(0.,3,31.,0,bs[1],12),CarState(0.,2,31.,0,bs[1],13),CarState(0.,1,31.,0,bs[1],14),CarState(6.,2,31.,0,bs[1],15)]

	dmodel = IDMMOBILModel(length(cs), pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,4)
	assert(length(nbhd) , 6)
	for jdx in nbhd[[1;2;3]]
		assert(jdx , 0)
	end
	assert(nbhd[4] , 3)
	assert(nbhd[5] , 2)
	assert(nbhd[6] , 1)
	#CASE: no one behind
	cs = CarState[CarState(12.,3,31.,0,bs[1],16),CarState(12.,2,31.,0,bs[1],17),CarState(12.,1,31.,0,bs[1],18),CarState(6.,2,31.,0,bs[1],19)]

	dmodel = IDMMOBILModel(length(cs), pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,4)
	assert(length(nbhd) , 6)
	for jdx in nbhd[[4;5;6]]
		assert(jdx , 0)
	end
	assert(nbhd[1] , 3)
	assert(nbhd[2] , 2)
	assert(nbhd[3] , 1)
	#CASE: full house
	cs = CarState[CarState(12.,2,31.,0,bs[1],20),CarState(0.,2,31.,0,bs[1],21),CarState(12.,1,31.,0,bs[1],22),CarState(0.,1,31.,0,bs[1],23),CarState(12.,3,35.,0,bs[1],24),CarState(0.,3,35.,0,bs[1],25),CarState(6.,2,31.,0,bs[1],26)]

	dmodel = IDMMOBILModel(length(cs), pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,7)
	assert(length(nbhd) , 6)
	for jdx in nbhd
		assert(jdx != 0)
	end
	assert(nbhd[1] , 3)
	assert(nbhd[2] , 1)
	assert(nbhd[3] , 5)
	assert(nbhd[4] , 4)
	assert(nbhd[5] , 2)
	assert(nbhd[6] , 6)
end
#TODO test_get_dv_ds(); FIX VVV

function test_is_lanechange_dangerous()
	println("\t\tTesting is_lanechange_dangerous")

	nb_lanes = 3
	_a = MLAction(0,0)
	pp = PhysicalParam(nb_lanes,lane_length=12.)
	bs = IDMMOBILBehavior[IDMMOBILBehavior(x[1],x[2],x[3],idx) for (idx,x) in enumerate(product(["cautious","normal","aggressive"],[pp.v_slow;pp.v_med;pp.v_fast],[pp.l_car]))]
	#CASE: no one, in the way
	cs = CarState[CarState(6.,2,35.,0,IDMMOBILBehavior("aggressive",35.,4.,1),27)]
	dmodel = IDMMOBILModel(length(cs),pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,1)
	assert(is_lanechange_dangerous(pp,s,nbhd,1,0),false)
	assert(is_lanechange_dangerous(pp,s,nbhd,1,1),false)
	assert(is_lanechange_dangerous(pp,s,nbhd,1,-1),false)

	#CASE: person close by, same velocity
	cs = CarState[CarState(6.,2,35.,0,IDMMOBILBehavior("aggressive",35.,4.,1),28),CarState(12.,3,35.,0,IDMMOBILBehavior("aggressive",35.,4.,1),29)]
	dmodel = IDMMOBILModel(length(cs),pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,1)
	assert(is_lanechange_dangerous(pp,s,nbhd,1,1),false)

	#slightly behind
	cs = CarState[CarState(6.,2,35.,0,IDMMOBILBehavior("aggressive",35.,4.,1),30),CarState(0.,3,35.,0,IDMMOBILBehavior("aggressive",35.,4.,1),31)]
	dmodel = IDMMOBILModel(length(cs),pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,1)
	assert(is_lanechange_dangerous(pp,s,nbhd,1,1),false)
	#CASE: person a little behind, slower
	cs = CarState[CarState(6.,2,35.,0,IDMMOBILBehavior("aggressive",35.,4.,1),32),CarState(1.,3,32.,0,IDMMOBILBehavior("aggressive",35.,4.,1),33)]
	dmodel = IDMMOBILModel(length(cs),pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,1)
	assert(is_lanechange_dangerous(pp,s,nbhd,1,1),false)
	#CASE: person a little ahead, faster
	cs = CarState[CarState(6.,2,35.,0,IDMMOBILBehavior("aggressive",35.,4.,1),34),CarState(1.,3,32.,0,IDMMOBILBehavior("aggressive",35.,4.,1),35)]
	dmodel = IDMMOBILModel(length(cs),pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,2)
	assert(is_lanechange_dangerous(pp,s,nbhd,2,-1),false)
	#CASE: person a little behind, faster
	cs = CarState[CarState(6.,2,32.,0,IDMMOBILBehavior("aggressive",35.,4.,1),36),CarState(1.,3,35.,0,IDMMOBILBehavior("aggressive",35.,4.,1),37)]
	dmodel = IDMMOBILModel(length(cs),pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,1)
	assert(is_lanechange_dangerous(pp,s,nbhd,1,1),true)
	#CASE: person a little ahead, slower
	cs = CarState[CarState(6.,2,32.,0,IDMMOBILBehavior("aggressive",35.,4.,1),38),CarState(1.,3,35.,0,IDMMOBILBehavior("aggressive",35.,4.,1),39)]
	dmodel = IDMMOBILModel(length(cs),pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,2)
	assert(is_lanechange_dangerous(pp,s,nbhd,2,-1),true)
	#CASE: next to each other
	cs = CarState[CarState(6.,2,32.,0,IDMMOBILBehavior("aggressive",35.,4.,1),40),CarState(6.,3,33.,0,IDMMOBILBehavior("aggressive",35.,4.,1),41)]
	dmodel = IDMMOBILModel(length(cs),pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,2)
	assert(is_lanechange_dangerous(pp,s,nbhd,2,-1),true)
end

function test_get_rear_accel()
	println("\t\tTesting get_rear_accel")
	nb_lanes = 3
	_a = MLAction(0,0)
	pp = PhysicalParam(nb_lanes,lane_length=12.)
	bs = IDMMOBILBehavior[IDMMOBILBehavior(x[1],x[2],x[3],idx) for (idx,x) in enumerate(product(["cautious","normal","aggressive"],[pp.v_slow;pp.v_med;pp.v_fast],[pp.l_car]))]
	#No one in that position
	cs = CarState[CarState(6.,2,35.,0,IDMMOBILBehavior("aggressive",35.,4.,1),42)]
	dmodel = IDMMOBILModel(length(cs),pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,1)
	a,a_ = get_rear_accel(pp,s,nbhd,1,0)
	assert(a,0.)
	assert(a_,0.)
	a,a_ = get_rear_accel(pp,s,nbhd,1,1)
	assert(a,0.)
	assert(a_,0.)
	a,a_ = get_rear_accel(pp,s,nbhd,1,-1)
	assert(a,0.)
	assert(a_,0.)
	#will go to free flow
	cs = CarState[CarState(6.,2,35.,0,IDMMOBILBehavior("aggressive",35.,4.,1),42),
								CarState(0.,3,35.,0,IDMMOBILBehavior("aggressive",35.,4.,1),43),
								CarState(0.,2,31.,0,IDMMOBILBehavior("aggressive",31.,4.,1),44),
								CarState(0.,1,27.,0,IDMMOBILBehavior("aggressive",27.,4.,1),45)]
	dmodel = IDMMOBILModel(length(cs),pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,1)
	a,a_ = get_rear_accel(pp,s,nbhd,1,0)
	assert(a_,a,>=)
	a,a_ = get_rear_accel(pp,s,nbhd,1,1)
	assert(a_,a,>=)
	a,a_ = get_rear_accel(pp,s,nbhd,1,-1)
	assert(a_,a,>=)
	#some synthetic case
end

function test_get_mobil_lane_change()
	println("\t\tTesting get_mobil_lane_change")
	nb_lanes = 2
	_a = MLAction(0,0)
	pp = PhysicalParam(nb_lanes,lane_length=12.)
	bs = IDMMOBILBehavior[IDMMOBILBehavior(x[1],x[2],x[3],idx) for (idx,x) in enumerate(product(["cautious","normal","aggressive"],[pp.v_slow;pp.v_med;pp.v_fast],[pp.l_car]))]
	#CASE: it's faster, but there's no space--is this even a real case?
	#CASE: it's faster and there is space
	cs = CarState[CarState(12.,2,35.,0,IDMMOBILBehavior("aggressive",35.,4.,1),1),
								CarState(11.,1,27.,0,IDMMOBILBehavior("cautious",27.,4.,1),2),
								CarState(6.,1,31.,0,IDMMOBILBehavior("aggressive",35.,4.,1),3)]
	dmodel = IDMMOBILModel(length(cs),pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,3)
	assert(get_mobil_lane_change(s.cars[3].behavior, pp,s,nbhd,3),1)
	#CASE: it's slower and there is space
	cs = CarState[CarState(12.,2,31.,0,IDMMOBILBehavior("cautious",31.,4.,1),6),CarState(11.,1,35.,0,IDMMOBILBehavior("cautious",35.,4.,1),5),CarState(6.,1,35.,0,IDMMOBILBehavior("aggressive",35.,4.,1),4)]
	dmodel = IDMMOBILModel(length(cs),pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,3)
	assert(get_mobil_lane_change(s.cars[3].behavior, pp,s,nbhd,3),0)
	#CASE: someone is going fast behind me and i'm slow
	cs = CarState[CarState(0.,1,35.,0,IDMMOBILBehavior("aggressive",35.,4.,1),7),CarState(6.,1,27.,0,IDMMOBILBehavior("cautious",27.,4.,1),8)]
	dmodel = IDMMOBILModel(length(cs),pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,2)
	assert(get_mobil_lane_change(s.cars[2].behavior,pp,s,nbhd,2),1)
	###repeat for other side
	#CASE: it's faster and there is space
	cs = CarState[CarState(12.,1,35.,0,IDMMOBILBehavior("aggressive",35.,4.,1),9),CarState(11.,2,27.,0,IDMMOBILBehavior("cautious",27.,4.,1),10),CarState(6.,2,31.,0,IDMMOBILBehavior("aggressive",35.,4.,1),11)]
	dmodel = IDMMOBILModel(length(cs),pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,3)
	assert(get_mobil_lane_change(s.cars[3].behavior,pp,s,nbhd,3),-1)
	#CASE: it's slower and there is space
	cs = CarState[CarState(12.,1,31.,0,IDMMOBILBehavior("cautious",31.,4.,1),0),CarState(11.,2,35.,0,IDMMOBILBehavior("cautious",35.,4.,1),0),CarState(6.,2,35.,0,IDMMOBILBehavior("aggressive",35.,4.,1),0)]
	dmodel = IDMMOBILModel(length(cs),pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,3)
	assert(get_mobil_lane_change(s.cars[3].behavior, pp,s,nbhd,3),0)
	#CASE: someone is going fast behind me and i'm slow
	cs = CarState[CarState(0.,2,35.,0,IDMMOBILBehavior("aggressive",35.,4.,1),0),CarState(6.,2,27.,0,IDMMOBILBehavior("cautious",27.,4.,1),0)]
	dmodel = IDMMOBILModel(length(cs),pp)
	s = MLState(0.0,0.0,cs)
	nbhd = get_neighborhood(pp,s,2)
	assert(get_mobil_lane_change(s.cars[2].behavior,pp,s,nbhd,2),-1)
end

function failure_2()
    println("\t\tTesting specific failure: car runs into ego")
    nb_lanes = 4
    pp = PhysicalParam(nb_lanes,lane_length=100.) #2.=>col_length=8
    _discount = 1.
    nb_cars=10
    rmodel = NoCrashRewardModel()
    dmodel = NoCrashIDMMOBILModel(nb_cars, pp)
    mdp = NoCrashMDP{NoCrashRewardModel, typeof(dmodel.behaviors)}(dmodel, rmodel, _discount, true);
    rng = MersenneTwister(2)
    s = Multilane.MLState(0.0,0.0,[Multilane.CarState(50.0,1.0,27.0,0.0,Multilane.NORMAL,0),Multilane.CarState(45.66875308699299,2.0,29.042803490064543,-0.6666666666666666,Multilane.IDMMOBILBehavior(Multilane.IDMParam(2.0,2.0,0.8,35.0,4.0,4.0),Multilane.MOBILParam(0.0,4.0,0.2),9),0),Multilane.CarState(85.41074952313767,3.0,29.831781823255017,0.0,Multilane.IDMMOBILBehavior(Multilane.IDMParam(2.0,2.0,0.8,31.0,4.0,4.0),Multilane.MOBILParam(0.0,4.0,0.2),6),0),Multilane.CarState(13.78112065279241,3.0,27.0,0.0,Multilane.IDMMOBILBehavior(Multilane.IDMParam(1.5,1.5,1.4,35.0,4.0,4.0),Multilane.MOBILParam(0.25,4.0,0.2),8),0),Multilane.CarState(5.897872040018292,1.0,27.39955827605183,0.0,Multilane.IDMMOBILBehavior(Multilane.IDMParam(1.5,1.5,1.4,35.0,4.0,4.0),Multilane.MOBILParam(0.25,4.0,0.2),8),0)])
    a = MLAction(0.0,0.0)
    sp, r = generate_sr(mdp, s, a, rng)
    @test !is_crash(mdp, s, sp)
end

#=
function failure_3()
    println("\t\tTesting specific failure: cars lane change into one another")
    mdp = Multilane.MLMDP{Multilane.MLState,Multilane.MLAction,Multilane.NoCrashIDMMOBILModel,Multilane.NoCrashRewardModel}(Multilane.NoCrashIDMMOBILModel(10,Multilane.PhysicalParam(0.75,2.0,4.0,31.0,4.0,35.0,27.0,31.0,35.0,27.0,20.0,4,100.0),Multilane.BehaviorModel[Multilane.IDMMOBILBehavior(Multilane.IDMParam(1.0,1.0,2.0,27.5,4.0,4.0),Multilane.MOBILParam(0.5,4.0,0.2),1),Multilane.IDMMOBILBehavior(Multilane.IDMParam(1.5,1.5,1.4,27.5,4.0,4.0),Multilane.MOBILParam(0.25,4.0,0.2),2),Multilane.IDMMOBILBehavior(Multilane.IDMParam(2.0,2.0,0.8,27.5,4.0,4.0),Multilane.MOBILParam(0.0,4.0,0.2),3),Multilane.IDMMOBILBehavior(Multilane.IDMParam(1.0,1.0,2.0,31.0,4.0,4.0),Multilane.MOBILParam(0.5,4.0,0.2),4),Multilane.IDMMOBILBehavior(Multilane.IDMParam(1.5,1.5,1.4,31.0,4.0,4.0),Multilane.MOBILParam(0.25,4.0,0.2),5),Multilane.IDMMOBILBehavior(Multilane.IDMParam(2.0,2.0,0.8,31.0,4.0,4.0),Multilane.MOBILParam(0.0,4.0,0.2),6),Multilane.IDMMOBILBehavior(Multilane.IDMParam(1.0,1.0,2.0,35.0,4.0,4.0),Multilane.MOBILParam(0.5,4.0,0.2),7),Multilane.IDMMOBILBehavior(Multilane.IDMParam(1.5,1.5,1.4,35.0,4.0,4.0),Multilane.MOBILParam(0.25,4.0,0.2),8),Multilane.IDMMOBILBehavior(Multilane.IDMParam(2.0,2.0,0.8,35.0,4.0,4.0),Multilane.MOBILParam(0.0,4.0,0.2),9)],StatsBase.WeightVec{Float64,Array{Float64,1}}([1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0],9.0),1.0,1.3333333333333333,0.5,20.0,0.5,[1.0,1.0,1.0,1.0],2.0,false),Multilane.NoCrashRewardModel(100.0,10.0,1.5,1),1.0)
    s = Multilane.MLState(false,0.0,0.0,[Multilane.CarState(50.0,2.5,30.75,0.0,Multilane.NORMAL,1),Multilane.CarState(99.56731292388635,3.0,34.36785842369873,0.0,Multilane.IDMMOBILBehavior(Multilane.IDMParam(2.0,2.0,0.8,35.0,4.0,4.0),Multilane.MOBILParam(0.0,4.0,0.2),9),4),Multilane.CarState(88.9857138446193,4.0,29.878676022074934,0.0,Multilane.IDMMOBILBehavior(Multilane.IDMParam(1.5,1.5,1.4,31.0,4.0,4.0),Multilane.MOBILParam(0.25,4.0,0.2),5),12),Multilane.CarState(83.55003341854245,1.0,27.45737686175689,0.0,Multilane.IDMMOBILBehavior(Multilane.IDMParam(1.5,1.5,1.4,27.5,4.0,4.0),Multilane.MOBILParam(0.25,4.0,0.2),2),13),Multilane.CarState(5.833404602969644,1.5,31.164522812071382,0.6666666666666666,Multilane.IDMMOBILBehavior(Multilane.IDMParam(2.0,2.0,0.8,35.0,4.0,4.0),Multilane.MOBILParam(0.0,4.0,0.2),9),14),Multilane.CarState(6.364105400258811,3.0,30.995876598243825,0.0,Multilane.IDMMOBILBehavior(Multilane.IDMParam(2.0,2.0,0.8,31.0,4.0,4.0),Multilane.MOBILParam(0.0,4.0,0.2),6),15),Multila
    sp, r = generate_sr(mdp, s, MLAction(0,0), MersenneTwister(5))
    @test !Multilane.occupation_overlap(sp.cars[4].y, sp.cars[5].y)
end
=#

function test_mobil()
	println("\tTesting MOBIL Units...")
	test_mobil_creation()
	test_mobil_equality()
	test_mobil_hashing()
	test_get_neighborhood()
	test_is_lanechange_dangerous()
	test_get_rear_accel()
	test_get_mobil_lane_change()
    failure_2()
    # failure_3() # removed because of behavior generator change
end


text:
using Distances
using LinearAlgebra
include("../src/kernels/kernel.jl")

#=
NOT important
compare two ways of compute the kernel matrix
compare my KernelMat function with "Distances.jl"

Result:
    x, y = 100 * rand(1000, 3)
    theta = [rand(), 2 * rand()]
    dist1:
        346.388 ms (6012 allocations: 137.65 MiB)
    dist2:
        189.920 ms (3000002 allocations: 328.06 MiB)

Conclusion:
    seems like my implementation is faster?
    so keep it for now.
 =#

function dist1(x, y, theta)
    r = pairwise(Euclidean(), x', y', dims=2)
    kernelfun = r -> Kernel_SE(r, theta)
    K = kernelfun.(r)
    return K
end

function dist2(x, y, theta)
    K = KernelMat(x, y, Kernel_SE, theta)
    return K
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    domain(object)

Return underlying domain of the `object`.
"""
function domain end

"""
    values(object)

Return the values of `object` as a table.
"""
function values end


text:
export solveMG,solveGMRES_MG,solveBiCGSTAB_MG,solveCG_MG,getAfun

function solveMG(param::MGparam{VAL,IND},b::Array{VAL},x::Array{VAL},verbose::Bool) where {VAL,IND}
#MGType = getMGType(param,b);

param = adjustMemoryForNumRHS(param,size(b,2));
tol = param.relativeTol;
numCores = param.numCores;
oneType = one(VAL);
zeroType = zero(VAL);
maxIter = param.maxOuterIter;
A = param.As[1];
r = param.memCycle[1].r;
r[:] = b;

if norm(x)==0
    res = norm(b);
    res_init = res;
else
	SpMatMul(-oneType,A,x,oneType,r,numCores)#  r -= A'*x;
    res = norm(r);
    res_init = res;
end

for count = 1:maxIter
	x = recursiveCycle(param,b,x,1);
	SpMatMul(-oneType,A,x,zeroType,r,numCores); #  r = -A'*x;
	addVectors(oneType,b,r); # r = r + b;
	res_prev = res;
	res = norm(r);
	if verbose
		println(string("Cycle ",count," done with relres: ",res/res_init,". Convergence factor: ",res/res_prev));
	end
	if res/res_init < tol
		break;
	end
end
return x,param;
end
####################################################################################################################
# this function checks if the memory allocated in param fits the nrhs and generates a function for applying the cycle.
# Also - it checks if the hierarchy needs to be transposed or not.
function getMultigridPreconditioner(param::MGparam{VAL,IND},B::Array,verbose::Bool=false) where {VAL,IND}
	n = size(B,1)
	nrhs = size(B,2);
	if hierarchyExists(param)==false
		println("You have to do a setup first.")
	end
	param = adjustMemoryForNumRHS(param,nrhs);
	z = param.memCycle[1].x;

	mixed_precision = VAL!=eltype(B);
	MMG = identity;
	if mixed_precision
		bl =  param.memCycle[1].b;
		z2 = zeros(eltype(B),size(B));
		MMG = (b) -> (z[:] .= 0.0;bl[:] .= b; recursiveCycle(param,bl,z,1); z2[:] .= z; return z2;);
	else
		MMG = (b) -> (z[:] .= 0.0;recursiveCycle(param,b,z,1); return z);
	end

	return MMG;
end

function getAfun(AT::SparseMatrixCSC{VAL,IND},Az::Array{VAL},numCores::Int64) where {VAL,IND}
	function Afun(z::Array{VAL})
		SpMatMul(AT,z,Az,numCores);
		return Az;
	end
	return Afun;
end

function solveBiCGSTAB_MG(AT::SparseMatrixCSC,param::MGparam{VAL,IND},b::Array,x0::Array,verbose::Bool = false) where {VAL,IND}
	return solveBiCGSTAB_MG(getAfun(AT,zeros(eltype(b),size(b)),param.numCores),param,b,x0,verbose);
end
function solveCG_MG(AT::SparseMatrixCSC,param::MGparam{VAL,IND},b::Array,x0::Array,verbose::Bool = false) where {VAL,IND}
	return solveCG_MG(getAfun(AT,zeros(eltype(b),size(b)),param.numCores),param,b,x0,verbose);
end
function solveGMRES_MG(AT::SparseMatrixCSC,param::MGparam{VAL,IND},b::Array,x0::Array,flexible::Bool,inner::Int64,verbose::Bool = false) where {VAL,IND}
	return solveGMRES_MG(getAfun(AT,zeros(eltype(b),size(b)),param.numCores),param,b,x0,flexible,inner,verbose);
end

function solveBiCGSTAB_MG(Afun::Function,param::MGparam{VAL,IND},b::Array,x0::Array,verbose::Bool = false) where {VAL,IND}
MMG = getMultigridPreconditioner(param,b,verbose);
out= -2;
if verbose
	out = 1;
end
if size(b,2)==1
	b = vec(b);
	x, flag,rnorm,iter = KrylovMethods.bicgstb(Afun,b,tol = param.relativeTol,maxIter = param.maxOuterIter,M1 = MMG,M2 = identity, x = x0,out=out);
else
	x, flag,rnorm,iter = KrylovMethods.blockBiCGSTB(Afun,b,tol = param.relativeTol,maxIter = param.maxOuterIter,M1 = MMG,M2 = identity, x = x0,out=out);
end
nprec = 2*iter*size(b,2) + (flag==-3)*size(b,2);
return x,param,iter,nprec
end

function solveCG_MG(Afun::Function,param::MGparam{VAL,IND},b::Array,x0::Array,verbose::Bool = false) where {VAL,IND}
MMG = getMultigridPreconditioner(param,b,verbose);
out = -2;
if verbose
	out = 1;
end
if size(b,2)==1
	b = vec(b);
	x, flag,rnorm,iter = KrylovMethods.cg(Afun,b,tol = param.relativeTol,maxIter = param.maxOuterIter,M = MMG, x = x0,out=out);
else
	x, flag,rnorm,iter = KrylovMethods.blockCG(Afun,b,tol = param.relativeTol,maxIter = param.maxOuterIter,M = MMG, X = x0,out=out);
end
return x,param,iter;
end

function solveGMRES_MG(Afun::Function,param::MGparam{VAL,IND},b::Array,x0::Array,flexible::Bool,inner::Int64,verbose::Bool = false) where {VAL,IND}
MMG = getMultigridPreconditioner(param,b,verbose);
out = -2;
if verbose
	out = 1;
end
if size(b,2)==1
	b = vec(b);
	x, flag,rnorm,iter,resvec = KrylovMethods.fgmres(Afun,b,inner,tol = param.relativeTol,maxIter = param.maxOuterIter,M = MMG, x = x0,out=out,flexible=flexible);
else
	x, flag,rnorm,iter,resvec = KrylovMethods.blockFGMRES(Afun,b,inner,tol = param.relativeTol,maxIter = param.maxOuterIter,M = MMG, X = x0,out=out,flexible=flexible);
end
return x,param,iter,resvec;
end