text:
using RayTracer
using Parameters

@with_kw mutable struct Particle{Dim}
    x⃗::Point{Dim}
    v⃗::Vec{Dim}
    m::Float = 1.0
end
@with_kw struct ExternalForces{Dim}
    F⃗::Vec{Dim}
    α::Float = 0.0
end

function step!(particle::Particle,forces::ExternalForces;dt = 1.0)
    @unpack x⃗,v⃗,m = particle
    @unpack F⃗,α = forces
    particle.x⃗ = x⃗ + dt .* v⃗
    particle.v⃗ = v⃗ + dt .* (F⃗/m - α .* v⃗)
    return particle.x⃗
end

using Plots

plt = plot(
    1,
    xlim = (-3, 3),
    ylim = (0, 3),
    title = "Universe",
    marker = 2,
    xlabel = "x",
    ylabel = "y",
)
projectile = Particle(x⃗ =  Point(-3, 0),v⃗ =  Vec(2, 3))
forces = ExternalForces(F⃗ = Vec(0, -1.0), α = .3)

# build an animated gif by pushing new points to the plot, saving every 2th frame
@gif for i=1:50
    push!(plt, (projectile.x⃗.coords)...)
    step!(projectile,forces,dt = 0.1)
end every 2


text:
# Raw memory management

function cualloc(T::Type, len::Integer)
    dptr_box = ptrbox(DevicePtr{Void})
    nbytes = Int(len) * sizeof(T)
    @cucall(:cuMemAlloc, (Ptr{DevicePtr{Void}}, Csize_t), dptr_box, nbytes)
    return ptrunbox(dptr_box)
end

cumemset(p::DevicePtr{Void}, value::Cuint, len::Integer) =
    @cucall(:cuMemsetD32, (DevicePtr{Void}, Cuint, Csize_t), p, value, len)

free(p::DevicePtr{Void}) = @cucall(:cuMemFree, (DevicePtr{Void},), p)

ext: jl


text:
# You notice a second question on the back of the homework assignment:
#
# What is the largest magnitude you can get from adding only two of the
# snailfish numbers?
#
# Note that snailfish addition is not commutative - that is, x + y and y + x
# can produce different results.
#
# Again considering the last example homework assignment above:
#
# [[[0,[5,8]],[[1,7],[9,6]]],[[4,[1,2]],[[1,4],2]]]
# [[[5,[2,8]],4],[5,[[9,9],0]]]
# [6,[[[6,2],[5,6]],[[7,6],[4,7]]]]
# [[[6,[0,7]],[0,9]],[4,[9,[9,0]]]]
# [[[7,[6,4]],[3,[1,3]]],[[[5,5],1],9]]
# [[6,[[7,3],[3,2]]],[[[3,8],[5,7]],4]]
# [[[[5,4],[7,7]],8],[[8,3],8]]
# [[9,3],[[9,9],[6,[4,9]]]]
# [[2,[[7,7],7]],[[5,8],[[9,3],[0,2]]]]
# [[[[5,2],5],[8,[3,7]]],[[5,[7,5]],[4,4]]]
#
# The largest magnitude of the sum of any two snailfish numbers in this list is
# 3993. This is the magnitude of [[2,[[7,7],7]],[[5,8],[[9,3],[0,2]]]] +
# [[[0,[5,8]],[[1,7],[9,6]]],[[4,[1,2]],[[1,4],2]]], which reduces to
# [[[[7,8],[6,6]],[[6,0],[7,7]]],[[[7,8],[8,8]],[[7,9],[0,6]]]].
#
# What is the largest magnitude of any sum of two different snailfish numbers
# from the homework assignment?

# read input
lines = map(eachline("input.txt")) do l
  r = []
  stack = [r]

  # first character will always be an opening bracket - skip it
  for c in l[2:end]
    # skip commas - numbers are always 1 digit in input, so we don't need to do
    # anything special with commas
    c == ',' && continue
    if c == '['
      nextpair = []
      push!(stack[end], nextpair)
      push!(stack, nextpair)
    elseif c == ']'
      pop!(stack)
    else
      # must be a number
      push!(stack[end], c - '0')
    end
  end
  return r
end

# utility function to handle explosions
function do_explosion(l)
  last_number = nothing
  add_to_next = nothing
  has_explosion = false
  exploder = function (p, depth)
    has_explosion && isnothing(add_to_next) && return p
    if depth == 5 && !has_explosion
      if !isnothing(last_number)
        last_number[1][last_number[2]] += p[1]
      end
      add_to_next = p[2]
      has_explosion = true
      return 0
    end

    if isa(p[1], Int)
      if isnothing(add_to_next)
        last_number = (p, 1)
      else
        p[1] += add_to_next
        add_to_next = nothing
        return p
      end
    else
      p[1] = exploder(p[1], depth + 1)
    end
    if isa(p[2], Int)
      if isnothing(add_to_next)
        last_number = (p, 2)
      else
        p[2] += add_to_next
        add_to_next = nothing
        return p
      end
    else
      p[2] = exploder(p[2], depth + 1)
    end
    return p
  end
  exploder(l, 1)
  return has_explosion
end

# utility function for splits
function do_split(l)
  had_split = false
  splitter = function (p)
    had_split && return p
    if isa(p, Int)
      if p >= 10
        (d,r) = divrem(p, 2)
        had_split = true
        return Any[d, d+r]
      end
      return p
    end
    p[1] = splitter(p[1])
    p[2] = splitter(p[2])
    return p
  end
  splitter(l)
  return had_split
end

# sum
function addition(a, b)
  r = Any[deepcopy(a), deepcopy(b)]
  while true
    do_explosion(r) || do_split(r) || break
  end
  return r
end

# find magnitude
function magnitude(p)
  isa(p, Int) && return p
  3 * magnitude(p[1]) + 2 * magnitude(p[2])
end

# do eet
largest_magnitude = 0
for i in 1:length(lines), j in 1:length(lines)
  i == j && continue
  m = magnitude(addition(lines[i], lines[j]))
  if m > largest_magnitude
    global largest_magnitude = m
  end
end

println("Result: ", largest_magnitude)


text:
using GeneralizedGenerated
using NamedTupleTools: merge_recursive

leaf_setter(::Type{Tx}, ::Type{Ty}) where {Tx, Ty} = leaf_setter(schema(Tx), schema(Ty))

@gg function leaf_setter(x::NamedTuple, y::NamedTuple)
    x = schema(x)
    y = schema(y)
    _leaf_setter(x, y)
end

@gg function leaf_setter(x::Tuple, y::Tuple)
    x = schema(x)
    y = schema(y)
    _leaf_setter(x,y)
end

function _leaf_setter(x, y)
    (xnames, xbody) = exprify(x)
    (ynames, ybody) = exprify(y)

    xargs = Expr(:tuple, xnames...)
    yargs = Expr(:tuple, ynames...)

    (_, body) = exprify(merge_recursive(
        leaf_setter(x)(xnames...),
        leaf_setter(y)(ynames...)); rename=false)

    q =  quote
        (x,y) -> begin
            $xargs = x
            $yargs = y
            return $body
        end
    end

    return q
end


text:
using Test
using Base.Iterators
using EndyIndexes

const ▶ = StartBasedIdx()
const 🛑 = EndBasedIdx()

# TODO tests for offset arrays
@testset "Endy Indexes" begin
    bv, dv1, ev1, ev2  = 2, 7, 7, 4
    N = 12
    a, b   = 🛑, 🛑-bv
    ai, bi =    N,    N-bv

    c, d, e    = 🛑:🛑, 🛑-dv1:🛑, 🛑-ev1:🛑-ev2
    ci, di, ei =    N:N,       N-dv1:N,       N-ev1:   N-ev2

    f, g, h    = 1:🛑, 2:🛑-bv,     🛑-dv1:11
    fi, gi, hi = 1:N,  2:N-bv,      N-dv1:11

    n, o   = ▶, ▶+bv
    ni, oi = 1,  1+bv

    p, q, r    = ▶:▶,   ▶:▶+dv1, ▶+ev2:▶+ev1
    pi, qi, ri =  1:1,    1:1+dv1,  1+ev2:1+ev1

    s, t, u    = ▶:N, 2:▶+dv1,     ▶+dv1:11
    si, ti, ui = 1:N,  2:1+dv1,      1+dv1:11

    endy_inds  = (a, b, c, d, e, f, g, h)
    iendy_inds = (ai, bi, ci, di, ei, fi, gi, hi)

    starty_inds  = (n, o, p, q, r, s, t, u)
    istarty_inds = (ni, oi, pi, qi, ri, si, ti, ui)

    @testset "display end" begin
        @test string(a) == "end"
        @test string(b) == "end-$bv"
        @test string(c) == "end:end"
        @test string(d) == "end-$dv1:end"
        @test string(e) == "end-$ev1:end-$ev2"
    end

    @testset "display start" begin
        @test string(n) == "start"
        @test string(o) == "start+$bv"
        @test string(p) == "start:start"
        @test string(q) == "start:start+$dv1"
        @test string(r) == "start+$ev2:start+$ev1"
    end

    @testset "endy types" begin
        @test a isa EndBasedIdx
        @test b isa EndBasedIdx
        @test c isa EndyRangeIdx
        @test d isa EndyRangeIdx
        @test e isa EndyRangeIdx
    end

    @testset "starty types" begin
        @test n isa StartBasedIdx
        @test o isa StartBasedIdx
        @test p isa EndyRangeIdx
        @test q isa EndyRangeIdx
        @test r isa EndyRangeIdx
    end

    @testset "single endy indexing" begin
        arr = [1:N;]
        @test arr[a] == arr[ai]
        @test arr[b] == arr[bi]
        @test arr[c] == arr[ci]
        @test arr[d] == arr[di]
        @test arr[e] == arr[ei]
    end

    @testset "single starty indexing" begin
        arr = [1:N;]
        @test arr[n] == arr[ni]
        @test arr[o] == arr[oi]
        @test arr[p] == arr[pi]
        @test arr[q] == arr[qi]
        @test arr[r] == arr[ri]
    end

    @testset "double endy indexing" begin
        arr = reshape([1:N^2;], N, N)

        all_endy_combos = product(endy_inds, endy_inds)
        all_inty_combos = product(iendy_inds, iendy_inds)
        for (inds, inds_ints) in zip(all_endy_combos, all_inty_combos)
            @show inds inds_ints
            @test arr[inds...] == arr[inds_ints...]
        end
    end

    @testset "double starty indexing" begin
        arr = reshape([1:N^2;], N, N)

        all_starty_combos = product(starty_inds, starty_inds)
        all_inty_combos = product(istarty_inds, istarty_inds)
        for (inds, inds_ints) in zip(all_starty_combos, all_inty_combos)
            @show inds inds_ints
            @test arr[inds...] == arr[inds_ints...]
        end
    end
end


text:
### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("FVheat1d");

useLog("FVheat1dlog")

# Configuration setup
domain(1)
solverType(FV)
timeStepper(RK4)

# Mesh
n = 40 # number of elements
mesh(LINEMESH, elsperdim=n)

# Variables and BCs
u = variable("u", location=CELL)
boundary(u, 1, FLUX, "0")

# Time interval and initial condition
T = 0.1;
timeInterval(T)
initial(u, "x<0.5 ? sin(pi*x)^4 : 0")

# The flux and source terms of the conservation equation
# F and S in the following equation:
# Dt(int(u dx)) = int(S dx) - int(F ds)
coefficient("D", 0.1) # Diffusion rate
flux(u, "-D*dot(grad(u),normal())") # in 1D -> -D * dx(u) * normal

exportCode("fvheat1dcode") # uncomment to export generated code to a file
#importCode("fvheat1dcode") # uncomment to import code from a file

solve(u)

finalize_finch()

##### Uncomment below to plot
# The initial condition
# u0 = zeros(n);
# x = Finch.fv_info.cellCenters[:]
# for i=1:n
#     u0[i] = x[i]<0.5 ? sin(pi*x[i])^4 : 0
# end

# using Plots
# pyplot();
# display(plot([x x], [u0 u.values[:]], markershape=:circle, label=["initial" "t="*string(T)]))

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
using DSP, FFTW, Test
using DSP: allocate_output
using Random: seed!

testfiles = ["dsp.jl", "util.jl", "windows.jl", "filter_conversion.jl",
    "filter_design.jl", "filter_response.jl", "filt.jl", "filt_stream.jl",
    "periodograms.jl", "multitaper.jl", "resample.jl", "lpc.jl", "estimation.jl", "unwrap.jl",
    "remez_fir.jl" ]

seed!(1776)

for testfile in testfiles
    eval(:(@testset $testfile begin include($testfile) end))
end


text:
module QOI

using FixedPointNumbers
using ColorTypes
using FileIO

#############
# Constants #
#############

const QOI_OP_INDEX = 0x00 # 00xxxxxx
const QOI_OP_DIFF  = 0x40 # 01xxxxxx
const QOI_OP_LUMA  = 0x80 # 10xxxxxx
const QOI_OP_RUN   = 0xc0 # 11xxxxxx
const QOI_OP_RGB   = 0xfe # 11111110
const QOI_OP_RGBA  = 0xff # 11111111
const QOI_MASK_2   = 0xc0 # 11000000
const QOI_MAGIC = UInt32('q') << 24 | UInt32('o') << 16 | UInt32('i') << 8 | UInt32('f')
const QOI_PIXELS_MAX = 400000000
const QOI_PADDING = (0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x01)

#########
# Pixel #
#########

struct Pixel
    r::UInt8
    g::UInt8
    b::UInt8
    a::UInt8
    function Pixel(r::UInt8, g::UInt8, b::UInt8, a::UInt8)
        new(r, g, b, a)
    end
end

@inline _qoi_color_hash(p::Pixel) = p.r*3 + p.g*5 + p.b*7 + p.a*11

