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

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

export BrickletMoisture
"""
Measures soil moisture
"""
mutable struct BrickletMoisture <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

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

        return new(ipcon, deviceInternal)
    end
end

export get_moisture_value
"""
    $(SIGNATURES)

Returns the current moisture value.
A small value corresponds to little moisture, a big
value corresponds to much moisture.

If you want to get the moisture value periodically, it is recommended
to use the :cb:`Moisture` callback and set the period with
:func:`Set Moisture Callback Period`.
"""
function get_moisture_value(device::BrickletMoisture)
    return device.deviceInternal.get_moisture_value()
end

export set_moisture_callback_period
"""
    $(SIGNATURES)

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

The :cb:`Moisture` callback is only triggered if the moisture value has changed
since the last triggering.
"""
function set_moisture_callback_period(device::BrickletMoisture, period)
    device.deviceInternal.set_moisture_callback_period(period)
end

export get_moisture_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Moisture Callback Period`.
"""
function get_moisture_callback_period(device::BrickletMoisture)
    return device.deviceInternal.get_moisture_callback_period()
end

export set_moisture_callback_threshold
"""
    $(SIGNATURES)

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

The following options are possible:

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

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

export get_moisture_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Moisture Callback Threshold`.
"""
function get_moisture_callback_threshold(device::BrickletMoisture)
    return device.deviceInternal.get_moisture_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callback

* :cb:`Moisture Reached`

is triggered, if the threshold

* :func:`Set Moisture Callback Threshold`

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

export get_debounce_period
"""
    $(SIGNATURES)

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

export set_moving_average
"""
    $(SIGNATURES)

Sets the length of a `moving averaging <https://en.wikipedia.org/wiki/Moving_average>`__
for the moisture value.

Setting the length to 0 will turn the averaging completely off. With less
averaging, there is more noise on the data.
"""
function set_moving_average(device::BrickletMoisture, average)
    device.deviceInternal.set_moving_average(average)
end

export get_moving_average
"""
    $(SIGNATURES)

Returns the length moving average as set by :func:`Set Moving Average`.
"""
function get_moving_average(device::BrickletMoisture)
    return device.deviceInternal.get_moving_average()
end

export get_identity
"""
    $(SIGNATURES)

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

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

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

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


text:
using Base.Test

tests = [
    "domain",
    "group",
    "waveletclass",
    "path",
    "meta",
    "spec",
    "spec1d",
    "spec2d",
    "filter",
    "fourierfilter",
    "weighting",
    "behavior",
    "morlet1d",
    "bank",
    "bank1d",
    "blob",
    "pointwise",
    "modulus",
    "log1p",
    "symbols",
    "layer",
    "inputlayer",
    "pointwiselayer",
    "fourierlayer",
    "waveletlayer",
    "inputlayerstate",
    "pointwiselayerstate",
    "fourierlayerstate",
    "net"
]

println("Running tests:")

for t in tests
    tfile = string(t, ".jl")
    println(" * ", tfile, " ...")
    include(tfile)
end

ext: jl


text:
export hausdorff_distance

"""
    hausdorff_distance(X::LazySet{N}, Y::LazySet{N}; [p]::N=N(Inf),
                       [ε]=N(1e-3)) where {N}

Compute the Hausdorff distance between two convex sets up to a given threshold.

### Input

- `X` -- convex set
- `Y` -- convex set
- `p` -- (optional, default: `Inf`) norm parameter of the Hausdorff distance
- `ε` -- (optional, default: `1e-3`) precision threshold; the true Hausdorff
         distance may diverge from the result by at most this value

### Output

A value from the ``ε``-neighborhood of the Hausdorff distance between ``X`` and
``Y``.

### Notes

Given a ``p``-norm, the Hausdorff distance ``d_H^p(X, Y)`` between sets ``X``
and ``Y`` is defined as follows:

```math
    d_H^p(X, Y) = \\inf\\{δ ≥ 0 \\mid Y ⊆ X ⊕ δ 𝐵_p^n \\text{ and } X ⊆ Y ⊕ δ 𝐵_p^n\\}
```

Here ``𝐵_p^n`` is the ``n``-dimensional unit ball in the ``p``-norm.

The implementation may internally rely on the support function of ``X`` and
``Y``; hence any imprecision in the implementation of the support function may
affect the result.
At the time of writing, the only set type with imprecise support function is the
lazy [`Intersection`](@ref).

### Algorithm

We perform binary search for bounding the Hausdorff distance in an interval
``[l, u]``, where initially ``l`` is ``0`` and ``u`` is described below.
The binary search terminates when ``u - l ≤ ε``, i.e., the interval becomes
sufficiently small.

To find an upper bound ``u``, we start with the heuristics of taking the biggest
distance in the axis-parallel directions.
As long as this bound does not work, we increase the bound by ``2``.

Given a value ``δ``, to check whether the sets are within Hausdorff distance
``δ``, we simply check the inclusions given above, where on the right-hand side
we use a lazy `MinkowskiSum` with a `Ballp` centered in the origin.
"""
function hausdorff_distance(X::LazySet{N}, Y::LazySet{N}; p::N=N(Inf),
                            ε=N(1e-3)) where {N}
    @assert ε > zero(N) "the value ε must be positive"
    @assert isbounded(X) && isbounded(Y) "the Hausdorff distance is only " *
        "defined for compact sets"

    n = dim(X)
    @assert dim(Y) == n "the Hausdorff distance is only defined between sets " *
                        "of the same dimension, but they had dimensions $n " *
                        "resp. $(dim(Y))"

    # phase 1: find a finite upper bound
    δ_upper = max(maximum(d -> abs(ρ(d, X) - ρ(d, Y)), BoxDirections{N}(n)),
                  N(1e-3))  # this initial bound should be strictly positive
    # verify that this is an upper bound
    while !_mutual_issubset_in_δ_bloating(X, Y, δ_upper, n, p)
        δ_upper *= N(2)
    end

    # phase 2: perform binary search between lower bound (initially 0) and upper
    # bound until convergence
    δ_lower = N(0)
    while δ_upper - δ_lower > ε
        δ = (δ_upper + δ_lower) / N(2)
        if _mutual_issubset_in_δ_bloating(X, Y, δ, n, p)
            δ_upper = δ
        else
            δ_lower = δ
        end
    end
    return δ_upper
end

function _mutual_issubset_in_δ_bloating(X, Y, δ, n, p)
    return _issubset_in_δ_bloating(X, Y, δ, n, p) &&
           _issubset_in_δ_bloating(Y, X, δ, n, p)
end

function _issubset_in_δ_bloating(X::LazySet{N}, Y, δ, n, p) where {N}
    return X ⊆ Y + Ballp(p, zeros(N, n), δ)
end

# for polytopes the default implementation of `⊆` requires membership in the rhs
# set, which will be a MinkowskiSum and hence not available; we use the
# alternative based on constraints_list on the right instead
function _issubset_in_δ_bloating(X::AbstractPolytope{N}, Y, δ, n, p) where {N}
    return LazySets._issubset_constraints_list(X, Y + Ballp(p, zeros(N, n), δ))
end


text:
using StanBase

stan_prog = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
  real empty[0];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

stanmodel = HelpModel( "help", stan_prog;
  method = StanBase.Help(help=:sample))

stan_help(stanmodel;n_chains=1)

run(`cat $(stanmodel.output_base*"_chain_1.log")`)


text:
export Map_Translate

# ****************************************************************

@doc raw"""
```julia
Map_Translate(;scale = 1)
```
Defines a translate transformation:

```math
X(\hat{X}) = \hat{X} + scale (θ - 1)
```

