text:
# file for functions related to calculating the sparsity pattern of a matrix

# can be generalized with some topology constants
function getSparsityCounts(mesh::PumiMeshDG, sparse_bnds::AbstractArray{Int32, 2}; getdofs=true)
# count how many local, remote nodes/dofs each dof is related to

  numNodesPerElement = mesh.numNodesPerElement
  if getdofs || mesh.numDofPerNode == 1
    numDofPerNode = mesh.numDofPerNode
    div_factor = 1
    add_factor = 0
  else
    numDofPerNode = 1
    div_factor = mesh.numDofPerNode
    add_factor = 1
  end

  edges = Array{Ptr{Void}}(4)
  part_nums = Array{Cint}(4)
  shared_edges = Array{Ptr{Void}}(4)
  for i=1:mesh.numEl
    el_ptr = mesh.elements[i]
    apf.getDownward(mesh.m_ptr, el_ptr, mesh.dim-1, edges)
    nremotes = 0
    for j=1:mesh.numFacesPerElement
      edge_j = edges[j]
      ncopies = apf.countCopies(mesh.shr_ptr, edge_j)
      apf.getCopies(part_nums, shared_edges)
      for k=1:ncopies
        if part_nums[k] != mesh.myrank
          nremotes += 1
        end
      end
    end

    nlocal = mesh.numFacesPerElement + 1  # 3 neighbors + 1
    nlocal = nlocal - nremotes

    # apply this to all dofs on this element
    for j=1:mesh.numNodesPerElement
      for k =1:numDofPerNode
        dof_k = mesh.dofs[k, j, i]
        dof_k = div(dof_k, div_factor) + add_factor
        sparse_bnds[1, dof_k] = nlocal*numNodesPerElement*numDofPerNode
        sparse_bnds[2, dof_k] = nremotes*numNodesPerElement*numDofPerNode
      end
    end
  end

  return nothing
end

# WIP: tighter sparsity bounds for CG
#=
function getDofConnectivity(mesh::PumiMesh2)

  # find the maximum number of elements a dof is related to
  max_el = 0
  up_cnt = 400 # maximum number of upward adjacencies
  up = Array{Ptr{Void}}(up_cnt)  # equivalent of apf::Up
  verts = Array{Ptr{Void}}(3)
  nel = 0
  els_all = Array{Ptr{Void}}(3*up_cnt)
  els_tmp = Array{Ptr{Void}}(3*up_cnt)
  for i=1:mesh.numEl
    apf.getDownward(mesh.m_ptr, 0, verts)
    pos = 1
    for j=1:3
      nel[j] = apf.countAdjacent(mesh.m_ptr, verts[j], 2)
      apf.getAdjacent(up)
      # copy element pointer into a single array
      for k=1:nel[j]
        els_all[pos] = up[k]
        pos += 1
      end
    end
    numel = pos - 1

    # count number of non duplicates
    els_used = sview(els_all, 1:numel);
    sort!(els_used)
    cnt = 1
    for j=2:length(els_used)
      if els_used[j] != els_used[j-1]
        cnt += 1
      end
    end
    if cnt > max_el
      max_el = cnt
    end
  end

  dof_elements = zeros(Int32, max_el, mesh.numNodes)
  # now get the element numbers
  for i=1:mesh.numEl
    apf.getDownward(mesh.m_ptr, 0, verts)
    pos = 1
    for j=1:3
      nel[j] = apf.countAdjacent(mesh.m_ptr, verts[j], 2)
      apf.getAdjacent(up)
      # copy element pointer into a single array
      for k=1:nel[j]
        els_all[pos] = up[k]
        pos += 1
      end
    end
    numel = pos - 1

    # de-duplicate whle copying into node_elements
    els_used = sview(els_all, 1:numel);
    sort!(els_used)
    els_tmp[1] = els_used[1]
    pos = 2
    for j=1:length(els_used)
      if els_used[j] != els_used[j-1]
        els_tmp[pos] = els_used[j]
      end
    end

  end
end
=#

    # copy the element pointers into a single array

# or maybe replace this with a tighter sparsity bound function
# can be generalized with a few if statements
function getSparsityBounds(mesh::PumiMesh, sparse_bnds::AbstractArray{Int32, 2}; getdofs=true)
# sparse_bnds : 2 by numDof array of the minimum, maximum dofs connected to
# a dof in the Jacobian
# or node numbers, if getdofs=false
# this works for high order elements despite not using mesh.elementNodeOffsets
# because it does not reference what element a particular entity belongs to
# getDofBounds also does not need it because only the minimum and maximum
# are important, not the order of them

#resetAllIts2(mesh.m_ptr)
# mesh iterator increment, retreval functions
#iterators_inc = [incrementVertIt, incrementEdgeIt, incrementFaceIt]
#iterators_inc = (VertIterator(), EdgeIterator(), FaceIterator(), ElIterator())

#iterators_get = [getVert, getEdge, getFace]
#iterators_get = (VertGetter(), EdgeGetter(), FaceGetter(), ElGetter)
num_entities = mesh.numEntitiesPerType
#num_entities = [mesh.numVert, mesh.numEdge, mesh.numEl]
num_nodes_entity = mesh.numNodesPerType  # number of nodes on each type
                                         # of mesh entity

if getdofs
  numbering_Nptr = mesh.dofnums_Nptr
  numDofPerNode = mesh.numDofPerNode
  @assert size(sparse_bnds, 2) == mesh.numDof
else
  numbering_Nptr = mesh.nodenums_Nptr
  numDofPerNode = 1
  @assert size(sparse_bnds, 2) == mesh.numNodes
end

for etype=1:(mesh.dim+1)  # loop over mesh entity types

  if (num_nodes_entity[etype] != 0)  # there are nodes herea
    it = apf.MeshIterator(mesh.m_ptr, etype-1)
    for entity = 1:num_entities[etype]  # loop over all entities of this type
#      entity_ptr = iterators_get[etype]()  # get pointer to mesh entity
      entity_ptr = apf.iterate(it)

      min, max = getDofBounds(mesh, entity_ptr, getdofs=getdofs)
      for node = 1:num_nodes_entity[etype]  # loop over nodes on each entity
	# get the minimum, maximum related dof numbers
	# use the same min, max for all dofs on this node

	for dof = 1:numDofPerNode
	  dofnum = apf.getNumberJ(numbering_Nptr, entity_ptr, node - 1, dof-1)
          sparse_bnds[1, dofnum] = min
	  sparse_bnds[2, dofnum] = max

	end  # end loop over dofs
      end  # end loop over nodes on entity
#      iterators_inc[etype]()
    end  # end loops over entities of this type
    apf.free(it)
  end  # end if statement
end  # end loop over entity types

return nothing

end

# rewrite this because it is terrible
# this could be generalized 3d?
function getDofBounds(mesh::PumiMesh2D, entity_ptr::Ptr{Void}; getdofs=true)
# gets the maximum, minimum dofs associated with the entity currently
# pointed to by the iterator specified by etype
# getDofs = true -> get dof numbers, false -> get node numbers
# this works for distance-0 and 1 colorings

#iterators_get = [getVert, getEdge, getFace]
#it = apf.MeshIterator(mesh.m_ptr, etype - 1)
#entity_ptr = iterators_get[etype]()

# get associated elements (distance-0 elements)
num_adj = apf.countAdjacent(mesh.m_ptr, entity_ptr, 2)
el_arr = apf.getAdjacent(num_adj)

# distance-1
#=
dofnums = zeros(Int32, mesh.numDofPerNode, mesh.numNodesPerElement, num_adj)

for i=1:num_adj
  el_i = apf.getNumberJ(mesh.el_Nptr, el_arr[i], 0, 0) + 1
  sub_arr = sub(dofnums, :, :, i)
  getGlobalNodeNumbers(mesh, el_i, sub_arr)
end

min, max = getMinandMax(dofnums)

return min, max
=#

  # get distance-2 elements
  # this is really distance-1
  if mesh.coloringDistance >= 2
    edge_arr = Array{Ptr{Void}}(num_adj*3)  # enough space for all edges, including repeats
    for i=1:num_adj  # get the edges
      sub_arr = sview(edge_arr, (3*(i-1) + 1):(3*i))
      apf.getDownward(mesh.m_ptr, el_arr[i], 1, sub_arr)
    end

    # edge_arr now populated with all edges

    # print the edge numbers
  #=
    for i=1:num_adj*3
      edge_num_i = apf.getNumberJ(mesh.edge_Nptr, edge_arr[i], 0, 0)
      println("edge ", i, " has number ", edge_num_i)
    end
  =#
    # now get the elements the edges belong to
    # count the number of elements first, then get them
    num_els = zeros(Int, length(edge_arr) + 1)  # count number of elements each edge has

  #  println("length(num_els) = ", length(num_els))
  #  println("length(edge_arr) = ", length(edge_arr))
    for i=1:length(edge_arr)
      num_els[i] = apf.countAdjacent(mesh.m_ptr, edge_arr[i], 2)
    end

    # now get the elements
    num_adj = sum(num_els)
    el_arr = Array{Ptr{Void}}(num_adj)  # rebind el_arr to new array
    start_idx = 1
    end_idx = num_els[1]
    for i=1:length(edge_arr)
      edge_i = edge_arr[i]
      sub_arr = sview(el_arr, start_idx:end_idx)
      apf.countAdjacent(mesh.m_ptr, edge_i, 2)
      apf.getAdjacent(sub_arr)

      # update indices
      start_idx = end_idx + 1
      end_idx += num_els[i+1]
    end

    # now el_arr has all the elements, including repeats
    # num_adj = legnth(el_arr)

  end  # end if distance-2

  # get dofnums for the elements
  if getdofs
    numDofPerNode = mesh.numDofPerNode
  else  # get node numbers
    numDofPerNode = 1
  end

  dofnums = zeros(Int32, numDofPerNode, mesh.numNodesPerElement, num_adj)

  for i=1:num_adj
    elnum_i = apf.getNumberJ(mesh.el_Nptr, el_arr[i], 0, 0) + 1
    getGlobalNodeNumbers(mesh, elnum_i, sview(dofnums, :, :, i), getdofs=getdofs)
  end

  min, max = getMinandMax(dofnums)

  return min, max

end


text:
using Distributions

function sample_q(intercept_mu::Float64, slope_mu::Float64, std::Float64, num::Int)
    intercepts = rand(Normal(intercept_mu, std), num)
    slopes = rand(Normal(slope_mu, std), num)
    (intercepts, slopes)
end

function log_q(intercept_mu::Float64, slope_mu::Float64, intercept::Float64, slope::Float64, std::Float64)
    intercept_dist = Normal(intercept_mu, std)
    slope_dist = Normal(slope_mu, std)
    logpdf(intercept_dist, intercept) + logpdf(slope_dist, slope)
end

function log_q_grad(intercept_mu::Float64, slope_mu::Float64, intercept::Float64, slope::Float64, std::Float64)
    intercept_grad = (intercept - intercept_mu) / (std * std)
    slope_grad = (slope - slope_mu) / (std * std)
    [intercept_grad, slope_grad]
end

function log_likelihood(x::Array{Float64,1}, y::Array{Float64,1}, noise_std::Float64, intercept::Float64, slope::Float64)
    line = x * slope + intercept
    ll = 0.0
    for (y_i, line_i) in zip(y, line)
        dist = Normal(line_i, noise_std)
        ll += logpdf(dist, y_i)
    end
    ll
end

function log_p_unnormalized(problem::Problem, intercept::Float64, slope::Float64)
    intercept_dist = Normal(0.0, problem.prior_std)
    slope_dist = Normal(0.0, problem.prior_std)
    log_prior = logpdf(intercept_dist, intercept)
    log_prior += logpdf(slope_dist, slope)
    log_lik = log_likelihood(problem.x, problem.y, problem.noise_std, intercept, slope)
    log_prior + log_lik
