text:
export TransfMatrix
export rotatex, rotatey, rotatez, rotate
export translate
export scalex, scaley, scalez, scale

import Base: *, ≈, show

################################################################################

"""
    struct TransfMatrix

A homogeneous transformation matrix. It contains the following fields:
- `m`
- `invm`

The multiplication `*` operation is defined on pairs of `TransMatrix`
objects. You can also use the approximate comparison operator `≈`
(useful for tests).
"""
struct TransfMatrix
    "4×4 matrix representing the transformation"
    m::StaticMatrix{4, 4, Float64}
    "4×4 matrix representing the inverse of the transformation"
    invm::StaticMatrix{4, 4, Float64}

    TransfMatrix() = new(SMatrix{4, 4}(Matrix{Float64}(I, 4, 4)),
                         SMatrix{4, 4}(Matrix{Float64}(I, 4, 4)))
    TransfMatrix(m) = new(m, inv(m))
    TransfMatrix(m, invm) = new(m, invm)
end

Base.:*(tr1::TransfMatrix, tr2::TransfMatrix) = TransfMatrix(tr1.m * tr2.m,
                                                             tr2.invm * tr1.invm)

Base.:≈(tr1::TransfMatrix, tr2::TransfMatrix) = (tr1.m ≈ tr2.m) && (tr1.invm ≈ tr2.invm)

Base.show(io::IO, tr::TransfMatrix) = println("TransfMatrix($(tr.m))")

################################################################################

"""
    identity_transform()

Return a `TransfMatrix` object that represents the identity transformation,
i.e., no transformation at all.
"""
function identity_transform()
    TransfMatrix()
end

################################################################################

"""
    translate(x, y, z)

Return a `TransfMatrix` object that represents a translation operation by
`x`, `y`, and `z` along the three axes.
"""
translate(x, y, z) = TransfMatrix(
    SMatrix{4, 4, Float64}([1 0 0 x;
                            0 1 0 y;
                            0 0 1 z;
                            0 0 0 1]),
    SMatrix{4, 4, Float64}([1 0 0 -x;
                            0 1 0 -y;
                            0 0 1 -z;
                            0 0 0 1])
)

"""
    translate(v::Vec3d)

Return a `TransfMatrix` object that represents a translation operator by
vector `v`.
"""
translate(v::Vec3d) = translate(v[1], v[2], v[3])

################################################################################

rotatex(angle) = TransfMatrix(
    SMatrix{4, 4, Float64}([1 0 0 0;
                            0 cos(angle) -sin(angle) 0;
                            0 sin(angle) cos(angle) 0;
                            0 0 0 1]),
    SMatrix{4, 4, Float64}([1 0 0 0;
                            0 cos(-angle) -sin(-angle) 0;
                            0 sin(-angle) cos(-angle) 0;
                            0 0 0 1])
)
rotatey(angle) = TransfMatrix(
    SMatrix{4, 4, Float64}([cos(angle) 0 sin(angle) 0;
                            0 1 0 0;
                            -sin(angle) 0 cos(angle) 0;
                            0 0 0 1]),
    SMatrix{4, 4, Float64}([cos(-angle) 0 sin(-angle) 0;
                            0 1 0 0;
                            -sin(-angle) 0 cos(-angle) 0;
                            0 0 0 1])
)
rotatez(angle) = TransfMatrix(
    SMatrix{4, 4, Float64}([cos(angle) -sin(angle) 0 0;
                            sin(angle) cos(angle) 0 0;
                            0 0 1 0;
                            0 0 0 1]),
    SMatrix{4, 4, Float64}([cos(-angle) -sin(-angle) 0 0;
                            sin(-angle) cos(-angle) 0 0;
                            0 0 1 0;
                            0 0 0 1])
)

"""
    rotatex(angle)
    rotatey(angle)
    rotatez(angle)

Return a `TransfMatrix` representing a rotation around one of the x, y, z axes.
The angle is expressed in radians.
"""
rotatex, rotatey, rotatez

"""
    rotate(angle, axis)

Return a `TransfMatrix` representing a rotation around an arbitrary angle. The
axis `axis` is assumed to have length one, and `angle` must be expressed in
radians.
"""
function rotate(angle, axis)
    cs, sn = cos(angle), sin(angle)
    ux, uy, uz = axis[1], axis[2], axis[3]
    m = SMatrix{4, 4, Float64}(
        cs + ux^2 * (1 - cs),
        uy * ux * (1 - cs) + uz * sn,
        uz * ux * (1 - cs) - uy * sn,
        0,
        ux * uy * (1 - cs) - uz * sn,
        cs + uy^2 * (1 - cs),
        uz * uy * (1 - cs) + ux * sn,
        0,
        ux * uz * (1 - cs) + uy * sn,
        uy * uz * (1 - cs) - ux * sn,
        cs + uz^2 * (1 - cs),
        0,
        0, 0, 0, 1,
    )

    TransfMatrix(m, inv(m))
end

################################################################################

scalex(f) = TransfMatrix(
    SMatrix{4, 4, Float64}([f 0 0 0;
                            0 1 0 0;
                            0 0 1 0;
                            0 0 0 1]),
    SMatrix{4, 4, Float64}([1/f 0 0 0;
                            0 1 0 0;
                            0 0 1 0;
                            0 0 0 1])
)

scaley(f) = TransfMatrix(
    SMatrix{4, 4, Float64}([1 0 0 0;
                            0 f 0 0;
                            0 0 1 0;
                            0 0 0 1]),
    SMatrix{4, 4, Float64}([1 0 0 0;
                            0 1/f 0 0;
                            0 0 1 0;
                            0 0 0 1])
)

scalez(f) = TransfMatrix(
    SMatrix{4, 4, Float64}([1 0 0 0;
                            0 1 0 0;
                            0 0 f 0;
                            0 0 0 1]),
    SMatrix{4, 4, Float64}([1 0 0 0;
                            0 1 0 0;
                            0 0 1/f 0;
                            0 0 0 1])
)

"""
    scalex(f)
    scaley(f)
    scalez(f)

Return a `TransfMatrix` object that represents a scaling operation along
one of the three axes `x`, `y`, or `z`. Use `scale` if you want a generic
scaling operation along the three axes.
"""
scalex, scaley, scalez

"""
    scale(fx, fy, fz)

Return a `TransfMatrix` object that represents a scaling operation along
the three axes `x`, `y`, and `z`. Use `scalex`, `scaley`, and `scalez` if
you are only interested in scaling along *one* of the three axes.
"""
scale(fx, fy, fz) = TransfMatrix(
    SMatrix{4, 4, Float64}([fx 0 0 0;
                            0 fy 0 0;
                            0 0 fz 0;
                            0 0 0 1]),
    SMatrix{4, 4, Float64}([1/fx 0 0 0;
                            0 1/fy 0 0;
                            0 0 1/fz 0;
                            0 0 0 1])
)

################################################################################

*(tr::TransfMatrix, v::Vec3d) = tr.m[1:3, 1:3] * v

function *(tr::TransfMatrix, r::Ray)

    # The origin point should be translated too
    neworig = vec3d(tr.invm[1, 1] * r.origin[1] +
                    tr.invm[1, 2] * r.origin[2] +
                    tr.invm[1, 3] * r.origin[3] +
                    tr.invm[1, 4],
                    tr.invm[2, 1] * r.origin[1] +
                    tr.invm[2, 2] * r.origin[2] +
                    tr.invm[2, 3] * r.origin[3] +
                    tr.invm[2, 4],
                    tr.invm[3, 1] * r.origin[1] +
                    tr.invm[3, 2] * r.origin[2] +
                    tr.invm[3, 3] * r.origin[3] +
                    tr.invm[3, 4])

    # Leave the translation off for the direction
    newdir = vec3d(tr.invm[1, 1] * r.dir[1] +
                    tr.invm[1, 2] * r.dir[2] +
                    tr.invm[1, 3] * r.dir[3],
                    tr.invm[2, 1] * r.dir[1] +
                    tr.invm[2, 2] * r.dir[2] +
                    tr.invm[2, 3] * r.dir[3],
                    tr.invm[3, 1] * r.dir[1] +
                    tr.invm[3, 2] * r.dir[2] +
                    tr.invm[3, 3] * r.dir[3])

    Ray(neworig, newdir)
end


text:
make_plots = true

using EchelleCCFs
using Statistics
using DataFrames, Query

include("read_soap_spectra_demo.jl")
#include("read_soap_spectra_gilbertson_etal_2020.jl")

λ_min = max(3950, minimum(λ) )  # Avoid Ca H & K lines
λ_max = min(6800, maximum(λ) )  # Avoid worst of telluric regions
idx_min = findfirst(x->x>=λ_min, λ)
idx_max = findlast(x->x<=λ_max, λ)
λ = λ[idx_min:idx_max]
flux = flux[idx_min:idx_max]

flux ./= mean(flux)

if make_plots
    using Plots
    plt = plot()
    #plt_idx = 1:length(λ)
    plt_idx = 1:2000
    plot!(plt, λ[plt_idx], flux[plt_idx], label=:none)
end

line_list_filename = joinpath(pkgdir(EchelleCCFs),"data","masks","G2.espresso.mas")
espresso_df = EchelleCCFs.read_linelist_espresso(line_list_filename)# # |>
                #@filter(λ_min < _.lambda < λ_max) |> DataFrame
line_list = EchelleCCFs.BasicLineList(espresso_df.lambda, espresso_df.weight)

resolution = 137000  # approximate R for EXPRESS
mask_width = EchelleCCFs.speed_of_light_mps / resolution  # m/s  Using a Gaussian CCF mask approximates effects of instrument resolution (neglects finite resolution of FTS)
mask_shape = GaussianCCFMask(mask_width)
#mask_shape = TopHatCCFMask(mask_width)   # In you wanted tophat mask for some reason
approx_v_offset_espresso_to_soap = -4.5e3  # just by eye so things aren't wildly shifted
max_bc = 30e3  #  ~ 2pi AU/year  in m/s
Δv_step = 400  # m/s arbitrary probably smaller/slower than you need
Δv_max = 30e3  # m/s arbitrary
ccf_plan = BasicCCFPlan(line_list=line_list, mask_shape=mask_shape, step=Δv_step, max=Δv_max, midpoint=approx_v_offset_espresso_to_soap,
    range_no_mask_change=max_bc)
v_grid = calc_ccf_v_grid(ccf_plan) .- approx_v_offset_espresso_to_soap
ccf = ccf_1D( λ, flux, ccf_plan)

if make_plots
    using Plots
    plot(v_grid,ccf)
    scatter(v_grid,ccf,markersize=1.2,label="")
    xlabel!("Δv (m/s)")
    ylabel!("CCF ")
end

ext: jl

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
using Plots; pyplot()
using DeSolver
using LinearAlgebra
using LaTeXStrings

# Number of collocation points
N = 50
h= pi/N

# create function
function f(u,t)
    du = zeros(length(u))
    for i in 1:length(u)
        if i ==1
            du[1] = (-2u[1] + 2u[2])/h^2
        elseif i == length(u)
            du[length(u)] = (2u[length(u)-1] - 2u[length(u)])/h^2
        else
            du[i] = (u[i-1] - 2u[i] + u[i+1])/h^2
        end
    end
    return du
end

xs=range(0, stop=pi, length=N)
u0=sin.(xs)

soltest1=odesolve(f, u0, (0.,10.), 1, method="RK4")
soltest2=odesolve(f, u0, (0.,10.), .1, method="RK4")
soltest3=odesolve(f, u0, (0.,10.), .01, method="RK4")


text:
ExprLike = Union{Expr,Number,Symbol}

## User should just use TransferFunction
immutable SisoGeneralized <: SisoTf
    expr::ExprLike
    function SisoGeneralized(expr::ExprLike)
        if isa(expr, Expr) && length(expr.args) == 3 && expr.args[1] == :(+) && expr.args[2] == 0
            #Get rid of the zero
            expr = expr.args[3]
        end
        new(expr)
    end
end

SisoGeneralized(str::AbstractString) = SisoGeneralized(parse(str))

function minreal(sys::SisoGeneralized, eps::Real=sqrt(eps()))
    warn("minreal is not implemented for generalized transferfunctions, returning same system")
    sys
end

function print_siso(io::IO, t::SisoGeneralized, var=:s)
    println(io, t.expr)
end

Base.promote_rule{T<:Real}(::Type{SisoGeneralized}, ::Type{T}) = SisoGeneralized
Base.convert(::Type{SisoGeneralized}, b::Real) = SisoGeneralized(b)

Base.zero(::Type{SisoGeneralized}) = SisoGeneralized(0)
Base.zero(::SisoGeneralized) = Base.zero(SisoGeneralized)

Base.length(t::SisoGeneralized) = error("length is not implemented for generalized transferfunctions")
Base.num(t::SisoGeneralized) = error("num is not implemented for generalized transferfunctions")
Base.den(t::SisoGeneralized) = error("den is not implemented for generalized transferfunctions")
pole(t::SisoGeneralized) = error("pole is not implemented for generalized transferfunctions")
tzero(t::SisoGeneralized) = error("tzero is not implemented for generalized transferfunctions")

#This makes sure that the function can compile once
function _preprocess_for_freqresp(sys::SisoGeneralized)
    _f = eval(:(s -> $(sys.expr)))
end

evalfr(f::Function, freq) = f(freq)
evalfr(sys::SisoGeneralized, freq) = _preprocess_for_freqresp(sys)(freq)

