text:
module ITK
using Cxx, Libdl

include("init.jl")

function __init__()
    loadcxx(libitk, dirname(libitk))
end

"Loads shared library and header file."
function loadcxx(libpath::String=libitk, headerdir::String=dirname(libitk))
    addHeaderDir(headerdir, kind=C_System)
    Libdl.dlopen(libpath, Libdl.RTLD_GLOBAL)
    cxxinclude("JuliaWrap.h")
end

"Sanity checker. Returns x if able to read C++ files."
function verifycxx(x::Int)
    sanity_check(n::Int) = @cxx sanity(n)
    return sanity_check(x)
end

"""
    registerframe(fixedImage, movingImage, outputImage, writeImage, optimizer)

Test method which registers movingImage in relation to fixedImage, writing the resulting image to outputImage if writeImage is true.
Uses Mattes Mutual Information Metric, and either Regular Step Gradient Descent or Amoeba optimizer.

# Arguments
- `fixedImage::String`: path to fixed image file
- `movingImage::String`: path to moving image file
- `outputImage::String`: path to desired output image file
- `writeImage::Bool=true`: whether or not to save the output image file
- `optimizer::String="Gradient"`: which optimizer to use, "Gradient" | "Amoeba"

# Returns
- (x translation, y translation, metric information)
"""
function registerframe(fixedImage::String, movingImage::String, outputImage::String, writeImage::Bool=true, optimizer::String="Gradient")
    register_amoeba_optimizer(fix::Ptr{UInt8}, moving::Ptr{UInt8}, output::Ptr{UInt8}, write::Bool) = @cxx test_registration1(fix, moving, output, write)
    register_gradient_optimizer(fix::Ptr{UInt8}, moving::Ptr{UInt8}, output::Ptr{UInt8}, write::Bool) = @cxx test_registration2(fix, moving, output, write)
    if optimizer == "Amoeba"
        result = register_amoeba_optimizer(pointer(fixedImage), pointer(movingImage), pointer(outputImage), writeImage)
    else
        result = register_gradient_optimizer(pointer(fixedImage), pointer(movingImage), pointer(outputImage), writeImage)
    end
    x, y, metric = unsafe_load(result,1), unsafe_load(result,2), unsafe_load(result,3)
    return x, y, metric
end

"""
    MMIGradientRegistration(fixedImage, movingImage, outputImage, writeImage, learningRate, minStepLength, maxIterations, relaxationFactor)

Register movingImage to fixedImage with Mattes Mutual Information metric and Regular Step Gradient Descent optimizer.

# Arguments
- `fixedImage::String`: path to fixed image file
- `movingImage::String`: path to moving image file
- `outputImage::String`: path to desired output image file
- `writeImage::Bool`: whether or not to save the output image file
- `learningRate::Float64`: higher values result in quicker but potentially unstable registration (recommended range 1.0:5.0)
- `minStepLength::Float64`: step size (recommended range 0.001:0.05)
- `maxIterations::Int64`: max number of iterations (recommended range 100:200)
- `relaxationFactor::Float64`: factor to slow down shift between iterations, higher values result in quicker, potentially unstable registration (recommended range 0.3:0.7)

# Returns
- (x translation, y translation, metric information)
"""
function MMIGradientRegistration(fixedImage::String, movingImage::String, outputImage::String, writeImage::Bool, learningRate::Float64, minStepLength::Float64, maxIterations::Int64, relaxationFactor::Float64)
    register(fix::Ptr{UInt8}, moving::Ptr{UInt8}, output::Ptr{UInt8}, write::Bool, LR::Float64, minStep::Float64, maxIter::Int64, relaxFactor::Float64) = @cxx MMIGradientTranslation(fix, moving, output, write, LR, minStep, maxIter, relaxFactor)
    result = register(pointer(fixedImage), pointer(movingImage), pointer(outputImage), writeImage, learningRate, minStepLength, maxIterations, relaxationFactor)
    x, y, metric = unsafe_load(result, 1), unsafe_load(result, 2), unsafe_load(result, 3)
    return x, y, metric
end

"""
    MMIAmoebaRegistration(fixedImage, movingImage, outputImage, writeImage, initialSimplex, pixelTolerance, metricTolerance, maxIterations)

Registers movingImage to fixedImage with MattesMutualInformation metric, with Amoeba optimizer with specified parameters.

# Arguments
- `fixedImage::String`: path to fixed image file
- `movingImage::String`: path to moving image file
- `outputImage::String`: path to desired output image file
- `writeImage::Bool`: whether or not to save the output image file
- `initialSimplex::Float64`: initial size of simplex moving along cost surface (recommended 5.0)
- `pixelTolerance::Float64`: convergence tolerance in pixels (recommended range 0.01:0.5)
- `metricTolerance::Float64`: function convergence tolerance (recommended 0.1)
- `maxIterations::Int64`: max number of iterations (recommended range 100:200)

# Returns
- (x translation, y translation, metric information)
"""
function MMIAmoebaRegistration(fixedImage::String, movingImage::String, outputImage::String, writeImage::Bool, initialSimplex::Float64, pixelTolerance::Float64, metricTolerance::Float64, maxIterations::Int64)
    register(fix::Ptr{UInt8}, moving::Ptr{UInt8}, output::Ptr{UInt8}, write::Bool, simplex::Float64, pixelT::Float64, metricT::Float64, maxIter::Int64) = @cxx MMIAmoebaTranslation(fix, moving, output, write, simplex, pixelT, metricT, maxIter)
    result = register(pointer(fixedImage), pointer(movingImage), pointer(outputImage), writeImage, initialSimplex, pixelTolerance, metricTolerance, maxIterations)
    x, y, metric = unsafe_load(result, 1), unsafe_load(result, 2), unsafe_load(result, 3)
    return x, y, metric
end

"""
    MMIOnePlusOneRegistration(fixedImage, movingImage, outputImage, writeImage, initialize, epsilon, maxIterations)

Registers movingImage to fixedImage with MattesMutualInformation metric, with 1+1 Evolutionary optimizer with specified parameters.

# Arguments
- `fixedImage::String`: path to fixed image file
- `movingImage::String`: path to moving image file
- `outputImage::String`: path to desired output image file
- `writeImage::Bool`: whether or not to save the output image file
- `initialRadius::Float64`: initial size of radius (recommended 6.25e-3)
- `epsilon::Float64`: minimum size of search radius (recommended 1.5e-6)
- `maxIterations::Int64`: max number of iterations (recommended range 100:200)

# Returns
- (x translation, y translation, metric information)
"""
function MMIOnePlusOneRegistration(fixedImage::String, movingImage::String, outputImage::String, writeImage::Bool, initialRadius::Float64, epsilon::Float64, maxIterations::Int64)
    register(fix::Ptr{UInt8}, moving::Ptr{UInt8}, output::Ptr{UInt8}, write::Bool, init::Float64, ep::Float64, maxIter::Int64) = @cxx MMIOnePlusOneTranslation(fix, moving, output, write, init, ep, maxIter)
    result = register(pointer(fixedImage), pointer(movingImage), pointer(outputImage), writeImage, initialRadius, epsilon, maxIterations)
    x, y, metric = unsafe_load(result, 1), unsafe_load(result, 2), unsafe_load(result, 3)
    return x, y, metric
end

"""
MeanSquaresGradientRegistration(fixedImage, movingImage, outputImage, writeImage, learningRate, minStepLength, maxIterations, relaxationFactor)

Register movingImage to fixedImage with Mean Squares metric and Regular Step Gradient Descent optimizer.

# Arguments
- `fixedImage::String`: path to fixed image file
- `movingImage::String`: path to moving image file
- `outputImage::String`: path to desired output image file
- `writeImage::Bool`: whether or not to save the output image file
- `learningRate::Float64`: higher values result in quicker but potentially unstable registration (recommended range 1.0:5.0)
- `minStepLength::Float64`: step size (recommended range 0.001:0.05)
- `maxIterations::Int64`: max number of iterations (recommended range 100:200)
- `relaxationFactor::Float64`: factor to slow down shift between iterations, higher values result in quicker, potentially unstable registration (recommended range 0.3:0.7)

# Returns
- (x translation, y translation, metric information)
"""
function MeanSquaresGradientRegistration(fixedImage::String, movingImage::String, outputImage::String, writeImage::Bool, learningRate::Float64, minStepLength::Float64, maxIterations::Int64, relaxationFactor::Float64)
    register(fix::Ptr{UInt8}, moving::Ptr{UInt8}, output::Ptr{UInt8}, write::Bool, LR::Float64, minStep::Float64, maxIter::Int64, relaxFactor::Float64) = @cxx MeanSquaresGradientTranslation(fix, moving, output, write, LR, minStep, maxIter, relaxFactor)
    result = register(pointer(fixedImage), pointer(movingImage), pointer(outputImage), writeImage, learningRate, minStepLength, maxIterations, relaxationFactor)
    x, y, metric = unsafe_load(result, 1), unsafe_load(result, 2), unsafe_load(result, 3)
    return x, y, metric
end

"""
NormalizedMutualInformationHistogramImageToImageMetric(fixedImage, movingImage, outputImage, writeImage, learningRate, minStepLength, maxIterations, relaxationFactor)

Register movingImage to fixedImage with mutual histogram and Regular Step Gradient Descent optimizer.

# Arguments
- `fixedImage::String`: path to fixed image file
- `movingImage::String`: path to moving image file
- `outputImage::String`: path to desired output image file
- `writeImage::Bool`: whether or not to save the output image file
- `learningRate::Float64`: higher values result in quicker but potentially unstable registration (recommended range 1.0:5.0)
- `minStepLength::Float64`: step size (recommended range 0.001:0.05)
- `maxIterations::Int64`: max number of iterations (recommended range 100:200)
- `relaxationFactor::Float64`: factor to slow down shift between iterations, higher values result in quicker, potentially unstable registration (recommended range 0.3:0.7)

# Returns
- (x translation, y translation, metric information)
"""
function NormalizedMIHImageToImageMetric(fixedImage::String, movingImage::String, outputImage::String, writeImage::Bool, learningRate::Float64, minStepLength::Float64, maxIterations::Int64, relaxationFactor::Float64)
    register(fix::Ptr{UInt8}, moving::Ptr{UInt8}, output::Ptr{UInt8}, write::Bool, LR::Float64, minStep::Float64, maxIter::Int64, relaxFactor::Float64) = @cxx NormalizedMutualInformationHistogramImageToImageMetric(fix, moving, output, write, LR, minStep, maxIter, relaxFactor)
    result = register(pointer(fixedImage), pointer(movingImage), pointer(outputImage), writeImage, learningRate, minStepLength, maxIterations, relaxationFactor)
    x, y, metric = unsafe_load(result, 1), unsafe_load(result, 2), unsafe_load(result, 3)
    return x, y, metric
end

export registerframe
export MMIGradientRegistration
export MMIAmoebaRegistration
export MMIOnePlusOneRegistration
export MeanSquaresGradientRegistration
export NormalizedMIHImageToImageMetric

end # module ITK


text:
L = 10
labels = rand(Bool, L);
scores = add_noise.(labels, 0.6)
curve = ROC(labels, scores)
plot(curve)

curve = PR(labels, scores)
plot(curve)

curve = DET(labels, scores)
plot(curve)

ext: jl


text:
module ContinuedFractions

import Base: start, done, next, length, eltype, collect

export ContinuedFraction, quotients, convergents, ConvergentIterator

immutable ContinuedFraction{T<:Integer}
	quotients::Vector{T}
end

quotients(cf::ContinuedFraction) = cf.quotients

immutable ConvergentIterator{T<:Integer}
    qs::Vector{T}
end

start(::ConvergentIterator) = 1
done(it::ConvergentIterator, state::Int) = state > length(it.qs)
length(it::ConvergentIterator) = length(it.qs)

function next(it::ConvergentIterator, state::Int)
    convergent = Rational(ContinuedFraction(it.qs[1:state]))
    convergent, state + 1
end

eltype(it::ConvergentIterator) = Rational{eltype(it.qs)}
collect(it::ConvergentIterator) = collect(eltype(it), it)

convergents(cf::ContinuedFraction) = convergents(quotients(cf))
convergents{T<:Integer}(qs::Vector{T}) = ConvergentIterator(qs)

function Base.Rational(cf::ContinuedFraction)
    qs = quotients(cf)
    isempty(qs) && return 0 // 1
    length(qs) == 1 && return qs[1] // 1

    remainder = qs[2:end]
    rat = Rational(ContinuedFraction(remainder))
    (qs[1] * rat.num + rat.den) // rat.num
end

