text:
# See run.jl for command-line example runs.
export run_trials, cont_var_result, print_cont_var_result, run_trial, writeheader, writerow, my_isdefined, check_parameters

@doc """ function run_trials(sim_record::ContVarEvolution.cont_var_result_type )
  Runs multiple trials of the simulation using the parameter file "\$(simname).jl" where  simname is the first command line argument.
  Trials can be run in parallel using the Julia parallel map (pmap) facility.
"""
#function run_trials( simname::AbstractString, sim_record::ContVarEvolution.cont_var_result_type )
function run_trials( simname::AbstractString, sim_record::cont_var_result_type )
  stream = open("$(simname).csv","w")
  println("stream: ",stream)
  # collect these parameter records into the list   sim_record_list_run
  sim_record_list_run = cont_var_result_type[]
  trial=1
  #if my_isdefined(:mutation_stddev_list)
  if length(sim_record.mutation_stddev_list) > 0
    for N in sim_record.N_list
      for num_attributes in sim_record.num_attributes_list
        for mutation_stddev in sim_record.mutation_stddev_list
          for trial = 1:sim_record.num_trials
            sim_record_run = deepcopy( sim_record )
            sim_record_run.N = N
            sim_record_run.num_attributes = num_attributes
            sim_record_run.mutation_stddev = mutation_stddev
            Base.push!(sim_record_list_run, sim_record_run )
          end
        end
      end
    end
  #elseif my_isdefined(:N_mut_list)
  elseif length(sim_record.N_mut_list) > 0
    for N in sim_record.N_list
      for num_attributes in sim_record.num_attributes_list
        for N_mut in sim_record.N_mut_list
          for trial = 1:sim_record.num_trials
            sim_record_run = deepcopy( sim_record )
            sim_record_run.N = N
            sim_record_run.num_attributes = num_attributes
            sim_record_run.mutation_stddev = N_mut/N
            #println("N:",sim_record_run.N,"  num_attr:",sim_record_run.num_attributes,"  mutation_stddev: ",sim_record_run.mutation_stddev)
            Base.push!(sim_record_list_run, sim_record_run )
          end
        end
      end
    end
  else
    error("Either mutation_stddev_list or N_mut_list must be defined in the configuration file.")
  end
  println("===================================")
  # Run the simulation function "cont_var_simulation" on each parameter record in parallel
  # For each trial, "cont_var_simulation" adds the trial results to the parameter/result record
  # You may want to change "pmap" to "map" for debugging purposes if you are getting complicated error messages.
  sim_record_list_result = pmap(cont_var_simulation, sim_record_list_run )
  #sim_record_list_result = map(cont_var_simulation, sim_record_list_run )
  trial = 1
  writeheader( stdout, sim_record )
  writeheader( stream, sim_record )
  for sim_record_result in sim_record_list_result    # write 1 record per trial
    writerow(stdout,trial,sim_record_result)
    writerow(stream,trial,sim_record_result)
    trial += 1
  end
end

function cont_var_result( N_list::Vector{Int64}, num_attributes_list::Vector{Int64}, mutation_stddev_list::Vector{Float64}, N_mut_list::Vector{Float64},
        num_trials, N::Int64, num_subpops::Int64, num_attributes::Int64, ngens::Int64, burn_in::Number,
        mutation_stddev::Float64, ideal::Float64, fit_slope::Float64, neutral::Bool=false )
  if typeof(burn_in) == Int64
    int_burn_in = burn_in
  else
    int_burn_in = Int(round(burn_in*N+50.0))
  end
  return cont_var_result_type( N_list, num_attributes_list, mutation_stddev_list, N_mut_list, num_trials, N, num_subpops, num_attributes, ngens, int_burn_in,
      mutation_stddev, ideal, fit_slope, neutral, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0,0,0,0 )
end

@doc """ function writeheader()
  Writes the parameters and the header line for the output CSV file.
  The parameters that do not vary by trial are first written as comment lines where the comment character is "#".
  Then the header line is written.  The headers of course must correspond to the values written by the writerows function.
"""
function writeheader( stream::IO, sim_record::cont_var_result_type )
  param_strings = [
    "# $(string(Dates.today()))",
    "# N_list=$(sim_record.N_list)",
    "# num_attributes_list=$(sim_record.num_attributes_list)",
    "# mutation_stddev_list=$(sim_record.mutation_stddev_list)",
    "# N_mut_list=$(sim_record.N_mut_list)",
    "# num_trials=$(sim_record.num_trials)",
    #"# num_attributes=$(sim_record.num_attributes)",
    "# ngens=$(sim_record.ngens)",
    "# neutral=$(sim_record.neutral)",
    "# ideal=$(sim_record.ideal)",
    "# fit_slope=$(sim_record.fit_slope)"]

  write(stream,join(param_strings,"\n"),"\n")
  heads = [
    "N",
    "N_mut",
    "mutation_stddev",
    "num_attributes",
    "int_burn_in",
    "attribute_mean",
    "attribute_median",
    "attribute_coef_var"
  ]
  fitness_heads = [   # not written if the simulation is neutral
    "fitness_mean",
    "fitness_median",
    "fitness_coef_var",
    "fit_diff_neg_fract",
    "fit_diff_neg_neutral",
    "fit_diff_pos_neutral",
    "fit_diff_pos_fract"
  ]
  if !sim_record.neutral
    heads = vcat(heads,fitness_heads)
  end
  write(stream,join(heads,","),"\n")
end

@doc """ function writerow()
  Writes one row of the CSV file, where a row corresponds to 1 trial of the simulation.
"""
function writerow( stream::IO, trial::Int64, sim_record::cont_var_result_type )
  sum_fitdiff = Float64(sum( (sim_record.neg_count, sim_record.neg_neutral, sim_record.pos_neutral, sim_record.pos_count) ))
  values = Any[
          sim_record.N,
          sim_record.N*sim_record.mutation_stddev,
          sim_record.mutation_stddev,
          sim_record.num_attributes,
          sim_record.int_burn_in,
          sim_record.attribute_mean,
          sim_record.attribute_median,
          sim_record.attribute_coef_var
  ]
  fitness_values = [
          sim_record.fitness_mean,
          sim_record.fitness_median,
          sim_record.fitness_coef_var,
          sim_record.neg_count/sum_fitdiff,
          sim_record.neg_neutral/sum_fitdiff,
          sim_record.pos_neutral/sum_fitdiff,
          sim_record.pos_count/sum_fitdiff
  ]
  if !sim_record.neutral
    values = vcat(values,fitness_values)
  end
  write(stream,join(values,","),"\n")
end


text:
struct LazyMulAdd{M,O,T<:NativeTypesV} <: Number
    data::T
    @inline LazyMulAdd{M,O,T}(data::T) where {M,O,T<:NativeTypesV} = new{M,O,T}(data)
end
# O for offset is kind of hard to read next to default of 0?
@inline LazyMulAdd{M,O}(data::T) where {M,O,T<:Union{Base.HWReal,AbstractSIMD}} = LazyMulAdd{M,O,T}(data)
@inline LazyMulAdd{M}(data::T) where {M,T<:NativeTypesV} = LazyMulAdd{M,0,T}(data)
@inline LazyMulAdd{0,O}(data::T) where {O,T<:Union{Base.HWReal,AbstractSIMD}} = @assert false#StaticInt{O}()
@inline LazyMulAdd{0}(data::T) where {T<:NativeTypesV} = @assert false#StaticInt{0}()
@inline LazyMulAdd(data::T, ::StaticInt{M}) where {M,T} = LazyMulAdd{M,0,T}(data)
@inline LazyMulAdd(data::T, ::StaticInt{M}, ::StaticInt{O}) where {M,O,T} = LazyMulAdd{M,O,T}(data)

@inline data(lm::LazyMulAdd) = data(getfield(lm, :data)) # calls data on inner for use with indexing (normally `data` only goes through one layer)

@inline _materialize(a::LazyMulAdd{M,O,I}) where {M,O,I} = vadd_nsw(vmul_nsw(StaticInt{M}(), getfield(a, :data)), StaticInt{O}())
@inline _materialize(x) = x
@inline Base.convert(::Type{T}, a::LazyMulAdd{M,O,I}) where {M,O,I,T<:Number} = convert(T, _materialize(a))
@inline Base.convert(::Type{LazyMulAdd{M,O,I}}, a::LazyMulAdd{M,O,I}) where {M,O,I} = a
# @inline Base.convert(::Type{LazyMulAdd{M,O,I}}, a::LazyMulAdd{M}) where {M,O,I} = a
# @inline Base.convert(::Type{LazyMulAdd{M,T,I}}, a::LazyMulAdd{M,StaticInt{O},I}) where {M,O,I,T} = a

Base.promote_rule(::Type{LazyMulAdd{M,O,I}}, ::Type{T}) where {M,O,I<:Number,T} = promote_type(I, T)
Base.promote_rule(::Type{LazyMulAdd{M,O,MM{W,X,I}}}, ::Type{T}) where {M,O,W,X,I,T} = promote_type(MM{W,X,I}, T)
Base.promote_rule(::Type{LazyMulAdd{M,O,Vec{W,I}}}, ::Type{T}) where {M,O,W,I,T} = promote_type(Vec{W,I}, T)

@inline lazymul(a, b) = vmul_nsw(a, b)
@inline lazymul(::StaticInt{M}, b) where {M} = LazyMulAdd{M}(b)
@inline lazymul(::StaticInt{1}, b) = b
@inline lazymul(::StaticInt{0}, b) = StaticInt{0}()
@inline lazymul(a, ::StaticInt{M}) where {M} = LazyMulAdd{M}(a)
@inline lazymul(a, ::StaticInt{1}) = a
@inline lazymul(a, ::StaticInt{0}) = StaticInt{0}()
@inline lazymul(a::LazyMulAdd, ::StaticInt{1}) = a
@inline lazymul(::StaticInt{1}, a::LazyMulAdd) = a
@inline lazymul(a::LazyMulAdd, ::StaticInt{0}) = StaticInt{0}()
@inline lazymul(::StaticInt{0}, a::LazyMulAdd) = StaticInt{0}()
# @inline lazymul(a::LazyMulAdd, ::StaticInt{0}) = StaticInt{0}()
# @inline lazymul(::StaticInt{M}, b::MM{W,X}) where {W,M,X} = LazyMulAdd{M}(MM{W}(getfield(b, :data), StaticInt{M}()*StaticInt{X}()))
# @inline lazymul(a::MM{W,X}, ::StaticInt{M}) where {W,M,X} = LazyMulAdd{M}(MM{W}(getfield(a, :data), StaticInt{M}()*StaticInt{X}()))
@inline lazymul(a::MM{W,X}, ::StaticInt{M}) where {W,M,X} = MM{W}(vmul_nsw(StaticInt{M}(), data(a)), StaticInt{X}() * StaticInt{M}())
@inline lazymul(::StaticInt{M}, b::MM{W,X}) where {W,M,X} = MM{W}(vmul_nsw(StaticInt{M}(), data(b)), StaticInt{X}() * StaticInt{M}())
@inline lazymul(a::MM{W,X}, ::StaticInt{1}) where {W,X} = a
@inline lazymul(::StaticInt{1}, a::MM{W,X}) where {W,X} = a
@inline lazymul(a::MM{W,X}, ::StaticInt{0}) where {W,X} = StaticInt{0}()
@inline lazymul(::StaticInt{0}, a::MM{W,X}) where {W,X} = StaticInt{0}()

@inline lazymul(::StaticInt{M}, ::StaticInt{N}) where {M, N} = StaticInt{M}()*StaticInt{N}()

@inline lazymul(::StaticInt{M}, ::StaticInt{1}) where {M} = StaticInt{M}()
@inline lazymul(::StaticInt,    ::StaticInt{0}) = StaticInt{0}()

@inline lazymul(::StaticInt{1}, ::StaticInt{M}) where {M} = StaticInt{M}()
@inline lazymul(::StaticInt{0}, ::StaticInt   ) = StaticInt{0}()

@inline lazymul(::StaticInt{0}, ::StaticInt{0}) = StaticInt{0}()
@inline lazymul(::StaticInt{0}, ::StaticInt{1}) = StaticInt{0}()
@inline lazymul(::StaticInt{1}, ::StaticInt{0}) = StaticInt{0}()
@inline lazymul(::StaticInt{1}, ::StaticInt{1}) = StaticInt{1}()

@inline function lazymul(a::LazyMulAdd{M,O}, ::StaticInt{N}) where {M,O,N}
    LazyMulAdd(getfield(a, :data), StaticInt{M}() * StaticInt{N}(), StaticInt{O}() * StaticInt{N}())
end
@inline function lazymul(::StaticInt{M}, b::LazyMulAdd{N,O}) where {M,O,N}
    LazyMulAdd(getfield(b, :data), StaticInt{M}() * StaticInt{N}(), StaticInt{O}() * StaticInt{M}())
end

@inline function lazymul(a::LazyMulAdd{M,<:MM{W,X}}, ::StaticInt{N}) where {M,N,W,X}
    LazyMulAdd(MM{W}(getfield(a, :data), StaticInt{M}()*StaticInt{N}()*StaticInt{X}()), StaticInt{M}()*StaticInt{N}())
end
@inline function lazymul(::StaticInt{M}, b::LazyMulAdd{N,<:MM{W,X}}) where {M,N,W,X}
    LazyMulAdd(MM{W}(getfield(b, :data), StaticInt{M}()*StaticInt{N}()*StaticInt{X}()), StaticInt{M}()*StaticInt{N}())
end

@inline lazymul(a::LazyMulAdd{M,<:MM{W,X}}, ::StaticInt{0}) where {M,W,X} = StaticInt{0}()
@inline lazymul(::StaticInt{0}, b::LazyMulAdd{N,<:MM{W,X}}) where {N,W,X} = StaticInt{0}()
@inline lazymul(a::LazyMulAdd{M,<:MM{W,X}}, ::StaticInt{1}) where {M,W,X} = a
@inline lazymul(::StaticInt{1}, b::LazyMulAdd{N,<:MM{W,X}}) where {N,W,X} = b
@inline function lazymul(a::LazyMulAdd{M}, b::LazyMulAdd{N}) where {M,N}
    LazyMulAdd(vmul_nsw(getfield(a, :data), getfield(b, :data)), StaticInt{M}()*StaticInt{N}())