function lsimabstract(sys::SisoGeneralized, uin, dt, Tend)
    #TODO make sure U and P are of equal length, fix input arguments, Tend results in half time, make sure u interp is using Tend
    N = round(Int, Tend/dt) + 1
    #T=N*dt
    T = Tend
    dw = pi/T
    omega = linspace(-pi/dt, pi/dt, 2N+1)
    u = [uin; zeros(N)]
    U = fft(u)
    Pf = _preprocess_for_freqresp(sys)
    P = Complex{Float64}[evalfr(Pf, omega[i]*im) for i in 1:2N]
    y = real(ifft(fftshift(P).*U))
    t = dt*(0:N-1)
    y[1:N]
end

function SisoRational(expr::Expr)
    if length(expr.args) == 1
        error("Unexpected operator $(expr.args[1]) in expression")
    end
    if expr.args[1] == :+
        return reduce(+, map(t -> SisoRational(t), expr.args[2:end]))
    elseif expr.args[1] == :*
        return reduce(*, map(t -> SisoRational(t), expr.args[2:end]))
    elseif length(expr.args) == 2 && expr.args[1] == :-
        return - SisoRational(expr.args[2])
    elseif length(expr.args) == 3
        if expr.args[1] == :-
            return SisoRational(expr.args[2]) - SisoRational(expr.args[3])
        elseif expr.args[1] == :/
            return SisoRational(expr.args[2]) / SisoRational(expr.args[3])
        elseif expr.args[1] == :^
            return SisoRational(expr.args[2]) ^ expr.args[3]
        else
            error("Unexpected operator $(expr.args[1]) in expression")
        end
    else
        error("Could not parse \"$(expr.args[1])\" in expression")
    end
end
SisoRational(expr::Symbol) = (expr == :s ? SisoRational([1, 0],[1]) : error("invalid symbol \"$exp\" only \"s\" is allowed"))
SisoRational(expr::Number) = isa(expr,Real) ? SisoRational([expr],[1]) : error("Only real numers are allowed in transferfunction")

==(t1::SisoGeneralized, t2::SisoGeneralized) = (t1.expr == t2.expr)

isapprox(t1::SisoGeneralized, t2::SisoGeneralized; kwargs...) = error("isapprox not implemented for generalized tf")

+(t1::SisoGeneralized, t2::SisoGeneralized) = SisoGeneralized(:($(t1.expr) + $(t2.expr)))
+(t::SisoGeneralized, n::Real) = SisoGeneralized(:($(t.expr) + $n))
+(n::Real, t::SisoGeneralized) = SisoGeneralized(:($n + $(t.expr)))

-(t1::SisoGeneralized, t2::SisoGeneralized) = SisoGeneralized(:($(t1.expr) - $(t2.expr)))
-(n::Real, t::SisoGeneralized) = SisoGeneralized(:($n - $(t.expr)))
-(t::SisoGeneralized, n::Real) = SisoGeneralized(:($(t.expr) - $n))

-(t::SisoGeneralized) = SisoGeneralized(:(- $(t.expr)))

*(t1::SisoGeneralized, t2::SisoGeneralized) = SisoGeneralized(:($(t1.expr) * $(t2.expr)))
*(t::SisoGeneralized, n::Real) = SisoGeneralized(:($(t.expr) * $n))
*(n::Real, t::SisoGeneralized) = SisoGeneralized(:($n * $(t.expr)))

/(n::Real, t::SisoGeneralized) = SisoGeneralized(:($n / $(t.expr)))
/(t::SisoGeneralized, n::Real) = SisoGeneralized(:($(t.expr) / $n))
/(t1::SisoGeneralized, t2::SisoGeneralized) = SisoGeneralized(:($(t1.expr) / $(t2.expr)))


text:
# This is code for computing a transit model and derivatives integrated over
# a time step, giving the fluence in units of time (since flux is normalized to unity).

include("transit_poly_struct.jl")

# Now the version with derivatives:
function integrate_timestep_gradient!(param::Array{T,1},trans::Transit_Struct{T},t1::T,t2::T,tol::T,maxdepth::Int64,fint::Array{T,1}) where {T <: Real}
  # Keep track of maximum number of depths for recursion:
  depthmax = 0
  # Define third to save time in division:
  third = inv(3)
  # Keep track of how many evaluations of the transit function are carried out:
  neval = 0

  # When transit computation is carried out, we need to fill in the vector for the flux and derivatives
  # at the midpoint of the current refinement time range:
  function fill_flux!(tmid::T,fmid0::T,trans_mid::Transit_Struct{T},fmid::Array{T,1}) where {T <: Real}
#  fmid = Array{T}(undef,6+trans_mid.n)
  fmid[1] = fmid0   # flux
  fmid[2] = trans_mid.dfdrb[1]  # r derivative
  binv = inv(trans_mid.b)
  fmid[3] = trans_mid.dfdrb[2]*binv*param[2]^2*(param[1]-tmid)  # t0 derivative
  fmid[4] = trans_mid.dfdrb[2]*param[2]*binv*(tmid-param[1])^2  # v derivative
  fmid[5] = trans_mid.dfdrb[2]*param[3]*binv  # b0 derivative
  @inbounds for i=1:trans_mid.n+1
    fmid[5+i] = trans_mid.dfdg[i]  # d_i derivatives
  end
  return fmid
  end

  # This routine gives a simple model for the impact parameter versus time:
  function solver(time::T) where {T <: Real}
  # This predicts the impact parameter as a function of time in
  # terms of the transit_struct.
  # For now, transit is approximated as a straight line, where params are given as:
  # param[1] = t0 (mid-transit time; units of days/JD)
  # param[2] = v/R_* (sky velocity in terms of radius of the star - units are day^{-1}
  # param[3] = b (impact parameter at time t0)
  return sqrt(param[3]^2+(param[2]*(time-param[1]))^2)
  end

  # Function to define the vector integration in integrator:
  function transit_flux_derivative(tmid::T,fmid::Array{T,1}) where {T <: Real}
#  function transit_flux_derivative(tmid::T) where {T <: Real}
#  if VERSION >=  v"0.7"
#    fmid = Array{T}(undef,6+trans.n)
#  else
#    fmid = Array{T}(6+trans.n)
#  end
  # Add the impact parameter to the transit structure:
  trans.b = solver(tmid)
  # Carry out the transit computation:
  fmid0 = transit_poly_d!(trans); neval += 1
  # Then, fill in the flux and derivatives in the fmid vector:
  fill_flux!(tmid,fmid0,trans,fmid)
#  println("b: ",trans.b," f/df: ",fmid)
  return fmid
  end

#  # Inner integrator which uses adaptive Simpson's method thanks to DFM:
#  function inner(func::Function,x0::T,dx::T,ym::Array{T,1},y0::Array{T,1},yp::Array{T,1},tol::T,max_depth::Int64,min_depth::Int64,depth::Int64) where {T <: Real}
#  # For each depth we need a left and a right copy of: x_m, val_m; x_p, val_p, int_mp, pred
#  x_m = x0 - 0.5*dx; val_m = func(x_m)
#  x_p = x0 + 0.5*dx; val_p = func(x_p)
#  int_mp = 0.5*dx * (2*y0+ ym + yp + 4*(val_m+val_p))*third
#  pred   =     dx * (4*y0+ ym + yp)*third
#  if (depth < min_depth || (depth < max_depth && maximum(abs,pred-int_mp)> tol))
#    return inner(func, x_m, 0.5*dx, ym, val_m, y0, tol, max_depth, min_depth, depth+1) +
#           inner(func, x_p, 0.5*dx, y0, val_p, yp, tol, max_depth, min_depth, depth+1)
#  end
#  return int_mp
#  end

  # Inner integrator (this version avoids broadcast operations; however, we don't need this for now).
  function inner(func::Function,x0::T,dx::T,ym::Array{T,1},y0::Array{T,1},yp::Array{T,1},tol::T,max_depth::Int64,min_depth::Int64,depth::Int64) where {T <: Real}
  if VERSION >= v"0.7"
    int_p = Array{T}(undef,6+trans.n); int_m = Array{T}(undef,6+trans.n); pred = Array{T}(undef,6+trans.n)
  else
    int_p = Array{T}(6+trans.n); int_m = Array{T}(6+trans.n); pred = Array{T}(6+trans.n)
  end
  x_m = x0 - 0.5*dx;
  val_m = Array{T}(6+trans.n)
#  val_m = func(x_m);
  func(x_m,val_m)
  x_p = x0 + 0.5*dx;
#  val_p = func(x_p);
  val_p = Array{T}(6+trans.n)
  func(x_p,val_p)
  maxdiff = zero(T)
  @inbounds for i=1:6+trans.n
    int_m[i] = 0.5*dx * (4*val_m[i] + ym[i] + y0[i])*third
    int_p[i] = 0.5*dx * (4*val_p[i] + yp[i] + y0[i])*third
    pred[i]  =     dx * (4*y0[i]    + ym[i] + yp[i])*third
    diffabs = abs(pred[i]-int_m[i]-int_p[i])
    if  diffabs > maxdiff
      maxdiff = diffabs
    end
  end
  if (depth < min_depth || (depth < max_depth && maxdiff > tol))
    int_m = inner(func, x_m, 0.5*dx, ym, val_m, y0, tol, max_depth, min_depth, depth+1);
    int_p = inner(func, x_p, 0.5*dx, y0, val_p, yp, tol, max_depth, min_depth, depth+1);
  end
  @inbounds for i=1:6+trans.n
     int_m[i] += int_p[i]
  end
  return int_m::Array{T,1}
  end

  # Outer integrators - this dispatch computes function at upper and lower limits:
  function outer!(func::Function,lower::T,upper::T,tol::T;max_depth=50,min_depth=0) where {T <: Real}
#  ym = func(lower)
  ym = Array{T}(6+trans.n)
  yp = Array{T}(6+trans.n)
  func(lower,ym)
  func(upper,yp)
  return outer!(func,lower,upper,ym,yp,tol,max_depth,min_depth)
  end

  # This version first compares trapezoid and Simpson's rule - if they are not precise
  # enough in comparison, then it begins refinement in inner recursion:
  function outer!(func::Function,lower::T,upper::T,ym::Array{T,1},yp::Array{T,1},tol::T,max_depth::Int64,min_depth::Int64) where {T <: Real}
  x0 = 0.5 * (upper + lower);
  dx = 0.5 * (upper - lower);
#  y0 = func(x0);
  y0 = Array{T}(6+trans.n)
  func(x0,y0);
  if VERSION >= v"0.7"
    int_trap = Array{T}(undef,6+trans.n)
    int_simp = Array{T}(undef,6+trans.n)
  else
    int_trap = Array{T}(6+trans.n)
    int_simp = Array{T}(6+trans.n)
  end
  maxdiff = zero(T)
  @inbounds for i=1:6+trans.n
    int_trap[i] = 0.5 * dx * (ym[i] + 2 * y0[i] + yp[i]);
    int_simp[i] = dx * (ym[i] + 4 * y0[i] + yp[i]) / 3;
    diffabs = abs(int_trap[i]-int_simp[i])
    if diffabs > maxdiff
      maxdiff = diffabs
    end
  end
  if (min_depth == 0 && maxdiff < tol)
    return int_simp;
  end
  return inner(func, x0, dx, ym, y0, yp, tol, max_depth, min_depth, 0);
  end

#  fint,ferr = hquadrature(trans.n+6,transit_flux_derivative,t1,t2,abstol=tol)
  # fint contains time-integrated flux (fluence) and derivatives:
  fint .= outer!(transit_flux_derivative,t1,t2,tol)
  # Return the number of evalutions and maximum depth for record-keeping:
return neval::Int64,depthmax::Int64
end


text:
# These oracles mimics behaviour of AlphaZero

struct AlphaRepresentation
end

struct AlphaDynamics
  gspec
end

struct AlphaPrediction
  net :: PredictionNetwork # PredictionNetwork
end

#AlphaZero-like behaviour
function forward(representation_h::AlphaRepresentation, state)
  # hidden_state = Flux.flatten(state)
  hidden_state = state
  return hidden_state
end

function (representation_h::AlphaRepresentation)(state)
  # hidden_state = state
  hidden_state = forward(representation_h, state)
  return hidden_state
end

#AlphaZero-like behaviour
function forward(dynamics_g::AlphaDynamics, hidden_state, action)
  game_g = GI.init(dynamics_g.gspec, hidden_state)
  GI.play!(game_g, action)
  reward = GI.white_reward(game_g)
  next_hidden_state = GI.current_state(game_g)
  return (reward, next_hidden_state)
end

function (dynamics_g::AlphaDynamics)(hidden_state, action)
  return forward(dynamics_g, hidden_state, action)
end

#AlphaZero-like behaviour
function forward(prediction_f::AlphaPrediction, hidden_state)
  x = Flux.flatten(hidden_state)
  forward(prediction_f.net, x)
end

function (prediction_f::AlphaPrediction)(hiddenstate)
  gspec = prediction_f.net.gspec
  x = GI.vectorize_state(gspec,hiddenstate)
  xnet = to_singletons(x)
  net_output = forward(prediction_f, xnet)
  p, v = from_singletons.(net_output)
  return (p, v[1])
end
Flux.@functor AlphaPrediction


text:
#以一个时层中所有的部件状况为一个子图，约束主要是子图间的时间关系
#也就是时层与时层的关系
using Plasmo
using Ipopt
using PlasmoTest
#using Plots

G = OptiGraph()
O = Ipopt.Optimizer

n1 = OptiGraph()
n2 = OptiGraph()
n3 = OptiGraph()
add_subgraph!(G,n1)
add_subgraph!(G,n2)
add_subgraph!(G,n3)

# “nxy” 表示 “第x个子图中的第y个结点”
n11 = @optinode(n1,n11)
n12 = @optinode(n1,n12)
n21 = @optinode(n2,n21)
n22 = @optinode(n2,n22)
n31 = @optinode(n3,n31)
n32 = @optinode(n3,n32)

