text:
import Observables
import Makie

export Fiber
export update

"""
    Represents a Hopf fiber.

fields: point, s2tos2map, s2tos3map, radius, top, bottom, s3rotation, config, segments, color,
observable and scale.
"""
mutable struct Fiber <: Sprite
    point::S²
    s2tos3map::Any
    s2tos2map::Any
    radius::Float64
    top::S¹
    bottom::S¹
    s3rotation::S³
    config::Biquaternion
    segments::Int
    color::Observables.Observable{Array{Makie.ColorTypes.RGBA{Float32},2}}
    observable::Tuple{Observables.Observable{Array{Float64,2}},
                      Observables.Observable{Array{Float64,2}},
                      Observables.Observable{Array{Float64,2}}}
    scale::Real
end

"""
    Fiber(scene,
          point,
          s2tos3map,
          s2tos2map,
          [radius, [top, [bottom, [s3rotation, [config, [segments, [color, [transparency
           [,scale]]]]]]]]])

Construct a Hopf fiber with the given `scene`, `points` in the base space, `s2tos3map`,
`s2tos2map`, `radius`, `top`, `bottom`, S³ rotation `s3rotation`, configuration `config`,
the number of `segments`, `color`, `transparency` and `scale`.
"""
function Fiber(scene::Makie.Scene,
               point::S²,
               s2tos3map,
               s2tos2map;
               radius::Float64 = 0.02,
               top::S¹ = U1(-pi),
               bottom::S¹ = U1(pi),
               s3rotation::S³ = Quaternion(1, 0, 0, 0),
               config::Biquaternion = Biquaternion(ℝ³(0, 0, 0)),
               segments::Int = 36,
               color::Makie.RGBAf = Makie.RGBAf(1.0, 0.2705, 0.0, 0.5),
               transparency::Bool = false,
               scale::Real = 1.0)
    fiber = constructfiber(point,
                           s2tos3map,
                           s2tos2map,
                           radius = radius,
                           top = top,
                           bottom = bottom,
                           s3rotation = s3rotation,
                           config = config,
                           segments = segments,
                           scale = scale)
    colorarray = Observables.Observable(fill(color, size(fiber)...))
    observable = buildsurface(scene, fiber, colorarray, transparency = transparency)
    Fiber(point, s2tos3map, s2tos2map, radius, top, bottom, s3rotation, config, segments,
          colorarray, observable, scale)
end

"""
    update(fiber, points)

Update a Fiber by changing its observable with the given `fiber` and `points` in the base
space.
"""
function update(fiber::Fiber, point::S²)
    fiber.point = point
    value = constructfiber(fiber.point,
                           fiber.s2tos3map,
                           fiber.s2tos2map,
                           radius = fiber.radius,
                           top = fiber.top,
                           bottom = fiber.bottom,
                           s3rotation = fiber.s3rotation,
                           config = fiber.config,
                           segments = fiber.segments,
                           scale = fiber.scale)
    updatesurface(value, fiber.observable)
end

"""
    update(fiber, s2tos3map, s2tos2map)

Update a Fiber by changing its observable with the given `fiber` and map `s2tos3map` from
the base space into the total space, f: S² → S³, and also the map `s2tos2map` from the base
space into itself, f: S² → S².
"""
function update(fiber::Fiber, s2tos3map::Any, s2tos2map::Any)
    fiber.s2tos3map = s2tos3map
    fiber.s2tos2map = s2tos2map
    value = constructfiber(fiber.point,
                           fiber.s2tos3map,
                           fiber.s2tos2map,
                           radius = fiber.radius,
                           top = fiber.top,
                           bottom = fiber.bottom,
                           s3rotation = fiber.s3rotation,
                           config = fiber.config,
                           segments = fiber.segments,
                           scale = fiber.scale)
    updatesurface(value, fiber.observable)
end

"""
    update(fiber, top, bottom)

Update a Fiber by changing its observable with the given `fiber`, `top` and `bottom`.
"""
function update(fiber::Fiber, top::S¹, bottom::S¹)
    fiber.top = top
    fiber.bottom = bottom
    value = constructfiber(fiber.point,
                           fiber.s2tos3map,
                           fiber.s2tos2map,
                           radius = fiber.radius,
                           top = fiber.top,
                           bottom = fiber.bottom,
                           s3rotation = fiber.s3rotation,
                           config = fiber.config,
                           segments = fiber.segments,
                           scale = fiber.scale)
    updatesurface(value, fiber.observable)
end

"""
    update(fiber, s3rotation)

Update a Fiber by changing its observable with the given `fiber` and S³ rotation
`s3rotation`.
"""
function update(fiber::Fiber, s3rotation::S³)
    fiber.s3rotation = s3rotation
    value = constructfiber(fiber.point,
                           fiber.s2tos3map,
                           fiber.s2tos2map,
                           radius = fiber.radius,
                           top = fiber.top,
                           bottom = fiber.bottom,
                           s3rotation = fiber.s3rotation,
                           config = fiber.config,
                           segments = fiber.segments,
                           scale = fiber.scale)
    updatesurface(value, fiber.observable)
end

"""
    update(fiber, config)

Update a Fiber by changing its observable with the given `fiber` and configuration `config`.
"""
function update(fiber::Fiber, config::Biquaternion)
    fiber.config = config
    value = constructfiber(fiber.point,
                           fiber.s2tos3map,
                           fiber.s2tos2map,
                           radius = fiber.radius,
                           top = fiber.top,
                           bottom = fiber.bottom,
                           s3rotation = fiber.s3rotation,
                           config = fiber.config,
                           segments = fiber.segments,
                           scale = fiber.scale)
    updatesurface(value, fiber.observable)
end

"""
    update(fiber, color)

Update a Fiber by changing its observable with the given `fiber` and `color`.
"""
function update(fiber::Fiber, color::Makie.RGBAf)
    fiber.color[] = fill(color, fiber.segments, 10)
end

"""
    update(fiber, point, s3rotation, radius, top, bottom)

Update a Fiber by constructing it from scratch with the given base space `point`, S³
rotation `s3rotation`, sectional `radius`, beginning U(1) action `top` and ending U(1) action `bottom`.
"""
function update(fiber::Fiber, point::S², s3rotation::S³, radius::Float64, top::S¹, bottom::S¹)
    fiber.point = point
    fiber.s3rotation = s3rotation
    fiber.radius = radius
    fiber.top = top
    fiber.bottom = bottom
    value = constructfiber(fiber.point,
                           fiber.s2tos3map,
                           fiber.s2tos2map,
                           radius = fiber.radius,
                           top = fiber.top,
                           bottom = fiber.bottom,
                           s3rotation = fiber.s3rotation,
                           config = fiber.config,
                           segments = fiber.segments,
                           scale = fiber.scale)
    updatesurface(value, fiber.observable)
end

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
"""
Google Pub/Sub API
"""
module _pubsub

export pubsub

using ..api
using ...root

"""
Google Pub/Sub API root.
"""
pubsub = APIRoot(
    "https://pubsub.googleapis.com/v1/projects/{project}",
    Dict(
        "cloud-platform" => "Full access to all resources and services in the specified Cloud Platform project.",
        "pubsub" => "View and manage Pub/Sub topics and subscriptions",
    );
    Subscription=APIResource("subscriptions";
        acknowledge=APIMethod(:POST, "{subscription}:acknowledge", "Acknowledges the messages associated with the ack_ids in the AcknowledgeRequest."),
        create=APIMethod(:PUT, "{subscription}", "Creates a subscription to a given topic."),
        delete=APIMethod(:DELETE, "{subscription}", "Deletes an existing subscription."),
        get=APIMethod(:GET, "{subscription}", "Gets the configuration details of a subscription."),
        getIamPolicy=APIMethod(:GET, "{subscription}:getIamPolicy", "Gets the access control policy for a resource."),
        list=APIMethod(:GET, "", "Lists matching subscriptions."),
        modifyAckDeadline=APIMethod(:POST, "{subscription}:modifyAckDeadline", "Modifies the ack deadline for a specific message."),
        modifyPushConfig=APIMethod(:POST, "{subscription}:modifyPushConfig", "Modifies the PushConfig for a specified subscription."),
        pull=APIMethod(:POST, "{subscription}:pull", "Pulls messages from the server."),
        setIamPolicy=APIMethod(:POST, "{subscription}:setIamPolicy", "Sets the access control policy on the specified resource."),
        testIamPermissions=APIMethod(:POST, "{subscription}:testIamPermissions", "Returns permissions that a caller has on the specified resource."),
    ),
    Topic=APIResource("topics";
        create=APIMethod(:PUT, "{topic}", "Creates the given topic with the given name."),
        delete=APIMethod(:DELETE, "{topic}", "Deletes the topic with the given name."),
        get=APIMethod(:GET, "{topic}", "Gets the configuration of a topic."),
        getIamPolicy=APIMethod(:GET, "{topic}:getIamPolicy", "Gets the access control policy for a resource."),
        list=APIMethod(:GET, "", "Lists matching topics."),
        publish=APIMethod(:POST, "{topic}:publish", "Adds one or more messages to the topic."),
        setIamPolicy=APIMethod(:POST, "{topic}:setIamPolicy", "Sets the access control policy on the specified resource."),
        testIamPermissions=APIMethod(:POST, "{topic}:testIamPermissions", "Returns permissions that a caller has on the specified resource."),
    ),
    TopicSubscription=APIResource("topics/{topic}/subscriptions";
        list=APIMethod(:GET, "", "Lists the name of the subscriptions for this topic."),
    ),
)

end

ext: jl


text:
# Posterior c.d.f.s for simulation example for variable selection using power likelihood
include("varsel.jl")

module RunBeta
using VarSel
using PyPlot
draw() = (pause(0.001); get_current_fig_manager()[:window][:raise_]())

# Settings
ns = [10000]  # sample sizes to use
alphas = [Inf,50]  # values of alpha to use
nreps = 1  # number of times to repeat each simulation
n_total = 5*10^4  # total number of MCMC samples
n_keep = n_total  # number of MCMC samples to record
nburn = int(n_keep/10)  # burn-in

# Data distribution
using Skew
a = [0.6,2.7,-3.3,-4.9,-2.5]
Q = [1.0 -0.89 0.93 -0.91 0.98
     -0.89 1.0 -0.94 0.97 -0.91
     0.93 -0.94 1.0 -0.96 0.97
     -0.91 0.97 -0.96 1.0 -0.93
     0.98 -0.91 0.97 -0.93 1.0]
sigma = 1
p = length(a)+1

g(x0) = -1 + 4*(x0 + (1/16)*x0.^2)

