text:
export polyhedral, PolyhedralTracker, PolyhedralStartSolutionsIterator, mixed_volume

using MixedSubdivisions: MixedCell, mixed_volume

"""
    PolyhedralStartSolutionsIterator

An iterator providing start solutions for the polyhedral homotopy.
"""
struct PolyhedralStartSolutionsIterator
    support::Vector{Matrix{Int32}}
    start_coefficients::Vector{Vector{ComplexF64}}
    lifting::Vector{Vector{Int32}}
    mixed_cells::Vector{MixedCell}
    BSS::BinomialSystemSolver
end

function PolyhedralStartSolutionsIterator(
    support::AbstractVector{<:AbstractMatrix{<:Integer}},
    coeffs::AbstractVector{<:AbstractVector{<:Number}},
    lifting = map(c -> zeros(Int32, length(c)), coeffs),
    mixed_cells = MixedCell[],
)
    BSS = BinomialSystemSolver(length(support))
    PolyhedralStartSolutionsIterator(
        convert(Vector{Matrix{Int32}}, support),
        coeffs,
        lifting,
        mixed_cells,
        BSS,
    )
end

Base.show(io::IO, C::PolyhedralStartSolutionsIterator) =
    print(io, "PolyhedralStartSolutionsIterator()")
Base.show(
    io::IO,
    ::MIME"application/prs.juno.inline",
    x::PolyhedralStartSolutionsIterator,
) = x
Base.IteratorSize(::Type{<:PolyhedralStartSolutionsIterator}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{<:PolyhedralStartSolutionsIterator}) = Base.HasEltype()
Base.eltype(iter::PolyhedralStartSolutionsIterator) = Tuple{MixedCell,Vector{ComplexF64}}

function compute_mixed_cells!(iter::PolyhedralStartSolutionsIterator)
    if isempty(iter.mixed_cells)
        res = MixedSubdivisions.fine_mixed_cells(iter.support)
        if isnothing(res) || isempty(res[1])
            throw(OverflowError("Cannot compute a start system."))
        end
        mixed_cells, lifting = res
        empty!(iter.mixed_cells)
        append!(iter.mixed_cells, mixed_cells)

        for (i, w) in enumerate(lifting)
            empty!(iter.lifting[i])
            append!(iter.lifting[i], w)
        end
    end
    iter
end

function Base.iterate(iter::PolyhedralStartSolutionsIterator)
    compute_mixed_cells!(iter)

    isempty(iter.mixed_cells) && return nothing

    cell = iter.mixed_cells[1]
    solve(iter.BSS, iter.support, iter.start_coefficients, cell)

    x = [iter.BSS.X[i, 1] for i = 1:length(iter.support)]
    el = (cell, x)
    next_state = size(iter.BSS.X, 2) == 1 ? (2, 1) : (1, 2)
    el, next_state
end

function Base.iterate(iter::PolyhedralStartSolutionsIterator, (i, j)::Tuple{Int,Int})
    i > length(iter.mixed_cells) && return nothing

    cell = iter.mixed_cells[i]
    j == 1 && solve(iter.BSS, iter.support, iter.start_coefficients, cell)
    x = [iter.BSS.X[i, j] for i = 1:length(iter.support)]
    el = (cell, x)
    next_state = size(iter.BSS.X, 2) == j ? (i + 1, 1) : (i, j + 1)
    return el, next_state
end

# Tracker
function polyhedral_system(support)
    m = 0
    p = Variable[]
    coeffs = map(support) do A
        c = variables(:c, m+1:m+size(A, 2))
        m += size(A, 2)
        append!(p, c)
        c
    end

    System(support, coeffs; variables = variables(:x, 1:length(support)), parameters = p)
end

"""
    PolyhedralTracker <: AbstractPathTracker

This tracker realises the two step approach of the polyhedral homotopy.
See also [`polyhedral`].
"""
struct PolyhedralTracker{H1<:ToricHomotopy,H2<:AbstractHomotopy,M} <: AbstractPathTracker
    toric_tracker::Tracker{H1,M}
    generic_tracker::EndgameTracker{H2,M}
    support::Vector{Matrix{Int32}}
    lifting::Vector{Vector{Int32}}
end

"""
    polyhedral(F::Union{System, AbstractSystem};
        only_non_zero = false,
        endgame_options = EndgameOptions(),
        tracker_options = TrackerOptions())

Solve the system `F` in two steps: first solve a generic system derived from the support
of `F` using a polyhedral homotopy as proposed in [^HS95], then perform a
coefficient-parameter homotopy towards `F`.
This returns a path tracker ([`PolyhedralTracker`](@ref) or [`OverdeterminedTracker`](@ref)) and an iterator to compute the start solutions.

If `only_non_zero` is `true`, then only the solutions with non-zero coordinates are computed.
In this case the number of paths to track is equal to the
mixed volume of the Newton polytopes of `F`.

If `only_non_zero` is `false`, then all isolated solutions of `F` are computed.
In this case the number of paths to track is equal to the
mixed volume of the convex hulls of ``supp(F_i) ∪ \\{0\\}`` where ``supp(F_i)`` is the support
of ``F_i``. See also [^LW96].

    function polyhedral(
        support::AbstractVector{<:AbstractMatrix},
        coefficientss::AbstractVector{<:AbstractVector{<:Number}};
        kwargs...,
    )

It is also possible to provide directly the support and coefficients of the system `F` to be solved.

[^HS95]: Birkett Huber and Bernd Sturmfels. “A Polyhedral Method for Solving Sparse Polynomial Systems.” Mathematics of Computation, vol. 64, no. 212, 1995, pp. 1541–1555
[^LW96]: T.Y. Li and Xiaoshen Wang. "The BKK root count in C^n". Math. Comput. 65, 216 (October 1996), 1477–1484.

### Example

We consider a system `f` which has in total 6 isolated solutions,
but only 3 where all coordinates are non-zero.
```julia
@var x y
f = System([2y + 3 * y^2 - x * y^3, x + 4 * x^2 - 2 * x^3 * y])
tracker, starts = polyhedral(f; only_non_zero = false)
# length(starts) == 8
count(is_success, track.(tracker, starts)) # 6

tracker, starts = polyhedral(f; only_non_zero = true)
# length(starts) == 3
count(is_success, track.(tracker, starts)) # 3
```
"""
function polyhedral(F::AbstractSystem; kwargs...)
    _, n = size(F)
    @var x[1:n]
    polyhedral(System(F(x), x); kwargs...)
end
function polyhedral(
    f::System;
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
    target_parameters = nothing,
    kwargs...,
)
    if target_parameters !== nothing
        return polyhedral(
            FixedParameterSystem(fixed(f; compile = compile), target_parameters);
            kwargs...,
        )
    end
    homogeneous = is_homogeneous(f)
    if homogeneous
        F = on_affine_chart(f; compile = compile)
        m, n = size(F)
        m ≥ n || throw(FiniteException(n - m))
        if m > n
            F = square_up(F; compile = compile)
        end
        @var x[1:n]
        support, target_coeffs = support_coefficients(System(F(x), x))
    else
        m, n = size(f)
        m ≥ n || throw(FiniteException(n - m))
        if m > n
            F = square_up(f; compile = compile)
            @var x[1:n]
            support, target_coeffs = support_coefficients(System(F(x), x))
        else
            support, target_coeffs = support_coefficients(f)
        end
    end
    tracker, starts = polyhedral(support, target_coeffs; compile = compile, kwargs...)
    if m > n
        tracker = OverdeterminedTracker(tracker, F)
    end
    tracker, starts
end

function has_zero_col(A)
    for j = 1:size(A, 2)
        has_non_zero = false
        for i = 1:size(A, 1)
            if !iszero(A[i, j])
                has_non_zero = true
                break
            end
        end
        has_non_zero || return true
    end
    false
end

paths_to_track(f::AbstractSystem, val::Val{:polyhedral}) = paths_to_track(System(f), val)
function paths_to_track(
    f::System,
    ::Val{:polyhedral};
    only_torus::Bool = false,
    only_non_zero::Bool = only_torus,
)
    supp, coeffs = begin
        if is_homogeneous(f)
            F = on_affine_chart(f)
            m, n = size(F)
            m ≥ n || throw(FiniteException(n - m))
            if m > n
                F = square_up(F)
            end
            @var x[1:n]
            support_coefficients(System(F(x), x))
        else
            m, n = size(f)
            m ≥ n || throw(FiniteException(n - m))
            if m > n
                F = square_up(f)
                @var x[1:n]
                support_coefficients(System(F(x), x))
            else
                support_coefficients(f)
            end
        end
    end

    if only_non_zero
        MixedSubdivisions.mixed_volume(supp)
    else
        supp′ = map(A -> has_zero_col(A) ? A : [A zeros(Int32, size(A, 1), 1)], supp)
        MixedSubdivisions.mixed_volume(supp′)
    end
end
MixedSubdivisions.mixed_volume(f::Union{System,AbstractSystem}) =
    paths_to_track(f; start_system = :polyhedral, only_torus = true)

function polyhedral(
    support::AbstractVector{<:AbstractMatrix},
    target_coeffs::AbstractVector;
    kwargs...,
)
    start_coeffs = map(c -> rand_approx_unit(length(c)) .* LA.norm(c, Inf), target_coeffs)
    polyhedral(support, start_coeffs, target_coeffs; kwargs...)
end

"""
    rand_approx_unit(n::Integer)::Vector{ComplexF64}

This samples uniformly from the rectangle ``[0.9,1.1] × [0,2π]`` and transforms the sampled
values with the map ``(r, φ) ↦ r * e^{i φ}``.
"""
rand_approx_unit(n) = map(_ -> (0.9 + 0.2 * rand()) * cis2pi(rand()), 1:n)

cis2pi(x) = complex(cospi(2x), sinpi(2x))

