text:
export SimpleSimulatedDAQ, SimpleSimulatedDAQParams, simulateLangevinInduced

Base.@kwdef struct SimpleSimulatedDAQParams <: DAQParams
  "All configured channels of this DAQ device."
  channels::Dict{String, DAQChannelParams}

  temperatureRise::Union{Dict{String, typeof(1.0u"K")}, Nothing} = nothing
  temperatureRiseSlope::Union{Dict{String, typeof(1.0u"s")}, Nothing} = nothing
  phaseChange::Union{Dict{String, typeof(1.0u"rad/K")}, Nothing} = nothing
  amplitudeChange::Union{Dict{String, typeof(1.0u"T/K")}, Nothing} = nothing
end

"Create the params struct from a dict. Typically called during scanner instantiation."
function SimpleSimulatedDAQParams(dict::Dict{String, Any})
  mainSectionFields = ["temperatureRise", "temperatureRiseSlope", "phaseChange", "amplitudeChange"]
  return createDAQParams(SimpleSimulatedDAQParams, dict)
end

Base.@kwdef mutable struct SimpleSimulatedDAQ <: AbstractDAQ
  "Unique device ID for this device as defined in the configuration."
  deviceID::String
  "Parameter struct for this devices read from the configuration."
  params::SimpleSimulatedDAQParams
  "Flag if the device is optional."
	optional::Bool = false
  "Flag if the device is present."
  present::Bool = false
  "Vector of dependencies for this device."
  dependencies::Dict{String, Union{Device, Missing}}

  # The following fields are only used for the simulation state!
  # The desired values for phase and amplitude are left within the sequence controller.
  "Base frequency to derive drive field frequencies."
  baseFrequency::typeof(1.0u"Hz") = 125.0u"MHz"
  "Divider of the baseFrequency to determine the drive field frequencies"
  divider::Union{Array{Int64, 2}, Missing} = missing
  "Applied drive field phase."
  phase::Union{Array{typeof(1.0u"rad"), 3}, Missing} = missing
  "Applied drive field voltage."
  amplitude::Union{Array{typeof(1.0u"T"), 3}, Array{typeof(1.0u"V"), 3}, Missing} = missing
  "Waveform type: sine, triangle or custom"
  waveform::Union{Array{Waveform, 2}, Missing} = missing

  txRunning::Bool = false

  "Number of periods within a frame."
  numPeriodsPerFrame::Int64 = 1
  "Number of samples per period."
  rxNumSamplingPoints::Int64 = 1

  rxChannelIDMapping::Dict{String, Int64} = Dict{String, Int64}()
  rxChanIDs::Vector{String} = []
  refChanIDs::Vector{String} = []

  currentFrame::Int64 = 1
  currentPeriod::Int64 = 1
end

function init(daq::SimpleSimulatedDAQ)
  @debug "Initializing simple simulated DAQ with ID `$(daq.deviceID)`."

  daq.present = true
end

neededDependencies(::SimpleSimulatedDAQ) = [SimulationController]
optionalDependencies(::SimpleSimulatedDAQ) = [TxDAQController, SurveillanceUnit]

Base.close(daq::SimpleSimulatedDAQ) = nothing

function setup(daq::SimpleSimulatedDAQ, sequence::Sequence)
  setupTx(daq, sequence)
  setupRx(daq, sequence)
end

function setupTx(daq::SimpleSimulatedDAQ, sequence::Sequence)
  daq.baseFrequency = txBaseFrequency(sequence)

  channels = electricalTxChannels(sequence)
  periodicChannels = [channel for channel in channels if channel isa PeriodicElectricalChannel]
  stepwiseChannels = [channel for channel in channels if channel isa StepwiseElectricalTxChannel]

  if !isempty(stepwiseChannels)
    @warn "The simple simulated DAQ can only process periodic channels. Other channels are ignored."
  end

  if any([length(component.amplitude) > 1 for channel in periodicChannels for component in channel.components])
    error("The simple simulated DAQ cannot work with more than one period in a frame or frequency sweeps yet.")
  end

  # Initialize fields
  numChannels_ = length(channels)
  numComponents_ = length([component for channel in periodicChannels for component in channel.components])
  daq.divider = zeros(Int64, (numComponents_, numChannels_))
  daq.phase = fill(0.0u"rad", (numComponents_, numChannels_, 1))
  daq.amplitude = fill(0.0u"T", (numComponents_, numChannels_, 1))
  daq.waveform = fill(WAVEFORM_SINE, (numComponents_, numChannels_))

  for channel in periodicChannels
    scannerChannel = daq.params.channels[channel.id]
    channelMapping = scannerChannel.channelIdx # Switch to getter?

    # Activate corresponding receive channels
    if !isnothing(scannerChannel.feedback)
      feedbackChannelID = scannerChannel.feedback.channelID
      scannerFeedbackChannel = daq.params.channels[feedbackChannelID]
      feedbackChannelIdx = scannerFeedbackChannel.channelIdx # Switch to getter?
      push!(daq.refChanIDs, feedbackChannelID)
      daq.rxChannelIDMapping[feedbackChannelID] = feedbackChannelIdx
    end

    componentIndex = 1
    for component in channel.components
      daq.divider[componentIndex, channelMapping, 1] = component.divider
      daq.phase[componentIndex, channelMapping, 1] = component.phase[1] # Only one period is allowed for now
      daq.waveform[componentIndex, channelMapping] = component.waveform # Only one period is allowed for now
      daq.amplitude[componentIndex, channelMapping, 1] = component.amplitude[1] # Only one period is allowed for now
      componentIndex += 1
    end
  end
end

function setupRx(daq::SimpleSimulatedDAQ, sequence::Sequence)
  daq.numPeriodsPerFrame = acqNumPeriodsPerFrame(sequence)
  daq.rxNumSamplingPoints = rxNumSamplesPerPeriod(sequence)

  for channel in rxChannels(sequence)
    scannerChannel = daq.params.channels[channel.id]
    push!(daq.rxChanIDs, channel.id)
    daq.rxChannelIDMapping[channel.id] = scannerChannel.channelIdx
  end
end

function startTx(daq::SimpleSimulatedDAQ)
  daq.txRunning = true
end

function stopTx(daq::SimpleSimulatedDAQ)
  daq.txRunning = false
end

function setTxParams(daq::SimpleSimulatedDAQ, Γ; postpone=false)
# Needs to update period and frame
end

function currentFrame(daq::SimpleSimulatedDAQ)
  return daq.currentFrame
end

function currentPeriod(daq::SimpleSimulatedDAQ)
  return daq.currentPeriod
end

function disconnect(daq::SimpleSimulatedDAQ)
end

enableSlowDAC(daq::SimpleSimulatedDAQ, enable::Bool, numFrames=0,
              ffRampUpTime=0.4, ffRampUpFraction=0.8) = 1

function readData(daq::SimpleSimulatedDAQ, startFrame::Integer, numFrames::Integer, numBlockAverages::Integer=1)
  startPeriod = (startFrame-1)*daq.numPeriodsPerFrame+1
  numPeriods = numFrames*daq.numPeriodsPerFrame

  uMeasPeriods, uRefPeriods, t = readDataPeriods(daq, startPeriod, numPeriods)

  measShape = (size(uMeasPeriods, 1), size(uMeasPeriods, 2), daq.numPeriodsPerFrame, numFrames)
  uMeas = reshape(uMeasPeriods, measShape)
  refShape = (size(uRefPeriods, 1), size(uRefPeriods, 2), daq.numPeriodsPerFrame, numFrames)
  uRef = reshape(uRefPeriods, refShape)
  tShape = (size(t, 1), daq.numPeriodsPerFrame, numFrames)
  t = reshape(t, tShape)

  return uMeas, uRef, t
end

function readDataPeriods(daq::SimpleSimulatedDAQ, startPeriod::Integer, numPeriods::Integer, numBlockAverages::Integer=1)
  simCont = dependencies(daq, SimulationController)[1]

  uMeas = zeros(daq.rxNumSamplingPoints, length(daq.rxChanIDs), numPeriods)u"V"
  uRef = zeros(daq.rxNumSamplingPoints, length(daq.refChanIDs), numPeriods)u"V"

  numSendChannels = size(daq.divider, 2)
  if numSendChannels != length(daq.rxChanIDs) || numSendChannels != length(daq.refChanIDs)
    error("For the simple simulated DAQ we assume a matching number of send, feedback and receive channels.")
  end

  cycle = lcm(daq.divider)/daq.baseFrequency
  startTime = upreferred(cycle*(startPeriod-1))
  t = collect(range(startTime, stop=startTime+numPeriods*cycle-1/daq.baseFrequency, length=daq.rxNumSamplingPoints*numPeriods))

  for (sendChannelID, sendChannel) in [(id, channel) for (id, channel) in daq.params.channels if channel isa DAQTxChannelParams]
    sendChannelIdx = sendChannel.channelIdx

    # Work with field for now
    if dimension(daq.amplitude[1]) == dimension(u"T")
      factor = 1.0u"T/T"
    else
      factor = 1/scannerChannel.calibration
    end

    temperatureRise = daq.params.temperatureRise[sendChannelID]
    temperatureRiseSlope = daq.params.temperatureRiseSlope[sendChannelID]
    phaseChange = daq.params.phaseChange[sendChannelID]
    amplitudeChange = daq.params.amplitudeChange[sendChannelID]

    # The temperature rises asymptotically to Tᵢₙᵢₜ+temperatureRise
    Tᵢₙᵢₜ = initialCoilTemperatures(simCont, sendChannelID)
    if t[1] == 0.0u"s"
      ΔT = fill(0.0u"K", size(t))
      ΔT[2:end] = temperatureRise*t[2:end]./(t[2:end].+temperatureRiseSlope)
    else
      ΔT = temperatureRise*t./(t.+temperatureRiseSlope)
    end
    T = Tᵢₙᵢₜ .+ ΔT
    ΔB = ΔT*amplitudeChange
    Δϕ = ΔT*phaseChange

    uₜₓ = zeros(length(t))u"V"
    uᵣₓ = zeros(length(t))u"V"
    uᵣₑ = zeros(length(t))u"V"
    for componentIdx in 1:size(daq.amplitude, 1)
      f = daq.baseFrequency/daq.divider[componentIdx, sendChannelIdx]
      Bₘₐₓ = daq.amplitude[componentIdx, sendChannelIdx, 1]*factor
      ϕ = daq.phase[componentIdx, sendChannelIdx, 1]

      Bᵢ = Bₘₐₓ.*sin.(2π*f*t.+ϕ) # Desired, ideal field without drift
      Bᵣ = (Bₘₐₓ.+ΔB).*sin.(2π*f*t.+ϕ.+Δϕ) # Field with drift of phase and amplitude

      uₜₓ .+= Bᵢ.*sendChannel.calibration
      uᵣₓ .+= simulateLangevinInduced(t, Bᵣ, f, ϕ.+Δϕ) # f is not completely correct due to the phase change, but this is only a rough approximation anyways

      # Assumes the same induced voltage from the field as given out with uₜₓ,
      # just with a slight change in phase and amplitude
      uᵣₑ .+= Bᵣ.*sendChannel.calibration
    end

    # Assumes one reference and one measurement channel for each send channel
    uMeas[:, sendChannelIdx, :] = reshape(uᵣₓ, (daq.rxNumSamplingPoints, 1, numPeriods))
    uRef[:, sendChannelIdx, :] = reshape(uᵣₑ, (daq.rxNumSamplingPoints, 1, numPeriods))

    # Update temperature in simulation controller for further use in e.g. simulated temperature sensors
    currentCoilTemperatures(simCont, sendChannelID, T[end])
  end

  totalPeriods = daq.currentPeriod+numPeriods
  daq.currentPeriod = totalPeriods % daq.numPeriodsPerFrame
  daq.currentFrame = totalPeriods ÷ daq.numPeriodsPerFrame

  return uMeas, uRef, reshape(t, (daq.rxNumSamplingPoints, numPeriods))
end

numTxChannelsTotal(daq::SimpleSimulatedDAQ) = 10 # Arbitrary number, since we are just simulating
numRxChannelsTotal(daq::SimpleSimulatedDAQ) = 10 # Arbitrary number, since we are just simulating
numTxChannelsActive(daq::SimpleSimulatedDAQ) = length([channel for (id, channel) in daq.params.channels if channel isa TxChannelParams])
numRxChannelsActive(daq::SimpleSimulatedDAQ) = numRxChannelsReference(daq)+numRxChannelsMeasurement(daq)
numRxChannelsReference(daq::SimpleSimulatedDAQ) = length(daq.refChanIDs)
numRxChannelsMeasurement(daq::SimpleSimulatedDAQ) = length(daq.rxChanIDs)
numComponentsMax(daq::SimpleSimulatedDAQ) = 1
canPostpone(daq::SimpleSimulatedDAQ) = false
canConvolute(daq::SimpleSimulatedDAQ) = false

"Very, very basic simulation of an MPI signal using the Langevin function."
function simulateLangevinInduced(t::Vector{typeof(1.0u"s")}, B::Vector{typeof(1.0u"T")}, f::typeof(1.0u"Hz"), ϕ::Union{typeof(1.0u"rad"), Vector{typeof(1.0u"rad")}})
  Bₘₐₓ = maximum(B)

  c = 0.5e-6 # mol/m^3
  ν = 22.459e3 # mol/m^3
  Mₛ = 477e3 # A/m (for magnetite)
  Dₖ = 30/1e9 # m
  μ₀ = 4π*1e-7 # Vs/Am

  k_B = 1.38064852*1e-23 # J/K (Boltzmann constant)
  Tₐ = 273.15+20 # K

  ξ = (π*Dₖ^3*Mₛ*ustrip.(u"T", B))/(3*k_B*Tₐ)
  ξ̇ = (2π^2*Dₖ^3*Mₛ*ustrip.(u"Hz", f)*ustrip.(u"T", Bₘₐₓ))/(3*k_B*Tₐ)*cos.(2π*ustrip.(u"Hz", f)*ustrip.(u"s", t).+ustrip.(u"rad", ϕ))
  Ṁ = c/(3*ν)*Mₛ*(-ξ̇.*csch.(ξ).^2 + ξ̇./(ξ.^2))

  u = -Ṁ; # Yes, this is wrong since it doesn't reflect the coil geometry or anything else, but I don't care at this point

  scalingFactor = 1/maximum(x->isnan(x) ? -Inf : x, u)
  u = u.*scalingFactor # Scale to ±1 for the following correction and leave it like that; we assume a nice LNA ;)

  # This at least makes the function work for most values of ϕ.
  # I don't think we have to care about all special cases here.
  if ϕ == 0.0
    indices = findall(x->isnan(x), u)
    for index in indices
      if index < length(u)
        u[index] = sign(u[index+1])
      else
        u[index] = sign(u[index-1])
      end
    end

    # Detect jumps, and yes this is crude...
    for i in 1:length(u)-1
      du = u[i]-u[i+1]
      if du > 0.2
        u[i+1] = sign(u[i+2])
      end
    end
  end

  return u*u"V"
