text:
# bookkeeping data for a screen

# Screen is unique also in that its instantiation requires a window. It's because often we'd simply
# supply the rootwin that has already been created.
function newTwScreen( win::Ptr{Void} )
    obj = TwObj( TwScreenData(), Val{ :Screen } )
    begy, begx = getwinbegyx( win )
    maxy, maxx = getwinmaxyx( win )

    obj.window = win
    obj.height = maxy
    obj.width = maxx
    obj.xpos = begx
    obj.ypos = begy

    obj.hasFocus = true
    obj.acceptsFocus = true
    obj.isVisible = true
    #wbkgdset( win, COLOR_PAIR( 15 ) | @compat UInt( '.'))
    obj
end

function registerTwObj( scr::TwObj{TwScreenData}, o::TwObj )
    if o.screen.value != nothing
        if o.screen.value == scr
            return
        end
        unregisterTwObj( o.screen.value, o )
    end
    push!( scr.data.objects, o )
    o.screenIndex = length(scr.data.objects)
    o.screen = WeakRef( scr )
end

function unregisterTwObj( scr::TwObj{TwScreenData}, o::TwObj )
    if o.screen.value != scr
        throw( "unregister obj that doesn't belong to the screen")
    end
    idx = o.screenIndex
    deleteat!( scr.data.objects, idx )
    if o.panel != nothing
        hide_panel( o.panel )
    end
    for i in idx:length( scr.data.objects )
        scr.data.objects[i].screenIndex = i
    end
    if scr.data.focus >= idx
        scr.data.focus -= 1
    end
    if scr.data.focus == idx
        setTwFocusNext( scr )
    end
    o.screen = WeakRef()
end

function setTwFocusNext( scr::TwObj{TwScreenData})
    n = scr.data.focus
    fst = n
    result = nothing
    while true
        n +=1
        if n > length(scr.data.objects)
            n = 1
        end
        o = scr.data.objects[ n ]
        if o.acceptsFocus
            result = o
            break
        end
        if n == fst
            break
        end
    end
    if result != nothing
        scr.data.focus = n
    else
        scr.data.focus = 0
    end
    return result
end

function setTwFocusPrevious( scr::TwObj{TwScreenData} )
    n = scr.data.focus
    fst = n
    result = nothing
    while true
        n -=1
        if n <= 0
            n = length( scr.data.objects )
        end
        o = scr.data.objects[ n ]
        if o.acceptsFocus
            result = o
            break
        end
        if n == fst
            break
        end
    end
    if result != nothing
        scr.data.focus = n
    else
        scr.data.focus = 0
    end
    return result
end

function swapTwObjIndices( scr::TwObj{TwScreenData}, n1::Int, n2::Int )
    if n1!=n2 && 1 <= n1 <= length(scr.data.objects) && 1 <= n2 <= length( scr.data.objects )
        o1 = scr.data.objects[n1]
        o2 = scr.data.objects[n2]
        scr.data.objects[n1] = o2
        scr.data.objects[n2] = o1
        o1.screenIndex = n2
        o2.screenIndex = n1
        if scr.data.focus == n1
            scr.data.focus = n2
            o1.hasFocus = false
            o2.hasFocus = true
        elseif scr.data.focus == n2
            scr.data.focus = n1
            o2.hasFocus = false
            o1.hasFocus = true
        end
    end
end

function raiseTwObject( o::TwObj )
    log( "Raise: " * string( o ) )
    scr = o.screen.value
    if scr != nothing
        si = o.screenIndex
        deleteat!( scr.data.objects, si )
        push!( scr.data.objects, o )
        for i in si:length( scr.data.objects )
            scr.data.objects[i].screenIndex = i
            scr.data.objects[i].hasFocus = false
        end
        o.hasFocus = true
        scr.data.focus = length( scr.data.objects )
        refresh( o )
        top_panel( o.panel )
    end
end

function lowerTwObject( o::TwObj )
    log( "Lower: " * string( o ) )
    scr = o.screen.value
    if scr != nothing
        si = o.screenIndex
        deleteat!( scr.data.objects, si )
        unshift!( scr.data.objects, o )
        for i in 1:length( scr.data.objects )
            scr.data.objects[i].screenIndex = i
            scr.data.objects[i].hasFocus = false
        end
        scr.data.objects[end].hasFocus = true
        scr.data.focus = length( scr.data.objects )
        refresh( o )
        bottom_panel( o.panel )
    end
end

#=
Blocking call. it doesn't use inject directly, because
    * it needs to wait for populating of widgets
    * once there are widgets, get the topmost focused widget (TFW) and
    * inject tokens into the TFW
    * if TFW returns exit, get the next in line that accepts focus
    * otherwise, exit

This is some UGLY code. To be streamlined...
=#
function activateTwObj( scr::TwObj{TwScreenData}, tokens::Any=nothing )
    # consume one token, this also makes sure the readtoken function is jitted
    #=
    if tokens == nothing
        readtoken( scr.window )
    end
    =#
    refresh(scr) # clear any potential wait message
    focusObj = nothing
    focusIdx = scr.data.focus
    retvalue = nothing

    if focusIdx > 0 && focusIdx <= length( scr.data.objects )
        focusObj = scr.data.objects[ focusIdx ]
    end

    handleStatus = (st, tok)->begin
        focusIdx = scr.data.focus
        if focusIdx > 0 && focusIdx <= length( scr.data.objects )
            focusObj = scr.data.objects[ focusIdx ]
        else
            focusObj = nothing
        end
        if st== :exit_ok || st== :exit_nothing
            if st== :exit_ok
                retvalue = scr.value
            end
            if focusObj != nothing
                unregisterTwObj( scr, focusObj )
                # find the next focusObj, if not found, just return
                # otherwise, make it the new focus and continue
                focusObj = nothing
                for i in length( scr.data.objects ) :-1: 1
                    o = scr.data.objects[i]
                    if o.isVisible && o.acceptsFocus
                        focusObj = o
                        o.hasFocus = true
                        scr.data.focus = i
                        refresh( o )
                        break
                    end
                end
            end
            if focusObj == nothing
                return :really_exit
            end
        elseif st == :pass
            if focusObj != nothing
                i = length( scr.data.objects )
                while( st == :pass && i >= 1 )
                    o = scr.data.objects[i]
                    if o!=focusObj && o.isVisible && o.grabUnusedKey
                        st = inject( o, tok )
                    end
                    i -= 1
                end
                if st == :exit_ok || st == :exit_nothing
                    return :really_exit
                end
            end
        end
    end

    if tokens == nothing
        while true
            update_panels()
            doupdate()
            token = readtoken( scr.window )
            status = inject( scr, token )

            for o in scr.data.objects
                global twGlobProgressData
                if objtype(o) == :Progress && isready( twGlobProgressData.statusChannel )
                    inject( o, :progressupdate )
                end
            end

            if handleStatus( status, token ) == :really_exit
                return retvalue
            end
            sleep(0.01) # this is important or @async task won't run
        end
    else
        for token in tokens
            update_panels()
            doupdate()
            token = readtoken( scr.window )
            status = inject( scr, token )
            if handleStatus( status, token ) == :really_exit
                return retvalue
            end
        end
        # exhausted all the tokens, no obvious response
        return nothing
    end
end

function inject( scr::TwObj{TwScreenData}, token )
    global rootTwScreen
    result = :pass
    if token == :KEY_MOUSE
        (mstate, x,y,bs ) = getmouse()
    end
    if scr.data.focus != 0
        result = inject( scr.data.objects[ scr.data.focus], token )
        if result != :pass
            return result
        end
        if token == :F1
            h = helptext( scr.data.objects[ scr.data.focus ] )
            if h != ""
                helper = newTwViewer( rootTwScreen, h, posy= :center, posx=:center, showHelp=false, showLineInfo=false, bottomText = "Esc to continue" )
                activateTwObj( helper )
                unregisterTwObj( rootTwScreen, helper )
                return :got_it
            end
        elseif token == :tab
            lowerTwObject(scr.data.objects[ scr.data.focus ])
        elseif token == :shift_tab
            raiseTwObject(scr.data.objects[1])
        end
    end

    for i in length( scr.data.objects ) : -1 : 1
        o = scr.data.objects[i]
        o.hasFocus = (i==scr.data.focus)
        if token == :KEY_MOUSE && o.isVisible
            rely,relx = screen_to_relative( o.window, y, x )
            log( @sprintf( "Screen: Test mouse click on %s (r%d,c%d)", string( o ), y, x ) )
            if 0<=relx<o.width && 0<=rely<o.height
                log( "  raising it" )
                raiseTwObject( o )
                result = :got_it
                break
            end
        elseif o.isVisible && !o.hasFocus && o.grabUnusedKey
            result = inject( scr.data.objects[i], token )
            if result != :pass
                break
            end
        end
    end
    return result
end

function draw( scr::TwObj{TwScreenData} )
    focused = 0
    i = length( scr.data.objects )
    while i>=1
        o = scr.data.objects[i]
        if o.isVisible
            if o.hasFocus && focused < 1
                focused = i
            end
        end
        i-=1
    end
    scr.data.focus = focused
    for (i,o) in enumerate( scr.data.objects )
        o.hasFocus = (i == focused)
        if o.isVisible
            if panel_hidden( o.panel )
                show_panel( o.panel )
            end
            top_panel( o.panel )
            draw( o )
        else
            if !panel_hidden( o.panel )
                hide_panel( o.panel )
            end
        end
    end
end

function refresh( scr::TwObj{TwScreenData} )
    focused = 0
    i = length( scr.data.objects )
    while i>=1
        o = scr.data.objects[i]
        if o.isVisible
            if o.hasFocus && focused < 1
                focused = i
            end
        else
            erase( o )
        end
        i -= 1
    end
    werase( scr.window )
    scr.data.focus = focused
    for (i,o) in enumerate( scr.data.objects )
        o.hasFocus = (i == focused)
        if o.panel == nothing
            throw( string( o ) * " has unset panel. Check its constructor." )
        end
        if o.isVisible
            if panel_hidden( o.panel )
                show_panel( o.panel )
            end
            top_panel( o.panel )
            draw( o )
        else
            if !panel_hidden( o.panel )
                hide_panel( o.panel )
            end
        end
    end
end

function erase( scr::TwObj{TwScreenData} )
    werase( scr.window )
    for o in scr.data.objects
        erase( o )
    end
end


text:
@testset "MPS Description" for (Q, C, A) in [(eye(2),eye(3,2),eye(4,2)),(speye(2),speye(3,2),speye(4,2))]
  q₁      = ones(2)
  q₂      = 0.
  c₁      = zeros(3)
  c₂      = ones(3)
  b       = ones(4)
  lb      = zeros(2)
  ub      = ones(2)
  vars    = ["x₁", "x₂"]
  qpname  = "Test Problem"
  @testset "Domain Errors" begin
    @test_throws DomainError MPSDescription(eye(3,2), q₁, q₂, c₁, C, c₂, A, b, lb, ub, vars, qpname)
    @test_throws DomainError MPSDescription(Q, ones(3), q₂, c₁, C, c₂, A, b, lb, ub, vars, qpname)
    @test_throws DomainError MPSDescription(Q, q₁, q₂, c₁, eye(3), c₂, A, b, lb, ub, vars, qpname)
    @test_throws DomainError MPSDescription(Q, q₁, q₂, zeros(2), C, c₂, A, b, lb, ub, vars, qpname)
    @test_throws DomainError MPSDescription(Q, q₁, q₂, c₁, C, ones(2), A, b, lb, ub, vars, qpname)
    @test_throws DomainError MPSDescription(Q, q₁, q₂, c₁, C, c₂, eye(4,3), b, lb, ub, vars, qpname)
    @test_throws DomainError MPSDescription(Q, q₁, q₂, c₁, C, c₂, A, ones(3), lb, ub, vars, qpname)
    @test_throws DomainError MPSDescription(Q, q₁, q₂, c₁, C, c₂, A, b, zeros(3), ub, vars, qpname)
    @test_throws DomainError MPSDescription(Q, q₁, q₂, c₁, C, c₂, A, b, lb, ones(3), vars, qpname)
    @test_throws DomainError MPSDescription(Q, q₁, q₂, c₁, C, c₂, A, b, lb, ub, ["x₁", "x₂", "x₃"], qpname)
  end
  @testset "Methods for MPSDescription{$T}" for T in subtypes(AbstractFloat)
    problem = MPSDescription(T, Q, q₁, q₂, c₁, C, c₂, A, b, lb, ub, vars, qpname)

    @test name(problem) == qpname

    @test variables(problem) == Symbol.(vars)

    QQ, qq1, qq2 = objective(problem)
    @test QQ ≈ Q && eltype(QQ) == T
    @test qq1 ≈ q₁ && eltype(qq1) == T
    @test qq2 ≈ q₂ && eltype(qq2) == T

    AA, bb = equalities(problem)
    @test AA ≈ A && eltype(AA) == T
    @test bb ≈ b && eltype(bb) == T

    CC, cc1, cc2 = inequalities(problem)
    @test CC ≈ C && eltype(CC) == T
    @test cc1 ≈ c₁ && eltype(cc1) == T
    @test cc2 ≈ c₂ && eltype(cc2) == T

    llb, uub = bounds(problem)
    @test llb ≈ lb && eltype(llb) == T
    @test uub ≈ ub && eltype(uub) == T
  end
