text:
include("shared.jl")

"""
Data type used in return value of the production kernel (V).
"""
struct NodeTypeAndXs
    node_type::Int
    xs::Vector{Float64}
end

"""
Data type for the return value of the aggregation kernel (W)
"""
struct CovFnAndMatrix
    node::Node
    cov_matrix::Matrix{Float64}
end

@gen (static) function production_kernel(xs::Vector{Float64})
    node_type = @trace(categorical(node_dist), :type)

    # to be passed to the corresponding aggregation kernel application
    v = NodeTypeAndXs(node_type, xs)

    # elements to be passed to each child production kernel application
    num_children = node_type_to_num_children[node_type]
    us = fill(xs, num_children)

    result = Production(v, us)
    return result
end

@gen function aggregation_kernel(node_type_and_xs::NodeTypeAndXs,
                                          child_outputs::Vector{CovFnAndMatrix})

    node_type::Int = node_type_and_xs.node_type
    xs::Vector{Float64} = node_type_and_xs.xs
    local node::Node
    local cov_matrix::Matrix{Float64}

    # constant kernel
    if node_type == CONSTANT
        @assert length(child_outputs) == 0
        param = @trace(uniform_continuous(0, 1), :param)
        node = Constant(param)
        cov_matrix = eval_cov_mat(node, xs)

    # linear kernel
    elseif node_type == LINEAR
        @assert length(child_outputs) == 0
        param = @trace(uniform_continuous(0, 1), :param)
        node = Linear(param)
        cov_matrix = eval_cov_mat(node, xs)

    # squared exponential kernel
    elseif node_type == SQUARED_EXP
        @assert length(child_outputs) == 0
        length_scale = 0.01 + @trace(uniform_continuous(0, 1), :length_scale)
        node = SquaredExponential(length_scale)
        cov_matrix = eval_cov_mat(node, xs)

    # periodic kernel
    elseif node_type == PERIODIC
        @assert length(child_outputs) == 0
        scale = 0.01 + @trace(uniform_continuous(0, 1), :scale)
        period = 0.01 + @trace(uniform_continuous(0, 1), :period)
        node = Periodic(scale, period)
        cov_matrix = eval_cov_mat(node, xs)

    # plus combinator
    elseif node_type == PLUS
        @assert length(child_outputs) == 2
        node = Plus(child_outputs[1].node, child_outputs[2].node)
        cov_matrix = child_outputs[1].cov_matrix .+ child_outputs[2].cov_matrix

    # times combinator
    elseif node_type == TIMES
        @assert length(child_outputs) == 2
        node = Times(child_outputs[1].node, child_outputs[2].node)
        cov_matrix = child_outputs[1].cov_matrix .* child_outputs[2].cov_matrix

    else
        error("unknown node type $node_type")
    end

    # to be passed to the parent aggregation kernel application
    w = CovFnAndMatrix(node, cov_matrix)

    return w
end

const covariance_prior_incremental = Recurse(
        production_kernel, aggregation_kernel,
        max_branch, # maximum number of children generated by production
        Vector{Float64}, # U (passed from production to its children)
        NodeTypeAndXs, # V (passed from production to aggregation)
        CovFnAndMatrix) # W (passed from aggregation to its parents, also Recurse's return type)

@gen (static) function model(xs::Vector{Float64})

    # sample covariance matrix
    cov_fn_and_matrix = @trace(covariance_prior_incremental(xs, 1), :tree)

    # sample diagonal noise
    noise = @trace(gamma(1, 1), :noise)

    # compute covariance matrix
    n = length(xs)
    cov_matrix = cov_fn_and_matrix.cov_matrix + Matrix((noise + 0.01) * LinearAlgebra.I, n, n)

    # sample from multivariate normal
    @trace(mvnormal(zeros(n), cov_matrix), :ys)

    node = cov_fn_and_matrix.node
    return node
end

Gen.load_generated_functions()

@gen function subtree_proposal_recursive(cur::Int)

    # base address for production kernel application 'cur'
    prod_addr = (cur, Val(:production))

    # base address for aggregation kernel application 'cur'
    agg_addr = (cur, Val(:aggregation))

    # sample node type
    node_type = @trace(categorical(node_dist), (cur, Val(:production)) => :type)

    # constant kernel
    if node_type == CONSTANT
        @trace(uniform_continuous(0, 1), agg_addr => :param)

    # linear kernel
    elseif node_type == LINEAR
        @trace(uniform_continuous(0, 1), agg_addr => :param)

    # squared exponential kernel
    elseif node_type == SQUARED_EXP
        @trace(uniform_continuous(0, 1), agg_addr => :length_scale)

    # periodic kernel
    elseif node_type == PERIODIC
        @trace(uniform_continuous(0, 1), agg_addr => :scale)
        @trace(uniform_continuous(0, 1), agg_addr => :period)

    # plus combinator
    elseif node_type == PLUS
        child1 = get_child(cur, 1, max_branch)
        child2 = get_child(cur, 2, max_branch)
        @trace(subtree_proposal_recursive(child1))
        @trace(subtree_proposal_recursive(child2))

    # times combinator
    elseif node_type == TIMES
        child1 = get_child(cur, 1, max_branch)
        child2 = get_child(cur, 2, max_branch)
        @trace(subtree_proposal_recursive(child1))
        @trace(subtree_proposal_recursive(child2))

    # unknown node type
    else
        error("Unknown node type: $node_type")
    end

    node_type
end

@gen function subtree_proposal(prev_trace)
    prev_subtree_node::Node = get_retval(prev_trace)
    (subtree_idx::Int, depth::Int) = @trace(pick_random_node(prev_subtree_node, 1, 0), :choose_subtree_root)
    new_subtree_node_type::Int = @trace(subtree_proposal_recursive(subtree_idx), :subtree)
    (subtree_idx, depth, new_subtree_node_type)
end

function subtree_involution(trace, fwd_choices::ChoiceMap, fwd_ret::Tuple, proposal_args::Tuple)
    (subtree_idx, subtree_depth, new_subtree_type) = fwd_ret
    model_args = get_args(trace)

    # populate constraints
    constraints = choicemap()
    set_submap!(constraints, :tree, get_submap(fwd_choices, :subtree))

    # obtain new trace and discard, which contains the previous subtree
    (new_trace, weight, _, discard) = update(trace, model_args, (NoChange(),), constraints)

    # populate backward assignment
    bwd_choices = choicemap()
    set_submap!(bwd_choices, :choose_subtree_root => :recurse_left,
        get_submap(fwd_choices, :choose_subtree_root => :recurse_left))
    for depth=0:subtree_depth-1
        bwd_choices[:choose_subtree_root => :done => depth] = false
    end
    if new_subtree_type in BINARY_OPS
        bwd_choices[:choose_subtree_root => :done => subtree_depth] = true
    end
    set_submap!(bwd_choices, :subtree, get_submap(discard, :tree))

    (new_trace, bwd_choices, weight)
end

@gen function noise_proposal(prev_trace)
    @trace(gamma(1, 1), :noise)
end

function inference(xs::Vector{Float64}, ys::Vector{Float64}, num_iters::Int, callback)

    # observed data
    constraints = choicemap()
    constraints[:ys] = ys

    # generate initial trace consistent with observed data
    (trace, _) = generate(model, (xs,), constraints)

    # do MCMC
    local covariance_fn::Node
    local noise::Float64
    for iter=1:num_iters

        covariance_fn = get_retval(trace)
        noise = trace[:noise]
        callback(covariance_fn, noise)

        # randomly pick a node to expand
        #root = pick_random_node(covariance_fn, 1, max_branch)

        # do MH move on the subtree
        (trace, _) = metropolis_hastings(trace, subtree_proposal, (), subtree_involution)

        # do MH move on the top-level white noise
        (trace, _) = metropolis_hastings(trace, noise_proposal, ())
    end
    return (covariance_fn, noise)
end

function experiment()

    # load and rescale the airline dataset
    (xs, ys) = get_airline_dataset()
    xs_train = xs[1:100]
    ys_train = ys[1:100]
    xs_test = xs[101:end]
    ys_test = ys[101:end]

    # set seed
    Random.seed!(0)

    # print MSE and predictive log likelihood on test data
    callback = (covariance_fn, noise) -> begin
        pred_ll = predictive_ll(covariance_fn, noise, xs_train, ys_train, xs_test, ys_test)
        mse = compute_mse(covariance_fn, noise, xs_train, ys_train, xs_test, ys_test)
        println("mse: $mse, predictive log likelihood: $pred_ll")
    end

    # do inference, time it
    @time (covariance_fn, noise) = inference(xs_train, ys_train, 1000, callback)

    # sample predictions
    pred_ys = predict_ys(covariance_fn, noise, xs_train, ys_train, xs_test)
end

experiment()


text:
@adjoint function evaluate(s::DotProduct, x::AbstractVector, y::AbstractVector)
  dot(x,y), Δ -> begin
    (nothing, Δ.*y, Δ.*x)
  end
end

ext: jl


text:
function f28825()
    eval(:(1+1))
end

max_stars_repo_name: musm/julia

max_issues_repo_name: musm/julia

max_forks_repo_name: musm/julia


text:
#Demo 1: Pulse generation
#-------------------------------------------------------------------------------

using CMDimCircuits
CMDimCircuits.@using_CData()

#Get a demo display:
include(CMDimCircuits.demoplotcfgscript)
#Normally use something like:
#CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
LBLAX_AMPV = "Amplitude (V)"
LBLAX_TIME = "Time (s)"
color1 = cons(:a, line = set(color=:red))
color2 = cons(:a, line = set(color=:blue))

#==Input data
===============================================================================#
tstart = 0
tstep = .01
tstop = 10

#==Computations
===============================================================================#
t = tstart:tstep:tstop
@show len = length(t)
tmax = maximum(t)

#Discrete-time data, ideal response:
tDT = DataTime(t) #Discrete-time time dataset
uideal = step(tDT, ndel=Index(len/4))
Πideal = pulse(tDT, ndel=Index(len/7), npw=Index(len/16))

#Continuous-time data, single-pole response:
tCT = DataF1(t) #Continuous-time time dataset
pulsewidth = tmax/16
u = step(tCT, Pole(3/pulsewidth,:rad), tdel=tmax/4)
Π = pulse(tCT, Pole(3/pulsewidth,:rad), tdel=tmax/7, tpw=pulsewidth)

#==Generate plot
===============================================================================#
plot = cons(:plot, nstrips=3,
	ystrip1 = set(axislabel=LBLAX_AMPV, striplabel="Step Response"),
	ystrip2 = set(axislabel=LBLAX_AMPV, striplabel="Pulse Response"),
	ystrip3 = set(axislabel=LBLAX_AMPV, striplabel="Sine wave"),
	xaxis = set(label=LBLAX_TIME),
)
push!(plot,
	cons(:wfrm, DataF1(uideal), color1, label="DT", strip=1),
	cons(:wfrm, u, color2, label="CT", strip=1),
	cons(:wfrm, DataF1(Πideal), color1, label="DT", strip=2),
	cons(:wfrm, Π, color2, label="CT", strip=2),
	cons(:wfrm, sin(tCT*(pi*10/tmax)), color2, strip=3),
)

pcoll = push!(cons(:plotcoll, title="Generating Simple Responses"), plot)
	pcoll.ncolumns = 1

#==Show results
===============================================================================#
PoleRad = SignalProcessing.PoleRad
f(x::Pole) = "Generic pole"
f(x::PoleRad{T}) where {T<:Number} = "Pole: $(value(x)) rad!"
phz=Pole(1,:Hz)
prad=Pole(1,:rad)
@show phz, prad
@show value(phz), value(prad)
@show value(:rad,phz), value(:Hz,prad)
@show f(phz), f(prad)

#==Display results in pcoll
===============================================================================#
EasyPlot.displaygui(pcoll)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll


text:
include("analyze_piece.jl")
include("simplify_set.jl")

function analyze_sequence(dyns1::AbstractVector, dyns2::AbstractVector,
                          times::AbstractVector, x0::AbstractVector, ε::Number,
                          m::Integer, ζ::Number;
                          log::Bool=false, return_sets::Bool=false)
    time = (0.0, times[1])
    P0 = □(x0, ε)
    for simplify_underapproximation in [true, false]
        if log
            println("* piece #1")
        end
        dyn1 = dyns1[1]
        dyn2 = dyns2[1]
        P1_ua, P1_oa, _, x1 =
            analyze_piece(dyn1, dyn2, x0, P0, ε, time, m, ζ, true; log=log)
        if return_sets
            P_uas = Vector{LazySet{Float64}}()
            push!(P_uas, P1_ua)
            P_oas = Vector{LazySet{Float64}}()
            push!(P_oas, P1_oa)
            xs = [x0, x1]
        end

        if isempty(P1_ua)
            answer = false
            is_definite_answer = false
            loop_break = true
            try_again = false
            n_loop = 0
            if log
                println("underapproximation is empty => give up")
            end
        else
            answer = true
            is_definite_answer = false
            loop_break = false
            try_again = false
            n_loop = length(dyns1)
        end

        for k in 2:n_loop
            if log
                println("* piece #$k")
            end

            dyn1 = dyns1[k]
            dyn2 = dyns2[k]
            time = (times[k-1], times[k])
            x0 = x1

            # overapproximation
            P0_oa = simplify_overapprox(P1_oa)
            P1_oa = analyze_piece(dyn1, dyn2, x0, P0_oa, ε, time, m, ζ, false;
                                  log=log)

            if isempty(P1_oa)
                if log
                    println("overapproximation is empty => no solution exists")
                end
                answer = false
                is_definite_answer = true
                loop_break = true
                break
            end

            # underapproximation
            P0_ua =
                simplify_underapprox(P1_ua; ignore=!simplify_underapproximation)
            P1_ua, _, _, x1 = analyze_piece(dyn1, dyn2, x0, P0_ua, ε, time, m,
                                            ζ, true; log=log)

            if isempty(P1_ua)
                answer = false
                is_definite_answer = false
                loop_break = true
                try_again = simplify_underapproximation
                if log
                    print("underapproximation is empty")
                    if try_again
                        println(" => trying again with more precise analysis")
                    else
                        println(" => give up")
                    end
                end
                break
            end

            if return_sets
                push!(P_uas, P1_ua)
                push!(P_oas, P1_oa)
                push!(xs, x1)
            end
        end

        if !loop_break
            answer = true
            is_definite_answer = !isempty(P1_ua)
        end

        if try_again
            if log
                println("trying analysis again without simplifying underapproximation")
            end
            continue
        end

        if return_sets
            return answer, is_definite_answer, P_uas, P_oas, xs
        end
        return answer, is_definite_answer
    end
end

max_stars_repo_name: HySynth/HySynth

max_issues_repo_name: HySynth/HySynth

max_forks_repo_name: HySynth/HySynth


text:
###############
# Prepare FMU #
###############

pathToFMU = joinpath(dirname(@__FILE__), "..", "model", ENV["EXPORTINGTOOL"], "SpringPendulum1D.fmu")

t_start = 0.0
t_stop = 1.0

myFMU = fmiLoad(pathToFMU)
@test myFMU.type == FMI.fmi2CoSimulation
@test fmi2IsCoSimulation(myFMU)
@test fmi2IsModelExchange(myFMU)
comp = fmiInstantiate!(myFMU; loggingOn=false)
@test comp != 0
# choose FMU or FMUComponent
fmuStruct = nothing
envFMUSTRUCT = ENV["FMUSTRUCT"]
if envFMUSTRUCT == "FMU"
    fmuStruct = myFMU
elseif envFMUSTRUCT == "FMUCOMPONENT"
    fmuStruct = comp
end
success = fmiSimulateCS(fmuStruct, t_start, t_stop)
@test success
sol, _ = fmiSimulateME(fmuStruct, t_start, t_stop)
@test sol.retcode == :Success
fmiUnload(myFMU)

myFMU = fmiLoad(pathToFMU; type=:ME)
@test myFMU.type == FMI.fmi2ModelExchange
comp = fmiInstantiate!(myFMU; loggingOn=false)
fmuStruct = nothing
envFMUSTRUCT = ENV["FMUSTRUCT"]
if envFMUSTRUCT == "FMU"
    fmuStruct = myFMU
elseif envFMUSTRUCT == "FMUCOMPONENT"
    fmuStruct = comp
end
sol, _ = fmiSimulate(fmuStruct, t_start, t_stop)
@test sol.retcode == :Success
fmiUnload(myFMU)

myFMU = fmiLoad(pathToFMU; type=:CS)
@test myFMU.type == FMI.fmi2CoSimulation
comp = fmiInstantiate!(myFMU; loggingOn=false)
fmuStruct = nothing
envFMUSTRUCT = ENV["FMUSTRUCT"]
if envFMUSTRUCT == "FMU"
    fmuStruct = myFMU
elseif envFMUSTRUCT == "FMUCOMPONENT"
    fmuStruct = comp
end
success = fmiSimulate(fmuStruct, t_start, t_stop)
@test success
fmiUnload(myFMU)


text:
using KernelEstimator
using Base.Test

# write your own tests here
include("testreg.jl")


text:
struct Queue{T, I, A <: AbstractAllocator{T, I}}
    alloc::A
    function Queue{T, I, A}(alloc::A) where {T, I, A <: AbstractAllocator{T, I}}
        new{T, I, A}(alloc)
    end
end

Base.push!(q::Queue{T, I, A}, t::T) where {T, I, A <: AbstractAllocator{T, I}} =
    allocateend!(q.alloc, t)

Base.isempty(q::Queue{T, I, A}) where {T, I, A <: AbstractAllocator{T, I}} =
    isempty(q.alloc)

function Base.pop!(q::Queue{T, I, A}) where {T, I, A <: AbstractAllocator{T, I}}
    t = q.alloc[q.alloc.first]
    deallocatebegin!(q.alloc)
    t
end

Base.empty!(q::Queue{T, I, A}) where {T, I, A <: AbstractAllocator{T, I}} =
    emptyend!(q.alloc)


text:
#!/usr/bin/julia
# github.com/atille/keeper.jl
# $ julia keeper.jl --url=https://... [--exts=jpg,gif...] [--dir=downloads] [--https]
function main(args)

	# We define globaly the variables that must be affected.
	url = nothing
	extensions = nothing
	output_dir = nothing
	ssl_forced = false
	protocol = "http"
	separator = "\""

	# For each argument submitted, we change the variable values.
	for arg in args
		if startswith(arg, "--url=")
			url = split(arg, "--url=")[2]
		end

		# Specific extensions, or default ones bellow.
		if startswith(arg, "--exts=")
			extensions = split(arg, "--exts=")[2]
			extensions = split(extensions, ",")
		end

		# Specify the output directory where will be downloaded the found files.
		if startswith(arg, "--dir=")
			output_dir = split(arg, "--dir=")[2]
		end

		# If declared, will change " as the separator to ' in HTML documents.
		is startswith(arg, "--separator=")
			separator = split(arg, "--separator=")[2]
		end

		# If specified, we force the HTTPS for download the files (if found).
		if isequal(arg, "--https")
			ssl_forced = true
			protocol = "https"
		end

		if isequal(arg, "--help")
			getHelp()
			quit()
		end
	end

	# If the script is launched without argument, we die it.
	if isequal(url, nothing)
		println("No argument submitted.")
		quit()
	end

	# If we didn't added --exts=, so we define few default extensions.
	if isequal(extensions, nothing)
		extensions = ["jpg", "png", "gif"]
	end

	# If we didn't added --dir=/path, we define the output directory as default.
	if isequal(output_dir, nothing)
		output_dir = "scrapped_content"
	end

	# If the specified folder isn't available, we create it, and all its subs.
	isdir(output_dir) || mkdir(output_dir)
	map(x -> mkdir(string(output_dir, "/", x)), extensions)

	# We download with wget/curl/fetch the HTML content.
	f = download(url)
	data = readstring(open(f, "r"))
	elements = split(data, separator)

	# For all the elements we got with the split function, we check if this is a URL that begins
	# with the right protocol (http/https) and ends with the wanted file extensions.
	for element in elements
		for extension in extensions
			if startswith(element, protocol) && endswith(element, extension)
				println(string("[$extension] : ", element))

				# We generate a new name, random.
				filename = randstring(6)

				# We launch the file download in background and asynchronous for more speed and
				# unprocedural process.
				@async run(`wget -O $output_dir/$extension/$filename.$extension $element -q`)
			end
		end
	end
end

function getHelp()

	println("Keeper | Content scrapper developed in Julia - POC")
	println(" - --url=\t define the url to scrap. Not optional.")
	println(" - --https\t files must be loaded with https. Optional.")
	println(" - --dir=\t directory of saved files.")
	println(" - --exts=\t separated with commas, target extensions. Default: jpg, gif, png.")
	println(" - --help\t display this page.")
end

@time main(ARGS)


text:
# Dates defines CompoundPeriod(t::Time)
Dates.CompoundPeriod(d::Date) =
    Year(d) + Month(d) + Day(d)

Dates.CompoundPeriod(dtm::DateTime) =
    Year(dtm) + Month(dtm) + Day(dtm) +
    Hour(dtm) + Minute(dtm) + Second(dtm) +
    Millisecond(dtm)

Dates.CompoundPeriod(nd::NanoDate) =
    Year(nd) + Month(nd) + Day(nd) +
    Hour(nd) + Minute(nd) + Second(nd) +
    Millisecond(nd) + Microsecond(nd) + Nanosecond(nd)

function Dates.Date(cperiod::CompoundPeriod)
    periods = canonicalize(cperiod).periods
    yr, mn, dy = 0, 1, 1
    qt, wk = 0, 0
    idx = findfirst(x->isa(x,Year), periods)
    if !isnothing(idx)
        yr = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Quarter), periods)
    if !isnothing(idx)
        qt = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Month), periods)
    if !isnothing(idx)
        mn = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Week), periods)
    if !isnothing(idx)
        wk = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Day), periods)
    if !isnothing(idx)
        dy = value(periods[idx])
    end
    mn = mn + 3*qt
    dy = dy + 7*wk
    Date(yr,mn,dy)
end

function Dates.Time(cperiod::CompoundPeriod)
    periods = canonicalize(cperiod).periods
    hr, mi, sc, ms, cs, ns = 0, 0, 0, 0, 0, 0
    idx = findfirst(x->isa(x,Hour), periods)
    if !isnothing(idx)
        hr = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Minute), periods)
    if !isnothing(idx)
        mi = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Second), periods)
    if !isnothing(idx)
        sc = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Millisecond), periods)
    if !isnothing(idx)
        ms = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Microsecond), periods)
    if !isnothing(idx)
        cs = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Nanosecond), periods)
    if !isnothing(idx)
        ns = value(periods[idx])
    end
    Time(hr,mi,sc,ms,cs,ns)
end

function Dates.DateTime(cperiod::CompoundPeriod)
    periods = canonicalize(cperiod).periods
    yr, mn, dy = 0, 1, 1
    qt, wk = 0, 0
    hr, mi, sc, ms = 0, 0, 0, 0
    idx = findfirst(x->isa(x,Year), periods)
    if !isnothing(idx)
        yr = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Quarter), periods)
    if !isnothing(idx)
        qt = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Month), periods)
    if !isnothing(idx)
        mn = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Week), periods)
    if !isnothing(idx)
        wk = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Day), periods)
    if !isnothing(idx)
        dy = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Hour), periods)
    if !isnothing(idx)
        hr = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Minute), periods)
    if !isnothing(idx)
        mi = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Second), periods)
    if !isnothing(idx)
        sc = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Millisecond), periods)
    if !isnothing(idx)
        ms = value(periods[idx])
    end
    mn = mn + 3*qt
    dy = dy + 7*wk
    DateTime(yr,mn,dy,hr,mi,sc,ms)
end

function NanoDate(cperiod::CompoundPeriod)
    periods = canonicalize(cperiod).periods
    yr, mn, dy = 0, 1, 1
    qt, wk = 0, 0
    hr, mi, sc, ms, cs, ns = 0, 0, 0, 0, 0, 0
    idx = findfirst(x->isa(x,Year), periods)
    if !isnothing(idx)
        yr = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Quarter), periods)
    if !isnothing(idx)
        qt = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Month), periods)
    if !isnothing(idx)
        mn = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Week), periods)
    if !isnothing(idx)
        wk = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Day), periods)
    if !isnothing(idx)
        dy = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Hour), periods)
    if !isnothing(idx)
        hr = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Minute), periods)
    if !isnothing(idx)
        mi = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Second), periods)
    if !isnothing(idx)
        sc = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Millisecond), periods)
    if !isnothing(idx)
        ms = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Microsecond), periods)
    if !isnothing(idx)
        cs = value(periods[idx])
    end
    idx = findfirst(x->isa(x,Nanosecond), periods)
    if !isnothing(idx)
        ns = value(periods[idx])
    end
    mn = mn + 3*qt
    dy = dy + 7*wk
    NanoDate(yr,mn,dy,hr,mi,sc,ms,cs,ns)
end

# length, iterate
Base.length(cperiod::CompoundPeriod) = length(cperiod.periods)

# reverse to process smaller period types before larger period types
Base.iterate(cperiod::CompoundPeriod) = iterate(reverse(cperiod.periods))
Base.iterate(cperiod::CompoundPeriod, state) = iterate(reverse(cperiod.periods), state)

function Base.:(+)(nd::NanoDate, cperiod::CompoundPeriod)
    for p in cperiod
        nd += p
    end
    nd
end
Base.:(+)(cperiod::CompoundPeriod, nd::NanoDate) = nd + cperiod

function Base.:(-)(nd::NanoDate, cperiod::CompoundPeriod)
    for p in cperiod
        nd -= p
    end
    nd
end

function Base.:(*)(a::Integer, b::CompoundPeriod)
    b = canonicalize(b)
    accum = similar(b.periods)
    for idx in eachindex(accum)
       accum[idx] = a * b.periods[idx]
    end
    CompoundPeriod(accum)
end

Base.:(*)(b::CompoundPeriod, a::Integer) = a * b

function Base.:(fld)(a::CompoundPeriod, b::Integer)
    b = canonicalize(a)
    accum = similar(a.periods)
    for idx in eachindex(accum)
       accum[idx] = fld(a.periods[idx], b)
    end
    CompoundPeriod(accum)
end

function Base.:(cld)(a::CompoundPeriod, b::Integer)
    b = canonicalize(a)
    accum = similar(a.periods)
    for idx in eachindex(accum)
       accum[idx] = cld(a.periods[idx], b)
    end
    CompoundPeriod(accum)
end

function Base.:(div)(a::CompoundPeriod, b::Integer)
    b = canonicalize(a)
    accum = similar(a.periods)
    for idx in eachindex(accum)
       accum[idx] = div(a.periods[idx], b)
    end
    CompoundPeriod(accum)
end


text:
#module TVVAR

# Code based on Forecasting with time-varying vector autoregressive models (2008), K. Triantafyllopoulos

# Constant
const FLOATVEC = Array{Fl, 1} where Fl <: AbstractFloat
const FLOATMAT = Array{Fl, 2} where Fl <: AbstractFloat
const FLOATARR = Array{Fl, 3} where Fl <: AbstractFloat

# Import
using LinearAlgebra: diag, kron, I, cholesky
using Distributions: Normal, MvNormal, MatrixNormal

# Include scripts
include("./src/types.jl")
include("./src/utils.jl")
include("./src/estimate.jl")
include("./src/simulate.jl")

# Exported types
#export TVVAR, KF

# Exported functions
#export estimate

#end


text:
export apply_back!, apply_back

const Rotor{D,T} = Union{RotationGate{D,T},PutBlock{D,<:Any,<:RotationGate{<:Any,T}}}

as_scalar(arr::AbstractArray{T,0}) where {T} = arr[]
as_scalar(arr) = arr

"""
    generator(rot::Rotor) -> AbstractBlock

Return the generator of rotation block.
"""
generator(rot::RotationGate) = rot.block
generator(rot::PutBlock) =
    PutBlock(rot.n, generator(rot |> content), rot |> occupied_locs)

function apply_back!(st, block::AbstractBlock, collector) #,AbstractContainer{<:PrimitiveBlock}
    out, outδ = st
    if nparameters(block) == 0
        adjblock = block'
        in = apply!(out, adjblock)
        inδ = apply!(outδ, adjblock)
        return (in, inδ)
    else
        throw(MethodError(apply_back, (st, block, collector)))
    end
end

function apply_back!(st, block::Subroutine{D}, collector) where {D}
    out, outδ = st
    focus!(out, block.locs)
    focus!(outδ, block.locs)
    apply_back!((out, outδ), content(block), collector)
    relax!(out, block.locs; to_nactive = nqudits(block))
    relax!(outδ, block.locs; to_nactive = nqudits(block))
    return (out, outδ)
end

function apply_back!(st, block::Rotor{D}, collector) where {D}
    out, outδ = st
    adjblock = block'
    backward_params!((out, outδ), block, collector)
    in = apply!(out, adjblock)
    inδ = apply!(outδ, adjblock)
    return (in, inδ)
end

function apply_back!(st, block::TimeEvolution, collector)
    out, outδ = st
    adjblock = block'

    out, outδ = st
    input = apply!(out, adjblock)
    for i=1:YaoArrayRegister._asint(nbatch(outδ))
        o = viewbatch(outδ, i)
        !all(x -> x ≈ 0.0im, o.state) && apply!(o, adjblock)
    end
    pushfirst!(collector, -sum(imag(input' * apply!(copy(outδ), block.H))))
    return (input, outδ)
end

function apply_back!(st, block::PutBlock, collector)
    out, outδ = st
    adjblock = block'
    in = apply!(out, adjblock)
    adjmat = outerprod(outδ, in)
    mat_back!(datatype(in), block, adjmat, collector)
    inδ = apply!(outδ, adjblock)
    return (in, inδ)
end

function apply_back!(st, block::KronBlock, collector)
    apply_back!(st, chain(nqudits(block), [put(loc => block[loc]) for loc in block.locs]), collector)
end

function apply_back!(st, block::ControlBlock, collector)
    out, outδ = st
    adjblock = block'
    in = apply!(out, adjblock)
    adjmat = outerprod(outδ, in)
    mat_back!(datatype(in), block, adjmat, collector)
    inδ = apply!(outδ, adjblock)
    return (in, inδ)
end

function apply_back!(st, block::Daggered, collector)
    out, outδ = st
    adjblock = block'
    in = apply!(out, adjblock)
    adjmat = outerprod(in, outδ)
    mat_back!(datatype(in), content(block), adjmat, collector)
    inδ = apply!(outδ, adjblock)
    return (in, inδ)
end

function apply_back!(st, block::Scale, collector)
    out, outδ = st
    # recoverting inputs
    if isunitary(block)
        in = apply!(out, block')
    elseif isunitary(content(block))
        in = apply!(regscale!(out, inv(factor(block))), content(block)')
    else
        @warn "rescaling a non-unitary block, call into inefficient back-propgation routine."
        in = apply!(out, matblock(inv(Matrix(block))))
    end
    # adjoint of matrix
    adjmat = outerprod(outδ, in)
    mat_back!(datatype(in), block, adjmat, collector)

    inδ = apply!(outδ, block')

    # differentiate scaling
    #if niparams(block) > 0
    #    g = dropdims(sum(conj.(state(in)) .* state(outδ), dims = (1, 2)), dims = (1, 2)) |> as_scalar
    #    pushfirst!(collector, g')
    #end
    return (in, inδ)
end

function apply_back!(st, circuit::ChainBlock, collector)
    for blk in Base.Iterators.reverse(subblocks(circuit))
        st = apply_back!(st, blk, collector)
    end
    return st
end

function apply_back!(st, circuit::Add, collector; in)
    out, outδ = st
    adjmat = outerprod(outδ, in)
    for blk in Base.Iterators.reverse(subblocks(circuit))
        mat_back!(datatype(in), blk, adjmat, collector)
    end
    inδ = apply!(outδ, circuit')
    (in, inδ)
end

function apply_back!(st, block::RepeatedBlock{D,C}, collector) where {D,C}
    if nparameters(content(block)) == 0
        return apply!.(st, Ref(block'))
    end
    res = Any[]
    st = apply_back!(st, chain(nqudits(block), [put(loc => content(block)) for loc in block.locs]), res)
    res = dropdims(sum(reshape(res, :, C), dims = 2), dims = 2) |> as_scalar
    prepend!(collector, res)
    return st
end

# TODO: concentrator, repeat, kron
apply_back!(st, block::Measure, collector) =
    throw(MethodError(apply_back!, (st, block, collector)))

function backward_params!(st, block::Rotor, collector)
    in, outδ = st
    Σ = generator(block)
    g =
        dropdims(sum(conj.(state(in |> Σ)) .* state(outδ), dims = (1, 2)), dims = (1, 2)) |>
        as_scalar
    pushfirst!(collector, -imag(g) / 2)
    in |> Σ
    nothing
end

"""
    apply_back(st::Tuple{<:AbstractArrayReg, <:AbstractArrayReg}, block::AbstractBlock; kwargs...) -> (out, outδ), paramsδ

The backward function of `apply!`. Returns a tuple of ((input register, gradient of input register), parameter gradients)
"""
function apply_back(st::Tuple{<:AbstractArrayReg,<:AbstractArrayReg}, block::AbstractBlock; kwargs...)
    col = []
    in, inδ = apply_back!(st, block, col; kwargs...)
    (in, inδ), col
end


text:
import Pkg
Pkg.activate(mktempdir())
Pkg.develop(Pkg.PackageSpec(path=ARGS[1]))
import Pluto
s = Pluto.ServerSession()
nb = Pluto.SessionActions.open_url(s, "https://raw.githubusercontent.com/fonsp/Pluto.jl/v0.12.16/sample/Basic.jl"; run_async=false)
@show [c.output_repr for c in nb.cells]
success = !any(c.errored for c in nb.cells)
exit(success ? 0 : 1)


text:
facts("Statistics") do

    fname = joinpath(dirname(@__FILE__), "..", "data", "test_Hz19.5-testing.bdf")

    s = read_SSR(fname)

    context("Global Field Power") do

        ~ = gfp(s.data)  # TODO: Check result
    end
end


text:
using DiffEqBase, DiffEqBase.InternalEuler
using IteratorInterfaceExtensions
using Test

f_1dlinear = (u,p,t) -> 1.01u
prob = ODEProblem(f_1dlinear,rand(),(0.0,1.0))
sol =solve(prob,InternalEuler.FwdEulerAlg();dt=1//2^(4))
array = collect(getiterator(sol))

@test length(array) == 17
@test collect(fieldnames(typeof(array[1]))) == [:timestamp, :value]

f_2dlinear = (du,u,p,t) -> du.=1.01u
prob = ODEProblem(f_2dlinear,rand(2),(0.0,1.0))
sol =solve(prob,InternalEuler.FwdEulerAlg();dt=1//2^(4))
array = collect(getiterator(sol))

@test length(array) == 17
@test collect(fieldnames(typeof(array[1]))) == [:timestamp, :value1, :value2]

struct LotkaVolterra
    syms::Vector{Symbol}
end
function (::LotkaVolterra)(du,u,p,t)
    du[1] = 1.5*u[1] - 1*u[1]*u[2]
    du[2] = -3*u[2] + 1*u[1]*u[2]
end
f_2dlinear_named = LotkaVolterra([:x,:y])
prob = ODEProblem(f_2dlinear_named,[1.0,1.0],(0.0,1.0))
sol =solve(prob,InternalEuler.FwdEulerAlg(),dt=0.1)
array = collect(getiterator(sol))

@test length(array) == 11
@test collect(fieldnames(typeof(array[1]))) == [:timestamp, :x, :y]


text:
# ESCB-2, "efficient sampling for combinatorial bandits".
# Based on https://arxiv.org/abs/1502.03475.

abstract type ESCB2OptimisationAlgorithm end

"""
  function optimise_linear_sqrtlinear(instance::CombinatorialInstance{T},
                                      algo::ESCB2OptimisationAlgorithm,
                                      linear::Dict{T, Float64},
                                      sqrtlinear::Dict{T, Float64},
                                      sqrtlinear_weight::Float64,
                                      bandit_round::Int;
                                      with_trace::Bool=false) where T

Optimise ESCB2's objective function, with a linear term (with coefficients in `linear`)
and a square root (with coefficients in `sqrtlinear`):

  max linear^T x + sqrtlinear_weight * sqrt(sqrtlinear^T x)
  s.t. x belongs to the combinatorial set defined by instance

Several implementations are provided, and can be chosen with the appropriate subtype of
`ESCB2OptimisationAlgorithm`.

Returns one (approximately) optimum solution. The exact guarantees depend on the
chosen algorithm. If `with_trace=true`, a second value is returned with comprehensive
details of the behaviour of the algorithm.
"""
function optimise_linear_sqrtlinear(::CombinatorialInstance{T},
                                    ::ESCB2OptimisationAlgorithm,
                                    ::Dict{T, Float64}, ::Dict{T, Float64}, ::Int,
                                    ::Float64) where T end

mutable struct ESCB2 <: Policy
  algo::ESCB2OptimisationAlgorithm
  confidence_bonus_scaling::Function
end

# Theoretical formula: log(n) + 4 m log(log(n)).
ESCB2(algo::ESCB2OptimisationAlgorithm) = ESCB2(algo, (state::State) -> log(state.round) / 2)
ESCB2(algo::ESCB2OptimisationAlgorithm, α::Float64) = ESCB2(algo, (state::State) -> α * log(state.round))

mutable struct ESCB2Details <: PolicyDetails
  n_iterations::Int
  best_objective::Float64
  solver_time::Float64

  function ESCB2Details()
    return new(0, 0.0, 0.0)
  end
end

function _escb2_confidence_values(policy::ESCB2, state::State{T}) where T
  weights = Dict(k => 1 / v for (k, v) in state.arm_counts)
  return weights, max(0, policy.confidence_bonus_scaling(state))
end

# TODO: shouldn't it go to a more generic function? Many algorithm work with indices, so maybe there is something to do... but only for black-box algorithms?
function escb2_index(policy::ESCB2, state::State{T}, solution::Vector{T}) where T
  w = state.arm_average_reward
  s2, f_n = _escb2_confidence_values(policy, state)
  return escb2_index(w, s2, f_n, solution)
end
function escb2_index(linear::Dict{T, Float64}, sqrtlinear::Dict{T, Float64}, sqrtlinear_weight::Float64, solution::Vector{T}) where T
  if length(solution) == 0
    return 0.0
  end

  avg_reward = sum(linear[arm] for arm in keys(linear) if arm in solution)
  confidence = sum(sqrtlinear[arm] for arm in keys(sqrtlinear) if arm in solution)
  return avg_reward + sqrt(sqrtlinear_weight * confidence)
end

function choose_action(instance::CombinatorialInstance{T}, policy::ESCB2, state::State{T}; with_trace::Bool=false) where T
  # For initialisation, if some arms have never been tried, force them to be tried.
  # Without this, the standard deviation term cannot exist.
  if any(v == 0 for v in values(state.arm_counts))
    weights = Dict(arm => iszero(state.arm_counts[arm]) ? 1.0 : 0.0 for arm in keys(state.arm_counts))

    t0 = time_ns()
    sol = solve_linear(instance, weights)
    t1 = time_ns()

    if with_trace
      run_details = ESCB2Details()
      run_details.n_iterations = 1
      run_details.solver_time = (t1 - t0) / 1_000_000_000
      return sol, run_details
    else
      return sol
    end
  end

  # If all arms have been tried, can compute the two parts of the objective function: the average reward and a standard deviation.
  w = state.arm_average_reward
  s2, f_n = _escb2_confidence_values(policy, state)

  # Solve the maximisation w^T x + \sqrt{f(n) s2^T x} with the required algorithm.
  return optimise_linear_sqrtlinear(instance, policy.algo, w, s2, f_n, state.round, with_trace=with_trace)
end


text:
include("MutableNamedTuples.jl")

import RecursiveArrayTools:
    recursive_unitless_eltype,
    recursive_bottom_eltype,
    recursive_unitless_bottom_eltype,
    recursivecopy

# Convenience closure for creating NamedTuple constructor
Struct(keys...) = (args...) -> (; zip(keys, args)...)
MutableStruct(keys...) = (args...) -> MutableNamedTuple((; zip(keys, args)...))

View{T} = SubArray{T,0,Array{T,1},Tuple{Int64},true}

struct NamedViewVector{T} <: AbstractVector{T}
    data::MutableNamedTuple
    vector::Array{View{T}, 1}
end
function NamedViewVector{T}(mnt::MutableNamedTuple) where {T}
    mnt = recursive_convert(T, mnt)
    return NamedViewVector{T}(mnt, attached_vect(mnt, T))
end
NamedViewVector{T}(nt::NamedTuple) where {T} = NamedViewVector{T}(MutableNamedTuple(nt))

# function attached_vect(mnt::MutableNamedTuple)
#     vector = []
#     for (key, val) in zip(keys(mnt), values(mnt))
#         if val[1] isa Number
#             push!(vector, view(val, 1))
#         elseif val[1] isa AbstractVector
#             for i in eachindex(val[1])
#                 push!(vector, view(val[1], i))
#             end
#         elseif val[1] isa MutableNamedTuple
#             push!(vector, attached_vect(val[1])...)
#         end
#     end
#     return vector
# end

# TODO: Make these all parameterized by ::Type{T}
# attached_vect(mnt, ::Type{T}) where {T} = recursive_convert(T, mnt) |> attached_vect
function attached_vect(mnt::MutableNamedTuple, ::Type{T}) where {T}
    mnt = recursive_convert(T, mnt)
    vector = View{T}[]
    for (key, val) in zip(keys(mnt), values(mnt))
        att_vect = attached_vect(val, T)
        push!(vector, att_vect...)
    end
    return vector
end
function attached_vect(val::AbstractVector{<:AbstractVector{<:T}}, ::Type{T}) where {T}
    vector = View{T}[]
    for i in eachindex(val[1])
        push!(vector, view(val[1], i))
    end
    return vector
end
attached_vect(val::AbstractVector{<:MutableNamedTuple}, ::Type{T}) where {T} =
    View{T}[attached_vect(val[1], T)...]
attached_vect(val, ::Type{T}) where {T} = View{T}[view(val, 1)]

recursive_convert(::Type{T}, num::Number) where {T} = T(num)
recursive_convert(::Type{T}, arr::AbstractArray) where {T} = recursive_convert.(T, arr)
function recursive_convert(::Type{T}, nt::NamedTuple) where {T}
    vals = []
    for (key, val) in zip(keys(nt), values(nt))
        push!(vals, key => recursive_convert(T, val))
    end
    return (; vals...)
end
recursive_convert(::Type{T}, mnt::MutableNamedTuple) where {T} =
    MutableNamedTuple(recursive_convert(T, namedtuple(mnt)))

Base.length(nvv::NamedViewVector) = length(getfield(nvv, :vector))

Base.size(nvv::NamedViewVector) = size(getfield(nvv, :vector))

Base.IndexStyle(nvv::NamedViewVector) = IndexLinear()

Base.getindex(nvv::NamedViewVector, i::Int) = getfield(nvv, :vector)[i][1][1]
# Base.getindex(nvv::NamedViewVector, I::Vararg{Int,N}) where {N} = getfield(nvv, :vector)[I...]
# Base.getindex(nvv::NamedViewVector, ::Colon) = getfield(nvv, :vector)[:]
Base.getindex(nvv::NamedViewVector, kr::AbstractRange) = getindex.(getfield(nvv, :vector)[kr], 1)
# Base.getindex(nvv::NamedViewVector, idx...) = getfield(nvv, :vector)[idx...]
Base.getindex(nvv::NamedViewVector, key::Symbol) = getfield(nvv, :data)[key]

function Base.setindex!(nvv::NamedViewVector, v, i::Int)
    getfield(nvv, :vector)[i][1] = v
    return nothing
end
# # Base.setindex!(nvv::NamedViewVector, v, I::Vararg{Int,N}) where {N} = (getfield(nvv, :vector)[I...] .= v)
function Base.setindex!(nvv::NamedViewVector, v, ::Colon)
    for i in eachindex(nvv)
        getfield(nvv, :vector)[i][1] = v
    end
    return nothing
end
function Base.setindex!(nvv::NamedViewVector, v, kr::AbstractRange)
    for i in kr
        getfield(nvv, :vector)[i][1] = v
    end
    return nothing
end
function Base.setindex!(nvv::NamedViewVector, v, key::Symbol)
    setproperty!(nvv, key, v)
    # getfield(nvv, :data)[key] = v
    return nothing
end

# Base.copyto!(nvv::NamedViewVector, src::AbstractArray) = (nvv[1:length(src)] = src)
# Base.copyto!(nvv::NamedViewVector, src::Base.Broadcast.Broadcasted) = (nvv[1:length(src)] = src)
# function Base.copyto!(nvv::NamedViewVector, doff, src, soff, N)
#     rng = 0:N-1
#     nvv[doff .+ rng] = src[soff .+ rng]
#     return nvv
# end

Base.getproperty(nvv::NamedViewVector, key::Symbol) = getfield(nvv, :data)[key]

function Base.setproperty!(nvv::NamedViewVector, key::Symbol, val)
    setproperty!(getfield(nvv, :data), key, val)
    # setindex!(getfield(getfield(getfield(nvv, :data), :data), key), val, 1)
    return nothing
end

Base.propertynames(nvv::NamedViewVector) = propertynames(getfield(nvv, :data))

Base.similar(nvv::NamedViewVector{T}) where {T} = NamedViewVector{T}(namedtuple(getfield(nvv, :data)))
Base.similar(nvv::NamedViewVector,::Type{T}) where T = NamedViewVector{T}(getfield(nvv, :data))

# Yikes. FIXME.
viewvect_bottom_type(vector) = map(x->x[1][1], vector) |> eltype
recursive_unitless_bottom_eltype(nvv::NamedViewVector) = eltype(nvv)
recursive_bottom_eltype(nvv::NamedViewVector) = recursive_unitless_bottom_eltype(nvv)
recursive_unitless_eltype(nvv::NamedViewVector) = recursive_unitless_bottom_eltype(nvv)
recursivecopy(nvv::NamedViewVector) = deepcopy(nvv)

Base.abs2(sa::SubArray) = abs2.(sa)
Base.:/(sa1::SubArray, sa2::SubArray) = sa1 ./ sa2
Base.isnan(sa::SubArray) = isnan.(sa)

vector(nvv::NamedViewVector) = getfield(nvv, :vector)

# Base.:+(nvv1::NamedViewVector, nvv2::NamedViewVector) = vector(nvv1) + vector(nvv2)
# Base.:*(x::Number, nvv::NamedViewVector) = x * vector(nvv)

function Base.show(io::IO, nvv::NamedViewVector)
    show(io, getfield(nvv, :data))
    return nothing
end
function Base.show(io::IO, ::MIME"text/plain", nvv::NamedViewVector)
    print(io, "NamedViewVector")
    show(io, getfield(nvv, :data))
    return nothing
end

max_stars_repo_name: jonniedie/ComposableDiffEqs

max_issues_repo_name: jonniedie/ComposableDiffEqs

max_forks_repo_name: jonniedie/ComposableDiffEqs


text:
using Turing
using Test
using Random

Random.seed!(125)

x = [1.5 2.0]

@model gdemo(x) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0,sqrt(s))
  for i in 1:length(x)
    x[i] ~ Normal(m, sqrt(s))
  end
  s, m
end

alg = IPMCMC(30, 500, 4)
chain = sample(gdemo(x), alg)
@test mean(chain[:s]) ≈ 49/24 atol=0.1
@test mean(chain[:m]) ≈ 7/6 atol=0.1


text:
mutable struct StandardSimplexe{T} <: AbstractSolver{T}
    M::Array{T, 2}
    Binv::Array{T, 2}
    xstar::Array{T, 1}
    vstar::T
    b_idx::Array{Int, 1}
    freelines::BitVector#index k is true if line k is linearly independent of the previous lines, dependent lines are found during phase 1.
    status::Status
    function StandardSimplexe{T}(M::Matrix{T}) where {T}
        @assert !(typeof(T) <: Integer)  "Type $T cannot be a subtype of Integer"
        m, n = size(M).-1
        ss = new{T}()
        ss.M = M
        ss.Binv = Array{T, 2}(I, m, m)
        ss.b_idx = -ones(T, m)
        ss.xstar = zeros(T, n)
        ss.status = Unknown
        ss.freelines = trues(m) #lines are inocent until proven guilty
        return ss
    end
end
function StandardSimplexe{T}(A::Matrix{T}, b::Vector{T}, c::Vector{T}) where {T}
    M = [A b; c' zero(T)]
    return  StandardSimplexe{T}(M)
end
function StandardSimplexe(A::Matrix{T}, b::Vector{T}, c::Vector{T}) where {T}
    M = [A b; c' zero(T)]
    return  StandardSimplexe{T}(M)
end
function StandardSimplexe(lp::AbstractLP{T}) where T
    A = getA(lp)
    b = getb(lp)
    c = getc(lp)
    return StandardSimplexe{T}(A, b, c)
end

function (ss::StandardSimplexe{T})(lp::AbstractLP{T}; verbose::Bool = false) where T
    m = size(ss.M, 1) - 1
    for k in 1:m
        if ss.M[k, end] <  0
            ss.Binv[k, :] *= -1
            ss.M[k, :] *= -1
        end
    end
    M = ss.M
    vM = @view M[1:end-1, :]
    iscanon, basis = iscanonical(vM)
    if iscanon
        ss.b_idx[:] = basis
        isOptimal(ss) && (ss.status = Optimal; lp.status = Optimal; lp.xstar[:] = ss.xstar; ss.vstar = lp.vstar; return ss.status)
        phase2!(ss, verbose = verbose)
    else
        phase1!(ss, verbose = verbose)
        isOptimal(ss) && (ss.status = Optimal; lp.status = Optimal; lp.xstar[:] = ss.xstar; ss.vstar = lp.vstar; return ss.status)
        (ss.status != Infeasible) && phase2!(ss, verbose = verbose)
    end
    lp.status = ss.status
    isOptimal(ss) && (lp.xstar[:] = ss.xstar; lp.vstar = ss.vstar)
    verbose && @show lp.vstar
    verbose && @show lp.xstar
    return ss.status
end

function println(ss::StandardSimplexe)
    m, n = size(ss.M)
    prd = "+" * prod("-------+" for _ in 1:n+1)
    println(prd)
    print("| V.B.  |")
    for k in 1:n-1
        s = ""
        xk = " x_$k "
        s*= xk
        s*= prod(" " for _ in 1:7-length(xk))
        print(s)
        print("|")
    end
    #@show ss.b_idx
    print(" T.D.  |")
    println()
    println(prd)
    for i in 1:m
        if i != m
            s = "| "#------"
            xk = ss.b_idx != [-1] ? "x_$(ss.b_idx[i])" : "   "
            s*= xk
            s*= prod(" " for _ in 1:6-length(xk))
            print(s)
        else
            print("| -z    ")
        end
        for j in 1:n
            print("| ")
            str = string(Float64(ss.M[i, j]))
            if length(str) >= 5
                str = str[1:5]
            end
            print(str)
            if length(1:(5 - length(str))) != 0
                print(prod(" "  for _ in 1:(5 - length(str))))
            end
            print(" ")
        end
        println("|")
        println(prd)
    end
end
function xstar(ss::StandardSimplexe{T})::Vector{T} where T
    ss.status != Optimal && @warn "probleme not optimal"
    return copy(ss.xstar)
end
function vstar(ss::StandardSimplexe{T})::T where T
    ss.status != Optimal && @warn "probleme not optimal"
    return ss.vstar
end
function isOptimal(ss::StandardSimplexe{T})::Bool where T
    return all(ss.M[end, 1:end-1] .>= -tol(T))
end

max_stars_repo_name: fbastin/Jasmin

max_issues_repo_name: fbastin/Jasmin

max_forks_repo_name: fbastin/Jasmin


text:
using CUDAnative: AS, DevicePtr

Base.@pure function datatype_align(::Type{T}) where {T}
    # typedef struct {
    #     uint32_t nfields;
    #     uint32_t alignment : 9;
    #     uint32_t haspadding : 1;
    #     uint32_t npointers : 20;
    #     uint32_t fielddesc_type : 2;
    # } jl_datatype_layout_t;
    field = T.layout + sizeof(UInt32)
    unsafe_load(convert(Ptr{UInt16}, field)) & convert(Int16, 2^9-1)
end

function kernel_vadd(a::DevicePtr{T,A}, b::DevicePtr{T,A}, c::DevicePtr{T,A}, i::Int) where {T,A}
    align = datatype_align(T)
    x = unsafe_load(a, i, Val(align)) + unsafe_load(b, i, Val(align))
    unsafe_store!(c, x, i, Val(align))
    return nothing
end

Base.code_llvm(stdout, kernel_vadd, Tuple{DevicePtr{Float32,AS.Global},
                                          DevicePtr{Float32,AS.Global},
                                          DevicePtr{Float32,AS.Global},
                                          Int})

using TinyCassette
struct Ctx end

Base.code_llvm(stdout, TinyCassette.execute,
               Tuple{Ctx, typeof(kernel_vadd),
                     DevicePtr{Float32,AS.Global},
                     DevicePtr{Float32,AS.Global},
                     DevicePtr{Float32,AS.Global},
                     Int})


text:
function day8()
    day8lines = filter(!isempty, strip.(readlines("AoCdata/AoC_2021_day8.txt")))
    part = [0, 0]

    patterns = [
        [String.(sort.(collect.(split(a)))) for a in split(line, "|")] for line in day8lines
    ]

    count1 = sum(count(length(x) == 2 for x in a[2]) for a in patterns)
    count4 = sum(count(length(x) == 4 for x in a[2]) for a in patterns)
    count7 = sum(count(length(x) == 3 for x in a[2]) for a in patterns)
    count8 = sum(count(length(x) == 7 for x in a[2]) for a in patterns)
    part[1] = sum([count1, count4, count7, count8])

    eachlinedig = Vector{String}[]
    for line in patterns
        dig = ["" for _ in 1:10]
        while any(dig .== "")
            for s in vcat(line...)
                length(s) == 2 && (dig[2] = s) # 1
                length(s) == 3 && (dig[8] = s) # 7
                length(s) == 4 && (dig[5] = s) # 4
                length(s) == 7 && (dig[9] = s) # 8
                if length(s) == 6 # 0, 6 or 9
                    if dig[5] != ""
                        if count(c -> c in s, dig[5]) == 4
                            dig[10] = s  # 9
                        elseif dig[2] != ""
                            if count(c -> c in s, dig[2]) == 2
                                dig[1] = s # 0
                            else
                                dig[7] = s # 6
                            end
                        end
                    end
                elseif length(s) == 5  # 2 or 3 or 5
                    if dig[2] != ""
                        if count(c -> c in s, dig[2]) == 2
                            dig[4] = s # 3
                        elseif dig[5] != ""
                            if count(c -> c in s, dig[5]) == 3
                                dig[6] = s # 5
                            else
                                dig[3] = s # 2
                            end
                        end
                    end
                end
            end
        end
        push!(eachlinedig, dig)
    end

    for (i, line) in enumerate(patterns)
        outputdigits = [findfirst(x -> x == s, eachlinedig[i]) - 1 for s in line[2]]
        part[2] += evalpoly(10, reverse(outputdigits))
    end

    println("Part 1:", part[1])
    println("Part 2:", part[2])
end

day8()


text:
function get_func_signature(fdef::Expr)
    @assert fdef.head == :function || fdef.head == :(=)
    sig , body = fdef.args
    @assert sig.head == :call && isa(sig,Expr)
    fname = sig.args[1]
    args = sig.args[2:end]

    args_parse = foldl(Tuple[],sig.args[2:end]) do acc, el
        if isa(el,Symbol)
            push!(acc,(el,true))
        else
            push!(acc,(el,false))
        end
    end
    calls = try
                open(string(fname),"r") do s
                    seekstart(s)
                    deserialize(s)
                end
            catch err
                if isa(err,EOFError) || isa(err,SystemError)
                    # empty, expected
                    Dict()
                else
                    println("Unexpected error reading file: ",err)
                    throw(err)
                end
            end
    println("Cache: for $fname: ",calls)
    Expr(:function,esc(sig),
         quote
            local args_tup = tuple($(args...))
            local key = foldl(Any[],enumerate($args_parse)) do acc, el
                local index = el[1]
                local symb = args_tup[index]
                if el[2][2]
                    if isempty(names(symb))
                        push!(acc,symb)
                    else
                        push!(acc,
                              [getfield(symb,n) for n in 1:length(names(symb))])
                    end
                else
                    push!(acc,symb)
                end
            end
            if haskey($calls,key)
                #println("Reading from file")
                $calls[key]
            else
                $calls[key] = let; $body; end
                open(string($fname) ,"w+") do s
                    serialize(s,$calls)
                end
                #println("Wrote to file")
                $calls[key]
            end
         end)
end

macro etcd_mock(func)
    get_func_signature(func)
end


text:
immutable BetaPrime <: ContinuousUnivariateDistribution
    betad::Beta

    BetaPrime(α::Float64, β::Float64) = new(Beta(α, β))
    BetaPrime(α::Float64) = new(Beta(α))
    BetaPrime() = new(Beta())
end

@distr_support BetaPrime 0.0 Inf

#### Parameters

params(d::BetaPrime) = params(d.betad)

#### Statistics

mean(d::BetaPrime) = ((α, β) = params(d); β > 1.0 ? α / (β - 1.0) : NaN)

mode(d::BetaPrime) = ((α, β) = params(d); α > 1.0 ? (α - 1.0) / (β + 1.0) : 0.0)

function var(d::BetaPrime)
    (α, β) = params(d)
    β > 2.0 ? α * (α + β - 1.0) / ((β - 2.0) * (β - 1.0)^2) : NaN
end

function skewness(d::BetaPrime)
    (α, β) = params(d)
    if β > 3.0
        s = α + β - 1.0
        2.0 * (α + s) / (β - 3.0) * sqrt((β - 2.0) / (α * s))
    else
        return NaN
    end
end

#### Evaluation

function logpdf(d::BetaPrime, x::Float64)
    (α, β) = params(d)
    (α - 1.0) * log(x) - (α + β) * log1p(x) - lbeta(α, β)
end

pdf(d::BetaPrime, x::Float64) = exp(logpdf(d, x))

cdf(d::BetaPrime, x::Float64) = cdf(d.betad, x / (1.0 + x))
ccdf(d::BetaPrime, x::Float64) = ccdf(d.betad, x / (1.0 + x))
logcdf(d::BetaPrime, x::Float64) = logcdf(d.betad, x / (1.0 + x))
logccdf(d::BetaPrime, x::Float64) = logccdf(d.betad, x / (1.0 + x))

quantile(d::BetaPrime, p::Float64) = (x = quantile(d.betad, p); x / (1.0 - x))
cquantile(d::BetaPrime, p::Float64) = (x = cquantile(d.betad, p); x / (1.0 - x))
invlogcdf(d::BetaPrime, p::Float64) = (x = invlogcdf(d.betad, p); x / (1.0 - x))
invlogccdf(d::BetaPrime, p::Float64) = (x = invlogccdf(d.betad, p); x / (1.0 - x))

#### Sampling

function rand(d::BetaPrime)
    (α, β) = params(d)
    rand(Gamma(α)) / rand(Gamma(β))
end


text:
export dense2padded!, padded2dense!

function dense2padded!{T}(backend::GPUBackend, dst::CuTensorBlob{T},
    src::CuTensorBlob{T}, pad_head::NTuple{2,Int}, pad_tail::NTuple{2,Int}, mirror::Bool)

    width, height, channels, num = size(src)
    chann_num = channels*num

    x_block = round(Int, ceil(float(chann_num)/CUDA.THREADS_PER_BLOCK_X))
    y_block = round(Int, ceil(float(width)/CUDA.THREADS_PER_BLOCK_Y))
    z_block = round(Int, ceil(float(height)/CUDA.THREADS_PER_BLOCK_Z))

    if T == Float32
      kernel = backend.mocha.dense_to_padded_float
    elseif T == Float64
      kernel = backend.mocha.dense_to_padded_double
    else
      error("Unsupported data type $T for padded copy")
    end
    mirror = convert(Int, mirror)
    CUDA.launch(kernel, (x_block,y_block,z_block),
        (CUDA.THREADS_PER_BLOCK_X,CUDA.THREADS_PER_BLOCK_Y,CUDA.THREADS_PER_BLOCK_Z),
        (dst.ptr.p,src.ptr.p,width,height,pad_head[1],pad_head[2],pad_tail[1],pad_tail[2],chann_num,mirror))
end

function padded2dense!{T}(backend::GPUBackend, dst::CuTensorBlob{T},
    src::CuTensorBlob{T}, pad_head::NTuple{2,Int}, pad_tail::NTuple{2,Int}, mirror::Bool)

    width, height, channels, num = size(dst)
    chann_num = channels*num

    x_block = round(Int, ceil(float(chann_num)/CUDA.THREADS_PER_BLOCK_X))
    y_block = round(Int, ceil(float(width)/CUDA.THREADS_PER_BLOCK_Y))
    z_block = round(Int, ceil(float(height)/CUDA.THREADS_PER_BLOCK_Z))

    if T == Float32
      kernel = backend.mocha.padded_to_dense_float
    elseif T == Float64
      kernel = backend.mocha.padded_to_dense_double
    else
      error("Unsupported data type $T for padded copy")
    end
    mirror = convert(Int, mirror)
    CUDA.launch(kernel, (x_block,y_block,z_block),
        (CUDA.THREADS_PER_BLOCK_X,CUDA.THREADS_PER_BLOCK_Y,CUDA.THREADS_PER_BLOCK_Z),
        (dst.ptr.p,src.ptr.p,width,height,pad_head[1],pad_head[2],pad_tail[1],pad_tail[2],chann_num, mirror))
end

function dense2padded!{T}(backend::GPUBackend, dst::CuTensorBlob{T},
    src::CuTensorBlob{T}, pad::NTuple{2,Int}, mirror=false)
  dense2padded!(backend, dst, src, pad, pad, mirror)
end
function padded2dense!{T}(backend::GPUBackend, dst::CuTensorBlob{T},
    src::CuTensorBlob{T}, pad::NTuple{2,Int}, mirror=false)
  padded2dense!(backend, dst, src, pad, pad, mirror)
end


text:
using DiffEqSensitivity, SafeTestsets
using Test, Pkg

const GROUP = get(ENV, "GROUP", "All")

function activate_gpu_env()
    Pkg.activate("gpu")
    Pkg.develop(PackageSpec(path=dirname(@__DIR__)))
    Pkg.instantiate()
end

@time begin

if GROUP == "All" || GROUP == "Core1" || GROUP == "Downstream"
    @time @safetestset "Forward Sensitivity" begin include("forward.jl") end
    @time @safetestset "Sparse Adjoint Sensitivity" begin include("sparse_adjoint.jl") end
    @time @safetestset "Second Order Sensitivity" begin include("second_order.jl") end
    @time @safetestset "Concrete Solve Derivatives" begin include("concrete_solve_derivatives.jl") end
    @time @safetestset "Branching Derivatives" begin include("branching_derivatives.jl") end
    @time @safetestset "Derivative Shapes" begin include("derivative_shapes.jl") end
    @time @safetestset "save_idxs" begin include("save_idxs.jl") end
    @time @safetestset "ArrayPartitions" begin include("array_partitions.jl") end
    @time @safetestset "Complex Adjoints" begin include("complex_adjoints.jl") end
    @time @safetestset "Forward Remake" begin include("forward_remake.jl") end
    @time @safetestset "Prob Kwargs" begin include("prob_kwargs.jl") end
    @time @safetestset "DiscreteProblem Adjoints" begin include("discrete.jl") end
    @time @safetestset "Time Type Mixing Adjoints" begin include("time_type_mixing.jl") end
end

if GROUP == "All" || GROUP == "Core2"
    @time @safetestset "hasbranching" begin include("hasbranching.jl") end
    @time @safetestset "Literal Adjoint" begin include("literal_adjoint.jl") end
    @time @safetestset "ForwardDiff Chunking Adjoints" begin include("forward_chunking.jl") end
    @time @safetestset "Stiff Adjoints" begin include("stiff_adjoints.jl") end
    @time @safetestset "Autodiff Events" begin include("autodiff_events.jl") end
    @time @safetestset "Null Parameters" begin include("null_parameters.jl") end
    @time @safetestset "Forward Mode Prob Kwargs" begin include("forward_prob_kwargs.jl") end
    @time @safetestset "Steady State Adjoint" begin include("steady_state.jl") end
    @time @safetestset "Concrete Solve Derivatives of Second Order ODEs" begin include("second_order_odes.jl") end
    @time @safetestset "Parameter Compatibility Errors" begin include("parameter_compatibility_errors.jl") end
end

if GROUP == "All" || GROUP == "Core3" || GROUP == "Downstream"
    @time @safetestset "Adjoint Sensitivity" begin include("adjoint.jl") end
end

if GROUP == "All" || GROUP == "Core4"
    @time @safetestset "Ensemble Tests" begin include("ensembles.jl") end
    @time @safetestset "GDP Regression Tests" begin include("gdp_regression_test.jl") end
    @time @safetestset "Layers Tests" begin include("layers.jl") end
    @time @safetestset "Layers SDE" begin include("layers_sde.jl") end
    @time @safetestset "Layers DDE" begin include("layers_dde.jl") end
    @time @safetestset "SDE - Neural" begin include("sde_neural.jl") end

    # No `@safetestset` since it requires running in Main
    @time @testset "Distributed" begin include("distributed.jl") end
end

if GROUP == "All" || GROUP == "Core5"
    @time @safetestset "Partial Neural Tests" begin include("partial_neural.jl") end
    @time @safetestset "Size Handling in Adjoint Tests" begin include("size_handling_adjoint.jl") end
    @time @safetestset "Callback - ReverseDiff" begin include("callback_reversediff.jl") end
    @time @safetestset "Hybrid DE" begin include("hybrid_de.jl") end
    @time @safetestset "HybridNODE" begin include("HybridNODE.jl") end
    @time @safetestset "ForwardDiff Sparsity Components" begin include("forwarddiffsensitivity_sparsity_components.jl") end
    @time @safetestset "Complex No u" begin include("complex_no_u.jl") end
end

if GROUP == "All" || GROUP == "SDE1"
    @time @safetestset "SDE Adjoint" begin include("sde_stratonovich.jl") end
    @time @safetestset "SDE Scalar Noise" begin include("sde_scalar_stratonovich.jl") end
    @time @safetestset "SDE Checkpointing" begin include("sde_checkpointing.jl") end
end

if GROUP == "All" || GROUP == "SDE2"
    @time @safetestset "SDE Non-Diagonal Noise" begin include("sde_nondiag_stratonovich.jl") end
end

if GROUP == "All" || GROUP == "SDE3"
    @time @safetestset "RODE Tests" begin include("rode.jl") end
    @time @safetestset "SDE Ito Conversion Tests" begin include("sde_transformation_test.jl") end
    @time @safetestset "SDE Ito Scalar Noise" begin include("sde_scalar_ito.jl") end
end

if GROUP == "Callbacks1"
    @time @safetestset "Discrete Callbacks with ForwardDiffSensitivity" begin include("callbacks/forward_sensitivity_callback.jl") end
    @time @safetestset "Discrete Callbacks with Adjoints" begin include("callbacks/discrete_callbacks.jl") end
    @time @safetestset "SDE Callbacks" begin include("callbacks/SDE_callbacks.jl") end
end

if GROUP == "Callbacks2"
    @time @safetestset "Continuous vs. discrete Callbacks" begin include("callbacks/continuous_vs_discrete.jl") end
    @time @safetestset "Continuous Callbacks with Adjoints" begin include("callbacks/continuous_callbacks.jl") end
    @time @safetestset "VectorContinuousCallbacks with Adjoints" begin include("callbacks/vector_continuous_callbacks.jl") end
end

if GROUP == "Shadowing"
    @time @safetestset "Shadowing Tests" begin include("shadowing.jl") end
end

if GROUP == "GPU"
    activate_gpu_env()
    @time @safetestset "Standard DiffEqFlux GPU" begin include("gpu/diffeqflux_standard_gpu.jl") end
    @time @safetestset "Mixed GPU/CPU" begin include("gpu/mixed_gpu_cpu_adjoint.jl") end
end
end


text:
module ADN_module_SIR

using Statistics

function infected_c(
    agent_state_1,
    agent_state_2,
    activity_good,
    activity_bad,
    m,
    mu,
    lambda,
    alpha,
    N_good,
    N_bad,
    # C_1,
    # C_2
)
    alpha_1=alpha
    alpha_2=.9
    N = N_good + N_bad
    AgentStateNext_1 = zeros(Int, N_good)
    AgentStateNext_2 = zeros(Int, N_bad)
    for k = 1:N_good
        AgentStateNext_1[k] = agent_state_1[k]
    end
    for k = 1:N_bad
        AgentStateNext_2[k] = agent_state_2[k]
    end

    for k = 1:N_good
        if agent_state_1[k] == 1
            if rand() < mu
                AgentStateNext_1[k] = 2
            end
        end
    end

    for k = 1:N_bad
        if agent_state_2[k] == 1
            if rand() < mu
                AgentStateNext_2[k] = 2
            end
        end
    end

    # loop over good agents
    for i = 1:N
        # if active
        # choose a random index
        if i <= N_good
            if rand() < activity_good[i]
                for j = 1:m
                    rInt = rand(1:N)

                    if rInt <= N_good

                        # the good active node gets the good neigbor sick
                        if agent_state_1[i] == 1
                            if (
                                (agent_state_1[rInt] == 0) &&
                                (rand() < alpha_1 *alpha_2 * lambda)
                            )
                                AgentStateNext_1[rInt] = 1
                            end
                        else
                            # agent_state_1[i] == 0
                            # the good active node gets sick from the good neighbor node
                            if (
                                (agent_state_1[rInt] == 1) &&
                                (rand() < alpha_1 *alpha_2* lambda)
                            )
                                AgentStateNext_1[i] = 1
                            end
                        end
                    else
                        rInt = rInt - N_good
                        if agent_state_1[i] == 1
                            if (
                                (agent_state_2[rInt] == 0) &&
                                (rand() < alpha_2 * lambda)
                            )
                                AgentStateNext_2[rInt] = 1
                            end
                        else
                            # agent_state_1[i] == 1
                            # the good active node gets sick from the good neighbor node
                            if (
                                (agent_state_2[rInt] == 1) &&
                                (rand() < alpha_2 * lambda)
                            )
                                AgentStateNext_1[i] = 1
                            end
                        end

                    end
                end
            end
        else
            i = i - N_good
            if rand() < activity_bad[i]
                for j = 1:m
                    rInt = rand(1:N)

                    if rInt <= N_good

                        # the good active node gets the good neigbor sick
                        if agent_state_2[i] == 1
                            if (
                                (agent_state_1[rInt] == 0) &&
                                (rand() < alpha_2 * lambda)
                            )
                                AgentStateNext_1[rInt] = 1
                            end
                        else
                            # agent_state_2[i] == 0
                            # the good active node gets sick from the good neighbor node
                            if (
                                (agent_state_1[rInt] == 1) &&
                                (rand() < alpha_2 * lambda)
                            )
                                AgentStateNext_2[i] = 1
                            end
                        end
                    else
                        rInt = rInt - N_good
                        if agent_state_2[i] == 1
                            if ((agent_state_2[rInt] == 0) && (rand() < lambda))
                                AgentStateNext_2[rInt] = 1
                            end
                        else
                            # agent_state_2[i] == 1
                            # the good active node gets sick from the good neighbor node
                            if ((agent_state_2[rInt] == 1) && (rand() < lambda))
                                AgentStateNext_2[i] = 1
                            end
                        end

                    end
                end
            end
        end
    end

    # if ((agent_state_2[rInt] == 0) && (rand() < alpha * lambda))
    #     AgentStateNext_2[rInt] = 1
    # end

    # C_1 = C_1 .| .!iszero.(AgentStateNext_1)
    # C_2 = C_2 .| .!iszero.(AgentStateNext_2)
    return AgentStateNext_1, AgentStateNext_2#, C_1, C_2
end

    function infected_totals(
    tmax,
    activity_good,
    activity_bad,
    m,
    mu,
    lambda,
    alpha,
    N_good,
    N_bad,
)
    agent_state_1 = zeros(Int, N_good)
    agent_state_2 = zeros(Int, N_bad)
    C_1 = .!iszero.(agent_state_1)
    C_2 = .!iszero.(agent_state_2)
    # Iend_1 = zeros(tmax)
    for k = 1:100
        agent_state_1[k] = 1
        agent_state_2[k] = 1
    end

    for i_time = 1:tmax
        agent_state_1, agent_state_2 = infected_c(
            agent_state_1,
            agent_state_2,
            activity_good,
            activity_bad,
            m,
            mu,
            lambda,
            alpha,
            N_good,
            N_bad,
            # C_1,
            # C_2
        )

    R=(sum(.!iszero.(agent_state_1)) + sum(.!iszero.(agent_state_2)))/(N_good+N_bad)
    T = sum([a==1 for a in agent_state_1]) + sum([a==1 for a in  agent_state_2])
    if R>.9 || T<4
        break
    end
    end

    return  (sum(.!iszero.(agent_state_1)) + sum(.!iszero.(agent_state_2)))/(N_good+N_bad)
end

   function parameter_sweep(n_trials, activity, grid)
    N = length(activity)
    proportion = grid[1]
    proportion = 1 - proportion
    # proportion =.1
    lambda = grid[2]
    alpha = .25
    N_bads = Int(ceil((proportion * N)))
    N_good = N - N_bads
    activity_bad = zeros(N_bads)
    activity_good = zeros(N_good)

    for k = 1:N_good
        activity_good[k] = activity[k]
    end
    for k = 1:N_bads
        activity_bad[k] = activity[k + N_good]
    end

    eta_ratio = 1

    # recovery rate
    mu = 0.05
    # infection rate
# *(1-.77 *proportion)

    m = 10
    tmax = 2000
    total_1 = 0.0
    total_2 = 0.0

    for k = 1:n_trials
        infected_1 = infected_totals(
            tmax,
            .8*activity_good,
            activity_bad,
            m,
            mu,
            lambda,
            alpha,
            N_good,
            N_bads,
        )
        total_1 = total_1 + infected_1
        # total_2 = total_2 + infected_2 / (N)

    end
    return total_1 / n_trials
end

end


text:
Base.append!(dest::TypedTables.Table, src::SingletonVector{<:NamedTuple}) =
    push!(dest, first(src))


text:
module Assignment03

export complement,
       ispurine,
       ispyrimidine,
       base_type,
       gc_content

include("assignment.jl")

end


text:
using Test
using Distributed
using Sockets
using Distributed: parse_machine, bind_client_port, SSHManager, LocalManager

@test parse_machine("127.0.0.1") == ("127.0.0.1", nothing)
@test parse_machine("127.0.0.1:80") == ("127.0.0.1", 80)
@test parse_machine("[2001:db8::1]") == ("2001:db8::1", nothing)
@test parse_machine("[2001:db8::1]:443") == ("2001:db8::1", 443)

@test parse_machine("127.0.0.1:90") == ("127.0.0.1", 90)
@test parse_machine("127.0.0.1:1") == ("127.0.0.1", 1)
@test parse_machine("127.0.0.1:65535") == ("127.0.0.1", 65535)
@test_throws ArgumentError parse_machine("127.0.0.1:-1")
@test_throws ArgumentError parse_machine("127.0.0.1:0")
@test_throws ArgumentError parse_machine("127.0.0.1:65536")

sock = bind_client_port(TCPSocket(), typeof(IPv4(0)))
addr, port = getsockname(sock)
@test addr == ip"0.0.0.0"

sock = bind_client_port(TCPSocket(), typeof(IPv6(0)))
addr, port = getsockname(sock)
@test addr == ip"::"

@test occursin(r"^SSHManager\(machines=.*\)$",
               sprint((t,x) -> show(t, "text/plain", x), SSHManager("127.0.0.1")))
@test sprint((t,x) -> show(t, "text/plain", x), LocalManager(1, true)) == "LocalManager()"

max_stars_repo_name: zierenberg/julia

max_issues_repo_name: zierenberg/julia

max_forks_repo_name: zierenberg/julia


text:
Base.@kwdef struct Grid{T}
    block::Block
    data::Matrix{T}
    cell_width::Int = 4
    cell_height::Int = 4

    function Grid(block, data::Matrix{T}, cell_width, cell_height) where T

        new{T}(block, data, 2 * (cell_width ÷ 2), 2 * (cell_height ÷ 2))

    end
end

function draw(grid::Grid, area::Rect, buf::Buffer)
    draw(grid.block, area, buf)
    grid_area = inner(grid.block, area)
    nrows, ncols = size(grid.data)

    ncols * grid.cell_width + 1 > width(grid_area) && return
    nrows * grid.cell_height + 1 > height(grid_area) && return

    for j in 1:ncols
        for i in 1:nrows
            x = left(grid_area) + grid.cell_width * 2 * (j-1) + grid.cell_width
            y = top(grid_area) + grid.cell_height * 1 * (i-1) + grid.cell_height ÷ 2

            set(buf, x, y, "$(grid.data[i, j])")

            for n in -grid.cell_width:grid.cell_width
                set(buf, x + n,  y - grid.cell_height ÷ 2, BOX_DRAWINGS_LIGHT_HORIZONTAL)
                set(buf, x + n,  y + grid.cell_height ÷ 2, BOX_DRAWINGS_LIGHT_HORIZONTAL)
            end

            for n in -grid.cell_height:grid.cell_height
                set(buf, x - grid.cell_width, y + n ÷ 2, BOX_DRAWINGS_LIGHT_VERTICAL)
                set(buf, x + grid.cell_width, y + n ÷ 2, BOX_DRAWINGS_LIGHT_VERTICAL)
            end
        end
    end

    for j in 1:ncols
        for i in 1:nrows
            x = left(grid_area) + grid.cell_width * 2 * (j-1) + grid.cell_width
            y = top(grid_area) + grid.cell_height * 1 * (i-1) + grid.cell_height ÷ 2
            set(buf, x + grid.cell_width, y + grid.cell_height ÷ 2, BOX_DRAWINGS_LIGHT_VERTICAL_AND_HORIZONTAL)
        end
    end

    for j in 1:ncols
        if j == ncols
            continue
        end

        x = left(grid_area) + grid.cell_width * 2 * (j-1) + grid.cell_width

        set(buf, x + grid.cell_width, top(grid_area), BOX_DRAWINGS_LIGHT_DOWN_AND_HORIZONTAL)
        set(buf, x + grid.cell_width, top(grid_area) + nrows * grid.cell_height, BOX_DRAWINGS_LIGHT_UP_AND_HORIZONTAL)

    end

    for i in 1:nrows
        if i == nrows
            continue
        end

        y = top(grid_area) + grid.cell_height * 1 * (i-1) + grid.cell_height ÷ 2

        set(buf, left(grid_area) + ncols * grid.cell_width * 2 , y + grid.cell_height ÷ 2, BOX_DRAWINGS_LIGHT_VERTICAL_AND_LEFT)
        set(buf, left(grid_area), y + grid.cell_height ÷ 2, BOX_DRAWINGS_LIGHT_VERTICAL_AND_RIGHT)

    end

    set(buf, left(grid_area), top(grid_area), BOX_DRAWINGS_LIGHT_DOWN_AND_RIGHT)
    set(buf, left(grid_area), top(grid_area) + nrows * grid.cell_height, BOX_DRAWINGS_LIGHT_UP_AND_RIGHT)

    set(buf, left(grid_area) + ncols * grid.cell_width * 2, top(grid_area), BOX_DRAWINGS_LIGHT_DOWN_AND_LEFT)
    set(buf, left(grid_area) + ncols * grid.cell_width * 2, top(grid_area) + nrows * grid.cell_height, BOX_DRAWINGS_LIGHT_UP_AND_LEFT)
end


text:
function keysuri(cli::Client, key::String)
    return "http://$(cli.host):$(cli.port)/$(cli.version)/keys$(key)"
end

"""
    get(cli, key, sort=false, recursive=false) -> Dict

Fetches the value for the specified `key` from the etcd cluster.
See the [etcd API](https://github.com/coreos/etcd/blob/master/Documentation/v2/api.md#get-the-value-of-a-key)
for more details.
"""
function Base.get(cli::Client, key::String; sort::Bool=false, recursive::Bool=false)
    return get(cli, key, Dict{String, Any}(); sort=sort, recursive=recursive)
end

function Base.get(cli::Client, key::String, opts::Dict{String, Any}; sort::Bool=false, recursive::Bool=false)
    opts["sorted"] = sort
    opts["recursive"] = recursive
    opts = filter!((key, val) -> !isa(val, Bool) || val, opts)
    return request(Requests.get, keysuri(cli, key), opts)
end

function put(cli::Client, key::String, opts::Dict{String, Any}; ttl::Int=-1)
    if ttl > 0
        opts["ttl"] = ttl
    end

    return request(Requests.put, keysuri(cli, key), opts)
end

function post(cli::Client, key::String, opts::Dict{String, Any}; ttl::Int=-1)
    if ttl > 0
        opts["ttl"] = ttl
    end

    return request(Requests.post, keysuri(cli, key), opts)
end

function delete(cli::Client, key::String, opts::Dict{String, Any})
    return request(Requests.delete, keysuri(cli, key), opts)
end

"""
    set(cli, key, value; ttl=-1, ordered=false) -> Dict

Sets the `value` for the specified `key` in the etcd cluster.
See the [etcd API](https://github.com/coreos/etcd/blob/master/Documentation/v2/api.md#changing-the-value-of-a-key)
for more details.
"""
function set(cli::Client, key::String, value::String; ttl::Int=-1, ordered=false)
    opts = Dict{String, Any}("value" => value)

    if ordered
        return post(cli, key, opts; ttl=ttl)
    else
        return put(cli, key, opts; ttl=ttl)
    end
end

"""
    setdir(cli, key; ttl=-1) -> Dict

Creates a directory with the specified `key` on the etcd cluster.

Reference: https://github.com/coreos/etcd/blob/master/Documentation/v2/api.md#creating-directories
"""
function setdir(cli::Client, key::String; ttl::Int=-1)
    opts = Dict{String, Any}("value" => "", "dir" => true)
    return put(cli, key, opts; ttl=ttl)
end

"""
    create(cli, key, value; ttl=-1) -> Dict

Creates a new key/value pair on the etcd cluster, asserting that the key does not already exist.
"""
function create(cli::Client, key::String, value::String; ttl::Int=-1)
    opts = Dict{String, Any}("value" => value, "prevExist" => false)
    return put(cli, key, opts; ttl=ttl)
end

"""
    update(cli, key, value; ttl=-1) -> Dict

Updates an existing key with a new value on the etcd cluster, asserting that the key must already exist.
"""
function update(cli::Client, key::String, value::String; ttl::Int=-1)
    opts = Dict{String, Any}("value" => value, "prevExist" => true)
    return put(cli, key, opts; ttl=ttl)
end

"""
    createdir(cli, key; ttl=-1) -> Dict

Creates a new directory on the etcd cluster, asserting that the directory does not already exist.
"""
function createdir(cli::Client, key::String; ttl::Int=-1)
    opts = Dict{String, Any}("value" => "", "prevExist" => false, "dir" => true)
    return put(cli, key, opts; ttl=ttl)
end

"""
    updatedir(cli, key; ttl=-1) -> Dict

Updates the name of a direcotry on the etcd cluster, asserting that the directory already exists.
"""
function updatedir(cli::Client, key::String; ttl::Int=-1)
    opts = Dict{String, Any}("value" => "", "prevExist" => true, "dir" => true)
    return put(cli, key, opts; ttl=ttl)
end

"""
    exists(cli, key) -> Bool

Returns whether the `key` exists in the etcd cluster.
"""
function exists(cli::Client, key::String)
    try
        Base.get(cli, key)
        return true
    catch err
        # expect "Key not found"(100) error
        if err.resp["errorCode"] == 100
            return false
        else
            rethrow()
        end
    end
end

Base.haskey(cli::Client, key::String) = exists(cli, key)

"""
    delete(cli, key) -> Dict

Removes the `key` from the etcd cluster.
"""
delete(cli::Client, key::String) = delete(cli, key, Dict{String,Any}())

"""
    deletedir(cli, key; recursive=false) -> Dict

Removes the directory from the etcd cluster.
"""
function deletedir(cli::Client, key::String; recursive::Bool=false)
    opts = Dict{String, Any}("dir" => true, "recursive" => recursive)
    delete(cli, key, filter!((k,v)->v, opts))
end

"""
    cad(cli, key; prev_value=nothing, prev_index=-1) -> Dict

Performs an atomic compare-and-delete on the `key` in the etcd cluster.
See the [etcd API](https://github.com/coreos/etcd/blob/master/Documentation/v2/api.md#atomic-compare-and-delete)
for more details.
"""
function cad(
    cli::Client, key::String;
    prev_value::Union{String, Void}=nothing, prev_index::Int=-1,
)
    opts = Dict{String, Any}()
    if prev_value !== nothing
       opts["prevValue"] = prev_value
    end

    if prev_index > 0
       opts["prevIndex"] = prev_index
    end

    return delete(cli, key, opts)
end

"""
    cas(cli, key, value; prev_value=nothing, prev_index=-1, ttl=-1) -> Dict

Performs an atomic compare-and-swap with the `key` and `value` on the etcd cluster.
See the [etcd API](https://github.com/coreos/etcd/blob/master/Documentation/v2/api.md#atomic-compare-and-swap)
for more details.
"""
function cas(
    cli::Client, key::String, value::String;
    prev_value::Union{String, Void}=nothing, prev_index::Int=-1, ttl::Int=-1
)
    opts = Dict{String, Any}("value" => value)

    if prev_value !== nothing
       opts["prevValue"] = prev_value
    end

    if prev_index > 0
       opts["prevIndex"] = prev_index
    end

    return put(cli, key, opts; ttl=ttl)
end


text:
###################################################################################
# Antony Schutz 2015, ANR - POLCA - 2016
###################################################################################
# Pre calculus, MATRIX Inversion - Method involving orthogonal matrices
# ---------------------------------------------------------------------------------
function invmat_par(F3D::Array,rho_y::Real,eta::Real,nw::Int)#,paral::Bool)
# step IV of [0]
# inverse of inner matrix of C^\lambda_n
# [F F^H + eta/rho I]^-1
# Compute the inverse matrix of size nb*nb for image reconstruction, pre calculus step
# F3D: non uniform 3D Fourier Transform Matrix
# rho_x: admm parameters
# eta: (nwvlt*OIDATA.rho_spat+OIDATA.rho_spec+OIDATA.rho_ps+OIDATA.epsilon)
    mat2inv = [(F3D[n], (rho_y / eta)) for n in 1:nw ]
    M = pmap(toinv, mat2inv)
    return M
end
# ---------------------------------------------------------------------------------
# for parallel calculus of inverse
function toinv(mat2inv)
   return inv((mat2inv[2] * mat2inv[1] * mat2inv[1]') + speye(size(mat2inv[1], 1)))
end
# ---------------------------------------------------------------------------------
# Phases To Phases difference Matrix
# ---------------------------------------------------------------------------------
# Method 1
function phasetophasediff(Closure_index::Matrix,nw::Integer,nb::Integer,T3::Integer,DP::Integer,dptype::ASCIIString,dpprm::Integer)
# create the matrix which link the phases to the phases difference
# Eq. 20 of [1] , 8 of [0]
# Closure phase and differential phase to phase matrix
# closure_index: index of phase closure (index related to phase vector index)
# nb is the number of base
# nw is the number of wavelength
# T3 and DP TBD for later, inform if TP and DP are present in data
# defined as in equation 14 of PAINTER
# dptype is the kind of differential phases: all, ref, diag, frame, sliding
# dpprm  is the parameters associated to ref (lambda ref) or frame and sliding (horizon)

    nbnw = (nb * nw)
    if(T3 == 1)
      NT3 = size(Closure_index, 1)
      rowt = repeat(collect(1:(NT3 * nw)), inner = [3])
      colt = repeat(vec(Closure_index'), outer = [nw]) + repeat(collect(0:(nw - 1)) * nb, inner = [3 * NT3])
      Valuet = repeat([1, 1, -1], outer=[NT3 * nw])
    end

# Differential phase to phase matrix
# defined as in equation 19 of PAINTER [1], 7 of [0]
# rank = Nbas*(Nwvl-1)
    if(DP == 1)
      if dptype == "all"
          # rank = Nbas*( Nwvl - 1 )
          HDP = alldp(nb, nw)
      elseif dptype == "ref"
          # rank = Nbas*( Nwvl - 1 )
          # dpprm is the reference channel
          HDP = alldpref(nb, nw, dpprm)
      elseif dptype == "diag"
          # rank = Nbas*( Nwvl - 1 )
          HDP = diagdp(nb, nw)
      elseif dptype == "frame"
          # rank = Nbas*floor( Nwvl - Nwvl/horizon )
          # dpprm is the horizon size of window
          HDP = framedp(nb, nw, dpprm)
      elseif dptype == "sliding"
          # rank = Nbas*( Nwvl - 1 )
          # dpprm is the horizon size of window
          HDP = slidingdp(nb, nw, dpprm)
      elseif dptype == "phase"
          # rank = Nbas*Nwvl (full)
          HDP = nodp(nb, nw)
      else
          error("wrong choice for Differential Phase to Phase Matrix: all (default), frame, ref, diag, sliding or 'phase' ")
      end
    end

# phase difference to phase matrix
# defined as in equation 20 of PAINTER [1], 6 of [0]
    if((T3 == 1)&&(DP == 0))
      HT3 = sparse(rowt, colt, Valuet, (NT3 * nw), nbnw)
      return HT3

    elseif((T3 == 0)&&(DP == 1))
      return HDP

    elseif((T3 == 1)&&(DP == 1))
      HT3 = sparse(rowt, colt, Valuet, NT3 * nw, nbnw)
      H = vcat(HT3, HDP)
      return H
    end
end
# All differential phase, Painter V1
function alldp(nb::Int64,nw::Int64)
    nbnw = (nb * nw)
    rowb = repeat(collect(1:(nb * (nw - 1))), inner = [2])
    colb = vec(vcat(repmat(collect(1:nb)', 1, (nw-1)),collect((nb + 1):nbnw)'))
    Valueb = repeat(vcat(1,-1), outer=[nb*(nw-1)])
    HDP = sparse(rowb, colb, Valueb, nb * (nw-1), nbnw)
    return HDP
end
# case visphi is phases of complexe visibilities
function nodp(nb::Int64,nw::Int64)
    nbnw = (nb * nw)
    HDP = speye(nbnw)
    return HDP
end
# All differential phase, with reference channel defined by lambdaref
function alldpref(nb::Int64,nw::Int64,lambdaref::Int64)
    tmp = (lambdaref-1)*nb + collect(1:nb)'
    nbnw = (nb * nw)
    rowb = repeat(collect(1:(nb * (nw - 1))), inner = [2])

    colb = vec( vcat( repmat(tmp, 1, nw - 1), hcat( collect( 1:(tmp[1] - 1) )', collect( (1 + tmp[end]):nbnw )' ) ))
    Valueb = repeat(vcat(1,-1), outer=[nb*(nw-1)])
    HDP = sparse(rowb, colb, Valueb, nb * (nw-1), nbnw)
   return HDP
end
# All differential phase, Painter V1
function diagdp(nb::Int64,nw::Int64)
    nbnw = (nb * nw)
    tmp = collect(1:(nb * (nw - 1)))
    rowb = vcat(tmp,tmp)
    colb = vcat(tmp,collect((nb + 1):nbnw))
    Valueb = repeat(vcat(1,-1), inner=[nb*(nw-1)])
    HDP = sparse(rowb, colb, Valueb, nb * (nw-1), nbnw)
    return HDP
end

# Non overlapping windows for DP, Painter SPIE 2016
function framedp(nb::Int64,nw::Int64,horiz::Int64)
    horiz = horiz-1;
    nbloc   = floor(nw/(horiz+1));
    rowb    = Int64[]
    Colb    = Int64[]
    Valueb  = Int64[]
    for n   in 1:nbloc
        rowbtm = (n-1)*nb*(horiz) + repeat(collect(1:(nb * horiz)), inner = [2])
        Colbtm = (n-1)*nb*(horiz+1) + vcat( repmat(collect(1:nb)',1,horiz), collect(nb+1:nb*(horiz+1))' )
        Valuebtm = repmat(vcat(1,-1),nb*horiz,1)
        rowb    = vcat(rowb, round(Int,rowbtm))
        Colb    = vcat(Colb, vec(round(Int,Colbtm)))
        Valueb  = vcat(Valueb, vec(round(Int,Valuebtm)))
    end
    # % rest
    horiz2  = round(Int, (nb*nw-maximum( vec(Colb) ))/nb-1 )
    if horiz2>0
        rowbtm  = round(Int, maximum( vec(rowb) )) +repeat(collect(1:(nb * horiz2)), inner = [2])
        Colbtm  = round(Int, maximum( vec(Colb) )) + vcat( repmat(collect(1:nb)',1,horiz2), collect(nb+1:nb*(horiz2+1))' )
        Valuebtm = repmat(vcat(1,-1),nb*horiz2,1)
        rowb    = vcat(rowb, rowbtm)
        Colb    = vcat(Colb, vec(Colbtm))
        Valueb  = vcat(Valueb, vec(Valuebtm))
    end
    H       = sparse(rowb, Colb, Valueb, maximum(rowb),nb*nw)
    return H
end

# sliding windows, overlap 1, for DP, Painter SPIE 2016
function slidingdp(nb::Int64,nw::Int64,horiz::Int64)
    horiz = horiz-1;
    nbnw = (nb * nw)
    row1 = repeat(collect(1:(nb*horiz)), inner = [2])
    col1 = vec(vcat(repmat(collect(1:nb)', 1, horiz),collect((nb + 1):(nb*(horiz+1)))'))
    value1 = repeat(vcat(1,-1), outer=[nb*horiz])

    row2    = Int64[]
    col2    = Int64[]
    value2  = Int64[]

    vnb = collect(1:nb)
    for n   in 1:(nw-1-horiz)
        rowbtm = (n-1)*nb + maximum(row1) + repeat(vnb, inner = [2])
        colbtm = vcat(n*nb + vnb' , (n-1)*nb + maximum(col1) + vnb')
        valuebtm = repmat(vcat(1,-1),nb,1)
        row2 = vcat(row2, round(Int,rowbtm))
        col2 = vcat(col2, vec(round(Int,colbtm)))
        value2 = vcat(value2, vec(round(Int,valuebtm)))
    end

    row = vcat(row1, row2)
    col = vcat(col1, vec(col2))
    value = vcat(value1, vec(value2))
    H = sparse(row, col, value, nb*(nw-1), nb*nw)
    return H
end

###################################################################################
# Tools for NUFFT
# ---------------------------------------------------------------------------------
function planarray_par(tab_u::Array,tab_v::Array,nx::Int,nw::Int)#,parral::Bool)
# create Array for plan for non uniform fft
# tab_u and tab_v are spatial frequencies
# nx is images side size
# nw is the number of wavelength
    Mtmp = [( hcat(tab_u[:,n] / nx, tab_v[:,n] / nx)', nx) for n in 1:nw]
    return pmap(toplan, Mtmp)
end
# ---------------------------------------------------------------------------------
# for parallel calculus
function toplan(Mtmp)
    NFFTPlan(Mtmp[1], (Mtmp[2], Mtmp[2]))
end
###################################################################################
# Matrices involved in data model
# ---------------------------------------------------------------------------------
# Non Uniform Discrete Fourier Matrices for Polychromatic data
# ---------------------------------------------------------------------------------
function nudft3d_par(tab_u::Array,tab_v::Array,nb::Int,nx::Int,nw::Int)
# create a list of non uniform dft Matrix per channel
# 3D nudft [number of bases (spatiale frequencies) per wavelength * number of pixels * number of wavelength]
# tab_u and tab_v must be of the form : [number of bases * number of wavelength]
# nx is images side size
# nw is the number of wavelength
# defined as in equation 3 of PAINTER [1]
    if(nx < 1)
        error("nudft3d:Image size must be positive")
    end
    nfu = length(tab_u)
    nfv = length(tab_v)

    if((nfu == nfv)&&(nfu > 0)&&(nfv > 0))

# test if U and V have same size and have same number of wavelength, otherwise PB
        if(size(tab_u,2) == size(tab_v,2))
            UVMAT = [(tab_u[:, n], tab_v[:, n], nx, nb) for n in 1:nw ]
	          F3D = pmap(non_uniform_dft_par, UVMAT)
            return F3D

        else
            error("nudft3d: U ad V vectors have not of same number of wavelength")
        end

    else
        error("nudft3d: U ad V vectors are not of same size")
    end
end
# ---------------------------------------------------------------------------------
function non_uniform_dft_par(UVMAT)
# Non Uniform DFT Matrix:
# UVMAT [u,v,nx]
# for a nx*nx pixels image
# u, v vector of normalized spatial frequencies (between -0.5 and 0.5)
# UVMAT[1] : U
# UVMAT[2] : V
# UVMAT[3] : nx
# UVMAT[4] : nb
    k1  	= repmat(collect((-UVMAT[3] / 2):(-1 + (UVMAT[3] / 2))), 1, UVMAT[3])
    k1m 	= repmat((vec(k1 )' / UVMAT[3]), UVMAT[4], 1)
    k2m 	= repmat((vec(k1')' / UVMAT[3]), UVMAT[4], 1)
    um 		= repmat(vec(UVMAT[1]), 1, (UVMAT[3] * UVMAT[3]))
    vm 		= repmat(vec(UVMAT[2]), 1, (UVMAT[3] * UVMAT[3]))
    return exp(-2 * im * pi * ((k1m .* um) + (k2m .* vm) )) ./ UVMAT[3]
end


text:
using FITSIO
using Base.Test

# -----------------------------------------------------------------------------
# Images

# Create a FITS instance and loop over supported types.
fname = tempname() * ".fits"
f = FITS(fname, "w")
for T in [Uint8, Int8, Uint16, Int16, Uint32, Int32, Int64,
          Float32, Float64]
    indata = reshape(T[1:100;], 5, 20)

    # Test writing the data to a new extension
    write(f, indata)

    # test reading the full array
    outdata = read(f[end])
    @test indata == outdata
    @test eltype(indata) == eltype(outdata)

    # test reading subsets of the array
    @test f[end][:, :] == indata
    @test f[end][4, 1:10] == indata[4, 1:10]  # 2-d array
    @test f[end][:, 4] == indata[:, 4]  # 1-d array
    @test f[end][2, 3] == indata[2, 3]  # scalar
    @test f[end][:, 1:2:10] == indata[:, 1:2:10]
    @test f[end][1:3, :] == indata[1:3, :]

    # test expected errors
    @test_throws DimensionMismatch f[end][:]
    @test_throws DimensionMismatch f[end][:, :, 1]
    @test_throws BoundsError f[end][1:6, :]
    @test_throws BoundsError f[end][1, 0]

end
close(f)
if isfile(fname)
    rm(fname)
end

# copy_section()
fname1 = tempname() * ".fits"
f1 = FITS(fname1, "w")
indata = reshape(Float32[1:400;], 20, 20)
write(f1, indata)

fname2 = tempname() * ".fits"
f2 = FITS(fname2, "w")
copy_section(f1[1], f2, 1:10, 1:10)
copy_section(f1[1], f2, 1:10, 1:2:20)
outdata = read(f2[1])
@test outdata == indata[1:10, 1:10]
outdata = read(f2[2])
@test outdata == indata[1:10, 1:2:20]
close(f1)
close(f2)
if isfile(fname1)
    rm(fname1)
end
if isfile(fname2)
    rm(fname2)
end

# -----------------------------------------------------------------------------
# FITSHeader

fname = tempname() * ".fits"
f = FITS(fname, "w")
inhdr = FITSHeader(["FLTKEY", "INTKEY", "BOOLKEY", "STRKEY", "COMMENT",
                    "HISTORY"],
                   [1.0, 1, true, "string value", nothing, nothing],
                   ["floating point keyword",
                    "",
                    "boolean keyword",
                    "string value",
                    "this is a comment",
                    "this is a history"])

inhdr["INTKEY"] = 2  # test setting by key
inhdr[1] = 2.0  # test settting by index
setcomment!(inhdr, "INTKEY", "integer keyword") # test setting a comment

indata = reshape(Float32[1:100;], 5, 20)
write(f, indata; header=inhdr)
outhdr = readheader(f[1])
@test outhdr["FLTKEY"] === 2.0
@test outhdr["INTKEY"] === 2
@test outhdr["BOOLKEY"] === true
@test outhdr["STRKEY"] == "string value"
@test getcomment(outhdr, 13) == "this is a comment"
@test getcomment(outhdr, 14) == "this is a history"
@test length(outhdr) == 14
@test haskey(outhdr, "FLTKEY")

# Read single keywords
@test readkey(f[1], 9) == ("FLTKEY", 2.0, "floating point keyword")
@test readkey(f[1], "FLTKEY") == (2.0, "floating point keyword")

# Test that show() works and that the beginning of output is what we expect.
io = IOBuffer()
show(io, f)
s = takebuf_string(io)
@test s[1:6] == "file: "

# Clean up from last test.
close(f)
if isfile(fname)
    rm(fname)
end

# -----------------------------------------------------------------------------
# parsing non-standard keyword records

# `create_test_file` : Create a simple FITS file for testing, with the
# given header string added after the required keywords. The length of
# `header` must be a multiple of 80.  The purpose of creating such
# files is to test the parsing of non-standard FITS keyword records
# (non-standard files can't be created with cfitsio).

function create_test_file(fname::String, header::ASCIIString)
    if length(header) % 80 != 0
        error("length of header must be multiple of 80")
    end

    f = open(fname, "w")

    stdhdr = "SIMPLE  =                    T / file does conform to FITS standard             BITPIX  =                  -64 / number of bits per data pixel                  NAXIS   =                    2 / number of data axes                            NAXIS1  =                   10 / length of data axis 1                          NAXIS2  =                   10 / length of data axis 2                          EXTEND  =                    T / FITS dataset may contain extensions            "
    endline = "END                                                                             "
    data = fill(0., (10, 10))  # 10x10 array of big-endian Float64 zeros

    # write header
    write(f, stdhdr)
    write(f, header)
    write(f, endline)

    # add padding
    block_position = (length(stdhdr) + length(header) + length(endline)) % 2880
    padding = (block_position == 0) ? 0 : 2880 - block_position
    write(f, " "^padding)

    # write data
    write(f, data)

    # add padding
    block_position = sizeof(data) % 2880
    padding = (block_position == 0) ? 0 : 2880 - block_position
    write(f, fill(0x00, (padding,)))

    close(f)
end

# Create a test file with a few non-standard keyword records
fname = tempname() * ".fits"
header = "        Warning: CROTA2 is inaccurate due to considerable skew                  SKEW    =  1.9511305786508E+00,  1.9234924037208E+00 /Measure of skew           "
create_test_file(fname, header)

# check that we can read the header (and data).
f = FITS(fname)
hdr = readheader(f[1])
data = read(f[1])
close(f)

# test that we can write it back out.
fname2 = tempname() * ".fits"
f2 = FITS(fname2, "w")
write(f2, data; header=hdr)
close(f2)

# clean up.
if isfile(fname)
    rm(fname)
end
if isfile(fname2)
    rm(fname2)
end

println("All tests passed.")


text:
function bw_digit(font::String, accent::String, digit, font_size::Int, img_size::Tuple)
    Drawing(img_size...)
    background("black")
    sethue("white")
    p = Point(6, 28)
    setfont(font*accent, font_size)
    settext("$digit", p)
    img = image_as_matrix()
    finish()
    img
end

function rgb_digit(font::String, accent::String, digit, font_size::Int, img_size::Tuple, col)
    Drawing(img_size...)
    background("black")
    sethue(col)
    p = Point(6, 28)
    setfont(font*accent, font_size)
    settext("$digit", p)
    img = image_as_matrix()
    finish()
    img
end

img_to_array_bw(img) = reshape(Float32.(channelview(Gray.(img))), size(img)...,1)
img_to_array_rgb(img) = permutedims(Float32.(channelview(RGB.(img))), (2,3,1))
array_to_img_bw(arr) = Gray.(reshape(permutedims(arr, (3,1,2)), size(arr)[1:2]...))
array_to_img_rgb(arr) = RGB.(arr[:,:,1]', arr[:,:,2]', arr[:,:,3]')

function transform(img, rot, shearx, sheary, zoom)
    pl = Rotate([rot]) |>
            ShearX([shearx]) |>
            ShearY([sheary]) |>
            CropSize(28, 28) |>
            Zoom([zoom])
   augment(img, pl)
end

function create_digit_bw(font_size, img_size, font, accent, dig, rot, shearx, sheary, zoom)
	orig = bw_digit(font, accent, dig, font_size, img_size)
	new = transform(orig, rot, shearx, sheary, zoom)
	img_to_array_bw(new)
end

function create_digit_rgb(font_size, img_size, font, accent, dig, col, rot, shearx, sheary, zoom)
    orig = rgb_digit(font, accent, dig, font_size, img_size, col)
    new = transform(orig, rot, shearx, sheary, zoom)
    img_to_array_rgb(new)
end

vectorize(x::AbstractArray{T,4}) where T = reshape(x, :, size(x,4))
devectorize(x::AbstractArray{T,2}, w, h) where T = reshape(x, w, h, 1, :)
devectorize(x::AbstractArray{T,2}, w, h, c) where T = reshape(x, w, h, c, :)
devectorize(x::AbstractArray{T,4}, args...; kwargs...) where T = x

function draw(x::AbstractArray{T,4}) where T
    if size(x,3) == 1
        return hcat(array_to_img_bw.([x[:,:,:,i] for i in 1:size(x,4)])...)
    elseif size(x,3) == 3
        return hcat(array_to_img_rgb.([x[:,:,:,i] for i in 1:size(x,4)])...)
    else
        error("only implemented for number of channels equal to 1 or 3")
    end
end
function draw(x::AbstractArray{T,3}) where T
    if size(x,3) == 1
        return array_to_img_bw(x)
    elseif size(x,3) == 3
        return array_to_img_rgb(x)
    else
        error("only implemented for number of channels equal to 1 or 3")
    end
end
function draw(x::AbstractArray{T,2}, s=(28,28)) where T
    draw(devectorize(x, s...))
end

max_stars_repo_name: vitskvara/HierarchicalAD

max_issues_repo_name: vitskvara/HierarchicalAD

max_forks_repo_name: vitskvara/HierarchicalAD