In peculiar, for ``(θ) = (1)`` the transformation is identity.
This is due to the presence of the ``θ- 1`` factor. This factor is introduced to stay consistent with other transformations where we
try to have identity transform for parameter vector of ones.
"""
struct Map_Translate{T} <: Abstract_Map
    _scale::T
end

function Map_Translate(;scale=1)
    Map_Translate{typeof(scale)}(scale)
end

parameter_size(m::Map_Translate) = 1

function eval_map(m::Map_Translate{T},
                  X_hat::AbstractVector{T_HAT},
                  θ::AbstractVector{T_θ}) where {T_HAT,T,T_θ}

    @assert length(θ) == parameter_size(m)

    θ = θ[1]

    X_T = promote_type(T_HAT,T,T_θ)
    X=similar(X_hat,X_T)

    @. X = X_hat + m._scale*(θ-1)

    X
end


text:
using CSV, JuliaDB, DungBase, Dates, UUIDs, Tables, TableOperations
import IntervalSets: width, (..), AbstractInterval, leftendpoint
import DataFrames: DataFrame
import Base.Threads: @spawn
# format2millisecond(_::Missing) = missing
# function format2millisecond(x)
#     m = match(r"^(\d*)\smilliseconds?$", x)
#     isnothing(m) && return missing
#     Millisecond(parse(Int, first(m.captures)))
# end
tons(_::Missing) = missing
tons(x) = Nanosecond(x)
function DungBase.AbstractTimeLine(x)
    files = Vector{VideoFile}(undef, length(x))
    for i in rows(x)
        files[i.index] = VideoFile(i.file_name, i.date_time, i.duration)
    end
    y = x[1]
    comment = ismissing(y.comment) ? "" : y.comment
    AbstractTimeLine(files, comment)
end
DungBase.AbstractPeriod(x, y::Missing) = AbstractPeriod(x)
function DungBase.Temporal(xs)
    x = xs[]
    video = x.AbstractTimeLine
    p1 = Millisecond(Nanosecond(x.start))
    p = ismissing(x.stop) ? p1 : p1..Millisecond(Nanosecond(x.stop))
    time = AbstractPeriod(p)
    try
        Temporal(video, time, ismissing(x.comment) ? "" : x.comment)
    catch
        Temporal(video, Prolonged(), ismissing(x.comment) ? "the time is made up because it was missing" : x.comment*" the time is made up because it was missing")
    end
end
function DungBase.Calibration(_x)
    x = _x[]
    board = Board(x.board, x.checker_width_cm, (x.checker_per_width, x.checker_per_height), x.board_description)
    comment = ismissing(x.comment) ? "" : x.comment
    Calibration(x.intrinsic, x.extrinsic, board, comment)
end
DungBase.POI(x) = POI(x.Calibration, x.Temporal)
function DungBase.Run(x)
    @assert allunique(x.type) "POI types are not unique: $(x.type)"
    pois = Dict(Symbol(i.type) => i.POI for i in x)
    tmp = select(table(x), Not(All(:date, :comment, :POI, :type)))
    # CSV.write("/home/yakir/coffeebeetlearticle/all/tmp.csv", tmp)
    setup = Dict(k => v for (k, v) in pairs(tmp[1]) if !ismissing(v))
    _x = x[1]
    date = _x.date
    comment = ismissing(_x.comment) ? "" : _x.comment
    metadata = Metadata(setup, comment, date)
    Run(pois, metadata)
end

function loaddeomcsv(source)

    t = CSV.File(joinpath(source, "video.csv")) |> TableOperations.transform(video = UUID)
    video = table(t, pkey = :video)
    t = CSV.File(joinpath(source, "videofile.csv")) |> TableOperations.transform(video = UUID, duration = tons)
    # t = CSV.File(joinpath(source, "videofile.csv")) |> TableOperations.transform(video = UUID, duration = format2millisecond)
    videofile = table(t, pkey = :file_name)
    x = join(video, videofile, rkey = :video)
    timeline = groupby(AbstractTimeLine,  x)
    t = CSV.File(joinpath(source, "interval.csv")) |> TableOperations.transform(interval = UUID, video = x -> ismissing(x) ? missing : UUID(x), start = tons, stop = tons)
    # t = CSV.File(joinpath(source, "interval.csv")) |> TableOperations.transform(interval = UUID, video = UUID, start = format2millisecond, stop = format2millisecond)
    interval = table(t, pkey = :interval)
    # interval = dropmissing(loadtable(joinpath(source, "interval.csv"), indexcols = :interval), :start)
    x = join(interval, timeline, lkey = :video)
    temporal = dropmissing(groupby(Temporal, x, :interval), :Temporal)

    board = loadtable(joinpath(source, "board.csv"), indexcols = :designation)
    t = CSV.File(joinpath(source, "calibration.csv")) |> TableOperations.transform(calibration = UUID, intrinsic = x -> ismissing(x) ? missing : UUID(x), extrinsic = UUID, board = String)
    calibration = table(t, pkey = :calibration)
    x = join(calibration, board, lkey = :board, rkey = :designation)
    x = join(x, temporal, rkey = :interval, lkey = :extrinsic)
    x = select(x, Not(:extrinsic))
    x = JuliaDB.rename(x, :Temporal => :extrinsic)
    x = join(x, temporal, rkey = :interval, lkey = :intrinsic, how = :left)
    x = select(x, Not(:intrinsic))
    x = JuliaDB.rename(x, :Temporal => :intrinsic)
    calibration = groupby(Calibration, x, :calibration)

    t = CSV.File(joinpath(source, "poi.csv")) |> TableOperations.transform(poi = UUID, run = UUID, calibration = x -> ismissing(x) ? missing : UUID(x), interval = x -> ismissing(x) ? missing : UUID(x))
    poi = table(t, pkey = :poi)
    coord = UUID.(first.(splitext.(readdir(joinpath(source, "pixel")))))
    poi = filter(in(coord), poi, select = :interval)
    x = join(poi, temporal, lkey = :interval)
    x = join(x, calibration, lkey = :calibration)
    x = select(x, Not(:calibration))
    x = JuliaDB.transform(x, :POI => POI.(x))

    temporal2pixel = Dict(Pair(Symbol(i.type), i.POI) => joinpath(joinpath(source, "pixel"), string(i.interval, ".csv")) for i in x)

    # JLSO.save(joinpath(source, "temporal2pixel.jlso"), "food" => "☕🥓🍳", "cost" => 11.95, "time" => Time(9, 0))

    poi = select(x, Not(All(:Temporal, :Calibration, :poi, :interval)))
    t = CSV.File(joinpath(source, "run.csv")) |> TableOperations.transform(run = UUID)
    runs = table(t, pkey = :run)
    # x = JuliaDB.transform(runs, :date => :date => Date)
    x = join(runs, poi, rkey = :run)
    experiment = loadtable(joinpath(source, "experiment.csv"), indexcols = :experiment)
    x = join(x, experiment, lkey = :experiment)

    data = groupby(x, :experiment, usekey = true) do k, xx
        runs = groupby(Run, select(table(xx), Not(All(:experiment_folder, :experiment_description))), :run)
        (Experiment = Experiment(select(runs, :Run), xx[1].experiment_description), name = k.experiment)
        # (Experiment = Experiment(select(runs, :Run), xx[1].experiment_description), name = xx[1].experiment_folder)
    end

    data = Dict(i.name => i.Experiment for i in data)

    return temporal2pixel, data

    # JLSO.save(joinpath(datadep"coffeebeetle", "data.jlso"), data)
end

goodvideo(file) = first(file) ≠ '.' && occursin(r"mts|mp4|avi|mpg|mov|mkv"i, last(splitext(file)))

function joinsources(sources; prefix = "source_")
    source = mktempdir(; prefix = prefix, cleanup = false)
    pixel = joinpath(source, "pixel")
    mkdir(pixel)
    @sync for s in sources
        @spawn begin
            pp = joinpath(s, "pixel")
            for p in readdir(pp)
                @assert p ∉ readdir(pixel) "duplicate pixel files"
                # mv(joinpath(pp, p), joinpath(pixel, p))
                cp(joinpath(pp, p), joinpath(pixel, p))
            end
            # rm(pp)
        end
    end
    @sync for x in ("board", "calibration", "experiment", "interval", "poi", "video", "videofile")
        @spawn begin

            open(joinpath(source, "$x.csv"), "w") do i
                open(joinpath(sources[1], "$x.csv")) do o
                    write(i, read(o))
                end
                for s in Iterators.drop(sources, 1)
                    open(joinpath(s, "$x.csv")) do o
                        readuntil(o, '\n')
                        write(i, read(o))
                    end
                end
            end

            # for s in sources
                # rm(joinpath(s, "$x.csv"))
            # end
        end
    end

    runs = vcat((DataFrame(CSV.File(joinpath(s, "run.csv"))) for s in sources)..., cols = :union)
    CSV.write(joinpath(source, "run.csv"), runs)

    @sync for s in sources
        @spawn begin
            for p in readdir(s)
                if goodvideo(p)
                    @assert p ∉ readdir(source) "duplicate video files"
                    cp(joinpath(s, p), joinpath(source, p))
                end
            end
        end
    end
    # for s in sources
    #     @assert isempty(readdir(s)) "some files remain?"
    #     rm(s)
    # end
    return source
end


text:
module RootsNewton

greet() = print("Hello World!")

using LinearAlgebra, ForwardDiff

#This function finds the roots of f(x) taking iv as initial point. The derivative is calculated using ForwardDiff
function newtonroot(f; iv, tolerance=1E-7, maxiter=1000)
    D(f) = x -> ForwardDiff.derivative(f, x)
    f_prime=D(f)
    return newtonroot(f,f_prime;iv=iv,tolerance=tolerance,maxiter=maxiter)
end

#This function finds the roots of f(x) taking iv as initial point
function newtonroot(f,f_prime; iv, tolerance=1E-7, maxiter=1000)
    # setup the algorithm
    x_old = iv
    normdiff = Inf
    iter = 1
    while normdiff > tolerance && iter <= maxiter
        x_new = x_old - f(x_old)/f_prime(x_old) # use the passed in map
        normdiff = norm(x_new - x_old)
        x_old = x_new
        iter = iter + 1
    end
	if iter>=maxiter || isnan(x_old)
		return (value = nothing, funcRes=nothing, iter=iter-1) # A named tuple
	else
		return (value = x_old, funcRes=f(x_old), iter=iter-1) # A named tuple
	end
end

export newtonroot

end # module


text:
module MAEBOL

import StatsBase: entropy, WeightVec, sample
using Distributions: VonMises, Normal, cdf
using PyPlot: imshow, plot, xlabel, ylabel, contour, figure, pause, hold
using Cubature
export figure

export SearchDomain, initial_belief, update_belief!, transition!
export Vehicle, VehicleSet
export O, observe, true_bearing
export fisher, EID
export plot_b, plot_eid, plot_mi
export plot_eid2, plot_eid3, meshgrid, plot_sim
export h_ij, hk2, update_phik!, phi_ij, fk
export ErgodicManager
export get_action, RandPolicy, ConstantPolicy, SMC, SMC2, FisherGreedy
export FisherGreedy2
export Simulation
export mutual_information

typealias Obs Int64
typealias Observation Obs
typealias ObsSet Vector{Obs}
typealias Action NTuple{2, Float64}
typealias ActionSet Vector{Action}
typealias Belief Matrix{Float64}

include("vehicleset.jl")
include("searchdomain.jl")
include("observations.jl")
include("ergodicity.jl")
include("eid.jl")
include("policy.jl")
include("simulation.jl")
include("infotheoretic.jl")
include("plotting.jl")

# TODO: add checks for staying within
"""
`transition(m, x, a)`