end

ext: jl


text:
import YaoBlocks: render_params, chsubblocks, subblocks,
        mat, dispatch!, niparams, getiparams, setiparams!,
        cache_key, print_block, apply!,
        ishermitian, isunitary, isreflexive,
        content, chcontent
import Base: ==, copy, hash

include("shortcuts.jl")
include("sequence.jl")
include("Diff.jl")
include("TrivilGate.jl")
include("Bag.jl")
include("Mod.jl")
include("ConditionBlock.jl")
include("RotBasis.jl")
include("reflect_gate.jl")
include("math_gate.jl")
include("pauli_strings.jl")
include("QFTBlock.jl")
include("FSimGate.jl")


text:
# Autogenerated wrapper script for Kaleido_jll for i686-w64-mingw32
export kaleido

JLLWrappers.@generate_wrapper_header("Kaleido")
JLLWrappers.@declare_executable_product(kaleido)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        kaleido,
        "bin\\kaleido.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Zygote, Optim

"""
    GRAPE(function, Array, Int, Int)

    Function using automatic differentiation to compute the gradient of a given
    functional

"""
function GRAPE(functional, controls; g_tol=1e-6)
    K, N = size(controls)
    # here we use an inplace gradient update
    function g!(G, x)
        _, back = Zygote.pullback(functional, reshape(x, (K, N)))
        G .= reshape(back(1)[1], size(G))
    end

    f = x -> functional(reshape(x, (K, N)))

    res = optimize(f, g!, controls, LBFGS(), Optim.Options(g_tol=g_tol, show_trace=true, store_trace = true))
    return res
end

function GRAPE_no_reshape(functional, controls; g_tol=1e-6)
    # here we use an inplace gradient update
    function g!(G, x)
		@views G .= Zygote.gradient(functional, x)[1]
    end
    res = optimize(functional, g!, controls, LBFGS(), Optim.Options(g_tol=g_tol, show_trace=true, store_trace = true))
    return res
end
"""
    CRAB
"""

function CRAB(functional, controls, n_SI)
    # N is the number of SI

end

function GROUP()
    println("use AD to optimise coefficients of dCRAB")
end

function Krotov()
    println("how do I work with AD?")
end


text:
export GenNeighbourhood
export GenNeighbourhoodIesimo

using Distributions
using Random
"""
    GenNeighbourhood(w::Vector{Real} ,µ::Real, var ::Real)
    Devuelve vecino del vector `w` de acorde a una normal
    de media `µ` y varianza `var`.
"""
function GenNeighbourhood(w::Vector{<:Real} ,µ::Real, var ::Real)::Vector{<:Real}
    distribution = Normal(µ, var)
    z=rand(distribution,size(w))
    w_new = map(
        x -> min(1, max(0,x)),
        w-z
    )
    return w_new
end

"""
GenNeighbourhoodIesimo(w::Vector{<:Real} ,µ::Real, var ::Real, index::Int)::Vector{<:Real}
    Devuelve vecino del vector `w` de acorde a una normal
    de media `µ` y varianza `var` donde solo se ha cambia la componente iéxima.
"""
function GenNeighbourhoodIesimo(w::Vector{<:Real} ,µ::Real, var ::Real, index::Int)::Vector{<:Real}
    distribution = Normal(µ, var)
    z=rand(distribution)
    w_new = w
    w_new[index] = min(1, max(0,w_new[index]-z))
    return w_new
end

"""
GenNeighbourhoodIesimo(w::Vector{<:Real} ,µ::Real, var ::Real, index::Int)::Vector{<:Real}
    Devuelve vecino del vector `w` de acorde a una normal
    de media `µ` y varianza `var` donde solo se ha cambia la componente iéxima.
"""
function GenNeighbourhoodIesimo(w::Vector{<:Real} ,µ::Real, var ::Real, n:: Int, coeficiente::Real )::Vector{<:Real}
    index = first(randperm(n), round(Int, coeficiente*n))
    distribution = Normal(µ, var)
    w_new = w
    for i in index
        w_new[i] = min(1, max(0,w_new[i]-rand(distribution)))
    end
    return w_new
end


text:
#math_graphics.jl:
#-------------------------------------------------------------------------------

#==Aliases
===============================================================================#
#Real data type used for display functionnality (low-level):
typealias DReal Float64

#==Data structures.
===============================================================================#
#TODO: Use library version?
abstract Point #So we can use as a function
immutable Point2D <: Point
	x::DReal
	y::DReal
end

#Could use BoundingBox, but cannot control type (and is technically not a BB).
immutable PExtents2D
	xmin::DReal
	xmax::DReal
	ymin::DReal
	ymax::DReal
end
PExtents2D(;xmin::Real=DNaN, xmax::Real=DNaN, ymin::Real=DNaN, ymax::Real=DNaN) =
	PExtents2D(xmin, xmax, ymin, ymax)

#Transform used for Canvas2D/CanvasF1 (only scale/offset - no rotation).
#xd = (xs + x0) * xu
#yd = (ys + y0) * yu
#where: xu/yu: user coordinates & xd/yd: device coordinates
#NOTE: Format used instead of xs*xu + x0 for improved numeric stability?
#TODO: Does immutable cause copy on function calls instead of using pointers??
immutable Transform2D
	xs::DReal
	x0::DReal
	ys::DReal
	y0::DReal
end

#==Constructor-like functions
===============================================================================#
#=Inputs
 -ext: Visible plot extents
 -inputb: Bounding box of display area.
NOTE:
 -Assumes device y-coordinate increases as we descend on the screen.
=#
function Transform2D(ext::PExtents2D, inputb::BoundingBox)
	xs = (inputb.xmax-inputb.xmin) / (ext.xmax-ext.xmin)
	ys = -(inputb.ymax-inputb.ymin) / (ext.ymax-ext.ymin)
	x0 = + inputb.xmin/xs - ext.xmin
	y0 = + inputb.ymax/ys - ext.ymin

	return Transform2D(xs, x0, ys, y0)
end

#==Basic operations
===============================================================================#
#Compute vector norm (assumes point refers to a vector):
function vecnorm(pt::Point2D)
	return sqrt(pt.x*pt.x+pt.y*pt.y)
end

#Simple scaling - no commutativity:
Base.(:*)(s::DReal, pt::Point2D) = Point2D(s*pt.x, s*pt.y)
Base.(:*)(s::Real, pt::Point2D) = DReal(s)*pt
Base.(:+)(p1::Point2D, p2::Point2D) = Point2D(p1.x+p2.x, p1.y+p2.y)

function union(e1::PExtents2D, e2::PExtents2D)
	return PExtents2D(
		min(e1.xmin, e2.xmin),
		max(e1.xmax, e2.xmax),
		min(e1.ymin, e2.ymin),
		max(e1.ymax, e2.ymax),
	)
end

#Overwrite with new extents, if defined
#-------------------------------------------------------------------------------
function Base.merge(base::PExtents2D, new::PExtents2D)
	#Pick maximum extents
	baseifnan(bv, nv) = isnan(nv)? bv: nv
	xmin = baseifnan(base.xmin, new.xmin)
	xmax = baseifnan(base.xmax, new.xmax)
	ymin = baseifnan(base.ymin, new.ymin)
	ymax = baseifnan(base.ymax, new.ymax)
	return PExtents2D(xmin, xmax, ymin, ymax)
end

Base.isfinite(ext::PExtents2D) =
	isfinite(ext.xmin) && isfinite(ext.xmax) &&
	isfinite(ext.ymin) && isfinite(ext.ymax)

#==Mapping/interpolation functions
===============================================================================#

#Apply transform that maps a data point to the canvas
#-------------------------------------------------------------------------------
function ptmap(xf::Transform2D, pt::Point2D)
	x = (pt.x + xf.x0)*xf.xs
	y = (pt.y + xf.y0)*xf.ys
	return Point2D(x, y)
end
function ptmap_rev(xf::Transform2D, pt::Point2D)
	x = pt.x/xf.xs - xf.x0
	y = pt.y/xf.ys - xf.y0
	return Point2D(x, y)
end

#Apply transform that maps a vector to the canvas (device) coordinates
#-------------------------------------------------------------------------------
function vecmap(xf::Transform2D, pt::Point2D)
	x = pt.x*xf.xs
	y = pt.y*xf.ys
	return Point2D(x, y)
end
function vecmap_rev(xf::Transform2D, pt::Point2D)
	x = pt.x/xf.xs
	y = pt.y/xf.ys
	return Point2D(x, y)
end

#Interpolate between two points.
#-------------------------------------------------------------------------------
function interpolate(p1::Point2D, p2::Point2D, x::DReal)
	m = (p2.y-p1.y) / (p2.x-p1.x)
	return m*(x-p1.x)+p1.y
end

#Last line


text:
using CSV, DataFrames, FeatherFiles

# Data from Gouda, The Netherlands

FILE_NAME = "./data/1_2008_Solar_Final_52.0 4.7.csv"
df1 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/2_2008_Solar_Final_52.0 4.7.csv"
df2 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/3_2008_Solar_Final_52.0 4.7.csv"
df3 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/4_2008_Solar_Final_52.0 4.7.csv"
df4 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/5_2008_Solar_Final_52.0 4.7.csv"
df5 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/6_2008_Solar_Final_52.0 4.7.csv"
df6 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/7_2008_Solar_Final_52.0 4.7.csv"
df7 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/8_2008_Solar_Final_52.0 4.7.csv"
df8 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/9_2008_Solar_Final_52.0 4.7.csv"
df9 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/10_2008_Solar_Final_52.0 4.7.csv"
df10 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/11_2008_Solar_Final_52.0 4.7.csv"
df11 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

FILE_NAME = "./data/12_2008_Solar_Final_52.0 4.7.csv"
df12 = CSV.File(FILE_NAME, normalizenames=true) |> DataFrame!

df = [df1;df2;df3;df4;df5;df6;df7;df8;df9;df10;df11;df12]
save("./data/solaryear.feather", df)


text:
"""
This file contains all the functions that are required to do basic mapping of
contour and contourf maps for climate data.

"""

function mollweide(clon::Float16)
    return subplot(projection=ccrs.Mollweide(central_longitude=clon));
end

function mollweide(gridbounds::Array)
    ax = plt.subplot(projection=ccrs.Mollweide()); ax.set_extent(gridbounds);
    return ax
end

function robinson(clon::Float16)
    ax = subplot(projection=ccrs.Robinson(central_longitude=clon));
    return ax
end

function robinson(gridbounds::Array)
    ax = subplot(projection=ccrs.Robinson()); ax.set_extent(gridbounds);
    return ax
end

function coastlines(ax)
    ax.coastlines()
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ cd6e8890-3632-11eb-0528-e1c9df81dc75
using Plots

