text:
__precompile__()
module Units

export regconv, toSI, fromSI, get_toSI_converter, get_fromSI_converter

"""
`Conversion` is the container for forward and backward conversion of units
of measure (`toSI` and `fromSI` fields)
"""
immutable Conversion
	toSI :: Function
	fromSI :: Function
end

"""
Dictionary of units conversion
"""
const conversions = Dict{AbstractString, Conversion}()

"""
Registers units conversion to and from SI units
"""
function regconv(units :: AbstractString, toSI, fromSI)
	conversions[units] = Conversion(toSI, fromSI)
	:done
end

function regconv{T <: AbstractString}(units :: Vector{T}, toSI, fromSI)
	for unit in units
		regconv(unit, toSI, fromSI)
	end
end

"""
Converts `value` of `units` of measure to SI units
"""
toSI(value, units :: AbstractString) = conversions[units].toSI(value)

"""
Converts `value` from SI units of measure to `units`
"""
fromSI(value, units :: AbstractString) = conversions[units].fromSI(value)

"""
Produces the function that converts the value from requested `units`
to SI units.

The `units` must be register prior the use of this function
(see `regconv` function)
"""
function get_toSI_converter(units, default = identity)
	c = get(conversions, units, :null)
	if c === :null
		return default
	else
		return c.toSI
	end
end

"""
Produces the function that converts the value from SI `units` to
requested units.

The `units` must be register prior the use of this function
(see `regconv` function)
"""
function get_fromSI_converter(units, default = identity)
	c = get(conversions, units, :null)
	if c === :null
		return default
	else
		return c.fromSI
	end
end

# Registering standard/useful conversions
# Length
regconv("m", identity, identity)
regconv("cm", x -> x / 100, x -> x * 100)
regconv("mm", x -> x / 1000, x -> x * 1000)
regconv("um", x -> x / 1e6, x -> x * 1e6)
regconv("μm", x -> x / 1e6, x -> x * 1e6)
regconv("km", x -> x * 1000, x -> x / 1000)
regconv("nm", x -> x / 1e9, x -> x * 1e9)

# Mass
regconv("kg", identity, identity)
regconv("g", x -> x / 1000, x -> x * 1000)
regconv(["t", "tonne"], x -> x * 1000, x -> x / 1000)

# Time
regconv(["s", "sec", "second", "seconds"], identity, identity)
regconv(["h", "hour", "hours", "hr", "hrs"], x -> x * 3600, x -> x / 3600)
regconv(["min", "minute", "minutes"], x -> x * 60, x -> x / 60)
regconv(["day", "d", "days"], x -> 24 * 3600, x -> x / 24 / 3600)

# Velocity
regconv(["m/s"], identity, identity)
regconv(["km/h", "kmph"], x -> x * 10 / 36, x -> x * 36 / 10)
regconv(["L/m2.h", "L/(m2.h)", "L/(m^2 h)", "L/(m^2.h)"],
	x -> x / 1000 / 3600,
	x -> x * 1000 * 3600)

# Diffusivity
regconv(["m2/s", "m^2/s"], identity, identity)
regconv(["cm2/s", "cm^2/s"], x -> x / 10000, x -> x * 10000)

# Area
regconv(["m2", "m^2"], identity, identity)
regconv(["cm2", "cm^2"], x -> x / 10000, x -> x * 10000)

# Volume
regconv(["m3", "m^3"], identity, identity)
regconv(["L", "l"], x -> x / 1000, x -> x * 1000)
regconv(["ml", "mL", "cm3", "cm^3"], x -> x / 1e6, x -> x * 1e6)

# Temperature
regconv("K", identity, identity)
regconv(["°C", "C"], x -> x + 273.15, x -> x - 273.15)

end


text:
####
#### EngineMode
####

# required
forwardstep!(p, ::EngineMode) = error("must implement")

supportsreverse(::EngineMode) = false
function reversestep!(p, m::EngineMode)
    supportsreverse(m) && error("supportsreverse was true but reversestep! undefined")
    return p
end

# optional
nameof(m::EngineMode) = string(Base.nameof(typeof(m)))
setup!(ui, p, ::EngineMode) = ui
teardown!(ui, p, ::EngineMode) = ui
reset!(p, ::EngineMode) = (reset!(p.model); p)
pausestep!(p, ::EngineMode) = pausestep!(p)
prepare!(ui, p, ::EngineMode) = ui
modeinfo(io1, io2, ui, p, ::EngineMode) = nothing
handlers(ui, p, ::EngineMode) = EventHandler[]

####
#### PassiveDynamics
####

struct PassiveDynamics <: EngineMode end
forwardstep!(p::PhysicsState, ::PassiveDynamics) = forwardstep!(p)

####
#### Controller
####

mutable struct Controller{F} <: EngineMode
    controller::F
    realtimefactor::Float64
end
Controller(controller) = Controller(controller, 1.0)

function setup!(ui::UIState, p::PhysicsState, x::Controller)
    dt = @elapsed x.controller(p.model)
    x.realtimefactor = timestep(p.model) / dt
    return ui
end

function teardown!(ui::UIState, p::PhysicsState, x::Controller)
    zerofullctrl!(getsim(p.model))
    return ui
end

function forwardstep!(p::PhysicsState, x::Controller)
    dt = @elapsed x.controller(p.model)
    x.realtimefactor = timestep(p.model) / dt
    # println("in ctrler: ",p.model.d.ctrl)
    return forwardstep!(p)
end

function modeinfo(io1, io2, ui::UIState, p::PhysicsState, x::Controller)
    println(io1, "Realtime Factor")
    @printf io2 "%.2fx\n" x.realtimefactor
    return nothing
end

####
#### Trajectory
####

mutable struct Trajectory{TR<:AbsVec{<:AbsMat}} <: EngineMode
    trajectories::TR
    k::Int
    t::Int

    burstmode::Bool
    bf_idx::Int
    bf_range::LinRange{Float64}
    bg_idx::Int
    bg_range::LinRange{Float64}
    doppler::Bool
    function Trajectory{TR}(trajectories) where {TR<:AbsVec{<:AbsMat}}
        new{TR}(trajectories, 1, 1, false, 1, LinRange(0, 1, 2), 1, LinRange(0, 1, 2), false)
    end
end
Trajectory(trajectories::AbsVec{<:AbsMat}) = Trajectory{typeof(trajectories)}(trajectories)
Trajectory(trajectories::AbsMat) = Trajectory([trajectories])

function setup!(ui::UIState, p::PhysicsState, m::Trajectory)
    setburstmodeparams!(m, p)
    setstate!(m, p)
    return ui
end

reset!(p::PhysicsState, m::Trajectory) = (m.t = 1; setstate!(m, p); p)

function forwardstep!(p::PhysicsState, m::Trajectory)
    m.t = inc(m.t, 1, getT(m))
    setstate!(m, p)
    return p
end

supportsreverse(::Trajectory) = true
function reversestep!(p::PhysicsState, m::Trajectory)
    m.t = dec(m.t, 1, getT(m))
    setstate!(m, p)
    return p
end

function prepare!(ui::UIState, p::PhysicsState, m::Trajectory)
    if m.burstmode
        bf = round(Int, m.bf_range[m.bf_idx])
        bg = m.bg_range[m.bg_idx]
        burst!(ui, p, gettraj(m), bf, m.t, gamma = bg, doppler = m.doppler)
    end
    return ui
end

function modeinfo(io1, io2, ui::UIState, p::PhysicsState, m::Trajectory)
    println(io1, "t")
    println(io2, "$(m.t)/$(getT(m))")
    println(io1, "k")
    println(io2, "$(m.k)/$(length(m.trajectories))")

    n = length(m.bf_range)
    println(io1, "Burst Mode")
    println(io2, m.burstmode ? "Enabled" : "Disabled")
    println(io1, "Burst Factor")
    println(io2, "$(m.bf_idx)/$n")
    println(io1, "Burst Gamma")
    println(io2, "$(m.bg_idx)/$n")
    println(io1, "Burst Doppler")
    println(io2, m.doppler ? "Enabled" : "Disabled")

    return nothing
end

function handlers(ui::UIState, p::PhysicsState, m::Trajectory)
    return let ui=ui, p=p, m=m
        [
            onscroll(MOD_CONTROL, what = "Change burst factor") do s, ev
                m.bf_idx = clamp(m.bf_idx + ev.dy, 1, length(m.bf_range))
            end,

            onscroll(MOD_SHIFT, what = "Change burst decay rate") do s, ev
                m.bg_idx = clamp(m.bg_idx + ev.dy, 1, length(m.bg_range))
            end,

            onkey(GLFW.KEY_B, MOD_CONTROL, what = "Toggle burst mode") do s, ev
                ispress_or_repeat(ev.action) && (m.burstmode = !m.burstmode)
            end,

            onkey(GLFW.KEY_D, MOD_CONTROL, what = "Toggle burst mode doppler effect") do s, ev
                ispress_or_repeat(ev.action) && (m.doppler = !m.doppler)
            end,

            onkey(GLFW.KEY_UP, what = "Cycle forwards through trajectories") do s, ev
                if ispress_or_repeat(ev.action)
                    m.k = inc(m.k, 1, length(m.trajectories))
                    ax = axes(gettraj(m), 2)
                    checkbounds(Bool, ax, m.t) || (m.t = last(ax))
                    setburstmodeparams!(m, p)
                    setstate!(m, p)
                end
            end,

            onkey(GLFW.KEY_DOWN, what = "Cycle backwards through trajectories") do s, ev
                if ispress_or_repeat(ev.action)
                    m.k = dec(m.k, 1, length(m.trajectories))
                    ax = axes(gettraj(m), 2)
                    checkbounds(Bool, ax, m.t) || (m.t = last(ax))
                    setburstmodeparams!(m, p)
                    setstate!(m, p)
                end
            end,
        ]
    end
end

function setburstmodeparams!(m::Trajectory, p::PhysicsState)
    steps = 40 # n scroll increments

    # Params calibrated on cartpole and scaled accordingly
    # Rendering up to bf0_max states for a len0 length trajectory
    # with a timestep of dt0 and gamma0 looks reasonalbe.
    len0 = 100
    bf0_max = 50
    gamma0 = 0.85
    dt0 = 0.01

    dt = timestep(p.model)
    len = size(gettraj(m), 2)

    bf_max = round(Int, bf0_max * (len / len0) * (dt / dt0))
    gamma = gamma0^(dt / dt0)

    m.bf_range = LinRange{Float64}(1, bf_max, steps) # render between 1 and bf_max states
    m.bg_range = LinRange{Float64}(gamma, 1, steps)

    m.bf_idx = round(Int, steps / 2)
    m.bg_idx = steps

    return m
end

@inline function setstate!(m::Trajectory, p::PhysicsState)
    LyceumMuJoCo.setstate!(p.model, view(gettraj(m), :, m.t))
    return m
end

