text:
import Base: +, -, *, /, ==, transpose, zero
import LinearAlgebra: dot, norm, normalize!
import SparseArrays: nnz
import SimilaritySearch: cosine_distance, angle_distance
export DVEC, compute_bow # add!

const DVEC{Ti,Tv<:Real} = Dict{Ti,Tv}
const BOW = DVEC{Symbol,Int}
const SVEC = DVEC{Int,Float64}

nnz(dvec::DVEC) = length(dvec)

"""
    compute_bow(tokenlist::AbstractVector{Symbol}, voc::DVEC)

Updates a DVEC using the given list of tokens; returns `voc`
"""
function compute_bow(tokenlist::AbstractVector{Symbol}, voc::BOW)
    maxfreq = 0
    for sym in tokenlist
        m = get(voc, sym, 0) + 1
        voc[sym] = m
        maxfreq = max(m, maxfreq)
    end

    voc, maxfreq
end

# these are needed to call `compute_bow` for symbol's list but also for simplicity of the API
compute_bow(tokenlist::AbstractVector{Symbol}) = compute_bow(tokenlist, BOW())

"""
    normalize!(bow::DVEC)

Inplace normalization of `bow`
"""
function normalize!(bow::DVEC{Ti,Tv}) where {Ti,Tv<:AbstractFloat}
    s = 1.0 / norm(bow)
    for (k, v) in bow
        bow[k] = convert(Tv, v * s)
    end

    bow
end

function normalize!(bow::DVEC{Ti,Tv}) where {Ti,Tv<:Integer}
    s = 1.0 / norm(bow)
    for (k, v) in bow
        bow[k] = round(T, v * s)
    end

    bow
end

function normalize!(matrix::AbstractVector{DVEC})
    for bow in matrix
        normalize!(bow)
    end
end

function dot(a::DVEC, b::DVEC)::Float64
    if length(b) < length(a)
        a, b = b, a  # a must be the smallest bow
    end

    s = 0.0
    for (k, v) in a
        w = convert(Float64, get(b, k, 0))
        s += v * w
    end

    s
end

function norm(a::DVEC)::Float64
    s = 0.0
    for w in values(a)
        s += w * w
    end

    sqrt(s)
end

function zero(::Type{DVEC{Ti,Tv}}) where {Ti,Tv<:Real}
    DVEC{Ti,Tv}()
end

## inplace sum
function add!(a::DVEC{Ti,Tv}, b::DVEC{Ti,Tv}) where {Ti,Tv<:Real}
    for (k, w) in b
        if w != 0
            a[k] = get(a, k, zero(Tv)) + w
        end
    end

    a
end

function add!(a::DVEC{Ti,Tv}, b::AbstractSparseArray) where {Ti,Tv<:Real}
    for (k, w) in zip(b.nzind, b.nzval)
        if w != 0
            a[k] = get(a, k, zero(Tv)) + w
        end
    end

    a
end

function add!(a::DVEC{Ti,Tv}, b::Pair{Ti,Tv}) where {Ti,Tv<:Real}
    k, w = b
    a[k] = get(a, k, zero(Tv)) + w
    a
end

function +(a::DVEC{Ti,Tv}, b::DVEC{Ti,Tv}) where {Ti,Tv<:Real}
    if length(a) < length(b)
        a, b = b, a  # a must be the largest bow
    end

    c = copy(a)
    for (k, w) in b
        if w != 0
            c[k] = get(c, k, zero(Tv)) + w
        end
    end

    c
end

function +(a::DVEC, b::Pair)
    c = copy(a)
    add!(c, b)
end

## definitions for substraction
function -(a::DVEC{Ti,Tv}, b::DVEC{Ti,Tv}) where {Ti,Tv<:Real}
    c = copy(a)
    for (k, w) in b
        if w != 0
            c[k] = get(c, k, zero(Tv)) - w
        end
    end

    c
end

## definitions for product

function *(a::DVEC{Ti,Tv}, b::DVEC{Ti,Tv}) where {Ti,Tv<:Real}
    if length(b) < length(a)
        a, b = b, a  # a must be the smallest bow
    end

    c = copy(a)
    for k in keys(a)
        w = get(b, k, zero(Tv))
        if w == 0
            delete!(c, k)
        else
            c[k] = convert(Tv, c[k] * w)
        end
    end

    c
end

function *(a::DVEC, b::F) where F<:Real
    c = copy(a)
    for (k, v) in a
        c[k] = convert(F, v * b)
    end

    c
end

function *(b::F, a::DVEC) where F<:Real
    a * b
end

function /(a::DVEC, b::F) where F<:Real
    c = copy(a)
    for (k, v) in a
        c[k] = convert(F, v / b)
    end

    c
end

"""
    cosine_distance(a::DVEC, b::DVEC)::Float64

Computes the cosine_distance between two weighted bags

It supposes that bags are normalized (see `normalize!` function)

"""
function cosine_distance(a::DVEC, b::DVEC)::Float64
    return 1.0 - dot(a, b)
end

const π_2 = π / 2
"""
    angle_distance(a::DVEC, b::DVEC)::Float64

Computes the angle  between two weighted bags

It supposes that all bags are normalized (see `normalize!` function)

"""
function angle_distance(a::DVEC, b::DVEC)::Float64
    d = dot(a, b)

    if d <= -1.0
        return π
    elseif d >= 1.0
        return 0.0
    elseif d == 0  # turn around for zero vectors, in particular for denominator=0
        return π_2
    else
        return acos(d)
    end
end

function cosine(a::DVEC, b::DVEC)::Float64
    return dot(a, b) # * a.invnorm * b.invnorm # it is already normalized
end


text:
"""
  get_linking_matrix(
    m::AbstractBlockNLPModel
  )

Returns the concatenated linking matrix ``A = [A_1, \\ldots, A_B]``

# Arguments

- `m::AbstractBlockNLPModel`: BlockNLPModel whose block matrix is to be extracted.
"""
function get_linking_matrix(m::AbstractBlockNLPModel)
    nb = m.problem_size.block_counter
    A = spzeros(m.problem_size.link_counter, m.problem_size.var_counter)
    for i = 1:nb
        for j = 1:length(m.linking_constraints)
            A[m.linking_constraints[j].idx, m.blocks[i].var_idx] =
                m.linking_constraints[j].linking_blocks[i]
        end
    end
    return A
end

"""
  get_linking_matrix_blocks(
    m::AbstractBlockNLPModel
  )

Returns a vector of linking matrix blocks ``[A_1, \\ldots, A_B]``

# Arguments

- `m::AbstractBlockNLPModel`: name of the BlockNLPModel whose block matrices are to be extracted.
"""
function get_linking_matrix_blocks(m::AbstractBlockNLPModel)
    nb = m.problem_size.block_counter
    A = [spzeros(m.problem_size.link_counter, m.blocks[i].meta.nvar) for i = 1:nb]
    for i = 1:nb
        for j = 1:length(m.linking_constraints)
            A[i][m.linking_constraints[j].idx, :] =
                m.linking_constraints[j].linking_blocks[i]
        end
    end
    return A
end

"""
    get_rhs_vector(
        m::AbstractBlockNLPModel
    )

Returns the concatenated RHS vector for all the linking constraints ``b``.

# Arguments

- `m::AbstractBlockNLPModel`: name of the BlockNLPModel whose RHS vector is to be extracted.
"""
function get_rhs_vector(m::AbstractBlockNLPModel)
    b = zeros(m.problem_size.link_counter)
    for j = 1:length(m.linking_constraints)
        b[m.linking_constraints[j].idx] = m.linking_constraints[j].rhs_vector
    end
    return b
end

"""
    n_constraints(
        m::AbstractBlockNLPModel
    )

Returns the total number of constraints in a BlockNLPModel.

# Arguments

- `m::AbstractBlockNLPModel`: name of the BlockNLPModel.
"""
function n_constraints(m::AbstractBlockNLPModel)
    return m.problem_size.link_counter + m.problem_size.con_counter
end

"""
    update_nnzh(
        meta::NLPModelMeta,
        new_nnzh::Int
    )

Updates the nnzh field of NLPModelMeta.

# Arguments

- `meta::NLPModelMeta`: name of the BlockNLPModel.
- `new_nnzh`: new nnzh value to update the meta with.
"""
function update_nnzh(meta::NLPModelMeta, new_nnzh::Int)
    field_names = [
        :nvar,
        :x0,
        :lvar,
        :uvar,
        :nlvb,
        :nlvo,
        :nlvc,
        :ncon,
        :y0,
        :lcon,
        :ucon,
        :nnzo,
        :nnzj,
        :lin_nnzj,
        :nln_nnzj,
        :nnzh,
        :lin,
        :minimize,
        :islp,
        :name,
    ]
    field_values = [getfield(meta, field) for field in field_names]

    new_meta = Dict(zip(field_names, field_values))
    new_meta[:nnzh] = new_nnzh

    nvar = new_meta[:nvar]
    pop!(new_meta, :nvar)

    return NLPModels.NLPModelMeta(nvar; new_meta...)
end

"""
    AugmentedHessianInfo

A data type to store information required to compute augmented hessian for a NLP block.
"""
mutable struct AugmentedHessianInfo
    augmented_hessian_struct::Tuple{AbstractVector,AbstractVector}
    block_hessian_struct::Tuple{AbstractVector,AbstractVector}
    ATA::Tuple{AbstractVector,AbstractVector,AbstractVector}
end

"""
    get_augmented_hessian_coord!(nlp::AbstractNLPModel,
        x::AbstractVector,
        vals::AbstractVector,
        obj_weight::Number;
        y::Union{AbstractVector, Nothing} = nothing,
    )

Returns the Hessian for an augmented subproblem as a sparse matrix.

# Arguments
- `m::AbstractNLPModel`: the subproblem
- `x::Union{AbstractVector, Nothing}`: current primal solution (optional)
- `obj_weight::Union{AbstractVector, Nothing}`: objective weight
- `y::Union{AbstractVector, Nothing}`: vector of dual variables
- `vals::Union{AbstractVector, Nothing}`: nonzero values of the Hessian matrix
"""
function get_augmented_hessian_coord!(
    nlp::AbstractNLPModel,
    x::AbstractVector,
    vals::AbstractVector,
    obj_weight::Number;
    y::Union{AbstractVector,Nothing} = nothing,
)
    # initialize vector indices
    main_idx = 1
    sub_idx1 = 1
    sub_idx2 = 1

    # assign pointers
    aug_hess = nlp.hess_info.augmented_hessian_struct
    blk_hess = nlp.hess_info.block_hessian_struct
    aug_term = nlp.hess_info.ATA
    if nlp.subproblem.meta.ncon > 0
        blk_hess_values =
            hess_coord(nlp.subproblem.problem_block, x, y, obj_weight = obj_weight)
    else
        blk_hess_values =
            hess_coord(nlp.subproblem.problem_block, x, obj_weight = obj_weight)
    end

    for (i, j) in zip(aug_hess[1], aug_hess[2])
        if (
            blk_hess[1][sub_idx1] == i &&
            blk_hess[2][sub_idx1] == j &&
            aug_term[1][sub_idx2] == i &&
            aug_term[2][sub_idx2] == j
        )
            vals[main_idx] = blk_hess_values[sub_idx1] + aug_term[3][sub_idx2]
            main_idx += 1
            sub_idx1 += 1
            sub_idx2 += 1
        elseif (blk_hess[1][sub_idx1] == i && blk_hess[2][sub_idx1] == j)
            vals[main_idx] = blk_hess_values[sub_idx1]
            main_idx += 1
            sub_idx1 += 1
        elseif (aug_term[1][sub_idx2] == i && aug_term[2][sub_idx2] == j)
            vals[main_idx] = aug_term[3][sub_idx2]
            main_idx += 1
            sub_idx2 += 1
        else # To-Do: probably not required, check and remove.
            vals[main_idx] = 0.0
            main_idx += 1
        end
    end
end

ext: jl


text:
# 2D CA using Agents.jl
module CA2D
using Agents

mutable struct Cell{T<:Integer, Y<:AbstractString} <: AbstractAgent
  id::T
  pos::Tuple{T, T}
  status::Y
end

"""
    build_model(;rules::Tuple, dims=(100,100), Moore=true)

Builds a 2D cellular automaton. `rules` is of type `Tuple{Integer,Integer,Integer}`. The numbers are DSR (Death, Survival, Reproduction). Cells die if the number of their living neighbors are <D, survive if the number of their living neighbors are <=S, come to life if their living neighbors are as many as R. `dims` is the x and y size a grid. `Moore` specifies whether cells should connect to their diagonal neighbors.
"""
function build_model(;rules::Tuple, dims=(100,100), Moore=true)
  space = Space(dims, moore=Moore)
  properties = Dict(:rules => rules, :Moore=>Moore)
  model = ABM(Cell, space; properties = properties, scheduler=by_id)
  nnodes = dims[1]*dims[2]
  for n in 1:nnodes
    add_agent!(Cell(n, vertex2coord(n, dims), "0"), (n,1), model)
  end
  return model
end

"""
    periodic_neighbors(pos::Tuple{Integer, Integer}, dims::Tuple{Integer, Integer})

Returns the the row and column numbers of the rows and columns before and after the pos.
"""
function periodic_neighbors(pos::Tuple{Integer, Integer}, dims::Tuple{Integer, Integer})
  if pos[1] == 1
    xbefore = dims[1]
    if dims[1] == 1
      xafter = 1
    else
      xafter = pos[1] + 1
    end
  elseif pos[1] == dims[1]
    xbefore = pos[1] - 1
    xafter = 1
  else
    xbefore = pos[1] - 1
    xafter = pos[1] + 1
  end
  if pos[2] == 1
    ybefore = dims[2]
    if dims[2] == 1
      yafter = 1
    else
      yafter = pos[2] + 1
    end
  elseif pos[2] == dims[2]
    ybefore = pos[2] - 1
    yafter = 1
  else
    ybefore = pos[2] - 1
    yafter = pos[2] + 1
  end
  return (xbefore, ybefore), (xafter, yafter)
end

