text:
@defcomp co2cycle begin
    b11     = Variable()                #Carbon cycle transition matrix atmosphere to atmosphere
	b21		= Variable()				#Carbon cycle transition matrix biosphere/shallow oceans to atmosphere
    b22     = Variable()                #Carbon cycle transition matrix shallow ocean to shallow oceans
    b32     = Variable()                #Carbon cycle transition matrix deep ocean to shallow ocean
    b33     = Variable()                #Carbon cycle transition matrix deep ocean to deep oceans
    MAT     = Variable(index=[time])    #Carbon concentration increase in atmosphere (GtC from 1750)
    ML      = Variable(index=[time])    #Carbon concentration increase in lower oceans (GtC from 1750)
    MU      = Variable(index=[time])    #Carbon concentration increase in shallow oceans (GtC from 1750)

    b12     = Parameter()               #Carbon cycle transition matrix atmosphere to shallow ocean
    b23     = Parameter()               #Carbon cycle transition matrix shallow to deep ocean
	mateq	= Parameter()				#Equilibrium concentration atmosphere (GtC)
	mueq 	= Parameter()				#Equilibrium concentration in upper strata (GtC)
	mleq	= Parameter()				#Equilibrium concentration in lower strata (GtC)
    E       = Parameter(index=[time])   #Total CO2 emissions (GtCO2 per year)
    mat0    = Parameter()               #Initial Concentration in atmosphere 2010 (GtC)
    ml0     = Parameter()               #Initial Concentration in lower strata 2010 (GtC)
    mu0     = Parameter()               #Initial Concentration in upper strata 2010 (GtC)

    function run_timestep(p, v, d, t)
		#Define function for b11
		v.b11 = 1 - p.b12

		#Define function for b21
		v.b21 = p.b12 * p.mateq/p.mueq

		#Define function for b22
		v.b22 = 1 - v.b21 - p.b23

		#Define function for b32
		v.b32 = p.b23 * p.mueq/p.mleq

		#Define function for b33
		v.b33 = 1 - v.b32

        #Define function for MAT
        if is_first(t)
            v.MAT[t] = p.mat0
        else
            v.MAT[t] = v.MAT[t-1] * v.b11 + v.MU[t-1] * v.b21 + (p.E[t-1]*(5/3.666))
        end

        #Define function for MU
        if is_first(t)
            v.MU[t] = p.mu0
        else
            v.MU[t] = v.MAT[t-1] * p.b12 + v.MU[t-1] * v.b22 + v.ML[t-1] * v.b32
        end

        #Define function for ML
        if is_first(t)
            v.ML[t] = p.ml0
        else
            v.ML[t] = v.ML[t-1] * v.b33 + v.MU[t-1] * p.b23
        end
    end
end


text:
using Polynomials

function cplxlagrange(tab::Array{Complex{T}}, k::CartesianIndex) where{T}
    p = Polynomial([one(T) + zero(T) * im, 0])
    for ind in CartesianIndices(tab)
        if tab[k] != tab[ind]
            p *= Polynomial([-tab[ind], one(Complex{T})]) / (tab[k] - tab[ind])
        end
    end
    return p
end

function getpoly(tab::Array{Complex{T},2}) where{T}
    p = Polynomial([zero(Complex{T}), zero(Complex{T})])
    for ind in CartesianIndices(tab)
        p += (T(ind.I[1]) + im * T(ind.I[2])) * cplxlagrange(tab, ind)
    end
    return p
end
function getpolyvec(tab::Array{Complex{T},2}) where{T}
    p = Polynomial([zero(Complex{T}), zero(Complex{T})])
    order = size(tab,1)-1
    origin = div(order,2)
    vec = [tab[:,origin] ; tab[origin,:] ]
    for ind in CartesianIndices(vec)
        v = if ind.I[1] <= order+1
            T(ind.I[1]) +im*T(origin)
        else
            T(origin) +im*T( ind.I[1])
        end
        p += v * cplxlagrange(vec, ind)
    end
    return p
end

ext: jl


text:
module PaperUtils

import Plots
import PlotUtils
import ColorTypes
import PyCall
  PyCall.@pyimport gc as pygc

export Normal, Full
export autofig, placeholder, annotate_color!, seqcolors,
       remove_ticklabels!

## Page layout settings
# To describe the extent of a figure we define these terms
@enum PlotSize Normal Full
const plotsizes = (Normal, Full)
"Indicates that a figure spans one column"
Normal
"Indicates that a figure spans the whole page"
Full

# DPI settings
const dpi_default = 100
const dpi_print = 300
const dpi_scale = 170

"Default plot aspect ratios"
const aspectratio = Dict(
  Normal => 1.5,
  Full => 1.6 )

"Widths of the margin, body text, and full page in inches"
const width_in = Dict(
  Normal => 3.3,
  Full => 6.9 )

"Sizes of figures in pixels"
const figsize = let
  pixelsizes = map(plotsizes) do s
    base_aspect = [1, 1/aspectratio[s]]
    dpi_scaled = base_aspect * dpi_print
    width_scaled = dpi_scaled * width_in[s]
    whole_numbers = round(Int, width_scaled)
  end

  Dict(zip(plotsizes, pixelsizes))
end

## Plot display options
# Font settings
const fontstyle = "serif"
const fontsize = Dict(
  Normal => 10,
  Full => 10 )
const font = Dict((s, Plots.font(fontstyle, fontsize[s])) for s in plotsizes)
const main_font_opts = Dict(
  :tickfont => font[Normal],
  :legendfont => font[Normal],
  :guidefont => font[Normal],
  :titlefont => font[Normal] )

# Plot appearance
const plotopts = Dict(
  :linewidth => 1.5,
  :grid => false,
  :dpi => dpi_scale )
Plots.pyplot(;plotopts...)

# Folder to save generated figures to
const figdir = "autofigs"

""" Make a plot to go in the thesis and save it

    plotfunc: The plotting function; should return a Plot object
    name: The file name to save the figure as
    s::PlotSize: One of `Normal` or `Full`. Indicates what size
                 the figure should be saved as

    vscale [optional, default=1]: Scale the plot to create more or less
                                  vertical space
    png [optional, default=1]: Save the plot as a png instead of a pdf """
function autofig(plotfunc, name, s::PlotSize; vscale=1, savepng=false)
  # workaround for python not closing files properly
  pygc.collect()

  # get the figure size and font size
  width = figsize[s][1]
  height = figsize[s][2] * vscale
  Plots.default(;main_font_opts...)
  if s == Full
    name *= "_big_fig"
  end

  # plot the figure and save it to file
  Plots.with(size=(width,height)) do
    p = plotfunc()
    figpath = joinpath(figdir, name)
    if savepng
      Plots.png(figpath)
    else
      Plots.svg(figpath)
      # bit of a workaround because the font doesn't embed
      run(`rsvg-convert -f pdf -o $figpath.pdf $figpath.svg`)
      rm("$figpath.svg")
    end
    p
  end
end

"Make a default plot"
function placeholder()
  Plots.plot([sin, cos], linspace(0, 2π), labels=["sin(x)" "cos(x)"],
    xlabel="This is the x-axis",
    ylabel="This is the y-axis",
    title="Placeholder plot")
end

typealias PlotType Union{Plots.Plot, Plots.Subplot}

"Annotate a plot with a colored label"
function annotate_color!(x, y, text, color;
                         position=:left, plotsize=Normal, rotation=0)
  basefont = font[plotsize]
  rot = deg2rad(rotation)
  formatted_font = Plots.font(basefont, color, position, rot)
  formatted_text = Plots.text(text, formatted_font)
  Plots.annotate!(x, y, formatted_text)
end
function annotate_color!(p::PlotType, x, y, text, color;
                         position=:left, plotsize=Normal, rotation=0)
  basefont = font[plotsize]
  rot = deg2rad(rotation)
  formatted_font = Plots.font(basefont, color, position, rot)
  formatted_text = Plots.text(text, formatted_font)
  Plots.annotate!(p, x, y, formatted_text)
end

"Remove the tick labels from a plot but leave the ticks and grid"
function remove_ticklabels!(p::PlotType; x=true, y=true)
  x && Plots.plot!(p; xformatter=_->"")
  y && Plots.plot!(p; yformatter=_->"")
  p
end

"Provide a palette of colors drawn sequentially from a color gradient"
function seqcolors(name, N, start=0, stop=1)
  grad = PlotUtils.cgrad(name)
  colors = ColorTypes.RGBA{Float64}[grad[n] for n in linspace(start, stop, N)]
end

end # module ThesisUtils


text:
"""
    UnstructuredQuadSortedBoundaryTypes{N, BCs<:NTuple{N, Any}}

General container to sort the boundary conditions by type for the unstructured quadrilateral solver.
It stores a set of global indices for each boundary condition type to expedite computation
during the call to `calc_boundary_flux!`. The original dictionary form of the boundary conditions
set by the user in the elixir file is also stored for printing.

!!! warning "Experimental code"
    This boundary condition container is experimental and can change any time.
"""
struct UnstructuredQuadSortedBoundaryTypes{N, BCs<:NTuple{N, Any}}
  boundary_condition_types::BCs # specific boundary condition type(s), e.g. BoundaryConditionWall
  boundary_indices::NTuple{N, Vector{Int}} # integer vectors containing global boundary indices
  boundary_dictionary::Dict{Symbol, Any} # boundary conditions as set by the user in the elixir file
end

# constructor that "eats" the original boundary condition dictionary and sorts the information
# from the `UnstructuredBoundaryContainer2D` in cache.boundaries according to the boundary types
# and stores the associated global boundary indexing in NTuple
function UnstructuredQuadSortedBoundaryTypes(boundary_conditions::Dict, cache)

  # extract the unique boundary function routines from the dictionary
  boundary_condition_types = Tuple(unique(collect(values(boundary_conditions))))
  n_boundary_types = length(boundary_condition_types)

  # pull and sort the indexing for each boundary type
  _boundary_indices = Vector{Any}(nothing, n_boundary_types)
  for j in 1:n_boundary_types
    indices_for_current_type = Int[]
    for (test_name, test_condition) in boundary_conditions
      temp_indices = findall(x->x===test_name, cache.boundaries.name)
      if test_condition === boundary_condition_types[j]
        indices_for_current_type = vcat(indices_for_current_type, temp_indices)
      end
    end
    _boundary_indices[j] = sort!(indices_for_current_type)
  end

  # convert the work array with the boundary indices into a tuple
  boundary_indices = Tuple(_boundary_indices)

  boundary_dictionary = boundary_conditions

  return UnstructuredQuadSortedBoundaryTypes{n_boundary_types, typeof(boundary_condition_types)}(
                                                                         boundary_condition_types,
                                                                         boundary_indices,
                                                                         boundary_dictionary)
end


text:
using Base.Test
using MIPVerify
using MIPVerify: UnrestrictedPerturbationFamily, BlurringPerturbationFamily
using MIPVerify: get_example_network_params, read_datasets, get_image
isdefined(:TestHelpers) || include("../../../TestHelpers.jl")
using TestHelpers: batch_test_adversarial_example

@testset "MNIST.n1" begin
    nn = get_example_network_params("MNIST.n1")
    mnist = read_datasets("mnist")

    pp_blur = BlurringPerturbationFamily((5, 5))
    pp_unrestricted = UnrestrictedPerturbationFamily()

    sample_index = 1
    @testset "sample index = $sample_index" begin
        x0 = get_image(mnist.test.images, sample_index)

        expected_objective_values = Dict(
            (1, pp_unrestricted, 1, 0) => 13.8219,
            (1, pp_unrestricted, Inf, 0) => 0.0964639,
            (3, pp_blur, 1, 0) => 88.0038,
            (3, pp_blur, Inf, 0) => 0.973011,
            (10, pp_unrestricted, 1, 0) => 4.64186,
            (10, pp_unrestricted, Inf, 0) => 0.0460847,
            (2, pp_blur, 1, 0) => NaN,
            (2, pp_blur, Inf, 0) => NaN,
            (5, pp_blur, 1, 0) => NaN,
            (5, pp_blur, Inf, 0) => NaN,
            (6, pp_blur, 1, 0) => NaN,
            (6, pp_blur, Inf, 0) => NaN
        )

        batch_test_adversarial_example(nn, x0, expected_objective_values)
    end

    sample_index = 3
    @testset "sample index = $sample_index" begin
        x0 = get_image(mnist.test.images, sample_index)

        expected_objective_values = Dict(
            (8, pp_unrestricted, 1, 0) => 0.582587,
            (8, pp_unrestricted, Inf, 0) => 0.00363594,
            (8, pp_blur, 1, 0) => 9.77784,
            (8, pp_blur, Inf, 0) => 0.277525
        )

        batch_test_adversarial_example(nn, x0, expected_objective_values)
    end

end


text:
function part1(input)
    return sum(count_score1.(split(read(input, String), "\n\n")))
end

function part2(input)
    return sum(count_score2.(split(read(input, String), "\n\n")))
end

function count_score1(s)
    s = replace(s, "\n" => "")
    return length(unique(split(s, "")))
end

function count_score2(s)
    answers = split(s, "\n", keepempty = false)
    counts = Dict{String, Int}()
    for answer in answers
        for question in unique(split(answer, ""))
            counts[question] = get(counts, question, 0) + 1
        end
    end
    return count(==(length(answers)), values(counts))
end


text:
# Simple buffered input that allows peeking an arbitrary number of characters
# ahead by maintaining a typically quite small buffer of a few characters.

mutable struct BufferedInput
    input::IO
    buffer::Vector{Char}
    offset::UInt64
    avail::UInt64

    function BufferedInput(input::IO)
        return new(input, Vector{Char}(undef, 0), 0, 0)
    end
end

# Read and buffer n more characters
function __fill(bi::BufferedInput, bi_input::IO, n::Integer)
    for i in 1:n
        c = eof(bi_input) ? '\0' : read(bi_input, Char)
        if bi.offset + bi.avail + 1 <= length(bi.buffer)
            bi.buffer[bi.offset + bi.avail + 1] = c
        else
            push!(bi.buffer, c)
        end
        bi.avail += 1
    end
end

_fill(bi::BufferedInput, n::Integer) = __fill(bi, bi.input, n)

# Peek the character in the i-th position relative to the current position.
# (0-based)
function peek(bi::BufferedInput, i::Integer=0)
    if bi.avail < i + 1
        _fill(bi, i + 1 - bi.avail)
    end
    return bi.buffer[bi.offset + i + 1]
end

# Return the string formed from the first n characters from the current position
# of the stream.
function prefix(bi::BufferedInput, n::Integer=1)
    if bi.avail < n + 1
        _fill(bi, n + 1 - bi.avail)
    end
    return string(bi.buffer[(bi.offset + 1):(bi.offset + n)]...)
end

# NOPE: This is wrong. What if n > bi.avail

# Advance the stream by n characters.
function forward!(bi::BufferedInput, n::Integer=1)
    if n < bi.avail
        bi.offset += n
        bi.avail -= n
    else
        n -= bi.avail
        bi.offset = 0
        bi.avail = 0
        while n > 0
            read(bi.input, Char)
            n -= 1
        end
    end
end

# Ugly hack to allow peeking of `StringDecoder`s
function peek(io::StringDecoder, ::Type{UInt8})
    c = read(io, UInt8)
    io.skip -= 1
    c
end

# The same but for Julia 1.3
peek(io::StringDecoder) = peek(io, UInt8)


text:
export set_conv

"""
    SetConv{T<:AV{<:Real}}

A set convolution layer.

# Fields
- `log_scales::T`: Natural logarithm of the length scales of every input channel.
- `density::Bool`: Include the density channel.
"""
struct SetConv{T<:AV{<:Real}}
    log_scales::T
    density::Bool
end

@Flux.functor SetConv

"""
    SetConvPD{T<:AV{<:Real}}

A set convolution layer for positive-definite-matrix outputs.

# Fields
- `log_scales::T`: Natural logarithm of the length scales of every input channel.
- `density::Bool`: Include the density channel.
"""
struct SetConvPD{T<:AV{<:Real}}
    log_scales::T
    density::Bool
end

@Flux.functor SetConvPD

"""
    set_conv(num_channels::Integer, scale::Float32; density::Bool=false)

Construct a set convolution layer.

# Arguments
- `num_channels::Integer`: Number of inputs channels, excluding the density channel.
- `scale::Float32`: Initialisation of the length scales.

# Keywords
- `density::Bool=false`: Include the density channel.
- `pd::Bool=false`: Generate positive-definite matrices as outputs.

# Returns
- `SetConv`: Corresponding set convolution layer.
"""
function set_conv(
    num_channels::Integer,
    scale::Float32;
    density::Bool=false,
    pd::Bool=false
)
    density && (num_channels += 1)
    scales = scale .* ones(Float32, num_channels)
    return (pd ? SetConvPD : SetConv)(log.(scales), density)
end

_get_scales(layer) = reshape(exp.(layer.log_scales), 1, 1, length(layer.log_scales), 1)

function _compute_weights(x, y, scales)
    dists² = compute_dists²(x, y)
    dists² = insert_dim(dists², pos=3)  # Add channel dimension.
    dists² = dists² ./ scales.^2
    return rbf(dists²)
end

function _prepend_density_channel(z)
    n, _, batch_size = size(z)
    density = ones_gpu(data_eltype(z), n, 1, batch_size)
    return cat(density, z, dims=2)
end

function _prepend_identity_channel(z)
    n, _, _, batch_size = size(z)
    identity = gpu(repeat(Matrix{Float32}(I, n, n), 1, 1, 1, batch_size))
    return cat(z, identity, dims=3)
end

function _normalise_by_first_channel(z)
    channels_dim = ndims(z) - 1  # Channels dimension is second to last.
    normaliser = slice_at(z, channels_dim, 1:1)
    others = slice_at(z, channels_dim, 2:size(z, channels_dim))
    return cat(normaliser, others ./ (normaliser .+ 1f-8), dims=channels_dim)
end

function code(layer::SetConv, xz::AA, z::AA, x::AA; kws...)
    weights = _compute_weights(x, xz, _get_scales(layer))
    layer.density && (z = _prepend_density_channel(z))
    if ndims(z) == 4
        z = permutedims(z, (1, 4, 2, 3))
        z = batched_mul(weights, z)
        z = permutedims(z, (1, 3, 4, 2))
    elseif ndims(z) == 3
        z = with_dummy(c -> batched_mul(weights, c), z)
    else
        error("Cannot deal with inputs of rank $(ndims(z)).")
    end
    layer.density && (z = _normalise_by_first_channel(z))
    return x, z
end

function code(layer::SetConv, xz::Nothing, z::Nothing, x::AA; kws...)
    return x, zeros_gpu(
        data_eltype(x),
        size(x, 1),               # Size of encoding
        length(layer.log_scales), # Number of z, including the density channel
        size(x, 3)                # Batch size
    )
end

function code(layer::SetConvPD, xz::AA, z::AA, x::AA; kws...)
    weights = _compute_weights(x, xz, _get_scales(layer))
    # TODO: Disengtangle the below.
    if layer.density
        z = insert_dim(_prepend_density_channel(z), pos=1)
        z = batched_mul(weights .* z, batched_transpose(weights))
        z = _normalise_by_first_channel(z)
        z = _prepend_identity_channel(z)
    else
        Ls = batched_mul(weights, z)
        z = batched_mul(Ls, batched_transpose(Ls))
    end
    return x, z
end

function code(layer::SetConvPD, xz::Nothing, z::Nothing, x::AA; kws...)
    z = zeros_gpu(
        data_eltype(x),
        size(x, 1),               # Size of encoding
        size(x, 1),               # Again size of encoding: encoding is square
        length(layer.log_scales), # Number of z, including the density channel
        size(x, 3)                # Batch size
    )
    # Also prepend density channel.
    layer.density && (z = _prepend_identity_channel(z))
    return x, z
end


text:
# Script to Read simul_teste.pdb and put its format in a suitable shape (for gmd.sh program)
# vinicius piccoli
########################################################

# Changing "0.00   " by "0.00   UBQ"
name = "simul_teste.pdb"
open("simul_alt.pdb","w") do action
  open(name, "r") do f
    for line in eachline(f)
      a =replace(line,"0.00            " => "0.00       UBQ",count = 1)
      write(action,"$a\n")
    end
  end
  close(action)
end

# Adjusting molecules names in the pdb file
for line in eachline(open("simul_alt.pdb","r+"))
  if  occursin("SOL",line)
    println(replace(line,"UBQ" => "SOL",count = 1))
  elseif occursin("BMI",line)
    println(replace(line,"UBQ" => "BMI",count = 1))
  elseif occursin("NC ",line)
    println(replace(line,"UBQ" => "NC ",count = 1))
  elseif occursin("Cl ",line)
    println(replace(line,"UBQ" => "Cl ",count = 1))
  elseif occursin("NO3",line)
    println(replace(line,"UBQ" => "NO3",count = 1))
  elseif occursin("Br ",line)
    println(replace(line,"UBQ" => "Br ",count = 1))
  elseif occursin("BF4",line)
    println(replace(line,"UBQ" => "BF4",count = 1))
  elseif occursin("MS ",line)
    println(replace(line,"UBQ" => "MS ",count = 1))
  elseif occursin("TCM",line)
    println(replace(line,"UBQ" => "TCM",count = 1))
  elseif occursin("EMI",line)
    println(replace(line,"UBQ" => "EMI",count = 1))
  elseif occursin("NA ",line)
    println(replace(line,"NA " => "OMI",count = 1))
  else
    println(line)
  end
end


text:
# Absolute power set

# Note: It might be more natural to define a Function for abspower and add
# constraints with the Interval set. But this would involve more effort to
# support it through JuMP, so we define a set instead.

"""
Set of (x,z) that satisfy the absolute power constraint:

    lhs ≤ sign(x + a) * abs(x + a)^n + c*z ≤ rhs

The exponent `n` should be greater than 1.
"""
struct AbsolutePowerSet{T <: Real} <: MOI.AbstractVectorSet
    n::T
    a::T
    c::T
    lhs::T
    rhs::T
end

"""
Shortcut constructor for the special case:

    z == sign(x) * abs(x)^n
"""
AbsolutePowerSet(n::Float64) = AbsolutePowerSet(n, 0.0, -1.0, 0.0, 0.0)

MOI.dimension(set::AbsolutePowerSet{<:Real}) = 2

# Absolute power constraint

const ABSPOWER = AbsolutePowerSet{Float64}

MOI.supports_constraint(::Optimizer, ::Type{VECTOR}, ::Type{ABSPOWER}) = true

function MOI.add_constraint(o::Optimizer, func::VECTOR, set::ABSPOWER)
    allow_modification(o)

    x, z = [VarRef(vi.value) for vi in func.variables]
    cr = add_abspower_constraint(o.mscip, x, set.a, set.n, z, set.c, set.lhs, set.rhs)
    ci = CI{VECTOR, ABSPOWER}(cr.val)
    register!(o, ci)
    register!(o, cons(o, ci), cr)
    return ci
end

function MOI.delete(o::Optimizer, ci::CI{VECTOR, ABSPOWER})
    allow_modification(o)
    delete!(o.constypes[VECTOR, ABSPOWER], ConsRef(ci.value))
    delete!(o.reference, cons(o, ci))
    delete(o.mscip, ConsRef(ci.value))
    return nothing
end

function MOI.get(o::Optimizer, ::MOI.ConstraintFunction, ci::CI{VECTOR, ABSPOWER})
    c = cons(o, ci)::Ptr{SCIP_CONS}
    lvar = SCIPgetLinearVarAbspower(o, c)
    nlvar = SCIPgetNonlinearVarAbspower(o, c)
    return VECTOR([VI(ref(o, nlvar).val), VI(ref(o, lvar).val)])
end

function MOI.get(o::Optimizer, ::MOI.ConstraintSet, ci::CI{VECTOR, ABSPOWER})
    c = cons(o, ci)::Ptr{SCIP_CONS}

    n = SCIPgetExponentAbspower(o, c)
    a = SCIPgetOffsetAbspower(o, c)
    coef = SCIPgetCoefLinearAbspower(o, c)
    lhs = SCIPgetLhsAbspower(o, c)
    rhs = SCIPgetRhsAbspower(o, c)

    return AbsolutePowerSet(n, a, coef, lhs, rhs)
end


text:
"""
Abstract atomistic or single-step solution.
"""
abstract type AtomicSolution{dType <: Number, tType <: Real} end

"Create AtomicSolution for ODE."
function AtomicSolution(equation::AbstractEquationODE{DT,TT}) where {DT,TT}
    AtomicSolutionODE{DT,TT}(ndims(equation))
end

"Create AtomicSolution for partitioned ODE."
function AtomicSolution(equation::AbstractEquationPODE{DT,TT}) where {DT,TT}
    AtomicSolutionPODE{DT,TT}(ndims(equation))
end

"Create AtomicSolution for DAE."
function AtomicSolution(equation::AbstractEquationDAE{DT,TT}) where {DT,TT}
    AtomicSolutionDAE{DT,TT}(ndims(equation))
end

"Create AtomicSolution for partitioned DAE."
function AtomicSolution(equation::AbstractEquationPDAE{DT,TT}) where {DT,TT}
    AtomicSolutionPDAE{DT,TT}(ndims(equation))
end

"Print error for AtomicSolutions of equations not implemented, yet."
function AtomicSolution(equation::Equation)
    error("No AtomicSolution found for equation ", equation)
end

"Copy solution from atomistic solution to solution object."
function copy_solution!(sol::Solution, asol::AtomicSolution, n, m)
    copy_solution!(sol, get_solution(asol)..., n, m)
end

function CommonFunctions.cut_periodic_solution!(asol::AtomicSolution{DT}, periodicity::Vector{DT}) where {DT}
    @assert length(asol.q) == length(periodicity)
    for k in eachindex(asol.q, asol.q̅, periodicity)
        if periodicity[k] ≠ 0
            while asol.q[k] < 0
                asol.q[k], asol.q̃[k] = compensated_summation(+periodicity[k], asol.q[k], asol.q̃[k])
                asol.q̅[k] += periodicity[k]
            end
            while asol.q[k] ≥ periodicity[k]
                asol.q[k], asol.q̃[k] = compensated_summation(-periodicity[k], asol.q[k], asol.q̃[k])
                asol.q̅[k] -= periodicity[k]
            end
        end
    end
end


text:
module AE

include("MNIST/MNIST.jl")
using Flux
using CUDA
using JLD2

function get_test_data_and_label()
    test_x, test_y  = MNIST.testdata(Float32)

    return test_x,  test_y
end

function get_data(batch_size)
    train_x, _ = MNIST.traindata(Float32)
    test_x, _ = MNIST.testdata(Float32)

    return (
        Flux.DataLoader((train_x, train_x), batchsize=batch_size, shuffle=true),
        Flux.DataLoader((test_x, test_x), batchsize=batch_size, shuffle=false)
    )
end

function update_model!(model_file_path, model)
    model = cpu(model)
    jldsave(model_file_path; model)
    @warn "model updated!"
end

function train(; η₀=1e-4, epoch=500, batch_size=100)
    if has_cuda()
        @info "CUDA is on"
        device = gpu
        CUDA.allowscalar(false)
    else
        device = cpu
    end

    m = Chain(
        flatten,
        Dense(28*28, 512, relu),
        Dense(512, 32, relu),

        Dense(32, 512, relu),
        Dense(512, 28*28, relu),
        x -> reshape(x, 28, 28, :)
    ) |> device

    loss(𝐱, 𝐲) = sum(abs2, 𝐲 .- m(𝐱)) / size(𝐱)[end]

    opt = Flux.ADAM(η₀)

    @info "gen data... "
    @time loader_train, loader_test = get_data(batch_size)

    losses = Float32[]
    function validate()
        validation_loss = sum(loss(device(𝐱), device(𝐲)) for (𝐱, 𝐲) in loader_test)/length(loader_test)
        @info "loss: $validation_loss"

        push!(losses, validation_loss)
        (losses[end] == minimum(losses)) && update_model!(joinpath(@__DIR__, "../model/model_h.jld2"), m)
    end

    data = [(𝐱, 𝐲) for (𝐱, 𝐲) in loader_train] |> device
    Flux.@epochs epoch @time begin
        Flux.train!(loss, params(m), data, opt)
        validate()
    end
end