function generate_sample(n)
    x = [ones(n) Skew.skewrndNormalized(n,Q,a)']
    x0 = vec(x[:,2])
    y = g(x0) + randn(n)*sigma
    return x,y
end

function interval(x,P)
    # Compute 100*P% interval from samples x
    N = length(x)
    xsort = sort(x)
    ai = (1-P)/2
    l,u = xsort[int(floor(ai*N))], xsort[int(ceil((1-ai)*N))]
    return l,u
end

# Run simulations
nns = length(ns)
for (i_a,alpha) in enumerate(alphas)
    for (i_n,n) in enumerate(ns)
        for rep in 1:nreps
            srand(n+rep) # Reset RNG

            # Sample data
            x,y = generate_sample(n)
            x0 = vec(x[:,2])

            # Run sampler
            zeta = (1/n)/(1/n + 1/alpha)
            beta_r,lambda_r,keepers = VarSel.run(x,y,n_total,n_keep,zeta)

            # Display c.d.f.s of coefficients
            xlims = Array[[-1.5,0.1],[2,4.5],[-.4,.4],[-.4,.4],[-.4,.4],[-.4,.4]]
            figure(i_a+2,figsize=(8,8)); clf(); hold(true)
            subplots_adjust(bottom=0.1,hspace=0.4)
            for j=1:p
                subplot(p,1,j)
                betas = vec(beta_r[j,nburn+1:end])
                n_use = n_keep - nburn
                plot(sort([betas,betas]),sort([0:n_use-1,1:n_use])/n_use,"b-",markersize=0,linewidth=2)
                ylim(0,1)
                xlim(xlims[j])
                xm,xM = xlim(); ym,yM = ylim()
                text((xM-xm)*0.92+xm,(yM-ym)*0.5+ym,"\$\\beta_$j\$",fontsize=18)
                l,u = interval(betas,0.95)
                plot([l,l],ylim(),"r-",markersize=0)
                plot([u,u],ylim(),"r-",markersize=0)
            end
            subplot(p,1,1)
            ylabel("\$\\\mathrm{c.d.f.}\$",fontsize=16)
            if alpha<Inf; title("\$\\mathrm{Coarsened}\$",fontsize=17)
            else; title("\$\\mathrm{Standard}\$",fontsize=17)
            end
            draw()
            savefig("varsel-beta-alpha=$alpha-n=$n-rep=$rep.png",dpi=150)
        end
    end
end

end # module

max_stars_repo_name: jwmi/CoarsenedPosterior

max_issues_repo_name: jwmi/CoarsenedPosterior

max_forks_repo_name: jwmi/CoarsenedPosterior


text:
function C3(t, fft_u, fft_v,
	    A1, A2,
	    matr, conjmatr,
	    sigma, llambda,
	    Ktaubis, epsilon, Ntaumm)

    h1u, h1v = ftau(t, fft_u, fft_v, A1, A2,
		    matr, conjmatr, sigma, llambda)

    h1u_fft = fft(h1u, dims=1)
    h1v_fft = fft(h1v, dims=1)

    dtu0u .= h1u_fft[1, :] / Ntaumm
    dtu0v .= h1v_fft[1, :] / Ntaumm
    h1u_fft[1, :] .= 0
    h1v_fft[1, :] .= 0
    h1u = epsilon * ifft(h1u_fft / (1im * Ktaubis), dims=1)
    h1v = epsilon * ifft(h1v_fft / (1im * Ktaubis), dims=1)

    champu1, champv1 = duftau(t, fft_u, fft_v,
			      dtu0u, dtu0v,
			      A1, A2,
			      matr, conjmatr,
			      sigma, llambda, Ntaumm)

    champu2, champv2 = dtftau(t, fft_u, fft_v,
			      A1, A2,
			      matr, conjmatr,
			      sigma, llambda, Ntaumm)

    champu = champu1 + champu2
    champv = champv1 + champv2

    champu_fft = fft(champu, dims=1)
    champv_fft = fft(champv, dims=1)

    dttu0u = champu_fft[1, :] / Ntaumm
    dttu0v = champv_fft[1, :] / Ntaumm

    champu_fft[1, :] .= 0
    champv_fft[1, :] .= 0

    dth1u = epsilon * ifft(champu_fft / (1im * Ktaubis), dims=1)
    dth1v = epsilon * ifft(champv_fft / (1im * Ktaubis), dims=1)

    champu, champv = ftau(t, fft_u .+ h1u, fft_v .+ h1v,
			  A1, A2,
			  matr, conjmatr,
			  sigma, llambda)

    h2u = champu - dth1u
    h2v = champv - dth1v

    h2u_fft = fft(h2u, dims=1)
    h2v_fft = fft(h2v, dims=1)

    champu_fft = fft(champu, dims=1)
    champv_fft = fft(champv, dims=1)

    dtu1u = champu_fft[1, :] / Ntaumm
    dtu1v = champv_fft[1, :] / Ntaumm

    h2u_fft[1, :] .= 0
    h2v_fft[1, :] .= 0

    h2u = epsilon * ifft(h2u_fft ./ (1im * Ktaubis), dims=1)
    h2v = epsilon * ifft(h2v_fft ./ (1im * Ktaubis), dims=1)

    champu1, champv1 = d2uftau(t, fft_u, fft_v,
			       dtu0u, dtu0v, dtu0u, dtu0v,
			       A1, A2, matr, conjmatr,
			       sigma, llambda, Ntaumm)

    champu2, champv2 = dtduftau(t, fft_u, fft_v,
				dtu0u, dtu0v,
				A1, A2, matr, conjmatr,
				sigma, llambda, Ntaumm)

    champu3, champv3 = duftau(t, fft_u, fft_v,
			      dttu0u, dttu0v,
			      A1, A2, matr, conjmatr,
			      sigma, llambda, Ntaumm)

    champu4, champv4 = d2tftau(t, fft_u, fft_v,
			       A1, A2, matr, conjmatr,
			       sigma, llambda, Ntaumm)

    champu = champu1 .+ 2 * champu2 .+ champu3 .+ champu4
    champv = champv1 .+ 2 * champv2 .+ champv3 .+ champv4

    champu_fft = fft(champu, dims=1)
    champv_fft = fft(champv, dims=1)

    champu_fft[1, :] .= 0
    champv_fft[1, :] .= 0

    dtth1u = epsilon * ifft(champu_fft / (1im * Ktaubis), dims=1)
    dtth1v = epsilon * ifft(champv_fft / (1im * Ktaubis), dims=1)

    champu1, champv1 = duftau(t, fft_u .+ h1u, fft_v .+ h1v,
                              dtu1u .+ dth1u, dtu1v .+ dth1v,
			      A1, A2, matr,
                              conjmatr, sigma, llambda, Ntaumm)

    champu2, champv2 = dtftau(t, fft_u .+ h1u, fft_v .+ h1v,
			      A1, A2, matr,
                              conjmatr, sigma, llambda, Ntaumm)

    champu = champu1 + champu2 - dtth1u
    champv = champv1 + champv2 - dtth1v

    champu_fft = fft(champu, dims=1)
    champv_fft = fft(champv, dims=1)

    champu_fft[1, :] .= 0
    champv_fft[1, :] .= 0

    dth2u = epsilon * ifft(champu_fft ./ (1im * Ktaubis), dims=1)
    dth2v = epsilon * ifft(champv_fft ./ (1im * Ktaubis), dims=1)

    h3u, h3v = ftau(t, fft_u .+ h2u, fft_v .+ h2v,
		    A1, A2, matr, conjmatr,
                    sigma, llambda)

    h3u = h3u - dth2u
    h3v = h3v - dth2v
    h3u_fft = fft(h3u, dims=1)
    h3v_fft = fft(h3v, dims=1)
    h3u_fft[1, :] .= 0
    h3v_fft[1, :] .= 0
    h3u = epsilon * ifft(h3u_fft / (1im * Ktaubis), dims=1)
    h3v = epsilon * ifft(h3v_fft / (1im * Ktaubis), dims=1)

    C3u = fft_u .+ h3u
    C3v = fft_v .+ h3v

    C3u, C3v

end


text:
using POMDPs
using ContinuousPOMDPTreeSearchExperiments
using POMDPToolbox
using POMCPOW
using VDPTag
using MCTS
using DataFrames
using ParticleFilters
using CPUTime

@everywhere using POMDPToolbox
@everywhere using CPUTime

N = 1000

pomdp = VDPTagPOMDP()

function create_pft(m)
    rng = MersenneTwister(13)
    rollout_policy = ToNextML(mdp(pomdp), rng)
    node_updater = ObsAdaptiveParticleFilter(pomdp, LowVarianceResampler(m), 0.05, rng)
    solver = DPWSolver(n_iterations=typemax(Int),
                       exploration_constant=40.0,
                       depth=10,
                       k_action=8.0,
                       alpha_action=1/20,
                       k_state=4.0,
                       alpha_state=1/20,
                       check_repeat_state=false,
                       check_repeat_action=false,
                       estimate_value=FORollout(rollout_policy),
                       # next_action=NextMLFirst(mdp(pomdp), rng),
                       next_action=RootToNextMLFirst(rng),
                       # default_action=TagAction(false,0.0),
                       rng=rng
                      )
    belief_mdp = GenerativeBeliefMDP(deepcopy(pomdp), node_updater)
    planner = solve(solver, belief_mdp)
    d = Dict(:cpu_us=>Int[],
             :n_nodes=>Int[],
            )
    return PolicyWrapper(planner, payload=d) do p, d, b
        start_us = CPUtime_us()
        a = action(p, b)
        cpu_us = CPUtime_us() - start_us
        push!(d[:cpu_us], cpu_us)
        tree = get(p.tree)
        push!(d[:n_nodes], length(tree.total_n))
        return a
    end
end

wrapped = Dict{String, Union{Solver,Policy}}(

    "pomcpow" => begin
        rng = MersenneTwister(13)
        rollout_policy = ToNextML(mdp(pomdp), rng)
        solver = POMCPOWSolver(tree_queries=10_000_000,
                               criterion=MaxUCB(40.0),
                               final_criterion=MaxTries(),
                               max_depth=10,
                               k_action=8.0,
                               alpha_action=1/20,
                               k_observation=4.0,
                               alpha_observation=1/20,
                               estimate_value=FORollout(rollout_policy),
                               check_repeat_act=false,
                               check_repeat_obs=false,
                               # next_action=NextMLFirst(mdp(pomdp), MersenneTwister(15)),
                               next_action=RootToNextMLFirst(rng),
                               default_action=TagAction(false,0.0),
                               rng=rng,
                              )
        planner = solve(solver, deepcopy(pomdp))
        d = Dict(:cpu_us=>Int[],
                 :n_nodes=>Int[],
                )
        PolicyWrapper(planner, payload=d) do p, d, b
            start_us = CPUtime_us()
            a = action(p, b)
            cpu_us = CPUtime_us() - start_us
            push!(d[:cpu_us], cpu_us)
            tree = get(p.tree)
            push!(d[:n_nodes], length(tree.total_n))
            return a
        end
    end,

    "pft_5" => create_pft(5),
    "pft_50" => create_pft(50),
    # "pft_100" => create_pft(100),
    # "pft_1000" => create_pft(1000)
)

alldata = DataFrame()
for t in logspace(-2,1,7)
# for t in logspace(-1,1,5)
# for t in [10.0]
    for (k, wrapper) in wrapped
        println("$k ($t)")
        wrapper.policy.solver.max_time = t
        sims = []
        for i in 1:N
            srand(wrapper, i+50_000)
            sim = Sim(deepcopy(pomdp),
                      deepcopy(wrapper),
                      ObsAdaptiveParticleFilter(deepcopy(pomdp),
                                                LowVarianceResampler(10_000),
                                                0.05, MersenneTwister(i+90_000)),
                      rng=MersenneTwister(i+70_000),
                      max_steps=100,
                      metadata=Dict(:solver=>k,
                                    :time=>t,
                                    :i=>i)
                     )
            push!(sims, sim)
        end

        data = run_parallel(sims) do sim, h
            stuff = sim.metadata
            if isa(sim.policy, PolicyWrapper)
                p = sim.policy
                if !isempty(p.payload[:cpu_us])
                    stuff[:mean_cpu_us] = mean(p.payload[:cpu_us])
                    stuff[:max_cpu_us] = maximum(p.payload[:cpu_us])
                    stuff[:min_cpu_us] = minimum(p.payload[:cpu_us])
                end
                if !isempty(p.payload[:n_nodes])
                    stuff[:mean_nodes] = mean(p.payload[:n_nodes])
                    stuff[:max_nodes] = maximum(p.payload[:n_nodes])
                    stuff[:min_nodes] = minimum(p.payload[:n_nodes])
                end
            end
            stuff[:reward] = discounted_reward(h)
            return stuff
        end

        rs = data[:reward]
        println(@sprintf("reward: %6.3f ± %6.3f", mean(rs), std(rs)/sqrt(length(rs))))
        alldata = vcat(alldata, data)
    end
end

filename = Pkg.dir("ContinuousPOMDPTreeSearchExperiments", "data", "vdp_trends_$(Dates.format(now(), "E_d_u_HH_MM")).csv")
println("saving to $filename...")
writetable(filename, alldata)
println("done.")


text:
using Base.Meta
export @cps

const compops = [symbol("=="), :>, :<, :>=, :<=]

Atom = Union(Number, Bool, Symbol, String, QuoteNode)

function cps(expr, cont)
    strippedexpr = striplineinfo(expr)
    cpsexpr = tc(strippedexpr, cont)
    # Trampolining function definitions prevent them from been added
    # to the calling scope.
    if isfuncdef(strippedexpr)
        cpsexpr
    else
        :(Stochy.trampoline(Stochy.Thunk(() -> $cpsexpr)))
    end
end

# This assumes line numbers have been stripped.
function isfuncdef(expr::Expr)
    expr.head == :function ||
    expr.head == :block && !isempty(expr.args) && isfuncdef(expr.args[1])
end

isfuncdef(expr) = false

macro cps(expr)
    cpsexpr = cps(expr, :identity)
    debug() && println(striplineinfo(cpsexpr))
    esc(cpsexpr)
end

# Distinguish thunks from functions returns by the transformed code.
immutable Thunk
    f::Function
    function Thunk(f)
        new(f)
    end
end

function trampoline(t::Thunk)
    thunkcount = 0
    while isa(t, Thunk)
        t = t.f()
        if debug()
            #println(t.f.code)
            thunkcount += 1
        end
    end
    debug() && println("trampolined thunks: ", thunkcount)
    t
end

# Convenience function used to kick-off trampolining in enum() &
# pmcmc().
trampoline(f::Function) = trampoline(Thunk(f))

# CPS transform based on:
# http://matt.might.net/articles/cps-conversion/
# Essentials of Programming Languages (3rd Ed.) [Friedman & Wand]

# | CPS                   | EoPL              | Might's Scheme CPS |                 |
# | tc                    | cps-of-exp        | T-c                |                 |
# | tk(::Function, expr)  | cps-of-exp/ctx    | T-k                | Ex. 6.30, p225. |
# | tk(::Function, exprs) | cps-of-exps       | T*-k               |                 |
# | m                     | cps-of-simple-exp | M                  |                 |

function tc(expr, cont)
    if isatomic(expr)
        :(Stochy.Thunk(() -> $cont($(m(expr)))))
    elseif expr.head == :block && length(expr.args) == 1
        tc(expr.args[1], cont)
    elseif expr.head == :block && length(expr.args) > 1
        tk(expr.args[1]) do _
            # Build block this way rather than with quoting so line
            # number expressions aren't added, preventing recursion
            # from ever bottoming out.
            tc(Expr(:block, expr.args[2:end]...), cont)
        end
    elseif expr.head == :local
        @assert expr.args[1].head == symbol("=")
        @assert length(expr.args) == 1
        val = expr.args[1].args[2]
        tk(e -> :($cont($e)), val)
    # elseif expr.head == :let
        # @assert all([arg.head == symbol("=") for arg in expr.args[2:end]])
        # body = expr.args[1]
        # vars = [arg.args[1]::Symbol for arg in expr.args[2:end]]
        # exprs = [arg.args[2] for arg in expr.args[2:end]]
        # tk(exprs) do e...
        #     :(let $(map((var, val) -> :($(var)=$(val)), vars, e)...); $(tc(body, cont)) end)
        # end
    elseif expr.head == :if
        @assert length(expr.args) == 3
        tk(expr.args[1]) do test
            :($test ? $(tc(expr.args[2], cont)) : $(tc(expr.args[3], cont)))
        end
    elseif expr.head == :function
        f = expr.args[1].args[1]
        args = expr.args[1].args[2:end]
        body = expr.args[2]
        v = gensym()
        :($cont(function $f($v, $(args...)); $(tc(body, v)); end))
    elseif expr.head == :||
        # Convert to if statement.
        @assert length(expr.args) == 2
        lhs, rhs = expr.args
        tc(:($lhs ? true : ($rhs ? true : false)), cont)
    elseif expr.head == :&&
        @assert length(expr.args) == 2
        lhs, rhs = expr.args
        tc(:($lhs ? ($rhs ? true : false) : false), cont)
    elseif isprimop(expr)
        f = expr.args[1]
        tk(expr.args[2:end]) do e...
            :($cont($f($(e...))))
        end
    elseif expr.head == :call
        tk(expr.args) do f, e...
            :($f($cont, $(e...)))
        end
    elseif expr.head == :.
        @assert isa(expr.args[2], QuoteNode) || isa(expr.args[2], Expr)
        @assert length(expr.args) == 2
        # The dot syntax turns in to two different Exprs depending on
        # whether you use the macro or :(<expr>). Handle both.
        if isa(expr.args[2], Expr)
            @assert expr.args[2].head == :quote
            @assert length(expr.args[2].args) == 1
            sym = expr.args[2].args[1]
        else
            @assert isa(expr.args[2].value, Symbol)
            sym = expr.args[2].value
        end
        tk(expr.args[1]) do e
            :($cont($(e).$(sym)))
        end
    elseif expr.head == :comparison && expr.args[2] in compops
        @assert length(expr.args) == 3
        tk([expr.args[1], expr.args[3]]) do x, y
            compexpr = Expr(expr.head, x, expr.args[2], y)
            :($cont($compexpr))
        end
    elseif expr.head == :vcat
        tk(expr.args) do e...
            :($cont([$(e...)]))
        end
    elseif expr.head == :ref
        @assert length(expr.args) == 2
        tk(expr.args) do x, y
            :($cont($x[$y]))
        end
    else
        error("$(expr.head) expression not recognized")
    end
end

function tk(k::Function, expr::Union(Expr, Atom))
    if isatomic(expr)
        k(m(expr))
    elseif expr.head == :local && isexpr(expr.args[1].args[2], :call)
        # This is an optimization.
        # When the result of a function call is bound to a local
        # variable, bind the variable by using its name as the
        # continuation's argument.
        # e.g. local a = f(x); <expr> => f(x, (a) -> <expr>)
        #
        # This optimization is contained wholly within this elseif.
        v = expr.args[1].args[1]
        callexpr = expr.args[1].args[2]
        tc(callexpr, :($v -> $(k(v))))
    elseif expr.head == :local
        v = expr.args[1].args[1]
        tc(expr, :($v -> $(k(v))))
    elseif expr.head == :...
        v = gensym()
        vsplat = Expr(:..., v)
        tc(expr.args[1], :($v -> $(k(vsplat))))
    else
        v = gensym()
        tc(expr, :($v -> Stochy.Thunk(() -> $(k(v)))))
    end
end

function tk(k::Function, exprs::Array)
    if isempty(exprs)
        k(()...)
    else
        tk(first(exprs)) do e
            tk(exprs[2:end]) do es...
                k(e, es...)
            end
        end
    end
end

function m(expr::Expr)
    @assert isatomic(expr)
    if expr.head == :->
        k = symbol("##k00")
        args = procargs(expr.args[1])
        # This type annotation is required to fix a failing test
        # caused by this issue:
        # https://github.com/JuliaLang/julia/issues/9770
        body = expr.args[2]::Expr
        :(($k, $(args...)) -> $(tc(body, k)))
    elseif expr.head == :quote
        expr
    else
        error("unreachable")
    end
end

m(expr::Atom) = expr

procargs(exp::Expr) = exp.args
procargs(exp) = [exp]

isatomic(exp::Atom) = true
isatomic(exp::Expr) = exp.head == :-> || exp.head == :quote
isatomic(exp) = false

isprimop(exp::Expr) = exp.head == :call && exp.args[1] in primitives

isline(exp::Expr) = exp.head == :line
isline(exp::LineNumberNode) = true
isline(exp) = false

not(pred) = x -> !(pred(x))
split(pos::Int, a::Array) = (a[1:pos-1],a[pos:end])

striplineinfo(exp::Expr) = Expr(exp.head, map(striplineinfo, filter(not(isline), exp.args))...)
striplineinfo(exp) = exp

function desugar(expr::Expr)
    if expr.head == :block
        # Expand local expressions of n variables into n local
        # expressions each introducing a single variable.
        newargs = map(expandlocals, expr.args) |> flatten
        Expr(:block, map(desugar, newargs)...)
    elseif expr.head == :local && length(expr.args) > 1
        # Add a top-level local expression to a block and recurse to
        # expand if required.
        desugar(Expr(:block, expr))
    elseif expr.head == :call && expr.args[1] == :~
        # Re-write calls to ~.
        Expr(:call, :sample, map(desugar, expr.args[2:end])...)
    else
        Expr(expr.head, map(desugar, expr.args)...)
    end
end

desugar(expr::Atom) = expr
desugar(expr::LineNumberNode) = expr

flatten(a) = vcat(a...)

function expandlocals(expr::Expr)
    if expr.head == :local
        map(arg->:(local $(arg.args[1])=$(arg.args[2])), expr.args)
    else
        expr
    end
end

expandlocals(expr::Atom) = expr
expandlocals(expr::LineNumberNode) = expr


text:
include("data.jl")
include("technical_indicators.jl")

function crossover(asset)
   """ creates crossover signal when the 26 moving average cross over the 52 moving average, returns either True or False """
   if asset.data[end-1].MA26 < asset.data[end-1].MA52 and asset.data[end].MA26 > asset.data[end].MA52
      return true
   else
      return false
   end
end

function crossunder(asset)
   """ creates a crossunder signal when the 26 moving average cross under the 52 moving average, return either True or False """
   if asset.data[end-1].MA26 > asset.data[end-1].MA52 and asset.data[end].MA26 < asset.data[end].MA52
      return True
   else
      return False
   end
end

function overRSI(asset)
   """ Returns True if RSI indicator is over a value of 60 """
   if asset.data[end].RSI > 60
      return true
   else
      return false
   end
end

function underRSI(asset)
   """ Returns True if RSI indicators is under a value of 40 """
   if asset.data[end].RSI < 40
      return true
   else
      return false
   end

function breaksupport(asset, support)
   """ Returns True when the last candlestick closes below a support """
   if asset.data[end].close < support
      return true
   else
      return false
   end
end

function breakresistance(asset)
   """ Returns True when the last candlestick closes above a resistance """
   if asset.data[end].close > resistance
      return true
   else
      return false
   end
end

max_stars_repo_name: SRelieved/JuliaBacktesting

max_issues_repo_name: MartinGenesse/JuliaBacktesting

max_forks_repo_name: MartinGenesse/JuliaBacktesting


text:
# Some extra LAPACK functions not defined in Base.

const liblapack = Base.liblapack_name

import Base.LinAlg: BlasFloat, BlasChar, BlasInt, blas_int, LAPACKException,
    DimensionMismatch, SingularException, PosDefException, chkstride1, chksquare

#Generic LAPACK error handlers
macro assertargsok() #Handle only negative info codes - use only if positive info code is useful!
    :(info[1]<0 && throw(ArgumentError("invalid argument #$(-info[1]) to LAPACK call")))
end
macro lapackerror() #Handle all nonzero info codes
    :(info[1]>0 ? throw(LAPACKException(info[1])) : @assertargsok )
end

macro assertnonsingular()
    :(info[1]>0 && throw(SingularException(info[1])))
end
macro assertposdef()
    :(info[1]>0 && throw(PosDefException(info[1])))
end

#Check that upper/lower (for special matrices) is correctly specified
macro chkuplo()
    :((uplo=='U' || uplo=='L') || throw(ArgumentError("""invalid uplo = $uplo

Valid choices are 'U' (upper) or 'L' (lower).""")))
end

num_packed(np) = div(-1 + isqrt(1+8*np), 2)

# (PO) positive-definite symmetric matrices,
# Cholesky decomposition, solvers (direct and factored) and inverse.
for (ppsv, pptrf, pptri, pptrs, elty, rtyp) in
    ((:dppsv_,:dpptrf_,:dpptri_,:dpptrs_,:Float64,:Float64),
     (:sppsv_,:spptrf_,:spptri_,:spptrs_,:Float32,:Float32),
     (:zppsv_,:zpptrf_,:zpptri_,:zpptrs_,:Complex128,:Float64),
     (:cppsv_,:cpptrf_,:cpptri_,:cpptrs_,:Complex64,:Float32))
    @eval begin
        #       SUBROUTINE DPPSV( UPLO, N, NRHS, AP, B, LDB, INFO )
        #       .. Scalar Arguments ..
        #       CHARACTER          UPLO
        #       INTEGER            INFO, LDB, N, NRHS
        #       ..
        #       .. Array Arguments ..
        #       DOUBLE PRECISION   AP( * ), B( LDB, * )
        function ppsv!(uplo::BlasChar, A::StridedVector{$elty}, B::StridedVecOrMat{$elty})
            chkstride1(A, B)
            n = size(B,1)
            @chkuplo
            if size(B,1) != n throw(DimensionMismatch("ppsv!")) end
            info    = Array(BlasInt, 1)
            ccall(($(string(ppsv)),liblapack), Void,
                  (Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
                   Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}),
                  &uplo, &n, &size(B,2), A, B, &max(1,stride(B,2)), info)
            @assertargsok
            @assertposdef
            A, B
        end
        #       SUBROUTINE DPPTRF( UPLO, N, AP, INFO )
        #
        #       .. Scalar Arguments ..
        #       CHARACTER          UPLO
        #       INTEGER            INFO, N
        #       ..
        #       .. Array Arguments ..
        #       DOUBLE PRECISION   AP( * )
        function pptrf!(uplo::BlasChar, A::StridedVector{$elty})
            chkstride1(A)
            @chkuplo
            info = Array(BlasInt, 1)
            n = Array(BlasInt,1)
            n[1] = num_packed(length(A))

            ccall(($(string(pptrf)),liblapack), Void,
                  (Ptr{BlasChar}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}),
                  &uplo, n, A, info)
            @assertargsok
            #info[1]>0 means the leading minor of order info[i] is not positive definite
            #ordinarily, throw Exception here, but return error code here
            #this simplifies isposdef! and factorize!
            return A, info[1]
        end
        #       SUBROUTINE DPPTRI( UPLO, N, AP, INFO )
        #
        #       .. Scalar Arguments ..
        #       CHARACTER          UPLO
        #       INTEGER            INFO, N
        #       ..
        #       .. Array Arguments ..
        #       DOUBLE PRECISION   AP( * )

        function pptri!(uplo::BlasChar, A::StridedVector{$elty})
            @chkuplo
            info = Array(BlasInt, 1)
            n = Array(BlasInt, 1)
            n[1] = num_packed(length(A))
            ccall(($(string(pptri)),liblapack), Void,
                  (Ptr{BlasChar}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}),
                  &uplo, n, A, info)
            @assertargsok
            @assertnonsingular
            A
        end
        #       SUBROUTINE DPPTRS( UPLO, N, NRHS, AP, B, LDB, INFO )
        #
        #       .. Scalar Arguments ..
        #       CHARACTER          UPLO
        #       INTEGER            INFO, LDB, N, NRHS
        #       ..
        #       .. Array Arguments ..
        #       DOUBLE PRECISION   AP( * ), B( LDB, * )
        #     SUBROUTINE DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
        #     .. Scalar Arguments ..
        #      CHARACTER          UPLO
        #      INTEGER            INFO, LDA, LDB, N, NRHS
        #     .. Array Arguments ..
        #      DOUBLE PRECISION   A( LDA, * ), B( LDB, * )
        function pptrs!(uplo::BlasChar, A::StridedVector{$elty}, B::StridedVecOrMat{$elty})
            chkstride1(A, B)
            n = Array(BlasInt,1)
            n[1] = num_packed(length(A))

            @chkuplo
            nrhs = size(B,2)
            if size(B,1) != n[1] throw(DimensionMismatch("left and right hand sides do not fit")) end
            ldb = max(1,stride(B,2))
            info = Array(BlasInt, 1)
            ccall(($(string(pptrs)),liblapack), Void,
                  (Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
                    Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}),
                   &uplo, n, &nrhs, A,
                   B, &ldb, info)
            @lapackerror
            return B
        end
    end