Returns vehicle state after taking action `a`. Not currently working.
"""
function transition!(m::SearchDomain, X::VehicleSet, A::ActionSet)
	for i = 1:length(X)
		X[i].x = min(max(X[i].x + A[i][1], 0), m.num_cells)
		X[i].y = min(max(X[i].y + A[i][2], 0), m.num_cells)
	end
end

# Finds the next closest cell you end up at
function nextcell(m::SearchDomain, v::Vehicle, a::Action)
	# Determine grid cell you belong to
	grid_x = round(Int, v.x + a[1], RoundDown)
	grid_y = round(Int, v.y + a[2], RoundDown)

	# Keep grid cell within bounds
	grid_x = min(max(grid_x, 0), m.num_cells-1)
	grid_y = min(max(grid_y, 0), m.num_cells-1)
	return grid_x, grid_y
end

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using DataFrames, CSV, Statistics
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing4.12"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	data = CSV.read(datadir("temperatures.csv"), copycols=true)
	brisT = data.Brisbane
	gcT = data.GoldCoast

	sigB = std(brisT)
	sigG = std(gcT)
	covBG = cov(brisT, gcT)

	meanVect = [mean(brisT) , mean(gcT)]
	covMat = [sigB^2  covBG
			  covBG   sigG^2]
end

# ╔═╡ 65fa5b06-246e-11eb-19b7-d77750e343fb
begin
	outfile = open(datadir("mvParams2.jl"), "w")
	write(outfile,"meanVect = $meanVect \ncovMat = $covMat")
	close(outfile)
	read(datadir("mvParams2.jl"), String)
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing4.12"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═65fa5b06-246e-11eb-19b7-d77750e343fb
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
mutable struct NetDecModel{T <: PM.AbstractPowerModel}
    node_partition::Vector{Set{Int64}}
    model_list::Vector{T}
    split_line_vars::Dict{Int64, Dict}
    function NetDecModel(
        node_partition::Vector{Set{Int64}},
        model_list::Vector{T} = T[],
        split_line_vars::Dict{Int64, Dict} = Dict{Int64, Dict}()
        ) where T <: PM.AbstractPowerModel
        return new{T}(node_partition, model_list, split_line_vars)
    end
end

get_partition(dm_model::NetDecModel)::Vector{Set{Int64}} = dm_model.node_partition

get_models(dm_model::NetDecModel)::Vector{<:PM.AbstractPowerModel} = dm_model.model_list

get_split_vars(dm_model::NetDecModel)::Dict{Int64, Dict} = dm_model.split_line_vars

function set_subnet_optimizer!(dm_model::NetDecModel, optimizer)
    for pm in get_models(dm_model)
        JuMP.set_optimizer(pm.model, optimizer)
    end
    return
end


text:
using OrdinaryDiffEq, StochasticDiffEq, Flux, DiffEqFlux

x = Float32[2.; 0.]
tspan = (0.0f0,25.0f0)
dudt = Chain(Dense(2,50,tanh),Dense(50,2))

neural_ode(dudt,x,tspan,Tsit5(),save_everystep=false,save_start=false)
neural_ode(dudt,x,tspan,Tsit5(),saveat=0.1)
neural_ode_rd(dudt,x,tspan,Tsit5(),saveat=0.1)

Flux.back!(sum(neural_ode(dudt,x,tspan,Tsit5(),saveat=0.0:0.1:10.0)))
Flux.back!(sum(neural_ode_rd(dudt,x,tspan,Tsit5(),saveat=0.1)))

mp = Float32[0.1,0.1]
neural_dmsde(dudt,x,mp,tspan,SOSRI(),saveat=0.1)
Flux.back!(sum(neural_dmsde(dudt,x,mp,tspan,SOSRI(),saveat=0.1)))


text:
function loadFortranData(dir)

    Nν, Nω = parse.(Int64, split(readlines(dir * "/ladderDGA.in")[4]))[1:2]
    Nq_x, NInt, Nk = parse.(Int64, split(readlines(dir * "/ladderDGA.in")[6]))
    Nq = floor(Int64,(Nq_x+2)*(Nq_x+1)*Nq_x/6)
    bubble =  Array{ComplexF64}(undef, 2*Nω+1, 2*Nν, Nq)
    χch =  Array{ComplexF64}(undef, 2*Nω+1, Nq)
    χsp =  Array{ComplexF64}(undef, 2*Nω+1, Nq)
    trilexch = Array{ComplexF64}(undef, 2*Nω+1, 2*Nν, Nq)
    trilexsp = Array{ComplexF64}(undef, 2*Nω+1, 2*Nν, Nq)
    Σ = Array{ComplexF64}(undef, 2*Nν, Nk)
    for ωi in 0:2*Nω
        istr = lpad(ωi,3,"0")
        f = FortranFile(String(dir * "/chi_bubble_raw/chi" * istr), "r")
        bubble[ωi+1,:,:] = read(f, (ComplexF64, size(bubble)[2:end]...))
        close(f)
        f = FortranFile(dir * "/chich_omega_raw/chi" * istr, "r")
        χch[ωi+1,:] = read(f, (ComplexF64, size(χch)[2]))
        close(f)
        f = FortranFile(dir * "/chisp_omega_raw/chi" * istr, "r")
        χsp[ωi+1,:] = read(f, (ComplexF64, size(χsp)[2]))
        close(f)
        f = FortranFile(dir * "/trilexch_omega_raw/tri" * istr, "r")
        trilexch[ωi+1,:,:] = read(f, (ComplexF64, size(trilexch)[2:end]...))
        close(f)
        f = FortranFile(dir * "/trilexsp_omega_raw/tri" * istr, "r")
        trilexsp[ωi+1,:,:] = read(f, (ComplexF64, size(trilexsp)[2:end]...))
        close(f)
    end
    bubble = permutedims(bubble, [1,3,2])
    trilexsp = permutedims(trilexsp, [1,3,2])
    trilexch = permutedims(trilexch, [1,3,2])
    f = FortranFile(dir * "/SelfE.dat", "r")
    Σ = read(f, (ComplexF32, size(Σ)...))
    close(f)
    return [bubble, χch, χsp, trilexch, trilexsp, Σ]
end


text:
using Revise
using CUDA
using Test
using BenchmarkTools
using LinearAlgebra
using NNlib

include("CUDA_getri.jl")

# Needs some warning if memory is getting too large !
FT = Float32
n = 32
nSpec = 10000

# Create CPU Matrices:
A_ = (randn(FT, n, n, nSpec));
B_ = (randn(FT, n, n, nSpec));
X_ = (randn(FT, n, n, nSpec));

orig_A = deepcopy(A_)
orig_B = deepcopy(B_)
orig_X = deepcopy(X_)

# And move to GPU as CuArray
A = CuArray(A_);
B = CuArray(B_);
X = CuArray(X_);

# batch Matrix inversion for CPU
function batch_solve!(X::CuArray{FT,3}, A::CuArray{FT,3}, B::CuArray{FT,3}) where {FT}
    temp = similar(A)
    pivot, info   = CUBLAS.getrf_strided_batched!(A, true);
    synchronize()
    getri_strided_batched!(A, temp, pivot); # inv stored in aux1
    NNlib.batched_mul!(X, temp, B)
    synchronize()
end

# batch Matrix inversion for CPU
function batch_solve!(X::AbstractArray{FT,3}, A::AbstractArray{FT,3}, B::AbstractArray{FT,3}) where {FT}
    for i = 1:size(A, 3)
        # @views X[:,:,i] = inv(A[:,:,i])
        @views ldiv!(X[:,:,i], qr!(A[:,:,i]), B[:,:,i])
    end
end

@time batch_solve!(X, A, B, );
@time batch_solve!(X_, A_, B_, );

# @test X_ ≈ Array(X) rtol = 1e-2

@time batch_solve!(X, A, B, );
@time batch_solve!(X_, A_, B_, );


text:
using JLD2, Test, RiskAdjustedLinearizations

# Use Wachter model with Disaster Risk to assess the constructors of a RiskAdjustedLinearization type
# for in-place and out-of-place functions
include(joinpath(dirname(@__FILE__), "../../examples/wachter_disaster_risk/wachter.jl"))
m = WachterDisasterRisk()

### In-place RiskAdjustedLinearization

## Deterministic steady state
detout = JLD2.jldopen(joinpath(dirname(@__FILE__), "../reference/det_ss_output.jld2"), "r")
z = vec(detout["z"])
y = vec(detout["y"])
Ψ = zeros(eltype(y), length(y), length(z))
ral = inplace_wachter_disaster_risk(m)

# Check outputs
update!(ral, z, y, Ψ)
nl = nonlinear_system(ral)
li = linearized_system(ral)
@testset "Evaluate WachterDisasterRisk in-place RiskAdjustedLinearization at deterministic steady state" begin
    @test nl.μ_sss ≈ detout["MU"]
    @test nl.Λ.cache ==  detout["LAM"]
    @test nl.Σ.cache.du ≈ detout["SIG"]
    @test nl.ξ_sss ≈ detout["XI"]
    @test nl.𝒱_sss ≈ detout["V"]
    @test li.Γ₁ ≈ detout["GAM1"]
    @test li.Γ₂ ≈ detout["GAM2"]
    @test li.Γ₃ ≈ detout["GAM3"]
    @test li.Γ₄ ≈ detout["GAM4"]
    @test li.Γ₅ ≈ detout["GAM5"]
    @test li.Γ₆ ≈ detout["GAM6"]
    @test li.JV ≈ detout["JV"]
end

## Stochastic steady state
sssout = JLD2.jldopen(joinpath(dirname(@__FILE__), "../reference/iterative_sss_output.jld2"), "r")
z = vec(sssout["z"])
y = vec(sssout["y"])
Ψ = sssout["Psi"]

# Check outputs
update!(ral, z, y, Ψ)
nl = nonlinear_system(ral)
li = linearized_system(ral)
@testset "Evaluate WachterDisasterRisk in-place RiskAdjustedLinearization at stochastic steady state" begin
    @test nl.μ_sss ≈ sssout["MU"]
    @test nl.Λ.cache ==  sssout["LAM"]
    @test nl.Σ.cache.du ≈ sssout["SIG"]
    @test nl.ξ_sss ≈ sssout["XI"]
    @test nl.𝒱_sss ≈ sssout["V"]
    @test li.Γ₁ ≈ sssout["GAM1"]
    @test li.Γ₂ ≈ sssout["GAM2"]
    @test li.Γ₃ ≈ sssout["GAM3"]
    @test li.Γ₄ ≈ sssout["GAM4"]
    @test li.Γ₅ ≈ sssout["GAM5"]
    @test li.Γ₆ ≈ sssout["GAM6"]
    @test li.JV ≈ sssout["JV"]
end

### Out-of-place RiskAdjustedLinearization

## Deterministic steady state
detout = JLD2.jldopen(joinpath(dirname(@__FILE__), "../reference/det_ss_output.jld2"), "r")
z = vec(detout["z"])
y = vec(detout["y"])
Ψ = zeros(eltype(y), length(y), length(z))
ral = outofplace_wachter_disaster_risk(m)

# Check outputs
update!(ral, z, y, Ψ)
nl = nonlinear_system(ral)
li = linearized_system(ral)
@testset "Evaluate WachterDisasterRisk out-of-place RiskAdjustedLinearization at deterministic steady state" begin
    @test nl.μ_sss ≈ detout["MU"]
    @test nl.Λ.cache ==  detout["LAM"]
    @test isnothing(nl.Σ.cache)
    @test nl.Σ(z) ≈ detout["SIG"]
    @test nl.ξ_sss ≈ detout["XI"]
    @test nl.𝒱_sss ≈ detout["V"]
    @test li.Γ₁ ≈ detout["GAM1"]
    @test li.Γ₂ ≈ detout["GAM2"]
    @test li.Γ₃ ≈ detout["GAM3"]
    @test li.Γ₄ ≈ detout["GAM4"]
    @test li.Γ₅ ≈ detout["GAM5"]
    @test li.Γ₆ ≈ detout["GAM6"]
    @test li.JV ≈ detout["JV"]
end

## Stochastic steady state
sssout = JLD2.jldopen(joinpath(dirname(@__FILE__), "../reference/iterative_sss_output.jld2"), "r")
z = vec(sssout["z"])
y = vec(sssout["y"])
Ψ = sssout["Psi"]

# Check outputs
update!(ral, z, y, Ψ)
nl = nonlinear_system(ral)
li = linearized_system(ral)
@testset "Evaluate WachterDisasterRisk out-of-place RiskAdjustedLinearization at stochastic steady state" begin
    @test nl.μ_sss ≈ sssout["MU"]
    @test nl.Λ.cache ==  sssout["LAM"]
    @test isnothing(nl.Σ.cache)
    @test nl.Σ(z) ≈ detout["SIG"]
    @test nl.ξ_sss ≈ sssout["XI"]
    @test nl.𝒱_sss ≈ sssout["V"]
    @test li.Γ₁ ≈ sssout["GAM1"]
    @test li.Γ₂ ≈ sssout["GAM2"]
    @test li.Γ₃ ≈ sssout["GAM3"]
    @test li.Γ₄ ≈ sssout["GAM4"]
    @test li.Γ₅ ≈ sssout["GAM5"]
    @test li.Γ₆ ≈ sssout["GAM6"]
    @test li.JV ≈ sssout["JV"]
end

nothing


text:
using LinearAlgebra
include("FEM1D.jl")
include("containers.jl")
include("Iterations.jl")

function ρ_error(dimension::Integer, k::Integer, iteration::Function, batch_size::Integer)
    x = randn((dimension, batch_size))
    y = x
    for i = 1:k
        y = iteration(y)
    end
    ρ = 0
    n_y, n_x = 0, 0
    for j = 1:batch_size
        ρ += (norm(y[:, j])/norm(x[:, j]))^(1/k)
    end
    ρ = ρ/batch_size
end

function ρ_F(dimension::Integer, k::Integer, iteration::Function, batch_size::Integer)
    x = rand((-1, 1), (dimension, batch_size))
    for l = 1:k
        x = iteration(x)
    end
    ρ = (norm(x)^2/batch_size)^(1/(2*k))
end

function prec_loss_function(parameters::AbstractArray{<:Real, 1}, prec::preconditioner, params::optimization_parameters)
    M = assemble_matrix(parameters, prec)
    iteration(x) = Richardson(x, prec.A, x -> M*x, parameters[end]; prec=prec.symmetry)
    dimension = size(prec.A)[1]
    if params.error == "Frobenius"
        loss = ρ_F(dimension, params.k_sweeps, iteration, params.k_batch)
    elseif params.error == "Spectral"
        loss = ρ_error(dimension, params.k_sweeps, iteration, params.k_batch)
    end
    return loss
end

function spectral_loss_function(parameters::AbstractArray{<:Real, 1}, prec::preconditioner, params::optimization_parameters)
    M = assemble_matrix(parameters, prec)
    iteration(x) = Richardson(x, prec.A, x -> M*x, 1/spectral_radius(parameters, prec, params); prec=prec.symmetry)
    dimension = size(prec.A)[1]
    if params.error == "Frobenius"
        loss = ρ_F(dimension, params.k_sweeps, iteration, params.k_batch)
    elseif params.error == "Spectral"
        loss = ρ_error(dimension, params.k_sweeps, iteration, params.k_batch)
    end
    return loss
end

function scalar_loss_function(parameters::AbstractArray{<:Real, 1}, M::SparseMatrixCSC{Float64,Int64}, prec::preconditioner, params::optimization_parameters)
    iteration(x) = Richardson(x, prec.A, x -> M*x, parameters[end]; prec=prec.symmetry)
    dimension = size(prec.A)[1]
    if params.error == "Frobenius"
        loss = ρ_F(dimension, params.k_sweeps, iteration, params.k_batch)
    elseif params.error == "Spectral"
        loss = ρ_error(dimension, params.k_sweeps, iteration, params.k_batch)
    end
    return loss
end

function ρ_error_T(dimension::Tuple{Integer, Integer}, k::Integer, iteration::Function, batch_size::Integer)
    n_y, n_x = 0, 0
    for j = 1:batch_size
        x = randn(dimension)
        n_x += norm(x)
        y = x
        for i = 1:k
            y = iteration(y)
        end
        n_y += norm(y)
    end
    ρ = (n_y/n_x)^(1/k)
end

function ρ_F_T(dimension::Tuple{Integer, Integer}, k::Integer, iteration::Function, batch_size::Integer)
    n_x = 0
    for j = 1:batch_size
        x = rand((-1, 1), dimension)
        for l = 1:k
            x = iteration(x)
        end
        n_x += norm(x)^2
    end
    ρ = (n_x/batch_size)^(1/(2*k))
end

function prec_loss_function(parameters::AbstractArray{<:Real, 1}, prec::T_preconditioner, params::optimization_parameters)
    N = 2^prec.L+1
    DOF_x, DOF_y = N - 2 + sum(prec.BCs[1:2] .== "N"), N - 2 + sum(prec.BCs[3:4] .== "N")
    dimension = (DOF_y, DOF_x)
    iteration(x) = Richardson_T(x, prec.C, x -> apply_T_preconditioner(x, parameters, prec), parameters[end]; prec=prec.symmetry)
    if params.error == "Frobenius"
        loss = ρ_F_T(dimension, params.k_sweeps, iteration, params.k_batch)
    elseif params.error == "Spectral"
        loss = ρ_error_T(dimension, params.k_sweeps, iteration, params.k_batch)
    end
    return loss
end

function spectral_loss_function(parameters::AbstractArray{<:Real, 1}, prec::T_preconditioner, params::optimization_parameters)
    N = 2^prec.L+1
    DOF_x, DOF_y = N - 2 + sum(prec.BCs[1:2] .== "N"), N - 2 + sum(prec.BCs[3:4] .== "N")
    dimension = (DOF_y, DOF_x)
    iteration(x) = Richardson_T(x, prec.C, x -> apply_T_preconditioner(x, parameters, prec), 1/spectral_radius(parameters, prec, params); prec=prec.symmetry)
    if params.error == "Frobenius"
        loss = ρ_F_T(dimension, params.k_sweeps, iteration, params.k_batch)
    elseif params.error == "Spectral"
        loss = ρ_error_T(dimension, params.k_sweeps, iteration, params.k_batch)
    end
    return loss
end

function Frobenius_loss_function(parameters::AbstractArray{<:Real, 1}, prec::T_preconditioner, params::optimization_parameters)
    N = 2^prec.L+1
    DOF_x, DOF_y = N - 2 + sum(prec.BCs[1:2] .== "N"), N - 2 + sum(prec.BCs[3:4] .== "N")
    dimension = (DOF_y, DOF_x)
    iteration(x) = x
    if prec.symmetry == "symmetric"
        iteration(x) = x - apply_T_preconditioner(reshape(prec.A*reshape(apply_T_preconditioner(x, parameters, prec), (DOF_y*DOF_x, )), dimension), parameters, prec)
    else
        iteration(x) = x - apply_T_preconditioner(reshape(prec.A*reshape(x, (DOF_y*DOF_x, )), dimension), parameters, prec)
    end
    loss = ρ_F_T(dimension, 1, iteration, params.k_batch)
end

function Frobenius_loss_function(parameters::AbstractArray{<:Real, 1}, prec::preconditioner, params::optimization_parameters)
    M = assemble_matrix(parameters, prec)
    iteration(x) = x
    if prec.symmetry == "left"
        iteration(x) = x - M*(prec.A*x)
    else
        iteration(x) = x - M*(prec.A*(M*x))
    end
    loss = ρ_F(size(prec.A)[1], 1, iteration, params.k_batch)
end

function scalar_loss_function(parameters::AbstractArray{<:Real, 1}, M::SparseMatrixCSC{Float64,Int64}, prec::T_preconditioner, params::optimization_parameters)
    iteration(x) = Richardson(x, prec.A, x -> M*x, parameters[end]; prec=prec.symmetry)
    dimension = size(prec.A)[1]
    if params.error == "Frobenius"
        loss = ρ_F(dimension, params.k_sweeps, iteration, params.k_batch)
    elseif params.error == "Spectral"
        loss = ρ_error(dimension, params.k_sweeps, iteration, params.k_batch)
    end
    return loss
end

function spectral_radius(parameters::AbstractArray{<:Real, 1}, prec::preconditioner, params::optimization_parameters)
    function iteration(x::AbstractArray{<:Real, 2}, M::SparseMatrixCSC{<:Real, <:Integer}, A::SparseMatrixCSC{<:Real, <:Integer}, symmetry::String)
        if symmetry == "symmetric"
            return M*A*M*x
        elseif symmetry == "left"
            return M*A*x
        end
    end
    M = assemble_matrix(parameters, prec)
    dimension = size(prec.A)[1]
    if params.error == "Frobenius"
        loss = ρ_F(dimension, params.k_sweeps, x -> iteration(x, M, prec.A, prec.symmetry), params.k_batch)
    elseif params.error == "Spectral"
        loss = ρ_error(dimension, params.k_sweeps, x -> iteration(x, M, prec.A, prec.symmetry), params.k_batch)
    end
    return loss
end

function spectral_radius(parameters::AbstractArray{<:Real, 1}, prec::T_preconditioner, params::optimization_parameters)
    function iteration(x::AbstractArray{<:Real, 2}, M::SparseMatrixCSC{<:Real, <:Integer}, A::SparseMatrixCSC{<:Real, <:Integer}, symmetry::String)
        if symmetry == "symmetric"
            return M*A*M*x
        elseif symmetry == "left"
            return M*A*x
        end
    end
    M = assemble_matrix(parameters, prec)
    dimension = size(prec.A)[1]
    if params.error == "Frobenius"
        loss = ρ_F(dimension, params.k_sweeps, x -> iteration(x, M, prec.A, prec.symmetry), params.k_batch)
    elseif params.error == "Spectral"
        loss = ρ_error(dimension, params.k_sweeps, x -> iteration(x, M, prec.A, prec.symmetry), params.k_batch)
    end
    return loss
end

function test_ρ_1(k::Integer, N::Integer, L::Integer)
    A, b, x = BVP_DD(L, x -> one(x), x -> zero(x), 1, 1, x -> x)
    ρ = maximum(abs.(eigvals(Array(A))))
    ρ_err = ρ_error(2^L-1, k, x -> A*x)
    ρ_Frob = ρ_F(2^L-1, k, N, x -> A*x)
    println("Spectral radius ", ρ)
    println("Estimation with error ", ρ_err)
    println("Estimation with Frobenius ", ρ_Frob)
end

function test_ρ_2(k::Integer, L::Integer)
    A, b, x = BVP_DD(L, x -> 1 + x^2, x -> exp(x), 1, 1, x -> x)
    ρ = maximum(abs.(eigvals(Array(A))))
    ρ_estimated = ρ_error(2^L-1, k, x -> A*x)
    println("Spectral radius ", ρ)
    println("Estimation ", ρ_estimated)
end

max_stars_repo_name: VLSF/neuralBPX

max_issues_repo_name: VLSF/neuralBPX

max_forks_repo_name: VLSF/neuralBPX


text:
export Poisson_Regression, predict, fit

type poissonRegression
    #array of weights
    w
    #w0 constnat
    w_0::Float64
    #learning rate
    eta::Float64
    #decay rate for the learning rate
    alpha::Float64
    #number of samples parsed through (will be incremnted automatically)
    n::Int64
    #number of features
    d::Int64
end

#constructor for convenience
function Poisson_Regression(;eta::Float64=0.0001, alpha::Float64=1.00, d::Int64=0)
    if d != 0
        model = poissonRegression(zeros(d), 0.0, eta, alpha, 0, d)
    else
        model = poissonRegression(nothing, 0.0, eta, alpha, 0, d)
    end

    #return the model
    model
end

#prediction using n examles (nxd matrix)
function predict(model::poissonRegression, x::Array{Float64,2})
    linear_prediction = x*model.w+model.w_0
    prediction = exp(linear_prediction)
end

#prediction using 1xd array
function predict(model::poissonRegression, x::Array{Float64,1})
    linear_prediction = x*model.w+model.w_0
    prediction = exp(linear_prediction)
end

#takes in nxd batch data as an input, conducts stochastic gradient descent
function fit(model::poissonRegression, y::Array{Float64, 1}, x::Array{Float64, 2})
    #checking if y and x match in size
    if length(y)!=size(x)[1]
        return nothing
    end

    #initiating weight array if necessary
    if model.d == 0
        model.d = size(x)[2]
        model.w = zeros(model.d)
    end

    #updating info (right now its just the number of examples parsed)
    model.n += length(y)
    num_data = length(y)

    #making prediction
    prediction = predict(model, x)

    #updating w_0
    model.w_0 = model.w_0 + model.eta*(sum(y-prediction, 1)[1]*1.0/num_data)

    #updating w
    model.w = model.w + model.eta*((transpose(x)*(y-prediction))/num_data)

    model
end


text:
"""
Data structure for particle filters
- effectiveSampleSize(pc :: ParticleContainer)
- normalise!(pc::ParticleContainer)
- consume(pc::ParticleContainer): return incremental likelihood
"""

typealias Particle Trace

type ParticleContainer{T<:Particle}
  model :: Function
  num_particles :: Int
  vals  :: Array{T,1}
  logWs :: Array{Float64,1}  # Log weights (Trace) or incremental likelihoods (ParticleContainer)
  logE  :: Float64           # Log model evidence
  conditional :: Union{Void,Conditional} # storing parameters, helpful for implementing rejuvenation steps
  n_consume :: Int # helpful for rejuvenation steps, e.g. in SMC2
  ParticleContainer(m::Function,n::Int) = new(m,n,Array{Particle,1}(),Array{Float64,1}(),0.0,nothing,0)
end

call{T}(::Type{ParticleContainer{T}}, m) = ParticleContainer{T}(m, 0)

Base.collect(pc :: ParticleContainer) = pc.vals # prev: Dict, now: Array
Base.length(pc :: ParticleContainer)  = length(pc.vals)
Base.similar(pc :: ParticleContainer) = ParticleContainer(pc.model, 0)
# pc[i] returns the i'th particle
Base.getindex(pc :: ParticleContainer, i :: Real) = pc.vals[i]

# registers a new x-particle in the container
function Base.push!(pc :: ParticleContainer, p :: Particle)
  pc.num_particles += 1
  push!(pc.vals, p)
  push!(pc.logWs, 0)
  pc
end
Base.push!(pc :: ParticleContainer) = Base.push!(pc, eltype(pc.vals)(pc.model))
function Base.push!(pc :: ParticleContainer, n :: Int)
  vals = Array{eltype(pc.vals), 1}(n)
  logWs = Array{eltype(pc.logWs), 1}(n)
  for i=1:n
    vals[i]  = eltype(pc.vals)(pc.model)
    logWs[i] = 0
  end
  append!(pc.vals, vals)
  append!(pc.logWs, logWs)
  pc.num_particles += n
  pc
end

# clears the container but keep params, logweight etc.
function Base.empty!(pc :: ParticleContainer)
  pc.num_particles = 0
  pc.vals  = Array{Particle,1}()
  pc.logWs = Array{Float64,1}()
  pc
end

# clones a theta-particle
function Base.copy(pc :: ParticleContainer)
  particles = collect(pc)
  newpc     = similar(pc)
  for p in particles
    newp = forkc(p)
    push!(newpc, newp)
  end
  newpc.logE        = pc.logE
  newpc.logWs       = deepcopy(pc.logWs)
  newpc.conditional = deepcopy(pc.conditional)
  newpc.n_consume   = pc.n_consume
  newpc
end

# run particle filter for one step, return incremental likelihood
function Base.consume(pc :: ParticleContainer)
  @assert pc.num_particles == length(pc)
  # normalisation factor: 1/N
  _, z1      = weights(pc)
  n = length(pc.vals)

  particles = collect(pc)
  num_done = 0
  for i=1:n
    p = pc.vals[i]
    score = consume(p)
    if isa(score, Real)
      increase_logweight(pc, i, Float64(score))
    elseif score == Val{:done}
      num_done += 1
    else
      "[consume]: error in running particle filter."
    end
  end

  if num_done == length(pc)
    res = Val{:done}
  elseif num_done != 0
    error("[consume]: mis-aligned execution traces, num_particles= $(n), num_done=$(num_done).")
  else
    # update incremental likelihoods
    _, z2      = weights(pc)
    res = increase_logevidence(pc, z2 - z1)
    pc.n_consume += 1
    # res = increase_loglikelihood(pc, z2 - z1)
  end

  res
end

function weights(pc :: ParticleContainer)
  @assert pc.num_particles == length(pc)
  logWs = pc.logWs
  Ws = exp(logWs-maximum(logWs))
  logZ = log(sum(Ws)) + maximum(logWs)
  Ws = Ws ./ sum(Ws)
  return Ws, logZ
end

function effectiveSampleSize(pc :: ParticleContainer)
  Ws, _ = weights(pc)
  ess = sum(Ws) ^ 2 / sum(Ws .^ 2)
end

increase_logweight(pc :: ParticleContainer, t :: Int, logw :: Float64) =
  (pc.logWs[t]  += logw)

increase_logevidence(pc :: ParticleContainer, logw :: Float64) =
  (pc.logE += logw)

function resample!( pc :: ParticleContainer,
                   randcat :: Function = Turing.resampleSystematic,
                   ref :: Union{Particle, Void} = nothing)
  n1, particles = pc.num_particles, collect(pc)
  @assert n1 == length(particles)

  # resample
  Ws, _ = weights(pc)
  n2    = isa(ref, Void) ? n1 : n1-1
  indx  = randcat(Ws, n2)

  # fork particles
  empty!(pc)
  num_children = zeros(Int,n1)
  map(i->num_children[i]+=1, indx)
  for i = 1:n1
    p = particles[i] == ref ? Traces.forkc(particles[i]) : particles[i]
    num_children[i] > 0 && push!(pc, p)
    for k=1:num_children[i]-1
      newp = Traces.forkc(p)
      push!(pc, newp)
    end
  end

  if isa(ref, Particle)
    # Insert the retained particle. This is based on the replaying trick for efficiency
    #  reasons. If we implement PG using task copying, we need to store Nx * T particles!
    push!(pc, ref)
  end

  pc
end


text:
facts("Trivial things are checkable") do

  context("Numbers can be floats") do
      @fact typeof(1.)    => Float64
      @pending typeof(1.) => Float32
  end

  context("Numbers can be ints") do
      @fact typeof(1) => Int
  end
end


text:
function calc_bias_potential(umbrella_center, data)
    spring_constant = 200.0 * (pi/180.0)^2
    return spring_constant * (data - umbrella_center)^2
end


text:
using Plots
pyplot(grid=false)
using DifferentialEquations
using ColorSchemes
using DataFrames
using CSV
using LaTeXStrings
using Statistics
using GLM
using Measures
using QuadGK
using Cairo, Fontconfig
import Gadfly
g = Gadfly
using Compose

fnt = Plots.font("Arial")
default(titlefont=fnt, guidefont=fnt, tickfont=fnt, legendfont=fnt)

################################################################################
# Model:
# A <-> R <-> I
# ka: R -> A
# kr: A -> R
# ki: R -> I
# kp: I -> R
################################################################################

function three_state_model(dc, c, k, t)
    dc[1] = k['a'] * c[2] - k['s'] * c[1]
    dc[2] = k['s'] * c[1] - k['a'] * c[2] - k['i'] * c[2] + k['p'] * c[3]
    dc[3] = k['i'] * c[2] - k['p'] * c[3]
end

function isoutofdomain(c, k, t)
    if (c[1]<=1) && (c[2]<=1) && (c[3]<=1) && (c[1]>=0) && (c[2]>=0) && (c[3]>=0)
        return false
    else
        return true
    end
end

function run_recruitment(k_time_pairs, get_irrev=true)
    t = fill(Float64, 0)
    u = fill(Float64[], 0)

    c0 = [1.0; 0.0; 0.0]

    for pair in k_time_pairs
        k = pair[1]
        times = pair[2]
        prob = ODEProblem(three_state_model, c0, times, k)
        sol = solve(prob, isoutofdomain=isoutofdomain, saveat=0.01)
        t = [t; sol.t]
        u = [u; sol.u]
        c0 = sol.u[end]
    end
    as = [x[1] for x in u]
    rs = [x[2] for x in u]
    is = [x[3] for x in u]
    if get_irrev
        data = [as rs is]
    else
        data = [as rs]
    end

    return [t, data]
end

function plot_recruitment_2state(k_time_pairs, grey_pairs, blue_pairs, title, filename)
    sol = run_recruitment(k_time_pairs, false)
    p1 = Plots.plot()
    for pair in grey_pairs
        p1 = vspan!(p1, pair, label="", color=parse(Colorant, "#F5EBCE"))
    end
    for pair in blue_pairs
        p1 = vspan!(p1, pair, label="", color=parse(Colorant, "#DEE6F0"))
    end

    p1 = plot!(p1,
               sol[1], sol[2],
               label=["Active" "Reversibly Silent"],
               color = [:green :grey],
               title=title,
               xlabel="Days",
               ylabel="Fraction of Cells",
               yticks = 0:0.5:1.0,
               xticks = [0, 5, 10, 15],
               linewidth=2,
               legendfont=12,
               legend=:topright,
               size=(500,250))

    savefig(p1, filename)
end

################################################################################
### FIGURE 1C
################################################################################

function plot_recruitment_2state_1c(k_time_pair_list, grey_pairs, blue_pairs, filename)
    sols = [run_recruitment(k, false) for k in k_time_pair_list]

    p1 = Plots.plot()
    for pair in grey_pairs
        p1 = vspan!(p1, pair, label="", color=parse(Colorant, "#F5EBCE"))
    end
    for pair in blue_pairs
        p1 = vspan!(p1, pair, label="", color=parse(Colorant, "#DEE6F0"))
    end

    colors = ["#009e73", "#0072b2"]
    labels = [L"k_s=10/\mathrm{d}, k_a=0.1/\mathrm{d}", L"k_s=0.1/\mathrm{d}, k_a=10/\mathrm{d}"]
    labels = [L"\mathrm{high} k_s\mathrm{, low} k_a", L"\mathrm{high} k_a\mathrm{, low} k_s"]
    labels = [L"high $k_s$, low $k_a$", L"low $k_s$, high $k_a$"]

    for i in 1:length(sols)
        sol = sols[i]
        p1 = plot!(p1,
                   sol[1], sol[2][:,2],
                   xlabel="Days",
                   ylabel="Fraction Silent",
                   color = colors[i],
                   label = labels[i],
                   xlims = [0, 12],
                   ylims = [0, 1],
                   yticks = 0:0.5:1.0,
                   xticks = [0, 6, 12],
                   linewidth=2,
                   legendfont=12,
                   legend=:bottomright,
                   size=(500,250))
    end

    savefig(p1, filename)
end

k_mem_rec = Dict('a' => 0,
                 's' => 0.1,
                 'i' => 0,
                 'p' => 0)
k_mem_rel = Dict('a' => 10,
                 's' => 0,
                 'i' => 0,
                 'p' => 0)
mem_time_pairs = [[k_mem_rec, (0.0, 6.0)],
                  [k_mem_rel, (6.0, 12.0)]]

k_sil_rec = Dict('a' => 0,
                 's' => 10,
                 'i' => 0,
                 'p' => 0)
k_sil_rel = Dict('a' => 0.1,
                 's' => 0,
                 'i' => 0,
                 'p' => 0)
sil_time_pairs = [[k_sil_rec, (0.0, 6.0)],
                  [k_sil_rel, (6.0, 12.0)]]

time_pairs = [sil_time_pairs, mem_time_pairs]

gp1 = [[0.0, 6.0]]
bp1 = [[6.0, 12.0]]
f1 = "fig_1c.pdf"

plot_recruitment_2state_1c(time_pairs, gp1, bp1, f1)

################################################################################
#### FIGURE 1D
################################################################################

c1 = parse(Colorant, "#1b9e77")
c2 = parse(Colorant, "#7570b3")
c3 = parse(Colorant, "#d95f02")

function cell_div_model(dc, c, k, t)
    ks = k['s']
    ka = k['a']
    kd = k['d']

    # number of on cells grows as gr, plus any reactivated off cells
    dc[1] = kd * c[1] + c[2] * ka - c[1] * ks
    # number of off cells grows as g, plus any silenced on cellsr
    dc[2] = kd * c[2] - c[2] * ka + c[1] * ks
end

function run_celldiv_model(ks, ka, kd, a0, tlim)
    c0 = [a0; 0.0]

    println(c0)

    k = Dict('s' => ks,
             'a' => ka,
             'd' => kd)

    times = [0.0, tlim]
    prob = ODEProblem(cell_div_model, c0, times, k)
    sol = solve(prob, saveat=0.01)
    t = sol.t
    u = sol.u
    ons = [x[1] for x in u]
    offs = [x[2] for x in u]
    data = [offs ons]
    return [t, data]
end

function plot_celldiv(ks, ka, kd, a0, tlim)
    t, data = run_celldiv_model(ks, ka, kd, a0, tlim)
    p1 = Plots.plot()
    p1 = plot!(p1,
               t, data,
               label = ["Reversibly Silent" "Active"],
               color = [c3 c2],
               xlabel = "Days",
               ylabel = "Number of Cells",
               ylims = [0, 700],
               yticks = [0, 350, 700],
               xlims = [0, 2],
               xticks = [0, 1, 2],
               linewidth = 2,
               legendfont = 12,
               size = (450, 300))
    return p1
end

p = plot_celldiv(1.0, 0.0, 1.0, 100, 2)
savefig(p, "fig_1d.pdf")

################################################################################
#### FIGURE 1E
################################################################################

signals = zeros(0)

function sinusoid(dc, c, k, t)
    dc[1] = k['a'] * c[2] - k['s'] * c[1]
    dc[2] = k['s'] * c[1] - k['a'] * c[2]
end

function isoutofdomain_sinus(c, k, t)
    for i in c
        if (i > 1) || (i < 0)
            return true
        end
    end
    return false
end

function run_sinus(k_time_pairs, get_irrev=true)
    t = fill(Float64, 0)
    u = fill(Float64[], 0)

    c0 = [0.75; 0.25]

    for pair in k_time_pairs
        k = pair[1]
        times = pair[2]
        prob = ODEProblem(sinusoid, c0, times, k)
        sol = solve(prob, isoutofdomain=isoutofdomain_sinus, saveat=0.001)
        t = [t; sol.t]
        u = [u; sol.u]
        c0 = sol.u[end]
    end
    as = [x[1] for x in u]
    rs = [x[2] for x in u]

    data = [as rs]
    return [t, data]
end

function get_plot_sinus(k_time_pairs, grey_pairs, blue_pairs, duration)
    sol = run_sinus(k_time_pairs)
    tlim = duration

    p1 = Plots.plot()

    crec = parse(Colorant, "#F5EBCE")
    crel = parse(Colorant, "#DEE6F0")

    for pair in grey_pairs
        ymin = 0.005
        ymax = 1.05
        p1 = plot!(p1,
                   [pair[1], pair[1], pair[2], pair[2]],
                   [ymin, ymax, ymax, ymin],
                   color = crec,
                   linealpha = 0,
                   label = "",
                   seriestype = :shape)
        #  p1 = vspan!(p1, pair, label="", linecolor=crec, fillcolor=crec, linealpha=0, fillalpha=0)
    end
    for pair in blue_pairs
        ymin = 0.005
        ymax = 1.05
        p1 = plot!(p1,
                   [pair[1], pair[1], pair[2], pair[2]],
                   [ymin, ymax, ymax, ymin],
                   color = crel,
                   linealpha = 0,
                   label = "",
                   seriestype = :shape)
        #  p1 = vspan!(p1, pair, label="", linecolor=crel, fillcolor=crel, linealpha=0)
    end

    p1 = plot!(p1,
               sol[1], sol[2][:,2],
               order = 10,
               label="",
               legend=:bottomright,
               xlabel="Days",
               #  xticks= [0, 1, 2, 3, 5, 10],
               xticks = [0, tlim/2.0, tlim],
               xlims = (0, tlim),
               ylabel="Fraction Silent",
               yticks = [0, 0.5, 1],
               ylims = (0, 1.05),
               linewidth=2,
               color=:green,
               bottom_margin=0.25cm,
               right_margin=0.25cm,
               titlefontfamily = "Arial",
               size=(450,125))

    #  tf2_vals = sol[2][:, 1]
    #  tf2_vals_trunc = tf2_vals[Int(floor(0.5*length(tf2_vals))):length(tf2_vals)]

    return p1
end

function get_sinus_timepairs(ka, ks, period, duration)
    k_rec = Dict('a' => ka,
                 's' => ks)
    k_rel = Dict('a' => ka,
                 's' => 0)
    l = [[[k_rec, (i, i+period)], [k_rel, (i+period, i+2*period)]] for i in 0:2*period:duration]
    return [l[i][j] for i in 1:length(l) for j in 1:length(l[i])]
end

function get_greypairs(period, duration)
    return [[i, i + period] for i in 0:2*period:duration]
end

function get_bluepairs(period, duration)
    return [[i, i + period] for i in period:2*period:duration]
end

function get_plotpairs(period, ks)
    duration = 12
    ka = 1.0
    return [get_sinus_timepairs(ka, ks, period, duration),
            get_greypairs(period, duration),
            get_bluepairs(period, duration),
            duration]
end

palette = ColorSchemes.Spectral_4
palette = [colorant"#e66101", colorant"#fdb863", colorant"#b2abd2", colorant"#5e3c99"]
palette = ["grey25", "brown4", "blue", "darkviolet"]
palette = [colorant"skyblue", colorant"purple", colorant"green", colorant"darkorange"]
palette = [colorant"grey25", colorant"darkorange", colorant"purple", colorant"green"]
p1 = get_plot_sinus(get_plotpairs(3.0, 1.0)...)
p1 = plot!(p1, title=L"Long Period, Low $k_s$", titlefontcolor = palette[1], titleweight = "bold")

p2 = get_plot_sinus(get_plotpairs(0.25, 1.0)...)
p2 = plot!(p2, title=L"Short Period, Low $k_s$", titlefontcolor = palette[2])

p3 = get_plot_sinus(get_plotpairs(3.0, 3.0)...)
p3 = plot!(p3, title=L"Long Period, High $k_s$", titlefontcolor = palette[3])

p4 = get_plot_sinus(get_plotpairs(0.25, 3.0)...)
p4 = plot!(p4, title=L"Short Period, High $k_s$", titlefontcolor = palette[4])

lay = @layout [a b ; c d]
p = plot(p4, p3, p2, p1, layout=lay, size=(600, 450), linewidth=2)
savefig(p, "fig_1e.pdf")

################################################################################
#### FIGURE 1F
################################################################################

function bot_val(t, ks, ka)
    num = -exp(ka*t)*ka + exp(t*(-ka-ks))*ka + exp(t*(-ka-ks))*ks - exp(ka*t + t*(-ka - ks))*ks
    den = (exp(ka*t) - exp(t*(-ka - ks)))*(ka + ks)
    # den = (exp(t*(2*ka+ks)) - 1)*(ka + ks)
    return -num/den
end

function top_val(t, ks, ka)
    num = exp(ka*t)*ka - exp(t*(-ka-ks))*ka - ks + exp(ka*t)*ks
    den = (exp(ka*t) - exp(t*(-ka - ks)))*(ka + ks)
    return num/den
end

t_vals = zeros(0)
k_vals = zeros(0)
d_vals = zeros(0)

for t in 0:0.01:10
    for k in 0:0.1:10
        d = top_val(t/2.0, k, 1.0) - bot_val(t/2.0, k, 1.0)
        if t == 0
            d = 0.0
        end
        push!(t_vals, t * 2)
        push!(k_vals, k)
        push!(d_vals, d)
    end
end

data = DataFrame(Dict("Period of Oscillation" => t_vals,
                      "Rate of Silencing (1/days)" => k_vals,
                      "Amplitude" => d_vals))

function get_pt_annotation(x, y, color="red")
    return g.Guide.annotation(compose(context(),
                                      Compose.circle(x, y, 0.24),
                                      Compose.stroke("white"),
                                      Compose.linewidth(0.5),
                                      fill(color),
                                     ))
end
palette = [colorant"#e66101", colorant"#fdb863", colorant"#b2abd2", colorant"#5e3c99"]
palette = ["grey25", "brown4", "blue", "darkviolet"]
palette = [colorant"skyblue", colorant"purple", colorant"green", colorant"darkorange"]
palette = [colorant"grey25", colorant"darkorange", colorant"green", colorant"purple"]

p = g.plot(data, x="Period of Oscillation", y="Rate of Silencing (1/days)", color="Amplitude",
           g.Geom.rectbin, g.Guide.xticks(ticks=collect(0:2:10)), g.Guide.yticks(ticks=collect(0:2:10)),
           g.Coord.cartesian(xmin=0, xmax=10, ymin=0, ymax=10, aspect_ratio=1),
           g.Scale.ContinuousColorScale(p -> get(ColorSchemes.viridis, p)),
           get_pt_annotation(6.0, 1.0, palette[1]),
           get_pt_annotation(0.5, 1.0, palette[2]),
           get_pt_annotation(6.0, 3.0, palette[4]),
           get_pt_annotation(0.5, 3.0, palette[3]),
           g.Theme(major_label_font="ArialMT", minor_label_font="ArialMT",
                   key_title_font="ArialMT", key_label_font="ArialMT", key_label_color="black",
                   grid_color="black",
                   minor_label_color="black", major_label_color="black", key_title_color="black"))

img = PDF("fig_1f.pdf")
draw(img, p)


text:
using OrderedCollections: OrderedDict

function ordered_dict_show(dict, padchar)
    if length(dict) == 0
        return "$(typeof(dict).name) with no entries"
    elseif length(dict) == 1
        return "$(typeof(dict).name) with 1 entry:\n" *
               "$padchar   └── $(dict.keys[1]) => $(typeof(dict.vals[1]).name)"
    else
        return string(typeof(dict).name, " with $(length(dict)) entries:\n",
                      Tuple("$padchar   ├── $name => $(typeof(dict[name]).name)\n" for name in dict.keys[1:end-1])...,
                            "$padchar   └── $(dict.keys[end]) => $(typeof(dict.vals[end]).name)")
    end
end


text:
include("euclidean-1-NN.jl")
include("../algoritmos-busqueda/busqueda-local.jl")
include("../utils/funcion-objetivo.jl")

"""
    BL_LearnerOneNN(data::Matrix{<:Real}, labels::Vector{<:Real})
Devuelve un clasificado 1-NN refinado por unos pesos en la distancia euclídea.
Ese vector de pesos se ha obtenido a partir de un algoritmo de búsqueda local
"""
function BL_LearnerOneNN(data::Matrix{<:Real}, labels)
    # Creamos función objetivo
    a = 0.5
    umbral_tasa_reduccion = 0.1
    F = CrearFuncionObjetivo(data, labels, a, umbral_tasa_reduccion )

    numero_atributos = size(data)[2] # los atributos son las columnas
    # Criterios de acorde a los requisitos
    numero_maximo_vecinos_sin_mejora = 20 * numero_atributos
    nuero_maximo_evaluaciones = 15000
    # Buscamos pesos de acorde a parámetros calculados
    w = PrimeroMejor(numero_maximo_vecinos_sin_mejora,
    nuero_maximo_evaluaciones,numero_atributos,umbral_tasa_reduccion, F)

    return WeightedLearnerEuclideanOneNN(w, data, labels), F(w), w
end


text:
using AtomicLevels: CSF, orbital_from_string

p_suffix(p::Parity) = isodd(p) ? "o" : ""
p_suffix(cfg::Configuration) = p_suffix(parity(cfg))

function parse_csf(filename)
    # It is assumed that all term symbols are on the form ML (and MLS
    # for intermediate terms) where M is multiplicity, and S is
    # seniority number.
    ref_csfs = NonRelativisticCSF[]
    open(filename) do file
        readline(file)
        core_cfg = close(fill(parse(Configuration{Orbital}, join(split(readline(file)), " "))))
        while true
            peel_cfg = readline(file)
            peel_cfg[1] == '*' && break
            peel_cfg = parse(Configuration{Orbital},
                             join(split(replace(replace(replace(peel_cfg, "( "=>""), "("=>""), ")"=>"")), " "))
            cfg = core_cfg + peel_cfg
            ts = split(readline(file))
            np = length(peel_cfg)
            its = map(enumerate(ts[1:np])) do (i,t)
                ip = p_suffix(parity(Configuration(peel_cfg[i]...)))
                IntermediateTerm(parse(Term, "$(t[1:2])$(ip)"), parse(Int, t[end]))
            end
            coupled_terms = map(enumerate(ts[vcat(1,np+1:end)])) do (i,t)
                parse(Term, "$(t[1:2])$(p_suffix(peel_cfg[1:i]))")
            end
            push!(ref_csfs, CSF(cfg, its, coupled_terms))
        end
    end
    ref_csfs
end


text:
#*******************************************************************************
#
#
#                Aplicação do split-operator em um sistema 2D
#               com condições de contorno tipo hardwall.
#
#
#*******************************************************************************

using PyPlot
using PyCall

@pyimport matplotlib.animation as anim # necessário para gerar vídeo da simulação

PyPlot.matplotlib[:rc]("text", usetex=true) # allow tex rendering
PyPlot.matplotlib[:rc]("font", family="serif")

# constantes físicas
Ry = 13.6056917253e3; # meV
a0 = 0.0529177208319; # nm
h2m = (Ry*a0^2); #  massa efetiva GaAs
hbar = 0.65821188926; # meV.ps

################################################################################
                    # Grid Bidimensional (Espaço real)
################################################################################
Lx = 1000; # Tamanho da dimensão "x" em nanometros (nm)
Ly = 200; # Tamanho da dimensão "y" em nanometros (nm)
Nx = div(Lx,3) # número de pontos na dimensão x
Ny = div(Ly,1) # número de pontos na dimensão y
# Nx = 2*Lx;
# Ny = 2*Ly;

x = linspace(-Lx/2, Lx/2, Nx+1)[1:end-1];
y = linspace(-Ly/2, Ly/2, Ny+1)[1:end-1];
dx = x[2] - x[1];
dy = y[2] - y[1];
X = kron(x', ones(Ny))
Y = kron(ones(Nx)', y)

################################################################################
                    # Espaço-k (na dimensão longitudinal à tira)
################################################################################
# A array "kx" deve ser 'invertida' via 'fftshift' para que todas as
# as funções de 'kx' possam ser multiplicadas elementwise sem confusão
# quando estivermos trabalhando com trasformadas de Fourier.

# Se houver interesse em observar como transformadas de Fourier e funções de
# 'kx' é interessante antes do 'plotting' reverter o 'shift'.

kx = fftshift(linspace(-pi/dx, pi/dx, Nx+1)[1:end]); # Valores possíveis de k_x
k0_x = pi/(2*dx);
k0_y = pi/(2*dy);

################################################################################
                    # Condição Inicial: Pacote Gaussiano
################################################################################

# Note que a condição inicial está representada no espaço de posição.

x_0 = -50; #  centro do pacote
y_0 = 0; #  centro do pacote
sigma_x = 30.0; # largura longitudinal
sigma_y = 30.0; # largura transversal
k0_x *= 0.5; # fração do kx_max
k0_y *= 0;   # fração do ky_max
# Norm = ((pi^2)*(sigma_x*sigma_y)^2)^(-1/4); # Fator para normalização
Norm = 1 ; # Fator para pacote NÃO normalizado

# pacote de ondas na representação de posições:
Ψ_0 = Norm * exp(-(X-x_0).^2/(2*sigma_x^2)-(Y-y_0).^2/(2*sigma_y^2)
                    + 1im*(k0_x*X + k0_y*Y))

################################################################################
                        # Hamiltoniano
################################################################################

#  O operador Hamiltoniano aqui é simplesmente aquele para uma partícula
# (Elétron) submetida a um potencial dependente da posição.

# Potencial:
α, β, γ = 10., .1, 50.;
V_degrau = α./(e.^(-β*(X-γ))+1);

# Cinética:
# matrizes Tridiagonais (Ny x Ny)
d2dy2 = 1/(dy^2) * SymTridiagonal(-2*ones(Ny), ones(Ny-1))

T_y  = - h2m * full(d2dy2)
T_kx =   h2m * kx.^2

################################################################################
                        # Funções para evolução temporal
################################################################################

function plot_simples(Ψ)
  pcolormesh(x,y,abs(Ψ).^2)
  ylim(-200,200)
  title("Gaussian distribution");
  xlabel(L"$x$");
  ylabel(L"$y$");
  return 0
end

function psi_dt(Ψ, Δt)
  # espaço de posições 'x'
  Ψ = [ exp(-1im*V_degrau[j,i]*Δt/(2*hbar)) * Ψ[j,i] for j=1:Ny, i=1:Nx ];

  # representação mista (kx,y)
  fft!(Ψ,2);
  for i = 1:Nx
    Ψ[:,i] = expm(-1im*(T_kx[i]*eye(Ny) + T_y)*Δt/hbar) * Ψ[:,i] ;
  end

  # espaço de posições 'x'
  ifft!(Ψ,2);
  Ψ = [exp(-1im*V_degrau[j,i]*Δt/(2*hbar)) * Ψ[j,i] for j=1:Ny, i=1:Nx ];

  return Ψ;
end

function psi_dt_animado(Δt)
  global Ψ_0
  if Δt == 0
    # Essa  condição existe para incluir a
    # cond. inicial no video da simulação
    clf()
    plot_simples(Ψ_0)
    return 0
  else
    # espaço de posições 'x'
    Ψ_0 = [ exp(-1im*V_degrau[j,i]*Δt/(2*hbar)) * Ψ_0[j,i] for j=1:Ny, i=1:Nx ];

    # representação mista (kx,y)
    fft!(Ψ_0,2);
    for i = 1:Nx
      Ψ_0[:,i] = expm(-1im*(T_kx[i]*eye(Ny) + T_y)*Δt/hbar) * Ψ_0[:,i] ;
    end

    # espaço de posições 'x'
    ifft!(Ψ_0,2);
    Ψ_0 = [exp(-1im*V_degrau[j,i]*Δt/(2*hbar)) * Ψ_0[j,i] for j=1:Ny, i=1:Nx ];

    clf()
    plot_simples(Ψ_0)
    return 0;
  end
end

function psit(Ψ, t_final, N_steps = 20)
  tspan = linspace(0, t_final, N_steps-1);
  delta_t = tspan[2]-tspan[1];

  for i = 1:N_steps
    psi_aux = psi_dt(Ψ,delta_t);
    Ψ = psi_aux;
  end
  Ψ
end

function psit_livre(Ψ, t)
  fft!(Ψ,2); # substitui Ψ por sua transformada de Fourier em "x"
  for i = 1:Nx
    Ψ[:,i] = expm(-1im*(T_kx[i]*eye(Ny) + T_y)* t/hbar)*Ψ[:,i] ;
  end
  ifft!(Ψ,2)
end

function simula_potencial(;t_final = 2.0, N_steps = 10, intervalo_frames = 100)

  Δt = t_final/N_steps;
  array_of_deltas = Δt*[0; ones(N_steps)];

  fig = figure();

  video = anim.FuncAnimation(fig,
                            psi_dt_animado,
                            frames = array_of_deltas,
                            interval = intervalo_frames);

  path = "/home/marcos/Dropbox/projetos/julia/split operator/";
  nome = string("simulação_com_potencial_degrau_em_x_50.mp4")

  video[:save](string(path,nome), extra_args=["-vcodec", "libx264"]);
  return 0
end

################################################################################
                        # Apresentação dos resultados
################################################################################

Ψ_1 = psit(Ψ_0, 2.0, 6)
# Ψ_1 = psit_livre(Ψ_0, 2)

figure(1)
pcolormesh(x,y,abs(Ψ_1).^2)
ylim(-200,200)
title("Gaussian distribution");
xlabel(L"$x$");
ylabel(L"$y$");

figure(2)
pcolormesh(X,Y,V_degrau)
ylim(-200,200)
title("Potential function");
xlabel(L"$x$");
ylabel(L"$y$");
colorbar()

# simula_potencial(t_final=10.0,N_steps=100,intervalo_frames=50)

max_stars_repo_name: mhlmedeiros/Split-Operator

max_issues_repo_name: mhlmedeiros/Split-Operator

max_forks_repo_name: mhlmedeiros/Split-Operator


text:
"""
$(SIGNATURES)
Request URLs for data blob download.

Args:
  navAbilityClient (NavAbilityClient): The NavAbility client.
  userId (String): The userId with access to the data.
  fileId (String): The unique file identifier of the data blob.
"""
function createDownloadEvent(
    navAbilityClient::NavAbilityClient,
    userId::AbstractString,
    fileId::AbstractString
  )
  #
  response = navAbilityClient.mutate(MutationOptions(
      "sdk_url_createdownload",
      GQL_CREATEDOWNLOAD,
      Dict(
          "userId" => userId,
          "fileId" => fileId
      )
  )) |> fetch
  rootData = JSON.parse(response.Data)
  if haskey(rootData, "errors")
      throw("Error: $(rootData["errors"])")
  end
  data = get(rootData,"data",nothing)
  if data === nothing return "Error" end
  urlMsg = get(data,"url","Error")
  # TODO: What about marshalling?
  return urlMsg
end

createDownload(w...) = @async createDownloadEvent(w...)

##

function getDataEvent(
    client::NavAbilityClient,
    userId::AbstractString,
    fileId::AbstractString
  )
  #
  url = createDownload(client, userId, fileId) |> fetch
  io = PipeBuffer()
  Downloads.download(url, io)
  io
end

getDataEvent(client::NavAbilityClient, context::Client, w...) = getDataEvent(client, context.userId, w...)
getData(w...) = @async getDataEvent(w...)

## =========================================================================
## Upload
## =========================================================================

"""
$(SIGNATURES)
Request URLs for data blob upload.