##############
# Exceptions #
##############

struct QOIException <: Exception
    msg::String
end
Base.showerror(io::IO, qoi::QOIException) = print(io, qoi.msg)

@noinline throw_magic_bytes_error(magic::UInt32)             = throw(QOIException("invalid magic bytes, got $(repr(magic)), expected $(repr(QOI_MAGIC))"))
@noinline throw_invalid_header_width(width::UInt32)          = throw(QOIException("invalid width in header, got $width"))
@noinline throw_invalid_header_height(height::UInt32)        = throw(QOIException("invalid height in header, got $height"))
@noinline throw_invalid_header_channels(channels::UInt8)     = throw(QOIException("invalid channels in header, got $channels"))
@noinline throw_invalid_header_colorspace(colorspace::UInt8) = throw(QOIException("invalid colorspace in header, got $colorspace"))
@noinline throw_unexpected_eof()                             = throw(QOIException("unexpected end of file"))

#############
# QOIHeader #
#############

@enum QOIChannel::UInt8 begin
    QOI_RGB  = 0x03
    QOI_RGBA = 0x04
end

@enum QOIColorSpace::UInt8 begin
    QOI_SRGB   = 0x00
    QOI_LINEAR = 0x01
end

struct QOIHeader
    width::UInt32
    height::UInt32
    channels::QOIChannel
    colorspace::QOIColorSpace
end
function QOIHeader(width::UInt32, height::UInt32, channels::UInt8, colorspace::UInt8)
    width == 0                   && throw_invalid_header_width(width)
    height == 0                  && throw_invalid_header_height(height)
    channels < 3 || channels > 4 && throw_invalid_header_channels(channels)
    colorspace > 1               && throw_invalid_header_colorspace(colorspace)
    return QOIHeader(width, height, QOIChannel(channels), QOIColorSpace(colorspace))
end

############
# Encoding #
############

mutable struct QOIWriter{V <: AbstractVecOrMat{UInt8}}
    v::V
    pos::Int
end
QOIWriter(v::AbstractVecOrMat{UInt8}) = QOIWriter(v, 0)

@inline function _qoi_write!(qoiw::QOIWriter, v::UInt8)
    qoiw.pos += 1
    qoiw.pos > length(qoiw.v) && resize!(qoiw.v, max(1, length(qoiw.v) * 2))
    @inbounds qoiw.v[qoiw.pos] = v
end

function _qoi_write_32!(qoiw::QOIWriter, v::UInt32)
    _qoi_write!(qoiw, ((0xff000000 & v) >> 24) % UInt8)
    _qoi_write!(qoiw, ((0x00ff0000 & v) >> 16) % UInt8)
    _qoi_write!(qoiw, ((0x0000ff00 & v) >> 8)  % UInt8)
    _qoi_write!(qoiw, ((0x000000ff & v))       % UInt8)
    return
end

qoi_encode_raw(image::AbstractVecOrMat{UInt8}, header::QOIHeader) =
    qoi_encode_raw!(Vector{UInt8}(undef, 256), image, header)

function qoi_encode_raw(io::IO, image::AbstractVecOrMat{UInt8}, header::QOIHeader)
    data = qoi_encode_raw(image, header)
    write(io, data)
end

function qoi_encode_raw!(data::AbstractVector{UInt8}, image::AbstractVecOrMat{UInt8}, header::QOIHeader)
    # TODO: Check size of data against QOI_PIXELS_MAX?

    qoiw = QOIWriter(data)

    # Header
    _qoi_write_32!(qoiw, QOI_MAGIC)
    _qoi_write_32!(qoiw, header.width)
    _qoi_write_32!(qoiw, header.height)
    _qoi_write!(qoiw, header.channels |> Integer)
    _qoi_write!(qoiw, header.colorspace |> Integer)

    index = fill(Pixel(0x00, 0x00, 0x00, 0x00), 64)
    run = 0x00
    px_prev = Pixel(0x00, 0x00, 0x00, 0xff)
    px = px_prev

    channels = Integer(header.channels)
    px_len = header.width * header.height * channels
    px_end = px_len - channels + 1

    # Data
    for px_pos in 1:channels:px_len
        r = image[px_pos + 0]
        g = image[px_pos + 1]
        b = image[px_pos + 2]
        a = header.channels == QOI_RGBA ? image[px_pos + 3] : 0xff
        px = Pixel(r, g, b, a)

        if px == px_prev
            run += 0x01
            if run == 62 || px_pos == px_end
                _qoi_write!(qoiw, QOI_OP_RUN |  (run-0x01))
                run = 0x00
            end
        else
            if run > 0
                _qoi_write!(qoiw, QOI_OP_RUN | (run-0x01))
                run = 0x00
            end

            index_pos = mod1(_qoi_color_hash(px)+1, 64) % UInt8
            if (@inbounds index[index_pos]) == px
                _qoi_write!(qoiw, QOI_OP_INDEX | (index_pos - 0x01))
            else
                @inbounds index[index_pos] = px
                if px.a == px_prev.a
                    vr = ((px.r) - (px_prev.r)) % Int8
                    vg = ((px.g) - (px_prev.g)) % Int8
                    vb = ((px.b) - (px_prev.b)) % Int8

                    vg_r = vr - vg
                    vg_b = vb - vg
                    if      vr > -3 && vr < 2 &&
                            vg > -3 && vg < 2 &&
                            vb > -3 && vb < 2
                        _qoi_write!(qoiw, QOI_OP_DIFF | ((vr + 0x02) % UInt8) << 4 | ((vg + 0x02) % UInt8) << 2 | (vb + 0x02)  % UInt8)
                    elseif  vg_r > -9 && vg_r < 8  &&
                            vg > -33  && vg   < 32 &&
                            vg_b > -9 && vg_b < 8
                        _qoi_write!(qoiw, QOI_OP_LUMA   | (vg + UInt8(32)) % UInt8)
                        _qoi_write!(qoiw, ((vg_r + 0x08) % UInt8) << 4 | (vg_b + 0x08) % UInt8)
                    else
                        _qoi_write!(qoiw, QOI_OP_RGB)
                        _qoi_write!(qoiw, px.r)
                        _qoi_write!(qoiw, px.g)
                        _qoi_write!(qoiw, px.b)
                    end
                else
                    _qoi_write!(qoiw, QOI_OP_RGBA)
                    _qoi_write!(qoiw, px.r)
                    _qoi_write!(qoiw, px.g)
                    _qoi_write!(qoiw, px.b)
                    _qoi_write!(qoiw, px.a)
                end
            end
        end
        px_prev = px
    end

    # Padding
    for x in QOI_PADDING
        _qoi_write!(qoiw, x)
    end

    sizehint!(data, qoiw.pos)
    resize!(data, qoiw.pos)
    return data
end

function qoi_encode(file::String, image::AbstractMatrix{T}) where T <: Colorant
    if T <: TransparentColor
        if T != RGBA{N0f8}
            image = RGBA{N0f8}.(image)
        end
        channel = QOI_RGBA
    else
        if T != RGB{N0f8}
            image = RGB{N0f8}.(image)
        end
        channel = QOI_RGB
    end
    header = QOIHeader(size(image, 2), size(image, 1), channel, QOI_SRGB)
    image = permutedims(image)
    open(file, "w") do io
        image_raw = reinterpret(UInt8, image)
        qoi_encode_raw(io, image_raw, header)
    end
end

############
# Decoding #
############

mutable struct QOIReader{V <: AbstractVecOrMat{UInt8}}
    v::V
    pos::Int
end
QOIReader(v::AbstractVecOrMat{UInt8}) = QOIReader(v, 0)

@inline _qoi_read!(qoir::QOIReader) = (qoir.pos+=1; @inbounds qoir.v[qoir.pos])

function _qoi_read_32!(qoir::QOIReader)
    a = UInt32(_qoi_read!(qoir))
    b = UInt32(_qoi_read!(qoir))
    c = UInt32(_qoi_read!(qoir))
    d = UInt32(_qoi_read!(qoir))
    return a << 24 | b << 16 | c << 8 | d
end

@inline _qoi_read_rgb!(qoir::QOIReader) = (_qoi_read!(qoir), _qoi_read!(qoir), _qoi_read!(qoir))
@inline _qoi_read_rgba!(qoir::QOIReader) = (_qoi_read!(qoir), _qoi_read!(qoir), _qoi_read!(qoir), _qoi_read!(qoir))

function qoi_decode_raw(v::AbstractVector{UInt8})
    qoir = QOIReader(v)

    if length(v) < 14 # 4 magic + 10 header
        throw_unexpected_eof()
    end

    # Magic
    magic = _qoi_read_32!(qoir)
    magic == QOI_MAGIC || throw_magic_bytes_error(magic)

    # Header
    width = _qoi_read_32!(qoir)
    height = _qoi_read_32!(qoir)
    channels = _qoi_read!(qoir)
    colorspace = _qoi_read!(qoir)
    header = QOIHeader(width, height, channels, colorspace)

    # Data
    n_pixels = header.width * header.height
    n_values = n_pixels * channels
    data = Vector{UInt8}(undef, n_values)
    index = fill(Pixel(0x00, 0x00, 0x00, 0x00), 64)
    px = Pixel(0x00, 0x00, 0x00, 0xFF)
    px_idx = 1
    run = 0x00
    lenv = length(v)

    for px_idx in 1:channels:n_values
        if qoir.pos > lenv - length(QOI_PADDING)
            throw_unexpected_eof()
        end
        if run > 0
            run -= 0x01
        else
            b1 = _qoi_read!(qoir)
            if b1 == QOI_OP_RGB
                r, g, b = _qoi_read_rgb!(qoir)
                px = Pixel(r, g, b, px.a)
            elseif b1 == QOI_OP_RGBA
                r, g, b, a = _qoi_read_rgba!(qoir)
                px = Pixel(r, g, b, a)
            elseif b1 & QOI_MASK_2 == QOI_OP_INDEX
                px = index[b1+0x01]
            elseif (b1 & QOI_MASK_2) == QOI_OP_DIFF
                r = px.r + ((b1 >> 0x04) & 0x03) - 0x02
                g = px.g + ((b1 >> 0x02) & 0x03) - 0x02
                b = px.b + ( b1          & 0x03) - 0x02
                px = Pixel(r, g, b, px.a)
            elseif ((b1 & QOI_MASK_2) == QOI_OP_LUMA)
                b2 = _qoi_read!(qoir)
                vg = (b1 & 0x3f) - UInt8(32)
                r = px.r + vg - 0x08 + ((b2 >> 4) & 0x0f)
                g = px.g + vg
                b = px.b + vg - 0x08 +  (b2       & 0x0f)
                px = Pixel(r, g, b, px.a)
            elseif (b1 & QOI_MASK_2) == QOI_OP_RUN
                run = (b1 & 0x3f)
            else
                error("unreachable")
            end
            @inbounds index[mod1(_qoi_color_hash(px)+1, 64)] = px
        end

        @inbounds data[px_idx+0] = px.r
        @inbounds data[px_idx+1] = px.g
        @inbounds data[px_idx+2] = px.b
        if header.channels == QOI_RGBA
            @inbounds data[px_idx+3] = px.a
        end
    end

    if qoir.pos != length(v) - length(QOI_PADDING)
        throw_unexpected_eof()
    end

    # Read padding
    for _ in 1:7
        x = _qoi_read!(qoir)
        x == 0 || throw_unexpected_eof()
    end
    x = _qoi_read!(qoir)
    x == 1 || throw_unexpected_eof()

    image = reshape(data, Int(header.width) * channels, Int(header.height))
    return header, image
end

function _to_colortype(header, raw_image)
    T = header.channels == QOI_RGBA ? RGBA{N0f8} : RGB{N0f8}
    return permutedims(reinterpret(T, raw_image))
end

function qoi_decode(v::AbstractVector{UInt8})
    header, raw_image = qoi_decode_raw(v)
    return _to_colortype(header, raw_image)
end

qoi_decode_raw(f::Union{String, IO}) = qoi_decode_raw(Base.read(f))
qoi_decode(f::Union{String, IO}) = qoi_decode(Base.read(f))

##########
# FileIO #
##########

load(f::File{format"QOI"}) = qoi_decode(f.filename)
save(f::File{format"QOI"}, image::AbstractMatrix{<:Colorant}) = qoi_encode(f.filename, image)

end


text:
_book_depot_dirname() = ".researchbook"
depotdir() = joinpath(bookdir(), _book_depot_dirname())
bookdir_relpath(path) = abspath(joinpath(bookdir(), path))


text:
include("../../mesh.jl")
include("../../utils.jl")
include("../../connectique.jl")
# include("../../visu.jl")

include("toyPbs.jl")

