text:
# common functors

abstract Functor{N}  # N is the number of arguments

typealias UnaryFunctor Functor{1}
typealias BinaryFunctor Functor{2}
typealias TernaryFunctor Functor{3}

#######################################
#
#  Macros to define functors
#
#######################################

macro functor1(F, T, f)
    quote
        global evaluate
        type $F <: Functor{1} end
        evaluate(::($F), x::($T)) = ($f)(x)
    end
end

macro functor2(F, T, f)
    quote
        global evaluate
        type $F <: Functor{2} end
        evaluate(::($F), x::($T), y::($T)) = ($f)(x, y)
    end
end

macro functor3(F, T, f)
    quote
        global evaluate
        type $F <: Functor{3} end
        evaluate(::($F), x::($T), y::($T), z::($T)) = ($f)(x, y, z)
    end
end

#######################################
#
#  Functors for operators
#
#######################################

export
    Negate, Add, Subtract, Multiply, Divide, Pow, Modulo,
    Greater, GreaterEqual, Less, LessEqual, Equal, NotEqual,
    Not, And, Or,
    BitwiseNot, BitwiseAnd, BitwiseOr, BitwiseXor

# arithmetic operators

type Negate <: Functor{1} end
evaluate(::Negate, x::Number) = -x

type Add <: Functor{2} end
evaluate(::Add, x::Number, y::Number) = x + y

type Subtract <: Functor{2} end
evaluate(::Subtract, x::Number, y::Number) = x - y

type Multiply <: Functor{2} end
evaluate(::Multiply, x::Number, y::Number) = x * y

type Divide <: Functor{2} end
evaluate(::Divide, x::Number, y::Number) = x / y

type Pow <: Functor{2} end
evaluate(::Pow, x::Number, y::Number) = x ^ y

type Modulo <: Functor{2} end
evaluate(::Modulo, x::Real, y::Real) = x % y

# comparison operators

type Greater <: Functor{2} end
evaluate(::Greater, x::Real, y::Real) = x > y

type GreaterEqual <: Functor{2} end
evaluate(::GreaterEqual, x::Real, y::Real) = x >= y

type Less <: Functor{2} end
evaluate(::Less, x::Real, y::Real) = x < y

type LessEqual <: Functor{2} end
evaluate(::LessEqual, x::Real, y::Real) = x <= y

type Equal <: Functor{2} end
evaluate(::Equal, x::Number, y::Number) = x == y

type NotEqual <: Functor{2} end
evaluate(::NotEqual, x::Number, y::Number) = x != y

# logical operators

type Not <: Functor{1} end
evaluate(::Not, x::Bool) = !x

type And <: Functor{2} end
evaluate(::And, x::Bool, y::Bool) = (x && y)

type Or <: Functor{2} end
evaluate(::Or, x::Bool, y::Bool) = (x || y)

type BitwiseNot <: Functor{1} end
evaluate(::BitwiseNot, x::Real) = (~x)

type BitwiseAnd <: Functor{2} end
evaluate(::BitwiseAnd, x::Real, y::Real) = (x & y)

type BitwiseOr <: Functor{2} end
evaluate(::BitwiseOr, x::Real, y::Real) = (x | y)

type BitwiseXor <: Functor{2} end
evaluate(::BitwiseXor, x::Real, y::Real) = (x $ y)

#######################################
#
#  Functors for elementary math
#
#######################################

export
    MaxFun, MinFun,
    AbsFun, Abs2Fun, SqrFun, RcpFun, SqrtFun, RsqrtFun,
    CbrtFun, RcbrtFun, HypotFun, FixPow, FixAbsPow,
    IsfiniteFun, IsinfFun, IsnanFun, IsequalFun,
    FloorFun, CeilFun, RoundFun, TruncFun,
    IfloorFun, IceilFun, IroundFun, ItruncFun,
    ModFun, FldFun, RemFun, DivFun,
    ExpFun, Exp2Fun, Exp10Fun, Expm1Fun,
    LogFun, Log2Fun, Log10Fun, Log1pFun, LdexpFun,
    SinFun, CosFun, TanFun, CotFun, SecFun, CscFun,
    SindFun, CosdFun, TandFun, CotdFun, SecdFun, CscdFun,
    SinhFun, CoshFun, TanhFun, CothFun, SechFun, CschFun,
    AsinFun, AcosFun, AtanFun, AcotFun, AsecFun, AcscFun,
    AsindFun, AcosdFun, AtandFun, AcotdFun, AsecdFun, AcscdFun,
    AsinhFun, AcoshFun, AtanhFun, AcothFun, AsechFun, AcschFun,
    Atan2Fun

# basic functors

@functor2 MaxFun Real max
@functor2 MinFun Real min

type _Max <: Functor{2} end
evaluate(::_Max, x::Real, y::Real) = ifelse(y > x, y, x)

type NonnegMax <: Functor{2} end
evaluate(::NonnegMax, x::Real, y::Real) = ifelse(y > x, y, x)

type _Min <: Functor{2} end
evaluate(::_Min, x::Real, y::Real) = ifelse(y < x, y, x)

# absolute value & power

@functor1 AbsFun   Number abs
@functor1 Abs2Fun  Number abs2
@functor1 SqrFun   Number sqr
@functor1 RcpFun   Number rcp
@functor1 SqrtFun  Number sqrt
@functor1 RsqrtFun Number rsqrt
@functor1 CbrtFun  Real   cbrt
@functor1 RcbrtFun Real   rcbrt
@functor2 HypotFun Real   hypot

immutable FixPow{T<:Real} <: Functor{1}
    p::T
end
evaluate(f::FixPow, x::Real) = (x ^ f.p)

immutable FixAbsPow{T<:Real} <: Functor{1}
    p::T
end
evaluate(f::FixAbsPow, x::Real) = (abs(x) ^ f.p)

# number classifying functors

@functor1 IsfiniteFun Real   isfinite
@functor1 IsnanFun    Real   isnan
@functor1 IsinfFun    Real   isinf
@functor2 IsequalFun  Number isequal

# rounding functors

@functor1 FloorFun Real floor
@functor1 CeilFun  Real ceil
@functor1 RoundFun Real round
@functor1 TruncFun Real trunc

@functor1 IfloorFun Real ifloor
@functor1 IceilFun  Real iceil
@functor1 IroundFun Real iround
@functor1 ItruncFun Real itrunc

# division & modulus functors

@functor2 ModFun Real mod
@functor2 FldFun Real fld
@functor2 RemFun Real rem
@functor2 DivFun Real div

# exponentiation & logarithm

@functor1 ExpFun   Number exp
@functor1 Exp2Fun  Number exp2
@functor1 Exp10Fun Number exp10

@functor1 LogFun   Number log
@functor1 Log2Fun  Number log2
@functor1 Log10Fun Number log10

@functor1 Expm1Fun Real expm1
@functor1 Log1pFun Real log1p

type LdexpFun <: Functor{2} end
evaluate(::LdexpFun, x::Real, n::Integer) = ldexp(x, n)

# trigonometric & hyperbolic

@functor1 SinFun Number sin
@functor1 CosFun Number cos
@functor1 TanFun Number tan
@functor1 CotFun Number cot
@functor1 SecFun Number sec
@functor1 CscFun Number csc

@functor1 SindFun Real sind
@functor1 CosdFun Real cosd
@functor1 TandFun Real tand
@functor1 CotdFun Real cotd
@functor1 SecdFun Real secd
@functor1 CscdFun Real cscd

@functor1 SinhFun Number sinh
@functor1 CoshFun Number cosh
@functor1 TanhFun Number tanh
@functor1 CothFun Number coth
@functor1 SechFun Number sech
@functor1 CschFun Number csch

@functor1 AsinFun Number asin
@functor1 AcosFun Number acos
@functor1 AtanFun Number atan
@functor1 AcotFun Number acot
@functor1 AsecFun Number asec
@functor1 AcscFun Number acsc

@functor1 AsindFun Real asind
@functor1 AcosdFun Real acosd
@functor1 AtandFun Real atand
@functor1 AcotdFun Real acotd
@functor1 AsecdFun Real asecd
@functor1 AcscdFun Real acscd

@functor1 AsinhFun Number asinh
@functor1 AcoshFun Number acosh
@functor1 AtanhFun Number atanh
@functor1 AcothFun Number acoth
@functor1 AsechFun Number asech
@functor1 AcschFun Number acsch

@functor2 Atan2Fun Real atan2

#######################################
#
#  Functors for special functions
#
#######################################

export ErfFun, ErfcFun, ErfInvFun, ErfcInvFun,
    GammaFun, LgammaFun, LfactFun, DigammaFun,
    BetaFun, LbetaFun, ZetaFun,
    AiryFun, AiryprimeFun, AiryaiFun, AiryaiprimeFun,
    AirybiFun, AirybiprimeFun,
    BesseljFun, Besselj0Fun, Besselj1Fun, BesseliFun, BesselkFun,
    LogitFun, LogisticFun, InvLogisticFun,
    XlogxFun, XlogyFun

# error functors

@functor1 ErfFun     Real erf
@functor1 ErfcFun    Real erfc
@functor1 ErfInvFun  Real erfinv
@functor1 ErfcInvFun Real erfcinv

# gamma functors

@functor1 GammaFun   Real gamma
@functor1 LgammaFun  Real lgamma
@functor1 LfactFun   Real lfact
@functor1 DigammaFun Real digamma

# beta functors

@functor2 BetaFun  Real beta
@functor2 LbetaFun Real lbeta
@functor1 EtaFun   Real eta
@functor1 ZetaFun  Real zeta

# airy functors

@functor1 AiryFun        Real airy
@functor1 AiryprimeFun   Real airyprime
@functor1 AiryaiFun      Real airyai
@functor1 AiryaiprimeFun Real airyaiprime
@functor1 AirybiFun      Real airybi
@functor1 AirybiprimeFun Real airybiprime

# bessel functors

@functor2 BesseljFun  Real besselj
@functor1 Besselj0Fun Real besselj0
@functor1 Besselj1Fun Real besselj1
@functor2 BesseliFun  Real besseli
@functor2 BesselkFun  Real besselk

# stats-related functors

@functor1 LogitFun       Real logit
@functor1 LogisticFun    Real logistic
@functor1 InvLogisticFun Real invlogistic

@functor1 XlogxFun Real xlogx
@functor2 XlogyFun Real xlogy

#######################################
#
#  Ternary functors
#
#######################################

export FMA, IfelseFun

type FMA <: Functor{3} end
evaluate(::FMA, x::Number, y::Number, z::Number) = (x + y * z)

type IfelseFun <: Functor{3} end
evaluate{T<:Number}(::IfelseFun, c::Bool, x::T, y::T) = ifelse(c, x, y)

#######################################
#
#  type inference facility
#
#######################################

fptype{T<:FloatingPoint}(::Type{T}) = T
fptype{T<:Integer}(::Type{T}) = Float64
fptype{T<:Union(Bool,Uint8,Int8,Uint16,Int16)}(::Type{T}) = Float32

typealias BasicBinaryFuns Union(Add,Subtract,Multiply,MaxFun,MinFun,_Max,_Min,NonnegMax)

result_type{T}(f::BasicBinaryFuns, ::Type{T}, ::Type{T}) = promote_type(T)
result_type{T1,T2}(f::BasicBinaryFuns, ::Type{T1}, ::Type{T2}) = promote_type(T1, T2)

result_type(f::Functor{1}, t1::Type) = typeof(evaluate(f, one(t1)))
result_type(f::Functor{2}, t1::Type, t2::Type) = typeof(evaluate(f, one(t1), one(t2)))
result_type(f::Functor{3}, t1::Type, t2::Type, t3::Type) = typeof(evaluate(f, one(t1), one(t2), one(t3)))

typealias BlasFP Union(Float32,Float64)


text:
module robot

struct Point
    x::Int
    y::Int
end

const Black = 0
const White = 1
const Color = Int

const Counterclockwise = 0
const Clockwise = 1
const Rotation = Int

const Up = Point(0, 1)
const Down = Point(0, -1)
const Left = Point(-1, 0)
const Right = Point(1, 0)

mutable struct Robot
    orientation::Array{Point,1}
    location::Point
    surface::Dict{Point,Color}
end

function deploy(surface::Dict{Point,Color})
    return Robot([Up, Left, Down, Right], Point(0, 0), surface)
end

function deploy(surface::Dict{Point,Color}, point::Point)
    return Robot([Up, Left, Down, Right], point, surface)
end

function color(robot::Robot)
    ## anything unpainted is black
    return haskey(robot.surface, robot.location) ?
        robot.surface[robot.location] :
        Black
end

function paint!(robot::Robot, color::Color)
    robot.surface[robot.location] = color
end

function move!(robot::Robot, rotation::Rotation)
    if rotation == Counterclockwise
        reverse!(robot.orientation)
        prepend!(robot.orientation, [pop!(robot.orientation)])
        reverse!(robot.orientation)
    else
        prepend!(robot.orientation, [pop!(robot.orientation)])
    end
    orientation = first(robot.orientation)
    robot.location = Point(robot.location.x + orientation.x, robot.location.y + orientation.y)
    l = robot.location
end

end

ext: jl


text:
# Programa de Calculadora
# ITO - 22 oct 2019

using Gtk
ENV["GTK_CSD"] = 0
cal = GtkWindow()
global status = 0
set_gtk_property!(cal,:height_request,500)
set_gtk_property!(cal,:width_request,500)
set_gtk_property!(cal,:title,"Calculadora 0.0.1")
set_gtk_property!(cal,:window_position, 3)

bt0 = GtkButton("0")
bt1 = GtkButton("1")
bt2 = GtkButton("2")
bt3 = GtkButton("3")
bt4 = GtkButton("4")
bt5 = GtkButton("5")
bt6 = GtkButton("6")
bt7 = GtkButton("7")
bt8 = GtkButton("8")
bt9 = GtkButton("9")
btOFF = GtkButton("OFF")
btMinus = GtkButton("-")
btPlus = GtkButton("+")
btMul = GtkButton("*")
btDiv = GtkButton("/")
btAC = GtkButton("AC")
btDot = GtkButton(".")
btEqual = GtkButton("=")
btOFF = GtkButton("OFF")

