text:
using Assignment05
using Test

@testset "Assignment05" begin

testpath = normpath(joinpath(@__DIR__, "..", "data", "test.txt"))
wordspath = normpath(joinpath(@__DIR__, "..", "data", "words.txt"))

@testset "Question 1" begin
    @test isreverse("racecar", "racecar")
    @test isreverse("GATA", "ATAG")
    @test !isreverse("GATA", "GTAG")
    @test !isreverse("1234", "1234")
end

@testset "Question 2" begin
    @test isreversecomplement("aaatttcg", "cgaaattt")
    @test !isreversecomplement("C", "A")
    @test_throws ErrorException isreversecomplement("TX", "AG")
    @test_throws ErrorException isreversecomplement("G", "CC")
end

@testset "Question 3" begin
    @test reverse_complement("AAATTT") ==  "AAATTT"
    @test reverse_complement("GCAT") == "ATGC"
    @test reverse_complement("TTGGG") == "CCCAA"

end

@testset "Question 4" begin
    q4a = string.(split(Assignment05.first15))[1:end-1]
    push!(q4a, wordspath)
    f15 = split(read(Cmd(q4a), String), '\n', keepempty=false)
    alllines = readlines(wordspath)
    @test f15 == alllines[1:15]
    q4b = string.(split(Assignment05.last15))[1:end-1]
    push!(q4b, wordspath)
    l15 = split(read(Cmd(q4b), String), '\n', keepempty=false)
    @test l15 == alllines[end-14:end]
end

@testset "Question 5" begin
    alllines = uppercase.(readlines(testpath))
    filter!(l-> all(c-> in(c, "ATGC"), l), alllines)
    @test length(alllines) == 3
    @test alllines == find_dna(testpath)
end

end # Assignment05


text:
# Switch to debug:
# macro dbgcore(x); esc(:(println(_dbg($x)))); end
macro dbgcore(x); end

# Contents:

# Here are the rough steps performed by g(x) where g=grad(f):
# 1. g is called with the same inputs as f.
# 2. g calls forward_pass which boxes x in a Value type and calls f(Value(x)).
# 3. If a primitive operator inside f gets a Value input, it records its action and returns a Value output.
# 4. g calls backward_pass which returns the gradient df/dx.

# And some background info:
# 5. How recording is done.
# 6. How new primitives and their gradients are defined.
# 7. How higher order gradients work.

# Details:

# 1. g is called with the same inputs as f.
# 1.1 g supports both regular and keyword args.
# 1.2 only one of the regular args is the gradient target, specified by the argnum argument of grad (defaults to 1).
# 1.3 in a typical model f would take parameters, return loss, with data kept in global variables.
# 1.4 to support multiple parameters, they can be grouped in a single arg using Array, Dict, or Tuple.

"""

grad(fun, argnum=1) -> gradfun

* fun: X->Y
* gradfun: X->dX

Returns a function which computes the gradient of `fun` with respect
to positional argument number `argnum`. The function `fun` should be
scalar-valued. The returned function `gradfun` takes the same
arguments as `fun`, but returns the gradient instead. The gradient has
the same type and size as the target argument which can be a Number,
Array, Tuple, or Dict.

"""
function grad(fun::Function, argnum::Int=1)
    @dbgcore((:grad,fun,argnum))
    function gradfun(args...; kwargs...)
        backward_pass(forward_pass(fun, args, kwargs, argnum)...)
    end
    return gradfun
end

# 2. g calls forward_pass which boxes argnum'th arg x in a Value type and calls f(Value(x))
# 2.1 f must be defined generically to accept Value arguments.
# 2.2 before the call a new Tape (tape1) is created (this is the only place a new tape is created)
# 2.3 we box the argument if it is not already (as may happen in higher order derivatives) and add it to the new tape
# 2.4 f is called with the boxed argument
# 2.5 the downstream operations on x are recorded on all its tapes
# 2.6 the output of f (end_value) could be a boxed Value or a regular value (if it does not depend on x)

function forward_pass(fun, args, kwargs, argnum)
    @dbgcore((:forw, argnum, fun, args..., kwargs...))
    tape = Tape()
    arg_wrt = args[argnum]
    if isa(arg_wrt,Value)
        Node(arg_wrt, tape)
        start_value = arg_wrt
    else
        start_value = Value(arg_wrt,tape)
    end
    args = Any[args...] # to make args writeable
    args[argnum] = start_value
    @dbgcore((:fcall, fun, args..., kwargs...))
    end_value = fun(args...; kwargs...)
    return start_value, end_value, tape
end

# forward_pass type: ((N(X)->N(Y)/Y),N(X)/X,K,I)->(N(X),N(Y)/Y,T)
# forward_pass deps: Tape, Value, Node

# 3. If a primitive operator inside f gets a Value input, it records its action and returns a Value output.

# 3.1 We implement this by dispatching f to r=recorder(f) if any of
# its arguments is a Value.  r unboxes the arguments, calls f, boxes
# and returns the result, recording the result and its dependencies on
# each boxed argument.

# We only need one recorder per function, but recorder(f) may be
# called many times for different methods.  To avoid duplication we
# hold recorders in a hash.

let fdict=ObjectIdDict()
global recorder

"""
recorder(fun) returns rfun, a recording version of fun.  It is used to
define primitive operations. rfun is defined with a generic signature
r(args...; kwargs...) and is intended to catch all invocations that
have at least one Value argument.
"""
function recorder(f)
r = get(fdict,f,0)
r != 0 && return r

function rfun(args...; kwargs...)
    @dbgcore((:call, f, args..., kwargs...))
    argvals = unbox(args)
    result = f(argvals...; kwargs...)
    for argnum = 1:length(args)
        arg = args[argnum]
        isa(arg,Value) || continue
        for t=1:length(arg.tapes)
            tape = arg.tapes[t]
            iscomplete(tape) && continue
            parent = arg.nodes[t]
            if !isa(result,Value)
                result = Value(result, tape; func=f, args=args, kwargs=kwargs)
                rnode = result.nodes[1]
            else
                s = findeq(result.tapes, tape)
                if s > 0
                    rnode = result.nodes[s]
                else
                    rnode = Node(result, tape)
                end
            end
            rnode.parents[argnum] = parent
        end
    end
    @dbgcore((:retn, f, result, args..., kwargs...))
    return result
end # function rfun
return (fdict[f] = rfun)
end # function recorder
end # let fdict

# recorder deps: Value, Node, iscomplete, findeq

"getval(x) unboxes x if it is a Value, otherwise returns x."
getval(x) = (if isa(x, Value); x.value; else; x; end)  # we never create Value(Value).

# this is much faster than map(getval,args)
function unbox(args)
    vals = Array(Any,length(args))
    for i=1:length(args)
        ai = args[i]
        if isa(ai,Value)
            vals[i] = ai.value
        else
            vals[i] = ai
        end
    end
    return vals
end

# findfirst uses == which is inefficient for tapes, so we define findeq with ===
function findeq(A,v)
    for i=1:length(A)
        if A[i] === v
            return i
        end
    end
    return 0
end

# 4. g calls backward_pass which returns the gradient df/dx.

# 4.1 backward_pass is called with start_value: Value(x), end_value:
# f(Value(x)) (which may or may not be a boxed Value), and the tape
# created by the corresponding forward_pass.  Note that Value(x) may
# point to more tapes in case of a higher order gradient.  It returns
# the gradient wrt the start_value.

function backward_pass(start_value, end_value, tape)
    @dbgcore((:back,:start,start_value,:end,end_value,:tape,tape,tape...))

# 4.2 If end_value is not a Value on the given tape, we return zero
# df/fx if x is a bits type, `nothing` otherwise.  end_value may not
# be a Value if the output of f does not depend on x.

    if !isa(end_value, Value) || 0==(tapeidx=findeq(end_value.tapes, tape))
        @dbgcore("Output seems independent of input. Returning zero gradient.")
        if isa(start_value,Number); return zero(start_value); else; return nothing; end
    end

    if !isa(end_value.value, Number)
        error("grad requires a scalar-valued function, got $(end_value.value)")
    end

# 4.3 backward_pass resets all node gradients except for the scalar
# output Value whose gradient is set to 1.0.

    end_value.nodes[tapeidx].outgrad = 1.0

# We need to complete!(tape) to prevent recording during backward_pass.

    complete!(tape)

# 4.4 the tape is read in reverse and for each node with a non-zero
# outgrad its ingrads are computed using the gradient methods.

    for node in tape[end-1:-1:1]  # note the end-1 because we pushed an eot marker
        node.outgrad == nothing && continue
        for i=1:length(node.parents)
            isassigned(node.parents,i) || continue
            parent = node.parents[i]
            v = node.value
            @dbgcore((:back,v.func,Grad{i},node.outgrad,v.value,v.args...,v.kwargs...))
            og = v.func(Grad{i},node.outgrad,v.value,v.args...;v.kwargs...)
            @dbgcore((:sum1,parent.outgrad,og))
            parent.outgrad = sum_outgrads(parent.outgrad, og)
            @dbgcore((:sum2,parent.outgrad))
        end
    end

# 4.5 tape[1].outgrad is returned.  How do we know this is the
# correct gradient df/dx?  Only x and its descendents are marked as
# Values and recorded on the tape. In the beginning the only non-empty
# outgrad is the one for the end_value.  Since the end_value is a
# boxed Value (otherwise we returned 0/nothing), it must depend on
# input x.  The input is the first thing recorded on tape by
# forward_pass, thus will be the last thing whose gradient is seen.
# If there are Values influenced by x but do not influence the
# end_value, their outgrad will remain empty, thus only the necessary
# gradients are computed.

    return tape[1].outgrad
end

# back deps: complete!, sum_outgrads

# 5. How recording is done.

# 5.1 Value: g=grad(f) calls forward_pass which calls f with one
# argument boxed in a Value type.  The primitives inside f call their
# recorder methods when one of their arguments is a Value.  The results
# of these recorder methods are also boxed in Value types, which will
# cause downstream primitives to be recorded as well.  The final
# output of f, if not independent of the input, will thus be a Value.

# 5.2 Node: Each result Value created by a primitive keeps track of
# the function and the arguments that created the Value.  Because a
# Value may need to be recorded in multiple tapes for higher order
# derivatives (see Sec. 7) these dependencies are kept in a separate
# data structure called a Node.  The parents field of a Node is an
# array that points to the Nodes of the arguments.  The gradient wrt
# the result is kept in the outgrad field.  outgrad is an array
# because a node can have multiple descendents each of which will push
# a gradient to outgrad to be summed.

# 5.3 Tape: When forward_pass is done, we have the computation graph
# (dependency tree) of the result recorded in Nodes.  However we also
# need the time order in which these Nodes were created for the
# backward_pass.  The gradient functions of all the children of a node
# need to be called before its own gradient function.  For example if
# z depends on x and y, and y depends on x, we want to compute the
# gradients in z-y-x order.  If we do it in z-x-y order, the gradient
# function of x will be called before its descendent y.  Thus we keep
# a Tape which is an array of Nodes in the order they were created.

if !isdefined(:Node)
type Node
    value
    outgrad
    parents::Vector{Node}
    Node(v) = new(v, nothing, Array(Node,length(v.args)))
end #type
end #if

typealias Tape Vector{Node}

if !isdefined(:Value)
type Value{T}
    value::T
    func::Function
    args::Tuple
    kwargs::Vector
    tapes::Vector{Tape}
    nodes::Vector{Node}
end # type Value{T}
end # if

function Value(value, tape::Tape=Tape(); func=rand, args=(), kwargs=[])
    self = Value(value,func,args,kwargs,Tape[tape],Array(Node,1))
    node = Node(self)
    push!(tape,node)
    self.nodes[1] = node
    return self
end

function Node(v::Value, t::Tape) # assumes v is not already in t
    n = Node(v)
    push!(t, n)
    push!(v.nodes, n)
    push!(v.tapes, t)
    return n
end

# Primitives with Value arguments may be called during the
# backward_pass. We do not want those primitives being recorded any
# more (at least on the tape created by the corresponding
# forward_pass, see Sec 7 for details).  We stop the recording on a
# Tape by calling its complete! method.

if !isdefined(:iscomplete)
let eot = Node(Value(nothing))
    global iscomplete, complete!
    iscomplete(a::Tape)=(!isempty(a) && a[end]===eot)
    complete!(a::Tape)=push!(a,eot)
end # let
end # if

# 6. How new primitives and their gradients are defined.

# 6.1 Primitives

# AutoGrad primitives record their actions when they are called with
# some arguments boxed in Values (However, see 6.3 for
# undifferentiable primitives). Julia supports multiple dispatch,
# i.e. a single function can have multiple methods with different arg
# types.  AutoGrad supports multiple dispatch for primitives and
# gradients, i.e. only some of the methods of a function can be
# defined as primitives and have gradients.  Calls to a particular
# method where some arguments are boxed in Values are directed to the
# recorder function. The following example makes `sin(x::Number)` a
# primitive, but says nothing about e.g. `sin(x::Array)`.