function validateToyProblem(nod2e::Array{Array{Int64, 1}, 1},
    e2neighs::Array{Array{Int64, 1}, 1}, nedg::Int64, e2edg::Array{Array{Int64, 1}, 1},
    freeedg::Array{Bool, 1}, freenod::Array{Bool, 1}, nod2e_ref::Array{Array{Int64, 1}, 1},
    e2neighs_ref::Array{Array{Int64, 1}, 1}, nedg_ref::Int64,
    e2edg_ref::Array{Array{Int64, 1}, 1}, freeedg_ref::Array{Bool, 1}, freenod_ref::Array{Bool, 1})
    # checking nod2e
    println("---------------------------------------------------")
    @assert length(nod2e)  == length(nod2e_ref)
    for n in eachindex(nod2e)
        println("Checking node: $n")
        @assert length(nod2e[n]) == length(nod2e_ref[n])
        for ne in eachindex(nod2e[n])
            println("  -> $(nod2e[n][ne])    -> $(nod2e_ref[n][ne])")
        end
    end
    # checking e2neighs
    println("---------------------------------------------------")
    @assert length(e2neighs) == length(e2neighs_ref)
    for e in eachindex(e2neighs)
        println("Checking element: $e")
        @assert length(e2neighs[e]) == length(e2neighs_ref[e])
        for ee in eachindex(e2neighs[e])
            println("  -> $(e2neighs[e][ee])    -> $(e2neighs_ref[e][ee])")
        end
    end
    # checking nedg
    println("---------------------------------------------------")
    println("Checking nedg: $nedg     $nedg_ref")
    # checking
    println("---------------------------------------------------")
    @assert length(e2edg) == length(e2edg_ref)
    for e in eachindex(e2edg)
        println("Checking element: $e")
        @assert length(e2edg[e]) == length(e2edg_ref[e])
        for ee in eachindex(e2edg[e])
            println("  -> $(e2edg[e][ee])    -> $(e2edg_ref[e][ee])")
        end
    end
    # check free edges
    println("---------------------------------------------------")
    println("Checking free edges")
    for e in eachindex(freeedg)
        println("  -> $(freeedg[e])    -> $(freeedg_ref[e])")
    end
    println("---------------------------------------------------")
    println("Checking free nodes")
    for n in eachindex(freenod)
        println("  -> $(freenod[n])    -> $(freenod_ref[n])")
    end
    nothing
end

function mainValidation2()
    vtx, elt, bnd = generateToyProblem2()
    # get mesh data
    e2neigh, nod2e = getNeighbours(vtx, elt)
    nedg           = getNbEdges(e2neigh)
    e2edg, freeedg, freenod = getMeshEdges(vtx, elt, nedg, e2neigh)
    # get correct mesh data
    nod2e_ref, e2neigh_ref, nedg_ref, e2edg_ref, freeedg_ref, freenod_ref =
        correctDataToyProblem2()
    # now check
    validateToyProblem(nod2e, e2neigh, nedg, e2edg, freeedg, freenod,
        nod2e_ref, e2neigh_ref, nedg_ref, e2edg_ref, freeedg_ref, freenod_ref)
    # connectivity
    msh = mesh(vtx, elt, e2edg, freeedg, freenod, bnd, nedg)
    conn = connectivity_k0(msh)
    conn_ref = connectivityToyProblem2()
    println("---------------------------------------------------")
    println("Checking connectivity")
    @assert length(conn.e2dof) == length(conn_ref)
    for e in eachindex(conn.e2dof)
        @assert length(conn.e2dof[e]) == length(conn_ref[e])
        for nd in eachindex(conn.e2dof[e])
            println("  -> $(conn.e2dof[e][nd])    -> $(conn_ref[e][nd])")
        end
    end
    #
    nothing
end

mainValidation2()


text:
# Routines for observations of light sources in particular images,
# rather than for sources in the abstract, in physical units,
# and rather than for images alone (that's image_model.jl).

"""
Attributes of the patch of sky surrounding a single
celestial object in a single image.

Attributes:
  - center: The approximate source location in world coordinates
  - radius_pix: The width of the influence of the object in pixel coordinates
  - psf: The point spread function in this region of the sky
  - wcs_jacobian: The jacobian of the WCS transform in this region of the
                  sky for each band
  - pixel_center: The pixel location of center in each band.
"""
struct SkyPatch
    center::Vector{Float64}
    radius_pix::Float64

    psf::Vector{PsfComponent}
    wcs_jacobian::Matrix{Float64}
    pixel_center::Vector{Float64}

    bitmap_offset::SVector{2, Int64}  # lower left corner index offset
    active_pixel_bitmap::Matrix{Bool}
end

"""
Choose a reasonable patch radius based on the catalog.

Args:
  - ce: The catalog entry for the object
  - img: The image
  - width_scale: A multiple of standard deviations to use
  - max_radius: The maximum radius in pixels.

Returns:
  - A radius in pixels chosen from the catalog entry.
"""
function choose_patch_radius(ce::CatalogEntry,
                             img::Image;
                             width_scale=1.0,
                             max_radius=25)
    psf_width = Model.get_psf_width(img.psf, width_scale=width_scale)

    # The galaxy scale is the point with half the light -- if the light
    # were entirely in a univariate normal, this would be at 0.67 standard
    # deviations.  We are being a bit conservative here.
    obj_width = ce.is_star ? 0.0 : width_scale * ce.gal_scale / 0.67
    obj_width += psf_width

    flux = ce.is_star ? ce.star_fluxes[img.b] : ce.gal_fluxes[img.b]
    if (!(flux > 0.)) @show ce end
    @assert flux > 0.

    # Choose enough pixels that the light is either 80% of the light
    # would be captured from a 1d gaussian or 5% of the sky noise,
    # whichever is a larger radius.
    epsilon = img.sky[div(img.H, 2), div(img.W, 2)]
    pdf_90 = exp(-0.5 * (1.64)^2) / (sqrt(2pi) * obj_width)
    pdf_target = min(pdf_90, epsilon / (20 * flux))
    rhs = log(pdf_target) + 0.5 * log(2pi) + log(obj_width)
    radius_req = sqrt(-2 * (obj_width ^ 2) * rhs)

    min(radius_req, max_radius)
end


text:
import Base: ==

abstract type Message end


text:
module SuffixArrays

export suffixsort

const CodeUnits = Union{UInt8,UInt16}
const IndexTypes = Union{Int8,Int16,Int32,Int64}
const IndexVector = AbstractVector{<:IndexTypes}

include("sais.jl")

function suffixsort(V::AbstractVector{U}, base::Integer=1) where {U<:CodeUnits}
    0 ≤ base || throw(ArgumentError("unsupported negative indexing base: $base"))
    n = length(V)
    # unsigned index type to return
    T = n+base-1 ≤ typemax(UInt8)  ? UInt8  :
        n+base-1 ≤ typemax(UInt16) ? UInt16 :
        n+base-1 ≤ typemax(UInt32) ? UInt32 : UInt64
    n ≤ 1 && return fill(T(base), n)
    # signed index type for algorithm
    S = n ≤ typemax(Int8)  ? Int8  :
        n ≤ typemax(Int16) ? Int16 :
        n ≤ typemax(Int32) ? Int32 : Int64
    if sizeof(T) == sizeof(S)
        I = zeros(T, n)
        sais(V, reinterpret(S, I), 0, n, Int(typemax(U))+1, false)
        base ≠ 0 && (I .+= base)
        return I
    else
        I = zeros(S, n)
        sais(V, I, 0, n, Int(typemax(U))+1, false)
        I′ = Vector{T}(undef, n)
        @inbounds for (i, x) in enumerate(I)
            I′[i] = (x + base) % T
        end
        return I′
    end
end

function suffixsort(s::AbstractString, base::Integer=1)
    return suffixsort(codeunits(s), base)
end

end # module


text:
using Documenter, MendelGeneDropping

ENV["DOCUMENTER_DEBUG"] = "true"

makedocs(
    format = :html,
    sitename = "MendelGeneDropping",
    modules = [MendelGeneDropping],
    authors = "Jeanette Papp",
    clean = true,
    debug = true,
    pages = [
        "index.md"
    ]
)

deploydocs(
    repo   = "github.com/OpenMendel/MendelGeneDropping.jl.git",
    target = "build",
    julia  = "1.0",
    deps   = nothing,
    make   = nothing
)


text:
using Random

Random.seed!(123)

# smooth case

f1 = SqrNormL2()
f2 = Translate(SqrNormL2(2.5), randn(10))
f = Sum(f1, f2)

predicates_test(f)

@test ProximalOperators.is_quadratic(f) == true
@test ProximalOperators.is_strongly_convex(f) == true
@test ProximalOperators.is_set(f) == false

xtest = randn(10)

result = f1(xtest) + f2(xtest)
@test f(xtest) ≈ result

grad1, val1 = gradient_test(f1, xtest)
grad2, val2 = gradient_test(f2, xtest)

gradsum, valsum = gradient_test(f, xtest)
@test gradsum ≈ grad1 + grad2

# nonsmooth case

g1 = NormL2()
g2 = Translate(SqrNormL2(2.5), randn(10))
g = Sum(g1, g2)

predicates_test(g)

@test ProximalOperators.is_smooth(g) == false
@test ProximalOperators.is_strongly_convex(g) == true
@test ProximalOperators.is_set(g) == false

xtest = randn(10)

result = g1(xtest) + g2(xtest)
@test g(xtest) ≈ result

grad1, val1 = gradient_test(g1, xtest)
grad2, val2 = gradient_test(g2, xtest)

gradsum, valsum = gradient_test(g, xtest)
@test gradsum ≈ grad1 + grad2


text:
using StructIO
using Base.Test

@struct immutable TwoUInts
    x::UInt
    y::UInt
end

buf = IOBuffer()
write(buf, UInt(1))
write(buf, UInt(2))
seekstart(buf)
@test unpack(buf, TwoUInts) == TwoUInts(1,2)

abstract SomeAbstractType
@struct immutable SomeConcreteType <: SomeAbstractType
    A::UInt32
    B::UInt16
    C::UInt32
    D::UInt8
end align_packed

@struct immutable ParametricType{S,T}
    A::S
    B::T
    C::T
end


text:
function  Bilderrekonstruktion(Auswahl,Aufloesung,prozent,p)

if Auswahl == "Bild1"
    Daten1 = load("../Computertomographie_Bilder/CTKrabbe.jpg");
elseif Auswahl == "Bild2"
    Daten1 = load("../Computertomographie_Bilder/CTEidechse.jpg")
elseif Auswahl == "Bild3"
    Daten1 = load("../Computertomographie_Bilder/CTFisch.jpg")
elseif Auswahl == "Bild4"
    Daten1 = load("../Computertomographie_Bilder/CTFledermaus.jpg")
elseif Auswahl == "Bild5"
    Daten1 = load("../Computertomographie_Bilder/CTChamaeleon.jpg")
elseif Auswahl == "Bild6"
    Daten1 = load("../Computertomographie_Bilder/CTSchildkroete.jpg")
elseif Auswahl == "Bild7"
    Daten1 = load("../Computertomographie_Bilder/CTSeepferdchen.jpg")
elseif Auswahl == "Bild8"
    Daten1 = load("../Computertomographie_Bilder/CTFrosch.jpg")
elseif Auswahl == "NaN"
    print("Ihr müsst das NaN durch Bild1, ... oder Bild8 ersetzen.");
end

Daten = Gray.(Daten1)
out = convert(Array{Float64,2}, Daten); #convert GrayMatrix into Float64, bei Vektor müsste die 2 durch eine 1ersetzt werden
Summe = zeros(Aufloesung,Aufloesung);
Durchschnitt = zeros(Aufloesung,Aufloesung);
Laenge = size(out,2)/Aufloesung;
SummeZeile = zeros(Int(Laenge),1);

for i = 1:Aufloesung #geht die Zeilen durch
    for j = 1:Aufloesung #geht die Spalten durch
        a = Int((i-1)*Laenge+1);
        b = Int(i*Laenge);
        c = Int((j-1)*Laenge+1);
        d = Int(j*Laenge);
        Matrix = out[a:b,c:d];
        SummeZeile = sum(Matrix, dims=2); # Jede Zeile wird aufaddiert und in die entsprechende Zeile des SummeZeilevektors geschrieben
        Summe = sum(SummeZeile,dims=1)/(Laenge^2);# Einträge von SummeZeile werden aufaddiertund normiert
        Durchschnitt[i,j] = Summe[1,1];
    end
end

nrows,ncols = size(Durchschnitt);
#f_nahBild = zeros(size(Durchschnitt,1),size(Durchschnitt,2));
DurchschnittVec = reshape(Durchschnitt,(Int(nrows)*Int(ncols),1));

x0,y0,N,JJ,nx,ny,x1,y1,verschx,verschy,a,b,A,theta1,counter = Dateneinlesen(Aufloesung);
A2 = erstelleA(x0,x1,y0,y1,N,JJ,nx,ny,verschx,verschy);
Messwerte =  A2*DurchschnittVec;

MesswerteMitFehler = erzeugeFehlerbehafteteDaten(Messwerte,prozent);

#f_nah = genTikhonov(A2,MesswerteMitFehler,p,Aufloesung);
#for i = 1:size(f_nah,1)
#    if f_nah[i] < 0
#        f_nah[i] = 0;
#    elseif f_nah[i] > 1
#        f_nah[i] = 1;
#    end
#end

f_nahfehler(p) = genTikhonov(A2,MesswerteMitFehler,p,Aufloesung);

if !isnan(p)
    Fehler = norm(DurchschnittVec-f_nahfehler(p))/norm(DurchschnittVec)*100;
    println("Der Fehler beträgt: ", Fehler)
    f_nahfehler = f_nahfehler(p);
    for i = 1:size(f_nahfehler,1)
        if f_nahfehler[i] < 0
            f_nahfehler[i] = 0;
        elseif f_nahfehler[i] > 1
            f_nahfehler[i] = 1;
        end
    end
    ObjektTikhonov = reshape(f_nahfehler,(nrows,ncols));
    f_nahBild = convert(Array{Gray{Normed{UInt8,8}},2},ObjektTikhonov);
else
    println("Ihr müsst einen Wert für p einsetzen.")
end
return f_nahBild;

end


text:
function finalgrade(exam::Number, projects::Number)::Int
    if exam > 90 || projects > 10
        100
    elseif exam > 75 && projects ≥ 5
        90
    elseif exam > 50 && projects ≥ 2
        75
    else
        0
     end
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test