function polyhedral(
    support::AbstractVector{<:AbstractMatrix},
    start_coeffs::AbstractVector,
    target_coeffs::AbstractVector;
    endgame_options = EndgameOptions(),
    tracker_options = TrackerOptions(),
    only_torus::Bool = false,
    only_non_zero::Bool = only_torus,
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
    kwargs...,
)
    unsupported_kwargs(kwargs)
    size.(support, 2) == length.(start_coeffs) == length.(target_coeffs) ||
        throw(ArgumentError("Number of terms do not coincide."))
    if only_non_zero
        min_vecs = minimum.(support, dims = 2)
        support = map((A, v) -> A .- v, support, min_vecs)
    elseif !all(has_zero_col, support)
        # Add 0 to each support
        starts = Vector{ComplexF64}[]
        targets = Vector{ComplexF64}[]
        supp = eltype(support)[]
        for (i, A) in enumerate(support)
            if has_zero_col(A)
                push!(starts, start_coeffs[i])
                push!(targets, target_coeffs[i])
                push!(supp, A)
            else
                push!(starts, vcat(start_coeffs[i], randn(ComplexF64)))
                push!(targets, vcat(target_coeffs[i], 0.0))
                push!(supp, [A zeros(Int32, size(A, 1), 1)])
            end
        end
        support = supp
        start_coeffs = starts
        target_coeffs = targets
    end

    F = fixed(polyhedral_system(support); compile = compile)
    H₁ = ToricHomotopy(F, start_coeffs)
    toric_tracker = Tracker(H₁; options = tracker_options)

    H₂ = begin
        p = reduce(append!, start_coeffs; init = ComplexF64[])
        q = reduce(append!, target_coeffs; init = ComplexF64[])
        CoefficientHomotopy(F, p, q)
    end
    generic_tracker =
        EndgameTracker(Tracker(H₂; options = tracker_options), options = endgame_options)

    S = PolyhedralStartSolutionsIterator(support, start_coeffs)
    tracker = PolyhedralTracker(toric_tracker, generic_tracker, S.support, S.lifting)

    tracker, S
end

function track(
    PT::PolyhedralTracker,
    start_solution::Tuple{MixedCell,<:AbstractVector{ComplexF64}};
    path_number::Union{Nothing,Int} = nothing,
    debug::Bool = false,
)
    cell, x∞ = start_solution
    H = PT.toric_tracker.homotopy
    # The tracker works in two stages
    # 1) Revert the totric degeneration by tracking 0 to 1
    #    Here we use the following strategy
    #    a) Reparameterize the path such that lowest power of t occuring is 1
    #    b)
    #      i) If maximal power is less than 10 track from 0 to 1
    #      ii) If maximal power is larger than 10 we can get some problems for values close to 1
    #          (since t^k is << 1 for t < 1 and k large)
    #          Therefore split in two stages:
    #           1) track from 0 to t₀ = clamp(0.1^(10 / max_weight), 0.9, 0.985)
    #           2) Second, reperamerize path s.t. max power of t is 10. This shifts the
    #              problem towards 0 again (where we have more accuracy available.)
    #              And then track to 1
    # 2) Perform coefficient homotopy with possible endgame
    min_weight, max_weight =
        update_weights!(H, PT.support, PT.lifting, cell, min_weight = 1.0)

    if debug
        println("Min-Max weight: ", min_weight, ", ", max_weight)
    end
    if max_weight < 10
        retcode = track!(
            PT.toric_tracker,
            x∞,
            0.0,
            1.0;
            ω = 20.0,
            μ = 1e-12,
            max_initial_step_size = 0.2,
            debug = debug,
        )
        @unpack μ, ω = PT.toric_tracker.state
    else
        t₀ = clamp(0.1^(10 / max_weight), 0.9, 1 - 1e-6)
        retcode = track!(
            PT.toric_tracker,
            x∞,
            0.0,
            t₀;
            ω = 20.0,
            μ = 1e-12,
            max_initial_step_size = 0.2,
            debug = debug,
        )

        @unpack μ, ω = PT.toric_tracker.state
        if is_success(retcode)
            min_weight, max_weight =
                update_weights!(H, PT.support, PT.lifting, cell, max_weight = 10.0)
            t_restart = t₀^(1 / min_weight)
            # set min_step_size to 0.0 to not accidentally loose a solution
            min_step_size = PT.toric_tracker.options.min_step_size
            PT.toric_tracker.options.min_step_size = 0.0

            retcode = track!(
                PT.toric_tracker,
                PT.toric_tracker.state.x,
                t_restart,
                1.0;
                ω = ω,
                μ = μ,
                τ = 0.1 * t_restart,
                keep_steps = true,
                debug = debug,
            )
            @unpack μ, ω = PT.toric_tracker.state
            PT.toric_tracker.options.min_step_size = min_step_size
        end
    end
    if !is_success(retcode)
        state = PT.toric_tracker.state
        return PathResult(
            return_code = :polyhedral_failed,
            solution = copy(state.x),
            start_solution = start_solution,
            t = real(state.t),
            accuracy = state.accuracy,
            singular = false,
            condition_jacobian = NaN,
            residual = NaN,
            winding_number = nothing,
            last_path_point = (copy(state.x), state.t),
            valuation = nothing,
            ω = state.ω,
            μ = state.μ,
            path_number = path_number,
            extended_precision = state.extended_prec,
            accepted_steps = state.accepted_steps,
            rejected_steps = state.rejected_steps,
            extended_precision_used = state.used_extended_prec,
        )
    end

    r = track(
        PT.generic_tracker,
        PT.toric_tracker.state.x;
        # Don't provide ω since this can be misleading a lead to a too large initial step
        # ω = ω,
        μ = μ,
        path_number = path_number,
        debug = debug,
    )
    r.start_solution = start_solution
    # report accurate steps
    r.accepted_steps += PT.toric_tracker.state.accepted_steps
    r.rejected_steps += PT.toric_tracker.state.rejected_steps

    r
end


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

using FEMBase, AsterReader, Test, LinearAlgebra, SparseArrays

include(joinpath("..", "docs", "make.jl"))

@testset "AsterReader.jl" begin
    include("test_read_aster_mesh.jl")
    include("test_read_aster_results.jl")
end

include(joinpath("..", "docs", "deploy.jl"))

ext: jl


text:
# ### m8.1t

using Pkg, DrWatson

@quickactivate "StatisticalRethinkingTuring"
using Turing
using StatisticalRethinking
Turing.setprogress!(false)

# Read in the `rugged` data as a DataFrame

df = CSV.read(sr_datadir("rugged.csv"), DataFrame)

# Apply log() to each element in rgdppc_2000 column and add it as a new column

df = hcat(df, map(log, df[!, Symbol("rgdppc_2000")]));

# Rename our col x1 => log_gdp

DataFrames.rename!(df, :x1 => :log_gdp);

# Now we need to drop every row where rgdppc_2000 == missing

# When this (https://github.com/JuliaData/DataFrames.jl/pull/1546) hits DataFrame it'll be conceptually easier: i.e., completecases!(d, :rgdppc_2000)

notisnan(e) = !ismissing(e)
dd = df[map(notisnan, df[:, :rgdppc_2000]), :];

# Define the Turing model

@model ppl8_1(y, x₁, x₂) = begin
    σ ~ truncated(Cauchy(0, 2), 0, Inf)
    βR ~ Normal(0, 10)
    βA ~ Normal(0, 10)
    βAR ~ Normal(0, 10)
    α ~ Normal(0, 100)

    for i ∈ 1:length(y)
        y[i] ~ Normal(α + βR * x₁[i] + βA * x₂[i] + βAR * x₁[i] * x₂[i], σ)
    end

    # alternative formulation using broadcasting
    # y .~ Normal.(α .+ βR * x₁ .+ βA * x₂ .+ βAR * x₁ .* x₂, σ)
end;

# Use Turing mcmc

m8_1t = ppl8_1(dd.log_gdp, dd.rugged, dd.cont_africa)
nchains = 4; sampler = NUTS(0.65); nsamples=2000
chns8_1t = mapreduce(c -> sample(m8_1t, sampler, nsamples), chainscat, 1:nchains)

# Here's the ulam() output from rethinking

m8_1s_cmdstan = "
Iterations = 1:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 1000

Empirical Posterior Estimates:
          Mean         SD       Naive SE       MCSE      ESS
    a  9.22360053 0.139119116 0.0021996664 0.0034632816 1000
   bR -0.20196346 0.076106388 0.0012033477 0.0018370185 1000
   bA -1.94430980 0.227080488 0.0035904578 0.0057840746 1000
  bAR  0.39071684 0.131889143 0.0020853505 0.0032749642 1000
sigma  0.95036370 0.052161768 0.0008247500 0.0009204073 1000

Quantiles:
          2.5%       25.0%       50.0%      75.0%        97.5%
    a  8.95307475  9.12719750  9.2237750  9.31974000  9.490234250
   bR -0.35217930 -0.25334425 -0.2012855 -0.15124725 -0.054216855
   bA -2.39010825 -2.09894500 -1.9432550 -1.78643000 -1.513974250
  bAR  0.13496995  0.30095575  0.3916590  0.47887625  0.650244475
sigma  0.85376115  0.91363250  0.9484920  0.98405750  1.058573750
";

# End of `08/m8.1t.jl`


text:
### Data Parser ###

# Data Read: Reads the structured data from the given instance file
# Input:  - filename: A string containing the file name with its relative path.
# Output: - n: An Integer containing the number of tasks;
#         - tasks_times: An (n) Array of Integers with each task's time;
#         - prec_relations: A (n,n) 2-dimensional Binary Array representing
#                           each task's Precedence Relation.
function parser_data_read(filename::String)
    # Gets every line in the instance file for processing
    lines = []
    open(filename,"r") do f
        lines = readlines(filename)
    end

    # The number n of tasks is given in the first line
    n = parse(Int64, lines[1])

    # Gets each task's time in the next n lines
    task_times = Int64[]
    for i = 2:(n+1)
        append!(task_times, parse(Int64, lines[i]))
    end

    # Gets the i,j Direct Precedence Relations until the "-1,-1" end mark
    # and transforms it into a Binary Matrix of the Relations
    prec_relations = zeros(Int64, (n,n))
    curr_line = n+2
    # Julia has no "Do While" so a first verification is necessary to avoid a break
    rs = split(lines[curr_line], ",")
    r = [parse(Int64, rs[1]), parse(Int64, rs[2])]
    while (r[1] > 0)
        prec_relations[r[1], r[2]] = 1

        curr_line += 1
        rs = split(lines[curr_line], ",")
        r = [parse(Int64, rs[1]), parse(Int64, rs[2])]
    end

    println("Data successfully read from file \"", filename, "\" containing ", n, " tasks.")
    println("")

    # Returns the structured data
    return n, task_times, prec_relations
end