# ╔═╡ e85c9b60-3632-11eb-22ba-9b4355ebac76
using PlutoUI

# ╔═╡ ff6eac80-3632-11eb-0cd1-5965f4ebeb11
using Unitful

# ╔═╡ f2500fd2-3632-11eb-14f1-8d0d9d355f8e
using OpticalPropagation

# ╔═╡ 154668e0-3633-11eb-169e-bd2dbc3c558d
gauss = MonoLightField2D(
    [exp(-(x^2+y^2)/3200) for x in -200:199, y in -200:199],
    wavelength=632.8u"nm",
    size=(1u"mm", 1u"mm")
);

# ╔═╡ a345c4f0-3634-11eb-17e3-f9096b32f28e
@bind d Slider(-20:20, default=0, show_value=true)

# ╔═╡ 34d8d670-3633-11eb-3f29-15e1181b491a
plot(plot(gauss), plot(fresnel2(gauss, d*u"cm")), size=(800,400))

# ╔═╡ Cell order:
# ╠═cd6e8890-3632-11eb-0528-e1c9df81dc75
# ╠═e85c9b60-3632-11eb-22ba-9b4355ebac76
# ╠═ff6eac80-3632-11eb-0cd1-5965f4ebeb11
# ╠═f2500fd2-3632-11eb-14f1-8d0d9d355f8e
# ╠═154668e0-3633-11eb-169e-bd2dbc3c558d
# ╠═a345c4f0-3634-11eb-17e3-f9096b32f28e
# ╠═34d8d670-3633-11eb-3f29-15e1181b491a


text:
macro d_str(str)
    # split the command by linebreaks
    s = split(str, r"\n|\r\n", limit=0, keepempty=false)
    # for each line, parse the string and construct the expression
    exs = Vector{Expr}()
    for i = 1:length(s)
        cmd = Douglass.parse(s[i])
        push!(exs, esc(Expr(:macrocall, Expr(:., :Douglass, QuoteNode(Symbol("@$(cmd.command)"))), @__LINE__,
            Douglass.active_df,
            cmd.by,
            cmd.sort,
            cmd.arguments,
            cmd.filter,
            cmd.use,
            cmd.options )) )
    end
    # For debugging:
    # @show exs
    # and return it to the REPL
    return Expr(:block, exs...)

end

macro use(t::Symbol)
    s = string(t)
    return esc(:( Douglass.set_active_df( Symbol($s) ) ))
end

function set_active_df(df::Symbol)
    global active_df
    active_df = df
end

# # the goal of this is just to return a string with the interpolated string
# macro m_str(str)
#     #esc(Meta.parse("\"$(escape_string(str))\""))
#     esc(Expr(:macrocall, Symbol("@d_exp_str"), @__LINE__, string(Meta.parse("\"$(escape_string(str))\"")) ) )
# end

# # experimental interface
# macro d_exp_str(str)
#     # parse the command
#     println("d_exp_str: String = $str")
#     cmd = Douglass.parse(str)
#     # return the call to the appropriate macro
#     # return esc(Expr(:macrocall, Expr(:., :Douglass, QuoteNode(Symbol("@$(cmd.command)"))), @__LINE__))
#     # return esc(Expr(:macrocall, Expr(:., :Douglass, QuoteNode(Symbol("@$(cmd.command)"))), @__LINE__, Douglass.active_df, cmd.arguments))
#     return esc(Expr(:macrocall, Expr(:., :Douglass, QuoteNode(Symbol("@$(cmd.command)"))), @__LINE__,
#         Douglass.active_df,
#         cmd.by,
#         cmd.sort,
#         cmd.arguments,
#         cmd.filter,
#         cmd.use,
#         cmd.options ))
# end


text:
using Base:Symbol
# utilities

_has_group(df::DataFrames.AbstractDataFrame, group::Any) = false
_has_group(df::DataFrames.AbstractDataFrame, group::Symbol) = hasproperty(df, group)
function _has_group(df::DataFrames.AbstractDataFrame, group::Vector{Symbol})
    all(x -> hasproperty(df, x), group)
end

_group_name(df::DataFrames.AbstractDataFrame, group::Symbol) = df[1, group]
function _group_name(df::DataFrames.AbstractDataFrame, groups::Vector{Symbol})
    join([df[1, g] for g in groups], ", ")
end

function _obtain_setindex_val(container::DataFrames.AbstractDataFrame, val::Vector{Symbol}, key::Symbol)
    if key == :dimensions
        out = []
        for v in val
            if hasproperty(container, v)
                push!(out, attr(label=v, values=container[!, v]))
            else
                push!(out, v)
            end
        end
        return out
    end
    val
end

# hook into setindex!
function _obtain_setindex_val(container::DataFrames.AbstractDataFrame, val::Symbol, key::Symbol)
    if hasproperty(container, val)
        if key == :dimensions
            return [attr(label=val, values=container[!, val])]
        end
        return container[!, val]
    end
    val
end

"""
$(SIGNATURES)
Build a trace of kind `kind`, using the columns of `df` where possible. In
particular for all keyword arguments, if the value of the keyword argument is a
Symbol and matches one of the column names of `df`, replace the value of the
keyword argument with the column of `df`

If `group` is passed and is a Symbol that is one of the column names of `df`,
then call `by(df, group)` and construct one trace per SubDataFrame, passing
all other keyword arguments. This means all keyword arguments are passed
applied to all traces

Also, when using this routine you can pass a function as a value for any
keyword argument. This function will be replaced by calling the function on the
DataFrame. For example, if I were to pass `name=(df) -> "Wage (average =
\$(mean(df[!, :X1])))"` then the `name` attribute on the trace would be replaced by
the  `Wage (average = XX)`, where `XX` is the average of the `X1` column in the
DataFrame.

The ability to pass functions as values for keyword arguments is particularly
useful when using the `group` keyword arugment, as the function will be applied
to each SubDataFrame. In the example above, the name attribute would set a
different mean for each group.
"""
function GenericTrace(df::DataFrames.AbstractDataFrame; group=missing, kind="scatter", kwargs...)
    if _has_group(df, group)
        _traces = []
        for dfg in collect(DataFrames.groupby(df, group))
            trace = GenericTrace(
                dfg;
                kind=kind,
                name=_group_name(dfg, group),
                legendgroup=_group_name(dfg, group),
                kwargs...
            )
            push!(_traces,  trace)
        end
        return GenericTrace[t for t in _traces]
    else
        if !ismissing(group)
            @warn "Unknown group $(group), skipping"
        end
    end

    attrs = attr()
    for (key, value) in kwargs
        # use `df` as `container` in setindex! for each property
        attrs[df, key] = value
    end
    if Symbol(kind) == :splom && :dimensions in keys(kwargs)
        attrs[:diagonal_visible] = false
        @assert attrs.dimensions isa Vector
        for d in attrs.dimensions
            d.axis_matches = true
        end
    end
    GenericTrace(kind; attrs...)
end

"""
$(SIGNATURES)
Pass the provided values of `x` and `y` as keyword arguments for constructing
the trace from `df`. See other method for more information
"""
function GenericTrace(df::DataFrames.AbstractDataFrame, x::Symbol, y::Symbol; kwargs...)
    GenericTrace(df; x=x, y=y, kwargs...)
end

"""
$(SIGNATURES)
Pass the provided value `y` as keyword argument for constructing the trace from
`df`. See other method for more information
"""
function GenericTrace(df::DataFrames.AbstractDataFrame, y::Symbol; kwargs...)
    GenericTrace(df; y=y, kwargs...)
end

"""
$(SIGNATURES)

lexographically sort obsvals to match orders, given keys into orders
"""
function _partial_val_sortperm(
        obsvals::Vector{<:Vector}, orders::Dict{Symbol}, order_keys::Vector{Symbol}
    )
    if length(obsvals) == 0 && length(order_keys) > 0
        error("obsvals is empty, but order_keys is not")
    end
    if length(obsvals) > 0 && (length(obsvals[1]) !== length(order_keys))
        error("obsvals[i] and order_keys must be same length")
    end

    to_sort = deepcopy(obsvals)
    for trace_i in 1:length(obsvals)
        for (group_i, colname) in enumerate(order_keys)
            if haskey(orders, Symbol(colname))
                want_order = orders[Symbol(colname)]
                group_val = obsvals[trace_i][group_i]
                to_sort[trace_i][group_i] = findfirst(isequal(group_val), want_order)
            end
        end
    end
    return sortperm(to_sort)
end

@inline function _ind2sub_rowmajor(nr::Int, nc::Int, ix::Int)
    out = CartesianIndices((nc, nr))[ix].I
    CartesianIndex(out[2], out[1])
end

function _faceted_subplot_titles(nr::Int, nc::Int, titles::Vector{String})
    subplot_titles = _Maybe{String}[missing for _ in 1:nr, _ in 1:nc]

    for title_ix in 1:length(titles)
        idx = _ind2sub_rowmajor(nr, nc, title_ix)
        subplot_titles[idx] = titles[title_ix]
    end

    # we need to transpose b/c column vs row major...
    permutedims(subplot_titles, (2, 1))
end

