text:
using ArgParse
function inferDelim(sampleLines)
	typicalDelim = [',','\t','|',' ']
	delim = typicalDelim[end]
	modFlag = 0
	chars = unique(sampleLines[1]) #We don't care about all possible chars in the sample. The delimiter should be in the first line (and any line)
	tempCounts = zeros(Int8, length(chars), 1)
	countDict = Dict(zip(chars,tempCounts))
	for line in sampleLines
		countChars!(line, countDict)
	end
	potentialDelims = Char[]
	for char in keys(countDict)
 		if countDict[char] % length(sampleLines) == 0
			push!(potentialDelims, char)
		end
	end
	if length(potentialDelims) > 1
		for candidate in potentialDelims
			if candidate in typicalDelim && findfirst(typicalDelim, candidate) <= findfirst(typicalDelim, delim)
				delim = candidate
				modFlag = 1
			end
		end
	end
	if modFlag == 0
		print("Unable to successfully identify the correct delimiter.")
		showDelim(sampleLines, potentialDelims)
		print("Which delimiter looks most appropriate?")
		delim = readline(STDIN)
	end
	return delim
end
function showDelim(sampleLines, potentialDelims)
	stringCandidates = arrayAsEngString(potentialDelims)
	print("The potential delimiter candidates are $stringCandidates")
	for delim in potentialDelims
		print("The split for $delim in this file is")
		splitHeader = truncSepSamples(sampleLines[1], delim)
		splitLine1 = truncSepSamples(sampleLines[2], delim)
		splitLine2 = truncSepSamples(sampleLines[3], delim)
	end
end
function truncSepSamples(sampleLine, delim)
	splitLine = split(sampleLine, ",")
	truncLine = ""
	for piece in splitLine
		pieceLength = length(piece)
		if pieceLength > 5
			truncLine *= piece[1:5] * "..."
		else
			truncLine *= piece * (" " ^ (7-pieceLength))
		end
	end
	return truncLine
end
function arrayAsEngString(inputArray::Array)
	str = ""
	for i in eachindex(inputArray)
		if i != size(inputArray)[1]
			str *= " $(inputArray[i]),"
		else
			str *= " and $(inputArray[i])"
		end
	end
	return str
end
function countChars!(line::String, charCountDict::Dict)
	for char in line
		try charCountDict[char] += 1
		catch
			continue
		end
	end
end
function main(args)
		#initialize settings for argparse library
	opts = ArgParseSettings(description = "Script to create live objects from rows in csv")
	@add_arg_table opts begin
		"--delimiter", "-d"
			help = "Specifies the delimiter for the csv. If none is specified, one is inferred"
		"file"
			help = "Argument for the csv file to read in. Can be either positional or indicated using a flag"
			required = true
	end

	parsed_args = parse_args(opts)
	#Read first line of csv in so we can figure out what the headings should be
	if length(parsed_args) > 1
		open(parsed_args["file"]) do csvInput
			src = readlines(csvInput)
			header = src[1]
			if(haskey(parsed_args, "--delimiter"))
				delim = parsed_args["--delimiter"]
			elseif(haskey(parsed_args, "-d"))
				delim = parsed_args["-d"]
			else
				lenFile = length(src)
				if(lenFile > 10)
					sampleLines = src[1:5:end]
				else
					sampleLines = src[:]
				end
				delim = inferDelim(sampleLines)
				print(delim)
			end
		end
	end
end
main(ARGS)


text:
# ============================================================================ #
# Queue
# ============================================================================ #

abstract type AbstractQueue{T}
end

mutable struct Queue{T} <:AbstractQueue{T}
    head_ptr::Int64
    tail_ptr::Int64
    capacity::Int64
    data::Vector{T}
    function Queue{T}(start_size::Int) where {T}
        head_ptr  = 0
        tail_ptr  = 0
        min_size  = minimum_queue_required_size()
        capacity  = start_size < 0 ? min_size : max(min_size, start_size)
        data      = Vector{T}(undef, capacity)
        new(head_ptr, tail_ptr, capacity, data)
    end
end

function Queue{T}() where {T}
    min_size = minimum_queue_required_size()
    return Queue{T}(min_size)
end

function minimum_queue_required_size()
    return 2
end

function Base.length(queue::Queue{T}; head = queue.head_ptr, tail = queue.tail_ptr) where {T}
    if head == 0 || tail == 0
        return 0
    else
        h = head
        t = tail
        L = queue.capacity
        return t - h + 1 + ((t < h) * L)
    end
end

function Base.size(queue::Queue{T}) where {T}
    return (length(queue),)
end

function Base.isempty(queue::Queue{T}) where {T}
    return length(queue) <= 0
end

function capacity(queue::Queue{T}) where {T}
    return queue.capacity
end

function next_index(queue::Queue{T}, index::Int) where {T}
    n = index + 1
    if n > queue.capacity
        n = 1
    end
    return n
end

function prev_index(queue::Queue{T}, index::Int) where {T}
    n = index - 1
    if n < 1
        n = queue.capacity
    end
    return n
end

function is_capacity_reached(queue::Queue{T}, add_one::Bool = false) where {T}
    zero_capa                   = queue.capacity == 0
    L = length(queue)
    exceeding_length_actual     = L >= queue.capacity
    exceeding_length_pretending = L + (1 * add_one) > queue.capacity

    return zero_capa || exceeding_length_actual || exceeding_length_pretending
end

function grow_queue!(queue::Queue{T}) where {T}
    growth_factor   = 2
    L               = queue.capacity
    new_capa        = convert(Int, ceil(queue.capacity * growth_factor))
    queue.capacity  = queue.capacity == 0 ? 1 : new_capa
    new_data        = similar( queue.data, T, queue.capacity )

    if queue.head_ptr >= queue.tail_ptr
        # 2 copies
        copyto!(new_data, 1,                      queue.data, queue.head_ptr, L - queue.head_ptr + 1)
        copyto!(new_data, L - queue.head_ptr + 2, queue.data, 1,              queue.tail_ptr)
        queue.head_ptr = 1
        queue.tail_ptr = L
    else
        # one copy
        copyto!(new_data, queue.data)
        queue.head_ptr = 1
        queue.tail_ptr = L
    end

    queue.data = new_data;

    return queue
end

function Base.push!(queue::Queue{T}, new_element::T) where {T}

    t_next = next_index(queue, queue.tail_ptr)

    if is_capacity_reached(queue, true)
        grow_queue!(queue)
        queue.tail_ptr += 1
    else
        queue.tail_ptr = t_next
    end

    queue.data[queue.tail_ptr] = new_element

    # first push
    if queue.head_ptr == 0 && queue.tail_ptr != 0
        queue.head_ptr = queue.tail_ptr
    end

    return queue
end

function Base.pop!(queue::Queue{T}) where {T}
    L = length(queue)
    if L <= 0
        underflow_error("empty queue")
    end

    elem   = queue.data[queue.head_ptr]

    if L == 1
        # reset queue
        queue.head_ptr = 0
        queue.tail_ptr = 0
    else
        queue.head_ptr = next_index(queue, queue.head_ptr)
    end

    return elem
end

function Base.peek(queue::Queue{T}) where {T}
    return queue.data[queue.head_ptr]
end

function Base.peek(queue::Queue{T}, n::Int = 0)::T where {T}
    idx = queue.head_ptr + n
    if idx > queue.capacity
        over = idx - queue.capacity
        idx = over
    end
    return queue.data[idx]
end

function Base.getindex(queue::Queue{T}, idx::Int) where {T}
    return queue.data[idx]
end

function Base.getindex(queue::Queue{T}, I::Vararg{Int, N}) where {T, N}
    return queue.data[I]
end

function Base.eltype(::Queue{T}) where {T}
    return T
end

function copy(queue::Queue{T}) where {T}
    new_queue           = Queue{T}()
    new_queue.data      = copy(queue.data)
    new_queue.capacity  = queue.capacity
    new_queue.head_ptr  = queue.head_ptr
    new_queue.tail_ptr  = queue.tail_ptr

    return new_queue
end

function Base.show(io::IO, queue::Queue{T}) where {T}
    println("Queue{$(T)} with properties:")
    println("\tcapacity:\t$(queue.capacity)")
    println("\thead_ptr:\t$(queue.head_ptr)")
    println("\ttail_ptr:\t$(queue.tail_ptr)")

    print("\tdata:\t\tVector{$(T)}")

    if queue.head_ptr > queue.tail_ptr
        print( view(queue.data, queue.head_ptr:queue.capacity) )
        print( view(queue.data, 1:queue.tail_ptr) )
    else
        if queue.head_ptr == 0 || queue.tail_ptr == 0
        else
            print( view(queue.data, queue.head_ptr:queue.tail_ptr) )
        end
    end

    print("\n\n")
end

# not implemented for now
# function shrink_to_fit!(stack::Stack{T}) where {T}
# end
# function merge!(stack::Stack{T}, vec::AbstractArray{T})
# end

ext: jl

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
@doc raw"""
    StochasticGradientProblem <: Problem

A stochastic gradient problem consists of
* a `Manifold M`
* a(n optional) cost function ``f(x) = \displaystyle\sum_{i=1}^n f_i(x)
* an array of gradients, i.e. a function that returns and array or an array of functions
``\{\operatorname{grad}f_i\}_{i=1}^n``, where both variants can be given in the allocating
variant and the array of function may also be provided as mutating functions `(X,x) -> X`.

# Constructors
    StochasticGradientProblem(M::Manifold, gradF::Function;
        cost=Missing(), evaluation=AllocatingEvaluation()
    )
    StochasticGradientProblem(M::Manifold, gradF::AbstractVector{<:Function};
        cost=Missing(), evaluation=AllocatingEvaluation()
    )

Create a Stochastic gradient problem with an optional `cost` and the gradient either as one
function (returning an array) or a vector of functions.
"""
struct StochasticGradientProblem{T,MT<:Manifold,TCost,TGradient} <:
       AbstractGradientProblem{T}
    M::MT
    cost::TCost
    gradient!!::TGradient
end
function StochasticGradientProblem(
    M::TM,
    gradF!!::G;
    cost::Union{Function,Missing}=Missing(),
    evaluation::AbstractEvaluationType=AllocatingEvaluation(),
) where {TM<:Manifold,G}
    return StochasticGradientProblem{typeof(evaluation),TM,typeof(cost),G}(M, cost, gradF!!)
end
function StochasticGradientProblem(
    M::TM,
    gradF!!::AbstractVector{<:Function};
    cost::Union{Function,Missing}=Missing(),
    evaluation::AbstractEvaluationType=AllocatingEvaluation(),
) where {TM<:Manifold}
    return StochasticGradientProblem{typeof(evaluation),TM,typeof(cost),typeof(gradF!!)}(
        M, cost, gradF!!
    )
end

@doc raw"""
    get_gradients(P::StochasticGradientProblem, x)
    get_gradients!(P::StochasticGradientProblem, Y, x)

Evaluate all summands gradients ``\{\operatorname{grad}f_i\}_{i=1}^n`` at `x` (in place of `Y`).

Note that for the [`MutatingEvaluation`](@ref) based problem and a single function for the
stochastic gradient, the allocating variant is not available.
"""
function get_gradients(
    p::StochasticGradientProblem{AllocatingEvaluation,<:Manifold,TC,<:Function}, x
) where {TC}
    return p.gradient!!(p.M, x)
end
function get_gradients(
    p::StochasticGradientProblem{AllocatingEvaluation,<:Manifold,TC,<:AbstractVector}, x
) where {TC}
    return [grad_i(p.M, x) for grad_i in p.gradient!!]
end
function get_gradients!(
    p::StochasticGradientProblem{AllocatingEvaluation,<:Manifold,TC,<:Function}, X, x
) where {TC}
    copyto!(X, p.gradient!!(p.M, x))
    return X
end
function get_gradients!(
    p::StochasticGradientProblem{AllocatingEvaluation,<:Manifold,TC,<:AbstractVector}, X, x
) where {TC}
    copyto!(X, [grad_i(p.M, x) for grad_i in p.gradient!!])
    return X
end
function get_gradients(
    ::StochasticGradientProblem{MutatingEvaluation,<:Manifold,TC,<:Function}, ::Any
) where {TC}
    return error(
        "For a mutating function type stochastic gradient, the allocating variant is not possible.",
    )
end
function get_gradients(
    p::StochasticGradientProblem{MutatingEvaluation,<:Manifold,TC,<:AbstractVector}, x
) where {TC}
    X = [zero_tangent_vector(p.M, x) for _ in 1:length(p.gradient!!)]
    return get_gradients!(p, X, x)
end
function get_gradients!(
    p::StochasticGradientProblem{MutatingEvaluation,<:Manifold,TC,<:Function}, X, x
) where {TC}
    return p.gradient!!(p.M, X, x)
end
function get_gradients!(
    p::StochasticGradientProblem{MutatingEvaluation,<:Manifold,TC,<:AbstractVector}, X, x
) where {TC}
    for i in 1:length(p.gradient!!)
        p.gradient!![i](p.M, X[i], x)
    end
    return X
end

@doc raw"""
    get_gradient(p::StochasticGradientProblem, k, x)
    get_gradient!(p::StochasticGradientProblem, Y, k, x)

Evaluate one of the summands gradients ``\operatorname{grad}f_k``, ``k∈\{1,…,n\}``, at `x` (in place of `Y`).

Note that for the [`MutatingEvaluation`](@ref) based problem and a single function for the
stochastic gradient mutating variant is not available, since it would require too many allocatins.
"""
function get_gradient(
    p::StochasticGradientProblem{AllocatingEvaluation,<:Manifold,TC,<:Function}, k, x
) where {TC}
    return p.gradient!!(p.M, x)[k]
end
function get_gradient(
    p::StochasticGradientProblem{AllocatingEvaluation,<:Manifold,TC,<:AbstractVector}, k, x
) where {TC}
    return p.gradient!![k](p.M, x)