end


text:
module PlayGround

# Version 0.1.0
include("styled-messages.jl")
include("Echidna.jl")
include("Newton-Raphson.jl")
include("Simulation.jl")

end # module

ext: jl


text:
"""
    BasisFunction{B <: AbstractBSplineBasis, T}

Describes a single basis function.

The basis function may belong to a [`BSplineBasis`](@ref) (in which case it's
effectively a B-spline), or to a basis derived from a B-spline basis (such as a
`RecombinedBSplineBasis`).

---

    BasisFunction(basis::AbstractBSplineBasis, i::Int, [T = Float64])

Construct i-th basis function of the given basis.

The constructed function can be evaluated as `b(x)`, returning a value of type
`T`.

---

    (b::BasisFunction)(x, [op::AbstractDifferentialOp])

Evaluate basis function at coordinate `x`.

To evaluate a derivative, pass `Derivative(n)` as the `op` argument, with `n`
the derivative order.

To evaluate multiple derivatives, pass a derivative range `Derivative(m:n)`.
In particular, `Derivative(m:n)` evaluates the basis function itself and its
first `n` derivatives.

More general differential operators, such as `Derivative(n) + λ Derivative(m)`,
are also supported.
"""
struct BasisFunction{Basis <: AbstractBSplineBasis, T}
    basis :: Basis
    i     :: Int
    @inline function BasisFunction(
            b::AbstractBSplineBasis, i, ::Type{T} = Float64) where {T}
        Basis = typeof(b)
        new{Basis, T}(b, i)
    end
end

basis(b::BasisFunction) = b.basis
knots(b::BasisFunction) = knots(basis(b))
order(b::BasisFunction) = order(basis(b))
Base.eltype(::Type{BasisFunction{B,T}}) where {B,T} = T

function Base.show(io::IO, b::BasisFunction)
    print(io, "Basis function i = ", b.i, "\n")
    print(io, "  from ")
    summary(io, basis(b))
    ind = support(b)
    i, j = first(ind), last(ind)
    ts = knots(b)
    ti, tj = map(n -> ts[n], (i, j))
    print(io, "\n  support: [", ti, ", ", tj, ") (knots ", ind, ")")
end

"""
    support(b::BasisFunction) -> UnitRange{Int}

Get range of knots supported by the basis function.

Returns the knot range `i:j` such that the basis function support is
``t ∈ [t_i, t_j)``.
"""
support(b::BasisFunction) = support(basis(b), b.i)

"""
    support(B::AbstractBSplineBasis, i::Integer) -> UnitRange{Int}

Get range of knots supported by the ``i``-th basis function.
"""
support(B::BSplineBasis, i::Integer) = i:(i + order(B))

"""
    nonzero_in_segment(B::AbstractBSplineBasis, n::Int) -> UnitRange{Int}

Returns the range of basis functions that are non-zero in a given knot segment.

The ``n``-th knot segment is defined by ``Ω_n = [t_n, t_{n + 1}]``.

For [`BSplineBasis`](@ref) and `RecombinedBSplineBasis`, the number of non-zero
functions in any given segment is generally equal to the B-spline order ``k``.
This number decreases near the borders, but this is not significant when
B-spline knots have multiplicity ``k`` there (as is the default).

For B-spline bases, and excepting the borders, the non-zero B-splines are
``\\left\\{ b_i \\right\\}_{i = n - k + 1}^{n}``.
This function thus returns `(n - k + 1):N` when `B` is a `BSplineBasis`.

See also [`support`](@ref) for the inverse operation.
"""
function nonzero_in_segment(B::BSplineBasis, n::Int; N = length(B))
    a = clamp(n - order(B) + 1, 1, typemax(Int))
    b = clamp(n, typemin(Int), N)
    a:b
end

"""
    common_support(b1::BasisFunction, b2::BasisFunction, ...) -> UnitRange{Int}

Get range of knots commonly supported by different basis functions.

If the supports don't intersect, an empty range is returned (e.g. `6:5`),
following the behaviour of `intersect`. The lack of intersection can be checked
using `isempty`, which returns `true` for such a range.
"""
common_support(bs::Vararg{BasisFunction}) = ∩(support.(bs)...)

(b::BasisFunction)(x, op=Derivative(0)) = evaluate(basis(b), b.i, x, op, eltype(b))

"""
    evaluate(B::AbstractBSplineBasis, i::Integer, x,
             [op::AbstractDifferentialOp], [T=Float64])

Evaluate ``i``-th basis function in the given basis at `x` (can be a coordinate or a
vector of coordinates).

To evaluate a derivative, pass `Derivative(n)` as the `op` argument, with `n`
the derivative order.

More general differential operators, such as `Derivative(n) + λ Derivative(m)`,
are also supported.

See also [`evaluate!`](@ref).
"""
function evaluate(B::BSplineBasis, i::Integer, x::Real,
                  op::AbstractDifferentialOp = Derivative(0),
                  ::Type{T} = Float64) where {T}
    N = length(B)
    if !(1 <= i <= N)
        throw(DomainError(i, "Basis function index must be in 1:$N"))
    end
    k = order(B)
    t = knots(B)
    evaluate_diff(op, BSplineOrder(k), t, i, x, T)
end

evaluate(B::BSplineBasis, i::Integer, x::Real, ::Type{T}) where {T} =
    evaluate(B, i, x, Derivative(0), T)

# No derivative
evaluate_diff(::Derivative{0}, ::BSplineOrder{k}, t, i, x,
              ::Type{T}) where {k,T} = _evaluate(BSplineOrder(k), t, i, x, T)

# N-th derivative
function evaluate_diff(::Derivative{N}, ::BSplineOrder{k}, t, i, x,
                       ::Type{T}) where {N,k,T}
    @assert N > 0
    y = zero(T)
    dt = t[i + k - 1] - t[i]
    local dy::T  # make sure `y` doesn't change type, e.g. if T = Float32
    if !iszero(dt)
        # Recursively evaluate derivative `N - 1` of B-spline of order `k - 1`.
        dy = evaluate_diff(Derivative(N - 1), BSplineOrder(k - 1), t, i, x, T) / dt
        y += dy
    end
    dt = t[i + k] - t[i + 1]
    if !iszero(dt)
        dy = evaluate_diff(Derivative(N - 1), BSplineOrder(k - 1), t, i + 1, x, T) / dt
        y -= dy
    end
    (y * (k - 1)) :: T
end

# TODO this can be optimised... it should be possible to compute all derivatives
# at once.
@generated function evaluate_diff(
        ::DerivativeUnitRange{m, n}, args...,
    ) where {m, n}
    ex = :(())
    for i = m:n
        push!(ex.args, :(evaluate_diff(Derivative($i), args...)))
    end
    ex
end

# More general diff operators
evaluate_diff(S::ScaledDerivative, etc...) =
    S.α * evaluate_diff(S.D, etc...)

function evaluate_diff(S::DifferentialOpSum, etc...)
    vals = map(D -> evaluate_diff(D, etc...), (S.a, S.b))
    sum(vals)
end

# For evaluation at multiple locations (allocates a vector).
evaluate(B::BSplineBasis, i, x::AbstractVector, args...) =
    evaluate.(B, i, x, args...)

"""
    evaluate!(b::AbstractVector, B::BSplineBasis, i::Integer,
              x::AbstractVector, args...)

Evaluate i-th basis function at positions `x` and write result to `b`.

See also [`evaluate`](@ref).
"""
function evaluate!(b::AbstractVector{T}, B::AbstractBSplineBasis, i,
                   x::AbstractVector, args...) where {T}
    broadcast!(x -> evaluate(B, i, x, args..., T), b, x)
end

# Specialisation for first order B-splines.
function _evaluate(::BSplineOrder{1}, t::AbstractVector, i::Integer,
                   x::Real, ::Type{T}) where {T}
    # Local support of the B-spline.
    @inbounds ta = t[i]
    @inbounds tb = t[i + 1]

    in_local_support(x, t, ta, tb) ? one(T) : zero(T)
end

# General case of order k >= 2.
function _evaluate(::BSplineOrder{k}, t::AbstractVector, i::Integer,
                   x::Real, ::Type{T}) where {T,k}
    k::Int
    @assert k >= 2

    # Local support of the B-spline.
    @inbounds ta = t[i]
    @inbounds tb = t[i + k]

    if !in_local_support(x, t, ta, tb)
        return zero(T)
    end

    @inbounds ta1 = t[i + 1]
    @inbounds tb1 = t[i + k - 1]

    # Recursively evaluate lower-order B-splines.
    y = zero(T)
    local dy::T  # make sure `y` doesn't change type, e.g. if T = Float32

    if tb1 != ta
        dy = _evaluate(BSplineOrder(k - 1), t, i, x, T) *
            (x - ta) / (tb1 - ta)
        y += dy
    end

    if ta1 != tb
        dy = _evaluate(BSplineOrder(k - 1), t, i + 1, x, T) *
            (tb - x) / (tb - ta1)
        y += dy
    end

    y :: T
end

# In principle, the support of a B-spline is [ta, tb[.
# The exception is when tb is the last point in the domain, in which case we
# include it in the interval (otherwise the right boundary always has zero value).
@inline in_local_support(x, t, ta, tb) = (ta <= x < tb) || (x == tb == t[end])

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
module UnitCellPlotter

using PyCall

"""
plot_udc_ruc()

# Args

- ruc_bbox::NTuple{4, Float64},
- inclusions::Dict{String, Matrix{Float64}};
- image_path::Union{String, Nothing} = nothing,
- matrix_color::Union{String, Symbol} = :grey,
- fibre_color::Union{String, Symbol} = :black,
- matrix_edge_color::Union{String, Symbol, Nothing} = nothing,
- fibre_edge_color::Union{String, Symbol, Nothing} = nothing,
- z_comp_in_data:: Bool = true,

"""
function plot_unit_cell(
    ruc_bbox::NTuple{4,Float64},
    inclusions::Dict{String,Matrix{Float64}};
    image_path::Union{String,Nothing} = nothing,
    matrix_color::Union{String,Symbol} = :grey,
    fibre_color::Union{String,Symbol} = :black,
    matrix_edge_color::Union{String,Symbol,Nothing} = nothing,
    fibre_edge_color::Union{String,Symbol,Nothing} = nothing,
    fibre_edge_thickness::Union{Float64,Nothing} = nothing,
    z_comp_in_data::Bool = true
)
    @pyinclude(joinpath(@__DIR__, "plot_2D_shapes.py"))
    py"plot_unit_cell"(
        ruc_bbox, inclusions,
        image_path = image_path,
        matrix_color = matrix_color,
        fibre_color = fibre_color,
        matrix_edge_color = matrix_edge_color,
        fibre_edge_color = fibre_edge_color,
        fibre_edge_thickness = fibre_edge_thickness,
        z_comp_in_data = z_comp_in_data,
    )

end

export plot_unit_cell

end


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

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

export BrickletLine
"""
Measures reflectivity of a surface
"""
mutable struct BrickletLine <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

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

        return new(ipcon, deviceInternal)
    end
end

export get_reflectivity
"""
    $(SIGNATURES)

Returns the currently measured reflectivity. The reflectivity is
a value between 0 (not reflective) and 4095 (very reflective).

Usually black has a low reflectivity while white has a high
reflectivity.

If you want to get the reflectivity periodically, it is recommended
to use the :cb:`Reflectivity` callback and set the period with
:func:`Set Reflectivity Callback Period`.
"""
function get_reflectivity(device::BrickletLine)
    return device.deviceInternal.get_reflectivity()
end

export set_reflectivity_callback_period
"""
    $(SIGNATURES)

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

The :cb:`Reflectivity` callback is only triggered if the reflectivity has
changed since the last triggering.
"""
function set_reflectivity_callback_period(device::BrickletLine, period)
    device.deviceInternal.set_reflectivity_callback_period(period)
end

export get_reflectivity_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Reflectivity Callback Period`.
"""
function get_reflectivity_callback_period(device::BrickletLine)
    return device.deviceInternal.get_reflectivity_callback_period()
end

export set_reflectivity_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Reflectivity 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 reflectivity is *outside* the min and max values"
 "'i'",    "Callback is triggered when the reflectivity is *inside* the min and max values"
 "'<'",    "Callback is triggered when the reflectivity is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the reflectivity is greater than the min value (max is ignored)"
"""
function set_reflectivity_callback_threshold(device::BrickletLine, option, min, max)
    device.deviceInternal.set_reflectivity_callback_threshold(option, min, max)
end

export get_reflectivity_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Reflectivity Callback Threshold`.
"""
function get_reflectivity_callback_threshold(device::BrickletLine)
    return device.deviceInternal.get_reflectivity_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callback

* :cb:`Reflectivity Reached`

is triggered, if the threshold

* :func:`Set Reflectivity Callback Threshold`

keeps being reached.
"""
function set_debounce_period(device::BrickletLine, 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::BrickletLine)
    return device.deviceInternal.get_debounce_period()
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::BrickletLine)
    return device.deviceInternal.get_identity()
end

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


text:
 simple_profile()=
        return Dict(
            "username"=> user_name(),
            "name"=> first_name(),
            "sex"=> rand(["M","F"]),
            "address"=> address(),
            "mail"=> free_email(),
            #"password"=>password(),
            "birthdate"=> date() )

     profile( fields="None")=
(
        fields == "None" && (fields = []);
        d = Dict(
            "job"=> job(),
            "company"=> company(),
            #"ssn"=> ssn(),
            "residence"=> address(),
            "current_location"=> (latitude(), longitude()),

            "blood_group"=> rand(["A", "B", "AB", "O"])*rand(["+", "-"]),

            "website"=> [url() for i=1:rand(2:5)]
        );
        if length(fields) >= 1
                e = Dict()
                for (k, v) in d
                    if (k in fields)
                        push!(e, k => v)
                    end
                end
                return e
        end;

        return d
)


text:
"""
    SemidirectProductOperation(action::AbstractGroupAction)

Group operation of a semidirect product group. The operation consists of the operation
`opN` on a normal subgroup `N`, the operation `opH` on a subgroup `H`, and an automorphism
`action` of elements of `H` on `N`. Only the action is stored.
"""
struct SemidirectProductOperation{A<:AbstractGroupAction} <: AbstractGroupOperation
    action::A
end

function show(io::IO, op::SemidirectProductOperation)
    print(io, "SemidirectProductOperation($(op.action))")
end

const SemidirectProductGroup{N,H,A} =
    GroupManifold{ProductManifold{Tuple{N,H}},SemidirectProductOperation{A}}

@doc raw"""
    SemidirectProductGroup(N::GroupManifold, H::GroupManifold, A::AbstractGroupAction)

A group that is the semidirect product of a normal group $\mathcal{N}$ and a subgroup
$\mathcal{H}$, written $\mathcal{G} = \mathcal{N} ⋊_θ \mathcal{H}$, where
$θ: \mathcal{H} × \mathcal{N} → \mathcal{N}$ is an automorphism action of $\mathcal{H}$ on
$\mathcal{N}$. The group $\mathcal{G}$ has the composition rule

````math
g \circ g' = (n, h) \circ (n', h') = (n \circ θ_h(n'), h \circ h')
````

and the inverse

````math
g^{-1} = (n, h)^{-1} = (θ_{h^{-1}}(n^{-1}), h^{-1}).
````
"""
function SemidirectProductGroup(N::GroupManifold, H::GroupManifold, A::AbstractGroupAction)
    N === g_manifold(A) || error("Subgroup $(N) must be the G-manifold of action $(A)")
    H === base_group(A) || error("Subgroup $(H) must be the base group of action $(A)")
    op = SemidirectProductOperation(A)
    M = ProductManifold(N, H)
    return GroupManifold(M, op)
end

function show(io::IO, G::SemidirectProductGroup)
    M = base_manifold(G)
    N, H = M.manifolds
    A = G.op.action
    print(io, "SemidirectProductGroup($(N), $(H), $(A))")
end

submanifold(G::SemidirectProductGroup, i) = submanifold(base_manifold(G), i)

_padpoint!(G::SemidirectProductGroup, q) = q

_padvector!(G::SemidirectProductGroup, X) = X

inv(G::GT, e::Identity{GT}) where {GT<:SemidirectProductGroup} = e

function inv!(G::SemidirectProductGroup, q, p)
    M = base_manifold(G)
    N, H = M.manifolds
    A = G.op.action
    np, hp = submanifold_components(G, p)
    nq, hq = submanifold_components(G, q)
    inv!(H, hq, hp)
    npinv = inv(N, np)
    apply!(A, nq, hq, npinv)
    @inbounds _padpoint!(G, q)
    return q
end
inv!(G::AG, p, e::Identity{AG}) where {AG<:SemidirectProductGroup} = identity!(G, p, e)

identity(G::GT, e::Identity{GT}) where {GT<:SemidirectProductGroup} = e

function identity!(G::SemidirectProductGroup, q, p)
    M = base_manifold(G)
    N, H = M.manifolds
    np, hp = submanifold_components(G, p)
    nq, hq = submanifold_components(G, q)
    identity!(N, nq, np)
    identity!(H, hq, hp)
    @inbounds _padpoint!(G, q)
    return q
end
identity!(G::GT, e::E, ::E) where {GT<:SemidirectProductGroup,E<:Identity{GT}} = e

compose(G::GT, p, e::Identity{GT}) where {GT<:SemidirectProductGroup} = p
compose(G::GT, e::Identity{GT}, p) where {GT<:SemidirectProductGroup} = p
compose(G::GT, e::E, ::E) where {GT<:SemidirectProductGroup,E<:Identity{GT}} = e

function compose!(G::SemidirectProductGroup, x, p, q)
    M = base_manifold(G)
    N, H = M.manifolds
    A = G.op.action
    np, hp = submanifold_components(G, p)
    nq, hq = submanifold_components(G, q)
    nx, hx = submanifold_components(G, x)
    compose!(H, hx, hp, hq)
    nxtmp = apply(A, hp, nq)
    compose!(N, nx, np, nxtmp)
    @inbounds _padpoint!(G, x)
    return x
end
compose!(G::GT, x, ::Identity{GT}, q) where {GT<:SemidirectProductGroup} = copyto!(x, q)
compose!(G::GT, x, p, ::Identity{GT}) where {GT<:SemidirectProductGroup} = copyto!(x, p)
function compose!(G::GT, x, e::E, ::E) where {GT<:SemidirectProductGroup,E<:Identity{GT}}
    return identity!(G, x, e)
end

function translate_diff!(G::SemidirectProductGroup, Y, p, q, X, conX::LeftAction)
    M = base_manifold(G)
    N, H = M.manifolds
    A = G.op.action
    np, hp = submanifold_components(G, p)
    nq, hq = submanifold_components(G, q)
    nX, hX = submanifold_components(G, X)
    nY, hY = submanifold_components(G, Y)
    translate_diff!(H, hY, hp, hq, hX, conX)
    nZ = apply_diff(A, hp, nq, nX)
    nr = apply(A, hp, nq)
    translate_diff!(N, nY, np, nr, nZ, conX)
    @inbounds _padvector!(G, Y)
    return Y
end

function get_vector!(G::SemidirectProductGroup, Y, p, X, B::VeeOrthogonalBasis)
    M = base_manifold(G)
    N, H = M.manifolds
    dimN = manifold_dimension(N)
    dimH = manifold_dimension(H)
    @assert length(X) == dimN + dimH
    np, hp = submanifold_components(G, p)
    nY, hY = submanifold_components(G, Y)
    get_vector!(N, nY, np, view(X, 1:dimN), B)
    get_vector!(H, hY, hp, view(X, dimN+1:dimN+dimH), B)
    @inbounds _padvector!(G, Y)
    return Y
end
eval(quote
    @invoke_maker 1 Manifold get_vector!(
        M::SemidirectProductGroup,
        Xⁱ,
        e::Identity,
        X,
        B::VeeOrthogonalBasis,
    )
end)

function get_coordinates!(G::SemidirectProductGroup, Y, p, X, B::VeeOrthogonalBasis)
    M = base_manifold(G)
    N, H = M.manifolds
    dimN = manifold_dimension(N)
    dimH = manifold_dimension(H)
    @assert length(Y) == dimN + dimH
    np, hp = submanifold_components(G, p)
    nY, hY = submanifold_components(G, X)
    get_coordinates!(N, view(Y, 1:dimN), np, nY, B)
    get_coordinates!(H, view(Y, dimN+1:dimN+dimH), hp, hY, B)
    return Y
end
eval(quote
    @invoke_maker 1 Manifold get_coordinates!(
        M::SemidirectProductGroup,
        Y,
        e::Identity,
        X,
        B::VeeOrthogonalBasis,
    )
end)

function zero_tangent_vector(G::SemidirectProductGroup, p)
    X = allocate_result(G, zero_tangent_vector, p)
    zero_tangent_vector!(G, X, p)
    return X
end

function zero_tangent_vector!(G::SemidirectProductGroup, X, p)
    M = base_manifold(G)
    N, H = M.manifolds
    np, hp = submanifold_components(G, p)
    nX, hX = submanifold_components(G, X)
    zero_tangent_vector!(N, nX, np)
    zero_tangent_vector!(H, hX, hp)
    return X
end

function isapprox(G::SemidirectProductGroup, p, q; kwargs...)
    M = base_manifold(G)
    N, H = M.manifolds
    np, hp = submanifold_components(G, p)
    nq, hq = submanifold_components(G, q)
    return isapprox(N, np, nq; kwargs...) && isapprox(H, hp, hq; kwargs...)
end
function isapprox(G::SemidirectProductGroup, p, X, Y; kwargs...)
    M = base_manifold(G)
    N, H = M.manifolds
    np, hp = submanifold_components(G, p)
    nX, hX = submanifold_components(G, X)
    nY, hY = submanifold_components(G, Y)
    return isapprox(N, np, nX, nY; kwargs...) && isapprox(H, hp, hX, hY; kwargs...)
end
function isapprox(G::GT, p, e::Identity{GT}; kwargs...) where {GT<:SemidirectProductGroup}
    return isapprox(G, e, p; kwargs...)
end
function isapprox(G::GT, e::Identity{GT}, p; kwargs...) where {GT<:SemidirectProductGroup}
    return isapprox(G, identity(G, p), p; kwargs...)
end
function isapprox(
    ::GT,
    ::E,
    ::E;
    kwargs...,
) where {GT<:SemidirectProductGroup,E<:Identity{GT}}
    return true
end


text:
struct Model
    vars::OffsetVector{Variable,Vector{Variable}}
    funcs::Vector{Func}
    loops::Vector{LoopSet}
    tracker::Dict{Symbol,Int} # maps symbol to varid
    mod::Module
    modsym::Symbol
    funcdict::Dict{Func,Int}
    # inputvars::Vector{Int}
    varid::Ref{UInt}
    tracked::Vector{Bool}
    gradmodel::Bool
    targets::Vector{Int}
    # targetinit::RefValue{Bool}
    function Model(mod::Module = ReverseDiffExpressions, gradmodel::Bool = false)
        target = Variable(Symbol("##TARGET##"), 0); target.tracked = true
        ein = Variable(Symbol("##ONE##"), 1)
        new(OffsetVector(Variable[target, ein], -1), Func[], LoopSet[], Dict{Symbol,Int}(), mod, Symbol(mod), Dict{Func,Int}(), Ref(UInt(2)), Bool[], gradmodel, [0])
    end
end
# struct DiffRuleVarFunc
    # diffrule::DiffRule
    # func::Func
    # variables::OffsetVector{Variable,Vector{Variable}}
# end

struct ∂Model # Houses meta data.
    m::Model
    mold::Model
    mapping::Vector{Int} # maps old vars to new ones
    ∂mapping::Vector{Int} # maps old vars to new ∂ones
    gradvar::Variable
end
function ∂Model(mold::Model)
    m = Model(mold.mod, true)
    gradvar = addvar!(m, Symbol("#∇#"))
    ∂Model(m, mold, Int[], Int[], gradvar)
end

# corresponding_func(∂m::∂Model, func::Func) = ∂m.m.funcs[∂m.mapping[∂m.mold.funcdict[func]]]

reset_funclowered!(m::Model) = foreach(f -> (f.lowered[] = false), m.funcs)
function reset_varinitialized!(m::Model)
    foreach(m.vars) do v
        v.lowered_count = 0
    end
end
reset_var_tracked!(m::Model) = foreach(v -> (v.tracked = false), m.vars)
onevar(m::Model) = @inbounds m.vars[1]
targetvar(m::Model) = @inbounds m.vars[0]

getparent(m::Model, v::Variable, i::Int = 1) = getfunc(m, parents(v)[i])
getfunc(m::Model, i::Integer) = m.funcs[i]
getvar(m::Model, i::Integer) = m.vars[i]
initialize_var!(m::Model, vid::Integer) = getvar(m, vid).initialized = true

function addvar!(m::Model, s::Symbol)::Variable
    addvar!(m, Variable(s, length(m.vars)))
end
function addvar!(m::Model, v::Variable)::Variable
    @unpack vars, tracker = m
    push!(vars, v)
    tracker[v.name] = v.varid
    v
end
function addvarref!(m::Model, ref)::Variable
    @assert !(ref isa Variable)
    @unpack vars, tracker = m
    s = gensym()
    v = Variable(s, length(m.vars))
    push!(vars, v)
    tracker[s] = v.varid
    v.ref = ref
    v
end
getvar!(m::Model, ref)::Variable = addvarref!(m, ref)

function getvar!(m::Model, s::Symbol)::Variable
    @unpack vars, tracker = m
    id = get(tracker, s, nothing)
    isnothing(id) ? addvar!(m, s) : vars[id]
end
getvar!(m::Model, s::Expr)::Variable = addvar!(m, s)

function getvar!(∂m::∂Model, vold::Variable)::Variable
    @unpack m = ∂m
    @unpack vars, tracker = m
    s = name(vold)
    vnew = if isref(vold)
        v = addvar!(m, s)
        v.ref = vold.ref
        v
    else
        id = get(tracker, s, nothing)
        isnothing(id) ? addvar!(m, s) : vars[id]
    end
    vnew.tracked = vold.tracked
    vnew
end

function addvar!(m::Model, s::Expr)::Variable # unpack LHS tuple
    @assert s.head === :tuple
    packedtuple = gensym(:packedtuple)
    pt = addvar!(m, packedtuple)
    for i ∈ eachindex(s.args)
        gi = Func(:getindex)
        v = getvar!(m, s.args[i])
        returns!(gi, v)
        uses!(gi, pt)
        uses!(gi, getvar!(m, i))
        addfunc!(m, gi)
    end
    pt
end

function Func(m::Model, instr::Symbol, args...)
    ins = Instruction(instr)
    ins = ins ∈ keys(LoopVectorization.COST) ? ins : Instruction(m.modsym, instr)
    Func(ins, args...)
end

function addfunc!(m::Model, f::Func)
    fid = get!(m.funcdict, f) do
        push!(m.funcs, f)
        fid = length(m.funcs)
        for p ∈ f.vparents
            push!(m.vars[p].useids, fid)
        end
        fid
    end
    ret = m.vars[m.funcs[fid].output[]]
    push!(parents(ret), fid)
    @assert !any(in(f.output[]), f.vparents)
    ret
end

uses!(func::Func, m::Model, x::Symbol) = uses!(func, getvar!(m, x))
# function uses!(f::Func, m::Model, x)
#     @unpack otherargs, vparents = f
#     push!(otherargs, x)
#     push!(vparents, -length(otherargs))
#     nothing
# end

function var_tracked_search!(m::Model, v::Variable)
    v.tracked && return true
    for fid ∈ parents(v), vid ∈ parents(m.funcs[fid])
        vp = m.vars[vid]
        var_tracked_search!(m, vp) && return (v.tracked = vp.tracked = true)
    end
    false
end

propagate_var_tracked!(m::Model) = foreach(v -> var_tracked_search!(m, v), m.vars)

function getconstindex!(m::Model, vout::Variable, vin::Variable, index::Int)
    if 0 < index < 10
        f = [:first, :second, :third, :fourth, :fifth, :sixth, :seventh, :eighth, :ninth][index]
        func = Func(Instruction(:ReverseDiffExpressionsBase, f), false)
        uses!(func, vin)
    else
        func = Func(Instruction(:Base, :getindex), false)
        vindex = addvar!(m, Symbol(""))
        vindex.ref = index
        uses!(func, vin)
        uses!(func, vindex)
    end
    returns!(func, vout)
    addfunc!(m, func)
    func
end

function nonindex_parent(m::Model, func::Func)
    if isindexfunc(func)
        parent = getvar(m, first(parents(func)))
        nonindex_parent(m, getfunc(m, first(parents(parent))))
    else
        func
    end
end


text:
# https://github.com/joshday/OnlineStats.jl/blob/master/notebooks/animations.jl
#using UnicodePlots
using Markdown
using InteractiveUtils
using JLD
using Pkg
using PyCall
using PyPlot

trace = load("JLD/PopulationScatter.jld", "trace")
evo_population = [t.metadata["pop"] for t in trace]
evo_loss = [t.value for t in trace]
evo_iteration = [t.iteration for t in trace]
points = last(evo_population)

using PyCall
@pyimport matplotlib.animation as anim
using PyPlot

function showanim(filename)
    base64_video = base64encode(open(filename))
    display("text/html", """<video controls src="data:video/x-m4v;base64,$base64_video">""")
end

fig = figure(figsize=(4,4))
ax = axes()

anim = @animate for (points, j, k) in zip(evo_population, evo_loss, evo_iteration)
    X = points[1,:]
    Y = points[2,:]
    Z = points[3,:]

    #scatter(X, Y, Z)
    show(PyPlot.scatter(X, Y, c=[j for i in 1:4]))#,markersize=j))#,markercolor=j))
    show(PyPlot.scatter(Z, Y, c=[j for i in 1:4]))#,markersize=j))#,markercolor=j))
    show(PyPlot.scatter(X, Z, c=[j for i in 1:4]))#,markersize=j))#,markercolor=j))

    #scatter(x,y, marker_z = z, markersize = 5*z,  color = :jet)
    #using VegaLite
    #display(@vlplot(:point, x=:X, y=:Y, color=:j))
    #cnt+=1
end
gif(anim, joinpath(directed, "evo.gif"), fps = 10)


text:
#=

Solves the price function for the Lucas tree in a continuous state
setting, using piecewise linear approximation for the sequence of
candidate price functions.  The consumption endownment follows the
log linear AR(1) process

    log y' = alpha log y + sigma epsilon

where y' is a next period y and epsilon is an iid standard normal
shock. Hence

    y' = y^alpha * xi   where xi = e^(sigma * epsilon)

The distribution phi of xi is

    phi = LN(0, sigma^2) where LN means lognormal

@authors : Spencer Lyon <spencer.lyon@nyu.edu>, John Stachurski

References
----------

http://quant-econ.net/jl/markov_asset.html

=#

using QuantEcon
using Distributions

"""
A function that takes two arrays and returns a function that approximates the
data using continuous piecewise linear interpolation.

"""
function lin_interp(x_vals::Vector{Float64}, y_vals::Vector{Float64})
    # == linear interpolation inside grid, constant values outside grid == #
    w = LinInterp(x_vals, y_vals)
    return w
end

"""
The Lucas asset pricing model --- parameters and grid data
"""
type LucasTree
    gamma::Real       # coefficient of risk aversion
    beta::Real        # Discount factor in (0, 1)
    alpha::Real       # Correlation coefficient in the shock process
    sigma::Real       # Volatility of shock process
    phi::Distribution # Distribution for shock process
    grid::Vector      # Grid of points on which to evaluate prices.
    shocks::Vector    # Draws of the shock
    h::Vector         # The h function represented as a vector
end

"""
Constructor for the Lucas asset pricing model
"""
function LucasTree(;gamma=2.0,
                beta=0.95,
                alpha=0.9,
                sigma=0.1,
                grid_size=100)

    phi = LogNormal(0.0, sigma)
    shocks = rand(phi, 500)

    # == build a grid with mass around stationary distribution == #
    ssd = sigma / sqrt(1 - alpha^2)
    grid_min, grid_max = exp(-4 * ssd), exp(4 * ssd)
    grid = collect(linspace(grid_min, grid_max, grid_size))

    # == set h(y) = beta * int u'(G(y,z)) G(y,z) phi(dz) == #
    h = similar(grid)
    for (i, y) in enumerate(grid)
        h[i] = beta * mean((y^alpha .* shocks).^(1 - gamma))
    end

    return LucasTree(gamma,
                    beta,
                    alpha,
                    sigma,
                    phi,
                    grid,
                    shocks,
                    h)
end

"""
The approximate Lucas operator, which computes and returns updated function
Tf on the grid points.
"""
function lucas_operator(lt::LucasTree, f::Vector{Float64})

    # == unpack names == #
    grid, alpha, beta, h = lt.grid, lt.alpha, lt.beta, lt.h
    z = lt.shocks

    Tf = similar(f)
    Af = lin_interp(grid, f)

    for (i, y) in enumerate(grid)
        Tf[i] = h[i] + beta * mean(Af.(y^alpha .* z))
    end
    return Tf
end

"""
Compute the equilibrium price function associated with Lucas tree `lt`
"""
function compute_lt_price(lt::LucasTree, max_iter=500)

    # == Simplify names == #
    grid = lt.grid
    alpha, beta, gamma = lt.alpha, lt.beta, lt.gamma

    # == Create suitable initial vector to iterate from == #
    f_init = zeros(grid)

    func(f_vec) = lucas_operator(lt, f_vec)
    f = compute_fixed_point(func, f_init;
                                    max_iter=max_iter,
                                    err_tol=1e-4,
                                    verbose=false)

    # p(y) = f(y) * y^gamma
    price = f .* grid.^(gamma)

    return price
end


text:
using Documenter
using VkComputeUtils

makedocs(
    sitename = "VkComputeUtils.jl",
    format = Documenter.HTML(),
    modules = [VkComputeUtils],
    authors = "The developers of VkComputeUtils.jl",
    pages = ["Home" => "index.md", "Reference" => "functions.md"],
)

deploydocs(
    repo = "github.com/LCSB-BioCore/VkComputeUtils.jl.git",
    target = "build",
    branch = "gh-pages",
    push_preview = true,
)


text:
"GRASP" in keys(ENV) || error("\$GRASP environment variable not defined -- unable to determine the location of GRASP")
isabspath(ENV["GRASP"]) || error("\$GRASP not set to an absolute path (GRASP=$(ENV["GRASP"]))")
isdir(ENV["GRASP"]) || error("\$GRASP not pointing to a directory (GRASP=$(ENV["GRASP"]))")

open("grasp-path.jl", "w") do io
    grasp_path = abspath(ENV["GRASP"])
    write(io, "const grasp = \"$(grasp_path)\"")
end
include("grasp-path.jl")

# Creating basic GRASP input/output files, using GRASP
# ./build-testfiles.sh cds into ../test/graspfiles
println("Running ./build-testfiles.sh")
run(`./build-testfiles.sh`)

# Build the Fortran shared library
run(`make`)


text:
using CayleyDickson
using Base.Test: @test, @test_throws

@test begin
    a = CayleyDickson.randomBigFloat()
    isreal(Exo4Real(Exo3Real(a)))
end

@test begin
    a = 1
    !isreal(Exo4Real(Exo3Real(a, a, a, a)))
end

@test begin
    z = random(Exo4Real{Int})
    z == +(z)
end

@test begin
    a = CayleyDickson.randomBigFloat()
    real(Exo4Real(Exo3Real(a))) == a
end

@test begin
    io = IOBuffer()
    show(io, Exo4Real(Exo3Real(1.0, 2, 3, 4, 5, 6, 7, 8)))
    l = String(take!(io))
    r = "[1: 1.0, a: 2.0, b: 3.0, ab: 4.0, c: 5.0, ac: 6.0, bc: 7.0, (ab)c: 8.0, d: 0.0, ad: 0.0, bd: 0.0, (ab)d: 0.0, cd: 0.0, (ac)d: 0.0, (bc)d: 0.0, ((ab)c)d: 0.0]"
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigInt()
    c = CayleyDickson.randomBigInt()
    d = CayleyDickson.randomBigInt()
    l = Exo4Real(Exo3Real(a, b, c, d))
    r = Exo4Real(Exo2Real(a, b, c, d))
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigInt()
    l = Exo4Real(Exo1Real(a, b))
    r = Exo4Real(Exo3Real(a, b))
    l == r
end