end

for (pbtrf, pbtrs, elty) in
    ((:dpbtrf_,:dpbtrs_,:Float64),
     (:spbtrf_,:spbtrs_,:Float32),
     (:zpbtrf_,:zpbtrs_,:Complex128),
     (:cpbtrf_,:cpbtrs_,:Complex64))
    @eval begin
        #       SUBROUTINE DPBTRF( UPLO, N, KD, AB, LDAB, INFO )
        #
        #       .. Scalar Arguments ..
        #       CHARACTER          UPLO
        #       INTEGER            INFO, KD, LDAB, N
        #       ..
        #       .. Array Arguments ..
        #       DOUBLE PRECISION   AB( LDAB, * )
        function pbtrf!(uplo::BlasChar, kd::Integer, AB::StridedMatrix{$elty})
            @chkuplo
            chkstride1(AB)
            info = Array(BlasInt, 1)
            n    = size(AB, 2)
            ccall(($(string(pbtrf)),liblapack), Void,
                  (Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt},
                   Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}),
                  &uplo, &n, &kd, AB, &max(1,stride(AB,2)), info)
            @lapackerror
            AB
        end

        #       SUBROUTINE DPBTRS( UPLO, N, KD, NRHS, AB, LDAB, B, LDB, INFO )
        #
        #       .. Scalar Arguments ..
        #       CHARACTER          UPLO
        #       INTEGER            INFO, KD, LDAB, LDB, N, NRHS
        #       ..
        #       .. Array Arguments ..
        #       DOUBLE PRECISION   AB( LDAB, * ), B( LDB, * )
        function pbtrs!(uplo::BlasChar, kd::Integer,
                        AB::StridedMatrix{$elty}, B::StridedVecOrMat{$elty})
            @chkuplo
            chkstride1(AB, B)
            info = Array(BlasInt, 1)
            n    = size(AB,2)
            if n != size(B,1) throw(DimensionMismatch("pbtrs!")) end
            ccall(($(string(pbtrs)),liblapack), Void,
                  (Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt},
                   Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty},   Ptr{BlasInt},
                   Ptr{BlasInt}),
                  &uplo, &n, &kd, &size(B,2), AB, &max(1,stride(AB,2)),
                  B, &max(1,stride(B,2)), info)
            @lapackerror
            B
        end
    end