#     local sin_r = recorder(sin)
#     sin{T<:Number}(x::Value{T}) = sin_r(x)

# With multiple arguments, things get a bit more complicated.  There
# is no easy way to say "at least one argument is a Value" in Julia.
# So one must define methods for all 2^N-1 combinations for
# boxed/unboxed arguments to be safe.  This example makes
# hypot(x1::Array,x2::Array) a primitive:

#     local hypot_r = recorder(hypot)
#     hypot{T<:Array,S<:Array}(x1::Value{T},x2::Value{S})=hypot_r(x1,x2)
#     hypot{T<:Array,S<:Array}(x1::Value{T},x2::S)=hypot_r(x1,x2)
#     hypot{T<:Array,S<:Array}(x1::T,x2::Value{S})=hypot_r(x1,x2)

# I wrote the @primitive macro in util.jl to automate this process.
# One restriction is the inability to target parametric methods such
# as `f{T<:Number}(AbstractArray{T})`.  Julia does not support
# `f{T<:Number,A<:AbstractArray{T}}(Value{A})` yet.

# One could also choose to be lazy and just say:

#     hypot(x...) = hypot_r(x...)

# This would send any argument combination not covered by regular
# hypot methods to the recorder function, which presumably includes
# calls with boxed arguments.  This is dangerous for several reasons:
# (1) the Julia base may contain a typeless method (e.g. it does for
# `vcat`) we are overwriting. (2) this catches boxed calls to hypot
# methods we may not support yet.  So generally I would not recommend
# it.

# 6.2 Gradients

if !isdefined(:Grad)
"Grad{N} creates a type used by AutoGrad to represent the gradient wrt N'th arg."
immutable Grad{N}; end
end

# In AutoGrad, gradients are defined using gradient methods that have
# the following signature:

#     f(Grad{i},dy,y,x...) => dx[i]

# Here `f` is the name of original function, Grad{i} is a Type
# constant that specifies the gradient wrt the i'th argument, `dy` is
# the gradient wrt the output `y`, and `x...` are the input arguments.
# In this case `f` was originally called with `f(x...)` and returned
# `y`.  Somebody handed us the gradient `dy` wrt the output and
# `f(Grad{i},...)` above is going to give us the gradient `dx[i]` wrt
# the i'th argument.

# Note that type declarations on the x's can be used to specialize
# this gradient to any method of the function `f`.  Here is the
# gradient for `sin`:

# `sin{T<:Number}(::Type{Grad{1}}, dy, y, x::Value{T})=dy*cos(x)

# For the second example a different gradient method is needed for
# each argument:

# `hypot{T<:Array,S<:Array}(::Type{Grad{1}},dy,y,x1::Value{T},x2::Value{S})=(dy.*x1./y)`
# `hypot{T<:Array,S<:Array}(::Type{Grad{2}},dy,y,x1::Value{T},x2::Value{S})=(dy.*x2./y)`

# And of course we need four more definitions for the other
# boxed/unboxed argument combinations, which the @primitive macro
# generates automatically.

# Finally, there are three cases of zero gradients that need to be
# handled:

# 6.3 Undifferentiable functions

# Piecewise constant functions such as `sign`, and non-numeric
# functions such as `size` are not differentiable wrt any of their
# arguments.  Unlike primitives, these functions do not need to record
# their action or return a boxed Value.  They can just unbox their
# arguments and return an unboxed value:

# `size(a::Value,i...)=size(a.value,i...)`

# The @zerograd macro defined in util.jl can be used to automate this.

# 6.4 Undifferentiable wrt unboxed arguments

# Methods such as `sum(a::Array,i::Int)` are only differentiable wrt
# some of their arguments (here `a` but not `i`).  These methods must
# record when their differentiable argument(s) are boxed and return
# boxed values.  If we are certain that their undifferentiable
# arguments are never going to be boxed, we can leave their gradients
# undefined:

#     local sum_r = recorder(sum)
#     sum{T<:Array}(a::Value{T},i::Int)=sum_r(a,i)
#     sum{T<:Array}(::Type{Grad{1}},dy,y,a::Value{T},i::Int)=dy.+zeros(a)

# 6.5 Undifferentiable wrt boxed arguments

# Finally, in the rare cases when an undifferentiable argument can be
# boxed, its gradient must be defined and must return `nothing`.  The
# utility function `ungetindex` in intefaces.jl which uses its first
# argument's shape as a template is one example of this rare class.

# 6.6 sum_outgrads

sum_outgrads(a::Number, b::Number)=a+b
sum_outgrads(a::Tuple, b::Tuple)=tuple([sum_outgrads(x,y) for (x,y) in zip(a,b)]...)
sum_outgrads(a::Associative, b::Associative) = (z=similar(a); for d in (a,b), (k,v) in d; z[k]=v+get(z,k,0); end; z)
sum_outgrads{T}(a::AbstractArray{T},b::AbstractArray{T})=(if isbits(T); (a+b); else; [sum_outgrads(x,y) for (x,y) in zip(a,b)]; end)
# sum_outgrads needs to be a primitive for higher order gradients:
sum_outgrads_r = recorder(sum_outgrads)
sum_outgrads(a::Value,b::Value)=sum_outgrads_r(a,b)
sum_outgrads(a::Value,b)=sum_outgrads_r(a,b)
sum_outgrads(a,b::Value)=sum_outgrads_r(a,b)
sum_outgrads{N}(::Type{Grad{N}},dy,y,x1,x2)=dy
# we use `nothing` to indicate zero gradients
sum_outgrads(::Void,::Void)=nothing
sum_outgrads(a::Value,::Void)=a   # to avoid ambiguity
sum_outgrads(::Void,a::Value)=a   # to avoid ambiguity
sum_outgrads(a,::Void)=a
sum_outgrads(::Void,a)=a

# 7. How higher order gradients work.

# Say g=grad(f) and h=grad(g) and we call h(x).
# h(x) calls forward_pass(g,x)
# merge_tapes in forward_pass(g,x) is a noop because x is not a Value.
# forward_pass(g,x) wraps x in v1=Value(x,t1:n1) with tape t1 and node n1 and calls g(v1).
# g(v1) calls forward_pass(f,v1), which creates v2=Value(x,t2:n2)
# merge_tapes in forward_pass(f,v1) creates v3=Value(x,[t1:n31,t2:n32]) with parents n31->n1, n32->n2.
# forward_pass(f,v1) calls f(v3)
# primitives in f(v3) push their result Nodes on both [t1,t2] and record parents on each tape separately.
# f(v3) returns v4=Value(y,[t1:n41,t2:n42]).
# g(v1) calls backward_pass(f)(v2,v4,t2).
# backward_pass(f) calls complete!(t2) and starts processing the nodes on t2 in reverse.
# the nodes on t2 only point to other nodes in t2, so backward_pass(f) fills outgrads on t2.
# backward_pass(f) calls gradient methods of the recorded primitives in f(v3).
# the operations of gradient methods are recorded only on t1, that is why we need iscomplete(t2) once we start backward_pass on t2.
# backward_pass(f) returns v5=Value(df/dx,t1:n5) which becomes the output of forward_pass(g,x)
# h(x) calls backward_pass(g)(v1,v5,t1).
# backward_pass(g) calls gradient methods recorded in t1.
# even though some inputs are Values again, nothing gets recorded and all primitives return values because t1 is complete.
# backward_pass(g) returns a regular value which becomes the output of h(x).

ext: jl


text:
import Base.isless

struct PixelKey{CT}
    val::CT
    time_step::Int
end
isless(a::PixelKey{T}, b::PixelKey{T}) where {T} = (a.val < b.val) || (a.val == b.val && a.time_step < b.time_step)

"""
```
segments                = watershed(img, markers)
```
Segments the image using watershed transform. Each basin formed by watershed transform corresponds to a segment.
If you are using image local minimas as markers, consider using [`hmin_transform`](@ref) to avoid oversegmentation.

Parameters:
-    img            = input grayscale image
-    markers        = An array (same size as img) with each region's marker assigned a index starting from 1. Zero means not a marker.
                      If two markers have the same index, their regions will be merged into a single region.
                      If you have markers as a boolean array, use `label_components`.

"""
function watershed(img::AbstractArray{T, N}, markers::AbstractArray{S,N}) where {T<:Images.NumberLike, S<:Integer, N}

    if axes(img) != axes(markers)
        error("image size doesn't match marker image size")
    end

    segments = copy(markers)
    pq = PriorityQueue{CartesianIndex{N}, PixelKey{T}}()
    time_step = 0

    R = CartesianIndices(axes(img))
    Istart, Iend = first(R), last(R)
    for i in R
        if markers[i] > 0
            for j in CartesianIndices(_colon(max(Istart,i-one(i)), min(i+one(i),Iend)))
                if segments[j] == 0
                    segments[j] = markers[i]
                    enqueue!(pq, j, PixelKey(img[i], time_step))
                    time_step += 1
                end
            end
        end
    end

    while !isempty(pq)
        current = dequeue!(pq)
        segments_current = segments[current]
        img_current = img[current]
        for j in CartesianIndices(_colon(max(Istart,current-one(current)), min(current+one(current),Iend)))
            if segments[j] == 0
                segments[j] = segments_current
                enqueue!(pq, j, PixelKey(img_current, time_step))
                time_step += 1
            end
        end
    end

    TM = meantype(T)
    num_segments        = Int(maximum(segments))
    labels              = Array(1:num_segments)
    region_means        = Dict{Int, TM}()
    region_pix_count    = Dict{Int, Int}()

    for i in R
        region_pix_count[segments[i]] = get(region_pix_count, segments[i], 0) + 1
        region_means[segments[i]] = get(region_means, segments[i], zero(TM)) + (img[i] - get(region_means, segments[i], zero(TM)))/region_pix_count[segments[i]]
    end
    return SegmentedImage(segments, labels, region_means, region_pix_count)
end

"""
```
out = hmin_transform(img, h)
```
Suppresses all minima in grayscale image whose depth is less than h.

H-minima transform is defined as the reconstruction by erosion of (img + h) by img. See Morphological image analysis by Soille pg 170-172.
"""
function hmin_transform(img::Array{T, N}, h::Real) where {T<:Images.NumberLike, N}
    out = img.+h
    while true
        temp = max.(img, erode(out))
        if temp == out
            break
        end
        out = temp
    end
    return out
end


text:
import Base.Broadcast: BroadcastStyle, DefaultArrayStyle, Style

# This is a `BroadcastStyle` for AbstractAbstractDimArray's
# It preserves the dimension names.
# `S` should be the `BroadcastStyle` of the wrapped type.
# Copied from NamedDims.jl (thanks @oxinabox).
struct DimensionalStyle{S <: BroadcastStyle} <: AbstractArrayStyle{Any} end
DimensionalStyle(::S) where {S} = DimensionalStyle{S}()
DimensionalStyle(::S, ::Val{N}) where {S,N} = DimensionalStyle(S(Val(N)))
DimensionalStyle(::Val{N}) where N = DimensionalStyle{DefaultArrayStyle{N}}()
function DimensionalStyle(a::BroadcastStyle, b::BroadcastStyle)
    inner_style = BroadcastStyle(a, b)
    # if the inner style is Unknown then so is the outer style
    if inner_style isa Unknown
        return Unknown()
    else
        return DimensionalStyle(inner_style)
    end
end

function BroadcastStyle(::Type{<:AbstractDimArray{T,N,D,A}}) where {T,N,D,A}
    inner_style = typeof(BroadcastStyle(A))
    return DimensionalStyle{inner_style}()
end

BroadcastStyle(::DimensionalStyle, ::Base.Broadcast.Unknown) = Unknown()
BroadcastStyle(::Base.Broadcast.Unknown, ::DimensionalStyle) = Unknown()
BroadcastStyle(::DimensionalStyle{A}, ::DimensionalStyle{B}) where {A, B} = DimensionalStyle(A(), B())
BroadcastStyle(::DimensionalStyle{A}, b::Style) where {A} = DimensionalStyle(A(), b)
BroadcastStyle(a::Style, ::DimensionalStyle{B}) where {B} = DimensionalStyle(a, B())
BroadcastStyle(::DimensionalStyle{A}, b::Style{Tuple}) where {A} = DimensionalStyle(A(), b)
BroadcastStyle(a::Style{Tuple}, ::DimensionalStyle{B}) where {B} = DimensionalStyle(a, B())
# We need to implement copy because if the wrapper array type does not
# support setindex then the `similar` based default method will not work
function Broadcast.copy(bc::Broadcasted{DimensionalStyle{S}}) where S
    _dims = _broadcasted_dims(bc)
    A = _firstdimarray(bc)
    data = copy(_unwrap_broadcasted(bc))
    return if A isa Nothing || _dims isa Nothing
        data
    else
        rebuild(A, data, _dims, refdims(A), Symbol(""))
    end
end