function ca_step!(model)
  agentnum = nagents(model)
  new_status = Array{String}(undef, agentnum)
  for agid in 1:agentnum
    agent = model.agents[agid]
    coord = agent.pos
    center = agent.status
    before, after = periodic_neighbors(coord, model.space.dimensions)
    right = model.agents[coord2vertex((after[1], coord[2]), model)].status
    left = model.agents[coord2vertex((before[1], coord[2]), model)].status
    top = model.agents[coord2vertex((coord[1], after[2]), model)].status
    bottom = model.agents[coord2vertex((coord[1], after[2]), model)].status
    topright = model.agents[coord2vertex((after[1], after[2]), model)].status
    topleft = model.agents[coord2vertex((before[1], after[2]), model)].status
    bottomright = model.agents[coord2vertex((after[1], before[2]), model)].status
    bottomleft = model.agents[coord2vertex((before[1], before[2]), model)].status

    if model.properties[:Moore]
      nstatus = [topleft, top, topright, left, right, bottomleft, bottom, bottomright]
      nlive = length(findall(a->a=="1", nstatus))
    else
      nstatus = [top, left, right, bottom]
      nlive = length(findall(a->a=="1", nstatus))
    end

    if agent.status == "1"
      if nlive < model.properties[:rules][1]
        new_status[agid] = "0"
      elseif nlive > model.properties[:rules][2]
        new_status[agid] = "0"
      else
        new_status[agid] = "1"
      end
    else
      if nlive == model.properties[:rules][3]
        new_status[agid] = "1"
      else
        new_status[agid] = "0"
      end
    end
  end
  for ss in 1:agentnum
    model.agents[ss].status = new_status[ss]
  end
end

"""
    ca_run(model::ABM, runs::Integer)

Runs a 2D cellular automaton.
"""
function ca_run(model::ABM, runs::Integer, plot_CA2Dgif::T; nodesize=2) where T<: Function
  data = step!(model, dummystep, ca_step!, 0, [:pos, :status], step0=true)
  anim = plot_CA2Dgif(data, nodesize=nodesize)
  for r in 1:runs
    data = step!(model, dummystep, ca_step!, 1, [:pos, :status], step0=false)
    anim = plot_CA2Dgif(data, anim=anim, nodesize=nodesize)
  end
  return anim
end

end  # module


text:
module GalleryWaveguide

using NonlinearEigenproblems.NEPCore
using NonlinearEigenproblems.NEPTypes
using NonlinearEigenproblems.LinSolvers
using NonlinearEigenproblems.Gallery
using SparseArrays
using LinearAlgebra
using Statistics
using FFTW
using IterativeSolvers
using LinearMaps

export wep_generate_preconditioner

# Specializalized NEPs
export WEP
export WEP_FD

# We overload these
import ..Gallery.nep_gallery
import ..NEPCore.compute_Mlincomb
import ..NEPCore.compute_Mder
import ..LinSolvers.lin_solve
import ..LinSolvers.create_linsolver

import Base.size
import SparseArrays.issparse
import Base.*
import Base.eltype
import LinearAlgebra.ldiv!

abstract type WEP <: NEP end

include("waveguide/Waveguide.jl")
include("waveguide/waveguide_FD.jl")
include("waveguide/waveguide_preconditioner.jl")

"""
    nep_gallery(WEP, kwargs...)

Create the NEP associated with finite difference discretization of the **Waveguide Eigenvalue Problem (WEP)** found in both
* E. Ringh, and G. Mele, and J. Karlsson, and E. Jarlebring, Sylvester-based preconditioning for the waveguide eigenvalue problem, Linear Algebra Appl., 2018
* E. Jarlebring, and G. Mele, and O. Runborg, The waveguide eigenvalue problem and the tensor infinite Arnoldi method, SIAM J. Sci. Comput., 2017

# keyword parameters:
* nx::Integer = 3 * 5 * 7,         disctretization points in x-direction
* nz::Integer = 3 * 5 * 7,         disctretization points in z-direction
* benchmark_problem = 'TAUSCH',    which waveguide (TAUSCH, JARLEBRING)
* neptype::String = 'WEP',         NEP-format (SPMF, SPMF_PRE, WEP) later format recommended
* delta = 0.1,                     Slack from the absorbing boundary conditions

See also: [`WEPLinSolverCreator`](@ref)
"""
function nep_gallery(::Type{T}; nx::Integer = 3*5*7, nz::Integer = 3*5*7, benchmark_problem::String = "TAUSCH", neptype::String = "WEP",  delta::Number = 0.1) where T<:WEP

    FFTW.set_num_threads(Sys.CPU_THREADS)

    waveguide = uppercase(benchmark_problem)
    neptype = uppercase(neptype)
    if !isodd(nz)
        error("Variable nz must be odd! You have used nz = ", nz, ".")
    end

    # Generate the matrices and formulate the problem is the sought format
    if (neptype == "SPMF") || (neptype == "SPMF_PRE")
        K, hx, hz, Km, Kp = generate_wavenumber_fd( nx, nz, waveguide, delta)
        Dxx, Dzz, Dz = generate_fd_interior_mat( nx, nz, hx, hz)
        C1, C2T = generate_fd_boundary_mat( nx, nz, hx, hz)
        if(neptype == "SPMF_PRE")
            pre_Schur_fact = true
        else
            pre_Schur_fact = false
        end
        nep = assemble_waveguide_spmf_fd(nx, nz, hx, Dxx, Dzz, Dz, C1, C2T, K, Km, Kp, pre_Schur_fact)

    elseif (neptype == "WEP")
        K, hx, hz, Km, Kp = generate_wavenumber_fd( nx, nz, waveguide, delta)
        Dxx, Dzz, Dz = generate_fd_interior_mat( nx, nz, hx, hz)
        C1, C2T = generate_fd_boundary_mat( nx, nz, hx, hz)
        nep = WEP_FD(nx, nz, hx, hz, Dxx, Dzz, Dz, C1, C2T, K, Km, Kp)

    else
        error("The NEP-type '", neptype, "' is not supported.")
    end

    return nep
end

end


text:
using StanSample, Test

ProjDir = @__DIR__
cd(ProjDir)

gq = "
  data {
    int<lower=0> N;
    int<lower=0> y[N];
  }
  parameters {
    real<lower=0.00001> Theta;
  }
  model {
    Theta ~ gamma(4, 1000);
    for (n in 1:N)
      y[n] ~ exponential(Theta);
  }
  generated quantities{
    real y_pred;
    y_pred = exponential_rng(Theta);
  }
";

gq_data = Dict(
  "N" => 3,
  "y" => [100, 950, 450]
);

stanmodel = SampleModel("Generate_quantities", gq);

(sample_file, log_file) = stan_sample(stanmodel; data=gq_data)

if !(sample_file == Nothing)
  stan_generate_quantities(stanmodel, 1)
end


text:
include("../src/PBCM.jl")

cell_general = cell()
cell_design_op = cylindrical_cell_designer(cell_general)
cost = cost_default()
cell_general, cost = convert_all(cell_general, cost, mult)

using PyPlot

    ###################################################### USD/kWh vs Thickness Cathode  ######################################################

param_x = []
param_y = []

cost.cell_costs.cathode.AM[1]   = 17.0
cell_general.cathode.por        = 0.4

for i in range(110, 200, length=1000)
    # cost.material.pos_AM[1] = i
    cell_general.cathode.th = converter(i, mult.pos_th)
    dollars_per_kWh, mfg_capacity = cost_calc(cell_general, cost, system="Cell", cost_verbosity=0)
    append!(param_y, dollars_per_kWh)
    append!(param_x, i)
end

# print(param_y)

clf()
plot(param_x, param_y)
xlabel(String("Cathode thickness (microns)"))
ylabel("Dollars per kWh")
title("Cost vs. Cathode Thickness")
figure(2)

max_stars_repo_name: BattModels/BatteryCost

max_issues_repo_name: battmodels/BatteryCost

max_forks_repo_name: battmodels/BatteryCost


text:
using SwiftObjectStores
using SwiftObjectStores.stat
using Base.Test
using JLD
import PyCall: PyError, @pyimport

@pyimport warnings
warnings.filterwarnings("ignore")

const serv = SwiftService(
	# Note that we will read settings from Enviroment Variables CI needs to set those
)

const verdata = "-$(VERSION.major)-$(VERSION.minor)-$(VERSION.patch)"
const osdata = string(ccall(:jl_get_UNAME, Any, ()))
const container_name = "swift.jl-testing-"*osdata*verdata

if haskey(list(serv), container_name)
    warn("Past test left mess. Deleting " * container_name * " for use in testing")
    delete(serv,container_name)
end

@testset "PreCreate" begin
    @test_throws PyError stat(serv, container_name)  # container should not exist
end

@testset "Stream User Path" begin
    const obj_name = "eg1buffer"

    eg_text = "The Julia Language: A fresh approach to technical computing. λ=x²"
    @test put_file(serv, container_name, obj_name, IOBuffer(eg_text)) != nothing
    get_stream(serv, container_name, obj_name) do fp
        @test readstring(fp) == eg_text
    end

    @test Dict(stat(serv, container_name)["items"])["Objects"] == "1"
    @test stat(serv, container_name, obj_name) !=nothing

    @test delete(serv,container_name, obj_name)!=nothing
    @test stat(serv, container_name, obj_name)["success"] == false
    @test stat(serv, container_name, "eg1buffer")["error"][:http_status] == 404
end;

@testset "File User Path" begin
    const obj_name = "eg2file"
    const file_name = joinpath(dirname(@__FILE__), "eg2file.txt")
    const eg_text= """
        Call me Ishmael. Some years ago- never mind how long precisely-
        having little or no money in my purse, and nothing particular to interest me on shore,
        I thought I would sail about a little and see the watery part of the world.
        It is a way I have of driving off the spleen and regulating the circulation."""

	open(file_name, "w") do fp
        print(fp, eg_text)
	end

    @test put_file(serv, container_name, obj_name, file_name) != nothing
    @test validate_file(serv, container_name, obj_name, file_name) #validate it as file

	@test get_file(serv, container_name, obj_name) do fn #Check reading it as a file
		open(readstring, fn, "r") == eg_text
	end

end

@testset "File User Path with PseudoDir" begin
    const obj_name = "eg3file"
    const file_name = joinpath(dirname(@__FILE__), "eg3file.txt")
    open(file_name, "w") do fp
        println(fp, """
        There now is your insular city of the Manhattoes,
        belted round by wharves as Indian isles by coral reefs
        - commerce surrounds it with her surf.
        Right and left, the streets take you waterward.
        Its extreme downtown is the battery,
        where that noble mole is washed by waves,
        and cooled by breezes,
        which a few hours previous were out of sight of land.
        Look at the crowds of water-gazers there
        """)
    end

    @test put_file(serv, container_name, "/psubdir/"*obj_name, file_name) != nothing
    @test validate_file(serv, container_name, "psubdir/"*obj_name, file_name)
end

immutable Mt
    m::String
    v::Int64
end

@testset "JLD" begin
    const obj_name = "eg2jld"
    val = Mt("Hello Dearest", 52)

    @test put_jld(serv, container_name, obj_name, mm=val) != nothing
    @test get_jld(serv, container_name, obj_name, "mm").m == val.m
    @test get_jld(serv, container_name, obj_name, "mm").v == val.v
end

@testset "Post_test" begin
    @test delete(serv,container_name)[end]["success"] == true
    @test_throws PyError stat(serv, container_name)
    @test haskey(SwiftObjectStores.list(serv), container_name) == false
end;


text:
print("gail-driver installer")
if VERSION != v"0.6.4"
  print("version is wrong. it must be 0.6.4")
  return
end

gail_root = ENV["GAIL_DRIVER_GITHUB"]

Pkg.clone(gail_root * "/Vec.jl")
Pkg.build("Vec.jl")
Pkg.clone("https://github.com/JuliaIO/VideoIO.jl")
Pkg.clone("https://github.com/JuliaGraphics/Cairo.jl")
Pkg.checkout("VideoIO.jl", "mindeps")
Pkg.checkout("Cairo.jl", "app")
Pkg.build("VideoIO.jl")
Pkg.build("Cairo.jl")
Pkg.clone(gail_root * "/AutomotiveDrivingModels.jl")
Pkg.build("AutomotiveDrivingModels.jl")
Pkg.clone(gail_root * "/AutoViz.jl")
Pkg.build("AutoViz.jl")
Pkg.clone(gail_root * "/ForwardNets.jl")
Pkg.build("ForwardNets.jl")
Pkg.clone(gail_root * "/NGSIM.jl")
Pkg.build("NGSIM.jl")
Pkg.clone(gail_root * "/Records.jl")
Pkg.build("Records.jl")
using Vec
using AutomotiveDrivingModels
using NGSIM
using AutoViz
using Records
using ForwardNets

max_stars_repo_name: Sawato/gail-driver

max_issues_repo_name: Sawato/gail-driver

max_forks_repo_name: Sawato/gail-driver


text:
using ITensors
using Test

using ChainRulesCore: rrule_via_ad

include("utils/chainrulestestutils.jl")

using Zygote: ZygoteRuleConfig