end

function gradient_est(params::Array{Float64,1}, problem::Problem, std::Float64, num::Int)
    @assert length(params) == 2
    intercept_mu = params[1]
    slope_mu = params[2]
    (intercepts, slopes) = sample_q(intercept_mu, slope_mu, std, num)
    elbo_est = 0.0
    grad_est = zeros(2)
    for (intercept, slope) in zip(intercepts, slopes)
        lj = log_p_unnormalized(problem, intercept, slope)
        lq = log_q(intercept_mu, slope_mu, intercept, slope, std)
        diff = lj - lq
        elbo_est += diff
        grad_est += diff * log_q_grad(intercept_mu, slope_mu, intercept, slope, std)
    end
    elbo_est /= num
    grad_est /= num
    (grad_est, elbo_est)
end

immutable BBVIResult
    intercept_mu::Float64
    slope_mu::Float64
    elbo_est::Float64
end

function linreg_fast_bbvi_foreign(x::Array{Float64,1}, y::Array{Float64,1}, prior_std::Float64, noise_std::Float64, iters::Int, std::Float64, num::Int, a::Float64, b::Float64)
    # num is the number of samples from q to use to approximate the gradient
    # (a + t)^(-b) is the step size at time t
    params = zeros(2)
    elbo_est = NaN
    problem = Problem(x, y, prior_std, noise_std)
    @assert iters > 0
    start = time() # in seconds
    for iter in 1:iters
        rho = (a + iter) ^ (-b)
        (grad_est, elbo_est) = gradient_est(params, problem, std, num)
        params += rho * grad_est
    end
    intercept_mu = params[1]
    slope_mu = params[2]
    BBVIResult(intercept_mu, slope_mu, elbo_est)
end

immutable BBVIParams
    # std of the intercept and slope normals of the variational family (constant)
    std::Float64

    # Number of samples to use for gradient estimates
    num::Int

    # Robbins Munro parameters
    a::Float64
    b::Float64

    # iterations of gradient ascent (no early stopping)
    iters::Int

    # model parameters
    prior_std::Float64
    noise_std::Float64
end

immutable BBVISampler <: Sampler
    intercept_dist::Normal
    slope_dist::Normal

    # when you construct the BBVI it does the optimization
    function BBVISampler(params::BBVIParams, x::Array{Float64,1}, y::Array{Float64,1})
        result = linreg_fast_bbvi_foreign(x, y, params.prior_std, params.noise_std,
                                          params.iters, params.std, params.num, params.a, params.b)

        intercept_dist = Normal(result.intercept_mu, params.std)
        slope_dist = Normal(result.slope_mu, params.std)
        new(intercept_dist, slope_dist)
    end
end

function simulate(sampler::BBVISampler)
    intercept = rand(sampler.intercept_dist)
    slope = rand(sampler.slope_dist)
    sample = [intercept, slope]
    log_weight = regenerate(sampler, sample)
    (sample, log_weight)
end

function regenerate(sampler::BBVISampler, intercept_and_slope::Array{Float64,1})
    @assert length(intercept_and_slope) == 2
    (intercept, slope) = intercept_and_slope
    logpdf(sampler.intercept_dist, intercept) + logpdf(sampler.slope_dist, slope)
end

ext: jl


text:
## NOTE:
# n-dim sparse array?
# each code word in CodingTheory is a point in an n-dim array

"""
    n_adjacencies(dim::Integer) -> Integer
    n_adjacencies(M::AbstractArray{T, N}) -> Integer
    n_adjacencies(I::AbstractIndexOrIndices{T, N}) -> Integer

Given a matrix or dimension, returns the number of elements adjacent to any given element in an infinite lattice/matrix (i.e., not including edges).  For edges, see `adjacencies`).
"""
n_adjacencies(dim::Integer) = 3^dim - 1
n_adjacencies(M::AbstractArray{T, N}) where {T, N} = n_adjacencies(ndims(M))
n_adjacencies(I::AbstractIndexOrIndices{N}) where {N} = n_adjacencies(length(first(I)))
# n_adjacencies(M::AbstractArray{T, N}) where {T, N} = n_adjacencies(N)

"""
Repeats a specified value n many times along the specified dimension.  If no `fill_elem` is given, and the element type of the matrix is a number, fills with zero.

```julia
julia> A = rand(Int8, 2, 2)
2×2 Array{Int8,2}:
  28  23
 -47  54

julia> append_n_times(A, 2, Int8(3), dims = 1) # repeat the value 3 twice along the first dimension
4×2 Array{Int8,2}:
  28  23
 -47  54
   3   3
   3   3
```
"""
function append_n_times(M::AbstractArray{T, N}, n::Integer, fill_elem::T; dims::Integer = 1) where {T, N}
    sz = ntuple(d -> d == dims ? n : size(M, d), max(N, dims))
    return cat(M, fill(fill_elem, sz); dims = dims)
end
append_n_times(M::AbstractArray{T, N}, n::Integer, dims::Integer = 1) where {T <: Number, N} = append_n_times(M, n, zero(T); dims = dims)

"""
```julia
julia> A = rand(Int8, 2, 2)
2×2 Array{Int8,2}:
  28  23
 -47  54

julia> append_n_times_backwards(A, 2, Int8(3), dims = 1) # repeat the value 3 twice along the first dimension
4×2 Array{Int8,2}:
  3   3
  3   3
  28  23
 -47  54
```
"""
function append_n_times_backwards(M::AbstractArray{T, N}, n::Integer, fill_elem::T; dims::Integer = 1) where {T, N}
    sz = ntuple(d -> d == dims ? n : size(M, d), max(N, dims))
    return cat(fill(fill_elem, sz), M; dims = dims)
end
append_n_times_backwards(M::AbstractArray{T, N}, n::Integer, dims::Integer = 1) where {T <: Number, N} = append_n_times(M, n, zero(T); dims = dims)

"""
    function promote_to_nD(M::AbstractArray{T, N}, n::Integer, fill_elem::T)

Assumes the given matrix M is an (n - 1) dimensional slice of an n-dimensional structure, and fills in the array to n dimensions with `fill_elem`.  If no `fill_elem` is given, and the element type of the matrix is a number, fills with zero.
"""
function promote_to_nD(M::AbstractArray{T, N}, n::Integer, fill_elem::T) where {T, N}
    ndims(M) == n && return M
    n < ndims(M) && throw(error("Cannot reduce the number of dimensions this array has.  See `resize`."))

    for d in (ndims(M) + 1):n
        M = append_n_times(M, 1, fill_elem, dims = d)
        M = append_n_times_backwards(M, 1, fill_elem, dims = d)
    end

    return M
end
promote_to_nD(M::AbstractArray{T, N}, n::Integer) where {T <: Number, N} = promote_to_nD(M, n, zero(T))

"""
    promote_to_3D(M::AbstractArray{T, N}, fill_elem::T)
    promote_to_3D(M::AbstractArray{T, N}) where T <: Number

A simple, self-evident wrapper for `promote_to_nD`.   No `fill_elem` assumes that the value is zero (if it can be).
"""
promote_to_3D(M::AbstractArray{T, N}, fill_elem::T) where {T, N} = promote_to_nD(M, 3, fill_elem)
promote_to_3D(M::AbstractArray{T, N}) where {T <: Number, N} = promote_to_3D(M, zero(T))

"""
    reduce_dim(A::AbstractArray{T, N}; reduce_by::Int = 1)

Reduce an `N`-dimensional array `A` to dimensions `N - reduce_by`.  Can reduce arbitrarily.  Currently unsafe (i.e., doesn't check that `reduce_by < ndims(A)`.

```julia
julia> A = zeros(Int, 3, 3, 3);

julia> reduce_dims(A)
3×3 Matrix{Vector{Int64}}:
 [0, 0, 0]  [0, 0, 0]  [0, 0, 0]
 [0, 0, 0]  [0, 0, 0]  [0, 0, 0]
 [0, 0, 0]  [0, 0, 0]  [0, 0, 0]
```
"""
function reduce_dims(A::AbstractArray{T, N}; reduce_by::Int = 1) where {T, N}
    sz = size(A)
    out_A = Array{Array{T, reduce_by}, N - reduce_by}(undef, sz[1:(end - reduce_by)])
    for i in CartesianIndices(sz[1:(end - reduce_by)])
        out_A[i] = A[i, (Colon() for _ in reduce_by)...]
    end
    return out_A
end


text:
import MinimalRLCore

using RecipesBase
import Images
using Images: permutedims, permutedims!

import ArcadeLearningEnvironment
const ALE = ArcadeLearningEnvironment

"""
    Atari(gamename; seed, frameskip, color_averaging, repeat_action_probabiliyt, gray_scale)

An interface adapted from https://github.com/JuliaML/AtariAlgos.jl/blob/master/src/AtariAlgos.jl with a backend
implemented by https://github.com/JuliaReinforcementLearning/ArcadeLearningEnvironment.jl. Because we want to have
some better fidelity with settings, reimplementing is easier than writting a wrapper around a wrapper around a wrapper...

Many current implementations of these utilities in Python use a compositional design pattern. I find this often
obfuscates what exactly is going on and forces users to travers several classes and types to understand the code.
Here we've decided to take a more procedural approach. The end user can still wrap this type to provide their own
settings or use this type as a template for their own implementation.

Use cases found in the literature not yet implemented:
- Terminate on life lost -> will implement w/ a boolean setting.
-

## Arguments
- `gamename::AbstractString`: Name of game (should be found in rom list.)
- `seed::Int32`: The seed passed to the ALE environment. This needs to be a 32 bit integer for ALE.
- `frameskip`: The number of frames skipped between actions. If set to 1, `color_averaging` must be set to `:none`
- `color_averaging`: Can be any of the set `(:none, :max, :average)`
- `gray_scale`: If true this will return the gray scale images provided by the ALE library. If false the full color images will be returned.
"""
mutable struct Atari <: MinimalRLCore.AbstractEnvironment
    # Pointer to ALE
    ale::ALE.ALEPtr

    # Details of current episode
    lives::Int
    died::Bool
    reward::Float64
    score::Float64

    # settings
    gamename::String
    frameskip::Int
    color_averaging::Symbol
    gray_scale::Bool

    # details needed for constructing the base state.
    width::Int
    height::Int
    rawscreen::Vector{Cuchar}  # raw screen data from the most recent frame
    state_buffer::Array{Array{UInt8}, 1}
end

function Atari(gamename::AbstractString,
               seed,
               frameskip,
               color_averaging,
               repeat_action_probability,
               gray_scale,
               terminal_on_life=false)

    @assert gamename ∈ ALE.getROMList()
    @assert color_averaging ∈ (:none, :max, :average)
    @assert frameskip >= 1
    @assert !(frameskip == 1 && color_averaging == :none)

    ale = ALE.ALE_new()
    ALE.setInt(ale, "random_seed", UInt16(seed))
    ALE.setFloat(ale, "repeat_action_probability", Float32(repeat_action_probability))
    ALE.loadROM(ale, gamename)

    w = ALE.getScreenWidth(ale)
    h = ALE.getScreenHeight(ale)

    rawscreen = if gray_scale
        Array{Cuchar}(undef, w * h)
    else
        Array{Cuchar}(undef, w * h * 3)
    end

    state = if gray_scale
        [fill(zero(UInt8), h, w), fill(zero(UInt8), h, w)]
    else
        [fill(zero(UInt8), h, w, 3), fill(zero(UInt8), h, w, 3)]
    end

    Atari(ale, 0, false, 0.0, 0.0, gamename, frameskip, color_averaging, gray_scale, w, h, rawscreen, state)
end

Base.show(io::IO, env::Atari) =
    println(io, "$(env.gamename)(frameskip=$(env.frameskip), color_averaging=$(env.color_averaging), gray_scale=$(env.gray_scale))")