#所有子图中的1结点和2结点分别相同，即n11、n21、n31相同
#Node n11 Model
@variable(n11,0 <= x11 <= 2)
@variable(n11,2 <= y11 <= 4)
@constraint(n11,x11 + y11 >= 5)

#Node n12 Model
@variable(n12,x12)
@NLnodeconstraint(n12,ref,exp(x12) >= 2)
@variable(n12,y12 >= 0)
@constraint(n12,x12 + y12 >= 3)

#Node n21 Model
@variable(n21,0 <= x21 <= 2)
@variable(n21,2 <= y21 <= 4)
@constraint(n21,x21 + y21 >= 5)

#Node n22 Model
@variable(n22,x22)
@NLnodeconstraint(n22,ref,exp(x22) >= 2)
@variable(n22,y22 >= 0)
@constraint(n22,x22 + y22 >= 3)

#Node n31 Model
@variable(n31,0 <= x31 <= 2)
@variable(n31,2 <= y31 <= 4)
@constraint(n31,x31 + y31 >= 5)

#Node n32 Model
@variable(n32,x32)
@NLnodeconstraint(n32,ref,exp(x32) >= 2)
@variable(n32,y32 >= 0)
@constraint(n32,x32 + y32 >= 3)

#Add linkconstraint
@linkconstraint(G,n11[:x11]==n12[:y12])
@linkconstraint(G,n12[:x12]==n21[:y21])
@linkconstraint(G,n12[:y12]==n22[:y22])
@linkconstraint(G,n22[:x22]==n31[:y31])
@linkconstraint(G,n22[:y22]==n32[:y32])
#Graph objective
@objective(G, Min, n12[:x12] + n22[:x22] + n32[:x32])

optimize!(G,O)
println("x11 = ",digit5(value(n11,n11[:x11])))
println("x12 = ",digit5(value(n12,n12[:x12])))
println("x21 = ",digit5(value(n21,n21[:x21])))
println("x22 = ",digit5(value(n22,n22[:x22])))
println("x31 = ",digit5(value(n31,n31[:x31])))
println("x32 = ",digit5(value(n32,n32[:x32])))
println("Objective = ", digit5(objective_value(G)))
#plt_G = plt_G = Plots.plot(G,node_labels = true, markersize=60, labelsize=30, linewidth=4, layout_options = Dict(:tol => 0.01, :iterations => 2))
#savefig("003.png")

max_stars_repo_name: BoyangSheng/Skill-studing

max_issues_repo_name: BoyangSheng/Skill-studying

max_forks_repo_name: BoyangSheng/Skill-studying


text:
struct BisectCostBottleneckSplitter{F, T}
    f::F
    ϵ::T
end

function partition_stripe(A::SparseMatrixCSC{Tv, Ti}, K, method::BisectCostBottleneckSplitter, args...; kwargs...) where {Tv, Ti}
    @inbounds begin
        (m, n) = size(A)
        f = oracle_stripe(SparseHint(), method.f, A, args...; kwargs...)
        ϵ = method.ϵ

        #=
            search returns the largest j′ such that f(j, j′) <= c, returns
            max(j, j′_lo) - 1 if no such j′ can be found
        =#
        @inline function search(j, j′_lo, j′_hi, k, c)
            j′_lo = max(j, j′_lo)
            while j′_lo <= j′_hi
                j′ = fld2(j′_lo + j′_hi)
                if f(j, j′, k) <= c
                    j′_lo = j′ + 1
                else
                    j′_hi = j′ - 1
                end
            end
            return j′_hi
        end

        spl_lo = ones(Int, K + 1)
        spl_lo[K + 1] = n + 1

        spl_hi = fill(n + 1, K + 1)
        spl_hi[1] = 1

        spl = undefs(Int, K + 1)
        spl[1] = 1
        spl[K + 1] = n + 1

        c_lo, c_hi = bound_stripe(A, K, args..., f) ./ 1

        while c_lo * (1 + ϵ) < c_hi
            c = (c_lo + c_hi) / 2
            spl[1] = 1
            chk = true
            for k = 1 : K - 1
                spl[k + 1] = search(spl[k], spl_lo[k + 1], spl_hi[k + 1], k, c)
                if spl[k + 1] < spl[k]
                    chk = false
                    spl[k + 1 : K] .= spl[k]
                    break
                end
            end
            if chk && f(spl[K], spl[K + 1], K) <= c
                c_hi = c
                spl_hi .= spl
            else
                c_lo = c
                spl_lo .= spl
            end
        end
        return SplitPartition(K, spl_hi)
    end
end

struct FlipBisectCostBottleneckSplitter{F, T}
    f::F
    ϵ::T
end

function partition_stripe(A::SparseMatrixCSC{Tv, Ti}, K, method::FlipBisectCostBottleneckSplitter, args...; kwargs...) where {Tv, Ti}
    @inbounds begin
        (m, n) = size(A)
        f = oracle_stripe(SparseHint(), method.f, A, args...; kwargs...)
        ϵ = method.ϵ

        #=
            search returns the smallest j′ such that f(j, j′) <= c, returns
            j′_hi + 1 if no such j′ can be found
        =#
        @inline function search(j, j′_lo, j′_hi, k, c)
            j′_lo = max(j, j′_lo)
            while j′_lo <= j′_hi
                j′ = fld2(j′_lo + j′_hi)
                if f(j, j′, k) <= c
                    j′_hi = j′ - 1
                else
                    j′_lo = j′ + 1
                end
            end
            return j′_lo
        end

        spl_lo = ones(Int, K + 1)
        spl_lo[K + 1] = n + 1

        spl_hi = fill(n + 1, K + 1)
        spl_hi[1] = 1

        spl = undefs(Int, K + 1)
        spl[1] = 1
        spl[K + 1] = n + 1

        c_lo, c_hi = bound_stripe(A, K, args..., f) ./ 1

        while c_lo * (1 + ϵ) < c_hi
            c = (c_lo + c_hi) / 2
            spl[1] = 1
            chk = true
            for k = 1 : K - 1
                spl[k + 1] = search(spl[k], spl_lo[k + 1], spl_hi[k + 1], k, c)
                if spl[k + 1] > n + 1
                    chk = false
                    spl[k + 1 : K] .= n + 1
                    break
                end
            end
            if chk && f(spl[K], spl[K + 1], K) <= c
                c_hi = c
                spl_lo .= spl
            else
                c_lo = c
                spl_hi .= spl
            end
        end
        return SplitPartition(K, spl_lo)
    end
end


text:
using BayesianOptimization, GaussianProcesses, Random, Statistics

# function and regret definitions
branin(x::Vector; kwargs...) = branin(x[1], x[2]; kwargs...)
branin(x1, x2; a = 1, b = 5.1/(4π^2), c = 5/π, r = 6, s = 10, t = 1/(8π),
       noiselevel = 0) =
    a * (x2 - b*x1^2 + c*x1 - r)^2 + s*(1 - t)*cos(x1) + s + noiselevel * randn()

minima(::typeof(branin)) = [[-π, 12.275], [π, 2.275], [9.42478, 2.475]], 0.397887

hartmann(x; α = [1., 1.2, 3., 3.2], A = [10 3 17 3.5 1.7 8;
                                         .05 10 17 .1 8 14;
                                         3 3.5 1.7 10 17 8;
                                         17 8 .05 10 .1 14],
         P = 1e-4 * [1312 1696 5569 124 8283 5886;
                     2329 4135 8307 3736 1004 9991;
                     2348 1451 3522 2883 3047 6650;
                     4047 8828 8732 5743 1091 381]) =
 -sum([α[i] * exp(-sum([A[i, j] * (x[j] - P[i, j])^2 for j in 1:6])) for i in 1:4])

minima(::typeof(hartmann)) = [[.20169, .150011, .476874, .275332, .311652, .6573]], -3.32237

euclidean(x, y) = √sum((x .- y).^2)
function regret(opt, func)
    mins, fmin = minima(func)
    (observed_dist = minimum(map(m -> euclidean(m, opt.observed_optimizer), mins)),
     observed_regret = abs(opt.observed_optimum - fmin),
     model_dist = minimum(map(m -> euclidean(m, opt.model_optimizer), mins)),
     model_regret = abs(Int(opt.sense) * opt.model_optimum - fmin))
end

# optimize noise-free branin

opt = BOpt(x -> branin(x, noiselevel = 0),
           ElasticGPE(2, mean = MeanConst(-10.), kernel = SEArd([0., 0.], 5.),
                      logNoise = -2., capacity = 3000),
           ExpectedImprovement(),
           MAPGPOptimizer(every = 50, noisebounds = [-4, 3],
                         kernbounds = [[-1, -1, 0], [4, 4, 10]],
                         maxeval = 40),
           [-5., 0.], [10., 15.], maxiterations = 200,
           sense = Min)
@time boptimize!(opt)
regret(opt, branin)

# optimize hartman

opt = BOpt(hartmann,
           ElasticGPE(6, mean = MeanConst(0.), kernel = Mat52Ard(zeros(6), 0.),
                      logNoise = -2., capacity = 3000),
           ExpectedImprovement(),
           MAPGPOptimizer(every = 20, noisebounds = [-4, 3],
                         kernbounds = [[-3*ones(6); -3], [4*ones(6); 3]],
                         maxeval = 100),
           zeros(6), ones(6), maxiterations = 300,
           sense = Min)
@time boptimize!(opt)
regret(opt, hartmann)

# compare model and observation optimizer on noisy branin

all_obs_regs = []
all_mod_regs = []
for _ in 1:10
    opt = BOpt(x -> branin(x, noiselevel = 1),
               ElasticGPE(2, mean = MeanConst(-50.), kernel = SEArd(zeros(2), 4.),
                          logNoise = 2., capacity = 1000),
               UpperConfidenceBound(),
               MAPGPOptimizer(every = 100, noisebounds = [-4, 3],
                              kernbounds = [[-1, -1, 0], [4, 4, 10]],
                              f_calls_limit = 40), repetitions = 5,
               [-5., 0.], [10., 15.], maxiterations = 20,
               sense = Min, verbosity = Silent)
    obs_regs = []
    mod_regs = []
    for i in 1:10
        res = boptimize!(opt);
        _, obs_reg, _, mod_reg = regret(opt, branin)
        push!(obs_regs, obs_reg)
        push!(mod_regs, mod_reg)
    end
    push!(all_obs_regs, obs_regs)
    push!(all_mod_regs, mod_regs)
end
using PGFPlotsX
x = (1:length(all_obs_regs[1])) * 20 * 5
@pgf Axis({ymode = "log", legend_entries = ["average observation regret",
                                            "average model regret"],
           legend_columns = 1, legend_pos = "north east",
           ylabel = "number of observations",
           title = "model optimizers become more accurate for noisy objectives"},
          Plot({red, very_thick, no_marks}, Coordinates(x, mean(all_obs_regs))),
          Plot({blue, very_thick, no_marks}, Coordinates(x, mean(all_mod_regs))),
          [Plot({red, thin, forget_plot}, Coordinates(x, y)) for y in all_obs_regs]...,
          [Plot({blue, thin, forget_plot}, Coordinates(x, y)) for y in all_mod_regs]...)


text:
# Plot of various Mathieu functions, inspired by https://dlmf.nist.gov/28.3
using Mathieu, PyPlot
const MPL = PyPlot.matplotlib

# get the same line colors
PyPlot.PyDict(MPL."rcParams")["axes.prop_cycle"] =
    MPL.cycler(color=["#2ca02c","#d62728","#1f77b4","#ff7f0e"])

# plot parameters
z = LinRange(0,π/2,101)
n = 0:3

for q in [1,10]
    f,ax = subplots(2,2,figsize=(10,8));
    f.suptitle("Mathieu Functions for \$q=$q\$");

    ax[1].set_title("Even π-Periodic Solutions");
    ax[1].plot(z, Mathieu.cep(n,q,z));
    ax[1].legend(string.("\$ce_",2n,"\$"));

    ax[2].set_title("Even π-Antiperiodic Solutions");
    ax[2].plot(z, Mathieu.cea(n,q,z));
    ax[2].legend(string.("\$ce_",2n.+1,"\$"));

    ax[3].set_title("Odd π-Antiperiodic Solutions");
    ax[3].plot(z, Mathieu.sea(n,q,z));
    ax[3].legend(string.("\$se_",2n.+1,"\$"));

    ax[4].set_title("Odd π-Periodic Solutions");
    ax[4].plot(z, Mathieu.sep(n,q,z));
    ax[4].legend(string.("\$se_",2n.+2,"\$"));
end


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("Parallel_Computing")

# # Distributed (or multi-core or multi-process) parallelism
#
# Julia has a built-in standard library — Distributed — that allows you to
# start and run multiple concurrent Julia processes. Imagine starting a slew
# of Julia instances and then having an easy way to run code on each and every
# one of them; that's what Distributed provides.
#
# ![](https://raw.githubusercontent.com/JuliaComputing/JuliaAcademyData.jl/master/courses/Parallel_Computing/images/Julia6x.png)

using Distributed
nprocs()

#-

addprocs(4)
@sync @everywhere workers() include("/opt/julia-1.0/etc/julia/startup.jl") # Needed just for JuliaBox
nprocs()

#-

myid()

# Now we can easily communicate with the other nodes:

r = @spawnat 2 (myid(), rand())

#-

fetch(r)

# This works kinda like an `@async` task!

@time r = @spawnat 2 (sleep(1), rand())
@time fetch(r)

# So we can repeat the same examples from tasks:

@time for i in 2:nprocs() # proc 1 is the controlling node
    @spawnat i sleep(1)
end

#-

@time @sync for i in 2:nprocs()
    @spawnat i sleep(1)
end

# Except unlike tasks, we're executing the code on a separate process — which
# can be performed on a different processor in parallel!

@everywhere function work(N)
    series = 1.0
    for i in 1:N
        series += (isodd(i) ? -1 : 1) / (i*2+1)
    end
    return 4*series
end

#-

@time work(1_000_000_000)
@time @sync for i in workers()
    @spawnat i work(1_000_000_000)
end

# Of course, this isn't very helpful. We're just performing exactly the same
# calculation on every worker... and then completely ignoring the result! Let's
# restructure our computation to be a bit more parallel friendly:

@everywhere function partial_pi(r)
    series = 0.0
    for i in r
        series += (isodd(i) ? -1 : 1) / (i*2+1)
    end
    return 4*series
end
a = partial_pi(0:999)
a, a-pi

#-

b = partial_pi(1000:9999)
(a + b), (a+b) - pi

# So now we can distribute this computation across our many workers!

r = 0:10_000_000_000
futures = Array{Future}(undef, nworkers())
@time begin
    for (i, id) in enumerate(workers())
        batch = 0:length(r)÷nworkers()-1
        futures[i] = @spawnat id partial_pi(batch .+ (i-1)*(length(r)÷nworkers()))
    end
    p = sum(fetch.(futures))
end
p - pi

# But that's rather annoying — needing to carefully divide up our workflow and
# manually collect all our results and such.  There's an easier way:

@time p = @distributed (+) for r in [(0:9999) .+ offset for offset in 0:10000:r[end]-1]
    partial_pi(r)
end
p - pi

# Why is this different from `@threads for` and `@simd for`? Why not just
# `@distributed for`?  Why the `@distributed (+) for`?

#-

# ## Data movement

#-

# Remember: Moving data is _expensive_!
#
# | System Event                   | Actual Latency | Scaled Latency |
# | ------------------------------ | -------------- | -------------- |
# | One CPU cycle                  |     0.4 ns     |     1 s        |
# | Level 1 cache access           |     0.9 ns     |     2 s        |
# | Level 2 cache access           |     2.8 ns     |     7 s        |
# | Level 3 cache access           |      28 ns     |     1 min      |
# | Main memory access (DDR DIMM)  |    ~100 ns     |     4 min      |
# | Intel Optane memory access     |     <10 μs     |     7 hrs      |
# | NVMe SSD I/O                   |     ~25 μs     |    17 hrs      |
# | SSD I/O                        |  50–150 μs     | 1.5–4 days     |
# | Rotational disk I/O            |    1–10 ms     |   1–9 months   |
# | Internet call: SF to NYC       |      65 ms     |     5 years    |
# | Internet call: SF to Hong Kong |     141 ms     |    11 years    |
#
# You really don't want to be taking a trip to the moon very frequently.
# Communication between processes can indeed be as expensive as hitting a disk —
# sometimes they're even implemented that way.
#
# So that's why Julia has special support for reductions built in to the
# `@distributed` macro: each worker can do its own (intermediate) reduction
# before returning just one value to our master node.

#-

# But sometimes you need to see those intermediate values. If you have a
# very expensive computation relative to the communication overhead, there are
# several ways to do this. The easiest is `pmap`:

@time pmap(partial_pi, [(0:99999) .+ offset for offset in 0:100000:r[end]-1])

# But if we have a large computation relative to the number of return values,
# pmap is great and easy.
#
# Increase the work on each worker by 100x and reduce the amount of communication by 100x:

@time pmap(partial_pi, [(0:9999999) .+ offset for offset in 0:10000000:r[end]-1])

# There are other ways of doing this, though, too — we'll get to them in a minute.
# But first, there's something else that I glossed over: the `@everywhere`s above.

#-

# ## Code movement

#-

# Each node is _completely_ independent; it's like starting brand new, separate
# Julia processes yourself. By default, `addprocs()` just launches the
# appropriate number of workers for the current workstation that you're on, but
# you can easily connect them to remote machines via SSH or even through cluster
# managers.

#-

# Those `@everywhere`s above are very important! They run the given expression
# on all workers to make sure the state between them is consistent. Without it,
# you'll see errors like this:

hello() = "hello world"
r = @spawnat 2 hello()

#-

fetch(r)

# Note that this applies to packages, too!

using Statistics # The Statistics stdlib defines mean
fetch(@spawnat 2 mean(rand(100_000)))

#-

@everywhere using Statistics
fetch(@spawnat 2 mean(rand(100_000)))

# # Other ways to structure and/or share data between processes
#
# Unlike `@threads`, we no longer have access to the same memory. While this
# does make expressing some algorithms a little more tricky, the "default"
# is much safer! There isn't any shared state to begin with, so it's harder
# to write an incorrect algorithm. It's also just harder to write some
# algorithms in the first place.
#
# So there are some special array types that can help bridge the gap between
# processes and make writing parallel code a bit easier.

#-

# ## The `SharedArray`
#
# If all workers are on the same physical machine, while they cannot share
# memory, they do all have shared access to the same hard drive(s)!
#
# The `SharedArray` makes use of this fact, allowing concurrent accesses to the
# same array — somewhat akin to threads default state.
#
# This is the prefix definition from the "thinking in parallel" course:
#
# ```
# using .Threads
# function prefix_threads!(y, ⊕)
#     l=length(y)
#     k=ceil(Int, log2(l))
#     for j=1:k
#         @threads for i=2^j:2^j:min(l, 2^k)       #"reduce"
#             @inbounds y[i] = y[i-2^(j-1)] ⊕ y[i]
#         end
#     end
#     for j=(k-1):-1:1
#         @threads for i=3*2^(j-1):2^j:min(l, 2^k) #"expand"
#             @inbounds y[i] = y[i-2^(j-1)] ⊕ y[i]
#         end
#     end
#     y
# end
# ```

using SharedArrays
function prefix!(y::SharedArray, ⊕)
    l=length(y)
    k=ceil(Int, log2(l))
    for j=1:k
        @distributed for i=2^j:2^j:min(l, 2^k)       #"reduce"
            @inbounds y[i] = y[i-2^(j-1)] ⊕ y[i]
        end
    end
    for j=(k-1):-1:1
        @distributed for i=3*2^(j-1):2^j:min(l, 2^k) #"expand"
            @inbounds y[i] = y[i-2^(j-1)] ⊕ y[i]
        end
    end
    y
end
data = rand(1_000_000);
A = SharedArray(data);

#-

prefix!(SharedArray(data), +) # compile
@time prefix!(A, +);

#-

A ≈ cumsum(data)

# What went wrong?

function prefix!(y::SharedArray, ⊕)
    l=length(y)
    k=ceil(Int, log2(l))
    for j=1:k
        @sync @distributed for i=2^j:2^j:min(l, 2^k)       #"reduce"
            @inbounds y[i] = y[i-2^(j-1)] ⊕ y[i]
        end
    end
    for j=(k-1):-1:1
        @sync @distributed for i=3*2^(j-1):2^j:min(l, 2^k) #"expand"
            @inbounds y[i] = y[i-2^(j-1)] ⊕ y[i]
        end
    end
    y
end
A = SharedArray(data)
@time prefix!(A, +)

#-

A ≈ cumsum(data)

# ## DistributedArrays
#
# We can, though, turn the problem on its head and allow the _data_ itself
# to determine how the problem gets split up. This can save us tons of indexing
# headaches.

@everywhere using Distributed
@everywhere using DistributedArrays
A = DArray(I->fill(myid(), length.(I)), (24, 24))

# The first argument takes a function that transforms the given set of indices
# to the _local portion_ of the distributed array.

A = DArray((24,24)) do I
    @show I
    fill(myid(), length.(I))
end

# Notice that none of the array actually lives on processor 1, but we can still
# display the contents — when we do we're requesting all workers give us their
# current data! While we've only talked about master-worker communcation so far,
# workers can communicate directly amongst themselves, too (by default).

@everywhere using BenchmarkTools
fetch(@spawnat 2 @benchmark $A[1,1])

#-

fetch(@spawnat 2 @benchmark $A[end,end])

# So it's fastest to work on a `DArray`'s "local" portion, but it's _possible_
# to grab other data if need be. This is perfect for any sort of tiled operation
# that works on neighboring values (like image filtering/convolution). Or Conway's
# game of life!

function life_step(d::DArray)
    DArray(size(d),procs(d)) do I
        # Compute the indices of the outside edge (that will come from other processors)
        top   = mod1(first(I[1])-1,size(d,1))
        bot   = mod1( last(I[1])+1,size(d,1))
        left  = mod1(first(I[2])-1,size(d,2))
        right = mod1( last(I[2])+1,size(d,2))
        # Create a new, temporary array that holds the local part + outside edge
        old = Array{Bool}(undef, length(I[1])+2, length(I[2])+2)
        # These accesses will pull data from other processors
        old[1      , 1      ] = d[top , left]
        old[2:end-1, 1      ] = d[I[1], left]   # left side (and corners)
        old[end    , 1      ] = d[bot , left]
        old[1      , end    ] = d[top , right]
        old[2:end-1, end    ] = d[I[1], right]  # right side (and corners)
        old[end    , end    ] = d[bot , right]
        old[1      , 2:end-1] = d[top , I[2]]   # top
        old[end    , 2:end-1] = d[bot , I[2]]   # bottom
        # But this big one is all local!
        old[2:end-1, 2:end-1] = d[I[1], I[2]]   # middle

        life_rule(old) # Compute the new segment!
    end
end
@everywhere function life_rule(old)
    # Now this part — the computational part — is entirely local and on Arrays!
    m, n = size(old)
    new = similar(old, m-2, n-2)
    for j = 2:n-1
        @inbounds for i = 2:m-1
            nc = (+)(old[i-1,j-1], old[i-1,j], old[i-1,j+1],
                     old[i  ,j-1],             old[i  ,j+1],
                     old[i+1,j-1], old[i+1,j], old[i+1,j+1])
            new[i-1,j-1] = (nc == 3 || nc == 2 && old[i,j])
        end
    end
    new
end

#-

A = DArray(I->rand(Bool, length.(I)), (20,20))
using Colors
Gray.(A)

#-

B = copy(A)

#-

B = Gray.(life_step(B))

# ## Clusters and more ways to distribute
#
# You can easily connect to completely separate machines with SSH access built in!
# But there are many other ways to connect to clusters:
#
# * [JuliaRun](https://juliacomputing.com/products/juliarun)
# * [Kubernetes](https://juliacomputing.com/blog/2018/12/15/kuber.html)
# * [MPI](https://github.com/JuliaParallel/MPI.jl)
# * [Cluster job queues with ClusterManagers](https://github.com/JuliaParallel/ClusterManagers.jl)
# * [Hadoop](https://github.com/JuliaParallel/Elly.jl)
# * [Spark](https://github.com/dfdx/Spark.jl)

#-

# # Multi-process parallelism is the heavy-duty workhorse in Julia
#
# It can tackle very large problems and distribute across a very large number
# of workers. Key things to remember
#
# * Each worker is a completely independent Julia process
#     * Data must move to them
#     * Code must move to them
# * Structure your algorithms and use a distributed mechanism that fits with the
#   time and memory parameters of your problem
#     * `@distributed` can be good for reductions and even relatively fast inner loops with limited (or no) explicit data transfer
#     * `pmap` is great for very expensive inner loops that return a value
#     * `SharedArray`s can be an easier drop-in replacement for threading-like behaviors (on a single machine)
#     * `DistributedArray`s can turn the problem on its head and let the data do the work splitting!

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
struct Camera{T<:Real}
    origin::Point3{T}
    rotation::RotMatrix3{T}
    horizontal::Vec3{T}
    vertical::Vec3{T}
    lower_left_corner::Vec3{T}
    aspect_ratio::T
end

origin(cam::Camera) = cam.origin
rotation(cam::Camera) = cam.rotation
horizontal(cam::Camera) = cam.horizontal
vertical(cam::Camera) = cam.vertical
lower_left_corner(cam::Camera) = cam.lower_left_corner
aspect_ratio(cam::Camera) = cam.aspect_ratio

function generate_ray(camera::Camera{T}, s::T, t::T)::Ray{T} where {T<:Real}
    Ray(origin(camera), normalize(lower_left_corner(camera) + s*horizontal(camera) + t*vertical(camera) - origin(camera)))
end

function make_camera(lookfrom::Point3{T}, lookat::Point3{T}, vup::Vec3{T}, vfov::T, aspect_ratio::T)::Camera{T} where {T<:Real}
    theta = deg2rad(vfov)
    h = tan(theta/2)
    viewport_height = 2 * h
    viewport_width = aspect_ratio * viewport_height

    w = normalize(Vec3{T}(lookfrom - lookat))
    u = normalize(cross(vup, w))
    v = cross(w, u)

    origin = lookfrom
    rotation = RotMatrix3{T}([w;u;v])
    horizontal = viewport_width * u
    vertical = viewport_height * v
    lower_left_corner = Vec3{T}(origin - horizontal/2 - vertical/2 - w)

    camera = Camera{T}(origin, rotation, horizontal, vertical, lower_left_corner, aspect_ratio)
end


text:
using .PyCall
using Unitful, UnitfulAtomic

export convert_from_ase_atoms
export convert_to_ase_atoms

const ase = PyNULL()
copy!(ase, pyimport_conda("ase", "ase", "conda-forge"))