# Parser Write: Parses the Data from the given Instance to fit GLPK pattern.
# Input:  - filename: A string with the name for the parser file that will be written;
#         - n, s, task_times, times_sum, prec_relations: Needed to write the Data.
# Output: None.
function parser_write(filename::String, n, s, task_times, prec_relations)
    open(filename, "w") do f
        print(f, "data;\r\n")
        print(f, "\r\n")
        print(f, "param n := ", n, ";\r\n")
        print(f, "param m := ", s, ";\r\n")
        print(f, "param times := ")
        for i = 1:n
            print(f, i, " ", task_times[i], "\r\n")
        end
        print(f, ";\r\n")
        print(f, "precedences : ")
        for i = 1:n
            print(f, i, " ")
        end
        print(f, ":=\r\n")
        for i = 1:n
            prec_line = string(prec_relations[i, :])
            prec_line = replace(prec_line, ",", "")
            prec_line = replace(prec_line, "[", "")
            prec_line = replace(prec_line, "]", "")
            print(f, i, " ", prec_line, "\r\n")
        end
        print(f, ";\r\n")
    end

    println("Solution saved in file: ", filename)
    println("")
end

### Main Function ###
function main()
    s = 3

    filename = ["instances/HAHN.IN2", "instances/LUTZ3.IN2", "instances/WEE-MAG.IN2"]
    out_file = ["instances/HAHN.DAT", "instances/LUTZ3.DAT", "instances/WEE-MAG.DAT"]

    for i = 1:(length(filename))
        n, task_times, prec_relations = parser_data_read(filename[i])
        parser_write(out_file[i], n, s, task_times, prec_relations)
    end

    println("Parsing finished.")
end

main()

max_stars_repo_name: AlexandreSSJr/SALBP-GA

max_issues_repo_name: AlexandreSSJr/SALBP-GA

max_forks_repo_name: AlexandreSSJr/SALBP-GA


text:
using CellularAutomata

#Rule 90 1-dim Cellular Automata run for 15 generations
function rule90example()
    #init=int(zeros(1000))
    #init[int(1000/2)]=1
    tic()
    ca = CellularAutomaton(90, 15)
    toc()

    ca
end

#Turbine pattern from Conway's game of life
function turbine_example()
    init = int(zeros(20,20))
    init[4, 4:9] = 1
    init[5, 4:9] = 1

    init[4:9, 11] = 1
    init[4:9, 12] = 1

    init[11, 7:12] = 1
    init[12, 7:12] = 1

    init[7:12, 4] = 1
    init[7:12, 5] = 1

    #game of life
    CA2d([3], [2,3], init, 5)
end


text:
function gencodeall(npm)
    # Code generation for all models

    gencodeelemface("Flux", npm, 1);
    gencodeelemface("Source", npm, 1);
    gencodeelemface("Tdfunc", npm, 1);

    gencodeelemface("Avfield", npm, 2);
    gencodeelemface("Output", npm, 2);

    gencodeelemface("Fbou", npm, 3);
    gencodeelemface("Ubou", npm, 3);

    gencodeelemface("Fhat", npm, 4);
    gencodeelemface("Uhat", npm, 4);
    gencodeelemface("Stab", npm, 4);

    gencodeelemface("Initu", npm, 5);
    gencodeelemface("Initq", npm, 5);
    gencodeelemface("Initudg", npm, 5);
    gencodeelemface("Initwdg", npm, 5);
    gencodeelemface("Initodg", npm, 5);

end

max_stars_repo_name: rloekvh/Exasim

max_issues_repo_name: rloekvh/Exasim

max_forks_repo_name: rloekvh/Exasim


text:
# Julia wrapper for header: types.h
# Automatically generated using Clang.jl

# Julia wrapper for header: cutensor.h
# Automatically generated using Clang.jl

function cutensorCreate(handle)
    @check @runtime_ccall((:cutensorCreate, libcutensor), cutensorStatus_t,
                 (Ptr{cutensorHandle_t},),
                 handle)
end

function cutensorDestroy(handle)
    @check @runtime_ccall((:cutensorDestroy, libcutensor), cutensorStatus_t,
                 (cutensorHandle_t,),
                 handle)
end

function cutensorCreateTensorDescriptor(desc, numModes, extent, stride, dataType, unaryOp,
                                        vectorWidth, vectorModeIndex)
    @check @runtime_ccall((:cutensorCreateTensorDescriptor, libcutensor), cutensorStatus_t,
                 (Ptr{cutensorTensorDescriptor_t}, UInt32, Ptr{Int64}, Ptr{Int64},
                  cudaDataType_t, cutensorOperator_t, UInt32, UInt32),
                 desc, numModes, extent, stride, dataType, unaryOp, vectorWidth,
                 vectorModeIndex)
end

function cutensorDestroyTensorDescriptor(desc)
    @check @runtime_ccall((:cutensorDestroyTensorDescriptor, libcutensor), cutensorStatus_t,
                 (cutensorTensorDescriptor_t,),
                 desc)
end

function cutensorElementwiseTrinary(handle, alpha, A, descA, modeA, beta, B, descB, modeB,
                                    gamma, C, descC, modeC, D, descD, modeD, opAB, opABC,
                                    typeCompute, stream)
    @check @runtime_ccall((:cutensorElementwiseTrinary, libcutensor), cutensorStatus_t,
                 (cutensorHandle_t, Ptr{Cvoid}, CuPtr{Cvoid}, cutensorTensorDescriptor_t,
                  Ptr{Int32}, Ptr{Cvoid}, CuPtr{Cvoid}, cutensorTensorDescriptor_t,
                  Ptr{Int32}, Ptr{Cvoid}, CuPtr{Cvoid}, cutensorTensorDescriptor_t,
                  Ptr{Int32}, CuPtr{Cvoid}, cutensorTensorDescriptor_t, Ptr{Int32},
                  cutensorOperator_t, cutensorOperator_t, cudaDataType_t, CUstream),
                 handle, alpha, A, descA, modeA, beta, B, descB, modeB, gamma, C, descC,
                 modeC, D, descD, modeD, opAB, opABC, typeCompute, stream)
end

function cutensorElementwiseBinary(handle, alpha, A, descA, modeA, gamma, C, descC, modeC,
                                   D, descD, modeD, opAC, typeCompute, stream)
    @check @runtime_ccall((:cutensorElementwiseBinary, libcutensor), cutensorStatus_t,
                 (cutensorHandle_t, Ptr{Cvoid}, CuPtr{Cvoid}, cutensorTensorDescriptor_t,
                  Ptr{Int32}, Ptr{Cvoid}, CuPtr{Cvoid}, cutensorTensorDescriptor_t,
                  Ptr{Int32}, CuPtr{Cvoid}, cutensorTensorDescriptor_t, Ptr{Int32},
                  cutensorOperator_t, cudaDataType_t, CUstream),
                 handle, alpha, A, descA, modeA, gamma, C, descC, modeC, D, descD, modeD,
                 opAC, typeCompute, stream)
end

function cutensorPermutation(handle, alpha, A, descA, modeA, B, descB, modeB, typeCompute,
                             stream)
    @check @runtime_ccall((:cutensorPermutation, libcutensor), cutensorStatus_t,
                 (cutensorHandle_t, Ptr{Cvoid}, CuPtr{Cvoid}, cutensorTensorDescriptor_t,
                  Ptr{Int32}, CuPtr{Cvoid}, cutensorTensorDescriptor_t, Ptr{Int32},
                  cudaDataType_t, CUstream),
                 handle, alpha, A, descA, modeA, B, descB, modeB, typeCompute, stream)
end

function cutensorContraction(handle, alpha, A, descA, modeA, B, descB, modeB, beta, C,
                             descC, modeC, D, descD, modeD, opOut, typeCompute, algo,
                             workspace, workspaceSize, stream)
    @check @runtime_ccall((:cutensorContraction, libcutensor), cutensorStatus_t,
                 (cutensorHandle_t, Ptr{Cvoid}, CuPtr{Cvoid}, cutensorTensorDescriptor_t,
                  Ptr{Int32}, CuPtr{Cvoid}, cutensorTensorDescriptor_t, Ptr{Int32},
                  Ptr{Cvoid}, CuPtr{Cvoid}, cutensorTensorDescriptor_t, Ptr{Int32},
                  CuPtr{Cvoid}, cutensorTensorDescriptor_t, Ptr{Int32}, cutensorOperator_t,
                  cudaDataType_t, cutensorAlgo_t, CuPtr{Cvoid}, UInt64, CUstream),
                 handle, alpha, A, descA, modeA, B, descB, modeB, beta, C, descC, modeC, D,
                 descD, modeD, opOut, typeCompute, algo, workspace, workspaceSize, stream)
end

function cutensorContractionGetWorkspace(handle, A, descA, modeA, B, descB, modeB, C,
                                         descC, modeC, D, descD, modeD, opOut, typeCompute,
                                         algo, pref, workspaceSize)
    @check @runtime_ccall((:cutensorContractionGetWorkspace, libcutensor), cutensorStatus_t,
                 (cutensorHandle_t, CuPtr{Cvoid}, cutensorTensorDescriptor_t, Ptr{Int32},
                  CuPtr{Cvoid}, cutensorTensorDescriptor_t, Ptr{Int32}, CuPtr{Cvoid},
                  cutensorTensorDescriptor_t, Ptr{Int32}, CuPtr{Cvoid},
                  cutensorTensorDescriptor_t, Ptr{Int32}, cutensorOperator_t,
                  cudaDataType_t, cutensorAlgo_t, cutensorWorksizePreference_t, Ptr{UInt64}),
                 handle, A, descA, modeA, B, descB, modeB, C, descC, modeC, D, descD,
                 modeD, opOut, typeCompute, algo, pref, workspaceSize)
end

function cutensorContractionMaxAlgos(maxNumAlgos)
    @check @runtime_ccall((:cutensorContractionMaxAlgos, libcutensor), cutensorStatus_t,
                 (Ptr{Int32},),
                 maxNumAlgos)
end

function cutensorReduction(handle, alpha, A, descA, modeA, beta, C, descC, modeC, D, descD,
                           modeD, opReduce, typeCompute, workspace, workspaceSize, stream)
    @check @runtime_ccall((:cutensorReduction, libcutensor), cutensorStatus_t,
                 (cutensorHandle_t, Ptr{Cvoid}, CuPtr{Cvoid}, cutensorTensorDescriptor_t,
                  Ptr{Int32}, Ptr{Cvoid}, CuPtr{Cvoid}, cutensorTensorDescriptor_t,
                  Ptr{Int32}, CuPtr{Cvoid}, cutensorTensorDescriptor_t, Ptr{Int32},
                  cutensorOperator_t, cudaDataType_t, CuPtr{Cvoid}, UInt64, CUstream),
                 handle, alpha, A, descA, modeA, beta, C, descC, modeC, D, descD, modeD,
                 opReduce, typeCompute, workspace, workspaceSize, stream)