end

vmul_nsw(::LazyMulAdd{M,O}, ::StaticInt{0}) where {M,O} = Zero()
vmul_nsw(::StaticInt{0}, ::LazyMulAdd{M,O}) where {M,O} = Zero()
vmul_nsw(a::LazyMulAdd{M,O}, ::StaticInt{1}) where {M,O} = a
vmul_nsw(::StaticInt{1}, a::LazyMulAdd{M,O}) where {M,O} = a
@inline function vmul_nsw(a::LazyMulAdd{M,O}, ::StaticInt{N}) where {M,N,O}
    LazyMulAdd(getfield(a, :data), StaticInt{M}() * StaticInt{N}(), StaticInt{O}() * StaticInt{N}())
end
@inline function vmul_nsw(::StaticInt{N}, a::LazyMulAdd{M,O}) where {M,N,O}
    LazyMulAdd(getfield(a, :data), StaticInt{M}() * StaticInt{N}(), StaticInt{O}() * StaticInt{N}())
end
@inline vmul_nsw(a::LazyMulAdd{M,0}, i::IntegerTypesHW) where {M} = LazyMulAdd{M,0}(vmul_nsw(getfield(a,:data), i))
@inline vmul_nsw(a::LazyMulAdd{M,O}, i::IntegerTypesHW) where {M,O} = vmul_nsw(_materialize(a), i)

@inline function Base.:(>>>)(a::LazyMulAdd{M,O}, ::StaticInt{N}) where {M,O,N}
    LazyMulAdd(getfield(a, :data), StaticInt{M}() >>> StaticInt{N}(), StaticInt{O}() >>> StaticInt{N}())
end

# The approach with `add_indices` is that we try and make `vadd_nsw` behave well
# but for `i` and `j` type combinations where it's difficult,
# we can add specific `add_indices` methods that increment the pointer.
@inline function add_indices(p::Ptr, i, j) # generic fallback
    p, vadd_nsw(i, j)
end
@inline vadd_nsw(i::LazyMulAdd, ::Zero) = i
@inline vadd_nsw(::Zero, i::LazyMulAdd) = i

# These following two definitions normally shouldn't be hit
@inline function vadd_nsw(a::LazyMulAdd{M,O,MM{W,X,I}}, b::IntegerTypesHW) where {M,O,W,X,I}
    MM{W}(vadd_nsw(vmul_nsw(StaticInt{M}(), data(a)), b), StaticInt{X}() * StaticInt{M}())
end
@inline function vadd_nsw(b::IntegerTypesHW, a::LazyMulAdd{M,O,MM{W,X,I}}) where {M,O,W,X,I}
    MM{W}(vadd_nsw(vmul_nsw(StaticInt{M}(), data(a)), b), StaticInt{X}() * StaticInt{M}())
end

@inline vadd_nsw(a::LazyMulAdd{M,O}, b) where {M,O} = vadd_nsw(_materialize(a), b)
@inline vadd_nsw(b, a::LazyMulAdd{M,O}) where {M,O} = vadd_nsw(b, _materialize(a))
@inline vsub_nsw(a::LazyMulAdd{M,O}, b) where {M,O} = vsub_nsw(_materialize(a), b)
@inline vsub_nsw(b, a::LazyMulAdd{M,O}) where {M,O} = vsub_nsw(b, _materialize(a))

@inline vadd_nsw(a::LazyMulAdd{M,O,MM{W,X,I}}, ::Zero) where {M,O,W,X,I} = a
@inline vadd_nsw(::Zero, a::LazyMulAdd{M,O,MM{W,X,I}}) where {M,O,W,X,I} = a
@inline function vsub_nsw(a::LazyMulAdd{M,O,MM{W,X,I}}, b::IntegerTypesHW) where {M,O,W,X,I}
    MM{W}(vsub_nsw(vmul_nsw(StaticInt{M}(), data(a)), b), StaticInt{X}() * StaticInt{M}())
end
@inline function vsub_nsw(b::IntegerTypesHW, a::LazyMulAdd{M,O,MM{W,X,I}}) where {M,O,W,X,I}
    MM{W}(vsub_nsw(b, vmul_nsw(StaticInt{M}(), data(a))), (StaticInt{-1}() * StaticInt{X}()) * StaticInt{M}())
end

# because we should hit this instead:
@inline add_indices(p::Ptr, b::Integer, a::LazyMulAdd{M,O}) where {M,O} = (p + b, a)
@inline add_indices(p::Ptr, a::LazyMulAdd{M,O}, b::Integer) where {M,O} = (p + b, a)
# but in the case of `VecUnroll`s, which skip the `add_indices`, it's useful to still have the former two definitions.
# However, this also forces us to write:
@inline add_indices(p::Ptr, ::StaticInt{N}, a::LazyMulAdd{M,O}) where {M,O,N} = (p, vadd_nsw(a, StaticInt{N}()))
@inline add_indices(p::Ptr, a::LazyMulAdd{M,O}, ::StaticInt{N}) where {M,O,N} = (p, vadd_nsw(a, StaticInt{N}()))

@inline function vadd_nsw(::StaticInt{N}, a::LazyMulAdd{M,O}) where {N,M,O}
    LazyMulAdd(getfield(a, :data), StaticInt{M}(), StaticInt{O}()+StaticInt{N}())
end
@inline function vadd_nsw(a::LazyMulAdd{M,O}, ::StaticInt{N}) where {N,M,O}
    LazyMulAdd(getfield(a, :data), StaticInt{M}(), StaticInt{O}()+StaticInt{N}())
end
@inline function vadd_nsw(::StaticInt{N}, a::LazyMulAdd{M,O,MM{W,X,I}}) where {N,M,O,W,X,I}
    LazyMulAdd(getfield(a, :data), StaticInt{M}(), StaticInt{O}()+StaticInt{N}())
end
@inline function vadd_nsw(a::LazyMulAdd{M,O,MM{W,X,I}}, ::StaticInt{N}) where {N,M,O,W,X,I}
    LazyMulAdd(getfield(a, :data), StaticInt{M}(), StaticInt{O}()+StaticInt{N}())
end
@inline function vadd_nsw(a::LazyMulAdd{M,O}, b::LazyMulAdd{M,A}) where {M,O,A}
    LazyMulAdd(vadd_nsw(getfield(a, :data), getfield(b, :data)), StaticInt{M}(), StaticInt{O}()+StaticInt{A}())
end

@inline function vsub_nsw(::StaticInt{N}, a::LazyMulAdd{M,O}) where {N,M,O}
    LazyMulAdd(getfield(a, :data), StaticInt{-1}()*StaticInt{M}(), StaticInt{N}()-StaticInt{O}())
end
@inline vsub_nsw(::Zero, a::LazyMulAdd{M,O}) where {M,O} = StaticInt{-1}()*a
@inline function vsub_nsw(a::LazyMulAdd{M,O}, ::StaticInt{N}) where {N,M,O}
    LazyMulAdd(getfield(a, :data), StaticInt{M}(), StaticInt{O}()-StaticInt{N}())
end
@inline vsub_nsw(a::LazyMulAdd{M,O}, ::Zero) where {M,O} = a
@inline function vsub_nsw(::StaticInt{N}, a::LazyMulAdd{M,O,MM{W,X,I}}) where {N,M,O,W,X,I}
    LazyMulAdd(getfield(a, :data), StaticInt{-1}()*StaticInt{M}(), StaticInt{N}()-StaticInt{O}())
end
@inline vsub_nsw(::Zero, a::LazyMulAdd{M,O,MM{W,X,I}}) where {M,O,W,X,I} = StaticInt{-1}()*a
@inline function vsub_nsw(a::LazyMulAdd{M,O,MM{W,X,I}}, ::StaticInt{N}) where {N,M,O,W,X,I}
    LazyMulAdd(getfield(a, :data), StaticInt{M}(), StaticInt{O}()-StaticInt{N}())
end
@inline vsub_nsw(a::LazyMulAdd{M,O,MM{W,X,I}}, ::Zero) where {M,O,W,X,I} = a
@inline function vsub_nsw(a::LazyMulAdd{M,O}, b::LazyMulAdd{M,A}) where {M,O,A}
    LazyMulAdd(vsub_nsw(getfield(a, :data), getfield(b, :data)), StaticInt{M}(), StaticInt{O}()-StaticInt{A}())
end

@inline add_indices(p::Ptr, a::LazyMulAdd{M,O,V}, b::LazyMulAdd{N,P,J}) where {M,O,V<:AbstractSIMDVector,N,P,J<:IntegerTypes} = (gep(p, b), a)
@inline add_indices(p::Ptr, b::LazyMulAdd{N,P,J}, a::LazyMulAdd{M,O,V}) where {M,O,V<:AbstractSIMDVector,N,P,J<:IntegerTypes} = (gep(p, b), a)
@inline add_indices(p::Ptr, a::LazyMulAdd{M,O,V}, b::LazyMulAdd{M,P,J}) where {M,O,V<:AbstractSIMDVector,P,J<:IntegerTypes} = (p, vadd_nsw(a, b))
@inline add_indices(p::Ptr, b::LazyMulAdd{M,P,J}, a::LazyMulAdd{M,O,V}) where {M,O,V<:AbstractSIMDVector,P,J<:IntegerTypes} = (p, vadd_nsw(a, b))

@inline add_indices(p::Ptr, a::AbstractSIMDVector, b::LazyMulAdd{M,O,I}) where {M,O,I<:IntegerTypes} = (gep(p, b), a)
@inline add_indices(p::Ptr, b::LazyMulAdd{M,O,I}, a::AbstractSIMDVector) where {M,O,I<:IntegerTypes} = (gep(p, b), a)
@inline function add_indices(p::Ptr, ::MM{W,X,StaticInt{A}}, a::LazyMulAdd{M,O,T}) where {M,O,T<:IntegerTypes,A,W,X}
    gep(p, a), MM{W,X}(StaticInt{A}())
end
@inline function add_indices(p::Ptr, a::LazyMulAdd{M,O,T}, ::MM{W,X,StaticInt{A}}) where {M,O,T<:IntegerTypes,A,W,X}
    gep(p, a), MM{W,X}(StaticInt{A}())
end

@generated function add_indices(p::Ptr, a::LazyMulAdd{M,O,MM{W,X,StaticInt{I}}}, b::LazyMulAdd{N,P,J}) where {M,O,W,X,I,N,P,J<:IntegerTypes}
    d, r = divrem(M, N)
    if iszero(r)
        quote
            $(Expr(:meta,:inline))
            p, VectorizationBase.LazyMulAdd{$N,$(I*M)}(MM{$W,$d}(getfield(b, :data)))
        end
    else
        quote
            $(Expr(:meta,:inline))
            gep(p, b), a
        end
    end
end
@inline add_indices(p::Ptr, b::LazyMulAdd{N,P,J}, a::LazyMulAdd{M,O,MM{W,X,StaticInt{I}}}) where {M,O,W,X,I,N,P,J<:IntegerTypes} = add_indices(p, a, b)
@generated function vadd_nsw(a::LazyMulAdd{M,O,MM{W,X,StaticInt{I}}}, b::LazyMulAdd{N,P,J}) where {M,O,W,X,I,N,P,J<:IntegerTypes}
    d, r = divrem(M, N)
    if iszero(r)
        quote
            $(Expr(:meta,:inline))
            VectorizationBase.LazyMulAdd{$N,$(I*M)}(MM{$W,$d}(getfield(b, :data)))
        end
    else
        quote
            $(Expr(:meta,:inline))
            vadd_nsw(a, _materialize(b))
        end
    end
end
@inline vadd_nsw(b::LazyMulAdd{N,P,J}, a::LazyMulAdd{M,O,MM{W,X,StaticInt{I}}}) where {M,O,W,X,I,N,P,J<:IntegerTypes} = vadd_nsw(a, b)
@inline vadd_nsw(a::VecUnroll, b::LazyMulAdd) = VecUnroll(fmap(vadd_nsw, getfield(a, :data), b))
@inline vadd_nsw(b::LazyMulAdd, a::VecUnroll) = VecUnroll(fmap(vadd_nsw, b, getfield(a, :data)))
@inline vadd_nsw(a::LazyMulAdd, b::LazyMulAdd) = vadd_nsw(_materialize(a), _materialize(b))
@inline vsub_nsw(a::LazyMulAdd, b::LazyMulAdd) = vsub_nsw(_materialize(a), _materialize(b))

@generated function vsub_nsw(a::LazyMulAdd{M,O,MM{W,X,StaticInt{I}}}, b::LazyMulAdd{N,P,J}) where {M,O,W,X,I,N,P,J<:IntegerTypes}
    d, r = divrem(M, N)
    if iszero(r)
        quote
            $(Expr(:meta,:inline))
            VectorizationBase.LazyMulAdd{$N,$(I*M)}(-MM{$W,$d}(getfield(b, :data)))
        end
    else
        quote
            $(Expr(:meta,:inline))
            vsub_nsw(a, _materialize(b))
        end
    end
end
@inline vsub_nsw(b::LazyMulAdd{N,P,J}, a::LazyMulAdd{M,O,MM{W,X,StaticInt{I}}}) where {M,O,W,X,I,N,P,J<:IntegerTypes} = vsub_nsw(a, b)

@inline vsub_nsw(a::VecUnroll, b::LazyMulAdd) = VecUnroll(fmap(vsub_nsw, getfield(a, :data), b))
@inline vsub_nsw(b::LazyMulAdd, a::VecUnroll) = VecUnroll(fmap(vsub_nsw, b, getfield(a, :data)))

ext: jl


text:
module ARPACK

import ..LinAlg: BlasInt, blas_int, ARPACKException

## aupd and eupd wrappers

