text:
# ----------------------------------------- #
# SyncPlot -- sync Plot object with display #
# ----------------------------------------- #
mutable struct SyncPlot
    plot::PlotlyBase.Plot
    scope::Scope
    options::Dict
    window::Union{Nothing,Blink.Window}
end

Base.getindex(p::SyncPlot, key) = p.scope[key] # look up Observables

WebIO.render(p::SyncPlot) = WebIO.render(p.scope)
Base.show(io::IO, mm::MIME"text/html", p::SyncPlot) = show(io, mm, p.scope)
Base.show(io::IO, mm::MIME"application/juno+plotpane", p::SyncPlot) = show(io, mm, p.scope)

function SyncPlot(
        p::Plot;
        options::AbstractDict=Dict("showLink"=> false),
        kwargs...
    )
    lowered = JSON.lower(p)
    id = string("#plot-", p.divid)

    # setup scope
    deps = [
        "Plotly" => joinpath(@__DIR__, "..", "assets", "plotly-latest.min.js"),
        joinpath(@__DIR__, "..", "assets", "plotly_webio.bundle.js")
    ]
    scope = Scope(imports=deps)
    scope.dom = dom"div"(id=string("plot-", p.divid))

    # INPUT: Observables for plot events
    scope["hover"] = Observable(Dict())
    scope["selected"] = Observable(Dict())
    scope["click"] = Observable(Dict())
    scope["relayout"] = Observable(Dict())
    scope["image"] = Observable("")
    scope["__gd_contents"] = Observable{Any}(Dict())  # for testing

    # OUTPUT: setup an observable which sends modify commands
    scope["_commands"] = Observable{Any}([])
    scope["_toImage"] = Observable(Dict())
    scope["_downloadImage"] = Observable(Dict())
    scope["__get_gd_contents"] = Observable("")

    onjs(scope["_toImage"], @js function (options)
        this.Plotly.toImage(this.plotElem, options).then(function(data)
            $(scope["image"])[] = data
        end)
    end)

    onjs(scope["__get_gd_contents"], @js function (prop)
        if prop == "data"
            $(scope["__gd_contents"])[] = this.plotElem.data
        end

        if prop == "layout"
            $(scope["__gd_contents"])[] = this.plotElem.layout
        end
    end)

    onjs(scope["_downloadImage"], @js function (options)
        this.Plotly.downloadImage(this.plotElem, options)
    end)

    # Do the respective action when _commands is triggered
    onjs(scope["_commands"], @js function (args)
        @var fn = args.shift()
        @var elem = this.plotElem
        @var Plotly = this.Plotly
        args.unshift(elem) # use div as first argument
        Plotly[fn].apply(this, args)
    end)

    onimport(scope, JSExpr.@js function (Plotly)

        # set up container element
        @var gd = this.dom.querySelector($id);

        # save some vars for later
        this.plotElem = gd
        this.Plotly = Plotly
        if (window.Blink !== undefined)
            # set css style for auto-resize
            gd.style.width = "100%";
            gd.style.height = "100vh";
            gd.style.marginLeft = "0%";
            gd.style.marginTop = "0vh";
        end

        window.onresize = () -> Plotly.Plots.resize(gd)

        # Draw plot in container
        Plotly.newPlot(
            gd, $(lowered[:data]), $(lowered[:layout]), $(options)
        )

        # hook into plotly events
        gd.on("plotly_hover", function (data)
            @var filtered_data = WebIO.CommandSets.Plotly.filterEventData(gd, data, "hover");
            if !(filtered_data.isnil)
                $(scope["hover"])[] = filtered_data.out
            end
        end)

        gd.on("plotly_unhover", () -> $(scope["hover"])[] = Dict())

        gd.on("plotly_selected", function (data)
            @var filtered_data = WebIO.CommandSets.Plotly.filterEventData(gd, data, "selected");
            if !(filtered_data.isnil)
                $(scope["selected"])[] = filtered_data.out
            end
        end)

        gd.on("plotly_deselect", () -> $(scope["selected"])[] = Dict())

        gd.on("plotly_relayout", function (data)
            @var filtered_data = WebIO.CommandSets.Plotly.filterEventData(gd, data, "relayout");
            if !(filtered_data.isnil)
                $(scope["relayout"])[] = filtered_data.out
            end
        end)

        gd.on("plotly_click", function (data)
            @var filtered_data = WebIO.CommandSets.Plotly.filterEventData(gd, data, "click");
            if !(filtered_data.isnil)
                $(scope["click"])[] = filtered_data.out
            end
        end)
    end)

    # create no-op `on` callback for image so it is _always_ sent
    # to us
    on(scope["image"]) do x end

    SyncPlot(p, scope, options, nothing)
end

function plot(args...; options=Dict(), kwargs...)
    SyncPlot(Plot(args...; kwargs...); options=options)
end

# Add some basic Julia API methods on SyncPlot that just forward onto the Plot
Base.size(sp::SyncPlot) = size(sp.plot)
Base.copy(sp::SyncPlot) = SyncPlot(copy(sp.plot), options=copy(sp.options))

Base.display(::PlotlyJSDisplay, p::SyncPlot) = display_blink(p::SyncPlot)

function display_blink(p::SyncPlot)
    p.window = Blink.Window()
    Blink.body!(p.window, p.scope)
end

function Base.close(p::SyncPlot)
    if p.window !== nothing && Blink.active(p.window)
        close(p.window)
    end
end

function send_command(scope, cmd, args...)
    # The handler for _commands is set up when plot is constructed
    scope["_commands"][] = [cmd, args...]
    nothing
end

# ----------------------- #
# Plotly.js api functions #
# ----------------------- #

# for each of these we first update the Julia object, then update the display

function restyle!(
        plt::SyncPlot, ind::Union{Int,AbstractVector{Int}},
        update::AbstractDict=Dict();
        kwargs...)
        restyle!(plt.plot, ind, update; kwargs...)
        send_command(plt.scope, :restyle, merge(update, prep_kwargs(kwargs)), ind .- 1)
end

function restyle!(plt::SyncPlot, update::AbstractDict=Dict(); kwargs...)
    restyle!(plt.plot, update; kwargs...)
    send_command(plt.scope, :restyle, merge(update, prep_kwargs(kwargs)))
end

function relayout!(plt::SyncPlot, update::AbstractDict=Dict(); kwargs...)
    relayout!(plt.plot, update; kwargs...)
    send_command(plt.scope, :relayout, merge(update, prep_kwargs(kwargs)))
end

function react!(plt::SyncPlot, data::AbstractVector{<:AbstractTrace}, layout::Layout)
    react!(plt.plot, data, layout)
    send_command(plt.scope, :react, data, layout)
end

function update!(
        plt::SyncPlot, ind::Union{Int,AbstractVector{Int}},
        update::AbstractDict=Dict();
        layout::Layout=Layout(),
        kwargs...)
    update!(plt.plot, ind, update; layout=layout, kwargs...)
    send_command(plt.scope, :update, merge(update, prep_kwargs(kwargs)), layout, ind .- 1)
end

function update!(
        plt::SyncPlot, update::AbstractDict=Dict(); layout::Layout=Layout(),
        kwargs...
    )
    update!(plt.plot, update; layout=layout, kwargs...)
    send_command(plt.scope, :update, merge(update, prep_kwargs(kwargs)), layout)
end

function addtraces!(plt::SyncPlot, traces::AbstractTrace...)
    addtraces!(plt.plot, traces...)
    send_command(plt.scope, :addTraces, traces)
end

function addtraces!(plt::SyncPlot, i::Int, traces::AbstractTrace...)
    addtraces!(plt.plot, i, traces...)
    send_command(plt.scope, :addTraces, traces, i-1)
end

function deletetraces!(plt::SyncPlot, inds::Int...)
    deletetraces!(plt.plot, inds...)
    send_command(plt.scope, :deleteTraces, collect(inds) .- 1)
end

function movetraces!(plt::SyncPlot, to_end::Int...)
    movetraces!(plt.plot, to_end...)
    send_command(plt.scope, :moveTraces, traces, collect(to_end) .- 1)
end

function movetraces!(
        plt::SyncPlot, src::AbstractVector{Int}, dest::AbstractVector{Int}
    )
    movetraces!(plt.plot, src, dest)
    send_command(plt.scope, :moveTraces, src .- 1, dest .- 1)
end

function redraw!(plt::SyncPlot)
    redraw!(plt.plot)
    send_command(plt.scope, :redraw)
end

function purge!(plt::SyncPlot)
    purge!(plt.plot)
    send_command(plt.scope, :purge)
end

function to_image(plt::SyncPlot; kwargs...)
    to_image(plt.plot)
    plt.scope["image"][] = ""  # reset
    plt.scope["_toImage"][] = Dict(kwargs)

    tries = 0
    while length(plt.scope["image"][]) == 0
        tries == 10 && error("Could not get image")
        sleep(0.25)
        tries += 1
    end
    return plt["image"].val
end

function download_image(plt::SyncPlot; kwargs...)
    download_image(plt.plot)
    plt.scope["_downloadImage"][] = Dict(kwargs)
    nothing
end

# unexported (by plotly.js) api methods
function extendtraces!(plt::SyncPlot, update::AbstractDict,
              indices::AbstractVector{Int}=[1], maxpoints=-1;)
    extendtraces!(plt.plot, update, indices, maxpoints)
    send_command(
        plt.scope, :extendTraces, prep_kwargs(update), indices .- 1, maxpoints
    )
end

function prependtraces!(plt::SyncPlot, update::AbstractDict,
               indices::AbstractVector{Int}=[1], maxpoints=-1;)
    prependtraces!(plt.plot, update, indices, maxpoints)
    send_command(
        plt.scope, :prependTraces, prep_kwargs(update), indices .- 1, maxpoints
    )
end

for f in [:restyle, :relayout, :update, :addtraces, :deletetraces,
          :movetraces, :redraw, :extendtraces, :prependtraces, :purge, :react]
    f! = Symbol(f, "!")
    @eval function $(f)(plt::SyncPlot, args...; kwargs...)
        out = SyncPlot(deepcopy(plt.plot))
        $(f!)(out, args...; kwargs...)
        out
    end
end

# add extra same extra methods we have on ::Plot for these functions
for f in (:extendtraces!, :prependtraces!)
    @eval begin
        function $(f)(p::SyncPlot, inds::Vector{Int}=[0],
                      maxpoints=-1; update...)
            d = Dict()
            for (k, v) in update
                d[k] = _tovec(v)
            end
            ($f)(p, d, inds, maxpoints)
        end

        function $(f)(p::SyncPlot, ind::Int, maxpoints=-1;
                      update...)
            ($f)(p, [ind], maxpoints; update...)
        end

        function $(f)(p::SyncPlot, update::AbstractDict, ind::Int,
                      maxpoints=-1)
            ($f)(p, update, [ind], maxpoints)
        end
    end
end

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

const js_default = Ref(:local)

function savehtml(io::IO, p::Plot, js::Symbol=js_default[])

    if js == :local
        script_txt = "<script src=\"$(_js_path)\"></script>"
    elseif js == :remote
        script_txt = "<script src=\"$(_js_cdn_path)\"></script>"
    elseif js == :embed
        script_txt = "<script>$(read(_js_path, String))</script>"
    else
        msg = """
        Unknown value for argument js: $js.
        Possible choices are `:local`, `:remote`, `:embed`
            """
        throw(ArgumentError(msg))
    end

    print(io, """
    <html>
    <head>
         $script_txt
    </head>
    <body>
         $(html_body(p))
    </body>
    </html>
    """)

end

function savehtml(p::Plot, fn::AbstractString, js::Symbol=js_default[])
    open(fn, "w") do f
        savehtml(f, p, js)
    end
end

"""
    savehtml(io::IO, p::Union{Plot,SyncPlot}, js::Symbol=js_default[])
    savehtml(p::Union{Plot,SyncPlot}, fn::AbstractString, js::Symbol=js_default[])

Save plot to standalone html file suitable for including in a website or
opening in a browser

Can either be written to an arbitrary IO stream, or saved to a file noted with
a string `fn`.

The `js` argument can be one of

- `:local`: Reference the local plotly.js file included in this Julia package
    Pros: small file size, offline viewing. Cons: Can't share with others or
    move to different machine..
- `:remote`: Reference plotly.js from a CDN. Pros small file size, move to
    other machine. Cons: need internet access to fetch from CDN
- `:embed`: Embed the entirety of your local copy of plotly.js in the
    outputted file. Pros: offline viewing, move to other machine. Con: large
    file size (adds about 2.7 MB)

The default is `:local`
"""
savehtml

for mime in ["text/plain", "application/vnd.plotly.v1+json"]
    function Base.show(io::IO, m::MIME{Symbol(mime)}, p::SyncPlot, args...)
        show(io, m, p.plot, args...)
    end
end

PlotlyBase.savejson(sp::SyncPlot, fn::String) = PlotlyBase.savejson(sp.plot, fn)
savehtml(io::IO, p::SyncPlot, js::Symbol=js_default[]) = savehtml(io, p.plot, js)
savehtml(p::SyncPlot, fn::AbstractString, js::Symbol=js_default[]) = savehtml(p.plot, fn, js)


text:
using PyPlot
using Colors

###############
# visualization
###############
export to_file, save_fig, make_tight_colorbar, visualize_pdf

to_file(s) = replace(lowercase(s), " " => "_")

function save_fig(fig_path, title)
  tikzplotlib.save("$(fig_path)/$(to_file(title)).tex",
                       axis_height="\\figureheight",
                       axis_width="\\figurewidth",
                       show_info=false)
end

# standard markers
station_center = "+"
sensor_center = (5, 2, 0)
selected_sensor = "*"
event_marker = "x"

agent_scale = 32.0

station_style = "--"

rgb_tuple(color::RGB) = (red(color), green(color), blue(color))

generate_colors(agents) = distinguishable_colors(length(agents) + 1)[2:end]