# Various default versions of Atari found in literature
"""
    RevisitingALEAtari(gamename, seed)

This is the settings used in "Revisiting the Arcade Learning Environment: Evaluation Protocols and Open Problems for General Agents." by Marlos Machado et. al. See https://arxiv.org/abs/1709.06009.

"""
RevisitingALEAtari(gamename, seed) =
    Atari(gamename, seed, 5, :max, 0.25f0, true)

"""
    DopamineALEAtari(gamename, seed)

This is the settings used in the Dopamine python package:
"""
DopamineALEAtari(gamename, seed) =
    Atari(gamename, seed, 4, :max, 0.25f0, true)

# Some utility functions
function Base.close(env::Atari)
    ALE.ALE_del(env.ale)
end

MinimalRLCore.get_actions(env::Atari) = ALE.getLegalActionSet(env.ale)
valid_action(env::Atari, action) = action in MinimalRLCore.get_actions(env)
get_minimal_actions(env::Atari) = ALE.getMinimalActionSet(env.ale)

function update_state!(env::Atari, t)
    # get the raw screen data

    if env.gray_scale
        ALE.getScreenGrayscale!(env.ale, env.rawscreen)
        permutedims!(env.state_buffer[t], reshape(env.rawscreen, (env.width, env.height)), (2,1))
    else
        ALE.getScreenRGB!(env.ale, env.rawscreen)
        permutedims!(env.state_buffer[t], reshape(env.rawscreen, (3, env.width, env.height)), (3,2,1))
    end
    env.lives = ALE.lives(env.ale)
    return env.state_buffer[t]
end

function _pool_state!(env::Atari)
    if env.color_averaging == :max
        env.state_buffer[2] .= max.(env.state_buffer[1], env.state_buffer[2])
    elseif env.color_averaging == :average
        env.state_buffer[2] .= UInt8.(round.(mean([env.state_buffer[1], env.state_buffer[2]])))
    end
end

# Set seed default to 0
function MinimalRLCore.start!(env::Atari)
    ALE.reset_game(env.ale)
    env.lives = 0
    env.died = false
    env.reward = 0.0f0
    env.score = 0
    fill!(env.state_buffer[1], 0)
    fill!(env.state_buffer[2], 0)
    update_state!(env, 2)
end

MinimalRLCore.start!(env::Atari, rng::AbstractRNG) = MinimalRLCore.start!(env)

function MinimalRLCore.environment_step!(env::Atari, action)

    env.reward = 0.0f0
    for i ∈ 1:env.frameskip
        env.reward += ALE.act(env.ale, action)
        if is_terminal(env)
            break
        elseif (env.frameskip - i) < 2
            update_state!(env, env.frameskip - i + 1)
        end
    end

    _pool_state!(env)
    env.score += env.reward
end

MinimalRLCore.get_reward(env::Atari) = env.reward
MinimalRLCore.is_terminal(env::Atari) = ALE.game_over(env.ale)
MinimalRLCore.get_state(env::Atari) = env.state_buffer[2]

"""
    image_manip_atari(s)

Helper functions for doing image pre-processing. This does the standard thing of turing the image gray scale
(if necessary) and resizes the image too 84*84.

As a note, Images.Gray does the conversion to Luminance as done in the original DQN paper (`0.299 * r + 0.587 * g + 0.114 * b`)
but then rounds to fit in a UInt8. It is unclear if this effects the final performance of the DQN or not... (which is a bit dissapointing).

"""
function image_manip_atari(s::Array{UInt8, 3})
    Images.colorview(Images.RGB{Images.N0f8}, permutedims(s, (3, 1, 2))) |> # Colorview of array
        (img)->Images.Gray{Images.N0f8}.(Images.imresize(img, (84,84))) |> # Resize image to 84x84
        (Images.rawview ∘ Images.channelview)
end

function image_manip_atari(s::Array{UInt8, 2})
    Images.imresize(Images.Gray{Images.N0f8}.(s./UInt8(255)), (84,84)) |>
        (Images.rawview ∘ Images.channelview)
end

image_norm(img) = img./255f0

#####
# Visualization through RecipeBase.jl
#####

function get_colorview(env::Atari, gray_scale::Val{true})
    Images.colorview(Images.Gray, permutedims(reshape(env.rawscreen, (env.width, env.height)), (2,1))./255)
end

get_colorview(env::Atari, gray_scale::Val{false}) =
    Images.colorview(Images.RGB, permutedims(reshape(env.rawscreen, (3, env.width, env.height)), (1,3,2))./255)

get_colorview_image_manip(env::Atari, gray_scale::Val{true}) =
    Images.Gray.(image_norm(image_manip_atari(MinimalRLCore.get_state(env))))

get_colorview_image_manip(env::Atari, gray_scale::Val{false}) =
    Images.Gray.(image_norm(Images.colorview(Images.Gray, image_manip_atari(MinimalRLCore.get_state(env)))))

@recipe function f(env::Atari, resize=false)
    ticks := nothing
    foreground_color_border := nothing
    grid := false
    legend := false
    aspect_ratio := 1
    if resize
        get_colorview_image_manip(env, Val(env.gray_scale))
    else
        get_colorview(env, Val(env.gray_scale))
    end
end


text:
using Geotherm.Geometry: Point2D
using Geotherm.Integrate: runge_kutta_iter, runge_kutta

@testset "Test `runge_kutta`" begin
    @test runge_kutta_iter(
        Point2D(1.0, 1.0),
        (x, y) -> x * 2 + y * 3,
    ) == Point2D(1.01, 1.05085856375)
    @test runge_kutta(
        Point2D(1.0, 1.0),
        (x, y) -> x * 2 + y * 3,
        0.01,
        4,
    ) == [
        Point2D(1.0, 1.0),
        Point2D(1.01, 1.05085856375),
        Point2D(1.02, 1.103469031571201),
        Point2D(1.03, 1.157884756885266),
    ]
end


text:
clearconsole()

using Dates
using DataFrames
using Documenter#, DocumenterTools
using Distributions, Statistics, StatsBase
using Gadfly
using Documenter,Test, DocumenterTools
using MCHammer, Compose, Cairo, Fontconfig

function DoIt()

    makedocs(
    sitename="MCHammer.jl",
    source  = "src",
    build   = "local",
    modules =[MCHammer],
    pages = Any[
    "Home" => "index.md",
    "User Manual" => Any[
        "Simulation Functions" => "manual/1_functions.md",
        "Charting & Analyzing" => "manual/2_charts.md",
        "Time-Series Functions" => "manual/3_time_series.md",
    ],
    "Tutorials" => Any[
        "My First Model" => "tutorials/1_first_model.md",
        "Correlating Inputs" => "tutorials/2_Correlated_Model.md",
        "Simulated CashFlow Model" => "tutorials/3_NPV_testmodel.md",
    ],
    "Import / Export Data" => "manual/4_moving_results.md"
    ])
end

# @testset "MCHammer" begin
#
#     doctest(MCHammer; manual = true)
#
# end


text:
module TestVanDerPol
include("preamble.jl")

using Bifurcations: LimitCycleProblem
using Bifurcations.Examples.DuffingVanDerPol
using Bifurcations.Continuations: as, ContinuationSolution, sweeps_as_vectors

using DiffEqBase: remake
param = DuffingVanDerPol.DuffingVanDerPolParam(
    d = 0.01,
)
ode = remake(
    DuffingVanDerPol.ode,
    p = param,
)

num_mesh = 20
degree = 3

# Start continuation of the limit cycle at the small damping limit.
xs0 = let
    n = num_mesh * degree
    dt = 2π / n
    ts = ((1:n) .- 1) .* dt
    hcat(
        2 .* cos.(ts),
        -2 .* sin.(ts),
    )'
end

prob = LimitCycleProblem(
    de_prob = ode;
    param_axis = DuffingVanDerPol.param_axis,
    num_mesh = num_mesh,
    degree = degree,
    xs0 = xs0,
    l0 = 2π,
    t_domain = (0.0001, 1.5), # bound it so that it works with above `num_mesh`
    t0 = get(ode.p, DuffingVanDerPol.param_axis),
)

@test size(prob.xs0) == (2, num_mesh * degree)

solver = init(
    prob;
    start_from_nearest_root = true,
    bidirectional_first_sweep = false,
)
solve!(solver)

#=
parameters, = let sol = as(solver.sol, ContinuationSolution)
    sweeps_as_vectors(sol, length(sol.sweeps[1].u[1]))
end
=#
parameters = let sol = as(solver.sol, ContinuationSolution)
    [u[end] for u in sol.sweeps[1].u]
end
# @test minimum(parameters) <= prob.t_domain[1]
@test maximum(parameters) >= prob.t_domain[end]

end  # module


text:
"""
	cracks2moduli(E₀, nu₀, α1, α3, β11, β13, β33)

Compute 6×6 stiffness matrix (in Voight notation) for a cracked material characterised by normalised crack density parameters α1, α3, β11, β13 and β33, with solid matrix moduli E₀ and nu₀. Crack density parameters are normalised by h = 3E₀(2-nu₀)/(32(1-nu₀²)). From Sayers and Kachanov 1995.
"""
function cracks2moduli(E₀, nu₀, α1, α3, β11, β13, β33)
    S011 = 1/E₀
    S012 = -nu₀/E₀
    h = cracknormfactor(E₀, nu₀)

    D = (S011 + α3/h + β33/h)*(S011 + S012 + α1/h + 4β11/h/3) -
        2*(S012 + β13/h)^2

    c11 = (1/2)*((S011+α3/h+β33/h)/D + 1/(S011-S012+α1/h+2β11/h/3))
    c33 = (S011 + S012 + α1/h + 4β11/h/3)/D
    c44 = 1/(2S011 - 2S012 + α1/h + α3/h + 4β13/h)
    c13 = -(S012 + β13/h)/D
    c66 = 1/(S011-S012+α1/h+2β11/h/3)/2

    return stiffmatrixTI(c11, c13, c33, c44, c66)
end

"""
	cracknormfactor(E₀, ν₀)

Compute normalisation factor h=3E₀*(2-nu₀)/(32(1-nu₀²)).
"""
function cracknormfactor(E₀, ν₀)
    return 3E₀*(2-ν₀)/(32(1-ν₀^2))
end

"""
	findcracks(p::Vector{VMeasure},
                    sv::Vector{VMeasure},
                    sh::Vector{VMeasure},
                    ρ::Real,
                    E₀::Real,
                    nu₀::Real,
                    α₁₀::Tuple{Real,Real},
                    α₃₀::Tuple{Real,Real},
                    β₁₁₀::Tuple{Real,Real},
                    β₁₃₀::Tuple{Real,Real},
                    β₃₃₀::Tuple{Real,Real})

Estimate normalised crack density parameters and phase angles from measurements of group P, SV ans SH wave velocities, using quasinewton method. A priori values ()₀ are tuples with value and standard deviation (a priori gaussian).

Return α₁₁, α₃₃, β₁₁₁₁, β₁₁₃₃, β₃₃₃₃, θp, θsv, θsh, CMpost.
"""
function findcracks(p::Vector{VMeasure},
                    sv::Vector{VMeasure},
                    sh::Vector{VMeasure},
                    ρ::Real,
                    E₀::Real,
                    nu₀::Real,
                    α₁₀::Tuple{Real,Real},
                    α₃₀::Tuple{Real,Real},
                    β₁₁₀::Tuple{Real,Real},
                    β₁₃₀::Tuple{Real,Real},
                    β₃₃₀::Tuple{Real,Real})

    dobs, CDi, np, nsv, nsh = inputdobs(p,sv,sh)

    #make vector of model parameters (first guess)
    m₀ = [α₁₀[1],α₃₀[1],β₁₁₀[1],β₁₃₀[1],β₃₃₀[1]]
    σm = [α₁₀[2],α₃₀[2],β₁₁₀[2],β₁₃₀[2],β₃₃₀[2]]

    C₀ = cracks2moduli(E₀, nu₀, m₀...)
    ϵ₀, δ₀, γ₀ = moduli2thomsen(C₀)
    vp,vs = velocities0(C₀,ρ) # here density does not matter

    for x in p
        push!(m₀, phaseangleP(x.angle, vp/vs, ϵ₀, δ₀))
        push!(σm, 100)
    end
    for x in sv
        push!(m₀, phaseangleSV(x.angle, vp/vs, ϵ₀, δ₀))
        push!(σm, 100)
    end
    for x in sh
        push!(m₀, phaseangleSH(x.angle, γ₀))
        push!(σm, 100)
    end

    CM = Diagonal(σm.^2)
    CMi = inv(CM)

    m, CMpost = quasinewton(dobs, _gf, m₀, CMi, CDi, (E₀, nu₀, ρ,np,nsv,nsh))

    return m[1:5], m[6:5+np], m[6+np:5+np+nsv], m[6+np+nsv:5+np+nsv+nsh], CMpost

end

function _gf(m,E₀,nu₀,ρ,np,nsv,nsh)
    C₀ = cracks2moduli(E₀, nu₀, m[1:5]...)
    ϵ, δ, γ = moduli2thomsen(C₀)
    vp0,vs0 = velocities0(C₀,ρ)
    mt = vcat([vp0,vs0,ϵ,δ,γ], m[6:end])
    return _gt(mt, np, nsv, nsh)
end


text:
module PandorasBoxLaserValve

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/laserValve" LaserValve(x::Integer, y::Integer, direction::String="Right", delay::Int=4)

const directionRotations = Dict{String, Number}(
    "Right" => 0,
    "Down" => pi / 2,
    "Left" => pi,
    "Up" => pi / 2 * 3
)

const placements = Ahorn.PlacementDict(
    "Laser Valve ($direction) (Pandora's Box)" => Ahorn.EntityPlacement(
        LaserValve,
        "point",
        Dict{String, Any}(
            "direction" => direction
        )
    ) for (direction, rot) in directionRotations
)

const offsets = Dict{String, Tuple{Int, Int}}(
    "Right" => (0, 0),
    "Down" => (16, 0),
    "Left" => (16, 16),
    "Up" => (0, 16)
)

Ahorn.editingOptions(entity::LaserValve) = Dict{String, Any}(
    "direction" => collect(keys(offsets))
)

sprite = "objects/pandorasBox/laser/valve/valve2"

function Ahorn.selection(entity::LaserValve)
    x, y = Ahorn.position(entity)

    direction = get(entity, "direction", "Right")
    rotation = get(directionRotations, direction, 0)

    return Ahorn.getSpriteRectangle(sprite, x, y, jx=0.5, jy=0.5)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::LaserValve, room::Maple.Room)
    direction = get(entity, "direction", "Right")
    rotation = get(directionRotations, direction, 0)

    Ahorn.drawSprite(ctx, sprite, offsets[direction]..., jx=0.5, jy=0.5, rot=rotation)
end

end

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
### Use the lasso path for the crime train example.

include("lasso.jl")
include("utils.jl")

using DataFrames
using DelimitedFiles: readdlm
using Plots
import UnicodePlots as Uplot

function DataSummary()
    P, H = readdlm("./L1_Lasso_Module/crime-train.txt", '\t'; header=true)
    Df = DataFrame(P, reshape(H, length(H)))
    return Df
end

PrintTitle("A Summarization of the Data: The Predictors")
TheData= DataSummary()
println(describe(TheData))
println("The size of the dataframe is: $(size(TheData))")
println("The predictant is:\"$(names(TheData)[1])\"")
PrintTitle("We are gonna prepare the into a matrix and vector. ")

function PrepareTheDataMatrix(data_frame::DataFrame)
    """
        Min, max, median and mean are assumed to be important predictor
    """
    FeatureMatrix = Matrix(data_frame[!, 2:end])
    Predictants = data_frame[!, 1]
    return FeatureMatrix, Predictants
end

A, b = PrepareTheDataMatrix(TheData)

function AnalyaisWithLasso(A::Matrix, b::Matrix)
    """
        Do the lasso path thing on the data matrix,
        features are columns of the matrix, instances are rows.

        Returns:
            tuple where the first element is the data matrix and the
            second element is the column vector of all the regularization
            parameter λ.
    """
    Instance = LassoSCOP(A, b)
    Results, λs =  LassoPath(Instance)
    VisualizeLassoPath(Instance, "cirme-data-lassopath-$(rand(1:999999999)).png")

end

AnalyaisWithLasso(A, reshape(b, (length(b), 1)))


text:
export prox_grad!, ProxGradParams

# solves minimize(f(x) + g(x))
# where f is differentiable and the prox of g is easy to evaluate
# prox_g(x, alpha) solves minimize(g(z) + 1/(2*alpha)||x-z||^2
function prox_grad!(x, # the initial value of the variable
	f::Function,
	grad_f::Function,
	g::Function,
	prox_g::Function,
	params::OptParams = SimpleParams(ConstantStepSize(1), 100);
	verbose::Bool = false)

	# initialize
	objective = x -> f(x) + g(x)
	objval = objective(x)
	take_step(x, stepsize) = prox_g(x - stepsize*grad_f(x), stepsize)

	if verbose
		@printf("%10s%12s\n", "iter", "obj")
		@printf("%10d%12.4e\n", 0, objval)
	end

	for iter = 1:params.maxiters
		oldobjval = copy(objval)
		x, objval = step(params.stepsizerule,
			             objective,
			             take_step,
			             x,
			             objval0 = objval)

		verbose && iter%1==0 && @printf("%10d%12.4e\n", iter, objval)

		if stop(params, objval, oldobjval)
			verbose && println("stopping criterion satisfied")
			break
		end
	end
	return x
end

type ProxGradParams<:OptParams
	stepsizerule::StepSizeRule
	maxiters::Int
end
ProxGradParams() = ProxGradParams(ConstantStepSize(.5), 100)

function step(s::HopefulStepSize, objective::Function, take_step::Function, x0;
	objval0 = objective(x0))
	stepsize = s.initial_stepsize
	x = take_step(x0, stepsize)
	objval = objective(x)
	while objval >= objval0
		@show stepsize *= s.decrease_by
		x = take_step(x0, stepsize)
		objval = objective(x)
	end
	# hope and change!
	if stepsize < s.initial_stepsize
		s.initial_stepsize = stepsize
	else
		s.initial_stepsize *= s.increase_by
	end
	return x, objval
end

max_stars_repo_name: madeleineudell/FirstOrderOptimization


text:
import ..LambdaCalculus: alpha_equivalent, freevars, substitute

const EqList = Vector{Tuple{Symbol, Symbol}}

function alpha_equivalent(t1::Var, t2::Var, equivalences::EqList)
    # see: https://cs.stackexchange.com/a/76621/28930
    # Either both bindings must appear in the same spot (they were bound by the same Abs) or
    # neither appears and they are equal (they are both free and they are equal).

    for (x, y) in equivalences
        if x == t1.name && y == t2.name
            return true
        elseif x != t1.name && y != t2.name
            continue
        else
            return false
        end
    end

    return t1.name == t2.name
end

alpha_equivalent(t1::App, t2::App, equivalences::EqList) =
    alpha_equivalent(t1.car, t2.car, equivalences) && alpha_equivalent(t1.cdr, t2.cdr, equivalences)
alpha_equivalent(t1::Abs, t2::Abs, equivalences::EqList) =
    alpha_equivalent(t1.body, t2.body, [(t1.boundname, t2.boundname); equivalences])
alpha_equivalent(t1::Term, t2::Term, equivalences::EqList) = false
alpha_equivalent(t1::Term, t2::Term) = alpha_equivalent(t1, t2, Tuple{Symbol, Symbol}[])

freevars(t::Var) = Set([t.name])
freevars(t::Abs) = filter(!isequal(t.boundname), freevars(t.body))
freevars(t::App) = freevars(t.car) ∪ freevars(t.cdr)

function substitute(name::Symbol, s::Term, t::Var)
    return name == t.name ? s : t
end

function substitute(name::Symbol, s::Term, t::App)
    return App(substitute(name, s, t.car), substitute(name, s, t.cdr))
end

function substitute(name::Symbol, s::Term, t::Abs)
    if name == t.boundname
        return t
    else
        fv = freevars(t.body)

        if name ∉ fv
            return Abs(t.boundname, substitute(name, s, t.body))
        else
            freshvar = Var(freshname(addprime(name), fv))
            t′ = Abs(freshvar.name, substitute(name, freshvar, t.body))
            return substitute(name, s, t′)
        end
    end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using Test

Test.@testset "optimization" begin
    include("test_mode_estimators.jl")
end


text:
include("vector.jl")
include("types.jl")
include("properties.jl")
import Images
using ImageView, FileIO, FixedPointNumbers,
			 Colors, ColorTypes, ImageMetadata, Base

			 println("Print some message!")
			 print("Print some message!")
someBoolValue = true
# Constants
AMBIENT = get(rayTracerProperties, "ambient", 0.1)

# Method to calculate the Intersection of a ray on a sphere
function intersection(s::Sphere, l::Ray)
	v = dot(l.direction, (l.origin - s.center))
	disc = v ^ 2 - ( (dot(l.origin - s.center, l.origin - s.center) - s.radius ^ 2) * dot(l.direction, l.direction))
	if disc < 0
		return Intersection(Vector(), -1, Vector(), s)
	else
		d = sqrt( disc)
		d1 = (-v + d)
		d2 = (-v - d)
		if 0 < d1 < d2
			point = l.origin + (l.direction * d1)
			distance = d1
		elseif 0 < d2 < d1
			point = l.origin + (l.direction * d2)
			distance = d2
		end
		return Intersection(point, distance, normal(point - s.center), s)
	end
end

function findIntersection(r::Ray, objects, ignore=nothing)
	intersect = Intersection(Vector(), -1, Vector(), nothing)
	for obj in objects
		if obj != ignore
			newIntersect = intersection(obj, r)
			if intersect.distance < 0 && newIntersect.distance > 0
				intersect = newIntersect
			elseif 0 < newIntersect.distance < intersect.distance
				intersect = newIntersect
			end
		end
	end
	return intersect
end

function blinnPhongShadeColor(l::Ray, intersect, lightSource, camera)
	n = intersect.normal
	l = normal(intersect.point - lightSource)
	v = normal(intersect.point - camera)
	h = (v + l) * (1 / magnitude(v + l))
	kSpecular = get(rayTracerProperties, "specular", Vector())
	kDiffuse = intersect.object.color
	I = 1 - AMBIENT

	return (kDiffuse * I * max(0, dot(n, l))) + (kSpecular * I * max(0, dot(n, h)))
end

function lambertShadeColor(l::Ray, intersect, lightSource)
	shade = dot(intersect.normal, normal(intersect.point - lightSource))
	if(shade <= 0)
		shade = 0
		return Vector(255, 255, 255) * AMBIENT
	end
	return intersect.object.color * (1 - AMBIENT) * shade
end

function traceRay(l::Ray, lightSource, objects, cameraPos)
	intersect = findIntersection(l, objects)
	if intersect.object == nothing
		return Vector(255, 255, 255) * AMBIENT
	end
	return blinnPhongShadeColor(l, intersect, lightSource, cameraPos)
end

function traceWorker(i, j, imageArray, cameraPos, lightSource, objects)
	ray = Ray(cameraPos, normal(Vector((i-0.5)/45, (j-0.5)/45, 0) - cameraPos))
	col = traceRay(ray, lightSource, objects, cameraPos)
	imageArray[j,i] = RGB(col.x/256, col.y/256, col.z/256)
end

function render()

	# Read from properties
	rtp = rayTracerProperties
	lightSource = get(rtp, "lightSource", Vector())
	cameraPos = get(rtp, "cameraPos", Vector())
	objects = get(rtp, "objects", nothing)
	imageWidth = get(rtp, "imageWidth", 0)
	imageHeight = get(rtp, "imageHeight", 0)
	pixelWidth = get(rtp, "pixelWidth", 0)
	#imageProperties = get(rtp, "imageProperties", nothing)

	# First set up image array
	image = fill(RGB(0,0.5,0), imageWidth, imageHeight)

	for x in 1:imageWidth
		for y in 1:imageHeight
			traceWorker(x, y, image, cameraPos, lightSource, objects)
		end
		end
	#image.data = imageArray
	return image

end
finalImage = render()
imshow(finalImage, pixelspacing = [1,1])
#display(finalImage)


text:
import DSP.Periodograms:
    PSDOnly,
    nextfastfft,
    Spectrogram,
    spectrogram,
    stft,
    fft2pow!,
    fft2oneortwosided!,
    compute_window,
    arraysplit,
    stfttype

function DSP.spectrogram(
    s::AbstractVector{T},
    estimator::Function,
    n::Int = length(s) >> 3,
    noverlap::Int = n >> 1;
    onesided::Bool = eltype(s) <: Real,
    nfft::Int = nextfastfft(n),
    fs::Real = 1,
    window::Union{Function,AbstractVector,Nothing} = nothing,
) where {T}

    out = stft(
        s,
        estimator,
        n,
        noverlap,
        PSDOnly();
        onesided = onesided,
        nfft = nfft,
        fs = fs,
        window = window,
    )
    Spectrogram(
        out,
        onesided ? DSP.rfftfreq(nfft, fs) : DSP.fftfreq(nfft, fs),
        (n/2:n-noverlap:(size(out, 2)-1)*(n-noverlap)+n/2) / fs,
    )

end

function stft(
    s::AbstractVector{T},
    estimator::Function,
    n::Int = length(s) >> 3,
    noverlap::Int = n >> 1,
    psdonly::Union{Nothing,PSDOnly} = nothing;
    onesided::Bool = eltype(s) <: Real,
    nfft::Int = nextfastfft(n),
    fs::Real = 1,
    window::Union{Function,AbstractVector,Nothing} = nothing,
) where {T}

    win, norm2 = compute_window(window, n)
    sig_split = arraysplit(s, n, noverlap, nfft, win)
    nout = onesided ? (nfft >> 1) + 1 : nfft
    out = zeros(stfttype(T, psdonly), nout, length(sig_split))

    freqs = onesided ? DSP.rfftfreq(nfft, fs) : DSP.fftfreq(nfft, fs)
    r = fs * norm2 / sqrt(length(freqs))
    offset = 0
    for (i, sig) in enumerate(sig_split)
        # mul!(tmp, plan, sig)
        tmp = estimator(sig, freqs)
        if isa(psdonly, PSDOnly)
            fft2pow!(out, tmp, nfft, r, onesided, offset)
        else
            fft2oneortwosided!(out, tmp, nfft, onesided, offset)
        end
        offset += nout
    end
    out
end

"""
    model_spectrum(f, h, args...; kwargs...)

DOCSTRING

# Arguments:
- `f`: the model-estimation function, e.g., `ar,arma`
- `h`: The sample time
- `args`: arguments to `f`
- `kwargs`: keyword arguments to `f`

# Example:
```
using ControlSystemIdentification, DSP
T = 1000
s = sin.((1:T) .* 2pi/10)
S1 = spectrogram(s,window=hanning)
estimator = model_spectrum(ar,1,2)
S2 = spectrogram(s,estimator,window=rect)
plot(plot(S1),plot(S2)) # Requires the package LPVSpectral.jl
```
"""
function model_spectrum(f, h, args...; kwargs...)
    function (s::AbstractArray{T}, freqs) where {T}
        d = iddata(s, h)
        model = f(d, args...; kwargs...)
        tmp = vec(Complex{T}.(freqresp(model, T(2pi) .* freqs)))
    end
end

# function rootspectrogram(model, fs)
#     roots = map(1:length(model)) do i
#         sys = tf(1,[1;-reverse(model.θ[:,i])],1)
#         (sort(pole(sys), by=imag, rev=true)[1:end÷2])
#         # log.(sort(complex.(eigvals(model.At[:,:,i])), by=imag, rev=true)[1:end÷2])
#     end
#     S = reduce(hcat,roots)
#     fs/(2pi) .* angle.(S)'
# end


text:
using BSONify, Dates
import Mongoc

# Define a custom parametric struct.
struct MyStruct{T, S}
    myDateTime::DateTime
    myArray::Array{S}
    myParametric::T
    myDict::Dict{T, S}
    myNamedTuple::NamedTuple{(:a, :b), Tuple{Int16, Float32}}
    myChild::Union{MyStruct{T, S}, S, Nothing} # Support for recursive structures and type unions.
end

# Define a custom type from the custom struct.
MyType = MyStruct{UInt8, String}

# Create some data.
myData = MyType(
    now(),
    ["The", "answer", "is"],
    42,
    Dict(4 => "is", 2 => "answer?!"),
    (a=1, b=2.3),
    MyType(
        now(),
        split("It's not a question, but an answer learned in time"),
        0,
        Dict(),
        (a=1, b=2),
        "sweet child o' mine"
    )
)

println("myData:\n", myData)

# Convert our data to BSON.
bson = Mongoc.BSON(myData)

println("\nbson:\n", bson)

# Convert BSON back to our custom type.
myReconstitutedData = as_type(MyType, bson)

println("\nmyReconstitutedData:\n", myReconstitutedData)

max_stars_repo_name: OliverColeman/BSONify

max_forks_repo_name: OliverColeman/BSONify


text:
abstract type AbstractSequentialLinearLeastSquares end

mutable struct InformationFormSequentialLinearLeastSquares{BFVT} <: AbstractSequentialLinearLeastSquares
    # Basis functions linear with respect to quantities to be estimated
    # h[i] ∀ i ∈ [0, n] should be a function which takes the independant
    # variable as an argument (i.e h[i](tₖ) where tₖ is the time at
    # the kth epoch.
    h::BFVT

    # Basis function matrix
    H::Matrix{Float64}

    # Estimate vector
    xhat::Vector{Float64}

    # Newest measurements
    y::Vector{Float64}

    # Newest measurement weights
    W::Symmetric{Float64, Matrix{Float64}}

    # Newest independant variables
    t::Vector{Float64}

    # Allocated information matricies
    Λkp1::Matrix{Float64}
    Λk::Matrix{Float64}
end

mutable struct CovarianceFormSequentialLinearLeastSquares{BFVT} <: AbstractSequentialLinearLeastSquares
    # Basis functions linear with respect to quantities to be estimated
    # h[i] ∀ i ∈ [0, n] should be a function which takes the independant
    # variable as an argument (i.e h[i](tₖ) where tₖ is the time at
    # the kth epoch.
    h::BFVT

    # Basis function matrix
    H::Matrix{Float64}

    # Estimate vector
    xhat::Vector{Float64}

    # Newest measurements
    y::Vector{Float64}

    # Newest measurement weights
    W::Symmetric{Float64, Matrix{Float64}}

    # Newest independant variables
    t::Vector{Float64}

    # Allocated covariance matricies
    Pkp1::Matrix{Float64}
    Pk::Matrix{Float64}
end

function SequentialLinearLeastSquares(h::AbstractVector, xhat0::AbstractVector, Λ0::AbstractMatrix, form::InformationForm)
    n = length(h)
    if length(xhat0) != n
        throw(ArgumentError("Initial estimate vector length is not the same length as the basis function vector."))
    end
    r, c = size(Λ0)
    if r != n || c != n
        throw(ArgumentError("Initial information matrix is not the correct size."))
    end

    InformationFormSequentialLinearLeastSquares(h, Matrix{Float64}(undef, (0, 0)),
        xhat0, Vector{Float64}(undef, 0), Symmetric(Matrix{Float64}(undef, (0,0))),
        Vector{Float64}(undef, 0), zeros(n,n), Λ0)
end

function SequentialLinearLeastSquares(h::AbstractVector, xhat0::AbstractVector, form::InformationForm)
    n = length(h)
    SequentialLinearLeastSquares(h, xhat0, zeros(n,n), form)
end

function SequentialLinearLeastSquares(h::AbstractVector, form::InformationForm)
    n = length(h)
    SequentialLinearLeastSquares(h, zeros(n), form)
end

function SequentialLinearLeastSquares(h::AbstractVector, xhat0::AbstractVector, P0::AbstractMatrix, form::CovarianceForm)
    n = length(h)
    if length(xhat0) != n
        throw(ArgumentError("Initial estimate vector length is not the same length as the basis function vector."))
    end
    r, c = size(P0)
    if r != n || c != n
        throw(ArgumentError("Initial covariance matrix is not the correct size."))
    end

    CovarianceFormSequentialLinearLeastSquares(h, Matrix{Float64}(undef, (0,0)),
        xhat0, Vector{Float64}(undef, 0), Symmetric(Matrix{Float64}(undef, (0,0))),
        Vector{Float64}(undef, 0), zeros(n,n), P0)
end

function FeedMeasurementBatch!(slls::AbstractSequentialLinearLeastSquares,
                               ys::AbstractVector, W::AbstractMatrix, ts::AbstractVector)
    m = length(ys)
    n = length(slls.h)
    if m != length(ts)
        throw(ArgumentError("Measurement vector and indepandant variable vector are not the same length."))
    end
    r, c = size(W)
    if r != c
        throw(ArgumentError("Weight matrix must be square!"))
    end

    # Set measurements and indepandant variable vectors
    slls.y = ys
    slls.W = Symmetric(W)
    slls.t = ts

    # Size and fill H matrix
    H = zeros(m, n)
    @inbounds for row in 1:m
        for col in 1:n
            H[row, col] = slls.h[col](slls.t[row])
        end
    end
    slls.H = H

    return nothing
end

function ComputeEstimate!(islls::InformationFormSequentialLinearLeastSquares)
    # Compute information update
    islls.Λkp1 .= islls.Λk .+ transpose(islls.H)*islls.W*islls.H

    # Compute gain matrix
    K    = transpose(islls.H)*islls.W
    Pkp1 = factorize(islls.Λkp1)
    ldiv!(Pkp1, K)

    # Compute estimate
    islls.xhat .+= K*(islls.y .- islls.H*islls.xhat)

    # Increment information matricies
    islls.Λk .= islls.Λkp1

    return nothing
end

function ComputeEstimate!(cslls::CovarianceFormSequentialLinearLeastSquares)
    # Compute gain
    fm = factorize(cslls.H*cslls.Pk*transpose(cslls.H) .+ inv(cslls.W))
    K  = cslls.Pk*transpose(cslls.H)
    rdiv!(K, fm)

    # Update estimate
    cslls.xhat .+= K*(cslls.y .- cslls.H*cslls.xhat)

    # Update covariance
    tm = -K*cslls.H
    map(i -> tm[i,i] += 1.0, 1:length(cslls.xhat))
    cslls.Pkp1 .= tm*cslls.Pk
    cslls.Pk   .= cslls.Pkp1

    return nothing
end

function GetEstimate(slls::AbstractSequentialLinearLeastSquares)
    return slls.xhat
end

function GetMeasurementEstimate(slls::AbstractSequentialLinearLeastSquares, t::Real)
    y = 0.0
    @inbounds for i in 1:length(slls.h)
        y += slls.h[i](t)*slls.xhat[i]
    end
    return y
end


text:
# # Simulations with TAMU data and [PowerSimulations.jl](https://github.com/NREL/PowerSimulations.jl)

