text:
using GAMSFiles

""" Turns GAMSFiles.GCall into an Expr. """
function eq_to_expr(eq::GAMSFiles.GCall, sets::Dict{String, Any})
    @assert(length(eq.args)==2)
    lhs, rhs, op = eq.args[1], eq.args[2], GAMSFiles.eqops[GAMSFiles.getname(eq)]
    lhsexpr = convert(Expr, lhs)
    rhsexpr = convert(Expr, rhs)
    GAMSFiles.replace_reductions!(lhsexpr, sets)
    GAMSFiles.replace_reductions!(rhsexpr, sets)
    if op in [:(=), :>]
        :($(lhsexpr) - $(rhsexpr))
    elseif op == :<
        :($(rhsexpr) - $(lhsexpr))
    else
        error("unexpected operator ", op)
    end
end

""" Returns whether a GAMSFiles.GCall is an equality. """
function is_equality(eq::GAMSFiles.GCall)
    @assert(length(eq.args)==2)
    GAMSFiles.eqops[GAMSFiles.getname(eq)] == :(=)
end

""" Finds and returns all varkeys in equation. """
function find_vars_in_eq(eq::GAMSFiles.GCall, vardict::Dict{Symbol, Any})
    vars = Symbol[]
    lhs, rhs = eq.args[1], eq.args[2]
    for (var, vinfo) in vardict
        if GAMSFiles.hasvar(lhs, string(var)) || GAMSFiles.hasvar(rhs, string(var))
            push!(vars, var)
        end
    end
    return vars
end

""" Takes gams and turns them into JuMP.Variables"""
function generate_variables!(model::JuMP.Model, gams::Dict{String, Any})
    gamsvars = GAMSFiles.allvars(gams)
    vardict = Dict{Symbol, Any}()
    constdict = Dict{Symbol, Any}()
    for (var, vinfo) in gamsvars
        if vinfo isa Union{Array, Real}
            model[Symbol(var)] = vinfo
            constdict[Symbol(var)] = vinfo
        elseif vinfo isa Base.RefValue
            model[Symbol(var)] = vinfo.x
            constdict[Symbol(var)] = vinfo.x
        else
            axs = vinfo.axs
            if axs == ()
                nv = @variable(model, base_name = var)
                model[Symbol(var)] = nv
                vardict[Symbol(var)] = nv
            else
                nv = JuMP.Containers.DenseAxisArray{JuMP.VariableRef}(undef, axs...)
                for idx in eachindex(nv)
                    nv[idx] = @variable(model)
                    set_name(nv[idx], "$(var)[$(join(Tuple(idx),","))]")
                end
                model[Symbol(var)] = nv
                vardict[Symbol(var)] = nv
            end
            if vinfo.typ != "free"
                if vinfo.typ == "positive"
                    JuMP.set_lower_bound.(model[Symbol(var)], 0)
                elseif vinfo.typ == "negative"
                    JuMP.set_upper_bound.(model[Symbol(var)], 0)
                elseif vinfo.typ == "binary"
                    JuMP.set_binary.(model[Symbol(var)])
                    JuMP.set_lower_bound.(model[Symbol(var)], 0)
                    JuMP.set_upper_bound.(model[Symbol(var)], 1)
                elseif vinfo.typ == "integer"
                    JuMP.set_integer.(model[Symbol(var)])
                else
                    throw(OCTHaGOnException("Type $(vinfo.typ) unknown for variable $(var)."))
                end
            end
            for (prop, val) in vinfo.assignments
                if hasproperty(prop, :indices)
                    inds = map(x->x.val, prop.indices)
                    nv = model[Symbol(var)][inds...]
                    c = val.val
                    if prop.name =="l"
                        JuMP.set_start_value(nv, c)
                    elseif prop.name == "fx"
                        JuMP.set_lower_bound(nv, c)
                        JuMP.set_upper_bound(nv, c)
                    elseif prop.name == "lo"
                        JuMP.set_lower_bound(nv, c)
                    elseif prop.name == "up"
                        JuMP.set_upper_bound(nv, c)
                    end
                else
                    nv = model[Symbol(var)]
                    c = val.val
                    if prop.text == "l"
                        JuMP.set_start_value(nv, c)
                    elseif prop.text == "fx"
                        JuMP.set_lower_bound(nv, c)
                        JuMP.set_upper_bound(nv, c)
                    elseif prop.text == "lo"
                        JuMP.set_lower_bound(nv, c)
                    elseif prop.text == "up"
                        JuMP.set_upper_bound(nv, c)
                    end
                end
            end
        end
    end
    return vardict, constdict
end

"""
    GAMS_to_GlobalModel(GAMS_DIR::String, filename::String)

Converts a GAMS optimization model to a GlobalModel.
GAMS_DIR is the directory to look in, while the filename
is the name of the .gms file.
"""
function GAMS_to_GlobalModel(GAMS_DIR::String, filename::String)
    model = JuMP.Model()
    # Parsing GAMS Files
    lexed = GAMSFiles.lex(GAMS_DIR * filename)
    gams = GAMSFiles.parsegams(GAMS_DIR * filename)
    GAMSFiles.parseconsts!(gams)

    vars = GAMSFiles.getvars(gams["variables"])
    sets = Dict{String, Any}()
    if haskey(gams, "sets")
        sets = gams["sets"]
    end
    preexprs, bodyexprs = Expr[], Expr[]
    if haskey(gams, "parameters") && haskey(gams, "assignments")
        GAMSFiles.parseassignments!(preexprs, gams["assignments"], gams["parameters"], sets)
    end

    # Getting variables
    vardict, constdict = generate_variables!(model, gams) # Actual JuMP variables
    obj_symb = :objvar
    # Getting objective
    if "minimizing" in keys(gams)
        if gams["minimizing"] isa String
            @objective(model, Min, JuMP.variable_by_name(model, gams["minimizing"]))
            obj_symb = Symbol(gams["minimizing"])
        else
            @objective(model, Min, sum([JuMP.variable_by_name(i) for i in gams["minimizing"]]))
        end
    elseif "maximizing" in keys(gams)
        @warn "$(filename) is a maximization. Make sure objvar is on LHS of constraints, with a positive coefficient and and equality."
        if gams["maximizing"] isa String
            @objective(model, Max, JuMP.variable_by_name(model, gams["maximizing"]))
            obj_symb = Symbol(gams["maximizing"])
        else
            @objective(model, Max, sum([JuMP.variable_by_name(i) for i in gams["maximizing"]]))
        end
    end

    # Creating GlobalModel
    gm = GlobalModel(model = model, name = filename)
    equations = [] # For debugging purposes...
    for (key, eq) in gams["equations"]
        push!(equations, key => eq)
    end
    for (key, eq) in equations
        if key isa GAMSFiles.GText
            constr_expr = eq_to_expr(eq, sets)
            # Substitute constant variables
            constkeys = find_vars_in_eq(eq, constdict)
            const_pairs = Dict(constkey => model[constkey] for constkey in constkeys)
            for (constkey, constval) in const_pairs
                constr_expr = OCTHaGOn.substitute(constr_expr, :($constkey) => constval)
            end
            # Designate free variables
            varkeys = find_vars_in_eq(eq, vardict)
            if !(obj_symb in varkeys)
                vars = Array{VariableRef}(flat([vardict[varkey] for varkey in varkeys]))
                input = Symbol.(varkeys)
                constr_fn = :(($(input...),) -> $(constr_expr))
                if length(input) == 1
                    constr_fn = :($(input...) -> $(constr_expr))
                end
                add_nonlinear_or_compatible(gm, constr_fn, vars = vars, expr_vars = [vardict[varkey] for varkey in varkeys],
                                        equality = is_equality(eq), name = gm.name * "_" * GAMSFiles.getname(key))
            else
                constr_expr = OCTHaGOn.substitute(constr_expr, :($(obj_symb)) => 0)
                # ASSUMPTION: objvar has positive coefficient, and is on the greater size.
                op = GAMSFiles.eqops[GAMSFiles.getname(eq)]
                # if !(op in [:<, :>])
                    # throw(OCTHaGOnException("Please make sure GAMS model has objvar on the greater than size of inequalities, " *
                    #                     " with a leading coefficient of 1."))
                # end
                varkeys = filter!(x -> x != obj_symb, varkeys)
                vars = Array{VariableRef}(flat([vardict[varkey] for varkey in varkeys]))
                input = Symbol.(varkeys)
                constr_fn = :(($(input...),) -> -$(constr_expr))
                if length(input) == 1
                    constr_fn = :($(input...) -> -$(constr_expr))
                end
                add_nonlinear_or_compatible(gm, constr_fn, vars = vars, expr_vars = [vardict[varkey] for varkey in varkeys],
                    dependent_var = vardict[obj_symb], equality = is_equality(eq), name = gm.name * "_" * GAMSFiles.getname(key))
            end
        elseif key isa GAMSFiles.GArray
            axs = GAMSFiles.getaxes(key.indices, sets)
            idxs = collect(Base.product([collect(ax) for ax in axs]...))
            names = [gm.name * "_" * key.name * string(idx) for idx in idxs]
            constr_expr = eq_to_expr(eq, sets)
                    # Substitute constant variables
            constkeys = find_vars_in_eq(eq, constdict)
            const_pairs = Dict(Symbol(constkey) => model[constkey] for constkey in constkeys)
            for (constkey, constval) in const_pairs
                constr_expr = OCTHaGOn.substitute(constr_expr, :($constkey) => constval)
            end
            # Designate free variables
            varkeys = find_vars_in_eq(eq, vardict)
            vars = Array{VariableRef}(flat([vardict[varkey] for varkey in varkeys]))
            input = Symbol.(varkeys)
            constr_fn = :(($(input...),) -> $(constr_expr))
            if length(input) == 1
                constr_fn = :($(input...) -> $(constr_expr))
            end
            constr_fns = []
            for idx in idxs
                new_fn = copy(constr_fn)
                for ax_number in 1:length(axs)
                    new_fn = OCTHaGOn.substitute(new_fn, Symbol(key.indices[ax_number].text) => idx[ax_number]);
                end
                push!(constr_fns, new_fn)
            end
            for i = 1:length(idxs)
                add_nonlinear_or_compatible(gm, constr_fns[i], vars = vars, expr_vars = [vardict[varkey] for varkey in varkeys],
                                         equality = is_equality(eq),
                                         name = names[i])
            end
        end
    end
    return gm
end

# """ Solver wrapper for GAMS benchmarking. """
# function nonlinear_solve(gm::GlobalModel; solver = IPOPT_SILENT)
#     nonlinearize!(gm)
#     set_optimizer(gm, solver)
#     optimize!(gm)
# end


text:
#=
    struct wrapped character constants for NanoDate formats
=#

# '⦰' is the reversed empty set 0x29b0
Base.isempty(ch::Char) = ch === EmptyChar

const CapitalT = 'T'      # ISO selected char, separates date from time
const SmallCapitalT = 'ᴛ' # latin letter small capital T 0x1d1b
const PunctDot = '.'      # separates fractions of a second, subseconds

const EmptyChar = '⦰'
const SingleSpace = ' '
const Underscore = '_'
const SmallWhiteCircle = '◦'
const SmallWhiteStar = '⭒'

#=
        make a string representation from a NanoDate
=#

Base.string(nd::NanoDate; sep::CharString=EmptyChar) =
    sep === EmptyChar ? nanodate_string(nd) : nanodate_string(nd, sep)

Base.string(nd::NanoDate, sep) = string(nd; sep=sep)

function nanodate_string(nd)
    str = string(nd.datetime)
    micronanos = value(nd.nanosecs)
    micronanos == 0 && return str
    millis = millisecond(nd.datetime)
    if millis == 0 # micronanos != 0
        str *= ".000"
    end
    micros, nanos = divrem(micronanos, NanosecondsPerMicrosecond)
    if nanos == 0 # micros != 0
        str *= lpad(micros, 3, '0')
    elseif micros == 0 # nanos != 0
        str *= lpad(nanos, 6, '0')
    else # micros != 0 && nanos != 0
        str *= lpad(micros, 3, '0')
        str *= lpad(nanos, 3, '0')
    end
    str