function visualize_agents(agents, colors)
  centers = map(get_center, agents);

  map(agents, colors) do agent, color
    center = get_center(agent)

#   scatter([center[1]], [center[2]], color = rgb_tuple(color),
#           s = 4*agent_scale, marker = station_center)

    map(get_block(agent)) do element
      plot_element(element; color = rgb_tuple(color))
    end

    plot_agent(agent, color=rgb_tuple(color))
  end

  nothing
end

function visualize_solution(p::PartitionProblem, X, agent_colors)

  elements = map(X) do x
    get_element(p.partition_matroid, x)
  end

  colors = [agent_colors[x[1]] for x in X]
  visualize_solution(elements, colors)
end

visualize_solution(p::PartitionProblem, s::Solution, agent_colors) =
  visualize_solution(p, s.elements, agent_colors)

# visualization
function circle_points(p; radius=1, scale=1, dth=0.1)
  scale*hcat(map(x->p+[radius*cos(x);radius*sin(x)], 0:dth:(2*pi)+dth)...)
end

function plot_circle(p; scale=1, radius=1, color="k", linestyle="-", linewidth=1.0, kwargs...)
  c = circle_points(p, scale=scale, radius=radius)
  plot(c[1,:][:], c[2,:][:], color=color, linestyle=linestyle, linewidth=linewidth, kwargs...)
end

function make_tight_colorbar(image)
  ax = gca()

  divider = ag.make_axes_locatable(ax)
  cax = divider.append_axes("right", "5%", pad="3%")
  colorbar(image, cax = cax)

  sca(ax)
end

function visualize_pdf(density::Matrix; limits = [0 1; 0 1], cmap = "viridis",
                       show_colorbar = true, kwargs...)
  ret = []

  xlim = limits[1,:]
  ylim = limits[2,:]

  image = imshow(density'; cmap=cmap, vmin=0.0, vmax=maximum(density[:]),
                 extent=[xlim[1], xlim[2], ylim[1], ylim[2]],
                 interpolation="nearest", origin="lower", kwargs...)
  push!(ret, image)

  if show_colorbar
    push!(ret, make_tight_colorbar(image))
  end

  ret
end

function visualize_pdf(fun; n = 1000, limits = [0 1; 0 1], kwargs...)
  xlim = limits[1,:]
  ylim = limits[2,:]

  volume = (xlim[2] - xlim[1]) * (ylim[2] - ylim[1])

  density = [evaluate_pdf(fun, [x,y]) for x in range(xlim[1], stop=xlim[2], length=n),
                                 y in range(ylim[1], stop=ylim[2], length=n)]

  # the density is normalized to be a proper probability density on the input
  # interval so values can be anywhere on real+, keep zero though
  density .= density .* length(density) / sum(density) * volume

  visualize_pdf(density; limits=limits, kwargs...)
end

ext: jl

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
using DataFrames
using Distributions
using NearestNeighbors
using StaticArrays
using StatsBase

struct ShapePrior
    std_values::MeanVec;
    std_value_stds::MeanVec;
    n_samples::Int;
end

distributions(prior::ShapePrior) = Normal.(prior.std_values, prior.std_value_stds)
rand(prior::ShapePrior) = rand.(distributions(prior))

var_posterior(prior::ShapePrior, eigen_values::T where T <: Union{MeanVec, StaticArray{Tuple{2},Float64,1}}; n_samples::TR where TR <: Real) =
    var_posterior.(prior.std_values, prior.std_value_stds, eigen_values; n_samples=n_samples, prior_n_samples=prior.n_samples)

# f(dx) = sign(dx) * sqrt(|dx/std|) * std
function var_posterior(prior_std::Float64, prior_std_std::Float64, eigen_value::Float64; n_samples::TR1 where TR1 <: Real, prior_n_samples::TR2 where TR2 <: Real)
    d_std = (sqrt(eigen_value) - prior_std)
    std_adj = (prior_std + sign(d_std) * (sqrt(abs(d_std) / prior_std_std + 1) - 1) * prior_std_std)

    return ((prior_n_samples * prior_std + n_samples * std_adj) / (prior_n_samples + n_samples)) ^ 2
end

function sample_var(d::Normal)
    @assert d.μ > 0
    while true
        v = rand(d)
        if v > 0
            return v.^2
        end
    end
end
sample_var(prior::ShapePrior) = sample_var.(distributions(prior))

mutable struct Component
    position_params::MvNormalF;
    composition_params::CategoricalSmoothed;
    n_samples::Int;
    prior_probability::Float64;

    shape_prior::Union{Nothing, ShapePrior};
    n_molecules_per_segment::Dict{Int, Int};

    guid::Int;

    Component(position_params::MvNormalF, composition_params::CategoricalSmoothed;
              n_samples::Int=0, shape_prior::Union{Nothing, ShapePrior}=nothing, guid::Int=-1) =
        new(position_params, composition_params, n_samples, 1.0, shape_prior, Dict{Int, Int}(), guid)
end

function maximize!(c::Component, pos_data::T1 where T1 <: AbstractMatrix{Float64}, comp_data::T2 where T2 <: AbstractVector{Int}, conf_data::T3 where T3 <: AbstractVector{Float64})
    c.n_samples = size(pos_data, 2)
    maximize!(c.composition_params, comp_data, conf_data);
    maximize!(c.position_params, pos_data, conf_data);

    if c.shape_prior !== nothing
        try
            adjust_cov_by_prior!(c.position_params.Σ, c.shape_prior; n_samples=sum(conf_data))
        catch
            @show c.position_params.Σ
            @show c.shape_prior
            @show sum(conf_data)
            rethrow()
        end
    end

    return c
end

pdf(params::Component, x::Float64, y::Float64, gene::Int64; use_smoothing::Bool=true)::Float64 =
    pdf(params.position_params, x, y) * pdf(params.composition_params, gene; use_smoothing=use_smoothing)

function adjust_cov_by_prior!(Σ::CovMat, prior::ShapePrior; n_samples::TR where TR <: Real)
    if (Σ[2, 1] / max(Σ[1, 1], Σ[2, 2])) < 1e-5 # temporary fix untill https://github.com/JuliaArrays/StaticArrays.jl/pull/694 is merged
        Σ[1, 2] = Σ[2, 1] = 0.0
    end

    fact = eigen(Σ)
    eigen_values_posterior = var_posterior(prior, fmax.(fact.values, 0.0); n_samples=n_samples)
    Σ .= fact.vectors * SMatrix{2, 2, Float64}(diagm(0 => eigen_values_posterior)) * inv(fact.vectors)
    Σ[1, 2] = Σ[2, 1]

    return adjust_cov_matrix!(Σ)
end

function normal_posterior!(μ::MeanVec, Σ::CovMat, μ_prior::MeanVec, Σ_prior::CovMat; n::TR1 where TR1 <: Real, n_prior::TR2 where TR2 <: Real)
    Σ = (n .* Σ .+ n_prior .* Σ_prior .+ (n * n_prior / (n + n_prior)) .* ((μ .- μ_prior) * (μ .- μ_prior)')) ./ (n + n_prior)
    μ .= (n .* μ .+ n_prior .* μ_prior) ./ (n + n_prior)
    adjust_cov_matrix!(Σ)
    return nothing
end

position(c::Component) = c.position_params.μ

eigen_values(c::Component) = eigen(shape(c.position_params)).values

max_stars_repo_name: YinuoJin/Baysor

max_issues_repo_name: YinuoJin/Baysor

max_forks_repo_name: YinuoJin/Baysor


text:
using Test
using GLM
using Statistics
using DataFrames
using Calculus

@testset "HWunconstrained.jl" begin

@testset "basics" begin

	@testset "Test Data Construction" begin

		# create data
		d = makeData(18)

		# test dataframe has 18 rows

		# test generated y vector has a mean less than 1
	end

	@testset "Test Return value of likelihood" begin

	end

	@testset "Test return value of gradient is nothing" begin
		# gradient should not return anything,
		# but modify a vector in place.

	end

	@testset "gradient vs finite difference" begin
	end
end

@testset "test maximization results" begin

	ttol = 2e-1

	@testset "maximize returns approximate result" begin
		m = HWunconstrained.maximize_like();
		d = makeData();
		@test maximum(abs,m.minimizer .- d["beta"]) < ttol
	end

	@testset "maximize_grad returns accurate result" begin
		m = HWunconstrained.maximize_like_grad();
		d = makeData();
		@test maximum(abs,m.minimizer .- d["beta"]) < ttol
	end

	@testset "maximize_grad_hess returns accurate result" begin
		m = HWunconstrained.maximize_like_grad_hess();
		d = makeData();
		@test maximum(abs,m.minimizer .- d["beta"]) < ttol
	end

	@testset "gradient is close to zero at max like estimate" begin
		m = HWunconstrained.maximize_like_grad();
		d = makeData()
		gradvec = ones(length(d["beta"]))
		r = HWunconstrained.grad!(gradvec,m.minimizer,d)

		@test r == nothing

		@test maximum(abs,gradvec) < 1e-6

	end

end

@testset "test against GLM" begin

	@testset "estimates vs GLM" begin

	end

	@testset "standard errors vs GLM" begin

	end

end

end


text:
###############################################################################
#
# Plot linear regression on axis
#
###############################################################################
"""
    plot_line_regress!(
                ax::PyObject,
                xs::Array,
                ys::Array;
                linestyle::String = "-",
                intercept::Bool = true,
                interval::Bool = false,
                color::String = "red",
                alpha::Number = 0.3
    )

Plor linear regression and confidence interval on the axis, given
- `ax` Given axis
- `xs` Array of x
- `ys` Array of y
- `linestyle` Optional. Line style for the regression curve ("-" by default)
- `intercept` Optional: if true, fit the data with an intercept
- `interval` Optional: if true, plot the confidence interval of fitted y
- `color` Color the fitted curve
- `alpha` Transparency of the confidence interval (same color as curve)
"""
function plot_line_regress!(
            ax::PyObject,
            xs::Array,
            ys::Array;
            linestyle::String = "-",
            intercept::Bool = true,
            interval::Bool = false,
            color::String = "red",
            alpha::Number = 0.3
)
    # make linear regression
    lr = line_regress(xs, ys, intercept=intercept, sorting=true);

    # plot the fittings
    ax.plot(lr.df.X, lr.df.predY, linestyle=linestyle, color=color);
    if interval
        ax.fill_between(lr.df.X,
                        lr.df.lowerY,
                        lr.df.upperY,
                        facecolor=color,
                        alpha=alpha);
    end

    return nothing
end


text:
function recvreq(socket)
  while true
    addr, msg = recvfrom(socket)
    println("RECEIVED")
    readval = ProtoBuf.readproto(IOBuffer(msg), Raft.RaftRPC.RPCRequest())
    return (addr, readval)
  end
end

function send_aereply(socket, ip, port, my_id, aer)
  reply = Raft.RaftRPC.RPCReply(_type=Int32(0),server_id=my_id)
  ib = IOBuffer()
  ProtoBuf.writeproto(ib, aer)
  ProtoBuf.set_field!(reply, :reply, takebuf_array(ib))

  println("\tSending AER reply")
  outb = IOBuffer()
  ProtoBuf.writeproto(outb, reply)
  send(socket, ip, port, takebuf_array(outb))
end

function run_follower(my_config::ServerConfig)
  println("Starting Listener")
  socket = UDPSocket()
  bind(socket, ip"127.0.0.1", my_config.port)
  while true
    addr, rpcreq = recvreq(socket)
    if !ProtoBuf.has_field(rpcreq, :_type) || !ProtoBuf.has_field(rpcreq, :request)
        println("\tDropped")
        continue
    elseif !ProtoBuf.has_field(rpcreq, :server_id)
      println("\tNo return address")
    end

    request_type = ProtoBuf.lookup(Raft.RaftRPC.ExchangeType, ProtoBuf.get_field(rpcreq, :_type))
    sender = my_config.servers[ProtoBuf.get_field(rpcreq, :server_id)]
    assert(sender.ipaddress == addr)
    if request_type == :APPENDENTRIES
      println("\tAER")
      aer = ProtoBuf.readproto(IOBuffer(ProtoBuf.get_field(rpcreq, :request)), Raft.RaftRPC.AppendEntriesRequest())
      send_aereply(socket, sender.ipaddress, sender.port, my_config.id, Raft.RaftRPC.AppendEntriesReply(term=UInt64(1),success=true))
    elseif request_type == :REQUESTVOTE
      rvr = ProtoBuf.readproto(IOBuffer(ProtoBuf.get_field(rpcreq, :request)), Raft.RaftRPC.RequestVoteRequest())
      println("\tRVR ", rvr)
    elseif request_type == :INSTALLSNAPSHOT
      isr = ProtoBuf.readproto(IOBuffer(ProtoBuf.get_field(rpcreq, :request)), Raft.RaftRPC.InstallSnapshotRequest())
      println("\tISR ", isr)
    else
      println("\tUNKNOWN REQUEST TYPE")
    end
  end
end


text:
# LLVM IR generation

function irgen(@nospecialize(job::CompilerJob), method_instance::Core.MethodInstance, world)
    entry, mod = @timeit_debug to "emission" compile_method_instance(job, method_instance, world)
    ctx = context(mod)

    # clean up incompatibilities
    @timeit_debug to "clean-up" begin
        for llvmf in functions(mod)
            # only occurs in debug builds
            delete!(function_attributes(llvmf), EnumAttribute("sspstrong", 0, ctx))

            if VERSION < v"1.5.0-DEV.393"
                # make function names safe for ptxas
                llvmfn = LLVM.name(llvmf)
                if !isdeclaration(llvmf)
                    llvmfn′ = safe_name(llvmfn)
                    if llvmfn != llvmfn′
                        LLVM.name!(llvmf, llvmfn′)
                        llvmfn = llvmfn′
                    end
                end
            end

            if Sys.iswindows()
                personality!(llvmf, nothing)
            end
        end

        # remove the exception-handling personality function
        if Sys.iswindows() && "__julia_personality" in functions(mod)
            llvmf = functions(mod)["__julia_personality"]
            @compiler_assert isempty(uses(llvmf)) job
            unsafe_delete!(mod, llvmf)
        end
    end

    # target-specific processing
    process_module!(job, mod)

    # sanitize function names
    # FIXME: Julia should do this, but apparently fails (see maleadt/LLVM.jl#201)
    for f in functions(mod)
        LLVM.isintrinsic(f) && continue
        llvmfn = LLVM.name(f)
        startswith(llvmfn, "julia.") && continue # Julia intrinsics
        startswith(llvmfn, "llvm.") && continue # unofficial LLVM intrinsics
        llvmfn′ = safe_name(llvmfn)
        if llvmfn != llvmfn′
            @assert !haskey(functions(mod), llvmfn′)
            LLVM.name!(f, llvmfn′)
        end
    end

    # rename and process the entry point
    if job.source.name !== nothing
        LLVM.name!(entry, safe_name(string("julia_", job.source.name)))
    end
    if job.source.kernel
        LLVM.name!(entry, mangle_call(entry, job.source.tt))
    end
    entry = process_entry!(job, mod, entry)

    # minimal required optimization
    @timeit_debug to "rewrite" ModulePassManager() do pm
        global current_job
        current_job = job

        linkage!(entry, LLVM.API.LLVMExternalLinkage)

        # internalize all functions, but keep exported global variables
        exports = String[LLVM.name(entry)]
        for gvar in globals(mod)
            push!(exports, LLVM.name(gvar))
        end
        internalize!(pm, exports)

        can_throw(job) || add!(pm, ModulePass("LowerThrow", lower_throw!))

        add_lowering_passes!(job, pm)

        run!(pm, mod)

        # NOTE: if an optimization is missing, try scheduling an entirely new optimization
        # to see which passes need to be added to the target-specific list
        #     LLVM.clopts("-print-after-all", "-filter-print-funcs=$(LLVM.name(entry))")
        #     ModulePassManager() do pm
        #         add_library_info!(pm, triple(mod))
        #         add_transform_info!(pm, tm)
        #         PassManagerBuilder() do pmb
        #             populate!(pm, pmb)
        #         end
        #         run!(pm, mod)
        #     end
    end

    return mod, entry
end

## name mangling

# we generate function names that look like C++ functions, because many NVIDIA tools
# support them, e.g., grouping different instantiations of the same kernel together.

function mangle_param(t, substitutions)
    t == Nothing && return "v"

    if isa(t, DataType) || isa(t, Core.Function)
        tn = safe_name(t)

        # handle substitutions
        sub = findfirst(isequal(tn), substitutions)
        if sub === nothing
            str = "$(length(tn))$tn"
            push!(substitutions, tn)
        elseif sub == 1
            str = "S_"
        else
            str = "S$(sub-2)_"
        end

        # encode typevars as template parameters
        if !isempty(t.parameters)
            str *= "I"
            for t in t.parameters
                str *= mangle_param(t, substitutions)
            end
            str *= "E"
        end

        str
    elseif isa(t, Integer)
        "Li$(t)E"
    else
        tn = safe_name(t)
        "$(length(tn))$tn"
    end
end

function mangle_call(f, tt)
    fn = safe_name(f)
    str = "_Z$(length(fn))$fn"

    substitutions = String[]
    for t in tt.parameters
        str *= mangle_param(t, substitutions)
    end

    return str
end

# make names safe for ptxas
safe_name(fn::String) = replace(fn, r"[^A-Za-z0-9_]"=>"_")
safe_name(f::Union{Core.Function,DataType}) = safe_name(String(nameof(f)))
safe_name(f::LLVM.Function) = safe_name(LLVM.name(f))
safe_name(x) = safe_name(repr(x))

## exception handling

# this pass lowers `jl_throw` and friends to GPU-compatible exceptions.
# this isn't strictly necessary, but has a couple of advantages:
# - we can kill off unused exception arguments that otherwise would allocate or invoke
# - we can fake debug information (lacking a stack unwinder)
#
# once we have thorough inference (ie. discarding `@nospecialize` and thus supporting
# exception arguments) and proper debug info to unwind the stack, this pass can go.
function lower_throw!(mod::LLVM.Module)
    job = current_job::CompilerJob
    ctx = context(mod)
    changed = false
    @timeit_debug to "lower throw" begin

    throw_functions = [
        # unsupported runtime functions that are used to throw specific exceptions
        "jl_throw"                      => "exception",
        "jl_error"                      => "error",
        "jl_too_few_args"               => "too few arguments exception",
        "jl_too_many_args"              => "too many arguments exception",
        "jl_type_error"                 => "type error",
        "jl_type_error_rt"              => "type error",
        "jl_undefined_var_error"        => "undefined variable error",
        "jl_bounds_error"               => "bounds error",
        "jl_bounds_error_v"             => "bounds error",
        "jl_bounds_error_int"           => "bounds error",
        "jl_bounds_error_tuple_int"     => "bounds error",
        "jl_bounds_error_unboxed_int"   => "bounds error",
        "jl_bounds_error_ints"          => "bounds error",
        "jl_eof_error"                  => "EOF error",
        # Julia-level exceptions that use unsupported inputs like interpolated strings
        r"julia_throw_exp_domainerror_\d+"      => "DomainError",
        r"julia_throw_complex_domainerror_\d+"  => "DomainError",
        r"julia_throw_domerr_powbysq_\d+"       => "DomainError",
        r"julia_throw_overflowerr_binaryop_\d+" => "OverflowError",
        r"julia_throw_overflowerr_negation_\d+" => "OverflowError",
        r"julia_throw_inexacterror_\d+"         => "InexactError",
        r"julia_throw_boundserror_\d+"          => "BoundsError",
    ]

    for f in functions(mod)
        fn = LLVM.name(f)
        for (throw_fn, name) in throw_functions
            occursin(throw_fn, fn) || continue

            for use in uses(f)
                call = user(use)::LLVM.CallInst

                # replace the throw with a PTX-compatible exception
                let builder = Builder(ctx)
                    position!(builder, call)
                    emit_exception!(builder, name, call)
                    dispose(builder)
                end

                # remove the call
                call_args = collect(operands(call))[1:end-1] # last arg is function itself
                unsafe_delete!(LLVM.parent(call), call)

                # HACK: kill the exceptions' unused arguments
                for arg in call_args
                    # peek through casts
                    if isa(arg, LLVM.AddrSpaceCastInst)
                        cast = arg
                        arg = first(operands(cast))
                        isempty(uses(cast)) && unsafe_delete!(LLVM.parent(cast), cast)
                    end

                    if isa(arg, LLVM.Instruction) && isempty(uses(arg))
                        unsafe_delete!(LLVM.parent(arg), arg)
                    end
                end

                changed = true
            end

            @compiler_assert isempty(uses(f)) job
            break
         end
     end

    end
    return changed
end

# report an exception in a GPU-compatible manner
#
# the exact behavior depends on the debug level. in all cases, a `trap` will be emitted, On
# debug level 1, the exception name will be printed, and on debug level 2 the individual
# stack frames (as recovered from the LLVM debug information) will be printed as well.
function emit_exception!(builder, name, inst)
    job = current_job::CompilerJob
    bb = position(builder)
    fun = LLVM.parent(bb)
    mod = LLVM.parent(fun)
    ctx = context(mod)

    # report the exception
    if Base.JLOptions().debug_level >= 1
        name = globalstring_ptr!(builder, name, "exception")
        if Base.JLOptions().debug_level == 1
            call!(builder, Runtime.get(:report_exception), [name])
        else
            call!(builder, Runtime.get(:report_exception_name), [name])
        end
    end

    # report each frame
    if Base.JLOptions().debug_level >= 2
        rt = Runtime.get(:report_exception_frame)
        ft = convert(LLVM.FunctionType, rt, ctx)
        bt = backtrace(inst)
        for (i,frame) in enumerate(bt)
            idx = ConstantInt(parameters(ft)[1], i)
            func = globalstring_ptr!(builder, String(frame.func), "di_func")
            file = globalstring_ptr!(builder, String(frame.file), "di_file")
            line = ConstantInt(parameters(ft)[4], frame.line)
            call!(builder, rt, [idx, func, file, line])
        end
    end

    # signal the exception
    call!(builder, Runtime.get(:signal_exception))

    emit_trap!(job, builder, mod, inst)
end

function emit_trap!(@nospecialize(job::CompilerJob), builder, mod, inst)
    ctx = context(mod)
    trap = if haskey(functions(mod), "llvm.trap")
        functions(mod)["llvm.trap"]
    else
        LLVM.Function(mod, "llvm.trap", LLVM.FunctionType(LLVM.VoidType(ctx)))
    end
    call!(builder, trap)
end

## kernel promotion

@enum ArgumentCC begin
    BITS_VALUE  # bitstype, passed as value
    BITS_REF    # bitstype, passed as pointer
    MUT_REF     # jl_value_t*, or the anonymous equivalent
    GHOST       # not passed
end

function classify_arguments(@nospecialize(job::CompilerJob), codegen_f::LLVM.Function)
    codegen_ft = eltype(llvmtype(codegen_f)::LLVM.PointerType)::LLVM.FunctionType
    source_sig = Base.signature_type(job.source.f, job.source.tt)::Type

    codegen_types = parameters(codegen_ft)
    source_types = [source_sig.parameters...]

    args = []
    codegen_i = 1
    for (source_i, source_typ) in enumerate(source_types)
        if isghosttype(source_typ) ||
           (VERSION >= v"1.5.0-DEV.581" && Core.Compiler.isconstType(source_typ))
            push!(args, (cc=GHOST, typ=source_typ))
            continue
        end

        codegen_typ = codegen_types[codegen_i]
        if codegen_typ isa LLVM.PointerType && !issized(eltype(codegen_typ))
            push!(args, (cc=MUT_REF, typ=source_typ,
                         codegen=(typ=codegen_typ, i=codegen_i)))
        elseif codegen_typ isa LLVM.PointerType && issized(eltype(codegen_typ)) &&
               !(source_typ <: Ptr) && !(VERSION >= v"1.5-" && source_typ <: Core.LLVMPtr)
            push!(args, (cc=BITS_REF, typ=source_typ,
                         codegen=(typ=codegen_typ, i=codegen_i)))
        else
            push!(args, (cc=BITS_VALUE, typ=source_typ,
                         codegen=(typ=codegen_typ, i=codegen_i)))
        end
        codegen_i += 1
    end

    return args
end

## byval lowering

# some back-ends don't support byval, or support it badly
# https://reviews.llvm.org/D79744

# generate a kernel wrapper to fix & improve argument passing
function lower_byval(@nospecialize(job::CompilerJob), mod::LLVM.Module, entry_f::LLVM.Function)
    ctx = context(mod)
    entry_ft = eltype(llvmtype(entry_f)::LLVM.PointerType)::LLVM.FunctionType
    @compiler_assert return_type(entry_ft) == LLVM.VoidType(ctx) job

    args = classify_arguments(job, entry_f)
    filter!(args) do arg
        arg.cc != GHOST
    end

    # generate the wrapper function type & definition
    wrapper_types = LLVM.LLVMType[]
    for arg in args
        typ = if arg.cc == BITS_REF
            eltype(arg.codegen.typ)
        else
            convert(LLVMType, arg.typ, ctx)
        end
        push!(wrapper_types, typ)
    end
    wrapper_fn = LLVM.name(entry_f)
    LLVM.name!(entry_f, wrapper_fn * ".inner")
    wrapper_ft = LLVM.FunctionType(LLVM.VoidType(ctx), wrapper_types)
    wrapper_f = LLVM.Function(mod, wrapper_fn, wrapper_ft)

    # emit IR performing the "conversions"
    let builder = Builder(ctx)
        entry = BasicBlock(wrapper_f, "entry", ctx)
        position!(builder, entry)

        wrapper_args = Vector{LLVM.Value}()

        # perform argument conversions
        for arg in args
            if arg.cc == BITS_REF
                # copy the argument value to a stack slot, and reference it.
                ptr = alloca!(builder, eltype(arg.codegen.typ))
                if LLVM.addrspace(arg.codegen.typ) != 0
                    ptr = addrspacecast!(builder, ptr, arg.codegen.typ)
                end
                store!(builder, parameters(wrapper_f)[arg.codegen.i], ptr)
                push!(wrapper_args, ptr)
            else
                push!(wrapper_args, parameters(wrapper_f)[arg.codegen.i])
                for attr in collect(parameter_attributes(entry_f, arg.codegen.i))
                    push!(parameter_attributes(wrapper_f, arg.codegen.i), attr)
                end
            end
        end

        call!(builder, entry_f, wrapper_args)

        ret!(builder)

        dispose(builder)
    end

    # early-inline the original entry function into the wrapper
    push!(function_attributes(entry_f), EnumAttribute("alwaysinline", 0, ctx))
    linkage!(entry_f, LLVM.API.LLVMInternalLinkage)

    fixup_metadata!(entry_f)
    ModulePassManager() do pm
        always_inliner!(pm)
        run!(pm, mod)
    end

    return wrapper_f
end

# HACK: get rid of invariant.load and const TBAA metadata on loads from pointer args,
#       since storing to a stack slot violates the semantics of those attributes.
# TODO: can we emit a wrapper that doesn't violate Julia's metadata?
function fixup_metadata!(f::LLVM.Function)
    for param in parameters(f)
        if isa(llvmtype(param), LLVM.PointerType)
            # collect all uses of the pointer
            worklist = Vector{LLVM.Instruction}(user.(collect(uses(param))))
            while !isempty(worklist)
                value = popfirst!(worklist)

                # remove the invariant.load attribute
                md = metadata(value)
                if haskey(md, LLVM.MD_invariant_load)
                    delete!(md, LLVM.MD_invariant_load)
                end
                if haskey(md, LLVM.MD_tbaa)
                    delete!(md, LLVM.MD_tbaa)
                end

                # recurse on the output of some instructions
                if isa(value, LLVM.BitCastInst) ||
                   isa(value, LLVM.GetElementPtrInst) ||
                   isa(value, LLVM.AddrSpaceCastInst)
                    append!(worklist, user.(collect(uses(value))))
                end

                # IMPORTANT NOTE: if we ever want to inline functions at the LLVM level,
                # we need to recurse into call instructions here, and strip metadata from
                # called functions (see CUDAnative.jl#238).
            end
        end
    end
end


text:
using StructuralCausalModels, Test

ProjDir = @__DIR__
cd(ProjDir) #do

d = OrderedDict(
  :w => :s,
  :d => [:a, :w, :m],
  :m => [:a, :s],
  :a => :s
);

dag = DAG("test_open_paths_01", d);
fname = joinpath(ProjDir, "test_open_paths_01.dot")
to_graphviz(dag, fname)
#Sys.isapple() && run(`open -a GraphViz.app $(fname)`)

to_dagitty(d)

bs = basis_set(dag)
ap  = all_paths(dag, :w, :d)
bp = backdoor_paths(dag, ap, :w)

cs = Vector{Symbol}[]
syms = syms_in_paths(bp, :w, :d)

function adjsets(d::DAG, paths, syms)
  lsyms = deepcopy(syms)
  adjustment_sets = Vector{Symbol}[]
  for s in lsyms
    if sym_in_all_paths(bp, s)
      if length(open_paths(dag, bp, [s])) == 0
        #println("$s closes all paths")
        append!(adjustment_sets, [[s]])
        setdiff!(lsyms, [s])
      end
    end
  end
  len = 2
  local csyms
  while length(lsyms) >= len
    csyms = collect(combinations(lsyms, len))
    for s in csyms
      if length(open_paths(dag, bp, s)) == 0
        append!(adjustment_sets, [s])
        setdiff!(lsyms, s)
      end
    end
    len += 1
  end
  adjustment_sets
end

@testset "Open_path_01" begin

  adjs = adjsets(dag, bp, syms)
  @test adjs == adjustment_sets(dag, :w, :d)

end


text:
@doc raw"""
    Oblique{N,M} <: AbstractPowerManifold

The oblique manifold $\mathcal{OB}(n,m)$ is the set of real-valued matrices with unit norm
column endowed with the metric from the embedding. This yields exactly the same metric as
considering the product metric of the unit norm vectors, i.e. [`PowerManifold`](@ref) of the
$(n-1)$-dimensional [`Sphere`](@ref).

The [`Sphere`](@ref) is stored internally within `M.manifold`, such that all functions of
[`AbstractPowerManifold`](@ref) can be used directly.

# Constructor

    Oblique(n,m)

Generate the manifold of matrices $\mathbb R^{n × m}$ such that the $m$ columns are unit
vectors, i.e. from the [`Sphere`](@ref)`(n-1)`.
"""
struct Oblique{N,M} <: AbstractPowerManifold{Sphere{N},ArrayPowerRepresentation}
    manifold::Sphere{N}
end

Oblique(n::Int, m::Int) = Oblique{n - 1,m}(Sphere(n - 1))

^(M::Sphere{N}, m::Int) where {N} = Oblique{N,m}(M)

@doc raw"""
    check_manifold_point(M::Oblique{n,m},p)

Checks whether `p` is a valid point on the [`Oblique`](@ref)`{m,n}` `M`, i.e. is a matrix
of `m` unit columns from $\mathbb R^{n+1}$, i.e. each column is a point from
[`Sphere`](@ref)`(n)`.
"""
check_manifold_point(::Oblique, ::Any)
function check_manifold_point(M::Oblique{Ns,Ms}, p; kwargs...) where {Ns,Ms}
    if size(p) != (Ns + 1, Ms)
        return DomainError(
            length(p),
            "The matrix in `p` ($(size(p))) does not match the dimension of Oblique $((Ns,Ms)).",
        )
    end
    return check_manifold_point(PowerManifold(M.manifold, Ms), p; kwargs...)
end
@doc raw"""
    check_tangent_vector(M::Oblique p, X; check_base_point = true, kwargs...)

Checks whether `X` is a valid tangent vector to `p` on the [`Oblique`](@ref) `M`.
This means, that `p` is valid, that `X` is of correct dimension and columnswise
a tangent vector to the columns of `p` on the [`Sphere`](@ref).
The optional parameter `check_base_point` indicates, whether to call [`check_manifold_point`](@ref)  for `p`.
"""
function check_tangent_vector(
    M::Oblique{Ns,Ms},
    p,
    X;
    check_base_point = true,
    kwargs...,
) where {Ns,Ms}
    if check_base_point && size(p) != (Ns + 1, Ms)
        return DomainError(
            length(p),
            "The matrix `p` ($(size(p))) does not match the dimension of Oblique $((Ns,Ms)).",
        )
    end
    if size(X) != (Ns + 1, Ms)
        return DomainError(
            length(X),
            "The matrix `X` ($(size(X))) does not match the dimension of Oblique $((Ns,Ms)).",
        )
    end
    return check_tangent_vector(
        PowerManifold(M.manifold, Ms),
        p,
        X;
        check_base_point = check_base_point,
        kwargs...,
    )
end

get_iterator(M::Oblique{Ns,Ms}) where {Ns,Ms} = 1:Ms

@generated manifold_dimension(::Oblique{N,M}) where {N,M} = (N) * M

@generated representation_size(::Oblique{N,M}) where {N,M} = (N + 1, M)

show(io::IO, ::Oblique{N,M}) where {N,M} = print(io, "Oblique($(N+1),$(M))")


text:
module Chess

include("magic.jl")
include("piece.jl")
include("square.jl")
include("squareset.jl")
include("move.jl")
include("chess960.jl")
include("board.jl")
include("game.jl")
include("san.jl")
include("pgn.jl")
include("uci.jl")
include("book.jl")
include("mime.jl")

end # module


text:
"""
Pollard's rho algorithm (https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm)
"""
module PollardRho

#=
According to <https://stackoverflow.com/a/48208795>, Pollard's rho algorithm does not always find
the factors of perfect powers.
=#

function squared_plus_1(n::Integer)
    function f(a::Integer)::Integer
        a2 = powermod(a, 2, n)
        mod(a2 + 1, n)
    end
end

function pollard_rho_factor(n::Integer, f::Function=squared_plus_1(n))::Integer
    # Brent's cycle detection algorithm
    l = 0
    tortoise = hare = 2  # the smallest prime
    cycle_len = 2  # prime gap
    factor = 1
    while factor == 1
        if l == cycle_len
            tortoise = hare
            cycle_len *= 2
            l = 0
        end
        hare = f(hare)
        l += 1
        factor = gcd(hare - tortoise, n)
    end
    factor
end

# Unit tests
using Test

function check_prime(n::Integer)::Bool
    factor = pollard_rho_factor(n)
    n == factor
end

function check_composite(n::Integer)::Bool
    factor = pollard_rho_factor(n)
    1 < factor < n && n % factor == 0
end

@testset "PollardRho" begin
    prime_list = [7, 4567, 15960005471, 1238926361552897]
    composite_list = [
        101 * 103, 1009 * 1013, 10007 * 10009,
        1234567,
        12389263_66_1552897,
    ]
    @test check_prime.(prime_list) |> all
    @test check_composite.(composite_list) |> all
end

end  # module

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
using Printf

function multifact(n::Integer, k::Integer)
    n > 0 && k > 0 || throw(DomainError())
    k > 1 || factorial(n)
    return prod(n:-k:2)
end

const khi = 5
const nhi = 10
println("Showing multifactorial for n in [1, $nhi] and k in [1, $khi].")
for k = 1:khi
    a = multifact.(1:nhi, k)
    lab = "n" * "!" ^ k
    @printf("  %-6s →  %s\n", lab, a)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using JuMP, Gurobi, LinearAlgebra, Suppressor, StatsBase, Compat, Mosek, DataFrames, CSV, Random
include("convex_penalties_edm.jl")
include("exact_edm.jl")

results_template = DataFrame(n=Int[], r=Real[], t_sdp=Real[], err_sdpdtr=Real[], nocorr_sdp=Int[], t_ex=Real[], err_exdtr=Real[], nocorr_ex=Int[], MIPGap=Real[], nodeCountEX=Int[], t_oa=Real[], err_oadtr=Real[], nocorr_oa=Int[], MIPGapOA=Real[], nodeCountOA=Int[], cutCountOA=Int[])

results_edmradio_ex=similar(results_template, 0)

# Gram matrix recovery test problem
r=2
d_radio=0.3 # You will need to vary this to reproduce all of the results in Table 6
n=30  # You will need to vary this to reproduce all of the results in Table 6
lambda=1.0*n^2
gamma=1.0/n

for ARG in ARGS #1-4
    d_radio=[0.1, 0.2, 0.3, 0.4][parse(Int, ARG)]

    xs=rand(n, r).-0.5 #generate n random points in r dimensional space uniformly distributed on [-0.5, 0.5]^r
    D=zeros(n,n) #the matrix we would like to recover.
    for i=1:n #in a loop is slightly slower, but unimportant
        for j=1:(i-1)
            D[i,j]=max(0.0, norm(xs[i,:].-xs[j,:])^2+0.01*(randn(1)[1])) #Corrupt the distance with some noise, in a symmetric fashion
            D[j,i]=D[i,j]
        end
    end

    # Get observed set of coordinates
    coords=[repeat(1:n, inner=size(1:n,1)) repeat(1:n, outer=size(1:n,1))];
    entries=Int[]
    for i=1:n^2
        if norm(D[coords[i,1], coords[i,2]])<=d_radio^2
            push!(entries, i) #Remark: D_ij=D_ji (including noise), so this ensures that if D_ij is in the radio range, so is D_ji; no need to account for symmetry directly
        end
    end

    numEntriesSampled=size(entries,1)

    # Solve the sdp relaxation
    t_sdp=@elapsed G_imputed, Y_relax=getMatrixTraceNorm_EDM_Gram(coords, entries, D, n, numEntriesSampled, d_radio, lambda, gamma, r);
    u,v,=svd(Y_relax)
    Y_ws=u[:,1:r]*u[:,1:r]'

    D_imputed=diag(G_imputed)*ones(n)'+ones(n)*diag(G_imputed)'.-2*G_imputed
    @show err_sdpDTR=norm(D_imputed.-D)/norm(D)
    @show sdp_corr=sum(abs.(D.-D_imputed).<1e-2)
    @show sdp_incorr =sum(abs.(D.-D_imputed).>=1e-2)

    # Commented out as the submitted version of the paper only includes multi tree
    t_oa=@elapsed G_imputed,nodecountOA, cutCountOA, miogapOA,=getMatrixFrobeniusNorm_cuttingplanesmultitree_exact_inout(coords, entries, D, n, numEntriesSampled, lambda, d_radio, r, gamma, Y_ws, Y_relax);
    D_imputed=diag(G_imputed)*ones(n)'+ones(n)*diag(G_imputed)'.-2*G_imputed
    @show err_oaDTR=norm(D_imputed.-D)/norm(D)
    @show oa_corr=sum(abs.(D.-D_imputed).<1e-2)
    @show oa_incorr =sum(abs.(D.-D_imputed).>=1e-2)

    # Commented out as the submitted version of the paper only includes multi tree
    t_oa=-1.0#@elapsed G_imputed,nodecountOA, cutCountOA, miogapOA,=getMatrixFrobeniusNorm_cuttingplanes_exact_inout(coords, entries, D, n, numEntriesSampled, lambda, d_radio, r, gamma, Y_ws, Y_relax);
    # D_imputed=diag(G_imputed)*ones(n)'+ones(n)*diag(G_imputed)'.-2*G_imputed
    @show err_oaDTR=-1.0#norm(D_imputed.-D)/norm(D)
    @show oa_corr=-1.0#sum(abs.(D.-D_imputed).<1e-2)
    @show oa_incorr =-1.0#sum(abs.(D.-D_imputed).>=1e-2)
    miogapOA=-1.0; nodeCountOA=-1; cutCountOA=-1

    t_ex=-1.0@elapsed G_imputed,miogap,nodeCountEX=getMatrixExact_EDM_Gram(coords, entries, D, n, r, numEntriesSampled, d_radio, lambda, gamma, D_imputed);
    # D_imputed=diag(G_imputed)*ones(n)'+ones(n)*diag(G_imputed)'.-2*G_imputed
    @show err_exDTR=-1.0#norm(D_imputed.-D)/norm(D)
    @show ex_corr=-1.0#sum(abs.(D.-D_imputed).<1e-2)
    @show ex_incorr =-1.0#sum(abs.(D.-D_imputed).>=1e-2)
    miogap=-1.0; nodeCountEX=-1

    push!(results_edmradio_ex, [n, r, t_sdp, err_sdpDTR, sdp_corr, t_ex, err_exDTR, ex_corr, miogap, nodeCountEX, t_oa, err_oaDTR, oa_corr, miogapOA, nodecountOA, cutCountOA])
end
CSV.write("oavsgurobi_edmgram.csv", results_edmradio_ex, append=true)

max_stars_repo_path: Euclidean Distance Matrix Completion/createTable6.jl

max_stars_repo_name: ryancorywright/MixedProjectionSoftware

max_issues_repo_name: ryancorywright/MixedProjectionSoftware

max_forks_repo_name: ryancorywright/MixedProjectionSoftware


text:
include(joinpath(@__DIR__, "pedagogy_graph.jl"))
using Random

dims = [15, 15] # if > (14, 14) the display becomes all black, don't know why !
periodic = false
H, J = 0.0, 0.01

ising = PRS.Ising(dims, J, H; periodic=periodic)

rng = Random.MersenneTwister(123)
@time sample = PRS.generate_sample_gibbs_perfect(rng, ising)

path = joinpath(@__DIR__, "output", "ising", "ising_gibbs_perfect.pdf")
p = plot(ising, sample, dims, path; nodelabel=nothing)


text:
function codegen(::Type{Val{:iterate}}, ::Type{RAM}, job::BasicMCJob)
  result::Expr
  update = []
  noupdate = []
  body = []
  dindex::Integer

  vform = variate_form(job.pstate)
  if vform != Univariate && vform != Multivariate
    error("Only univariate or multivariate parameter states allowed in MH code generation")
  end

  push!(body, :(_job.sstate.count += 1))

  if job.tuner.verbose
    push!(body, :(_job.sstate.tune.proposed += 1))
  end

  if vform == Univariate
    push!(body, :(_job.sstate.randnsample = randn()))
    push!(body, :(_job.sstate.pstate.value = _job.pstate.value+_job.sstate.S*_job.sstate.randnsample))
  elseif vform == Multivariate
    push!(body, :(_job.sstate.randnsample[:] = randn(_job.pstate.size)))
    push!(body, :(_job.sstate.pstate.value[:] = _job.pstate.value+_job.sstate.S*_job.sstate.randnsample))
  end

  push!(body, :(_job.parameter.logtarget!(_job.sstate.pstate)))

  push!(body, :(_job.sstate.ratio = _job.sstate.pstate.logtarget-_job.pstate.logtarget))

  if vform == Univariate
    push!(update, :(_job.pstate.value = _job.sstate.pstate.value))
  elseif vform == Multivariate
    push!(update, :(_job.pstate.value = copy(_job.sstate.pstate.value)))
  end
  push!(update, :(_job.pstate.logtarget = _job.sstate.pstate.logtarget))
  if in(:loglikelihood, job.outopts[:monitor]) && job.parameter.loglikelihood! != nothing
    push!(update, :(_job.pstate.loglikelihood = _job.sstate.pstate.loglikelihood))
  end
  if in(:logprior, job.outopts[:monitor]) && job.parameter.logprior! != nothing
    push!(update, :(_job.pstate.logprior = _job.sstate.pstate.logprior))
  end
  dindex = findfirst(job.outopts[:diagnostics], :accept)
  if dindex != 0
    push!(update, :(_job.pstate.diagnosticvalues[$dindex] = true))
    push!(noupdate, :(_job.pstate.diagnosticvalues[$dindex] = false))
  end
  if job.tuner.verbose
    push!(update, :(_job.sstate.tune.accepted += 1))
  end

  push!(body, Expr(:if, :(_job.sstate.ratio > 0 || (_job.sstate.ratio > log(rand()))), Expr(:block, update...), noupdate...))

  if job.tuner.verbose
    fmt_iter = format_iteration(ndigits(job.range.burnin))
    fmt_perc = format_percentage()

    push!(body, :(
      if _job.sstate.tune.totproposed <= _job.range.burnin && mod(_job.sstate.tune.proposed, _job.tuner.period) == 0
        rate!(_job.sstate.tune)
        println(
          "Burnin iteration ",
          $(fmt_iter)(_job.sstate.tune.totproposed),
          " of ",
          _job.range.burnin,
          ": ",
          $(fmt_perc)(100*_job.sstate.tune.rate),
          " % acceptance rate"
        )
        reset_burnin!(_job.sstate.tune)
      end
    ))
  end

  if vform == Univariate
    push!(body, :(_job.sstate.η = min(1, _job.sstate.count^(-_job.sampler.γ))))

    push!(body, :(_job.sstate.SST = _job.sstate.η*(min(1, exp(_job.sstate.ratio))-_job.sampler.targetrate)))

    push!(body, :(_job.sstate.SST = abs2(_job.sstate.S)*(1+_job.sstate.SST)))

    push!(body, :(_job.sstate.S = chol(_job.sstate.SST)))
  elseif vform == Multivariate
    push!(body, :(_job.sstate.η = min(1, _job.pstate.size*_job.sstate.count^(-_job.sampler.γ))))

    push!(
      body,
      :(
        _job.sstate.SST[:, :] = (
          _job.sstate.randnsample*_job.sstate.randnsample'/dot(_job.sstate.randnsample, _job.sstate.randnsample)*
          _job.sstate.η*(min(1, exp(_job.sstate.ratio))-_job.sampler.targetrate)
        )
      )
    )

    push!(body, :(_job.sstate.SST[:, :] = _job.sstate.S*(eye(_job.pstate.size)+_job.sstate.SST)*_job.sstate.S'))

    push!(body, :(_job.sstate.S[:, :] = ctranspose(chol(_job.sstate.SST))))
  end

  if !job.plain
    push!(body, :(produce()))
  end

  @gensym _iterate

  result = quote
    function $_iterate(_job::BasicMCJob)
      $(body...)
    end
  end

  result
end


text:
hassym(lib, sym) = Libdl.dlsym_e(lib, sym) != C_NULL

# call dlsym_e on a sequence of symbols and return the symbol that gives
# the first non-null result
function findsym(lib, syms...)
    for sym in syms
        if hassym(lib, sym)
            return sym
        end
    end
    error("no symbol found from: ", syms)
end

# Need to set PythonHome before calling GetVersion to avoid warning (#299).
# Unfortunately, this poses something of a chicken-and-egg problem because
# we need to know the Python version to set PythonHome via the API.  Note
# that the string (or array) passed to Py_SetPythonHome needs to be a
# constant that lasts for the lifetime of the program, which is why we
# can't use Cwstring here (since that creates a temporary copy).
function Py_SetPythonHome(libpy, PYTHONHOME, wPYTHONHOME, pyversion)
    if !isempty(PYTHONHOME)
        if pyversion.major < 3
            ccall(Libdl.dlsym(libpy, :Py_SetPythonHome), Void, (Cstring,), PYTHONHOME)
        else
            ccall(Libdl.dlsym(libpy, :Py_SetPythonHome), Void, (Ptr{Cwchar_t},), wPYTHONHOME)
        end
    end
end

# need to be able to get the version before Python is initialized
Py_GetVersion(libpy) = unsafe_string(ccall(Libdl.dlsym(libpy, :Py_GetVersion), Ptr{UInt8}, ()))


text:
import BenchmarkTools as BT
import DataStructures as DS
import Dates
import GridWorlds as GW
import ReinforcementLearningBase as RLBase

function run_random_policy!(env, num_resets, steps_per_reset)
    for reset_number in 1 : num_resets
        RLBase.reset!(env)
        for step_number in 1 : steps_per_reset
            state = RLBase.state(env)
            action = rand(RLBase.action_space(env))
            env(action)
            is_terminated = RLBase.is_terminated(env)
            reward = RLBase.reward(env)
        end
    end

    return nothing
end

function benchmark(Env, num_resets, steps_per_reset)
    benchmark = DS.OrderedDict()

    env = GW.RLBaseEnv(Env())

    benchmark[:random_policy] = BT.@benchmark run_random_policy!($(Ref(env))[], $(Ref(num_resets))[], $(Ref(steps_per_reset))[])
    benchmark[:reset] = BT.@benchmark RLBase.reset!($(Ref(env))[])
    benchmark[:state] = BT.@benchmark RLBase.state($(Ref(env))[])

    action_names = GW.get_action_names(env.env)

    for action in RLBase.action_space(env)
        benchmark[action_names[action]] = BT.@benchmark $(Ref(env))[]($(Ref(action))[])
    end

    benchmark[:is_terminated] = BT.@benchmark RLBase.is_terminated($(Ref(env))[])
    benchmark[:reward] = BT.@benchmark RLBase.reward($(Ref(env))[])

    @info "$(nameof(Env)) benchmark complete"

    return benchmark
end

function benchmark(Envs::Vector, num_resets, steps_per_reset)
    benchmarks = DS.OrderedDict()

    for Env in Envs
        benchmarks[nameof(Env)] = benchmark(Env, num_resets, steps_per_reset)
    end

    @info "Benchmarks complete"

    return benchmarks
end

function get_summary(trial::BT.Trial)
    median_trial = BT.median(trial)
    memory = BT.prettymemory(median_trial.memory)
    median_time = BT.prettytime(median_trial.time)
    return memory, median_time
end

function get_table(benchmark)
    title = "|"
    separator = "|"
    data = "|"

    for key in keys(benchmark)
        title = title * String(key) * "|"
        separator = separator * ":---:|"
        memory, median_time = get_summary(benchmark[key])
        data = data * "$(memory)<br>$(median_time)|"
    end

    return title, separator, data
end

function generate_benchmark_file(; Envs = GW.ENVS, num_resets = 100, steps_per_reset = 100, file_name = nothing)
    date = Dates.format(Dates.now(), "yyyy_mm_dd_HH_MM_SS")

    if isnothing(file_name)
        file_name = date * ".md"
    end

    io = open(file_name, "w")

    benchmarks = benchmark(Envs, num_resets, steps_per_reset)

    println(io, "Date: $(date) (yyyy_mm_dd_HH_MM_SS)")
    println(io)
    println(io, "**Note:** The time in benchmarks is the median time. While benchmarking the random policy, each environment was run for $(num_resets) episodes of $(steps_per_reset) steps each, that is, a total of $(num_resets * steps_per_reset) steps.")
    println(io)
    println(io, "## List of Environments")

    for Env in Envs
        name_string = String(nameof(Env))
        println(io, "  1. [$(name_string)](#$(lowercase(name_string)))")
    end

    println(io)

    for key in keys(benchmarks)
        println(io, "### " * String(key))
        title, separator, data = get_table(benchmarks[key])
        println(io, title)
        println(io, separator)
        println(io, data)
        println(io)
    end

    close(io)

    return nothing
end


text:
using DualNumbers
using Interpolations
using Test

@testset "Extrapolation" begin

    f(x) = sin((x-3)*2pi/9 - 1)
    xmax = 10
    A = Float64[f(x) for x in 1:xmax]

    itpg = interpolate(A, BSpline(Linear()))

    etpg = extrapolate(itpg, Flat())
    @test typeof(etpg) <: AbstractExtrapolation

    @test etpg(-3) == etpg(-4.5) == etpg(0.9) == etpg(1.0) == A[1]
    @test etpg(10.1) == etpg(11) == etpg(148.298452) == A[end]

    etpf = @inferred(extrapolate(itpg, NaN))
    @test typeof(etpf) <: Interpolations.FilledExtrapolation
    @test parent(etpf) === itpg

    @test @inferred(size(etpf)) == (xmax,)
    @test isnan(@inferred(etpf(-2.5)))
    @test isnan(etpf(0.999))
    @test @inferred(etpf(1)) ≈ f(1)
    @test etpf(10) ≈ f(10)
    @test isnan(@inferred(etpf(10.001)))

    @test etpf(2.5,1) == etpf(2.5)   # for show method
    @test_throws BoundsError etpf(2.5,2)
    @test_throws BoundsError etpf(2.5,2,1)

    x =  @inferred(etpf(dual(-2.5,1)))
    @test isa(x, Dual)

    etpl = extrapolate(itpg, Line())
    k_lo = A[2] - A[1]
    x_lo = -3.2
    @test etpl(x_lo) ≈ A[1] + k_lo * (x_lo - 1)
    k_hi = A[end] - A[end-1]
    x_hi = xmax + 5.7
    @test etpl(x_hi) ≈ A[end] + k_hi * (x_hi - xmax)

    xmax, ymax = 8,8
    g(x, y) = (x^2 + 3x - 8) * (-2y^2 + y + 1)

    itp2g = interpolate(Float64[g(x,y) for x in 1:xmax, y in 1:ymax], (BSpline(Quadratic(Free(OnGrid()))), BSpline(Linear())))
    etp2g = extrapolate(itp2g, (Line(), Flat()))

    @test @inferred(etp2g(-0.5,4)) ≈ itp2g(1,4) - 1.5 * epsilon(etp2g(dual(1,1),4))
    @test @inferred(etp2g(5,100)) ≈ itp2g(5,ymax)

    etp2ud = extrapolate(itp2g, ((Line(), Flat()), Flat()))
    @test @inferred(etp2ud(-0.5,4)) ≈ itp2g(1,4) - 1.5 * epsilon(etp2g(dual(1,1),4))
    @test @inferred(etp2ud(5, -4)) == etp2ud(5,1)
    @test @inferred(etp2ud(100, 4)) == etp2ud(8,4)
    @test @inferred(etp2ud(-.5, 100)) == itp2g(1,8) - 1.5 * epsilon(etp2g(dual(1,1),8))

    etp2ll = extrapolate(itp2g, Line())
    @test @inferred(etp2ll(-0.5,100)) ≈ (itp2g(1,8) - 1.5 * epsilon(etp2ll(dual(1,1),8))) + (100 - 8) * epsilon(etp2ll(1,dual(8,1)))

    # Allow element types that don't support conversion to Int (#87):
    etp87g = extrapolate(interpolate([1.0im, 2.0im, 3.0im], BSpline(Linear())), 0.0im)
    @test @inferred(etp87g(1)) == 1.0im
    @test @inferred(etp87g(1.5)) == 1.5im
    @test @inferred(etp87g(0.75)) == 0.0im
    @test @inferred(etp87g(3.25)) == 0.0im

    # Make sure it works with Gridded too
    etp100g = extrapolate(interpolate(([10;20],),[100;110], Gridded(Linear())), Flat())
    @test @inferred(etp100g(5)) == 100
    @test @inferred(etp100g(15)) == 105
    @test @inferred(etp100g(25)) == 110
    # issue #178
    a = randn(10,10) + im*rand(10,10)
    etp = @inferred(extrapolate(interpolate((1:10, 1:10), a, Gridded(Linear())), 0.0))
    @test @inferred(etp(-1,0)) === 0.0+0.0im

    # check all extrapolations work with vectorized indexing
    for E in [0,Flat(),Line(),Periodic(),Reflect()]
        @test (@inferred(extrapolate(interpolate([0,0],BSpline(Linear())),E))([1.2, 1.8, 3.1])) == [0,0,0]
    end

    # Issue #156
    F     = *(collect(1.0:10.0), collect(1:4)')
    itp   = interpolate(F, (BSpline(Linear()), NoInterp()));
    itps   = scale(itp, 1:10, 1:4)
    itpe   = extrapolate(itps, (Line(), Throw()))
    @test itpe(10.1, 1) ≈ 10.1
    @test_throws BoundsError itpe(9.9, 0)

    # Issue #232
    targeterr = ArgumentError("cannot create a filled extrapolation with a type Line, use a value of this type instead (e.g., Line())")
    @test_throws targeterr extrapolate(itp, Line)

    # Issue #244
    xs = range(1e-2, stop = 8.3, length = 3)
    ys = sort(rand(3))
    itp = LinearInterpolation(xs, ys, extrapolation_bc = Flat())
    @test itp(8.3) ≈ ys[end]

    # Issue #288
    etp = extrapolate(interpolate(collect(1:12), NoInterp()) , 0)
    @test etp(0:23) == [0:12; zeros(Int, 11)]

    include("type-stability.jl")
    include("non1.jl")
end


text:
# Belief Propagation                #
# --------------------------------- #
@rule NormalMeanVariance(:v, Marginalisation) (m_out::PointMass, m_μ::UnivariateNormalDistributionsFamily) = begin
    m_out_mean        = mean(m_out)
    m_μ_mean, m_μ_cov = mean_cov(m_μ)

    return ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(m_μ_cov + x) / 2 - 1 / (2x) * (m_out_mean - m_μ_mean) ^ 2)
end

@rule NormalMeanVariance(:v, Marginalisation) (m_out::UnivariateNormalDistributionsFamily, m_μ::PointMass) = begin
    m_out_mean, m_out_cov = mean_cov(m_out)
    m_μ_mean              = mean(m_μ)

    return ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(m_out_cov + x) / 2 - 1 / (2x) * (m_out_mean - m_μ_mean) ^ 2)
end

@rule NormalMeanVariance(:v, Marginalisation) (m_out::UnivariateNormalDistributionsFamily, m_μ::UnivariateNormalDistributionsFamily) = begin
    m_out_mean, m_out_cov = mean_cov(m_out)
    m_μ_mean, m_μ_cov     = mean_cov(m_μ)

    return ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(m_out_cov + m_μ_cov + x) / 2 - 1 / (2x) * (m_out_mean - m_μ_mean) ^ 2)
end


text:
module CodecBGZF

export
    BGZFCompressorStream,
    BGZFDecompressorStream,
    BGZFError,
    VirtualOffset,
    offsets,
    gzi

using LibDeflate
import LibDeflate: GzipExtraField
using TranscodingStreams

import TranscodingStreams:
    TranscodingStreams,
    TranscodingStream,
    Memory,
    Error

import Base.Threads.@spawn

const NO_GZIP_EXTRA_FIELD_BSIZE = "No GZIP extra field \"BSIZE\""

const DE_COMPRESSOR = Union{Compressor, Decompressor}

"""
    BGZFError(message::String)

BGZF de/compressor errored with `message` when processing data."
"""
struct BGZFError <: Exception
    message::String
end

function bitload(T::Type{<:Base.BitInteger}, data::Vector{UInt8}, p::Integer)
    ltoh(unsafe_load(Ptr{T}(pointer(data, p))))
end

function bitstore(v::Base.BitInteger, data::Vector{UInt8}, p::Integer)
    unsafe_store!(Ptr{typeof(v)}(pointer(data, p)), htol(v))
end

@noinline bgzferror(s::String) = throw(BGZFError(s))

include("virtualoffset.jl")
include("block.jl")
include("bgzfstream.jl")

end # module


text:
using LinearAlgebra
include("../ellell.jl")
include("../expmap.jl")

x1 = [0.; 0.; 0.];
a1 = 1.;
b1 = 1.;
c1 = 1.;
R1 = Matrix{Float64}(I, 3, 3)

x2 = [0.; 0.; 2.0];
a2 = 1.;
b2 = 1.;
c2 = 1.1;
R2 = Matrix{Float64}(I, 3, 3)

(conpnt, normal, depth, A, B) = ellipsoid_ellipsoid_contact(x1, a1, b1, c1, R1, x2, a2, b2, c2, R2)

println("-------------------------------------------------------------")
println("for x1 = ", x1, "; a1, b1, c1 = ", a1, ", ", b1, ", ", c1, "; R1 = ", R1)
println("and x2 = ", x2, "; a2, b2, c2 = ", a2, ", ", b2, ", ", c2, "; R2 = ", R2)
println("conpnt = ", conpnt, ", normal = ", normal, ", depth = ", depth, ", A, B = ", A, ", ", B)

max_stars_repo_name: tkoziara/edyna

max_issues_repo_name: tkoziara/edyna

max_forks_repo_name: tkoziara/edyna


text:
@testset "Differential evolution optimizer" begin

ss = symmetric_search_space(1, (0.0, 10.0))
fake_problem = FunctionBasedProblem(x -> 0.0, "test_problem", MinimizingFitnessScheme, ss)
DE = de_rand_1_bin(fake_problem, ParamsDict(
    :Population => reshape(collect(1.0:10.0), (1, 10)),
    :f => 0.4, :cr => 0.5, :NumParents => 3))

@testset "ask()/tell!()" begin
    for i in 1:NumTestRepetitions
        res = BlackBoxOptim.ask(DE)
        @test BlackBoxOptim.candi_pool_size(BlackBoxOptim.population(DE)) == 0 # no candidates in the pool as we just exhausted it
        @test length(res) == 2

        trial, target = res

        @test ndims(trial.params) == 1
        @test (1 <= trial.index <= popsize(DE.population))
        @test in(trial.params, DE.embed.searchSpace)

        @test ndims(target.params) == 1
        @test (1 <= target.index <= popsize(DE.population))
        @test in(target.params, DE.embed.searchSpace)

        @test trial.index == target.index

        BlackBoxOptim.tell!(DE, res)
        @test BlackBoxOptim.candi_pool_size(BlackBoxOptim.population(DE)) == 2 # test that all candidates returned to the pool
    end
end

end


text:
using SciMLBase

@testset "DiffEqOperator" begin
    A = rand(10,10);
    @test eachindex(A) === eachindex(SciMLBase.DiffEqArrayOperator(A))
    @test eachindex(A') === eachindex(SciMLBase.DiffEqArrayOperator(A'))
end


text:
@testset "Creating a canvas" begin
    canvas = Canvas{Int}(10, 20)
    @test size(canvas, 1) == 10
    @test size(canvas, 2) == 20
    @test iszero(canvas.pixels)
end
@testset "Writing pixels to a canvas" begin
    c = Canvas{Int}(10, 20)
    red = Color(1, 0, 0)
    c[2, 3] = red
    @test c[2, 3] == red
end


text:
include(joinpath(dirname(@__FILE__),"..", "src/WebBLAST/", "hit.jl"))
using BioSeq

# check_threshold tests

test_hit = Hit(1,"gi|18412314|ref|NP_565204.1|","tTest Hit","NP_565204",166,
	[Hsp(1,132.88,2.93692e-37,1,62,0,0,62,62,0,62,
	aminoacid("NQLQQLQNPGESPPVHPFVAPLSYLLGTWRGQGEGEYPTIPSFRYGEEIRFSHSGKPVIAY"),
	"MNQLQQLQNPGESPPVHPFVAPLSYLLGTWRGQGEGEYPTIPSFRYGEEIRFSHSGKPVIAY",
	aminoacid("MNQLQQLQNPGESPPVHPFVAPLSYLLGTWRGQGEGEYPTIPSFRYGEEIRFSHSGKPVIAY"),
	aminoacid("MNQLQQLQNPGESPPVHPFVAPLSYLLGTWRGQGEGEYPTIPSFRYGEEIRFSHSGKPVIAY"))])

@test check_threshold(0.01, test_hit) == true
@test check_threshold(2.9e-38, test_hit) == false

# fastarepresentation test

@test fastarepresentation(test_hit) == ("gi|18412314|ref|NP_565204.1|","MNQLQQLQNPGESPPVHPFVAPLSYLLGTWRGQGEGEYPTIPSFRYGEEIRFSHSGKPVIAY")


text:
nt = 1000
seed = 99
ssm = acqr(0.8,1.5,2,3,0.5,2)
m = acqr_simula(ssm,nt,seed)
y = m["y"]
a = ssm.a
c = ssm.c
q = ssm.q
r = ssm.r
x10 = ssm.m1
P10 = ssm.P1

###########################
print("Testing acqr_kalman: ")
start = acqr(0.5,0.5,1,1,0,0)
me = acqr_em_s(y,start,max_iter=10)
mk = acqr_kalman(y,me.fit,s=true)
if sum(mk.prediction.^2) > 1e-6
    println("OK")
end


text:
include("run.jl")
include(joinpath(SLIDE_DIR,"src","build","share_cfs.jl"))
include(joinpath(SLIDE_DIR,"src","build","share_gsp.jl"))
include(joinpath(SLIDE_DIR,"src","build","share_pce.jl"))
include(joinpath(SLIDE_DIR,"src","build","share_sgf.jl"))
include(joinpath(SLIDE_DIR,"src","build","share_utd.jl"))

bshr_int = read_to_check(joinpath(SLIDE_DIR,"data","readfiles","benchmark","shr","shr_int.yml"));
bshr_out = read_to_check(joinpath(SLIDE_DIR,"data","readfiles","benchmark","shr","shr_out.yml"));
bshr = merge(bshr_int, bshr_out)

# files_share = write_yaml(READ_DIR, XLSXInput("generate_yaml.xlsx", "share", "B1:G150", "share"))
# y = [read_file(files_share[ii]) for ii in 1:length(files_share)]
# files_share = run_yaml(ensurearray(files_share))
# shr_read = Dict(Symbol(y[ii]["PathOut"][end][1:end-4]) =>
#     read_file(joinpath(y[ii]["PathOut"]...)) for ii in 1:length(y))

# Filter data and extrapolate values as appropriate.
shr = copy(shr_read)
shr = Dict(k => sort(filter_with(df, set; extrapolate = true)) for (k, df) in shr)
shr[:va0] = edit_with(io[:va0], Rename(:j,:s))

share_pce!(shr)
share_sgf!(shr)
share_utd!(shr, set)
share_region!(shr, set)
share_labor!(shr, set)
share_rpc!(shr, set)

# Read files to benchmark.
shr_bench = Dict()
[benchmark!(shr_bench, k, bshr_out, shr) for k in intersect(keys(bshr_out),keys(shr))]

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
# This file was generated, do not modify it. # hide
coerce!(df, :zipcode => Multiclass)
df.isrenovated  = @. !iszero(df.yr_renovated)
df.has_basement = @. !iszero(df.sqft_basement)
schema(df)


text:
using Documenter, GridapGeosciences

makedocs(;
    modules=[GridapGeosciences],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/santiagobadia/GridapGeosciences.jl/blob/{commit}{path}#L{line}",
    sitename="GridapGeosciences.jl",
    authors="Santiago Badia",
    assets=String[],
)

deploydocs(;
    repo="github.com/santiagobadia/GridapGeosciences.jl",
)

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
#=ДАНО: На ограниченном внешней прямоугольной рамкой поле имеется ровно одна внутренняя перегородка в форме прямоугольника.
 Робот - в произвольной клетке поля между внешней и внутренней перегородками.

 РЕЗУЛЬТАТ: Робот - в исходном положении и по всему периметру внутренней перегородки поставлены маркеры.=#
 function mark_between_border(r::Robot)
    num_steps_1=moves!(r,Sud)
    num_steps_2=moves!(r,West)
    num_steps_3=moves!(r,Sud)
#учитывается возможность изменения траектории из-за внутренней прямоугольной рамки

    find_border!(r,Ost,Nord)
    putmarker!(r)
    for side in (Nord,Ost,Sud,West)
        while isborder(r,HorizonSide(mod(Int(side)-1,4)))==true
            move!(r,side)
            putmarker!(r)
        end
        move!(r,HorizonSide(mod(Int(side)-1,4)))
        putmarker!(r)
    end
 #проставляем маркеры по периметру прямоугольной рамки

    back(r,Sud)
    back(r,West)
#возвращает в угол

    moves!(r,Nord,num_steps_3)
    moves!(r,Ost,num_steps_2)
    moves!(r,Nord,num_steps_1)
#возвращает в исходное положение
end
function find_border!(r::Robot,side_to_border::HorizonSide, side_of_move::HorizonSide)
    while isborder(r,side_to_border)==false
        if isborder(r,side_of_move)==false
            move!(r,side_of_move)
        else
            move!(r,side_to_border)
            side_of_move=inverse(side_of_move)
        end
    end
end
#=функция,при которой находится перегородка;входныепеременные:
робот,по какой линии смотреть(вертикальной или горизонтальной),какие перегородки ищем=#

function moves!(r::Robot, side::HorizonSide)
    num_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_steps+=1
    end
    return num_steps
end

function moves!(r::Robot,side::HorizonSide,num_steps::Int)
    for _ in 1:num_steps
        move!(r,side)
    end
end

function back(r::Robot,side::HorizonSide)
    while isborder(r,side)==false
        move!(r,side)
    end
end
#просто путь до границы

inverse(side::HorizonSide) = HorizonSide(mod(Int(side)+2, 4))

max_stars_repo_name: ilya-design/mirea-progs

max_issues_repo_name: ilya-design/mirea-progs

max_forks_repo_name: ilya-design/mirea-progs


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# ╔═╡ fab940e4-260e-11eb-02fb-017189d5dd7e
begin
	using Pkg; Pkg.activate(".") #sets the Julia package environment
	using CSV #For processing csv files
	using DataFrames
	using StatsBase
	using Query #Enables SQL type query of DataFrames
	using MLDataUtils #Tools for splitting data
	path = pwd();
	data_path = path*"/../data/cp_data_cleaned_jl.csv"

	#Set random seed for reproducibility
	using Random
	const rng_seed = 42
	Random.seed!(rng_seed)
end;

# ╔═╡ d5f6a828-260c-11eb-360a-7fdf63568dd6
md"
# Splitting data into the train/validation/test dataset

It is important to split your full dataset into train/validation/test datasets, and reliably use the same datasets for your modeling tasks later.

An example of using different splitting for train/validation/test is shown in the figure in the [python notebook](https://github.com/anthony-wang/BestPractices/blob/master/notebooks/2-data_splitting.ipynb)

Lets load the Julia packages we need and the preprocessed data generated from [1-data\_loading\_cleanup\_processing.ipynb]()
"

# ╔═╡ b038f004-260f-11eb-021d-cd976cca1f98
dataframe = CSV.File(data_path) |> DataFrame

# ╔═╡ e632ec3c-260f-11eb-1432-43dad5c863e0
md"
Now we need to split the data into `inputs` and `output`. Here the input is the chemical composition and temperature and the output is the heat capacity. The standard notation is to describe inputs as `features` or `descriptors` with variable `X` and output as the `target` with variable `y`.

Note that the use of `!` in the dataframe type indicates all rows.

"

# ╔═╡ 32c4846e-2610-11eb-0a57-5f7bdb29198f
begin
	X = dataframe[!,[:formula,:T]]
	y = dataframe[!,[:Cₚ]]
	md"Size of input: $(size(X))
	   Size of output: $(size(y))
	"
end

# ╔═╡ 8f072bd2-2612-11eb-383f-c1688c98898a
md"
The next step is to utilize standard ML tools that enable random shuffling of the dataset and then splitting it into training and testing. Its more common to split this into training, validation, and testing where the validation is used to help assist with model suitability and hyperparameter evaluation, but in materials science the datasets usually are too small that such splittign would limit valuable training data.

As mentioned in the orignal jupyter notebook, a rule of thumb one can go by is:

  --             | Training Split | Validation Split | Testing Split
:--------------  | :------------  | :--------------: | --------:
Portion to keep  | 50-70%         | 20-30%           | 10-20%

Note the random number generator used for MLDataUtils uses the `Random.GLOBAL_RNG` value which is set by the seed.
"

# ╔═╡ f7e7cba6-2611-11eb-0ddd-51add42ce94f
begin
	Xshuffle,yshuffle = shuffleobs((X,y))
	(Xtrain,ytrain),(Xtest,ytest) = splitobs((Xshuffle,yshuffle),at=0.80)
	md"Size of train input: $(size(Xtrain))
	   SIze of test input: $(size(Xtest))
	"
end

# ╔═╡ e92ea6d0-2635-11eb-2649-27d8ce85ea43
md"

One thing we didn't consider is the fact that the input and output features in the test set shouldn't also be in the training set. In other words, the data has many entries for the same chemical formula given that the data was taken at different temperatures.

Therefore we have to do some additional data processing to ensure that we properly manage this.
"

# ╔═╡ 4619c080-2638-11eb-048a-571e53d02fef
begin
	nrows = nrow(Xtrain);
	nunique_formulae = length(unique(Xtrain[:formula]));
end;

# ╔═╡ 4de8b85e-2636-11eb-2679-c9b1ea4fd78e
md"
There are a total of $(nrows) in the feature training set. But the number of unique formulae is $(nunique_formulae). Below is the unique formula and the counts for each.
"

# ╔═╡ 5f7430e4-2636-11eb-061b-79d0c6cb107d
counts_train = sort(combine(groupby(Xtrain,[:formula]),nrow=>:count))

# ╔═╡ 6651a89a-2638-11eb-204e-51f5c3182c50
counts_test = sort(combine(groupby(Xtest,[:formula]),nrow=>:count))

# ╔═╡ 9d0c5eca-2638-11eb-38be-1b343afbdb6e
md"
Compare the two unique count databases for train and test we can will see there is overlap which is not good practice.

What we need to do is get the list of unique chemical formula and then split the training and test datasets from that list.
"

# ╔═╡ 51aabf24-263a-11eb-0edc-a56f41717300
X_unique_formulae = unique(X[:formula]);

# ╔═╡ de263442-263a-11eb-110a-43b9ac45b2b0
md"

$(X_unique_formulae)

\

Now from the unique list we will split into training, validation, and test data. Theen we can use these seperate label arrays to split the actual dataframe. We will use the Query.jl package to allow for macros that can make constructing the different data sets easier.
"

# ╔═╡ b05f5720-2638-11eb-0259-eb0d51fad3c0
begin
	val_size = 0.20;
	test_size = 0.10;
	train_size = 1.00-val_size-test_size;
unique_train,unique_val,unique_test = splitobs(shuffleobs(X_unique_formulae),at=(train_size,val_size));
end;

# ╔═╡ 71020cfe-2646-11eb-0a9a-f1a5e3bed706
md"

The number of unique formulae for each dataset is:

Training: $(length(unique_train)) \
Validation: $(length(unique_val)) \
Testing: $(length(unique_test))

Note that this shuffling and splitting won't be the same each time which is different from the original jupyter implementation.
"

# ╔═╡ a69242c4-2643-11eb-36ca-275312dd2739
begin
	df_train = filter(row-> row[:formula] ∈ unique_train, dataframe)
	# Query.jl approach below
	#dataframe |> @filter(_.formula ∈ unique_test) |> DataFrame
	df_val = filter(row -> row[:formula] ∈ unique_val, dataframe)
	df_test = filter(row -> row[:formula] ∈ unique_test, dataframe)
end;

# ╔═╡ b10eb172-264a-11eb-054b-63a9436f86eb
md"

Below we implement the original manual approach from the jupyter notebook. The steps are:

1. Get number of entries based on fraction split.
2. Randomly sample from the list of unique formulae.
3. Update the list by removing those selected.
4. Slice the data from the training, validation, and test list.

"

# ╔═╡ 627b84b2-26cd-11eb-3289-85202be68813
begin
	Random.seed!(rng_seed)

	#Get the size of the split data
	all_formulae = copy(X_unique_formulae);
	num_unique = length(all_formulae);
	num_val_samples = round(Int,val_size*num_unique);
	num_test_samples = round(Int,test_size*num_unique);
	num_train_samples = round(Int,train_size*num_unique);

	#Select formula from unique list for each data set
	val_formulae = sample(all_formulae,num_val_samples,replace=false);
	all_formulae = [formula for formula in all_formulae if formula ∉ val_formulae];
	test_formulae = sample(all_formulae,num_test_samples,replace=false);
	all_formulae = [formula for formula in all_formulae if formula ∉ test_formulae];
	train_formulae = copy(all_formulae);

end	;

# ╔═╡ c1e61032-26d2-11eb-02b0-39feb496857c
begin
	df_train_org = filter(row-> row[:formula] ∈ train_formulae, dataframe)
	df_val_org = filter(row -> row[:formula] ∈ val_formulae, dataframe)
	df_test_org = filter(row -> row[:formula] ∈ test_formulae, dataframe)
end;

# ╔═╡ ebdccc4e-264a-11eb-2b6c-6de724bb411c
begin
	name_df_pair =  zip(("cp_train","cp_val","cp_test"),
				        (df_train_org,df_val_org,df_test_org));
	for (name,df) in name_df_pair
		datafile = path*"/../data/$(name)_jl.csv"
		df |> CSV.write(datafile);
	end
end;

# ╔═╡ Cell order:
# ╟─d5f6a828-260c-11eb-360a-7fdf63568dd6
# ╠═fab940e4-260e-11eb-02fb-017189d5dd7e
# ╠═b038f004-260f-11eb-021d-cd976cca1f98
# ╟─e632ec3c-260f-11eb-1432-43dad5c863e0
# ╠═32c4846e-2610-11eb-0a57-5f7bdb29198f
# ╟─8f072bd2-2612-11eb-383f-c1688c98898a
# ╠═f7e7cba6-2611-11eb-0ddd-51add42ce94f
# ╟─e92ea6d0-2635-11eb-2649-27d8ce85ea43
# ╠═4619c080-2638-11eb-048a-571e53d02fef
# ╟─4de8b85e-2636-11eb-2679-c9b1ea4fd78e
# ╠═5f7430e4-2636-11eb-061b-79d0c6cb107d
# ╠═6651a89a-2638-11eb-204e-51f5c3182c50
# ╟─9d0c5eca-2638-11eb-38be-1b343afbdb6e
# ╠═51aabf24-263a-11eb-0edc-a56f41717300
# ╟─de263442-263a-11eb-110a-43b9ac45b2b0
# ╠═b05f5720-2638-11eb-0259-eb0d51fad3c0
# ╟─71020cfe-2646-11eb-0a9a-f1a5e3bed706
# ╠═a69242c4-2643-11eb-36ca-275312dd2739
# ╟─b10eb172-264a-11eb-054b-63a9436f86eb
# ╠═627b84b2-26cd-11eb-3289-85202be68813
# ╠═c1e61032-26d2-11eb-02b0-39feb496857c
# ╠═ebdccc4e-264a-11eb-2b6c-6de724bb411c


text:
function Initi_model_struct(H::T,Hc::T,N::T,Wf::T,B::T,αf::T,β₁::T,β₂::T,
    Nₛ::T,rₘ::T,a_Jmax::T,b_Jmax::T,i::T,Kₓₗ₀::T) where {T<:Float64}
    cons = Constants()
    env = EnvironmentStruct()
    kinetic = PhotoKineticRates()
    photo = PhotoPar(kinetic,env.Tₐ)
    hydPar = HydraulicsPar(;Kₓₗ₀=Kₓₗ₀,i=i)
    treepar = TreePar(;αf=αf,β₁=β₁,β₂=β₂,rₘ=rₘ,Nₛ=Nₛ,a_Jmax=a_Jmax,b_Jmax=b_Jmax)

    Hs = H-Hc
    As = Wf/treepar.αf
    Ww = treepar.ρw*As*(treepar.β₁*H+treepar.β₂*Hs)
    treesize = TreeSize(Wf,Ww,H,Hs,As,B,N)

    model = CCPHStruct(cons,env,treepar,treesize,photo,hydPar)

    return model,kinetic
end

function Initi_model_struct(j::Integer,data::RoData,
    αf::T,β₁::T,β₂::T,Nₛ::T,rₘ::T,a_Jmax::T,b_Jmax::T,i::T,Kₓₗ₀::T) where {T<:Float64}
    data_ind = Find_data_ind(j)
    H = data.H[data_ind]
    N = data.N[data_ind]
    Wf = data.Wf[data_ind]
    B = data.B[data_ind]
    Hc = data.Hc
    model,kinetic = Initi_model_struct(H,Hc,N,Wf,B,αf,β₁,β₂,Nₛ,rₘ,a_Jmax,b_Jmax,i,Kₓₗ₀)

    return model,kinetic
end

function Run_RO_C_F_CCPH(par::Array{Float64,1},weatherts_F::WeatherTS,
    weatherts_C::WeatherTS,data_F::RoData,data_C::RoData;sim_steps::Integer=83)

    αf,β₁,β₂,b_Jmax = 460.0,1.27,-0.27,0.0
    Nₛ,rₘ,a_Jmax,Kₓₗ₀,Nₛ_C,i = par[1:6]

    #Fertilized stand
    GPP_model_F = zeros(sim_steps+1)
    gₛ_model_F = zeros(sim_steps+1)
    for j = 1:sim_steps+1

        model,kinetic = Initi_model_struct(j,data_F,αf,β₁,β₂,Nₛ,rₘ,a_Jmax,b_Jmax,i,Kₓₗ₀)

        GPP_model_F[j],gₛ_model_F[j] = CalcGPP!(j,model,weatherts_F,kinetic)
    end

    #Control stand
    GPP_model_C = zeros(sim_steps+1)
    gₛ_model_C = zeros(sim_steps+1)
    for j = 1:sim_steps+1

        model,kinetic = Initi_model_struct(j,data_C,αf,β₁,β₂,Nₛ_C,rₘ,a_Jmax,b_Jmax,i,Kₓₗ₀)

        GPP_model_C[j],gₛ_model_C[j] = CalcGPP!(j,model,weatherts_C,kinetic)
    end

    return GPP_model_F,GPP_model_C,gₛ_model_F,gₛ_model_C
end

function RO_C_F_CCPH_ModelOutput(par::Array{Float64,1},weatherts_F::WeatherTS,
    weatherts_C::WeatherTS,data_F::RoData,data_C::RoData;sim_steps::Integer=83)

    αf,β₁,β₂,b_Jmax = 460.0,1.27,-0.27,0.0
    Nₛ,rₘ,a_Jmax,Kₓₗ₀,Nₛ_C,i = par[1:6]

    #Fertilized stand
    ModelOutput_F = CCPHOutput[]
    k_cost_max_F = Float64[]
    gₛ_F = Float64[]
    gₛ_crit_F = Float64[]
    Nₘ_f_F = Float64[]
    for j = 1:sim_steps+1

        model,kinetic = Initi_model_struct(j,data_F,αf,β₁,β₂,Nₛ,rₘ,a_Jmax,b_Jmax,i,Kₓₗ₀)
        modeloutput,gₛ_opt,Nₘ_f_opt = OptimCCPH!(j,model,weatherts_F,kinetic)

        push!(ModelOutput_F,modeloutput)
        push!(k_cost_max_F,Max_k_cost(model))
        push!(gₛ_F,gₛ_opt)
        push!(gₛ_crit_F,CCPH.Calc_K_costᵢₙᵥ(0.12,model))
        push!(Nₘ_f_F,Nₘ_f_opt)
    end

    #Control stand
    ModelOutput_C = CCPHOutput[]
    k_cost_max_C = Float64[]
    gₛ_C = Float64[]
    gₛ_crit_C = Float64[]
    Nₘ_f_C = Float64[]
    for j = 1:sim_steps+1

        model,kinetic = Initi_model_struct(j,data_C,αf,β₁,β₂,Nₛ_C,rₘ,a_Jmax,b_Jmax,i,Kₓₗ₀)
        modeloutput,gₛ_opt,Nₘ_f_opt = OptimCCPH!(j,model,weatherts_C,kinetic)

        push!(ModelOutput_C,modeloutput)
        push!(k_cost_max_C,Max_k_cost(model))
        push!(gₛ_C,gₛ_opt)
        push!(gₛ_crit_C,CCPH.Calc_K_costᵢₙᵥ(0.12,model))
        push!(Nₘ_f_C,Nₘ_f_opt)
    end

    return ModelOutput_F,ModelOutput_C,k_cost_max_F,k_cost_max_C,gₛ_F,gₛ_C,gₛ_crit_F,gₛ_crit_C,Nₘ_f_F,Nₘ_f_C
end

function Calc_error_RO_C_F_CCPH(
    GPP_model_F::Array{Float64,1},gₛ_model_F::Array{Float64,1},data_F::RoData,
    GPP_model_C::Array{Float64,1},gₛ_model_C::Array{Float64,1},data_C::RoData)

    #[21,44,64,84]
    ind_start = [1,22,45,65]
    ind_end = [21,44,64,84]
    #Fertilized stand

    error_F = 0.0
    for i = 1:4
        error_F += 1-calcR²(data_F.GPP[ind_start[i]:ind_end[i]],GPP_model_F[ind_start[i]:ind_end[i]])
        error_F += 1-calcR²(data_F.gₛ[ind_start[i]:ind_end[i]],gₛ_model_F[ind_start[i]:ind_end[i]])
    end

    #Control stand
    error_C = 0.0
    for i = 1:4
        error_C += 1-calcR²(data_C.GPP[ind_start[i]:ind_end[i]],GPP_model_C[ind_start[i]:ind_end[i]])
        error_C += 1-calcR²(data_C.gₛ[ind_start[i]:ind_end[i]],gₛ_model_C[ind_start[i]:ind_end[i]])
    end

    return  error_C#error_F+error_C
end

function Check_para_RO_C_F_CCPH(par::Array{Float64,1})
    #Nₛ,rₘ,a_Jmax,Kₓₗ₀,Nₛ_C,i,X0,τ,Smax,τ_C
    lo = [0.001,10.0,0.01,0.001,0.001,0.1,-8.0,1.0,10.0,1.0]

    up = [0.1,40.0,1.0,0.1,0.1,6.0,-1.0,15.0,25.0,15.0]

    return any(par.<lo)||any(par.>up)
end

function Post_distri_RO_C_F_CCPH(par::Array{Float64,1},RO_data::RO_raw_data)
    if Check_para_RO_C_F_CCPH(par)
        return post = -Inf
    else
        try
            X0,τ,Smax,τ_C = par[7:10]
            weatherts_F,GPP_data_F = create_weather_struct_RO(RO_data;stand_type="Fertilized",X0=X0,τ=τ,Smax=Smax)
            weatherts_C,GPP_data_C = create_weather_struct_RO(RO_data;stand_type="Control",X0=X0,τ=τ_C,Smax=Smax)

            data_F,data_C = Create_RoData_C_F(GPP_data_F,GPP_data_C,weatherts_F,weatherts_C)

            GPP_model_F,GPP_model_C,gₛ_model_F,gₛ_model_C = Run_RO_C_F_CCPH(par,weatherts_F,weatherts_C,data_F,data_C)

            error_val = Calc_error_RO_C_F_CCPH(
            GPP_model_F,gₛ_model_F,data_F,
            GPP_model_C,gₛ_model_C,data_C)

            return post = -error_val
        catch err
            println("Parameters Error: ", err)
            return post = -Inf
        end
    end
end

function run_simple_trait_model_C_F_tuning_RO_ts()
    #Init Ro data
    RO_data = Load_RO_data()

    nsamples = 5000
    nchains = 12

    #par_guess = [0.036975489745880365, 17.40739827323768, 0.5, 0.01,
    #0.025675473753997178, -4.230533359595931, 12.328105896346644, 19.079489364590533,12.328105896346644]

    par_guess = [0.011976551772617215, 25.42466609569836, 0.5009242298864488, 0.023503724920524185,
    0.010452800544803813,1.0,-1.312709470268466, 3.3784814235281933, 17.510523050796454, 9.566445194647766]

    x_current = [par_guess for i = 1:nchains]
    P_current = [Post_distri_RO_C_F_CCPH(par_guess,RO_data) for i = 1:nchains]

    println(P_current)

    q_vec = [0.04,24.0,0.033,0.01,0.04,1.0,-4.0,7.0,16.0,7.0]/25.0
    metropolis_mcmc(x_current,P_current,
    x::Array{Float64,1}->Post_distri_RO_C_F_CCPH(x::Array{Float64,1},RO_data),q_vec,"RO_C_F_GPP_Only_Few_Para_NoRoot_20220201";
    n_chains=nchains,n_samples=nsamples,burn_in=7000,sample_freq=5)
end

function run_C_F_ts_mean()

    RO_data = Load_RO_data()

    file_name = "RO_C_F_GPP_Only_Few_Para_NoRoot_20220201"

    samples = load("./output/"*file_name*".jld","samples_container")
    P_samples = load("./output/"*file_name*".jld","P_samples_container")
    Draw_trace_plot(samples,file_name)
    Draw_histogram(samples,file_name)

    par = Find_max_sample(samples,P_samples)

    println(par)

    X0,τ,Smax = par[7:10]
    weatherts_F,GPP_data_F = create_weather_struct_RO(RO_data;stand_type="Fertilized",X0=X0,τ=τ,Smax=Smax)
    weatherts_C,GPP_data_C = create_weather_struct_RO(RO_data;stand_type="Control",X0=X0,τ=τ,Smax=Smax)

    data_F,data_C = Create_RoData_C_F(GPP_data_F,GPP_data_C,weatherts_F,weatherts_C)

    GPP_model_F,GPP_model_C,gₛ_model_F,gₛ_model_C = Run_RO_C_F_CCPH(par,weatherts_F,weatherts_C,data_F,data_C)

    @show GPP_R²_F =  calcR²(data_F.GPP,GPP_model_F)
    @show GPP_cor_F = cor(data_F.GPP,GPP_model_F)
    @show gₛ_R²_F =  calcR²(data_F.gₛ,gₛ_model_F)
    @show gₛ_cor_F = cor(data_F.gₛ,gₛ_model_F)

    @show GPP_R²_C =  calcR²(data_C.GPP,GPP_model_C)
    @show GPP_cor_C = cor(data_C.GPP,GPP_model_C)
    @show gₛ_R²_C =  calcR²(data_C.gₛ,gₛ_model_C)
    @show gₛ_cor_C = cor(data_C.gₛ,gₛ_model_C)

    println("Fertilized")
    clac_GPP_R²_annual(data_C.GPP,GPP_model_F,weatherts_F.date)
    println("Control")
    clac_GPP_R²_annual(data_C.GPP,GPP_model_C,weatherts_C.date)

    plot(weatherts_F.date,data_F.GPP,label="Data")
    pl1 = plot!(weatherts_F.date,GPP_model_F,label="Model")
    plot(weatherts_C.date,data_C.GPP,label="Data")
    pl2 = plot!(weatherts_C.date,GPP_model_C,label="Model")

    plot(weatherts_F.date,data_F.gₛ,label="Data")
    pl3 = plot!(weatherts_F.date,gₛ_model_F,label="Model")
    plot(weatherts_C.date,data_C.gₛ,label="Data")
    pl4 = plot!(weatherts_C.date,gₛ_model_C,label="Model")

    plot(pl1,pl3,pl2,pl4,layout=(2,2),legends=false)
end

function test_adaptive_rwm()

    file_name = "adaptive_rwm_RO_C_GPP_Only_20220214"
    n_chains = 4
    n_iter = 40000

    RO_data = Load_RO_data()

    par_guess = [0.011976551772617215, 25.42466609569836, 0.5009242298864488, 0.023503724920524185,
    0.010452800544803813,1.0,-1.312709470268466, 3.3784814235281933, 17.510523050796454, 9.566445194647766]

    out_mat_tot = Array{Matrix{Float64},1}(undef,n_chains)

    Threads.@threads for j = 1:n_chains
        # Run n_iter iterations of the Adaptive Metropolis method:
        out = AdaptiveMCMC.adaptive_rwm(par_guess,
        x::Array{Float64,1}->Post_distri_RO_C_F_CCPH(x::Array{Float64,1},RO_data),
        n_iter; algorithm=:am)

        out_mat_tot[j] = transpose(out.X)
    end

    try
        save("./output/"*file_name*".jld","samples",out_mat_tot)

        @show typeof(out_mat_tot)
        @show size(out_mat_tot)

        n_features = size(out_mat_tot[1])[2]

        histogram(out_mat_tot[1],layout=n_features,legend=false)
        savefig("./plots/"*file_name*".svg")

        plot(;layout=n_features,legends=false)
        for i = 1:n_features
            for j=1:n_chains
                plot!(out_mat_tot[j][:,i],subplot=i)
            end
        end
        savefig("./plots/"*file_name*"_trace_plot.svg")

        # Calculate '95% credible intervals':
        para_stat = mapslices(x->"$(mean(x)) ± $(1.96std(x))", out_mat_tot[1], dims=1)
        for line in para_stat
            println(line)
        end
    catch err
        println("Error: ",err)
    end
end

function find_max_adaptive_rwm()

    n_chains = 12
    file_name = "adaptive_rwm_RO_C_F_GPP_Only_20220203"
    out_mat_tot = load("./output/"*file_name*".jld","samples")
    out_mat = out_mat_tot[1]
    n_features = size(out_mat)[2]
    n_samples = size(out_mat)[1]

    RO_data = Load_RO_data()

    #=
    log_p_vec = Array{Float64,1}(undef,n_samples)

    Threads.@threads for j = 1:n_samples
        println("Done: $(j)/$(n_samples)")
        para_val = out_mat[j,:]
        log_p_vec[j] = Post_distri_RO_C_F_CCPH(para_val,RO_data)
    end

    save("./output/"*file_name*"_log_p.jld","log_p",log_p_vec)
    =#

    log_p_vec = load("./output/"*file_name*"_log_p.jld","log_p")

    @show min_ind = argmax(log_p_vec)

    @show size(out_mat)
    @show typeof(out_mat)

    @show par = reshape(mean(out_mat, dims=1),n_features)
    @show par = out_mat[min_ind,:]

    X0,τ,Smax,τ_C = par[7:10]
    weatherts_F,GPP_data_F = create_weather_struct_RO(RO_data;stand_type="Fertilized",X0=X0,τ=τ,Smax=Smax)
    weatherts_C,GPP_data_C = create_weather_struct_RO(RO_data;stand_type="Control",X0=X0,τ=τ_C,Smax=Smax)

    data_F,data_C = Create_RoData_C_F(GPP_data_F,GPP_data_C,weatherts_F,weatherts_C)

    GPP_model_F,GPP_model_C,gₛ_model_F,gₛ_model_C = Run_RO_C_F_CCPH(par,weatherts_F,weatherts_C,data_F,data_C)

    @show GPP_R²_F =  calcR²(data_F.GPP,GPP_model_F)
    @show GPP_cor_F = cor(data_F.GPP,GPP_model_F)
    @show gₛ_R²_F =  calcR²(data_F.gₛ,gₛ_model_F)
    @show gₛ_cor_F = cor(data_F.gₛ,gₛ_model_F)

    @show GPP_R²_C =  calcR²(data_C.GPP,GPP_model_C)
    @show GPP_cor_C = cor(data_C.GPP,GPP_model_C)
    @show gₛ_R²_C =  calcR²(data_C.gₛ,gₛ_model_C)
    @show gₛ_cor_C = cor(data_C.gₛ,gₛ_model_C)

    println("Fertilized")
    clac_GPP_R²_annual(data_C.GPP,GPP_model_F,weatherts_F.date)
    println("Control")
    clac_GPP_R²_annual(data_C.GPP,GPP_model_C,weatherts_C.date)

    plot(weatherts_F.date,data_F.GPP,label="Data",ylabel="GPP")
    pl1 = plot!(weatherts_F.date,GPP_model_F,label="Model")
    plot(weatherts_C.date,data_C.GPP,label="Data",ylabel="GPP")
    pl2 = plot!(weatherts_C.date,GPP_model_C,label="Model")

    plot(weatherts_F.date,data_F.gₛ,label="Data",ylabel="gₛ")
    pl3 = plot!(weatherts_F.date,gₛ_model_F,label="Model")
    plot(weatherts_C.date,data_C.gₛ,label="Data",ylabel="gₛ")
    pl4 = plot!(weatherts_C.date,gₛ_model_C,label="Model")

    plot(pl1,pl3,pl2,pl4,layout=(2,2),legends=false)
    savefig("./plots/"*file_name*"_result.svg")

    # Calculate '95% credible intervals':
    para_stat = mapslices(x->"$(mean(x)) ± $(1.96std(x))", out_mat, dims=1)
    for line in para_stat
        println(line)
    end
end

#run_simple_trait_model_C_F_tuning_RO_ts()
#run_C_F_ts_mean()
#test_adaptive_rwm()
find_max_adaptive_rwm()


text:
"""
Particles struct.
fields: positon x, velocity v, weight w
"""
struct ParticleList{T}
    x::Vector{T}
    v::Vector{T}
    w::Vector{T}
end

Base.length(pl::ParticleList) = length(pl.w)