"""
$(SIGNATURES)
Construct a plot using the columns of `df` if possible. For each keyword
argument, if the value of the argument is a Symbol and the `df` has a column
whose name matches the value, replace the value with the column of the `df`.

If `group` is passed and is a Symbol that is one of the column names of `df`,
then call `by(df, group)` and construct one trace per SubDataFrame, passing
all other keyword arguments. This means all keyword arguments are passed
applied to all traces

All keyword arguments are passed directly to the trace constructor, with the
following exceptions:

- `category_orders`: A mapping from DataFrame column name to a vector of values that
  appear in that column. Wherever the column is used in the chart (facets, x, y, etc.)
  the chart elements tied to the column will be sorted according to the vector of values
  specified here. This ordering can be partial, in which case all elements that appear
  in the DataFrame, but not in the vector of values specified here, will be included
  *after* values explicitly listed here.
- `labels`: A mapping from column names (symbols) to text that should appear on the chart.
  This can be used to customize the title of axes, title of legend, etc.
- `facet_row`: A symbol denoting which DataFrame column should be used to create
  facetted subplots in the vertical direction.
- `facet_col`: A symbol denoting which DataFrame column should be used to create
  facetted subplots in the horizontal direction.
- `symbol`: A symbol noting which column should be used for the marker symbol. The values
  of this column are not used directly, but rather a new trace will be created for each
  unique value in this column. Each trace will have a different marker symbol, such
  as circle, diamond, square, etc.
- `color`: Similar to `symbol`, but for the color of a trace
- `line_dash`: Similar to `symbol`, but for the line style of a trace.
"""
function Plot(
        df::DataFrames.AbstractDataFrame,
        l::Layout=Layout();
        category_orders::_Maybe{Union{PlotlyAttribute,<:Dict{Symbol}}}=missing,
        labels::_Maybe{Union{PlotlyAttribute,<: Dict}}=missing,
        facet_row::_Maybe{Symbol}=missing,
        facet_row_wrap::_Maybe{Int}=missing,
        facet_col::_Maybe{Symbol}=missing,
        facet_col_wrap::_Maybe{Int}=missing,
        group::_Maybe{Union{Symbol,Vector{Symbol}}}=missing,
        symbol::_Maybe{Symbol}=missing,
        color::_Maybe{Symbol}=missing,
        line_dash::_Maybe{Symbol}=missing,
        kw...
    )
    groupby_cols = Symbol[]
    facet_cols = Symbol[]
    kwargs = Dict(kw)
    kw_keys = keys(kwargs)

    # prep groupings (symbol, color)
    _group_args = (;symbol, color, line_dash)
    _group_maps = Dict{Symbol,Dict{Any,String}}()
    _group_cols = Dict{Symbol,Symbol}()
    for (_group_arg, _group_col) in pairs(_group_args)
        if !ismissing(_group_col)
            if !_has_group(df, _group_col)
                error("DataFrame is missing $(_group_col) column, cannot set as $(_group_arg)")
            else
                # mark that we have this group and add to list of groupby_cols
                push!(groupby_cols, _group_col)

                # make note of the column that we will use for the groupby
                _group_cols[_group_arg] = _group_col

                # prepare map from group observation in data to the value in trace
                _group_map = Dict{Any,String}()
                defaults = _get_default_seq(l, _group_arg)
                for (i, val) in enumerate(unique(df[!, _group_col]))
                    _group_map[val] = defaults[i]
                end
                _group_maps[_group_arg] = _group_map
            end
        end
    end
    group_attr_pairs = [(:symbol, :marker_symbol), (:color, :marker_color), (:line_dash, :line_dash)]

    if length(groupby_cols) > 0 && !ismissing(group)
        @warn "One of color or symbol present AND group -- group will be ignored"
    end

    # handle labels
    label_map = Dict{Symbol,Any}(pairs((;facet_row, facet_col, color, symbol, line_dash)))
    for v in values(label_map)
        !ismissing(v) && setindex!(label_map, v, v)
    end
    if !ismissing(labels)
        for (k, v) in pairs(_symbol_dict(labels))
            label_map[k] = v
        end
    end

    # axis titles
    for ax in [:x, :y, :z]
        if ax in kw_keys
            ax_val = kwargs[ax]
            if !(ax_val isa Symbol)
                continue
            end
            # use `get` below because if something was passed into labels it
            # will already  be set in `label_map` and we should use the Given
            # label. Otherwise use hte value directly.
            # e.g.  in plot(x=:total_bill, labels=Dict(:total_bill=>"Total Bill"))
            # label_map[:total_bill] will be "Total Bill", so we should set
            # label_map[:x] = "Total Bill".
            label_map[ax] = get(label_map, ax_val, ax_val)
            label_map[ax_val] = ax_val
        end
    end

    # Handle category orders
    orders = Dict{Symbol,Any}()
    if !ismissing(category_orders)
        for (k, v) in pairs(category_orders)
            if !_has_group(df, k)
                @warn "Unknown category $k (dataframe has no column $k). Skipping"
                continue
            end
            unique_vals = unique(df[!, k])
            cat_order = deepcopy(v)
            for v in unique_vals
                if !(v in cat_order)
                    push!(cat_order, v)
                end
            end
            orders[k] = cat_order
        end
    end

    # handle facets
    rows = []
    cols = []
    if !ismissing(facet_row)
        if !_has_group(df, facet_row)
            error("DataFrame is missing $(facet_row), cannot set as facet_row")
        else
            rows = get(orders, facet_row, unique(df[:, facet_row]))
        end
    end

    if !ismissing(facet_col)
        if !_has_group(df, facet_col)
            error("DataFrame is missing $(facet_col), cannot set as facet_col")
        else
            cols = get(orders, facet_col, unique(df[:, facet_col]))
        end
    end

    Nrows = max(1, length(rows))
    Ncols = max(1, length(cols))

    # handle facets here
    if !ismissing(facet_row) && !ismissing(facet_col_wrap)
        @warn "cannot set facet_row and facet_col_wrap -- setting facet_col_wrap to missing"
        facet_col_wrap = missing
    end

    if !ismissing(facet_col) && !ismissing(facet_row_wrap)
        @warn "cannot set facet_col and facet_row_wrap -- setting facet_row_wrap to missing"
        facet_row_wrap = missing
    end

    out_layout = l
    if Nrows > 1 || Ncols > 1
        subplot_kw = Dict{Symbol,Any}()
        if Ncols > 1

            if !ismissing(facet_col_wrap) && Ncols > facet_col_wrap
                nr = ceil(Int, Ncols / facet_col_wrap)
                nc = facet_col_wrap
                subplot_kw[:cols] = nc
                subplot_kw[:rows] = nr
                subplot_kw[:subplot_titles] = _faceted_subplot_titles(
                    nr, nc, map(x -> "$(label_map[facet_col])=$(x)", cols)
                )
                subplot_kw[:shared_yaxes] = "all"
            else
                subplot_kw[:cols] = Ncols
                subplot_kw[:column_titles] = map(x -> "$(label_map[facet_col])=$(x)", cols)
                subplot_kw[:shared_yaxes] = true
            end

            push!(facet_cols, facet_col)
        end

        if Nrows > 1
            if !ismissing(facet_row_wrap) && Nrows > facet_row_wrap
                nr, nc = facet_row_wrap, ceil(Int, Nrows / facet_row_wrap)
                subplot_kw[:rows] = nr
                subplot_kw[:cols] = nc
                subplot_kw[:subplot_titles] = _faceted_subplot_titles(
                    nr, nc,  map(x -> "$(label_map[facet_row])=$(x)", rows)
                )
                subplot_kw[:shared_xaxes] = "all"
            else
                subplot_kw[:rows] = Nrows
                subplot_kw[:row_titles] = map(x -> "$(label_map[facet_row])=$(x)", rows)
                subplot_kw[:shared_xaxes] = true
            end

            push!(facet_cols, facet_row)
        end

        subplots = Subplots(; subplot_kw...)
        out_layout = Layout(subplots; out_layout.fields...)
    end

    out = Plot(out_layout)

    legend_order = Vector{<:Any}[]
    for dfg in collect(DataFrames.groupby(df, facet_cols))
        row_ix = Nrows == 1 ? 1 : findfirst(isequal(dfg[1, facet_row]), rows)
        col_ix = Ncols == 1 ? 1 : findfirst(isequal(dfg[1, facet_col]), cols)

        # @show row_ix, col_ix, :before, dfg[1, facet_col]
        if !ismissing(facet_col_wrap) || !ismissing(facet_row_wrap)
            ## ind2sub
            # when making facets with wrap, we need to convert from the
            # non-wrapped row_ix (=1 when facet_col_wrap) and col_ix (potentially
            # greater than facet_col_wrap) to the wrapped row_ix and col_ix
            # we do this by converting from a linear index in the faceted + wrapped
            # dimenmsion only into a cartesian index in both dimensions
            _facet_ix = ismissing(facet_col_wrap) ? row_ix : col_ix
            row_ix, col_ix = _ind2sub_rowmajor(size(out.layout.subplots.grid_ref)..., _facet_ix).I
        end

        for sub_dfg in collect(DataFrames.groupby(dfg, groupby_cols))
            all_kw = attr(;kwargs...)
            group_obs = Any[sub_dfg[1, g] for g in groupby_cols]
            group_name = join(string.(group_obs), ",")
            if length(group_name) > 0
                all_kw.name = group_name
                all_kw.legendgroup = group_name
                all_kw.showlegend = false
            end
            for (groupname, trace_attr) in group_attr_pairs
                if !ismissing(getfield(_group_args, groupname))
                    obs = sub_dfg[1, _group_cols[groupname]]
                    all_kw[trace_attr] = _group_maps[groupname][obs]
                end
            end
            traces = GenericTrace(sub_dfg; group=group, all_kw...,)

            traces_vec = traces isa GenericTrace ? [traces] : traces
            for trace in traces_vec
                push!(legend_order, group_obs)
                add_trace!(out, trace, row=row_ix, col=col_ix)
            end
        end
    end

    # turn on legend once for each
    if length(groupby_cols) > 0
        group_labels = getindex.(Ref(label_map), groupby_cols)
        out.layout[:legend_title_text] = join(string.(group_labels), ", ")
        seen = Set()
        for trace in out.data
            grp = trace.legendgroup
            if !(grp in seen)
                push!(seen, grp)
                trace.showlegend = true
            end
        end
    end
    out.layout.legend_tracegroupgap = 0

    # fix order of traces according to category_orders
    trace_order = _partial_val_sortperm(legend_order, orders, groupby_cols)
    out.data = out.data[trace_order]

    # add x axis title to bottom row of subplots and y axis title to left column
    for gr in out.layout.subplots.grid_ref[:, 1]
        yname = gr[1].layout_keys[2]
        if :y in keys(label_map)
            setifempty!(out.layout, Symbol("$(yname)_title_text"), label_map[:y])

            if haskey(orders, kwargs[:y])
                # order the yticks
                setifempty!(out.layout, Symbol("$(yname)_categoryorder"), "array")
                setifempty!(out.layout, Symbol("$(yname)_categoryarray"), orders[kwargs[:y]])
            end
        end
    end
    for gr in out.layout.subplots.grid_ref[end, :]
        xname = gr[1].layout_keys[1]
        if :x in keys(label_map)
            setifempty!(out.layout, Symbol("$(xname)_title_text"), label_map[:x])

            if haskey(orders, kwargs[:x])
                # order the xticks
                setifempty!(out.layout, Symbol("$(xname)_categoryorder"), "array")
                setifempty!(out.layout, Symbol("$(xname)_categoryarray"), orders[kwargs[:x]])
            end
        end
    end

    if Symbol(get(kwargs, :kind, "")) == :splom
        out.layout.dragmode = "select"
    end

    return out
end

"""
$(SIGNATURES)
Construct a plot from `df`, passing the provided values of x and y as keyword
arguments. See docstring for other method for more information.
"""
function Plot(d::DataFrames.AbstractDataFrame, x::Symbol, y::Symbol, l::Layout=Layout();
              kwargs...)
    Plot(d, l; x=x, y=y, kwargs...)
end

"""
$(SIGNATURES)
Construct a plot from `df`, passing the provided value y as a keyword argument.
See docstring for other method for more information.
"""
function Plot(d::DataFrames.AbstractDataFrame, y::Symbol, l::Layout=Layout();
              kwargs...)
    Plot(d, l; y=y, kwargs...)
end

for t in _TRACE_TYPES
    str_t = string(t)
    @eval $t(df::DataFrames.AbstractDataFrame; kwargs...) = GenericTrace(df; kind=$(str_t), kwargs...)
end


text:
# conversions and promotions from integer and rational types to algebraic number types
import Base.convert
import Base.promote_rule

# Algebraic number from integer
#convert{T<:Integer}(::Type{AlgebraicNumber}, x::T) =
#    AlgebraicNumber(BigInt[-x,one(T)], Complex{BigFloat}(x))
# Algebraic number from rational
#convert{T<:Integer}(::Type{AlgebraicNumber}, x::Rational{T}) =
#    AlgebraicNumber(BigInt[-num(x), den(x)], Complex{BigFloat}(x))

# promotions
promote_rule(x::Type{T},          y::Type{AlgebraicNumber{S,F}}) where {T<:Integer,S,F} = AlgebraicNumber
promote_rule(x::Type{Rational{T}},y::Type{AlgebraicNumber{S,F}}) where {T<:Integer,S,F} = AlgebraicNumber

# conversions back
function convert(::Type{Int64},an::AlgebraicNumber)
	c = an.coeff
	if length(c)==2 && abs(c[2])==1
		return convert(Int64, -c[1]*c[2])
	else
		throw(InexactError())
	end
end


text:
numTests = 5
maxDepth = 50
numColors = 6
target = [100.0, 100.0, 100.0]
fnEvals=0
# alphaScoreMatrix = 10 .* rand(numColors+1,numColors+1) .- 5
# betaScoreMatrix = 10 .* rand(numColors+1,numColors+1) .- 5
# gammaScoreMatrix = 10 .* rand(numColors+1,numColors+1) .- 5

using LinearAlgebra
using Plotly
include("EvaluationFunctions.jl")
include("GenerationFunctions.jl")
#generate = GenerateRandomly
generate = ExhaustiveHillClimbing

scores = Array{Float64, 2}(undef, numTests, maxDepth)
for i=1:numTests

    global    alphaScoreMatrix = 10 .* rand(numColors+1,numColors+1) .- 5
    global    betaScoreMatrix = 10 .* rand(numColors+1,numColors+1) .- 5
    global    gammaScoreMatrix = 10 .* rand(numColors+1,numColors+1) .- 5
    testTree1 = generate(1)
    scores[i,1]= GetCombinedScore(testTree1)
    for j=2:maxDepth
        testTree1 = generate(1, cand=testTree1)
        scores[i,j]= GetCombinedScore(testTree1)
    end