function Base.copyto!(dest::AbstractArray, bc::Broadcasted{DimensionalStyle{S}}) where S
    _dims = comparedims(dims(dest), _broadcasted_dims(bc))
    copyto!(dest, _unwrap_broadcasted(bc))
    A = _firstdimarray(bc)
    if A isa Nothing || _dims isa Nothing
        dest
    else
        rebuild(A, dest, _dims, refdims(A))
    end
end
function Base.copyto!(dest::AbstractDimArray, bc::Broadcasted{DimensionalStyle{S}}) where S
    _dims = comparedims(dims(dest), _broadcasted_dims(bc))
    copyto!(parent(dest), _unwrap_broadcasted(bc))
    A = _firstdimarray(bc)
    if A isa Nothing || _dims isa Nothing
        dest
    else
        rebuild(A, parent(dest), _dims, refdims(A))
    end
end

function Base.similar(bc::Broadcast.Broadcasted{DimensionalStyle{S}}, ::Type{T}) where {S,T}
    A = _firstdimarray(bc)
    rebuildsliced(A, similar(_unwrap_broadcasted(bc), T, axes(bc)...), axes(bc), Symbol(""))
end

# Recursively unwraps `AbstractDimArray`s and `DimensionalStyle`s.
# replacing the `AbstractDimArray`s with the wrapped array,
# and `DimensionalStyle` with the wrapped `BroadcastStyle`.
function _unwrap_broadcasted(bc::Broadcasted{DimensionalStyle{S}}) where S
    innerargs = map(_unwrap_broadcasted, bc.args)
    return Broadcasted{S}(bc.f, innerargs)
end
_unwrap_broadcasted(x) = x
_unwrap_broadcasted(nda::AbstractDimArray) = parent(nda)

# Get the first dimensional array in the broadcast
_firstdimarray(x::Broadcasted) = _firstdimarray(x.args)
_firstdimarray(x::Tuple{<:AbstractDimArray,Vararg}) = x[1]
_fistdimarray(ext::Base.Broadcast.Extruded) = _firstdimarray(ext.x)
function _firstdimarray(x::Tuple{<:Broadcasted,Vararg})
    found = _firstdimarray(x[1])
    if found isa Nothing
        _firstdimarray(tail(x))
    else
        found
    end
end
_firstdimarray(x::Tuple) = _firstdimarray(tail(x))
_firstdimarray(x::Tuple{}) = nothing

# Make sure all arrays have the same dims, and return them
_broadcasted_dims(bc::Broadcasted) = _broadcasted_dims(bc.args...)
_broadcasted_dims(a, bs...) = comparedims(_broadcasted_dims(a), _broadcasted_dims(bs...))
_broadcasted_dims(a::AbstractDimArray) = dims(a)
_broadcasted_dims(a) = nothing


text:
module ExtremeFloats
export ExtremeFloat

struct ExtremeFloat <: AbstractFloat
    negative::Bool
    mantissa::UInt128
    exponent::Int64

    @inline ExtremeFloat(negative, mantissa, exponent) = new(negative, mantissa, exponent)
end

function normalize(ef::ExtremeFloat)
    tz = trailing_zeros(ef.mantissa)
    ExtremeFloat(ef.negative, ef.mantissa >> tz, ef.exponent + tz)
end

Base.Float64(ef::ExtremeFloat) = ifelse(ef.negative, -1, 1) * Float64(ef.mantissa) * 2.0 ^ ef.exponent

function ExtremeFloat(f::Float64)
    if iszero(f)
        return ExtremeFloat(false, UInt128(0), Int64(0))
    end

    negative = false
    if f < zero(f)
        f *= -1
        negative = true
    end

    mantissa = round(UInt128, significand(f) * 2.0 ^ 52)
    exponent = Base.exponent(f) - 52
    ExtremeFloat(negative, mantissa, exponent) |> normalize
end

function ExtremeFloat(i::Integer)
    negative = false
    if i < zero(i)
        negative = true
        i *= -1
    end

    ExtremeFloat(negative, UInt128(i), Int64(0))
end

function Base.:*(ef1::ExtremeFloat, ef2::ExtremeFloat)
    exponent = ef1.exponent + ef2.exponent
    negative = ef1.negative ⊻ ef2.negative

    mantissa1 = ef1.mantissa
    mantissa2 = ef2.mantissa

    # shift mantissas such that their product has at most one leading zero
    num_digits = 128
    remaining_zeros = leading_zeros(mantissa1) + leading_zeros(mantissa2) - num_digits
    shift1 = remaining_zeros ÷ 2
    shift2 = remaining_zeros - shift1
    mantissa1 <<= shift1
    mantissa2 <<= shift2
    mantissa = mantissa1 * mantissa2
    exponent -= shift1 + shift2

    ExtremeFloat(negative, mantissa, exponent) |> normalize
end

function Base.:*(other::Number, ef::ExtremeFloat)
    ExtremeFloat(other) * ef
end

function Base.:+(ef1::ExtremeFloat, ef2::ExtremeFloat)
    if iszero(ef1.mantissa)
        return ef2
    elseif iszero(ef2.mantissa)
        return ef1
    end

    ef1, ef2 = if ef1.exponent > ef2.exponent
        (ef1, ef2)
    else
        (ef2, ef1)
    end
    exponent_diff = ef1.exponent - ef2.exponent
    shift1 = min(leading_zeros(ef1.mantissa) - 1, exponent_diff)
    shift2 = max(0, exponent_diff - shift1)

    mantissa1 = ef1.mantissa << shift1
    mantissa2 = ef2.mantissa >> shift2

    negative = if mantissa1 > mantissa2
        ef1.negative
    else
        ef2.negative
    end

    mantissa =
        ifelse(ef1.negative, -1, 1) * mantissa1 +
        ifelse(ef2.negative, -1, 1) * mantissa2
    if negative
        mantissa *= -1
    end
    exponent = ef1.exponent - shift1

    ExtremeFloat(negative, mantissa, exponent) |> normalize
end

function Base.:-(ef1::ExtremeFloat, ef2::ExtremeFloat)
    ef1 + ExtremeFloat(!ef2.negative, ef2.mantissa, ef2.exponent)
end

function Base.inv(ef::ExtremeFloat)
    mantissa_inv = ef.mantissa |> inv |> ExtremeFloat
    pow2 = ExtremeFloat(ef.negative, UInt128(1), -ef.exponent) |> normalize
    mantissa_inv * pow2
end

const LN2 = ExtremeFloat(log(2.))

function Base.log(ef::ExtremeFloat)
    ExtremeFloat(log(ef.mantissa)) + ExtremeFloat(ef.exponent) * LN2
end

const E = ℯ |> Float64 |> ExtremeFloat
const C1 = ExtremeFloat(362880)
const C2 = ExtremeFloat(362880)
const C3 = ExtremeFloat(181440)
const C4 = ExtremeFloat(60480)
const C5 = ExtremeFloat(15120)
const C6 = ExtremeFloat(3024)
const C7 = ExtremeFloat(504)
const C8 = ExtremeFloat(72)
const C9 = ExtremeFloat(9)
const C0 = ExtremeFloat(2.75573192e-6)

function Base.exp(ef::ExtremeFloat)
    non_frac = ef.mantissa << ef.exponent
    exp_non_frac = if iszero(non_frac)
        ExtremeFloat(1)
    else
        E ^ non_frac
    end

    frac = ef - ExtremeFloat(ef.negative, non_frac, 0)
    # from https://stackoverflow.com/a/10552567
    exp_frac = (C1+frac*(C2+frac*(C3+frac*(C4+frac*(C5+frac*(C6+frac*(C7+frac*(C8+frac*(C9+frac)))))))))*C0

    result = exp_non_frac * exp_frac
    if ef.negative
        result = inv(result)
    end

    result
end

end # module


text:
#QE calls these flags
struct QEFlagInfo{T}
    name::Symbol
    # default::Union{T,Nothing,Symbol}  #check again v0.7 Some
    description::String
end
QEFlagInfo() = QEFlagInfo{Nothing}(:error, "")
Base.eltype(x::QEFlagInfo{T}) where {T} = T

abstract type AbstractBlockInfo end

struct QEControlBlockInfo <: AbstractBlockInfo
    name::Symbol
    flags::Vector{<:QEFlagInfo}
end
Base.in(f::Symbol, i::QEControlBlockInfo) = findfirst(x -> x.name == f, i.flags) !== nothing

#TODO rewrite this this is garbage
function qe_flaginfo(block::AbstractBlockInfo, variable_name::Symbol)
    varstr1 = string(variable_name)
    for var in block.flags
        varstr2 = string(var.name)
        if occursin(varstr2, varstr1) &&
           (length(varstr1) == length(varstr2) || length(varstr1) == length(varstr2) + 2)
            return var
        end
    end
    return QEFlagInfo()
end

struct QEDataBlockInfo <: AbstractBlockInfo
    name                :: Symbol
    description         :: String
    options             :: Vector{Symbol}
    options_description :: String
    flags               :: Vector{<:QEFlagInfo}
end

struct QECalculationInfo
    exec::String
    control::Vector{QEControlBlockInfo}
    data::Vector{QEDataBlockInfo}
end

function allflags(info::QECalculationInfo)
    return vcat([[i.flags for i in info.control]; [i.flags for i in info.data]]...)
end

include(joinpath(DFC.DEPS_DIR, "qeflags.jl"))
const QECalculationInfos = _QEINPUTINFOS()
push!(QECalculationInfos,
      QECalculationInfo("pw2wannier90.x",
                        [QEControlBlockInfo(:inputpp,
                                            [QEFlagInfo{String}(:outdir,
                                                                "location of temporary output files"),
                                             QEFlagInfo{String}(:prefix,
                                                                "pwscf filename prefix"),
                                             QEFlagInfo{String}(:seedname,
                                                                "wannier90 calculation/output filename prefix"),
                                             QEFlagInfo{String}(:wan_mode,
                                                                "'standalone' or 'library'"),
                                             QEFlagInfo{String}(:spin_component,
                                                                "'none', 'up' or 'down'"),
                                             QEFlagInfo{Bool}(:write_spn,
                                                              "Write .spn matrix elements."),
                                             QEFlagInfo{Bool}(:write_mmn,
                                                              "compute M_mn matrix"),
                                             QEFlagInfo{Bool}(:write_amn,
                                                              "compute A_mn matrix"),
                                             QEFlagInfo{Bool}(:write_unk,
                                                              "write wavefunctions to file"),
                                             QEFlagInfo{Bool}(:write_uHu,
                                                              "write the hamiltonian elements between different k-values"),
                                             QEFlagInfo{Bool}(:wvfn_formatted,
                                                              "formatted or unformatted output for wavefunctions"),
                                             QEFlagInfo{Bool}(:reduce_unk,
                                                              "output wavefunctions on a coarse grid to save memory")])],
                        QEDataBlockInfo[]))

function qe_calculation_info(calculation::Calculation{QE})
    return getfirst(x -> occursin(x.exec, calculation.exec.exec), QECalculationInfos)
end
function qe_calculation_info(exec::AbstractString)
    return getfirst(x -> occursin(x.exec, exec), QECalculationInfos)
end
qe_calculation_flags(exec::AbstractString) = allflags(qe_calculation_info(exec))

function qe_flaginfo(calculation_info::QECalculationInfo, variable_name::Symbol)
    for block in vcat(calculation_info.control, calculation_info.data)
        var = qe_flaginfo(block, variable_name)
        if eltype(var) != Nothing
            return var
        end
    end
    return QEFlagInfo()
end

function qe_flaginfo(variable_name::Symbol)
    for info in QECalculationInfos
        var = qe_flaginfo(info, variable_name)
        if eltype(var) != Nothing
            return var
        end
    end
    return QEFlagInfo()
end

function qe_block_variable(calculation_info::QECalculationInfo, variable_name)
    for block in vcat(calculation_info.control, calculation_info.data)
        var = qe_flaginfo(block, variable_name)
        if eltype(var) != Nothing
            return block, var
        end
    end
    return :error, QEFlagInfo()
end

function qe_flaginfo(exec::Exec, varname)
    for calculation_info in QECalculationInfos
        if occursin(calculation_info.exec, exec.exec)
            return qe_flaginfo(calculation_info, varname)
        end
    end
    return QEFlagInfo()
end

function qe_block_info(block_name::Symbol)
    for calculation_info in QECalculationInfos
        for block in [calculation_info.control; calculation_info.data]
            if block.name == block_name
                return block
            end
        end
    end
end

function qe_all_block_flags(calculation::Calculation{QE}, block_name)
    return getfirst(x -> x.name == block, qe_calculation_info(calculation).control).flags
end
function qe_all_block_flags(exec::AbstractString, block_name)
    return getfirst(x -> x.name == block_name, qe_calculation_info(exec).control).flags
end

function qe_block_variable(exec::AbstractString, flagname)
    for calculation_info in QECalculationInfos
        if occursin(calculation_info.exec, exec)
            return qe_block_variable(calculation_info, flagname)
        end
    end
    return :error, QEFlagInfo()
end