function train_dim_reduction(; η₀=1e-4, epoch=500, batch_size=100)
    if has_cuda()
        @info "CUDA is on"
        device = gpu
        CUDA.allowscalar(false)
    else
        device = cpu
    end

    m = get_model("model_h")
    m = Chain(
        m[1:3],
        Dense(32, 3),
        Dense(3, 32),
        m[4:6]
    ) |> device

    loss(𝐱, 𝐲) = sum(abs2, 𝐲 .- m(𝐱)) / size(𝐱)[end]

    opt = Flux.ADAM(η₀)

    @info "gen data... "
    @time loader_train, loader_test = get_data(batch_size)

    losses = Float32[]
    function validate()
        validation_loss = sum(loss(device(𝐱), device(𝐲)) for (𝐱, 𝐲) in loader_test)/length(loader_test)
        @info "loss: $validation_loss"

        push!(losses, validation_loss)
        (losses[end] == minimum(losses)) && update_model!(joinpath(@__DIR__, "../model/model.jld2"), m)
    end

    data = [(𝐱, 𝐲) for (𝐱, 𝐲) in loader_train] |> device
    Flux.@epochs epoch @time begin
        Flux.train!(loss, params(m), data, opt)
        validate()
    end
end

function get_model(name="model")
    f = jldopen(joinpath(@__DIR__, "../model/$name.jld2"))
    model = f["model"]
    close(f)

    return model
end

end


text:
using NMRTools
using Test
using SafeTestsets

@safetestset "NMRData" begin
    include("nmrdata_test.jl")
end


text:
# Question
# Find the sum of all the multiples of 3 or 5 below 1000.

# Time
# O(1)

function ap(n, d, a)
    return n*(2*a + (n-1)*d)/2
end
function main()
    limit = 1000
    n3 =  floor(1000/3)
    n5 = floor((1000-1)/5)
    n15 = floor(1000/15)

    sum = ap(n3, 3, 3) + ap(n5, 5, 5) - ap(n15, 15, 15)
    println(convert(Int, trunc(sum)))
end

@time main()


text:
using OrdinaryDiffEq
using ModelingToolkit
using DataDrivenDiffEq
using LinearAlgebra, DiffEqSensitivity, Optim
using DiffEqFlux, Flux
using Plots
gr()

function triNN!(du,u,p,t,dens,cons)
    # unpack rates and constants
    nᵣ,nₓ,n₃ = u
    k₁,k₋₁,k₂,k₋₂ = cons
    mᵣ,mₗ,mₓ,A = dens
    z = L(u,p)
    # model
    du[1] = dnᵣ = A*k₁*mᵣ*mₗ - k₋₁*nᵣ - z[1] + z[2]
    du[2] = dnₓ = A*k₂*mₓ*mₗ - k₋₂*nₓ - z[3] + z[4]
    du[3] = dn₃ = z[1] + z[3] - z[2] - z[4] + z[5]

end
function triNN!(du,u,p,t,dens,cons)
    # unpack rates and constants
    nᵣ,nₓ,n₃ = u
    k₁,k₋₁,k₂,k₋₂ = cons
    mᵣ,mₗ,mₓ,A = dens
    z = L(u,p)
    # model
    du[1] = dnᵣ = A*k₁*mᵣ*mₗ - k₋₁*nᵣ + z[1] - z[2]
    du[2] = dnₓ = A*k₂*mₓ*mₗ - k₋₂*nₓ + z[3] - z[4]
    du[3] = dn₃ = z[5] - z[6]

end

L = FastChain(FastDense(3, 50, tanh),
                FastDense(50, 100, tanh),
                FastDense(100, 100, tanh),
                FastDense(100, 50, tanh),
                FastDense(50, 6))
#L = FastChain(FastDense(3, 20, tanh),FastDense(20, 20, tanh), FastDense(20, 2))
p = initial_params(L)

# Define the experimental parameter
tspan = (0.0,16.1) # A little bit longer because there are two values to save at final t?
u0 = Float32[0.0,0.0,0.0]

# Xₙ = n_tc1 #load data
Xₙ = [0.02020270731751947
 0.0
 0.06187540371808745
 0.10536051565782635
 0.083381608939051
 0.19845093872383823
 0.3011050927839216
 0.3011050927839216
 0.3566749439387324
 0.38566248081198473
 0.5108256237659907
 0.6931471805599453
 0.8209805520698303
 0.7339691750802005
 0.8209805520698303
 0.7339691750802005]
 t = [0.25
  0.25
  0.5
  0.5
  1.0
  1.0
  2.0
  2.0
  4.0
  4.0
  6.0
  6.0
  8.0
  8.0
 16.0
 16.0]
densities_tc1 = [15.0,38.0,10.0,1.0]
cons_tc = [ 0.0006692541890287495,0.8662361534770547,1.169965568192585e-6,0.4]

f = (du,u,p,t) -> triNN!(du,u, p,t,densities_tc1,cons_tc)
prob_nn = ODEProblem(f,u0, tspan, p)
sol_nn = solve(prob_nn, Tsit5(), saveat = t)

# plot(solution)
plot(sol_nn)
# summ = reduce(vcat,sum(sol_nn,dims=1))
# h = plot!(sol_nn.t,summ,linecolor=:black)

function predict(θ)

    tmp_prob = remake(prob_nn,u0=u0,p=θ)
    tmp_sol =  solve(tmp_prob, VCABM(), saveat = t,
                  abstol=1e-5, reltol=1e-5)
                  # backsolve=true)
    Σ_sol = sum(Array(tmp_sol),dims=1) # Note: this returns a row vector!
end