getT(m::Trajectory) = size(m.trajectories[m.k], 2)
gettraj(m::Trajectory) = m.trajectories[m.k]

####
#### Util
####

function burst!(
    ui::UIState,
    p::PhysicsState,
    states::AbstractMatrix,
    n::Integer,
    t::Integer;
    gamma::Real = 0.9995,
    alphamin::Real = 0.05,
    alphamax::Real = 0.55,
    doppler::Bool = true,
)
    T = size(states, 2)

    T >= n > 0 || error("n must be in range [1, size(states, 2)]")
    0 < t || error("t must be > 0")
    0 < gamma <= 1|| error("gamma must be in range (0, 1)")
    0 < alphamin <= 1 || error("alphamin must be in range (0, 1]")
    0 < alphamax <= 1 || error("alphamin must be in range (0, 1]")

    scn = ui.scn
    sim = getsim(p.model)
    n = min(n, fld(MAXGEOM, sim.m.ngeom))
    geoms = unsafe_wrap(Array, scn[].geoms, scn[].maxgeom)

    function color!(tprime, from)
        for i = from:scn[].ngeom
            geom = @inbounds geoms[i]
            if geom.category == Int(MJCore.mjCAT_DYNAMIC)
                dist = abs(tprime - t)
                r, g, b, alpha0 = geom.rgba

                alpha = alpha0 * gamma^dist

                if doppler
                    g = 0
                    if tprime < t
                        beta = (dist + 1) / t / 2 + 0.5
                        r = beta * r
                        b = (1 - beta) * b
                    else
                        beta = dist / (T - t) / 2 + 0.5
                        r = (1 - beta) * r
                        b = beta * b
                    end
                    r = clamp(r, 0, 1)
                    b = clamp(b, 0, 1)
                end

                geoms[i] = @set!! geom.rgba = SVector{4,Cfloat}(r, g, b, alpha)
            end
        end
    end

    LyceumMuJoCo.setstate!(p.model, view(states, :, t))
    mjv_updateScene(
        sim.m,
        sim.d,
        ui.vopt,
        p.pert,
        ui.cam,
        MJCore.mjCAT_ALL,
        scn,
    )

    if n > 1
        fromidx = scn[].ngeom + 1
        from = scn[].ngeom + 1
        for tprime in Iterators.map(x -> round(Int, x), LinRange(1, T, n))
            if tprime != t
                LyceumMuJoCo.setstate!(p.model, view(states, :, tprime))
                mjv_addGeoms(
                    sim.m,
                    sim.d,
                    ui.vopt,
                    p.pert,
                    MJCore.mjCAT_DYNAMIC,
                    scn,
                )
                color!(tprime, from)
            end
            from = scn[].ngeom + 1
        end
    end

    # reset the model to the state it had when it was passed in
    LyceumMuJoCo.setstate!(p.model, view(states, :, t))

    return ui
end

function pausestep!(p::PhysicsState)
    sim = getsim(p.model)
    mjv_applyPerturbPose(sim.m, sim.d, p.pert, 1)
    forward!(sim)
    return p
end

function forwardstep!(p::PhysicsState)
    sim = getsim(p.model)
    fill!(sim.d.xfrc_applied, 0)
    mjv_applyPerturbPose(sim.m, sim.d, p.pert, 0)
    mjv_applyPerturbForce(sim.m, sim.d, p.pert)
    step!(p.model)
    return p
end

ext: jl


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

export BrickletAnalogOut
"""
Generates configurable DC voltage between 0V and 5V
"""
mutable struct BrickletAnalogOut <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

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

        return new(ipcon, deviceInternal)
    end
end

export set_voltage
"""
    $(SIGNATURES)

Sets the voltage. Calling this function will set
the mode to 0 (see :func:`Set Mode`).
"""
function set_voltage(device::BrickletAnalogOut, voltage)
    device.deviceInternal.set_voltage(voltage)
end

export get_voltage
"""
    $(SIGNATURES)

Returns the voltage as set by :func:`Set Voltage`.
"""
function get_voltage(device::BrickletAnalogOut)
    return device.deviceInternal.get_voltage()
end

export set_mode
"""
    $(SIGNATURES)

Sets the mode of the analog value. Possible modes:

* 0: Normal Mode (Analog value as set by :func:`Set Voltage` is applied)
* 1: 1k Ohm resistor to ground
* 2: 100k Ohm resistor to ground
* 3: 500k Ohm resistor to ground

Setting the mode to 0 will result in an output voltage of 0 V. You can jump
to a higher output voltage directly by calling :func:`Set Voltage`.
"""
function set_mode(device::BrickletAnalogOut, mode)
    device.deviceInternal.set_mode(mode)
end

export get_mode
"""
    $(SIGNATURES)

Returns the mode as set by :func:`Set Mode`.
"""
function get_mode(device::BrickletAnalogOut)
    return device.deviceInternal.get_mode()
end

export get_identity
"""
    $(SIGNATURES)

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

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

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


text:
module Kuber

using JSON
using Swagger
using HTTP

include("api/Kubernetes.jl")
using .Kubernetes
import .Kubernetes: getAPIVersions

import Base: convert, get, put!, delete!, show
import Swagger: SwaggerModel

include("typealiases.jl")
include("apialiases.jl")
include("helpers.jl")
include("simpleapi.jl")

export KuberContext, set_server, set_ns, get_server, get_ns, kuber_type, kuber_obj, @K_str
export get, put!, update!, delete!, sel, get_logs

end # module


text:
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

compute(n::Int)::Int = n * (n + 1) * (3n ^ 2 - n - 2) ÷ 12

compute(10)

compute(100)

@benchmark compute(100)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using Flux
using Flux:@adjoint

function kan_rv_prod_inner(phi::Array{AbstractFloat,2}, v::NTuple)
	vl = collect(v)
	h = 0.5 .- vl
	hKh =  ( 0.5 * dot(h, phi * h))^(size(phi)[1] / 2)
	(-1.0)^(sum(vl)) * hKh
end
"""
calculates the expectation of products of Gaussian RVs, with
exponents s and covariance matrix phi, using the identity in proposition 1
of Kan 2008.
"""
function kan_rv_prod(phi)
	st = size(phi)[1]
	mapreduce(v -> kan_rv_prod_inner(phi, v), +, Iterators.product(fill(0:1, st)...)) / factorial(st ÷ 2)
end

function kan_rv_prod2(phi)
    kan_rv_prod(phi)
end

function kan_rv_prod2_adj(phi)
    st = size(phi)[1]
    g = v -> gradient(kan_rv_prod_inner, phi, v)[1]
    sum(g, Iterators.product(fill(0:1, st)...)) / factorial(st ÷ 2)
    # gradient(kan_rv_prod, phi)[1]
end

@adjoint function kan_rv_prod2(phi)
    kan_rv_prod2(phi), x -> (x * kan_rv_prod2_adj(phi),)
end

println(kan_rv_prod2_adj(3 * ones(5, 5)))
# gradient(kan_rv_prod_inner, ones(2, 2), (1, 0))
# print(gradient(kan_rv_prod, ones((2, 2)))[1] - gradient(kan_rv_prod2, ones((2, 2)))[1])
"""
@btime VolterraGP.kan_rv_prod(ones((5, 5)))

@btime gradient(kan_rv_prod, ones((2, 2)))

@btime gradient(kan_rv_prod2, ones((2, 2)))