end


text:
using Base.Test
using FastGaussQuadrature

println("Chebyshev tests")
include("test_gausschebyshev.jl")
println("Legendre tests")
include("test_gausslegendre.jl")
println("Jacobi tests")
include("test_gaussjacobi.jl")
println("Radau tests")
include("test_gaussradau.jl")
println("Lobatto tests")
include("test_gausslobatto.jl")
println("Bessel roots tests")
include("test_besselroots.jl")
println("Laguerre tests")
include("test_gausslaguerre.jl")
println("Hermite tests")
include("test_gausshermite.jl")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, BinaryBuilderBase

name = "ripgrep"
version = v"13.0.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/BurntSushi/ripgrep/archive/$(version).tar.gz",
                  "0fb17aaf285b3eee8ddab17b833af1e190d73de317ff9648751ab0660d763ed2"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/ripgrep-*/
cargo build --release --features 'pcre2'
mkdir -p "${bindir}"
cp "target/${rust_target}/release/rg${exeext}" "${bindir}/."
install_license COPYING LICENSE-MIT UNLICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)
# Rust toolchain for i686 Windows is unusable
filter!(p -> !Sys.iswindows(p) || arch(p) != "i686", platforms)
# All Musl-based ARM platforms fail with this error:
#
#     = note: /workspace/srcdir/ripgrep-13.0.0/target/aarch64-unknown-linux-musl/release/deps/libpcre2_sys-1c7f8ac14084107f.rlib(pcre2_jit_compile.o): In function `sljit_generate_code':
#             /opt/x86_64-linux-musl/registry/src/github.com-1ecc6299db9ec823/pcre2-sys-0.2.5/pcre2/src/sljit/sljitExecAllocator.c:178: undefined reference to `__clear_cache'
#             collect2: error: ld returned 1 exit status
filter!(p -> libc(p) != "musl" || proc_family(p) != "arm", platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("rg", :rg),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat="1.6", compilers=[:c, :rust], lock_microarchitecture=false)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
    DarioHitsFile(file::Union{IOStream, AbstractString}; batch_size::Integer=10)

represents a `.root.hits` file given by MaGe's Dario output scheme. A `DarioHitsFile` can be
iterated or read to yield events, which contain data on energy depositions, i.e. hits.

A `DarioHitsFile` is also a Tables.jl compatible row table of events. `batch_size`
determines the number of events grouped into partitions when `Tables.partitions` is used.
To change the default of 10, see `DARIO_HITS_BATCH_SIZE`.
"""
struct DarioHitsFile
    stream::IOBuffer
    batch_size::Int

    function DarioHitsFile(stream::IOStream; batch_size=DARIO_HITS_BATCH_SIZE[])
        new(IOBuffer(mmap(stream)), batch_size)
    end
end

function DarioHitsFile(path::AbstractString; batch_size=DARIO_HITS_BATCH_SIZE[])
    if occursin(r".root.hits$", path)
        return DarioHitsFile(open(path); batch_size=batch_size)
    else
        throw(ArgumentError("$path is not a .root.hits file"))
    end
end

const DARIO_HITS_BATCH_SIZE = Ref(10)

const DarioHitsEventTuple = NamedTuple{
    (:eventnum, :primcount, :pos, :E, :time, :particleID, :trkID, :trkparentID, :volumeID),
    Tuple{
        Int32, Int32, Vector{SVector{3, Float32}}, Vector{Float32}, Vector{Float32},
        Vector{Int32}, Vector{Int32}, Vector{Int32}, Vector{String}
    }
}

# TODO: Add documentation for this ^

function Base.read(f::DarioHitsFile)
    eventnum  = Parsers.parse(Int32, f.stream)
    hitcount  = Parsers.parse(Int32, f.stream)
    primcount = Parsers.parse(Int32, f.stream)

    # skip newline
    skip(f.stream, 1)

    pos         = Vector{SVector{3, Float32}}(undef, hitcount)
    E           = Vector{            Float32}(undef, hitcount)
    time        = Vector{            Float32}(undef, hitcount)
    particleID  = Vector{              Int32}(undef, hitcount)
    trkID       = Vector{              Int32}(undef, hitcount)
    trkparentID = Vector{              Int32}(undef, hitcount)
    volumeID    = Vector{             String}(undef, hitcount)

    @inbounds for i in 1:hitcount
        pos[i] = SVector{3, Float32}(
            Parsers.parse(Float32, f.stream),
            Parsers.parse(Float32, f.stream),
            Parsers.parse(Float32, f.stream)
        )
        E[i]           = Parsers.parse(Float32, f.stream)
        time[i]        = Parsers.parse(Float32, f.stream)
        particleID[i]  = Parsers.parse(Int32, f.stream)
        trkID[i]       = Parsers.parse(Int32, f.stream)
        trkparentID[i] = Parsers.parse(Int32, f.stream)
        volumeID[i]    = String(readuntil(f.stream, UInt8('\n')))
    end

    return DarioHitsEventTuple((
        eventnum, primcount, pos, E, time,
        particleID, trkID, trkparentID, volumeID
    ))
end

Base.eof(f::DarioHitsFile) = eof(f.stream)

function Base.iterate(f::DarioHitsFile, state = nothing)
    eof(f) && return nothing
    return read(f), nothing
end

Base.IteratorSize(::Type{DarioHitsFile}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{DarioHitsFile}) = Base.HasEltype()
Base.eltype(::Type{DarioHitsFile}) = DarioHitsEventTuple

Tables.isrowtable(::Type{DarioHitsFile}) = true
Tables.schema(f::DarioHitsFile) = Tables.Schema(eltype(f))

Tables.partitions(f::DarioHitsFile) = Iterators.partition(f, f.batch_size)


text:
struct EvolutionaryOptimization <: QuerySynthesisOptimizer
    method::Symbol
    params::Dict{Symbol, Any}
    function EvolutionaryOptimization(method::Symbol; params...)
        new(method, Dict(params))
    end
end

function query_synthesis_optimize(f::Function, optimizer::EvolutionaryOptimization, data::Array{T, 2}, labels::Vector{Symbol})::Array{T, 2} where T <: Real
    (x_opt, fitness, iterations) = eval(optimizer.method)(x -> -first(f(reshape(x, length(x), 1))), size(data, 1); optimizer.params...)
    return reshape(x_opt, size(data, 1), 1)
end


text:
using AQFED, Test
import AQFED.PLN: EtoreGobetPLNPricer, LeFlochLehmanPLNPricer, priceEuropean
import AQFED.TermStructure: CapitalizedDividend, Dividend, futureValue
import AQFED.Basket: DeelstraBasketPricer, DeelstraLBBasketPricer

@testset "EtoreGobetSingle" begin
    spot = 100.0
    vol = 0.3
    discountRate = 0.0
    divAmount = 7.0
    tte = 1.0
    ttd = 0.9 * tte
    ttp = tte
    isCall = true
    discountFactor = exp(-discountRate * ttp)
    dividends = Vector{CapitalizedDividend}(undef, 1)
    #amount::T        exDate::Float64        payDate::Float64        isProportional::bool        isKnown::bool
    dividends[1] =
        CapitalizedDividend(Dividend(divAmount, ttd, ttd, false, false), exp((tte - ttd) * discountRate))
    println("S Model Price Error")
    eg2 = EtoreGobetPLNPricer(2)
    eg3 = EtoreGobetPLNPricer(3)

    ll2 = LeFlochLehmanPLNPricer(2)
    ll3 = LeFlochLehmanPLNPricer(3)
    d = DeelstraBasketPricer(1, 3)
    dlb = DeelstraLBBasketPricer(1, 3)
    refHHL = [43.24845580, 9.07480013, 1.06252880, 0.10473887]
    refEg2 = [43.24846889, 9.07479039, 1.06253441, 0.10474035]
    refEg3 = [43.24845582, 9.07480026, 1.06252875, 0.10473885]
    for (i, strike) in enumerate(LinRange(50.0, 200.0, 4))
        price = priceEuropean(
            eg2,
            isCall,
            strike,
            spot * exp(tte * discountRate),
            vol * vol * tte,
            tte,
            discountFactor,
            dividends,
        )
        println(strike, " EG-2 ", price, " ", price - refHHL[i])
        @test isapprox(refEg2[i], price, atol = 1e-8)
        price = priceEuropean(
            eg3,
            isCall,
            strike,
            spot * exp(tte * discountRate),
            vol * vol * tte,
            tte,
            discountFactor,
            dividends,
        )
        println(strike, " EG-3 ", price, " ", price - refHHL[i])
        @test isapprox(refEg3[i], price, atol = 1e-8)
        price = priceEuropean(
            ll2,
            isCall,
            strike,
            spot * exp(tte * discountRate),
            vol * vol * tte,
            tte,
            discountFactor,
            dividends,
        )
        println(strike, " LL-2 ", price, " ", price - refHHL[i])
        price = priceEuropean(
            ll3,
            isCall,
            strike,
            spot * exp(tte * discountRate),
            vol * vol * tte,
            tte,
            discountFactor,
            dividends,
        )
        println(strike, " LL-3 ", price, " ", price - refHHL[i])
        price = AQFED.Basket.priceEuropean(
            d,
            isCall,
            strike,
            spot * exp(tte * discountRate),
            vol * vol * tte,
            tte,
            discountFactor,
            dividends,
        )
        println(strike, " Deelstra ", price, " ", price - refHHL[i])
        price = AQFED.Basket.priceEuropean(
            dlb,
            isCall,
            strike,
            spot * exp(tte * discountRate),
            vol * vol * tte,
            tte,
            discountFactor,
            dividends,
        )
        println(strike, " Deelstra-LB ", price, " ", price - refHHL[i])

    end

end

@testset "MultipleVellekoop" begin
    spot = 100.0
    σ = 0.25
    discountRate = 0.06
    tte = 7.0
    isCall = true
    ttd = 0.9
    dividends = [CapitalizedDividend(Dividend(6.0, ttd, ttd, false, false), exp((tte - ttd) * discountRate)),
    CapitalizedDividend(Dividend(6.5, ttd + 1, ttd + 1, false, false), exp((tte - ttd - 1) * discountRate)),
    CapitalizedDividend(Dividend(7.0, ttd + 2, ttd + 2, false, false), exp((tte - ttd - 2) * discountRate)),
    CapitalizedDividend(Dividend(7.5, ttd + 3, ttd + 3, false, false), exp((tte - ttd - 3) * discountRate)),
    CapitalizedDividend(Dividend(8.0, ttd + 4, ttd + 4, false, false), exp((tte - ttd - 4) * discountRate)),
    CapitalizedDividend(Dividend(8.0, ttd + 5, ttd + 5, false, false), exp((tte - ttd - 5) * discountRate)),
    CapitalizedDividend(Dividend(8.0, ttd + 6, ttd + 6, false, false), exp((tte - ttd - 6) * discountRate))]
    rawForward = spot * exp(discountRate * tte)
    df = exp(-discountRate * tte)
    ll3 = LeFlochLehmanPLNPricer(3)
    d = DeelstraBasketPricer(3, 3)
    dlb = DeelstraLBBasketPricer(3, 3)
    refPrices = [
        34.19664044,
        30.49396147,
        27.21393776,
        24.31297892,
        21.74885502,
        19.48228598,
        17.47762235,
        15.70300783,
        14.13025908,
        12.73460326,
        11.49435473,
    ]
    f = spot * exp(discountRate * tte) - sum(futureValue(cd) for cd in dividends)

    for (i, strike) in enumerate(LinRange(50.0, 150.0, 11))
        refVol = Black.impliedVolatility(isCall, refPrices[i], f, strike, tte, df)
        price = priceEuropean(ll3, isCall, strike, rawForward, σ^2 * tte, tte, df, dividends)
        vol = Black.impliedVolatility(isCall, price, f, strike, tte, df)
        println(strike, " LL-3 ", price, " ", price / refPrices[i] - 1, " ", vol - refVol)
        price = AQFED.Basket.priceEuropean(d, isCall, strike, rawForward, σ^2 * tte, tte, df, dividends)
        vol = Black.impliedVolatility(isCall, price, f, strike, tte, df)
        println(strike, " Deelstra ", price, " ", price / refPrices[i] - 1, " ", vol - refVol)
        price = AQFED.Basket.priceEuropean(dlb, isCall, strike, rawForward, σ^2 * tte, tte, df, dividends)
        vol = Black.impliedVolatility(isCall, price, f, strike, tte, df)
        println(strike, " Deelstra-LB ", price, " ", price / refPrices[i] - 1, " ", vol - refVol)

    end
end

@testset "MultipleGocsei" begin
    spot = 100.0
    σ = 0.25
    discountRate = 0.03
    divAmount = 2.0
    tte = 10.0
    isCall = true
    ttd0 = 1.0 / 365
    dividends = Vector{CapitalizedDividend}(undef, 20)
    for i = 1:20
        ttd = ttd0 + (i - 1) / 2
        dividends[i] =
            CapitalizedDividend(Dividend(divAmount, ttd, ttd, false, false), exp((tte - ttd) * discountRate))
    end
    rawForward = spot * exp(discountRate * tte)
    df = exp(-discountRate * tte)
    ll3 = LeFlochLehmanPLNPricer(3)
    d = DeelstraBasketPricer(1, 3) #1,3 minimizes relative error in price, 3,3 minimizes volatility error
    dlb = DeelstraLBBasketPricer(1, 3)
    refPrices = [37.25772649, 28.59516740, 22.28192756, 17.61042951, 14.09731856, 11.41453970, 9.33695198]
    f = spot * exp(discountRate * tte) - sum(futureValue(cd) for cd in dividends)

    for (i, strike) in enumerate(LinRange(50.0, 200.0, 7))
        refVol = Black.impliedVolatility(isCall, refPrices[i], f, strike, tte, df)
        price = priceEuropean(ll3, isCall, strike, rawForward, σ^2 * tte, tte, df, dividends)
        vol = Black.impliedVolatility(isCall, price, f, strike, tte, df)
        println(strike, " LL-3 ", price, " ", price / refPrices[i] - 1, " ", vol - refVol)
        price = AQFED.Basket.priceEuropean(d, isCall, strike, rawForward, σ^2 * tte, tte, df, dividends)
        vol = Black.impliedVolatility(isCall, price, f, strike, tte, df)
        println(strike, " Deelstra ", price, " ", price / refPrices[i] - 1, " ", vol - refVol)
        price = AQFED.Basket.priceEuropean(dlb, isCall, strike, rawForward, σ^2 * tte, tte, df, dividends)
        vol = Black.impliedVolatility(isCall, price, f, strike, tte, df)
        println(strike, " Deelstra-LB ", price, " ", price / refPrices[i] - 1, " ", vol - refVol)

    end
end


text:
# This example mirrors http://www.dabeaz.com/ply/example.html

using PyLexYacc

module rules
	using PyLexYacc

	tokens = (
		"NAME","NUMBER",
		"PLUS","MINUS","TIMES","DIVIDE","EQUALS",
		"LPAREN","RPAREN"
		)

	# Tokens

	t_PLUS    = "\\+"
	t_MINUS   = "-"
	t_TIMES   = "\\*"
	t_DIVIDE  = "/"
	t_EQUALS  = "="
	t_LPAREN  = "\\("
	t_RPAREN  = "\\)"
	t_NAME    = "[a-zA-Z_][a-zA-Z0-9_]*"

	t_NUMBER = rule("\\d+") do t, lexer
			t["value"] = int(t["value"])
			return t
		end

	t_ignore  = " \t"

	t_newline = rule("[\\r\\n]+") do t, lexer
			lexer["lineno"] += count((x)->(x == '\n'),t["value"])
			return nothing
		end

	t_error = rule() do t, lexer
			@printf("Illegal character '%s'\n", t["value"][1])
			skip(lexer, 1)
			return t
		end

	precedence = (
			("left", "PLUS", "MINUS"),
			("left", "TIMES", "DIVIDE"),
			("right", "UMINUS")
		)

	vars = Dict()  # symbol table

	p_statement_assign = parserule("statement : NAME EQUALS expression") do t
			vars[t[2]] = t[4]
		end

	p_statement_expr = parserule("statement : expression") do t
			println(t[2])
		end

	p_expression_binop = parserule(
			"""expression : expression PLUS expression
			              | expression MINUS expression
			              | expression TIMES expression
			              | expression DIVIDE expression""") do t
			if t[3] == "+"
				t[1] = t[2] + t[4]
			elseif t[3] == "-"
				t[1] = t[2] - t[4]
			elseif t[3] == "*"
				t[1] = t[2] * t[4]
			elseif t[3] == "/"
				t[1] = t[2] / t[4]
			end
		end

	p_expression_uminus = parserule("expression : MINUS expression %prec UMINUS") do t
			t[1] = -t[3]
		end

	p_expression_group = parserule("expression : LPAREN expression RPAREN") do t
			t[1] = t[3]
		end

	p_expression_number = parserule("expression : NUMBER") do t
			t[1] = t[2]
		end

	p_expression_name = parserule("expression : NAME") do t
			try
				t[1] = vars[t[2]]
			catch
				@printf("Undefined name '%s'\n", t[2])
				t[1] = 0
			end
		end

	p_error = rule() do t, lexer  # the error rule actually takes a LexToken so make it a rule()
			@printf("Syntax error at '%s'\n", t["value"][1])
		end
end

l = lexer(rules)
p = parser(rules, "statement")  # specifiy start rule here or through start variable in module

# replicates the input function using in calc.py
function input(str::String)
	print(str)
	return readline(STDIN)
end

println("\\q to quit")
line = input("calc > ")
while chomp(line) != "\\q"
	parse(p, l, line)
	line = input("calc > ")
end


text:
using Tsys
using Test

# These J1939-6342 (aka 1934-638) flux models are from different epochs, so
# returned flux is expected to be similar, but not identical.
model_1e6 = [-30.7667, 26.4908, -7.0977, 0.605334] # From MeerKAT, ν1=1e6
model_1e9 = [1.186E+00, 1.191E-01, -1.191E+00]     # From ATCA, ν1=1e9

@testset "Tsys.jl" begin
    @test model_flux(model_1e6) == 10^model_1e6[1]
    @test model_flux(model_1e9) == 10^model_1e9[1]
end


text:
using Temporal
const YAHOO_URL = "http://real-chart.finance.yahoo.com/table.csv"
const QUANDL_URL = "https://www.quandl.com/api/v3/datasets"

# haspkg(pkg::String)::Bool = isdir("$(Pkg.dir())/$pkg")

function matchcount(s::String, c::Char=',')
    x = 0
    @inbounds for i = 1:length(s)
        if s[i] == c
            x += 1
        end
    end
    return x
end

function tsread(file::String; dlm::Char=',', header::Bool=true, eol::Char='\n',
                indextype::Type=Date, format::String="yyyy-mm-dd")
    @assert indextype == Date || indextype == DateTime "Argument `indextype` must be either `Date` or `DateTime`."
    csv = Vector{String}(split(readstring(file), eol))
    if csv[end] == ""
        pop!(csv)  # remove final blank line
    end
    if header
        fields = Vector{String}(split(shift!(csv), dlm)[2:end])
        k = length(fields)
        n = length(csv)
    else
        k = matchcount(csv[1], dlm)
        n = length(csv)
        fields = autocolname(1:k)
    end
    # Fill data
    arr = zeros(Float64, (n,k))
    idx = fill("", n)::Vector{String}
    for i = 1:n
        s = Vector{String}(split(csv[i], dlm))
        idx[i] = shift!(s)
        s[s.==""] = "NaN"
        arr[i,:] = float(s)
    end
    return ts(arr, indextype(idx), fields)
end

function tswrite(x::TS, file::String; dlm::Char=',', header::Bool=true, eol::Char='\n')
    outfile = open(file, "w")
    if header
        write(outfile, "Index$(dlm)$(join(x.fields, dlm))$(eol)")
    end
    arr = x.values
    idx = x.index
    for i = 1:length(idx)
        write(outfile, "$(idx[i])$(dlm)$(join(arr[i,:],dlm))$(eol)")
    end
    close(outfile)
end

# ==============================================================================
# WEB INTERFACE ================================================================
# ==============================================================================

function dateconv(s::String)
    Dates.datetime2unix(Dates.DateTime(s))
end

function isdate(t::Vector{DateTime})
    h = Dates.hour(t)
    m = Dates.minute(t)
    s = Dates.second(t)
    ms = Dates.millisecond(t)
    return all(h.==h[1]) && all(m.==m[1]) && all(s.==s[1]) && all(ms.==ms[1])
end

function csvresp(resp; sort::String="des")
    @assert resp.status == 200 "Error in download request."
    rowdata = Vector{String}(split(readall(resp), '\n'))
    header = Vector{String}(split(shift!(rowdata), ','))
    pop!(rowdata)
    if sort == "des"
        reverse!(rowdata)
    end
    N = length(rowdata)
    k = length(header)
    data = zeros(Float64, (N,k-1))
    v = map(s -> Array{String}(split(s, ',')), rowdata)
    t = map(s -> Dates.DateTime(s[1]), v)
    isdate(t) ? t = Date(t) : nothing
    @inbounds for i = 1:N
        j = (v[i] .== "")
        v[i][find(j)] = "NaN"
        data[i,:] = float(v[i][2:k])
    end
    return (data, t, header)
end

# ==============================================================================
# QUANDL INTERFACE =============================================================
# ==============================================================================

function quandl_auth{T<:String}(key::T="")
    authfile = "$(Pkg.dir())/quandl-auth"
    if key == ""
        if isfile(authfile)
            key = open(readstring, authfile)
        end
    else
        f = open(authfile, "w")
        write(f, key)
        close(f)
    end
    return key
end

"""
`quandl(code::String;
        from::String="",
        thru::String="",
        freq::String="daily",
        calc::String="none",
        sort::String="asc",
        rows::Int=0,
        auth::String=quandl_auth())`