@testset "ChainRules rrules" begin
  i = Index(2, "i")
  A = randomITensor(i', dag(i))
  Ac = randomITensor(ComplexF64, i', dag(i))
  B = randomITensor(i', dag(i))
  C = ITensor(3.4)

  test_rrule(getindex, ITensor(3.4); check_inferred=false)
  test_rrule(getindex, A, 1, 2; check_inferred=false)
  test_rrule(*, A', A; check_inferred=false)
  test_rrule(*, 3.2, A; check_inferred=false)
  test_rrule(*, A, 4.3; check_inferred=false)
  test_rrule(+, A, B; check_inferred=false)
  test_rrule(prime, A; check_inferred=false)
  test_rrule(prime, A, 2; check_inferred=false)
  test_rrule(prime, A; fkwargs=(; tags="i"), check_inferred=false)
  test_rrule(prime, A; fkwargs=(; tags="x"), check_inferred=false)
  test_rrule(replaceprime, A, 1 => 2; check_inferred=false)
  test_rrule(swapprime, A, 0 => 1; check_inferred=false)
  test_rrule(addtags, A, "x"; check_inferred=false)
  test_rrule(addtags, A, "x"; fkwargs=(; plev=1), check_inferred=false)
  test_rrule(removetags, A, "i"; check_inferred=false)
  test_rrule(replacetags, A, "i" => "j"; check_inferred=false)
  test_rrule(
    swaptags, randomITensor(Index(2, "i"), Index(2, "j")), "i" => "j"; check_inferred=false
  )
  test_rrule(replaceind, A, i' => sim(i); check_inferred=false)
  test_rrule(replaceinds, A, (i, i') => (sim(i), sim(i)); check_inferred=false)
  test_rrule(swapind, A, i', i; check_inferred=false)
  test_rrule(swapinds, A, (i',), (i,); check_inferred=false)
  test_rrule(itensor, randn(2, 2), i', i; check_inferred=false)
  test_rrule(ITensor, randn(2, 2), i', i; check_inferred=false)
  test_rrule(ITensor, 2.3; check_inferred=false)
  test_rrule(dag, A; check_inferred=false)
  test_rrule(permute, A, reverse(inds(A)); check_inferred=false)

  f = function (x)
    j = Index(2, "j")
    b = itensor([0, 0, 1, 1], i, j)
    k = itensor([0, 1, 0, 0], i, j)
    T = itensor([0 x x^2 1; 0 0 sin(x) 0; 0 cos(x) 0 exp(x); x 0 0 0], i', j', i, j)
    return x * real((b' * T * k)[])
  end
  args = (0.3,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)

  f = function (x)
    j = Index(2, "j")
    b = itensor([0, 0, 1, 1], i, j)
    k = itensor([0, 1, 0, 0], i, j)
    T = itensor([0 x x^2 1; 0 0 sin(x) 0; 0 cos(x) 0 exp(x); x 0 0 0], i, j, i', j')
    return x * real((b' * T * k)[])
  end
  args = (0.3,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)

  f = x -> sin(scalar(x)^3)
  args = (C,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> sin(x[]^3)
  args = (C,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = adjoint
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = (x, y) -> (x * y)[1, 1]
  args = (A', A)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> prime(x, 2)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> x'[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> addtags(x, "x")[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x' * x)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (prime(x) * x)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> ((x'' * x') * x)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x'' * (x' * x))[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = (x, y, z) -> (x * y * z)[1, 1]
  args = (A'', A', A)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x'' * x' * x)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x''' * x'' * x' * x)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x''' * x'' * x' * x)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = (x, y) -> (x + y)[1, 1]
  args = (A, B)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x + x)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (2x)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x + 2x)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x + 2 * mapprime(x' * x, 2 => 1))[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = (x, y) -> (x * y)[]
  args = (A, δ(dag(inds(A))))
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x * x)[]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x * δ(dag(inds(x))))[]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = function (x)
    y = x' * x
    tr = δ(dag(inds(y)))
    return (y * tr)[]
  end
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = function (x)
    y = x'' * x' * x
    tr = δ(dag(inds(y)))
    return (y * tr)[]
  end
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x^2 * δ((i', i)))[1, 1]
  args = (6.2,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x^2 * δ(i', i))[1, 1]
  args = (5.2,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> itensor([x^2 x; x^3 x^4], i', i)
  args = (2.54,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> ITensor([x^2 x; x^3 x^4], i', i)
  args = (2.1,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> ITensor(x)
  args = (2.12,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = function (x)
    j = Index(2)
    T = itensor([x^2 sin(x); x^2 exp(-2x)], j', dag(j))
    return real((dag(T) * T)[])
  end
  args = (2.8,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  args = (2.8 + 3.1im,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = function f(x)
    j = Index(2)
    v = itensor([exp(-3.2x), cos(2x^2)], j)
    T = itensor([x^2 sin(x); x^2 exp(-2x)], j', dag(j))
    return real((dag(v') * T * v)[])
  end
  args = (2.8,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  args = (2.8 + 3.1im,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = function (x)
    j = Index(2)
    return real((x^3 * ITensor([sin(x) exp(-2x); 3x^3 x+x^2], j', dag(j)))[1, 1])
  end
  args = (3.4 + 2.3im,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> prime(permute(x, reverse(inds(x))))[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> prime(x; plev=1)[1, 1]
  args = (A,)
  @test_throws ErrorException f'(args...)
end


text:
x, y = rand(1000), rand(1000)

import HypothesisTests: OneSampleTTest

@test jdd(x, y) isa Vector
@test jdd(OneSampleTTest, x, y) isa OneSampleTTest


text:
#------------------------------------------------------------------------------
#  Example.jl
#
#  Test of LSTAR code used in "Carry Trade" paper.
#
#
#
#
#
#
#
#
#  Paul.Soderlind@unisg.ch   18 January 2013, to Julia Nov 2015
#------------------------------------------------------------------------------

using Statistics, Dates, DelimitedFiles, LinearAlgebra, Optim, ForwardDiff

include("jlFiles/FindNoNaNPs.jl")
include("jlFiles/NewEst3Ps.jl")
include("jlFiles/OlsPs.jl")
include("jlFiles/OlsLStar3Ps.jl")
include("jlFiles/printmat.jl")
#------------------------------------------------------------------------------

x     = readdlm("CTData.csv",',')      #date, FXV, Return_CT, SP, Ty

FXV   = x[:,2]                         #FX volatility
Re_CT = x[:,3]                         #carry trade return
SP    = x[:,4]                         #S&P return
Ty    = x[:,5]                         #treasury return

T = size(Re_CT,1)

gM = range(1,stop=4,length=15)
cM = range(-2,stop=3,length=21)

z = (FXV .- mean(FXV))/std(FXV)        #standardised regime variable
#------------------------------------------------------------------------------

gcKeep = [NaN;NaN]                       #(gamma,c) set to NaN if estimated in NLS, otherwise imposed
#gcKeep = [3.0;NaN]                      #try to uncomment this to restrict gamma to 3
#gcKeep = [NaN;1.3]                      #or this

w = zeros(T,0)                           #no regressors without "regimes": empty vector
(theta,Stdtheta,fnOutput) = OlsLStar3Ps(Re_CT,[ones(T) SP Ty],w,true,z,gM,cM,gcKeep)

if all(isnan.(gcKeep))
  println("\n","theta is [g;c;b_low;b_high;(slopes without regimes, if any)]")
  rowNames = ["γ";"c";string.("b",1:3,"_L");string.("b",1:3,"_H")]
elseif isnan.(gcKeep) == [true;false]
  println("\n","theta is [g;b_low;b_high;(slopes without regimes, if any)]")
  rowNames = ["γ";string.("b",1:3,"_L"); string.("b",1:3,"_H")]
elseif isnan.(gcKeep) == [false;true]
  println("\n","theta is [c;b_low;b_high;(slopes without regimes, if any)]")
  rowNames = ["c";string.("b",1:3,"_L");string.("b",1:3,"_H")]
end

println("Estimates and t-stats")
printTable([theta Stdtheta],["coef","t-stat"],rowNames)

println("\n","difference of slope (b, high minus low state)")
printTable(fnOutput.slopeDiff,["coef","t-stat"],string.("b",1:3))
#------------------------------------------------------------------------------

#Plot of the loss function. Comment out this if you do not have PyPlot installed.
using PyPlot
#PyPlot.svg(true)           #for ipynb notebooks

using3D()                   #active 3D plotting
fig = figure(figsize=(12,8))
ax = gca(projection="3d")
  surf(gM,cM,fnOutput.sseM',rstride=2,cstride=2,cmap=ColorMap("summer"),alpha=0.8)
  scatter(fnOutput.gcHat[1],fnOutput.gcHat[2],zs=minimum(fnOutput.sseM),s=200,color="k")
  ax.view_init(elev=20.0, azim=10)
  title("Sum of squared residuals (optimum at the dot)")
  xlabel(L"$\gamma$")
  ylabel(L"$c$")
  #display(gcf())          #uncomment in Atom/Juno
#------------------------------------------------------------------------------

max_stars_repo_name: PaulSoderlind/PaulSoderlindCode

max_issues_repo_name: PaulSoderlind/PaulSoderlindCode

max_forks_repo_name: PaulSoderlind/PaulSoderlindCode


text:
using EmpiricalRisks
using Base.Test

## Auxiliary functions

function verify_multipred(pred::PredictionModel{1,0}, θ, X::DenseMatrix)
    n = size(X, 2)
    @test ninputs(pred, X) == n
    rr = zeros(n)
    for i = 1:n
        rr[i] = predict(pred, θ, X[:,i])
    end
    @test_approx_eq predict(pred, θ, X) rr

    buffer = zeros(n)
    @test_approx_eq predict!(buffer, pred, θ, X) rr
end

function verify_multipred(pred::PredictionModel{1,1}, θ, X::DenseMatrix)
    n = size(X, 2)
    @test ninputs(pred, X) == n
    p = length(predict(pred, θ, X[:,1]))
    rr = zeros(p, n)
    for i = 1:n
        rr[:,i] = predict(pred, θ, X[:,i])
    end
    @test_approx_eq predict(pred, θ, X) rr

    buffer = zeros(p, n)
    @test_approx_eq predict!(buffer, pred, θ, X) rr
end

## Data preparation

d = 5
k = 3
n = 12

X = randn(d, n)

## Predictors

# LinearPred

a = randn()
w = randn(d)
wa = [w; a]
b = 2.5

pred = LinearPred(d)
@test inputlen(pred) == d
@test inputsize(pred) == (d,)
@test outputlen(pred) == 1
@test outputsize(pred) == ()
@test paramlen(pred) == d
@test paramsize(pred) == (d,)
@test isvalidparam(pred, w)

for i = 1:n
    x_i = X[:,i]
    @test_approx_eq predict(pred, w, x_i) dot(w, x_i)
end
verify_multipred(pred, w, X)

# AffinePred

pred = AffinePred(d, b)
@test inputlen(pred) == d
@test inputsize(pred) == (d,)
@test outputlen(pred) == 1
@test outputsize(pred) == ()
@test paramlen(pred) == d + 1
@test paramsize(pred) == (d+1,)
@test isvalidparam(pred, wa)

for i = 1:n
    x_i = X[:,i]
    @test_approx_eq predict(pred, wa, x_i) dot(w, x_i) + a * b
end
verify_multipred(pred, wa, X)

# MvLinearPred

a = randn(k)
W = randn(k, d)
Wa = [W a]
b = 2.5

pred = MvLinearPred(d, k)
@test inputlen(pred) == d
@test inputsize(pred) == (d,)
@test outputlen(pred) == k
@test outputsize(pred) == (k,)
@test paramlen(pred) == d * k
@test paramsize(pred) == (k, d)
@test isvalidparam(pred, W)

for i = 1:n
    x_i = X[:,i]
    @test_approx_eq predict(pred, W, x_i) W * x_i
end
verify_multipred(pred, W, X)

# MvAffinePred

pred = MvAffinePred(d, k, b)
@test inputlen(pred) == d
@test inputsize(pred) == (d,)
@test outputlen(pred) == k
@test outputsize(pred) == (k,)
@test paramlen(pred) == (d+1) * k
@test paramsize(pred) == (k, d+1)
@test isvalidparam(pred, Wa)

for i = 1:n
    x_i = X[:,i]
    @test_approx_eq predict(pred, Wa, x_i)  W * x_i + a * b
end
verify_multipred(pred, Wa, X)


text:
module Runtests_withPlot

using TinyModia
using Test

const  test_title = "Test TinyModia with simulation (with " * currentPlotPackage() * ")"

@testset "$test_title" begin
    include("include_all.jl")
end

end


text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.reduce
# Date: 12/15/2014

module EncounterDBN

export
    AbstractEncounterDBN,
    AddObserver,

    getInitialState,
    initialize,
    update,
    get,

    CorrAEMDBN,
    StarDBN,
    SideOnDBN,
    PairwiseCorrAEMDBN

import CommonInterfaces: initialize, update

using AbstractEncounterDBNImpl
import AbstractEncounterDBNInterfaces: getInitialState, get

using CorrAEMDBNImpl
using StarDBNImpl
#using PairwiseCorrAEMDBNImpl
using SideOnDBNImpl

end


text:
const OFFSET = 100.0 # offset in degrees!!!!
touch(string(getpid(),".pid"))
push!(LOAD_PATH, pwd())
using UnicodePlots, PyCall, Reactive, DataFrames, Light
@pyimport seabreeze
seabreeze.use("pyseabreeze")
@pyimport seabreeze.spectrometers as sb
devices = sb.list_devices()
s = sb.Spectrometer(devices[1])
maxy = s[:_dev][:interface][:_MAX_PIXEL_VALUE]
wl = s[:wavelengths]()
nwl = length(wl)
function watchit()
	watch_file("integrationtime")
	sleep(.1)
	put!(it, parse(Int,readchomp("integrationtime")))
end
function watchnr()
	watch_file("nreplicates")
	sleep(.1)
	take!(nr)
	put!(nr, parse(Int,readchomp("nreplicates")))
end
function watchr()
	watch_file("rotation")
	sleep(.1)
	put!(sp, readchomp("rotation"))
end
it = Channel{Int}(1)
nr = Channel{Int}(1)
put!(nr,1)
sp = Channel{AbstractString}(1)
yy = Input(wl)
X = DataFrame(wl = wl, l315c0 = zeros(nwl), l0c0 = zeros(nwl), l45c0 = zeros(nwl), l90c0 = zeros(nwl), l0c315 = zeros(nwl), l0c45 = zeros(nwl), dark = zeros(nwl), I = zeros(nwl), dolp = zeros(nwl), aop = zeros(nwl), docp = zeros(nwl))
writetable("data.csv", X)

@async map(readall, [`pdflatex plot.tex`, `evince plot.pdf`])

function plotit(y)
	p = scatterplot(wl,y,xlim=[300,800],ylim=[0,maxy])
	display(p)
	sleep(0.07)
end
lift(plotit,yy)
function normalize(I, dolp, aop, docp)
	I = I < 0 ? 0.0 : I > 10maxy ? 10maxy : I
	dolp = dolp < 0 ? 0.0 : dolp > 2 ? 2.0 : dolp
	docp = docp < -2 ? -2.0 : docp > 2 ? 2.0 : docp
	return (I, dolp, aop, docp)
end
function draw()
	if any(map(Bool,X[:l0c315])) && any(map(Bool,X[:l0c45]))
		s0, s1, s2, s3 = measurements2stokes(X[:l315c0] - X[:dark], X[:l0c0] - X[:dark], X[:l45c0] - X[:dark], X[:l90c0] - X[:dark], X[:l0c315] - X[:dark], X[:l0c45] - X[:dark])
	else
		s0, s1, s2, s3 = measurements2stokes(X[:l315c0] - X[:dark], X[:l0c0] - X[:dark], X[:l45c0] - X[:dark], X[:l90c0] - X[:dark])
	end
	for i = 1:nwl
		I, dolp, aop, docp = normalize(stokes2polar(s0[i], s1[i], s2[i], s3[i])...)
		X[i,:I] = I
		X[i,:dolp] = dolp
		X[i,:aop] = rad2deg(aop) + OFFSET
		X[i,:docp] = docp
	end
	writetable("data.csv", X)
	readall(`pdflatex plot.tex`)
end
function iter()
	wait(nr)
	nreps = fetch(nr)
	if isready(it)
		s[:integration_time_micros](take!(it))
	end
	y = zeros(nwl)
	for i = 1:nreps
		y .+= s[:intensities]()
	end
	y /= nreps
	push!(yy,y)
	if isready(sp)
		X[symbol(take!(sp))] = y
		draw()
	end
end
sleep(0)
s[:intensities]()
sleep(0)

@schedule begin
	while true
		watchit()
	end
end
@schedule begin
	while true
		watchr()
	end
end
@schedule begin
	while true
		watchnr()
	end
end
@schedule begin
	while true
		iter()
	end
end

sleep(60*60*24)


text:
function eftSqr{T<:StdFloat}(a::T)
    hi = a * a
    lo = fma(a, a, -hi)
    hi,lo
end

function eftIncr{T<:StdFloat}(a::T)
    b = one(T)
    hi = a + b
    t = hi - a
    lo = (a - (hi - t)) + (b - t)
    hi,lo
end

function eftDecr{T<:StdFloat}(a::T)
    b = one(T)
    hi = a - b
    t = hi - a
    lo = (a - (hi - t)) - (b + t)
    hi,lo
end

function accInv{T<:StdFloat}(b::T)
    hi = one(T)/b
    r = fma(-b,hi,one(T))
    lo = r/b
    r = fma(-b,lo,r)
    lo += r/b
    hi,lo
end

# for use when the sign of the part is all that is needed of the lo part
function accInvForSign{T<:StdFloat}(b::T)
    hi = one(T)/b          # hi has the sign of b
    loapprox = fma(-b,hi,one(T))
    hi, loapprox
end

function accSqrt{T<:StdFloat}(a::T)
    hi = sqrt(a)
    lo = fma(-hi,hi,a)
    lo = lo / (2*hi)
    hi, lo
end

# for use when the sign of the part is all that is needed of the lo part
function accSqrtForSign(a::StdFloat)
     hi = sqrt(a)
     loApprox = fma(-hi,hi,a)
     hi,loApprox
end

for fn in (:eftSqr, :eftIncr, :eftDecr, :accInv, :accSqrt, :accSqrtForSign)
  @eval begin
    ($fn){T<:Real}(a::T) = ($fn)(Float64(a))
  end
end


text:
export  KalmanFilter,
        kalmanfilter,
        kalmansmoother

struct KalmanFilter
    means::Matrix{Float}
    covs::Array{Float,3}
    # Storage for recycling in kalman smoother
    means_::Matrix{Float}
    covs_::Array{Float,3}
end

struct KalmanSmoother
    means::Matrix{Float}
    covs::Array{Float,3}
end

function kalmanfilter(lg::LinearGaussian, observations::Matrix{Float},
                      mu0::Vector{Float}, C0::Matrix{Float}
                      )::KalmanFilter
    @assert issymmetric(C0) && isposdef(C0) "Cov mat must be sym + pos def"

    nsteps = size(observations, 2)

    kf_means = zeros(lg.dimx, nsteps)
    kf_covs  = zeros(lg.dimx, lg.dimx, nsteps)

    kf_means[:,1]  = mu0
    kf_covs[:,:,1] = C0

    kf_means_ = copy(kf_means)
    kf_covs_  = copy(kf_covs)

    for k = 2:nsteps
        # intermediate (k|k-1)
        mu_ = lg.A * kf_means[:,k-1]
        C_  = lg.A * kf_covs[:,:,k-1] * lg.A' + lg.Q
        G_  = (C_ * lg.B') / (lg.R + lg.B*C_*lg.B')
        # update (k|k)
        kf_means[:,k]  = mu_ + G_*(observations[:,k]-lg.B*mu_)
        kf_covs[:,:,k] = (eye(lg.dimx) - G_*lg.B)*C_
        # storage for smoothing
        kf_means_[:,k]  = mu_
        kf_covs_[:,:,k] = C_
    end
    KalmanFilter(kf_means, kf_covs, kf_means_, kf_covs_)
end

"""
    kalmansmoother(lg, observations, kf)

Two-filter smoother in linear gaussian case (everything is tractable).
"""
function kalmansmoother(lg::LinearGaussian, observations::Matrix{Float},
                        kf::KalmanFilter)::KalmanSmoother
    ks_means = similar(kf.means)
    ks_covs  = similar(kf.covs)
    # Nsteps
    K = size(observations,2)
    # Pre-computations
    iQA, iRB = lg.Q\lg.A, lg.R\lg.B
    Pbi      = lg.B'*iRB
    Pbi_K    = copy(Pbi)
    c        = iRB'*observations[:,K]
    # Initialisation
    ks_covs[:,:,K] = inv( inv(kf.covs_[:,:,K]) + Pbi_K )
    ks_means[:,K]  = ks_covs[:,:,K]*( kf.covs_[:,:,K]\kf.means_[:,K] + c )
    # main loop
    for k = (K-1):-1:1
        # Kalman Smoother with 2 Filter Smoother update
        G   = (eye(lg.dimx) + lg.Q*Pbi)\lg.A
        Pbi = Pbi_K+iQA'*(lg.A-G)
        c   = iRB'*observations[:,k] + G'*c
        #
        ks_covs[:,:,k] = inv( inv(kf.covs_[:,:,k]) + Pbi )
        ks_means[:,k]  = ks_covs[:,:,k]*( kf.covs_[:,:,k]\kf.means_[:,k] + c )
    end
    KalmanSmoother(ks_means, ks_covs)
end


text:
if basename(pwd()) == "aoc"
    cd("2016/21")
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
module TestDay06

import AOC2021.Day06.part1, AOC2021.Day06.part2
using Test

test_input_raw = """
3,4,3,1,2"""

test_input = test_input_raw # no transformation needed her

@testset "Day 06" begin
    @testset "part 1" begin
        @test part1(test_input) == 5934
    end
    @testset "part 2" begin
        @test part2(test_input) == 26984457539
    end
end

end # module


text:
using Documenter, XScales

makedocs(;
    modules=[XScales],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/zahachtah/XScales.jl/blob/{commit}{path}#L{line}",
    sitename="XScales.jl",
    authors="Jon Norberg, Stockholm Resilience Centre, Stockholm University",
    assets=String[],
)

deploydocs(;
    repo="github.com/zahachtah/XScales.jl",
)


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
project.jl defines methods for the project(x, y) function
"""

# --- Exports

export project

# --- Method definitions

using LinearAlgebra: I

"""
    project(M::AbstractMultivector, B::AbstractBlade)::AbstractMultivector

Compute the projection of `M` onto the subspace represented by `B`.

    project(B::AbstractBlade, C::AbstractBlade;
            return_blade=true)::Union{AbstractBlade, AbstractFloat,
                                      Vector, Matrix, LinearAlgebra.I}

Compute the projection of `B` onto the subspace represented by `C`.

When `return_blade` is true, the return value is an AbstractBlade. Otherwise,
the return value is an AbstractFloat if the result is a scalar, a Vector if the
result is a vector, a Matrix if the result is a blade with 1 < grade < `dim`,
and a multiple of LinearAlgebra.I if the result is a pseudoscalar.
"""
function project end

# ------ Specializations involving an AbstractMultivector instance

# M::AbstractMultivector, B::Zero
project(M::AbstractMultivector, B::Zero) = B

# ------ Specializations involving an AbstractBlade instance

# B::AbstractBlade, C::AbstractScalar
# B::AbstractScalar, C::AbstractBlade
project(B::AbstractScalar, C::AbstractBlade; return_blade::Bool=true) =
    return_blade ? B : value(B)
project(B::AbstractBlade, C::AbstractScalar; return_blade::Bool=true) =
    return_blade ? zero(B) : 0

# B::AbstractBlade, C::Zero
project(B::AbstractBlade, C::Zero; return_blade::Bool=true) =
    return_blade ? zero(B) : 0

# B::AbstractBlade, C::Real
# B::Real, C::AbstractBlade
project(B::AbstractBlade, C::Real; return_blade::Bool=true) =
    return_blade ? zero(B) : 0
project(B::Real, C::AbstractBlade; return_blade::Bool=true) =
    return_blade ? Scalar(B) : B

# ------ Specializations involving a Blade instance

# B::Blade, C::Blade
function project(B::Blade, C::Blade; return_blade::Bool=true)
    # --- Check arguments

    assert_dim_equal(B, C)

    # --- Handle edge cases

    # When grade(B) > grade(C), the projection is zero.
    if grade(B) > grade(C)
        return return_blade ? zero(B) : 0
    end

    # Compute the projections of basis(B) onto basis(C)
    projections = Matrix{typeof(volume(B))}(undef, dim(B), grade(B))
    for i in 1:grade(B)
        projections[:, i] = project(basis(B)[:, i], C, return_blade=false)
    end

    # Incorporate volume(B) into the norm of the first projection vector
    projections[:, 1] *= volume(B)

    # Compute projection (using fact that projection is an outermorphism)
    return_blade ? Blade(projections) : projections
end

# B::Blade, C::Pseudoscalar
# B::Pseudoscalar, C::Blade
function project(B::Blade, C::Pseudoscalar; return_blade::Bool=true)
    assert_dim_equal(B, C)
    if return_blade
        return B
    else
        projections = copy(basis(B))
        projections[:, 1] *= volume(B)
        return projections
    end
end

function project(B::Pseudoscalar, C::Blade; return_blade::Bool=true)
    assert_dim_equal(B, C)
    return_blade ? zero(B) : 0
end

# B::Blade, v::Vector
# v::Vector, B::Blade
function project(B::Blade, v::Vector{<:Real}; return_blade::Bool=true)
    # Check arguments
    assert_dim_equal(v, B)

    # Compute projection
    projection = (grade(B) == 1) ?
        v * LinearAlgebra.dot(v, basis(B)) : 0

    return_blade ? Blade(projection) : projection
end

function project(v::Vector{<:Real}, B::Blade; return_blade::Bool=true)
    # Check arguments
    assert_dim_equal(v, B)

    # Compute projection
    projection = (grade(B) == 1) ?
        basis(B) * LinearAlgebra.dot(v, basis(B)) :
        basis(B) * transpose(transpose(v) * basis(B))

    return_blade ? Blade(projection) : projection
end

# ------ Specializations involving a Pseudoscalar instance

# B::Pseudoscalar, C::Pseudoscalar
function project(B::Pseudoscalar, C::Pseudoscalar; return_blade=true)
    assert_dim_equal(B, C)
    return_blade ? B : value(B) * I
end

# B::Pseudoscalar, v::AbstractScalar
project(B::Pseudoscalar, C::AbstractScalar; return_blade::Bool=true) =
    return_blade ? zero(B) : 0

# B::Pseudoscalar, C::Zero
project(B::Pseudoscalar, C::Zero; return_blade::Bool=true) =
    return_blade ? zero(B) : 0

# B::Pseudoscalar, v::Vector
# v::Vector, B::Pseudoscalar
function project(B::Pseudoscalar, v::Vector{<:Real}; return_blade::Bool=true)
    assert_dim_equal(v, B)
    return_blade ? zero(B) : 0
end

function project(v::Vector{<:Real}, B::Pseudoscalar; return_blade::Bool=true)
    assert_dim_equal(v, B)
    return_blade ? Blade(v) : v
end

# ------ Specializations involving an AbstractScalar instance

# B::AbstractScalar, C::AbstractScalar
project(B::AbstractScalar, C::AbstractScalar; return_blade::Bool=true) =
    return_blade ? B : value(B)

# B::AbstractScalar, C::Zero
project(B::AbstractScalar, C::Zero; return_blade::Bool=true) =
    return_blade ? zero(B) : 0

# B::AbstractScalar, x::Real
# x::Real, B::AbstractScalar
project(B::AbstractScalar, x::Real; return_blade::Bool=true) =
    return_blade ? B : value(B)

project(x::Real, B::AbstractScalar; return_blade::Bool=true) =
    return_blade ? Scalar{typeof(value(B))}(x) : x

# B::AbstractScalar, v::Vector
# v::Vector, B::AbstractScalar
project(B::AbstractScalar, v::Vector{<:Real}; return_blade::Bool=true) =
    return_blade ? B : value(B)

project(v::Vector{<:Real}, B::AbstractScalar; return_blade::Bool=true) =
    return_blade ? zero(B) : 0

# ------ Specializations involving a Zero instance

# B::Zero, x::Real
# x::Real, B::Zero
project(B::Zero, x::Real; return_blade::Bool=true) =
    return_blade ? zero(B) : 0

project(x::Real, B::Zero; return_blade::Bool=true) =
    return_blade ? zero(B) : 0

# B::Zero, v::Vector
# v::Vector, B::Zero
project(B::Zero, v::Vector{<:Real}) = B
project(v::Vector{<:Real}, B::Zero) = B


text:
using Test
using StatsPlots
using MCMCChains

n_iter = 500
n_name = 3
n_chain = 2

val = randn(n_iter, n_name, n_chain) .+ [1, 2, 3]'
val = hcat(val, rand(1:2, n_iter, 1, n_chain))

chn = Chains(val)

@testset "Plotting tests" begin

    # plotting singe plotting types
    ps_trace = traceplot(chn, 1)
    @test isa(ps_trace, Plots.Plot)

    ps_mean = meanplot(chn, 1)
    @test isa(ps_mean, Plots.Plot)

    ps_density = density(chn, 1)
    @test isa(ps_density, Plots.Plot)

    ps_density = density(chn, 1, append_chains=true)
    @test isa(ps_density, Plots.Plot)

    ps_autocor = autocorplot(chn, 1)
    @test isa(ps_autocor, Plots.Plot)

    #ps_contour = plot(chn, :contour)

    ps_hist = histogram(chn, 1)
    @test isa(ps_hist, Plots.Plot)

    ps_mixed = mixeddensity(chn, 1)
    @test isa(ps_mixed, Plots.Plot)

    # plotting combinations
    ps_trace_mean = plot(chn)
    @test isa(ps_trace_mean, Plots.Plot)

    ps_trace_mean = plot(chn, append_chains=true)
    @test isa(ps_trace_mean, Plots.Plot)

    savefig("demo-plot.png")

    ps_mixed_auto = plot(chn, seriestype = (:mixeddensity, :autocorplot))
    @test isa(ps_mixed_auto, Plots.Plot)

    # Test plotting using colordim keyword
    p_colordim = plot(chn, colordim = :parameter)
    @test isa(p_colordim, Plots.Plot)

    # Test if plotting a sub-set work.s
    p_subset = plot(chn, 2)
    @test isa(p_subset, Plots.Plot)

    p_subset_colordim = plot(chn, 2, colordim = :parameter)
    @test isa(p_subset_colordim, Plots.Plot)

    rm("demo-plot.png")
end


text:
function rugplot(; x::AbstractVector = Int[])
    v = VegaVisualization(width = 400, height = 20, padding = "auto")

    add_data!(v, x = x, y = zeros(Int,length(x)))

    #Get unique table name
    table = v.data[1].name

    default_scales!(v)

    #Creates axis and light gray lines
    v.axes = Array{VegaAxis}(1)
    v.axes[1] = VegaAxis(_type = "x",
                                  scale = "x",
                                  properties = Dict{Any, Any}("stroke" => Dict{Any, Any}("value" => "#000000"),
                                                              "opacity" => Dict{Any, Any}("value" => 0.08)
                                                             ),
                                  layer = "back",
                                  ticks = 5,
                                  titleOffset = 38,
                                  grid = true,
                              )

    #Add data
    v.marks = Array{VegaMark}(1)
    v.marks[1] = VegaMark(_type = "rect",
                           from = VegaMarkFrom(data = table),
                           properties = VegaMarkProperties(enter = VegaMarkPropertySet(x = VegaValueRef(scale = "x", field = "x"),
                                                                                y =  VegaValueRef(value = 0),
                                                                                width =  VegaValueRef(value = 1),
                                                                                height =  VegaValueRef(value = 14),
                                                                                fill = VegaValueRef(scale = "group", field = "group"),
                                                                                opacity = VegaValueRef(value = 0.7)
                                                                                )
                                                          )
                          )

    return v
end


text:
struct Simulation{T, TM, TI, TA, TO}
    model :: TM
    q_integ :: TI
    analyzer :: TA
    output :: TO
    dz0 :: T
    phik :: T
    phip :: T
    phimax :: T
    Icrit :: T
end

function Simulation(
    model, n2, N0;
    prefix="", dz0, dzout=dz0, phimax=pi/100, Icrit=Inf, alg=RK4(),
)
    (; grid, field, response, zu, z) = model

    # integrator:
    # use initial J=0 as initial condition to avoid creating dummy array
    q_prob = ODEIntegrators.Problem(q_func!, response.J, (model,))
    q_integ = ODEIntegrators.Integrator(q_prob, alg)

    analyzer = FieldAnalyzer(grid, z)

    output = Output(grid, field, N0, zu, z, dzout; prefix)

    phik = phi_kerr(field, n2, zu)
    phip = phi_plasma(field, N0, zu)

    T = float_type(grid.arch)
    TM = typeof(model)
    TI = typeof(q_integ)
    TA = typeof(analyzer)
    TO = typeof(output)
    return Simulation{T, TM, TI, TA, TO}(
        model, q_integ, analyzer, output, dz0, phik, phip, phimax, Icrit,
    )
end

function run!(simulation)
    (; model, q_integ, analyzer, output) = simulation
    (; dz0, phik, phip, phimax, Icrit) = simulation
    (; grid, field, response, z, zmax, nonlinearity) = model

    # calculate response to find the initial electron density:
    if nonlinearity
        calculate_response!(response, field.E, z)
    end

    analyze!(analyzer, grid, field, response, z)
    write_output(output, field, analyzer, z)

    dz = dz0
    zfirst = true

    while abs(z - zmax) > dz / 2
        if nonlinearity
            dzk = phimax / (phik * analyzer.Imax)
            dzp = phimax / (phip * analyzer.nemax)
            dz = min(dz0, dzk, dzp)
        end
        z = z + dz

        time2frequency!(field.E, field.FFT)
        if nonlinearity
            @timeit "Q step" ODEIntegrators.step!(q_integ, field.E, z, dz)
        end
        @timeit "K step" K_step!(field, model.KK, dz)
        frequency2time!(field.E, field.FFT)
        apply_guard_real_domain!(field.E, model.guard)

        @timeit "output" begin
            analyze!(analyzer, grid, field, response, z)
            write_output(output, field, analyzer, z)
        end

        # Exclude the first step from the timings:
        if zfirst
            reset_timer!(get_defaulttimer())
            zfirst = false
        end

        if analyzer.Imax >= Icrit
            @warn "Imax >= Icrit"
            break
        end
    end

    println(get_defaulttimer())
end

function K_step!(field::FieldT, KK, dz)
    @. field.E = field.E * exp(1im * KK * dz)
    return nothing
end

function K_step!(field::FieldRT, KK, dz)
    space2frequency!(field.E, field.DHT)
    @. field.E = field.E * exp(1im * KK * dz)
    frequency2space!(field.E, field.DHT)
    return nothing
end

function phi_kerr(field, n2, zu)
    (; w0, Eu) = field

    mu = permeability(w0)
    k0 = k_func(w0)
    chi3 = chi3_func(w0, n2)

    QQ0 = MU0 * mu * w0^2 / (2 * k0) * zu / Eu
    Rk0 = EPS0 * chi3 * 3 / 4 * Eu^3
    return QQ0 * abs(real(Rk0))
end

function phi_plasma(field, neu, zu)
    (; w0, Eu) = field

    mu = permeability(w0)
    k0 = k_func(w0)

    QQ0 = MU0 * mu * w0^2 / (2 * k0) * zu / Eu
    Rp0 = 1im / w0 * QE^2 / ME / (-1im * w0) * neu * Eu
    return QQ0 * abs(real(Rp0))
end


text:
module CollabUtils2MiniHeartDoor

using ..Ahorn, Maple

@mapdef Entity "CollabUtils2/MiniHeartDoor" MiniHeartDoor(x::Integer, y::Integer, width::Integer=40, height::Integer=Maple.defaultBlockHeight,
    requires::Integer=0, startHidden::Bool=false, levelSet::String="SpringCollab2020/1-Beginner", color::String="beginner", doorID::String="")

const placements = Ahorn.PlacementDict(
    "Mini Heart Door (Collab Utils 2 / READ DOCS)" => Ahorn.EntityPlacement(
        MiniHeartDoor,
        "rectangle"
    )
)

Ahorn.editingOptions(entity::MiniHeartDoor) = Dict{String, Any}(
  "color" => String["beginner", "intermediate", "advanced", "expert", "grandmaster"]
)

Ahorn.minimumSize(entity::MiniHeartDoor) = 40, 8
Ahorn.resizable(entity::MiniHeartDoor) = true, true

Ahorn.nodeLimits(entity::MiniHeartDoor) = 0, 1

const heartPadding = 4
const wallColor = (47, 187, 255, 255) ./ 255

function Ahorn.selection(entity::MiniHeartDoor, room::Maple.Room)
    nodes = get(entity.data, "nodes", ())

    if isempty(nodes)
        return Ahorn.getEntityRectangle(entity)
    else
        nx, ny = Int.(nodes[1])
        width = get(entity.data, "width", 40)

        return [Ahorn.getEntityRectangle(entity), Ahorn.Rectangle(nx - 8, ny, width + 16, 8)]
    end
end

function heartsWidth(heartSprite::Ahorn.Sprite, hearts::Integer)
    return hearts * (heartSprite.width + heartPadding) - heartPadding
end

function heartsPossible(width::Integer, edgeSprite::Ahorn.Sprite, heartSprite::Ahorn.Sprite, required::Integer)
    rowWidth = width - 2 * edgeSprite.width

    for i in 0:required
        if heartsWidth(heartSprite, i) > rowWidth
            return i - 1
        end
    end

    return required
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::MiniHeartDoor, room::Maple.Room)
    x, y = Ahorn.position(entity)
    nodes = get(entity.data, "nodes", ())

    width = Int(get(entity.data, "width", 8))

    if !isempty(nodes)
        nx, ny = Int.(nodes[1])
        dy = ny - y

        Ahorn.drawRectangle(ctx, x, ny, width, 1, (1.0, 0.0, 0.0, 1.0), (0.0, 0.0, 0.0, 0.0))
        Ahorn.drawRectangle(ctx, x, 2 * y - ny, width, 1, (1.0, 0.0, 0.0, 1.0), (0.0, 0.0, 0.0, 0.0))

        Ahorn.drawRectangle(ctx, nx - 8, ny, width + 16, 8, (1.0, 0.0, 0.0, 1.0), (0.0, 0.0, 0.0, 0.0))
    end
end

# Not completely accurate on heart positions, but good enough
function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::MiniHeartDoor, room::Maple.Room)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 40))
    height = get(entity.data, "height", Maple.defaultBlockHeight)
    hearts = Int(get(entity.data, "requires", 0))

    edgeSprite = Ahorn.getSprite("objects/heartdoor/edge", "Gameplay")
    heartSprite = Ahorn.getSprite("objects/heartdoor/icon00", "Gameplay")

    Ahorn.drawRectangle(ctx, x, y - height, width, height * 2, wallColor, (0.0, 0.0, 0.0, 0.0))

    for i in y-height:edgeSprite.height:y+height-1
        Ahorn.Cairo.save(ctx)

        Ahorn.drawImage(ctx, edgeSprite, x + width - edgeSprite.width, i)
        Ahorn.scale(ctx, -1, 1)
        Ahorn.drawImage(ctx, edgeSprite, -x - edgeSprite.width, i)

        Ahorn.Cairo.restore(ctx)
    end

    if hearts > 0
        fits = heartsPossible(width, edgeSprite, heartSprite, hearts)
        rows = ceil(Int, hearts / fits)

        for row in 1:rows
            fits = heartsPossible(width, edgeSprite, heartSprite, hearts)
            drawWidth = heartsWidth(heartSprite, fits)

            startX = x + round(Int, (width - drawWidth) / 2) + edgeSprite.width - 1
            startY = y - round(Int, rows / 2 * (heartSprite.height + heartPadding)) - heartPadding - 1

            for col in 1:fits
                drawX = (col - 1) * (heartSprite.width + heartPadding) - heartPadding
                drawY = row * (heartSprite.height + heartPadding) - heartPadding

                Ahorn.drawImage(ctx, heartSprite, startX + drawX, startY + drawY)
            end

            hearts -= fits
        end
    end
end

end


text:
import NLPModels: obj, grad, grad!, hess, objgrad

export LineModel
export obj, grad, derivative, grad!, derivative!, hess, redirect!, objgrad

"""A type to represent the restriction of a function to a direction.
Given f : R → Rⁿ, x ∈ Rⁿ and a nonzero direction d ∈ Rⁿ,

    ϕ = LineModel(nlp, x, d)

represents the function ϕ : R → R defined by

    ϕ(t) := f(x + td).
"""
mutable struct LineModel <: AbstractNLPModel
  meta :: NLPModelMeta
  counters :: Counters
  nlp :: AbstractNLPModel
  x :: Vector
  d :: Vector
  ∇ft :: Vector
  ## obj_grad :: Any ### On a dit qu'on voulait le obj_grad dans LineModel
  ## mais qu'est-ce qu'on veut? une méthode ? un compteur?
  ## la méthode obj_grad existe déjà et modifie les compteur de obj et grad
end

function LineModel(nlp :: AbstractNLPModel,
                   x   :: Vector,
                   d   :: Vector)
  meta = NLPModelMeta(1, x0=zeros(1), name="LineModel to $(nlp.meta.name))")
  g = Vector(undef, size(x)[1])
  return LineModel(meta, Counters(), nlp, x, d, g)
end

"""`redirect!(ϕ, x, d)`

Change the values of x and d of the LineModel ϕ, but retains the counters.
"""
function redirect!(ϕ :: LineModel, x :: Vector, d :: Vector)
  ϕ.x, ϕ.d = x, d
  return ϕ
end

"""`obj(f, t)` evaluates the objective of the `LineModel`

    ϕ(t) := f(x + td).
"""
function obj(f :: LineModel, t :: Number)
  NLPModels.increment!(f, :neval_obj)
  return obj(f.nlp, f.x + t * f.d)
end

"""`grad(f, t)` evaluates the first derivative of the `LineModel`

    ϕ(t) := f(x + td),

i.e.,

    ϕ'(t) = ∇f(x + td)ᵀd.
"""
function grad(f :: LineModel, t :: Number)
  NLPModels.increment!(f, :neval_grad)
  f.∇ft = grad(f.nlp, f.x + t * f.d)
  return dot(f.∇ft, f.d)
end
derivative(f :: LineModel, t :: Number) = grad(f, t)

"""`grad!(f, t, g)` evaluates the first derivative of the `LineModel`

    ϕ(t) := f(x + td),

i.e.,

    ϕ'(t) = ∇f(x + td)ᵀd.

The gradient ∇f(x + td) is stored in `g`.
"""
function grad!(f :: LineModel, t :: Number, g :: Vector)
  NLPModels.increment!(f, :neval_grad)
  f.∇ft = grad(f.nlp, f.x + t * f.d)
  g = dot(f.∇ft, f.d)
  return g
end
derivative!(f :: LineModel, t :: Number, g :: Vector) = grad!(f, t, g)

"""Evaluate the second derivative of the `LineModel`

    ϕ(t) := f(x + td),

i.e.,

    ϕ"(t) = dᵀ∇²f(x + td)d.
"""
function hess(f :: LineModel, t :: Number)
  NLPModels.increment!(f, :neval_hess)
  return dot(f.d, hprod(f.nlp, f.x + t * f.d, f.d))
end

"""Evaluate the objective function and the derivative at the same time"""
function objgrad(f :: LineModel, t :: Number)
	NLPModels.increment!(f, :neval_obj); NLPModels.increment!(f, :neval_grad)
  	f.∇ft = grad(f.nlp, f.x + t * f.d)
	obj_f = obj(f.nlp, f.x + t * f.d)
	return obj_f, dot(f.∇ft, f.d)
end

#end of module
# end


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

"""
    FoldingMethod

A method for creating folds out of a spatial object for
cross-validatory error estimation methods.
"""
abstract type FoldingMethod end

"""
    folds(object, method)

Return iterator of folds of `object` according to `method`.
"""
function folds end

# ----------------
# IMPLEMENTATIONS
# ----------------

include("folding/one.jl")
include("folding/uniform.jl")
include("folding/block.jl")
include("folding/ball.jl")


text:
"""
    Buffer(xs, ...)

`Buffer` is an array-like type which is mutable when taking gradients. You can
construct a `Buffer` with the same syntax as `similar` (e.g. `Buffer(xs, 5)`)
and then use normal indexing. Finally, use `copy` to get back a normal array.

For example:

```julia
julia> function vstack(xs)
           buf = Buffer(xs, length(xs), 5)
           for i = 1:5
             buf[:, i] = xs
           end
           return copy(buf)
         end
vstack (generic function with 1 method)

julia> vstack([1, 2, 3])
3×5 Array{Int64,2}:
 1  1  1  1  1
 2  2  2  2  2
 3  3  3  3  3

julia> gradient(x -> sum(vstack(x)), [1, 2, 3])
([5.0, 5.0, 5.0],)
```

`Buffer` is not an `AbstractArray` and can't be used for linear algebra
operations like matrix multiplication. This prevents it from being captured by
pullbacks.

`copy` is a semantic copy, but does not allocate memory. Instead the `Buffer`
is made immutable after copying.
"""
mutable struct Buffer{T,A<:AbstractArray{T}}
  data::A
  freeze::Bool
end

Buffer(xs::AbstractArray, args...) =
  Buffer(similar(xs, args...), false)

bufferfrom(xs::AbstractArray) = Buffer(xs, false)

Base.getindex(b::Buffer, i...) = b.data[i...]

function Base.setindex!(b::Buffer, v, i...)
  b.freeze && error("Buffer is frozen")
  b.data[i...] = v
end

function Base.copy(b::Buffer)
  b.freeze = true
  return b.data
end

@forward Buffer.data Base.eltype, Base.length, Base.ndims, Base.size, Base.axes, Base.eachindex, Base.stride, Base.strides

grad_mut(b::Buffer) = fill!(similar(b.data, Any), nothing)
grad_mut(b::Buffer{T}) where T<:Number = fill!(similar(b.data, float(T)), 0)

@nograd Buffer

@adjoint function getindex(b::Buffer, i...)
  b[i...], function (Δ)
    grad = grad_mut(__context__, b)
    grad[i...] = accum(grad[i...], Δ)
    return
  end
end

@adjoint! function setindex!(b::Buffer, v, i...)
  setindex!(b, v, i...), function (_)
    grad = grad_mut(__context__, b)
    v̄ = grad[i...]
    zero = eltype(grad) <: Number ? 0 : nothing
    if i isa NTuple{N,Integer} where N
      grad[i...] = zero
    else
      grad[i...] .= zero
    end
    (nothing, v̄, map(_->nothing, i)...)
  end
end

@adjoint function copy(b::Buffer)
  copy(b), function (b̄)
    grad_mut(__context__, b)[:] = b̄
    return
  end
end


text:
include("data.jl")
include("model.jl")
include("train.jl")
include("exp.jl")

function MNIST_pretrain_fn(m)
    function train_fn(model)
        ds_fn = () -> load_MNIST_ds(batch_size=50)
        ds, test_ds = ds_fn();
        x, y = next_batch!(ds) |> gpu;
        model(x)

        opt = ADAM(1e-3);
        train!(model, opt, ds, train_steps=2000, print_steps=100)
        model
    end
    maybe_train("trained/pretrain-MNIST.bson", m, train_fn)
end

function MNIST_pretrained_fn()
    m = get_Madry_model()
    MNIST_pretrain_fn(m)
end

function MNIST_AE_pretrain_fn(m)
    ds_fn = () -> load_MNIST_ds(batch_size=50)
    function train_fn(model)
        ds, test_ds = ds_fn();
        x, y = next_batch!(ds) |> gpu;
        model(x)

        opt = ADAM(1e-3);
        aetrain!(model, opt, ds, train_steps=2000, print_steps=100)
    end
    maybe_train("trained/pretrain-MNIST-ae.bson", m, train_fn)
end

function MNIST_exp_helper(expID, lr, total_steps, λ; pretrain=false)
    # This put log and saved model into MNIST subfolder
    expID = "MNIST/" * expID

    model_fn = () -> get_Madry_model()
    ds_fn = () -> load_MNIST_ds(batch_size=50)

    adv_exp_helper(expID, lr, total_steps, λ,
                   model_fn, ds_fn,
                   if pretrain MNIST_pretrain_fn else (a)->a end,
                   print_steps=20, save_steps=40,
                   test_per_steps=100, test_run_steps=20,
                   attack_fn=attack_PGD_k(40))
end

function MNIST_free_exp_helper(expID, lr, total_steps)
    expID = "MNIST-free/" * expID

    model_fn = () -> get_Madry_model()
    ds_fn = () -> load_MNIST_ds(batch_size=50)
    free_exp_helper(expID, lr, total_steps, 0.3,
                    model_fn, ds_fn,
                    (a)->a,
                    print_steps=20, save_steps=40,
                    test_per_steps=100, test_run_steps=20,
                    test_attack_fn=attack_PGD_k(40))
end

function exp_free()
    MNIST_free_exp_helper("test-$(now)", 1e-3, 1000)
end

function MNIST_dyattack_exp_helper(expID, lr, attack_fn, total_steps)
    # This put log and saved model into MNIST subfolder
    expID = "MNIST-dyattack/" * expID

    model_fn = () -> get_Madry_model()
    ds_fn = () -> load_MNIST_ds(batch_size=50)

    adv_exp_helper(expID, lr, total_steps, 0,
                   model_fn, ds_fn,
                   (a)->a,
                   print_steps=20, save_steps=40,
                   test_per_steps=100, test_run_steps=20,
                   attack_fn=attack_fn,
                   test_attack_fn=attack_PGD_k(40))
end

function exp_dyattack(schedule)
    expID = replace("$schedule", " "=>"")
    @show expID
    for m in schedule
        @show m
        # TODO lr schedule, reduce lr should increase acc further
        # TODO pretrain? Probably not.
        MNIST_dyattack_exp_helper(expID, 1e-3, attack_PGD_k(m[1]), m[2])
    end
end

function exp_dymix(schedule)
    expID = replace("$schedule", " "=>"")
    @show expID
    for m in schedule
        @show m
        steps = m[2]
        λ = m[1]
        MNIST_exp_helper(expID, 1e-3, steps, λ)
    end
end

function tmp()
    # testing recording test_attack_fn using PGD-40, this would be the default one
    # TODO tune schedule hyper-parameters
    exp_dyattack((5=>400, 10=>800, 15=>1200, 20=>1600, 30=>2000, 40=>3000))
    # exp_dyattack("test-$(now())")
    exp_dymix((5=>1000, 4=>2000, 3=>3000, 2=>4000, 1=>5000, 0=>6000))
end

function MNIST_advae_exp_helper(expID, lr, total_steps; λ=0, γ=0, β=1, pretrain=false)
    # This put log and saved model into MNIST subfolder
    expID = "MNIST-advae/" * expID

    ae_model_fn = CNN_AE
    ds_fn = () -> load_MNIST_ds(batch_size=50)

    advae_exp_helper(expID, lr, total_steps,
                     ae_model_fn, ds_fn,
                     if pretrain MNIST_AE_pretrain_fn else (a)->a end,
                     MNIST_pretrained_fn,
                     λ=λ, γ=γ, β=β,
                     print_steps=20, save_steps=40,
                     test_per_steps=100, test_run_steps=20,
                     attack_fn=attack_PGD_k(40))
end

# CAUTION the function names are the same for CIFAR

function exp_itadv(lr, total_steps)
    expID = "f0-$lr"
    MNIST_exp_helper(expID, lr, total_steps, 0)
end

function exp_pretrain(lr, total_steps)
    expID = "pretrain-$lr"
    MNIST_exp_helper(expID, lr, total_steps, 0, pretrain=true)
end

function exp_f1(lr, total_steps)
    expID = "f1-$lr"
    MNIST_exp_helper(expID, lr, total_steps, 1)
end

function exp_itadv()
    # FIXME should I use learning rate decay at the same time?
    #
    # NOTE: the steps must devide all metric steps, especially save steps,
    # otherwise it won't be saved correctly.
    exp_itadv(1e-1, 600)
    exp_itadv(5e-2, 600)
    exp_itadv(1e-2, 600)
    exp_itadv(5e-3, 600)
    exp_itadv(1e-3, 600)
    exp_itadv(5e-4, 1000)
    exp_itadv(4e-4, 1000)

    # converging from 3e-4, (HEBI: this is the border line)
    exp_itadv(3e-4, 6000)
    exp_itadv(1e-4, 8000)
    exp_itadv(5e-5, 4000)

    exp_itadv(1e-5, 3000)
end

function exp_pretrain()
    exp_pretrain(1e-2, 1000)
    # This does not converge, and I would expect nat acc to graduallly reduce 0.1
    exp_pretrain(1e-3, 1000)
    exp_pretrain(8e-4, 2000)
    # TODO and this is important because it is next to border-line. FIXME It
    # also does work
    exp_pretrain(7e-4, 2000)

    # This converges, and (HEBI: this is the border line)
    exp_pretrain(6e-4, 5000)
    # it is working here, but struggled
    exp_pretrain(5e-4, 5000)
    # TODO and I want to show how the worked one perform with a pretrained start
    exp_pretrain(3e-4, 3000)
end

function exp_f1()
    # TODO what about starting from pretrained?
    #
    # (HEBI: I hope this not to reach high accuracy)
    #
    # FIXME what if we just use a simple lr decay? The key point should be, no
    # matter how the lr change, the accuracy should still not reach the high
    # value. This might make more sense on CIFAR10 than MNIST. I'll need to stop
    # here and (HEBI: move to CIFAR NOW).
    #
    # this defnintely does not converge, this lr may not converge even for clean
    # train, I didn't try though.
    exp_f1(1e-2, 2000)
    # this converges, but end acc is not high, as expected
    exp_f1(8e-3, 3000)
    exp_f1(5e-3, 3000)
    # TODO what is the borderline of fast+acc
    exp_f1(3e-3, 5000)
    exp_f1(2e-3, 5000)
    # this should be the most promising results for this exp setting
    exp_f1(1e-3, 5000)
    exp_f1(5e-4, 1000)
end

function exp()
    # itadv train with different learning rate
    exp_itadv()
    # TODO pretrain CNN with different learning rate
    exp_pretrain()
    # TODO nat+acc 1:1 with different learning rate
    exp_f1()
    # TODO mixing data with schedule
    # TODO dynamic attacking strength
    # FIXME I will need to monitor the acc of which attack?
end

function main()
    test()
    exp()
end

##############################
## AdvAE experiments
##############################

function exp_advae_f0(lr, total_steps)
    expID = "f0-$lr"
    MNIST_advae_exp_helper(expID, lr, total_steps, λ=0)
end

function exp_advae_pretrain(lr, total_steps)
    expID = "pretrain-$lr"
    MNIST_advae_exp_helper(expID, lr, total_steps, λ=0, pretrain=true)
end

function exp_advae_f1(lr, total_steps)
    expID = "f1-$lr"
    MNIST_advae_exp_helper(expID, lr, total_steps, λ=1)
end

function exp_advae_f01(lr, total_steps)
    expID = "f01-$lr"
    MNIST_advae_exp_helper(expID, lr, total_steps, λ=0, γ=1)
end

function exp_advae_test(lr, total_steps; expID="test-$(now())")
    @show expID
    MNIST_advae_exp_helper(expID, lr, total_steps, λ=1, γ=10, β=1)
end

function test()
    exp_advae_test(2e-3, 2000, expID="test-2019-11-17T22:46:27.503")
    exp_advae_test(2e-3, 1000)

    @load "trained/pretrain-MNIST.bson" model
    # gpu(model)(x)
    exp_advae_f0(1e-3, 2000)
    exp_advae_f0(1e-4, 2000)
    exp_advae_pretrain(1e-4, 2000)
    # this works pretty good. Two things:
    #
    # 1. pretraining seems to be important here. I probably want to enable
    # pretrain by default, that does not hurt.
    #
    # 2. TODO IMPORTANT lr decay might help
    # 3. what is the performance gap (both acc and time) with itadv?
    exp_advae_pretrain(1e-3, 3000)
    # surprisingly this does not work smoothly, it struggled until a point clean
    # acc jumps
    exp_advae_f1(1e-3, 2000)
    exp_advae_f1(2e-3, 2000)

    # rec loss
    exp_advae_f01(1e-3, 2000)
end

max_stars_repo_name: lihebi/AdvAE

max_issues_repo_name: lihebi/AdvAE

max_forks_repo_name: lihebi/AdvAE


text:
"""
   interp_chunk_to_grid_gp_brute_force!( flux_out, var_out, chunk_of_spectrum, wavelengths )
Return spectra interpolated onto a grid of points using Gaussian Process interpolation.
# Arguments:
- flux_out: (results stored into this array)
- var_out: (results stored into this array)
- chunk_of_spectrum
- wavelengths: AbstractRange or AbstractArray of locations where chunk is to be interpolated to
# Returns
- flux_out

NOTE:  Using own GP code for now, since include predicting derivatives and can minimize unnecessary dependancies.  We may need to
revisit this if we want improved speed
"""

function interp_chunk_to_shifted_grid_gp_brute_force!( flux_out::AA1, var_out::AA2, chunk::AC, grid::AR, boost_factor::Real ) where { T1<:Real, AA1<:AbstractArray{T1,1}, T2<:Real, AA2<:AbstractArray{T2,1}, AC<:AbstractChunkOfSpectrum, AR<:Union{AbstractRange,AbstractArray{T2,1}} }
    @assert size(flux_out) == size(var_out)
    @assert size(flux_out) == size(grid)
    rho = 2*5000/speed_of_light_mps * mean(chunk.λ)
    sigmasq_kernel = 0.25 # Float64(mean(chunk.var))
    println(" rho = ", rho, "  σ²_kernel = ", sigmasq_kernel)
    flux_out .= GPInterpolation.predict_mean(chunk.λ./boost_factor, chunk.flux, grid, sigmasq_obs = chunk.var, kernel = GPs.matern52_sparse_kernel, rho=rho, sigmasq_cor=sigmasq_kernel ) # 	sigmasq_cor=1.0, rho=1
    # TODO: Update var_out to actually use the right GP or at least do something more sensible
    #var_out .= GPInterpolation.predict_mean(chunk.λ, chunk.var, grid, sigmasq_obs = chunk.var, kernel = GPs.matern52_sparse_kernel, rho=5000/speed_of_light_mps) # 	sigmasq_cor=1.0, rho=1
    return flux_out
end

function interp_chunk_to_shifted_grid_gp_brute_force( chunk::AC, grid::AR, boost_factor::Real ) where {  AC<:AbstractChunkOfSpectrum, T2<:Real, AR<:Union{AbstractRange,AbstractArray{T2,1}} }
    flux_out = Array{Float64,1}(undef,length(grid))
    var_out = Array{Float64,1}(undef,length(grid))
    interp_chunk_to_grid_gp_brute_force!(flux_out, var_out, chunk, grid, boost_factor=boost_factor)
    return (flux=flux_out, var=var_out)
end

function interp_chunk_to_grid_gp_brute_force!( flux_out::AA1, var_out::AA2, chunk::AC, grid::AR ) where { T1<:Real, AA1<:AbstractArray{T1,1}, T2<:Real, AA2<:AbstractArray{T2,1}, AC<:AbstractChunkOfSpectrum, AR<:Union{AbstractRange,AbstractArray{T2,1}} }
    @assert size(flux_out) == size(var_out)
    @assert size(flux_out) == size(grid)
    rho = 2*5000/speed_of_light_mps * mean(chunk.λ)
    sigmasq_kernel = 0.25 # Float64(mean(chunk.var))
    println(" rho = ", rho, "  σ²_kernel = ", sigmasq_kernel)
    flux_out .= GPInterpolation.predict_mean(chunk.λ, chunk.flux, grid, sigmasq_obs = chunk.var, kernel = GPs.matern52_sparse_kernel, rho=rho, sigmasq_cor=sigmasq_kernel ) # 	sigmasq_cor=1.0, rho=1
    # TODO: Update var_out to actually use the right GP or at least do something more sensible
    #var_out .= GPInterpolation.predict_mean(chunk.λ, chunk.var, grid, sigmasq_obs = chunk.var, kernel = GPs.matern52_sparse_kernel, rho=5000/speed_of_light_mps) # 	sigmasq_cor=1.0, rho=1
    return flux_out
end

function interp_chunk_to_grid_gp_brute_force( chunk::AC, grid::AR ) where {  AC<:AbstractChunkOfSpectrum, T2<:Real, AR<:Union{AbstractRange,AbstractArray{T2,1}} }
    flux_out = Array{Float64,1}(undef,length(grid))
    var_out = Array{Float64,1}(undef,length(grid))
    interp_chunk_to_grid_gp_brute_force!(flux_out, var_out, chunk, grid, )
    return (flux=flux_out, var=var_out)
end


text:
import RowEchelon

#=
Homework 5
Section 2.1
Problem 10
=#

# Initialize Matrices
A = [3 -6 ; -1 2]
B = [-1 1 ; 3 4]
C = [-3 -5 ; 2 1]

print("AB = ")
println(A*B)

print("AC = ")
println(A*C)

println("Determine whether B == C:")
if (B == C)
	println("True, B == C")
else
	println("False, B != C")
end


text:
using Revise # reduce recompilation time
using Plots
# using Documenter
using LinearAlgebra
using SparseArrays
using BenchmarkTools
using UnPack

push!(LOAD_PATH, "./src")
using CommonUtils
using Basis1D

using SetupDG

push!(LOAD_PATH, "./examples/EntropyStableEuler.jl/src")
using EntropyStableEuler
using EntropyStableEuler.Fluxes1D

function wavespeed_1D(rho,rhou,E)
    p = pfun_nd(rho,(rhou,),E)
    cvel = @. sqrt(γ*p/rho)
    return @. abs(rhou/rho) + cvel
end
unorm(U) = sum(map((x->x.^2),U))
function pfun_nd(rho, rhoU, E)
    rhoUnorm2 = unorm(rhoU)./rho
    return @. (γ-1)*(E - .5*rhoUnorm2)
end

function primitive_to_conservative_hardcode(rho,U,p)
    rhoU = rho.*U
    Unorm = unorm(U)
    E = @. p/(γ-1) + .5*rho*Unorm
    return (rho,rhoU,E)
end

function euler_fluxes_1D(rhoL,uL,betaL,rhologL,betalogL,
    rhoR,uR,betaR,rhologR,betalogR)

    rholog = logmean.(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean.(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = (@. .5*(rhoL+rhoR))
    uavg   = (@. .5*(uL+uR))

    unorm = (@. uL*uR)
    pa    = (@. rhoavg/(betaL+betaR))
    f4aux = (@. rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm)

    FxS1 = (@. rholog*uavg)
    FxS2 = (@. FxS1*uavg + pa)
    FxS3 = (@. f4aux*uavg)

    return (FxS1,FxS2,FxS3)
end

function euler_fluxes(rhoL,uL,betaL,rhoR,uR,betaR)
    rhologL,betalogL,rhologR,betalogR = map(x->log.(x),(rhoL,betaL,rhoR,betaR))
    return euler_fluxes_1D(rhoL,uL,betaL,rhologL,betalogL,
                           rhoR,uR,betaR,rhologR,betalogR)
end

const TOL = 1e-16
"Approximation parameters"
N = 2 # The order of approximation
K = 100
T = 3.0

# Becker viscous shocktube
const γ = 1.4
const M_0 = 3.0
const mu = 0.01
const lambda = 2/3*mu
const Pr = 3/4
const cp = γ/(γ-1)
const cv = 1/(γ-1)
const kappa = mu*cp/Pr

const v_inf = 0.2
const rho_0 = 1.0
const v_0 = 1.0
const m_0 = rho_0*v_0
const v_1 = (γ-1+2/M_0^2)/(γ+1)
const v_01 = sqrt(v_0*v_1)

const uL = v_0+v_inf
const uR = v_1+v_inf
const rhoL = m_0/v_0
const rhoR = m_0/v_1
const eL = 1/(2*γ)*((γ+1)/(γ-1)*v_01^2-v_0^2)
const eR = 1/(2*γ)*((γ+1)/(γ-1)*v_01^2-v_1^2)
const pL = (γ-1)*rhoL*eL
const pR = (γ-1)*rhoR*eR
const EL = pL/(γ-1)+0.5*rhoL*uL^2
const ER = pR/(γ-1)+0.5*rhoR*uR^2

const Bl = -1.0
const Br = 1.5
gr(size=(300,300),ylims=(0,5.0),legend=false,markerstrokewidth=1,markersize=2)
plot()

# # Sod shocktube
# const γ = 1.4
# const Bl = -0.5
# const Br = 0.5
# const rhoL = 1.0
# const rhoR = 0.125
# const pL = 1.0
# const pR = 0.1
# const uL = 0.0
# const ur = 0.0
# const xC = 0.0
# const EL = pL/(γ-1)+0.5*rhoL*uL^2
# const ER = pR/(γ-1)+0.5*rhoR*uR^2
# T = 0.2
# t0 = 0.0
# gr(size=(300,300),ylims=(0,1.2),legend=false,markerstrokewidth=1,markersize=2)
# plot()

# "Viscous parameters"
# const Re = 10000
# const Ma = 0.3
# const mu = 1/Re
# const lambda = -2/3*mu
# const Pr = .71
# const cv = 1/γ/(γ-1)/Ma^2
# const kappa = γ*cv*mu/Pr

"Mesh related variables"
VX = LinRange(Bl,Br,K+1)
EToV = transpose(reshape(sort([1:K; 2:K+1]),2,K))

"Initialize reference element"
r,_ = gauss_lobatto_quad(0,0,N)   # Reference nodes
VDM = vandermonde_1D(N,r)         # modal to nodal
Dr = grad_vandermonde_1D(N,r)/VDM # nodal differentiation
V1 = vandermonde_1D(1,r)/vandermonde_1D(1,[-1;1]) # nodal linear interpolation

Nq = N
rq,wq = gauss_lobatto_quad(0,0,Nq)
Vq = vandermonde_1D(N,rq)/VDM
M = Vq'*diagm(wq)*Vq
Mlump = zeros(size(M))
Mlump_inv = zeros(size(M))
for i = 1:Nq+1
    Mlump[i,i] = sum(M[i,:])
    Mlump_inv[i,i] = 1.0/Mlump[i,i]
end
rf = [-1.0;1.0]
nrJ = [-1.0;1.0]
Vf = vandermonde_1D(N,rf)/VDM
LIFT = Mlump_inv*Vf'

# operators
Qr = M*Dr
B = zeros(N+1,N+1)
B[1,1] = -1
B[end,end] = 1
L = Array(spdiagm(0=>-2*ones(N+1), 1=>ones(N), -1=>ones(N)))
L[1,1] = -1
L[end,end] = -1
psi = pinv(L)*-1/2*B*ones(N+1)
S0 = zeros(N+1,N+1)
for i = 1:N+1
    for j = 1:N+1
        if L[i,j] != 0
            S0[i,j] = psi[j] - psi[i]
        end
    end
end
Qr0 = S0+1/2*B

# Drop zeros
Qr = Matrix(droptol!(sparse(Qr),TOL))
Qr0 = Matrix(droptol!(sparse(Qr0),TOL))
B = Matrix(droptol!(sparse(B),TOL))
L = Matrix(droptol!(sparse(L),TOL))
S0 = Matrix(droptol!(sparse(S0),TOL))
S = Matrix(droptol!(sparse((Qr-Qr')/2),TOL))

"""High order mesh"""
x = V1*VX[transpose(EToV)]
xf = Vf*x
mapM = reshape(1:2*K,2,K)
mapP = copy(mapM)
mapP[1,2:end] .= mapM[2,1:end-1]
mapP[2,1:end-1] .= mapM[1,2:end]

# """Periodic"""
# mapP[1] = mapM[end]
# mapP[end] = mapP[1]

"""Geometric factors"""
J = repeat(transpose(diff(VX)/2),N+1,1)
nxJ = repeat([-1;1],1,K)
rxJ = 1.0

"""Geometric factors"""
J = repeat(transpose(diff(VX)/2),N+1,1)
nxJ = repeat([-1;1],1,K)
rxJ = 1.0

"""Initial condition"""
function bisection_solve_velocity(x,max_iter,tol)
    v_L = v_1
    v_R = v_0
    num_iter = 0

    L_k = kappa/m_0/cv
    f(v) = -x+2*L_k/(γ+1)*(v_0/(v_0-v_1)*log((v_0-v)/(v_0-v_01))-v_1/(v_0-v_1)*log((v-v_1)/(v_01-v_1)))

    v_new = (v_L+v_R)/2
    while num_iter < max_iter
        v_new = (v_L+v_R)/2

        if abs(f(v_new)) < tol
            return v_new
        elseif sign(f(v_L)) == sign(f(v_new))
            v_L = v_new
        else
            v_R = v_new
        end
        num_iter += 1
    end

    return v_new
end

const max_iter = 100
const tol = 1e-14

function exact_sol_viscous_shocktube(x,t)
    u   = bisection_solve_velocity(x-v_inf*t,max_iter,tol)
    rho = m_0/u
    e   = 1/(2*γ)*((γ+1)/(γ-1)*v_01^2-u^2)
    return rho, rho*(v_inf+u), rho*(e+1/2*(v_inf+u)^2)
end

U = exact_sol_viscous_shocktube.(x,0.0)
U = ([x[1] for x in U], [x[2] for x in U], [x[3] for x in U])

function zhang_wavespd(rhoL,rhouL,EL,sigma2L,sigma3L,pL,rhoR,rhouR,ER,sigma2R,sigma3R,pR)
    uL   = rhouL/rhoL
    eL   = (EL - .5*rhoL*uL^2)/rhoL
    tauL = sigma2L
    qL   = uL*tauL-sigma3L
    uR   = rhouR/rhoR
    eR   = (ER - .5*rhoR*uR^2)/rhoR
    tauR = sigma2R
    qR   = uR*tauR-sigma3R
    wavespdL = abs(uL)+1/(2*rhoL^2*eL)*(sqrt(rhoL^2*qL^2+2*rhoL^2*eL*abs(tauL-pL)^2)+rhoL*abs(qL))
    wavespdR = abs(uR)+1/(2*rhoR^2*eR)*(sqrt(rhoR^2*qR^2+2*rhoR^2*eR*abs(tauR-pR)^2)+rhoR*abs(qR))
    return max(wavespdL,wavespdR)
end

function zhang_wavespd(rhoL,rhouL,EL,sigma2L,sigma3L,pL)
    uL   = rhouL/rhoL
    eL   = (EL - .5*rhoL*uL^2)/rhoL
    tauL = sigma2L
    qL   = uL*tauL-sigma3L
    wavespdL = abs(uL)+1/(2*rhoL^2*eL)*(sqrt(rhoL^2*qL^2+2*rhoL^2*eL*abs(tauL-pL)^2)+rhoL*abs(qL))
    return wavespdL
end

function flux_lowIDP(U_i,U_j,f_i,f_j,c_ij,wavespd)
    return c_ij*(f_i+f_j)-abs(c_ij)*wavespd*(U_j-U_i)
end

function flux_viscous(sigma_i,sigma_j,c_ij)
    return c_ij*(sigma_i+sigma_j)
end

function rhs_low_viscous(U,K,N)
    J = (Br-Bl)/K/2 # assume uniform interval
    Nc = 3
    rhsU  = [zeros(N+1,K) for i = 1:Nc]
    theta = [zeros(N+1,K) for i = 1:Nc]
    sigma = [zeros(N+1,K) for i = 1:Nc]

    p = pfun_nd.(U[1],U[2],U[3])

    VU = v_ufun(U...)
    VUf = (x->Vf*x).(VU)
    VUP = (x->x[mapP]).(VUf)
    VL = v_ufun(rhoL,rhoL*uL,EL)
    VUP[1][1] = VL[1]
    VUP[2][1] = VL[2]
    VUP[3][1] = VL[3]
    VUP[1][end] = VUf[1][end]
    VUP[2][end] = VUf[2][end]
    VUP[3][end] = VUf[3][end]

    surfx(uP,uf) = LIFT*(@. .5*(uP-uf)*nxJ)
    # Strong form, dv/dx
    VUx = (x->Dr*x).(VU)
    VUx = VUx .+ surfx.(VUP,VUf)
    VUx = VUx./J

    # σ = K dv/dx
    for k = 1:K
        for i = 1:N+1
            Kx = zeros(3,3)
            v1 = VU[1][i,k]
            v2 = VU[2][i,k]
            v4 = VU[3][i,k]
            Kx[2,2] = -(2*mu-lambda)/v4
            Kx[2,3] = (2*mu-lambda)*v2/v4^2
            Kx[3,2] = Kx[2,3]
            Kx[3,3] = -(2*mu-lambda)*v2^2/v4^3+kappa/cv/v4^2

            sigma[2][i,k] += Kx[2,2]*VUx[2][i,k] + Kx[2,3]*VUx[3][i,k]
            sigma[3][i,k] += Kx[3,2]*VUx[2][i,k] + Kx[3,3]*VUx[3][i,k]
        end
    end

    sxf = (x->Vf*x).(sigma)
    sxP = (x->x[mapP]).(sxf)
    sxP[1][1] = sxf[1][1]
    sxP[2][1] = sxf[2][1]
    sxP[3][1] = sxf[3][1]
    sxP[1][end] = sxf[1][end]
    sxP[2][end] = sxf[2][end]
    sxP[3][end] = sxf[3][end]
    # strong form, dσ/dx
    penalization(uP,uf) = LIFT*(@. -.5*(uP-uf))
    sigmax = (x->Dr*x).(sigma)
    sigmax = sigmax .+ surfx.(sxP,sxf) #.+ penalization.(VUP,VUf)
    sigmax = sigmax./J

    return sigmax,sigma
end

function rhs_IDP(U,K,N,Mlump_inv,S0,S,wq,dt)
    p = pfun_nd.(U[1],U[2],U[3])
    flux = zero.(U)
    @. flux[1] = U[2]
    @. flux[2] = U[2]^2/U[1]+p
    @. flux[3] = U[3]*U[2]/U[1]+p*U[2]/U[1]

    Ub = zero.(U)
    @. Ub[1] = U[1]
    @. Ub[2] = U[2]/U[1]
    @. Ub[3] = U[1]/(2*p)

    J = (Br-Bl)/K/2 # assume uniform interval

    # Low order and high order algebraic fluxes
    F_low       = [zeros(N+1,N+1),zeros(N+1,N+1),zeros(N+1,N+1)]
    F_low_P     = [zeros(2),zeros(2),zeros(2)] # 1: left boundary, 2: right boundary
    L           = zeros(N+1,N+1) # Array of limiting params
    rhsU        = [zeros(N+1,K),zeros(N+1,K),zeros(N+1,K)]

    # TODO: redundant!
    _,sigma = rhs_low_viscous(U,K,N)
    wavespd_arr = zeros(N+1,K)
    beta_arr = zeros(N+1,K)
    for k = 1:K
        for i = 1:N+1
            wavespd_arr[i,k] = wavespeed_1D(U[1][i,k],U[2][i,k],U[3][i,k])
            beta_arr[i,k] = zhang_wavespd(U[1][i,k],U[2][i,k],U[3][i,k],sigma[2][i,k],sigma[3][i,k],p[i,k])
        end
    end

    L_plot = zeros(N+1,K)

    for k = 1:K
        L = zeros(N+1,N+1)
        for i = 1:N+1
            for j = 1:N+1
                if i != j
                    wavespd = max(wavespd_arr[i,k],wavespd_arr[j,k],beta_arr[i,k],beta_arr[j,k])
                    for c = 1:3
                        F_low[c][i,j]  = -flux_viscous(sigma[c][i,k],sigma[c][j,k],S[i,j])+flux_lowIDP(U[c][i,k],U[c][j,k],flux[c][i,k],flux[c][j,k],S0[i,j],wavespd)
                    end
                end
            end
        end

        U_left   = (k == 1) ? [rhoL; rhoL*uL; EL]                 : [U[1][end,k-1]; U[2][end,k-1]; U[3][end,k-1]]
        f_left   = (k == 1) ? [rhoL*uL; rhoL*uL^2+pL; uL*(pL+EL)] : [flux[1][end,k-1]; flux[2][end,k-1]; flux[3][end,k-1]]
        U_right  = (k == K) ? [rhoR; rhoR*uR; ER]                 : [U[1][1,k+1]; U[2][1,k+1]; U[3][1,k+1]]
        f_right  = (k == K) ? [rhoR*uR; rhoR*uR^2+pR; uR*(pR+ER)] : [flux[1][1,k+1]; flux[2][1,k+1]; flux[3][1,k+1]]
        p_left   = (k == 1) ? pL : p[end,k-1]
        p_right  = (k == K) ? pR : p[1,k+1]
        # Assume velocity doesn't change at boundary, isothermal
        TL = pL/rhoL/(γ-1)/cv
        TR = pR/rhoR/(γ-1)/cv
        sigma_left  = (k == 1) ? [sigma[1][1,1];sigma[2][1,1];sigma[3][1,1]]       : [sigma[1][end,k-1]; sigma[2][end,k-1]; sigma[3][end,k-1]]
        sigma_right = (k == K) ? [sigma[1][end,k];sigma[2][end,k];sigma[3][end,k]] : [sigma[1][1,k+1]; sigma[2][1,k+1]; sigma[3][1,k+1]]
        wavespd_l = max(wavespd_arr[1,k],wavespeed_1D(U_left[1],U_left[2],U_left[3]),beta_arr[1,k],zhang_wavespd(U_left[1],U_left[2],U_left[3],sigma_left[2],sigma_left[3],p_left))
        wavespd_r = max(wavespd_arr[end,k],wavespeed_1D(U_right[1],U_right[2],U_right[3]),beta_arr[end,k],zhang_wavespd(U_right[1],U_right[2],U_right[3],sigma_right[2],sigma_right[3],p_right))

        for c = 1:3
            F_low_P[c][1] = -flux_viscous(sigma[c][1,k],sigma_left[c],-1/2)+flux_lowIDP(U[c][1,k],U_left[c],flux[c][1,k],f_left[c],-0.5,wavespd_l)
            F_low_P[c][2] = -flux_viscous(sigma[c][end,k],sigma_right[c],1/2)+flux_lowIDP(U[c][end,k],U_right[c],flux[c][end,k],f_right[c],0.5,wavespd_r)
        end

        for c = 1:3
            # With limiting
            rhsU[c][:,k] = -sum(F_low[c],dims=2)

            rhsU[c][1,k] += -F_low_P[c][1]
            rhsU[c][N+1,k] += -F_low_P[c][2]

            rhsU[c][:,k] .= 1/J*Mlump_inv*rhsU[c][:,k]
        end
    end

    return rhsU
end

# Time stepping
"Time integration"
t = 0.0
U = collect(U)
resU = [zeros(size(x)),zeros(size(x)),zeros(size(x))]
resW = [zeros(size(x)),zeros(size(x)),zeros(size(x))]
resZ = [zeros(size(x)),zeros(size(x)),zeros(size(x))]

Vp = vandermonde_1D(N,LinRange(-1,1,10))/VDM

anim = Animation()
const ptL = Bl+(Br-Bl)/K/(N+1)/2
const ptR = Br-(Br-Bl)/K/(N+1)/2
const hplot = (Br-Bl)/K/(N+1)
i = 1

while t < T

    # SSPRK(3,3)
    # rhsU,dt,_ = rhs_IDP(U,K,N,Mlump_inv,S0,S,wq)
    # dt = min(dt,T-t)
    dt = min(1e-4,T-t)
    rhsU = rhs_IDP(U,K,N,Mlump_inv,S0,S,wq,dt)
    @. resW = U + dt*rhsU
    rhsU = rhs_IDP(resW,K,N,Mlump_inv,S0,S,wq,dt)
    @. resZ = resW+dt*rhsU
    @. resW = 3/4*U+1/4*resZ
    rhsU = rhs_IDP(resW,K,N,Mlump_inv,S0,S,wq,dt)
    @. resZ = resW+dt*rhsU
    @. U = 1/3*U+2/3*resZ

    global t = t + dt
    global i = i + 1
    println("Current time $t with time step size $dt, and final time $T")
    # if mod(i,1000) == 1
    #     plot(x,U[1])
    # end
end

plot(x,U[1])
#gif(anim,"~/Desktop/tmp.gif",fps=15)

exact_U = @. exact_sol_viscous_shocktube.(x,T)
exact_rho = [x[1] for x in exact_U]
exact_u = [x[2] for x in exact_U]
exact_u = exact_u./exact_rho
exact_E = [x[3] for x in exact_U]

rho = U[1]
u = U[2]./U[1]
E = U[3]
p = pfun_nd.(U[1],U[2],U[3])
J = (Br-Bl)/K/2

Linferr = maximum(abs.(exact_rho-rho))/maximum(abs.(exact_rho)) +
          maximum(abs.(exact_u-u))/maximum(abs.(exact_u)) +
          maximum(abs.(exact_E-E))/maximum(abs.(exact_E))

L1err = sum(J*wq.*abs.(exact_rho-rho))/sum(J*wq.*abs.(rho)) +
        sum(J*wq.*abs.(exact_u-u))/sum(J*wq.*abs.(u)) +
        sum(J*wq.*abs.(exact_E-E))/sum(J*wq.*abs.(E))

println("N = $N, K = $K")
println("L1 error is $L1err")
println("Linf error is $Linferr")

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
function options_ASB07()
    𝑂spec = Vector{OptionSpec}()

    push!(𝑂spec, OptionSpec(:δ, 1e-2, domain=Float64, aliases=[:sampling_time],
                            domain_check=(v  ->  v > 0.), info="time step"))

    push!(𝑂spec, OptionSpec(:max_order, 10, domain=Int,
                            info="maximum allowed order of zonotopes"))

    push!(𝑂spec, OptionSpec(:order_discretization, 2, domain=Int,
                            info="order of Taylor approximation in " *
                                 "discretization"))

    push!(𝑂spec, OptionSpec(:set_operations_discretization, "zonotope",
                            domain=String,
                            domain_check=(v  ->  v ∈ ["zonotope", "lazy"]),
                            info="type of set operations applied during " *
                            "discretization"))

    return 𝑂spec
end


text:
"""
    function integratorEuler_(p, platform)

Adapter for the Euler integration step method.

```math
x(t+Δt) = x(t) + a*Δt +b*ΔW
```

The algorithm can be reduced to three kernel invocations performing the following actions:
 - Clean the interaction parameters array
 - Compute interaction parameters at original position v₀ (if there are interaction parameters)
 - Compute final position v₁
"""
function addIntegratorEuler_!(p::Program_, platform::String)

    if "Equation" in keys(p.agent.declaredUpdates)

        code = p.agent.declaredUpdates["Equation"]

        #Create interaction parameter kernel if there is any interaction parameter updated
        if "UpdateInteraction" in keys(p.agent.declaredUpdates)

            k1 = loop_[p.neighbors](p,p.agent.declaredUpdates["UpdateInteraction"],platform)
            k1 = vectorize_(p.agent,k1,p,interaction=true)
            k1 = wrapInFunction_(:interactionCompute_!,k1)
            push!(p.declareF.args,k1)

        end

        #Create integration step function
        for (i,j) in enumerate(p.agent.declaredSymbols["Local"])
            code = postwalk(x -> @capture(x,g_(s_)=v__) && g == DIFFSYMBOL && s == j ? :($j = $j + $(v...)) : x, code)
            code = postwalk(x -> @capture(x,dW) ? :(Normal(0.,sqrt(dt))) : x, code)
        end
        code = vectorize_(p.agent,code,p)
        f = simpleFirstLoopWrapInFunction_(platform,:integrationStep1_!,code)
        push!(p.declareF.args,f)

        #Create wrapped integration step function
        if "UpdateInteraction" in keys(p.agent.declaredUpdates)
            addInteraction = [:(@platformAdapt cleanInteraction_!(ARGS_);@platformAdapt interactionCompute_!(ARGS_))]
            push!(p.execInit.args, :(@platformAdapt cleanInteraction_!(ARGS_); @platformAdapt interactionCompute_!(ARGS_)))
        else
            addInteraction = []
        end
        push!(p.declareF.args,
            :(begin
                function integrationStep_!(ARGS_)
                    $(addInteraction...)
                    @platformAdapt integrationStep1_!(ARGS_)
                    #println(localVCopy[1,:])

                    return
                end
            end)
        )

        push!(p.execInloop.args,
            :(integrationStep_!(ARGS_))
        )

    end

    return
end