@test begin
    x = random(Exo4Real{BigFloat})
    y = random(Exo4Real{BigFloat})
    x != y
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(Exo4Real{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = random(Exo3Real{BigInt})
    z = random(Exo4Real{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(Exo4Real{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    a = random(Exo3Real{BigInt})
    z = random(Exo4Real{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    x = random(Exo4Real{BigInt})
    y = random(Exo4Real{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(Exo4Real{BigInt})
    y = random(Exo4Real{BigInt})
    l = x * y
    r = y * x
    l != r
end

@test begin
    x = random(Exo4Real{BigInt})
    y = random(Exo4Real{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(Exo4Real{BigInt})
    y = random(Exo4Real{BigInt})
    z = random(Exo4Real{BigInt})
    l = (x + y) + z
    r = x + (y + z)
    l == r
end

@test begin
    x = random(Exo4Real{BigInt})
    y = random(Exo4Real{BigInt})
    z = random(Exo4Real{BigInt})
    l = (x * y) * z
    r = x * (y * z)
    l != r
end

@test begin
    x = random(Exo4Real{BigInt})
    y = random(Exo4Real{BigInt})
    l = alternatorL(x, y)
    r = zero(Exo4Real{BigInt})
    l != r
end

@test begin
    x = random(Exo4Real{BigInt})
    y = random(Exo4Real{BigInt})
    l = alternatorR(x, y)
    r = zero(Exo4Real{BigInt})
    l != r
end

@test begin
    x = random(Exo4Real{BigInt})
    conj(conj(x)) == x
end

@test begin
    x = random(Exo4Real{Rational{BigInt}})
    inv(inv(x)) == x
end

@test begin
    x = random(Exo4Real{BigInt})
    y = random(Exo4Real{BigInt})
    l = conj(x * y)
    r = conj(y) * conj(x)
    l == r
end

@test begin
    x = random(Exo4Real{BigInt})
    y = random(Exo4Real{BigInt})
    l = conj(x + y)
    r = conj(x) + conj(y)
    l == r
end

@test begin
    x = random(Exo4Real{BigInt})
    y = random(Exo4Real{BigInt})
    l = conj(x - y)
    r = conj(x) - conj(y)
    l == r
end

# @test begin
#     x = random(Exo4Real{Rational{BigInt}})
#     y = random(Exo4Real{Rational{BigInt}})
#     z = x * y
#     z / y == x
# end

# @test begin
#     x = random(Exo4Real{Rational{BigInt}})
#     y = random(Exo4Real{Rational{BigInt}})
#     z = x * y
#     x \ z == y
# end

# @test begin
#     x = random(Exo4Real{Rational{BigInt}})
#     y = random(Exo4Real{Rational{BigInt}})
#     l = inv(x * y)
#     r = inv(y) * inv(x)
#     l == r
# end

# @test begin
#     x = random(Exo4Real{Rational{BigInt}})
#     y = random(Exo4Real{Rational{BigInt}})
#     l = inv(x / y)
#     r = y / x
#     l == r
# end

# @test begin
#     x = random(Exo4Real{Rational{BigInt}})
#     y = random(Exo4Real{Rational{BigInt}})
#     l = inv(x \ y)
#     r = y \ x
#     l == r
# end

@test begin
    x = random(Exo4Real{BigInt})
    y = random(Exo4Real{BigInt})
    z = random(Exo4Real{BigInt})
    l = x * (y + z)
    r = (x * y) + (x * z)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(Exo4Real{BigInt})
    y = random(Exo4Real{BigInt})
    l = a * (x + y)
    r = (a * x) + (a * y)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(Exo4Real{BigInt})
    l = a * x
    r = x * a
    l == r
end

@test begin
    a = random(Exo3Real{BigInt})
    x = random(Exo4Real{BigInt})
    l = commutator(a, x)
    r = zero(Exo4Real{BigInt})
    l != r
end

@test begin
    x = random(Exo4Real{BigInt})
    y = random(Exo4Real{BigInt})
    z = random(Exo4Real{BigInt})
    l = (x + y) * z
    r = (x * z) + (y * z)
    l == r
end

@test begin
    x = random(Exo4Real{BigInt})
    y = random(Exo4Real{BigInt})
    z = random(Exo4Real{BigInt})
    l = x * (y - z)
    r = (x * y) - (x * z)
    l == r
end

@test begin
    x = random(Exo4Real{BigInt})
    y = random(Exo4Real{BigInt})
    z = random(Exo4Real{BigInt})
    l = (x - y) * z
    r = (x * z) - (y * z)
    l == r
end

@test begin
    x = random(Exo4Real{BigInt})
    quadrance(x) > big(0) || quadrance(x) == big(0)
end

@test begin
    x = random(Exo4Real{BigInt})
    y = random(Exo4Real{BigInt})
    l = quadrance(x * y)
    r = quadrance(x) * quadrance(y)
    l == r
end

@test_throws ErrorException begin
    inv(Exo4Real(Exo3Real(0)))
end

@test_throws ErrorException begin
    random(Exo4Real{Int}) / 0
end

@test_throws ErrorException begin
    0 \ random(Exo4Real{Int})
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(Exo4Real{Rational{BigInt}})
    l = z / a
    r = a \ z
    l == r
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(Exo4Real{Rational{BigInt}})
    l = a / z
    r = z \ a
    l == r
end

@test begin
    z = random(Exo4Real{BigInt})
    z == z + zero(Exo4Real{BigInt})
end

@test begin
    z = random(Exo4Real{BigInt})
    z == z * one(Exo4Real{BigInt})
end

@test begin
    z = random(Exo4Real{BigInt})
    z == z + zero(z)
end

@test begin
    z = random(Exo4Real{BigInt})
    z == z * one(z)
end

@test begin
    x = random(Exo4Real{BigInt})
    cloak(cloak(x)) == x
end

@test begin
    x = random(Exo4Real{BigInt})
    dagger(dagger(x)) == x
end

@test begin
    x = random(Exo4Real{BigInt})
    star(star(x)) == dagger(x)
end

@test begin
    x = random(Exo4Real{Rational{BigInt}})
    selfstar(selfstar(x)) == selfstar(x)
end

@test begin
    x = random(Exo4Real{Rational{BigInt}})
    anti_selfstar(anti_selfstar(x)) == anti_selfstar(x)
end

@test begin
    l = Exo4Real(1,2,3,4)
    r = Exo4Real(Exo2Real(1,2,3,4))
    l == r
end

@test begin
    l = Exo4Real(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17)
    r = Exo4Real(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16)
    l == r
end


text:
# EXPERIMENTAL: random braid generation
# These functions are the subject of my ongoing research. They are far from
# finalized. Use caution.

"""
Draw a uniform random braid.
"""
function randombraid(maxlen::Int, g::BraidGroup)
    l = convert(Int, floor(rand() * maxlen)) + 1
    x = Tuple{Int,Int}[]
    for i = 1:l
        idx = convert(Int, floor((g.n - 1) * rand())) + 1
        expt = rand() > 0.5 ? 1 : -1
        push!(x, (idx, expt))
    end
    simplifybraid(BraidElement(g, x))
end

"""
Sample a random braid given the draw function `f`, where `f` should return
null to terminate the braid, and a nonempty `Nullable{BraidElement}` otherwise.
"""
function randombraid(f::Function, g::BraidGroup)
    c = one(g)
    x = convert(Nullable{BraidElement}, f())
    while !isnull(x)
        c = c * get(x)
        x = convert(Nullable{BraidElement}, f())
    end
    c
end

"""
Get a sampler based on a random walk over strands with equal probability of
passing over and under. Stops with probability `dprob` and switches
strands with probability `sprob`.
"""
function getmarkovsampler(g::BraidGroup, sprob::Float64, dprob::Float64)
    cur = 1
    () -> begin
        if rand() < dprob
            return Nullable{BraidElement}()
        end
        if rand() < sprob
            if cur == (g.n - 1)
                ir = -1 * cur
                cur = cur - 1
                Nullable{BraidElement}(braidfromintrep(g, [ir]))
            elseif cur == 1
                ir = cur
                cur = cur + 1
                Nullable{BraidElement}(braidfromintrep(g, [ir]))
            elseif rand() < 0.5
                ir = cur
                cur = cur + 1
                Nullable{BraidElement}(braidfromintrep(g, [ir]))
            else
                ir = -1 * cur
                cur = cur - 1
                Nullable{BraidElement}(braidfromintrep(g, [ir]))
            end
        else
            Nullable{BraidElement}(one(g))
        end
    end
end

"""
Get a sampler based on a crude discretization of the Gaussian (for
experimentation only).
"""
function getgsampler(g::BraidGroup, sprob::Float64)
    ct = 0
    () -> begin
      if rand() < sprob
          Nullable{BraidElement}()
      else
          ir = convert(Int64, max(min(floor(randn() * (g.n / 2)), g.n - 1), -1*g.n + 1))
          if ir == 0
              Nullable{BraidElement}(one(g))
          else
              Nullable{BraidElement}(braidfromintrep(g, [ir]))
          end
      end
    end
end


text:
@inline function abs_db_db(x::DoubleFloat{T}) where {T<:IEEEFloat}
    return DoubleFloat{T}(abs_dd_dd(HILO(x)))
end

@inline function neg_db_db(x::DoubleFloat{T}) where {T<:IEEEFloat}
    return DoubleFloat{T}(neg_dd_dd(HILO(x)))
end

@inline function negabs_db_db(x::DoubleFloat{T}) where {T<:IEEEFloat}
    return DoubleFloat{T}(negabs_dd_dd(HILO(x)))
end

@inline function inv_db_db(x::DoubleFloat{T}) where {T<:IEEEFloat}
    return DoubleFloat{T}(inv_dd_dd(HILO(x)))
end

@inline function square_db_db(x::DoubleFloat{T}) where {T<:IEEEFloat}
    return DoubleFloat{T}(square_dd_dd(HILO(x)))
end

@inline function cube_db_db(x::DoubleFloat{T}) where {T<:IEEEFloat}
    return DoubleFloat{T}(cube_dd_dd(HILO(x)))
end

@inline function sqrt_db_db(x::DoubleFloat{T}) where {T<:IEEEFloat}
    return DoubleFloat{T}(sqrt_dd_dd(HILO(x)))
end

@inline function cbrt_db_db(x::DoubleFloat{T}) where {T<:IEEEFloat}
    return DoubleFloat{T}(cbrt_dd_dd(HILO(x)))
end


text:
function MCTS.RandomBandit(mdp::Walk1DMDP.Walk1D; kwargs...)
    RandomBandit(;
        enable_action_pw=true,
        check_repeat_action=true,
        kwargs...)
end
function MCTS.DPWBandit(mdp::Walk1DMDP.Walk1D; kwargs...)
    DPWBandit(;
        enable_action_pw=true,
        check_repeat_action=true,
        exploration_constant=0.25,
        k_action=2.0,
        alpha_action=0.9,
        kwargs...)
end
function MCTS.ModularSolver(mdp::Walk1DMDP.Walk1D, b::ModularBandit, seed=0;
                            kwargs...)
    rng = MersenneTwister(seed)
    ModularSolver(;
        n_iterations=100,
        bandit=b,
        depth=mdp.p.t_max,
        k_state=0.1,
        alpha_state=0.0,
        keep_tree=true,
        check_repeat_state=true,
        rng=rng,
        init_Q=0.0,
        init_N=0,
        next_action=RandomGaussian(rng),
        estimate_value=RolloutEstimator(RandomGaussian(rng)),
        #estimate_value=TrackedRolloutEstimator(rng, RandomSolver(rng)),
        track_best_path=true,
        kwargs...
       )
end
function MCTS.CBTSBandit(mdp::Walk1DMDP.Walk1D; kwargs...)
    CBTSBandit(;
        enable_action_pw=true,
        check_repeat_action=true,
        exploration_constant=1.0,
        A_max=20,
        n_proposes=100,
        log_length_scale=0.0,
        log_signal_sigma=0.0,
        log_obs_noise=-1.0,
        fit_qs=false,
        action_dims=1,
        n_sig=2.0,  #number of standard deviations for GP-UCB
        kwargs...)
end
function MCTS.CBTSDPWBandit(mdp::Walk1DMDP.Walk1D; kwargs...)
    CBTSDPWBandit(;
        enable_action_pw=true,
        check_repeat_action=true,
        exploration_constant=1.0,
        k_action=1.0,
        alpha_action=0.5,
        n_proposes=100,
        log_length_scale=0.0,
        log_signal_sigma=0.0,
        log_obs_noise=-1.0,
        fit_qs=false,
        action_dims=1,
        n_sig=2.0,  #number of standard deviations for GP-UCB
        kwargs...)
end
function MCTS.ModularSolver(mdp::Walk1DMDP.Walk1D, b::CBTSBandit, seed=0;
                            kwargs...)
    rng = MersenneTwister(seed)
    ModularSolver(;
        n_iterations=100,
        bandit=b,
        depth=mdp.p.t_max,
        k_state=0.1,
        alpha_state=0.0,
        keep_tree=true,
        check_repeat_state=true,
        rng=rng,
        init_Q=0.0,
        init_N=0,
        next_action=RandomGaussian(rng),
        estimate_value=RolloutEstimator(RandomGaussian(rng)),
        #estimate_value=TrackedRolloutEstimator(rng, RandomSolver(rng)),
        track_best_path=true,
        kwargs...
       )
end
Base.string(::Walk1D) = "Walk1D"


text:
# include("../src/integration/basisFuncTriangle.jl")
# include("../src/mesh/master2D.jl")
# include("../src/mesh/mesh2D.jl")
# include("../src/mesh/master3D.jl")

using luteos
using Base.Test

println("BASIS FUNCTIONS TEST")
include("unit/testBasisFunc.jl")

println("MESH TESTS")
println(" 2D")
include("unit/testMesh2D.jl")
println(" 3D")
include("unit/testMesh3D.jl")

println("2D MASTER TESTS")
# include("unit/testMaster2D.jl")

println("3D MASTER TESTS")
# include("unit/testMaster3D.jl")


text:
# First passage time density (and stationary rate, and ISI statistics)
# for a single LIF neuron without refractoriness!

push!(LOAD_PATH, abspath(@__DIR__,"..","src"))
using LinearAlgebra,Statistics,StatsBase
using Plots,NamedColors ; theme(:dark)
using SparseArrays
using SpikingRNNs; const global S = SpikingRNNs
using FFTW
using ProgressMeter

function onesparsemat(w::Real)
  mat=Matrix{Float64}(undef,1,1) ; mat[1,1]=w
  return sparse(mat)
end

# Time parameters
dt = 1E-3
Ttot = 6.
ntimes = floor(Int64,Ttot/dt)+1
# LIF neuron parameters
myτ = 0.5
v_th = 10.
v_r = -5.0
τrefr = 0.0 # refractoriness
τpcd = 0.2 # post synaptic current decay
myinput = 9.5
myσ = 1.0
mynt = S.NTLIF(myτ,v_th,v_r,τrefr,τpcd)
myps = S.PSLIF(mynt,1)

## inputs here
# one periodic input (freq is the inverse of period)
ωin = inv(2.0)
infun = function(t)
  α = 1.
  return myinput + α*sin(ωin*t*2π)
end
in_state = S.PSSimpleInput(S.InputSomeFunction(infun))
conn_in = S.BaseConnection(onesparsemat(1.0))
# ... plus noise
in_state_noise = S.PSSimpleInput(S.InputIndependentNormal(myσ))
##
mypop = S.Population(myps,(conn_in,S.FakeConnection()),
    (in_state,in_state_noise))

## that's it, let's make the network

myntw = S.RecurrentNetwork(dt,(mypop,))

## Recorder
krec = 10
rec_state = S.RecStateNow(myps,krec,dt,Ttot)
rec_spikes = S.RecCountSpikes(myps,dt)

## Run

times = (0:myntw.dt:Ttot)
nt = length(times)
# initial conditions
myps.state_now[1] = v_r

S.reset!.([rec_state,rec_spikes])

for (k,t) in enumerate(times)
  rec_state(t,k,myntw)
  rec_spikes(t,k,myntw)
  S.dynamics_step!(t,myntw)
end
r0 = S.get_mean_rates(rec_spikes,dt,Ttot)[1]
##
_ = let plt=plot(;xlabel="time",ylabel="voltage",leg=false)
  plot!(plt,rec_state.times,rec_state.state_now[1,:],linewidth=2)
end
##

function simulation_linear_response_coef(ω,E;dt=1E-3,Ttot::Float64=1000.0)
  # network without extra input
  in_noise = S.PSSimpleInput(S.InputIndependentNormal(myσ))
  in_constant = S.PSSimpleInput(S.InputSimpleOffset(myinput))
  mypop = S.Population(myps,
      (S.FakeConnection(),S.FakeConnection()),
      (in_constant,in_noise))

  rec_spikes = S.RecCountSpikes(myps,dt)
  times = (0:dt:Ttot)
  # initial conditions
  myps.state_now[1] = v_r
  ntw = S.RecurrentNetwork(dt,(mypop,))
  S.reset!(rec_spikes)
  @showprogress "1/2 --> " for (k,t) in enumerate(times)
    rec_spikes(t,k,ntw)
    S.dynamics_step!(t,ntw)
  end
  r0 = S.get_mean_rates(rec_spikes,dt,Ttot)[1]
  # network with extra input
  infun = function(t)
    return myinput + E*sin(ω*t*2π)
  end
  in_func = S.PSSimpleInput(S.InputSomeFunction(infun))
  conn_in = S.BaseConnection(onesparsemat(1.0))
  mypop = S.Population(myps,
      (conn_in,S.FakeConnection()),
      (in_func,in_noise))

  myps.state_now[1] = v_r
  ntw = S.RecurrentNetwork(dt,(mypop,))
  S.reset!(rec_spikes)
  @showprogress "2/2 --> " for (k,t) in enumerate(times)
    rec_spikes(t,k,ntw)
    S.dynamics_step!(t,ntw)
  end
  r1 = S.get_mean_rates(rec_spikes,dt,Ttot)[1]
  return (r1-r0)/E
end

##

t_period = range(0.1,2;length=20)
lrs = map(t_period) do _t
  simulation_linear_response_coef(inv(_t),10.0;Ttot=1500.0)
end
lrs2 = map(t_period) do _t
  simulation_linear_response_coef(inv(_t),30.0;Ttot=1500.0)
end

plot(t_period,[lrs lrs2];label=["E=10" "E=30"], marker=:star)

############################
### now follow the approach of Richardson et al
##
# start with stationary case

dv = 1E-3
v_lb = v_r - 3.0
ndv = round(Int64,(v_th - v_lb)/dv)
k_re = round(Int64,(v_r - v_lb)/dv)
ks = collect(0:ndv)
Vks = @. v_lb + ks * dv

j0s = fill(NaN,ndv+1)
p0s = fill(NaN,ndv+1)

ksrev = collect(ndv+1:-1:2)

j0s[ndv+1]=1.0
p0s[ndv+1]=0.0

for k in ksrev
  krond = k==k_re ? 1.0 : 0.0
  j0s[k-1] = j0s[k] - krond
  Gk = (Vks[k]-myinput)/(myσ^2)
  Hk = myτ*j0s[k]/(myσ^2)
  Ak = exp(dv*Gk)
  Bk = Gk==0 ? dv*Hk : Hk/Gk*(Ak-1)
  p0s[k-1] = p0s[k]*Ak + Bk
end

r0 = inv(dv*sum(p0s))
P0s = r0 .* p0s

## now the non-stationary linear addendum

function theory_linear_response!(jr,je,pr,pe,ω::Real;E=1.0)
  jr[ndv+1] = 1.0
  je[ndv+1] = 0.
  pr[ndv+1] = 0.
  pe[ndv+1] = 0.
  for k in ksrev
    krond = k==k_re ? 1.0 : 0.0
    jr[k-1] = jr[k] + im*2π*ω*pr[k] - krond
    je[k-1] = je[k] + im*2π*ω*pe[k]
    G = (Vks[k]-myinput)/(myσ^2)
    Hr = myτ*jr[k]/(myσ^2)
    He = (myτ*je[k]-P0s[k])/(myσ^2)
    Ak = exp(dv*G)
    pr[k-1] = G==0 ? pr[k] + Hr*dv : pr[k]*Ak + Hr/G*(Ak-1)
    pe[k-1] = G==0 ? pe[k] + He*dv : pe[k]*Ak + He/G*(Ak-1)
    #@show pe[k]
    #@show pe[k-1]
    @assert isfinite(pe[k-1])
  end
  #return je
  AEω = (-je[1]/jr[1])
  Pfull = @. r0*p0s + E*AEω + E*pe
  return AEω ,Pfull
end

_jr,_je,_pr,_pe =ntuple(_->similar(p0s,ComplexF64),4)

_Aω,_Pfullω = theory_linear_response!(_jr,_je,_pr,_pe,inv(6.0))

_ = let y = real.(_Pfullω),
  x=Vks
  plot(x,y)
end


text:
using Blender

Blender.render("simple")


text:
function deckreader(file)
    lines = readlines(file)
    filter!(x -> occursin(r"^# \d+x", x), lines)
    matches = match.(Ref(r"^# (\d+)x"), lines)
    counts = map(x -> parse(Int, x.captures[1]), matches)
    matches = match.(Ref(r"\(\d+\) (.+)$"), lines)
    names = String.(strip.(map(x -> convert(String, x.captures[1]), matches)))
    collect(zip(names, counts))
end

max_stars_repo_name: asieev/HearthstoneSim

max_issues_repo_name: asieev/HearthstoneSim

max_forks_repo_name: asieev/HearthstoneSim


text:
using DistributionsLH
using GLM, StatsBase, Random, Test

d = DistributionsLH;

function pretty_print(x)
    show(stdout, MIME("text/plain"), x);
    println("");
end

# The numbers here need to be nice. Otherwise it takes very large samples to get accurate cov matrices.
function make_mean_std(n :: Integer)
    meanV = collect(range(-2.0, 3.0, length = n));
    stdV  = collect(range(0.4, 1.2, length = n));
    # if length(stdV) > 2
    #     stdV[2] = 0.0;
    # end
    return meanV, stdV
end

# Make lower triangular weight matrix
function make_weight_matrix(stdV)
   n = length(stdV);
   rng = MersenneTwister(123);
   wtM = randn(rng, n,n) .- 0.5;
   for i1 = 1 : n
      if i1 < n
         wtM[i1, (i1+1) : n] .= 0.0;
      end
      if stdV[i1] == 0.0
         wtM[i1, :] .= 0.0;
         wtM[:, i1] .= 0.0;
      end
      # Diagonal elements are 1
      wtM[i1, i1] = 1.0;
   end
   return wtM
end

function mvn_basic_test(n :: Integer)
    @testset "Basics" begin
        meanV, stdV = make_mean_std(n);
        switches = MvNormalLHSwitches(meanV, stdV);
        m = MvNormalLH(switches);
        @test isapprox(d.means(m), meanV)
        println(m);

        wtM = make_weight_matrix(stdV);
        @test check_weight_matrix(m, wtM)
        covM = cov_matrix(m, wtM);
        @test check_cov_matrix(m, covM)
    end
end

function drawing_test(n :: Integer)
    meanV, stdV = make_mean_std(n);
    switches = MvNormalLHSwitches(meanV, stdV);
    m = MvNormalLH(switches);
    rng = MersenneTwister(43);
    nObs = 50_000;
    wtM = make_weight_matrix(stdV);
    drawM = draw_from_weights(m, wtM, nObs, rng);
    @test isequal(size(drawM), (nObs, n_vars(m)))

    mean2V = vec(mean(drawM, dims = 1));
    @test isapprox(mean2V, meanV, atol = 0.03)
    std2V = vec(StatsBase.std(drawM, dims = 1))
    @test isapprox(std2V, stdV, rtol = 0.02)

    covM = Statistics.cov(drawM, dims = 1);
    cov2M = cov_matrix(m, wtM);

    maxDiff = maximum((abs.(covM .- cov2M) ./ max.(0.2, cov2M)));
    @test maxDiff < 0.05
    if maxDiff > 0.05
        pretty_print(covM)
        pretty_print(cov2M)
        pretty_print((covM .- cov2M) ./ max.(0.2, cov2M))
    end
end

function conditional_test()
    @testset "Conditional distribution" begin
        n = 5;
        meanV, stdV = make_mean_std(n);
        switches = MvNormalLHSwitches(meanV, stdV);
        m = MvNormalLH(switches);
        wtM = make_weight_matrix(stdV);
        covM = cov_matrix(m, wtM);

        idx2V = [2, 4];
        value2V = [0.5 -0.3];
        idx1V = [j  for j = 1 : n  if !(j ∈ idx2V)];
        n1 = length(idx1V);
        conditional_distrib(m, idx2V, value2V, covM);

        # Draw a sample, by [variable, observation]
        rng = MersenneTwister(43);
        nObs = 10_000;
        # Transpose so we have draws by [obs, var]
        drawM = draw_from_weights(m, wtM, nObs, rng);

        condMeanM = zeros(nObs, n1);
        condStdM  = zeros(nObs, n1);
        for i1 = 1 : nObs
            # Columns are variables
            value2V = drawM[i1 : i1, idx2V];
            condMeanM[i1, :], condStdM[i1,:], _ =
                conditional_distrib(m, idx2V, value2V, covM);
        end

        # Test with multiple inputs (tests matrix expansion)
        condMean2M, _ = conditional_distrib(m, idx2V, drawM[:, idx2V], covM);
        @test isapprox(condMean2M, condMeanM,  atol = 0.0001);

        # ----  Check by simulation + regression
        for iVar = 1 : n1
            # Regress y against its conditional mean
            idx1 = idx1V[iVar];
            linModel = lm(hcat(ones(nObs), condMeanM[:, iVar]),  drawM[:, idx1]);
            betaV = coef(linModel);
            seBetaV = stderror(linModel);

            # Test that intercept is 0 and slope is 1
            @test !any(abs.(betaV .- [0.0, 1.0]) .> 2.0 * seBetaV)

            # Residuals
            residV = drawM[:, idx1] .- condMeanM[:, iVar];
            std2 = std(residV);
            @test abs(std2 - mean(condStdM[:, iVar])) < 0.01
        end
	end
end

# Conditional mean weights
function cond_mean_weights_test()
    @testset "Conditional mean weights" begin
        rng = MersenneTwister(42);
        n = 7;
        # With means of 0 and conditioning values of 1, we can compare
        # conditional means and weights for conditional means
        meanV = zeros(n);
        stdV  = rand(n) * 2.0;

        m = MvNormalLH(meanV, stdV);
        wtM = make_weight_matrix(stdV);
        covM = cov_matrix(m, wtM);

        # Conditioning indices
        idx2V = [2, 4, 5];
        value2V = ones(1, length(idx2V));
        condMeanV, _ = conditional_distrib(m, idx2V, value2V, covM);
        # Compare with conditioning weights
        wtM = cond_mean_weights(m, idx2V, covM);

        @test isapprox(vec(condMeanV), vec(sum(wtM, dims = 2)), atol = 1e-4);
    end
end

@testset "MvNormal" begin
    for n ∈ [2, 5]
        mvn_basic_test(n);
        drawing_test(n);
    end
    conditional_test();
    cond_mean_weights_test();
end

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


text:
# Factor Analysis

"""
This type contains factor analysis model parameters.
"""
struct FactorAnalysis{T<:Real} <: LatentVariableDimensionalityReduction
    mean::Vector{T}       # sample mean: of length d
    W::Matrix{T}          # factor loadings matrix: of size d x p
    Ψ::Vector{T}          # noise covariance: diagonal of size d x d
end

"""
    size(M::PPCA)

Returns a tuple with values of the input dimension ``d``, *i.e* the dimension of
the observation space, and the output dimension ``p``, *i.e* the dimension of
the principal subspace.
"""
size(M::FactorAnalysis) = size(M.W)

"""
    mean(M::PPCA)

Get the mean vector (of length ``d``).
"""
mean(M::FactorAnalysis) = fullmean(size(M)[1], M.mean)

"""
    projection(M::FactorAnalysis)

Recovers principle components from the weight matrix of the model `M`.
"""
projection(M::FactorAnalysis) = svd(M.W).U

"""
    cov(M::FactorAnalysis)

Returns the covariance of the model `M`.
"""
cov(M::FactorAnalysis) = M.W*M.W' + Diagonal(M.Ψ)

"""
    var(M::FactorAnalysis)

Returns the variance of the model `M`.
"""
var(M::FactorAnalysis) = M.Ψ

"""
    loadings(M::FactorAnalysis)

Returns the factor loadings matrix of the model `M`.
"""
loadings(M::FactorAnalysis) = M.W

## use

"""
    predict(M::FactorAnalysis, x)

Transform observations `x` into latent variables. Here, `x` can be either a vector
of length `d` or a matrix where each column is an observation.
"""
function predict(m::FactorAnalysis, x::AbstractVecOrMat{T}) where {T<:Real}
    xn = centralize(x, mean(m))
    W = m.W
    WᵀΨ⁻¹ = W'*diagm(0 => 1 ./ m.Ψ)  # (q x d) * (d x d) = (q x d)
    return inv(I+WᵀΨ⁻¹*W)*(WᵀΨ⁻¹*xn)  # (q x q) * (q x d) * (d x 1) = (q x 1)
end

"""
    reconstruct(M::FactorAnalysis, z)

Approximately reconstruct observations from the latent variable given in `z`.
Here, `z` can be either a vector of length ``p`` or a matrix where each column gives
the latent variables for an observation.
"""
function reconstruct(m::FactorAnalysis, z::AbstractVecOrMat{T}) where {T<:Real}
    W  = m.W
    # ΣW(W'W)⁻¹z+μ = ΣW(W'W)⁻¹W'Σ⁻¹(x-μ)+μ = Σ(WW⁻¹)((W')⁻¹W')Σ⁻¹(x-μ)+μ
    # = ΣΣ⁻¹(x-μ)+μ = (x-μ)+μ = x
    return cov(m)*W*inv(W'W)*z .+ mean(m)
end

## show

function show(io::IO, M::FactorAnalysis)
    i,o = size(M)
    print(io, "Factor Analysis(indim = $i, outdim = $o)")
end

## core algorithms

"""
    faem(S, mean, n; ...)

Performs factor analysis using an expectation-maximization algorithm for a given sample covariance matrix `S`[^2].

**Parameters**
- `S`: The sample covariance matrix.
- `mean`: The mean vector of original samples, which can be a vector of length ``d``,
or an empty vector indicating a zero mean.
- `n`: The number of observations.

Returns the resultant [`FactorAnalysis`](@ref) model.

**Note:** This function accepts two keyword arguments: `maxoutdim`,`tol`, and `maxiter`.
"""
function faem(S::AbstractMatrix{T}, mv::Vector{T}, n::Int;
             maxoutdim::Int=size(X,1)-1,
             tol::Real=1.0e-6,   # convergence tolerance
             maxiter::Integer=1000) where T<:Real

    d = size(S,1)
    q = maxoutdim
    W = Matrix{T}(I,d,q)
    Ψ = fill(T(0.01),d)

    L_old = 0.
    for c in 1:maxiter
        # EM-steps
        Ψ⁻¹W = diagm(0 => 1 ./ Ψ)*W
        SΨ⁻¹W = S*Ψ⁻¹W
        H = SΨ⁻¹W*inv(I + W'*Ψ⁻¹W)
        W⁺ = SΨ⁻¹W*inv(I + H'*Ψ⁻¹W)
        Ψ⁺ = diag(S - H*(W⁺)')
        # new parameters
        W = W⁺
        Ψ = Ψ⁺
        # log likelihood
        Ψ⁻¹ = diagm(0 => 1 ./ Ψ)
        WᵀΨ⁻¹ = W'*Ψ⁻¹
        logdetΣ = sum(log, Ψ) + logabsdet(I + WᵀΨ⁻¹*W)[1]
        Σ⁻¹ = Ψ⁻¹ - Ψ⁻¹*W*inv(I + WᵀΨ⁻¹*W)*WᵀΨ⁻¹
        L = (-n/2)*(d*log(2π) + logdetΣ + tr(Σ⁻¹*S))
        @debug "Likelihood" iter=c L=L ΔL=abs(L_old - L)
        if abs(L_old - L) < tol
            break
        end
        L_old = L
    end

    return FactorAnalysis(mv, W, Ψ)
end

"""
    facm(S, mean, n; ...)

Performs factor analysis using a fast conditional maximization algorithm for a given sample covariance matrix `S`[^3].

**Parameters**
- `S`: The sample covariance matrix.
- `mean`: The mean vector of original samples, which can be a vector of length ``d``,
or an empty vector indicating a zero mean.
- `n`: The number of observations.

Returns the resultant [`FactorAnalysis`](@ref) model.

**Note:** This function accepts two keyword arguments: `maxoutdim`,`tol`, `maxiter`, and `η`.
"""
function facm(S::AbstractMatrix{T}, mv::Vector{T}, n::Int;
             maxoutdim::Int=size(X,1)-1,
             tol::Real=1.0e-6,   # convergence tolerance
             η = tol,            # variance low bound
             maxiter::Integer=1000) where T<:Real

    d = size(S,1)

    q = maxoutdim
    W = Matrix{T}(I,d,q)
    Ψ = fill(T(0.01),d)
    V = zeros(T,q)
    eᵢeᵢ = zeros(T,d,d)
    Bᵢ⁻¹ = zeros(T,d,d)
    addconst = d*log(2π)

    L_old = 0.
    for c in 1:maxiter
        # CM-step 1
        Ψ⁻ʰ = diagm(0 => 1 ./ sqrt.(Ψ))
        S⁺ = Ψ⁻ʰ*S*Ψ⁻ʰ

        F = eigen(S⁺)
        λ = real(F.values)
        ord = sortperm(λ, rev=true)
        λ = λ[ord]

        q′ = λ[q] > 1 ? q : findlast(λ .> 1)
        q′ === nothing && error("Eigenvalues less then 1")

        λq = 0.0
        @inbounds for i in 1:q′
            λq += log(λ[i]) - λ[i] + 1.
            V[i] = sqrt(λ[i] - 1.)
        end

        L = (addconst + log(prod(Ψ)) + tr(S⁺) + λq)*(-n/2)

        U = convert(Matrix{T}, F.vectors[:,ord[1:q′]])
        W = U*diagm(0 => V[1:q′]) # set new parameter

        # CM-step 2
        @inbounds for i in 1:q′
            V[i] = 1 ./ λ[i] - 1.
        end
        eBe⁻¹ = ωᵢᵗ⁺¹ = 0.0
        for i in 1:d
            if i == 1
                Bᵢ⁻¹ = U*diagm(0 => V[1:q′])*U' + I
            else
                eᵢeᵢ[i-1,i-1] = 1.
                Bᵢ⁻¹ = Bᵢ⁻¹ -  ωᵢᵗ⁺¹*Bᵢ⁻¹*eᵢeᵢ*Bᵢ⁻¹ / (1. + ωᵢᵗ⁺¹ * Bᵢ⁻¹[i-1,i-1])
                eᵢeᵢ[i-1,i-1] = 0.
            end
            eBe⁻¹ = 1 ./ Bᵢ⁻¹[i,i]
            ωᵢᵗ⁺¹ = (Bᵢ⁻¹*S⁺*Bᵢ⁻¹)[i,i] * eBe⁻¹ * eBe⁻¹ - eBe⁻¹
            Ψ[i] = max(η, (ωᵢᵗ⁺¹ + 1.)*Ψ[i]) # set new parameter
        end

        @debug "Likelihood" iter=c L=L ΔL=abs(L_old - L)
        if abs(L_old - L) < tol
            break
        end
        L_old = L
    end

    return FactorAnalysis(mv, Diagonal(sqrt.(Ψ))*W, Ψ)
end

## interface functions
"""
    fit(FactorAnalysis, X; ...)

Perform factor analysis over the data given in a matrix `X`.
Each column of `X` is an observation.
This method returns an instance of [`FactorAnalysis`](@ref).

**Keyword arguments:**

Let `(d, n) = size(X)` be respectively the input dimension and the number of observations:

- `method`: The choice of methods:
    - `:em`: use EM version of factor analysis
    - `:cm`: use CM version of factor analysis (*default*)
- `maxoutdim`: Maximum output dimension (*default* `d-1`)
- `mean`: The mean vector, which can be either of:
    - `0`: the input data has already been centralized
    - `nothing`: this function will compute the mean (*default*)
    - a pre-computed mean vector
- `tol`: Convergence tolerance (*default* `1.0e-6`)
- `maxiter`: Maximum number of iterations (*default* `1000`)
- `η`: Variance low bound (*default* `1.0e-6`)

**Notes:** This function calls [`facm`](@ref) or [`faem`](@ref) internally, depending on the choice of method.
"""
function fit(::Type{FactorAnalysis}, X::AbstractMatrix{T};
             method::Symbol=:cm,
             maxoutdim::Int=size(X,1)-1,
             mean=nothing,
             tol::Real=1.0e-6,   # convergence tolerance
             η = tol,            # variance low bound
             maxiter::Integer=1000) where T<:Real

    d, n = size(X)

    # process mean
    mv = preprocess_mean(X, mean)
    S = covm(X, isempty(mv) ? 0 : mv, 2)
    if method == :em
        M = faem(S, mv, n, maxoutdim=maxoutdim, tol=tol, maxiter=maxiter)
    elseif method == :cm
        M = facm(S, mv, n, maxoutdim=maxoutdim, tol=tol, maxiter=maxiter, η = η)
    else
        throw(ArgumentError("Invalid method name $(method)"))
    end

    return M::FactorAnalysis
end


text:
using Reinforce:CartPole, actions, reset!, Episode, finished
import Reinforce.action
using Flux, StatsBase, Plots

gr()

#Define custom policy for choosing action
mutable struct CartPolePolicy <: Reinforce.AbstractPolicy end

#Load game environment
env = CartPole()

#Parameters
EPISODES = 3000
STATE_SIZE = length(env.state)
ACTION_SIZE = length(actions(env, env.state))
MEM_SIZE = 2000
BATCH_SIZE = 32
γ = 0.95    # discount rate
ϵ = 1.0  # exploration rate
ϵ_min = 0.01
ϵ_decay = 0.995
η = 0.001   #learning rate

memory = [] #used to remember past results

#Model Architecture
model = Chain(Dense(STATE_SIZE, 24, σ), Dense(24, 24, σ), Dense(24, ACTION_SIZE))
loss(x, y) = Flux.mse(model(x), y)
opt = ADAM(params(model), η)
fit_model(dataset) = Flux.train!(loss, dataset, opt)

function remember(state, action, reward, next_state, done)
    if length(memory) == MEM_SIZE
        deleteat!(memory, 1)
    end
    push!(memory, (state, action, reward, next_state, done))
end

function action(policy::CartPolePolicy, reward, state, action)
    if rand() <= ϵ
        return rand(1:ACTION_SIZE)
    end
    act_values = model(state)
    return Flux.argmax(act_values)[1]  # returns action
end

function replay()
    global ϵ
    minibatch = sample(memory, BATCH_SIZE, replace = false)

    for (state, action, reward, next_state, done) in minibatch
        target = reward
        if !done
            target += γ * maximum(model(next_state))
        end
        target_f = model(state).data
        target_f[action, 1] = target
        dataset = zip(state, target_f)
        fit_model(dataset)
    end

    if ϵ > ϵ_min
        ϵ *= ϵ_decay
    end
end

#Render the environment
on_step(env::CartPole, niter, sars) = gui(plot(env))

function episode!(env, policy = RandomPolicy(); stepfunc = on_step, kw...)
    ep = Episode(env, policy)
    for sars in ep
        stepfunc(ep.env, ep.niter, sars)
        state, action, reward, next_state = sars
        state = reshape(state, STATE_SIZE, 1)
        next_state = reshape(next_state, STATE_SIZE, 1)
        done = finished(ep.env, next_state) #check of game is over)
        reward = !done ? reward : -1 #Penalty of -10 if game is over
        remember(state, action, reward, next_state, done)
    end
    ep.total_reward
end

for e=1:EPISODES
    reset!(env)
    total_reward = episode!(env, CartPolePolicy())
    println("Episode: $e/$EPISODES | Score: $total_reward | ϵ: $ϵ")
    if length(memory) >= BATCH_SIZE
        replay()
    end
end

max_stars_repo_name: GantZA/model-zoo

max_issues_repo_name: GantZA/model-zoo

max_forks_repo_name: GantZA/model-zoo


text:
# TODO: input and document
__DEFAULT_LBV_PRE(nlevels::Integer) = ntuple(nlevels-1) do l
    Smoothers(
        RedBlackGaussSeidel(ForwardSweep(), iter=2^(l-1)),
        RedBlackGaussSeidel(BoundarySweep(ForwardSweep()), iter=2^l),
    )
end

__DEFAULT_LBV_COARSE() = Smoothers(
    RedBlackGaussSeidel(ForwardSweep(), iter=512),
    RedBlackGaussSeidel(BackwardSweep(), iter=512),
)

__DEFAULT_LBV_POST(nlevels::Integer) = ntuple(nlevels-1) do l
    Smoothers(
        RedBlackGaussSeidel(BoundarySweep(BackwardSweep()), iter=2^l),
        RedBlackGaussSeidel(BackwardSweep(), iter=2^(l-1)),
    )
end

"""
    lbv(
        f::AbstractArray{<:AbstractFloat, N ∈ (3, 4)},
        mask::AbstractArray{Bool, 3},
        vsz::NTuple{3, Real};
        atol::Real = sqrt(eps(Float64)),
        rtol::Real = sqrt(eps(Float64)),
        maxit::Integer = maximum(size(f)),
        verbose::Bool = false
    ) -> typeof(similar(f))

Laplacian boundary value problem (LBV) [1].

The Laplacian is computed using second order central finite differences.
The resulting Poisson's equation is then solved inside an ROI (`mask`) with
homogenous Dirichlet boundary condition (BC) using a multigrid-preconditioned
conjugate gradient method. The boundary of the ROI is set such that values
outside of it (`mask = 0`) are taken as boundary points and values inside of it
(`mask = 1`) as interior points, ie. BC: `fl[!mask] = 0`.

### Arguments
- `f::AbstractArray{<:AbstractFloat, N ∈ (3, 4)}`: unwrapped (multi-echo) field/phase
- `mask::AbstractArray{Bool, 3}`: binary mask of region of interest
- `vsz::NTuple{3, Real}`: voxel size

### Keywords
- `atol::Real = sqrt(eps(Float64))`: absolute stopping tolerance
- `rtol::Real = sqrt(eps(Float64))`: relative stopping tolerance
- `maxit::Integer = maximum(size(f))`: maximum number of cg iterations
- `verbose::Bool = false`: print convergence information

### Returns
- `typeof(similar(f))`: background corrected local field/phase

### References
[1] Zhou D, Liu T, Spincemaille P, Wang Y. Background field removal by solving
    the Laplacian boundary value problem. NMR in Biomedicine. 2014 Mar;27(3):312-9.
"""
function lbv(
    f::AbstractArray{<:AbstractFloat, N},
    mask::AbstractArray{Bool, 3},
    vsz::NTuple{3, Real};
    atol::Real = sqrt(eps(Float64)),
    rtol::Real = sqrt(eps(Float64)),
    maxit::Integer = maximum(size(f)),
    verbose::Bool = false
) where {N}
    N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
    return _lbv!(tzero(f), f, mask, vsz, atol, rtol, maxit, verbose)
end

function _lbv!(
    fl::AbstractArray{<:AbstractFloat, N},
    f::AbstractArray{T, N},
    mask::AbstractArray{Bool, 3},
    vsz::NTuple{3, Real},
    atol::Real,
    rtol::Real,
    maxit::Integer,
    verbose::Bool,
) where {T<:AbstractFloat, N}
    N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
    size(fl) == size(f) || throw(DimensionMismatch())
    size(mask) == size(f)[1:3] || throw(DimensionMismatch())

    # crop to avoid unnecessary work. fl[!mask] = 0
    sz = size(mask)
    Rc = crop_indices(mask)

    # pad for finite diffy stencil
    I1, I2 = first(Rc), last(Rc)
    Rc = CartesianIndices(ntuple(ndims(Rc)) do n
        max(1, I1[n]-1) : min(sz[n], I2[n]+1)
    end)

    # TODO: performance test proper cutoff
    szc = size(Rc)
    _crop = maximum(sz .- szc) > 10

    if !_crop
        mc = mask
        fc = f
        flc = fl
    else
        mc = similar(mask, szc)
        mc = _tcopyto!(mc, @view(mask[Rc]))

        if N == 3
            fc = @view(f[Rc])
            flc = similar(fl, szc)
            flc = _tcopyto!(flc, @view(fl[Rc]))
        else
            fc = @view(f[Rc,:])
            flc = similar(fl, (szc..., size(fl, 4)))
            flc = _tcopyto!(flc, @view(fl[Rc,:]))
        end
    end

    # set boundaries
    @inbounds for t in axes(f, 4)
        fct = @view(f[Rc,t])
        @batch for I in eachindex(fct, mc)
            fct[I] *= mc[I]
        end
    end

    # Laplacian
    A = Poisson2(T, mc, vsz)
    d2fc = A*fc

    # MGPCG options
    nlevels = __DEFAULT_MGPCG_DEPTH(mc) # stc/utils/poisson_solver/mgpcg.jl

    opts = (
        presmoother  = __DEFAULT_LBV_PRE(nlevels),
        coarsesolver = __DEFAULT_LBV_COARSE(),
        postsmoother = __DEFAULT_LBV_POST(nlevels),
        atol = atol,
        rtol = rtol,
        maxit = maxit,
        verbose = verbose,
    )

    # solve
    flc = mgpcg!(flc, A, d2fc; opts...)

    if _crop
        if N == 3
            _tcopyto!(@view(fl[Rc]), flc)
        else
            _tcopyto!(@view(fl[Rc,:]), flc)
        end
    end

    return fl
end


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

using Markdown

# ╔═╡ ec524cbb-6bf3-4905-9bcf-b12b132f5a73
md"""
# Julia es rápido

Frecuentemente, se usan benchmarks para comparar lenguajes. Éstos benchmarks pueden llevar a largas discusiones; primero, para saber lo que se está midiendo y segundo para explicar sus diferencias. Estas preguntas sencillas a veces son mucho más complicadas de lo que uno se imagina.

El propósito de este notebook es para que tú puedas hacer un benchmark simple. Uno puede leer el notebook y ver que sucedió en la Macbook Pro 4-core Intel Core i7 del autor, o correrlo uno mismo.

(Este material empezó como parte de una clase que dió Steven Johnson en MIT: https://github.com/stevengj/18S096/blob/master/lectures/lecture1/Boxes-and-registers.ipynb.)
"""

# ╔═╡ 7c937466-d3e3-40a2-a783-7aee0593194c
md"""
# Esquema de este notebook
- Definir la función suma
- Implementación y benchmark de la función suma en ...
    - C
    - python (interno)
    - python (numpy)
    - python (hecho a mano)
    - Julia (interno)
    - Julia (hecho a mano)
- Resúmenes de resultados
"""

# ╔═╡ 7b4560ee-5b96-4839-b177-e241d2b0a590
md"""
# `sum`: una función fácil de entender
"""

# ╔═╡ 1a6dc528-c4f0-4728-9584-7aa2603750e7
md"""
Considera la función  **sum** `sum(a)`, la cual calcula
$$
\mathrm{sum}(a) = \sum_{i=1}^n a_i,
$$
con $n$ longitud de `a`.
"""

# ╔═╡ 88ae83e6-fde5-4d8c-aa1a-ee8c20c4b4b6
a = rand(10^7) # Vector 1D uniforme en [0,1)

# ╔═╡ f7904227-b760-4756-8bce-00e4bfa4b8ac
sum(a)

# ╔═╡ c10e3205-7dcc-4d5a-9c3a-9d5ad0ea64ff
md"""
El resultado esperado es 0.5 * 10^7, pues el promedio de cada entrada es .5
"""

# ╔═╡ 91a933be-f773-4953-a445-dc2705337cd6
md"""
# Benchmarks de algunas manera en algunos lenguajes
"""

# ╔═╡ 15a815af-8cad-44ba-a7a9-8d1dfe46c473
md"""
Julia cuenta con el paquete `BenchmarkTools.jl` para hacer benchmarks fácil y rápidos
"""

# ╔═╡ 6d50f534-c4ec-4a52-89c4-2c9e99bccc7d
#Pkg.add("BenchmarkTools")

# ╔═╡ 0b34c831-76a6-4f64-8dbf-99141c69a334
using BenchmarkTools

# ╔═╡ c5d36d38-bb99-48de-9131-0e5cb99bc9c5
md"""
#  1. El lenguaje C

C es considerado el estándar dorado: difícil para el usuario, fácil para la máquina. Estar dentro de un factor de 2 de C puede ser muy satisfactorio. Sin emargo, aún dentro de C, existen muchos tipos de optimizaciones posibles que un usuario de C promedio puedo o no aprovechar.

Su autor no habla C, entonces no va a leer la celda siguiente, pero será feliz en saber que puedes poner código de C en una sesión de Julia, compilarlo, y correrlo. Nota que `"""` denota una cadena de varias líneas.
"""

# ╔═╡ 8fc6e2d0-9547-4c76-baa3-e57fa0d1fd8f
C_code = """
#include <stddef.h>
double c_sum(size_t n, double *X) {
    double s = 0.0;
    for (size_t i = 0; i < n; ++i) {
        s += X[i];
    }
    return s;
}
"""

const Clib = tempname()   # Haz un directorio temporario

# compila a una biblioteca compartida pipeando C_code a gcc
# (funciona sólo con gcc instalado):

open(`gcc -fPIC -O3 -msse3 -xc -shared -o $(Clib * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code)
end

# define una funcion de Julia que llama a la función de C:
c_sum(X::Array{Float64}) = ccall(("c_sum", Clib), Float64, (Csize_t, Ptr{Float64}), length(X), X)

# ╔═╡ 33d262a0-1ff7-4acd-87b1-3040f5a185fd
c_sum(a)

# ╔═╡ 94f2bead-9d31-4d39-9f84-582c496b6baa
c_sum(a) ≈ sum(a) # teclea \approx y luego <TAB> para obtener el símbolo ≈

# ╔═╡ fea09a3c-ea4d-47a0-a21b-64f9ebebc7e3
≈  # alias para la función `isapprox`

# ╔═╡ 3adf64d4-51d4-4b09-b6cf-99ec92f1fbf2
?isapprox

# ╔═╡ 0e21a579-0f38-4e4c-bcaf-dcd2f0ef5ffe
md"""
¡Ahora podemos correr el benchmark directo desde Julia!
"""

# ╔═╡ e122b029-882e-4ba0-9fd8-f52bf8b9f355
c_bench = @benchmark c_sum($a)

# ╔═╡ 5fb9fd14-211d-4c0b-915c-33f56dc55a46
println("C: Tiempo más rápido fue $(minimum(c_bench.times) / 1e6) msec")

# ╔═╡ 831a1618-5a4a-420d-a10f-c1a79519b8db
d = Dict()  # un diccionario a.k.a un arreglo asociativo
d["C"] = minimum(c_bench.times) / 1e6  # en milisegundos
d

# ╔═╡ d176bf1d-f920-43ec-a0f5-258f66939c6b
using Plots
gr()

# ╔═╡ 5ef530cd-ee1a-43da-a62b-a4d413c6d1f3
t = c_bench.times / 1e6 # tiempos en milisegundos
m, σ = minimum(t), std(t)

histogram(t, bins=500,
    xlim=(m - 0.01, m + σ),
    xlabel="milliseconds", ylabel="count", label="")

# ╔═╡ 7eedc060-e5f4-47a5-9430-8f12784f8afa
md"""
# 2. Python y `sum` interno
"""

# ╔═╡ d3456c29-d8f3-45d2-9147-dfd228dc5779
md"""
El paquete `PyCall` provee un interfaz de Julia a Python
"""

# ╔═╡ 32033011-8f48-45a7-bb3d-5977bff9f9b7
#Pkg.add("PyCall")

# ╔═╡ 4c6dbf22-bdb6-46c1-816c-de7e5d43d60d
using PyCall

# ╔═╡ 60e548a1-0224-4570-bdae-3361fb5cd828
# Llama una función de bajo nivel de PyCall para obtener una lista de Python
# porque por default PyCall convertirá a un arreglo de NumPy en vez (benchmarkeamos NumPy más abajo)

apy_list = PyCall.array2py(a, 1, 1)

# obtener el "sum" interno de Python
pysum = pybuiltin("sum")

# ╔═╡ c8a7a4d1-5c23-41ce-866b-ee519c411068
pysum(a)

# ╔═╡ 690c3eb5-b885-4c8b-970e-d53335be6e0f
pysum(a) ≈ sum(a)

# ╔═╡ c09c9a8e-556b-42da-8fcb-ffd4d682fda7
py_list_bench = @benchmark $pysum($apy_list)

# ╔═╡ 0f2c481d-32e5-4b5c-a885-b6ab9748fb7a
d["Python interno"] = minimum(py_list_bench.times) / 1e6
d

# ╔═╡ 4a64a3ce-f2fd-4d1f-b66e-b888ed2ce2f9
md"""
# 3. Python: `numpy`

## Aprovechar arquitectura "SIMD" pero sólo cuando funciona

`numpy` es una biblioteca de C optimizada, llamable desde Python, que se puede instalar en Julia haciendo:
"""

# ╔═╡ 874ebdb9-a746-4547-bdba-34d8d94ae9d8
using Conda
#Conda.add("numpy")

# ╔═╡ 92d796e9-1a66-40e3-b9fb-6dba5b322428
numpy_sum = pyimport("numpy")["sum"]
apy_numpy = PyObject(a) # convierte a un arreglo de NumPy por default

py_numpy_bench = @benchmark $numpy_sum($apy_numpy)

# ╔═╡ b68015f2-04ba-4233-8442-4a67b3231b98
numpy_sum(apy_list) # python thing

# ╔═╡ ac78ba96-cb1b-4acb-b866-e965acd3a240
numpy_sum(apy_list) ≈ sum(a)

# ╔═╡ c694f08f-f21d-49c6-9584-3c56a559256c
d["Python numpy"] = minimum(py_numpy_bench.times) / 1e6
d

# ╔═╡ da9cea03-a889-44cb-9f9d-fe62e1f7f2eb
md"""
# 4. Python, hecho a mano
"""

# ╔═╡ 7c70c208-bf85-4b53-8e68-a3cbf7bdaad8
py"""
def py_sum(a):
    s = 0.0
    for x in a:
        s = s + x
    return s
"""

sum_py = py"py_sum"

# ╔═╡ ffea97e6-c458-4887-af2f-1614e772b12a
py_hand = @benchmark $sum_py($apy_list)

# ╔═╡ 6c821697-2bff-48ef-82ea-afbaf2d72fed
sum_py(apy_list)

# ╔═╡ fd6f2717-54f4-4642-8899-ecfc61b38408
sum_py(apy_list) ≈ sum(a)

# ╔═╡ 40bce87d-2c98-4b99-8c6f-e24ae1bf89f6
d["Python hecho a mano"] = minimum(py_hand.times) / 1e6
d

# ╔═╡ 5d21c3fb-06dc-4605-a161-5b2ae0661e58
md"""
# 5. Julia (interno)

   ## Escrito directo en Julia, ¡no en C!
"""

# ╔═╡ 622cda9b-81e6-4a9f-999a-3340b5e38d58
@which sum(a)

# ╔═╡ c2316f62-e28b-42ea-8fbf-141842a8ad3c
j_bench = @benchmark sum($a)

# ╔═╡ 5fc5090a-8e35-4b16-81d6-61c72f2eed9f
d["Julia interno"] = minimum(j_bench.times) / 1e6
d

# ╔═╡ 09dfaec6-229c-43c7-b09d-ff3e523f5117
md"""
# 6. Julia (hecho a mano)
"""

# ╔═╡ adfa1a38-1835-49fe-8e7d-d15fba49be87
function mysum(A)
    s = 0.0  # s = zero(eltype(A))
    for a in A
        s += a
    end
    s
end

# ╔═╡ f88824c1-6f59-478c-accd-a6ec70d09861
j_bench_hand = @benchmark mysum($a)

# ╔═╡ d9aad37c-2a57-44aa-aafc-e0d9ffffa185
d["Julia hecho a mano"] = minimum(j_bench_hand.times) / 1e6
d

# ╔═╡ 3c676aa8-8e94-4745-abca-4402452f684b
md"""
# Summary
"""

# ╔═╡ 13a1f41f-24a4-4660-a112-3718b1ae1265
for (key, value) in sort(collect(d))
    println(rpad(key, 20, "."), lpad(round(value, 1), 8, "."))
end

# ╔═╡ bc697e8a-63a2-4e2f-bf4d-f4e7a34910eb
for (key, value) in sort(collect(d), by=x->x[2])
    println(rpad(key, 20, "."), lpad(round(value, 2), 10, "."))
end

# ╔═╡ Cell order:
# ╟─ec524cbb-6bf3-4905-9bcf-b12b132f5a73
# ╟─7c937466-d3e3-40a2-a783-7aee0593194c
# ╟─7b4560ee-5b96-4839-b177-e241d2b0a590
# ╟─1a6dc528-c4f0-4728-9584-7aa2603750e7
# ╠═88ae83e6-fde5-4d8c-aa1a-ee8c20c4b4b6
# ╠═f7904227-b760-4756-8bce-00e4bfa4b8ac
# ╟─c10e3205-7dcc-4d5a-9c3a-9d5ad0ea64ff
# ╟─91a933be-f773-4953-a445-dc2705337cd6
# ╟─15a815af-8cad-44ba-a7a9-8d1dfe46c473
# ╠═6d50f534-c4ec-4a52-89c4-2c9e99bccc7d
# ╠═0b34c831-76a6-4f64-8dbf-99141c69a334
# ╟─c5d36d38-bb99-48de-9131-0e5cb99bc9c5
# ╠═8fc6e2d0-9547-4c76-baa3-e57fa0d1fd8f
# ╠═33d262a0-1ff7-4acd-87b1-3040f5a185fd
# ╠═94f2bead-9d31-4d39-9f84-582c496b6baa
# ╠═fea09a3c-ea4d-47a0-a21b-64f9ebebc7e3
# ╠═3adf64d4-51d4-4b09-b6cf-99ec92f1fbf2
# ╟─0e21a579-0f38-4e4c-bcaf-dcd2f0ef5ffe
# ╠═e122b029-882e-4ba0-9fd8-f52bf8b9f355
# ╠═5fb9fd14-211d-4c0b-915c-33f56dc55a46
# ╠═831a1618-5a4a-420d-a10f-c1a79519b8db
# ╠═d176bf1d-f920-43ec-a0f5-258f66939c6b
# ╠═5ef530cd-ee1a-43da-a62b-a4d413c6d1f3
# ╟─7eedc060-e5f4-47a5-9430-8f12784f8afa
# ╟─d3456c29-d8f3-45d2-9147-dfd228dc5779
# ╠═32033011-8f48-45a7-bb3d-5977bff9f9b7
# ╠═4c6dbf22-bdb6-46c1-816c-de7e5d43d60d
# ╠═60e548a1-0224-4570-bdae-3361fb5cd828
# ╠═c8a7a4d1-5c23-41ce-866b-ee519c411068
# ╠═690c3eb5-b885-4c8b-970e-d53335be6e0f
# ╠═c09c9a8e-556b-42da-8fcb-ffd4d682fda7
# ╠═0f2c481d-32e5-4b5c-a885-b6ab9748fb7a
# ╟─4a64a3ce-f2fd-4d1f-b66e-b888ed2ce2f9
# ╠═874ebdb9-a746-4547-bdba-34d8d94ae9d8
# ╠═92d796e9-1a66-40e3-b9fb-6dba5b322428
# ╠═b68015f2-04ba-4233-8442-4a67b3231b98
# ╠═ac78ba96-cb1b-4acb-b866-e965acd3a240
# ╠═c694f08f-f21d-49c6-9584-3c56a559256c
# ╟─da9cea03-a889-44cb-9f9d-fe62e1f7f2eb
# ╠═7c70c208-bf85-4b53-8e68-a3cbf7bdaad8
# ╠═ffea97e6-c458-4887-af2f-1614e772b12a
# ╠═6c821697-2bff-48ef-82ea-afbaf2d72fed
# ╠═fd6f2717-54f4-4642-8899-ecfc61b38408
# ╠═40bce87d-2c98-4b99-8c6f-e24ae1bf89f6
# ╟─5d21c3fb-06dc-4605-a161-5b2ae0661e58
# ╠═622cda9b-81e6-4a9f-999a-3340b5e38d58
# ╠═c2316f62-e28b-42ea-8fbf-141842a8ad3c
# ╠═5fc5090a-8e35-4b16-81d6-61c72f2eed9f
# ╟─09dfaec6-229c-43c7-b09d-ff3e523f5117
# ╠═adfa1a38-1835-49fe-8e7d-d15fba49be87
# ╠═f88824c1-6f59-478c-accd-a6ec70d09861
# ╠═d9aad37c-2a57-44aa-aafc-e0d9ffffa185
# ╟─3c676aa8-8e94-4745-abca-4402452f684b
# ╠═13a1f41f-24a4-4660-a112-3718b1ae1265
# ╠═bc697e8a-63a2-4e2f-bf4d-f4e7a34910eb

max_stars_repo_name: miguelraz/JuliaTutorials

max_issues_repo_name: miguelraz/JuliaTutorials

max_forks_repo_name: miguelraz/JuliaTutorials


text:
#!/usr/bin/env julia

using NarrativeTest

# Workaround for https://github.com/JuliaLang/julia/pull/43787.
function Base.display(d::TextDisplay, M::MIME"text/plain", @nospecialize x)
    show(d.io, M, x)
end

# Handle filesystem paths under Windows.
subs = NarrativeTest.common_subs()
if Sys.iswindows()
    drive = splitdrive(pwd())[1]
    pushfirst!(subs, r" /?…/" => SubstitutionString(" $drive\\\\…\\\\"))
end

NarrativeTest.runtests(subs=subs)