end

function nanodate_string(nd::NanoDate, sep::CharString)
    dtm = nd.datetime
    str = string(nd.datetime)
    datepart, timepart = split(str, CapitalT)
    nanos = value(nd.nanosecs)
    nanos === 0 && return str
    micros, nanos = fldmod(nanos, 1_000)

    millis = millisecond(nd.datetime)
    if millis === 0
        str = str * ".000"
    end

    if nanos === 0
        padded = sep * lpad(micros, 3, '0')
    elseif micros === 0
        padded = sep * "000" * sep * lpad(nanos, 3, '0')
    else
        padded = sep * lpad(micros, 3, '0') * sep * lpad(nanos, 3, '0')
    end
    str * padded
end

Dates.format(nd::NanoDate, df::DateFormat=NANODATE_FORMAT;
             sep::CharString=EmptyChar) =
    sep === EmptyChar ? nanodate_format(nd, df) :
                        nanodate_format(nd, df, sep)

function nanodate_format(nd, df)
    datetime = nd.datetime
    str = Dates.format(datetime, df)
    value(nd.nanosecs) == 0 && return str

    millis = millisecond(datetime)
    millistr = PunctDot * lpad(millis, 3, '0')
    str = split(str, PunctDot)[1]

    nsubsecfields = 0
    lasttoken = df.tokens[end]
    if isa(lasttoken, Dates.DatePart) &&
       typeof(lasttoken).parameters[1] === 's'
        nsubsecfields = lasttoken.width
    end
    nsubsecfields == 0 && return str
    str = str * millistr
    nsubsecfields == 1 && return str
    nanos = value(nd.nanosecs)
    cs, ns = divrem(nanos, 1_000)
    str = str * lpad(cs, 3, '0')
    nsubsecfields == 2 && return str
    str = str * lpad(ns, 3, '0')
    str
end

function nanodate_format(nd, df, sep)
    datetime = nd.datetime
    str = Dates.format(datetime, df)
    value(nd.nanosecs) == 0 && return str

    millis = millisecond(datetime)
    millistr = PunctDot * lpad(millis, 3, '0')
    str = split(str, PunctDot)[1]

    nsubsecfields = 0
    lasttoken = df.tokens[end]
    if isa(lasttoken, Dates.DatePart) &&
       typeof(lasttoken).parameters[1] === 's'
        nsubsecfields = lasttoken.width
    end
    nsubsecfields == 0 && return str
    str = str * millistr
    nsubsecfields == 1 && return str
    nanos = value(nd.nanosecs)
    cs, ns = divrem(nanos, 1_000)
    str = str * sep * lpad(cs, 3, '0')
    nsubsecfields == 2 && return str
    str = str * sep * lpad(ns, 3, '0')
    str
end

function Time(s::String)
    n = length(s)
    n <= 12 && return Dates.Time(s)
    timefromstring(s)
end

function timefromstring(s::String)
    dotidx = findlast(x -> x == '.', s)
    if isnothing(dotidx) || length(s) - dotidx <= 3
        return parse(Time, s)
    else
        secs, subsecs = split(s, '.')
        secs *= "." * subsecs[1:3]
        subsecs = subsecs[4:end]
        tm = parse(Time, secs)
        isubsecs = Meta.parse("1" * subsecs)
        micros, nanos = divrem(isubsecs, 1_000)
        micros = mod(micros, 1_000)
        return tm + Microsecond(micros) + Nanosecond(nanos)
    end
end

ext: jl


text:
using JuliaSandbox

using Test

function foo(m::Module)
    m.foo()
end

foo(JuliaSandbox.InternalModule)

@testset "JuliaSandbox.jl" begin
    println("TODO")
    # Write your tests here.
end

max_stars_repo_name: erwanlecarpentier/JuliaSandbox

max_issues_repo_name: erwanlecarpentier/JuliaSandbox

max_forks_repo_name: erwanlecarpentier/JuliaSandbox


text:
println("Fitting Volterra series model")
n,k,λ= 19,5,logspace(0,6,20)[15]
M3 = RegularizedVolterraModel(n,N1,k,λ)

volterra = estfun(M3,p1,Q1)


text:
using RegionCovarianceDescriptor
using Test

@testset "RegionCovarianceDescriptor.jl" begin
    # Write your own tests here.
end


text:
### A Pluto.jl notebook ###
# v0.12.7

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 9bc349e8-efb0-11ea-1392-cd735828af33
# packaged loads, etc
begin
	using DelimitedFiles
	using NIfTI
	using Images
	using ImageView
	using StatsBase
	using Plots
	using PlutoUI
	using Printf
	using LinearAlgebra
	using CSV
end

# ╔═╡ df5232e0-04e2-11eb-0c26-f19b9d9c7b14
begin
	using GLM
	using DataFrames
end

# ╔═╡ 8029070e-efb0-11ea-1d0a-7f402f2e9989
md"""# Images for moodle quiz """

# ╔═╡ 4df55a88-f98a-11ea-1c9f-cd52a4b5062f
md"""## General idea

Set of N different plots with voxel timecourses that look
"""

# ╔═╡ 2b2f6be4-f07f-11ea-1459-775749ca949f
md"""## Load in data """

# ╔═╡ 368e2640-f8e0-11ea-1ddf-b37069ba7588
d, h= readdlm("data/stimulus_timing.csv", ',', Float64, header = true )

# ╔═╡ 51594b8e-0741-11eb-0c03-338b971e945f
df = DataFrame(d,[h...]) # read as dataframe

# ╔═╡ bf362532-0741-11eb-137e-67c2924f4f92
fmri_data = niread("data/filtered_func_data.nii");

# ╔═╡ db903b82-0741-11eb-2df7-0dcebbb93c01
# this is how you get the header info!
fmri_data.header

# ╔═╡ 19bc9784-0742-11eb-18a4-edca29549ab4

# ╔═╡ cac0339a-f8e4-11ea-0876-f91dc7db1454
md"""## Have a quick look through data"""

# ╔═╡ 1058934e-f8e1-11ea-1896-35d75c4b6a8e
begin
md"""
*Make a slider that controls which column to plot*

**column:** $(@bind whichCol Slider(2:size(d,2), default=2, show_value=true))
"""
end

# ╔═╡ c6cab78c-f836-11ea-2b5b-23f0006f8e93
begin
	plot(d[:,1], d[:,whichCol], xlabel="Time (volumes)",
		ylabel="fMRI response", label="$(h[whichCol])", xlims=(0,160), ylims=(-5,5))
end

# ╔═╡ da4ac85c-f8e4-11ea-192e-03d5b6862d5a
md"""## Scaling = multiplying"""

# ╔═╡ 39d465a6-f8e2-11ea-2d71-75cbe586428c
function plotAtCoords(theData, theX, theY, theZ)
	p = plot(
		plot(1:160, view(theData,theX, theY, theZ, :),
		xlim=(0,160),
			c=:black,
			lw=2,
			grid_lw = 1,
		xlabel="Time (volumes)",
		ylabel="fMRI response (%signal change)",
			xticks=0:16:160,
			xminorticks=2,
			label=string("data at [$theX, $theY, $theZ]"))
		)
	p
end

# ╔═╡ a233a2e4-07c7-11eb-07d1-5561b4cc64b2
function plotWithStimuli(fmri_data, xs, ys, zs)
	p = plotAtCoords(fmri_data, xs, ys, zs)
	plot(p)
	xdiff = 16

	xtloc1 = xdiff.+(0:32:144)
	xtloc2 = xdiff.+(16:32:160)

	yl = ylims()
	ydiff = (yl[2]-yl[1])./40
	ann = [ (xtloc1, yl[1]+ydiff, Plots.text("face", 10, :red, :right, :bottom )),
			(xtloc2, yl[1]+ydiff, Plots.text("obj", 10, :blue, :right, :bottom ))]

	plot!(; annotation=ann)
	#annotate!([(0, 20, Plots.text("this is #5", 16, :red, :center))])
	#annotate!(0, 20, "bla") #GR.annotate!(0, 20, "my text", :color)
end