end

plot(transpose(scores))

max_stars_repo_name: DesignEngrLab/ColorGraph

max_issues_repo_name: DesignEngrLab/ColorGraph

max_forks_repo_name: DesignEngrLab/ColorGraph


text:
include("../imports.jl")

@testset "keypointpreprocessing.jl" begin

    ks = [
        SVector{2, Float32}(10, 10),
        SVector{2, Float32}(50, 80),
    ]
    sz = (100, 100)

    block = Keypoints{2}(2)
    enc = KeypointPreprocessing(sz)
    ctx = Training()

    testencoding(enc, block, ks)
    y = encode(enc, ctx, block, ks)
    ks_ = decode(enc, ctx, encodedblock(enc, block), y)
    @test ks ≈ ks_

end


text:
function ueberpruefeLaengeStrahl(A,Strahl,Laenge1, Laenge2, Laenge3,Laenge4, Laenge5, Laenge6,Laenge7,Laenge8,Laenge9)

if round(A[Strahl,1];digits=1) == round(Laenge1;digits=1)
    if round(A[Strahl,2];digits=1) == round(Laenge2;digits=1)
        if round(A[Strahl,3];digits=1) == round(Laenge3;digits=1)
            if round(A[Strahl,4];digits=1) == round(Laenge4;digits=1)
                if round(A[Strahl,5];digits=1) == round(Laenge5;digits=1)
                    if round(A[Strahl,6];digits=1) == round(Laenge6;digits=1)
                        if round(A[Strahl,7];digits=1) == round(Laenge7;digits=1)
                           if round(A[Strahl,8];digits=1) == round(Laenge8;digits=1)
                               if round(A[Strahl,9];digits=1) == round(Laenge9;digits=1)
                                    println("Ihr habt die Zeile der Matrix korrekt interpretiert.")
                               else
                                    println("Überprüft noch einmal den Eintrag für Laenge9.")
                               end
                           else
                                println("Überprüft noch einmal den Eintrag für Laenge8.")
                           end
                        else
                            println("Überprüft noch einmal den Eintrag für Laenge7.")
                        end
                    else
                        println("Überprüft noch einmal den Eintrag für Laenge6.")
                    end
                else
                    println("Überprüft noch einmal den Eintrag für Laenge5.")
                end
            else
                println("Überprüft noch einmal den Eintrag für Laenge4.")
            end
        else
            println("Überprüft noch einmal den Eintrag für Laenge3.")
        end
    else
         println("Überprüft noch einmal den Eintrag für Laenge2.")
    end
else
    println("Überprüft noch einmal den Eintrag für Laenge1.")
end

end


text:
# This file is a part of LegendHDF5IO.jl, licensed under the MIT License (MIT).

__precompile__(true)

module LegendHDF5IO

using ArraysOfArrays
using ElasticArrays
using EncodedArrays
using LegendDataTypes
using RadiationDetectorSignals
using StaticArrays
using StructArrays
using Tables
using Unitful

import HDF5
import TypedTables

using LegendDataTypes: readdata, writedata, getunits, setunits!,
    units_from_string, units_to_string,
    read_from_properties, write_to_properties!

using RadiationDetectorSignals: RealQuantity, ArrayOfDims, AosAOfDims, SArrayOfDims,
    recursive_ndims

include("generic_io.jl")
include("radsig_io.jl")
include("geant4_hdf5.jl")

const _datatype_dict = Dict{String,Type}()

function __init__()
    _datatype_dict[datatype_to_string(EventType)] = EventType
    _datatype_dict["table{t0,dt,values}"] = Vector{<:RDWaveform}
end

end # module


text:
function filter_CIR_keep_fixed_number(δ, γ, σ, λ, data, fixed_number::Int64; silence = false)
    # println("filter_WF_mem2")

    function prune_keeping_fixed_number(Λ_of_t, wms_of_t)
        Λ_of_t_kept, wms_of_t_kept = keep_fixed_number_of_weights(Λ_of_t, wms_of_t, fixed_number)
        return Λ_of_t_kept, normalise(wms_of_t_kept)
    end

    filter_CIR_pruning(δ, γ, σ, λ, data, prune_keeping_fixed_number; silence = silence)

end

function filter_CIR_keep_above_threshold(δ, γ, σ, λ, data, ε::Float64; silence = false)
    # println("filter_WF_mem2")

    function prune_keeping_above_threshold(Λ_of_t, wms_of_t)
        Λ_of_t_kept, wms_of_t_kept = keep_above_threshold(Λ_of_t, wms_of_t, ε)
        return Λ_of_t_kept, normalise(wms_of_t_kept)
    end

    filter_CIR_pruning(δ, γ, σ, λ, data, prune_keeping_above_threshold; silence = silence)

end

function filter_CIR_fixed_fraction(δ, γ, σ, λ, data, fraction::Float64; silence = false)
    # println("filter_WF_mem2")
    function prune_keeping_fixed_fraction(Λ_of_t, wms_of_t)
        Λ_of_t_kept, wms_of_t_kept = keep_fixed_fraction(Λ_of_t, wms_of_t, fraction)
        return Λ_of_t_kept, normalise(wms_of_t_kept)
    end

    filter_CIR_pruning(δ, γ, σ, λ, data, prune_keeping_fixed_fraction; silence = silence)

end

function filter_CIR_pruning(δ, γ, σ, λ, data, do_the_pruning::Function; silence = false, trim0 = true)

    times = keys(data) |> collect |> sort
    Λ_of_t = Dict{Float64, Array{Int64,1}}()
    wms_of_t = Dict{Float64, Array{Float64,1}}()
    θ_of_t = Dict{Float64, Float64}()

    filtered_θ, filtered_Λ, filtered_wms = update_CIR_params([1.], δ, γ/σ^2, λ, [0], data[0])
    pruned_Λ, pruned_wms = do_the_pruning(filtered_Λ, filtered_wms)

    Λ_of_t[0] = filtered_Λ
    wms_of_t[0] = filtered_wms # = 1.
    θ_of_t[0] = filtered_θ

    for k in 1:(length(times)-1)
        if (!silence)
            println("Step index: $k")
            println("Number of components before pruning : $(length(filtered_Λ))")
            println("Number of components after pruning: $(length(pruned_Λ))")
        end
        filtered_θ, filtered_Λ, filtered_wms = get_next_filtering_distribution(pruned_Λ, pruned_wms, filtered_θ, times[k], times[k+1], δ, γ, σ, λ, data[times[k+1]])
        pruned_Λ, pruned_wms = do_the_pruning(filtered_Λ, filtered_wms)
        if trim0
            pruned_Λ, pruned_wms = pruned_Λ[pruned_wms.>0], pruned_wms[pruned_wms.>0]
        end
        Λ_of_t[times[k+1]] = filtered_Λ
        wms_of_t[times[k+1]] = filtered_wms
        θ_of_t[times[k+1]] = filtered_θ
    end

    return Λ_of_t, wms_of_t, θ_of_t

end

function filter_CIR_pruning_logweights(δ, γ, σ, λ, data, do_the_pruning::Function; silence = false)
    Λ_of_t, wms_of_t, θ_of_t = filter_CIR_pruning(δ, γ, σ, λ, data, do_the_pruning::Function; silence = silence)
    logwms_of_t = wms_of_t |> convert_weights_to_logweights
    return Λ_of_t, logwms_of_t, θ_of_t
end

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
@testset "Structural model tests" begin
    @testset "Constant signal with basic structural model" begin
        y = ones(30)
        model = structural(y, 2)

        @test isa(model, StateSpaceModel)
        @test model.mode == "time-invariant"

        ss = statespace(model)

        @test ss.filter_type <: KalmanFilter
        @test isa(ss, StateSpace)

        @test all(ss.model.H .< 1e-6)
        @test all(ss.model.Q .< 1e-6)
        compare_forecast_simulation(ss, 20, 1000, 1e-3)
    end

    @testset "Constant signal with exogenous variables" begin
        y = ones(15)
        X = randn(30, 2)

        model = structural(y, 2; X = X)

        @test isa(model, StateSpaceModel)
        @test model.mode == "time-variant"

        ss = statespace(model)

        @test ss.filter_type <: KalmanFilter
        @test isa(ss, StateSpace)
        @test all(ss.model.H .< 1e-6)
        @test all(ss.model.Q .< 1e-6)

        compare_forecast_simulation(ss, 10, 1000, 1e-3)
    end

    @testset "Error tests" begin
        y = ones(15, 1)
        @test_throws ErrorException structural(y, 2; X = ones(10, 2))

        Z = Array{Float64, 3}(undef, 1, 2, 20)
        for t = 1:20
            Z[:, :, t] = [1 0]
        end
        T = [1. 1; 0 1]
        R = [1. 0; 0 1]

        model = StateSpaceModel(y, Z, T, R)
        ss = statespace(model)
        @test_throws ErrorException sim = simulate(ss, 10, 1000)
    end

    @testset "Model with 32-bit precision" begin
        model = structural(ones(Float32, 50), 3)
        @test size(model) == (50, 1, 4, 3)
    end
end


text:
# This file was generated, do not modify it. # hide
Dists.logpdf(Dists.Normal(), 1.0)

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
# FYI: TagBot (github.com/apps/julia-tagbot) depends on the "Repository", "Version",
# "Commit", and "Release notes" fields. If you're going to change the format here,
# please ping @christopher-dG or open an issue on TagBot!
function pull_request_contents(;
    registration_type::AbstractString,
    package::AbstractString,
    repo::AbstractString,
    user::AbstractString,
    version::VersionNumber,
    commit::AbstractString,
    release_notes::AbstractString,
    gitref::AbstractString="",
    reviewer::AbstractString="",
    reference::AbstractString="",
    meta::AbstractString="",
)
    title = if isempty(registration_type)
        "Registering: $package v$version"
    else
        "$registration_type: $package v$version"
    end

    # Build the PR body one line at a time.
    lines = [
        "- Registering package: $package",
        "- Repository: $repo",
        "- Created by: $user",
        "- Version: v$version",
        "- Commit: $commit",
    ]

    isempty(gitref) || push!(lines, "- Git reference: $gitref")
    isempty(reviewer) || push!(lines, "- Reviewed by: $reviewer")
    isempty(reference) || push!(lines, "- Reference: $reference")
    isempty(release_notes) || push!(
        lines,
        "- Release notes:",
        "<!-- BEGIN RELEASE NOTES -->",
        join(map(line -> "> $line", split(release_notes, "\n")), "\n"),
        "<!-- END RELEASE NOTES -->",
        ""
    )
    isempty(meta) || push!(lines, meta)

    return title, strip(join(lines, "\n"))
end


text:
"""
Magnetic levitation system
q = (z), u = (z1dot), extra_state = i
g(q) = (q-0.05, 0.20-q)
W = del g / del q = [1 -1]
M = m*Matrix(1.0f0*I, 1, 1)
h = m*g*ones(Float32, 1)
"""

using LinearAlgebra
using PyPlot
using TimeStepping

par = Dict{Symbol, Float32}(
    :m=>0.01187, :g=>9.81, :C=>1.4e-4, :L1=>0.65, :R=>28.7, :qd=>0.1, :tf=>4.0,
    :Δt=>1e-3, :ε=>0.1
)

δ = convert.(Float32, [0.05, 0.20])       # Hard stops of the maglev system.
q0 = convert.(Float32, [δ[1]])            # Initial ball position
u0 = convert.(Float32, [0.0])             # Initial ball velocity
extra_state0 = convert.(Float32, [0.0])   # Initial magnet current    (5.784 to lift-off)

function gap(q::AbstractArray, u::AbstractArray)
    return [q[1].-δ[1], δ[2].-q[1]], convert.(Float32, [1.0 -1])
end

function dynamics(q::AbstractArray, u::AbstractArray, extra_state::AbstractArray)
    return par[:m]*Matrix(1.0f0*I, 1, 1),
        (par[:m]*par[:g] .- par[:C]*(extra_state[1]./q[1]).^2)*ones(Float32, 1)
end

function fl_controller(extra_state::AbstractArray, q::AbstractArray, u::AbstractArray)
    m = par[:m]
    g = par[:g]
    C = par[:C]
    L1 = par[:L1]
    R = par[:R]
    qd = par[:qd]

    x = (q[1], u[1], extra_state[1])
    ξ = (x[1]-qd, x[2], g - C/m*(x[3]/x[1])^2)
    L = L1 + 2*C / x[1]

    fx = -4*C*C/m/L*x[2]*(x[3]^2)/(x[1]^4) + 2*R*C/m/L*(x[3]^2)/(x[1]^2) +
        2*C/m*x[2]*(x[3]^2)/(x[1]^3)
    gx = -2*C/m/L*x[3]/(x[1]^2)
    p = 4.6022850                # p = 4.6022850 just grazes the pedestal with
                                 #          voltage_limit = 165.54182708519832
    # p = 2.5
    k = [p^3, 3*p^2, 3*p]
    w = dot(-k, ξ)
    return 1/gx*(w - fx)
end

function curr_cmd_controller(extra_state::AbstractArray, q::AbstractArray, u::AbstractArray)
    m = par[:m]
    g = par[:g]
    C = par[:C]
    L1 = par[:L1]
    R = par[:R]
    qd = par[:qd]

    x = (q[1], u[1], extra_state[1])

    v = -25.0*(x[1]-qd) - 10.0*x[2]
    i_ref = x[1]*sqrt( m/C*(g - v) )
    return R * i_ref - 2 * C * (x[2]*x[3] / x[1] / x[1]) # second term cancels out the nonlinear dynamics
end

function extradynamics(extra_state::AbstractArray, q::AbstractArray, u::AbstractArray)
    m = par[:m]
    g = par[:g]
    C = par[:C]
    L1 = par[:L1]
    R = par[:R]
    qd = par[:qd]

    L = L1 + 2*C / q[1]
    voltage_limit = 165.54182708519832      # This is the minimum voltage needed to pull the ball from δ[2] = 0.2

    control_input = clamp(fl_controller(extra_state, q, u), -voltage_limit, voltage_limit)
    # control_input = clamp(curr_cmd_controller(extra_state, q, u), -voltage_limit, voltage_limit)

    return (-R/L * extra_state[1] + 2 * C / L * (u[1]*extra_state[1] / q[1] / q[1])
        + 1/L*control_input) * ones(Float32, 1)
end

maglev = Integrator(gap, dynamics, q0, u0, extra_state0, Δt=par[:Δt], ε=par[:ε])
integrate(maglev, extradynamics, par[:tf])

fig = figure(1, clear=true, figsize=(12.80,12.80))
fig.add_subplot(2,2,1)
fig.add_subplot(2,2,2)
fig.add_subplot(2,2,3)
fig.add_subplot(2,2,4)
fig.axes[1].clear()
fig.axes[1].plot(maglev.t, getindex.(maglev.q,1), "b", linewidth=2, label="ball position")
fig.axes[2].plot(maglev.t, getindex.(maglev.u,1), "k", linewidth=2, label="ball velocity")
fig.axes[3].plot(maglev.t, getindex.(maglev.extra_state,1), "m", linewidth=2, label="current")
# fig.axes[4].plot(maglev.t[2:end], getindex.(maglev.Λ,1)./maglev.Δt, "m", linewidth=2, label="contact force top")
# fig.axes[4].plot(maglev.t[2:end], getindex.(maglev.Λ,2)./maglev.Δt, "g:", linewidth=2, label="current force bottom")
fig.axes[4].plot(maglev.t, par[:g] .- par[:C]/par[:m]*(getindex.(maglev.extra_state,1)./getindex.(maglev.q,1)).^2,
    "g", linewidth=2, label="contact force top")

fig.axes[1].set_ylabel(L"$z(t)$", fontsize=18)
fig.axes[2].set_ylabel(L"$\dot{z}(t)$", fontsize=18)
fig.axes[3].set_ylabel(L"$i(t)$", fontsize=18)
# fig.axes[4].set_ylabel(L"$\Lambda(t)$", fontsize=18)
fig.axes[4].set_ylabel(L"$\xi_3(t)$", fontsize=18)

# fig.axes[4].legend(loc="best", fontsize="x-large")

for i = 1:4
    fig.axes[i].set_xlabel(L"$t$", fontsize=18)
    # fig.axes[i].autoscale_view()
    fig.axes[i].grid()
end
fig.suptitle("Magnetic levitation", fontsize=18)
fig.tight_layout()


text:
using TensorNetworkCodes
using Test

@testset "plot_code" begin
    code = TensorNetworkCode(five_qubit_code())
    plot_code(code)  # no error thrown
end

@testset "plot_operator" begin
    code = TensorNetworkCode(steane_code())
    plot_operator(code, code.stabilizers[1])  # no error thrown
end


text:
module QuantumRelay

using SymPy
using JuMP
using GLPK
using IterTools
using LLLplus
using GSL
using PyPlot
using Distributions
using LinearAlgebra
using SmithNormalForm

export
	qrelay_op,
	op_mat,
	scan_maker,
        nonneg,
	QRelaySampler

include("utility.jl")
include("operateurs.jl")
include("scan.jl")
include("distributions.jl")

struct QRelaySampler
	prob::Function #return the probablility
	psetproposal::Function #return next combination of sample

    #input parameters:
    #mat: the matrix p_ij in the note (10)
    #coef: the coefficient c in the note (10)
    #omega: the coefficient A in the note (10)
    #pdet0: the probability of detection for each detector
	function QRelaySampler(mat::Array{T, 2}, coef, omega, pdet0) where T <: Int
		F=smith(mat)
                U=F.S
                V=F.T
                S=diagm(F)
        #the PolyLib returns P=USV. Inverse the matrices so Ui/Vi is the same as U/V in the note (18)
		Ui = inv(U)
		Vi = inv(V)
		s = diag(S)
		r = count(!iszero, s)
		s0 = s[1:r]
		@assert s0 == ones(r)
		ui1 = Ui[1:r, :]
		ui2 = Ui[r+1:end, :]
		vi1 = Vi[:, 1:r]
		vi2 = Vi[:, r+1:end]
		vi2 = lll(vi2)[1] #Lenstra–Lenstra–Lovász lattice basis reduction
		T0 = vi1*ui1
		ui2oc = orthocomp(ui2) #orhogonal complement

		setc, scan = scan_maker(vi2) #make the scanner for the algorithm1 in the note

        #compute the probability for an ideal system
        #na: the photon numbers in a output mode
		function prob(na)
		    @assert count(!iszero, ui2*na) == 0
		    b = T0*na
		    setc(-b)
		    total = 0.0
		    for x in Channel(scan)
		        nab = vi2*x + b #the photon numbers for each item in the sum in the note (10)
		        total += prod([c.^complex(n)/factorial(n) for (c, n) in zip(coef, nab)])
		    end
		    return abs2.(total*omega)
		end

        #compute the probability of detection
        #q: the number of photons detectors report
        #na: the number of photons arrived at detector
        #mask: if there is no detector in this channel, mask=0
		function prob(q, na, mask)
		    q0 = round.(Int, q.>0)
		    m0 = round.(Int, mask)
		    return prod((q0 + (ones(length(q0))-2q0).*pdet0(na)).^m0)
		end

		psetproposal(x::Vector) = QuantumRelay.OrthoNNDist(x, ui2oc)

		new(prob, psetproposal)
	end
end

end


text:
# Interface with openBuildNet node

# A callback = a function + optional arguments in a tuple
typealias OBNCallback Tuple{Function,Tuple}
#Base.call(f::OBNCallback) = f[1](f[2]...) # Make the callback callable
(f::OBNCallback)() = f[1](f[2]...) # Make the callback callable
OBNCallback() = OBNCallback((() -> nothing, ()))

typealias InputPortDict Dict{Cint,OBNInputAbstract}

# An openBuildNet node
type OBNNode
  node_id::Csize_t
  valid::Bool

  # Dictionary of input ports
  input_ports::InputPortDict

  # Callbacks
  block_output_cb::Dict{Int64,OBNCallback}
  block_state_cb::Dict{Int64,OBNCallback}

  init_cb::OBNCallback
  restart_cb::OBNCallback
  term_cb::OBNCallback

  function OBNNode(name::String, workspace::String = "", server::String = "")
    myid = Ref{Csize_t}(0)  # To receive the node id
    result = ccall(_api_createOBNNode, Cint, (Cstring, Cstring, Cstring, Ref{Csize_t}), name, workspace, server, myid)
    if result != 0
      error("Error creating openBuildNet node [$result]: ", lastErrorMessage())
    end

    obj = new(myid[],                       # id
              true,                         # valid
              InputPortDict(),              # dictionary of input ports
              Dict{Int64,OBNCallback}(),    # output Callbacks
              Dict{Int64,OBNCallback}(),    # state callbacks
              OBNCallback(),                # init callback
              OBNCallback(),                # restart callback
              OBNCallback())                # terminate callback
    finalizer(obj, x -> delete(x))
    obj
  end
end

# Delete / destroy a node
function delete(node::OBNNode)
  if node.valid
    ccall(_api_deleteOBNNode, Cint, (Csize_t,), node.node_id);
    node.valid = false
    true
  else
    false
  end
end

# Set callback for computing outputs of a block (UPDATE_Y)
function on_block_output(f::Function, node::OBNNode, blkid::Integer, optargs...)
  @assert node.valid "Node is not valid."

  if blkid < 0 || blkid > max_blockid()
    error("Invalid computation block ID.")
  end

  node.block_output_cb[blkid] = OBNCallback( (f, optargs) )
end

function on_block_state(f::Function, node::OBNNode, blkid::Integer, optargs...)
  @assert node.valid "Node is not valid."

  if blkid < 0 || blkid > max_blockid()
    error("Invalid computation block ID.")
  end

  node.block_state_cb[blkid] = OBNCallback( (f, optargs) )
end

function on_init(f::Function, node::OBNNode, optargs...)
  @assert node.valid "Node is not valid."

  node.init_cb = OBNCallback( (f, optargs) )
end

function on_restart(f::Function, node::OBNNode, optargs...)
  @assert node.valid "Node is not valid."

  node.restart_cb = OBNCallback( (f, optargs) )
end

function on_term(f::Function, node::OBNNode, optargs...)
  if !node.valid
    error("Node is not valid.")
  end

  node.term_cb = OBNCallback( (f, optargs) )
end

function do_updatey(node::OBNNode, mask::OBNUpdateMask)
  for (id,f) in node.block_output_cb
    if mask == 0
      break;
    end
    if (mask & (1 << id)) != 0
      # If the id is in the mask, run the callback
      f()
      mask $= (1 << id)   # reset that bit
    end
  end
end

function do_updatex(node::OBNNode, mask::OBNUpdateMask)
  for (id,f) in node.block_state_cb
    if mask == 0
      break;
    end
    if (mask & (1 << id)) != 0
      # If the id is in the mask, run the callback
      f()
      mask $= (1 << id)   # reset that bit
    end
  end
end

import Base.run

"""
    run(OBNNode, timeout, stopIfTimeout)

Run simulation

If stopIfTimeout is true [default], the node will automatically  stop if there is a timeout error

Returns: 1 if timeout; 2 if the simulation stopped properly; 3 if stopped with an error
"""
function run(node::OBNNode, timeout::Number = -1.0, stopIfTimeout::Bool = true)
  @assert node.valid "Node is not valid."

  if iserror(node)
    error("Node currently has error and cannot run; clear the error first by calling stopSimulation.")
  end

  # If the node is running --> may be not right, give a warning
  if isrunning(node)
    warn("Node is currently running; may be an error but we will run it anyway.")
  end

  # timeoutStart = tic;

  event_type = Ref{OBNEI_EventType}()
  event_args = Ref{OBNEI_EventArg}()

  result = 0

  while result == 0
    result = ccall(_api_simRunStep, Cint, (Csize_t, Cdouble, Ref{OBNEI_EventType}, Ref{OBNEI_EventArg}), node.node_id, timeout, event_type, event_args)

    if result == 0
      # Got an event
      if event_type[] == OBNEI_Event_Y
        do_updatey(node, event_args[].mask)
      elseif event_type[] == OBNEI_Event_X
        do_updatex(node, event_args[].mask)
      elseif event_type[] == OBNEI_Event_INIT
        event_result = node.init_cb()
        if isa(event_result, Integer)
          # Set the event result to OBN
          ccall(_api_simSetEventResult, Cint, (Csize_t, Int64), node.node_id, event_result)
        end
      elseif event_type[] == OBNEI_Event_TERM
        node.term_cb()
      elseif event_type[] == OBNEI_Event_RESTART
        event_result = node.restart_cb()
        if isa(event_result, Integer)
          # Set the event result to OBN
          ccall(_api_simSetEventResult, Cint, (Csize_t, Int64), node.node_id, event_result)
        end
      elseif event_type[] == OBNEI_Event_RCV
        # Port's RCV event
        portid = Cint(event_args[].index)
        if haskey(node.input_ports, portid)
          node.input_ports[portid].rcv_cb()
        else
          warn("Input port $portid does not exist.")
        end
      else
        error("Internal error: Unknown event type.")
      end
    elseif result == 1
      # Timeout
      warn("Simulation has timed out.")
      if stopIfTimeout
        # Stop the simulation immediately
        stop(node, true)
      end
      # if toc(timeoutStart) > timeout
      #     % Real timeout error occurred
      #     if stopIfTimeout
      #         this(k).stopSimulation();  % stop the simulation immediately
      #     end
      #     b(k) = false;
      # else
      #     % We can continue running the simulation
      #     status(k) = 0;
      # end
    elseif result == 2
      # Stop properly
      # disp(['Simulation of node ' this(k).obnnode_nodeName ' has stopped properly.']);
    elseif result == 3
      # Stop with error
      warn("Simulation has stopped due to an error: ", lastErrorMessage())
    else
      error("Internal error: Unknown running state ", result)
    end
  end
  result
end

# Similar to run(node) but it catches any exception (error) from Julia and will try to stop the simulation properly before re-throwing the exception.
# This could be useful for scripts (vs. the REPL) because any runtime exception will cause the script to stop abruptly and the simulation may stuck (of course the SMN may terminate the simulation manually).
# Alternatively the user can wrap run(node) inside try...catch or try...finally to achieve the same result.
function runsafe(node::OBNNode, timeout::Number = -1.0, stopIfTimeout::Bool = true)
  @assert node.valid "Node is not valid."

  result = 3
  try
    result = run(node, timeout, stopIfTimeout)
  catch
    stop(node, true)
    rethrow()
  end
  result
end

# Wait for and process the next port event (by executing its callback), up until a given
# timeout (in seconds, can be non-positive for no timeout).
# Returns true if an event has been processed; false otherwise (timeout).
# This method is useful for event-triggered processing inside another callback.
function process_port_events(node::OBNNode, timeout::Number = -1.0)
  @assert node.valid "Node is not valid."

  # Get the next port event
  event_type = Ref{OBNEI_EventType}()
  portid = Ref{Csize_t}()
  result = ccall(_api_simGetPortEvent, Cint, (Csize_t, Cdouble, Ref{OBNEI_EventType}, Ref{Csize_t}), node.node_id, timeout, event_type, portid)

  if result != 0
    return false
  end

  if event_type[] == OBNEI_Event_RCV
    # A message received at a port
    if haskey(node.input_ports, portid[])
      node.input_ports[portid[]].rcv_cb()
    else
      warn("Input port $(portid[]) does not exist.")
      return false
    end
  else
    # An unexpected event is returned -> error
    error("Internal error: unexpected event type $(event_type[]).")
  end
  return true
end

# Stop the simulation.
# If stopnow is false [default], the node will only request the SMN to stop simulation.
# If the SMN accepts the request, it will send a TERM message to all nodes, including this one.
# If stopnow is true, the node will stop its simulation immediately regardless (after sending the request to the SMN).
# Returns true if successful
function stop(node::OBNNode, stopnow::Bool = false)
  @assert node.valid "Node is not valid."

  result = stopnow?ccall(_api_nodeStopSimulation, Cint, (Csize_t,), node.node_id):ccall(_api_nodeRequestStopSimulation, Cint, (Csize_t,), node.node_id)
  result == 0
end

# Converts a given time value in a given unit (default: seconds) to clock ticks.
# Valid units are :s, :m, :h, :ms, :us, :t (clock ticks -> returns the given time value immediately)
# Return type is OBNSimTimeType
function to_ticks(node::OBNNode, t::Number, tu::Symbol = :s)
  if tu == :t
    return t
  end

  atomictu = timeunit(node)

  # Convert t to clock ticks
  const tu_scale = Dict{Symbol,Float64}(:s => 1e6, :ms => 1e3, :us => 1, :m => 60*1e6, :h => 60*60*1e6)
  ceil(OBNSimTimeType, t*tu_scale[tu] / atomictu)
end

# Requests a triggering (i.e., update) of certain blocks at a future time (in clock ticks)
# Here the time is in the clock ticks, i.e., the number of atomic time units from the beginning of the simulation.
# Returns the status of the request:
# 0 if successful (accepted), -1 if timeout (failed), -2 if request is invalid, >0 if other errors (failed, see OBN documents for details).
import Base.schedule
function schedule(node::OBNNode, blks::OBNUpdateMask, t::OBNSimTimeType, timeout::Number = -1.0)
  @assert node.valid "Node is not valid."

  ccall(_api_simRequestFutureUpdate, Cint, (Csize_t, OBNSimTimeType, OBNUpdateMask, Cdouble), node.node_id, t, blks, timeout)
end

# Requests a triggering (i.e., update) of certain blocks at a future time.
# Here, the time is given in simulation time with a given unit (default = seconds) from the beginning of the simulation.
# Valid units are :s, :m, :h, :ms, :us
# This function converts the future time to the atomic clock ticks and call the default method
function schedule(node::OBNNode, blks::OBNUpdateMask, t::Number, tu::Symbol = :s, timeout::Number = -1.0)
  schedule(node, blks, to_ticks(node, t, tu), timeout)
end

# Requests a triggering (i.e., update) of certain blocks at a future time.
# Here, the time is given in wallclock time.
# This function converts the future time to the atomic clock ticks and call the default method
function schedule(node::OBNNode, blks::OBNUpdateMask, t::DateTime, timeout::Number = -1.0)
  # Duration until the requested future time, in milliseconds
  dt = t - wallclock_time(node)
  @assert dt.value>0 "Requested time must be strictly in the future."

  schedule(node, blks, sim_time(node, :ms) + dt.value, :ms, timeout)
end

# Get current simulation time in given unit
# Possible units as symbols: us, ms, s, m, h, t (for atomic ticks)
function sim_time(node::OBNNode, unit::Symbol = :s)
  @assert node.valid "Node is not valid."

  if unit == :t
    return sim_time_ticks(node)
  end

  # Time unit: 0 = second, -1 = millisecond, -2 = microsecond, 1 = minute, 2 = hour
  const timeunits = Dict{Symbol,Cint}(:s => 0, :ms => -1, :us => -2, :m => 1, :h => 2)
  @assert haskey(timeunits, unit) "Invalid time unit."

  curtime = Ref{Cdouble}(0.0)
  result = ccall(_api_nodeSimulationTime, Cint, (Csize_t, Cint, Ref{Cdouble}), node.node_id, timeunits[unit], curtime)
  if result != 0
    error("Error querying the simulation time [$result]: ", lastErrorMessage())
  end
  curtime[]
end

# Get current simulation time in atomic time ticks
function sim_time_ticks(node::OBNNode)
  @assert node.valid "Node is not valid."

  curtime = Ref{OBNSimTimeType}(0)
  result = ccall(_api_nodeSimulationTimeTicks, Cint, (Csize_t, Ref{OBNSimTimeType}), node.node_id, curtime)
  if result != 0
    error("Error querying the simulation time ticks [$result]: ", lastErrorMessage())
  end
  curtime[]
end

# Returns the atomic time unit as an integer in microseconds
function timeunit(node::OBNNode)
  @assert node.valid "Node is not valid."

  tu = Ref{OBNSimTimeType}(0)
  result = ccall(_api_nodeTimeUnit, Cint, (Csize_t, Ref{OBNSimTimeType}), node.node_id, tu)
  if result != 0
    error("Error querying the time unit [$result]: ", lastErrorMessage())
  end
  tu[]
end

# Get current wallclock time as Julia's DateTime
function wallclock_time(node::OBNNode)
  @assert node.valid "Node is not valid."

  unixtime = Ref{Int64}(0)
  result = ccall(_api_nodeWallClockTime, Cint, (Csize_t, Ref{Int64}), node.node_id, unixtime)
  if result != 0
    error("Error querying the wallclock time [$result]: ", lastErrorMessage())
  end

  Dates.unix2datetime(unixtime[])
end

# Internal function to check state of node
function _is_state(node::OBNNode, api::Ptr{Void})
  @assert node.valid "Node is not valid."

  result = ccall(api, Cint, (Csize_t,), node.node_id)
  (result<0)?error("Error querying node's state."):(result>0)
end

# Check node's state
isstopped(node::OBNNode) = _is_state(node, _api_nodeIsStopped)
isrunning(node::OBNNode) = _is_state(node, _api_nodeIsRunning)
iserror(node::OBNNode) = _is_state(node, _api_nodeIsError)


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _kepler_solver{T<:AbstractFloat}(M::T, e::T)
    if 0 <= e <= 1
        # M is assumed to be in the range [-pi, pi], see Markley (1995), page 2.
        # First restrict it to [0, 2pi], then move values above pi to [-pi, 0].
        M = cirrange(M, 2.0*pi)
        M > pi && (M -= 2.0*pi)
        if M == 0
            return 0.0
        else
            # equation (20)
            α = (3.0*pi*pi + 1.6*pi*(pi - abs(M))/(1.0 + e))/(pi*pi - 6.0)
            # equation (5)
            d = 3.0*(1.0 - e) + α*e
            # equation (9)
            q = 2.0*α*d*(1.0 - e) - M*M
            # equation (10)
            r = 3.0*α*d*(d - 1.0 + e)*M + M*M*M
            # equation (14)
            w = (abs(r) + sqrt(q*q*q + r*r))^(2.0/3.0)
            # equation (15)
            E1 = (2.0*r*w/(w*w + w*q + q*q) + M)/d
            # equation (21)
            f0 = E1 - e*sin(E1) - M
            # equation (25)
            f1 = 1.0 - e*cos(E1)
            # equation (26)
            f2 = e*sin(E1)
            # equation (27)
            f3 = 1.0 - f1
            # equation (28)
            f4 = -f2
            # equation (22)
            δ3 = -f0/(f1 - f0*f2/(2.0*f1))
            # equation (23)
            δ4 = -f0/(f1 + δ3*f2/2.0 + δ3*δ3*f3/6.0)
            # equation (24)
            δ5 = -f0/(f1 + δ4*f2/2.0 + δ4*δ4*f3/6.0 + δ4*δ4*δ4*f4/24.0)
            return E1 + δ5 # equation 29
        end
    else
        error("eccentricity must be in the range [0, 1]")
    end
end

"""
    kepler_solver(M, e) -> E

### Purpose ###

Solve Kepler's equation in the elliptic motion regime (\$0 \\leq e \\leq 1\$)
and return eccentric anomaly \$E\$.

### Explanation ###

In order to find the position of a body in elliptic motion (e.g., in the
two-body problem) at a given time \$t\$, one has to solve the [Kepler's
equation](https://en.wikipedia.org/wiki/Kepler%27s_equation)

\$\$ M(t) = E(t) - e\\sin E(t)\$\$

where \$M(t) = (t - t_{0})/P\$ is the mean anomaly, \$E(t)\$ the eccentric
anomaly, \$e\$ the eccentricity of the orbit, \$t_0\$ is the time of periapsis
passage, and \$P\$ is the period of the orbit.  Usually the eccentricity is
given and one wants to find the eccentric anomaly \$E(t)\$ at a specific time
\$t\$, so that also the mean anomaly \$M(t)\$ is known.

### Arguments ###

* `M`: mean anomaly.  This can be either a scalar or an array
* `e`: eccentricity, in the elliptic motion regime (\$0 \\leq e \\leq 1\$)

### Output ###

The eccentric anomaly \$E\$, restricted to the range \$[-\\pi, \\pi]\$.  If an
array of mean anomalies is provided in input, an array of the same length as `M`
is returned.

### Method ###

Many different numerical methods exist to solve Kepler's equation.  This
function implements the algorithm proposed in Markley (1995) Celestial Mechanics
and Dynamical Astronomy, 63, 101
(DOI:[10.1007/BF00691917](http://dx.doi.org/10.1007/BF00691917)).  This method
is not iterative, requires only four transcendental function evaluations, and
has been proved to be fast and efficient over the entire range of elliptic
motion \$0 \\leq e \\leq 1\$.

### Example ###

(1) Find the eccentric anomaly for an orbit with eccentricity \$e = 0.7\$ and
for \$M(t) = 8\\pi/3\$.

``` julia
ecc = 0.7;
E = kepler_solver(8pi/3, ecc)
# => 2.5085279492864223
```

(2) Plot the eccentric anomaly as a function of mean anomaly for eccentricity
\$e = 0\$, \$0.5\$, \$0.9\$.  Recall that `kepler_solver` gives \$E \\in [-\\pi,
\\pi]\$, use `cirrange` to have it in \$[0, 2\\pi]\$.  Use
[PyPlot.jl](https://github.com/stevengj/PyPlot.jl) for plotting.

``` julia
using PyPlot
M=linspace(0, 2pi, 1001)[1:end-1];
for ecc in (0, 0.5, 0.9); plot(M, cirrange(kepler_solver(M, ecc), 2pi)); end
```

### Notes ###

The true anomaly can be calculated with `trueanom` function.
"""
kepler_solver(M::Real, e::Real) =
    _kepler_solver(promote(float(M), float(e))...)

function kepler_solver{R1<:Real,R2<:Real}(M::AbstractArray{R1}, e::R2)
    typee = promote_type(float(R1), float(R2))
    E = similar(M, typee)
    for i in eachindex(M)
        E[i] = kepler_solver(M[i], e)
    end
    return E
end


text:
# An illustration of invariance in Julia Sets.

abstract type Animal end
struct Cat <: Animal name end
struct Dog <: Animal name end

dogs = Set([Dog("rex"), Dog("spot")])
@assert typeof(dogs) == Set{Dog}      # Julia infers this
@assert isa(dogs, Set{Dog})
@assert !isa(dogs, Set{Animal})       # Surprised?

function something(animals::Set{Animal}) end

@assert (try
  something(dogs)                     # Not allowed!
catch e
  1
end) == 1

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
using FreeParameters, Documenter, LaTeXStrings

makedocs(
  sitename = "FreeParameters.jl",
  doctest = false,
  strict = false,
  format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        mathengine = MathJax(Dict(
            :TeX => Dict(
                :equationNumbers => Dict(:autoNumber => "AMS"),
                :Macros => Dict()
            )
        ))
  ),
  clean = false,
  modules = [Documenter, FreeParameters],
  pages = Any[
    "Home" => "index.md",
  ],
)

deploydocs(
           repo = "github.com/CliMA/FreeParameters.jl.git",
           target = "build",
          )


text:
module FluxMPI

include("mpi_extensions.jl")

using .MPIExtensions: Iallreduce!, Ibcast!, JuliaTaskRequest
using CUDA, MPI
using Dates: now
using Flux: params
using Flux.Optimise: AbstractOptimiser
using MPI: Request, Waitall!, Allreduce!, Bcast!
using Zygote: @nograd, Params
using LearnBase: ObsDim
import MLDataUtils: nobs, getobs
import Flux.Optimise: update!, apply!
import Base: getproperty, setproperty!

include("common.jl")
include("synchronize.jl")
include("optimiser.jl")
include("data.jl")

export MPIExtensions, MPI
export local_rank, total_workers, DistributedOptimiser, broadcast_parameters, clean_print, clean_println,
       DistributedDataContainer

end


text:
using Luxor, Test

Drawing(700, 350, "dimensioning-test.svg")
origin()
background("white")
pentagon = ngonside(O, 120, 5, vertices=true)
sethue("black")
setline(.2)
poly(pentagon, :stroke, close=true)

@testset "dimensioning" begin
    d, t = dimension(O, pentagon[4],
    fromextension = [0, 0])

    @test isapprox(d, 102.0780970022448)
    @test t == "102.0780970022448"

    d, t = dimension(pentagon[1], pentagon[2],
        offset        = -60,
        fromextension = [20, 50],
        toextension   = [20, 50],
        textrotation  = 2π/5,
        textgap       = 20,
        format        = (d) -> string(round(d, digits=4), "pts"))

    @test isapprox(d, 119.99999999999997)
    @test t == "120.0pts"

    d, t = dimension(pentagon[2], pentagon[3],
         offset        = -40,
         format        =  string)

    @test isapprox(d, 120.0)
    @test t =="120.0"

    d, t = dimension(pentagon[5], Point(pentagon[5].x, pentagon[4].y),
        offset        = 60,
        format        = (d) -> string("approximately ", round(d, digits=4)),
        fromextension = [5, 5],
        toextension   = [80, 5])

    @test isapprox(d, 97.08203932499367)
    @test t ==  "approximately 97.082"

    d, t = dimension(pentagon[1], midpoint(pentagon[1], pentagon[5]),
        offset               = 70,
        fromextension        = [65, -5],
        toextension          = [65, -5],
        texthorizontaloffset = -5,
        arrowheadlength      = 5,
        format               = (d) ->
            begin
                if isapprox(d, 60.0)
                    string("exactly ", round(d, digits=4), "pts")
                else
                    string("≈ ", round(d, digits=4), "pts")
                end
            end)

    @test isapprox(d, 60.00000000000001)
    @test t == "exactly 60.0pts"

    d, t = dimension(pentagon[1], pentagon[5],
        offset               = 120,
        fromextension        = [5, 5],
        toextension          = [115, 5],
        textverticaloffset   = 0.5,
        texthorizontaloffset = 0,
        textgap              = 5)

    @test isapprox(d, 120.00000000000001)
    @test t == "120.00000000000001"

    finish()
    nothing
end


text:
using CartesianGP
using Base.Test

# Half adder

c = chromosome_half_adder()
g0 = goal_half_adder()
g1 = Goal(2, (0b0110, 0b1001))

@test_approx_eq fitness(c, g0) 1.0
@test_approx_eq fitness(c, g1) (1.0 / 2.0)


text:
using SolitonBEC
using Printf
using Plots
using Elliptic
hbar = 6.62607015e-34 / (2 * pi)
width = 1.42e-6
mass = 6.941 * 1.660539e-27
omega_perp = 2 * pi * 710
l_perp = sqrt(hbar / (mass * omega_perp))
N = 4e3
as = -0.21e-9

interaction_g = abs(2*hbar^2 * as / mass / l_perp^2)
ggg = 2*hbar*omega_perp*abs(as)

print("\ninteraction_g: ", interaction_g)

# mass as total mass or mass per atom??
time_unit = hbar^2/(mass * ggg * N)
space_unit = hbar^3/(mass * ggg^2 * N^2)
velocity_unit = ggg * N / hbar

# --------- Simulation ---------
khaykovich_gpe = Simulation(
  "GPE",
  "barrier",
)

# --------- Apparata ---------
std_apparatus = Apparatus(
  mass, #m (conversion AMU -> kg)
  as, # as
  omega_perp, # ω_perp
  N, #N
  1.77e-11, #γ
  2 * pi * 4,# ω_z
)

# Energy = 1 / 2 * (std_apparatus.as^2 / (width^2) + std_apparatus.as^2 / (2 * l_perp^2) + width^2 / std_apparatus.as^2 - interaction_g * std_apparatus.as^3 / (l_perp * width^2)) * hbar * omega_perp
L = 2 * space_unit

velocity_unit = ggg*N/hbar
print("\n\tvelocity unit: ", velocity_unit, " m/s")
v0 = 1*velocity_unit
x0 =-L/4
# --------- InitialStates ---------
InitialState1 = InitialState(
  "sech", #type
  width*1.8, # width
  v0, #v0
  x0,
)

# --------- Numerics ---------
T = abs(InitialState1.x0)/ v0 * 2
print("\n\ttime_unit: ", time_unit)
print("\n\tspace_unit: ", space_unit)
num = Numerics(
  T, #T
  T * 1e-3, #dt
  L, #S
  L * 1e-3, #ds
)

#--------- Potential ---------

function barrier_height(energy::Float64)
  r = Potential(
    "barrier", # type
    L * 1e-4, #width
    0, #position
    energy, #energy
    0, #ϵ
    10e-6 #a
  )
  return r
end

energy_unit = mass * N^2 * ggg^2 / hbar^2
print("\n\tenergy unit: ", energy_unit, " J")
print("\n==> desired group velocity: ", v0, " m/s\n\n")

GSEnergy = energy_unit * (-N/24)
#print("\nground state energy: ", GSEnergy/hbar, " hbar\n")
phys_vel = 4.918e-3
normd_vel = phys_vel*hbar/ggg/N
Energy = GSEnergy + energy_unit * normd_vel^2*N/2
#print("\ntraveling state energy: ", Energy/hbar, " hbar\n")

## Configurations
configs = []
for barrier_energy in LinRange(0, 1, 3)
  potential = barrier_height(barrier_energy * energy_unit)
  push!(configs, (num, khaykovich_gpe, std_apparatus, potential, InitialState1))
end

mem_limit = 15000000000 #byte
cnt = 1

pyplot()
p = Plots.palette(:rainbow_bgyr_35_85_c72_n256, length(configs) + 3)

## Soliton - barrier collision --------------------------------------------
run_dynamics(configs[2]...)

fig1 = plot(title="|ψ|^2 after collision with barrier",
  xlabel="space [mm]",
  ylabel="|ψ|^2",
  reuse=false,
  size=(800, 400),
  legend=:topleft)

for (num, sim, app, pot, state) in configs

  global cnt
  @printf("\n------Running simulation # %3i \n", cnt)
  estimate = mem_estimate(num)

  if estimate < mem_limit
    #@printf("Estimated memory consumption: %4.1f MiB\n", estimate / (1024^2))

    coeffs = get_coefficients(sim, app, pot, state)

    time, space, ψ, ψ_spect = @time ssfm_solve(num, coeffs)
    #plot_dynamics(time, space, ψ, ψ_spect)
    peak, idx = findmax(abs.(ψ[:, end]).^2)

    # plot!(space*1e3, abs.(ψ[:, end]).^2)

    print("peak: ", peak)
    group_velocity = (space[idx]-state.x0)/time[length(time)]
    print("\n==> measured group velocity: ", group_velocity, " m/s\n")

    print("\n ratio group velocity: ", state.v0/group_velocity, "\n")

    plot!(space * 1e3,
      abs.(ψ[:, end]) .^ 2,
      label="peak potential = $(pot.energy / (sqrt(2 * pi) * pot.width) /1.602176634e-19)   [eV]",
      lw=1,
      color=p[cnt])

    # plot!(space * 1e3,
    #   abs.(map(s -> potential(sim, app, pot, s), space)) / hbar,
    #   lw=1,
    #   ls=:dot,
    #   label=nothing,
    #   color=p[cnt])

    else
    @printf("Estimated memory consumption (%4.1f MiB) exceed maximum!\n", estimate / (1024^2))
  end

  cnt += 1

end
display(fig1)

max_stars_repo_name: lorenzifrancesco/soliton-BEC

max_issues_repo_name: lorenzifrancesco/soliton-BEC

max_forks_repo_name: lorenzifrancesco/soliton-BEC