text:
using .RationalFunctionFields: interpolate_rational_function,
        interpolate_multivariate_rational_function, polynomial_reconstruction,
        random_linear_shift, decompose_by_degrees, interpolate_polynomial,
        generate_kronecker_points, backward_kronecker,
        decompose_by_good_degrees, unpackfrac

function test_rational_function_interpolation(f)

    R = parent(denominator(f))

    # hope we have multivariate ring here
    vs = collect(AA.gens(R))
    nvs = length(vs)

    npoints = 100

    varidx = 1

    points = [ [ one(AA.QQ) for _ in 1:nvs ] for _ in 1:npoints ]
    for i in 1:npoints
        points[i][varidx] = AA.QQ(i)
    end

    xs = []
    ys = []
    for i in 1:npoints
        push!(xs, points[i][varidx])
        push!(ys, AA.evaluate(numerator(f), points[i]) // AA.evaluate(denominator(f), points[i]))
    end

    uniring, x = AA.PolynomialRing(AA.QQ, "x")

    inter = interpolate_rational_function(uniring, xs, ys)

    degf = AA.degree(numerator(f), vs[varidx]), AA.degree(denominator(f), vs[varidx])
    degi = AA.degree(numerator(inter)), AA.degree(denominator(inter))

    @test degi == degf
end

@testset "Polynomial Reconstruction tests" begin
    R, a = Nemo.PolynomialRing(Nemo.QQ, "a")

    ms = [
        (a - 2)*(a - 4)*(a + 5//3)*(a - 10//21)*(a - 1//8),
        (a^2 + a + 1)*(a + 10)^3*(a - 1)*(a - 2)
    ]
    fs = [
        a^2 + a + 1,
        (a - 3)^2,
        a,
        R(1)
    ]

    for f in fs
        for m in ms
            q = polynomial_reconstruction(f, m)
            n, d = unpackfrac(q)

            @test mod(n, m) == mod(d*f, m)
            @test AA.degree(n) + AA.degree(d) <= AA.degree(m)
        end
    end

    #=
    # should throw as long as the factor in f vanishes
    f = (a - 2)*(a - 1)*(a - 4)
    @test_throws DomainError polynomial_reconstruction(f, ms[1])
    @test_throws DomainError polynomial_reconstruction(f, ms[2])
    =#

end

#=
@testset "Rational Interpolation tests" begin

    R, (x1, x2, x3) = AA.PolynomialRing(AA.QQ, ["x1", "x2", "x3"])

    fs = [
      (x1 + x2) // x2,
      x1^5 // (x2 + x1 + x3),
      (x2 + x3) // (x1^3 + x1*x2 + x1*x3),
      (x1^4*x2^3 + x1^4*x3 + x2) // x1,
      (x1^6 + x2) // (x1^9 + x3),
      (x1^36 + x3^12 + x1^23 - x2^2 - x2*x3^2 + x1^2*x2^3+(2*x3^4)) // (x1^2 + x2 + x1*x2 + x3)
    ]

    for f in fs
        test_rational_function_interpolation(f)
    end

end
=#

###############################################################################

#=
@testset "Interpolation + shift tests" begin

    R, (x1, x2, x3) = AA.PolynomialRing(AA.QQ, ["x1", "x2", "x3"])

    f = (x1^2*x2 + x3) // (x2 + 3)

    npoints = 100
    ndeg = 3
    nvs = length(AA.gens(R))

    points = [ [ BigInt(j)^i for i in [ (ndeg+1)^k for k in 1:nvs ] ] for j in 1:npoints ]

    xs = []
    ys = []

    shifts, inverses = random_linear_shift(AA.QQ, length(points[1]))
    fp = [ fj(pj) for (fj, pj) in zip(shifts, AA.gens(R)) ]
    g = AA.evaluate(numerator(f), fp) // AA.evaluate(denominator(f), fp)

    for i in 1:npoints
        p = points[i]

        push!(xs, p[1])
        push!(ys,
            AA.evaluate(numerator(g), p) // AA.evaluate(denominator(g), p)
        )
    end

    uniring, x = AA.PolynomialRing(AA.QQ, "x")

    inter = interpolate_multivariate_rational_function(R, uniring, xs, ys)

    lol = [R(0), R(0)]
    for (i, D) in [ [1, numerator(inter)], [2, denominator(inter)] ]
        polybuilder = AA.MPolyBuildCtx(R)
        for (e, c) in zip(0:AA.degree(D), AA.coefficients(D))
            if c == 0
                continue
            end
            AA.push_term!(polybuilder, c, decompose_by_degrees(e, ndeg+1, nvs))
        end
        lol[i] = AA.finish(polybuilder)
    end

    F = lol[1] // lol[2]

    shiftedvariables = [ fj(pj) for (fj, pj) in zip(inverses, AA.gens(R)) ]

    F = AA.evaluate(numerator(F), shiftedvariables) // AA.evaluate(denominator(F), shiftedvariables)

    F = AA.leading_coefficient(numerator(F)) *inv(AA.leading_coefficient(denominator(F))) * (numerator(F) * inv(AA.leading_coefficient(numerator(F)))) // (denominator(F) * inv(AA.leading_coefficient(denominator(F))))

    @test isone(denominator(F // f)) && isone(numerator(F // f))

end
=#

#=
@testset "Degree prediction tests" begin

    fs0 = [
        (x) -> 2,
        (x) -> x,
        (x) -> 1 // x,
        (x) -> (x^2 + 11) // (x + 111)
    ]
    answers0 = [
        (0, 0),
        (1, 0),
        (0, 1),
        (2, 1)
    ]

    fs = [
        (x) -> (x^3 + 1),
        (x) -> (x^16 - 23x^11 + 123),
        (x) -> (x^2 // (x + 1)^4),
        (x) -> ((x - 1) // (x^18 + x^17 - 123)),
        (x) -> ((5x^113 + 2x^55 + 3x^33 + 1) // (x^11 + x + 22))
    ]

    answers = [
        (3, 0),
        (16, 0),
        (2, 4),
        (1, 18),
        (113, 11)
    ]

    FF = AA.GF(2^31-1)

    for (i, f) in enumerate(fs0)
        d = predict_degrees(f, AA.QQ)
        @test answers0[i] == d
    end

    for (i, f) in enumerate(fs)
        d = predict_degrees(f, FF)
        @test answers[i] == d
    end

end
=#

#################################################################

#=
@testset "Kronecker points generation tests" begin
    ground = Sing.GF(2^31 - 1)
    R, (a, b, c) = AA.PolynomialRing(ground, ["a", "b", "c"])

    maxexp = 10
    nvariables = 3

    polys = [
        a,
        a + b + c,
        a^2 + b^2 + c^3,
        a*c + b*c + c^3,
        a*b*c + a^4 + b,
        a^1*b^2*c^3 - 2a - 3b - 4c,
        b^10
    ]

    xs, points = generate_kronecker_points(ground, maxexp, nvariables)

    for poly in polys
        for (point, x) in zip(points, xs)
            symbolic_sub = [a, a^(maxexp+1), a^((maxexp+1)^2)]
            ans = AA.evaluate(AA.evaluate(poly, symbolic_sub), [x, ground(0), ground(0)])
            @test ans == AA.evaluate(poly, point)
        end
    end

end
=#

#=
@testset "Good Kronecker tests" begin
    ns = [1, 5, 37]
    prods = [[1, 2], [1, 2, 6], [1, 5, 25]]

    anses = [
        [[1, 0],  [1, 0, 0], [1, 0, 0]],
        [[1, 2],  [1, 2, 0], [0, 1, 0]],
        [[1, 18], [1, 0, 6], [2, 2, 1]]
    ]

    for (n, ans) in zip(ns, anses)
        hmm = map(p -> decompose_by_good_degrees(n, p, length(p)), prods)

        @test hmm == ans
    end
end
=#

#=
@testset "Kronecker backward substitution tests" begin
    FF = Sing.QQ
    R, (x1, x2, x3) = AA.PolynomialRing(FF, ["x1", "x2", "x3"])
    uni, x = AA.PolynomialRing(FF, "x")

    truefs = [
        x1,
        x2,
        x3,
        x1^2 + x1*x2^2 + x1*x2*x3 + x3 + 1,
        x1 + x2^5 + x3,
    ]
    fs = [
        x,
        x^2,
        x^4,
        x^2 + x*x^8 + x*x^4*x^16 + x^16 + 1,
        x + x^30 + x^36
    ]
    maxexps = [
        1,
        1,
        1,
        3,
        5
    ]

    for (truef, f, maxexp) in zip(truefs, fs, maxexps)
        @test backward_kronecker(f, R, maxexp) == truef
    end
end
=#

@testset "Polynomial interpolation tests" begin

    # over Singular finite field
    FF = Sing.GF(2^31 - 1)

    uniring, x = AA.PolynomialRing(FF, "x")
    test_cases = [
        uniring(1),
        x,
        x^2 + x + 1,
        x^4 - x,
        x^18 + x^17 + 22x^10 + 11
    ]

    for case in test_cases
        numpoints = AA.degree(case) + 1
        xs = [rand(FF) for _ in 1:numpoints]
        ys = [AA.evaluate(case, t) for t in xs]
        ans = interpolate_polynomial(uniring, xs, ys)
        @test ans == case
    end

    # over Singular QQ
    FF = Sing.QQ

    uniring, x = AA.PolynomialRing(FF, "x")
    test_cases = [
        uniring(1),
        x,
        x^2 + x + 1,
        x^4 - x,
        x^18 + x^17 + 22x^10 + 11
    ]

    for case in test_cases
        numpoints = AA.degree(case) + 5
        xs = [rand(FF) for _ in 1:numpoints]
        ys = [AA.evaluate(case, t) for t in xs]
        ans = interpolate_polynomial(uniring, xs, ys)
        @test ans == case
    end
end

@testset "Interpolation tests" begin

    FF = AA.GF(2^31 - 1)
    uniring, x = AA.PolynomialRing(FF, "x")
    test_cases = [
        1 // x,
        (x + 1) // (x + 2),
        (x^3 - 7) // (3 * x^2 + 1),
        (x^18 - 16) // (x + 1),
        (x + 1)^100 // (x + 2)^100,
    ]

    for case in test_cases
        n, d = unpackfrac(case)
        numpoints = 2 * (AA.degree(n) + AA.degree(d))
        xs = [rand(FF) for _ in 1:numpoints]
        ys = [AA.evaluate(n, xval) // AA.evaluate(d, xval) for xval in xs]
        ans = interpolate_rational_function(uniring, xs, ys)
        @test ans == case
    end

    xs = [FF(2), FF(4), FF(8), FF(16)]
    ys = [FF(1), FF(1), FF(1), FF(1)]
    ans = interpolate_rational_function(uniring, xs, ys)
    @test ans == FF(1)

    xs = [FF(2), FF(3)]
    ys = [FF(1), FF(1)]
    ans = interpolate_rational_function(uniring, xs, ys)
    @test ans == FF(1)

end

#=
# for future

using Logging

logger = Logging.SimpleLogger(stderr, Logging.Debug)

Logging.global_logger(logger)
@debug "lol"

uniring, x = AA.PolynomialRing(AA.QQ, "x")
f = (x) -> (x + 1) // (x + 2),

ans = interpolate_rational_function(uniring, xs, ys)
println(ans)

=#


text:
module Fleck

include("prefixsearch.jl")
include("sample/neverdist.jl")
include("sample/interface.jl")
include("sample/firstreaction.jl")
include("sample/naive.jl")
include("sample/fixeddirect.jl")
include("sample/direct.jl")
include("sample/anderson.jl")
include("vas.jl")

export Never

end

ext: jl


text:
module MaxHelpingHandReskinnableCrushBlock

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/ReskinnableCrushBlock" ReskinnableCrushBlock(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight,
	axes::String="both", chillout::Bool=false, spriteDirectory::String="objects/crushblock", fillColor::String="62222b", crushParticleColor1::String="ff66e2", crushParticleColor2::String="68fcff",
    activateParticleColor1::String="5fcde4", activateParticleColor2::String="ffffff")

const placements = Ahorn.PlacementDict(
    "Kevin (Both, Reskinnable) (max480's Helping Hand)" => Ahorn.EntityPlacement(
        ReskinnableCrushBlock,
        "rectangle"
    ),
    "Kevin (Vertical, Reskinnable) (max480's Helping Hand)" => Ahorn.EntityPlacement(
        ReskinnableCrushBlock,
        "rectangle",
        Dict{String, Any}(
            "axes" => "vertical"
        )
    ),
    "Kevin (Horizontal, Reskinnable) (max480's Helping Hand)" => Ahorn.EntityPlacement(
        ReskinnableCrushBlock,
        "rectangle",
        Dict{String, Any}(
            "axes" => "horizontal"
        )
    ),
)

Ahorn.editingOptions(entity::ReskinnableCrushBlock) = Dict{String, Any}(
    "axes" => Maple.kevin_axes
)

Ahorn.minimumSize(entity::ReskinnableCrushBlock) = 24, 24
Ahorn.resizable(entity::ReskinnableCrushBlock) = true, true

Ahorn.selection(entity::ReskinnableCrushBlock) = Ahorn.getEntityRectangle(entity)

# Todo - Use randomness to decide on Kevin border
function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::ReskinnableCrushBlock, room::Maple.Room)
    kevinColor = Ahorn.argb32ToRGBATuple(parse(Int, get(entity.data, "fillColor", "62222b"), base=16))[1:3] ./ 255
    spriteDirectory = get(entity.data, "spriteDirectory", "objects/crushblock")

    frameImage = Dict{String, String}(
        "none" => "$(spriteDirectory)/block00",
        "horizontal" => "$(spriteDirectory)/block01",
        "vertical" => "$(spriteDirectory)/block02",
        "both" => "$(spriteDirectory)/block03"
    )

    smallFace = "$(spriteDirectory)/idle_face"
    giantFace = "$(spriteDirectory)/giant_block00"

    axes = lowercase(get(entity.data, "axes", "both"))
    chillout = get(entity.data, "chillout", false)

    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    giant = height >= 48 && width >= 48 && chillout
    face = giant ? giantFace : smallFace
    frame = frameImage[lowercase(axes)]
    faceSprite = Ahorn.getSprite(face, "Gameplay")

    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    Ahorn.drawRectangle(ctx, 2, 2, width - 4, height - 4, kevinColor)
    Ahorn.drawImage(ctx, faceSprite, div(width - faceSprite.width, 2), div(height - faceSprite.height, 2))

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, frame, (i - 1) * 8, 0, 8, 0, 8, 8)
        Ahorn.drawImage(ctx, frame, (i - 1) * 8, height - 8, 8, 24, 8, 8)
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, 0, (i - 1) * 8, 0, 8, 8, 8)
        Ahorn.drawImage(ctx, frame, width - 8, (i - 1) * 8, 24, 8, 8, 8)
    end

    Ahorn.drawImage(ctx, frame, 0, 0, 0, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, width - 8, 0, 24, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, 0, height - 8, 0, 24, 8, 8)
    Ahorn.drawImage(ctx, frame, width - 8, height - 8, 24, 24, 8, 8)
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
module Day10

using LinearAlgebra

function readdata(s)
    sidelen = s |> length |> sqrt |> Int
    Set((j - 1, i) for i in 0:sidelen - 1, j in 1:sidelen if s[i * sidelen + j] == '#')
end

asteroids = read("day10.e", String) |> x -> replace(x, "\n" => "") |> readdata

angle(pos, other) = atan(other[2] - pos[2], other[1] - pos[1])

function numseen(asteroids, pos)
    slopes = Set()
    count(asteroids) do other
        (other == pos || angle(pos, other) in slopes) && return false
        push!(slopes, angle(pos, other))
        return true
    end
end

solvea(asteroids) = argmax(Dict(x => numseen(asteroids, x) for x in asteroids))

function solveb(asteroids) #incomplete
    start = solvea(asteroids)
    asts = [(x, angle(start, x)) for x in asteroids if x != start]
    sort!(asts, by = x -> x[2])
    ang = -pi / 2
    nextval = findfirst(x -> x[2] >= ang, asts) # straight up
    for _ in 1:199
        idxs = findall(x -> x[2] == asts[nextval][2], asts)
        idx = argmin([x -> (x, norm(start .- asts[x][1])) for x in idxs])
        idx == nothing && (idx = 1)
        ang = asts[idx[1]][2]
        @info asts[idx]
        splice!(asts, idx)
    end

    asts[idx]
end

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
AbstractString
println("all string types are subtypes of the abstract type AbstractString")
println(typeof('x'))
str = """asdfkasdf"quote" """
str = "hello, world\n"
for i in 1:length(str)
	println(str[i])
end
println(str[end-1])
println(str[2:4])

println("hello" * "world")  # string concatenation is noncommutative
println("interpolation $str")

println(occursin("a", "asdf"))
println(occursin("x", "sdf"))  # 'a' in 'sdf'
println(length(str))
println(repeat(str, 10))
println(join([str, str]))

println("Julia's regex is perl-compatible (PCRE)")
occursin(r"^\s*(?:#|$)", "not a comment")  # false
match(r"^\s*(?:#|$)", "# a comment")  # RegexMatch("#")
m = match(r"#", "#")
if m === nothing
	println("nothing")
else
	println(m)
end

println(replace("Programming Python", "Python" => "Julia"))
println("perlre(1)")

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
using JSON3

module API
using StructTypes
include("graphqlgen_types.jl")
include("graphqlgen_functions.jl")
end;

NODE_QUERY = """
id
"""
f = API.node(NODE_QUERY)

# Returns a GraphQL query and the variables required to execute it.
x = f("10")

# You can serialize the variables to JSON (as an example) and send use it
# in an HTTP request.
j = JSON3.write(x.variables)
# output: "{\"id\":\"10\"}"


text:
module StokesTaylorHoodTests

using Test
using Gridap
import Gridap: ∇

using LinearAlgebra: tr, ⋅

# Using automatic differentiation
u(x) = VectorValue( x[1]^2 + 2*x[2]^2, -x[1]^2 )
p(x) = x[1] + 3*x[2]
f(x) = -Δ(u)(x) + ∇(p)(x)
g(x) = (∇⋅u)(x)
∇u(x) = ∇(u)(x)

#u(x) = VectorValue( x[1]^2 + 2*x[2]^2, -x[1]^2 )
#∇u(x) = TensorValue( 2*x[1], 4*x[2], -2*x[1], zero(x[1]) )
#Δu(x) = VectorValue( 6, -2 )
#
#p(x) = x[1] + 3*x[2]
#∇p(x) = VectorValue(1,3)
#
#f(x) = -Δu(x) + ∇p(x)
#g(x) = tr(∇u(x))
#
#∇(::typeof(u)) = ∇u
#∇(::typeof(p)) = ∇p

domain = (0,2,0,2)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

labels = get_face_labeling(model)
add_tag_from_tags!(labels,"dirichlet",[1,2,5])
add_tag_from_tags!(labels,"neumann",[6,7,8])

order = 2

ref_style = [:reference,:physical]

for ref_st in ref_style
  V = TestFESpace(
  model=model,
  order=order,
  reffe=:Lagrangian,
  labels=labels,
  valuetype=VectorValue{2,Float64},
  dirichlet_tags="dirichlet",
  dof_space=ref_st,
  conformity=:H1)

  Q = TestFESpace(
  model=model,
  order=order-1,
  reffe=:Lagrangian,
  valuetype=Float64,
  dof_space=ref_st,
  conformity=:H1)

  U = TrialFESpace(V,u)
  P = TrialFESpace(Q)

  Y = MultiFieldFESpace([V,Q])
  X = MultiFieldFESpace([U,P])

  trian = get_triangulation(model)
  degree = order
  quad = CellQuadrature(trian,degree)

  btrian = BoundaryTriangulation(model,labels,"neumann")
  bdegree = order
  bquad = CellQuadrature(btrian,bdegree)
  n = get_normal_vector(btrian)

  function a(x,y)
    u,p = x
    v,q = y
    ∇(v)⊙∇(u) - (∇⋅v)*p + q*(∇⋅u)
  end

  function l(y)
    v,q = y
    v⋅f + q*g
  end

  function l_Γb(y)
    v,q = y
    v⋅(n⋅∇u) - (n⋅v)*p
  end

  t_Ω = AffineFETerm(a,l,trian,quad)
  t_Γb = FESource(l_Γb,btrian,bquad)

  op = AffineFEOperator(X,Y,t_Ω,t_Γb)

  uh, ph = solve(op)

  eu = u - uh
  ep = p - ph

  l2(v) = v⋅v
  h1(v) = v⋅v + ∇(v)⊙∇(v)

  eu_l2 = sqrt(sum(integrate(l2(eu),trian,quad)))
  eu_h1 = sqrt(sum(integrate(h1(eu),trian,quad)))
  ep_l2 = sqrt(sum(integrate(l2(ep),trian,quad)))

  tol = 1.0e-9
  @test eu_l2 < tol
  @test eu_h1 < tol
  @test ep_l2 < tol
end

end # module


text:
#
# Multi connection server implementation. Works through TCP
# protocol for connection with clients. Fully asynchronous and
# one threaded. It implements RPC-similar logic, where clients
# may run functions with parameters on server and obtain results.
#
# It uses green threads (or coroutines) inside. See this link
# http://julia.readthedocs.org/en/latest/manual/control-flow/#man-tasks
# for details regarding coroutines. Because of green threads,
# it uses "magic" inside: there are two places in this module,
# where this magic happens. First place is some main code or
# code, which creates the server and run it. It may be any
# other module or code. Second place is internal server's event
# loop machine, which is used for accepting clients connections
# and obtaining data (see @async macro in code). So, if you create
# a server with create() method, it will return Server's data
# object imediately. Calling run() method it also returns control
# back at that moment, but in reality server starts to listen
# clients in some loop (in our example this loop is out the run()
# code - in a parent code). After running, server creates one task
# (green thread) per one connection. To run server out you have to
# call yield() function in your (parent) code all the time.
# See example below for details:
#
#     #
#     # This callback function will be called as many times
#     # as server will obtain client's request (command)
#     #
#     function onCommand(sock::Base.TCPSocket, cmd::Connection.Command, ans::Connection.Answer)
#       # This is how we create the answer for client.
#       # ans.data has type Any.
#       ans.data = "answer"
#     end
#     #
#     # Creation of server doesn't run it
#     #
#     connection = Server.create(ip"127.0.0.1", 2016)
#     #
#     # Before running we have to bind command event listeners
#     #
#     Event.on(connection.observer, Connection.EVENT_BEFORE_RESPONSE, onCommand)
#     #
#     # This is how our server run itself
#     #
#     Server.run(connection)
#
#     ...
#
#     #
#     # This is client's (parent) code, which was mentioned before
#     #
#     while <condition>
#       # ...your stuff is here...
#
#       #
#       # This call is used for switching between Tasks (green threads),
#       # checks new connections, obtaining commands from clients,...
#       yield()
#     end
#     ...
#     #
#     # This is how we stop our server
#     #
#     Server.stop(connection)
#
# You have two possibilities how to run main code: you may have
# a loop (like in example above) or your code may wait for
# connection in background (without blocking loop). For example,
# you may run this example in a REPL without infinite loop at the
# end. In this case, you don't need to call yield() manually.
#
# Events:
#     command{Connection.Command, Connection.Answer} Fires if new
#            command was obtained from client. Contains command
#            itself and special answer object, where you may set your
#            custom data.
#
# @author DeadbraiN
# TODO: add EVENT_AFTER_REQUEST logic description
# TODO: add description of fast part of this module
module Server
  import Config.@if_debug
  import Event
  import Helper
  using Connection

  export create
  export run
  export stop
  export isOk
  export ServerConnection
  #
  # Describes a server. It contains clients sockets, tasks, server object
  # and it's observer.
  #
  type ServerConnection
    tasks   ::Array{Task, 1}
    socks   ::Array{Base.TCPSocket, 1}
    server  ::Base.TCPServer
    observer::Event.Observer
    host    ::Base.IPAddr
    port    ::Int
    fast    ::Bool
  end
  #
  # Creates a server. Returns special server's data object, which identifies
  # this server and takes an ability to listen it events. It also contains
  # server's related tasks. See Server.ServerConnection type for
  # details. It doesn't run the server (use run() method for this), but
  # it start to listen specified host and port using Base.listen() method.
  # Setting port to zero you may create "empty" connection. In fact, in
  # this case, server will not be created and will not work.
  # @param host Host we are listening to
  # @param port Port we are listening to
  # @param fast Turns on serialization mode (slow speed mode).
  # @return {Server.ServerConnection} Server's related data object
  #
  function create(host::Base.IPAddr, port::Int, fast::Bool = false)
    local tasks::Array{Task, 1} = Task[]
    local socks::Array{Base.TCPSocket, 1} = Base.TCPSocket[]
    local obs::Event.Observer = Event.create()
    local con::ServerConnection

    if port > 0
      try
        local server::Base.TCPServer = listen(host, port)
        con = ServerConnection(tasks, socks, server, obs, host, port, fast)
        Helper.info(string(fast ? "Fast" : "Slow", " Server created on ", con.host, ":", con.port))
        return con
      catch e
        Helper.warn("Server.create(): $e")
        @if_debug showerror(STDOUT, e, catch_backtrace())
      end
    end

    ServerConnection(tasks, socks, Base.TCPServer(), obs, host, port, fast)
  end
  #
  # Runs the server. Starts listening clients connections
  # and starts answering on requests. This method implements
  # main asynchronous client-server communication logic. Here
  # all green threads are used. See this link for details:
  # http://julia.readthedocs.org/en/latest/manual/control-flow/#man-tasks
  # Don't remember to call yield() in your main code (or loop).
  # It also calls _update() for removing closed connections
  # ans failed tasks. It means that
  # serialize()/deserialize() functions will be used for sending and
  # receiving data. Otherwise native read()/write() functions will be
  # used.
  # @param con Server connection object returned by Server.create()
  # @return {Bool} Run status
  #
  function run(con::Server.ServerConnection)
    if !isOk(con)
      Helper.warn("Server.run(): Server wasn\'t created correctly. Try to change Server.create() arguments.")
      return false
    end

    Helper.info(string(con.fast ? "Fast" : "Slow", " server has run on ", con.host, ":", con.port))
    @async begin
      while true
        try
          #
          # This line handles new connections
          #
          local sock::Base.TCPSocket = accept(con.server)
          push!(con.socks, sock)
        catch e
          @if_debug showerror(STDOUT, e, catch_backtrace())
          #
          # Possibly Server.stop() was called.
          #
          if !isOk(con) _close(con); break end
          Helper.warn("Server.run(): $e")
        end
        push!(con.tasks, @async while Helper.isopen(sock)
          yield()
          con.fast ? fastAnswer(sock, con.observer, _onAnswerException) : answer(sock, con.observer, _onAnswerException)
          _update(con)
        end)
      end
    end
    #
    # This yield() prevents server from error:
    # ArgumentError("server not connected, make sure \"listen\" has been called")
    #
    yield()
    true
  end
  #
  # Makes request to client. This method is not blocking. It returns
  # just after the call. Answer will be obtained in run() method
  # async loop.
  # @param sock Client socket returned by accept() function
  # @param fn Remote function id
  # @param args Custom fn arguments
  # @return true - request was sent, false wasn't
  #
  function request(sock::Base.TCPSocket, fn::Int, args...)
    Connection.request(sock, fn, args...)
  end
  #
  # Fast version of request function. It uses native read()/write()
  # function for data transfer. This function is non blocking one.
  # TODO: describe dataIndex and data
  # @param sock Client socket returned by accept() function
  # @param dataIndex Index in api array
  # @param args Data to transfer in request
  #
  function request(sock::Base.TCPSocket, dataIndex::UInt8, args...)
    fastRequest(sock, dataIndex, args...)
  end
  #
  # Returns server's state. true means - created and run.
  # @param con Client connection state
  # @return {Bool}
  #
  function isOk(con::Server.ServerConnection)
    Helper.isopen(con.server)
  end
  #
  # Stops the server. Stops listening all connections and drops
  # existing if exist.
  # @param con Server object returned by create() method.
  #
  function stop(con::Server.ServerConnection)
    try
      _close(con)
      Helper.info(string("Server has stopped: ", con.host, ":", con.port))
    catch e
      Helper.warn("Server.stop(): $e")
      @if_debug showerror(STDOUT, e, catch_backtrace())
    finally
      yield()
    end
  end

  #
  # Closes all available clients connections and server's one
  # @param con Server object returned by create() method.
  #
  function _close(con::Server.ServerConnection)
    for i::Int = 1:length(con.socks)
      close(con.socks[i])
      yield()
    end
    close(con.server)
    yield()
  end
  #
  # This method should be called in main server's loop or outside
  # code for removing stopped tasks and sockets (connections)
  # between clients and this server. See an example at the beginning
  # of this module for details.
  # @param con Connection object returned by create() method
  #
  function _update(con::Server.ServerConnection)
    i::Int = 1

    while i <= length(con.socks)
      if Helper.isopen(con.socks[i])
       	i += 1
      else
        deleteat!(con.socks, i)
        deleteat!(con.tasks, i)
      end
    end
  end
  #
  # Is called when answer() or fastAnswer() catch network
  # related exception
  # @param sock Client's socket we are working with
  # @param e Exception object
  # @return true Always true
  #
  function _onAnswerException(sock::Base.TCPSocket, e::Exception)
    #
    # This error means that client has disconnected
    #
    if isa(e, EOFError)
      Helper.info("Client has disconnected")
    else
      Helper.warn(string("Server._onAnswerException(): ", e))
    end
    if Helper.isopen(sock) close(sock) end

    true
  end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
@testset "Manopt Cyclic Proximal Point" begin
using Dates
n = 100
N = Power(Circle(),(n,))
f = PowPoint(artificialS1Signal(n))
F = x -> costL2TV(N,f,0.5,x)
proxes = [ (λ,x) -> proxDistance(N,λ,f,x), (λ,x) -> proxTV(N,0.5*λ,x) ]
fR,rec = cyclicProximalPoint(N,F,proxes, f;
    λ = i -> π/(2*i),
    stoppingCriterion = stopWhenAll( stopAfter(Second(10)), stopAfterIteration(5000) ),
    debug = [DebugIterate()," ",DebugCost()," ",DebugProximalParameter(),"\n",10000],
    record = [RecordProximalParameter(), RecordIterate(f), RecordCost()]
    )
@test F(f) > F(fR)
#
o = CyclicProximalPointOptions(f, stopAfterIteration(1), i -> π/(2*i))
p = ProximalProblem(N,F,proxes,[1,2])
@test_throws ErrorException getProximalMap(p,1.,f,3)
end


text:
using DataFrames, PyCall, PyPlot, Query
import Mortchartgen, CSV, Requests, ZipFile, DataStructures

rffiles = Dict(
	       :chol => "IHME_GBD_2015_COVARIATES_1980_2015_MEAN_CHOLESTEROL_Y2016M10D07",
	       :sbp => "IHME_GBD_2015_COVARIATES_1980_2015_MEAN_SBP_Y2016M10D07",
	       :smok => "IHME_GBD_2015_COVARIATES_1980_2015_SMOKING_PREV_Y2016M10D07"
	       )

function download_gbd(filename)
	mkpath("data")
	preurl = "https://cloud.ihme.washington.edu/index.php/s/SotC3shZAcGVFch/download?path=%2FGBD%202015%20Covariates&files="
	filefullname = "$filename.zip"
	fileurl = "$preurl$filefullname"
	dlpath = joinpath("data", filefullname)
	Requests.save(Requests.get(fileurl), dlpath)
	r = ZipFile.Reader(dlpath)
	for f in r.files
		write(joinpath("data", f.name), read(f))
	end
	rm(dlpath)
end

for fname in values(rffiles)
	if !(isfile(joinpath("data", "$fname.CSV")))
		download_gbd(fname)
	end
end

PyDict(matplotlib["rcParams"])["axes.formatter.use_locale"] = true

if !(isdefined(:frames))
	frames = Mortchartgen.load_frames()
end

gbdlid2who = Dict(67=>3160, 78=>4050, 79=>4070, 90=>4220, 93=>4290, 95=>4310, 102=>2450)

rfs = Dict(:chol=>:mean, :sbp=>:mean, :smok=>:prev)

## RRS from:
#  - doi:10.1371/journal.pone.0065174 (chol, sbp)
#  - doi:10.1161/CIRCULATIONAHA.104.521708 (smok)
ages = DataStructures.OrderedDict(
	:a35_44 => Dict(:sage=>13, :eage=>14, :gbd_ageid=>12,
			:chol=>2.20, :sbp=>1.68^(1/10), :smok=>Dict(1=>5.51, 2=>2.26)),
	:a45_54 => Dict(:sage=>15, :eage=>16, :gbd_ageid=>14,
			:chol=>1.82, :sbp=>1.56^(1/10),:smok=>Dict(1=>3.04, 2=>3.78)),
	:a55_64 => Dict(:sage=>17, :eage=>18, :gbd_ageid=>16,
			:chol=>1.44, :sbp=>1.45^(1/10),
			:smok=>Dict(1=>mean([3.04;1.88]), 2=>mean([3.78;2.53]))),
	:a65_74 => Dict(:sage=>19, :eage=>20, :gbd_ageid=>18,
			:chol=>1.27, :sbp=>1.33^(1/10),
			:smok=>Dict(1=>mean([1.88;1.44]), 2=>mean([2.53;1.68]))),
	:a75_84 => Dict(:sage=>21, :eage=>22, :gbd_ageid=>20,
			:chol=>1.18, :sbp=>1.26^(1/10),
			:smok=>Dict(1=>mean([1.44;1.05]), 2=>mean([1.68;1.38])))
	)
ageints(age) = map(x->parse(x), split("$age"[2:end], '_'))
agelabel(age) = replace("$age", "_", "\u2013")[2:end]
rfframes = Dict(
	:chol => CSV.read(joinpath("data", "$(rffiles[:chol]).CSV")),
	:sbp => CSV.read(joinpath("data", "$(rffiles[:sbp]).CSV")),
	:smok => CSV.read(joinpath("data", "$(rffiles[:smok]).CSV"), datarow=14330),
	)
rfalsy(rfframe, gbd_ageid, gbd_locid, sex, syr, eyr) = @from i in rfframe begin
	@where i.age_group_id==gbd_ageid && i.location_id==gbd_locid &&
		i.sex_id==sex && i.year_id>=syr && i.year_id<=eyr
	@select i
	@collect DataFrame
end

function rrs(rf, age, gbd_locid, sex, syr, eyr, ind)
	if (rf==:chol || rf==:sbp)
		gbd_ageid=22
	else
		gbd_ageid=ages[age][:gbd_ageid]
	end
	frame = rfalsy(rfframes[rf], gbd_ageid, gbd_locid, sex, syr, eyr)
	agerr = ages[age][rf]
	if typeof(agerr)<:Dict
		agesexrr = agerr[sex]
	else
		agesexrr = agerr
	end
	if rfs[rf]==:mean
		rfdiff = frame[:val][ind].-frame[:val]
		agesexrr.^rfdiff
	elseif rfs[rf]==:prev
		rfr1 = frame[:val][ind].*agesexrr
		rfr = frame[:val].*agesexrr
		(rfr1.+(1.-frame[:val][ind])) ./ (rfr.+(1.-frame[:val]))
	end
end

transforms = Dict(
		  :lin=>Dict(:func=>identity, :label=>"dödstal"),
		  :log=>Dict(:func=>log, :label=>"log(dödstal)")
		  )

function agsrrs(age, gbd_locid, sex, syr, eyr)
	yrs = syr:eyr
	indyr = syr
	ind = 1
	who_country = gbdlid2who[gbd_locid]
	ihd = Mortchartgen.propplot_sexesyrs("ihd", "pop", [sex], who_country,
		ages[age][:sage], ages[age][:eage], syr:eyr, true, frames, "sv",
		joinpath(tempdir(), "mout.html"), false)
	rates = ihd[sex][:propframe][:value]
	rates_sm = ihd[sex][:propsm]
	rrs_chol = rrs(:chol, age, gbd_locid, sex, syr, eyr, ind)
	rrs_sbp = rrs(:sbp, age, gbd_locid, sex, syr, eyr, ind)
	rrs_smok = rrs(:smok, age, gbd_locid, sex, syr, eyr, ind)
	rs_c = rates_sm .* rrs_chol
	rs_cb = rs_c .* rrs_sbp
	rs_cbs = rs_cb .* rrs_smok
	ratediff = rates_sm[1]-rates_sm[end]
	rdiffs_cbs = rdiffs(rates_sm, ratediff, [rrs_chol rrs_sbp rrs_smok], ind)
	Dict(
	     :rrs_chol=>rrs_chol, :rrs_sbp=>rrs_sbp, :rrs_smok=>rrs_smok,
	     :rates=>rates, :rates_sm=>rates_sm, :rs_c=>rs_c, :rs_cb=>rs_cb, :rs_cbs=>rs_cbs,
	     :ratediff=>ratediff, :rdiffs_cbs=>rdiffs_cbs, :age=>age, :sex=>sex, :indyr=>indyr,
	     :yrs=>yrs, :who_country=>who_country
	     )
end

function plotsex(agsd, trans=:log)
	yrs = agsd[:yrs]
	indyr = agsd[:indyr]
	ind = 1
	clab = Mortchartgen.conf["countries"]["$(agsd[:who_country])"]["alias"]["sv"]
	slab = Mortchartgen.conf["sexes"]["$(agsd[:sex])"]["alias"]["sv"]
	plot(yrs, transforms[trans][:func].(agsd[:rates]), "*", label=L"r", color="blue")
	plot(yrs, transforms[trans][:func].(agsd[:rates_sm]), "-", label=L"r_{sm}", color="blue")
	plot(yrs, transforms[trans][:func].(agsd[:rs_c]), "-",
		label=latexstring("r_{sm} \\times RR_{TC($indyr)}"))
	plot(yrs, transforms[trans][:func].(agsd[:rs_cb]), "-",
		label=latexstring("r_{sm} \\times RR_{TC($indyr)} \\times RR_{SBP($indyr)}"))
	plot(yrs, transforms[trans][:func].(agsd[:rs_cbs]), "-",
		label=latexstring("r_{sm} \\times RR_{TC($indyr)} \\times RR_{SBP($indyr)} \\times RR_{SMOK($indyr)}"))
	plot(yrs, transforms[trans][:func].(agsd[:rdiffs_cbs][:rs_nobasech]), "-",
	     label=latexstring("r_{sm}($indyr) \\times 1/ \\prod RR_{i(x)}"))
	legend(loc=3)
	grid(1)
	xlabel("År")
	ylabel(transforms[trans][:label])
	title("Dödstal IHD vs kontrafaktiska riskfaktormönster, $slab $clab $(agelabel(agsd[:age]))")
end

function rdiffs(rates_sm, ratediff, rfs, ind)
	rs_nobasech = rates_sm[1] .* 1./prod(rfs, 2)
	rs_rfs = prod([rates_sm rfs], 2)
	ratediff_rf_st = rates_sm[1]-rs_nobasech[end]
	ratediff_rf_stp = ratediff_rf_st/ratediff
	ratediff_rf_end = rs_rfs[end]-rates_sm[end]
	ratediff_rf_endp = ratediff_rf_end/ratediff
	Dict(
	     :rs_nobasech=>rs_nobasech,
	     :rs_rfs=>rs_rfs, :ratediff_rf_st=>ratediff_rf_st, :ratediff_rf_stp=>ratediff_rf_stp,
	     :ratediff_rf_end=>ratediff_rf_end, :ratediff_rf_endp=>ratediff_rf_endp,
	 )
end

function plotrdiffps(ages, agsds, rdiffds)
	sages = map(a->ageints(a)[1], ages)
	agsd = agsds[1]
	clab = Mortchartgen.conf["countries"]["$(agsd[:who_country])"]["alias"]["sv"]
	slab = Mortchartgen.conf["sexes"]["$(agsd[:sex])"]["alias"]["sv"]
	stps = map(d->d[:ratediff_rf_stp], rdiffds)
	endps = map(d->d[:ratediff_rf_endp], rdiffds)
	plot(sages, stps, "-*", label="Andel startår")
	plot(sages, endps, "-*", label="Andel slutår")
	legend(loc=3)
	grid(1)
	xlabel("Ålder")
	ylabel("Förklarad andel")
	title("Förklarad andel av förändring IHD $(agsd[:yrs][1])\u2013$(agsd[:yrs][end]), $slab $clab")
end

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
struct PlotlinesModel{T₁ <: AbstractVector} <: AbstractModel
    data::T₁
end

mutable struct PlotlinesControl <: AbstractControl
    isenabled::Bool
end

struct PlotContext{T₁ <: AbstractControl,   T₂ <: AbstractModel,  T₃ <: AbstractDisplayProperties} <: AbstractPlotContext
    control::T₁
    model::T₂
    display_properties::T₃
end

function (context::PlotContext{<: PlotlinesControl,   <: PlotlinesModel,  <: PlotlinesDisplayProperties})()
        control = context.control
        model = context.model
        display_properties = context.display_properties
        id = get_id(display_properties)
        captioner = get_captioner(display_properties)
        caption = captioner(nothing)
        is_new_window(display_properties) ? CImGui.Begin(caption*id,C_NULL, CImGui.ImGuiWindowFlags_NoBringToFrontOnFocus) : nothing
        isenabled(control) ? control(model, display_properties) : nothing
        is_new_window(display_properties) ? CImGui.End() : nothing
end

function get_data(model::PlotlinesModel)
    model.data
end

function (control::PlotlinesControl)(model::PlotlinesModel, properties::PlotlinesDisplayProperties)
    id = get_id(properties)
    captioner = get_captioner(properties)
    caption = captioner(nothing)
    col = get_col(properties)
    rectangle = get_layout(properties)
    xaxis = get_xaxis(properties)
    yaxis = get_yaxis(properties)
    pos = CImGui.GetCursorScreenPos()
    totalwidth = get_width(rectangle)
    totalheight = get_height(rectangle)
    padding = get_padding(properties)
    data = get_data(model)
    CImGui.Dummy(ImVec2(0, padding[1]))
    CImGui.Indent(padding[2])
    width = totalwidth - padding[2]
    height = totalheight
    CImGui.PlotLines(id, data , length(data), 0 , caption, minimum(data), maximum(data), (width, height))
    isenabled(xaxis) ? draw_horizontal_ticks(data, width, height, get_ytick(properties)) : nothing
    isenabled(yaxis) ? draw_vertical_ticks(data, width, height, get_xtick(properties)) : nothing
    CImGui.Unindent(padding[2])
end

function draw_vertical_ticks(data, width, height, tick::Tickmark)
    if isenabled(tick)
        # Draw the x-axis
        draw_list = CImGui.GetWindowDrawList()
        black = Base.convert(ImU32, ImVec4(0,0, 0, 1))
        pos = CImGui.GetCursorScreenPos()
        CImGui.AddLine(draw_list, pos, ImVec2(pos.x + width, pos.y), black, get_thickness(tick));
        # Draw the concomitant tick marks.
        x = pos.x
        y = pos.y
        interpret = get_interpreter(tick)
        # for xₙ in range(x, stop = x + width; length = floor(Int, width / get_spacing(tick)) )
        #     # This line represents the tick mark.
        #     CImGui.AddLine(draw_list, ImVec2(xₙ, y), ImVec2(xₙ, y + get_length(tick)), black, get_thickness(tick)) #TODO add tick direction
        #     # Display value associated with that tickmark
        #     index = round(Int,stretch_linearly(xₙ, x,  x + width, 1, length(data)))
        #     # Convert the raw data to a text description after applying the
        #     # transformation associated with the value2text function.
        #     str = string(interpret(index))
        #     xoffset = div(length(str) * CImGui.GetFontSize(), 4)
        #     yoffset = get_length(tick) + 2
        #     CImGui.AddText(draw_list, ImVec2(xₙ - xoffset, y + yoffset), black, "$str",);
        # end
        gap = floor(Int, width / get_spacing(tick))
        increment = round(Int, stretch_linearly(gap , 1, width, 1, length(data)))
        for index in 1:increment:length(data)
            xₙ = Cfloat(stretch_linearly(index , 1, length(data), x,  x + width))
            # This line represents the tick mark.
            CImGui.AddLine(draw_list, ImVec2(xₙ, y), ImVec2(xₙ, y + get_length(tick)), black, get_thickness(tick)) #TODO add tick direction
            # Display value associated with that tickmark
            # Convert the raw data to a text description after applying the
            # transformation associated with the value2text function.
            str = string(interpret(index))
            xoffset = div(length(str) * CImGui.GetFontSize(), 4)
            yoffset = get_length(tick) + 2
            CImGui.AddText(draw_list, ImVec2(xₙ - xoffset, y + yoffset), black, "$str",);
        end
    end
end

function draw_horizontal_ticks(data, width, height, tick::Tickmark)

    if isenabled(tick)
        draw_list = CImGui.GetWindowDrawList()
        black = Base.convert(ImU32, ImVec4(0,0, 0, 1))
        col = Base.convert(ImU32, ImVec4(0.95,0.1, 0.1, 0.1))
        pos = CImGui.GetCursorScreenPos()
        # Draw the y-axis
        CImGui.AddLine(draw_list, pos, ImVec2(pos.x, pos.y - height), black, get_thickness(tick));
        minval = minimum(data)
        maxval = maximum(data)
        # Draw the concomitant tick marks.
        x = pos.x
        y = pos.y
        interpret = get_interpreter(tick)
        # TODO: Perhaps better to iterate over bonafide values and round pixel coordinates?
        span = range(y - height, stop = y; length = floor(Int, height / get_spacing(tick)))
        count = 1
        for yₙ in span
            # This line represents the tick mark.
            CImGui.AddLine(draw_list, ImVec2(x, yₙ), ImVec2(x - get_length(tick), yₙ), black, get_thickness(tick)); #TODO add tick direction
            # Display value associated with that tickmark
            val = round(stretch_linearly(yₙ, y,  y - height, minval, maxval); digits = 3)
            # Convert the raw data to a text description after applying the
            # transformation associated with the value2text function.
            str = string(interpret(val))
            xoffset = get_length(tick) + 3 + div(length(str) * CImGui.GetFontSize(), 2)
            yoffset = div(CImGui.GetFontSize(), 2)
            CImGui.AddText(draw_list, ImVec2(x - xoffset, yₙ - yoffset), black, "$str",);
            # Draw a faint filled rectangle as a visual guide for every second row
            if iseven(count)
                 CImGui.AddRectFilled(draw_list, ImVec2(x, yₙ), ImVec2(x + width, yₙ + step(span)), col);
            end
            count += 1
        end
    end
end


text:
# -*- coding: utf-8 -*-
# # Univariate Polynomial Regression
#
# This example demonstrates how to use manipulate a univariate regression model in an algebraically precise way. It builds on the example `pseudo_polynomial_regression.jl`.

# Let `m` be a statistical model like
# ```julia
# f(x) = βx
# X = rand(Normal(0,1), n)
# target[i] = f(X[i]) + rand(Normal(0,1), 1)
# min_{a,b} sum_i (f(X[i]) - target[i])^2
# ```
#
# So it is a linear least squares regression model. The parameter $\beta$ is the regression coefficients that are learned from data.
#
# Define the set of transformations $T$ as the free monoid generated by
#
# ```julia
# T_x = f(x)->x*f(x)
# T_1 = f(x)->f(x) + C
# ```
#
# Namely, the set of strings over the alphabet $\{T_x, T_1\}$
#
# The monoid generated by this set of transformations acts on `m` to create all polynomial univariate regressions.
#
# ```julia
# f(x) = sum(β_i*x^i)
# X = rand(Normal(0,1), n)
# target[i] = f(X[i]) + rand(Normal(0,1), 1)
# min_{β} sum_i (f_β(x[i]) - target[i])^2
# ```
#
# for all polynomials `f`. The ring of polynomials is well studied and we should be able to apply our knowledge of polynomial algebras to say something about the model transformations.
#
# We will  replicate the feature selection process used in standard polynomial regression models using our model augmentation tools.
#
# Here is an example implementation of this theory.

# ## Example Model

#Our working example of a multivariate nonlinear regression model
expr = quote
    module Regression
    using Random
    using LsqFit
    using LinearAlgebra

    function f(x, β)
        # This .+ node is added so that we have something to grab onto
        # in the metaprogramming. It is the ∀a .+(a) == a.
        return .+(β[1].* x.^0)
    end
    # Generate the data for this example
    function sample(g::Function, n)
        x = randn(Float64, n)
        target = g(x) .+ randn(Float64, n[1])./1600
        return x, target
    end

    # Compute the regression statistics
    function describe(fit)
        if !fit.converged
            error("Did not converge")
        end
        return (β = fit.param, r=norm(fit.resid,2), n=length(fit.resid))
    end
    #setup

    # Random.seed!(42)
    β = (1/2, 1/2)
    n = (1000)
    g(x) = β[1].*x.^3 + β[2].*x.^2
    X, target = sample(g, n)
    # x, y = X[:,1], X[:,2]
    # @show size(X), size(target)
    # loss(a) = sum((f.(a, x, y) .- target).^2)
    # @show loss.([-1,-1/2,-1/4, 0, 1/4,1/3,1/2,2/3, 1])

    a₀ = [1.5]
    try
        ŷ₀ = f(X, a₀)
        catch except
        @show except
        error("Could not execute f on the initial data")
    end

    #solving
    fit = curve_fit(f, X, target, a₀)#; autodiff=:forwarddiff)
    result = describe(fit)
    end
end
eval(expr.args[2])
Regression.fit.param

# ## Implementation Details
# The following code is the implementation details for representing the models as an `AbstractProblem` and representing the transformations as sequences of operations and applying the transformations onto the models.
#
# You can skip these if you are not interested in how this is implemented.

# +
using LinearAlgebra
import Base: show

using SemanticModels
using SemanticModels.ExprModels
import SemanticModels: model, AbstractModel
import SemanticModels.ExprModels: isexpr
import SemanticModels.ExprModels.Parsers: findfunc, findassign
import SemanticModels.ExprModels.Transformations: Pow
# -

ExprModels.Transformations.Pow

# +
"""    Lsq

A program that solves min_β || f(X,β) - y ||_2

Example:

`f(X, β) = β[1].*X.^p .+ β[2].*X.^q`

See also [`(t::Pow)(m::Lsq)`](@ref)
"""
struct Lsq <: AbstractModel
    expr
    f
    coefficient
    p₀
end

function show(io::IO, m::Lsq)
    write(io, "Lsq(\n  f=$(repr(m.f)),\n  coefficient=$(repr(m.coefficient)),\n  p₀=$(repr(m.p₀))\n)")
end

function model(::Type{Lsq}, expr::Expr)
    if expr.head == :block
        return model(Lsq, expr.args[2])
    end
    objective = :l2norm
    f = callsites(expr, :curve_fit)[end].args[2]
    coeff = callsites(expr, f)[1].args[end]
    p₀ = callsites(expr, :curve_fit)[end].args[end]
    return Lsq(expr, f, coeff, p₀)
end

"""    poly(m::Lsq)::Expr

find the part of the model that implements the polynomial model for regression.
"""
function poly(m::Lsq)
    func = findfunc(m.expr, m.f)[1]
    poly = func.args[2].args[end].args[1]
    return poly
end

"""    (t::Pow)(m::Lsq)

Example:

If `m` is a program implementing `f(X, β) = β[1]*X^p + β[2]*X^q`

a) and `t = Pow(2)` then `t(m)` is the model implementing
`f(X, β) = β[1]*X^p+2 + β[2]*X^q+q`.

"""
function (t::Pow)(m::Lsq)
    p = poly(m)
    for i in 2:length(p.args)
        slot = p.args[i]
        pow = callsites(slot, :(.^))
        pow[end].args[3] += t.inc
    end
    return m
end

struct AddConst <: ExprModels.Transformations.Transformation end

"""    (c::AddConst)(m::Lsq)

Example:

If `m` is a program implementing `f(X, β) = β[1]*X^p + β[2]*X^q`

a) and `c = AddConst()` then `c(m)` is the model implementing
`f(X, β) = β[1]*X^p + β[2]*X^0`.

"""
function (c::AddConst)(m::Lsq)
    p = poly(m)
    ix = map(t->t.args[2].args[2], p.args[2:end])
    i = maximum(ix)+1
    @show p
    push!(p.args, :(β[$i].*x.^0))
    assigns = findassign(m.expr, m.p₀)
    @show assigns
    b = assigns[end].args[2].args
    push!(b, 1)
    return m
end
# -

m = model(Lsq, deepcopy(expr))

# ## Model Augmentation with Group Actions
#
# Now we have all the machinery in place to build novel models from old models.

# Let's build an instance of the model object from the code snippet expr
m = model(Lsq, deepcopy(expr))
@show m
poly(m)

# Some *generator elements* will come in handy for building elements of the transformation group.
# $T_x,T_1$ are *generators* for our group of transformations $T = \langle T_x, T_1 \rangle$. Application of $T_1$ adds a constant to our polynomial while application of $T_x$ increments all the powers of the terms by 1. Any polynomial can be generated by these two operations.

@show Tₓ = Pow(1)
@show T₁ = AddConst()
Tₓ, T₁

# TODO: make these tests

m = model(Lsq, deepcopy(expr))
@show poly(Tₓ(m))
@show poly(T₁(m))
@show poly(Tₓ(m))
@show poly(Tₓ(m))
@show poly(T₁(m))

# Theorem: Any polynomial regression model can be produced by repeated application of $T_x,T_1$.
#
# Proof: Let $p(x) = \sum_i^d \beta_i x^i$ Starting from $q_0(x)= \beta[1] = \beta[1] x^0$,
#
# ```
# for i in 1:d
#    if \beta_{d-i} != 0 then
#       apply T_1
#    end
#    apply T_x
# end
# ```
#
# This algorithm is Horner's rule for evaluating $p(x)$ on the symbolic coeffients $\beta_i$

# ## Exploring the Orbits
#
# We use the Levenberg Marquardt algorithm for general least squares problems in data analysis. The goal is to fix a model class and find the best coefficients for that class. Our algebraic representation allows us to have a similar treatment of model class (in this case the exponents $i,j$ in our formula). We can sample from the orbits of the group when applied to the model and solve for the best coefficients in order to find the best model class.
#
# In this case we have reduced the code search space to a model search space to an algebraic formulation.

m = model(Lsq, deepcopy(expr))
results = []
for i in 1:5
    p = poly(m)
    M = eval(m.expr)
    push!(results, (i=i, M.result..., p=deepcopy(p)))
    T₁(Tₓ(m))
end

# It turns out that this method recovers the *true* model order $2,3$

for r in results
    # \t$(map(x->string(x.args[3]), r.p.args[2:end]))
    β̂ = map(βᵢ -> round(βᵢ, digits=4), r.β)
    s = join(map(x->string(x.args[3].args[3]), r.p.args[2:end]),", ")
    dof = length(r.p.args)-1
    println("DoF: $(dof)\tPoly: x^{$s}")
    println("Residual: $(round(r.r, digits=4))\tβ̂: $(β̂))")
    ℓ1 = norm(β̂, 1)
    println("ℓ₁ norm of β̂: $ℓ1")
    println("--------------------------------------------------------\n\n")
end

# ## Appied Algebra
#
# We know that the data generating function is an odd polynomial and that any polynomial in $\langle T_x\dot T_x, T_1\rangle$ will be even. The following metamodel shows that no even polynomial fits our data well.

# +
m = model(Lsq, deepcopy(expr))
results = []
for i in 1:5
    p = poly(m)
    M = eval(m.expr)
    push!(results, (i=i, M.result..., p=deepcopy(p)))
    T₁(Tₓ(Tₓ(m)))
end

for r in results
    # \t$(map(x->string(x.args[3]), r.p.args[2:end]))
    β̂ = map(βᵢ -> round(βᵢ, digits=4), r.β)
    s = join(map(x->string(x.args[3].args[3]), r.p.args[2:end]),", ")
    dof = length(r.p.args)-1
    println("DoF: $(dof)\tPoly: x^{$s}")
    println("Residual: $(round(r.r, digits=4))\tβ̂: $(β̂))")
    ℓ1 = norm(β̂, 1)
    println("ℓ₁ norm of β̂: $ℓ1")
    println("--------------------------------------------------------\n\n")
end
# -

# ## Conclusions
#
# We have seen how abstract algebra can be applied to the category of models to build a systematic treatment of model augmentation. This proves that the model transformations can be arranged into a simple algebraic structure that can act on a model to build new models. The structure of the transformations are easier to analyze than the changes to the models themselves.


text:
# Ruiz Omar, auditing the course
function q1()
@eval using Random
@eval using Distributions
@eval using LinearAlgebra
# a) We set the seed 1234
Random.seed!(1234)
# i. A 10x7 matrix  with random numbers U[-5,10]
A=rand(Uniform(-5,10),10,7)
# ii. B 10x7 matrix  with random numbers N(-2,15)
B=rand(Normal(-2.0,15.0),10,7)
# iii. C 5x7 matrix with first five rows and columns of A and first five rows and two last columns of B
C=[A[1:5,1:5] B[1:5,6:7]]
# iv. D a matrix with values of A but if they are positive, convert them to zeros.
D=A
D[D.>0].=0
D
# b) List the number of elements of A
prod(size(A))
# c) List unique elements of D
prod(size(unique(D)))
# d) Vectorize B and call it E with reshape
E=reshape(B,70,1)
# Alternative solution, not easier than wih reshape
E=[B[1:10,1]; B[1:10,2]; B[1:10,3] ; B[1:10,4]; B[1:10,5]; B[1:10,6]; B[1:10,7]]
# e) F, an array with A and B in the first and second columnsof the third dim, respectively
F=reshape([A B], (10,7,2))
# f) Reshape F so it is now a 2x7x10 array
F=permutedims(F, [3,2,1])
# g) G computes the Kronecker product of B and C
G=kron(B,C)
# When trying the same with C and F a problem occurs due to the fact that F is not a matrix, but an
# array containing 10 different so the Kronecker product is not computed.
# h) Save A,B,C,D,E,F and G as matrixpractice.jld
@eval using JLD
save("/matrixpractice.jld","A",A,"B",B,"C",C,"D", D,"E",E,"F",F,"G",G)
# i) Save matrices A-D as firstmatrix.jld
save("/firstmatrix.jld","A",A,"B",B,"C",C,"D", D)
# j) Export C as a csv file
@eval using CSV
@eval using DataFrames
C_1=DataFrame(C)
CSV.write("Cmatrix.csv", C_1)
# k) Export D as a tab-delimited file, called Dmatrix.dat
@eval using DelimitedFiles
writedlm("Dmatrix.dat", D)
return A, B, C, D
end
A, B, C, D = q1()

# 2 Loops and comprehensions
function q2(A,B,C)
    # a) Write a loop or a comprehension for the element-by-element product of A and B.
    AB=map((A,B) -> A.*B, A, B)
    # a) Create a matrix called AB2 that accomplishes the same without a comprehension or a loop
    AB2=A.*B
    # b) Write a loop that creates a column vector Cprime which contains the elements of C between -5 and 5
    Cprime=Float64[]
    for i in C
     if 5.0>=i>=-5.0
        append!(Cprime, i)
     end
    end
    # b) Create a vector Cprime2 that does the same without a loop
    Cprime2=filter(x -> -5.0<=x<=5.0,reshape(C,35,1))
    # c) Crate a 3-dim array called X of dimension NxKxT 15169x6x5
    n=15169
    T=5
    for i in 1:T
        local onesx=ones(n)
        local DM=rand(n)
        for p in eachindex(DM)
            if DM[p] > 0.75*(6-i)/5
                DM[p]=0
            else DM[p]=1
            end
        end
        local NV=rand(Normal(15.0+i-1,5.0*(i-1)),n,1)
        local NEV=rand(Normal(π*(6-i)/3,1/ℯ),n,1)
    # A discrete normal distribution is the limiting case when a binomial distribution has N increasing to infinity
    # For the case in which the mean is 12 and its standard deviation is 2.19^3 we can use the approximation of
    # a N->∞ binomial which mean is μ=Np and sd is σ^2=Np(1-p). So if we take as given mean at 12 and solve for p and then substitute in
    # σ^2=Np(1-p) and solve for N, we get that N=96.2219 and p=0.12471173 and q=1-p=0.875288266
    # Using these parameters (rounded) to approximate a discrete normal distribution through a binomial distribution:
        local DN=rand(Binomial(96,0.125),n,1)
        local BN=rand(Binomial(20,0.5),n,1)
        if i==1
        global    H=[onesx DM NV NEV DN BN]
        else
        global    H=[H;[onesx DM NV NEV DN BN]]
        end
    end
    X=reshape(H, (n,6,5))

    # d) β KxT matrix with evolution through time is as given. Use comprehensions.
    β=[[1+(i-1)*0.25 for i in 1:5] [log(i) for i in 1:5] [-√i for i in 1:5] [ℯ^i - ℯ^(i+1)  for i in 1:5] [i  for i in 1:5] [i/3  for i in 1:5]]

    # e) Y is a NxT defined by Y_t=X_tB_t+e_t where e_t is N(0,σ=0.36)
    for i in 1:5
        if i==1
        global Y=X[:,:,1]*β[1,:]+rand(Normal(0,0.36),n,1)
        else
        global Y=[Y X[:,:,i]*β[i,:]+rand(Normal(0,0.36),n,1) ]
        end
    end
    end

q2(A,B,C)

function q3()
  # a) Clear workspace and read file,
   @eval using DataFrames
   @eval using CSV
   @eval using HTTP
   # I correct for misidentified data types
   nls= CSV.read(download("https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS1-julia-intro/nlsw88.csv"), missingstrings=["NA", "na", "n/a", "nothing","missing"],types=[Int,Int,String,Int,Int,Int,Int,Int,Int,Int,String,String,Int,Float64,Int,Float64,Float64], silencewarnings=true)
   DataFrame(nls)
   describe(nls)
   save("/nls88.jld","nls88",nls)
   # b) Percentage of sample never been married - no information is provided about marriage status
   mean(nls.never_married)*100
   # b) To compute percentage of college graduates, we use the variable collgrad
   mean(nls.collgrad)*100
   # c) What percentage of sample is in each race category?
   c=nrow(nls)
   count=by(nls,:race,nrow)
   count.percent=count.nrow/c*100
   count
   # d) summarystats is a matrix with mean, median, sd, min, max, unique, interquantile range of nls
   # I define the dataframe to work with
   summarystats=describe(nls, :mean, :median, :std, :min, :max, :nunique, :q75, :q25 )
   DataFrame(summarystats)
   # Since interquantile range is not provided directly by describe, we compute it skipping those non-numeric
   # variables
   c=nrow(summarystats)
   summarystats[:intqtile]=zeros(c)
   for i=1:c
    if summarystats[i,:q75]==nothing
    else summarystats[i,:intqtile]=summarystats[i,:q75]-summarystats[i,:q25]
    end
   end
  # e) cross-tabulation of industry and occupation
  count=by(nls,[:occupation, :industry],nrow)
  # count shows the cross-tabulation with levels of industrya in rows and levels of occupation as column names
  count=unstack(count, :occupation, :industry, :nrow)
  # f) Mean wage over industry and occupation categories
  nlsn=select!(nls,:industry, :occupation, :wage)
  # Mean wage over industry categories
  nlsna=groupby(nlsn,:industry)
  combine(nlsna, :wage => mean)
  # Mean wage over occupation categories
  nlsna=groupby(nlsn,:occupation)
  combine(nlsna, :wage => mean)
  # Mean wage over occupation-industry categories
  nlsna=groupby(nlsn,[:occupation,:industry])
  combine(nlsna, :wage => mean)
end
q3()

function q4()
    # a) Load firstmatrix.jld
   @eval using JLD
   @eval using LinearAlgebra
    fm=load("/firstmatrix.jld")
    # b) write a function called matrixops
    function matrixops(A,B)
    # c) This function computes and prints element by element product of both matrices AB, product of A'B (ApB) and the sum of each entry in A+B
        sa=size(A)
        sb=size(B)
     # e) Error if not the same size
        if sa!=sb
            println("inputs must have the same size")
        else
        AB=map((A,B) -> A.*B, A, B)
        # product A'B
        ApB=transpose(A)*B
        # sum of all elements of A+B
        sAB=sum(A+B)
        return AB, ApB, sAB
        end
    end
    # d) Evaluate using A and B
    matrixops(A,B)
    # f) Evaluate using C and D
    matrixops(C,D)
    # g) Evaluate matrixops with ttl_exp and wage
    @eval using CSV
    nls= CSV.read(download("https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS1-julia-intro/nlsw88.csv"), missingstrings=["NA", "na", "n/a", "nothing","missing"],types=[Int,Int,String,Int,Int,Int,Int,Int,Int,Int,String,String,Int,Float64,Int,Float64,Float64], silencewarnings=true)
    J=convert(Array,nls.ttl_exp)
    L=convert(Array,nls.wage)
    matrixops(J,L)
end
q4()


text:
module IsApprox

export AbstractApprox, Equal, EachApprox, Approx, UpToPhase
export ispossemidef, isunitary, isinvolution, isidempotent, isnormal, commutes, anticommutes

include("core.jl")
include("base_applications.jl")
include("other_applications.jl")

end # module


text:
using Text
using Base.Test, Stage, Ollam, DataStructures

const logger = Log(STDERR)

all       = lazy_map(l -> split(chomp(l), '\t')[2], filelines("data/saic-sms/all.tsv.gz"))
all_truth = map(l -> split(chomp(l), '\t')[1], filelines("data/saic-sms/all.tsv.gz"))

# -------------------------------------------------------------------------------------------------------------------------
# representative split
# -------------------------------------------------------------------------------------------------------------------------
xsection = DefaultDict(String, Vector{(String, String)}, () -> ((String, String))[])
for fvt in zip(all, all_truth)
  push!(xsection[fvt[2]], fvt)
end

total = 0
for k in keys(xsection)
  @info logger @sprintf("%-30s %10d", k, length(xsection[k]))
  total += length(xsection[k])
end
@sep logger
@info logger @sprintf("%-30s %10d", "total", total)

test        = vcat([ map(fvt -> fvt[1], xsection[k][1:int(length(xsection[k])*0.15)]) for k in keys(xsection) ]...)
test_truth  = vcat([ map(fvt -> fvt[2], xsection[k][1:int(length(xsection[k])*0.15)]) for k in keys(xsection) ]...)
train       = vcat([ map(fvt -> fvt[1], xsection[k][int(length(xsection[k])*0.15)+1:end]) for k in keys(xsection) ]...)
train_truth = vcat([ map(fvt -> fvt[2], xsection[k][int(length(xsection[k])*0.15)+1:end]) for k in keys(xsection) ]...)
@info logger "train: $(length(train)), test: $(length(test))"

# -------------------------------------------------------------------------------------------------------------------------
# LID
# -------------------------------------------------------------------------------------------------------------------------
bkgmodel, fextractor, model = lid_train(train, train_truth, lid_tokenizer,
                                        trainer = (fvs, truth, init_model) -> train_mira(fvs, truth, init_model, iterations = 20, k = 4, C = 0.01, average = true),
                                        iteration_method = :eager)
confmat = DefaultDict(String, DefaultDict{String, Int32}, () -> DefaultDict(String, Int32, 0))
@info @sprintf("mira test set error rate: %7.3f", test_classification(model, lazy_map(fextractor, test), test_truth, record = (t, h) -> confmat[t][h] += 1) * 100.0)
print_confusion_matrix(confmat)

bkgmodel, fextractor, model = lid_train(collect(train), collect(train_truth), lid_tokenizer,
                                        trainer = (fvs, truth, init_model) -> train_libsvm(fvs, truth, C = 1.0),
                                        iteration_method = :eager)
confmat = DefaultDict(String, DefaultDict{String, Int32}, () -> DefaultDict(String, Int32, 0))
@info @sprintf("svm test set error rate: %7.3f", test_classification(model, lazy_map(fextractor, test), test_truth, record = (t, h) -> confmat[t][h] += 1) * 100.0)
print_confusion_matrix(confmat)


text:
module CUSOLVER

import CUDAdrv
using CUDAdrv: CuContext, CuDevice
using CUDAnative

using ..CuArrays
const CuStream_t = Ptr{Nothing}

using ..CuArrays: libcusolver, active_context, _getindex

using LinearAlgebra

import Base.one
import Base.zero

include("libcusolver_types.jl")

const _dense_handles = Dict{CuContext,cusolverDnHandle_t}()
const _dense_handle = Ref{cusolverDnHandle_t}(C_NULL)

function dense_handle()
    if _dense_handle[] == C_NULL
        @assert isassigned(active_context) # some other call should have initialized CUDA
        _dense_handle[] = get!(_dense_handles, active_context[]) do
            context = active_context[]
            handle = cusolverDnCreate()
            atexit(()->CUDAdrv.isvalid(context) && cusolverDnDestroy(handle))
            handle
        end
    end

    return _dense_handle[]
end

include("error.jl")
include("libcusolver.jl")
include("dense.jl")
include("highlevel.jl")

end


text:
# crude method for binarization and being with limit
# they work great!
# need some crude  filters also :D

export boxconstraint!
export gradualbinarize!

function boxconstraint! end
function (sim::Simulation)(::typeof(boxconstraint!),m::AbstractArray,rng,maskF = EVERYWHERE ,gridtype::GridType = p̂)
	mask = maskF.(Coordinates(sim.grid,gridtype))
    for i in eachindex(m)
        if real(m[i]) < rng[1] && mask[i]  m[i] = rng[1] end
        if real(m[i]) > rng[2] && mask[i]  m[i] = rng[2] end
    end
end

function gradualbinarize! end

function (sim::Simulation)(::typeof(gradualbinarize!),m::AbstractArray,rng,δ,maskF = EVERYWHERE ,gridtype::GridType = p̂)
	mask = maskF.(Coordinates(sim.grid,gridtype))
    for i in eachindex(m)
        if real(m[i] - rng[1]) < δ && mask[i]  m[i] = rng[1] end
        if real(rng[2] - m[i]) < δ && mask[i]  m[i] = rng[2] end
    end
end


text:
function f()
	for x = 0:1
		y = 0
	end
	println(y)
end
f()

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
# biowulf.jl
# functions for use on the NIH Biowulf super computer

"""
makeswarm(;G::Int=2,cell="HCT116",swarmfile::String="fit",label="label",inlabel=label,timestamp="",nsets=1,datafolder::String="HCT116_testdata",fish= false,cycle=true,thresholdlow::Float64=0.,thresholdhigh::Float64=1e8,conds::String="MOCK",resultfolder::String= "fit_result",infolder=resultfolder,batchsize=1000,maxtime = 60.,nchains::Int=2,nthreads::Int=1,transient::Bool=false,fittedparam=collect(1:2*G-1),fixedeffects=(),juliafile::String="fitscript",root=".",samplesteps::Int=100000,warmupsteps=20000,annealsteps=0,temp=1.,tempanneal=100.,cv = 0.02,yieldprior=0.05)

makeswarm(genes::Vector;G::Int=2,cell="HCT116",swarmfile::String="fit",label="label",inlabel=label,timestamp="",nsets=1,datafolder::String="HCT116_testdata",fish=false,cycle=true,conds::String="MOCK",resultfolder::String="fit_result",infolder=resultfolder,batchsize=1000,maxtime=60.,nchains::Int=2,nthreads::Int=1,transient::Bool=false,fittedparam=collect(1:2*G-1),fixedeffects=(),juliafile::String="fitscript",root=".",samplesteps::Int=100000,warmupsteps=20000,annealsteps=0,temp=1.,tempanneal=100.,cv=0.02,yieldprior=0.05)

Arguments
    - `G`: number of gene states
    - `cell': cell type for halflives and allele numbers
    - `infolder`: name of folder for initial parameters
    - `swarmfile`: name of swarmfile to be executed by swarm
    - `label`: label of output files produced
    - `inlabel`: label of files used for initial conditions
    - `nsets`: number of histograms to be fit (e.g. one for wild type and one for perturbation)
    - `datafolder`: folder holding histograms, if two folders use `-` (hyphen) to separate, e.g.  "data\folder1-data\folder2"
    - `fish`: Data file type, set to true for FISH and false to scRNA (FISH type assumes different folder structure)
    - `thresholdlow`: lower threshold for halflife for genes to be fit
    - `threhsoldhigh`: upper threshold
    - `conds`: string describing conditions to be fit with `-` to separate if two conditions, e.g. "WT-AUXIN"
    - `result`: folder for results
    - `batchsize`: number of jobs per swarmfile, default = 1000
    - `maxtime`: maximum wall time for run, default = 2 hrs
    - `nchains`: number of MCMC chains = number of processors called by Julia, default = 2
    - 'nthreads`: number of Julia threads per processesor, default = 1
    - `transient::Bool`: true means fit transient model (T0, T30, T120)
    - `fittedparam`: vector of rate indices to be fit, e.g. [1,2,3,5,6,7]
    - `fixedeffects`: tuple of vectors of rates that are fixed between control and treatment where first index is fit and others are fixed to first, e.g. ([3,8],) means  index 8 is fixed to index 3
         (each vector in tuple is a fixed rate set)
    - `juliafile`: name of file to be called by julia in swarmfile
    - `root`: name of root directory for project, e.g. "scRNA\"
    - `samplesteps`: number of MCMC sampling steps
    - `warmupsteps`: number of MCMC warmup steps to find proposal distribution covariance
    - `annealsteps`: number of annealing steps (during annealing temperature is dropped from tempanneal to temp)
    - `temp`: MCMC temperature
    - `tempanneal`: annealing temperature
    - `cv`: coefficient of variation (mean/std) of proposal distribution
    - `yieldprior`: prior for yield, default = .05, (set to 1 for FISH if using scRNA data format for FISH data)

returns swarmfile that calls a julia file that is executed on biowulf

"""

function makeswarm(;G::Int=2,cell="HCT116",swarmfile::String="fit",label="label",inlabel=label,timestamp="",nsets=1,datafolder::String="HCT116_testdata",fish= false,cycle=true,thresholdlow::Float64=0.,thresholdhigh::Float64=1e8,conds::String="MOCK",resultfolder::String= "fit_result",infolder=resultfolder,batchsize=1000,maxtime = 60.,nchains::Int=2,nthreads::Int=1,transient::Bool=false,fittedparam=collect(1:2*G-1),fixedeffects=(),juliafile::String="fitscript",root=".",samplesteps::Int=100000,warmupsteps=20000,annealsteps=0,temp=1.,tempanneal=100.,cv = 0.02,yieldprior=0.05)
    if occursin.("-",conds)
        cond = string.(split(conds,"-"))
    else
        cond = conds
    end
    genes = checkgenes(root,cond,datafolder,cell,thresholdlow,thresholdhigh)
    makeswarm(genes,G=G,cell=cell,infolder=infolder,swarmfile=swarmfile,label=label,inlabel=inlabel,timestamp=timestamp,nsets=nsets,datafolder=datafolder,fish=fish,cycle=cycle,conds=conds,resultfolder=resultfolder,batchsize=batchsize,maxtime=maxtime,nchains=nchains,nthreads=nthreads,transient=transient,fittedparam=fittedparam,fixedeffects=fixedeffects,juliafile=juliafile,root=root,samplesteps=samplesteps,warmupsteps=warmupsteps,annealsteps=annealsteps,temp=temp,tempanneal=tempanneal,cv=cv,yieldprior=yieldprior)
end

function makeswarm(genes::Vector;G::Int=2,cell="HCT116",swarmfile::String="fit",label="label",inlabel=label,timestamp="",nsets=1,datafolder::String="HCT116_testdata",fish=false,cycle=true,conds::String="MOCK",resultfolder::String="fit_result",infolder=resultfolder,batchsize=1000,maxtime=60.,nchains::Int=2,nthreads::Int=1,transient::Bool=false,fittedparam=collect(1:2*G-1),fixedeffects=(),juliafile::String="fitscript",root=".",samplesteps::Int=100000,warmupsteps=20000,annealsteps=0,temp=1.,tempanneal=100.,cv=0.02,yieldprior=0.05)
    if label == "label"
        if fish
            label = "FISH-ss"
        else
            label = "scRNA-ss"
        end
        label = label * "-" * cell
        if ~isempty(timestamp)
            label = label * "-" * timestamp
        end
        label = label * "_" * conds
        inlabel = label
    end
    ngenes = length(genes)
    println("number of genes: ",ngenes)
    juliafile = juliafile * "_" * label * "_" * "$G" * ".jl"
    if ngenes > batchsize
        batches = getbatches(genes,ngenes,batchsize)
        for batch in eachindex(batches)
            sfile = swarmfile * "_" * label * "_" * "$G" * "_" * "$batch" * ".swarm"
            write_swarmfile(sfile,nchains,nthreads,juliafile,batches[batch])
        end
    else
        sfile = swarmfile * "_" * label * "_" * "$G" * ".swarm"
        write_swarmfile(sfile,nchains,nthreads,juliafile,genes)
    end
    write_fitfile(juliafile,nchains,cell,conds,G,maxtime,fittedparam,fixedeffects,infolder,resultfolder,datafolder,fish,cycle,inlabel,label,nsets,transient,samplesteps,warmupsteps,annealsteps,temp,tempanneal,root,cv,yieldprior)
end

"""
    fix(folder)

    Finds jobs that failed and writes a swarmfile for those genes

"""
fix(folder) = writeruns(fixruns(findjobs(folder)))

function fix_filenames(folder,old="scRNA-ss-",new="scRNA-ss_")
    files = readdir(folder)
    for file in files
        if occursin(old,file)
            nfile = replace(file,old => new)
            mv(joinpath(folder,file),joinpath(folder,nfile),force=true)
        end
    end
end

"""
    setup(rootfolder = "scRNA")

    Sets up the folder system prior to use
    Defaults to "scRNA"

"""
function rna_setup(root = ".")

    data = joinpath(root,"data")
    results = joinpath(root,"results")
    alleles = joinpath(data,"alleles")
    halflives = joinpath(data,"halflives")
    testdata = joinpath(data,"HCT116_testdata")

    if ~ispath(data)
        mkpath(data)
    end
    if ~ispath(results)
        mkpath(results)
    end
    if ~ispath(alleles)
        mkpath(alleles)
    end
    Downloads.download("https://raw.githubusercontent.com/nih-niddk-mbs/StochasticGene.jl/master/data/alleles/CH12_alleles.csv","$alleles/CH12_alleles.csv")
    Downloads.download("https://raw.githubusercontent.com/nih-niddk-mbs/StochasticGene.jl/master/data/alleles/HCT116_alleles.csv","$alleles/HCT116_alleles.csv")
    if ~ispath(halflives)
        mkpath(halflives)
    end
    Downloads.download("https://raw.githubusercontent.com/nih-niddk-mbs/StochasticGene.jl/master/data/halflives/ESC_halflife.csv","$halflives/ESC_halflife.csv")
    Downloads.download("https://raw.githubusercontent.com/nih-niddk-mbs/StochasticGene.jl/master/data/halflives/CH12_halflife.csv","$halflives/CH12_halflife.csv")
    Downloads.download("https://raw.githubusercontent.com/nih-niddk-mbs/StochasticGene.jl/master/data/halflives/HCT116_halflife.csv","$halflives/HCT116_halflife.csv")
    Downloads.download("https://raw.githubusercontent.com/nih-niddk-mbs/StochasticGene.jl/master/data/halflives/OcaB_halflife.csv","$halflives/OcaB_halflife.csv")
    Downloads.download("https://raw.githubusercontent.com/nih-niddk-mbs/StochasticGene.jl/master/data/halflives/aB_halflife.csv","$halflives/aB_halflife.csv")
    if ~ispath(testdata)
        mkpath(testdata)
    end
    Downloads.download("https://raw.githubusercontent.com/nih-niddk-mbs/StochasticGene.jl/master/data/HCT116_testdata/CENPL_MOCK.txt","$testdata/CENPL_MOCK.txt")
    Downloads.download("https://raw.githubusercontent.com/nih-niddk-mbs/StochasticGene.jl/master/data/HCT116_testdata/MYC_MOCK.txt","$testdata/MYC_MOCK.txt")

end

"""
make_fitfile(fitfile,fittedparam,fixedeffects)

make the file the swarm file calls to execute julia code

"""
function write_fitfile(fitfile,nchains,cell,datacond,G,maxtime,fittedparam,fixedeffects,infolder,resultfolder,datafolder,fish,cycle,inlabel,label,nsets,runcycle,samplesteps,warmupsteps,annealsteps,temp,tempanneal,root,cv,yieldprior)
        f = open(fitfile,"w")
        s =   '"'
        write(f,"@everywhere using StochasticGene\n")
        write(f,"@time fit_rna($nchains,ARGS[1],$s$cell$s,$fittedparam,$fixedeffects,$s$datacond$s,$G,$maxtime,$s$infolder$s,$s$resultfolder$s,$s$datafolder$s,$fish,$cycle,$s$inlabel$s,$s$label$s,$nsets,$cv,$runcycle,$samplesteps,$warmupsteps,$annealsteps,$temp,$tempanneal,$s$root$s,$yieldprior)\n")
        close(f)
end

function write_fitfile_include(fitfile,nchains,cell,datacond,G,maxtime,fittedparam,fixedeffects,infolder,resultfolder,datafolder,fish,cycle,inlabel,label,nsets,runcycle,samplesteps,warmupsteps,annealsteps,temp,tempanneal,root,modulepath,cv,yieldprior)
        f = open(fitfile,"w")
        s =   '"'
        write(f,"@everywhere include($s$modulepath$s)\n")
        # write(f,"@everywhere using StochasticGene\n")
        write(f,"@time StochasticGene.fit_rna($nchains,ARGS[1],$s$cell$s,$fittedparam,$fixedeffects,$s$datacond$s,$G,$maxtime,$s$infolder$s,$s$resultfolder$s,$s$datafolder$s,$fish,$cycle,$s$inlabel$s,$s$label$s,$nsets,$cv,$runcycle,$samplesteps,$warmupsteps,$annealsteps,$temp,$tempanneal,$s$root$s,$yieldprior)\n")
        close(f)
end

function getbatches(genes,ngenes,batchsize)
    nbatches = div(ngenes,batchsize)
    batches = Vector{Vector{String}}(undef,nbatches+1)
    println(batchsize," ",nbatches+1)
    for i in 1:nbatches
        batches[i] = genes[batchsize*(i-1)+1:batchsize*(i)]
    end
    batches[end] = genes[batchsize*nbatches+1:end]
    return batches
end

function write_swarmfile(sfile,nchains,nthreads,juliafile,genes::Vector)
    f = open(sfile,"w")
    for gene in genes
        gene = check_genename(gene,"(")
        writedlm(f,["julia -t $nthreads -p" nchains juliafile gene])
        # writedlm(f,["julia -p" nchains juliafile nchains gene cell cond G maxtime infolder resultfolder datafolder fish inlabel label nsets runcycle transient fittedparam fixedeffects])
    end
    close(f)
end

function checkgenes(root,conds::Vector,datafolder,celltype::String,thresholdlow::Float64,thresholdhigh::Float64)
    genes = Vector{Vector}(undef,2)
    if occursin.("-",datafolder)
        datafolder = string.(split(datafolder,"-"))
        for i in 1:2
            genes[i] = checkgenes(root,conds[i],datafolder[i],celltype,thresholdlow,thresholdhigh)
        end
    else
        for i in 1:2
            genes[i] = checkgenes(root,conds[i],datafolder,celltype,thresholdlow,thresholdhigh)
        end
    end
    intersect(genes[1],genes[2])
end

function checkgenes(root,cond::String,datafolder,cell::String,thresholdlow::Float64,thresholdhigh::Float64)
    datafolder = joinpath("data",datafolder)
    genes = intersect(get_halflives(root,cell,thresholdlow,thresholdhigh), get_genes(root,cond,datafolder))
    alleles = get_alleles(root,cell)
    if ~isnothing(alleles)
        return intersect(genes,alleles)
    else
        return genes
    end
end

function checkgenes(cond::String,datafolder,cell::String,thresholdlow::Float64,thresholdhigh::Float64)
    genes = intersect(get_halflives(".",cell,thresholdlow,thresholdhigh), get_genes(cond,datafolder))
    alleles = get_alleles(".",cell)
    if ~isnothing(alleles)
        return intersect(genes,alleles)
    else
        return genes
    end
end

function get_halflives(root,cell,thresholdlow::Float64,thresholdhigh::Float64)
    file = get_file(root,"data/halflives",cell,"csv")
    get_halflives(file,thresholdlow,thresholdhigh)
end

function get_halflives(file,thresholdlow::Float64,thresholdhigh::Float64)
    genes = Vector{String}(undef,0)
    halflives = readdlm(file,',')
    for row in eachrow(halflives)
        if typeof(row[2]) <: Number
            if thresholdlow <= row[2] < thresholdhigh
                push!(genes,string(row[1]))
            end
        end
    end
    return genes
end

function get_alleles(root,cell)
    file = get_file(root,"data/alleles",cell,"csv")
    if ~isnothing(file)
        return readdlm(file,',')[2:end,1]
    else
        return nothing
    end
end

function get_file(root,folder,type,suffix)
    folder = joinpath(root,folder)
    if ispath(folder)
        files = readdir(folder)
        for file in files
            name = split(file,"_")[1]
            if occursin(suffix,file) && name == type
                path = joinpath(folder,file)
                return path
            end
        end
    else
        throw("folder $folder does not exist")
    end
    nothing
end

function findjobs(folder)
    files = readdir(folder)
    files = files[occursin.("swarm_",files)]
    for (i,file) in enumerate(files)
        files[i] = split(file,"_")[2]
    end
    unique(files)
end

function fixruns(jobs,message="FAILED")
    runlist = Vector{String}(undef,0)
    for job in jobs
        if occursin(message,read(`jobhist $job`,String))
            swarmfile = findswarm(job)
            list = readdlm(swarmfile,',')
            runs =  chomp(read(pipeline(`jobhist $job`, `grep $message`),String))
            runs = split(runs,'\n')
            println(job)
            for run in runs
                linenumber = parse(Int,split(split(run," ")[1],"_")[2]) + 1
                while linenumber < length(list)
                    a = String(list[linenumber])
                    linenumber += 1000
                    println(a)
                    push!(runlist,a)
                end
            end
        end
    end
    return runlist
end

function writeruns(runs,outfile="fitfix.swarm")

    f = open(outfile,"w")
    for run in runs
        writedlm(f,[run],quotes=false)
    end
    close(f)

end

function findswarm(job)
    sc = "Swarm Command"
    line = read(pipeline(`jobhist $job`, `grep $sc`),String)
    list = split(line," ")
    list[occursin.(".swarm",list)][1]
end

function collate(folders,type="rates_scRNA_T120")
    genes = Vector{Vector}(undef,0)
    for folder in folders
        files = readdir(folder)
        for file in files
            if occursin(type,file)
                contents=readdlm(file,',')
                for row in eachrow(contents)
                    row[1] .== genes
                end
            end
        end
    end
end

"""
    get_missing_genes(datafolder::String,folder::String,cell,type,label,cond,model)
"""

function get_missing_genes(datafolder::String,folder::String,cell,type,label,cond,model,root=".")
    genes = checkgenes(root,cond,datafolder,cell,0.,1e8)
    genes1=get_genes(folder,type,label,cond,model)
    union(setdiff(genes1,genes),setdiff(genes,genes1))
end

function get_missing_genes(genes::Vector,folder,type,label,cond,model)
    genes1=get_genes(folder,type,label,cond,model)
    union(setdiff(genes1,genes),setdiff(genes,genes1))
end

get_missing_genes(genes,genes1) = union(setdiff(genes1,genes),setdiff(genes,genes1))

function scan_swarmfiles(jobid,folder=".")
    if ~(typeof(jobid) <: String)
        jobid = string(jobid)
    end
    genes = Array{String,1}(undef,0)
    files = readdir(folder)
    files = files[occursin.(jobid,files)]
    for file in files
        genes = vcat(genes,scan_swarmfile(file))
    end
    return genes
end

function scan_swarmfile(file)
    genes = Array{String,1}(undef,0)
    contents = readdlm(file,'\t')
    lines = contents[occursin.("[\"",string.(contents))]
    for line in lines
        push!(genes,split.(string(line)," ")[1])
    end
    return genes
end

function scan_fitfile(file,folder=".")
    genes = Array{String,1}(undef,0)
    joinpath(folder,file)
    file = readdlm(file,'\t')
    for line in eachrow(file)
        push!(genes,line[4])
    end
    return genes
end


text:
module Requests

using Genie, Genie.Router
using HTTP

export jsonpayload, rawpayload, filespayload, postpayload, getpayload, getrequest

function jsonpayload()
  @params(Genie.PARAMS_JSON_PAYLOAD)
end

function rawpayload()
  @params(Genie.PARAMS_RAW_PAYLOAD)
end

function filespayload()
  @params(Genie.PARAMS_FILES)
end
function filespayload(filename::String)
  @params(Genie.PARAMS_FILES)[filename]
end

function postpayload()
  @params(Genie.PARAMS_POST_KEY)
end

function getpayload()
  @params(Genie.PARAMS_GET_KEY)
end

function getrequest()
  @params(Genie.PARAMS_REQUEST_KEY)
end

end


text:
export RedPitayaCluster, master, readDataPeriods, numChan, readDataSlow

import Base: length

struct RedPitayaCluster
  rp::Vector{RedPitaya}
end

#TODO: set first RP to master
function RedPitayaCluster(hosts::Vector{String}, port=5025)
  # the first RP is the master
  rps = RedPitaya[ RedPitaya(host, port, i==1) for (i,host) in enumerate(hosts) ]

  for rp in rps
    triggerMode(rp, "EXTERNAL")
  end

  return RedPitayaCluster(rps)
end

length(rpc::RedPitayaCluster) = length(rpc.rp)
numChan(rpc::RedPitayaCluster) = 2*length(rpc)
master(rpc::RedPitayaCluster) = rpc.rp[1]

const _currentFrame = Ref(0)

function currentFrame(rpc::RedPitayaCluster)
  _currentFrame[] = currentFrame(rpc.rp[1]) #[ currentFrame(rp) for rp in rpc.rp ]
  @debug "Current frame: $(_currentFrame[])"
  #return minimum(currentFrames)
  return _currentFrame[]
end

function currentPeriod(rpc::RedPitayaCluster)
  currentPeriods = currentPeriod(rpc.rp[1])  #[ currentPeriod(rp) for rp in rpc.rp ]
  @debug "Current period: $currentPeriods"
  #return minimum(currentPeriods)
  return currentPeriods
end

for op in [:periodsPerFrame,  :samplesPerPeriod, :decimation, :keepAliveReset,
           :triggerMode, :slowDACPeriodsPerFrame]
  @eval $op(rpc::RedPitayaCluster) = $op(master(rpc))
  @eval begin
    function $op(rpc::RedPitayaCluster, value)
      for rp in rpc.rp
        $op(rp, value)
      end
    end
  end
end

for op in [:connectADC, :stopADC, :disconnect, :connect]
  @eval begin
    function $op(rpc::RedPitayaCluster)
      for rp in rpc.rp
        $op(rp)
      end
    end
  end
end

function startADC(rpc::RedPitayaCluster)
  wp = currentWP(rpc.rp[1])
  for rp in rpc.rp
    startADC(rp, wp)
  end
end

function masterTrigger(rpc::RedPitayaCluster, val::Bool)
    if val
        masterTrigger(master(rpc), val)
    else
        keepAliveReset(rpc, true)
        masterTrigger(master(rpc), false)
        keepAliveReset(rpc, false)
    end
end
masterTrigger(rpc::RedPitayaCluster) = masterTrigger(master(rpc))
bufferSize(rpc::RedPitayaCluster) = bufferSize(master(rpc))

# "TRIGGERED" or "CONTINUOUS"
function ramWriterMode(rpc::RedPitayaCluster, mode::String)
  for rp in rpc.rp
    ramWriterMode(rp, mode)
  end
end

for op in [:amplitudeDAC, :amplitudeDACNext, :frequencyDAC, :phaseDAC, :modulusFactorDAC, :modulusDAC]
  @eval function $op(rpc::RedPitayaCluster, chan::Integer, component::Integer)
    idxRP = div(chan-1, 2) + 1
    chanRP = mod1(chan, 2)
    return $op(rpc.rp[idxRP], chanRP, component)
  end
  @eval function $op(rpc::RedPitayaCluster, chan::Integer, component::Integer, value)
    idxRP = div(chan-1, 2) + 1
    chanRP = mod1(chan, 2)
    return $op(rpc.rp[idxRP], chanRP, component, value)
  end
end

for op in [:signalTypeDAC,  :DCSignDAC, :jumpSharpnessDAC]
  @eval function $op(rpc::RedPitayaCluster, chan::Integer)
    idxRP = div(chan-1, 2) + 1
    chanRP = mod1(chan, 2)
    return $op(rpc.rp[idxRP], chanRP)
  end
  @eval function $op(rpc::RedPitayaCluster, chan::Integer, value)
    idxRP = div(chan-1, 2) + 1
    chanRP = mod1(chan, 2)
    return $op(rpc.rp[idxRP], chanRP, value)
  end
end

function setSlowDAC(rpc::RedPitayaCluster, chan, value)
  idxRP = div(chan-1, 2) + 1
  chanRP = mod(chan-1, 2)
  setSlowDAC(rpc.rp[idxRP], chanRP, value)
end

function getSlowADC(rpc::RedPitayaCluster, chan::Integer)
  idxRP = div(chan-1, 2) + 1
  chanRP = mod(chan-1, 2)
  getSlowADC(rpc.rp[idxRP], chanRP)
end

function numSlowADCChan(rpc::RedPitayaCluster)
  tmp = [ numSlowADCChan(rp) for rp in rpc.rp]
  return sum(tmp)
end

function numSlowADCChan(rpc::RedPitayaCluster, num)
  for rp in rpc.rp
    numSlowADCChan(rp, num)
  end
  return
end

function passPDMToFastDAC(rpc::RedPitayaCluster, val::Vector{Bool})
  for (d,rp) in enumerate(rpc.rp)
    passPDMToFastDAC(rp, val[d])
  end
end

function passPDMToFastDAC(rpc::RedPitayaCluster)
  return [ passPDMToFastDAC(rp) for rp in rpc.rp]
end

modeDAC(rpc::RedPitayaCluster) = modeDAC(master(rpc))

function modeDAC(rpc::RedPitayaCluster, mode::String)
  for rp in rpc.rp
    modeDAC(rp, mode)
  end
end

function enableSlowDAC(rpc::RedPitayaCluster, enable::Bool, numFrames::Int64=0,
                       ffRampUpTime::Float64=0.4, ffRampUpFraction::Float64=0.8)
  # We just use the first rp currently
  #res = [enableSlowDAC(rp, enable) for rp in rpc.rp]
  #return maximum(res)
  return enableSlowDAC(rpc.rp[1], enable, numFrames, ffRampUpTime, ffRampUpFraction)
end

function slowDACInterpolation(rpc::RedPitayaCluster, enable::Bool)
  for rp in rpc.rp
    slowDACInterpolation(rp, enable)
  end
end

# High level read. numFrames can adress a future frame. Data is read in
# chunks
function readData(rpc::RedPitayaCluster, startFrame, numFrames, numBlockAverages=1, numPeriodsPerPatch=1)
  numSampPerPeriod = master(rpc).samplesPerPeriod
  numSamp = numSampPerPeriod * numFrames
  numPeriods = master(rpc).periodsPerFrame
  numSampPerFrame = numSampPerPeriod * numPeriods

  if rem(numSampPerPeriod,numBlockAverages) != 0
    error("block averages has to be a divider of numSampPerPeriod")
  end

  numTrueSampPerPeriod = div(numSampPerPeriod,numBlockAverages*numPeriodsPerPatch)

  data = zeros(Float32, numTrueSampPerPeriod, numChan(rpc), numPeriods*numPeriodsPerPatch, numFrames)
  wpRead = startFrame
  l=1

  numFramesInMemoryBuffer = bufferSize(master(rpc)) / numSampPerFrame #numSamp
  @debug "numFramesInMemoryBuffer = $numFramesInMemoryBuffer"

  # This is a wild guess for a good chunk size
  chunkSize = max(1,  round(Int, 1000000 / numSampPerFrame)  )
  @debug "chunkSize = $chunkSize"
  while l<=numFrames
    wpWrite = currentFrame(rpc)
    while wpRead >= wpWrite # Wait that startFrame is reached
      wpWrite = currentFrame(rpc)
      @debug wpWrite
    end
    chunk = min(wpWrite-wpRead,chunkSize) # Determine how many frames to read
    @debug chunk
    if l+chunk > numFrames
      chunk = numFrames - l + 1
    end

    if wpWrite - numFramesInMemoryBuffer >= wpRead
      @error("WARNING: We have lost data !!!!!!!!!!")
    end

    @debug "Read from $wpRead until $(wpRead+chunk-1), WpWrite $(wpWrite), chunk=$(chunk)"

    #@sync begin
    done = zeros(Bool, length(rpc.rp))
    @async for (d,rp) in enumerate(rpc.rp)

      u = readData_(rp, Int64(wpRead), Int64(chunk))
      utmp1 = reshape(u,2,numTrueSampPerPeriod,numBlockAverages,
                        size(u,3)*numPeriodsPerPatch,size(u,4))
      utmp2 = numBlockAverages > 1 ? mean(utmp1,dims=3) : utmp1

      data[:,2*d-1,:,l:(l+chunk-1)] = utmp2[1,:,1,:,:]
      data[:,2*d,:,l:(l+chunk-1)] = utmp2[2,:,1,:,:]
      done[d] = true
    end
    timeout = 10
    t = Timer(timeout)
    while isopen(t)
      sleep(0.001)
      if all(done)
        break
      end
    end
    if !all(done)
      error("Timout reached when reading from sockets")
    end
    #end

    l += chunk
    wpRead += chunk
  end

  return data
end

function readDataPeriods(rpc::RedPitayaCluster, startPeriod, numPeriods, numBlockAverages=1)
  numSampPerPeriod = master(rpc).samplesPerPeriod
  numSamp = numSampPerPeriod * numPeriods

  if rem(numSampPerPeriod,numBlockAverages) != 0
    error("block averages has to be a divider of numSampPerPeriod")
  end

  numAveragedSampPerPeriod = div(numSampPerPeriod,numBlockAverages)

  data = zeros(Float32, numAveragedSampPerPeriod, numChan(rpc), numPeriods)
  wpRead = startPeriod
  l=1

  # This is a wild guess for a good chunk size
  chunkSize = max(1,  round(Int, 1000000 / numSampPerPeriod)  )
  println("chunkSize = $chunkSize; numPeriods = $numPeriods")
  while l<=numPeriods
    wpWrite = currentPeriod(rpc)
    while wpRead >= wpWrite # Wait that startPeriod is reached
      wpWrite = currentPeriod(rpc)
      @debug wpWrite
    end
    chunk = min(wpWrite-wpRead,chunkSize) # Determine how many periods to read
    @debug chunk
    if l+chunk > numPeriods
      chunk = numPeriods - l + 1
    end

    @debug "Read from $wpRead until $(wpRead+chunk-1), WpWrite $(wpWrite), chunk=$(chunk)"

    for (d,rp) in enumerate(rpc.rp)
    # @sync   @async begin
        u = readDataPeriods_(rp, Int64(wpRead), Int64(chunk))

        utmp1 = reshape(u,2,numAveragedSampPerPeriod,numBlockAverages,size(u,3))
        utmp2 = numBlockAverages > 1 ? mean(utmp1,dims=3) : utmp1

        data[:,2*d-1,l:(l+chunk-1)] = utmp2[1,:,1,:]
        data[:,2*d,l:(l+chunk-1)] = utmp2[2,:,1,:]
    #  end
    end

    l += chunk
    wpRead += chunk
  end

  return data
end

# High level read. numFrames can adress a future frame.
function readDataSlow(rpc::RedPitayaCluster, startFrame, numFrames)
  numPeriods = master(rpc).periodsPerFrame
  numChanTotal = numSlowADCChan(rpc)
  numSampPerFrame = numPeriods * numChanTotal

  data = zeros(Float32, numChanTotal, numPeriods, numFrames)
  wpRead = startFrame
  l=1

  # This is a wild guess for a good chunk size
  chunkSize = max(1,  round(Int, 1000000 / numSampPerFrame)  )
  @debug "chunkSize = $chunkSize"
  while l<=numFrames
    wpWrite = currentFrame(rpc)
    while wpRead >= wpWrite # Wait that startFrame is reached
      wpWrite = currentFrame(rpc)
      @debug wpWrite
    end
    chunk = min(wpWrite-wpRead,chunkSize) # Determine how many frames to read
    @debug chunk
    if l+chunk > numFrames
      chunk = numFrames - l + 1
    end

    @debug "Read from $wpRead until $(wpRead+chunk-1), WpWrite $(wpWrite), chunk=$(chunk)"

    p = 1
    for (d,rp) in enumerate(rpc.rp)
      u = readDataSlow_(rp, Int64(wpRead), Int64(chunk))
      numChan = size(u,1)

      #utmp1 = reshape(u, numChan, numBlockAverages,
      #                    div(size(u,2),numBlockAverages),size(u,3))
      #utmp2 = numBlockAverages > 1 ? mean(utmp1,dims=2) : utmp1

      data[p:(p+numChan-1),:,l:(l+chunk-1)] = u #utmp2[:,1,:,:]
      p += numChan
    end

    l += chunk
    wpRead += chunk
  end

  return data
end

max_stars_repo_name: fredrikekre/RedPitayaDAQServer

max_issues_repo_name: fuerpy/RedPitayaDAQServer

max_forks_repo_name: fuerpy/RedPitayaDAQServer


text:
using Dates
using Hospitalplanning
HP = Hospitalplanning

path = "./test/Sample data/PatientOverview_test.xlsx"
sheet = "Sheet1"
startdate = Date("2019-01-01")
enddate = Date("2019-02-01")
mastercalendar = HP.MasterCalendar(startdate,enddate)

columns = Dict(:Consultation => "Consultation" , :Telefon =>"Telephone" , :TTE => "TTE", :AEKG => "AEKG",  :MR => "MR", :Holter=>"Holter")
patients, visits = HP.readPatientTable(path,sheet,columns,mastercalendar)

path_resourceOverview = "./test/Sample data/GUCHamb_Timeslot_test.xlsx"
sheet_amb = "GUCH AMB"
sheet_resources = "External"
GUCHAmb_resources,timeslots,workpattern = Hospitalplanning.readWorkPattern(path_resourceOverview,sheet_amb)
external_resources = Hospitalplanning.readWorkPattern!(GUCHAmb_resources,timeslots,workpattern,path_resourceOverview,sheet_resources,)
path_TimeDelta = "./test/Sample data/GUCHamb_timeDelta.csv"
timeDelta = HP.loadTimeDelta(path_TimeDelta)
HP.generateCalendarFromPattern!(timeslots,GUCHAmb_resources,workpattern,mastercalendar,0.5)

mp, subs,sets = HP.columngeneration(patients,visits,GUCHAmb_resources,timeslots,mastercalendar,timeDelta; multithreading = false,setuponly=false)

solution = HP.extractsolution(mp,sets,timeslots,visits)
HP.save(solution,"./solution.csv")

import Pkg
const PT = Pkg.Types

Pkg.activate(pwd())             # current directory as the project
ctx = PT.Context()
pkg = ctx.env.pkg
if pkg ≡ nothing
    @error "Not in a package, I won't generate REQUIRE."
    exit(1)
else
    @info "found package" pkg = pkg
end

deps = PT.get_deps(ctx)
non_std_deps = sort(collect(setdiff(keys(deps), values(ctx.stdlibs))))

open("REQUIRE", "w") do io
    println(io, "julia 0.7")
    for d in non_std_deps
        println(io, d)
        @info "listing $d"
    end
end

#=-----------------------1--------------------

=#