end

function cutensorReductionGetWorkspace(handle, A, descA_, modeA, C, descC_, modeC, D,
                                       descD_, modeD, opReduce, typeCompute, workspaceSize)
    @check @runtime_ccall((:cutensorReductionGetWorkspace, libcutensor), cutensorStatus_t,
                 (cutensorHandle_t, CuPtr{Cvoid}, cutensorTensorDescriptor_t, Ptr{Int32},
                  CuPtr{Cvoid}, cutensorTensorDescriptor_t, Ptr{Int32}, CuPtr{Cvoid},
                  cutensorTensorDescriptor_t, Ptr{Int32}, cutensorOperator_t,
                  cudaDataType_t, Ptr{UInt64}),
                 handle, A, descA_, modeA, C, descC_, modeC, D, descD_, modeD, opReduce,
                 typeCompute, workspaceSize)
end

function cutensorGetErrorString(error)
    @runtime_ccall((:cutensorGetErrorString, libcutensor), Cstring,
          (cutensorStatus_t,),
          error)
end

function cutensorGetVersion()
    @runtime_ccall((:cutensorGetVersion, libcutensor), Csize_t, ())
end

function cutensorGetCudartVersion()
    @runtime_ccall((:cutensorGetCudartVersion, libcutensor), Csize_t, ())
end


text:
# ================================================
# Concrete overapproximation with a hyperrectangle
# ================================================

"""
    box_approximation(S::LazySet{N}) where {N}

Overapproximate a set by a tight hyperrectangle.

### Input

- `S` -- set

### Output

A tight hyperrectangle.

### Algorithm

The center of the hyperrectangle is obtained by averaging the support function
of the given set in the canonical directions, and the lengths of the sides can
be recovered from the distance among support functions in the same directions.
"""
function box_approximation(S::LazySet{N}) where {N}
    c, r = box_approximation_helper(S)
    if r[1] < 0
        return EmptySet{N}(dim(S))
    end
    return Hyperrectangle(c, r)
end

"""
    interval_hull

Alias for `box_approximation`.
"""
interval_hull = box_approximation

"""
    box_approximation_helper(S::LazySet{N}) where {N}

Common code of `box_approximation` and `box_approximation_symmetric`.

### Input

- `S` -- set

### Output

A tuple containing the data that is needed to construct a tightly
overapproximating hyperrectangle.

- `c` -- center
- `r` -- radius

### Algorithm

The center of the hyperrectangle is obtained by averaging the support function
of the given set in the canonical directions.
The lengths of the sides can be recovered from the distance among support
functions in the same directions.
"""
@inline function box_approximation_helper(S::LazySet{N}) where {N}
    n = dim(S)
    c = Vector{N}(undef, n)
    r = Vector{N}(undef, n)
    @inbounds for i in 1:n
        htop = ρ(SingleEntryVector(i, n, one(N)), S)
        hbottom = -ρ(SingleEntryVector(i, n, -one(N)), S)
        c[i] = (htop + hbottom) / 2
        r[i] = (htop - hbottom) / 2
        if r[i] < 0
            # contradicting bounds => set is empty
            # terminate with first radius entry being negative
            r[1] = r[i]
            break
        end
    end
    return c, r
end

# ===============
# Specializations
# ===============

# hyperrectangle specialization
function box_approximation(S::AbstractHyperrectangle)
    return Hyperrectangle(center(S), radius_hyperrectangle(S))
end

# empty set specialization
box_approximation(∅::EmptySet) = ∅

"""
    box_approximation(S::CartesianProductArray{N, <:AbstractHyperrectangle}) where {N}

Return a tight overapproximation of the Cartesian product array of a finite
number of convex sets with and hyperrectangle.

### Input

- `S`              -- Cartesian product array of a finite number of convex set
- `Hyperrectangle` -- type for dispatch

### Output

A hyperrectangle.

### Algorithm

This method falls back to the corresponding `convert` method. Since the sets wrapped
by the Cartesian product array are hyperrectangles, it can be done efficiently
without overapproximation.
"""
function box_approximation(S::CartesianProductArray{N, <:AbstractHyperrectangle}) where {N}
    return convert(Hyperrectangle, S)
end

"""
    box_approximation(S::CartesianProduct{N, <:AbstractHyperrectangle, <:AbstractHyperrectangle}) where {N}

Return a tight overapproximation of the Cartesian product of two
hyperrectangles by a new hyperrectangle.

### Input

- `S`              -- Cartesian product of two hyperrectangular sets
- `Hyperrectangle` -- type for dispatch

### Output

A hyperrectangle.

### Algorithm

This method falls back to the corresponding `convert` method. Since the sets wrapped
by the Cartesian product are hyperrectangles, it can be done efficiently
without overapproximation.
"""
function box_approximation(S::CartesianProduct{N, <:AbstractHyperrectangle, <:AbstractHyperrectangle}) where {N}
    return convert(Hyperrectangle, S)
end

"""
    box_approximation(lm::LinearMap{N, <:AbstractHyperrectangle}) where {N}

Return a tight overapproximation of the linear map of a hyperrectangular set
using a hyperrectangle.

### Input

- `S`              -- linear map of a hyperrectangular set
- `Hyperrectangle` -- type for dispatch

### Output

A hyperrectangle.

### Algorithm

If `c` and `r` denote the center and vector radius of a hyperrectangle `H`,
a tight hyperrectangular overapproximation of `M * H` is obtained by transforming
`c ↦ M*c` and `r ↦ abs.(M) * r`, where `abs.(⋅)` denotes the element-wise absolute
value operator.
"""
function box_approximation(lm::LinearMap{N, <:AbstractHyperrectangle}) where {N}
    M, X = lm.M, lm.X
    center_MX = M * center(X)
    radius_MX = abs.(M) * radius_hyperrectangle(X)
    return Hyperrectangle(center_MX, radius_MX)
end

"""
    box_approximation(r::Rectification{N}) where {N}

Overapproximate the rectification of a convex set by a tight hyperrectangle.

### Input

- `r`              -- rectification of a convex set
- `Hyperrectangle` -- type for dispatch

### Output

A hyperrectangle.

### Algorithm

Box approximation and rectification distribute.
Hence we first check whether the wrapped set is empty.
If so, we return the empty set.
Otherwise, we compute the box approximation of the wrapped set, rectify the
resulting box (which is simple), and finally convert the resulting set to a box.
"""
function box_approximation(r::Rectification{N}) where {N}
    if isempty(r.X)
        return EmptySet{N}(dim(r))
    end
    return convert(Hyperrectangle, Rectification(box_approximation(r.X)))
end

"""
    box_approximation(Z::AbstractZonotope)

Return a tight overapproximation of a zonotope with an axis-aligned box.

### Input

- `Z`              -- zonotope
- `Hyperrectangle` -- type for dispatch

### Output

A hyperrectangle.

### Algorithm

This function implements the method in [Section 5.1.2, 1]. A zonotope
``Z = ⟨c, G⟩`` can be overapproximated tightly by an axis-aligned box
(i.e. a `Hyperrectangle`) such that its center is ``c`` and the radius along
dimension ``i`` is the column-sum of the absolute values of the ``i``-th row
of ``G`` for ``i = 1,…, p``, where ``p`` is the number of generators of ``Z``.

[1] *Althoff, M., Stursberg, O., & Buss, M. (2010). Computing reachable sets of
hybrid systems using a combination of zonotopes and polytopes. Nonlinear analysis:
hybrid systems, 4(2), 233-249.*
"""
function box_approximation(Z::AbstractZonotope)
    r = sum(abs.(genmat(Z)), dims=2)[:]
    return Hyperrectangle(center(Z), r)
end

"""
    box_approximation(am::AbstractAffineMap{N, <:AbstractHyperrectangle}) where {N}

Overapproximate the affine map of a hyperrectangular set
using a hyperrectangle.

### Input

- `am`             -- affine map of a hyperrectangular set
- `Hyperrectangle` -- type for dispatch

### Output

A hyperrectangle.

### Algorithm

If `c` and `r` denote the center and vector radius of a hyperrectangle `H`
and `v` the translation vector, a tight hyperrectangular overapproximation of
`M * H + v` is obtained by transforming `c ↦ M*c+v` and `r ↦ abs.(M) * r`, where
`abs.(⋅)` denotes the element-wise absolute value operator.
"""
function box_approximation(am::AbstractAffineMap{N, <:AbstractHyperrectangle}) where {N}
    M, X, v = matrix(am), set(am), vector(am)
    center_MXv = M * center(X) + v
    radius_MX = abs.(M) * radius_hyperrectangle(X)
    return Hyperrectangle(center_MXv, radius_MX)
end

function box_approximation(P::Union{VPolytope, VPolygon})
    n = dim(P)
    vlist = vertices_list(P)
    @assert !isempty(vlist) "cannot overapproximate an empty polytope"

    @inbounds v1 = vlist[1]
    center = similar(v1)
    radius = similar(v1)

    @inbounds for i in 1:n
        low_i = v1[i]
        high_i = v1[i]
        for v in vlist
            if v[i] > high_i
                high_i = v[i]
            elseif v[i] < low_i
                low_i = v[i]
            end
        end
        radius_i = (high_i - low_i) / 2
        center[i] = low_i + radius_i
        radius[i] = radius_i
    end
    return Hyperrectangle(center, radius)
end


text:
## spectral.jl : Gaussian random field generator based on the spectral (eigenvalue) decomposition of the covariance matrix

"""
    Spectral()

Returns a [`GaussianRandomFieldGenerator`](@ref) based on a spectral (eigenvalue) decomposition of the covariance matrix.

# Optional Arguments for `GaussianRandomField`
- `n::Integer`: the number of eigenvalues to compute. By default, we compute all eigenvalues.
- `eigensolver::EigenSolver`: which method to use for the eigenvalue decomposition (see [`AbstractEigenSolver`](@ref)). The default is `EigenSolver()`.

# Examples
```jldoctest
julia> cov = CovarianceFunction(2, Matern(.3, 1))
2d Matérn covariance function (λ=0.3, ν=1.0, σ=1.0, p=2.0)

julia> pts = range(0, stop=1, length=51)
0.0:0.02:1.0

julia> grf = GaussianRandomField(cov, Spectral(), pts, pts)
Gaussian random field with 2d Matérn covariance function (λ=0.3, ν=1.0, σ=1.0, p=2.0) on a 51×51 structured grid, using a spectral decomposition

julia> heatmap(grf)
[...]

```
!!! tip

    Try using the Karhunen-Loève expansion if evaluating the covariance matrix is too expensive.

This is also useful when computing Gaussian random fields on a Finite Element mesh using a truncated KL expansion. Here's an example that computes the first 10 eigenfunctions on an L-shaped domain.

```jldoctest
julia> cov = CovarianceFunction(2, Matern(.3, 1))
2d Matérn covariance function (λ=0.3, ν=1.0, σ=1.0, p=2.0)

julia> p, t = Lshape()
[...]

julia> grf = GaussianRandomField(cov, Spectral(), p, t, n=10)
Gaussian random field with 2d Matérn covariance function (λ=0.3, ν=1.0, σ=1.0, p=2.0) on a mesh with 998 points and 1861 elements, using a spectral decomposition

```
See also: [`Cholesky`](@ref), [`KarhunenLoeve`](@ref), [`CirculantEmbedding`](@ref)
"""
struct Spectral <: GaussianRandomFieldGenerator end

# container type for eigenvalues and eigenfunctions
struct SpectralData{X,Y}
    eigenval::X
    eigenfunc::Y
end

function _GaussianRandomField(mean, cov::CovarianceFunction, method::Spectral, pts...;
                              n::Integer=0, eigensolver=EigenSolver())
    C = apply(cov, pts...)

    # compute eigenvalue decomposition
    if n == 0
        F = eigen(C)
        idx = sortperm(F.values,rev=true)
        Λ = F.values[idx]
        U = F.vectors[:,idx]
    else # thanks to #4
        eigenval, eigenfunc = compute(C, n, eigensolver)
        idx = sortperm(eigenval,rev=true)[1:n]
        Λ = eigenval[idx]
        U = eigenfunc[:,idx]
    end

    # if negative eigenvalues detected, remove them
    m = findfirst(x -> x < 0, Λ)
    if m != nothing
        m -= 1
        @warn "$(length(Λ) - m) negative eigenvalues ≥ $(Λ[end]) detected, Gaussian
            random field will be approximated (ignoring all negative eigenvalues)"

        resize!(Λ, m)
        U = U[:, 1:m]
    end

    # store eigenvalues and eigenfunctions
    Λ .= sqrt.(Λ) # note: store sqrt of eigenval for more efficient sampling
    data = SpectralData(Λ, U)

    GaussianRandomField{Spectral,typeof(cov),typeof(pts),typeof(mean),typeof(data)}(mean, cov, pts, data)
end

# returns the required dimension of the random points
randdim(grf::GaussianRandomField{Spectral}) = length(grf.data.eigenval)

# see KarhunenLoeve.jl for sample function

Base.show(io::IO, ::Spectral) = print(io, "spectral decomposition")


text:
@testset "complex" begin
	# Just catch that it doesn't error for now
	x0 = randn(ComplexF64, 3)
	@test_nowarn OnceDifferentiable(x -> sum(abs2, x), x0)
	@test_nowarn  NLSolversBase.value!(OnceDifferentiable(x -> sum(abs2, x), x0), x0)
	@test_nowarn  NLSolversBase.gradient!(OnceDifferentiable(x -> sum(abs2, x), x0), x0)
	@test_nowarn  NLSolversBase.value_gradient!(OnceDifferentiable(x -> sum(abs2, x), x0), x0)
end


text:
module KNNGraphs

using DataStructures
using Distances
using LightGraphs
using SparseArrays

include("abstract.jl")
include("heap_edge.jl")
include("threaded_heap_graph.jl")
include("heap_graph.jl")
include("heap_utils.jl")

# export public interface
export ApproximateKNNGraph, HeapKNNGraph, LockHeapKNNGraph, HeapKNNGraphEdge
export knn_diameter, knn_matrices
export flag, weight, edge_indices, node_edge, node_edges, update_flag!

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 2dd190ae-3888-11eb-2e68-090f885617d6
using DataStructures

# ╔═╡ 47c147e6-3888-11eb-1574-cd2d6aea7b1c
file = open(f->read(f, String), "day_10.txt")

# ╔═╡ 019df412-388a-11eb-17c1-2bfe6376975b
#first match get command, second get value of comand
regex = r"(\d+)"

# ╔═╡ f875f352-3952-11eb-3d4a-979e6ccb162f
to_int(x) = parse(Int64,x)

# ╔═╡ 913adf24-3955-11eb-2cf3-5bbfeee310db
input =append!([0],[ to_int(x.captures[1]) for x in eachmatch(regex,file)])

# ╔═╡ 5959a5a6-3ae8-11eb-304e-93efdbb1449d
sort!(input)

# ╔═╡ 604c7afc-3ae8-11eb-0b3e-ad3c795d8d39
difference= input[2:end] - input[1:end-1]

# ╔═╡ f2da1cce-3ae8-11eb-3823-93b721546763
dif1 =sum(map(x -> (x==1), difference))

# ╔═╡ 0cb01a04-3ae9-11eb-1503-713cd2879c56
#your device has always a difference of 3 in the last slot, thats why the +1
dif3 = sum(map(x -> (x==3), difference)) +  1

# ╔═╡ 4980f462-3ae9-11eb-0258-6fb53fe670f1
dif1 * dif3

# ╔═╡ aab19fbe-3956-11eb-0e2e-4d1994109154
md"# Question 2"

# ╔═╡ 907408d6-3aea-11eb-2ff8-650fd94ee053
md"Consider the input already in order:
.we have to include the pairs with difference 3.
If not, we will make a 4+ jump.

So divide the input in continuos '1' difference

Calculate with dp for each group
"

# ╔═╡ 60c7816c-3af9-11eb-01c7-a9830ae2db73

# ╔═╡ 6a3b0a72-3aed-11eb-037d-cfdd3bb9eb2c
begin
	group_size = 0
	size_continuos_groups = []

	for x in difference

		group_size +=1

		if(x==3)

			if(group_size > 2)
				push!(size_continuos_groups,group_size)
			end
			group_size = 0
		end
	end

	push!(size_continuos_groups,group_size + 1)

end

# ╔═╡ 03af7394-3afa-11eb-1a11-c77bce80ccda
size_continuos_groups

# ╔═╡ e0c5fc98-3af1-11eb-3935-0b3da67afc71
input

# ╔═╡ c9c035fa-3af5-11eb-2372-a5b6156890ce
difference

# ╔═╡ 130b05ce-3a23-11eb-17bb-218a402f4383
number_of_solutions = Dict([(2,1), (3,2), (1,1), (0,1)])

# ╔═╡ 39c08774-3af1-11eb-3fa2-7f55584e157f
function find_solutions(n)
	if haskey(number_of_solutions, n)
		return number_of_solutions[n]
	else
		result = find_solutions(n-1) + find_solutions(n-2) + find_solutions(n-3)
		number_of_solutions[n] = result
		return result
	end
end

# ╔═╡ 552d820c-3af4-11eb-29a2-e3868c3f3ae5
begin
	awn2 = 1
	for x in size_continuos_groups
		awn2 *= find_solutions(x)
	end
end

# ╔═╡ 6c9fe808-3af4-11eb-25db-838ef16b38de
awn2

# ╔═╡ Cell order:
# ╠═2dd190ae-3888-11eb-2e68-090f885617d6
# ╠═47c147e6-3888-11eb-1574-cd2d6aea7b1c
# ╠═019df412-388a-11eb-17c1-2bfe6376975b
# ╠═f875f352-3952-11eb-3d4a-979e6ccb162f
# ╠═913adf24-3955-11eb-2cf3-5bbfeee310db
# ╠═5959a5a6-3ae8-11eb-304e-93efdbb1449d
# ╠═604c7afc-3ae8-11eb-0b3e-ad3c795d8d39
# ╠═f2da1cce-3ae8-11eb-3823-93b721546763
# ╠═0cb01a04-3ae9-11eb-1503-713cd2879c56
# ╠═4980f462-3ae9-11eb-0258-6fb53fe670f1
# ╠═aab19fbe-3956-11eb-0e2e-4d1994109154
# ╟─907408d6-3aea-11eb-2ff8-650fd94ee053
# ╠═60c7816c-3af9-11eb-01c7-a9830ae2db73
# ╠═6a3b0a72-3aed-11eb-037d-cfdd3bb9eb2c
# ╠═03af7394-3afa-11eb-1a11-c77bce80ccda
# ╠═e0c5fc98-3af1-11eb-3935-0b3da67afc71
# ╠═c9c035fa-3af5-11eb-2372-a5b6156890ce
# ╠═130b05ce-3a23-11eb-17bb-218a402f4383
# ╠═39c08774-3af1-11eb-3fa2-7f55584e157f
# ╠═552d820c-3af4-11eb-29a2-e3868c3f3ae5
# ╠═6c9fe808-3af4-11eb-25db-838ef16b38de


text:
#!/bin/bash
#=
exec julia --project=. -e 'include(popfirst!(ARGS))' \
    "${BASH_SOURCE[0]}" "$@"
=#

# Script to extract twix headers into a new temporary directory

using ArgParse

argdef = ArgParseSettings(exc_handler=isinteractive() ? ArgParse.debug_handler : ArgParse.default_handler)
@add_arg_table argdef begin
    "in_file"
        help = "path to the input twix file (Siemens meas_*.dat RAID format)"
        required = true
    "out_dir"
        help = "path to the output directory"
        required = true
end

parsed_args = parse_args(ARGS, argdef)

using MagneticResonanceSignals
out_dir = parsed_args["out_dir"]
if !isdir(out_dir)
    mkpath(out_dir)