end
function get_gradient!(
    p::StochasticGradientProblem{AllocatingEvaluation,<:Manifold,TC,<:Function}, X, k, x
) where {TC}
    copyto!(X, p.gradient!!(p.M, x)[k])
    return X
end
function get_gradient!(
    p::StochasticGradientProblem{AllocatingEvaluation,<:Manifold,TC,<:AbstractVector},
    X,
    k,
    x,
) where {TC}
    copyto!(X, p.gradient!![k](p.M, x))
    return X
end
function get_gradient(
    p::StochasticGradientProblem{MutatingEvaluation,<:Manifold,TC}, k, x
) where {TC}
    X = zero_tangent_vector(p.M, x)
    return get_gradient!(p, X, k, x)
end
function get_gradient!(
    ::StochasticGradientProblem{MutatingEvaluation,<:Manifold,TC,<:Function},
    ::Any,
    ::Any,
    ::Any,
) where {TC}
    return error(
        "A mutating variant of the stochastic gradient as a single function is not implemented.",
    )
end
function get_gradient!(
    p::StochasticGradientProblem{MutatingEvaluation,<:Manifold,TC,<:AbstractVector}, X, k, x
) where {TC}
    return p.gradient!![k](p.M, X, x)
end

"""
    AbstractStochasticGradientDescentOptions <: Options

A generic type for all options related to stochastic gradient descent methods
"""
abstract type AbstractStochasticGradientProcessor <: DirectionUpdateRule end

"""
    StochasticGradientDescentOptions <: AbstractStochasticGradientDescentOptions

Store the following fields for a default stochastic gradient descent algorithm,
see also [`StochasticGradientProblem`](@ref) and [`stochastic_gradient_descent`](@ref).

# fields

# Fields
* `x` the current iterate
* `stopping_criterion` ([`StopAfterIteration`](@ref)`(1000)`)– a [`StoppingCriterion`](@ref)
* `stepsize` ([`ConstantStepsize`](@ref)`(1.0)`) a [`Stepsize`](@ref)
* `evaluation_order` – (`:Random`) – whether
  to use a randomly permuted sequence (`:FixedRandom`), a per
  cycle permuted sequence (`:Linear`) or the default `:Random` one.
* `order` the current permutation
* `retraction_method` – (`ExponentialRetraction()`) a `retraction(M,x,ξ)` to use.

# Constructor
    StochasticGradientDescentOptions(x)

Create a [`StochasticGradientDescentOptions`](@ref) with start point `x`.
all other fields are optional keyword arguments.
"""
mutable struct StochasticGradientDescentOptions{
    TX,
    TV,
    D<:DirectionUpdateRule,
    TStop<:StoppingCriterion,
    TStep<:Stepsize,
    RM<:AbstractRetractionMethod,
} <: AbstractGradientOptions
    x::TX
    gradient::TV
    direction::D
    stop::TStop
    stepsize::TStep
    order_type::Symbol
    order::Vector{<:Int}
    retraction_method::RM
    k::Int # current iterate
end
function StochasticGradientDescentOptions(
    x,
    X,
    direction::DirectionUpdateRule;
    order_type::Symbol=:RandomOrder,
    order::Vector{<:Int}=Int[],
    retraction_method::AbstractRetractionMethod=ExponentialRetraction(),
    stoping_criterion::StoppingCriterion=StopAfterIteration(1000),
    stepsize::Stepsize=ConstantStepsize(1.0),
)
    return StochasticGradientDescentOptions{
        typeof(x),
        typeof(X),
        typeof(direction),
        typeof(stoping_criterion),
        typeof(stepsize),
        typeof(retraction_method),
    }(
        x,
        X,
        direction,
        stoping_criterion,
        stepsize,
        order_type,
        order,
        retraction_method,
        0,
    )
end

"""
    StochasticGradient <: DirectionUpdateRule

The default gradient processor, which just evaluates the (stochastic) gradient or a subset
thereof.
"""
struct StochasticGradient{T} <: AbstractStochasticGradientProcessor
    dir::T
end

function (s::StochasticGradient)(
    p::StochasticGradientProblem, o::StochasticGradientDescentOptions, iter
)
    # for each new epoche choose new order if we are at random order
    ((o.k == 1) && (o.order_type == :Random)) && shuffle!(o.order)
    # i is the gradient to choose, either from the order or completely random
    j = o.order_type == :Random ? rand(1:length(o.order)) : o.order[o.k]
    return o.stepsize(p, o, iter), get_gradient!(p, s.dir, j, o.x)
end
function MomentumGradient(
    p::StochasticGradientProblem,
    x0::P,
    s::DirectionUpdateRule=StochasticGradient(zero_tangent_vector(p.M, x0));
    gradient=zero_tangent_vector(p.M, x0),
    momentum=0.2,
    vector_transport_method::VTM=ParallelTransport(),
) where {P,VTM<:AbstractVectorTransportMethod}
    return MomentumGradient{P,typeof(gradient),typeof(momentum),VTM}(
        deepcopy(x0), gradient, momentum, s, vector_transport_method
    )
end
function AverageGradient(
    p::StochasticGradientProblem,
    x0::P,
    n::Int=10,
    s::DirectionUpdateRule=StochasticGradient(zero_tangent_vector(p.M, x0));
    gradients=fill(zero_tangent_vector(p.M, x0), n),
    vector_transport_method::VTM=ParallelTransport(),
) where {P,VTM}
    return AverageGradient{P,eltype(gradients),VTM}(
        gradients, deepcopy(x0), s, vector_transport_method
    )
end


text:
# Function to create matrix from array of arrays. This version is faster than vcat
function fast_mat(S)
    X = Array{Float32, 2}(undef, length(S[1]), length(S))
    for i=1:length(S)
        X[:,i] = S[i]
    end
    X
end

# This function samples episodes from the policy using importance sampling to weight them
function sample_episodes(mdp, policy::Policy, Nsamps::Int64; max_steps::Int64 = 1000,
                         a_and_p = action_and_probability, ϵ = 0)
    S, G, avg_ep_return = [], Float32[], 0.
    Neps = 0
    while length(G) < Nsamps
        Neps += 1
        s, steps = rand(initialstate(mdp)), 0
        Si, Ri, ρi = [], [], []
        push!(Si, convert_s(Array{Float64,1}, s, mdp))
        try
            push!(Ri, reward(mdp, s))
        catch
            push!(Ri, 0.)
        end
        while !isterminal(mdp, s)
            if  rand() < ϵ
                a, prob = rand(actions(mdp)), 1. / length(actions(mdp))
            else
                a, prob = a_and_p(policy, s)
            end
            push!(ρi, policy.action_probability(mdp, s, a) / prob)
            s, r = gen(mdp, s, a)
            push!(Si, convert_s(Array{Float64,1}, s, mdp))
            push!(Ri, r)
            steps += 1
            steps >= max_steps && break
        end
        push!(ρi, 1.)
        steps >= max_steps && println("Episode timeout at ", max_steps, " steps")
        weighted_returns = reverse(cumsum(reverse(Ri))) .* reverse(cumprod(reverse(ρi)))
        push!(S, Si...)
        push!(G, weighted_returns...)
        avg_ep_return += sum(Ri)
    end
    S, G, avg_ep_return = fast_mat(S), G', avg_ep_return/Neps
    @assert size(S,2) == size(G,2)
    @assert size(G,1) == 1
    S, G, avg_ep_return
end


text:
load ("Mongrel2")
using Mongrel2

t = run_server("6DFF1523-C091-49B8-B635-598640E864B3", "tcp://127.0.0.1:9997", "tcp://127.0.0.1:9996")
while true
     (conn, req) = consume (t)
     response = "Sender: $(req.sender_id)
ConnectionId: $(req.connection_id)
Path: $(req.path)
Headers: $(string(req.headers))
 Body: $(req.body)"

     if is_disconnected(req);continue; end
     reply_http(conn, req, response); disconnect_client(conn,req);
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct ActivePowerPI <: ActivePowerControl
        Kp_p::Float64
        Ki_p::Float64
        ωz::Float64
        P_ref::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
    end

Parameters of a Proportional-Integral Active Power controller for a specified power reference

# Arguments
- `Kp_p::Float64`: Proportional Gain, validation range: `(0, nothing)`
- `Ki_p::Float64`: Integral Gain, validation range: `(0, nothing)`
- `ωz::Float64`: filter frequency cutoff, validation range: `(0, nothing)`
- `P_ref::Float64`: Reference Power Set-point, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states of the ActivePowerPI model are:
	σp_oc: Integrator state of the PI Controller,
	p_oc: Measured active power of the inverter model
- `n_states::Int`: ActivePowerPI has two states
"""
mutable struct ActivePowerPI <: ActivePowerControl
    "Proportional Gain"
    Kp_p::Float64
    "Integral Gain"
    Ki_p::Float64
    "filter frequency cutoff"
    ωz::Float64
    "Reference Power Set-point"
    P_ref::Float64
    ext::Dict{String, Any}
    "The states of the ActivePowerPI model are:
	σp_oc: Integrator state of the PI Controller,
	p_oc: Measured active power of the inverter model"
    states::Vector{Symbol}
    "ActivePowerPI has two states"
    n_states::Int
end

function ActivePowerPI(Kp_p, Ki_p, ωz, P_ref=1.0, ext=Dict{String, Any}(), )
    ActivePowerPI(Kp_p, Ki_p, ωz, P_ref, ext, [:σp_oc, :p_oc], 2, )
end

function ActivePowerPI(; Kp_p, Ki_p, ωz, P_ref=1.0, ext=Dict{String, Any}(), states=[:σp_oc, :p_oc], n_states=2, )
    ActivePowerPI(Kp_p, Ki_p, ωz, P_ref, ext, states, n_states, )
end

# Constructor for demo purposes; non-functional.
function ActivePowerPI(::Nothing)
    ActivePowerPI(;
        Kp_p=0,
        Ki_p=0,
        ωz=0,
        P_ref=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`ActivePowerPI`](@ref) `Kp_p`."""
get_Kp_p(value::ActivePowerPI) = value.Kp_p
"""Get [`ActivePowerPI`](@ref) `Ki_p`."""
get_Ki_p(value::ActivePowerPI) = value.Ki_p
"""Get [`ActivePowerPI`](@ref) `ωz`."""
get_ωz(value::ActivePowerPI) = value.ωz
"""Get [`ActivePowerPI`](@ref) `P_ref`."""
get_P_ref(value::ActivePowerPI) = value.P_ref
"""Get [`ActivePowerPI`](@ref) `ext`."""
get_ext(value::ActivePowerPI) = value.ext
"""Get [`ActivePowerPI`](@ref) `states`."""
get_states(value::ActivePowerPI) = value.states
"""Get [`ActivePowerPI`](@ref) `n_states`."""
get_n_states(value::ActivePowerPI) = value.n_states

"""Set [`ActivePowerPI`](@ref) `Kp_p`."""
set_Kp_p!(value::ActivePowerPI, val) = value.Kp_p = val
"""Set [`ActivePowerPI`](@ref) `Ki_p`."""
set_Ki_p!(value::ActivePowerPI, val) = value.Ki_p = val
"""Set [`ActivePowerPI`](@ref) `ωz`."""
set_ωz!(value::ActivePowerPI, val) = value.ωz = val
"""Set [`ActivePowerPI`](@ref) `P_ref`."""
set_P_ref!(value::ActivePowerPI, val) = value.P_ref = val
"""Set [`ActivePowerPI`](@ref) `ext`."""
set_ext!(value::ActivePowerPI, val) = value.ext = val


text:
include("../src/wiedemannchik.jl")

#------------------------------------------------------------------------------

using Plots
using Nemo

#------------------------------------------------------------------------------

sizes = 50 : 50 : 2500
ZZ = GF(2^31 - 1)

iters = 1:10

timeit(start) = (time_ns() - start) * 1e-9

function compare_det_and_rand()

    for _ in iters
        for (idx, n) in enumerate(sizes)
            println("n = ", n)
            nnz_per_row = 4
            nnz_in_b    = 1

            S = MatrixSpace(Nemo.ZZ, n, n)
            s = MatrixSpace(Nemo.ZZ, n, 1)

            while true
                dataA = zeros(Nemo.ZZ, n, n)
                dataB = zeros(Nemo.ZZ, n, 1)

                for m in 1:n
                    for z in 1:nnz_per_row
                        i, j = m, rand(1:n)
                        dataA[i, j] = rand(Nemo.ZZ, 1:2^31-1)
                    end
                    # so, we have nnz diagonal!
                    # the matrix will be nonsingular almost surely
                    dataA[m, m] = 1
                end

                for z in 1:nnz_in_b
                    i = rand(1:n)
                    dataB[i] = rand(Nemo.ZZ, 1:2^31-1)
                end

                dataB_ = reshape(dataB, (n))
                A = from_dense(dataA, ZZ)
                b = from_dense(dataB_, ZZ)
                A_ = S(dataA)
                b_ = s(dataB)

                try
                    println("W. rand")
                    start = time_ns()
                    @time y_W = wiedemann_solve(A, b, proved=false)
                    res_W_rand = timeit(start)

                    println("W. det")
                    start = time_ns()
                    @time y_W = wiedemann_solve(A, b, proved=true)
                    res_W_det = timeit(start)

                    times_W_rand[idx] += res_W_rand
                    times_W_det[idx]  += res_W_det

                    break

                catch e
                    if isa(e, InterruptException)
                        rethrow(e)
                    end
                    @info e
                end
            end
            println("-------------")
        end
    end
end

times_W_rand = zeros(length(sizes))
times_W_det  = zeros(length(sizes))

compare_det_and_rand()

for i in 1 : length(times_W_det)
    times_W_rand[i] = times_W_rand[i] / length(iters)
    times_W_det[i]  = times_W_det[i] / length(iters)
end

plot(sizes, times_W_det, label="W. det", ylabel="time, s",  xlabel="order, n", title="rand vs det, 4 nnz per row")
plot!(sizes, times_W_rand, label="W. rand")

