text:
"""
    step!(integ::DEIntegrator [, dt [, stop_at_tdt]])
Perform one (successful) step on the integrator.

Alternative, if a `dt` is given, then `step!` the integrator until
there is a temporal difference `≥ dt` in `integ.t`.  When `true` is
passed to the optional third argument, the integrator advances exactly
`dt`.
"""
function step!(d::DEIntegrator) error("Integrator stepping is not implemented") end

"""
    resize!(integrator::DEIntegrator,k::Int)

Resizes the DE to a size `k`. This chops off the end of the array, or adds blank values at the end, depending on whether
`k > length(integrator.u)`.
"""
Base.resize!(i::DEIntegrator,ii::Int) = error("resize!: method has not been implemented for the integrator")

"""
    deleteat!(integrator::DEIntegrator,idxs)

Shrinks the ODE by deleting the `idxs` components.
"""
Base.deleteat!(i::DEIntegrator,ii) = error("deleteat!: method has not been implemented for the integrator")

"""
    addat!(integrator::DEIntegrator,idxs,val)

Grows the ODE by adding the `idxs` components. Must be contiguous indices.
"""
addat!(i::DEIntegrator,idxs,val=zeros(length(idxs))) = error("addat!: method has not been implemented for the integrator")

"""
    get_tmp_cache(i::DEIntegrator)

Returns a tuple of internal cache vectors which are safe to use as temporary arrays. This should be used
for integrator interface and callbacks which need arrays to write into in order to be non-allocating.
The length of the tuple is dependent on the method.
"""
get_tmp_cache(i::DEIntegrator) = error("get_tmp_cache!: method has not been implemented for the integrator")
user_cache(i::DEIntegrator) = error("user_cache: method has not been implemented for the integrator")
u_cache(i::DEIntegrator) = error("u_cache: method has not been implemented for the integrator")
du_cache(i::DEIntegrator) = error("du_cache: method has not been implemented for the integrator")
ratenoise_cache(i::DEIntegrator) = ()
rand_cache(i::DEIntegrator) = ()

"""
    full_cache(i::DEIntegrator)

Returns an iterator over the cache arrays of the method. This can be used to change internal values as needed.
"""
full_cache(i::DEIntegrator) = error("full_cache: method has not been implemented for the integrator")

"""
    resize_non_user_cache!(integrator::DEIntegrator,k::Int)

Resizes the non-user facing caches to be compatible with a DE of size `k`. This includes resizing Jacobian caches.

!!! note
    In many cases, [`resize!`](@ref) simply resizes [`full_cache`](@ref) variables and then
    calls this function. This finer control is required for some `AbstractArray`
    operations.
"""
resize_non_user_cache!(i::DEIntegrator,ii::Int) = error("resize_non_user_cache!: method has not been implemented for the integrator")

"""
    deleteat_non_user_cache!(integrator::DEIntegrator,idxs)

[`deleteat!`](@ref)s the non-user facing caches at indices `idxs`. This includes resizing Jacobian caches.

!!! note
    In many cases, `deleteat!` simply `deleteat!`s [`full_cache`](@ref) variables and then
    calls this function. This finer control is required for some `AbstractArray`
    operations.
"""
deleteat_non_user_cache!(i::DEIntegrator,idxs) = error("deleteat_non_user_cache!: method has not been implemented for the integrator")

"""
    addat_non_user_cache!(i::DEIntegrator,idxs)

[`addat!`](@ref)s the non-user facing caches at indices `idxs`. This includes resizing Jacobian caches.
!!! note
    In many cases, `addat!` simply `addat!`s [`full_cache`](@ref) variables and then
    calls this function. This finer control is required for some `AbstractArray`
    operations.
"""
addat_non_user_cache!(i::DEIntegrator,idxs) = error("addat_non_user_cache!: method has not been implemented for the integrator")

"""
    terminate!(i::DEIntegrator[, retcode = :Terminated])

Terminates the integrator by emptying `tstops`. This can be used in events and callbacks to immediately
end the solution process.  Optionally, `retcode` may be specified (see: [Return Codes (RetCodes)](@ref retcodes)).
"""
terminate!(i::DEIntegrator) = error("terminate!: method has not been implemented for the integrator")

"""
    get_du(i::DEIntegrator)

Returns the derivative at `t`.
"""
get_du(i::DEIntegrator) = error("get_du: method has not been implemented for the integrator")

"""
    get_du!(out,i::DEIntegrator)

Write the current derivative at `t` into `out`.
"""
get_du!(out,i::DEIntegrator) = error("get_du: method has not been implemented for the integrator")
get_dt(i::DEIntegrator) = error("get_dt: method has not been implemented for the integrator")

"""
    get_proposed_dt(i::DEIntegrator)

Gets the proposed `dt` for the next timestep.
"""
get_proposed_dt(i::DEIntegrator) = error("get_proposed_dt: method has not been implemented for the integrator")

"""
    set_proposed_dt(i::DEIntegrator,dt)
    set_proposed_dt(i::DEIntegrator,i2::DEIntegrator)

Sets the proposed `dt` for the next timestep. If second argument isa `DEIntegrator` then it sets the timestepping of
first argument to match that of second one. Note that due to PI control and step acceleration this is more than matching
the factors in most cases.
"""
set_proposed_dt!(i::DEIntegrator,dt) = error("set_proposed_dt!: method has not been implemented for the integrator")

"""
    savevalues!(integrator::DEIntegrator,
      force_save=false) -> Tuple{Bool, Bool}

Try to save the state and time variables at the current time point, or the
`saveat` point by using interpolation when appropriate. It returns a tuple that
is `(saved, savedexactly)`. If `savevalues!` saved value, then `saved` is true,
and if `savevalues!` saved at the current time point, then `savedexactly` is
true.

The saving priority/order is as follows:
  - `save_on`
    - `saveat`
    - `force_save`
    - `save_everystep`
"""
savevalues!(i::DEIntegrator) = error("savevalues!: method has not been implemented for the integrator")

"""
    u_modified!(i::DEIntegrator,bool)

Sets `bool` which states whether a change to `u` occurred, allowing the solver to handle the discontinuity. By default,
this is assumed to be true if a callback is used. This will result in the re-calculation of the derivative at
`t+dt`, which is not necessary if the algorithm is FSAL and `u` does not experience a discontinuous change at the
end of the interval. Thus if `u` is unmodified in a callback, a single call to the derivative calculation can be
eliminated by `u_modified!(integrator,false)`.
"""
u_modified!(i::DEIntegrator,bool) = error("u_modified!: method has not been implemented for the integrator")

"""
    add_tstop!(i::DEIntegrator,t)

Adds a `tstop` at time `t`.
"""
add_tstop!(i::DEIntegrator,t) = error("add_tstop!: method has not been implemented for the integrator")

"""
    add_saveat!(i::DEIntegrator,t)

Adds a `saveat` time point at `t`.
"""
add_saveat!(i::DEIntegrator,t) = error("add_saveat!: method has not been implemented for the integrator")

set_abstol!(i::DEIntegrator,t) = error("set_abstol!: method has not been implemented for the integrator")
set_reltol!(i::DEIntegrator,t) = error("set_reltol!: method has not been implemented for the integrator")

"""
    reinit!(integrator::DEIntegrator,args...; kwargs...)

The reinit function lets you restart the integration at a new value.

# Arguments

- `u0`: Value of `u` to start at. Default value is `integrator.sol.prob.u0`

# Keyword Arguments
- `t0`: Starting timepoint. Default value is `integrator.sol.prob.tspan[1]`
- `tf`: Ending timepoint. Default value is `integrator.sol.prob.tspan[2]`
- `erase_sol=true`: Whether to start with no other values in the solution, or keep the previous solution.
- `tstops`, `d_discontinuities`, & `saveat`: Cache where these are stored. Default is the original cache.
- `reset_dt`: Set whether to reset the current value of `dt` using the automatic `dt` determination algorithm. Default is
  `(integrator.dtcache == zero(integrator.dt)) && integrator.opts.adaptive`
- `reinit_callbacks`: Set whether to run the callback initializations again (and `initialize_save` is for that). Default is `true`.
- `reinit_cache`: Set whether to re-run the cache initialization function (i.e. resetting FSAL, not allocating vectors)
  which should usually be true for correctness. Default is `true`.

Additionally, once can access [`auto_dt_reset!`](@ref) which will run the auto `dt` initialization algorithm.
"""
reinit!(integrator::DEIntegrator,args...; kwargs...) =
       error("reinit!: method has not been implemented for the integrator")

"""
initialize_dae!(integrator::DEIntegrator,initializealg = integrator.initializealg)

Runs the DAE initialization to find a consistent state vector. The optional
argument `initializealg` can be used to specify a different initialization
algorithm to use.
"""
initialize_dae!(integrator::DEIntegrator) =
       error("initialize_dae!: method has not been implemented for the integrator")

"""
    auto_dt_reset!(integrator::DEIntegrator)

Run the auto `dt` initialization algorithm.
"""
auto_dt_reset!(integrator::DEIntegrator) = error("auto_dt_reset!: method has not been implemented for the integrator")

"""
    change_t_via_interpolation!(integrator::DEIntegrator,t,modify_save_endpoint=Val{false})

Modifies the current `t` and changes all of the corresponding values using the local interpolation. If the current solution
has already been saved, one can provide the optional value `modify_save_endpoint` to also modify the endpoint of `sol` in the
same manner.
"""
change_t_via_interpolation!(i::DEIntegrator,args...) = error("change_t_via_interpolation!: method has not been implemented for the integrator")

addsteps!(i::DEIntegrator,args...) = nothing

"""
    reeval_internals_due_to_modification!(integrator::DDEIntegrator)

Recalculate interpolation data and update ODE integrator after changes by callbacks.
"""
reeval_internals_due_to_modification!(integrator::DEIntegrator) = nothing

"""
    set_t!(integrator::DEIntegrator, t)

Set current time point of the `integrator` to `t`.
"""
set_t!(integrator::DEIntegrator, t) =
    error("set_t!: method has not been implemented for the integrator")

"""
    set_u!(integrator::DEIntegrator, u)

Set current state of the `integrator` to `u`.
"""
set_u!(integrator::DEIntegrator, u) =
    error("set_u!: method has not been implemented for the integrator")

"""
    set_ut!(integrator::DEIntegrator, u, t)

Set current state of the `integrator` to `u` and `t`
"""
function set_ut!(integrator::DEIntegrator, u, t)
  set_u!(integrator, u)
  set_t!(integrator, t)
end

### Addat isn't a real thing. Let's make it a real thing Gretchen

function addat!(a::AbstractArray,idxs,val=nothing)
  if val === nothing
    resize!(a,length(a)+length(idxs))
  else
    error("real addat! on arrays isn't supported yet")
    #=
    flip_range = last(idxs):-1:idxs.start
    @show idxs,flip_range
    splice!(a,flip_range,val)
    =#
  end
end

### Integrator traits

has_reinit(i::DEIntegrator) = false

### Display

Base.summary(I::DEIntegrator) = string(
                  TYPE_COLOR, nameof(typeof(I)),
                  NO_COLOR, " with uType ",
                  TYPE_COLOR, typeof(I.u),
                  NO_COLOR, " and tType ",
                  TYPE_COLOR, typeof(I.t), NO_COLOR)

function Base.show(io::IO, A::DEIntegrator)
  println(io,string("t: ",A.t))
  print(io,"u: ")
  show(io, A.u)
end
function Base.show(io::IO, m::MIME"text/plain", A::DEIntegrator)
  println(io,string("t: ",A.t))
  print(io,"u: ")
  show(io,m,A.u)
end
TreeViews.hastreeview(x::DEIntegrator) = true
function TreeViews.treelabel(io::IO,x::DEIntegrator,
                             mime::MIME"text/plain" = MIME"text/plain"())
  show(io,mime,Base.Text(Base.summary(x)))
end

### Error check (retcode)

last_step_failed(integrator::DEIntegrator) = false

"""
    check_error(integrator)

Check state of `integrator` and return one of the
[Return Codes](http://docs.juliadiffeq.org/dev/basics/solution.html#Return-Codes-(RetCodes)-1)
"""
function check_error(integrator::DEIntegrator)
  # This implementation is intended to be used for ODEIntegrator and
  # SDEIntegrator.
  if isnan(integrator.dt)
    if integrator.opts.verbose
      @warn("NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.")
    end
    return :DtNaN
  end
  if integrator.iter > integrator.opts.maxiters
    if integrator.opts.verbose
      @warn("Interrupted. Larger maxiters is needed.")
    end
    return :MaxIters
  end

  # The last part:
  # If you are close to the end, don't exit: let the user hit the end!
  # However, if we try that and the step fails, exit instead of infinite loop
  if !integrator.opts.force_dtmin && integrator.opts.adaptive &&
     abs(integrator.dt) <= abs(integrator.opts.dtmin) &&
     (((hasproperty(integrator,:opts) && hasproperty(integrator.opts,:tstops)) ?
     integrator.t + integrator.dt < integrator.tdir*first(integrator.opts.tstops) :
     true) || (hasproperty(integrator,:accept_step) && !integrator.accept_step))
    if integrator.opts.verbose
      @warn("dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.")
    end
    return :DtLessThanMin
  end
  if integrator.opts.unstable_check(integrator.dt,integrator.u,integrator.p,integrator.t)
    if integrator.opts.verbose
      @warn("Instability detected. Aborting")
    end
    return :Unstable
  end
  if last_step_failed(integrator)
    if integrator.opts.verbose
      @warn("Newton steps could not converge and algorithm is not adaptive. Use a lower dt.")
    end
    return :ConvergenceFailure
  end
  return :Success
end

function postamble! end

"""
    check_error!(integrator)

Same as `check_error` but also set solution's return code
(`integrator.sol.retcode`) and run `postamble!`.
"""
function check_error!(integrator::DEIntegrator)
  code = check_error(integrator)
  if code != :Success
    integrator.sol = solution_new_retcode(integrator.sol, code)
    postamble!(integrator)
  end
  return code
end