Quandl data download
"""
function quandl(code::String;
                from::String="",
                thru::String="",
                freq::String="daily",
                calc::String="none",
                sort::String="asc",
                rows::Int=0,
                auth::String=quandl_auth())
    # Check arguments =========================================================
    @assert from=="" || (from[5]=='-' && from[8]=='-') "Argument `from` has invlalid format."
    @assert thru=="" || (thru[5]=='-' && thru[8]=='-') "Argument `thru` has invlalid format."
    @assert freq in ["daily","weekly","monthyl","quarterly","annual"] "Invalid `freq` argument."
    @assert calc in ["none","diff","rdiff","cumul","normalize"] "Invalid `calc` argument."
    @assert sort  == "asc" || sort == "des" "Argument `sort` must be either \"asc\" or \"des\"."
    if rows != 0 && (from != "" || thru != "")
        error("Cannot specify `rows` and date range (`from` or `thru`).")
    end
    # Format URL ===============================================================
    if rows == 0
        fromstr = from == "" ? "" : "&start_date=$from"
        thrustr = thru == "" ? "" : "&end_date=$thru"
        url = "$QUANDL_URL/$code.csv?$(fromstr)$(thrustr)&order=$sort&collapse=$freq&transform=$calc&api_key=$auth"
    else
        url = "$QUANDL_URL/$code.csv?&rows=$rows&order=$sort&collapse=$freq&transform=$calc&api_key=$auth"
    end
    indata = csvresp(get(url), sort=sort)
    return ts(indata[1], indata[2], indata[3][2:end])
end

"""
quandl_meta(database::String, dataset::String)
Quandl dataset metadata downloaded into a Julia Dict
"""
function quandl_meta(database::String, dataset::String)
    resp = get("$QUANDL_URL/$database/$dataset/metadata.json")
    @assert resp.status == 200 "Error downloading metadata from Quandl."
    return JSON.parse(readall(resp))["dataset"]
end

"""
quandl_search(;db::String="", qry::String="", perpage::Int=1, pagenum::Int=1)
Search Quandl for data in a given database, `db`, or matching a given query, `qry`.
"""
function quandl_search(;db::String="", qry::String="", perpage::Int=1, pagenum::Int=1)
    @assert db!="" || qry!="" "Must enter a database or a search query."
    dbstr = db   == "" ? "" : "database_code=$db&"
    qrystr = qry  == "" ? "" : "query=$(replace(qry, ' ', '+'))&"
    resp = get("$QUANDL_URL.json?$(dbstr)$(qrystr)per_page=$perpage&page=$pagenum")
    @assert resp.status == 200 "Error retrieving search results from Quandl"
    return JSON.parse(readall(resp))
end

# ==============================================================================
# YAHOO INTERFACE ==============================================================
# ==============================================================================

function yahoo(symb::String;
               from::String="1900-01-01",
               thru::String=string(Dates.today()),
               freq::Char='d')
    @assert freq in ['d','w','m','v'] "Argument `freq` must be in ['d','w','m','v']"
    @assert from[5] == '-' && from[8] == '-' "Argument `from` has invalid date format."
    @assert thru[5] == '-' && thru[8] == '-' "Argument `thru` has invalid date format."
    m = Base.parse(Int, from[6:7]) - 1
    a = m<10 ? string("0",m) : string(m)
    b = from[9:10]
    c = from[1:4]
    m = Base.parse(Int, thru[6:7]) - 1
    d = m<10 ? string("0",m) : string(m)
    e = thru[9:10]
    f = thru[1:4]
    indata = csvresp(get("$YAHOO_URL?s=$symb&a=$a&b=$b&c=$c&d=$d&e=$e&f=$f&g=$freq&ignore=.csv"))
    return ts(indata[1], indata[2], indata[3][2:end])
end

function yahoo(syms::Vector{String};
               from::String="1900-01-01",
               thru::String=string(Dates.today()),
               freq::Char='d')
    out = Dict()
    for s = syms
        out[s] = yahoo(s, from=from, thru=thru, freq=freq)
    end
    return out
end


text:
struct UnixDateCovariate <: AbstractCovariate{Float64}
    name::String
    datetimevar::AbstractDateTimeVariate
end

Base.length(var::UnixDateCovariate) = length(var.datetimevar)

function UnixDateCovariate(datetimevar::AbstractDateTimeVariate)
    UnixDateCovariate(getname(datetimevar), datetimevar)
end

function slice(unixdtcov::UnixDateCovariate, fromobs::Integer, toobs::Integer, slicelength::Integer)
    basevar = unixdtcov.datetimevar
    f = (tick::Int64) -> tick == zero(Int64) ? NaN64 : Dates.datetime2unix(Dates.epochms2datetime(tick))
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slices = slice(basevar, fromobs, toobs, slicelength)
    mapslice(f, slices, slicelength, Float64)
end

function Base.map(unixdtcov::UnixDateCovariate, dataframe::AbstractDataFrame)
    UnixDateCovariate(unixdtcov.name, map(unixdtcov.datetimevar, dataframe))
end

function covariate(datetimevar::AbstractDateTimeVariate)
    UnixDateCovariate(getname(datetimevar), datetimevar)
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
using Documenter
using MultiDimEquations

push!(LOAD_PATH,"../src/")
makedocs(sitename="MultiDimEquations.jl Documentation",
         pages = [
            "Index" => "index.md",
         ],
         format = Documenter.HTML(prettyurls = false)
)
# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/sylvaticus/MultiDimEquations.jl.git",
    devbranch = "master"
)

max_forks_repo_name: sylvaticus/MultiDimEquations


text:
using Gaugefields
using MPI
using LinearAlgebra

npe_arg = [1 1 1 1]
npe_arg[1] = parse(Int,　ARGS[1]) # npe_x
npe_arg[2] = parse(Int,　ARGS[2]) # npe_y
npe_arg[3] = parse(Int,　ARGS[3]) # npe_z
npe_arg[4] = parse(Int,　ARGS[4]) # npe_t

function MDtest!(snet,U,Dim,mpi=false)
    p = initialize_TA_Gaugefields(U)
    Uold = similar(U)
    substitute_U!(Uold,U)
    MDsteps = 100
    temp1 = similar(U[1])
    temp2 = similar(U[1])
    comb = 6
    factor = 1/(comb*U[1].NV*U[1].NC)
    numaccepted = 0

    numtrj = 100
    for itrj = 1:numtrj
        @time accepted = MDstep!(snet,U,p,MDsteps,Dim,Uold)
        numaccepted += ifelse(accepted,1,0)

        plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
        println("$itrj plaq_t = $plaq_t")
        println("acceptance ratio ",numaccepted/itrj)
    end
end

function calc_action(snet,U,p)
    NC = U[1].NC
    Sg = -evaluate_GaugeAction(snet,U)/NC
    Sp = p*p/2
    S = Sp + Sg
    return real(S)
end

function MDstep!(snet,U,p,MDsteps,Dim,Uold)
    Δτ = 1/MDsteps
    gauss_distribution!(p)
    Sold = calc_action(snet,U,p)
    substitute_U!(Uold,U)

    for itrj=1:MDsteps
        U_update!(U,p,0.5,Δτ,Dim,snet)

        P_update!(U,p,1.0,Δτ,Dim,snet)

        U_update!(U,p,0.5,Δτ,Dim,snet)

    end
    Snew = calc_action(snet,U,p)
    println("Sold = $Sold, Snew = $Snew")
    println("Snew - Sold = $(Snew-Sold)")
    ratio = min(1,exp(Snew-Sold))
    if rand() > ratio
        substitute_U!(U,Uold)
        return false
    else
        return true
    end
end

function U_update!(U,p,ϵ,Δτ,Dim,snet)
    temps = get_temporary_gaugefields(snet)
    temp1 = temps[1]
    temp2 = temps[2]
    expU = temps[3]
    W = temps[4]

    for μ=1:Dim
        exptU!(expU,ϵ*Δτ,p[μ],[temp1,temp2])
        mul!(W,expU,U[μ])
        substitute_U!(U[μ],W)

    end
end

function P_update!(U,p,ϵ,Δτ,Dim,snet) # p -> p +factor*U*dSdUμ
    NC = U[1].NC
    temps = get_temporary_gaugefields(snet)
    dSdUμ = temps[end]
    factor =  -ϵ*Δτ/(NC)

    for μ=1:Dim
        calc_dSdUμ!(dSdUμ,snet,μ,U)
        mul!(temps[1],U[μ],dSdUμ) # U*dSdUμ
        Traceless_antihermitian_add!(p[μ],factor,temps[1])
    end
end

function test1()
    NX = 24
    NY = 24
    NZ = 24
    NT = 48
    Nwing = 1
    Dim = 4
    NC = 3

    mpi = true
    #mpi = false
    if mpi
        println("mpi=$mpi npe=$(npe_arg)")
    else
	println("mpi=$mpi")
    end
    if mpi
        PEs = npe_arg # (1,1,1,2)

        u1 = IdentityGauges(NC,Nwing,NX,NY,NZ,NT,mpi=true,PEs = PEs,mpiinit = false)
        U = Array{typeof(u1),1}(undef,Dim)
        U[1] = u1
        for μ=2:Dim
            U[μ] = IdentityGauges(NC,Nwing,NX,NY,NZ,NT,mpi=true,PEs = PEs,mpiinit = false)
        end
    else

        u1 = IdentityGauges(NC,Nwing,NX,NY,NZ,NT)
        U = Array{typeof(u1),1}(undef,Dim)
        U[1] = u1
        for μ=2:Dim
            U[μ] = IdentityGauges(NC,Nwing,NX,NY,NZ,NT)
        end
    end

    snet = GaugeAction(U)
    plaqloop = make_loops_fromname("plaquette")
    append!(plaqloop,plaqloop')
    β = 5.7/2
    push!(snet,β,plaqloop)

    #show(snet)

    MDtest!(snet,U,Dim,mpi)

end

test1()


text:
include("filtering.jl");
##
list = ["PreVehicle",
    "Altanserin",
    "SB242084",
    "Way_100135",
    "Methysergide",
    "Citalopram",
    "SB242084_opt",
    "Saline"]
s = filter(r->r.Treatment in list &&
    r.Trial < 61 &&
    r.Trial_Travel_to < 10 &&
    r.MouseID != "pc7",
    streaks)
s[s.Treatment .== "PreVehicle",:Treatment] .= "Control"
s[s.Treatment .== "Saline",:Treatment] = [o ? "Optogenetic" : "Control" for o in s[s.Treatment .== "Saline",:Stim]]
s[s.Treatment .== "SB242084_opt",:Phase] .=  "SB242084_opt"
s[s.Treatment .== "SB242084_opt",:Treatment] = [o ? "SB242084_opt" : "Control" for o in s[s.Treatment .== "SB242084_opt",:Stim]]
##Webers Law

release = filter(r -> r.Treatment in ["Citalopram","Optogenetic"],s)
WebersLaw(release,:AfterLast,:Treatment)
savefig(joinpath(figs_loc,"WebersLaw/ReleaseManipulation.pdf"))

selective = filter(r -> r.Treatment in ["SB242084","Altanserin"],s)
WebersLaw(selective,:AfterLast,:Treatment)
savefig(joinpath(figs_loc,"WebersLaw/SelectiveManipulation.pdf"))

both = filter(r -> r.Treatment in ["Citalopram","Optogenetic","SB242084","Altanserin"],s)
WebersLaw(both,:AfterLast,:Treatment)
savefig(joinpath(figs_loc,"WebersLaw/EffectiveManipulation.pdf"))

manipulations = filter(r -> r.Treatment != "Control",s)
WebersLaw(manipulations,:AfterLast,:Treatment)
savefig(joinpath(figs_loc,"WebersLaw/AllManipulation.pdf"))

WebersLaw(s,:AfterLast,:Treatment)
savefig(joinpath(figs_loc,"WebersLaw/WithControl.pdf"))
##
Drug_colors!(s)

WebersLaw(s,:AfterLast,:Treatment)

################ 2c opto comparison across days ##############
optopharma = filter(r ->
    r.Phase in ["Optogenetic", "SB242084_opt"], s)
gd = groupby(optopharma,[:MouseID,:Phase,:Stim])
df1 = combine([:Num_pokes,:AfterLast,:Trial_Travel_to] => (n,a,t) ->
    (Num_pokes = mean(n), AfterLast = mean(a), Travel = mean(t))
    ,gd)
df2 = select(df1, Not([:Travel,:Num_pokes]))
df3 = combine(groupby(df2,:Stim)) do dd
    subdf = unstack(dd,:Phase,:AfterLast)
    # current_drug = Symbol(subdf[1,:Phase])
    # rename!(subdf, current_drug => :Drug)
end
open_html_table(df3)
df4 = combine(groupby(df3,[:Stim])) do dd
    wilcoxon(dd,:SB242084_opt, :Optogenetic; f = x -> mean(skipmissing(x)))
end

plot_wilcoxon(df4)
############ Find the exact day before SB manipulation as control ###############
using Dates
stim_drug_days = [Date(2016,07,22),
    Date(2016,07,23),
    Date(2016,07,26),
    Date(2016,07,28),
    Date(2016,07,30),
    Date(2016,08,01)]
optopharma = filter(r ->
    r.Phase in ["Optogenetic", "SB242084_opt"], s)
open_html_table(optopharma)
union(filter(r -> r.Treatment == "SB242084_opt",optopharma).Day)
gd = groupby(optopharma,[:MouseID,:Phase,:Stim])
df1 = combine([:Num_pokes,:AfterLast,:Trial_Travel_to] => (n,a,t) ->
    (Num_pokes = mean(n), AfterLast = mean(a), Travel = mean(t))
    ,gd)
open_html_table(df1)
stim = filter(r -> r.Phase == "Optogenetic" && r.Stim, df1).AfterLast
nostim = filter(r -> r.Phase == "Optogenetic" && !r.Stim, df1).AfterLast
sb_stim = filter(r -> r.Phase == "SB242084_opt" && r.Stim, df1).AfterLast
sb_nostim = filter(r -> r.Phase == "SB242084_opt" && !r.Stim, df1).AfterLast

SignedRankTest(nostim,sb_nostim)


text:
module Transducers

export Transducer, Map, Filter, Cat, MapCat, Take, PartitionBy, Scan, Zip,
    Replace, TakeWhile, TakeNth, Drop, DropLast, DropWhile, Keep, Unique,
    Interpose, Dedupe, Partition, Iterated, Count, GroupBy, ReduceIf,
    TakeLast, FlagFirst, MapSplat, ScanEmit, Enumerate, NotA, OfType,
    transduce, eduction, setinput, Reduced, reduced, unreduced, ifunreduced,
    Completing, OnInit, CopyInit, right, reducingfunction, dreduce, dtransduce,
    withprogress, AdHocFoldable

# Deprecated:
export Distinct

using Base.Broadcast: Broadcasted

using ArgCheck
using BangBang: BangBang, Empty, append!!, empty!!, push!!, setindex!!
using Distributed: Distributed, @everywhere
using Logging: LogLevel, @logmsg
using Requires
using InitialValues: InitialValues, InitialValue, SpecificInitialValue, Init,
    hasinitialvalue

import Setfield
using Setfield: @lens, @set, set

# Dummy `ConstructionBase` module for supporting older `Setfield`:
module ConstructionBase
    using Setfield
    @static if isdefined(Setfield, :constructorof)
        const constructorof = Setfield.constructorof
    else
        const constructorof = Setfield.constructor_of
    end
end

@static if VERSION >= v"1.3-alpha"
    using Base.Threads: @spawn
else
    # Mock `@spawn` using `@async`:
    @eval const $(Symbol("@spawn")) = $(Symbol("@async"))
end

include("showutils.jl")
include("basics.jl")
include("core.jl")
include("library.jl")
include("simd.jl")
include("processes.jl")
include("reduce.jl")
include("dreduce.jl")
include("air.jl")
include("lister.jl")
include("show.jl")
include("comprehensions.jl")
include("progress.jl")
include("deprecated.jl")

include("interop/blockarrays.jl")
include("interop/lazyarrays.jl")

function __init__()
    @require BlockArrays="8e7c35d0-a365-5155-bbbb-fb81a777f24e" begin
        __foldl__(rf, acc, coll::BlockArrays.BlockArray) =
            _foldl_blockarray(rf, acc, coll)
    end
    @require LazyArrays="5078a376-72f3-5289-bfd5-ec5146d43c02" begin
        __foldl__(rf, acc, coll::LazyArrays.Hcat) =
            _foldl_lazy_hcat(rf, acc, coll)
        __foldl__(rf, acc, coll::LazyArrays.Vcat) =
            _foldl_lazy_vcat(rf, acc, coll)
    end
    @require OnlineStatsBase="925886fa-5bf2-5e8e-b522-a9147a512338" begin
        include("interop/onlinestats.jl")
    end
end

end # module


text:
using egads

# const libegads_include = "/home/mopg/code/EngSketchPad/lib/libegads.so"
# include the above in the precompilation step. use libdl.find_library to look for common places where this library file is located

egO   = Ref{ego}() # Pointer to a pointer essentially
egO_2 = Ref{ego}()
egO_3 = Ref{ego}()
egO_4 = Ref{ego}()

# the other way to do this is in the function
# function blabla(cid)
#            ptrobj = Ref{ego}()
#            call c function bullshit
#            return ptrobj[]
# end

# if you want to initialize a pointer to an object, you need to use for example
#       bla = ego( 0 ) # initialize to something

width = 5.0
rad   = 1.0

status = EG_open(egO)
if (status < 0) error("Can't open, failure code: %i", status) end
context = egO[] # now dereference the pointer to obtain the
@show context

# build cylinder
datacyl = [0.0, 0.0,  width/2.0,
           0.0, 0.0, -width/2.0,
           rad]
status = EG_makeSolidBody(context, CYLINDER, datacyl, egO_2)
ebody1 = egO_2[] # dereference the pointer to obtain the pointer to egObject
if (status != EGADS_SUCCESS) cleanup(status, context) end

# build end points of sphere
datasph = [0.0, 0.0, -width/2.0, rad]
status = EG_makeSolidBody(context, SPHERE, datasph, egO_3)
ebody2 = egO_3[] # dereference the pointer to obtain the pointer to egObject
if (status != EGADS_SUCCESS) cleanup(status, context) end

# boolean operations
status = EG_solidBoolean(ebody1, ebody2, FUSION, egO_4)
emodel = egO_4[]
if (status != EGADS_SUCCESS) cleanup(status, context) end

# save model
status = EG_saveModel(emodel, "test.egads")
if (status != EGADS_SUCCESS) cleanup(status, context) end
status = EG_saveModel(emodel, "test.step")
if (status != EGADS_SUCCESS) cleanup(status, context) end

# delete
status = EG_deleteObject(ebody1)
if (status != EGADS_SUCCESS) cleanup(status, context) end
status = EG_deleteObject(ebody2)
if (status != EGADS_SUCCESS) cleanup(status, context) end
status = EG_deleteObject(emodel)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# close everything
cleanup( context )


text:
using Static, Aqua
using Test

@testset "Static.jl" begin
    Aqua.test_all(Static)

    @testset "StaticInt" begin
        @test static(UInt(8)) === StaticInt(UInt(8)) === StaticInt{8}()
        @test iszero(StaticInt(0))
        @test !iszero(StaticInt(1))
        @test !isone(StaticInt(0))
        @test isone(StaticInt(1))
        @test @inferred(one(StaticInt(1))) === StaticInt(1)
        @test @inferred(zero(StaticInt(1))) === StaticInt(0)
        @test @inferred(one(StaticInt)) === StaticInt(1)
        @test @inferred(zero(StaticInt)) === StaticInt(0) === StaticInt(StaticInt(Val(0)))
        @test eltype(one(StaticInt)) <: Int

        x = StaticInt(1)
        @test @inferred(Bool(x)) isa Bool
        @test @inferred(BigInt(x)) isa BigInt
        @test @inferred(Integer(x)) === x
        @test @inferred(%(x, Integer)) === 1
        # test for ambiguities and correctness
        for i ∈ Any[StaticInt(0), StaticInt(1), StaticInt(2), 3]
            for j ∈ Any[StaticInt(0), StaticInt(1), StaticInt(2), 3]
                i === j === 3 && continue
                for f ∈ [+, -, *, ÷, %, <<, >>, >>>, &, |, ⊻, ==, ≤, ≥]
                    (iszero(j) && ((f === ÷) || (f === %))) && continue # integer division error
                    @test convert(Int, @inferred(f(i,j))) == f(convert(Int, i), convert(Int, j))
                end
            end
            i == 3 && break
            for f ∈ [+, -, *, /, ÷, %, ==, ≤, ≥]
                w = f(convert(Int, i), 1.4)
                x = f(1.4, convert(Int, i))
                @test convert(typeof(w), @inferred(f(i, 1.4))) === w
                @test convert(typeof(x), @inferred(f(1.4, i))) === x # if f is division and i === StaticInt(0), returns `NaN`; hence use of ==== in check.
                (((f === ÷) || (f === %)) && (i === StaticInt(0))) && continue
                y = f(convert(Int, i), 2 // 7)
                z = f(2 // 7, convert(Int, i))
                @test convert(typeof(y), @inferred(f(i, 2 // 7))) === y
                @test convert(typeof(z), @inferred(f(2 // 7, i))) === z
            end
        end

        @test UnitRange{Int16}(StaticInt(-9), 17) === Int16(-9):Int16(17)
        @test UnitRange{Int16}(-7, StaticInt(19)) === Int16(-7):Int16(19)
        @test UnitRange(-11, StaticInt(15)) === -11:15
        @test UnitRange(StaticInt(-11), 15) === -11:15
        @test UnitRange{Int}(StaticInt(-11), StaticInt(15)) === -11:15
        @test UnitRange(StaticInt(-11), StaticInt(15)) === -11:15
        @test float(StaticInt(8)) === static(8.0)

        # test specific promote rules to ensure we don't cause ambiguities
        SI = StaticInt{1}
        IR = typeof(1//1)
        PI = typeof(pi)
        @test @inferred(convert(SI, SI())) === SI()
        @test @inferred(promote_rule(SI, PI)) <: promote_type(Int, PI)
        @test @inferred(promote_rule(SI, IR)) <: promote_type(Int, IR)
        @test @inferred(promote_rule(SI, SI)) <: Int
        @test @inferred(promote_rule(Missing, SI)) <: promote_type(Missing, Int)
        @test @inferred(promote_rule(Nothing, SI)) <: promote_type(Nothing, Int)
        @test @inferred(promote_rule(SI, Missing)) <: promote_type(Int, Missing)
        @test @inferred(promote_rule(SI, Nothing)) <: promote_type(Int, Nothing)
        @test @inferred(promote_rule(Union{Missing,Int}, SI)) <: promote_type(Union{Missing,Int}, Int)
        @test @inferred(promote_rule(Union{Nothing,Int}, SI)) <: promote_type(Union{Nothing,Int}, Int)
        @test @inferred(promote_rule(Union{Nothing,Missing,Int}, SI)) <: Union{Nothing,Missing,Int}
        @test @inferred(promote_rule(Union{Nothing,Missing}, SI)) <: promote_type(Union{Nothing,Missing}, Int)
        @test @inferred(promote_rule(SI, Missing)) <: promote_type(Int, Missing)
        @test @inferred(promote_rule(Base.TwicePrecision{Int}, StaticInt{1})) <: Base.TwicePrecision{Int}

        @test static(Int8(-18)) === static(-18)
        @test static(0xef) === static(239)
        @test static(Int16(-18)) === static(-18)
        @test static(0xffef) === static(65519)
        if sizeof(Int) == 8
            @test static(Int32(-18)) === static(-18)
            @test static(0xffffffef) === static(4294967279)
        end
    end

    @testset "StaticBool" begin
        t = static(static(true))
        f = StaticBool(static(false))

        @test StaticBool{true}() === t
        @test StaticBool{false}() === f

        @test @inferred(StaticInt(t)) === StaticInt(1)
        @test @inferred(StaticInt(f)) === StaticInt(0)

        @test @inferred(~t) === f
        @test @inferred(~f) === t
        @test @inferred(!t) === f
        @test @inferred(!f) === t
        @test @inferred(+t) === StaticInt(1)
        @test @inferred(+f) === StaticInt(0)
        @test @inferred(-t) === StaticInt(-1)
        @test @inferred(-f) === StaticInt(0)
        @test @inferred(sign(t)) === t
        @test @inferred(abs(t)) === t
        @test @inferred(abs2(t)) === t
        @test @inferred(iszero(t)) === f
        @test @inferred(isone(t)) === t
        @test @inferred(iszero(f)) === t
        @test @inferred(isone(f)) === f

        @test @inferred(xor(true, f))
        @test @inferred(xor(f, true))
        @test @inferred(xor(f, f)) === f
        @test @inferred(xor(f, t)) === t
        @test @inferred(xor(t, f)) === t
        @test @inferred(xor(t, t)) === f

        @test @inferred(|(true, f))
        @test @inferred(|(true, t)) === t
        @test @inferred(|(f, true))
        @test @inferred(|(t, true)) === t
        @test @inferred(|(f, f)) === f
        @test @inferred(|(f, t)) === t
        @test @inferred(|(t, f)) === t
        @test @inferred(|(t, t)) === t

        @test @inferred(Base.:(&)(true, f)) === f
        @test @inferred(Base.:(&)(true, t))
        @test @inferred(Base.:(&)(f, true)) === f
        @test @inferred(Base.:(&)(t, true))
        @test @inferred(Base.:(&)(f, f)) === f
        @test @inferred(Base.:(&)(f, t)) === f
        @test @inferred(Base.:(&)(t, f)) === f
        @test @inferred(Base.:(&)(t, t)) === t

        @test @inferred(<(f, f)) === false
        @test @inferred(<(f, t)) === true
        @test @inferred(<(t, f)) === false
        @test @inferred(<(t, t)) === false

        @test @inferred(<=(f, f)) === true
        @test @inferred(<=(f, t)) === true
        @test @inferred(<=(t, f)) === false
        @test @inferred(<=(t, t)) === true

        @test @inferred(==(f, f)) === true
        @test @inferred(==(f, t)) === false
        @test @inferred(==(t, f)) === false
        @test @inferred(==(t, t)) === true

        @test @inferred(*(f, t)) === t & f
        @test @inferred(-(f, t)) === StaticInt(f) - StaticInt(t)
        @test @inferred(+(f, t)) === StaticInt(f) + StaticInt(t)

        @test @inferred(^(t, f)) == ^(true, false)
        @test @inferred(^(t, t)) == ^(true, true)

        @test @inferred(^(2, f)) == 1
        @test @inferred(^(2, t)) == 2

        @test @inferred(^(BigInt(2), f)) == 1
        @test @inferred(^(BigInt(2), t)) == 2

        @test @inferred(div(t, t)) === t
        @test_throws DivideError div(t, f)

        @test @inferred(rem(t, t)) === f
        @test_throws DivideError rem(t, f)
        @test @inferred(mod(t, t)) === f

        @test @inferred(all((t, t, t)))
        @test !@inferred(all((t, f, t)))
        @test !@inferred(all((f, f, f)))

        @test @inferred(any((t, t, t)))
        @test @inferred(any((t, f, t)))
        @test !@inferred(any((f, f, f)))

        x = StaticInt(1)
        y = StaticInt(0)
        z = StaticInt(-1)
        @test @inferred(Static.eq(y)(x)) === f
        @test @inferred(Static.eq(x, x)) === t

        @test @inferred(Static.ne(y)(x)) === t
        @test @inferred(Static.ne(x, x)) === f

        @test @inferred(Static.gt(y)(x)) === t
        @test @inferred(Static.gt(y, x)) === f

        @test @inferred(Static.ge(y)(x)) === t
        @test @inferred(Static.ge(y, x)) === f

        @test @inferred(Static.lt(x)(y)) === t
        @test @inferred(Static.lt(x, y)) === f

        @test @inferred(Static.le(x)(y)) === t
        @test @inferred(Static.le(x, y)) === f

        @test @inferred(Static.ifelse(t, x, y)) === x
        @test @inferred(Static.ifelse(f, x, y)) === y

        @test @inferred(promote_rule(True, True)) <: StaticBool
        @test @inferred(promote_rule(True, Bool)) <: Bool
        @test @inferred(promote_rule(Bool, True)) <: Bool
    end

    @testset "StaticSymbol" begin
        x = StaticSymbol(:x)
        y = StaticSymbol("y")
        z = StaticSymbol(1)
        @test y === StaticSymbol(:y)
        @test z === StaticSymbol(Symbol(1))
        @test @inferred(StaticSymbol(x)) === x
        @test @inferred(StaticSymbol(x, y)) === StaticSymbol(:x, :y)
        @test @inferred(StaticSymbol(x, y, z)) === static(:xy1)
    end

    @testset "static interface" begin
        v = Val((:a, 1, true))

        @test static(1) === StaticInt(1)
        @test static(true) === True()
        @test static(:a) === StaticSymbol{:a}()
        @test Symbol(static(:a)) === :a
        @test static((:a, 1, true)) === (static(:a), static(1), static(true))
        @test @inferred(static(v)) === (static(:a), static(1), static(true))
        @test_throws ErrorException static("a")

        @test @inferred(Static.is_static(v)) === True()
        @test @inferred(Static.is_static(typeof(v))) === True()
        @test @inferred(Static.is_static(typeof(static(true)))) === True()
        @test @inferred(Static.is_static(typeof(static(1)))) === True()
        @test @inferred(Static.is_static(typeof(static(:x)))) === True()
        @test @inferred(Static.is_static(typeof(1))) === False()
        @test @inferred(Static.is_static(typeof((static(:x),static(:x))))) === True()
        @test @inferred(Static.is_static(typeof((static(:x),:x)))) === False()
        @test @inferred(Static.is_static(typeof(static(1.0)))) === True()

        @test @inferred(Static.known(v)) === (:a, 1, true)
        @test @inferred(Static.known(typeof(v))) === (:a, 1, true)
        @test @inferred(Static.known(typeof(static(true))))
        @test @inferred(Static.known(typeof(static(false)))) === false
        @test @inferred(Static.known(typeof(static(1.0)))) === 1.0
        @test @inferred(Static.known(typeof(static(1)))) === 1
        @test @inferred(Static.known(typeof(static(:x)))) === :x
        @test @inferred(Static.known(typeof(1))) === missing
        @test @inferred(Static.known(typeof((static(:x),static(:x))))) === (:x, :x)
        @test @inferred(Static.known(typeof((static(:x),:x)))) === (:x, missing)

        @test @inferred(Static.dynamic((static(:a), static(1), true))) === (:a, 1, true)
    end

    @testset "tuple utilities" begin
        x = (static(1), static(2), static(3))
        y = (static(3), static(2), static(1))
        z = (static(1), static(2), static(3), static(4))
        T = Tuple{Int,Float64,String}
        @test @inferred(Static.invariant_permutation(x, x)) === True()
        @test @inferred(Static.invariant_permutation(x, y)) === False()
        @test @inferred(Static.invariant_permutation(x, z)) === False()

        @test @inferred(Static.permute(x, Val(x))) === x
        @test @inferred(Static.permute(x, (static(1), static(2)))) === (static(1), static(2))
        @test @inferred(Static.permute(x, x)) === x
        @test @inferred(Static.permute(x, y)) === y
        @test @inferred(Static.eachop(getindex, x)) === x

        get_tuple_add(::Type{T}, ::Type{X}, dim::StaticInt) where {T,X} = Tuple{Static._get_tuple(T, dim),X}
        @test @inferred(Static.eachop_tuple(Static._get_tuple, y, T)) === Tuple{String,Float64,Int}
        @test @inferred(Static.eachop_tuple(get_tuple_add, y, T, String)) === Tuple{Tuple{String,String},Tuple{Float64,String},Tuple{Int,String}}
        @test @inferred(Static.find_first_eq(static(1), y)) === static(3)
        # inferred is Union{Int,Nothing}
        @test Static.find_first_eq(1, map(Int, y)) === 3
    end

    @testset "NDIndex" begin
        x = NDIndex((1,2,3))
        y = NDIndex((1,static(2),3))
        z = NDIndex(static(3), static(3), static(3))

        @testset "constructors" begin
            @test static(CartesianIndex(3, 3, 3)) === z == Base.setindex(Base.setindex(x, 3, 1), 3, 2)
            @test @inferred(CartesianIndex(z)) === @inferred(Static.dynamic(z)) === CartesianIndex(3, 3, 3)
            @test @inferred(Static.known(z)) === (3, 3, 3)
            @test Tuple(@inferred(NDIndex{0}())) === ()
            @test @inferred(NDIndex{3}(1, static(2), 3)) === y
            @test @inferred(NDIndex{3}((1, static(2), 3))) === y
            @test @inferred(NDIndex{3}((1, static(2)))) === NDIndex(1, static(2), static(1))
            @test @inferred(NDIndex(x, y)) === NDIndex(1, 2, 3, 1, static(2), 3)
        end

        @test @inferred(Base.IteratorsMD.split(x, Val(2))) === (NDIndex(1, 2), NDIndex(3,))
        @test @inferred(length(x)) === 3
        @test @inferred(length(typeof(x))) === 3
        @test @inferred(y[2]) === 2
        @test @inferred(y[static(2)]) === static(2)

        @test @inferred(-y) === NDIndex((-1,-static(2),-3))
        @test @inferred(y + y) === NDIndex((2,static(4),6))
        @test @inferred(y - y) === NDIndex((0,static(0),0))
        @test @inferred(zero(x)) === NDIndex(static(0),static(0),static(0))
        @test @inferred(oneunit(x)) === NDIndex(static(1),static(1),static(1))
        @test @inferred(x * 3) === NDIndex((3,6,9))
        @test @inferred(3 * x) === NDIndex((3,6,9))

        @test @inferred(min(x, z)) === x
        @test @inferred(max(x, z)) === NDIndex(3, 3, 3)
        @test !@inferred(isless(y, x))
        @test @inferred(isless(x, z))
        @test @inferred(Static.lt(oneunit(z), z)) === static(true)

        A = rand(3,3,3);
        @test @inferred(to_indices(A, axes(A), (x,))) === (1, 2, 3)
        @test @inferred(to_indices(A, axes(A), ([y,y],))) == ([y, y],)
    end

    @test repr(static(float(1))) == "static($(float(1)))"
    @test repr(static(1)) == "static(1)"
    @test repr(static(:x)) == "static(:x)"
    @test repr(static(true)) == "static(true)"
    @test repr(static(CartesianIndex(1,1))) == "NDIndex(static(1), static(1))"
end

# for some reason this can't be inferred when in the "Static.jl" test set
known_length(x) = known_length(typeof(x))
known_length(::Type{T}) where {N,T<:Tuple{Vararg{Any,N}}} = N
known_length(::Type{T}) where {T} = missing
maybe_static_length(x) = Static.maybe_static(known_length, length, x)
x = ntuple(+, 10)
y = 1:10
@test @inferred(maybe_static_length(x)) === StaticInt(10)
@test @inferred(maybe_static_length(y)) === 10

include("float.jl")

#=
A = rand(3,4);

offset1(x::Base.OneTo) = static(1)
offset1(x::AbstractUnitRange) = first(x)

offsets(x) = Static._eachop(offset1, (axes(x),), Static.nstatic(Val(ndims(x))))

=#


text:
if basename(pwd()) == "aoc"
    cd("2016/12")
end

mutable struct AssembunnyVM
    a::Int
    b::Int
    c::Int
    d::Int
end
AssembunnyVM() = AssembunnyVM(0, 0, 0, 0)

function Base.getindex(vm::AssembunnyVM, i::AbstractString)
    if isdigit(last(i))
        parse(Int, i)
    else
        getfield(vm, Symbol(i))
    end
end
function Base.setindex!(vm::AssembunnyVM, val::AbstractString, i::AbstractString)
    setfield!(vm, Symbol(i), vm[val])
end
function Base.setindex!(vm::AssembunnyVM, val::Integer, i::AbstractString)
    setfield!(vm, Symbol(i), val)
end

function assembunny(filename::AbstractString; vm = AssembunnyVM())
    lines = split.(readlines(filename))
    pc = 1
    while firstindex(lines) <= pc <= lastindex(lines)
        line = lines[pc]
        try
            if line[1] == "cpy"
                vm[line[3]] = line[2]
            elseif line[1] == "inc"
                vm[line[2]] += 1
            elseif line[1] == "dec"
                vm[line[2]] -= 1
            elseif line[1] == "jnz"
                if !iszero(vm[line[2]])
                    pc += vm[line[3]]
                    continue
                end
            end
        catch err
            println("Crashing on line ", pc, ": ", join(line))
            rethrow(err)
        end
        pc += 1
    end
    vm
end

function part1()
    @assert assembunny("example.txt").a == 42
    assembunny("input.txt").a
end
part2() = assembunny("input.txt"; vm = AssembunnyVM(0, 0, 1, 0)).a

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
abstract type AbstractIsolate end

"""
    CellComposition structure

