text:
using ZMQ, Mass2
context=Context()
sub=Socket(context, SUB)
ZMQ.connect(sub, "tcp://localhost:5555")
ZMQ.set_subscribe(sub, "0")
t0 = time()
function plot_hist(coadded_hist,n)
	fig = figure(1)
	fig[:clf]()
	ax = gca()
	ax[:plot](midpoints(coadded_hist), counts(coadded_hist),label="coadded $n ch")
	xlabel("energy (eV)")
	ylabel("counts per $(binsize(coadded_hist)) bin")
	title("using $n channels, $(counted(coadded_hist)) total counts")
	legend()
	show()
end

try
	while true
		msg = ZMQ.recv(sub)
		println("recieved message $(time()-t0) after starting")
		out=seek(convert(IOStream, msg),1)
		coadded_hist,n = Base.deserialize(out)
		plot_hist(coadded_hist,n)
	end
finally
	close(sub)
	close(context)
end


text:
module APCA
# Write your package code here.
using DataFrames
using Wavelets
using StatsBase
using ShiftedArrays

function runAPCA(oC::Vector, M::Int)
     #K. Chakrabarti et al.
     #Table V. An Algorithm to Produce the APCA
     #Algorithm Compute APCA(C,M )2
     #begin
     #1. if length(C) is not a power of two, pad it with zeros to make it so.
     #2. Perform the Haar Discrete Wavelet Transform on C.
     #3. Sort coefficients in order of decreasing normalized magnitude, truncate after M.
     #4. Reconstruct approximation (APCA representation) of C from retained coeffs.
     #5. If C was padded with zeros, truncate it to the original length.
     #6. Replace approximate segment mean values with exact mean values.
     #7. while the number of segments is greater than M
     #8. Merge the pair of segments that can be merged with least rise in error
     #9. endwhile
     #end
    # test values
    #M  = 3
    #oC = [7.,5.,5.,3.,3.,3.,4.,6.] # example in paper
    #oC = vcat(rand(40).+3,rand(40).+2,rand(40).-1)

    # 1. If length of oC is not a power of two, pad it with zeros
    newexp = ceil(log(size(oC)[1])/ log(2) )
    padsize = Int32(2^newexp - size(oC)[1])
    oCp = vcat(oC,fill(0.,padsize))
    fvlength = log2(size(oCp)[1]) # length power

    # 2. Perform the Haar Discrete Wavelet Transform on C.
    Haar = dwt(oCp,wavelet(WT.haar))

    #  Transform into weighting and structure in paper
    transHaar = []
    for i in [1:1:Int32(fvlength);]
        j = i-1
        addvec = -Haar[2^j+1:2^i] ./ sqrt(2^(fvlength - i + 1)) # get diff coeffs
        addvec2 = addvec ./ 2^((i-1)/2) # normalize
        append!(transHaar,addvec2)
    end
    #add final mean to normalized and non-normalized
    transHaar=append!(first(Haar,1)/sqrt(2^fvlength),transHaar)

    #3. Sort coefficients in order of decreasing normalized magnitude, truncate after M
    # go to float32 because of rounding
    if size(transHaar)[1] > M
        inds = sort(DataFrame(t=abs.(Float32.(transHaar)), i=1:length(transHaar)),:t,rev=true)[!,:i][1:M,:]
    else
        inds = sort(DataFrame(t=abs.(Float32.(transHaar)), i=1:length(transHaar)),:t,rev=true)[!,:i]
    end
    # get non-normalized coefs
    #4. Reconstruct approximation (APCA representation) of C from retained coeffs.
    Haar[Not(inds)] .= 0
    newrep =  idwt(Haar,wavelet(WT.haar))

    #5. If C was padded with zeros, truncate it to the original length.
    newrep = newrep[1:size(oC)[1]]

    #6. Replace approximate segment mean values with exact mean values.
    repeatlengths = vcat(1,rle(newrep)[2]) # get repeat lengths

    exactrep=[]
    for k in 2:size(repeatlengths)[1]
        s =sum(repeatlengths[1:k-1])
        e =repeatlengths[k] + s - 1
        um =mean(oCp[s:e])

        append!(exactrep,fill(um,repeatlengths[k]))
    end

    #7. while the number of segments is greater than M
    while size(rle(exactrep)[2])[1] > M
        z=vcat(rle(exactrep)[1] .- lag(rle(exactrep)[1]))
        replace!(z, missing =>-999.)
        z=abs.(convert(Vector{Float64}, z))
        l=findfirst(minimum(z) .== z)

        newreplengths = vcat(1,rle(exactrep)[2])
        exactrep
        s =sum(newreplengths[1:l-1])
        e =newreplengths[l]+newreplengths[l+1] + s - 1
        um =mean(exactrep[s:e])
        exactrep[s:e] .= um

    end

    return exactrep

  end
export runAPCA

end

ext: jl


text:
#=
Model: Heat (200 variables, 1 input)
=#
using Reachability, MAT, SparseArrays

heat(o::Pair{Symbol, <:Any}...) = heat(Options(Dict{Symbol, Any}(o)))

function heat(input_options::Options)
    # =====================
    # Problem specification
    # =====================
    file = matopen(@relpath "heat.mat")
    A = read(file, "A")

    # initial set
    # - x1-x300 are 0.0,
    # - the rest is in [0.002, 0.0015]
    X0 = Hyperrectangle([fill(0.6125, 2); zeros(198)],
                        [fill(0.0125, 2); zeros(198)])

    # input set
    B = sparse([67], [1], [1.0], size(A, 1), 1)
    U = BallInf([0.0], 0.5)

    # instantiate continuous LTI system
    S = InitialValueProblem(
        ConstrainedLinearControlContinuousSystem(A, B, nothing, U), X0)

    # property: x133 < 0.1
    property = SafeStatesProperty(HalfSpace(sparsevec([133], [1.0], 200), 0.1))

    # =======================
    # Problem default options
    # =======================
    partition = [(2*i-1:2*i) for i in 1:100] # 2D blocks

    if input_options[:mode] == "reach"
        problem_options = Options(:vars => [133],
                                  :partition => partition,
                                  :plot_vars => [0, 133])
    elseif input_options[:mode] == "check"
        problem_options = Options(:vars => [133],
                                  :partition => partition,
                                  :property => property)
    end

    return (S, merge(problem_options, input_options))
end

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
include("../envs/RLEnvironments.jl")
using .RLEnvironments
using Plots

function act(Q, s, ϵ)
    rand() <= ϵ && return rand(1:size(Q, 1))
    a = @view Q[:, s]
    return rand(findall(a .== maximum(a)))
end

function episode!(e, env, maxt, Q, ϵ)
    reset!(env)
    done = false
    t = 0
    while !done && t < maxt
        t += 1
        s = env.current
        a = act(Q, s, ϵ)
        _, r, done = step!(env, a)
        e[t] = (s, a, r)
    end
    return t
end

function run!(Q, env; epochs=20, maxt=100, ϵ=0.1, α=0.5, γ=0.99)
    rewards = zeros(epochs)
    steps = zeros(Int64, epochs)
    e = Vector{Tuple{Int64, Int64, Float64}}(undef, maxt)

    for i in 1:epochs
        t = episode!(e, env, maxt, Q, ϵ)
        rdisc = 0.0
        for step in Iterators.reverse(view(e, 1:t))
            rdisc = step[3] + γ * rdisc
            Q[step[2], step[1]] += α * (rdisc - Q[step[2], step[1]])

            rewards[i] += step[3]
        end
        steps[i] = t
        ϵ = max(ϵ * 0.995, 0.01)
    end
    return rewards, steps
end

env = SimpleRooms{Int64, Int64, Float64}()
Q = zeros(length(env.actionspace), length(env.observationspace))

r, s = run!(Q, env)
plt = plot([r s], labels=["Reward" "Steps"], color=[:red :blue], layout=(2, 1))
display(plt)
#savefig(plt, "mccontrol")


text:
# curly A{a,b...} in the context of *using* a parametric type
# declaring a parametric function or parametric type are separately
# considered in lintfunction and linttype, respectively.

# contracts for common collections / type parametrized types
# TODO: Can we be more specific here? What about detecting contract?
const CURLY_CONTRACTS = Dict{Symbol, Any}(
    :Array  => (Type, Integer),
    :Dict   => (Type, Type),
    :Matrix => (Type,),
    :Set    => (Type,),
    :Type   => (Type,),
    :Val    => (Any,),
    :Vector => (Type,))

function lintcurly(ex::Expr, ctx::LintContext)
    head = ex.args[1]
    if head == :Ptr && length(ex.args) == 2 && ex.args[2] == :Void
        return
    end
    contract = get(CURLY_CONTRACTS, head, nothing)
    for i = 2:length(ex.args)
        a = ex.args[i]
        if isexpr(a, :parameters) # only used for Traits.jl, AFAIK
            continue # grandfathered. We worry about linting this later
        elseif isexpr(a, :($))
            continue # grandfathered
        else
            t = guesstype(a, ctx)
            if !(t <: Type || t == Symbol || isbits(t) || t == Any)
                msg(ctx, :W441, a, "probably illegal use inside curly")
            elseif contract != nothing
                if i - 1 > length(contract)
                    msg(ctx, :W446, head, "too many type parameters")
                elseif !(t <: contract[i - 1] || t == Any)
                    msg(ctx, :W447, t, "can't be #$(i-1) type parameter for $head;" *
                        "it should be of type $(contract[i-1])")
                end
            end
            lintexpr(a, ctx)
        end
    end
end


text:
@testset "ZeroOne" begin

	@testset "Score Fitness" begin
		Random.seed!(SEED)

		#Dataset
		dset = BitsMNIST.Datasets.mnist()
		sx, sy = BitsMNIST.ZeroOne.sample(dset["train_x"], dset["train_y"], 0.01) #600 examples

		#Model
		model = Chain(BitDense(784, 800), BitDense(800, 2))

		#Error t0
		error0 = BitsMNIST.Statistics.error(model, sx, sy)

		#Genetic: it will have to be updated in the next major release
		fitness = BitsMNIST.ZeroOne.generate_score_fitness(sx, sy)
		tset = Genetic.TrainingSet(model, model.layers, fitness, mutationRate=0.05)
		Genetic.train!(tset, genNumber=5)

		#Error t1
		error1 = BitsMNIST.Statistics.error(model, sx, sy)

		@test error1 < error0
	end

	@testset "MCC Fitness" begin
		Random.seed!(SEED)

		#Dataset
		dset = BitsMNIST.Datasets.mnist()
		sx, sy = BitsMNIST.ZeroOne.sample(dset["train_x"], dset["train_y"], 0.01) #600 examples

		#Model
		model = Chain(BitDense(784, 800), BitDense(800, 2))

		#Error t0
		error0 = BitsMNIST.Statistics.error(model, sx, sy)

		#Genetic: it will have to be updated in the next major release
		fitness = BitsMNIST.ZeroOne.generate_mcc_fitness(sx, sy)
		tset = Genetic.TrainingSet(model, model.layers, fitness, mutationRate=0.05)
		Genetic.train!(tset, genNumber=5)

		#Error t1
		error1 = BitsMNIST.Statistics.error(model, sx, sy)

		@test error1 < error0
	end

end


text:
"""
    inverse_order!(s, var1, var2, p)

For solution `s`, this function reverses the order of the values between `var1` and `var2`. For example, the inverse order of [1.2, 3.4, 6.0, 8.9] is [8.9, 6.0, 3.4, 1.2].

The status of the solution (e.g., feasibility) is updated using the information in problem `p`.

# Example
See examples in tests.
"""
function inverse_order!(s::FixedLengthArray{T},
					    firstVariable::Int,
					    lastVariable::Int,
					    p::Problem) where {T<:Real}
	local numberOfElements::Int = lastVariable - firstVariable;
	local counter = 0
	for i in 0:2:numberOfElements
		swap_values!(s, firstVariable+counter, lastVariable-counter, p)
		counter += 1
	end
end

"""
    undo_inverse!(s, var1, var2, p)

This function undoes the last called of function `inverse_order`.

# Example
See examples in tests.
"""
function undo_inverse!(s::FixedLengthArray{T},
					   firstVariable::Int,
					   lastVariable::Int,
					   p::Problem) where {T<:Real}
	inverse_order!(s, firstVariable, lastVariable, p);
end

"""
    swap_values!(s, var1, var2, p)

For solution `s`, this function changes the values of variables `var1` and `var2`, so the `var1` now takes the value of `var2` and `var2` the value of `var1`.

The status of the solution (e.g., feasibility) is updated using the information in problem `p`.

# Example
See examples in tests.
"""
function swap_values!(s::FixedLengthArray{T},
					  variable1::Int,
					  variable2::Int,
					  p::Problem) where {T<:Real}
	local value1::T = get_solution(s, variable1)
	local value2::T = get_solution(s, variable2)
	add_solution_and_update_status!(s, variable2, value1, p);
	add_solution_and_update_status!(s, variable1, value2, p);
end

"""
    undo_swap!(s, var1, var2, p)

This function undoes the last called of function `swap_values`.

# Example
See examples in tests.
"""
function undo_swap!(s::FixedLengthArray{T},
					variable1::Int,
					variable2::Int,
					p::Problem) where {T<:Real}
	swap_values!(s, variable1, variable2, p); # doing the swap again undoes the swap
end

"""
    flip_value!(s, var, p)

For solution `s`, this function puts `var` to 0 if it has value 1 and viceversa. If the value of the variable is different than 1 or 0, then nothing happens.

The status of the solution (e.g., feasibility) is updated using the information in problem `p`.

# Example
See examples in tests.
"""
function flip_value!(s::FixedLengthArray{T},
					 variable::Int,
					 p::Problem) where {T<:Real}
	local value::T = get_solution(s, variable);
	if value == zero(T)
		add_solution_and_update_status!(s, variable, one(T), p);
	elseif value == one(T)
		add_solution_and_update_status!(s, variable, zero(T), p);
	end
end

"""
    undo_flip!(s, var, p)

This function undoes the last called of function `flip_value`.

# Example
See examples in tests.
"""
function undo_flip!(s::FixedLengthArray{T},
					variable::Int,
					p::Problem) where {T<:Real}
	flip_value!(s, variable, p); # if we flip again, we undo the change
end

"""
    mirror_value!(s, var, p)

If variable `var` in solution `s` has a value `x[var]=z`, and the variable is bounded to the interval (it may be discrete) `[lb, ub]`, then the variable is assigned value `ub+lb-z`.

For example, if `x[var]=4` and `[lb, ub] = [1,5]`, then, after calling the function, `x[var] = 2`. Notice that if you draw the bounds in the real line: 1-2-3-4-5, value 2 is the mirror of 4. Meaning that if you fold the interval by its middle point (which in the example is 3), value 2 and 4 would touch.

If the value of the variable is not within its bounds, then nothing happens.

The status of the solution (e.g., feasibility) is updated using the information in problem `p`.

# Example
See examples in tests.
"""
function mirror_value!(s::FixedLengthArray{T},
					   variable::Int,
					   p::Problem) where {T<:Real}
	local value::T = get_solution(s, variable);
	if is_value_within_the_domain(p, variable, value)
		local sumBounds::T = T(get_ub_variable(p, variable) + get_lb_variable(p, variable));
		local newValue::T =  sumBounds - value;
		add_solution_and_update_status!(s, variable, newValue, p);
	end
end

"""
    undo_mirror!(s, var, p)

This function undoes the last called of function `mirror_value`.

# Example
See examples in tests.
"""
function undo_mirror!(s::FixedLengthArray{T},
					  variable::Int,
					  p::Problem) where {T<:Real}
	mirror_value!(s, variable, p); # if we mirror again, we undo the change
end


text:
abstract type Field end
abstract type Vortex end
abstract type CoreShape end

"""
    Torus <: Field

Specify toroidal boundary conditions: doubly periodic in 2D.

# Arguments
- `ψ::Array{Complex{Float64},2}`: wavefunction.
- `x::Vector{Float64}`: `x` coordinates.
- `y::Vector{Float64}`: `y` coordinates.

See also: [`Sphere`](@def), [`Field`](@ref)
"""
@with_kw mutable struct Torus <: Field
    ψ::Array{Complex{Float64},2}
    x::Vector{Float64}
    y::Vector{Float64}
end

"""
    Sphere <: Field

Specify spherical boundary conditions in 2D.

# Arguments
- `ψ::Array{Complex{Float64},2}`: wavefunction.
- `x::Vector{Float64}`: `x` coordinates.
- `y::Vector{Float64}`: `y` coordinates.

See also: [`Torus`](@def), [`Field`](@ref)
"""
@with_kw mutable struct Sphere <: Field
    ψ::Array{Complex{Float64},2}
    x::Vector{Float64}
    y::Vector{Float64}
end

@with_kw mutable struct PointVortex <: Vortex
    xv::Float64
    yv::Float64
    qv::Int64
end

@with_kw mutable struct Ansatz <: CoreShape
    f::Interpolations.GriddedInterpolation{Float64,1,Float64,Gridded{Linear},Tuple{Array{Float64,1}}}
    ξ::Float64
    Λ::Float64
end

@with_kw mutable struct Exact <: CoreShape
    f::Interpolations.GriddedInterpolation{Float64,1,Float64,Gridded{Linear},Tuple{Array{Float64,1}}}
    ξ::Float64
end

@with_kw mutable struct ScalarVortex{T<:CoreShape} <: Vortex
    core::T
    vort::PointVortex
end


text:
using Plots
# physical coordinates
inner_radius = 6378100 # [m]
outer_radius = 6378100 + 100000 # [m] adding 100km for the atmosphere [about 1.6%] (its 3km for the ocean)
r1 = (inner_radius / outer_radius)
#

nr = 5
nh = 10
ξ¹ = reshape(collect(range(-1,1, length = nh)), (1,nh))
ξ² = reshape(collect(range(-1,1, length = nr)) , (nr,1))

r1 = 0.5 # inner radius
r2 = 1.0 # outer radius
npolygon = 5 # polygon approximation
φ = 2π/(2 * npolygon)
x¹ = @. ( (r2-r1)*sin(φ) * (ξ² + 1)/2 + r1 * sin(φ) ) * ξ¹
x² = @. (r2-r1)*cos(φ) * (ξ²+1)/2 + 0 * ξ¹ + r1*cos(φ)
p1 = scatter(x¹[:], x²[:], xlims = (-r2*1.1,r2*1.1), ylims = (-r2*1.1,r2*1.1), legend = false,aspectratio = 1)

##
for j in 1:1:(npolygon-1)
    θ = 2π/npolygon * j
    xrot = @. cos(θ) * x¹ - sin(θ) * x²
    yrot = @. sin(θ) * x¹ + cos(θ) * x²
    p1 = scatter!(xrot[:], yrot[:], xlims = (-r2*1.1,r2*1.1), ylims = (-r2*1.1,r2*1.1))
    display(p1)
    sleep(0.1)
end

plot(p1)

##
j=npolygon
θ = 2π/npolygon * j
xrot = @. cos(θ) * x¹ - sin(θ) * x²
yrot = @. sin(θ) * x¹ + cos(θ) * x²
p1 = scatter(xrot[:], yrot[:], xlims = (-r2*1.1,r2*1.1), ylims = (-r2*1.1,r2*1.1), aspectratio = 1, legend = false)
display(p1)
sleep(0.1)
scatter!(x¹[:], x²[:], xlims = (-r2*1.1,r2*1.1), ylims = (-r2*1.1,r2*1.1), legend = false,aspectratio = 1)

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
using DFTK
using ForwardDiff
using Test
include("testcases.jl")

@testset "Force derivatives using ForwardDiff" begin
    function compute_force(ε1, ε2; metal=false)
        T = promote_type(typeof(ε1), typeof(ε2))
        pos = [[1.01, 1.02, 1.03] / 8, -ones(3) / 8 + ε1 * [1., 0, 0] + ε2 * [0, 1., 0]]
        if metal
            # Silicon reduced HF is metallic
            model = model_DFT(Matrix{T}(silicon.lattice), silicon.atoms, pos, [];
                              temperature=1e-3)
        else
            model = model_LDA(Matrix{T}(silicon.lattice), silicon.atoms, pos)
        end
        basis = PlaneWaveBasis(model; Ecut=5, kgrid=[2, 2, 2], kshift=[0, 0, 0])

        is_converged = DFTK.ScfConvergenceDensity(1e-10)
        scfres = self_consistent_field(basis; is_converged,
                                       response=ResponseOptions(verbose=true))
        compute_forces_cart(scfres)
    end

    F = compute_force(0.0, 0.0)
    derivative_ε1_fd = let ε1 = 1e-5
        (compute_force(ε1, 0.0) - F) / ε1
    end
    derivative_ε1 = ForwardDiff.derivative(ε1 -> compute_force(ε1, 0.0), 0.0)
    @test norm(derivative_ε1 - derivative_ε1_fd) < 1e-4

    derivative_ε2_fd = let ε2 = 1e-5
        (compute_force(0.0, ε2) - F) / ε2
    end
    derivative_ε2 = ForwardDiff.derivative(ε2 -> compute_force(0.0, ε2), 0.0)
    @test norm(derivative_ε2 - derivative_ε2_fd) < 1e-4

    @testset "Multiple partials" begin
        grad = ForwardDiff.gradient(v -> compute_force(v...)[1][1], [0.0, 0.0])
        @test abs(grad[1] - derivative_ε1[1][1]) < 1e-4
        @test abs(grad[2] - derivative_ε2[1][1]) < 1e-4

        jac = ForwardDiff.jacobian(v -> compute_force(v...)[1], [0.0, 0.0])
        @test norm(grad - jac[1, :]) < 1e-10
    end

    @testset "Derivative for metals" begin
        F = compute_force(0.0, 0.0; metal=true)
        derivative_ε1_fd = let ε1 = 1e-5
            (compute_force(ε1, 0.0; metal=true) - F) / ε1
        end
        derivative_ε1 = ForwardDiff.derivative(ε1 -> compute_force(ε1, 0.0; metal=true), 0.0)
        @test norm(derivative_ε1 - derivative_ε1_fd) < 1e-4
    end
end

@testset "Force pseudo-sensitivity using ForwardDiff" begin
    # TODO Maybe later change to a band-energy test once eigenvalue derivatives work.
    function compute_force(ε::T) where {T}
        psp = load_psp("hgh/lda/si-q4")
        rloc = convert(T, psp.rloc)

        pspmod = PspHgh(psp.Zion, rloc,
                        psp.cloc, psp.rp .+ [0, ε], psp.h;
                        psp.identifier, psp.description)
        elem  = ElementPsp(silicon.atnum, psp=pspmod)
        atoms = [elem, elem]
        positions = [[1.01, 1.02, 1.03] / 8, -ones(3) / 8]
        model = model_LDA(Matrix{T}(silicon.lattice), atoms, positions)
        basis = PlaneWaveBasis(model; Ecut=5, kgrid=[2, 2, 2], kshift=[0, 0, 0])

        is_converged = DFTK.ScfConvergenceDensity(1e-10)
        scfres = self_consistent_field(basis; is_converged,
                                       response=ResponseOptions(verbose=true))
        compute_forces_cart(scfres)
    end

    derivative_fd = let ε = 1e-5
        (compute_force(ε) - compute_force(-ε)) / 2ε
    end
    derivative_ε = ForwardDiff.derivative(compute_force, 0.0)
    @test norm(derivative_fd - derivative_ε) < 1e-4
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: workmailmessageflow

using AWS.Compat
using AWS.UUIDs
"""
    GetRawMessageContent()

Retrieves the raw content of an in-transit email message, in MIME format.

# Required Parameters
- `messageId`: The identifier of the email message to retrieve.

"""