function aupd_wrapper(T, matvecA::Function, matvecB::Function, solveSI::Function, n::Integer,
                      sym::Bool, cmplx::Bool, bmat::ASCIIString,
                      nev::Integer, ncv::Integer, which::ASCIIString,
                      tol::Real, maxiter::Integer, mode::Integer, v0::Vector)

    lworkl = cmplx ? ncv * (3*ncv + 5) : (sym ? ncv * (ncv + 8) :  ncv * (3*ncv + 6) )
    TR = cmplx ? T.types[1] : T
    TOL = Array(TR, 1)
    TOL[1] = tol

    v      = Array(T, n, ncv)
    workd  = Array(T, 3*n)
    workl  = Array(T, lworkl)
    rwork  = cmplx ? Array(TR, ncv) : Array(TR, 0)

    if isempty(v0)
    	resid  = Array(T, n)
    	info   = zeros(BlasInt, 1)
    else
	resid  = deepcopy(v0)
    	info   = ones(BlasInt, 1)
    end
    iparam = zeros(BlasInt, 11)
    ipntr  = zeros(BlasInt, (sym && !cmplx) ? 11 : 14)
    ido    = zeros(BlasInt, 1)

    iparam[1] = blas_int(1)       # ishifts
    iparam[3] = blas_int(maxiter) # maxiter
    iparam[7] = blas_int(mode)    # mode

    zernm1 = 0:(n-1)

    while true
        if cmplx
            naupd(ido, bmat, n, which, nev, TOL, resid, ncv, v, n,
                  iparam, ipntr, workd, workl, lworkl, rwork, info)
        elseif sym
            saupd(ido, bmat, n, which, nev, TOL, resid, ncv, v, n,
                  iparam, ipntr, workd, workl, lworkl, info)
        else
            naupd(ido, bmat, n, which, nev, TOL, resid, ncv, v, n,
                  iparam, ipntr, workd, workl, lworkl, info)
        end

        # Check for warnings and errors
        # Refer to ex-*.doc files in ARPACK/DOCUMENTS for calling sequence
        if info[1] == 3; warn("try eigs/svds with a larger value for ncv"); end
        if info[1] == 1; warn("maximum number of iterations reached; check nconv for number of converged eigenvalues"); end
        if info[1] < 0; throw(ARPACKException(info[1])); end

        load_idx = ipntr[1]+zernm1
        store_idx = ipntr[2]+zernm1
        x = workd[load_idx]
        if ido[1] == -1
            if mode == 1
                workd[store_idx] = matvecA(x)
            elseif mode == 2
                if sym
                    temp = matvecA(x)
                    workd[load_idx] = temp    # overwrite as per Remark 5 in dsaupd.f
                    workd[store_idx] = solveSI(temp)
                else
                    workd[store_idx] = solveSI(matvecA(x))
                end
            elseif mode == 3
                if bmat == "I"
                    workd[store_idx] = solveSI(x)
                elseif bmat == "G"
                    workd[store_idx] = solveSI(matvecB(x))
                end
            end
        elseif ido[1] == 1
            if mode == 1
                workd[store_idx] = matvecA(x)
            elseif mode == 2
                if sym
                    temp = matvecA(x)
                    workd[load_idx] = temp
                    workd[store_idx] = solveSI(temp)
                else
                    workd[store_idx] = solveSI(matvecA(x))
                end
            elseif mode == 3
                workd[store_idx] = solveSI(workd[ipntr[3]+zernm1])
            end
        elseif ido[1] == 2
            workd[store_idx] = matvecB(x)
        elseif ido[1] == 99
            break
        else
            error("Internal ARPACK error")
        end
    end

    return (resid, v, n, iparam, ipntr, workd, workl, lworkl, rwork, TOL)
end

function eupd_wrapper(T, n::Integer, sym::Bool, cmplx::Bool, bmat::ASCIIString,
                      nev::Integer, which::ASCIIString, ritzvec::Bool,
                      TOL::Array, resid, ncv::Integer, v, ldv, sigma, iparam, ipntr,
                      workd, workl, lworkl, rwork)

    howmny = "A"
    select = Array(BlasInt, ncv)
    info   = zeros(BlasInt, 1)

    if cmplx

        d = Array(T, nev+1)
        sigma = ones(T, 1)*sigma
        workev = Array(T, 2ncv)
        neupd(ritzvec, howmny, select, d, v, ldv, sigma, workev,
              bmat, n, which, nev, TOL, resid, ncv, v, ldv,
              iparam, ipntr, workd, workl, lworkl, rwork, info)
        if info[1] != 0; throw(ARPACKException(info[1])); end
        return ritzvec ? (d[1:nev], v[1:n, 1:nev],iparam[5],iparam[3],iparam[9],resid) : (d[1:nev],iparam[5],iparam[3],iparam[9],resid)

    elseif sym

        d = Array(T, nev)
        sigma = ones(T, 1)*sigma
        seupd(ritzvec, howmny, select, d, v, ldv, sigma,
              bmat, n, which, nev, TOL, resid, ncv, v, ldv,
              iparam, ipntr, workd, workl, lworkl, info)
        if info[1] != 0; throw(ARPACKException(info[1])); end
        return ritzvec ? (d, v[1:n, 1:nev],iparam[5],iparam[3],iparam[9],resid) : (d,iparam[5],iparam[3],iparam[9],resid)

    else

        dr     = Array(T, nev+1)
        di     = Array(T, nev+1)
        sigmar = ones(T, 1)*real(sigma)
        sigmai = ones(T, 1)*imag(sigma)
        workev = Array(T, 3*ncv)
        neupd(ritzvec, howmny, select, dr, di, v, ldv, sigmar, sigmai,
              workev, bmat, n, which, nev, TOL, resid, ncv, v, ldv,
              iparam, ipntr, workd, workl, lworkl, info)
        if info[1] != 0; throw(ARPACKException(info[1])); end
        evec = complex(zeros(T, n, nev+1), zeros(T, n, nev+1))
        j = 1
        while j <= nev
            if di[j] == 0.0
                evec[:,j] = v[:,j]
            else
                evec[:,j]   = v[:,j] + im*v[:,j+1]
                evec[:,j+1] = v[:,j] - im*v[:,j+1]
                j += 1
            end
            j += 1
        end
        d = complex(dr[1:nev],di[1:nev])
        return ritzvec ? (d, evec[1:n, 1:nev],iparam[5],iparam[3],iparam[9],resid) : (d,iparam[5],iparam[3],iparam[9],resid)
    end

end

for (T, saupd_name, seupd_name, naupd_name, neupd_name) in
    ((:Float64, :dsaupd_, :dseupd_, :dnaupd_, :dneupd_),
     (:Float32, :ssaupd_, :sseupd_, :snaupd_, :sneupd_))
    @eval begin

        function naupd(ido, bmat, n, evtype, nev, TOL::Array{$T}, resid::Array{$T}, ncv, v::Array{$T}, ldv,
                       iparam, ipntr, workd::Array{$T}, workl::Array{$T}, lworkl, info)

            ccall(($(string(naupd_name)), :libarpack), Void,
                  (Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt},
                   Ptr{$T}, Ptr{$T}, Ptr{BlasInt}, Ptr{$T}, Ptr{BlasInt},
                   Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{BlasInt}, Ptr{BlasInt}, Clong, Clong),
                  ido, bmat, &n, evtype, &nev, TOL, resid, &ncv, v, &ldv,
                  iparam, ipntr, workd, workl, &lworkl, info, sizeof(bmat), sizeof(evtype))
        end

        function neupd(rvec, howmny, select, dr, di, z, ldz, sigmar, sigmai,
                  workev::Array{$T}, bmat, n, evtype, nev, TOL::Array{$T}, resid::Array{$T}, ncv, v, ldv,
                  iparam, ipntr, workd::Array{$T}, workl::Array{$T}, lworkl, info)

            ccall(($(string(neupd_name)), :libarpack), Void,
                  (Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{$T},
                   Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{$T}, Ptr{Uint8}, Ptr{BlasInt},
                   Ptr{Uint8}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{BlasInt}, Ptr{$T},
                   Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T},
                   Ptr{BlasInt}, Ptr{BlasInt}, Clong, Clong, Clong),
                  &rvec, howmny, select, dr, di, z, &ldz, sigmar, sigmai,
                  workev, bmat, &n, evtype, &nev, TOL, resid, &ncv, v, &ldv,
                  iparam, ipntr, workd, workl, &lworkl, info,
                  sizeof(howmny), sizeof(bmat), sizeof(evtype))
        end

        function saupd(ido, bmat, n, which, nev, TOL::Array{$T}, resid::Array{$T}, ncv, v::Array{$T}, ldv,
                       iparam, ipntr, workd::Array{$T}, workl::Array{$T}, lworkl, info)

            ccall(($(string(saupd_name)), :libarpack), Void,
                  (Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt},
                   Ptr{$T}, Ptr{$T}, Ptr{BlasInt}, Ptr{$T}, Ptr{BlasInt},
                   Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{BlasInt}, Ptr{BlasInt}, Clong, Clong),
                  ido, bmat, &n, which, &nev, TOL, resid, &ncv, v, &ldv,
                  iparam, ipntr, workd, workl, &lworkl, info, sizeof(bmat), sizeof(which))

        end

        function seupd(rvec, howmny, select, d, z, ldz, sigma,
                       bmat, n, evtype, nev, TOL::Array{$T}, resid::Array{$T}, ncv, v::Array{$T}, ldv,
                       iparam, ipntr, workd::Array{$T}, workl::Array{$T}, lworkl, info)

            ccall(($(string(seupd_name)), :libarpack), Void,
                  (Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{BlasInt}, Ptr{$T},
                   Ptr{Uint8}, Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt},
                   Ptr{$T}, Ptr{$T}, Ptr{BlasInt}, Ptr{$T}, Ptr{BlasInt}, Ptr{BlasInt},
                   Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{BlasInt}, Ptr{BlasInt}, Clong, Clong, Clong),
                  &rvec, howmny, select, d, z, &ldz, sigma,
                  bmat, &n, evtype, &nev, TOL, resid, &ncv, v, &ldv,
                  iparam, ipntr, workd, workl, &lworkl, info, sizeof(howmny), sizeof(bmat), sizeof(evtype))
        end

    end
end

for (T, TR, naupd_name, neupd_name) in
    ((:Complex128, :Float64, :znaupd_, :zneupd_),
     (:Complex64,  :Float32, :cnaupd_, :cneupd_))
    @eval begin

        function naupd(ido, bmat, n, evtype, nev, TOL::Array{$TR}, resid::Array{$T}, ncv, v::Array{$T}, ldv,
                       iparam, ipntr, workd::Array{$T}, workl::Array{$T}, lworkl,
                       rwork::Array{$TR}, info)

            ccall(($(string(naupd_name)), :libarpack), Void,
                  (Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt},
                   Ptr{$TR}, Ptr{$T}, Ptr{BlasInt}, Ptr{$T}, Ptr{BlasInt},
                   Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{BlasInt},
                   Ptr{$TR}, Ptr{BlasInt}),
                  ido, bmat, &n, evtype, &nev, TOL, resid, &ncv, v, &ldv,
                  iparam, ipntr, workd, workl, &lworkl, rwork, info)

        end

        function neupd(rvec, howmny, select, d, z, ldz, sigma, workev::Array{$T},
                       bmat, n, evtype, nev, TOL::Array{$TR}, resid::Array{$T}, ncv, v::Array{$T}, ldv,
                       iparam, ipntr, workd::Array{$T}, workl::Array{$T}, lworkl,
                       rwork::Array{$TR}, info)

            ccall(($(string(neupd_name)), :libarpack), Void,
                  (Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{BlasInt},
                   Ptr{$T}, Ptr{$T}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt},
                   Ptr{$TR}, Ptr{$T}, Ptr{BlasInt}, Ptr{$T}, Ptr{BlasInt}, Ptr{BlasInt},
                   Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{BlasInt}, Ptr{$TR}, Ptr{BlasInt}),
                  &rvec, howmny, select, d, z, &ldz, sigma, workev,
                  bmat, &n, evtype, &nev, TOL, resid, &ncv, v, &ldv,
                  iparam, ipntr, workd, workl, &lworkl, rwork, info)

        end

    end
end

end # module ARPACK

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
using Gadfly, Compat

plot(layer((x,y) -> x*exp(-(x-(@compat round(Int, x)))^2-y^2), -8., 8, -2., 2),
     layer((x, y) -> sqrt(hypot(x, y)), -2, 2, -1, 1),
     layer(x=[3], y=[1.5], Geom.point))


text:
function setup_dual_infeasibility()
        options = Dict(:verbose => false,
                       :eps_abs => 1e-05,
                       :eps_rel => 1e-05,
               :eps_prim_inf => 1e-15,
               :check_termination => 1)
    return options
end

tol = 1e-5

@testset "dual_infeasibility" begin

    @testset "dual_infeasible_lp" begin
        P = spzeros(2, 2)
        q = [2.; -1.]
        A = sparse(I, 2, 2)
        u = Inf * ones(2)
        l = [0.; 0.]
        options = setup_dual_infeasibility()

        model = OSQP.Model()
        OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)
        results = OSQP.solve!(model)

        @test results.info.status == :Dual_infeasible

    end

    @testset "dual_infeasible_qp" begin
        P = sparse(Diagonal([4.; 0.]))
        q = [0.; 2]
        A = sparse([1. 1.; -1. 1])
        u = [2.; 3]
        l = -Inf * ones(2)
        options = setup_dual_infeasibility()

        model = OSQP.Model()
        OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)
        results = OSQP.solve!(model)

        @test results.info.status == :Dual_infeasible

    end

    @testset "primal_dual_infeasible" begin
        P = spzeros(2, 2)
        q = [-1.; -1.]
        A = sparse([1. -1.; -1. 1; 1. 0; 0. 1])
        u = Inf * ones(4)
        l = [1., 1., 0., 0.]
        options = setup_dual_infeasibility()

        model = OSQP.Model()
        OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)

        # Warm start to avoid infeasibility detection at first step
        OSQP.warm_start!(model, x=[50.; 30.], y=[-2.; -2; -2; -2])

        results = OSQP.solve!(model)

        @test results.info.status == :Dual_infeasible

    end