# Función para boton OFF
signal_connect(btOFF, :clicked) do widget
    destroy(cal)
end

# Función para boton 1
signal_connect(bt1, :clicked) do widget
    if status == 0
        set_gtk_property!(pantalla, :label, "1")
        global status = 1
    else
        msg = get_gtk_property(pantalla, :label, String)
        newMsg = string(msg,"1")
        set_gtk_property!(pantalla, :label, newMsg)
        global status = 1
    end
end

signal_connect(cal, "key-press-event") do widget, event
    if event.keyval == 49 || event.keyval == 65457
        if status == 0
            set_gtk_property!(pantalla, :label, "1")
            global status = 1
        else
            msg = get_gtk_property(pantalla, :label, String)
            newMsg = string(msg,"1")
            set_gtk_property!(pantalla, :label, newMsg)
            global status = 1
        end
    end
end

# Función para boton 2
signal_connect(bt2, :clicked) do widget
    if status == 0
        set_gtk_property!(pantalla, :label, "2")
        global status = 1
    else
        msg = get_gtk_property(pantalla, :label, String)
        newMsg = string(msg,"2")
        set_gtk_property!(pantalla, :label, newMsg)
        global status = 1
    end
end

# Función para boton 3
signal_connect(bt3, :clicked) do widget
    if status == 0
        set_gtk_property!(pantalla, :label, "3")
        global status = 1
    else
        msg = get_gtk_property(pantalla, :label, String)
        newMsg = string(msg,"3")
        set_gtk_property!(pantalla, :label, newMsg)
        global status = 1
    end
end

# Función para boton 4
signal_connect(bt4, :clicked) do widget
    if status == 0
        set_gtk_property!(pantalla, :label, "4")
        global status = 1
    else
        msg = get_gtk_property(pantalla, :label, String)
        newMsg = string(msg,"4")
        set_gtk_property!(pantalla, :label, newMsg)
        global status = 1
    end
end

# Función para boton 5
signal_connect(bt5, :clicked) do widget
    if status == 0
        set_gtk_property!(pantalla, :label, "5")
        global status = 1
    else
        msg = get_gtk_property(pantalla, :label, String)
        newMsg = string(msg,"5")
        set_gtk_property!(pantalla, :label, newMsg)
        global status = 1
    end
end

# Función para boton 6
signal_connect(bt6, :clicked) do widget
    if status == 0
        set_gtk_property!(pantalla, :label, "6")
        global status = 1
    else
        msg = get_gtk_property(pantalla, :label, String)
        newMsg = string(msg,"6")
        set_gtk_property!(pantalla, :label, newMsg)
        global status = 1
    end
end

# Función para boton 7
signal_connect(bt7, :clicked) do widget
    if status == 0
        set_gtk_property!(pantalla, :label, "7")
        global status = 1
    else
        msg = get_gtk_property(pantalla, :label, String)
        newMsg = string(msg,"7")
        set_gtk_property!(pantalla, :label, newMsg)
        global status = 1
    end
end

# Función para boton 8
signal_connect(bt8, :clicked) do widget
    if status == 0
        set_gtk_property!(pantalla, :label, "8")
        global status = 1
    else
        msg = get_gtk_property(pantalla, :label, String)
        newMsg = string(msg,"8")
        set_gtk_property!(pantalla, :label, newMsg)
        global status = 1
    end
end

# Función para boton 9
signal_connect(bt9, :clicked) do widget
    if status == 0
        set_gtk_property!(pantalla, :label, "9")
        global status = 1
    else
        msg = get_gtk_property(pantalla, :label, String)
        newMsg = string(msg,"9")
        set_gtk_property!(pantalla, :label, newMsg)
        global status = 1
    end
end

# Función para boton 0
signal_connect(bt0, :clicked) do widget
    msg = get_gtk_property(pantalla, :label, String)
    newMsg = string(msg,"0")
    set_gtk_property!(pantalla, :label, newMsg)
    global status = 1
end

# Función para botonAC
signal_connect(btAC, :clicked) do widget
    set_gtk_property!(pantalla, :label, 0)
    global status = 0
end

signal_connect(btPlus, :clicked) do widget
    msg = get_gtk_property(pantalla, :label, String)
    v1S = msg
    global v1N = parse(Float64, v1S)
    newMsg = string(msg,"+")
    global numVal = length(newMsg)
    set_gtk_property!(pantalla, :label, newMsg)
end

signal_connect(btEqual, :clicked) do widget
    global v1N, numVal
    msg = get_gtk_property(pantalla, :label, String)
    global v2S = msg[(numVal+1):end]
    global v2N = parse(Float64, v2S)
    SUM = v1N + v2N
    set_gtk_property!(pantalla, :label, SUM)
end
g1 = GtkGrid()
set_gtk_property!(g1,:column_homogeneous,true)
set_gtk_property!(g1,:row_homogeneous,true)
set_gtk_property!(g1, :column_spacing, 15)
set_gtk_property!(g1, :row_spacing, 15)
set_gtk_property!(g1, :margin_top, 20)
set_gtk_property!(g1, :margin_left, 20)
set_gtk_property!(g1, :margin_right, 20)
set_gtk_property!(g1, :margin_bottom, 20)

borde = GtkFrame()
pantalla = GtkLabel("0")
set_gtk_property!(pantalla, :tooltip_text, "Resultado")
set_gtk_property!(pantalla, :xalign, 1.0)

g1[1:4,1] = borde
g1[1:2,2] = btAC
g1[3:4,2] = btOFF
g1[3,3] = bt9
g1[1,3] = bt7
g1[2,3] = bt8
g1[3,3] = bt9
g1[1,4] = bt4
g1[2,4] = bt5
g1[3,4] = bt6
g1[1,5] = bt1
g1[2,5] = bt2
g1[3,5] = bt3
g1[1,6] = bt0
g1[2,6] = btDot
g1[3,6] = btEqual

g1[4,3] = btMinus
g1[4,4] = btPlus
g1[4,5] = btMul
g1[4,6] = btDiv

push!(borde, pantalla)
push!(cal, g1)
showall(cal)


text:
module TestUpdating

using Test
using BenchmarkCI: GitUtils

function setup_dummy_user()
    run(`git config user.email "DUMMY@users.noreply.github.com"`)
    run(`git config user.name DUMMY`)
end

function init_random_repo(dir, branch)
    mkpath(dir)
    cd(dir) do
        run(`git init .`)
        setup_dummy_user()

        run(`git checkout -b $branch`)
        write("README.txt", "hello")
        run(`git add .`)
        run(`git commit --message "First commit"`)
        run(`git checkout -b $branch.0`)
    end
end

@testset "updating" begin
    mktempdir() do dir
        branch = "somebranch"
        origin = joinpath(dir, "origin")
        init_random_repo(origin, branch)
        url = "file://" * abspath(origin)

        GitUtils.updating(url, branch) do ctx
            write("spam.txt", "hello")
        end
        cd(origin) do
            run(`git checkout $branch`)
            @test read("spam.txt", String) == "hello"
            run(`git checkout -`)
        end

        GitUtils.updating(url, branch) do ctx
            write("spam.txt", "hello hello")
        end
        cd(origin) do
            run(`git checkout $branch`)
            @test read("spam.txt", String) == "hello hello"
            run(`git checkout -`)
        end
    end
end

@testset "prepare_ssh_command" begin
    mktempdir() do dir
        @test occursin(dir, GitUtils.prepare_ssh_command(dir, "dummy key"))
    end
end

@testset "_push_with_retry" begin
    mktempdir() do dir
        branch = "somebranch"
        origin = joinpath(dir, "origin")
        init_random_repo(origin, branch)
        url = "file://" * abspath(origin)

        workdir = joinpath(dir, "workdir")
        run(`git clone $url $workdir`)
        cd(workdir) do
            setup_dummy_user()
        end

        # Advance `origin`
        cd(origin) do
            run(`git checkout $branch`)
            write("file-1", "hello")
            run(`git add .`)
            run(`git commit --message "Add file-1"`)
            run(`git checkout -`)
        end

        cd(workdir) do
            # Advance local `workdir`:
            write("file-2", "hello")
            run(`git add .`)
            run(`git commit --message "Add file-2"`)

            logs, = Test.collect_test_logs() do
                GitUtils._push_with_retry(args -> `git $args`, branch)
            end
            all_messages = join((l.message for l in logs), "\n")
            @test occursin("retrying", all_messages)
        end
    end
end

end  # module


text:
#=
  False coin in Turing.jl

  From Peter Winkler
  Celebration of Mind (Gathering For Gardner):
  "CoM Apr 2021 - Drawing from Urns - Peter Winkler"
  https://www.youtube.com/watch?v=A0GxvYPTKDk&list=PL5_D5QQxvWKVtX4Vklvhsf5KTYZBLDSwO&index=3

  @3:43
  """
  You have two quarters.
  Both are sitting on a table heads up.
  One of them is an ordinary quarter, but the
  other has a head on the other side as well.

  You get two flips to help you guess which is which.

  Should you flip each coin once or one coin twice.
  """

  According to this model it's better to flip the same coin twice
  than once for each coin. But that's under the condition that
  all the flips shows head. If either flip shows tail then it's
  settled that it's the ordinary coin.

  """
  Distributions of variable coinA (num:0)
  1.00000 =>    8066  (0.806600)
  2.00000 =>    1934  (0.193400)

  Distributions of variable sideA (num:0)
  1.00000 =>   10000  (1.000000)

  Distributions of variable coinB (num:0)
  2.00000 =>    8066  (0.806600)
  1.00000 =>    1934  (0.193400)

  Distributions of variable sideB (num:0)
  1.00000 =>   10000  (1.000000)

  Distributions of variable coinSameIsFake (num:0)
  1.00000 =>    8066  (0.806600)
  0.00000 =>    1934  (0.193400)

  Distributions of variable coinBothIsFake (num:0)
  0.00000 =>    8066  (0.806600)
  1.00000 =>    1934  (0.193400)

  """

  (See false_coin.jl for an easier variant of the problem.)

=#

using Turing
include("jl_utils.jl")

@model function false_coin2()
    #  The two coins
    coin1 = 1 #  fake coin: head,head
    coin2 = 2 #  proper coin: head,tail
    coins = [coin1,coin2];

    #  Sides of the coins
    head = 1
    tail = 2

    function flip_coin(cond)
        return cond ? Categorical([1.0,0.0]) :  Categorical([0.5,0.5])
    end

    #  The setup: both coins are head (no real flipping)
    coinA ~ Categorical([0.5,0.5])
    sideA ~ flip_coin(coinA == coin1)
    true ~ Dirac(sideA == head)

    coinB ~ Dirac(coinA == coin1 ? coin2 : coin1)
    sideB ~ flip_coin(coinB == coin1);
    true ~ Dirac(sideB == head)

    #  condition(coinA == coin1 || coinB == coin1)

    #  Flip same coin
    coinSame_1 ~ Dirac(coinA)
    coinSame_1_flip ~ flip_coin(coinSame_1 == coin1);
    coinSame_2 ~ Dirac(coinSame_1)
    coinSame_2_flip ~ flip_coin(coinSame_2 == coin1);
    true ~ Dirac(coinSame_1_flip == head)
    true ~ Dirac(coinSame_2_flip == head)

    #  Flip both coins
    coinBoth_1 ~ Dirac(coinA)
    coinBoth_1_flip ~ flip_coin(coinBoth_1 == coin1);
    coinBoth_2 ~ Dirac(coinB)
    coinBoth_2_flip ~ flip_coin(coinBoth_2 == coin1);
    true ~ Dirac(coinBoth_1_flip == head)
    true ~ Dirac(coinBoth_2_flip == head)

    coinSameIsFake ~ Dirac(coinA ==coin1)
    coinBothIsFake ~ Dirac(coinB ==coin1)
end

model = false_coin2()
# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, IS(), 10_000)
# chns = sample(model, SMC(), 10_000)

display(chns)
show_var_dist_pct(chns,:coinA)
show_var_dist_pct(chns,:sideA)
show_var_dist_pct(chns,:coinB)
show_var_dist_pct(chns,:sideB)

show_var_dist_pct(chns,:coinSameIsFake)
show_var_dist_pct(chns,:coinBothIsFake)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
import Term.Renderables: Renderable, RenderableText, AbstractRenderable
import Term.Segments: Segment
import Term: fillin

@testset "\e[34mSegment" begin
    seg = Segment("test", "default")
    @test seg.text == "\e[22mtest\e[22m"
    @test seg.measure.w == 4

    @test_nothrow println(seg)

    seg = Segment("test", "red")
    @test seg.text == "\e[31mtest\e[39m"
    @test seg.measure.w == 4

    seg = Segment("aa\n{blue}123{/blue}")
    @test size(seg.measure) == (3, 2)

    seg = Segment("test")
    @test seg * "test" isa Segment
    @test (seg * "t2").text == "testt2"
    @test ("t2" * seg).text == "t2test"
    @test (seg * seg).text == "testtest"
end

@testset "\e[34mRenderables - Renderable" begin
    lorem = "Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore"
    r = RenderableText(lorem)

    @test size(r.measure) == (length(lorem), 1)

    r2 = RenderableText(r)
    @test size(r.measure) == (length(lorem), 1)

    width = 22
    r = RenderableText(lorem; width = width)
    @test string(r) ==
          "Lorem ipsum dolor sit \namet, consectetur     \nadipiscing elit, sed  \ndo eiusmod tempor     \nincididunt ut labore  "
    @test r.measure.w == width

    r = RenderableText(lorem; width = width, style = "red")
    @test string(r) ==
          "\e[31mLorem ipsum dolor sit \e[39m\n\e[31mamet, consectetur     \e[39m\n\e[31madipiscing elit, sed  \e[39m\n\e[31mdo eiusmod tempor     \e[39m\n\e[31mincididunt ut labore  \e[39m"
    @test r.measure.w == width

    @test string(RenderableText("a string")) == "a string"
    @test string(RenderableText("a\nstring")) == "a     \nstring"
