text:
@enum(CbSolStrat, CPXCALLBACKSOLUTION_CHECKFEAS, CPXCALLBACKSOLUTION_PROPAGATE)

@enum(CbInfo,
    CPXCALLBACKINFO_THREADID, CPXCALLBACKINFO_NODECOUNT,
    CPXCALLBACKINFO_ITCOUNT, CPXCALLBACKINFO_BEST_SOL, CPXCALLBACKINFO_BEST_BND,
    CPXCALLBACKINFO_THREADS, CPXCALLBACKINFO_FEASIBLE, CPXCALLBACKINFO_TIME,
    CPXCALLBACKINFO_DETTIME)

mutable struct CallbackContext
    model::Model
    context::Ptr{Cvoid}
end

function callback_wrapper(context::Ptr{Cvoid},
                          context_id::Clong,
                          user_handle::Ptr{Cvoid})
    (model, callback_func) = unsafe_pointer_to_objref(user_handle)::Tuple{Model, Function}
    # Run the `callback_func`.
    callback_func(CallbackContext(model, context), context_id)
    # Exit normally.
    return Cint(0)
end

function return_status_or_throw(model, status)
    if status != 0
        throw(CplexError(model.env, status))
    else
        return status
    end
end

 """
    cpx_callbacksetfunc(model::Model, context_mask::Clong, callback_func::Function)

Set the general callback function of model to `callback_func`.

`callback_func` is a function that takes two arguments. The first is a
`CallbackContext` instance, and the second is a `Clong` indicating where the
callback is being called from.

### Examples

     function my_callback(cb_context::CPLEX.CallbackContext, context_id::Clong)
        CPLEX.cpx_callbackabort(cb_context)
    end
    CPLEX.cpx_callbacksetfunc(model, Clong(0), my_callback)
"""
function cbsetfunc(model::Model, context_mask::Clong, callback_func::Function)
    c_callback = @cfunction(callback_wrapper, Cint, (Ptr{Cvoid}, Clong, Ptr{Cvoid}))
    user_handle = (model, callback_func)::Tuple{Model, Function}
    return_status = @cpx_ccall(callbacksetfunc,
                    Cint,
                    (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Ptr{Cvoid}, Any),
                    model.env, model.lp, context_mask, c_callback, user_handle)
    # We need to keep a reference to the callback function so that it isn't
    # garbage collected.
    model.callback = user_handle
    return return_status_or_throw(model, return_status)
end

function cbabort(callback_data::CallbackContext)
    return_status = @cpx_ccall(
        callbackabort, Cint, (Ptr{Cvoid},), callback_data.context)
    return return_status_or_throw(callback_data.model, return_status)
end

function cbgetrelaxationpoint(
        callback_data::CallbackContext, x::Vector{Float64}, begin_index::Int,
        end_index::Int, obj::Ref{Float64})
    return_status = @cpx_ccall(callbackgetrelaxationpoint,
        Cint,
        (Ptr{Cvoid}, Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}),
        callback_data.context, x, begin_index - 1, end_index - 1, obj)
    return return_status_or_throw(callback_data.model, return_status)
end

function cbpostheursoln(
        callback_data::CallbackContext, num_entries::Int, columns::Vector{Int},
        values::Vector{Cdouble}, objective_value::Cdouble, strat::CbSolStrat)
    return_status = @cpx_ccall(callbackpostheursoln,
        Cint,
        (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}, Cdouble, Cint),
        callback_data.context, num_entries, columns .- 1, values,
            objective_value, strat)
    return return_status_or_throw(callback_data.model, return_status)
end

function cbaddusercuts(
        callback_data::CallbackContext, reject_count::Int, nonzero_count::Int,
        rhs::Cdouble, sense::Char, sparse_index_begin::Vector{Int},
        sparse_index::Vector{Int}, sparse_index_coeff::Vector{Cdouble},
        purgeable::Int, lcl::Int)
    return_status = @cpx_ccall(callbackaddusercuts,
        Cint,
        (Ptr{Cvoid}, Cint, Cint, Ref{Cdouble}, Ptr{UInt8}, Ptr{Cint}, Ptr{Cint},
            Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}),
        callback_data.context, reject_count, nonzero_count, rhs, string(sense),
            sparse_index_begin .- 1, sparse_index .- 1, sparse_index_coeff,
            purgeable, lcl)
    return return_status_or_throw(callback_data.model, return_status)
end

# Check if the callback for the context CPX_CALLBACKCONTEXT_CANDIDATE is invoked by a feasible integer point
function cbcandidateispoint(callback_data::CallbackContext, is_point::Int)
    return_status = @cpx_ccall(callbackcandidateispoint,
        Cint, (Ptr{Cvoid}, Ptr{Cint}), callback_data.context, is_point)
    return return_status_or_throw(callback_data.model, return_status)
end

# Check if the callback for the context CPX_CALLBACKCONTEXT_CANDIDATE is invoked by an unbounded relaxation
function cbcandidateisray(callback_data::CallbackContext, is_ray::Int)
    return_status = @cpx_ccall(callbackcandidateisray,
        Cint, (Ptr{Cvoid}, Ptr{Cint}), callback_data.context, is_ray)
    return return_status_or_throw(callback_data.model, return_status)
end

function cbgetcandidatepoint(
        callback_data::CallbackContext, x::Vector{Float64}, begin_index::Int,
        end_index::Int, obj::Ref{Float64})
    return_status = @cpx_ccall(callbackgetcandidatepoint,
        Cint,
        (Ptr{Cvoid}, Ptr{Cdouble}, Cint, Cint, Ref{Cdouble}),
        callback_data.context, x, begin_index - 1, end_index - 1, obj)
    return return_status_or_throw(callback_data.model, return_status)
end

function cbgetcandidateray(
        callback_data::CallbackContext, x::Vector{Cdouble}, begin_index::Int,
        end_index::Int)
    return_status = @cpx_ccall(callbackgetcandidateray,
        Cint,
        (Ptr{Cvoid}, Ptr{Cdouble}, Cint, Cint),
        callback_data.context, x, begin_index - 1, end_index - 1)
    return return_status_or_throw(callback_data.model, return_status)
end

# This function may have bug, should be tested.
# function cbgetfunc(model::Model, context_id::Clong, callback_func::Function)
#     c_callback = @cfunction(callback_wrapper, Cint, (Ptr{Cvoid}, Clong, Ptr{Cvoid}))
#     user_handle = (model, callback_func)::Tuple{Model, Function}
#     return_status = @cpx_ccall(callbackgetfunc,
#                     Cint,
#                     (Ptr{Cvoid}, Ptr{Cvoid}, Ref{Clong}, Ptr{Cvoid}, Any),
#                     model.env, model.lp, context_mask, c_callback, user_handle)
#     if return_status != 0
#         throw(CplexError(model.env, return_status))
#     end
#     model.callback = user_handle
#     return return_status
# end

function cbgetincumbent(
        callback_data::CallbackContext, x::Vector{Cdouble}, begin_index::Int,
        end_index::Int, obj::Ref{Float64})
    return_status = @cpx_ccall(callbackgetincumbent,
        Cint,
        (Ptr{Cvoid}, Ptr{Cdouble}, Cint, Cint, Ref{Cdouble}),
        callback_data.context, x, begin_index - 1, end_index - 1, obj)
    return return_status_or_throw(callback_data.model, return_status)
end

function cbgetinfodbl(
        callback_data::CallbackContext, cbinfo::CbInfo, data_p::Cdouble)
    return_status = @cpx_ccall(callbackgetinfodbl,
        Cint,
        (Ptr{Cvoid}, Cint, Ref{Cdouble}),
        callback_data.context, cbinfo, data_p)
    return return_status_or_throw(callback_data.model, return_status)
end

function cbgetinfoint(
        callback_data::CallbackContext, cbinfo::CbInfo, data_p::Int)
    return_status = @cpx_ccall(callbackgetinfoint,
        Cint,
        (Ptr{Cvoid}, Cint, Ref{Cint}),
        callback_data.context, cbinfo, data_p)
    return return_status_or_throw(callback_data.model, return_status)
end

function cbgetinfolong(
        callback_data::CallbackContext, cbinfo::CbInfo, data_p::Clong)
    return_status = @cpx_ccall(callbackgetinfolong,
        Cint,
        (Ptr{Cvoid}, Cint, Ref{Clong}),
        callback_data.context, cbinfo, data_p)
    return return_status_or_throw(callback_data.model, return_status)
end

function cbrejectcandidate(
        callback_data::CallbackContext, reject_count::Int, nonzero_count::Int,
        rhs::Cdouble, sense::Char, sparse_index_begin::Vector{Int},
        sparse_index::Vector{Int}, sparse_index_coeff::Vector{Cdouble})
    return_status = @cpx_ccall(callbackrejectcandidate,
        Cint,
        (Ptr{Cvoid}, Cint, Cint, Ref{Cdouble}, Ptr{UInt8}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}),
        callback_data.context, reject_count, nonzero_count, rhs, string(sense),
            sparse_index_begin .- 1, sparse_index .- 1, sparse_index_coeff)
    return return_status_or_throw(callback_data.model, return_status)
end

function cbgetlocallb(
        callback_data::CallbackContext, lb::Vector{Cdouble}, begin_index::Int,
        end_index::Int)
    return_status = @cpx_ccall(callbackgetlocallb,
        Cint,
        (Ptr{Cvoid}, Ptr{Cdouble}, Cint, Cint),
        callback_data.context, lb, begin_index - 1, end_index - 1)
    return return_status_or_throw(callback_data.model, return_status)
end

function cbgetlocalub(
        callback_data::CallbackContext, ub::Vector{Cdouble}, begin_index::Int,
        end_index::Int)
    return_status = @cpx_ccall(callbackgetlocalub,
        Cint,
        (Ptr{Cvoid}, Ptr{Cdouble}, Cint, Cint),
        callback_data.context, ub, begin_index - 1, end_index - 1)
    return return_status_or_throw(callback_data.model, return_status)
end

function cbgetgloballb(
        callback_data::CallbackContext, lb::Vector{Cdouble}, begin_index::Int,
        end_index::Int)
    return_status = @cpx_ccall(callbackgetgloballb,
        Cint,
        (Ptr{Cvoid}, Ptr{Cdouble}, Cint, Cint),
        callback_data.context, lb, begin_index - 1, end_index - 1)
    return return_status_or_throw(callback_data.model, return_status)
end

function cbgetglobalub(
        callback_data::CallbackContext, ub::Vector{Cdouble}, begin_index::Int,
        end_index::Int)
    return_status = @cpx_ccall(callbackgetglobalub,
        Cint,
        (Ptr{Cvoid}, Ptr{Cdouble}, Cint, Cint),
        callback_data.context, ub, begin_index - 1, end_index - 1)
    return return_status_or_throw(callback_data.model, return_status)
end


text:
@testset "Polynomial equality constraint with domain with $(factory.constructor)" for factory in sdp_factories
    @polyvar(x, y)

    m = SOSModel(factory)

    @constraint(m, x^3 + x*y^2 == x, domain=(@set x^2 + y^2 >= 1 && x^2 + y^2 <= 1))

    JuMP.optimize!(m)

    @test JuMP.primal_status(m) == MOI.FeasiblePoint
end

ext: jl


text:
using LinearAlgebra

mutable struct NesterovAcceleration{R <: Real, C <: Union{R, Complex{R}}, T <: AbstractArray{C}}
    k::Int
    s::T
end

function NesterovAcceleration(x::T) where {R <: Real, C <: Union{R, Complex{R}}, T <: AbstractArray{C}}
    NesterovAcceleration{R, C, T}(0, zero(x))
end

function update!(L::NesterovAcceleration{R, C, T}, s, y) where {R, C, T}
    L.s .= s
    L.k += 1
    return L
end

function reset!(L::NesterovAcceleration{R, C, T}) where {R, C, T}
    L.k = 0
end

import Base: *

function (*)(L::NesterovAcceleration, v)
    w = similar(v)
    mul!(w, L, v)
end

import LinearAlgebra: mul!

function mul!(d::T, L::NesterovAcceleration{R, C, T}, v::T) where {R, C, T}
    if L.k == 0
        d .= 0
    elseif L.k == 1
        d .= v
    else
        d .= v .- (L.k - 1)/(L.k + 2) .* L.s
    end
end


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_ps_m7_opt_rho6_0(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_ps_m7_opt_rho6_0!(A_copy)
end

@inline function exp_ps_m7_opt_rho6_0!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=10
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    memslots9 = similar(A,T)
    memslots10 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 Ba7 Bb7 B7 Ba8 Bb8 B8 T2k10
    # Computing Bb2 = x*I+x*A
    coeff1=0.40493955401846876
    coeff2=0.4037226736739943
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=0.40493955401846876
    coeff2=0.4037226736739943
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=0.05431906428723421
    coeff2=0.9217267072847067
    coeff3=0.253861883536394
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=0.7792960542315325
    coeff2=0.4594622982985042
    coeff3=0.085965940313993
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.07101339405931287
    coeff2=0.2183449794974525
    coeff3=0.6325401240010657
    coeff4=0.5077334065861822
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=0.7611160284006683
    coeff2=0.14582635451246748
    coeff3=0.01346186908673783
    coeff4=0.005570337889186752
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.01237456812732672
    coeff2=-0.007417121502402069
    coeff3=0.36461901499271315
    coeff4=0.7318998265697871
    coeff5=0.42311906620920164
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.5304750742911883
    coeff2=0.5804911103105471
    coeff3=0.15347975441496023
    coeff4=0.02847544670127242
    coeff5=0.017703603503626326
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-1.1142498138003632e-9
    coeff2=-0.00044073573263142074
    coeff3=-0.009944137968525259
    coeff4=0.05904148896574505
    coeff5=0.16671013087805137
    coeff6=0.9891460760326148
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=5.222157680049785e-8
    coeff2=-2.3658628414902608e-9
    coeff3=-2.6488325501827903e-10
    coeff4=6.060309152649078e-11
    coeff5=1.0661525964254732e-12
    coeff6=8.967326795327424e-15
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing Ba7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=0.00030937511267045
    coeff2=-0.004992129336289725
    coeff3=0.10003129254607988
    coeff4=-0.21265617105393061
    coeff5=-0.2535622236995746
    coeff6=0.9892639607059927
    coeff7=1.0457211702940734e-7
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=0.006312989983639836
    coeff2=3.05349059290055e-5
    coeff3=1.2830919585643768e-5
    coeff4=3.7913704630241052e-6
    coeff5=2.3225691022612054e-6
    coeff6=9.290484746562112e-8
    coeff7=1.0000000001559715
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B7 with operation: mult
    mul!(memslots9,memslots2,memslots3)
    # Deallocating Ba7 in slot 2
    # Deallocating Bb7 in slot 3
    # Computing Ba8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=0.990632734847538
    coeff2=-0.037157343305066634
    coeff3=0.09692046328345978
    coeff4=0.4186089963980845
    coeff5=0.6631289790906699
    coeff6=0.7614234275617816
    coeff7=-1.005632757710081e-8
    coeff8=-0.0713903017943822
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=0.1078899138880161
    coeff2=0.03205343917230924
    coeff3=0.005656891024698354
    coeff4=0.010962105218802998
    coeff5=0.0013745451070322613
    coeff6=-0.004668205384343422
    coeff7=-1.152223044748429e-5
    coeff8=1.0009187967337483
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B8 with operation: mult
    mul!(memslots10,memslots2,memslots3)
    # Deallocating Ba8 in slot 2
    # Deallocating Bb8 in slot 3
    # Computing T2k10 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7+x*B8
    coeff1=0.7754984117964325
    coeff2=0.3475417761570095
    coeff3=0.4359396287990991
    coeff4=0.007610711188155261
    coeff5=0.43170388767357853
    coeff6=0.037684845865425066
    coeff7=1.5754616856674562e-7
    coeff8=0.0007760687742691415
    coeff9=1.0058204367206869
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9 .+ coeff9.*memslots10
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    # Deallocating B7 in slot 9
    # Deallocating B8 in slot 10
    return memslots1 # Returning T2k10
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
"""
    module ProcessBNF

Special processing routines for back-and-forth pumping tissue slice culture experiments.

`function process_spectra`
`function generate_protocol`
"""
module ProcessBNF

using DataFrames
import CSV
using Dates
import NMR
using LsqFit
import Dates
import Plots
using Distributions
using Statistics

export process_spectra, generate_protocol,rate_plot, linearRegP

"""
    function process_fid(n,dataStore,dataLocation)

process the FID contained in the file `dataStore`*`dataLocation`/`n`.
Returns a `Data1D` object with the processed spectrum.
"""
function process_fid(n,dataStore,dataLocation)
        SW = 19.8269803808675;
        SW_h= 600*SW;
        d=NMR.readBrukerFID(dataStore * dataLocation * "/$(n)/fid")[77:end];
        fid=NMR.Data1D(d,0.,length(d)/SW_h);
        spect=NMR.FourierTransform(fid,SI=16384,PPM=600.,CTR=4.82,LB=1.0π)
        spect=NMR.PhaseCorrect(spect,Ph0=2.95,Ph1=-0)
        spect=NMR.cut(spect,0.,10.)
        spect=NMR.BaseLineCorrect(spect,kfactor=4.5,regions=256,order=31)
        #spect-=NMR.medianBaseline(spect)
        return spect
end

"""
    function generate_protocol(fname)

Generates a protocol csv file, which each row corresponding to a single
NMR spectrum, in the directory designated by `fname`.
"""
function generate_protocol(fname)
    # cd(fname)

    d=readdir(fname)

    boxf=findfirst(x->match(r"^box-[0-9]*\.txt",x)!=nothing,d)
    protof=findfirst(x->match(r"^labsmith-[0-9]*\.txt",x)!=nothing,d)

    # @show d[[boxf,protof]]

    # clean up labsmith protocol file

    infile=open(fname*d[protof],"r")
    outfile=open(fname*"protocol-raw.csv","w")
    write(outfile,"Date,Time,PumpName,FlowRate,CurrVol,TargetVol,Push,Pull,BackForthActive,OtherTarget\n")
    while !eof(infile)
        line=readline(infile)
        write(outfile,replace(line,r"\0"=>""))
        write(outfile,"\n")
    end
    close(infile)
    close(outfile)

    # read in temperature log

    tlograw=CSV.read(fname*d[boxf], DataFrame,types=[String,Float64,Float64,Float64,Float64,Float64,Float64,String])

    tlograw.DateTime=[try DateTime(x,"dd/mm/yyyy HH:MM:SS") catch end for x in tlograw.Time]
    tlog=tlograw[tlograw.DateTime.!=nothing,:]

    # read in protocol

    protocol=CSV.read(fname*"protocol-raw.csv",DataFrame)
    protocol.DateTime=Date.(protocol[!,"Date"],"dd/mm/yyyy").+protocol[!,"Time"]

    # look up temperature data from tlog

    protocol.Temperature = [  begin
            dt=f.DateTime ;
            tlog[findmin(abs.(tlog.DateTime.-dt))[2],:Temperature]
        end
    for f in eachrow(protocol) ]

    # select only the rows in which back-and-forth pumping is active

    protocol_select=protocol[(protocol.BackForthActive.==true) .& ( (protocol.Push.==1) .| (protocol.Pull.==1) ),:]

    protocol_final=DataFrame()
    lasttime=protocol_select[1,:DateTime]
    for k in eachrow(protocol_select)
        if((k.DateTime-lasttime)>Millisecond(2000)) push!(protocol_final,k) end
        lasttime=k.DateTime
    end

    # time differentials. Long duration (more than 4min) is indicative of medium refresh

    deltaT=[protocol_final[k,:DateTime]-protocol_final[k-1,:DateTime] for k in 2:nrow(protocol_final)]
    prepend!(deltaT,[Minute(9)])
    protocol_final.duration = deltaT

    CSV.write(fname*"protocol-final.csv",protocol_final)

end

"""
    function rate_plot(protocol::DataFrame,compound::String; stride=50, model=(t,p)->p[1].+p[2].*t)

compute linear fits to segments of the data of length `stride`, and return a tuple `(fits,c)`, where
`c` is a plot of the concentration of `compound` vs time, and `fits` is a vector of `LsqFit` objects
containing information on the fit in each of the segments.
"""
function rate_plot(protocol::DataFrame,compound::String; stride=50, model=(t,p)->p[1].+p[2].*t)
    ydata=protocol[!,compound]
    xdata=Dates.value.(protocol[!,"DateTime"].-protocol[1,"DateTime"]) ./ 3.6e6

    p0=[1.0,-1.0]
    c=Plots.plot(layout=[1,1],size=(500,600),dpi=200,legend=false)
    a=Plots.plot!(c,xdata,ydata,xlabel="Time [h]",ylabel="Concentration [mM]",xlims=[0,xdata[end]],label=compound,subplot=1,link=:x)
    fits=Array{Any,1}()

    for k=1:stride:(length(ydata))
        rng=k:(k+stride-1)
        f1=curve_fit(model,xdata[rng].-xdata[k],ydata[rng],p0)
        Plots.plot!(c,xdata[rng],model(xdata[rng].-xdata[k],f1.param),linewidth=2,linecolor=:red,label=false,link=:x,subplot=1)
        push!(fits,f1)
    end

    b=Plots.plot!(c,xdata[(stride>>1):stride:end],[f.param[2] for f in fits],
        seriestype=:scatter,err=[margin_error(f,0.1)[2] for f in fits],
        xlims=[0,xdata[end]],
        xlabel="Time [h]",
        ylabel="Rate [mM/h]",
        label=compound,link=:x,subplot=2)

    return (fits,c)

end

"""
        function linearRegP(f)

compute p-value for linear regression fit `f`, which must be a result of a linear
        regression through `LsqFit.jl`, with parameter set intercept,slope.
"""
function linearRegP(f)
        F=((var(f.jacobian*f.param.+f.resid)-var(f.resid))/var(f.resid))*(length(f.resid)-2)
        p=1-cdf(FDist(1,length(f.resid)-2),F)
end

end

max_stars_repo_name: marcel-utz/NMR

max_issues_repo_name: marcel-utz/NMR

max_forks_repo_name: marcel-utz/NMR


text:
# Memcopy based on AvX instructions.
"""
    _mov(src::Ptr, dest::Ptr, ::Val{N}) where {N}

Unroll move `N` elements from `src` to `ptr`.
"""
@generated function _mov(
        ::Type{SIMD.Vec{N,T}},
        dest::Ptr{UInt8},
        src::Ptr{UInt8},
        ::Val{U},
    ) where {N,T,U}

    loads = load_impl(SIMD.Vec{N,T}, U)
    stores = store_impl(SIMD.Vec{N,T}, U)
    return quote
        $(Expr(:meta, :inline))
        _src = convert(Ptr{$T}, src)
        _dest = convert(Ptr{$T}, dest)
        $(loads...)
        $(stores...)
    end
end

function load_impl(::Type{T}, U::Integer) where {T <: SIMD.Vec}
    return map(0:U-1) do j
        x = Symbol("i_$j")
        # Vector Load
        # Val(true) implies aligned load
        :($x = SIMD.vload($T, _src + $(sizeof(T)) * $j, nothing, Val(true)))
    end
end

function store_impl(::Type{T}, U::Integer) where {T <: SIMD.Vec}
    return map(0:U-1) do j
        x = Symbol("i_$j")
        # Vector Store
        # first Val(true) implies aligned store
        # second Val(true) implies nontemporal store
        #
        # See SIMD.jl for documentation.
        :(SIMD.vstore($x, _dest + $(sizeof(T)) * $j, nothing, Val(true), Val(true)))
    end
end

function sfence()
    Base.llvmcall(
        raw"""
        tail call void asm sideeffect "sfence", "~{memory},~{dirflag},~{fpsr},~{flags}"()
        ret void
        """,
        Nothing,
        Tuple{},
    )
    return nothing
end

# Chunk up an array into N views so the start of each array is aligned
# on a cache boundary.
function aligned_chunk(nbytes, nchunks)
    isone(nchunks) && return (0, nbytes)

    # Round down `nbytes` to a multiple of 64
    extra = mod(nbytes, 64)
    _nbytes = nbytes - extra

    bytes_per_chunk, tail = divrem(_nbytes, nchunks)
    chunk_remainder = mod(bytes_per_chunk, 64)
    bytes_per_chunk = bytes_per_chunk - chunk_remainder

    last_chunk = bytes_per_chunk + extra + nchunks * chunk_remainder + tail
    return (bytes_per_chunk, last_chunk)
end

_isaligned(x::Integer) = iszero(mod(x, 64))
_isaligned(x::Ptr) = _isaligned(convert(Int, x))

# Top level entry point
#
#
# NOTE: Dest and Src must not alias!

"""
    memcpy!(dest::AbstractArray, src::AbstractArray, [toremote = false]; [forcesingle])

Copy the contents from `src` to `dest` using non-temporal AVX store instructions.
Pass `toremote = true` if `dest` lives in Remote for better performance.
Set keyword `forcesingle = true` to force only one thread to do the copy operations.

LIMITATIONS
-----------
* `dest` and `src` must not alias at all.
* The base pointers for `dest` and `src` must be 64-byte aligned.
"""
function memcpy!(dest::AbstractArray{T}, src::AbstractArray{T}; kw...) where {T}
    return _memcpy!(dest, src; kw...)
end

# Specialize `memcpy` for CachedArrays.
# TODO: Figure out maximum number of threads for copying from Remote to Remote
function memcpy!(dest::CachedArray{T}, src::AbstractArray{T}; nthreads = nothing) where {T}
    # Perform standard thread chosing logic
    if nthreads === nothing
        maxthreads = Threads.nthreads()
        nthreads = pool(dest) == Remote ? min(maxthreads, 4) : maxthreads
    end

    return _memcpy!(dest, src; nthreads = nthreads)
end

function _memcpy!(dest::AbstractArray{T}, src::AbstractArray{T}; nthreads = nothing) where {T}
    # Do type check
    if !isbitstype(T)
        throw(ArgumentError("Can only memcpy isbits types."))
    end

    # Do length checks
    if length(dest) != length(src)
        throw(ArgumentError("Source and Destination must have the same lengths."))
    end

    dest_ptr = pointer(dest)
    src_ptr = pointer(src)

    # Do alignment checks
    if !_isaligned(dest_ptr) || !_isaligned(src_ptr)
        @warn "Copying unaligned"
        copyto!(dest, src)
        return nothing
    end

    return _memcpy!(
        dest_ptr,
        src_ptr,
        length(dest);
        nthreads = nthreads === nothing ? Threads.nthreads() : nthreads
    )
end

function _memcpy!(dest::Ptr{T}, src::Ptr{T}, len; kw...) where {T}
    iszero(sizeof(T)) && error("Cannot move items of type $T")
    return _memcpy!(
        convert(Ptr{UInt8}, dest),
        convert(Ptr{UInt8}, src),
        sizeof(T) * len;
        kw...
    )
end

function _memcpy!(dest::Ptr{Nothing}, src::Ptr{Nothing}, len; kw...)
    return _memcpy!(convert(Ptr{UInt8}, dest), convert(Ptr{UInt8}, src), len; kw...)
end

function _memcpy!(dest::Ptr{UInt8}, src::Ptr{UInt8}, bytes; nthreads = nothing)
    # Handle `isnothing` case first.
    if nthreads === nothing
        unsafe_memcpy!(+, dest, src, bytes)
        sfence()
        return nothing
    end

    # Now, deal with the specified nmber of threads
    bytes_per_chunk, last_chunk = aligned_chunk(bytes, nthreads)
    @timeit "Memcpy" Polyester.@batch per=core for i in Base.OneTo(nthreads)
        offset = bytes_per_chunk * (i-1)
        copybytes = (i == nthreads) ? last_chunk : bytes_per_chunk
        unsafe_memcpy!(+, dest + offset, src + offset, copybytes)
        sfence()
    end
    return nothing
end

const AVX512BYTES = 64
_adjust(::typeof(+), x::Ptr{UInt8}, len) = x
inv(::typeof(+)) = -

"""
    unsafe_memcpy!(f, dest::Ptr{UInt8}, src::Ptr{UInt8}, len)

Copy `len` bytes from `src` to `dest`.
If `f == Base.+`, the copy will be done from low addresses to high addresses.
If `f == Base.-`, the copy will be done from low addresses to high addresses.
"""
function unsafe_memcpy!(f, dest::Ptr{UInt8}, src::Ptr{UInt8}, len)
    # By this point, assume that `dest`  and `src` are aligned.
    dest = _adjust(f, dest, len)
    src = _adjust(f, src, len)

    # While we can, manually unroll the AVX512 move operation by 32 times.
    while (len >= 32 * AVX512BYTES)
        _mov(SIMD.Vec{8,Int64}, dest, src, Val{32}())
        len = inv(f)(len, 32 * AVX512BYTES)
        dest = f(dest, 32 * AVX512BYTES)
        src = f(src, 32 * AVX512BYTES)
    end

    if (len >= 16 * AVX512BYTES)
        _mov(SIMD.Vec{8,Int64}, dest, src, Val{16}())
        len = inv(f)(len, 16 * AVX512BYTES)
        dest = f(dest, 16 * AVX512BYTES)
        src = f(src, 16 * AVX512BYTES)
    end

    if (len >= 8 * AVX512BYTES)
        _mov(SIMD.Vec{8,Int64}, dest, src, Val{8}())
        len = inv(f)(len, 8 * AVX512BYTES)
        dest = f(dest, 8 * AVX512BYTES)
        src = f(src, 8 * AVX512BYTES)
    end

    if (len >= 4 * AVX512BYTES)
        _mov(SIMD.Vec{8,Int64}, dest, src, Val{4}())
        len = inv(f)(len, 4 * AVX512BYTES)
        dest = f(dest, 4 * AVX512BYTES)
        src = f(src, 4 * AVX512BYTES)
    end

    if (len >= 2 * AVX512BYTES)
        _mov(SIMD.Vec{8,Int64}, dest, src, Val{2}())
        len = inv(f)(len, 2 * AVX512BYTES)
        dest = f(dest, 2 * AVX512BYTES)
        src = f(src, 2 * AVX512BYTES)
    end

    if (len >= 1 * AVX512BYTES)
        _mov(SIMD.Vec{8,Int64}, dest, src, Val{1}())
        len = inv(f)(len, 1 * AVX512BYTES)
        dest = f(dest, 1 * AVX512BYTES)
        src = f(src, 1 * AVX512BYTES)
    end

    # Normal MEMCPY the rest.
    iszero(len) || unsafe_copyto!(dest, src, len)
    return nothing
end


text:
"""
wordcount(sentence)

Program to find word count in the given sentence. The program will return count with the word.
A full description of the algorithm can be found on [exercism](https://exercism.org/tracks/julia/exercises/word-count)

# Arguments:
- `sentence`: The sentence to find the word count.

# Examples/Tests
```julia
julia> wordcount(The quick brown fox jumps over the lazy dog)
Dict{Any, Any}("jumps" => 1, "the" => 2, "brown" => 1, "over" => 1, "quick" => 1, "lazy" => 1, "dog" => 1, "fox" => 1)

julia> wordcount(the sky is blue and beautiful)
Dict{Any, Any}("and" => 1, "the" => 1, "sky" => 1, "blue" => 1, "is" => 1, "beautiful" => 1)

```

# Algorithm:

```julia

for word in eachmatch(reg_expression, sentence)
    if !haskey(counts, word.match)
        counts[word.match] = 1
    else
        counts[word.match] += 1
    end
end

```

# References:
(https://exercism.org/tracks/julia/exercises/word-count)

```

# Contributed by:- [Ihjass Thasbekha](https://github.com/Ihjass)
"""
function word_count(sentence)
    lowercase_sentence = lowercase(sentence)
    counts = Dict()
    regex = r"(\w+('\w)?)"
    for w in eachmatch(regex, lowercase_sentence)
        if !haskey(counts, w.match)
            counts[w.match] = 1
        else
            counts[w.match] += 1
        end
    end
    return counts
end

max_stars_repo_name: standardgalactic/Julia-algorithms


text:
# Model and data files don't exist

using Turing, ContinuousBenchmarks
using Mamba: describe
using FileIO, JLD2

include(joinpath(ContinuousBenchmarks.STAN_MODELS_DIR, "sv.model.jl"))

sv_data = load(joinpath(ContinuousBenchmarks.DATA_DIR, "/nips-2017/sv-data.jld2"))["data"]

# setadbackend(:forward_diff)
# setchunksize(1000)
# chain_nuts = sample(model_f, NUTS(sample_n, 0.65))
# describe(chain_nuts)

# setchunksize(5)
# chain_gibbs = sample(model_f, Gibbs(sample_n, PG(50,1,:h), NUTS(1000,0.65,:ϕ,:σ,:μ)))
# describe(chain_gibbs)

turnprogress(false)

mf = sv_model(data=sv_data[1])
chain_nuts = sample(mf, HMC(2000, 0.05, 10))


text:
using Combinatorics

l = ["iced", "jam", "plain"]
println("List: ", l, "\nCombinations:")
for c in with_replacement_combinations(l, 2)
    println(c)
end

@show length(with_replacement_combinations(1:10, 3))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module Juliaup

include("Tar/src/Tar.jl")

import Downloads, .Tar, CodecZlib, JSON

include("versions_database.jl")

function julia_main()
    try
        real_main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

function get_juliauphome_path()
	# TODO Handle JULIA_DEPOT env var here
	return joinpath(homedir(), ".julia", "juliaup")
end

function get_juliaupconfig_path()
	return joinpath(get_juliauphome_path(), "juliaup.json")
end

function tryparse_channel(value::AbstractString)
	parts = split(value, '.')
	if length(parts)==2 && !any(i->tryparse(Int, i)===nothing, parts)
		return value
	elseif length(parts)==1 && tryparse(Int, value)!==nothing
		return value
	else
		return nothing
	end
end

function tryparse_versionstring(value::String)
	parts = split(value, '~')

	version = nothing
	platform = nothing

	if length(parts)==1
		version = tryparse(VersionNumber, value)
		platform=Int===Int64 ? "x64" : "x86"
	elseif length(parts)==2 && parts[2] in ("x64", "x86")
		version = tryparse(VersionNumber, parts[1])
		platform=parts[2]
	end

	if version===nothing
		return nothing
	else
		return (version=version, platform=platform)
	end
end

function get_julia_versions()
	version_db = get_version_db()

	relevant_versions = filter(version_db) do i
		return i[2]["stable"] && any(i[2]["files"]) do j
			return j["os"]=="winnt" && j["kind"]=="archive"
		end
	end

	typed_versions = map(relevant_versions) do i
		return VersionNumber(i[1])
	end

	sort!(typed_versions)

	return typed_versions
end

function getJuliaVersionsThatMatchChannel(channelString)
	parts = split(channelString, '.')

	versionsThatWeCouldUse = VersionNumber[]

	# Collect all the known versions of Julia that exist that match our channel into a vector
	for currVersion in reverse(get_julia_versions())
		if length(parts) == 1 && parts[1] == string(currVersion.major)
			push!(versionsThatWeCouldUse, currVersion)
		elseif length(parts) == 2 && parts[1] == string(currVersion.major) && parts[2] == string(currVersion.minor)
			push!(versionsThatWeCouldUse, currVersion)
		end
	end

	return versionsThatWeCouldUse
end

function install_version(version::String, config_data::Dict{String,Any})
	if haskey(config_data["InstalledVersions"], version)
		return
	else
		version_db = get_version_db()

		download_url = version_db["AvailableVersions"][version]["Url"]

		version_split = tryparse_versionstring(version)

		target_path = joinpath(get_juliauphome_path(), "julia-$version")

		println("Installing Julia $(version_split.version) ($(version_split.platform)).")

		temp_file = Downloads.download(download_url)

		try
			open(temp_file) do tar_gz
				tar = CodecZlib.GzipDecompressorStream(tar_gz)
				try
					mktempdir() do extract_temp_path
						Tar.extract(tar, extract_temp_path, same_permissions=false)
						folders = readdir(extract_temp_path, join=true)
						if length(folders)!=1
							error("The archive for this version has a folder structure that juliaup does not understand.")
						end
						mv(folders[1], target_path, force=true)
					end
				finally
					close(tar)
				end
			end

			config_data["InstalledVersions"][version] = Dict("Path" => joinpath(".", "julia-$version"))

			println("New version successfully installed.")
		finally
			rm(temp_file, force=true)
		end
	end
end

const g_version_db = Ref{Dict}()

function download_version_db()
	# Disable for now because we are now reliably distributing this file via the store
	# download_url = "https://www.david-anthoff.com/juliaup-versionsdb-winnt-$(Int===Int64 ? "x64" : "x86").json"
	# Downloads.download(download_url, joinpath(get_juliauphome_path(), "juliaup-versionsdb-winnt-$(Int===Int64 ? "x64" : "x86").json"))
end

function get_version_db()
	if !isassigned(g_version_db)
		# For now make the MSIX deployed version file the priority as the store deployment works reliably
		version_db_search_paths = [
			joinpath(Sys.BINDIR, "..", "..", "VersionsDB", "juliaup-versionsdb-winnt-$(Int===Int64 ? "x64" : "x86").json"), # This only exists when MSIX deployed
			joinpath(get_juliauphome_path(), "juliaup-versionsdb-winnt-$(Int===Int64 ? "x64" : "x86").json")
		]
		for i in version_db_search_paths
			if isfile(i)
				g_version_db[] = JSON.parsefile(i, use_mmap=false)
				return g_version_db[]
			end
		end

		error("No version database found.")
	else
		return g_version_db[]
	end
end

function load_config_db()
	juliaup_config_file_path = get_juliaupconfig_path()
	if isfile(juliaup_config_file_path)
		return JSON.parsefile(juliaup_config_file_path, use_mmap=false)
	else
		return Dict{String,Any}("Default"=>"release", "InstalledVersions"=>Dict{String,Any}(), "InstalledChannels"=>Dict{String,Any}())
	end
end

function save_config_db(config_db)
	juliaup_config_file_path = get_juliaupconfig_path()

	open(juliaup_config_file_path, "w") do f
		JSON.print(f, config_db, 4)
	end
end

function is_valid_channel(version_db::Dict{String,Any}, channel::String)
	return haskey(version_db["AvailableChannels"], channel)
end

function get_latest_version_for_channel(channel::String)
	version_db = get_version_db()
	return version_db["AvailableChannels"][channel]["Version"]
end

function garbage_collect_versions(config_data::Dict{String,Any})
	default_channel = config_data["Default"]
	versions_to_uninstall = filter(config_data["InstalledVersions"]) do i
		version = i[1]
		return default_channel!=version && all(config_data["InstalledChannels"]) do j
			if haskey(j[2], "Version")
				return j[2]["Version"] != version
			else
				return true
			end
		end
	end

	for i in versions_to_uninstall
		path_to_delete = joinpath(get_juliauphome_path(), i[2]["Path"])
		try
			rm(path_to_delete, force=true, recursive=true)

			delete!(config_data["InstalledVersions"], i[1])
		catch
			println("WARNING: Failed to delete $path_to_delete. You can try to delete at a later point by running `juliaup gc`.")
		end
	end
end

function update_channel(config_db::Dict{String,Any}, channel::String)
	version_db = get_version_db()

	if version_db["AvailableChannels"][channel]["Version"]!=config_db["InstalledChannels"][channel]["Version"]
		install_version(version_db["AvailableChannels"][channel]["Version"], config_db)

		config_db["InstalledChannels"][channel]["Version"] = version_db["AvailableChannels"][channel]["Version"]
	end
end

function real_main()
    if length(ARGS)==0
        println("Julia Version Manager Preview")
		println()
		println("juliaup command line utility enables configuration of the default Julia version from the command line.")
		println()
		println("usage: juliaup [<command>] [<options>]")
		println()
		println("The following commands are available:")
		println()
		println("  default       Set the default Julia version")
		println("  add           Add a specific Julia version to your system")
		println("  link          Link an existing Julia binary")
		println("  update        Update all or a specific channel to the latest Julia version")
		println("  status        Show all installed Julia versions")
		println("  remove        Remove a Julia version from your system")
		println()
		println("For more details on a specific command, pass it the help argument. [-?] (not yet implemented)")
		println()
		println("The following options are available:")
		println("  -v,--version  Display the version of the tool")
		println("  --info        Display general info of the tool")
		println()
    elseif length(ARGS)>0

		if ARGS[1] == "-v" || ARGS[1] == "--version"
			if length(ARGS)==1
				println("v", JULIA_APP_VERSION)
			else
				println("ERROR: The ", ARGS[1], " argument does not accept any additional arguments.")
			end
		elseif ARGS[1] == "--info"
			if length(ARGS)==1
				println("Julia Version Manager Preview")
				println("Copyright (c) David Anthoff")
			else
				println("ERROR: The --info argument does not accept any additional arguments.")
			end
		elseif ARGS[1] == "default"
			if length(ARGS)==2
				full_channel = ARGS[2]

				data = load_config_db()

				if haskey(data["InstalledChannels"], full_channel)
					data["Default"] = full_channel

					save_config_db(data)

					println("Configured the default Julia version to be '", ARGS[2], "'.")
				else
					println("ERROR: '", ARGS[2], "' is not a valid Julia version.")
				end
			else
				println("ERROR: The default command only accepts one additional argument.")
			end
		elseif ARGS[1] == "add"
			if length(ARGS)==2
				full_channel = ARGS[2]

				download_version_db()

				if is_valid_channel(get_version_db(), full_channel)
					config_data = load_config_db()

					if !haskey(config_data["InstalledChannels"], full_channel)
						required_version = get_latest_version_for_channel(full_channel)

						install_version(required_version, config_data)

						config_data["InstalledChannels"][full_channel] = Dict{String,Any}("Version"=>required_version)

						save_config_db(config_data)
					else
						println("ERROR: '", ARGS[2], "' is already installed.")
					end
				else
					println("ERROR: '", ARGS[2], "' is not a valid Julia version.")
				end
			else
				println("ERROR: The add command only accepts one additional argument.")
			end
		elseif ARGS[1] == "link"
			if length(ARGS)>=3
				channel_name = ARGS[2]
				destination_command = join(ARGS[3:end], ' ')

				config_db = load_config_db()
				version_db = get_version_db()

				if !haskey(config_db["InstalledChannels"], channel_name)
					if haskey(version_db["AvailableChannels"], channel_name)
						println("WARNING: The channel name `$channel_name` is also a system channel. By linking your custom binary to this channel you are hiding this system channel.")
					end

					config_db["InstalledChannels"][channel_name] = Dict{String,Any}("Command"=>destination_command)
				else
					println("ERROR: Channel name `$channel_name` is already used.")
				end

				save_config_db(config_db)
			else
				println("ERROR: The link command requires at least two additional arguments.")
			end
		elseif ARGS[1] == "update" || ARGS[1] == "up"
			if length(ARGS)==1
				download_version_db()

				config_db = load_config_db()
				version_db = get_version_db()

				for i in config_db["InstalledChannels"]
					if haskey(i[2], "Version")
						update_channel(config_db, i[1])
					end
				end

				garbage_collect_versions(config_db)

				save_config_db(config_db)
			elseif length(ARGS)==2
				full_channel = ARGS[2]
				config_db = load_config_db()
				version_db = get_version_db()

				if haskey(config_db["InstalledChannels"], full_channel)
					if haskey(config_db["InstalledChannels"][full_channel], "Version")
						update_channel(config_db, full_channel)

						garbage_collect_versions(config_db)
					else
						println("ERROR: `$full_channel` is a linked channel that cannot be updated.")
					end
				else
					println("Julia $full_channel cannot be updated because it is currently not installed.")
				end

				save_config_db(config_db)
			else
				println("ERROR: The update command accepts at most one additional argument.")
			end
		elseif ARGS[1] == "remove" || ARGS[1] == "rm"

			if length(ARGS)==2
				full_channel = ARGS[2]
				config_data = load_config_db()

				if haskey(config_data["InstalledChannels"], full_channel)
					if full_channel!=config_data["Default"]
						delete!(config_data["InstalledChannels"], full_channel)

						garbage_collect_versions(config_data)

						save_config_db(config_data)

						println("Julia $full_channel successfully removed.")
					else
						println("ERROR: `$(full_channel)` cannot be removed because it is currently configured as the default channel.")
					end
				else
					println("Julia $full_channel cannot be removed because it is currently not installed.")
				end
			else
				println("ERROR: The remove command only accepts one additional argument.")
			end
		elseif ARGS[1] == "status" || ARGS[1] == "st"
			if length(ARGS)==1
				config_data = load_config_db()

				download_version_db()

				version_db = get_version_db()

				defaultJulia = config_data["Default"]

				println("Installed Julia channels (default marked with *):")

				for i in config_data["InstalledChannels"]
					if i[1] == defaultJulia
						print("  * ")
					else
						print("    ")
					end
					print("$(i[1])")

					if haskey(i[2], "Command")
						print(" (linked to `$(i[2]["Command"])`)")
					elseif (version_db["AvailableChannels"][i[1]]["Version"]!=i[2]["Version"])
						print(" (Update from $(i[2]["Version"]) to $(version_db["AvailableChannels"][i[1]]["Version"]) available)")
					end

					println()
				end
			else
				println("ERROR: The status command does not accept any additional arguments.")
			end
		elseif ARGS[1] == "gc"
			if length(ARGS)==1
				config_data = load_config_db()
				garbage_collect_versions(config_data)
			else
				println("ERROR: The gc command accepts no additional argument.")
			end
		else
			println("ERROR: '", ARGS[1], "' is not a recognized command.")
		end
	else
		println("Internal error.")
	end

    return
end

if abspath(PROGRAM_FILE) == @__FILE__
    real_main()
end

end

max_stars_repo_name: piever/juliaup

max_issues_repo_name: piever/juliaup

max_forks_repo_name: piever/juliaup


text:
using Fmincon

# --- function to optimize ----
file = "barnes.jl"
func = "barnes"
gradients = true

# -------- starting point and bounds --------------
x0 = [10.0, 10.0]
ub = [65.0, 70.0]
lb = [0.0, 0.0]

# ---- set options ----
options = Dict(
    "Algorithm" => "active-set",
    "AlwaysHonorConstraints" => "bounds",
    "display" => "iter-detailed",
    "MaxIter" => 1000,
    "MaxFunEvals" => 10000,
    "TolCon" => 1e-6,
    "TolFun" => 1e-6,
    "Diagnostics" => "on")

printfile = "fmincon_summary.out"

# --- run optimization ----
xopt, fopt, exitflag, output = Fmincon.fmincon(file, func, x0, lb, ub,
    options=options, gradients=gradients, printfile=printfile)

# --- print results
@show xopt
@show fopt
@show exitflag
@show output


text:
# Tests defines

# ≈   \approx
# x ≉ y is equivalent to !isapprox(x,y)
# isapprox(x, y; rtol::Real=atol>0 ? 0 : √eps, atol::Real=0, nans::Bool=false, norm::Function)

# but I want my own defs where I know what they are doing

# comon definitions for testing
machine_precision = eps(Float64)
accuracy_goal1 = 1e-12 # default low-precision goal
accuracy_goal2 = 1e-15 # default high-precision goal
relerr(z, x) = z == x ? 0.0 : abs(z - x) / abs(x)
relerrc(z, x) = max(relerr(real(z),real(x)), relerr(imag(z),imag(x)))
function near_equal(x, z, accuracy_goal)
    # for most cases we want small relative errors
    # but near zero, you run the risk of (1.0e-20 - 0.0)/0.0, or (0.0 - 1.0e-20)/1.0e-20
    # also, relative error in (a-b) might be bigger than the actual relative error in the terms
    if abs(z) < 1.0e-4 && abs(x - z) ≤ accuracy_goal
        return true
    else
        return relerr(x,z) ≤ accuracy_goal
    end
end
≈(a,b) = near_equal(a,b,accuracy_goal1)
≅(a,b) = near_equal(a,b,accuracy_goal2)
≒(a,b) = a === b


text:
using Plots, NipponColors

D = dict_color()

x = 0:0.02:1 |> collect
plot(x, sin.(2π .* x), color = D["mizu"])
plot!(x, cos.(2π .* x), color = D["梔子"])


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Replace(old₁ => new₁, old₂ => new₂, ..., oldₙ => newₙ)

Replaces `oldᵢ` value with `newᵢ` value in the table.

# Examples

```julia
Replace(1 => -1, 5 => -5)
Replace(1 => 1.5, 5 => 5.5, 4 => true)
```
"""
struct Replace{K,V} <: Colwise
  pairs::IdDict{K,V}
end

Replace() = throw(ArgumentError("Cannot create a Replace object without arguments."))

Replace(pairs::Pair...) = Replace(IdDict(values(pairs)))

isrevertible(::Type{<:Replace}) = true

function colcache(transform::Replace, x)
  olds = keys(transform.pairs)
  inds = [findall(v -> v === old, x) .=> old for old in olds]
  Dict(reduce(vcat, inds))
end

colapply(transform::Replace, x, c) = [get(transform.pairs, xᵢ, xᵢ) for xᵢ in x]

colrevert(transform::Replace, x, c) = [get(c, i, x[i]) for i in 1:length(x)]


text:
using Distributed

include("misctests_singleprocess.jl")
include("productsplit.jl")
include("paralleltests.jl")

for workersused in [1, 2, 4, 8]
	addprocs(workersused)

	try
		include("paralleltests.jl")
	finally
		rmprocs(workers())
	end
end


text:
using Oceananigans.Architectures: device_event
using Oceananigans.BoundaryConditions: OBC

using Adapt
using KernelAbstractions: @kernel, @index
using Base: @propagate_inbounds

import Oceananigans.BoundaryConditions: fill_halo_regions!
import Statistics: norm, mean, mean!

struct Field{LX, LY, LZ, O, G, T, D, B, S} <: AbstractField{LX, LY, LZ, G, T, 3}
    grid :: G
    data :: D
    boundary_conditions :: B
    operand :: O
    status :: S

    # Inner constructor that does not validate _anything_!
    function Field{LX, LY, LZ}(grid::G, data::D, bcs::B, op::O, status::S) where {LX, LY, LZ, G, D, B, O, S}
        T = eltype(data)
        return new{LX, LY, LZ, O, G, T, D, B, S}(grid, data, bcs, op, status)
    end
end

function validate_field_data(loc, data, grid)
    Tx, Ty, Tz = total_size(loc, grid)

    if size(data) != (Tx, Ty, Tz)
        LX, LY, LZ = loc
        e = "Cannot construct field at ($LX, $LY, $LZ) with size(data)=$(size(data)). " *
            "`data` must have size ($Tx, $Ty, $Tz)."
        throw(ArgumentError(e))
    end

    return nothing
end

validate_boundary_condition_location(bc, ::Center, side) = nothing                  # anything goes for centers
validate_boundary_condition_location(::Union{OBC, Nothing}, ::Face, side) = nothing # only open or nothing on faces
validate_boundary_condition_location(::Nothing, ::Nothing, side) = nothing          # its nothing or nothing
validate_boundary_condition_location(bc, loc, side) = # everything else is wrong!
    throw(ArgumentError("Cannot specify $side boundary condition $bc on a field at $(loc)!"))

validate_boundary_conditions(loc, grid, ::Missing) = nothing
validate_boundary_conditions(loc, grid, ::Nothing) = nothing

function validate_boundary_conditions(loc, grid, bcs)
    sides = (:east, :west, :north, :south, :bottom, :top)
    directions = (1, 1, 2, 2, 3, 3)

    for (side, dir) in zip(sides, directions)
        topo = topology(grid, dir)()
        ℓ = loc[dir]()
        bc = getproperty(bcs, side)

        # Check that boundary condition jives with the grid topology
        validate_boundary_condition_topology(bc, topo, side)

        # Check that boundary condition is valid given field location
        topo isa Bounded && validate_boundary_condition_location(bc, ℓ, side)

        # Check that boundary condition arrays, if used, are on the right architecture
        validate_boundary_condition_architecture(bc, architecture(grid), side)
    end

    return nothing
end

# Common outer constructor for all field flavors that validates data and boundary conditions
function Field(loc::Tuple, grid::AbstractGrid, data, bcs, op, status)
    validate_field_data(loc, data, grid)
    validate_boundary_conditions(loc, grid, bcs)
    LX, LY, LZ = loc
    return Field{LX, LY, LZ}(grid, data, bcs, op, status)
end

#####
##### Vanilla, non-computed Field
#####

"""
    Field{LX, LY, LZ}(grid::AbstractGrid,
                      T::DataType=eltype(grid); kw...) where {LX, LY, LZ}

Construct a `Field` on `grid` with data type `T` at the location `(LX, LY, LZ)`.
Each of `(LX, LY, LZ)` is either `Center` or `Face` and determines the field's
location in `(x, y, z)` respectively.

Keyword arguments
=================

- `data :: OffsetArray`: An offset array with the fields data. If nothing is providet the
  field is filled with zeros.
- `boundary_conditions`: If nothing is provided, then field is created using the default
  boundary conditions via [`FieldBoundaryConditions`](@ref).

Example
=======

```jldoctest
julia> using Oceananigans

julia> ω = Field{Face, Face, Center}(RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)))
1×1×1 Field{Face, Face, Center} on RectilinearGrid on CPU
├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo
├── boundary conditions: west=Periodic, east=Periodic, south=Periodic, north=Periodic, bottom=ZeroFlux, top=ZeroFlux, immersed=ZeroFlux
└── data: 3×3×3 OffsetArray(::Array{Float64, 3}, 0:2, 0:2, 0:2) with eltype Float64 with indices 0:2×0:2×0:2
    └── max=0.0, min=0.0, mean=0.0
```
"""
function Field{LX, LY, LZ}(grid::AbstractGrid,
                           T::DataType=eltype(grid);
                           kw...) where {LX, LY, LZ}

    return Field((LX, LY, LZ), grid, T; kw...)
end

function Field(loc::Tuple,
               grid::AbstractGrid,
               T::DataType = eltype(grid);
               data = new_data(T, grid, loc),
               boundary_conditions = FieldBoundaryConditions(grid, loc))

    return Field(loc, grid, data, boundary_conditions, nothing, nothing)
end

Field(f::Field) = f # indeed
Field(z::ZeroField) = z

#####
##### Field utils
#####

# Canonical `similar` for Field (doesn't transfer boundary conditions)
function Base.similar(f::Field, grid=f.grid)
    loc = location(f)
    return Field(loc,
                 grid,
                 new_data(eltype(parent(f)), grid, loc),
                 FieldBoundaryConditions(grid, loc),
                 f.operand,
                 deepcopy(f.status))
end

# Fallback: cannot infer boundary conditions.
boundary_conditions(field) = nothing
boundary_conditions(f::Field) = f.boundary_conditions

function interior(a::Union{Field, OffsetArray}, (LX, LY, LZ), grid)
    TX, TY, TZ = topology(grid)
    ii = interior_parent_indices(LX, TX, grid.Nx, grid.Hx)
    jj = interior_parent_indices(LY, TY, grid.Ny, grid.Hy)
    kk = interior_parent_indices(LZ, TZ, grid.Nz, grid.Hz)
    return view(parent(a), ii, jj, kk)
end

"Returns a view of `f` that excludes halo points."
interior(f::Field) = interior(f, location(f), f.grid)

function interior_copy(f::Field)
    LX, LY, LZ = location(f)
    return parent(f)[interior_parent_indices(LX, topology(f, 1), f.grid.Nx, f.grid.Hx),
                     interior_parent_indices(LY, topology(f, 2), f.grid.Ny, f.grid.Hy),
                     interior_parent_indices(LZ, topology(f, 3), f.grid.Nz, f.grid.Hz)]
end

# Don't use axes(f) to checkbounds; use axes(f.data)
Base.checkbounds(f::Field, I...) = Base.checkbounds(f.data, I...)

@propagate_inbounds Base.getindex(f::Field, inds...) = getindex(f.data, inds...)
@propagate_inbounds Base.getindex(f::Field, i::Int)  = parent(f)[i]
@propagate_inbounds Base.setindex!(f::Field, val, i, j, k) = setindex!(f.data, val, i, j, k)
@propagate_inbounds Base.lastindex(f::Field) = lastindex(f.data)
@propagate_inbounds Base.lastindex(f::Field, dim) = lastindex(f.data, dim)

Base.fill!(f::Field, val) = fill!(parent(f), val)
Base.parent(f::Field) = parent(f.data)
Adapt.adapt_structure(to, f::Field) = Adapt.adapt(to, f.data)

data(f::Field) = f.data
cpudata(a::Field) = arch_array(CPU(), a.data)

#####
##### Special constructors for tracers and velocity fields
#####

"""
    CenterField(grid; kwargs...)

Returns `Field{Center, Center, Center}` on `arch`itecture and `grid`.
Additional keyword arguments are passed to the `Field` constructor.
"""
CenterField(grid::AbstractGrid, T::DataType=eltype(grid); kw...) = Field{Center, Center, Center}(grid, T; kw...)

"""
    XFaceField(grid; kw...)

Returns `Field{Face, Center, Center}` on `grid`.
Additional keyword arguments are passed to the `Field` constructor.
"""
XFaceField(grid::AbstractGrid, T::DataType=eltype(grid); kw...) = Field{Face, Center, Center}(grid, T; kw...)

"""
    YFaceField(grid; kw...)

Returns `Field{Center, Face, Center}` on `grid`.
Additional keyword arguments are passed to the `Field` constructor.
"""
YFaceField(grid::AbstractGrid, T::DataType=eltype(grid); kw...) = Field{Center, Face, Center}(grid, T; kw...)

"""
    ZFaceField(grid; kw...)

Returns `Field{Center, Center, Face}` on `grid`.
Additional keyword arguments are passed to the `Field` constructor.
"""
ZFaceField(grid::AbstractGrid, T::DataType=eltype(grid); kw...) = Field{Center, Center, Face}(grid, T; kw...)

#####
##### Interface for field computations
#####

"""
    compute!(field)

Computes `field.data` from `field.operand`.
"""
compute!(field, time=nothing) = nothing # fallback

"""
    @compute(exprs...)

Call `compute!` on fields after defining them.
"""
macro compute(def)
    expr = Expr(:block)
    field = def.args[1]
    push!(expr.args, :($(esc(def))))
    push!(expr.args, :(compute!($(esc(field)))))
    return expr
end

# Computation "status" for avoiding unnecessary recomputation
mutable struct FieldStatus{T}
    time :: T
end

FieldStatus() = FieldStatus(0.0)
Adapt.adapt_structure(to, status::FieldStatus) = (; time = status.time)

"""
    compute_at!(field, time)

Computes `field.data` at `time`. Falls back to compute!(field).
"""
compute_at!(field, time) = compute!(field)

"""
    compute_at!(field, time)

Computes `field.data` if `time != field.status.time`.
"""
function compute_at!(field::Field, time)
    if isnothing(field.status) # then always compute:
        compute!(field, time)

    # Otherwise, compute only on initialization or if field.status.time is not current,
    elseif time == zero(time) || time != field.status.time
        compute!(field, time)
        field.status.time = time
    end

    return field
end

# This edge case occurs if `fetch_output` is called with `model::Nothing`.
# We do the safe thing here and always compute.
compute_at!(field::Field, ::Nothing) = compute!(field, nothing)

#####
##### Fields that are reduced along one or more dimensions
#####

const XReducedField = Field{Nothing}
const YReducedField = Field{<:Any, Nothing}
const ZReducedField = Field{<:Any, <:Any, Nothing}

const YZReducedField = Field{<:Any, Nothing, Nothing}
const XZReducedField = Field{Nothing, <:Any, Nothing}
const XYReducedField = Field{Nothing, Nothing, <:Any}

const XYZReducedField = Field{Nothing, Nothing, Nothing}

const ReducedField = Union{XReducedField, YReducedField, ZReducedField,
                           YZReducedField, XZReducedField, XYReducedField,
                           XYZReducedField}

reduced_dimensions(field::Field) = ()

reduced_dimensions(field::XReducedField) = tuple(1)
reduced_dimensions(field::YReducedField) = tuple(2)
reduced_dimensions(field::ZReducedField) = tuple(3)

reduced_dimensions(field::YZReducedField) = (2, 3)
reduced_dimensions(field::XZReducedField) = (1, 3)
reduced_dimensions(field::XYReducedField) = (1, 2)

reduced_dimensions(field::XYZReducedField) = (1, 2, 3)

function fill_halo_regions!(field::Field, arch, args...; kwargs...)
    reduced_dims = reduced_dimensions(field)
    if reduced_dimensions === () # the field is not reduced!
        return fill_halo_regions!(field.data,
                                  field.boundary_conditions,
                                  architecture(field),
                                  field.grid,
                                  args...;
                                  kwargs...)
    else
        return fill_halo_regions!(field.data,
                                  field.boundary_conditions,
                                  architecture(field),
                                  field.grid,
                                  args...;
                                  reduced_dimensions=reduced_dims,
                                  kwargs...)
    end
end

@propagate_inbounds Base.getindex(r::XReducedField, i, j, k) = getindex(r.data, 1, j, k)
@propagate_inbounds Base.getindex(r::YReducedField, i, j, k) = getindex(r.data, i, 1, k)
@propagate_inbounds Base.getindex(r::ZReducedField, i, j, k) = getindex(r.data, i, j, 1)

@propagate_inbounds Base.setindex!(r::XReducedField, v, i, j, k) = setindex!(r.data, v, 1, j, k)
@propagate_inbounds Base.setindex!(r::YReducedField, v, i, j, k) = setindex!(r.data, v, i, 1, k)
@propagate_inbounds Base.setindex!(r::ZReducedField, v, i, j, k) = setindex!(r.data, v, i, j, 1)

@propagate_inbounds Base.getindex(r::YZReducedField, i, j, k) = getindex(r.data, i, 1, 1)
@propagate_inbounds Base.getindex(r::XZReducedField, i, j, k) = getindex(r.data, 1, j, 1)
@propagate_inbounds Base.getindex(r::XYReducedField, i, j, k) = getindex(r.data, 1, 1, k)

@propagate_inbounds Base.setindex!(r::YZReducedField, v, i, j, k) = setindex!(r.data, v, i, 1, 1)
@propagate_inbounds Base.setindex!(r::XZReducedField, v, i, j, k) = setindex!(r.data, v, 1, j, 1)
@propagate_inbounds Base.setindex!(r::XYReducedField, v, i, j, k) = setindex!(r.data, v, 1, 1, k)

@propagate_inbounds Base.getindex(r::XYZReducedField, i, j, k) = getindex(r.data, 1, 1, 1)
@propagate_inbounds Base.setindex!(r::XYZReducedField, v, i, j, k) = setindex!(r.data, v, 1, 1, 1)

# Preserve location when adapting fields reduced on one or more dimensions
function Adapt.adapt_structure(to, reduced_field::ReducedField)
    LX, LY, LZ = location(reduced_field)
    return Field{LX, LY, LZ}(nothing,
                             adapt(to, reduced_field.data),
                             nothing,
                             nothing,
                             nothing)
end

#####
##### Field reductions
#####

# TODO: needs test
Statistics.dot(a::Field, b::Field) = mapreduce((x, y) -> x * y, +, interior(a), interior(b))

# TODO: in-place allocations with function mappings need to be fixed in Julia Base...
const SumReduction     = typeof(Base.sum!)
const MeanReduction    = typeof(Statistics.mean!)
const ProdReduction    = typeof(Base.prod!)
const MaximumReduction = typeof(Base.maximum!)
const MinimumReduction = typeof(Base.minimum!)
const AllReduction     = typeof(Base.all!)
const AnyReduction     = typeof(Base.any!)

initialize_reduced_field!(::SumReduction,  f, r::ReducedField, c) = Base.initarray!(interior(r), Base.add_sum, true, interior(c))
initialize_reduced_field!(::ProdReduction, f, r::ReducedField, c) = Base.initarray!(interior(r), Base.mul_prod, true, interior(c))
initialize_reduced_field!(::AllReduction,  f, r::ReducedField, c) = Base.initarray!(interior(r), &, true, interior(c))
initialize_reduced_field!(::AnyReduction,  f, r::ReducedField, c) = Base.initarray!(interior(r), |, true, interior(c))

initialize_reduced_field!(::MaximumReduction, f, r::ReducedField, c) = Base.mapfirst!(f, interior(r), c)
initialize_reduced_field!(::MinimumReduction, f, r::ReducedField, c) = Base.mapfirst!(f, interior(r), c)

filltype(f, c) = eltype(c)
filltype(::Union{AllReduction, AnyReduction}, grid) = Bool

function reduced_location(loc; dims)
    if dims isa Colon
        return (Nothing, Nothing, Nothing)
    else
        return Tuple(i ∈ dims ? Nothing : loc[i] for i in 1:3)
    end
end

function reduced_dimension(loc)
    dims = ()
    for i in 1:3
        loc[i] == Nothing ? dims = (dims..., i) : dims
    end
    return dims
end

## Allow support for ConditionalOperation

get_neutral_mask(::Union{AllReduction, AnyReduction})  = true
get_neutral_mask(::Union{SumReduction, MeanReduction}) =   0
get_neutral_mask(::MinimumReduction) =   Inf
get_neutral_mask(::MaximumReduction) = - Inf
get_neutral_mask(::ProdReduction)    =   1

# If func = identity and condition = nothing, nothing happens
@inline condition_operand(f::typeof(identity), operand::AbstractField, ::Nothing, mask) = operand

@inline condition_operand(operand, condition, mask) = condition_operand(identity, operand, condition, mask)

@inline conditional_length(c::AbstractField)        = length(c)
@inline conditional_length(c::AbstractField, dims)  = mapreduce(i -> size(c, i), *, unique(dims); init=1)

# Allocating and in-place reductions
for reduction in (:sum, :maximum, :minimum, :all, :any, :prod)

    reduction! = Symbol(reduction, '!')

    @eval begin

        # In-place
        Base.$(reduction!)(f::Function, r::ReducedField, a::AbstractArray;
                           condition = nothing, mask = get_neutral_mask(Base.$(reduction!)), kwargs...) =
            Base.$(reduction!)(identity, interior(r), condition_operand(f, a, condition, mask); kwargs...)

        Base.$(reduction!)(r::ReducedField, a::AbstractArray;
                           condition = nothing, mask = get_neutral_mask(Base.$(reduction!)), kwargs...) =
            Base.$(reduction!)(identity, interior(r), condition_operand(a, condition, mask); kwargs...)

        # Allocating
        function Base.$(reduction)(f::Function, c::AbstractField;
                                   condition = nothing, mask = get_neutral_mask(Base.$(reduction!)),
                                   dims=:)
            if dims isa Colon
                r = zeros(architecture(c), c.grid, 1, 1, 1)
                Base.$(reduction!)(identity, r, condition_operand(f, c, condition, mask))
                return CUDA.@allowscalar r[1, 1, 1]
            else
                T = filltype(Base.$(reduction!), c)
                loc = reduced_location(location(c); dims)
                r = Field(loc, c.grid, T)
                initialize_reduced_field!(Base.$(reduction!), identity, r, condition_operand(f, c, condition, mask))
                Base.$(reduction!)(identity, r, condition_operand(f, c, condition, mask), init=false)
                return r
            end
        end

        Base.$(reduction)(c::AbstractField; kwargs...) = Base.$(reduction)(identity, c; kwargs...)
    end
end

function Statistics._mean(f, c::AbstractField, ::Colon; condition = nothing, mask = 0)
    operator = condition_operand(f, c, condition, mask)
    return sum(operator) / conditional_length(operator)
end

function Statistics._mean(f, c::AbstractField, dims; condition = nothing, mask = 0)
    operator = condition_operand(f, c, condition, mask)
    r = sum(operator; dims)
    n = conditional_length(operator, dims)
    r ./= n
    return r
end

Statistics.mean(f::Function, c::AbstractField; condition = nothing, dims=:) = Statistics._mean(f, c, dims; condition)
Statistics.mean(c::AbstractField; condition = nothing, dims=:) = Statistics._mean(identity, c, dims; condition)

function Statistics.mean!(f::Function, r::ReducedField, a::AbstractArray; condition = nothing, mask = 0)
    sum!(f, r, a; condition, mask, init=true)
    dims = reduced_dimension(location(r))
    n = conditional_length(condition_operand(f, a, condition, mask), dims)
    r ./= n
    return r
end

Statistics.mean!(r::ReducedField, a::AbstractArray; kwargs...) = Statistics.mean!(identity, r, a; kwargs...)

function Statistics.norm(a::AbstractField; condition = nothing)
    r = zeros(a.grid, 1)
    Base.mapreducedim!(x -> x * x, +, r, condition_operand(a, condition, 0))
    return CUDA.@allowscalar sqrt(r[1])
end

function Base.isapprox(a::AbstractField, b::AbstractField; kw...)
    conditioned_a = condition_operand(a, nothing, one(eltype(a)))
    conditioned_b = condition_operand(b, nothing, one(eltype(b)))
    # TODO: Make this non-allocating?
    return all(isapprox.(conditioned_a, conditioned_b; kw...))
end


text:
module MyExample
using ForwardDiff

# Write your package code here.
include("extra_file.jl")
export my_f, derivative_of_my_f

end


text:
"""
This library provides modal functions and pattern matching in Julia.

The design is based on the Ivo programming language.
Unlike Ivo, the compilation strategy follows the JMatch design
using resumable functions (aka iterators) rather than a LogicT-like monad.
This turned out to be simpler to implement using Julia macros
(especially since there's already a resumable function package).

The pattern matching implementation is based on Rematch.jl.
"""
module Logic

import MacroTools
import MacroTools: @capture

# Error when a formula has no solutions.
struct MatchFailure <: Exception
  value
end

# Error when a formula has more than one solution.
struct MatchAmbiguous <: Exception
  value
end

# How to translate a function definition with patterns?
#
#

# How to translate a relation definition?
#
# @rel +(in x, out y, in z)
#   y = z - x
#   @yield
# end
#
# -->
#
# @resumable
# function +(x, y::Out, z)
#   y = z - x
#   @yield y
# end
#
# Another example:
#
# @rel +(out x, out y, in z::T)
#   x = typemin(T)
#
#   while true
#     y = z - x
#     @yield
#
#     if x == typemax(T)
#       return
#     end
#
#     x += 1
#   end
# end
#
# @let 1+x = 2
# ->
# x = +(1, Out, 2)

abstract type Mode end
struct In <: Mode end
struct Out <: Mode end

include("prims.jl")

# Returns true if the pattern has unknowns.
# If not, it can just be evaluated in forward mode.
function has_unknowns(pat, bound)::Bool
  # println("pat $pat")
  # println("bound $bound")

  if @capture(pat, p_ = e_)
    has_unknowns(p, bound)
  elseif @capture(pat, p_ in e_)
    has_unknowns(p, bound)
  elseif pat == :(_)
    true
  elseif !(pat isa Expr || pat isa Symbol)
    # constants (including strings, true and false)
    false
  elseif pat == :nothing || pat == :missing
    # nothing
    false
  elseif @capture(pat, _quote_macrocall)
    false
  elseif @capture(pat, Symbol(_))
    false
  elseif @capture(pat, x_Symbol)
    ! (x in bound)
  elseif @capture(pat, p1_ || p2_)
    has_unknowns(p1, bound) || has_unknowns(p2, bound)
  elseif @capture(pat, p1_ && p2_)
    has_unknowns(p1, bound) || has_unknowns(p2, bound)
  elseif @capture(pat, (p1_ where p2_))
    has_unknowns(p1, bound) || has_unknowns(p2, bound)
  elseif @capture(pat, T_Symbol(ps__))
    # call (on a symbol)
    # there are some special cases
    # if T is + or *, 1+2+3+4+5 has the tree +(1,2,3,4,5)
    # 5x is *(5,:x)
    # T might be a struct name. Then the pattern is Foo(x=2, y=3)
    Base.any(p -> has_unknowns(p, bound), ps)
  elseif @capture(pat, (ps__,))
    # tuple
    Base.any(p -> has_unknowns(p, bound), ps)
  elseif @capture(pat, [ps__,])
    # array
    Base.any(p -> has_unknowns(p, bound), ps)
  elseif @capture(pat, [ps__;])
    # array
    Base.any(p -> has_unknowns(p, bound), ps)
  elseif @capture(pat, p_::T_)
    # type assert
    has_unknowns(p, bound)
  else
    # anything else must be evaluated in forward mode
    # T{a}
    # a[1,2,3]
    # if, while, for, etc
    # ? :
    # call on a lambda
    # lambda
    false
  end
end

function match_tuple(ps, v, i, body, bound, rename)
  println("ps $ps")
  if i > length(ps)
    body
  else
    @assert 1 <= i && i <= length(ps)
    # FIXME binds variables in the wrong order
    rest = match_tuple(ps, v, i+1, body, bound, rename)
    translate_match(ps[i], :($v[$i]), rest, bound, rename)
  end
end

# Some of this code is based on Rematch.jl.
# But rather than return a boolean, we eval the body for each
# solution.
function translate_match(pat, value, body, bound::Set{Symbol}, rename::Dict{Symbol,Symbol})
  if !has_unknowns(pat, bound)
    # constants and other forward-mode expressions
    quote
      if $pat == $value
        $body
      end
    end
  elseif @capture(pat, p_ = e_)
    # the pattern is a bind formula. the value should be true
    # note, unlike in a true relational language, we only bind
    # unknowns on the left; the term on the right is not a pattern
    body1 = translate_match(p, e, body, bound, rename)

    if value == true
      body1
    elseif value == false
      :(throw(MatchFailure(value)))
    else
      quote
        if $value
          $body1
        end
      end
    end
  elseif @capture(pat, p_ in e_)
    # the pattern is an `in` pattern (`<-` in Ivo).
    tmp = gensym()
    quote
      if $value
        for $tmp in $e
          $(translate_match(p, tmp, body, bound, rename))
        end
      else
        nothing
      end
    end
  elseif pat == :(_)
    tmp = gensym()
    quote
      let $tmp = $value;
        $body
      end
    end
  elseif @capture(pat, x_Symbol)
    # variable

    # if the pattern doesn't match, we don't want to set the variable,
    # so use a temporary

    # FIXME: don't rename yet
    tmp = x # get!(rename, x, gensym(x))

    if x in bound
      # already bound; eval in forward mode
      quote
        if $tmp == $value
          $body
        end
      end
    else
      # println("$pat :: bound $bound + $x")
      push!(bound, x)
      quote
        let $tmp = $value;
          $body
        end
      end
    end
  elseif @capture(pat, p1_ || p2_)
    matched = gensym()
    v = gensym()

    bound1 = copy(bound)
    bound2 = copy(bound)

    body1 = translate_match(p1, v, :($matched = true), bound1, rename)
    body2 = translate_match(p2, v, :($matched = true), bound2, rename)

    union!(bound, intersect(bound1, bound2))

    quote
      let $v = $value;
        let $matched = false;
          begin
            $body1;
            if !$matched
              $body2
            end
            if $matched
              $body
            end
          end
        end
      end
    end
  elseif @capture(pat, p1_ && p2_)
    v = gensym()

    # HACK: compute the bound from p1 and use that for p2, but then redo with the original bound
    # have to do this because we pass the body in and bind in p2 first
    bound1 = copy(bound)
    translate_match(p1, v, :(), bound1, rename)

    body1 = translate_match(p1, v, translate_match(p2, v, body, bound1, rename), bound, rename)
    union!(bound, bound1)

    quote
      let $v = $value;
        $body1
      end
    end
  elseif @capture(pat, (p_ where guard_))
    # both p and guard can bind variables
    #
    # similar to &&, but matches the guard against true

    v = gensym()

    # HACK: compute the bound from p1 and use that for p2, but then redo with the original bound
    # have to do this because we pass the body in and bind in p2 first
    bound1 = copy(bound)
    translate_match(p, v, :(), bound1, rename)

    body1 = translate_match(p, v, translate_match(guard, true, body, bound1, rename), bound, rename)
    union!(bound, bound1)

    quote
      let $v = $value;
        $body1
      end
    end
  elseif @capture(pat, f_(ps__))
    result = gensym()

    begin
      # FIXME bound is wrong because of unknowns bound in earlier patterns are *not* unknowns in later patterns
      outputs = filter(p -> has_unknowns(p, bound), ps)

      inputs = map(p -> has_unknowns(p, bound) ? :(Out()) : p, ps)
      newcall = Expr(:call, f, value, inputs...)

      # println("pat ", pat)
      # println("value ", value)
      # println("f ", f)
      # println("ps ", ps)
      # println("inputs ", inputs)
      # println("outputs ", outputs)
      # println("newcall ", newcall)
      # println("newcall.head ", newcall.head)
      # println("newcall.args ", newcall.args)
      # for (i,a) in enumerate(newcall.args)
      #   println("newcall.args[$i] $a")
      # end

      if length(outputs) == 1
        quote
          for $result in $newcall
            $(translate_match(outputs[1], result, body, bound, rename))
          end
        end
      else
        results_pat = Expr(:tuple, outputs...)
        quote
          for $result in $newcall
            $(translate_match(results_pat, result, body, bound, rename))
          end
        end
      end
    end

    #=
    # struct or call
    # TODO TODO TODO
    #
    # eval the function, inject into monad if not overloaded
    # eval the input arguments
    # eval $value
    # select-thunk(apply-amb (f >>= mode-select [invert-mode ps]))

    (apply-amb u x) = u >>= \g -> g x   [g fresh]
    (apply-amb u x1 x2..xn) = (apply-amb (apply-amb u x1)
                                         x2..xn)

    select-thunk (amb ...)
        -- filter: if exactly one unique, done
                   if all normal, done
                   if all default, done
                   if normal and unique, error
                   if default and unique, throw out the default
                   if default and normal, throw out the default

        -- if there's one, choose it
        -- if there's more than one, MatchAmbiguous
        -- if there's zero, MatchFailure
        --

        =#
  elseif @capture(pat, (ps__,))
    # tuple
    v = gensym()

    quote
      let $v = $value
        if isa($v, Tuple)
          if length($v) == $(length(ps))
            $(match_tuple(ps, v, 1, body, bound, rename))
          end
        end
      end
    end
  elseif @capture(pat, [ps__,])
    # array
    v = gensym()
    w = gensym()
    quote
      let $v = $value
        if isa($v, AbstractArray)
          if length($v) == $(length(ps))
            $(match_tuple(ps, v, 1, body, bound, rename))
          end
        end
      end
    end
  elseif @capture(pat, p_::T_)
    # type assert
    v = gensym()
    quote
      let $v = $value
        if isa($v, $(esc(T)))
          $(translate_match(p, v, body, bound, rename))
        end
      end
    end
  else
    error("Unrecognized pattern syntax: $pat")
  end
end

"""
    @match value begin
        pattern1 => result1
        pattern2 => result2
        ...
    end

Return `result` for the first matching `pattern`.
If there are no matches, throw `MatchFailure`.
If there are multiple matches, throw `MatchAmbiguous`.
"""
macro match(value, body)
  println(value)
  println(body)
  bound1 = Set{Symbol}()
  bound2 = Set{Symbol}()
  rename1 = Dict{Symbol,Symbol}()
  rename2 = Dict{Symbol,Symbol}()

  if @capture(body, begin cases__ end)
    input = gensym("value")
    output = gensym("result")
    count = gensym("count")

    accum = []
    counter = []

    for case in cases
      if @capture(case, pat_ => exp_)
        result = translate_match(pat, input, :true, bound1, rename1)
        push!(counter, result)

        result = translate_match(pat, input, exp, bound2, rename2)
        push!(accum, result)
      else
        @error("bad match case")
      end
    end

    quote
      begin
        $input = $value
        $output = nothing
        $count = 0
        begin
            $((:($count = $count + (isnothing($s) ? 0 : 1)) for s in counter)...)
        end
        if $count == 0
          throw(MatchFailure($input))
        elseif $count > 1
          throw(MatchAmbiguous($input))
        else
          $((:($output = $s) for s in accum)...)
        end
        $output
      end
    end
  else
    @error("bad match statement")
  end
end

"""
    @foreach(formula)
        ...
    end

Iterate over all solutions of the formula. Unknowns in the formula are bound in the body.
"""
macro foreach(formula, body)
  println(formula)
  println(body)
  bound = Set{Symbol}()
  rename = Dict{Symbol,Symbol}()
  x = translate_match(formula, :true, body, bound, rename)
  println(x)
  x
end

"""
    @once(formula)
        ...
    end

Iterate over all solutions of the formula. Unknowns in the formula are bound in the body.
"""
macro once(formula, body)
  println(formula)
  println(body)
  bound = Set{Symbol}()
  rename = Dict{Symbol,Symbol}()
  x = translate_match(formula, :true, :(begin; $body; break; end), bound, rename)
  println(x)
  quote
    for _ in [()]
      $x
    end
  end
end

"""
    @define(formula)

    @let formula begin
        body
    end

If `formula` has a solution, bind variables and return `true`.
Throws an error if there is not exactly one solution.
If there are no solutions, throw `MatchFailure`.
If there is more than one solutions, throw `MatchAmbiguous`.
"""
macro define(formula)
  #=
    quote
        begin
            apply_thunk(unwrap($translate_match(formula, :true, body)))
        end
    end
    =#
end

export MatchFailure, MatchAmbiguous
export @match
export @foreach, @once
export Out

end #module


text:
# Module
module EngThermBase

# Imports
using Reexport
@reexport using Unitful
@reexport using Measurements

# Includes - abstract supertypes
include("abstract.jl")

# Includes - default settings
include("settings.jl")

# Includes - type exactness
include("exactness.jl")

# Includes - concrete types
include("amounts.jl")
include("operations.jl")
include("constants.jl")
include("states.jl")

# Includes - functionalities
include("utils.jl")

# Module
end


text:
# dynamic solution
function maxpathsum(t::Array{Array{I, 1}, 1}) where I
    T = deepcopy(t)
    for r in length(T)-1:-1:1
        for c in linearindices(T[r])
            T[r][c] += max(T[r+1][c], T[r+1][c+1])
        end
    end
    return T[1][1]
end

test = [[55],
        [94, 48],
        [95, 30, 96],
        [77, 71, 26, 67],
        [97, 13, 76, 38, 45],
        [07, 36, 79, 16, 37, 68],
        [48, 07, 09, 18, 70, 26, 06],
        [18, 72, 79, 46, 59, 79, 29, 90],
        [20, 76, 87, 11, 32, 07, 07, 49, 18],
        [27, 83, 58, 35, 71, 11, 25, 57, 29, 85],
        [14, 64, 36, 96, 27, 11, 58, 56, 92, 18, 55],
        [02, 90, 03, 60, 48, 49, 41, 46, 33, 36, 47, 23],
        [92, 50, 48, 02, 36, 59, 42, 79, 72, 20, 82, 77, 42],
        [56, 78, 38, 80, 39, 75, 02, 71, 66, 66, 01, 03, 55, 72],
        [44, 25, 67, 84, 71, 67, 11, 61, 40, 57, 58, 89, 40, 56, 36],
        [85, 32, 25, 85, 57, 48, 84, 35, 47, 62, 17, 01, 01, 99, 89, 52],
        [06, 71, 28, 75, 94, 48, 37, 10, 23, 51, 06, 48, 53, 18, 74, 98, 15],
        [27, 02, 92, 23, 08, 71, 76, 84, 15, 52, 92, 63, 81, 10, 44, 10, 69, 93]]

@show maxpathsum(test)


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

using Markdown
using InteractiveUtils

# ╔═╡ 47b82739-35ba-496c-b5e0-2e73e59cf1c4
begin
	import Pkg
	Pkg.activate(mktempdir())
end

# ╔═╡ 1b7a3fe9-8716-427f-b796-ca94f7e48a09
begin
	Pkg.add(url="https://github.com/fyzycyst/WebMSISE00.jl.git")
	using WebMSISE00
end

# ╔═╡ d8961618-08aa-4424-9cba-960c5023d1fb
begin
	Pkg.add("Plots")
	Pkg.add("StatPlots")
	Pkg.add("DataFrames")
	using DataFrames
	using Plots
	using StatPlots
end

# ╔═╡ 28dc8709-edc3-4754-a2dd-fc98d0890469
using Plots.PlotMeasures

# ╔═╡ 9ba46bb7-66de-4756-ba4e-613cb1e6f751
md"""
# Example use of _WebMSISE00.jl_

This notebook just has some quick & dirty examples of generating an atmosphere and
plotting the density profile(s) vs. altitude.
"""

# ╔═╡ db0d3958-a55a-4d49-b27e-c586007c8676
md"""

_First, we load all the packages needed._

In addition to `WebMSISE00`, we load in `Plots`, `StatPlots`, and `DataFrames`, since
we're going to plot from the returned dataframe.
"""

# ╔═╡ dabf432d-434a-4526-a61e-e4d2d61d77d3
plotly()

# ╔═╡ 0d766211-5fd1-4810-b799-5c0f28e088e8
md"""

_Next, we set up and run an atmosphere profile._

"""

# ╔═╡ 177978c2-935d-432e-b5e3-c623d5f6b492
md"""
In order to do a run, we need to first discuss the inputs to the model.
"""

# ╔═╡ 234718cf-381a-4d06-8403-67d993922b50
md"""
The inputs are all entered as strings, even numerical inputs. The important inputs (entered as a Dict) are:\
_[NOTE: valid dates are 1960/02/14 - 2020/09/16, with the end date updated monthly]_\
_[NOTE 2: If none of the solar activity indices (F10.7 or Ap) are entered, the real value for the date specified is automatically used.]_
- `"year" => "2020"`,              (4 digits)
- `"month" => "01"`,               (2 digits)
- `"day" => "01"`,                 (2 digits)
- `"time_flag" => "0"`,            ('0' for Universal Time; '1' for local time)
- `"hour" => "1.5"`,               (hours after midnight)
- `"latitude" => "55."`,           (from -90. to 90.)
- `"longitude" => "45."`,          (from 0. to 360.)
- `"profile" => "1"`,              (defines the parameter sweep performed)
    * '1' = Height (from 0. to 1000.)
    * '2' = Latitude (from -90. to 90.)
    * '3' = Longitude (from 0. to 360.)
    * '4' = Year (from 1961. to 2007.)
    * '5' = Month (from 1. to 12.)
    * '6' = Day of Month (from 1. to 31.)
    * '7' = Day of Year (from 1. to 366.)
    * '8' = Hour of Day (from 1. to 24.)
- `"start" => "0."`,                (start of the parameter sweep)
- `"stop" => "1000."`,              (end of the parameter sweep)
- `"step" => "50."`,                (step size for the parameter sweep)
- `"f10_7" => ""`,                  (F10.7 daily value)
- `"f10_7_3" => ""`,                (F10.7 three month average)
- `"ap" => ""`,                     (Ap daily value)
- `"vars" => ["05", "08", "09"]`    (lists the independent & dependent variables)

"vars" determines the columns written to the output datafile which ends up loaded into the DataFrame output.\
 \
The `vars` Array values that are exposed in `WebMSISE00` are:\
_Independent variables:_ \
00 = Year; 01 = Month; 02 = Day of Month; 03 = Day of Year;\
04 = Hour of Day; 05 = Height (km); 06 = Latitude; 07 = Longitude\
_Calculated Parameters:_ \
08 = Atomic Oxygen number density (cm^-3); 09 = N2 number density;\
10 = O2 number density; 11 = total mass density (g/cm^3);\
12 = Neutral Temp (K); 13 = Exospheric Temp (K); 14 = He number density;\
15 = Ar number density; 16 = H number density; 17 = N number density;\
18 = Anomalous Oxygen number density\
_Indices taken from user input vs. data base_ \
19 = F10.7 (daily); 20 = F10.7 (3 month avg); 21 = ap (daily)\
 \
Finally, the values listed in the above example Dict entries are the defaults used by `WebMSISE00` for any not explicitly set by the user in the input Dict.
"""

# ╔═╡ d56a1b81-1d17-463b-b5d1-751cf531da76
md"""
As an initial example, let's choose a time of low solar activity... 1 Jan 2020, near dawn (hour 5, local time), for a mid latitude of $40^{o}$, and an F10.3 three month average value of 70. In addition, let's get a profile for the total density *and* all the individual constituents, from ground to 1000 km (this range is a default setting), at a step size of 25 km.
"""

# ╔═╡ 50fc3dd4-fc7b-4283-804f-053c29edbe72
input_low = Dict(
	"year" => "2020",
	"month" => "01",
	"day" => "01",
	"time_flag" => "1",
	"hour" => "5",
	"latitude" => "40.",
	"step" => "25.",
	"f10_7_3" => "70.",
	"vars" => ["05", "08", "09","10","11","14","15","16","17","18"]
	)

# ╔═╡ a9c5f4df-6c71-449b-8ae9-8195c5b62357
output_low = ReadMSISE(input_low)

# ╔═╡ 150bacf2-b05e-42cf-957a-46df242d4d3a
begin
	@df output_low plot(cols(1), cols(2:4),
		yaxis=:log10,
		size=(650,600),
		ylims=(1e-19,1e20),
		label=["O" "N2" "O2"],
		legend=(0.8,0.95),
		title="Density of atmospheric constituents with altitude",
		xlabel="Altitude in km",
		ylabel="Number density in cm-3"
	)
	@df output_low plot!(cols(1), cols(6:10),
		label=["He" "Ar" "H" "N" "Anom. O"]
	)
end

# ╔═╡ 189a1664-bc02-40f0-8347-c5725426faa4
md"""

_Now, compare two days at different parts of the solar cycle to exhibit how the
atmosphere gets 'puffy' when the Sun is particularly active._

"""

# ╔═╡ b4e055c8-cedc-4aef-b170-f51f5afe17d3
md"""
We already have the total atmosphere density (for the low activity portion of the solar cycle) in the `output_low` dataframe, so now we generate the density for high solar activity. We'll just generate the total density without the number densities of all the constituents.
"""

# ╔═╡ 1d49508e-9a41-4835-9499-ba0f08082d87
md"""
For this case, we choose: 1 Jul 1989, afternoot (hour 15, local time), for a mid latitude of $40^{o}$, and an F10.3 three month average value of 250. All other settings remain the same, except for the calculated parameters to be output.
"""

# ╔═╡ a4a3407d-080c-4ac2-b9c6-0a96e24cc78c
input_high = Dict(
	"year" => "1989",
	"month" => "07",
	"day" => "01",
	"time_flag" => "1",
	"hour" => "15",
	"latitude" => "40.",
	"step" => "25.",
	"f10_7_3" => "250.",
	"vars" => ["05", "11"]
	)

# ╔═╡ d15dc97c-f349-4dfe-be7e-489acca266df
output_high = ReadMSISE(input_high)

# ╔═╡ 11b1301f-92fc-4425-af9e-557036377384
begin
	@df output_low plot(cols(1), cols(5),
		yaxis=:log10,
		size=(625,600),
		ylims=(1e-18,1e-3),
		legend=(0.75,0.9),
		label="Low Activity",
		title="Atmospheric Density Changes with Solar Activity",
		xlabel="Altitude in km",
		ylabel="Mass density in g/cm-3"
	)
	@df output_high plot!(cols(1), cols(2),
		label="High Activity")
end

# ╔═╡ 603f9f54-a9f0-4320-8223-d23617b61224
md"""
The immediate idea I had for application of this in an instructional environment would be to illustrate the effects of high solar activity making the atmosphere 'puffy' by feeding these density profiles into some sort of orbit decay model. That, however, is the proverbial exercise-left-to-reader.
"""

# ╔═╡ Cell order:
# ╟─9ba46bb7-66de-4756-ba4e-613cb1e6f751
# ╟─db0d3958-a55a-4d49-b27e-c586007c8676
# ╠═47b82739-35ba-496c-b5e0-2e73e59cf1c4
# ╠═1b7a3fe9-8716-427f-b796-ca94f7e48a09
# ╠═d8961618-08aa-4424-9cba-960c5023d1fb
# ╠═28dc8709-edc3-4754-a2dd-fc98d0890469
# ╠═dabf432d-434a-4526-a61e-e4d2d61d77d3
# ╟─0d766211-5fd1-4810-b799-5c0f28e088e8
# ╟─177978c2-935d-432e-b5e3-c623d5f6b492
# ╟─234718cf-381a-4d06-8403-67d993922b50
# ╟─d56a1b81-1d17-463b-b5d1-751cf531da76
# ╠═50fc3dd4-fc7b-4283-804f-053c29edbe72
# ╠═a9c5f4df-6c71-449b-8ae9-8195c5b62357
# ╠═150bacf2-b05e-42cf-957a-46df242d4d3a
# ╟─189a1664-bc02-40f0-8347-c5725426faa4
# ╟─b4e055c8-cedc-4aef-b170-f51f5afe17d3
# ╟─1d49508e-9a41-4835-9499-ba0f08082d87
# ╠═a4a3407d-080c-4ac2-b9c6-0a96e24cc78c
# ╠═d15dc97c-f349-4dfe-be7e-489acca266df
# ╠═11b1301f-92fc-4425-af9e-557036377384
# ╟─603f9f54-a9f0-4320-8223-d23617b61224


text:
@testset "Genus representatives" begin
  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x-1;
  K, a = number_field(f)
  D = matrix(K, 3, 3, [-717210130, 0, 0, 0, -55, 0, 0, 0, -1173298395869600]);
  gens = [[1, 0, 0], [1, 0, 0], [4//5, 2//5, 0], [4//5, 2//5, 0], [1493459//5517001, 0, 1//15778622860], [1493459//5517001, 0, 1//15778622860]]
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  @test length(genus_representatives(L)) == 6

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x-1;
  K, a = number_field(f)
  D = matrix(K, 3, 3, [-34, 0, 0, 0, -17, 0, 0, 0, -17192032]);
  gens = [[26, 0, 0], [234, 0, 0], [177//17, 2//17, 0], [177//17, 2//17, 0], [27//2, 0, 1//1768], [27//2, 0, 1//1768]]
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  @test length(genus_representatives(L)) == 93

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x-1;
  K, a = number_field(f)
  D = matrix(K, 3, 3, [-98, 0, 0, 0, -2, 0, 0, 0, -5829824]);
  gens = [[8//7, 0, 0], [8//7, 0, 0], [5//14, 13//2, 0], [45//14, 117//2, 0], [29//28, 3//4, 1//1456], [29//28, 3//4, 1//1456]]
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  @test length(genus_representatives(L)) == 114

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x^2-2;
  K, a = number_field(f)
  D = matrix(K, 3, 3, [2, 0, 0, 0, -1, 0, 0, 0, -7436]);
  gens = [[13, 0, 0], [156*a+143, 0, 0], [3//2*a+5, 1, 0], [3//2*a+5, 1, 0], [21//2*a, 0, 1//26], [21//2*a, 0, 1//26]]
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  @test length(genus_representatives(L)) == 1

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x - 1
  K, a = number_field(f)
  M = matrix(QQ, 2, 2, [47, 80, 80, 560])
  w = Int[]
  for i in 1:10
    M[2, 2] = 560 + i
    V = quadratic_space(K, M)
    L = lattice(V, identity_matrix(K, 2))
    g = length(Hecke._genus_representatives_binary_quadratic_definite(L))
    push!(w, g)
  end
  @test w == Int[37, 16, 19, 11, 6, 3, 52, 7, 13, 21]

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x^2 - 2;
  K, a = number_field(f)
  D = matrix(K, 2, 2, [15, 2, 2, 32]);
  gens = [[1, 0], [1, 0], [0, 1], [0, 1]]
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  @test length(genus_representatives(L)) == 6

  # Local isometry test

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x - 1;
  K, a = number_field(f)
  D = matrix(K, 7, 7, [8, -4, 3, 4, 0, 1, 1, -4, 8, 1, 0, 4, 1, 1, 3, 1, 8, 4, 0, 1, 1, 4, 0, 4, 8, 3, 0, 4, 0, 4, 0, 3, 8, 4, 0, 1, 1, 1, 0, 4, 8, -4, 1, 1, 1, 4, 0, -4, 8]);
  gens = [[1, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 0, 1], [0, 0, 0, 0, 0, 0, 1]]
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)

  D = matrix(K, 7, 7, [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 14, 0, 0, 0, 0, 0, 0, 0, 56]);
  gens = [[1, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 0, 1], [0, 0, 0, 0, 0, 0, 1]]
  LL = quadratic_lattice(K, generators = gens, gram_ambient_space = D)

  p = prime_decomposition(maximal_order(K), 2)[1][1]
  @test islocally_isometric(L, LL, p)

  # Rank 2 case
  # This is the Zlattice with basis [1 2; 3 4]

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x - 1;
  K, a = number_field(f)
  D = matrix(K, 2, 2, [1, 0, 0, 1]);
  gens = [[1, 2], [1, 2], [3, 4], [3, 4]]
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  @test length(genus_representatives(L)) == 1

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x - 1
  K, a = number_field(f)
  D = matrix(K, 3, 3, [-18, -6, -9, -6, -3, -3, -9, -3, -6])
  gens = [[1, 0, 0], [1, 0, 0], [0, 1, 0], [0, 1, 0], [0, 0, 1], [0, 0, 1]]
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  @test length(genus_representatives(L)) == 1

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x - 1;
  K, a = number_field(f)
  D = matrix(K, 2, 2, [2, 0, 0, 3]);
  gens = [[1, 0], [1, 0], [0, 1], [0, 1]]
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  @test length(genus_representatives(L)) == 1

  B = QQ[2 0; 3//2 1//2]
  G = QQ[1 0; 0 23]
  V = quadratic_space(QQ, G)
  L = lattice(V, B)
  @test length(genus_representatives(L)) == 2
  @test length(genus_representatives(Zlattice(gram = gram_matrix(L)))) == 2

  L = Zlattice(ZZ[4 3; 3 8])
  @test length(genus_representatives(L)) == 4
end


text:
module _TestIndirectImportsUpstream

using IndirectImports

"""
    fun

Documentation at the indirect function declaration.
"""
@indirect function fun end

# This does not work:
#
# """
#     fun(x::Complex)
#
# Documentation at the indirect function definition.
# """
@indirect function fun(x::Complex)
    return x + 1 + 1im
end

@indirect function op end

"""
    op

"Off-site" documentation for an indirect function.
"""
op

function reduceop(config, acc, xs)
    for x in xs
        acc = op(config, acc, x)
    end
    return acc
end

end # module


text:
using Weave

function weaveit(name)
    @info "Building $name"
    weave(name, out_path=joinpath("..","docs","src","$(splitext(name)[1]).md"), doctype="github")
end

let start_dir = pwd()
    cd(@__DIR__)
    outpath = normpath(joinpath(@__DIR__, "..", "docs", "src"))
    @show outpath
    if !isdirpath(outpath)
        mkpath(outpath)
    end

    weaveit.(("errorbars.jmd",  "k412refs.jmd", "K412fit.jmd",  "XRFspectra.jmd",  "K412quick.jmd",  "continuummodel.jmd" ))

    cd(start_dir)
end

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
@testset "740.delete-and-earn.jl" begin
    @test delete_and_earn([3, 4, 2]) == 6
    @test delete_and_earn([2, 2, 3, 3, 3, 4]) == 9
end


text:
#definition of particles
struct particle
    pos::Array{Float64,1} #postitions
    pgrid::Array{Float64,1} #positions in grid
    vel::Array{Float64,1} #velocities
    mass::Float64 #mass
    indbox::Int64 #the box where the particle is
    tp::Int64
    function particle(dim::Array{Int64,1}, m::Float64, tp::Int64)
        this = new()
        this.pos = rand(2) .* dim
        this.pgrid = this.pos[:]
        this.vel = [rand()*rand([-1,1]),rand()*rand([-1,1])]
        this.mass = m
        this.indbox = ceil(this.pos[1]) + dim[1] * (ceil(this.pos[2])-1)
        this.tp = tp
        return this
    end
    function particle(x::Int64,dim::Array{Int64,1}, m::Float64, tp::Int64)
        this = new()
        this.pos = [rand()+(x-1),rand()* dim[2]]
        this.pgrid = this.pos[:]
        this.vel = [rand()*rand([-1,1]),rand()*rand([-1,1])]
        this.mass = m
        this.indbox = ceil(this.pos[1]) + dim[1] * (ceil(this.pos[2])-1)
        this.tp = tp
        return this
    end
end
#definition of a box
struct box
    vel::Array{Float64,1} # the velocity of the box.
    np::Array{Int64,1}#  not really sure if is neccesary to have the number of particles of each box.
    ind::Int64
    function box(m::Int64, x::Int64)
        this = new()
        this.vel = zeros(2)
        this.np = zeros(Int64,m)
        this.ind = x
        return this
    end
end


text:
function FieldNodeDetector(arg0::FieldOrbit, arg1::Frame)
    return FieldNodeDetector((FieldOrbit, Frame), arg0, arg1)
end

function FieldNodeDetector(arg0::RealFieldElement, arg1::FieldOrbit, arg2::Frame)
    return FieldNodeDetector((RealFieldElement, FieldOrbit, Frame), arg0, arg1, arg2)
end

function g(obj::FieldNodeDetector, arg0::FieldSpacecraftState)
    return jcall(obj, "g", RealFieldElement, (FieldSpacecraftState,), arg0)
end

function get_frame(obj::FieldNodeDetector)
    return jcall(obj, "getFrame", Frame, ())
end