#------------------------------------------------------------------------------

#=

magical code

using Nemo;
Base.one(gfp_elem) = RR(1)

Nemo.gfp_elem(x::gfp_elem) = RR(x)

Base.zero(gfp_elem) = RR(0)

import Base: /
/(x::gfp_elem, y::gfp_elem) = RR(x) // RR(y)

Base.abs(x::gfp_elem) = RR(x)

Base.isless(x::gfp_elem, y::gfp_elem) = x.data < y.data

RR = GF(7)

A = [
    RR(50) RR(-10);
    RR(3) RR(49);
]

b = [
    RR(5); RR(64);
]

y = A \ b
=#


text:
using FromFile
@from "Core.jl" import Options, Dataset, RecordType
@from "EquationUtils.jl" import stringTree, countNodes
@from "LossFunctions.jl" import scoreFunc
@from "MutationFunctions.jl" import genRandomTree
@from "PopMember.jl" import PopMember
# A list of members of the population, with easy constructors,
#  which allow for random generation of new populations
mutable struct Population{T<:Real}
    members::Array{PopMember{T}, 1}
    n::Int
end
"""
    Population(pop::Array{PopMember{T}, 1})

Create population from list of PopMembers.
"""
Population(pop::Array{PopMember{T}, 1}) where {T<:Real} = Population{T}(pop, size(pop, 1))
"""
    Population(dataset::Dataset{T}, baseline::T;
               npop::Int, nlength::Int=3, options::Options,
               nfeatures::Int)

Create random population and score them on the dataset.
"""
Population(dataset::Dataset{T}, baseline::T;
           npop::Int, nlength::Int=3,
           options::Options,
           nfeatures::Int) where {T<:Real} = Population([PopMember(dataset, baseline, genRandomTree(nlength, options, nfeatures), options) for i=1:npop], npop)
"""
    Population(X::AbstractMatrix{T}, y::AbstractVector{T},
               baseline::T; npop::Int, nlength::Int=3,
               options::Options, nfeatures::Int)

Create random population and score them on the dataset.
"""
Population(X::AbstractMatrix{T}, y::AbstractVector{T}, baseline::T;
           npop::Int, nlength::Int=3,
           options::Options,
           nfeatures::Int) where {T<:Real} = Population(Dataset(X, y), baseline, npop=npop, options=options, nfeatures=nfeatures)

# Sample 10 random members of the population, and make a new one
function samplePop(pop::Population, options::Options)::Population
    idx = rand(1:pop.n, options.ns)
    return Population(pop.members[idx])
end

# Sample the population, and get the best member from that sample
function bestOfSample(pop::Population, options::Options)::PopMember
    sample = samplePop(pop, options)
    if options.probPickFirst == 1.0
        best_idx = argmin([sample.members[member].score for member=1:options.ns])
        return sample.members[best_idx]
    else
        sort_idx = sortperm([sample.members[member].score for member=1:options.ns])
        # Lowest comes first
        k = range(0.0, stop=options.ns-1, step=1.0) |> collect
        p = options.probPickFirst

        # Weighted choice:
        prob_each = p * (1 - p) .^ k
        prob_each /= sum(prob_each)
        cumprob = cumsum(prob_each)
        chosen_idx = findfirst(cumprob .> rand(Float32))

        return sample.members[chosen_idx]
    end
end

function finalizeScores(dataset::Dataset{T},
                        baseline::T, pop::Population,
                        options::Options)::Population where {T<:Real}
    need_recalculate = options.batching
    if need_recalculate
        @inbounds @simd for member=1:pop.n
            pop.members[member].score = scoreFunc(dataset, baseline,
                                                  pop.members[member].tree,
                                                  options)
        end
    end
    return pop
end

# Return best 10 examples
function bestSubPop(pop::Population; topn::Int=10)::Population
    best_idx = sortperm([pop.members[member].score for member=1:pop.n])
    return Population(pop.members[best_idx[1:topn]])
end

function record_population(pop::Population{T}, options::Options)::RecordType where {T<:Real}
    RecordType("population"=>[RecordType("tree"=>stringTree(member.tree, options),
                                         "loss"=>member.score,
                                         "complexity"=>countNodes(member.tree),
                                         "birth"=>member.birth,
                                         "ref"=>member.ref,
                                         "parent"=>member.parent)
                             for member in pop.members],
               "time"=>time()
    )
end


text:
subclones = [
    SomaticEvolution.CloneTracker(
        1, 1, 1.4437821887595856, [2], 2, 1, 1.0, 9
    ),
    SomaticEvolution.CloneTracker(
        1, 1, 1.7344043713478965, [1, 4], 3, 2, 1.6666666666666667, 90
    )
]
N = 100
@test SomaticEvolution.getclonesize(N, subclones) == [1, 9, 90]

#check birth and death rates are calculated correctly
b, d, selection = 1, 0, [1, 2]
brates, drates = SomaticEvolution.set_branching_birthdeath_rates(b, d, selection)
@test brates == [1, 2, 3]
@test drates == [0, 0, 0]


text:
module LoggerTests

using VarianceReducedSG
using Random
using Test

function test_loggers()
	mt = MersenneTwister(741)

	x0 = randn(mt,1)
	n = 10

	grad!(out,x,i) = out .= -1.0
	vrg = VRGradient(grad!, [0.0], n)
	VarianceReducedSG.store_from_point!(vrg, x0)

	alg = SAGA(vrg,x0,1.0,MersenneTwister(800))
	alg.x .= x0

	f1(x,y) = x[1]
	f2(x,y) = length(y.y)

	f1_cache = CacheFuncVal(f1)
	f2_cache = CacheFuncVal(f2,Int)

	nlogs = 23
	niter = 97
	extralogs = 10
	logiters = zeros(Int,nlogs + extralogs)
	f1_vals = zeros(nlogs + extralogs)
	f2_vals = zeros(Int,nlogs + extralogs)
	x0_vals = zeros(Float64,nlogs + extralogs)

	loggers = (
		ShowNewLine("Begining of log:"),
		NoLog(),
		ShowIterations(),
		ShowTime(),
		ShowNewLine(),
		f1_cache,
		f2_cache,
		ShowFuncVal(f1_cache, "Iterate"),
		ShowFuncVal(f2_cache, "Number of functions"),
		ShowAlgState(alg->alg.stepsize, "Stepsize"),
		ShowAlgState(alg->alg.x0, "Initial point"),
		StoreAlgState(alg->alg.x0[], x0_vals),
		ShowNewLine("\nEnd of log.\n"),
		StoreFuncVal(f1_cache, f1_vals),
		StoreFuncVal(f2_cache, f2_vals),
		StoreLogIterations(logiters),
		)

	iterate!(alg, niter, loggers, nlogs, warmstart=true)

	loginterval = Int(ceil(niter/nlogs))
	correct_logiters = range(0,step=loginterval,length=nlogs)
	@test all(logiters[1:nlogs] .== correct_logiters)
	@test all(logiters[(nlogs+1):end] .== 0)
	@test all(isapprox.(f1_vals[1:nlogs], x0 .+ correct_logiters))
	@test all(isapprox.(f1_vals[nlogs+1:end], 0.0))
	@test all(f2_vals[1:nlogs] .== n)
	@test all(f2_vals[nlogs+1:end] .== 0)
	@test all(x0_vals[1:nlogs] .== x0)
	@test all(x0_vals[nlogs+1:end] .== 0)
end

test_loggers()

end


text:
# Copyright 2018 Martin Holters
# See accompanying license file.

import Compat
using Compat: @warn, BitSet, argmax, argmin, copyto!, dropdims, findall, rmul!, undef
using Compat.Markdown: @doc_str

if isdefined(Base, :NamedTuple) # 0.7.0-DEV.2738 to 0.7.0-DEV.3226
    kwargs_pairs(kwargs::NamedTuple) = pairs(kwargs)
end
kwargs_pairs(kwargs) = kwargs

@static if !isdefined(Compat.SparseArrays, :blockdiag) # prior to 0.7.0-DEV.4498
    const blockdiag = Compat.SparseArrays.blkdiag
else
    using Compat.SparseArrays: blockdiag
end


text:
is_constraint(::typeof(power_balance)) = true

function _power_balance!(
    F, v_m, v_a, pinj, qinj, ybus_re, ybus_im, pv, pq, ref, nbus,
)
    npv = length(pv)
    npq = length(pq)
    if isa(F, Array)
        device = KA.CPU()
        kernel! = residual_kernel!(device)
    else
        device = KA.CUDADevice()
        kernel! = residual_kernel!(device)
    end
    ev = kernel!(
        F, v_m, v_a,
        ybus_re.colptr, ybus_re.rowval,
        ybus_re.nzval, ybus_im.nzval,
        pinj, qinj, pv, pq, nbus,
        ndrange=npv+npq,
        dependencies=Event(device)
    )
    wait(ev)
end

function power_balance(polar::PolarForm, cons, vm, va, pbus, qbus)
    nbus = get(polar, PS.NumberOfBuses())
    ref = polar.indexing.index_ref
    pv = polar.indexing.index_pv
    pq = polar.indexing.index_pq
    ybus_re, ybus_im = get(polar.topology, PS.BusAdmittanceMatrix())

    fill!(cons, 0.0)
    _power_balance!(
        cons, vm, va, pbus, qbus,
        ybus_re, ybus_im,
        pv, pq, ref, nbus
    )
end

function power_balance(polar::PolarForm, cons, buffer::PolarNetworkState)
    Vm, Va, pbus, qbus = buffer.vmag, buffer.vang, buffer.pinj, buffer.qinj
    power_balance(polar, cons, Vm, Va, pbus, qbus)
end

function size_constraint(polar::PolarForm, ::typeof(power_balance))
    npv = PS.get(polar.network, PS.NumberOfPVBuses())
    npq = PS.get(polar.network, PS.NumberOfPQBuses())
    return 2 * npq + npv
end

function bounds(polar::PolarForm{T, IT, VT, MT}, ::typeof(power_balance)) where {T, IT, VT, MT}
    m = size_constraint(polar, power_balance)
    return xzeros(VT, m) , xzeros(VT, m)
end

# Adjoint
function adjoint!(
    polar::PolarForm,
    pbm::AutoDiff.TapeMemory{F, S, I},
    cons, ∂cons,
    vm, ∂vm,
    va, ∂va,
    pinj, ∂pinj,
) where {F<:typeof(power_balance), S, I}
    nbus = get(polar, PS.NumberOfBuses())
    ref = polar.indexing.index_ref
    pv = polar.indexing.index_pv
    pq = polar.indexing.index_pq
    qinj = polar.reactive_load
    ybus_re, ybus_im = get(polar.topology, PS.BusAdmittanceMatrix())

    fill!(pbm.intermediate.∂edge_vm_fr , 0.0)
    fill!(pbm.intermediate.∂edge_vm_to , 0.0)
    fill!(pbm.intermediate.∂edge_va_fr , 0.0)
    fill!(pbm.intermediate.∂edge_va_to , 0.0)

    adj_residual_polar!(
        cons, ∂cons,
        vm, ∂vm,
        va, ∂va,
        ybus_re, ybus_im, polar.topology.sortperm,
        pinj, ∂pinj, qinj,
        pbm.intermediate.∂edge_vm_fr,
        pbm.intermediate.∂edge_vm_to,
        pbm.intermediate.∂edge_va_fr,
        pbm.intermediate.∂edge_va_to,
        pv, pq, nbus,
    )
end

function matpower_jacobian(polar::PolarForm, X::Union{State, Control}, ::typeof(power_balance), V)
    nbus = get(polar, PS.NumberOfBuses())
    pf = polar.network
    ref = pf.ref ; nref = length(ref)
    pv = pf.pv ; npv = length(pv)
    pq = pf.pq ; npq = length(pq)
    Ybus = pf.Ybus

    dSbus_dVm, dSbus_dVa = PS.matpower_residual_jacobian(V, Ybus)

    if isa(X, State)
        j11 = real(dSbus_dVa[[pv; pq], [pv; pq]])
        j12 = real(dSbus_dVm[[pv; pq], pq])
        j21 = imag(dSbus_dVa[pq, [pv; pq]])
        j22 = imag(dSbus_dVm[pq, pq])
        return [j11 j12; j21 j22]
    elseif isa(X, Control)
        j11 = real(dSbus_dVm[[pv; pq], [ref; pv]])
        j12 = sparse(I, npv + npq, npv)
        j21 = imag(dSbus_dVm[pq, [ref; pv]])
        j22 = spzeros(npq, npv)
        return [j11 -j12; j21 j22]
    end
end

# Hessian
function matpower_hessian(
    polar::PolarForm,
    ::typeof(power_balance),
    buffer::PolarNetworkState,
    λ::AbstractVector,
)
    ref = polar.indexing.index_ref
    pv = polar.indexing.index_pv
    pq = polar.indexing.index_pq
    V = buffer.vmag .* exp.(im .* buffer.vang)
    hxx, hxu, huu = PS.residual_hessian(V, polar.network.Ybus, λ, pv, pq, ref)
    return FullSpaceHessian(hxx, hxu, huu)
end


text:
#!/bin/julia

# This file must be ran inside the build/ directory of the project. It is assumed that
# the upper directory i.e. ../ is the main project directory where the Project.toml and
# Manifest.toml files reside.
const PROJECT = "\"Garamond\""
const BUILD_DIR = "/" * joinpath(split(@__FILE__, "/")[1:end-1]...)       # root project directory
const PROJECT_DIR = abspath(joinpath(BUILD_DIR, ".."))
const BASE_ENV_PATH = joinpath(ENV["HOME"], ".julia", "environments",
                               "v" * string(Int(Base.VERSION.major)) * "." *
                                     string(Int(Base.VERSION.minor))
                              )
const COMPILED_DIR = "compiled"
const APPS_DIR = "apps"
const APPS_PATH = abspath(joinpath(PROJECT_DIR, APPS_DIR))
const COMPILED_PATH = abspath(joinpath(BUILD_DIR, COMPILED_DIR))        # build directory
const TARGETS = ["gars", "garc", "garw"]                                # files to be compiles
const TARGETS_PATH = abspath.(joinpath.(APPS_PATH, TARGETS))            # full path to the targes
const REQUIRED_PACKAGES = ["PackageCompiler"]                           # required packages for compilation
const PACKAGE_REVS = Dict(                                              # Packag revisions/branches/etc
    #"PackageCompiler"=>("https://github.com/JuliaLang/PackageCompiler.jl", "master"),
    "Languages"=>("https://github.com/JuliaText/Languages.jl", "master")
   )

##############
### Checks ###
##############

# Check that the upper directory has a Project.toml and the project is Garamond
project_file = joinpath(PROJECT_DIR, "Project.toml")
if isfile(project_file)
    local key_name, project_name
    try
        key_name, project_name = strip.(split(readline(project_file), "="))
    catch
        @warn "Something went wrong while parsing the first line of $project_file. Exiting..."
        exit()
    end
    if key_name != "name" || project_name != PROJECT
        @warn "Malformed first line of Project.toml (expected: name=$PROJECT). Exiting..."
        exit()
    end
else
    @warn "$project_file for $PROJECT does not exist. Exiting..."
    exit()
end

# Check for a manifest file, warn if not present
manifest_file = joinpath(PROJECT_DIR, "Manifest.toml")
if !isfile(manifest_file)
    @warn "$manifest_file for $PROJECT does not exist. Will continue..."
end

@info "Pre-checks complete."

# Build or cleanup directory
if !isdir(COMPILED_PATH)
    mkdir(COMPILED_PATH)
    @info "Created $COMPILED_PATH."
else
    rm(COMPILED_PATH, recursive=true, force=true)
    @info "Cleaned up $COMPILED_PATH."
end

####################
### Dependencies ###
####################

using Pkg
try
    # Add project depedencies to the required ones
    Pkg.activate(PROJECT_DIR)
    Pkg.instantiate()
    #for pkg in keys(Pkg.installed())
    #    push!(REQUIRED_PACKAGES, pkg)
    #end
catch e
    @warn "Could not read $PROJECT package dependencies! ($e)"
end

@info "Dependencies are: $REQUIRED_PACKAGES"
Pkg.activate(BASE_ENV_PATH)  # reactivate default environment
Pkg.update()
for pkg in REQUIRED_PACKAGES
    if  !(pkg in keys(Pkg.installed()))
        if pkg in keys(PACKAGE_REVS)
            # Add custom revision, branch etc
            url, rev = PACKAGE_REVS[pkg]
            Pkg.add(PackageSpec(url=url, rev=rev))
        else
            # Add registered version
            Pkg.add(pkg)
        end
    end
end
@info "Installed dependencies."

########################
### Build Executable ###
########################

if length(ARGS) != 0 && ARGS[1] == "--deps-only"
    @info "Skipping build (deps-only)."
else
    # Check that all targets exist
    for target in TARGETS_PATH
        if !isdir(target)
            @warn "Could not find $target. Exiting..."
            exit()
        end
    end
    @info "Checks complete."

    using PackageCompiler
    cd(APPS_PATH)
    Pkg.activate(PROJECT_DIR)
    for (i, target) in enumerate(TARGETS)
        @info "*** Building $(uppercase(target)) ***"
        create_app(target, joinpath(COMPILED_PATH, target), force=true)
    end
    @info "Build complete."
end


text:
@testset "LookupTable" begin
    table = FinancialDSL.Core.Compiler.OptimizingIR.LookupTable{FinancialDSL.Core.RiskFactor}()
    rf1 = FinancialDSL.Core.DiscountFactor(:PRE, Date(2019, 2, 1))
    @test rf1 ∉ table
    i = FinancialDSL.Core.Compiler.OptimizingIR.addentry!(table, rf1)
    @test rf1 ∈ table
    @test table[i] == rf1

    for rf in table
        @test rf == rf1
    end

    rf2 = FinancialDSL.Core.SpotCurrency(USD)
    @test rf2 ∉ table
    i = FinancialDSL.Core.Compiler.OptimizingIR.addentry!(table, rf2)
    @test rf2 ∈ table
    @test table[i] == rf2
    @test rf1 ∈ table

    for rf in table
        @test rf == rf1 || rf == rf2
    end

    let
        filtered_risk_factors = filter( x -> FinancialDSL.Core.risk_factor_symbol(x) === :USD, table)
        @test filtered_risk_factors[1] == FinancialDSL.Core.SpotCurrency(USD)
        @test length(filtered_risk_factors) == 1
        @test isa(filtered_risk_factors, AbstractVector)
    end
end

@testset "compile" begin

    pricing_date = Date(2018, 5, 29)
    currency_to_curves_map = Dict( "onshore" => Dict( :BRL => :PRE, :USD => :cpUSD ))
    static_model = FinancialDSL.Core.StaticHedgingModel(BRL, currency_to_curves_map)
    attr = FinancialDSL.Core.ContractAttributes("riskfree_curves" => "onshore", "carry_type" => "none")
    empty_provider = FinancialDSL.MarketData.EmptyMarketDataProvider()

    scenario = FinancialDSL.Core.FixedScenario()
    scenario[FinancialDSL.Core.SpotCurrency(USD)] = 3.0BRL
    scenario[FinancialDSL.Core.DiscountFactor(:PRE, Date(2019, 1, 2))] = 0.7
    scenario[FinancialDSL.Core.DiscountFactor(:cpUSD, Date(2019, 2, 1))] = 0.8

    let
        # Test for ReduceObs
        obs_scale = FinancialDSL.Core.ReduceObs(+, Float64)
        FinancialDSL.Core.push_reduce_obs!(obs_scale, FinancialDSL.Core.Konst(5.0))
        FinancialDSL.Core.push_reduce_obs!(obs_scale, FinancialDSL.Core.LiftObs2(*, FinancialDSL.Core.Konst(2.0), FinancialDSL.Core.Konst(3.0)))
        FinancialDSL.Core.push_reduce_obs!(obs_scale, FinancialDSL.Core.LiftObs(-, FinancialDSL.Core.Konst(1.0)))

        contract = FinancialDSL.Core.Scale(obs_scale, FinancialDSL.Core.Unit(FinancialDSL.Core.SpotCurrency(BRL)))
        pricer = FinancialDSL.Core.compile_pricer(empty_provider, pricing_date, static_model, contract, attr)
        p = FinancialDSL.Core.price(pricer, scenario)
        @test p ≈ 10.0
    end

    let
        # Test for ReduceObs with a single element
        obs_scale = FinancialDSL.Core.ReduceObs(+, Float64)
        FinancialDSL.Core.push_reduce_obs!(obs_scale, FinancialDSL.Core.Konst(1))
        contract = FinancialDSL.Core.Scale(obs_scale, FinancialDSL.Core.Unit(FinancialDSL.Core.SpotCurrency(BRL)))
        pricer = FinancialDSL.Core.compile_pricer(empty_provider, pricing_date, static_model, contract, attr)
        p = FinancialDSL.Core.price(pricer, scenario)
        @test p ≈ 1.0
    end

    let
        contract = FinancialDSL.Core.Both(
            FinancialDSL.Core.Both(FinancialDSL.Core.Amount(10.0, BRL), FinancialDSL.Core.Amount(10.0, BRL)),
            FinancialDSL.Core.Amount(10.0, USD))

        pricer = FinancialDSL.Core.compile_pricer(empty_provider, pricing_date, static_model, contract, attr)
        println("Printing pricer for contract")
        println(contract)
        println(pricer)
        p = FinancialDSL.Core.price(pricer, scenario)
        @test p ≈ 10.0 + 10.0 + 10.0*3.0
        @test FinancialDSL.Core.get_functional_currency(pricer) == BRL

        ex = FinancialDSL.Core.exposures(FinancialDSL.Core.DeltaNormalExposuresMethod(), pricer, scenario)
        @test ex[FinancialDSL.Core.SpotCurrency(USD)] ≈ 10.0*3.0
    end

    let
        contract = FinancialDSL.Core.WhenAt(Date(2019, 2, 1), FinancialDSL.Core.Amount(10.0, USD))

        pricer = FinancialDSL.Core.compile_pricer(empty_provider, pricing_date, static_model, contract, attr)
        p = FinancialDSL.Core.price(pricer, scenario)
        @test p ≈ 10.0 * 3.0 * 0.8

        ex = FinancialDSL.Core.exposures(FinancialDSL.Core.DeltaNormalExposuresMethod(), pricer, scenario)
        @test ex[FinancialDSL.Core.SpotCurrency(USD)] ≈ 10.0 * 3.0 * 0.8
        @test ex[FinancialDSL.Core.DiscountFactor(:cpUSD, Date(2019, 2, 1))] ≈ 10.0 * 3.0 * 0.8
    end

    let
        contract = FinancialDSL.Core.WhenAt(Date(2019, 2, 1), FinancialDSL.Core.Amount(10.0, USD))
        program = FinancialDSL.Core.compile_pricing_program(empty_provider, pricing_date, static_model, contract, attr)

        mem_buff = Vector{Any}(undef, 1_000_000)
        input_values_buff = Vector{Any}(undef, 10_000)

        pricer = FinancialDSL.Core.create_pricer(program, memory_buffer=mem_buff, input_values_buffer=input_values_buff)
        p = FinancialDSL.Core.price(pricer, scenario)
        @test p ≈ 10.0 * 3.0 * 0.8

        @test length(mem_buff) == 1_000_000
        @test length(input_values_buff) == 10_000

        ex = FinancialDSL.Core.exposures(FinancialDSL.Core.DeltaNormalExposuresMethod(), pricer, scenario)
        @test ex[FinancialDSL.Core.SpotCurrency(USD)] ≈ 10.0 * 3.0 * 0.8
        @test ex[FinancialDSL.Core.DiscountFactor(:cpUSD, Date(2019, 2, 1))] ≈ 10.0 * 3.0 * 0.8

        native_pricer = FinancialDSL.Core.create_pricer(program, compiler=:native)
        p = FinancialDSL.Core.price(native_pricer, scenario)
        @test p ≈ 10.0 * 3.0 * 0.8

        ex = FinancialDSL.Core.exposures(FinancialDSL.Core.DeltaNormalExposuresMethod(), native_pricer, scenario)
        @test ex[FinancialDSL.Core.SpotCurrency(USD)] ≈ 10.0 * 3.0 * 0.8
        @test ex[FinancialDSL.Core.DiscountFactor(:cpUSD, Date(2019, 2, 1))] ≈ 10.0 * 3.0 * 0.8
    end

    let
        contract = FinancialDSL.Core.WhenAt(Date(2019, 2, 1), FinancialDSL.Core.Amount(10.0, USD))
        program = FinancialDSL.Core.compile_pricing_program(empty_provider, pricing_date, static_model, contract, attr)

        mem_buff = Vector{Any}(undef, 1)
        input_values_buff = Vector{Any}()

        @test length(mem_buff) != FinancialDSL.OptimizingIR.required_memory_size(program.program)
        @test length(input_values_buff) != FinancialDSL.OptimizingIR.required_input_values_size(program.program)

        pricer = FinancialDSL.Core.create_pricer(program, memory_buffer=mem_buff, input_values_buffer=input_values_buff)
        p = FinancialDSL.Core.price(pricer, scenario)
        @test p ≈ 10.0 * 3.0 * 0.8

        @test length(mem_buff) == FinancialDSL.OptimizingIR.required_memory_size(program.program)
        @test length(input_values_buff) == FinancialDSL.OptimizingIR.required_input_values_size(program.program)

        ex = FinancialDSL.Core.exposures(FinancialDSL.Core.DeltaNormalExposuresMethod(), pricer, scenario)
        @test ex[FinancialDSL.Core.SpotCurrency(USD)] ≈ 10.0 * 3.0 * 0.8
        @test ex[FinancialDSL.Core.DiscountFactor(:cpUSD, Date(2019, 2, 1))] ≈ 10.0 * 3.0 * 0.8

        native_pricer = FinancialDSL.Core.create_pricer(program, compiler=:native)
        p = FinancialDSL.Core.price(native_pricer, scenario)
        @test p ≈ 10.0 * 3.0 * 0.8

        ex = FinancialDSL.Core.exposures(FinancialDSL.Core.DeltaNormalExposuresMethod(), native_pricer, scenario)
        @test ex[FinancialDSL.Core.SpotCurrency(USD)] ≈ 10.0 * 3.0 * 0.8
        @test ex[FinancialDSL.Core.DiscountFactor(:cpUSD, Date(2019, 2, 1))] ≈ 10.0 * 3.0 * 0.8
    end
end


text:
# validate a previous solve

# using Distributed
# addprocs(8)

using Revise

using DistributedFactorGraphs
using Caesar, RoME
# @everywhere using Caesar, RoME
using LCMCore, BotCoreLCMTypes
using Dates
using DataStructures
using Logging
using JSON2
using DSP

# few local application specific functions
include(joinpath(@__DIR__, "SandsharkUtils.jl"))
include(joinpath(@__DIR__, "Plotting.jl"))

logpath = "/media/dehann/temp2/caesar/2020-01-05T13:22:16.359/"

fg = initfg()
loadDFG(logpath*"fg_final", Main, fg)

dontMarginalizeVariablesAll!(fg)
plotLocalProduct(fg, :x533)

Juno.@enter plotLocalProduct(fg, :x533)

ls(fg, :x533)

# # draw plots in getLogPath(fg)
# include(joinpath(@__DIR__, "GenerateResults.jl"))

## BATCH SOLVE

# dontMarginalizeVariablesAll!(fg)
# tree, smt, hist = solveTree!(fg)

#


text:
using Plots
using LinearAlgebra
using SparseArrays
using Arpack
# 1. initialize rho, uvec, f_i and feq_i
# 2. Streaming: move f_i -> f*_i in the direction of evec_i
# 3. Compute macroscopic rho and uvec from f*_i using:
# rho(xvec,t) = sum from i=0:8 of f_i(xvec,t), and
# uvec(xvec,t) = 1/rho * sum from i=0:8 of c*f_i*evec_i
# where c = delta x/delta t is the lattice speed
# 4. compute feq_i using feq_i(xvec,t) = w_i*rho + rho*s_i(uvec(xvec,t))
# 5. calculate updated distribution function f_i = f*_i - 1/tau * (f*_i - feq_i)

# make u, e three dimensional to contain index and 2D vector info?

struct StreamTensor

    E :: SparseMatrixCSC
    SE :: SparseMatrixCSC
    S :: SparseMatrixCSC
    SW :: SparseMatrixCSC
    W :: SparseMatrixCSC
    NW :: SparseMatrixCSC
    N :: SparseMatrixCSC
    NE :: SparseMatrixCSC

end

function PeriodicStreamTensor(Nx, Ny)

    Fx = spdiagm(-1 => ones(Nx - 1))
    Fx[1,Nx] = 1
    Bx = spdiagm(1 => ones(Nx - 1))
    Bx[Nx,1] = 1

    Fy = spdiagm(-1 => ones(Ny - 1))
    Fy[1,Ny] = 1
    By = spdiagm(1 => ones(Ny - 1))
    By[Ny,1] = 1

    Ix = sparse(I, Nx, Nx)
    Iy = sparse(I, Ny, Ny)

    E = kron(Iy, Fx)
    SE = kron(By, Fx)
    S = kron(By, Ix)
    SW = kron(By, Bx)
    W = kron(Iy, Bx)
    NW = kron(Fy, Bx)
    N = kron(Fy, Ix)
    NE = kron(Fy, Fx)

    return StreamTensor(E, SE, S, SW, W, NW, N, NE)
end

function MirrorStreamTensor(Nx, Ny)

    Fx = spdiagm(-1 => ones(Nx - 1))
    Fx[1,1] = 1
    Bx = spdiagm(1 => ones(Nx - 1))
    Bx[Nx,Nx] = 1

    Fy = spdiagm(-1 => ones(Ny - 1))
    Fy[1,1] = 1
    By = spdiagm(1 => ones(Ny - 1))
    By[Ny,Ny] = 1

    Ix = sparse(I, Nx, Nx)
    Iy = sparse(I, Ny, Ny)

    E = kron(Iy, Fx)
    SE = kron(By, Fx)
    S = kron(By, Ix)
    SW = kron(By, Bx)
    W = kron(Iy, Bx)
    NW = kron(Fy, Bx)
    N = kron(Fy, Ix)
    NE = kron(Fy, Fx)

    return StreamTensor(E, SE, S, SW, W, NW, N, NE)
end

function Base.:*(S::StreamTensor, f::Matrix{Float64})

    out = zeros(size(f))

    out[:,1] = f[:,1]
    out[:,2] = S.E * f[:,2]
    out[:,3] = S.N * f[:,3]
    out[:,4] = S.W * f[:,4]
    out[:,5] = S.S * f[:,5]
    out[:,6] = S.NE * f[:,6]
    out[:,7] = S.NW * f[:,7]
    out[:,8] = S.SW * f[:,8]
    out[:,9] = S.SE * f[:,9]

    return out
end

function lattice_boltzmann_method(xres, yres, f, N, tau, deltax, deltat, count=0)
    # vector at each point with nine dimensions corresponding to each direction 0-8
    c = deltax/deltat

    ftransform = zeros(xres*yres, 9)
    for j=1:yres
        for i=1:xres
            ftransform[(j-1)*xres+i,:] = f[i,j,:]
        end
    end

    stensor = PeriodicStreamTensor(xres, yres)
    fnew = stensor * ftransform

    rho = zeros((xres*yres))
    for i=1:xres*yres
        rho[i] = sum(fnew[i,:])
    end

    u = zeros(xres*yres, 2)
    for i=1:xres*yres
        ux = (fnew[i,2] - fnew[i,4] + 1/sqrt(2) * (fnew[i,6] + fnew[i,9] - fnew[i,7] - fnew[i,8])) / rho[i]
        uy = (fnew[i,3] - fnew[i,5] + 1/sqrt(2) * (fnew[i,6] + fnew[i,8] - fnew[i,9] - fnew[i,8])) / rho[i]
        u[i,1] = ux
        u[i,2] = uy
    end

    unitvecs = zeros(9,2)
    unitvecs[:,1] = [0, c, 0, -c, 0, 1/sqrt(2)*c, -1/sqrt(2)*c, -1/sqrt(2)*c, 1/sqrt(2)*c]
    unitvecs[:,2] = [0, 0, c, 0, -c, 1/sqrt(2)*c, 1/sqrt(2)*c, -1/sqrt(2)*c, -1/sqrt(2)*c]

    w = zeros(9)
    w[1] = 4/9
    w[2] = 1/9
    w[3] = 1/9
    w[4] = 1/9
    w[5] = 1/9
    w[6] = 1/36
    w[7] = 1/36
    w[8] = 1/36
    w[9] = 1/36

    s = zeros((xres*yres, 9))
    usquared = zeros((xres*yres))
    for i=1:xres*yres
        usquared[i] = dot(u[i,:],u[i,:])
        s[i,:] = w[:] .* (3*(unitvecs*u[i,:])/c^2 .+ 9/2*((unitvecs*u[i,:])).^2/c^3 .- 3/2*usquared[i]/c^2)
    end

    feq = zeros((xres*yres, 9))
    for i=1:xres*yres
        feq[i,:] = w[:].*rho[i] .+ rho[i].*s[i,:]
    end

    fnext = zeros((xres*yres, 9))
    for i=1:xres*yres
        fnext[i,:] = fnew[i,:] .- (deltat/(tau+.5*deltat)).*(fnew[i,:] - feq[i,:])
    end

    """for i=1:xres
        for j=1:yres
            for k=1:9
                # fnext[i,j,k] = 1/(1 - 1/tau) * (fnew[i,j,k] - 1/tau*feq[i,j,k])
                # fnext[i,j,k] = fnew[i,j,k] - (1/tau)*(fnew[i,j,k]-feq[i,j,k])
                fnext[i,j,k] = fnew[i,j,k] - (deltat/(tau+.5*deltat))*(fnew[i,j,k]-feq[i,j,k])
            end
        end
    end"""

    fn = zeros(xres, yres, 9)
    for j=1:yres
        for i=1:xres
            fn[i,j,:] = fnext[(j-1)*xres+i,:]
        end
    end

    count = count + 1
    if count < N
        lattice_boltzmann_method(xres, yres, fn, N, tau, deltax, deltat, count)
    else
        return fn, rho
    end
end

xres = 100
yres = 100
# gaussian distribution in x and y?
fgauss = zeros(xres, yres, 9)
for i=1:xres
    for j=1:yres
        for k=1:9
            fgauss[i,j,k] = 1/(50^2*sqrt(2*pi)) * exp(-(i-xres/2)^2/(2*50^2)) * exp(-(j-yres/2)^2/(2*50^2))
        end
    end
end

fsine = zeros(xres, yres, 9)
for i=1:xres
    for j=1:yres
        for k=1:9
            fsine[i,j,k] = sin(i*pi/xres)*sin(j*pi/yres)
        end
    end
end

N = 100
tau = 10000
deltax = 1
deltat = 1
c = deltax/deltat

#=rhofull = zeros(xres*yres)
anim = @animate for i=1:N
    fsine, rho = lattice_boltzmann_method(xres, yres, fsine, 1, tau, deltax, deltat)
    fplot = zeros(xres, yres)
    if i==1
        rhofull = rho
    else
        rhofull = hcat(rhofull, rho)
    end
    rhot = zeros(xres, yres)
    for j=1:yres
        for i=1:xres
            rhot[i,j] = rho[(j-1)*xres+i]
            for k=1:9
                fplot[i,j] = fplot[i,j] + fsine[i,j,k]
            end
        end
    end
    println(fplot[Int32(xres/2),Int32(yres/2)])
    heatmap(rhot)
end
# println(size(rhofull))
p = size(rhofull)[2]
covarmatrix = 1/(p-1) * rhofull * transpose(rhofull)

"""evals = eigvals(covarmatrix)
evecs = eigvecs(covarmatrix)
# println(size(evals))
# println(abs(evals[1]))
absevals = zeros(size(evals)[1])
for i=1:size(evals)[1]
    absevals[i] = Float64(abs(evals[i]))
end
# println(absevals[1])
sortedevals = sort(absevals)
# print(sortedevals[23])
print(sortedevals)"""

evals, evecs = eigs(covarmatrix, nev=20, which=:LM)
println(evals)

evecstransform = zeros(xres, yres)
for j=1:yres
    for i=1:xres
        evecstransform[i,j] = evecs[(j-1)*xres+i,3]
    end
end

gif(anim, fps=10)
heatmap(evecstransform)=#


text:
export reconstruction_simple, reconstruction_multiEcho, reconstruction_multiCoil, reconstruction_multiCoilMultiEcho, reconstruction_lowRank, RecoParameters

"""
Performs iterative image reconstruction independently for the data of all coils,
contrasts and slices

# Arguments
* `acqData::AcquisitionData`            - AcquisitionData object
* `reconSize::NTuple{2,Int64}`              - size of image to reconstruct
* `reg::Regularization`                 - Regularization to be used
* `sparseTrafo::AbstractLinearOperator` - sparsifying transformation
* `weights::Vector{Vector{Complex{<:AbstractFloat}}}` - sampling density of the trajectories in acqData
* `solvername::String`                  - name of the solver to use
* (`normalize::Bool=false`)             - adjust regularization parameter according to the size of k-space data
* (`params::Dict{Symbol,Any}`)          - Dict with additional parameters
"""
function reconstruction_simple( acqData::AcquisitionData{T}
                              , reconSize::NTuple{D,Int64}
                              , reg::Vector{Regularization}
                              , sparseTrafo::Trafo
                              , weights::Vector{Vector{Complex{T}}}
                              , solvername::String
                              , normalize::Bool=false
                              , encodingOps=nothing
                              , params::Dict{Symbol,Any}=Dict{Symbol,Any}()) where {D, T <: AbstractFloat}

  encDims = ndims(trajectory(acqData))
  if encDims!=length(reconSize)
    error("reco-dimensionality $(length(reconSize)) and encoding-dimensionality $(encDims) do not match")
  end

  numContr, numChan, numSl = numContrasts(acqData), numChannels(acqData), numSlices(acqData)

  encParams = getEncodingOperatorParams(;params...)

  # set sparse trafo in reg
  reg[1].params[:sparseTrafo] = sparseTrafo

  # reconstruction
  Ireco = zeros(Complex{T}, prod(reconSize), numSl, numContr, numChan)
  @sync for k = 1:numSl
    Threads.@spawn begin
      if encodingOps!=nothing
        F = encodingOps[:,k]
      else
        F = encodingOps_simple(acqData, reconSize, slice=k; encParams...)
      end
      for j = 1:numContr
        W = WeightingOp(weights[j])
        for i = 1:numChan
          kdata = kData(acqData,j,i,k).* weights[j]
          solver = createLinearSolver(solvername, W∘F[j]; reg=reg, params...)

          I = solve(solver, kdata, startVector=get(params,:startVector,Complex{T}[]),
                                 solverInfo=get(params,:solverInfo,nothing))

          if isCircular( trajectory(acqData, j) )
            circularShutter!(reshape(I, reconSize), 1.0)
          end
          Ireco[:,k,j,i] = I
        end
      end
    end
  end

  Ireco = reshape(Ireco, volumeSize(reconSize, numSl)..., numContr, numChan)

  return makeAxisArray(Ireco, acqData)
end

"""
Performs a iterative image reconstruction jointly for all contrasts. Different slices and coil images
are reconstructed independently.

# Arguments
* `acqData::AcquisitionData`            - AcquisitionData object
* `reconSize::NTuple{2,Int64}`              - size of image to reconstruct
* `reg::Regularization`                 - Regularization to be used
* `sparseTrafo::AbstractLinearOperator` - sparsifying transformation
* `weights::Vector{Vector{Complex{<:AbstractFloat}}}` - sampling density of the trajectories in acqData
* `solvername::String`                  - name of the solver to use
* (`normalize::Bool=false`)             - adjust regularization parameter according to the size of k-space data
* (`params::Dict{Symbol,Any}`)          - Dict with additional parameters
"""
function reconstruction_multiEcho(acqData::AcquisitionData{Complex{T}}
                              , reconSize::NTuple{D,Int64}
                              , reg::Vector{Regularization}
                              , sparseTrafo::Trafo
                              , weights::Vector{Vector{Complex{T}}}
                              , solvername::String
                              , normalize::Bool=false
                              , encodingOps=nothing
                              , params::Dict{Symbol,Any}=Dict{Symbol,Any}()) where {D , T <: AbstractFloat}

  encDims = ndims(trajectory(acqData))
  if encDims!=length(reconSize)
    error("reco-dimensionality $(length(reconSize)) and encoding-dimensionality $(encDims) do not match")
  end

  numContr, numChan, numSl = numContrasts(acqData), numChannels(acqData), numSlices(acqData)
  encParams = getEncodingOperatorParams(;params...)

  # set sparse trafo in reg
  reg[1].params[:sparseTrafo] = diagOp( repeat([sparseTrafo],numContr)... )

  W = WeightingOp( vcat(weights...) )

  # reconstruction
  Ireco = zeros(Complex{T}, prod(reconSize)*numContr, numChan, numSl)
  @sync for i = 1:numSl
    Threads.@spawn begin
      if encodingOps != nothing
        F = encodingOps[i]
      else
        F = encodingOp_multiEcho(acqData, reconSize, slice=i; encParams...)
      end
      for j = 1:numChan
        kdata = multiEchoData(acqData, j, i) .* vcat(weights...)
        solver = createLinearSolver(solvername, W∘F; reg=reg, params...)

        Ireco[:,j,i] = solve(solver,kdata; params...)
        # TODO circular shutter
      end
    end
  end

  if encDims==2
    # 2d reconstruction
    Ireco = reshape(Ireco, reconSize[1], reconSize[2], numContr, numChan, numSl)
    Ireco = permutedims(Ireco, [1,2,5,3,4])
  else
    # 3d reconstruction
    Ireco = reshape(Ireco, reconSize[1], reconSize[2], reconSize[3], numContr, numChan)
  end

  return makeAxisArray(permutedims(Ireco,[1,2,5,3,4]), acqData)
end

"""
Performs a SENSE-type iterative image reconstruction. Different slices and contrasts images
are reconstructed independently.

# Arguments
* `acqData::AcquisitionData`            - AcquisitionData object
* `reconSize::NTuple{2,Int64}`              - size of image to reconstruct
* `reg::Regularization`                 - Regularization to be used
* `sparseTrafo::AbstractLinearOperator` - sparsifying transformation
* `weights::Vector{Vector{Complex{<:AbstractFloat}}}` - sampling density of the trajectories in acqData
* `solvername::String`                  - name of the solver to use
* `senseMaps::Array{Complex{<:AbstractFloat}}`        - coil sensitivities
* (`normalize::Bool=false`)             - adjust regularization parameter according to the size of k-space data
* (`params::Dict{Symbol,Any}`)          - Dict with additional parameters
"""
function reconstruction_multiCoil(acqData::AcquisitionData{T}
                              , reconSize::NTuple{D,Int64}
                              , reg::Vector{Regularization}
                              , sparseTrafo::Trafo
                              , weights::Vector{Vector{Complex{T}}}
                              , solvername::String
                              , senseMaps::Array{Complex{T}}
                              , normalize::Bool=false
                              , encodingOps=nothing
                              , params::Dict{Symbol,Any}=Dict{Symbol,Any}()) where {D , T}

  encDims = ndims(trajectory(acqData))
  if encDims!=length(reconSize)
    error("reco-dimensionality $(length(reconSize)) and encoding-dimensionality $(encDims) do not match")
  end

  numContr, numChan, numSl = numContrasts(acqData), numChannels(acqData), numSlices(acqData)
  encParams = getEncodingOperatorParams(;params...)

  # set sparse trafo in reg
  reg[1].params[:sparseTrafo] = sparseTrafo

  # solve optimization problem
  Ireco = zeros(Complex{T}, prod(reconSize), numSl, numContr, 1)
  @sync for k = 1:numSl
    Threads.@spawn begin
      if encodingOps != nothing
        E = encodingOps[:,k]
      else
        E = encodingOps_parallel(acqData, reconSize, senseMaps; slice=k, encParams...)
      end

      for j = 1:numContr
        W = WeightingOp(weights[j],numChan)
        kdata = multiCoilData(acqData, j, k) .* repeat(weights[j], numChan)

        EFull = ∘(W, E[j], isWeighting=true)
        solver = createLinearSolver(solvername, EFull; reg=reg, params...)
        I = solve(solver, kdata; params...)

        if isCircular( trajectory(acqData, j) )
          circularShutter!(reshape(I, reconSize), 1.0)
        end
        Ireco[:,k,j] = I
      end
    end
  end

  Ireco = reshape(Ireco, volumeSize(reconSize, numSl)..., numContr, 1)

  return makeAxisArray(Ireco, acqData)
end

"""
Performs a SENSE-type iterative image reconstruction which reconstructs all contrasts jointly.
Different slices are reconstructed independently.

# Arguments
* `acqData::AcquisitionData`            - AcquisitionData object
* `reconSize::NTuple{2,Int64}`              - size of image to reconstruct
* `reg::Regularization`                 - Regularization to be used
* `sparseTrafo::AbstractLinearOperator` - sparsifying transformation
* `weights::Vector{Vector{Complex{<:AbstractFloat}}}` - sampling density of the trajectories in acqData
* `solvername::String`                  - name of the solver to use
* `senseMaps::Array{Complex{<:AbstractFloat}}`        - coil sensitivities
* (`normalize::Bool=false`)             - adjust regularization parameter according to the size of k-space data
* (`params::Dict{Symbol,Any}`)          - Dict with additional parameters
"""
function reconstruction_multiCoilMultiEcho(acqData::AcquisitionData{T}
                              , reconSize::NTuple{D,Int64}
                              , reg::Vector{Regularization}
                              , sparseTrafo::Trafo
                              , weights::Vector{Vector{Complex{T}}}
                              , solvername::String
                              , senseMaps::Array{Complex{T}}
                              , normalize::Bool=false
                              , encodingOps=nothing
                              , params::Dict{Symbol,Any}=Dict{Symbol,Any}()) where {D, T}

  encDims = ndims(trajectory(acqData))
  if encDims!=length(reconSize)
    error("reco-dimensionality $(length(reconSize)) and encoding-dimensionality $(encDims) do not match")
  end

  numContr, numChan, numSl = numContrasts(acqData), numChannels(acqData), numSlices(acqData)
  encParams = getEncodingOperatorParams(;params...)

  # set sparse trafo in reg
  reg[1].params[:sparseTrafo] = diagOp( repeat([sparseTrafo],numContr)... )

  W = WeightingOp( vcat(weights...), numChan )

  Ireco = zeros(Complex{T}, prod(reconSize)*numContr, numSl)
  @sync for i = 1:numSl
    Threads.@spawn begin
      if encodingOps != nothing
        E = encodingOps[i]
      else
        E = encodingOp_multiEcho_parallel(acqData, reconSize, senseMaps; slice=i, encParams...)
      end

      kdata = multiCoilMultiEchoData(acqData, i) .* repeat(vcat(weights...), numChan)
      solver = createLinearSolver(solvername, W∘E; reg=reg, params...)

      Ireco[:,i] = solve(solver, kdata; params...)
    end
  end

  if encDims==2
    # 2d reconstruction
    Ireco = reshape(Ireco, reconSize[1], reconSize[2], numContr, numSl, 1)
    Ireco = permutedims(Ireco, [1,2,4,3,5])
  else
    # 3d reconstruction
    Ireco = reshape(Ireco, reconSize[1], reconSize[2], reconSize[3], numContr, 1)
  end

  return makeAxisArray(Ireco, acqData)
end


text:
"""
So, ReX = oscillatory integral + hyperbolic integral. I expanded the first integral ages ago, but it was not an arbitrary precision algorithm. I tested solving ReX with this old expansion for the oscillatory integral added to the hyperbolic integral solved by arbitrary precision quadgk, and it still diverged when digits beyond Float64 became important. Thus, the oscillatory integral needs to be rewritten in the arbitrary precision summation algorithm (the algorithm I used for BesselI-StruveL). I will do this here and also test it against brute forcing the integral with arb. prec. quadgk provided quadgk manages to solve this!
"""

include("bessel_minus_struve.jl")
import .bessel_minus_struve

using QuadGK
using Plots
using ArbNumerics
plotly()
Plots.PlotlyBackend()

function oscillatory_integral(Ω, β, α, v, w)

    R = (v^2 - w^2) / (w^2 * v)
    a = sqrt(β^2 / 4 + R * β * coth(β * v / 2))
    b = R * β / sinh(β * v / 2)

    coefficient = 2 * α * β^(3 / 2) * v^3 * sinh(Ω * β / 2) / (3 * √π * w^3 * sinh(β / 2))

    integrand(x) = (sin(Ω * x) * cos(x) / (x^2 + a^2 - b * cos(v * x))^(3 / 2))

    integral = QuadGK.quadgk(x -> integrand(x), 0, Inf)[1]
    @show(coefficient * integral)
    BigFloat(coefficient * integral)
end

function arb_binomial(x, y)
    x = ArbReal("$x")
    y = ArbReal("$y")
    one = ArbReal("1")
    ArbNumerics.gamma(x + one) / (ArbNumerics.gamma(y + one) * ArbNumerics.gamma(x - y + one))
end

function oscillatory_integral_expansion(Ω, β, α, v, w; prec = 64)

    # Initialise precision of ArbReal to prec.
    setextrabits(0)
    setprecision(ArbReal, prec + 8)

    Ω = ArbReal("$Ω")
    β = ArbReal("$β")
    α = ArbReal("$α")
    v = ArbReal("$v")
    w = ArbReal("$w")

    R = ArbReal((v^2 - w^2) / (w^2 * v))
    a = ArbReal(sqrt(β^2 / 4 + R * β * coth(β * v / 2)))
    b = ArbReal(R * β / sinh(β * v / 2))

    coefficient = ArbReal(α * β^(3/2) * v^3 / (3 * a * √π * w^3 * sinh(β / 2)))

    θ(z, n) = ArbReal(√π * ArbNumerics.gamma(-n - 1/2) * abs(z)^n * z * bessel_minus_struve.BesselI_minus_StruveL(n + 1, a * abs(z); prec = prec) / 2)

    s = ArbReal(sinh(Ω * β / 2))

    n = ArbReal("0")
    result = ArbReal("0.0")
    err = eps(ArbReal(0, bits = prec + 8))  # Machine accuracy of specified precision prec.

    while true

        previous_result = ArbReal("$result")

        if mod(n, 2) == 0
            even_term = ArbReal("0.0")
            @fastmath @inbounds @simd for z in [ArbReal(Ω + 1), ArbReal(Ω - 1)]
                even_term += ArbReal(arb_binomial(n, n/2) * s * θ(z, n) / 2)
            end
        else
            even_term = ArbReal("0.0")
        end

        k_term = ArbReal("0.0")
        @fastmath @inbounds @simd for k in ArbReal("0"):ArbReal(floor(n / 2 - 1 / 2))
            k_coeff = ArbReal(arb_binomial(n, k))
            @fastmath @inbounds @simd for z in [ArbReal(Ω + 1 + v * (n - 2 * k)), ArbReal(Ω - 1 + v * (n - 2 * k)), ArbReal(Ω + 1 - v * (n - 2 * k)), ArbReal(Ω - 1 - v * (n - 2 * k))]
                k_term += ArbReal(k_coeff * s * θ(z, n) / 2)
            end
        end

        term = ArbReal(arb_binomial(-3/2, n) * (-b / (4 * a))^n * (even_term
         + k_term))

        # Break loop if term smaller than accuracy of result.
        if abs(term) < err
            break
        end

        result += ArbReal(term)
        println("term: n = ", n, "\nterm value: ", coefficient * term, "\ncurrent result: ", coefficient * result, "\n")

        n += ArbReal("1")

        # Double precision if rounding error in result exceeds accuracy specified by prec.
        if ball(result)[2] > err
            setprecision(ArbReal, precision(result) + 8)

            println("Not precise enough. Required error < ", err, ". Increasing precision to ", precision(result) + 8, " bits.\n")

            Ω = ArbReal("$Ω")
            β = ArbReal("$β")
            α = ArbReal("$α")
            v = ArbReal("$v")
            w = ArbReal("$w")

            R = ArbReal((v^2 - w^2) / (w^2 * v))
            a = ArbReal(sqrt(β^2 / 4 + R * β * coth(β * v / 2)))
            b = ArbReal(R * β / sinh(β * v / 2))

            coefficient = ArbReal(α * β^(3/2) * v^3 / (3 * a * √π * w^3 * sinh(β / 2)))

            θ(z, n) = ArbReal(√π * ArbNumerics.gamma(-n - 1/2) * abs(z)^n * z * bessel_minus_struve.BesselI_minus_StruveL(n + 1, a * abs(z); prec = precision(result)) / 2)

            s = ArbReal(sinh(Ω * β / 2))

            n = ArbReal("$(n - 1)")
            result = ArbReal("$previous_result")
        end
    end

    result *= ArbReal(coefficient)
    println("Frequency: ", ArbReal(Ω, bits = prec + 8), ". Final result: ", ArbReal(result, bits = prec + 8))
    ArbReal(result, bits = prec + 8)
end

Ω_range = 0.01:0.5:20.01
β = 2.0
α = 7.0
v = 5.8
w = 1.6
# For Ω_range = 3.1:
# 9.828299 seconds (46.17 M allocations: 1.859 GiB, 12.89% gc time)
@time oscill_expansion = [oscillatory_integral_expansion(Ω, β, α, v, w; prec = 64) for Ω in Ω_range]
# oscill_expansion = fill(36.009791933858445279331705890513148915488)
# @show(oscill_expansion)
# 273.411207 seconds (931.57 M allocations: 32.002 GiB, 5.19% gc time)
@time oscill_integral = [oscillatory_integral(Ω, β, α, v, w) for Ω in Ω_range]

diff = oscill_expansion .- oscill_integral
@show(diff)
# oscill_integral = fill(-5.81034244862553995092)
# @show(oscill_integral)
# p = plot(Ω_range, abs.(oscill_expansion), yaxis = :log, label = "Oscillatory", xlabel = "Ω", ylabel = "ReX")
# plot!(Ω_range, abs.(oscill_integral), label = "Integral")
# display(p)


text:
using JobSchedulers
using Base.Threads

scheduler_start()
sleep(2)
scheduler_status()

job = Job(@task(begin
    sleep(2)
    println("highpriority")
end), name = "high_priority", priority = 0)
display(job)
submit!(job)

max_stars_repo_name: GTrunSec/DevSecOps-Cells-Lab

max_issues_repo_name: GTrunSec/DevSecOps-cells

max_forks_repo_name: GTrunSec/DevSecOps-Cells-Lab


text:
using Test, LinearAlgebra
using WaveProp
using WaveProp.Nystrom
using WaveProp.Geometry
using WaveProp.Integration
using WaveProp.Mesh

@testset "Greens identity test" begin
    # construct interior solution
    Geometry.clear!()
    pde  = Helmholtz(dim=2,k=1)
    xout = SVector(3,3)
    u    = (x)   -> SingleLayerKernel(pde)(xout,x)
    dudn = (x,n) -> DoubleLayerKernel(pde)(xout,x,n)
    geo = Circle()
    Ω   = Domain(geo)
    M   = meshgen(Ω,h=0.1)
    mesh = NystromMesh(view(M,boundary(Ω));order=5)
    γ₀u   = γ₀(u,mesh)
    γ₁u   = γ₁(dudn,mesh)
    𝐒     = SingleLayerOperator(pde,mesh)
    𝐃     = DoubleLayerOperator(pde,mesh)
    e0    = WaveProp.Nystrom.error_interior_green_identity(𝐒,𝐃,γ₀u,γ₁u)
    norm(e0,Inf)/norm(γ₀u,Inf)
    shand = Kress()
    q     = SingularQuadratureRule(GaussLegendre(8),shand)
    δ𝐒    = singular_weights(𝐒,q)
    δ𝐃    = singular_weights(𝐃,q)
    SS    = 𝐒 + δ𝐒
    DD    = 𝐃 + δ𝐃
    e0    = WaveProp.Nystrom.error_interior_green_identity(SS,DD,γ₀u,γ₁u)
    norm(e0,Inf)/norm(γ₀u,Inf)
    # δS    = singular_weights_dim(𝐒)
    # δD    = singular_weights_dim(𝐃)
    # Sfull = Matrix(𝐒) + δS
    # Dfull = Matrix(𝐃) + δD
    # e1 = WaveProp.Nystrom.error_interior_green_identity(Sfull,Dfull,γ₀u,γ₁u)
    # @test norm(e1,Inf) < norm(e0,Inf)
    # S,D = Nystrom.single_double_layer(pde,mesh)
    # e2 = WaveProp.Nystrom.error_interior_green_identity(S,D,γ₀u,γ₁u)
    # @test e1 ≈ e2
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
function format(package_name)
    ctx = Pkg.Types.Context()
    haskey(ctx.env.project.deps, package_name) || error("Unknown package $package_name.")
    pkg_uuid = ctx.env.project.deps[package_name]
    pkg_path = ctx.env.manifest[pkg_uuid].path
    pkg_path===nothing && error("Package must be deved to be formatted.")

    DocumentFormat.format(Path(pkg_path))
end


text:
using IOCapture
using Test, Random

# hasfield was added in Julia 1.2. This definition borrowed from Compat.jl (MIT)
# Note: this can not be inside the testset
(VERSION < v"1.2.0-DEV.272") && (hasfield(::Type{T}, name::Symbol) where T = Base.fieldindex(T, name, false) > 0)

hascolor(io) = VERSION >= v"1.6.0-DEV.481" && get(io, :color, false)
has_escapecodes(s) = occursin(r"\e\[[^m]*m", s)
strip_escapecodes(s) = replace(s, r"\e\[[^m]*m" => "")

# Callable object for testing
struct Foo
    x
end
(foo::Foo)() = println(foo.x)

@testset "IOCapture.jl" begin
    # Capturing standard output
    c = IOCapture.capture() do
        println("test")
    end
    @test !c.error
    @test c.output == "test\n"
    @test c.value === nothing
    @test c.backtrace isa Vector
    @test isempty(c.backtrace)

    # Capturing standard error
    c = IOCapture.capture() do
        println(stderr, "test")
    end
    @test !c.error
    @test c.output == "test\n"
    @test c.value === nothing
    @test c.backtrace isa Vector
    @test isempty(c.backtrace)

    # Return values
    c = IOCapture.capture() do
        println("test")
        return 42
    end
    @test !c.error
    @test c.output == "test\n"
    @test c.value === 42
    @test c.backtrace isa Vector
    @test isempty(c.backtrace)

    c = IOCapture.capture() do
        println("test")
        println(stderr, "test")
        return rand(5,5)
    end
    @test !c.error
    @test c.output == "test\ntest\n"
    @test c.value isa Matrix{Float64}
    @test c.backtrace isa Vector
    @test isempty(c.backtrace)

    # Callable objects
    c = IOCapture.capture(Foo("callable test"))
    @test !c.error
    @test c.output == "callable test\n"
    @test c.value === nothing

    # Colors get discarded
    c = IOCapture.capture() do
        printstyled("foo", color=:red)
    end
    @test !c.error
    @test c.output == "foo"
    @test c.value === nothing

    # Colors are preserved if it's supported
    c = IOCapture.capture(color=true) do
        printstyled("foo", color=:red)
    end
    @test !c.error
    if hascolor(stdout)
        @test c.output == "\e[31mfoo\e[39m"
    else
        @test c.output == "foo"
    end
    @test c.value === nothing

    # This test checks that deprecation warnings are captured correctly
    c = IOCapture.capture(color=true) do
        println("println")
        @info "@info"
        f() = (Base.depwarn("depwarn", :f); nothing)
        f()
    end
    @test !c.error
    @test c.value === nothing
    # The output is dependent on whether the user is running tests with deprecation
    # warnings enabled or not. To figure out whether that is the case or not, we can
    # look at the .depwarn field of the undocumented Base.JLOptions object.
    @test isdefined(Base, :JLOptions)
    @test hasfield(Base.JLOptions, :depwarn)
    if Base.JLOptions().depwarn == 0 # --depwarn=no, default on Julia >= 1.5
        @test has_escapecodes(c.output) === hascolor(stderr)
        @test strip_escapecodes(c.output) == "println\n[ Info: @info\n"
    else # --depwarn=yes
        @test has_escapecodes(c.output) === hascolor(stderr)
        output_nocol = strip_escapecodes(c.output)
        @test startswith(output_nocol, "println\n[ Info: @info\n┌ Warning: depwarn\n")
    end

    # Exceptions -- normally rethrown
    @test_throws ErrorException IOCapture.capture() do
        println("test")
        error("error")
        return 42
    end

    # .. but can be controlled with rethrow
    c = IOCapture.capture(rethrow=Union{}) do
        println("test")
        error("error")
        return 42
    end
    @test c.error
    @test c.output == "test\n"
    @test c.value isa ErrorException
    @test c.value.msg == "error"

    c = IOCapture.capture(rethrow=Union{}) do
        error("error")
        println("test")
        return 42
    end
    @test c.error
    @test c.output == ""
    @test c.value isa ErrorException
    @test c.value.msg == "error"

    # .. including interrupts
    c = IOCapture.capture(rethrow=Union{}) do
        println("test")
        throw(InterruptException())
        return 42
    end
    @test c.error
    @test c.output == "test\n"
    @test c.value isa InterruptException

    # .. or setting rethrow = InterruptException
    @test_throws InterruptException IOCapture.capture(rethrow=InterruptException) do
        println("test")
        throw(InterruptException())
        return 42
    end

    # .. or a union of exception types
    @test_throws DivideError IOCapture.capture(rethrow=Union{DivideError,InterruptException}) do
        println("test")
        div(1, 0)
        return 42
    end
    @test_throws InterruptException IOCapture.capture(rethrow=Union{DivideError,InterruptException}) do
        println("test")
        throw(InterruptException())
        return 42
    end

    # don't throw on errors that don't match rethrow
    c = IOCapture.capture(rethrow=Union{DivideError,InterruptException}) do
        println("test")
        three = "1" + "2"
        return 42
    end
    @test c.error
    @test c.output == "test\n"
    @test c.value isa MethodError

    # Invalid rethrow values
    @test_throws TypeError IOCapture.capture(()->nothing, rethrow=:foo)
    @test_throws TypeError IOCapture.capture(()->nothing, rethrow=42)
    @test_throws TypeError IOCapture.capture(()->nothing, rethrow=true)
    @test_throws TypeError IOCapture.capture(()->nothing, rethrow=false)

    # Make sure that IOCapture does not stall if we are printing _a lot_ of bytes into
    # stdout. X-ref: https://github.com/fredrikekre/Literate.jl/issues/138
    @testset "Buffer filling" begin
        for nrows = 2 .^ (0:20)
            c = IOCapture.capture() do
                for _ in 1:nrows; print("="^80); end
            end
            @test length(c.output) == 80 * nrows
        end
    end

    # Make sure the global rng isn't affected (JuliaLang/julia#41184).
    Random.seed!(1)
    r = rand()
    Random.seed!(1)
    c = IOCapture.capture(() -> rand())
    @test r == c.value
end


text:
using PiecewiseAffineTransforms
using MultivariateStats
using ActiveAppearanceModels
using Images
using ImageView
using Colors
using FaceDatasets
using Boltzmann
using HDF5
using JLD
using Clustering

const DATA_DIR = expanduser("~/data")
const DATA_DIR_PUT = joinpath(DATA_DIR, "PUT")
const DATA_DIR_CK = joinpath(DATA_DIR, "CK")
const MODEL_DIR = expanduser("~/data/models")
const RESIZERATIO = .25

include("view.jl")
include("utils.jl")


text:
import Base: *
import LinearAlgebra: mul!

abstract type FFTAPlan{T,N} <: Plan{T} end

struct FFTAInvPlan{T,N} <: FFTAPlan{T,N} end

struct FFTAPlan_cx{T,N} <: FFTAPlan{T,N}
    callgraph::NTuple{N, CallGraph{T}}
    region::Union{Int,AbstractVector{<:Int}}
    dir::Direction
    pinv::FFTAInvPlan{T}
end

struct FFTAPlan_re{T,N} <: FFTAPlan{T,N}
    callgraph::NTuple{N, CallGraph{T}}
    region::Union{Int,AbstractVector{<:Int}}
    dir::Direction
    pinv::FFTAInvPlan{T}
    flen::Int
end

function AbstractFFTs.plan_fft(x::AbstractArray{T}, region; kwargs...)::FFTAPlan_cx{T} where {T <: Complex}
    N = length(region)
    @assert N <= 2 "Only supports vectors and matrices"
    if N == 1
        g = CallGraph{T}(size(x,region[]))
        pinv = FFTAInvPlan{T,N}()
        return FFTAPlan_cx{T,N}((g,), region, FFT_FORWARD, pinv)
    else
        sort!(region)
        g1 = CallGraph{T}(size(x,region[1]))
        g2 = CallGraph{T}(size(x,region[2]))
        pinv = FFTAInvPlan{T,N}()
        return FFTAPlan_cx{T,N}((g1,g2), region, FFT_FORWARD, pinv)
    end
end

function AbstractFFTs.plan_bfft(x::AbstractArray{T}, region; kwargs...)::FFTAPlan_cx{T} where {T <: Complex}
    N = length(region)
    @assert N <= 2 "Only supports vectors and matrices"
    if N == 1
        g = CallGraph{T}(size(x,region[]))
        pinv = FFTAInvPlan{T,N}()
        return FFTAPlan_cx{T,N}((g,), region, FFT_BACKWARD, pinv)
    else
        sort!(region)
        g1 = CallGraph{T}(size(x,region[1]))
        g2 = CallGraph{T}(size(x,region[2]))
        pinv = FFTAInvPlan{T,N}()
        return FFTAPlan_cx{T,N}((g1,g2), region, FFT_BACKWARD, pinv)
    end
end

function AbstractFFTs.plan_rfft(x::AbstractArray{T}, region; kwargs...)::FFTAPlan_re{Complex{T}} where {T <: Real}
    N = length(region)
    @assert N <= 2 "Only supports vectors and matrices"
    if N == 1
        g = CallGraph{Complex{T}}(size(x,region[]))
        pinv = FFTAInvPlan{Complex{T},N}()
        return FFTAPlan_re{Complex{T},N}(tuple(g), region, FFT_FORWARD, pinv, size(x,region[]))
    else
        sort!(region)
        g1 = CallGraph{Complex{T}}(size(x,region[1]))
        g2 = CallGraph{Complex{T}}(size(x,region[2]))
        pinv = FFTAInvPlan{Complex{T},N}()
        return FFTAPlan_re{Complex{T},N}(tuple(g1,g2), region, FFT_FORWARD, pinv, size(x,region[1]))
    end
end

function AbstractFFTs.plan_brfft(x::AbstractArray{T}, len, region; kwargs...)::FFTAPlan_re{T} where {T}
    N = length(region)
    @assert N <= 2 "Only supports vectors and matrices"
    if N == 1
        g = CallGraph{T}(len)
        pinv = FFTAInvPlan{T,N}()
        return FFTAPlan_re{T,N}((g,), region, FFT_BACKWARD, pinv, len)
    else
        sort!(region)
        g1 = CallGraph{T}(len)
        g2 = CallGraph{T}(size(x,region[2]))
        pinv = FFTAInvPlan{T,N}()
        return FFTAPlan_re{T,N}((g1,g2), region, FFT_BACKWARD, pinv, len)
    end
end

function AbstractFFTs.plan_bfft(p::FFTAPlan_cx{T,N}) where {T,N}
    return FFTAPlan_cx{T,N}(p.callgraph, p.region, -p.dir, p.pinv)
end

function AbstractFFTs.plan_brfft(p::FFTAPlan_re{T,N}) where {T,N}
    return FFTAPlan_re{T,N}(p.callgraph, p.region, -p.dir, p.pinv, p.flen)
end

function LinearAlgebra.mul!(y::AbstractVector{U}, p::FFTAPlan{T,1}, x::AbstractVector{T}) where {T,U}
    fft!(y, x, 1, 1, p.dir, p.callgraph[1][1].type, p.callgraph[1], 1)
end

function LinearAlgebra.mul!(y::AbstractArray{U,N}, p::FFTAPlan{T,1}, x::AbstractArray{T,N}) where {T,U,N}
    Rpre = CartesianIndices(size(x)[1:p.region-1])
    Rpost = CartesianIndices(size(x)[p.region+1:end])
    for Ipre in Rpre
        for Ipost in Rpost
            @views fft!(y[Ipre,:,Ipost], x[Ipre,:,Ipost], 1, 1, p.dir, p.callgraph[1][1].type, p.callgraph[1], 1)
        end
    end