end


text:
"""
The TinStandard type contains the state variables (swe), the inputs
(prec, tair) for one time step, the parameters (param) and the time step
length (tstep) for a standard temperature index model.
"""
type TinStandard <: Snow

  swe::Array{Float64,1}
  lw::Array{Float64,1}
  prec::Array{Float64,1}
  tair::Array{Float64,1}
  q_sim::Array{Float64,1}
  param::Array{Float64,1}
  frac::Array{Float64,1}
  tstep::Float64
  time::DateTime

end

"""
    TinStandard(tstep, time, frac)

Constructor for TinStandard with predefined state variables, parameters and inputs.
The time step (tstep) is given as a fraction of one day. Thus, for hourly input
data tstep should be set to 1/24. The fraction of elevation bands should sum
up to unity.
"""
function TinStandard(tstep::Float64, time::DateTime, frac::Array{Float64,1})

  nzones = length(frac)
  swe    = zeros(Float64, nzones)
  lw     = zeros(Float64, nzones)
  prec   = zeros(Float64, nzones)
  tair   = zeros(Float64, nzones)
  q_sim  = zeros(Float64, nzones)
  param  = zeros(Float64, 5)

  TinStandard(swe, lw, prec, tair, q_sim, param, frac, tstep, time)

end

"""
    TinStandard(tstep, time, param, frac)

Constructor for TinStandard with predefined state variables and inputs.
The time step (tstep) is given as a fraction of one day. Thus, for hourly input
data tstep should be set to 1/24. The fraction of elevation bands should sum
up to unity.
"""
function TinStandard(tstep::Float64, time::DateTime, param::Array{Float64,1}, frac::Array{Float64,1})

  nzones = length(frac)
  swe    = zeros(Float64, nzones)
  lw     = zeros(Float64, nzones)
  prec   = zeros(Float64, nzones)
  tair   = zeros(Float64, nzones)
  q_sim  = zeros(Float64, nzones)

  TinStandard(swe, lw, prec, tair, q_sim, param, frac, tstep, time)

end

"""
    init_states(mdata::TinStandard)

Initilize the state variables of the model.
"""
function init_states(mdata::TinStandard)

  for ireg in eachindex(mdata.swe)
    mdata.swe[ireg] = 0.
    mdata.lw[ireg] = 0.
  end

end

"""
    get_param_range(mdata::TinStandard)

Get allowed parameter ranges for the calibration of the model.
"""
function get_param_range(mdata::TinStandard)

  param_range_snow = [(-0.5, 0.5), (0.5, 3.0), (0.5, 4.0), (0.01, 0.10), (0.5, 2.0)]

end

"""
    get_states(mdata::TinStandard)

Get state variables for computing penelty during calibration.

"""
function get_states(mdata::TinStandard)

  return [mean(mdata.swe); mean(mdata.lw)]

end

"""
    assign_param(mdata::TinStandard, param::Array{Float64,1})

Assign parameter values to the TinBasic type.
"""
function assign_param(mdata::TinStandard, param::Array{Float64,1})

  for ireg in eachindex(mdata.param)
    mdata.param[ireg] = param[ireg]
  end

end

"""
    enkf_snow(mdata::Array{TinStandard}, obs_ens, q_sim)

Implementation of the ensemble Kalman filter for the standard temperature index
snow model.
"""
function enkf_snow(mdata::Array{TinStandard}, obs_ens, q_sim)

  nens = length(mdata)

  # Allocate arrays

  swe = zeros(Float64, length(mdata[1].swe), nens)
  lw = zeros(Float64, length(mdata[1].lw), nens)

  # Add states to arrays

  for iens = 1:nens

    swe[:, iens] = mdata[iens].swe
    lw[:, iens] = mdata[iens].lw

  end

  # Run ensemble kalman filter

  swe = enkf(swe, obs_ens, q_sim)
  lw = enkf(lw, obs_ens, q_sim)

  # Check limits of states

  swe[swe .< 0] = 0.

  lwmax = swe * mdata[1].param[4]

  lw[lw .< 0] = 0.
  lw[lw .> lwmax] = 0.

  # Add arrays to states

  for iens = 1:nens

    mdata[iens].swe = swe[:, iens]
    mdata[iens].lw = lw[:, iens]

  end

  nothing

end

"""
    run_timestep(mdata::TinStandard)

Propagate the model one time step and compute simulated snowpack dischage.
"""
function run_timestep(mdata::TinStandard)

    # Parameters

    tth     = mdata.param[1]
    ddf_min = mdata.param[2]
    ddf_max = mdata.param[3]
    whcap   = mdata.param[4]
    pcorr   = mdata.param[5]

    ddf = compute_ddf(mdata.time, ddf_min, ddf_max)

    for ireg in eachindex(mdata.swe)

        swe = mdata.swe[ireg]
        lw = mdata.lw[ireg]

        # Compute solid and liquid precipitation

        psolid, pliquid = split_prec(mdata.prec[ireg], mdata.tair[ireg], tth)

        # Apply snowfall correction factor

        psolid = pcorr * psolid

        # Compute frozen water in snowpack

        fw  = swe - lw

        # Constrain fw and lw

        if fw < 0.0
            fw = 0.0
        end

        if lw < 0.0
            lw = 0.0
        end

        # Compute melt rates

        melt = compute_melt(mdata.tair[ireg], fw, ddf, tth)

        # Compute refreezing rates

        refr = compute_refreeze(mdata.tair[ireg], lw, ddf, tth)

        # Massbalance for frozen water in snowpack

        dfw = psolid + refr - melt
        fw  = fw + dfw

        # Massbalance for liquid water in snowpack

        dlw = pliquid + melt - refr
        lw  = lw + dlw

        lwmax = fw * whcap / (1 - whcap)

        q_sim = lw - lwmax

        if q_sim < 0.
            q_sim = 0.
        end

        lw = lw - q_sim

        # Massbalance for snowpack

        swe = fw + lw

        # Assign final results

        mdata.swe[ireg] = swe
        mdata.lw[ireg] = lw
        mdata.q_sim[ireg] = q_sim

    end

    mdata.time += Dates.Hour(mdata.tstep)

    return nothing

end

# Compute snowmelt

function compute_melt(tair::Float64, fw::Float64, ddf::Float64, tth_melt::Float64)

    m_melt = 0.5

    t_m = (tair - tth_melt) / m_melt
    melt = ddf * m_melt * (t_m + log(1. + exp(-t_m)))

    melt = min(max(0.,fw), melt)

    return melt

end

# Compute refreeze

function compute_refreeze(tair::Float64, lw::Float64, ddf::Float64, tth_refr::Float64)

    m_refr = 0.5

    t_m = (tth_refr - tair) / m_refr
    refr = ddf * m_refr * (t_m + log(1. + exp(-t_m)))

    if refr > lw
        refr = lw
    end

    refr = min(lw,refr)

end

max_stars_repo_name: jmgnve/Vann

max_forks_repo_name: jmgnve/Vann


text:
"""
    vocab_count(corpus, vocab; max_vocab=100_000, min_count=10, verbose=2)

Extracts unigram counts.

# Arguments
  * `corpus::AbstractString` input corpus file path
  * `vocab::AbstractString` output vocabulary file path

# Keyword arguments
  * `verbose::Int` set verbosity: 0, 1 or 2 (default)
  * `max_vocab::Int` upper bound on vocabulary size, i.e. keep the <int> most frequent words. The minimum frequency words are randomly sampled so as to obtain an even distribution over the alphabet.
  * `min_count::Int` lower limit such that words which occur fewer than <int> times are discarded.

# Examples
```
julia> vocab_count("corpus.txt", "vocab.txt", verbose=2, max_vocab=100_000, min_count=10)
```
"""
function vocab_count(corpus::AbstractString,
                     vocab::AbstractString;
                     max_vocab::Int=100_000,
                     min_count::Int=10,
                     verbose::Int=2)
    command = joinpath(dirname(@__FILE__), "..",
                       "deps", "src", "GloVe-c", "build", "./vocab_count")
    parameters = String[]
    args = ["-max-vocab", "-min-count", "-verbose"]
    values = [max_vocab, min_count, verbose]

    for (arg, value) in zip(args, values)
        push!(parameters, arg)
        push!(parameters, string(value))
    end

    run(pipeline(pipeline(corpus, `$command $parameters`), vocab))
    return nothing
end

"""
    cooccur(corpus, vocab, cooccurrences; verbose=2, symmetric=0, window_size=15, memory=4.0, max_product=nothing, overflow_length=nothing, overflow_file="overflow", distance_weighting=1)

Calculates word-word cooccurrence statistics.

# Arguments
  * `corpus::AbstractString` input corpus file path
  * `vocab::AbstractString` input vocabulary file path (the vocabulary contains truncated unigram counts, produced by `vocab_count`)
  * `cooccurrences::AbstractString` output cooccurrences file path

# Keyword arguments
  * `verbose::Int` set verbosity: 0, 1, 2 (default) or 3
  * `symmetric::Int` if <int> = 0, only use left context; if <int> = 1 (default), use left and right
  * `window_size::Int` number of context words to the left (and to the right, if symmetric = 1); default 15
  * `memory::Float64` soft limit for memory consumption, in GB -- based on simple heuristic, so not extremely accurate; default 4.0
  * `max_product::Union{Nothing, Int}` limit the size of dense cooccurrence array by specifying the max product <int> of the frequency counts of the two cooccurring words. This value overrides that which is automatically produced by `memory`. Typically only needs adjustment for use with very large corpora
  * `overflow_length::Union{Nothing, Int}` limit to length <int> the sparse overflow array, which buffers cooccurrence data that does not fit in the dense array, before writing to disk. This value overrides that which is automatically produced by `memory`. Typically only needs adjustment for use with very large corpora
  * `overflow_file::String` filename, excluding extension, for temporary files; default "overflow"
  * `distance_weighting::Int` if <int> = 0, do not weight cooccurrence count by distance between words; if <int> = 1 (default), weight the cooccurrence count by inverse of distance between words"

# Examples
```
# It is assumed that vocab.txt exists and has been created by `vocab_count`
julia> cooccur("corpus.txt", "vocab.txt", "cooccurrences.bin", verbose=2, symmetric=0, window_size=10, memory=8.0, overflow_file="tempoverflow")
```
"""
function cooccur(corpus::AbstractString,
                 vocab::AbstractString,
                 cooccurrences::AbstractString;
                 verbose::Int=2,
                 symmetric::Int=1,
                 window_size::Int=15,
                 memory::Float64=4.0,
                 max_product::Union{Nothing, Int}=nothing,
                 overflow_length::Union{Nothing, Int}=nothing,
                 overflow_file::String="overflow",
                 distance_weighting::Int=1)
    command = joinpath(dirname(@__FILE__), "..",
                       "deps", "src", "GloVe-c", "build", "./cooccur")
    parameters = String[]
    args = ["-verbose", "-symmetric", "-window-size", "-memory",
            "-max-product", "-overflow-length", "-overflow-file",
            "-distance-weighting"]
    values = [verbose, symmetric, window_size, memory,
              max_product, overflow_length, overflow_file,
              distance_weighting]

    for (arg, value) in zip(args, values)
        if value != nothing
            push!(parameters, arg)
            push!(parameters, string(value))
        end
    end
    # Push '-vocab-file' separately
    push!(parameters, "-vocab-file")
    push!(parameters, vocab)

    run(pipeline(pipeline(corpus, `$command $parameters`), cooccurrences))
    return nothing
end

"""
    shuffle(cooccurences, shuffled; verbose=2, memory=4.0, array_size=nothing, temp_file="temp_shuffle")

Shuffles entries of word-word cooccurrence files.

# Arguments
  * `cooccurences::AbstractString` input cooccurrences file path
  * `shuffled::AbstractString` output shuffled cooccurences file path

# Keyword arguments
  * `verbose::Int` set verbosity: 0, 1, or 2 (default)
  * `memory::Float64` soft limit for memory consumption, in GB; default 4.0
  * `array_size::Union{Nothing, Int}` limit to length <int> the buffer which stores chunks of data to shuffle before writing to disk. This value overrides that which is automatically produced by `memory`
  * `temp_file`::String filename, excluding extension, for temporary files; default "temp_shuffle"

# Examples
```
# It is assumed that cooccurences.bin exists and has been created by `cooccur`
julia> shuffle("cooccurences.bin", "shuffled.bin", verbose=2, memory=8.0)
```
"""
function shuffle(cooccurrences::AbstractString,
                 shuffled::AbstractString;
                 verbose::Int=2,
                 memory::Float64=4.0,
                 array_size::Union{Nothing, Int}=nothing,
                 temp_file::String="temp_shuffle")
    command = joinpath(dirname(@__FILE__), "..",
                       "deps", "src", "GloVe-c", "build", "./shuffle")
    parameters = String[]
    args = ["-verbose", "-memory", "-array-size", "-temp-file"]
    values = [verbose, memory, array_size, temp_file]

    for (arg, value) in zip(args, values)
        if value != nothing
            push!(parameters, arg)
            push!(parameters, string(value))
        end
    end

    run(pipeline(pipeline(cooccurrences, `$command $parameters`), shuffled))
    return nothing
end