Structure that contains the cellular composition of isolates.
"""
@units @description mutable struct IsolateComposition <: AbstractIsolate
    Cell_volume::Vector{Float64}        |  m^3   | "Cell volume"
    Protein_volume::Vector{Float64}     |  m^3   | "Protein volume"
    Ribosome_volume::Vector{Float64}    |  m^3   | "Ribosome volume"
    Envelope_volume::Vector{Float64}    |  m^3   | "Cell envelope volume"
    CNP::Matrix{Float64}                | _      | "Cell stoichiometry"
    YEV::Matrix{Float64}                | mol/mol| "Yield of structure on CNP-reserve"
end

function IsolateComposition(Genome_size, Min_gen_time, Gram_stain)
    V_cell = genome_size_to_cell_volume(collect(skipmissing(Genome_size)))
    gmax   = log(2)./collect(skipmissing(Min_gen_time))
    V_protein = cell_volume_to_protein_volume(V_cell)
    V_r = cell_volume_to_ribosome_volume(V_cell, gmax)
    V_env = cell_volume_to_envelope_volume(V_cell, collect(skipmissing(Gram_stain)))
    CNP = zeros(3, size(V_cell,1))
    for i in 1:size(V_cell,1)
        CNP[:,i] = cell_volume_to_stoichiometry([V_cell[i]], [gmax[i]])
    end
    YEV = zeros(3, size(V_cell,1))
    for i in 1:size(V_cell,1)
        YEV[:,i] = cell_volume_to_yields([V_cell[i]], [gmax[i]])
    end
    return IsolateComposition(V_cell, V_protein, V_r, V_env, CNP, YEV)
end

function genome_size_to_cell_volume(L_DNA::Vector{Float64})
    # Kempes et al. (2016), Eq. 8
    v_N = 1.47e-27
    V_DNA = v_N*L_DNA
    D_0 = 3e-17
    β_D = 0.21
    V_c = (V_DNA/D_0).^(1/β_D) # m^3
end

function cell_volume_to_genome_size(V_c::Vector{Float64})
    # Kempes et al. (2016), Eq. 8
    D_0 = 3e-17
    β_D = 0.21
    V_DNA = D_0*V_c.^(β_D) # m^3
    v_N = 1.47e-27
    L_DNA = V_DNA./v_N
end

function cell_volume_to_protein_volume(V_c::Vector{Float64})
    # Kempes et al. (2016), Eq. 12
    P_0 = 3.42e-7
    β_p = 0.70
    V_p = P_0*V_c.^β_p # m^3
end

function cell_volume_to_ribosome_volume(V_c::Vector{Float64}, gmax::Vector{Float64})
    # Kempes et al. (2016), Eq. 13
    l_p = 975
    V_p = cell_volume_to_protein_volume(V_c)
    m_p = 5.81e-20
    d_p = 1.37e6
    v_p = m_p/d_p
    N_p = V_p/v_p
    ϕ = 6.20e-5
    η = 6.20e-5
    l_r = 4566
    r_r = 63
    μ = gmax./60^2
    # μ_0 = 4e7
    # β_B = 1.64
    # μ = μ_0*V_c.^(β_B-1)
    N_r = @. l_p*N_p*(ϕ/μ + 1)/(r_r/μ - l_r*(η/μ + 1))
    v_r = 3.04e-24
    V_r = N_r*v_r # m^3
end

function cell_volume_to_mRNA_volume(V_c::Vector{Float64}, gmax::Vector{Float64})
    # Kempes et al. (2016), Eq. 18
    V_r = cell_volume_to_ribosome_volume(V_c, gmax./60^2)
    v_r = 3.04e-24
    N_r = V_r/v_r
    v_mRNA = 1.43e-24
    n_mRNA = 1.08
    V_mRNA = v_mRNA*n_mRNA*N_r # m^3
end

function cell_volume_to_tRNA_volume(V_c::Vector{Float64}, gmax::Vector{Float64})
    # Kempes et al. (2016), Eq. 17
    V_r = cell_volume_to_ribosome_volume(V_c, gmax./60^2)
    v_r = 3.04e-24
    N_r = V_r/v_r
    v_tRNA = 3.10e-26
    n_tRNA = 9.3
    V_tRNA = v_tRNA*n_tRNA*N_r # m^3
end

function cell_volume_to_envelope_volume(V_c::Vector{Float64}, Gram_stain::Vector{String})
    # Kempes et al. (2016), Eq. 19
    r_c = (3*V_c/(4*π)).^(1/3)
    p_p = 0.149
    V_env = zeros(size(r_c,1))
    for i in 1:size(r_c, 1)
        if Gram_stain[i] == "(-)"
            r_env = 40.5e-9      # 27.5 if gram +
            V_env[i] = (V_c[i] - 4/3*π*(r_c[i] - r_env).^3)*(1-p_p) # m^3
        else
            r_env = 27.5e-9
            V_env[i] = (V_c[i] - 4/3*π*(r_c[i] - r_env).^3)*(1-p_p) # m^3
        end
    end
    return V_env
end

function cell_volume_to_dry_mass(V_c::Vector{Float64}, gmax::Vector{Float64}, Gram_stain::Vector{String})
    v_0 = 3.0e-17
    β_D = 0.21
    V_DNA = v_0*V_c.^β_D
    V_p = cell_volume_to_protein_volume(V_c)
    V_r = cell_volume_to_ribosome_volume(V_c, gmax)
    V_mRNA = cell_volume_to_mRNA_volume(V_c, gmax)
    V_tRNA = cell_volume_to_tRNA_volume(V_c, gmax)
    V_env = cell_volume_to_envelope_volume(V_c, Gram_stain)
    d_DNA = 2e6
    d_p = 1.37e6
    d_r = 1.79e6
    d_mem = 1.05e6
    d_RNA = 2e6
    m_d = d_DNA*V_DNA + d_p*V_p + d_r*V_r + d_mem*V_env + d_RNA*V_mRNA + d_RNA*V_tRNA  # g
end

function cell_volume_to_cell_number_density(V_c::Vector{Float64}, gmax::Vector{Float64}, Gram_stain::Vector{String})
    # assuming 47% of dry mass is C (BNID 100649)
    m_d = cell_volume_to_dry_mass(V_c, gmax, Gram_stain)
    M_C = 12.0111
    N_A = 6.022e23
    λ_B = M_C./(0.47*m_d*N_A)    # mol cells/mol C
end

function cell_volume_to_cellular_density(V_c::Vector{Float64}, gmax::Vector{Float64}, Gram_stain::Vector{String})
    # Kempes et al. (2016), Eq. S53
    v_0 = 3.0e-17
    β_D = 0.21
    V_DNA = v_0*V_c.^β_D
    V_p = cell_volume_to_protein_volume(V_c)
    V_r = cell_volume_to_ribosome_volume(V_c, gmax)
    V_mRNA = cell_volume_to_mRNA_volume(V_c, gmax)
    V_tRNA = cell_volume_to_tRNA_volume(V_c, gmax)
    V_env = cell_volume_to_envelope_volume(V_c, Gram_stain)
    V_w = V_c - (V_DNA + V_p + V_r + V_env + V_tRNA + V_mRNA)
    d_DNA = 2e6
    d_p = 1.37e6
    d_r = 1.79e6
    d_mem = 1.05e6
    d_RNA = 2e6
    d_w = 1e6
    d_c = @. (d_DNA*V_DNA + d_p*V_p + d_r*V_r + d_mem*V_env + d_RNA*V_mRNA + d_RNA*V_tRNA + d_w*V_w)/V_c    #g/m^3
end

function genome_size_to_rRNA_copy_number(L_DNA::Vector{Float64})
    # Roller et al. (2016) - Supp. Fig. 1
    rRNA  = @. 2^(L_DNA./1e6 - 2^2.8)/0.66 + 1.0
    return round.(rRNA)
end

function surface_area_volume_ratio(V_c::Vector{Float64})
    A_c = @. π^(1/3)*(6*V_c)^(2/3)
    SAV = @. A_c/V_c
end

function sav_from_growth_rate(gmax::Vector{Float64})
    SAV = @. -2.8*gmax + 9.3
end

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait