text:
const DISTMSG = "Right-hand side of a ~ must be subtype of Distribution or a vector of " *
    "Distributions."

const INTERNALNAMES = (:_model, :_sampler, :_context, :_varinfo)

"""
    isassumption(expr)

Return an expression that can be evaluated to check if `expr` is an assumption in the
model.

Let `expr` be `:(x[1])`. It is an assumption in the following cases:
    1. `x` is not among the input data to the model,
    2. `x` is among the input data to the model but with a value `missing`, or
    3. `x` is among the input data to the model with a value other than missing,
       but `x[1] === missing`.

When `expr` is not an expression or symbol (i.e., a literal), this expands to `false`.
"""
function isassumption(expr::Union{Symbol, Expr})
    vn = gensym(:vn)

    return quote
        let $vn = $(varname(expr))
            # This branch should compile nicely in all cases except for partial missing data
            # For example, when `expr` is `:(x[i])` and `x isa Vector{Union{Missing, Float64}}`
            if !$(DynamicPPL.inargnames)($vn, _model) || $(DynamicPPL.inmissings)($vn, _model)
                true
            else
                # Evaluate the LHS
                $expr === missing
            end
        end
    end
end

# failsafe: a literal is never an assumption
isassumption(expr) = :(false)

#################
# Main Compiler #
#################

"""
    @model(expr[, warn = true])

Macro to specify a probabilistic model.

If `warn` is `true`, a warning is displayed if internal variable names are used in the model
definition.

# Example

Model definition:

```julia
@model function model_generator(x = default_x, y)
    ...
end
```

To generate a `Model`, call `model_generator(x_value)`.
"""
macro model(expr, warn=true)
    esc(model(expr, warn))
end

function model(expr, warn)
    modelinfo = build_model_info(expr)

    # Generate main body
    modelinfo[:main_body] = generate_mainbody(modelinfo[:main_body], modelinfo[:args], warn)

    return build_output(modelinfo)
end

"""
    build_model_info(input_expr)

Builds the `model_info` dictionary from the model's expression.
"""
function build_model_info(input_expr)
    # Extract model name (:name), arguments (:args), (:kwargs) and definition (:body)
    modeldef = ExprTools.splitdef(input_expr)
    modeldef[:whereparams] = get(modeldef,:whereparams, [])
    modeldef[:args] = get(modeldef,:args, [])
    # Function body of the model is empty
    warn_empty(modeldef[:body])

    # Extracting the argument symbols from the model definition
    arg_syms = map(modeldef[:args]) do arg
        # @model demo(x)
        if arg isa Symbol
            arg
        elseif arg == Expr(:kw, IsEqual(x->true), IsEqual(x->true))
            expr, default = arg.args
            # @model demo(::Type{T}=Float64) where {T}
            if !isnothing(get_type(expr))
                # TODO support t::Type{T}
                !Meta.isexpr(expr, :(::), 2) || throw(ArgumentError("The syntax `t::Type{T}` is currently unsupported please use `::Type{T}` instead."))
                get_type(expr)
            # @model demo(x::Int = 1) or demo(x::Vector{Int} = [1,2,3])
            elseif !isnothing(get_symbol(expr))
                get_symbol(expr)
            # @model demo(x = 1)
            elseif expr isa Symbol
                expr
            else
                throw(ArgumentError("Unsupported default argument $arg to the `@model` macro."))
            end
        else
            throw(ArgumentError("Unsupported argument $arg to the `@model` macro."))
        end
    end
    if length(arg_syms) == 0
        args_nt = :(NamedTuple())
    else
        nt_type = Expr(:curly, :NamedTuple,
            Expr(:tuple, QuoteNode.(arg_syms)...),
            Expr(:curly, :Tuple, [:(Core.Typeof($x)) for x in arg_syms]...)
        )
        args_nt = Expr(:call, :($namedtuple), nt_type, Expr(:tuple, arg_syms...))
    end
    args = map(get(modeldef, :args, [])) do arg
        if (arg isa Symbol)
            arg
        elseif !isnothing(get_type(arg.args[1]))
            Texpr, Tval = arg.args
            T = get_type(Texpr)
            S = nothing
            for x in modeldef[:whereparams]
                if x == T
                    S = :Any
                elseif x == Expr(:<:, T, IsEqual(x->true))
                    S = x.args[2]
                end
            end
            !isnothing(S) || throw(ArgumentError("Please make sure type parameters are properly used. Every `Type{T}` argument need to have `T` in the a `where` clause"))
            Expr(:kw, :($T::Type{<:$S}), Tval)
        else
            arg
        end
    end
    args_nt = to_namedtuple_expr(arg_syms)

    default_syms = []
    default_vals = []
    foreach(get(modeldef, :args, [])) do arg
        if arg == Expr(:kw, IsEqual(x->true), IsEqual(x->true))
            var, val = arg.args
            sym = if var isa Symbol
                var
            elseif var == Expr(:(::), IsEqual(issymbol), IsEqual(x->true))
                var.args[1]
            elseif !isnothing(get_type(var))
                get_type(var)
            else
                error("Unsupported keyword argument given $arg")
            end
            push!(default_syms, sym)
            push!(default_vals, val)
        end
    end
    defaults_nt = to_namedtuple_expr(default_syms, default_vals)

    model_info = Dict(
        :name => modeldef[:name],
        :main_body => modeldef[:body],
        :arg_syms => arg_syms,
        :args_nt => args_nt,
        :defaults_nt => defaults_nt,
        :args => args,
        :whereparams => modeldef[:whereparams]
    )

    return model_info
end

"""
    generate_mainbody(expr, args, warn)

Generate the body of the main evaluation function from expression `expr` and arguments
`args`.

If `warn` is true, a warning is displayed if internal variables are used in the model
definition.
"""
generate_mainbody(expr, args, warn) = generate_mainbody!(Symbol[], expr, args, warn)

generate_mainbody!(found, x, args, warn) = x
function generate_mainbody!(found, sym::Symbol, args, warn)
    if warn && sym in INTERNALNAMES && sym ∉ found
        @warn "you are using the internal variable `$(sym)`"
        push!(found, sym)
    end
    return sym
end
function generate_mainbody!(found, expr::Expr, args, warn)
    # Do not touch interpolated expressions
    expr.head === :$ && return expr.args[1]

    # Apply the `@.` macro first.
    if Meta.isexpr(expr, :macrocall) && length(expr.args) > 1 &&
        expr.args[1] === Symbol("@__dot__")
        return generate_mainbody!(found, Base.Broadcast.__dot__(expr.args[end]), args, warn)
    end

    # Modify dotted tilde operators.
    args_dottilde = getargs_dottilde(expr)
    if args_dottilde !== nothing
        L, R = args_dottilde
        return generate_dot_tilde(generate_mainbody!(found, L, args, warn),
                                  generate_mainbody!(found, R, args, warn),
                                  args) |> Base.remove_linenums!
    end

    # Modify tilde operators.
    args_tilde = getargs_tilde(expr)
    if args_tilde !== nothing
        L, R = args_tilde
        return generate_tilde(generate_mainbody!(found, L, args, warn),
                              generate_mainbody!(found, R, args, warn),
                              args) |> Base.remove_linenums!
    end

    return Expr(expr.head, map(x -> generate_mainbody!(found, x, args, warn), expr.args)...)
end

"""
    generate_tilde(left, right, args)

Generate an `observe` expression for data variables and `assume` expression for parameter
variables.
"""
function generate_tilde(left, right, args)
    @gensym tmpright
    top = [:($tmpright = $right),
           :($tmpright isa Union{$Distribution,AbstractVector{<:$Distribution}}
             || throw(ArgumentError($DISTMSG)))]

    if left isa Symbol || left isa Expr
        @gensym out vn inds
        push!(top, :($vn = $(varname(left))), :($inds = $(vinds(left))))

        # It can only be an observation if the LHS is an argument of the model
        if vsym(left) in args
            @gensym isassumption
            return quote
                $(top...)
                $isassumption = $(DynamicPPL.isassumption(left))
                if $isassumption
                    $left = $(DynamicPPL.tilde_assume)(
                        _context, _sampler, $tmpright, $vn, $inds, _varinfo)
                else
                    $(DynamicPPL.tilde_observe)(
                        _context, _sampler, $tmpright, $left, $vn, $inds, _varinfo)
                end
            end
        end

        return quote
            $(top...)
            $left = $(DynamicPPL.tilde_assume)(_context, _sampler, $tmpright, $vn, $inds,
                                               _varinfo)
        end
    end

    # If the LHS is a literal, it is always an observation
    return quote
        $(top...)
        $(DynamicPPL.tilde_observe)(_context, _sampler, $tmpright, $left, _varinfo)
    end
end

"""
    generate_dot_tilde(left, right, args)

Generate the expression that replaces `left .~ right` in the model body.
"""
function generate_dot_tilde(left, right, args)
    @gensym tmpright
    top = [:($tmpright = $right),
           :($tmpright isa Union{$Distribution,AbstractVector{<:$Distribution}}
             || throw(ArgumentError($DISTMSG)))]

    if left isa Symbol || left isa Expr
        @gensym out vn inds
        push!(top, :($vn = $(varname(left))), :($inds = $(vinds(left))))

        # It can only be an observation if the LHS is an argument of the model
        if vsym(left) in args
            @gensym isassumption
            return quote
                $(top...)
                $isassumption = $(DynamicPPL.isassumption(left))
                if $isassumption
                    $left .= $(DynamicPPL.dot_tilde_assume)(
                        _context, _sampler, $tmpright, $left, $vn, $inds, _varinfo)
                else
                    $(DynamicPPL.dot_tilde_observe)(
                        _context, _sampler, $tmpright, $left, $vn, $inds, _varinfo)
                end
            end
        end

        return quote
            $(top...)
            $left .= $(DynamicPPL.dot_tilde_assume)(
                _context, _sampler, $tmpright, $left, $vn, $inds, _varinfo)
        end
    end

    # If the LHS is a literal, it is always an observation
    return quote
        $(top...)
        $(DynamicPPL.dot_tilde_observe)(_context, _sampler, $tmpright, $left, _varinfo)
    end
end

const FloatOrArrayType = Type{<:Union{AbstractFloat, AbstractArray}}
hasmissing(T::Type{<:AbstractArray{TA}}) where {TA <: AbstractArray} = hasmissing(TA)
hasmissing(T::Type{<:AbstractArray{>:Missing}}) = true
hasmissing(T::Type) = false

"""
    build_output(model_info)

Builds the output expression.
"""
function build_output(model_info)
    # Arguments with default values
    args = model_info[:args]
    # Argument symbols without default values
    arg_syms = model_info[:arg_syms]
    # Arguments namedtuple
    args_nt = model_info[:args_nt]
    # Default values of the arguments
    # Arguments namedtuple
    defaults_nt = model_info[:defaults_nt]
    # Where parameters
    whereparams = model_info[:whereparams]
    # Model generator name
    model_gen = model_info[:name]
    # Main body of the model
    main_body = model_info[:main_body]

    unwrap_data_expr = Expr(:block)
    for var in arg_syms
        push!(unwrap_data_expr.args,
              :($var = $(DynamicPPL.matchingvalue)(_sampler, _varinfo, _model.args.$var)))
    end

    @gensym(evaluator, generator)
    generator_kw_form = isempty(args) ? () : (:($generator(;$(args...)) = $generator($(arg_syms...))),)
    model_gen_constructor = :($(DynamicPPL.ModelGen){$(Tuple(arg_syms))}($generator, $defaults_nt))

    return quote
        function $evaluator(
            _model::$(DynamicPPL.Model),
            _varinfo::$(DynamicPPL.AbstractVarInfo),
            _sampler::$(DynamicPPL.AbstractSampler),
            _context::$(DynamicPPL.AbstractContext),
        )
            $unwrap_data_expr
            $main_body
        end

        $generator($(args...)) = $(DynamicPPL.Model)($evaluator, $args_nt, $model_gen_constructor)
        $(generator_kw_form...)

        $(Base).@__doc__ $model_gen = $model_gen_constructor
    end
end

function warn_empty(body)
    if all(l -> isa(l, LineNumberNode), body.args)
        @warn("Model definition seems empty, still continue.")
    end
    return
end

"""
    matchingvalue(sampler, vi, value)

Convert the `value` to the correct type for the `sampler` and the `vi` object.
"""
function matchingvalue(sampler, vi, value)
    T = typeof(value)
    if hasmissing(T)
        return get_matching_type(sampler, vi, T)(value)
    else
        return value
    end
end
matchingvalue(sampler, vi, value::FloatOrArrayType) = get_matching_type(sampler, vi, value)

"""
    get_matching_type(spl, vi, ::Type{T}) where {T}
Get the specialized version of type `T` for sampler `spl`. For example,
if `T === Float64` and `spl::Hamiltonian`, the matching type is `eltype(vi[spl])`.
"""
function get_matching_type end


text:
module netgengen

abstract type CSGObject end
abstract type LineObject end
abstract type CurveObject <: LineObject end
abstract type CSGCompositeObject <: CSGObject end
abstract type CSGPrimitive <: CSGObject end

export torus, plane, CSGObject, brick, csgunion, csgstring, declare, not, cylinder, sphere, tlo, intersection, curve2d, LineObject, revolution, CurveObject

function ComposeModifiers(flags...)
  output = ""
  for flag in flags
    if(typeof(flag[2])==Bool)
      output = string(output, " -", string(flag[1]))
    else
      output = string(output, " -", string(flag[1]), "=", string(flag[2]))
    end
  end
  return output