"""
    glove(shuffled, vocab, vectors="vectors"; verbose=2, write_header=0, vector_size=50, threads=8, iter=25, eta=0.05, alpha=0.75, x_max=100.0, binary=0, model=2, gradsq_file="gradsq", save_gradsq=0, checkpoint_every=0)

GloVe: Global Vectors for Word Representation, v0.2

# Arguments
  * `shuffled::AbstractString` binary input file of shuffled cooccurrence data (produced by `cooccur` and `shuffle`)
  * `vocab::AbstractString` file containing vocabulary (truncated unigram counts, produced by `vocab_count`)
  * `vectors::AbstractString` filename, excluding extension, for word vector output; default "vectors"

# Keyword arguments
  * `verbose::Int` set verbosity: 0, 1, or 2 (default)
  * `write_header::Int` if 1, write vocab_size/vector_size as first line. Do nothing if 0 (default)
  * `vector_size::Int` dimension of word vector representations (excluding bias term); default 50
  * `threads::Int` number of threads; default 8
  * `iter::Int` number of training iterations; default 25
  * `eta::Float64` initial learning rate; default 0.05
  * `alpha::Float64` parameter in exponent of weighting function; default 0.75
  * `x_max::Float64` parameter specifying cutoff in weighting function; default 100.0
  * `binary::Int` save output in binary format (0: text, 1: binary, 2: both); default 0
  * `model::Int` model for word vector output (for text output only); default 2
        0: output all data, for both word and context word vectors, including bias terms
        1: output word vectors, excluding bias terms
        2: output word vectors + context word vectors, excluding bias terms
  * `gradsq_file::String` filename, excluding extension, for squared gradient output; default "gradsq"
  * `save_gradsq::Int` save accumulated squared gradients; default 0 (off); ignored if `gradsq_file`
  * `checkpoint_every::Int` checkpoint a  model every <int> iterations; default 0 (off)

# Examples
```
# It is assumed that:
#  - shuffled.bin exists and has been created by `shuffled`
#  - vocab.txt exists and has been created by `vocab_count`
julia> glove("shuffled.bin", "vocab.txt", "vectors", gradsq_file="gradsq", verbose=2, vector_size=100, threads=16, alpha=0.75, x_max=100.0, eta=0.05, binary=2, model=2)
```
"""

#printf("./glove -input-file cooccurrence.shuf.bin -vocab-file vocab.txt -save-file vectors -gradsq-file gradsq -verbose 2 -vector-size 100 -threads 16 -alpha 0.75 -x-max 100.0 -eta 0.05 -binary 2 -model 2\n\n");
function glove(shuffled::AbstractString,
               vocab::AbstractString,
               vectors::AbstractString="vectors";
               verbose::Int=2,
               write_header::Int=0,
               vector_size::Int=50,
               threads::Int=8,
               iter::Int=25,
               eta::Float64=0.05,
               alpha::Float64=0.75,
               x_max::Float64=100.0,
               binary::Int=0,
               model::Int=2,
               gradsq_file::Union{Nothing, String}=nothing,
               save_gradsq::Int=0,
               checkpoint_every::Int=0)
    command = joinpath(dirname(@__FILE__), "..",
                       "deps", "src", "GloVe-c", "build", "./glove")
    parameters = String[]
    args = ["-verbose", "-write-header", "-vector-size", "-threads", "-iter",
            "-eta", "-alpha", "-x-max", "-binary", "-model", "-save-gradsq",
            "-checkpoint-every"]
    values = [verbose, write_header, vector_size, threads, iter,
              eta, alpha, x_max, binary, model, save_gradsq,
              checkpoint_every]

    for (arg, value) in zip(args, values)
        push!(parameters, arg)
        push!(parameters, string(value))
    end

    # Handle gradsq_file and save_gradsq
    if gradsq_file != nothing
        push!(parameters, "-gradsq-file")
        push!(parameters, string(gradsq_file))
    end
    if gradsq_file == nothing && save_gradsq != 0
        push!(parameters, "-save-gradsq")
        push!(parameters, string(save_gradsq))
        push!(parameters, "-gradsq-file")
        push!(parameters, "gradsq")
    end
    # Handle input and output arguments
    push!(parameters, "-input-file")
    push!(parameters, shuffled)
    push!(parameters, "-vocab-file")
    push!(parameters, vocab)
    push!(parameters, "-save-file")
    push!(parameters, vectors)

    run(`$command $parameters`)
    return nothing
end


text:
using ForwardDiff
#=
Functions to compute chemical potential-related characteristics
=#

function exc_helmholtz(
    mix::MBWREoSMixture,
    nmol,
    vol,
    RT;
    buf = thermo_buffer(mix, nmol)
)
    Tx = RT / GAS_CONSTANT_SI
    total_mol = sum(nmol)
    ρ = total_mol / vol
    fx, hx = scaling_coeffs(mix, nmol, Tx; buf)
    Tref = Tx / fx
    ρref = ρ * hx

    ρmoll = ρref * 1e-3
    G1, G2, G3, G4, G5 = (
        -0.2804337729e-3,
        0.1180666107,
        -0.3756325860e+1,
        0.5624374521e3,
        -0.9354759605e5,
    )
    G6, G7, G8, G9 = (
        -0.4557405505e-4,
        0.1530044332,
        -0.1078107476e+3,
        0.2218072099e+5,
    )
    G10, G11, G12 = (
        0.6629473971e-5,
        -0.6199354447e-2,
        0.6754207966e+1,
    )
    G13 = 0.6472837570e-3
    G14, G15 = -0.6804325262e-1, -0.9726162355e+1
    G16 = 0.5097956459e-2
    G17, G18 = -0.1004655900e-3, 0.4363693352e-1
    G19 = -0.1249351947e-2
    G20, G21 = 0.2644755879e+5, -0.7944237270e+7
    G22, G23 = -0.7299920845e+3, 0.5381095003e+8
    G24, G25 = 0.3450217377e+1, 0.9936666689e+3
    G26, G27 = -0.2166699036, -0.1612103424e+5
    G28, G29 = -0.3633126990e-3, 0.1108612343e+1
    G30, G31, G32 = -0.1330932838e-4, -0.3157701101e-2, 0.1423083811
    iTsq = inv(Tref^2)
    iT = inv(Tref)
    a2 = G1 * Tref + G2 * sqrt(Tref) + G3 + (G4 * Tref + G5) * iTsq
    a3 = G6 * Tref + G7 + (G8 * Tref + G9) * iTsq
    a4 = G10 * Tref + G11 + G12 * iT
    a5 = G13
    a6 = (G14 * Tref + G15) * iTsq
    a7 = G16 * iT
    a8 = (G17 * Tref + G18) * iTsq
    a9 = G19 * iTsq

    b3 = (G20 * Tref + G21) * iT * iTsq
    b5 = (G22 * Tref^2 + G23) * iTsq^2
    b7 = (G24 * Tref + G25) * iT * iTsq
    b9 = (G26 * Tref^2 + G27) * iTsq^2
    b11 = (G28 * Tref + G29) * iT * iTsq
    b13 = (G32 + Tref * (G31 + Tref * G30)) * iTsq^2

    a = (a2, a3 / 2, a4 / 3, a5 / 4, a6 / 5, a7 / 6, a8 / 7, a9 / 8)
    Aexc = evalpoly(ρmoll, a) * ρmoll * 1e3
    ρr = C3_VC * ρref
    ρcmoll = 1e-3 / C3_VC
    In = 1 - exp(-ρr^2)
    bfactor = ρcmoll^2
    Ifactor = ρr^2
    Ar = 0.5 * In * b3 * bfactor
    for (n, b) in enumerate((b5, b7, b9, b11, b13))
        bfactor *= ρcmoll^2
        In = n * In - Ifactor * exp(-ρr^2)
        Ar += 0.5 * In * b * bfactor
        Ifactor *= ρr^2
    end
    Aexc += Ar * 1e3
    return Aexc * fx * total_mol
end

"""
    log_c_activity(mixture, nmol, volume, RT[; buf])

Return vector of ln(c_a) - logarithm of activity coefficients
for components of `mixture` at given `nmol`, `volume`, `RT`.
If buffers are provided as keyword arguments, their contents
are modified during the intermediate calculations.

# Arguments

- `mix`: mixture
- `nmol::AbstractVector`: amount of each component (mol)
- `volume`: volume of the mixture (m³)
- `RT`: thermal energy (J/mol)

# Keywords

- `buf::Union{NamedTuple, AbstractDict, BrusilovskyThermoBuffer}`: buffers for intermediate
    calculations (see [`pressure`](@ref))

# Returns

- `AbstractVector`: the logarithms of activity coefficients of the components at given
    number of moles, volume and temperature

See also: [`log_c_activity!`](@ref), [`log_c_activity_wj`](@ref), [`log_c_activity_wj!`](@ref)
"""
function log_c_activity(
    mix::MBWREoSMixture,
    nmol::AbstractVector,
    volume::Real,
    RT::Real;
    buf = thermo_buffer(mix, nmol),
)
    Aexc(nmol) = exc_helmholtz(mix, nmol, volume, RT; buf)
    log_ca = ForwardDiff.gradient(Aexc, nmol)
    log_ca ./= RT

    return log_ca
end

"""
    log_c_activity!(log_ca, mixture, nmol, volume, RT[; buf])

Return vector of ln(c_a) - logarithm of activity coefficients
for components of `mixture` at given `nmol`, `volume`, `RT`.
If buffers are provided as keyword arguments, their contents
are modified during the intermediate calculations. The answer is stored
in `log_ca`.

# Arguments

- `log_ca::AbstractVector`: buffer to store the result
- `mix::BrusilovskyEoSMixture`: mixture
- `nmol::AbstractVector`: amount of each component (mol)
- `volume`: volume of the mixture (m³)
- `RT`: thermal energy (J/mol)

# Keywords

- `buf::Union{BrusilovskyThermoBuffer, NamedTuple, AbstractDict}`: buffers for intermediate
    calculations (see [`pressure`](@ref))

# Returns

- `AbstractVector`: the logarithms of activity coefficients of the components at given
    number of moles, volume and temperature

See (Jirí Mikyska, Abbas Firoozabadi // 10.1002/aic.12387)

See also: [`log_c_activity`](@ref), [`log_c_activity_wj`](@ref), [`log_c_activity_wj!`](@ref)
"""
function log_c_activity!(
    log_ca::AbstractVector,
    mix::MBWREoSMixture,
    nmol::AbstractVector,
    volume,
    RT;
    buf = thermo_buffer(mix, nmol),
)
    Aexc(nmol) = exc_helmholtz(mix, nmol, volume, RT; buf)
    log_ca .= ForwardDiff.gradient(Aexc, nmol) ./ RT

    return log_ca
end

"""
    log_c_activity_wj(mixture, nmol, volume, RT[; buf])

Return vector of ln(c_a) - logarithm of activity coefficient
for components of `mixture` at given `nmol`, `volume`, `RT` -
and the jacobian ∂ln(c_a[i]) / ∂n[j].
If buffers are provided as keyword arguments, their contents
are modified during the intermediate calculations.

# Arguments

- `mix`: mixture
- `nmol::AbstractVector`: amount of each component (mol)
- `volume`: volume of the mixture (m³)
- `RT`: thermal energy (J/mol)

# Keywords

- `buf::BrusilovskyThermoBuffer`: buffers for intermediate calculations
    (see [`thermo_buffer`](@ref))

# Returns

- `Tuple{AbstractVector,AbstractMatrix}`: the logarithms of activity coefficients
of the components at given number of moles, volume and temperature, and the jacobian
matrix ∂ln(c_a[i]) / ∂n[j].

See also: [`log_c_activity`](@ref), [`log_c_activity!`](@ref), [`log_c_activity_wj!`](@ref)
"""
function log_c_activity_wj(
    mix::MBWREoSMixture{T},
    nmol::AbstractVector,
    volume::Real,
    RT::Real;
    buf::MBWRThermoBuffer = thermo_buffer(mix, nmol),
) where {T}
    Aexc(nmol) = exc_helmholtz(mix, nmol, volume, RT; buf)
    log_ca = ForwardDiff.gradient(Aexc, nmol)
    log_ca ./= RT
    jacobian = ForwardDiff.hessian(Aexc, nmol)
    jacobian ./= RT
    return log_ca, jacobian
end

"""
    log_activity_wj!(log_ca, jacobian, mixture, nmol, volume, RT[; buf])

Return vector of ln(c_a) - logarithm of activity coefficient
for components of `mixture` at given `nmol`, `volume`, `RT` -
and the jacobian ∂ln(c_a[i]) / ∂n[j]. The first two arguments get overwritten
by the result.
If buffers are provided as keyword arguments, their contents
are modified during the intermediate calculations.

# Arguments

- `log_ca::AbstractVector`: a vector to store the activity coefficients
- `jacobian::AbstractMatrix`: a matrix to store ∂ln(c_a[i]) / ∂n[j]
- `mix`: mixture
- `nmol::AbstractVector`: amount of each component (mol)
- `volume`: volume of the mixture (m³)
- `RT`: thermal energy (J/mol)

# Keywords

- `buf::BrusilovskyThermoBuffer`: buffers for intermediate calculations
    (see [`thermo_buffer`](@ref))

# Returns

- `Tuple{AbstractVector,AbstractMatrix}`: the logarithms of activity coefficients
of the components at given number of moles, volume and temperature, and the jacobian
matrix ∂ln(c_a[i]) / ∂n[j]. The values are aliases for the first two arguments.

See also: [`log_c_activity`](@ref), [`log_c_activity!`](@ref), [`log_c_activity_wj`](@ref)
"""
function log_c_activity_wj!(
    log_ca::AbstractVector,
    jacobian::AbstractMatrix,
    mix::MBWREoSMixture,
    nmol::AbstractVector,
    volume::Real,
    RT::Real;
    buf::MBWRThermoBuffer=thermo_buffer(mix, nmol),
)
    Aexc(nmol) = exc_helmholtz(mix, nmol, volume, RT; buf)
    log_ca .= ForwardDiff.gradient(Aexc, nmol) ./ RT
    jacobian .= ForwardDiff.hessian(Aexc, nmol) ./ RT
    return log_ca, jacobian
end