end


text:
using Pkg

# This is Pkg.precompile without the try/catch
function precompile_fatal()
    ctx = Pkg.API.Context()
    pkgids = [Base.PkgId(uuid, name) for (name, uuid) in ctx.env.project.deps if !Pkg.API.is_stdlib(uuid)]
    if ctx.env.pkg !== nothing && isfile( joinpath( dirname(ctx.env.project_file), "src", ctx.env.pkg.name * ".jl") )
        push!(pkgids, Base.PkgId(ctx.env.pkg.uuid, ctx.env.pkg.name))
    end

    for pkg in pkgids
        paths = Base.find_all_in_cache_path(pkg)
        sourcepath = Base.locate_package(pkg)
        sourcepath === nothing && continue
        # Heuristic for when precompilation is disabled
        occursin(r"\b__precompile__\(\s*false\s*\)", read(sourcepath, String)) && continue
        stale = true
        for path_to_try in paths::Vector{String}
            staledeps = Base.stale_cachefile(sourcepath, path_to_try)
            staledeps === true && continue
            # TODO: else, this returns a list of packages that may be loaded to make this valid (the topological list)
            stale = false
            break
        end
        if stale
            Base.compilecache(pkg, sourcepath)
        end
    end
    return nothing
end

try
    Pkg.instantiate(;verbose=true)
    precompile_fatal()
catch
    try
        @info "Precompilation faild, trying Pkg.build()"
        Pkg.build(;verbose=true)
        precompile_fatal()
    catch
        @info "Pkg.build() failed, starting from clean depot"
        rm(DEPOT_PATH[1]; recursive=true)
        Pkg.instantiate(;verbose=true)
        precompile_fatal()
    end
end

Pkg.status()

max_stars_repo_name: CliMA/slurm-buildkite

max_issues_repo_name: CliMA/slurm-buildkite

max_forks_repo_name: CliMA/slurm-buildkite


text:
using Flux
using Zygote
using MAT
#using Plots
using Statistics
using Distributions
using LinearAlgebra
using Random
using BSON
using JLD
using PyPlot
using CUDA
#using Base.Iterators: partition

gpu_id = 1
CUDA.device!(gpu_id)

#=
mf = matopen("/home/naxos2-raid1/maxiao/Projects/pTXDL/JuliaPrograms/ImageBasedDLForpTX.jl/Data/Data_processed_71Vol.mat")
Xtrain = read(mf, "Xtrain")
Ytrain = read(mf, "Ytrain")
close(mf)

for idx_subj = 1:size(Xtrain,5)
    Xtrain[:,:,:,:,idx_subj] = Xtrain[:,:,:,:,idx_subj]./maximum(Xtrain[:,:,:,:,idx_subj])
    Ytrain[:,:,:,:,idx_subj] = Ytrain[:,:,:,:,idx_subj]./maximum(Ytrain[:,:,:,:,idx_subj])
end

Xall = Xtrain
Yall = Ytrain

Xtrain = Xall[:,:,:,:,1:4]
Ytrain = Yall[:,:,:,:,1:4]

Sall = size(Xtrain)
Xtrain = reshape(Xtrain,Sall[1],Sall[2],1,Sall[3]*Sall[4]*Sall[5])

Sall = size(Ytrain)
Ytrain = reshape(Ytrain,Sall[1],Sall[2],1,Sall[3]*Sall[4]*Sall[5])

XValid = Xall[:,:,:,:,5]
YValid = Yall[:,:,:,:,5]

Xall = nothing;
Yall = nothing;

mf = matopen("/home/naxos2-raid1/maxiao/Projects/pTXDL/JuliaPrograms/ImageBasedDLForpTX.jl/Data/T1wBrain_processed.mat")
T1wBrain_resize_ds = read(mf, "T1wBrain_resize_ds")
close(mf)

for idx_subj = 1:size(T1wBrain_resize_ds,4)
    T1wBrain_resize_ds[:,:,:,idx_subj] = T1wBrain_resize_ds[:,:,:,idx_subj]./maximum(T1wBrain_resize_ds[:,:,:,idx_subj])
end

Xtrain_new = zeros(Float32,size(Xtrain,1),size(Xtrain,1),2,size(Xtrain,4))
Xtrain_new[:,:,1,:] = Xtrain
Xtrain_new[:,:,2,:] = reshape(repeat(T1wBrain_resize_ds[:,:,:,1:4],1,1,71,1),size(Xtrain,1),size(Xtrain,1),1,size(Xtrain,4))

Xtrain = Xtrain_new
Xtrain_new=nothing

Sall = size(XValid)
XValid = reshape(XValid,Sall[1],Sall[2],1,Sall[3]*Sall[4])
YValid = reshape(YValid,Sall[1],Sall[2],1,Sall[3]*Sall[4])

XValid_new = zeros(Float32,size(XValid,1),size(XValid,1),2,size(XValid,4))
XValid_new[:,:,1,:] = XValid
XValid_new[:,:,2,:] = reshape(repeat(T1wBrain_resize_ds[:,:,:,5],1,1,71),size(XValid,1),size(XValid,1),1,size(XValid,4))

XValid = XValid_new
XValid_new=nothing

JLD.save("Data/XYtrain_XYvalid_processed_71Vol_T1w.jld","Xtrain",Xtrain,"Ytrain",Ytrain,"XValid",XValid,"YValid",YValid)
=#

Xtrain = JLD.load("Data/XYtrain_XYvalid_processed_71Vol_T1w.jld","Xtrain")
Ytrain = JLD.load("Data/XYtrain_XYvalid_processed_71Vol_T1w.jld","Ytrain")
XValid = JLD.load("Data/XYtrain_XYvalid_processed_71Vol_T1w.jld","XValid")
YValid = JLD.load("Data/XYtrain_XYvalid_processed_71Vol_T1w.jld","YValid")

idx_shuffle = shuffle(collect(1:size(Xtrain,4)))
Xtrain = Xtrain[:,:,1,idx_shuffle]
Ytrain = Ytrain[:,:,:,idx_shuffle]

Xtrain = reshape(Xtrain,size(Xtrain,1),size(Xtrain,1),1,size(Xtrain,3))

Xtrain = gpu(Xtrain)
Ytrain = gpu(Ytrain)

batch_size = 20

train_loader = Flux.Data.DataLoader(Xtrain,Ytrain, batchsize=batch_size,shuffle=true)

XValid = XValid[:,:,1,:]
XValid = reshape(XValid,size(XValid,1),size(XValid,1),1,size(XValid,3))

XValid = gpu(XValid)
YValid = gpu(YValid)

channels_in=1
channels_out=1

# Unet define
function _random_normal(shape...)
  return Float32.(rand(Normal(0.f0,0.02f0),shape...))
end

BatchNormWrap(out_ch) = BatchNorm(out_ch)

ResNetConvBlock(chs, kernel = (3, 3)) =
    Chain(Conv(kernel, chs=>chs,pad = (1, 1);init=_random_normal),
          BatchNorm(chs),
          x->relu.(x),
          Conv(kernel, chs=>chs,pad = (1, 1);init=_random_normal),
          BatchNorm(chs))

r = Chain(
            Conv((4,4),channels_in=>64,pad=(1,1),stride=(2,2);init=_random_normal),
            x->relu.(x),
            SkipConnection(ResNetConvBlock(64),(mx,x)->mx+x),
            SkipConnection(ResNetConvBlock(64),(mx,x)->mx+x),
            SkipConnection(ResNetConvBlock(64),(mx,x)->mx+x),
            SkipConnection(ResNetConvBlock(64),(mx,x)->mx+x),
            ConvTranspose((2,2),64=>channels_out,stride=(2,2);init=_random_normal),
            BatchNorm(channels_out)
        )

r = gpu(r)

#loss(x,y) = mean(abs,(u(x).-y))
loss(x,y) = mean((r(x).-y).^2)
validationloss() = mean(Flux.mse(r(XValid[:,:,:,(idx-1)*100+1:idx*100]),YValid[:,:,:,(idx-1)*100+1:idx*100]) for idx=1:round(Int,size(XValid,4)/100))
#validationloss() = mean(mean(abs,(u(XValid[:,:,:,(idx-1)*100+1:idx*100]).-YValid[:,:,:,(idx-1)*100+1:idx*100])) for idx=1:round(Int,size(XValid,4)/100))

opt = ADAM(0.001)

function my_custom_train!(loss, ps, data, opt)
    # training_loss is declared local so it will be available for logging outside the gradient calculation.
    local training_loss
    ps = Params(ps)
    for d in data
      if d isa AbstractArray{<:Number}
        gs = Flux.gradient(ps) do
          training_loss = loss(d)
          return training_loss
        end
      else
        gs = Flux.gradient(ps) do
          training_loss = loss(d...)
          return training_loss
        end
      end
      # Insert whatever code you want here that needs training_loss, e.g. logging.
      # logging_callback(training_loss)
      # Insert what ever code you want here that needs gradient.
      # E.g. logging with TensorBoardLogger.jl as histogram so you can see if it is becoming huge.
      #println("training_loss = ",training_loss,"; gradient[1] = ",gs[ps[1]],"; W1 = ",ps.order.data[1][1],"; b1 = ",ps.order.data[2][1])
      println("training_loss = ",training_loss)
      Flux.update!(opt, ps, gs)
      # Here you might like to check validation set accuracy, and break out to do early stopping.
    end
  end

loss_valid = zeros(100)
for epoch in 1:100
  learn_rate = 0.001*exp(log(0.1)/99*(epoch-1))
  opt = ADAM(learn_rate)
  @time Flux.train!(loss, Flux.params(r), train_loader, opt)
  #@time my_custom_train!(loss, Flux.params(r), train_loader, opt)
  #@time Flux.train!(loss, Flux.params(u), train_set, opt)
  #model = cpu(u)
  #BSON.@save string("Results/model_resnet_ptx_processed_7dir_gpu_epoch",string(epoch),".bson") model
  #u = gpu(model)
  loss_valid[epoch] = validationloss()
  loss_valid[epoch]
end

model = cpu(r)
BSON.@save string("Results/model_resnet_ptx_processed_71Vol_gpu_batch20_epoch100.bson") model

#test

BSON.@load string("Results/model_resnet_ptx_processed_71Vol_gpu_batch20_epoch100.bson") model
slc_list = collect(1:100)
dir_list = collect(1:7)
XtestOutput4write = zeros(size(XValid,1),size(XValid,2),length(slc_list),length(dir_list))

Xtest = zeros(Float32,size(XValid,1),size(XValid,2),1,1)
Ytest = zeros(Float32,size(XValid,1),size(XValid,2),1,1)
for idx_dir_test in 1:length(dir_list)
  for idx_slc_test in 1:length(slc_list)
    idx_test = idx_slc_test + (idx_dir_test-1)*100
    Xtest[:,:,:,1] = XValid[:,:,:,idx_test]
    #Ytest[:,:,1,1] = YValid[:,:,idx_slc_test,idx_dir_test]
    XtestOutput4write[:,:,idx_slc_test,idx_dir_test] = model(Xtest)
  end
end

file = matopen("Results/XValidOutput_subj5_processed_71dir_gpu_resnet_batch20_epoch100.mat", "w")
write(file, "XtestOutput4write", XtestOutput4write)
#write(file, "XtestOutput4write_slc10_25_dir1_2", XtestOutput4write_slc10_25_dir1_2)
close(file)

# testing on training data
BSON.@load string("model_unet_ptx_processed_7dir_epoch",string(10),".bson") u
XValid = Xtrain[:,:,:,1:700]
XValid = reshape(XValid,size(XValid,1),size(XValid,2),100,7)
XtestOutput4write = zeros(size(XValid))

Xtest = zeros(Float32,size(XValid,1),size(XValid,2),1,1)
Ytest = zeros(Float32,size(XValid,1),size(XValid,2),1,1)
for idx_dir_test in 1:size(XValid,4)
  for idx_slc_test in 1:size(XValid,3)
    Xtest[:,:,1,1] = XValid[:,:,idx_slc_test,idx_dir_test]
    #Ytest[:,:,1,1] = YValid[:,:,idx_slc_test,idx_dir_test]
    XtestOutput4write[:,:,idx_slc_test,idx_dir_test] = u(Xtest)
  end
end

file = matopen("XValidOutput_traingsubj1_processed_7dir_epoch10.mat", "w")
write(file, "XtestOutput4write", XtestOutput4write)
close(file)


text:
@testset "Routes.jl" begin
    filename = "data\\test.routes"
    dc = DDR2import.Routes.read(filename)
    @test dc["ABESI8TLIME"].type == "DP"
    @test dc["ABDIL1ALFMD"].route[2].wp == "TUPOX"
    @test dc["ABDIL1ALFTZ"].route[3].location_type == "SP"
    @test dc["ABIRO2DGMTT"].route[2].wp == "ABIRO"
end


text:
module Test

using Test
using MathOptInterface
using ConstraintProgrammingExtensions
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const CP = ConstraintProgrammingExtensions

const Config{T <: Real} = MOIT.Config
const setup_test = MOIT.setup_test

"""
[See MOI docs](https://jump.dev/MathOptInterface.jl/stable/submodules/Test/reference/#MathOptInterface.Test.runtests).
"""
function runtests(
    model::MOI.ModelLike,
    config::Config;
    include::Vector{String} = String[],
    exclude::Vector{String} = String[],
    warn_unsupported::Bool = false,
)
    # Copy-paste from MOI.
    for name_sym in names(@__MODULE__; all = true)
        name = string(name_sym)
        if !startswith(name, "test_")
            continue  # All test functions start with test_
        elseif !isempty(include) && !any(s -> occursin(s, name), include)
            continue
        elseif !isempty(exclude) && any(s -> occursin(s, name), exclude)
            continue
        end
        @testset "$(name)" begin
            test_function = getfield(@__MODULE__, name_sym)
            c = copy(config)
            tear_down = MOIT.setup_test(test_function, model, c)
            # Make sure to empty the model before every test!
            MOI.empty!(model)
            try
                test_function(model, c)
            catch err
                MOIT._error_handler(err, name, warn_unsupported)
            end
            if tear_down !== nothing
                tear_down()
            end
        end
    end
    return
end

# Include all test files in CP.

for file in readdir(@__DIR__)
    if startswith(file, "test_")
        include(file)
    end
end

end


text:
function countdistinct_vectorofvariables_test(
    model::MOI.ModelLike,
    config::MOIT.Config,
)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.EqualTo{Int})
    @test MOI.supports_constraint(
        model,
        MOI.VectorOfVariables,
        CP.CountDistinct,
    )

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x3, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x4, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, x1, MOI.EqualTo(1))
    c2 = MOI.add_constraint(model, x2, MOI.EqualTo(1))
    c3 = MOI.add_constraint(model, x3, MOI.EqualTo(2))

    c4 = MOI.add_constraint(model, MOI.VectorOfVariables([x4, x1, x2, x3]), CP.CountDistinct(3))

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, x3)
    @test MOI.is_valid(model, x4)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)
    @test MOI.is_valid(model, c3)
    @test MOI.is_valid(model, c4)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x2) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x3) == 2
        @test MOI.get(model, MOI.VariablePrimal(), x4) == 2
    end
end

function countdistinct_vectoraffinefunction_test(
    model::MOI.ModelLike,
    config::MOIT.Config,
)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{Int},
        MOI.EqualTo{Int},
    )
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{Int},
        CP.CountDistinct,
    )

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x3, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x4, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, 1 * x1, MOI.EqualTo(1))
    c2 = MOI.add_constraint(model, 1 * x2, MOI.EqualTo(1))
    c3 = MOI.add_constraint(model, 1 * x3, MOI.EqualTo(2))

    c4 = MOI.add_constraint(model, MOIU.vectorize([x4, x1, x2, x3]), CP.CountDistinct(3))

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, x3)
    @test MOI.is_valid(model, x4)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)
    @test MOI.is_valid(model, c3)
    @test MOI.is_valid(model, c4)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x2) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x3) == 2
        @test MOI.get(model, MOI.VariablePrimal(), x4) == 2
    end
end

const countdistincttests = Dict(
    "countdistinct_vectorofvariables" =>
        countdistinct_vectorofvariables_test,
    "countdistinct_vectoraffinefunction" =>
        countdistinct_vectoraffinefunction_test,
)

MOIT.@moitestset countdistinct


text:
@testset "API communication" begin
    @test 1 == 1
end


text:
using Plots

function plot_decision_boundary(elm::T₁,
                                samples::T₂,
                                targets::T₃) where {T₁ <: ELM,
                                                    T₂ <: AbstractMatrix,
                                                    T₃ <: AbstractVector}
    x, y = samples[1,:], samples[2,:]
    colours = [label == 1 ? "#F8766D" : "#00BFC4" for label in targets]
    rng = -10:0.1:10
    plt = scatter(x, y,
            markersize = 4,
            markerstrokewidth = 0,
            markercolor = colours,
            xaxis = false,
            yaxis = false,
            grid = false,
            legend = false)
    plt = contour!(rng, rng, (x, y) -> predict(elm, [x, y]),
            levels = [0],
            linewidth = 1,
            linecolor = :black,
            legend = false)
end


text:
# Tests SerializedMsgFormat in async mode operation.
# Both server and client run in the same process.
# In this mode the transport can be in-memory.
using JuliaWebAPIPlugins

include("srvr.jl")
include("clnt.jl")

function test_amqp_transport_ser()
    stport = JuliaWebAPIPlugins.AMQPTransport("juliawebapi", :server)
    run_srvr(JuliaWebAPI.SerializedMsgFormat(), stport, true)
    run_clnt(JuliaWebAPI.SerializedMsgFormat(), JuliaWebAPIPlugins.AMQPTransport("juliawebapi", :client))
    close(stport)
end

function test_amqp_transport_json()
    stport = JuliaWebAPIPlugins.AMQPTransport("juliawebapi", :server)
    run_srvr(JuliaWebAPI.JSONMsgFormat(), stport, true)
    run_clnt(JuliaWebAPI.JSONMsgFormat(), JuliaWebAPIPlugins.AMQPTransport("juliawebapi", :client))
    close(stport)
end

# run tests
!isempty(ARGS) && (ARGS[1] == "--ser") && test_amqp_transport_ser()
!isempty(ARGS) && (ARGS[1] == "--json") && test_amqp_transport_json()

max_stars_repo_name: tanmaykm/JuliaWebAPIPlugins

max_issues_repo_name: tanmaykm/JuliaWebAPIPlugins

max_forks_repo_name: tanmaykm/JuliaWebAPIPlugins


text:
module Murmur3

export hash32

@inline function rotl32(var::Integer, value::Integer)
    (var << value) | ((var & 0xffffffff) >> (32-value))
end

@inline function fmix32(value::Integer)
    value ⊻= ((value & 0xffffffff) >> 16)
    value *= 0x85ebca6b
    value ⊻= ((value & 0xffffffff) >> 13)
    value *= 0xc2b2ae35
    value ⊻= ((value & 0xffffffff) >> 16)
end

function hash32(data::Vector{UInt8}, seed::UInt32=zero(UInt32))
    c1 = 0xcc9e2d51
    c2 = 0x1b873593
    h = seed

    len = length(data)
    blocks = div(len, 4)
    p = convert(Ptr{UInt32}, pointer(data))

    # Body
    for next_block ∈ 1:blocks
        k = big(unsafe_load(p, next_block))
        k *= c1
        k = rotl32(k, 15)
        k *= c2

        h ⊻= k
        h = rotl32(h, 13)
        h = h * 5 + 0xe6546b64
    end

    # Tail
    k = zero(BigInt)
    remainder = len & 3
    last_ = blocks * 4
    if remainder == 3
        k |= (UInt32(data[last_ + 3]) & 0xff) << 16
    end
    if remainder ∈ 2:3
        k |= (UInt32(data[last_ + 2]) & 0xff) << 8
    end
    if remainder ∈ 1:3
        k |= UInt32(data[last_ + 1]) & 0xff
        k *= c1
        k = rotl32(k, 15)
        k *= c2
        h ⊻= k
    end

    # Finalization
    h ⊻= len
    h = fmix32(h)

    h & 0xffffffff

end

hash32(data::Vector{UInt8}, seed::Int) = hash32(data, UInt32(seed))
hash32(data::AbstractString, seed::UInt32) = hash32(UInt8.(collect(data)), seed)
hash32(data::AbstractString, seed::Int) = hash32(UInt8.(collect(data)), UInt32(seed))
hash32(data::AbstractString) = hash32(UInt8.(collect(data)), UInt32(0))

end


text:
struct PulseSummary
  pretrig_mean      ::Float64
  pretrig_rms       ::Float64
  pulse_average     ::Float64
  pulse_rms         ::Float64
  rise_time         ::Float64
  postpeak_deriv    ::Float64
  peak_index        ::Int16
  peak_value        ::UInt16
  min_value         ::UInt16
end

function summary(data::LJH.LJHRecord, npresamples, nsamples, average_pulse_peak_index, frametime)
  summarize(data.data, npresamples, nsamples, average_pulse_peak_index, frametime)
end

function summarize(data::Vector, npresamples, nsamples, average_pulse_peak_index, frametime)
  length(data) == nsamples || error("wrong length data. nsamples: $nsamples\n data: $data")
  s = s2 = 0
  min_idx = 0
  peak_idx = 0
  peak_val = 0
  min_idx = 0
  min_val = typemax(Int)
  npostsamples = length(data)-npresamples
  for j = 1:npresamples
      d=Int(data[j])
      if d > peak_val
          peak_idx, peak_val = j, d
      elseif d < min_val
          min_idx, min_val = j,d
      end
      s+=d
      s2+=d*d
  end
  ptm = s/npresamples # pre trigger mean
  pretrig_rms = sqrt(abs(s2/npresamples - ptm*ptm))

  # Now post-trigger calculations
  s = s2 = 0
  for j = npresamples+1:length(data)
      d=Int(data[j])
      if d > peak_val
          peak_idx, peak_val = j, d
      elseif d < min_val
          min_idx, min_val = j,d
      end
      # d2=d-ptm
      s+=d
      s2+=d*d
  end

  rise_time::Float64 = estimate_rise_time(data, npresamples+1:peak_idx,
                                 peak_val, ptm, frametime)

  postpeak_deriv = max_timeseries_deriv_mass(data, average_pulse_peak_index)
  # @show length(data)
  # @show average_pulse_peak_index, npresamples, nsamples
  # @show deriv_range = average_pulse_peak_index:length(data)
  # deriv_work = zeros(length(deriv_range))
  # deriv_data = data[deriv_range]
  # @show typeof(deriv_data), typeof(deriv_work)
  # postpeak_deriv = max_timeseries_deriv!(deriv_work, deriv_data, true)

  if 0<=peak_val-ptm<=typemax(UInt16)
    peak_val_ret = round(Int,peak_val-ptm)
  else
    peak_val_ret=0
  end

  # Copy results into the PulseSummaries object
  pulse_average = s/npostsamples-ptm
  pulse_rms = sqrt(abs(s2/npostsamples - ptm*(ptm+2*pulse_average)))
  PulseSummary(ptm, pretrig_rms, pulse_average, pulse_rms, rise_time,
   postpeak_deriv, peak_idx, peak_val_ret, min_val)
end

"estimate_rise_time(pulserecord, searchrange, peakval, ptm, frametime)
uses the slope between the 10 and 90 percent points to extrapolate rise time to 100 percent
`pulserecord` = Vector
`searchrange` = range to look for risetime in, should go from base of pulse to peak, last(searchrange) is peakindex
`peakval` = maximum(pulserecord[postrig_region], has not had ptm subtracted off)
`ptm` = pretrigger mean
`frametime` = time spacing between points"
function estimate_rise_time(pulserecord, searchrange, peakval,ptm,frametime)
    idx10 = first(searchrange)
    peakindex = last(searchrange)
    (peakindex > length(pulserecord) || length(searchrange)==0) && (return Float64(length(pulserecord)))

    idx90 = peakindex
    thresh10 = 0.1*(peakval-ptm)+ptm
    thresh90 = 0.9*(peakval-ptm)+ptm
    for j = idx10:peakindex
        if pulserecord[j] > thresh10
            idx10 = j-1
            break
        end
    end
    for j = idx10+1:peakindex
        if pulserecord[j] > thresh90
            idx90 = j-1
            break
        end
    end
    #divide out fraction rise to extrapolate to 10% rise time
    fracrise = (pulserecord[idx90]-pulserecord[idx10])/(peakval-ptm)
    rise_time = (idx90-idx10)*frametime/fracrise
    return rise_time
end

function max_timeseries_deriv_mass(p, s)
  #following mass.analysis_algorithms.compute_max_deriv
  # use filter=SG and spike reject = true options from mass, not optional
  N = length(p)
  Nk = 5
  k1,k2,k3,k4,k5=2,1,0,-1,-2 # default kernel from mass
  t_max_deriv=0

  t0 = k5 * p[s+1] + k4 * p[s+2] + k3 * p[s+3] + k2 * p[s+4] + k1 * p[s+5]
  t1 = k5 * p[s+2] + k4 * p[s+3] + k3 * p[s+4] + k2 * p[s+5] + k1 * p[s+6]
  t2 = k5 * p[s+3] + k4 * p[s+4] + k3 * p[s+5] + k2 * p[s+6] + k1 * p[s+7]
  t_max_deriv = min(t2,t0)

  for j=s+8:N
    t3 = k5 * p[j-4] + k4*p[j-3] + k3*p[j-2] + k2*p[j-1] + k1*p[j]
          # t4 = t3 if t3 < t1 else t1
    t4 = t3<t1 ? t3 : t1
    t_max_deriv = max(t4,t_max_deriv)
    t0, t1, t2 = t1, t2, t3
  end
  t_max_deriv/10
end


text:
export project_cardinality_parallel

function project_cardinality_parallel(x::Vector,
                             k ::Int64,
                             n_proc::Int64
                             )
"""
Project the vector x onto the set of vectors with cardinality (l0 'norm') less then or equal to k.
project m onto {m | card(m)<= k} : x = argmin_x 1/2 ||x-m||2^2 s.t. card(x)<=k
A parallel version
"""

#sort_ind = sortperm( a, by=abs, rev=true)
#val_max=x[sort_ind[1:k]]
#x=zeros(length(x))
#x[sort_ind[1:k]]=val_max

ns = convert(Int64,floor((length(x)/n_proc)))
sort_ind=Vector{Vector{Int64}}(n_proc)
sort_ind[1] = sortperm( x[1:ns], by=abs, rev=true)
sort_ind[2] = sortperm( x[ns+1:end], by=abs, rev=true)

val_max=Vector{Float64}(2*k)
val_max[1:k]=x[sort_ind[1][1:k]]
val_max[1+k:k*2]=x[sort_ind[2][1:k]]

sort_final_ind = sortperm(val_max,by=abs,rev=true  )
#sort_ind[1]=sort_ind[1][1:k]
#sort_ind[2]=sort_ind[2][1:k]
#sort_ind_total=[sort_ind[1];sort_ind[2]]
x[sort_ind[1][k+1:end]]=0.0;
x[sort_ind[2][k+1:end]+ns]=0;.0

for i=(k+1):length(sort_final_ind)
  print(i)
  if sort_final_ind[i]<=k
    x[sort_ind[1][sort_final_ind[i]]]=0.0
  else
    x[sort_ind[2][sort_final_ind[i]-k]]=0.0
  end
end

  return x
end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
type OCCAOnDevice
end

type OCCAOnHost
end

ondev = OCCAOnDevice();
onhost = OCCAOnHost();