### Default Iterator Interface
function done(integrator::DEIntegrator)
  if ! (integrator.sol.retcode in (:Default, :Success))
    return true
  elseif isempty(integrator.opts.tstops)
    postamble!(integrator)
    return true
  elseif integrator.just_hit_tstop
    integrator.just_hit_tstop = false
    if integrator.opts.stop_at_next_tstop
      postamble!(integrator)
      return true
    end
  end
  false
end
function Base.iterate(integrator::DEIntegrator,state=0)
  done(integrator) && return nothing
  state += 1
  step!(integrator) # Iter updated in the step! header
  # Next is callbacks -> iterator  -> top
  return integrator,state
end

Base.eltype(::Type{T}) where {T<:DEIntegrator} = T
Base.IteratorSize(::Type{<:DEIntegrator}) = Base.SizeUnknown()

### Other Iterators

struct IntegratorTuples{I}
 integrator::I
end

function Base.iterate(tup::IntegratorTuples, state=0)
  done(tup.integrator) && return nothing
  step!(tup.integrator) # Iter updated in the step! header
  state += 1
  # Next is callbacks -> iterator  -> top
  return (tup.integrator.u,tup.integrator.t),state
end

Base.eltype(::Type{IntegratorTuples{I}}) where {U, T, I<:DEIntegrator{<:Any, <:Any, U, T}} = Tuple{U, T}
Base.IteratorSize(::Type{<:IntegratorTuples}) = Base.SizeUnknown()

RecursiveArrayTools.tuples(integrator::DEIntegrator) = IntegratorTuples(integrator)

"""
$(TYPEDEF)
"""
struct IntegratorIntervals{I}
 integrator::I
end

function Base.iterate(tup::IntegratorIntervals,state=0)
  done(tup.integrator) && return nothing
  state += 1
  step!(tup.integrator) # Iter updated in the step! header
  # Next is callbacks -> iterator  -> top
  return (tup.integrator.uprev,tup.integrator.tprev,tup.integrator.u,tup.integrator.t),state
end

Base.eltype(::Type{IntegratorIntervals{I}}) where {U, T, I<:DEIntegrator{<:Any, <:Any, U, T}} = Tuple{U, T, U, T}
Base.IteratorSize(::Type{<:IntegratorIntervals}) = Base.SizeUnknown()

intervals(integrator::DEIntegrator) = IntegratorIntervals(integrator)

struct TimeChoiceIterator{T,T2}
  integrator::T
  ts::T2
end

function Base.iterate(iter::TimeChoiceIterator,state=1)
  state > length(iter.ts) && return nothing
  t = iter.ts[state]
  integrator = iter.integrator
  if isinplace(integrator.sol.prob)
    tmp = first(get_tmp_cache(integrator))
    if t == integrator.t
      tmp .= integrator.u
    elseif t < integrator.t
      integrator(tmp,t)
    else
      step!(integrator,t-integrator.t)
      integrator(tmp,t)
    end
    return (tmp,t),state+1
  else
    if t == integrator.t
      tmp = integrator.u
    elseif t < integrator.t
      tmp = integrator(t)
    else
      step!(integrator,t-integrator.t)
      tmp = integrator(t)
    end
    return (tmp,t),state+1
  end
end

Base.length(iter::TimeChoiceIterator) = length(iter.ts)

@recipe function f(integrator::DEIntegrator;
                    denseplot=(integrator.opts.calck || typeof(integrator) <: AbstractSDEIntegrator)  && integrator.iter>0,
                    plotdensity =10,
                    plot_analytic=false,vars=nothing)

  int_vars = interpret_vars(vars,integrator.sol)

  if denseplot
    # Generate the points from the plot from dense function
    plott = collect(range(integrator.tprev;step=integrator.t,length=plotdensity))
    plot_timeseries = integrator(plott)
    if plot_analytic
      plot_analytic_timeseries = [integrator.sol.prob.f(Val{:analytic},t,integrator.sol.prob.u0) for t in plott]
    end
  end # if not denseplot, we'll just get the values right from the integrator.

  dims = length(int_vars[1])
  for var in int_vars
    @assert length(var) == dims
  end
  # Should check that all have the same dims!

  plot_vecs = []
  for i in 2:dims
    push!(plot_vecs,[])
  end

  labels = String[]# Array{String, 2}(1, length(int_vars)*(1+plot_analytic))
  for x in int_vars
    for j in 2:dims
      if denseplot
        push!(plot_vecs[j-1], u_n(plot_timeseries, x[j],integrator.sol,plott,plot_timeseries))
      else # just get values
        if x[j] == 0
          push!(plot_vecs[j-1], integrator.t)
        elseif x[j]==1 && !(typeof(integrator.u) <: AbstractArray)
          push!(plot_vecs[j-1], integrator.u)
        else
          push!(plot_vecs[j-1], integrator.u[x[j]])
        end
      end
    end
    add_labels!(labels,x,dims,integrator.sol)
  end

  if plot_analytic
    for x in int_vars
      for j in 1:dims
        if denseplot
          push!(plot_vecs[j], u_n(plot_timeseries, x[j],sol,plott,plot_timeseries))
        else # Just get values
          if x[j] == 0
            push!(plot_vecs[j], integrator.t)
          elseif x[j]==1 && !(typeof(integrator.u) <: AbstractArray)
            push!(plot_vecs[j], integrator.sol.prob.f(Val{:analytic},integrator.t,integrator.sol[1]))
          else
            push!(plot_vecs[j], integrator.sol.prob.f(Val{:analytic},integrator.t,integrator.sol[1])[x[j]])
          end
        end
      end
      add_labels!(labels,x,dims,integrator.sol)
    end
  end

  xflip --> integrator.tdir < 0

  if denseplot
    seriestype --> :path
  else
    seriestype --> :scatter
  end

  # Special case labels when vars = (:x,:y,:z) or (:x) or [:x,:y] ...
  if typeof(vars) <: Tuple && (typeof(vars[1]) == Symbol && typeof(vars[2]) == Symbol)
    xlabel --> vars[1]
    ylabel --> vars[2]
    if length(vars) > 2
      zlabel --> vars[3]
    end
  end
  if getindex.(int_vars,1) == zeros(length(int_vars)) || getindex.(int_vars,2) == zeros(length(int_vars))
    xlabel --> "t"
  end

  linewidth --> 3
  #xtickfont --> font(11)
  #ytickfont --> font(11)
  #legendfont --> font(11)
  #guidefont  --> font(11)
  label --> reshape(labels,1,length(labels))
  (plot_vecs...,)
end

function step!(integ::DEIntegrator, dt, stop_at_tdt = false)
    (dt * integ.tdir) < 0 * oneunit(dt) && error("Cannot step backward.")
    t = integ.t
    next_t = t+dt
    stop_at_tdt && add_tstop!(integ,next_t)
    while integ.t*integ.tdir < next_t*integ.tdir
        step!(integ)
        integ.sol.retcode in (:Default, :Success) || break
    end
end

has_destats(i::DEIntegrator) = false

"""
    is_integrator_adaptive(i::DEIntegrator)

Checks if the integrator is adaptive
"""
isadaptive(integrator::DEIntegrator) =
        isdefined(integrator.opts, :adaptive) ? integrator.opts.adaptive : false


text:
# function distoration_indices(choose_distortion::Vector{Symbol} = [:Y,:Z,:ψ,:θ,:α,:β])
#     all_distortion_order = [:Y,:Z,:ψ,:θ,:φ,:α,:β]
#     inds = [findfirst(all_distortion_order .== d) for d in choose_distortion]
#     return inds
# end

function camera_distoration(camera_reference::VideoCamera, camera::VideoCamera)

    keys = [:X,:Y,:Z,:ψ,:θ,:φ,:f]
    δ = [- (camera.xyz - camera_reference.xyz);
        camera.ψθφ - camera_reference.ψθφ;
        focal_length(camera) - focal_length(camera_reference)
    ]

    return Dict(keys .=> δ)
end

function rail_uvs_to_distortion(left_track_uvs::AbstractVector{S}, right_track_uvs::AbstractVector{S}, camera_reference::VideoCamera{T}, trackproperties::TrackProperties{T};
        iterations::Int = 1, choose_distortions::Vector{Symbol} = [:Y,:Z,:ψ,:θ,:α,:β]
    ) where {T<:AbstractFloat,S<:AbstractVector{T}}

    uLs = [uv[1] for uv in left_track_uvs];
    vLs = [uv[2] for uv in left_track_uvs];

    uRs = [uv[1] for uv in right_track_uvs];
    vRs = [uv[2] for uv in right_track_uvs];

    camera = camera_reference

    line_distortions = [:Y,:Z,:ψ,:θ,:φ,:ψT,:θT,:φT,:f]
    w = if length(intersect(line_distortions, choose_distortions)) <= 4
        zero(T)
    else
        eps(T)^T(2/3)
    end

    L = length(choose_distortions)
    δ = zeros(T,L)

    for j = 1:iterations
        ΔuL = left_track_Δu(camera, trackproperties, vLs; choose_distortions = choose_distortions);
        ΔuR = right_track_Δu(camera, trackproperties, vRs;  choose_distortions = choose_distortions);

        M = transpose(ΔuL) * ΔuL + transpose(ΔuR) * ΔuR + w .* I;

        uLs0 = left_track_image_u(camera, trackproperties, vLs);
        uRs0 = right_track_image_u(camera, trackproperties, vRs);
        duLs = uLs - uLs0;
        duRs = uRs - uRs0;

        A = transpose(ΔuL) * duLs + transpose(ΔuR) * duRs;
        δ = inv(M) * A
        #
        #
        # inv(transpose(ΔuL) * ΔuL) * transpose(ΔuL) * duLs
        #
        # inv(transpose(ΔuL) * ΔuL) * transpose(ΔuL) * (ΔuL * δs[i])
        #
        # inv(transpose(ΔuL) * ΔuL + transpose(ΔuR) * ΔuR)

        camera = VideoCamera(camera, Dict(choose_distortions .=> δ))
    end

    dict = camera_distoration(camera_reference, camera)
    for s in [:α,:β,:ψT,:θT,:φT]
        k = findfirst(choose_distortions .== s)
        if !isnothing(k)
            dict[s] = δ[k]
        end
    end

   return dict
end

"""
    Δu(vs::AbstractArray, Y::T, camera::VideoCamera)

Returns an array where Δu[i,j] corresponds to ∂_j u_i, and u_i = u(vs[i]).
"""
function Δu(camera::VideoCamera{T}, Y::T, vs::AbstractArray{T}; choose_distortions::Vector{Symbol} = [:Y,:Z,:ψ,:θ,:α,:β]) where T
    # get the basis functions corresponding to the choose distorations
    fs = [eval(Symbol("v_to_dud",s)) for s in choose_distortions]

    δus = map(fs) do f
        δu = f(camera, Y)
        δu.(vs)
    end

    return hcat(δus...)
end

left_track_Δu(camera::VideoCamera{T}, trackproperties::TrackProperties{T}, vs::AbstractArray{T}; kws...) where T = Δu(camera, - camera.xyz[2] - trackproperties.track_gauge / T(2), vs; kws...)

right_track_Δu(camera::VideoCamera{T}, trackproperties::TrackProperties{T}, vs::AbstractArray{T}; kws...) where T = Δu(camera, - camera.xyz[2] + trackproperties.track_gauge / T(2), vs; kws...)

function v_to_dudf(camera::VideoCamera{T}, Y::T; Z::T = - camera.xyz[3]) where T
    ψ, θ, φ = camera.ψθφ

    # does not depend on v
    δu(v::T) = (Y*sin(θ) + Z*cos(θ)*sin(φ)) / (Y*cos(θ)*sin(ψ) - Z*(cos(φ)*cos(ψ) + sin(θ)*sin(φ)*sin(ψ)))
    return δu
end

function v_to_dudψ(camera::VideoCamera{T}, Y::T; Z::T = - camera.xyz[3]) where T
    ψ, θ, φ = camera.ψθφ

    f = camera.opticalproperties.focal_length * camera.opticalproperties.pixelspermeter

    denom = (Z * cos(φ) * cos(ψ) + (-Y * cos(θ) + Z * sin(θ) * sin(φ)) * sin(ψ))^2

    dau = (Z^2 * cos(φ)^2 + (Y * cos(θ) -  Z * sin(θ) * sin(φ))^2) / denom
    dbu = - f * (Y * sin(θ) + Z * cos(θ) * sin(φ)) * (Y * cos(θ) * cos(ψ) - Z * cos(ψ) * sin(θ) * sin(φ) + Z * cos(φ) * sin(ψ)) / denom

    δu(v::T) = dau * v + dbu

    return δu
end

function v_to_dudθ(camera::VideoCamera{T}, Y::T; Z::T = - camera.xyz[3]) where T
    ψ, θ, φ = camera.ψθφ

    f = camera.opticalproperties.focal_length * camera.opticalproperties.pixelspermeter
    denom = (Z * cos(φ) * cos(ψ) + (-Y * cos(θ) + Z * sin(θ) * sin(φ)) * sin(ψ))^2

    dau = - Z * cos(φ) * (Y * sin(θ) + Z * cos(θ) * sin(φ)) / denom
    # dbu = (Z * cos(ψ) * (-2Y * cos(θ) * cos(φ) + Z * sin(θ) * sin(2φ)) + 2 * (Y^2 + Z^2 * sin(φ)^2) * sin(ψ)) / 2

    dbu = Z * cos(φ) * cos(ψ) * (-Y * cos(θ) + Z * sin(θ) * sin(φ)) + (Y^2 + Z^2 * sin(φ)^2) * sin(ψ)
    dbu = f * dbu / denom

    δu(v::T) = dau * v + dbu

    return δu
end

function v_to_dudφ(camera::VideoCamera{T}, Y::T; Z::T = - camera.xyz[3]) where T
    ψ, θ, φ = camera.ψθφ

    f = camera.opticalproperties.focal_length * camera.opticalproperties.pixelspermeter
    denom = (Z * cos(φ) * cos(ψ) + (-Y * cos(θ) + Z * sin(θ) * sin(φ)) * sin(ψ))^2

    dau = Z * (-Z * sin(θ) + Y * cos(θ) * sin(φ)) / denom
    dbu = f * Z * (-cos(ψ) * (Z * cos(θ) + Y * sin(θ) * sin(φ)) + Y * cos(φ) * sin(ψ)) / denom

    δu(v::T) = dau * v + dbu

    return δu
end

function v_to_dudY(camera::VideoCamera{T}, Y::T; Z::T = - camera.xyz[3]) where T
    ψ, θ, φ = camera.ψθφ

    f = camera.opticalproperties.focal_length * camera.opticalproperties.pixelspermeter
    denom = (Z * cos(φ) * cos(ψ) + (-Y * cos(θ) + Z * sin(θ) * sin(φ)) * sin(ψ))^2

    dau = Z * cos(θ) * cos(φ) / denom
    dbu = - f * Z * (cos(φ) * cos(ψ) * sin(θ) + sin(φ) * sin(ψ)) / denom

    δu(v::T) = dau * v + dbu

    return δu
end

function v_to_dudZ(camera::VideoCamera{T}, Y::T; Z::T = - camera.xyz[3]) where T
    ψ, θ, φ = camera.ψθφ

    f = camera.opticalproperties.focal_length * camera.opticalproperties.pixelspermeter

    denom = (Z * cos(φ) * cos(ψ) + (-Y * cos(θ) + Z * sin(θ) * sin(φ)) * sin(ψ))^2;

    dau = - (Y * cos(θ) * cos(φ)) / denom
    dbu = f * Y * (cos(φ) * cos(ψ) * sin(θ) + sin(φ) * sin(ψ)) / denom

    δu(v::T) = dau * v + dbu

    return δu
end

function v_to_X(camera::VideoCamera{T}, Y::T; Z::T = - camera.xyz[3]) where T
    ψ, θ, φ = camera.ψθφ

    f = camera.opticalproperties.focal_length * camera.opticalproperties.pixelspermeter

    top = [-f * cos(ψ) * (Z * cos(θ) + Y * sin(θ) * sin(φ)) + f * Y * cos(φ) * sin(ψ),
        -Z * sin(θ) + Y * cos(θ) * sin(φ)]

    bottom = [f * (cos(φ) * cos(ψ) * sin(θ) + sin(φ) * sin(ψ)),
    -cos(θ) * cos(φ)]

    X(v::T) = (top[1] + top[2] * v) / (bottom[1] + bottom[2] * v)

    return X
end

function v_to_dudβ(camera::VideoCamera{T}, Y::T; Z::T = - camera.xyz[3]) where T
    ψ, θ, φ = camera.ψθφ

    f = camera.opticalproperties.focal_length * camera.opticalproperties.pixelspermeter

    dudY = v_to_dudY(camera, Y; Z = Z)
    X = v_to_X(camera, Y; Z = Z)

    δu(v::T) = dudY(v) * X(v)^2

    return δu
end

function v_to_dudα(camera::VideoCamera{T}, Y::T; Z::T = - camera.xyz[3]) where T
    ψ, θ, φ = camera.ψθφ

    f = camera.opticalproperties.focal_length * camera.opticalproperties.pixelspermeter

    dudZ = v_to_dudZ(camera, Y; Z = Z)
    X = v_to_X(camera, Y;  Z = Z)

    δu(v::T) = dudZ(v) * X(v)^2

    return δu
end

function v_to_dudψT(camera::VideoCamera{T}, Y::T; kws...) where T
    ψ, θ, φ = camera.ψθφ

    dudψ = v_to_dudψ(camera, Y; kws...)
    dudθ = v_to_dudθ(camera, Y; kws...)
    dudφ = v_to_dudφ(camera, Y; kws...)

    δu(v::T) = dudψ(v) * cos(φ) / cos(θ) - dudθ(v) * sin(φ) + dudφ(v) * cos(φ) * tan(θ)

    return δu
end

function v_to_dudθT(camera::VideoCamera{T}, Y::T; kws...) where T
    ψ, θ, φ = camera.ψθφ

    dudψ = v_to_dudψ(camera, Y; kws...)
    dudθ = v_to_dudθ(camera, Y; kws...)
    dudφ = v_to_dudφ(camera, Y; kws...)

    δu(v::T) = dudψ(v) * sec(θ) * sin(φ) + dudθ(v) * cos(φ) + dudφ(v) * sin(φ) * tan(θ)

    return δu
end

v_to_dudφT(camera::VideoCamera{T}, Y::T; kws...) where T = v_to_dudφ(camera, Y; kws...)

ext: jl


text:
f(a, b) = sum(a' * b + a * b')

@rdcmp f a=rand(3,3) b=rand(3,3)

#-----------------

function main_1093()
    ex = :(sum(a' * b + a * b'))
    a = rand(3,3)
    b = rand(3,3)
    inputs = [:a => a, :b => b]
    ctx = Dict()
    dex = xdiff(ex; ctx=ctx, inputs...)
end


text:
using MPI
using StaticArrays
using CLIMA
using CLIMA.VariableTemplates
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.MPIStateArrays
using CLIMA.DGmethods
using CLIMA.DGmethods.NumericalFluxes
using Printf
using LinearAlgebra
using Logging
using GPUifyLoops

const ArrayType = CLIMA.array_type()

import CLIMA.DGmethods: BalanceLaw, vars_aux, vars_state, vars_gradient,
                        vars_diffusive, vars_integrals, integrate_aux!,
                        flux_nondiffusive!, flux_diffusive!, source!, wavespeed,
                        update_aux!, indefinite_stack_integral!,
                        reverse_indefinite_stack_integral!,  boundary_state!,
                        gradvariables!, init_aux!, init_state!,
                        init_ode_state, LocalGeometry

struct IntegralTestSphereModel{T} <: BalanceLaw
  Rinner::T
  Router::T
end

function update_aux!(dg::DGModel, m::IntegralTestSphereModel, Q::MPIStateArray, t::Real)
  indefinite_stack_integral!(dg, m, Q, dg.auxstate, t)
  reverse_indefinite_stack_integral!(dg, m, dg.auxstate, t)

  return true
end

vars_integrals(::IntegralTestSphereModel, T) = @vars(v::T)
vars_aux(m::IntegralTestSphereModel,T) = @vars(int::vars_integrals(m,T), rev_int::vars_integrals(m,T), r::T, a::T)

vars_state(::IntegralTestSphereModel, T) = @vars()
vars_diffusive(::IntegralTestSphereModel, T) = @vars()

flux_nondiffusive!(::IntegralTestSphereModel, _...) = nothing
flux_diffusive!(::IntegralTestSphereModel, _...) = nothing
source!(::IntegralTestSphereModel, _...) = nothing
boundary_state!(_, ::IntegralTestSphereModel, _...) = nothing
init_state!(::IntegralTestSphereModel, _...) = nothing
wavespeed(::IntegralTestSphereModel,_...) = 1

function init_aux!(m::IntegralTestSphereModel, aux::Vars, g::LocalGeometry)

  x,y,z = g.coord
  aux.r = hypot(x, y, z)
  θ = atan(y, x)
  ϕ = asin(z / aux.r)
  # Exact integral
  aux.a = 1 + cos(ϕ)^2 * sin(θ)^2 + sin(ϕ)^2
  aux.int.v = exp(-aux.a * aux.r^2) - exp(-aux.a * m.Rinner^2)
  aux.rev_int.v = exp(-aux.a * m.Router^2) - exp(-aux.a * aux.r^2)
end

@inline function integrate_aux!(m::IntegralTestSphereModel, integrand::Vars, state::Vars, aux::Vars)
  integrand.v = -2aux.r * aux.a * exp(-aux.a * aux.r^2)
end

if !@isdefined integration_testing
  const integration_testing =
    parse(Bool, lowercase(get(ENV,"JULIA_CLIMA_INTEGRATION_TESTING","false")))
end

using Test
function run(mpicomm, topl, ArrayType, N, FT, Rinner, Router)
  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                          meshwarp = Topologies.cubedshellwarp,
                                         )
  dg = DGModel(IntegralTestSphereModel(Rinner, Router),
               grid,
               Rusanov(),
               CentralNumericalFluxDiffusive(),
               CentralGradPenalty())

  Q = init_ode_state(dg, FT(0))
  dQdt = similar(Q)

  exact_aux = copy(dg.auxstate)
  dg(dQdt, Q, nothing, 0.0)
  euclidean_distance(exact_aux, dg.auxstate)
end

let
  CLIMA.init()
  mpicomm = MPI.COMM_WORLD
  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = ll == "DEBUG" ? Logging.Debug :
  ll == "WARN"  ? Logging.Warn  :
  ll == "ERROR" ? Logging.Error : Logging.Info
  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))

  polynomialorder = 4

  base_Nhorz = 4
  base_Nvert = 2
  Rinner = 1//2
  Router = 1

  polynomialorder = 4

  expected_result = [4.662884229467401e-7,
                     7.218989778540723e-9,
                     1.1258613174916711e-10,
                     1.7587739986848968e-12]
  lvls = integration_testing ? length(expected_result) : 1

  for FT in (Float64,) # Float32)
    err = zeros(FT, lvls)
    for l = 1:lvls
      @info (ArrayType, FT, "sphere", l)
      Nhorz = 2^(l-1) * base_Nhorz
      Nvert = 2^(l-1) * base_Nvert
      Rrange = grid1d(FT(Rinner), FT(Router); nelem=Nvert)
      topl = StackedCubedSphereTopology(mpicomm, Nhorz, Rrange)
      err[l] = run(mpicomm, topl, ArrayType, polynomialorder, FT,
                    FT(Rinner), FT(Router))
      @test expected_result[l] ≈ err[l] rtol=1e-3 atol=eps(FT)
    end
    if integration_testing
      @info begin
        msg = ""
        for l = 1:lvls-1
          rate = log2(err[l]) - log2(err[l+1])
          msg *= @sprintf("\n  rate for level %d = %e\n", l, rate)
        end
        msg
      end
    end
  end
end

nothing

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
# Module for simulating the kite control unit (KCU).

module KitePodModels

using Reexport
@reexport using KiteUtils
using Parameters, DocStringExtensions

export calc_alpha_depower, init_kcu, set_depower_steering, get_depower, get_steering
export KCU

if ! isfile("./data/settings.yaml")
    copy_settings()
end

"""
    mutable struct KCU

Stuct that is storing the settings and the state of the kite control unit.
"""
@with_kw mutable struct KCU
    set::KiteUtils.Settings = se()
    set_depower::Float64 =         0.0
    set_steering::Float64 =        0.0
    depower::Float64 =             0.0                         #    0 .. 1.0
    steering::Float64 =            0.0                         # -1.0 .. 1.0
end

"""
    KCU(set::KiteUtils.Settings)

Constructor for KCU struct.
"""
function KCU(set::KiteUtils.Settings)
    kcu = KCU()
    init_kcu(kcu, set)
    return kcu
end

"""
    init_kcu(kcu::KCU, set::KiteUtils.Settings)

Inititalze the model of the kite control unit (KCU).
The actual and the set values of depower are initialized to set.depower_offset * 0.01,
the actual and the set values of steering to zero.
"""
function init_kcu(kcu::KCU, set::KiteUtils.Settings)
    kcu.set = set
    kcu.set_depower =         set.depower_offset * 0.01
    kcu.set_steering =        0.0
    kcu.depower =             set.depower_offset * 0.01   #    0 .. 1.0
    kcu.steering =            0.0                         # -1.0 .. 1.0
    nothing
end

# Calculate the length increase of the depower line [m] as function of the relative depower
# setting [0..1].
function calc_delta_l(kcu::KCU, rel_depower)
    u = kcu.set.depower_drum_diameter
    l_ro = 0.0
    rotations = (rel_depower - 0.01 * kcu.set.depower_offset) * 10.0 * 11.0 / 3.0 * (3918.8 - 230.8) / 4096.
    while rotations > 0.0
         l_ro += u * π
         rotations -= 1.0
         u -= kcu.set.tape_thickness
    end
    if rotations < 0.0
         l_ro += (-(u + kcu.set.tape_thickness) * rotations + u * (rotations + 1.0)) * π * rotations
    end
    return l_ro
end

"""
    calc_alpha_depower(kcu::KCU, rel_depower)

Calculate the change of the angle between the kite and the last tether segment [rad] as function of the
actual rel_depower value.

Returns `nothing` in case of error.
"""
function calc_alpha_depower(kcu::KCU, rel_depower)
    a   =  kcu.set.power2steer_dist
    b_0 = kcu.set.h_bridle + 0.5 * kcu.set.height_k
    b   = b_0 + 0.5 * calc_delta_l(kcu, rel_depower) # factor 0.5 due to the pulleys

    c = sqrt(a * a + b_0 * b_0)
    # print 'a, b, c:', a, b, c, rel_depower
    if c >= a + b
         println("ERROR in calc_alpha_depower!")
         return nothing
    else
        tmp = 1/(2*a*b)*(a*a+b*b-c*c)
        if tmp > 1.0
            println("-->>> WARNING: tmp > 1.0: $tmp")
            tmp = 1.0
        elseif tmp < -1.0
            println("-->>> WARNING: tmp < 1.0: $tmp")
            tmp = -1.0
        end
        return pi/2.0 - acos(tmp)
    end
end

"""
    set_depower_steering(kcu::KCU, depower, steering)

Set the values of depower and steering. The value for depower must be between 0.0 and 1.0,
the value for steering between -1.0 and +1.0 .
"""
function set_depower_steering(kcu::KCU, depower, steering)
    kcu.set_depower  = depower
    kcu.set_steering = steering
    nothing
end

"""
    get_depower(kcu::KCU)

Read the current depower value. Result will be between 0.0 and 1.0.
"""
function get_depower(kcu::KCU);  return kcu.depower;  end

"""
    get_steering(kcu::KCU)

Read the current depower value. Result will be between -1.0 and 1.0.
"""
function get_steering(kcu::KCU); return kcu.steering; end

"""
    on_timer(kcu::KCU, dt = 0.0)

Must be called at each clock tick. Parameter: Δt in seconds
Updates the current values of steering and depower depending
on the set values and the last value.
If dt == 0.0, then it will be set to 1.0 / kcu.set.sample_freq .
"""
function on_timer(kcu::KCU, dt = 0.0)
    if dt == 0.0
        dt = 1.0 / kcu.set.sample_freq
    end
    # calculate the depower motor velocity
    vel_depower = (kcu.set_depower - kcu.depower) * kcu.set.depower_gain
    # println("vel_depower: $(vel_depower)")
    if vel_depower > kcu.set.v_depower
        vel_depower = kcu.set.v_depower
    elseif vel_depower < -kcu.set.v_depower
        vel_depower = -kcu.set.v_depower
    end
    # update the position
    kcu.depower += vel_depower * dt
    # calculate the steering motor velocity
    vel_steering = (kcu.set_steering - kcu.steering) * kcu.set.steering_gain
    if vel_steering > kcu.set.v_steering
        vel_steering = kcu.set.v_steering
    elseif vel_steering < -kcu.set.v_steering
        vel_steering = -kcu.set.v_steering
    end
    kcu.steering += vel_steering * dt
    nothing
end

end


text:
using ArrayInterface: dimnames

@testset "dimensions" begin

###
### define wrapper with dimnames
###

struct NamedDimsWrapper{L,T,N,P<:AbstractArray{T,N}} <: ArrayInterface.AbstractArray2{T,N}
    parent::P
    NamedDimsWrapper{L}(p) where {L} = new{L,eltype(p),ndims(p),typeof(p)}(p)
end
ArrayInterface.parent_type(::Type{T}) where {P,T<:NamedDimsWrapper{<:Any,<:Any,<:Any,P}} = P
ArrayInterface.has_dimnames(::Type{T}) where {T<:NamedDimsWrapper} = true
ArrayInterface.dimnames(::Type{T}) where {L,T<:NamedDimsWrapper{L}} = static(Val(L))
function ArrayInterface.dimnames(::Type{T}, dim) where {L,T<:NamedDimsWrapper{L}}
    if ndims(T) < dim
        return static(:_)
    else
        return static(L[dim])
    end
end
ArrayInterface.has_dimnames(::Type{T}) where {T<:NamedDimsWrapper} = true
Base.parent(x::NamedDimsWrapper) = x.parent

@testset "dimension permutations" begin
    a = ones(2, 2, 2)
    perm = PermutedDimsArray(a, (3, 1, 2))
    mview = view(perm, :, 1, :)
    madj = mview'
    vview = view(madj, 1, :)
    vadj = vview'

    @test @inferred(ArrayInterface.to_parent_dims(typeof(a))) == (1, 2, 3)
    @test @inferred(ArrayInterface.to_parent_dims(typeof(perm))) == (3, 1, 2)
    @test @inferred(ArrayInterface.to_parent_dims(typeof(mview))) == (1, 3)
    @test @inferred(ArrayInterface.to_parent_dims(typeof(madj))) == (2, 1)
    @test @inferred(ArrayInterface.to_parent_dims(typeof(vview))) == (2,)
    @test @inferred(ArrayInterface.to_parent_dims(typeof(vadj))) == (2, 1)
    @test @inferred(ArrayInterface.to_parent_dims(typeof(vadj), static(1))) == 2
    @test @inferred(ArrayInterface.to_parent_dims(typeof(vadj), 1)) == 2
    @test @inferred(ArrayInterface.to_parent_dims(typeof(vadj), static(3))) == 2
    @test @inferred(ArrayInterface.to_parent_dims(typeof(vadj), 3)) == 2

    @test @inferred(ArrayInterface.from_parent_dims(typeof(a))) == (1, 2, 3)
    @test @inferred(ArrayInterface.from_parent_dims(typeof(perm))) == (2, 3, 1)
    @test @inferred(ArrayInterface.from_parent_dims(typeof(mview))) == (1, 0, 2)
    @test @inferred(ArrayInterface.from_parent_dims(typeof(madj))) == (2, 1)
    @test @inferred(ArrayInterface.from_parent_dims(typeof(vview))) == (0, 1)
    @test @inferred(ArrayInterface.from_parent_dims(typeof(vadj))) == (2,)
    @test @inferred(ArrayInterface.from_parent_dims(typeof(vadj), static(1))) == 2
    @test @inferred(ArrayInterface.from_parent_dims(typeof(vadj), 1)) == 2

    @test_throws DimensionMismatch ArrayInterface.to_parent_dims(typeof(vadj), 0)
    @test_throws DimensionMismatch ArrayInterface.to_parent_dims(typeof(vadj), static(0))

    @test_throws DimensionMismatch ArrayInterface.from_parent_dims(typeof(vadj), 0)
    @test_throws DimensionMismatch ArrayInterface.from_parent_dims(typeof(vadj), static(0))

    if VERSION ≥ v"1.6.0-DEV.1581"
        colormat = reinterpret(reshape, Float64, [(R = rand(), G = rand(), B = rand()) for i ∈ 1:100])
        @test @inferred(ArrayInterface.from_parent_dims(typeof(colormat))) === (static(2),)
        @test @inferred(ArrayInterface.to_parent_dims(typeof(colormat))) === (static(0), static(1),)

        Rr = reinterpret(reshape, Int32, ones(4))
        @test @inferred(ArrayInterface.from_parent_dims(typeof(Rr))) === (static(2),)
        @test @inferred(ArrayInterface.to_parent_dims(typeof(Rr))) === (static(0), static(1),)

        Rr = reinterpret(reshape, Int64, ones(4))
        @test @inferred(ArrayInterface.from_parent_dims(typeof(Rr))) === (static(1),)
        @test @inferred(ArrayInterface.to_parent_dims(typeof(Rr))) === (static(1),)

        Sr = reinterpret(reshape, Complex{Int64}, zeros(2, 3, 4))
        @test @inferred(ArrayInterface.from_parent_dims(typeof(Sr))) === (static(0), static(1), static(2))
        @test @inferred(ArrayInterface.to_parent_dims(typeof(Sr))) === (static(2), static(3))
    end
end

@testset "order_named_inds" begin
    n1 = (static(:x),)
    n2 = (n1..., static(:y))
    n3 = (n2..., static(:z))
    @test @inferred(ArrayInterface.order_named_inds(n1, NamedTuple{(),Tuple{}}(()) )) == ()
    @test @inferred(ArrayInterface.order_named_inds(n1, (x=2,))) == (2,)
    @test @inferred(ArrayInterface.order_named_inds(n2, (x=2,))) == (2, :)
    @test @inferred(ArrayInterface.order_named_inds(n2, (y=2,))) == (:, 2)
    @test @inferred(ArrayInterface.order_named_inds(n2, (y=20, x=30))) == (30, 20)
    @test @inferred(ArrayInterface.order_named_inds(n2, (x=30, y=20))) == (30, 20)
    @test @inferred(ArrayInterface.order_named_inds(n3, (x=30, y=20))) == (30, 20, :)

    @test_throws ErrorException ArrayInterface.order_named_inds(n2, (x=30, y=20, z=40))
end

val_has_dimnames(x) = Val(ArrayInterface.has_dimnames(x))

@testset "dimnames" begin
    d = (static(:x), static(:y))
    x = NamedDimsWrapper{d}(ones(2,2));
    y = NamedDimsWrapper{(:x,)}(ones(2));
    dnums = ntuple(+, length(d))
    @test @inferred(val_has_dimnames(x)) === Val(true)
    @test @inferred(ArrayInterface.has_dimnames(ones(2,2))) === false
    @test @inferred(ArrayInterface.has_dimnames(Array{Int,2})) === false
    @test @inferred(val_has_dimnames(typeof(x))) === Val(true)
    @test @inferred(val_has_dimnames(typeof(view(x, :, 1, :)))) === Val(true)
    @test @inferred(dimnames(x)) === d
    @test @inferred(dimnames(parent(x))) === (static(:_), static(:_))
    @test @inferred(dimnames(x')) === reverse(d)
    @test @inferred(dimnames(y')) === (static(:_), static(:x))
    @test @inferred(dimnames(PermutedDimsArray(x, (2, 1)))) === reverse(d)
    @test @inferred(dimnames(PermutedDimsArray(x', (2, 1)))) === d
    @test @inferred(dimnames(view(x, :, 1))) === (static(:x),)
    @test @inferred(dimnames(view(x, :, :, :))) === (static(:x),static(:y), static(:_))
    @test @inferred(dimnames(view(x, :, 1, :))) === (static(:x), static(:_))
    @test @inferred(dimnames(x, ArrayInterface.One())) === static(:x)
    @test @inferred(dimnames(parent(x), ArrayInterface.One())) === static(:_)
end

@testset "to_dims" begin
    x = NamedDimsWrapper{(:x, :y)}(ones(2,2));
    y = NamedDimsWrapper{(:x, :y, :a, :b, :c, :d)}(ones(6));

    @test @inferred(ArrayInterface.to_dims(x, :)) == Colon()
    @test @inferred(ArrayInterface.to_dims(x, 1)) == 1
    @testset "small case" begin
        @test @inferred(ArrayInterface.to_dims(x, (:x, :y))) == (1, 2)
        @test @inferred(ArrayInterface.to_dims(x, (:y, :x))) == (2, 1)
        @test @inferred(ArrayInterface.to_dims(x, :x)) == 1
        @test @inferred(ArrayInterface.to_dims(x, :y)) == 2
        @test_throws DimensionMismatch ArrayInterface.to_dims(x, static(:z))  # not found
        @test_throws DimensionMismatch ArrayInterface.to_dims(x, :z)  # not found
    end

    @testset "large case" begin
        @test @inferred(ArrayInterface.to_dims(y, :x)) == 1
        @test @inferred(ArrayInterface.to_dims(y, :a)) == 3
        @test @inferred(ArrayInterface.to_dims(y, :d)) == 6
        @test_throws DimensionMismatch ArrayInterface.to_dims(y, :z) # not found
    end
end

@testset "methods accepting dimnames" begin
    d = (static(:x), static(:y))
    x = NamedDimsWrapper{d}(ones(2,2));
    y = NamedDimsWrapper{(:x,)}(ones(2));
    @test @inferred(size(x, first(d))) == size(parent(x), 1)
    @test @inferred(ArrayInterface.size(y')) == (1, size(parent(x), 1))
    @test @inferred(axes(x, first(d))) == axes(parent(x), 1)
    @test strides(x, :x) == ArrayInterface.strides(parent(x))[1]
    @test @inferred(ArrayInterface.axes_types(x, static(:x))) <: Base.OneTo{Int}
    @test ArrayInterface.axes_types(x, :x) <: Base.OneTo{Int}
    @test @inferred(ArrayInterface.axes_types(LinearIndices{2,NTuple{2,Base.OneTo{Int}}})) <: NTuple{2,Base.OneTo{Int}}
    CI = CartesianIndices{2,Tuple{Base.OneTo{Int},UnitRange{Int}}}
    @test @inferred(ArrayInterface.axes_types(CI, static(1))) <: Base.OneTo{Int}

    x[x = 1] = [2, 3]
    @test @inferred(getindex(x, x = 1)) == [2, 3]
end

end


text:
function simplified_matvec(Ξ::AbstractArray{T, 2}, basis) where T <: Real
    eqs = Operation[]
    for i=1:size(Ξ, 2)
        eq = nothing
        for j = 1:size(Ξ, 1)
            if !iszero(Ξ[j,i])
                if eq === nothing
                    eq = basis[j]*Ξ[j,i]
                else
                    eq += basis[j]*Ξ[j,i]
                end
            end
        end
        if eq != nothing
            push!(eqs, eq)
        end
    end
    eqs
end

function simplified_matvec(Ξ::AbstractArray{T,1}, basis) where T <: Real
    eq = nothing
    @inbounds for i in 1:size(Ξ, 1)
        if !iszero(Ξ[i])
            if eq === nothing
                eq = basis[i]*Ξ[i]
            else
                eq += basis[i]*Ξ[i]
            end
        end

    end
    eq
end

function normalize_theta!(scales::AbstractArray, θ::AbstractArray)
    @assert length(scales) == size(θ, 1)
    @inbounds for (i, ti) in enumerate(eachrow(θ))
        scales[i] = norm(ti, 2)
        normalize!(ti, 2)
    end
    return
end

function rescale_xi!(scales::AbstractArray, Ξ::AbstractArray)
    @inbounds for (si, ti) in zip(scales, eachrow(Ξ))
        ti .= ti / si
    end
    return
end

function sparse_regression(X::AbstractArray, Ẋ::AbstractArray, Ψ::Basis, p::AbstractArray , maxiter::Int64 , opt::T, denoise::Bool, normalize::Bool ) where T <: Optimise.AbstractOptimiser
    @assert size(X)[end] == size(Ẋ)[end]
    nx, nm = size(X)
    ny, nm = size(Ẋ)

    Ξ = zeros(eltype(X), length(Ψ), ny)
    scales = ones(eltype(X), length(Ψ))
    θ = Ψ(X, p = p)

    denoise ? optimal_shrinkage!(θ') : nothing
    normalize ? normalize_theta!(scales, θ) : nothing

    Optimise.init!(Ξ, opt, θ', Ẋ')
    Optimise.fit!(Ξ, θ', Ẋ', opt, maxiter = maxiter)

    normalize ? rescale_xi!(scales, Ξ) : nothing

    return Ξ
end

function sparse_regression!(Ξ::AbstractArray, X::AbstractArray, Ẋ::AbstractArray, Ψ::Basis, p::AbstractArray , maxiter::Int64 , opt::T, denoise::Bool, normalize::Bool ) where T <: Optimise.AbstractOptimiser
    @assert size(X)[end] == size(Ẋ)[end]
    nx, nm = size(X)
    ny, nm = size(Ẋ)
    @assert size(Ξ) == (length(Ψ), ny)

    scales = ones(eltype(X), length(Ψ))
    θ = Ψ(X, p = p)

    denoise ? optimal_shrinkage!(θ') : nothing
    normalize ? normalize_theta!(scales, θ) : nothing

    Optimise.init!(Ξ, opt, θ', Ẋ')
    Optimise.fit!(Ξ, θ', Ẋ', opt, maxiter = maxiter)

    normalize ? rescale_xi!(scales, Ξ) : nothing

    return
end

# For pareto
function sparse_regression!(Ξ::AbstractArray, θ::AbstractArray, Ẋ::AbstractArray, maxiter::Int64 , opt::T, denoise::Bool, normalize::Bool) where T <: Optimise.AbstractOptimiser

    scales = ones(eltype(Ξ), size(θ, 1))

    denoise ? optimal_shrinkage!(θ') : nothing
    normalize ? normalize_theta!(scales, θ) : nothing

    Optimise.init!(Ξ, opt, θ', Ẋ')'
    Optimise.fit!(Ξ, θ', Ẋ', opt, maxiter = maxiter)

    normalize ? rescale_xi!(scales, Ξ) : nothing
    normalize ? rescale_xi!(scales, θ) : nothing

    return
end

# One Variable on multiple derivatives
function SInDy(X::AbstractArray{S, 1}, Ẋ::AbstractArray, Ψ::Basis; kwargs...) where S <: Number
    return SInDy(X', Ẋ, Ψ; kwargs...)
end

# Multiple on one
function SInDy(X::AbstractArray{S, 2}, Ẋ::AbstractArray{S, 1}, Ψ::Basis; kwargs...) where S <: Number
    return SInDy(X, Ẋ', Ψ; kwargs...)
end

# General
function SInDy(X::AbstractArray{S, 2}, Ẋ::AbstractArray{S, 2}, Ψ::Basis; p::AbstractArray = [], maxiter::Int64 = 10, opt::T = Optimise.STRRidge(), denoise::Bool = false, normalize::Bool = true) where {T <: Optimise.AbstractOptimiser, S <: Number}
    Ξ = sparse_regression(X, Ẋ, Ψ, p, maxiter, opt, denoise, normalize)
    Basis(simplified_matvec(Ξ, Ψ.basis), variables(Ψ), parameters = p)
end

function SInDy(X::AbstractArray{S, 1}, Ẋ::AbstractArray, Ψ::Basis, thresholds::AbstractArray; kwargs...) where S <: Number
    return SInDy(X', Ẋ, Ψ, thresholds; kwargs...)
end

function SInDy(X::AbstractArray{S, 2}, Ẋ::AbstractArray{S, 1}, Ψ::Basis, thresholds::AbstractArray; kwargs...) where S <: Number
    return SInDy(X, Ẋ', Ψ, thresholds; kwargs...)
end

function SInDy(X::AbstractArray{S, 2}, Ẋ::AbstractArray{S, 2}, Ψ::Basis, thresholds::AbstractArray ; p::AbstractArray = [], maxiter::Int64 = 10, opt::T = Optimise.STRRidge(),denoise::Bool = false, normalize::Bool = true) where {T <: Optimise.AbstractOptimiser, S <: Number}
    @assert size(X)[end] == size(Ẋ)[end]
    nx, nm = size(X)
    ny, nm = size(Ẋ)

    θ = Ψ(X, p = p)

    ξ = zeros(eltype(X), length(Ψ), ny)
    Ξ_opt = zeros(eltype(X), length(Ψ), ny)
    Ξ = zeros(eltype(X), length(thresholds), ny, length(Ψ))
    x = zeros(eltype(X), length(thresholds), ny, 2)
    pareto = zeros(eltype(X),  ny, length(thresholds))
    scales = ones(eltype(X), length(Ψ))

    denoise ? optimal_shrinkage!(θ') : nothing
    normalize ? normalize_theta!(scales, θ) : nothing

    @inbounds for (j, threshold) in enumerate(thresholds)
        set_threshold!(opt, threshold)

        sparse_regression!(ξ, θ, Ẋ, maxiter, opt, false, false)

        [x[j, i, :] = [norm(xi, 0); norm(view(Ẋ , i, :) - θ'*xi, 2)] for (i, xi) in enumerate(eachcol(ξ))]

        normalize ? rescale_xi!(scales, ξ) : nothing
        Ξ[j, :, :] = ξ[:, :]'
    end

    # Create the evaluation
    @inbounds for i in 1:ny
        x[:, i, 2] .= x[:, i, 2]./maximum(x[:, i, 2])
        pareto[i, :] = [norm(x[j, i, :], 2) for j in 1:length(thresholds)]
        _, indx = findmin(pareto[i, :])
        Ξ_opt[:, i] = Ξ[indx, i, :]
    end

    return Basis(simplified_matvec(Ξ_opt, Ψ.basis), variables(Ψ), parameters = p)
end


text:
export RESP_wont_read

function resp_unit_split(buf::Array{UInt8,1}, L::Int64)
  i = 1
  while i+2 < L
    if buf[i] == 0x20 && buf[i+1] == 0x2d && buf[i+2] == 0x20
      return String(buf[1:i-1])
    end
    i = i+1
  end
  return String(buf[1:L])
end

# *** FIX ME ==> version in dataless.jl is better
function parse_resp_date!(buf::Array{UInt8,1}, L::Int64, T::Array{Int16,1})
  is_u8_digit(buf[1]) || return typemax(Int64)
  fill!(T, zero(Int16))
  i = 0
  j = 1 # counter to start position in buf
  k = 1 # counter to y,j,h,m,s
  while i ≤ L
    i += 1
    if i > L
      T[k] = buf_to_int(buf, i-1, j)
      T[k] = parse(Int16, String(buf[j:i-1]))
      break
    elseif buf[i] in (0x2c, 0x2e, 0x3a)
      T[k] = buf_to_int(buf, i-1, j)
      k += 1
      k > 6 && break
      j = i+1
    end
  end
  return mktime(T)
end

function close_resp_channel!(S::SeisData, C::SeisChannel, ts::Int64, te::Int64, fname::String, ns::Int64)
  L = length(C.resp.fs)
  for f in fieldnames(MultiStageResp)
    deleteat!(getfield(C.resp, f), (ns+1):L)
  end

  i = findid(S, C.id)
  if i != 0
    t0 = isempty(S.t[i]) ? ts : S.t[i][1,2]
    if ts ≤ t0 < te
      S.resp[i] = C.resp
      S.gain[i] = C.gain
      note!(S, i, string("seed_resp, ", fname, ", overwrote :gain, :resp"))
    end
  else
    # ObsPy does this but I don't think it's safe; stage fs is not stored uniformly.
    if length(C.resp.fs) > 0
      fac = C.resp.fac[end]
      C.fs = C.resp.fs[end] * (fac > 0 ? 1/fac : 1.0)
    end
    push!(S, C)
  end

  # fill!(delay, zero(Float64))
  # fill!(fg, zero(Float64))
  # fill!(fs_in, zero(Float64))
  # fill!(gain, zero(Float64))
  # fill!(corr, zero(Float64))
  # fill!(fac, zero(Int64))
  # fill!(os, zero(Int64))

  return nothing
end

function add_stage!(C::SeisChannel, n::Int64, tfc::UInt8)
  if C.resp.stage[n] == nothing
    C.resp.stage[n] =
      if tfc == 0x41
        PZResp64()
      elseif tfc == 0x44
        CoeffResp()
      else
        GenResp()
      end
    return nothing
  end
end

function read_seed_resp!(S::GphysData, fpat::String;
  units::Bool=false)

  buf   = BUF.buf
  C     = SeisChannel(resp = MultiStageResp(12))
  R     = C.resp
  blk   = zero(UInt64)

  # Containers for parts of the response
  # nstg    = 12
  # delay = zeros(Float64, nstg)        # Estimated delay (seconds) (047F08, 057F07)
  # corr  = zeros(Float64, nstg)        # Correction applied (seconds) (047F09, 057F08)
  # fg    = zeros(Float64, nstg)        # Frequency of gain (058F05)
  # fs_in = zeros(Float64, nstg)        # Input sample rate (047F05, 057F04)
  # gain  = zeros(Float64, nstg)        # Gain (058F04)
  # fac   = zeros(Int64,   nstg)        # Decimation factor (047F06, 057F05)
  # os    = zeros(Int64,   nstg)        # Decimation offset (047F07, 057F06)
  seq_n = zero(Int64)                 # Stage sequence number (053F04, 054F04, 057F03, 058F03, 060F04, 061F03, 062F04)
  seq_n_old = zero(Int64)             # Last sequence number
  d_arr = zeros(Int16, 6)             # Date
  nmax = zero(Int64)

  # Array containers
  X = Float64[]
  D = Float64[]
  Z = ComplexF64[]
  P = ComplexF64[]

  # Unit strings containers
  units_out = ""            # Response out units lookup (041F07, 043F07, 044F07, 053F06, 054F06, 061F07, 062F06)

  if safe_isfile(fpat)
    files = [fpat]
  else
    files = ls(fpat)
  end
  for (nf, file) in enumerate(files)

    # Counters, etc.
    c = 0x00
    i = 1
    read_state = 0x00
    tfc = 0x00
    ts = zero(Int64)          # Start date (052F22)
    te = zero(Int64)          # End date (052F23)
    ND = zero(Int64)          # Number of denominators (044F11, 054F10)
    NN = zero(Int64)          # Number of coefficients (041F08, 044F08, 054F07, 061F08, 062F14)
    NZ = zero(Int64)          # Number of zeroes (043F10, 053F09)
    NP = zero(Int64)          # Number of poles (043F15, 053F14)
    a0 = one(Float64)         # A0 normalization factor (043F08, 048F05, 053F07)
    f0 = one(Float64)         # Normalization frequency (043F09, 048F06, 053F08)
    id = UInt8[]
    units_in = ""             # Response in units lookup (041F06, 043F06, 044F06, 053F05, 054F05, 061F06, 062F05)

    io = open(file, "r")
    chip = false
    while !eof(io)
      c = read(io, UInt8)

      # comment line ------------------------------
      if c == 0x23
        while c != 0x0a
          c = read(io, UInt8)
        end
      end

      # start of blockette line -------------------
      if c == 0x42
        #= format: %c%3i%c%2i, 'B', blockette_#, 'F', field_#

        Here, I map the blockette and field numbers (along with 'F') to a UInt64.
        This mapping should be unique; for any blockette blk, I expect that
        String(reinterpret(UInt8, [blk]))[1:6] = bbbFff, where bbb is
        blockette # and ff is field #
        =#
        blk  = UInt64(read(io, UInt8))
        blk |= UInt64(read(io, UInt8)) << 8
        blk |= UInt64(read(io, UInt8)) << 16
        blk |= UInt64(read(io, UInt8)) << 24
        blk |= UInt64(read(io, UInt8)) << 32
        blk |= UInt64(read(io, UInt8)) << 40
        read_state = 0x01

        # Coefficient parsing ==================================================
        if blk in (0x0000393046313630, 0x0000393046313430)
          # Numerator coefficients (061F09, 041F09)
          X = Array{Float64,1}(undef, NN)
          for n = 1:NN
            # To first non-whitespace
            if n > 1
              skip(io, 7)
            end
            c = read(io, UInt8)
            while c == 0x20
              c = read(io, UInt8)
            end
            j = get_coeff_n(io, c, buf)
            c = skip_whitespace(io, c)
            # To next newline
            k = 1
            while c != 0x0a
              buf[k] = c
              k += 1
              c = read(io, UInt8)
            end
            store_dbl!(X, buf, k, j)
          end

        elseif blk in (0x0000313146333430, 0x0000303146333530, 0x0000363146333430, 0x0000353146333530)
          # Complex zeros (043F11, 053F10) / Complex poles (043F16, 053F15)

          if blk in (0x0000313146333430, 0x0000303146333530)
            N = NZ
            Z = Array{ComplexF64,1}(undef, N)
            is_p = false
          else
            N = NP
            P = Array{ComplexF64,1}(undef, N)
            is_p = true
          end

          for n = 1:N
            # coefficient number
            if n > 1
              skip(io, 10)
            else
              skip(io, 3)
            end
            c = skip_whitespace(io, c)
            j = get_coeff_n(io, c, buf)

            # real part
            c = skip_whitespace(io, c)
            k = 1
            while c != 0x20
              buf[k] = c
              k += 1
              c = read(io, UInt8)
            end
            xr = buf_to_double(buf, k)

            # imaginary part
            c = skip_whitespace(io, c)
            k = 1
            while c != 0x20
              buf[k] = c
              k += 1
              c = read(io, UInt8)
            end
            xi = buf_to_double(buf, k)

            # store
            if is_p
              P[j] = complex(xr, xi)
            else
              Z[j] = complex(xr, xi)
            end

            # rest irrelevant
            c = to_newline(io, c)
          end

        elseif blk in (0x0000383046343530, 0x0000313146343530, 0x0000353146323630)
          # Numerator coefficients (054F08) / Denominator coefficients (054F11) / Polynomial coefficients (B062F15)
          if blk == 0x0000383046343530
            N = NN
            X = Array{Float64,1}(undef, N)
            is_n = true
          else
            N = ND
            D = Array{Float64,1}(undef, N)
            is_n = false
          end
          for n = 1:N
            # To first non-whitespace
            if n > 1
              skip(io, 10)
            else
              skip(io, 3)
            end
            c = skip_whitespace(io, c)
            j = get_coeff_n(io, c, buf)
            c = skip_whitespace(io, c)

            # To next whitesoace
            k = 1
            while c != 0x20
              buf[k] = c
              k += 1
              c = read(io, UInt8)
            end

            # Store
            if is_n
              store_dbl!(X, buf, k, j)
            else
              store_dbl!(D, buf, k, j)
            end
            c = to_newline(io, c)
          end

        else
          # To separator
          while c != 0x3a
            c = read(io, UInt8)
          end

          # To first non-whitespace
          c = skip_whitespace(io, c)
        end
      end

      # header info --------------------------------
      if read_state == 0x01
        i = 1
        while c != 0x0a
          buf[i] = c
          i += 1
          eof(io) && break
          c = read(io, UInt8)
        end
        i -= 1

        # header parsing =======================================================
        if blk in (0x0000383046333430, 0x0000373046333530)
          # A0 normalization factor (043F08, 053F07)
          a0 = buf_to_double(buf, i)

        elseif blk == 0x0000343046323530
          # Channel (052F04)
          ni = min(i,3)
          L = length(id)
          append!(id, zeros(UInt8, ni+1))
          id[L+1] = 0x2e
          unsafe_copyto!(id, L+2, buf, 1, ni)
          setfield!(C, :id, String(id))
          # println("Starting channel ", C.id)
          chip = true

        elseif blk in (0x0000383046373530, 0x0000393046373430)
          # Correction applied (seconds) (057F08) / Response correction (047F09)
          store_dbl!(R.corr, buf, i, seq_n)

        elseif blk in (0x0000353046373530, 0x0000363046373430)
          # Decimation factor (057F05) / Response decimation factor (047F06)
          store_int!(R.fac, buf, i, seq_n)

        elseif blk in (0x0000363046373530, 0x0000373046373430)
          # Decimation offset (057F06) /  Response decimation offset (047F07)
          store_int!(R.os, buf, i, seq_n)

        elseif blk == 0x0000333246323530
          # End date (052F23)
          te = parse_resp_date!(buf, i, d_arr)

        elseif blk in (0x0000373046373530, 0x0000383046373430)
          # Estimated delay (seconds) (057F07) /  Response delay (047F08)
          store_dbl!(R.delay, buf, i, seq_n)

        elseif blk in (0x0000353046383530, 0x0000363046383430)
          # Frequency of gain (058F05) / Frequency of sensitivity (048F06)
          # println("blk = ", String(reinterpret(UInt8, [blk])), ", buf = ", String(buf[1:i]))
          store_dbl!(R.fg, buf, i, seq_n)

        elseif blk in (0x0000343046383530, 0x0000353046383430)
          # Gain (058F04) / Sensitivity (048F05)
          if seq_n > 0
            store_dbl!(R.gain, buf, i, seq_n)
          else
            C.gain = buf_to_double(buf, i)
          end

        elseif blk in (0x0000343046373530, 0x0000353046373430)
          # Input sample rate (057F04) / Response input sample rate (047F05)
          store_dbl!(R.fs, buf, i, seq_n)

        elseif blk == 0x0000333046323530
          # Location (052F03)
          push!(id, 0x2e)
          for j = 1:min(i,2)
            if buf[j] != 0x3f
              push!(id, buf[j])
            end
          end

        elseif blk == 0x0000363146303530
          # Network (050F16)
          ni = min(i,2)
          prepend!(id, zeros(UInt8, ni))
          copyto!(id, 1, buf, 1, ni)

        elseif blk in (0x0000393046333430, 0x0000383046333530)
          # Normalization frequency (043F09, 053F08)
          f0 = buf_to_double(buf, i)

        elseif blk in (0x0000313146343430, 0x0000303146343530)
          # Number of denominators (044F11, 054F10)
          ND = buf_to_int(buf, i)

        elseif blk in (0x0000383046313430, 0x0000383046343430, 0x0000373046343530, 0x0000383046313630, 0x0000343146323630)
          # Number of numerators (041F08, 044F08, 054F07, 061F08) / Number of coefficients (062F14)
          NN = buf_to_int(buf, i)

        elseif blk in (0x0000353146333430, 0x0000343146333530)
          # Number of poles (043F15, 053F14)
          NP = buf_to_int(buf, i)

        elseif blk in (0x0000303146333430, 0x0000393046333530)
          # Number of zeroes (043F10, 053F09)
          NZ = buf_to_int(buf, i)

        elseif blk in (0x0000363046313430, 0x0000363046333430, 0x0000363046343430, 0x0000353046333530, 0x0000353046343530, 0x0000363046313630, 0x0000353046323630)
          # Response in units lookup (041F06, 043F06, 044F06, 053F05, 054F05, 061F06, 062F05)
          if units || isempty(C.units)
            units_in = fix_units(resp_unit_split(buf, i))
            if seq_n > 0 && units == true
              C.resp.i[seq_n] = units_in
            end
          end

        elseif blk in (0x0000373046313430, 0x0000373046333430, 0x0000373046343430, 0x0000363046333530, 0x0000363046343530, 0x0000373046313630, 0x0000363046323630)
          # Response out units lookup (041F07, 043F07, 044F07, 053F06, 054F06, 061F07, 062F06)
          if units
            units_out = fix_units(resp_unit_split(buf, i))
            if seq_n > 0
              C.resp.o[seq_n] = units_out
            end
          end

        elseif blk in (0x0000343046333530, 0x0000343046343530, 0x0000333046373530, 0x0000333046383530, 0x0000343046303630, 0x0000333046313630, 0x0000343046323630)
          # Stage sequence number (053F04, 054F04, 057F03, 058F03, 060F04, 061F03, 062F04)
          seq_n = buf_to_int(buf, i)
          nmax = max(nmax, seq_n)

          #= Here is where we dump everything to seq_n_old.

          Fortunately RESP files list stage 0 last, which makes stupid tricks
          like this viable.
          =#
          if seq_n_old != seq_n

            # initialize stage seq_n
            if seq_n > 0
              if length(C.resp.stage) < seq_n
                append!(C.resp, MultiStageResp(6))
              end
            end

            # dump to seq_n_old
            if seq_n_old > 0
              if isa(C.resp.stage[seq_n_old], CoeffResp)
                C.resp.stage[seq_n_old].b = X
                C.resp.stage[seq_n_old].a = D
                X = Float64[]
                D = Float64[]
              elseif isa(C.resp.stage[seq_n_old], PZResp64)
                C.resp.stage[seq_n_old].a0 = a0
                C.resp.stage[seq_n_old].f0 = f0
                C.resp.stage[seq_n_old].z = Z
                C.resp.stage[seq_n_old].p = P
                Z = ComplexF64[]
                P = ComplexF64[]
                a0 = one(Float64)
                f0 = one(Float64)
              end
            end

            seq_n_old = seq_n
          end

          # trigger new channel on seq_n
          if blk in (0x0000343046333530, 0x0000343046343530, 0x0000343046323630)
            if seq_n > 0 && tfc != 0x00
              add_stage!(C, seq_n, tfc)
              tfc = 0x00
            end
          end

        elseif blk == 0x0000323246323530
          # Start date (052F22)
          ts = parse_resp_date!(buf, i, d_arr)

        elseif blk == 0x0000333046303530
          # Station (050F03)
          if chip
            # close_resp_channel!(S, C, ts, te, file, corr, delay, fg, fs_in, gain, fac, os, length(C.resp.stage))
            close_resp_channel!(S, C, ts, te, file, nmax)
            C = SeisChannel(resp = MultiStageResp(12))
            R = C.resp
            chip = false
            seq_n = zero(Int64)
            seq_n_old = zero(Int64)
            tfc = 0x00
            units_in = ""
            nmax = zero(Int64)
          end
          ni = min(i,5)
          id = zeros(UInt8, ni+1)
          unsafe_copyto!(id, 2, buf, 1, ni)
          id[1] = 0x2e

        elseif blk in (0x0000353046333430, 0x0000353046343430, 0x0000333046333530, 0x0000333046343530, 0x0000333046323630)
          # Response type (043F05, 044F05) / Transfer function type (053F03, 054F03, 062F03)
          tfc = buf[1]

          # trigger on response type
          if blk in (0x0000353046333430, 0x0000353046343430)
            if seq_n > 0
              add_stage!(C, seq_n, tfc)
              tfc = 0x00
            end
          end

        elseif blk in (0x0000353046313430, 0x0000353046313630)
          # Symmetry type (041F05, 061F05) tells us we have a digital filter
          tfc = 0x44
          if seq_n > 0
            add_stage!(C, seq_n, tfc)
            tfc = 0x00
          end

        end

        if seq_n == 1 && isempty(C.units) && !isempty(units_in)
          C.units = units_in
          units_in = ""
        end

        read_state = 0x00
      end
    end
    close(io)
    close_resp_channel!(S, C, ts, te, file, nmax)
    if nf != lastindex(files)
      C = SeisChannel(resp = MultiStageResp(12))
      R = C.resp
      chip = false
      seq_n = zero(Int64)
      seq_n_old = zero(Int64)
      tfc = 0x00
      units_in = ""
      nmax = zero(Int64)
    end
  end
  return S
end

function read_seed_resp(fpat::String;
  units::Bool=false)
  S = SeisData()
  read_seed_resp!(S, fpat, units=units)
  return S
end

"""
    RESP_wont_read()

The following is a list of breaking SEED RESP issues that we've encountered
in real data. Files with these issues don't read correctly into any known
program (e.g., ObsPy, SAC, SeisIO).

| Network | Station(s)    | Problem(s)
| :----   | :----         | :----
| CN      | (broadbands)  | B058F05-06 contain units; should be B053F05-06
"""
function RESP_wont_read()
  return nothing
end


text:
randuniform(dims::Base.Dims) = randuniform(Float64, dims)
randuniform(::Type{T}, dims::Base.Dims) where {T<:Number} = rand(T, dims)

randnormal(dims::Base.Dims) = randnormal(Float64, dims)
randnormal(::Type{T}, dims::Base.Dims) where {T<:Number} = randn(T, dims)

randisometry(dims::Base.Dims{2}) = randisometry(Float64, dims)
randisometry(::Type{T}, dims::Base.Dims{2}) where {T<:Number} = _leftorth!(randnormal(T, dims), QRpos())[1]


text:
using Revise, PyPlot, LinearAlgebra, Sobol, ForwardDiff, Random
using MaximumEntropyBasisFunctions

nBasis = 10; nEval = 1000;

xData = 1.1*sort(reshape(2*rand(MersenneTwister(0), nBasis) .- 1,1,nBasis),dims=2);
# xData = reshape(1.05*collect(range(-1,1,length=nBasis)),1,nBasis);

xEval = reshape(collect(range(-1,1,length=nEval)),1,nEval);

d = size(xData,1);
# Generate mexent basis functions
Phi = zeros(nEval,nBasis);
derPhi = zeros(nEval,nBasis);
m = zeros(nEval,nBasis);

ϕ = MaximumEntropyBasisFunctions.CreateMaxEntBasis(xData,nn=5);

@time for i ∈ 1:nEval
	Phi[i,:] = MaximumEntropyBasisFunctions.EvaluateMaxEntBasis(ϕ,xEval[:,i])
	derPhi[i,:] = MaximumEntropyBasisFunctions.getBasisDerivative(ϕ,xEval[:,i],Phi[i,:]);
	for j in 1:nBasis
		m[i,j] = MaximumEntropyBasisFunctions.priorFunction(xEval[:,i]-ϕ.x[:,j],ϕ.β);
	end
end

# Plot ϕ and ϕ_x
nCol = 5;
if  mod(nBasis,nCol) == 0
	nRow = div(nBasis,nCol);
else
	nRow = div(nBasis,nCol)+1;
end

# Find maximum
vals, inds = findmax(Phi,dims=1);

figure(1,figsize=(10,4)); clf();
pygui(true)
rc("text", usetex=true);
rc("font", serif="Helvetica");
rc("font", size=10);
for i in 1:nBasis
	subplot(nRow,nCol,i);
	# plot(xEval,Phi[:,i],linewidth=1.0); grid("on");
	plot(xEval',Phi[:,i],linewidth=1.0,color="b");
	plot(xEval',derPhi[:,i],linewidth=1.0,color="r");
	plot(xEval',m[:,i],linewidth=1.0,color="g");
	grid("on");

	# # Plots maximum point
	# X = xEval[1,inds[i][1]];
	# plot([X;X],[0;vals[i]],color="k",lw=0.5);

	local str = "\$" * "\\phi_" * "{$i}" * "\$";
	title(LaTeXString(str));
	# axis("off");
end
tight_layout();
legend((L"\phi",L"\phi_x",L"e^{-\beta x^2}"));

max_stars_repo_name: isrlab/MaximumEntropyBasisFunctions

max_issues_repo_name: isrlab/MaximumEntropyBasisFunctions

max_forks_repo_name: isrlab/MaximumEntropyBasisFunctions


text:
module TestComponents
	using Roxygen

	srcpath = "test/example.jl"

	chunks = Roxygen.extractdocs(srcpath)

	lexemes = Roxygen.lex(chunks[1])
	doc = Roxygen.parse_plus_semantics(lexemes)
	summary = Roxygen.markdown(doc)

	# For debugging, set io = STDOUT
	io = IOBuffer()
	print(io, CharString(chunks[1]))
	print(io, summary)

	lexemes = Roxygen.lex(chunks[2])
	doc = Roxygen.parse_plus_semantics(lexemes)
	summary = Roxygen.markdown(doc)

	# For debugging, set io = STDOUT
	io = IOBuffer()
	print(io, CharString(chunks[2]))
	print(io, summary)
end


text:
import .DiffEqCallbacks: SavingCallback, SavedValues

function SavingCallback(semi::AbstractSemidiscretization; kwargs...)
  T = eltype(semi.derivative)

  save_func = (u,t,integrator) -> integrate(u->ScalarIntegralQuantities(u,u^2),
                                              u, integrator.f.f)
  saved_values = SavedValues(T, ScalarIntegralQuantities{T})
  SavingCallback(save_func, saved_values; kwargs...)
end


text:
macro makeshow(CV, fields)
    fields = fields.args
    CVstr = string(CV)
    CVesc = esc(CV)
    objfields = Array(Expr, length(fields))
    for i = 1:length(fields)
        if isa(fields[i], Expr)
            sym = fields[i].args[2].args[1]
            objfields[i] = :(c.c.$sym)
        else
            objfields[i] = :(c.$(fields[i]))
        end
    end
    exs = [:(showcompact(io, $(fn))) for fn in objfields]
    exc = [d < length(fields) ? (:(print(io, ','))) : (:(print(io, ')'))) for d = 1:length(fields)]
    exboth = hcat(exs, exc)'
    ex = Expr(:block, exboth...)
    quote
        function Base.show{T,f}(io::IO, c::$CVesc{FixedPointNumbers.UfixedBase{T,f}})
            print(io, "$($CVstr){Ufixed", f, "}(")
            $ex
        end
    end
end

@makeshow RGB (r, g, b)
@makeshow XYZ (x, y, z)
@makeshow RGBA (c.r, c.g, c.b, alpha)

# Displaying color swatches
# -------------------------

function writemime(io::IO, ::MIME"image/svg+xml", c::ColorValue)
    write(io,
        """
        <?xml version"1.0" encoding="UTF-8"?>
        <!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"
         "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
        <svg xmlns="http://www.w3.org/2000/svg" version="1.1"
             width="25mm" height="25mm" viewBox="0 0 1 1">
             <rect width="1" height="1"
                   fill="#$(hex(c))" stroke="none"/>
        </svg>
        """)
end

function writemime{T <: ColorValue}(io::IO, ::MIME"image/svg+xml",
                                    cs::AbstractVecOrMat{T})
    m,n = ndims(cs) == 2 ? size(cs) : (1,length(cs))

    xsize,xpad = n > 50 ? (250/n,0) : n > 18 ? (5.,1) : n > 12 ? (10.,1) : n > 1 ? (15.,1) : (25.,0)
    ysize,ypad = m > 28 ? (150/m,0) : m > 14 ? (5.,1) : m > 9 ? (10.,1) : m > 1 ? (15.,1) : (25.,0)

    write(io,
        """
        <?xml version="1.0" encoding="UTF-8"?>
        <!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"
         "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
        <svg xmlns="http://www.w3.org/2000/svg" version="1.1"
             width="$(n*xsize)mm" height="$(m*ysize)mm"
             shape-rendering="crispEdges">
        """)

    for i in 1:m, j in 1:n
        c = ndims(cs) == 2 ? cs[i, j] : cs[j]
        write(io,
            """
            <rect x="$((j-1)*xsize)mm" y="$((i-1)*ysize)mm"
                  width="$(xsize - xpad)mm" height="$(ysize - ypad)mm"
                  fill="#$(hex(c))" stroke="none" />
            """)
    end

    write(io, "</svg>")
end


text:
@testset "Reduce" begin
    let rng = MersenneTwister(123456)
        import Nabla.fmad

        # Check that `mapreduce`, `mapfoldl`and `mapfoldr` work as expected with all unary
        # functions, some composite functions which use FMAD under both `+` and `*`.
        let N = 3
            for functional in (mapreduce, mapfoldl, mapfoldr)

                # Sensitivities implemented in Base.
                for (package, f) in Nabla.unary_sensitivities

                    # Generate some data and get the function to be mapped.
                    f = eval(f)
                    domain = domain1(f)
                    domain === nothing && error("Could not determine domain for $f.")
                    lb, ub = domain
                    x = rand(rng, N) .* (ub - lb) .+ lb

                    # Test +.
                    x_ = Leaf(Tape(), x)
                    s = functional(f, +, x_)
                    @test ∇(s)[x_] ≈ ∇.(f, Arg{1}, x)
                end

                # Some composite sensitivities.
                composite_functions = (x->5x, x->1 / (1 + x), x->10+x)
                for f in composite_functions

                    # Generate some data.
                    x = randn(rng, 100)

                    # Test +.
                    x_ = Leaf(Tape(), x)
                    s = functional(f, +, x_)
                    @test unbox(s) ≈ functional(f, +, x)
                    @test ∇(s)[x_] ≈ map(x->fmad(f, (x,), Val{1}), x)
                end
            end
        end

        # Check that something sensible (if rather slow) happens for some arbitrary `op`s.
        let

        end

        # Check that `reduce`, `foldl` and `foldr` work as expected for `+` and `*`.
        let
            for functional in (reduce, foldl, foldr)

                # Test `+`.
                x = randn(rng, 100)
                x_ = Leaf(Tape(), x)
                s = functional(+, x_)
                @test unbox(s) == functional(+, x)
                @test ∇(s)[x_] ≈ oneslike(100)
            end
        end

        # Check that `sum` and `prod` work as expected under all implemented unary functions
        # and some composite functions which use FMAD.
        let N = 5
            # Sensitivities implemented in Base.
            for (package, f) in Nabla.unary_sensitivities

                # Generate some data and get the function to be mapped.
                f = eval(f)
                domain = domain1(f)
                domain === nothing && error("Could not determine domain for $f.")
                lb, ub = domain
                x = rand(rng, N) .* (ub .- lb) .+ lb

                # Test +.
                x_ = Leaf(Tape(), x)
                s = sum(f, x_)
                @test unbox(s) == sum(f, x)
                @test ∇(s)[x_] ≈ ∇.(f, Arg{1}, x)
            end

            # Some composite functions.
            composite_functions = (x->5x, x->1 / (1 + x), x->10+x)
            for f in composite_functions

                # Generate some data.
                x = randn(rng, N)

                # Test +.
                x_ = Leaf(Tape(), x)
                s = sum(f, x_)
                @test unbox(s) == sum(f, x)
                @test ∇(s)[x_] ≈ map(x->fmad(f, (x,), Val{1}), x)
            end
        end
    end
end


text:
"""
    interpret_odds(odds::Number; rules::Union{String, Rules}="chen2010", logodds::Bool=false)

Interpret odds ratio with rules of thumb.

The odds ratio (OR) is probably one of the most widely used index of effect size.

# Arguments
- `odds`: Odds ratio value.
- `rules`: Can be 'chen2010' or custom set of [Rules](@ref).
- `logodds`: Are these log odds ratio (*e.g.*, extracted from a logistic regression)?

# Rules of thumb
- **Chen et al., (2010)**: Thee authors demonstrate that OR = 1.68, 3.47, and 6.71 are equivalent to Cohen's d = 0.2 (small), 0.5 (medium), and 0.8 (large) and Cohen's d < 0.2 when OR < 1.5, and Cohen's d > 0.8 when OR > 5.

# Examples
```jldoctest
interpret_odds(2, rules="chen2010")

# output

"small"

# References
- Chen, H., Cohen, P., & Chen, S. (2010). How big is a big odds ratio? Interpreting the magnitudes of odds ratios in epidemiological studies. Communications in Statistics—Simulation and Computation, 39(4), 860-864.
```
"""
function interpret_odds(odds::Number; rules::Union{String, Rules}="chen2010", logodds::Bool=false)

    if logodds == true
        odds = exp(abs(odds))
    end

    if isa(rules, String)
        if rules == "chen2010"
            rules = Rules(
                [1.68, 3.47, 6.71],
                ["very small", "small", "medium", "large"],
                true
                )
        else
            throw(ArgumentError("`rules` must be a String ('chen2010') or custom Rules."))
        end
    end

    interpretation = interpret(abs(odds), rules)
    return interpretation
end


text:
@recipe function f(bath::OhmicBath, func_symbol, ω)
    if (func_symbol == :γ) || (func_symbol == :spectrum)
        yguide --> "\$\\frac{\\gamma(\\omega)}{\\eta g^2 \\omega_c}\$"
    elseif func_symbol == :S
        yguide --> "\$\\frac{S(\\omega)}{\\eta g^2 \\omega_c}\$"
    else
        throw(ArgumentError("Must be :γ or :S for Ohmic bath."))
    end
    y = [eval(Expr(:call, func_symbol, x, bath)) for x in ω * bath.ωc]
    xguide --> "\$\\omega / \\omega_c\$"
    (ω, y ./ bath.η / bath.ωc)
end

@recipe function f(bath::EnsembleFluctuator, func_symbol, x_axis)
    y = [eval(Expr(:call, func_symbol, x, bath)) for x in x_axis]
    (x_axis, y)
end


text:
using GeoDataFrames; const GDF = GeoDataFrames
using Test
using DataFrames
using GeoFormatTypes; const GFT = GeoFormatTypes
using Pkg.PlatformEngines

# Use ArchGDAL datasets to test with
probe_platform_engines!()  # for download

const testdatadir = @__DIR__

REPO_URL = "https://github.com/yeesian/ArchGDALDatasets/blob/master/"

remotefiles = [
    ("ospy/data1/sites.dbf", "7df95edea06c46418287ae3430887f44f9116b29715783f7d1a11b2b931d6e7d"),
    ("ospy/data1/sites.prj", "81fb1a246728609a446b25b0df9ede41c3e7b6a133ce78f10edbd2647fc38ce1"),
    ("ospy/data1/sites.sbn", "198d9d695f3e7a0a0ac0ebfd6afbe044b78db3e685fffd241a32396e8b341ed3"),
    ("ospy/data1/sites.sbx", "49bbe1942b899d52cf1d1b01ea10bd481ec40bdc4c94ff866aece5e81f2261f6"),
    ("ospy/data1/sites.shp", "69af5a6184053f0b71f266dc54c944f1ec02013fb66dbb33412d8b1976d5ea2b"),
    ("ospy/data1/sites.shx", "1f3da459ccb151958743171e41e6a01810b2a007305d55666e01d680da7bbf08"),
]
@info "Downloading test files..."
for (f, sha) in remotefiles
    localfn = joinpath(testdatadir, basename(f))
    url = REPO_URL * f * "?raw=true"
    PlatformEngines.download_verify(url, sha, localfn; force=true)
end

@testset "GeoDataFrames.jl" begin
    fn = "sites.shp"
    coords = zip(rand(10), rand(10))

    @testset "Read actual shapefile" begin
        # Read large file
        t = GDF.read(fn)
        @test nrow(t) == 42
        @test "ID" in names(t)
    end

    @testset "Read self written file" begin
        # Save table with a few random points
        table = DataFrame(geom=createpoint.(coords), name="test")
        GDF.write("test_points.shp", table)
        GDF.write("test_points.gpkg", table, "test_points")
        GDF.write("test_points.geojson", table, "test_points")

        ntable = GDF.read("test_points.shp")
        @test nrow(ntable) == 10
        ntable = GDF.read("test_points.gpkg")
        @test nrow(ntable) == 10
        ntable = GDF.read("test_points.geojson")
        @test nrow(ntable) == 10
    end

    @testset "Write actual shapefile" begin

        t = GDF.read(fn)

        # Save table from reading
        GDF.write("test_read.shp", t, "test_coastline")
        GDF.write("test_read.gpkg", t, "test_coastline")
        GDF.write("test_read.geojson", t, "test_coastline")

    end

    @testset "Spatial operations" begin
        table = DataFrame(geom=createpoint.(coords), name="test")

        # Buffer to also write polygons
        table.geom = buffer(table.geom, 10)
        GDF.write("test_polygons.shp", table)
        GDF.write("test_polygons.gpkg", table)
        GDF.write("test_polygons.geojson", table)
    end

    @testset "Reproject" begin
        table = DataFrame(geom=createpoint.([[0,0,0]]), name="test")
        reproject(table.geom, GFT.EPSG(4326), GFT.EPSG(28992))
        @test GDF.AG.getpoint(table.geom[1], 0)[1] ≈ -587791.596556932
    end
end


text:
# SIF.jl, adds an imread function for Andor .sif images
# 2013 Ronald S. Rock, Jr.

import Images.imread
function imread{S<:IO}(stream::S, ::Type{Images.AndorSIF})
    # line 1
    seek(stream, 0)
    l = strip(readline(stream))
    l == "Andor Technology Multi-Channel File" || error("Not an Andor file: " * l)

    # line 2
    l = strip(readline(stream))
    l == "65538 1" || error("Unknown Andor version number at line 2: " * l)

    # line 3: TInstaImage thru "Head model"
    l = strip(readline(stream))
    fields = split(l)
    fields[1] == "65547" || fields[1] == "65558" ||
        error("Unknown TInstaImage version number at line 3: " * fields[1])

    ixon = {"data_type" => int(fields[2])}
    ixon["active"] = int(fields[3])
    ixon["structure_vers"] = int(fields[4]) # (== 1)
    # date is recored as seconds counted from 1970.1.1 00:00:00
    ixon["date"] = int(fields[5]) # need to convert to actual date
    ixon["temperature"] = max(int(fields[6]), int(fields[48]))
    ixon["temperature_stable"] = int(fields[6]) != -999
    ixon["head"] = fields[7]
    ixon["store_type"] = fields[8]
    ixon["data_type"] = fields[9]
    ixon["mode"] = fields[10]
    ixon["trigger_source"] = fields[11]
    ixon["trigger_level"] = fields[12]
    ixon["exposure_time"] = float(fields[13])
    ixon["frame_delay"] = float(fields[14])
    ixon["integration_cycle_time"] = float(fields[15])
    ixon["no_integrations"] = int(fields[16])
    ixon["sync"] = fields[17]
    ixon["kin_cycle_time"] = float(fields[18])
    ixon["pixel_readout_time"] = float(fields[19])
    ixon["no_points"] = int(fields[20])
    ixon["fast_track_height"] = int(fields[21])
    ixon["gain"] = int(fields[22])
    ixon["gate_delay"] = float(fields[23])
    ixon["gate_width"] = float(fields[24])
    ixon["gate_step"] = float(fields[25])
    ixon["track_height"] = int(fields[26])
    ixon["series_length"] = int(fields[27])
    ixon["read_pattern"] = fields[28]
    ixon["shutter_delay"] = fields[29]
    ixon["st_center_row"] = int(fields[30])
    ixon["mt_offset"] = int(fields[31])
    ixon["operation_mode"] = fields[32]
    ixon["flipx"] = fields[33]
    ixon["flipy"] = fields[34]
    ixon["clock"] = fields[35]
    ixon["aclock"] = fields[36]
    ixon["MCP"] = fields[37]
    ixon["prop"] = fields[38]
    ixon["IOC"] = fields[39]
    ixon["freq"] = fields[40]
    ixon["vert_clock_amp"] = fields[41]
    ixon["data_v_shift_speed"] = float(fields[42])
    ixon["output_amp"] = fields[43]
    ixon["pre_amp_gain"] = float(fields[44])
    ixon["serial"] = int(fields[45])
    ixon["num_pulses"] = int(fields[46])
    ixon["m_frame_transfer_acq_mode"] = int(fields[47])
    ixon["unstabilized_temperature"] = int(fields[48])
    ixon["m_baseline_clamp"] = int(fields[49])
    ixon["m_pre_scan"] = int(fields[50])
    ixon["m_em_real_gain"] = int(fields[51])
    ixon["m_baseline_offset"] = int(fields[52])
    _ = fields[53]
    _ = fields[54]
    ixon["sw_vers1"] = int(fields[55])
    ixon["sw_vers2"] = int(fields[56])
    ixon["sw_vers3"] = int(fields[57])
    ixon["sw_vers4"] = int(fields[58])

    # line 4
    ixon["camera_model"] = strip(readline(stream))

    # line 5, something like camera dimensions??
    _ = readline(stream)

    # line 6
    ixon["original_filename"] = strip(readline(stream))

    # line 7
    l = strip(readline(stream))
    fields = split(l)
    fields[1] == "65538" || error("Unknown TUserText version number in line 7: $fields[1]")
    usertextlen = int(fields[2]) # don't need?

    # line 8
    usertext = strip(readline(stream))
    # ixon["usertext"] = usertext # Not useful

    # line 9 TShutter
    l = strip(readline(stream)) # Weird!

    # line 10 TCalibImage
    l = strip(readline(stream))

    # lines 11-22
    _ = strip(readline(stream))
    _ = strip(readline(stream))
    _ = strip(readline(stream))
    _ = strip(readline(stream))
    _ = strip(readline(stream))
    _ = strip(readline(stream))
    _ = strip(readline(stream))
    _ = strip(readline(stream))
    _ = strip(readline(stream))
    _ = strip(readline(stream))
    _ = strip(readline(stream))

    # what a bizarre file format here
    # length of the next string is in this line
    next_str_len = int(strip(readline(stream)))
    # and here is the next string, followed by the length
    # of the following string, with no delimeter in between!
    l = strip(readline(stream))
    next_str_len = int(l[(next_str_len + 1):])
    # lather, rinse, repeat...
    l = strip(readline(stream))
    next_str_len = int(l[(next_str_len + 1):])
    l = strip(readline(stream))
    l = l[(next_str_len + 1):]
    fields = split(l)
    fields[1] == "65538" || error("Unknown version number at image dims record")
    ixon["image_format_left"] = int(fields[2])
    ixon["image_format_top"] = int(fields[3])
    ixon["image_format_right"] = int(fields[4])
    ixon["image_format_bottom"] = int(fields[5])
    frames = int(fields[6])
    ixon["frames"] = frames
    ixon["num_subimages"] = int(fields[7])
    ixon["total_length"] = int(fields[8]) # in pixels across all frames
    ixon["single_frame_length"] = int(fields[9])

    # Now at the first (and only) subimage
    l = strip(readline(stream))
    fields = split(l)
    fields[1] == "65538" || error("unknown TSubImage version number: " * fields[1])
    left = int(fields[2])
    top = int(fields[3])
    right = int(fields[4])
    bottom = int(fields[5])
    vertical_bin = int(fields[6])
    horizontal_bin = int(fields[7])
    subimage_offset = int(fields[8])

    # calculate frame width, height, with binning
    width = right - left + 1
    mod = width%horizontal_bin
    width = int((width - mod)/horizontal_bin)
    height = top - bottom + 1
    mod = height%vertical_bin
    height = int((height - mod)/vertical_bin)

    ixon["left"] = left
    ixon["top"] = top
    ixon["right"] = right
    ixon["bottom"] = bottom
    ixon["vertical_bin"] = vertical_bin
    ixon["horizontal_bin"] = horizontal_bin
    ixon["subimage_offset"] = subimage_offset

    # rest of the header is a timestamp for each frame
    # (actually, just a bunch of zeros). Skip
    for i = 1:frames
        _ = readline(stream)
    end
    offset = position(stream) # start of the actual pixel data, 32-bit float, little-endian

    pixels = read(stream, Float32, width, height, frames)
    prop = {"colorspace" => "Gray",
            "spatialorder" => ["y", "x"],
            "ixon" => ixon,
            "suppress" => Set("ixon"),
            "pixelspacing" => [1, 1]}
    if frames > 1
        prop["timedim"] = 3
    end
    Image(float64(pixels), prop)
end


text:
"""functions to grab the physiome model repository.

to just get a folder of cellml files, use `cellml_models`

"""
module CellMLModelRepository

using JSON3, JSONTables, CSV, DataFrames
using Base.Threads, Downloads

using HTTP, Cascadia, Gumbo, EzXML
using CellMLToolkit, OrdinaryDiffEq

const datadir = joinpath(@__DIR__, "../data/")
mkpath(datadir)

include("curl.jl")
include("repos.jl")

export cellml_models, cellml_metadata
export cellml_repo_table, cellml_workspaces, clone_physiome, run_all_repos

end


text:
import Pkg; Pkg.add("HDF5")
import Pkg; Pkg.add("ArgParse")
import Pkg; Pkg.add("HypergeometricFunctions")
import Pkg; Pkg.add("SpecialFunctions")
import Pkg; Pkg.add("StaticArrays")
import Pkg; Pkg.add("Interpolations")
import Pkg; Pkg.add("PolynomialRoots")

max_stars_repo_name: KerwannTEP/CAT

max_issues_repo_name: KerwannTEP/CAT

max_forks_repo_name: KerwannTEP/CAT


text:
"""
Generate the low-level `src/AWSServices.jl` file with all definitions.
"""
function _generate_low_level_wrappers(services::AbstractArray, repo_name::String, auth::GitHub.OAuth2)
    service_definitions = sort(_generate_low_level_definitions(services, repo_name, auth))

    template = """
    # This file is auto-generated by AWSMetadata.jl

    module AWSServices

    using AWS
    using OrderedCollections: LittleDict

    $(join(service_definitions, "\n"))

    end
    """

    open(services_path, "w") do f
        print(f, template)
    end

    return template
end

"""
Get the low-level definitions for all AWS Services.
"""
function _generate_low_level_definitions(services::AbstractArray, repo_name::String, auth::GitHub.OAuth2)
    service_definitions = String[]

    for service in services
        service_name = service["path"]
        @info "Generating low level wrapper for $service_name"

        # Get the contents of the ${Service}.normal.json file
        service_blob = blob(repo_name, service["sha"]; auth=auth)
        service = JSON.parse(String(base64decode(service_blob.content)))
        service_metadata = service["metadata"]

        definition = _generate_low_level_definition(service_metadata)
        push!(service_definitions, definition)
    end

    return service_definitions
end

"""
Get the low-level definition for an AWS Service.
"""
function _generate_low_level_definition(service::AbstractDict)
    protocol = service["protocol"]
    service_name = haskey(service, "signingName") ? service["signingName"] : service["endpointPrefix"]
    service_id = replace(lowercase(service["serviceId"]), ' ' => '_')
    api_version = service["apiVersion"]

    service_specifics = LittleDict{String, String}()

    if service_id == "glacier"
        service_specifics[service_id] = "LittleDict(\"x-amz-glacier-version\" => \"$(service["apiVersion"])\")"
    end

    if protocol == "rest-xml"
        return "const $service_id = AWS.RestXMLService(\"$service_name\", \"$api_version\")"
    elseif protocol in ("ec2", "query")
        return "const $service_id = AWS.QueryService(\"$service_name\", \"$api_version\")"
    elseif protocol == "rest-json" && haskey(service_specifics, service_id)
        return "const $service_id = AWS.RestJSONService(\"$service_name\", \"$api_version\", $(service_specifics[service_id]))"
    elseif protocol == "rest-json"
        return "const $service_id = AWS.RestJSONService(\"$service_name\", \"$api_version\")"
    elseif protocol == "json"
        json_version = service["jsonVersion"]
        target = service["targetPrefix"]
        return "const $service_id = AWS.JSONService(\"$service_name\", \"$api_version\", \"$json_version\", \"$target\")"
    else
        throw(ProtocolNotDefined("$service_id is using a new protocol; $protocol which is not supported."))
    end
end


text:
# TODO: Optimize for fg! calls
# TODO: Implement safeguards

function interpolating_linesearch!{T}(d::Union{DifferentiableFunction,
                                               TwiceDifferentiableFunction},
                                      x::Vector,
                                      p::Vector,
                                      x_new::Vector,
                                      gr_new::Vector,
                                      lsr::LineSearchResults{T},
                                      c::Real,
                                      mayterminate::Bool;
                                      c1::Real = 1e-4,
                                      c2::Real = 0.9,
                                      rho::Real = 2.0)

    # Parameter space
    n = length(x)

    # Count function calls
    f_calls, g_calls = 0, 0

    # Step-sizes
    a_0 = 0.0
    a_iminus1 = a_0
    a_i = 1.0
    a_max = 65536.0

    # phi(alpha) = f(x + alpha * p)
    phi_0 = d.f(x)
    f_calls += 1
    phi_a_iminus1 = phi_0
    phi_a_i = NaN

    # phi'(alpha) = vecdot(g(x + alpha * p), p)
    d.g!(x, gr_new)
    g_calls += 1
    phiprime_0 = vecdot(gr_new, p)
    phiprime_a_i = NaN

    # Iteration counter
    i = 1

    while a_i < a_max
        # Update x_new
        @simd for index in 1:n
            @inbounds x_new[index] = x[index] + a_i * p[index]
        end

        # Evaluate phi(a_i)
        phi_a_i = d.f(x_new)
        f_calls += 1

        # Test Wolfe conditions
        if (phi_a_i > phi_0 + c1 * a_i * phiprime_0) ||
             (phi_a_i >= phi_a_iminus1 && i > 1)
            a_star, f_up, g_up = zoom(a_iminus1, a_i,
                                      phiprime_0, phi_0,
                                      d.f, d.g!, x, p, x_new, gr_new)
            return a_star, f_calls + f_up, g_calls + g_up
        end

        # Evaluate phi'(a_i)
        d.g!(x_new, gr_new)
        g_calls += 1
        phiprime_a_i = vecdot(gr_new, p)

        # Check condition 2
        if abs(phiprime_a_i) <= -c2 * phiprime_0
            return a_i, f_calls, g_calls
        end

        # Check condition 3
        if phiprime_a_i >= 0.0
            a_star, f_up, g_up = zoom(a_i, a_iminus1,
                                      phiprime_0, phi_0,
                                      d.f, d.g!, x, p, x_new, gr_new)
            return a_star, f_calls + f_up, g_calls + g_up
        end

        # Choose a_iplus1 from the interval (a_i, a_max)
        a_iminus1 = a_i
        a_i *= rho

        # Update phi_a_iminus1
        phi_a_iminus1 = phi_a_i

        # Update iteration count
        i += 1
    end

    # Quasi-error response
    return a_max, f_calls, g_calls
end

function zoom(a_lo::Real,
              a_hi::Real,
              phiprime_0::Real,
              phi_0::Real,
              f::Function,
              g!::Function,
              x::Vector,
              p::Vector,
              x_new::Vector,
              gr_new::Vector;
              c1::Real = 1e-4,
              c2::Real = 0.9)

    # Parameter space
    n = length(x)

    # Count function and gradient calls
    f_calls, g_calls = 0, 0

    # Step-size
    a_j = NaN

    # Count iterations
    iteration = 0
    max_iterations = 10

    # Shrink bracket
    while iteration < max_iterations
        iteration += 1

        # Cache phi_a_lo
        @simd for index in 1:n
            @inbounds x_new[index] = x[index] + a_lo * p[index]
        end
        phi_a_lo = f(x_new)
        g!(x_new, gr_new)
        f_calls += 1
        g_calls += 1
        phiprime_a_lo = vecdot(gr_new, p)

        # Cache phi_a_hi
        @simd for index in 1:n
            @inbounds x_new[index] = x[index] + a_hi * p[index]
        end
        phi_a_hi = f(x_new)
        g!(x_new, gr_new)
        f_calls += 1
        g_calls += 1
        phiprime_a_hi = vecdot(gr_new, p)

        # Interpolate a_j
        if a_lo < a_hi
            a_j = interpolate(a_lo, a_hi,
                              phi_a_lo, phi_a_hi,
                              phiprime_a_lo, phiprime_a_hi)
        else
            # TODO: Check if this is needed
            a_j = interpolate(a_hi, a_lo,
                              phi_a_hi, phi_a_lo,
                              phiprime_a_hi, phiprime_a_lo)
        end

        # Update x_new
        @simd for index in 1:n
            @inbounds x_new[index] = x[index] + a_j * p[index]
        end

        # Evaluate phi(a_j)
        phi_a_j = f(x_new)
        f_calls += 1

        # Check Armijo
        if (phi_a_j > phi_0 + c1 * a_j * phiprime_0) ||
             (phi_a_j > phi_a_lo)
            a_hi = a_j
        else
            # Evaluate phiprime(a_j)
            g!(x_new, gr_new)
            g_calls += 1
            phiprime_a_j = vecdot(gr_new, p)

            if abs(phiprime_a_j) <= -c2 * phiprime_0
                return a_j, f_calls, g_calls
            end

            if phiprime_a_j * (a_hi - a_lo) >= 0.0
                a_hi = a_lo
            end

            a_lo = a_j
        end
    end

    # Quasi-error response
    return a_j, f_calls, g_calls
end

# a_lo = a_{i - 1}
# a_hi = a_{i}
function interpolate(a_i1::Real, a_i::Real,
                     phi_a_i1::Real, phi_a_i::Real,
                     phiprime_a_i1::Real, phiprime_a_i::Real)
    d1 = phiprime_a_i1 + phiprime_a_i -
           3.0 * (phi_a_i1 - phi_a_i) / (a_i1 - a_i)
    d2 = sqrt(d1 * d1 - phiprime_a_i1 * phiprime_a_i)
    return a_i - (a_i - a_i1) *
            ((phiprime_a_i + d2 - d1) /
             (phiprime_a_i - phiprime_a_i1 + 2.0 * d2))
end