get_raw_message_content(messageId; aws_config::AWSConfig=global_aws_config()) = workmailmessageflow("GET", "/messages/$(messageId)"; aws_config=aws_config)
get_raw_message_content(messageId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = workmailmessageflow("GET", "/messages/$(messageId)", args; aws_config=aws_config)


text: include(joinpath(@__DIR__, "example_data.jl"))


text:
function reof(
    array::AbstractArray{<:Any,3};
    nmodes = nothing,
    maxrotations::Int = 50,
    nsims::Int = 10,
)

    eof = EmpiricalOrthogonalFunction(array)

    if isnothing(nmodes)
        centered = eof.dataset[:, eof.valididx]
        eigenvals = eigenvalues(eof)
        sigmodes = montecarlo_significance(centered, eigenvals; nsims = nsims)
        # nmodes = max(sigmodes...)
        nmodes = length(sigmodes)
    else
        sigmodes = 1:nmodes
    end

    @info "Found $nmodes significant modes"

    eof_slice = index_eof(eof,sigmodes)

    reof = orthorotation(eof_slice)

    temporal_modes = convert(Array{Float32}, pcs(reof))
    spatial_modes = reshape(eofs(reof), (size(array)[1:2]..., nmodes))
    center = reshape(eof.center, size(array)[1:2])

    varfrac = variancefraction(reof)

    return center, spatial_modes, temporal_modes, varfrac
end

function montecarlo_significance(
    centered::AbstractArray{<:Any,2},
    eigenvals::AbstractArray{<:Any,1};
    nsims::Int = 10,
)
    nt = length(eigenvals)

    norm = nt - 1
    mc_lamb = zeros(Float64, nt, nsims)

    obstemp = zeros(size(centered))
    obstemp[:] = centered[:]

    nobs = size(obstemp,2)[1]

    @simd for i = 1:nsims
        # ----- Randomize the observation (space-wise randomization)
        # @threads for j = 1:size(obstemp, 1)
        #     seed = now().instant.periods.value
        #     rng = MersenneTwister(seed)
        #     shuffle!(rng, @view obstemp[i, :])
        # end
        seed = now().instant.periods.value
        rng = MersenneTwister(seed)
        rndidx = mapslices(sortperm, rand(rng, size(obstemp)...), dims=2)
        # obstemp =  obstemp[:,rndidx]

        # run eof to get eigen vales of randomized (in space) observations
        rnd_eof = EmpiricalOrthogonalFunction(@view obstemp[rndidx])

        mc_lamb[:, i] = eigenvalues(rnd_eof)
    end

    mean_mc_lamb = Statistics.mean(mc_lamb, dims = 2)
    std_mc_lamb = Statistics.std(mc_lamb, dims = 2)

    # eigendist = Normal.(mean_mc_lamb[:,1], std_mc_lamb[:,1])
    # mc_ci = quantile.(eigendist,α)
    # eigengt = mc_ci .< eigenvals

    eigengt = (eigenvals .> mean_mc_lamb) .& (eigenvals / sum(eigenvals) .> 1e-6)

    sig_modes = findall(eigengt[:, 1])

    return sig_modes
end

function index_eof(eof::EmpiricalOrthogonalFunction, I)
    eof_ = eof.eofs[:,I]
    pcs_ = eof.pcs[:,I]
    eigenvalues_ = eof.eigenvals[I]

    return EmpiricalOrthogonalFunction(
        eof.dataset,
        eof.center,
        eof.valididx,
        eof_,
        eigenvalues_,
        pcs_,
    )
end

function filter_dates(
    df::DataFrame,
    dates::Vector{DateTime};
    timecol::AbstractString = "Date",
)
    # TODO: check to make sure that dates range is within df[:,"datetime"] range
    # produces inconsistent results if not ... not sure why...
    mask = sum(map(x -> x .== df[:, timecol], dates), dims = 1)[1]
    indices = findall(!iszero, mask)

    df[indices, :]
end

function fitmodes(
    temporal_modes::AbstractArray{<:Real,2},
    dfs::Vector{DataFrame};
    datacol::AbstractString = "H",
    removeoutliers::Bool = true,
)

    nmodes = size(temporal_modes, 2)
    cormat = zeros(length(dfs), nmodes)
    outcoeffs = zeros(3, nmodes)

    @inbounds for i = 1:nmodes
        models = Polynomial[]

        @inbounds for (j, df) in enumerate(dfs)
            x = df[!, datacol]
            if removeoutliers
                x_z = StatsBase.zscore(x)
                idx = findall(x -> abs(x) <= 3.0, x_z)
            else
                idx = findall(x -> !ismissing(x), x_z)
            end

            @views cormat[j, i] = cor(x[idx], temporal_modes[idx, i])

            @views model = Polynomials.fit(x[idx], temporal_modes[idx, i], 2)

            push!(models, model)
        end

        m, best_idx = findmax(abs.(cormat[:, i]))
        best_model = models[best_idx]

        outcoeffs[:, i] = coeffs(best_model)
    end

    return outcoeffs, cormat
end

max_stars_repo_name: SERVIR/fier-cli

max_issues_repo_name: SERVIR/fier-cli

max_forks_repo_name: SERVIR/fier-cli


text:
@testset "Paired read datastores" begin
    function get_fastq_seqs(file)
        seqs = map(open(FASTQ.Reader, file) do rdr
            collect(rdr)
        end) do rec
            FASTQ.sequence(LongDNASeq, rec)
        end
        for seq in seqs
            if length(seq) > 300
                resize!(seq, 300)
            end
        end
        return seqs
    end

    function check_round_trip(R1, R2)
        r1_seqs = get_fastq_seqs(R1)
        r2_seqs = get_fastq_seqs(R2)

        fqa = open(FASTQ.Reader, R1)
        fqb = open(FASTQ.Reader, R2)

        ds = PairedReads(fqa, fqb,
                         "ecoli-pe.prds", "ecoli-pe",
                         UInt64(250), UInt64(300), UInt64(0), FwRv)

        ds_seqs = collect(ds)

        return ds_seqs[1:2:end] == r1_seqs && ds_seqs[2:2:end] == r2_seqs
    end

    function check_show(ds, msg)
        buf = IOBuffer()
        show(buf, ds)
        return String(take!(buf)) == msg
    end

    @test check_round_trip("ecoli_tester_R1.fastq", "ecoli_tester_R2.fastq")

    ds = open(PairedReads, "ecoli-pe.prds")
    @test ReadDatastores.name(ds) == "ecoli-pe"
    @test ReadDatastores.maxseqlen(ds) == 300
    @test ReadDatastores.orientation(ds) == FwRv
    @test check_show(ds, "Paired Read Datastore 'ecoli-pe': 20 reads")
    @test firstindex(ds) == 1
    @test lastindex(ds) == 20
    @test Base.IteratorSize(ds) == Base.HasLength()
    @test Base.IteratorEltype(ds) == Base.HasEltype()
    @test Base.eltype(ds) == LongSequence{DNAAlphabet{4}}
    @test_throws BoundsError ds[100]
    @test_throws BoundsError buffer(ds)[100]
    @test_throws BoundsError load_sequence!(ds, 100, dna"")
    @test collect(ds) == collect(buffer(ds)) == open(PairedReads, "ecoli-pe.prds") do ds
        collect(ds)
    end

end


text:
include("../src/RichVehicleRoutingProblem.jl")
const RVRP = RichVehicleRoutingProblem

using Test

include("unit_tests/unit_tests.jl")

unit_tests()


text:
include("lib_for_practic_10.jl")

function count_markers_in_labirint!(robot::CoordRobot)
    coord = (0,0)
    function recurs_for_labirint()
        coord = get_coord(get_coord(robot))
        if !(coord in set_coords)
            push!(set_coords, coord)
            if ismarker(robot)
                num_markers += 1
            end
            for side in (Nord, West, Sud, Ost)
                if try_move!(robot, side)
                    recurs_for_labirint()
                    move!(robot, inverse(side))
                end
            end
        end
    end

    set_coords = Set{NTuple{2,Int}}()
    num_markers = 0
    recurs_for_labirint()
    return num_markers
end

count_markers_in_labirint!(r)

max_stars_repo_name: AntonDudDm/mirea-progs

max_issues_repo_name: AntonDudDm/mirea-progs

max_forks_repo_name: AntonDudDm/mirea-progs


text:
function par2lpc(state::SpectralParamState{PartialAutoCorrelation};
                 kargs...)
    @assert !has_loggain(state)
    pardef = paramdef(state)
    data = rawdata(state)
    data = SPTK.par2lpc(data; kargs...)
    newdef = LinearPredictionCoef(param_order(pardef))
    SpectralParamState(newdef, data,
                       has_loggain(state), gain_normalized(state))
end


text:
# structure for iterating over a csv file
# no automatic type inference is done, but types are allowed to be passed
# for as many columns as desired; `CSV.detect(row, i)` can also be used to
# use the same inference logic used in `CSV.File` for determing a cell's typed value
struct Rows{transpose, O, IO}
    name::String
    names::Vector{Symbol} # only includes "select"ed columns
    types::Vector{Type} # only includes "select"ed columns
    columnmap::Vector{Int} # maps "select"ed column index to actual file column index
    typecodes::Vector{TypeCode} # includes *all* columns (whether or not selected)
    cols::Int64
    e::UInt8
    buf::IO
    datapos::Int64
    len::Int
    limit::Int64
    options::O # Parsers.Options
    coloptions::Union{Nothing, Vector{Parsers.Options}}
    positions::Vector{Int64}
    reusebuffer::Bool
    tapes::Vector{Vector{UInt64}}
    intsentinels::Vector{Int64}
    lookup::Dict{Symbol, Int}
end

function Base.show(io::IO, r::Rows)
    println(io, "CSV.Rows(\"$(r.name)\"):")
    println(io, "Size: $(r.cols)")
    show(io, Tables.schema(r))
end

"""
    CSV.Rows(source; kwargs...) => CSV.Rows

Read a csv input (a filename given as a String or FilePaths.jl type, or any other IO source), returning a `CSV.Rows` object.

While similar to [`CSV.File`](@ref), `CSV.Rows` provides a slightly different interface, the tradeoffs including:
  * Very minimal memory footprint; while iterating, only the current row values are buffered
  * Only provides row access via iteration; to access columns, one can stream the rows into a table type
  * Performs no type inference; each column/cell is essentially treated as `Union{String, Missing}`, users can utilize the performant `Parsers.parse(T, str)` to convert values to a  more specific type if needed

Opens the file and uses passed arguments to detect the number of columns, ***but not*** column types.
The returned `CSV.Rows` object supports the [Tables.jl](https://github.com/JuliaData/Tables.jl) interface
and can iterate rows. Each row object supports `propertynames`, `getproperty`, and `getindex` to access individual row values.
Note that duplicate column names will be detected and adjusted to ensure uniqueness (duplicate column name `a` will become `a_1`).
For example, one could iterate over a csv file with column names `a`, `b`, and `c` by doing:

```julia
for row in CSV.Rows(file)
    println("a=\$(row.a), b=\$(row.b), c=\$(row.c)")
end
```

Supported keyword arguments include:
* File layout options:
  * `header=1`: the `header` argument can be an `Int`, indicating the row to parse for column names; or a `Range`, indicating a span of rows to be concatenated together as column names; or an entire `Vector{Symbol}` or `Vector{String}` to use as column names; if a file doesn't have column names, either provide them as a `Vector`, or set `header=0` or `header=false` and column names will be auto-generated (`Column1`, `Column2`, etc.)
  * `normalizenames=false`: whether column names should be "normalized" into valid Julia identifier symbols; useful when iterating rows and accessing column values of a row via `getproperty` (e.g. `row.col1`)
  * `datarow`: an `Int` argument to specify the row where the data starts in the csv file; by default, the next row after the `header` row is used. If `header=0`, then the 1st row is assumed to be the start of data
  * `skipto::Int`: similar to `datarow`, specifies the number of rows to skip before starting to read data
  * `limit`: an `Int` to indicate a limited number of rows to parse in a csv file; use in combination with `skipto` to read a specific, contiguous chunk within a file
  * `transpose::Bool`: read a csv file "transposed", i.e. each column is parsed as a row
  * `comment`: rows that begin with this `String` will be skipped while parsing
  * `use_mmap::Bool=!Sys.iswindows()`: whether the file should be mmapped for reading, which in some cases can be faster
  * `ignoreemptylines::Bool=false`: whether empty rows/lines in a file should be ignored (if `false`, each column will be assigned `missing` for that empty row)
* Parsing options:
  * `missingstrings`, `missingstring`: either a `String`, or `Vector{String}` to use as sentinel values that will be parsed as `missing`; by default, only an empty field (two consecutive delimiters) is considered `missing`
  * `delim=','`: a `Char` or `String` that indicates how columns are delimited in a file; if no argument is provided, parsing will try to detect the most consistent delimiter on the first 10 rows of the file
  * `ignorerepeated::Bool=false`: whether repeated (consecutive) delimiters should be ignored while parsing; useful for fixed-width files with delimiter padding between cells
  * `quotechar='"'`, `openquotechar`, `closequotechar`: a `Char` (or different start and end characters) that indicate a quoted field which may contain textual delimiters or newline characters
  * `escapechar='"'`: the `Char` used to escape quote characters in a quoted field
  * `strict::Bool=false`: whether invalid values should throw a parsing error or be replaced with `missing`
  * `silencewarnings::Bool=false`: if `strict=false`, whether warnings should be silenced
* Iteration options:
  * `reusebuffer=false`: while iterating, whether a single row buffer should be allocated and reused on each iteration; only use if each row will be iterated once and not re-used (e.g. it's not safe to use this option if doing `collect(CSV.Rows(file))`)
"""
function Rows(source;
    # file options
    # header can be a row number, range of rows, or actual string vector
    header::Union{Integer, Vector{Symbol}, Vector{String}, AbstractVector{<:Integer}}=1,
    normalizenames::Bool=false,
    # by default, data starts immediately after header or start of file
    datarow::Integer=-1,
    skipto::Union{Nothing, Integer}=nothing,
    footerskip::Integer=0,
    limit::Integer=typemax(Int64),
    transpose::Bool=false,
    comment::Union{String, Nothing}=nothing,
    use_mmap::Bool=!Sys.iswindows(),
    ignoreemptylines::Bool=false,
    select=nothing,
    drop=nothing,
    # parsing options
    missingstrings=String[],
    missingstring="",
    delim::Union{Nothing, Char, String}=nothing,
    ignorerepeated::Bool=false,
    quotechar::Union{UInt8, Char}='"',
    openquotechar::Union{UInt8, Char, Nothing}=nothing,
    closequotechar::Union{UInt8, Char, Nothing}=nothing,
    escapechar::Union{UInt8, Char}='"',
    dateformat::Union{String, Dates.DateFormat, Nothing}=nothing,
    dateformats::Union{AbstractDict, Nothing}=nothing,
    decimal::Union{UInt8, Char}=UInt8('.'),
    truestrings::Union{Vector{String}, Nothing}=["true", "True", "TRUE"],
    falsestrings::Union{Vector{String}, Nothing}=["false", "False", "FALSE"],
    # type options
    type=nothing,
    types=nothing,
    typemap::Dict=Dict{TypeCode, TypeCode}(),
    categorical::Union{Bool, Real}=false,
    pool::Union{Bool, Real}=0.1,
    strict::Bool=false,
    silencewarnings::Bool=false,
    debug::Bool=false,
    parsingdebug::Bool=false,
    reusebuffer::Bool=false,
    kw...)

    h = Header(source, header, normalizenames, datarow, skipto, footerskip, limit, transpose, comment, use_mmap, ignoreemptylines, false, select, drop, missingstrings, missingstring, delim, ignorerepeated, quotechar, openquotechar, closequotechar, escapechar, dateformat, dateformats, decimal, truestrings, falsestrings, type, types, typemap, categorical, pool, strict, silencewarnings, debug, parsingdebug, true)
    tapes = [Vector{UInt64}(undef, usermissing(h.typecodes[i]) ? 0 : 1) for i = 1:h.cols]
    types = Type[gettype(T) for T in h.typecodes]
    columnmap = [i for i = 1:h.cols]
    deleteat!(h.names, h.todrop)
    deleteat!(types, h.todrop)
    deleteat!(columnmap, h.todrop)
    lookup = Dict(nm=>i for (i, nm) in enumerate(h.names))
    return Rows{transpose, typeof(h.options), typeof(h.buf)}(
        h.name,
        h.names,
        types,
        columnmap,
        h.typecodes,
        h.cols,
        h.e,
        h.buf,
        h.datapos,
        h.len,
        limit,
        h.options,
        h.coloptions,
        h.positions,
        reusebuffer,
        tapes,
        fill(INT_SENTINEL, h.cols),
        lookup
    )
end

Tables.rowaccess(::Type{<:Rows}) = true
Tables.rows(r::Rows) = r
Tables.schema(r::Rows) = Tables.Schema(r.names, r.types)
Base.eltype(r::Rows) = Row2
Base.IteratorSize(::Type{<:Rows}) = Base.SizeUnknown()

const EMPTY_TYPEMAP = Dict{TypeCode, TypeCode}()
const EMPTY_POSLENS = Vector{Vector{UInt64}}()
const EMPTY_REFS = Vector{Dict{String, UInt64}}()
const EMPTY_LASTREFS = UInt64[]

@inline function Base.iterate(r::Rows{transpose}, (pos, len, row)=(r.datapos, r.len, 1)) where {transpose}
    (pos > len || row > r.limit) && return nothing
    pos > len && return nothing
    tapes = r.reusebuffer ? r.tapes : [Vector{UInt64}(undef, usermissing(r.typecodes[i]) ? 0 : 1) for i = 1:r.cols]
    pos = parserow(1, Val(transpose), r.cols, EMPTY_TYPEMAP, tapes, EMPTY_POSLENS, r.buf, pos, len, r.limit, r.positions, 0.0, EMPTY_REFS, EMPTY_LASTREFS, 0, r.typecodes, r.intsentinels, false, r.options, r.coloptions)
    intsentinels = r.reusebuffer ? r.intsentinels : copy(r.intsentinels)
    return Row2(r.names, r.types, r.columnmap, r.typecodes, r.lookup, tapes, r.buf, r.e, r.options, r.coloptions, intsentinels), (pos, len, row + 1)
end

struct Row2{O} <: Tables.AbstractRow
    names::Vector{Symbol}
    types::Vector{Type}
    columnmap::Vector{Int}
    typecodes::Vector{TypeCode}
    lookup::Dict{Symbol, Int}
    tapes::Vector{Vector{UInt64}}
    buf::Vector{UInt8}
    e::UInt8
    options::O
    coloptions::Union{Nothing, Vector{Parsers.Options}}
    intsentinels::Vector{Int64}
end

getnames(r::Row2) = getfield(r, :names)
gettypes(r::Row2) = getfield(r, :types)
getcolumnmap(r::Row2) = getfield(r, :columnmap)
gettypecodes(r::Row2) = getfield(r, :typecodes)
getlookup(r::Row2) = getfield(r, :lookup)
gettapes(r::Row2) = getfield(r, :tapes)
getbuf(r::Row2) = getfield(r, :buf)
gete(r::Row2) = getfield(r, :e)
getoptions(r::Row2) = getfield(r, :options)
getcoloptions(r::Row2) = getfield(r, :coloptions)
getintsentinels(r::Row2) = getfield(r, :intsentinels)

Tables.columnnames(r::Row2) = getnames(r)

Base.checkbounds(r::Row2, i) = 0 < i < length(r)

Tables.getcolumn(r::Row2, nm::Symbol) = Tables.getcolumn(r, getlookup(r)[nm])
Tables.getcolumn(r::Row2, i::Int) = Tables.getcolumn(r, gettypes(r)[i], i, getnames(r)[i])

Base.@propagate_inbounds function Tables.getcolumn(r::Row2, ::Type{Missing}, i::Int, nm::Symbol)
    @boundscheck checkbounds(r, i)
    return missing
end

Base.@propagate_inbounds function Tables.getcolumn(r::Row2, ::Type{T}, i::Int, nm::Symbol) where {T}
    @boundscheck checkbounds(r, i)
    j = getcolumnmap(r)[i]
    @inbounds x = reinterp_func(T)(gettapes(r)[j][1])
    return x
end

Base.@propagate_inbounds function Tables.getcolumn(r::Row2, ::Type{Union{Missing, T}}, i::Int, nm::Symbol) where {T}
    @boundscheck checkbounds(r, i)
    j = getcolumnmap(r)[i]
    @inbounds x = gettapes(r)[j][1]
    return ifelse(x === sentinelvalue(T), missing, reinterp_func(T)(x))
end

Base.@propagate_inbounds function Tables.getcolumn(r::Row2, ::Type{Union{Missing, Int64}}, i::Int, nm::Symbol)
    @boundscheck checkbounds(r, i)
    j = getcolumnmap(r)[i]
    @inbounds x = reinterp_func(Int64)(gettapes(r)[j][1])
    return ifelse(x === getintsentinels(r)[j], missing, x)
end

Base.@propagate_inbounds function Tables.getcolumn(r::Row2, ::Type{String}, i::Int, nm::Symbol)
    @boundscheck checkbounds(r, i)
    j = getcolumnmap(r)[i]
    @inbounds offlen = gettapes(r)[j][1]
    s = PointerString(pointer(getbuf(r), getpos(offlen)), getlen(offlen))
    return escapedvalue(offlen) ? unescape(s, gete(r)) : String(s)
end

Base.@propagate_inbounds function Tables.getcolumn(r::Row2, ::Type{Union{Missing, String}}, i::Int, nm::Symbol)
    @boundscheck checkbounds(r, i)
    j = getcolumnmap(r)[i]
    @inbounds offlen = gettapes(r)[j][1]
    missingvalue(offlen) && return missing
    s = PointerString(pointer(getbuf(r), getpos(offlen)), getlen(offlen))
    return escapedvalue(offlen) ? unescape(s, gete(r)) : String(s)
end

@noinline stringsonly() = error("Parsers.parse only allowed on String column types")

Base.@propagate_inbounds function Parsers.parse(::Type{T}, r::Row2, i::Int) where {T}
    @boundscheck checkbounds(r, i)
    j = getcolumnmap(r)[i]
    typecode = gettypecodes(r)[j]
    (typecode == STRING || typecode == (STRING | MISSING)) || stringsonly()
    @inbounds offlen = gettapes(r)[j][1]
    missingvalue(offlen) && return missing
    pos = getpos(offlen)
    colopts = getcoloptions(r)
    opts = colopts === nothing ? getoptions(r) : colopts[j]
    x, code, vpos, vlen, tlen = Parsers.xparse(T, getbuf(r), pos, pos + getlen(offlen), opts)
    return Parsers.ok(code) ? x : missing
end

Base.@propagate_inbounds function detect(r::Row2, i::Int)
    @boundscheck checkbounds(r, i)
    j = getcolumnmap(r)[i]
    typecode = gettypecodes(r)[j]
    (typecode == STRING || typecode == (STRING | MISSING)) || stringsonly()
    @inbounds offlen = gettapes(r)[j][1]
    missingvalue(offlen) && return missing
    pos = getpos(offlen)
    colopts = getcoloptions(r)
    opts = colopts === nothing ? getoptions(r) : colopts[j]
    x = detect(getbuf(r), pos, pos + getlen(offlen) - 1, opts)
    return x === nothing ? r[i] : x
end

function Parsers.parse(::Type{T}, r::Row2, nm::Symbol) where {T}
    @inbounds x = Parsers.parse(T, r, getlookup(r)[nm])
    return x
end

function detect(r::Row2, nm::Symbol)
    @inbounds x = detect(r, getlookup(r)[nm])
    return x
end


text:
using GHCNData
using Test

@testset "GHCNData.jl" begin
end


text:
# To test the prediction with the previous implementation
# You need to write oraclepred, which predicts similar to oracleloss,
# Do not use invalid actions as well as gold actions
include("test_header.jl")

function load_data()
    d = load(language_model); v = create_vocab(d); wmodel = makewmodel(d);
    corpus = load_conllu(small_data_file, v)
    # newly added
    dev    = load_conllu(small_dev_file, v)
    corpora = []; push!(corpora, corpus); push!(corpora, dev); cc=vcat(corpora...);
    fillvecs!(wmodel, cc, v)

    # ppl = fillvecs!(wmodel, corpus, v) old
    unk = unkrate(corpus)
    info("Caching lm vectors...")
    map(cachelmvec!,corpora)
    #cachelmvec!(corpora[1]);  old

    info("Model initialization...")
    pmodel, poptim = initmodel1(odict, corpus[1])
    #p2 = odict[:arctype](corpus[2901]); pvecs = postagv(old_pmodel);
    return corpora, pmodel, poptim
end

function test_pred()
    corpora, model, poptims = load_data(); corpus=corpora[1];
    pdrop = (0.5, 0.6); batchsize = 16;

    @msg "Drop out probabilities $pdrop"
    info("calculating initial accuracies")
    acc1 = oracletest(model, corpora[1], odict[:arctype], odict[:feats], batchsize; pdrop=(0.0, 0.0))
    acc2 = oracletest(model, corpora[2], odict[:arctype], odict[:feats], batchsize; pdrop=(0.0, 0.0))
    @msg "Initial tracc $acc1 devacc $acc2"
    sentbatches = minibatch(corpus, batchsize, maxlen=64, minlen=2, shuf=false) # shuffling imp. acc
    nsent = sum(map(length,sentbatches)); nsent0 = length(corpus)
    nword = sum(map(length,vcat(sentbatches...))); nword0 = sum(map(length,corpus))
    for epoch=1:20
        @msg("nsent=$nsent/$nsent0 nword=$nword/$nword0")
        nwords = StopWatch()
        losses = Any[0,0,0]
        @time for sentences in sentbatches
            grads = oraclegrad(model, sentences, odict[:arctype], odict[:feats]; losses=losses, pdrop=pdrop)
            update!(model, grads, poptims)
            nw = sum(map(length,sentences))
            if (speed = inc(nwords, nw)) != nothing
                date("$(nwords.ncurr) words $(round(Int,speed)) wps $(losses[3]) avgloss")
                gc();Knet.gc();gc();
                #free_KnetArray()
            end

        end
        empty_parses!(corpus)
        acc = oracletest(model, corpus, odict[:arctype], odict[:feats], batchsize; pdrop=(0.0, 0.0))
        acc2 = oracletest(model, corpora[2], odict[:arctype], odict[:feats], batchsize; pdrop=(0.0, 0.0))

        println()
        println("epoch $epoch train acc $acc dev acc $acc2");
    end
#    p2 = odict[:arctype](corpus[2901]); pvecs = postagv(model);
#    idx = [ rand(1:length(corpus)) for i in 1:10 ];
end
!isinteractive() && test_pred()


text:
"""
    FieldBoundaryConditions

An alias for `NamedTuple{(:x, :y, :z)}` that represents a set of three `CoordinateBoundaryCondition`s
applied to a field along x, y, and z.
"""
const FieldBoundaryConditions = NamedTuple{(:x, :y, :z)}

"""
    FieldBoundaryConditions(x, y, z)

Construct a `FieldBoundaryConditions` using a `CoordinateBoundaryCondition` for each of the
`x`, `y`, and `z` coordinates.
"""
FieldBoundaryConditions(x, y, z) = FieldBoundaryConditions((x, y, z))

DefaultBoundaryCondition(::Union{Type{Grids.Periodic}, Type{Flat}}, loc) = PeriodicBoundaryCondition()

DefaultBoundaryCondition(::Type{Bounded}, ::Type{Cell}) = NoFluxBoundaryCondition()
DefaultBoundaryCondition(::Type{Bounded}, ::Type{Face}) = ImpenetrableBoundaryCondition()

function validate_bcs(topology, left_bc, right_bc, default_bc, left_name, right_name, dir)
    if topology isa Periodic && (left_bc != default_bc || right_bc != default_bc)
        e = "Cannot specify $left_name or $right_name boundary conditions with " *
            "$topology topology in $dir-direction."
        throw(ArgumentError(e))
    end
    return true
end

"""
    FieldBoundaryConditions(grid, loc;
          east = DefaultBoundaryCondition(topology(grid)[1], loc[1]),
          west = DefaultBoundaryCondition(topology(grid)[1], loc[1]),
         south = DefaultBoundaryCondition(topology(grid)[2], loc[2]),
         north = DefaultBoundaryCondition(topology(grid)[2], loc[2]),
        bottom = DefaultBoundaryCondition(topology(grid)[3], loc[3]),
           top = DefaultBoundaryCondition(topology(grid)[3], loc[3]))

Construct `FieldBoundaryConditions` for a field with location `loc` (a 3-tuple of `Face` or `Cell`)
defined on `grid` (the grid's topology is what defined the default boundary conditions that are
imposed).

Specific boundary conditions can be applied along the x dimension with the `west` and `east` kwargs,
along the y-dimension with the `south` and `north` kwargs, and along the z-dimension with the `bottom`
and `top` kwargs.
"""
function FieldBoundaryConditions(grid, loc;
      east = DefaultBoundaryCondition(topology(grid)[1], loc[1]),
      west = DefaultBoundaryCondition(topology(grid)[1], loc[1]),
     south = DefaultBoundaryCondition(topology(grid)[2], loc[2]),
     north = DefaultBoundaryCondition(topology(grid)[2], loc[2]),
    bottom = DefaultBoundaryCondition(topology(grid)[3], loc[3]),
       top = DefaultBoundaryCondition(topology(grid)[3], loc[3]))

    TX, TY, TZ = topology(grid)
    x_default_bc = DefaultBoundaryCondition(topology(grid)[1], loc[1])
    y_default_bc = DefaultBoundaryCondition(topology(grid)[2], loc[2])
    z_default_bc = DefaultBoundaryCondition(topology(grid)[3], loc[3])

    validate_bcs(TX, west,   east, x_default_bc, :west,   :east, :x)
    validate_bcs(TY, south, north, y_default_bc, :south, :north, :y)
    validate_bcs(TZ, bottom,  top, z_default_bc, :bottom,  :top, :z)

    x = CoordinateBoundaryConditions(west, east)
    y = CoordinateBoundaryConditions(south, north)
    z = CoordinateBoundaryConditions(bottom, top)

    return FieldBoundaryConditions(x, y, z)
end

  UVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Face, Cell, Cell); user_defined_bcs...)
  VVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell, Face, Cell); user_defined_bcs...)
  WVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell, Cell, Face); user_defined_bcs...)
     TracerBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell, Cell, Cell); user_defined_bcs...)
   PressureBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell, Cell, Cell); user_defined_bcs...)
DiffusivityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell, Cell, Cell); user_defined_bcs...)


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

using Markdown
using InteractiveUtils

# ╔═╡ 96d4f44e-49d1-4398-b7f8-4daa7290b287
begin
	using PlutoUI
	PlutoUI.TableOfContents()
end

# ╔═╡ 915d5707-daf4-4127-abf8-68a5528d9209
begin
	using ImageUtils
	using Plots

	p = shepp_logan(256)

	heatmap(reverse(p,dims=1), c=:viridis)
end

# ╔═╡ 9a1941ff-3d6c-40b3-b47b-254e859c5b6c
begin
	using FFTW

	pFT = fftshift( fft(p) )
	heatmap(reverse(log.(abs.(pFT).+1.1),dims=1), c=:viridis)
end

# ╔═╡ 71d11567-6098-4a0e-b1ad-5d0319d80ecf
begin
	using Wavelets

	pW = dwt(p, wavelet(WT.db2))
	heatmap(reverse(log.(abs.(pW).+1.1),dims=1), c=:viridis)
end

# ╔═╡ cd73e41b-283e-4fe8-84c8-ee86ce7968ff
begin
	using SparseArrays

	function myfunc(a::AbstractMatrix{T}) where T<:Number
  		# generic fallback
	end

	function myfunc(a::StridedMatrix{T}) where T<:Number
  		# call BLAS routines
	end

	function myfunc(a::SparseMatrixCSC{T}) where T<:Number
 		# exploit sparsity
	end
