text:
using HDF5, Compat, KernelDensity

"""
    cal_data(data::Array; bit = 1, nqubits = 1, num_repeats = 2, rm_cals=true)

Normalize data using calibration points.  Note the data should not be complex.

# Arguments
  - `data::Array{Float64}`: data to normalize with cals at the end
  - `bit::Int`: 1, 2, 3... from LSB to MSB
  - `nqubits::Int`: number of qubits
  - `num_repeats::Int`: number of calibration points per computational state
  - `rm_cals::boolean`: remove calibration points from the end of the sequences

# Examples
```julia-repl
julia> data = vcat(rand(10,1), [0,0,1,1])
julia> Qlab.cal_data(data)
-0.037458  -0.763034  -0.830473  -1.69414  …  -1.25552  -2.45348  -2.45213
```
"""
function cal_data(data::Array;
                  bit = 1,
                  nqubits = 1,
                  num_repeats = 2,
                  rm_cals = true)
	zero_cal = mean(data[end-num_repeats*(2^nqubits)+1:end-num_repeats*(2^nqubits-1)])
	one_cal = mean(data[end-num_repeats*(2^nqubits-2^(bit-1))+1:end-num_repeats*(2^nqubits-2^(bit-1)-1)])
	scale_factor = -(one_cal - zero_cal)/2;
	if rm_cals
		data = data[1:end-(2^nqubits)*num_repeats]
	end
	data = (data .- zero_cal)/scale_factor .+ 1
end

"""
  cal_data(data; qubit="", cal0="0", cal1="1")

Normalize data with reference measurements defined in metadata. Format used in
Auspex.  Assume data with structure: [data, data, ...., cal0, cal0, ..., cal1,
cal1, cal1, ..., data, ..., cal0, ...].  The number of data and cals can be
different for each set, as long as they are contiguous.

# Arguments
  - `data`: dictionary (Auspex format)
  - `qubit::string`: qubit name
  - `cal0/1::string`: reference measurement for qubit in 0/1
"""

function cal_data(data::Tuple{Dict{String,Array{Any,N} where N},Dict{String,Any}};
                  qubit::String = "",
                  cal0::String = "0",
                  cal1::String = "1",
                  quad = :real)
    if length(collect(keys(data[1]))) == 1
        qubit = collect(keys(data[1]))[1]
    elseif isempty(qubit)
        error("More than one qubit. Specify qubit name and calibration labels")
    elseif ~(qubit in keys(data[1]))
        error("Qubit not found")
    end

    metadata = data[2][qubit]["meta_data"]
    if(length(metadata))>1
       error("Invalid metadata")
    end

    metadata_values = data[2][qubit]["meta_data"][(collect(keys(metadata)))[1]]

    data_out = []
    ind0 = findall(x -> x==cal0, metadata_values)
    ind1 = findall(x -> x==cal1, metadata_values)
    ind_data = findall(x -> x==metadata_values[1], metadata_values)
    ind0_edge = push!(filter(x -> ind0[x] != ind0[x+1]-1, 1:length(ind0)-1), length(ind0)) #find consecutive cals
    ind1_edge = push!(filter(x -> ind1[x] != ind1[x+1]-1, 1:length(ind1)-1), length(ind1))
    ind_data_edge = push!(filter(x -> ind_data[x] != ind_data[x+1]-1, 1:length(ind_data)-1), length(ind_data)) #find consecutive data. Assume that every sweep starts with data
    for s in 1:length(ind0_edge)
        ind0_s =  s>1 ? ind0[ind0_edge[s-1]+1:ind0_edge[s]] : ind0[1:ind0_edge[1]]
        ind1_s =  s>1 ? ind1[ind1_edge[s-1]+1:ind1_edge[s]] : ind1[1:ind1_edge[1]]
        ind_data_s =  s>1 ? ind_data[ind_data_edge[s-1]+1:ind_data_edge[s]] : ind_data[1:ind_data_edge[1]]
        data_s = eval(Symbol(quad)).(data[1][qubit][ind_data_s])
        zero_cal = mean(eval(Symbol(quad)).(data[1][qubit][ind0_s]))
        one_cal = mean(eval(Symbol(quad)).(data[1][qubit][ind1_s]))
        scale_factor = -(one_cal - zero_cal)/2
        data_s = (data_s .- zero_cal)/scale_factor .+ 1
        push!(data_out, data_s)
    end
      return data_out
end

"""
    get_fidelity(shots_0, shots_1; nbins=51, showPlot=false)

Get readout fidelity from single-shot measurements
"""
function get_fidelity(shots_0, shots_1, nbins = 51, showPlot = false)
  bins = range(min(minimum(shots_0),minimum(shots_1)), stop=max(maximum(shots_0),length=maximum(shots_1)),nbins)
  hist_0 = kde(shots_0[:],bins)
  hist_1 = kde(shots_1[:],bins)
  cdf_0 = cumsum(hist_0.density)/cumsum(hist_0.density)[end]
  cdf_1 = cumsum(hist_1.density)/cumsum(hist_1.density)[end]
  if showPlot
    plot(hist_0.x, cdf_0)
    plot(hist_1.x, cdf_1)
    ylabel("Cumulative counts")
    xlabel("Homodyne voltage (a.u.)")
  end
  fidelity = 1-(1-maximum(abs.(cdf_0-cdf_1)))/2
  threshold = bins[argmax(abs.(cdf_0-cdf_1))]
  return fidelity, threshold
end

"""
    nanmean(vec)

Calculate the mean of vec, ignoring nan values.
"""
function nanmean(vec)
    return sum(.~isnan.(vec).*vec)/length(findall(.~isnan.(vec)))
end

"""
    get_extr_loc(M, dim; getmax=true)

Return index and value of max/min (getmax = True/False) in each column/row
(dim=1/2) of M
"""
function get_extr_loc(M, dim, getmax = true)
  #return index and value of max/min in each column (1)/row (2)
  if getmax
    _, aindx = findmax(M, dim)
  else
    _, aindx = findmin(M, dim)
  end
  msize=size(M);
  col = zeros(msize[-dim+3], 1);
  for i=1:msize[-dim+3]
    if dim == 1
      col[i], _ = ind2sub(msize, aindx[i])
    elseif dim == 2
      _, col[i] = ind2sub(msize, aindx[i])
    end
  end
  return convert(Array{Int8,2}, col), M[aindx]
end

"""
    get3pops(data)

get transmon populations after an experiment alternating Id, pi_01, pi_12 as tomography pulses
"""
function get3pops(data)
  calRepeat = 2;
  pvec3 = zeros(3,Int((length(data)-3*calRepeat)/3));
  pvec2 = zeros(3,Int((length(data)-3*calRepeat)/3));
  pvec2temp = zeros(2,Int((length(data)-3*calRepeat)/3));
  v0 = mean(data[end-3*calRepeat+1:end-2*calRepeat]);
  v1 = mean(data[end-2*calRepeat+1:end-calRepeat]);
  v2 = mean(data[end-calRepeat+1:end]);
  m0 = data[1:3:end]
  m1 = data[2:3:end]
  m2 = data[3:3:end]

  M = [v0 v1 v2; v1 v0 v2; v0 v2 v1];
  for c = 1:Int((length(data)-3*calRepeat)/3)
    pvec3[:,c] = M\[m0[c];m1[c];m2[c]]
  end

  #fixed p0+p1+p2 = 1
  M2 = [v0-v2 v1-v2; v1-v2 v0-v2];
  for c = 1:Int((length(data)-3*calRepeat)/3)
    pvec2temp[:,c] = M2\[m0[c]-v2;m1[c]-v2];
  end
  pvec2 = [pvec2temp; 1-sum(pvec2temp,1)];
  return (pvec3, pvec2)
end

"""
    get_feedback_fidelity(data::Int,
                          nqubits::Int,
                          nrounds::Int,
                          bit;
                          cal_repeats=2,
                          target_state=0)

Get fidelity for resetting 'bit' in a 'nqubits' register. Note that this
estimate assumes no measurement crosstalk.

# Arguments
  - `nrounds::Int`: number of feedback rounds
  - `bit`: 1, 2, 3... from LSB to MSB
  - `cal_repeats`: number of calibration points per computational state
  - `target_state`: 0/1 to target either qubit state
"""
function get_feedback_fidelity(data, nqubits, nrounds, bit; cal_repeat = 2, target_state = 0)
    zero_cal = mean(data[end-cal_repeat*(2^nqubits)+1:end-cal_repeat*(2^nqubits-1)])
    one_cal = mean(data[end-cal_repeat*(2^nqubits-2^(nqubits-bit))+1:end-cal_repeat*(2^nqubits-1-2^(nqubits-bit))])
    scale_factor = -(one_cal - zero_cal)/2;
    data = data[1:end-cal_repeat*2^nqubits]
    normdata = (data - zero_cal)./scale_factor + 1
    fidmat = zeros(2^nqubits, nrounds)
    for state_ind = 1:size(fidmat,1)
        for round_ind = 1:size(fidmat,2)
            fidmat[state_ind, round_ind] = (1+(-1)^target_state*normdata[(state_ind-1)*(nrounds+1)+round_ind+1])/2
        end
    end
    return fidmat
end

"""Unwrap angles by changing deltas between values to 2π complement.

Unwrap radian phase ϕ by changing absolute jumps greater than `discont` to
their 2π complement.

Args:
  ϕ: Input array.
  discont: Maximum discontinuity between values. Optional, default is π.

Returns:
  out: Output array of unwrapped phases.

See also:
  numpy.unwrap()
"""
function unwrap!(ϕ::Array{T}; discont=π) where T <: AbstractFloat
  Δ = diff(ϕ)
  Δmod = mod.(Δ .+ π, 2 * π) .- π
  Δmod[(Δmod .== -π) .& (Δ .> 0)] .= π
  ϕcorr = Δmod .- Δ
  ϕcorr[abs.(Δ) .< discont] .= 0
  return ϕ .+ vcat(0, cumsum(ϕcorr))
end

""" In-place version of unwrap. """
function unwrap(ϕ::Array{T}; discont=π) where T <: AbstractFloat
  return unwrap!(copy(ϕ), discont=discont)
end

function CR_hamiltonian(datapath,
                        filelist,
                        subdir=Dates.format(Dates.today(),"yymmdd"),
                        params_estimate = [5,5,0,0,0,1])
    IX = zeros(length(filelist))
    IY = zeros(length(filelist))
    IZ = zeros(length(filelist))
    ZX = zeros(length(filelist))
    ZY = zeros(length(filelist))
    ZZ = zeros(length(filelist))
    fit_curves = Dict()

    for (ct, filenum) in enumerate(filelist)
        data = load_data(datapath, filenum, datestr);

        norm_data = Qlab.cal_data(data[1], qubit="main")[1]
        data_len = Int(length(norm_data)/6)
        xvec = reshape(norm_data[1:3:end], data_len, 2)
        yvec = reshape(norm_data[2:3:end], data_len, 2)
        zvec = reshape(norm_data[3:3:end], data_len, 2);
        rvec = sqrt.(sum(xvec,2).^2+sum(yvec,2).^2+sum(zvec,2).^2)/2 # missing factor of 2 in PRA
        xpoints = data[2]["q1-main"][1]["points"][1:3:end-4][1:data_len]

        if ct>1
            params_estimate = tau_vec[1,:]
        end
        tau_vec = zeros(2,6)
        for ind=1:2
            result = optimize(x-> calculate_residueCR(x, xpoints, [xvec[:,ind],yvec[:,ind],zvec[:,ind]])[1], params_estimate)
            τ = Optim.minimizer(result)
            tau_vec[ind, :] = τ
            fit_curves[ind-1] = calculate_residueCR(τ, xpoints,  [xvec[:,ind],yvec[:,ind],zvec[:,ind]])[2]
        end
        IX[ct] = (tau_vec[1,1]+tau_vec[2,1])/2/(2*pi)
        IY[ct] = (tau_vec[1,2]+tau_vec[2,2])/2/(2*pi)
        IZ[ct] = (tau_vec[1,3]+tau_vec[2,3])/2/(2*pi)
        ZX[ct] = (tau_vec[1,1]-tau_vec[2,1])/2/(2*pi)
        ZY[ct] = (tau_vec[1,2]-tau_vec[2,2])/2/(2*pi)
        ZZ[ct] = (tau_vec[1,3]-tau_vec[2,3])/2/(2*pi)
        rfit = sqrt.(sum(hcat(fit_curves[0]+fit_curves[1]...).^2,2))/2
    end
end


text:
# Autogenerated wrapper script for RDKit_jll for aarch64-linux-gnu-cxx11
export librdkitcffi

using FreeType2_jll
using boost_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("RDKit")
JLLWrappers.@declare_library_product(librdkitcffi, "librdkitcffi.so.1")
function __init__()
    JLLWrappers.@generate_init_header(FreeType2_jll, boost_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        librdkitcffi,
        "lib/librdkitcffi.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
!isdir("src") && mkdir("src")
!isdir("src/rr") && cd("src") do
    run(`git clone https://github.com/mozilla/rr`)
end
!isdir("build/rr") && run(`mkdir -p build/rr`)
DEPS_DIR=pwd()
!isfile("build/rr/build.ninja") && cd("build/rr") do
    run(`cmake -G Ninja -DRR_BUILD_SHARED=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=$(DEPS_DIR)/usr ../../src/rr`)
end
cd("build/rr") do
    run(`ninja`)
    run(`ninja install`)
end


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

"""
    Ptr{T}

A memory address referring to data of type `T`.  However, there is no guarantee that the
memory is actually valid, or that it actually represents data of the specified type.
"""
Ptr

## converting pointers to an appropriate unsigned ##

"""
    C_NULL

The C null pointer constant, sometimes used when calling external code.
"""
const C_NULL = bitcast(Ptr{Cvoid}, 0)

# TODO: deprecate these conversions. C doesn't even allow them.

# pointer to integer
convert(::Type{T}, x::Ptr) where {T<:Integer} = T(UInt(x))

# integer to pointer
convert(::Type{Ptr{T}}, x::Union{Int,UInt}) where {T} = Ptr{T}(x)

# pointer to pointer
convert(::Type{Ptr{T}}, p::Ptr{T}) where {T} = p
convert(::Type{Ptr{T}}, p::Ptr) where {T} = bitcast(Ptr{T}, p)

# object to pointer (when used with ccall)

"""
    unsafe_convert(T, x)

Convert `x` to a C argument of type `T`
where the input `x` must be the return value of `cconvert(T, ...)`.

In cases where [`convert`](@ref) would need to take a Julia object
and turn it into a `Ptr`, this function should be used to define and perform
that conversion.

Be careful to ensure that a Julia reference to `x` exists as long as the result of this
function will be used. Accordingly, the argument `x` to this function should never be an
expression, only a variable name or field reference. For example, `x=a.b.c` is acceptable,
but `x=[a,b,c]` is not.

The `unsafe` prefix on this function indicates that using the result of this function after
the `x` argument to this function is no longer accessible to the program may cause undefined
behavior, including program corruption or segfaults, at any later time.

See also [`cconvert`](@ref)
"""
function unsafe_convert end

unsafe_convert(::Type{Ptr{UInt8}}, x::Symbol) = ccall(:jl_symbol_name, Ptr{UInt8}, (Any,), x)
unsafe_convert(::Type{Ptr{Int8}}, x::Symbol) = ccall(:jl_symbol_name, Ptr{Int8}, (Any,), x)
unsafe_convert(::Type{Ptr{UInt8}}, s::String) = convert(Ptr{UInt8}, pointer_from_objref(s)+sizeof(Int))
unsafe_convert(::Type{Ptr{Int8}}, s::String) = convert(Ptr{Int8}, pointer_from_objref(s)+sizeof(Int))
# convert strings to String etc. to pass as pointers
cconvert(::Type{Ptr{UInt8}}, s::AbstractString) = String(s)
cconvert(::Type{Ptr{Int8}}, s::AbstractString) = String(s)

unsafe_convert(::Type{Ptr{T}}, a::Array{T}) where {T} = ccall(:jl_array_ptr, Ptr{T}, (Any,), a)
unsafe_convert(::Type{Ptr{S}}, a::AbstractArray{T}) where {S,T} = convert(Ptr{S}, unsafe_convert(Ptr{T}, a))
unsafe_convert(::Type{Ptr{T}}, a::AbstractArray{T}) where {T} = error("conversion to pointer not defined for $(typeof(a))")

# unsafe pointer to array conversions
"""
    unsafe_wrap(Array, pointer::Ptr{T}, dims; own = false)

Wrap a Julia `Array` object around the data at the address given by `pointer`,
without making a copy.  The pointer element type `T` determines the array
element type. `dims` is either an integer (for a 1d array) or a tuple of the array dimensions.
`own` optionally specifies whether Julia should take ownership of the memory,
calling `free` on the pointer when the array is no longer referenced.

This function is labeled "unsafe" because it will crash if `pointer` is not
a valid memory address to data of the requested length.
"""
function unsafe_wrap(::Union{Type{Array},Type{Array{T}},Type{Array{T,N}}},
                     p::Ptr{T}, dims::NTuple{N,Int}; own::Bool = false) where {T,N}
    ccall(:jl_ptr_to_array, Array{T,N}, (Any, Ptr{Cvoid}, Any, Int32),
          Array{T,N}, p, dims, own)
end
function unsafe_wrap(::Union{Type{Array},Type{Array{T}},Type{Array{T,1}}},
                     p::Ptr{T}, d::Integer; own::Bool = false) where {T}
    ccall(:jl_ptr_to_array_1d, Array{T,1},
          (Any, Ptr{Cvoid}, Csize_t, Cint), Array{T,1}, p, d, own)
end
unsafe_wrap(Atype::Type, p::Ptr, dims::NTuple{N,<:Integer}; own::Bool = false) where {N} =
    unsafe_wrap(Atype, p, convert(Tuple{Vararg{Int}}, dims), own = own)

"""
    unsafe_load(p::Ptr{T}, i::Integer=1)

Load a value of type `T` from the address of the `i`th element (1-indexed) starting at `p`.
This is equivalent to the C expression `p[i-1]`.

The `unsafe` prefix on this function indicates that no validation is performed on the
pointer `p` to ensure that it is valid. Incorrect usage may segfault your program or return
garbage answers, in the same manner as C.
"""
unsafe_load(p::Ptr, i::Integer=1) = pointerref(p, Int(i), 1)

"""
    unsafe_store!(p::Ptr{T}, x, i::Integer=1)

Store a value of type `T` to the address of the `i`th element (1-indexed) starting at `p`.
This is equivalent to the C expression `p[i-1] = x`.

The `unsafe` prefix on this function indicates that no validation is performed on the
pointer `p` to ensure that it is valid. Incorrect usage may corrupt or segfault your
program, in the same manner as C.
"""
unsafe_store!(p::Ptr{Any}, @nospecialize(x), i::Integer=1) = pointerset(p, x, Int(i), 1)
unsafe_store!(p::Ptr{T}, x, i::Integer=1) where {T} = pointerset(p, convert(T,x), Int(i), 1)

# convert a raw Ptr to an object reference, and vice-versa
"""
    unsafe_pointer_to_objref(p::Ptr)

Convert a `Ptr` to an object reference. Assumes the pointer refers to a valid heap-allocated
Julia object. If this is not the case, undefined behavior results, hence this function is
considered "unsafe" and should be used with care.

See also: [`pointer_from_objref`](@ref).
"""
unsafe_pointer_to_objref(x::Ptr) = ccall(:jl_value_ptr, Any, (Ptr{Cvoid},), x)

"""
    pointer_from_objref(x)

Get the memory address of a Julia object as a `Ptr`. The existence of the resulting `Ptr`
will not protect the object from garbage collection, so you must ensure that the object
remains referenced for the whole time that the `Ptr` will be used.

This function may not be called on immutable objects, since they do not have
stable memory addresses.

See also: [`unsafe_pointer_to_objref`](@ref).
"""
function pointer_from_objref(@nospecialize(x))
    @_inline_meta
    typeof(x).mutable || error("pointer_from_objref cannot be used on immutable objects")
    ccall(:jl_value_ptr, Ptr{Cvoid}, (Any,), x)
end

## limited pointer arithmetic & comparison ##

isequal(x::Ptr, y::Ptr) = (x === y)
isless(x::Ptr{T}, y::Ptr{T}) where {T} = x < y

==(x::Ptr, y::Ptr) = UInt(x) == UInt(y)
<(x::Ptr,  y::Ptr) = UInt(x) < UInt(y)
-(x::Ptr,  y::Ptr) = UInt(x) - UInt(y)

+(x::Ptr, y::Integer) = oftype(x, add_ptr(UInt(x), (y % UInt) % UInt))
-(x::Ptr, y::Integer) = oftype(x, sub_ptr(UInt(x), (y % UInt) % UInt))
+(x::Integer, y::Ptr) = y + x

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
Base.signbit(x::Double64) = signbit(Hi(x))
signbits(x::Double64) = (signbit(Hi(x), signbit(Lo(x))))
Base.sign(x::Double64) = sign(Hi(x))
signs(x::Double64) = (sign(Hi(x), sign(Lo(x))))
Base.significand(x::Double64) = (significand(Hi(x)), significand(Lo(x)))
Base.exponent(x::Double64) = (exponent(Hi(x)), exponent(Lo(x)))

significant_bits(::Type{Float64})    =   53
significant_bits(::Type{Double64})   = 2*53
significant_bits(::Type{ComplexF64}) =   53
significant_bits(::Type{ComplexD64}) = 2*53

Base.frexp(x::Double64) = (frexp(Hi(x)), frexp(Lo(x)))

function Base.ldexp(x::Tuple{Tuple{Float64,Int64},Tuple{Float64,Int64}})
    hi = ldexp(x[1]...)
    lo = ldexp(x[2]...)
    return Double64((hi, lo))
end

# like frexp, returns an integer-valued significand
function ntexp(x::Float64)
    fr, ex = frexp(x)
    nt = Int128(ldexp(fr, 53))
    ex -= 53
    return nt, ex
end

#=
    num, xp2, den = Base.decompose(x)
    if isone(den)
       BigFloat(x) == ldexp(BigFloat(num), xp2)
    else
       BigFloat(x) == ldexp(BigFloat(num), xp2) / den
    end
=#
function Base.decompose(x::Double64)
    hi = Hi(x)
    isnan(hi)  && return 0, 0, 0
    isinf(hi)  && return sign(x), 0, 0
    iszero(hi) && return 0, 0, sign(x)

    nt_hi, ex_hi = ntexp(Hi(x))
    nt_lo, ex_lo = ntexp(Lo(x))
    num = (nt_hi << abs(ex_hi-ex_lo)) + nt_lo
    xp2 = ex_lo
    return num, xp2, 1
end

Base.copysign(x::Double64, y) = signbit(y) ? -abs(x) : abs(x)
Base.flipsign(x::Double64, y) = signbit(y) ? -x : x

Base.:(-)(x::Double64) = Double64((-Hi(x), -Lo(x)))
Base.abs(x::Double64) = signbit(x) ? -x : x
Base.abs2(x::Double64) = x*x
fastabs(x::Double64) = abs(Hi(x))
fastabs2(x::Double64) = abs2(Hi(x))

Base.:(-)(x::ComplexD64) = ComplexD64((-Hi(x), -Lo(x)))
fastabs(x::ComplexD64) = abs(Hi(x))
fastabs2(x::ComplexD64) = abs2(Hi(x))

fastabs(x::Float64) = abs(x)
fastabs(x::ComplexF64) = abs(x)
fastabs2(x::Float64) = abs2(x)
fastabs2(x::ComplexF64) = abs2(x)

function Base.ceil(x::Double64)
    xhi = Hi(x)
    hi = ceil(Hi(x))
    if hi === xhi
        t = ceil(Lo(x))
        hi = hi + t
        lo = t - (hi - xhi)
    else
        lo = 0.0
    end
    return Double64((hi, lo))
end

Base.ceil(::Type{T}, x::Double64) where {T<:Integer} = T(ceil(x))

function Base.floor(x::Double64)
    xhi = Hi(x)
    hi = floor(Hi(x))
    if hi === xhi
        t = floor(Lo(x))
        hi = hi + t
        lo = t - (hi - xhi)
    else
        lo = 0.0
    end
    return Double64((hi, lo))
end

Base.floor(::Type{T}, x::Double64) where {T<:Integer} = T(floor(x))

function Base.trunc(x::Double64)
    return signbit(x) ? ceil(x) : floor(x)
end

Base.trunc(::Type{T}, x::Double64) where {T<:Integer} = T(trunc(x))

function Base.div(x::Double64, y::Double64)
    (!isfinite(x) || isnan(y)) && return NaND64
    !isfinite(y) && return zero(Double64)
    return trunc(x / y)
end

function Base.fld(x::Double64, y::Double64)
    (!isfinite(x) || isnan(y)) && return NaND64
    !isfinite(y) && return zero(Double64)
    return floor(x / y)
end

function Base.cld(x::Double64, y::Double64)
    (!isfinite(x) || isnan(y)) && return NaND64
    !isfinite(y) && return zero(Double64)
    return ceil(x / y)
end

function Base.rem(x::Double64, y::Double64)
    (!isfinite(x) || isnan(y)) && return NaND64
    !isfinite(y) && return x
    return x - div(x,y) * y
end

function Base.mod(x::Double64, y::Double64)
    (!isfinite(x) || isnan(y)) && return NaND64
    !isfinite(y) && return x
    return x - fld(x,y) * y
end

function Base.divrem(x::Double64, y::Double64)
    (!isfinite(x) || isnan(y)) && return NaND64
    !isfinite(y) && return zero(Double64)
    dv = trunc(x / y)
    rm = x - dv * y
    return dv, rm
end

function Base.fldmod(x::Double64, y::Double64)
    (!isfinite(x) || isnan(y)) && return NaND64
    !isfinite(y) && return zero(Double64)
    fr = floor(x / y)
    md = x - fr * y
    return fr, md
end


text:
#Define some quantities and a very simple TargetModel for all tests below
rinitseed1 = 879
rnburnin1 = 2
rniter1 = 4
smhnprops1 = 1 #number of test proposals for standard mh
gmhnprops1 = 2 #number of test proposals for generalized mh
rtime1 = linspace(0.0,10.0,100)
rlower1 = [2.5]
rupper1 = [3.5]
rdefault1 = [3.0]
rcov1 = [0.01]
rparas1 = parameters([:a],rlower1,rupper1,rdefault1)
rdata1 = data(:array,rtime1,sin(3rtime1))
rnoise1 = noise(:gaussian,rcov1)
rmodel1 = model(:target,rparas1,rdata1,rnoise1,(p,t)->sin(p[1]*t),rtime1;name="RunnersTestModel")
rnparas1 = numparas(rmodel1)
rsampler1 = sampler(:mh,:normal,0.1,eye(1))
rtuner1 = tuner(:monitor,1)
rtuner2 = tuner(:scale,2,0.5,:logistic)


text:
"""
    Objective{F <: Function}
A structure to handle objectives in a solver.
```
struct Objective{F <: Function}
    name::String
    f::F
end
````
"""
struct Objective{F <: Function} <: FunctionContainer
    name::String
    f::F
end

"""
    Objective(F, o::Objective{F2}) where {F2 <: Function}
Constructor used in specializing a solver. Should never be called externally.
"""
function Objective(F, o::Objective{F2}) where {F2 <: Function}
    return Objective{F}(o.name, o.f)
end

"""
    objective(func, name)
Construct an objective with a function `func` that should be applied to a collection of variables.
"""
function objective(func, name)
    Objective(name, func)
end


text:
@testset "select_region" begin

@testset "select_region_null" begin
    chrom, start, stop = "01", 5000, 50000
    b = Bgen(example_16bits)
    @test b.idx !== nothing
    @test length(select_region(b, "02")) == 0
end

@testset "select_whole_chrom" begin
    chrom, start, stop = "01", 5000, 50000
    b = Bgen(example_16bits)
    lt = (x, y) -> isless((x.chrom, x.pos), (y.chrom, y.pos))
    variants = collect(select_region(b, chrom))
    for (x, y) in zip(sort(variants; lt=lt), sort(gen_data; lt=lt))
        @test (x.rsid, x.chrom, x.pos) == (y.rsid, y.chrom, y.pos)
    end
end

@testset "select_after_position" begin
    chrom, start, stop = "01", 5000, 50000
    b = Bgen(example_16bits)
    lt = (x, y) -> isless((x.chrom, x.pos), (y.chrom, y.pos))
    variants = collect(select_region(b, chrom; start=start))
    gen_data_f = filter(x -> x.pos >= start, gen_data)
    for (x, y) in zip(sort(variants; lt=lt), sort(gen_data_f; lt=lt))
        @test (x.rsid, x.chrom, x.pos) == (y.rsid, y.chrom, y.pos)
    end
end

@testset "select_in_region" begin
    chrom, start, stop = "01", 5000, 50000
    b = Bgen(example_16bits)
    lt = (x, y) -> isless((x.chrom, x.pos), (y.chrom, y.pos))
    variants = collect(select_region(b, chrom; start=start, stop=stop))
    gen_data_f = filter(x -> start <= x.pos <= stop, gen_data)
    for (x, y) in zip(sort(variants; lt=lt), sort(gen_data_f; lt=lt))
        @test (x.rsid, x.chrom, x.pos) == (y.rsid, y.chrom, y.pos)
    end
end
end


text:
# Adapted from https://discourse.julialang.org/t/expression-parser/41880/7
# code by Alan R. Rogers, Professor of Anthropology, University of Utah

import Nemo: fmpq_mpoly

function myeval(e::Union{Expr,Symbol,Number}, map::Dict{Symbol,fmpq_mpoly})
    try
        return _myeval(e, map)
    catch ex
        println("Can't parse \"$e\"")
        rethrow(ex)
    end
end

function _myeval(s::Symbol, map::Dict{Symbol,fmpq_mpoly})
    if haskey(map, s)
        return map[s]
    else
        @info "Can not find $s in $map while parsing.."
        throw(ParseException("$s"))
    end
end

# Numbers are converted to fmpq.
function _myeval(x::Number, map::Dict{Symbol,fmpq_mpoly})
    return fmpq(x)
end

# a helper definition for floats
function _myeval(x::Float64, map::Dict{Symbol,fmpq_mpoly})
    result = fmpq(0)

    # Getting the result from the string representation in order
    # to avoid approximations caused by the float representation
    s = string(x)
    denom = 1
    extra_num = 1
    if occursin(r"[eE]", s)
        s, exp = split(s, r"[eE]")
        if exp[1] == "-"
            denom = fmpq(10)^(-parse(Int, exp))
        else
            extra_num = fmpq(10)^(parse(Int, exp))
        end
    end
    frac = split(s, ".")
    if length(frac) == 1
        result = fmpq(parse(fmpz, s)) * extra_num // denom
    else
        result = fmpq(parse(fmpz, frac[1] * frac[2])) * extra_num // (denom * 10^(length(frac[2])))
    end

    # too verbose for now
    # @warn "a possibility of inexact float conversion" from=x to=result
    return result
end

# To parse an expression, convert the head to a singleton
# type, so that Julia can dispatch on that type.
function _myeval(e::Expr, map::Dict{Symbol,fmpq_mpoly})
    return _myeval(Val(e.head), e.args, map)
end

# Call the function named in args[1]
function _myeval(::Val{:call}, args, map::Dict{Symbol,fmpq_mpoly})
    return _myeval(Val(args[1]), args[2:end], map)
end

# Addition
function _myeval(::Val{:+}, args, map::Dict{Symbol,fmpq_mpoly})
    x = 0
    for arg ∈ args
        x += _myeval(arg, map)
    end
    return x
end

# Subtraction and negation
function _myeval(::Val{:-}, args, map::Dict{Symbol,fmpq_mpoly})
    len = length(args)
    if len == 1
        return -_myeval(args[1], map)
    else
        return _myeval(args[1], map) - _myeval(args[2], map)
    end
end

# Multiplication
function _myeval(::Val{:*}, args, map::Dict{Symbol,fmpq_mpoly})
    x = 1
    for arg ∈ args
        x *= _myeval(arg, map)
    end
    return x
end

# Division
function _myeval(::Val{:/}, args, map::Dict{Symbol,fmpq_mpoly})
    if typeof(_myeval(args[2], map)) <: fmpq_mpoly
        @warn "We can not parse polynomial fractions, sorry"
        throw(ParseException("Polynomial fractions are not supported"))
    end

    return _myeval(args[1], map) // _myeval(args[2], map)
end

# Absolute value
function _myeval(::Val{:abs}, args, map::Dict{Symbol,fmpq_mpoly})
    @warn "Absolute values of polynomials are ambigous"
    throw(ParseException("abs($(args[1])) encountered, exiting.."))
end

# Exponentiation
function _myeval(::Val{:^}, args, map::Dict{Symbol, fmpq_mpoly})
    if typeof( _myeval(args[2], map)) <: fmpq_mpoly
        @warn "We can not parse polynomial fractions, sorry"
        throw(ParseException("Polynomial fractions are not supported"))
    end

    if _myeval(args[2], map) < 0
        @warn "Negative exponent encountered while parsing"
        throw(ParseException("Negative exponents are not supported"))
    end

    return _myeval(args[1], map) ^ Int(numerator(_myeval(args[2], map)))
end


text:
__precompile__()

module SpinMonteCarlo

using Random
using Printf
using Markdown
using Statistics
using LinearAlgebra

using DataStructures

include("observables/MCObservables.jl")
include("API/api.jl")
include("model/model.jl")
include("lattice/Lattices.jl")
include("runMC.jl")
include("snapshot.jl")
end # of module


text:
const rules = rstrip(read("day07.txt", String), '\n')

function parse_line(line::AbstractString)
    color, contained = split(line, " bags contain ")
    pattern = r"(\d+) (.*?) bag"
    return color => [(parse(Int, m.captures[1]), m.captures[2]) for m in eachmatch(pattern, contained)]
end

BagsColors = Vector{Tuple{Int,String}}

function create_assoc(rules::AbstractString)
    return Dict{String,BagsColors}(parse_line.(line) for line in split(rules, '\n'))
end

function find_possibles(rules::AbstractString)
    to_check = ["shiny gold"]
    assoc = create_assoc(rules)
    possible = Set{String}()
    for color in to_check
        for (rule_color, rule_contained) in assoc
            if color in [r[2] for r in rule_contained]
                push!(possible, rule_color)
                if !(rule_color in to_check)
                    push!(to_check, rule_color)
                end
            end
        end
    end
    return possible
end

function part1(rules::AbstractString)
    possibles = find_possibles(rules)
    return length(possibles)
end

println("Part 1: $(part1(rules))")

function count_bags(assoc::Dict{String,BagsColors}, key::String)::Int
    if isempty(assoc[key])
        return 0
    end
    return sum(c + c * count_bags(assoc, new_key) for (c, new_key) in assoc[key])
end

function part2(rules::AbstractString)
    assoc = create_assoc(rules)
    return count_bags(assoc, "shiny gold")
end

println("Part 2: $(part2(rules))")


text:
# # Benchmarking flags for FFTW plans

# We test planning the FFTW with different flags.

# Here are the packages we need.

using FFTW
using Plots
using LinearAlgebra: mul!
using Test
using BenchmarkTools
using Random
using BenchmarkPlots
using StatsPlots

@info "Threads: $(FFTW.nthreads())"

# ## Preparing the benchmarks

# We prepare a suite of benchmarks, using [BenchmarkTools.BenchmarkGroup()](https://juliaci.github.io/BenchmarkTools.jl/stable/manual/#The-BenchmarkGroup-type).

# The best plan usually depends on the size of the array, which in this case is $N \times N$. In particular, it depens on the factorization of $N$. So, we benchmark different values of $N$.

# We choose the following values:

Ns = [2^6, 3^4, 2 * 3 * 5 * 7, 2^8, 2^2 * 3^2 * 5^2, 2^10, 2^11, 2^2 * 3 * 5^2 * 7]

# We must also prepare some other variables.

# List of flags for planning the FFTs:

flags = ["ESTIMATE", "MEASURE", "PATIENT", "EXHAUSTIVE"]
ext_flags = ["NO PLAN", flags...]
nothing

# Physical space:
L = 2π
κ₀ = 2π/L
nothing

# Excited modes for testing field:
rng = Xoshiro(123)
num_modes = 24
max_mode = 16

kxs = rand(rng, 1:max_mode, num_modes)
kys = rand(rng, 1:max_mode, num_modes)
ars = 10*randn(rng, num_modes)
ais = 10*randn(rng, num_modes)
nothing

# The suite of benchmarks:

suite = BenchmarkGroup()
for flag in ext_flags
    suite[flag] = BenchmarkGroup()
end

plan_stats = Dict{String, Dict{Int, Float64}}(flag => Dict() for flag in flags)

# Now we are ready to prepare the suite of benchmarks. Keep in mind that preparing the suite includes planning the transforms, and plans with a `PATIENT` or `EXAUSTIVE` flag take some time for large $N$.

for N in Ns
    x = y = (L/N):(L/N):L
    vort = sum(
        [
            2κ₀^2 * (kx^2 + ky^2) * (
                ar * cos.(κ₀ * (kx * one.(y) * x' + ky * y * one.(x)'))
                - ai * sin.(κ₀ * (kx * one.(y) * x' + ky * y * one.(x)'))
            )
            for (kx, ky, ar, ai) in zip(kxs, kys, ars, ais)
        ]
    )

    vort_hat = rfft(vort)

    flag = "NO PLAN"
    @info "N = $N; flag: $flag"
    suite[flag][N] = @benchmarkable rfft(w) setup = (w = copy($vort));

    for flag in flags
        @info "N = $N; flag: $flag"
        planed, pstats... = @timed plan_rfft(vort, flags = eval(Meta.parse("FFTW.$flag")))
        plan_stats[flag][N] = pstats.time

        suite[flag][N] = @benchmarkable mul!(w_hat, p, w) setup = (
            w_hat = copy($vort_hat);
            p = $planed;
            w = copy($vort)
        );
    end
end

suite

# Looking at `plan_stats`, we can see the time spent in planning the transforms.

for N in Ns
    @info "N = $N"
    for flag in flags
        @info "$flag: \t$(BenchmarkTools.prettytime(plan_stats[flag][N] * 1.0e+9))"
    end
end

# ## Running the benchmark

# This should take some time as well.

results, stats... = @timed run(suite)
nothing

# Here are the stats of the run:

stats

# Let's take a look at the results. The time shown is the minimum time of the trial runs for each benchmark.

results

# ## Analysis of the benchmark

# We start plotting the minimum and median times of the benchmark trials.

# First with the minimum time for low values of `N`.

plt = plot(
    title="Minimum times for different plans with vector fields of different sizes",
    xlabel = "N",
    ylabel = "time (ns)",
    xticks = Ns[1:4],
    rotation = 90,
    titlefont=10,
    legend=:topleft
)

for flag in ext_flags
    plot!(plt, Ns[1:4], N -> minimum(values(results[flag][N]).times),
    linestyle = :dash,
    markershape = :rect,
    label="$flag"
    )
end

plt

# Now, the minimum time for all values of `N`:

plt = plot(
    title="Minimum times for different plans with vector fields of different sizes",
    xlabel = "N",
    ylabel = "time (ns)",
    xticks = Ns,
    rotation = 90,
    titlefont=10,
    legend=:topleft
)

for flag in ext_flags
    plot!(plt, Ns, N -> minimum(values(results[flag][N]).times),
    linestyle = :dash,
    markershape = :rect,
    label="$flag"
    )
end

plt

# Next with the median time, starting with low values of `N`:

plt = plot(
    title="Median times for different plans with vector fields of different sizes",
    xlabel = "N",
    ylabel = "time (ms)",
    xticks = Ns[1:4],
    rotation = 90,
    titlefont=10,
    legend=:topleft
)

for flag in ext_flags
    plot!(plt, Ns[1:4], N -> median(values(results[flag][N]).times),
    linestyle = :solid,
    markershape = :circle,
    label="$flag"
    )
end
plt

# Now with all values of `N`:

plt = plot(
    title="Median times for different plans with vector fields of different sizes",
    xlabel = "N",
    ylabel = "time (ms)",
    xticks = Ns,
    rotation = 90,
    titlefont=10,
    legend=:topleft
)

for flag in ext_flags
    plot!(plt, Ns, N -> median(values(results[flag][N]).times)./ 1.0e+6,
    linestyle = :solid,
    markershape = :circle,
    label="$flag"
    )
end

plt

# Next we have a look at the set of trials, with violin plots.

plts = []

for N in Ns
    push!(
        plts,
        violin(
            [results[flag][N].times for flag in flags],
            title = "Trials with N = $N",
            titlefont = 12,
            xticks = (1:length(flags), string.(flags)),
            yaxis = "time (ns)",
            legend = nothing
        )
    )
end

if isodd(length(Ns))
    push!(plts, plot(border = :none))
end

plt = plot(plts..., layout = (div(length(plts), 2), 2), size = (800, 1200))

# We may look at the plot recipe build for the results of running a benchmark suite. But need to work on sorting the values of `N`.

# Let us take a final closer look at the median values.

for N in Ns
    @info "N = $N"
    for flag in flags
        @info "median time for flag $flag: $(BenchmarkTools.prettytime(median(values(results[flag][N]).times)))"
    end
end

# **Notice that `N` with powers of 2 are a little better.** The timings for the non-powers-of-2 are skewed up.

# We can have a better look at that by fitting the expected order of complexity. The discrete Fourier transform is of order $K\log(K)$, but this is a two-dimensional problem, so $K = N^2$.

# Let us fit a single plan, say `PATIENT`, with the median times.

flag = "PATIENT"

# We only fit the data for the powers of two.

twos = filter(N -> isinteger(log2(N)), Ns)

# Here is the Vandermonde matrix for the fit.

mat = [ones(length(twos)) [N^2 * log(N^2) for N in twos]]

# We get
a, b = mat \ [median(values(results[flag][N].times)) for N in twos]

# Now we plot the resulting fit agains the data

plt = plot(
    title="Median times for plan $flag with vector fields of different sizes",
    xlabel = "N (stretched out as N²)",
    ylabel = "time (ms)",
    xticks = (Ns.^2, string.(Ns)),
    rotation = 90,
    titlefont=10,
    legend=:topleft
)

plot!(plt, Ns.^2, map(N -> median(values(results[flag][isqrt(N)]).times)./ 1.0e+6, Ns.^2),
linestyle = :solid,
markershape = :circle,
label="$flag"
)

plot!(plt, Ns.^2, K -> (a + b * K * log(K) )./ 1.0e+6, linestyle = :dash, markershape = :square, label="N²log(N²) fit")

plt

# Now we take a closer look at the results for low `N`. Because the values are so discrepant for low and high `Ns` (they scale with $N^2$), instead of zooming into the result for low `N`, we do another fit, for these values. One option is to do a weighted least square fit, but that didn't turn out to be the best, so we we just restrict the fit to the low `N`.

# We look at the powers of two that have low `N`:

lowtwos = filter(N -> isinteger(log2(N)), Ns[1:4])

# And fit them with the Vandermonde matrix

lowmat = [ones(length(lowtwos)) [N^2 * log(N^2) for N in lowtwos]]

# We get
c, d = lowmat \ [median(values(results[flag][N].times)) for N in lowtwos]

# Now we can visualize the result.

plt = plot(
    title="Median times for plan $flag with vector fields of different sizes",
    xlabel = "N (stretched out as N²)",
    ylabel = "time (ms)",
    xticks = (Ns[1:4].^2, string.(Ns[1:4])),
    rotation = 90,
    titlefont=10,
    legend=:topleft
)

plot!(plt, Ns[1:4].^2, map(N -> median(values(results[flag][isqrt(N)]).times)./ 1.0e+6, Ns[1:4].^2),
linestyle = :solid,
markershape = :circle,
label="$flag"
)

plot!(plt, Ns[1:4].^2, K -> (c + d * K * log(K) )./ 1.0e+6, linestyle = :dash, markershape = :square, label="N²log(N²) fit")

plt

# Hmm, notice there are just two powers of two in this low range of `N`, so we could have just draw a straight line joint these two, but we leave it here in this generic way in case we change the values in the vector `Ns` of values of `N`.

# ## Conclusions

# Well, no plan is consistently worse, with all the others being closer together. In general, "PATIENT" and "EXHAUSTIVE" perform better, but planning them is costly, especially "EXHAUSTIVE", which doesn't seem to be worth it.

# Hence, for a quick plan, "MEASURE" seems to be a good choice. Otherwise, "PATIENT" seems like the best option.

# Moreover, powers of single primes, and in particular powers of 2, tend also to be faster, relatively speaking.

# Hence, prefer `N` with powers of 2, plan with `PATIENT`, and apply the plan with `mul!`.


text:
using RBM
using Base.Test

# test initialisation
rbm = rbm_create(11, 12, 13, 14, 0.1, 0.2, 0.3, 100, 15)
@test rbm.n          == 11
@test rbm.m          == 12
@test rbm.k          == 13
@test rbm.uditer     == 14
@test rbm.alpha      == 0.1
@test rbm.momentum   == 0.2
@test rbm.weightcost == 0.3
@test rbm.numepochs  == 100
@test rbm.batchsize  == 15

# test standard initialisation
rbm = rbm_create_with_standard_values(1, 2, 3)
@test rbm.n          == 1
@test rbm.m          == 2
@test rbm.k          == 3
@test rbm.uditer     == 10
@test rbm.alpha      == 0.1
@test rbm.momentum   == 0.5
@test rbm.weightcost == 0.001
@test rbm.numepochs  == 10000
@test rbm.batchsize  == 50

# test random initialisation functions
rbm = rbm_create_with_standard_values(10, 10, 10)
@test rbm.W == zeros(rbm.m, rbm.n)
@test rbm.V == zeros(rbm.m, rbm.k)
@test rbm.b == zeros(rbm.n)
@test rbm.c == zeros(rbm.m)

rbm_init_weights_random!(rbm)

@test rbm.W != zeros(rbm.m, rbm.n)
@test rbm.V != zeros(rbm.m, rbm.k)
@test size(rbm.W) == (rbm.m,rbm.n)
@test size(rbm.V) == (rbm.m,rbm.k)
@test rbm.b == zeros(rbm.n)
@test rbm.c == zeros(rbm.m)

rbm_init_output_bias_random!(rbm)
@test rbm.c != zeros(rbm.m)
@test size(rbm.c) == (rbm.m,)

rbm_init_hidden_bias_random!(rbm)
@test rbm.b != zeros(rbm.n)
@test size(rbm.b) == (rbm.n,)

#function rbm_init_hidden_bias_random!(rbm::RBM_t)

# test visible bias initialisation
r = [ [0, 0, 0, 0, 0, 0, 0, 0, 0, 0] [1, 0, 0, 0, 0, 0, 0, 0, 0, 0] [1, 1, 0, 0, 0, 0, 0, 0, 0, 0] [1, 1, 1, 0, 0, 0, 0, 0, 0, 0] [1, 1, 1, 1, 0, 0, 0, 0, 0, 0] [1, 1, 1, 1, 1, 0, 0, 0, 0, 0] [1, 1, 1, 1, 1, 1, 0, 0, 0, 0] [1, 1, 1, 1, 1, 1, 1, 0, 0, 0] [1, 1, 1, 1, 1, 1, 1, 1, 0, 0] [1, 1, 1, 1, 1, 1, 1, 1, 1, 0] [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]]

rbm = rbm_create_with_standard_values(11, 1, 1)
rbm_init_visible_bias!(rbm, r)

for i=1:11
  if i == 1
    s = 0.0
  elseif i == 11
    p = (10 - 1.0) / 10.0
    s = log2( p / ( 1 - p))
  else
    p = (i - 1.0) / 10.0
    s = log2( p / ( 1 - p))
  end
  @test rbm.b[i] == s
end

# test rbm copying
src   = rbm_create_with_standard_values(10, 10, 10)
src.W = rand(10,10)
src.V = rand(10,10)
src.b = rand(10)
src.c = rand(10)
dest  = rbm_copy(src)

@test dest.n  == src.n
@test dest.m  == src.m
@test dest.k  == src.k
@test dest.W  == src.W
@test dest.V  == src.V
@test dest.b  == src.b
@test dest.c  == src.c
@test dest.vW == src.vW
@test dest.vV == src.vV
@test dest.vb == src.vb
@test dest.vc == src.vc

dest.W[1,1] = dest.W[1,1] * 2
@test dest.W != src.W

# test weight scaling
rbm = rbm_create_with_standard_values(10, 10, 10)
rbm.W = rand(10,10) * 100
rbm.V = rand(10,10) * 100
cpy = rbm_copy(rbm)
rbm_rescale_weights!(rbm, 10.0)
mw = maximum(abs(rbm.W))
mv = maximum(abs(rbm.V))
if mw > mv
  m = mw
else
  m = mv
end
@test m == 10.0

mw = maximum(abs(cpy.W))
mv = maximum(abs(cpy.V))
if mw > mv
  mx = mw
else
  mx = mv
end

factor = m / mx

for i=1:size(rbm.W)[1]
  for j=1:size(rbm.W)[2]
    rbm.W[i,j]  = cpy.W[i,j] * factor
    rbm.W[i,j] != cpy.W[i,j]
  end
end

for i=1:size(rbm.V)[1]
  for j=1:size(rbm.V)[2]
    rbm.V[i,j]  = cpy.V[i,j] * factor
    rbm.V[i,j] != cpy.V[i,j]
  end
end

# test of assertion
#= rbm = rbm_create(11, 12, 13, 14, 0.1, 0.2, 100, 15) =#


text:
module FileStructs
export fileStructs

struct filePathAndModuleName
    filePath::String
    moduleName::String
    end

end # module


text:
using SimplTraits
@traitdef IsExpFam{D}
@traitimpl IsExpFam{D} <- isExpFam(D)

isExpFam(D) = false

@traitimpl IsExpFam{Normal{T} where T} = true

@traitfn f(d::D) where {D; IsExpFam{D}}


text:
using SparseArrays
using Random
using SpatialEcology
using Distances
using Test

@testset "ComMatrix" begin
    Random.seed!(1337)

    datb = sprand(Bool, 12,8,0.9)

    dati = sprand(8,13,0.9).*100
    dati = round.(Int, dati)

    datf = sprand(11,9,0.9)

    cmb = ComMatrix(datb)
    @test cmb isa ComMatrix{Bool}
    cmi = ComMatrix(dati,
        speciesnames = ["sp1", "sp2", "sp3", "sp4", "sp5", "sp6", "sp7", "sp8"],
        sitenames = [:a, :b, :c, :d, :e, :f, :g, :h, :i, :j, :k, :l, :m])
    @test cmi isa ComMatrix{<:Union{Int64, Int32}}
    cmf = ComMatrix(datf, sitecolumns = false)

    @test ComMatrix(Matrix(datf)).occurrences == ComMatrix(datf).occurrences

    @test occupancy(cmb) == [6, 6, 8, 7, 7, 8, 8, 8, 6, 6, 6, 7]
    @test occupancy(cmi) == [11, 13, 10, 12, 12, 13, 9, 13]
    @test occupancy(cmf) == [10, 10, 9, 10, 11, 8, 8, 10, 9]

    @test richness(cmb) == [12, 11, 10, 10, 11, 10, 12, 7]
    @test richness(cmi) == [7, 7, 8, 8, 5, 8, 8, 7, 8, 6, 6, 8, 7]
    @test richness(cmf) == [8, 8, 8, 9, 8, 8, 8, 7, 7, 7, 7]

    @test nsites(cmb) == 8
    @test nsites(cmi) == 13
    @test nsites(cmf) == 11

    @test nspecies(cmb) == 12
    @test nspecies(cmi) == 8
    @test nspecies(cmf) == 9

    @test occurring(cmb) == 1:12
    @test occurring(cmi) == 1:8
    @test occurring(cmf) == 1:9

    @test occupied(cmb) == 1:8
    @test occupied(cmi) == 1:13
    @test occupied(cmf) == 1:11

    @test noccurring(cmb) == 12
    @test noccurring(cmi) == 8
    @test noccurring(cmf) == 9

    @test noccupied(cmb) == 8
    @test noccupied(cmi) == 13
    @test noccupied(cmf) == 11

    @test getspecies(cmb, 4) == [true, true, true, true, true, true, true, false]
    @test getspecies(cmi, 2) == [48, 89, 34, 39, 29, 73, 49, 23, 62, 71, 60, 57, 19]
    gcm = getspecies(cmf, 6)
    @test gcm[3] == 0.37980641711782925
    @test gcm isa SubArray

    @test speciesnames(cmi)[8] == "sp8"
    @test speciesnames(cmf)[3] == "species3"

    @test sitenames(cmi)[4] == "d"
    @test sitenames(cmb)[5] == "site5"

    @test length(sitenames(cmf)) == nsites(cmf)
    @test length(sitenames(cmi)) == nsites(cmi)
    @test length(sitenames(cmb)) == nsites(cmb)

    @test length(speciesnames(cmf)) == nspecies(cmf)
    @test length(speciesnames(cmi)) == nspecies(cmi)
    @test length(speciesnames(cmb)) == nspecies(cmb)

    @test sitetotals(cmb) == [12, 11, 10, 10, 11, 10, 12, 7]
    @test sitetotals(cmb) == richness(cmb)
    @test sitetotals(cmi) == [334, 352, 466, 443, 139, 339, 314, 448, 414, 303, 473, 411, 275]
    @test length(sitetotals(cmf)) == nsites(cmf)
    @test sitetotals(cmf)[5] ≈ 4.52174342959108

    @test speciestotals(cmb) == occupancy(cmb)
    @test speciestotals(cmi) == [474, 653, 501, 577, 777, 521, 464, 744]
    @test length(speciestotals(cmf)) == nspecies(cmf)
    @test speciestotals(cmf)[8] ≈ 6.542348823677926

    @test size(cmb) == (12, 8)
    @test size(cmi, 1) == 8

    @test cooccurring(cmb, 1, 3) == [true, true, true, true, true, false, true, false]
    @test cooccurring(cmf, [8, 3]) == [false, true, false, true, true, true, false, true, true, true, true]

    dist = pairwise(BrayCurtis(), view(cmi, sites=1:2))
    @test round.(dist, digits=1) == [0.0 0.5; 0.5 0.0]
  # getindex and setindex are to do

end


text:
@testset verbose=true "indexing" begin
    import SphericalFunctions: WignerHsize, WignerHrange, WignerHindex
    import SphericalFunctions: WignerDsize, WignerDrange, WignerDindex
    import SphericalFunctions: Ysize, Yrange, Yindex, deduce_limits
    import SphericalFunctions: theta_phi

    ell_max = 16
    ell_max_slow = ell_max ÷ 2

    @testset "WignerHrange" begin
        r1(mp_max, ell_max) = hcat([
            [ell, mp, m] for ell in 0:ell_max
            for mp in -min(ell, mp_max):min(ell, mp_max)
            for m in abs(mp):ell
        ]...)'
        for ell_max in 0:ell_max
            a = WignerHrange(ell_max)  # Implicitly, mp_max=ell_max
            b = r1(ell_max, ell_max)
            @test a == b
            for mp_max in 0:ell_max
                a = WignerHrange(ell_max, mp_max)
                b = r1(mp_max, ell_max)
                @test a == b
            end
        end
    end

    @testset "WignerHsize" begin
        for i in -5:-1
            @test WignerHsize(i) == 0
            for j in -5:5
                @test WignerHsize(i, j) == 0
            end
        end
        for ell_max in 0:ell_max
            @test WignerHsize(ell_max) == size(WignerHrange(ell_max, ell_max), 1)
            for mp_max in ell_max
                @test WignerHsize(ell_max, mp_max) == size(WignerHrange(ell_max, mp_max), 1)
            end
        end
    end

    @testset "WignerHindex" begin
        function fold_H_indices(ell, mp, m)
            if m < -mp
                if m < mp
                    return [ell, -mp, -m]
                else
                    return [ell, -m, -mp]
                end
            else
                if m < mp
                    return [ell, m, mp]
                else
                    return [ell, mp, m]
                end
            end
        end

        for ell_max_i in 0:ell_max_slow
            r = WignerHrange(ell_max_i)
            for ell in 0:ell_max_i
                for mp in -ell:ell
                    for m in -ell:ell
                        i = WignerHindex(ell, mp, m)
                        @test r[i, :] == fold_H_indices(ell, mp, m)
                    end
                end
            end
            for mp_max in 0:ell_max_i
                r = WignerHrange(ell_max_i, mp_max)
                for ell in 0:ell_max_i
                    for mp in -ell:ell
                        for m in -ell:ell
                            if abs(mp) > mp_max && abs(m) > mp_max
                                continue
                            end
                            i = WignerHindex(ell, mp, m, mp_max)
                            @test r[i, :] == fold_H_indices(ell, mp, m)
                        end
                    end
                end
            end
        end
    end

    @testset "WignerDrange" begin
        function r2(ell_min, mp_max, ell_max)
            a = [
                [ℓ, mp, m]
                for ℓ in ell_min:ell_max
                for mp in -min(ℓ, mp_max):min(ℓ, mp_max)
                for m in -ℓ:ℓ
            ]
            collect(transpose(reshape(collect(Iterators.flatten(a)), (3, length(a)))))
        end

        for ell_max in 0:ell_max_slow÷2
            let ell_min = 0
                a = WignerDrange(ell_max)  # Implicitly, mp_max=ell_max
                b = r2(ell_min, ell_max, ell_max)
                @test a == b
                for mp_max in 0:ell_max
                    a = WignerDrange(ell_max, mp_max)
                    b = r2(ell_min, mp_max, ell_max)
                    @test a == b
                end
            end
        end
    end

    @testset "WignerDsize" begin
        for ell_max in 0:ell_max
            let ell_min = 0
                for mp_max in 0:ell_max
                    a = WignerDsize(ell_max, mp_max)
                    b = size(WignerDrange(ell_max, mp_max), 1)
                    @test a == b
                end
            end
        end
        for ell_max in 0:ell_max
            let ell_min = 0
                for mp_max in 0:ell_max
                    a = WignerDsize(ell_max, mp_max)
                    b = size(WignerDrange(ell_max, mp_max), 1)
                    @test a == b
                end
            end
        end
        for ell_max in 0:ell_max
            let ell_min = 0
                for mp_max in [ell_max]
                    a = WignerDsize(ell_max, mp_max)
                    # a = WignerDsize_ellmin(ell_min, ell_max)
                    b = size(WignerDrange(ell_max, mp_max), 1)
                    @test a == b
                end
            end
        end
        for ell_max in 0:ell_max
            let ell_min = 0
                for mp_max in [ell_max]
                    a = WignerDsize(ell_max, mp_max)
                    # a = WignerDsize(ell_max)
                    b = size(WignerDrange(ell_max, mp_max), 1)
                    @test a == b
                end
            end
        end
    end

    @testset "WignerDindex" begin
        for ellmax in 0:ell_max_slow
            r = WignerDrange(ellmax)
            for ell in 0:ellmax
                for mp in -ell:ell
                    for m in -ell:ell
                        i = WignerDindex(ell, mp, m)
                        @test r[i, :] == [ell, mp, m]
                    end
                end
            end
            let ell_min = 0
                r = WignerDrange(ellmax)
                for ell in ell_min:ellmax
                    for mp in -ell:ell
                        for m in -ell:ell
                            i = WignerDindex(ell, mp, m)
                            @test r[i, :] == [ell, mp, m]
                        end
                    end
                end
                for mp_max in 0:ellmax
                    r = WignerDrange(ellmax, mp_max)
                    for ell in ell_min:ellmax
                        for mp in -min(ell, mp_max):min(ell, mp_max)
                            for m in -ell:ell
                                i = WignerDindex(ell, mp, m, mp_max)
                                @test r[i, :] == [ell, mp, m]
                            end
                        end
                    end
                end
            end
        end
    end

    @testset "Yrange" begin
        function r3(ell_min, ell_max)
            a = [
                [ℓ, m]
                for ℓ in ell_min:ell_max
                for m in -ℓ:ℓ
            ]
            collect(transpose(reshape(collect(Iterators.flatten(a)), (2, length(a)))))
        end
        for ell_max in 0:ell_max
            for ell_min in 0:ell_max
                a = Yrange(ell_min, ell_max)
                b = r3(ell_min, ell_max)
                @test a == b
            end
        end
    end

    @testset "Ysize" begin
        for ell_max in 0:ell_max
            for ell_min in 0:ell_max
                a = Ysize(ell_min, ell_max)
                b = size(Yrange(ell_min, ell_max))[1]
                @test a == b
            end
        end
    end

    @testset "deduce_limits" begin
        for ℓmax in 0:4096
            for ℓmin in 0:min(2, ℓmax)
                deduced = deduce_limits(Ysize(ℓmin, ℓmax), nothing)
                @test deduced == (ℓmin, ℓmax)
            end
        end
        for ℓmax in 0:ell_max
            for ℓmin in [0]
                deduced = deduce_limits(Ysize(ℓmin, ℓmax))
                @test deduced == (ℓmin, ℓmax)
            end
        end
        for ℓmax in 0:ell_max
            for ℓmin in 0:ℓmax
                deduced = deduce_limits(Ysize(ℓmin, ℓmax), ℓmin)
                @test deduced == (ℓmin, ℓmax)
            end
        end
        for ℓmax in 1:ell_max
            prev_size = Ysize(ℓmax-1)
            this_size = Ysize(ℓmax)
            if abs(prev_size - this_size) > 1
                mid_size = (prev_size + this_size) ÷ 2
                @test_throws ErrorException deduce_limits(mid_size)
            end
        end
    end

    @testset "Yindex" begin
        for ell_max in 0:ell_max
            let ell_min = 0
                r = Yrange(ell_min, ell_max)
                for ell in ell_min:ell_max
                    for m in -ell:ell
                        i = Yindex(ell, m)
                        @test r[i, :] == [ell, m]
                    end
                end
                s = Yrange(ell_max)
                @test array_equal(r, s)
            end
            for ell_min in 0:ell_max
                r = Yrange(ell_min, ell_max)
                for ell in ell_min:ell_max
                    for m in -ell:ell
                        i = Yindex(ell, m, ell_min)
                        @test r[i, :] == [ell, m]
                    end
                end
            end
        end
    end

    @testset "phi_theta <-> theta_phi" for T in [Float64, Float32, BigFloat]
        for nθ in 2:(2ell_max+1)
            for nϕ in 2:(2ell_max+1)
                θϕ = theta_phi(nθ, nϕ, T)
                ϕθ = permutedims(phi_theta(nϕ, nθ, T), [2, 1, 3])[:, :, 2:-1:1]
                @test array_equal(θϕ, ϕθ)
            end
        end
    end

end


text:
using MAT, JuMP

function read_lp(name::String)
    lp_data = matopen("../netlib/$(name).mat")

    A = read(lp_data,"A")
    b = read(lp_data,"b")[:]
    c = read(lp_data,"c")[:]

    lbounds = nothing
    ubounds = nothing
    try
      lbounds = read(lp_data,"lbounds")[:]
      ubounds = read(lp_data,"ubounds")[:]
    catch(e)
      lbounds = zeros(length(c))
      ubounds = Inf * ones(length(c))
    end

    return A, b, c, lbounds, ubounds
end

function read_lp_into_Class_QP(name::String)
  A, b, c, lbounds, ubounds = read_lp(name)

  threshold = 1e20
  li = (lbounds .> -threshold)
  ui = (ubounds .< threshold)

  n = length(c)
  Ahat = [A; -A; speye(n)[li,:]; -speye(n)[ui,:]];
  bhat = [b; -b; lbounds[li]; -ubounds[ui]];

  return Class_QP(bhat, Ahat, c)
end
function read_lp_into_JuMP(name::String, perturb::Float64=0.0, pretend_is_nlp::Bool=false)
  A, b, c, lbounds, ubounds = read_lp(name)

  read_lp_into_JuMP(A, b, c, lbounds, ubounds, perturb, pretend_is_nlp)
end

function read_lp_into_JuMP(A, b, c, lbounds, ubounds, perturb::Float64=0.0, pretend_is_nlp::Bool=false)
  threshold = 1e20
  li = (lbounds .> -threshold)
  ui = (ubounds .< threshold)

  mod = Model()

  n = length(c)
  x = Array{JuMP.Variable,1}(n)
  for i = 1:n
    if true == li[i] && true == ui[i]
      x[i] = @variable(mod, lowerbound = lbounds[i] - perturb, upperbound = ubounds[i] + perturb)
    elseif true == li[i]
      x[i] = @variable(mod, lowerbound = lbounds[i] - perturb)
    elseif true == ui[i]
      x[i] = @variable(mod, upperbound = ubounds[i] + perturb)
    else
      x[i] = @variable(mod)
    end
  end

  if pretend_is_nlp
    @NLobjective(mod, Min, sum(c[i] * x[i] for i = 1:n))
  else
    @objective(mod, Min, sum(c[i] * x[i] for i = 1:n))
  end
  if length(b) > 0
    @constraint(mod, A * x .== b)
  end

  return mod
end

function run_IPOPT(A, b, c, lbounds, ubounds, max_iter, bound_relax_factor,mehrotra_algorithm)
    m = read_lp_into_JuMP(A, b, c, lbounds, ubounds, 0.0, true)
    nlp_raw = MathProgNLPModel(m)
    return run_IPOPT(nlp_raw, max_iter, bound_relax_factor,mehrotra_algorithm)
end

function run_IPOPT(nlp_raw, max_iter, bound_relax_factor, mehrotra_algorithm)
  t = 0
  hist = Array{ipopt_alg_history,1}()

  println(pd("t"), pd("comp"), pd("prm"), pd("dual"), pd("|x|"), pd("|y|"))
  A = jac(nlp_raw,nlp_raw.meta.x0)
  c = grad(nlp_raw,nlp_raw.meta.x0)

  ucon = nlp_raw.meta.ucon
  lcon = nlp_raw.meta.lcon
  #@show lbounds, ubounds
  infty = 1e30
  lbounds = max(-infty,nlp_raw.meta.lvar)
  ubounds = min(infty,nlp_raw.meta.uvar)
  #@show lbounds, ubounds

  for t = 1:max_iter

    solver = IpoptSolver(print_level=2, acceptable_iter=t, tol=1e-6, acceptable_tol=Inf, acceptable_compl_inf_tol=Inf, acceptable_constr_viol_tol=Inf, nlp_scaling_method="none", bound_relax_factor=bound_relax_factor, mehrotra_algorithm=mehrotra_algorithm)
    #setsolver(m, solver)
    #solve(m)
    #inner = m.internalModel.inner

    mp = NLPModels.NLPtoMPB(nlp_raw, solver)
    MathProgBase.optimize!(mp)
    inner = mp.inner

    hist_it = ipopt_alg_history()

    hist_it.t = t
    bounds_vio = max(0.0,maximum(lbounds - inner.x), max(0.0,maximum(inner.x -  ubounds )))
    a = cons(nlp_raw, inner.x)
    vio_A = max(maximum(a -  lcon), maximum(ucon - a))
    #norm(A * inner.x - b, Inf)

    hist_it.con_vio = max(vio_A, bounds_vio)
    #hist_it.primal_residual = max(norm(A * inner.x - b, Inf), max(0.0,maximum(lbounds - inner.x)), max(0.0,maximum(inner.x -  ubounds )))
    v1 = inner.mult_x_U - inner.mult_x_L
    v2 = A' * inner.mult_g
    v3 = v1 + v2

    hist_it.norm_grad_lag = norm(c + v3,Inf)
    hist_it.comp = max(maximum(abs(inner.x - lbounds) .* inner.mult_x_L),maximum(abs(inner.x - ubounds) .* inner.mult_x_U))
    hist_it.fval = dot(c, inner.x)
    hist_it.x_norm = norm(inner.x,Inf)
    hist_it.y_norm = max(norm(inner.mult_g,Inf),norm(inner.mult_x_U,Inf), norm(inner.mult_x_U,Inf))

    if t > 1 && hist_it.comp == hist[end].comp && hist_it.norm_grad_lag == hist[end].norm_grad_lag && hist_it.con_vio == hist[end].con_vio && hist_it.norm_grad_lag / (1.0 + norm(hist_it.y_norm,Inf)) < 1e-4 && hist[end].con_vio < 1e-4 && hist_it.comp < 1e-4
        break
    else
      println(pd(t), rd(hist_it.comp), rd(hist_it.con_vio), rd(hist_it.norm_grad_lag), rd(hist_it.x_norm), rd(hist_it.y_norm))
      push!(hist, hist_it)
    end
  end

  return hist
end

function load_netlib(num_nz::Int64)
    file_list = readdir("../netlib")
    problem_list = Array{String,1}()
    for file_name in file_list
        if file_name[(end-3):end] == ".mat" #&& filesize("netlib/$file_name") < file_size
          problem_name = file_name[1:(end-4)]
          println("loading...", problem_name)
          A, b, c, lbounds, ubounds = read_lp(problem_name)
          if nnz(A) <= num_nz
            push!(problem_list, problem_name)
          end
        end
    end

    return problem_list
end


text:
# This file is part of RLinearAlgebra.jl

module ProceduralTestLSVRDetermCycle

using Test, RLinearAlgebra, Random

Random.seed!(1010)

@testset "LSVR Deterministic Cyclic -- Procedural" begin

    # Verify appropriate super type
    @test supertype(LinSysVecRowDetermCyclic) == LinSysVecRowSampler

    # Test for appropriate cycling of rows
    @test let
        A = rand(10,3)
        b = rand(10)
        x = rand(3)

        cyc = LinSysVecRowDetermCyclic()

        flag = true

        for i = 11:20
            α, β = RLinearAlgebra.sample(cyc, A, b, x, i)
            flag = flag & (α == A[i-10,:])
            flag = flag & (β == b[i-10])
        end

        flag
    end

end

end


text:
using Oceananigans.Architectures: device
using Oceananigans.Operators: div_xyᶜᶜᵃ, Δzᵃᵃᶠ

"""
Compute the vertical velocity w by integrating the continuity equation from the bottom upwards

    `w^{n+1} = -∫ [∂/∂x (u^{n+1}) + ∂/∂y (v^{n+1})] dz`
"""
compute_w_from_continuity!(model) = compute_w_from_continuity!(model.velocities, model.architecture, model.grid)

function compute_w_from_continuity!(velocities, arch, grid)

    event = launch!(arch,
                    grid,
                    :xy,
                    _compute_w_from_continuity!,
                    velocities,
                    grid,
                    dependencies=Event(device(arch)))

    wait(device(arch), event)

    return nothing
end

@kernel function _compute_w_from_continuity!(U, grid)
    i, j = @index(Global, NTuple)
    # U.w[i, j, 1] = 0 is enforced via halo regions.
    @unroll for k in 2:grid.Nz+1
        @inbounds U.w[i, j, k] = U.w[i, j, k-1] - Δzᵃᵃᶠ(i, j, k, grid) * div_xyᶜᶜᵃ(i, j, k-1, grid, U.u, U.v)
    end
end


text:
module DrugImpurityAnalysis

using Base
using Plots
using CSV
using DataFrames
using LightXML
using JSON3
using Dates
using Statistics

export create_impurity_profiles_batch, create_impurity_profile, analyse_spectrum
export compound_percentage_occurrence, compound_trends, find_similar_samples
export mz_integrals_to_csv, manual_integration, visualize_peak_range, visualize_peak_determination_steps, sort_spectra, retrieve_spectrum
export batch_import, filter_XIC

include("helpers.jl")
include("analysis.jl")
include("manual_inspection.jl")
include("MS_Import_light.jl")

#------ Global Variables ------#
const MAIN_RT_DEVIATION = 1 # Absolute maximum shift of cocaine and IS peak
const MAX_RT_MODIFIER_DIFFERENCE = 0.015 # Difference RT_modifier cocaine and IS must be lower or equal

const MAX_MZ_DEVIATION = 0.5 # Maximum deviation from given mz value in XIC spectrum
const MAX_RT_SHIFT = 0.05 # Maximum amount that RT can shift from determined RT

# 40 scans is 0.1 minute
const MAX_SCANS_PEAK_SEARCH = 16 # Used for differentiating peak from noise
const MAX_SCANS_PEAK_LEFT = 40 # Bounds for peak integration (left)
const MAX_SCANS_PEAK_RIGHT = 80 # Bounds for peak integration (right)

# If fraction of median/mean crossings is above either value, peak attributed to noise
const NOISE_MEDIAN_CROSSINGS_FRACTION = 0.50
const NOISE_MEAN_CROSSINGS_FRACTION = 0.25

# Amount of scans for determining baseline
const BASELINE_SCANS = 20

# Used for determining overlap
const OVERLAP_CONSECUTIVE_BELOW_MEDIAN = 2
const OVERLAP_CONSECUTIVE_ABOVE_MEDIAN = 5

const COMPOUNDS_CSV_LOCATION = dirname(@__DIR__)
const SETTINGS_JSON_LOCATION = dirname(@__DIR__)
#------ Global Variables End------#

"""
	create_impurity_profiles_batch(pathin; pathout=pathin, start_at=1, append=false)
Creates impurity profiles from all folders in pathin and writes to one csv file at pathout
start_at: which folder to start (ordered by folder name), skips all folders prior
append: set to true if appending to already existing csv file
"""
function create_impurity_profiles_batch(pathin::String; pathout::String=pathin, start_at::Int=1, append::Bool=false)
	csvout = joinpath(pathout, "impurity_profile.csv")
	subdirs = [dir for dir in readdir(pathin) if isdir(joinpath(pathin, dir)) == true]

	for i=start_at:length(subdirs)
		impurity_profile = create_impurity_profile(joinpath(pathin, subdirs[i]))
		CSV.write(csvout, impurity_profile, append=append)
		println("\n---Processed directory $(subdirs[i]) ($i of $(length(subdirs)))---\n")
		append = true
	end
end

"""
	create_impurity_profile(pathin; csvout=nothing)
Creates and returns impurity profile from all samples in directory pathin
Writes impurity profile to csvout if given
Set use_filters to false to ignore filters and include all spectra in impurity profile
"""
function create_impurity_profile(pathin; csvout=nothing, use_filters=true)
	# Read settings.json
	json_string = read(joinpath(SETTINGS_JSON_LOCATION, "settings.json"), String)
	settings_json = JSON3.read(json_string)
	main_compound_name = settings_json[:main_settings]["main_compound"]
	IS_name = settings_json[:main_settings]["internal_standard"]
	main_IS_min_ratio = settings_json[:main_settings]["main/IS min ratio"]

	# Import spectra
	spectra, metadata_headers = batch_import(pathin, settings_json, use_filters=use_filters)

	# Import RT and mz info of valid compounds into DataFrame
	compounds_csv = CSV.read(joinpath(COMPOUNDS_CSV_LOCATION, "compounds.csv"), DataFrame)
	filter!(row -> !(any(ismissing, (row.RT, row.mz)) || any((row.RT, row.mz) .== 0)), compounds_csv)
	main_compound = filter(row -> row.compound == main_compound_name, compounds_csv)[1, :]
	internal_standard = filter(row -> row.compound == IS_name, compounds_csv)[1, :]

	# Create DataFrame for storing impurity profile (output)
	exclude_list = settings_json[:exclude_from_impurity_profile]
	compounds_in_profile = filter(row -> !(row.compound in exclude_list), compounds_csv).compound
	impurity_profiles = DataFrame()
	for header in metadata_headers
		insertcols!(impurity_profiles, Symbol(header) => String[])
	end
	for compound_name in compounds_in_profile
		insertcols!(impurity_profiles, Symbol(compound_name) => Float32[])
	end

	# Analyse all spectra
	for i=1:length(spectra)
		@info "Analysing spectrum $i..."
		spectrum = spectra[i]["MS1"]

		compound_mz_integrals = analyse_spectrum(spectrum, compounds_csv, main_compound, internal_standard)
		# Create impurity profile and add to DataFrame
		sample_metadata, sample_profile = calculate_impurity_profile_values(
							 spectrum, compound_mz_integrals, compounds_csv,
							 compounds_in_profile, metadata_headers,
							 main_compound, IS_name, main_IS_min_ratio)
		if sample_profile[1] > 0 #TODO excludes samples not containing cocaine, should be done more dynamically
			push!(impurity_profiles, append!(sample_metadata, sample_profile))
		end
	end

	# Write to csv if csvout is given
	if !isnothing(csvout)
		CSV.write(csvout, impurity_profiles)
	end

	return impurity_profiles
end

"""
	analyse_spectrum(spectrum, compounds_csv, main_compound, internal_standard)
Analyses all compounds in spectrum and returns dictionary with mz values for each compound
"""
function analyse_spectrum(spectrum, compounds_to_analyse, main_compound, internal_standard)
	RT_modifier = determine_RT_modifier(spectrum, main_compound, internal_standard)

	# Create Dict for storing integral for each compound and mz value
	compound_mz_integrals = Dict(Symbol(compound) => Dict{Symbol, Int}() for compound in compounds_to_analyse.compound)

	# No cocaine peak found
	if RT_modifier == -1
		return 0
	end

	# Iterate over compounds in compounds.csv
	for row in eachrow(compounds_to_analyse)
		RT = row.RT * RT_modifier
		mz_vals, overlap_RT_vals = parse_data(row)
		overlap_RT_vals .*= RT_modifier

		# Iterate over each mz value per compound
		for mz in mz_vals
			spectrum_XIC = filter_XIC(spectrum, mz)
			peak_exists, max_scan = search_peak(spectrum_XIC, RT, spectrum)

			if peak_exists == false
				compound_mz_integrals[Symbol(row.compound)][Symbol(mz)] = 0
				continue
			end

			# Determine baseline and bounds
			pre_bounds = max_scan - MAX_SCANS_PEAK_LEFT : max_scan + MAX_SCANS_PEAK_RIGHT
			overlap_max_left, overlap_max_right = determine_overlap(spectrum_XIC, max_scan, overlap_RT_vals, spectrum)
			baseline = determine_baseline(spectrum_XIC, max_scan, overlap_max_left, overlap_max_right, pre_bounds)
			peak_bounds = determine_bounds(spectrum_XIC, max_scan, overlap_max_left, overlap_max_right, baseline, pre_bounds)

			# Calculate and store integral
			peak_integral = integrate_peak(spectrum_XIC, peak_bounds, baseline, pre_bounds)
			compound_mz_integrals[Symbol(row.compound)][Symbol(mz)] = peak_integral
		end
	end

	return compound_mz_integrals
end

"""
	determine_RT_modifier(spectrum, main_compound, internal_standard)
main_compound and internal_standard should be input as DataFrameRow
Determines retention time shift of cocaine and internal standard
and returns mean shift
"""
function determine_RT_modifier(spectrum, main_compound, internal_standard)

	RT_modifiers = Vector{Float32}()
	for compound in (main_compound, internal_standard)
		# Read predicted RT and highest mz value
		RT_predicted = compound.RT
		highest_mz_val = maximum(parse_data(compound)[1])

		# Search for maximum around predicted RT
		RT_range = (RT_predicted - MAIN_RT_DEVIATION, RT_predicted + MAIN_RT_DEVIATION)
		RT_index_range = RT_to_scans(spectrum, RT_range)
		spectrum_XIC_part = filter_XIC(spectrum, highest_mz_val)[RT_index_range[1] : RT_index_range[2]]

		# No peak at all
		if maximum(spectrum_XIC_part) <= 0
			@info "No internal standard and/or cocaine found"
			return -1
		end

		max_scan = findmax(spectrum_XIC_part)[2] + RT_index_range[1] - 1
		RT_actual = spectrum["Rt"][max_scan]
		push!(RT_modifiers, round(RT_actual / RT_predicted, digits=3))
	end

	if abs(RT_modifiers[1] - RT_modifiers[2]) > MAX_RT_MODIFIER_DIFFERENCE
		@info "No internal standard and/or cocaine found"
		return -1
	end

	return round(mean(RT_modifiers), digits=3)
end

"""
	search_peak(spectrum_XIC, RT, spectrum)
Searches for peak in ion-extracted chromatogram
Returns true or false, and scan number of maximum if true
"""
function search_peak(spectrum_XIC, RT, spectrum, left_scan=0, right_scan=0)
	# Range to search peak in
	RT_range = (RT - MAX_RT_SHIFT, RT + MAX_RT_SHIFT)
	scan_range = collect(RT_to_scans(spectrum, RT_range))

	# If scans pre-defined
	scan_range[1] = left_scan == 0 ? scan_range[1] : left_scan
	scan_range[2] = right_scan == 0 ? scan_range[2] : right_scan

	max_intensity, max_scan_number = findmax(spectrum_XIC[scan_range[1] : scan_range[2]])
	max_scan_number += scan_range[1] - 1

	# No peak at all
	if max_intensity <= 0
		return false, -1
	end

	# Bounds used for determining whether peak is noise or not
	left_scan = left_scan == 0 ? max_scan_number - round(Int, MAX_SCANS_PEAK_SEARCH / 2) : left_scan
	right_scan = right_scan == 0 ? max_scan_number + round(Int, MAX_SCANS_PEAK_SEARCH / 2) : right_scan

	# Count number of median crosses
	spectrum_part_median = median(spectrum_XIC[left_scan:right_scan])
	median_crossings_count = 0
	for i in (left_scan + 1):right_scan
		if (spectrum_XIC[i] - spectrum_part_median) * (spectrum_XIC[i - 1] - spectrum_part_median) <= 0
			median_crossings_count += 1
		end
	end

	# Count number of mean crosses
	spectrum_part_mean = mean(spectrum_XIC[left_scan:right_scan])
	mean_crossings_count = 0
	for i in (left_scan + 1):right_scan
		if (spectrum_XIC[i] - spectrum_part_mean) * (spectrum_XIC[i - 1] - spectrum_part_mean) <= 0
			mean_crossings_count += 1
		end
	end

	# More crossings than allowed, noise instead of peak
	if median_crossings_count / length(spectrum_XIC[left_scan:right_scan]) > NOISE_MEDIAN_CROSSINGS_FRACTION ||
	   		mean_crossings_count / length(spectrum_XIC[left_scan:right_scan]) > NOISE_MEAN_CROSSINGS_FRACTION
		return false, -1
	else
		return true, max_scan_number
	end
end

"""
	determine_overlap(spectrum_XIC, max_scan, RT_overlap_vals, spectrum)
Determines scan numbers of maxima of overlapping peaks
For side at which no overlap is defined, takes a simple and
non-thorough approach for determining whether overlap is present
"""
function determine_overlap(spectrum_XIC, max_scan, RT_overlap_vals, spectrum)
	overlap_max_scans = zeros(Int, 2)
	for (i, RT_overlap) in enumerate(RT_overlap_vals)

		left_scan, right_scan = 0, 0
		# Overlap not defined
		if RT_overlap == 0
			# Determine iteration sequence
			direction  = i == 1 ? -1 : 1
			scan_upper_bound = i == 1 ? -MAX_SCANS_PEAK_LEFT : MAX_SCANS_PEAK_RIGHT
			spectrum_scan_range = max_scan:direction:(scan_upper_bound + max_scan)

			# Determine mean of range and initialize variables
			range_mean = mean(spectrum_XIC[spectrum_scan_range])
			below_consecutive, above_consecutive = 0, 0
			reached_below = false
			overlap_max_scan = 0

			# Iterate through sequence and search for peak separate from main peak
			for scan=spectrum_scan_range
				if spectrum_XIC[scan] < range_mean
					below_consecutive += 1
					above_consecutive = 0
				else
					below_consecutive = 0
					above_consecutive += 1
				end

				# Consecutively below median
				if below_consecutive >= OVERLAP_CONSECUTIVE_BELOW_MEDIAN
					reached_below = true
				# Previously consecutively below median and now consecutively above median
				elseif reached_below && above_consecutive >= OVERLAP_CONSECUTIVE_ABOVE_MEDIAN
					overlap_range = (scan - OVERLAP_CONSECUTIVE_ABOVE_MEDIAN * direction):direction:spectrum_scan_range[end]

					# Bounds for search_peak
					left_scan = direction == -1 ? overlap_range[end] : overlap_range[begin]
					right_scan = direction == -1 ? overlap_range[begin] : overlap_range[end]
					break
				end
			end

			# Set RT_overlap to 0 if not found, else to RT of maximum
			RT_overlap = overlap_max_scan == 0 ? 0 : spectrum["Rt"][overlap_max_scan]
		else
			overlap_scan = RT_to_scans(spectrum, RT_overlap)[1]
			scan_between = round(Int, (overlap_scan + max_scan) / 2)
			left_scan = i == 1 ? 0 : scan_between
			right_scan = i == 1 ? scan_between : 0
		end

		if RT_overlap > 0
			# Final check whether peak exists
			overlap_max_scans[i] = search_peak(spectrum_XIC, RT_overlap, spectrum, left_scan, right_scan)[2]
		end
	end

	return overlap_max_scans
end

"""
	determine_baseline(spectrum_XIC, max_scan, overlap_max_left, overlap_max_right, bounds)
Determines baseline based on range left of peak
"""
function determine_baseline(spectrum_XIC, max_scan, overlap_max_left, overlap_max_right, bounds)
	#TODO take section on right side if overlap on left
	if overlap_max_left > 0
		@warn "Overlap on left side" overlap_max_left, overlap_max_right
	end
	# Noise end as minimum closest to peak
	noise_end = max_scan - findmin(spectrum_XIC[max_scan:-1:max_scan - 30])[2] + 1
	spectrum_part = spectrum_XIC[noise_end - BASELINE_SCANS:noise_end]

	# Calculate median and set median as baseline
	spectrum_part_median = round(Int, median(spectrum_part))
	baseline = fill(spectrum_part_median, length(bounds))

	return baseline
end

"""
Determine bounds of peak
Set bounds to baseline crossing closest to peak, or minimum if it does not cross
"""
function determine_bounds(spectrum_XIC, max_scan, overlap_max_left, overlap_max_right, baseline, pre_bounds)
	# Set bound to max of overlapping peak if it exists
	pre_left_bound = overlap_max_left > 0 && overlap_max_left > pre_bounds[begin] ? overlap_max_left : pre_bounds[begin]
	pre_right_bound = overlap_max_right > 0 && overlap_max_right < pre_bounds[end] ? overlap_max_right : pre_bounds[end]

	spectrum_part = spectrum_XIC[pre_left_bound:pre_right_bound]
	part_max_scan = max_scan - pre_left_bound + 1

	# Subtract baseline from spectrum, negative values set to zero
	spectrum_part = map((x_spectr, x_base) -> (x_spectr < x_base ? 0 : x_spectr - x_base), spectrum_part, baseline)

	# Set bounds to zero closest to maximum or minimum if no zero value
	left_bound  = max_scan - findmin(spectrum_part[part_max_scan:-1:begin])[2] + 1
	right_bound = max_scan + findmin(spectrum_part[part_max_scan:end])[2] - 1

	return left_bound:right_bound
end

"""
	integrate_peak(spectrum_XIC, bounds, baseline, pre_bounds)
Integrates peak within bounds with baseline subtracted
"""
function integrate_peak(spectrum_XIC, bounds, baseline, pre_bounds)
	spectrum_part = spectrum_XIC[pre_bounds]
	# spectrum_part .-= baseline
	bounds = bounds .- (pre_bounds[begin] - 1)

	# Subtract baseline from spectrum, negative values set to zero
	spectrum_part = map((x_spectr, x_base) -> (x_spectr < x_base ? 0 : x_spectr - x_base), spectrum_part, baseline)

	# Calculate integral within bounds
	peak_integral = sum(spectrum_part[bounds])

	return peak_integral
end

"""
	calculate_impurity_profile_values(spectrum, compound_mz_integrals, compounds_csv,
									  compounds_in_profile, metadata_headers,
									  main_compound, IS_name, main_IS_min_ratio)
									  Returns metadata, impurity_profile_values
Calculates values to use in impurity profile from analyzed spectra and retrieves
metadata to store.
"""
function calculate_impurity_profile_values(spectrum, compound_mz_integrals, compounds_csv,
										   compounds_in_profile, metadata_headers,
										   main_compound, IS_name, main_IS_min_ratio)

	# Retrieve metadata
	metadata = Vector{Any}([spectrum[header] for header in metadata_headers])

	# No cocaine present
	if compound_mz_integrals == 0
		return metadata, zeros(length(compounds_in_profile))
	end

	# Determine integral of highest mz value of main compound
	main_highest_mz_value = maximum(parse_data(main_compound)[1])
	main_highest_mz_integral = compound_mz_integrals[Symbol(main_compound.compound)][Symbol(main_highest_mz_value)]

	# Sum integrals of all mz values of internal standard
	IS_integral = sum(values(compound_mz_integrals[Symbol(IS_name)]))

	impurity_profile_values = zeros(length(compounds_in_profile))
	# Ratio too low, no significant amount of main compound present
	if main_highest_mz_integral / IS_integral < main_IS_min_ratio || main_highest_mz_integral == 0
		return metadata, zeros(length(compounds_in_profile))
	elseif IS_integral == 0 && main_highest_mz_integral != 0
		@warn "main compound detected, but internal standard not detected"
		return metadata, zeros(length(compounds_in_profile))
	end

	# Sum all integrals of cocaine
	main_integral_sum = sum(values(compound_mz_integrals[Symbol(main_compound.compound)]))

	# Calculate percent ratio for all compounds
	for (i, compound_name) in enumerate(compounds_in_profile)
		mz_dict = compound_mz_integrals[Symbol(compound_name)]

		# Check if first mass has no intensity
		first_mz = parse_data(filter(row -> row.compound == compound_name, compounds_csv)[1, :])[1][1]
		if mz_dict[Symbol(first_mz)] <= 0
			impurity_profile_values[i] = 0
			continue
		end

		# Calculate percent ratio
		integral_sum = sum(values(mz_dict))
		percent_ratio = round(integral_sum / main_integral_sum * 100, digits=2)
		impurity_profile_values[i] = percent_ratio
	end

	return metadata, impurity_profile_values
end

end # Module


text:
pages = [
        "Home" => "index.md",
        "Tutorials" => tutorials,
        "Problems" => "problems.md",
        "Solvers" => Any[
          "solvers/common.md",
          "solvers/koopman.md",
          "solvers/optimization.md",
          "solvers/symbolic_regression.md"
          ],
        "Basis" => "basis.md",
        "Solutions" => "solutions.md",
        "Utilities" => "utils.md",
        "Contributing" => "contributions.md",
        "Citing" => "citations.md"
        ]


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function readLAMMPS(filename,datamode)

if datamode==0 # binary
    tmp = reinterpret(Float64,read(filename));
else # text
    tmp = Main.DelimitedFiles.readdlm(filename);
end

config = CONFIGStruct();

config.dim = 3;
config.ncq = 0;
config.ncv = 0;
config.nce = 1;
config.ncf = 3;
config.ncx = 3;

k = 1;
sz1, sz2 = size(tmp);
while (k < sz1)
    n = tmp[k,1]; # number of atoms
    if k ==1
        config.natom = reshape([n],1,1);
    else
        config.natom = [config.natom n];
    end
    b = transpose(reshape(tmp[(k+1):(k + n),:], n, 7));
    if k==1
        config.t = reshape(b[1,:],1,n);
        config.x = b[2:4,:];
        config.f = b[5:7,:];
    else
        config.t = [config.t reshape(b[1,:],1,n)];
        config.x = [config.x b[2:4,:]];
        config.f = [config.f b[5:7,:]];
    end
    b = tmp[(k+n+1):(k+n+2),:];
    if k == 1
        config.e = reshape([b[1,1]],1,1);
        config.a = reshape([b[2,1] 0 0],3,1);
        config.b = reshape([0 b[2,2] 0],3,1);
        config.c = reshape([0 0 b[2,3]],3,1);
    else
        config.e = [config.e reshape([b[1,1]],1,1)];
        config.a = [config.a [b[2,1] 0 0]'];
        config.b = [config.b [0 b[2,2] 0]'];
        config.c = [config.c [0 0 b[2,3]]'];
    end
    k = k + n + 3;
end

config.q = reshape([],0,2);
config.v = reshape([],0,2);

return config

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    Atom, QueryAtom,
    Bond, QueryBond

abstract type Atom <: AbstractNode end
abstract type QueryAtom <: AbstractNode end

abstract type Bond <: UndirectedEdge end
abstract type QueryBond <: UndirectedEdge end

# Aliases

getatom(mol, i) = nodeattr(mol, i)
getbond(mol, i) = edgeattr(mol, i)
getbond(mol, u, v) = edgeattr(mol, u, v)
hasbond(mol, u, v) = hasedge(mol, u, v)

setatom!(mol, i, attr) = setnodeattr!(mol, i, attr)
setbond!(mol, i, attr) = setedgeattr!(mol, i, attr)
addatom!(mol, attr) = addnode!(mol, attr)
addbond!(mol, u, v, attr) = addedge!(mol, u, v, attr)

atomcount(mol) = nodecount(mol)
bondcount(mol) = edgecount(mol)


text:
using PolyChaos, Test, LinearAlgebra

# values taken from https://en.wikipedia.org/wiki/Hermite_polynomials#Definition
c =    [
       [-0.0] ,
       [-1.0, -0.0]   ,
       [0.0, -3.0, -0.0]       ,
       [3.0, 0.0, -6.0, -0.0]     ,
       [-0.0, 15.0, 0.0, -10.0, -0.0]  ,
       [-15.0, -0.0, 45.0, 0.0, -15.0, -0.0]    ,
       [0.0, -105.0, -0.0, 105.0, 0.0, -21.0, -0.0]  ,
       [105.0, 0.0, -420.0, -0.0, 210.0, 0.0, -28.0, -0.0] ,
       [-0.0, 945.0, 0.0, -1260.0, -0.0, 378.0, 0.0, -36.0, -0.0]  ,
       [-945.0, -0.0, 4725.0, 0.0, -3150.0, -0.0, 630.0, 0.0, -45.0, -0.0]   ,
       [0.0, -10395.0, -0.0, 17325.0, 0.0, -6930.0, -0.0, 990.0, 0.0, -55.0, -0.0]
       ]
α, β = rm_hermite_prob(length(c))
cc = rec2coeff(α,β)

@testset "showpoly" begin
    for (i, c_) in enumerate(c)
        @test isapprox(norm(c_ - cc[i]), 0)
    end
end


text:
include("./jtalk.jl")
tts = Jtalk.JTalk(voicePath="mei_normal")

tts = Jtalk.JTalk()
println("available voices:")
for item in Jtalk.getVoices(tts)
    println(item[:path])
end
Jtalk.setVoice(tts,rand(Jtalk.getVoices(tts)))
println("current voice: $(Jtalk.getVoice(tts)[:name])")
println("dic= $(Jtalk.getDic(tts))")

println("s  = $(Jtalk.getS(tts))")
println("p  = $(Jtalk.getP(tts))")
println("a  = $(Jtalk.getA(tts))")
println("b  = $(Jtalk.getB(tts))")
println("r  = $(Jtalk.getR(tts))")
println("fm = $(Jtalk.getFm(tts))")
println("u  = $(Jtalk.getU(tts))")
println("jm = $(Jtalk.getJm(tts))")
println("jf = $(Jtalk.getJf(tts))")
println("g  = $(Jtalk.getG(tts))")

Jtalk.speakAsync(tts,"聞こえますか")
Jtalk.waitUntilDone(tts)


text:
module SubpixelRegistration

using AbstractFFTs
using Compat
using FFTW
using Statistics

export phase_offset, register, coregister, coregister!, fourier_shift

"""
    phase_offset(source, target; upsample_factor=1)

Return the shift between `source` and `target` along each axis by measuring the maximum in the cross-correlation between the images. This algorithm can achieve `1/upsample_factor` precision by locally upsampling the cross-correlation via a matrix-multiplication DFT.[^1]

# Examples

```jldoctest
julia> image = reshape(1.0:100.0, 10, 10);

julia> shift = (-1.6, 2.8)
(-1.6, 2.8)

julia> target = fourier_shift(image, shift);

julia> phase_offset(image, target)
(shift = (2.0, -3.0), error = 0.9999999999997532, phasediff = 0.0)

julia> phase_offset(image, target; upsample_factor=5)
(shift = (1.8, -2.8), error = 0.9999999971025185, phasediff = 0.0)

julia> @. isapprox(ans.shift, -1 * shift, atol=1/5)
(true, true)
```

[^1]: Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup, ["Efficient subpixel image registration algorithms,"](http://www.opticsinfobase.org/ol/fulltext.cfm?uri=ol-33-2-156&id=148843) Opt. Lett. 33, 156-158 (2008)
"""
function phase_offset(source::AbstractArray, target::AbstractArray; kwargs...)
    plan = plan_fft(source)
    return phase_offset(plan, plan * source, plan * target; kwargs...)
end

"""
    phase_offset(plan, source_freq, target_freq; upsample_factor=1)

Returns the phase shift between the two images which have already been Fourier transformed with the given `plan`.
"""
function phase_offset(plan, source_freq::AbstractMatrix{<:Complex{T}}, target_freq; upsample_factor=1) where T
    # whole-pixel shift
    # compute cross-correlation via iFFT
    image_product = transpose(source_freq) .* target_freq'
    # phase normalization
    @. image_product /= max(abs(image_product), 100 * eps(T))
    # ifft to calculate cross correlation
    if isone(upsample_factor)
        # no upsampling means we can modify this array
        cross_correlation = ifft!(image_product)
    else
        cross_correlation = plan \ image_product
    end

    # locate maximums
    maxima, maxidx = @compat findmax(abs, cross_correlation)
    shape = size(source_freq)
    midpoints = map(ax -> (first(ax) + last(ax)) / 2, axes(source_freq))
    idxoffset = map(first, axes(cross_correlation))
    shift = @. float(ifelse(maxidx.I > midpoints, maxidx.I - shape, maxidx.I) - idxoffset)

    isone(upsample_factor) && return (;shift, calculate_stats(maxima, source_freq, target_freq)...)

    # upsample with matrix-multiply DFT
    shift = @. round(shift * upsample_factor) / upsample_factor
    upsample_region_size = ceil(upsample_factor * 1.5)
    # center of output array at dftshift + 1
    dftshift = div(upsample_region_size, 2)
    # matmul DFT
    sample_region_offset = @. dftshift - shift * upsample_factor
    cross_correlation = upsampled_dft(image_product, upsample_region_size, upsample_factor, sample_region_offset)
    maxima, maxidx = @compat findmax(abs, cross_correlation)
    shift = @. shift + (maxidx.I - dftshift - idxoffset) / upsample_factor

    stats = calculate_stats(maxima, source_freq, target_freq)
    return (;shift, stats...)
end

"""
    SubpixelRegistration.upsampled_dft(data, region_size, upsample_factor, offsets)

Calculate the cross-correlation in a region of size `region_size` via an upsampled DFT. The DFT uses matrix-multiplication to super-sample the input by `upsample_factor`. The frequencies will be shifted and centered around `offsets`.
"""
function upsampled_dft(data::AbstractMatrix{T}, region_size, upsample_factor, offsets) where {T<:Complex}
    shiftrange = 1:region_size
    idxoffset = map(first, axes(data))
    sample_rate = inv(upsample_factor)
    freqs = fftfreq(size(data, 2), sample_rate)
    kernel = @. cis(-2π * (shiftrange - offsets[2] - idxoffset[2]) * freqs')

    _data = kernel * data'

    freqs = fftfreq(size(data, 1), sample_rate)
    kernel = @. cis(2π * (shiftrange - offsets[1] - idxoffset[1]) * freqs')
    _data = kernel * _data'
    return _data
end

"""
    SubpixelRegistration.calculate_stats(crosscor_maxima, source_freq, target_freq)

Calculate the normalized root-mean-square error (NRMSE) and total phase difference between the two complex arrays, `source_freq` and `target_freq`, with maximum cross-correlation value `crosscor_maxima`
"""
function calculate_stats(crosscor_maxima, source_freq, target_freq)
    source_amp = mean(abs2, source_freq)
    target_amp = mean(abs2, target_freq)
    error = 1 - abs2(crosscor_maxima) / (source_amp * target_amp)
    phasediff = atan(imag(crosscor_maxima), real(crosscor_maxima))
    return (; error, phasediff)
end

"""
    fourier_shift(image, shift)

Shift the given `image` by `shift` along each axis, using the Fourier phase information.
"""
function fourier_shift(image, shift, args...)
    FT = plan_fft(image)
    shifted = fourier_shift!(FT * image, shift, args...)
    return real(FT \ shifted)
end

"""
    fourier_shift!(image_freq::AbstractMatrix{<:Complex}, shift, phasediff=0)

Shift the given image, which is already in frequency-space, by `shift` along each axis. Modifies `image_freq` inplace.
"""
function fourier_shift!(image_freq::AbstractMatrix{<:Complex}, shift, phasediff=0)
    shape = size(image_freq)

    freqs1 = fftfreq(shape[1])'
    freqs2 = fftfreq(shape[2])
    @. image_freq *= cis(-2π * (freqs1 * shift[1] + freqs2 * shift[2]) + phasediff)
    return image_freq
end

"""
    register(source, target; upsample_factor=1)

Register `target` image to `source` image by first finding the phase offset ([`phase_offset`](@ref)), and then Fourier shifting `target` with [`fourier_shift`](@ref).

# Examples

```jldoctest
julia> image = reshape(1.0:100.0, 10, 10);

julia> shift = (-1.6, 2.8)
(-1.6, 2.8)

julia> target = fourier_shift(image, shift);

julia> target_shift = register(image, target; upsample_factor=5);
```

# See also
[`phase_offset`](@ref)
"""
function register(source, target; kwargs...)
    plan = plan_fft(source)
    target_freq = plan * target
    result = phase_offset(plan, plan * source, target_freq; kwargs...)
    shifted = fourier_shift!(target_freq, result.shift, result.phasediff)
    return real(plan \ shifted)
end

"""
    coregister(cube; dims, refidx=firstindex(cube, dims), upsample_factor=1)

Coregister a cube of data along `dims`, using the `refidx` slice as the source frame. Other keyword arguments will be passed to [`phase_offset`](@ref)

# See also
[`coregister!`](@ref)
"""
coregister(cube::AbstractArray; kwargs...) = coregister!(copy(cube); kwargs...)

"""
    coregister!(cube; dims, kwargs...)

Coregister slices in `cube`, modifyingn it inplace.

# See also
[`coregister`](@ref)
"""
function coregister!(cube::AbstractArray; dims, refidx=firstindex(cube, dims), kwargs...)
    source = selectdim(cube, dims, refidx)
    plan = plan_fft(source)
    source_freq = plan * source
    @inbounds for idx in axes(cube, dims)[begin + 1:end]
        target = selectdim(cube, dims, idx)
        target_freq = plan * target
        # measure offset and fourier shift
        result = phase_offset(plan, source_freq, target_freq; kwargs...)
        fourier_shift!(target_freq, result.shift, result.phasediff)
        # target is a view, update inplace
        target .= real.(plan \ target_freq)
    end
    return cube
end

end


text:
Base.one(A::HA) where {HA<:HybridMatrix} = HA(one(parent(A)))

# This should make one(sized subarray) return SArray
@inline function StaticArrays._construct_sametype(a::Type{<:SSubArray{Tuple{S,S},T}}, elements) where {S,T}
    return SMatrix{S,S,T}(elements)
end
@inline function StaticArrays._construct_sametype(a::SSubArray{Tuple{S,S},T}, elements) where {S,T}
    return SMatrix{S,S,T}(elements)
end

Base.fill!(A::HybridArray, x) = fill!(parent(A), x)

@inline function Base.zero(a::HybridArray{S}) where {S}
    return HybridArray{S}(zero(parent(a)))
end

@inline function Base.zero(a::SSubArray{S,T}) where {S,T}
    return StaticArrays.zeros(SArray{S,T})
end


text:
abstract type Stencil end

###
struct Star <: Stencil end

const star = Star()


text:
module LinearAgents

import ..JuliaRL: step!, start!, get_action

import ..AbstractAgent

import ..FeatureCreators
import ..Learning

import ..FeatureCreators: AbstractFeatureCreator, create_features, feature_size
import ..Learning.LinearRL: AbstractQFunction, update!, Optimizer, LearningUpdate, feature_type
import ..AbstractQPolicy
import ..get

using Random

export LinearQAgent, TileCoderAgent, get_action, start!, step!

mutable struct LinearQAgent{QT<:AbstractQFunction, FC<:AbstractFeatureCreator, P<:AbstractQPolicy, T<:Number} <: AbstractAgent
    Q::QT
    fc::FC
    π::P
    γ::Float64
    lu::LearningUpdate
    ϕ_t::Array{T, 1}
    ϕ_tp1::Array{T, 1}
    action::Int64
    LinearQAgent(Q::QT, fc::FC, π::P, γ::Float64, lu::LearningUpdate) where{T<:Number, QT<:AbstractQFunction, FC<:AbstractFeatureCreator, P<:AbstractQPolicy} =
        new{QT, FC, P, feature_type(Q)}(Q, fc, π, γ, lu, zeros(feature_type(Q), feature_size(fc)), zeros(feature_type(Q), feature_size(fc)), 0)
end

get_action(agent::LinearQAgent, ϕ; rng=Random.GLOBAL_RNG, kwargs...) = get(agent.π, [agent.Q(ϕ, a) for a = agent.π.actions]; rng=rng)

function start!(agent::LinearQAgent, env_s_tp1; rng=Random.GLOBAL_RNG, kwargs...)

    agent.ϕ_t .= create_features(agent.fc, env_s_tp1)
    agent.action = get_action(agent, agent.ϕ_t; rng=rng)

    return agent.action
end

function step!(agent::LinearQAgent, env_s_tp1, r, terminal; rng=Random.GLOBAL_RNG, kwargs...)

    agent.ϕ_tp1 .= create_features(agent.fc, env_s_tp1)

    update!(
        agent.Q,
        agent.lu,
        agent.ϕ_t,
        agent.ϕ_tp1,
        r,
        agent.γ,
        1.0,
        terminal,
        agent.action)

    agent.ϕ_t .= agent.ϕ_tp1
    agent.action = get_action(agent, agent.ϕ_t; rng=rng)
    return agent.action
end

function TileCoderAgent(opt::LearningUpdate,
                        size_env_state::Integer,
                        num_actions::Integer,
                        tilings::Integer,
                        tiles::Integer,
                        γ::Float64,
                        policy::AbstractQPolicy)

    fc = FeatureCreators.TileCoder(tilings, tiles, size_env_state; wrap=false, wrapwidths=0.0)
    num_features_per_action = (tilings*(tiles+1)^size_env_state)
    num_features = num_features_per_action*num_actions
    Q = Learning.LinearRL.ActionSparseQFunction(
        num_features_per_action,
        num_actions)
    return LinearQAgent(Q, fc, policy, γ, opt)
end

end

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
xnorm(x::AbstractVector) = norm(x, 2)
xnorm(x::CuVector) = CUBLAS.nrm2(x)

@enum(
    SolveStatus,
    Unsolved,
    MaxIterations,
    NotANumber,
    Converged,
    Diverged,
)

"""
    bicgstab(A, b, P, xi;
             tol=1e-8,
             maxiter=size(A, 1),
             verbose=false,
             maxtol=1e20)

BiCGSTAB implementation according to

> Van der Vorst, Henk A.
> "Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems."
> SIAM Journal on scientific and Statistical Computing 13, no. 2 (1992): 631-644.

"""
function bicgstab(A, b, P, xi;
                  tol=1e-8, maxiter=size(A, 1), verbose=false, maxtol=1e20)
    # parameters
    n    = size(b, 1)
    mul!(xi, P, b)
    ri   = b - A * xi
    br0  = copy(ri)
    rho0 = 1.0
    alpha = 1.0
    omega0 = 1.0
    vi = similar(xi)
    pi = similar(xi)
    fill!(vi, 0)
    fill!(pi, 0)

    rhoi   = copy(rho0)
    omegai = copy(omega0)
    residual = copy(b)

    y = similar(pi)
    s = similar(pi)
    z = similar(pi)
    t1 = similar(pi)
    t2 = similar(pi)
    pi1 = similar(pi)
    vi1 = similar(pi)
    t = similar(pi)

    go = true
    status = Unsolved
    iter = 1
    restarts = 0
    while go
        rhoi1 = dot(br0, ri) ;
        if abs(rhoi1) < 1e-20
            restarts += 1
            ri .= b
            mul!(ri, A, xi, -1.0, 1.0)
            br0 .= ri
            residual .= b
            rho0 = 1.0
            rhoi = rho0
            rhoi1 = dot(br0,ri)
            alpha = 1.0
            omega0 = 1.0
            omegai = 1.0
            fill!(vi, 0.0)
            fill!(pi, 0.0)
        end
        beta = (rhoi1/rhoi) * (alpha / omegai)
        pi1 .= ri .+ beta .* (pi .- omegai .* vi)
        mul!(y, P, pi1)
        mul!(vi1, A, y)
        alpha = rhoi1 / dot(br0, vi1)
        s .= ri .- (alpha .* vi1)

        mul!(z, P, s)
        mul!(t, A, z)
        mul!(t1, P, t)
        mul!(t2, P, s)
        omegai1 = dot(t1, t2) / dot(t1, t1)
        xi .= xi .+ alpha .* y .+ omegai1 .* z

        mul!(residual, A, xi, 1.0, -1.0)
        anorm = xnorm(residual)

        if verbose
            @printf("%4d %10.4e\n", iter, anorm)
        end

        if isnan(anorm)
            go = false
            status = NotANumber
        end
        if anorm < tol
            go = false
            status = Converged
            verbose && println("Tolerance reached at iteration $iter")
        elseif anorm > maxtol
            go = false
            status = Diverged
        elseif maxiter == iter
            go = false
            status = MaxIterations
            verbose && println("Not converged")
        end

        ri     .= s .- omegai1 .* t
        rhoi   = rhoi1
        pi     .= pi1
        vi     .= vi1
        omegai = omegai1
        iter   += 1
    end

    return xi, iter, status
end