end

"""
`function tlo(E::CSGObject, buf=STDOUT; flags...)`

Prints out netgen string of `E` to buffer `buf` and appends `flags` to the end of the string

"""
function tlo(E::CSGObject, buf=STDOUT; flags...)
  if !(E.declared)
    declare(E,buf)
  end
  print(buf, "tlo ", csgstring(E), ComposeModifiers(flags...), ";\n")
end

"""
`type not`

`not(E)` constructs complement of `E`

`not(name, E; modifiers...)` constructs complement of `E`, assigns it with String `name` and `modifiers...`
"""
type not <: CSGCompositeObject
  name::String
  E::Array{CSGObject,1}
  expr::String
  modifiers::String
  declared::Bool
  not(E) = (x = new();
    x.E = [E];
    x.expr = string("not ", csgstring(E));
    x.name = "";
    x.modifiers = "";
    x.declared = false;
    return x)
  not(name, E; modifiers...) = (x = not(E); x.name=name; x.modifiers = ComposeModifiers(modifiers...); return x)
end

"""
`type torus`

`torus(p,q,m,n)` constructs torus from length 3 vectors `p`,`q` and scalars `m` and `n` where

* `p` is the center of the torus
* `q` is the normal vector to the torus
* `m` is the distance from center to the center of the torus ring
* `n` is the radius of the torus ring

`torus(name,p,q,m,n;modifiers...)` same as above but assigns it with `name` and `modifiers...`
"""
type torus <: CSGObject
  name::String
  p::Array{AbstractFloat, 1}
  q::Array{AbstractFloat, 1}
  m::AbstractFloat
  n::AbstractFloat
  expr::String
  modifiers::String
  declared::Bool
  torus(p,q,m,n) = (
    x = new();
    x.p = p;
    x.q = q;
    x.m = m;
    x.n = n;
    x.expr = string("torus ($(x.p[1]), $(x.p[2]), $(x.p[3]); $(x.q[1]), $(x.q[2]), $(x.q[3]); $(x.m); $(x.n))");
    x.modifiers = "";
    x.declared = false;
    x.name = "";
    return x)
  torus(name, p, q, m,n; modifiers...) = (
    x = (torus(p,q,m,n));
    x.name = name;
    x.modifiers= ComposeModifiers(modifiers...);
    return x)
end

"""
`type cylinder`

`cylinder(p,q,r)` constructs cylinder where

* `p` and `q` are length 3 vectors for endpoints of the cylinder and `r` is the radius of the cylinder

`cylinder(name,p,q,r;modifiers...)`

* same as above but adds `name` and `modifiers...`
"""
type cylinder <: CSGObject
  name::String
  p::Array{AbstractFloat, 1}
  q::Array{AbstractFloat, 1}
  r::AbstractFloat
  expr::String
  modifiers::String
  declared::Bool
  cylinder(p,q,r) = (
    x = new();
    x.p = p;
    x.q = q;
    x.r = r;
    x.expr = string("cylinder ($(x.p[1]), $(x.p[2]), $(x.p[3]); $(x.q[1]), $(x.q[2]), $(x.q[3]); $(x.r))");
    x.modifiers = "";
    x.declared = false;
    x.name = "";
    return x)
  cylinder(name, p, q, r; modifiers...) = (
    x = cylinder(p,q,r);
    x.name = name;
    x.modifiers= ComposeModifiers(modifiers...);
    return x)
end

"""
`type plane`

`plane(p,n)` constructs half space where

* `p` is a point on its boundary and
* `n` outside normal vector to the plane
"""
type plane <: CSGObject
  name::String
  p::Array{AbstractFloat, 1}
  n::Array{AbstractFloat, 1}
  expr::String
  modifiers::String
  declared::Bool
  plane(p,n;modifiers...) = (
    x = new();
    x.p = p;
    x.n = n;
    x.expr = string("plane ($(x.p[1]), $(x.p[2]),$(x.p[3]); $(x.n[1]), $(x.n[2]), $(x.n[3]))");
    x.modifiers = ComposeModifiers(modifiers...);
    x.declared = false;
    x.name = "";
    return x)
  plane(name::String, p, n; modifiers...) = (
    x = plane(p,n);
    x.name=name;
    x.modifiers = ComposeModifiers(modifiers...);
    x.declared = false;
    return x)
end

"""
`type brick`

`brick(p,q)` constructs 3D interval where with coordinate minima in vector `p` and maxima in `q`
"""
type brick <: CSGObject
  name::String
  p::Array{AbstractFloat, 1}
  q::Array{AbstractFloat, 1}
  expr::String
  modifiers::String
  declared::Bool
  brick(p,q) = (
    x = new();
    x.p = p;
    x.q = q;
    x.expr = string("orthobrick ($(x.p[1]), $(x.p[2]), $(x.p[3]); $(x.q[1]), $(x.q[2]), $(x.q[3]))");
    x.modifiers = "";
    x.declared = false;
    x.name = "";
    return x)
  brick(name::String, p, q; modifiers...) = (
    x = brick(p,q);
    x.name = name;
    x.modifiers = ComposeModifiers(modifiers...);
    x.declared = false;
    return x)
end

type curve3d <: CurveObject
  name::String
  np::Integer
  ns::Integer
  p::Array{Float64, 2}
  t::Array{Array{Integer,1},1}
  expr::String
  modifiers::String
  declared::Bool
  curve3d(name, p,t) = (
    x = new();
    x.p = p;
    x.t = t;
    x.np = size(p,2);
    x.ns = length(t);
    x.name = name;
    evalexpr(x);
    return x;
    )
end

type curve2d <: CurveObject
  name::String
  np::Integer
  ns::Integer
  p::Array{Float64, 2}
  t::Array{Array{Integer,1},1}
  expr::String
  modifiers::String
  declared::Bool
  curve2d(name, p,t) = (
    x = new();
    x.p = p;
    x.t = t;
    x.np = size(p,2);
    x.ns = size(t,1);
    x.name = name;
    x.expr = evalexpr(x);
    return x;
    )
end

function evalexpr(x::curve2d)
  expr = string("curve2d $(x.name) = ($(x.np);");
  for k = 1:x.np
    expr = string(expr, "\n\t\t$(x.p[1,k]), $(x.p[2,k]);")
  end
  expr = string(expr, "\n\t\t$(x.ns);")
  for k = 1:length(x.t)-1
    if length(x.t[k]) == 2
      expr = string(expr, "\n\t\t2, $(x.t[k][1]), $(x.t[k][2]);")
    elseif length(x.t[k]) == 3
      expr = string(expr, "\n\t\t3, $(x.t[k][1]), $(x.t[k][2]), $(x.t[k][3]);")
    else
      error("Wrong size of segment $(x.t[k])")
    end
  end
  seg = x.t[end]
  if length(seg) == 2
    expr = string(expr, "\n\t\t2, $(seg[1]), $(seg[2])")
  elseif length(seg) == 3
    expr = string(expr, "\n\t\t3, $(seg[1]), $(seg[2]), $(seg[3])")
  else
    error("Wrong size of segment $(seg)")
  end
  expr = string(expr, ");\n")
  return expr
end

function evalexpr(x::curve3d)
  expr = string("curve3d $(x.name) = ($(x.np),");
  for k = 1:x.np
    expr = string(expr, "\n\t\t$(x.p[1,k]), $(x.p[2,k]), $(x.p[3,k]);")
  end
  expr = string(expr, "\n\t\t$(x.ns);")
  for k = 1:length(x.t)-1
    if length(x.t[k]) == 2
      expr = string(expr, "\n\t\t2, $(x.t[k][1]), $(x.t[k][2]);")
    elseif length(x.t[k]) == 3
      expr = string(expr, "\n\t\t3, $(x.t[k][1]), $(x.t[k][2]), $(x.t[k][3]);")
    else
      error("Wrong size of segment $(x.t[k])")
    end
  end
  seg = x.t[end]
  if length(seg) == 2
    expr = string(expr, "\n\t\t2, $(seg[1]), $(seg[2])")
  elseif length(seg) == 3
    expr = string(expr, "\n\t\t3, $(seg[1]), $(seg[2]), $(seg[3])")
  else
    error("Wrong size of segment $(seg)")
  end
  expr = string(expr, ");\n")
  return expr
end

function MakeExpression!{T}(object::T)
  error("Not done!")
end

"""
`type sphere`

`sphere(c,R)` construct sphere with center `c` and radius `R`
"""
type sphere <: CSGObject
  name::String
  c::Array{AbstractFloat, 1}
  R::AbstractFloat
  expr::String
  modifiers::String
  declared::Bool
  sphere(c,R) = (
    x = new();
    x.c = c;
    x.R = R;
    x.expr = string("sphere ($(x.c[1]), $(x.c[2]), $(x.c[3]); $(R))");
    x.modifiers = "";
    x.declared = false;
    x.name = "";
    return x)
  sphere(name, c, R; modifiers...) = (
    x = (sphere(c,R));
    x.name = name;
    x.modifiers = ComposeModifiers(modifiers...);
    x.declared = false;
    return x)
end

type revolution <: CSGObject
  name::String
  p1::Array{AbstractFloat,1}
  p2::Array{AbstractFloat,1}
  curvename::String
  expr::String
  modifiers::String
  declared::Bool
  revolution(p1,p2,curvename) = (
    x = new();
    x.p1 = p1;
    x.p2 = p2;
    x.expr = string("revolution ($(x.p1[1]), $(x.p1[2]), $(x.p1[3]); $(x.p2[1]), $(x.p2[2]), $(x.p2[3]); $(curvename))");
    x.modifiers = "";
    x.declared = false;
    x.name = "";
    return x)
  revolution(name, p1, p2, curvename; modifiers...) = (
    x = (revolution(p1, p2, curvename));
    x.name = name;
    x.modifiers = ComposeModifiers(modifiers...);
    x.declared = false;
    return x)
end

"""
`type csgunion`

`csgunion(E)` constructs union of objects in array `E`

`csgunion(name,E,modifiers...)` constructs union of objects in array `E` and assigns it with `name` and `modifiers...`
"""
type csgunion <: CSGCompositeObject
  name::String
  E::Array{CSGObject,1}
  expr::String
  modifiers::String
  declared::Bool
  csgunion(E) = (
    x = new();
    x.E = E;
    x.expr = evalexpr(x);
    x.name = "";
    x.modifiers = "";
    x.declared = false;
    return x)
  csgunion(name, E; modifiers...) = (
    x = csgunion(E);
    x.name = name;
    x.modifiers = ComposeModifiers(modifiers...);
    x.declared = false;
    return x)
end

"""
`type intersection`

`intersection(E)` constructs intersection of objects in array `E`
`intersection(name,E,modifiers...)` constructs intersection of objects in array `E` and assigns it with `name` and `modifiers...`
"""
type intersection <: CSGCompositeObject
  name::String
  E::Array{CSGObject,1}
  expr::String
  modifiers::String
  declared::Bool
  intersection(E) = (
    x = new();
    x.E = E;
    x.expr = evalexpr(x);
    x.name = "";
    x.modifiers = "";
    x.declared = false;
    return x)
  intersection(name, E; modifiers...) = (
    x = intersection(E);
    x.name = name;
    x.modifiers = ComposeModifiers(modifiers...);
    x.declared = false;
    return x)
end

function csgstring(A::CSGObject)
  return A.name == "" ? A.expr : A.name;
end

function declare(A::CSGCompositeObject, buf=STDOUT)
  for sub_obj in A.E
    if(!sub_obj.declared)
      declare(sub_obj, buf)
    end
  end
  A.declared = true
  if A.name == ""
    return true
  end
  println(buf, "solid ", A.name, " = ", A.expr, A.modifiers, ";")
  return true
end

function declare(A::CSGPrimitive, buf=STDOUT)
  if(!(A.declared))
    A.declared = true
  else
    return true
  end
  if A.name == ""
    return true
  end
  println(buf, "solid ", A.name, " = ", A.expr, A.modifiers, ";")
  return true
end

function declare(A::CSGObject, buf=STDOUT)
  if(!(A.declared))
    A.declared = true
  else
    return true
  end
  if A.name == ""
    return true
  end
  println(buf, "solid ", A.name, " = ", A.expr, A.modifiers, ";")
  return true
end

function declare(A::CurveObject, buf=STDOUT)
  println("warning: declare(A::CurveObject) might not do what wanted")
  return A.expr;
end

function evalexpr(A::csgunion)
  S = ""
  for m = 1:length(A.E)-1
    S = string(S, csgstring(A.E[m]), " or\n\t")
  end
  S = string(S, csgstring(A.E[length(A.E)]))
end

function evalexpr(A::intersection)
  S = ""
  for m = 1:length(A.E)-1
    S = string(S, csgstring(A.E[m]), " and\n\t")
  end
  S = string(S, csgstring(A.E[length(A.E)]))
end

# Utility functions
function finitecylinder(name, p,q,r; modifiers...)
  return intersection(name, [
      cylinder(p,q,r),
      plane(p,p-q; bc=a),
      plane(q,q-p; bc=b)]; modifiers...)
end

function finitecylinder(name, p,q,r, a, b; modifiers...)
  plane_a = plane(string(name, "_", a), p,p-q; bc=a)
  plane_b = plane(string(name, "_", b), q,q-p; bc=b)
  return intersection(name, [
      cylinder(p,q,r),
      plane_a, plane_b]; modifiers...)
end

end

ext: jl


text:
# See https://en.wikipedia.org/wiki/List_of_fractals_by_Hausdorff_dimension

const Fractal{D} = Vector{Segment{D}}

function unit_circle(npoints::Int)::Fractal{2}
    @assert npoints > 2
    θ = linspace(0, 2*pi, npoints)
    x = cos.(θ)
    y = sin.(θ)
    [Segment(Point(x[i], y[i]), Point(x[i+1], y[i+1])) for i in 1:npoints-1]
end

#"""Step a segment into
#   /\
#--/  \--
#"""
function koch_line(line::Segment{2})::Fractal{2}
    p0, p1, p2, p3 = line.O, line(1./3), line(2./3), line.B

    v = 0.5*(p1 + p2)
    v += Point(-line.t[2], line.t[1])*sqrt(3)/2*line.length/3

    [Segment(p0, p1), Segment(p1, v), Segment(v, p2), Segment(p2, p3)]
end

"""Grow by stepiing from straigt line. FD = log(4)/log(3)"""
function koch_line(depth::Int)::Fractal{2}
    @assert depth > 0
    lines = [Segment(Point(0, 0), Point(1, 0))]
    while depth > 0
        lines = vcat(map(koch_line, lines)...)
        depth -= 1
    end
    lines
end

"""Grow (outward) by stepping from quilateral triangle. FD = log(4)/log(3)"""
function koch_flake(depth::Int)::Fractal{2}
    @assert depth > 0
    lines = [Segment(Point(0, 0), Point(0.5, sqrt(3)/2)),
             Segment(Point(0.5, sqrt(3)/2), Point(1, 0)),
             Segment(Point(1, 0), Point(0, 0))]
    while depth > 0
        lines = vcat(map(koch_line, lines)...)
        depth -= 1
    end
    lines
end

"""Like Kock flake but grows inward. Triflake in wiki?. FD = log(4)/log(3)"""
function triflake(depth::Int)::Fractal{2}
    @assert depth > 0
    lines = [Segment(Point(0, 0), Point(1, 0)),
             Segment(Point(1, 0), Point(0.5, sqrt(3)/2)),
             Segment(Point(0.5, sqrt(3)/2), Point(0, 0))]
    while depth > 0
        lines = vcat(map(koch_line, lines)...)
        depth -= 1
    end
    lines
end

# --------------------------------------------------------------------

#"""                  |-|
#--------  becomes   -| |-
#"""
function koch_quadratic_1(line::Segment{2})
    v1, v2, v3, v4 = line.O, line(1./3), line(2./3), line.B
    shift = Point(-line.t[2], line.t[1])*line.length/3

    [Segment(v1, v2),
     Segment(v2, v2 + shift),
     Segment(v2 + shift, v3 + shift),
     Segment(v3 + shift, v3),
     Segment(v3, v4)]
end

"""Grow by stepping from (0, 0) - (1, 0). FD = log(5)/log(3)"""
function koch_quadratic_1(depth::Int)::Fractal{2}
    @assert depth > 0
    lines = [Segment(Point(0, 0), Point(1, 0))]
    while depth > 0
        lines = vcat(map(koch_quadratic_1, lines)...)
        depth -= 1
    end
    lines
end

# --------------------------------------------------------------------

#"""                  |-|
#--------  becomes   -| | |-
#                       |-|
#"""
function koch_quadratic_2(line::Segment{2})
    v1, v2, v3, v4, v5 = line.O, line(0.25), line(0.5), line(0.75), line.B
    shift = Point(-line.t[2], line.t[1])*line.length/4

    [Segment(v1, v2),
     Segment(v2, v2 + shift),
     Segment(v2 + shift, v3 + shift),
     Segment(v3 + shift, v3),
     Segment(v3, v3 - shift),
     Segment(v3 - shift, v4 - shift),
     Segment(v4 - shift, v4),
     Segment(v4, v5)]
end

"""Grow by stepping from (0, 0) - (1, 0). FD = 1.5"""
function koch_quadratic_2(depth::Int)::Fractal{2}
    @assert depth > 0
    lines = [Segment(Point(0, 0), Point(1, 0))]
    while depth > 0
        lines = vcat(map(koch_quadratic_2, lines)...)
        depth -= 1
    end
    lines
end

# --------------------------------------------------------------------

#"""               _
#\  / becomes   | |
# \/            |_|
#"""
function dragon_curve(lines::NTuple{2, Segment{2}})
    line = first(lines)
    v1, v2 = line.O, line.B
    shift = Point(line.t[2], -line.t[1])*line.length/2
    m = 0.5*(v1 + v2) + shift

    l1 = (Segment(v1, m), Segment(m, v2))

    line = last(lines)
    v1, v2 = line.O, line.B
    shift = Point(-line.t[2], line.t[1])*line.length/2
    m = 0.5*(v1 + v2) + shift

    [l1, (Segment(v1, m), Segment(m, v2))]
end

"""Grow by stepping from (0, 0) - (1, -1) - (2, 0). FD = 1.526"""
function dragon_curve(depth::Int)::Fractal{2}
    @assert depth > 0
    pairs = [(Segment(Point(0, 0), Point(1, -1)),
              Segment(Point(1, -1), Point(2, 0)))]
    while depth > 0
        pairs = vcat(map(dragon_curve, pairs)...)
        depth -= 1
    end
    # Saw them together
    [l for pair in pairs for l in pair]
end

max_stars_repo_name: MiroK/fractional-dim


text:
# This file is part of the TaylorIntegration.jl package; MIT licensed
@auto_hash_equals struct TaylorInterpolant{T,U,N}
    t0::T
    t::AbstractVector{T}
    x::AbstractArray{Taylor1{U},N}
    #Inner constructor
    function TaylorInterpolant{T,U,N}(
        t0::T,
        t::AbstractVector{T},
        x::AbstractArray{Taylor1{U},N}
    ) where {T<:Real, U<:Number, N}
        @assert size(x)[1] == length(t)-1
        @assert issorted(t) || issorted(t, rev=true)
        return new{T,U,N}(t0, t, x)
    end
end

#outer constructor
function TaylorInterpolant(t0::T, t::AbstractVector{T},
        x::AbstractArray{Taylor1{U},N}) where {T<:Real, U<:Number, N}
    return TaylorInterpolant{T,U,N}(t0, t, x)
end

# return time vector index corresponding to interpolation range
function getinterpindex(tinterp::TaylorInterpolant{T,U,N}, t::V) where {T<:Real, U<:Number, V<:Number, N}
    t00 = constant_term(constant_term(t))
    tmin, tmax = minmax(tinterp.t[end], tinterp.t[1])
    Δt = t-tinterp.t0
    Δt00 = t00-tinterp.t0
    @assert tmin ≤ Δt00 ≤ tmax "Evaluation time outside range of interpolation"
    if Δt00 == tinterp.t[end] # compute solution at final time from last step expansion
        ind = lastindex(tinterp.t)-1
    elseif issorted(tinterp.t) # forward integration
        ind = searchsortedlast(tinterp.t, Δt00)
    elseif issorted(tinterp.t, rev=true) # backward integration
        ind = searchsortedlast(tinterp.t, Δt00, rev=true)
    end
    return ind, Δt
end

# function-like (callability) methods
function (tinterp::TaylorInterpolant{T,U,1})(t::V) where {T<:Real, U<:Number, V<:Number}
    ind, Δt = getinterpindex(tinterp, t)
    δt = Δt-tinterp.t[ind]
    return tinterp.x[ind](δt)
end

function (tinterp::TaylorInterpolant{T,U,2})(t::V) where {T<:Real, U<:Number, V<:Number}
    ind, Δt = getinterpindex(tinterp, t)
    δt = Δt-tinterp.t[ind]
    return tinterp.x[ind,:](δt)
end

# reverse independent variable function
function reverse(tinterp::TaylorInterpolant{T,U,N}) where {T<:Real, U<:Number, N}
    tinterp_rev_t0 = tinterp.t[end]
    tinterp_rev_t = tinterp.t[end:-1:1] .- tinterp_rev_t0
    tinterp_rev_x = vcat(tinterp(tinterp.t[end]+Taylor1(tinterp.x[1].order))', tinterp.x[end:-1:2,:])
    return TaylorInterpolant(tinterp_rev_t0, tinterp_rev_t, tinterp_rev_x)
end


text:
#=
   Interpolation module

    Interpolate two and three dimensional arrays.
=#

module Interpolate

export interpolateCalculateSpace2D!, interpolateCalculateSpace3D!, interpolateArray2D!, interpolateArray3D!

#------------------------------------------------------------------
#
# interpolateCalculateSpace2D -calculating the lower and upper
# bounds of an interpolated 2D array with regard to interpFactor.
#
# Parameters:
# - odims, 4 elements
# - interpFactor, interpolation factor
# - fdims, 4 elements
#
#------------------------------------------------------------------
function interpolateCalculateSpace2D!(odims, interpfactor, fdims)

  fdims[1] = odims[1]
  fdims[3] = odims[3]
  fdims[2] = ((odims[2] - odims[1] + 1) + (odims[2] - odims[1]) * interpfactor) + (odims[1] - 1)
  fdims[4] = ((odims[4] - odims[3] + 1) + (odims[4] - odims[3]) * interpfactor) + (odims[3] - 1)
end # function interpolatecalculatespace2d

#------------------------------------------------------------------
#
# interpolateCalculateSpace3D
#
#
# Subroutine for calculating the lower and upper bounds of an
# interpolated 3D array with regard to interpFactor.  The r is
# expected to this routine, allocate the interpolated array,
#  function: interpolateArray3D to fill the
# interpolated array.
#
#------------------------------------------------------------------
function interpolateCalculateSpace3D!(odims, interpfactor, fdims)

  fdims[1] = odims[1]
  fdims[3] = odims[3]
  fdims[5] = odims[5]
  fdims[2] = ((odims[2] - odims[1] + 1) + (odims[2] - odims[1]) * interpfactor) + (odims[1] - 1)
  fdims[4] = ((odims[4] - odims[3] + 1) + (odims[4] - odims[3]) * interpfactor) + (odims[3] - 1)
  fdims[6] = odims[6] # The 3rd dimension is unchanged.
end # function interpolatecalculatespace3d

#------------------------------------------------------------------
#
# interpolateArray2D
#
# Bilinear interpolation of f with respect to oa.
#
# f"s size has been calculated in function interpolateCalculateSpace2D
# f has been allocated in the tree above this function
#
# Subroutine Parameters:
# oa - original array, 2D
# odims - low / high subscripts for oa
# f - interpolated array, 2D
# fdims - low / high subscripts for f
# interpFactor - interpolation factor

# Here are some illustrations of interpFactor.
# We are assuming that the grid points are equally spaced.
#
# An interpolation factor of 1 means 1 new interpolated element.
# So for example, a 3x3 matrix with interpFactor = 1 becomes a 5x5 matrix
#
# x"s are data from oa, the o"s are points to be interpolated
#
# [ x x x ]     == > [x o x o x]
# [ x x x ]          [o o o o o]
# [ x x x ]          [x o x o x]
#                    [o o o o o]
#                    [x o x o x]
#
# And an interpFactor = 2, means a 3x3 matrix becomes a 7x7 matrix.
# [ x x x ]     == > [x o o x o o x]
# [ x x x ]          [o o o o o o o]
# [ x x x ]          [o o o o o o o]
#                    [x o o x o o x]
#                    [o o o o o o o]
#                    [o o o o o o o]
#                    [x o o x o o x]
#
#------------------------------------------------------------------
function interpolateArray2D!(oa, odims, f, fdims, interpfactor::Int64)

  # Unlock the subscript bounds of the originalArray.
  ol1 = odims[1]
  ou1 = odims[2]
  ol2 = odims[3]
  ou2 = odims[4]

  # Unlock the subscript bounds of the interpolatedArray.
  fl1 = fdims[1]
  fu1 = fdims[2]
  fl2 = fdims[3]
  fu2 = fdims[4]

  # Intersperse the oa points into f.
  for j = ol2: ou2
    for i = ol1: ou1
      f[ol1 + (i - ol1) * (interpfactor + 1), ol2 + (j - ol2) * (interpfactor + 1)] = oa[i, j]
    end
  end

  # Loop over all the squares.
  for j = fl2: interpfactor + 1 : fu2 - 1
    for i = fl1: interpfactor + 1 : fu1 - 1
      # Find the indices of the corner points
      x1 = i
      x2 = i + interpfactor + 1
      y1 = j
      y2 = j + interpfactor + 1
      # Find the value of the corner points
      fq11 = f[x1, y1]
      fq12 = f[x1, y2]
      fq21 = f[x2, y1]
      fq22 = f[x2, y2]
      # Loop over all the interpolated points
      for iy = j: j + interpfactor + 1
        for ix = i: i + interpfactor + 1
          # Skip the corner points
          if ((ix == x1) && (iy == y1)) ||
             ((ix == x1) && (iy == y2)) ||
             ((ix == x2) && (iy == y1)) ||
             ((ix == x2) && (iy == y2))
            continue
          end
          # Get the weights for the x direction
          wx1 = real(x2 - ix) / real(x2 - x1)
          wx2 = real(ix - x1) / real(x2 - x1)
          # Get the weights for the y direction
          wy1 = real(y2 - iy) / real(y2 - y1)
          wy2 = real(iy - y1) / real(y2 - y1)
          # interpolate
          f[ix, iy] = wy1 * (wx1 * fq11 + wx2 * fq21) + wy2 * (wx1 * fq12 + wx2 * fq22)
        end # x loop
      end   # y loop
    end     # i loop
  end       # j loop
end # function interpolatearray2d

#------------------------------------------------------------------
#
# interpolateArray3D
#
# Bilinear interpolation of f with respect to oa.  The interpolation
# ignores the 3rd dimension.
#
#------------------------------------------------------------------
function interpolateArray3D!(oa, odims, f, fdims, interpfactor::Int64)
  for i = odims[5]: odims[6]
    interpolateArray2D!(view(oa, :,:,i), odims, view(f, :,:,i), fdims, interpfactor)
  end
end # function interpolatearray3d

end # module interpolate


text:
module ScoringTest

import PlayingCards: Card
using  Test

import Cribbage
import Cribbage.Scoring

@testset "Scoring Hands" begin

    hand = parse(Array{Card,1}, "8♣ 7♥ 7♠ 6♦")
    starter = parse(Card, "2♠")

    @test Cribbage.Scoring.scorefifteen(hand, starter) == 8
    @test Cribbage.Scoring.scorepair(hand, starter) == 2
    @test Cribbage.Scoring.scorerun(hand, starter) == 6
    @test Cribbage.Scoring.scoreflush(hand, starter) == 0
    @test Cribbage.Scoring.scorenobs(hand, starter) == 0

    @test Cribbage.Scoring.scorehand(hand, starter) == 16

    # Perfect Hand
    hand = parse(Array{Card,1}, "5♣ 5♥ J♠ 5♦")
    starter = parse(Card, "5♠")

    @test Cribbage.Scoring.scorefifteen(hand, starter) == 16
    @test Cribbage.Scoring.scorepair(hand, starter) == 12
    @test Cribbage.Scoring.scorerun(hand, starter) == 0
    @test Cribbage.Scoring.scoreflush(hand, starter) == 0
    @test Cribbage.Scoring.scorenobs(hand, starter) == 1

    @test Cribbage.Scoring.scorehand(hand, starter) == 29
end

@testset "Scoring Play" begin

    # Scoring the play
    prior = parse(Array{Card,1}, "5♣")
    card = parse(Card, "10♠")
    @test Cribbage.Scoring.scoreplay(prior, card) == 2
    push!(prior, card)

    card = parse(Card, "10♦")
    @test Cribbage.Scoring.scoreplay(prior, card) == 2
    push!(prior, card)

    card = parse(Card, "6♠")
    @test Cribbage.Scoring.scoreplay(prior, card) == 2
    push!(prior, card)

    prior = parse(Array{Card,1}, "9♣ 6♦ 8♥")
    card = parse(Card, "7♣")
    @test Cribbage.Scoring.scoreplay(prior, card) == 4

    # Score the starter card
    @test Cribbage.Scoring.scoreheels(parse(Card, "J♠")) == 2
    @test Cribbage.Scoring.scoreheels(parse(Card, "A♠")) == 0

end;

end  # module ScoringTest


text:
module OptimTests

using Photon, Test
using Photon.Layers
using Photon.Optimisers: Adam, SGD, Momentum
using Photon.Callbacks: SilentMeter
using Photon.Losses

function get_model()
    Sequential(
        Conv2D(16, (3,3), :relu),
        Conv2D(32, (5,5), :relu),
        MaxPool2D(),
        Dense(64, :relu),
        Dense(10)
    )
end

function testOptim(opt)
    model = get_model()

    X = [randn(50,50,3,4) for _=1:5]
    Y = [randn(10,4) for _=1:5]

    workout = Workout(model, MSELoss(), opt)
    train!(workout, zip(X,Y), epochs=2, cb=SilentMeter())
    @test hasmetric(workout, :loss)
end

@testset "Optim" begin
    testOptim(Adam())
    testOptim(SGD())
    testOptim(Momentum())
end

end


text:
# Read input
readarray(x) = "[$x]" |> Meta.parse |> eval
firstline, rest = Base.Iterators.peel(readlines("input.txt"))
nums = readarray(firstline)
boards = map(readarray ∘ x -> join(x, '\n'), Base.Iterators.partition(rest, 6))

# Bingo functions
bingo_rows = vcat(
    [[CartesianIndex(i, j) for i in 1:5] for j in 1:5], # column indices
    [[CartesianIndex(i, j) for j in 1:5] for i in 1:5], # row indices
)
mark_board(board, i) = reduce(.|, nums[j] .∈ board for j in 1:i)
check_win(marks) = map(all ∘ x -> marks[x], bingo_rows) |> any
win_time(board) = map(check_win ∘ i -> mark_board(board, i), eachindex(nums)) |> findfirst
score_board(board) = (i = win_time(board); nums[i] * sum(board[.!mark_board(board, i)]))
win_times = map(win_time, boards)

# Part 1
first_time, i = findmin(win_times)
@info "Score of first winning board = $(score_board(boards[i]))" first_time

# Part 2
last_time, i = findmax(win_times)
@info "Score of last winning board = $(score_board(boards[i]))" last_time


text:
"""
    KaleidoLens <: Lens

Internal abstract type for Kaleido.jl.
"""
abstract type KaleidoLens <: Lens end

_getfields(obj) = map(n -> getfield(obj, n), fieldnames(typeof(obj))) :: Tuple

struct Prefixed
    prefix::String
    name::Symbol
end

prefixof(f) = join(fullname(parentmodule(f)), '.')
prefixof(f::Prefixed) = f.prefix
Base.nameof(f::Prefixed) = f.name

function print_apply(io, f, args)
    if !get(io, :limit, false)
        # Don't show full name in REPL etc.:
        print(io, prefixof(f), '.')
    end
    print(io, nameof(f))
    if length(args) == 1
        print(io, '(')
        show(io, args[1])
        print(io, ')')
    else
        show(io, args)
    end
    return
end

_default_show(io, obj) = print_apply(io, typeof(obj), _getfields(obj))

Base.show(io::IO, lens::KaleidoLens) = _default_show(io, lens)
Base.show(io::IO, ::MIME"text/plain", lens::KaleidoLens) = prettylens(io, lens)

# TODO: Remove custom `_constructor_of` definitions.  These are
# different from ConstructionBase spec.
@generated _constructor_of(::Type{T}) where T =
    getfield(parentmodule(T), nameof(T))
_constructor_of(::Type{<:NamedTuple{names}}) where names = NamedTuple{names}
_constructor_of(::Type{<:Tuple}) = Tuple

"""
    prefer_singleton_callable(f)

Convert `f` to an callable singleton object if possible.  Useful if `f`
is a `Type`.

# Examples
```jldoctest
julia> using Kaleido: prefer_singleton_callable

julia> sizeof((Int,))
8

julia> sizeof((prefer_singleton_callable(Int),))
0

julia> prefer_singleton_callable(Int)(1.0)
1
```
"""
prefer_singleton_callable(::Type{T}) where T = SingletonCallable{T}()
prefer_singleton_callable(f) = f

struct SingletonCallable{T} end
(::SingletonCallable{T})(x) where T = T(x)

@nospecialize
_singleton_callable(::SingletonCallable{T}) where T = T
_singleton_callable(f) = f
@specialize

## Specialized foldl

_tail(t) = Base.tail(t)
_tail(t::NamedTuple{names}) where names = NamedTuple{Base.tail(names)}(t)

struct _Zip{T1, T2}
    it1::T1
    it2::T2
end

const EmptyTuple = Union{Tuple{}, NamedTuple{(),Tuple{}}}
const EmptyItr = Union{EmptyTuple, _Zip{<:EmptyTuple, <:EmptyTuple}}
const AnyItr = Union{Tuple, NamedTuple, _Zip}
const _zip = _Zip

@inline _tail(it::_Zip) = _Zip(_tail(it.it1), _tail(it.it2))
@inline Base.getindex(it::_Zip, i) = (it.it1[i], it.it2[i])

@inline _mapfoldl(::Any, ::Any, ::EmptyItr, init) = init
@inline _mapfoldl(f, op, xs::AnyItr, init) =
    _mapfoldl(f, op, _tail(xs), op(init, f(xs[1])))

@inline _foldl(op, xs, init) = _mapfoldl(identity, op, xs, init)

_enumerate(xs) = _zip(ntuple(identity, length(xs)), xs)

_push(xs, x) = (xs..., x)
_map(f, xs) = _mapfoldl(f, _push, xs, ())

_cat() = ()
_cat(xs, tuples...) = (xs..., _cat(tuples...)...)

headtail(xs) = _headtail(xs...)
_headtail(x1, xs...) = (x1, xs)

newpartition(x, by) = (
    key = by(x),
    values = (x,),
    by = by,
)

function partitionby(values::Tuple, by)
    x, rest = headtail(values)
    return _foldl(_partitionstep, rest, (newpartition(x, by),))
end

_partitionstep(partitions, x) =
    if partitions[end].by(x) == partitions[end].key
        modify(partitions, @lens _[length(partitions)].values) do values
            _push(values, x)
        end
    else
        _push(partitions, newpartition(x, partitions[end].by))
    end


text:
## Example showing usage of the CUFFT wrappers
# Author: Stefan Kroboth

load("cumatrix.jl")

print("CUFFT version: $(cufftGetVersion())\n\n")

sx, sy = 16, 16
#sx, sy = 512, 512
#sx, sy = 1024, 1024

# create CUFFT plan
plan = cufftPlan2d(sx, sy, CUFFT_C2C)

# create data
h_data = ones(Complex64, sx, sy)
d_data = CuMatrix(h_data)  # from host to device
o_data = copy(d_data)
d_data_inv = copy(d_data)

# FFT
print("FFT:\n")
@time cufftExec(plan, d_data, o_data, CUFFT_FORWARD)
print(o_data)
println()

# IFFT
print("IFFT:\n")
@time cufftExec(plan, o_data, d_data_inv, CUFFT_INVERSE)
print(d_data_inv*float32(1.0/(sx*sy)))

max_stars_repo_name: stefan-k/julia-CuMatrix

max_issues_repo_name: stefan-k/julia-CuMatrix

max_forks_repo_name: stefan-k/julia-CuMatrix


text:
export PooledString
struct PooledString <: AbstractString end

struct PointerString <: AbstractString
    ptr::Ptr{UInt8}
    len::Int
end

function Base.hash(s::PointerString, h::UInt)
    h += Base.memhash_seed
    ccall(Base.memhash, UInt, (Ptr{UInt8}, Csize_t, UInt32), s.ptr, s.len, h % UInt32) + h
end

import Base: ==
function ==(x::String, y::PointerString)
    sizeof(x) == y.len && ccall(:memcmp, Cint, (Ptr{UInt8}, Ptr{UInt8}, Csize_t), pointer(x), y.ptr, y.len) == 0
end
==(y::PointerString, x::String) = x == y

Base.ncodeunits(s::PointerString) = s.len
@inline function Base.codeunit(s::PointerString, i::Integer)
    @boundscheck checkbounds(s, i)
    GC.@preserve s unsafe_load(s.ptr + i - 1)
end
Base.String(x::PointerString) = unsafe_string(x.ptr, x.len)

const TypeCode = Int8

# default value to signal that parsing should try to detect a type
const EMPTY       = 0b00000000 % TypeCode

# MISSING is a mask that can be combined w/ any other TypeCode for Union{T, Missing}
const MISSING     = 0b10000000 % TypeCode
missingtype(x::TypeCode) = (x & MISSING) === MISSING

# MISSINGTYPE is for a column like `Vector{Missing}`
# if we're trying to detect a column type and the 1st value of a column is `missing`
# we basically want to still treat it like EMPTY and try parsing other types on each row
const MISSINGTYPE = 0b00000001 % TypeCode

# enum-like type codes for basic supported types
const INT         = 0b00000010 % TypeCode
const FLOAT       = 0b00000011 % TypeCode
const DATE        = 0b00000100 % TypeCode
const DATETIME    = 0b00000101 % TypeCode
const TIME        = 0b00000110 % TypeCode
const BOOL        = 0b00000111 % TypeCode
const STRING      = 0b00001000 % TypeCode
const POOL        = 0b00010000 % TypeCode
pooled(x::TypeCode) = (x & POOL) == POOL

# a user-provided type; a mask that can be combined w/ basic types
const USER     = 0b00100000 % TypeCode
user(x::TypeCode) = (x & USER) === USER

const TYPEBITS = 0b00011111 % TypeCode
typebits(x::TypeCode) = x & TYPEBITS

typecode(::Type{Missing}) = MISSINGTYPE
typecode(::Type{<:Integer}) = INT
typecode(::Type{<:AbstractFloat}) = FLOAT
typecode(::Type{Date}) = DATE
typecode(::Type{DateTime}) = DATETIME
typecode(::Type{Time}) = TIME
typecode(::Type{Bool}) = BOOL
typecode(::Type{<:AbstractString}) = STRING
typecode(::Type{PooledString}) = POOL
typecode(::Type{CategoricalString{UInt32}}) = POOL
typecode(::Type{Union{}}) = EMPTY
typecode(::Type{Union{T, Missing}}) where {T} = typecode(T) | MISSING
typecode(::Type{T}) where {T} = EMPTY
typecode(x::T) where {T} = typecode(T)

const TYPECODES = Dict(
    EMPTY => Missing,
    MISSINGTYPE => Missing,
    INT => Int64,
    FLOAT => Float64,
    DATE => Date,
    DATETIME => DateTime,
    TIME => Time,
    BOOL => Bool,
    POOL => PooledString,
    STRING => String,
    INT | MISSING => Union{Int64, Missing},
    FLOAT | MISSING => Union{Float64, Missing},
    DATE | MISSING => Union{Date, Missing},
    DATETIME | MISSING => Union{DateTime, Missing},
    BOOL | MISSING => Union{Bool, Missing},
    POOL | MISSING => Union{PooledString, Missing},
    STRING | MISSING => Union{String, Missing}
)

gettypecodes(x::Dict) = Dict(typecode(k)=>typecode(v) for (k, v) in x)
gettypecodes(x::Dict{TypeCode, TypeCode}) = x

# bit patterns for missing value, int value, escaped string, position and len in tape parsing
const MISSING_BIT = 0x8000000000000000
missingvalue(x::UInt64) = (x & MISSING_BIT) == MISSING_BIT

const INT_BIT = 0x4000000000000000
intvalue(x::UInt64) = (x & INT_BIT) == INT_BIT

const ESCAPE_BIT = 0x2000000000000000
escapedvalue(x::UInt64) = (x & ESCAPE_BIT) == ESCAPE_BIT

getpos(x::UInt64) = (x & 0x1fffffffffff0000) >> 16
getlen(x::UInt64) = x & 0x000000000000ffff

# utilities to convert values to raw UInt64 and back for tape writing
int64(x::UInt64) = Core.bitcast(Int64, x)
float64(x::UInt64) = Core.bitcast(Float64, x)
bool(x::UInt64) = x == 0x0000000000000001
date(x::UInt64) = Date(Dates.UTD(int64(x)))
datetime(x::UInt64) = DateTime(Dates.UTM(int64(x)))
time(x::UInt64) = Time(Nanosecond(int64(x)))
ref(x::UInt64) = unsafe_trunc(UInt32, x)

uint64(x::Int64) = Core.bitcast(UInt64, x)
uint64(x::Float64) = Core.bitcast(UInt64, x)
uint64(x::Bool) = UInt64(x)
uint64(x::Union{Date, DateTime, Time}) = uint64(Dates.value(x))
uint64(x::UInt32) = UInt64(x)

reinterp_func(::Type{Int64}) = int64
reinterp_func(::Type{Float64}) = float64
reinterp_func(::Type{Date}) = date
reinterp_func(::Type{DateTime}) = datetime
reinterp_func(::Type{Time}) = time
reinterp_func(::Type{Bool}) = bool

@noinline function consumeBOM!(source)
    # BOM character detection
    startpos = pos = 1
    len = length(source)
    if pos <= len && source[pos] == 0xef
        pos += 1
        if pos <= len && source[pos] == 0xbb
            pos += 1
        else
            pos = startpos
        end
        if pos <= len && source[pos] == 0xbf
            pos += 1
        else
            pos = startpos
        end
    end
    return pos
end

function getsource(source, use_mmap)
    if source isa Vector{UInt8}
        return source
    elseif source isa Cmd
        return Base.read(source)
    elseif use_mmap && !isa(source, IO)
        return Mmap.mmap(source)
    elseif !isa(source, IO)
        m = Mmap.mmap(source)
        m2 = Mmap.mmap(Vector{UInt8}, length(m))
        copyto!(m2, 1, m, 1, length(m))
        finalize(m)
        return m2
    else
        iosource = source isa IO ? source : open(String(source))
        io = IOBuffer()
        while !eof(iosource)
            Base.write(io, iosource)
        end
        A = Mmap.mmap(Vector{UInt8}, io.size)
        copyto!(A, 1, io.data, 1, io.size)
        return A
    end
end

getname(buf::Vector{UInt8}) = "<raw buffer>"
getname(str) = String(str)
getname(io::I) where {I <: IO} = string("<", I, ">")

const RESERVED = Set(["local", "global", "export", "let",
    "for", "struct", "while", "const", "continue", "import",
    "function", "if", "else", "try", "begin", "break", "catch",
    "return", "using", "baremodule", "macro", "finally",
    "module", "elseif", "end", "quote", "do"])

normalizename(name::Symbol) = name
function normalizename(name::String)::Symbol
    uname = strip(Unicode.normalize(name))
    id = Base.isidentifier(uname) ? uname : map(c->Base.is_id_char(c) ? c : '_', uname)
    cleansed = string((isempty(id) || !Base.is_id_start_char(id[1]) || id in RESERVED) ? "_" : "", id)
    return Symbol(replace(cleansed, r"(_)\1+"=>"_"))
end

function makeunique(names)
    set = Set(names)
    length(set) == length(names) && return Symbol[Symbol(x) for x in names]
    nms = Symbol[]
    for nm in names
        if nm in nms
            k = 1
            newnm = Symbol("$(nm)_$k")
            while newnm in set || newnm in nms
                k += 1
                newnm = Symbol("$(nm)_$k")
            end
            nm = newnm
        end
        push!(nms, nm)
    end
    return nms
end

initialtypes(T, x::AbstractDict{String}, names) = TypeCode[haskey(x, string(nm)) ? typecode(x[string(nm)]) | USER : T for nm in names]
initialtypes(T, x::AbstractDict{Symbol}, names) = TypeCode[haskey(x, nm) ? typecode(x[nm]) | USER : T for nm in names]
initialtypes(T, x::AbstractDict{Int}, names)    = TypeCode[haskey(x, i) ? typecode(x[i]) | USER : T for i = 1:length(names)]

function timetype(df::Dates.DateFormat)
    date = false
    time = false
    for token in df.tokens
        T = typeof(token)
        if T == Dates.DatePart{'H'}
            time = true
        elseif T == Dates.DatePart{'y'} || T == Dates.DatePart{'Y'}
            date = true
        end
    end
    return ifelse(date & time, DateTime, ifelse(time, Time, Date))
end

roundup(a, n) = (a + (n - 1)) & ~(n - 1)

function unescape(s, e)
    n = ncodeunits(s)
    buf = Base.StringVector(n)
    len = 1
    i = 1
    @inbounds begin
        while i <= n
            b = codeunit(s, i)
            if b == e
                i += 1
                b = codeunit(s, i)
            end
            @inbounds buf[len] = b
            len += 1
            i += 1
        end
    end
    resize!(buf, len - 1)
    return String(buf)
end

"""
    CSV.detect(str::String)

Use the same logic used by `CSV.File` to detect column types, to parse a value from a plain string.
This can be useful in conjunction with the `CSV.Rows` type, which returns each cell of a file as a String.
The order of types attempted is: `Int64`, `Float64`, `Date`, `DateTime`, `Bool`, and if all fail, the input String is returned.
No errors are thrown.
For advanced usage, you can pass your own `Parsers.Options` type as a keyword argument `option=ops` for sentinel value detection.
"""
function detect end

function detect(str::String; options=Parsers.OPTIONS)
    buf = codeunits(str)
    pos = 1
    len = sizeof(str)
    int, code, vpos, vlen, tlen = Parsers.xparse(Int64, buf, pos, len, options)
    if Parsers.sentinel(code) && code > 0
        return missing
    end
    if Parsers.ok(code) && vlen == len
        return int
    end
    float, code, vpos, vlen, tlen = Parsers.xparse(Float64, buf, pos, len, options)
    if Parsers.ok(code) && vlen == len
        return float
    end
    if options.dateformat === nothing
        try
            date, code, vpos, vlen, tlen = Parsers.xparse(Date, buf, pos, len, options)
            if Parsers.ok(code) && vlen == len
                return date
            end
        catch e
        end
        try
            datetime, code, vpos, vlen, tlen = Parsers.xparse(DateTime, buf, pos, len, options)
            if Parsers.ok(code) && vlen == len
                return datetime
            end
        catch e
        end
    else
        try
            # use user-provided dateformat
            DT = timetype(options.dateformat)
            dt, code, vpos, vlen, tlen = Parsers.xparse(DT, buf, pos, len, options)
            if Parsers.ok(code) && vlen == len
                return dt
            end
        catch e
        end
    end
    bool, code, vpos, vlen, tlen = Parsers.xparse(Bool, buf, pos, len, options)
    if Parsers.ok(code) && vlen == len
        return bool
    end
    return str
end


text:
using Distributed
@everywhere using Test
@everywhere using FileTrees
@everywhere using FileTrees: Thunk
@everywhere using Dates
@everywhere import FilePathsBase: /

@testset "maketree" begin
    global t = maketree(["a" => ["b" => ["a"],
                                 "c" => ["b", "a", "c"=>[]]]])

    @test strip(repr(t)) == """
                            ./
                            └─ a/
                               ├─ b/
                               │  └─ a
                               └─ c/
                                  ├─ b
                                  ├─ a
                                  └─ c/"""
end

@testset "indexing" begin
    @test name(t) == "."
    @test t["a"]["b"]["a"] isa File
    @test t["a"]["c"]["c"] isa FileTree

    @test path(t["a"]["b"]["a"])  == p"." /"a"/"b"/"a"

    t1 = FileTrees.rename(t, "foo")
    @test path(t1["a"]["b"]["a"])  == p"foo" /"a"/"b"/"a"

    @test isequal(t[r"a|b"], t)

    @test isempty(t[r"^c"])

    @test isequal(t["a"]["c"][["b", r"a|d"]],
                  maketree("c" => ["b","a"]))

    @test_throws ErrorException t["c"]
end

@testset "filter" begin
    @test isequal(filter(f->f isa FileTree || name(f) == "c", t),
                  maketree(["a"=>["b"=>[],
                                  "c"=>["c"=>[]]]]))
end

@testset "merge" begin
    @test isequal(merge(t,t), t)
end

@testset "diff" begin
    @test isempty(FileTrees.diff(t,t))
    @test isequal(FileTrees.diff(t,maketree([])), t)
end

import FileTrees: attach

@testset "touch mv and cp" begin
    global t1 = maketree([])
    for yr = 1:9
        for mo = 1:6
            t1 = attach(t1, joinpath(string("0", yr), string("0", mo)),
                        File(nothing, "data.csv", (yr, mo)) )
        end
    end
    data1 = reducevalues(vcat, t1) |> exec

    t4 = mv(t1, r"^(.*)/(.*)/data.csv$", s"\1/\2.csv")

    t5 = maketree([])
    for yr = 1:9
        for mo = 1:6
            t5 = attach(t5, string("0", yr),
                        File(nothing, string("0", mo) * ".csv", (yr, mo)))
        end
    end

    @test isequal(t4, t5)

    data2 = reducevalues(vcat, t4) |> exec

    @test data1 == data2

    t7 = mapsubtrees(t1, glob"*") do subtree
        reducevalues(vcat, subtree)
    end

    @test reducevalues(hcat, t7) == [(j, i) for i=1:6, j=1:9]

    t6 = cp(t1, r"^(.*)/(.*)/data.csv$", s"\1/\2.csv")

    @test isequal(t6, merge(t1, t5))
end

@testset "values" begin
    t1 = FileTrees.load(x->string(path(x)), t)
    if isdir("test_dir")
        rm("test_dir", recursive=true)
    end

    @test t1["a/b/a"][] == string(p"." / "a" / "b" / "a")

    @test reducevalues(*, mapvalues(lowercase, t1)) == lowercase(reducevalues(*, t1))

    FileTrees.save(maketree("test_dir" => [t1])) do f
        @test f isa File
        open(path(f), "w") do io
            print(io, f[])
        end
    end

    t2 = FileTree("test_dir")
    t3 = FileTrees.load(t2) do f
        open(path(f), "r") do io
            String(read(io))
        end
    end

    t4 = filter(!isempty, t1)

    @test isequal(t3, FileTrees.rename(t4, "test_dir"))
    if isdir("test_dir")
        rm("test_dir", recursive=true)
    end
end

@testset "lazy-exec" begin

    if isdir("test_dir_lazy")
        rm("test_dir_lazy", recursive=true)
    end

    t1 = FileTrees.load(x->uppercase(string(path(x))), t, lazy=true)

    @test t1["a/b/a"][] isa Thunk
    @test exec(t1)["a/b/a"][] == string(p"."/"A"/"B"/"A")

    @test exec(reducevalues(*, mapvalues(lowercase, t1))) == lowercase(exec(reducevalues(*, t1)))

    s = FileTrees.save(maketree("test_dir_lazy" => [t1])) do f
        open(path(f), "w") do io
            print(io, f[])
        end
    end

    @test isdir("test_dir_lazy")
    @test isfile("test_dir_lazy/a/b/a")

    t2 = FileTree("test_dir_lazy")
    t3 = FileTrees.load(t2; lazy=true) do f
        open(path(f), "r") do io
            (String(read(io)), now())
        end
    end
    toc = now()
    sleep(0.01)
    tic = exec(reducevalues((x,y)->x, mapvalues(last, t3)))

    @test tic > toc

    t4 = filter(!isempty, t1) |> exec

    t5 = mapvalues(first, t3) |> exec
    @test isequal(t5, FileTrees.rename(t4, "test_dir_lazy"))

    if isdir("test_dir_lazy")
        rm("test_dir_lazy", recursive=true)
    end
end


text:
struct Switch{C, N, K, T} <: GenerativeFunction{T, Trace}
    branches::NTuple{N, GenerativeFunction{T}}
    cases::Dict{C, Int}
    function Switch(gen_fns::GenerativeFunction...)
        @assert !isempty(gen_fns)
        rettype = get_return_type(getindex(gen_fns, 1))
        new{Int, length(gen_fns), typeof(gen_fns), rettype}(gen_fns, Dict{Int, Int}())
    end
    function Switch(d::Dict{C, Int}, gen_fns::GenerativeFunction...) where C
        @assert !isempty(gen_fns)
        rettype = get_return_type(getindex(gen_fns, 1))
        new{C, length(gen_fns), typeof(gen_fns), rettype}(gen_fns, d)
    end
end
export Switch

has_argument_grads(switch_fn::Switch) = map(zip(map(has_argument_grads, switch_fn.branches)...)) do as
    all(as)
end
accepts_output_grad(switch_fn::Switch) = all(accepts_output_grad, switch_fn.branches)

function (gen_fn::Switch)(index::Int, args...)
    (_, _, retval) = propose(gen_fn, (index, args...))
    retval
end

function (gen_fn::Switch{C})(index::C, args...) where C
    (_, _, retval) = propose(gen_fn, (gen_fn.cases[index], args...))
    retval
end

include("trace.jl")
include("assess.jl")
include("propose.jl")
include("simulate.jl")
include("generate.jl")
include("update.jl")
include("regenerate.jl")
include("backprop.jl")

@doc(
"""
    gen_fn = Switch(gen_fns::GenerativeFunction...)

Returns a new generative function that accepts an argument tuple of type `Tuple{Int, ...}` where the first index indicates which branch to call.

    gen_fn = Switch(d::Dict{T, Int}, gen_fns::GenerativeFunction...) where T

Returns a new generative function that accepts an argument tuple of type `Tuple{Int, ...}` or an argument tuple of type `Tuple{T, ...}` where the first index either indicates which branch to call, or indicates an index into `d` which maps to the selected branch. This form is meant for convenience - it allows the programmer to use `d` like if-else or case statements.

`Switch` is designed to allow for the expression of patterns of if-else control flow. `gen_fns` must satisfy a few requirements:

1. Each `gen_fn` in `gen_fns` must accept the same argument types.
2. Each `gen_fn` in `gen_fns` must return the same return type.

Otherwise, each `gen_fn` can come from different modeling languages, possess different traces, etc.
""", Switch)


text:
module ElmerResultReader

export elmervar, readelmervars, readelmerdat, readlinearsystem

type elmervar
  perm::Array{Int64,2}
  permsize::Array{Int64,1}
  values::Array{Float64,2}
  name
end

function readmatch(file::IOStream, matcher; error_msg="Invalid parse")
  if (!(matcher(readline(file))))
    close(file)
    throw(ParseError(error_msg))
  end
end

function readmatch(data::AbstractString, matcher; error_msg="Invalid parse")
  if (!(matcher(data)))
    throw(ParseError(error_msg))
  end
end

function readnlines(data::IOStream, n::Integer)
  outbuf = IOBuffer()
  #=print("reading $(n) lines\n")=#
  for i = 1:n
    line = readline(data)
    write(outbuf, line)
  end
  #=print("ok\n")=#
  return outbuf
end

@doc """
function readelmervars(filename; verbose=1)

Populates elmervar array from results given in ```filename```.
Works only with 1 timestep currently.
""" ->
function readelmervars(filename; verbose=1)
  verbose > 1 ? print("Reading variables from $(filename)\n") : Union{}
  totaldofmatch(x) = x[1:12] == " Total DOFs:"

  charmatch(x,y) = x[1] == y
  commentmatch(x) = charmatch(x,'!')

  ms_t1 = open(filename)
  readmatch(ms_t1, x -> x == " ASCII 3\n")
  readmatch(ms_t1, commentmatch)
  readmatch(ms_t1, x-> x == " Degrees of freedom: \n")

  # read dof table
  line = readline(ms_t1)
  dofdata = []
  while(!totaldofmatch(line))
    readmatch(line, x->ismatch(r".+:.+:", x), error_msg= "Couldn't get dof line")
    dofline = readdlm(IOBuffer(line), ':', ASCIIString)
    dofs = readdlm(IOBuffer(rstrip(dofline[2])), Int64)
    push!(dofdata, [rstrip(dofline[1]) dofs lstrip(dofline[3])])
    line = readline(ms_t1)
  end

  nvars = parse(Int64, match(r".+:\s*(\d*)", line).captures[1])
  line = readline(ms_t1)

  readmatch(line, x -> " Number Of Nodes" == x[1:16])
  nnodes = parse(Int64, match(r".+:\s*(\d*)", line).captures[1])

  t_elmervars = []

  while (true)
    try
      line = readline(ms_t1)
      readmatch(line, x -> "Time:" == x[1:5])
    catch readmatch_error
      if isa(readmatch_error, ParseError)
        print(readmatch_error.msg)
        break
      else
        if isa(readmatch_error, BoundsError)
          verbose > 2 ? print("Did not find next timestep, exiting\n") : Union{}
          break
        end
      end
    end
    timestep_n = parse(Int64, match(r"^Time:\s*(\d+)\s*(\d+)\s*", line).captures[1])
    verbose>0 ? print("Reading timestep $(timestep_n)") : Union{}

    elmervars = []
    perm = []
    prev_ndof = 0
    for dof in dofdata
      line = readline(ms_t1)
      readmatch(line, x -> ismatch(Regex(dof[1]), line), error_msg="Error: $(line)")
      line = readline(ms_t1)
      readmatch(line, x-> x[1:5] == "Perm:")
      if (!ismatch(r"use previous", line))
        prev_ndof = parse(Int64, match(r".+:\s*\d+\s+(\d+)",line).captures[1])
        dof[2] = prev_ndof
        lines = readnlines(ms_t1, dof[2])
        seekstart(lines)
        perm = readdlm(lines, Int64)
      else
        dof[2] = prev_ndof
      end
      lines = readnlines(ms_t1, dof[2])
      seekstart(lines)
      vals = readdlm(lines, Float64)
      push!(elmervars, elmervar(perm, dof[2:4], vals, dof[1]))
      verbose > 0 ? print(".") : Union{}
    end

    push!(t_elmervars, elmervars)
    verbose > 0 ? print(" Read $(length(elmervars)) variables.\n") : print(".\n")
  end

  return t_elmervars

end

@doc """
function readelmerdat(name, NPROC=1, DICT=false)
Reads .dat file generated with savescalars
""" ->
function readelmerdat(name, NPROC=1, DICT=false)
  datas = Array[]
  if NPROC > 1
    for N = 1:NPROC
      push!(datas, readdlm("$name.$(N-1)"))
    end
  else
    push!(datas, readdlm(name))
  end
  desc = open(readlines, "$name.names")[4:end]
  desc = map(strip, desc)
  if DICT
    datasum = sum(datas)
    datadict = Dict{ASCIIString, Array{Float64,1}}()
    k = 1
    for s in desc
      datadict[s] = datasum[:,k]
      k = k + 1
    end
    return datasum, desc, datadict
  end
  return sum(datas), desc
end

@doc """
function readlinearsystem(prefix; suffix="")

Loads elmer linear system that is saved with SaveLinearSystem.

* Output: A::SparseMatrixCSC{Float64, Int64}, b::Array{Float64,1}
* Input: read data from <prefix>_a.dat and <prefix>_b.dat
""" ->
function readlinearsystem(prefix; suffix="")
  ri = x -> round(Int64, x)
  Alist = readdlm(string(prefix, "_a.dat", suffix))
  blist = readdlm(string(prefix, "_b.dat", suffix))
  A = sparse(ri(Alist[:,1]), ri(Alist[:,2]), Alist[:,3])
  b = blist[:,2]
  return A,b
end

end # module ElmerResultReader


text:
using Test
using JustSayIt

functions = [:internet]

@testset "$(basename(@__FILE__))" begin
    @testset "Module Internet available" begin
        @test @isdefined Internet
    end
    @testset "function $f available" for f in functions
        @test f in names(Internet; all=true)
    end;
end;


text:
"""
Abstract type for coordinates in three dimensions

All sub types have x, y, z coordinates.
And conversion is available between subtypes using the convert function.

```julia
bv_coord = (0.3, 2, 3.1)
tal_coord = convert(Talairach, mni)
```
"""
abstract type Coordinate end

"""
Type for BrainVision coordinate system.
"""
mutable struct BrainVision <: Coordinate
    x::typeof(1.0u"m")
    y::typeof(1.0u"m")
    z::typeof(1.0u"m")
    BrainVision(x::AbstractQuantity, y::AbstractQuantity, z::AbstractQuantity) =
        new(x, y, z)
    BrainVision(x::Number, y::Number, z::Number) = new(x * u"m", y * u"m", z * u"m")
end

"""
Type for Talairach coordinate system.
"""
mutable struct Talairach <: Coordinate
    x::typeof(1.0u"m")
    y::typeof(1.0u"m")
    z::typeof(1.0u"m")
    Talairach(x::AbstractQuantity, y::AbstractQuantity, z::AbstractQuantity) = new(x, y, z)
    Talairach(x::Number, y::Number, z::Number) = new(x * u"m", y * u"m", z * u"m")
end

"""
Type for SPM coordinate system.
"""
mutable struct SPM <: Coordinate
    x::typeof(1.0u"m")
    y::typeof(1.0u"m")
    z::typeof(1.0u"m")
    SPM(x::AbstractQuantity, y::AbstractQuantity, z::AbstractQuantity) = new(x, y, z)
    SPM(x::Number, y::Number, z::Number) = new(x * u"m", y * u"m", z * u"m")
end

"""
Type to be used when the coordinate system is unknown.
"""
mutable struct UnknownCoordinate <: Coordinate
    x::typeof(1.0u"m")
    y::typeof(1.0u"m")
    z::typeof(1.0u"m")
    UnknownCoordinate(x::AbstractQuantity, y::AbstractQuantity, z::AbstractQuantity) =
        new(x, y, z)
    UnknownCoordinate(x::Number, y::Number, z::Number) = new(x * u"m", y * u"m", z * u"m")
end

import Base.show
function show(c::S) where {S<:Coordinate}
    println(
        "Coordinate: $(typeof(c)) - ($(c.x |> u"mm"), $(c.y |> u"mm"), $(c.z |> u"mm"))",
    )
end


text:
struct CFLPFacility
    id::Int64
    capacity::Int64
    fixed_cost::Float64
    costs::Vector{Float64}
end

function Base.show(io::IO, facility::CFLPFacility)
    print(io, "F($(facility.id))")
end

struct CFLPCustomer
    id::Int64
    demand::Int64
    costs::Vector{Float64}
end

function Base.show(io::IO, customer::CFLPCustomer)
    print(io, "C($(customer.id))")
end

struct CFLPData
    name::String
    facilities::Vector{CFLPFacility}
    customers::Vector{CFLPCustomer}

    lb::Float64
    ub::Float64
end

function Base.show(io::IO, data::CFLPData)
    print(io, "CFLP Data $(data.name)")
    print(io, " ($(length(data.facilities)) facilities,")
    print(io, " $(length(data.customers)) customers)")
    print(io, " [$(data.lb), $(data.ub)]")
end


text:
# This file was generated, do not modify it. # hide
model = @pipeline(Standardizer(),
                     OneHotEncoder(),
                     LinearRegressor())

pipe  = machine(model, Xc, y)
fit!(pipe, rows=train)
ŷ = predict(pipe, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)


text: @deprecate axes_keys(x) keys.(axes(x))


text:
function render!(
    rendermodel::RenderModel,
    scene::Scene;
    car_color::Colorant=COLOR_CAR_OTHER,
    car_colors::Dict{Int,Colorant}=Dict{Int,Colorant}(), #  id -> color
    )

    for veh in scene
        render!(rendermodel, veh, get(car_colors, veh.def.id, car_color))
    end

    rendermodel
end

function render(ctx::CairoContext, scene::Scene, roadway::Roadway;
    rendermodel::RenderModel=RenderModel(),
    cam::Camera=SceneFollowCamera(),
    car_colors::Dict{Int,Colorant}=Dict{Int,Colorant}(),
    )

    canvas_width = floor(Int, Cairo.width(ctx))
    canvas_height = floor(Int, Cairo.height(ctx))

    clear_setup!(rendermodel)

    render!(rendermodel, roadway)
    render!(rendermodel, scene, car_colors=car_colors)

    camera_set!(rendermodel, cam, scene, roadway, canvas_width, canvas_height)

    render(rendermodel, ctx, canvas_width, canvas_height)
    ctx
end
function render(scene::Scene, roadway::Roadway;
    canvas_width::Int=DEFAULT_CANVAS_WIDTH,
    canvas_height::Int=DEFAULT_CANVAS_HEIGHT,
    rendermodel::RenderModel=RenderModel(),
    cam::Camera=SceneFollowCamera(),
    car_colors::Dict{Int,Colorant}=Dict{Int,Colorant}(), # id
    )

    s, ctx = get_surface_and_context(canvas_width, canvas_height)
    render(ctx, scene, roadway, rendermodel=rendermodel, cam=cam, car_colors=car_colors)
    s
end


text:
using Documenter
using IndirectLikelihood
using ContinuousTransformations

"""
Documentation generation wrapped into a single function for development.

Call `mkdocs()` repeatedly when working on the documentation interactively.
"""
mkdocs() = makedocs(format = :html,
                    sitename = "IndirectLikelihood",
                    pages = ["index.md",
                             "general.md",
                             "auxiliary.md",
                             "internals.md"])

mkdocs()

deploydocs(repo = "github.com/tpapp/IndirectLikelihood.jl.git",
           target = "build",
           deps = nothing,
           make = nothing,
           julia = "0.6")


text:
################################################################################
##                                 CALLBACKS                                  ##
################################################################################

# librdkafka provides several callback (CB) functions such as delivery report CB and error CB
# we defne fixed @cfunction callbacks to pass them during compilation
# these @cfunction's then call corresponding Julia functions from XXX_CALLBACKS dictionaries
const DELIVERY_CALLBACKS = Dict{Ptr{Cvoid}, Function}()
function delivery_callback(rk::Ptr{Cvoid}, rkmessage::Ptr{CKafkaMessage}, opaque::Ptr{Cvoid})::Cvoid
    msg = Message(unsafe_load(rkmessage))
    if haskey(DELIVERY_CALLBACKS, rk)
        cb = DELIVERY_CALLBACKS[rk]
        cb(msg)
    end
end
delivery_callback_c = @cfunction(delivery_callback, Cvoid, (Ptr{Cvoid}, Ptr{CKafkaMessage}, Ptr{Cvoid}))

function error_callback(rk::Ptr{Cvoid}, err::Cint, reason::Ptr{Cchar}, opaque::Ptr{Cvoid})::Cvoid
    # TODO: implement
    @warn "Error callback is not implemented"
end
error_callback_c = @cfunction(error_callback, Cvoid, (Ptr{Cvoid}, Cint, Ptr{Cvoid}, Ptr{Cvoid}))
const ERROR_CALLBACKS = Dict{Ptr, Function}()

mutable struct KafkaClient
    conf::Dict{Any, Any}
    typ::Cint
    rk::Ptr{Cvoid}
end

function KafkaClient(typ::Integer, conf::Dict=Dict(); dr_cb=nothing, err_cb=nothing)
    # TODO: add options for callbacks, creata a task calling kafka_poll(rk, timeout) every 1 second
    c_conf = kafka_conf_new()
    for (k, v) in conf
        kafka_conf_set(c_conf, string(k), string(v))
    end
    if dr_cb !== nothing
        # set callback in config before creating rk
        kafka_conf_set_dr_msg_cb(c_conf, delivery_callback_c)
    end
    if err_cb !== nothing
        kafka_conf_set_error_cb(c_conf, error_callback_c)
    end
    rk = kafka_new(c_conf, Cint(typ))
    client = KafkaClient(conf, typ, rk)
    if dr_cb !== nothing
        # set Julia callback after rk is created
        DELIVERY_CALLBACKS[rk] = dr_cb
    end
    if err_cb !== nothing
        ERROR_CALLBACKS[rk] = err_cb
    end
    return client
end

Base.show(io::IO, kc::KafkaClient) = print(io, "KafkaClient($(kc.typ))")

mutable struct KafkaTopic
    conf::Dict{Any, Any}
    topic::String
    rkt::Ptr{Cvoid}
end

function KafkaTopic(kc::KafkaClient, topic::String, conf::Dict=Dict())
    c_conf = kafka_topic_conf_new()
    for (k, v) in conf
        kafka_topic_conf_set(c_conf, string(k), string(v))
    end
    rkt = kafka_topic_new(kc.rk, topic, c_conf)
    topic = KafkaTopic(conf, topic, rkt)
    return topic
end

function Base.close(topic::KafkaTopic)
    kafka_topic_destroy(topic.rkt)
end

Base.show(io::IO, kt::KafkaTopic) = print(io, "KafkaTopic($(kt.topic))")

struct Message{K,P}
    err::Int
    topic::KafkaTopic
    partition::Int32
    key::Union{K, Cvoid}
    payload::Union{P, Cvoid}
    offset::Int64
end

Base.convert(::Type{String}, data::Vector{UInt8}) = String(data)

function Message{K,P}(c_msg::CKafkaMessage) where {K,P}
    topic = KafkaTopic(Dict(), "<unkown>", c_msg.rkt)
    if c_msg.err == 0
        key = convert(K, unsafe_load_array(c_msg.key, c_msg.key_len))
        payload = convert(P, unsafe_load_array(c_msg.payload, c_msg.len))
        return Message(Int(c_msg.err), topic, c_msg.partition, key, payload, c_msg.offset)
    else
        return Message{K,P}(Int(c_msg.err), topic, c_msg.partition, nothing, nothing, c_msg.offset)
    end
end
Message(c_msg::CKafkaMessage) = Message{Vector{UInt8}, Vector{UInt8}}(c_msg)

Base.show(io::IO, msg::Message) = print(io, "Message($(msg.key): $(msg.payload))")


text:
clear()
model = Chain(
IndRNN(2,128,relu),
IndRNN(128,64,relu),
Dense(64,1,relu)
)

params = paramsof(model)
optimi = Adam(params;lr=1e-4)

B = 32
T = 100
x, s = addingProblemData(B,T)

epochs = 1000
lossval = zeros(epochs,1)
for e=1:epochs
    resethidden(model)
    for t = 1:T-1
        tmp = forward(model, Variable( reshape(x[:,t,:], 2,B) ) )
        zeroDelta(tmp)
    end
    y = forward(model, Variable( reshape(x[:,T,:], 2,B) ) )
    COST = mseCost(y, Variable( reshape(s,1,B) ) )
    backward()
    update!(optimi, params; clipvalue=1e4)
    zerograds!(params)
    lossval[e] = COST/B
end

lineplot(vec(lossval))

function addingProblemData(N::Int,T::Int)
    #= 两个长度为 T 的序列组成输入：
    第一个序列在 (0,1)范围内均匀采样，
    第二个序列的前两个为 1，其余都为 0.
    例如：
    x1 = 0.90  0.93  0.90  0.56  0.46
    x2 = 1.0   1.0   0.0   0.0   0.0
    而输出为 sum(x1 .* x2) = 0.90 + 0.93
    =#
    @assert (T>1) "The sequence length should lager than 1"
    x1 =  rand(1,T,N)./3
    x2 = zeros(1,T,N)
    x2[1,1,:] .= 1.0
    x2[1,2,:] .= 1.0
    y = sum(x1 .* x2, dims=2)
    return [x1;x2],y
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
# My first attempt at comming up with efficient connected components

function groupTRUES( mask::Array{Bool,2} )

	h, w = size( mask );

	finalGroups   = [];
	prevColGroups = [];

	for col in 1:w

		# Finding connected TRUES ( groups ) in the current column, 1D

	    colGroups = groupsInColumn( view( mask, :, col ), col, length(finalGroups) );
	    numGroups = length( colGroups );

	    # For each group in colGroups, we have to check if it overlaps with any group in the previous column
	    # Since the groups are ordered from "top" to "bottom", we use the variable "start" to skip all groups
	    # in the previous column that are above the currently considered group in colGroups

	    start = 1;

	    for idx1 in 1:numGroups

	        thisGroup = colGroups[idx1];
	        overlaps  = false;

	        for idx2 in start:length(prevColGroups)

	            prevGroup = prevColGroups[idx2]

	            if     thisGroup[1] > prevGroup[2] # thisGroup is above prevGroup

	                start = idx2+1;

	            elseif thisGroup[2] < prevGroup[1] # thisGroup is below prevGroup
	                break;

	            else                               # thisGroup and prevGroup overlap

	                # adding coordinates of thisGroup to prevGroup's parent
	                push!( finalGroups[ prevGroup[4] ], thisGroup[1:3] );

	                # changing thisGroup's parent to prevGroup's parent
	                colGroups[idx1][4] = prevGroup[4];
	                diditoverlap = true;
	            end
	        end

	        # If thisGroup does not overlap, it will be added to the list of groups
	        # and have the change to nucleate its own group

	        if !overlaps
	            push!( finalGroups, [ thisGroup, ] )
	        end
	    end
	    prevColGroups = colGroups;
	end

	return finalGroups
end

function groupsInColumn( column::AbstractArray{Bool,1}, colIdx, last )

	groups    = [];
	prevPixel = false;

    for row in 1:length(column)
        pixel = column[ row ]
        if pixel
            if prevPixel
                groups[end][2] = row;
            else
                push!( groups, [ row, row, colIdx, last+1 ] );
            end
        end
        prevPixel = pixel;
    end

    return groups;
end


text:
__precompile__()

module Collective

import Base: size, getindex, isless, show
using Iterators: subsets
using DataStructures: OrderedDict

export Corpus,
       wordlist,
       analyze,
       best_feature,
       best_cluster,
       best_clusters,
       common_features

cleanup(s::String) = replace(lowercase(strip(s)), r"[^a-z]", "")

wordlist(data::IO) = cleanup.(vec(readdlm(data, ',', String)))
wordlist(s::String) = wordlist(IOBuffer(s))

include("bitstally.jl")
import .BitsTallies: BitsTally, isanagram, istransaddition
include("feature_expressions.jl")
include("features.jl")

type Corpus{F}
    features::FeatureSet{F}
    frequencies::Vector{Float64}
end

show(io::IO, c::Corpus) = print(io, "Corpus with $(length(c.features.descriptions)) features")

function Corpus(words::AbstractArray{String}, features=allfeatures())
    featureset = FeatureSet(features)
    frequencies = sum(featureset.evaluate(lowercase(word)) for word in words) / length(words)
    frequencies = max.(frequencies, 1/length(words))
    Corpus(featureset, frequencies)
end

immutable FeatureResult
    description::String
    evaluate::Function
    satisfied::BitArray{1}
    probability::Float64
end

(r::FeatureResult)(word::String) = r.evaluate(word)

show(io::IO, r::FeatureResult) = print(io, "Feature: $(sum(r.satisfied))/$(length(r.satisfied)) matches, p=$(r.probability): $(r.description)")

isless(f1::FeatureResult, f2::FeatureResult) = f1.probability < f2.probability

binomial_probability(successes, trials, frequency) = binomial(trials, successes) * frequency^successes * (1 - frequency)^(trials - successes)

function binomial_pvalue(successes, trials, frequency)
    p = binomial_probability(successes, trials, frequency)
    expected = frequency * trials
    if successes > expected
        for k in (successes + 1):trials
            p += binomial_probability(k, trials, frequency)
        end
    else
        for k in 1:(successes - 1)
            p += binomial_probability(k, trials, frequency)
        end
    end
    p
end

function evaluate(corpus::Corpus, words::AbstractArray{String})
    corpus.features.evaluate.(lowercase.(words))
end

function analyze(corpus::Corpus, words::AbstractArray{String}, allowed_misses=length(words))
    match = evaluate(corpus, words)
    num_matches = sum(match)
    probabilities = binomial_pvalue.(num_matches, length(words), corpus.frequencies)
    results = FeatureResult.(corpus.features.descriptions, corpus.features.evaluators, [[m[i] for m in match] for i in 1:length(corpus.features.descriptions)], probabilities)
    required_matches = length(words) - allowed_misses
    [r for r in results if sum(r.satisfied) >= required_matches]
end

function best_feature(corpus::Corpus, words::AbstractArray{String})
    best_feature(corpus, evaluate(corpus, words))
end

function best_feature(corpus::Corpus, vals::AbstractArray{BitArray{1}})
    matches = sum(vals)
    n = length(vals)
    num_features = length(corpus.features.descriptions)
    p, i = findmin(Collective.binomial_pvalue(matches[i], n, corpus.frequencies[i]) for i in 1:num_features)
    FeatureResult(corpus.features.descriptions[i],
                  corpus.features.evaluators[i],
                  [m[i] for m in vals],
                  p)
end

immutable Cluster
    words::Vector{String}
    feature::FeatureResult
end

isless(c1::Cluster, c2::Cluster) = c1.feature < c2.feature

function best_cluster(corpus::Corpus, words::AbstractArray{String}, n::Integer)
    vals = evaluate(corpus, words)
    x = Inf
    best_subset = Int[]
    for subset in subsets(collect(1:length(vals)), n)
        f = Collective.best_feature(corpus, vals[subset])
        p = f.probability
        if p < x
            x = p
            best_subset = subset
        end
    end
    Cluster(words[best_subset], best_feature(corpus, words[best_subset]))
end

best_clusters(corpus::Corpus, words::AbstractArray{String}) =
    [best_cluster(corpus, words, n) for n in 2:div(length(words), 2) if length(words) % n == 0]

best_cluster(corpus::Corpus, words::AbstractArray{String}) = minimum(best_clusters(corpus, words))

function common_features(featureset::FeatureSet, words::AbstractArray{String})
    vals = featureset.evaluate.(lowercase.(words))
    common_indices = [i for i in 1:length(vals[1]) if all(v[i] for v in vals)]
    FeatureResult[
        FeatureResult(featureset.descriptions[i],
                      featureset.evaluators[i],
                      [m[i] for m in vals],
                      NaN) for i in common_indices]
end

common_features(corpus::Corpus, words::AbstractArray{String}) =
    common_features(corpus.features, words)

end


text:
import Base: getindex, setindex!, enumerate, call, length
import DataStructures: SortedDict, pastendtoken, deref, deref_value, startof, advance

export Chromosome, DenseChromosome, IntervalChromosome, setoffspring!, inherit

abstract Chromosome{Gene}

type DenseChromosome{Gene} <: Chromosome{Gene}
    data::Vector{Gene}
end
call{T}(::Type{DenseChromosome{T}}, len::Integer) = DenseChromosome{T}(Vector{T}(len))

type IntervalChromosome{Gene, Tk<:Real} <: Chromosome{Gene}
    data::SortedDict{Tk, Gene}
    n::Tk
end
call{T, Tk<:Real}(::Type{IntervalChromosome{T, Tk}}, len::Tk) =
    IntervalChromosome{T,Tk}(SortedDict(Dict{Tk,T}()), len)

getindex(chr::DenseChromosome, i::Integer) = chr.data[i]
function getindex(chr::IntervalChromosome, i::Real)
    (i < 0 || chr.n < i) && throw(BoundsError("Outside of a valid chromosomal range."))
    val = find(chr.data, i)
    if val == pastendtoken(chr.data)
        Nullable{typeof(i)}()
    else
        Nullable(deref_value(val))
    end
end
setindex!{Gene}(chr::DenseChromosome{Gene}, g::Gene, i::Integer) = chr.data[i] = g
setindex!{Gene}(chr::IntervalChromosome{Gene}, g::Gene, i::Real) = chr.data[i] = g

enumerate(chr::DenseChromosome) = enumerate(chr.data)

enumerate(chr::IntervalChromosome) = chr.data

length(chr::DenseChromosome) = length(chr.data)
length(chr::IntervalChromosome) = chr.n

function _writedata!{C<:Chromosome}(d::C, p::C)
    # write new data
    for (i, j) in enumerate(p)
        d.data[i] = j
    end
end

writedata!{C<:DenseChromosome}(d::C, p::C) = _writedata!(d, p)

function writedata!{C<:IntervalChromosome}(d::C, p::C)
    empty!(d.data)
    _writedata!(d, p)
    d.n = p.n
end

setoffspring!{C<:Chromosome}(d::C, p::C) = writedata!(d, p)

function setoffspring!{C<:DenseChromosome, I<:Integer}(d::C, ps::(C, C), selected::Int, recs::Vector{I})
    i0 = 1
    for recblock = 1:length(recs)
        for i = i0:recs[recblock]
            d.data[i] = inherit(ps[selected].data[i])
        end
        selected = 3 - selected
        i0 = recs[recblock] + 1
    end
    for i = i0:length(d)
        d.data[i] = inherit(ps[selected].data[i])
    end
end

function setoffspring!{C<:IntervalChromosome, I<:Real}(d::C, ps::(C, C), selected::Int, recs::Vector{I})
    empty!(d.data)
    # Assuming position is always positive.
    pos0 = -1.0
    iends = (pastendtoken(ps[1].data), pastendtoken(ps[2].data))
    is = [startof(ps[1].data), startof(ps[2].data)]
    for recsite in [recs; length(ps[1])]
        i = is[selected]
        while i != iends[selected]
            pos, val = deref(i)
            if pos > recsite
                # If the current position is past the next recombination site, switch to another parantal chromosome.

                # Cache the current state
                is[selected] = i
                selected = 3 - selected
                i = is[selected]
                pos0 = recsite
                break
            elseif pos <= pos0
                i = advance(i)
            else
                d.data[pos] = inherit(val)
                i = advance(i)
            end
        end
    end
end

inherit(i::Number) = i


text:
using RecyclePermanent
using RNGPool
using JLD2
using BenchmarkTools

include("../models/bpModel.jl")

@load "data/bpData.jld" ys θ
model = makeBPModel(θ, ys)

N = 1024*1024
xs = Vector{Int64Particle}(undef, N)
for i in 1:N
  xs[i] = Int64Particle()
end
gs = Vector{Float64}(undef, N)

function fillμ(μ!::F, xs::Vector{Int64Particle}) where F<:Function
  rng = getRNG()
  for i in 1:length(xs)
    @inbounds μ!(xs[i], rng)
  end
end

function fillG(G::F, xs::Vector{Int64Particle}, gs::Vector{Float64}) where
  F<:Function
  for i in 1:length(xs)
    @inbounds gs[i] = G(1, xs[i])
  end
end

@btime fillμ($model.μ!, $xs)
@btime fillG($model.G, $xs, $gs)


text:
function LaplaceExpansion(perm_mat, full_mat)

    function collect_sub_mat_perm(A)

        if length(size(A)) == 1
            v = [(A...)...]
            return [v[setdiff(1:end)] for i in 1:length(v)]
        else
            sub_mat = [remove_row_col(A, [], [i]) for i in 1:size(A)[2]]
            return [fast_glynn_perm(a) for a in sub_mat]
        end

    end

    res = []
    v_perms = collect_sub_mat_perm(perm_mat)

    for i in 1:size(full_mat)[1]
        rowA = full_mat[i,:]
        ans = 0
        for j in 1:length(v_perms)
            ans += rowA[j] * v_perms[j]
        end
        push!(res, ans)
    end

    return [abs(r)^2 for r in res]

end

function total_variation_distance(p, q)

    if length(p) > length(q)
        vcat(q, zeros(length(p)-length(q)))
    elseif length(q) > length(p)
        vcat(p, zeros(length(q)-length(p)))
    end

    return sum(abs(p[i]-q[i]) for i = 1:length(p))

end

function collect_sub_mat_perm(U)
    sub_mat = [remove_row_col(U, [], [i]) for i in 1:size(U)[2]]
    return [fast_glynn_perm(m) for m in sub_mat]
end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the Lattice-Boltzmann equations for the D2Q9 scheme

equations = LatticeBoltzmannEquations2D(Ma=0.1, Re=1000)

initial_condition = initial_condition_lid_driven_cavity
boundary_conditions = (
                       x_neg=boundary_condition_wall_noslip,
                       x_pos=boundary_condition_wall_noslip,
                       y_neg=boundary_condition_wall_noslip,
                       y_pos=boundary_condition_lid_driven_cavity,
                      )

solver = DGSEM(polydeg=5, surface_flux=flux_godunov)

coordinates_min = (0, 0)
coordinates_max = (1, 1)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=4,
                periodicity=false,
                n_cells_max=10_000,)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 400.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 1000
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_restart = SaveRestartCallback(interval=1000,
                                   save_final_restart=true)

save_solution = SaveSolutionCallback(interval=1000,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2macroscopic)

stepsize_callback = StepsizeCallback(cfl=1.0)

collision_callback = LBMCollisionCallback()

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_restart, save_solution,
                        stepsize_callback,
                        collision_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
using Test

using ChipSort
using SIMD

include("testutils.jl")

@testset "Transpose" begin

# types_to_test = [UInt8, UInt16, UInt32, UInt64, Float32, Float64]
types_to_test = [UInt16, Float64]

function test_transpose(T, V, J, K)
    data = randa(T, V*J*K)
    sol = permutedims(reshape(data, V, K, J), [1,3,2])[:]
    # seed=find_transpose_cycles(Val(J), Val(K))
    # transpose!(aa, Val(V), Val(J), Val(K), seed)
    transpose!(data, Val(V), Val(J), Val(K))
    @test sol == data
end

@testset for T in types_to_test
    for V in 2 .^ (1:2:5)
        for J in 2 .^ (1:2:5)
            for K in 2 .^ (2:4:8)
                test_transpose(T, V, J, K)
            end
        end
    end
end

function test_transpose_vecs(T, N, L)
    uu = ntuple(l->Vec(tuple([convert(T, mod(l*N+n-N, 2^(8*T.size-2))) for n in 1:N]...)), L)
    ttref = ntuple(l->Vec(tuple([convert(T, mod(n*N+l-N, 2^(8*T.size-2))) for n in 1:L]...)), N)
    tthat = transpose_vecs(uu...)
    @test all(all(j) for j in (tthat.==ttref))
end

@testset for T in types_to_test
    for N in 2 .^ (1:2:5)
        for L in 2 .^ (1:2:5)
            test_transpose_vecs(T, N, L)
        end
    end
end

end


text:
####### Payoffs definition

## Engines
include("engines/general_european_engine.jl")
include("engines/general_american_engine.jl")
include("engines/general_path_dependent_engine.jl")
include("engines/general_bermudan_engine.jl")


text:
include("featurescaling.jl")

X = [randn(100); 5+10*randn(100); 50+2*randn(100)]

X = reshape(X,(100,3))'

X_stand = standardization(X)
mean(X_stand,2)
std(X_stand,2)

standardization!(X)
mean(X,2)
std(X,2)

X = [randn(100); 5+10*randn(100); 50+2*randn(100)]

X = reshape(X,(100,3))

X_stand = standardization(X)
mean(X_stand,2)
std(X_stand,2)

standardization!(X)
mean(X,1)
std(X,1)

X = 5+10*randn(100)

X_stand = standardization(X)
mean(X_stand)
std(X_stand)

standardization!(X)
mean(X)
std(X)

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


text:
using DataFrames
using CSV
using Plots
#%%
# run(`./pi2`)
file1 = open("error_empirical.dat")
file2 = open("error_theo.dat")
df1=CSV.read(file1, header=false, delim=" ")
df2=CSV.read(file2, header=false, delim=" ")

plot(df1[1],df1[2], label="empirical", title="Monte Carlo Pi error")
plot!(df2[1],df2[2], label="theoretical")
xaxis!("N", xscale= :log10)
yaxis!("error", ylims=(1e-5,10), yscale= :log10)

max_stars_repo_name: aurelio-amerio/corso-cpp


text:
export iauJdcalf
"""
Julian Date to Gregorian Calendar, expressed in a form convenient
for formatting messages:  rounded to a specified precision.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   ndp       int      number of decimal places of days in fraction
   dj1,dj2   double   dj1+dj2 = Julian Date (Note 1)

Returned:
   iymdf     int[4]   year, month, day, fraction in Gregorian
                     calendar

Returned (function value):
            int      status:
                        -1 = date out of range
                        0 = OK
                        +1 = NDP not 0-9 (interpreted as 0)

Notes:

   1. The Julian Date is apportioned in any convenient way between
   the arguments dj1 and dj2.  For example, JD=2450123.7 could
   be expressed in any of these ways, among others:

            dj1            dj2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   2. In early eras the conversion is from the "Proleptic Gregorian
   Calendar";  no account is taken of the date(s) of adoption of
   the Gregorian Calendar, nor is the AD/BC numbering convention
   observed.

   3. Refer to the function iauJd2cal.

   4. NDP should be 4 or less if internal overflows are to be
   avoided on machines which use 16-bit integers.

Called:
   iauJd2cal    JD to Gregorian calendar

Reference:

   Explanatory Supplement to the Astronomical Almanac,
   P. Kenneth Seidelmann (ed), University Science Books (1992),
   Section 12.92 (p604).

This revision:  2016 December 2

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauJdcalf(int ndp, double dj1, double dj2, int iymdf[4])
function iauJdcalf(ndp::Int, dj1::Real, dj2::Real)
   iymdf = zeros(Int32, 4)

   status = ccall((:iauJdcalf, libsofa_c), Cint,
       (Cint, Cdouble, Cdouble, Ptr{Cint}),
       convert(Int32, ndp), convert(Float64, dj1), convert(Float64, dj2), iymdf)

   if status != 0
       @warn "iauJdcalf return non-zero exit code"
   end

   return status, SVector{4}(iymdf)
end


text:
# using Pkg
# Pkg.activate(joinpath(@__DIR__, ".."))
# Pkg.instantiate()

module LabClimateDashboard

using Dash
using DataFrames, CSV
using Dates
using Statistics
import Sockets: getipaddr

include("utils.jl")
include("layout.jl")

# default lower limit of the time axis of the plot
const START_DATE = now() - Dates.Day(3)
# refresh interval in s -> the Plot will update in this interval
const INTERVAL = 5
# when did the last update of the plot happen -> we need this to evaluate what
# is new data and should be appended to the data already in the plot
const last_update = Ref{DateTime}(now()) # use Ref to make a global of fixed type
# Have the dataframe that holds all the values global
const df = Ref{DataFrame}()
global LOGDATA_PATH = ""
# initialize the app
const app = dash(external_stylesheets=["https://codepen.io/chriddyp/pen/bWLwgP.css"])
layout!(app)

include("callbacks.jl")

"""
    run(LOGDATA_PATH; port=8050, debug=false)

Start the server. `LOGDATA_PATH` should be the path the the folder where the
log data resides.
"""
function run(logdata_path; port=8050, debug=false)
    global LOGDATA_PATH = logdata_path

    println("Loading in the log data...")
    load_logdata!(df, LOGDATA_PATH)
    sort!(df[])

    isempty(df[]) && @warn "Did not find any data!"

    ip4 = getipaddr()
    println("""
    Starting the server...
    In your browser visit http://127.0.0.1:$port/ from this machine to get to
    the dashboard.

    From other machines in the same network as this machine visit
    http://$ip4:$port/ to view the dashboard.

    Terminate the process by pressing 'Ctrl+c'
    """)

    println("Starting the server...")
    run_server(app, "0.0.0.0"; debug)

    nothing
end

end


text:
using DataStructures

function readmatrix(filename)
  lines = readlines(filename)
  nums = [[parse(Int, x) for x in line] for line in lines]
  M = zeros(Int, (length(nums),length(nums[1])))
  for r in 1:length(nums)
    M[r,:] = nums[r]
  end
  M
end

const M = readmatrix("inputs/11-15/11.txt")

neighbors((x,y)) = [(x+i,y+j) for i in -1:1, j in -1:1 if (i,j) != (0,0) && 0 < x+i <= 10 && 0 < y+j <= 10]

function day11(M)
  for s in 1:500
    f = 0
    M.+=1
    S = Stack{Tuple{Int64, Int64}}()
    for i in 1:10, j in 1:10
      if M[i,j] > 9
        push!(S, (i,j))
      end
    end
    while !isempty(S)
      (i,j) = pop!(S)
      f += 1
      for (ni, nj) in neighbors((i,j))
        M[ni,nj] += 1
        if M[ni,nj] == 10
          push!(S, (ni,nj))
        end
      end
    end
    f == 100 && return s
    M = map(x -> x > 9 ? 0 : x, M)
  end
  return -1
end

using BenchmarkTools

@show day11(M)
@btime day11(M) setup=(M=readmatrix("inputs/11-15/11.txt"))