Args:
  navAbilityClient (NavAbilityClient): The NavAbility client.
  filename (String): file/blob name.
  filesize (Int): total number of bytes to upload.
  parts (Int): Split upload into multiple blob parts, FIXME currently only supports parts=1.
"""
function createUploadEvent(
    navAbilityClient::NavAbilityClient,
    filename::AbstractString,
    filesize::Int,
    parts::Int=1
  )
  #
  response = navAbilityClient.mutate(MutationOptions(
    "sdk_url_createupload",
    GQL_CREATEUPLOAD,
    Dict(
      "filename" => filename,
      "filesize" => filesize,
      "parts" => parts
    )
  )) |> fetch
  rootData = JSON.parse(response.Data)
  if haskey(rootData, "errors")
    throw("Error: $(rootData["errors"])")
  end
  data = get(rootData,"data",nothing)
  if data === nothing return "Error" end
  uploadResp = get(data,"createUpload","Error")
  # TODO: What about marshalling?
  return uploadResp
end

createUpload(w...) = @async createUploadEvent(w...)

## Complete the upload

function completeUploadSingleEvent(
    navAbilityClient::NavAbilityClient,
    fileId::AbstractString,
    uploadId::AbstractString,
    eTag::AbstractString,
  )
  response = navAbilityClient.mutate(MutationOptions(
    "completeUpload",
    GQL_COMPLETEUPLOAD_SINGLE,
    Dict(
      "fileId" => fileId,
      "uploadId" => uploadId,
      "eTag" => eTag
    )
  )) |> fetch
  rootData = JSON.parse(response.Data)
  if haskey(rootData, "errors")
    throw("Error: $(rootData["errors"])")
  end
  data = get(rootData,"data",nothing)
  if data === nothing return "Error" end
  uploadResp = get(data,"completeUpload","Error")
  return uploadResp
end

completeUploadSingle(w...) = @async completeUploadSingleEvent(w...)

##

function addDataEvent(
    client::NavAbilityClient,
    blobname::AbstractString,
    blob::AbstractVector{UInt8}
  )
  #
  io = IOBuffer(blob)

  filesize = io.size
  np = 1
  # create the upload url destination
  d = NVA.createUploadEvent(client, blobname, filesize, np)

  url = d["parts"][1]["url"]
  uploadId = d["uploadId"]
  fileId = d["file"]["id"]

  # custom header for pushing the file up
  headers = [
    "Content-Length" => filesize,
    "Accept" => "application/json, text/plain, */*",
    "Accept-Encoding" => "gzip, deflate, br",
    "Sec-Fetch-Dest" => "empty",
    "Sec-Fetch-Mode" => "cors",
    "Sec-Fetch-Site" => "cross-site",
    "Sec-GPC" => 1,
    "Connection" => "keep-alive"
  ]
  #

  # fid = open(filepath,"r")
  resp = HTTP.put(url, headers, io)
  # close(fid)

  # Extract eTag
  eTag = match(r"[a-zA-Z0-9]+",resp["eTag"]).match

  # close out the upload
  res = NVA.completeUploadSingleEvent(client, fileId, uploadId, eTag)

  res == "Accepted" ? nothing : @error("Unable to upload blob, $res")

  fileId
end

addData(w...) = @async addDataEvent(w...)

##

function addDataEntryEvent(
    navAbilityClient::NavAbilityClient,
    userId::AbstractString,
    robotId::AbstractString,
    sessionId::AbstractString,
    variableLabel::AbstractString,
    dataId::AbstractString,
    dataLabel::AbstractString,
    mimeType::AbstractString="",
  )
  response = navAbilityClient.mutate(MutationOptions(
    "sdk_adddataentry",
    GQL_ADDDATAENTRY,
    Dict(
      "userId" => userId,
      "robotId" => robotId,
      "sessionId" => sessionId,
      "variableLabel" => variableLabel,
      "dataId" => dataId,
      "dataLabel" => dataLabel,
      "mimeType" => mimeType,
    )
  )) |> fetch
  rootData = JSON.parse(response.Data)
  if haskey(rootData, "errors")
    throw("Error: $(rootData["errors"])")
  end
  data = get(rootData,"data",nothing)
  if data === nothing return "Error" end
  addentryresp = get(data,"addDataEntry","Error")
  return addentryresp
end

addDataEntryEvent(client::NavAbilityClient,
                  context::Client,
                  w...) = addDataEntryEvent(client,
                                            context.userId,
                                            context.robotId,
                                            context.sessionId,
                                            w...)
#

addDataEntry(w...) = @async addDataEntryEvent(w...)

##

function listDataEntriesEvent(
    navAbilityClient::NavAbilityClient,
    userId::AbstractString,
    robotId::AbstractString,
    sessionId::AbstractString,
    variableLabel::AbstractString
  )
  #
  response = navAbilityClient.mutate(MutationOptions(
    "sdk_listdataentries",
    GQL_LISTDATAENTRIES,
    Dict(
      "userId" => userId,
      "robotId" => robotId,
      "sessionId" => sessionId,
      "variableLabel" => variableLabel
    )
  )) |> fetch
  rootData = JSON.parse(response.Data)
  if haskey(rootData, "errors")
    throw("Error: $(rootData["errors"])")
  end
  data = get(rootData,"data",nothing)
  if data === nothing return "Error" end

  listdata = data["users"][1]["robots"][1]["sessions"][1]["variables"][1]["data"]
  ret = []
  for d in listdata
    tupk = Tuple(Symbol.(keys(d)))
    nt = NamedTuple{tupk}( values(d) )
    push!(ret,
      nt
    )
  end

  return ret
end

listDataEntriesEvent(client::NavAbilityClient,
                      context::Client,
                      w...) = listDataEntriesEvent(client,
                                                    context.userId,
                                                    context.robotId,
                                                    context.sessionId,
                                                    w...)
#

listDataEntries(w...) = @async listDataEntriesEvent(w...)

##


text:
pushfirst!(LOAD_PATH, @__DIR__)
pushfirst!(LOAD_PATH, joinpath(@__DIR__, ".."))
using Documenter, MATDaemon

makedocs(;
    modules = [MATDaemon],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
    ),
    sitename = "MATDaemon.jl",
    authors = "Jonathan Doucette",
    pages = [
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo = "github.com/jondeuce/MATDaemon.jl.git",
    push_preview = true,
    deploy_config = Documenter.GitHubActions(),
)


text:
include("abstract_action_selector.jl")
include("epsilon_greedy_selector.jl")
include("ucb_selector.jl")
include("weighted_selector.jl")
include("alternate_selector.jl")


text:
# This file is a part of SimilaritySearch.jl

export IHCSearch

# Iterated Hill Climbing Search

"""
    IHCSearch(; restarts=length(hints), localimprovements=false)
    IHCSearch(restarts::Integer=20; localimprovements=false)