function ContinuedFraction{T<:Integer}(rat::Rational{T})
    a = div(rat.num, rat.den)
    a * rat.den == rat.num && return ContinuedFraction(T[a])  # Exact!

    cf = ContinuedFraction(rat.den//(rat.num - a*rat.den))
    unshift!(quotients(cf), a) # insert at index 1
    cf
end

# let rationalize handle conversion from floating point
ContinuedFraction(x::BigFloat) = ContinuedFraction(rationalize(BigInt, x))
ContinuedFraction(x::AbstractFloat) = ContinuedFraction(rationalize(x))

ContinuedFraction(x::Integer) = ContinuedFraction([x])

end # module ContinuedFractions


text:
"""
    product(f, a, b)

Takes the Cartesian outer product of two containers and evaluates `f` on all pairings of
elements.

For example, if `a` and `b` are vectors, this returns a matrix `out` such that
`out[i,j] = f(a[i], b[j])` for `i in keys(a)` and `j in keys(b)`. See also `productview`.

# Example

```julia
julia> product(+, [1,2], [1,2,3])
2×3 Array{Int64,2}:
 2  3  4
 3  4  5
```
"""
function product(f::Callable, a, b)
	out = similar(a, promote_op(f, eltype(a), eltype(b)), (axes(a)..., axes(b)...))
	return product!(f, out, a, b)
end

@inline function product!(f::Callable, out, a, b)
	axes_a = axes(a)
	axes_b = axes(b)
	axes_out = axes(out)

	@boundscheck if axes_out != (axes_a..., axes_b...)
	    throw(DimensionMismatch("Output does not match input dimensions"))
	end

	@inbounds for i_a in CartesianIndices(axes_a)
		for i_b in CartesianIndices(axes_b)
			i_out = CartesianIndex(i_a.I..., i_b.I...)
			out[i_out] = f(a[i_a], b[i_b])
		end
	end

	return out
end

struct ProductArray{T, N, F, A <: AbstractArray, B <: AbstractArray} <: AbstractArray{T, N}
	f::F
	a::A
	b::B
end

Base.IndexStyle(::Type{<:ProductArray}) = IndexCartesian()
Base.axes(p::ProductArray) = (axes(p.a)..., axes(p.b)...)
Base.size(p::ProductArray) = (size(p.a)..., size(p.b)...)
@propagate_inbounds function Base.getindex(p::ProductArray{T}, i::Int...) where {T}
	ndims_a = ndims(p.a)
	ndims_b = ndims(p.b)
	i_a = ntuple(j -> @inbounds(return i[j]), Val(ndims_a))
	i_b = ntuple(j -> @inbounds(return i[j + ndims_a]), Val(ndims_b))
	return p.f(p.a[i_a...], p.b[i_b...])::T
end

"""
    productview(f, a, b)

Return a view of a Cartesian product of `a` and `b` where the output elements are `f`
evaluated those of `a` and `b`. See also `product` and `ProductArray`

# Example

```julia
julia> productview(+, [1,2], [1,2,3])
2×3 ProductArray{Int64,2,typeof(+),Array{Int64,1},Array{Int64,1}}:
 2  3  4
 3  4  5
```
"""
function productview(f::Callable, a, b)
	T = promote_op(f, eltype(a), eltype(b))
	N = ndims(a) + ndims(b)
	return ProductArray{T, N, typeof(f), typeof(a), typeof(b)}(f, a, b)
end


text:
#
# V.0.1 (2016.11.04)
#
# ROC curve for the real data.
# Run experiment with the real data set.
# Plot ROC curve.
# Corresponding tests are in: test4_concept_proof_real.R, .jl
#
include("./lib/metadata.jl")
include("./lib/pccf.jl")
include("./2bayes_detector_replica.jl")
include("./lib/SignalGenerator.jl")
include("./lib/roc.jl")

using SignalGenerator
using BayesChpOrig
using Metadata
using Roc
using Pccf

"""
In this experiment we measure perfoemance of the Online Bayesian Detector
1) without PCCF
2) with PCCF
and plot correspondint ROC curves (exp1.R).
"""
function run_exp1(inputDataFolder, outputfolder)
    #work_path = "./data/activity-recognition-sample/"
    #signal = readdlm(work_path * "sample-sig-10.dat", ',')
    println("... exp1.jl.run_exp > read input data from ", inputDataFolder)
    println("... exp1.jl.run_exp > save results into ", outputfolder)

    signal      = readdlm(inputDataFolder * "sample-sig-10.dat", ',')
    changes     = readdlm(inputDataFolder * "changes-sig-active-10.dat")
    changes     = reshape(changes, length(changes))
    mid_changes = readdlm(inputDataFolder * "mid-changes-sig-active-10.dat")
    #signal      = readdlm(projectdir * "src_julia/data/activity-recognition-sample/sample-sig-10.dat", ',')
    #changes     = readdlm(projectdir * "src_julia/data/activity-recognition-sample/changes-sig-active-10.dat")
    #changes     = reshape(changes, length(changes))
    #mid_changes = readdlm(projectdir * "src_julia/data/activity-recognition-sample/mid-changes-sig-active-10.dat")

    # BEFORE ADDING MID_CHANGES
    mu0 = mean(diff(changes)) # ~84.6
    sd0 = std(diff(changes))  # ~9.16
    @printf("... exp1.jl > mu = %0.2f sigma = %0.2f\n", mu0, sd0)
    # end BEFORE ADDING MID_CHANGES

    if 1 == 1
        # Add middle changes
        mid_changes = reshape(mid_changes, length(mid_changes))
        changes = vcat(changes, mid_changes)
        sort!(changes)
    end

    x = signal[:,1]
    y = signal[:,2]
    n = length(x)

    max_delay = 10

    h_from    = 80   # target values (~ test4_concept_proof_real.jl) l_common = 100
    h_to      = 300  # 300 - on ROC points become too close
    h_by      = 10
    threshold = 0.3

    tp_rates = []
    fa_rates = []
    tp_rates_pccf = []
    fa_rates_pccf = []
    h_vals = []

    probs = scale01(pccf(n, mu0, sd0)) # usefull only at the beginning

    # Make input signal shorter so that PCCF effect is more visible (also below in a loop)
    shortenInput = false # true
    if shortenInput
        max_t = 500
        x = x[1:max_t]
        changes = filter(x -> x <= max_t, changes)
        probs = probs[1:max_t]
    end
    # end Make input signal shorter so that PCCF effect is more visible (also below in a loop)

    for h in h_from : h_by : h_to
        push!(h_vals, h)

        l_common = h
        if mod(h, 10)==0
            # print(h/h_to,"\n")
        end

        maxes, mat = bayes_detector(x, lambda = l_common)
        detections = extract_changes(maxes)

        pccf_maxes, pccf_mat = pccf_bayes_detector(x, probs, l_common + 100, l_common, threshold, h_common = l_common, max_d = 500)
        pccf_detections = extract_changes(pccf_maxes)

        # Make input signal shorter so that PCCF effect is more visible
        if shortenInput
            detections      = filter(x-> x<= max_t, detections)
            pccf_detections = filter(x-> x<= max_t, pccf_detections)
        end
        # END Make input signal shorter so that PCCF effect is more visible

        tpr      = tpRate(detections, changes, max_delay)
        tpr_pccf = tpRate(pccf_detections, changes, max_delay)

        fpr      = faRateScaled(detections     , changes, max_delay, n, tn_price_scale = 0.1 * 0.5)
        fpr_pccf = faRateScaled(pccf_detections, changes, max_delay, n, tn_price_scale = 0.1 * 0.5)

        # To find optimal h values to demonstrate PCCF effect
        if fpr <= 0.30 && tpr >= 0.6
            @printf("... exp1.jl > TP      = %0.3f / FP      = %0.3f /  H = %0.3f \n"  , tpr, fpr, h)
            @printf("... exp1.jl > TP_PCCF = %0.3f / FP_PCCF = %0.3f /  H = %0.3f \n\n", tpr_pccf, fpr_pccf, h)
        end

        push!(tp_rates, tpr)
        push!(fa_rates, fpr)
        push!(tp_rates_pccf, tpr_pccf)
        push!(fa_rates_pccf, fpr_pccf)

    end

    writedlm(outputfolder * "tp_rates.dat"     , tp_rates)
    writedlm(outputfolder * "fa_rates.dat"     , fa_rates)
    writedlm(outputfolder * "tp_rates_pccf.dat", tp_rates_pccf)
    writedlm(outputfolder * "fa_rates_pccf.dat", fa_rates_pccf)
    writedlm(outputfolder * "h_vals.dat"       , h_vals)
    # writedlm(outdir * "exp1/tp_rates.dat", tp_rates)
    # writedlm(outdir * "exp1/fa_rates.dat", fa_rates)
    # writedlm(outdir * "exp1/tp_rates_pccf.dat", tp_rates_pccf)
    # writedlm(outdir * "exp1/fa_rates_pccf.dat", fa_rates_pccf)
    # writedlm(outdir * "exp1/h_vals.dat", h_vals)

    println("... exp1.jl > Results are saved into ", outdir * "exp1/")

end

run_exp1(projectdir * "src_julia/data/activity-recognition-sample/",
         outdir * "exp1/")

# changes =
# max_delay = 10
# tpr = tpRate(detections, chps, max_delay)
# fpr = faRateScaled(detections, changes, max_delay, length(sig), tn_price_scale=0.1)
# println("... test_bd > TPrate = ", tpr, " / FArate = ", fpr)

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
include("trace.jl")

struct GenFunction <: Generator{Any,GFTrace}
    params_grad::Dict{Symbol,Any}
    params::Dict{Symbol,Any}
    arg_types::Vector{Type}
    julia_function::Function
    has_argument_grads::Vector{Bool}
    accepts_output_grad::Bool
end

function GenFunction(arg_types::Vector{Type}, julia_function::Function,
                     has_argument_grads, accepts_output_grad::Bool)
    params_grad = Dict{Symbol,Any}()
    params = Dict{Symbol,Any}()
    GenFunction(params_grad, params, arg_types, julia_function,
                has_argument_grads, accepts_output_grad)
end

function (g::GenFunction)(args...)
    trace = simulate(g, args)
    call = get_call_record(trace)
    call.retval
end

exec(gf::GenFunction, state, args::Tuple) = gf.julia_function(state, args...)

get_static_argument_types(gen::GenFunction) = gen.arg_types

# whether there is a gradient of score with respect to each argument
# it returns 'nothing' for those arguemnts that don't have a derivatice
has_argument_grads(gen::GenFunction) = gen.has_argument_grads

# if it is true, then it expects a value, otherwise it expects 'nothing'
accepts_output_grad(gen::GenFunction) = gen.accepts_output_grad

const self = gensym("self")

function parse_arg_types(args)
    types = Vector{Any}()
    for arg in args
        if isa(arg, Expr) && arg.head == :(::)
            typ = arg.args[2]
            push!(types, typ)
        elseif isa(arg, Symbol)
            push!(types, :Any)
        else
            error("Error parsing argument to generative function: $arg")
        end
    end
    types
end

# TODO handle @compiled @ad @gen function and @ad compiled @gen function
macro ad(ast)
    if ast.head != :macrocall || ast.args[1] != Symbol("@ad") || length(ast.args) != 2
        error("Syntax error in @ad $ast")
    end
    gen_ast = ast.args[2]
    if (!isa(gen_ast, Expr)
        || gen_ast.head != :macrocall
        || gen_ast.args[1] != Symbol("@gen"))
        error("Syntax error in @ad $ast")
    end
    parse_gen_function(gen_ast, true)
end

macro addr(expr::Expr, addr, args_change)
    if expr.head != :call
        error("syntax error in @addr at $(expr)")
    end
    fn = esc(expr.args[1])
    args = map(esc, expr.args[2:end])
    Expr(:call, :addr, esc(state), fn, Expr(:tuple, args...), esc(addr), esc(args_change))
end

macro gen(ast)
    parse_gen_function(ast, false)
end

function args_for_gf_julia_fn(args, has_argument_grads)
    # if an argument is marked for AD, we change remove the type annotation for
    # the Julia function implementation, because we use reverse mode AD using
    # boxed tracked values. we don't give any type information about the argument.
    escaped_args = []
    for (arg, has_grad) in zip(args, has_argument_grads)
        if has_grad
            if isa(arg, Expr) && arg.head == :(::)
                push!(escaped_args, esc(arg.args[1]))
            elseif isa(arg, Symbol)
                push!(escaped_args, esc(arg))
            end
        else
            push!(escaped_args, esc(arg))
        end
    end
    escaped_args
end

function parse_gen_function(ast, ad_annotation::Bool)
    if ast.head != :function
        error("syntax error at $(ast) in $(ast.head)")
    end
    if length(ast.args) != 2
        error("syntax error at $(ast) in $(ast.args)")
    end
    signature = ast.args[1]
    body = ast.args[2]
    if signature.head != :call
        error("syntax error at $(ast) in $(signature)")
    end
    args = signature.args[2:end]
    has_argument_grads = map(marked_for_ad, args)
    args = map(strip_marked_for_ad, args)

    # julia function definition
    escaped_args = args_for_gf_julia_fn(args, has_argument_grads)
    gf_args = [esc(state), escaped_args...]
    julia_fn_defn = Expr(:function,
        Expr(:tuple, gf_args...),
        esc(body))

    # create generator and assign it a name
    generator_name = signature.args[1]
    arg_types = map(esc, parse_arg_types(args))
    Expr(:block,
        Expr(:(=),
            esc(generator_name),
            Expr(:call, :GenFunction,
                quote Type[$(arg_types...)] end,
                julia_fn_defn, has_argument_grads, ad_annotation)))
end

function address_not_found_error_msg(addr)
    "Address $addr not found"
end

function read_param(state, name::Symbol)
    state.params[name]
end

#####################
# static parameters #
#####################

macro param(expr_or_symbol)
    local name::Symbol
    if isa(expr_or_symbol, Symbol)
        name = expr_or_symbol
    elseif expr_or_symbol.head == :(::)
        name = expr_or_symbol.args[1]
    else
        error("Syntax in error in @static at $(expr_or_symbol)")
    end
    Expr(:(=), esc(name), Expr(:call, :read_param, esc(state), QuoteNode(name)))
end

function set_param!(gf::GenFunction, name::Symbol, value)
    gf.params[name] = value
end

function get_param(gf::GenFunction, name::Symbol)
    gf.params[name]
end

function get_param_grad(gf::GenFunction, name::Symbol)
    gf.params_grad[name]
end

function zero_param_grad!(gf::GenFunction, name::Symbol)
    gf.params_grad[name] = zero(gf.params[name])
end

function init_param!(gf::GenFunction, name::Symbol, value)
    set_param!(gf, name, value)
    zero_param_grad!(gf, name)
end

#########
# delta #
#########

macro argschange()
    Expr(:call, :get_args_change, esc(state))
end

macro change(addr)
    Expr(:call, :get_addr_change, esc(state), esc(addr))
end

macro retchange!(value)
    Expr(:call, :set_ret_change!, esc(state), esc(value))
end

function addr(state, gen::Generator{T,U}, args, address) where {T,U}
    addr(state, gen, args, address, nothing)
end

function addr(state, dist::Distribution{T}, args, address) where {T}
    addr(state, dist, args, address, nothing)
end

##################
# AddressVisitor #
##################

struct AddressVisitor
    visited::DynamicAddressSet
end

AddressVisitor() = AddressVisitor(DynamicAddressSet())

function visit!(visitor::AddressVisitor, addr)
    push_leaf_node!(visitor.visited, addr)
end

function _diff(trie, addrs::AddressSet)
    diff_trie = HomogenousTrie{Any,Any}()

    for (key_field, value) in get_leaf_nodes(trie)
        if !has_leaf_node(addrs, key_field)
            set_leaf_node!(diff_trie, key_field, value)
        end
    end

    for (key_field, node) in get_internal_nodes(trie)
        if !has_leaf_node(addrs, key_field)
            if has_internal_node(addrs, key_field)
                sub_addrs = get_internal_node(addrs, key_field)
            else
                sub_addrs = EmptyAddressSet()
            end
            diff_node = _diff(node, sub_addrs)
            if !isempty(diff_node)
                set_internal_node!(diff_trie, key_field, diff_node)
            end
        end
    end

    diff_trie
end

function get_unvisited(visitor::AddressVisitor, choices)
    _diff(choices, visitor.visited)
end

function lightweight_got_internal_node_err(addr)
    error("Expected a value at address $addr but trace contained an assignment")
end

function lightweight_got_leaf_node_err(addr)
    error("Expected an assignment at address $addr but trace contained a value")
end

function lightweight_retchange_already_set_err()
    error("@retchange! was already used")
end

include("simulate.jl")
include("assess.jl")
include("generate.jl")
include("fix_update.jl")
include("update.jl")
include("extend.jl")
include("project.jl")
include("regenerate.jl")
include("ungenerate.jl")
include("backprop_params.jl")
include("backprop_trace.jl")

export set_param!, get_param, get_param_grad, zero_param_grad!, init_param!
export @ad
export @delta
export @param
export @addr
export @gen
export @change
export @argschange
export @retchange!

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
using Documenter, MonteCarloX

DocMeta.setdocmeta!(MonteCarloX, :DocTestSetup, :(using MonteCarloX); recursive = true)

makedocs(;
    modules = [MonteCarloX],
    format = Documenter.HTML(),
    pages = [
        "Getting Started" => "index.md",
        "Equilibrium Tools" => "equilibrium.md",
        "Non-Equilibrium Tools" => "non_equilibrium.md",
        "Helper" => "helper.md",
    ],
    sitename = "MonteCarloX",
    authors = "Johannes Zierenberg",)

deploydocs(;
    deps = Deps.pip("pygments", "mkdocs", "python-markdown-math"),
    repo = "github.com/zierenberg/MonteCarloX.jl",
    target="build"
)


text:
using ModelingToolkit
using Test

struct Lorenz <: AbstractComponent
    x::Variable
    y::Variable
    z::Variable
    σ::Variable
    ρ::Variable
    β::Variable
    eqs::Vector{Equation}
end
function generate_lorenz_eqs(t,x,y,z,σ,ρ,β)
    D = Differential(t)
    [D(x) ~ σ*(y-x)
     D(y) ~ x*(ρ-z)-y
     D(z) ~ x*y - β*z]
end
function Lorenz(t)
    @Unknown x(t) y(t) z(t)
    @Param σ ρ β
    Lorenz(x, y, z, σ, ρ, β, generate_lorenz_eqs(t, x, y, z, σ, ρ, β))
end

@Param t
lz1 = Lorenz(t)
lz2 = Lorenz(t)
lz1.x ~ lz2.x


text:
# A chain of layers used to build models
struct Chain
    layers::Array
    Chain(layers...) = new(collect(layers))
end
(c::Chain)(x) = (for l in c.layers; x = l(x); end; x)
(c::Chain)(x,y) = nll(c(x),y)
(c::Chain)(d::Data) = mean(c(x,y) for (x,y) in d)

# Linear layer
mutable struct Linear; w; b; end
Linear(i::Int,o::Int) = Linear(param(o,i),param0(o))
(l::Linear)(x) = l.w * x .+ l.b

# Softmax classifier layer
mutable struct SoftmaxCls; w; b; end
SoftmaxCls(i::Int,o::Int) = SoftmaxCls(param(o,i),param0(o))
(l::SoftmaxCls)(x) = softmax(l.w * x .+ l.b)

# Dense (normal) layer
mutable struct Dense; w; b; f; pdrop; end
Dense(i::Int,o::Int,f=relu; pdrop=0) = Dense(param(o,i),param0(o),f,pdrop)
(l::Dense)(x) = l.f.(l.w * dropout(x,l.pdrop) .+ l.b)

# Flatten layer. I had to make it a struct because otherwise it does not save to file
struct Flatten; n::Int end
(r::Flatten)(x) = reshape(x, (r.n,:))

# Convolutional + pooling layer
mutable struct ConvPool; w; b; f; p; end
(c::ConvPool)(x) = c.f.(pool(conv4(c.w, dropout(x,c.p)) .+ c.b))
ConvPool(w1::Int,w2::Int,cx::Int,cy::Int,f=relu;pdrop=0) = ConvPool(param(w1,w2,cx,cy), param0(1,1,cy,1), f, pdrop)

# Convolutional layer
mutable struct Conv; w; b; f; padding; stride; end
function (c::Conv)(x)
    c.f.(conv4(c.w, x, padding=c.padding, stride=c.stride) .+ c.b)
end
function Conv(w1::Int,w2::Int,cx::Int,cy::Int,f=relu;padding=0,stride=1)
    return Conv(param(w1,w2,cx,cy), param0(1,1,cy,1), f, padding, stride)
end

function my_bn(x, bn_moments, bn_params)
    eps = 1e-5
    ivar = 1 ./ sqrt.(bn_moments.var .+ eps)

    g = reshape(bn_params[1:size(x, 3)], (1, 1, size(x, 3), 1))
    b = reshape(bn_params[size(x, 3)+1:end], (1, 1, size(x, 3), 1))
    return g .* (x .- bn_moments.mean) .* ivar .+ b
end

# Convolutional + Batchnorm layer
mutable struct ConvBN; w; b; f; padding; stride; bn_params; bn_moments; end
function (c::ConvBN)(x)
    y_conv = conv4(c.w, x, padding=c.padding, stride=c.stride) .+ c.b
    c.f.(batchnorm(y_conv, c.bn_moments, c.bn_params))
end
function (c::ConvBN)(x, skip_bn::Bool)
    y_conv = conv4(c.w, x, padding=c.padding, stride=c.stride) .+ c.b
    if skip_bn
        return c.f.(y_conv)
    end
    return c.f.(batchnorm(y_conv, c.bn_moments, c.bn_params))
end
function ConvBN(w1::Int,w2::Int,cx::Int,cy::Int,f=relu;padding=0,stride=1)
    return ConvBN(param(w1,w2,cx,cy), param0(1,1,cy,1), f, padding, stride, Param(convert(atype(), bnparams(cy))), bnmoments())
end

# Pooling layer
struct Pool
    window
    stride
    padding
    mode
    Pool(window=2, stride=2, padding=0, mode=0) = new(window, stride, padding, mode)
end
(p::Pool)(x) = pool(x, window=p.window, stride=p.stride, padding=p.padding, mode=p.mode)

include("inception.jl")

"Builds a CNN model consisting of 3x3 filters and 2x2 max pools"
function create_cnn_model(num_channels::Int, num_classes::Int, use_bn::Bool=false)
    C = use_bn ? ConvBN : Conv
    Chain(
        C(3, 3, num_channels, 32),
        C(3, 3, 32, 48),
        Pool(), # 14x14

        C(3, 3, 48, 64),
        C(3, 3, 64, 96),
        Pool(), # 5x5

        C(3, 3, 96, 128),
        C(3, 3, 128, 192),

        Flatten(192),
        Linear(192, num_classes)
    )
end

"Builds an MLP model with any number of hidden layers of any given unit sizes"
function create_mlp_model(i::Int, o::Int, h_units::Int ... ; f=relu)
    weight_dims = [i]
    for u in h_units
        push!(weight_dims, u, u)
    end
    push!(weight_dims, o)

    # Unit counts are used as such: [(i, h1), (h1, h2), (h2, h3), ..., (hn, o)]
    layer_params = []
    for i::Int in 1:(length(weight_dims)/2)
        d1 = weight_dims[2*i - 1]
        d2 = weight_dims[2*i]
        push!(layer_params, [d1, d2, f])
    end

    # Set the last function to be identity
    layer_params[end][3] = identity

    @show layer_params

    layers = (Dense(lp...) for lp in layer_params)
    return Chain(layers...)
end

# Utility function that calculates the output channel size of Inception modules
function out_channel_count(i::InceptionA)
    return size(i.c1_alone.w, 4) + size(i.c3.w, 4) + size(i.cd3_2.w, 4) + size(i.c1_after_pool.w, 4)
end
function out_channel_count(i::InceptionB)
    return size(i.c3.w, 4) + size(i.cd3_2.w, 4) + size(i.c1_before_3.w, 3)
end


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

using Test, Tensors

let dtime = 0.25,
    parameters = ChabocheParameterState(E=200.0e3,
                                        nu=0.3,
                                        R0=100.0,
                                        Kn=100.0,
                                        nn=10.0,
                                        C1=10000.0,
                                        D1=100.0,
                                        C2=50000.0,
                                        D2=1000.0,
                                        Q=50.0,
                                        b=0.1),
    mat = Chaboche(parameters = parameters),
    dstrain11 = 1e-3*dtime,
    dstrain12 = 1e-3*dtime,
    dtimes = dtime*[1.0, 1.0, 1.0, 1.0, 4.0],
    dstrains11 = dstrain11*[1.0, 1.0, 1.0, -1.0, -4.0],
    dstrains12 = dstrain12*[1.0, 1.0, 1.0, -1.0, -4.0]

    plastic_flow_occurred = zeros(Bool, length(dtimes) - 1)
    for i in 1:length(dtimes)
        dstrain11 = dstrains11[i]
        dstrain12 = dstrains12[i]
        dtime = dtimes[i]
        biaxial_increment!(mat, dstrain11, dstrain12, dtime)
        update_material!(mat)
        if i > 1
            plastic_flow_occurred[i-1] = (mat.variables.cumeq > 0.0)
        end
        @test !iszero(mat.variables.stress[1,1]) && !iszero(mat.variables.stress[1,2])
        @test isapprox(tovoigt(mat.variables.stress)[2:5], zeros(4); atol=1e-8)
    end
    @test any(plastic_flow_occurred)
end


text:
using Parallel
using Documenter

DocMeta.setdocmeta!(Parallel, :DocTestSetup, :(using Parallel); recursive=true)

makedocs(;
    modules=[Parallel],
    authors="William S. Moses <gh@wsmoses.com> and contributors",
    repo="https://github.com/wsmoses/Parallel.jl/blob/{commit}{path}#{line}",
    sitename="Parallel.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://wsmoses.github.io/Parallel.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/wsmoses/Parallel.jl",
    devbranch="main",
)


text:
using ITensors
using LinearAlgebra
using Random

include(joinpath(ITensors.examples_dir(), "src", "electronk.jl"))
include(joinpath(ITensors.examples_dir(), "src", "hubbard.jl"))

function main(; Nx::Int = 6, Ny::Int = 3, U::Float64 = 4.0, t::Float64 = 1.0,
                maxdim::Int = 3000, conserve_ky = true, use_splitblocks = true,
                nsweeps = 10, blas_num_threads = 1, strided_num_threads = 1,
                use_threaded_blocksparse = true, outputlevel = 1,
                seed = 1234)
  Random.seed!(seed)
  ITensors.Strided.set_num_threads(strided_num_threads)
  BLAS.set_num_threads(blas_num_threads)
  if use_threaded_blocksparse
    ITensors.enable_threaded_blocksparse()
  else
    ITensors.disable_threaded_blocksparse()
  end

  if outputlevel > 0
    @show Threads.nthreads()
    @show Sys.CPU_THREADS
    @show ITensors.blas_get_num_threads()
    @show ITensors.Strided.get_num_threads()
    @show ITensors.using_threaded_blocksparse()
    println()
  end

  N = Nx * Ny

  sweeps = Sweeps(nsweeps)
  maxdims = min.([100, 200, 400, 800, 2000, 3000, maxdim], maxdim)
  maxdim!(sweeps, maxdims...)
  cutoff!(sweeps, 1e-6)
  noise!(sweeps, 1e-6, 1e-7, 1e-8, 0.0)

  if outputlevel > 0
    @show sweeps
  end

  sites = siteinds("ElecK", N; conserve_qns = true,
                   conserve_ky = conserve_ky, modulus_ky = Ny)

  ampo = hubbard(Nx = Nx, Ny = Ny, t = t, U = U, ky = true)
  H = MPO(ampo, sites)

  if outputlevel > 0
    @show use_splitblocks
  end

  # This step makes the MPO more sparse but also
  # introduces more blocks.
  # It generally improves DMRG performance
  # at large bond dimensions.
  if use_splitblocks
    H = splitblocks(linkinds, H)
  end

  # Number of structural nonzero elements in a bulk
  # Hamiltonian MPO tensor
  if outputlevel > 0
    @show nnz(H[end÷2])
    @show nnzblocks(H[end÷2])
  end

  # Create start state
  state = Vector{String}(undef, N)
  for i in 1:N
    x = (i - 1) ÷ Ny
    y = (i - 1) % Ny
    if x % 2 == 0
      if y % 2 == 0
        state[i] = "Up"
      else
        state[i] = "Dn"
      end
    else
      if y % 2 == 0
        state[i] = "Dn"
      else
        state[i] = "Up"
      end
    end
  end

  psi0 = randomMPS(sites, state, 10)

  energy, psi = @time dmrg(H, psi0, sweeps; outputlevel = outputlevel)

  if outputlevel > 0
    @show Nx, Ny
    @show t, U
    @show flux(psi)
    @show maxlinkdim(psi)
    @show energy
  end
  return nothing
end

println("################################")
println("Compilation")
println("################################")
println("Without threaded block sparse:\n")
main(nsweeps = 2, use_threaded_blocksparse = false, outputlevel = 0)
println()
println("With threaded block sparse:\n")
main(nsweeps = 2, use_threaded_blocksparse = true, outputlevel = 0)
println()

println("################################")
println("Runtime")
println("################################")
println()
println("Without threaded block sparse:\n")
main(nsweeps = 10, use_threaded_blocksparse = false)
println()
println("With threaded block sparse:\n")
main(nsweeps = 10, use_threaded_blocksparse = true)
println()


text:
### Fitting Routine through (penalized MLE)

struct EFMarginalDistribution{EFD,EB,S} <: Distribution{Univariate,S}
    efd::EFD
    Z::EB
end

function Base.show(io::IO, ef::EFMarginalDistribution)
    println(io, "MarginalDistribution")
    println(io, " Prior")
    println(io, ef.efd)
    println(io, " EBayes Sample")
    return print(io, ef.Z)
end

function EFMarginalDistribution(
    efd::EFD, Z::EB
) where {EFD,EB<:Empirikos.ContinuousEBayesSample}
    return EFMarginalDistribution{EFD,EB,Continuous}(efd, Z)
end

function EFMarginalDistribution(
    efd::EFD, Z::EB
) where {EFD,EB<:Empirikos.DiscreteEBayesSample}
    return EFMarginalDistribution{EFD,EB,Discrete}(efd, Z)
end

function Distributions.pdf(efm::EFMarginalDistribution, x::Real)
    efd = efm.efd
    Z = efm.Z
    _f =
        μ ->
            pdf(efd, μ; include_basemeasure=false) *
            pdf(Empirikos.likelihood_distribution(Z, μ), x)
    return efd.integrator(_f)
end

function Distributions.logpdf(efm::EFMarginalDistribution, x::Real)
    return log(Distributions.pdf(efm, x))
end

function Distributions.cdf(efm::EFMarginalDistribution, x::Real)
    efd = efm.efd
    Z = efm.Z
    _f =
        μ ->
            pdf(efd, μ; include_basemeasure=false) *
            cdf(Empirikos.likelihood_distribution(Z, μ), x)
    return efd.integrator(_f)
end

function Empirikos.marginalize(Z::EBayesSample, efd::ExponentialFamilyDistribution)
    return EFMarginalDistribution(efd, Z)
end

Base.@kwdef struct PenalizedMLE{EF<:ExponentialFamily,T<:Real}
    ef::EF
    c0::T = 0.01
    solver = NewtonTrustRegion()
    optim_options = Optim.Options(; show_trace=true, show_every=1, g_reltol=1e-8)
    initializer = LindseyMethod(; ef=ef)
end

Base.@kwdef mutable struct FittedPenalizedMLE{P<:PenalizedMLE,T,VT<:AbstractVector{T},EFD}
    pen::P
    α_opt::VT
    α_bias = zero(α_opt)
    efd::EFD = pen.ef(α_opt)
    α_covmat = nothing
    nll_hessian = nothing
    nll_gradient = nothing
    regularizer_hessian = nothing
    regularizer_gradient = nothing
    fitter = nothing
end

Base.broadcastable(fitted::FittedPenalizedMLE) = Ref(fitted)

StatsBase.fit(fitted::FittedPenalizedMLE, args...; kwargs...) = fitted

function StatsBase.fit(pen::PenalizedMLE, Zs::Empirikos.VectorOrSummary)
    return _fit(pen, Zs)
end

function _fit(pen::PenalizedMLE, Zs)
    @unpack ef, c0, solver, optim_options, initializer = pen
    n = nobs(Zs)

    if isa(initializer, LindseyMethod)
        # initialize method through Lindsey's method
        # this mostly makes sense for normal Zs
        fit_init = StatsBase.fit(initializer, response.(Zs), StatsBase.weights(Zs))
        α_init = fit_init.α
    else
        α_init = initializer
    end

    # set up objective
    function _nll(α)
        efd = ef(α)
        return -loglikelihood(Zs, efd) / n
    end
    _s(α) = c0 * norm(α) / n
    _penalized_nll(α) = _nll(α) + _s(α)

    # ready to descend
    optim_res = optimize(_penalized_nll, α_init, solver, optim_options; autodiff=:forward)

    α_opt = Optim.minimizer(optim_res)
    @show α_opt
    #---------------------------------------
    # Extract curvature and bias information
    #---------------------------------------
    hessian_storage_nll = DiffResults.HessianResult(α_opt)
    ForwardDiff.hessian!(hessian_storage_nll, _nll, α_opt)

    nll_hessian = DiffResults.hessian(hessian_storage_nll)
    nll_gradient = DiffResults.gradient(hessian_storage_nll)

    # project onto psd cone
    nll_hessian_eigen = eigen(Symmetric(nll_hessian))
    nll_hessian_eigen.values .= max.(nll_hessian_eigen.values, 0.0)
    nll_hessian = Matrix(nll_hessian_eigen)

    hessian_storage_s = DiffResults.HessianResult(α_opt)
    ForwardDiff.hessian!(hessian_storage_s, _s, α_opt)

    regularizer_hessian = DiffResults.hessian(hessian_storage_s)
    regularizer_gradient = DiffResults.gradient(hessian_storage_s)

    inv_mat = inv((nll_hessian + regularizer_hessian) .* sqrt(n))
    α_covmat = inv_mat * nll_hessian * inv_mat

    α_bias = -inv_mat * regularizer_gradient .* sqrt(n)

    return FittedPenalizedMLE(;
        pen=pen,
        α_opt=α_opt,
        α_covmat=α_covmat,
        α_bias=α_bias,
        nll_hessian=nll_hessian,
        nll_gradient=nll_gradient,
        regularizer_hessian=regularizer_hessian,
        regularizer_gradient=regularizer_gradient,
        fitter=optim_res,
    )
end

function target_bias_std(
    fcef::FittedPenalizedMLE, target::Empirikos.EBayesTarget; bias_corrected=true, clip=true
)
    _fun(α) = target(fcef.pen.ef(α))
    α_opt = fcef.α_opt
    target_gradient = ForwardDiff.gradient(_fun, α_opt)
    target_bias = LinearAlgebra.dot(target_gradient, fcef.α_bias)
    target_variance = LinearAlgebra.dot(target_gradient, fcef.α_covmat, target_gradient)
    target_value = bias_corrected ? _fun(α_opt) - target_bias : _fun(α_opt)

    if clip
        target_value = clamp(target_value, extrema(target)...)
    end

    return (
        estimated_target=target_value,
        estimated_bias=target_bias,
        estimated_std=sqrt(target_variance),
    )
end

function confint(
    fcef::FittedPenalizedMLE,
    target::Empirikos.EBayesTarget;
    level::Real=0.95,
    clip=true,
    bias_corrected=true,
)
    α = 1 - level
    res = target_bias_std(fcef, target; clip=false, bias_corrected=bias_corrected)
    _estimate = res[:estimated_target]
    _std = res[:estimated_std]

    q_mult = quantile(Normal(), 1 - α / 2) * _std

    L, U = _estimate .+ (-1, 1) .* q_mult

    if clip
        L, U = clamp.((L, U), extrema(target)...)
    end

    return Empirikos.BiasVarianceConfidenceInterval(;
        target=target,
        method=nothing,
        α=α,
        se=_std,
        estimate=_estimate,
        lower=L,
        upper=U,
        halflength=q_mult,
    )
end

function confint(
    pmle::Union{PenalizedMLE,FittedPenalizedMLE},
    target::Empirikos.EBayesTarget,
    args...;
    kwargs...,
)
    _fit = StatsBase.fit(pmle, args...)
    return confint(_fit, target; kwargs...)
end

function Base.broadcasted_kwsyntax(
    ::typeof(confint),
    pmle::PenalizedMLE,
    targets::AbstractArray,
    args...;
    level=0.95,
    kwargs...,
)
    _fit = StatsBase.fit(pmle, args...)
    _first_ci = confint(_fit, targets[1]; level=level, kwargs...)
    confint_vec = fill(_first_ci, length(targets))
    for (index, target) in enumerate(targets[2:end])
        confint_vec[index + 1] = StatsBase.confint(_fit, target; level=level)
    end
    return confint_vec
end

function Base.broadcasted(::typeof(confint), pmle::PenalizedMLE, targets, args...)
    return confint.(pmle, targets, args...; level=0.95)
end


text:
using Plasmo
using Plots, PlasmoPlots
gr()

graph = OptiGraph()

@optinode(graph,n1)
@variable(n1, y >= 2)
@variable(n1,x >= 0)
@constraint(n1,x + y >= 3)
@objective(n1, Min, y)

@optinode(graph,n2)
@variable(n2, y)
@variable(n2,x >= 0)
@constraint(n2,x + y >= 3)
@objective(n2, Min, y)

@optinode(graph,n3)
@variable(n3, y )
@variable(n3,x >= 0)
@constraint(n3,x + y >= 3)
@objective(n3, Min, y)

@linkconstraint(graph, n1[:x] + n2[:x] + n3[:x] == 3)

plt_graph1 = layout_plot(graph,node_labels = true,markersize = 60,labelsize = 30,
linewidth = 4,layout_options = Dict(:tol => 0.01,:iterations => 2));

plt_matrix1 = matrix_plot(graph,node_labels = true,markersize = 30);


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 10b19f7a-1f8b-11eb-128e-11743ce72ef5
using PlutoUI, Plots

# ╔═╡ e5a688ac-1f91-11eb-104a-2753c9263f14
md"""
 # Floating Point Arithmetic
"""

# ╔═╡ 2df31082-1f8b-11eb-145a-93fa6e01e345
@bind f Slider(1:100)

# ╔═╡ 37c0afe8-1f8b-11eb-1154-31f5954ca4ff
f

# ╔═╡ 8d2203e8-1f8a-11eb-0332-7d512d831ec1
bitstring(f)

# ╔═╡ 4aa7bc3c-1f8b-11eb-30dd-89d64f42fc27
length(bitstring(f))

# ╔═╡ c90f7b06-1f8a-11eb-17d7-35a09af858f1
float(f)

# ╔═╡ d4cd33f2-1f8a-11eb-1af7-bdbc1fb236da
bitstring(float(f))

# ╔═╡ 4f354440-1f8b-11eb-225a-f55552bc9e29
length(bitstring(float(f)))

# ╔═╡ b7eb6210-1f9e-11eb-2a57-95b4fab4187f
md"""
IEEE 64 bit representation (double precision) (64 bits0
"""

# ╔═╡ 2ee2a8ec-1f8e-11eb-23e6-9d5e2bb00d1e
ieee0(x::Float64) = [bitstring(x)[1:1], bitstring(x)[2:12], bitstring(x)[13:64]]

# ╔═╡ d6b5d1f6-1f9e-11eb-0c32-9b9cb51352a0
md"""
## Let me work with numbers from 1 to 2
"""

# ╔═╡ e0cd119a-1f9e-11eb-1181-3d3fa6e0653f
 nextfloat(1.0)-1.0

# ╔═╡ f3357c5a-1f9e-11eb-270d-d3b5fd6bd542
2.0 ^(-52)

# ╔═╡ fbe68df8-1f9e-11eb-1278-a102df99c88c

# ╔═╡ 0bf5f86a-1f8f-11eb-17fd-75f4b45f16d5
@bind j Slider(1: 2^(52), show_value = true)

# ╔═╡ 2ff48c68-1f8f-11eb-0bd0-cde93f2cc18b
g = 1 + j/2^52

# ╔═╡ 316b7dd0-1f8e-11eb-0da9-9fc5b63a4e80
ieee0(g)

# ╔═╡ 56e8258a-1f8f-11eb-3249-6592feed84b3
ieee0(1.0)

# ╔═╡ 5d7dcbe0-1f8f-11eb-05e9-1b2ae35cf681
ieee0(-1.0)

# ╔═╡ bf3bd13c-1f95-11eb-3187-a7311125be18
nextfloat(1.0)-1

# ╔═╡ d3a6bbaa-1f95-11eb-27de-e99b0e230d21
eps(1.0)

# ╔═╡ c6d813ba-1f95-11eb-1be6-55cb3ffeb358
nextfloat(1.0f0)-1

# ╔═╡ d66f2a16-1f95-11eb-3edc-4182c1a4b8ad
eps(1.0f0)

# ╔═╡ d31baed4-1f95-11eb-19df-15ff5c378487
eps(1.0f0) * 2^23

# ╔═╡ 81a7921c-1f8f-11eb-33c6-397d38b57a75
ieee(x) = [ parse.(Int,ieee0(x),base=2)      ieee0(x)]

# ╔═╡ 85c8f822-1f8f-11eb-15a1-7956d269779a
ieee(g)

# ╔═╡ dd3031d8-1f9f-11eb-39b2-cb7ac57c6bdf
ieee(1.0)

# ╔═╡ e4b2733a-1f9f-11eb-30ac-0dc5305b8938
ieee(-1.0)

# ╔═╡ ff15233a-1f9f-11eb-0012-f198e05abaf9
md"""
# s, e, m  (sign, exponent, mantissa)  =  (-1)^s  * 2^(e-1023) * (1 + m/2^52)
"""

# ╔═╡ 3f36e4a8-1fa0-11eb-1c46-97ff95905b41
ieee(prevfloat(1.0))

# ╔═╡ f8612fa6-1f91-11eb-36d6-cbd0f9b5b0ea
float2ieee(x::Float64) =  parse.(Int,ieee0(x),base=2)

# ╔═╡ 66f85a94-1fa0-11eb-32bd-178184d8aef9
2^52-1

# ╔═╡ 0a44c782-1f92-11eb-1b67-eb93f1dfe1ff
float2ieee(1.0)

# ╔═╡ 37039e06-1f92-11eb-192e-e59ffb46a368
begin
	i =-2:2
	x =  2.0 .^i
	plot(yticks=[],xticks=round.(x,digits=2),size=(500,100))
	for xx∈x

		z = (1+1/16):(1/16):(2-1/16)
		for zz∈ xx.*z

	  plot!( [zz, zz], [-.1, 0.0], color=:red, ratio=10, legend=false, ylims=(-.1,.1),tickfontsize=1)

		end
		  plot!( [xx, xx], [-.1, .1], color=:blue, ratio=5, legend=false, ylims=(-.1,.1),tickfontsize=4)
	end

	plot!()
end

# ╔═╡ c51ad4de-1fa1-11eb-3e01-fb46ed72f493
(1 + 1.000001*eps(1.0)/2)-1

# ╔═╡ 79e5e5ca-1fa0-11eb-087f-9fad526c1b0e
ieee(0.0) ## does not follow the same rule

# ╔═╡ 8fd89a82-1fa0-11eb-2921-0d3c31c93cbc
ieee(-0.0) ## is this the same as 0.0?

# ╔═╡ a207d416-1fa0-11eb-3bf0-bdc40925736e
1.0 / 0.0

# ╔═╡ c7473096-1fa0-11eb-39fd-31c4c4b1eb04
ieee(Inf)

# ╔═╡ d08a02dc-1fa0-11eb-3369-a7c84dad581a
ieee(1.0/-0.0)

# ╔═╡ f072cf66-1fa0-11eb-0789-733f904d41d0
0.0 / 0.0

# ╔═╡ ff8975cc-1fa0-11eb-1fa6-514db7898e52
md"""
# Not a Number
"""

# ╔═╡ 0b340f86-1fa1-11eb-37bd-839bd7525dc9
NaN == NaN

# ╔═╡ 0d1c1902-1f96-11eb-0734-1bb3dda2aa7d
md"""
Why do we need this floating point representation?
"""

# ╔═╡ a40c4352-1f94-11eb-263e-632828cb6079
float2ieee(0.0)

# ╔═╡ b7b6a198-1f94-11eb-0093-8f3ab007d770
float2ieee(-0.0)

# ╔═╡ bb04eb2c-1f94-11eb-0a4a-3d2765531e94
float2ieee(Inf)

# ╔═╡ c0e4f8f2-1f94-11eb-1305-d9028395ab2f
float2ieee(-Inf)

# ╔═╡ c500ec66-1f94-11eb-079a-1b5b4befade2
1/-0.0

# ╔═╡ 1e37c20a-1f95-11eb-37f9-11c5a6041e1b
Inf - Inf

# ╔═╡ 25be9fbc-1f95-11eb-294a-1dd62a135bfa
0.0/0.0

# ╔═╡ 7345c08e-1f96-11eb-0bf1-c50d0108facd
md"""
$$\begin{matrix}   \text{exponent\\mantissa} & 0 & ≠0  \\ \hline 0 & ±0  &  \text{denormal} \\
\text{regular} &  2^k & \text{general} \\ \text{max} &  \pm\infty &  \text{NaN} \end{matrix}$$
"""

# ╔═╡ b22dc266-1f9a-11eb-199b-170ef50b702c
md"""
Rules of IEEE arithmetic
"""

# ╔═╡ 3f2943c6-1f95-11eb-0342-6569169f32ca
Float32(1.0)

# ╔═╡ 47153ebe-1f95-11eb-2224-716d2e152db0
bitstring(Float32(1.0))

# ╔═╡ 60135b2e-1f95-11eb-293d-4f8a96a677b2
float32ieee0(x::Float32) = [bitstring(x)[1:1], bitstring(x)[2:9], bitstring(x)[10:32]]

# ╔═╡ 77de76f8-1f95-11eb-00bc-915dd7977c75
float32ieee0(1.0f0)

# ╔═╡ 65ef756c-1f96-11eb-2546-35b37cb7f54c
## The rules of +, -, / , * and √

# ╔═╡ f707de2e-1fa1-11eb-0147-3fca8c03ab5f

# ╔═╡ 2b0265d4-1f96-11eb-1553-ef2ecb7d451b
md"""
What precision do you need?
 - Small absolute precision?
 - Small relative precision?
 - Neither?
"""

# ╔═╡ c5f5c97c-1f97-11eb-3683-55c16be8274e
md"""
What is displayed? What can you represent?
"""

# ╔═╡ b86c366c-1f95-11eb-0a18-a52cdf00da4c
1/5

# ╔═╡ d400e394-1f97-11eb-0f74-476c32bdd9d0
bitstring(1/5)

# ╔═╡ f5e3b162-1f97-11eb-2879-f7c7954143a0
big(1/5)

# ╔═╡ 000c656c-1f98-11eb-0f0e-657b9651f0bd
big(1)/5

# ╔═╡ 149b4c9e-1f98-11eb-19b7-eb9988c9333a
BigFloat(1//5,precision=1000)

# ╔═╡ 3aae6bac-1f98-11eb-22a5-bdc5e5de3f4c
BigFloat(π,precision=1000)

# ╔═╡ 43bed31c-1f98-11eb-1131-f5dd8dd27086
BigFloat(float(pi),precision=1000)

# ╔═╡ 56fb9c94-1f98-11eb-216e-f9e6e04d6f45
md"""
Note that every IEEE number terminates (is not a repeating expansion)
"""

# ╔═╡ 32cf8afa-1f99-11eb-0193-e969cca12f86
md"""
Find the first floating point number bigger than 1 such that x * (1/x) ≠ 1
"""

# ╔═╡ 51baa858-1f99-11eb-3329-4d9bf850b1c0
begin
	t = 1+rand()
	t,ieee(t*1/t)
end

# ╔═╡ 64c9235e-1f99-11eb-24f3-d940be06c7c8
let
	 α = 1.0
     while( α * (1/α) == 1.0)
		α = nextfloat(α)
	 end
	 Int((α-1)/eps(1.0))
end

# ╔═╡ ecc5c55a-1f99-11eb-3f44-4161850426cc
α= 1 + 257736490*eps()

# ╔═╡ 40575ec2-1f9a-11eb-2d29-5b517d1f3b18
α*(1/α)

# ╔═╡ 0b0719c2-1f9c-11eb-3ed6-6927c927145d

# ╔═╡ 11ebb920-1f9c-11eb-061a-8df3899e8945

# ╔═╡ Cell order:
# ╟─e5a688ac-1f91-11eb-104a-2753c9263f14
# ╠═2df31082-1f8b-11eb-145a-93fa6e01e345
# ╠═37c0afe8-1f8b-11eb-1154-31f5954ca4ff
# ╠═8d2203e8-1f8a-11eb-0332-7d512d831ec1
# ╠═4aa7bc3c-1f8b-11eb-30dd-89d64f42fc27
# ╠═c90f7b06-1f8a-11eb-17d7-35a09af858f1
# ╠═d4cd33f2-1f8a-11eb-1af7-bdbc1fb236da
# ╠═4f354440-1f8b-11eb-225a-f55552bc9e29
# ╟─b7eb6210-1f9e-11eb-2a57-95b4fab4187f
# ╠═2ee2a8ec-1f8e-11eb-23e6-9d5e2bb00d1e
# ╟─d6b5d1f6-1f9e-11eb-0c32-9b9cb51352a0
# ╠═e0cd119a-1f9e-11eb-1181-3d3fa6e0653f
# ╠═f3357c5a-1f9e-11eb-270d-d3b5fd6bd542
# ╠═fbe68df8-1f9e-11eb-1278-a102df99c88c
# ╠═0bf5f86a-1f8f-11eb-17fd-75f4b45f16d5
# ╠═2ff48c68-1f8f-11eb-0bd0-cde93f2cc18b
# ╠═316b7dd0-1f8e-11eb-0da9-9fc5b63a4e80
# ╠═56e8258a-1f8f-11eb-3249-6592feed84b3
# ╠═5d7dcbe0-1f8f-11eb-05e9-1b2ae35cf681
# ╠═bf3bd13c-1f95-11eb-3187-a7311125be18
# ╠═d3a6bbaa-1f95-11eb-27de-e99b0e230d21
# ╠═c6d813ba-1f95-11eb-1be6-55cb3ffeb358
# ╠═d66f2a16-1f95-11eb-3edc-4182c1a4b8ad
# ╠═d31baed4-1f95-11eb-19df-15ff5c378487
# ╠═81a7921c-1f8f-11eb-33c6-397d38b57a75
# ╠═85c8f822-1f8f-11eb-15a1-7956d269779a
# ╠═dd3031d8-1f9f-11eb-39b2-cb7ac57c6bdf
# ╠═e4b2733a-1f9f-11eb-30ac-0dc5305b8938
# ╠═ff15233a-1f9f-11eb-0012-f198e05abaf9
# ╠═3f36e4a8-1fa0-11eb-1c46-97ff95905b41
# ╠═f8612fa6-1f91-11eb-36d6-cbd0f9b5b0ea
# ╠═66f85a94-1fa0-11eb-32bd-178184d8aef9
# ╠═0a44c782-1f92-11eb-1b67-eb93f1dfe1ff
# ╠═37039e06-1f92-11eb-192e-e59ffb46a368
# ╠═c51ad4de-1fa1-11eb-3e01-fb46ed72f493
# ╠═79e5e5ca-1fa0-11eb-087f-9fad526c1b0e
# ╠═8fd89a82-1fa0-11eb-2921-0d3c31c93cbc
# ╠═a207d416-1fa0-11eb-3bf0-bdc40925736e
# ╠═c7473096-1fa0-11eb-39fd-31c4c4b1eb04
# ╠═d08a02dc-1fa0-11eb-3369-a7c84dad581a
# ╠═f072cf66-1fa0-11eb-0789-733f904d41d0
# ╟─ff8975cc-1fa0-11eb-1fa6-514db7898e52
# ╠═0b340f86-1fa1-11eb-37bd-839bd7525dc9
# ╟─0d1c1902-1f96-11eb-0734-1bb3dda2aa7d
# ╠═a40c4352-1f94-11eb-263e-632828cb6079
# ╠═b7b6a198-1f94-11eb-0093-8f3ab007d770
# ╠═bb04eb2c-1f94-11eb-0a4a-3d2765531e94
# ╠═c0e4f8f2-1f94-11eb-1305-d9028395ab2f
# ╠═c500ec66-1f94-11eb-079a-1b5b4befade2
# ╠═1e37c20a-1f95-11eb-37f9-11c5a6041e1b
# ╠═25be9fbc-1f95-11eb-294a-1dd62a135bfa
# ╠═7345c08e-1f96-11eb-0bf1-c50d0108facd
# ╟─b22dc266-1f9a-11eb-199b-170ef50b702c
# ╠═3f2943c6-1f95-11eb-0342-6569169f32ca
# ╠═47153ebe-1f95-11eb-2224-716d2e152db0
# ╠═60135b2e-1f95-11eb-293d-4f8a96a677b2
# ╠═77de76f8-1f95-11eb-00bc-915dd7977c75
# ╠═65ef756c-1f96-11eb-2546-35b37cb7f54c
# ╠═f707de2e-1fa1-11eb-0147-3fca8c03ab5f
# ╠═2b0265d4-1f96-11eb-1553-ef2ecb7d451b
# ╟─c5f5c97c-1f97-11eb-3683-55c16be8274e
# ╠═b86c366c-1f95-11eb-0a18-a52cdf00da4c
# ╠═d400e394-1f97-11eb-0f74-476c32bdd9d0
# ╠═f5e3b162-1f97-11eb-2879-f7c7954143a0
# ╠═000c656c-1f98-11eb-0f0e-657b9651f0bd
# ╠═149b4c9e-1f98-11eb-19b7-eb9988c9333a
# ╠═3aae6bac-1f98-11eb-22a5-bdc5e5de3f4c
# ╠═43bed31c-1f98-11eb-1131-f5dd8dd27086
# ╟─56fb9c94-1f98-11eb-216e-f9e6e04d6f45
# ╟─32cf8afa-1f99-11eb-0193-e969cca12f86
# ╠═51baa858-1f99-11eb-3329-4d9bf850b1c0
# ╠═64c9235e-1f99-11eb-24f3-d940be06c7c8
# ╠═ecc5c55a-1f99-11eb-3f44-4161850426cc
# ╠═40575ec2-1f9a-11eb-2d29-5b517d1f3b18
# ╠═0b0719c2-1f9c-11eb-3ed6-6927c927145d
# ╠═11ebb920-1f9c-11eb-061a-8df3899e8945
# ╠═10b19f7a-1f8b-11eb-128e-11743ce72ef5

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
## Produces the logstic plot
# Requires PyPlot

function logisticPlot()
  tic()
  r = 2.9:.00005:4; numAttract = 100;
  steady = ones(length(r),1)*.25;
  for i=1:400 ## Get to steady state
    @devec steady = r.*steady.*(1-steady);
  end
  x = zeros(length(steady),numAttract);
  x[:,1] = steady;
  @simd for i=2:numAttract ## Grab values at the attractor
    @inbounds @fastmath x[:,i] = r.*x[:,i-1].*(1-x[:,i-1]);
  end
  toc()
  fig = figure(figsize=(20,10));
  plot(collect(r),x,"b.",markersize=.06)
  savefig("plot.png",dpi=300);
end
using PyPlot
@time logisticPlot()


text:
# This file is a part of EncodedArrays.jl, licensed under the MIT License (MIT).

import Test
Test.@testset "Package EncodedArrays" begin

include("test_encoded_array.jl")
include("test_varlen_io.jl")
include("test_varlen_diff_codec.jl")

end # testset


text:
## DOC STRING GENERATION

function detailed_doc_string(M; header="", example="", body="")
    ret = "    $header"
    ret *= "\n\n"
    ret *= "An iteration control, as in, `$example`. "
    ret *= "\n\n"
    ret *= body
    return ret
end

_err_create_docs() = error(
"@create_docs syntax error. Usage: \n"*
    "@create_docs(ControlSubType, header=..., example=..., body=...)")

macro create_docs(M_ex, exs...)
    M_ex isa Symbol || _err_create_docs()
    h = ""
    e = ""
    b = ""
    for ex in exs
        ex.head == :(=) || _err_create_docs()
        ex.args[1] == :header &&  (h = ex.args[2])
        ex.args[1] == :example && (e = ex.args[2])
        ex.args[1] == :body &&    (b = ex.args[2])
    end
    esc(quote
        "$(IterationControl.detailed_doc_string($M_ex, header=$h, example=$e, body=$b))"
        function $M_ex end
        end)
end

## FOR EVALUATION OF FUNCTIONS PASSED TO LOGGING CONTROLS

_log_eval(f, model) = f(expose(model))
_log_eval(f::String, model) = f


text:
## NLEVP
"""
    Term{T}

Single term of a linear operator family.

# Fields
- `coeff::T`: matrix coefficient of the term
- `func::Tuple`: tuple of scalar-valued functions multiplying the matrix coefficient
- `symbol::String`: character string for displaying the function(s)
- `params::Tuple`: tuple of tuples containing function symbols for each function
- `operator:String`: String for displaying the matrix coefficient
- `varlist::Array{Symbol,1}`: Array containing all symbols of fucntion arguments
that are used at least once
"""
struct Term{T}
  coeff::T
  func::Tuple
  symbol::String
  params::Tuple
  operator::String
  varlist::Array{Symbol,1} #TODO: this is redundant information
end

#constructor
function Term(coeff,func::Tuple,params::Tuple,symbol::String,operator::String)
  varlist=Symbol[]
  for par in params
    for var in par
      push!(varlist,var)
    end
  end
  unique!(varlist)
  return Term(coeff,func,symbol,params,operator, varlist)
end

"""
    term=Term(coeff,func::Tuple,params::Tuple,operator::String)

Standard constructor for type `Term`.

# Arguments
- `coeff`: matrix coefficient of the term
- `func::Tuple`: tuple of scalar-valued functions multiplying the matrix coefficient
- `params::Tuple`: tuple of tuples containing function symbols for each function
- `operator:String`: String for displaying the matrix coefficient

# Notes

The rendering of the functions of `term` ist automotized. If the passed functions
implement a method `f(z::Symbol)` then this method is used, otherwise the
functions will be simply displayed with the string `f`. You can overwrite this
behavior by explicitly passing a string for the function display using the syntax

  term =Term(coeff,func::Tuple,params::Tuple,symbol::String,operator::String)

where the extra argument `symbol`is the string used to display the function.

See also: [LinearOperatorFamily](@ref)
"""
function Term(coeff,func::Tuple,params::Tuple,operator::String)
  symbol=""
  for (f,p) in zip(func,params)
    if applicable(f,p...)
      symbol*=f(p...)
    else
      symbol*="f("
      for idx = 1:length(p)-1
        symbol*="$(p[idx]),"
      end
      symbol*="$(p[end]))"
    end
  end
  return Term(coeff,func,params,symbol,operator)
end

#Solution object
"""
    Solution

Type for the solution of an iterative eigensolver.

# Fields
- `params`: Dictionary mapping parameter symbols to their values
- `v`: right (direct) eigenvector
- `v_adj`: left (adjoint) eigenvector
- `eigval`: the symbol of the parameter that is the eigenvalue
- `eigval_pert`: extra data for asymptotic series expansion of the eigenvalue
- `v_pert`: extra data for asymptotic series expansion of the eigenvector
- `auxval`: the symbol of the parameter that is the auxiliary eigenvalue

See also: [LinearOperatorFamily](@ref)
"""
mutable struct Solution #IDEA: make immutable and parametrized type
  params
  v
  v_adj
  eigval
  eigval_pert
  v_pert
  auxval
end
#constructor
"""
    sol=Solution(params,v,v_adj,eigval,auxval=Symbol())

Standard constructor for a solution object. See [Solution](@ref) for details.
"""
function Solution(params,v,v_adj,eigval,auxval=Symbol())
  return Solution(deepcopy(params),v,v_adj,eigval,Dict{Symbol,Any}(),Dict{Symbol,Any}(),auxval) #TODO: copy params?
end

"""
    LinearOperatorFamily

Type for a linear operator family. The type is mutable.

# Fields
- `terms`: Array of terms forming the operator family
- `params`: Dictionary mapping parameter symbols to their values
- `eigval`: the symbol of the parameter that is the eigenvalue
- `auxval`: the symbol of the parameter that is the auxiliary eigenvalue
- `active`: list of symbols that can be actively change when using an instance of
the family as a function.
- `mode`: mode that is controlling which terms are evalauted when calling an
instance of the family. This field is not to be modified by the user.

See also: [Solution](@ref), [Term](@ref)
"""
mutable struct LinearOperatorFamily #TODO: add example to the doc
  terms::Array{Term,1}
  params::Dict{Symbol,ComplexF64}
  eigval::Symbol
  auxval::Symbol
  active::Array{Symbol,1}
  mode::Symbol
end

#constructors
"""
    L=LinearOperatorFamily(params,values)

Standard constructor for an empty Linear operator family.

# Arguments
- `params`: List of parameter symbols
- `values`: List of corresponding parameter values

# Note
The first parameter appearing in `params` will be assigned as eigenvalue. If
there are more than 1 parameters in the list, the last parameter will be designated
as auxiliary eigenvalue. (These choices can be changed after construction)
If `values`is omitted, then all parameters will be initialized with `NaN+NaN*im`.
If also `params` is omitted the family will be initialized with `:λ` as its
eigenvalue an no auxiliary eigenvalue.

Terms can be added to the family using the `+` operator or the more memory efficient
`push!` function. For instance `L+=term` or `push!(L,term)` both add `term`
to the list of terms.

See also: [Solution](@ref), [Term](@ref)
"""
# standard constructors
function LinearOperatorFamily(params,values)
  terms=Term[]
  eigval=Symbol(params[1])
  if length(params)>1
    auxval=Symbol(params[end])
  else
    auxval=Symbol("")
  end
  active=[eigval]
  pars=Dict{Symbol,ComplexF64}()
  for (p,v) in zip(params,values) #Implement alphabetical sorting
    pars[Symbol(p)]=v
  end
  mode=:all
  return LinearOperatorFamily(terms,pars,eigval,auxval,active,mode)
end
function LinearOperatorFamily()
  return LinearOperatorFamily(["λ",],[NaN+NaN*1im,])
end
function LinearOperatorFamily(params)
  return LinearOperatorFamily(params,[NaN+NaN*1im for a in params])
end

#loader
"""
    L=LinearOperatorFamily(fname::String)

Load and construct `LinearOperatorFamily` from file `fname`.

See also: [`save`](@ref)
"""
function LinearOperatorFamily(fname::String)
  D=read_toml(fname)
  eigval=D["eigval"]
  auxval=D["auxval"]
  params=[]
  vals=[]
  for (par,val) in D["params"]
    push!(params,par)
    push!(vals,val)
  end
  L=LinearOperatorFamily(params,vals)
  L.eigval=eigval
  L.auxval=auxval
  L.active=[eigval]
  for idx =1:length(D["/terms"])
    term=D["/terms"]["/"*string(idx)]
    I=term["/sparse_matrix"]["I"]
    J=term["/sparse_matrix"]["J"]
    V=term["/sparse_matrix"]["V"]
    m, n = term["size"]
    M=sparse(I,J,V,m,n)
    push!(L,Term(M,term["functions"],term["params"],term["symbol"],term["operator"]))
  end
  return L
end

import Dates
#saver
"""
    save(fname::String,L::LinearOperatorFamily)

Save `L` to file `fname`. The file is utf8 encoded and adheres to a Julia-enriched TOML standard.

See also: [`LinearOperatorFamily`](@ref)
"""
function save(fname,L::LinearOperatorFamily)
  date=string(Dates.now(Dates.UTC))

  eq=""
  for term in L.terms
    if term.operator[1]=='_' #TODO: put this into signature?
      continue
    end
    eq*="+"*string(term)
  end

  open(fname,"w") do f
    write(f,"# LinearOperatorFamily version 0\n")
    write(f,"#"*date*"\n")
    write(f,"#"*eq*"\n")
    write(f,"params=[")
    for (key,value) in L.params
        write(f,"(:$(string(key)),$value),\n")
    end
    write(f,"]\n")

    write(f,"eigval=:$(L.eigval)\n")
    write(f,"auxval=:$(L.auxval)\n")
    #TODO activitiy and terms
    write(f,"[terms]\n")
    for (idx,term) in enumerate(L.terms)
      write(f,"\t[terms.$idx]\n")
      write(f,"\tfunctions=(")
      for func in term.func
        write(f,"$func,")
      end
      write(f,")\n")
      write(f,"\tsymbol=\"$(term.symbol)\"\n")
      write(f,"\tparams=$(term.params)\n")
      # for param in term.params
      #   write(f,"[")
      #   for par in param
      #     write(f,":$(String(par)),")
      #   end
      #   write(f,"],")
      # end
      # write(f,"]\n")
      write(f,"\toperator=\"$(term.operator)\"\n")
      m,n=size(term.coeff)
      write(f,"\tsize=[$m,$n]\n")
      #write(f,"\tis_sparse=$(SparseArrays.issparse(term.coeff))\n")
      write(f,"\t\t[terms.$idx.sparse_matrix]\n")
      I,J,V=SparseArrays.findnz(term.coeff)
      write(f,"\t\tI=$I\n")
      write(f,"\t\tJ=$J\n")
      #write(f,"\t\tV=$V\n")
      write(f,"\t\tV=Complex{Float64}[")
      for v in V
        if imag(v)>=0
          write(f,"$(real(v))+$(imag(v))im,")
        else
          write(f,"$(real(v))$(imag(v))im,")
        end
      end
      write(f,"]\n")
      write(f,"\n")
    end
  end
end

#TODO lesen
#eval funktion nutzen
#beginnt ist mit :  ????
#beginnt es mit (
#beginnt es mit [   ???
#dann eval
# es gibt nur drei arten: tags, variablen, listen/tuple rest ist eval

import Base.push!
function push!(L::LinearOperatorFamily,T::Term)
  d=Dict()
  for (idx,term) in enumerate(L.terms)
    d[(term.func, term.params)]=idx
  end
  signature=(T.func,T.params)
  if signature in keys(d) #change existing term if signature known
    idx=d[signature]
    coeff=L.terms[idx].coeff+T.coeff
    if LinearAlgebra.norm(coeff)==0
      deleteat!(L.terms,idx) #delte if resulting coeff is 0
      #check for unbound variables and delete
      vars=[]
      for term in L.terms
        for pars in term.params
          for par in pars
            if par ∉ vars
              push!(vars,par)
            end
          end
        end
      end
      for par in keys(T.params)
        if par ∉ vars
          delete!(L.params,par)
        end
      end

    else
      L.terms[idx]=Term(coeff,L.terms[idx].func,L.terms[idx].symbol,L.terms[idx].params,L.terms[idx].operator,L.terms[idx].varlist) #overwrite term if coeff is non-zero
    end
  else #add term if signature is new
    for pars in T.params
      for par in pars
        if par ∉ keys(L.params)
          L.params[par]=NaN+NaN*1im #initialize variable if its new
        end
      end
    end
    push!(L.terms,T)
  end
end
# function push!(a::LinearOperatorFamily,b::LinearOperatorFamily)
#    push!(a.terms,b.terms)
# end
#
#
import Base.(+)
function (+)(a::LinearOperatorFamily,b::Term)
  L=deepcopy(a)
  push!(L,b)
  return L
end
function (+)(b::Term,a::LinearOperatorFamily,)
  L=deepcopy(a)
  push!(L,b)
  return L
end
import Base.(-)
function (-)(a::LinearOperatorFamily,b::Term)
  L=deepcopy(a)
  push!(L,Term(-b.coeff,b.func,b.symbol,b.params,b.operator,b.varlist))
  return L
end
function (-)(b::Term,a::LinearOperatorFamily)
  L=deepcopy(a)
  for i=1:length(L.terms)
    L.terms[i].coeff*=-1
  end
  push!(L,b)
  return L
end

import Base.size
"""
  d=size(L::LinearOperatorFamily)

Get dimension `d` of `d×d` linear operator family `L`.
"""
function size(L::LinearOperatorFamily)
  if length(L.terms)!=0
    shape=size(L.terms[1].coeff)
  else
    shape=(0,0)
  end

  return shape[1]
end

#nice display of LinearOperatorFamilies in interactive and other modes
import Base.show
function show(io::IO,L::LinearOperatorFamily)
  if !isempty(L.terms)
    shape=size(L.terms[1].coeff)
    if length(shape)==2
      txt="$(shape[1])×$(shape[2])-dimensional operator family: \n\n"
    else
      txt="$(shape[1])-dimensional vector family: \n\n"
    end
  else
    txt="empty operator family\n\n"
  end

  eq=""
  for term in L.terms
    if term.operator[1]=='_'
      continue
    end
    eq*="+"*string(term)
  end

  parameter_list="\n\nParameters\n----------\n"
  for (par,val) in L.params
    parameter_list*=string(par)*"\t"*string(val)*"\n"
  end
  print(io, txt*eq[2:end]*parameter_list)
end

import Base.string
function string(T::Term)
  if T.symbol==""
    txt=""
  else
    txt=string(T.symbol)*"*"
  end
  txt*=T.operator
end
function show(io::IO,T::Term)
  print(io,string(T))
end

#make solution showable
function string(sol::Solution)
  txt="""####Solution####
  eigval:
  $(sol.eigval) = $(sol.params[sol.eigval])

  Parameters:
  """
  for (key,val) in sol.params
    if key!=sol.eigval && key!=sol.auxval
      txt*="$key = $val\n"
    end
  end

  if sol.auxval in keys(sol.params)
    txt*="""

    Residual:
    abs($(sol.auxval)) = $(abs(sol.params[sol.auxval]))
    """
  end
  return txt
end
function show(io::IO,sol::Solution)
  print(io,string(sol))
end

#make terms callable
function (term::Term)(dict::Dict{Symbol,Tuple{ComplexF64,Int64}})
  coeff::ComplexF64=1 #TODO parametrize type
  for (func,pars) in zip(term.func, term.params)
    args=[]
    dargs=[]
    for par in pars
      push!(args,dict[par][1])
      push!(dargs,dict[par][2])
    end

    coeff*=func(args...,dargs...)
  end
  return coeff*term.coeff
end

#make LinearOperatorFamily callable
function(L::LinearOperatorFamily)(args...;oplist=[],in_or_ex=false)
  if L.mode==:all # if mode is all the first n args correspond to the values of the n active variables
    for (var,val) in zip(L.active,args)
      L.params[var]=val #change the active variables
    end
  end
  if L.mode==:all && length(args)==length(L.active)
    derivs=zeros(Int64,length(L.active))
  else
    derivs=args[end-length(L.active)+1:end] #TODO implement sanity check on length of args
  end

  derivDict=Dict{Symbol,Int64}()
  for (var, drv) in zip(L.active,derivs)
    derivDict[var]=drv #create a dictionairy for the active variables
  end

  coeff=spzeros(size(L.terms[1].coeff)...)  #TODO: improve this to handle more flexible matrices

  for term in L.terms
    if (!in_or_ex && term.operator in oplist) || (in_or_ex && !(term.operator in oplist)) || (L.mode!=:householder && term.operator=="__aux__")  #TODO: consider deprecating this feature together with nicoud and picard
      continue
    end

    #check whether term is constant w.r.t. to some parameter then deriv is 0 thus continue
    skip=false
    for (var,d) in zip(L.active, derivs)
      if d>0 && !(var in term.varlist)
        skip=true
        break
      end
    end
    if skip
      continue
    end
    dict=Dict{Symbol,Tuple{Complex{Float64},Int64}}()
    for var in term.varlist
      dict[var]=(L.params[var], var in keys(derivDict) ? derivDict[var] : 0)
    end
    coeff+=term(dict) #
  end

  if L.mode in [:compact,:householder]
    coeff/=prod(factorial.(float.(args[end-length(L.active)+1:end])))
  end

  return coeff
end

#wrapper to perturbation theory
#TODO: functional programming renders this obsolete, no?
"""
    perturb!(sol::Solution,L::LinearOperatorFamily,param::Symbol,N::Int; <keyword arguments>)

Compute the `N`th order power series perturbation coefficients for the solution `sol` of the nonlineaer eigenvalue problem given by the operator family `L` with respect to the parameter `param`. The coefficients will be stored in the field `sol.eigval_pert` and `sol.v_pert` for the eigenvalue and the eignevector, respectively.

# Keyword Arguments
- `mode = :compact`: parameter controlling internal programm flow. Use the default, unless you know what you are doing.

# Notes
For large perturbation orders `N` the method might be slow.

See also: [`perturb_fast!`](@ref)
"""
 function perturb!(sol::Solution,L::LinearOperatorFamily,param::Symbol,N::Int;mode=:compact)
   active=L.active #TODO: copy?
   params=L.params
   L.params=sol.params
   current_mode=L.mode
   L.active=[sol.eigval, param]
   L.mode=mode
   pade_symbol=Symbol("$(string(param))/Taylor")
   sol.eigval_pert[pade_symbol],sol.v_pert[pade_symbol]=perturb(L,N,sol.v,sol.v_adj)
   sol.eigval_pert[pade_symbol][1]=sol.params[sol.eigval]
   L.active=active
   L.mode=current_mode
   L.params=params
   return
 end

"""
    perturb_fast!(sol::Solution,L::LinearOperatorFamily,param::Symbol,N::Int; <keyword arguments>)

Compute the `N`th order power series perturbation coefficients for the solution `sol` of the nonlineaer eigenvalue problem given by the operator family `L` with respect to the parameter `param`. The coefficients will be stored in the field `sol.eigval_pert` and `sol.v_pert` for the eigenvalue and the eignevector, respectively.

# Keyword Arguments
- `mode = :compact`: parameter controlling internal programm flow. Use the default, unless you know what you are doing.

# Notes
This method reads multi-indeces for the computation of the power series coefficients from disk. Make sure that JulHoltz is properly installed to use this fast method.

See also: [`perturb!`](@ref)
"""
 function perturb_fast!(sol::Solution,L::LinearOperatorFamily,param::Symbol,N::Int;mode=:compact)
   active=L.active #TODO: copy?
   params=L.params
   L.params=sol.params
   current_mode=L.mode
   L.active=[sol.eigval, param]
   L.mode=mode
   pade_symbol=Symbol("$(string(param))/Taylor")
   sol.eigval_pert[pade_symbol],sol.v_pert[pade_symbol]=perturb_disk(L,N,sol.v,sol.v_adj)
   sol.eigval_pert[pade_symbol][1]=sol.params[sol.eigval]
   L.active=active
   L.mode=current_mode
   L.params=params
   return
 end

"""
    perturb_norm!(sol::Solution,L::LinearOperatorFamily,param::Symbol,N::Int; <keyword arguments>)

 Compute the `N`th order power series perturbation coefficients for the solution `sol` of the nonlineaer eigenvalue problem given by the operator family `L` with respect to the parameter `param`. The coefficients will be stored in the field `sol.eigval_pert` and `sol.v_pert` for the eigenvalue and the eignevector, respectively.

 # Keyword Arguments
 - `mode = :compact`: parameter controlling internal programm flow. Use the default, unless you know what you are doing.

 # Notes
 This method reads multi-indeces for the computation of the power series coefficients from disk. Make sure that JulHoltz is properly installed to use this fast method.

 See also: [`perturb!`](@ref)
 """
  function perturb_norm!(sol::Solution,L::LinearOperatorFamily,param::Symbol,N::Int;mode=:compact)
    active=L.active #TODO: copy?
    params=L.params
    L.params=sol.params
    current_mode=L.mode
    L.active=[sol.eigval, param]
    L.mode=mode
    pade_symbol=Symbol("$(string(param))/Taylor")
    sol.eigval_pert[pade_symbol],sol.v_pert[pade_symbol]=perturb_norm(L,N,sol.v,sol.v_adj)
    sol.eigval_pert[pade_symbol][1]=sol.params[sol.eigval]
    L.active=active
    L.mode=current_mode
    L.params=params
    return
  end

 #TODO: Implement solution class

function pade(ω,L,M) #TODO: harmonize symbols for eigenvalues, modes etc
  A=zeros(ComplexF64,M,M)
  for i=1:M
    for j=1:M
      if L+i-j>=0
        A[i,j]=ω[L+i-j+1] #+1 is for 1-based indexing
      end
    end
  end
  b=A\(-ω[L+2:L+M+1])
  b=[1;b]
  a=zeros(ComplexF64,L+1)
  for l=0:L
    for m=0:l
      if m<=M
        a[l+1]+=ω[l-m+1]*b[m+1] #+1 is for zero-based indexing
      end
    end
  end
  return a,b
end

function pade!(sol::Solution,param,L::Int64,M::Int64;vector=false)
  #TODO: implement sanity check whether L+M+1=length(sol.eigval_pert[:Taylor])
  pade_symbol=Symbol("$(string(param))/[$L/$M]")
  taylor_symbol=Symbol("$(string(param))/Taylor")
  sol.eigval_pert[pade_symbol]=pade(sol.eigval_pert[taylor_symbol],L,M)

  #TODO Degeneracy
  if vector
    D=length(sol.v)
    #preallocation
    A=Array{Array{ComplexF64}}(undef,L+1)
    B=Array{Array{ComplexF64}}(undef,M+1)
    for idx in 1:length(A)
      A[idx]=Array{ComplexF64}(undef,D)
    end
    for idx in 1:length(B)
      B[idx]=Array{ComplexF64}(undef,D)
    end
    for idx in 1:D
      v=Array{ComplexF64}(undef,L+M+1)
      for ord = 1:L+M+1
        v[ord]=sol.v_pert[taylor_symbol][ord][idx]
      end
      a,b=pade(v,L,M)
      for ord=1:length(A)
        A[ord][idx]=a[ord]
      end
      for ord=1:length(B)
        B[ord][idx]=b[ord]
      end
    end
    sol.v_pert[pade_symbol]=A,B
  end
  return
end

#make solution object callable
function (sol::Solution)(param::Symbol,ε,L=0,M=0; vector=false) #Todo not really performant, consider syntax that allows for vectorization in eigenvector
    pade_symbol=Symbol("$(string(param))/[$L/$M]")
    if pade_symbol ∉ keys(sol.eigval_pert) || (vector && pade_symbol ∉ keys(sol.v_pert))
      pade!(sol,param,L,M,vector=vector)
    end
    a,b=sol.eigval_pert[pade_symbol]
    Δε=ε-sol.params[param]
    eigval=polyval(a,Δε)/polyval(b,Δε)
    if !vector
      return eigval
    else
      A, B = sol.v_pert[pade_symbol]
      eigvec = polyval(A,Δε) ./ polyval(B,Δε)
      return eigval, eigvec
    end
end
#TODO: implement parameter activity in solution type

# function polyval(A,z)
#     f=zeros(eltype(A[1]),length(A[1]))
#     for (idx,a) in enumerate(A)
#         f.+=a*(z^(idx-1))
#     end
#     return f
# end

#polyval(p,z)= sum(p.*(z.^(0:length(p)-1)))  #TODO: Although, this is not performance critical. Consider some performance aspects
"""
    f=polyval(p,z)
Evaluate polynomial f(z)=∑_i p[i]z^1 at z, using Horner's scheme.
"""
function polyval(p,z)
  f=ones(size(z))*p[end]
  for i = length(p)-1:-1:1
     f.*=z
     f.+=p[i]
  end
  return f
end
function polyval(p,z::Number)
  f=p[end]
  for i = length(p)-1:-1:1
     f*=z
     f+=p[i]
  end
  return f
end

function estimate_pol(ω::Array{Complex{Float64},1})
  N=length(ω)
  Δε=zeros(ComplexF64,N-2)
  k=zeros(ComplexF64,N-2)
  for j =2:length(ω)-1
    i=j-1 # index shift to account for 1-based indexing
    denom=(i+1)*ω[j+1]*ω[j-1]-i*ω[j]^2
    Δε[i]=ω[j]*ω[j-1]/denom
    k[i]=(i^2-1)*ω[j+1]*ω[j-1]-(i*ω[j])^2
  end
  return Δε,k
end

function estimate_pol(sol::Solution,param::Symbol)
  pade_symbol=Symbol("$(string(param))/Taylor")
  return estimate_pol(sol.eigval_pert[pade_symbol])
end

function conv_radius(a::Array{Complex{Float64},1})
  N=length(a)
  r=zeros(Float64,N-1)
  for n=1:N-1
    r[n]=abs(a[n]/a[n+1])
  end
  return r
end

function conv_radius(sol::Solution,param::Symbol)
  pade_symbol=Symbol("$(string(param))/Taylor")
  return conv_radius(sol.eigval_pert[pade_symbol])
end


text:
using Match

import Base: convert

"""
Specifies a particular electronic level within a species.

That is, all the quantum numbers to specify a state, except for the projection
of the total angular momentum on the quantization axis.
"""
abstract type LevelSpec end

"""
Specifies a particular electronic state using the relevant low-field
quantum numbers.

The total spin quantum number is implicitly assumed to be :math:`1 / 2` (cf.
`NoHyperfineOneElectronSpecies`).
"""
immutable NoHyperfineNumberSpec <: LevelSpec
    "Orbital angular momentum quantum number."
    l::Integer

    "Total angular momentum quantum number."
    j::Rational
end

"""
Specifies a particular electronic state using spectroscopic notation.
"""
immutable SpectroscopicSpec <: LevelSpec
    val::String
end

convert{T<:LevelSpec}(::Type{T}, s::String) = convert(T, SpectroscopicSpec(s))

function convert(::Type{NoHyperfineNumberSpec}, s::SpectroscopicSpec)
    val = s.val

    l_sym = val[1]
    l = findfirst(['S', 'P', 'D', 'F', 'G'], l_sym) - 1
    if l == -1
        throw(ArgumentError("Unknown symbol '$(s[1])' for orbital angular momentum quantum number: $(s.val)"))
    end

    val = lstrip(val[2 : end], '_')

    if val[1] == '{'
        if val[end] != '}'
            throw(ArgumentError("Invalid total angular momentum quantum number specification '$val': $(s.val)"))
        end
        val = val[2 : end - 1]
    end

    # Parse numerator.
    fraction_idx = findfirst(val, '/')
    if fraction_idx < 2
        throw(ArgumentError("Expected fractional total angular momentum quantum number, not '$val': $(s.val)"))
    end

    j_num_str = val[1 : fraction_idx - 1]
    j_num = try
        parse(Int64, j_num_str)
    catch
        throw(ArgumentError("Invalid numerator '$j_num_str' for total angular momentum quantum number: $(s.val)"))
    end

    # Skip fraction slash.
    val = val[fraction_idx + 1 : end]

    if length(val) == 0
        throw(ArgumentError("Expected denominator for total angular momentum quantum number: $(s.val)"))
    end

    # Skip optional second slash (Julia's fraction notation).
    val = lstrip(val, '/')

    # Parse denominator.
    j_den = try
        parse(Int64, val)
    catch
        throw(ArgumentError("Invalid denominator '$val' for total angular momentum quantum number: $(s.val)"))
    end

    NoHyperfineNumberSpec(l, j_num // j_den)
end

immutable StateSpec
    level::LevelSpec

    """
    Total momentum projection on the quantisation (z) axis (usually denoted
    :math:`m_J` or :math:`m_F` depending on the species).
    """
    m::Rational
end

export LevelSpec, NoHyperfineNumberSpec, convert, StateSpec, SpectroscopicSpec


text:
using Test
using SequencerJ
using StatsBase
using Random
using LightGraphs
using Images

@testset "SequencerJ.jl" begin

    # Random.seed!(2732);

    # declare a loss function to measure the error in calculations
    loss(A,B) = L2(A, B)

    imgsmall = Gray.(load(joinpath(@__DIR__,"..","resources","colony.png")))
    SMALL = convert(Matrix{Float32}, imgsmall)

    imgmed = Gray.(load(joinpath(@__DIR__,"..","resources","Aged Whisky.jpg")));
    MED = convert(Matrix{Float32}, imgmed)

    imgbig = Gray.(load(joinpath(@__DIR__,"..","resources","Hummingbird!.jpeg")))
    BIG = convert(Matrix{Float32}, imgbig)

    include("algotests.jl")

    @testset "autoscale produces one integer result" begin
        A = rand(100,100)
        @test autoscale(A) isa Int
    end

    @testset "autoscale produces one integer result" begin
        A = rand(100,100)
        @test autoscale(A) isa Int
    end

    @testset "autoscale produces a Fibonacci number" begin
        A = rand(100,100)
        @test autoscale(A) ∈  SequencerJ.FIB
    end

    @testset "autoscale p=0.2" begin
        A = rand(100,100)
        @test autoscale(A; p=0.2) ∈ SequencerJ.FIB
    end

    @testset "autoscale should fail on too few rows" begin
        A = rand(5,100)
        @test_throws AssertionError autoscale(A)
    end

    @testset "autoscale should fail on too few columns" begin
        A = rand(100,5)
        @test_throws AssertionError autoscale(A)
    end

    @testset "real image does not fail." begin
        @test sequence(SMALL; scales=[2]) !== nothing
    end

    @testset "small image does not fail." begin
        @test sequence(SMALL; scales=[2]) !== nothing
    end

    @testset "small image, multiple scales." begin
        @test sequence(SMALL; scales=[1,2,8]) !== nothing
    end

    @testset "large image does not fail." begin
        @test sequence(BIG; scales=[2], metrics=(KLD, WASS1D)) !== nothing
    end

    @testset "reorder small image, scale=1" begin
        Idx = shuffle(axes(SMALL,2))
        imgshuff = SMALL[:,Idx]
        seqres = sequence(imgshuff; scales=(1)) #all metrics
        ind = order(seqres)
        res = imgshuff[:, ind]
        # This way the test will also fail if the
        # reordering begins to work perfectly, i.e. loss ~ 0!
        # @test loss(SMALL,res) ≈ 0
        @test loss(SMALL,res) < 500
    end

    @testset "reorder small image, autoscale" begin
        Idx = shuffle(axes(SMALL,2))
        imgshuff = SMALL[:,Idx]
        seqres = sequence(imgshuff) #autoscale, all metrics
        ind = order(seqres)
        res = imgshuff[:, ind]
        # This way the test will also fail if the
        # reordering begins to work perfectly, i.e. loss ~ 0!
        # @test loss(SMALL,res) ≈ 0
        @test loss(SMALL,res) < 500
    end

    @testset "reorder med image, orig. portrait, scales=(1,2)" begin
        Idx = shuffle(axes(MED,2))
        imgshuff = MED[:,Idx]
        seqres = sequence(imgshuff; scales=(1,2), metrics=(KLD, L2))
        ind = order(seqres)
        res = imgshuff[:, ind]
        # This way the test will also fail if the
        # reordering begins to work perfectly, i.e. loss ~ 0!
        # @test loss(MED,res) ≈ 0
        @test loss(MED,res) < 700
    end

    @testset "reorder med image, orig. portrait, autoscale" begin
        Idx = shuffle(axes(MED,2))
        imgshuff = MED[:,Idx]
        seqres = sequence(imgshuff; metrics=(KLD, L2))
        ind = order(seqres)
        res = imgshuff[:, ind]
        # This way the test will also fail if the
        # reordering begins to work perfectly, i.e. loss ~ 0!
        # @test loss(MED,res) ≈ 0
        @test loss(MED,res) < 500
    end

    @testset "reorder med image, landscape, scales=(1,2)" begin
        MED = permutedims(MED)
        Idx = shuffle(axes(MED,2))
        imgshuff = MED[:,Idx]
        seqres = sequence(imgshuff; scales=(1,2), metrics=(KLD, L2))
        ind = order(seqres)
        res = imgshuff[:, ind]
        # This way the test will also fail if the
        # reordering begins to work perfectly, i.e. loss ~ 0!
        # @test loss(MED,res) ≈ 0
        @test loss(MED,res) < 1100
    end

    @testset "reorder med image, landscape, autoscale" begin
        MED = permutedims(MED)
        Idx = shuffle(axes(MED,2))
        imgshuff = MED[:,Idx]
        seqres = sequence(imgshuff; metrics=(KLD, L2))
        ind = order(seqres)
        res = imgshuff[:, ind]
        # This way the test will also fail if the
        # reordering begins to work perfectly, i.e. loss ~ 0!
        # @test loss(MED,res) ≈ 0
        @test loss(MED,res) < 1100
    end

    # This test fails on some combinations of random numbers
    @testset "reorder large image, orig. landscape, scales=(1,2,4)" begin
        Idx = shuffle(axes(BIG,2))
        imgshuff = BIG[:,Idx]
        seqres = sequence(imgshuff; scales=(1,2,4), metrics=(L2, KLD))
        ind = order(seqres)
        res = imgshuff[:, ind]
        @test loss(BIG,res) < 100
    end

    # This test fails on certain combinations of random numbers
    @testset "reorder large image, orig. landscape, autoscale" begin
        Idx = shuffle(axes(BIG,2))
        imgshuff = BIG[:,Idx]
        seqres = sequence(imgshuff; scales=(1,2,4), metrics=(L2, KLD))
        ind = order(seqres)
        res = imgshuff[:, ind]
        @test loss(BIG,res) < 100
    end

    # So far, this is the only case that systematically produces ~ 0 loss
    @testset "reorder large image, portrait, scales=(1,2,4)" begin
        BIG = permutedims(BIG)
        Idx = shuffle(axes(BIG,2))
        imgshuff = BIG[:,Idx]
        seqres = sequence(imgshuff; scales=(1,2,4), metrics=(L2, KLD))
        ind = order(seqres)
        res = imgshuff[:, ind]
        @test loss(BIG,res) ≈ 0
    end

    # So far, this is the only case that systematically produces ~ 0 loss
    @testset "reorder large image, portrait, autoscale" begin
        BIG = permutedims(BIG)
        Idx = shuffle(axes(BIG,2))
        imgshuff = BIG[:,Idx]
        seqres = sequence(imgshuff, metrics=(L2, KLD))
        ind = order(seqres)
        res = imgshuff[:, ind]
        @test loss(BIG,res) ≈ 0
    end

    @testset "PeriodicEuclidean small image scale=1" begin
        Idx = shuffle(axes(SMALL,2))
        imgshuff = SMALL[:,Idx]
        ρ = rand()
        # @show ρ
        P = vec(fill(ρ, size(imgshuff,1))) # period = 0.5 per dimension
        # @show P
        seqres = sequence(imgshuff, metrics=(PeriodicEuclidean(P),), scales=(1,)) #autoscale, all metrics
        ind = order(seqres)
        res = imgshuff[:, ind]
        # This way the test will also fail if the
        # reordering begins to work perfectly, i.e. loss ~ 0!
        # @test loss(SMALL,res) ≈ 0
        @test loss(SMALL,res) < 25
    end

    @testset "PeriodicEuclidean med image scale=1" begin
        Idx = shuffle(axes(MED,2))
        imgshuff = MED[:,Idx]
        @show size(imgshuff)
        ρ = 0.5
        @show ρ
        P = vec(fill(ρ, size(imgshuff,1))) # period = 0.5 per dimension
        @show P
        seqres = sequence(imgshuff, metrics=(PeriodicEuclidean(P),), scales=(1,)) #autoscale, all metrics
        ind = order(seqres)
        res = imgshuff[:, ind]
        # This way the test will also fail if the
        # reordering begins to work perfectly, i.e. loss ~ 0!
        # @test loss(MED,res) ≈ 0
        @test loss(MED,res) < 30
    end

    @testset "PeriodicEuclidean large image scale=1" begin
        Idx = shuffle(axes(BIG,2))
        imgshuff = BIG[:,Idx]
        @show size(imgshuff)
        ρ = 0.5
        @show ρ
        P = vec(fill(ρ, size(imgshuff,1))) # period = 0.5 per dimension
        seqres = sequence(imgshuff, metrics=(PeriodicEuclidean(P),), scales=(1,))
        ind = order(seqres)
        res = imgshuff[:, ind]
        # This way the test will also fail if the
        # reordering begins to work perfectly, i.e. loss ~ 0!
        # @test loss(BIG,res) ≈ 0
        @test loss(BIG,res) < 5
    end

    @testset "Accessor functions" begin
        seqres = sequence(SMALL; scales=(2,), metrics=(L2,))
        for f in (:elong, :order, :D, :mst)
            @testset "$f" begin
                @eval @test $f($seqres) !==nothing
            end
        end
    end

    @testset "show" begin
        seqres = sequence(SMALL; scales=(2,), metrics=(L2,))
        @test occursin("Sequencer Result", "$(seqres)")
    end

    @testset "prettyp of a vector" begin
        v = collect(1:20)
        r = "1,2,3...18,19,20"
        @test r == prettyp(v)
    end

    @testset "prettyp of a vector, more digits" begin
        v = collect(1:20)
        r = "1,2,3,4,5...16,17,18,19,20"
        @test r == prettyp(v, 5)
    end

    @testset "prettyp of a vector, small" begin
        v = collect(1:5)
        r = "1,2,3,4,5"
        @test r == prettyp(v)
    end
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

using Base.Test

using Base.LinAlg: BlasComplex, BlasFloat, BlasReal

n = 10

# Split n into 2 parts for tests needing two matrices
n1 = div(n, 2)
n2 = 2*n1

srand(1234321)

areal = randn(n,n)/2
aimg  = randn(n,n)/2
a2real = randn(n,n)/2
a2img  = randn(n,n)/2
breal = randn(n,2)/2
bimg  = randn(n,2)/2

@testset for eltya in (Float32, Float64, Complex64, Complex128)
    a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal)
    a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(a2real, a2img) : a2real)
    asym = a'+a                  # symmetric indefinite
    apd  = a'*a                 # symmetric positive-definite
    ε = εa = eps(abs(float(one(eltya))))

    @testset for eltyb in (Float32, Float64, Complex64, Complex128, Int)
        b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex(breal, bimg) : breal)
        εb = eps(abs(float(one(eltyb))))
        ε = max(εa,εb)

        α = rand(eltya)
        aα = fill(α,1,1)
        @test lqfact(α)[:L]*lqfact(α)[:Q] ≈ lqfact(aα)[:L]*lqfact(aα)[:Q]
        @test lq(α)[1]*lq(α)[2] ≈ lqfact(aα)[:L]*lqfact(aα)[:Q]
        @test abs(lqfact(α)[:Q][1,1]) ≈ one(eltya)
        tab = promote_type(eltya,eltyb)

        for i = 1:2
            let a = i == 1 ? a : view(a, 1:n - 1, 1:n - 1), b = i == 1 ? b : view(b, 1:n - 1), n = i == 1 ? n : n - 1
                lqa   = lqfact(a)
                l,q   = lqa[:L], lqa[:Q]
                qra   = qrfact(a)
                @testset "Basic ops" begin
                    @test size(lqa,1) == size(a,1)
                    @test size(lqa,3) == 1
                    @test size(lqa[:Q],3) == 1
                    @test Base.LinAlg.getq(lqa) == lqa[:Q]
                    @test_throws KeyError lqa[:Z]
                    @test full(lqa') ≈ a'
                    @test lqa * lqa' ≈ a * a'
                    @test lqa' * lqa ≈ a' * a
                    @test q*full(q, thin = false)' ≈ eye(eltya,n)
                    @test l*q ≈ a
                    @test full(lqa) ≈ a
                    @test full(copy(lqa)) ≈ a
                end
                @testset "Binary ops" begin
                    @test_approx_eq_eps a*(lqa\b) b 3000ε
                    @test_approx_eq_eps lqa*b qra[:Q]*qra[:R]*b 3000ε
                    @test_approx_eq_eps A_mul_Bc(eye(eltyb,size(q.factors,2)),q)*full(q, thin=false) eye(n) 5000ε
                    if eltya != Int
                        @test eye(eltyb,n)*q ≈ convert(AbstractMatrix{tab},q)
                    end
                    @test_approx_eq_eps q*b full(q, thin=false)*b 100ε
                    @test_approx_eq_eps q.'*b full(q, thin=false).'*b 100ε
                    @test_approx_eq_eps q'*b full(q, thin=false)'*b 100ε
                    @test_approx_eq_eps a*q a*full(q, thin=false) 100ε
                    @test_approx_eq_eps a*q.' a*full(q, thin=false).' 100ε
                    @test_approx_eq_eps a*q' a*full(q, thin=false)' 100ε
                    @test_throws DimensionMismatch q*b[1:n1 + 1]
                    @test_throws DimensionMismatch Ac_mul_B(q,ones(eltya,n+2,n+2))
                    @test_throws DimensionMismatch ones(eltyb,n+2,n+2)*q
                end
            end
        end

        @testset "Matmul with LQ factorizations" begin
            lqa = lqfact(a[:,1:n1])
            l,q = lqa[:L], lqa[:Q]
            @test full(q)*full(q)' ≈ eye(eltya,n1)
            @test (full(q,thin=false)'*full(q,thin=false))[1:n1,:] ≈ eye(eltya,n1,n)
            @test_throws DimensionMismatch A_mul_B!(eye(eltya,n+1),q)
            @test Ac_mul_B!(q,full(q)) ≈ eye(eltya,n1)
            @test_throws DimensionMismatch A_mul_Bc!(eye(eltya,n+1),q)
            @test_throws BoundsError size(q,-1)
        end
    end
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
struct FileFactor{T<:Unsigned} <: AbstractFactor{T}
    name::String
    length::Int64
    levels::AbstractVector{<:AbstractString}
    datapath::String
end

Base.length(factor::FileFactor{T}) where {T<:Unsigned} = factor.length

function slice(factor::FileFactor{T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:Unsigned}
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    sizeT = sizeof(T)
    iostream = open(factor.datapath, "r")
    seek(iostream, sizeT * (fromobs - 1))
    buffer = Vector{T}(undef, slicelength)
    map(Seq(Vector{T}, (buffer, iostream, fromobs, toobs, slicelength), nextdatachunk), SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}) do slice
        view(slice, 1:length(slice))
    end
end

function Base.map(factor::FileFactor, dataframe::AbstractDataFrame; permute::Bool = false)
    if permute
        PermuteFactor(factor, dataframe[getname(factor)])
    else
        dataframe[getname(factor)]
    end
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
iy = ix
iz = ix
incx = 1
incy = 1
incz = 1


text:
using StringDistances, Base.Test

@test_approx_eq_eps compare(Winkler(Jaro(), 0.1, 0.0), "martha", "marhta") 0.9611 1e-4
@test_approx_eq_eps compare(Winkler(Jaro(), 0.1, 0.0), "dwayne", "duane") 0.84 1e-4
@test_approx_eq_eps compare(Winkler(Jaro(), 0.1, 0.0), "dixon", "dicksonx") 0.81333 1e-4
@test_approx_eq_eps compare(Winkler(Jaro(), 0.1, 0.0), "william", "williams") 0.975 1e-4
@test_approx_eq_eps compare(Winkler(Jaro(), 0.1, 0.0), "", "foo") 0.0 1e-4
@test_approx_eq_eps compare(Winkler(Jaro(), 0.1, 0.0), "a", "a") 1.0 1e-4
@test_approx_eq_eps compare(Winkler(Jaro(), 0.1, 0.0), "abc", "xyz") 0.0 1e-4

strings = [
("martha", "marhta"),
("dwayne", "duane") ,
("dixon", "dicksonx"),
("william", "williams"),
("", "foo"),
("a", "a"),
("abc", "xyz"),
("abc", "ccc"),
("kitten", "sitting"),
("saturday", "sunday"),
("hi, my name is", "my name is"),
("alborgów", "amoniak"),
("cape sand recycling ", "edith ann graham"),
( "jellyifhs", "jellyfish"),
("ifhs", "fish"),
("leia", "leela"),
]
solutions = [0.03888889 0.16000000 0.18666667 0.02500000 1.00000000 0.00000000 1.00000000 0.44444444 0.25396825 0.22250000 0.16190476 0.43928571 0.49166667 0.04444444 0.16666667 0.17333333]
for i in 1:length(solutions)
	@test_approx_eq_eps compare(Winkler(Jaro(), 0.1, 0.0), strings[i]...) (1 - solutions[i]) 1e-4
end

@test_approx_eq_eps compare(Hamming(), "", "abc") 0.0 1e-4
@test_approx_eq_eps compare(Hamming(), "acc", "abc") 2/3 1e-4
@test_approx_eq_eps compare(Hamming(), "saturday", "sunday") 1/8 1e-4

@test_approx_eq_eps compare(QGram(1), "", "abc") 0.0 1e-4
@test_approx_eq_eps compare(QGram(1), "abc", "cba") 1.0 1e-4
@test_approx_eq_eps compare(QGram(1), "abc", "ccc") 1/3 1e-4

@test_approx_eq compare(Partial(RatcliffObershelp()), "New York Yankees",  "Yankees") 1.0
@test_approx_eq compare(Partial(RatcliffObershelp()), "New York Yankees",  "") 0.0
@test_approx_eq compare(Partial(RatcliffObershelp()),"mariners vs angels", "los angeles angels at seattle mariners") 0.444444444444

s = "HSINCHUANG"
@test_approx_eq compare(Partial(RatcliffObershelp()), s, "SINJHUAN") 0.875
@test_approx_eq compare(Partial(RatcliffObershelp()), s, "LSINJHUANG DISTRIC") 0.8
@test_approx_eq compare(Partial(RatcliffObershelp()), s, "SINJHUANG DISTRICT") 0.8
@test_approx_eq compare(Partial(RatcliffObershelp()), s,  "SINJHUANG") 0.8888888888888

@test_approx_eq compare(Partial(Hamming()), "New York Yankees",  "Yankees") 1
@test_approx_eq compare(Partial(Hamming()), "New York Yankees",  "") 1

@test_approx_eq compare(TokenSort(RatcliffObershelp()), "New York Mets vs Atlanta Braves", "Atlanta Braves vs New York Mets")  1.0
@test_approx_eq compare(TokenSet(RatcliffObershelp()),"mariners vs angels", "los angeles angels of anaheim at seattle mariners") 1.0 - 0.09090909090909094

@test_approx_eq compare(TokenSort(RatcliffObershelp()), "New York Mets vs Atlanta Braves", "")  0.0
@test_approx_eq compare(TokenSet(RatcliffObershelp()),"mariners vs angels", "") 0.0

@test_approx_eq compare(TokenMax(RatcliffObershelp()),"mariners vs angels", "") 0.0

#@test_approx_eq compare(TokenSort(RatcliffObershelp()), graphemeiterator("New York Mets vs Atlanta Braves"), graphemeiterator("Atlanta Braves vs New York Mets"))  1.0
#@test_approx_eq compare(TokenSet(RatcliffObershelp()),graphemeiterator("mariners vs angels"), graphemeiterator("los angeles angels of anaheim at seattle mariners")) 1.0 - 0.09090909090909094
#@test_approx_eq compare(TokenSort(RatcliffObershelp()), graphemeiterator("New York Mets vs Atlanta Braves"), graphemeiterator(""))  0.0
#@test_approx_eq compare(TokenSet(RatcliffObershelp()),graphemeiterator("mariners vs angels"), graphemeiterator("")) 0.0

@test_approx_eq compare(Winkler(Partial(Jaro())),"mariners vs angels", "los angeles angels at seattle mariners") 0.7378917378917379
@test_approx_eq compare(TokenSet(Partial(RatcliffObershelp())),"mariners vs angels", "los angeles angels at seattle mariners") 1.0