immutable OCCADevFlag
    flag::Uint8;
    function OCCADevFlag(ondev::OCCAOnDevice)
        return new(convert(1,Uint8));
    end
    function OCCADevFlag(onhost::OCCAOnHost)
        return new(convert(2,Uint8));
    end
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using Tricks

# %%
?static_hasmethod

# %%
struct A end
struct B end
abstract type FooBar end
struct Foo <: FooBar x::Int end
foomethod(foo::Foo) = nothing
struct Bar <: FooBar x::Int end

# %%
foo = Foo(23)
hasmethod(foomethod, Tuple{Foo})

# %%
attr_dynamic(fb::T) where T <: FooBar = hasmethod(foomethod, Tuple{T}) ? A() : B()
f_dynamic(fb::FooBar) = f(attr_dynamic(fb), fb)
f(::A, fb) = fb.x + 100
f(::B, fb) = fb.x - 100

# %%
f_dynamic(Foo(23))

# %%
f_dynamic(Bar(-23))

# %%
@code_typed f_dynamic(Foo(23))

# %%
foo = Foo(23)
static_hasmethod(foomethod, Tuple{Foo})

# %%
attr_static(fb::T) where T <: FooBar = static_hasmethod(foomethod, Tuple{T}) ? A() : B()
f_static(fb::FooBar) = f(attr_static(fb), fb)

# %%
f_static(Foo(23))

# %%
f_static(Bar(-23))

# %%
@code_typed f_static(Foo(23))

# %%
@code_typed f_static(Bar(-23))

# %%
foomethod(fb::Bar) = nothing

# %%
f_static(Bar(-23))

# %%
@code_typed f_static(Bar(-23))

# %%

max_stars_repo_path: 0005/static_hasmethod example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0005/static_hasmethod example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0005/static_hasmethod example.jl

max_forks_repo_name: genkuroki/public


text:
"Problem 10 in the Hock-Schittkowski suite"
function hs10()

  nlp = Model()

  @variable(nlp, x[i=1:2])
  setvalue(x[1], -10)
  setvalue(x[2],  10)

  @NLobjective(
    nlp,
    Min,
    x[1] - x[2]
  )

  @NLconstraint(
    nlp,
    -3 * x[1]^2 + 2 * x[1] * x[2] - x[2]^2 ≥ -1
  )

  return nlp
end


text:
using Base.Test
using DataFrames

let
	dv = dataones(3)
	push!(dv, 3.0)
	push!(dv, NA)

	@assert isequal(dv, DataVector[1.0, 1.0, 1.0, 3.0, NA])

	a, b = pop!(dv), pop!(dv)
	@assert isna(a)
	@assert b == 3.0

	enqueue!(dv, 3.0)
	enqueue!(dv, NA)

	@assert isequal(dv, DataVector[NA, 3.0, 1.0, 1.0, 1.0])

	a, b = shift(dv), shift(dv)
	@assert isna(a)
	@assert b == 3.0
end


text:
import FractionalTimeDG
using OffsetArrays
using PyPlot

include("../FPDE_problem.jl")
α = 1.0

Nx = 8
Nt = 4
rx = 4
rt = 3
Mt = 2rt

x = OffsetVector(range(0, L, length=Nx+1), 0:Nx)
grid = FractionalTimeDG.SpatialGrid(x, rx, ("essential", "essential"), rx)

U0 = FractionalTimeDG.spatial_L2_projection(u0, grid)
t, U = FractionalTimeDG.PDEDG(1.0, f, U0, grid, max_t, Nt, rt, Mt)

pts_per_time_interval = 5
pts_per_space_interval = 5
xvals, pcwise_t, pcwise_U = FractionalTimeDG.evaluate_pcwise_poly(
    U, t, grid, pts_per_time_interval, pts_per_space_interval)

figure(1)
Nxvals = length(xvals)
Ntvals = pts_per_time_interval
X = Float64[ xvals[i] for i = 1:Nxvals, j = 1:Ntvals ]
for n = 1:Nt
    T = Float64[ pcwise_t[j,n] for i = 1:Nxvals, j = 1:Ntvals ]
    Z = Float64[ pcwise_U[i,j,n] for i = 1:Nxvals, j = 1:Ntvals ]
    surf(X, T, Z)
end
xlabel(L"$x$")
ylabel(L"$t$")

t[1:Nt] = [ max_t * (n/Nt)^2 for n = 1:Nt ]
U = FractionalTimeDG.PDEDG(1.0, f, U0, grid, t, rt, Mt)
xvals, pcwise_t, pcwise_U = FractionalTimeDG.evaluate_pcwise_poly(
    U, t, grid, pts_per_time_interval, pts_per_space_interval)

figure(2)
Nxvals = length(xvals)
Ntvals = pts_per_time_interval
X = Float64[ xvals[i] for i = 1:Nxvals, j = 1:Ntvals ]
for n = 1:Nt
    T = Float64[ pcwise_t[j,n] for i = 1:Nxvals, j = 1:Ntvals ]
    Z = Float64[ pcwise_U[i,j,n] for i = 1:Nxvals, j = 1:Ntvals ]
    surf(X, T, Z)
end
xlabel(L"$x$")
ylabel(L"$t$")


text:
function update_psu_conf!(psu_conf, psu, refresh_cnt, base)
     # Meas Voltage more often
     refresh_cnt==base*1 && (psu_conf.C1.volt_meas = string(get_meas(psu, "C1", "voltage")))
     refresh_cnt==base*2 && (psu_conf.C2.volt_meas = string(get_meas(psu, "C2", "voltage")))
     refresh_cnt==base*3 && (psu_conf.C3.volt_meas = string(get_meas(psu, "C3", "voltage")))
     # # Meas Current more often
     refresh_cnt==base*4 && (psu_conf.C1.curr_meas = string(get_meas(psu, "C1", "current")))
     refresh_cnt==base*5 && (psu_conf.C2.curr_meas = string(get_meas(psu, "C2", "current")))
     refresh_cnt==base*6 && (psu_conf.C3.curr_meas = string(get_meas(psu, "C3", "current")))
     # # Update settings
     # # Voltage
     refresh_cnt==base*7 && (psu_conf.C1.volt_set = string(get_source_lev(psu, "C1")))
     refresh_cnt==base*8 && (psu_conf.C2.volt_set = string(get_source_lev(psu, "C2")))
     refresh_cnt==base*9 && (psu_conf.C3.volt_set = string(get_source_lev(psu, "C3")))
     # # Max Current
     refresh_cnt==base*10 && (psu_conf.C1.curr_set = string(get_max_curr(psu, "C1")))
     refresh_cnt==base*11 && (psu_conf.C2.curr_set = string(get_max_curr(psu, "C2")))
     refresh_cnt==base*12 && (psu_conf.C3.curr_set = string(get_max_curr(psu, "C3")))
     # # OVP
     refresh_cnt==base*13 && (psu_conf.C1.ovp_set = string(get_volt_protection(psu, "C1")))
     refresh_cnt==base*14 && (psu_conf.C2.ovp_set = string(get_volt_protection(psu, "C2")))
     refresh_cnt==base*15 && (psu_conf.C3.ovp_set = string(get_volt_protection(psu, "C3")))
     # # OCP
     refresh_cnt==base*16 && (psu_conf.C1.ocp_set = get_curr_protection(psu, "C1"))
     refresh_cnt==base*17 && (psu_conf.C2.ocp_set = get_curr_protection(psu, "C2"))
     refresh_cnt==base*18 && (psu_conf.C3.ocp_set = get_curr_protection(psu, "C3"))
	 # OUTPUT
	 (refresh_cnt==base*10 || refresh_cnt==base*19)&& (psu_conf.output_state = get_outp(psu, psu_conf.crt_chan))
	 # OUTPUT TOGGLE BUTTON STATE
	 if psu_conf.output_state == "off"
		psu_conf.output_btn.state = false
	 else
		psu_conf.output_btn.state = true
	 end
end

max_stars_repo_name: iuliancioarca/TIVM

max_issues_repo_name: iuliancioarca/TIVM

max_forks_repo_name: iuliancioarca/TIVM


text:
@shorthands meanplot
@shorthands autocorplot
@shorthands mixeddensity
@shorthands pooleddensity
@shorthands traceplot
@shorthands corner

struct _TracePlot; c; val; end
struct _MeanPlot; c; val;  end
struct _DensityPlot; c; val;  end
struct _HistogramPlot; c; val;  end
struct _AutocorPlot; lags; val;  end

# define alias functions for old syntax
const translationdict = Dict(
                        :traceplot => _TracePlot,
                        :meanplot => _MeanPlot,
                        :density => _DensityPlot,
                        :histogram => _HistogramPlot,
                        :autocorplot => _AutocorPlot,
                        :pooleddensity => _DensityPlot
                      )

const supportedplots = push!(collect(keys(translationdict)), :mixeddensity, :corner)

@recipe f(c::AbstractChains, s::Symbol) = c, [s]

@recipe function f(c::AbstractChains, i::Int;
    colordim = :chain,
    barbounds = (-Inf, Inf),
    maxlag = nothing,
    section = :parameters,
    append_chains = false)
    st = get(plotattributes, :seriestype, :traceplot)
    c = append_chains || st == :pooleddensity ? pool_chain(c) : c

    if colordim == :parameter
        title --> "Chain $(chains(c)[i])"
        label --> names(c)
        val = c.value[:, :, i]
    elseif colordim == :chain
        title --> names(c)[i]
        label --> map(k -> "Chain $(chains(c)[k])", 1:size(c)[3])
        val = c.value[:, i, :]
    else
        throw(ArgumentError("`colordim` must be one of `:chain` or `:parameter`"))
    end

    if st == :mixeddensity || st == :pooleddensity
        discrete = MCMCChains.indiscretesupport(c, barbounds)
        st = if colordim == :chain
            discrete[i] ? :histogram : :density
        else
            # NOTE: It might make sense to overlay histograms and density plots here.
            :density
        end
        seriestype := st
    end

    if st == :autocorplot
        lags = 0:(maxlag === nothing ? round(Int, 10 * log10(length(range(c)))) : maxlag)
        ac = MCMCChains.autocor(c, lags=collect(lags); showall=true)
        ac_mat = convert(Matrix, ac)
        val = colordim == :parameter ? ac_mat[:, :, i]' : ac_mat[i, :, :]
        _AutocorPlot(lags, val)
    elseif st ∈ supportedplots
        translationdict[st](c, val)
    else
        range(c), val
    end
end

@recipe function f(p::_DensityPlot)
    xaxis --> "Sample value"
    yaxis --> "Density"
    [collect(skipmissing(p.val[:,k])) for k in 1:size(p.val, 2)]
end

@recipe function f(p::_HistogramPlot)
    xaxis --> "Sample value"
    yaxis --> "Frequency"
    fillalpha --> 0.7
    bins --> 25
    [collect(skipmissing(p.val[:,k])) for k in 1:size(p.val, 2)]
end

@recipe function f(p::_MeanPlot)
    seriestype := :path
    xaxis --> "Iteration"
    yaxis --> "Mean"
    range(p.c), MCMCChains.cummean(p.val)
end

@recipe function f(p::_AutocorPlot)
    seriestype := :path
    xaxis --> "Lag"
    yaxis --> "Autocorrelation"
    p.lags, p.val
end

@recipe function f(p::_TracePlot)
    seriestype := :path
    xaxis --> "Iteration"
    yaxis --> "Sample value"
    range(p.c), p.val
end

@recipe function f(chn::MCMCChains.AbstractChains, parameters::AbstractVector{Symbol};
        colordim = :chain, section = :parameters, append_chains = false)
    c = Chains(chn, section)
    c = append_chains ? pool_chain(c) : c
    colordim != :chain && error("Symbol names are interpreted as parameter names, only compatible with `colordim = :chain`")
    ret = indexin(parameters, Symbol.(keys(c)))
    any(y -> y == nothing, ret) && error("Parameter not found")
    c, Int.(ret)
end

@recipe function f(chn::MCMCChains.AbstractChains,
                   parameters::AbstractVector{<:Integer} = Int[];
                   width = 500,
                   height = 250,
                   colordim = :chain,
                   section = :parameters,
                   append_chains = false
                  )
    c = isempty(parameters) ? Chains(chn, section; sorted=true) : sort(chn)
    c = append_chains ? pool_chain(c) : c
    ptypes = get(plotattributes, :seriestype, (:traceplot, :mixeddensity))
    ptypes = ptypes isa AbstractVector || ptypes isa Tuple ? ptypes : (ptypes,)
    @assert all(map(ptype -> ptype ∈ supportedplots, ptypes))
    ntypes = length(ptypes)
    nrows, nvars, nchains = size(c)
    isempty(parameters) && (parameters = colordim == :chain ? (1:nvars) : (1:nchains))
    N = length(parameters)

    if !(:corner ∈ ptypes)
        layout := (N, ntypes)
        size --> (ntypes*width, N*height)
        indices = reshape(1:N*ntypes, ntypes, N)'
        legend --> false
        for (j, ptype) in enumerate(ptypes)
            for (i, par) in enumerate(parameters)
                @series begin
                    subplot := indices[i, j]
                    colordim := colordim
                    seriestype := ptype
                    c, par
                end
            end
        end
    else
        length(ptypes) > 1 && error(":corner is not compatible with multiple seriestypes")
        Corner(c, names(c)[parameters])
    end
end

struct Corner
    c
    parameters
end

@recipe function f(corner::Corner)
    label --> permutedims(corner.parameters)
    compact --> true
    size --> (600, 600)
    ar = collect(Array(corner.c.value[:, corner.parameters,i]) for i in chains(corner.c))
    RecipesBase.recipetype(:cornerplot, vcat(ar...))
end


text:
# Modified from the file automatically generated by ProtoBuf.jl
# See https://github.com/pszufe/OpenStreetMapX.jl/pull/52/
# syntax: proto2
using ProtoBuf
import ProtoBuf.meta

mutable struct Blob <: SimpleProtoType
    raw_size::Int32
    raw::Vector{UInt8}
    zlib_data::Vector{UInt8}
    lzma_data::Vector{UInt8}
    OBSOLETE_bzip2_data::Vector{UInt8}
    lz4_data::Vector{UInt8}
    zstd_data::Vector{UInt8}
    function Blob()
        return new()
    end
end # mutable struct Blob
const __meta_Blob = Ref{ProtoMeta}()
function meta(::Type{Blob})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Blob)
            __meta_Blob[] = target = ProtoMeta(Blob)
            fnum = Int[2,1,3,4,5,6,7]
            allflds = Pair{Symbol,Union{Type,String}}[:raw_size => Int32, :raw => Vector{UInt8}, :zlib_data => Vector{UInt8}, :lzma_data => Vector{UInt8}, :OBSOLETE_bzip2_data => Vector{UInt8}, :lz4_data => Vector{UInt8}, :zstd_data => Vector{UInt8}]
            oneofs = Int[0,1,1,1,1,1,1]
            oneof_names = Symbol[Symbol("data")]
            meta(target, Blob, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_Blob[]
    end
end

mutable struct BlobHeader <: SimpleProtoType
    _type::AbstractString
    indexdata::Vector{UInt8}
    datasize::Int32
    function BlobHeader()
        return new()
    end
end # mutable struct BlobHeader
const __meta_BlobHeader = Ref{ProtoMeta}()
function meta(::Type{BlobHeader})
    ProtoBuf.metalock() do
        if !isassigned(__meta_BlobHeader)
            __meta_BlobHeader[] = target = ProtoMeta(BlobHeader)
            req = Symbol[:_type,:datasize]
            allflds = Pair{Symbol,Union{Type,String}}[:_type => AbstractString, :indexdata => Vector{UInt8}, :datasize => Int32]
            meta(target, BlobHeader, allflds, req, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_BlobHeader[]
    end
end

export Blob, BlobHeader


text:
using POMDPModels
using Test

let
    problem = MountainCar()
    policy = RandomPolicy(problem)
    sim = RolloutSimulator(rng=MersenneTwister(1), max_steps=1000)

    r = simulate(sim, problem, policy)
    @test r < 0.0

    sv = convert_s(Array{Float64}, (0.5, 0.25), problem)
    @test sv == [0.5, 0.25]
    s = convert_s(Tuple{Float64,Float64}, sv, problem)
    @test s == (0.5, 0.25)

    problem = MountainCar(discount=1.0, cost=-0.1, jackpot=100.0)
    policy = RandomPolicy(problem)
    r = simulate(sim, problem, policy)
    @test r < 0.0
    policy = Energize()
    r = simulate(sim, problem, policy)
    @test r > 0.0
end


text:
# -----------------------------------------------------------
#
# Construction of multidimensional DG derivative matrices
# based on the (precomputed) 1-D derivative matrix
#
# -----------------------------------------------------------

# Efficiency criticality: HIGH
# Matrix computations are only performed once,
# but this can be the main bottleneck if not done right

# Accuracy criticality: HIGH
# Critical for accurate PDE evolution

function D_matrix{D}(i::Int, k::Int, n::Int,
							srefVD::Array{NTuple{3,CartesianIndex{D}},1},
							srefDV::Dict{NTuple{3,CartesianIndex{D}},Int};
							scheme="sparse")

	cutoff = get_cutoff(scheme, D, n)

	V2D_1D = V2Dref(1,k,n)
	D2V_1D = D2Vref(1,k,n)
	Dmat_1D = periodic_DLF_matrix(k, n)

	len = length(srefVD[:,1])
	I = Int[]
	J = Int[]
	V = Float64[]

	for c1 in 1:len
		lpf = srefVD[c1]
		l = lpf[1][i]
		p = lpf[2][i]
		f = lpf[3][i]
		vc1 = D2V_1D[(CartesianIndex(l),CartesianIndex(p),CartesianIndex(f))]
		dvc1s = view(Dmat_1D,:, vc1)

		for j in 1:length(dvc1s)
			lpf2 = V2D_1D[j]
			# using tuple constructor independent of
			# Julia's metaprogramming:
			level2 = make_cartesian_index(i, lpf2[1][1], lpf[1])
			cutoff(level2) && continue

			cell2 = make_cartesian_index(i, lpf2[2][1], lpf[2])
			mode2 = make_cartesian_index(i, lpf2[3][1], lpf[3])
			c2 = srefDV[(level2, cell2, mode2)]
			push!(I, c2)
			push!(J, c1)
			push!(V, dvc1s[j])
		end
	end
	# dropzeros! does not seem helpful for this matrix:
	return sparse(I, J, V, len, len, +)
end

function D_matrix(D::Int, i::Int, k::Int, n::Int; scheme="sparse")
	# Precompute the 1D derivatve matrix elements as global variables
	# if they are not yet formed
	precompute_diffs()
	VD = V2Dref(D, k, n; scheme=scheme)
	DV = D2Vref(D, k, n; scheme=scheme)
	return D_matrix(i, k, n, VD, DV; scheme=scheme)
end

function grad_matrix(D::Int, k::Int, n::Int; scheme="sparse")
	return [D_matrix(D, i, k, n; scheme=scheme) for i in 1:D]
end

function laplacian_matrix(D::Int, k::Int, n::Int; scheme="sparse")
	len = get_size(D, k, n; scheme=scheme)
	laplac = spzeros(len, len)
	for i in 1:D
		D_op = D_matrix(D, i, k, n; scheme=scheme)
		laplac += *(D_op, D_op)
	end
	return laplac
end


text:
include("common.jl")

function f(x, y)
    0.0001*(10/(1+x^2) + x * y + 10*y^2)
end

Δt = 0.01
mmesh = Mesh(30, 30, 1/30)
bdnode = bcnode(mmesh)
bdnode = [bdnode; bdnode .+ mmesh.ndof]
# nu = constant(eval_f_on_gauss_pts(f, mmesh))

nu = Variable(ones(get_ngauss(mmesh))) * 0.0001

nodes = fem_nodes(mmesh)
u = @. sin(2π * nodes[:,1])
v = @. sin(2π * nodes[:,2])
u0 = [u;v]
u0[bdnode] .= 0.0
us = solve_burgers(u0, 10, nu)

U = matread("fenics/fwd2.mat")["U"]
loss = sum((us - U)^2)

_loss = Float64[]
cb = (vs, iter, loss)->begin
    global _loss
    push!(_loss, loss)
    printstyled("[#iter $iter] loss=$loss\n", color=:green)
    if mod(iter, 10)==1
        @info iter
        close("all")
        visualize_scalar_on_gauss_points(vs[1], mmesh)
        matwrite("fenics/bwd3-$iter.mat", Dict("iter"=>iter,"loss"=>_loss, "nu"=>vs[1]))
        savefig("fenics/bwd3_pixel$iter.png")
    end
end

# @info run(sess, loss)
sess = Session(); init(sess)
BFGS!(sess, loss, vars = [nu], callback = cb)


text:
using QuantumOptics
using Test

@testset "bloch-redfield" begin

Δ = 0.2 * 2*π
ϵ0 = 1.0 * 2*π
γ1 = 0.5

b = SpinBasis(1//2)
sx = sigmax(b)
sz = sigmaz(b)

H = -Δ/2.0 * sx - ϵ0/2.0 * sz

function ohmic_spectrum(ω)
    ω == 0.0 ? (return γ1) : (return γ1/2 * (ω / (2*π)) * (ω > 0.0))
end

R, ekets = timeevolution.bloch_redfield_tensor(H, [[sx, ohmic_spectrum]])

known_result =  [0.0+0.0im        0.0+0.0im            0.0+0.0im       0.245145+0.0im
 0.0+0.0im  -0.161034-6.40762im        0.0+0.0im            0.0+0.0im
 0.0+0.0im        0.0+0.0im      -0.161034+6.40762im        0.0+0.0im
 0.0+0.0im        0.0+0.0im            0.0+0.0im      -0.245145+0.0im]
@test isapprox(dense(R).data, known_result, atol=1e-5)

psi0 = spindown(b)
tout, ρt = timeevolution.master_bloch_redfield([0.0:0.1:2.0;], psi0, R, H)
rho_end = [0.38206-0.0im 0.0466443+0.0175017im
            0.0466443-0.0175017im 0.61794+0.0im]
@test isapprox(ρt[end].data, rho_end, atol=1e-5)
@test ρt[end] != ρt[end-1]
@test isa(ρt, Vector{<:DenseOperator})

# Test fout
fout(t,rho) = copy(rho)
tout, ρt2 = timeevolution.master_bloch_redfield([0.0:0.1:2.0;], psi0, R, H; fout=fout)
@test all(ρt .== ρt2)

fout2(t,rho) = expect(sz,rho)
tout, z = timeevolution.master_bloch_redfield([0.0:0.1:2.0;], psi0, R, H; fout=fout2)
@test length(z) == length(ρt)
@test isa(z, Vector{ComplexF64})

end # testset


text:
# Helper function

function zip_url_to_df(url)
    zip_response = HTTP.get(url)
    zipb = zip_response.body
    ziparchive = ZipFile.Reader(IOBuffer(zipb, read=true, write=true))
    return df = CSV.File(ziparchive.files[1]) |> DataFrame
end

"""
    fundamentals(api_key)

Download the CORE US Fundamentals dataset from Sharadar via Quandl.

Equity fundamentals data. Equity fundamentals are financial statement
information (cash flow, earnings statements, balance sheets).

To see an example, [here](https://uk.finance.yahoo.com/quote/MANU/financials?p=MANU) you can find the fundamentals data for Manchester United
football club (MANU).

# Examples
```julia-repl
julia> api_key = "-31swyRz8v4SoNL_zyYc"
julia> fundamentals(api_key)

```

# Related data

The dataset generated by `fundamentals` lists firms by their stock tickers
(`AAPL` instead of Apple), and uses short codes to indicate the fundamentals
(`rnd` instead of Research and Development).

To download descriptions of the indicator codes, use
```julia-repl
julia> fundamental_indicators(api_key)

```
Or if you just need to know a couple of indicators, go directly to the table
[here](https://www.quandl.com/tables/SF1/SHARADAR-INDICATORS)

To download a database of the tickers, company names and other information
including industry and listing status, use
```julia-repl
julia> tickers(api_key)

```
"""
function fundamentals(api_key)
    url  = "https://www.quandl.com/api/v3/datatables/SHARADAR/SF1.csv?&qopts.export=true&dimension=MRY&api_key="
    println("Downloading the following data:")
    println("Dataset: Sharadar CORE US Equities Bundle")
    println("Table:   CORE US Fundamentals")
    println("Source:  Quandl")
    println("URL:     https://www.quandl.com/tables/SFA/SHARADAR-SF1")
    try
        println("Generating download link")
        df_  = CSV.File(HTTP.get(url*api_key).body) |> DataFrame
        url_ = df_[!,Symbol("file.link")][1]
        println("Downloading the dataset")
        df = zip_url_to_df(url_)
        return df
    catch
        println("I couldn't download the dataset. This could be due to one of the ")
        println("following issues:")
        println("")
        println(" 1. Make sure that you provide your API key, with double quote")
        println("    marks as follows:")
        println("          fundamentals(\"api_key_goes_here\")")
        println("    Your API key will look something like -31swyRz8v4SoNL_zyYc ")
        println("    If you haven't received an API key for Quandl, inform your ")
        println("    instructor.")
        println(" 2. Make sure that you're connected to the internet.")
        println(" 3. Make sure that you can login to Quandl.com. ")
        println("    Check by going to the following URL:")
        println("    https://www.quandl.com/databases/SFA/data")
        println("    and check that you can open the datasets in your browser.")
    end
end

function fundamentals()
    println("Please provide your Quandl API key as an argument to the function:")
    println("      fundamentals(\"api_key_goes_here\")")
    println("Your API key will look something like -31swyRz8v4SoNL_zyYc ")
    println("Make sure that you include the double quote marks as above.")
end

"""
    events(api_key)

Download the CORE US Fundamentals Events dataset from Sharadar via Quandl.

Equity events data. Equity events relate to SEC filings, where companies
publicly announce major events including financial restatements, missed
debt payments and so on.

# Examples
```julia-repl
julia> api_key = "-31swyRz8v4SoNL_zyYc"
julia> events(api_key)

```

# Related data

The dataset generated by `events` lists firms by their stock tickers
(`AAPL` instead of Apple), and uses short codes to indicate the events
(`13` instead of Bankruptcy or Receivership).

A full list of event codes can be found [here](https://www.sec.gov/fast-answers/answersform8khtm.html).

To download a database of the tickers, company names and other information
including industry and listing status, use
```julia-repl
julia> tickers(api_key)

```
"""
function events(api_key)
    url  = "https://www.quandl.com/api/v3/datatables/SHARADAR/EVENTS.csv?&qopts.export=true&api_key="
    println("Downloading the following data:")
    println("Dataset: Sharadar CORE US Equities Bundle")
    println("Table:   CORE US Fundamentals Events")
    println("Source:  Quandl")
    println("URL:     https://www.quandl.com/tables/SFA/SHARADAR-EVENTS")
    try
        println("Generating download link")
        df_  = CSV.File(HTTP.get(url*api_key).body) |> DataFrame
        url_ = df_[!,Symbol("file.link")][1]
        println("Downloading the dataset")
        df = zip_url_to_df(url_)
        df[!, :eventcodes] = split.(df[:,:eventcodes], '|');
        df = flatten(df, :eventcodes)
        return df
    catch
        println("I couldn't download the dataset. This could be due to one of the ")
        println("following issues:")
        println("")
        println(" 1. Make sure that you provide your API key, with double quote")
        println("    marks as follows:")
        println("          events(\"api_key_goes_here\")")
        println("    Your API key will look something like -31swyRz8v4SoNL_zyYc ")
        println("    If you haven't received an API key for Quandl, inform your ")
        println("    instructor.")
        println(" 2. Make sure that you're connected to the internet.")
        println(" 3. Make sure that you can login to Quandl.com. ")
        println("    Check by going to the following URL:")
        println("    https://www.quandl.com/tables/SFA/SHARADAR-EVENTS")
        println("    and check that you can open the datasets in your browser.")
    end
end

function events()
    println("Please provide your Quandl API key as an argument to the function:")
    println("      events(\"api_key_goes_here\")")
    println("Your API key will look something like -31swyRz8v4SoNL_zyYc ")
    println("Make sure that you include the double quote marks as above.")
end

function splitrows(df_::DataFrame,col::Symbol,delim::Char;sortcols=names(df_))
    df = deepcopy(df_)
    while .!isnothing.(findfirst.(delim,df[:,col])) |> sum > 0
        filter  = isnothing.(findfirst.(delim,df[:,col]))
        df_temp = df[.!filter,:]
        df      = df[filter,:]
        n = size(df_temp)[1]
        for i in 1:n
            values_ = df_temp[i,col]
            delim_index = findfirst(delim,values_)
            df_temp[i,col] = values_[1:(delim_index-1)]
            push!(df_temp,
                (ticker     = df_temp[i,:ticker],
                 date       = df_temp[i,:date],
                 eventcodes = values_[(delim_index+1):end]))
        end
        df = vcat(df,df_temp)
    end
    sort!(df,sortcols)
    return df
end

function splitrows_(df_::DataFrame,col::Symbol,delim::Char;sortcols=names(df_))
    df = deepcopy(df_)
    df[!, col] = split.(df[:,col], delim);
    df = flatten(df, col)
    return df
end

"""
    metadata(api_key)

Download the Tickers and Metadata dataset from Sharadar via Quandl.

Metadata includes stock tickers, company names, industry codes, and related
information.

# Examples
```julia-repl
julia> api_key = "-31swyRz8v4SoNL_zyYc"
julia> metadata(api_key)

```
"""
function metadata(api_key)

    url  = "https://www.quandl.com/api/v3/datatables/SHARADAR/TICKERS.csv?&qopts.export=true&api_key="
    println("Downloading the following data:")
    println("Dataset: Sharadar CORE US Equities Bundle")
    println("Table:   Tickers and Metadata")
    println("Source:  Quandl")
    println("URL:     https://data.nasdaq.com/tables/SFA/SHARADAR-TICKERS")
    try
        println("Generating download link")
        df_  = CSV.File(HTTP.get(url*api_key).body) |> DataFrame
        url_ = df_[!,Symbol("file.link")][1]
        println("Downloading the dataset")
        df = zip_url_to_df(url_)
        return df
    catch
        println("I couldn't download the dataset. This could be due to one of the ")
        println("following issues:")
        println("")
        println(" 1. Make sure that you provide your API key, with double quote")
        println("    marks as follows:")
        println("          metadata(\"api_key_goes_here\")")
        println("    Your API key will look something like -31swyRz8v4SoNL_zyYc ")
        println("    If you haven't received an API key for Quandl, inform your ")
        println("    instructor.")
        println(" 2. Make sure that you're connected to the internet.")
        println(" 3. Make sure that you can login to Quandl.com. ")
        println("    Check by going to the following URL:")
        println("    https://www.quandl.com/databases/SFA/data")
        println("    and check that you can open the datasets in your browser.")
    end
end

function metadata()
    println("Please provide your Quandl API key as an argument to the function:")
    println("      metadata(\"api_key_goes_here\")")
    println("Your API key will look something like -31swyRz8v4SoNL_zyYc ")
    println("Make sure that you include the double quote marks as above.")
end

"""
    stockprices(api_key)

Download the Equity Prices dataset from Sharadar via Quandl.

# Examples
```julia-repl
julia> api_key = "-31swyRz8v4SoNL_zyYc"
julia> stockprices(api_key)

```
"""

function stockprices(
    api_key::String;ticker="",
    start_date=Dates.today() - Dates.Year(1)
    )

    url = "https://www.quandl.com/api/v3/datatables/SHARADAR/SEP.csv?&qopts.export=true"*
            (ticker == "" ? "" : "&ticker="*ticker)*
            "&date.gte="*Dates.format(start_date,"yyyy-mm-dd")*
            "&api_key="*
            api_key

    println("api call = "*url)
    println("Downloading the following data:")
    println("Dataset: Sharadar CORE US Equities Bundle")
    println("Table:   Equity Prices")
    println("Source:  Quandl")
    println("URL:     https://data.nasdaq.com/tables/SFA/SHARADAR-SEP")
    println("This will take a few minutes.")
    try
        println("Generating download link")
        df_  = CSV.File(HTTP.get(url).body) |> DataFrame
        url_ = df_[!,Symbol("file.link")][1]
        println("Downloading the dataset")
        df = zip_url_to_df(url_)
        sort!(df,[:ticker,:date])
        return df
    catch
        println("I couldn't download the dataset. This could be due to one of the ")
        println("following issues:")
        println("")
        println(" 1. Make sure that you provide your API key, with double quote")
        println("    marks as follows:")
        println("          stockprices(\"api_key_goes_here\")")
        println("    Your API key will look something like -31swyRz8v4SoNL_zyYc ")
        println("    If you haven't received an API key for Quandl, inform your ")
        println("    instructor.")
        println(" 2. Make sure that you're connected to the internet.")
        println(" 3. Make sure that you can login to Quandl.com. ")
        println("    Check by going to the following URL:")
        println("    https://www.quandl.com/databases/SFA/data")
        println("    and check that you can open the datasets in your browser.")
    end
end

function stockprices()
    println("Please provide your Quandl API key as an argument to the function:")
    println("      stockprices(\"api_key_goes_here\")")
    println("Your API key will look something like -31swyRz8v4SoNL_zyYc ")
    println("Make sure that you include the double quote marks as above.")
end

"""
    fundprices(api_key)

Download the Equity Prices dataset from Sharadar via Quandl.

# Examples
```julia-repl
julia> api_key = "-31swyRz8v4SoNL_zyYc"
julia> fundprices(api_key)

```
"""

function fundprices(
    api_key::String;
    ticker="",
    start_date=Dates.today() - Dates.Year(1)
    )

    url = "https://www.quandl.com/api/v3/datatables/SHARADAR/SFP.csv?&qopts.export=true"*
            (ticker == "" ? "" : "&ticker="*ticker)*
            "&date.gte="*Dates.format(start_date,"yyyy-mm-dd")*
            "&api_key="*
            api_key

    println("api call = "*url)
    println("Downloading the following data:")
    println("Dataset: Sharadar CORE US Equities Bundle")
    println("Table:   Equity Prices")
    println("Source:  Quandl")
    println("URL:     https://data.nasdaq.com/tables/SFA/SHARADAR-SFP")
    println("This will take a few minutes.")
    try
        println("Generating download link")
        df_  = CSV.File(HTTP.get(url).body) |> DataFrame
        url_ = df_[!,Symbol("file.link")][1]
        println("Downloading the dataset")
        df = zip_url_to_df(url_)
        sort!(df,[:ticker,:date])
        return df
    catch
        println("I couldn't download the dataset. This could be due to one of the ")
        println("following issues:")
        println("")
        println(" 1. Make sure that you provide your API key, with double quote")
        println("    marks as follows:")
        println("          fundprices(\"api_key_goes_here\")")
        println("    Your API key will look something like -31swyRz8v4SoNL_zyYc ")
        println("    If you haven't received an API key for Quandl, inform your ")
        println("    instructor.")
        println(" 2. Make sure that you're connected to the internet.")
        println(" 3. Make sure that you can login to Quandl.com. ")
        println("    Check by going to the following URL:")
        println("    https://www.quandl.com/databases/SFA/data")
        println("    and check that you can open the datasets in your browser.")
    end
end

function fundprices()
    println("Please provide your Quandl API key as an argument to the function:")
    println("      fundprices(\"api_key_goes_here\")")
    println("Your API key will look something like -31swyRz8v4SoNL_zyYc ")
    println("Make sure that you include the double quote marks as above.")
end

max_stars_repo_name: alfredjmduncan/KentFinancialEconomics

max_issues_repo_name: alfredjmduncan/KentFinancialEconomics

max_forks_repo_name: alfredjmduncan/KentFinancialEconomics


text:
"Data type for candidate solutions from master problem."
struct CandSol
    zsol::Array{Bool,2}
    qsol::Vector{Float64}
    ϕsol::Vector{Float64}
end

"Result data from GBD algorithm"
struct Result
    status::MOI.TerminationStatusCode
    solution
    primalbound::Float64
    dualbound::Float64
    niter::Int
end

"""
Extract sub-topology from active arcs and incident nodes.

The option `keep_nodes` allows to also include the isolated nodes.
"""
function topology_from_candsol(topo::Topology, ysol::Vector{Float64},
                               keep_nodes=false)
    actarcs = topo.arcs[ysol .> 0.5]
    if keep_nodes
        return Topology(topo.nodes, actarcs)
    end
    tails = [arc.tail for arc in actarcs]
    heads = [arc.head for arc in actarcs]
    nodeidx = tails ∪ heads
    nodemap = fill(0, length(topo.nodes))
    for (target, idx) in enumerate(nodeidx)
        nodemap[idx] = target
    end
    Topology(topo.nodes[nodeidx],
             [Arc(nodemap[arc.tail], nodemap[arc.head]) for arc in actarcs])
end

"Generic no-good cut"
function nogood(model::JuMP.Model, vars::AbstractArray{JuMP.VariableRef},
                sol::AbstractArray{T}) where T<:Real
    @assert size(vars) == size(sol)
    nvars = length(vars)
    active = (sol .> 0.5)
    nactive = sum(active)
    coef = 2.0*active .- 1.0
    @constraint(model, sum(coef[i]*vars[i] for i=1:nvars) <= nactive - 1)
    return 1
end

"Cut off all y values for given, undirected topology"
function avoid_topo_cut(model::JuMP.Model, y::AbstractArray{JuMP.VariableRef},
                        topo::Topology, edges::Vector{Arc})
    arcidx = arcindex(topo)
    antidx = antiparallelindex(topo)

    fwd = [arcidx[e] for e in edges]
    bwd = [antidx[a] for a in fwd]
    arcs = vcat(fwd, bwd)
    @constraint(model, sum(y[a] for a in arcs) ≤ length(edges) - 1)
end


text:
# Enable CUDA if the CUDA package is installed (enables to use the package for CPU-only without requiring the CUDA package installed if the installation procedure allows it).
const CUDA_IS_INSTALLED = (Base.find_package("CUDA")!==nothing)
const ENABLE_CUDA = true # NOTE: Can be set to CUDA_IS_INSTALLED, or to true or false independent of it.
const PKG_CUDA = :CUDA
const PKG_THREADS = :Threads
const PKG_NONE = :PKG_NONE
@static if ENABLE_CUDA
    using CUDA
    const SUPPORTED_PACKAGES = [PKG_THREADS, PKG_CUDA]
else
    const SUPPORTED_PACKAGES = [PKG_THREADS]
end
using MacroTools
import MacroTools: postwalk, isexpr, inexpr

## CONSTANTS

const GENSYM_SEPARATOR = ", "
gensym_world(tag::String, generator::Module) = gensym(string(tag, GENSYM_SEPARATOR, generator)) #NOTE: this function needs to be defind before constants using it.
gensym_world(tag::Symbol, generator::Module) = gensym(string(tag, GENSYM_SEPARATOR, generator))

const INT_CUDA = Int64
const NTHREADS_MAX = 256
const RANGES_VARNAME = gensym_world("ranges", @__MODULE__)
const RANGELENGTHS_VARNAMES = (gensym_world("rangelength_x", @__MODULE__), gensym_world("rangelength_y", @__MODULE__), gensym_world("rangelength_z", @__MODULE__))
const THREADIDS_VARNAMES = (gensym_world("tx", @__MODULE__), gensym_world("ty", @__MODULE__), gensym_world("tz", @__MODULE__))
const RANGES_TYPE_1D = UnitRange{}
const RANGES_TYPE_1D_TUPLE = Tuple{UnitRange{}}
const RANGES_TYPE_2D = Tuple{UnitRange{},UnitRange{}}
const RANGES_TYPE = Tuple{UnitRange{},UnitRange{},UnitRange{}}
const MAXSIZE_TYPE_1D = Integer
const MAXSIZE_TYPE_1D_TUPLE = Tuple{T} where T <: Integer
const MAXSIZE_TYPE_2D = Tuple{T, T} where T <: Integer
const MAXSIZE_TYPE = Tuple{T, T, T} where T <: Integer
const BOUNDARY_WIDTH_TYPE_1D = Integer
const BOUNDARY_WIDTH_TYPE_1D_TUPLE = Tuple{T} where T <: Integer
const BOUNDARY_WIDTH_TYPE_2D = Tuple{T, T} where T <: Integer
const BOUNDARY_WIDTH_TYPE = Tuple{T, T, T} where T <: Integer
const OPERATORS = [:-, :+, :*, :/, :%, :!, :&&, :||] #NOTE: ^ is not contained as causes an error.
const SUPPORTED_LITERALTYPES = [Float16, Float32, Float64, Complex{Float16}, Complex{Float32}, Complex{Float64}, Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8] # NOTE: Not isbitstype as required for CUDA: BigFloat, BigInt, Complex{BigFloat}, Complex{BigInt}
const SUPPORTED_NUMBERTYPES  = [Float16, Float32, Float64, Complex{Float16}, Complex{Float32}, Complex{Float64}]
const NUMBERTYPE_NONE = DataType
const ERRMSG_UNSUPPORTED_PACKAGE = "unsupported package for parallelization"
const ERRMSG_CHECK_PACKAGE  = "package has to be one of the following: $(join(SUPPORTED_PACKAGES,", "))"
const ERRMSG_CHECK_NUMBERTYPE = "numbertype has to be one of the following: $(join(SUPPORTED_NUMBERTYPES,", "))"
const ERRMSG_CHECK_LITERALTYPES = "the type given to 'literaltype' must be one of the following: $(join(SUPPORTED_LITERALTYPES,", "))"

## FUNCTIONS TO DEAL WITH FUNCTION DEFINITIONS: SIGNATURES, BODY AND RETURN STATEMENT

function push_to_signature!(kernel::Expr, arg::Expr)
    push!(kernel.args[1].args, arg)
    return kernel
end

function get_body(kernel::Expr)
    return kernel.args[2]
end

function set_body!(kernel::Expr, body::Expr)
    kernel.args[2] = body
end

function remove_return(body::Expr)
    if !(body.args[end] in [:(return), :(return nothing), :(nothing)])
        @ArgumentError("invalid kernel in @parallel kernel definition: the last statement must be a `return nothing` statement ('return' or 'return nothing' or 'nothing') as required for any CUDA kernels.")
    end
    remainder = copy(body)
    remainder.args = body.args[1:end-2]
    if inexpr_walk(remainder, :return) @ArgumentError("invalid kernel in @parallel kernel definition: only one return statement is allowed in the kernel and it must return nothing and be the last statement (required to ensure equal behaviour with different packages for parallellization).") end
    return remainder
end

function add_return(body::Expr)
    quote
        $body
        return nothing
    end
end

## FUNCTIONS TO DEAL WITH FUNCTION/MACRO CALLS: CHECK IF DEFINITION/CALL, EXTRACT, SPLIT AND EVALUATE ARGUMENTS

is_function(arg)    = ( isa(arg, Expr) && ( (arg.head == :function) || (arg.head == :(=) && isa(arg.args[1], Expr) && arg.args[1].head == :call) ) )
is_call(arg)        = ( isa(arg, Expr) && (arg.head == :call) )
is_block(arg)       = ( isa(arg, Expr) && (arg.head == :block) )
is_parallel_call(x) = isexpr(x, :macrocall) && (x.args[1] == Symbol("@parallel") || x.args[1] == :(@parallel))

macro get_args(args...) return args end
extract_args(call::Expr, macroname::Symbol) = eval(substitute(deepcopy(call), macroname, Symbol("@get_args")))

function split_args(args)
    posargs   = [x for x in args[1:end-1] if !(isa(x,Expr) && x.head == :(=))]
    kwargs    = [x for x in args[1:end-1] if isa(x,Expr) && x.head == :(=)]
    kernelarg = args[end]
    if any([x.args[1] in [:blocks, :threads] for x in kwargs]) @KeywordArgumentError("Invalid keyword argument in @parallel call: blocks / threads. They must be passed as positional arguments or been omited.") end
    return posargs, kwargs, kernelarg
end

function eval_arg(caller::Module, arg)
    try
        return @eval(caller, $arg)
    catch e
        @ArgumentEvaluationError("argument $arg could not be evaluated at parse time (in module $caller).")
    end
end

## FUNCTIONS FOR COMMON MANIPULATIONS ON EXPRESSIONS

function substitute(expr::Expr, old, new)
    postwalk(expr) do x
        if x == old
            return new
        else
            return x;
        end
    end
end

function inexpr_walk(expr::Expr, s::Symbol)
    found = false
    postwalk(expr) do x
        if ((isa(x,Expr) && x.head==s) || (isa(x,Symbol) && x==s)) found = true end
        return x
    end
    return found
end

## FUNCTIONS FOR ERROR HANDLING

check_package(P)                = ( if !isa(P, Symbol) || !(P in SUPPORTED_PACKAGES)  @ArgumentError("$ERRMSG_CHECK_PACKAGE (obtained: $P)." ) end )
check_numbertype(T::DataType)   = ( if !(T in SUPPORTED_NUMBERTYPES) @ArgumentError("$ERRMSG_CHECK_NUMBERTYPE (obtained: $T)." ) end )
check_literaltype(T::DataType)  = ( if !(T in SUPPORTED_LITERALTYPES) @ArgumentError("$ERRMSG_CHECK_LITERALTYPES (obtained: $T)." ) end )
check_numbertype(datatypes...)  = check_numbertype.(datatypes)
check_literaltype(datatypes...) = check_literaltype.(datatypes)

## FUNCTIONS AND MACROS FOR UNIT TESTS

symbols(eval_mod::Union{Symbol,Module}, mod::Union{Symbol,Module}) = @eval(eval_mod, names($mod, all=true, imported=true))
prettystring(expr::Expr)                                           = string(remove_linenumbernodes!(expr))
gorgeousstring(expr::Expr)                                         = string(simplify_varnames!(remove_linenumbernodes!(expr)))
longnameof(f)                                                      = "$(parentmodule(f)).$(nameof(f))"
macro require(condition)               condition_str = string(condition); esc(:( if !($condition) error("pre-test requirement not met: $($condition_str).") end )) end  # Verify a condition required for a unit test (in the unit test results, this should not be treated as a unit test).
macro symbols(eval_mod, mod)           symbols(eval_mod, mod) end
macro macroexpandn(n::Integer, expr)   return QuoteNode(macroexpandn(__module__, expr, n)) end
macro prettyexpand(n::Integer, expr)   return QuoteNode(remove_linenumbernodes!(macroexpandn(__module__, expr, n))) end
macro gorgeousexpand(n::Integer, expr) return QuoteNode(simplify_varnames!(remove_linenumbernodes!(macroexpandn(__module__, expr, n)))) end
macro prettyexpand(expr)               return QuoteNode(remove_linenumbernodes!(macroexpand(__module__, expr; recursive=true))) end
macro gorgeousexpand(expr)             return QuoteNode(simplify_varnames!(remove_linenumbernodes!(macroexpand(__module__, expr; recursive=true)))) end

function macroexpandn(m::Module, expr, n::Integer)
    for i = 1:n
        expr = macroexpand(m, expr; recursive=false)
    end
    return expr
end

function remove_linenumbernodes!(expr::Expr)
    expr = Base.remove_linenums!(expr)
    args = expr.args
    for i=1:length(args)
        if isa(args[i], LineNumberNode)
             args[i] = nothing
        elseif typeof(args[i]) == Expr
            args[i] = remove_linenumbernodes!(args[i])
        end
    end
    return expr
end

function simplify_varnames!(expr::Expr)
    args = expr.args
    for i=1:length(args)
        if isa(args[i], Symbol)
            varname = string(args[i]);
            if startswith(varname, "##")
                varname = replace(varname, "##" => "")
                varname = replace(varname, r"#\d*" => "")
                varname = split(varname, GENSYM_SEPARATOR)[1]
                args[i] = Symbol(varname)
            end
        elseif isa(args[i], Expr)
            args[i] = simplify_varnames!(args[i])
        end
    end
    return expr
end


text:
function simulate_exp_gpreal(t,alpha,beta,ndev)
#
# Function for simulating a GP based upon an exponential correlation
# matrix using an analytic form for the Cholesky decomposition and
# a computation in O(N) operations.  The autocorrelation function is:
#
#    K_ij = alpha * exp(-beta*(|t_i-t_j|))
#
# Requirements:
#  - alpha should be a real positive number.
#  - The times t *must* be sorted in order from least to greatest.
#  - beta may be complex
#  - ndev: normal deviates drawn from N(0,1) with the same length as time vector t.
# Output:
#  - data is a GP drawn from this correlation function with length nt
#
nt = length(t)
data = zeros(eltype(beta),nt)
data[nt] = sqrt(alpha)*ndev[nt]
gamma = zero(eltype(beta))
for i=nt-1:-1:1
  gamma = exp(-beta*(t[i+1]-t[i]))
#  println(i," ",abs(sqrt(1.0-gamma^2)))
  data[i] = sqrt(1.0-gamma^2)*sqrt(alpha)*ndev[i]+gamma*data[i+1]
end
return data
end


text:
using AbstractPlotting
using MakieRipserer
using Ripserer
using Suppressor
using Test

@testset "Nothing throws" begin
    data = [(sin(t), cos(t)) for t in range(0, 2π; length=13)[1:12]]
    time = Observable(0.0)
    @suppress begin
        @test begin
            MakieRipserer.app(data, Rips; time)
            time[] = 1
            true
        end
    end
end


text:
using Colors: colormap
using RecipesBase

"shift and scale numbers to interval [lb=0, ub=1]"
function shiftnscale(data::AbstractArray{T}; lb=zero(T), ub=one(T)) where T<:Real
    @assert ub > lb
    mi, ma = extrema(data)
    range = ma == mi ? 1.0 : ma - mi
    (data - mi) / range * (ub - lb) + lb
end

"map numbers to colors"
function map2color(data::AbstractArray{T}; ncolors::Int=20, cmap="Blues") where T<:Real
    data = shiftnscale(data)

    # take two extra colors, one to leave out at each extremity, and another one
    # because the normalized data is in [0, 1] inclusive
    colors = colormap(cmap, ncolors + 3)
    indices = round.(Int, ncolors*data + 2, RoundDown)
    colors[indices]
end

# a recipe to display topologies with arcs and nodes
@recipe function f(nodes::AbstractArray{PipeLayout.Node},
                   arcs::AbstractArray{PipeLayout.Arc})
    # prepare node data
    nodepos = reshape(reinterpret(Float64, nodes), (2, length(nodes)))'
    x = nodepos[:, 1]
    y = nodepos[:, 2]

    # prepare arc data
    arcidx = reshape(reinterpret(Int64, arcs), (2, length(arcs)))

    # first draw the arcs
    @series begin
        color --> :black
        x[arcidx], y[arcidx]
    end

    # then the nodes
    @series begin
        seriestype := :scatter
        aspect_ratio --> :equal
        markercolor --> :white
        markersize --> 3
        grid --> false
        foreground_color_border	--> :white
        legend --> false
        ticks --> nothing
        x, y
    end
end

# a recipe to display only nodes
@recipe function f(nodes::AbstractArray{PipeLayout.Node})
    markeralpha --> 0.5
    markercolor --> :red
    markersize --> 6
    nodes, Arc[]
end

# a recipe for topologies: just nodes and arcs
@recipe function f(topo::PipeLayout.Topology)
    topo.nodes, topo.arcs
end

# a recipe to display an instance: nodes with demand
@recipe function f(instance::PipeLayout.Instance)
    nodes = instance.nodes
    demand = instance.demand

    colors = fill(:gray, length(nodes))
    colors[demand .> 0] = :red
    colors[demand .< 0] = :blue

    sizes = sqrt.(abs.(demand))
    sizes = 10 * sizes / maximum(sizes) + 4

    markercolor --> colors
    markersize --> sizes
    nodes, Arc[]
end

# a recipe to display a candidate solution: active arcs with diameter choices
@recipe function f(topo::PipeLayout.Topology,
                   cand::PipeLayout.GndStr.CandSol)
    narcs, ndiams = size(cand.zsol)
    arcidx, diamidx = findn(cand.zsol)
    arcs = topo.arcs[arcidx]

    linewidth --> 3 * diamidx'
    linecolor --> map2color(diamidx', ncolors=ndiams, cmap="Greens")

    # need to redraw nodes, because arcs reference them
    markersize --> 0

    topo.nodes, arcs
end

# a recipe to display a pipe dimensioning solution: diameter choices
@recipe function f(inst::PipeLayout.Instance,
                   topo::PipeLayout.Topology,
                   sol::PipeLayout.PipeDim.Solution)
    narcs, ndiams = size(sol.lsol)
    equiv = [serialmerge(inst.diameters, sol.lsol[a,:]) for a=1:narcs]

    linewidth --> 8 * equiv'
    linecolor --> map2color(equiv', ncolors=ndiams, cmap="Greens")

    # need to redraw nodes, because arcs reference them
    markersize --> 0

    topo.nodes, topo.arcs
end

# a recipe to display a junction location solution: nodes & diameter choices
@recipe function f(inst::PipeLayout.Instance,
                   topo::PipeLayout.Topology,
                   sol::PipeLayout.JuncLoc.Solution)
    nodes = sol.nodes
    narcs, ndiams = size(sol.lsol)
    equiv = [serialmerge(inst.diameters, sol.lsol[a,:]) for a=1:narcs]

    linewidth --> 8 * equiv'
    linecolor --> map2color(equiv', ncolors=ndiams, cmap="Greens")

    # need to redraw nodes, because arcs reference them.
    # also: show junction nodes (not in Instance)
    markersize --> 8

    sol.nodes, topo.arcs
end


text:
using DeepGetfield, MacroTools
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

struct A
    x::Int
    y
    z
end

struct B
    xx
    yy
end

struct C
    zz
end

a = A(2, 1., B(10.,[C(1),C(2)]))

ops = findfield(a, :yy)
@test @capture(ops[1], x->getfield(x, arg_))
@test @capture(ops[2], x->getfield(x, arg_))

ops_b = DeepGetfield.tobroadcast(ops)
@test @capture(ops_b[1], x->getfield(x, arg_))
@test @capture(ops_b[2], x->getfield(x, arg_))

fun = getter(a, :yy)
@test fun(a) isa Vector{C}

ops = findfield(a, :zz)
@test @capture(ops[3], x->getindex(x, 1))
@test @capture(ops[4], x->getfield(x, arg_))

ops_b = DeepGetfield.tobroadcast(ops)
@test length(ops_b) == length(ops) - 1 == 3
@test @capture(ops_b[1], x->getfield(x, arg_))
@test @capture(ops_b[2], x->getfield(x, arg_))
@test @capture(ops_b[3], x->getfield.(x, arg_))

fun2 = getter(a, :zz)
@test fun2(a) == [1,2] == a.z.yy |> x -> getfield.(x,:zz)

@deep(a.yy) == fun(a)
@deep(a.zz) == [1,2]

ops = [:(x-> getindex(x,1))]
b = [1,2]
@test DeepGetfield.follow(b,ops) == 1

f = eval(ops[1])
@test Base.invokelatest(f,b) == 1

x = [1, 2]
y = [x, x, x]
z = [y, y, y, y]

@test DeepGetfield.calcdims(x) == [length(x)]
@test DeepGetfield.calcdims(y) == [length(y), length(x)]
@test DeepGetfield.calcdims(z) == [length(z), length(y), length(x)]
@test denest(x) == x
@test denest(y) == [x x x]
@test denest(z) == cat(3, denest(y), denest(y), denest(y), denest(y))