end
MagneticResonanceSignals.dump_twix_headers(parsed_args["in_file"], out_dir)

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
# Here we solve the equation on a unit square, with a uniform internal source.
# The strong form of the (linear) heat equation is given by
#
# ```math
#  -\nabla \cdot (k \nabla u) = f  \quad x \in \Omega,
# ```
#
# where $u$ is the unknown temperature field, $k$ the heat conductivity,
# $f$ the heat source and $\Omega$ the domain. For simplicity we set $f = 1$
# and $k = 1$. We will consider homogeneous Dirichlet boundary conditions such that
# ```math
# u(x) = 0 \quad x \in \partial \Omega,
# ```
# where $\partial \Omega$ denotes the boundary of $\Omega$.
#
# The HDG weak form is given by
# ```math
#(\vec{\sigma}_{h},\vec{v_{h}})_{\tau_{h}}-(u_{h},\nabla_{h}\cdot\vec{v}_{h})_{\tau_{h}}+\left\langle \hat{u}_{h},\vec{v}_{h}\cdot\vec{n}_{h}\right\rangle _{\partial\tau_{h}} & =0\quad\forall\vec{v}_{h}\in\vec{V}_{h},\\
#(w_{h},\nabla_{h}\cdot\vec{\sigma}_{h})_{\tau_{h}}+\left\langle (\hat{\sigma}_{h}-\vec{\sigma}_{h})\cdot\vec{n},w_{h}\right\rangle _{\partial\tau_{h}} & =(f,w_{h})_{\tau_{h}}\quad\forall w_{h}\in W_{h},\\
#\left\langle \vec{\sigma}_{h},\xi_{h}\right\rangle _{\partial\tau_{h}\backslash\Gamma} & =0\quad\forall\xi_{h}\in M_{h}^{o},\\
#\left\langle \hat{u}_{h},\xi_{h}\right\rangle _{\Gamma} & =\left\langle g,\xi_{h}\right\rangle _{\Gamma}\quad\forall\xi_{h}\in M_{h}^{\partial},\\
#\hat{\sigma}_{h}\cdot\vec{n} & =\vec{\sigma}_{h}\cdot\vec{n}+\tau(u_{h}-\hat{u}_{h}).
# ```
#-

using HDiscontinuousGalerkin
using Tensors
using BlockArrays
using LinearAlgebra

# Load mesh
#To load mesh from triangle
#root_file = "mesh/figure.1"
#@time mesh = parse_mesh_triangle(root_file)

#or use internal mesh
mesh = rectangle_mesh(TriangleCell, (10,10), Vec{2}((0.0,0.0)), Vec{2}((1.0,1.0)))

# ### Trial and test functions
dim = 2
finiteElement = GenericFiniteElement(Dubiner{dim,RefTetrahedron,1}())
Wh = ScalarFunctionSpace(mesh, finiteElement)
Vh = VectorFunctionSpace(mesh, finiteElement)
Mh = ScalarTraceFunctionSpace(Wh, GenericFiniteElement(Legendre{dim-1,RefTetrahedron,1}()))

# Variables
û_h = TrialFunction(Mh)
σ_h = TrialFunction(Vh)
u_h = TrialFunction(Wh)

# ### Boundary conditions
dbc = Dirichlet(û_h, mesh, "boundary", x -> 0)