IHCSearch is an iterated hill climbing algorithma, a local search algorithm. It greedily navigates the search graph
and restart the search `restarts` times.
Multithreading applications must have copies of this object due to shared cache objects.

- `restarts`: The number of restarts.
- `localimprovements`: An experimental technique that if it is true it will achieve very high quality results, at cost of increasing searching time.
"""
@with_kw mutable struct IHCSearch <: LocalSearchAlgorithm
    restarts::Int32 = 16
    localimprovements::Bool = false
end

Base.copy(ihc::IHCSearch; restarts=ihc.restarts, localimprovements=ihc.localimprovements) =
    IHCSearch(; restarts, localimprovements)

function Base.copy!(dst::IHCSearch, src::IHCSearch)
    dst.restarts = src.restarts
    dst.localimprovements = src.localimprovements
end

"""
    hill_climbing(isearch::IHCSearch, index::SearchGraph, q, res::KnnResult, nodeID::Integer)

Runs a single hill climbing search process starting in vertex `nodeID`
"""
function hill_climbing(isearch::IHCSearch, index::SearchGraph, q, res::KnnResult, nodeID::Integer, vstate)
    omin::Int32 = -1
    dmin::Float32 = typemax(Float32)
    localimprovements = isearch.localimprovements
    while true
        dmin = typemax(Float32)
        omin = -1
        @inbounds for childID in keys(index.links[nodeID])
            visited(vstate, childID) && continue
            visit!(vstate, childID)
            d = convert(Float32, evaluate(index.dist, index[childID], q))
            if localimprovements  ## this yields to better quality but has no early stopping
                push!(res, childID, d)
                if d < dmin
                    dmin = d
                    omin = childID
                end
            else
                if push!(res, childID, d) && d < dmin
                    dmin = d
                    omin = childID
                end
            end
        end

        if omin < 0
            break
        else
            nodeID = omin
        end
    end
end

function searchat(isearch::IHCSearch, index::SearchGraph, q, res, startpoint, vstate)
    if !visited(vstate, startpoint)
        visit!(vstate, startpoint)
        d = convert(Float32, evaluate(index.dist, q, index[startpoint]))
        push!(res, startpoint, d)
        hill_climbing(isearch, index, q, res, startpoint, vstate)
    end
end

"""
    search(isearch::IHCSearch, index::SearchGraph, q, res::KnnResult, hints, vstate)

