text:
"""
    LaplaceOutput
Return type for Laplace computations.
Fields are `nothing` on return unless requested.
See individual FMM/direct computation function
documentation for specifics.
"""
mutable struct LaplaceOutput <: FMMVals
    pot
    grad
    hess

    pottarg
    gradtarg
    hesstarg

    pgt
    pg
end
function LaplaceOutput(pgt,pg=nothing)
    LaplaceOutput(nothing,nothing,nothing,
                  nothing,nothing,nothing,
                  pgt,pg)
end

"""
```julia
    vals = lfmm3d(eps,sources;charges=nothing,dipvecs=nothing,
                  targets=nothing,pg=0,pgt=0,nd=1)
```
This function computes the N-body Laplace interactions
in three dimensions where the interaction kernel is given
by ``\$1/r\$`` and its gradients.  This is the
``O(N)`` fast multipole code which computes the interactions
to the requested precision.
```math
u(x) = \\sum_{j=1}^{N} c_{j} / \\|x-x_{j}\\| + v_{j}
\\cdot \\nabla( 1/\\|x-x_{j}\\|)  \\, ,
u(x) = \\sum_{j=1}^{N} c_{j} \\frac{1}{\\|x-x_{j}\\|}
- v_{j} \\cdot \\nabla \\left( \\frac{1}{\\|x-x_{j}\\|}
\\right)  \\, ,
```

where ``c_{j}`` are the charge densities,
      ``v_{j}`` are the dipole orientation vectors, and
      ``x_{j}`` are the source locations.
      When ``x=x_{m}``, the term corresponding to
      ``x_{m}`` is dropped from the sum
# Input
# Input
* `eps::Float64` precision requested
* `sources::Array{Float64}` size (3,n) source locations (``x_{j}``)
* `charges::Array{Float64}` size (nd,n) or (n) charge densities (c_{j})
* `dipvecs::Array{Float64}` size (nd,3,n) or (3,n)) dipole orientation vectors (``v_{j}``)
* `targets::Array{Float64}` size (3,nt) target locations (``x``)
* `pg::Integer` source eval flag.
    + Potential (``u``) at sources evaluated if `pg == 1`.
    + Potential and gradient (``\\nabla u``) at sources evaluated if `pg == 2`
* `pgt::Integer` target eval flag.
    + Potential at targets evaluated if `pgt == 1`.
    + Potential and gradient at targets evaluated if `pgt == 2`
* `nd::Integer` number of densities
# Output

`vals::LaplaceOutput` with the fields
Note that the Hessian is returned as a length
6 vector at each point with the second derivatives
in the order: ``\\partial_{xx}``, ``\\partial_{yy}``,
``\\partial_{zz}``, ``\\partial_{xy}``,
``\\partial_{xz}``, ``\\partial_{yz}``.
* `vals.pot::Array{Float64}` size (nd,n) potential at source locations if requested
* `vals.grad::Array{Float64}` size (nd,3,n) gradient at source locations if requested
* `vals.hess::Array{Float64}` size (nd,6,n) Hessian at source locations if requested
* `vals.pottarg::Array{Float64}` size (nd,nt) potential at target locations if requested
* `vals.gradtarg::Array{Float64}` size (nd,3,nt) gradient at target locations if requested
* `vals.hesstarg::Array{Float64}` size (nd,6,nt) Hessian at target locations if requested
"""
function lfmm3d(eps::Float64,sources::Array{Float64};charges::TFN=nothing,
                dipvecs::TFN=nothing,targets::TFN=nothing,
                pg::Integer=0,pgt::Integer=0,nd::Integer=1)

    # check inputs

    @assert size(sources,1) == 3
    @assert nd >= 0
    @assert (0 <= pg  && pg  <= 3) "flag pg not in expected range"
    @assert (0 <= pgt && pgt <= 3) "flag pgt not in expected range"

    # default values

    vals = LaplaceOutput(pgt,pg)

    ifcharge = 0
    ifdipole = 0

    zero = Float64(0)

    pot = zero
    grad = zero
    hess = zero
    pottarg = zero
    gradtarg = zero
    hesstarg = zero

    n = div(length(sources),3)

    nt = 0

    if targets !== nothing
        @assert size(targets,1) == 3
        nt = div(length(targets),3)
    else
        targets = 0.0
    end

    if charges !== nothing
        @assert div(length(charges),nd) == n
        ifcharge = 1
    else
        charges = zero
        ifcharge = 0
    end

    if dipvecs !== nothing
        @assert div(length(dipvecs),nd) == n*3
        ifdipole = 1
    else
        dipvecs = zero
        ifdipole = 0
    end

    if ifcharge == 0 && ifdipole == 0
        @warn "no charges or dipoles provided, doing nothing"
        return vals
    end

    if ((pg != 1 && pg != 2 && pg !=3) &&
        (pgt != 1 && pgt != 2 && pgt != 3))
        @warn "no output requested, doing nothing"
        return vals
    end

    if (pgt == 1 || pgt == 2 || pgt == 3) && targets === nothing
        @warn "target values requested but no targets provided"
    end

    # allocate memory for return values

    if pg == 1 || pg == 2 || pg == 3
        if nd > 1
            pot = zeros(Float64,nd,n)
        else
            pot = zeros(Float64,n)
        end
    end

    if pg == 2 || pg == 3
        if nd > 1
            grad = zeros(Float64,nd,3,n)
        else
            grad = zeros(Float64,3,n)
        end
    end

    if pg == 3
        if nd > 1
            hess = zeros(Float64,nd,6,n)
        else
            hess = zeros(Float64,6,n)
        end
    end

    if pgt == 1 || pgt == 2 || pgt == 3
        if nd > 1
            pottarg = zeros(Float64,nd,nt)
        else
            pottarg = zeros(Float64,nt)
        end
    end

    if pgt == 2 || pgt == 3
        if nd > 1
            gradtarg = zeros(Float64,nd,3,nt)
        else
            gradtarg = zeros(Float64,3,nt)
        end
    end

    if pgt == 3
        if nd > 1
            hesstarg = zeros(Float64,nd,6,nt)
        else
            hesstarg = zeros(Float64,6,nt)
        end
    end

    # actually call the function
    #
    # fortran calling sequence:
    # subroutine lfmm3d(nd,eps,nsource,source,ifcharge,
    #     $    charge,ifdipole,dipvec,ifpgh,pot,grad,hess,ntarg,
    #     $    targ,ifpghtarg,pottarg,gradtarg,hesstarg)
    #

    ccall((:lfmm3d_,libfmm3d),Cvoid,(Fi,Fd,Fi,Fd,Fi,Fd,
                                    Fi,Fd,Fi,Fd,Fd,Fd,Fi,
                                    Fd,Fi,Fd,Fd,Fd),
          nd,eps,n,sources,ifcharge,charges,ifdipole,
          dipvecs,pg,pot,grad,hess,nt,targets,pgt,
          pottarg,gradtarg,hesstarg)

    # load requested values

    if pg == 1 || pg == 2 || pg == 3; vals.pot = pot end
    if pg == 2 || pg == 3; vals.grad = grad end
    if pg == 3; vals.hess = hess end
    if pgt == 1 || pgt == 2 || pgt == 3; vals.pottarg = pottarg end
    if pgt == 2 || pgt == 3; vals.gradtarg = gradtarg end
    if pgt == 3; vals.hesstarg = hesstarg end

    return vals

end

"""
```julia
    vals = l3ddir(sources,targets;charges=nothing,
                dipvecs=nothing,pgt=0,nd=1,
                thresh=1e-16)
```
This function computes the N-body Laplace interactions
in three dimensions where the interaction kernel is given
by ``1/r`` and its gradients.  This is the
``O(N^2)`` direct evaluation code.
```math
 u(x) = \\sum_{j=1}^{N} c_{j} \\frac{1}{\\|x-x_{j}\\|}
- v_{j} \\cdot \\nabla \\left( \\frac{1}{\\|x-x_{j}\\|}
\\right)  \\, ,
```

where ``c_{j}`` are the charge densities,
      ``v_{j}`` are the dipole orientation vectors, and
      ``x_{j}`` are the source locations.
      When ``x=x_{m}``, the term corresponding to
      ``x_{m}`` is dropped from the sum
# Input
* `sources::Array{Float64}` size (3,n) source locations (``x_{j}``)
* `targets::Array{Float64}` size (3,nt) target locations (``x``)
* `charges::Array{ComplexF64}` size (nd,n) or (n) charge densities (c_{j})
* `dipvecs::Array{ComplexF64}` size (nd,3,n) or (3,n)) dipole orientation vectors (``v_{j}``)
* `pgt::Integer` target eval flag.
    + Potential at targets evaluated if `pgt == 1`.
    + Potential and gradient at targets evaluated if `pgt == 2`
* `nd::Integer` number of densities
* `thresh::Float64` threshold for ignoring interactions when ``\\|x-x_{j}\\| \\leq thresh``
# Output

`vals::LaplaceOutput` with the fields
* `vals.pottarg::Array{ComplexF64}` size (nd,nt) or (nt) potential at target locations if requested
* `vals.gradtarg::Array{ComplexF64}` size (nd,3,nt) or (3,nt) gradient at target locations if requested
"""
function l3ddir(sources::Array{Float64},targets::Array{Float64};
                charges::TFN=nothing,dipvecs::TFN=nothing,
                pgt::Integer=0,nd::Integer=1,
                thresh::Float64=1e-16)

    # check inputs

    @assert size(sources,1) == 3
    @assert size(targets,1) == 3
    @assert nd >= 0
    @assert (0 <= pgt && pgt <= 2) "flag pgt not in expected range"

    # default values

    vals = LaplaceOutput(pgt)

    ifcharge = 0
    ifdipole = 0

    zero = Float64(0)

    pottarg = zero
    gradtarg = zero
    hesstarg = zero

    n = div(length(sources),3)
    nt = div(length(targets),3)

    if charges !== nothing
        @assert div(length(charges),nd) == n
        ifcharge = 1
    else
        charges = zero
        ifcharge = 0
    end

    if dipvecs !== nothing
        @assert div(length(dipvecs),nd) == n*3
        ifdipole = 1
    else
        dipvecs = zero
        ifdipole = 0
    end

    if ifcharge == 0 && ifdipole == 0
        @warn "no charges or dipoles provided, doing nothing"
        return vals
    end

    if (pgt != 1 && pgt != 2)
        @warn "no output requested, doing nothing"
        return vals
    end

    # allocate memory for return values

    if pgt == 1 || pgt == 2
        if nd > 1
            pottarg = zeros(Float64,nd,nt)
        else
            pottarg = zeros(Float64,nt)
        end
    end

    if pgt == 2
        if nd > 1
            gradtarg = zeros(Float64,nd,3,nt)
        else
            gradtarg = zeros(Float64,3,nt)
        end
    end

    # dispatch to appropriate wrapper

    if pgt == 1
        if ifcharge == 1
            if ifdipole == 1
                l3ddirectcdp!(nd,sources,charges,
                              dipvecs,n,targets,nt,
                              pottarg,thresh)
                vals.pottarg = pottarg
            else
                l3ddirectcp!(nd,sources,charges,
                             n,targets,nt,pottarg,
                             thresh)
                vals.pottarg = pottarg
            end
        else
            if ifdipole == 1
                l3ddirectdp!(nd,sources,
                             dipvecs,n,targets,nt,
                             pottarg,thresh)
                vals.pottarg = pottarg
            end
        end
    elseif pgt == 2
        if ifcharge == 1
            if ifdipole == 1
                l3ddirectcdg!(nd,sources,charges,
                              dipvecs,n,targets,nt,
                              pottarg,gradtarg,thresh)
                vals.pottarg = pottarg
                vals.gradtarg = gradtarg
            else
                l3ddirectcg!(nd,sources,charges,
                             n,targets,nt,pottarg,
                             gradtarg,thresh)
                vals.pottarg = pottarg
                vals.gradtarg = gradtarg
            end
        else
            if ifdipole == 1
                l3ddirectdg!(nd,sources,
                              dipvecs,n,targets,nt,
                             pottarg,gradtarg,thresh)
                vals.pottarg = pottarg
                vals.gradtarg = gradtarg
            end
        end
    end

    return vals

end

function l3ddirectcp!(nd,sources,charges,n,targets,nt,
                      pot,thresh)
    ccall((:l3ddirectcp_,libfmm3d),Cvoid,(Fi,Fd,Fd,Fi,
                                          Fd,Fi,Fd,Fd),
          nd,sources,charges,n,targets,nt,pot,thresh)
    return
end

function l3ddirectdp!(nd,sources,dipvecs,n,targets,nt,
                      pot,thresh)
    ccall((:l3ddirectdp_,libfmm3d),Cvoid,(Fi,Fd,Fd,Fi,
                                          Fd,Fi,Fd,Fd),
          nd,sources,dipvecs,n,targets,nt,pot,thresh)
    return