function runtests(name, isolate=true; seed=nothing)
    old_print_setting = Test.TESTSET_PRINT_ENABLE[]
    Test.TESTSET_PRINT_ENABLE[] = false
    try
        if isolate
            # Simple enough to type and random enough so that no one will hard
            # code it in the test
            mod_name = Symbol("Test", rand(1:100), "Main_", replace(name, '/', '_'))
            m = @eval(Main, module $mod_name end)
        else
            m = Main
        end
        @eval(m, using Test)
        ex = quote
            @timed @testset $"$name" begin
                srand($seed)
                include($"$name.jl")
            end
        end
        res_and_time_data = eval(m, ex)
        rss = Sys.maxrss()
        #res_and_time_data[1] is the testset
        passes,fails,error,broken,c_passes,c_fails,c_errors,c_broken = Test.get_test_counts(res_and_time_data[1])
        if res_and_time_data[1].anynonpass == false
            res_and_time_data = (
                                 (passes+c_passes,broken+c_broken),
                                 res_and_time_data[2],
                                 res_and_time_data[3],
                                 res_and_time_data[4],
                                 res_and_time_data[5])
        end
        vcat(collect(res_and_time_data), rss)
    finally
        Test.TESTSET_PRINT_ENABLE[] = old_print_setting
    end
end

# looking in . messes things up badly
filter!(x->x!=".", LOAD_PATH)

nothing # File is loaded via a remotecall to "include". Ensure it returns "nothing".

max_stars_repo_name: anowacki/julia

max_issues_repo_name: anowacki/julia

max_forks_repo_name: anowacki/julia


text:
"""
PENBBScan2D()
Performs a scan moving in x-y axis
...
Function to perform an automate scan in 2D (x,y axis)
This function requires as input the starting point (x_start,y_start), and the end point (x_ends,y_ends) as well
as the step sizes in both axis(step_x,step_y). All values are in mm. The range specified muss be in the interval x[0.,10.], y[0.,10.]
In addition you can specify the name of the holder or sample and time of data taking for each position
PENBBScan2D(x_start, y_start, step_x, step_y, x_ends, y_ends, HolderName::String, time_per_point)
Example: PENBBScan2D(0.0,0.0,20.0,20.0,40.0,40.0,"small",2)
...
# Arguments
- x_start: intial point in x
- y_start: intial point in y
- step_x: step size in x
- step_y: step size in y
- x_ends: final point in x
- y_ends: final point in y
- HolderName::String: name of the holder of piece you are scanning
- time_per_point: time of data taking in each point
...
"""
function PENBBScan2D(x_start, y_start, step_x, step_y, x_ends, y_ends, HolderName::String, time_per_point)
    if x_start < 0.0 || y_start < 0.0 || x_ends > 100.0 || y_ends > 100.
        @info("Error: value out of range: you muss use values in the range x[0.,10.], y[0.,10.]")
    else
	motor = mymotor()
	Calibrate(motor)
        for i in  collect(x_start:step_x:x_ends)
            XMoveMM(i,motor)
            for j in collect(y_start:step_y:y_ends)
                @info("position: ",i,j)
                YMoveMM(j,motor)
                pos_x = string(PosX(motor))
                pos_y = string(PosY(motor))
                name_file = string("2D_PEN_Scan_Holder_",HolderName,"_x_",pos_x,"_y_",pos_y)
                settings = (fadc = "gelab-fadc08",
                    output_basename = name_file,
                    data_dir = "./",
                    conv_data_dir = "./",
                    measurement_time = time_per_point,
                    number_of_measurements = 1,
                    channels = [1,2,3,4,5,6,7],
                    trigger_threshold = 45,
                    trigger_pmt = 5,
                    peakTime = 2,
                    gapTime = 2,
                    nPreTrig = 70,
                    nSamples = 128,
                    saveEnergy = false,
                    delete_dat = true,
                    h5_filesize_limit = 1000,
                    filter_faulty_events = false,
                    coincidence_interval = 4e-9
                );
                println(name_file)
                take_struck_data(settings, calibration_data=false)
                sleep(2)
            end
        end
        @info("PEN BB 2D scan completed, see you soon!")
    end

end

"""
PENBBScan1DY()
Performs a scan moving in y axis
...
Function to perform an automate scan in 1D (y axis) keeping x constant
This function requires as input the starting point (x_start,y_start), and the end point (x_start,y_ends) as well
as the step size(step_y). All values are in mm. The range specified muss be in the interval x[0.,10.], y[0.,10.]
In addition you can specify the name of the holder or sample and time of data taking for each position
...
# Arguments
- x_start: intial point in x
- y_start: intial point in y
- step_y: step size in y
- y_ends: final point in y
- HolderName::String: name of the holder of piece you are scanning
- time_per_point: time of data taking in each point
...
"""
function PENBBScan1DY(x_start, y_start, step_y, y_ends, HolderName::String, time_per_point)
    if x_start < 0.0 || y_start < 0.0 || x_start > 100.0 || y_ends > 100.
        @info("Error: value out of range: you muss use soem value in the range x[0.,10.], y[0.,10.]")
    else
	motor = mymotor()
	Calibrate(motor)
        XMoveMM(x_start,motor)
        for j in collect(y_start:step_y:y_ends)
            @info("position: ",x_start,j)
            YMoveMM(j,motor)
            pos_x = string(PosX(motor))
            pos_y = string(PosY(motor))
            name_file = string("2D_PEN_Scan_Holder_",HolderName,"_x_",pos_x,"_y_",pos_y)
            settings = (fadc = "gelab-fadc08",
                    output_basename = name_file,
                    data_dir = "./",
                    conv_data_dir = "./",
                    measurement_time = time_per_point,
                    number_of_measurements = 1,
                    channels = [1,2,3,4,5,6,7],
                    trigger_threshold = 45,
                    trigger_pmt = 5,
                    peakTime = 2,
                    gapTime = 2,
                    nPreTrig = 70,
                    nSamples = 128,
                    saveEnergy = false,
                    delete_dat = true,
                    h5_filesize_limit = 1000,
                    filter_faulty_events = false,
                    coincidence_interval = 4e-9
            );
            println(name_file)
            take_struck_data(settings, calibration_data=false)
            sleep(2)
        end
        @info("PEN 1D BB scan completed, see you soon!")
    end

end

"""
PENBBScan1DX()
Performs a scan moving in x axis
...
Function to perform an automate scan in 1D (x axis) keeping y constant
This function requires as input the starting point (y_start,x_start), and the end point (y_start,x_ends) as well
as the step size(step_x). All values are in mm. The range specified muss be in the interval x[0.,10.], y[0.,10.]
In addition you can specify the name of the holder or sample and time of data taking for each position
...
# Arguments
- x_start: intial point in x
- y_start: intial point in y
- step_x: step size in y
- x_ends: final point in y
- HolderName::String: name of the holder of piece you are scanning
- time_per_point: time of data taking in each point
...
"""

function PENBBScan1DX(y_start, x_start, step_x, x_ends, HolderName::String, time_per_point)
    if x_start < 0.0 || y_start < 0.0 || x_start > 100.0 || y_start > 100.
        @info("Error: value out of range: you muss use soem value in the range x[0.,10.], y[0.,10.]")
    else
	motor = mymotor()
	Calibrate(motor)
        YMoveMM(y_start,motor)
        for i in collect(x_start:step_x:x_ends)
            @info("position: ",i,y_start)
            XMoveMM(i,motor)
            pos_x = string(PosX(motor))
            pos_y = string(PosY(motor))
            name_file = string("2D_PEN_Scan_Holder_",HolderName,"_x_",pos_x,"_y_",pos_y)
            settings = (fadc = "gelab-fadc08",
                    output_basename = name_file,
                    data_dir = "./",
                    conv_data_dir = "./",
                    measurement_time = time_per_point,
                    number_of_measurements = 1,
                    channels = [1,2,3,4,5,6,7],
                    trigger_threshold = 45,
                    trigger_pmt = 5,
                    peakTime = 2,
                    gapTime = 2,
                    nPreTrig = 70,
                    nSamples = 128,
                    saveEnergy = false,
                    delete_dat = true,
                    h5_filesize_limit = 1000,
                    filter_faulty_events = false,
                    coincidence_interval = 4e-9
            );
            println(name_file)
            take_struck_data(settings, calibration_data=false)
            sleep(2)
        end
        @info("PEN BB 1D scan completed, see you soon!")
    end

end


text:
# see: http://stackoverflow.com/a/37732384/5075246

@userplot PortfolioComposition

# this shows the shifting composition of a basket of something over a variable
# - "returns" are the dependent variable
# - "weights" are a matrix where the ith column is the composition for returns[i]
# - since each polygon is its own series, you can assign labels easily
@recipe function f(pc::PortfolioComposition)
    weights, returns = pc.args
    n = length(returns)
    weights = cumsum(weights,2)
    seriestype := :shape

	# create a filled polygon for each item
    for c=1:size(weights,2)
        sx = vcat(weights[:,c], c==1 ? zeros(n) : reverse(weights[:,c-1]))
        sy = vcat(returns, reverse(returns))
        @series Plots.isvertical(d) ? (sx, sy) : (sy, sx)
    end
end


text:
function sph_zero_spurious_modes!(A::AbstractMatrix)
    M, N = size(A)
    n = N÷2
    @inbounds for j = 1:n-1
        @simd for i = M-j+1:M
            A[i,2j] = 0
            A[i,2j+1] = 0
        end
    end
    @inbounds @simd for i = M-n+1:M
        A[i,2n] = 0
        2n < N && (A[i,2n+1] = 0)
    end
    A
end

function sphrand(::Type{T}, m::Int, n::Int) where T
    A = zeros(T, m, 2n-1)
    for i = 1:m
        A[i,1] = rand(T)
    end
    for j = 1:n-1
        for i = 1:m-j
            A[i,2j] = rand(T)
            A[i,2j+1] = rand(T)
        end
    end
    A
end

function sphrandn(::Type{T}, m::Int, n::Int) where T
    A = zeros(T, m, 2n-1)
    for i = 1:m
        A[i,1] = randn(T)
    end
    for j = 1:n-1
        for i = 1:m-j
            A[i,2j] = randn(T)
            A[i,2j+1] = randn(T)
        end
    end
    A
end

function sphones(::Type{T}, m::Int, n::Int) where T
    A = zeros(T, m, 2n-1)
    for i = 1:m
        A[i,1] = one(T)
    end
    for j = 1:n-1
        for i = 1:m-j
            A[i,2j] = one(T)
            A[i,2j+1] = one(T)
        end
    end
    A
end

sphzeros(::Type{T}, m::Int, n::Int) where T = zeros(T, m, 2n-1)

function normalizecolumns!(A::AbstractMatrix)
    m, n = size(A)
    @inbounds for j = 1:n
        nrm = zero(eltype(A))
        for i = 1:m
            nrm += abs2(A[i,j])
        end
        nrm = sqrt(nrm)
        for i = 1:m
            A[i,j] /= nrm
        end
    end
    A
end

function maxcolnorm(A::AbstractMatrix)
    m, n = size(A)
    nrm = zeros(n)
    @inbounds for j = 1:n
        nrm[j] = 0
        for i = 1:m
            nrm[j] += abs2(A[i,j])
        end
        nrm[j] = sqrt(nrm[j])
    end
    norm(nrm, Inf)
end

"""
Pointwise evaluation of real orthonormal spherical harmonic:

```math
Y_\\ell^m(\\theta,\\varphi) = (-1)^{|m|}\\sqrt{(\\ell+\\frac{1}{2})\\frac{(\\ell-|m|)!}{(\\ell+|m|)!}} P_\\ell^{|m|}(\\cos\\theta) \\sqrt{\\frac{2-\\delta_{m,0}}{2\\pi}} \\left\\{\\begin{array}{ccc} \\cos m\\varphi & {\\rm for} & m \\ge 0,\\\\ \\sin(-m\\varphi) & {\\rm for} & m < 0.\\end{array}\\right.
```
"""
sphevaluate(θ, φ, L, M) = sphevaluatepi(θ/π, φ/π, L, M)

sphevaluatepi(θ::Number, φ::Number, L::Integer, M::Integer) = sphevaluatepi(θ, L, M)*sphevaluatepi(φ, M)

function sphevaluatepi(θ::Number, L::Integer, M::Integer)
    ret = one(θ)/sqrt(two(θ))
    if M < 0 M = -M end
    c, s = cospi(θ), sinpi(θ)
    for m = 1:M
        ret *= sqrt((m+half(θ))/m)*s
    end
    tc = two(c)*c

    if L == M
        return ret
    elseif L == M+1
        return sqrt(two(θ)*M+3)*c*ret
    else
        temp = ret
        ret *= sqrt(two(θ)*M+3)*c
        for l = M+1:L-1
            ret, temp = (sqrt(l+half(θ))*tc*ret - sqrt((l-M)*(l+M)/(l-half(θ)))*temp)/sqrt((l-M+1)*(l+M+1)/(l+3half(θ))), ret
        end
        return ret
    end
end

sphevaluatepi(φ::Number, M::Integer) = sqrt((two(φ)-δ(M, 0))/(two(φ)*π))*(M ≥ 0 ? cospi(M*φ) : sinpi(-M*φ))


text:
module ReHorde

import LinearAlgebra: dot

using LoopVectorization

using Tullio

struct FeatureCumulant
    idx::Int
end
(c::FeatureCumulant)(o, x::Vector{Int}) = c.idx ∈ x

struct ObservationCumulant
    idx::Int
end
(c::ObservationCumulant)(o, x) = o[c.idx]

struct OnPolicy end

# What are we doing.
# Massively parallel GVF learning
struct GVF{W, Z, C, Π, Γ}
    w::W
    z::Z
    c::C
    π::Π
    γ::Γ
end

(gvf::GVF{<:AbstractVector})(x) = dot(gvf.w, x)
(gvf::GVF{<:AbstractVector})(x::AbstractVector{Int}) = begin; w = gvf.w; @tullio ret := w[x[i]]; end #sum(@view gvf.w[x])

struct TDλ
    α::Float32
    λ::Float32
end