end

function LinearAlgebra.mul!(y::AbstractArray{U,N}, p::FFTAPlan{T,2}, x::AbstractArray{T,N}) where {T,U,N}
    R1 = CartesianIndices(size(x)[1:p.region[1]-1])
    R2 = CartesianIndices(size(x)[p.region[1]+1:p.region[2]-1])
    R3 = CartesianIndices(size(x)[p.region[2]+1:end])
    y_tmp = similar(y, axes(y)[p.region])
    rows,cols = size(x)[p.region]
    for I1 in R1
        for I2 in R2
            for I3 in R3
                for k in 1:cols
                    @views fft!(y_tmp[:,k],  x[I1,:,I2,k,I3], 1, 1, p.dir, p.callgraph[1][1].type, p.callgraph[1], 1)
                end

                for k in 1:rows
                    @views fft!(y[I1,k,I2,:,I3], y_tmp[k,:], 1, 1, p.dir, p.callgraph[2][1].type, p.callgraph[2], 1)
                end
            end
        end
    end
end

function *(p::FFTAPlan{T,1}, x::AbstractVector{T}) where {T<:Union{Real,Complex}}
    y = similar(x, T <: Real ? Complex{T} : T)
    LinearAlgebra.mul!(y, p, x)
    y
end

function *(p::FFTAPlan{T,N1}, x::AbstractArray{T,N2}) where {T<:Union{Real, Complex}, N1, N2}
    y = similar(x, T <: Real ? Complex{T} : T)
    LinearAlgebra.mul!(y, p, x)
    y
end

function *(p::FFTAPlan_re{T,1}, x::AbstractVector{T}) where {T<:Union{Real, Complex}}
    if p.dir == FFT_FORWARD
        y = similar(x, T <: Real ? Complex{T} : T)
        LinearAlgebra.mul!(y, p, x)
        return y[1:end÷2 + 1]
    else
        x_tmp = similar(x, p.flen)
        x_tmp[1:end÷2 + 1] .= x
        x_tmp[end÷2 + 2:end] .= iseven(p.flen) ? conj.(x[end-1:-1:2]) : conj.(x[end:-1:2])
        y = similar(x_tmp)
        LinearAlgebra.mul!(y, p, x_tmp)
        return y
    end
end

function *(p::FFTAPlan_re{T,N}, x::AbstractArray{T,2}) where {T<:Union{Real, Complex}, N}
    if p.dir == FFT_FORWARD
        y = similar(x, T <: Real ? Complex{T} : T)
        LinearAlgebra.mul!(y, p, x)
        return y[1:end÷2 + 1,:]
    else
        x_tmp = similar(x, p.flen, size(x)[2])
        x_tmp[1:end÷2 + 1,:] .= x
        x_tmp[end÷2 + 2:end,:] .= iseven(p.flen) ? conj.(x[end-1:-1:2,:]) : conj.(x[end:-1:2,:])
        y = similar(x_tmp)
        LinearAlgebra.mul!(y, p, x_tmp)
        return y
    end
end


text:
module Expr2LaTeX
using LaTeXStrings
using DocStringExtensions

export to_math, to_latex, render, renderval, @render

include("render.jl")
include("mathast.jl")
include("expr2math.jl")
include("math2latex.jl")

end # module


text:
# This file contains types and structures for economic DEA models
"""
    AbstractEconomicDEAlModel
An abstract type representing an economic DEA model.
"""
abstract type AbstractEconomicDEAModel <: AbstractDEAModel end

"""
    AbstractCostDEAModel
An abstract type representing a cost DEA model.
"""
abstract type AbstractCostDEAModel <: AbstractEconomicDEAModel end

"""
    AbstractRevenueDEAModel
An abstract type representing a revenue DEA model.
"""
abstract type AbstractRevenueDEAModel <: AbstractEconomicDEAModel end

"""
    AbstractProfitDEAModel
An abstract type representing a revenue DEA model.
"""
abstract type AbstractProfitDEAModel <: AbstractEconomicDEAModel end

"""
    AbstractProfitabilityDEAModel
An abstract type representing a revenue DEA model.
"""
abstract type AbstractProfitabilityDEAModel <: AbstractEconomicDEAModel end

"""
    efficiency(model::AbstractEconomicDEAModel)
Return efficiency scores of an economic DEA model.

# Optional Arguments
- `type=Economic`: type of efficiency scores to return.

Type specification:
- `:Economic`: returns economic efficiency of the model.
- `:Technical`: returns technical efficiency.
- `:Allocative`: returns allocative efficiency.

Some models also allow these types:
- `:CRS`: returns technical efficiency under constant returns to scale.
- `:VRS`: returns technical efficiency under variable returns to scale.
- `:Scale`: returns scale efficiency.

# Examples
```jldoctest
julia> X = [5 3; 2 4; 4 2; 4 8; 7 9.0];
julia> Y = [7 4; 10 8; 8 10; 5 4; 3 6.0];
julia> W = [2 1; 2 1; 2 1; 2 1; 2 1.0];
julia> P = [3 2; 3 2; 3 2; 3 2; 3 2.0];
julia> profitbl = deaprofitability(X, Y, W, P)
julia> efficiency(profitbl)
5-element Vector{Float64}:
 0.38795983677810825
 0.9999999082180037
 0.7652173234322985
 0.24999998462128403
 0.15879016408241559

julia> efficiency(profitbl, :Allocative)
5-element Vector{Float64}:
 0.6096511887087086
 0.9999999387832732
 0.765217344390068
 0.9999999256504392
 0.608695614904087
```
"""
function efficiency(model::AbstractEconomicDEAModel, type::Symbol = :Economic)::Vector

    if type == :Economic
        return model.eff
    end

    if type == :Technical
        if isdefined(model, :techeff)
            return model.techeff
        end
    end

    if type == :CRS
        if isdefined(model, :crseff)
            return model.crseff
        end
    end

    if type == :VRS
        if isdefined(model, :vrseff)
            return model.vrseff
        end
    end

    if type == :Scale
        if isdefined(model, :scaleff)
            return model.scaleff
        end
    end

    if type == :Allocative
        if isdefined(model, :alloceff)
            return model.alloceff
        end
    end

    throw(ArgumentError("$(typeof(model)) has no efficiency $(type)"));

end

"""
    targets(model::AbstractEconomicDEAModel, target::Symbol)
Return targets of an economic DEA model.
# Examples
```jldoctest
julia> X = [5 3; 2 4; 4 2; 4 8; 7 9.0];
julia> Y = [7 4; 10 8; 8 10; 5 4; 3 6.0];
julia> W = [2 1; 2 1; 2 1; 2 1; 2 1.0];
julia> P = [3 2; 3 2; 3 2; 3 2; 3 2.0];
julia> profit = deaprofit(X, Y, W, P, Gx = :Monetary, Gy = :Monetary);
julia> targets(profit, :X)
5×2 Matrix{Float64}:
 2.0  4.0
 2.0  4.0
 2.0  4.0
 2.0  4.0
 2.0  4.0

julia> targets(profit, :Y)
5×2 Matrix{Float64}:
 10.0  8.0
 10.0  8.0
 10.0  8.0
 10.0  8.0
 10.0  8.0

```
"""
function targets(model::AbstractEconomicDEAModel, target::Symbol)::Matrix

    if target == :X
        return model.Xtarget
    elseif target == :Y
        return model.Ytarget
    end

    throw(ArgumentError("`target` must be :X or :Y"));

end

"""
    normfactor(model::AbstractEconomicDEAModel)
Return the normalization factor of an economic DEA model.

# Examples
```jldoctest
julia> X = [5 3; 2 4; 4 2; 4 8; 7 9.0];
julia> Y = [7 4; 10 8; 8 10; 5 4; 3 6.0];
julia> W = [2 1; 2 1; 2 1; 2 1; 2 1.0];
julia> P = [3 2; 3 2; 3 2; 3 2; 3 2.0];
julia> profit = deaprofit(X, Y, W, P, Gx = :Ones, Gy = :Ones)
julia> normfactor(profit)
5-element Vector{Float64}:
 8.0
 8.0
 8.0
 8.0
 8.0
```
"""
function normfactor(model::AbstractEconomicDEAModel)::Vector
    if isdefined(model, :normalization)
        return model.normalization
    else
        throw(ArgumentError("$(typeof(model)) has no normalization factor"));
    end
end

"""
    ismonetary(model::AbstractEconomicDEAModel)
Indicate whether inefficiency is in monetary units.

# Examples
```jldoctest
julia> X = [5 3; 2 4; 4 2; 4 8; 7 9.0];
julia> Y = [7 4; 10 8; 8 10; 5 4; 3 6.0];
julia> W = [2 1; 2 1; 2 1; 2 1; 2 1.0];
julia> P = [3 2; 3 2; 3 2; 3 2; 3 2.0];
julia> profit = deaprofit(X, Y, W, P, Gx = :Ones, Gy = :Ones, monetary = true);
julia> ismonetary(profit)
true
```
"""
function ismonetary(model::AbstractEconomicDEAModel)::Bool
    if isdefined(model, :monetary)
        return model.monetary
    else
        throw(ArgumentError("$(typeof(model)) has no monetary identifier"));
    end
end


text:
# module Bukdu

module Logger

import ..Bukdu: ApplicationController, ApplicationLayout, ApplicationView, LayoutDivision

const level_false =   0
const level_fatal = 100
const level_error = 200
const level_warn  = 300
const level_info  = 400
const level_debug = 500

"""
    Logger.level::Union{Bool,Symbol}

Get the current Logger level.
"""
level = :info::Union{Bool,Symbol}

const levels = Dict{Union{Bool,Symbol},Int}(
    false  => level_false,
    :fatal => level_fatal,
    :error => level_error,
    :warn  => level_warn,
    :info  => level_info,
    :debug => level_debug
)

settings = Dict(
    :level => level_info,
    :have_color => Base.have_color,
    :have_datetime => false,
    :info_prefix => "INFO",
    :info_sub => ""
)

fatal(block::Function) = settings[:level] >= level_fatal && fatal(block())
error(block::Function) = settings[:level] >= level_error && error(block())
warn(block::Function)  = settings[:level] >= level_warn && warn(block())
info(block::Function)  = settings[:level] >= level_info && info(block())
debug(block::Function) = settings[:level] >= level_debug && debug(block())

fatal(args...; kw...) = settings[:level] >= level_fatal && print_log(:magenta, "FATAL", "", args...; kw...)
error(args...; kw...) = settings[:level] >= level_error && print_log(:red, "ERROR", "", args...; kw...)
warn(args...; kw...)  = settings[:level] >= level_warn && print_log(:yellow, "WARN ", "", args...; kw...)
info(args...; kw...)  = settings[:level] >= level_info && print_info(args...; kw...)
debug(args...; kw...) = settings[:level] >= level_debug && print_log(:cyan, "DEBUG", "", args...; kw...)

"""
    Logger.set_level(lvl::Union{Symbol,Bool})

Set the log level.
Options: `:debug`, `:info`, `:error`, `:warn`, `:fatal`, `false`

```julia
julia> Logger.set_level(false)
false
```
"""
function set_level(lvl::Union{Symbol,Bool})
    if haskey(levels, lvl)
        global level
        settings[:level] = levels[lvl]
        level = lvl
    else
        valids = join(map(repr, keys(levels)), ", ")
        throw(ArgumentError("invalid argument for level"))
    end
end

function print_info(args...; kw...)
    prefix = settings[:info_prefix]
    sub = string(settings[:info_sub])
    print_log(:green, prefix, sub, args...; kw...)
end

function inner_callstack(stackframe::Vector{StackFrame})
    !settings[:have_color] && return stackframe
    pat = r"(.* at )(?P<file>.*.jl):(?P<lineno>\d*)"
    string('\n', join(map(stackframe) do frame
        str = string(frame)
        m = match(pat, str)
        isa(m, RegexMatch) ? string(m[1], with_color(:bold, m[:file]), ':', with_color(:bold, m[:lineno])) : nothing
    end, '\n'))
end

function inner_contents(args...)
    if 1 == length(args)
        arg = first(args)
        if any(x->isa(arg,x), [Array, Tuple])
            if isa(arg, Vector{StackFrame})
                inner_callstack(arg)
            else
                inner_contents(arg...)
            end
        else
            arg
        end
    else
        join(map(inner_contents, args), ' ')
    end
end

function print_log(color::Symbol, prefix::String, sub::String, args...; LF=true)
    print(string(
        settings[:have_datetime] ? rpad(string(Dates.now()), 24) : "",
        with_color(color, prefix),
        ' ',
        isempty(sub) ? "" : string(sub, ' '),
        inner_contents(args...),
        (LF ? "\n" : "")))
end

function have_datetime(enabled::Bool)
    settings[:have_datetime] = enabled
end

function have_color(enabled::Bool)
    settings[:have_color] = enabled
end

function with_color(color::Symbol, text)::String
    if settings[:have_color]
        string(Base.text_colors[color], text, Base.color_normal)
    else
        string(text)
    end
end

## log_message
function log_message(prefix::String)
    settings[:info_prefix] = prefix
end

function log_message{AL<:ApplicationLayout}(D::LayoutDivision{AL})
    settings[:info_sub] = D
end

function log_message{AV<:ApplicationView}(::Type{AV})
    settings[:info_sub] = AV.name.name
end

function log_message(modul::Module)
    settings[:info_sub] = Base.module_name(modul)
end

function log_message{AC<:ApplicationController}(c::AC)
    controller = AC
    action = Base.function_name(c[:action])
    settings[:info_sub] = string(AC, '.', action)
end

end # module Bukdu.Logger


text:
export ft_reload

function ft_reload(vid::Integer, pid::Integer)
  ft_status = ccall((:FT_Reload, d2xx),
                     Cuint,
                     (Cuint, Cuint),
                     vid, pid)
  checkstatus(ft_status)
  return nothing
end