end

function l3ddirectcdp!(nd,sources,charges,dipvecs,n,targets,
                       nt,pot,thresh)
    ccall((:l3ddirectcdp_,libfmm3d),Cvoid,(Fi,Fd,Fd,Fd,Fi,
                                          Fd,Fi,Fd,Fd),
          nd,sources,charges,dipvecs,n,targets,nt,pot,thresh)
    return
end

function l3ddirectcg!(nd,sources,charges,n,targets,nt,
                      pot,grad,thresh)
    ccall((:l3ddirectcg_,libfmm3d),Cvoid,(Fi,Fd,Fd,Fi,
                                          Fd,Fi,Fd,Fd,Fd),
          nd,sources,charges,n,targets,nt,pot,grad,thresh)
    return
end

function l3ddirectdg!(nd,sources,dipvecs,n,targets,nt,
                      pot,grad,thresh)
    ccall((:l3ddirectdg_,libfmm3d),Cvoid,(Fi,Fd,Fd,Fi,
                                          Fd,Fi,Fd,Fd,Fd),
          nd,sources,dipvecs,n,targets,nt,pot,grad,thresh)
    return
end

function l3ddirectcdg!(nd,sources,charges,dipvecs,n,targets,
                       nt,pot,grad,thresh)
    ccall((:l3ddirectcdg_,libfmm3d),Cvoid,(Fi,Fd,Fd,Fd,Fi,
                                          Fd,Fi,Fd,Fd,Fd),
          nd,sources,charges,dipvecs,n,targets,nt,pot,grad,thresh)
    return
end


text:
# using Revise
using ADCME
using ADSeismic
using PyPlot
using DelimitedFiles
include("database.jl")

scale = 201

if length(ARGS)==1
    scale = parse(Int64, ARGS[1])
end

num_repeat = 3

param = AcousticPropagatorParams(NX=scale, NY=scale,
    NSTEP=100, DELTAT=1e-4,  DELTAX=1.0, DELTAY=1.0,
    PropagatorKernel = 1)

rc = Gauss(param, 1/pi, 100, 1e6)
srci = [div(param.NX,2)]
srcj = [div(param.NY,2)]
srcv = reshape(rc, :, 1)
src = AcousticSource(srci, srcj, srcv)

C = placeholder(3300*ones(param.NX+2, param.NY+2))
model = AcousticPropagatorSolver(param, src, C)

loss = sum(model.u^2)
g = gradients(loss, C)
sess = Session(); init(sess)
u = run(sess, model.u)
_ = run(sess, g)

time = 0
for i = 1:num_repeat
    init(sess)
    global time += (@timed run(sess, model.u))[2]
end
time /= num_repeat
t1 = time

time = 0
for i = 1:num_repeat
    init(sess)
    global time += (@timed run(sess, g))[2]
end
time /= num_repeat
t2 = time

insert_record("gpu_simulation", scale, t1, t2)

ext: jl


text:
using Distributions
using SpinGlassNetworks

function ramp(t::T, τ::T, α::T, pi::T, pf::T) where T <: Real
    p = (pf + pi) + (pf - pi) * tanh(α * (2.0 * t / τ - 1.0))
    p / 2.0
end

@testset "Simple Coherent Ising Machine simulator for small Ising instances." begin
    L = 4
    ig = ising_graph("$(@__DIR__)/instances/basic/$(L)_001.txt")

    scale = 0.2
    noise = Normal(0.0, 0.1)

    x0 = zeros(L)
    sat = 1.0
    time = 1000.
    pi, pf, α = -15.0, 0.0, 2.0
    momentum = 0.4

    pump = [ramp(t, time, α, pi, pf) for t ∈ 0:time]

    opo = OpticalOscillators{Float64}(ig, scale, noise)
    dyn = OPODynamics{Float64}(x0, sat, pump, momentum)

    @testset  "OpticalOscillators and OPODynamics work properly." begin
        @test opo.scale ≈ scale
        @test opo.noise ≈ noise

        @test dyn.initial_state ≈ x0
        @test dyn.saturation ≈ sat
        @test dyn.pump ≈ pump
    end

    N = 500
    states = Vector{Vector{Int}}(undef, N)
    Threads.@threads for i ∈ 1:N
        states[i] = evolve_optical_oscillators(opo, dyn)
    end

    @test minimum(energy(states, ig)) ≈ brute_force(ig, :CPU, num_states=1).energies[1]
end


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

using Markdown
using InteractiveUtils

# ╔═╡ c6144a10-47a0-11eb-22b9-dfd100f82e1f
begin
	using Pkg
	Pkg.add.(["Images","ImageIO","ImageMagick"])
	using Images
end

# ╔═╡ 27ba2730-47a1-11eb-1b96-cf281ec5864c
url = "https://i.imgur.com/VGPeJ6s.jpg"

# ╔═╡ 6d7e90d0-47a1-11eb-23b0-b956b56c1b22
download(url,"philip.jpg")

# ╔═╡ 0b54a4b0-47a3-11eb-2fda-c12dded73b7b
RGBX(1,0.5,0.5)

# ╔═╡ 1998d59e-47a3-11eb-24c7-0df767077078
size(philip)

# ╔═╡ 2b246960-47a3-11eb-0f4d-534a1fd3d56b
philip[1:1000,1:400]

# ╔═╡ 6cb402f0-47a3-11eb-0fe2-1349ec0f89c9
begin
	(h, w) = size(philip)
	head = philip[(h ÷ 2):h, (w ÷ 10):(9w ÷ 10)]
end

# ╔═╡ 9ce7561e-47a3-11eb-303c-f14193f1052f
size(head)

# ╔═╡ a01ae2d0-47a3-11eb-361a-37a0b8256bc6
size(philip)

# ╔═╡ a2af2240-47a3-11eb-0de5-034785d18c92
[head head;head head]

# ╔═╡ c0170eae-47a3-11eb-3635-7900777c0748
[
	head reverse(head,dims=2)
	reverse(head, dims=1) reverse(reverse(head,dims=2),dims=1)
]

# ╔═╡ 0a659450-47a4-11eb-3db3-dbe84bf0939b
new_phil = copy(head)

# ╔═╡ 4dd8baf0-47a4-11eb-1683-85215ab7f5bc
new_phil[100:200,1:100] .= RGB(0,1,0)

# ╔═╡ 5f3a7f8e-47a4-11eb-207c-658d802f61be
new_phil

# ╔═╡ 9884f0f0-47a4-11eb-376e-ed8838083417
function redify(color)
	return RGB(color.r,0,0)
end

# ╔═╡ a4d408a0-47a4-11eb-22ec-ab85c6551f41
new_phil2 = copy(new_phil)

# ╔═╡ b2c97d52-47a4-11eb-0dfe-4115fcd1f59f
begin
	poor_phil = decimate(new_phil,3)
	size(poor_phil)
	poor_phil
end

# ╔═╡ 33ab4e20-47a6-11eb-07a8-8bb46671145c
convolve(head, blur(2))

# ╔═╡ a14b3332-47a8-11eb-14e2-196bfbfeecfe

# ╔═╡ a1199be0-47a8-11eb-1971-fdb89c561105
element = 1//1

# ╔═╡ a578b040-47a8-11eb-1fbd-f9e5eb680df1
fill(element,3,4)

# ╔═╡ cb7ff2d0-47a8-11eb-3bc9-f361ff13bbf0
typeof(element)