text:
# ---
# title: 1588. Sum of All Odd Length Subarrays
# id: problem1588
# author: Indigo
# date: 2021-08-30
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/sum-of-all-odd-length-subarrays/description/>
# hidden: true
# ---
#
# Given an array of positive integers `arr`, calculate the sum of all possible
# odd-length subarrays.
#
# A subarray is a contiguous subsequence of the array.
#
# Return  _the sum of all odd-length subarrays of  _`arr`.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [1,4,2,5,3]
#     Output: 58
#     Explanation: The odd-length subarrays of arr and their sums are:
#     [1] = 1
#     [4] = 4
#     [2] = 2
#     [5] = 5
#     [3] = 3
#     [1,4,2] = 7
#     [4,2,5] = 11
#     [2,5,3] = 10
#     [1,4,2,5,3] = 15
#     If we add all these together we get 1 + 4 + 2 + 5 + 3 + 7 + 11 + 10 + 15 = 58
#
# **Example 2:**
#
#
#
#     Input: arr = [1,2]
#     Output: 3
#     Explanation: There are only 2 subarrays of odd length, [1] and [2]. Their sum is 3.
#
# **Example 3:**
#
#
#
#     Input: arr = [10,11,12]
#     Output: 66
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 100`
#   * `1 <= arr[i] <= 1000`
#
#
## @lc code=start
using LeetCode

function sum_odd_length_subarrays(arr::Vector{Int})::Int
    res = 0
    len = length(arr)
    for i in 1:len
        l, r = i - 1, len - i
        lodd, rodd = (l + 1) ÷ 2, (r + 1) ÷ 2
        leven, reven = l ÷ 2 + 1, r ÷ 2 + 1
        res += arr[i] * (lodd * rodd + leven * reven)
    end
    res
end
## @lc code=end


text:
function build_tree(p::AdaOPSPlanner, b0)
    D = AdaOPSTree(p, b0)
    b = 1
    trial = 1
    start = CPUtime_us()

    Depth = sizehint!(Int[], 2000)
    sol = solver(p)
    while D.u[1]-D.l[1] > sol.epsilon_0 &&
          CPUtime_us()-start < sol.T_max*1e6 &&
          trial <= sol.max_trials
        push!(Depth, explore!(D, 1, p))
        trial += 1
    end
    if (CPUtime_us()-start)*1e-6 > sol.T_max*sol.timeout_warning_threshold
        @warn(@sprintf("Surpass the time limit. The actual runtime is %3.1fs.
        Hyperparameters: delta=%4.2f, m_min=%3d, m_max=%3d, zeta=%4.2f, grid=%s, bounds=%s",
        (CPUtime_us()-start)*1e-6, sol.delta, sol.m_min, sol.m_max, sol.zeta, typeof(sol.grid), typeof(sol.bounds)))
        info_analysis(Dict(:tree=>D, :depth=>Depth))
    end
    return D, Depth
end

function explore!(D::AdaOPSTree, b::Int, p::AdaOPSPlanner)
    sol = solver(p)
    while D.Delta[b] < sol.max_depth
        if isempty(D.children[b]) # a leaf
            Δl, Δu = expand!(D, b, p)
            if backup!(D, b, p, Δl, Δu) || excess_uncertainty(D, b, sol.xi, p) <= 0.0
                break
            end
        end
        b = next_best(D, b, p)
    end
    if D.Delta[b] == sol.max_depth
        backup!(D, b, p, -D.l[b], -D.u[b])
    end

    return D.Delta[b]
end

function backup!(D::AdaOPSTree, b::Int, p::AdaOPSPlanner, Δl::Float64, Δu::Float64)
    D.u[b] += Δu
    D.l[b] += Δl
    disc = discount(p.pomdp)
    best_a_change = false
    while b != 1
        bp = b
        ba = D.parent[bp]
        b = D.ba_parent[ba]

        # Update u
        D.ba_u[ba] += disc * D.obs_prob[bp] * Δu
        largest_u = -Inf
        largest_ba = 0
        for ba in D.children[b]
            if D.ba_u[ba] > largest_u
                largest_u = D.ba_u[ba]
                largest_ba = ba
            end
        end
        if largest_ba != ba
            best_a_change = true
        end
        Δu = largest_u - D.u[b]
        D.u[b] = largest_u

        # Update l
        if Δl != 0.0
            D.ba_l[ba] += disc * D.obs_prob[bp] * Δl
            if D.l[b] < D.ba_l[ba]
                Δl = D.ba_l[ba] - D.l[b]
                D.l[b] = D.ba_l[ba]
            else
                Δl = 0.0
            end
        end
    end
    return best_a_change
end

function next_best(D::AdaOPSTree, b::Int, p::AdaOPSPlanner)
    max_u = -Inf
    best_ba = first(D.children[b])
    @inbounds for ba in D.children[b]
        if D.ba_u[ba] > max_u
            max_u = D.ba_u[ba]
            best_ba = ba
        end
    end

    max_eu = -Inf
    best_bp = first(D.ba_children[best_ba])
    tolerated_gap = solver(p).xi * max(D.u[1]-D.l[1], 0.0) / p.discounts[D.Delta[best_bp]+1]
    @inbounds for bp in D.ba_children[best_ba]
        eu = D.obs_prob[bp] * (D.u[bp]-D.l[bp] - tolerated_gap)
        if eu > max_eu
            max_eu = eu
            best_bp = bp
        end
    end

    return best_bp
end

function excess_uncertainty(D::AdaOPSTree, b::Int, ξ::Float64, p::AdaOPSPlanner)
    return D.obs_prob[b] * (D.u[b]-D.l[b] - ξ * max(D.u[1]-D.l[1], 0.0) / p.discounts[D.Delta[b]+1])
end


text:
using Plots
import Plots: plot
export plot

plot(a::AbstractExpression) = plot(to_expr(a), shape = :circle, fontsize=10, shorten=0.01, axis_buffer=0.15)
plot(a::Equation) = plot(to_expr(a), shape = :circle, fontsize=10, shorten=0.01, axis_buffer=0.15)


text:
mutable struct ExaTronProblem{VI, VD, TM <: AbstractTronMatrix}
    n::Cint                 # number of variables
    nnz::Integer            # number of Hessian entries
    nnz_a::Integer          # number of Hessian entries in the strict lower
    A::TM
    B::TM
    L::TM
    indfree::VI   # a working array of dimension n
    iwa::VI       # a working array of dimension 3*n
    g::VD         # gradient
    xc::VD        # a working array of dimension n
    s::VD         # a working array of dimension n
    wa::VD        # a working array of dimension 7*n

    x::VD
    x_l::VD
    x_u::VD
    rows::VI
    cols::VI
    values::VD

    eval_f::Function
    eval_grad_f::Function
    eval_h::Function

    options::Dict{String, Any}

    # Statistics
    delta::Cdouble
    f::Cdouble
    gnorm_two::Float64
    gnorm_inf::Float64
    nfev::Integer
    ngev::Integer
    nhev::Integer
    minor_iter::Integer
    status::Symbol

    ExaTronProblem{VI, VD, TM}() where {VI, VD, TM} = new{VI, VD, TM}()
end

function gpnorm(n, x, x_l, x_u, g)
    two_norm = 0.0
    inf_norm = 0.0

    for i=1:n
        if x_l[i] != x_u[i]
            if x[i] == x_l[i]
                val = (min(g[i], 0.0))^2
            elseif x[i] == x_u[i]
                val = (max(g[i], 0.0))^2
            else
                val = g[i]^2
            end

            two_norm += val
            val = sqrt(val)
            if inf_norm < val
                inf_norm = val
            end
        end
    end

    two_norm = sqrt(two_norm)

    return two_norm, inf_norm
end

function set_default_options!(prob::ExaTronProblem)
    prob.options = Dict{String,Any}()
    set_default_options!(prob.options)
    return
end

function set_default_options!(options::Dict{String, Any})
    options["max_feval"] = 500
    options["max_minor"] = 200
    options["p"] = 5
    options["verbose"] = 0
    options["tol"] = 1e-6
    options["fatol"] = 0.0
    options["frtol"] = 1e-12
    options["fmin"] = -1e32
    options["cgtol"] = 0.1
    options["tron_code"] = :Julia
    options["matrix_type"] = :Sparse

    return
end

function getOption(prob::ExaTronProblem, keyword::String)
    if !haskey(prob.options, keyword)
        error("ExaTron does not have option with name $(keyword).")
    end
    return prob.options[keyword]
end

function setOption(prob::ExaTronProblem, keyword::String, value::Symbol)
    if !haskey(prob.options, keyword)
        error("ExaTron does not have option with name $(keyword).")
    end
    prob.options[keyword] = value
end

function setOption(prob::ExaTronProblem, keyword::String, value::Integer)
    if !haskey(prob.options, keyword)
        error("ExaTron does not have option with name $(keyword).")
    end
    prob.options[keyword] = value
    return
end

function setOption(prob::ExaTronProblem, keyword::String, value::Float64)
    if !haskey(prob.options, keyword)
        error("ExaTron does not have option with name $(keyword).")
    end
    prob.options[keyword] = value
    return
end

function instantiate_memory!(tron::ExaTronProblem{VI,VD}, n, nele_hess) where {VI,VD}
    tron.n = convert(Cint, n)
    tron.nnz = nele_hess
    tron.indfree = VI(undef, n)
    tron.iwa = VI(undef, 3*n)

    tron.g = tron_zeros(VD, n)
    tron.xc = tron_zeros(VD, n)
    tron.s = tron_zeros(VD, n)
    tron.wa = tron_zeros(VD, 7*n)

    tron.x = tron_zeros(VD, n)
    tron.x_l = tron_zeros(VD, n)
    tron.x_u = tron_zeros(VD, n)
    tron.rows = VI(undef, nele_hess)
    tron.cols = VI(undef, nele_hess)
    tron.values = tron_zeros(VD, nele_hess)
end

function createProblem(n::Integer, x_l::AbstractVector{Float64}, x_u::AbstractVector{Float64},
                       nele_hess::Integer, eval_f_cb, eval_grad_f_cb, eval_h_cb; options...)
    @assert n == length(x_l) == length(x_u)
    @assert typeof(x_l) == typeof(x_u)

    # Load options first, as we need to know :matrix_type for type inference
    options_dict = Dict{String, Any}()
    set_default_options!(options_dict)
    for (name, value) in options
        options_dict[string(name)] = value
    end

    VI = Vector{Cint}
    VD = typeof(x_l)

    if options_dict["matrix_type"] == :Sparse
        TM = TronSparseMatrixCSC{VI, VD}
    else
        AT = isa(x_l, Array) ? Array{Float64, 2} : CuArray{Float64, 2}
        TM = TronDenseMatrix{AT}
    end

    tron = ExaTronProblem{VI, VD, TM}()
    tron.options = options_dict
    instantiate_memory!(tron, n, Int64(nele_hess))
    copyto!(tron.x_l, 1, x_l, 1, n)
    copyto!(tron.x_u, 1, x_u, 1, n)

    tron.eval_f = eval_f_cb
    tron.eval_grad_f = eval_grad_f_cb
    tron.eval_h = eval_h_cb

    tron.eval_h(tron.x, :Structure, tron.rows, tron.cols, 1.0, Float64[], Float64[])
    # Instantiate sparse matrix
    p = tron.options["p"]

    if tron.options["matrix_type"] == :Sparse
        tron.A = TronSparseMatrixCSC(tron.rows, tron.cols, tron.values, n)
        tron.B = TronSparseMatrixCSC{VI, VD}(n, Int64(nele_hess))
        tron.L = TronSparseMatrixCSC{VI, VD}(n, Int64(nele_hess + n*p))
        tron.nnz_a = tron.A.nnz
    else
        tron.A = TronDenseMatrix(tron.rows, tron.cols, tron.values, n)
        if isa(x_l, Array)
            tron.B = TronDenseMatrix{Array{Float64, 2}}(n)
            tron.L = TronDenseMatrix{Array{Float64, 2}}(n)
        else
            tron.B = TronDenseMatrix{CuArray{Float64, 2}}(n)
            tron.L = TronDenseMatrix{CuArray{Float64, 2}}(n)
        end
        tron.nnz_a = n*n
    end
    tron.status = :NotSolved

    return tron
end

function solveProblem(tron::ExaTronProblem)
    task = Vector{UInt8}(undef, 60)
    for (i,s) in enumerate("START")
        task[i] = UInt8(s)
    end

    VI = Vector{Cint}
    VD = typeof(tron.x)

    isave = VI(undef, 3)
    dsave = tron_zeros(VD, 3)
    fatol = tron.options["fatol"]::Float64
    frtol = tron.options["frtol"]::Float64
    fmin = tron.options["fmin"]::Float64
    cgtol = tron.options["cgtol"]::Float64
    gtol = tron.options["tol"]::Float64
    max_feval = tron.options["max_feval"]::Int
    tcode = tron.options["tron_code"]::Symbol
    max_minor = tron.options["max_minor"]::Int

    # Project x into its bound. Tron requires x to be feasible.
    tron.x .= max.(tron.x_l, min.(tron.x_u, tron.x))

    itermax = tron.n
    tron.minor_iter = 0
    tron.nfev = tron.ngev = tron.nhev = 0
    tron.status = :NotSolved
    search = true

    while (search)

        # Evaluate the function.

        if Char(task[1]) == 'F' || unsafe_string(pointer(task), 5) == "START"
            tron.f = tron.eval_f(tron.x)
            tron.nfev += 1
            if tron.nfev >= max_feval
                tron.status = :Maximum_Iterations_Exceeded
                search = false
            end
        end

        # Evaluate the gradient and Hessian.

        if (Char(task[1]) == 'G' && Char(task[2]) == 'H') || unsafe_string(pointer(task), 5) == "START"
            tron.eval_grad_f(tron.x, tron.g)
            tron.eval_h(tron.x, :Values, Int[], Int[], 1.0, Float64[], tron.values)
            tron.ngev += 1
            tron.nhev += 1
            tron.minor_iter += 1

            # Copy values in the CSC matrix.
            fill!(tron.A, 0.0)
            transfer!(tron.A, tron.rows, tron.cols, tron.values, tron.nnz)
        end

        # Initialize the trust region bound.

        if unsafe_string(pointer(task), 5) == "START"
            gnorm0 = norm(tron.g)
            tron.delta = gnorm0
        end

        # Call Tron.

        if search
            tron.delta = ExaTronKernels.dtron(tron.n, tron.x, tron.x_l, tron.x_u,
                            tron.f, tron.g, tron.A,
                            frtol, fatol,
                            fmin, cgtol, itermax, tron.delta,
                            task, tron.B, tron.L,
                            tron.xc, tron.s, tron.indfree,
                            isave, dsave, tron.wa, tron.iwa)
        end

        if unsafe_string(pointer(task), 4) == "NEWX"
            tron.gnorm_two, tron.gnorm_inf = gpnorm(tron.n, tron.x, tron.x_l, tron.x_u, tron.g)
            if tron.gnorm_inf <= gtol
                for (i,s) in enumerate("CONVERGENCE: GTOL TEST SATISFIED")
                    task[i] = UInt8(s)
                end
            end

            if tron.minor_iter >= max_minor
                tron.status = :Maximum_Iterations_Exceeded
                search = false
            end
        end

        if unsafe_string(pointer(task), 4) == "CONV"
            tron.status = :Solve_Succeeded
            search = false
        end

        if unsafe_string(pointer(task), 4) == "WARN"
            tron.status = :TooSmall
            search = false
        end
    end

    return tron.status
end


text:
include("tools.jl")
using Main.tools
using LinearAlgebra

n = 200
println("Using n=",n)
xi = get_xi(n)

# for JIT Compiler to do the work
rhs = get_rhs_norm(50)
rhs = get_rhs_sin(50)
a1, s1 = hosvd(rhs, 1e-4, true)

println()
println("First example with sinus")

rhs = get_rhs_sin(xi)
println("Starting hosvd")
# tic()
elapsed_time = @elapsed a1, s1 = hosvd(rhs, 1e-4, true)
# toc()
println("Elapsed time: ", elapsed_time)
println("The core tensor has size ",size(a1.core))
residual = rhs - fullten(a1)
err = norm(residual[:])
println("The error ||rhs - full(a1)||_F is: ", err)
println("10⁻⁴ ⋅ ||B₁|| ", 1e-4*norm(rhs[:]))
println("Relative error is: ", err/norm(rhs[:]))

println()
println("Second example with norm")

rhs = get_rhs_norm(xi)
println("Starting hosvd")
# tic()
elapsed_time = @elapsed a2, s2 = hosvd(rhs, 1e-4, true)
# toc()
println("Elapsed time: ", elapsed_time)
println("The core tensor has size ", size(a2.core))
residual = rhs - fullten(a2)
err = norm(residual[:])
println("The error ||rhs - full(a2)||_F is: ", err)
println("10⁻⁴ ⋅ ||B₂|| ", 1e-4*norm(rhs[:]))
println("Relative Error is: ", err/norm(rhs[:]))

### Plotting
# println("creating Plots")
# ENV["GKSwstype"]="png"
# using Plots;

# p1 = scatter(s1[1],  yaxis=(:log))
# # savefig(p1, "report/eigvalb1.png")
# savefig(p1, "eigvalb1.png")

# p2 = scatter(s2[1], yaxis=(:log))
# # savefig(p2, "report/eigvalb2.png")
# savefig(p2, "eigvalb2.png")


text:
using InfiniteLinearAlgebra, DSP, FillArrays, Test

@testset "conv" begin
    @test conv(1:∞, [2]) ≡ conv([2], 1:∞) ≡ 2:2:∞
    @test conv(1:2:∞, [2]) ≡ conv([2], 1:2:∞) ≡ 2:4:∞
    @test conv(1:∞, Ones(∞))[1:5] == conv(Ones(∞),1:∞)[1:5] == [1,3,6,10,15]
    @test conv(Ones(∞), Ones(∞)) ≡ 1.0:1.0:∞
    @test conv(Ones{Int}(∞), Ones{Int}(∞)) ≡ oneto(∞)
    @test conv(Ones{Bool}(∞), Ones{Bool}(∞)) ≡ oneto(∞)
    @test conv(Fill{Int}(2,∞), Fill{Int}(1,∞)) ≡ conv(Fill{Int}(2,∞), Ones{Int}(∞)) ≡
                conv(Ones{Int}(∞), Fill{Int}(2,∞)) ≡ 2:2:∞
    @test conv(Ones{Int}(∞), [1,5,8])[1:10] == conv([1,5,8], Ones{Int}(∞))[1:10] ==
                    conv(fill(1,100),[1,5,8])[1:10] == conv([1,5,8], fill(1,100))[1:10]
    @test conv(Ones{Int}(∞), 1:4)[1:10] == conv(1:4, Ones{Int}(∞))[1:10] ==
                    conv(fill(1,100),collect(1:4))[1:10] == conv(collect(1:4), fill(1,100))[1:10]
end


text:
using Base64

_open_device(format::Symbol, plot_file; kargs...) = rcall(format, plot_file; kargs...)

_close_device() = rcall(Symbol("dev.off"))

_DEFAULT_FORMAT = :png

function _open_temp_plot(format)
	temp_plot_file = tempname()
	if format == :svg
		_open_device(format, temp_plot_file)
	elseif format == :png
		_open_device(format, temp_plot_file, width=672, height=672, units="px", res=100)
	else
		throw(ErrorException(
				"`$format` is not supported, try with `png` or `svg`"
				))
	end
	temp_plot_file
end

_open_temp_plot() = _open_temp_plot(_DEFAULT_FORMAT)

function _close_and_show_temp_plot(temp_plot_file, format)
	_close_device()
	if isfile(temp_plot_file)
		if format == :png
			plot_base64 = Base64.base64encode(read(temp_plot_file))
			HTML("<img src='data:image/png;base64, $plot_base64'/>")
		elseif format == :svg
			svg_str = read(temp_plot_file, String)
			HTML(_fix_svg_ids(svg_str))
		end
	end
end

function _close_and_show_temp_plot(temp_plot_file)
	_close_and_show_temp_plot(temp_plot_file, _DEFAULT_FORMAT)
end

"""
It executes the R code and returns an `HTML` object containing the plot.
This macro is handy to inline R plots inside Pluto notebooks.

`@html_plot` can take two arguments; the first is the code block to evaluate
to produce the R plot. The second argument is the plot format (optional);
it could be `png` (the default) or `svg`. For example:

```julia
@html_plot R"plot(c(1,2,5,3,4), type='o', col='blue')"
```

To use the `svg` format:

```julia
@html_plot R"plot(c(1,2,5,3,4), type='o', col='blue')" svg
```

**NOTE:** To see a *ggplot2* plot, it is necessary to `print` it explicitly:

```julia
@html_plot R\"""
library(ggplot2)

plt <- ggplot(data = diamonds) +
  geom_bar(
    mapping = aes(x = cut, fill = clarity),
    position = "fill"
  )

print(plt)
\"""
```

"""
macro html_plot(r_code, format...)
	quote
		temp_plot_file = _open_temp_plot($format...)
		$r_code
		_close_and_show_temp_plot(temp_plot_file, $format...)
	end