# RHS function
f(x::Vec{dim}) = 2*π^2*sin(π*x[1])*sin(π*x[2])
# ### Assembling the linear system
# Now we have all the pieces needed to assemble the linear system, $K û = f$.
function doassemble(Vh, Wh, Mh, τ = 1.0)
    # Allocate Matrices
    n_basefuncs = getnbasefunctions(Vh)
    n_basefuncs_s = getnbasefunctions(Wh)
    n_basefuncs_t = getnbasefunctions(Mh)
    Ae = zeros(n_basefuncs, n_basefuncs)
    Be = zeros(n_basefuncs, n_basefuncs_s)
    Ce = zeros(n_basefuncs_s, n_basefuncs_s)
    Ee = zeros(n_basefuncs,3*n_basefuncs_t)
    Fe = zeros(n_basefuncs_s, 3*n_basefuncs_t)
    be = zeros(n_basefuncs_s,1)
    He = zeros(3*n_basefuncs_t, 3*n_basefuncs_t)

    # create a matrix assembler and rhs vector
    assembler = start_assemble(getnfaces(mesh)*n_basefuncs_t)
    rhs = zeros(Float64, getnfaces(mesh)*n_basefuncs_t)
    K_element = Array{AbstractMatrix{Float64}}(undef, getncells(mesh))
    b_element = Array{AbstractVector{Float64}}(undef, getncells(mesh))

    @inbounds for (cell_idx, cell) in enumerate(CellIterator(mesh))
        reinit!(Wh, cell)
        reinit!(Vh, cell)
        fill!(Ae, 0)
        fill!(Be, 0)
        fill!(Ce, 0)
        fill!(Ee, 0)
        fill!(Fe, 0)
        fill!(He, 0)
        fill!(be, 0)
        #Cell integrals
        for q_point in 1:getnquadpoints(Vh)
            dΩ = getdetJdV(Vh, q_point)
            for i in 1:n_basefuncs
                vh  = shape_value(Vh, q_point, i)
                div_vh = shape_divergence(Vh, q_point, i)
                for j in 1:n_basefuncs
                    σ = shape_value(Vh, q_point, j)
                    # Integral σ⋅v dΩ
                    Ae[i, j] += (σ ⋅ vh) * dΩ
                end
                for j in 1:n_basefuncs_s
                    u = shape_value(Wh, q_point, j)
                    # Integral u*∇⋅v dΩ
                    Be[i,j] += (u*div_vh) * dΩ
                end
            end
        end
        #RHS
        for q_point in 1:getnquadpoints(Wh)
            dΩ = getdetJdV(Wh, q_point)
            for i in 1:n_basefuncs_s
                w  = shape_value(Wh, q_point, i)
                fh = function_value(f, Wh, cell, q_point)
                # Integral f*u dΩ
                be[i] += fh*w*dΩ
            end
        end
        #Face integrals
        for face_idx in 1:getnfaces(cell)
            for q_point in 1:getnfacequadpoints(Wh)
                dS = getfacedetJdS(Wh, face_idx, q_point)
                orientation = face_orientation(cell, face_idx)
                for i in 1:n_basefuncs_s
                    w = shape_value(Wh, face_idx, q_point, i)
                    for j in 1:n_basefuncs_s
                        u = shape_value(Wh, face_idx, q_point, j)
                        # Integral_∂T τ*u ⋅ w dS
                        Ce[i,j] += τ*(u*w) * dS
                    end
                    w = shape_value(Wh, face_idx, q_point, i, orientation)
                    for j in 1:n_basefuncs_t
                        û = shape_value(Mh, q_point, j)
                        # Integral_∂T τ*û*w  dS
                        Fe[i,n_basefuncs_t*(face_idx-1)+j] += (τ*(û*w)) * dS
                    end
                end
                dS = getfacedetJdS(Mh, face_idx, q_point)
                n = get_normal(Vh, face_idx)
                for i in 1:n_basefuncs
                    v = shape_value(Vh, face_idx, q_point, i, orientation)
                    for j in 1:n_basefuncs_t
                        û = shape_value(Mh, q_point, j)
                        # Integral_∂T û(v⋅n)  dS
                        Ee[i,n_basefuncs_t*(face_idx-1)+j] += (û*(v⋅n)) * dS
                    end
                end
                for i in 1:n_basefuncs_t
                    μ = shape_value(Mh, q_point, i)
                    for j in 1:n_basefuncs_t
                        û = shape_value(Mh, q_point, j)
                        # Integral_∂T û*μ  dS
                        He[n_basefuncs_t*(face_idx-1)+i,n_basefuncs_t*(face_idx-1)+j] += (û*μ) * dS
                    end
                end
            end
        end
        #Assamble Ke
        Me = BlockArray{Float64}(undef, [n_basefuncs,n_basefuncs_s], [n_basefuncs,n_basefuncs_s])
        setblock!(Me, Ae, 1, 1)
        setblock!(Me, Be', 2, 1)
        setblock!(Me, -Be, 1, 2)
        setblock!(Me, Ce, 2, 2)
        Mei = factorize(Array(Me))
        EFe = BlockArray{Float64}(undef, [n_basefuncs,n_basefuncs_s], [3*n_basefuncs_t])
        setblock!(EFe, -Ee, 1, 1)
        setblock!(EFe, Fe, 2, 1)
        Ge = BlockArray{Float64}(undef, [n_basefuncs,n_basefuncs_s], [3*n_basefuncs_t])
        setblock!(Ge, Ee, 1, 1)
        setblock!(Ge, Fe, 2, 1)
        Ge = Ge'
        Bte = BlockArray{Float64}(undef, [n_basefuncs,n_basefuncs_s],[1])
        setblock!(Bte, zeros(n_basefuncs,1),1,1)
        setblock!(Bte, be,2,1)
        K_element[cell_idx]=Mei\EFe
        Ate = Ge*K_element[cell_idx]-He
        b_element[cell_idx]=(Mei\Bte)[:,1]
        bte = -Ge*b_element[cell_idx]
        #Assemble
        gdof = Vector{Int}()
        for fidx in mesh.cells[cell_idx].faces
            for j in 1:n_basefuncs_t
                push!(gdof, fidx*n_basefuncs_t-(n_basefuncs_t-j))
            end
        end
        assemble!(assembler, gdof, Ate)
        assemble!(rhs, gdof, bte[:,1])
    end
    return end_assemble(assembler), rhs, K_element, b_element
end

### Solution of the system
K, b, K_e, b_e = doassemble(Vh,Wh,Mh);

# To account for the boundary conditions we use the `apply!` function.
# This modifies elements in `K` and `f` respectively, such that
# we can get the correct solution vector `u` by using `\`.
apply!(K,b,dbc)
û = K \ b;
#Now we recover original variables from skeleton û
function get_uσ!(σ_h,u_h,û_h,û, K_e, b_e, mesh)
    n_cells = getncells(mesh)
    n_basefuncs = getnbasefunctions(σ_h)
    n_basefuncs_t = getnbasefunctions(Mh)
    for cell_idx in 1:getncells(mesh)
        #get dofs
        û_e = Vector{Float64}()
        for (k,face) in enumerate(mesh.cells[cell_idx].faces)
            push!(û_e, û[n_basefuncs_t*face-1:n_basefuncs_t*face]...)
            û_h.m_values[cell_idx,:,k] = û[n_basefuncs_t*face-1:n_basefuncs_t*face]
        end
        dof = K_e[cell_idx]*û_e+b_e[cell_idx]
        σ_h.m_values[cell_idx,:] = dof[1:n_basefuncs]
        u_h.m_values[cell_idx,:] = dof[n_basefuncs+1:end]
    end
end

get_uσ!(σ_h, u_h,û_h,û, K_e, b_e, mesh)
#Compute errors
u_ex(x::Vec{dim}) = sin(π*x[1])*sin(π*x[2])
Etu_h = errornorm(u_h, u_ex)
Etu_h <= 0.00006

### Plot Solution
# TODO: Turn into a recipe
using Makie, CairoMakie
#AbstractPlotting.__init__()
using ColorSchemes
connectivity = getcells_matrix(mesh)
coordinates = get_vertices_matrix(mesh)

# We avg since u_h is discontinuous
color = nodal_avg(u_h)
scene = Scene(resolution = (400, 200), colorrange = (0.0, 1.0), colormap = ColorSchemes.viridis.colors)
#poly!(scene, coordinates, connectivity, color = color, strokecolor = (:black, 0.6), strokewidth = 1)
for row in 1:size(connectivity)[1]
	#read coordinates
	points = node(:poly, Point2f0[coordinates[i,:] for i in connectivity[row,:]])
	colors = get(ColorSchemes.viridis, [color[i] for i in connectivity[row,:]])
    #colors = get(ColorSchemes.viridis, [value(u_h, node, row) for node in connectivity[row,:]])
	poly!(scene, points, strokewidth = 1, strokecolor = (:black,0.6), color = colors, show_axis = true, scale_plot = false)
end
display(scene)


text:
mpow2(a::AbstractArray) = a .^ 2

"""
    svd_back(U, S, V, dU, dS, dV)

adjoint for SVD decomposition.

References:
    https://j-towns.github.io/papers/svd-derivative.pdf
    https://giggleliu.github.io/2019/04/02/einsumbp.html
"""
function svd_back(U::AbstractArray, S::AbstractArray{T}, V, dU, dS, dV; η::Real=1e-40) where T
    all(x -> x isa Nothing, (dU, dS, dV)) && return nothing
    η = T(η)
    NS = length(S)
    S2 = mpow2(S)
    Sinv = @. S/(S2+η)
    F = S2' .- S2
    F ./= (mpow2(F) .+ η)

    res = ZeroAdder()
    if !(dU isa Nothing)
        UdU = U'*dU
        J = F.*(UdU)
        res += (J+J')*Diagonal(S) + Diagonal(1im*imag(diag(UdU)) .* Sinv)
    end
    if !(dV isa Nothing)
        VdV = V'*dV
        K = F.*(VdV)
        res += Diagonal(S) * (K+K')
    end
    if !(dS isa Nothing)
        res += Diagonal(dS)
    end

    res = U*res*V'

    if !(dU isa Nothing) && size(U, 1) != size(U, 2)
        res += (dU - U* (U'*dU)) * Diagonal(Sinv) * V'
    end

    if !(dV isa Nothing) && size(V, 1) != size(V, 2)
        res = res + U * Diagonal(Sinv) * (dV' - (dV'*V)*V')
    end
    res
end


text:
using Base.Meta
using BenchmarkTools

struct Variable
    index::Int
end

# Represents sum(coefficients[i]*variables[i] for i in ...) + offset
mutable struct AffineExpression
    variables::Vector{Variable}
    coefficients::Vector{Float64}
    offset::Float64
end

Base.zero(::Type{AffineExpression}) = AffineExpression(Variable[],Float64[],0.0)
Base.:*(c::Float64,v::Variable) = AffineExpression([v],[c],0.0)
function Base.:+(aff1::AffineExpression,aff2::AffineExpression)
    return AffineExpression([aff1.variables;aff2.variables],[aff1.coefficients;aff2.coefficients],aff1.offset+aff2.offset)
end

function build_generator(n)
    return sum(i^1.5*Variable(i) for i in 1:n)
end

function build_loop(n)
    aff = zero(AffineExpression)
    for i in 1:n
        aff += i^1.5*Variable(i)
    end
    return aff
end

function build_manual(n)
    aff = zero(AffineExpression)
    sizehint!(aff.variables, n)
    sizehint!(aff.coefficients, n)
    for i in 1:n
        push!(aff.variables, Variable(i))
        push!(aff.coefficients, i^1.5)
    end
    return aff
end

@enum NodeClass CONST VAR MUL ADD

struct OperationNode
    nc::NodeClass
    children::Vector{OperationNode}
    metadata::Float64
end

OperationNode(v::Variable) = OperationNode(VAR, [], v.index)
OperationNode(n::Float64) = OperationNode(CONST, [], n)

Base.:*(o1::OperationNode,o2::OperationNode) = OperationNode(MUL, [o1,o2], 0.0)
Base.:+(o1::OperationNode,o2::OperationNode) = OperationNode(ADD, [o1,o2], 0.0)

# special structure:
#          ADD
#         /   \
#       MUL    ADD
#      /   \     \
#  CONST   VAR   ...

function graph_to_aff(g::OperationNode, aff::AffineExpression = zero(AffineExpression))
    if g.nc == ADD
        for c in g.children
            graph_to_aff(c, aff)
        end
    elseif g.nc == MUL
        # we know it is CONST*VAR
        coef = 0.0
        varidx = 0
        @assert length(g.children) == 2
        for c in g.children
            if c.nc == CONST
                coef = c.metadata
            elseif c.nc == VAR
                varidx = Int(c.metadata)
            end
        end
        push!(aff.variables, Variable(varidx))
        push!(aff.coefficients, coef)
    else
        error()
    end
    return aff
end

function build_graph(n)

    g = sum(OperationNode(i^1.5)*OperationNode(Variable(i)) for i in 1:n)

    return graph_to_aff(g)
end

@benchmark build_generator(100)
@benchmark build_loop(100)
@benchmark build_manual(100)
@benchmark build_graph(100)

@benchmark build_generator(1000)
@benchmark build_loop(1000)
@benchmark build_manual(1000)
@benchmark build_graph(1000)

ex = :(sum(i^1.5*Variable(i) for i in 1:n))
dump(ex)
ex.args[2]
ex.args[2].args

function rewrite_generator(ex)
    @assert isexpr(ex,:call)
    @assert ex.args[1] == :sum
    @assert isexpr(ex.args[2], :generator)
    return quote
        aff = zero(AffineExpression)
        len = length($(ex.args[2].args[2].args[2]))
        sizehint!(aff.variables, len)
        sizehint!(aff.coefficients, len)
        for $(ex.args[2].args[2].args[1]) in $(ex.args[2].args[2].args[2])
            push!(aff.variables, $(ex.args[2].args[1].args[3]))
            push!(aff.coefficients, $(ex.args[2].args[1].args[2]))
        end
        aff
    end
end

rewrite_generator(:(sum(i^1.5*Variable(i) for i in 1:n)))

macro build_expression(ex)
    return rewrite_generator(ex)
end

function build_macro(n)
    return @build_expression sum(i^1.5*Variable(i) for i in 1:n)
end

@benchmark build_macro(100)
@benchmark build_macro(1000)


text:
"""
Partial out variables in a Dataframe

### Arguments
* `df`: A table
* `formula::FormulaTerm`: A formula created using `@formula`
* `add_mean::Bool`: Should the initial mean added to the returned variable?
* `method::Symbol`: A symbol for the method. Default is :cpu. Alternatively,  :gpu requires `CuArrays`. In this case, use the option `double_precision = false` to use `Float32`.
* `maxiter::Integer`: Maximum number of iterations
* `double_precision::Bool`: Should the demeaning operation use Float64 rather than Float32? Default to true.
* `tol::Real`: Tolerance
* `align::Bool`: Should the returned DataFrame align with the original DataFrame in case of missing values? Default to true

### Returns
* `::DataFrame`: a dataframe with as many columns as there are dependent variables and as many rows as the original dataframe.
* `::Vector{Int}`: a vector of iterations for each column
* `::Vector{Bool}`: a vector of success for each column

### Details
`partial_out` returns the residuals of a set of variables after regressing them on a set of regressors. The syntax is similar to `reg` - but it accepts multiple dependent variables. It returns a dataframe with as many columns as there are dependent variables and as many rows as the original dataframe.
The regression model is estimated on only the rows where *none* of the dependent variables is missing. With the option `add_mean = true`, the mean of the initial variable is added to the residuals.

### Examples
```julia
using  RDatasets, DataFrames, FixedEffectModels, Gadfly
df = dataset("datasets", "iris")
result = partial_out(df, @formula(SepalWidth + SepalLength ~ fe(Species)), add_mean = true)
plot(layer(result[1], x="SepalWidth", y="SepalLength", Stat.binmean(n=10), Geom.point),
   layer(result[1], x="SepalWidth", y="SepalLength", Geom.smooth(method=:lm)))
```
"""
function partial_out(
    @nospecialize(df),
    @nospecialize(f::FormulaTerm);
    @nospecialize(weights::Union{Symbol, Expr, Nothing} = nothing),
    @nospecialize(add_mean = false),
    @nospecialize(maxiter::Integer = 10000),
    @nospecialize(contrasts::Dict = Dict{Symbol, Any}()),
    @nospecialize(method::Symbol = :cpu),
    @nospecialize(double_precision::Bool = true),
    @nospecialize(tol::Real = double_precision ? 1e-8 : 1e-6),
    @nospecialize(align = true))

    if  (ConstantTerm(0) ∉ eachterm(f.rhs)) & (ConstantTerm(1) ∉ eachterm(f.rhs))
        f = FormulaTerm(f.lhs, tuple(ConstantTerm(1), eachterm(f.rhs)...))
    end
    formula, formula_endo, formula_iv = parse_iv(f)
    has_iv = formula_iv != nothing
    has_iv && throw("partial_out does not support instrumental variables")
    has_weights = weights != nothing

    # create a dataframe without missing values & negative weights
    all_vars = StatsModels.termvars(formula)
    esample = completecases(df[!, all_vars])
    if has_weights
        esample .&= BitArray(!ismissing(x) & (x > 0) for x in df[!, weights])
    end

    # initialize iterations & converged
    iterations = Int[]
    convergeds = Bool[]

    # Build fixedeffects, an array of AbtractFixedEffects
    fes, ids, ids_fes, formula = parse_fixedeffect(df, formula)
    has_fes = !isempty(fes)

    nobs = sum(esample)
    (nobs > 0) || throw("sample is empty")
    # Compute weights
    if has_weights
        weights = Weights(convert(Vector{Float64}, view(df, esample, weights)))
        all(isfinite, weights) || throw("Weights are not finite")
    else
        weights = uweights(sum(esample))
    end

    if has_fes
        # in case some FixedEffect does not have interaction, remove the intercept
        if any(isa(fe.interaction, UnitWeights) for fe in fes)
            formula = FormulaTerm(formula.lhs, tuple(ConstantTerm(0), (t for t in eachterm(formula.rhs) if t!= ConstantTerm(1))...))
            has_fes_intercept = true
        end
        fes = FixedEffect[fe[esample] for fe in fes]
        feM = AbstractFixedEffectSolver{double_precision ? Float64 : Float32}(fes, weights, Val{method})
    end

    # Compute residualized Y
    vars = unique(StatsModels.termvars(formula))
    subdf = Tables.columntable(disallowmissing!(df[esample, vars]))
    formula_y = FormulaTerm(ConstantTerm(0), (ConstantTerm(0), eachterm(formula.lhs)...))
    formula_y_schema = apply_schema(formula_y, schema(formula_y, subdf, contrasts), StatisticalModel)
    Y = convert(Matrix{Float64}, modelmatrix(formula_y_schema, subdf))

    ynames = coefnames(formula_y_schema)[2]
    if !isa(ynames, Vector)
        ynames = [ynames]
    end
    if add_mean
        m = mean(Y, weights, dims = 1)
    end
    if has_fes
        Y, b, c = solve_residuals!(Y, feM; maxiter = maxiter, tol = tol, progress_bar = false)
        append!(iterations, b)
        append!(convergeds, c)
    end

    # Compute residualized X
    formula_x = FormulaTerm(ConstantTerm(0), formula.rhs)
    formula_x_schema = apply_schema(formula_x, schema(formula_x, subdf, contrasts), StatisticalModel)
    X = convert(Matrix{Float64}, modelmatrix(formula_x_schema, subdf))
    if has_fes
        X, b, c = solve_residuals!(X, feM; maxiter = maxiter, tol = tol, progress_bar = false)
        append!(iterations, b)
        append!(convergeds, c)
    end
    if has_weights
        Y .= Y .* sqrt.(weights)
        X .= X .* sqrt.(weights)
    end
    # Compute residuals
    if size(X, 2) > 0
        residuals = Y .- X * (X \ Y)
    else
        residuals = Y
    end

    # rescale residuals
    if has_weights
        residuals .= residuals ./ sqrt.(weights)
    end
    if add_mean
        residuals .= residuals .+ m
    end

    # Return a dataframe
    out = DataFrame()
    j = 0

    for y in ynames
        j += 1
        if align & (nobs < length(esample))
            out[!, Symbol(y)] = Vector{Union{Float64, Missing}}(missing, size(df, 1))
            out[esample, Symbol(y)] = residuals[:, j]
        else
            out[!, Symbol(y)] = residuals[:, j]
        end
    end
    return out, esample, iterations, convergeds
end


text:
#!/usr/bin/env julia
# current command for Stripe 82 Experiment
# julia run_celeste_on_field_mcmc.jl \
#   --ais-output-dir ais-output-s82 \
#   --use-full-initialization \
#   --initialization-catalog ~/Proj/Celeste.jl/benchmark/accuracy/output/sdss_4263_5_119_primary_b97a8fda22.csv \
#   --target-source-range $SOURCE_RANGE
import Celeste.ArgumentParse
parser = Celeste.ArgumentParse.ArgumentParser()
ArgumentParse.add_argument(parser, "--joint", help="Use joint inference", action=:store_true)
ArgumentParse.add_argument(
    parser,
    "--use-full-initialization",
    help="Use all information from initialization catalog (otherwise will just use noisy position",
    action=:store_true,
)
ArgumentParse.add_argument(
    parser,
    "--target-source-range",
    help="Target only the given number of sources, for quicker testing",
    default="5:10",
)
ArgumentParse.add_argument(
    parser,
    "--images-jld",
    help="FITS file containing synthetic imagery; if not specified, will use SDSS (primary) images",
)
ArgumentParse.add_argument(
    parser,
    "--initialization-catalog",
    help="CSV catalog for initialization. Default is to initialize by source detection.",
)
ArgumentParse.add_argument(
    parser,
    "--ais-output-dir",
    help="location to save AIS+MCMC samples",
)
ArgumentParse.add_argument(
    parser,
    "--num_procs",
    help="number of processers to use --- a node has 64"
)
parsed_args = ArgumentParse.parse_args(parser, ARGS)
println(parsed_args)

# Set up AIS-MCMC sample output directory
ais_output_dir = "ais-output"
if haskey(parsed_args, "ais-output-dir")
    ais_output_dir = parsed_args["ais-output-dir"]
end
if !isdir(ais_output_dir)
    mkdir(ais_output_dir)
end

# handle multiple cores (alternatively --- run `julia -p <n_procs>`)
num_procs = 63
#if haskey(parsed_args, "num-procs")
#    num_procs = parse(Int, parsed_args["num_procs"])
#end
println("------ using ", num_procs, " cores -----------")
addprocs(num_procs)

# this insane syntax broadcasts the local variable parsed_args to every
# every process that we added above
@eval @everywhere parsed_args = $parsed_args
@eval @everywhere ais_output_dir = $ais_output_dir

@everywhere begin

############################
# do imports everywhere    #
############################
using DataFrames
import JLD
import Celeste
import Celeste: Config, Model, MCMC
import Celeste.ParallelRun
import Celeste.SDSSIO
import Celeste.Log
import Celeste.AccuracyBenchmark
const OUTPUT_DIRECTORY = joinpath(splitdir(Base.source_path())[1], "output")

###################################################################
# Parse input images --- if none, assume reference STRIPE82_RCF   #
###################################################################
if haskey(parsed_args, "images-jld")
    images = JLD.load(parsed_args["images-jld"], "images")
    catalog_label = splitext(basename(parsed_args["images-jld"]))[1]
else
    dataset = SDSSIO.SDSSDataSet(AccuracyBenchmark.SDSS_DATA_DIR)
    strategy = SDSSIO.PlainFITSStrategy(AccuracyBenchmark.SDSS_DATA_DIR)
    images = SDSSIO.load_field_images(dataset, AccuracyBenchmark.STRIPE82_RCF)
    catalog_label = @sprintf("sdss_%s_%s_%s", rcf.run, rcf.camcol, rcf.field)
end
@assert length(images) == 5

###########################################################################
# parse the initialization catalog
#  For stripe 82 --- use full initialization --- for
#  For synthetic experiment, use initialization catalog, but do not use
#   full init
###########################################################################
if !haskey(parsed_args, "initialization-catalog")
    parsed_args["initialization-catalog"] = expanduser(
      "~/Proj/Celeste.jl/benchmark/accuracy/output/sdss_4263_5_119_primary_b97a8fda22.csv")
end
if haskey(parsed_args, "initialization-catalog")
    catalog_data = AccuracyBenchmark.read_catalog(parsed_args["initialization-catalog"])
    if parsed_args["use-full-initialization"]
        println("Using full initialization from ", parsed_args["initialization-catalog"])
    end
    catalog_entries = AccuracyBenchmark.make_initialization_catalog(
        catalog_data, parsed_args["use-full-initialization"])
    target_sources = collect(1:length(catalog_entries))
end

###########################################
# Generate sky patches for all sources    #
###########################################
function match_entry_to_patches(entry, detected_patches; detected_locs=nothing)
    if detected_locs == nothing
       detected_locs = hcat([p.world_center
                            for p in detected_patches[:, 1]]...)
    end
    dists = [ sum((detected_locs[:,i]-entry.pos).^2)
              for i in 1:size(detected_locs,2) ]
    min_dist, min_idx = findmin(dists)
    return detected_patches[min_idx, :], min_idx
end

function match_catalog_sources_to_detected_patches(catalog, detected_patches)
    # center of each detected patch
    detected_locs = hcat([p.world_center
                          for p in detected_patches[:, 1]]...)
    matched_idxs = []
    for entry in catalog
        _, m_idx = match_entry_to_patches(entry, detected_patches; detected_locs=detected_locs)
        push!(matched_idxs, m_idx)
    end
    return detected_patches[matched_idxs,:]
end

# if "use-full-initialization" flag is caught, then use the Model sky patches function
if parsed_args["use-full-initialization"]
    println("Using full initialization --- including Model Patches")
    catalog_patches = Model.get_sky_patches(images, catalog_entries)
else
    println("Not using full initialization, inferring patches from image")
    # detect patches with new function --- match to primary catalog to make comparison
    detected_catalog, detected_patches = Celeste.detect_sources(images)
    catalog_patches = match_catalog_sources_to_detected_patches(catalog_entries, detected_patches)
end

# generate neighboring source map
box = ParallelRun.BoundingBox(-1000.0, 1000.0, -1000.0, 1000.0)
entry_in_range = entry->((box.ramin < entry.pos[1] < box.ramax) &&
                         (box.decmin < entry.pos[2] < box.decmax))
target_ids = find(entry_in_range, catalog_entries)

# find objects with patches overlapping (for this specific patch)
neighbor_map = Dict(id=>Model.find_neighbors(catalog_patches, id)
                    for id in target_ids)

#######################################
# define process function for pmap    #
#######################################
function proc(ts)
    println("========== source ", ts, " of ", length(target_sources))
    # config = px_radius, num_temperatures, num_ais_samples, num_vi_iterations
    config = Config(25.0, 200, 25, 0)

    neighbor_ids = neighbor_map[ts]
    res = ParallelRun.process_source_mcmc(
        config, ts, catalog_entries, catalog_patches, neighbor_ids, images)

    # save source samples
    println(" ... saving source ", ts)
    fname = joinpath(ais_output_dir, "source-$(lpad(ts, 4, 0)).jld")
    JLD.save(fname, "res", res)
    return res
end

end ## everywhere end

#####################################################################
# Start Script                                                      #
#####################################################################
@printf("Loaded %d sources...\n", length(catalog_entries))

if haskey(parsed_args, "target-source-range")
    src_start, src_end = split(parsed_args["target-source-range"], ":")
    target_sources = parse(src_start):parse(src_end)
    println("====== inferring sources ", target_sources, " =========")
    target_sources = collect(target_sources)
end


text:
tests = ["conjugates",
         "conjugates_normal",
         "conjugates_mvnormal"]

for t in tests
    fpath = "$t.jl"
    println("running $fpath ...")
    include(fpath)
end


text:
function gen_rho( dr, σ1, σ2 )
    Npoints = size(dr)[1]
    rho = Array{Float64}(undef,Npoints)
    c1 = 2*σ1^2
    c2 = 2*σ2^2
    cc1 = sqrt(2*pi*σ1^2)^3
    cc2 = sqrt(2*pi*σ2^2)^3
    for ip=1:Npoints
        g1 = exp(-dr[ip]^2/c1)/cc1
        g2 = exp(-dr[ip]^2/c2)/cc2
        rho[ip] = g2 - g1
    end
    return rho
end