# ╔═╡ ebdcb53e-47a8-11eb-1660-353aa67bfa6e
keeptrack = [typeof(1), typeof(1.0), typeof("one"), typeof(1//1)]

# ╔═╡ fc2d09e0-47a8-11eb-22c2-013caff40755
typeof(keeptrack)

# ╔═╡ 03bf5cd0-47a9-11eb-3a9a-cfaf895e6cbf
string

# ╔═╡ Cell order:
# ╠═c6144a10-47a0-11eb-22b9-dfd100f82e1f
# ╠═27ba2730-47a1-11eb-1b96-cf281ec5864c
# ╠═6d7e90d0-47a1-11eb-23b0-b956b56c1b22
# ╠═0b54a4b0-47a3-11eb-2fda-c12dded73b7b
# ╠═1998d59e-47a3-11eb-24c7-0df767077078
# ╠═2b246960-47a3-11eb-0f4d-534a1fd3d56b
# ╠═6cb402f0-47a3-11eb-0fe2-1349ec0f89c9
# ╠═9ce7561e-47a3-11eb-303c-f14193f1052f
# ╠═a01ae2d0-47a3-11eb-361a-37a0b8256bc6
# ╠═a2af2240-47a3-11eb-0de5-034785d18c92
# ╠═c0170eae-47a3-11eb-3635-7900777c0748
# ╠═0a659450-47a4-11eb-3db3-dbe84bf0939b
# ╠═4dd8baf0-47a4-11eb-1683-85215ab7f5bc
# ╠═5f3a7f8e-47a4-11eb-207c-658d802f61be
# ╠═9884f0f0-47a4-11eb-376e-ed8838083417
# ╠═a4d408a0-47a4-11eb-22ec-ab85c6551f41
# ╠═b2c97d52-47a4-11eb-0dfe-4115fcd1f59f
# ╠═33ab4e20-47a6-11eb-07a8-8bb46671145c
# ╠═a14b3332-47a8-11eb-14e2-196bfbfeecfe
# ╠═a1199be0-47a8-11eb-1971-fdb89c561105
# ╠═a578b040-47a8-11eb-1fbd-f9e5eb680df1
# ╠═cb7ff2d0-47a8-11eb-3bc9-f361ff13bbf0
# ╠═ebdcb53e-47a8-11eb-1660-353aa67bfa6e
# ╠═fc2d09e0-47a8-11eb-22c2-013caff40755
# ╠═03bf5cd0-47a9-11eb-3a9a-cfaf895e6cbf


text:
#' % FIR filter design with Julia
#' % Matti Pastell
#' % 21th April 2016

#' # Introduction

#' This an example of a julia script that can be published using
#' [Weave](http://mpastell.github.io/Weave.jl/latest/usage/).
#' The script can be executed normally using Julia
#' or published to HTML or pdf with Weave.
#' Text is written in markdown in lines starting with "`#'` " and code
#' is executed and results are included in the published document.

#' Notice that you don't need to define chunk options, but you can using
#' `#+`. just before code e.g. `#+ term=True, caption='Fancy plots.'`.
#' If you're viewing the published version have a look at the
#' [source](FIR_design.jl) to see the markup.

#' # FIR Filter Design

#' We'll implement lowpass, highpass and ' bandpass FIR filters. If
#' you want to read more about DSP I highly recommend [The Scientist
#' and Engineer's Guide to Digital Signal
#' Processing](http://www.dspguide.com/) which is freely available
#' online.

#' ## Calculating frequency response

#' DSP.jl package doesn't (yet) have a method to calculate the
#' the frequency response of a FIR filter so we define it:

using Gadfly, DSP

function FIRfreqz(b::Array, w = linspace(0, π, 1024))
    n = length(w)
    h = Array{Complex64}(n)
    sw = 0
    for i = 1:n
      for j = 1:length(b)
        sw += b[j]*exp(-im*w[i])^-j
      end
      h[i] = sw
      sw = 0
    end
    return h
end

#' ## Design Lowpass FIR filter

#' Designing a lowpass FIR filter is very simple to do with DSP.jl, all you
#' need to do is to define the window length, cut off frequency and the
#' window. We will define a lowpass filter with cut off frequency at 5Hz for a signal
#' sampled at 20 Hz.
#' We will use the Hamming window, which is defined as:
#' $w(n) = \alpha - \beta\cos\frac{2\pi n}{N-1}$, where $\alpha=0.54$ and $\beta=0.46$

fs = 20
f = digitalfilter(Lowpass(5, fs = fs), FIRWindow(hamming(61)))
w = linspace(0, pi, 1024)
h = FIRfreqz(f, w)

#' ## Plot the frequency and impulse response

#' The next code chunk is executed in term mode, see the [script](FIR_design.jl) for syntax.
#+ term=true

h_db = log10(abs(h))
ws = w/pi*(fs/2)
plot(y = h_db, x = ws, Geom.line,
      Guide.xlabel("Frequency (Hz)"), Guide.ylabel("Magnitude (db)"))

#' And again with default options

h_phase = unwrap(-atan2(imag(h),real(h)))
plot(y = h_phase, x = ws, Geom.line,
    Guide.xlabel("Frequency (Hz)"), Guide.ylabel("Phase (radians)"))


text:
# # Introduction

# Broadly the Wombling algorithm traverses a 'landscape' (for example species
# richness at different degree squares) and describes the the landscape in terms
# of the *Rate of Change (`m`)* (think slope/gradient) and *Direction of Change
# (`θ`)* (direction of slope). This process is outlined in the figure below.

# ![image info](assets/fig_concept.png)

# Both rate and direction of change are calculated by interpolating points in a
# given search window which (depending on how the points are arranged in space),
# can be done in two ways.

# **Lattice Wombling:** For points that are regularly arranged in space the
# search window is 2 × 2 points.

# **Triangulation Wombling:** For points that irregularly arranged in space
# nearest neigbours are found using delaunay triangualtion and interpolation is
# between 3 points.

# In this example, we will see how the `SpatialBoundaries.jl` package works (as
# well as interpreting the rates and direction of change), by taking a random
# mid-point displacement landscape, and measuring its rate and direction of
# change.

using SpatialBoundaries
using NeutralLandscapes
using StatsPlots

# We will set a few options for the default plots:

default(; dpi=500, size=(600, 600), aspectratio=1, c=:davos, frame=:box)

# The landscape generation is done using the `NeutralLandscapes` package, and we
# will pick a 500x500 grid:

landscape_size = (500, 500)
landscape = rand(MidpointDisplacement(0.75), landscape_size...)

# By default, lattice wombling will assume that the cells have the same size,
# which is 1/n (where n is the number of cells on each side), but  you can
# specify your own `x` and `y` arguments.

# We can take a quick peek at the landscape:

heatmap(landscape)

# Getting the lattice wombling is done with

W = wombling(landscape);

# The resulting `LatticeWomble` object has fields for the rate of change (`m`),
# the direction of change in degrees (`θ`), and the values of the centers of the
# cells at `x` and `y` (note that the grid of rates of change is one cell
# smaller than the original grid!). These points are in a grid so the **Lattice
# Wombling** function was used - note that `wombling()` will select the
# appropriate algorithm based on data input.

# Let's have a look at the rate of change:

heatmap(W.m, c=:bilbao, clim=(0, maximum(W.m)))

# The rate of change informs us on the potential for there to be a boundary
# (zone of change) within a window. Cells with a high rate of change are
# indicative of large differences (changes) in the landscape 'topology' and are
# suggestive of a boundary as we shift from one 'state' to another.

# The *direction* of change is also given, and is expressed a *wind* direction;
# for instance, an angle of 180° means that the value is smaller in the South,
# and larger in the North:

heatmap(W.θ, c=:brocO, clim=(0., 360.))

# The direction of change is *not* the direction the boundary would be if you
# were to draw it on the landscape but rather the direction the rate of change
# is 'moving in'. This means it is possible to think of and use the direction of
# change independently of caluclating boundaries *per se* and can be used to
# inform how the landscape is behaving/changing in a more 'continuous' way as
# opposed to discrete zones/boundaries. For example if changes in species
# richness are more gradual (rate of change is near constant) but the the
# direction of change is consistently South-North (i.e. 180°) we can still infer
# that species richness is 'uniformly' increasing in a South-North direction.

# A note on outputs: The new $x$ and $y$ co-ordinates correspond to latitude and
# longitude respectively.


text:
include("help.jl")
function full_field_in_marks!(r::Robot)
    num_vert = get_num_steps_movements!(r, Sud)
    num_hor = get_num_steps_movements!(r, West)
    #Мы в левом нижнем углу и помним наше число шагов

    side = Ost
    mark_row!(r,side)
    putmarker!(r)

    movements!(r,Sud)
    movements!(r, West)

    #Возвращение в первоначальную позицию
    movements!(r,Ost,num_hor)
    movements!(r,Nord,num_vert)
end

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

function mark_row!(r::Robot,side::HorizonSide) #Движение змейкой
    while isborder(r,side) == false #По горизонтали
        putmarker!(r)
        move!(r,side)
    end
    if (isborder(r,Nord) == false) #По вертикали
        putmarker!(r)
        move!(r,Nord)
        side = invers(side::HorizonSide)
        mark_row!(r,side)
    end
end

function get_num_steps_movements!(r::Robot, side::HorizonSide) #пока нет стенки идем в направлении side и запоминаем кол-во шагов
    num_steps = 0
    while isborder(r, side) == false
        move!(r,side)
        num_steps += 1
    end
    return num_steps
end

function movements!(r::Robot,side::HorizonSide) #Движение до стенки
    while isborder(r,side) == false
        move!(r,side)
    end
end
function movements!(r::Robot,side::HorizonSide,num_steps::Int) #возвращаемся на num_steps шагов
    for _ in 1:num_steps
        move!(r,side)
    end
end

max_stars_repo_name: MCPavlo/mirea-progs

max_issues_repo_name: MCPavlo/mirea-progs

max_forks_repo_name: MCPavlo/mirea-progs


text:
# piecewise_schedule.jl
# original: https://github.com/denizyuret/Knet.jl/blob/master/examples/reinforcement-learning/dqn/piecewise_schedule.jl

linear_interpolation(l, r, α) = l + α * (r - l)

struct PiecewiseSchedule
    endpoints::Vector{Tuple{Int, Float64}}
end

function value(schedule::PiecewiseSchedule, t::Int)
    for ((l_t, l), (r_t, r)) in zip(schedule.endpoints[1:end-1], schedule.endpoints[2:end])
        if l_t <= t < r_t
            α = (t - l_t) / (r_t - l_t)
            return linear_interpolation(l, r, α)
        end
    end
    return schedule.endpoints[end][2]
end

max_stars_repo_name: JuliaTagBot/JuliaGymDemo

max_issues_repo_name: JuliaTagBot/JuliaGymDemo

max_forks_repo_name: JuliaTagBot/JuliaGymDemo


text:
function adaptivesimpsons_inner(
    f::Function,
    a::Real,
    b::Real,
    ε::Real,
    S::Real,
    fa::Real,
    fb::Real,
    fc::Real,
    bottom::Integer,
)
    c = (a + b) / 2
    h = b - a
    d = (a + c) / 2
    g = (c + b) / 2
    fd = f(d)
    fe = f(g)
    Sleft = (h / 12) * (fa + 4 * fd + fc)
    Sright = (h / 12) * (fc + 4 * fe + fb)
    S2 = Sleft + Sright
    if bottom <= 0 || abs(S2 - S) <= 15 * ε
        return S2 + (S2 - S) / 15
    end
    res1 = adaptivesimpsons_inner(
        f,
        a,
        c,
        ε / 2,
        Sleft,
        fa,
        fc,
        fd,
        bottom - 1
    )
    res2 = adaptivesimpsons_inner(
        f,
        c,
        b,
        ε / 2,
        Sright,
        fc,
        fb,
        fe,
        bottom - 1
    )
    return res1 + res2
end

function adaptivesimpsons(
    f::Function,
    a::Real,
    b::Real,
    accuracy::Real = 10e-10,
    max_iterations::Integer = 50,
)
    c = (a + b) / 2
    h = b - a
    fa = f(a)
    fb = f(b)
    fc = f(c)
    S = (h / 6) * (fa + 4 * fc + fb)
    return adaptivesimpsons_inner(
        f,
        a,
        b,
        accuracy,
        S,
        fa,
        fb,
        fc,
        max_iterations
    )
end


text:
# This file was generated, do not modify it. # hide
ŷ = predict(mach, X)
round(rms(ŷ, y), sigdigits=4)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""
    ObjsResources(;
        excluded_ids=nothing,
        ids=nothing,
        wildcard=nothing,
    )

    - excluded_ids::Vector{Any}
    - ids::Vector{Any}
    - wildcard::Bool
"""
mutable struct ObjsResources <: SwaggerModel
    excluded_ids::Any # spec type: Union{ Nothing, Vector{Any} } # spec name: excluded_ids
    ids::Any # spec type: Union{ Nothing, Vector{Any} } # spec name: ids
    wildcard::Any # spec type: Union{ Nothing, Bool } # spec name: wildcard

    function ObjsResources(;excluded_ids=nothing, ids=nothing, wildcard=nothing)
        o = new()
        validate_property(ObjsResources, Symbol("excluded_ids"), excluded_ids)
        setfield!(o, Symbol("excluded_ids"), excluded_ids)
        validate_property(ObjsResources, Symbol("ids"), ids)
        setfield!(o, Symbol("ids"), ids)
        validate_property(ObjsResources, Symbol("wildcard"), wildcard)
        setfield!(o, Symbol("wildcard"), wildcard)
        o
    end
end # type ObjsResources

const _property_map_ObjsResources = Dict{Symbol,Symbol}(Symbol("excluded_ids")=>Symbol("excluded_ids"), Symbol("ids")=>Symbol("ids"), Symbol("wildcard")=>Symbol("wildcard"))
const _property_types_ObjsResources = Dict{Symbol,String}(Symbol("excluded_ids")=>"Vector{Any}", Symbol("ids")=>"Vector{Any}", Symbol("wildcard")=>"Bool")
Base.propertynames(::Type{ ObjsResources }) = collect(keys(_property_map_ObjsResources))
Swagger.property_type(::Type{ ObjsResources }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ObjsResources[name]))}
Swagger.field_name(::Type{ ObjsResources }, property_name::Symbol) =  _property_map_ObjsResources[property_name]

function check_required(o::ObjsResources)
    (getproperty(o, Symbol("ids")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ObjsResources }, name::Symbol, val)
end


text:
using Test
using MLDatasets
using ImageCore

tests = [
    "tst_iris.jl",
    "tst_boston_housing.jl",
    "tst_cifar10.jl",
    "tst_cifar100.jl",
    "tst_mnist.jl",
    "tst_fashion_mnist.jl",
    "tst_svhn2.jl",
]

for t in tests
    @testset "$t" begin
        include(t)
    end
end

#temporary to not stress CI
if !parse(Bool, get(ENV, "CI", "false"))
    @testset "other tests" begin
        # PTBLM
        x, y = PTBLM.traindata()
        x, y = PTBLM.testdata()

        # UD_English
        x = UD_English.traindata()
        x = UD_English.devdata()
        x = UD_English.testdata()
    end
end
nothing


text:
# Cours      : Metaheuristiques
# Algorithme : Resolution par AG de la fonction de Shubert
# Auteur     : Xavier Gandibleux
# Date       : 2017 - rev aout 2020 pour V1.5

@static if VERSION < v"1.5-"
           error("NOT COMPLIANT WITH JULIA < v1.5.0")
        end

using Printf
using Statistics

# Parametres de l'algorithme genetique; ces parametres DOIVENT etre calibres!
nPop         = 10   # nombre d'individus (attention : nombre obligatoirement pair)
nGenerations = 2   # nombre de generations
probCros     = 0.5   # probabilite de crossover
probMut      = 0.5   # probabilite de mutation

# Parametres du codage binaire des variables reelles
lVar         =  15   # nbre de bits codant une variable (=> chromosome de 30 bits)
xinf         = -10   # valeur inferieure des variables
delta        =  20   # range des valeurs des variables
epsilon      =   3   # precision (digits) apres la virgule

# Variables d'environnement
struct t_environment
  # for display's options: false => mute, true => verbose
  displayPrint   :: Bool # display all the prints to follow the algorithm's activity
  displayGraphic :: Bool # display the results on a graphic
end

# Experimentation
nRun         =   1   # nombre de repetition de la resolution du probleme

# Fonction de Shubert
function shubert(x1::Float64,x2::Float64)
    sum1=0.0
    sum2=0.0
    for i=1:5
        sum1 = sum1 + i * cos((i+1)*x1+i)
        sum2 = sum2 + i * cos((i+1)*x2+i)
    end
    return sum1 * sum2
end

# Decoder des variables x1, x2 d'un chromosome
function decoderVariables(x)
    # codage binaire avec bit de poids faible a gauche du mot
    x1 = 0.0
    x2 = 0.0
    for i=1:lVar
        x1=x1+x[i]*2^(i-1)
        x2=x2+x[lVar+i]*2^(i-1)
    end
    x1 = xinf + x1 * delta / (2^lVar - 1)
    x2 = xinf + x2 * delta / (2^lVar - 1)
    return x1, x2
end

# Evaluation de la fonction sur les 2 variables
function evaluerFonction(individu)
    x1,x2 = decoderVariables(individu)
    return shubert(x1,x2)
end

# Generation de la population initiale
function genererPopulationInitiale(env::t_environment, nPop, lVar)
    # les 2 variables sont encodees consecutivement dans un chromosome
    pop = Matrix{Int8}(undef,nPop,2*lVar)
    for i =1:nPop
        pop[i,:]= rand(0:1,2*lVar)
    end

    if env.displayPrint == true
        println("Population initiale ")
        for i=1:nPop
            x1,x2 = decoderVariables(pop[i,:])
            fx1x2 = evaluerFonction(pop[i,:])
            @printf("[%3d]  x1 = % 3.3f  x2 = % 2.3f  |  f(x1,x2) = %8.3f \n",i,x1,x2,fx1x2)
        end
    end

    return pop
end

# Affichage des informations (indice; variables x1, x2; valeur de f(x1,x2) ) d'un individu
function afficherIndividu(env::t_environment, num, individu)
    if env.displayPrint == true
        x1,x2 = decoderVariables(individu)
        findividu = evaluerFonction(individu)
        @printf("Individu  %1d       : x1 = % 3.3f  x2 = % 2.3f  |  f(x1,x2) = %8.3f \n",num,x1,x2,findividu)
    end
end

# Application d'un crossover a deux points
function crossover(env::t_environment, xParent1, xParent2)
    if rand() <= probCros
        # application d'un point de coupe aleatoire sur chaque variable
        coupe1=rand(1:lVar)
        coupe2=rand(lVar+1:2*lVar)

        a1=xParent1[1:coupe1];       b1=xParent1[coupe1+1:lVar]
        a2=xParent2[1:coupe1];       b2=xParent2[coupe1+1:lVar]

        c1=xParent1[lVar+1:coupe2];  d1=xParent1[coupe2+1:2*lVar]
        c2=xParent2[lVar+1:coupe2];  d2=xParent2[coupe2+1:2*lVar]

        xEnfant1= vcat(a1,b2,c1,d2)
        xEnfant2= vcat(a2,b1,c2,d1)

        if env.displayPrint == true
            println("crossover")
            println("  points de coupe   : ", coupe1, " ", coupe2)
            println("  Parents 1 (x1,x2) : ", xParent1[1:lVar], " ", xParent1[lVar+1:2*lVar])
            println("  Parents 2 (x1,x2) : ", xParent2[1:lVar], " ", xParent2[lVar+1:2*lVar])
            println("  Enfants 1 (x1,x2) : ", xEnfant1[1:lVar], " ", xEnfant1[lVar+1:2*lVar])
            println("  Enfants 2 (x1,x2) : ", xEnfant2[1:lVar], " ", xEnfant2[lVar+1:2*lVar])
        end

        return xEnfant1, xEnfant2
    else
        return xParent1, xParent2
    end
end

# Application d'une mutation
function mutation(env::t_environment, individu)
    if rand() <= probMut
        # application d'un flip aleatoire d'un bit sur chaque variable
        ibit1 = rand(1:lVar) ;                     ibit2 = rand(lVar+1:2*lVar)
        individu[ibit1] = (individu[ibit1]+1)%2 ;  individu[ibit2] = (individu[ibit2]+1)%2
        if env.displayPrint == true
            println("Mutation")
            println("  bits mutes  : ", ibit1, " ", ibit2)
            println("  Individu (x1,x2)  : ", individu[1:lVar], " ", individu[lVar+1:2*lVar])
        end
    end
    return individu
end

# Selection d'un parent par tournoi binaire
function selectionTournoiPopulation(env::t_environment, pop)
    i1 = rand(1:size(pop,1)) ;                 i2 = rand(1:size(pop,1))
    xParent1 = pop[i1,:] ;                     xParent2 = pop[i2,:]
    fparent1 = evaluerFonction(xParent1) ;     fparent2 = evaluerFonction(xParent2)

    if env.displayPrint == true
        @printf("Selection : [%3d] %8.3f  <->  [%3d] %8.3f \n", i1, fparent1, i2, fparent2)
    end

    if fparent1 < fparent2
        return i1, xParent1, fparent1
    else
        return i2, xParent2, fparent2
    end
end

# Selection d'un individu par tournoi binaire
function selectionTournoiSurvivant(env::t_environment, individu1, individu2)
    findividu1 = evaluerFonction(individu1) ;   findividu2 = evaluerFonction(individu2)
    if env.displayPrint == true
        @printf("Selection : %8.3f  <->  %8.3f \n", findividu1, findividu2)
    end

    if findividu1 < findividu2
        return individu1
    else
        return individu2
    end
end

# Calcul du fitness moyen de l'ensemble de la population
function fitnessAvgPopulation(pop)
    fitnessAvg = 0.0
    for i=1:size(pop,1)
      fitnessAvg = fitnessAvg + evaluerFonction(pop[i,:])
    end
    return fitnessAvg / size(pop,1)
end

# Calcul du fitness min de l'ensemble de la population
function fitnessMinPopulation(pop)
    fitnessMin = 0.0
    for i=1:size(pop,1)
      fitnessMin = min(fitnessMin, evaluerFonction(pop[i,:]))
    end
    return fitnessMin
end

function plotAvg(allMinPop, allAvgPop)
    figure("bilan valeurs moyennes toutes generations",figsize=(6,6)) # Create a new figure
    title("GA-Shubert | f(x1,x2) moyen")
    xlabel("generations")
    ylabel("f(x1,x2)")

    nPoint = length(allAvgPop)
    x=collect(0:div(nPoint,10):nPoint)
    xticks(x)
    plot(allAvgPop,linestyle="--", lw=0.5, marker="o", ms=2, color="blue", label="f(x1,x2) avg")
    plot(allMinPop,linestyle="--", lw=0.5, marker="_", ms=4, color="red", label="f(x1,x2) min")
    legend(loc=1, fontsize ="small")
end

function plotfct()
    figure("Fonction à minimiser",figsize=(6,6)) # Create a new figure
    zlabel("shubert(x1,x2)")

    n=1000
    x=range(-10, stop=10, length=n) #linspace(-10,10,n)
    xgrid = repeat(x', outer=(n,1)) # repmat(x',n,1)
    ygrid = repeat(x, outer=(1,n))  # repmat(x,1,n)
    z = zeros(n,n)
    for i=1:n
        for j=1:n
            z[i,j]=shubert(xgrid[i,j],ygrid[i,j])
        end
    end
    plot_surface(xgrid, ygrid, z, cmap="nipy_spectral", linewidth=1)
end

function initplotpop(pop)
    fig = figure()
    n=1000
    x=range(-10, stop=10, length=n) #linspace(-10,10,n)
    xgrid = repeat(x', outer=(n,1)) # repmat(x',n,1)
    ygrid = repeat(x, outer=(1,n))  # repmat(x,1,n)
    z = zeros(n,n)
    for i=1:n
        for j=1:n
            z[i,j]=shubert(xgrid[i,j],ygrid[i,j])
        end
    end
    plot_wireframe(xgrid, ygrid, z, linewidth=0.1)

    ax = gca(projection="3d")
    vx1=[]
    vx2=[]
    vf=[]
    for ind=1:nPop
        x1,x2 = decoderVariables(pop[ind,:])
        f = shubert(x1,x2)
        push!(vx1,x1)
        push!(vx2,x2)
        push!(vf,f)
    end
    plot3D(vx1,vx2,vf,".",c="black",markersize=3)
end

function plotpop(couleur, pop)
    vx1=[]
    vx2=[]
    vf=[]
    for ind=1:nPop
        x1,x2 = decoderVariables(pop[ind,:])
        f = shubert(x1,x2)
        push!(vx1,x1)
        push!(vx2,x2)
        push!(vf,f)
    end
    plot3D(vx1,vx2,vf,".",c=couleur,markersize=3)
end

function AlgorithmeGenetique(env)

    # Initialisations
    pop = genererPopulationInitiale(env, nPop, lVar)

    avgPop = fitnessAvgPopulation(pop)
    allAvgPop = [] # liste des fitness moyens pour chaque generation
    push!(allAvgPop, avgPop)

    minGlobal = fitnessMinPopulation(pop)
    allMinPop = [] # liste des fitness min pour chaque generation
    push!(allMinPop, minGlobal)

    println("Fitness pop init = ", avgPop)

    if env.displayGraphic == true
        initplotpop(pop)
    end

    # Generations de l'AG
    for generation = 1:nGenerations

        newPop = zeros(Int8,nPop,2*lVar) # individus de la nouvelle generation

        for couple = 1:div(nPop,2)

            # Selection des parents 1 et 2 dans la populaton par tournoi binaire
            iParent1, parent1, fparent1 = selectionTournoiPopulation(env, pop)
            iParent2, parent2, fparent2 = selectionTournoiPopulation(env, pop)
            if env.displayPrint == true
                @printf("Parents   : [%3d] %8.3f   |   [%3d] %8.3f \n", iParent1, fparent1, iParent2, fparent2)
            end

            # Operateur evolutionnaire : procreation des enfants par crossover
            enfant1, enfant2 = crossover(env, parent1,parent2)
            afficherIndividu(env, 1,enfant1) ;  afficherIndividu(env, 2,enfant2)

            # Operateur evolutionnaire : perturbation des enfants par mutation
            enfant1 = mutation(env, enfant1)
            afficherIndividu(env, 1,enfant1)

            enfant2 = mutation(env, enfant2)
            afficherIndividu(env, 2,enfant2)

            fenfant1 = evaluerFonction(enfant1)
            fenfant2 = evaluerFonction(enfant2)

            newPop[2*couple-1,:] = selectionTournoiSurvivant(env, parent1,parent2)
            newPop[2*couple,:]   = selectionTournoiSurvivant(env, enfant1,enfant2)

            # Archive la meilleure valeur trouvee jusqu'a present
            minGlobal = min(minGlobal, fparent1, fparent2, fenfant1, fenfant2)
            if env.displayPrint == true
                println(" ")
            end
        end
        pop = deepcopy(newPop)

        avgPop = fitnessAvgPopulation(pop)
        push!(allAvgPop, avgPop)
        println("Fitness pop = ", avgPop)
        push!(allMinPop, minGlobal)

        if env.displayGraphic == true
            if generation==5
                plotpop("blue", pop)
            elseif generation==10
                plotpop("orange", pop)
            elseif generation==nGenerations
                plotpop("red", pop)
            end
        end

    end
    return minGlobal, allMinPop, allAvgPop
end

# =============================================================================

# Activation ou pas des sorties ecran
displayPrint=false; displayGraphic=true
env = t_environment(displayPrint,displayGraphic)

allAvgRun = [] # liste des fitness moyens pour chaque run
if env.displayGraphic == true
    using PyPlot
    plotfct()
end

# Go pour l'experimentation numerique
for run = 1:nRun
    minGlobal, allMinPop, allAvgPop = AlgorithmeGenetique(env)
    println("Minimum global trouve = ", minGlobal)
    push!(allAvgRun, allAvgPop[end])

    if env.displayGraphic == true
        plotAvg(allMinPop, allAvgPop)
    end
end

# Resultats
println("Fitness Min / Moy / Max population moyenne sur tous les runs : ")
println(" Min = ", minimum(allAvgRun))
println(" Moy = ", mean(allAvgRun))
println(" Max = ", maximum(allAvgRun))

println("That's all folk!")

max_stars_repo_name: xgandibleux/metaNLOstu

max_issues_repo_name: xgandibleux/metaNLOstu

max_forks_repo_name: xgandibleux/metaNLOstu


text:
using BlackBoxOptim

const opt_methods = (:adaptive_de_rand_1_bin, :adaptive_de_rand_1_bin_radiuslimited, :separable_nes, :xnes, :de_rand_1_bin, :de_rand_2_bin, :de_rand_1_bin_radiuslimited, :de_rand_2_bin_radiuslimited, :random_search, :generating_set_search, :probabilistic_descent)

# opt_fields = [(:r,(0.01,2.0)),(:volfrac,(0.,0.14)),(:ρ,(0.0,100.))]
# fix_fields = [(:c,[1.0+0.0im,1.0+0.0im])]
# optf = f_slow1
# medium = Medium{Float64}(1.0,1.)
# ωs = 0.01:0.01:1.0
# num_species=2
# MaxTime=5
# method = :adaptive_de_rand_1_bin

function optimal_species(optf, medium, ωs;
                opt_fields = [(:r,(0.1,2.0)), (:volfrac,(0.01,0.1))],
                fix_fields = [(:ρ,[0.0,Inf]), (:c,[1.0+0.0im,1.0+0.0im])],
                num_species=2, MaxTime=100., method=:choose_method,
                opt_kws...
)

  function F(; ρ=[0.,Inf], r=[0.5,0.5], c=[1.0+0.0im,1.0+0.0im], volfrac=[0.1,0.1])
    N = num_species
    sps = [Specie(ρ[i], r[i], c[i]; volfrac=volfrac[i]) for i=1:N]
    kTs = wavenumber_low_volfrac(ωs, medium, sps)
    # kTs = [wavenumber_low_volfrac(ω, medium, sps) for ω in ωs]
    optf(sps, medium, ωs, kTs)
  end
  function optF(x)
    var_fields = map(1:length(opt_fields)) do n
      (opt_fields[n][1], x[(1+(n-1)*num_species):(n*num_species)])
    end
    F(; var_fields..., fix_fields...)
  end
  if method == :choose_method
    fits =[
      best_fitness(bboptimize(optF; MaxTime = min(100.0,MaxTime), Method = m,
        SearchRange =  repeat([o[2] for o in opt_fields], inner=[num_species])))
    for m in opt_methods]
    method = opt_methods[findmin(fits)[2]]
    println("For MaxTime=$MaxTime, the optimal method was = $method.")
  end

  res = bboptimize(optF; MaxTime = MaxTime, Method=method,
    SearchRange =  repeat([o[2] for o in opt_fields], inner=[num_species]))
  parameters = reshape(best_candidate(res), num_species, length(opt_fields))
  # include the fixed paramters
  parameters = hcat(parameters, hcat([f[2] for f in fix_fields]...))
  fields = [f[1] for f in (opt_fields..., fix_fields...)]
  [ Specie(Float64; zip(fields, parameters[n,:])...) for n =1:num_species]
end


text:
@doc raw"""
    ps_cost(canonical_model::CanonicalModel,
                variable::JuMP.Containers.DenseAxisArray{JV},
                cost_component::Float64,
                dt::Float64,
                sign::Float64) where {JV <: JuMP.AbstractVariableRef}

Returns linear cost terms for sum of variables with common factor to be used for cost expression for canonical model.

# Equation

``` gen_cost = sum(variable)*cost_component ```

# LaTeX

`` cost = dt\times sign\sum_{i\in I} c x_i ``

Returns:

``` sign*gen_cost*dt ```

# Arguments

* canonical_model::CanonicalModel : the canonical model built in PowerSimulations
* variable::JuMP.Containers.DenseAxisArray{JV} : variable array
* cost_component::Float64 : cost to be associated with variable
* dt::Float64 : fraction of hour
* sign::Float64 : positive or negative sign to be associated cost term
"""
function ps_cost(canonical_model::CanonicalModel,
                variable::JuMP.Containers.DenseAxisArray{JV},
                cost_component::Float64,
                dt::Float64,
                sign::Float64) where {JV<:JuMP.AbstractVariableRef}

    gen_cost = sum(variable)*cost_component

    return sign*gen_cost*dt

end

@doc raw"""
    ps_cost(canonical_model::CanonicalModel,
                variable::JuMP.Containers.DenseAxisArray{JV},
                cost_component::PSY.VariableCost{Float64},
                dt::Float64,
                sign::Float64) where {JV <: JuMP.AbstractVariableRef}

Returns linear cost terms for sum of variables with common factor to be used for cost expression for canonical model.
Does this by calling ```ps_cost``` that has Float64 cost component input.

Returns:

``` ps_cost(canonical_model, variable, PSY.get_cost(cost_component), dt, sign) ```

# Arguments

* canonical_model::CanonicalModel : the canonical model built in PowerSimulations
* variable::JuMP.Containers.DenseAxisArray{JV} : variable array
* cost_component::PSY.VariableCost{Float64} : container for cost to be associated with variable
* dt::Float64 : fraction of hour
* sign::Float64 : positive or negative sign to be associated cost term
"""
function ps_cost(canonical_model::CanonicalModel,
                variable::JuMP.Containers.DenseAxisArray{JV},
                cost_component::PSY.VariableCost{Float64},
                dt::Float64,
                sign::Float64) where {JV<:JuMP.AbstractVariableRef}

    return  ps_cost(canonical_model, variable, PSY.get_cost(cost_component), dt, sign)

end

@doc raw"""
    ps_cost(canonical_model::CanonicalModel,
                variable::JuMP.Containers.DenseAxisArray{JV},
                cost_component::PSY.VariableCost{NTuple{2, Float64}}
                dt::Float64,
                sign::Float64) where {JV <: JuMP.AbstractVariableRef}

Returns quadratic cost terms for sum of variables with common factor to be used for cost expression for canonical model.

# Equation

``` gen_cost = dt*sign*(sum(variable.^2)*cost_component[1] + sum(variable)*cost_component[2]) ```

# LaTeX

`` cost = dt\times sign (sum_{i\in I} c_1 v_i^2 + sum_{i\in I} c_2 v_i ) ``

for quadratic factor large enough. Otherwise

``` return ps_cost(canonical_model, variable, cost_component[2], dt, 1.0) ```

Returns ```gen_cost```

# Arguments

* canonical_model::CanonicalModel : the canonical model built in PowerSimulations
* variable::JuMP.Containers.DenseAxisArray{JV} : variable array
* cost_component::PSY.VariableCost{NTuple{2, Float64}} : container for quadratic and linear factors
* sign::Float64 : positive or negative sign to be associated cost term
"""
function ps_cost(canonical_model::CanonicalModel,
                 variable::JuMP.Containers.DenseAxisArray{JV},
                 cost_component::PSY.VariableCost{NTuple{2, Float64}},
                 dt::Float64,
                 sign::Float64) where {JV<:JuMP.AbstractVariableRef}

    if cost_component[1] >= eps()
        gen_cost = sum(variable.^2)*cost_component[1] + sum(variable)*cost_component[2]
        return sign*gen_cost*dt
    else
        return ps_cost(canonical_model, variable, cost_component[2], dt, 1.0)
    end

end

@doc raw"""
    _pwlparamcheck(cost_)

Returns True/False depending on compatibility of the cost data with the linear implementation method

Returns ```flag```

# Arguments

* cost_::PSY.VariableCost{NTuple{2, Float64}} : container for quadratic and linear factors
"""
function _pwlparamcheck(cost_)
    flag = true

    if abs(cost_[1][1]/cost_[1][2] - ((cost_[2][1] - cost_[1][1])/(cost_[2][2] - cost_[1][2]))) > 1.0
        flag =false
    end
    l = length(cost_)
    for i in 1:(l-2)
        if ((cost_[i+1][1] - cost_[i][1])/(cost_[i+1][2] - cost_[i][2])) > ((cost_[i+2][1] - cost_[i+1][1])/(cost_[i+2][2] - cost_[i+1][2]))
            flag = false
        end
    end
    return flag
end

@doc raw"""
    _pwlgencost_sos(canonical_model::CanonicalModel,
                variable::JuMP.Containers.DenseAxisArray{JV},
                cost_component::PSY.VariableCost{NTuple{2, Float64}}) where {JV <: JuMP.AbstractVariableRef}

Returns piecewise cost expression using SOS Type-2 implementation for canonical model.

# Equations

``` variable = sum(sos_var[i]*cost_component[2][i])```

``` gen_cost = sum(sos_var[i]*cost_component[1][i]) ```

# LaTeX

`` variable = (sum_{i\in I} c_{2, i} sos_i) ``

`` gen_cost = (sum_{i\in I} c_{1, i} sos_i) ``

Returns ```gen_cost```

# Arguments

* canonical_model::CanonicalModel : the canonical model built in PowerSimulations
* variable::JuMP.Containers.DenseAxisArray{JV} : variable array
* cost_component::PSY.VariableCost{NTuple{2, Float64}} : container for quadratic and linear factors
"""
function _pwlgencost_sos(canonical_model::CanonicalModel,
        variable::JV,
        cost_component::Vector{NTuple{2, Float64}}) where {JV<:JuMP.AbstractVariableRef}

    gen_cost = JuMP.GenericAffExpr{Float64, _variable_type(canonical_model)}()
    pwlvars = JuMP.@variable(canonical_model.JuMPmodel, [i = 1:length(cost_component)],
                            base_name = "{$(variable)}_{sos}",
                            start = 0.0, lower_bound = 0.0, upper_bound = 1.0)

    JuMP.@constraint(canonical_model.JuMPmodel, sum(pwlvars)  == 1.0)
    JuMP.@constraint(canonical_model.JuMPmodel, pwlvars in MOI.SOS2(collect(1:length(pwlvars))))
    for (ix, var) in enumerate(pwlvars)
        JuMP.add_to_expression!(gen_cost, cost_component[ix][1] * var)
    end

    JuMP.@constraint(canonical_model.JuMPmodel, variable ==
                        sum([var*cost_component[ix][2] for (ix, var) in enumerate(pwlvars) ]) )

    return gen_cost

end

@doc raw"""
    _pwlgencost_sos(canonical_model::CanonicalModel,
                variable::JuMP.Containers.DenseAxisArray{JV},
                cost_component::PSY.VariableCost{NTuple{2, Float64}}) where {JV <: JuMP.AbstractVariableRef}

Returns piecewise cost expression using linear implementation for canonical model.

# Equations

``` 0 <= pwl_var[i] <= (cost_component[2][i] - cost_component[2][i-1])```

``` variable = sum(pwl_var[i])```

``` gen_cost = sum(pwl_var[i]*cost_component[1][i]/cost_component[2][i]) ```

# LaTeX
`` 0 <= pwl_i <= (c_{2, i} - c_{2, i-1})``

`` variable = (sum_{i\in I} pwl_i) ``

`` gen_cost = (sum_{i\in I}  pwl_i) c_{1, i}/c_{2, i} ``

Returns ```gen_cost```

# Arguments

* canonical_model::CanonicalModel : the canonical model built in PowerSimulations
* variable::JuMP.Containers.DenseAxisArray{JV} : variable array
* cost_component::PSY.VariableCost{NTuple{2, Float64}} : container for quadratic and linear factors
"""
function _pwlgencost_linear(canonical_model::CanonicalModel,
        variable::JV,
        cost_component::Vector{NTuple{2, Float64}}) where {JV<:JuMP.AbstractVariableRef}

    gen_cost = JuMP.GenericAffExpr{Float64, _variable_type(canonical_model)}()
    upperbound(i) = (i == 1 ? cost_component[i][2] : (cost_component[i][2] - cost_component[i-1][2]))
    pwlvars = JuMP.@variable(canonical_model.JuMPmodel, [i = 1:length(cost_component)],
                            base_name = "{$(variable)}_{pwl}", start = 0.0,
                            lower_bound = 0.0, upper_bound = upperbound(i))

    for (ix, pwlvar) in enumerate(pwlvars)
        if ix == 1
            JuMP.add_to_expression!(gen_cost, cost_component[ix][1] * (pwlvar / cost_component[ix][2]))
        else
            JuMP.add_to_expression!(gen_cost, (cost_component[ix][1] - cost_component[ix-1][1]) *
                                            (pwlvar/(cost_component[ix][2] - cost_component[ix-1][2])))
        end
    end

    c = JuMP.@constraint(canonical_model.JuMPmodel, variable == sum([pwlvar for (ix, pwlvar) in enumerate(pwlvars) ]) )

    return gen_cost

end

@doc raw"""
    _gen_cost(cost_)

Returns JuMP expression for a piecewise linear cost function depending on the data compatibility.

Returns ```gen_cost```

# Arguments

* canonical_model::CanonicalModel : the canonical model built in PowerSimulations
* variable::JuMP.Containers.DenseAxisArray{JV} : variable array
* cost_component::PSY.VariableCost{NTuple{2, Float64}} : container for quadratic and linear factors
"""
function _pwl_cost(canonical_model::CanonicalModel,
                    variable::JV,
                    cost_component::Vector{NTuple{2, Float64}}) where {JV<:JuMP.AbstractVariableRef}

    # If array is full of tuples with zeros return 0.0
    all(iszero.(last.(cost_component))) && return 0.0

    if !_pwlparamcheck(cost_component)
        @warn("The cost function provided for $(variable) device is not compatible with a linear PWL cost function.
        An SOS-2 formulation will be added to the model.
        This will result in additional binary variables added to the model.")
        gen_cost = _pwlgencost_sos(canonical_model, variable, cost_component)
    else
        gen_cost = _pwlgencost_linear(canonical_model, variable, cost_component)
    end
    return gen_cost
end

@doc raw"""
    ps_cost(canonical_model::CanonicalModel,
                 variable::JuMP.Containers.DenseAxisArray{JV},
                 cost_component::PSY.VariableCost{Vector{NTuple{2, Float64}}},
                 dt::Float64,
                 sign::Float64) where {JV<:JuMP.AbstractVariableRef}

Creates piecewise linear cost function using a sum of variables and expression with sign and time step included.

# Expression

```JuMP.add_to_expression!(gen_cost, c)```

Returns sign*gen_cost*dt

# LaTeX

``cost = sign\times dt \sum_{v\in V} c_v``

where ``c_v`` is given by

`` c_v = \sum_{i\in Ix} \frac{y_i - y_{i-1}}{x_i - x_{i-1}} v^{p.w.}_i ``

# Arguments

* canonical_model::CanonicalModel : the canonical model built in PowerSimulations
* variable::JuMP.Containers.DenseAxisArray{JV} : variable array
* cost_component::PSY.VariableCost{Vector{NTuple{2, Float64}}}
* dt::Float64 : fraction of hour
* sign::Float64 : positive or negative sign to be associated cost term
"""
function ps_cost(canonical_model::CanonicalModel,
                 variable::JuMP.Containers.DenseAxisArray{JV},
                 cost_component::PSY.VariableCost{Vector{NTuple{2, Float64}}},
                 dt::Float64,
                 sign::Float64) where {JV<:JuMP.AbstractVariableRef}

    gen_cost = JuMP.GenericAffExpr{Float64, _variable_type(canonical_model)}()
    cost_array = cost_component.cost
    for var in variable
        c = _pwl_cost(canonical_model, var, cost_array)
        JuMP.add_to_expression!(gen_cost, c)
    end

    return sign*gen_cost*dt

end

@doc raw"""
    add_to_cost(canonical_model::CanonicalModel,
                     devices::D,
                     var_name::Symbol,
                     cost_symbol::Symbol,
                     sign::Float64 = 1.0) where {D<:IS.FlattenIteratorWrapper{<:PSY.Device}}

Adds cost expression for each device using appropriate call to ```ps_cost```.

# Expression

for d in devices

```    cost_expression = ps_cost(canonical_model,
                              variable[PSY.get_name(d), :],
                              getfield(PSY.get_op_cost(d), cost_symbol),
                              dt,
                              sign) ```
``` canonical_model.cost_function += cost_expression ```

# LaTeX

`` COST = \sum_{d\in D} cost_d ``

# Arguments

* canonical_model::CanonicalModel : the canonical model built in PowerSimulations
* devices::D : set of devices
* var_name::Symbol : name of variable
* cost_symbol::Symbol : symbol associated with costx
"""
function add_to_cost(canonical_model::CanonicalModel,
                     devices::D,
                     var_name::Symbol,
                     cost_symbol::Symbol,
                     sign::Float64 = 1.0) where {D<:IS.FlattenIteratorWrapper{<:PSY.Device}}

    resolution = model_resolution(canonical_model)
    dt = Dates.value(Dates.Minute(resolution))/60
    variable = var(canonical_model, var_name)

    for d in devices
        cost_component = getfield(PSY.get_op_cost(d), cost_symbol)
        cost_expression = ps_cost(canonical_model,
                                  variable[PSY.get_name(d), :],
                                  cost_component,
                                  dt,
                                  sign)
        T_ce = typeof(cost_expression)
        T_cf = typeof(canonical_model.cost_function)
        if  T_cf<:JuMP.GenericAffExpr && T_ce<:JuMP.GenericQuadExpr
            canonical_model.cost_function += cost_expression
        else
            JuMP.add_to_expression!(canonical_model.cost_function, cost_expression)
        end
    end

    return

end


text:
# Add local directory to path
push!(LOAD_PATH, ".")

# Run local _init.jl if present
isfile("_init.jl") && include(joinpath(pwd(), "_init.jl"))

# Revise and run without reloading manually
@schedule begin
    sleep(0.1)
    @eval using Revise
end

max_stars_repo_name: briochemc/dotfiles

max_issues_repo_name: briochemc/dotfiles

max_forks_repo_name: briochemc/dotfiles


text:
module huobi

using HTTP
using JSON
using Nettle
using Base64
using Dates
using Printf
using DataStructures

global const BASE_URL = "https://api.huobi.pro"

publickey() = ENV["HUOBI-PUBLIC-KEY"]
secretkey() = ENV["HUOBI-SECRET-KEY"]

function signature(
    method::AbstractString,
    host::AbstractString,
    path::AbstractString,
    query::AbstractDict,
    secret::AbstractString
)
    pairs = join(string.(keys(query),"=", HTTP.URIs.escapeuri.(values(query))), "&")
    str_query = join([method, host, path, pairs], "\n")

    hmac = HMACState("sha256", secret)
    Nettle.update!(hmac, str_query)

    Base64.base64encode(Nettle.digest!(hmac))
end

function request(
    method::String,
    path::String,
    query::AbstractDict,
    public::Union{Missing,String},
    secret::Union{Missing,String},
)

    @debug "request" method path (;zip(keys(query) .|> Symbol, values(query))...)...

    !(method ∈ ["GET", "POST"]) &&
        error("unknown method $method")

    timestamp = Dates.format(now(UTC), "yyyy-mm-ddTHH:MM:SS")

    params = [
        "AccessKeyId" => public,
        "SignatureMethod" => "HmacSHA256",
        "SignatureVersion" => 2,
        "Timestamp" => timestamp
    ] |> OrderedDict

    path = "/" * path
    host = HTTP.URI(BASE_URL).host
    url = string(BASE_URL, path)

    isequal(method, "GET") &&
        merge!(params, query)

    params["Signature"] =
        signature(method, host, path, params, secret)

    response = if method == "GET"
        HTTP.request(method, url, query=params, status_exception=false)
    else
        jsn = isempty(query) ? "" : JSON.json(query)
        HTTP.request(method, url, query=params, jsn, status_exception=false)
    end

    json = JSON.parse(String(response.body))

    if haskey(json, "status")
        if json["status"] == "error"
            error(string(
                json["err-msg"], " ",
                JSON.json(query)
            ))
        end
    end

    if haskey(json, "code")
        if json["code"] != 200
            error(string(
                JSON.json(query)
            ))
        end
    end

    @debug "response" method path (;zip(keys(json) .|> Symbol, values(json))...)...

    return json
end

request(
    method::String,
    path::String,
    query::AbstractDict = OrderedDict()
) = request(method, path, query, publickey(), secretkey())

#Signed requests

accounts() =
    request("GET", "v1/account/accounts", OrderedDict())

function balances(type::String)
    response = accounts()
    filter!(x-> x["type"] == "$type", response["data"])
    account_id = map(x -> x["id"], response["data"])[1]
    request("GET", "v1/account/accounts/$account_id/balance", OrderedDict())
end

function balance(currency::String, type::String)
    response = balances("$type")
    for cur in response["data"]["list"]
        if currency == cur["currency"]
            return cur["balance"]
        else
            nothing
        end
    end
end

#Public requests
function details()
    endpoint = "/v1/common/symbols"
    response = HTTP.request("GET", string(BASE_URL, endpoint))
    json = JSON.parse(String(response.body))
    symbols = map(
        x -> x["symbol"],
        filter(
            x -> x["state"] == "online" && x["api-trading"] == "enabled",
            json["data"],
        ),
    )
    base_currency = map(
        x -> x["base-currency"],
        filter(
            x -> x["state"] == "online" && x["api-trading"] == "enabled",
            json["data"],
        ),
    )
    quote_currency = map(
        x -> x["quote-currency"],
        filter(
            x -> x["state"] == "online" && x["api-trading"] == "enabled",
            json["data"],
        ),
    )
    return (
        symbols = symbols,
        base_currency = base_currency,
        quote_currency = quote_currency,
    )
end

function stats_24hr(symbol::String)
    path = "market/detail/merged"
    url = join([BASE_URL, path], "/")
    params = "symbol=$(lowercase(symbol))"
    curl = join([url, params], "?")
    response = HTTP.request("GET", curl)
    jsn = JSON.parse(String(response.body))
    return jsn["tick"]["vol"]
end

function order_book(symbol::String, type::String="step0")
    _path = "market/depth"
    url = join([BASE_URL, _path], "/")
    params = "?symbol=$(lowercase(symbol))&type=$type"
    response = HTTP.request("GET", url * params)
    json = JSON.parse(String(response.body))
    return json["tick"]
end

function klines(symbol::String, period::String, size::Int)
    _path = "market/history/kline"
    point = join([BASE_URL, _path], "/")
    params = "period=$period&size=$size&symbol=$(lowercase(symbol))"
    curl = join([point, params], "?")
    response = HTTP.request("GET", curl)
    json = JSON.parse(String(response.body))
    return json["data"]
end

function historical_trades(symbol::String, size::Int)
    _path = "market/history/trade"
    point = join([BASE_URL, _path], "/")
    params = "symbol=$(lowercase(symbol))&size=$size"
    curl = join([point, params], "?")
    response = HTTP.request("GET", curl)
    json = JSON.parse(String(response.body))
    return json["data"]
end

end


text:
instr = readlines("input.txt")
reg = filter(x -> x != "1", unique(map(x -> x[2], split.(instr))))
cmd = Dict("set" => "=", "add" => "+=", "mul" => "*=", "mod" => "%=")
get(x) = x in reg ? eval(Symbol(x)) : parse(Int, x)

foreach(x -> eval(parse("$x = 0")), reg)
pointer = 1
last_sound = 0
while pointer <= length(instr)
   d = split(instr[pointer])
   if d[1] == "snd"
       last_sound = get(d[2])
   elseif d[1] == "jgz" && get(d[2]) > 0
       pointer += get(d[3])
       continue
   elseif d[1] == "rcv" && get(d[2]) != 0
       println(last_sound)
       if last_sound != 0
           break
       end
   elseif d[1] != "jgz" && d[2] != "rcv"
       eval(parse("$(d[2]) $(cmd[d[1]]) $(d[3])"))
   end
   pointer += 1
end


text:
module ModuleWithOneInlineTest
using Base.Test

let
  @test true
end

end


text:
type RenderTexture
    ptr::Ptr{Void}
    _view::View

    function RenderTexture(ptr::Ptr{Void})
        r = new(ptr)
        finalizer(r, destroy)
        r
    end
end

function RenderTexture(width::Integer, height::Integer, depth_buffer::Bool=false)
    RenderTexture(ccall((:sfRenderTexture_create, libcsfml_graphics), Ptr{Void}, (UInt32, UInt32, Int32,), width, height, depth_buffer))
end

function destroy(texture::RenderTexture)
    ccall((:sfRenderTexture_destroy, libcsfml_graphics), Void, (Ptr{Void},), texture.ptr)
end

function get_size(texture::RenderTexture)
    return to_vec2i(ccall((:sfRenderTexture_getSize, libcsfml_graphics), Vector2u, (Ptr{Void},), texture.ptr))
end

function set_active(texture::RenderTexture, active::Bool)
    ccall((:sfRenderTexture_setActive, libcsfml_graphics), Void, (Ptr{Void}, Int32,), texture.ptr, active)
end

function display(texture::RenderTexture)
    ccall((:sfRenderTexture_display, libcsfml_graphics), Void, (Ptr{Void},), texture.ptr)
end

function clear(texture::RenderTexture, color::Color)
    ccall((:sfRenderTexture_clear, libcsfml_graphics), Void, (Ptr{Void}, Color,), texture.ptr, color)
end

function set_view(texture::RenderTexture, view::View)
    ccall((:sfRenderTexture_setView, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Void},), texture.ptr, view.ptr)
    texture._view = view
end

function get_view(texture::RenderTexture)
    View(ccall((:sfRenderTexture_getView, libcsfml_graphics), Ptr{Void}, (Ptr{Void},), texture.ptr))
end

function get_default_view(texture::RenderTexture)
    View(ccall((:sfRenderTexture_getDefaultView, libcsfml_graphics), Ptr{Void}, (Ptr{Void},), texture.ptr))
end

function get_viewport(texture::RenderTexture)
    ccall((:sfRenderTexture_getViewport, libcsfml_graphics), IntRect, (Ptr{Void},), texture.ptr)
end

function pixel2coords(texture::RenderTexture, point::Vector2, targetview::View=get_view(texture))
    point = to_vec2i(point)
    ccall((:sfRenderTexture_mapPixelToCoords, libcsfml_graphics), Vector2f, (Ptr{Void}, Vector2i, Ptr{Void},), texture.ptr, point, targetview.ptr)
end

function coords2pixel(texture::RenderTexture, point::Vector2, targetview=get_view(texture))
    point = to_vec2f(point)
    ccall((:sfRenderTexture_mapCoordsToPixel, libcsfml_graphics), Vector2i, (Ptr{Void}, Vector2f, Ptr{Void},), texture.ptr, point, targetview.ptr)
end

function draw(texture::RenderTexture, object::Drawable, shader::Shader)
    state = RenderStates(shader)
    draw(texture, object, state)
end

function draw(texture::RenderTexture, object::Drawable, state::RenderStates)
    draw(texture, object, state.ptr)
end

function draw(texture::RenderTexture, object::Drawable)
    draw(texture, object, C_NULL)
end

function draw(texture::RenderTexture, object::CircleShape, renderStates::Ptr{Void})
    ccall((:sfRenderTexture_drawCircleShape, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void},), texture.ptr, object.ptr, renderStates)
end

function draw(texture::RenderTexture, object::RectangleShape, renderStates::Ptr{Void})
    ccall((:sfRenderTexture_drawRectangleShape, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void},), texture.ptr, object.ptr, renderStates)
end

function draw(texture::RenderTexture, object::Sprite, renderStates::Ptr{Void})
    ccall((:sfRenderTexture_drawSprite, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void},), texture.ptr, object.ptr, renderStates)
end

function draw(texture::RenderTexture, object::RenderText, renderStates::Ptr{Void})
    ccall((:sfRenderTexture_drawText, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void},), texture.ptr, object.ptr, renderStates)
end

function draw(texture::RenderTexture, object::ConvexShape, renderStates::Ptr{Void})
    ccall((:sfRenderTexture_drawConvexShape, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void},), texture.ptr, object.ptr, renderStates)
end

function draw(texture::RenderTexture, object::VertexArray, renderStates::Ptr{Void})
    ccall((:sfRenderTexture_drawVertexArray, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void},), texture.ptr, object.ptr, renderStates)
end

function get_texture(texture::RenderTexture)
    return Texture(ccall((:sfRenderTexture_getTexture, libcsfml_graphics), Ptr{Void}, (Ptr{Void},), texture.ptr))
end

function set_smooth(texture::RenderTexture, smooth::Bool)
    ccall((:sfRenderTexture_setSmooth, libcsfml_graphics), Void, (Ptr{Void}, Int32,), texture.ptr, smooth)
end

function is_smooth(texture::RenderTexture)
    Bool(ccall((:sfRenderText_isSmooth, libcsfml_graphics), Int32, (Ptr{Void},), texture.ptr))
end

function set_repeated(texture::RenderTexture, repeated::Bool)
    ccall((:sfRenderTexture_setRepeated, libcsfml_graphics), Void, (Ptr{Void}, Int32,), texture.ptr, repeated)
end

function is_repeated(texture::RenderTexture)
    Bool(ccall((:sfRenderTexture_isRepeated, libcsfml_graphics), Int32, (Ptr{Void},), texture.ptr))
end


text:
"""
    PreImputingMatrix{T <: Number, U <: AbstractMatrix{T}, V <: AbstractVector{T}} <: AbstractMatrix{T}

A parametrized matrix that fills in elements from a default vector of parameters whenever a `missing`
element is encountered during multiplication.

# Examples
```jldoctest
julia> A = PreImputingMatrix(ones(2, 2), -ones(2))
2×2 PreImputingMatrix{Float64, Matrix{Float64}, Vector{Float64}}:
W:
 1.0  1.0
 1.0  1.0

ψ:
 -1.0  -1.0

julia> A * [0 1; missing -1]
2×2 Matrix{Float64}:
 -1.0  0.0
 -1.0  0.0
```

See also: [`PreImputingMatrix`](@ref).
"""
struct PreImputingMatrix{T <: Number, U <: AbstractMatrix{T}, V <: AbstractVector{T}} <: AbstractMatrix{T}
    W::U
    ψ::V
end

Flux.@functor PreImputingMatrix

"""
    PreImputingMatrix(W::AbstractMatrix{T}, ψ=zeros(T, size(W, 2))) where T

Construct a [`PreImputingMatrix`](@ref) with multiplication parameters `W` and default parameters `ψ`.

# Examples
```jldoctest
julia> PreImputingMatrix([1 2; 3 4])
2×2 PreImputingMatrix{Int64, Matrix{Int64}, Vector{Int64}}:
W:
 1  2
 3  4

ψ:
 0  0
```

See also: [`PostImputingMatrix`](@ref).
"""
PreImputingMatrix(W::AbstractMatrix{T}) where T = PreImputingMatrix(W, zeros(T, size(W, 2)))

Flux.@forward PreImputingMatrix.W Base.size, Base.getindex, Base.setindex!, Base.firstindex, Base.lastindex

Base.hcat(As::PreImputingMatrix...) = PreImputingMatrix(hcat((A.W for A in As)...), vcat((A.ψ for A in As)...))
function Base.vcat(As::PreImputingMatrix...)
    ArgumentError("It doesn't make sense to vcat PreImputingMatrices") |> throw
end

A::PreImputingMatrix * b::AbstractVector = (_check_mul(A, b); _mul(A, b))
Zygote.@adjoint A::PreImputingMatrix * b::AbstractVector = (_check_mul(A, b); Zygote.pullback(_mul, A, b))
A::PreImputingMatrix * B::AbstractMatrix = (_check_mul(A, B); _mul(A, B))
Zygote.@adjoint A::PreImputingMatrix * B::AbstractMatrix = (_check_mul(A, B); Zygote.pullback(_mul, A, B))

_mul(A::PreImputingMatrix, B::AbstractVecOrMat) = A.W * B
_mul(A::PreImputingMatrix, ::AbstractVector{Missing}) = A.W * A.ψ
_mul(A::PreImputingMatrix, B::AbstractMatrix{Missing}) = repeat(A.W * A.ψ, 1, size(B, 2))
_mul(A::PreImputingMatrix, B::AbstractVecOrMat{Maybe{T}}) where {T <: Number} = A.W * _mul_maybe(A.ψ, B)

_mul_maybe(ψ, B) = _impute_row(ψ, B)[1]
function ChainRulesCore.rrule(::typeof(_mul_maybe), ψ, B)
    X, m = _impute_row(ψ, B)
    X, Δ -> (NoTangent(), @thunk(vec(sum(.!m .* Δ, dims=2))), @thunk(m .* Δ))
end
function _impute_row(ψ, B)
    m = .!ismissing.(B)
    X = similar(ψ, size(B))
    X .= ψ
    X[m] = @view B[m]
    X, m
end


text:
@doc """
Bellman operator.

# Arguments

* `Q::Array{V,2}`: The q-function array. Will be modified in-place.
* `value::Array{V,1}`: The previous value function.
* `transition::Array{P,3}`: Transition probability matrices.
* `reward::Array{R,2}`: The state-based rewards.
* `δ::Float64`: Discount fatcor.

# Returns

* `Q::Array{V,2}`: The modified q-function. Use it to calculate the new
  value-function and current policy.

# WARNING

Will not check that transition and reward specify a valid MDP. This method is
meant to be called by other methods that check the validity of the transition
and reward arrays. Consider using one of the MDP types instead.

""" ->
function bellman!(
    Q::AbstractQFunction,
    value_prev::Vector,
    P::AbstractTransitionProbability,
    R::AbstractReward,
    δ,
)
    0 < δ <= 1 || error("δ musy be in the interval (0, 1].")
    S = num_states(P)
    A = num_actions(P)
    S == num_states(R) == num_states(Q) == length(value_prev) || throw(DimensionMismatch())
    V = valuetype(Q)
    @inbounds for a = 1:A
        for s1 = 1:S
            Σ = zero(V)
            @simd for s0 = 1:S
                Σ += convert(V, probability(P, s0, s1, a)*value_prev[s0])
            end
            setvalue!(Q, convert(V, reward(R, s1, a) + δ*Σ), s1, a)
        end
    end
    return Q
end

function bellman!(Q, value_prev, P, R, δ)
    qfunction = bellman!(QFunction(Q), value_prev, TransitionProbability(P), Reward(R), δ)
    return qfunction.array
end

bellman(Q, value_prev, P, R, δ) =
    bellman!(copy(Q), value_prev, P, R, δ)

function bellman!(value, policy, value_prev, P, R, δ)
    qfunction = bellman!(QFunction(value, policy), value_prev, TransitionProbability(P), Reward(R), δ)
    return qfunction.value, qfunction.policy
end

bellman(value, policy, value_prev, P, R, δ) =
    bellman!(copy(value), copy(policy), value_prev, P, R, δ)


text:
Base.:(==)(a::IndicesArray, b::IndicesArray) = parent(a) == parent(b)
Base.:(==)(a::AbstractArray, b::IndicesArray) = a == parent(b)
Base.:(==)(a::IndicesArray, b::AbstractArray) = parent(a) == b

Base.isequal(a::IndicesArray, b::IndicesArray) = isequal(parent(a), parent(b))
Base.isequal(a::AbstractArray, b::IndicesArray) = isequal(a, parent(b))
Base.isequal(a::IndicesArray, b::AbstractArray) = isequal(parent(a), b)


text:
# Dominique Orban, <dominique.orban@gerad.ca>
# Montreal, QC, November 2016-January 2017.

export lslq

"""
    lslq(A, b, λ=0.0)

Solve the regularized linear least-squares problem

    minimize ‖b - Ax‖₂² + λ² ‖x‖₂²

using the LSLQ method, where λ ≥ 0 is a regularization parameter.
LSLQ is formally equivalent to applying SYMMLQ to the normal equations

    (A'A + λ² I) x = A'b

but is more stable.

#### Main features

* the solution estimate is updated along orthogonal directions
* the norm of the solution estimate ‖xᴸₖ‖₂ is increasing
* the error ‖eₖ‖₂ := ‖xᴸₖ - x*‖₂ is decreasing
* it is possible to transition cheaply from the LSLQ iterate to the LSQR iterate if there is an advantage (there always is in terms of error)
* if `A` is rank deficient, identify the minimum least-squares solution

#### Input arguments

* `A::AbstractLinearOperator`
* `b::Vector{Float64}`

#### Optional arguments

* `M::AbstractLinearOperator=opEye(size(A,1))`: a symmetric and positive definite dual preconditioner
* `N::AbstractLinearOperator=opEye(size(A,2))`: a symmetric and positive definite primal preconditioner
* `sqd::Bool=false` indicates whether or not we are solving a symmetric and quasi-definite augmented system
  If `sqd = true`, we solve the symmetric and quasi-definite system

      [ E   A' ] [ r ]   [ b ]
      [ A  -F  ] [ x ] = [ 0 ],

  where E = M⁻¹  and F = N⁻¹.

  If `sqd = false`, we solve the symmetric and indefinite system

      [ E   A' ] [ r ]   [ b ]
      [ A   0  ] [ x ] = [ 0 ].

  In this case, `N` can still be specified and indicates the norm in which `x` and the forward error should be measured.
* `λ::Float64=0.0` is a regularization parameter (see the problem statement above)
* `σ::Float64=0.0` is an underestimate of the smallest nonzero singular value of `A`---setting `σ` too large will result in an error in the course of the iterations
* `atol::Float64=1.0e-8` is a stopping tolerance based on the residual
* `btol::Float64=1.0e-8` is a stopping tolerance used to detect zero-residual problems
* `etol::Float64=1.0e-8` is a stopping tolerance based on the lower bound on the error
* `window::Int=5` is the number of iterations used to accumulate a lower bound on the error
* `utol::Float64=1.0e-8` is a stopping tolerance based on the upper bound on the error
* `itmax::Int=0` is the maximum number of iterations (0 means no imposed limit)
* `conlim::Float64=1.0e+8` is the limit on the estimated condition number of `A` beyond which the solution will be abandoned
* `verbose::Bool=false` determines verbosity.

#### Return values

`lslq()` returns the tuple `(x_lq, x_cg, err_lbnds, err_ubnds_lq, err_ubnds_cg, stats)` where

* `x_lq::Vector{Float64}` is the LQ solution estimate
* `x_cg::Vector{Float64}` is the CG solution estimate (i.e., the LSQR point)
* `err_lbnds::Vector{Float64}` is a vector of lower bounds on the LQ error---the vector is empty if `window` is set to zero
* `err_ubnds_lq::Vector{Float64}` is a vector of upper bounds on the LQ error---the vector is empty if `σ == 0` is left at zero
* `err_ubnds_cg::Vector{Float64}` is a vector of upper bounds on the CG error---the vector is empty if `σ == 0` is left at zero
* `stats::SimpleStats` collects other statistics on the run.

#### Stopping conditions

The iterations stop as soon as one of the following conditions holds true:

* the optimality residual is sufficiently small (`stats.status = "found approximate minimum least-squares solution"`) in the sense that either
  * ‖Aᵀr‖ / (‖A‖ ‖r‖) ≤ atol, or
  * 1 + ‖Aᵀr‖ / (‖A‖ ‖r‖) ≤ 1
* an approximate zero-residual solution has been found (`stats.status = "found approximate zero-residual solution"`) in the sense that either
  * ‖r‖ / ‖b‖ ≤ btol + atol ‖A‖ * ‖xᴸ‖ / ‖b‖, or
  * 1 + ‖r‖ / ‖b‖ ≤ 1
* the estimated condition number of `A` is too large in the sense that either
  * 1/cond(A) ≤ 1/conlim (`stats.status = "condition number exceeds tolerance"`), or
  * 1 + 1/cond(A) ≤ 1 (`stats.status = "condition number seems too large for this machine"`)
* the lower bound on the LQ forward error is less than etol * ‖xᴸ‖
* the upper bound on the CG forward error is less than utol * ‖xᶜ‖

#### References

* R. Estrin, D. Orban and M. A. Saunders, *Estimates of the 2-Norm Forward Error for SYMMLQ and CG*, Cahier du GERAD G-2016-70, GERAD, Montreal, 2016. DOI http://dx.doi.org/10.13140/RG.2.2.19581.77288.
* R. Estrin, D. Orban and M. A. Saunders, *LSLQ: An Iterative Method for Linear Least-Squares with an Error Minimization Property*, Cahier du GERAD G-2017-xx, GERAD, Montreal, 2017.
"""
function lslq(A :: AbstractLinearOperator, b :: AbstractVector{T};
              M :: AbstractLinearOperator=opEye(size(A,1)),
              N :: AbstractLinearOperator=opEye(size(A,2)),
              sqd :: Bool=false, λ :: Float64=0.0, σ :: Float64=0.0,
              atol :: Float64=1.0e-8, btol :: Float64=1.0e-8, etol :: Float64=1.0e-8,
              window :: Int=5, utol :: Float64=1.0e-8, itmax :: Int=0,
              conlim :: Float64=1.0e+8, verbose :: Bool=false) where T <: Number

  m, n = size(A)
  size(b, 1) == m || error("Inconsistent problem size")
  verbose && @printf("LSLQ: system of %d equations in %d variables\n", m, n)

  # If solving an SQD system, set regularization to 1.
  sqd && (λ = 1.0)
  λ² = λ * λ
  ctol = conlim > 0.0 ? 1/conlim : 0.0

  x_lq = zeros(T, n)    # LSLQ point
  x_cg = zeros(T, n)    # LSQR point
  err_lbnds = T[]
  err_ubnds_lq = T[]
  err_ubnds_cg = T[]

  # Initialize Golub-Kahan process.
  # β₁ M u₁ = b.
  Mu = copy(b)
  u = M * Mu
  β₁ = sqrt(@kdot(m, u, Mu))
  β₁ == 0.0 && return (x_lq, x_cg, err_lbnds, err_ubnds_lq, err_ubnds_cg,
                       SimpleStats(true, false, [0.0], [0.0], "x = 0 is a zero-residual solution"))
  β = β₁

  @kscal!(m, 1.0/β₁, u)
  @kscal!(m, 1.0/β₁, Mu)
  Nv = copy(A' * u)
  v = N * Nv
  α = sqrt(@kdot(n, v, Nv))  # = α₁

  # A'b = 0 so x = 0 is a minimum least-squares solution
  α == 0.0 && return (x_lq, x_cg, err_lbnds, err_ubnds_lq, err_ubnds_cg,
                      SimpleStats(true, false, [β₁], [0.0], "x = 0 is a minimum least-squares solution"))
  @kscal!(n, 1.0/α, v)
  @kscal!(n, 1.0/α, Nv)

  Anorm = α
  Anorm² = α * α

  # condition number estimate
  σmax = 0.0
  σmin = Inf
  Acond  = 0.0

  xlqNorm  = 0.0
  xlqNorm² = 0.0
  xcgNorm  = 0.0
  xcgNorm² = 0.0

  w = zeros(T, n)       # = w₀
  w̄ = copy(v)        # = w̄₁ = v₁

  err_lbnd = 0.0
  err_vec = zeros(T, window)

  # Initialize other constants.
  αL = α
  βL = β
  ρ̄ = -σ
  γ̄ = α
  ss = β₁
  c = -1.0
  s = 0.0
  δ = -1.0
  τ = α * β₁
  ζ = 0.0
  ζ̃  = T(0)
  csig = -1.0

  rNorm = β₁
  rNorms = [rNorm]
  ArNorm = α * β
  ArNorms = [ArNorm]

  verbose && @printf("%5s  %7s  %7s  %7s  %7s  %8s  %8s  %7s  %7s  %7s\n",
                     "Aprod", "‖r‖", "‖A'r‖", "β", "α", "cos", "sin", "‖A‖²", "κ(A)", "‖xL‖")
  verbose && @printf("%5d  %7.1e  %7.1e  %7.1e  %7.1e  %8.1e  %8.1e  %7.1e  %7.1e  %7.1e\n",
                     1, rNorm, ArNorm, β, α, c, s, Anorm², Acond, xlqNorm)

  iter = 0
  itmax == 0 && (itmax = m + n)

  status = "unknown"
  solved = solved_mach = solved_lim = (rNorm ≤ atol)
  tired  = iter ≥ itmax
  ill_cond = ill_cond_mach = ill_cond_lim = false
  zero_resid = zero_resid_mach = zero_resid_lim = false
  fwd_err_lbnd = false
  fwd_err_ubnd = false

  while ! (solved || tired || ill_cond)

    # Generate next Golub-Kahan vectors.
    # 1. βu = Av - αu
    @kscal!(m, -α, Mu)
    @kaxpy!(m, 1.0, A * v, Mu)
    u = M * Mu
    β = sqrt(@kdot(m, u, Mu))
    if β != 0.0
      @kscal!(m, 1.0/β, u)
      @kscal!(m, 1.0/β, Mu)

      # 2. αv = A'u - βv
      @kscal!(n, -β, Nv)
      @kaxpy!(n, 1.0, A' * u, Nv)
      v = N * Nv
      α = sqrt(@kdot(n, v, Nv))
      if α != 0.0
        @kscal!(n, 1.0/α, v)
        @kscal!(n, 1.0/α, Nv)
      end

      # rotate out regularization term if present
      αL = α
      βL = β
      if λ != 0.0
        (cL, sL, βL) = sym_givens(β, λ)
        αL = cL * α

        # the rotation updates the next regularization parameter
        λ = sqrt(λ² + (sL * α)^2)
      end
      Anorm² = Anorm² + αL * αL + βL * βL;  # = ‖Lₖ‖²
      Anorm = sqrt(Anorm²)
    end

    # Continue QR factorization of Bₖ
    #
    #       k   k+1     k   k+1      k  k+1
    # k   [ c'   s' ] [ γ̄      ] = [ γ   δ  ]
    # k+1 [ s'  -c' ] [ β   α⁺ ]   [     γ̄ ]
    (cp, sp, γ) = sym_givens(γ̄, βL)
    τ = -τ * δ / γ  # forward substitution for t
    δ = sp * αL
    γ̄ = -cp * αL

    if σ > 0.0
      # Continue QR factorization for error estimate
      μ̄ = -csig * γ
      (csig, ssig, ρ) = sym_givens(ρ̄, γ)
      ρ̄ = ssig * μ̄ + csig * σ
      μ̄ = -csig * δ

      # determine component of eigenvector and Gauss-Radau parameter
      h = δ * csig / ρ̄
      ω = sqrt(σ * (σ - δ * h))
      (csig, ssig, ρ) = sym_givens(ρ̄, δ)
      ρ̄ = ssig * μ̄ + csig * σ
    end

    # Continue LQ factorization of Rₖ
    ϵ̄ = -γ * c
    η = γ * s
    (c, s, ϵ) = sym_givens(ϵ̄, δ)

    # condition number estimate
    # the QLP factorization suggests that the diagonal of M̄ approximates
    # the singular values of B.
    σmax = max(σmax, ϵ, abs(ϵ̄))
    σmin = min(σmin, ϵ, abs(ϵ̄))
    Acond = σmax / σmin

    # forward substitution for z, ζ̄
    ζold = ζ
    ζ = (τ - ζ * η) / ϵ
    ζ̄ = ζ / c

    # residual norm estimate
    rNorm = sqrt((ss * cp - ζ * η)^2 + (ss * sp)^2)
    push!(rNorms, rNorm)

    ArNorm = sqrt((γ * ϵ * ζ)^2 + (δ * η * ζold)^2)
    push!(ArNorms, ArNorm)

    # compute LSQR point
    x_cg = x_lq + ζ̄ * w̄
    xcgNorm² = xlqNorm² + ζ̄ * ζ̄

    if σ > 0.0 && iter > 0
      err_ubnd_cg = sqrt(ζ̃ * ζ̃ - ζ̄  * ζ̄ )
      push!(err_ubnds_cg, err_ubnd_cg)
      fwd_err_ubnd = err_ubnd_cg ≤ utol * sqrt(xcgNorm²)
    end

    test1 = rNorm / β₁
    test2 = ArNorm / (Anorm * rNorm)
    test3 = 1 / Acond
    t1    = test1 / (1.0 + Anorm * xlqNorm / β₁)
    rtol  = btol + atol * Anorm * xlqNorm / β₁

    verbose && @printf("%5d  %7.1e  %7.1e  %7.1e  %7.1e  %8.1e  %8.1e  %7.1e  %7.1e  %7.1e\n",
                       1 + 2 * iter, rNorm, ArNorm, β, α, c, s, Anorm, Acond, xlqNorm)

    # update LSLQ point for next iteration
    w = c * w̄ + s * v
    w̄ = s * w̄ - c * v
    x_lq = x_lq + ζ * w
    xlqNorm² += ζ * ζ
    xlqNorm = sqrt(xlqNorm²)

    # check stopping condition based on forward error lower bound
    err_vec[mod(iter, window) + 1] = ζ
    if iter ≥ window
      err_lbnd = norm(err_vec)
      push!(err_lbnds, err_lbnd)
      fwd_err_lbnd = err_lbnd ≤ etol * xlqNorm
    end

    # compute LQ forward error upper bound
    if σ > 0.0
      η̃ = ω * s
      ϵ̃ = -ω * c
      τ̃ = -τ * δ / ω
      ζ̃ = (τ̃ - ζ * η̃) / ϵ̃
      push!(err_ubnds_lq, abs(ζ̃ ))
    end

    # Stopping conditions that do not depend on user input.
    # This is to guard against tolerances that are unreasonably small.
    ill_cond_mach = (1.0 + test3 ≤ 1.0)
    solved_mach = (1.0 + test2 ≤ 1.0)
    zero_resid_mach = (1.0 + t1 ≤ 1.0)

    # Stopping conditions based on user-provided tolerances.
    tired  = iter ≥ itmax
    ill_cond_lim = (test3 ≤ ctol)
    solved_lim = (test2 ≤ atol)
    zero_resid_lim = (test1 ≤ rtol)

    ill_cond = ill_cond_mach | ill_cond_lim
    solved = solved_mach | solved_lim | zero_resid_mach | zero_resid_lim | fwd_err_lbnd | fwd_err_ubnd

    iter = iter + 1
  end

  tired         && (status = "maximum number of iterations exceeded")
  ill_cond_mach && (status = "condition number seems too large for this machine")
  ill_cond_lim  && (status = "condition number exceeds tolerance")
  solved        && (status = "found approximate minimum least-squares solution")
  zero_resid    && (status = "found approximate zero-residual solution")
  fwd_err_lbnd  && (status = "forward error lower bound small enough")
  fwd_err_ubnd  && (status = "forward error upper bound small enough")

  stats = SimpleStats(solved, !zero_resid, rNorms, ArNorms, status)
  return (x_lq, x_cg, err_lbnds, err_ubnds_lq, err_ubnds_cg, stats)
end


text:
######### Stan program example  ###########

using StanSample, MCMCChains
using StatsPlots

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

observed_data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

# Default for tmpdir is to create a new tmpdir location
# To prevent recompilation of a Stan progam, choose a fixed location,
tmpdir=mktempdir()

sm = SampleModel("bernoulli", bernoullimodel,
  tmpdir=tmpdir);

rc = stan_sample(sm, data=observed_data);

if success(rc)
  chns = read_samples(sm; output_format=:mcmcchains)

  # Describe the results
  chns

  # Optionally, read samples as a a DataFrame
  df=read_samples(sm, output_format=:dataframe)
  first(df, 5)
  println()

  # Look at effective sample saize
  ess(chns) |> display
  println()

  # Check if StatsPlots is available and show basic MCMCChains plots
  if isdefined(Main, :StatsPlots)
    cd(@__DIR__) do
      p1 = plot(chns)
      savefig(p1, joinpath(tmpdir, "traceplot.pdf"))
      p2 = pooleddensity(chns)
      savefig(p2, joinpath(tmpdir, "pooleddensity.pdf"))
    end
  end

  df = read_summary(sm)
  df[df.parameters .== :theta, [:mean, :ess]]
end


text:
#############################################################################
# union.jl
# Handles unions of combinatorial sets.
#############################################################################

import Base: union

export union, get_elements

mutable struct UnionAtom <: CombiSet
  head::Symbol
  id_hash::UInt64
  children::Tuple
  elements::ValOrNothing
  baseset::AbstractArray
  cardinality::Int

  function UnionAtom(elements::ValOrNothing, children::Tuple, baseset::AbstractArray)
    this = new(:union, 0, children, elements, baseset, length(baseset))
    this.id_hash = object_id(this)
    return this
  end
end

function union(set::AllCombiSet...)
  sets = collect(set)
  baseset = Set([])
  n = length(sets)
  if n < 2
    error("Cannot compute the union of one set.")
  else
    if isa(sets, Array{Number})
      elements = []
      for i = 1:n
        elements = vcat(elements, sets[i])
      end
      elements = collect(Set(elements))
      return elements
    else
      for i = 1:n
        if isa(sets[i], CombiSet)
          baseset = union(baseset, Set(sets[i].baseset))
        else
          seti = Set(sets[i])
          baseset = union(baseset, seti)
        end
      end
      baseset = collect(baseset)
      newset = UnionAtom(nothing, set, baseset)
      return newset
    end
  end
end

function get_elements(x::UnionAtom)
  elements = []
  n = length(x.children)
  for i = 1:n
    if isa(x.children[i], CombiSet)
      elements = union(elements, Set(get_elements(x.children[i])))
    else
      seti = Set(x.children[i])
      elements = union(elements, seti)
    end
  end
  x.elements = collect(elements)
  return x.elements
end