function qe_exec(calculation::Calculation{QE})
    if !(calculation.exec.exec ∈ QE_EXECS)
        error("Calculation $calculation does not have a valid QE executable, please set it first.")
    end
    return calculation.exec
end

function qe_block_variable(calculation::Calculation, flagname)
    return qe_block_variable(qe_exec(calculation).exec, flagname)
end

function flagtype(calculation::Calculation{QE}, flag)
    return eltype(qe_flaginfo(qe_exec(calculation), flag))
end
flagtype(::Type{QE}, exec, flag) = eltype(qe_flaginfo(exec, flag))

isbands(c::Calculation{QE})   = get(c, :calculation, nothing) == "bands"
isnscf(c::Calculation{QE})    = get(c, :calculation, nothing) == "nscf"
isscf(c::Calculation{QE})     = get(c, :calculation, nothing) == "scf"
isvcrelax(c::Calculation{QE}) = get(c, :calculation, nothing) == "vc-relax"
isrelax(c::Calculation{QE})   = get(c, :calculation, nothing) == "relax"
isprojwfc(c::Calculation{QE}) = c.exec.exec == "projwfc.x"
ishp(c::Calculation{QE}) = c.exec.exec == "hp.x"

function ispw(c::Calculation{QE})
    return isbands(c) || isnscf(c) || isscf(c) || isvcrelax(c) || isrelax(c)
end

issoc(c::Calculation{QE}) = get(c, :lspinorb, false)

function ismagnetic(c::Calculation{QE})
    return get(c, :nspin, 0.0) > 0.0 || get(c, :total_magnetization, 0.0)
end

function outfiles(c::Calculation{QE})
    files = [c.outfile]
    for (is, fuzzies) in zip(("projwfc.x", "hp.x"), (("pdos",), ("Hubbard_parameters",)))
        if c.exec.exec == is
            append!(files, fuzzies)
        end
    end
    return unique(files)
end

ψ_cutoff_flag(::Calculation{QE}) = :ecutwfc
ρ_cutoff_flag(::Calculation{QE}) = :ecutrho

function kgrid(na, nb, nc, ::Type{QE})
    return reshape([(a, b, c, 1 / (na * nb * nc))
                    for a in collect(range(0; stop = 1, length = na + 1))[1:end-1],
                        b in collect(range(0; stop = 1, length = nb + 1))[1:end-1],
                        c in collect(range(0; stop = 1, length = nc + 1))[1:end-1]],
                   (na * nb * nc))
end

function set_kpoints!(c::Calculation{QE}, k_grid::NTuple{3,Int}; print = true) #nscf
    print && !isnscf(c) && (@warn "Expected calculation to be 'nscf'.\nGot $c.")
    d = data(c, :k_points)
    if d !== nothing
        d.data = kgrid(k_grid..., c)
        d.option = :crystal
    else
        push!(c.data, InputData(:k_points, :crystal, kgrid(k_grid..., c)))
    end
    prod(k_grid) > 100 && set_flags!(c, :verbosity => "high"; print = print)
    return c
end

function set_kpoints!(c::Calculation{QE}, k_grid::NTuple{6,Int}; print = true) #scf
    print &&
        !(isscf(c) || isvcrelax(c) || isrelax(c)) &&
        (@warn "Expected calculation to be scf, vc-relax, relax.\nGot $calc.")
    d = data(c, :k_points)
    if d !== nothing
        d.data = [k_grid...]
        d.option = :automatic
    else
        push!(c.data, InputData(:k_points, :automatic, [k_grid...]))
    end
    prod(k_grid[1:3]) > 100 && set_flags!(c, :verbosity => "high"; print = print)
    return c
end

function set_kpoints!(c::Calculation{QE}, k_grid::Vector{<:NTuple{4}}; print = true,
                      k_option = :crystal_b)
    print &&
        isbands(c) != "bands" &&
        (@warn "Expected calculation to be bands, got $(c[:calculation]).")
    @assert in(k_option, [:tpiba_b, :crystal_b, :tpiba_c, :crystal_c]) error("Only $([:tpiba_b, :crystal_b, :tpiba_c, :crystal_c]...) are allowed as a k_option, got $k_option.")
    if k_option in [:tpiba_c, :crystal_c]
        @assert length(k_grid) == 3 error("If $([:tpiba_c, :crystal_c]...) is selected the length of the k_points needs to be 3, got length: $(length(k_grid)).")
    end
    num_k = 0.0
    for k in k_grid
        num_k += k[4]
    end
    if num_k > 100.0
        set_flags!(c, :verbosity => "high"; print = print)
        if print
            @info "Verbosity is set to high because num_kpoints > 100,\n
                       otherwise bands won't get printed."
        end
    end
    d = data(c, :k_points)
    if d !== nothing
        data(c, :k_points).data = k_grid
        data(c, :k_points).option = k_option
    else
        push!(c.data, InputData(:k_points, k_option, k_grid))
    end
    return c
end

"""
    gencalc_scf(template::Calculation{QE}, kpoints::NTuple{6, Int}, newflags...; name="scf")

Uses the information from the template and `supplied` kpoints to generate an scf calculation.
Extra flags can be supplied which will be set for the generated calculation.
"""
function gencalc_scf(template::Calculation{QE}, kpoints::NTuple{6,Int}, newflags...;
                     name = "scf")
    return calculation_from_kpoints(template, name, kpoints, :calculation => "scf",
                                    newflags...)
end

"""
    gencalc_vcrelax(template::Calculation{QE}, kpoints::NTuple{6, Int}, newflags...; name="scf")

Uses the information from the template and supplied `kpoints` to generate a vcrelax calculation.
Extra flags can be supplied which will be set for the generated calculation.
"""
function gencalc_vcrelax(template::Calculation{QE}, kpoints::NTuple{6,Int}, newflags...;
                         name = "vcrelax")
    return calculation_from_kpoints(template, name, kpoints, :calculation => "vc-relax",
                                    newflags...)
end

"""
    gencalc_bands(template::Calculation{QE}, kpoints::Vector{NTuple{4}}, newflags...; name="bands")

Uses the information from the template and supplied `kpoints` to generate a bands calculation.
Extra flags can be supplied which will be set for the generated calculation.
"""
function gencalc_bands(template::Calculation{QE}, kpoints::Vector{<:NTuple{4}}, newflags...;
                       name = "bands")
    return calculation_from_kpoints(template, name, kpoints, :calculation => "bands",
                                    newflags...)
end

"""
    gencalc_nscf(template::Calculation{QE}, kpoints::NTuple{3, Int}, newflags...; name="nscf")

Uses the information from the template and supplied `kpoints` to generate an nscf calculation.
Extra flags can be supplied which will be set for the generated calculation.
"""
function gencalc_nscf(template::Calculation{QE}, kpoints::NTuple{3,Int}, newflags...;
                      name = "nscf")
    return calculation_from_kpoints(template, name, kpoints, :calculation => "nscf",
                                    newflags...)
end

"""
    gencalc_projwfc(template::Calculation{QE}, Emin, Emax, DeltaE, newflags...; name="projwfc")

Uses the information from the template and supplied `kpoints` to generate a `projwfc.x` calculation.
Extra flags can be supplied which will be set for the generated calculation.
"""
function gencalc_projwfc(template::Calculation{QE}, Emin, Emax, DeltaE, extraflags...;
                         name = "projwfc")
    occflag = get(template, :occupations, "fixed")
    ngauss  = 0
    if occflag == "smearing"
        smearingflag = get(template, :smearing,"gaussian")
        if smearingflag ∈ ("methfessel-paxton", "m-p", "mp")
            ngauss = 1
        elseif smearingflag ∈ ("marzari-vanderbilt", "cold", "m-v", "mv")
            ngauss = -1
        elseif smearingflag ∈ ("fermi-dirac", "f-d", "fd")
            ngauss = -99
        end
    end
    tdegaussflag = get(template, :degauss, nothing)
    degauss = tdegaussflag !== nothing ? tdegaussflag : 0.0
    exec = Exec(dir = template.exec.dir, exec = "projwfc.x", modules = template.exec.modules, flags = template.exec.flags)
    empty!(exec.flags)
    out = Calculation(deepcopy(template); name = name, exec = exec, data = InputData[])
    set_name!(out, "projwfc")
    empty!(out.flags)
    set_flags!(out, :Emin => Emin, :Emax => Emax, :DeltaE => DeltaE, :ngauss => ngauss,
               :degauss => degauss; print = false)
    set_flags!(out, extraflags...)
    return out
end


text:
g = 9.81
#Insira o valor de L
L = 1000 #Em metros
T = 2pi*sqrt(L/g)

max_stars_repo_name: PedroMoraesFranco/fisica-computacional

max_issues_repo_name: PedroMoraesFranco/fisica-computacional

max_forks_repo_name: PedroMoraesFranco/fisica-computacional


text:
function make_file(path::AbstractString, size)
    open(path, "w") do io
        written = 0
        while written < size
            batch = min(2*1024*1024, size - written)
            written += write(io, zeros(UInt8, batch))
        end
    end
end

@testset "get_disk_freespace" begin
    mktempdir() do dir
        # Note that other disk activity can cause this test to be a little noisy,
        # so we do approximate math and use large file sizes in an attempt to be
        # more robust against random fluctuations.
        freespace = get_disk_freespace(dir)
        make_file(joinpath(dir, "big"), 50*1024*1024)
        @test abs(Int64(freespace - get_disk_freespace(dir)) - 50*1024*1024) < 20*1024*1024
        rm(joinpath(dir, "big"); force=true)
        @test abs(Int64(freespace - get_disk_freespace(dir))) < 20*1024*1024
    end
end


text:
"""
The `HeatmapCanvas` is also Unicode-based.
It has a half the resolution of the `BlockCanvas`.
This canvas effectively turns every character
into two pixels (top and bottom).
"""
struct HeatmapCanvas <: LookupCanvas
    grid::Array{UInt8,2}
    colors::Array{ColorType,2}
    pixel_width::Int
    pixel_height::Int
    origin_x::Float64
    origin_y::Float64
    width::Float64
    height::Float64
end

const HALF_BLOCK = '▄'
@inline x_pixel_per_char(::Type{HeatmapCanvas}) = 1
@inline y_pixel_per_char(::Type{HeatmapCanvas}) = 2

@inline lookup_encode(c::HeatmapCanvas) = [0 0; 1 1]
@inline lookup_decode(c::HeatmapCanvas) = [HALF_BLOCK; HALF_BLOCK]

@inline nrows(c::HeatmapCanvas) = div(size(grid(c), 2) + 1, 2)

function HeatmapCanvas(args...; kwargs...)
    CreateLookupCanvas(HeatmapCanvas, args...; min_char_width=1, min_char_height=1, kwargs...)
end

_toCrayon(c) = c === nothing ? 0 : (c isa Unsigned ? Int(c) : c)

function printrow(io::IO, c::HeatmapCanvas, row::Int)
    0 < row <= nrows(c) || throw(ArgumentError("Argument row out of bounds: $row"))
    y = 2*row
    # extend the plot upwards by half a row
    if isodd(size(grid(c), 2))
        y -= 1
    end
    iscolor = get(io, :color, false)
    for x in 1:ncols(c)
        if iscolor
            fgcol = _toCrayon(c.colors[x, y])
            if (y - 1) > 0
                bgcol = _toCrayon(c.colors[x, y - 1])
                print(io, Crayon(foreground=fgcol, background=bgcol), HALF_BLOCK)
            # for odd numbers of rows, only print the foreground for the top row
            else
                print(io, Crayon(foreground=fgcol), HALF_BLOCK)
            end
        else
            print(io, HALF_BLOCK)
        end
    end
    if iscolor
        print(io, Crayon(reset=true))
    end
end

function printcolorbarrow(io::IO, c::HeatmapCanvas, row::Int, colormap::Any, border::Symbol, lim, plot_padding, zlabel)
    b = bordermap[border]
    min_z = lim[1]
    max_z = lim[2]
    if row == 1
        # print top border and maximum z value
        print_color(:light_black, io, b[:tl], b[:t], b[:t], b[:tr])
        max_z_str = isinteger(max_z) ? max_z : float_round_log10(max_z)
        print(io, plot_padding)
        print_color(:light_black, io, max_z_str)
    elseif row == nrows(c)
        # print bottom border and minimum z value
        print_color(:light_black, io, b[:bl], b[:b], b[:b], b[:br])
        min_z_str = isinteger(min_z) ? min_z : float_round_log10(min_z)
        print(io, plot_padding)
        print_color(:light_black, io, min_z_str)
    else
        # print gradient
        print_color(:light_black, io, b[:l])
        # if min and max are the same, single color
        if min_z == max_z
            bgcol = colormap(1, 1, 1)
            fgcol = bgcol
        # otherwise, blend from min to max
        else
            n = 2*(nrows(c) - 2)
            r = row - 2
            bgcol = colormap(n - (2*r),     1, n)
            fgcol = colormap(n - (2*r + 1), 1, n)
        end
        print(io, Crayon(foreground=fgcol, background=bgcol), HALF_BLOCK)
        print(io, HALF_BLOCK)
        print(io, Crayon(reset=true))
        print_color(:light_black, io, b[:r])

        # print z label
        if row == div(nrows(c), 2) + 1
            print(io, plot_padding)
            print(io, zlabel)
        end
    end