function update!(lu::TDλ, gvf::GVF, o_t, a_t, o_tp1, x_t, x_tp1, r, μ_t)

    z = gvf.z
    w = gvf.w

    λ = lu.λ
    ρ_t = if gvf.π isa OnPolicy
        1f0
    else
        gvf.π(o_t, a_t)/μ_t
    end
    γ_t = if gvf.γ isa AbstractFloat
        gvf.γ
    else
        gvf.γ(o_t)
    end
    c = gvf.c(o_tp1, x_tp1)
    δ = c + γ_t*gvf(x_tp1) - gvf(x_t)

    if eltype(x_t) <: Integer
        z .*= γ_t*λ
        # view(z, x_t) .+= 1
        @tullio z[x_t[i]] = z[x_t[i]] + 1
        z .*= ρ_t
        w .+= (lu.α * δ) .* z
    else
        gvf.z .= ρ_t .* ((γ_t*λ) * gvf.z .+ x_t)
        gvf.w .+= lu.α * δ * gvf.z
    end
end

function update!(lu, gvfs::Vector{G}, o_t, a_t, o_tp1, x_t, x_tp1, r, μ_t) where G<:GVF
    Threads.@threads for i in 1:length(gvfs)
	update!(lu, gvfs[i], o_t, a_t, o_tp1, x_t, x_tp1, r, μ_t)
    end
end

# Write your package code here.
include("critterbot_utils.jl")
include("Critterbot.jl")

import ProgressMeter: @showprogress
const CBU = ReHorde.CritterbotUtils

function main_experiment()

    tiled_features_raw = CBU.loadTiles()
    tiled_features = [tiled_features_raw[i, :] for i in 1:size(tiled_features_raw, 1)]
    sensors_raw = Float32.(CBU.relevant_sensors())
    sensors = [sensors_raw[i, :] for i in 1:size(sensors_raw, 1)]

    horde = vcat([[
        [GVF(
	    zeros(Float32, CBU.numFeatures()),
	    zeros(Float32, CBU.numFeatures()),
	    ReHorde.ObservationCumulant(i),
	    ReHorde.OnPolicy(),
	    γ) for i in 1:length(sensors[1])];
        [GVF(
	    zeros(Float32, CBU.numFeatures()),
	    zeros(Float32, CBU.numFeatures()),
	    ReHorde.FeatureCumulant(i),
	    ReHorde.OnPolicy(),
	    γ) for i in rand(1:CBU.numFeatures(), 487)]]
                  for γ in [0.0, 0.8, 0.95, 0.9875]]...)

    @show length(horde)

    lu = TDλ(0.1/length(tiled_features[1]), 0.9)

    @showprogress 0.1 "Time Step: " for t = 2:length(tiled_features)
        x_t = tiled_features[t-1]
        x_tp1 = tiled_features[t]
        o_t = sensors[t-1]
        o_tp1 = sensors[t]
        a = nothing
        μ = nothing
        r = nothing

        update!(lu, horde, o_t, a, o_tp1, x_t, x_tp1, r, μ)
    end

end

end


text:
export save_as_dot, save_circuit, save_as_sdd, save_as_tex, save_as_dot2tex

#####################
# Save lines
#####################

function save_lines(name::String, lines::CircuitFormatLines)
    open(name, "w") do f
        for line in lines
            println(f, line)
        end
    end
end

#####################
# decompile for nodes
#####################

"Decompile for sdd circuit, used during saving of circuits to file"
decompile(n::Union{SddLiteralNode,PlainStructLiteralNode}, node2id, vtree2id) =
    UnweightedLiteralLine(node2id[n], vtree2id[n.vtree], literal(n), false)

decompile(n::Union{SddConstantNode,PlainStructConstantNode}, node2id, _) =
    AnonymousConstantLine(node2id[n], constant(n), false)

decompile(n::Union{Sdd⋁Node,PlainStruct⋁Node}, node2id, vtree2id) =
    DecisionLine(node2id[n], vtree2id[n.vtree], UInt32(num_children(n)),
                 map(c -> make_element(c, node2id), children(n)))

make_element(n::Union{Sdd⋀Node,PlainStruct⋀Node}, node2id) =
    SDDElement(node2id[children(n)[1]],  node2id[children(n)[2]])

make_element(_::StructLogicCircuit, _) =
    error("Given circuit is not an SDD, its decision node elements are not conjunctions.")

# TODO: decompile for logic circuit to some file format

#####################
# build maping
#####################

function get_node2id(circuit::LogicCircuit)
    node2id = Dict{LogicCircuit, ID}()
    outnodes = filter(n -> !is⋀gate(n), circuit)
    sizehint!(node2id, length(outnodes))
    index = ID(0) # node id start from 0
    for n in outnodes
        node2id[n] = index
        index += ID(1)
    end
    node2id
end

function get_vtree2id(vtree::Vtree):: Dict{Vtree, ID}
    vtree2id = Dict{Vtree, ID}()
    sizehint!(vtree2id, num_nodes(vtree))
    index = ID(0) # vtree id start from 0
    foreach(vtree) do n
        vtree2id[n] = index
        index += ID(1)
    end
    vtree2id
end

#####################
# saver for circuits
#####################

"Returns header for SDD file format"
function sdd_header()
    """
    c ids of sdd nodes start at 0
    c sdd nodes appear bottom-up, children before parents
    c
    c file syntax:
    c sdd count-of-sdd-nodes
    c F id-of-false-sdd-node
    c T id-of-true-sdd-node
    c L id-of-literal-sdd-node id-of-vtree literal
    c D id-of-decomposition-sdd-node id-of-vtree number-of-elements {id-of-prime id-of-sub}*
    c
    c File generated by Juice.jl
    c"""
end

"""
    save_as_sdd(name::String, circuit::LogicCircuit, vtree)

Save a SDD circuit to file. File name should end with '.sdd'.
"""
function save_as_sdd(name::String, circuit, vtree)
    @assert endswith(name, ".sdd")
    node2id = get_node2id(circuit)
    vtree2id = get_vtree2id(vtree)
    formatlines = Vector{CircuitFormatLine}()
    append!(formatlines, parse_sdd_file(IOBuffer(sdd_header())))
    push!(formatlines, SddHeaderLine(num_nodes(circuit)))
    for n in filter(n -> !is⋀gate(n), circuit)
        push!(formatlines, decompile(n, node2id, vtree2id))
    end
    save_lines(name, formatlines)
end

"""
    save_circuit(name::String, circuit, vtree)

Save a circuit to file. Currently defaults to and only supports saving '.sdd' format.
"""
save_circuit(name::String, circuit, vtree) =
    save_as_sdd(name, circuit, vtree)

"Rank nodes in the same layer left to right"
function get_nodes_level(circuit::Node)
    levels = Vector{Vector{Node}}()
    current = Vector{Node}()
    next = Vector{Node}()

    push!(next, circuit)
    push!(levels, Base.copy(next))
    while !isempty(next)
        current, next = next, current
        while !isempty(current)
            n = popfirst!(current)
            if isinner(n)
                for c in children(n)
                    if !(c in next) push!(next, c); end
                end
            end
        end
        push!(levels, Base.copy(next))
    end

    return levels
end

"""Effectively returns the inverse of `get_nodes_level`, returning a dictionary where keys are the
nodes and values are their levels, and the highest level in the circuit."""
function get_level_nodes(circuit::LogicCircuit)::Tuple{Dict{LogicCircuit, Int}, Int}
    levels = Dict{LogicCircuit, Int}()
    current = Vector{Node}()
    next = Vector{Node}()

    push!(next, circuit)
    i = 1
    while !isempty(next)
        current, next = next, current
        while !isempty(current)
            n = popfirst!(current)
            if isinner(n)
                for c in children(n)
                    if c ∉ next
                        push!(next, c)
                        levels[c] = i
                    end
                end
            end
        end
        i += 1
    end

    return levels, i
end

"""
    save_as_dot(file::String, circuit::LogicCircuit)

Save logic circuit in 'dot' file format.
"""
function save_as_dot(file::String, circuit::LogicCircuit)
    circuit_nodes = linearize(circuit)
    node_cache = Dict{LogicCircuit, Int64}()
    for (i, n) in enumerate(circuit_nodes)
        node_cache[n] = i
    end

    levels = get_nodes_level(circuit)

    f = open(file, "w")
    write(f,"digraph Circuit {\nsplines=false\nedge[arrowhead=\"none\",fontsize=6]\n")

    for level in levels
        if length(level) > 1
            write(f,"{rank=\"same\";newrank=\"true\";rankdir=\"LR\";")
            rank = ""
            foreach(x->rank*="$(node_cache[x])->",level)
            rank = rank[1:end-2]
            write(f, rank)
            write(f,"[style=invis]}\n")
        end
    end

    for n in reverse(circuit_nodes)
        if is⋀gate(n)
            write(f, "$(node_cache[n]) [label=\"*$(node_cache[n])\"]\n")
        elseif is⋁gate(n)
            write(f, "$(node_cache[n]) [label=\"+$(node_cache[n])\"]\n")
        elseif isliteralgate(n) && ispositive(n)
            write(f, "$(node_cache[n]) [label=\"+$(variable(n))\"]\n")
        elseif isliteralgate(n)  && isnegative(n)
            write(f, "$(node_cache[n]) [label=\"-$(variable(n))\"]\n")
        elseif isfalse(n)
            write(f, "$(node_cache[n]) [label=\"F\"]\n")
        elseif istrue(n)
            write(f, "$(node_cache[n]) [label=\"T\"]\n")
        else
            throw("unknown node type")
        end
    end

    for n in reverse(circuit_nodes)
        if isinnergate(n)
            for c in children(n)
                write(f, "$(node_cache[n]) -> $(node_cache[c])\n")
            end
        end
    end

    write(f, "}\n")
    flush(f)
    close(f)
end

"""
    save_as_dot2tex(file::String, circuit::LogicCircuit; **other_args)

Save logic circuit as a LaTeX TikZ .tex file using the dot2tex engine.

Argument `V` is a label map for each vertex; `E` is a label map for each edge; `⋀_style`,
`⋁_style`, `lit_style`, `⊤_style`, and `⊥_style` are styles for each node type; `edge_turn` is how
far down should edges "bend"; `node_xsep` and `node_ysep` are x and y-axis separation sizes for
terminal nodes.

Note: `save_as_tex` does not work well on large circuits, and the resulting LaTeX file should serve
as a first sketch in (serious) need for adjustments.
"""
function save_as_dot2tex(file::String, circuit::LogicCircuit,
                     V::Dict{LogicCircuit, Any} = Dict{LogicCircuit, Any}(),
                     E::Dict{Tuple{LogicCircuit, Int}, Any} = Dict{Tuple{LogicCircuit, Int}, Any}(),
                     ⋀_style::String = "and gate,fill=blue!50!red!30",
                     ⋁_style::String = "or gate,fill=blue!50!green!30",
                     ⊤_style::String = "", ⊥_style::String = "", lit_style::String = "",
                     edge_turn::Float64 = 0.25, ortho_splines::Bool = true,
                     node_xsep::Float64 = 1.0, node_ysep::Float64 = 1.25)
    preamble = """
    \\documentclass{standalone}

    \\usepackage{mathtools}
    \\usepackage{amssymb}
    \\usepackage{amsfonts}
    \\usepackage{tikz}
    \\usepackage{xcolor}
    \\usepackage{dot2texi}
    \\usetikzlibrary{shapes,arrows,positioning,fit,circuits.logic.US}

    \\begin{document}

    \\begin{tikzpicture}[>=latex',circuit logic US,every circuit symbol/.style={thick,point up}]
        \\tikzstyle{and} = [$(⋀_style)]
        \\tikzstyle{or} = [$(⋁_style)]
        \\tikzstyle{top} = [$(⊤_style)]
        \\tikzstyle{bot} = [$(⊥_style)]
        \\tikzstyle{lit} = [$(lit_style)]
        \\begin{dot2tex}[dot,tikz,codeonly,styleonly,options=-s -tmath]
            digraph G {
                graph [pad="0.75", ranksep="0.25", nodesep="0.25"];
                node [label=""];
                $(ortho_splines ? "edge [style=\"invis\"];" : "")
    """
    postamble = """
    \\end{tikzpicture}
    \\end{document}
    """

    nodes = linearize(circuit)
    cache = Dict{LogicCircuit, Int}()
    for (i, n) in enumerate(nodes) cache[n] = i end
    L = get_nodes_level(circuit)

    function label(v::LogicCircuit)::String
        if haskey(V, v)
            return string(V[v])
        elseif isleaf(v)
            if isliteralgate(v)
                l = literal(v)
                return l < 0 ? "\$\\neg $(-l)\$" : "\$$l\$"
            end
            return istrue(v) ? "\$\\top\$" : "\$\\bot\$"
        end
        return ""
    end
    label(v::LogicCircuit, c::Int)::String = (e = (v, c); haskey(E, e) ? string(E[e]) : "")
    compute_xsep(m::Int, i::Int, sep::Float64)::String = string(sep*(i-1-((m-1)/2)))
    compute_edge_turn(m::Int, i::Int)::String = string(-edge_turn*(0.75*m-abs((m+1)/2 - i)))
    draw_edge(pid::Int, cid::Int, c::Int, m::Int, l::String)::String =
        "    \\draw (n$pid.input $c) -- ++(0,$(compute_edge_turn(m, c))) -| $l (n$cid);\n"
    leaf_type(t::LogicCircuit)::String = istrue(t) ? "top" : isfalse(t) ? "bot" : "lit"
    leaf_tikz(id::Int, pid::Int, m::Int, c::Int, label::String)::String =
        "    \\node (n$id) at (\$(n$pid.input $c) + ($(compute_xsep(m, c, node_xsep/2)),-$(node_ysep/2))\$) {$label};\n"

    f = open(file, "w")
    write(f, preamble)

    # Define nodes, except leaves which are treated as a special case.
    for N in Iterators.reverse(nodes)
        if isinner(N)
            write(f, "            n$(cache[N]) [style=\"$(is⋀gate(N) ? "and" : "or"),inputs=$('n'^num_children(N))\", label=\"$(label(N))\"];\n")
        end
    end

    # Write ranks.
    for l in L
        if length(l) > 1
            write(f, "            {rank=\"same\";newrank=\"true\";rankdir=\"LR\";")
            rank = ""
            foreach(x -> rank *= "n$(cache[x])->", l)
            rank = rank[1:end-2]
            write(f, rank)
            write(f, "[style=\"invis\"]}\n")
        end
    end

    # Draw "invisible" edges top-down to force layout.
    for N in Iterators.reverse(nodes)
        if isleaf(N) continue end
        for C in children(N)
            write(f, "            n$(cache[N]) -> n$(cache[C]);\n")
        end
    end

    write(f, "        }\n    \\end{dot2tex}\n")

    # Draw actual edges and terminal nodes.
    if ortho_splines
        for N in Iterators.reverse(nodes)
            if isleaf(N) continue end
            N_id = cache[N]
            m = num_children(N)
            for (i, C) in enumerate(children(N))
                C_id = cache[C]
                if isleaf(C) write(f, leaf_tikz(C_id, N_id, m, i, label(C))) end
                l = label(N, i)
                write(f, draw_edge(N_id, C_id, i, m, isempty(l) ? "" : "node {\\small $l}"))
            end
        end
    end

    write(f, postamble)
    flush(f)
    close(f)