convert_to_ase_atoms(atoms::Atoms, R::Matrix) =
    ase.Atoms(positions=ustrip.(u"Å", R'u"bohr"), symbols=string.(atoms.types))

convert_to_ase_atoms(atoms::Atoms, R::Matrix, ::InfiniteCell) =
    convert_to_ase_atoms(atoms, R)

function convert_to_ase_atoms(atoms::Atoms, R::Matrix, cell::PeriodicCell)
    ase.Atoms(
        positions=ustrip.(u"Å", R'u"bohr"),
        cell=ustrip.(u"Å", cell.vectors'u"bohr"),
        symbols=string.(atoms.types),
        pbc=cell.periodicity)
end

function convert_to_ase_atoms(atoms::Atoms, R::Vector{<:Matrix}, cell::AbstractCell)
    convert_to_ase_atoms.(Ref(atoms), R, Ref(cell))
end

convert_from_ase_atoms(ase_atoms::PyObject) =
    Atoms(ase_atoms), positions(ase_atoms), Cell(ase_atoms)

Atoms(ase_atoms::PyObject) = Atoms{Float64}(Symbol.(ase_atoms.get_chemical_symbols()))

positions(ase_atoms::PyObject) = austrip.(ase_atoms.get_positions()'u"Å")

function Cell(ase_atoms::PyObject)
    if all(ase_atoms.cell.array .== 0)
        return InfiniteCell()
    else
        return PeriodicCell{Float64}(austrip.(ase_atoms.cell.array'u"Å"), ase_atoms.pbc)
    end
end


text:
module DifferentiableBackwardEuler

import DifferentialEquations
import ChainRulesCore
import LinearAlgebra
import NLsolve

#solve M * (y1 - y0) - deltah * f(y1, p, t) = 0
function dae_step(y0, deltah, f, f_y, f_p, f_t, p, t, M; kwargs...)
	function f!(F, y)
		myF = M * (y .- y0) .- deltah * f(y, p, t)
		copy!(F, myF)
	end
	function j!(J, y)
		thisJ = M - deltah * f_y(y, p, t)
		copyto!(J, thisJ)
	end
	J0 = f_y(y0, p, t)
	j!(J0, y0)
	F0 = f(y0, p, t)
	f!(F0, y0)
	df = NLsolve.OnceDifferentiable(f!, j!, y0, F0, J0)
	soln = NLsolve.nlsolve(df, y0; kwargs...)
	if !NLsolve.converged(soln)
		@show soln
		@show t
		error("solution did not converge")
	end
	return soln.zero

end

#solve y1 - y0 - deltah * f(y1, p, t) = 0
function step(y0, deltah, f, f_y, f_p, f_t, p, t; kwargs...)
	return dae_step(y0, deltah, f, f_y, f_p, f_t, p, t, LinearAlgebra.I; kwargs...)
end

function make_dae_step_pullback(y1, y0, deltah, f, f_y, f_p, f_t, p, t, M)
	function dae_step_pullback(delta)
		step_y0 = -M
		step_deltah = -f(y1, p, t)
		step_p = -deltah * f_p(y1, p, t)
		step_t = -deltah * f_t(y1, p, t)
		lambda = transpose(M - deltah * f_y(y1, p, t)) \ delta
		return (ChainRulesCore.NoTangent(),#step function
				@ChainRulesCore.thunk(-transpose(step_y0) * lambda),#y0
				@ChainRulesCore.thunk(-transpose(step_deltah) * lambda),#deltah
				ChainRulesCore.NoTangent(), ChainRulesCore.NoTangent(), ChainRulesCore.NoTangent(), ChainRulesCore.NoTangent(),#f through f_t
				@ChainRulesCore.thunk(-transpose(step_p) * lambda),#p
				@ChainRulesCore.thunk(-transpose(step_t) * lambda),#t
				ChainRulesCore.NoTangent())#M -- this is wrong but we're not supporting differentation w.r.t. M
	end
	return dae_step_pullback
end

function ChainRulesCore.rrule(::typeof(dae_step), y0, deltah, f, f_y, f_p, f_t, p, t, M; kwargs...)
	y1 = dae_step(y0, deltah, f, f_y, f_p, f_t, p, t, M; kwargs...)
	pullback = make_dae_step_pullback(y1, y0, deltah, f, f_y, f_p, f_t, p, t, M)
	return y1, pullback
end

function steps(y0, f, f_y, f_p, f_t, p, ts; kwargs...)
	return dae_steps(y0, f, f_y, f_p, f_t, p, ts, LinearAlgebra.I; kwargs...)
end

function dae_steps(y0, f, f_y, f_p, f_t, p, ts, M; kwargs...)
	ys = zeros(eltype(y0), length(y0), length(ts))
	ys[:, 1] = y0
	for i = 1:length(ts) - 1
		ys[:, i + 1] = dae_step(ys[:, i], ts[i + 1] - ts[i], f, f_y, f_p, f_t, p, ts[i + 1], M; kwargs...)
	end
	return ys
end

function make_dae_steps_pullback(y, y0, f, f_y, f_p, f_t, p, ts, M)
	dae_step_pullbacks = [make_dae_step_pullback(y[:, i], y[:, i - 1], ts[i] - ts[i - 1], f, f_y, f_p, f_t, p, ts[i], M) for i = 2:size(y, 2)]
	function dae_steps_pullback(delta)
		_, thisdy0, _, _, _, _, _, thisdp, _ = dae_step_pullbacks[size(y, 2) - 1](delta[:, size(y, 2)])
		dy0 = ChainRulesCore.unthunk(thisdy0) .+ delta[:, size(y, 2) - 1]
		dp = ChainRulesCore.unthunk(thisdp)
		for i = reverse(1:size(y, 2) - 2)
			_, thisdy0, _, _, _, _, _, thisdp, _ = dae_step_pullbacks[i](dy0)
			dy0 .= ChainRulesCore.unthunk(thisdy0) .+ delta[:, i]
			dp .= ChainRulesCore.unthunk(thisdp) .+ dp
		end
		return (ChainRulesCore.NoTangent(),#steps
				dy0,#y0
				ChainRulesCore.NoTangent(), ChainRulesCore.NoTangent(), ChainRulesCore.NoTangent(), ChainRulesCore.NoTangent(),#f through f_t
				dp,#p
				ChainRulesCore.ZeroTangent(),#ts -- just make it zero for now
				ChainRulesCore.NoTangent())#M -- again this is wrong, but we basically don't support this for now
	end
	return dae_steps_pullback
end

function ChainRulesCore.rrule(::typeof(dae_steps), y0, f, f_y, f_p, f_t, p, ts, M; kwargs...)
	y = dae_steps(y0, f, f_y, f_p, f_t, p, ts, M; kwargs...)
	dae_steps_pullback = make_dae_steps_pullback(y, y0, f, f_y, f_p, f_t, p, ts, M)
	return y, dae_steps_pullback
end

function steps_diffeq(y0, f, f_y, f_p, f_t, p, t0, tfinal; kwargs...)
	return dae_steps_diffeq(y0, f, f_y, f_p, f_t, p, t0, tfinal, LinearAlgebra.I; kwargs...)
end
function dae_steps_diffeq(y0, f, f_y, f_p, f_t, p, t0, tfinal, M; kwargs...)
	odef = DifferentialEquations.ODEFunction(f; jac=f_y, jac_prototype=f_y(y0, p, t0), mass_matrix=M)
	prob = DifferentialEquations.ODEProblem(odef, y0, (t0, tfinal), p)
	soln = DifferentialEquations.solve(prob, DifferentialEquations.ImplicitEuler(); kwargs...)
	return soln
end
function steps(y0, f, f_y, f_p, f_t, p, t0, tfinal; kwargs...)
	return dae_steps(y0, f, f_y, f_p, f_t, p, t0, tfinal, LinearAlgebra.I; kwargs...)
end
function dae_steps(y0, f, f_y, f_p, f_t, p, t0, tfinal, M; kwargs...)
	return dae_steps_diffeq(y0, f, f_y, f_p, f_t, p, t0, tfinal, M; kwargs...)[:, :]
end

function make_dae_steps_pullback(y, y0, f, f_y, f_p, f_t, p, t0, tfinal, ts, M)
	pullback1 = make_dae_steps_pullback(y, y0, f, f_y, f_p, f_t, p, ts, M)
	function dae_steps_pullback(delta)
		_, dy0, _, _, _, _, dp, _ = pullback1(delta)
		return (ChainRulesCore.NoTangent(),#steps
				dy0,#y0
				ChainRulesCore.NoTangent(), ChainRulesCore.NoTangent(), ChainRulesCore.NoTangent(), ChainRulesCore.NoTangent(),#f through f_t
				dp,#p
				ChainRulesCore.ZeroTangent(),#t0 -- just make it zero for now
				ChainRulesCore.ZeroTangent(),#tfinal -- just make it zero for now
				ChainRulesCore.NoTangent())#M -- this is incorrect but differentiating w.r.t. M is no supported
	end
	return dae_steps_pullback
end

function ChainRulesCore.rrule(::typeof(dae_steps), y0, f, f_y, f_p, f_t, p, t0, tfinal, M; kwargs...)
	soln = dae_steps_diffeq(y0, f, f_y, f_p, f_t, p, t0, tfinal, M; kwargs...)
	y = soln[:, :]
	steps_pullback = make_dae_steps_pullback(y, y0, f, f_y, f_p, f_t, p, t0, tfinal, soln.t, M)
	return y, steps_pullback
end

end


text:
using ProbNumDiffEq
using Test
using LinearAlgebra

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems;
importodeproblems();
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_fitzhughnagumo, prob_ode_vanstiff

@testset "Smoothing with small constant steps" begin
    prob = ProbNumDiffEq.remake_prob_with_jac(prob_ode_fitzhughnagumo)
    @test solve_ieks(prob, IEKS(order=4, diffusionmodel=:fixed)) isa
          ProbNumDiffEq.ProbODESolution
end


text:
module FluidFluidQueues

using DiscretisedFluidQueues
using LinearAlgebra
using Random
using SparseArrays

include("queue.jl")
include("operators.jl")
include("partition.jl")
include("simulate.jl")

export Rates, FluidFluidQueue
export augment_model
export AbstractSign, Plus, Minus, Zero, NotPlus, NotMinus, NotZero
export index
export FluidFluidOperator, FluidFluidGenerator, RatesOperator, InOutGenerator
export build_psi, build_xi, build_limit_dist_operators
export simulate, first_exit_y

end


text:
module Proxy
include("handler_proxy.jl")
include("client_protocol_proxy.jl")
include("writer_proxy.jl")
end


text:
export convolve
export marchonintime

function convolve(Z,x,j,k0)
    T = promote_type(eltype(Z), eltype(x))
    M,N,K = size(Z)
    @assert M == size(x,1)
y = zeros(T,M)
    for m in 1:M
        for n in 1:N
            for k in k0:min(j,K)
                i = j - k + 1
                y[m] += Z[m,n,k] * x[n,i]
            end
        end
    end
    return y
end

"""
    marchonintime(Z0,Z,B,I)

Solve by marching-on-in-time the causal convolution problem defined by `(W0,Z,B)`
up to timestep `I`. Here, `Z` is an array of order 3 that contains a discretisation
of a time translation invariant retarded potential operator. `W0` is the inverse of
the slice `Z[:,:,1]`.
"""
function marchonintime(W0,Z,B,I)
    T = eltype(Z)
    M,N,K = size(Z)
    @assert M == size(B,1)
    x = zeros(T,N,I)
    for i in 1:I
        b = B[:,i] - convolve(Z,x,i,2)
        x[:,i] = W0 * b
        (i % 10 == 0) && print(i, "[", I, "] - ")
    end
    return x
end

function companion(Z)
    M, N, K = size(Z)
end


text:
using DrWatson
@quickactivate "PECUZAL_Julia"

using DynamicalSystems
using DelayEmbeddings
using DelimitedFiles
using BenchmarkTools

include("../../src/pecuzal_method.jl")
include("../../src/data_analysis_functions.jl")

## We analyze the computational complexity of the proposed PECUZAL method in
# comparison to TDE, the method from G&A and MDOP. We feed in the y-component
# of the Rössler system for different time series length and compute the average
# computation time.

roe = Systems.roessler()

# set time interval for integration
N = 20000 # number of samples
dt_tra = 0.05
t = 0:dt_tra:(dt_tra*N)

tr = trajectory(roe, (N*dt_tra); dt = dt_tra, Ttr = (2000*dt_tra))
tr = regularize(tr)

# define function which wrap the whole embedding process
function tde_time(data, w, taus)
    # use cao's method
    Y_tde, τ_tde, _ = optimal_traditional_de(data; w = w, τs = taus)
end

function GA_time(data, w, taus)
    Y_GA, τ_vals_GA, ts_vals_GA, FNNs_GA , ns_GA = garcia_almeida_embedding(data;
                                                        τs = taus , w = w, T = w)
end

function mdop_time(data, w)
    tw1 = mdop_maximum_delay(data)
    lm1 = DelayEmbeddings.findlocalminima(tw1[2])
    if lm1[1] ≤ 2
        lmm1 = try lm1[2]
        catch
            lm1[1]
        end
    else
        lmm1 = lm1[1]
    end

    taus_mdop = 0:lmm1

    Y_mdop, τ_vals_mdop, ts_vals_mdop, FNNs_mdop , βs_mdop = mdop_embedding(data;
                                                        τs = taus_mdop , w = w)
end

function pec_time(data, w, taus)
    Y_pec, τ_vals_pec, ts_vals_pec, Ls_pec , εs_pec = pecuzal_embedding_update(data;
                                                                τs = taus , w = w)
end

step = 1000
Nss = 1000:step:N

times_tde = zeros(length(Nss))
times_GA = zeros(length(Nss))
times_mdop = zeros(length(Nss))
times_pec = zeros(length(Nss))

for (i,Ns) in enumerate(Nss)
    display("run: $i")
    data = tr[1:Ns,2]
    # compute Theiler window
    w = estimate_delay(data, "mi_min", 0:150)
    # and range of encountered time delays
    taus1 = 1:(4*w)
    taus2 = 0:(4*w)

    b_tde = @benchmark tde_time($data, $w, $taus1)
    b_GA = @benchmark GA_time($data, $w, $taus2)
    b_mdop = @benchmark mdop_time($data, $w)
    b_pec = @benchmark pec_time($data, $w, $taus2)

    times_tde[i] = median(b_tde).time
    times_GA[i] = median(b_GA).time
    times_mdop[i] = median(b_mdop).time
    times_pec[i] = median(b_pec).time

end

writedlm("./scripts/performance/results/times_tde.csv", times_tde)
writedlm("./scripts/performance/results/times_GA.csv", times_GA)
writedlm("./scripts/performance/results/times_mdop.csv", times_mdop)
writedlm("./scripts/performance/results/times_pec.csv", times_pec)
writedlm("./scripts/performance/results/times.csv", Nss)


text:
function build_list_of_approx_tensors_test(triangles::Array{Tuple{Tuple{Int,Int,Int},Tuple{F,F,F}},1},
    eps::F=1.0;avg_bins=false,test_mode=false) where {F <: AbstractFloat}
    """-------------------------------------------------------------------------
        Builds a dictionary linking approximate angles to that each triangle
      closest maps to and returns a list of pairs containing the indices and the
      approximated triangle angles for a given threshold eps.

    Inputs:
    -------
    * triangles - (List of Int 3-tuple, Float64 3-tuple, pairs ):
      Dictionary linking the indices of the hyperedges to the angles in the
      triangle.
    * eps - (optional float):
      The resolution for approximating angles (angle bin size).
      TODO: change var to bin_size
    * avg_bins - (optional bool):
      Indicates whether or not to alter the approximated angles to the average
      of each point within the bin, as opposed to the mid-point.

    Outputs:
    --------
    * Tensors - (Dict of 3-tuple -> list of 3-tuples):
      Dictionary linking approximate angles to the list of indices associated
      with the triangles approximated by each approximate angles.
    -------------------------------------------------------------------------"""
    #could probably improve performance using 3 tuples
    Tensors = Dict{Array{F,1},Array{Tuple{Int,Int,Int},1}}()

    if avg_bins
        new_approx_angle_tensors = Dict{Array{F,1},Array{Tuple{Int,Int,Int},1}}()
        bin_averages = Dict{Array{F,1},Array{F,1}}()
    end

    for (indices,og_angles) in triangles
        for (perm_indices,angles) in zip(permutations(indices),permutations(og_angles))

            aprox_angle1 = eps * (floor(angles[1] / eps) + .5)
            aprox_angle2 = eps * (floor(angles[2] / eps) + .5)
            aprox_angle3 = 180 - aprox_angle1 - aprox_angle2
            approx_angles = [aprox_angle1,aprox_angle2,aprox_angle3]

            indices = tuple(perm_indices[1],perm_indices[2],perm_indices[3])

            if avg_bins
                if haskey(bin_averages,approx_angles)
                    bin_averages[approx_angles] += [angles[1],angles[2],angles[3]]
                else
                    bin_averages[approx_angles] = [angles[1],angles[2],angles[3]]
                end
            end

            if haskey(Tensors,approx_angles)
                push!(Tensors[approx_angles],indices)
            else
                Tensors[approx_angles] = [indices]
            end
        end
    end

    if avg_bins

        for (approx_angles,indices) in Tensors
            avg = bin_averages[approx_angles]/length(indices)
            new_approx_angle_tensors[avg] = indices
        end

        return new_approx_angle_tensors
    end

    return Tensors

end

function Make_HOFASM_tensor_pairs_test(H_indices::Array{Array{Int64,2},1},B_indices::Array{Array{Int64,2},1},
                            B_vals::Array{Array{Float64,1},1},n::Int,m::Int)
    marginalized_ten_pairs = Array{Tuple{SparseMatrixCSC{Float64,Int64},SparseMatrixCSC{Float64,Int64}},1}(undef,0)
    for (H_idx,B_idx,v) in zip(H_indices,B_indices,B_vals)
        push!(marginalized_ten_pairs,
              (perm_mode1_marginalization(H_idx,ones(size(H_idx,1)),[1,2,3],n),
               perm_mode1_marginalization(B_idx,v,[1,2,3],m))
             )
    end
    return marginalized_ten_pairs
end

function HOFASM_contraction!_test(tensor_pairs::Array{NTuple{2,SparseMatrixCSC{Float64,Int64}},1},
                             X_k::Array{Float64,2},X_k_1::Array{Float64,2}) #where {F <: AbstractFloat}

    n,m = size(X_k_1)

    for i=1:n
       for j=1:m
           X_k_1[i,j] = 0.0
       end
    end

    for (Hn,Bn) in tensor_pairs
         X_k_1 .+= (Bn*(Hn*X_k)' )'
    end

end


text:
# allele_freq.jl
# ==============
#
# Compute allele frequencies with BioJulia data types.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/GeneticVariation.jl/blob/master/LICENSE.md

"""
    gene_frequencies(seqcounts::Composition{T}) where T <: Sequence

Compute gene_frequencies from a `BioSequences.Composition` variable that contains
unique sequence counts.
"""
function gene_frequencies(seqcounts::Composition{T}) where T <: Sequence
    n = sum(values(seqcounts))
    frequencies = Dict{T, Float64}()
    @inbounds for (seq, count) in seqcounts
        frequencies[seq] = count / n
    end
    return frequencies
end

"""
    gene_frequencies(iterable)

Compute the gene frequencies for any iterable with an `eltype` which is a
concrete subtype of the abstract `Sequence` type.
"""
function gene_frequencies(iterable)
    return _gene_frequencies(iterable, eltype(iterable), Base.IteratorSize(iterable))
end

# Default for most iterables, throws an error.
_gene_frequencies(iterable, eltype, is) = error("Iterable not supported.")

# Action to take for a sequence iterable which has no known size.
function _gene_frequencies(iterable, ::Type{<:Sequence}, is::Base.SizeUnknown)
    composition = BioSequences.composition(iterable)
    return gene_frequencies(composition)
end

# Action to take for a sequence iterable which has a known size.
# This version computes frequencies directly as n is known in advance.
# The other method first computes composition, and computes frequencies from
# that.
function _gene_frequencies(sequences, ::Type{T}, is::Union{Base.HasLength,Base.HasShape}) where T<:Sequence
    inc = 1 / length(sequences)
    frequencies = Dict{T, Float64}()
    @inbounds for seq in sequences
        old = get(frequencies, seq, 0.0)
        frequencies[seq] = old + inc
    end
    return frequencies
end


text:
using Stheno
using TemporalGPs: to_sde, smooth, SArrayStorage
using Optim
using Zygote: gradient
using Random

using GPARatScale
export create_lgssm, get_sde_predictions

"""
Create LGSSM for the latent locations and the given observation noise.
The user can specify a noise vector, in which case it overrides the noise from
the parameters.
"""
function create_lgssm(
    latent_locations,
    l,
    process_var,
    noise_sigma,
    kernel_structure::Kernel;
    sde_storage::SArrayStorage = SArrayStorage(Float64),
    noise_vector = nothing,
    debug::Bool = false,
)
    # Transform the gp into a  Linear Gaussian State Space model (LGSSM)
    # This is done by indexing into the LTISDE.
    if debug
        println("Creating LGSSM with parameters\n\tl=$(l)\n\tprocess_var=$(process_var)\n\tnoise_sigma=$(noise_sigma)")
    end
    # Generate the GP
    gp_kernel = kernel(kernel_structure, l = l, s = process_var^2)
    gp = GP(gp_kernel, GPC())
    # Transform the gp into an LTISDE
    gp_sde = to_sde(gp, sde_storage)
    if isnothing(noise_vector)
        lgssm = gp_sde(latent_locations, noise_sigma^2)
        return lgssm
    else
        lgssm = gp_sde(latent_locations, noise_vector)
        return lgssm
    end
end

"""
Create an LGSSM conditioned on the observed data and output the marginals of
the observations at the output locations
"""
function get_sde_predictions(
    data_locations,
    data_outputs,
    output_locations;
    kernel_structure::Kernel = Matern52(),
    sde_storage::SArrayStorage = SArrayStorage(Float64),
    debug::Bool = true,
)
    # Concatanete all the input locations
    latent_locations = vcat(data_locations, output_locations)
    # Concatenate the knwon outputs with the outputs we want to infer
    outputs = vcat(data_outputs, repeat([0], length(output_locations)))
    # Get the permutation for sorting the latent vars since TemporalGPs expects
    # a timeseries. Use the permutation for later sorting
    sorting_perm = sortperm(latent_locations)
    reverse_perm = sortperm(sorting_perm)  # used to reverse the sorting
    # Sort the inputs and output so that the input locations are in
    # ascending order. This is needed in order to model a timeseries.
    s_latent_locations = latent_locations[sorting_perm]
    s_outputs = outputs[sorting_perm]

    # Helper function used to compute negative log marignal lik for params
    function nlml(params)
        l, process_var, noise_sigma = unpack_gp(params)
        lgssm = create_lgssm(
            data_locations,
            l,
            process_var,
            noise_sigma,
            kernel_structure,
        )
        return -logpdf(lgssm, data_outputs)
    end
    # Optimize the parameters of the GP
    params = randn(3)
    results = Optim.optimize(nlml, params, NelderMead())
    opt_l, opt_process_var, opt_noise_sigma = unpack_gp(results.minimizer)
    # Create the noise vector passed into the creation of the LGSSM.
    # We assume (almost) infinite noise for the output locations
    noise_vector = vcat(
        repeat([opt_noise_sigma^2], length(data_locations)),
        repeat([1e10], length(output_locations)),
    )
    s_noise_vector = noise_vector[sorting_perm]

    opt_lgssm = create_lgssm(
        s_latent_locations,
        opt_l,
        opt_process_var,
        opt_noise_sigma,
        kernel_structure,
        noise_vector = s_noise_vector,
        debug = debug,
    )
    # Get the marginal posteriors for all observations
    _, y_smooth, _ = smooth(opt_lgssm, s_outputs)

    output_observations =
        y_smooth[reverse_perm][length(data_outputs)+1:length(y_smooth)]
    return opt_lgssm, output_observations
end

max_stars_repo_name: willtebbutt/GPAR-at-scale

max_issues_repo_name: willtebbutt/GPAR-at-scale

max_forks_repo_name: willtebbutt/GPAR-at-scale


text:
import Base.remotecall
import Base.remotecall_fetch
import Base.remotecall_wait

const WS_MODE_RR = 1    # Round Robin
const WS_MODE_FF = 2    # First Free

export WS_MODE_RR, WS_MODE_FF
export WorkerSet

type WorkerSet
    idset::Array{Int,1}
    freeset::Array{Int,1}
    nextidx::Int
    mode::Int
    req_wait::RemoteRef
    worker_wait::RemoteRef

    WorkerSet(w, mode) = begin
        ws = new(w, w, 1, mode, RemoteRef(), RemoteRef())
        if (mode == WS_MODE_FF)
            # start a new task that will manage the request Q
            @schedule begin
                manage_q(ws)
            end
        end
        ws
    end
end

WorkerSet(w) = WorkerSet(w, WS_MODE_RR)

function remotecall(ws::WorkerSet, f, args...)
    if (ws.mode != WS_MODE_RR)
        error("remotecall(WorkerSet) only supported in WS_MODE_RR mode")
    end

    remotecall(chooseproc(ws), f, args...)
end

remotecall_fetch(ws::WorkerSet, f, args...) = remotecall_internal(ws, :FETCH, f, args...)
remotecall_wait(ws::WorkerSet, f, args...) = remotecall_internal(ws, :WAIT, f, args...)

function remotecall_internal(ws::WorkerSet, rrtype, f, args...)
    if (ws.mode == WS_MODE_RR)
        if (rrtype == :FETCH)
            return remotecall_fetch(chooseproc(ws), f, args...)
        else
            return remotecall_wait(chooseproc(ws), f, args...)
        end
    elseif (ws.mode == WS_MODE_FF)
        if length(ws.freeset) > 0
            id = shift!(ws.freeset)
            if rrtype == :FETCH
                resp = remotecall_fetch(id, f, args...)
            else
                resp = remotecall_wait(id, f, args...)
            end
            push_and_notify(ws, id)
            return resp
        else
            rr = RemoteRef()
            put(ws.req_wait, (rr, () -> f(args...), rrtype))
            return take(rr)
        end
    end
    error("Unsupported WorkerSet mode $(ws.mode)")
end

function chooseproc(ws::WorkerSet)
    nextidx = ws.nextidx
    if nextidx > length(ws.idset)
        p = ws.idset[1]
        ws.nextidx = 2
    else
        p = ws.idset[nextidx]
        ws.nextidx += 1
    end
    p
end

function manage_q(ws::WorkerSet)
#    println("Started manage_q")
    while (true)

        # Tuple of remoteref, remotefunction and type.
        (rr, c, t) = take(ws.req_wait)
#        println("Got a Queued request")
        if length(ws.freeset) == 0
#            println("Waiting for a free worker")
            take(ws.worker_wait)
        end
        if length(ws.freeset) == 0
            println("FIXME: No Free Worker found")
            continue
        end

        id = shift!(ws.freeset)
        @async let id=id,c=c,rr=rr
            if t == :FETCH
                resp = remotecall_fetch(id, c)
            else
                resp = remotecall_wait(id, c)
            end
            push_and_notify(ws, id)
            put(rr, resp)
        end
    end
end

function push_and_notify(ws, id)
    push!(ws.freeset, id)
    if (!isready(ws.worker_wait))
        put(ws.worker_wait, true)
    end
end


text:
module StreamStats
    import StatsBase, Distributions

    export update!, state, nobs, replicates

    include("streamstat.jl")
    include("mean.jl")
    include("var.jl")
    include("moments.jl")
    include("cov.jl")
    include("min.jl")
    include("max.jl")
    # include("approx_quantile.jl")
    # include("approx_cdf.jl")
    # include("approx_pmf.jl")
    include("hyper_log_log.jl")
    include("approx_distinct.jl")
    include("bootstrap.jl")
    include("sample.jl")
    include("approx_ols.jl")
    include("approx_logit.jl")
    include("approx_ridge.jl")
    include("approx_l2_logit.jl")
end


text:
@testset "primitive" begin

    arrtype = KnetLayers.arrtype

    @testset "Multiply" begin
        inputDim = 10; outputDim = 2;
        m = Multiply(input=inputDim, output=outputDim, winit=randn)
        x = arrtype(zeros(inputDim, 1))
        y1 = m(x)
        y2 = arrtype(zeros(outputDim, 1))
        @test y1==y2
    end

    @testset "Embed" begin
        m  = Embed(input=10, output=2, winit=randn)
        x  = arrtype(zeros(10,1))
        ind = rand(1:10)
        x[ind,1] = 1.0
        y1 = m(x)
        y2 = m([ind])
        @test y1==y2
    end

    @testset "Bias" begin
        m  = Bias(10)
        @test size(m.b) == (10,)
        @test sum(m.b) == 0
     end

    @testset "Linear" begin
        inputDim=10; outputDim=3
        m = Linear(input=inputDim, output=outputDim, winit=randn,binit=zeros)
        x = arrtype(zeros(10,2))
        y = m(x)
        @test true
    end

    @testset "Dense" begin
        inputDim=10; outputDim=3
        m = Dense(input=inputDim, output=outputDim, winit=randn,binit=zeros)
        x = arrtype(zeros(10,2))
        y = m(x)
        @test true
    end
end


text:
module BreaskCancerExample
using IntervalCensoredEventTimes
using CSV, DataFrames, Ipopt, NLopt, SparseArrays, Random, Test, UnicodePlots

@testset "Breast Cancer" begin
bcfile = joinpath(@__DIR__, "..", "data", "breast_cancer.csv")
df  = CSV.File(bcfile) |> DataFrame
g   = Float64.(df[!, :Group] .== "RCT")
L   = Float64.(df[!, :L])
R   = map(x -> ismissing(x) ? Inf : Float64(x), df[!, :R])
icm = IntervalCensoredModel(reshape(g, length(g), 1), L, R)
display([g icm.L icm.R icm.C]); println()
@info "NPMLE of survival functin S(t)"
npmle!(icm)
plt = lineplot([0; icm.ts], [1; icm.S₀], xlabel = "Time by Months", ylabel = "Survival Functions")
display(plt); println()
@info "Log-likelihood at starting point"
initialize_uniform!(icm)
display(DataFrame(t=icm.ts, St=icm.S₀, Λ₀t=icm.Λ₀, λ₀t=icm.λ₀)); println()
println("logl = $(loglikelihood!(icm))")
@info "MLE of PH model"
# solver = NLopt.NLoptSolver(algorithm = :LN_BOBYQA, ftol_rel = 1e-12, ftol_abs = 1e-8, maxeval = 10000)
solver = Ipopt.IpoptSolver(print_level=5)
# solver = NLopt.NLoptSolver(algorithm = :LN_BOBYQA)
# solver = NLopt.NLoptSolver(algorithm = :LN_COBYLA, maxeval = 10000)
# solver = NLopt.NLoptSolver(algorithm = :LD_MMA)
# solver = NLopt.NLoptSolver(algorithm = :LD_SLSQP)
# solver = NLopt.NLoptSolver(algorithm = :LD_LBFGS)
fit!(icm, solver, init = initialize_uniform!(icm))
display(icm)
@info "IHT"
Random.seed!(123)
X = randn(length(icm.C), 100)
k = 5
icsm = IntervalCensoredSparseModel(icm, X, k)
iht!(icsm, verbose=true, debiasing=true, printfreq=1)
βx = sparsevec(icsm.βx)
display(DataFrame(j=βx.nzind, βj=βx.nzval)); println()
end

end


text:
"""
`Kernel` is a module implementing filtering (correlation) kernels of full
dimensionality. The following kernels are supported:

  - `sobel`
  - `prewitt`
  - `ando3`, `ando4`, and `ando5`
  - `scharr`
  - `bickley`
  - `gaussian`
  - `DoG` (Difference-of-Gaussian)
  - `LoG` (Laplacian-of-Gaussian)
  - `Laplacian`
  - `gabor`
  - `moffat`

See also: [`KernelFactors`](@ref).
"""
module Kernel

using StaticArrays, OffsetArrays
using ..ImageFiltering
using ..ImageFiltering.KernelFactors
import ..ImageFiltering: _reshape, IdentityUnitRange

# We would like to do `using ..ImageFiltering.imgradients` so that that
# Documenter.jl (the documentation system) can parse a reference such as `See
# also: [`ImageFiltering.imgradients`](@ref)`. However, imgradients is not yet
# in scope because of the order in which include files are included into
# ImageFiltering.jl. With the more general `using ImageFiltering`, we seem to
# sidestep the scope problem, although I don't actually understand the mechanism
# form why this works. - ZS

function product2d(kf)
    k1, k2 = kf
    k1[1].*k1[2], k2[1].*k2[2]
end

"""
```julia
    diff1, diff2 = sobel()
```

Return ``3 \\times 3`` correlation kernels for two-dimensional gradient compution
using the Sobel operator. The `diff1` kernel computes the gradient along the
y-axis (first dimension), and the `diff2` kernel computes the gradient along the
x-axis (second dimension). `diff1 == rotr90(diff2)`

```julia
    (diff,) = sobel(extended::NTuple{N,Bool}, d)
```
Return (a tuple of) the N-dimensional correlation kernel for gradient compution
along the dimension `d` using the Sobel operator. If `extended[dim]` is false,
`diff` will have size 1 along that dimension.

# Citation
P.-E. Danielsson and O. Seger, "Generalized and separable sobel operators," in  *Machine Vision for Three-Dimensional Scenes*,  H. Freeman, Ed.  Academic Press, 1990,  pp. 347–379. [doi:10.1016/b978-0-12-266722-0.50016-6](https://doi.org/doi:10.1016/b978-0-12-266722-0.50016-6)

See also: [`KernelFactors.sobel`](@ref), [`Kernel.prewitt`](@ref),
[`Kernel.ando3`](@ref), [`Kernel.scharr`](@ref), [`Kernel.bickley`](@ref) and
[`imgradients`](@ref).
"""
sobel() = product2d(KernelFactors.sobel())

sobel(extended, d) = (broadcast(*, KernelFactors.sobel(extended, d)...),)

"""
```julia
    diff1, diff2 = prewitt()
```

Return ``3 \\times 3`` correlation kernels for two-dimensional gradient compution
using the Prewitt operator. The `diff1` kernel computes the gradient along the
y-axis (first dimension), and the `diff2` kernel computes the gradient along the
x-axis (second dimension). `diff1 == rotr90(diff2)`

```julia
    (diff,) = prewitt(extended::NTuple{N,Bool}, d)
```
Return (a tuple of) the N-dimensional correlation kernel for gradient compution
along the dimension `d` using the Prewitt operator. If `extended[dim]` is false,
`diff` will have size 1 along that dimension.

# Citation
J. M. Prewitt, "Object enhancement and extraction," *Picture processing and Psychopictorics*, vol. 10, no. 1, pp. 15–19, 1970.

See also: [`KernelFactors.prewitt`](@ref), [`Kernel.sobel`](@ref),
[`Kernel.ando3`](@ref), [`Kernel.scharr`](@ref),[`Kernel.bickley`](@ref) and
[`ImageFiltering.imgradients`](@ref).
"""
prewitt() = product2d(KernelFactors.prewitt())

prewitt(extended, d) = (broadcast(*, KernelFactors.prewitt(extended, d)...),)

"""
```julia
    diff1, diff2 = ando3()
```

Return ``3 \\times 3`` correlation kernels for two-dimensional gradient
compution using Ando's "optimal" filters. The `diff1` kernel computes the
gradient along the y-axis (first dimension), and the `diff2` kernel computes the
gradient along the x-axis (second dimension). `diff1 == rotr90(diff2)`

```julia
    (diff,) = ando3(extended::NTuple{N,Bool}, d)
```
Return (a tuple of) the N-dimensional correlation kernel for gradient compution
along the dimension `d` using Ando's "optimal" filters of size 3. If
`extended[dim]` is false, `diff` will have size 1 along that dimension.

# Citation
S. Ando, "Consistent gradient operators," *IEEE Transactions on Pattern Analysis and Machine Intelligence*, vol. 22, no.3, pp. 252–265, 2000. [doi:10.1109/34.841757](https://doi.org/doi:10.1109/34.841757)

See also: [`KernelFactors.ando3`](@ref), [`Kernel.ando4`](@ref),
[`Kernel.ando5`](@ref) and  [`ImageFiltering.imgradients`](@ref).
"""
ando3() = product2d(KernelFactors.ando3())

ando3(extended, d) = (broadcast(*, KernelFactors.ando3(extended, d)...),)

"""
```julia
    diff1, diff2 = ando4()
```

Return ``4 \\times 4`` correlation  kernels for two-dimensional gradient
compution using Ando's "optimal" filters.  The `diff1` kernel computes the
gradient along the y-axis (first dimension), and  the `diff2` kernel computes
the gradient along the x-axis (second dimension). `diff1 == rotr90(diff2)`

```julia
    (diff,) = ando4(extended::NTuple{N,Bool}, d)
```
Return (a tuple of) the N-dimensional correlation kernel for gradient compution
along the dimension `d` using Ando's "optimal" filters of size 4. If
`extended[dim]` is false, `diff` will have size 1 along that dimension.

# Citation
S. Ando, "Consistent gradient operators," *IEEE Transactions on Pattern Analysis and Machine Intelligence*, vol. 22, no.3, pp. 252–265, 2000. [doi:10.1109/34.841757](https://doi.org/doi:10.1109/34.841757)

See also: [`KernelFactors.ando4`](@ref), [`Kernel.ando3`](@ref),
[`Kernel.ando5`](@ref) and [`ImageFiltering.imgradients`](@ref).
"""
function ando4()
    f = centered(@SMatrix [ -0.022116 -0.025526  0.025526  0.022116
                            -0.098381 -0.112984  0.112984  0.098381
                            -0.098381 -0.112984  0.112984  0.098381
                            -0.022116 -0.025526  0.025526  0.022116 ])
    return f', f
end

function ando4(extended::Tuple{Bool,Bool}, d)
    all(extended) || error("all dimensions must be extended")
    (ando4()[d],)
end

"""
```julia
    diff1, diff2 = ando5()
```

Return ``5 \\times 5`` correlation  kernels for two-dimensional gradient
compution using Ando's "optimal" filters.  The `diff1` kernel computes the
gradient along the y-axis (first dimension), and  the `diff2` kernel computes
the gradient along the x-axis (second dimension). `diff1 == rotr90(diff2)`

```julia
    (diff,) = ando5(extended::NTuple{N,Bool}, d)
```
Return (a tuple of) the N-dimensional correlation kernel for gradient compution
along the dimension `d` using Ando's "optimal" filters of size 5. If
`extended[dim]` is false, `diff` will have size 1 along that dimension.

# Citation
S. Ando, "Consistent gradient operators," *IEEE Transactions on Pattern Analysis and Machine Intelligence*, vol. 22, no.3, pp. 252–265, 2000. [doi:10.1109/34.841757](https://doi.org/doi:10.1109/34.841757)

See also: [`KernelFactors.ando5`](@ref), [`Kernel.ando3`](@ref),
[`Kernel.ando4`](@ref) and  [`ImageFiltering.imgradients`](@ref).
"""
function ando5()
    f = centered(@SMatrix [ -0.003776 -0.010199  0.0  0.010199  0.003776
                            -0.026786 -0.070844  0.0  0.070844  0.026786
                            -0.046548 -0.122572  0.0  0.122572  0.046548
                            -0.026786 -0.070844  0.0  0.070844  0.026786
                            -0.003776 -0.010199  0.0  0.010199  0.003776 ])
    return f', f
end

function ando5(extended::Tuple{Bool,Bool}, d)
    all(extended) || error("all dimensions must be extended")
    (ando5()[d],)
end

"""
```julia
    diff1, diff2 = scharr()
```

Return ``3 \\times 3`` correlation kernels for two-dimensional gradient
compution using the Scharr operator. The `diff1` kernel computes the gradient
along the y-axis (first dimension), and the `diff2` kernel  computes the
gradient along the x-axis (second dimension). `diff1 == rotr90(diff2)`

```julia
    (diff,) = scharr(extended::NTuple{N,Bool}, d)
```
Return (a tuple of) the N-dimensional correlation kernel for gradient compution
along the dimension `d` using the Scharr operator. If `extended[dim]` is false,
`diff` will have size 1 along that dimension.

# Citation
H. Scharr and  J. Weickert, "An anisotropic diffusion algorithm with optimized rotation invariance," *Mustererkennung 2000*, pp. 460–467, 2000. [doi:10.1007/978-3-642-59802-9_58](https://doi.org/doi:10.1007/978-3-642-59802-9_58)

See also: [`KernelFactors.scharr`](@ref), [`Kernel.prewitt`](@ref),
[`Kernel.ando3`](@ref), [`Kernel.bickley`](@ref) and
[`ImageFiltering.imgradients`](@ref).
"""
scharr() = product2d(KernelFactors.scharr())

scharr(extended, d) = (broadcast(*, KernelFactors.scharr(extended, d)...),)

"""
```julia
    diff1, diff2 = bickley()
```

Return ``3 \\times 3`` correlation kernels for two-dimensional gradient
compution using the Bickley operator. The `diff1` kernel computes the gradient
along the y-axis (first dimension), and the `diff2` kernel computes the gradient
along the x-axis (second dimension). `diff1 == rotr90(diff2)`

```julia
    (diff,) = bickley(extended::NTuple{N,Bool}, d)
```
Return (a tuple of) the N-dimensional correlation kernel for gradient compution
along the dimension `d` using the Bickley operator. If `extended[dim]` is false,
`diff` will have size 1 along that dimension.

# Citation
W. G. Bickley, "Finite difference formulae for the square lattice," *The Quarterly Journal of Mechanics and Applied Mathematics*, vol. 1, no. 1, pp. 35–42, 1948.  [doi:10.1093/qjmam/1.1.35](https://doi.org/doi:10.1137/12087092x)

See also: [`KernelFactors.bickley`](@ref), [`Kernel.prewitt`](@ref),
[`Kernel.ando3`](@ref),  [`Kernel.scharr`](@ref) and
[`ImageFiltering.imgradients`](@ref).
"""
bickley() = product2d(KernelFactors.bickley())

bickley(extended, d) = (broadcast(*, KernelFactors.bickley(extended, d)...),)

"""
    gaussian((σ1, σ2, ...), [(l1, l2, ...)]) -> g
    gaussian(σ)                  -> g

Construct a multidimensional gaussian filter, with standard deviation
`σd` along dimension `d`. Optionally provide the kernel length `l`,
which must be a tuple of the same length.

If `σ` is supplied as a single number, a symmetric 2d kernel is
constructed.

See also: [`KernelFactors.gaussian`](@ref).
"""
@inline gaussian(σs::NTuple{N,Real}, ls::NTuple{N,Integer}) where {N} =
    broadcast(*, KernelFactors.gaussian(σs, ls)...)
gaussian(σ::Tuple{Real}, l::Tuple{Integer}) = KernelFactors.gaussian(σ[1], l[1])
gaussian(σ::Tuple{}, l::Tuple{}) = reshape([1])  # 0d
gaussian(σs::AbstractVector{T}, ls::AbstractVector{I}) where {T<:Real,I<:Integer} =
    gaussian((σs...,), (ls...,))

@inline gaussian(σs::NTuple{N,Real}) where {N} = broadcast(*, KernelFactors.gaussian(σs)...)
gaussian(σs::AbstractVector{T}) where {T<:Real} = gaussian((σs...,))
gaussian(σ::Tuple{Real}) = KernelFactors.gaussian(σ[1])
gaussian(σ::Tuple{}) = reshape([1])

gaussian(σ::Real) = gaussian((σ, σ))

"""
    DoG((σp1, σp2, ...), (σm1, σm2, ...), [l1, l2, ...]) -> k
    DoG((σ1, σ2, ...))                                   -> k
    DoG(σ::Real)                                         -> k

Construct a multidimensional difference-of-gaussian kernel `k`, equal
to `gaussian(σp, l)-gaussian(σm, l)`.  When only a single `σ` is
supplied, the default is to choose `σp = σ, σm = √2 σ`. Optionally
provide the kernel length `l`; the default is to extend by two
`max(σp,σm)` in each direction from the center. `l` must be odd.

If `σ` is provided as a single number, a symmetric 2d DoG kernel is
returned.

See also: [`KernelFactors.IIRGaussian`](@ref).
"""
DoG(σps::NTuple{N,Real}, σms::NTuple{N,Real}, ls::NTuple{N,Integer}) where {N} =
    gaussian(σps, ls) - gaussian(σms, ls)
function DoG(σps::NTuple{N,Real}) where N
    σms = map(s->s*√2, σps)
    neg = gaussian(σms)
    l = map(length, axes(neg))
    gaussian(σps, l) - neg
end
DoG(σ::Real) = DoG((σ,σ))

"""
    LoG((σ1, σ2, ...)) -> k
    LoG(σ)             -> k

Construct a Laplacian-of-Gaussian kernel `k`. `σd` is the gaussian width
along dimension `d`.  If `σ` is supplied as a single number, a
symmetric 2d kernel is returned.

See also: [`KernelFactors.IIRGaussian`](@ref) and [`Kernel.Laplacian`](@ref).
"""
function LoG(σs::NTuple{N}) where N
    ws = map(n->(ceil(Int,8.5*n)>>1), σs)
    R = CartesianIndices(map(w->IdentityUnitRange(-w:w), ws))
    σ = SVector(σs)
    C = 1/(prod(σ)*(2π)^(N/2))
    σ2 = σ.^2
    σ2i = sum(1 ./ σ2)
    function df(I::CartesianIndex, σ2, σ2i)
        x = SVector(Tuple(I))
        xσ = x.^2 ./ σ2
        (sum(xσ./σ2) - σ2i) * exp(-sum(xσ)/2)
    end
    [C*df(I, σ2, σ2i) for I in R]
end
LoG(σ::Real) = LoG((σ,σ))

struct Laplacian{N}
    flags::NTuple{N,Bool}
    offsets::Vector{CartesianIndex{N}}

    function Laplacian{N}(flags::NTuple{N,Bool}) where {N}
        offsets = Array{CartesianIndex{N}}(undef, 0)
        for i = 1:N
            if flags[i]
                push!(offsets,
                      CartesianIndex{N}((ntuple(d->0, i-1)..., 1, ntuple(d->0, N-i)...)))
            end
        end
        new{N}(flags, offsets)
    end
end

"""
    Laplacian((true,true,false,...))
    Laplacian(dims, N)
    Laplacian()

Laplacian kernel in `N` dimensions, taking derivatives along the
directions marked as `true` in the supplied tuple. Alternatively, one
can pass `dims`, a listing of the dimensions for
differentiation. (However, this variant is not inferrable.)

`Laplacian()` is the 2d laplacian, equivalent to `Laplacian((true,true))`.

The kernel is represented as an opaque type, but you can use
`convert(AbstractArray, L)` to convert it into array format.
"""
Laplacian(flags::NTuple{N,Bool}) where {N} = Laplacian{N}(flags)
Laplacian() = Laplacian((true,true))

function Laplacian(dims, N::Int)
    flags = falses(N)
    flags[[dims...]] .= true
    Laplacian((flags...,))
end

Base.axes(L::Laplacian) = map(f->f ? (-1:1) : (0:0), L.flags)
Base.isempty(L::Laplacian) = false
function Base.convert(::Type{AbstractArray}, L::Laplacian{N}) where N
    A = fill!(OffsetArray{Int}(undef, axes(L)), 0)
    for I in L.offsets
        A[I] = A[-I] = 1
    end
    A[ntuple(d->0, Val(N))...] = -2*length(L.offsets)
    A
end
_reshape(L::Laplacian{N}, ::Val{N}) where {N} = L

"""
    gabor(size_x,size_y,σ,θ,λ,γ,ψ) -> (k_real,k_complex)

Returns a 2 Dimensional Complex Gabor kernel contained in a tuple where

  - `size_x`, `size_y` denote the size of the kernel
  - `σ` denotes the standard deviation of the Gaussian envelope
  - `θ` represents the orientation of the normal to the parallel stripes of a Gabor function
  - `λ` represents the wavelength of the sinusoidal factor
  - `γ` is the spatial aspect ratio, and specifies the ellipticity of the support of the Gabor function
  - `ψ` is the phase offset

#Citation
N. Petkov and P. Kruizinga, “Computational models of visual neurons specialised in the detection of periodic and aperiodic oriented visual stimuli: bar and grating cells,” Biological Cybernetics, vol. 76, no. 2, pp. 83–96, Feb. 1997. doi.org/10.1007/s004220050323
"""
function gabor(size_x::Integer, size_y::Integer, σ::Real, θ::Real, λ::Real, γ::Real, ψ::Real)

    σx = σ
    σy = σ/γ
    nstds = 3
    c = cos(θ)
    s = sin(θ)

    validate_gabor(σ,λ,γ)

    if(size_x > 0)
        xmax = floor(Int64,size_x/2)
    else
        @warn "The input parameter size_x should be positive. Using size_x = 6 * σx + 1 (Default value)"
        xmax = round(Int64,max(abs(nstds*σx*c),abs(nstds*σy*s),1))
    end

    if(size_y > 0)
        ymax = floor(Int64,size_y/2)
    else
        @warn "The input parameter size_y should be positive. Using size_y = 6 * σy + 1 (Default value)"
        ymax = round(Int64,max(abs(nstds*σx*s),abs(nstds*σy*c),1))
    end

    xmin = -xmax
    ymin = -ymax

    x = [j for i in xmin:xmax,j in ymin:ymax]
    y = [i for i in xmin:xmax,j in ymin:ymax]
    xr = x*c + y*s
    yr = -x*s + y*c

    kernel_real = (exp.(-0.5*(((xr.*xr)/σx^2) + ((yr.*yr)/σy^2))).*cos.(2*(π/λ)*xr .+ ψ))
    kernel_imag = (exp.(-0.5*(((xr.*xr)/σx^2) + ((yr.*yr)/σy^2))).*sin.(2*(π/λ)*xr .+ ψ))

    kernel = (kernel_real,kernel_imag)
    return kernel
end

function validate_gabor(σ::Real,λ::Real,γ::Real)
    if !(σ>0 && λ>0 && γ>0)
        throw(ArgumentError("The parameters σ, λ and γ must be positive numbers."))
    end
end

"""
    moffat(α, β, ls) -> k

Constructs a 2D, symmetric Moffat kernel `k` with core width, `α`, and power, `β`.
Size of kernel defaults to 4 * full-width-half-max or as specified in `ls`.
See [this notebook](https://nbviewer.jupyter.org/github/ysbach/AO_2017/blob/master/04_Ground_Based_Concept.ipynb#1.2.-Moffat) for details.

# Citation
Moffat, A. F. J. "A theoretical investigation of focal stellar images in the photographic emulsion and application to photographic photometry." Astronomy and Astrophysics 3 (1969): 455.
"""
function moffat(α::Real, β::Real, ls::Tuple{Integer, Integer})
    ws = map(n->(ceil(Int,n)>>1), ls)
    R = CartesianIndices(map(w->IdentityUnitRange(-w:w), ws))
    α2 = α^2
    amp = (β - 1)/(π * α2)
    @. amp*((1+df(R)/α2)^-β)
end
moffat(α::Real, β::Real, ls::Integer)    = moffat(α, β, (ls,ls))
moffat(α::Real, β::Real)                 = moffat(α, β, ceil(Int, (α*2*sqrt(2^(1/β) - 1))*4))

@inline function df(I::CartesianIndex)
    x = SVector(Tuple(I))
    sum(x.^2)
end

"""
    reflect(kernel) --> reflectedkernel

Compute the pointwise reflection around 0, 0, ... of the kernel
`kernel`.  Using `imfilter` with a `reflectedkernel` performs convolution,
rather than correlation, with respect to the original `kernel`.
"""
function reflect(kernel::AbstractArray)
    inds = map(reflectind, axes(kernel))
    out = similar(kernel, inds)
    for I in CartesianIndices(axes(kernel))
        out[-I] = kernel[I]
    end
    out
end

reflectind(r::AbstractUnitRange) = -last(r):-first(r)

end


text:
struct SingleVariate <: MCTrait end

mutable struct MonteCarloModel{P <: AbstractPathPricer, RSG <: AbstractRandomSequenceGenerator, S <: StochasticProcess1D} <: AbstractMonteCarloModel
  pathGenerator::PathGenerator{RSG, S}
  pathPricer::P
  sampleAccumulator::RiskStatistics
  isAntitheticVariate::Bool
end

function add_samples!(mcmodel::MonteCarloModel, samples::Int, idx::Int=1)
  # re-init the risk data
  adding_data!(mcmodel.sampleAccumulator, samples)
  @inbounds @simd for j = 1:samples
    path = get_next!(mcmodel.pathGenerator)
    price = mcmodel.pathPricer(path.value)

    # TODO Control Variate
    if mcmodel.isAntitheticVariate
      path2 = get_antithetic!(mcmodel.pathGenerator)
      price2 = mcmodel.pathPricer(path2.value)
      add_sample!(mcmodel.sampleAccumulator, (price+price2)/2.0, path.weight, idx)
    else
      add_sample!(mcmodel.sampleAccumulator, price, path.weight, idx)
    end
    idx += 1
  end

  return mcmodel
end


text:
using Turing
import Turing.translate_tilde!

model_info = Dict(:name => "model", :main_body_names => Dict(:model => :model, :vi => :vi, :sampler => :sampler), :arg_syms => [], :tent_pvars_list => [])

ex = :(y ~ Normal(1,1))
model_info[:main_body] = ex
translate_tilde!(model_info)
res = model_info[:main_body]
Base.@assert res.head == :block

ex = quote
  x = 1
  y = rand()
  y ~ Normal(0,1)
end
model_info[:main_body] = ex
translate_tilde!(model_info)
res = model_info[:main_body]
Base.@assert res.head == :block


text:
function setup_unconstrained()
    options = Dict(:verbose => false,
                   :eps_abs => 1e-08,
                   :eps_rel => 1e-08,
                   :eps_dual_inf => 1e-18)
    return options
end

tol = 1e-5

@testset "unconstrained" begin

    @testset "unconstrained_problem" begin
        srand(1)

        n = 30
        m = 0
        P = sparse(Diagonal(rand(n)) + 0.2 * sparse(I, n, n))
        q = randn(n)
        A = spzeros(m, n)
        u = Float64[]
        l = Float64[]

        # Explicit solution
        invP = inv(Array(P))
        x_test = - invP * q
        y_test = zeros(m)
        obj_test = - .5 * q' * invP * q

        options = setup_unconstrained()

        model = OSQP.Model()
        OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)
        results = OSQP.solve!(model)

        @test isapprox(results.x, x_test, atol=tol)
        @test isapprox(results.y, y_test, atol=tol)
        @test isapprox(results.info.obj_val, obj_test, atol=tol)
        @test results.info.status == :Solved

    end

end