"""
# function sum2(op, arr)
# 	return sum(op, arr)
# end

# function sum2adj(Δ, op, arr)
# 	n = length(arr)
# 	g = x -> Δ * gradient(op, x)[1]
# 	return (nothing, map(g, arr))
# end

# @adjoint function sum2(op, arr)
# 	return sum2(op, arr), Δ -> sum2adj(Δ, op, arr)
# end

# function t1(arr, op)
# 	g = sum2(op, arr)
# 	return g
# end

# function t0(arr, op)
# 	g = sum(op, arr)
# 	return g
# end

# function testPerf()
	# nb = 100000
# 	arr = rand(nb)
# 	op = x -> x + x

# 	res = t1(arr, op)
# 	println(res)

# 	df = x -> gradient(x -> t1(x, op), x)
# 	res1 = df(arr)[1]

# 	df0 = x -> gradient(x -> t0(x, op), x)
# 	res0 = df0(arr)[1]

# 	diff = sum((res0 - res1).^2)
# 	println("diff : ")
# 	println(diff)

# 	println("Benchmarking t1")
# 	@time for i = 1:100
# 		t1(arr, op)
# 	end

# 	println("Benchmarking gradient")
# 	@time for i = 1:100
# 		df(arr)
# 	end

# end

# testPerf()


text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4"]);
gb = SignatureGB.f5([x1, x2, x3, x4]);F = [3136*x1^4+6944*x1^2*x2^2-5823*x1^2*x3^2+8176*x1^2*x4^2-2414*x1*x2^2*x3-10650*x1*x2*x3^2+4133*x2^4+2550*x2^3*x3-6403*x2^2*x3^2+9052*x2^2*x4^2+9409*x3^4-14162*x3^2*x4^2+5329*x4^4-9744*x1^3+11360*x1^2*x3-13508*x1*x2^2-12000*x1*x2*x3+16878*x1*x3^2-6248*x1*x3*x4-12702*x1*x4^2+1496*x2^2*x4+6600*x2*x3*x4+13969*x1^2-11644*x1*x3-7040*x1*x4+2788*x2^2+12300*x2*x3+1936*x4^2-13120*x1+7216*x4+6724, 13888*x1^2*x2-4828*x1*x2*x3-10650*x1*x3^2+16532*x2^3+7650*x2^2*x3-12806*x2*x3^2+18104*x2*x4^2-27016*x1*x2-12000*x1*x3+2992*x2*x4+6600*x3*x4+5576*x2+12300*x3, -11646*x1^2*x3-2414*x1*x2^2-21300*x1*x2*x3+2550*x2^3-12806*x2^2*x3+37636*x3^3-28324*x3*x4^2+11360*x1^2-12000*x1*x2+33756*x1*x3-6248*x1*x4+6600*x2*x4-11644*x1+12300*x2, 16352*x1^2*x4+18104*x2^2*x4-28324*x3^2*x4+21316*x4^3-6248*x1*x3-25404*x1*x4+1496*x2^2+6600*x2*x3-7040*x1+3872*x4+7216];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
function stemleaf{T<:Real}(a::Array{T,1}, leafsize=1)
    ls = 10^int(log10(leafsize))
    (stem, leaf) = divrem(sort(int(a/ls)), 10)
    leaf[sign(stem) .== -1] *= -1
    negzero = leaf .< 0
    if any(negzero)
        leaf[negzero] *= -1
        nz = @sprintf "%10s | " "-0"
        nz *= join(map(string, leaf[negzero]), " ")
        nz *= "\n"
        stem = stem[!negzero]
        leaf = leaf[!negzero]
    else
        nz = ""
    end
    slp = ""
    for i in stem[1]:stem[end]
        i != 0 || (slp *= nz)
        slp *= @sprintf "%10d | " i
        slp *= join(map(string, leaf[stem .== i]), " ")
        slp *= "\n"
    end
    slp *= " Leaf Unit = " * string(convert(T, ls)) * "\n"
    return slp
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
immutable PhotoReceptorModel{P<:GeneralizedMetropolisHastings.AbstractParameter} <: AbstractModel

    #Generic model specs
    name::AbstractString
    parameters::Vector

    photons::GeneralizedMetropolisHastings.DataArray
    measurements::GeneralizedMetropolisHastings.DataArray
    noisemodel::GeneralizedMetropolisHastings.AbstractNoiseModel

    policy::PhotoReceptorPolicy
    receptor::AbstractPhotoReceptor

    ###variable holding the result of a model evaluation
    lightinducedcurrent::Vector{Float32}

    ###temp location to store model data
    fixedbumpvals::Vector{Float32}
    latency::Vector{Float32}
    refractory::Vector{Float32}

    ###Inner constructor
    function PhotoReceptorModel(parameters::Vector{P},photons::GeneralizedMetropolisHastings.DataArray,measurements::GeneralizedMetropolisHastings.DataArray,
                                noisemodel::GeneralizedMetropolisHastings.AbstractNoiseModel,policy::PhotoReceptorPolicy,receptor::AbstractPhotoReceptor,fixedbumpvals)
        @assert (traitvalue(policy.bump) == :fixed && length(fixedbumpvals) > 0) || (traitvalue(policy.bump) == :sample && length(fixedbumpvals) == 0)
        n = "PhotoReceptorModel"
        nsteps = numtimesteps(receptor)
        nvilli = numvilli(receptor)
        nphotons = numphotons(receptor)
        new(n,parameters,photons,measurements,noisemodel,policy,receptor,zeros(Float32,nsteps),fixedbumpvals,zeros(Float32,nphotons),zeros(Float32,nphotons))
    end
end

###Factory functions
_paramprior(::Type{Val{:uniform}},args...) = Distributions.Uniform(args...)
_paramprior(::Type{Val{:normal}},args...) = Distributions.Normal(args...)
_paramprior(::Type{Val{:lognormal}},args...) = Distributions.LogNormal(args...)

function _parameters(::Type{Val{:photoreceptor}},policy::PhotoReceptorPolicy,args...)
    latencykeys = paramkeys(policy.latency)
    refractorykeys = paramkeys(policy.refractory)
    bumpkeys = paramkeys(policy.bump)
    latencypriors = parampriors(policy.latency)
    refractorypriors = parampriors(policy.refractory)
    bumppriors = parampriors(policy.bump)
    numparams = 4 + length(bumpkeys)
    params = Vector{GeneralizedMetropolisHastings.AbstractParameter}(numparams)
    pcounter = 0
    acounter = 0
    for i=1:length(latencykeys)
        params[pcounter+=1] = GeneralizedMetropolisHastings.parameter(latencykeys[i],_paramprior(Val{latencypriors[i]},args[acounter+=1]...))
    end
    if length(latencykeys) == 1
        params[pcounter+=1] = GeneralizedMetropolisHastings.parameter(:unused,0.0f0)
    end
    for i=1:length(refractorykeys)
        params[pcounter+=1] = GeneralizedMetropolisHastings.parameter(refractorykeys[i],_paramprior(Val{refractorypriors[i]},args[acounter+=1]...))
    end
    if length(refractorykeys) == 1
        params[pcounter+=1] = GeneralizedMetropolisHastings.parameter(:unused,0.0f0)
    end
    for i=1:length(bumpkeys)
        params[pcounter+=1] = GeneralizedMetropolisHastings.parameter(bumpkeys[i],_paramprior(Val{bumppriors[i]},args[acounter+=1]...))
    end
    params
end

function _model(::Type{Val{:photoreceptor}},parameters::Vector,photons::GeneralizedMetropolisHastings.DataArray,measurements::GeneralizedMetropolisHastings.DataArray,
                variance::Vector,nvilli::Int,policy::PhotoReceptorPolicy)
    @assert numvalues(photons) == numvalues(measurements)
    @assert (traitvalue(policy.bump) == :fixed && length(parameters) == 4) || (traitvalue(policy.bump) == :sample && length(parameters) == 7)
    n = noise(:gaussian,variance)
    r = photoreceptor(:cell,datavalues(photons),nvilli,nvilli<typemax(UInt16)?UInt16:UInt32)
    b = traitvalue(policy.bump)==:fixed?bump(fixedbumpvalues(policy.fixedbump)...):Float32[]
    PhotoReceptorModel{eltype(parameters)}(parameters,photons,measurements,n,policy,r,b)
end

###evaluate! function

function evaluate!(model::PhotoReceptorModel,vals::AbstractVector,times::Int)
    nvals = numvalues(model.photons)
    result = zeros(Float32,nvals,times)
    for i=1:times
        l = evaluate!(model,vals)
        copy!(result,(i-1)*nvals+1,l,1,nvals)
    end
    result
end

function evaluate!(model::PhotoReceptorModel,vals::AbstractVector)
    _resettemp(model)
    bumpvals = _bumpvals(traittype(model.policy.bump),model,vals)
    _latency!(traittype(model.policy.latency),model,vals)
    _refractory!(traittype(model.policy.refractory),model,vals)
    lightinducedcurrent!(model.lightinducedcurrent,model.receptor,model.latency,model.refractory,bumpvals)
end

###Utility functions used in generic implementations in AbstractModel
@inline dataindex(model::PhotoReceptorModel) = dataindex(model.measurements)
@inline measurements(model::PhotoReceptorModel) = datavalues(model.measurements)
@inline noisemodel(model::PhotoReceptorModel) = model.noisemodel
lightinducedcurrent(model::PhotoReceptorModel) = model.lightinducedcurrent

@inline _setval!{N<:Number}(vals::Vector{N},v::N) = @simd for i=1:length(vals) @inbounds vals[i] = v end

@inline _resettemp(model::PhotoReceptorModel) = _setval!(model.lightinducedcurrent,0.0f0)

@inline _latency!(::Type{Val{:deterministic}},m::PhotoReceptorModel,vals::AbstractVector) = _setval!(m.latency,Float32(vals[1]))
@inline function _latency!(::Type{Val{:lognormal}},m::PhotoReceptorModel,vals::AbstractVector)
    randlognormal!(m.latency,Float32(vals[1]),Float32(vals[2]))
end

@inline _refractory!(::Type{Val{:deterministic}},m::PhotoReceptorModel,vals::AbstractVector) = _setval!(m.refractory,Float32(vals[3]))
@inline function _refractory!(::Type{Val{:lognormal}},m::PhotoReceptorModel,vals::AbstractVector)
    randlognormal!(m.refractory,Float32(vals[3]),Float32(vals[4]))
end

@inline _bumpvals(::Type{Val{:fixed}},m::PhotoReceptorModel,vals::AbstractVector) = m.fixedbumpvals
@inline function _bumpvals(::Type{Val{:sample}},m::PhotoReceptorModel,vals::AbstractVector)
    bump(Float32(vals[5]),Float32(vals[6]),Float32(vals[7]),m.policy.fixedbump.cutoff)
end

###Base functionality
function show(io::IO,m::PhotoReceptorModel)
    println(io,"Model ",m.name)
    print(io,"parameters: ") ; show(io,m.parameters)
    print(io,"photons: ") ; show(io,m.photons)
    print(io,"measurements: ") ; show(io,m.measurements)
    print(io,"noisemodel: ") ; show(io,m.noisemodel)
    print(io,"policy: ") ; show(io,m.policy)
    print(io,"receptor: ") ; show(io,m.receptor)
    println(io)
    nothing
end


text:
if VERSION < v"0.5.0-dev"
    using BaseTestNext
else
    using Base.Test
end
using ExpressionUtils

include("functions.jl")
include("utils.jl")


text:
module PerfTests

using BenchmarkTools
using PBC

using PBC.Curve: EP, EP2, FP, BN
using PBC.Shamir: LagrangeCoeffGenerator, BarycentricWeightGenerator, lagrange_interpolate_c0

BenchmarkTools.DEFAULT_PARAMETERS.samples = 1000
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10
BenchmarkTools.DEFAULT_PARAMETERS.evals = 1
BenchmarkTools.DEFAULT_PARAMETERS.gctrial = true
BenchmarkTools.DEFAULT_PARAMETERS.gcsample = false
@show BenchmarkTools.DEFAULT_PARAMETERS

Base.rand(::Type{PublicKey}) = PublicKey(rand(PBC.Config.@EP2))
Base.rand(::Type{Signature}) = Signature(rand(PBC.Config.@EP))
Base.rand(::Type{Hash}) = Hash(rand(PBC.Config.@EP))
Base.rand(::Type{Identity}) = Identity(rand(PublicKey))
Base.rand(::Type{PublicKeyPoly}; n=500) = PublicKeyPoly(PrivateKeyPoly(n))
function Base.rand(::Type{SignatureShares}; n=501)
    sigshares = SignatureShares()
    for i in 1:n
        sigshares[rand(Identity)] = Signature(rand(PBC.Config.@EP))
    end
    return sigshares
end

const x = 1:1000

suite = BenchmarkGroup()
suite["PBC"] = BenchmarkGroup()
# suite["PBC"]["rand(PrivateKey)"] = @benchmarkable rand(PrivateKey)
# suite["PBC"]["PrivateKey(::Vector{UInt8})"] = @benchmarkable PrivateKey($(rand(UInt8, PBC.Config.PRIVATE_KEY_SIZE)))
# suite["PBC"]["PublicKey(::PrivateKey)"] = @benchmarkable PublicKey($(rand(PrivateKey)))
# suite["PBC"]["PublicKey(::Vector{UInt8})"] = @benchmarkable PublicKey($(Vector{UInt8}(rand(PublicKey))))
# suite["PBC"]["Vector{UInt8}(::PrivateKey)"] = @benchmarkable Vector{UInt8}($(rand(PrivateKey)))
# suite["PBC"]["Vector{UInt8}(::PublicKey)"] = @benchmarkable Vector{UInt8}($(rand(PublicKey)))
# suite["PBC"]["PBC.sign(::PrivateKey, ::Hash)"] = @benchmarkable PBC.sign($(rand(PrivateKey)), $(rand(Hash)))
suite["PBC"]["PBC.verify(::Signature, ::PublicKey, ::Hash)"] = @benchmarkable PBC.verify($(rand(Signature)), $(rand(PublicKey)), $(rand(Hash)))
# suite["PBC"]["PBC.verify_share(::PublicKeyPoly, ::Identity, ::PrivateKey)"] = @benchmarkable PBC.verify_share($(rand(PublicKeyPoly, n=500)), $(rand(Identity)), $(rand(PrivateKey)))
# suite["PBC"]["::PublicKeyPoly + ::PublicKeyPoly"] = @benchmarkable $(rand(PublicKeyPoly, n=500)) + $(rand(PublicKeyPoly, n=500))
# suite["PBC"]["Signature(::SignatureShares)"] = @benchmarkable Signature($(rand(SignatureShares, n=501)))
# suite["PBC"]["isvalid(::PrivateKey)"] = @benchmarkable isvalid($(rand(PrivateKey)))

# suite["Shamir"] = BenchmarkGroup()
# suite["Shamir"]["lagrange_interpolate_c0(::LagrangeCoeffGenerator, ::Vector{BN})"] = @benchmarkable lagrange_interpolate_c0($(LagrangeCoeffGenerator(x)), $(Dict(i=>rand(BN) for i in x)))
# suite["Shamir"]["lagrange_interpolate_c0(::LagrangeCoeffGenerator, ::Vector{EP})"] = @benchmarkable lagrange_interpolate_c0($(LagrangeCoeffGenerator(x)), $(Dict(i=>rand(EP) for i in x)))
# suite["Shamir"]["lagrange_interpolate_c0(::LagrangeCoeffGenerator, ::Vector{EP2})"] = @benchmarkable lagrange_interpolate_c0($(LagrangeCoeffGenerator(x)), $(Dict(i=>rand(EP2) for i in x)))

# suite["Shamir"]["lagrange_interpolate_c0(::BarycentricWeightGenerator, ::Vector{BN}))"] =  @benchmarkable lagrange_interpolate_c0($(BarycentricWeightGenerator(x)), $(Dict(i=>rand(BN) for i in x)))
# suite["Shamir"]["lagrange_interpolate_c0(::BarycentricWeightGenerator, ::Vector{EP}))"] =  @benchmarkable lagrange_interpolate_c0($(BarycentricWeightGenerator(x)), $(Dict(i=>rand(EP) for i in x)))
# suite["Shamir"]["lagrange_interpolate_c0(::BarycentricWeightGenerator, ::Vector{EP2}))"] =  @benchmarkable lagrange_interpolate_c0($(BarycentricWeightGenerator(x)), $(Dict(i=>rand(EP2) for i in x)))

# suite["Shamir"]["LagrangeCoeffGenerator(1:501)"] = @benchmarkable LagrangeCoeffGenerator($(1:501))
# suite["Shamir"]["BarycentricWeightGenerator(1:501)"] = @benchmarkable BarycentricWeightGenerator($(1:501))
# suite["Shamir"]["LagrangeCoeffGenerator(1:1000) - delete 499 keys"] = @benchmarkable for i in ids; delete!(lc, i); end setup=(ids=1:501; lc=LagrangeCoeffGenerator(ids))
# suite["Shamir"]["BarycentricWeightGenerator(1:1000) - delete 499 keys"] = @benchmarkable for i in ids; delete!(lc, i); end setup=(ids=1:501; lc=BarycentricWeightGenerator(x))

function format_trial(suite, group, res)
    a = allocs(res)
    gct = BenchmarkTools.prettytime(gctime(res))
    t = BenchmarkTools.prettytime(time(res))
    m = BenchmarkTools.prettymemory(memory(res))
    return "[$suite] $group: $t (alloc: $a, mem: $m, gc: $gct)"
end

# If a cache of tuned parameters already exists, use it, otherwise, tune and cache
# the benchmark parameters. Reusing cached parameters is faster and more reliable
# than re-tuning `suite` every time the file is included.
paramspath = joinpath(@__DIR__, "params.json")
if isfile(paramspath)
    loadparams!(suite, BenchmarkTools.load(paramspath)[1], :evals);
else
    println("First run - tuning params (please be patient) ...")
    tune!(suite)
    BenchmarkTools.save(paramspath, params(suite));
end

# print the results
results = run(suite, verbose = true)
for suiteres in results
    for groupres in suiteres.second
        msg = format_trial(suiteres.first, groupres.first, groupres.second)
        println(msg)
    end
end

end


text:
function get_num_markers_lab(robot)
    #УТВ: Робот - в произвольной клетке лабиринта

    num_steps, side_steps = move_to_start!(robot, (Sud,West))
    num_steps=reverse(num_steps)
    side_steps=reverse(side_steps)
    #УТВ: Робот - в юго-западном углу

    side = Ost
    num_markers = get_num_markers!(robot, side)
    #ИНВАРИАНТ: Робот - в "начале" ряда
    while is_move_posible!(robot, Nord, side) == true
        move!(robot, Nord)
        movements!(robot,side)
        side = inverse(side)
        num_markers += get_num_markers!(robot, side)
    end
    #УТВ: Маркеры посчитаны, Робот где-то у северной границы в углу

    move_to_start!(robot, (Sud,West))
    #УТВ: Робот - снова в юго-западном углу

    for (i,n) in enumerate(num_steps)
        movements!(robot, side_steps[i], n)
    end
    #УТВ: Робот - в исходном положении
    return num_markers
end

function get_num_markers!(robot, side)
#=
    ДАНО: Робот - в начале ряда, side - направление перемещения
    РЕЗУЛЬТАТ: Робот - в противоположном конце ряда и возвращено число маркеров, имеющихся в этом ряду
=#
    nun_markers = ismarker(robot) ? 1 : 0
    while !isborder(robot,side)
        move!(robot,side)
        ismarker(robot) && num_markers += 1 # - это тоже самое, что и if ismarker(robot) num_markers += 1 end
    end
    return num_markers
end

function is_move_posible!(robot, move_side, find_side)
#=
    ДАНО: Робот - в одном из концов ряда
        -- find_side - направление к противоположному концу ряда
        -- move_side - направление для перемещения в следующий ряда
    РЕЗУЛЬТАТ:
        Проверена возможность переместить Робота (но не перемещет его) в направлении move_side,
        при необходимости выполнив смещение Робота в направлении find_side на некоторое количество шагов

        Если перемещение в направлении move_side возможно, то Робот остановится возле первого возможного прохода,
        и тогда будет возвращено - true;
        в противном случае будет возвращено - false, и Робот окажется у противоположного конца ряда.
=#
    while isborder(robot, move_side)
        if !isborder(robot, find_side)
            move!(robot, find_side)
        else
            return false
        end
    end
    return true
end

function move_to_start(robot, start_side::Tuple)
#=
    Перемещает Робота из произвольной клетки в угол, заданный кортежем из двух ортогональных направлений start_side,
    и возвращает два массива равной длины:
    первый массив содержит числа сделанных шагов в направлениях, содержащихся во втором массиве, в естественном порядке
=#
    num_steps = Int[]
    side_steps = HorizoSide[]
    while !isborder(robot,start_side[2]) || is_posible!(robot, start_side[1], inverse(start_side[2]), num_steps, side_steps) # <=> Робот - не в стартовом углу
        for i in 1:2
            push!(num_steps, get_num_movements!(robot,start_side[i]))
            push!(side_steps, start_side[i])
        end
    end
    return num_steps, side_steps
end

function is_move_posible!(robot, move_side, find_side, num_steps::Vector, side_steps::Vector)
#=
    ДАНО: Робот - в одном из концов ряда
        -- find_side - направление к противоположному концу ряда
        -- move_side - направление для перемещения в следующий ряда
        -- num_steps, side_steps - внешние массивы, которые додны быть пополнены новыми элементами
    РЕЗУЛЬТАТ:
        Проверена возможность переместиться в направлении move_side, при необходимости сместившись
        в направлении find_side на некотое количество шагов
        Если перемещение возможно, то Робот остановится возле первого возможного прохода, и будет возвращено -  true
        В противном случае будет возвращено - false, и Робот окажется у противоположного конца ряда

        При этом во внешние массивы num_steps, side_steps добавлены число сделанных при поиске прохода шагов и
        направление поиска, соответственно
=#
    num=0
    while isborder(robot, move_side)
        if !isborder(robot, find_side)
            num += 1
            move!(robot,find_side)
        else
            push!(num_steps, num)
            push!(side_steps, find_side)
            return false
        end
    end
    push!(num_steps, num)
    push!(side_steps, find_side)

    return true
end


text:
#Tests dispersion calculations

f = linspace(100THz, 200THz, 20)
neff = collect(linspace(1, 2, 20).^2)
aeff = OpticalAttr(1.0e-12m^-2, "Aeff")
attr = OpticalAttr(f, neff, "neff")
mode = Mode(attr, aeff)
@test attr(f[10]) ≈ neff[10] # fit function should be close to data
@test NonlinearOptics.der(attr, f[end]; order=0) ≈ 4.0
# =neff + 2*neff in this special case
@test get_groupindex(mode, f[1]) ≈ 3.0
attr_nondispersive = OpticalAttr(1.0, "nondispersive")
mode_nondispersive = Mode(attr_nondispersive, aeff)
@test_throws Exception get_beta(mode_nondispersive, f[1], 1)


text:
using DirectSearch
using Test

include("./test_OrthoMADS.jl")

include("./test_LTMADS.jl")
include("./test_Cache.jl")
include("./test_Search.jl")
include("./test_Constraints.jl")
include("./test_Core.jl")
include("./solve_tests.jl")


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Clingcon_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Clingcon")
JLLWrappers.@generate_main_file("Clingcon", UUID("f3fadb3f-5422-5a6c-be27-a20f6ffb29ff"))
end  # module Clingcon_jll


text:
export BatchExplorer

using Flux

"""
    BatchExplorer(explorer::AbstractExplorer)