end

"""
    save_as_tex(file::String, circuit::LogicCircuit; **other_args)

Save logic circuit as a LaTeX TikZ .tex file.
Argument `V` is a label map for each vertex, `E` is a label map for each edge, `⋀_style` is the
TikZ conjunction node style to use, `⋁_style` is the TikZ disjunction node style to use,
`node_ysep` is the y-axis separation distance for nodes, `node_xsep` for the x-axis, `edge_turn` is
how far down should edges "bend", `verbose` if set to true prints comments to the LaTeX file for
better readability.

Note: `save_as_tex` does not work well on large circuits, and the resulting LaTeX file should serve
as a first sketch in (serious) need for adjustments.
"""
function save_as_tex(file::String, circuit::LogicCircuit;
                     V::Dict{LogicCircuit, Any} = Dict{LogicCircuit, Any}(),
                     E::Dict{Tuple{LogicCircuit, Int}, Any} = Dict{Tuple{LogicCircuit, Int}, Any}(),
                     ⋀_style::String = "and gate,fill=blue!50!red!30",
                     ⋁_style::String = "or gate,fill=blue!50!green!30",
                     node_ysep::Float64 = 1.25, node_xsep::Float64 = 1.0, edge_turn::Float64 = 0.25,
                     verbose::Bool = false)
    preamble = """
    \\documentclass{standalone}

    \\usepackage{mathtools}
    \\usepackage{amssymb}
    \\usepackage{amsfonts}
    \\usepackage{tikz}
    \\usepackage{xcolor}
    \\usetikzlibrary{shapes,arrows,positioning,fit,circuits.logic.US}

    \\newcommand{\\andNode}[4]{\\node[style={thick},#4,$(⋀_style)] (#1) at #2 {\\rotatebox{-90}{#3}}}
    \\newcommand{\\orNode}[4]{\\node[style={thick},#4,$(⋁_style)] (#1) at #2 {\\rotatebox{-90}{#3}}}

    \\begin{document}

    \\begin{tikzpicture}[circuit logic US,every circuit symbol/.style={thick,point up}]
    """
    postamble = """
    \\end{tikzpicture}
    \\end{document}
    """

    nodes = linearize(circuit)
    cache = Dict{LogicCircuit, Int}()
    for (i, n) in enumerate(nodes) cache[n] = i end
    L, max_level = get_level_nodes(circuit)

    function label(v::LogicCircuit)::String
        if haskey(V, v)
            return string(V[v])
        elseif isleaf(v)
            if isliteralgate(v)
                l = literal(v)
                return l < 0 ? "\$\\neg $(-l)\$" : "\$$l\$"
            end
            return istrue(v) ? "\$\\top\$" : "\$\\bot\$"
        end
        return ""
    end
    label(v::LogicCircuit, c::Int)::String = (e = (v, c); haskey(E, e) ? string(E[e]) : "")

    compute_xsep(m::Int, i::Int, sep::Float64, l::Int)::String = string(((max_level-l)*sep)*(i-1-((m-1)/2)))
    compute_edge_turn(m::Int, i::Int)::String = string(-edge_turn*(0.75*m-abs((m+1)/2 - i)))

    # Draws a TikZ root logical node. Arguments are: node type `t`, number of children `n`,
    # node ID `id`, and node label.
    root_tikz(t::String, n::Int, id::Int, label::String)::String =
        "  \\$(t)Node{n$id}{(0,0)}{$label}{inputs=$('n'^n)};\n"

    # Draws a TikZ inner logical node. Arguments are: node type `t`, number of children `n`, node
    # ID `id`, parent node ID `pid`, number of siblings `m`, child index `c`, node level `l`, and
    # node label.
    inner_tikz(t::String, n::Int, id::Int, pid::Int, m::Int, c::Int, l::Int, label::String)::String =
        "  \\$(t)Node{n$id}{(\$(n$pid.input $c) + ($(compute_xsep(m, c, node_xsep, l)),-$node_ysep)\$)}{$label}{inputs=$('n'^n)};\n"

    # Draws a TikZ leaf node. Arguments are: node ID `id`, parent node ID `pid`, number of siblings
    # `n`, child index `c`, and a terminal label.
    leaf_tikz(id::Int, pid::Int, m::Int, c::Int, label::String)::String =
        "  \\node (n$id) at (\$(n$pid.input $c) + ($(compute_xsep(m, c, node_xsep/2, max_level-1)),-$(node_ysep/2))\$) {$label};\n"

    # Draws a TikZ edge. Arguments are: outgoing edge node ID `pid`, incoming edge node ID `cid`,
    # number of "sibling" edges `m`, child index `c`, and edge label.
    edge_tikz(pid::Int, cid::Int, m::Int, c::Int, label::String)::String =
        (l = isempty(label) ? "" : "node {\\small $label}";
         "  \\draw (n$pid.input $c) -- ++(0,$(compute_edge_turn(m, c))) -| $l (n$cid);\n")

    inner_rep(N::LogicCircuit)::String = is⋀gate(N) ? "and" : "or"

    f = open(file, "w")

    write(f, preamble)

    Q = LogicCircuit[circuit]
    U = Dict{LogicCircuit, Bool}(circuit => true)
    first = true
    while !isempty(Q)
        N = popfirst!(Q)
        N_id = cache[N]
        m = num_children(N)
        if first
            verbose && write(f, "  % Root of circuit with ID: $N_id\n")
            write(f, root_tikz(inner_rep(N), m, N_id, label(N)))
            first = false
        end
        for (i, C) in enumerate(children(N))
            if haskey(U, C) continue end
            C_id = cache[C]
            if isleaf(C)
                verbose && write(f, "  % Leaf node with ID: $C_id, child of $N_id\n")
                write(f, leaf_tikz(C_id, N_id, m, i, label(C)))
            else
                t = inner_rep(C)
                verbose && write(f, "  % Inner node ($t gate) with ID: $C_id, child of $N_id\n")
                write(f, inner_tikz(t, num_children(C), C_id, N_id, m, i, L[C], label(C)))
                push!(Q, C); U[C] = true
            end
            verbose && write(f, "  % Edge from $N_id to $C_id\n")
            write(f, edge_tikz(N_id, C_id, m, i, label(N, i)))
        end
    end

    write(f, postamble)
    flush(f)
    close(f)
end


text:
searchdir(path, key) = filter(x -> occursin(key, x), readdir(path))
checkpoints = searchdir(pwd(), filename_1 * "_checkpoint_iteration")
if length(checkpoints) > 0
    checkpointed = true
    checkpoint_filepath = joinpath(pwd(), checkpoints[end])
    @info "Restoring from checkpoint: $checkpoint_filepath"
    model = restore_from_checkpoint(checkpoint_filepath, boundary_conditions = bcs, forcing = forcings, advection = advection_scheme)
else
    checkpointed = false
	model = IncompressibleModel(
	           architecture = arch,
	             float_type = FT,
	                   grid = grid,
	               coriolis = coriolis,
	               buoyancy = buoyancy,
	                closure = closure,
	                tracers = (:b,),
        boundary_conditions = bcs,
                    forcing = forcings,
                   advection = advection_scheme,
                   timestepper = timestepping_scheme,
	)
end

if !checkpointed
    if !geostrophic_balance
        set!(model, b=B₀)
    else
        set!(model, u = U₀, b=B₀)
    end
end

if !checkpointed
	Δt= 120.0
else
	Δt = 300.0
end

max_stars_repo_name: sandreza/Mesoscale

max_issues_repo_name: sandreza/Mesoscale

max_forks_repo_name: sandreza/Mesoscale


text:
module ILMs

using DataFrames, Distances, Distributions, Gadfly, Interact, Mamba

export
  create_event_db,
  create_pop_db_univariates,
  create_pop_db_bivariate,
  infect_recover_loop,
  state_timeseries,
  state_animation,
  spatial_states,
  create_loglikelihood

include("population.jl")
include("events.jl")
include("visuals.jl")
include("inference.jl")

end # module


text:
@testset "disrelation.jl" begin
    x, A = demodata()
    Ã = A .- mean(A, dims=1)

    Λ = disrelation(Ã)
    @plottest begin
        heatmap(x, x, Λ)
        plot!(xflip=true, yflip=true)
        xlabel!("Peak position")
        ylabel!("Peak position")
        title!("Disrelation")
    end "disrelation.png"

    # Can this be prettier?
    @test .≈(diag(Λ), 0) |> all

    Ψ = asynchronous(Ã)
    @test maximum(sqrt.((Λ .- Ψ).^2)) ≈ 0 atol=2e-3
end


text:
# This file includes the plugin for calculation of fast fourier transform of data
using FFTW

"""
    Fft(dims::Int)

Constructs an `Fft` plugin. The  [`process(plg::Fft, x)`](@ref) function performes an `fft` operatinon along `dims` of `x`. See also: [`fft`](https://juliamath.github.io/AbstractFFTs.jl/stable/api/#AbstractFFTs.fft)
"""
struct Fft <: AbstractPlugin
    dims::Int
end
Fft(;dims::Int=1) = Fft(dims)

show(io::IO, plg::Fft) = print(io, "Fft(dims:$(plg.dims))")

"""
    process(plg::Fft, x)

Performes an `fft` transformation for the input data `x`.

# Example
```julia
julia> x = collect(reshape(1:16, 4,4))
4×4 Array{Int64,2}:
 1  5   9  13
 2  6  10  14
 3  7  11  15
 4  8  12  16

julia> plg = Plugins.Fft(dims=1)
Fft(dims:1)

julia> process(plg, x)
4×4 Array{Complex{Float64},2}:
 10.0+0.0im  26.0+0.0im  42.0+0.0im  58.0+0.0im
 -2.0+2.0im  -2.0+2.0im  -2.0+2.0im  -2.0+2.0im
 -2.0+0.0im  -2.0+0.0im  -2.0+0.0im  -2.0+0.0im
 -2.0-2.0im  -2.0-2.0im  -2.0-2.0im  -2.0-2.0im
```
"""
process(plg::Fft, x) = fft(x, plg.dims)


text:
println("Testing crossover function.")
using SymbolicRegression
using Test
using SymbolicRegression: crossover_trees
options = SymbolicRegression.Options(;
    default_params...,
    binary_operators=(+, *, /, -),
    unary_operators=(cos, exp),
    npopulations=8,
)
tree1 = cos(Node("x1")) + (3.0f0 + Node("x2"))
tree2 = exp(Node("x1") - Node("x2") * Node("x2")) + 10.0f0 * Node("x3")

# See if we can observe operators flipping sides:
cos_flip_to_tree2 = false
exp_flip_to_tree1 = false
swapped_cos_with_exp = false
for i in 1:1000
    child_tree1, child_tree2 = crossover_trees(tree1, tree2)
    if occursin("cos", repr(child_tree2))
        # Moved cosine to tree2
        global cos_flip_to_tree2 = true
    end
    if occursin("exp", repr(child_tree1))
        # Moved exp to tree1
        global exp_flip_to_tree1 = true
    end
    if occursin("cos", repr(child_tree2)) && occursin("exp", repr(child_tree1))
        global swapped_cos_with_exp = true
        # Moved exp with cos
        @assert !occursin("cos", repr(child_tree1))
        @assert !occursin("exp", repr(child_tree2))
    end

    # Check that exact same operators, variables, numbers before and after:
    rep_tree_final = sort([a for a in repr(child_tree1) * repr(child_tree2)])
    rep_tree_final = strip(String(rep_tree_final), ['(', ')', ' '])
    rep_tree_initial = sort([a for a in repr(tree1) * repr(tree2)])
    rep_tree_initial = strip(String(rep_tree_initial), ['(', ')', ' '])
    @test rep_tree_final == rep_tree_initial
end

@test cos_flip_to_tree2
@test exp_flip_to_tree1
@test swapped_cos_with_exp
println("Passed.")


text:
using Clang.cindex
using Clang.wrap_c

JULIAHOME=EnvHash()["JULIAHOME"]

clang_includes = map(x->joinpath(JULIAHOME, x), [
  "deps/llvm-3.2/build/Release/lib/clang/3.2/include",
  "deps/llvm-3.2/include",
  "deps/llvm-3.2/include",
  "deps/llvm-3.2/build/include/",
  "deps/llvm-3.2/include/"
  ])
clang_extraargs = ["-D", "__STDC_LIMIT_MACROS", "-D", "__STDC_CONSTANT_MACROS"]

wrap_hdrs = map( x-> joinpath("/usr/include/curl", x), [ "curl.h", "easy.h", "multi.h" ])

wc = wrap_c.init(".", "../src/lC_common_h.jl", clang_includes, clang_extraargs, (th, h) -> contains(wrap_hdrs, h) , h -> "libcurl", h -> "../src/lC_" * replace(last(split(h, "/")), ".", "_")  * ".jl" )
wc.options.wrap_structs = true

wrap_c.wrap_c_headers(wc, ["/usr/include/curl/curl.h"])

# generate export statements.....
fe = open("../src/lC_exports_h.jl", "w+")
println(fe, "#   Generating exports")

fc = open("../src/lC_curl_h.jl", "r")
curljl = split(readall(fc), "\n")
close(fc)

for e in curljl
  m = match(r"^\s*\@c\s+[\w\:\{\}\_]+\s+(\w+)", e)
  if (m != nothing)
#    println (m.captures[1])
    @printf fe "export %s\n"  m.captures[1]
  end
end

fc = open("../src/lC_common_h.jl", "r")
curljl = split(readall(fc), "\n")
close(fc)

for e in curljl
  m = match(r"^\s*\@ctypedef\s+(\w+)", e)
  if (m != nothing)
#    println(m.captures[1])
    @printf fe "export %s\n"  m.captures[1]
  else
    m = match(r"^\s*const\s+(\w+)", e)
    if (m != nothing)
#        println(m.captures[1])
        @printf fe "export %s\n"  m.captures[1]
    end
  end
end

# #defines generated
const ign_defs = [
  "_PULL_SYS_TYPES_H",
  "_PULL_STDINT_H",
  "_PULL_INTTYPES_H",
  "_PULL_SYS_SOCKET_H",
  "_SIZEOF_LONG",
  "_TYPEOF_CURL_SOCKLEN_T",
  "_SIZEOF_CURL_SOCKLEN_T",
  "_TYPEOF_CURL_OFF_T",
  "_FORMAT_CURL_OFF_T",
  "_FORMAT_CURL_OFF_TU",
  "_FORMAT_OFF_T",
  "_SIZEOF_CURL_OFF_T",
  "_SUFFIX_CURL_OFF_T",
  "_SUFFIX_CURL_OFF_TU"
]

f = open("../src/lC_defines_h.jl", "w+")
println(f, "#   Generating #define constants")

hashdefs = split(readall(`gcc -E -dD -P /usr/include/curl/curl.h`), "\n")
for e in hashdefs
  m = match(r"^\s*#define\s+CURL(\w+)\s+(.+)", e)
  if (m != nothing) && !contains(ign_defs, strip(m.captures[1]))
    c2 = replace (m.captures[2], "(unsigned long)", "")
    @printf f "const CURL%-30s = %s\n"  m.captures[1]  c2
    @printf fe "export CURL%s\n"  m.captures[1]
  end
end

close(f)
close(fe)

# Generate the File extension to MIME type mapping.
f = open("../src/Mime_ext.jl", "w+")
write(f, "module Mime_ext\n\nexport MimeExt\n\nMimeExt = {\n")

mimes = split(open(readall, "/etc/mime.types"), "\n")
for e in mimes
  m = match(r"^\s*#", e)
  if (m == nothing)
    tokens = split(e)
    if (length(tokens) > 1)
        for i = 2:length(tokens)
            write (f, "  \"$(tokens[i])\" => \"$(tokens[1])\",\n")
        end
    end
  end
end

write(f, "  \"\" => \"\"\n}\n\nend\n")
close(f)


text:
export ScannerCoordinateSystem, toScannerCoords
abstract type DeviceCoords{T<:Unitful.Length} <: AbstractVector{T} end

struct ScannerCoords{T<:Unitful.Length} <: AbstractVector{T}
  data::Vector{T}
end

struct ScannerCoordinateSystem
  """Matrix defining the axes of the scanner coordinate system, the first column is the scanner x axis in robot coordinates"""
  axes::Matrix{Float64}
  """The origin should be set, so that the objGeometry is centered in the crosssection of the bore, the zero of the bore axis is not important, as the available boreAxis range is set by minMaxBoreAxis"""
  origin::Vector{typeof(1.0u"mm")}
  function ScannerCoordinateSystem(dof::Integer)
    return new(Matrix(1.0LinearAlgebra.I, dof, dof), zeros(dof)*u"mm")
  end
  function ScannerCoordinateSystem(axes::Matrix{<:Real}, origin::Vector{<:Unitful.Length})
    dim = length(origin)
    @assert all(isequal(dim), size(axes))
    return new(axes, origin)
  end
end

Base.Broadcast.broadcastable(val::ScannerCoordinateSystem) = Ref(val) # Otherwise tryuparse.(val) returns [val]

ScannerCoordinateSystem(axes::Matrix{<:Real}, origin::Nothing) = ScannerCoordinateSystem(axes, zeros(size(axes, 1))*u"mm")
ScannerCoordinateSystem(axes::Nothing, origin::Vector{<:Unitful.Length}) = ScannerCoordinateSystem(Matrix(1.0LinearAlgebra.I, length(origin), length(origin)), origin)
function ScannerCoordinateSystem(axes::AbstractString, origin)
  axes_split = split(lowercase(axes),",")
  cs_matrix = zeros(length(axes_split),length(axes_split))
  for (i,ax) in enumerate(axes_split)
    cs_matrix[i,ax[end]-'w'] = (ax[1]=='-' ? -1. : 1.)
  end
  return ScannerCoordinateSystem(cs_matrix, origin)
end
ScannerCoordinateSystem(axes::AbstractVector{<:AbstractVector{<:Real}}, origin) = ScannerCoordinateSystem(Matrix{Float64}(hcat(axes...)), origin)

function toScannerCoords(sys::ScannerCoordinateSystem, coords::DeviceCoords)
  shifted = coords.data - sys.origin
  return ScannerCoords(sys.axes * shifted)
end
toScannerCoords(sys::ScannerCoordinateSystem, coords::ScannerCoords) = coords

function toDeviceCoords(sys::ScannerCoordinateSystem, coords::ScannerCoords, T::Type{<:DeviceCoords})
  rotated = inv(sys.axes) * coords.data
  return T(rotated + sys.origin)
end


text:
function round_messages{T,F,R,K}(messages::Vector{Float64},
                        S::SparseMatrixCSC{T,Int64},w::Vector{F},
                        alpha::R,beta::K,
                        rp::Vector{Int64},ci::Vector{Int64},
                        tripi::Vector{Int64},n::Int64,
                        m::Int64,perm::Vector{Int64},
                        li::Vector{Int64},lj::Vector{Int64})
  ai = zeros(Float64,length(tripi))
  ai[tripi.>0] = messages[perm]
  M_output = MatrixNetworks.bipartite_matching_primal_dual(rp,ci,ai,m,n)
  mi = MatrixNetworks.edge_indicator(M_output,li,lj)
  matchweight = sum(w[mi.>0])#M_output.weight
  cardinality = sum(mi) #M_output.cardinality
  overlap = dot(mi,(S*mi)/2)
  f = alpha*matchweight + beta*overlap
  info = [f matchweight cardinality overlap]
  return info
end

max_stars_repo_name: nassarhuda/NetworkAlign

max_issues_repo_name: nassarhuda/NetworkAlign

max_forks_repo_name: nassarhuda/NetworkAlign


text:
abstract type AbstractDesign end

# make designs iterable
Base.iterate(design::AbstractDesign, state = 0) = state > 0 ? nothing : (design, state + 1)
Base.length(design::AbstractDesign) = 1

Base.show(io::IO, design::AbstractDesign) = print(io, string(design))
Base.show(io::IO, ::MIME"application/prs.juno.inline", design::AbstractDesign) = print(io, string(design))

function string(design::AbstractDesign)
    nn1         = n1(design)
    cont_region = Int(max(0, early_futility(design) + 1)):Int(min(nn1, early_efficacy(design) - 1))
    if length(cont_region) == 0 # one stage design
        return @sprintf "%s<n=%i;c=%i>" string(typeof(design)) nn1 early_futility(design)
    end
    n2_cont     = n2.(design, cont_region)
    if minimum(n2_cont) == maximum(n2_cont) # group sequential design
        return @sprintf "%s<n1=%i;n2:[%i,%i]->%i>" string(typeof(design)) nn1 cont_region[1] cont_region[end] n2_cont[1]
    end
    # generic two stage design
    return @sprintf "%s<n1=%i;n2:[%i,%i]->[%i,%i]>" string(typeof(design)) nn1 cont_region[1] cont_region[end] minimum(n2_cont) maximum(n2_cont)
end

n1(design::AbstractDesign)             = length(design.n2) - 1
n2(design::AbstractDesign, x1::Int)    = valid(design, x1) ? design.n2[x1 + 1] : error("0 <= x1 <= n1 violated")
n(design::AbstractDesign, x1::Int)     = n1(design) + n2(design, x1)
n(design::AbstractDesign)              = n1(design) .+ design.n2

early_futility(design::AbstractDesign) = any(design.c2 .== Inf) ? findlast(design.c2 .== Inf) - 1 : -Inf
futility_region(design::AbstractDesign) = collect(0:n1(design))[(n2.(design, 0:n1(design)) .== 0) .& (c2.(design, 0:n1(design)) .== Inf)]
early_efficacy(design::AbstractDesign) = any(design.c2 .== -Inf) ? findfirst(design.c2 .== -Inf) - 1 : Inf
efficacy_region(design::AbstractDesign) = collect(0:n1(design))[(n2.(design, 0:n1(design)) .== 0) .& (c2.(design, 0:n1(design)) .== -Inf)]
continuation_region(design::AbstractDesign) = collect(0:n1(design))[n2.(design, 0:n1(design)) .> 0]
early_stop_region(design::AbstractDesign) = vcat(futility_region(design), efficacy_region(design))

c2(design::AbstractDesign, x1::Int)    = valid(design, x1) ? design.c2[x1 + 1] : error("0 <= x1 <= n1 violated")

function as_table(design::AbstractDesign)
    DataFrames.DataFrame(
        n1 = repeat([n1(design)], n1(design) + 1),
        x1 = 0:n1(design),
      phat = (0:n1(design)) ./ repeat([n1(design)], n1(design) + 1),
        n2 = design.n2,
        n  = n(design),
        c2 = design.c2
    )
end

function sample_space(design::TD) where {TD<:AbstractDesign}

    return [ [x1, x2] for x1 in 0:n1(design) for x2 in 0:n2(design, x1) ] |>
        x -> hcat(x...)' |>
        x -> convert(Array{eltype(x),2}, x)
end

function sample_space(design::TD, x1::TI) where {TI<:Integer,TD<:AbstractDesign}
    return collect(0:n2(design, x1))
end

reject(x2::Int, c2::Real) = x2 > c2

function reject(x1::Int, x2::Int, design::AbstractDesign)
    !valid(design, x1, x2) ? error("invalid x1 / x2 for given design") : nothing
    return reject(x2, c2(design, x1))
end

valid(design::AbstractDesign, x1::Int) = 0 <= x1 <= n1(design)
valid(design::AbstractDesign, x1::Int, x2::Int) = valid(design, x1) & (0 <= x2 <= n2(design, x1))

function fisher_information_integrand(p::TR, x::TI, n::TI)::TR where {TR<:Real,TI<:Integer}
    ( x/p - (n - x)/(1 - p) )^2
end

function fisher_information(p::T, design::TD)::T where {T<:Real,TD<:AbstractDesign}
    # 0 < p < 1 ? nothing : error("p must be in (0, 1)")
    (p ≈ 0) | (p ≈ 1) ? (return Inf) : nothing
    XX        = sample_space(design)
    x1, x2    = XX[:,1], XX[:,2]
    log_prob  = log.(
        pmf.(x2, n2.(design, x1), p) .*
        pmf.(x1, n1(design), p)
    )
    x = x1 .+ x2
    return sum( exp.( log_prob .+ log.(
        ( x./p .- (n.(design, x1) .- x)./(1 - p) ).^2
    ) ) )
end

mutable struct Design <: AbstractDesign
    n2::Vector{Int}
    c2::Vector{Real}
    function Design(n2::Vector{N2}, c2::Vector{C2}) where {N2<:Integer, C2<:Real}
        !all(0 .<= n2) ? error("n must be positive") : nothing
        # ToDo: check for ineffective stopping (n2 > 0 and early stopping)
        length(n2) != length(c2) ? error("n and c must be of equal length") : nothing
        new(n2, c2)
    end
end
Design(n2, c2) = Design(convert(Vector{Integer}, n2), convert(Vector{Real}, c2))


text:
module TestCountMin

using CountMin

if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

include("testutil.jl")