end


text:
using PlotRecipes
using Base.Test

# write your own tests here
@test 1 == 1

# TODO: actual tests

# -----------------------------------------
# marginalhist

# using Distributions
# n = 1000
# x = rand(Gamma(2), n)
# y = -0.5x + randn(n)
# marginalhist(x, y)

# -----------------------------------------
# portfolio composition map

# # fake data
# tickers = ["IBM", "Google", "Apple", "Intel"]
# N = 10
# D = length(tickers)
# weights = rand(N,D)
# weights ./= sum(weights, 2)
# returns = sort!((1:N) + D*randn(N))

# # plot it
# portfoliocomposition(weights, returns, labels = tickers')

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


text:
using FileIO
using ColorTypes
using ColorTypes.FixedPointNumbers
using HTTP
using Pkg
using Test

@testset "Integration" begin
    if haskey(ENV, "CI") && Base.VERSION >= v"1.3"  # ImageIO can't be loaded on Julia < 1.3, so putting it as test dep is problematic
        Pkg.add("ImageIO")
        img = rand(RGB{N0f8}, 50, 50)
        io = IOBuffer()
        save(Stream{format"PNG"}(io), img)
        buf = take!(io)
        io2 = IOBuffer(buf)
        img2 = load(io2)
        @test img2 == img

        uri = HTTP.URI("https://upload.wikimedia.org/wikipedia/commons/thumb/b/b3/Wikipedia-logo-v2-en.svg/135px-Wikipedia-logo-v2-en.svg.png")
        @test isa(load(uri), Matrix)
    end
end


text:
using Revise
using AdFem
using PyPlot

function ricker(;dt = 0.005, f0 = 3.0)
    nw = 2/(f0*dt)
    nc = floor(Int, nw/2)
    t = dt*collect(-nc:1:nc)
    b = (π*f0*t).^2
    w = @. (1 - 2b)*exp(-b)
end

n = 50
mmesh = Mesh(n, n, 1/n)
NT = 400
Δt = 2/NT

xy = gauss_nodes(mmesh)
xg, yg = xy[:,1], xy[:,2]
xy = fem_nodes(mmesh)
x, y = xy[:,1], xy[:,2]

β = zeros(mmesh.ndof)
β_g = zeros(get_ngauss(mmesh))
βprime = zeros(get_ngauss(mmesh))
λ = constant(ones(get_ngauss(mmesh)))
μ = constant(0.5ones(get_ngauss(mmesh)))

nv = [ones(get_ngauss(mmesh)) zeros(get_ngauss(mmesh))]

R = 5000
pf = n->R*n^3
pfprime = n->3*R*n^2

for i = 1:get_ngauss(mmesh)
    n = 0.0

    if yg[i]<0.2
        nv[i,:] = [0.0;-1.0]
        n = max(n, abs(0.2-yg[i]))
    elseif yg[i]>0.8
        nv[i,:] = [0.0;1.0]
        n = max(abs(0.8-yg[i]), n)
    end

    if xg[i]<0.2
        nv[i,:] = [-1.0;0.0]
        n = max(n, abs(0.2-xg[i]))
    elseif xg[i]>0.8
        nv[i,:] = [1.0;0.0]
        n = max(n, abs(0.8-xg[i]))
    end

    if xg[i]<0.2 && yg[i]<0.2
        nv[i,:] = [xg[i]-0.2 yg[i]-0.2]
        ns = norm(nv[i,:])
        nv[i,:] = nv[i,:]/ns
    end
    if xg[i]<0.2 && yg[i]>0.8
        nv[i,:] = [xg[i]-0.2 yg[i]-0.8]
        ns = norm(nv[i,:])
        nv[i,:] = nv[i,:]/ns
    end
    if xg[i]>0.8 && yg[i]>0.8
        nv[i,:] = [xg[i]-0.8 yg[i]-0.8]
        ns = norm(nv[i,:])
        nv[i,:] = nv[i,:]/ns
    end
    if xg[i]>0.8 && yg[i]<0.2
        nv[i,:] = [xg[i]-0.8 yg[i]-0.2]
        ns = norm(nv[i,:])
        nv[i,:] = nv[i,:]/ns
    end

    if (0.2<xg[i]<0.8) && (0.2<yg[i]<0.8)
        continue
    end
    β_g[i] = pf(n)
    βprime[i] = pfprime(n)
end

for i = 1:mmesh.ndof
    n = 0.0

    if x[i]<0.2
        n = abs(0.2-x[i])
    elseif x[i]>0.8
        n = abs(0.8-x[i])
    end

    if y[i]<0.2
        n = max(n, abs(0.2-y[i]))
    elseif y[i]>0.8
        n = max(abs(0.8-y[i]), n)
    end

    if (0.2<x[i]<0.8) && (0.2<y[i]<0.8)
        continue
    end

    β[i] = pf(n)
end

nv = constant(nv)
β_g = constant(β_g)
βprime = constant(βprime)
β = constant(β)
β = [β;β]


text:
"""
    CalculateSampling(in)

Calculate the sampling operator of an n-dimension input. The output has the
same size as the input.

"""
function CalculateSampling(in)

cutoff = 1e-10
	itrace = 1

	n=size(in)
	in=reshape(in,n[1],:)
        wd = zeros(Float32,size(in))
	n2=size(in,2)
	for itrace = 1 : n2
		a = sum(in[:,itrace].*in[:,itrace])
		if (a > cutoff)
			wd[:,itrace] .= 1.
		end
	end
	wd=reshape(wd,n)
	return wd;

end


text:
"Meta data to attach to ω::Ω"
const Tags{K, V} = NamedTuple{K, V}

"Does tag type contain `t`, forall t in `tags`?"
hastags(::Type{Tags{K, V}}, tags::Symbol...) where {K, V} = all([t in K for t in tags])

"""Combine tags for a particular tag type

`combinetags(::Type{Val{:tagtype}}, a, b)`

Tag types should implement this
"""
function combinetags end

"Sample space tagged with meta-data.  Enables `iid`, `replace`, `trackerror`"
struct TaggedΩ{I, TAGS <: Tags, ΩT} <: Ω{I}
  taggedω::ΩT
  tags::TAGS
  TaggedΩ(ω::ΩT, tags::TAGS) where {I, ΩT <: Ω{I}, TAGS} =
    new{I, TAGS, ΩT}(ω, tags)
end

hastags(::Type{TaggedΩ{I, TAGS, ΩT}}, tags...) where {I, TAGS, ΩT} = hastags(TAGS, tags...)

tag(ω, tag_::Tags) = TaggedΩ(ω, tag_)
tag(tω::TaggedΩ, tag_::Tags) = TaggedΩ(tω.taggedω, merge(combinetags, tag_, tω.tags))

"Returns ω1 with tags from ω2 transferred over"
transfertags(ω1::TaggedΩ, ω2::Ω) = ω1                       # ω1 has no tags to transfer
transfertags(ω1::TaggedΩ, ω2::TaggedΩ) = tag(ω1, ω2.tags)  # merge tags
transfertags(ω1::Ω, ω2::TaggedΩ) = tag(ω1, ω2.tags)  # merge tags
transfertags(ω1::Ω, ω2::Ω) = ω1                             # Neither have tags

# Pass-throughs (tω::TaggedΩ should work like its tω.taggedω, but preserve tags)

Base.values(tω::TaggedΩ) = values(tω.taggedω)
Base.keys(tω::TaggedΩ) = keys(tω.taggedω)

@inline Base.getindex(tω::TaggedΩ, i) = TaggedΩ(getindex(tω.taggedω, i), tω.tags)
# Base.rand(tω::TaggedΩ, args...) = rand(tω.taggedω, args...; rng = rng(tω))

# Just pass through rand to tω.taggedω (have to repeat due to ambiguities with Random)
@inline Base.rand(tω::TaggedΩ) = rand(tω.taggedω; rng = rng(tω))
@inline Base.rand(tω::TaggedΩ, dims::Dims) = rand(tω.taggedω, dims; rng = rng(tω))
@inline Base.rand(tω::TaggedΩ, dim::Integer, dims::Integer...) = rand(tω.taggedω, dim, dims...; rng = rng(tω))

@inline Base.rand(tω::TaggedΩ, ::Type{T}) where T = rand(tω.taggedω, T; rng = rng(tω))
@inline Base.rand(tω::TaggedΩ, ::Type{T}, dims::Dims) where T = rand(tω.taggedω, T, dims; rng = rng(tω))
@inline Base.rand(tω::TaggedΩ, ::Type{T}, dim::Integer, dims::Integer...) where T= rand(tω.taggedω, T, dim, dims...; rng = rng(tω))

@inline Base.rand(tω::TaggedΩ, arr::Array) = rand(tω.taggedω, arr; rng = rng(tω))
@inline Base.rand(tω::TaggedΩ, arr::Array, dims::Dims) = rand(tω.taggedω, arr, dims; rng = rng(tω))
@inline Base.rand(tω::TaggedΩ, arr::Array, dim::Integer, dims::Integer...) = rand(tω.taggedω, arr, dim, dims...; rng = rng(tω))

@inline Base.rand(tω::TaggedΩ, ur::UnitRange) = rand(tω.taggedω, ur; rng = rng(tω))
@inline Base.rand(tω::TaggedΩ, ur::UnitRange, dims::Dims) = rand(tω.taggedω, ur, dims; rng = rng(tω))
@inline Base.rand(tω::TaggedΩ, ur::UnitRange, dim::Integer, dims::Integer...) = rand(tω.taggedω, ur, dim, dims...; rng = rng(tω))

@inline Base.rand(tω::TaggedΩ, D::Distributions.Normal, args...)  =
  rand(tω.taggedω, D, args...; rng = rng(tω))

@inline parentω(tω::TaggedΩ) = TaggedΩ(parentω(tω), tω.tags)


text:
"""
Definition of `DynamicalSystem` as well as all integrators
used in the ecosystem of DynamicalSystems.jl.

Also contains predefined well-known systems.
"""
module DynamicalSystemsBase

using DelayEmbeddings
import SciMLBase
# Import some core functions that are extended a lot:
import SciMLBase: init, step!, isinplace, reinit!, u_modified!

include("deprecations.jl")

include("core/dynamicalsystem.jl")
include("core/create_dynamics.jl")
include("core/api_docstrings.jl")
include("core/discrete.jl")
include("core/continuous.jl")
include("core/pretty_printing.jl")

include("advanced_integrators/projected_integrator.jl")
include("advanced_integrators/stroboscopic_map.jl")

include("predefined/famous_systems.jl") # includes both discrete and continuous

export GeneralizedDynamicalSystem, DynamicalSystem
export ContinuousDynamicalSystem, DiscreteDynamicalSystem
export get_rule_for_print, isdiscretetime
export dimension, get_state, get_states
export set_parameter!, trajectory, jacobian
export integrator, tangent_integrator, parallel_integrator
export set_state!, get_state, get_deviations, set_deviations!, current_time
export SciMLBase, init, step!, isinplace, reinit!
export Systems
export SVector, SMatrix, @SVector, @SMatrix
export projected_integrator, stroboscopicmap

end


text:
# Translated from PyTorch version, see cnn_minst.py

using Avalon
import Avalon: accuracy
using MLDatasets

# include("../src/core.jl")
# __init__()

mutable struct Net
    conv1::Conv2d
    conv2::Conv2d
    fc1::Linear
    fc2::Linear
end

Net() = Net(
    Conv2d(1, 20, 5),
    Conv2d(20, 50, 5),
    Linear(4 * 4 * 50, 500),
    Linear(500, 10)
)

function (m::Net)(x::AbstractArray)
    x = maxpool2d(relu.(m.conv1(x)), (2, 2))
    x = maxpool2d(relu.(m.conv2(x)), (2, 2))
    # prod(size(x)[1:3])
    x = reshape(x, 4*4*50, :)
    x = relu.(m.fc1(x))
    x = logsoftmax(m.fc2(x))
    return x
end

function get_mnist_data(train::Bool; device=CPU())
    X, Y = train ? MNIST.traindata() : MNIST.testdata()
    X = convert(Array{Float64}, reshape(X, 28, 28, 1, :)) |> device
    # replace class label like "0" with its position like "1"
    Y = Y .+ 1 |> device
    return X, Y
end

function main()
    # choose device: if CUDA is available on the system, GPU() will be used, otherwise - CPU()
    device = best_available_device()
    # instantiate the model
    m = Net() |> device
    # load training data
    X_trn, Y_trn = get_mnist_data(true);
    # set loss function and optimizer, then fit the model
    loss_fn = NLLLoss()
    opt = SGD(1e-2; momentum=0)
    @time fit!(m, X_trn, Y_trn, loss_fn; n_epochs=10, opt=opt, batch_size=100, device=device)

    # load test data
    X_tst, Y_tst = get_mnist_data(false)
    # convert to device
    X_tst = X_tst |> device; Y_tst = Y_tst |> device
    # predict log probabilities and calculate accuracy
    Ŷ = m(X_tst)
    @info accuracy(Y_tst, Ŷ)