# **Originally Contributed by**: Clayton Barrows

# ## Introduction

# This is a basic simulation example using the [TAMU Cases](https://electricgrids.engr.tamu.edu/).

# ## Dependencies
using SIIPExamples
pkgpath = dirname(dirname(pathof(SIIPExamples)))
using Dates

# ### Modeling Packages
using PowerSystems
using PowerSimulations

# ### Optimization packages
# For this simple example, we can use the Cbc solver with a relatively relaxed tolerance.
using Cbc #solver
solver = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 1, "ratioGap" => 0.5)

# ### Create a `System` from TAMU data
# We have included some of the TAMU cases (with truncated time series info)
# in the [PowerSystemsTestData](https://github.com/nrel-siip/powersystemstestdata)
# repository for testing, so we can just use that.
PowerSystems.download(PowerSystems.TestData; branch = "master") # *note* add `force=true` to get a fresh copy
base_dir = dirname(dirname(pathof(PowerSystems)));

# The TAMU data format relies on a folder containing `.m` or `.raw` files and `.csv`
# files for the time series data. We have provided a parser for the TAMU data format with
# the `TamuSystem()` fuction.

TAMU_DIR = joinpath(base_dir, "data", "ACTIVSg2000");
sys = TamuSystem(TAMU_DIR)
horizon = 24;
interval = Dates.Hour(24);
transform_single_time_series!(sys, horizon, interval);

# ## Run a PCM
# note that the TAMU data doesn't contain startup and shutdown costs, or minimum up/down
# time limits, so a UC problem merely respects minmum generation levels.

sim_folder = mkpath(joinpath(pkgpath, "TAMU-sim"))
stages_definition = Dict(
    "UC" => Stage(
        GenericOpProblem,
        template_unit_commitment(network = CopperPlatePowerModel),
        sys,
        solver,
    ),
)
order = Dict(1 => "UC")
horizons = Dict("UC" => 24)
intervals = Dict("UC" => (Hour(24), Consecutive()))
DA_sequence = SimulationSequence(
    step_resolution = Hour(24),
    order = order,
    horizons = horizons,
    intervals = intervals,
    ini_cond_chronology = IntraStageChronology(),
)

# ### Define and build a simulation
sim = Simulation(
    name = "TAMU-test",
    steps = 3,
    stages = stages_definition,
    stages_sequence = DA_sequence,
    simulation_folder = sim_folder,
)

build!(sim)

# ### Execute the simulation
#nb execute!(sim)

# ### Load and analyze results
#nb sim_results = SimulationResults(sim);
#nb uc_results = get_stage_results(sim_results, "UC")
#nb read_realized_variables(uc_results, names = [:On__ThermalStandard])[:On__ThermalStandard]


text:
using Pages,Blink

using Compat; import Compat.readstring

Pages.start()

Endpoint("/examples/pages") do request::HTTP.Request
    readstring(joinpath(dirname(@__FILE__),"pages.html"))
end

function example_pages()
    nwin = 3
    for i = 1:nwin
        w = Window(Dict(:title => "Pages: Window #$i", :url => "http://localhost:8000/examples/pages"))
        tools(w)
        wait(Pages.conditions["connected"])
    end

    Pages.broadcast("say","Hello everyone!")
    for (sid,s) in Pages.pages["/examples/pages"].sessions
        Pages.message(s,"say","You are connection #$(s.id).")
    end
end
# example_pages()


text:
# DESCRIPTION
# This file contains the helicopter parameters.. Always re-run DI_design.jl
# after modifying the aircraft parameters.

# ------------------------- MISCELLANEOUS CONSTANTS ----------------------------

# air density [slugs-ft³]
RHO=0.002378
# acceleration of gravity [ft/s²]
g=32.17

# ----------------------------- MASS PROPERTIES --------------------------------

# gross weight [lb]
W=16000.0
# inertias [sl-ft²]
Ix=4659.0*W/16000.0
Iy=38512.0*W/16000.0
Iz=36796.0*W/16000.0
Ixz=1882.0*W/16000.0
# CG location (station, water line, butt line) [in]
FSCG=358.0
BLCG=0.0
WLCG=248.2

# --------------------------------- FUSELAGE -----------------------------------

# fuselage drag [ft²]
CDSF=35.04
# location of fuselage aerodynamic center (station, water line, butt line) [in]
FSfus=345.5
BLfus=0.0
WLfus=234.0
# number of blades
NB=4;
# single blade weight [lb]
WBLADE=256.9
# aircraft gross weight (rotor blades not included) [lb]
WEIGHTNR = W-NB*WBLADE;
# aircraft mass (rotor blades not included) [slug]
mass=WEIGHTNR/g;
# main rotor location (station, water line, butt line) [in]
FSMR=341.25
BLMR=0.0
WLMR=315.0
# fuselage CG location (without rotor) (station, water line, butt line) [in]
FSCGB=((W*FSCG)-(NB*WBLADE*FSMR))/WEIGHTNR
WLCGB=((W*WLCG)-(NB*WBLADE*WLMR))/WEIGHTNR
BLCGB=((W*BLCG)-(NB*WBLADE*BLMR))/WEIGHTNR
# moment arms from fuselage CG to fuselage aero center [ft]
FWT = (FSCG-FSfus)/12
WWT = (WLCG-WLfus)/12
BWT = (BLCG-BLfus)/12
# location of fuselage aero center relative to CG [ft]
rFS=[FWT; BWT; WWT];

# -------------------------------- MAIN ROTOR ----------------------------------

# rotor Radius [ft]
RADIUS=26.83
# rotor speed [rad/s]
OMEGA=27.0
# solidity
SOL=0.0821
# blade lift slope [1/rad]
A0=5.73
# blade twist [rad]
TWIST=-13.0*pi/180.0
# lock number = rho*a0*c*R^4/Ibeta
GAMMA=8.0755
# moment arms from fuselage CG to main rotor hub [ft]
Xh=(FSCGB-FSMR)/12.0
Yh=(BLCGB-BLMR)/12.0
Zh=(WLCGB-WLMR)/12.0
# main rotor location
rMR=[Xh; Yh; Zh]
# rotor blade profile CD parameters
DELTA0=0.01
DELTA2=250.0

# -------------------------------- TAIL ROTOR ----------------------------------

# tail rotor speed [rad/s]
OMEGTR=124.62
# tail rotor radius [ft]
RTR=5.5
# tail rotor blade lift slope [1/rad]
A0TR=5.73
# tail rotor solidity
SOLTR=0.1875
# tail rotor twist [rad]
TWISTTR=-17.2*pi/180.0
# location (station, water line, butt line) [in]
FSTR=732.0
BLTR=-14.0
WLTR=324.7
# moment arms from fuselage CG to tail rotor hub
XTR=-(FSTR-FSCGB)/12.0
YTR=(BLTR-BLCGB)/12.0
ZTR=-(WLTR-WLCGB)/12.0
# location of tail rotor relative to CG [ft]
rTR=[XTR; YTR; ZTR]

# -------------------------- HORIZONTAL STABILIZER -----------------------------

# area [ft²]
SHT=45.0
# lift slope [1/rad]
AHT=2.3
# horizontal stabilizer location (station, water line, butt line) [in]
FSHT=700.0
BLHT=0.0
WLHT=244.0
# moment arms from fuselage CG to horizontal stabilizer aero center
XHT=-(FSHT-FSCG)/12.0
YHT=(BLHT-BLCG)/12.0
ZHT=-(WLHT-WLCG)/12.0
# location of tail rotor relative to CG [ft]
rHT=[XHT; YHT; ZHT]

# --------------------------- VERTICAL STABILIZER ------------------------------

# area [ft²]
SVT=32.3
# lift slope [1/rad]
AVT=2.3
# horizontal stabilizer location (station, water line, butt line) [in]
FSVT=695.0
BLVT=0.0
WLVT=273.0
# moment arms from fuselage CG to vertical stabilizer aero center
XVT=-(FSVT-FSCG)/12.0
YVT=(BLVT-BLCG)/12.0
ZVT=-(WLVT-WLCG)/12.0
# location of tail rotor relative to CG [ft]
rVT=[XVT; YVT; ZVT]

# ------------------------------ CONTROL MIXING --------------------------------

# conversion from stick percentages to stick inches
perc2in = [0.1 0   0   0
           0   0.1 0   0
           0   0   0.1 0
           0   0   0   5.38/100]
# mixer gain matrix
LNKGAIN = [0.2062 0      0      0
           0      0.2172 0      0
           0      0      0.2025 0
           0      0      0      0.3780]
# convert stick inputs to servo commands
MIXGAIN = [0.6800  0      0.9930   0
           0       1.2020 0.9050  -0.4288
           0      -1.2020 1.2990   0.4288
           0       0      -0.8554  1.6043]
# gain for conversion servo to blade deflections
SWASHGAIN = [11.3413 -5.6706 -5.6706
              0      -5.6706  5.6706
              0       3.5034  3.5034]
TRGAIN = -9.142
# bias for conversion servo to blade deflections
SWASHBIAS = [-7.9738, 9.7280, 10.1111]
TRBIAS = 21.98

# --------------------- LINEARIZATION AND TRIM PARAMETERS ----------------------

# number of states
NSTATES=12
# number of control inputs
NCTRLS=4
# number of output variables
NOUT=6
# state vector perturbations for linearization
DELXLIN=[0.1; 0.1; 0.1; pi/180*0.1*ones(6); 0.1; 0.1; 0.1]
# control vector perturbations for linearization
DELCLIN=[0.1; 0.1; 0.1; 0.1]
# trim targets
TRIMTARG=collect(1:12)
# trim variables
TRIMVARS=[collect(1:8); collect(13:16)]
# tolerance on trim error
TOL=0.001*[0.001; 0.001; 0.001; pi/180*0.001*ones(6); 0.001; 0.001; 0.001;
           0.001; 0.001; 0.001]

max_stars_repo_name: umbertosaetti/J-SimpleHel

max_issues_repo_name: umbertosaetti/J-SimpleHel

max_forks_repo_name: umbertosaetti/J-SimpleHel


text:
function inter(mach::Vector{Float64},coeff::Vector{Float64})
    coeff_inter = LinearInterpolation(mach, coeff, extrapolation_bc=Flat())
    return coeff_inter

end

function inter_matrix(mach::Vector{Float64},coeff::Array{Float64, 2},α::Vector{Float64})
    coeff_inter = LinearInterpolation((mach,α), coeff,extrapolation_bc=Flat())
    return coeff_inter

end

function drag(CD0,CDδ2,αₜ,mach)
    CD = CD0(mach) + CDδ2(mach) * (sin(αₜ))^2
    return CD
end
function drag(CX0,CX2,αₜ,mach,CNa)
    CX = CX0(mach) + CX2(mach) * (sin(αₜ))^2
    #CD = CNa(mach) *(sin(αₜ))^2 - CX*cos(αₜ)
    CD = CNa(mach) *sin(αₜ) + CX*cos(αₜ)
    return CD
end
function drag(CX0,CX2,αₜ, α,β,mach,CNa)
    CX = CX0(mach) + CX2(mach) * (sin(αₜ))^2
    #CD = CNa(mach) *(sin(αₜ))^2 - CX*cos(αₜ)
    CD = CNa(mach) *sin(α)*cos(β) + CX*cos(α)*cos(β)
    return CD
end

function lift(CLα,CLα³,αₜ,mach)
    CL = CLα + CLα³ * (sin(αₜ))^2
    return CL
end
function lift(CX0,CX2,αₜ,mach,CNa)
    CX = CX0(mach) + CX2(mach) * (sin(αₜ))^2
    #CL = CNa(mach) *cos(αₜ) + CX
    CL =-CNa(mach)*cos(αₜ)+CX*sin(αₜ)
    return CL
end

function matrix(m,args...)

for arg in args
    m=hcat(m,arg)
end
return m
end


text:
"""
    $SIGNATURES

A lazy representation of a OMETIFF file. This custom type is needed since TIFF
files are laid out noncontiguously and nonregularly. It uses an internal index
to determine the mapping from indices to the locations of data slices on disk.
These slices are generally XY slices and are usually loaded in all at once so it
is quickly loaded into an internal cache to speed up the process. Externally,
this type should behave very similarly to an in-memory array, albeit with a
higher cost of accessing an element.

$(FIELDS)
"""
mutable struct DiskOMETaggedImage{T, N1, N2, O, AA <: AbstractArray} <: AbstractDenseTIFF{T, N2}
    """
    A map of dimensions (sans XY) to the corresponding [`TiffFile`](@ref) and [`IFD`](@ref)
    """
    ifds::OrderedDict{NTuple{N1, Int}, Tuple{TiffFile{O}, IFD{O}}}

    """
    The full set of dimensions of the TIFF file, including XY
    """
    dims::NTuple{N2, Int}

    """
    An internal cache to fill when reading from disk
    """
    cache::AA

    """
    The dimension indices corresponding to the slice currently in the cache
    """
    cache_index::NTuple{N1, Int}

    function DiskOMETaggedImage(ifds::OrderedDict{NTuple{N1, Int}, Tuple{TiffFile{O}, IFD{O}}}, dims::NTuple{N2, Int}) where {N1, O, N2}
        if N2 - 2 != N1
            error("$N2 dimensions given, but the IFDs are indexed on $N1 dimensions instead of "*
                  "expected $(N2-2).")
        end
        ifd = first(values(ifds))[2]
        cache = getcache(ifd)
        new{eltype(cache), N1, N2, O, typeof(cache)}(ifds, dims, cache, (-1, -1, -1, -1))
    end
end

Base.size(A::DiskOMETaggedImage) = A.dims

function Base.getindex(A::DiskOMETaggedImage{T, N1, N2, O, AA}, i1::Int, i2::Int, i::Vararg{Int, N1}) where {T, N1, N2, O, AA}
    # check the loaded cache is already the correct slice
    if A.cache_index == i
        return A.cache[i2, i1]
    end

    file, ifd = A.ifds[i]

    # if the file isn't open, lets open a handle and update it
    if !isopen(file.io)
        path = file.filepath
        file.io = getstream(format"OMETIFF", open(path), path)
    end

    read!(A.cache, file, ifd)

    A.cache_index = i

    return A.cache[i2, i1]
end

function Base.setindex!(A::T, X, I...) where {T <: DiskOMETaggedImage}
    error("This array is on disk and is read only. Convert to a mutable in-memory version by running "*
          "`copy(arr)`. \n\n𝗡𝗼𝘁𝗲: For large files this can be quite expensive. A future PR will add "*
          "support for reading and writing to/from disk.")
end


text:
## Calculus

Base.sum{C<:Chebyshev,DD}(f::Fun{JacobiWeight{C,DD}})=sum(setdomain(f,canonicaldomain(f))*fromcanonicalD(f))
linesum{C<:Chebyshev,DD}(f::Fun{JacobiWeight{C,DD}})=linesum(setdomain(f,canonicaldomain(f))*abs(fromcanonicalD(f)))

for (Func,Len) in ((:(Base.sum),:complexlength),(:linesum,:arclength))
    @eval begin
        function $Func{C<:Chebyshev,DD<:Interval}(f::Fun{JacobiWeight{C,DD}})
            tol=1e-10
            d,α,β,n=domain(f),f.space.α,f.space.β,ncoefficients(f)
            g=Fun(f.coefficients,space(f).space)
            if α ≤ -1.0 && abs(first(g))≤tol
                $Func(increase_jacobi_parameter(-1,f))
            elseif β ≤ -1.0 && abs(last(g))≤tol
                $Func(increase_jacobi_parameter(+1,f))
            elseif α ≤ -1.0 || β ≤ -1.0
                fs = Fun(f.coefficients,f.space.space)
                return Inf*0.5*$Len(d)*(sign(fs(d.a))+sign(fs(d.b)))/2
            elseif α == β == -0.5
                return 0.5*$Len(d)*f.coefficients[1]*π
            elseif α == β == 0.5
                return 0.5*$Len(d)*(n ≤ 2 ? f.coefficients[1]/2 : f.coefficients[1]/2 - f.coefficients[3]/4)*π
            elseif α == 0.5 && β == -0.5
                return 0.5*$Len(d)*(n == 1 ? f.coefficients[1] : f.coefficients[1] + f.coefficients[2]/2)*π
            elseif α == -0.5 && β == 0.5
                return 0.5*$Len(d)*(n == 1 ? f.coefficients[1] : f.coefficients[1] - f.coefficients[2]/2)*π
            else
                c = zeros(eltype(f),n)
                c[1] = 2.^(α+β+1)*gamma(α+1)*gamma(β+1)/gamma(α+β+2)
                if n > 1
                    c[2] = c[1]*(α-β)/(α+β+2)
                    for i=1:n-2
                        c[i+2] = (2(α-β)*c[i+1]-(α+β-i+2)*c[i])/(α+β+i+2)
                    end
                end
                return 0.5*$Len(d)*dotu(f.coefficients,c)
            end
        end
        $Func{PS<:PolynomialSpace,DD}(f::Fun{JacobiWeight{PS,DD}}) = $Func(Fun(f,
                                                                       JacobiWeight(space(f).α,
                                                                                space(f).β,
                                                                                Chebyshev(domain(f)))))
    end
end

function differentiate{J<:JacobiWeight,DD<:Interval}(f::Fun{J,DD})
    S=f.space
    d=domain(f)
    ff=Fun(f.coefficients,S.space)
    if S.α==S.β==0
        u=differentiate(ff)
        Fun(u.coefficients,JacobiWeight(0.,0.,space(u)))
    elseif S.α==0
        x=Fun(identity,d)
        M=tocanonical(d,x)
        Mp=tocanonicalD(d,d.a)
        u=-Mp*S.β*ff +(1-M).*differentiate(ff)
        Fun(u.coefficients,JacobiWeight(0.,S.β-1,space(u)))
    elseif S.β==0
        x=Fun(identity,d)
        M=tocanonical(d,x)
        Mp=tocanonicalD(d,d.a)
        u=Mp*S.α*ff +(1+M).*differentiate(ff)
        Fun(u.coefficients,JacobiWeight(S.α-1,0.,space(u)))
    else
        x=Fun(identity,d)
        M=tocanonical(d,x)
        Mp=tocanonicalD(d,d.a)
        u=(Mp*S.α)*(1-M).*ff- (Mp*S.β)*(1+M).*ff +(1-M.^2).*differentiate(ff)
        Fun(u.coefficients,JacobiWeight(S.α-1,S.β-1,space(u)))
    end
end

function integrate{SS,DD<:Interval}(f::Fun{JacobiWeight{SS,DD}})
    S=space(f)
    # we integrate by solving u'=f
    D=Derivative(S)
    tol=1e-10
    g=Fun(f.coefficients,S.space)
    if isapprox(S.α,0.) && isapprox(S.β,0.)
        integrate(g)
    elseif S.α ≤ -1.0 && abs(first(g))≤tol
        integrate(increase_jacobi_parameter(-1,f))
    elseif S.β ≤ -1.0 && abs(last(g))≤tol
        integrate(increase_jacobi_parameter(+1,f))
    elseif isapprox(S.α,-1) && isapprox(S.β,-1)
        error("Implement")
    elseif isapprox(S.α,-1) && isapprox(S.β,0)
        p=first(g)  # first value without weight
        fp = Fun(f-Fun([p],S),S.space)  # Subtract out right value and divide singularity via conversion
        d=domain(f)
        Mp=tocanonicalD(d,d.a)
        integrate(fp)⊕Fun([p/Mp],LogWeight(1.,0.,S.space))
    elseif isapprox(S.α,-1) && S.β > 0 && isapproxinteger(S.β)
        # convert to zero case and integrate
        integrate(Fun(f,JacobiWeight(S.α,0.,S.space)))
    elseif isapprox(S.β,-1) && isapprox(S.α,0.)
        p=last(g)  # last value without weight
        fp = Fun(f-Fun([p],S),S.space)  # Subtract out right value and divide singularity via conversion
        d=domain(f)
        Mp=tocanonicalD(d,d.a)
        integrate(fp)⊕Fun([-p/Mp],LogWeight(0.,1.,S.space))
    elseif isapprox(S.β,-1) && S.α > 0 && isapproxinteger(S.α)
        # convert to zero case and integrate
        integrate(Fun(f,JacobiWeight(0.,S.β,S.space)))
    elseif isapprox(S.α,0) || isapprox(S.β,0)
        D\f   # this happens to pick out a smooth solution
    else
        s=sum(f)
        if abs(s)<1E-14
            linsolve(D,f;tolerance=1E-14)  # if the sum is 0 we don't get step-like behaviour
        else
            # we normalized so it sums to zero, and so backslash works
            w=Fun(x->exp(-40x^2),81)
            w1=Fun(coefficients(w),S)
            w2=Fun(x->w1(x),domain(w1))
            c=s/sum(w1)
            v=f-w1*c
            (c*integrate(w2))⊕linsolve(D,v;tolerance=100eps())
        end
    end
end

function Base.cumsum{SS,DD<:Interval}(f::Fun{JacobiWeight{SS,DD}})
    g=integrate(f)
    S=space(f)

    if (S.α==0 && S.β==0) || S.α>-1
        g-first(g)
    else
        warn("Function is not integrable at left endpoint.  Returning a non-normalized indefinite integral.")
        g
    end
end

## Operators

function jacobiweightDerivative{SS,DDD<:Interval}(S::JacobiWeight{SS,DDD})
    d=domain(S)

    if d!=Interval()
        # map to canonical
        Mp=fromcanonicalD(d,d.a)
        DD=jacobiweightDerivative(setdomain(S,Interval()))

        return DerivativeWrapper(SpaceOperator(DD.op.op,S,setdomain(rangespace(DD),d))/Mp,1)
    end

    if S.α==S.β==0
        DerivativeWrapper(SpaceOperator(Derivative(S.space),S,JacobiWeight(0.,0.,rangespace(Derivative(S.space)))),1)
    elseif S.α==0
        w=Fun([1.],JacobiWeight(0,1,ConstantSpace(d)))

        DD=-S.β + w*Derivative(S.space)
        rs=S.β==1?rangespace(DD):JacobiWeight(0.,S.β-1,rangespace(DD))
        DerivativeWrapper(SpaceOperator(DD,S,rs),1)
    elseif S.β==0
        w=Fun([1.],JacobiWeight(1,0,ConstantSpace(d)))

        DD=S.α + w*Derivative(S.space)
        rs=S.α==1?rangespace(DD):JacobiWeight(S.α-1,0.,rangespace(DD))
        DerivativeWrapper(SpaceOperator(DD,S,rs),1)
    else
        w=Fun([1.],JacobiWeight(1,1,ConstantSpace(d)))
        x=Fun()

        DD=S.α*(1-x) - S.β*(1+x) + w*Derivative(S.space)
        rs=S.α==1&&s.β==1?rangespace(DD):JacobiWeight(S.α-1,S.β-1,rangespace(DD))
        DerivativeWrapper(SpaceOperator(DD,S,rs),1)
    end
end

Derivative{SS,DDD<:Interval}(S::JacobiWeight{SS,DDD})=jacobiweightDerivative(S)

function Derivative{SS,DD<:Interval}(S::JacobiWeight{SS,DD},k::Integer)
    if k==1
        Derivative(S)
    else
        D=Derivative(S)
        DerivativeWrapper(TimesOperator(Derivative(rangespace(D),k-1),D),k)
    end
end

## Multiplication

#Left multiplication. Here, S is considered the domainspace and we determine rangespace accordingly.

function Multiplication{S1,S2,DD<:IntervalDomain,T}(f::Fun{JacobiWeight{S1,DD},T},S::JacobiWeight{S2,DD})
    M=Multiplication(Fun(f.coefficients,space(f).space),S.space)
    if space(f).α+S.α==space(f).β+S.β==0
        rsp=rangespace(M)
    else
        rsp=JacobiWeight(space(f).α+S.α,space(f).β+S.β,rangespace(M))
    end
    MultiplicationWrapper(f,SpaceOperator(M,S,rsp))
end

function Multiplication{D,T,SS,DD<:IntervalDomain}(f::Fun{D,T},S::JacobiWeight{SS,DD})
    M=Multiplication(f,S.space)
    rsp=JacobiWeight(S.α,S.β,rangespace(M))
    MultiplicationWrapper(f,SpaceOperator(M,S,rsp))
end

function Multiplication{SS,T,ID<:IntervalDomain}(f::Fun{JacobiWeight{SS,ID},T},S::PolynomialSpace{ID})
    M=Multiplication(Fun(f.coefficients,space(f).space),S)
    rsp=JacobiWeight(space(f).α,space(f).β,rangespace(M))
    MultiplicationWrapper(f,SpaceOperator(M,S,rsp))
end

#Right multiplication. Here, S is considered the rangespace and we determine domainspace accordingly.

function Multiplication{DD<:IntervalDomain,SS,S2,T}(S::JacobiWeight{SS,DD},f::Fun{JacobiWeight{S2,DD},T})
    M=Multiplication(Fun(f.coefficients,space(f).space),S.space)
    dsp=canonicalspace(JacobiWeight(S.α-space(f).α,S.β-space(f).β,rangespace(M)))
    MultiplicationWrapper(f,SpaceOperator(M,dsp,S))
end

function Multiplication{D,SS,DD<:IntervalDomain,T}(S::JacobiWeight{SS,DD},f::Fun{D,T})
    M=Multiplication(f,S.space)
    dsp=JacobiWeight(S.α,S.β,rangespace(M))
    MultiplicationWrapper(f,SpaceOperator(M,dsp,S))
end

# function Multiplication{D<:JacobiWeight,T,V,ID<:IntervalDomain}(S::Space{V,D},f::Fun{ID,T})
#     M=Multiplication(Fun(f.coefficients,space(f).space),S)
#     dsp=JacobiWeight(-space(f).α,-space(f).β,rangespace(M))
#     MultiplicationWrapper(f,SpaceOperator(M,dsp,S))
# end

## Conversion

# The second case handles zero
isapproxinteger(x) = isapprox(x,round(Int,x))  || isapprox(x+1,round(Int,x+1))

for (OPrule,OP) in ((:maxspace_rule,:maxspace),(:union_rule,:union))
    @eval begin
        function $OPrule(A::JacobiWeight,B::JacobiWeight)
            if domainscompatible(A,B) && isapproxinteger(A.α-B.α) && isapproxinteger(A.β-B.β)
                ms=$OP(A.space,B.space)
                if min(A.α,B.α)==0.0 && min(A.β,B.β) == 0.0
                    return ms
                else
                    return JacobiWeight(min(A.α,B.α),min(A.β,B.β),ms)
                end
            end
            NoSpace()
        end
        $OPrule{T,D<:IntervalDomain}(A::JacobiWeight,B::Space{T,D})=$OPrule(A,JacobiWeight(0.,0.,B))
    end
end

for FUNC in (:hasconversion,:isconvertible)
    @eval begin
        $FUNC{S1,S2,D<:IntervalDomain}(A::JacobiWeight{S1,D},B::JacobiWeight{S2,D})=isapproxinteger(A.α-B.α) &&
            isapproxinteger(A.β-B.β) && A.α ≥ B.α && A.β ≥ B.β && $FUNC(A.space,B.space)

        $FUNC{T,S,D<:IntervalDomain}(A::JacobiWeight{S,D},B::Space{T,D})=$FUNC(A,JacobiWeight(0.,0.,B))
        $FUNC{T,S,D<:IntervalDomain}(B::Space{T,D},A::JacobiWeight{S,D})=$FUNC(JacobiWeight(0.,0.,B),A)
    end
end

# return the space that has banded Conversion to the other, or NoSpace

function conversion_rule(A::JacobiWeight,B::JacobiWeight)
    if isapproxinteger(A.α-B.α) && isapproxinteger(A.β-B.β)
        ct=conversion_type(A.space,B.space)
        ct==NoSpace()?NoSpace():JacobiWeight(max(A.α,B.α),max(A.β,B.β),ct)
    else
        NoSpace()
    end
end

conversion_rule{T,D<:IntervalDomain}(A::JacobiWeight,B::UnivariateSpace{T,D})=conversion_type(A,JacobiWeight(0,0,B))

# override defaultConversion instead of Conversion to avoid ambiguity errors
function defaultConversion{JS1,JS2,DD<:IntervalDomain}(A::JacobiWeight{JS1,DD},B::JacobiWeight{JS2,DD})
    @assert isapproxinteger(A.α-B.α) && isapproxinteger(A.β-B.β)

    if isapprox(A.α,B.α) && isapprox(A.β,B.β)
        ConversionWrapper(SpaceOperator(Conversion(A.space,B.space),A,B))
    else
        @assert A.α≥B.α&&A.β≥B.β
        # first check if a multiplication by JacobiWeight times B.space is overloaded
        # this is currently designed for Jacobi multiplied by (1-x), etc.
        αdif,βdif=round(Int,A.α-B.α),round(Int,A.β-B.β)
        d=domain(A)
        M=Multiplication(jacobiweight(αdif,βdif,d),
                         A.space)

        if rangespace(M)==JacobiWeight(αdif,βdif,A.space)
            # M is the default, so we should use multiplication by polynomials instead
            x=Fun(identity,d)
            y=mobius(d,x)   # we use mobius instead of tocanonical so that it works for Funs
            m=(1+y).^αdif.*(1-y).^βdif
            MC=promoterangespace(Multiplication(m,A.space),B.space)

            ConversionWrapper(SpaceOperator(MC,A,B))# Wrap the operator with the correct spaces
        else
            ConversionWrapper(SpaceOperator(promoterangespace(M,B.space),
                                            A,B))
        end
    end
end

defaultConversion{JS,D<:IntervalDomain}(A::RealUnivariateSpace{D},B::JacobiWeight{JS,D})=ConversionWrapper(
    SpaceOperator(
        Conversion(JacobiWeight(0,0,A),B),
        A,B))
defaultConversion{JS,D<:IntervalDomain}(A::JacobiWeight{JS,D},B::RealUnivariateSpace{D})=ConversionWrapper(
    SpaceOperator(
        Conversion(A,JacobiWeight(0,0,B)),
        A,B))

## Evaluation

function  Base.getindex{J<:JacobiWeight}(op::ConcreteEvaluation{J,Bool},kr::Range)
    S=op.space
    @assert op.order<=1
    d=domain(op)

    if op.x
        @assert S.β>=0
        if S.β==0
            if op.order==0
                2^S.α*getindex(Evaluation(S.space,op.x),kr)
            else #op.order ===1
                @assert isa(d,IntervalDomain)
                2^S.α*getindex(Evaluation(S.space,op.x,1),kr)+
                    (tocanonicalD(d,d.a)*S.α*2^(S.α-1))*getindex(Evaluation(S.space,op.x),kr)
            end
        else
            @assert op.order==0
            zeros(kr)
        end
    else
        @assert S.α>=0
        if S.α==0
            if op.order==0
                2^S.β*getindex(Evaluation(S.space,op.x),kr)
            else #op.order ===1
                @assert isa(d,IntervalDomain)
                2^S.β*getindex(Evaluation(S.space,op.x,1),kr)-(tocanonicalD(d,d.a)*S.β*2^(S.β-1))*getindex(Evaluation(S.space,op.x),kr)
            end
        else
            @assert op.order==0
            zeros(kr)
        end
    end
end

## Definite Integral

for (Func,Len) in ((:DefiniteIntegral,:complexlength),(:DefiniteLineIntegral,:arclength))
    ConcFunc = parse("Concrete"*string(Func))

    @eval begin
        function getindex{LT,D<:Interval,T}(Σ::$ConcFunc{JacobiWeight{Ultraspherical{LT,D},D},T},k::Integer)
            λ = order(domainspace(Σ))
            dsp = domainspace(Σ)
            d = domain(Σ)
            C = $Len(d)/2

            if dsp.α==dsp.β==λ-0.5
                k == 1? C*gamma(λ+one(T)/2)*gamma(one(T)/2)/gamma(λ+one(T)) : zero(T)
            else
                sum(Fun([zeros(k-1);1],dsp))
            end
        end

        function getindex{LT,D<:Interval,T}(Σ::$ConcFunc{JacobiWeight{Ultraspherical{LT,D},D},T},kr::Range)
            λ = order(domainspace(Σ))
            dsp = domainspace(Σ)
            d = domain(Σ)
            C = $Len(d)/2

            if dsp.α==dsp.β==λ-0.5
                promote_type(T,typeof(C))[k == 1? C*gamma(λ+one(T)/2)*gamma(one(T)/2)/gamma(λ+one(T)) : zero(T) for k=kr]
            else
                promote_type(T,typeof(C))[sum(Fun([zeros(k-1);1],dsp)) for k=kr]
            end
        end

        function bandinds{LT,D<:Interval}(Σ::$ConcFunc{JacobiWeight{Ultraspherical{LT,D},D}})
            λ = order(domainspace(Σ))
            α,β = domainspace(Σ).α,domainspace(Σ).β
            if α==β && isapproxinteger(α-0.5-λ) && λ ≤ ceil(Int,α)
                0,2*(ceil(Int,α)-λ)
            else
                0,∞
            end
        end

        function getindex{D<:Interval,T}(Σ::$ConcFunc{JacobiWeight{Chebyshev{D},D},T},k::Integer)
            dsp = domainspace(Σ)
            d = domain(Σ)
            C = $Len(d)/2

            if dsp.α==dsp.β==-0.5
                k == 1? C*π : zero(T)
            else
                sum(Fun([zeros(k-1);1],dsp))
            end
        end

        function getindex{D<:Interval,T}(Σ::$ConcFunc{JacobiWeight{Chebyshev{D},D},T},kr::Range)
            dsp = domainspace(Σ)
            d = domain(Σ)
            C = $Len(d)/2

            if dsp.α==dsp.β==-0.5
                promote_type(T,typeof(C))[k == 1? C*π : zero(T) for k=kr]
            else
                promote_type(T,typeof(C))[sum(Fun([zeros(k-1);1],dsp)) for k=kr]
            end
        end

        function bandinds{D<:Interval}(Σ::$ConcFunc{JacobiWeight{Chebyshev{D},D}})
            α,β = domainspace(Σ).α,domainspace(Σ).β
            if α==β && isapproxinteger(α-0.5) && 0 ≤ ceil(Int,α)
                0,2ceil(Int,α)
            else
                0,∞
            end
        end
    end
end


text:
"""
    count_nucleotides(strand)

The count of each nucleotide within `strand` as a dictionary.

Invalid strands raise a `DomainError`.

"""
function count_nucleotides(strand)
    valid_chars = Set(['A', 'C', 'G', 'T'])
    counts = Dict(c=>0 for c in valid_chars)
    foreach(c -> counts[c] = get(counts, c, 0) + 1, strand)
    issubset(keys(counts), valid_chars) ? counts : throw(DomainError("Invalid strand"))
end

max_stars_repo_name: koladilip/julia-exercism

max_issues_repo_name: koladilip/julia-exercism

max_forks_repo_name: koladilip/julia-exercism