@testset "CountMinSketch tests" begin
    @testset "CountMinSketch empty constructor" begin
        cms = CountMinSketch()
        @test eltype(cms.sketch) == UInt8
        @test cms.sketch == zeros(UInt8, 1000, 4)
        @test cms.tables == 4
        @test cms.tablesize == 1000
    end

    @testset "CountMinSketch normal constructor" begin
        cms = CountMinSketch{UInt16}(4, 100)
        @test eltype(cms.sketch) == UInt16
        @test size(cms.sketch) == (100, 4)
        @test cms.tables == 4
        @test cms.tablesize == 100

        # table too small @test_throws Exception CountMinSketch{UInt8}(4,0)
        # too few tables
        @test_throws Exception CountMinSketch{UInt8}(0,100)
        @test_throws Exception CountMinSketch{UInt8}(100,100)
    end

    @testset "CountMinSketch size/eltype" begin
        cms = CountMinSketch{UInt16}(4, 100)
        @test eltype(cms.sketch) == UInt16
        @test eltype(cms) == UInt16
        @test size(cms.sketch) == (100, 4)
        @test size(cms) == (4, 100)
    end

    @testset "CountMinSketch push/pop/add" begin
        cms = CountMinSketch{UInt16}(4, 100)
        @test cms[hash(1)] == UInt16(0)
        push!(cms, hash(1))
        @test cms[hash(1)] == UInt16(1)
        push!(cms, hash(1))
        @test cms[hash(1)] == UInt16(2)
        pop!(cms, hash(1))
        @test cms[hash(1)] == UInt16(1)
        add!(cms, hash(1), 100)
        @test cms[hash(1)] == UInt16(101)
        add!(cms, hash(1), -10)
        @test cms[hash(1)] == UInt16(91)
    end

    @testset "CountMinSketch getindex" begin
        cms = CountMinSketch{UInt16}(4, 100)
        item = "ABCD"
        @test cms[item] == UInt16(0)
        push!(cms, item)
        @test cms[item] == UInt16(1)
        push!(cms, item)
        @test cms[item] == UInt16(2)
        pop!(cms, item)
        @test cms[item] == UInt16(1)
        add!(cms, item, 100)
        @test cms[item] == UInt16(101)
        add!(cms, item, -10)
        @test cms[item] == UInt16(91)
    end

    @testset "CountMinSketch over/underflow" begin
        cms = CountMinSketch{UInt8}(4, 100)
        item = "ABCD"

        @test cms[item] == 0
        add!(cms, item, 256)
        @test cms[item] == 255
        add!(cms, item, -256)
        @test cms[item] == 0
    end

    @testset "CountMinSketch Read/Write" begin
        intempdir() do
            cms = CountMinSketch{UInt16}(4, 100)
            item = "ABCD"

            # Add item
            @test cms[item] == UInt16(0)
            push!(cms, item)
            @test cms[item] == UInt16(1)

            # Save and reload
            writecms("cms.h5", cms)

            newcms = readcms("cms.h5")
            # Shape & Type
            @test size(cms) == size(newcms)
            @test eltype(cms) == eltype(newcms)

            # Ensure value is preserved
            @test newcms[item] == UInt16(1)
        end
    end

    @testset "collisionrate(CountMinSketch)" begin
        cms = CountMinSketch{UInt16}(4, 100)

        # Fill the first half of each table. This is super dodgy
        for table in 1:4
            for i in 1:50
                cms.sketch[i, table] = 1
            end
        end
        # The above is to satisfy the following
        # @test sum(cms.sketch, 1) == UInt64[50 50 50 50]

        @test collisionrate(cms) == 0.5^4
    end
end # testset CountMinSketch

end # module TestCountMinSketch


text:
using BiomolecularStructures.WebBLAST
using BiomolecularStructures.Mafft
using Requests
using BiomolecularStructures.PDB
# 1. search for S_0 in PDB

#seq = "MNQLQQLQNPGESPPVHPFVAPLSYLLGTWRGQGEGEYPTIPSFRYGEEIRFSHSGKPVIAY"

#=seq = "MKKQLKYCFF SLFVSLSSIL SSCGSTTFVL ANFESYISPL LLERVQEKHP LTFLTYPSNE
KLINGFANNT YSVAVASTYA VSELIERDLL SPIDWSQFNL KKSSSSSDKV NNASDAKDLF
IDSIKEISQQ TKDSKNNELL HWAVPYFLQN LVFVYRGEKI SELEQENVSW TDVIKAIVKH
KDRFNDNRLV FIDDARTIFS LANIVNTNNN SADVNPKEDG IGYFTNVYES FQRLGLTKSN
LDSIFVNSDS NIVINELASG RRQGGIVYNG DAVYAALGGD LRDELSEEQI PDGNNFHIVQ
PKISPVALDL LVINKQQSNF QKEAHEIIFD LALDGADQTK EQLIKTDEEL GTDDEDFYLK
GAMQNFSYVN YVSPLKVISD PSTGIVSSKK NNAEMKSKQM STDQMTSEKE FDYYTETLKA
LLEKEDSAEL NENEKKLVET IKKAYTIEKD SSIRWNQLVE KPISPLQRSN LSLSWLDFKL
HWW" =#

seq = "QIKDLLVSSSTDLDTTLVLVNAIYFKGMWKTAFNAEDTREMPFHVTKQESKPVQMMCMNNSFNVATLPAE
KMKILELPFASGDLSMLVLLPDEVSDLERIEKTINFEKLTEWTNPNTMEKRRVKVYLPQMKIEEKYNLTS
VLMALGMTDLFIPSANLTGISSAESLKISQAVHGAFMELSEDGIEMAGSTGVIEDIKHSPESEQFRADHP
FLFLIKHNPTNTIVYFGRYWSP"

threshold = 0.005

results = webblast(seq, threshold, :ncbi, true)

# PDB IDs and Chain IDs
pdbs = (AbstractString,AbstractString)[]

fastastring = ""
for result in results
	fasta = fastarepresentation(result)

	accession = split(result.accession, "_")

	push!(pdbs, (convert(AbstractString,accession[1]), convert(AbstractString,accession[2])))

	fastastring = string(fastastring, ">", fasta[1], "\n", fasta[2], "\n")
end

# 2. Run MAFFT on results
if length(results) > 0
	print_fasta(mafft_from_string(fastastring))
end

# 3. Get PDBs

structures = Any[]

for pdb in pdbs
  push!(structures, get_structure(get_remote_pdb(pdb[1])))
end

# 4. Do something with the PDBs (Superimpose, etc)
for struc in structures
	m = structure_to_matrix(struc)
	println(size(m))
end


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Behavior.Gherkin
using Behavior:
    ResultAccumulator, accumulateresult!, issuccess, featureresults,
    FeatureResult, ScenarioResult, Given, SuccessfulStepExecution, Scenario, StepFailed

@testset "Result Accumulator   " begin
    @testset "Accumulate results; One feature with a successful scenario; Total result is success" begin
        accumulator = ResultAccumulator()

        feature = Gherkin.Feature(FeatureHeader("", [], []), Scenario[])
        given = Given("some precondition")
        successfulscenario = ScenarioResult(
            [SuccessfulStepExecution()],
            Scenario("Some scenario", String[], ScenarioStep[given]),
            Background(),
            ScenarioStep[])
        featureresult = FeatureResult(feature, [successfulscenario])

        accumulateresult!(accumulator, featureresult)

        @test issuccess(accumulator)
    end

    @testset "Accumulate results; One feature with a failing scenario; Total result is fail" begin
        accumulator = ResultAccumulator()

        feature = Gherkin.Feature(FeatureHeader("", [], []), Scenario[])
        given = Given("some precondition")
        scenario = ScenarioResult(
            [StepFailed("")],
            Scenario("Some scenario", String[], ScenarioStep[given]),
            Background(),
            ScenarioStep[])
        featureresult = FeatureResult(feature, [scenario])

        accumulateresult!(accumulator, featureresult)

        @test issuccess(accumulator) == false
    end

    @testset "Accumulate results; One scenario with one successful and one failing step; Total result is fail" begin
        accumulator = ResultAccumulator()

        feature = Gherkin.Feature(FeatureHeader("", [], []), Scenario[])
        given = Given("some precondition")
        when = When("some action")
        scenario = ScenarioResult(
            [SuccessfulStepExecution(), StepFailed("")],
            Scenario("Some scenario", String[], ScenarioStep[given, when]),
            Background(),
            ScenarioStep[])
        featureresult = FeatureResult(feature, [scenario])

        accumulateresult!(accumulator, featureresult)

        @test issuccess(accumulator) == false
    end

    @testset "Accumulate results; Two scenarios, one failing; Total result is fail" begin
        accumulator = ResultAccumulator()

        feature = Gherkin.Feature(FeatureHeader("", [], []), Scenario[])
        given = Given("some precondition")
        when = When("some action")
        successfulscenario = ScenarioResult(
            [SuccessfulStepExecution()],
            Scenario("Some scenario", String[], ScenarioStep[given]),
            Background(),
            ScenarioStep[])
        failingscenario = ScenarioResult(
            [StepFailed("")],
            Scenario("Some other scenario", String[],  ScenarioStep[when]),
            Background(),
            ScenarioStep[])
        featureresult = FeatureResult(feature, [successfulscenario, failingscenario])

        accumulateresult!(accumulator, featureresult)

        @test issuccess(accumulator) == false
    end

    @testset "Accumulate results; One feature, one passing Scenario; One success and zero failures" begin
        accumulator = ResultAccumulator()

        feature = Gherkin.Feature(FeatureHeader("", [], []), Scenario[])
        given = Given("some precondition")
        successfulscenario = ScenarioResult(
            [SuccessfulStepExecution()],
            Scenario("Some scenario", String[], ScenarioStep[given]),
            Background(),
            ScenarioStep[])
        featureresult = FeatureResult(feature, [successfulscenario])

        accumulateresult!(accumulator, featureresult)

        @test featureresults(accumulator)[1].n_success == 1
        @test featureresults(accumulator)[1].n_failure == 0
    end

    @testset "Accumulate results; Two scenarios, one failing; One success and one failure" begin
        accumulator = ResultAccumulator()

        feature = Gherkin.Feature(FeatureHeader("", [], []), Scenario[])
        given = Given("some precondition")
        when = When("some action")
        successfulscenario = ScenarioResult(
            [SuccessfulStepExecution()],
            Scenario("Some scenario", String[], ScenarioStep[given]),
            Background(),
            ScenarioStep[])
        failingscenario = ScenarioResult(
            [StepFailed("")],
            Scenario("Some other scenario", String[],  ScenarioStep[when]),
            Background(),
            ScenarioStep[])
        featureresult = FeatureResult(feature, [successfulscenario, failingscenario])

        accumulateresult!(accumulator, featureresult)

        @test featureresults(accumulator)[1].n_success == 1
        @test featureresults(accumulator)[1].n_failure == 1
    end

    @testset "Accumulate results; Seven scenarios, two failing; 5 success and 2 failures" begin
        accumulator = ResultAccumulator()

        feature = Gherkin.Feature(FeatureHeader("", [], []), Scenario[])
        given = Given("some precondition")
        when = When("some action")
        successfulscenario = ScenarioResult(
            [SuccessfulStepExecution()],
            Scenario("Some scenario", String[], ScenarioStep[given]),
            Background(),
            ScenarioStep[])
        failingscenario = ScenarioResult(
            [StepFailed("")],
            Scenario("Some other scenario", String[],  ScenarioStep[when]),
            Background(),
            ScenarioStep[])
        featureresult = FeatureResult(feature,
            [successfulscenario, successfulscenario, successfulscenario, successfulscenario, successfulscenario,
             failingscenario, failingscenario])

        accumulateresult!(accumulator, featureresult)

        @test featureresults(accumulator)[1].n_success == 5
        @test featureresults(accumulator)[1].n_failure == 2
    end

    @testset "Accumulate results; Two successful features; Both have no failures" begin
        accumulator = ResultAccumulator()

        feature1 = Gherkin.Feature(FeatureHeader("1", [], []), Scenario[])
        feature2 = Gherkin.Feature(FeatureHeader("1", [], []), Scenario[])
        given = Given("some precondition")
        successfulscenario = ScenarioResult(
            [SuccessfulStepExecution()],
            Scenario("Some scenario", String[], ScenarioStep[given]),
            Background(),
            ScenarioStep[])
        featureresult1 = FeatureResult(feature1, [successfulscenario])
        featureresult2 = FeatureResult(feature1, [successfulscenario])

        accumulateresult!(accumulator, featureresult1)
        accumulateresult!(accumulator, featureresult2)

        @test featureresults(accumulator)[1].n_success > 0
        @test featureresults(accumulator)[1].n_failure == 0
        @test featureresults(accumulator)[2].n_success > 0
        @test featureresults(accumulator)[2].n_failure == 0
    end

    @testset "Accumulate results; Failing Background with a successful Scenario; Total result is failed" begin
        accumulator = ResultAccumulator()

        feature = Gherkin.Feature(FeatureHeader("", [], []), Scenario[])
        given = Given("some precondition")
        successfulscenario = ScenarioResult(
            [SuccessfulStepExecution()],
            Scenario("Some scenario", String[], ScenarioStep[given]),
            Background("Failing background", ScenarioStep[Given("some failing background step")]),
            StepExecutionResult[
                StepFailed(""),
            ])
        featureresult = FeatureResult(feature, [successfulscenario])

        accumulateresult!(accumulator, featureresult)

        @test issuccess(accumulator) == false
    end

    @testset "Accumulate results; One feature with a successful scenario; Results accumulator is not empty" begin
        accumulator = ResultAccumulator()

        feature = Gherkin.Feature(FeatureHeader("", [], []), Scenario[])
        given = Given("some precondition")
        successfulscenario = ScenarioResult(
            [SuccessfulStepExecution()],
            Scenario("Some scenario", String[], ScenarioStep[given]),
            Background(),
            ScenarioStep[])
        featureresult = FeatureResult(feature, [successfulscenario])

        accumulateresult!(accumulator, featureresult)

        @test !isempty(accumulator)
    end

    @testset "Accumulate results; No features; Results accumulator is empty" begin
        accumulator = ResultAccumulator()

        @test isempty(accumulator)
    end

    @testset "Accumulate results; One feature with syntax error; Total result is failure" begin
        accumulator = ResultAccumulator()

        parseresult = Gherkin.BadParseResult{Feature}(:somereason, :someexpected, :someactual, 0, "Some line")

        accumulateresult!(accumulator, parseresult, "features/some/path/to/my.feature")

        @test issuccess(accumulator) == false
    end
end


text:
@testset "get_column_id" begin
    table = init_log_table(
        (id=:open_nodes, name="#Open"),
        (id=:closed_nodes, name="#Closed"),
    )
    TableLogger.get_column_id(table, :closed_nodes) == 2
end