text:
@with_kw mutable struct Reactor <: AbstractReactor

  T_bar::AbstractSymbol

  n_bar::AbstractSymbol = symbols(:n_bar)
  I_P::AbstractSymbol = symbols(:I_P)

  R_0::AbstractSymbol = symbols(:R_0)
  B_0::AbstractSymbol = symbols(:B_0)

  sigma_v::AbstractSymbol = symbols(:sigma_v)

  mode_scaling::Dict = H_mode_scaling

  deck::AbstractKey = nothing

  is_pulsed::Bool = true
  is_symbolic::Bool = false
  is_consistent::Bool = false

  is_solved::Bool = false
  is_valid::Bool = true
  is_good::Bool = true

  H::AbstractSymbol = 1.1
  Q::AbstractSymbol = 39.86

  epsilon::AbstractSymbol = 0.3226
  kappa_95::AbstractSymbol = 1.590
  delta_95::AbstractSymbol = 0.333

  nu_n::AbstractSymbol = 0.27
  nu_T::AbstractSymbol = 1.094

  Z_eff::AbstractSymbol = 2.584
  f_D::AbstractSymbol = 0.7753

  A::AbstractSymbol = 2.735

  l_i::AbstractSymbol = 1.155

  rho_m::AbstractSymbol = symbols(:rho_m)
  gamma::AbstractSymbol = symbols(:gamma)

  N_G::AbstractSymbol = 1.2

  eta_T::AbstractSymbol = 0.3531
  eta_RF::AbstractSymbol = 0.4

  tau_FT::AbstractSymbol = 7200.0
  C_ejima::AbstractSymbol = 0.3

  pi::AbstractSymbol = AbstractFloat(Base.pi)

  speed_of_light::AbstractSymbol = 3e8

  B_CS::AbstractSymbol = 12.77

  k::AbstractSymbol = 0.333

  sigma_CS::AbstractSymbol = 660.0
  sigma_TF::AbstractSymbol = 660.0

  J_CS::AbstractSymbol = 50.0
  J_TF::AbstractSymbol = 461.1

  max_beta_N::AbstractSymbol = 0.02591
  max_q_95::AbstractSymbol = 2.5
  max_q_DV::AbstractSymbol = 10.0
  max_P_E::AbstractSymbol = 1500.0
  max_P_W::AbstractSymbol = 8.0

  eta_CD::AbstractSymbol = 0.2721

  wave_theta::AbstractSymbol = 135.0
  eta_LH::AbstractSymbol = 0.6

  tau_E::AbstractCalculated = nothing
  p_bar::AbstractCalculated = nothing
  P_F::AbstractCalculated = nothing

  beta_N::AbstractCalculated = nothing
  q_95::AbstractCalculated = nothing
  q_DV::AbstractCalculated = nothing
  P_E::AbstractCalculated = nothing
  P_W::AbstractCalculated = nothing

  norm_beta_N::AbstractCalculated = nothing
  norm_q_95::AbstractCalculated = nothing
  norm_q_DV::AbstractCalculated = nothing
  norm_P_E::AbstractCalculated = nothing
  norm_P_W::AbstractCalculated = nothing

  f_BS::AbstractCalculated = nothing
  f_CD::AbstractCalculated = nothing
  f_IN::AbstractCalculated = nothing

  resistance::AbstractCalculated = nothing
  voltage::AbstractCalculated = nothing
  inductance::AbstractCalculated = nothing

  volume::AbstractCalculated = nothing

  W_M::AbstractCalculated = nothing
  cost::AbstractCalculated = nothing

  a::AbstractCalculated = nothing
  b::AbstractCalculated = nothing
  c::AbstractCalculated = nothing
  d::AbstractCalculated = nothing

  h_CS::AbstractCalculated = nothing
  R_CS::AbstractCalculated = nothing

  constraint::AbstractKey = nothing

  skipped_limits::Vector{Symbol} = [:pcap]
  ignored_limits::Vector{Symbol} = [:heat]

  branch_id::Int = 1

end

function _Reactor!(cur_reactor::AbstractReactor, cur_kwargs::Dict)
  @assert !( cur_reactor.is_pulsed && cur_reactor.is_consistent )

  for (cur_key, cur_value) in cur_kwargs
    cur_field = getfield(cur_reactor, cur_key)

    value_type = eltype(cur_value)
    field_type = eltype(cur_field)

    if value_type <: Int && !(field_type <: Int)
      cur_value = float(cur_value)
    end

    setfield!(cur_reactor, cur_key, cur_value)
  end

  if !cur_reactor.is_pulsed
    cur_reactor.tau_FT = 1.6e9
  end

  if cur_reactor.constraint == nothing
    if cur_reactor.is_pulsed
      cur_reactor.constraint = "kink"
    else
      cur_reactor.constraint = "beta"
    end
  end

  is_bad_reactor = false

  is_bad_reactor |= ( isa(cur_reactor.f_D, AbstractFloat) && ( cur_reactor.f_D <= 0 ) )

  is_bad_reactor |= ( isa(cur_reactor.f_D, AbstractFloat) && ( cur_reactor.f_D > 1 ) )

  is_bad_reactor |= ( isa(cur_reactor.Z_eff, AbstractFloat) && ( cur_reactor.Z_eff < 1 ) )

  is_bad_reactor |= ( isa(cur_reactor.N_G, AbstractFloat) && ( cur_reactor.N_G <= 0 ) )

  is_bad_reactor |= ( isa(cur_reactor.epsilon, AbstractFloat) && ( cur_reactor.epsilon <= 0 ) )

  is_bad_reactor |= ( isa(cur_reactor.kappa_95, AbstractFloat) && ( cur_reactor.kappa_95 < 1 ) )

  is_bad_reactor |= ( isa(cur_reactor.rho_m, AbstractFloat) && ( cur_reactor.rho_m > 1 ) )

  is_bad_reactor |= ( isa(cur_reactor.rho_m, AbstractFloat) && ( cur_reactor.rho_m < 0 ) )

  is_bad_reactor |= ( isa(cur_reactor.max_beta_N, AbstractFloat) && ( cur_reactor.max_beta_N <= 0 ) )

  is_bad_reactor |= ( isa(cur_reactor.max_P_W, AbstractFloat) && ( cur_reactor.max_P_W <= 0 ) )

  is_bad_reactor |= ( isa(cur_reactor.max_q_95, AbstractFloat) && ( cur_reactor.max_q_95 <= 0 ) )

  if is_bad_reactor
    cur_reactor.is_good = false
    return
  end

  isa(cur_reactor.kappa_95, AbstractFloat) || return

  @assert (
    !isa(cur_reactor.l_i, SymEngine.Basic) ||
    !isa(cur_reactor.gamma, SymEngine.Basic) ||
    !isa(cur_reactor.rho_m, SymEngine.Basic)
  )

  @assert (
    ( isa(cur_reactor.l_i, SymEngine.Basic) && isa(cur_reactor.gamma, SymEngine.Basic) ) ||
    ( isa(cur_reactor.rho_m, SymEngine.Basic) && isa(cur_reactor.l_i, SymEngine.Basic) ) ||
    ( isa(cur_reactor.gamma, SymEngine.Basic) && isa(cur_reactor.rho_m, SymEngine.Basic) ) ||
    !( isa(cur_reactor.l_i, SymEngine.Basic) || isa(cur_reactor.gamma, SymEngine.Basic) || isa(cur_reactor.rho_m, SymEngine.Basic) )
  )

  if isa(cur_reactor.l_i, SymEngine.Basic)
    if isa(cur_reactor.rho_m, AbstractFloat)
      isa(cur_reactor.gamma, SymEngine.Basic) &&
        ( cur_reactor.gamma = 1 / ( 1 - cur_reactor.rho_m ^ 2 ) )
    end

    cur_reactor.l_i = int_b_p(cur_reactor.gamma)
    cur_reactor.l_i /= 1 + cur_reactor.kappa_95 ^ 2
    cur_reactor.l_i *= 4 * cur_reactor.kappa_95
  end

  if isa(cur_reactor.gamma, SymEngine.Basic)
    cur_lhs = cur_reactor.l_i
    cur_lhs *= 1 + cur_reactor.kappa_95 ^ 2
    cur_lhs /= 4 * cur_reactor.kappa_95

    cur_gamma = NaN
    for cur_attempt in 1:100
      try
        cur_gamma = fzero(
          tmp_gamma -> int_b_p(tmp_gamma) - cur_lhs,
          2 * rand() - 1
        )
      catch
        continue
      end

      iszero(cur_gamma) || break
    end

    cur_reactor.gamma = cur_gamma
  end

  if isa(cur_reactor.rho_m, SymEngine.Basic)
    cur_gamma = cur_reactor.gamma

    if cur_gamma < 1
      cur_reactor.rho_m = 0.0
    else
      cur_reactor.rho_m = sqrt( cur_gamma - 1 )
      cur_reactor.rho_m /= sqrt( cur_gamma )
    end
  end

  is_bad_reactor |= !isfinite(cur_reactor.rho_m)

  is_bad_reactor |= !isfinite(cur_reactor.gamma)

  is_bad_reactor |= !isfinite(cur_reactor.l_i)

  is_bad_reactor |= ( cur_reactor.rho_m > 1 )

  is_bad_reactor |= ( cur_reactor.rho_m < 0 )

  is_bad_reactor && ( cur_reactor.is_good = false )

end

function Reactor(cur_temp::Symbol; cur_kwargs...)
  cur_dict = merge!(Dict(), Dict(cur_kwargs))

  cur_reactor = Reactor(SymEngine.symbols(cur_temp), cur_dict)

  cur_reactor
end

function Reactor(cur_temp::AbstractSymbol; cur_kwargs...)
  cur_dict = merge!(Dict(), Dict(cur_kwargs))

  cur_reactor = Reactor(cur_temp, cur_dict)

  cur_reactor
end

function Reactor(cur_temp::AbstractSymbol, tmp_dict::Dict)
  cur_dict = deepcopy(tmp_dict)

  if haskey(cur_dict, :deck) && cur_dict[:deck] != nothing
    cur_deck_symbol = Symbol("$(cur_dict[:deck])_deck")
    cur_deck_func = getfield(Fussy, cur_deck_symbol)

    cur_reactor = cur_deck_func(cur_temp)

    cur_redundant_keys = [:rho_m, :gamma, :l_i]

    for cur_key in cur_redundant_keys
      haskey(cur_dict, cur_key) || continue
      for cur_sub_key in cur_redundant_keys
        ( cur_key == cur_sub_key ) && continue
        cur_dict[cur_sub_key] = symbols(cur_sub_key)
      end
      break
    end
  else
    cur_reactor = Reactor(T_bar=cur_temp)
  end

  _Reactor!(cur_reactor, cur_dict)

  cur_reactor
end

function _SymbolizeReactor!(cur_reactor::AbstractReactor)
  for cur_field_name in fieldnames(cur_reactor)
    ( cur_field_name == :T_bar ) && continue

    cur_field = getfield(cur_reactor, cur_field_name)
    isa(cur_field, AbstractFloat) || continue

    cur_field = symbols(cur_field_name)
    setfield!(cur_reactor, cur_field_name, cur_field)
  end
end

function _NanizeReactor!(cur_reactor::AbstractReactor)
  for cur_field_name in fieldnames(cur_reactor)
    ( cur_field_name == :T_bar ) && continue
    ( cur_field_name == :constraint ) && continue
    ( cur_field_name == :deck ) && continue

    cur_field = getfield(cur_reactor, cur_field_name)

    is_nan_field = isa(cur_field, AbstractFloat)
    is_nan_field |= isa(cur_field, SymEngine.Basic)
    is_nan_field |= cur_field == nothing

    is_nan_field || continue

    setfield!(cur_reactor, cur_field_name, NaN)
  end
end

function NanReactor(cur_temp::AbstractSymbol, cur_dict::Dict)
  cur_reactor = Reactor(T_bar = cur_temp)

  _NanizeReactor!(cur_reactor)

  cur_reactor
end

function NanReactor(cur_temp::AbstractSymbol=symbols(:T_bar); cur_kwargs...)
  cur_dict = merge!(Dict(), Dict(cur_kwargs))

  cur_reactor = NanReactor(cur_temp, cur_dict)

  cur_reactor
end

function SymbolicReactor(cur_temp::AbstractSymbol=symbols(:T_bar); cur_kwargs...)
  cur_dict = merge!(Dict(), Dict(cur_kwargs))

  cur_scaling = haskey(cur_dict, :mode_scaling) ?
    cur_dict[:mode_scaling] : symbol_scaling

  delete!(cur_dict, :mode_scaling)

  cur_reactor = Reactor(
    T_bar = cur_temp,
    is_symbolic = true,
    mode_scaling = cur_scaling
  )

  _SymbolizeReactor!(cur_reactor)
  _Reactor!(cur_reactor, cur_dict)

  cur_reactor
end

function BaseReactor(cur_dict::Dict)
  @assert haskey(cur_dict, :T_bar)

  cur_T_bar = cur_dict[:T_bar]

  cur_reactor = Reactor(T_bar = cur_T_bar)

  _SymbolizeReactor!(cur_reactor)

  cur_reactor.pi = AbstractFloat(Base.pi)

  cur_reactor.speed_of_light = 3e8

  _Reactor!(cur_reactor, cur_dict)

  cur_reactor
end

function BaseReactor(cur_temp::AbstractSymbol=symbols(:T_bar); cur_kwargs...)
  cur_dict = merge!(Dict(), Dict(cur_kwargs))

  BaseReactor(cur_dict)
end

Reactor(cur_temp::Number; cur_kwargs...) =
  Reactor(real(float(cur_temp)); cur_kwargs...)

SymbolicReactor(cur_temp::Number; cur_kwargs...) =
  SymbolicReactor(real(float(cur_temp)); cur_kwargs...)

BaseReactor(cur_temp::Number; cur_kwargs...) =
  BaseReactor(real(float(cur_temp)); cur_kwargs...)

NanReactor(cur_temp::Number; cur_kwargs...) =
  NanReactor(real(float(cur_temp)); cur_kwargs...)


text:
import JSON

# export readmaster, readcompounds, readids

function readexpansion(f::String)
    nexp = JSON.parsefile(f)
    x = Vector{Int}(nexp["x"])
    t = Vector{Int}(nexp["t"])
    compounds = IDs(nexp["compounds"])
    reactions = IDs(nexp["reactions"])
    X = Vector{Vector{Int}}(nexp["X"])
    Y = Vector{Vector{Int}}(nexp["Y"])
    Expansion(x,t,compounds,reactions,X,Y)
end

## The point of this is to verify the structure of the
## master file or object.
"""
    readmaster(master)

Return a `Dict{Any,Reaction}` of reactions coerced into Reaction struct.
When applied to a string, look for a file to read in.
"""
function readmaster end

function readmaster(f::String)
    master = JSON.parsefile(f)
    readmaster(master)
end

function readmaster(master::Dict)
    reactions = Reactions()
    for (k,v) in master["reactions"]
        reactions[k] = Reaction(id=k,left=v["left"],right=v["right"],metadata=v["metadata"])
    end
    return reactions
end

"""
    readcompounds(dir)

Returns a `Dict{Symbol,Compound}` of compounds coerced into Compound struct.
"""
function readcompounds(dir::String)
    compounds = Compounds()
    for f in readdir(dir)
        if endswith(f,".json")
            compound = JSON.parsefile(joinpath(dir,f))[1]
            compound = symbolizekeys(compound)
            compounds[compound[:entry_id]] = Compound(;compound...)
        end
    end
    return compounds
end

"""
    readids(f)

Return a `Vector{String}` from a vector of reactions or compounds.
When applied to a string, look for a file to read in.
"""
function readids end

function readids(f::String)
    ids = JSON.parsefile(f)
    readids(ids)
end

function readids(ids::Vector{<:Any})
    return convert(IDs,ids)
end

"""
    readkeyedids(f)

Return a `Dict{String,IDs}` from a dict-like with values of reactions or compounds.
When applied to a string, look for a file to read in.
"""
function readkeyedids(f::String)
    ids = JSON.parsefile(f)
    readkeyedids(ids)
end

function readkeyedids(ids::Dict)
    return convert(Dict{String,IDs},ids)
end

"""
    symbolizekeys(d)

Converts the keys of a dict to the type `Symbol`.
"""
function symbolizekeys(d::Dict)
    newdict = Dict{Symbol,Any}()
    for (k,v) in d
        newdict[Symbol(k)] = v
    end
    return newdict
end

ext: jl

max_stars_repo_name: hbsmith/BioXP

max_issues_repo_name: hbsmith/BioXP

max_forks_repo_name: hbsmith/BioXP


text:
# Utilities for working with LLVM source code

## cloning

export clone_into!, clone

type_mapper_callback(typ, type_mapper) =
    Base.unsafe_convert(API.LLVMTypeRef, type_mapper[](LLVMType(typ)))
materializer_callback(val, materializer) =
    Base.unsafe_convert(API.LLVMValueRef, materializer[](Value(val)))

function clone_into!(new::Function, old::Function;
                     value_map::Dict{Value,Value}=Dict{Value,Value}(),
                     changes=API.LLVMCloneFunctionChangeTypeLocalChangesOnly,
                     suffix::String="", type_mapper=nothing, materializer=nothing)
    value_map_array = Value[]
    for (src, dest) in value_map
        push!(value_map_array, src)
        push!(value_map_array, dest)
    end
    if type_mapper === nothing
        type_mapper_ptr = C_NULL
        type_mapper_data = C_NULL
    else
        type_mapper_ptr = @cfunction(type_mapper_callback, API.LLVMTypeRef, (API.LLVMTypeRef,Any))
        type_mapper_data = Ref(type_mapper)
    end
    if materializer === nothing
        materializer_ptr = C_NULL
        materializer_data = C_NULL
    else
        materializer_ptr = @cfunction(materializer_callback, API.LLVMValueRef, (API.LLVMValueRef,Any))
        materializer_data = Ref(materializer)
    end
    API.LLVMCloneFunctionInto(new, old, value_map_array, length(value_map), changes, suffix,
                              type_mapper_ptr, type_mapper_data,
                              materializer_ptr, materializer_data)
end

function clone(f::Function; value_map::Dict{Value,Value}=Dict{Value,Value}())
    argtypes = LLVMType[]

    # The user might be deleting arguments to the function by specifying them in
    # the VMap. If so, we need to not add the arguments to the arg ty vector
    for arg in parameters(f)
        if !in(arg, keys(value_map))    # Haven't mapped the argument to anything yet?
            push!(argtypes, llvmtype(arg))
        end
    end

    # Create a new function type...
    vararg = isvararg(eltype(llvmtype(f)))
    fty = FunctionType(return_type(eltype(llvmtype(f))), argtypes; vararg)

    # Create the new function...
    new_f = Function(parent(f), name(f), fty)
    linkage!(new_f, linkage(f))
    # TODO: address space

    # Loop over the arguments, copying the names of the mapped arguments over...
    for (arg, new_arg) in zip(parameters(f), parameters(new_f))
        if !in(arg, keys(value_map))    # Is this argument preserved?
            name!(new_arg, name(arg))   # Copy the name over...
            value_map[arg] = new_arg    # Add mapping to VMap
        end
    end

    clone_into!(new_f, f;
                value_map, changes=API.LLVMCloneFunctionChangeTypeLocalChangesOnly)

   return new_f
end


text:
set_id() = "ID_" * randstring(10)

#function deal_with_NA(data)
#    res = Array{Any,0}
#    for i = 1:nrow(data)
#        tmp = isna(data[i,1])? '-' : data[i,1]
#        push!(res,tmp)
#    end
#    return res
#end

function ec_init()
    opt = Dict(
        "title" => Dict{ASCIIString,Any}(),
        "legend" => Dict{ASCIIString,Any}(),
        "toolbox" =>
            Dict{ASCIIString,Any}(
                "show" => true,
                "feature" =>
                    Dict{ASCIIString,Any}(
                        "mark" => Dict{ASCIIString,Bool}("show"=>true),
                        "dataView" => Dict{ASCIIString,Bool}("show"=>true,"readOnly" => false),
                        "restore" =>  Dict{ASCIIString,Bool}("show"=>true),
                        "saveAsImage" => Dict{ASCIIString,Bool}("show"=>true)
                    )
            ),
        "tooltip" => Dict{ASCIIString,ASCIIString}(),
        "dataZoom" => Dict{ASCIIString,Any}(),
        "series" => Array(Dict{ASCIIString,Any}, 0)
        )

    return opt
end

function ec_init_cartesian()
    opt = ec_init()
    opt["xAxis"] = Array(Dict{ASCIIString,Any}, 0)
    opt["yAxis"] = Array(Dict{ASCIIString,Any}, 0)

    return opt
end

function ec_set_visible_toolbox(opt::Dict, visible::Bool)
    opt["toolbox"]["show"] = visible
end

function ec_set_color(opt::Dict, colorlist::Vector{ASCIIString})
    opt["color"] = colorlist
end

function ec_set_title{S<:String}(opt::Dict, title::S)
    opt["title"]["text"] = title
end


text:
const expander = Core.atdoc
const setter! = Core.atdoc!

"""
$(:SIGNATURES)

Set the docstring expander function to first call `func` before calling the default expander.

To remove a hook that has been applied using this method call [`hook!()`](@ref).
"""
hook!(func) = setter!((args...) -> expander(func(args...)...))

"""
$(:SIGNATURES)

Reset the docstring expander to only call the default expander function. This clears any
'hook' that has been set using [`hook!(func)`](@ref).
"""
hook!() = setter!(expander)

"""
$(:SIGNATURES)

Defines a docstring template that will be applied to all docstrings in a module that match
the specified category or tuple of categories.

# Examples

```julia
@template DEFAULT =
    \"""
    \$(SIGNATURES)
    \$(DOCSTRING)
    \"""
```

`DEFAULT` is the default template that is applied to a docstring if no other template
definitions match the documented expression. The `DOCSTRING` abbreviation is used to mark
the location in the template where the actual docstring body will be spliced into each
docstring.

```julia
@template (FUNCTIONS, METHODS, MACROS) =
    \"""
    \$(SIGNATURES)
    \$(DOCSTRING)
    \$(METHODLIST)
    \"""
```

A tuple of categories can be specified when a docstring template should be used for several
different categories.

```julia
@template MODULES = ModName
```

The template definition above will define a template for module docstrings based on the
template for modules found in module `ModName`.

!!! note

    Supported categories are `DEFAULT`, `FUNCTIONS`, `METHODS`, `MACROS`, `TYPES`,
    `MODULES`, and `CONSTANTS`.

"""
macro template(ex)
    template(__source__, __module__, ex)
end

const TEMP_SYM = gensym("templates")

function template(src::LineNumberNode, mod::Module, ex::Expr)
    Meta.isexpr(ex, :(=), 2) || error("invalid `@template` syntax.")
    template(src, mod, ex.args[1], ex.args[2])
end

function template(source::LineNumberNode, mod::Module, tuple::Expr, docstr::Union{Symbol, Expr})
    Meta.isexpr(tuple, :tuple) || error("invalid `@template` syntax on LHS.")
    isdefined(mod, TEMP_SYM) || Core.eval(mod, :(const $(TEMP_SYM) = $(Dict{Symbol, Vector}())))
    local block = Expr(:block)
    for category in tuple.args
        local key = Meta.quot(category)
        local vec = Meta.isexpr(docstr, :string) ?
            Expr(:vect, docstr.args...) : :($(docstr).$(TEMP_SYM)[$(key)])
        push!(block.args, :($(TEMP_SYM)[$(key)] = $(vec)))
    end
    push!(block.args, nothing)
    return esc(block)
end

function template(src::LineNumberNode, mod::Module, sym::Symbol, docstr::Union{Symbol, Expr})
    template(src, mod, Expr(:tuple, sym), docstr)
end

# The signature for the atdocs() calls changed in v0.7
# On v0.6 and below it seems it was assumed to be (docstr::String, expr::Expr), but on v0.7
# it is (source::LineNumberNode, mod::Module, docstr::String, expr::Expr)
function template_hook(source::LineNumberNode, mod::Module, docstr, expr::Expr)
    # During macro expansion we only need to wrap docstrings in special
    # abbreviations that later print out what was before and after the
    # docstring in it's specific template. This is only done when the module
    # actually defines templates.
    if isdefined(mod, TEMP_SYM)
        dict = getfield(mod, TEMP_SYM)
        # We unwrap interpolated strings so that we can add the `:before` and
        # `:after` abbreviations. Otherwise they're just left as is.
        unwrapped = Meta.isexpr(docstr, :string) ? docstr.args : [docstr]
        before, after = Template{:before}(dict), Template{:after}(dict)
        # Rebuild the original docstring, but with the template abbreviations
        # surrounding it.
        docstr = Expr(:string, before, unwrapped..., after)
    end
    return (source, mod, docstr, expr)
end

function template_hook(docstr, expr::Expr)
    source, mod, docstr, expr::Expr = template_hook(LineNumberNode(0), current_module(), docstr, expr)
    docstr, expr
end

template_hook(args...) = args

get_template(t::Dict, k::Symbol) = haskey(t, k) ? t[k] : get(t, :DEFAULT, Any[DOCSTRING])


text:
module Wrappers

using CommonRLInterface

export
    AbstractWrapper,
    QuickWrapper

"""
    AbstractWrapper

Abstract base class for environment wrappers. For a subtype of `AbstractWrapper`, all CommonRLInterface functions will be forwarded to the wrapped environment defined by `wrapped_env`.

Interface functions can be selectively overridden for the new wrapper type. `provided` and optional functions will be handled correctly by default.

# Example

```julia
struct MyActionWrapper{E} <: AbstractWrapper
    env::E
end

# Any subclass of AbstractWrapper MUST implement wrapped_env
CommonRLInterface.wrapped_env(w::MyActionWrapper) = w.env

# Now all CommonRLFunctions functions are forwarded
w = MyActionWrapper(env)
observe(w) # will return an observation from env
actions(w) # will return the action space from env

# The actions function for the wrapper can be overridden
CommonRLInterface.actions(w::MyActionWrapper) = [-1, 1]
actions(w) # will now return [-1, 1]
"""
abstract type AbstractWrapper end

"""
Indicate what the wrapped environment is for an AbstractWrapper (see the AbstractWrapper docstring)
"""
function wrapped_env end

macro forward_to_wrapped(f)
    return :($f(w::AbstractWrapper, args...; kwargs...) = $f(wrapped_env(w), args...; kwargs...))
end

@forward_to_wrapped CommonRLInterface.reset!
@forward_to_wrapped CommonRLInterface.actions
@forward_to_wrapped CommonRLInterface.observe
@forward_to_wrapped CommonRLInterface.act!
@forward_to_wrapped CommonRLInterface.terminated

@forward_to_wrapped CommonRLInterface.player

@forward_to_wrapped CommonRLInterface.render
@forward_to_wrapped CommonRLInterface.state
@forward_to_wrapped CommonRLInterface.setstate!
@forward_to_wrapped CommonRLInterface.valid_actions
@forward_to_wrapped CommonRLInterface.valid_action_mask
@forward_to_wrapped CommonRLInterface.observations
# not straightforward to provide clone

CommonRLInterface.provided(f::Function, w::AbstractWrapper, args...) = provided(f, wrapped_env(w), args...)
CommonRLInterface.provided(::typeof(clone), w::AbstractWrapper, args...) = false

include("quick_wrapper.jl")

end


text:
module ProtonFunctions

using QuadGK
using Trapz

export adiabatic_energy_loss, interaction_rate, hubbles_law, survival_probability

"""
    `adiabatic_energy_loss(lorentz_initial, z_initial, z_final)`

Calculates the Lorentz factor of a particle after it propagates from a redshift z_initial to a redshift z_final.
The particle loses energy due to the adiabatic expansion of the universe.

Based on eq. (2-4) from arXiv:1505.01347 https://arxiv.org/abs/1505.01347
"""
function adiabatic_energy_loss(lorentz_initial::Real, z_initial::Real, z_final::Real)
    integral, err = quadgk(z -> 1/(1+z), z_initial, z_final) # numerical integration using QuadGK.jl
    energy_loss = ℯ^(-integral)
    return lorentz_initial/energy_loss
end

"""
    `interaction_rate(Γ, ϵ, σ)`
Calculates the interaction rate of protons with photons for a given proton Lorentz factor, redshift, photon energy, and interaction cross section.

Inputs:
- Γ: proton Lorentz factor
- z: redshift
- ϵ: photon energy (in the proton rest frame) (eV)
- σ: pγ interaction cross-section (micrometer^2)

Source: Eq. (10.9) from Cosmic Rays and Particle Physics, 2nd Edition, by Thomas K. Gaisser, Ralph Engel, and Elisa Resconi.
"""
function interaction_rate(Γ::Real, z::Real, ϵ::Vector{T1}, σ::Vector{T2}) where {T1<:Number, T2<:Number}
    c = (2.99792e8)*1e9 # speed of light, nanometer/s
    kT = 0.234822*1e-9 # boltzmann constant * z=0 CMB temperature, MeV
    ħc = .197327*1e-3 # reduced planck constant * speed of light, MeV * nanometer
    M = @. ϵ*σ*(-log(1-ℯ^(-ϵ/(2*Γ*kT))))
    # trapezoidal integration of M over the array of ϵ
    I = trapz(ϵ, M)
    @assert I >= 0.0 # An interaction rate will never be negative
    rate_at_zero_redshift = (c*kT)/((2*π^2)*(ħc^3)*(Γ^2)) * I
    return rate_at_zero_redshift*(1+z)^3
end

"""
    'hubbles_law(z)'

Calculates dt/dz at the input redshift z. Source: eq. (3) from arXiv:1505.01347 https://arxiv.org/abs/1505.01347
"""
function hubbles_law(z::Real)
    H0 = 2.198e-18 # Hubble constant in s^-1
    Ωm = 0.3 # matter density parameter
    ΩΛ = 0.7 # dark energy density parameter
    return -1/(H0*(1+z)*sqrt((1+z^3)*Ωm + ΩΛ))
end

"""
    'survival_probability(p_prev, rate_prev, z_prev, rate, z)'

The probability that the proton has not undergone a reaction at each redshift step, approximated
using the trapezoidal rule. Source: Eq. (7) arXiv:1505.01347 https://arxiv.org/abs/1505.01347

Inputs:
- p_prev = the survival probability at the previous redshift step
- rate_prev = the interaction rate at the previous redshift step
- z_prev = the redshift at the previous redshift step
- rate = the interaction rate at the current redshift step
- z = the redshift at the current redshift step
"""
function survival_probability(p_prev::Real, rate_prev::Real, z_prev::Real, rate::Real, z::Real)
    @assert p_prev >= 0 # probabilites cannot be negative
    logp = log(p_prev) - (rate_prev*hubbles_law(z_prev) + rate*hubbles_law(z))*(0.5)*(z - z_prev)
    return ℯ^logp
end

end


text:
# test LaMEM i/o routines
using Test, GeophysicalModelGenerator

# Load LaMEM input file:
Grid        =   ReadLaMEM_InputFile("test_files/SaltModels.dat")
@test Grid.X[10] ≈ -2.40625

# Transfer into ParaviewData struct:
Phases      =   zeros(Int32,   size(Grid.X));
Temp        =   zeros(Float64, size(Grid.X));
Model3D     =   CartData(Grid, (Phases=Grid.Z,));
@test  Model3D.y.val[100]==-1.9375km

# Read Partitioning file:
PartitioningFile = "test_files/ProcessorPartitioning_4cpu_1.2.2.bin"
Nprocx,Nprocy,Nprocz, xc,yc,zc, nNodeX,nNodeY,nNodeZ = GetProcessorPartitioning(PartitioningFile)
@test Nprocz==2
@test yc[2]==0.0

# Save serial output
Save_LaMEMMarkersParallel(Model3D, verbose=false)

# Save parallel output
Save_LaMEMMarkersParallel(Model3D, PartitioningFile=PartitioningFile, verbose=false)

# Test creating model setups
Grid        =   ReadLaMEM_InputFile("test_files/Subduction2D_FreeSlip_Particles_Linear_DirectSolver.dat")
Phases      =   zeros(Int32,   size(Grid.X));

# constant T
Temp        =   ones(Float64, size(Grid.X))*1350;
AddBox!(Phases,Temp,Grid, xlim=(0,500), zlim=(-50,0), phase=ConstantPhase(3), DipAngle=10, T=ConstantTemp(1000))
@test sum(Temp) == 1.1905107e9

# Add a layer above the slab with a different phase but no thermal structure
AddBox!(Phases,Temp,Grid, xlim=(0,500), zlim=(0,20), phase=ConstantPhase(1), DipAngle=10, Origin=(0,0,0))
@test sum(Temp) == 1.1905107e9

# Linear T
Temp        =   ones(Float64, size(Grid.X))*1350;
AddBox!(Phases,Temp,Grid, xlim=(0,500), zlim=(-50,0), phase=ConstantPhase(3), DipAngle=10, T=LinearTemp(Tbot=1350, Ttop=200))
@test sum(Temp) == 1.1881296265169694e9

# Halfspace cooling T structure
Phases      =   zeros(Int32,   size(Grid.X));
Temp        =   ones(Float64, size(Grid.X))*1350;
AddBox!(Phases,Temp,Grid, xlim=(0,500), zlim=(-500,0), phase=LithosphericPhases(Layers=[15 15 250], Phases=[1 2 3 0], Tlab=1250), DipAngle=10, T=HalfspaceCoolingTemp(Age=20, Adiabat=0.3))
@test sum(Temp) == 1.1942982365477426e9

# Mid-oceanic ridge cooling temperature structure
Phases      =   zeros(Int32,   size(Grid.X));
Temp        =   ones(Float64, size(Grid.X))*1350;
AddBox!(Phases,Temp,Grid, xlim=(0,500), zlim=(-500,-20), phase=LithosphericPhases(Layers=[15 15 250], Phases=[1 2 3 0],Tlab=1250), DipAngle=10, T=SpreadingRateTemp(MORside="right", SpreadingVel=3))
@test sum(Temp) == 1.189394358568891e9

Model3D     =   CartData(Grid, (Phases=Phases,Temp=Temp));
Write_Paraview(Model3D,"LaMEM_ModelSetup")                  # Save model to paraview

# Test writing a LaMEM topography file
X,Y,Z = XYZGrid(-20:20,-10:10,0);
Z = cos.(2*pi.*X./5).*cos.(2*pi.*Y./10)

Topo = CartData(X,Y,Z,(Topography=Z,))
@test Save_LaMEMTopography(Topo, "test_topo.dat")==nothing
rm("test_topo.dat")

# Test adding geometric primitives
Grid    = ReadLaMEM_InputFile("test_files/GeometricPrimitives.dat")
Phases  = zeros(Int32,size(Grid.X));
Temp    = zeros(Float64,size(Grid.X));
AddSphere!(Phases,Temp,Grid, cen=(0,0,-6), radius=2, phase=ConstantPhase(1), T=ConstantTemp(800))
@test Phases[55,55,55] == 1
@test Phases[56,56,56] == 0
@test Temp[44,52,21]   == 800.0
@test Temp[44,52,20]   == 0.0

AddEllipsoid!(Phases,Temp,Grid, cen=(-2,-1,-7), axes=(1,2,3), StrikeAngle=90, DipAngle=45, phase=ConstantPhase(2), T=ConstantTemp(600))
@test Phases[11,37,28] == 2
@test Phases[10,37,28] == 0
@test Temp[31,58,18]   == 600.0
@test Temp[31,59,18]   == 0.0

AddCylinder!(Phases,Temp,Grid, base=(0,0,-5), cap=(3,3,-2), radius=1.5, phase=ConstantPhase(3), T=ConstantTemp(400))
@test Phases[55,65,75] == 3
@test Phases[54,65,75] == 0
@test Temp[55,46,45]   == 400.0
@test Temp[55,45,45]   == 800.0


text:
# Julia wrapper for header: /usr/include/scip/pub_benderscut.h
# Automatically generated using Clang.jl wrap_c

function SCIPbenderscutComp(elem1, elem2)
    ccall((:SCIPbenderscutComp, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPbenderscutCompName(elem1, elem2)
    ccall((:SCIPbenderscutCompName, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPbenderscutGetData(benderscut)
    ccall((:SCIPbenderscutGetData, libscip), Ptr{SCIP_BENDERSCUTDATA}, (Ptr{SCIP_BENDERSCUT},), benderscut)
end

function SCIPbenderscutSetData(benderscut, benderscutdata)
    ccall((:SCIPbenderscutSetData, libscip), Cvoid, (Ptr{SCIP_BENDERSCUT}, Ptr{SCIP_BENDERSCUTDATA}), benderscut, benderscutdata)
end

function SCIPbenderscutGetName(benderscut)
    ccall((:SCIPbenderscutGetName, libscip), Cstring, (Ptr{SCIP_BENDERSCUT},), benderscut)
end

function SCIPbenderscutGetDesc(benderscut)
    ccall((:SCIPbenderscutGetDesc, libscip), Cstring, (Ptr{SCIP_BENDERSCUT},), benderscut)
end

function SCIPbenderscutGetPriority(benderscut)
    ccall((:SCIPbenderscutGetPriority, libscip), Cint, (Ptr{SCIP_BENDERSCUT},), benderscut)
end

function SCIPbenderscutGetNCalls(benderscut)
    ccall((:SCIPbenderscutGetNCalls, libscip), Clonglong, (Ptr{SCIP_BENDERSCUT},), benderscut)
end

function SCIPbenderscutGetNFound(benderscut)
    ccall((:SCIPbenderscutGetNFound, libscip), Clonglong, (Ptr{SCIP_BENDERSCUT},), benderscut)
end

function SCIPbenderscutIsInitialized(benderscut)
    ccall((:SCIPbenderscutIsInitialized, libscip), UInt32, (Ptr{SCIP_BENDERSCUT},), benderscut)
end

function SCIPbenderscutGetSetupTime(benderscut)
    ccall((:SCIPbenderscutGetSetupTime, libscip), Cdouble, (Ptr{SCIP_BENDERSCUT},), benderscut)
end

function SCIPbenderscutGetTime(benderscut)
    ccall((:SCIPbenderscutGetTime, libscip), Cdouble, (Ptr{SCIP_BENDERSCUT},), benderscut)
end

function SCIPbenderscutGetAddedConss(benderscut, addedconss, naddedconss)
    ccall((:SCIPbenderscutGetAddedConss, libscip), SCIP_RETCODE, (Ptr{SCIP_BENDERSCUT}, Ptr{Ptr{Ptr{SCIP_CONS}}}, Ptr{Cint}), benderscut, addedconss, naddedconss)
end

function SCIPbenderscutGetAddedCuts(benderscut, addedcuts, naddedcuts)
    ccall((:SCIPbenderscutGetAddedCuts, libscip), SCIP_RETCODE, (Ptr{SCIP_BENDERSCUT}, Ptr{Ptr{Ptr{SCIP_ROW}}}, Ptr{Cint}), benderscut, addedcuts, naddedcuts)
end

function SCIPbenderscutIsLPCut(benderscut)
    ccall((:SCIPbenderscutIsLPCut, libscip), UInt32, (Ptr{SCIP_BENDERSCUT},), benderscut)
end

function SCIPbenderscutSetEnabled(benderscut, enabled)
    ccall((:SCIPbenderscutSetEnabled, libscip), Cvoid, (Ptr{SCIP_BENDERSCUT}, UInt32), benderscut, enabled)
end


text:
# # MNIST Classification using Neural ODEs

# ## Package Imports
using Lux
using Pkg #hide
Pkg.activate(joinpath(dirname(pathof(Lux)), "..", "examples")) #hide
using ComponentArrays, CUDA, DiffEqSensitivity, NNlib, Optimisers, OrdinaryDiffEq, Random,
      Statistics, Zygote
import MLDatasets: MNIST
import MLDataUtils: convertlabel, LabelEnc
import MLUtils: DataLoader, splitobs
CUDA.allowscalar(false)

# ## Loading MNIST
## Use MLDataUtils LabelEnc for natural onehot conversion
function onehot(labels_raw)
    convertlabel(LabelEnc.OneOfK, labels_raw, LabelEnc.NativeLabels(collect(0:9)))
end

function loadmnist(batchsize, train_split)
    ## Load MNIST: Only 1500 for demonstration purposes
    N = 1500
    imgs = MNIST.traintensor(1:N)
    labels_raw = MNIST.trainlabels(1:N)

    ## Process images into (H,W,C,BS) batches
    x_data = Float32.(reshape(imgs, size(imgs, 1), size(imgs, 2), 1, size(imgs, 3)))
    y_data = onehot(labels_raw)
    (x_train, y_train), (x_test, y_test) = splitobs((x_data, y_data); at=train_split)

    return (
            ## Use DataLoader to automatically minibatch and shuffle the data
            DataLoader(collect.((x_train, y_train)); batchsize=batchsize, shuffle=true),
            ## Don't shuffle the test data
            DataLoader(collect.((x_test, y_test)); batchsize=batchsize, shuffle=false))
end

# ## Define the Neural ODE Layer
#
# The NeuralODE is a ContainerLayer. It stores a `model` and the parameters and states of the NeuralODE is
# same as that of the underlying model.
struct NeuralODE{M <: Lux.AbstractExplicitLayer, So, Se, T, K} <:
       Lux.AbstractExplicitContainerLayer{(:model,)}
    model::M
    solver::So
    sensealg::Se
    tspan::T
    kwargs::K
end

function NeuralODE(model::Lux.AbstractExplicitLayer;
                   solver=Tsit5(),
                   sensealg=InterpolatingAdjoint(; autojacvec=ZygoteVJP()),
                   tspan=(0.0f0, 1.0f0),
                   kwargs...)
    return NeuralODE(model, solver, sensealg, tspan, kwargs)
end

function (n::NeuralODE)(x, ps, st)
    function dudt(u, p, t)
        u_, st = n.model(u, p, st)
        return u_
    end
    prob = ODEProblem{false}(ODEFunction{false}(dudt), x, n.tspan, ps)
    return solve(prob, n.solver; sensealg=n.sensealg, n.kwargs...), st
end

function diffeqsol_to_array(x::ODESolution{T, N, <:AbstractVector{<:CuArray}}) where {T, N}
    dropdims(gpu(x); dims=3)
end
diffeqsol_to_array(x::ODESolution) = dropdims(Array(x); dims=3)

# ## Create and Initialize the Neural ODE Layer
function create_model()
    ## Construct the Neural ODE Model
    model = Chain(FlattenLayer(),
                  Dense(784, 20, tanh),
                  NeuralODE(Chain(Dense(20, 10, tanh), Dense(10, 10, tanh),
                                  Dense(10, 20, tanh));
                            save_everystep=false,
                            reltol=1.0f-3,
                            abstol=1.0f-3,
                            save_start=false),
                  diffeqsol_to_array,
                  Dense(20, 10))

    rng = Random.default_rng()
    Random.seed!(rng, 0)

    ps, st = Lux.setup(rng, model)
    ps = ComponentArray(ps) |> gpu
    st = st |> gpu

    return model, ps, st
end

# ## Define Utility Functions
get_class(x) = argmax.(eachcol(x))

logitcrossentropy(ŷ, y) = mean(-sum(y .* logsoftmax(ŷ); dims=1))

function loss(x, y, model, ps, st)
    ŷ, st = model(x, ps, st)
    return logitcrossentropy(ŷ, y), st
end

function accuracy(model, ps, st, dataloader)
    total_correct, total = 0, 0
    st = Lux.testmode(st)
    iterator = CUDA.functional() ? CuIterator(dataloader) : dataloader
    for (x, y) in iterator
        target_class = get_class(cpu(y))
        predicted_class = get_class(cpu(model(x, ps, st)[1]))
        total_correct += sum(target_class .== predicted_class)
        total += length(target_class)
    end
    return total_correct / total
end

# ## Training
function train()
    model, ps, st = create_model()

    ## Training
    train_dataloader, test_dataloader = loadmnist(128, 0.9)

    opt = Optimisers.ADAM(0.001f0)
    st_opt = Optimisers.setup(opt, ps)

    ### Warmup the Model
    img, lab = gpu(train_dataloader.data[1][:, :, :, 1:1]),
               gpu(train_dataloader.data[2][:, 1:1])
    loss(img, lab, model, ps, st)
    (l, _), back = pullback(p -> loss(img, lab, model, p, st), ps)
    back((one(l), nothing))

    ### Lets train the model
    nepochs = 9
    for epoch in 1:nepochs
        stime = time()
        iterator = CUDA.functional() ? CuIterator(train_dataloader) : train_dataloader
        for (x, y) in iterator
            (l, _), back = pullback(p -> loss(x, y, model, p, st), ps)
            gs = back((one(l), nothing))[1]
            st_opt, ps = Optimisers.update(st_opt, ps, gs)
        end
        ttime = time() - stime

        println("[$epoch/$nepochs] \t Time $(round(ttime; digits=2))s \t Training Accuracy: " *
                "$(round(accuracy(model, ps, st, train_dataloader) * 100; digits=2))% \t " *
                "Test Accuracy: $(round(accuracy(model, ps, st, test_dataloader) * 100; digits=2))%")
    end
end

train()


text:
using FluidFields
using Test


text:
# # Power Flow

# `PowerSystems.jl` provides the capability to run a power flow with the intention of
# providing a valid initial AC operating point to the system.

# The power flow tool is not meant for analytics where the principal goal is to determine
# if the system has settings that produce a feasible AC solution. This power flow routine
# does not check for reactive power limits or other limiting mechanisms in the grid, and can
# therefore be used to check for solver convergence - making no guarantees
# of the solution feasibility.

# The power flow solver uses [NLsolve.jl](https://github.com/JuliaNLSolvers/NLsolve.jl) under
# the hood and takes any keyword argument accepted by NLsolve. The solver uses the current
# operating point in the buses to provide the initial guess.

# **Limitations**: The PowerFlow solver doesn't support systems with HVDC lines or
# Phase Shifting transformers yet. The power flow solver can't handle systems with islands.

# Check section [Power Flow](@ref pf) for detailed usage instructions

using PowerSystems
const PSY = PowerSystems

DATA_DIR = "../../../data" #hide
system_data = System(joinpath(DATA_DIR, "matpower/case14.m"))

# `PowerSystems.jl` has two modes of using the power flow solver.

# 1. Solving the power flow for the current operating point in the system.
#    Takes the data in the buses, the `active_power` and `reactive_power` fields
#    in the static injection devices. Returns a dictionary with results in a DataFrame that
#    can be exported or manipulated as needed.
#
# 2. Solves the power flow and updated the devices in the system to the operating condition.
#    This model will update the values of magnitudes and angles in the system's buses. It
#    also updates the active and reactive power flows in the branches and devices connected
#    to PV buses. It also updates the active and reactive power of the injection devices
#    connected to the Slack bus, and updates only the reactive power of the injection devices
#    connected to PV buses. If multiple devices are connected to the same bus, the power is
#    divided proportional to the base power.
#    This utility is useful to initialize systems before serializing or checking the
#    addition of new devices is still AC feasible.

# Solving the power  flow with mode 1:

results = solve_powerflow(system_data)
results["bus_results"]

# Solving the power  flow with mode 2:

# Before running the power flow command these are the values of the
# voltages:

for b in get_components(Bus, system_data)
    println("$(get_name(b)) - Magnitude $(get_magnitude(b)) - Angle (rad) $(get_angle(b))")
end

# [`solve_powerflow!`](@ref) return true or false to signal the successful result of the power
# flow. This enables the integration of a power flow into functions and use the return as check.
# For instance, initializing dynamic simulations. Also, because [`solve_powerflow!`](@ref) uses
# [NLsolve.jl](https://github.com/JuliaNLSolvers/NLsolve.jl) all the parameters used for NLsolve
# are also available for [`solve_powerflow!`](@ref)

solve_powerflow!(system_data; finite_diff = true, method = :newton)

# After running the power flow command this are the values of the
# voltages:

for b in get_components(Bus, system_data)
    println("$(get_name(b)) - Magnitude $(get_magnitude(b)) - Angle (rad) $(get_angle(b))")
end


text:
using Pkg
Pkg.activate(".")

verbose = true
using Revise
if verbose   println("# Loading RvSpecML")    end
using RvSpectML
if verbose   println("# Loading other packages")    end
using DataFrames, Query

expres_data_path = "."
 target_subdir = "101501"

 # TODO: USER:  Either the paths that specify where datafiles are stored here or in examples/data_paths.jl
 if verbose   println("# Creating manifest of files to process.")    end
 if isdir("examples")
    cd("examples")
 end
 if isfile("data_paths.jl")
      include("data_paths.jl")
 end
 expres_data_path

num_spectra_to_use = 20
 #bjd_first_good = 2458745.1296134139
 #bjd_last_good = 2458745.283
 df_files = EXPRES.make_manifest(joinpath(expres_data_path,target_subdir))
 df_files_use = df_files |>
   @filter( _.target == "101501" ) |>
   # @take(num_spectra_to_use) |>
   DataFrame

if verbose println("# Reading in FITS files.")  end
@time expres_data = map(EXPRES.read_data,eachrow(df_files_use))


text:
function add_noise!(layout::RQCLayout{VT, PT, SCGate}, v, d) where {VT, PT}
    if gates(layout, v)[d] !== FSim
        gates(layout, v)[d] = Noise
        return true
    end
    return false
end

function update_noise!(layout::RQCLayout{VT, PT, SCGate}; approx = false) where {VT, PT}
    to_update = true
    while to_update
        to_update = false
        for v in vertices(layout)
            gates_v = gates(layout, v)
            for i in 1:length(gates_v)
                gates_v[i] === Noise || continue
                for j in (i-1, i+1)
                    1 <= j <= length(gates_v) || continue
                    if gates_v[j] === FSim
                        u = partner(layout, v, depth_to_cycle(j))
                        gates_u = gates(layout, u)
                        if Noise in (approx ? (gates_u[i], gates_v[2j-i], gates_u[2j-i]) : (gates_u[i], gates_v[2j-i]))
                            gates_v[j] = Noise
                            gates_u[j] = Noise
                            # u == 53 && println("FSim @ $(depth_to_cycle(j)) from $v")
                            # v == 53 && println("FSim @ $(depth_to_cycle(j))")
                            to_update = true
                        end
                    elseif gates_v[j] in (Id, Single)
                        gates_v[j] = Noise
                        # v == 53 && println("Single @ $(depth_to_cycle(j))")
                        to_update = true
                    end
                end
            end
        end
    end
    return layout
end

function update_id!(layout::RQCLayout{VT, PT, SCGate}) where {VT, PT}
    for v in vertices(layout)
        gates_v = gates(layout, v)
        for i = 1:(length(gates_v)-1)
            if gates_v[i] === Noise && gates_v[i+1] === Noise
                i > 1 && (gates_v[i] = Id)
            end
        end
    end
    for v in vertices(layout)
        gates_v = gates(layout, v)
        if gates_v[1] === Noise && gates_v[2] === Noise
            gates_v[2] = Id
        end
    end
    return layout
end

simplify!(layout::RQCLayout{VT, PT, SCGate}; approx = false) where {VT, PT} = update_noise!(layout; approx = approx)
function simplify!(layout::RQCLayout{VT, PT, SCGate}, cuts; approx = false) where {VT, PT}
    for c in cuts
        add_noise!(layout, c[1], c[2])
    end
    return simplify!(layout; approx = approx)
end

function add_noise_greedy!(layout::RQCLayout, vs = collect(vertices(layout)),
    flux_min = 15, out_min = 25)
    update_noise!(layout)
    to_update = true
    while to_update
        to_update = false
        for d in 1:depth(layout), v in vs
            if gates(layout, v)[d] in (Single, Id)
                lo_new = copy(layout)
                gates(lo_new, v)[d] = Noise
                update_noise!(lo_new)
                if flux(lo_new, vs) >= flux_min &&
                        sum(gates(lo_new, v)[end] !== Noise for v in vs) >= out_min
                    add_noise!(layout, v, d)
                    update_noise!(layout)
                    to_update = true
                end
            end
        end
    end
    return layout
end


text: push(xs::StaticArrays.StaticVector, x) = vcat(xs, StaticArrays.SVector(x))


text:
#include("src/book.jl")
#include("src/journal.jl")

include("book.jl")
include("journal.jl")
include("journal.jl")

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
# 1:2*length(layerconfig) -> forward lstm
# 2*lenghth(layerconfig)+1:4*length(layerconfig) -> bacward lstm
# final 3 or 4 is embedding and final layer prediction parameters
function initparams(atype, layerconfig, embedsize, vocab, winit; single_embedding=false)
    if !single_embedding
        parameters = Array(Any, 4*length(layerconfig)+4)
        parameters[end] = winit * randn(vocab, embedsize) # backward lstm embedding
        parameters[end-1] = winit * randn(vocab, embedsize) # forward lstm embedding
        parameters[end-2] =  winit * randn(layerconfig[end]*2, vocab)# final layer weight
        parameters[end-3] = zeros(1, vocab) # final layer bias
    else
        parameters = Array(Any, 4*length(layerconfig)+3)
        parameters[end] = winit * randn(vocab, embedsize) # embedding shared among forward and backward lstm
        parameters[end-1] =  winit * randn(layerconfig[end], vocab)# final layer weight
        parameters[end-2] = zeros(1, vocab) # final layer bias
    end

    len = length(layerconfig)
    # forward layer parameters
    input = embedsize
    for k=1:len
        parameters[2k-1] = winit * randn(input+layerconfig[k], 4*layerconfig[k])
        parameters[2k] = zeros(1, 4*layerconfig[k])
        parameters[2k][1:layerconfig[k]] = 1 # forget gate bias
        input = layerconfig[k]
    end
    # bacward layer parameters
    input = embedsize
    for k=1:len
        parameters[2k-1+2len] = winit * randn(input+layerconfig[k], 4*layerconfig[k])
        parameters[2k+2len] = zeros(1, 4*layerconfig[k])
        parameters[2k+2len][1:layerconfig[k]] = 1 # forget gate bias
        input = layerconfig[k]
    end
    return map(p->convert(atype, p), parameters)
end

# 1:2*length(layerconfig) -> forward lstm
# 2*length(layerconfig)+1:4*length(layerconfig) -> backward lstm
function initstates(atype, layerconfig, batchsize)
    states = Array(Any, 4*length(layerconfig))
    len = length(layerconfig)

    # forward layer initialization
    for k=1:len
        states[2k-1] = zeros(batchsize, layerconfig[k])
        states[2k] = zeros(batchsize, layerconfig[k])
    end

    # backward layer initialization
    for k=1:len
        states[2k-1+2len] = zeros(batchsize, layerconfig[k])
        states[2k+2len] = zeros(batchsize, layerconfig[k])
    end
    return map(s->convert(atype, s), states)
end

function resetstate!(states, atype)
    for k=1:length(states)
        states[k] = zeros(size(states))
        convert(atype, states[k])
    end
end

function lstm(weight, bias, hidden, cell, input)
    gates   = hcat(input,hidden) * weight .+ bias
    hsize   = size(hidden,2)
    forget  = sigm(gates[:,1:hsize])
    ingate  = sigm(gates[:,1+hsize:2hsize])
    outgate = sigm(gates[:,1+2hsize:3hsize])
    change  = tanh(gates[:,1+3hsize:end])
    cell    = cell .* forget + ingate .* change
    hidden  = outgate .* tanh(cell)
    return (hidden,cell)
end

# multi layer lstm forward function
function forward(parameters, states, input)
    x = input
    for i=1:2:length(states)
        (states[i], states[i+1]) = lstm(parameters[i], parameters[i+1], states[i], states[i+1], x)
        x = states[i]
    end
    return x
end

# forward reads from 1 to n-1, backward reads from n to 2
# that loss assumes double embedding
# sequence corresponds the n minitbatches in the literature
function loss(parameters, states, sequence)
    total = 0.0
    count = 0
    atype = typeof(AutoGrad.getval(parameters[1]))
    hlayers = length(states) / 4
    hlayers = convert(Int, hlayers)

    # forward lstm
    fhiddens = Array(Any, length(sequence))
    t1 = states[1:2*hlayers]
    for i=1:length(sequence)-1
        input = convert(atype, sequence[i])
        x = input * parameters[end-1]
        fhiddens[i+1] = forward(parameters[1:2*hlayers], t1, x)
        #(prev_1 == t1) && println("State is not being modified")
        #(fhiddens[i+1] == t1[1]) && println("I get the correct hidden")
        #prev_1 = t1[1:end]
    end
    fhiddens[1] = convert(atype, zeros(size(fhiddens[2])))

    # backward lstm
    bhiddens = Array(Any, length(sequence))
    t2 = states[2*hlayers+1:4*hlayers]
    for i=length(sequence):-1:2
        input = convert(atype, sequence[i])
        x = input * parameters[end]
        bhiddens[i-1] = forward(parameters[2*hlayers+1:4*hlayers], t2, x)
    end
    bhiddens[end] = convert(atype, zeros(size(bhiddens[2])))

    # merge layer
    for i=1:length(fhiddens)
        ypred = hcat(fhiddens[i], bhiddens[i]) * parameters[end-2] .+ parameters[end-3]
        ynorm = logp(ypred, 2)
        ygold = convert(atype, sequence[i])
        count += size(ygold, 1)
        total += sum(ygold .* ynorm)
    end
    return - total / count
end

lossgradient = grad(loss)


text:
module Day2

export get_inputs, get_solution1, get_solution2

## Input getting
# Functions
function parse_line(str)
    policy, password = split(str, ": ")
    lohi, letter = split(policy, " ")
    letter = only(letter)
    lo, hi = parse.(Int, split(lohi, "-"))
    return (; lo, hi, letter, password)
end

read_input(io) = parse_line.(readlines(io))

# Inputs
function get_inputs()
    test_input1 = test_input2 = read_input(
        IOBuffer(
            """
            1-3 a: abcde
            1-3 b: cdefg
            2-9 c: ccccccccc
            """
        )
    )
    test_output1 = 2
    test_output2 = 1

    data = read_input(joinpath(@__DIR__, "input.txt"))

    return (; test_input1, test_input2, test_output1, test_output2, data)
end

## Solution functions (change these)
# Password validators
validator1(line) = line.lo ≤ count(c->c==line.letter, line.password) ≤ line.hi
function validator2(line)
    lo = line.password[line.lo]
    hi = line.password[line.hi]
    letter = line.letter
    return (lo==letter && !(hi==letter)) || (!(lo==letter) && hi==letter)
end

# Solution-getting closure
get_solution_with(password_validator) = data -> sum(password_validator.(data))

# Part 1
get_solution1 = get_solution_with(validator1)

# Part 2
get_solution2 = get_solution_with(validator2)

end


text:
@time include("total_variation_test.jl")
@time include("fourier_pade_test.jl")


text:
using Oceananigans.Fields: AbstractField

#####
##### The turbulence closure proposed by Leith
#####

struct TwoDimensionalLeith{FT, CR, GM, M} <: AbstractTurbulenceClosure{ExplicitTimeDiscretization}
                  C :: FT
             C_Redi :: CR
               C_GM :: GM
    isopycnal_model :: M

    function TwoDimensionalLeith{FT}(C, C_Redi, C_GM, isopycnal_model) where FT
        C_Redi = convert_diffusivity(FT, C_Redi, Val(false))
        C_GM = convert_diffusivity(FT, C_GM, Val(false))
        return new{FT, typeof(C_Redi), typeof(C_GM), typeof(isopycnal_model)}(C, C_Redi, C_GM)
    end
end

"""
    TwoDimensionalLeith(FT=Float64;
                        C=0.3, C_Redi=1, C_GM=1,
                        isopycnal_model=SmallSlopeIsopycnalTensor())

Return a `TwoDimensionalLeith` type associated with the turbulence closure proposed by
Leith (1965) and Fox-Kemper & Menemenlis (2008) which has an eddy viscosity of the form

```julia
νₑ = (C * Δᶠ)³ * √(|∇ₕ ζ|² + |∇ₕ ∂w/∂z|²)
```

and an eddy diffusivity of the form...

where `Δᶠ` is the filter width, `ζ = ∂v/∂x - ∂u/∂y` is the vertical vorticity,
and `C` is a model constant.

Keyword arguments
=================
  - `C`: Model constant
  - `C_Redi`: Coefficient for down-gradient tracer diffusivity for each tracer.
              Either a constant applied to every tracer, or a `NamedTuple` with fields
              for each tracer individually.
  - `C_GM`: Coefficient for down-gradient tracer diffusivity for each tracer.
            Either a constant applied to every tracer, or a `NamedTuple` with fields
            for each tracer individually.

References
==========
Leith, C. E. (1968). "Diffusion Approximation for Two‐Dimensional Turbulence", The Physics of
    Fluids 11, 671. doi: 10.1063/1.1691968

Fox‐Kemper, B., & D. Menemenlis (2008), "Can large eddy simulation techniques improve mesoscale rich
    ocean models?", in Ocean Modeling in an Eddying Regime, Geophys. Monogr. Ser., vol. 177, pp. 319–337.
    doi: 10.1029/177GM19

Pearson, B. et al. (2017) , "Evaluation of scale-aware subgrid mesoscale eddy models in a global eddy
    rich model", Ocean Modelling 115, 42-58. doi: 10.1016/j.ocemod.2017.05.007
"""
TwoDimensionalLeith(FT=Float64; C=0.3, C_Redi=1, C_GM=1, isopycnal_model=SmallSlopeIsopycnalTensor()) =
    TwoDimensionalLeith{FT}(C, C_Redi, C_GM, isopycnal_model)

function with_tracers(tracers, closure::TwoDimensionalLeith{FT}) where FT
    C_Redi = tracer_diffusivities(tracers, closure.C_Redi)
    C_GM = tracer_diffusivities(tracers, closure.C_GM)

    return TwoDimensionalLeith{FT}(closure.C, C_Redi, C_GM, closure.isopycnal_model)
end

@inline function abs²_∇h_ζ(i, j, k, grid, U)
    vxx = ℑyᵃᶜᵃ(i, j, k, grid, ∂²xᶜᶠᶜ, U.v)
    uyy = ℑxᶜᵃᵃ(i, j, k, grid, ∂²yᶠᶜᶜ, U.u)
    uxy = ℑyᵃᶜᵃ(i, j, k, grid, ∂xᶜᶠᶜ, ∂yᶠᶠᶜ, U.u)
    vxy = ℑxᶜᵃᵃ(i, j, k, grid, ∂xᶠᶜᶜ, ∂yᶜᶜᶜ, U.v)

    return (vxx - uxy)^2 + (vxy - uyy)^2
end

const ArrayOrField = Union{AbstractArray, AbstractField}

@inline ψ²(i, j, k, grid, ψ::Function, args...) = ψ(i, j, k, grid, args...)^2
@inline ψ²(i, j, k, grid, ψ::ArrayOrField, args...) = @inbounds ψ[i, j, k]^2

@inline function abs²_∇h_wz(i, j, k, grid, w)
    wxz² = ℑxᶜᵃᵃ(i, j, k, grid, ψ², ∂xᶠᶜᶜ, ∂zᶜᶜᶜ, w)
    wyz² = ℑyᵃᶜᵃ(i, j, k, grid, ψ², ∂yᶜᶠᶜ, ∂zᶜᶜᶜ, w)
    return wxz² + wyz²
end

@inline νᶜᶜᶜ(i, j, k, grid, closure::TwoDimensionalLeith{FT}, buoyancy, U, C) where FT =
    (closure.C * Δᶠ(i, j, k, grid, closure))^3 * sqrt(  abs²_∇h_ζ(i, j, k, grid, U)
                                              + abs²_∇h_wz(i, j, k, grid, U.w))

#####
##### Abstract Smagorinsky functionality
#####

# Diffusive fluxes for Leith diffusivities

@inline function diffusive_flux_x(i, j, k, grid, closure::TwoDimensionalLeith,
                                  c, ::Val{tracer_index}, clock, diffusivities, C, buoyancy) where tracer_index

    νₑ = diffusivities.νₑ

    @inbounds C_Redi = closure.C_Redi[tracer_index]
    @inbounds C_GM = closure.C_GM[tracer_index]

    νₑⁱʲᵏ = ℑxᶠᵃᵃ(i, j, k, grid, νₑ)

    ∂x_c = ∂xᶠᶜᶜ(i, j, k, grid, c)
    ∂z_c = ℑxzᶠᵃᶜ(i, j, k, grid, ∂zᶜᶜᶠ, c)

    R₁₃ = isopycnal_rotation_tensor_xz_fcc(i, j, k, grid, buoyancy, C, closure.isopycnal_model)

    return - νₑⁱʲᵏ * (                 C_Redi * ∂x_c
                      + (C_Redi - C_GM) * R₁₃ * ∂z_c)
end

@inline function diffusive_flux_y(i, j, k, grid, closure::TwoDimensionalLeith,
                                  c, ::Val{tracer_index}, clock, diffusivities, C, buoyancy) where tracer_index

    νₑ = diffusivities.νₑ

    @inbounds C_Redi = closure.C_Redi[tracer_index]
    @inbounds C_GM = closure.C_GM[tracer_index]

    νₑⁱʲᵏ = ℑyᵃᶠᵃ(i, j, k, grid, νₑ)

    ∂y_c = ∂yᶜᶠᶜ(i, j, k, grid, c)
    ∂z_c = ℑyzᵃᶠᶜ(i, j, k, grid, ∂zᶜᶜᶠ, c)

    R₂₃ = isopycnal_rotation_tensor_yz_cfc(i, j, k, grid, buoyancy, C, closure.isopycnal_model)
    return - νₑⁱʲᵏ * (                  C_Redi * ∂y_c
                             + (C_Redi - C_GM) * R₂₃ * ∂z_c)
end

@inline function diffusive_flux_z(i, j, k, grid, closure::TwoDimensionalLeith,
                                  c, ::Val{tracer_index}, clock, diffusivities, C, buoyancy) where tracer_index

    νₑ = diffusivities.νₑ

    @inbounds C_Redi = closure.C_Redi[tracer_index]
    @inbounds C_GM = closure.C_GM[tracer_index]

    νₑⁱʲᵏ = ℑzᵃᵃᶠ(i, j, k, grid, νₑ)

    ∂x_c = ℑxzᶜᵃᶠ(i, j, k, grid, ∂xᶠᶜᶜ, c)
    ∂y_c = ℑyzᵃᶜᶠ(i, j, k, grid, ∂yᶜᶠᶜ, c)
    ∂z_c = ∂zᶜᶜᶠ(i, j, k, grid, c)

    R₃₁ = isopycnal_rotation_tensor_xz_ccf(i, j, k, grid, buoyancy, C, closure.isopycnal_model)
    R₃₂ = isopycnal_rotation_tensor_yz_ccf(i, j, k, grid, buoyancy, C, closure.isopycnal_model)
    R₃₃ = isopycnal_rotation_tensor_zz_ccf(i, j, k, grid, buoyancy, C, closure.isopycnal_model)

    return - νₑⁱʲᵏ * (
          (C_Redi + C_GM) * R₃₁ * ∂x_c
        + (C_Redi + C_GM) * R₃₂ * ∂y_c
                 + C_Redi * R₃₃ * ∂z_c)
end

const L2D = TwoDimensionalLeith

@inline viscous_flux_ux(i, j, k, grid, closure::L2D, clock, U, K) = - 2 * ν_σᶜᶜᶜ(i, j, k, grid, closure, K.νₑ, Σ₁₁, U.u, U.v, U.w)
@inline viscous_flux_uy(i, j, k, grid, closure::L2D, clock, U, K) = - 2 * ν_σᶠᶠᶜ(i, j, k, grid, closure, K.νₑ, Σ₁₂, U.u, U.v, U.w)
@inline viscous_flux_uz(i, j, k, grid, closure::L2D, clock, U, K) = - 2 * ν_σᶠᶜᶠ(i, j, k, grid, closure, K.νₑ, Σ₁₃, U.u, U.v, U.w)

@inline viscous_flux_vx(i, j, k, grid, closure::L2D, clock, U, K) = - 2 * ν_σᶠᶠᶜ(i, j, k, grid, closure, K.νₑ, Σ₂₁, U.u, U.v, U.w)
@inline viscous_flux_vy(i, j, k, grid, closure::L2D, clock, U, K) = - 2 * ν_σᶜᶜᶜ(i, j, k, grid, closure, K.νₑ, Σ₂₂, U.u, U.v, U.w)
@inline viscous_flux_vz(i, j, k, grid, closure::L2D, clock, U, K) = - 2 * ν_σᶜᶠᶠ(i, j, k, grid, closure, K.νₑ, Σ₂₃, U.u, U.v, U.w)

@inline viscous_flux_wx(i, j, k, grid, closure::L2D, clock, U, K) = - 2 * ν_σᶠᶜᶠ(i, j, k, grid, closure, K.νₑ, Σ₃₁, U.u, U.v, U.w)
@inline viscous_flux_wy(i, j, k, grid, closure::L2D, clock, U, K) = - 2 * ν_σᶜᶠᶠ(i, j, k, grid, closure, K.νₑ, Σ₃₂, U.u, U.v, U.w)
@inline viscous_flux_wz(i, j, k, grid, closure::L2D, clock, U, K) = - 2 * ν_σᶜᶜᶜ(i, j, k, grid, closure, K.νₑ, Σ₃₃, U.u, U.v, U.w)

function calculate_diffusivities!(diffusivity_fields, closure::TwoDimensionalLeith, model)
    arch = model.architecture
    grid = model.grid
    velocities = model.velocities
    tracers = model.tracers
    buoyancy = model.buoyancy

    event = launch!(arch, grid, :xyz,
                    calculate_nonlinear_viscosity!,
                    diffusivity_fields.νₑ, grid, closure, buoyancy, velocities, tracers,
                    dependencies = device_event(arch))

    wait(device(arch), event)

    return nothing
end

"Return the filter width for a Leith Diffusivity on a regular rectilinear grid."
@inline Δᶠ(i, j, k, grid::RectilinearGrid, ::TwoDimensionalLeith) = sqrt(Δxᶜᶜᶜ(i, j, k, grid) * Δyᶜᶜᶜ(i, j, k, grid))

function DiffusivityFields(grid, tracer_names, bcs, ::L2D)
    default_eddy_viscosity_bcs = (; νₑ = FieldBoundaryConditions(grid, (Center, Center, Center)))
    bcs = merge(default_eddy_viscosity_bcs, bcs)
    νₑ = CenterField(grid, boundary_conditions=bcs.νₑ)

    return (; νₑ)
end


text:
seed!(1234)

@testset "LossFunctions.jl - binary" begin
    y = categorical(["yes", "yes", "no", "yes"])
    yes, no = y[1], y[3]
    dyes = MLJBase.UnivariateFinite([yes, no], [0.6, 0.4])
    dno =  MLJBase.UnivariateFinite([yes, no], [0.3, 0.7])
    yhat = [dno, dno, dyes, dyes]
    X = nothing
    w = [1, 2, 3, 4]

    @test MLJBase.value(ZeroOneLoss(), yhat, X, y, nothing) ≈ [1, 1, 1, 0]
    @test MLJBase.value(ZeroOneLoss(), yhat, X, y, w) ≈ [1, 2, 3, 0] ./10 .* 4

    N = 10
    y = categorical(rand(["yes", "no"], N), ordered=true)
    levels!(y, ["no", "yes"])
    no, yes = MLJBase.classes(y[1])
    @test MLJBase.pm1([yes, no]) in [[+1, -1], [-1, +1]]
    ym = MLJBase.pm1(y) # observations for raw LossFunctions measure
    p_vec = rand(N) # probabilities of yes
    yhat  = map(p_vec) do p
        MLJBase.UnivariateFinite([yes, no], [p, 1 - p])
    end
    yhatm = MLJBase._scale.(p_vec) # predictions for raw LossFunctions measure
    w = rand(N)
    X = nothing

    for m in [ZeroOneLoss(), L1HingeLoss(), L2HingeLoss(), LogitMarginLoss(),
              ModifiedHuberLoss(), PerceptronLoss(), SmoothedL1HingeLoss(0.9),
              L2MarginLoss(), ExpLoss(), SigmoidLoss(), DWDMarginLoss(0.9)]
        @test MLJBase.value(m, yhat, X, y, nothing) ≈ LossFunctions.value(m, yhatm, ym)
        @test mean(MLJBase.value(m, yhat, X, y, w)) ≈
                LossFunctions.value(m, yhatm, ym, AggMode.WeightedMean(w))
    end
end

@testset "LossFunctions.jl - continuous" begin
    # losses for continuous targets:
    N    = 10
    y    = randn(N)
    yhat = randn(N)
    X    = nothing
    w    = rand(N)

    for m in [LPDistLoss(0.5), L1DistLoss(), L2DistLoss(),
              HuberLoss(0.9), EpsilonInsLoss(0.9), L1EpsilonInsLoss(0.9),
              L2EpsilonInsLoss(0.9), LogitDistLoss(), QuantileLoss(0.7)]

        @test MLJBase.value(m, yhat, X, y, nothing) ≈ LossFunctions.value(m, yhat, y)
        @test mean(MLJBase.value(m, yhat, X, y, w)) ≈
                LossFunctions.value(m, yhat, y, AggMode.WeightedMean(w))
    end
end


text:
module SparseRegressionAlgorithms

import SweepOperator
import StatsBase

typealias AVec{T} AbstractVector{T}
typealias AMat{T} AbstractMatrix{T}

include("sweep.jl")
end

max_stars_repo_name: joshday/SparseRegressionAlgorithms

max_issues_repo_name: joshday/SparseRegressionAlgorithms

max_forks_repo_name: joshday/SparseRegressionAlgorithms


text:
function vgg19_layers()
  ws = weights("vgg19.bson")
  ls = Chain(
    Conv(ws[:conv1_1_w_0][end:-1:1,:,:,:][:,end:-1:1,:,:], ws[:conv1_1_b_0], relu, pad = (1,1), stride = (1,1), dilation = (1,1)),
    Conv(ws[:conv1_2_w_0][end:-1:1,:,:,:][:,end:-1:1,:,:], ws[:conv1_2_b_0], relu, pad = (1,1), stride = (1,1), dilation = (1,1)),
    x -> maxpool(x, (2,2)),
    Conv(ws[:conv2_1_w_0][end:-1:1,:,:,:][:,end:-1:1,:,:], ws[:conv2_1_b_0], relu, pad = (1,1), stride = (1,1), dilation = (1,1)),
    Conv(ws[:conv2_2_w_0][end:-1:1,:,:,:][:,end:-1:1,:,:], ws[:conv2_2_b_0], relu, pad = (1,1), stride = (1,1), dilation = (1,1)),
    x -> maxpool(x, (2,2)),
    Conv(ws[:conv3_1_w_0][end:-1:1,:,:,:][:,end:-1:1,:,:], ws[:conv3_1_b_0], relu, pad = (1,1), stride = (1,1), dilation = (1,1)),
    Conv(ws[:conv3_2_w_0][end:-1:1,:,:,:][:,end:-1:1,:,:], ws[:conv3_2_b_0], relu, pad = (1,1), stride = (1,1), dilation = (1,1)),
    Conv(ws[:conv3_3_w_0][end:-1:1,:,:,:][:,end:-1:1,:,:], ws[:conv3_3_b_0], relu, pad = (1,1), stride = (1,1), dilation = (1,1)),
    Conv(ws[:conv3_4_w_0][end:-1:1,:,:,:][:,end:-1:1,:,:], ws[:conv3_4_b_0], relu, pad = (1,1), stride = (1,1), dilation = (1,1)),
    x -> maxpool(x, (2,2)),
    Conv(ws[:conv4_1_w_0][end:-1:1,:,:,:][:,end:-1:1,:,:], ws[:conv4_1_b_0], relu, pad = (1,1), stride = (1,1), dilation = (1,1)),
    Conv(ws[:conv4_2_w_0][end:-1:1,:,:,:][:,end:-1:1,:,:], ws[:conv4_2_b_0], relu, pad = (1,1), stride = (1,1), dilation = (1,1)),
    Conv(ws[:conv4_3_w_0][end:-1:1,:,:,:][:,end:-1:1,:,:], ws[:conv4_3_b_0], relu, pad = (1,1), stride = (1,1), dilation = (1,1)),
    Conv(ws[:conv4_4_w_0][end:-1:1,:,:,:][:,end:-1:1,:,:], ws[:conv4_4_b_0], relu, pad = (1,1), stride = (1,1), dilation = (1,1)),
    x -> maxpool(x, (2,2)),
    Conv(ws[:conv5_1_w_0][end:-1:1,:,:,:][:,end:-1:1,:,:], ws[:conv5_1_b_0], relu, pad = (1,1), stride = (1,1), dilation = (1,1)),
    Conv(ws[:conv5_2_w_0][end:-1:1,:,:,:][:,end:-1:1,:,:], ws[:conv5_2_b_0], relu, pad = (1,1), stride = (1,1), dilation = (1,1)),
    Conv(ws[:conv5_3_w_0][end:-1:1,:,:,:][:,end:-1:1,:,:], ws[:conv5_3_b_0], relu, pad = (1,1), stride = (1,1), dilation = (1,1)),
    Conv(ws[:conv5_4_w_0][end:-1:1,:,:,:][:,end:-1:1,:,:], ws[:conv5_4_b_0], relu, pad = (1,1), stride = (1,1), dilation = (1,1)),
    x -> maxpool(x, (2,2)),
    x -> reshape(x, :, size(x, 4)),
    Dense(ws[:fc6_w_0]', ws[:fc6_b_0], relu),
    Dropout(0.5f0),
    Dense(ws[:fc7_w_0]', ws[:fc7_b_0], relu),
    Dropout(0.5f0),
    Dense(ws[:fc8_w_0]', ws[:fc8_b_0]),
    softmax)
  Flux.testmode!(ls)
  return ls
end

struct VGG19 <: ClassificationModel{ImageNet.ImageNet1k}
  layers::Chain
end

VGG19() = VGG19(vgg19_layers())

Base.show(io::IO, ::VGG19) = print(io, "VGG19()")

@treelike VGG19

(m::VGG19)(x) = m.layers(x)


text:
import Base: iterate, getindex, length

abstract type AbstractNormCentralMoment end

struct stdNormCentralMoment <: AbstractNormCentralMoment
    n::Int
end

struct NormCentralMoment{T<:Real} <: AbstractNormCentralMoment
    n::Int
    σ::T
end

Base.length(iter::AbstractNormCentralMoment) = iter.n

σ(iter::stdNormCentralMoment) = 1.0
σ(iter::NormCentralMoment) = iter.σ

Base.getindex(iter::AbstractNormCentralMoment, k::Integer) = iseven(k) ? 0.0 : 2^(-k/2) * factorial(k) / factorial(k/2) * σ(iter)^k

function Base.iterate(iter::AbstractNormCentralMoment, state=(1, 1))
    # Even moment p is σᵖ(p-1)!! where `(p-1)!! = (p-1)(p-3)...(3)(1)`
    #    Note: k!! is a double factorial (product of numbers 1 to p with same parity (even/odd) as p)
    # i:  keeps track of which moment (1st, 2nd, ..., pth, etc...)
    # nf: keeps track of the double-factorial part
    i, nf = state

    # If done
    if i > iter.n
        return nothing

    # Odd moments are centered around 0
    elseif isodd(i)
        return (0.0, (i+1, nf))

    else
        nf *= (i-1)
        return (nf*σ(iter)^i, (i+1, nf))
    end
end


text:
using SeisRequests, Test
using Seis
using Dates: DateTime, Minute, Second

@testset "get_data" begin
    server = "http://service.example.com"
    SUCCESS = SeisRequests.CODE_SUCCESS
    mseed_file = joinpath(dirname(pathof(Seis)), "..", "test", "test_data", "io",
        "miniseed_GB.CWF.single_sample_gaps.mseed")
    data = read(mseed_file)

    @testset "parse_data_response" begin
        @testset "No response" begin
            request = FDSNDataSelect()
            response = HTTP.Messages.Response(request.nodata)
            @test SeisRequests.parse_data_response(Float64, request, response, server) == []

            @testset "Eltype $T" for T in (Float32, Float64)
                @test SeisRequests.parse_data_response(T, request, response, server) isa
                    Vector{Seis.Trace{T, Vector{T}, Seis.Geographic{T}}}
            end

            response.body = [1]
            @test (@test_logs (:warn, "unexpected data in response reporting no data"
                ) SeisRequests.parse_data_response(Float64, request, response, server) == [])

            request = FDSNDataSelect(nodata=404)
            response = HTTP.Messages.Response(request.nodata)
            @test SeisRequests.parse_data_response(Float64, request, response, server) == []
        end
    end

    @testset "'Success' but no data" begin
        request = FDSNDataSelect()
        response = HTTP.Messages.Response(SUCCESS)
        @test_throws ErrorException SeisRequests.parse_data_response(Float64, request,
            response, server)
    end

    @testset "Content-Type" begin
        @testset "Miniseed" begin
            request = FDSNDataSelect()

            @testset "Wrong" begin
                response_correct = HTTP.Messages.Response(SUCCESS,
                    Dict("Content-Type"=>"application/vnd.fdsn.mseed"), body=data)
                response_type_wrong = HTTP.Messages.Response(SUCCESS,
                    Dict("Content-Type"=>"WRONG"), body=data)
                @test (@test_logs (:warn, "content type of response is \"WRONG\", " *
                                          "not \"application/vnd.fdsn.mseed\" as expected"
                        ) SeisRequests.parse_data_response(Float64, request, response_type_wrong, server)) ==
                    SeisRequests.parse_data_response(Float64, request, response_correct, server)
            end

            @testset "Empty" begin
                response_correct = HTTP.Messages.Response(SUCCESS,
                    Dict("Content-Type"=>"application/vnd.fdsn.mseed"), body=data)
                response_type_empty = HTTP.Messages.Response(SUCCESS,
                    Dict(), body=data)
                @test SeisRequests.parse_data_response(Float64, request,
                        response_type_empty, server) ==
                    SeisRequests.parse_data_response(Float64, request, response_correct, server)
            end
        end
    end

    @testset "Format" begin
        @testset "Miniseed" begin
            request = FDSNDataSelect()
            @testset "Eltype $T" for T in (Float32, Float64)
                response = HTTP.Messages.Response(SUCCESS,
                    Dict("Content-Type"=>"application/vnd.fdsn.mseed"), body=data)
                out = SeisRequests.parse_data_response(T, request, response, server)
                @test out isa Vector{Seis.Trace{T, Vector{T}, Seis.Geographic{T}}}
                @test length(out) == 2
                @test channel_code(out[1]) == "GB.CWF..BHZ"
                @test channel_code(out[2]) == "GB.CWF..HHZ"
                @test all(out.meta.server .== server)
                t = Seis.read_mseed(Seis.Trace{T, Vector{T}, Seis.Geographic{T}},
                    mseed_file)
                t.meta.mseed_file = missing
                t.meta.server = server
                @test out == t
            end
        end
    end

    @testset "Real requests" begin
        @testset "Verbosity" begin
            @testset "No data" begin
                @test (@test_logs (:info, "Request status: Request was properly " *
                    "formatted and submitted but no data matches the selection"
                    ) get_data(code="XX.XXX.XX.BHZ", starttime="2000-01-01",
                               endtime="2000-01-02")) == []
                @test_logs get_data(code="XX.XXX.XX.BHZ", starttime="2000-01-01",
                               endtime="2000-01-02", verbose=false)
            end
            @testset "Data" begin
                @test_logs (:info, "Request status: Successful request, results follow"
                    ) get_data(code="GB.JSA..BHZ",
                    starttime="2015-01-01", endtime="2015-01-01T00:00:01")
            end
        end
    end

    @testset "Simple" begin
        code = "IU.ANMO..BHZ"
        stime = DateTime(1990)
        etime = stime + Second(1)
        @testset "Eltype $T" for T in (Float32, Float64)
            t = get_data(code=code, starttime=stime, endtime=etime,
                T=T, verbose=false)
            @test t isa Vector{Seis.Trace{T, Vector{T}, Seis.Geographic{T}}}
        end
        @testset "Default eltype" begin
            t = get_data(code=code, starttime=stime, endtime=etime, verbose=false)
            @test t isa Vector{Seis.Trace{Float64, Vector{Float64}, Seis.Geographic{Float64}}}
            # These should never change, but rely on the IRIS data centre
            @test length(t) == 1
            @test all(x -> channel_code(x) == "IU.ANMO..BHZ", t)
            # Should be 1 second of data
            @test nsamples(first(t)) ≈ 1/first(t).delta
        end
    end

    @testset "Station" begin
        @testset "Trace eltype $Tt" for Tt in (Float32, Float64)
            @testset "Station eltype $Ts" for Ts in (Float32, Float64)
                stas = [Station{Ts}(net="G", sta="HYB", loc="", cha="BHZ",
                                    lon=78.553, lat=17.417, elev=510, azi=0, inc=0),
                        Station{Ts}(net="G", sta="SCZ", loc="", cha="BHZ",
                                    lon=-121.403, lat=36.598, elev=261, azi=0, inc=0)]
                t = get_data(stas, DateTime(2001), DateTime(2001)+Second(1), T=Tt,
                    verbose=false)
                @test t isa Vector{Seis.Trace{Tt,Vector{Tt},Seis.Geographic{Tt}}}
                # Assume data is returned; will fail if datacentre fails
                if isempty(t)
                    @warn("No data returned for Trace $Tt Station $Ts")
                else
                    @test all(x -> channel_code(x) in channel_code.(stas), t)
                    for tt in t
                        @test tt.sta in convert.(GeogStation{Tt}, stas)
                    end
                end
            end
        end
    end

    @testset "Event and stations" begin
        types = (Float32,  Float64)
        @testset "Trace eltype $Tt" for Tt in types
            @testset "Station eltype $Ts" for Ts in types
                stas = [Station{Ts}(net="IU", sta="FURI", loc="00", cha=comp,
                                    lon=-62.35, lat=82.5033, elev=60, azi=azi, inc=inc)
                        for (comp, azi, inc) in zip(
                            ("BHE", "BHN", "BHZ"), (90, 0, 0), (90, 90, 0))]
                @testset "Event eltype $Te" for Te in types
                    evt = Event{Te}(lon=-88.729, lat=-12.997, dep=82.9,
                                    time=DateTime("2001-01-13T17:33:34.58"))
                    t = get_data(evt, stas, Second(5), Minute(1) + Second(1),
                        verbose=false, T=Tt)
                    @test t isa Vector{Trace{Tt,Vector{Tt},Seis.Geographic{Tt}}}
                    if isempty(t)
                        @warn("No data returned for Event $Te Station $Ts Trace $Tt")
                    else
                        # Ensure correction conversion so that comparison with
                        # trace fields isn't wrong due to loss of precision
                        mintype = any((Tt, Ts, Tt) .== Float32) ? Float32 : Float64
                        for tt in t
                            @test channel_code(tt) in channel_code.(stas)
                            @test convert(GeogStation{mintype}, tt.sta) in
                                convert.(GeogStation{mintype}, stas)
                            @test convert(GeogEvent{mintype}, tt.evt) ==
                                convert(GeogEvent{mintype}, evt)
                            @test starttime(tt) ≈ 5 atol=0.1
                            @test endtime(tt) ≈ 61 atol=0.1
                        end
                    end
                    t′ = get_data(evt, stas, 5, 61, verbose=false, T=Tt)
                    @test t == t′
                end
            end
        end
    end
end


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

using Actors, Test

struct MyMsg <: Msg
    x::Int
end

struct MySource <: Msg
    x::Int
    from::Link
end

function writeLk(lk::Link, src=true, slp=false)
    for i in 1:4
        put!(lk.chn, MyMsg(i))
        yield()
        slp && sleep(rand()*0.1)
    end
    if src
        put!(lk.chn, MySource(10, lk2))
        yield()
        slp && sleep(rand()*0.1)
    end
    for i in 5:8
        put!(lk.chn, MyMsg(i))
        yield()
        slp && sleep(rand()*0.1)
    end
end

function readLk(lk::Link)
    buf = Int[]
    while isready(lk.chn)
        push!(buf, take!(lk.chn).x)
    end
    buf
end

lk1 = Actors.newLink()
lk2 = Actors.newLink()

@test Actors._match(MyMsg(1), nothing, nothing)
@test Actors._match(MySource(1,lk2), nothing, nothing)

@test Actors._match(MyMsg(1), MyMsg, nothing)
@test !Actors._match(MyMsg(1), Response, nothing)

@test !Actors._match(MyMsg(1), nothing, lk2)
@test Actors._match(MySource(1,lk2), nothing, lk2)
@test !Actors._match(MySource(1,lk2), nothing, lk1)

@test !Actors._match(MyMsg(1), MyMsg, lk2)
@test !Actors._match(MySource(1,lk2), MyMsg, lk2)
@test Actors._match(MySource(1,lk2), MySource, lk2)
@test !Actors._match(MySource(1,lk2), MySource, lk1)

writeLk(lk1)
@test length(lk1.chn.data) == 9
msg = receive(lk1, MySource, lk2)
@test msg == MySource(10, lk2)
@test length(lk1.chn.data) == 8
@test readLk(lk1) == collect(1:8)
@test length(lk1.chn.data) == 0

msg = receive(lk1, MySource, lk2, timeout=1)
@test msg == Actors.Timeout()

writeLk(lk1, false)
@test length(lk1.chn.data) == 8
msg = receive(lk1, MySource, lk2, timeout=0)
@test msg == Actors.Timeout()
@test length(lk1.chn.data) == 8
@test readLk(lk1) == collect(1:8)
@test length(lk1.chn.data) == 0

@async writeLk(lk1, true, true)
msg = receive(lk1, MySource, lk2)
@test msg == MySource(10, lk2)
sleep(1)
@test length(lk1.chn.data) == 8
@test readLk(lk1) == collect(1:8)
@test length(lk1.chn.data) == 0

comtest(msg::MySource) = nothing
function comtest(msg::Request)
    res = (1, (2), (3,4,5), [6,7,8])
    if msg.x in 1:4
        send(msg.from, Response(res[msg.x], self()))
    else
        send(msg.from, Response("test", self()))
    end
    return nothing
end

A = Actors.spawn(Bhv(comtest))
send(A, MySource(1, lk2))
msg = receive(lk2, timeout=1)
@test msg == Actors.Timeout()
send(A, Request(10, lk2))
msg = receive(lk2, timeout=1)
@test msg == Response("test", A)

res = request(A, MySource(1, lk2), timeout=1)
@test res == Actors.Timeout()
res = request(A, Request(1, lk2), full=true, timeout=1)
@test res == Response(1, A)
@test request(A, Request(1, lk2), timeout=1) == 1
@test request(A, Request(2, lk2), timeout=1) == 2
@test request(A, Request(3, lk2), timeout=1) == (3,4,5)
@test request(A, Request(4, lk2), timeout=1) == [6,7,8]
@test request(A, Request(99, lk2), timeout=1) == "test"

# test send_after
counter = Ref(0)
tick(c, t, n) = n > 0 && (c[]+=1; send_after(self(), t, c, t, n-1))
ticker = Actors.spawn(tick)
send_after(ticker, 0.1, counter, 0.1, 5)
sleep(0.6)
@test counter[] == 5
@test Actors.diag(ticker) == :ok
exit!(ticker)


text:
"""
    write_metadata_NDE_training(FILE_PATH, train_files, train_epochs, train_tranges, opts, uw_NN, vw_NN, wT_NN)
"""
function write_metadata_NDE_training(FILE_PATH, train_files, train_epochs, train_tranges, train_parameters, opts, uw_NN, vw_NN, wT_NN)
    jldopen(FILE_PATH, "w") do file
        training_info = JLD2.Group(file, "training_info")
        training_info["train_files"] = train_files
        training_info["train_epochs"] = train_epochs
        training_info["train_tranges"] = train_tranges
        training_info["optimizers"] = opts
        training_info["parameters"] = train_parameters
        training_info["uw_neural_network"] = uw_NN
        training_info["vw_neural_network"] = vw_NN
        training_info["wT_neural_network"] = wT_NN

        training_data = JLD2.Group(file, "training_data")
        loss = JLD2.Group(training_data, "loss")
        neural_network = JLD2.Group(training_data, "neural_network")
        uw = JLD2.Group(neural_network, "uw")
        vw = JLD2.Group(neural_network, "vw")
        wT = JLD2.Group(neural_network, "wT")
        η = JLD2.Group(training_data, "η")
        β = JLD2.Group(training_data, "β")
        state = JLD2.Group(training_data, "state")
    end
end

function write_data_NDE_training(FILE_PATH, losses, loss_scalings, uw_NN, vw_NN, wT_NN, stage, optimizer)
    jldopen(FILE_PATH, "a") do file
        profile_loss = losses.u + losses.v + losses.T
        gradient_loss = losses.∂u∂z + losses.∂v∂z + losses.∂T∂z
        total_loss = profile_loss + gradient_loss

        if !haskey(file, "training_data/loss/total/$stage")
            file["training_info/loss_scalings"] = loss_scalings
            file["training_data/loss/total/$stage/1"] = total_loss
            file["training_data/loss/profile/$stage/1"] = profile_loss
            file["training_data/loss/gradient/$stage/1"] = gradient_loss

            file["training_data/loss/u/$stage/1"] = losses.u
            file["training_data/loss/v/$stage/1"] = losses.v
            file["training_data/loss/T/$stage/1"] = losses.T

            file["training_data/loss/∂u∂z/$stage/1"] = losses.∂u∂z
            file["training_data/loss/∂v∂z/$stage/1"] = losses.∂v∂z
            file["training_data/loss/∂T∂z/$stage/1"] = losses.∂T∂z

            file["training_data/neural_network/uw/$stage/1"] = uw_NN
            file["training_data/neural_network/vw/$stage/1"] = vw_NN
            file["training_data/neural_network/wT/$stage/1"] = wT_NN

            file["training_data/optimizer/η/$stage/1"] = optimizer.eta
            file["training_data/optimizer/β/$stage/1"] = optimizer.beta
            file["training_data/optimizer/state/$stage/1"] = optimizer.state
        else
            count = length(keys(file["training_data/loss/total/$stage"])) + 1
            file["training_data/loss/total/$stage/$count"] = total_loss
            file["training_data/loss/profile/$stage/$count"] = profile_loss
            file["training_data/loss/gradient/$stage/$count"] = gradient_loss

            file["training_data/loss/u/$stage/$count"] = losses.u
            file["training_data/loss/v/$stage/$count"] = losses.v
            file["training_data/loss/T/$stage/$count"] = losses.T

            file["training_data/loss/∂u∂z/$stage/$count"] = losses.∂u∂z
            file["training_data/loss/∂v∂z/$stage/$count"] = losses.∂v∂z
            file["training_data/loss/∂T∂z/$stage/$count"] = losses.∂T∂z

            file["training_data/neural_network/uw/$stage/$count"] = uw_NN
            file["training_data/neural_network/vw/$stage/$count"] = vw_NN
            file["training_data/neural_network/wT/$stage/$count"] = wT_NN

            file["training_data/optimizer/η/$stage/$count"] = optimizer.eta
            file["training_data/optimizer/β/$stage/$count"] = optimizer.beta
            file["training_data/optimizer/state/$stage/$count"] = optimizer.state
        end
    end
end

function write_metadata_NN_training(FILE_PATH, train_files, train_parameters, train_epochs, opts, NN, NN_type)
    jldopen(FILE_PATH, "w") do file
        training_info = JLD2.Group(file, "training_info")
        training_info["train_files"] = train_files
        training_info["train_epochs"] = train_epochs
        training_info["optimizers"] = opts
        training_info["parameters"] = train_parameters
        training_info["$(NN_type)_neural_network"] = NN
        training_data = JLD2.Group(file, "training_data")
    end
end

function write_data_NN_training(FILE_PATH, loss, NN)
    jldopen(FILE_PATH, "a") do file
        if !haskey(file, "training_data/loss")
            file["training_data/loss/1"] = loss
        else
            count = length(keys(file["training_data/loss"])) + 1
            file["training_data/loss/$count"] = loss
        end

        if !haskey(file, "training_data/neural_network")
            file["training_data/neural_network/1"] = NN
        else
            count = length(keys(file["training_data/neural_network"])) + 1
            file["training_data/neural_network/$count"] = NN
        end

    end
end

function write_data_NN(FILE_PATH, uw_NN, vw_NN, wT_NN)
    jldopen(FILE_PATH, "w") do file
        file["neural_network/uw"] = uw_NN
        file["neural_network/vw"] = vw_NN
        file["neural_network/wT"] = wT_NN
    end
end

function write_metadata_modified_pacanowski_philander_optimisation(FILE_PATH, train_files, train_epochs, train_tranges, train_parameters, loss_scalings, opts)
    jldopen(FILE_PATH, "w") do file
        training_info = JLD2.Group(file, "training_info")
        training_info["train_files"] = train_files
        training_info["train_epochs"] = train_epochs
        training_info["train_tranges"] = train_tranges
        training_info["optimizers"] = opts
        training_info["parameters"] = train_parameters
        training_info["loss_scalings"] = loss_scalings

        training_data = JLD2.Group(file, "training_data")
        loss = JLD2.Group(training_data, "loss")
        parameters = JLD2.Group(training_data, "parameters")
    end
end

function write_data_modified_pacanowski_philander_optimisation(FILE_PATH, losses, parameters)
    profile_loss = losses.u + losses.v + losses.T
    gradient_loss = losses.∂u∂z + losses.∂v∂z + losses.∂T∂z
    total_loss = profile_loss + gradient_loss

    jldopen(FILE_PATH, "a") do file
        if !haskey(file, "training_data/loss/total")
            file["training_data/loss/total/1"] = total_loss
            file["training_data/loss/profile/1"] = profile_loss
            file["training_data/loss/gradient/1"] = gradient_loss

            file["training_data/loss/u/1"] = losses.u
            file["training_data/loss/v/1"] = losses.v
            file["training_data/loss/T/1"] = losses.T

            file["training_data/loss/∂u∂z/1"] = losses.∂u∂z
            file["training_data/loss/∂v∂z/1"] = losses.∂v∂z
            file["training_data/loss/∂T∂z/1"] = losses.∂T∂z

            file["training_data/parameters/1"] = parameters
        else
            count = length(keys(file["training_data/loss/total"])) + 1

            file["training_data/loss/total/$count"] = total_loss
            file["training_data/loss/profile/$count"] = profile_loss
            file["training_data/loss/gradient/$count"] = gradient_loss

            file["training_data/loss/u/$count"] = losses.u
            file["training_data/loss/v/$count"] = losses.v
            file["training_data/loss/T/$count"] = losses.T

            file["training_data/loss/∂u∂z/$count"] = losses.∂u∂z
            file["training_data/loss/∂v∂z/$count"] = losses.∂v∂z
            file["training_data/loss/∂T∂z/$count"] = losses.∂T∂z

            file["training_data/parameters/$count"] = parameters
        end
    end
end


text:
using Base.Test
include("src/Robots.jl")

test1 = [1 1 ; 2 3]
test2 = [0 0 ; 0 0]
scrap_robots!(test1,test2)
@test test1 == [1 1 ; 0 0]
@test test2 == [0 0 ; 1 1]

a = [1 1 ; 2 2]
set_indices_to!(a, [1 1;2 2],3)
@test a == [3 1; 2 3]
a = [1 1 ; 2 2]
set_indices_to!(a, [1 1],3)
@test a == [3 1; 2 2]

@test towards(5,5) == 5
@test towards(5,4) == 4
@test towards(8,2) == 7

g = Game(Human())

@test has_robots(g.board)

win = TermWin.initscr()
TermWin.noecho()
TermWin.curs_set(0)

play(g,g.player,win);

TermWin.endwin()

println("You made it to level ",g.level, " with a score of ",score(g),"!")

max_stars_repo_name: RatanRSur/robots-from-scratch

max_issues_repo_name: RatanRSur/robots-from-scratch

max_forks_repo_name: RatanRSur/robots-from-scratch


text:
module Preconditioners

using JuLIP, LinearAlgebra

@doc raw"""
`AdjacancyPrecon` :

WARNING: While the code accounts for PBC, the following description does not.

If `(r0, r1, k)` is in the `intervals` list and if
`r0 <= rij <= r1` then the 3 x 3 block
```math
   P_{ij} = - k \big( (1-c) \hat{R}_{ij} \otimes \hat{R}_{ij} + c I \big),
```
where ``\hat{R}_{ij}`` is the bond direction and `c = innerstab`. The diagonal
is then obtained by
```math
   P_{ii} = \sum_{j \neq i} P_{ij} + c_g I
```
where `cg = stab` is a second, global, stabilisation constant.

Regarding the choice of `innerstab`:
* Choose `innerstab = 0` to get full infinitesimal rotation-invariance
* Choose `innerstab = 1` to get full resistance against infinitesimal rotations
* In optimisation we often found that `innerstab = 0.1` is a decent choice.
"""
struct AdjacancyPrecon
   intervals::Vector{NamedTuple{(:r0, :r1, :k), Tuple{Float64, Float64, Float64}}}
   stab::Float64
   innerstab::Float64
end

import JuLIP.FIO: read_dict, write_dict

AdjacancyPrecon(intervals; stab = 1e-4, innerstab = 0.1) =
   AdjacancyPrecon(intervals, stab, innerstab)

_cutoff(P::AdjacancyPrecon) =
   maximum( i.r1 for i in P.intervals )

function (P::AdjacancyPrecon)(R::JVecF)
   r = norm(R)
   R̂ = R / r
   Π = (1 - P.innerstab) * R̂ * R̂' + P.innerstab * one(JMat{Float64})
   for i in P.intervals
      if i.r0 <= r <= i.r1
         return i.k * Π
      end
   end
   return 0.0 * Π
end

function (precon::AdjacancyPrecon)(at::Atoms)
   P = zeros(3*length(at), 3*length(at))
   nlist = neighbourlist(at, _cutoff(precon))
   for i = 1:length(at)
      Js, Rs = neigs(nlist, i)
      for (j, Rij) in zip(Js, Rs)
         if j < i; continue; end
         Pij = precon(Rij)
         inds = 3*(i-1) .+ (1:3)
         jnds = 3*(j-1) .+ (1:3)
         for a = 1:3, b = 1:3
            P[inds[a], jnds[b]] -= Pij[a,b]
            P[jnds[b], inds[a]] -= Pij[a,b]
            P[inds[a], inds[b]] += Pij[a,b]
            P[jnds[a], jnds[b]] += Pij[a,b]
         end
      end
   end
   P += precon.stab * I
   return P
end

end


text:
@testset "189.rotate-array.jl" begin
    @test rotate_arr([1,2,3,4,5,6,7], 3) == [5,6,7,1,2,3,4]
    @test rotate_arr([-1,-100,3,99], 2) == [3,99,-1,-100]
end


text:
export get_sources

"""
    get_sources(con::SeisCon)

Returns an array of the source location coordinate pairs, NOT the minimum and maximum values.
Unlike `get_headers`, `get_sources` checks to make sure that SourceX and SourceY are consistant
throughout each block, that is `min == max`.

Column 1 of the returned array is SourceX, and Column 2 is SourceY.

# Example

sx = get_sources(s)

"""
function get_sources(con::SeisCon)

    sx_v = get_header(con, "SourceX")
    sy_v = get_header(con, "SourceY")

    sx_const = all([sx_v[i,1] == sx_v[i,2] for i in 1:size(sx_v)[1]])
    sy_const = all([sy_v[i,1] == sy_v[i,2] for i in 1:size(sy_v)[1]])

    if sx_const && sy_const
        sx = [sx_v[i,1] for i in 1:size(sx_v)[1]]
        sy = [sy_v[i,1] for i in 1:size(sy_v)[1]]
    else
        @error "Source locations are not constistant for all blocks."
    end

    return hcat(sx,sy)

end


text:
# # # Code:
# using PyPlot
# include("types.jl")
# #
# # colors = Dict(
# #                "red"    =>  "r",
# #                "blue"   =>  "b",
# #                "black"  =>  "k",
# #              )
# #
# function plot_frame(pose::Pose2D, color="black")
#
#    plt[:figure]()
#    ax = plt[:gca]()
#
#    #ax[:arrow](x,y,dx,dy, fc="k", ec="k", head_width=0.1, head_length=0.1 )
#    ax[:arrow](pose.p.x, pose.p.y, pose.ξ[1,1], pose.ξ[2,1], fc=colors[color], ec=colors[color], head_width=0.1, head_length=0.1)
#    ax[:arrow](pose.p.x, pose.p.y, pose.ξ[1,2], pose.ξ[2,2], fc=colors[color], ec=colors[color], head_width=0.1, head_length=0.1)
#
#    # TODO(elsuizo): ver para que quede centrado alrededor del punto p con un offset
#    ax[:set_xlim]([-5 - pose.p.x, 5 + pose.p.x])
#    xlabel(L"x")
#    ax[:set_ylim]([-5 - pose.p.x, 5 + pose.p.x])
#    ylabel(L"y")
#
#    title("Frame plot")
#    grid("on")
#    plt[:draw]()
#    plt[:show]()
# end
# #
# # function plot_frame(pose_vector::Array{Pose2D,1})
# #
# #    plt[:figure]()
# #    ax = plt[:gca]()
# #    color = "black"
# #
# #    x_position_sum = 0
# #    y_position_sum = 0
# #    number_of_pose = length(pose_vector)
# #    #ax[:arrow](x,y,dx,dy, fc="k", ec="k", head_width=0.1, head_length=0.1 )
# #    for pose in pose_vector
# #       ax[:arrow](pose.p.x, pose.p.y, pose.ξ[1,1], pose.ξ[2,1], fc=colors[color], ec=colors[color], head_width=0.1, head_length=0.1)
# #       ax[:arrow](pose.p.x, pose.p.y, pose.ξ[1,2], pose.ξ[2,2], fc=colors[color], ec=colors[color], head_width=0.1, head_length=0.1)
# #    end
# #
# #    for pose in pose_vector
# #       x_position_sum += pose.p.x
# #       y_position_sum += pose.p.y
# #    end
# #    x_mean = x_position_sum / number_of_pose
# #    y_mean = x_position_sum / number_of_pose
# #
# #    # TODO(elsuizo): ver para que quede centrado alrededor del punto p con un offset
# #    ax[:set_xlim]([-5 - x_mean, 5 + x_mean])
# #    xlabel(L"x")
# #    ax[:set_ylim]([-5 - y_mean, 5 + y_mean])
# #    ylabel(L"y")
# #
# #
# #    title("Frame plots")
# #    grid("on")
# #    plt[:draw]()
# #    plt[:show]()
# #
# # end
# #
# # #
# # # graphics.jl ends here

max_stars_repo_name: elsuizo/Robot

max_issues_repo_name: elsuizo/Robot

max_forks_repo_name: elsuizo/Robot


text:
include("../interfaces.jl")

import PhysicalConstants.CODATA2018 as Constants

### Boundary conditions
Temperature = Model(
    T = 293.15u"K",
    port = HeatPort,
    equations = :[port.T=T]
)

# Base Classes
PartialResistance = Model(
    port_a = HeatPort,
    port_b = HeatPort,
    partialEquations  = :[
        0 = port_a.Q_flow + port_b.Q_flow
        port_a.Q_flow = (port_a.T - port_b.T)/R
    ]
)

### Resistance Models
Sphere = PartialResistance | Model(
    r_inner = 0.1u"m",
    r_outer = 1.0u"m",
    lambda = 1.0u"W/(m*K)",
    R = Var(start=1.0u"K/W"),
    π=π,
    equations = :[
        R = 1/(4*π*lambda)*(1/r_inner - 1/r_outer)
    ]
)

Convection = PartialResistance | Model(
    surfaceArea = 0.01u"m^2",
    alpha = 1000.0u"W/(m^2*K)",
    R = Var(start=1.0u"K/W"),
    equations = :[
        R = 1/(alpha*max(eps()u"m^2",surfaceArea))
    ]
)

Contact = PartialResistance | Model(
    surfaceArea = 0.01u"m^2",
    Rc_pp = 1.0u"m^2*K/W",
    R = Var(start=1.0u"K/W"),
    equations = :[
        R = Rc_pp/max(eps()u"m^2",surfaceArea)
    ]
)

Radiation = PartialResistance | Model(
    surfaceArea = 0.01u"m^2",
    epsilon = 1.0u"m/m", # no units... don't know how to represent that
    R = Var(start=1.0u"K/W"),
    useExact=true,
    sigma=Constants.StefanBoltzmannConstant, # how to use this within a model without specifying as a parameter/variable?
    equations = :[
        R = if useExact
                1/(surfaceArea*sigma*epsilon*(port_a.T^2+port_b.T^2)*(port_a.T + port_b.T))
            else
                1/(4*surfaceArea*sigma*epsilon*(0.5*(port_a.T + port_b.T))^3)
            end
    ]
)


text:
#=
Retieves data to replicate Chiu et al (2011) and Ghysels (2016)

From 1947 to 2011:

* Real GDP, quarterly (BEA)
* Industrial Production, monthly (FRED)
* CPI Inflation, monthly (BLS)
* Unemployment rate, monthly (BLS)

=#
using BeaData, FredData
using DataFrames

#--------- BEA Data
bea = Bea()
startyear = 1947
endyear = 2011
# Real GDP
nipa116 = get_nipa_table(bea, 6, "Q", startyear, endyear)
qdata = nipa116.df[[:date, :line1]]
rename!(qdata, :line1, :GDPlev)

#--------- FRED Data
fred = Fred()
startdate = "1947-01-01"
enddate = "2011-12-31"
# Unemployment rate
UNEMP_id = "UNRATE"
unemp_data = FredData.get_data(fred, UNEMP_id, observation_start = startdate, observation_end = enddate)
UNEMP = unemp_data.df[[:date, :value]]
rename!(UNEMP, :value, :UNEMP)
# CPI
CPI_id = "CPIAUCSL"
cpi_data = FredData.get_data(fred, CPI_id, observation_start = startdate, observation_end = enddate)
CPI = cpi_data.df[[:date, :value]]
rename!(CPI, :value, :CPI)
# Industrial Production
IP_id = "INDPRO"
ip_data = FredData.get_data(fred, IP_id, observation_start = startdate, observation_end = enddate)
INDPRO = ip_data.df[[:date, :value]]
rename!(INDPRO, :value, :INDPRO)

mdata = join(INDPRO, CPI,  on=:date, kind=:left)
mdata = join(mdata, UNEMP, on=:date, kind=:left)
sort!(mdata, cols = :date)

# Save Data
writetable("original_data/mdata.csv", mdata)
writetable("original_data/qdata.csv", qdata)

max_stars_repo_name: lnsongxf/MF-VAR

max_issues_repo_name: lnsongxf/MF-VAR

max_forks_repo_name: lnsongxf/MF-VAR