Performs an iterated hill climbing search for `q`.
"""
function search(isearch::IHCSearch, index::SearchGraph, q, res::KnnResult, hints, vstate)
    n = length(index)

    for startpoint in hints
        searchat(isearch, index, q, res, startpoint, vstate)
    end

    if length(res) == 0
        _range = 1:n
        for i in 1:min(isearch.restarts, n)
            startpoint = rand(_range)
            searchat(isearch, index, q, res, startpoint, vstate)
        end
    end

    res
end


text:
module DataKnots4Postgres

import DataKnots:
    AbstractShape,
    BlockOf,
    BlockVector,
    DataKnot,
    DataShape,
    IsLabeled,
    Pipeline,
    Runtime,
    TableCell,
    TupleOf,
    TupleVector,
    ValueOf,
    as_tuples,
    block_cardinality,
    chain_of,
    column,
    cover,
    designate,
    fits,
    flatten,
    lookup,
    quoteof,
    quoteof_inner,
    render_cell,
    shapeof,
    syntaxof,
    target,
    with_elements,
    x0to1,
    x1to1,
    x1toN

using Tables

using LibPQ:
    _DEFAULT_TYPE_MAP,
    PQ_SYSTEM_TYPES,
    Connection,
    execute

using PostgresCatalog:
    PGCatalog,
    PGColumn,
    PGForeignKey,
    PGSchema,
    PGTable,
    PGType,
    introspect

struct Options
end

struct Handle{E}
    conn::Connection
    ety::E
    opt::Options
end

Handle(conn, ety::E) where {E} =
    Handle{E}(conn, ety, Options())

Handle(conn) =
    Handle(conn, introspect(conn))

struct VectorWithHandle{H,V,T} <: AbstractVector{T}
    h::H
    v::V
end

VectorWithHandle(h::H, v::V) where {H,T,V<:AbstractVector{T}} =
    VectorWithHandle{H,V,T}(h, v)

handle(hv::VectorWithHandle) =
    hv.h

output(hv::VectorWithHandle) =
    hv.v

Base.getindex(hv::VectorWithHandle) =
    hv.v

Base.getindex(hv::VectorWithHandle, k::Integer) =
    hv.v[k]

Base.getindex(hv::VectorWithHandle, ks::AbstractVector) =
    VectorWithHandle(hv.h, hv.v[ks])

Base.size(hv::VectorWithHandle) =
    size(hv.v)

Base.eltype(hv::VectorWithHandle) =
    eltype(hv.v)

Base.IndexStyle(::Type{<:VectorWithHandle{H,V}}) where {H,V<:AbstractVector} =
    IndexStyle(V)

struct EntityShape{E} <: DataShape
    ety::E
    opt::Options
    out::AbstractShape
end

entity(shp::EntityShape) = shp.ety

options(shp::EntityShape) = shp.opt

output(shp::EntityShape) = shp.out

replace_output(shp::EntityShape, f) =
    EntityShape(shp.ety, shp.opt, f isa AbstractShape ? f : f(shp.out))

Base.eltype(shp::EntityShape) =
    eltype(shp.out)

quoteof(shp::EntityShape) =
    Expr(:call, nameof(EntityShape), quoteof_inner(shp.ety), quoteof_inner(shp.opt), quoteof_inner(shp.out))

syntaxof(shp::EntityShape) =
    Symbol(string(shp.ety))

shapeof(hv::VectorWithHandle) =
    EntityShape(hv.h.ety, hv.h.opt, shapeof(hv.v))

fits(shp1::EntityShape, shp2::EntityShape) =
    shp1.ety == shp2.ety && shp1.opt == shp2.opt && fits(shp1.out, shp2.out)

DataKnot(conn::Connection) =
    DataKnot(Any, VectorWithHandle(Handle(conn), TupleVector(1)), x1to1)

render_cell(shp::EntityShape, vals::VectorWithHandle, idx::Int, avail::Int, depth::Int=0) =
    if fits(output(shp), TupleOf())
        TableCell(string(entity(shp)))
    else
        render_cell(output(shp), output(vals), idx, avail, depth)
    end

output() = Pipeline(output)

output(::Runtime, input) =
    output(input)

postgres_name(name::AbstractString) =
    "\"$(replace(name, "\"" => "\"\""))\""

postgres_name(qname::Tuple) =
    join(postgres_name.(qname), '.')

postgres_name(names::AbstractVector) =
    join(postgres_name.(names), ", ")

postgres_name(name::Integer) =
    "\$$name"

load_postgres_table(tbl_name, col_names, col_types) =
    Pipeline(load_postgres_table, tbl_name, col_names, col_types)

load_postgres_table(tbl_name, col_names, col_types, icol_names) =
    Pipeline(load_postgres_table, tbl_name, col_names, col_types, icol_names)

function load_postgres_table(::Runtime, input::AbstractVector, tbl_name, col_names, col_types, icol_names=String[])
    @assert input isa VectorWithHandle
    sql = "SELECT $(postgres_name(col_names)) FROM $(postgres_name(tbl_name))"
    if !isempty(icol_names)
        condition = join(["$(postgres_name(icol_name)) = $(postgres_name(i))"
                          for (i, icol_name) in enumerate(icol_names)], " AND ")
        sql = "$sql WHERE $condition"
    end
    results = []
    for row in input
        res = execute(handle(input).conn, sql, values(row))
        push!(results, (length(res), columntable(res)))
    end
    offs = Vector{Int}(undef, length(results)+1)
    offs[1] = top = 1
    lbls = Symbol.(col_names)
    cols = AbstractVector[col_type[] for col_type in col_types]
    for k = eachindex(results)
        sz, res = results[k]
        top += sz
        offs[k+1] = top
        for j = 1:length(cols)
            append!(cols[j], res[j])
        end
    end
    tv = TupleVector(lbls, top-1, cols)
    h = handle(input)
    ety = h.ety
    ety′ = get_catalog(ety)[tbl_name[1]][tbl_name[2]]
    h′ = Handle(h.conn, ety′, h.opt)
    hv = VectorWithHandle(h′, tv)
    BlockVector(offs, hv)
end

get_catalog(cat::PGCatalog) = cat

get_catalog(tbl::PGTable) = tbl.schema.catalog

function get_type(col::PGColumn)
    T = get_type(col.type)
    if !col.not_null
        T = Union{T, Missing}
    end
    T
end

function get_type(typ::PGType)
    key = Symbol(typ.name)
    if typ.schema.name == "pg_catalog"
        key = Symbol(typ.name)
        if haskey(PQ_SYSTEM_TYPES, key)
            key_oid = PQ_SYSTEM_TYPES[key]
            if haskey(_DEFAULT_TYPE_MAP, key_oid)
                return _DEFAULT_TYPE_MAP[key_oid]
            end
        end
    end
    String
end

guess_name(fk::PGForeignKey) =
    fk.name

guess_referring_name(fk::PGForeignKey) =
    "$(fk.table.name)_via_$(fk.name)"

function lookup(src::EntityShape{PGCatalog}, name::Symbol)
    p = lookup(EntityShape(src.ety["public"], src.opt, src.out), name)
    p !== nothing || return p
    p |> designate(src, target(p))
end

function lookup(src::EntityShape{PGSchema}, name::Symbol)
    scm = src.ety
    tbl = get(scm, string(name), nothing)
    tbl !== nothing || return nothing
    @assert tbl.primary_key !== nothing
    tbl_name = (tbl.schema.name, tbl.name)
    col_names = [col.name for col in tbl.primary_key.columns]
    col_types = Type[get_type(col) for col in tbl.primary_key.columns]
    p = load_postgres_table(tbl_name, col_names, col_types)
    tgt = EntityShape(tbl, src.opt, TupleOf(Symbol.(col_names), ValueOf.(col_types))) |> IsLabeled(name) |> BlockOf
    p |> designate(src, tgt)
end

function lookup(src::EntityShape{PGTable}, name::Symbol)
    tbl = src.ety
    for col in tbl
        if col.name == string(name)
            col_type = get_type(col)
            tbl_name = (tbl.schema.name, tbl.name)
            col_names = [col.name]
            col_types = Type[col_type]
            icol_names = [col.name for col in tbl.primary_key.columns]
            c = cover(ValueOf(col_type) |> IsLabeled(name))
            return chain_of(
                    load_postgres_table(tbl_name, col_names, col_types, icol_names),
                    block_cardinality(x1to1),
                    with_elements(chain_of(output(), column(1), c)),
                    flatten(),
            ) |> designate(src, target(c))
        end
    end
    for fk in tbl.foreign_keys
        if guess_name(fk) == string(name)
            ttbl = fk.target_table
            @assert ttbl.primary_key !== nothing
            tbl_name = (tbl.schema.name, tbl.name)
            col_names = [col.name for col in fk.columns]
            col_types = Type[get_type(col) for col in fk.columns]
            icol_names = [col.name for col in tbl.primary_key.columns]
            p0 = load_postgres_table(tbl_name, col_names, col_types, icol_names)
            tbl_name = (ttbl.schema.name, ttbl.name)
            col_names = [col.name for col in ttbl.primary_key.columns]
            col_types = Type[get_type(col) for col in ttbl.primary_key.columns]
            icol_names = [col.name for col in fk.target_columns]
            p1 = load_postgres_table(tbl_name, col_names, col_types, icol_names)
            card = x1to1
            if any(col -> !col.not_null, fk.columns)
                card |= x0to1
            end
            tgt = BlockOf(EntityShape(ttbl, src.opt, TupleOf(Symbol.(col_names), ValueOf.(col_types))) |> IsLabeled(name), card)
            return chain_of(
                    p0,
                    block_cardinality(x1to1),
                    with_elements(chain_of(p1, block_cardinality(card))),
                    flatten(),
            ) |> designate(src, tgt)
        end
    end
    for fk in tbl.referring_foreign_keys
        if guess_referring_name(fk) == string(name)
            ttbl = fk.table
            @assert ttbl.primary_key !== nothing
            tbl_name = (tbl.schema.name, tbl.name)
            col_names = [col.name for col in fk.target_columns]
            col_types = Type[get_type(col) for col in fk.target_columns]
            icol_names = [col.name for col in tbl.primary_key.columns]
            p0 = load_postgres_table(tbl_name, col_names, col_types, icol_names)
            tbl_name = (ttbl.schema.name, ttbl.name)
            col_names = [col.name for col in ttbl.primary_key.columns]
            col_types = Type[get_type(col) for col in ttbl.primary_key.columns]
            icol_names = [col.name for col in fk.columns]
            p1 = load_postgres_table(tbl_name, col_names, col_types, icol_names)
            card = x0to1
            if !any(uk -> uk.columns == fk.columns, ttbl.unique_keys)
                card |= x1toN
            end
            tgt = BlockOf(EntityShape(ttbl, src.opt, TupleOf(Symbol.(col_names), ValueOf.(col_types))) |> IsLabeled(name), card)
            return chain_of(
                    p0,
                    block_cardinality(x1to1),
                    with_elements(chain_of(p1, block_cardinality(card))),
                    flatten(),
            ) |> designate(src, tgt)
        end
    end
    nothing
end

lookup(src::EntityShape, name::Symbol) =
    error("not implemented")

end