function loss(θ)
    pred = predict(θ)
    # sum(abs2, (Xₙ .- pred') .* t), pred
    sum(abs2, (Xₙ .- pred')), pred
end

# Test
@time loss(p)

losses = []

callback(θ,l,pred) = begin
    push!(losses, l)
    @show l
    if length(losses)%50==0
        println("Current loss after $(length(losses)) iterations: $(losses[end])")
    end
    p = plot(t, pred')
    scatter!(p, t, Xₙ)
    display(p)
    false
end

# First train with ADAM for better convergence
@time res1 = DiffEqFlux.sciml_train(loss, p, ADAM(0.005), cb=callback, maxiters = 10)

# Train with BFGS
res2 = DiffEqFlux.sciml_train(loss, res1.minimizer, BFGS(initial_stepnorm=0.01),
                                cb=callback, maxiters = 1000)

res2 = DiffEqFlux.sciml_train(loss, res1.minimizer, NelderMead(),
                                cb=callback)

println("Final training loss after $(length(losses)) iterations: $(losses[end])")

# Plot the losses
plot(losses, yaxis = :log, xaxis = :log, xlabel = "Iterations", ylabel = "Loss")

# Plot the data and the approximation
NNsolution = predict(res2.minimizer)
# Trained on noisy data vs real solution
plot(t, NNsolution')
scatter!(t, Xₙ)


text:
using CoolPropExtensions
using Test

@testset "CoolPropExtensions.jl" begin
    atm = 101325
    constant_pressure_specific_heat = CoolPropExtensions.Property("C")
    pressuretemperature = CoolPropExtensions.InputPair("PT_INPUTS")
    water = CoolPropExtensions.WorkingFluid("HEOS","Water")
    CoolPropExtensions.update(water, pressuretemperature, 101325, 300)
    c_p = CoolPropExtensions.keyed_output(water,
        constant_pressure_specific_heat)
    @test c_p≈4180.6357765560715
    mass_density = CoolPropExtensions.Property("Dmass")
    pressure = CoolPropExtensions.Property("P")
    temperature = CoolPropExtensions.Property("T")
    pd = CoolPropExtensions.first_partial_deriv(water,
        mass_density, pressure, temperature)
    @test pd≈4.489626212601196e-7
    let
        work = CoolPropExtensions.isothermalwork(WorkingFluid("HEOS","Water"),
            300.0, 10atm, 100atm, rtol=1e-8)
        @test work ≈ 9128.426693091658
    end
end


text:
using PDFmerger
using Test
import Base.Filesystem

test_files = [joinpath(pkgdir(PDFmerger), "test", "file_1.pdf"),
              joinpath(pkgdir(PDFmerger), "test", "file_2.pdf")]

function make_n_copies(file, n)
    for i in 1:n
        cp(file, replace(file, r"\.pdf$" => "$(i).pdf"), force=true)
    end
end

@testset "merge_pdfs" begin

    # -- no cleanup
    # setup test directory
    test_dir = mktempdir()
    single_files = joinpath.(test_dir, ["file_1.pdf", "file_2.pdf"])
    cp.(test_files, single_files)

    out_file = joinpath(test_dir, "out.pdf")

    merge_pdfs(single_files, out_file)
    @test PDFmerger.n_pages(out_file) == 2

    merge_pdfs(single_files, out_file)
    @test PDFmerger.n_pages(out_file) == 2

    # test merging existing file with itself
    merge_pdfs([single_files..., out_file], out_file)
    @test PDFmerger.n_pages(out_file) == 2 + 2

    @test readdir(test_dir) |> length == 3
    @test all(PDFmerger.n_pages.(single_files) .== 1)

    # -- with cleanup
    # setup test directory
    test_dir = mktempdir()
    single_files = joinpath.(test_dir, ["file_1.pdf", "file_2.pdf"])
    cp.(test_files, single_files)

    out_file = joinpath(test_dir, "out.pdf")

    merge_pdfs(single_files, out_file, cleanup=true)
    @test PDFmerger.n_pages(out_file) == 2
    @test readdir(test_dir) |> length == 1

    # -- large number of files
    # setup test directory
    test_dir = mktempdir()
    single_file = joinpath(test_dir, "file.pdf")
    cp(test_files[1], single_file)

    out_file = joinpath(test_dir, "out.pdf")

    n = 1500
    make_n_copies(single_file, n)
    readdir(test_dir)
    merge_pdfs([joinpath(test_dir, "file$i.pdf") for i in 1:n], out_file, cleanup=true)
    @test PDFmerger.n_pages(out_file) == n
    @test readdir(test_dir) |> length == 2
end

@testset "append_pdf" begin

    # -- no cleanup
    # setup test directory
    test_dir = mktempdir()
    single_files = joinpath.(test_dir, ["file_1.pdf", "file_2.pdf"])
    cp.(test_files, single_files)

    out_file = joinpath(test_dir, "out.pdf")

    @test_throws ErrorException append_pdf!(out_file, single_files[1], create=false)

    for k in 1:10
        append_pdf!(out_file, single_files[1])
        @test PDFmerger.n_pages(out_file) == k
    end

    # test existing file with itself
    append_pdf!(out_file, out_file)
    @test PDFmerger.n_pages(out_file) == 2*10
    @test all(PDFmerger.n_pages.(single_files) .== 1)

    @test readdir(test_dir) |> length == 3

    # -- with cleanup
    # setup test directory
    test_dir = mktempdir()
    single_files = joinpath.(test_dir, ["file_1.pdf", "file_2.pdf"])
    cp.(test_files, single_files)

    out_file = joinpath(test_dir, "out.pdf")

    readdir(test_dir)
    @test_throws ErrorException append_pdf!(out_file, single_files[1], create=false, cleanup=true)

    readdir(test_dir)

    append_pdf!(out_file, single_files[1], cleanup=true)

    @test PDFmerger.n_pages(out_file) == 1
    @test "out.pdf" ∈ readdir(test_dir)
    @test "file_1.pdf" ∉ readdir(test_dir)

    append_pdf!(out_file, single_files[2], cleanup=true)
    @test PDFmerger.n_pages(out_file) == 2
    @test "out.pdf" ∈ readdir(test_dir)
    @test "file_2.pdf" ∉ readdir(test_dir)

    # test existing file with itself
    append_pdf!(out_file, out_file, cleanup=true)
    @test PDFmerger.n_pages(out_file) == 2*2

    @test readdir(test_dir) |> length == 1

end


text:
export Train, TrainStatusPrint

function Train(HN,logs,TrOpts,train_data,val_data,train_labels,val_labels,P,image_weights_train,image_weights_val,active_z_slice::Union{Array{Any,1},Int}=[])

  for j=1:TrOpts.maxiter
    for k=1:TrOpts.batchsize
      rand_ind = randperm(length(train_data))[1]
      # Evaluate objective and gradients
      fval, dc2val = LossTotal(HN,TrOpts,train_data[rand_ind],train_labels[rand_ind],P[rand_ind],image_weights_train[rand_ind],active_z_slice)
    end
    for p in get_params(HN)
        update!(TrOpts.opt, p.data, p.grad)
        #update!(lr_decay_fn, p.data, p.grad)
    end
    clear_grad!(HN)

    if mod(j, TrOpts.eval_every) == 0 || j==1
      if isempty(train_labels[1])==false
        ioupos_train,iouneg_train = IoU(HN,train_data,train_labels)
        ioupos_train   = ioupos_train[ioupos_train.>0.0]
        logs.IoU_train = cat(logs.IoU_train,[mean(ioupos_train) mean(iouneg_train)],dims=1)
      end
      fvalepoch_train = 0.0
      dvalepoch_train = 0.0
      for i=1:length(train_data)
        f,d = LossTotal(HN,TrOpts,train_data[i],train_labels[i],P[i],image_weights_train[i],active_z_slice)
        fvalepoch_train = fvalepoch_train + f
        dvalepoch_train = dvalepoch_train + d
      end
      logs.train     = vcat(logs.train,fvalepoch_train/length(train_data))
      logs.dc2_train = vcat(logs.dc2_train,dvalepoch_train/length(train_data))

      #validation data/labels
      if isempty(val_labels[1])==false
        ioupos_val,iouneg_val = IoU(HN,val_data,val_labels)
        ioupos_val    = ioupos_val[ioupos_val.>0.0]
        logs.IoU_val  = cat(logs.IoU_val,[mean(ioupos_val) mean(iouneg_val)],dims=1)

        fvalepoch_val = 0.0
        dvalepoch_val = 0.0
        for i=1:length(val_data)
          f,d = LossTotal(HN,TrOpts,val_data[i],val_labels[i],P[i],image_weights_val[i],active_z_slice)
          fvalepoch_val = fvalepoch_val + f
          dvalepoch_val = dvalepoch_val + d
        end
        logs.val     = vcat(logs.val,fvalepoch_val/length(val_data))
        logs.dc2_val = vcat(logs.dc2_val,dvalepoch_val/length(val_data))
      end

      TrainStatusPrint(j,logs)
      #print("Iteration: ", j, "; ftrain = ", logs.train[end], "; dtrain = ", logs.dc2_train[end], "; fval = ", logs.val[end], "; dval = ", logs.dc2_val[end], ";  IoUtrain:", logs.IoU_train[end,:] , ";  IoUval:" , logs.IoU_val[end,:], "\n")
      clear_grad!(HN)
    end

  end

  return logs
end

function Train(HN,logs,TrOpts,train_data,val_data,train_labels,val_labels,P,image_weights_train,image_weights_val,active_z_slice,P_mode::String,TD_OP,P_sub,alpha_CQ)

  for j=1:TrOpts.maxiter
    for k=1:TrOpts.batchsize
      rand_ind = randperm(length(train_data))[1]
      # Evaluate objective and gradients
      fval, dc2val = LossTotal(HN,TrOpts,train_data[rand_ind],train_labels[rand_ind],P[rand_ind],image_weights_train[rand_ind],active_z_slice,P_mode,TD_OP,P_sub[rand_ind],alpha_CQ)
    end
    for p in get_params(HN)
        update!(TrOpts.opt, p.data, p.grad)
        #update!(lr_decay_fn, p.data, p.grad)
    end
    clear_grad!(HN)

    if mod(j, TrOpts.eval_every) == 0 || j==1
      if isempty(train_labels[1])==false
        ioupos_train,iouneg_train = IoU(HN,train_data,train_labels)
        ioupos_train   = ioupos_train[ioupos_train.>0.0]
        logs.IoU_train = cat(logs.IoU_train,[mean(ioupos_train) mean(iouneg_train)],dims=1)
      end
      fvalepoch_train = 0.0
      dvalepoch_train = 0.0
      for i=1:length(train_data)
        f,d = LossTotal(HN,TrOpts,train_data[i],train_labels[i],P[i],image_weights_train[i],active_z_slice,P_mode,TD_OP,P_sub[i],alpha_CQ)
        fvalepoch_train = fvalepoch_train + f
        dvalepoch_train = dvalepoch_train + d
      end
      logs.train       = vcat(logs.train,fvalepoch_train/length(train_data))
      logs.dc2_train   = vcat(logs.dc2_train,dvalepoch_train/length(train_data))

      #validation data/labels
      if isempty(val_labels[1])==false
        ioupos_val,iouneg_val = IoU(HN,val_data,val_labels)
        ioupos_val=ioupos_val[ioupos_val.>0.0]
        logs.IoU_val  = cat(logs.IoU_val,[mean(ioupos_val) mean(iouneg_val)],dims=1)

        fvalepoch_val = 0.0
        dvalepoch_val = 0.0
        for i=1:length(val_data)
          PI = x -> x#no projector for validation
          f,d = LossTotal(HN,TrOpts,val_data[i],val_labels[i],PI,image_weights_val[i],active_z_slice)
          fvalepoch_val = fvalepoch_val + f
          dvalepoch_val = dvalepoch_val + d
        end
        logs.val     = vcat(logs.val,fvalepoch_val/length(val_data))
        logs.dc2_val = vcat(logs.dc2_val,dvalepoch_val/length(val_data))
      end

      TrainStatusPrint(j,logs)
      #print("Iteration: ", j, "; ftrain = ", logs.train[end], "; dtrain = ", logs.dc2_train[end], "; fval = ", logs.val[end], "; dval = ", logs.dc2_val[end], ";  IoUtrain:", logs.IoU_train[end,:] , ";  IoUval:" , logs.IoU_val[end,:], "\n")
      clear_grad!(HN)
    end

  end

  return logs
end

function TrainStatusPrint(iteration,logs)
  printstring=string("Iteration: ", iteration)
  if isempty(logs.train)==false
    printstring = printstring * string("; ftrain = ", logs.train[end])
  end
  if isempty(logs.dc2_train)==false
    printstring = printstring * string("; dtrain = ", logs.dc2_train[end])
  end
  if isempty(logs.val)==false
    printstring = printstring * string("; fval = ", logs.val[end])
  end
  if isempty(logs.dc2_val)==false
     printstring = printstring * string("; dval = ", logs.dc2_val[end])
   end
   if isempty(logs.IoU_train)==false
     printstring = printstring * string(";  IoUtrain:", logs.IoU_train[end,:])
   end
   if isempty(logs.IoU_val)==false
     printstring = printstring * string(";  IoUval:" , logs.IoU_val[end,:])
   end
   println(printstring)
end


text:
using example
using Test

@testset "example.jl" begin
    # Write your tests here.
end


text:
using Flux
using PyCall
const DTYPE = Float32

function reversedims(a::AbstractArray{T,N}) where {T<:AbstractFloat,N}
    permutedims(a, N:-1:1)
end

function ch2dense(link, σ = Flux.identity)
    W = link.W.array
    b = link.b.array
    Dense(W, b, σ)
end

function ch2conv(link, σ = Flux.identity)
    # get weight W and bias b
    W = reversedims(link.W.array)
    # flip kernel data
    W = W[end:-1:1, end:-1:1, :, :]
    if !isa(link.b, Nothing)
        b = reversedims(link.b.array)
    else
        b = zeros(DTYPE, size(W)[4])
    end
    pad = link.pad
    s = link.stride
    Conv(W, b, σ, pad = pad, stride = s)
end

function ch2dwconv(link, σ = Flux.identity)
    # get weight W and bias b
    W = permutedims(link.W.array, (4, 3, 1, 2))
    # flip kernel data
    W = W[end:-1:1, end:-1:1, :, :]
    if !isa(link.b, Nothing)
        b = reversedims(link.b.array)
    else
        b = zeros(DTYPE, size(W)[4])
    end
    pad = link.pad
    s = link.stride
    DepthwiseConv(W, b, σ, pad = pad, stride = s)
end

function ch2bn(link, λ = Flux.identity)
    β = link.beta.array
    γ = link.gamma.array
    ϵ = DTYPE(link.eps)
    μ = link.avg_mean
    σ² = link.avg_var
    sz = size(μ)[1]
    momentum = DTYPE(0.1) #unused parameter

    #set `false` so as to BatchNorm works as test mode
    activate = false
    bn = BatchNorm(λ, β, γ, μ, σ², ϵ, momentum, false)
    return bn
end


text:
# Run a script that supposedly creates a new fake disk here.

# Now, synthesize a new set of channel maps
# First, create the necessary input files
run(`DJ_initialize.jl`)

# make plots of the structure
run(`DJ_plot_structure.jl`)

# run the synthesis
run(`DJ_synthesize_model.jl`)

# plot the channel maps
run(`DJ_plot_chmaps.jl`)

# plot the moments
# run(`plot_moments.jl`)

# Try downsampling the model to the baseline locations
run(`DJ_write_model.jl`)

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
module sbp

include("config.jl")
include("grid.jl")
include("operator.jl")
include("metrics.jl")
include("sparse.jl")
include("vtk.jl")
include("test.jl")
include("source.jl")
include("table.jl")
include("lmv.jl")
include("mms.jl")
include("boundary_operator.jl")

# Time integrators
include("time_integrators/lsrk4.jl")
include("time_integrators/staggered_rk4.jl")

# Operators
include("../operators/oreilly_petersson_2019/operators.jl")
include("../operators/strand_1994/operators.jl")

# Staggered operators and discretizations
include("staggered.jl")
include("staggered/acoustic.jl")

# Collocated operators and discretizations
include("collocated.jl")
include("collocated/acoustic.jl")

end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
using Plasm
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using DataStructures

(V, EV) = ([0.4147 0.3513 1.0727 1.0207 0.9985 0.9819 0.9589 0.9551 0.2472 0.3382 0.412 0.6323 0.6883 0.6905 0.7134 0.7493 0.9692 0.96975 0.97765 0.98495 0.98675 1.0061 1.0118 1.0196 1.0212 0.7513 0.8803 0.9053 1.1325 0.8074 0.7588 0.8152 1.0346 1.071 1.2059 1.145 1.12855 1.1198 1.0817 1.0793 0.9135 0.894 0.8909 0.5633 0.5419 0.5337 0.8759 0.8777 0.8805 0.6816 0.678 0.67025 0.6533 0.6396 0.6011 0.5156 0.4937 0.4369 0.3842 0.1423 0.188 0.3885 0.3111 0.3739 0.5279 0.6445 0.4494 0.4046 0.3577 0.5681 0.6194 0.6245 0.6312 0.6457 0.7953 0.7642 0.4687 0.5068 0.4716 0.4184 0.9952 0.9689 0.8672 0.8052 0.8873 1.0581 0.7346 0.6812 0.6802 0.6416 0.5996 0.2089 0.1971 0.1637 0.5901 0.6027 0.6054 0.6599 0.6653 0.6873 0.8794 0.8735 0.8456 0.8338 1.0287 0.9359 0.9698 0.9867 1.0316 1.0844 1.1472 1.192 1.1772 1.1558 1.116 0.9926 1.0564 0.8367 0.6039 0.5565 0.7989 0.868 0.8747 0.8885 0.8912 1.1651 1.2045 1.2104 1.0543 1.039 1.0099 0.4669 0.6515 0.6552 0.6721 0.6753 0.7547 0.7621 0.995 1.0185 1.1216 1.1358 0.1423 0.1606 0.1941 0.5786 0.4388 0.8556 0.8585 0.9516 1.0371 1.1562 0.396 0.3752 1.0705 1.09785 1.1442 1.1637 0.3416 0.2369 0.0492 0.0 1.2505 1.0215 0.9192 0.944 0.9664 1.0162 0.4171 0.3993 0.3606 0.254 0.2095 0.0984 1.2328 1.1039 1.033 1.0084 0.3604 0.3398 0.2462 0.2247 0.2226 0.042 0.5514 0.5684 0.5978 0.66255 0.664 0.6678 0.742 0.4943 0.5545 0.579 0.5919 0.225 0.2006 0.194 0.2391 0.2201 0.7876 0.7645 0.756 0.7482 0.5815 0.5367 0.5016 0.4217 0.3506 0.3548 0.3662 1.0513 0.965 0.9079 0.7189 0.7069 0.6592 0.6005 0.4101 0.4085 0.3773 1.0137 1.00855 0.8911 0.7751 0.7739 0.7613 0.5301 0.4518 0.3103 1.0211 1.1883 1.2607 1.0758 1.0546 1.0327 0.89135 0.8844 0.9508 0.9783 0.9904 1.0214 0.8913 0.8323 0.9119 0.7358 0.7209 0.5202 0.9087 0.8283 1.2488 0.9538 1.1125 0.9496 0.8673 0.4435 0.688 0.7065 0.9402 1.0538 1.0825 1.33 0.6925 0.7525 0.8367 0.3249 0.2557 0.1853 0.9102 1.00805 1.1588 0.6854 0.6691 0.555 0.2412 0.2551 0.2636 0.4074 0.5142 0.3093 0.1165 1.0201 0.7699 0.7295 0.7237 0.8186 0.6238 0.3453 0.3965 0.7716 0.704 0.8378 0.7839 0.6726 0.6944 0.6354 0.8668 0.8745 0.8808 0.9726 1.1682 0.4632 0.419 0.6471 0.6784 0.8065 0.4551 0.7609 0.859 0.7564 1.1209 0.8271 0.54 0.715 0.5082 0.4986 0.3978 0.1423 0.3351 0.4001 0.6971 0.2674 0.2182 0.685 0.6432 0.6226 0.7084 0.6228 0.7452; 0.4904 0.5168 0.9814 0.9485 0.9344 0.9239 0.9094 0.907 1.0989 0.9581 0.8439 0.5094 0.4607 0.4581 0.438 0.4064 0.9969 0.9937 0.9482 0.9061 0.8957 0.7843 0.7513 0.7061 0.6971 1.0753 1.0703 1.0693 1.0605 0.126 0.139 0.0947 0.1359 0.1427 0.4411 0.5693 0.604 0.62245 0.7027 0.7077 0.3983 0.3189 0.3063 1.0409 0.9443 0.9075 0.6632 0.6004 0.4992 0.7729 0.7751 0.78 0.7907 0.7992 0.8234 0.87695 0.8907 0.9263 0.9593 0.7807 0.8221 1.0037 0.3333 0.328 0.3149 1.0586 0.8413 0.7914 0.7392 0.8981 0.7818 0.7703 0.7551 0.7223 0.4622 0.4731 0.577 0.4283 0.4861 0.5735 0.7424 0.741 0.7355 0.7322 0.7214 0.7624 0.9219 0.83545 0.8342 0.7718 0.7041 0.7098 0.7733 0.9528 0.2515 0.2782 0.2841 0.40015 0.4117 0.4586 0.1956 0.2002 0.222 0.2312 0.7575 0.8673 0.3942 0.3897 0.3778 0.3638 0.3471 0.3352 0.2386 0.1788 0.0674 0.8776 0.7764 0.969 0.7684 0.7276 0.8568 0.921 0.9272 0.94 0.4834 0.6226 0.6426 0.6456 0.2169 0.154 0.0345 0.2752 0.1923 0.1906 0.183 0.1816 0.1459 0.1426 0.1651 0.1592 0.1333 0.1297 0.5157 0.4583 0.3529 0.4549 0.4957 0.2896 0.2918 0.3629 0.4282 0.5192 0.7772 0.7094 0.9229 0.8917 0.8389 0.8167 0.0061 0.07755 0.20555 0.2391 0.3048 0.0 0.9417 0.9705 0.9966 1.0545 0.6364 0.8568 0.4165 0.4388 0.4481 0.4713 0.2841 0.379 0.4312 0.4493 0.3815 0.3698 0.3169 0.3048 0.3036 0.2015 0.3637 0.3353 0.2863 0.17805 0.17595 0.17035 0.0459 0.4395 0.3708 0.3428 0.3281 0.2774 0.5455 0.6184 0.0357 0.3895 0.8019 0.7491 0.7296 0.7117 0.0774 0.1127 0.1404 0.2034 0.2004 0.1761 0.1098 0.3943 0.322 0.2742 0.1159 1.0629 1.0605 1.0575 0.2903 0.2855 0.1907 0.6278 0.6221 0.4932 0.366 0.5464 0.4516 0.1362 0.2323 0.406 0.4165 0.6226 0.7119 0.8832 0.8972 0.9118 1.00515 1.01 0.3652 0.2815 0.2447 0.1504 0.5897 0.6359 1.0997 0.2456 0.25 0.3092 0.4978 0.477 0.6228 0.622 0.3856 0.2593 0.1954 0.3299 0.4201 0.4269 1.0199 0.9191 0.8937 0.6741 0.7285 0.7296 0.7311 0.661 0.5967 0.5312 0.869 0.9032 0.9558 0.7758 0.7743 0.764 0.2396 0.45575 0.5878 0.3142 0.1172 0.8137 0.7443 0.2341 0.3234 0.3378 0.3399 1.0191 0.8927 0.1943 0.0961 0.7605 0.6513 0.3507 0.3443 0.3311 0.902 0.604 0.7244 0.9202 1.0818 0.7086 0.698 0.7478 1.0476 0.1394 0.1863 0.3782 0.6345 0.5357 1.0131 0.9593 0.8888 0.9262 0.3152 0.4392 0.8675 1.1366 1.0586 0.7229 0.5083 0.0418 0.196 0.4482 0.4784 0.7245 0.6409 0.2299 0.1181 0.4672 0.4508], Array{Int64,1}[[1, 2], [3, 4], [4, 5], [4, 18], [4, 260], [5, 6], [5, 19], [5, 236], [6, 7], [6, 19], [6, 20], [7, 8], [7, 20], [7, 298], [9, 10], [10, 11], [10, 61], [10, 62], [12, 13], [13, 14], [13, 100], [14, 15], [14, 100], [14, 329], [15, 16], [15, 257], [15, 314], [17, 18], [18, 19], [18, 167], [19, 166], [20, 21], [20, 270], [21, 22], [21, 269], [21, 270], [22, 23], [22, 105], [22, 106], [23, 24], [23, 82], [23, 86], [24, 25], [24, 39], [24, 300], [26, 27], [27, 28], [27, 123], [27, 299], [28, 29], [28, 237], [28, 245], [30, 31], [32, 33], [33, 34], [33, 130], [33, 131], [35, 36], [36, 37], [36, 177], [36, 232], [37, 38], [37, 125], [37, 126], [38, 39], [38, 126], [38, 223], [39, 40], [39, 301], [41, 42], [42, 43], [42, 149], [42, 150], [44, 45], [45, 46], [45, 57], [45, 66], [47, 48], [48, 49], [48, 243], [48, 244], [50, 51], [51, 52], [51, 272], [51, 273], [52, 53], [52, 88], [52, 273], [53, 54], [53, 89], [53, 90], [54, 55], [54, 71], [54, 89], [55, 56], [55, 70], [55, 71], [56, 57], [56, 217], [56, 315], [57, 58], [57, 67], [58, 59], [58, 67], [58, 303], [60, 61], [61, 93], [61, 94], [63, 64], [64, 65], [64, 180], [64, 220], [67, 68], [67, 302], [68, 69], [68, 169], [68, 170], [71, 72], [71, 119], [72, 73], [72, 90], [72, 119], [73, 74], [73, 90], [73, 91], [75, 76], [76, 77], [76, 226], [76, 227], [78, 79], [79, 80], [79, 146], [79, 147], [81, 82], [82, 83], [82, 85], [83, 84], [83, 297], [83, 298], [87, 88], [88, 89], [88, 118], [88, 295], [90, 273], [92, 93], [93, 280], [93, 281], [95, 96], [96, 97], [96, 134], [96, 187], [97, 98], [97, 187], [97, 247], [98, 99], [98, 194], [98, 257], [99, 100], [99, 193], [99, 257], [100, 328], [101, 102], [102, 103], [102, 254], [102, 255], [103, 104], [103, 137], [103, 214], [107, 108], [108, 109], [108, 150], [108, 151], [109, 110], [109, 212], [109, 213], [110, 111], [110, 176], [110, 240], [111, 112], [111, 175], [111, 176], [113, 114], [114, 115], [114, 141], [114, 163], [116, 117], [119, 120], [119, 274], [121, 122], [122, 123], [122, 298], [122, 312], [123, 124], [123, 298], [126, 127], [126, 232], [127, 128], [127, 232], [127, 233], [129, 130], [130, 140], [130, 141], [132, 133], [133, 134], [133, 188], [133, 326], [134, 135], [134, 188], [135, 136], [135, 188], [135, 305], [136, 137], [136, 190], [136, 305], [137, 138], [137, 215], [139, 140], [140, 241], [140, 242], [141, 142], [141, 164], [143, 144], [144, 145], [144, 173], [144, 174], [148, 149], [149, 214], [149, 283], [150, 213], [150, 239], [151, 152], [151, 176], [151, 177], [153, 154], [155, 156], [156, 157], [156, 261], [156, 311], [157, 158], [157, 261], [157, 262], [159, 160], [160, 161], [160, 182], [160, 199], [161, 162], [161, 183], [161, 184], [165, 166], [166, 167], [166, 237], [167, 168], [167, 259], [171, 172], [172, 173], [172, 181], [172, 276], [173, 183], [173, 197], [176, 253], [177, 178], [177, 231], [179, 180], [180, 181], [180, 230], [181, 182], [181, 275], [182, 183], [182, 200], [183, 196], [185, 186], [186, 187], [186, 194], [186, 313], [187, 248], [188, 189], [188, 206], [189, 190], [190, 191], [190, 304], [190, 327], [192, 193], [193, 194], [193, 256], [194, 195], [197, 198], [197, 267], [197, 268], [201, 202], [202, 203], [202, 264], [202, 290], [203, 204], [203, 264], [203, 265], [205, 206], [206, 207], [206, 320], [207, 208], [207, 229], [207, 279], [209, 210], [210, 211], [210, 288], [210, 289], [213, 214], [213, 240], [214, 254], [216, 217], [217, 218], [219, 220], [220, 221], [220, 229], [222, 223], [223, 224], [223, 252], [224, 225], [224, 249], [224, 250], [228, 229], [229, 278], [232, 251], [234, 235], [235, 236], [235, 261], [235, 270], [236, 260], [236, 270], [237, 238], [240, 241], [240, 254], [241, 254], [241, 282], [246, 247], [247, 283], [247, 305], [257, 258], [260, 261], [260, 271], [263, 264], [264, 291], [266, 267], [267, 318], [267, 319], [273, 296], [276, 277], [276, 322], [276, 323], [283, 284], [283, 293], [284, 285], [284, 293], [284, 294], [286, 287], [292, 293], [293, 306], [305, 321], [307, 308], [309, 310], [316, 317], [324, 325]])

model = V,EV
V,EVs = Lar.biconnectedComponent(model)
EW = convert(Lar.Cells, cat(EVs))
VV = [[k] for k=1:size(V,2)]

hpc = Plasm.lar2hpc(V,EW)
cyan = Plasm.color("cyan")
out = cyan(hpc)
Plasm.view( out )

colors = ["orange", "red", "green", "blue", "cyan", "magenta", "yellow",
"white", "black", "purple", "gray", "brown"]

hpccolors = [Plasm.color(key) for key in colors]
comps = [Plasm.lar2hpc(V,EV) for EV in EVs]::Array{Plasm.Hpc,1}

using PyCall
p = PyCall.pyimport("pyplasm")

Plasm.view(p["STRUCT"]([hpccolors[k](comps[k]) for k=1:length(comps)]))


text:
#=
complex_number:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-07
=#

module CompNum
    # expressions
    expr1 = 1 + 2im
    expr2 = 2 - 3im

    # methods
    function get_expr1()
        return expr1
    end

    function get_expr2()
        return expr2
    end

    function multiply(expr1, expr2)
        return expr1 * expr2
    end

    function get_real(expr)
        return real(expr)
    end

    function get_imag(expr)
        return imag(expr)
    end

    function get_conj(expr)
        return conj(expr)
    end

    function get_abs(expr)
        return abs(expr)
    end

    function main()
        expr_1 = get_expr1()
        expr_2 = get_expr2()
        println("Expr1 = ", expr_1)
        println("Expr2 = ", expr_2)
        println("Expr1 * Expr2 = ", multiply(expr_1, expr_2))
        println("Real of expr1 = ", get_real(expr_1))
        println("Imag of expr1 = ", get_imag(expr_1))
        println("Conj of expr1 = ", get_conj(expr_1))
        println("Abs of expr1 = ", get_abs(expr_1))
        println("Real of expr2 = ", get_real(expr_2))
        println("Imag of expr2 = ", get_imag(expr_2))
        println("Conj of expr2 = ", get_conj(expr_2))
        println("Abs of expr2 = ", get_abs(expr_2))
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .CompNum
    CompNum.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
using RunTests
using Base.Test

@testmodule BoundaryReflection begin
	# @parameterize 1:4 function (x)
	# 	@test
	# end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

module PolicySetDefinitionsClient

using Random
using Dates
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_PolicySetDefinitionsApi.jl")

# export models
export convert, CloudError
export convert, ErrorAdditionalInfo
export convert, ErrorResponse
export convert, ParameterDefinitions
export convert, ParameterDefinitionsValue
export convert, ParameterValues
export convert, ParameterValuesValue
export convert, PolicyDefinitionGroup
export convert, PolicyDefinitionReference
export convert, PolicySetDefinition
export convert, PolicySetDefinitionListResult
export convert, PolicySetDefinitionProperties

# export operations
export convert, PolicySetDefinitionsApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
__precompile__()

module Requires

include("init.jl")
include("require.jl")

function __init__()
    push!(package_callbacks, loadpkg)
end

end # module


text:
#=
Additional methods to allow Reproject to work.
=#

using Reproject

"""
    img_proj, mask = reproject(img_in::AstroImageMat, img_out::AstroImageMat)

Reprojects the AstroImageMat `img_in` to the coordinates of `img_out`
according to the WCS information/headers using interpolation.
"""
function Reproject.reproject(img_in::AstroImageMat, img_out::AstroImageMat)
    data_out, mask = reproject(img_in, img_out)
    # TODO: should copy the WCS headers from img_out and the remaining
    # headers from img_in.
    return copyheaders(img_in, data_out)
end

function Reproject.parse_input_data(input_data::AstroImageMat, hdu)
    input_data, input_data.wcs
end
function Reproject.parse_output_projection(output_data::AstroImageMat, hdu)
    output_data.wcs, size(output_data)
end
function Reproject.pad_edges(array_in::AstroImageMat{T}) where {T}
    image = Matrix{T}(undef, size(array_in)[1] + 2, size(array_in)[2] + 2)
    image[2:end-1,2:end-1] = array_in
    image[2:end-1,1] = array_in[:,1]
    image[2:end-1,end] = array_in[:,end]
    image[1,:] = image[2,:]
    image[end,:] = image[end-1,:]
    return AstroImageMat(image, headers(array_in), wcs(array_in))
end


text:
# ---
# title: 224. Basic Calculator
# id: problem224
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math, Stack
# link: <https://leetcode.com/problems/basic-calculator/description/>
# hidden: true
# ---
#
# Implement a basic calculator to evaluate a simple expression string.
#
# The expression string may contain open `(` and closing parentheses `)`, the
# plus `+` or minus sign `-`, **non-negative** integers and empty spaces ` `.
#
# **Example 1:**
#
#
#
#     Input: "1 + 1"
#     Output: 2
#
#
# **Example 2:**
#
#
#
#     Input: " 2-1 + 2 "
#     Output: 3
#
# **Example 3:**
#
#
#
#     Input: "(1+(4+5+2)-3)+(6+8)"
#     Output: 23
#
# **Note:**
#
#   * You may assume that the given expression is always valid.
#   * **Do not** use the `eval` built-in library function.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Base: String

include("Calculator.jl")
include("asset.jl")

x = interest(10, 0.1)
println("Interest: ", x)
sum_res = sum(x, 10)
println("Sum: ", sum_res)


text:
using Plots; pyplot()
using Random

using Pkg
using Distributions

using GLM
using DataFrames

# The regression function
m(x) = 3*x.-3*x.^2+x.^3-0.1*x.^4

Random.seed!(79015)
nn = 30
xx = range(0, stop = 3, length = nn)
yy = m(xx)+rand(Normal(0, 0.15), nn)

plot(xx, yy,
    seriestype=:scatter,
    label = "",
    m = (0.9, :RebeccaPurple, 5),
    background_color = :Lavender
    )
plot!(m, 0, 3,
    line = (:tomato, 1, 2),
    label = "Regression function")

# We adjust a first order linear estimator
model1 = lm(hcat(repeat([1], outer = 30), xx), yy)

# applies a function to each value of an array and returns a new array
# containing the resulting values:
roundmap(vect, b) = map(x -> round(x, digits = b), vect)

roundmap(coef(model1), 3)

#0.582
#-0.019
#the regression function is 3*x.-3*x.^2+x.^3-0.1*x.^4

coef(model1)

plot(xx, yy,
    seriestype=:scatter,
    label = "",
    m = (0.9, :RebeccaPurple, 5),
    background_color = :Lavender
    )
plot!(m, 0, 3,
    line = (:tomato, 1, 2),
    label = "Regression function")
plot!(x -> sum([1, x].*coef(model1)), 0, 3,
    label = "Linear fit",
    line = (:RebeccaPurple, 1, 2))

#A quadratic fit
model2 = lm(hcat(repeat([1], outer = nn), xx, map(x -> x^2, xx)), yy)

roundmap(coef(model2), 3)

#0.599
#-0.055
#0.012
#the regression function is 3*x.-3*x.^2+x.^3-0.1*x.^4

plot(xx, yy,
    seriestype=:scatter,
    label = "",
    m = (0.9, :RebeccaPurple, 5),
    background_color = :Lavender
    )
plot!(m, 0, 3,
    line = (:tomato, 1, 2),
    label = "Regression function")
plot!(x -> sum([1, x, x^2].*coef(model2)), 0, 3,
    label = "Quadratic fit",
    line = (:RebeccaPurple, 1, 2))

#A cubic fit
model3 = lm(hcat(repeat([1], outer = nn), xx, map(x -> x^2, xx),
    map(x -> x^3, xx)), yy)

roundmap(coef(model3), 3)

#0.014
#2.507
#-2.16
#0.483
#the regression function is 3*x.-3*x.^2+x.^3-0.1*x.^4

plot(xx, yy,
    seriestype=:scatter,
    label = "",
    m = (0.9, :RebeccaPurple, 5),
    background_color = :Lavender
    )
plot!(m, 0, 3,
    line = (:tomato, 1, 2),
    label = "Regression function")
plot!(x -> sum([1, x, x^2, x^3].*coef(model3)), 0, 3,
    label = "Cubic fit",
    line = (:RebeccaPurple, 1, 2))

model4 = lm(hcat(repeat([1], outer = nn), xx, map(x -> x^2, xx),
    map(x -> x^3, xx), map(x -> x^4, xx)), yy)

roundmap(coef(model4), 3)

#-0.069
#3.16
#-3.178
#1.016
#-0.089
#the regression function is 3*x.-3*x.^2+x.^3-0.1*x.^4

plot(xx, yy,
    seriestype=:scatter,
    label = "",
    m = (0.9, :RebeccaPurple, 5),
    background_color = :Lavender
    )
plot!(m, 0, 3,
    line = (:tomato, 1, 2),
    label = "Regression function")
plot!(x -> sum([1, x, x^2, x^3, x^4].*coef(model4)), 0, 3,
    label = "4 degree fit",
    line = (:RebeccaPurple, 1, 2))

A = [repeat([1], outer = 10) xx[1:10]]
B = [repeat([1], outer = 10) xx[11:20]]
C = [repeat([1], outer = 10) xx[21:30]]
Z = [repeat([0], outer = 10) repeat([0], outer = 10)]

X = [A Z Z;
    Z B Z;
    Z Z C]

bets = X'X\X'yy

plot(xx, yy,
    seriestype=:scatter,
    label = "",
    m = (0.9, :RebeccaPurple, 5),
    background_color = :Lavender
    )
plot!(m, 0, 3,
    line = (:tomato, 1, 2),
    label = "Regression function")
plot!(x -> sum([1, x].*bets[1:2]),
    0, 1,
    label = "",
    line = (:RebeccaPurple, 1, 2))
plot!(x -> sum([1, x].*bets[3:4]),
    1, 2,
    label = "",
    line = (:RebeccaPurple, 1, 2))
plot!(x -> sum([1, x].*bets[5:6]),
    2, 3,
    label = "",
    line = (:RebeccaPurple, 1, 2))

k1 = 1
k2 = 2
Kt =[1 k1 -1 -k1 0 0;
    0 0 1 k2 -1 -k2]

aux1 = X'X\Kt'
aux2 = (Kt*((X'X)\Kt'))\(Kt*bets)
tilbets = bets-aux1*aux2

plot(xx, yy,
    seriestype=:scatter,
    label = "",
    m = (0.9, :RebeccaPurple, 5),
    background_color = :Lavender
    )
plot!(m, 0, 3,
    line = (:tomato, 1, 2),
    label = "Regression function")
plot!(x -> sum([1, x].*tilbets[1:2]),
    0, 1,
    label = "",
    line = (:RebeccaPurple, 1, 2))
plot!(x -> sum([1, x].*tilbets[3:4]),
    1, 2,
    label = "",
    line = (:RebeccaPurple, 1, 2))
plot!(x -> sum([1, x].*tilbets[5:6]),
    2, 3,
    label = "",
    line = (:RebeccaPurple, 1, 2))

A = hcat(repeat([1], outer = 10), xx[1:10],
        map(x -> x^2, xx[1:10]), map(x -> x^3, xx[1:10]))
B = hcat(repeat([1], outer = 10), xx[11:20],
        map(x -> x^2, xx[11:20]), map(x -> x^3, xx[11:20]))
C = hcat(repeat([1], outer = 10), xx[21:30],
        map(x -> x^2, xx[21:30]), map(x -> x^3, xx[21:30]))
Z = reshape(repeat(repeat([0], outer = 10), outer = 4), (10, 4))

X = [A Z Z;
    Z B Z;
    Z Z C]

bets = X'X\X'yy

Kt =[1 k1 k1^2 k1^3 -1 -k1 -k1^2 -k1^3  0   0     0     0;
     0  0    0    0  1  k2  k2^2  k2^3 -1 -k2 -k2^2 -k2^3]

aux1 = X'X\Kt'
aux2 = (Kt*((X'X)\Kt'))\(Kt*bets)
tilbets = bets-aux1*aux2

plot(xx, yy,
    seriestype=:scatter,
    label = "",
    m = (0.9, :RebeccaPurple, 5),
    background_color = :Lavender
    )
plot!(m, 0, 3,
    line = (:tomato, 1, 2),
    label = "Regression function")
plot!(x -> sum([1, x, x^2, x^3].*tilbets[1:4]),
    0, 1,
    label = "",
    line = (:RebeccaPurple, 1, 2))
plot!(x -> sum([1, x, x^2, x^3].*tilbets[5:8]),
    1, 2,
    label = "",
    line = (:RebeccaPurple, 1, 2))
plot!(x -> sum([1, x, x^2, x^3].*tilbets[9:12]),
    2, 3,
    label = "",
    line = (:RebeccaPurple, 1, 2))

Kt =[Kt;
     0  1  2k1  3k1^2  0  -1  -2k1 -3k1^2  0   0     0      0;
     0  0    0      0  0   1   2k2  3k2^2  0  -1  -2k2 -3k2^2]

aux1 = X'X\Kt'
aux2 = (Kt*((X'X)\Kt'))\(Kt*bets)
tilbets = bets-aux1*aux2

plot(xx, yy,
    seriestype=:scatter,
    label = "",
    m = (0.9, :RebeccaPurple, 5),
    background_color = :Lavender
    )
plot!(m, 0, 3,
    line = (:tomato, 1, 2),
    label = "Regression function")
plot!(x -> sum([1, x, x^2, x^3].*tilbets[1:4]),
    0, 1,
    label = "",
    line = (:RebeccaPurple, 1, 2))
plot!(x -> sum([1, x, x^2, x^3].*tilbets[5:8]),
    1, 2,
    label = "",
    line = (:RebeccaPurple, 1, 2))
plot!(x -> sum([1, x, x^2, x^3].*tilbets[9:12]),
    2, 3,
    label = "",
    line = (:RebeccaPurple, 1, 2))

Kt =[Kt;
    0 0 2 6k1 0 0 -2 -6k1 0 0  0    0;
    0 0 0   0 0 0  2  6k2 0 0 -2 -6k2]

aux1 = X'X\Kt'
aux2 = (Kt*((X'X)\Kt'))\(Kt*bets)
tilbets = bets-aux1*aux2

plot(xx, yy,
    seriestype=:scatter,
    label = "",
    m = (0.9, :RebeccaPurple, 5),
    background_color = :Lavender
    )
plot!(m, 0, 3,
    line = (:tomato, 1, 2),
    label = "Regression function")
plot!(x -> sum([1, x, x^2, x^3].*tilbets[1:4]),
    0, 1,
    label = "",
    line = (:RebeccaPurple, 1, 2))
plot!(x -> sum([1, x, x^2, x^3].*tilbets[5:8]),
    1, 2,
    label = "",
    line = (:RebeccaPurple, 1, 2))
plot!(x -> sum([1, x, x^2, x^3].*tilbets[9:12]),
    2, 3,
    label = "",
    line = (:RebeccaPurple, 1, 2),
    ylims = (-0.5, 1.5))

using Dierckx
spl = Spline1D(xx, yy, [k1, k2], k = 3, bc = "extrapolate")
spl_evaluated(x) = evaluate(spl,x)

plot(xx, yy,
    seriestype=:scatter,
    label = "",
    m = (0.9, :RebeccaPurple, 5),
    background_color = :Lavender
    )
plot!(m, 0, 3,
    line = (:tomato, 1, 2),
    label = "Regression function")
plot!(spl_evaluated, 0, 3,
    label = "",
    line = (:RebeccaPurple, 1, 2),
    background_color = :Lavender,
    ylims = (-0.5, 1.5)
    )

max_stars_repo_name: IrvingGomez/academic-hugo

max_issues_repo_name: IrvingGomez/academic-hugo

max_forks_repo_name: IrvingGomez/academic-hugo


text:
using DimReduc
using LazySets
#import DimReduc: forward_network, forward_affine_map, ishull

nnet = read_nnet("nnet/84442.nnet")

delta = 0.1

solver1 = MaxSens(delta)
solver2 = FastGrid(delta)
solver3 = SpeGuid(delta)
solver4 = FastTree(delta)

in_hyper = Hyperrectangle(fill(0.5, 8), fill(0.5, 8))
out_hyper = Hyperrectangle(fill(0.0, 2), fill(0.2, 2))
problem = Problem(nnet, in_hyper, out_hyper)

#reach = forward_network(solver3, problem.network, in_hyper)

file = open("results/group1.txt", "a")
print(file, "Test Result of Group 1: delta = $(delta)\n\n")
close(file)
#solver1

#solver2

#solver3

time3 = 0

#solve(solver3, problem)
for i = 1:2
    timed_result =@timed solve(solver3, problem)
    file = open("results/group1.txt", "a")
    print(file, "SpeGuid - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(file, " - Output: " * string(timed_result.value) * "\n")
    close(file)
    global time3 += timed_result.time
end

print(file, "Average time: " * string(time3/2) * " s\n\n")

#solver4

time4 = 0

#solve(solver4, problem)
for i = 1:2
    timed_result =@timed solve(solver4, problem)
    file = open("results/group1.txt", "a")
    print(file, "FastTree - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(file, " - Output: " * string(timed_result.value) * "\n")
    close(file)
    global time4 += timed_result.time
end
file = open("results/group1.txt", "a")
print(file, "Average time: " * string(time4/2) * " s\n\n")

close(file)


text:
"""
    JuMPNorm

Abstract type for norms to a JuMP model.
"""
abstract type JuMPNorm end

"""
    L1NormLinear
    L1NormLinear()

Add a set of linear constraints to a model to map an expression to a variable which is the L1 norm of that expression.
"""
struct L1NormLinear <: JuMPNorm end
"""
    MaxNormLinear
    MaxNormLinear()

Add a set of linear constraints to a model to map an expression to a variable which is the max norm of that expression.
"""
struct MaxNormLinear <: JuMPNorm end

"""
    ScaleNorm{S<:Real,N} <: JuMPNorm
    ScaleNorm(scale, n)

Scales result from `n` with a factor `scale`.
"""
struct ScaleNorm{S<:Real,N<:JuMPNorm} <: JuMPNorm
    scale::S
    n::N
end

"""
    SumNorm{N<:JuMPNorm} <: JuMPNorm

Sum of `ns`.
"""
struct SumNorm{N<:JuMPNorm} <: JuMPNorm
    ns::Vector{N}
end
SumNorm(ns::JuMPNorm...) = SumNorm(collect(ns))
SumNorm(sns::Pair{<:Real, <:JuMPNorm}...) = SumNorm(ScaleNorm.(first.(sns), last.(sns))...)

"""
    norm!(s::L1NormLinear, model, X)

Add a set of linear constraints to a model to map `X` to an expression `X′` which is the L1 norm of `X`.

Note that it only works for the objective function and only for minimization.
"""
function norm!(s::L1NormLinear, model, X, denom=1)
    # Use trick from http://lpsolve.sourceforge.net/5.1/absolute.htm to make min abs(expression) linear
    X′ = @variable(model, [1:length(X)])
    @constraint(model,  X .<= X′ .* denom)
    @constraint(model, -X .<= X′ .* denom)
    return @expression(model, sum(X′))
end

"""
    norm!(s::L1NormLinear, model, X)

Add a set of linear constraints to a model to map `X` to a variable `X′` which is the max norm of `X`.

Note that it only works for the objective function and only for minimization.
"""
function norm!(s::MaxNormLinear, model, X, denom=1)
    # Use trick from https://math.stackexchange.com/questions/2589887/how-can-the-infinity-norm-minimization-problem-be-rewritten-as-a-linear-program to make min abs(expression) linear
    X′ = @variable(model)
    @constraint(model,  X .<= X′ .* denom)
    @constraint(model, -X .<= X′ .* denom)
    return X′
end

function norm!(s::ScaleNorm, model, X, denom=1)
    X′ = norm!(s.n, model, X, denom)
    return @expression(model, s.scale * X′)
end

norm!(s::SumNorm, model, X, denom=1) = mapfoldl(n -> norm!(n, model, X, denom), (X′,X″) -> @expression(model, X′+X″), s.ns, init=@expression(model, 0))


text:
#=

  Set covering problem in Julia ConstraintSolver.jl

  Problem from
  Katta G. Murty: "Optimization Models for Decision Making", page 302f
  http://ioe.engin.umich.edu/people/fac/books/murty/opti_model/junior-7.pdf

  10 senators making a committee, where there must at least be one
  representative from each group:
  group:        senators:
  southern      1 2 3 4 5
  northern      6 7 8 9 10
  liberals      2 3 8 9 10
  conservative  1 5 6 7
  democrats     3 4 5 6 7 9
  republicans   1 2 8 10

  The objective is to minimize the number of senators.

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function set_covering3(problem,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   # "all_solutions"=> all_solutions,
                                                            "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 100, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))
    belongs = problem[:belongs]
    num_groups,num_senators = size(belongs)

    # Which senator to select
    @variable(model, x[1:num_senators], Bin)
    @variable(model, 0 <= min_val <= num_senators, Int)

    # cover all groups with senators
    for i in 1:num_groups
        @constraint(model, sum(x.*belongs[i,:])>=1)
    end

    @constraint(model, min_val == sum(x))

    @objective(model,Min,min_val)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                min_val_val = convert.(Integer,JuMP.value.(min_val; result=sol))
                println("min_val:$min_val_val")
                println("x:$x_val")
                println("senators selected:$([i for i in 1:num_senators if x_val[i] == 1])")
                println()
            end
        end
    else
        println("status:$status")
    end

    return status
end

problem = Dict(
    :belongs => resize_matrix([
        [1, 1, 1, 1, 1, 0, 0, 0, 0, 0],   # 1 southern
        [0, 0, 0, 0, 0, 1, 1, 1, 1, 1],   # 2 northern
        [0, 1, 1, 0, 0, 0, 0, 1, 1, 1],   # 3 liberals
        [1, 0, 0, 0, 1, 1, 1, 0, 0, 0],   # 4 conservative
        [0, 0, 1, 1, 1, 1, 1, 0, 1, 0],   # 5 democrats
        [1, 1, 0, 0, 0, 0, 0, 1, 0, 1]    # 6 republicans
    ]                                   )
)

@time set_covering3(problem,true,true)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using Documenter, Tricks

makedocs(;
    modules=[Tricks],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/oxinabox/Tricks.jl/blob/{commit}{path}#L{line}",
    sitename="Tricks.jl",
    authors="Lyndon White",
    assets=String[],
)

deploydocs(;
    repo="github.com/oxinabox/Tricks.jl",
)


text:
# Author: Gabrio Rizzuti, rizzuti.gabrio@gatech.edu
# Date: February, 2020
# Copyright: Georgia Institute of Technology, 2020
#
# Test forward modeling

using LinearAlgebra, PyPlot
push!(LOAD_PATH, "./src/"); using WaveOpsAVP

# Model
nz = 256
dz = R(5)
m_geom = ModelGeom(nz, dz)
nb = 8
m = R(1/2000^2)*ones(R, nz, nb)
m[102:end, :] .= R(1/2000^2)
# m[52:end, 2] .= R(1/2100^2)
# m[152:end, 3] .= R(1/2200^2)
M = Model(m_geom, m)

# Data parameters
nt = 1000
dt = R(0.001)
θ = range(R(0), stop = R(sqrt(3)/2), length = 10)
p = sin.(θ)*sqrt(max(m...))
dat_pars = DataParams(nt, dt, p)

# Source parameters
wav = cfreq_wavelet(50, nt, dt, (R(5), R(10), R(30), R(60)))
src_pars = SrcParams(wav, dt)

# ABC parameters
abc_size = (20, 40)
abc_fact = R(0.0001)
abc_geom = ABC_Geom(nz, abc_size, abc_fact)

# Modeling operator
acc_mod = 4
F = seismod_fun(dat_pars, src_pars; abc_geom = abc_geom, acc = acc_mod)

# Modeling
@time dat = F.eval(M; store_extracomp = true);

# # Plot
# figure(); imshow(dat.d[:, :, 1], aspect = "auto")


text:
"""
    Integer, Integer, Boolean -> Array{UInt8,1}

Convert Integer to an Array{UInt8}
"""
function int2bytes(x::Integer, l::Integer=0, little_endian::Bool=false)
    result = reinterpret(UInt8, [hton(x)])
    i = findfirst(x -> x != 0x00, result)
    if l != 0
        i = length(result) - l + 1
    end
    result = result[i:end]
    if little_endian
        reverse!(result)
    end
    return result
end

"""
    BigInt -> Array{UInt8,1}

Convert BigInt to an Array{UInt8}
"""
function int2bytes(x::BigInt)
    n_bytes_with_zeros = x.size * sizeof(Sys.WORD_SIZE)
    uint8_ptr = convert(Ptr{UInt8}, x.d)
    n_bytes_without_zeros = 1

    if ENDIAN_BOM == 0x04030201
        # the minimum should be 1, else the result array will be of
        # length 0
        for i in n_bytes_with_zeros:-1:1
            if unsafe_load(uint8_ptr, i) != 0x00
                n_bytes_without_zeros = i
                break
            end
        end

        result = Array{UInt8}(undef, n_bytes_without_zeros)

        for i in 1:n_bytes_without_zeros
            @inbounds result[n_bytes_without_zeros + 1 - i] = unsafe_load(uint8_ptr, i)
        end
    else
        for i in 1:n_bytes_with_zeros
            if unsafe_load(uint8_ptr, i) != 0x00
                n_bytes_without_zeros = i
                break
            end
        end

        result = Array{UInt8}(undef, n_bytes_without_zeros)

        for i in 1:n_bytes_without_zeros
            @inbounds result[i] = unsafe_load(uint8_ptr, i)
        end
    end
    return result
end

"""
    Array{UInt8,1}, Boolean -> Integer

Convert Array{UInt8,1} to an Integer
"""
function bytes2int(x::Array{UInt8,1}, little_endian::Bool=false)
    if length(x) > 8
        bytes2big(x)
    else
        missing_zeros = div(Sys.WORD_SIZE, 8) -  length(x)
        if missing_zeros > 0
            if little_endian
                for i in 1:missing_zeros
                    push!(x,0x00)
                end
            else
                for i in 1:missing_zeros
                    pushfirst!(x,0x00)
                end
            end
        end
        if ENDIAN_BOM == 0x04030201 && little_endian
        elseif ENDIAN_BOM == 0x04030201 || little_endian
            reverse!(x)
        end
        return reinterpret(Int, x)[1]
    end
end

"""
    Array{UInt8,1} -> BigInt

Convert Array{UInt8,1} to an Integer
"""
function bytes2big(x::Array{UInt8,1})
    hex = bytes2hex(x)
    return parse(BigInt, hex, base=16)
end

# Alternative implementation

# function bytes2int2(x::Array{UInt8,1})
#     if isempty(x)
#         return 0
#     end
#     if length(x) > div(Sys.WORD_SIZE, 8)
#         T = BigInt
#     else
#         T = Int
#     end
#     result = zero(T)
#     if ENDIAN_BOM == 0x01020304
#         reverse!(x)
#     end
#     for c in x
#         result <<= 8
#         result += c
#     end
#     return result
# end
#
# function faster_bytes2big(x::Array{UInt8,1})
#     xsize = cld(length(x), Base.GMP.BITS_PER_LIMB / 8)
#     if ENDIAN_BOM == 0x04030201
#         reverse!(x)
#     end
#     result = Base.GMP.MPZ.realloc2(xsize * Base.GMP.BITS_PER_LIMB)
#     result.size = xsize
#     unsafe_copyto!(result.d, convert(Ptr{Base.GMP.Limb}, pointer(x)), xsize)
#     return result
# end


text:
using DiffPrecTest
using JLD

coverage = 0
coverage_stud = 0

fName = ARGS[1]
alpha = parse(Float64, ARGS[2])

NUM_REP = 1000

res = load(fName, "results")

for rep=1:NUM_REP
    global coverage
    global coverage_stud

    CI = simulCI(res[rep], alpha)
    coverage += all(0 .<= CI[:, 2]) * all(0 .>= CI[:, 1]) ? 1 : 0

    CI = simulCIstudentized(res[rep], alpha)
    coverage_stud += all(0 .<= CI[:, 2]) * all(0 .>= CI[:, 1]) ? 1 : 0
end

@show coverage / NUM_REP
@show coverage_stud / NUM_REP


text:
using Test
using BKTrees

@testset "Constructors" begin
        T = Int
        val::T = 1
        vals = T[1,2,3]
        foo(x,y)=abs(x-y)

        node = Node{T}()
        @test node isa Node{T}
        @test node.item == nothing

        node = Node(val)
        @test node isa Node{T} && node.item==val &&
            isempty(node.children)

        bkt = BKTree{T}()
        @test bkt isa BKTree{T}
        @test bkt.root isa Node{T}
        @test bkt.root.item == nothing

        bkt = BKTree{T}(foo)
        @test bkt isa BKTree{T}
        @test bkt.root isa Node{T}
        @test bkt.root.item == nothing
        @test bkt.f == foo
        @test isempty(bkt.root.children)

        node = Node(val)
        bkt = BKTree(node)
        @test bkt isa BKTree{T}
        @test bkt.root isa Node{T}
        @test bkt.root.item == val
        @test isempty(bkt.root.children)

        bkt = BKTree(val)
        @test bkt isa BKTree{T}
        @test bkt.root isa Node{T}
        @test bkt.root.item == val
        @test isempty(bkt.root.children)

        bkt = BKTree(vals)
        @test bkt isa BKTree{T}
        @test bkt.root isa Node{T}
        @test bkt.f == BKTrees.hamming_distance
        @test !isempty(bkt.root.children)

        bkt = BKTree(foo, vals)
        @test bkt isa BKTree{T}
        @test bkt.root isa Node{T}
        @test bkt.f == foo
        @test !isempty(bkt.root.children)
end

@testset "add!, find" begin
    bkt = BKTree{Int}()
    for val in [0,1,2,5]
        add!(bkt, val)
    end
    target = 3
    result = find(bkt, target, 10, k=5)
	for (i, (dist, item)) in enumerate(result)
		i==1 && @test dist == 1.0 && item==1
		i==2 && @test dist == 1.0 && item==2
		i==3 && @test dist == 2.0 && item==0
		i==4 && @test dist == 2.0 && item==5
	end
end

# show methods
@testset "Show methods" begin
    buf = IOBuffer()
    T = Int
    try
        node = Node{T}()
        show(buf, node)
        @test true
    catch
        @test false
    end
    try
        tree = BKTree{T}()
        show(buf, tree)
        @test true
    catch
        @test false
    end
end


text:
# subsets(xs) = IterTools.subsets(xs)
# subsets(xs, k) = IterTools.subsets(xs, k)


text:
module Ionisation
import SpecialFunctions: gamma
import GSL: hypergeom
import HDF5
import Pidfile: mkpidlock
import Logging: @info
import Luna.PhysData: c, ħ, electron, m_e, au_energy, au_time, au_Efield, wlfreq
import Luna.PhysData: ionisation_potential, quantum_numbers
import Luna: Maths
import Luna: @hlock

"""
    ionrate_fun!_ADK(ionpot::Float64, threshold=true)
    ionrate_fun!_ADK(material::Symbol)

Return a closure `ionrate!(out, E)` which calculates the ADK ionisation rate for the electric
field `E` and places the result in `out`. If `threshold` is true, use [`ADK_threshold`](@ref)
to avoid calculation below floating-point precision. If `cycle_average` is `true`, calculate
the cycle-averaged ADK ionisation rate instead.
"""
function ionrate_fun!_ADK(ionpot::Float64, threshold=true; cycle_average=false)
    nstar = sqrt(0.5/(ionpot/au_energy))
    cn_sq = 2^(2*nstar)/(nstar*gamma(nstar+1)*gamma(nstar))
    ω_p = ionpot/ħ
    ω_t_prefac = electron/sqrt(2*m_e*ionpot)

    if threshold
        thr = ADK_threshold(ionpot)
    else
        thr = 0
    end

    # Zenghu Chang: Fundamentals of Attosecond Optics (2011) p. 184
    # Section 4.2.3.1 Cycle-Averaged Rate
    # ̄w_ADK(Fₐ) = √(3/π) √(Fₐ/F₀) w_ADK(Fₐ) where Fₐ is the field amplitude
    Ip_au = ionpot / au_energy
    F0_au = (2Ip_au)^(3/2)
    F0 = F0_au*au_Efield
    avfac = sqrt.(3/(π*F0))

    ionrate! = let nstar=nstar, cn_sq=cn_sq, ω_p=ω_p, ω_t_prefac=ω_t_prefac, thr=thr
        function ir(E)
            if abs(E) >= thr
                r = (ω_p*cn_sq*
                    (4*ω_p/(ω_t_prefac*abs(E)))^(2*nstar-1)
                    *exp(-4/3*ω_p/(ω_t_prefac*abs(E))))
                if cycle_average
                    r *= avfac*sqrt(abs(E))
                end
                return r
            else
                return zero(E)
            end
        end
        function ionrate!(out, E)
            out .= ir.(E)
        end
    end

    return ionrate!
end

function ionrate_fun!_ADK(material::Symbol; kwargs...)
    return ionrate_fun!_ADK(ionisation_potential(material); kwargs...)
end

function ionrate_ADK(IP_or_material, E; kwargs...)
    out = zero(E)
    ionrate_fun!_ADK(IP_or_material; kwargs...)(out, E)
    return out
end

function ionrate_ADK(IP_or_material, E::Number; kwargs...)
    out = [zero(E)]
    ionrate_fun!_ADK(IP_or_material; kwargs...)(out, [E])
    return out[1]
end

"""
    ADK_threshold(ionpot)

Determine the lowest electric field strength at which the ADK ionisation rate for the
ionisation potential `ionpot` is non-zero to within 64-bit floating-point precision.
"""
function ADK_threshold(ionpot)
    out = [0.0]
    ADKfun = ionrate_fun!_ADK(ionpot, false)
    E = 1e3
    while out[1] == 0
        E *= 1.01
        ADKfun(out, [E])
    end
    return E
end

"""
    ionrate_fun!_PPTaccel(material::Symbol, λ0; kwargs...)
    ionrate_fun!_PPTaccel(ionpot::Float64, λ0, Z, l; kwargs...)

Create an accelerated (interpolated) PPT ionisation rate function.
"""
function ionrate_fun!_PPTaccel(material::Symbol, λ0; kwargs...)
    n, l, Z = quantum_numbers(material)
    ip = ionisation_potential(material)
    ionrate_fun!_PPTaccel(ip, λ0, Z, l; kwargs...)
end

function ionrate_fun!_PPTaccel(ionpot::Float64, λ0, Z, l; kwargs...)
    E, rate = makePPTcache(ionpot, λ0, Z, l; kwargs...)
    return makePPTaccel(E, rate)
end

"""
    ionrate_fun!_PPTcached(material::Symbol, λ0; kwargs...)
    ionrate_fun!_PPTcached(ionpot::Float64, λ0, Z, l; kwargs...)

Create a cached (saved) interpolated PPT ionisation rate function. If a saved lookup table
exists, load this rather than recalculate.

# Keyword arguments
- `N::Int`: Number of samples with which to create the `CSpline` interpolant.
- `Emax::Number`: Maximum field strength to include in the interpolant.
- `cachedir::String`: Path to the directory where the cache should be stored and loaded from.
    Defaults to \$HOME/.luna/pptcache

Other keyword arguments are passed on to [`ionrate_fun_PPT`](@ref)
"""
function ionrate_fun!_PPTcached(material::Symbol, λ0; kwargs...)
    n, l, Z = quantum_numbers(material)
    ip = ionisation_potential(material)
    ionrate_fun!_PPTcached(ip, λ0, Z, l; kwargs...)
end

function ionrate_fun!_PPTcached(ionpot::Float64, λ0, Z, l;
                                sum_tol=1e-4, cycle_average=false, N=2^16, Emax=nothing,
                                cachedir=joinpath(homedir(), ".luna", "pptcache"))
    h = hash((ionpot, λ0, Z, l, sum_tol, cycle_average, N, Emax))
    fname = string(h, base=16)*".h5"
    fpath = joinpath(cachedir, fname)
    lockpath = joinpath(cachedir, "pptlock")
    isdir(cachedir) || mkpath(cachedir)
    if isfile(fpath)
        @info "Found cached PPT rate for $(ionpot/electron) eV, $(λ0*1e9) nm"
        pidlock = mkpidlock(lockpath)
        rate = loadPPTaccel(fpath)
        close(pidlock)
        return rate
    else
        E, rate = makePPTcache(ionpot::Float64, λ0, Z, l;
                               sum_tol=sum_tol, cycle_average, N=N, Emax=Emax)
        @info "Saving PPT rate cache for $(ionpot/electron) eV, $(λ0*1e9) nm in $cachedir"
        pidlock = mkpidlock(lockpath)
        if isfile(fpath) # makePPTcache takes a while - has another process saved first?
            rate = loadPPTaccel(fpath)
            close(pidlock)
            return rate
        end
        @hlock HDF5.h5open(fpath, "cw") do file
            file["E"] = E
            file["rate"] = rate
        end
        close(pidlock)
        return makePPTaccel(E, rate)
    end
end

function loadPPTaccel(fpath)
    isfile(fpath) || error("PPT cache file $fpath not found!")
    E, rate = @hlock HDF5.h5open(fpath, "r") do file
        (read(file["E"]), read(file["rate"]))
    end
    makePPTaccel(E, rate)
end

function makePPTcache(ionpot::Float64, λ0, Z, l;
                      sum_tol=1e-4, cycle_average=false, N=2^16, Emax=nothing)
    Emax = isnothing(Emax) ? 2*barrier_suppression(ionpot, Z) : Emax

    # ω0 = 2π*c/λ0
    # Emin = ω0*sqrt(2m_e*ionpot)/electron/0.5 # Keldysh parameter of 0.5
    Emin = Emax/5000

    E = collect(range(Emin, stop=Emax, length=N));
    @info "Pre-calculating PPT rate for $(ionpot/electron) eV, $(λ0*1e9) nm"
    rate = ionrate_PPT(ionpot, λ0, Z, l, E; sum_tol=sum_tol, cycle_average);
    @info "PPT pre-calcuation done"
    return E, rate
end

"""
    barrier_suppression(ionpot, Z)

Calculate the barrier-suppresion **field strength** for the ionisation potential `ionpot`
and charge state `Z`.
"""
function barrier_suppression(ionpot, Z)
    Ip_au = ionpot / au_energy
    ns = Z/sqrt(2*Ip_au)
    Z^3/(16*ns^4) * au_Efield
end

"""
    keldysh(material, λ, E)

Calculate the Keldysh parameter for the given `material` at wavelength `λ` and electric field
strength `E`.
"""
function keldysh(material, λ, E)
    Ip_au = ionisation_potential(material)/au_energy
    E_au = E/au_Efield
    ω0_au = wlfreq(λ)*au_time
    ω0_au*sqrt(2Ip_au)/E_au
end

"""
    ionfrac(rate, E, δt)

Given an ionisation rate function `rate` and an electric field array `E` sampled with time
spacing `δt`, calculate the ionisation fraction as a function of time on the same time axis.

The function `rate` should have the signature `rate!(out, E)` and place its results into
`out`, like the functions returned by e.g. `ionrate_fun!_ADK` or `ionrate_fun!_PPTcached`.
"""
function ionfrac(rate, E, δt)
    r = similar(E)
    frac = similar(E)
    rate(r, E)
    Maths.cumtrapz!(frac, r, δt)
    @. frac = 1 - exp(-frac)
end

function makePPTaccel(E, rate)
    cspl = Maths.CSpline(E, log.(rate); bounds_error=true)
    Emin = minimum(E)
    # Interpolating the log and re-exponentiating makes the spline more accurate
    ir(E) = abs(E) <= Emin ? 0.0 : exp(cspl(abs(E)))
    function ionrate!(out, E)
        out .= ir.(E)
    end
end

function ionrate_fun!_PPT(args...)
    ir = ionrate_fun_PPT(args...)
    function ionrate!(out, E)
        out .= ir.(E)
    end
    return ionrate!
end

"""
    ionrate_fun_PPT(ionpot::Float64, λ0, Z, l; sum_tol=1e-4, cycle_average=false)

Create closure to calculate PPT ionisation rate.

# Keyword arguments
- `sum_tol::Number`: Relative tolerance used to truncate the infinite sum.
- `cycle_average::Bool`: If `false` (default), calculate the cycle-averaged rate

# References
[1] Ilkov, F. A., Decker, J. E. & Chin, S. L.
Ionization of atoms in the tunnelling regime with experimental evidence
using Hg atoms. Journal of Physics B: Atomic, Molecular and Optical
Physics 25, 4005–4020 (1992)

[2] 1.Bergé, L., Skupin, S., Nuter, R., Kasparian, J. & Wolf, J.-P.
Ultrashort filaments of light in weakly ionized, optically transparent
media. Rep. Prog. Phys. 70, 1633–1713 (2007)
(Appendix A)
"""
function ionrate_fun_PPT(ionpot::Float64, λ0, Z, l; sum_tol=1e-4, cycle_average=false)
    Ip_au = ionpot / au_energy
    ns = Z/sqrt(2Ip_au)
    ls = ns-1
    Cnl2 = 2^(2ns)/(ns*gamma(ns + ls + 1)*gamma(ns - ls))

    ω0 = 2π*c/λ0
    ω0_au = au_time*ω0
    E0_au = (2*Ip_au)^(3/2)

    ionrate = let ω0_au=ω0_au, Cnl2=Cnl2, ns=ns, sum_tol=sum_tol
        function ionrate(E)
            E_au = abs(E)/au_Efield
            γ = ω0_au*sqrt(2Ip_au)/E_au
            γ2 = γ*γ
            β = 2γ/sqrt(1 + γ2)
            α = 2*(asinh(γ) - γ/sqrt(1+γ2))
            Up_au = E_au^2/(4*ω0_au^2)
            Uit_au = Ip_au + Up_au
            v = Uit_au/ω0_au
            ret = 0
            divider = 0
            for m = -l:l
                divider += 1
                mabs = abs(m)
                flm = ((2l + 1)*factorial(l + mabs)
                    / (2^mabs*factorial(mabs)*factorial(l - mabs)))
                # Following 5 lines are [1] eq. 8 and lead to identical results:
                # G = 3/(2γ)*((1 + 1/(2γ2))*asinh(γ) - sqrt(1 + γ2)/(2γ))
                # Am = 4/(sqrt(3π)*factorial(mabs))*γ2/(1 + γ2)
                # lret = sqrt(3/(2π))*Cnl2*flm*Ip_au
                # lret *= (2*E0_au/(E_au*sqrt(1 + γ2))) ^ (2ns - mabs - 3/2)
                # lret *= Am*exp(-2*E0_au*G/(3E_au))
                # [2] eq. (A14)
                lret = 4sqrt(2)/π*Cnl2
                lret *= (2*E0_au/(E_au*sqrt(1 + γ2))) ^ (2ns - mabs - 3/2)
                lret *= flm/factorial(mabs)
                lret *= exp(-2v*(asinh(γ) - γ*sqrt(1+γ2)/(1+2γ2)))
                lret *= Ip_au * γ2/(1+γ2)
                # Remove cycle average factor, see eq. (2) of [1]
                if !cycle_average
                    lret *= sqrt(π*E0_au/(3E_au))
                end
                k = ceil(v)
                n0 = ceil(v)
                sumfunc = let k=k, β=β, m=m
                    function sumfunc(x, n)
                        diff = n-v
                        return x + exp(-α*diff)*φ(m, sqrt(β*diff))
                    end
                end
                # s, success, steps = Maths.aitken_accelerate(
                #     sumfunc, 0, n0=n0, rtol=sum_tol, maxiter=Inf)
                s, success, steps = Maths.converge_series(
                    sumfunc, 0, n0=n0, rtol=sum_tol, maxiter=Inf)
                lret *= s
                ret += lret
            end
            return ret/(au_time*divider)
        end
    end
    return ionrate
end

"""
    φ(m, x)

Calculate the φ function for the PPT ionisation rate.

Note that w_m(x) in [1] and φ_m(x) in [2] look slightly different but
are in fact identical.
"""
function φ(m, x)
    mabs = abs(m)
    return (exp(-x^2)
            * sqrt(π)
            * x^(mabs+1)
            * gamma(mabs+1)
            * hypergeom(1/2, 3/2 + mabs, x^2)
            / (2*gamma(3/2 + mabs)))
end

function ionrate_fun_PPT(material::Symbol, λ0; kwargs...)
    n, l, Z = quantum_numbers(material)
    ip = ionisation_potential(material)
    return ionrate_fun_PPT(ip, λ0, Z, l; kwargs...)
end

function ionrate_PPT(ionpot, λ0, Z, l, E; kwargs...)
    return ionrate_fun_PPT(ionpot, λ0, Z, l; kwargs...).(E)
end

function ionrate_PPT(material::Symbol, λ0, E; kwargs...)
    n, l, Z = quantum_numbers(material)
    ip = ionisation_potential(material)
    return ionrate_PPT(ip, λ0, Z, l, E; kwargs...)
end

end


text:
module MiniKalman

using GaussianDistributions, FillArrays
using GaussianDistributions: dim, logpdf
using StaticArrays
using Statistics, Random, LinearAlgebra
using Unitful: ustrip, unit   # could be @required ?
using QuickTypes: roottypeof  # TODO: get rid of this dependency
using Random: GLOBAL_RNG
using DocStringExtensions
@template (FUNCTIONS, METHODS) =
    """
    $(SIGNATURES)
    $(DOCSTRING)
    """

export kfilter, kalman_filter, white_noise1, white_noise2,
    kalman_smoother, kalman_sample, no_noise, log_likelihood, cumulative_log_likelihood,
    sample_and_recover, optimize, marginal, kalman_predicted, kalman_filtered

include("utils.jl")
include("filtering.jl")
include("smoothing.jl")
include("optim.jl")
include("sampling.jl")

end  # module


text:
using RegNeuralDE, OrdinaryDiffEq, Flux, CUDA, Tracker
CUDA.allowscalar(false)

mlp_dynamics = TDChain(Dense(3, 10, CUDA.tanh), Dense(11, 2)) |> track |> gpu

x_ = rand(2, 1) |> gpu

# Unregularized Neural ODE
node_unreg = TrackedNeuralODE(
    mlp_dynamics,
    [0.0f0, 1.0f0],
    true,
    false,
    Tsit5(),
    save_everystep = false,
    reltol = 1.4f-8,
    abstol = 1.4f-8,
    save_start = false,
)

## Forward Pass
@code_warntype node_unreg(x_ |> track)

## Backward Pass
@code_warntype Tracker.gradient(p -> sum(node_unreg(x_ |> track, p)[1]), node_unreg.p)

# Regularized Neural ODE -- Error Estimates
node_errreg = TrackedNeuralODE(
    mlp_dynamics,
    [0.0f0, 1.0f0],
    true,
    true,
    Tsit5(),
    save_everystep = false,
    reltol = 1.4f-8,
    abstol = 1.4f-8,
    save_start = false,
)

## Forward Pass
@code_warntype node_errreg(
    x_ |> track;
    func = (u, t, integrator) -> integrator.EEst * integrator.dt,
)

## Backward Pass
@code_warntype Tracker.gradient(
    p -> begin
        res, _, sv = node_errreg(
            x_ |> track,
            p;
            func = (u, t, integrator) -> integrator.EEst * integrator.dt,
        )
        sum(res) + sum(sv.saveval)
    end,
    node_errreg.p,
)

# Regularized Neural ODE -- Stiffness Estimates
node_stiffreg = TrackedNeuralODE(
    mlp_dynamics,
    [0.0f0, 1.0f0],
    true,
    true,
    AutoTsit5(Tsit5()),
    save_everystep = false,
    reltol = 1.4f-8,
    abstol = 1.4f-8,
    save_start = false,
)

## Forward Pass
@code_warntype node_stiffreg(
    x_ |> track;
    func = (u, t, integrator) -> abs(integrator.eigen_est * integrator.dt),
)

## Backward Pass
@code_warntype Tracker.gradient(
    p -> begin
        res, _, sv = node_stiffreg(
            x_ |> track,
            p;
            func = (u, t, integrator) -> abs(integrator.eigen_est * integrator.dt),
        )
        sum(res) + sum(filter(x -> !iszero(x), sv.saveval))
    end,
    node_stiffreg.p,
)


text:
@testset "Mass of hermitian lattices" begin
  Qx, x = PolynomialRing(FlintQQ, "x")
  f = x-1
  K, a = NumberField(f, "a", cached = false)
  Kt, t = PolynomialRing(K, "t")
  g = t^2+(2)
  E, b = NumberField(g, "b", cached = false)
  D = matrix(E, 3, 3, [(1), 0, 0, 0, (1), 0, 0, 0, (1)])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [(-2)*b+(-2), b+(6), 0]), map(E, [0, (1), (1)]), map(E, [b+(-6), (-6)*b+(6), 0])]
  L = hermitian_lattice(E, generators = gens, gram_ambient_space = D)
  m = @inferred mass(L)
  @test m == 1//32

  Qx, x = PolynomialRing(FlintQQ, "x")
  f = x^2-3
  K, a = NumberField(f, "a", cached = false)
  Kt, t = PolynomialRing(K, "t")
  g = t^2+(1)
  E, b = NumberField(g, "b", cached = false)
  D = matrix(E, 3, 3, [(1), 0, 0, 0, (1), 0, 0, 0, (1)])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [(-1//2*a-1//2)*b+(1//2*a+5//2), (-1), 0]), map(E, [(a-9)*b+(15*a+3), 0, (-a)*b+(8*a+13)]), map(E, [(-1610*a+2990)*b+(2070*a-2990), 0, (-68*a+393)*b+(19*a+462)])]
  L = hermitian_lattice(E, generators = gens, gram_ambient_space = D)
  m = @inferred mass(L)
  @test m == 1//108

  Qx, x = PolynomialRing(FlintQQ, "x")
  f = x^4-x^3-4*x^2+4*x+1
  K, a = NumberField(f, "a", cached = false)
  Kt, t = PolynomialRing(K, "t")
  g = t^2+(a^3 - 1*a^2 - 4*a + 5)
  E, b = NumberField(g, "b", cached = false)
  D = matrix(E, 3, 3, [(1), 0, 0, 0, (1), 0, 0, 0, (1)])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [(1), 0, 0]), map(E, [0, (1), 0]), map(E, [0, 0, (1)])]
  L = hermitian_lattice(E, generators = gens, gram_ambient_space = D)
  m = @inferred mass(L)
  @test m == 1//10125
end