"""
struct BatchExplorer{E} <: AbstractExplorer
    explorer::E
end

"""
    (x::BatchExplorer)(values::AbstractMatrix)

Apply inner explorer to each column of `values`.
"""
(x::BatchExplorer)(values::AbstractMatrix) = [x.explorer(v) for v in eachcol(values)]

(x::BatchExplorer)(values::AbstractMatrix, mask::AbstractMatrix) =
    [x.explorer(v, m) for (v, m) in zip(eachcol(values), eachcol(mask))]

(x::BatchExplorer)(v::AbstractVector) = x.explorer(v)
(x::BatchExplorer)(v::AbstractVector, m::AbstractVector) = x.explorer(v, m)

Flux.testmode!(x::BatchExplorer, mode = true) = testmode!(x.explorer, mode)


text:
using Oceananigans: AbstractModel
using Oceananigans.Grids
using Oceananigans.Utils: tupleit

struct KernelComputedField{X, Y, Z, S, A, G, K, F, P} <: AbstractField{X, Y, Z, A, G}
                  data :: A
                  grid :: G
                kernel :: K
    field_dependencies :: F
            parameters :: P
                status :: S

    """
        KernelComputedField(loc, kernel, grid)

    Example
    =======

    ```julia
    using KernelAbstractions: @index, @kernel
    using Oceananigans.Fields: AveragedField, KernelComputedField, compute!
    using Oceananigans.Grids: Center, Face

    @inline ψ²(i, j, k, grid, ψ, Ψ) = @inbounds (ψ[i, j, k] - Ψ[i, j, k])^2
    @kernel function compute_variance!(var, grid, ϕ, Φ)
        i, j, k = @index(Global, NTuple)

        @inbounds var[i, j, k] = ψ′²(i, j, k, grid, ϕ, Φ)
    end

    u, v, w = model.velocities

    U = AveragedField(u, dims=(1, 2))
    V = AveragedField(V, dims=(1, 2))

    u′² = KernelComputedField(Face, Center, Center, compute_variance!, model; field_dependencies=(u, U,))
    v′² = KernelComputedField(Center, Face, Center, compute_variance!, model; field_dependencies=(v, V,))
    w′² = KernelComputedField(Center, Center, Face, compute_variance!, model; field_dependencies=(w, 0,))

    compute!(u′²)
    compute!(v′²)
    compute!(w′²)
    ```
    """
    function KernelComputedField{X, Y, Z}(kernel::K, arch, grid;
                                          field_dependencies = (),
                                          parameters::P = nothing,
                                          data = nothing,
                                          recompute_safely = true) where {X, Y, Z, K, P}

        field_dependencies = tupleit(field_dependencies)

        if isnothing(data)
            data = new_data(arch, grid, (X, Y, Z))
        end

        # Use FieldStatus if we want to avoid always recomputing
        status = recompute_safely ? nothing : FieldStatus(0.0)

        G = typeof(grid)
        A = typeof(data)
        S = typeof(status)
        F = typeof(field_dependencies)

        return new{X, Y, Z, S,
                   A, G, K, F, P}(data, grid, kernel, field_dependencies, parameters)
    end
end

KernelComputedField(X, Y, Z, kernel, model::AbstractModel; kwargs...) =
    KernelComputedField{X, Y, Z}(kernel, model.architecture, model.grid; kwargs...)

KernelComputedField(X, Y, Z, kernel, arch::AbstractArchitecture, grid::AbstractGrid; kwargs...) =
    KernelComputedField{X, Y, Z}(kernel, arch, grid; kwargs...)

function compute!(kcf::KernelComputedField{X, Y, Z}) where {X, Y, Z}

    for dependency in kcf.field_dependencies
        compute!(dependency)
    end

    arch = architecture(kcf.data)

    workgroup, worksize = work_layout(kcf.grid,
                                      :xyz,
                                      location=(X, Y, Z),
                                      include_right_boundaries=true)

    compute_kernel! = kcf.kernel(device(arch), workgroup, worksize)

    event = isnothing(kcf.parameters) ?
        compute_kernel!(kcf.data, kcf.grid, kcf.field_dependencies...) :
        compute_kernel!(kcf.data, kcf.grid, kcf.field_dependencies..., kcf.parameters)

    wait(device(arch), event)

    return nothing
end

compute!(field::KernelComputedField{X, Y, Z, <:FieldStatus}, time) where {X, Y, Z} =
    conditional_compute!(field, time)

Adapt.adapt_structure(to, kcf::KernelComputedField) = Adapt.adapt(to, kcf.data)


text:
module Day18

export get_inputs, get_solution1, get_solution2

## Input getting
function get_inputs()
    test_input1 = test_input2 = readlines(IOBuffer(
        """
        1 + 2 * 3 + 4 * 5 + 6
        1 + (2 * 3) + (4 * (5 + 6))
        2 * 3 + (4 * 5)
        5 + (8 * 3 + 9 + 3 * 4 * 3)
        5 * 9 * (7 * 3 * 3 + 9 * 3 + (8 + 6 * 4))
        ((2 + 4 * 9) * (6 + 9 * 8 + 6) + 6) + 2 + 4 * 2
        """))
    test_output1 = 71 + 51 + 26 + 437 + 12240 + 13632
    test_output2 = 231 + 51 + 46 + 1445 + 669060 + 23340
    data = readlines(joinpath(@__DIR__, "input.txt"))
    return (; test_input1, test_input2, test_output1, test_output2, data)
end

## Solution functions
# Define new addition operator from same precedence pool as *
const ⊠ = +

# And new multiplication operator from same precence pool as +
const ⊞ = *

# Replace the desired operation rules and parse as an Expr object to be evaled
parse_expr(str, rules...) = Meta.parse(reduce(replace, rules, init=str))

# Create expressions and eval them
get_solution(data, rules...) = sum(eval.(parse_expr.(data, rules...)))

# Part 1
get_solution1(data) = get_solution(data, '+'=>'⊠')

# Part 2
get_solution2(data) = get_solution(data, '+'=>'⊠', '*'=>'⊞')

end


text:
"""
    StaticSymbol