# ╔═╡ 1c3b6c56-f8e2-11ea-3969-f59871025c5e
begin
	sz = convert.(Int64,size(fmri_data).//2)
md"""
How can we *scale* the idealised prediction to fit the data better?

**x:** $(@bind xs Slider(1:size(fmri_data,1), default=sz[1], show_value=true))

**y:** $(@bind ys Slider(1:size(fmri_data,2), default=sz[2], show_value=true))

**z:** $(@bind zs Slider(1:size(fmri_data,3), default=1, show_value=true))

	"""
end

# ╔═╡ 45e986f0-074c-11eb-0d49-dbce692ade1a
begin
	heatmap(view(fmri_data, :, :, zs, 1),
		aspect_ratio=1,
		c=:greys, colorbar=:none)
	xaxis!(showaxis=:no)
	plot!([xs], [ys], m=:square,
		msw=4, msc=:red, ms=4, c=:red)
end

# ╔═╡ 665f1302-07b8-11eb-0848-8fb424cae65c
p = plotAtCoords(fmri_data, xs, ys, zs)

# ╔═╡ 2522ee1a-07ca-11eb-1d15-5f6b9e8f3aa6
illustration = plotWithStimuli(fmri_data, xs, ys, zs)

# ╔═╡ bae1b4e6-07ae-11eb-1fb9-c56c0cf7298a
md"""
$(@bind fileroot TextField(;default="testfile"))

$(@bind test Button("clicked"))
"""

# ╔═╡ e8716624-07b8-11eb-0219-81e10bd83100
filename = @sprintf("%s%+02d%+02d%+02d.pdf", fileroot, xs, ys, zs)

# ╔═╡ fe67814c-2299-11eb-2080-79f06d510355
save(filename, illustration)

# ╔═╡ 5176ef4e-07ba-11eb-09e7-9b6894fcc592
begin
	test
	allFilenames = [];
	for (xx,yy,zz) in [(18,14,4), (17,13,4), (19,9,4), (22,23,2), (20,21,2), (40,17,2), (36,1,5),(46,11,1),(19,11,2)]
		sfilename = @sprintf("%s%+02d%+02d%+02d.pdf", fileroot, xx+1, yy+1, zz+1)
		#q = plotAtCoords(fmri_data, xx+1, yy+1, zz+1)
		q = plotWithStimuli(fmri_data, xx+1, yy+1, zz+1)
		save(sfilename, q)
		append!(allFilenames, [sfilename])
	end
	with_terminal(print, allFilenames)

end

# ╔═╡ 63820b34-07af-11eb-21c3-75356e38ec43
begin
	test
	save(filename, p)
	md""" because `test` hooked up to `Button` this is reactive, **$(filename)**"""
end

# ╔═╡ 417c2dba-07b1-11eb-18cc-61d23caecbc6
rand(1:9999, 1,1)

# ╔═╡ 7e11e66e-07b3-11eb-1186-fba413764c97
with_terminal(dump, [1,2,[3,4]])

# ╔═╡ 7c2633f6-2297-11eb-16f1-7b4f1b1a283c
# 4 each for now F, O, Both / NB coords here in FSL 0-indexed form!, filenames will be 1-indexd!
begin
quizVoxels = [
	(17,13,4),
	(18,14,5),
	(18,8,4),
	(18,14,4),

	(40,20,3),
	(40,19,2),
	(22,23,2),
	(40,17,1),

	(36,10,1),
	(37,9,1),
	(37,9,2),
	(25,7,2)];

	quizVoxels3 = quizVoxels[[12,1,6]]
end

# ╔═╡ 6a1c2364-2297-11eb-36b5-7f100ff97687
# for the 2020/21 lab report:
begin
	allFilenamesAssess = [];
	for (xx,yy,zz) in quizVoxels
		sfilename = @sprintf("%s%+02d%+02d%+02d.pdf", fileroot, xx+1, yy+1, zz+1)
		#q = plotAtCoords(fmri_data, xx+1, yy+1, zz+1)
		q = plotWithStimuli(fmri_data, xx+1, yy+1, zz+1)
		save(sfilename, q)
		append!(allFilenamesAssess, [sfilename])
	end
	with_terminal(print, allFilenamesAssess)

end

# ╔═╡ 7d28a38a-229a-11eb-2cec-2153fbbe4e63
# for the 2020/21 lab report:
begin
	if true
	allFilenamesCSV = [];
	for (xx,yy,zz) in quizVoxels3 show
		with_terminal(print, (xx,yy,zz))
		sfilename = @sprintf("%s%+02d%+02d%+02d.csv", "lab-report" , xx+1, yy+1, zz+1)
		# writeVoxelDataToFile(fmri_data, xx+1, yy+1, zz+1, sfilename)
		append!(allFilenamesCSV, [sfilename])
	end
	with_terminal(print, allFilenamesCSV)
	end

end

# ╔═╡ 946112fc-229b-11eb-3027-bd10d12cb59b
normalizeTimecourse = function(tc)
	tc = 100.0 .* (tc./mean(tc) .- 1)
end

# ╔═╡ 32372d56-22c7-11eb-0923-4db4ed166e26
@bind quizPlot Slider(1:12, show_value=true)

# ╔═╡ 47b7c70c-229b-11eb-303a-41674232ea90
begin
	tm = normalizeTimecourse(fmri_data[1 .+ quizVoxels[quizPlot]...,:])
	plot(tm, label="quiz: $quizPlot, $(1 .+ quizVoxels[quizPlot])")
end

# ╔═╡ da7f81ca-229a-11eb-138f-43dd54c21ad5
writeVoxelDataToFile = function(fmri_data, xx, yy, zz, sfilename)

	t = [collect(0.0:159.0).*1., normalizeTimecourse( fmri_data[xx+1, yy+1, zz+1,:] ) ]
	df = DataFrame(t, [:Time_seconds, :fMRI_response])
	CSV.write(sfilename, df)
end

# ╔═╡ fc1de00a-229b-11eb-22f3-13b7f379ed7d
writeVoxelDataToFile(fmri_data, quizVoxels[1]..., "bla")

# ╔═╡ 2315ce8c-fb57-11ea-31a4-eb68d1d096d3
md""" ### Determining what fits best

Plot data, fit and errors together to understand what number `m` gives the best fit"""

# ╔═╡ d9f6c1e8-fa89-11ea-05b3-07daaa06326f
begin
md"""
How can we *scale* the idealised prediction to fit the data better?

**multiplier:** $(@bind m Slider(-10:0.5:10, default=1, show_value=true))
"""
end

# ╔═╡ a1b031a4-fc08-11ea-2314-13e2e79c2ca3
md"""### Now: more than one component to fit"""

# ╔═╡ e6309486-fc08-11ea-3ade-c93a24b6866b
begin
	plot(d[:,1], d[:,3],
		xlim=(0,160), ylim=(-1, 1),
		xlabel="Time (volumes)",
		ylabel="fMRI response (%signal change)",
		label="response (FACES)",
		color=:red, lw = 2, layout = (2, 1), subplot=1)
	plot!(d[:,1], d[:,4],
		color=:blue, lw = 2,
		label="response (OBJECTS)", subplot=2,xlim=(0,160), ylim=(-1, 1))

end

# ╔═╡ 9c3b0258-fc08-11ea-1922-fd1bec66b4ff
begin
md"""
How can we *scale* the idealised prediction to fit the data better?

**multiplier (faces):** $(@bind mf Slider(-10:0.5:10, default=1, show_value=true))

**multiplier (objects):** $(@bind mo Slider(-10:0.5:10, default=1, show_value=true))

"""
end

# ╔═╡ 975860d2-fc08-11ea-2def-ad1c45826911
begin
	plot(d[:,1], d[:,5],
		xlim=(0,160), ylim=(-5, 5),
		xlabel="Time (volumes)",
	ylabel="fMRI response (%signal change)", label="data")
	plot!(d[:,1], d[:,[3,4]]*[mf; mo], lw=2,
		label="model * [w_1; w_2]")

end

# ╔═╡ e7fb832e-04e0-11eb-388d-f7520b672a6e
md"""

## Calculating the weights for the best fit

*Linear regression* can be used to find/calculate the weights that will give the **best** fit of the mixture of the two curves.

**Best** here means the settings that will result in the *least squares fit*, the fit that makes the *sum of squared errors* between each measurement in the timeseries ($\text{data}_i$) and the corresponding model description ($\text{model}_i$) smallest.

$\sum_i (\text{data}_i - \text{model}_i)^2$

"""

# ╔═╡ 2b070720-04e2-11eb-3101-cbc40b17682b
# X, the design matrix
designMatrix = d[:,3:4]

# ╔═╡ 4b968440-04e2-11eb-35ce-5f050dcd4b6b
y_1 = d[:,5] #

# ╔═╡ 6464f806-04e2-11eb-37cc-7363fa76c80e
β_1 = designMatrix \ y_1

# ╔═╡ 30b5258c-04e5-11eb-3cdc-11d36a2f35bc
html"""<br>"""

# ╔═╡ 111aea9a-04e5-11eb-066d-4d245b04dc5e
md"""
##

turn data into a `dataframe` and use `GLM` module to get a linear model fit

"""

# ╔═╡ 1e370abc-04e3-11eb-064a-21fbb9fafc7c
theData = DataFrame(d[:, 3:6], [:face_response, :object_response, :Voxel1, :Voxel2 ])

# ╔═╡ 1d60545a-04e4-11eb-1720-b5f04c5da239
ols1 = lm(@formula(Voxel1 ~ face_response + object_response), theData)

# ╔═╡ 67b724d0-04e4-11eb-1254-f59bf78b1345
ols2 = lm(@formula(Voxel2 ~ face_response + object_response), theData)

# ╔═╡ 86d10676-04e4-11eb-2121-2f0c8945cc1a
coefnames(ols1)

# ╔═╡ 979732dc-04e4-11eb-3803-7302dc91e777
coef(ols1)

# ╔═╡ a33a18f2-04e4-11eb-0587-6117c036c641
GLM.confint(ols1)

# ╔═╡ 31071f6a-081a-11eb-2eda-f99e2073028c

# ╔═╡ c3899100-04e4-11eb-37aa-394d3f60bd53
import Pkg; Pkg.add("StatsPlots")

# ╔═╡ d5c0b236-f98b-11ea-0444-6fe340642707
md"""

## Some settings and helper functions

"""

# ╔═╡ 28cb2832-0748-11eb-0c14-832f1df303ca

# ╔═╡ 79d08b08-fa8b-11ea-0b4b-01ceef03660a
# calculate sum of squared errors for different betas
begin

	betaVals = -10:0.5:10 ;#range(-10, 10, length=nPoints)
	e = zeros(size(betaVals))
	X = Vector(d[:,3])
	theY = Vector(d[:,5])
	for (i, theB) in enumerate(betaVals)
		evec = X * theB .- theY
		e[i] = transpose(evec)* evec
	end
	#e
	md"""calculate squared errors"""
end

# ╔═╡ 58c63ca8-fa8b-11ea-296d-0b7f613e8976
begin
	# plot error curve and current choice of M
function plotError(betas, es, currentM, subplotNum)
	plot!(betas, es, lw=2, color=:black,
		label="error", subplot=subplotNum,)
	idx = betas .== currentM
	if es[idx]... ≈ minimum(es)
		theColor = :green
		theLabel = "min SSE"
	else
		theColor = :red
		theLabel = "SSE"

	end
	plot!(betas[idx], es[idx],
		marker=:circle,
		ms=10,
		color= theColor,
		subplot=subplotNum,
		xlabel= "β value",
		ylabel = "SSE",
		label = theLabel)
end
	md"""a function that plots error for a given choice of m `plotError` """
end

# ╔═╡ 96844d74-fa87-11ea-0966-89bbbc8a4353
begin
function plotWithBeta(xData,yData,yModel, m)
	l = @layout [ grid(2,1) a{0.30w}
					 ]
	plot(xData, yData,
		xlim=(0,160), ylim=(-5, 5),
		ylabel="fMRI response (%change)",
		label="data",
		layout = l,
	subplot=1)
	plot!(xData, yModel*m, lw=2,
		subplot=1,
		label="model * $(@sprintf "%.1f" m )")
	sticks!(xData, yData .- yModel*m,
			color=:gray,
			label="data-model",
				xlabel="Time (volumes)",
		subplot=2,
	xlim=(0,160), ylim=(-5, 5))

	# betaVals from other cell
	# ditto for e
	plotError(betaVals, e, m, 3)
end
	md"""a function that plots data + model with a
	particular choice of m `plotWithBeta` """
end

# ╔═╡ bfdd7d56-f98a-11ea-3064-a9033ba994a5
#plotly()
gr()

# ╔═╡ eda9e2dc-f98b-11ea-3eca-539ac938107c
nPoints = 20; # number of points in plot above (not used)

# ╔═╡ f6e4754c-f98b-11ea-2a42-034d82f52ca1
addNoise = 0.5 .* randn(nPoints);

# ╔═╡ Cell order:
# ╠═8029070e-efb0-11ea-1d0a-7f402f2e9989
# ╠═9bc349e8-efb0-11ea-1392-cd735828af33
# ╠═4df55a88-f98a-11ea-1c9f-cd52a4b5062f
# ╠═2b2f6be4-f07f-11ea-1459-775749ca949f
# ╠═368e2640-f8e0-11ea-1ddf-b37069ba7588
# ╠═51594b8e-0741-11eb-0c03-338b971e945f
# ╠═bf362532-0741-11eb-137e-67c2924f4f92
# ╠═db903b82-0741-11eb-2df7-0dcebbb93c01
# ╠═19bc9784-0742-11eb-18a4-edca29549ab4
# ╟─cac0339a-f8e4-11ea-0876-f91dc7db1454
# ╠═1058934e-f8e1-11ea-1896-35d75c4b6a8e
# ╟─c6cab78c-f836-11ea-2b5b-23f0006f8e93
# ╟─da4ac85c-f8e4-11ea-192e-03d5b6862d5a
# ╟─45e986f0-074c-11eb-0d49-dbce692ade1a
# ╠═665f1302-07b8-11eb-0848-8fb424cae65c
# ╠═a233a2e4-07c7-11eb-07d1-5561b4cc64b2
# ╠═39d465a6-f8e2-11ea-2d71-75cbe586428c
# ╠═2522ee1a-07ca-11eb-1d15-5f6b9e8f3aa6
# ╠═fe67814c-2299-11eb-2080-79f06d510355
# ╟─1c3b6c56-f8e2-11ea-3969-f59871025c5e
# ╟─bae1b4e6-07ae-11eb-1fb9-c56c0cf7298a
# ╠═e8716624-07b8-11eb-0219-81e10bd83100
# ╠═5176ef4e-07ba-11eb-09e7-9b6894fcc592
# ╟─63820b34-07af-11eb-21c3-75356e38ec43
# ╠═417c2dba-07b1-11eb-18cc-61d23caecbc6
# ╠═7e11e66e-07b3-11eb-1186-fba413764c97
# ╠═7c2633f6-2297-11eb-16f1-7b4f1b1a283c
# ╠═6a1c2364-2297-11eb-36b5-7f100ff97687
# ╠═7d28a38a-229a-11eb-2cec-2153fbbe4e63
# ╠═946112fc-229b-11eb-3027-bd10d12cb59b
# ╠═47b7c70c-229b-11eb-303a-41674232ea90
# ╠═32372d56-22c7-11eb-0923-4db4ed166e26
# ╠═da7f81ca-229a-11eb-138f-43dd54c21ad5
# ╠═fc1de00a-229b-11eb-22f3-13b7f379ed7d
# ╟─2315ce8c-fb57-11ea-31a4-eb68d1d096d3
# ╟─d9f6c1e8-fa89-11ea-05b3-07daaa06326f
# ╟─a1b031a4-fc08-11ea-2314-13e2e79c2ca3
# ╠═e6309486-fc08-11ea-3ade-c93a24b6866b
# ╠═9c3b0258-fc08-11ea-1922-fd1bec66b4ff
# ╠═975860d2-fc08-11ea-2def-ad1c45826911
# ╟─e7fb832e-04e0-11eb-388d-f7520b672a6e
# ╠═2b070720-04e2-11eb-3101-cbc40b17682b
# ╠═4b968440-04e2-11eb-35ce-5f050dcd4b6b
# ╠═6464f806-04e2-11eb-37cc-7363fa76c80e
# ╠═df5232e0-04e2-11eb-0c26-f19b9d9c7b14
# ╠═30b5258c-04e5-11eb-3cdc-11d36a2f35bc
# ╟─111aea9a-04e5-11eb-066d-4d245b04dc5e
# ╟─1e370abc-04e3-11eb-064a-21fbb9fafc7c
# ╠═1d60545a-04e4-11eb-1720-b5f04c5da239
# ╠═67b724d0-04e4-11eb-1254-f59bf78b1345
# ╠═86d10676-04e4-11eb-2121-2f0c8945cc1a
# ╠═979732dc-04e4-11eb-3803-7302dc91e777
# ╠═a33a18f2-04e4-11eb-0587-6117c036c641
# ╠═31071f6a-081a-11eb-2eda-f99e2073028c
# ╠═c3899100-04e4-11eb-37aa-394d3f60bd53
# ╠═d5c0b236-f98b-11ea-0444-6fe340642707
# ╠═28cb2832-0748-11eb-0c14-832f1df303ca
# ╟─79d08b08-fa8b-11ea-0b4b-01ceef03660a
# ╟─58c63ca8-fa8b-11ea-296d-0b7f613e8976
# ╟─96844d74-fa87-11ea-0966-89bbbc8a4353
# ╠═bfdd7d56-f98a-11ea-3064-a9033ba994a5
# ╠═eda9e2dc-f98b-11ea-3eca-539ac938107c
# ╠═f6e4754c-f98b-11ea-2a42-034d82f52ca1

max_stars_repo_name: schluppeck/hands-on-brain-data

max_issues_repo_name: schluppeck/hands-on-brain-data

max_forks_repo_name: schluppeck/hands-on-brain-data


text:
import LinearAlgebra.dot

function ismatrix(x::AbstractExpr)
    return (s = size(x); length(s) == 2 && s[1] > 1 && s[2] > 1)
end
ismatrix(::AbstractMatrix) = true
ismatrix(::Any) = false

# NOTE: Using asvec avoids broadcast-specific behaviors that we want to avoid, such
# as extending singleton dimensions. We need to ensure that the inputs have the same
# length, which broadcast will check for us if both inputs are vectors.
asvec(x) = convert(AbstractExpr, ismatrix(x) ? vec(x) : x)
_vecdot(x, y) = sum(broadcast(*, asvec(x), asvec(y)))

dot(x::AbstractExpr, y::AbstractExpr) = _vecdot(x, y)
dot(x::Value, y::AbstractExpr) = _vecdot(x, y)
dot(x::AbstractExpr, y::Value) = _vecdot(x, y)


text:
__precompile__(true)
module Pkg3

import Random
import REPL
using REPL.TerminalMenus

depots() = Base.DEPOT_PATH
logdir() = joinpath(depots()[1], "logs")
devdir() = get(ENV, "JULIA_PKG_DEVDIR", joinpath(homedir(), ".julia", "dev"))

have_warned_session = false
function print_first_command_header()
    global have_warned_session
    have_warned_session && return
    isinteractive() || return
    if !PKG3_IS_PRECOMPILED && !haskey(ENV, "JULIA_PKG3_DISABLE_PRECOMPILE_WARNING")
        @info """
        Pkg3 is running without precompile statements, first action will be slow.
        Rebuild julia with the environment variable `JULIA_PKG3_PRECOMPILE` set to enable precompilation of Pkg3.
        This message can be disabled by setting the env variable `JULIA_PKG3_DISABLE_PRECOMPILE_WARNING`.
        """
    end
    have_warned_session = true
end

# load snapshotted dependencies
include("../ext/TOML/src/TOML.jl")

include("GitTools.jl")
include("PlatformEngines.jl")
include("Types.jl")
include("Display.jl")
include("Pkg2/Pkg2.jl")
include("GraphType.jl")
include("Resolve.jl")
include("Operations.jl")
include("API.jl")
include("REPLMode.jl")

import .API: add, rm, up, test, gc, init, build, installed, pin, free, checkout, develop, generate
const update = up
import .REPLMode: @pkg_str
export @pkg_str

function __init__()
    if isdefined(Base, :active_repl)
        REPLMode.repl_init(Base.active_repl)
    else
        atreplinit() do repl
            if isinteractive() && repl isa REPL.LineEditREPL
                repl.interface = REPL.setup_interface(repl)
                REPLMode.repl_init(repl)
            end
        end
    end
end

using .Types
using UUIDs
import LibGit2
import Dates
# This crashes low memory systems and some of Julia's CI
# so keep it disabled by default for now.
if haskey(ENV, "JULIA_PKG3_PRECOMPILE")
    const PKG3_IS_PRECOMPILED = true
    include("precompile.jl")
else
    const PKG3_IS_PRECOMPILED = false
end

end # module

max_stars_repo_name: NymanLauri/julia

max_issues_repo_name: NymanLauri/julia

max_forks_repo_name: NymanLauri/julia


text:
Base.@kwdef struct Range
    lb::Float64 = 0.0 # to represent the normal lowerbound
    ub::Float64 = MAXNUMBER # to represent the normal upperbound
end
Range(v::Real) = Range(v, v)

Base.@kwdef struct Flexibility
    flexibility_level::Int = 0 # in level zero the nominal value is mandatory, in level k the nominal value can unsatisfyied if there was no feasbile solutions to the constraints of level 0 to k-1 that satisfy this level k contraint; a negative level means that the nominal value is optional, i.e. it must be statisfyied only it improves the economic value of the solution.
    fixed_price::Float64 = 0.0 # fixed_cost for not satisfying the nominal value, or fixed reward for satisfying it if it was optional
    unit_price::Float64 = 0.0 # while respecting the hard value, this represent a cost per unit away from the moninal value, or reward per unit away from the moninal value if the nominal value was optional
end

Base.@kwdef struct FlexibleRange
    soft_range::Range = Range()
    hard_range::Range = Range()
    lb_flex::Flexibility = Flexibility()
    ub_flex::Flexibility = Flexibility()
end

simple_flex_range(hard_lb::Float64,
                  soft_lb::Float64,
                  soft_ub::Float64,
                  hard_ub::Float64,
                  soft_violation_price::Float64) =
    FlexibleRange(Range(soft_lb, soft_ub),
                  Range(hard_lb, hard_ub),
                  Flexibility(0,0.0,soft_violation_price),
                  Flexibility(0,0.0,soft_violation_price))

Base.@kwdef mutable struct Location # Location where can be a Depot, Pickup, Delivery, Recharging, ..., or a combination of those services
    id::String = ""
    index::Int = -1 # used for matrices such as travel distance, travel time ...
    lat_y::Float64 = -1.0
    long_x::Float64 = -1.0
    opening_time_windows::Vector{Range} = [Range()]
    energy_fixed_cost::Float64 = 0.0 # an entry fee
    energy_unit_cost::Float64 = 0.0 # recharging cost per unit of energy
    energy_recharging_speeds::Vector{Float64} = Float64[] # if recharging in this location: the i-th speep is associted to the i-th energy interval defined for the vehicle
end

Base.@kwdef mutable struct LocationGroup # optionally defined to identify a set of locations with some commonalities, such as all possible pickups for a request
    id::String = ""
    location_ids::Vector{String} = String[]
end

Base.@kwdef mutable struct ProductCompatibilityClass # To define preceedence or conflict restriction between requested products
    id::String = ""
    conflict_compatib_class_ids::Vector{String} = String[]
    prohibited_predecessor_compatib_class_ids::Vector{String} = String[]
end

Base.@kwdef mutable struct ProductSharingClass # To define global availabitily restrictions for a product that is shared between different requests
    id::String = ""
    restricted_pickup_availabitilies::Bool = true
    restricted_delivery_capacities::Bool = true
    pickup_availabitilies_at_location_ids::Dict{String,Float64} = Dict{String,Float64}() # defined only if pickup locations have a restricted capacity; provides capcity for each pickup location where the product is avaiblable in restricted capacity
    delivery_capacities_at_location_ids::Dict{String,Float64} = Dict{String,Float64}() # defined only if delivery locations have a restricted capacity; provides capcity for each delivery location where the product can be delivered in restricted capacity
end

Base.@kwdef mutable struct ProductSpecificationClass # To define capacity consumption of a requested product
    id::String = ""
    capacity_consumptions::Dict{String,Tuple{Float64,Float64}} = Dict{String,Tuple{Float64,Float64}}() # to quantify the vehicle/compartment capacity that is used for accomodating lot-sizes of the request along several independant capacity measures whose string id key are in the dictionary: as weight, value, volume; for each such key, the capacity used is the float coef 2 * roundup(quantity / shipment_lot_size = float coef 1)
    property_requirements::Dict{String,Float64} = Dict{String,Float64}() # to check if the vehicle has the property of accomodating the request: yes if request requirement <= vehicle property capacity for each string id referenced requirement
end

Base.@kwdef mutable struct Request # can be
    # a shipment from a depot to a delivery location, or
    # a shipment from a pickup location to a depot, or
    # a delivery of a product that is shared by several requests, some of which are supplying the product while others are demanding the product, or
    # a pickup of a product that is shared by several requests, some of which are supplying the product while others are demanding the product, or
    # a shipment from a given pickup location to a given delivery location of a product that is specific to the request, or
    # a shipment from a given pickup location to any location of a group delivery locations of a product that is specific to the request, or
    # a shipment from any location of a group pickup locations to a given delivery location of a product that is specific to the request, or
    # a shipment from any location of a group pickup locations to any location of a group delivery locations of a product that is specific to the request.

    id::String = ""
    request_type::Int = 0 # 0 : pickup and delivery, 1: pickup only (delivery info is ignored) 2: delivery only (pickup info is ignored)
    product_compatibility_class_id::String = "default_id"
    product_sharing_class_id::String = "default_id"
    product_specification_class_id::String = "default_id"
    split_fulfillment::Bool = false # true if split delivery/pickup is allowed, default is false
    request_flexibility::Flexibility = Flexibility()
    precedence_status::Int = 0 # default = 0 = product predecessor restrictions;  1 = after all pickups, 2 =  after all deliveries.
    product_quantity_range::Range = Range(1) # of the request
    pickup_location_group_id::String = "default_id" # empty string for delivery-only requests. LocationGroup representing alternatives for pickup, otherwise.
    delivery_location_group_id::String = "default_id" # empty string for pickup-only requests. LocationGroup representing alternatives for delivery, otherwise.
    pickup_service_time::Float64 = 0.0 # used to measure pre-cleaning or loading time for instance
    delivery_service_time::Float64 = 0.0 # used to measure post-cleaning or unloading time for instance
    max_duration::Float64 = MAXNUMBER # to enforce a max duration between pickup and delivery
    duration_unit_cost::Float64 = 0.0 # to measure the cost of the time spent between pickup and delivery
    pickup_time_windows::Vector{FlexibleRange} = [FlexibleRange()]
    delivery_time_windows::Vector{FlexibleRange} = [FlexibleRange()]
end

Base.@kwdef mutable struct VehicleCharacteristics
    of_vehicle::Dict{String,Float64} = Dict{String,Float64}() # defined only if measured at the vehicle level; for string id key associated with properties or capacity measures that need to be checked on the vehicle, as for instance weight, value, volume
    of_compartments::Dict{String,Dict{String,Float64}} = Dict{String,Float64}() # defined only if measured at the compartment level; for string id key associated with properties or capacity measures that need to be checked on the vehicle, as for instance weight, value, volume, ... For each such property, the Dictionary specifies the capacity for each compartment id key.
end

Base.@kwdef mutable struct VehicleCategory
    id::String = ""
    capacities::VehicleCharacteristics = VehicleCharacteristics()
    properties::VehicleCharacteristics = VehicleCharacteristics()
    loading_option::Int = 0 # 0 = no restriction (=default), 1 = one request per compartment, 2 = removable compartment separation (note that product conflicts are measured within a compartment)
    energy_interval_lengths::Vector{Float64} = Float64[] # at index i, the length of the i-th energy interval. empty if no recharging.

end

Base.@kwdef mutable struct CostPeriod
    period::Range = Range()
    # CostSpecification
    fixed_cost::Float64 = 0.0
    travel_distance_unit_cost::Float64 = 0.0 # may depend on both driver and vehicle
    travel_time_unit_cost::Float64 = 0.0 # may depend on both driver and vehicle
    service_time_unit_cost::Float64 = 0.0
    waiting_time_unit_cost::Float64 = 0.0
end

Base.@kwdef mutable struct HomogeneousVehicleSet # vehicle type in optimization instance.
    id::String = ""
    route_mode::Int = 0 # 0 closed at departure and arrival, 1 open at arrival, 2 open at departure, 3 open at both depature and arrival
    vehicle_category_id::String = "default_id"
    cost_periods::Vector{CostPeriod} = [CostPeriod()]
    departure_location_group_id::String = "" # Vehicle routes start from one of the depot locations in the group
    arrival_location_group_id::String = "" # Vehicle routes end at one of the depot locations in the group
    work_periods::Vector{FlexibleRange} = [FlexibleRange()] # Define the work periods for which vehicles can be used, with some flexibility, does not have to be contiguous
    initial_energy_charge::Float64 = MAXNUMBER
    max_working_time::Float64 = MAXNUMBER # within each work period
    max_travel_distance::Float64 = MAXNUMBER # within each work period
    allow_shipment_over_multiple_work_periods::Bool = false # true if the vehicles do not need to complete all their requests by the end of each time period of the planning
    nb_of_vehicles_range::FlexibleRange = FlexibleRange()
end

Base.@kwdef mutable struct TravelSpecification
   id::String = ""
   travel_time_matrix::Array{Float64,2} = Array{Float64,2}(undef, 0, 0)
   travel_distance_matrix::Array{Float64,2} = Array{Float64,2}(undef, 0, 0)
   energy_consumption_matrix::Array{Float64,2} = Array{Float64,2}(undef, 0, 0)
end

Base.@kwdef mutable struct TravelPeriod
   period::Range = Range()
   travel_specification_id::String = ""
end

Base.@kwdef mutable struct RvrpInstance
    id::String = ""
    distance_mode::Int = 0 # 0 (default) - long and lat, 1 - x and y
    coordinate_mode::Int = 0 # 0 (default) user-defined, 1 - euc, 2 - manhattan, 3 - gps
    travel_specifications::Vector{TravelSpecification} = TravelSpecification[]
    travel_periods::Vector{TravelPeriod} = TravelPeriod[]
    locations::Vector{Location} = Location[]
    location_groups::Vector{LocationGroup} = LocationGroup[]
    product_compatibility_classes::Vector{ProductCompatibilityClass} = ProductCompatibilityClass[]
    product_sharing_classes::Vector{ProductSharingClass} = ProductSharingClass[]
    product_specification_classes::Vector{ProductSpecificationClass} = ProductSpecificationClass[]
    requests::Vector{Request} = Request[]
    vehicle_categories::Vector{VehicleCategory} = VehicleCategory[]
    vehicle_sets::Vector{HomogeneousVehicleSet} = HomogeneousVehicleSet[]
end


text:
@doc """
      Heun's Method (with iteration)
      Supported by:
      -PhD Kelvyn B. Sánchez kelvyn.baruc@gmail.com
 """ ->
function heun(f, tspan, y0; h = 1, AbsTol = 0.00001)

t = tspan[1]:h:tspan[end]
n = length(t)

      if t[n] < tspan[end]
            t = [t, tspan[end]]
            n += 1
      end

y1 = f(t[1],y0)
eq = length(y1) # Test the number of equations in f
y = zeros(Int(n),eq)
y[1,:] = y0
y1 = zeros(Int(n),eq)
y2 = zeros(Int(n),eq)
y3 = zeros(Int(n),eq)

	for i=1:n-1 # Solve for each x
            for j=1:eq # Solve for each equation in f
                  s1 = f(t[i],y[i,:])
                  y1[i,j] = y[i,j] + s1[j] * (t[i+1] - t[i])

                  s2 = f(t[i+1],y1[i,:])
                  y2[i,j] = y[i,j] + ((s1[j] + s2[j])/2) * (t[i+1] - t[i])

                  err = 100
                  while err >= AbsTol
                        s1 = f(t[i],y[i,:])
                        s2 = f(t[i+1],y2[i,:])
                        y3[i,j] = y[i,j] + ((s1[j] + s2[j]) / 2) * (t[i+1] - t[i])
                        err = abs((y3[i,j] - y2[i,j])/y3[i,j]) * 100
                        y2[i,j] = y3[i,j]
                        y[i+1,j] = y3[i,j]
                  end
            end
	end
	return t, y
end


text:
#!/usr/bin/env julia
using DocOpt, HDF5

doc = """
Pope (Pass one pulse examiner)
Usage:
  popeonce.jl <ljhpath> <preknowledge> <output>
  popeonce.jl --overwriteoutput <ljhpath> <preknowledge> <output>
  popeonce.jl -h | --help

Options:
  --overwriteoutput   Overwrite the output file if it already exists.
Help exposition:
  Where <ljhpath> points to a single ljh file, and <preknowledge> points to a valid HDF5 file in the pope preknowledge format. <output> points to a location where the output hdf5 file will be written.

  If you have succesfully run `Pkg.test("Pope")` from within Julia, this should work:
  ./popeonce.jl --overwriteoutput ~/.julia/v0.6/ReferenceMicrocalFiles/ljh/20150707_D_chan13.ljh ~/.julia/v0.6/Pope/test/preknowledge.h5 output.h5
"""

arguments = docopt(doc, version=v"0.0.1")
preknowledge_filename = expanduser(arguments["<preknowledge>"])
ljhpath = expanduser(arguments["<ljhpath>"])
output_filename = expanduser(arguments["<output>"])

using Pope.LJH

pkfile = h5open(preknowledge_filename,"r")
if !isfile(output_filename) || arguments["--overwriteoutput"]
  output_file = Pope.h5create(output_filename)
else
  println("ERROR: $output_filename already exists, cannot use it as output file")
  exit()
end

println("starting Pope analysis:")
@show ljhpath
@show preknowledge_filename
@show output_filename
println("")

Pope.init_for_zmqdatasink(Pope.ZMQ_PORT,verbose=true)
readers = Pope.Readers()
channels = Int[]
for name in names(pkfile)
  channel_number = startswith(name,"chan") ? parse(Int,name[5:end]) : parse(Int,name) # handle chan1 or 1
  ljh_filename = LJH.fnames(ljhpath,channel_number)
  if !isfile(ljh_filename)
    println("Channel $channel_number: exists in preknowledge file, but ljh does not exist")
    continue
  end
  analyzer = try
    analyzer = Pope.analyzer_from_preknowledge(pkfile[name])
  catch
    println("Channel $channel_number: failed to generate analyzer from preknowledge file")
    continue
  end
  product_writer = Pope.make_buffered_hdf5_and_zmq_multisink(output_file, channel_number, analyzer)
  reader = Pope.make_reader(ljh_filename, analyzer, product_writer,progress_meter=true)
  push!(readers, reader)
  push!(channels, channel_number)
end
println("Channels: ",repr(channels))
schedule(readers)
Pope.stop(readers) # stop readers immediatley, they will analyze all pulses in each ljh, but not look for new pulses
wait(readers)

println("done")


text:
import StaticArrays: SVector

function verified_prediction_lags(lags)
    ηs = sort(lags[lags .!= 0])

    if length(ηs) % 2 != 0
        throw(ArgumentError("Need an even number of lags (as many negative as positive prediction lags)"))
    end

    ηs_neg = ηs[ηs .< 0]
    ηs_pos = ηs[ηs .> 0]

    if length(ηs_neg) != length(ηs_pos)
        throw(ArgumentError("There must be as many negative as positive prediction lags. Got $ηs_neg and $ηs_pos"))
    end

    ηs_neg_sort = sort(abs.(ηs_neg))
    ηs_pos_sort = sort(ηs_pos)
    if !(all(ηs_neg_sort .== ηs_pos_sort))
        throw(ArgumentError("Negative amd positive prediction lags must be symmetric. Got $ηs_neg and $ηs_pos"))
    end

    return ηs
end

function verified_prediction_lags(lag::Int)
    # If only one η is provided (either positive or negative), just make a prediction
    # for that lag.
    ηs = [-abs(lag), abs(lag)]
end

"""
    PredictiveAsymmetryTest(predictive_test::CausalityTest)

The parameters for a predictive asymmetry causality test [1]. For a
normalised version, see [`NormalisedPredictiveAsymmetryTest`](@ref).

## Mandatory keywords

- **`predictive_test`**: An instance of a predictive causality test that explicitly
    uses prediction lags (e.g. [`VisitationFrequencyTest`](@ref) or
    [`TransferOperatorGridTest`](@ref)).

## About the prediction lags

The prediction lags in the predictive causality test must consist of `n` negative
integers and `n` positive integers that are symmetric around zero.

In other words, negative lags  must exactly match the positive lags but with opposite
sign. The zero lag can be included, but will be ignored, so it is possible to give
ranges too.

## Examples

```julia
bin = RectangularBinning(5) # divide each axis into 5 equal-length intervals
ηs = [-10:-1; 1:10] # exclude the zero lag (it is irrelevant for the asymmetry)

test_visitfreq = VisitationFrequencyTest(ηs = ηs)
test_transferoperator = TransferOperatorGridTest(ηs = ηs)

# Note that `predictive_test` is a *mandatory* keyword.
PredictiveAsymmetryTest(predictive_test = test_visitfreq)
PredictiveAsymmetryTest(predictive_test = test_transferoperator)
```

## References

1. Diego, David, Kristian Agasøster Haaga, Jo Brendryen, and Bjarte Hannisdal.
    A simple test for causal asymmetry in complex systems. In prep.
"""
Base.@kwdef mutable struct PredictiveAsymmetryTest{TEST, N} <: AbstractPredictiveAsymmetryTest where TEST
    predictive_test::TEST

    function PredictiveAsymmetryTest(test::T) where {T <: TransferEntropyCausalityTest}
        # Check that prediction lags are okay
        verified_prediction_lags(test.ηs)
        N = length(test.ηs[test.ηs .> 0])
        return new{T, N}(test)
    end
end

function causalbalance(lags, tes)
    posinds = findall(lags .> 0)
    neginds = findall(lags .< 0)
    sum(tes[posinds]) - sum(tes[neginds])
end

function update_ηs(test::VisitationFrequencyTest)

    VisitationFrequencyTest(
        k = test.k, l = test.l, m = test.m, n = test.n,
        τ = test.τ,
        binning_summary_statistic = test.binning_summary_statistic,
        binning = test.binning,
        ηs = verified_prediction_lags(test.ηs))
end

function update_ηs(test::TransferOperatorGridTest)

    TransferOperatorGridTest(
        k = test.k, l = test.l, m = test.m, n = test.n,
        τ = test.τ,
        binning_summary_statistic = test.binning_summary_statistic,
        binning = test.binning,
        ηs = verified_prediction_lags(test.ηs))
end

function update_ηs(test::NearestNeighbourMITest{N}) where N

    NearestNeighbourMITest(
        k = test.k, l = test.l, m = test.m, n = test.n,
        τ = test.τ,
        ηs = verified_prediction_lags(test.ηs))
end

function return_predictive_asymmetry(ηs, As, N)
    #T = typeof(As[1])
    #SVector{N, T}(As)
    As
end

return_predictive_asymmetry(η::Int, As, N) = As[1]

"""
    predictive_asymmetry(source, target, test::PredictiveAsymmetryTest{T, N}) where {T <: TransferEntropyCausalityTest, N}

Compute the predictive asymmetry from `source` to `target` using the provided predictive test `p`.
The test can be a [`VisitationFrequencyTest`](@ref), [`TransferOperatorGridTest`](@ref) or
[`NearestNeighbourMITest`](@ref) and has to be defined for prediction lags symmetrically
around zero.

## Example

We'll use the [`nonlinear3d`](@ref) system, which consists of three variables
and is unidirectionally coupled ``x \\to y \\to z``.

```julia
# Example data
sys = nonlinear3d()
npts = 300
orbit = trajectory(sys, npts, Ttr = 300);
x, y, z = columns(orbit);

# Test setup
η_max = 10
ηs = [-η_max:-1; 1:η_max]
bin = RectangularBinning(floor(Int, npts^(1/4)))
test = VisitationFrequencyTest(ηs = ηs, binning = bin)
pa_test = PredictiveAsymmetryTest(test)

# Analysis
pas_xy = causality(x, y, pa_test)
pas_yx = causality(y, x, pa_test)
pas_yz = causality(y, z, pa_test)
pas_zy = causality(z, y, pa_test)

[pas_xy pas_yx pas_yz pas_zy]
```

If the test correctly picked up the correct directionality, then the values of
of `asymmetries` corresponding to the causal interactions (1st and 3rd column)
should be mostly positive, and those corresponding to non-causal interactions
should be mostly negative.
"""
function predictive_asymmetry(source, target,
        p::PredictiveAsymmetryTest{T, N}) where {T <: TransferEntropyCausalityTest, N}

    # Update the test parameters so that we have symmetric prediction lags
    test = update_ηs(p.predictive_test)

    # Predictions from source to target.
    preds = causality(source, target, test)

    # The number of predictive asymmetries is half the number of prediction lags,
    # which is encoded in the type parameter `N`
    ηs = test.ηs
    As = zeros(Float64, N)

    for (i, η) in enumerate(ηs[ηs .> 0])
        lag_idxs = findfirst(ηs .== -η):findfirst(ηs .== η)
        As[i] = causalbalance(ηs[lag_idxs], preds[lag_idxs])
    end

    return return_predictive_asymmetry(p.predictive_test.ηs, As, N)
end

function predictive_asymmetry(source, target, cond,
    p::PredictiveAsymmetryTest{T, N}) where {T <: TransferEntropyCausalityTest, N}

    # Update the test parameters so that we have symmetric prediction lags
    test = update_ηs(p.predictive_test)

    # Predictions from source to target.
    preds = causality(source, target, cond, test)

    # The number of predictive asymmetries is half the number of prediction lags,
    # which is encoded in the type parameter `N`
    ηs = test.ηs
    As = zeros(Float64, N)

    for (i, η) in enumerate(ηs[ηs .> 0])
        lag_idxs = findfirst(ηs .== -η):findfirst(ηs .== η)
        As[i] = causalbalance(ηs[lag_idxs], preds[lag_idxs])
    end

    return return_predictive_asymmetry(p.predictive_test.ηs, As, N)
end

function causality(source::AbstractVector{T}, target::AbstractVector{T}, p::PredictiveAsymmetryTest{CT}) where {T<:Real, CT}
    predictive_asymmetry(source, target, p)
end

function causality(source::AbstractVector{T}, target::AbstractVector{T}, cond::AbstractVector{T},
        p::PredictiveAsymmetryTest{CT}) where {T<:Real, CT}
    predictive_asymmetry(source, target, cond, p)
end

# there is no need to define custom causality method for a uncertain data
# as long as the predictive tests that it uses supports them

export
    PredictiveAsymmetryTest,
    predictive_asymmetry,
    update_ηs,
    verified_prediction_lags


text:
using ApproxFun, ODE, Plots

import ODE: ode23, ode45, ode78, ode4, ode4ms, ode23s, ode4s

for ode in (:ode23, :ode45, :ode78, :ode4, :ode4ms, :ode23s, :ode4s)
    @eval function $ode(fn, t0, y0, d::Space; reltol=eps(), kwargs...)
        t0added=false
        for logn = 4:20
            t = points(d, 2^logn)
            t0 ∉ t && (unshift!(t,t0); t0added=true)
            t, u = $ode(fn, y0, t; points=:specified, reltol=reltol, kwargs...)
            cf = t0added ? Fun(ApproxFun.transform(d,u[2:end]), d) : Fun(ApproxFun.transform(d,u), d)

            maxabsc=maxabs(cf.coefficients)
            if maxabsc==0 && maxabsfr==0
                return(zeros(d))
            end

            if length(cf) > 8 && maxabs(cf.coefficients[end-8:end]) < reltol
                return chop!(cf,reltol/10)
            end
        end
        warn("Maximum length "*string(2^20+1)*" reached")

        t = points(d, 2^21)
        t0 ∉ t && (unshift!(t,t0); t0added=true)
        t, u = $ode(fn, y0, t; points=:specified, kwargs...)

        t0added ? Fun(ApproxFun.transform(d,u[2:end]), d) : Fun(ApproxFun.transform(d,u), d)
    end
end

sp = Chebyshev(Interval(2,0))

u0 = ode45((t,u)->-u^2/(2-.1^2*u),2.0,-3.0,sp;reltol=eps(Float32),abstol=realmin())

pyplot()

plot(u0)

norm(u0'+u0^2/(2-.1^2*u0))

N=u->[u(2)+3.0;u'+u^2/(2-.1^2*u)]

u = newton(N,u0)

norm(u'+u^2/(2-.1^2*u))


text:
# (4) Retail sales, excl. food services
retail08 = get_data(api_fred, "RSXFS"; observation_start = "2008-06-01", observation_end = "2011-03-01")
retail20 = get_data(api_fred, "RSXFS"; observation_start = "2020-02-01")

idx08 = retail08.data[1, :value]
retail08.data[:, :value] = retail08.data[:, :value] ./ idx08 * 100

idx20 = retail20.data[1, :value]
retail20.data[:, :value] = retail20.data[:, :value] ./ idx20 * 100

len08 = length(retail08.data[:, :value])
len20 = length(retail20.data[:, :value])

plot_retail08 = layer(retail08.data, x = 1:len08, y = :value, Geom.line, Theme(default_color=colorant"red"));
plot_retail20 = layer(retail20.data, x = 1:len20, y = :value, Geom.line, Theme(default_color=colorant"deepskyblue"));
plt_retail = plot(plot_retail08, plot_retail20,
        Guide.manual_color_key("Retail sales, excluding food services",
                                ["2008 Recession",
                                 "2020 Pandemic ($(maximum(retail20.data[:date])): $(round(last(retail20.data[:value]), digits = 1)))"],
                                ["red", "deepskyblue"],
                                pos = [0.13w,0.38h]
                               ),
         Guide.xlabel("Month(s) since the start of the recession"), Guide.ylabel("Index (100: pre-decline peak)"),
         Coord.cartesian(xmin=1, xmax=len08)
    )
pngout = SVG("plt_retail.svg", 5inch, 4inch)
draw(pngout, plt_retail)


text:
using OMEinsum: tsetdiff, tunique

@testset "utils" begin
    @test tsetdiff((1,2,3), (2,)) == [1,3]
    @test tunique((1,2,3,3,)) == [1,2,3]
    @test asarray(3) isa Array
    @test asarray(3, randn(3,3)) isa Array
    @test asarray(randn(3,3)) isa Array
end

@testset "allunique" begin
    @test allunique(())
    @test allunique(('i',))
    @test allunique((1,2))
    @test !allunique((1,2,1))
end

@testset "conditioned_permutedims" begin
    a = randn(100, 100)
    OMEinsum.conditioned_permutedims(a, [1,2])
    @test (@allocated OMEinsum.conditioned_permutedims(a, (1,2))) < 100
end

@testset "align_types" begin
    a = randn(100, 100)
    b = randn(ComplexF64, 100, 100)
    res = OMEinsum.align_eltypes(a, b)
    @test res[2] === b
    @test res == (ComplexF64.(a), b)
end


text:
using BinaryProvider

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libnnpack"], :libnnpack),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaPackaging/Yggdrasil/releases/download/NNPACK-v2018.06.22-0"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc) => ("$bin_prefix/NNPACK.v2018.6.22.aarch64-linux-gnu.tar.gz", "e0c6e21ba4c47acfd5a3d3e3510e8786474080f654338f4583b88860296c1437"),
    Linux(:i686, libc=:glibc) => ("$bin_prefix/NNPACK.v2018.6.22.i686-linux-gnu.tar.gz", "e9b6685001bc5a5d17acef15f3f6ffeb7beb6081926300f23ed4a442beac71ca"),
    Linux(:i686, libc=:musl) => ("$bin_prefix/NNPACK.v2018.6.22.i686-linux-musl.tar.gz", "36c1d3c30b3bc3e0b34f215945bb46319f88e28f011fc758f21ba888b1fd9e25"),
    MacOS(:x86_64) => ("$bin_prefix/NNPACK.v2018.6.22.x86_64-apple-darwin14.tar.gz", "b30046223a11470b15a2ceb0d0df6f7d8a43260fe52f4a2f8ebe5f0b2df822ca"),
    Linux(:x86_64, libc=:glibc) => ("$bin_prefix/NNPACK.v2018.6.22.x86_64-linux-gnu.tar.gz", "150d5b6ca81fa72bfdc8bbda2428f0d3483fd11a5813724646c6d6c6a7ef969f"),
    Linux(:x86_64, libc=:musl) => ("$bin_prefix/NNPACK.v2018.6.22.x86_64-linux-musl.tar.gz", "d961a104f814ec5b356519a82746a70a1df193ae37fc8130f38ffb61336def16"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    @warn "Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by NNPACK!
          You will only be able to use only the default NNlib backend."
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
# Download and install binaries
use_nnpack = get(ENV, "NNLIB_USE_NNPACK", "false") == "true"
os_support = Sys.islinux() || Sys.isapple()
if use_nnpack && os_support
  if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    install(dl_info...; prefix=prefix, force=true, verbose=verbose)
  end
  # Write out a deps.jl file that will contain mappings for our products
  write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)
else
  open(joinpath(@__DIR__, "deps.jl"), "w") do io
    write(io, "check_deps() = false")
  end
end


text:
using DelimitedFiles
using Logging
using Plots
#using XfromProjections.curve_utils
#using XfromProjections.snake_forward
#using XfromProjections
using Printf
using LinearAlgebra
include("./utils.jl")
include("../../phantoms/sperm_phantom.jl")
include("../snake_forward.jl")
include("../snake.jl")
include("../curve_utils.jl")

cwd = @__DIR__
images, tracks = get_sperm_phantom(11,grid_size=0.1)

#Radius function
r(s) = 1.0
bins = collect(-38.0:0.125:38.0)
frame_nr = 1
nangles = 1
angles = [π/2]

#determine outline from skeleton
outline, normals = get_outline(tracks[1], r)
sinogram = parallel_forward(outline, angles, bins)

num_points = 30

template = get_straight_template(sinogram[:,1], r, [0.0 0.0], angles[1], num_points,bins)

plt = plot(template[:,1], template[:,2], aspect_ratio=:equal, label="template", framestyle=:none, color=:black,  size=(400,400), legend=false)
plot!(outline[:,1], outline[:,2], label = "target", color=:green)

n = 5000
recon = deepcopy(template)

w_u = ones(num_points+2)
w_u[1] = 0.0
w_u[2] = 0.0
w_l = ones(num_points+2)
w_l[1] = 0.0
w_l[2] = 0.0

recon = recon2d_tail(recon,r,angles,bins,sinogram,1, 0.0, 0.1, 1, w_u, w_l, plot=true)

plot!(recon[:,1],recon[:,2], label="reconstruction", color=:orange, linewidth=2)

outline, normals = get_outline(recon, r)
s = parallel_forward(outline, angles, bins)

y_vals = cat(bins,s,dims=2)
y_vals = rotate_points(translate_points(y_vals, [0.0 20.0]), angles[1])
plot_vals = y_vals[10:end-250,:]
plot!(plot_vals[:,1], plot_vals[:,2],label="FP", fill = (0, 0.2, :orange), color=:white)

y_vals = cat(bins,sinogram,dims=2)
y_vals = rotate_points(translate_points(y_vals, [0.0 20.0]), angles[1])
plot_vals = y_vals[10:end-250,:]
plot!(plot_vals[:,1], plot_vals[:,2],label="projection", color=:green)

cwd = @__DIR__
savepath = normpath(joinpath(@__DIR__, "result"))
!isdir(savepath) && mkdir(savepath)
cd(savepath)
savefig("stalemate")

cd(cwd)

plot!()

max_stars_repo_name: QPImaging/xfromProjections

max_issues_repo_name: QPImaging/xfromProjections

max_forks_repo_name: QPImaging/xfromProjections


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/UMAT.jl/blob/master/LICENSE

using UMAT
using Test
using Suppressor

@testset "UMAT.jl" begin
    @testset "Binary dependencies" begin
        include("test_binary_dependencies.jl")
    end
    @testset "libelastic.so" begin
        # include("test_elastic.jl")
    end
    @testset "UmatMaterial" begin
        # include("test_umatmaterial.jl")
    end
    @testset "DruckerPragerMaterial" begin
        include("test_druckerpragermaterial.jl")
    end
    @testset "GursonMaterial" begin
        # include("test_gurson.jl")
    end
end


text:
using CUDA
using KernelAbstractions
using LinearAlgebra
using SparseArrays
using Random
using Printf
using TimerOutputs

import ConjGrad
import IterativeSolvers
import Krylov
import KrylovKit
import Solvent

struct BenchmarkSetup{AT, T}
    A::AT
    b::Vector{T}
    xexact::Vector{T}
    xinit::Vector{T}

    function BenchmarkSetup(A::AbstractArray{<:Number})
        T = eltype(A)
        _, ncols = size(A)
        xexact = rand(T, ncols)
        xinit = rand(T, ncols)
        return new{typeof(A), T}(A, A * xexact, xexact, xinit)
    end
end

mutable struct StatTracker{T<:Number}
    min::T
    max::T
    tot::T

    function StatTracker(value::T) where T<:Number
        return new{T}(value, value, value)
    end
end

function addstat!(st::StatTracker{T}, value::T) where T<:Number
    st.tot += value
    if value < st.min
        st.min = value
    elseif value > st.max
        st.max = value
    end
end

mutable struct BenchmarkTracker
    runs::Int
    time::StatTracker
    alloc::StatTracker
    iter::StatTracker
    res::StatTracker
    err::StatTracker

    function BenchmarkTracker(time::AbstractFloat, allocs::AbstractFloat,
            iters::Integer, res::AbstractFloat, err::AbstractFloat)
        return new(1, StatTracker(time), StatTracker(allocs),
            StatTracker(iters), StatTracker(res), StatTracker(err))
    end
end

function addbenchmark!(bt::BenchmarkTracker, time::AbstractFloat,
        allocs::AbstractFloat, iters::Integer, res::AbstractFloat,
        err::AbstractFloat)
    bt.runs += 1
    addstat!(bt.time, time)
    addstat!(bt.alloc, allocs)
    addstat!(bt.iter, iters)
    addstat!(bt.res, res)
    addstat!(bt.err, err)
end

mutable struct Benchmarker
    storage::Dict{String, BenchmarkTracker}
    enabled::Bool

    function Benchmarker(enabled::Bool = true)
        return new(Dict{String, BenchmarkTracker}(), enabled)
    end
end

function enable!(b::Benchmarker)
    b.enabled = true
end

function disable!(b::Benchmarker)
    b.enabled = false
end

function recordbenchmark!(b::Benchmarker, name::String, time::AbstractFloat,
        allocs::AbstractFloat, iters::Integer, res::AbstractFloat,
        err::AbstractFloat)
    if b.enabled
        if haskey(b.storage, name)
            addbenchmark!(b.storage[name], time, allocs, iters, res, err)
        else
            b.storage[name] = BenchmarkTracker(time, allocs, iters, res, err)
        end
    end
end

function printbenchmarks(b::Benchmarker)
    storage = b.storage
    maxnamelen = maximum([length(name) for name in keys(storage)])
    print(
        rpad("Benchmark Name", maxnamelen),
        "| Runs |       Time (ms)      | Allocations (MiB) | Iterations  |",
        " Relative Residual |  Absolute Error  \n",
        " " ^ maxnamelen,
        "|      |   Min    Max    Avg  |  Min   Max   Avg  | Min Max Avg |",
        "  Min   Max   Avg  |  Min   Max   Avg \n",
        "-" ^ (maxnamelen + 103),
        "\n")
    for name in sort!(collect(keys(storage)))
        bt = storage[name]
        print(rpad(name, maxnamelen + 3))
        @printf("%-3d   %6.1f %6.1f %6.1f   %5.1f %5.1f %5.1f   %3d %3d %3d   ",
            bt.runs,
            bt.time.min, bt.time.max, bt.time.tot / bt.runs,
            bt.alloc.min, bt.alloc.max, bt.alloc.tot / bt.runs,
            bt.iter.min, bt.iter.max, bt.iter.tot / bt.runs)
        @printf("%.0e %.0e %.0e   %.0e %.0e %.0e\n",
            bt.res.min, bt.res.max, bt.res.tot / bt.runs,
            bt.err.min, bt.err.max, bt.err.tot / bt.runs)
    end
end

# Constructs an m×n random sparse matrix with a relatively low condition number
# and a density no less than d.
function simplesparsesetup(::Type{T}, m::Integer, α::Number=0.01,
        d::AbstractFloat=0.05, n::Integer=m) where T
    return BenchmarkSetup(I + T(α) * sprandn(T, m, n, d))
end

# Constructs a m×n random sparse matrix with reciprocal condition number rc and
# a density roughly equal to (though slightly above) d. Based on the code from
# https://searchcode.com/file/124101202/scripts/sparse/private/__sprand__.m/
function sparsesetup(::Type{T}, m::Integer, rc::Real=0.9,
        d::AbstractFloat=0.05, n::Integer=m) where T
    if rc < 0 || rc > 1
        error("reciprocal condition number lies outside the interval [0, 1]")
    end

    singularvalues = rand(min(m, n)) .* (1 - rc) .+ rc
    singularvalues[1] = 1
    singularvalues[end] = rc
    S = SparseMatrixCSC{T, Int}(spdiagm(m, n, 0=>singularvalues))

    U = SparseMatrixCSC{T, Int}(sparse(I, m, m))
    V = SparseMatrixCSC{T, Int}(sparse(I, n, n))
    k = m * n * d
    while nnz(S) < k
        if m > 1
            angle = rand() * 2 * pi
            c = cos(angle)
            s = sin(angle)
            i, j = randperm(m)[1:2]
            U[i, i] = c
            U[i, j] = -s
            U[j, i] = s
            U[j, j] = c
            S = U * S
            U[i, i] = U[j, j] = 1
            U[i, j] = U[j, i] = 0
            dropzeros!(U)
        end
        if n > 1
            angle = rand() * 2 * pi
            c = cos(angle)
            s = sin(angle)
            i, j = randperm(n)[1:2]
            V[i, i] = c
            V[i, j] = s
            V[j, i] = -s
            V[j, j] = c
            S = S * V
            V[i, i] = V[j, j] = 1
            V[i, j] = V[j, i] = 0
            dropzeros!(V)
        end
    end

    return BenchmarkSetup(S)
end

# Constructs the matrix approximation of -∇² for an lx×ly box of nx×ny points
# with Dirichlet boundary conditions. Based on the code from
# http://math.mit.edu/~stevenj/18.303/lecture-10.html
function laplacesetup(::Type{T}, nx::Integer, ny::Integer=nx, lx::Integer=1,
        ly::Integer=lx) where T
    dx = T(lx) / T(nx + 1)
    dy = T(ly) / T(ny + 1)
    Dx = [[one(T) spzeros(T, 1, nx - 1)]; spdiagm(1=>ones(T, nx - 1)) - I] / dx
    Dy = [[one(T) spzeros(T, 1, ny - 1)]; spdiagm(1=>ones(T, ny - 1)) - I] / dy
    Ax = Dx' * Dx
    Ay = Dy' * Dy
    A = kron(sparse(I, ny, ny), Ax) + kron(Ay, sparse(I, nx, nx))
    return BenchmarkSetup(A)
end

function benchmarkcg!(benchmarker::Benchmarker, bs::BenchmarkSetup,
        setupname::String, maxiters::Integer, tol::AbstractFloat)
    A = bs.A
    b = bs.b
    xexact = bs.xexact
    xinit = bs.xinit
    algnames = ["ConjGrad", "IterativeSolvers", "Krylov", "KrylovKit",
        "Solvent"]
    times = Array{Float64, 1}(undef, 5)
    allocs = Array{Int, 1}(undef, 5)
    iters = Array{Int, 1}(undef, 5)
    xfinals = [copy(xinit) for _ in 1:5]
    result, times[1], allocs[1], _, _ = @timed begin
        ConjGrad.cg!((x) -> (A * x), b, xfinals[1]; tol = Float64(tol),
            maxIter = maxiters)
    end
    iters[1] = result[2]
    result, times[2], allocs[2], _, _ = @timed begin
        IterativeSolvers.cg!(xfinals[2], A, b, tol = tol, maxiter = maxiters,
            log = true)
    end
    iters[2] = result[2].iters
    result, times[3], allocs[3], _, _ = @timed begin
        Krylov.cg(A, b; atol = tol, rtol = tol, itmax = maxiters)
    end
    iters[3] = length(result[2].residuals) - 1 # 1st residual found before loop.
    xfinals[3] = result[1]
    result, times[4], allocs[4], _, _ = @timed begin
        KrylovKit.linsolve(A, b, xinit; isposdef = true, krylovdim = 1,
            maxiter = maxiters, atol = tol, rtol = tol)
    end
    iters[4] = result[2].numops - 1 # Linear operator applied once before loop.
    xfinals[4] = result[1]
    iters[5], times[5], allocs[5], _, _ = @timed begin
        linearsolver = Solvent.LinearSolver(
            (y, x) -> (y .= A * x),
            Solvent.ConjugateGradientMethod(M = maxiters),
            xinit;
            pc_alg = Solvent.Identity(pc_side=Solvent.PCleft()),
            rtol = tol,
            atol = tol,
        )
        Solvent.linearsolve!(linearsolver, xfinals[5], b)
    end
    prefix = join(["CG", setupname, eltype(xexact), size(A)[1], maxiters], '-')
    for i in 1:5
        recordbenchmark!(benchmarker, string(prefix, '-', algnames[i]),
            times[i] * 1000, allocs[i] / 2^20, iters[i],
            norm(A * xfinals[i] - b) / norm(A * xinit - b),
            norm(xfinals[i] - xexact))
    end
end

function benchmarkgmres!(benchmarker::Benchmarker, bs::BenchmarkSetup,
        setupname::String, dims::Integer, maxiters::Integer, tol::AbstractFloat)
    A = bs.A
    b = bs.b
    xexact = bs.xexact
    xinit = bs.xinit
    algnames = ["IterativeSolvers", "Krylov", "KrylovKit", "Solvent"]
    times = Array{Float64, 1}(undef, 4)
    allocs = Array{Int, 1}(undef, 4)
    iters = Array{Int, 1}(undef, 4)
    xfinals = [copy(xinit) for _ in 1:4]
    result, times[1], allocs[1], _, _ = @timed begin
        IterativeSolvers.gmres!(xfinals[1], A, b, tol = tol, restart = dims,
            maxiter = dims * maxiters, log = true)
    end
    iters[1] = result[2].iters
    result, times[2], allocs[2], _, _ = @timed begin
        Krylov.dqgmres(A, b; atol = tol, rtol = tol, itmax = dims * maxiters)
    end
    iters[2] = length(result[2].residuals) - 1 # 1st residual found before loop.
    xfinals[2] = result[1]
    result, times[3], allocs[3], _, _ = @timed begin
        KrylovKit.linsolve(A, b, xinit; isposdef = false, krylovdim = dims,
            maxiter = maxiters, atol = tol, rtol = tol)
    end
    iters[3] = result[2].numops - 2 # Linear operator applied twice before loop.
    xfinals[3] = result[1]
    iters[4], times[4], allocs[4], _, _ = @timed begin
        linearsolver = Solvent.LinearSolver(
            (y, x) -> (y .= A * x),
            Solvent.GeneralizedMinimalResidualMethod(M = dims, K = maxiters),
            xinit;
            pc_alg = Solvent.Identity(pc_side=Solvent.PCright()),
            rtol = tol,
            atol = tol,
        )
        Solvent.linearsolve!(linearsolver, xfinals[4], b)
    end
    prefix = join(["GMRES", setupname, eltype(xexact), size(A)[1], dims,
        maxiters], '-')
    for i in 1:4
        recordbenchmark!(benchmarker, string(prefix, '-', algnames[i]),
            times[i] * 1000, allocs[i] / 2^20, iters[i],
            norm(A * xfinals[i] - b) / norm(A * xinit - b),
            norm(xfinals[i] - xexact))
    end
end

const b = Benchmarker()
for T in [Float32, Float64]
    tol = sqrt(eps(T))
    for i in 0:100
        i == 0 && disable!(b)
        benchmarkcg!(b, laplacesetup(T, 100), "Laplace", 5000, tol)
        benchmarkgmres!(b, laplacesetup(T, 100), "Laplace", 500, 10, tol)
        benchmarkgmres!(b, simplesparsesetup(T, 10000), "SimpleSparse", 14, 1, tol)
        benchmarkgmres!(b, simplesparsesetup(T, 10000), "SimpleSparse", 12, 3, tol)
        benchmarkgmres!(b, simplesparsesetup(T, 10000), "SimpleSparse", 11, 5, tol)
        benchmarkgmres!(b, simplesparsesetup(T, 10000), "SimpleSparse", 10, 8, tol)
        # Solvent does not converge for SimpleSparse with M = 9, even if K = 10000.
        benchmarkgmres!(b, simplesparsesetup(T, 10123), "SimpleSparse", 10, 90, tol)
        # Packages stop converging for Sparse 1, 2, 5, and Inf when M is lowered.
        benchmarkgmres!(b, sparsesetup(T, 500, 1), "Sparse1", 460, 2, tol)
        benchmarkgmres!(b, sparsesetup(T, 500, 0.5), "Sparse2", 475, 2, tol)
        benchmarkgmres!(b, sparsesetup(T, 500, 0.2), "Sparse5", 490, 2, tol)
        benchmarkgmres!(b, sparsesetup(T, 500, 0), "SparseInf", 499, 2, tol)
        i == 0 && enable!(b)
    end
end
printbenchmarks(b)

max_stars_repo_name: CliMA/SystemSolvers

max_issues_repo_name: CliMA/SystemSolvers


text:
"""
    Execution time comparison.

    As for the first experimental scenario, we fixed
    random thresholds for both nodes and hyperedges.

    For each heuristic, we consider the average time (in seconds)
    required to complete the task.
"""
using Pkg, Distributed
#addprocs(4)
Pkg.activate(".")
@everywhere using Distributed, Pkg
@everywhere Pkg.activate(".")
using LTMSim, SimpleHypergraphs, Random, Serialization
@everywhere using LTMSim, SimpleHypergraphs, Random, Serialization

#
#
#
project_path = dirname(pathof(LTMSim))
data_path = joinpath(project_path, "..", "data", "hgs")
res_path = joinpath(project_path, "..", "res", "journal", "benchmark_randV_randE.data")

hg_files = readdir(data_path)
hgs = [prune_hypergraph!(hg_load(joinpath(data_path, hg_file))) for hg_file in hg_files]

runs = 50
data = Dict{String, Vector{Vector{Float64}}}()

push!(data, "BinarySearch(H)" => Vector{Float64}[])
push!(data, "Greedy(H)" => Vector{Float64}[])
push!(data, "Greedy([H]₂)" => Vector{Float64}[])
push!(data, "SubTSS(H)" => Vector{Float64}[])

push!(data, "BinarySearch(H)-noOpt" => Vector{Float64}[])
push!(data, "Greedy(H)-noOpt" => Vector{Float64}[])
push!(data, "Greedy([H]₂)-noOpt" => Vector{Float64}[])
push!(data, "SubTSS(H)-noOpt" => Vector{Float64}[])

for index in 1:length(hgs)
    println("Index=$index, $(hg_files[index])")

    h = hgs[index]

    results = @distributed (append!) for run=1:runs
		println("run=$run at proc $(myid())")
        Random.seed!(run)

        metaV = randMetaV(h)
        metaE = randMetaE(h)

        r1 = @belapsed greedy_tss_2section($h, $metaV, $metaE; opt=true)
        r2 = @belapsed bisect($h, $metaV, $metaE; opt=true)
        r3 = @belapsed greedy_tss($h, $metaV, $metaE; opt=true)
        r4 = @belapsed sub_tss($h, $metaV, $metaE; opt=true)

        #no-opt
        r5 = @belapsed greedy_tss_2section($h, $metaV, $metaE; opt=false)
        r6 = @belapsed bisect($h, $metaV, $metaE; opt=false)
        r7 = @belapsed greedy_tss($h, $metaV, $metaE; opt=false)
        r8 = @belapsed sub_tss($h, $metaV, $metaE; opt=false)

        [(r1,r2,r3,r4,r5,r6,r7,r8)]
    end

    push!(data["Greedy([H]₂)"], [r[1] for r in results])
    push!(data["BinarySearch(H)"], [r[2] for r in results])
    push!(data["Greedy(H)"], [r[3] for r in results])
    push!(data["SubTSS(H)"], [r[4] for r in results])

    #no-opt
    push!(data["Greedy([H]₂)-noOpt"], [r[5] for r in results])
    push!(data["BinarySearch(H)-noOpt"], [r[6] for r in results])
    push!(data["Greedy(H)-noOpt"], [r[7] for r in results])
    push!(data["SubTSS(H)-noOpt"], [r[8] for r in results])
end

serialize(res_path, data)


text:
export Polygonic
export Polygon, Rectangle

export GeoParametric
export Sphere

export isIntersect, isInside
export show

include("polygonic.jl")
include("geoparametric.jl")


text:
using StableRNGs
@testset "Manifolds" begin
    rng = StableRNG(213)

    # Test case: find eigenbasis for first two eigenvalues of a symmetric matrix by minimizing the Rayleigh quotient under orthogonality constraints
    n = 4
    m = 2
    A = Diagonal(range(1, stop=2, length=n))
    fmanif(x) = real(dot(x,A*x)/2)
    gmanif(x) = A*x
    gmanif!(stor,x) = copyto!(stor,gmanif(x))
    # A[2,2] /= 10 #optional: reduce the gap to make the problem artificially harder
    x0 = randn(rng,n,m)+im*randn(rng,n,m)

    @testset "Stiefel, retraction = $retraction" for retraction in [:SVD,:CholQR]
        manif = Optim.Stiefel(retraction)

        # AcceleratedGradientDescent should be compatible also, but I haven't been able to make it converge
        for ls in (Optim.BackTracking,Optim.HagerZhang,Optim.StrongWolfe,Optim.MoreThuente)
            for method in (Optim.GradientDescent, Optim.ConjugateGradient, Optim.LBFGS, Optim.BFGS,
                           Optim.NGMRES, Optim.OACCEL)
                debug_printing && printstyled("Solver: $(summary(method())), linesearch: $(summary(ls()))\n", color=:green)
                res = Optim.optimize(fmanif, gmanif!, x0, method(manifold=manif,linesearch=ls()), Optim.Options(allow_f_increases=true,g_tol=1e-6))
                debug_printing && printstyled("Iter\tf-calls\tg-calls\n", color=:green)
                debug_printing && printstyled("$(Optim.iterations(res))\t$(Optim.f_calls(res))\t$(Optim.g_calls(res))\n", color=:red)
                @test Optim.converged(res)
            end
        end
        res = Optim.optimize(fmanif, gmanif!, x0, Optim.MomentumGradientDescent(mu=0.0, manifold=manif))
        @test Optim.converged(res)
    end

    @testset "Power and Product" begin
        # Power
        @views fprod(x) = fmanif(x[:,1]) + fmanif(x[:,2])
        @views gprod!(stor,x) = (gmanif!(stor[:, 1],x[:, 1]);gmanif!(stor[:, 2],x[:, 2]);stor)
        m1 = Optim.PowerManifold(Optim.Sphere(), (n,), (2,))
        rng = StableRNG(0)
        x0 = randn(rng, n, 2) + im*randn(rng, n, 2)
        res = Optim.optimize(fprod, gprod!, x0, Optim.ConjugateGradient(manifold=m1))
        @test Optim.converged(res)
        minpow = Optim.minimizer(res)

        # Product
        @views fprod(x) = fmanif(x[1:n]) + fmanif(x[n+1:2n])
        @views gprod!(stor,x) = (gmanif!(stor[1:n],x[1:n]);gmanif!(stor[n+1:2n],x[n+1:2n]);stor)
        m2 = Optim.ProductManifold(Optim.Sphere(), Optim.Sphere(), (n,), (n,))
        rng = StableRNG(0)
        x0 = randn(rng,2n) + im*randn(rng,2n)
        res = Optim.optimize(fprod, gprod!, x0, Optim.ConjugateGradient(manifold=m2))
        @test Optim.converged(res)
        minprod = Optim.minimizer(res)

        # results should be exactly equal: same initial guess, same sequence of operations
        @test minpow[:,1] == minprod[1:n]
        @test minpow[:,2] == minprod[n+1:2n]
    end
end


text:
module FastArrayOperations

function map(f::Function, array::Array{T} where T <:Number)
	# Implementation will be added shortly
	# For now it is just a template
	return array
end

end

max_stars_repo_name: PCwarlord/Julia-GPU-CPU-balancer

max_issues_repo_name: PCwarlord/Julia-GPU-CPU-balancer

max_forks_repo_name: PCwarlord/Julia-GPU-CPU-balancer


text:
"""
timedev(data, rate; [frequency=false], [overlapping=true], [taus=Octave])
Calculates the time deviation

#parameters:
* `<data>`:			The data array to calculate the deviation from either as as phases or frequencies.
* `<rate>`:			The rate of the data given.
* `[frequency]`:		True if `data` contains frequency data otherwise (default) phase data is assumed.
* `[overlapping]`:	True (default) to calculate overlapping deviation, false otherwise.
* `[taus]`:			Taus to calculate the deviation at. This can either be an AllanTauDescriptor type `AllTaus, Decadade, HalfDecade, Octave, HalfOctave, QuarterOctave`, an array of taus to calculate at, a float number to build a custom log-scale on or an integer to build a specific number of log spaced points.

#returns: named tupple (tau, deviation, error, count)
* `tau`:		Taus which where used.
* `deviation`:	Deviations calculated.
* `error`:		Respective errors.
* `count`:		Number of contributing terms for each deviation.
"""
function timedev(
		data::Array{T, 1},
		rate::AbstractFloat;
		frequency::Bool = false,
		overlapping::Bool = true,
		taus::Union{Type{U}, Integer, AbstractFloat, Array{Float64}} = 192) where {T, U <: AllanTauDescriptor}

	n = length(data)
	if n < 4
		error("Length for `data` in timedev must be at least 4 or greater")
		#we check this here, so that we can output the right function name in case of the error
	end

	(mdtaus, mddeviation, mderror, mdcount) = mallandev(data, rate, frequency = frequency, overlapping = overlapping, taus = taus)
	mdm = mdtaus ./ sqrt(3)

	(tau = mdtaus, deviation = mdm .* mddeviation, error = mdm .* mderror, count = mdcount)
end


text:
using Documenter

makedocs(;
    authors = "Julia Astro",
    repo = "https://github.com/juliaastro/juliaastro.github.io/blob/{commit}{path}#L{line}",
    sitename = "JuliaAstro",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://juliaastro.github.io/",
        assets = String[],
    ),
    pages = [
        "JuliaAstro" => "index.md",
    ],
    strict = true
)

deploydocs(;
    repo = "github.com/JuliaAstro/JuliaAstro.github.io",
    push_preview = true,
    branch = "master",
    devbranch = "source"
)