end


text:
module TestCentrality
  using Test
  using EcologicalNetworks

  A = UnipartiteNetwork(map(Bool, [0.0 1.0 0.0 0.0; 0.0 0.0 0.0 1.0; 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0]))

  c_1_1 = centrality_katz(A, a=0.1, k=1)
  c_2_1 = centrality_katz(A, a=0.2, k=1)
  c_1_2 = centrality_katz(A, a=0.1, k=2)

  e_c_1_1 = [0.0 2.0/3.0 0.0 1.0/3.0]
  e_c_2_1 = [0.0 2.0/3.0 0.0 1.0/3.0]
  e_c_1_2 = [0.0 0.625 0.0 0.374]

  for i = eachindex(species(A))
    si = species(A)[i]
    @test e_c_1_1[i] ≈ c_1_1[si] atol=0.01
    @test e_c_2_1[i] ≈ c_2_1[si] atol=0.01
    @test e_c_1_2[i] ≈ c_1_2[si] atol=0.01
  end

  # Test network

  N = UnipartiteNetwork([
    false true true false false;
    false false false true false;
    false false false true false;
    false false false false true;
    false false false false false
  ])

  @test centrality_degree(N)["s1"] ≈ 2/3
  @test centrality_degree(N)["s2"] ≈ 2/3
  @test centrality_degree(N)["s3"] ≈ 2/3
  @test centrality_degree(N)["s4"] ≈ 1.0
  @test centrality_degree(N)["s5"] ≈ 1/3

  @test centrality_closeness(N)["s1"] ≈ 4 / 7
  @test centrality_closeness(N)["s2"] ≈ 1.3333333333333
  @test centrality_closeness(N)["s3"] ≈ 1.3333333333333
  @test centrality_closeness(N)["s4"] ≈ 4 / 1
  @test centrality_closeness(N)["s5"] ≈ 0.0

end


text:
"""
    Quantizer{T <: Unsigned}

Apply value quantization by scaling it down from one type to another, which
requires less bits than the original type.
"""
struct Quantizer{T <: Unsigned}
    minval::Float32
    maxval::Float32
    scale::Float32
    zero_point::Float32
end

"""
    Quantizer{T}(minval, maxval)

Create a new `Quantizer` which can quantize values from the range
`[minval, maxval]`.
"""
Quantizer{T}(minval, maxval) where T <: Unsigned = begin
    scale = (maxval - minval) / typemax(T)
    zero_point = -minval / scale
    return Quantizer{T}(minval, maxval, scale, zero_point)
end

"""
    Quantizer{T}(data::Vector{Float32})

Create a new `Quantizer` for the given `data`.
"""
Quantizer{T}(data::Vector{Float32}) where T <: Unsigned = begin
    minval = minimum(data)
    maxval = maximum(data)
    return Quantizer{T}(minval, maxval)
end

"""
    quantize(q::Quantizer{T}, x) where T <: Unsigned

Quantize `x` with the Quantizer `q`.
"""
quantize(q::Quantizer{T}, x) where T <: Unsigned = round(T, x / q.scale + q.zero_point)

"""
    dequantize(q::Quantizer, x)

Dequantize `x` with the Quantizer `q`.
"""
dequantize(q::Quantizer, x) = q.scale * (x - q.zero_point)


text:
import Base: ==

# syntactic equality of systems
function ==(sys1::AbstractSystem, sys2::AbstractSystem)
    if typeof(sys1) != typeof(sys2)
        return false
    end
    for field in fieldnames(typeof(sys1))
        if getfield(sys1, field) != getfield(sys2, field)
            return false
        end
    end
    return true
end

# approximate equality of systems
function Base.isapprox(sys1::AbstractSystem, sys2::AbstractSystem; kwargs...)
    if typeof(sys1) != typeof(sys2)
        return false
    end
    for field in fieldnames(typeof(sys1))
        if !isapprox(getfield(sys1, field), getfield(sys2, field); kwargs...)
            return false
        end
    end
    return true
end

"""
    typename(system::AbstractSystem)

Returns the base type of `system` without parameter information.

### Input

- `system` -- `AbstractSystem`

### Output

The base type of `system`.

"""
function typename(system::AbstractSystem)
    return Base.typename(typeof(system)).wrapper
end

"""
    _complementary_type(system_type::Type{<:AbstractSystem})

Return the complementary type of a system type `system_type`.

### Input

- `system_type` -- type of `AbstractSystem`

### Ouput

Return complementary type of `system_type`.

### Notes

There are two main subclasses of abstract types: continuous types and discrete
types. A complementary type of `system_type` has the same fields as `system_type`
but belongs to the other subclass, e.g. for a `LinearContinuousSystem` which is
a subtype of `AbstractContinuousSystem` and has the field `:A`, the subtype of
`AbstractDiscreteSystem` with the field `:A`, i.e. `LinearDiscreteSystem`,
is returned.

To get the complementary type of system type, use
`_complementary_type(typename(system))`.
"""
@generated function _complementary_type(type::Type{<:AbstractSystem})
    system_type = type.parameters[1]
    type_string = string(system_type)
    if supertype(system_type) == AbstractDiscreteSystem
        type_string = replace(type_string, "Discrete"=>"Continuous")
    elseif supertype(system_type) == AbstractContinuousSystem
        type_string = replace(type_string, "Continuous"=>"Discrete")
    else
        error("$system_type <: $(supertype(system_type)) is neither discrete nor continuous")
    end
    return Meta.parse(type_string)
end


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

JLLWrappers.@generate_main_file_header("SuiteSparse_LBT")
JLLWrappers.@generate_main_file("SuiteSparse_LBT", UUID("8a7bcfce-81ee-5a07-a07c-c8520ac69669"))
end  # module SuiteSparse_LBT_jll


text:
"""
    function updates_!(p,abm)

Function that checks the variables in the model that are modified at each step in order to make appropiate copy vectors and add them to program.
"""
function updates_!(p::Program_,abm::Agent,space::SimulationSpace)
    syms = []
    for i in keys(abm.declaredUpdates)
        if !(emptyquote_(abm.declaredUpdates[i])) && !(i in ["UpdateInteraction","UpdateLocalInteraction"])
            s = symbols_(abm,abm.declaredUpdates[i])
            append!(syms,s[Bool.(s[:,"updated"] .+ s[:,"assigned"]),"Symbol"])
        end
    end

    unique!(syms)

    ## Assign updates of variable types
    for t in keys(abm.declaredSymbols)
        dict = Dict{Symbol,Int}()
        counter = 1
        for i in syms
            if i in abm.declaredSymbols[t]
                dict[i] = counter
                counter += 1
            end
        end
        p.update[t] = dict
    end

    ## Check equations and their update variables
    dict = Dict{Symbol,Int}()
    counter = 1
    if "Equation" in keys(abm.declaredUpdates)
        s = symbols_(abm,abm.declaredUpdates["Equation"]).Symbol
        for i in abm.declaredSymbols["Local"]
            ss = Meta.parse(string(EQUATIONSYMBOL,i))
            if ss in s
                dict[i] = counter
                counter += 1
                if !(i in keys(p.update["Local"]))
                    if isempty(p.update["Local"])
                        p.update["Local"][i] = 1
                    else
                        p.update["Local"][i] = maximum(values(p.update["Local"])) + 1
                    end
                end
            end
        end
    end
    p.update["Variables"] = dict

    ## Check division and their update variables
    dict = Dict{Symbol,Int}()
    counter = 1
    if "EventDivision" in keys(abm.declaredUpdates)
        s = symbols_(abm,abm.declaredUpdates["EventDivision"]).Symbol
        for place in ["Local","Identity"]
            for i in abm.declaredSymbols[place]
                for j in DIVISIONSYMBOLS
                    ss = Meta.parse(string(i,j))
                    if ss in s && !(i in keys(dict))
                        dict[i] = counter
                        counter += 1
                        if !(i in keys(p.update[place]))
                            if isempty(p.update[place])
                                p.update[place][i] = 1
                            else
                                p.update[place][i] = maximum(values(p.update[place])) + 1
                            end
                        end
                    end
                end
            end
        end
    end
    p.update["EventDivision"] = dict

    ## Check space and bound updates
    for i in space.box
        if i.s in keys(p.update["Local"])
            for j in keys(i.addSymbols)
                for k in i.addSymbols[j]
                    if !(k in keys(p.update["Local"]))
                        p.update["Local"][k] = maximum(values(p.update["Local"])) + 1
                    end
                end
            end
        end
    end

    ## Check equations and their update variables
    dict = Dict{Symbol,Int}()
    counter = 1
    if "UpdateMedium" in keys(abm.declaredUpdates)
        s = symbols_(abm,abm.declaredUpdates["UpdateMedium"]).Symbol
        for i in abm.declaredSymbols["Medium"]
            ss = Meta.parse(string(MEDIUMSYMBOL,i))
            if ss in s
                dict[i] = counter
                counter += 1
                if !(i in keys(p.update["Medium"]))
                    if isempty(p.update["Medium"])
                        p.update["Medium"][i] = 1
                    else
                        p.update["Medium"][i] = maximum(values(p.update["Medium"])) + 1
                    end
                end
            end
        end
    end
    p.update["Medium"] = dict

    return
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
using DrWatson
@quickactivate "EPExperiments"

using EPExperiments
using Random
using Logging
using LoggingExtras
using JLD2
using FileIO
using CSV
using DataFrames
using LinearAlgebra: I

const LOG_DIR_NAME = "posterior_predictive_mh_baseline"
const LOG_FILENAME = "out.log"

@with_kw struct PostPredMHBaseline
    experiment_name::String = "test"
    num_samples::Int = 100
    num_chains::Int = 4
    seed::Int = 1234
end

function post_pred_mh_baseline(experiment_config)
    @unpack seed, num_samples, num_chains = experiment_config

    result_folder = mkpath(joinpath(
        datadir(LOG_DIR_NAME),
        savename(experiment_config)
    ))

    logger = TeeLogger(
        ConsoleLogger(),
        FileLogger(joinpath(result_folder, LOG_FILENAME))
    )

    Random.seed!(seed)

    post_pred_turing = post_pred(POST_PRED_Y_OBSERVED).gamma2

    @time begin
    nuts_samples = sample(
        post_pred_turing,
        Turing.MH(Matrix{Float64}(I, POST_PRED_DIMENSION, POST_PRED_DIMENSION)),
        MCMCThreads(),
        num_samples,
        num_chains;
        discard_adapt=false
    )
    end

    @tagsave(joinpath(result_folder, "results.jld2"), Dict(
        "nuts_samples" => nuts_samples,
        "experiment_config" => experiment_config
    ))

    with_logger(logger) do
        @info "MCMC samples:" nuts_samples
    end
end

post_pred_mh_baseline(PostPredMHBaseline(
    experiment_name="mh_baseline",
    num_samples=10000,
    num_chains=28,
    seed=1390
))


text:
# ## Scatter plot
# ### Scatter plot 2d

#md # ![FILE_NAME.png](images/FILE_NAME.png)

using Random
Random.seed!(2018)

n = 50
x = rand(n)
y = rand(n)
ms = rand(50) * 30
scatter(x, y, markersize=ms)


text:
module SpaceAnalysis

include("utility.jl")
include("constants.jl")
include("signal.jl")
include("MVA.jl")
include("MFA.jl")

end


text:
function get_partition_dirs(table_name)
  meta_path = joinpath("data", table_name, "_meta")
  meta_file = open(meta_path, "r")
  meta_data = read(meta_file, String)
  # Could parse JSON but this is faster + easier
  dir_regex = r"\"dir\":\s*\"(.*)\""
  sort([m.captures[1] for m in eachmatch(dir_regex, meta_data)])
end

function get_column_type(col_ext)
  if col_ext == ".u64"
    UInt64
  else
    println("Unknown column type ", col_ext)
    exit(2)
  end
end

function get_type_column(col_type)
  if col_type == UInt32
    ".u32"
  else
    println("Unknown column type ", col_type)
    exit(3)
  end
end

function convert_column(partition_dirs, col_name, to_type)
  col_files = readdir(partition_dirs[1])
  col_basenames = filter(d -> startswith(d, col_name * "."), col_files)
  if length(col_basenames) > 1
    println("Found multiple columns ", col_basenames)
    exit(1)
  end
  (col_basename, col_ext) = splitext(col_basenames[1])
  for p in partition_dirs
    col_path = joinpath(p, col_basename * col_ext)
    println(col_path)
    io_in = open(col_path, "r")
    data = UInt8[]
    readbytes!(io_in, data, Inf)
    close(io_in)
    data = reinterpret(get_column_type(col_ext), data)
    data = convert(Array{to_type}, data)
    println(length(data))
    mv(col_path, col_path * ".bak")
    new_path = joinpath(p, col_basename * get_type_column(to_type))
    println(new_path)
    io_out = open(new_path, "w")
    write(io_out, data)
    close(io_out)
  end
end

function convert_table_column(table_name, col_name, to_type)
  dirs = get_partition_dirs(table_name)
  convert_column(dirs, col_name, to_type)
end

convert_table_column("agg1m", "volume", UInt32)

max_stars_repo_name: clickingbuttons/zdb

max_issues_repo_name: clickingbuttons/zdb

max_forks_repo_name: clickingbuttons/zdb


text:
using Documenter, FFTW

makedocs(
    modules = [FFTW],
    clean = false,
    format = :html,
    sitename = "FFTW.jl",
    pages = Any[
        "Home" => "index.md",
        "API" => "fft.md",
    ],
)

deploydocs(
    julia = "nightly",
    repo = "github.com/JuliaMath/FFTW.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
import Base: merge, hcat, vcat, map

###### merge ####################

function _merge_outer(::Type{IndexType}, ta1::TimeArray{T,N,D}, ta2::TimeArray{T,M,D}, padvalue, meta) where {IndexType,T,N,M,D}
    timestamps, new_idx1, new_idx2 = sorted_unique_merge(IndexType, timestamp(ta1), timestamp(ta2))
    vals = fill(convert(T, padvalue), (length(timestamps), length(colnames(ta1)) + length(colnames(ta2))))
    insertbyidx!(vals, values(ta1), new_idx1)
    insertbyidx!(vals, values(ta2), new_idx2, size(values(ta1), 2))
    TimeArray(timestamps, vals, [colnames(ta1); colnames(ta2)], meta; unchecked = true)
end

function merge(ta1::TimeArray{T,N,D}, ta2::TimeArray{T,M,D}, method::Symbol = :inner;
               colnames::Vector = Symbol[], meta = nothing,
               padvalue=NaN) where {T,N,M,D}

    if colnames isa Vector{<:AbstractString}
        @warn "`merge(...; colname::Vector{<:AbstractString})` is deprecated, " *
              "use `merge(...; colnames=Symbol.(colnames))` instead."
        colnames = Symbol.(colnames)
    end

    meta = if _meta(ta1) == _meta(ta2) && meta ≡ nothing
        _meta(ta1)
    elseif typeof(_meta(ta1)) <: AbstractString && typeof(_meta(ta2)) <: AbstractString && meta ≡ nothing
        string(_meta(ta1), "_", _meta(ta2))
    else
        meta
    end

    if method == :inner

        idx1, idx2 = overlap(timestamp(ta1), timestamp(ta2))
        vals = [values(ta1[idx1]) values(ta2[idx2])]
        ta = TimeArray(timestamp(ta1[idx1]), vals, [_colnames(ta1); _colnames(ta2)], meta; unchecked = true)

    elseif method == :left

        new_idx2, old_idx2 = overlap(timestamp(ta1), timestamp(ta2))
        right_vals = fill(convert(T, padvalue), (length(ta1), length(_colnames(ta2))))
        insertbyidx!(right_vals, values(ta2), new_idx2, old_idx2)
        ta = TimeArray(timestamp(ta1), [values(ta1) right_vals], [_colnames(ta1); _colnames(ta2)], meta; unchecked = true)

    elseif method == :right

        ta = merge(ta2, ta1, :left; padvalue = padvalue)
        ncol2 = length(_colnames(ta2))
        vals = [values(ta)[:, (ncol2+1):end] values(ta)[:, 1:ncol2]]
        ta = TimeArray(timestamp(ta), vals, [_colnames(ta1); _colnames(ta2)], meta; unchecked = true)

    elseif method == :outer

        ta = if (length(timestamp(ta1)) + length(timestamp(ta2))) > typemax(Int32)
            _merge_outer(Int64, ta1, ta2, padvalue, meta)
        else
            _merge_outer(Int32, ta1, ta2, padvalue, meta)
        end

    else
        throw(ArgumentError(
            "merge method must be one of :inner, :left, :right, :outer"))
    end

    return setcolnames!(ta, colnames)

end

# hcat ##########################

function hcat(x::TimeArray, y::TimeArray)
    tsx = timestamp(x)
    tsy = timestamp(y)

    if length(tsx) != length(tsx) || tsx != tsy
        throw(DimensionMismatch(
            "timestamps not consistent, please checkout `merge`."))
    end

    m = ifelse(meta(x) == meta(y), meta(x), nothing)

    TimeArray(tsx, [values(x) values(y)], [colnames(x); colnames(y)], m)
end

hcat(x::TimeArray, y::TimeArray, zs::Vararg{TimeArray}) =
    hcat(hcat(x, y), zs...)

# collapse ######################

function collapse(ta::TimeArray{T,N,D}, period::Function, timestamp::Function,
                  value::Function = timestamp) where {T,N,D}

    length(ta) == 0 && return ta

    ncols = length(colnames(ta))
    collapsed_tstamps = D[]
    collapsed_values = values(ta)[1:0, :]

    tstamp = _timestamp(ta)[1]
    mapped_tstamp = period(tstamp)
    cluster_startrow = 1

    for i in 1:length(ta)-1

        next_tstamp = _timestamp(ta)[i+1]
        next_mapped_tstamp = period(next_tstamp)

        if mapped_tstamp != next_mapped_tstamp
          push!(collapsed_tstamps, timestamp(_timestamp(ta)[cluster_startrow:i]))
          collapsed_values = [collapsed_values; T[value(values(ta)[cluster_startrow:i, j]) for j in 1:ncols]']
          cluster_startrow = i+1
        end #if

        tstamp = next_tstamp
        mapped_tstamp = next_mapped_tstamp

    end #for

    push!(collapsed_tstamps, timestamp(_timestamp(ta)[cluster_startrow:end]))
    collapsed_values = [collapsed_values; T[value(values(ta)[cluster_startrow:end, j]) for j in 1:ncols]']

    N == 1 && (collapsed_values = vec(collapsed_values))
    return TimeArray(collapsed_tstamps, collapsed_values, colnames(ta), meta(ta))

end

# vcat ######################

function vcat(TA::TimeArray...)
    # Check all meta fields are identical.
    prev_meta = meta(TA[1])
    for ta in TA
        if meta(ta) != prev_meta
            throw(ArgumentError("metadata doesn't match"))
        end
    end

    # Check column names are identical.
    prev_colnames = colnames(TA[1])
    for ta in TA
        if colnames(ta) != prev_colnames
            throw(ArgumentError("column names don't match"))
        end
    end

    # Concatenate the contents.
    ts  = vcat([timestamp(ta) for ta in TA]...)
    val = vcat([values(ta) for ta in TA]...)

    order = sortperm(ts)
    if ndims(TA[1]) == 1 # Check for 1D to ensure values remains a 1D vector.
        return TimeArray(ts[order], val[order], prev_colnames, prev_meta)
    else
        return TimeArray(ts[order], val[order, :], prev_colnames, prev_meta)
    end
end

# map ######################

@generated function map(f, ta::TimeArray{T,N}) where {T,N}
    input_val  = (N == 1) ? :(values(ta)[i]) : :(vec(values(ta)[i, :]))
    output_val = (N == 1) ? :(vals[i]) : :(vals[i, :])

    output_vals = (N == 1) ? :(vals[order]) : :(vals[order, :])

    quote
        ts   = similar(timestamp(ta))
        vals = similar(values(ta))

        for i in eachindex(ta)
            @inbounds ts[i], $output_val = f(timestamp(ta)[i], $input_val)
        end

        order = sortperm(ts)
        TimeArray(ts[order], $output_vals, colnames(ta), meta(ta))
    end
end


text:
using Luxor, Colors

chord = [0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0];
adjacency = hcat([circshift(chord, p) for p in 0:11] ...);
laplacian = sum(chord) * eye(12) - adjacency;
eigvs = eigvecs(laplacian);

radius = 20
Drawing(16*2*radius, 10*2*radius)

eigv = eigvs[:, 4]
bluered = diverging_palette(240, 10, 200)

origin()
background(bluered[101])
setopacity(1.0)

h = sqrt(3)

xs = [0,  3,  2,  2,  0,  4,  2,  1,  1,  3,  3,  1]
ys = [0,  1,  0,  2,  1,  2,  1,  0,  2,  0,  2,  1]

for p in 1:length(xs)
    c = clamp(imag(eigv[p] * 2), -1, 1)
    c = Int(round(c * 99)) + 101
    setcolor(bluered[c])
    for oct_x in -1:1
        for oct_y in -1:1
            x = xs[p]
            y = ys[p]

            ngon(
                (x + mod(y, 2)/2 + 4oct_x + 1.5*oct_y - 2) * radius * h,
                (y + 3oct_y - 1) * radius * 1.5,
                radius, 6, pi/2, :fill)
            text("C")
        end
    end
end

finish()
preview()


text:
###### Periodic domains

abstract type PeriodicDomain{T} <: Domain{T} end

isperiodic(::PeriodicDomain) = true

canonicaldomain(::PeriodicDomain) = PeriodicSegment()

points(d::PeriodicDomain{T},n::Integer) where {T} =
    fromcanonical.(Ref(d), fourierpoints(real(eltype(T)),n))

fourierpoints(n::Integer) = fourierpoints(Float64,n)
fourierpoints(::Type{T},n::Integer) where {T<:Number} = convert(T,π)*collect(0:2:2n-2)/n

function indomain(x, d::PeriodicDomain{T}) where T
    y=tocanonical(d,x)
    if !isapprox(fromcanonical(d,y),x)
        return false
    end

    l=arclength(d)
    if isinf(l)
        abs(imag(y))<20eps(T)
    else
        abs(imag(y))/l<20eps(T)
    end
end

first(d::PeriodicDomain) = fromcanonical(d,0)
last(d::PeriodicDomain) = fromcanonical(d,2π)

rand(d::PeriodicDomain,k...) = fromcanonical.(Ref(d),2π*rand(k...)-π)
checkpoints(d::PeriodicDomain) = fromcanonical.(Ref(d),[1.223972,3.14,5.83273484])
boundary(d::PeriodicDomain) = EmptyDomain()

for op in (:rdirichlet,:ldirichlet,:lneumann,:rneumann,:ivp,:bvp)
    @eval $op(::PeriodicDomain) = error("Periodic domains do not have boundaries")
end

struct AnyPeriodicDomain <: PeriodicDomain{UnsetNumber} end
isambiguous(::AnyPeriodicDomain)=true

convert(::Type{D},::AnyDomain) where {D<:PeriodicDomain} = AnyPeriodicDomain()

include("PeriodicSegment.jl")
include("Circle.jl")
include("PeriodicLine.jl")
include("PeriodicCurve.jl")
include("Disk.jl")


text:
module Pendulum3D_1

println("... Pendulum3D_1.jl: Simple pendulum modeled with multibody components")

using Modia
@usingModiaPlot

Pendulum = Model(
    world     = Object3D(feature=Scene()),
    body      = Object3D(feature=Solid(massProperties=MassProperties(mass=1.0))),
    bodyFrame = Object3D(parent=:body, translation=[-0.5, 0.0, 0.0]),
    rev       = Revolute(obj1=:world, obj2=:bodyFrame)
)

pendulum = @instantiateModel(buildModia3D(Pendulum))
simulate!(pendulum, stopTime=3.0)
plot(pendulum, "rev.phi")
end


text:
export
    KnitroProblem,
    createProblem, freeProblem,
    initializeProblem,
    solveProblem,
    restartProblem,
    setFuncCallback, setGradCallback, setHessCallback,
    setCallbacks, setMIPCallback,
    loadOptionsFile,
    loadTunerFile,
    setOption, getOption,
    applicationReturnStatus

mutable struct KnitroProblem
    # For KNITRO
    env::Ptr{Nothing} # pointer to KTR_context
    eval_status::Int32 # scalar input used only for reverse comms
    status::Int32  # Final status
    mip::Bool # whether it is a Mixed Integer Problem

    # For MathProgBase
    x::Vector{Float64}  # Starting and final solution
    lambda::Vector{Float64}
    g::Vector{Float64}  # Final constraint values
    obj_val::Vector{Float64}  # (length 1) Final objective

    # Callbacks
    eval_f::Function
    eval_g::Function
    eval_grad_f::Function
    eval_jac_g::Function
    eval_h::Function
    eval_hv::Function
    eval_mip_node::Function

    function KnitroProblem()
        kp = new(newcontext(),
                    0,
                    100, # Code for :Uninitialized
                    false)
        finalizer(freeProblem, kp)
        kp
    end
end

createProblem() = KnitroProblem()

function freeProblem(kp::KnitroProblem)
    kp.env == C_NULL && return
    return_code = @ktr_ccall(free, Int32, (Ptr{Nothing},), [kp.env])
    if return_code != 0
        error("KNITRO: Error freeing memory")
    end
    kp.env = C_NULL
end

function initializeKP(kp, x0, lambda0, g; mip = false)
    kp.status = 101 # code for :Initialized
    kp.mip = mip
    kp.x = x0
    kp.lambda = lambda0
    kp.g = g
    kp.obj_val = zeros(Float64, 1)
end

function initializeProblem(kp, objGoal, objType, x_l, x_u, c_Type, g_lb,
                            g_ub, jac_var, jac_con, hess_row, hess_col;
                            initial_x = C_NULL, initial_lambda = C_NULL)
    initializeKP(kp, (initial_x != C_NULL) ? initial_x :
                        zeros(Float64, length(x_l)),
                        (initial_lambda != C_NULL) ? initial_lambda :
                        zeros(Float64, length(x_l) + length(g_lb)),
                        zeros(Float64, length(g_lb)))
    init_problem(kp, objGoal, objType, x_l, x_u, c_Type, g_lb, g_ub,
                    jac_var, jac_con, hess_row, hess_col; initial_x=kp.x,
                    initial_lambda=kp.lambda)
end

function initializeProblem(kp, objGoal, objType, x_l, x_u, c_Type, g_lb,
                            g_ub, jac_var, jac_con; initial_x = C_NULL,
                            initial_lambda = C_NULL)
    hessopt = zeros(Int32 , 1)
    getOption(kp, "hessopt", hessopt)
    @assert hessopt[1] != KTR_HESSOPT_EXACT
    # KNITRO documentation:
    # If user option hessopt is not set to KTR_HESSOPT_EXACT, then Hessian
    # nonzeros will not be used. In this case, set nnzH=0, and pass NULL
    # pointers for hessIndexRows and hessIndexCols.
    initializeKP(kp, (initial_x != C_NULL) ? initial_x :
                        zeros(Float64, length(x_l)),
                        (initial_lambda != C_NULL) ? initial_lambda :
                        zeros(Float64, length(x_l) + length(g_lb)),
                        zeros(Float64, length(g_lb)))
    init_problem(kp, objGoal, objType, x_l, x_u, c_Type, g_lb, g_ub,
                    jac_var, jac_con; initial_x=kp.x,initial_lambda=kp.lambda)
end

# Initialization for MIP
function initializeProblem(kp, objGoal, objType, objFnType,
                            x_Type, x_l, x_u, c_Type, c_FnType, g_lb,
                            g_ub, jac_var, jac_con, hess_row, hess_col;
                            initial_x = C_NULL, initial_lambda = C_NULL)
    initializeKP(kp, (initial_x != C_NULL) ? initial_x :
                        zeros(Float64, length(x_l)),
                        (initial_lambda != C_NULL) ? initial_lambda :
                        zeros(Float64, length(x_l) + length(g_lb)),
                        zeros(Float64, length(g_lb)),
                        mip=true)
    mip_init_problem(kp, objGoal, objType, objFnType, x_Type, x_l, x_u,
                        c_Type, c_FnType, g_lb, g_ub, jac_var, jac_con,
                        hess_row, hess_col; initial_x=kp.x,
                        initial_lambda=kp.lambda)
end

function initializeProblem(kp, objGoal, objType, objFnType,
                            x_Type, x_l, x_u, c_Type, c_FnType, g_lb,
                            g_ub, jac_var, jac_con; initial_x = C_NULL,
                            initial_lambda = C_NULL)
    hessopt = Array(Int32, 1)
    getOption(kp, "hessopt", hessopt)
    @assert hessopt[1] != KTR_HESSOPT_EXACT
    # KNITRO documentation:
    # If user option hessopt is not set to KTR_HESSOPT_EXACT, then Hessian
    # nonzeros will not be used. In this case, set nnzH=0, and pass NULL
    # pointers for hessIndexRows and hessIndexCols.
    initializeKP(kp, (initial_x != C_NULL) ? initial_x :
                        zeros(Float64, length(x_l)),
                        (initial_lambda != C_NULL) ? initial_lambda :
                        zeros(Float64, length(x_l) + length(g_lb)),
                        zeros(Float64, length(g_lb)),
                        mip=true)
    mip_init_problem(kp, objGoal, objType, objFnType, x_Type, x_l, x_u,
                        c_Type, c_FnType, g_lb, g_ub, jac_var, jac_con;
                        initial_x=kp.x, initial_lambda=kp.lambda)
end

function solveProblem(kp::KnitroProblem)
    # Ensure that Knitro will not segfault when using Julia callbacks.
    setOption(kp, "par_numthreads", 1)
    setOption(kp, "par_msnumthreads", 1)
    if kp.mip
        kp.status = mip_solve_problem(kp, kp.x, kp.lambda, kp.eval_status,
                                        kp.obj_val)
    else
        kp.status = solve_problem(kp, kp.x, kp.lambda, kp.eval_status,
                                    kp.obj_val)
    end
end

function solveProblem(kp::KnitroProblem,
                        cons::Vector{Float64},
                        objGrad::Vector{Float64},
                        jac::Vector{Float64},
                        hess::Vector{Float64},
                        hessVector::Vector{Float64})
    # Ensure that Knitro will not segfault when using Julia callbacks.
    setOption(kp, "par_numthreads", 1)
    setOption(kp, "par_msnumthreads", 1)
    if kp.mip
        kp.status = mip_solve_problem(kp, kp.x, kp.lambda, kp.eval_status,
                                        kp.obj_val, cons, objGrad, jac, hess,
                                        hessVector)
    else
        kp.status = solve_problem(kp, kp.x, kp.lambda, kp.eval_status,
                                    kp.obj_val, cons, objGrad, jac, hess,
                                    hessVector)
    end
end

function restartProblem(kp, x0, lambda0)
    kp.status = 101 # code for :Initialized
    kp.eval_status = 0
    restart_problem(kp, x0, lambda0)
end

function restartProblem(kp)
    kp.status = 101 # code for :Initialized
    kp.eval_status = 0
    restart_problem(kp)
end

# -----
# Callback Wrappers
# -----
function eval_fc_wrapper(evalRequestCode::Cint,
                            n::Cint,
                            m::Cint,
                            nnzJ::Cint,
                            nnzH::Cint,
                            x_::Ptr{Cdouble},
                            lambda_::Ptr{Cdouble},
                            obj_::Ptr{Cdouble},
                            c_::Ptr{Cdouble},
                            g_::Ptr{Cdouble},
                            J_::Ptr{Cdouble},
                            H_::Ptr{Cdouble},
                            HV_::Ptr{Cdouble},
                            userParams_::Ptr{Nothing})
    if evalRequestCode != KTR_RC_EVALFC
        return KTR_RC_CALLBACK_ERR
    end
    kp = unsafe_pointer_to_objref(userParams_)::KnitroProblem
    x = unsafe_wrap(Array,x_,n)

    # calculate the new objective function value
    unsafe_store!(obj_, kp.eval_f(x))
    # calculate the new constraint values
    kp.eval_g(x,unsafe_wrap(Array,c_,m))

    Int32(0)
end

function eval_ga_wrapper(evalRequestCode::Cint,
                            n::Cint,
                            m::Cint,
                            nnzJ::Cint,
                            nnzH::Cint,
                            x_::Ptr{Cdouble},
                            lambda_::Ptr{Cdouble},
                            obj_::Ptr{Cdouble},
                            c_::Ptr{Cdouble},
                            g_::Ptr{Cdouble},
                            J_::Ptr{Cdouble},
                            H_::Ptr{Cdouble},
                            HV_::Ptr{Cdouble},
                            userParams_::Ptr{Nothing})
    if evalRequestCode != KTR_RC_EVALGA
        return KTR_RC_CALLBACK_ERR
    end
    kp = unsafe_pointer_to_objref(userParams_)::KnitroProblem
    x = unsafe_wrap(Array,x_,n)

    # evaluate the gradient
    kp.eval_grad_f(x,unsafe_wrap(Array,g_,n))
    # evaluate the jacobian
    if m > 0
        kp.eval_jac_g(x,unsafe_wrap(Array,J_,nnzJ))
    end

    Int32(0)
end

function eval_hess_wrapper(evalRequestCode::Cint,
                            n::Cint,
                            m::Cint,
                            nnzJ::Cint,
                            nnzH::Cint,
                            x_::Ptr{Cdouble},
                            lambda_::Ptr{Cdouble},
                            obj_::Ptr{Cdouble},
                            c_::Ptr{Cdouble},
                            g_::Ptr{Cdouble},
                            J_::Ptr{Cdouble},
                            H_::Ptr{Cdouble},
                            HV_::Ptr{Cdouble},
                            userParams_::Ptr{Nothing})
    kp = unsafe_pointer_to_objref(userParams_)::KnitroProblem
    x = unsafe_wrap(Array,x_, n)
    lambda = unsafe_wrap(Array,lambda_, m+n)

    if evalRequestCode == KTR_RC_EVALH
        kp.eval_h(x, lambda, 1.0, unsafe_wrap(Array,H_, nnzH))
    elseif evalRequestCode == KTR_RC_EVALH_NO_F
        kp.eval_h(x, lambda, 0.0, unsafe_wrap(Array,H_, nnzH))
    elseif evalRequestCode == KTR_RC_EVALHV
        kp.eval_hv(x, lambda, 1.0, unsafe_wrap(Array,HV_, n))
    elseif evalRequestCode == KTR_RC_EVALHV_NO_F
        kp.eval_hv(x, lambda, 0.0, unsafe_wrap(Array,HV_, n))
    else
        return KTR_RC_CALLBACK_ERR
    end
    Int32(0)
end

function eval_mip_node_wrapper(evalRequestCode::Cint,
                            n::Cint,
                            m::Cint,
                            nnzJ::Cint,
                            nnzH::Cint,
                            x_::Ptr{Cdouble},
                            lambda_::Ptr{Cdouble},
                            obj_::Ptr{Cdouble},
                            c_::Ptr{Cdouble},
                            g_::Ptr{Cdouble},
                            J_::Ptr{Cdouble},
                            H_::Ptr{Cdouble},
                            HV_::Ptr{Cdouble},
                            userParams_::Ptr{Nothing})
    kp = unsafe_pointer_to_objref(userParams_)::KnitroProblem
    obj = unsafe_load(obj_)
    kp.eval_mip_node(kp,obj)
    Int32(0)
end

function setFuncCallback(kp::KnitroProblem,
                            eval_f::Function,
                            eval_g::Function)
    kp.eval_f = eval_f
    kp.eval_g = eval_g
    set_func_callback(kp,eval_fc_wrapper)
end

function setGradCallback(kp::KnitroProblem,
                            eval_grad_f::Function,
                            eval_jac_g::Function)
    kp.eval_grad_f = eval_grad_f
    kp.eval_jac_g = eval_jac_g
    set_grad_callback(kp,eval_ga_wrapper)
end

function setHessCallback(kp::KnitroProblem,
                            eval_h::Function,
                            eval_hv::Function)
    kp.eval_h = eval_h
    kp.eval_hv = eval_hv
    set_hess_callback(kp,eval_hess_wrapper)
end

function setCallbacks(kp::KnitroProblem,
                        eval_f::Function,
                        eval_g::Function,
                        eval_grad_f::Function,
                        eval_jac_g::Function,
                        eval_h::Function,
                        eval_hv::Function)
    setFuncCallback(kp, eval_f, eval_g)
    setGradCallback(kp, eval_grad_f, eval_jac_g)
    setHessCallback(kp, eval_h, eval_hv)
end

function setMIPCallback(kp::KnitroProblem, eval_mip_node::Function)
    kp.eval_mip_node = eval_mip_node
    set_mip_node_callback(kp,eval_mip_node_wrapper)
end

# Getters and Setters for Parameters/Options
loadOptionsFile(kp, filename) = load_param_file(kp, filename)
loadTunerFile(kp, filename) = load_tuner_file(kp, filename)
setOption(args...) = set_param(args...)
getOption(args...) = get_param(args...)

function applicationReturnStatus(kp::KnitroProblem)
    if kp.status == 100
        # chosen not to clash with any of the KTR_RC_* codes
        return :Uninitialized
    elseif kp.status == 101
        # chosen not to clash with any of the KTR_RC_* codes
        return :Initialized
    elseif kp.status == 0
        return :Optimal
    elseif 1 <= kp.status <= 11
        return :ReverseComms
    elseif -199 <= kp.status <= -100
        return :FeasibleApproximate
    elseif -299 <= kp.status <= -200
        return :Infeasible
    elseif kp.status == -300
        return :Unbounded
    elseif -499 <= kp.status <= -400
        return :UserLimit
    elseif -599 <= kp.status <= -500
        return :KnitroError
    else
        return :Undefined
    end
end


text:
using Documenter, GridCreation

makedocs(modules = [GridCreation], doctest = false, sitename = "GridCreation")
deploydocs(repo = "github.com/JuliaOptics/GridCreation.jl.git")

max_stars_repo_name: JuliaOptics/Gridding

max_issues_repo_name: JuliaOptics/Gridding

max_forks_repo_name: JuliaOptics/Gridding