A statically typed `Symbol`.
"""
struct StaticSymbol{s}
    StaticSymbol{s}() where {s} = new{s::Symbol}()
    StaticSymbol(s::Symbol) = new{s}()
    StaticSymbol(x::StaticSymbol) = x
    StaticSymbol(x) = StaticSymbol(Symbol(x))
end
StaticSymbol(x, y) = StaticSymbol(Symbol(x, y))
StaticSymbol(x::StaticSymbol, y::StaticSymbol) = _cat_syms(x, y)
@generated function _cat_syms(::StaticSymbol{x}, ::StaticSymbol{y}) where {x,y}
    return :(StaticSymbol{$(QuoteNode(Symbol(x, y)))}())
end
StaticSymbol(x, y, z...) = StaticSymbol(StaticSymbol(x, y), z...)

Base.Symbol(::StaticSymbol{s}) where {s} = s::Symbol

Base.show(io::IO, ::StaticSymbol{s}) where {s} = print(io, "static(:$s)")


text:
import Base.-
import Base.*
import Base.+

mutable struct Conserved
    mass::Float64
    momentum::Float64
    energy::Float64
end

function +(argl::Conserved, argr::Conserved)
    return Conserved(argl.mass+argr.mass,
                     argl.momentum+argr.momentum,
                     argl.energy+argr.energy)
end

function *(s::Float64, c::Conserved)
    return Conserved(s*c.mass,
                     s*c.momentum,
                     s*c.energy)
end

function *(s::Int64, c::Conserved)
    return convert(Float64, s)*c
end

function -(argl::Conserved, argr::Conserved)
    return argl+(-1)*argr
end

max_stars_repo_name: bolverk/Sirion

max_issues_repo_name: bolverk/Sirion

max_forks_repo_name: bolverk/Sirion


text:
import Base: size, eltype
import RegularizedLeastSquares: initkaczmarz, dot_with_matrix_row, kaczmarz_update!

export reconstructionMultiPatch, MultiPatchOperator

# Necessary for Multi-System-Matrix FF reconstruction
voxelSize(bSF::MultiMPIFile) = voxelSize(bSF[1])
sfGradient(bSF::MultiMPIFile,dim) = sfGradient(bSF[1],dim)
generateHeaderDict(bSF::MultiMPIFile,bMeas::MPIFile) =
   generateHeaderDict(bSF[1],bMeas)

function reconstructionMultiPatch(bSF, bMeas::MPIFile;
  minFreq=0, maxFreq=1.25e6, SNRThresh=-1,maxMixingOrder=-1, numUsedFreqs=-1, sortBySNR=false, recChannels=1:numReceivers(bMeas), kargs...)

  freq = filterFrequencies(bSF,minFreq=minFreq, maxFreq=maxFreq,recChannels=recChannels, SNRThresh=SNRThresh, numUsedFreqs=numUsedFreqs, sortBySNR=sortBySNR)

  @debug "selecting $(length(freq)) frequencies"

  return reconstructionMultiPatch(bSF, bMeas, freq; kargs...)
end

function reconstructionMultiPatch(bSF, bMeas::MPIFile, freq;
            frames=nothing, bEmpty=nothing, nAverages=1, numAverages=nAverages,
	    spectralLeakageCorrection=true, voxelsize=voxelSize(bSF), kargs...)

  bgCorrection = (bEmpty != nothing)

  FFOp = MultiPatchOperatorHighLevel(bSF, bMeas, freq, bgCorrection;
                    kargs... )

  L = acqNumFGFrames(bMeas)
  (frames==nothing) && (frames=collect(1:L))
  nFrames=length(frames)

  uTotal_ = getMeasurementsFD(bMeas,frequencies=freq, frames=frames, numAverages=numAverages,
                             spectralLeakageCorrection=spectralLeakageCorrection)

  periodsSortedbyFFPos = unflattenOffsetFieldShift(ffPos(bMeas))
  uTotal = similar(uTotal_,size(uTotal_,1),length(periodsSortedbyFFPos),size(uTotal_,3))

  for k=1:length(periodsSortedbyFFPos)
      uTotal[:,k,:] = mean(uTotal_[:,periodsSortedbyFFPos[k],:], dims=2)
  end

  # Here we call a regular reconstruction function
  c_ = reconstruction(FFOp, uTotal; kargs...)
  c = reshape(c_, shape(FFOp.grid)..., :)

  # calculate axis
  shp = size(c)
  pixspacing = (voxelsize ./ sfGradient(bMeas,3) .* sfGradient(bSF,3)) * 1000u"mm"
  offset = (fieldOfViewCenter(FFOp.grid) .- 0.5.*fieldOfView(FFOp.grid) .+ 0.5.*spacing(FFOp.grid)) * 1000u"mm"

  # TODO does this provide the correct value in the multi-patch case?
  dtframes = acqNumAverages(bMeas)*dfCycle(bMeas)*numAverages*1u"s"
  # create image
  c = reshape(c,1,size(c)...)
  im = makeAxisArray(c, pixspacing, offset, dtframes)
  imMeta = ImageMeta(im,generateHeaderDict(bSF,bMeas))
  return imMeta
end

# MultiPatchOperator is a type that acts as the MPI system matrix but exploits
# its sparse structure.
# Its very important to keep this type typestable
mutable struct MultiPatchOperator{V<:AbstractMatrix, U<:Positions}
  S::Vector{V}
  grid::U
  N::Int
  M::Int
  RowToPatch::Vector{Int}
  xcc::Vector{Vector{Int}}
  xss::Vector{Vector{Int}}
  sign::Matrix{Int}
  nPatches::Int
  patchToSMIdx::Vector{Int}
end

eltype(FFOp::MultiPatchOperator) = eltype(FFOp.S[1])

function MultiPatchOperatorHighLevel(SF::MPIFile, bMeas, freq, bgCorrection::Bool; kargs...)
  return MultiPatchOperatorHighLevel(MultiMPIFile([SF]), bMeas, freq, bgCorrection; kargs...)
end

function MultiPatchOperatorHighLevel(bSF::MultiMPIFile, bMeas, freq, bgCorrection::Bool;
        FFPos = zeros(0,0), FFPosSF = zeros(0,0), kargs...)

  FFPos_ = ffPos(bMeas)

  periodsSortedbyFFPos = unflattenOffsetFieldShift(FFPos_)
  idxFirstPeriod = getindex.(periodsSortedbyFFPos,1)
  FFPos_ = FFPos_[:,idxFirstPeriod]

  if length(FFPos) > 0
    FFPos_[:] = FFPos
  end

  if length(FFPosSF) == 0
    L = length(ffPos(bSF[1]))
    FFPosSF_ = [vec(ffPos(SF))[l] for l=1:L, SF in bSF] #[vec(ffPos(SF)) for SF in bSF]
  else
    FFPosSF_ = FFPosSF #[vec(FFPosSF[:,l]) for l=1:size(FFPosSF,2)]
  end

  gradient = acqGradient(bMeas)[:,:,1,idxFirstPeriod]

  FFOp = MultiPatchOperator(bSF, bMeas, freq, bgCorrection,
                  FFPos = FFPos_,
                  gradient = gradient,
                  FFPosSF = FFPosSF_; kargs...)
  return FFOp
end

function MultiPatchOperator(SF::MPIFile, bMeas, freq, bgCorrection::Bool; kargs...)
  return MultiPatchOperator(MultiMPIFile([SF]), bMeas, freq, bgCorrection; kargs...)
end

function findNearestPatch(ffPosSF, FFPos, gradientSF, gradient)
  idx = -1
  minDist = 1e20
  for l = 1:size(ffPosSF,2)
    if gradientSF[l][:,:,1,1] == gradient
      dist = norm(ffPosSF[:,l].-FFPos)
      if dist < minDist
        minDist = dist
        idx = l
      end
    end
  end
  if idx < 0
    error("Something went wrong")
  end
  return idx
end

function MultiPatchOperator(SFs::MultiMPIFile, bMeas, freq, bgCorrection::Bool;
        mapping=zeros(0), kargs...)
  if length(mapping) > 0
    return MultiPatchOperatorExpliciteMapping(SFs,bMeas,freq,bgCorrection; mapping=mapping, kargs...)
  else
    return MultiPatchOperatorRegular(SFs,bMeas,freq,bgCorrection; kargs...)
  end
end

function MultiPatchOperatorExpliciteMapping(SFs::MultiMPIFile, bMeas, freq, bgCorrection::Bool;
                    denoiseWeight=0, FFPos=zeros(0,0), FFPosSF=zeros(0,0),
                    gradient=zeros(0,0,0),
                    roundPatches = false,
                    SFGridCenter = zeros(0,0),
                    systemMatrices = nothing,
                    mapping=zeros(0),
		    calibsize = nothing,
                    grid = nothing, kargs...)

  @debug "Loading System matrix"
  numPatches = size(FFPos,2)
  M = length(freq)
  RowToPatch = kron(collect(1:numPatches), ones(Int,M))

  if systemMatrices == nothing
    S = [getSF(SF,freq,nothing,"kaczmarz", bgCorrection=bgCorrection)[1] for SF in SFs]
  else
    S = systemMatrices
  end

  if length(SFGridCenter) == 0
    SFGridCenter = zeros(3,length(SFs))
    for l=1:length(SFs)
      SFGridCenter[:,l] = calibFovCenter(SFs[l])
    end
  end

  gradientSF = [acqGradient(SF) for SF in SFs]

  grids = RegularGridPositions[]
  patchToSMIdx = mapping

  sign = ones(Int, M, numPatches)

  # We first check which system matrix fits best to each patch. Here we use only
  # those system matrices where the gradient matches. If the gradient matches, we take
  # the system matrix with the closes focus field shift
  for k=1:numPatches
    idx = mapping[k]
    SF = SFs[idx]

    diffFFPos = FFPosSF[:,idx] .- FFPos[:,k]

    calibsizeTmp = (calibsize == nothing) ? calibSize(SF) : calibsize

    push!(grids, RegularGridPositions(calibsizeTmp,calibFov(SF),SFGridCenter[:,idx].-diffFFPos))
  end

  # We now know all the subgrids for each patch, if the corresponding system matrix would be taken as is
  # and if a possible focus field missmatch has been taken into account (by changing the center)
  @debug "Calculate Reconstruction Grid"
  if grid == nothing
    recoGrid = RegularGridPositions(grids)
  else
    recoGrid = grid
    # Test whether given grid covers all subgrids
    recoGridTest = RegularGridPositions(grids)
    testMin = round.(recoGridTest.center,digits=3) - round.(recoGridTest.fov./2,digits=3)
    testMax = round.(recoGridTest.center,digits=3) + round.(recoGridTest.fov./2,digits=3)
    gridMin = round.(recoGrid.center,digits=3) - round.(recoGrid.fov./2,digits=3)
    gridMax = round.(recoGrid.center,digits=3) + round.(recoGrid.fov./2,digits=3)
    if false in ((gridMin .<= testMin) .* (gridMax .>= testMax))
        @warn "A larger grid is used for reconstruction, since the given grid does not cover all subgrids: $recoGrid ⊅  $recoGridTest"
        recoGrid = recoGridTest
    end
  end

  # Within the next loop we will refine our grid since we now know our reconstruction grid
  for k=1:numPatches
    idx = mapping[k]
    SF = SFs[idx]

    issubgrid = isSubgrid(recoGrid,grids[k])
    if !issubgrid
      grids[k] = deriveSubgrid(recoGrid, grids[k])
    end
  end
  @debug "Use $(length(S)) patches"

  @debug "Calculate LUT"
  # now that we have all grids we can calculate the indices within the recoGrid
  xcc, xss = calculateLUT(grids, recoGrid)

  return MultiPatchOperator(S, recoGrid, length(recoGrid), M*numPatches,
             RowToPatch, xcc, xss, sign, numPatches, patchToSMIdx)
end

function MultiPatchOperatorRegular(SFs::MultiMPIFile, bMeas, freq, bgCorrection::Bool;
                    denoiseWeight=0, gradient=zeros(0,0,0),
                    roundPatches = false, FFPos=zeros(0,0), FFPosSF=zeros(0,0),
                    kargs...)

  @debug "Loading System matrix"
  numPatches = size(FFPos,2)
  M = length(freq)
  RowToPatch = kron(collect(1:numPatches), ones(Int,M))

  S = AbstractMatrix[]
  SOrigIdx = Int[]
  SIsPlain = Bool[]

  gradientSF = [acqGradient(SF) for SF in SFs]

  grids = RegularGridPositions[]
  matchingSMIdx = zeros(Int,numPatches)
  patchToSMIdx = zeros(Int,numPatches)

  # We first check which system matrix fits best to each patch. Here we use only
  # those system matrices where the gradient matches. If the gradient matches, we take
  # the system matrix with the closes focus field shift
  for k=1:numPatches
    idx = findNearestPatch(FFPosSF, FFPos[:,k], gradientSF, gradient[:,:,k])

    SF = SFs[idx]

    if isapprox(FFPosSF[:,idx],FFPos[:,k])
      diffFFPos = zeros(3)
    else
      diffFFPos = FFPosSF[:,idx] .- FFPos[:,k]
    end
    push!(grids, RegularGridPositions(calibSize(SF),calibFov(SF),calibFovCenter(SF).-diffFFPos))
    matchingSMIdx[k] = idx
  end

  # We now know all the subgrids for each patch, if the corresponding system matrix would be taken as is
  # and if a possible focus field missmatch has been taken into account (by changing the center)
  @debug "Calculate Reconstruction Grid"
  recoGrid = RegularGridPositions(grids)

  # Within the next loop we will refine our grid since we now know our reconstruction grid
  for k=1:numPatches
    idx = matchingSMIdx[k]
    SF = SFs[idx]

    issubgrid = isSubgrid(recoGrid,grids[k])
    if !issubgrid &&
       roundPatches &&
       spacing(recoGrid) == spacing(grids[k])

      issubgrid = true
      grids[k] = deriveSubgrid(recoGrid, grids[k])

    end

    # if the patch is a true subgrid we don't need to apply interpolation and can load the
    # matrix as is.
    if issubgrid
      # we first check if the matrix is already in memory
      u = -1
      for l=1:length(SOrigIdx)
        if SOrigIdx[l] == idx && SIsPlain[l]
          u = l
          break
        end
      end
      if u > 0 # its already in memory
        patchToSMIdx[k] = u
      else     # not yet in memory  -> load it
        S_, grid = getSF(SF,freq,nothing,"kaczmarz", bgCorrection=bgCorrection)
        push!(S,S_)
        push!(SOrigIdx,idx)
        push!(SIsPlain,true) # mark this as a plain system matrix (without interpolation)
        patchToSMIdx[k] = length(S)
      end
    else
      # in this case the patch grid does not fit onto the reco grid. Lets derive a subgrid
      # that is very similar to grids[k]
      newGrid = deriveSubgrid(recoGrid, grids[k])

      # load the matrix on the new subgrid
      S_, grid = getSF(SF,freq,nothing,"kaczmarz", bgCorrection=bgCorrection,
                   gridsize=shape(newGrid),
                   fov=fieldOfView(newGrid),
                   center=fieldOfViewCenter(newGrid).-fieldOfViewCenter(grids[k]))
                   # @TODO: I don't know the sign of aboves statement

      grids[k] = newGrid # we need to change the stored Grid since we now have a true subgrid
      push!(S,S_)
      push!(SOrigIdx,idx)
      push!(SIsPlain,false)
      patchToSMIdx[k] = length(S)
    end
  end
  @debug "Use $(length(S)) patches"

  @debug "Calculate LUT"
  # now that we have all grids we can calculate the indices within the recoGrid
  xcc, xss = calculateLUT(grids, recoGrid)

  sign = ones(Int, M, numPatches)

  return MultiPatchOperator(S, recoGrid, length(recoGrid), M*numPatches,
             RowToPatch, xcc, xss, sign, numPatches, patchToSMIdx)
end

function calculateLUT(grids, recoGrid)
  xss = Vector{Int}[]
  xcc = Vector{Int}[]
  for k=1:length(grids)
    N = length(grids[k])
    push!(xss, collect(1:N))
    xc = zeros(Int64,N)
    for n=1:N
      xc[n] = posToLinIdx(recoGrid,grids[k][n])
    end
    push!(xcc, xc)
  end
  return xcc, xss
end

function size(FFOp::MultiPatchOperator,i::Int)
  if i==2
    return FFOp.N
  elseif i==1
    return FFOp.M
  else
    error("bounds error")
  end
end

size(FFOp::MultiPatchOperator) = (FFOp.M,FFOp.N)

length(FFOp::MultiPatchOperator) = size(FFOp,1)*size(FFOp,2)

### The following is intended to use the standard kaczmarz method ###

function calculateTraceOfNormalMatrix(Op::MultiPatchOperator, weights)
  if length(Op.S) == 1
    trace = calculateTraceOfNormalMatrix(Op.S[1],weights)
    trace *= Op.nPatches #*prod(Op.PixelSizeSF)/prod(Op.PixelSizeC)
  else
    trace = sum([calculateTraceOfNormalMatrix(S,weights) for S in Op.S])
    #trace *= prod(Op.PixelSizeSF)/prod(Op.PixelSizeC)
  end
  return trace
end

setlambda(::MultiPatchOperator, ::Any) = nothing

function dot_with_matrix_row(Op::MultiPatchOperator, x::AbstractArray{T}, k::Integer) where T
  p = Op.RowToPatch[k]
  xs = Op.xss[p]
  xc = Op.xcc[p]

  j = mod1(k,div(Op.M,Op.nPatches))
  A = Op.S[Op.patchToSMIdx[p]]
  sign = Op.sign[j,Op.patchToSMIdx[p]]

  return dot_with_matrix_row_(A,x,xs,xc,j,sign)
end

function dot_with_matrix_row_(A::AbstractArray{T},x,xs,xc,j,sign) where T
  tmp = zero(T)
  @simd  for i = 1:length(xs)
     @inbounds tmp += sign*A[j,xs[i]]*x[xc[i]]
  end
  tmp
end

function kaczmarz_update!(Op::MultiPatchOperator, x::AbstractArray, k::Integer, beta)
  p = Op.RowToPatch[k]
  xs = Op.xss[p]
  xc = Op.xcc[p]

  j = mod1(k,div(Op.M,Op.nPatches))
  A = Op.S[Op.patchToSMIdx[p]]
  sign = Op.sign[j,Op.patchToSMIdx[p]]

  kaczmarz_update_!(A,x,beta,xs,xc,j,sign)
end

function kaczmarz_update_!(A,x,beta,xs,xc,j,sign)
  @simd for i = 1:length(xs)
    @inbounds x[xc[i]] += beta* conj(sign*A[j,xs[i]])
  end
end

function initkaczmarz(Op::MultiPatchOperator,λ,weights::Vector)
  T = typeof(real(Op.S[1][1]))
  denom = T[] #zeros(T,Op.M)
  rowindex = Int64[] #zeros(Int64,Op.M)

  MSub = div(Op.M,Op.nPatches)

  if length(Op.S) == 1
    for i=1:MSub
      s² = rownorm²(Op.S[1],i)*weights[i]^2
      if s²>0
        for l=1:Op.nPatches
          k = i+MSub*(l-1)
          push!(denom,weights[i]^2/(s²+λ)) #denom[k] = weights[i]^2/(s²+λ)
          push!(rowindex,k) #rowindex[k] = k
        end
      end
    end
  else
    for l=1:Op.nPatches
      for i=1:MSub
        s² = rownorm²(Op.S[Op.patchToSMIdx[l]],i)*weights[i]^2
        if s²>0
          k = i+MSub*(l-1)
          push!(denom,weights[i]^2/(s²+λ)) #denom[k] = weights[i]^2/(s²+λ)
          push!(rowindex,k) #rowindex[k] = k
        end
      end
    end
  end

  denom, rowindex
end


text:
module Smearing
    using LinearAlgebra, SpecialFunctions
    import ..LTK_universe:Universe,calc_gaugeforce!,expF_U!
    import ..LieAlgebrafields:LieAlgebraFields,clear!,add!
    import ..Gaugefields:GaugeFields,normalize!,SU,evaluate_wilson_loops!,TA!,set_wing!
    import ..Wilsonloops:Wilson_loop_set,calc_coordinate,make_plaq_staple_prime,calc_shift,make_plaq,make_plaq_staple
    import ..Actions:SmearingParam_single,SmearingParam_multi

    function add!(a::Array{N,1},α,b::Array{N,1}) where N <: LieAlgebraFields
        for μ=1:4
            add!(a[μ],α,b[μ]) # a_\mu+= \beta b_\mu
        end
    end

    function gradientflow!(U::Array{N,1},univ::Universe,tempW1,tempW2,Nflow::Int = 1,eps::Float64 = 0.01) where N <: GaugeFields
        # Here we use definition in 1006.4518 except for the hermiticity.
        # ref https://www2.physik.uni-bielefeld.de/fileadmin/user_upload/theory_e6/Master_Theses/Masterarbeit_LukasMazur.pdf
        NC = univ.NC
        NX = univ.NX
        NY = univ.NY
        NZ = univ.NZ
        NT = univ.NT

        W1 = tempW1
        W2 = tempW2

        F0 = Array{LieAlgebraFields,1}(undef,4)
        for μ=1:4
            F0[μ] = LieAlgebraFields(NC,NX,NY,NZ,NT)
        end
        F1 = Array{LieAlgebraFields,1}(undef,4)
        for μ=1:4
            F1[μ] = LieAlgebraFields(NC,NX,NY,NZ,NT)
        end
        F2 = Array{LieAlgebraFields,1}(undef,4)
        for μ=1:4
            F2[μ] = LieAlgebraFields(NC,NX,NY,NZ,NT)
        end
        Ftmp = Array{LieAlgebraFields,1}(undef,4)
        for μ=1:4
            Ftmp[μ] = LieAlgebraFields(NC,NX,NY,NZ,NT)
        end

        for istep=1:Nflow #RK4 integrator
            calc_gaugeforce!(F0,U,univ) #F
            expF_U!(W1,F0,-eps*(1/4),univ)  #exp(eps*F)*U
            #
            calc_gaugeforce!(F1,W1,univ) #F
            clear!(Ftmp)
            add!(Ftmp,-(8/9*eps),F1)
            add!(Ftmp,(17/36*eps),F0)
            expF_U!(W2,Ftmp,1,univ)
            #
            calc_gaugeforce!(F2,W2,univ) #F
            clear!(Ftmp)
            add!(Ftmp,-(3/4*eps),F2)
            add!(Ftmp,(8/9*eps),F1)
            add!(Ftmp,-(17/36*eps),F0)
            expF_U!(U,Ftmp,1,univ)
            #
        end
        #return W1 #test
        #return W0
    end

    function gradientflow!(U::Array{N,1},univ::Universe,Nflow::Int = 1,eps::Float64 = 0.01) where N <: GaugeFields
        W1 = deepcopy(U)
        W2 = deepcopy(U)
        gradientflow!(U,univ,W1,W2,Nflow,eps)
        return
    end
    #=
    function gradientflow_euler!(univ::Universe) #backup
        eps = 0.01
        Nflow = 1
        NC = univ.NC
        NX = univ.NX
        NY = univ.NY
        NZ = univ.NZ
        NT = univ.NT

        F = Array{LieAlgebraFields,1}(undef,4)
        for μ=1:4
            F[μ] = LieAlgebraFields(NC,NX,NY,NZ,NT)
        end

        U = similar(univ.U)

        for istep=1:Nflow #Euler integrator
            calc_gaugeforce!(F,univ) #F
            expF_U!(U,F,eps,univ)  #exp(eps*F)*U
        end
    end
    =#

    function calc_stout(U::Array{GaugeFields{SU{NC}},1},ρ) where NC
        Uout = similar(U)
        calc_stout!(Uout,U,ρ)
        return Uout
    end

    function calc_stout!(Uout::Array{GaugeFields{SU{NC}},1},ρ) where NC
        Uin = deepcopy(Uout)
        calc_stout!(Uout,Uin,ρ)
    end

    function calc_stout!(Uout::Array{GaugeFields{SU{NC}},1},U::Array{GaugeFields{SU{NC}},1},ρ) where NC
        @assert Uout != U "input U and output U should not be same!"

        NX = U[1].NX
        NY = U[1].NY
        NZ = U[1].NZ
        NT = U[1].NT
        Q = zeros(ComplexF64,NC,NC)
        C = zeros(ComplexF64,NC,NC)

        for μ=1:4
            loops = make_plaq_staple_prime(μ)
            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        for ix=1:NX
                            C .= 0
                            evaluate_wilson_loops!(C,loops,U,ix,iy,iz,it)
                            #display(C)
                            #println("\t")

                            Ω = ρ*C[:,:]*U[μ][:,:,ix,iy,iz,it]'
                            Q = (im/2)*(Ω' .- Ω) .- (im/(2NC))*tr(Ω' .- Ω)
                            Uout[μ][:,:,ix,iy,iz,it] = exp(im*Q)*U[μ][:,:,ix,iy,iz,it]
                            set_wing!(Uout[μ],ix,iy,iz,it)

                        end
                    end
                end
            end
        end
        #display(Uout[1][:,:,1,1,1,1])
        return
    end

    function calc_multihit!(Uout::Array{GaugeFields{SU{NC}},1},U::Array{GaugeFields{SU{NC}},1},β) where NC
        @assert Uout != U "input U and output U should not be same!"
        @assert U[1].NC == 2 "NC=2 is required for now but $(U[1].NC)"

        NX = U[1].NX
        NY = U[1].NY
        NZ = U[1].NZ
        NT = U[1].NT
        Q = zeros(ComplexF64,NC,NC)
        A = zeros(ComplexF64,NC,NC)

        for μ=1:4
            loops = make_plaq_staple_prime(μ)
            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        for ix=1:NX
                            A .= 0
                            evaluate_wilson_loops!(A,loops,U,ix,iy,iz,it)
                            K = sqrt(abs(det(A)))

                            #display(C)
                            #println("\t")

                            #Ω = ρ*C[:,:]*U[μ][:,:,ix,iy,iz,it]'
                            #Q = (im/2)*(Ω' .- Ω) .- (im/(2NC))*tr(Ω' .- Ω)
                            Uout[μ][:,:,ix,iy,iz,it] = K*inv(A)*besseli(2, β*K)/besseli(1, β*K)
                            #0  #exp(im*Q)*U[μ][:,:,ix,iy,iz,it]
                            set_wing!(Uout[μ],ix,iy,iz,it)

                        end
                    end
                end
            end
        end
        #display(Uout[1][:,:,1,1,1,1])
        return
    end

    function calc_fatlink_APE!(Uout::Array{GaugeFields{SU{NC}},1},α,β;normalize_method= "special unitary",temporal_dir_smear=true) where NC
        Uin = deepcopy(Uout)
        calc_fatlink_APE!(Uout,Uin,α,β,normalize_method=normalize_method,temporal_dir_smear=temporal_dir_smear)
    end

    function calc_fatlink_APE!(Uout::Array{GaugeFields{SU{NC}},1},U::Array{GaugeFields{SU{NC}},1},α,β;normalize_method= "special unitary",temporal_dir_smear=true) where NC
        @assert Uout != U "input U and output U should not be same!"
        NX = U[1].NX
        NY = U[1].NY
        NZ = U[1].NZ
        NT = U[1].NT
        Vtemp = zeros(ComplexF64,NC,NC)
        nmethod(A) = if normalize_method == "unitary"
                        A*(A'*A)^(-1/2)#inv(sqrt(A'*A))
                    elseif normalize_method == "special unitary"
                        normalize!(A)
                        A
                    else
                        error("normalize_method should be unitary or special unitary. Now, $normalize_method")
                    end

        #println(normalize_method)
        d=3
        if temporal_dir_smear
            d=4
        end
        for μ=1:d
            loops = make_plaq_staple_prime(μ)
            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        for ix=1:NX
                            Vtemp .= 0
                            evaluate_wilson_loops!(Vtemp,loops,U,ix,iy,iz,it)
                            #Vtemp = Vtemp'
                            Vtemp[:,:] = (1-α)*U[μ][:,:,ix,iy,iz,it] .+ (β/6)*Vtemp[:,:]
                            #nmethod!(Vtemp)
                            Uout[μ][:,:,ix,iy,iz,it] = nmethod(Vtemp)#Vtemp[:,:]
                            set_wing!(Uout[μ],ix,iy,iz,it)
                        end
                    end
                end
            end
        end
        return
    end

    function calc_fatlink_APE(U::Array{GaugeFields{SU{NC}},1},α,β;normalize_method= "special unitary",temporal_dir_smear=true) where NC
        Uout = similar(U)
        calc_fatlink_APE!(Uout,U,α,β,normalize_method=normalize_method,temporal_dir_smear=temporal_dir_smear)
        return Uout
    end

    function calc_multihit!(Uout::Array{GaugeFields{SU{NC}},1},U::Array{GaugeFields{SU{NC}},1},β) where NC
        @assert Uout != U "input U and output U should not be same!"
        @assert U[1].NC == 2 "NC=2 is required for now but $(U[1].NC)"

        NX = U[1].NX
        NY = U[1].NY
        NZ = U[1].NZ
        NT = U[1].NT
        Q = zeros(ComplexF64,NC,NC)
        A = zeros(ComplexF64,NC,NC)

        for μ=1:4
            loops = make_plaq_staple_prime(μ)
            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        for ix=1:NX
                            A .= 0
                            evaluate_wilson_loops!(A,loops,U,ix,iy,iz,it)
                            K = sqrt(abs(det(A)))

                            #display(C)
                            #println("\t")

                            #Ω = ρ*C[:,:]*U[μ][:,:,ix,iy,iz,it]'
                            #Q = (im/2)*(Ω' .- Ω) .- (im/(2NC))*tr(Ω' .- Ω)
                            Uout[μ][:,:,ix,iy,iz,it] = K*inv(A)*besseli(2, β*K)/besseli(1, β*K)
                            #0  #exp(im*Q)*U[μ][:,:,ix,iy,iz,it]
                            set_wing!(Uout[μ],ix,iy,iz,it)

                        end
                    end
                end
            end
        end
        #display(Uout[1][:,:,1,1,1,1])
        return
    end
end


text:
@testset "DisjointSet" begin

    @testset "IntDisjointSet" begin
        for T in [Int, UInt8, Int8, UInt16, Int16, UInt32, Int32, UInt64]
            @testset "eltype = $(T)" begin
                s = IntDisjointSet(T(10))
                s2 = IntDisjointSet{T}(10)

                @testset "basic tests" begin
                    @test length(s) == 10
                    @test length(s2) == 10
                    @test eltype(s) == T
                    @test eltype(s2) == T
                    @test eltype(typeof(s)) == T
                    @test eltype(typeof(s2)) == T
                    @test num_groups(s) == T(10)
                    @test num_groups(s2) == T(10)

                    for i = 1:10
                        @test find_root!(s, T(i)) == T(i)
                    end
                    @test_throws BoundsError find_root!(s, T(11))

                    @test !in_same_set(s, T(2), T(3))
                end

                @testset "union!" begin
                    union!(s, T(2), T(3))
                    @test num_groups(s) == T(9)
                    @test in_same_set(s, T(2), T(3))
                    @test find_root!(s, T(3)) == T(2)
                    union!(s, T(3), T(2))
                    @test num_groups(s) == T(9)
                    @test in_same_set(s, T(2), T(3))
                    @test find_root!(s, T(3)) == T(2)
                end

                @testset "more tests" begin
                    # We cannot support arbitrary indexing and still use @inbounds with IntDisjointSet
                    # (and it's not useful anyway)
                    @test_throws MethodError push!(s, T(22))

                    @test push!(s) == T(11)
                    @test num_groups(s) == T(10)

                    @test union!(s, T(8), T(7)) == T(8)
                    @test union!(s, T(5), T(6)) == T(5)
                    @test union!(s, T(8), T(5)) == T(8)
                    @test num_groups(s) == T(7)
                    @test find_root!(s, T(6)) == T(8)
                    union!(s, T(2), T(6))
                    @test find_root!(s, T(2)) == T(8)
                    root1 = find_root!(s, T(6))
                    root2 = find_root!(s, T(2))
                    @test root_union!(s, T(root1), T(root2)) == T(8)
                    @test union!(s, T(5), T(6)) == T(8)
                end
            end
        end
    end

    @testset "IntDisjointSet overflow" begin
        for T in [UInt8, Int8]
            s = IntDisjointSet(T(typemax(T)-1))
            push!(s)
            @test_throws ArgumentError push!(s)
        end
    end

    @testset "DisjointSet" begin
        # DisjointSet supports arbitrary indices
        s = DisjointSet{Int}(1:10)

        @testset "constructor" begin
            @test DisjointSet() isa DisjointSet{Any}
            @test DisjointSet{Int}(1:10) isa DisjointSet{Int}
            @test DisjointSet(collect(1:10)) isa DisjointSet{Int}
            @test DisjointSet(x*im for x = 1:10) isa DisjointSet{Complex{Int}}
            g = (x % 2 == 0 ? x+1//x : x*im for x = 1:10)
            @test DisjointSet(g) isa DisjointSet{Number}
            @test DisjointSet(Any[1,2,3]) isa DisjointSet
            @test DisjointSet{Any}(Any[1,2,3]) isa DisjointSet{Any}
            @test DisjointSet{Int}(Any[1,2,3]) isa DisjointSet{Int}
        end

        @testset "basic tests" begin
            @test length(s) == 10
            @test num_groups(s) == 10
            @test eltype(s) == Int
            @test eltype(typeof(s)) == Int
            @test length(empty(s)) == num_groups(empty(s)) == 0
            @test empty(s) isa DisjointSet{eltype(s)}
            @test collect(s) == collect(1:10)
            g = (x % 2 == 0 ? x+1//x : x*im for x = 1:10)
            s1 = DisjointSet(g)
            @test length(s1) == 10
            @test sizehint!(s1, 100) === s1

            r = [find_root!(s, i) for i in 1 : 10]
            @test isequal(r, collect(1:10))
        end

        @testset "union!" begin
            for i = 1 : 5
                x = 2 * i - 1
                y = 2 * i
                union!(s, x, y)
                @test find_root!(s, x) == find_root!(s, y)
            end

            @test union!(s, 1, 4) == find_root!(s, 1)
            @test union!(s, 3, 5) == find_root!(s, 1)
            @test union!(s, 7, 9) == find_root!(s, 7)

            @test length(s) == 10
            @test num_groups(s) == 2
        end

        @testset "r0" begin
            r0 = [ find_root!(s,i) for i in 1:10 ]
            # Since this is a DisjointSet (not IntDisjointSet), the root for 17 will be 17, not 11
            push!(s, 17)

            @test length(s) == 11
            @test num_groups(s) == 3

            r0 = [ r0 ; 17]
            r = [find_root!(s, i) for i in [1 : 10; 17] ]
            @test isequal(r, r0)
        end

        @testset "root_union!" begin
            root1 = find_root!(s, 7)
            root2 = find_root!(s, 3)
            @test root1 != root2
            root_union!(s, 7, 3)
            @test find_root!(s, 7) == find_root!(s, 3)
        end

        @testset "Some tests using non-integer disjoint sets" begin
            elems = ["a", "b", "c", "d"]
            a = DisjointSet{AbstractString}(elems)
            @test collect(a) == ["a", "b", "c", "d"]
            union!(a, "a", "b")
            @test in_same_set(a,"a","b")
            @test find_root!(a,"a") == find_root!(a,"b")
            @test find_root!(a,"a") in elems
            @test !in_same_set(a, "c", "d")
            # union returns new root
            @test find_root!(a,"a") == union!(a,"b","c")
            union!(a,"c","d")
            # Now they should be in same set, and a is transitively connected to d
            @test in_same_set(a,"a", "d")
            # Root element should thus be same for all:
            @test all(find_root!(a,first(elems)) .== map(x->find_root!(a,x),elems))

            #@test_throws KeyError find_root!(a,"f")

            push!(a, "f")
            @test find_root!(a,"a") != find_root!(a,"f")
        end
    end

    @testset "Test deprecations from #700" begin
        @test DisjointSets(1:10) isa DisjointSet
        @test IntDisjointSets(10) isa IntDisjointSet
        @test DisjointSets(1:10...) == DisjointSet(1:10)
        @test typeof(DisjointSet(1:10)) <: DisjointSets
    end
end # @testset DisjointSet


text:
files = [
    "001.jl",
    "002.jl",
    "003.jl",
    "004.jl",
    "005.jl",
    "006.jl",
    "007.jl",
    "008.jl",
    "009.jl",
    "010.jl",
    "011.jl",
    "012.jl",
]

for file in files
    path = joinpath("julia", file)
    cmd = `julia-0.4 $path`
    println(cmd)
    run(cmd)
end


text:
using Pkg
# change to course environment
Pkg.activate(dirname(@__FILE__))
# download packages if required
Pkg.instantiate()


text:
# MIT License

# Copyright (c) Microsoft Corporation.

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE

# small offset along ray direction to prevent ray from reintersecting the surface it just intersected. Optical pathlength is measured in mm so this is much smaller than a wavelength of light.
const RAY_OFFSET = 1e-6
# if the power of a ray drops below this value return nothing and do not continue the ray tracing recursion
const POWER_THRESHOLD = 1e-6
# if the ray is traced recursively more than this many times then the trace will be terminated and the ray ignored (mostly happens with TIR)
const TRACE_RECURSION_LIMIT = 256

# triangulated surfaces are imprecise mostly due to floating point precision, particularly if using trig functions
# expanding every triangle by a fraction resolves missed intersections
const TRIANGULATION_EXPANSION = 2 * eps(Float64)
# we also need to extend the bounds of the surface to ensure we don't miss anythgin at the edges of the surface
const TRIANGULATION_EXTENSION = 0.02

# for accelerated surface the triangulated surface can often be very imprecise, this will mean rays which start close
# to the surface might miss intersection because the ray origin lies above the triangulated surface but below the true surface
# (where for this example the ray direction is straight upwards), offsetting the ray fixes this, and as we check the alpha
# of the intersection point afterwards there are no problems
const ACCEL_SURF_RAY_OFFSET = 0.1

# settings for visualization
const MIN_VIS_TRI_SIZE = 1e-3
const MESH_PRECISION = 1e-3
const VIS_RECURSION_LIMIT = 64

# number of M/N to precompute for various variables involved in the QType calculations
const QTYPE_PRECOMP = 16

# maximum number of orders for a thin grating interface
const GRATING_MAX_ORDERS = 10

# LensAssemblyN methods will be generated at compile time (faster on first run) for N up to this number
const PREGENERATED_LENS_ASSEMBLY_SIZE = 10