end


text:
ACTIVATION_FUNCTIONS = [σ, relu, leakyrelu, elu, swish, selu, softplus, softsign];

function test_value_float_precision_preserving(a)
    @testset "$(a): " begin
        for T in [Float32, Float64]
            for val in [-10, -1, 0, 1, 10]
                val = @inferred a(T(val))
                @test typeof(val) == T
            end
        end
    end
end

function test_value_int_input_forces_float64(a)
    @testset "$(a): " begin
        for T in [Int32, Int64]
            for val in [-10, -1, 0, 1, 10]
                val = @inferred a(T(val))
                @test typeof(val) == Float64
            end
        end
    end
end

if Base.find_in_path("ForwardDiff") ≠ nothing
    using ForwardDiff
    function test_value_duals(a)
        @testset "$(a): " begin
        for T in [Float32, Float64, Int32, Int64]
          val = @inferred a(ForwardDiff.Dual(float(T(1)), one(float(T))))
          @test typeof(val) == ForwardDiff.Dual{Void,float(T),1}
        end
        end
    end

    test_value_duals.(ACTIVATION_FUNCTIONS)
end

@testset "Activation Functions" begin

  @test σ(0.0) == 0.5
  @test relu(0.0) == 0.0
  @test leakyrelu(0.0) == 0.0
  @test elu(0.0) == 0.0
  @test swish(0.0) == 0.0
  @test softplus(0.0) ≈ log(2.0)
  @test softsign(0.0) == 0.0
  @test selu(0.0) == 0.0

  @test σ(1.0) == 1.0 / (1.0 + exp(-1.0))
  @test relu(1.0) == 1.0
  @test leakyrelu(1.0) == 1.0
  @test elu(1.0) == 1.0
  @test swish(1.0) == 1.0 / (1.0 + exp(-1.0))
  @test softplus(1.0) ≈ log(exp(1.0) + 1.0)
  @test softsign(1.0) == 0.5
  @test selu(1.0) == 1.0507009873554804934193349852946

  @test σ(-1.0) == 1.0 / (1.0 + exp(1.0))
  @test relu(-1.0) == 0.0
  @test leakyrelu(-1.0) == -0.01
  @test elu(-1.0) == exp(-1.0) - 1.0
  @test swish(-1.0) == -1.0 / (1.0 + exp(1.0))
  @test softplus(-1.0) ≈ log(exp(-1.0) + 1.0)
  @test softsign(-1.0) == -0.5
  @test selu(-1.0) == 1.0507009873554804934193349852946 * 1.6732632423543772848170429916717 * (exp(-1.0) - 1.0)

  @testset "Float inference" begin
    test_value_float_precision_preserving.(ACTIVATION_FUNCTIONS)
  end

  @testset "Test Integer64 and Integer32 inputs will force Float64 outputs" begin
    test_value_int_input_forces_float64.(filter(x -> x != relu, ACTIVATION_FUNCTIONS))

    @testset "relu: " begin
      # relu doesn't have to force floating point outputs
      @test typeof(relu(Int64(1))) == Int64
      @test typeof(relu(Int32(1))) == Int32
    end
  end

  xs = rand(5,5)

  @test all(sum(softmax(xs), 1) .≈ 1)

  @test sum(softmax(vec(xs))) ≈ 1

  @testset "elu" begin
      @test elu(42) == 42
      @test elu(42.) == 42.

      @test elu(-4) ≈ (exp(-4) - 1)
  end

  @test leakyrelu( 0.4,0.3) ≈  0.4
  @test leakyrelu(-0.4,0.3) ≈ -0.12

end


text:
using Base64

# This is adapted from the Documenter.jl 'git_push' function.
function update_homepage(repo, branch, target)
    println("Updating homepage...")

    dirname = ""

    mktempdir() do temp
        # Get current directory.
        root = Base.source_dir()
        root === nothing ? pwd() : root

        # Other variables.
        sha = cd(root) do
            # We'll make sure we run the git commands in the source directory (root), in case
            # the working directory has been changed (e.g. if the makedocs' build argument is
            # outside root).
            try
                readchomp(`git rev-parse --short HEAD`)
            catch
                # git rev-parse will throw an error and return code 128 if it is not being
                # run in a git repository, which will make run/readchomp throw an exception.
                # We'll assume that if readchomp fails it is due to this and set the sha
                # variable accordingly.
                "(not-git-repo)"
            end
        end

        key = ""
        target_dir = ""
        upstream = ""
        keyfile = ""

        cd(root) do
            # Grab the push key.
            key = get(ENV, "DOCUMENTER_KEY", "")

            dirname = isempty(dirname) ? temp : joinpath(temp, dirname)
            isdir(dirname) || mkpath(dirname)

            keyfile = abspath(joinpath(root, ".documenter"))
            target_dir = abspath(target)

            # The upstream URL to which we push new content and the ssh decryption commands.
            upstream = "git@$(replace(repo, "github.com/" => "github.com:"))"

            write(keyfile, String(base64decode(key)))
            chmod(keyfile, 0o600)
        end

        try
            # Use a custom SSH config file to avoid overwriting the default user config.
            withfile(joinpath(homedir(), ".ssh", "config"),
                """
                Host github.com
                    StrictHostKeyChecking no
                    HostName github.com
                    IdentityFile $keyfile
                """
            ) do
                # cd(temp) do
                # end

                cd(temp) do
                    println("Setting up git...")
                    # Setup git.
                    run(`git init`)
                    run(`git config user.name "autodocs"`)
                    run(`git config user.email "autodocs"`)

                    # Fetch from remote and checkout the branch.
                    run(`git remote add upstream $upstream`)
                    run(`git fetch upstream`)

                    try
                        run(`git checkout -b $branch upstream/$branch`)
                    catch e
                        run(`git checkout --orphan $branch`)
                        run(`git commit --allow-empty -m "Initial empty commit for homepage"`)
                    end

                    # Copy the target directory to the specified directory.
                    filecopy(target_dir, dirname)

                    # Add, commit, and push the docs to the remote.
                    run(`git add -A .`)
                    if !success(`git diff --cached --exit-code`)
                        run(`git commit -m "homepage update based on $sha"`)
                        run(`git push -q upstream HEAD:$branch`)
                    end
                end
            end
        catch e
            println("Errored out building homepage, $e")
        finally
            # Remove the unencrypted private key.
            isfile(keyfile) && rm(keyfile)
        end
    end
end

function withfile(func, file::AbstractString, contents::AbstractString)
    hasfile = isfile(file)
    original = hasfile ? read(file, String) : ""
    open(file, "w") do stream
        print(stream, contents)
        flush(stream) # Make sure file is written before continuing.
    end
    try
        func()
    finally
        if hasfile
            open(file, "w") do stream
                print(stream, original)
            end
        else
            #rm(file)
        end
    end
end

function filecopy(src, dst, force_copy = true)
    for item in readdir(src)
        cp(joinpath(src, item), joinpath(dst, item), force = force_copy)
    end
end

function filecopy_deep(src, dst)
    for (root, dirs, files) in walkdir(src)
        for file in files
            source = joinpath(root, file)
            target = replace(source, src => dst)
            isdir(dirname(target)) || mkpath(dirname(target))
            cp(source, target, force=true)
        end
    end
end

function preprocess_markdown(folder)
    yaml_dict = Dict()

    try
        for (root, dirs, files) in walkdir(folder)
            for file in files
				if endswith(file, ".md")
	                full_path = joinpath(root, file)
	                yaml_dict[full_path] = remove_yaml(full_path)
				end
            end
        end
    catch e
        # println("Markdown copy error: $e")
        rethrow(e)
    end

    return yaml_dict