end

# ╔═╡ 1257912e-7eea-495b-9668-6b952cbd1b3b
begin
	using MRIReco

	# File loading
	f = BrukerFile("data/brukerfileCart")
	raw = RawAcquisitionData(f);
end

# ╔═╡ f12894c1-578d-47b4-8c71-5ccf31a5b84c
using LinearOperators, DelimitedFiles, LinearAlgebra

# ╔═╡ 2a2d075a-accb-11eb-36a3-5bacf8e2498d
md"""
# MRIReco.jl
#### Julia package for MRI reconstruction

by Tobias Knopp and Mirco Grosser
"""

# ╔═╡ 82453453-f2bb-423e-a7bb-099d1ed9dade
md"""
## Getting Started

This tutorial can be used and run in two ways:
* You can directly use the Pluto notebook you are reading right now. This needs some preparation discussed next.
* You can also copy/paste into a Julia REPL. This is in particular a good way if you are reading an HTML or PDF version of this tutorial.
* Keep in mind that this notebook will run all cells, which needs some time since we also showcase some more computationally intensive methods.
"""

# ╔═╡ 12dcb593-838f-4cf3-81e0-ca035dee56b5
md"""
### Preparation

* Go to [https://julialang.org/downloads/](https://julialang.org/downloads/) and download the most recent stable release (version 1.6 when this notebook was prepared).
* Checkout [https://github.com/MagneticResonanceImaging/MRIRecoTutorial](https://github.com/MagneticResonanceImaging/MRIRecoTutorial).
* Start Julia from within a terminal within the folder you checked out the repository.
* We need to install some packages. With ] we can enter the Pkg mode and enter
  ```julia
  (@v1.6) pkg> activate .
  ```
  which will install all necessary packages.
* To run this notebook use the command `import Pluto; Pluto.run()` and open the notebook in the appearing file browser.
* In the box *Open from file:* enter the filename of the notebook (MRIRecoTutorial.jl) and press enter.
* We are using [Pluto.jl](https://github.com/fonsp/Pluto.jl), which is similar to Jupyter but has some advantages, in particular the notebook is stored as plain Julia code that you can access with a simple text editor.
* Instead of running the notebook in the browser you can also run it in the REPL
  ```julia
  julia> include("MRIRecoTutorial.jl")
  ```
  Note that this, however will only show the plot of the last reconstruction since this notebook is tailored for usage in Pluto.
"""

# ╔═╡ c8b39aa6-d85a-4e18-822a-5dccc3b7e5c7
md"""
## Introduction to Julia

* High level programming language with Matlab like syntax
* Provides a powerful REPL for interactive use
* Strongly influenced by C++ and Python
* Generates efficient machine code:
  no need to write performance critical things in C/C++

In more detail:
* Julia uses JIT-compilation to generate efficient machine code
* Julia has a sophisticated type system $\rightarrow$ allows the compiler to identify statically typed code fragments $\rightarrow$ efficient machine code
* Dynamic program style is fine where convenience is more important than performance
* Package manager allows efficient modularization
"""

# ╔═╡ b6a623f0-499c-4bb7-81c6-679bd491452f
md"""
In a nutshell (somewhat biased...):

**Julia = (C++ $\cup$ Python $\cup$ Matlab) $\cap$ OnlyTheNiceThings**
"""

# ╔═╡ 056d2ff9-cc18-41bb-b35f-07936844cc60
md"""
### Julia Tooling

Commonly Julia is developed in

* Packages: Several source files / modules grouped together
* Scripts: Usually load packages and perform a dedicated operation

Let us explore some of the available functionality
"""

# ╔═╡ 56c864db-9b71-4d9d-81ee-3df9c839bfb4
md"""
Let's do some image processing
"""

# ╔═╡ 65ce8d0b-d4b7-4ce9-8954-803f5628d536
let
	pfilt = imfilter(p, Kernel.gaussian(3))
	heatmap(reverse(pfilt,dims=1).-0.1, c=:viridis)
end

# ╔═╡ ae8a1021-f85f-40a0-ad35-e571e974d3d2
md"Let's switch into frequency space"

# ╔═╡ bca47852-7812-4b3b-8f6a-b90dd727f31f
md"Or let's look at the wavelet transform"

# ╔═╡ f42ade9a-b8bc-499a-93f5-2cc077218e5b

# ╔═╡ 3ef6ac29-bcc6-4648-8c0c-e758ddd8f818
md"""
### Multiple Dispatch

Julia uses multiple dispatch to decide, which method is called
"""

# ╔═╡ e7d573d2-1bac-40c5-b7e6-2abdf94806c7
md"""
### How low level is Julia?

You can directly call C code with zero overhead
"""

# ╔═╡ 4adcb409-4632-4ddb-9e47-75e8d2dc4d83
t = ccall(:clock, Int32, ())

# ╔═╡ aeb4be94-bb33-44b0-9532-2d752fe2739b
md"""
One can look at the assembly code
"""

# ╔═╡ c1cea40b-6de6-4a83-a937-ece75dfca1a8
double(x) = 2*x

# ╔═╡ e1a3e63e-e3ce-4291-b54b-41db9d099c1d
with_terminal() do
  @code_native double(2)
end

# ╔═╡ 41e3e109-d1ad-4f51-835b-b4fa0f2c51ba
with_terminal() do
  @code_native double(2.0)
end

# ╔═╡ 7fac2740-1036-4cf6-88c9-9049f9ad7e41
md"""
## MRIReco: Architecture

* Julia package intended to provide a powerful MRI reconstruction toolbox
* In development (not as stable / feature complete as Gadgetron / BART)
* Focus on reconstruction researchers providing a hackable package
* Modular concept using existing goodies from Julialand
"""

# ╔═╡ 30660c21-4cce-4ed9-8b73-f68e60cbec92
md"""
### Dependency Graph
"""

# ╔═╡ f0cbdf18-6870-48be-8b42-a17305d6d7c1
md"The following shows a dependency graph of MRIReco and its sister package MPIReco (used for Magnetic Particle Imaging Reconstruction). This shows the philosophy to share common functionality into dedicated packages."

# ╔═╡ bfa0e160-a18f-4856-86c2-28485fa2e5e6
PlutoUI.LocalResource("./img/JuliaImagingPackages.png", :width => 700)

# ╔═╡ 2b7d6b3a-c4b9-4fb5-be8f-bc9ad4bd7723
md"""
### Data Types and Flow

Basic overview of the data flow during reconstruction
"""

# ╔═╡ 583111a4-2a79-4fe7-9122-294456c7478a
PlutoUI.LocalResource("./img/overview.svg", :width => 600)

# ╔═╡ f600f9b5-b759-4a54-abd9-cb3f14095583
md"## MRIReco: Examples"

# ╔═╡ bca790f3-da86-4be4-a97a-335741b259b5
md"""
### Example 1: Reconstruction Pipeline

Lets first get some data to perform a simple image reconstruction
"""

# ╔═╡ c3164c4a-8553-4188-aa15-2a7d9e9da106
md"""
The following code shows all the layers how to perform an image reconstruction and store the result:
```julia
using MRIReco

# File loading
f = BrukerFile("data/brukerfileCart")
raw = RawAcquisitionData(f)
acq = AcquisitionData(raw)

# Reconstruction
params = Dict{Symbol, Any}()
params[:reco] = "direct"
params[:reconSize] = (acq.encodingSize[1],acq.encodingSize[2])
Ireco = reconstruction(acq, params)

# Image handling
saveImage("reco.nii", Ireco, true);

```

Let's go through this step by step
"""

# ╔═╡ 984c46f2-e62c-4dac-97a6-eb9895d148a4
md"""
The `RawAcquisitionData` has a similar design as an ISMRMRD data file.
It has two members:
* `raw.params` is a dict that holds general parameter of the file
* `raw.profiles` is a vector with all measured profiles
"""

# ╔═╡ e9b2124f-5cb9-4705-80df-4e2f1b85d846
raw.params

# ╔═╡ 6b62ccaa-5ae8-43bf-a620-989af8754c6c
length(raw.profiles)

# ╔═╡ c948a18c-43c7-4bdc-a51b-591abcb57497
raw.profiles[1]

# ╔═╡ b572b96f-2bc8-4fc4-9113-d0a052c2bce4
md"""
#### Reconstruction

Next we will reconstruct the data. To do so we first convert the raw acquisition data
object into an acquisition data object. This will reorder the data in such a way, that
the reconstruction knows how to deal with it.
"""

# ╔═╡ 470557df-b938-4312-ba85-88e9de4a8564
acq = AcquisitionData(raw);

# ╔═╡ 78614fe5-e980-45fc-873a-105c8dba80e4
md" Now we can perform the reconstruction. It can be controlled using a dictionary of parameters"

# ╔═╡ 46ff4011-29c1-4c94-b794-0f655ec67daa
begin
	local params = Dict{Symbol, Any}()
	params[:reco] = "direct"
	params[:reconSize] = (acq.encodingSize[1],acq.encodingSize[2]);

	Ireco = reconstruction(acq, params);
end

# ╔═╡ 1e9452da-7710-4f3d-a389-2e0bc221acc7
md"The reconstructed image is stored as a 5-dimensional `AxisArray`. The latter is simply an array with named axis.

Upon inspection of the reconstructed data we see that the first three dimensions correspond to the three spatial dimensions, whereas the fourth dimension encodes different contrasts and the fifth dimension encodes different receive channels."

# ╔═╡ 1f762afd-12bf-43f1-9ea8-1aba7a0a41b0
md"The image data can be stored as a NIfTI file by calling"

# ╔═╡ dc0545e8-e309-480c-8f6e-3a517b0c7c4f
saveImage("reco.nii", Ireco, true);

# ╔═╡ 5f1d779b-7796-4103-b512-b06c8dfc8776
md"To load this image back just call"

# ╔═╡ 52f2408a-1cde-451f-8422-da5f4791e668
Ireco2 = loadImage("reco.nii")

# ╔═╡ 35c98353-b1fd-4916-b143-e27a23174716
md"Lets plot the reconstructed data. We will be using Plots.jl for this"

# ╔═╡ da39843a-d344-4509-b7c7-13ef9ed37a65
heatmap(reverse( circshift(abs.(Ireco[:,:,10,1,1]),(0,5)),dims=1), c=:viridis)

# ╔═╡ 44b4863b-cd5f-4bb1-80ed-fe11458adf84
md"Bruker does provide reconstruction data as well:"

# ╔═╡ a7c32c4b-0299-4bd6-b82f-f5ffb2ed9cbc
begin
	b = BrukerFile("data/brukerfileCart")
	Ibruker = recoData(b)

	heatmap(reverse( Ibruker[:,:,10,1,1]',dims=1), c=:viridis)
end

# ╔═╡ 617c5a93-36d8-424c-ab54-fa261be07439
md"""
#### Conversion

As we have loaded a `RawAcquistionData` from file, it is also possible to save it.
Currently, it is only possible to save into the ISMRMRD data format, which seems
to be the most sensible choice anyway.
"""

# ╔═╡ 023e6072-b02d-43d2-a958-711e989eaed7
begin
	fout = ISMRMRDFile("outputfile.h5")

	save(fout, raw);
end

# ╔═╡ 079bb385-df0a-47d1-ac9c-b352ae94f07f
md"We may load the data to check that it is still the same:"

# ╔═╡ fa80020c-1056-400d-a753-81993ace124e
begin
	rawCopy = RawAcquisitionData(fout)

	rawCopy.params
end

# ╔═╡ ed50e3d3-a228-4c51-a081-e98db113feb1
md"""
### Example 2: Simulation & Reconstruction

The following code shows a simple MRI simulation and reconstruction considering offresonance effects
"""

# ╔═╡ 107c4794-95a8-4e85-b789-2139091dc5de
begin
	N = 128

	I = MRIReco.shepp_logan(N)
	I = circularShutterFreq!(I,1)
	cmap = 1im*quadraticFieldmap(N,N,125*2pi)

	# simulation parameters
	local params = Dict{Symbol, Any}()
	params[:simulation] = "fast"
	params[:trajName] = "Spiral"
	params[:numProfiles] = 10
	params[:numSamplingPerProfile] = div(N*N,8)
	params[:windings] = 8
	params[:AQ] = 10.0e-3
	params[:correctionMap] = cmap[:,:,1]

	# do simulation
	acqData = simulation(I, params)
end

# ╔═╡ fc564eb0-ee21-43ea-8c65-f07677428ef2
begin
	tr = trajectory(acqData)
	nodes = kspaceNodes(tr)
	plot(nodes[1,1:2048],nodes[2,1:2048], label="trajectory")
end

# ╔═╡ 3ccafef9-5cb4-488f-a697-1c950399ee7c
begin
	# reco parameters
	local params = Dict{Symbol, Any}()
	params[:reco] = "standard"
	params[:regularization] = "L2"
	params[:iterations] = 3
	params[:solver] = "cgnr"
	params[:reconSize] = (N,N)

	INoCorr = reconstruction(acqData, params)

	params[:correctionMap] = cmap
	ICorr = reconstruction(acqData, params)

	heatmap(reverse( abs.(INoCorr[:,:,1,1,1]),dims=1), c=:viridis)
end

# ╔═╡ 6500d0af-d0f4-4dd7-a84a-e70b3424ee2f
heatmap(reverse( abs.(ICorr[:,:,1,1,1]),dims=1), c=:viridis)

# ╔═╡ 1c5c2818-fdc5-4214-bd3a-7c9a35e27153
md"""
### Example 3: Iterative Reconstruction
Iterative reconstructions can be performed using a high-level interface. All the relevant reconstruction parameters are collected in a `Dict{Symbol,Any}`. The reconstruction methods use those to formulate the image reconstruction problem and solve it.

To illustrate that, lets consider the following dataset, which contains one slice of a 3dFSE scan obtained from [mridata.org](mridata.org).
"""

# ╔═╡ f331e51f-a1e0-4959-8bc1-ef3b62d2d787
begin
	# load fully sampled data
	fKnee = ISMRMRDFile("data/knee_3dFSE_slice170.h5")
	acqDataKnee = AcquisitionData(fKnee);
end

# ╔═╡ f71ca727-1c64-41fd-a928-33eea4328598
md"For reference, let us reconstruct the fully sampled data"

# ╔═╡ 34f7fc8e-fd22-46e3-984a-940d504ce235
begin
	# reconstruct
	local params = Dict{Symbol, Any}()
	params[:reco] = "direct"
	params[:reconSize] = (320,320) # this size is also contained in acqData.encodingSize

	img = reconstruction(acqDataKnee, params)
	img = sqrt.(sum(img.^2,dims=5))

	# show image
	heatmap(reverse( abs.(img[:,:,1,1,1]),dims=1), c=:viridis)
end

# ╔═╡ c66c94b3-e5a2-4b87-a943-e69cb49092bd
md"To simulate an undersampled reconstruction, we retrospectively undersample the data using a Poisson disk pattern."

# ╔═╡ 96ed1677-3042-4aca-bb03-1719ebc0cba8
begin
	# sampling
	redFac= 4.0
	acqDataSub = sample_kspace(acqDataKnee,redFac,"poisson",calsize=30,profiles=false);

	# show sampling pattern
	msk = zeros(acqDataSub.encodingSize[1],acqDataSub.encodingSize[2])
	msk[acqDataSub.subsampleIndices[1]] .= 1

	heatmap(reverse( msk,dims=1), c=:grays)
end

# ╔═╡ 9592634e-ea3d-42eb-b013-cea1ff98601e
md"Estimate the coil sensitivities using ESPIRiT and reconstruct using SENSE"

# ╔═╡ 9aec589f-3615-49f6-9a92-c87afb7568a9
begin
	# coil sensitivities
	smaps = espirit(acqDataSub,(6,6),30,eigThresh_1=0.035,eigThresh_2=0.98)

	# SENSE reconstruction
	local params = Dict{Symbol, Any}()
	params[:reco] = "multiCoil"
	params[:reconSize] = (320,320)
	params[:senseMaps] = smaps

	params[:solver] = "cgnr"
	params[:regularization] = "L2"
	params[:λ] = 1.e1
	params[:iterations] = 5

	img_l2 = reconstruction(acqDataSub, params)

	# show image
	heatmap(reverse( abs.(img_l2[:,:,1,1,1]),dims=1), c=:viridis)
end

# ╔═╡ 936a7154-f601-4f06-ba41-e684e40677f3
md"Using ℓ¹-Wavelet-regularization recquires us to change some parameters."

# ╔═╡ d21a0949-6256-4325-84ca-3f78ab382a51
begin
	local params = Dict{Symbol, Any}()
	params[:reco] = "multiCoil"
	params[:reconSize] = (320,320)
	params[:senseMaps] = smaps

	params[:solver] = "admm"
	params[:regularization] = "L1"
	params[:sparseTrafo] = "Wavelet"
	params[:λ] = 2.e1
	params[:iterations] = 50
	params[:ρ] = 1.0
	params[:absTol] = 1.e-4
	params[:relTol] = 1.e-2
	params[:tolInner] = 1.e-2

	img_wavelet = reconstruction(acqDataSub, params)

# show image
heatmap(reverse( abs.(img_wavelet[:,:,1,1,1]),dims=1), c=:viridis)
end

# ╔═╡ 63d00efb-c10e-4c8a-8160-d6156a6fcfd2
md"""
### Example 4: Extending MRIReco

Reconstructions parameters, such as sparsifying transforms or the form of regularization can be incorporated in two ways:
* specify name of of an already implemented transform or regularizer
* directly pass the `Regularization` or `LinearOperator` object to be used

The second option yields a simple way to
* modify reconstructions algorithms
* implement and test new reconstruction methods

#### Learned Sparsifying Transform

Let us implement a simple data-driven sparsifying transform (see: S. Ravishankar and Y. Bresler, IEEE Trans. Med. Imaging, 30 (5), 2011).
* Learn a dictionary from a reference image (e.g. adjacent slice) using KSVD  $\color{green}\checkmark$
* Implement sparsifying transform which analyses the input image in terms of the dictionary
"""

# ╔═╡ 873a78e3-f4cd-425a-8507-0d24cef8b771
PlutoUI.LocalResource("img/dictTrafo.png", :width => 500)

# ╔═╡ 484b7d08-6af4-4795-b7d6-e741ff3fb503
function analyzeImage(x::Vector{T}, D::Matrix{T}, xsize::NTuple{2,Int64},
		              psize::NTuple{2,Int64}; t0::Int64=size(D,2),
		              tol=1e-3, nmax::Int64=size(D,1)) where T
	nx,ny = xsize
	px,py = psize
	x = reshape(x,nx,ny)
	x_pad = repeat(x,2,2)[1:nx+px-1,1:ny+py-1] # pad image using periodic boundary conditions
	α = zeros(T,size(D,2),nx,ny)
	patch = zeros(T,px*py)
	for j=1:ny, i=1:nx
		patch[:] .= vec(x_pad[i:i+px-1,j:j+py-1])
		norm(patch)==0 && continue
		# matchingpursuit is contained in Wavelets.jl
		α[:,i,j] .= matchingpursuit(patch, x->D*x, x->transpose(D)*x, tol, nmax)
	end
	return vec(α)
end

# ╔═╡ 316f47e9-58a8-45b7-89d4-375a8c0024a2
function synthesizeImage(α::Vector{T}, D::Matrix{T}, xsize::NTuple{2,Int64},
		                 psize::NTuple{2,Int64}) where T
	nx,ny = xsize
	px,py = psize
	x = zeros(T,nx+px-1,ny+py-1)
	α = reshape(α,:,nx,ny)
	for j=1:ny, i=1:nx
		x[i:i+px-1,j:j+py-1] .+= reshape(D*α[:,i,j],px,py)
	end
	return vec(x[1:nx,1:ny])/(px*py)
end

# ╔═╡ 41b190f7-f9f6-414b-b260-3c5316c08619
function dictOp(D::Matrix{T},xsize::NTuple{2,Int64},psize::NTuple{2,Int64},
		        tol::Float64=1.e-3,nmax::Int64=size(D,1)) where T
	produ = x->analyzeImage(x,D,xsize,psize,tol=tol,nmax=nmax)
	ctprodu = x->synthesizeImage(x,D,xsize,psize)
	return LinearOperator(prod(xsize)*size(D,2),prod(xsize),false,false
			, produ
			, nothing
			, ctprodu )
end

# ╔═╡ 38b1f34c-c71a-4f8f-a672-27e745898deb
md"Next, let us Load the Dictionary, build the sparsifying transform and perform reconstruct the knee data"

# ╔═╡ e85a04b1-7a0d-4b51-a18c-48e1cce6a09e
begin
	# load dict
	D = ComplexF64.(readdlm("data/kneeDict165.tsv"))

	local params = Dict{Symbol, Any}()
	params[:reco] = "multiCoil"
	params[:reconSize] = (320,320)
	params[:senseMaps] = smaps

	params[:solver] = "admm"
	params[:regularization] = "L1"
	params[:sparseTrafo] = dictOp(D,(320,320),(6,6),1.e-3,12)
	params[:λ] = 2.e1
	params[:iterations] = 30
	params[:ρ] = 1.0
	params[:absTol] = 1.e-4
	params[:relTol] = 1.e-2
	params[:tolInner] = 1.e-2

	img_dict = reconstruction(acqDataSub, params)

	# show image
	heatmap(reverse( abs.(img_dict[:,:,1,1,1]),dims=1), c=:viridis)
end

# ╔═╡ 2514b8f8-2492-41b1-8c85-0debe400af12
md"""
#### Developing MRIReco.jl

Changing MRIReco.jl is pretty easy. Go into the package manager and develop the package

```julia
] dev MRIReco
```
After developing the package, the code is located in

   `~/.julia/dev/MRIReco`

and can be directly modified.

Note that Julia requires a restart, if you change code. This limitation
can be circumvented by using the `Revise.jl` package (look at the Github page).

"""

# ╔═╡ f78b55aa-7f5b-4412-a293-0eccac11dd91
md"## Summary
MRIReco is a MRI reconstruction package written in julia, which aims to:
* be fast
* accessible
* easily extendable

The package is still relatively young. Thus, while being usable, further updates and extensions can be expected.

To further improve the package, we welcome feedback, bug reports and ideas for new features to be added.

For further details on *MRIReco.jl*, have a look at:
* The packages github repository:  [https://github.com/MagneticResonanceImaging/MRIRecoTutorial](https://github.com/MagneticResonanceImaging/MRIRecoTutorial)
* The accompanying paper: [https://onlinelibrary.wiley.com/doi/full/10.1002/mrm.28792](https://onlinelibrary.wiley.com/doi/full/10.1002/mrm.28792)
"

# ╔═╡ Cell order:
# ╟─2a2d075a-accb-11eb-36a3-5bacf8e2498d
# ╟─96d4f44e-49d1-4398-b7f8-4daa7290b287
# ╟─82453453-f2bb-423e-a7bb-099d1ed9dade
# ╟─12dcb593-838f-4cf3-81e0-ca035dee56b5
# ╟─c8b39aa6-d85a-4e18-822a-5dccc3b7e5c7
# ╟─b6a623f0-499c-4bb7-81c6-679bd491452f
# ╟─056d2ff9-cc18-41bb-b35f-07936844cc60
# ╠═915d5707-daf4-4127-abf8-68a5528d9209
# ╟─56c864db-9b71-4d9d-81ee-3df9c839bfb4
# ╠═65ce8d0b-d4b7-4ce9-8954-803f5628d536
# ╟─ae8a1021-f85f-40a0-ad35-e571e974d3d2
# ╠═9a1941ff-3d6c-40b3-b47b-254e859c5b6c
# ╟─bca47852-7812-4b3b-8f6a-b90dd727f31f
# ╠═71d11567-6098-4a0e-b1ad-5d0319d80ecf
# ╠═f42ade9a-b8bc-499a-93f5-2cc077218e5b
# ╟─3ef6ac29-bcc6-4648-8c0c-e758ddd8f818
# ╠═cd73e41b-283e-4fe8-84c8-ee86ce7968ff
# ╟─e7d573d2-1bac-40c5-b7e6-2abdf94806c7
# ╠═4adcb409-4632-4ddb-9e47-75e8d2dc4d83
# ╟─aeb4be94-bb33-44b0-9532-2d752fe2739b
# ╠═c1cea40b-6de6-4a83-a937-ece75dfca1a8
# ╠═e1a3e63e-e3ce-4291-b54b-41db9d099c1d
# ╠═41e3e109-d1ad-4f51-835b-b4fa0f2c51ba
# ╟─7fac2740-1036-4cf6-88c9-9049f9ad7e41
# ╟─30660c21-4cce-4ed9-8b73-f68e60cbec92
# ╟─f0cbdf18-6870-48be-8b42-a17305d6d7c1
# ╟─bfa0e160-a18f-4856-86c2-28485fa2e5e6
# ╟─2b7d6b3a-c4b9-4fb5-be8f-bc9ad4bd7723
# ╟─583111a4-2a79-4fe7-9122-294456c7478a
# ╟─f600f9b5-b759-4a54-abd9-cb3f14095583
# ╟─bca790f3-da86-4be4-a97a-335741b259b5
# ╟─c3164c4a-8553-4188-aa15-2a7d9e9da106
# ╠═1257912e-7eea-495b-9668-6b952cbd1b3b
# ╟─984c46f2-e62c-4dac-97a6-eb9895d148a4
# ╠═e9b2124f-5cb9-4705-80df-4e2f1b85d846
# ╠═6b62ccaa-5ae8-43bf-a620-989af8754c6c
# ╠═c948a18c-43c7-4bdc-a51b-591abcb57497
# ╟─b572b96f-2bc8-4fc4-9113-d0a052c2bce4
# ╠═470557df-b938-4312-ba85-88e9de4a8564
# ╟─78614fe5-e980-45fc-873a-105c8dba80e4
# ╠═46ff4011-29c1-4c94-b794-0f655ec67daa
# ╟─1e9452da-7710-4f3d-a389-2e0bc221acc7
# ╟─1f762afd-12bf-43f1-9ea8-1aba7a0a41b0
# ╠═dc0545e8-e309-480c-8f6e-3a517b0c7c4f
# ╟─5f1d779b-7796-4103-b512-b06c8dfc8776
# ╠═52f2408a-1cde-451f-8422-da5f4791e668
# ╟─35c98353-b1fd-4916-b143-e27a23174716
# ╠═da39843a-d344-4509-b7c7-13ef9ed37a65
# ╟─44b4863b-cd5f-4bb1-80ed-fe11458adf84
# ╠═a7c32c4b-0299-4bd6-b82f-f5ffb2ed9cbc
# ╟─617c5a93-36d8-424c-ab54-fa261be07439
# ╠═023e6072-b02d-43d2-a958-711e989eaed7
# ╟─079bb385-df0a-47d1-ac9c-b352ae94f07f
# ╠═fa80020c-1056-400d-a753-81993ace124e
# ╟─ed50e3d3-a228-4c51-a081-e98db113feb1
# ╠═107c4794-95a8-4e85-b789-2139091dc5de
# ╠═fc564eb0-ee21-43ea-8c65-f07677428ef2
# ╠═3ccafef9-5cb4-488f-a697-1c950399ee7c
# ╠═6500d0af-d0f4-4dd7-a84a-e70b3424ee2f
# ╟─1c5c2818-fdc5-4214-bd3a-7c9a35e27153
# ╠═f331e51f-a1e0-4959-8bc1-ef3b62d2d787
# ╟─f71ca727-1c64-41fd-a928-33eea4328598
# ╠═34f7fc8e-fd22-46e3-984a-940d504ce235
# ╟─c66c94b3-e5a2-4b87-a943-e69cb49092bd
# ╠═96ed1677-3042-4aca-bb03-1719ebc0cba8
# ╟─9592634e-ea3d-42eb-b013-cea1ff98601e
# ╠═9aec589f-3615-49f6-9a92-c87afb7568a9
# ╟─936a7154-f601-4f06-ba41-e684e40677f3
# ╠═d21a0949-6256-4325-84ca-3f78ab382a51
# ╟─63d00efb-c10e-4c8a-8160-d6156a6fcfd2
# ╟─873a78e3-f4cd-425a-8507-0d24cef8b771
# ╠═f12894c1-578d-47b4-8c71-5ccf31a5b84c
# ╠═484b7d08-6af4-4795-b7d6-e741ff3fb503
# ╠═316f47e9-58a8-45b7-89d4-375a8c0024a2
# ╠═41b190f7-f9f6-414b-b260-3c5316c08619
# ╟─38b1f34c-c71a-4f8f-a672-27e745898deb
# ╠═e85a04b1-7a0d-4b51-a18c-48e1cce6a09e
# ╟─2514b8f8-2492-41b1-8c85-0debe400af12
# ╟─f78b55aa-7f5b-4412-a293-0eccac11dd91


text:
using Bezier
using Base.Test

# linear interpolation
P11 = [0., 1]
P12 = [2., 3]
@test [0,     1] == bezier(0,    P11, P12)
@test [0.5, 1.5] == bezier(0.25, P11, P12)
@test [1,     2] == bezier(0.5,  P11, P12)
@test [1.5, 2.5] == bezier(0.75, P11, P12)
@test [2,     3] == bezier(1,    P11, P12)

P21 = [0. 1]
P22 = [2. 3]
@test [0     1] == bezier(0,    P21, P22)
@test [0.5 1.5] == bezier(0.25, P21, P22)
@test [1     2] == bezier(0.5,  P21, P22)
@test [1.5 2.5] == bezier(0.75, P21, P22)
@test [2     3] == bezier(1,    P21, P22)

# quadratic interpolation with vector of vectors
@test Array[[10, 5]] == bezier(0.5, Vector{Float64}[[0., 0], [10, 10], [20, 0]])

# quadratic interpolation with matrix
@test [[10 5];] == bezier(0.5, [[0. 0]; [10 10]; [20 0]])
@test [[10 5];] == bezier(0.5, [0. 0; 10 10; 20 0])

# boundary checks
@test_throws DomainError bezier(prevfloat(0.), [0., 1], [2., 3])
@test_throws DomainError bezier(nextfloat(1.), [0., 1], [2., 3])
@test_throws DomainError bezier(prevfloat(0.), Vector{Float64}[[0., 1], [2, 3]])
@test_throws DomainError bezier(nextfloat(1.), Vector{Float64}[[0., 1], [2, 3]])
@test_throws DomainError bezier(prevfloat(0.), [0. 1; 2 3])
@test_throws DomainError bezier(nextfloat(1.), [0. 1; 2 3])


text:
module LeetCode

const LC = LeetCode
export LC

using Markdown
using DataStructures
include("common.jl")
include("problems/problems.jl")

# export all
for n in names(@__MODULE__(); all=true)
    if Base.isidentifier(n) &&
       !startswith(String(n), "_") &&
       n ∉ (Symbol(@__MODULE__()), :eval, :include)
        @eval export $n
    end
end

end


text:
using TupleLabXD
using Documenter

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

makedocs(;
    modules=[TupleLabXD],
    authors="PhysicsCodesLab",
    repo="https://github.com/PhysicsCodesLab/TupleLabXD.jl/blob/{commit}{path}#{line}",
    sitename="TupleLabXD.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://PhysicsCodesLab.github.io/TupleLabXD.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/PhysicsCodesLab/TupleLabXD.jl",
    devbranch="master",
)


text:
using Revise
using POMDPs
using GridGoalMDPs
using DiscreteValueIteration
using POMDPGifs

m = TwoAgentMDP(size=(25,25))
solver = SparseValueIterationSolver(verbose=true)
p = solve(solver, m)

gif = makegif(m, p, max_steps=100, show_progress=true)
fname = gif.filename
@info("Attempting to open $(fname) with google chrome.")
if Sys.iswindows()
    run(`cmd /C start chrome "$fname"`)
elseif Sys.isapple()
    run(`open -a "Google Chrome" $fname`)
else
    run(`google-chrome $fname`)
end


text:
if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

using Boltzmann
using Benchmark
using DataFrames    # Needed for organization of benchmark results

#=
Available test groups are:
* :unit,
* :integration,
* :acceptance,
* :benchmark,
=#
const DEFAULT_TEST_GROUPS = Set([:unit, :integration])
TEST_GROUPS = DEFAULT_TEST_GROUPS

# If the JULIA_TEST_GROUPS ENV variable is set
# then use the inte
if haskey(ENV, "JULIA_TEST_GROUPS")
    TEST_GROUPS = Set(
        Symbol[split(ENV["JULIA_TEST_GROUPS"])...]
    )
end

tests = [
    "utils",
    "rbm",
    "nets",
    "conditional",
]

for t in tests
    println("\n\nRunning $t tests...\n")
    include(string("test", t, ".jl"))
end


text:
"Mean for circular data"
circmean(α::AbstractVector,w=ones(length(α)))=sum(w.*exp.(im*α))
"""
Mean Resultant Vector Length for circular data
    α: sample of angles in radians
    w: number of incidences in case of binned angle data
    d: spacing of bin centers for binned data, if supplied
       correction factor is used to correct for bias
"""
function circr(α::AbstractVector,w=ones(length(α)),d=nothing)
    r = abs(sum(w.*exp.(im*α)))/sum(abs.(w))
    if d!=nothing
        r *= d/2/sin(d/2)
    end
    r
end

"Circular Variance"
circvar(α::AbstractVector,w=ones(length(α)),d=nothing) = 1-circr(α,w,d)

"""
Transforms p-axial data to a common scale
  α: sample of angles in radians
  p: number of modes
"""
circaxial(α,p=2)=(p*α)%(2pi)

"Mean direction for circular data with axial"
circaxialmean(α::AbstractVector,p=2)=angle(sum(exp.(im*α*p)))/p

"""
Hotelling's T-Squared test for one multivariate sample.

Performs Hotelling's T^2 test on multivariate samples X to determine
if the data have mean MU.  X should be a NxP matrix with N observations
of P-dimensional data, and the mean MU to be tested should be 1xP.
the significance level α is 0.05 by default.

H is 1 if the null hypothesis (that the mean of X is equal to MU) can be
rejected at significance level α.  P is the actual P value.

The code is based on HotellingT2.m by A. Trujillo-Ortiz and
R. Hernandez-Walls.
"""
function hotellingt2test(X,mu,α = 0.05)
	n=size(X,1); p=size(X,2); m = mean(X,dims=1); S = cov(X);
	T2 = n*(m-mu)*inv(S)*(m-mu)'  # T-square statistics

	if n>= 50 # Chi-square approximation
		P = 1-chisqcdf(T2[1],p)
	else  # F approximation
		F = (n-p)/((n-1)*p)*T2[1]
		P = ccdf(FDist(p, n-p), F)
	end
	H = P<α
	return P
end

"""
Direction tuning significance using dot product with empirical orientation preference
This function calculates the probability that the "true" direction tuning
vector of a neuron has non-zero length. It performs this by empirically
determing the unit orientation vector, and then computing the dot
product of the direction vector for each trial onto the overall orientation
vector, and then looking to see if the average is non-zero.

Inputs:  ANGLES is a vector of direction angles at which the response has
       been measured.
       RATES is the response of the neuron in to each angle; each row
       should contain responses from a different trial.
Output:  P the probability that the "true" direction tuning vector is
       non-zero.

See: Mazurek, Kagan, Van Hooser 2014;  Frontiers in Neural Circuits
"""
function dirsigtest(diraxis, dirvec)
	 #  compute the trial by trial dot products (really the square of the dot product)
	dot_prods = [real(dirvec) imag(dirvec)] * [real(diraxis) imag(diraxis)]'
	 # now compute p value
	p=pvalue(OneSampleTTest(dot_prods[:], 0))

	return p
end

"""
This function calculates the ROC curve, which represents the 1-specificity and sensitivity of two classes of data, (i.e., class_1 and class_2).

The function also returns all the needed quantitative parameters: threshold position, distance to the optimum point, sensitivity,
	specificity, accuracy, area under curve (AROC), positive and negative predicted values (PPV, NPV), false negative and
	positive rates (FNR, FPR), false discovery rate (FDR), false omission rate (FOR), F1 score, Matthews correlation coefficient
	(MCC), Informedness (BM) and Markedness; as well as the number of true positives (TP), true negatives (TN), false positives (FP),
	and false negatives (FN).

Rewrite from Víctor Martínez-Cagigal (2020). ROC Curve, MATLAB Central File Exchange. Retrieved February 11, 2020.
"""
function roccurve(class_1, class_2)

	# Calculating the threshold values between the data points
	s_data = unique(sort(hcat(class_1, class_2),dims=1))  # Sorted data points
	# s_data(isnan(s_data)) = []                 # Delete NaN values
	d_data = diff(s_data)                      # Difference between consecutive points
	if isempty(d_data)
		AROC = 0
		return AROC
	    # @warn "Both class data are the same!"
	end
	d_data = vcat(d_data,d_data[end]) # Last point
	thres=[];
	thres = append!(thres,s_data[1] - d_data[1])                 # First point  # Threshold values
	thres = vcat(thres,s_data + d_data./2)
	# Calculating the sensibility and specificity of each threshold
	curve = zeros(size(thres,1),2)
	distance = zeros(size(thres,1),1)
	a=1
	for i = 1:length(thres)
	    TP = sum(class_2 .>= thres[i])    # True positives
	    FP = sum(class_1 .>= thres[i])    # False positives
	    FN = sum(class_2 .< thres[i])     # False negatives
	    TN = sum(class_1 .< thres[i])     # True negatives

	    curve[i,1] = TP/(TP + FN)   # Sensitivity
	    curve[i,2] = TN/(TN + FP)	# Specificity

	    # Distance between each point and the optimum point (0,1)
	    distance[i]= sqrt((1-curve[i,1])^2+(curve[i,2]-1)^2)
	end

	# Optimum threshold and parameters
	opt = findmin(distance)[2][1]
	TP = sum(class_2 .>= thres[opt])    # No. true positives
	FP = sum(class_1 .>= thres[opt])    # No. false positives
	FN = sum(class_2 .< thres[opt])     # No. false negatives
	TN = sum(class_1 .< thres[opt])     # No. true negatives

	# Output parameters
	# param = DataFrame()
	# param.Threshold = thres[opt]                #Optimum threshold position
	# param.Sensi = curve[opt,1]                 # Sensitivity
	# param.Speci = curve[opt,2]                 # Specificity
	# param.AROC  = abs(trapz(1 .- curve[:,2], curve[:,1])) # Area under curve
	# param.Accuracy = (TP+TN)/(TP+TN+FP+FN)     # Aaccuracy
	# param.PPV   = TP/(TP+FP)                   # Positive predictive value
	# param.NPV   = TN/(TN+FN)                   # Negative predictive value
	# param.FNR   = FN/(FN+TP)                   # False negative rate
	# param.FPR   = FP/(FP+TN)                   # False positive rate
	# param.FDR   = FP/(FP+TP)                   # False discovery rate
	# param.FOR   = FN/(FN+TN)                   # False omission rate
	# param.F1_score = 2*TP/(2*TP+FP+FN)         # F1 score
	# param.MCC   = (TP*TN-FP*FN)/sqrt((TP+FP)*(TP+FN)*(TN+FP)*(TN+FN))  # Matthews correlation coefficient
	# param.BM    = curve[opt,1] + curve[opt,2] - 1    # Informedness
	# param.MK    = TP/(TP+FP) + TN/(TN+FN) - 1        # Markedness
	#
	# param.TP = TP    # No. true positives
	# param.FP = FP    # No. false positives
	# param.FN = FN    # No. false negatives
	# param.TN = TN    # No. true negatives
    AROC  = abs(trapz(1 .- curve[:,2], curve[:,1]))  # so far , only need AUC
	return AROC
end


text:
#
#  InterMonomerForces.jl
#  CollagenPolymer
#
#  Created by Christopher Revell on 30/03/2020.
#
#

module InterMonomerForces

# Import Julia packages
using LinearAlgebra
using StaticArrays

# Import local program modules
using LennardJones

function interMonomerForces!(pairsList,pos,F,nDomains,ϵ,σ,dx,WCAthreshSq,intrctnThrshld)

    for (ii,jj) in pairsList
        if floor(Int8,(ii-1)/nDomains)==floor(Int8,(jj-1)/nDomains) && abs(ii-jj)<=1
            # Skip adjacent particles in same monomer
        else
            dx = pos[jj] - pos[ii]
            dxMag = norm(dx)
            if dxMag < intrctnThrshld
                dx = dx*lennardJones(dxMag,ϵ,σ)/dxMag
                F[ii] += dx
                F[jj] -= dx
            end

        	# if (ii+3)%nDomains == (jj-1)%nDomains && floor(Int8,(ii-1)/nDomains)!=floor(Int8,(jj-1)/nDomains)
            # 	# Apply adhesive van der waals force in stepped fashion between monomers
            #     dx = dx*lennardJones(dx,ϵ,σ)
            #     F[ii] += dx
            #     F[jj] -= dx
            # elseif dxmag_sq < WCAthreshSq
            #     # For all other particles, apply WCA potential (truncated repulsive Lennard-Jones)
            #     dx = dx*lennardJones(dx,ϵ,σ)
            #     F[ii] += dx
            #     F[jj] -= dx
            # else
            #     # Skip any pairs within interaction range, beyond WCA range, and without specified adhesive rule
            # end
        end
    end

    return nothing

end

export interMonomerForces!

end

max_stars_repo_name: chris-revell/collagen-model-jl

max_issues_repo_name: chris-revell/collagen-model-jl

max_forks_repo_name: chris-revell/collagen-model-jl


text:
using Test
import Pluto: Notebook, ServerSession, ClientSession, Cell, update_caches!, updated_topology, is_just_text

@testset "Analysis" begin
    notebook = Notebook([
        Cell(""),
        Cell("md\"a\""),
        Cell("html\"a\""),
        Cell("md\"a \$b\$\""),
        Cell("""
        let
            x = md"a"
            md"r \$x"
        end
        """),
        Cell("html\"a \$b\""),

        Cell("md\"a \$b\""),
        Cell("@a md\"asdf\""),
        Cell("x()"),
        Cell("x() = y()"),
        Cell("1 + 1"),
        Cell("import Dates"),
        Cell("import Dates"),
        Cell("while false end"),
        Cell("for i in [1,2]; end"),
        Cell("[i for i in [1,2]]"),

    ])

    update_caches!(notebook, notebook.cells)
    old = notebook.topology
    new = notebook.topology = updated_topology(old, notebook, notebook.cells)

    @testset "Only-text detection" begin
        @test is_just_text(new, notebook.cells[1])
        @test is_just_text(new, notebook.cells[2])
        @test is_just_text(new, notebook.cells[3])
        @test is_just_text(new, notebook.cells[4])
        @test is_just_text(new, notebook.cells[5])
        @test is_just_text(new, notebook.cells[6])

        @test !is_just_text(new, notebook.cells[7])
        @test !is_just_text(new, notebook.cells[8])
        @test !is_just_text(new, notebook.cells[9])
        @test !is_just_text(new, notebook.cells[10])
        @test !is_just_text(new, notebook.cells[11])
        @test !is_just_text(new, notebook.cells[12])
        @test !is_just_text(new, notebook.cells[13])
        @test !is_just_text(new, notebook.cells[14])
        @test !is_just_text(new, notebook.cells[15])
        @test !is_just_text(new, notebook.cells[16])
    end
end


text:
# RUN update-readme.sh MANUALLY

using Documenter, DocumenterMarkdown, HDF5Logging

makedocs(
    modules = [HDF5Logging],
    format = DocumenterMarkdown.Markdown(),
    authors = "Tamas K. Papp",
    sitename = "HDF5Logging.jl",
    pages = Any["index.md"],
    # strict = true,
    # clean = true,
    checkdocs = :exports,
)


text:
#---
#julia: 1.2.3
#---


text:
# check if a list of indices specifies a tensor contraction in ncon style
function isnconstyle(network)
    allindices = Vector{Int}()
    for ind in network
        all(i->isa(i, Integer), ind) || return false
        append!(allindices, ind)
    end
    while length(allindices) > 0
        i = pop!(allindices)
        if i > 0 # positive labels represent contractions or traces and should appear twice
            k = findfirst(isequal(i), allindices)
            k === nothing && return false
            l = findnext(isequal(i), allindices, k+1)
            l !== nothing && return false
            deleteat!(allindices, k)
        elseif i < 0 # negative labels represent open indices and should appear once
            findfirst(isequal(i), allindices) === nothing || return false
        else # i == 0
            return false
        end
    end
    return true
end

function ncontree(network)
    contractionindices = Vector{Vector{Int}}(undef, length(network))
    for k = 1:length(network)
        indices = network[k]
        # trace indices have already been removed, remove open indices by filtering on positive values
        contractionindices[k] = filter(i->i>0, indices)
    end
    partialtrees = collect(Any, 1:length(network))
    _ncontree!(partialtrees, contractionindices)
end

function _ncontree!(partialtrees, contractionindices)
    if length(partialtrees) == 1
        return partialtrees[1]
    end
    if all(isempty, contractionindices) # disconnected network
        partialtrees[end-1] = Any[partialtrees[end-1], partialtrees[end]]
        pop!(partialtrees)
        pop!(contractionindices)
    else
        let firstind = minimum(vcat(contractionindices...))
            i1 = findfirst(x->in(firstind,x), contractionindices)
            i2 = findnext(x->in(firstind,x), contractionindices, i1+1)
            @assert i1 !== nothing && i2 !== nothing
            newindices = unique2(vcat(contractionindices[i1], contractionindices[i2]))
            newtree = Any[partialtrees[i1], partialtrees[i2]]
            partialtrees[i1] = newtree
            deleteat!(partialtrees, i2)
            contractionindices[i1] = newindices
            deleteat!(contractionindices, i2)
        end
    end
    _ncontree!(partialtrees, contractionindices)
end

function nconindexcompletion(ex)
    if isassignment(ex) || isdefinition(ex)
        lhs, rhs = getlhs(ex), getrhs(ex)
        # process left hand side
        if istensor(lhs) && istensorexpr(rhs)
            indices = getindices(rhs)

            if lhs.head == :ref && length(lhs.args) == 2 && lhs.args[2] == :(:)
                if all(isa(i, Integer) && i < 0 for i in indices)
                    lhs = Expr(:ref, lhs.args[1], sort(indices, rev=true)...)
                else
                    error("cannot automatically infer index order of left hand side")
                end
            end
            return Expr(ex.head, lhs, rhs)
        else
            return ex
        end
    elseif ex isa Expr
        return Expr(ex.head, map(nconindexcompletion, ex.args)...)
    else
        return ex
    end
end