end

function postprocess_markdown(folder, yaml_dict; original = "")
    try
        for (root, dirs, files) in walkdir(folder)
            for file in files
                full_path = joinpath(root, file)
                original_path = full_path

                if length(original) > 0
                    # original_path = abspath(original, file)
                    original_path = replace(full_path, folder => original)
                    # original_path = replace(
                    #     full_path,
                    #     joinpath("docs", "site") => joinpath("docs", "src")
                    # )
                end

                if haskey(yaml_dict, original_path)
                    # println("Original: $original_path => Full path: $full_path")

                    txt = open(f -> read(f, String), full_path)
                    open(full_path, "w+") do f
                        # Add in the yaml block.
                        for line in yaml_dict[original_path]
                            write(f, line)
                        end

						txt = replace(txt, "api.md" => "{{site.baseurl}}/docs/library/")

                        # Add the rest of the text.
						if original_path == full_path
							write(f, txt)
						else
							write(f, replace(txt, "![](figures" => "![](/{{site.baseurl}}/tutorials/figures"))
						end
                    end
                elseif endswith(file, ".md")
                    println("Original: $original_path")
                    println("Full:     $full_path \n")
                end

				# Make specific api items headers.
				if file == "api.md" && original_path != full_path
					tidy_api(full_path)
				end
            end
        end
    catch e
        # println("Markdown copy error: $e")
        rethrow(e)
    end

    return yaml_dict
end

function remove_yaml(file)
    lines = readlines(file, keep=true)
    yaml = []

    yaml_start = 1
    yaml_stop = length(lines)

    # Read out YAML block.
    if replace(lines[1], "\n" => "") == "---"
        for i = 2:length(lines)
            if replace(lines[i], "\n" => "") == "---"
                yaml = lines[1:i]
                yaml_stop = i
                break
            end
        end

        # Write non-YAML back to file.
        open(file, "w+") do f
            for line in lines[yaml_stop + 1 : end]
                write(f, line)
            end
        end
    end

    return yaml
end

function tidy_api(file)
	lines = readlines(file, keep=true)

    # Find the ID sections.
    for i = 1:(length(lines)-1)
		first = lines[i]
		second = lines[i+1]
        if startswith(first, "<a id=") && startswith(second, "**")
            first = replace(first, "\n" => "")
			second = replace(second, "\n" => "")
			# final_line = replace(first, ">" => ">$second", count = 1)

			lines[i] = "### $first $second"
			lines[i+1] = ""
        end
    end

    # Write lines back.
    open(file, "w+") do f
        for line in lines
            write(f, line)
        end
    end
end


text:
# This file is a part of NeighborhoodApproximationIndex.jl
import SimilaritySearch: search
using InvertedFiles, Intersections, KCenters, StatsBase, Parameters, LinearAlgebra
export DeloneInvertedFile, search

const GlobalEncodeKnnResult = [KnnResult(10)]
const GlobalEncodeVector = [SVEC()]

function __init__()
    for i in 2:Threads.nthreads()
        push!(GlobalEncodeKnnResult, KnnResult(10))
        push!(GlobalEncodeVector, SVEC())
    end
end

@inline getencodeknnresult() = @inbounds GlobalEncodeKnnResult[Threads.threadid()]
@inline getencodevector() = @inbounds GlobalEncodeVector[Threads.threadid()]

@with_kw struct DeloneInvertedFile{
        DistType<:PreMetric,
        DataType<:AbstractVector,
        InvertedFileType<:InvertedFile,
        CentersType<:AbstractSearchContext
        } <: AbstractSearchContext
    dist::DistType = SqL2Distance()
    db::DataType
    centers::CentersType
    invfile::InvertedFileType
    k::Int32 = 7
    t::Int32 = 1
end

@inline Base.getindex(D::DeloneInvertedFile, i) = @inbounds D.db[i]

Base.copy(D::DeloneInvertedFile;
    dist=D.dist,
    db=D.db,
    centers=D.centers,
    invfile=D.invfile,
    k=D.k,
    t=D.t
    ) = DeloneInvertedFile(; dist, db, centers, invfile, k, t)

function encode!(D::DeloneInvertedFile, obj::T, v::DVEC) where T
    empty!(v)
    encres = getencodeknnresult()
    empty!(encres, D.k)

    for (id_, d_) in search(D.centers, obj, encres)
        v[id_] = d_
    end

    normalize!(v)
end

function Base.push!(D::DeloneInvertedFile, obj, v::DVEC)
    push!(D.db, obj)
    push!(D.invfile, length(D.db) => v)
end

function Base.append!(D::DeloneInvertedFile, db; parallel_block=1, encode=true, verbose=true)
    encode = encode || !(db[1] isa DVEC)

    if !encode
        for v in db
            push!(D, v, v)
        end
    elseif parallel_block > 1
        nt = Threads.nthreads()
        sp = 1
        n = length(db)
        venc = [SVEC() for i in 1:parallel_block]
        while sp < n
            ep = min(n, sp + parallel_block - 1)
            verbose && println(stderr, "$(typeof(D)) appending chunk ", (sp=sp, ep=ep, n=n), " ", Dates.now())

            Threads.@threads for i in sp:ep
                @inbounds encode!(D, db[i], venc[i-sp+1])
            end

            for i in sp:ep
                @inbounds push!(D, db[i], venc[i-sp+1])
            end

            sp = ep + 1
        end
    else
        for v in db
            push!(D, v, encode!(D, v, getencodevector()))
        end
    end

    D
end

"""
    DeloneInvertedFile(
        dist::PreMetric,
        db::AbstractVector;
        numcenters=ceil(Int, sqrt(length(db))),
        refs=:rand,
        centers=:delone,
        k=5,
        t=1,
        parallel_block=Threads.nthreads() > 1 ? 1024 : 1,
        verbose=true
    )

Creates a `DeloneInvertedFile`, high level interface.

- `dist`: Distance object (a `PreMetric` object, see `Distances.jl`)
- `db`: Objects to be indexed
- `centers`: An index on a set of references used to create the index. It can also be a symbol that describes how to create the index. Valid symbols are:
  - `:delone` creates a DeloneInvertedFile
  - `:graph` creates a SearchGraph
  - `:exhaustive` creates an ExhaustiveSearch
  Please note that these indexes will be created with generic parameters
- `numcenters` number of centers to insert, only used if `centers` is not an index
- `refs`: References, it can be raw objects or symbols indicating how to compute them `:rand`, `:dnet`, `:fft` (see `KCenters.jl`). Only used if `centers` is not an index
- `k` the number of references to be used
- `t` the number of occurrences in the inverted file required to accept an object as candidate in the result set
- `parallel_block` Parallel construction works on batches, this is the size of these blocks
- `verbose` true if you want to see messages
"""
function DeloneInvertedFile(
        dist::PreMetric,
        db::AbstractVector;
        numcenters=ceil(Int, sqrt(length(db))),
        refs=:rand,
        centers=:delone,
        k=5,
        t=1,
        parallel_block=Threads.nthreads() > 1 ? 1024 : 1,
        verbose=true
    )

    if centers === nothing || centers in (:delone, :graph, :exhaustive)
        if refs === :rand
            refs = db[unique(rand(1:length(db), numcenters))]
        elseif refs in (:dnet, :fft)
            train = db[unique(rand(1:length(db), 3*numcenters))]
            C = kcenters(dist, train, numcenters; initial=refs, maxiters=0)
            refs = C.centers[C.dmax .> 0.0]
        end

        if centers === :delone
            m = ceil(Int, length(refs) / log2(length(refs)))
            centers = DeloneInvertedFile(dist, refs; numcenters=m, centers=:exhaustive, k=3, t=1, verbose)
        elseif centers === :exhaustive
            centers = ExhaustiveSearch(dist, refs)
        elseif centers === :graph
            centers = SearchGraph(; dist, db=eltype(refs)[])
            append!(centers, refs; parallel_block=1024)
        end
    end

    D = DeloneInvertedFile(;
        dist, k, t,
        db=eltype(db)[],
        invfile=InvertedFile(length(centers)),
        centers=centers
    )

    append!(D, db; parallel_block, verbose)
    D
end

"""
    search(D::DeloneInvertedFile, q::T, res::KnnResult)

Searches nearest neighbors of `q` inside the `index` under the distance function `dist`.
"""
function search(D::DeloneInvertedFile, q, res::KnnResult; t=D.t)
    Q = prepare_posting_lists_for_querying(D.invfile, encode!(D, q, getencodevector()))

	umerge(Q, t) do L, P, m
        @inbounds begin
            id = L[1].I[P[1]]
		    push!(res, id, evaluate(D.dist, D.db[id], q))
        end
	end

    res
end

#=
"""
    optimize!(perf::Performance, index::DeloneInvertedFile, recall=0.9; numqueries=128, verbose=false, k=10)

Tries to configure `index` to achieve the specified recall for fetching `k` nearest neighbors.
"""
function optimize!(perf::Performance, index::DeloneInvertedFile, recall=0.9; numqueries=128, verbose=false, k=10)
    verbose && println("$(string(index)) optimizing for recall=$(recall)")
    index.opts.ksearch = 1
    p = probe(perf, index)

    while p.macrorecall < recall && index.opts.ksearch < length(index.lists)
        index.opts.ksearch += 1
        verbose && println(stderr, "$(string(index)) optimize! step ksearch=$(index.opts.ksearch), performance ", p)
        p = probe(perf, index)
    end

    verbose && println("$(string(index)) reached performance ", p)
    index
end=#


text:
import Term.Boxes: ROUNDED, get_row, get_title_row, fit
using Term: Term
println("get row")
@time get_row(ROUNDED, [3], :top)
@time get_row(ROUNDED, 3, :top)

println("get title row")
@time get_title_row(:top, Term.box.ROUNDED, nothing, width = 12, style = "red")

for justify in (:left, :center, :right)
    @time get_title_row(
        :top,
        Term.box.ROUNDED,
        "test",
        width = 12,
        justify = justify,
        title_style = "blue",
        style = "red",
    )
end


text:
using Compose, Gadfly, RDatasets

set_default_plot_size(21cm, 16cm)

mpg = dataset("ggplot2","mpg")
p1 = plot(mpg, x=:Class, y=:Hwy, color=:Class, Geom.boxplot, Theme(boxplot_spacing=0.3cx) )

p2 = plot(mpg, x=:Class, y=:Hwy, color=:Drv, Geom.boxplot,
    Guide.colorkey(labels=["front", "4-wheel\t\t","rear"]),
    Theme(boxplot_spacing=0.1cx) )

vstack(p1, p2)


text:
using Oceananigans.Utils: prettytime
using Oceananigans: short_show

"""Show the innards of a `Model` in the REPL."""
Base.show(io::IO, model::ShallowWaterModel{G, A, T}) where {G, A, T} =
    print(io, "ShallowWaterModel{$(Base.typename(A)), $T}",
        "(time = $(prettytime(model.clock.time)), iteration = $(model.clock.iteration)) \n",
        "├── grid: $(short_show(model.grid))\n",
        "├── tracers: $(tracernames(model.tracers))\n",
        "└── coriolis: $(typeof(model.coriolis))")


text:
export Up, Down, Left, Right
export UP, DOWN, LEFT, RIGHT, LRUD

#####
# Directions
#####

struct Up end
const UP = Up()
(x::Up)(p::CartesianIndex{2}) = p + CartesianIndex(-1, 0)

struct Down end
const DOWN = Down()
(x::Down)(p::CartesianIndex{2}) = p + CartesianIndex(1, 0)

struct Left end
const LEFT = Left()
(x::Left)(p::CartesianIndex{2}) = p + CartesianIndex(0, -1)

struct Right end
const RIGHT = Right()
(x::Right)(p::CartesianIndex{2}) = p + CartesianIndex(0, 1)

const LRUD = Union{Left, Right, Up, Down}


text:
using DataFrames
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = page_model()
    include("../src/components/NonMarketDamages.jl")

    nonmarketdamages = addnonmarketdamages(m)

    set_param!(m, :NonMarketDamages, :rtl_realizedtemperature, readpagedata(m, "test/validationdata/$valdir/rtl_realizedtemperature.csv"))
    set_param!(m, :NonMarketDamages, :rcons_per_cap_MarketRemainConsumption, readpagedata(m, "test/validationdata/$valdir/rcons_per_cap_MarketRemainConsumption.csv"))
    set_param!(m, :NonMarketDamages, :rgdp_per_cap_MarketRemainGDP, readpagedata(m, "test/validationdata/$valdir/rgdp_per_cap_MarketRemainGDP.csv"))
    set_param!(m, :NonMarketDamages, :atl_adjustedtolerableleveloftemprise, readpagedata(m, "test/validationdata/atl_adjustedtolerableleveloftemprise_nonmarket.csv"))
    set_param!(m, :NonMarketDamages, :imp_actualreduction, readpagedata(m, "test/validationdata/$valdir/imp_actualreduction_nonmarket.csv"))
    set_param!(m, :NonMarketDamages, :isatg_impactfxnsaturation, 17.0)

    p = load_parameters(m)
    p["y_year_0"] = 2015.
    p["y_year"] = Mimi.dim_keys(m.md, :time)
    set_leftover_params!(m, p)

    run(m)
    rcons_per_cap = m[:NonMarketDamages, :rcons_per_cap_NonMarketRemainConsumption]
    rcons_per_cap_compare = readpagedata(m, "test/validationdata/$valdir/rcons_per_cap_NonMarketRemainConsumption.csv")
    @test rcons_per_cap ≈ rcons_per_cap_compare rtol=1e-2

    rgdp_per_cap = m[:NonMarketDamages, :rgdp_per_cap_NonMarketRemainGDP]
    rgdp_per_cap_compare = readpagedata(m, "test/validationdata/$valdir/rgdp_per_cap_NonMarketRemainGDP.csv")
    @test rgdp_per_cap ≈ rgdp_per_cap_compare rtol=1e-2
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE