text:
# Triangulations of point sets
# ===========================
#
#
# These examples can be loaded into Julia (Revise.jl recommended)
# and run by calling one of the methods with the optional arguments "Plotter=PyPlot".
# Alternatively, you can download a [jupyter notebook](example_pointsets.ipynb) created
# from this source.
#
#
#---
#
# Set  up environment

using Triangulate
using Test

injupyter()=(isdefined(Main, :IJulia) && Main.IJulia.inited)

if injupyter()
    import PyPlot
end

# ### Delaunay triangulation of point set
#
# Create a set of random points in the plane and calculate
# the Delaunay triangulation of this set of points. It is a triangulation
# where for each triangle, the interior of its circumcircle  does not contain any points
# of the trianglation.
#
# The Delaunay triangulation of a set of points in general position
# (no 4 of them on a circle) is unique. At the same time, it is a
# triangulation of the convex hull of these points.
#
# Given an input list of points, without any further flags, Triangle creates
# just this triangulation (the "Q" flag suppresses the text output of Triangle).
# For this and the next examples, the input list of points is created randomly,
# but on a raster, preventing the appearance of too close points.
#
function example_convex_hull(;Plotter=nothing,n=10,raster=10)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=hcat(unique([ Cdouble[rand(1:raster)/raster, rand(1:raster)/raster] for i in 1:n])...)
    display(triin)
    (triout, vorout)=triangulate("Q", triin)
    display(triout)
    plot_in_out(Plotter,triin,triout,title="Convex hull")
    @test numberofpoints(triin)==numberofpoints(triout)
    @test numberoftriangles(triout)>0
end
#
injupyter()&& example_convex_hull(Plotter=PyPlot,n=10,raster=10);

# ### Delaunay triangulation of point set with boundary
#
# Same as the previous example, but in addition specify the "c" flag
# In this case, Triangle outputs an additional list of segments
# describing the boundary of the convex hull. In fact this is a constrained
# Delaunay triangulation (CDT) where the boundary segments
# are the seen constraining edges which must appear in the output.
function example_convex_hull_with_boundary(;Plotter=nothing,n=10,raster=10)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=hcat(unique([ Cdouble[rand(1:raster)/raster, rand(1:raster)/raster] for i in 1:n])...)
    display(triin)
    (triout, vorout)=triangulate("cQ", triin)
    display(triout)
    @test numberofpoints(triin)==numberofpoints(triout)
    @test numberoftriangles(triout)>0
    @test numberofsegments(triout)>0
    plot_in_out(Plotter,triin,triout,title="Convex hull with boundary")
end
#
injupyter()&&  example_convex_hull_with_boundary(Plotter=PyPlot,n=10,raster=10);

# ### Delaunay triangulation of point set with Voronoi diagram
#
# Same as the previous example, but instead of "c" specify the "v" flag
# In this case, Triangle outputs information about the Voronoi diagram
# of the point set which is a structure dual to  the Delaunay triangulation.
#
# The Voronoi cell around a point $p$ a point set $S$ is defined
# as the set of points $x$ such that $|x-p|<|x-q|$ for all $q\in S$
# such that $p\neq q$.
# The Voronoi cells of boundary points
# of the convex hull are of infinite size. The corners of the Voronoi
# cells are the circumcenters of the triangles. They can be far
# outside of the triangulated domain.
function example_convex_hull_voronoi(;Plotter=nothing,n=10,raster=10)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=hcat(unique([ Cdouble[rand(1:raster)/raster, rand(1:raster)/raster] for i in 1:n])...)
    display(triin)
    (triout, vorout)=triangulate("vQ", triin)
    display(triout)
    display(vorout)
    plot_in_out(Plotter,triin,triout,voronoi=vorout,title="Convex hull with Voronoi diagram")
    @test numberofpoints(triin)==numberofpoints(triout)
    @test numberoftriangles(triout)>0
    @test numberofpoints(vorout)>0
end
#
injupyter()&&  example_convex_hull_voronoi(Plotter=PyPlot,n=10,raster=10);

# ### Boundary conforming Delaunay triangulation of point set
# Specify "c" flag for convex hull segments, "v" flag for Voronoi
# and "D" flag for creating a boundary conforming Delaunay triangulation of
# the point set. In this case additional points are created which split
# the boundary segments and ensure that all triangle circumcenters
# lie within the convex hull.
# Due to random input, there may be situations where Triangle fails with this task,
# so we check for the corresponding exception.
function example_convex_hull_voronoi_delaunay(;Plotter=nothing,n=10,raster=10)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=hcat(unique([ Cdouble[rand(1:raster)/raster, rand(1:raster)/raster] for i in 1:n])...)
    display(triin)
    try
        (triout, vorout)=triangulate("vcDQ", triin)
        plot_in_out(Plotter,triin,triout,voronoi=vorout,title="Convex hull with Voronoi diagram")
        display(triout)
        display(vorout)
        @test numberofpoints(triin)<=numberofpoints(triout)
        @test numberoftriangles(triout)>0
        @test numberofpoints(vorout)>0
    catch err
        if typeof(err)==TriangulateError
            println("Triangle had some problem.")
            return true
        end
    end
end
#
injupyter()&&  example_convex_hull_voronoi_delaunay(Plotter=PyPlot,n=10,raster=10);

# ### Constrained Delaunay triangulation (CDT) of a point set with edges
# Constrained Delaunay triangulation (CDT) of a point set with
# additional constraints given a priori. This is obtained when
# specifying the "p" flag and an additional list of segments each described
# by two points
# which should become edges of the triangulation. Note that
# the resulting triangulation is not Delaunay in the sense
# given above.
function example_cdt(;Plotter=nothing,n=10,raster=10)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=hcat(unique([ Cdouble[rand(1:raster)/raster, rand(1:raster)/raster] for i in 1:n])...)
    npt=size(triin.pointlist,2)
    triin.segmentlist=Matrix{Cint}([1 2; npt-1 npt-2;  1 npt;]')
    triin.segmentmarkerlist=Vector{Cint}([2,3,4])
    display(triin)
    (triout, vorout)=triangulate("pcQ", triin)
    display(triout)
    plot_in_out(Plotter,triin,triout,title="CDT")
    @test numberofpoints(triin)<=numberofpoints(triout)
    @test numberofsegments(triout)>=numberofsegments(triin)
    @test numberoftriangles(triout)>0
end
#
injupyter()&&  example_cdt(Plotter=PyPlot,n=10,raster=10);


text:
using StochasticIntegrators
using GeometricProblems.KuboOscillator
using Test

Δt = .1
t0 = 0.
x0 = rand(2)
q0 = rand(1)
p0 = q0.^2
λ0 = rand(1)
λ1 = rand(1)
v0 = rand(2)
y0 = rand(1)
z0 = rand(1)
ΔW = rand(3)
ΔZ = rand(3)

@testset "$(rpad("Atomic Solution Constructors",80))" begin
    sde   = kubo_oscillator_sde_1()
    psde  = kubo_oscillator_psde_1()
    spsde = kubo_oscillator_spsde_1()

    @test typeof(AtomicSolution(sde))   <: AtomicSolutionSDE
    @test typeof(AtomicSolution(psde))  <: AtomicSolutionPSDE
    @test typeof(AtomicSolution(spsde)) <: AtomicSolutionPSDE
end

@testset "$(rpad("Atomic SDE Solution",80))" begin
    asol = AtomicSolutionSDE(t0, x0, ΔW, ΔZ)
    @test get_solution(asol) == (zero(t0), zero(x0))

    set_solution!(asol, (t0, x0))
    @test get_solution(asol) == (t0, x0)
    @test asol.t  == t0
    @test asol.q  == x0
    @test asol.t̅ == zero(t0)
    @test asol.q̅ == zero(x0)

    set_increments!(asol, (ΔW, ΔZ))
    δW = zero(ΔW)
    δZ = zero(ΔZ)
    get_increments!(asol, δW, δZ)
    @test get_increments(asol) == (ΔW, ΔZ)
    @test asol.ΔW == ΔW
    @test asol.ΔZ == ΔZ
    @test δW == ΔW
    @test δZ == ΔZ

    reset!(asol, Δt)
    @test asol.t̅ == t0
    @test asol.q̅ == x0

    update!(asol, v0)
    @test asol.t == t0  + Δt
    @test asol.q == x0 .+ v0
end

@testset "$(rpad("Atomic PSDE Solution",80))" begin
    asol = AtomicSolutionPSDE(t0, q0, p0, ΔW, ΔZ)
    @test get_solution(asol) == (zero(t0), zero(q0), zero(p0))
    @test get_increments(asol) == (zero(ΔW), zero(ΔZ))

    set_solution!(asol, (t0, q0, p0))
    @test get_solution(asol) == (t0, q0, p0)
    @test asol.t  == t0
    @test asol.q  == q0
    @test asol.p  == p0
    @test asol.t̅  == zero(t0)
    @test asol.q̅  == zero(q0)
    @test asol.p̅  == zero(p0)

    set_increments!(asol, (ΔW, ΔZ))
    δW = zero(ΔW)
    δZ = zero(ΔZ)
    get_increments!(asol, δW, δZ)
    @test get_increments(asol) == (ΔW, ΔZ)
    @test asol.ΔW == ΔW
    @test asol.ΔZ == ΔZ
    @test δW == ΔW
    @test δZ == ΔZ

    reset!(asol, Δt)
    @test asol.t̅ == t0
    @test asol.q̅ == q0
    @test asol.p̅ == p0

    update!(asol, y0, z0)
    @test asol.t == t0  + Δt
    @test asol.q == q0 .+ y0
    @test asol.p == p0 .+ z0
end

ext: jl


text:
__precompile__()

module PlotlyJS

using Base.Iterators
using JSON
using Blink
using Colors
using DocStringExtensions
using DataFrames

import Base: ==

# import LaTeXStrings and export the handy macros
using LaTeXStrings
export @L_mstr, @L_str

# export some names from JSON
export json

# globals for this package
const _js_path = joinpath(dirname(dirname(@__FILE__)),
                          "deps", "plotly-latest.min.js")
const _js_cdn_path = "https://cdn.plot.ly/plotly-latest.min.js"
const _mathjax_cdn_path =
    "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS-MML_SVG"

const _autoresize = [true]
autoresize(b::Bool) = (_autoresize[1] = b; b)
autoresize() = _autoresize[1]

_isijulia() = isdefined(Main, :IJulia) && Main.IJulia.inited

_symbol_dict(x) = x
_symbol_dict(d::Associative) =
    Dict{Symbol,Any}([(Symbol(k), _symbol_dict(v)) for (k, v) in d])

# include these here because they are used below
include("traces_layouts.jl")
include("styles.jl")
abstract type AbstractPlotlyDisplay end

# core plot object
mutable struct Plot{TT<:AbstractTrace}
    data::Vector{TT}
    layout::AbstractLayout
    divid::Base.Random.UUID
    style::Style
end

# include the rest of the core parts of the package
include("display.jl")
include("util.jl")
include("json.jl")
include("subplots.jl")
include("api.jl")
include("savefig.jl")
include("convenience_api.jl")
include("recession_bands.jl")

# Set some defaults for constructing `Plot`s
function Plot(;style::Style=DEFAULT_STYLE[1])
    Plot(GenericTrace{Dict{Symbol,Any}}[], Layout(), Base.Random.uuid4(), style)
end

function Plot{T<:AbstractTrace}(data::AbstractVector{T}, layout=Layout();
                                style::Style=DEFAULT_STYLE[1])
    Plot(data, layout, Base.Random.uuid4(), style)
end

function Plot(data::AbstractTrace, layout=Layout();
              style::Style=DEFAULT_STYLE[1])
    Plot([data], layout; style=style)
end

function docs()
    schema_path = joinpath(dirname(dirname(@__FILE__)), "deps", "schema.html")
    if !isfile(schema_path)
        msg = "schema docs not built. Run `Pkg.build(\"PlotlyJS\")` to generate"
        error(msg)
    end
    w = Blink.Window()
    for f in [
        "http://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/css/bootstrap.min.css",
        "https://ajax.googleapis.com/ajax/libs/jquery/1.12.4/jquery.min.js",
        "http://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/js/bootstrap.min.js"
        ]
        Blink.load!(w.content, f)
    end
    Blink.content!(w, "html", open(readstring, schema_path), fade=false)
end

# NOTE: we export trace constructing types from inside api.jl
# NOTE: we export names of shapes from traces_layouts.jl
export

    # core types
    Plot, GenericTrace, Layout, Shape, ElectronDisplay, JupyterDisplay,
    ElectronPlot, JupyterPlot, AbstractTrace, AbstractLayout,

    # other methods
    savefig, svg_data, png_data, jpeg_data, webp_data, autoresize,

    # plotly.js api methods
    restyle!, relayout!, addtraces!, deletetraces!, movetraces!, redraw!,
    extendtraces!, prependtraces!, purge!, to_image, download_image,

    # non-!-versions (forks, then applies, then returns fork)
    restyle, relayout, addtraces, deletetraces, movetraces, redraw,
    extendtraces, prependtraces,

    # helper methods
    plot, fork, vline, hline, attr,

    # new trace types
    stem,

    # convenience stuff
    add_recession_bands!,

    # frontend methods
    init_notebook,

    # styles
    use_style!, style, Style, Cycler

## borrowed from Compose.jl
function isinstalled(pkg, ge=v"0.0.0-")
    try
        # Pkg.installed might throw an error,
        # we need to account for it to be able to precompile
        ver = Pkg.installed(pkg)
        ver == nothing && try
            # Assume the version is new enough if the package is in LOAD_PATH
            ex = Expr(:import, Symbol(pkg))
            @eval $ex
            return true
        catch
            return false
        end
        return ver >= ge
    catch
        return false
    end
end

@static if isinstalled("Cairo") && isinstalled("Rsvg")
    include("savefig_cairo.jl")
else
    function _savefig_cairo(x...)
        lib_fn = joinpath(Base.LOAD_CACHE_PATH[1], "PlotlyJS.ji")
        msg = """
        Cairo and Rsvg need to be installed in order to save in this format

        Use `Pkg.add("Rsvg")` to install both of them.

        You then need to delete $(lib_fn) and restart your Julia session
        """
        error(msg)
    end
end

function __init__()
    # --------------------------------------------- #
    # Code to run once when the notebook starts up! #
    # --------------------------------------------- #
    if !isfile(_js_path)
        info("plotly.js javascript libary not found -- downloading now")
        include(joinpath(dirname(_js_path), "build.jl"))
    end

    if _isijulia()

        init_notebook()

        @eval begin
            import IJulia
            import IJulia.CommManager: Comm, send_comm
        end

        # set up the comms we will use to send js messages to be executed
        global const _ijulia_eval_comm = Ref(Comm(:plotlyjs_eval))
        global const _ijulia_return_comms = ObjectIdDict()

        @eval begin
            function IJulia.display_dict(p::JupyterPlot)
                if p.view.displayed
                    Dict()
                else
                    p.view.displayed = true
                    Dict("text/html" => html_body(p),
                         "application/vnd.plotly.v1+json" => json(p))
                end
            end

            SyncPlot(p::Plot) = SyncPlot(p, JupyterDisplay(p))

            IJulia.display_dict(p::Plot) =
                Dict("text/plain" => sprint(show, "text/plain", p))

        end
    else
        @eval SyncPlot(p::Plot) = SyncPlot(p, ElectronDisplay(p))
    end

    global const DEFAULT_STYLE = [_default_style()]

end

end # module


text:
for (gemm, elt) in NNlib.gemm_datatype_mappings
    @eval begin
        @inline function NNlib.gemm!(transA::Val, transB::Val,
                               M::Int, N::Int, K::Int,
                               alpha::$(elt), A::Ptr{$elt}, B::Ptr{$elt},
                               beta::$(elt), C::Ptr{$elt})
            # Convert our compile-time transpose marker to a char for BLAS
            convtrans(V::Val{false}) = 'N'
            convtrans(V::Val{true})  = 'T'

            if transA == Val(false)
                lda = M
            else
                lda = K
            end
            if transB == Val(false)
                ldb = K
            else
                ldb = N
            end
            ldc = M
            ccall((@blasfunc($(gemm)), libblas), Nothing,
                  (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                   Ref{BlasInt}, Ref{$elt}, Ptr{$elt}, Ref{BlasInt},
                   Ptr{$elt}, Ref{BlasInt}, Ref{$elt}, Ptr{$elt},
                   Ref{BlasInt}),
                  convtrans(transA), convtrans(transB), M, N, K,
                  alpha, A, lda, B, ldb, beta, C, ldc)
        end
    end
end


text:
using Documenter
using SequencerJ

makedocs(
    sitename = "SequencerJ",
    format = Documenter.HTML(),
    modules = [SequencerJ, Distances],
    pages = [
        "ABCs of SequencerJ" => "index.md",
        "Metrics" => "algorithms.md",
        "Examples" => "examples.md",
        "API" => "api.md"
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#


text:
module CompoundPeriods

export Period, CompoundPeriod, ReverseCompoundPeriod,
    typesof, canonical,
    Years, Months, Weeks, Days, Hours, Minutes, Seconds,
    Milliseconds, Microseconds, Nanoseconds, TimeUnits,
    years, months, weeks, days, hours, minutes, seconds,
    milliseconds, microseconds, nanoseconds

import Base: convert,
    iterate, getindex, lastindex, eltype, length,
    min, max, minmax, reverseind, string, show, reverse, fldmod,
    (==), (!=), (<=), (<), (>=), (>), isequal, isless,
    sign, signbit,
    (+), (-), (*)

import Dates: AbstractTime, Period, CompoundPeriod, canonicalize,
    Year, Month, Week, Day, Hour, Minute, Second, Millisecond, Microsecond, Nanosecond,
    year, month, week, day, hour, minute, second, millisecond, microsecond, nanosecond,
    Time, Date, DateTime

using Dates

include("timeunitper.jl")
include("reversecompound.jl")

const Diurnal   = Union{Nanosecond, Microsecond, Millisecond, Second, Minute, Hour, Day, Week}
const Compounds = Union{CompoundPeriod, ReverseCompoundPeriod}

include("iterable.jl")
include("indexable.jl")
include("misc.jl")
include("selector.jl")
include("canonical.jl")
include("compare.jl")
include("pluralperiods.jl")

end # Compound Periods


text:
const download = MLDatasets.MNIST.download

function traindata(; dir = nothing)
    X, T = MLDatasets.MNIST.traindata(dir = dir)
    reshape_data(X, T)
end

function traindata(n_samples::Int; dir = nothing)
    X, T = MLDatasets.MNIST.traindata(dir = dir)
    X₀, T₀ = random_sample(n_samples, X, T)
    reshape_data(X₀, T₀)
end

function traindata(filtering_labels::Vector{Int}; dir = nothing)
    X, T = filter_data(MLDatasets.MNIST.traindata(dir = dir)..., filtering_labels)
    reshape_data(X, T)
end

function traindata(n_samples::Int, filtering_labels::Vector{Int}; dir = nothing)
    X, T = filter_data(MLDatasets.MNIST.traindata(dir = dir)..., filtering_labels)
    X₀, T₀ = random_sample(n_samples, X, T)
    reshape_data(X₀, T₀)
end

function testdata(; dir = nothing)
    X, T = MLDatasets.MNIST.testdata(dir = dir)
    reshape_data(X, T)
end

function testdata(n_samples::Int; dir = nothing)
    X, T = MLDatasets.MNIST.testdata(dir = dir)
    X₀, T₀ = random_sample(n_samples, X, T)
    reshape_data(X₀, T₀)
end

function testdata(filtering_labels::Vector{Int}; dir = nothing)
    X, T = filter_data(MLDatasets.MNIST.testdata(dir = dir)..., filtering_labels)
    reshape_data(X, T)
end

function testdata(n_samples::Int, filtering_labels::Vector{Int}; dir = nothing)
    X, T = filter_data(MLDatasets.MNIST.testdata(dir = dir)..., filtering_labels)
    X₀, T₀ = random_sample(n_samples, X, T)
    reshape_data(X₀, T₀)
end

function filter_data(samples::T1, labels::T2,
                     filtering_labels::Vector{Int}) where {T1 <: AbstractArray, T2 <: AbstractVector}
    indices = findall(x -> x ∈ filtering_labels, labels)
    (samples[:,:,indices], labels[indices])
end

function reshape_data(samples::T1, targets::T2) where {T1 <: AbstractArray, T2 <: AbstractVector}
    r, c, N = size(samples)
    (reshape(samples, (r*c, N)), vector_to_1hot(targets))
end

function random_sample(n_samples::Int, samples::T1, labels::T2) where {T1 <: AbstractArray, T2 <: AbstractVector}
    @assert n_samples <= length(labels) "Blah."
    range = length(labels)
    indices = Random.randperm(range)[1:n_samples]
    (samples[:,:,indices], labels[indices])
end

function vector_to_1hot(targets::T) where {T <: AbstractVector}
    labels = Set(targets)
    result = zeros(eltype(targets), length(labels), length(targets))
    for (i, x) in enumerate(labels)
        result[i,:] .= targets .== x
    end
    result
end


text:
#find row based on the code, to improve this function
# and the find_column could do a hi_lo funtion
# and then call that for both find_row and find_column

function find_row(code)
  row_code = code[1:7]
  lo_n = 0
  hi_n = 127

  count = 1

  for i in row_code
     #find the bounds of the range depending if we can the upper or lower half
     if i =='F'
        hi_n = lo_n + div(hi_n+1 - lo_n,2) - 1
        if count == 7
           return hi_n
        end
     else
        lo_n = lo_n + div(hi_n+1 - lo_n,2)
        if count == 7
           return lo_n
        end
     end

    count += 1

  end

  return 1
end

function find_column(code)
  col_code = code[8:end]

  lo_n = 0
  hi_n = 7

  count = 1

  for i in col_code

     if i =='L'
        hi_n = lo_n + div(hi_n+1 - lo_n,2) - 1
        if count == 3
           return hi_n
        end
     else
        lo_n = lo_n + div(hi_n+1 - lo_n,2)
        if count == 3
           return lo_n
        end
     end

    count += 1

  end

  return 1
end

using DelimitedFiles

# reading the input
data = readdlm("input_day5")

max_id = 0

for temp_code in data
   r = find_row(temp_code)
   c = find_column(temp_code)

   #calculate id
   id = r*8 + c

   #find max id
   if id > max_id
     global max_id = id
   end

end

println(max_id)


text:
using Plots; contourf, plot

@doc raw"""
Represents a continuous wavelet coherence analysis.

# Fields:
- `coh`: The actual coherence as a complex value. The layout of this array follows the layout of
  `ContinuousWaveletTransformed`. That is, each voice is stored in a column, starting with the
  smallest scale.
- `scales`: The vector of scales at which the wavelet transform is performed.
- `wavelet`: The wavelet being used.
- `α`: The point-wise empirical signigicance level of the coherence analysis obtained by means of
  surrogate data. Holds `NaN` if no skipped.
"""
struct ContinuousWaveletCoherence <: AbstractArray{Complex{Float64}, 2}
    coh::Array{Complex{Float64}, 2}
    scales::Array{Float64, 1}
    wavelet::GenericContinuousWavelet
    α::Float64
end

@doc raw"""
    coherence(wt::ContinuousWaveletTransform,
              X::AbstractArray{Float64, 2};
              nsurrogate=0, α=0.05)

Implements the *self-coherence* of a set of time-series. The *self-coherence* is defined as
```math
 C(X)(b,a) = \frac{\left\vert E\left[\mathcal{W}[X](b,a)\right]\right\vert^2}{E\left[\left\vert\mathcal{W}[X](b,a)\right\vert^2\right]} =
  \frac{\left\vert\mathcal{W}[E\left[X\right]](b,a)\right\vert^2}{E\left[\left\vert\mathcal{W}[X](b,a)\right\vert^2\right]}
```
This measure can be interpreted as how coherent a process is from realization to realization.
Typical applications are the analyses of EEG event-related potentials and similar time-series of
repeated experiments.

# Arguments
- `wt`: Specifies the wavelet transform to use.
- `X`: Specifies the time-series to analyze. Each single time-series (trial) is stored as a column in `X`.
- 'nsurrogate': A keyword argument specifying the number of surrogate samples to generate from the
  given time-series for the estimation of the point-wise significant test. By default it is 0 and
  no surrogate samples are generated.
- `α`: A keyword arguement specifying the significance level for the point-wise significance test.
  By default set to `α=0.05`.
"""
function coherence(wt::ContinuousWaveletTransform, X::AbstractArray{Float64, 2}; nsurrogate=0, α=0.05)
    N,M = size(X)
    C = zeros(Complex{Float64}, N, length(wt.scales))
    P = zeros(Float64, N, length(wt.scales))
    for i in 1:M
        wX = transform(wt, X[:,i])
        C[:,:] += wX.wX
        P[:,:] += map(z->z^2, map(abs, wX.wX))
    end
    C = map(z->z^2, map(abs, C)) ./ M;
    C ./= P

    threshold = NaN;
    if nsurrogate > 0
        empC = Array{Float64, 1}[]
        # generate nsurrogate surrogate datasets and obtain coherence
        for i in 1:nsurrogate
            Cs = coherence(wt, surrogate(X); nsurrogate=0).coh
            empC = vcat(empC, map(abs, reshape(Cs, (prod(size(Cs)),)) ) )
        end
        # get threshold corresponding to the given α-level
        sort!(empC)
        i = Int32(floor((1.0-α)*length(empC)))
        threshold = empC[(i<1) ? 1 : i]
    end
    ContinuousWaveletCoherence(C, wt.scales, wt.wavelet, threshold)
end

@doc raw"""
    coherence(wt::ContinuousWaveletTransform,
              X::AbstractArray{Float64, 2}, Y::AbstractArray{Float64, 2};
              nsurrogate=0, α=0.05)

Canonical wavelet coherence as defined by
```math
 C(X)(b,a) = \frac{E\left[\mathcal{W}[X](b,a)\,\mathcal{W}[Y]^{*}(b,a)\right]}{
  E\left[\left\vert\mathcal{W}[X](b,a)\right\vert\,\left\vert\mathcal{W}[Y](b,a)\right\vert\right]}
```
That is the normalized expectation value of the product of the two wavelet transformed, where
`z^*` is the complex conjugate.

In contrast to the *self-coherence*, this analysis does estimate the coherence of a process to
itself but the coherence between two observed processes. It is therefore necessary that the two
sets of time-series `X` and `Y` are of the same shape. That is, each trial of `X` and `Y` must
consists of the same number of samples (rows) and both sets must contain the same number of trials
(columns).

# Arguments
- `wt`: Specifies the wavelet transform to use.
- `X`: Specifies one time-series to analyze. Each single time-series (trial) is stored as a column
    in `X`. Must be of same shape as `Y`.
- `Y`: Specifies one time-series to analyze. Each single time-series (trial) is stored as a column
    in `Y`. Must be of same shape as `X`.
- 'nsurrogate': A keyword argument specifying the number of surrogate samples to generate from the
    given time-series for the estimation of the point-wise significant test. By default it is 0
    and no surrogate samples are generated.
- `α`: A keyword arguement specifying the significance level for the point-wise significance test.
    By default set to `α=0.05`.
"""
function coherence(wt::ContinuousWaveletTransform, X::AbstractArray{Float64, 2}, Y::AbstractArray{Float64, 2}; nsurrogate=0, α=0.05)
    (size(X) != size(Y)) && error("Arrays X & Y must be of same size!");
    N,M = size(X)
    C = zeros(Complex{Float64}, N, length(wt.scales))
    P = zeros(Float64, N, length(wt.scales))
    for i in 1:M
        wX = transform(wt, X[:,i])
        wY = transform(wt, Y[:,i])
        C[:,:] += wX.wX .* conj.(wY.wX)
        P[:,:] += map(abs, wX.wX) .* map(abs, wX.wX)
    end
    C ./= P

    threshold = NaN;
    if nsurrogate > 0
        empC = Array{Float64, 1}[]
        # generate nsurrogate surrogate datasets and obtain coherence
        for i in 1:nsurrogate
            Cs = coherence(wt, surrogate(X); nsurrogate=0).coh
            empC = vcat(empC, map(abs, reshape(Cs, (prod(size(Cs)),)) ) )
        end
        # get threshold corresponding to the given α-level
        sort!(empC)
        i = Int32(floor((1.0-α)*length(empC)))
        threshold = empC[(i<1) ? 1 : i]
    end
    ContinuousWaveletCoherence(C, wt.scales, wt.wavelet, threshold)
end

Base.size(A::ContinuousWaveletCoherence) = size(A.coh)
Base.getindex(A::ContinuousWaveletCoherence, I::Vararg{Int, N}) where {N} = getindex(A.coh, I...);
Base.setindex!(A::ContinuousWaveletCoherence, v, I::Vararg{Int, N}) where {N} = setindex!(A.coh, v, I...)

@doc raw"""
    contourf(A::ContinuousWaveletCoherence;  drawvalid=true, shadevalid=true, kw...)

Implements a plotting helper function to plot wavelet-coherence analyses using a filled contour
plot.

If a significance test was performed using surrogate data during the analyses, the
point-wise significant areas are shown using a single (thicker) contour line.

If the optional keyword argument `drawvalid` is `true` (default), the area of valid wavelet
coefficients is marked by a thick black line. If the optional keyword arguement `shadevalid` is
also `true` (default) the invalid wavelet coefficients are shaded gray. A valid wavelet coefficient
is one, where the accociated wavelet does not overlapp with the edges of the time-series.

Additional keyword arguments are passed to the default implementation of `Plots.contourf()`.
"""
function Plots.contourf(A::ContinuousWaveletCoherence;  drawvalid=true, shadevalid=true, kw...)
    N,M = size(A)
    X = 1:N; Y = A.scales
    Z = transpose(map(abs, A.coh))
    cn = contourf(X, Y, Z;
                  linewidth=0, linetype=:heatmap,
                  xlims=(1,N), ylims=(minimum(A.scales),maximum(A.scales)),
                  kw...);
    # draw point-wise significance level if set
    if isfinite(A.α)
        contour!(X, Y, Z; levels=[A.α], linewidth=3, linecolor=:green)
    end
    # draw valid range
    if drawvalid
        x = [1.0]; y = Float64[A.scales[1]]
        for i in 1:M
            tc = cutoff_time(A.wavelet)*A.scales[i]
            if (tc < N÷2)
                push!(y, A.scales[i])
                push!(x, tc)
            end
        end
        append!(x, reverse(N.-x))
        append!(y, reverse(y))
        if shadevalid
            plot!(x,y;
                linecolor=:black, linewidth=3, label="",
                fillrange=A.scales[end], fillalpha=0.5, fillcolor=:black)
        else
            plot!(x,y;
                linecolor=:black, linewidth=3, label="")
        end
    end
    return cn
end


text:
#==========================================================================================#

# TYPE

mutable struct Mullahy <: GMM

    method::String
    sample::Microdata
    β::Vector{Float64}
    V::AbstractMatrix{Float64}
    W::AbstractMatrix{Float64}

    Mullahy() = new()
end

#==========================================================================================#

# CONSTRUCTOR

function Mullahy(MD::Microdata, W::AbstractMatrix{Float64}, method::String)
    obj        = Mullahy()
    obj.sample = MD
    obj.W      = W
    obj.method = method
    return obj
end

#==========================================================================================#

# INTERFACE

function fit(
        ::Type{Mullahy},
        MD::Microdata;
        novar::Bool = false,
        method::String = "One-step GMM"
    )

    if method == "Poisson"

        FSD                   = Microdata(MD)
        FSD.mapping[:control] = vcat(MD.mapping[:treatment], MD.mapping[:control])

        pop!(FSD.mapping, :treatment)
        pop!(FSD.mapping, :instrument)

        obj = Poisson(FSD)

        _fit!(obj, getweights(obj))

    elseif method == "Reduced form"

        FSD                   = Microdata(MD)
        FSD.mapping[:control] = vcat(MD.mapping[:instrument], MD.mapping[:control])

        pop!(FSD.mapping, :treatment)

        obj = Poisson(FSD)

        _fit!(obj, getweights(obj))

    else

        if length(MD.mapping[:treatment]) == length(MD.mapping[:instrument])
            k   = length(MD.mapping[:instrument]) + length(MD.mapping[:control])
            obj = Mullahy(MD, Matrix{Float64}(I, k, k), "Method of moments")
        elseif method == "One-step GMM"
            k   = length(MD.mapping[:instrument]) + length(MD.mapping[:control])
            obj = Mullahy(MD, Matrix{Float64}(I, k, k), "One-step GMM")
        elseif (method == "TSLS") | (method == "2SLS")
            W   = crossprod(getmatrix(MD, :instrument, :control), getweights(MD))
            obj = Mullahy(MD, W, "Two-step GMM")
        elseif (method == "Two-step GMM") | (method == "Optimal GMM")
            W     = crossprod(getmatrix(MD, :instrument, :control), getweights(MD))
            obj   = Mullahy(MD, W, method) ; _fit!(obj, getweights(obj))
            obj.W = wmatrix(obj, getcorr(obj), getweights(obj))
        else
            throw("unknown method")
        end

        _fit!(obj, getweights(obj))
    end

    novar || _vcov!(obj, getcorr(obj), getweights(obj))

    return obj
end

#==========================================================================================#

# ESTIMATION

function _fit!(obj::Mullahy, ::UnitWeights)

    has_offset = haskey(obj.sample.mapping, :offset)

    y  = getvector(obj, :response)
    x  = getmatrix(obj, :treatment, :control)
    z  = getmatrix(obj, :instrument, :control)
    W  = obj.W * nobs(obj)^2
    μ  = Array{Float64}(undef, length(y))
    xx = Array{Float64}(undef, size(x)...)

    if isdefined(obj, :β)
        β₀ = obj.β
    else
        β₀ = vcat(fill(0.0, size(x, 2) - 1), log(mean(y)))
    end

    has_offset && (xo = getmatrix(obj, :offset, :treatment, :control))

    function L(β::Vector)

        has_offset ? mul!(μ, xo, vcat(1.0, β)) : mul!(μ, x, β)

        μ .= muladd.(y, exp.(- μ), - 1.0)
        m  = z' * μ
        mw = W \ m

        return 0.5 * dot(m, mw)
    end

    function G!(g::Vector, β::Vector)

        has_offset ? mul!(μ, xo, vcat(1.0, β)) : mul!(μ, x, β)

        μ  .= y .* exp.(- μ)
        xx .= x .* μ
        d   = - xx' * z
        μ  .= μ .- 1.0
        mw  = W \ (z' * μ)

        mul!(g, d, mw)
    end

    function LG!(g::Vector, β::Vector)

        has_offset ? mul!(μ, xo, vcat(1.0, β)) : mul!(μ, x, β)

        μ  .= y .* exp.(- μ)
        xx .= x .* μ
        d   = - xx' * z
        μ  .= μ .- 1.0
        m   = z' * μ
        mw  = W \ m

        mul!(g, d, mw)

        return 0.5 * dot(m, mw)
    end

    function H!(h::AbstractMatrix, β::Vector)

        has_offset ? mul!(μ, xo, vcat(1.0, β)) : mul!(μ, x, β)

        μ  .= y .* exp.(- μ)
        xx .= x .* μ
        d   = z' * xx
        dw  = W \ d

        mul!(h, transpose(d), dw)
    end

    res = optimize(TwiceDifferentiable(L, G!, LG!, H!, β₀), β₀, Newton())

    if Optim.converged(res)
        obj.β = Optim.minimizer(res)
    else
        throw("minimization did not converge")
    end
end

function _fit!(obj::Mullahy, w::AbstractWeights)

    has_offset = haskey(obj.sample.mapping, :offset)

    y  = getvector(obj, :response)
    x  = getmatrix(obj, :treatment, :control)
    z  = getmatrix(obj, :instrument, :control)
    W  = obj.W * nobs(obj)^2
    μ  = Array{Float64}(undef, length(y))
    xx = Array{Float64}(undef, size(x)...)

    if isdefined(obj, :β)
        β₀ = obj.β
    else
        β₀ = vcat(fill(0.0, size(x, 2) - 1), log(mean(y)))
    end

    has_offset && (xo = getmatrix(obj, :offset, :treatment, :control))

    function L(β::Vector)

        has_offset ? mul!(μ, xo, vcat(1.0, β)) : mul!(μ, x, β)

        μ .= w .* muladd.(y, exp.(- μ), - 1.0)
        m  = z' * μ
        mw = W \ m

        return 0.5 * dot(m, mw)
    end

    function G!(g::Vector, β::Vector)

        has_offset ? mul!(μ, xo, vcat(1.0, β)) : mul!(μ, x, β)

        μ  .= y .* exp.(- μ)
        xx .= x .* μ .* w
        d   = - xx' * z
        μ  .= w .* (μ .- 1.0)
        mw  = W \ (z' * μ)

        mul!(g, d, mw)
    end

    function LG!(g::Vector, β::Vector)

        has_offset ? mul!(μ, xo, vcat(1.0, β)) : mul!(μ, x, β)

        μ  .= y .* exp.(- μ)
        xx .= x .* μ .* w
        d   = - xx' * z
        μ  .= w .* (μ .- 1.0)
        m   = z' * μ
        mw  = W \ m

        mul!(g, d, mw)

        return 0.5 * dot(m, mw)
    end

    function H!(h::AbstractMatrix, β::Vector)

        has_offset ? mul!(μ, xo, vcat(1.0, β)) : mul!(μ, x, β)

        μ  .= w .* y .* exp.(- μ)
        xx .= x .* μ
        d   = z' * xx
        dw  = W \ d

        mul!(h, transpose(d), dw)
    end

    res = optimize(TwiceDifferentiable(L, G!, LG!, H!, β₀), β₀, Newton())

    if Optim.converged(res)
        obj.β = Optim.minimizer(res)
    else
        throw("minimization did not converge")
    end
end

#==========================================================================================#

# SCORE (MOMENT CONDITIONS)

score(obj::Mullahy) = Diagonal(residuals(obj)) * getmatrix(obj, :instrument, :control)

# EXPECTED JACOBIAN OF SCORE × NUMBER OF OBSERVATIONS

function jacobian(obj::Mullahy, ::UnitWeights)

    y = getvector(obj, :response)
    x = getmatrix(obj, :treatment, :control)
    z = getmatrix(obj, :instrument, :control)

    if haskey(obj.sample.mapping, :offset)
        x = getmatrix(obj, :offset, :treatment, :control)
        v = xo * vcat(1.0, obj.β)
    else
        x = getmatrix(obj, :treatment, :control)
        v = x * obj.β
    end

    v .= y .* exp.(- v)

    return - z' * (x .* v)
end

function jacobian(obj::Mullahy, w::AbstractWeights)

    y = getvector(obj, :response)
    x = getmatrix(obj, :treatment, :control)
    z = getmatrix(obj, :instrument, :control)

    if haskey(obj.sample.mapping, :offset)
        x = getmatrix(obj, :offset, :treatment, :control)
        v = xo * vcat(1.0, obj.β)
    else
        x = getmatrix(obj, :treatment, :control)
        v = x * obj.β
    end

    v .= w .* y .* exp.(- v)

    return - z' * (x .* v)
end

#==========================================================================================#

# LINEAR PREDICTOR

function linear_predictor(obj::Mullahy, MD::Microdata)
    if getnames(obj, :treatment, :control) != getnames(MD, :treatment, :control)
        throw("missing variables")
    end
    if haskey(obj.sample.mapping, :offset)
        return getmatrix(MD, :offset, :treatment, :control) * vcat(1.0, obj.β)
    else
        return getmatrix(MD, :treatment, :control) * obj.β
    end
end

# FITTED VALUES

predict(obj::Mullahy, MD::Microdata) = exp.(linear_predictor(obj, MD))

# DERIVATIVE OF FITTED VALUES

function jacobexp(obj::Mullahy)

    if haskey(obj.sample.mapping, :offset)
        v = getmatrix(obj, :offset, :treatment, :control) * vcat(1.0, obj.β)
    else
        v = getmatrix(obj, :treatment, :control) * obj.β
    end

    v .= exp.(v)

    return Diagonal(v) * getmatrix(obj, :treatment, :control)
end

# RESIDUALS

function residuals(obj::Mullahy, MD::Microdata)
    r  = predict(obj, MD)
    r .= response(obj) ./ r .- 1.0
    return r
end

#==========================================================================================#

# UTILITIES

coefnames(obj::Mullahy) = getnames(obj, :treatment, :control)

max_stars_repo_name: lbittarello/Microeconometrics

max_issues_repo_name: lbittarello/Microeconometrics

max_forks_repo_name: lbittarello/Microeconometrics


text:
# automatically generated by the FlatBuffers compiler, do not modify

# module: Example2

module Example2
    include("Monster.jl")
end


text:
export AlazarDSPGetModules, AlazarDSPGetInfo, AlazarDSPGenerateWindowFunction
export AlazarDSPGetBuffer, AlazarDSPAbortCapture, AlazarFFTSetWindowFunction
export AlazarFFTVerificationMode, AlazarFFTSetup

AlazarDSPGetModules(handle::U32, numEntries, modules, numModules) =
    ccall((:AlazarDSPGetModules,ats), U32,
        (U32, U32, Ptr{dsp_module_handle}, Ptr{U32}),
        handle, numEntries, modules, numModules)

AlazarDSPGetInfo(handle::dsp_module_handle, dspModuleId,
        versionMajor, versionMinor, maxLength, reserved0, reserved1) =
    ccall((:AlazarDSPGetInfo,ats), U32, (dsp_module_handle, Ptr{U32},
        Ptr{U16}, Ptr{U16}, Ptr{U32}, Ptr{U32}, Ptr{U32}), handle, dspModuleId,
        versionMajor, versionMinor, maxLength, reserved0, reserved1)

AlazarDSPGenerateWindowFunction(windowType, window,
        windowLength_samples, paddingLength_samples) =
    ccall((:AlazarDSPGenerateWindowFunction,ats), U32,
        (U32, Ptr{Cfloat}, U32, U32), windowType, window,
        windowLength_samples, paddingLength_samples)

AlazarFFTSetWindowFunction(handle::dsp_module_handle,
        samplesPerRecord, reArray, imArray) =
    ccall((:AlazarFFTSetWindowFunction,ats), U32,
        (dsp_module_handle, U32, Ptr{Cfloat}, Ptr{Cfloat}),
        handle, samplesPerRecord, reArray, imArray)

# Undocumented in API.
AlazarFFTVerificationMode(handle::dsp_module_handle, enable,
        reArray, imArray, recordLength_samples) =
    ccall((:AlazarFFTVerificationMode,ats), U32,
        (dsp_module_handle, Bool, Ptr{S16}, Ptr{S16}, Csize_t),
        handle, enable, reArray, imArray, recordLength_samples)

AlazarFFTSetup(handle::dsp_module_handle, channel, recordLength_samples,
        fftLength_samples, outputFormat, reserved, footer, bytesPerOutputRecord) =
    ccall((:AlazarFFTSetup,ats), U32,
        (dsp_module_handle, U16, U32, U32, U32, U32, U32, Ptr{U32}),
        handle, CHANNEL_A, recordLength_samples, fftLength_samples,
        outputFormat, footer, reserved, bytesPerOutputRecord)

AlazarDSPGetBuffer(handle::U32, buffer, timeout_ms) =
    ccall((:AlazarDSPGetBuffer,ats), U32, (U32, Ptr{Void}, U32),
    handle, buffer, timeout_ms)

AlazarDSPAbortCapture(handle::U32) =
    ccall((:AlazarDSPAbortCapture,ats), U32, (U32,), handle)


text:
module pfnet

# Dependencies
if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("PFNET is not properly installed. Please run Pkg.build(\"pfnet\")")
end

# Vector
function Vector(vec::Ptr{Void}, own::Bool)
    jvec = unsafe_wrap(Array,
                       ccall((:VEC_get_data, libpfnet), Ptr{Float64}, (Ptr{Void},), vec),
                       ccall((:VEC_get_size, libpfnet), Int, (Ptr{Void},), vec),
                       own)
    if own
        Base.Libc.free(vec) # delete container
    end
    return jvec
end

function pfnet_vec_from_Array{T}(ar::Array{T,1})
    return ccall((:VEC_new_from_array, libpfnet),
                 Ptr{Void},
                 (Ptr{T}, Int,),
                 Ref(ar, 1),
                 size(ar)[1])
end

# Matrix
function SparseMatrixCOO(mat::Ptr{Void}, own::Bool)
    nnz = m = ccall((:MAT_get_nnz, libpfnet), Int, (Ptr{Void},), mat)
    I = unsafe_wrap(Array,
                    ccall((:MAT_get_row_array, libpfnet), Ptr{Int32}, (Ptr{Void},), mat),
                    nnz,
                    own)+1
    J = unsafe_wrap(Array,
                    ccall((:MAT_get_col_array, libpfnet), Ptr{Int32}, (Ptr{Void},), mat),
                    nnz,
                    own)+1
    V = unsafe_wrap(Array,
                    ccall((:MAT_get_data_array, libpfnet), Ptr{Float64}, (Ptr{Void},), mat),
                    nnz,
                    own)
    if own
        Base.Libc.free(mat) # delete container
    end
    return (I,J,V)
end

function SparseMatrixCSC(mat::Ptr{Void}, own::Bool)
    m = ccall((:MAT_get_size1, libpfnet), Int, (Ptr{Void},), mat)
    n = ccall((:MAT_get_size2, libpfnet), Int, (Ptr{Void},), mat)
    I, J, V = SparseMatrixCOO(mat, own)
    return sparse(I, J, V, m, n)
end

# Includes
include("strings.jl")
include("net.jl")
include("parser.jl")
include("function.jl")
include("constraint.jl")
include("problem.jl")

# Symbols
for name in names(pfnet, true)
    @eval export $(Symbol(name))
end

end


text:
export chanspec, seis_www, track_on!, track_off!

# Generic handler for getting data by HTTP
# Returns:
#   R::Array{UInt8,1}, either request body or error data
#   parsable::Bool, whether or not R is parsable
function get_http_req(url::String, req_info_str::String, to::Int; status_exception::Bool=false)
  (R::Array{UInt8,1}, parsable::Bool) = try
    req = request(  "GET", url, webhdr,
                    readtimeout = to,
                    status_exception = status_exception  )
    if req.status == 200
      (req.body, true)

    else
      @warn(string("Request failed", req_info_str,
      "\nRESPONSE = ", req.status, " (", statustext(req.status), ")",
      "\n\nHTTP response is in misc[\"data\"]"))
      (Array{UInt8,1}(string(req)), false)
    end

  catch err
    T = typeof(err)
    @warn(  string( "Error thrown", req_info_str,
                    "\n\nERROR = ", T,
                    "\n\nTrying to store error message in misc[\"data\"]"
                    )
          )
    msg_data::Array{UInt8,1} = Array{UInt8,1}( try; string(getfield(err, :response)); catch; try; string(getfield(err, :msg));  catch; ""; end; end )
    (msg_data, false)
  end

  return R, parsable
end

function get_http_post(url::String, body::String, to::Int; status_exception::Bool=false)
  try
    req = request(  "POST", url, webhdr, body,
                    readtimeout = to,
                    status_exception = status_exception)
    if req.status == 200
      return (req.body, true)
    else
      @warn(string( "Request failed!\nURL: ", url, "\nPOST BODY: \n", body, "\n",
                    "RESPONSE: ", req.status, " (", statustext(req.status), ")\n" ) )
      return (Array{UInt8,1}(string(req)), false)
    end

  catch err
    @warn(string( "Error thrown:\nURL: ", url, "\nPOST BODY: \n", body, "\n",
                  "ERROR TYPE: ", typeof(err), "\n" ) )
    msg_data::Array{UInt8,1} = Array{UInt8,1}( try; string(getfield(err, :response)); catch; try; string(getfield(err, :msg));  catch; ""; end; end )
    return (msg_data, false)
  end
end

datareq_summ(src::String, ID::String, d0::String, d1::String) = ("\n" * src *
  " query:\nID = " * ID * "\nSTART = " * d0 * "\nEND = " * d1)

# ============================================================================
# Utility functions not for export
hashfname(str::Array{String,1}, ext::String) = string(hash(str), ".", ext)

# Start tracking channel IDs and data lengths
"""
    track_on!(S::SeisData)

Track changes to S.id, changes to channel structure of S, and the sizes of data
vectors in S.x. Does not track data processing operations to any channel i
unless length(S.x[i]) changes for channel i.

**Warning**: If you have or suspect gapped data in any channel, do not use
ungap! while tracking is active.
"""
function track_on!(S::SeisData)
  if S.n > 0
    ids = unique(getfield(S, :id))
    nxs = zeros(Int64, S.n)
    for i = 1:S.n
      nxs[i] = length(S.x[i])
    end
    S.misc[1]["track"] = (ids, nxs)
  end
  return nothing
end

# Stop tracking channel IDs and data lengths; report which have changed
"""
    u = track_off!(S::SeisData)

Turn off tracking in S and return a boolean vector of which channels have
been added or altered significantly.
"""
function track_off!(S::SeisData)
  k = findfirst([haskey(S.misc[i],"track") for i = 1:S.n])
  if S.n == 0
    return nothing
  elseif k == nothing
    return trues(S.n)
  end
  u = falses(S.n)
  (ids, nxs) = S.misc[k]["track"]
  for (n, id) in enumerate(S.id)
    j = findfirst(x -> x == id, ids)
    if j == nothing
      u[n] = true
    else
      if nxs[j] != length(S.x[n])
        u[n] = true
      end
    end
  end
  delete!(S.misc[k], "track")
  return u
end

function savereq(D::Array{UInt8,1}, ext::String, id::String, s::String)
  if ext == "miniseed"
    ext = "mseed"
  elseif occursin("sac", ext)
    ext = "SAC"
  end
  s_str = int2tstr(tstr2int(s))
  yy = s_str[1:4]
  mm = s_str[6:7]
  dd = s_str[9:10]
  HH = s_str[12:13]
  MM = s_str[15:16]
  SS = s_str[18:19]
  nn = lpad(div(parse(Int64, s_str[21:26]), 1000), 3, '0')
  jj = lpad(md2j(yy, mm, dd), 3, '0')
  fname = join([yy, jj, HH, MM, SS, nn, id, "R", ext], '.')
  safe_isfile(fname) && @warn(string("File ", fname, " contains an identical request; overwriting."))
  f = open(fname, "w")
  write(f, D)
  close(f)
  return nothing
end

"""
| String | Source             |
|:------:|:-------------------|
|BGR    | http://eida.bgr.de |
|EMSC   | http://www.seismicportal.eu |
|ETH    | http://eida.ethz.ch |
|GEONET | http://service.geonet.org.nz |
|GFZ    | http://geofon.gfz-potsdam.de |
|ICGC   | http://ws.icgc.cat |
|INGV   | http://webservices.ingv.it |
|IPGP   | http://eida.ipgp.fr |
|IRIS   | http://service.iris.edu |
|ISC    | http://isc-mirror.iris.washington.edu |
|KOERI  | http://eida.koeri.boun.edu.tr |
|LMU    | http://erde.geophysik.uni-muenchen.de |
|NCEDC  | http://service.ncedc.org |
|NIEP   | http://eida-sc3.infp.ro |
|NOA    | http://eida.gein.noa.gr |
|ORFEUS | http://www.orfeus-eu.org |
|RESIF  | http://ws.resif.fr |
|SCEDC  | http://service.scedc.caltech.edu |
|TEXNET | http://rtserve.beg.utexas.edu |
|USGS   | http://earthquake.usgs.gov |
|USP    | http://sismo.iag.usp.br |
"""
seis_www = Dict("BGR" => "http://eida.bgr.de",
                "EMSC" => "http://www.seismicportal.eu",
                "ETH" => "http://eida.ethz.ch",
                "GEONET" => "http://service.geonet.org.nz",
                "GFZ" => "http://geofon.gfz-potsdam.de",
                "ICGC" => "http://ws.icgc.cat",
                "INGV" => "http://webservices.ingv.it",
                "IPGP" => "http://eida.ipgp.fr",
                "IRIS" => "http://service.iris.edu",
                "ISC" => "http://isc-mirror.iris.washington.edu",
                "KOERI" => "http://eida.koeri.boun.edu.tr",
                "LMU" => "http://erde.geophysik.uni-muenchen.de",
                "NCEDC" => "http://service.ncedc.org",
                "NIEP" => "http://eida-sc3.infp.ro",
                "NOA" => "http://eida.gein.noa.gr",
                "ODC" => "http://www.orfeus-eu.org",
                "ORFEUS" => "http://www.orfeus-eu.org",
                "RESIF" => "http://ws.resif.fr",
                "SCEDC" => "http://service.scedc.caltech.edu",
                "TEXNET" => "http://rtserve.beg.utexas.edu",
                "USGS" => "http://earthquake.usgs.gov",
                "USP" => "http://sismo.iag.usp.br")

fdsn_uhead(src::String) = haskey(seis_www, src) ? seis_www[src] * "/fdsnws/" : src

"""
## CHANNEL ID SPECIFICATION
Channel ID data can be passed to SeisIO web functions in three ways:

1. String: a comma-delineated list of IDs formatted `"NET.STA.LOC.CHA"` (e.g. `"PB.B004.01.BS1,PB.B004.01.BS2"`)
2. Array{String,1}: one ID per entry, formatted `"NET.STA.LOC.CHA"` (e.g. `["PB.B004.01.BS1","PB.B004.01.BS2"]`)
3. Array{String,2}: one ID per row, formatted `["NET" "STA" "LOC" "CHA"]` (e.g. `["PB" "B004" "01" "BS?"; "PB" "B001" "01" "BS?"]`)

The `LOC` field can be left blank (e.g. `"UW.ELK..EHZ", ["UW" "ELK" "" "EHZ"]`).

The allowed subfield widths before channel IDs break is identical to the FDSN
standard: NN.SSSSS.LL.CCC (network name length ≤ 2 chars, etc.)

#### SEEDLINK ONLY
For SeedLink functions (`seedlink!`, `has_stream`, etc.), channel IDs can
include a fifth field (i.e. NET.STA.LOC.CHA.T) to set the "type" flag (one of
DECOTL, for Data, Event, Calibration, blOckette, Timing, or Logs). Note that
SeedLink calibration, timing, and logs are not supported by SeisIO.
"""
function chanspec()
  return nothing
end


text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.edu
# Date: 11/23/2014

module ACASXSensorImpl

export
    initialize,
    update,

    updateSensor,

    ACASXSensor,
    ACASXSensorState,
    ACASXSensorInput,
    ACASXSensorOutput

using AbstractSensorImpl
using AbstractSensorInterfaces
using CommonInterfaces

using CASInterface

import CommonInterfaces.initialize
import CommonInterfaces.update
import AbstractSensorInterfaces.updateSensor

typealias ACASXSensorOutput CASInterface.Input #STM Input

type ACASXSensor <: AbstractSensor
  my_id::Int64 #my id
  max_intruders::Int64
  output::ACASXSensorOutput
end

function ACASXSensor(my_id::Int, max_intruders::Int)
  return ACASXSensor(my_id, max_intruders, ACASXSensorOutput(max_intruders))
end

type ACASXSensorState
  x::Float64 #East
  y::Float64 #North
  z::Float64 #Up
  vx::Float64
  vy::Float64
  vz::Float64
end

type ACASXSensorInput
  states::Vector{ACASXSensorState}
end

getListId(list_owner_id::Integer,id::Integer) = id < list_owner_id ? id : id - 1

function updateSensor(sensor::ACASXSensor, input::ACASXSensorInput)

  own_state = input.states[sensor.my_id]
  ownInput = sensor.output.ownInput

  ownInput.dz = own_state.vz
  ownInput.z = own_state.z #baro alt
  ownInput.psi = atan2(own_state.vx, own_state.vy) #zero when aligned with y-axis / north
  ownInput.h = own_state.z #agl alt
  ownInput.modes = UInt32(sensor.my_id)

  intruders = sensor.output.intruders

  for i = 1:endof(input.states)
    if i != sensor.my_id
      intr_i = getListId(sensor.my_id, i)
      intr_state = input.states[i]
      intruders[intr_i].valid = true
      intruders[intr_i].id = UInt32(i)
      intruders[intr_i].modes = UInt32(i)
      intruders[intr_i].sr = norm([own_state.x, own_state.y, own_state.z]-
                                    [intr_state.x, intr_state.y, intr_state.z]) #slant range (feet)
      intr_psi = atan2(intr_state.x - own_state.x, intr_state.y - own_state.y)
      intruders[intr_i].chi = to_plusminus_pi(intr_psi - ownInput.psi) #bearing (radians, nose is 0, clockwise is positive)
      intruders[intr_i].z = intr_state.z #altitude (feet)  #quantized

      ## These are the defaults.  Modified through coordination
      intruders[intr_i].cvc = 0x0
      intruders[intr_i].vrc = 0x0
      intruders[intr_i].vsb = 0x0
      intruders[intr_i].equipage = EQUIPAGE_TCAS
      intruders[intr_i].quant = intruders[intr_i].equipage == EQUIPAGE_TCAS ? 25 : 100
      intruders[intr_i].sensitivity_index = 0x0
      intruders[intr_i].protection_mode = 0x0
    end
  end

  return sensor.output
end

update(sensor::ACASXSensor, input) = update(sensor, convert(ACASXSensorInput, input))
update(sensor::ACASXSensor, input::ACASXSensorInput) = updateSensor(sensor, input)

function initialize(sensor::ACASXSensor)
  reset!(sensor.output)
end

#mods x to the range [-b, b]
function to_plusminus_b(x::AbstractFloat, b::AbstractFloat)
  z = mod(x, 2 * b)
  return (z > b) ? (z - 2 * b) : z
end

to_plusminus_pi(x::AbstractFloat) = to_plusminus_b(x, Float64(pi))

end


text:
### A Pluto.jl notebook ###
# v0.14.3

using Markdown
using InteractiveUtils

# ╔═╡ 5c00529c-54cb-4cae-9154-5bad2818314d
using Pkg, DrWatson

# ╔═╡ 63def520-35a3-4bdf-a333-a8f86caba00f
begin
	@quickactivate "NumericalMethodsforEngineers"
	using NumericalMethodsforEngineers
	using Test
end

# ╔═╡ 9c686554-5b34-4929-ae8a-309c33a0cc48
md" ### f-02-03.jl"

# ╔═╡ 2b6b8cd8-f5d1-4ee5-abaa-1eb25369b633
begin
	a = Float64[16 4 8;4 5 -4;8 -4 22]
	b = Float64[4, 2, 5]
	d = zeros(3)
end;

# ╔═╡ 634ecc96-cf0d-44ce-b1d3-10c99e3684d5
begin
	c = a \ b
end

# ╔═╡ b893aebf-2b1c-4e30-a47c-9095f5a8fee8
at, dt = ldlt(a, d)

# ╔═╡ 7ac63cd7-245f-45e3-ab91-0b4ed6fde3df
dt

# ╔═╡ c85a1d6c-8fa6-4992-ac93-ef28eea94828
bt = ldlfor(at, b)

# ╔═╡ 87de1def-6f4a-43c0-b793-0e9562186445
for i in 1:size(a, 1)
  at[i, :] = at[i, :] / dt[i]
end

# ╔═╡ bd666706-7a4d-4330-9061-481fe02110d0
at

# ╔═╡ 119e89db-6e66-413e-9341-063ce1280235
btt = subbac(at, bt)

# ╔═╡ 25cca764-5472-4354-9a86-6bccfe4547c9
@test btt == c

# ╔═╡ a7de228a-f876-485a-9c44-09df62672f41
md" ### End of j-02-03.jl"

# ╔═╡ Cell order:
# ╟─9c686554-5b34-4929-ae8a-309c33a0cc48
# ╠═5c00529c-54cb-4cae-9154-5bad2818314d
# ╠═63def520-35a3-4bdf-a333-a8f86caba00f
# ╠═2b6b8cd8-f5d1-4ee5-abaa-1eb25369b633
# ╠═634ecc96-cf0d-44ce-b1d3-10c99e3684d5
# ╠═b893aebf-2b1c-4e30-a47c-9095f5a8fee8
# ╠═7ac63cd7-245f-45e3-ab91-0b4ed6fde3df
# ╠═c85a1d6c-8fa6-4992-ac93-ef28eea94828
# ╠═87de1def-6f4a-43c0-b793-0e9562186445
# ╠═bd666706-7a4d-4330-9061-481fe02110d0
# ╠═119e89db-6e66-413e-9341-063ce1280235
# ╠═25cca764-5472-4354-9a86-6bccfe4547c9
# ╟─a7de228a-f876-485a-9c44-09df62672f41


text:
"""$(TYPEDEF)
Augmented Lagrangian Trajectory Optimizer (ALTRO) is a solver developed by the Robotic Exploration Lab at Stanford University.
    The solver is special-cased to solve Markov Decision Processes by leveraging the internal problem structure.

ALTRO consists of two "phases":
1) AL-iLQR: iLQR is used with an Augmented Lagrangian framework to solve the problem quickly to rough constraint satisfaction
2) Projected Newton: A collocation-flavored active-set solver projects the solution from AL-iLQR onto the feasible subspace to achieve machine-precision constraint satisfaction.

# Constructor

    ALTROSolver(prob::Problem, opts::SolverOptions; [infeasible, R_inf, kwarg_opts...])

The `infeasible` keyword is a boolean flag that specifies whether the solver should
augment the controls to make it artificially fully actuated, allowing state initialization.
The `R_inf` is the weight on these augmented controls. Any solver options can be passed
as additional keyword arguments and will be set in the solver.

# Getters
* `get_model`
* `get_objective`
* `get_constraints`
* `Altro.get_ilqr`
* `TO.get_initial_state`

# Other methods
* `Base.size`: returns `(n,m,N)`
* `TO.integration`
* `TO.is_constrained`
"""
struct ALTROSolver{T,S} <: ConstrainedSolver{T}
    opts::SolverOptions{T}
    stats::SolverStats{T}
    solver_al::AugmentedLagrangianSolver{T,S}
    solver_pn::ProjectedNewtonSolver{T}
end

function ALTROSolver(prob::Problem{Q,T}, opts::SolverOptions=SolverOptions();
        infeasible::Bool=false,
        R_inf::Real=1.0,
        solver_uncon=iLQRSolver,
        kwarg_opts...
    ) where {Q,T}
    if infeasible
        # Convert to an infeasible problem
        prob = InfeasibleProblem(prob, prob.Z, R_inf/prob.Z[1].dt)

        # Set infeasible constraint parameters
        # con_inf = get_constraints(prob).constraints[end]
        # con_inf::ConstraintVals{T,Control,<:InfeasibleConstraint}
        # con_inf.params.μ0 = opts.penalty_initial_infeasible
        # con_inf.params.ϕ = opts.penalty_scaling_infeasible
    end
    set_options!(opts; kwarg_opts...)
    stats = SolverStats{T}(parent=solvername(ALTROSolver))
    solver_al = AugmentedLagrangianSolver(prob, opts, stats; solver_uncon=solver_uncon)
    solver_pn = ProjectedNewtonSolver(prob, opts, stats)
    link_constraints!(get_constraints(solver_pn), get_constraints(solver_al))
    S = typeof(solver_al.solver_uncon)
    solver = ALTROSolver{T,S}(opts, stats, solver_al, solver_pn)
    reset!(solver)
    # set_options!(solver; opts...)
    solver
end

# Getters
@inline Base.size(solver::ALTROSolver) = size(solver.solver_pn)
@inline TO.get_trajectory(solver::ALTROSolver)::Traj = get_trajectory(solver.solver_al)
@inline TO.get_objective(solver::ALTROSolver) = get_objective(solver.solver_al)
@inline TO.get_model(solver::ALTROSolver) = get_model(solver.solver_al)
@inline get_initial_state(solver::ALTROSolver) = get_initial_state(solver.solver_al)
solvername(::Type{<:ALTROSolver}) = :ALTRO
TrajectoryOptimization.integration(solver::ALTROSolver) = integration(solver.solver_al)
is_constrained(solver::ALTROSolver) = !isempty(get_constraints(solver.solver_al))
@inline get_ilqr(solver::ALTROSolver) = solver.solver_al.solver_uncon

function TO.get_constraints(solver::ALTROSolver)
    if solver.opts.projected_newton
        get_constraints(solver.solver_pn)
    else
        get_constraints(solver.solver_al)
    end
end

# Solve Methods
function solve!(solver::ALTROSolver)
    reset!(solver)
    conSet = get_constraints(solver.solver_al)

    if isempty(conSet)
        ilqr = solver.solver_al.solver_uncon
        solve!(ilqr)
        terminate!(solver)
        return solver
    end

    # Set terminal condition if using projected newton
    opts = solver.opts
    ϵ_con = solver.opts.constraint_tolerance
    if opts.projected_newton
        opts_al = solver.solver_al.opts
        if opts.projected_newton_tolerance >= 0
            opts_al.constraint_tolerance = opts.projected_newton_tolerance
        else
            opts_al.constraint_tolerance = 0
            opts_al.kickout_max_penalty = true
        end
    end

    # Solve with AL
    solve!(solver.solver_al)

    if status(solver) <= SOLVE_SUCCEEDED || opts.force_pn
        # Check convergence
        i = solver.solver_al.stats.iterations
        if i > 1
            c_max = solver.solver_al.stats.c_max[i]
        else
            c_max = TO.max_violation(solver.solver_al)
        end

        opts.constraint_tolerance = ϵ_con
        if (opts.projected_newton && c_max > opts.constraint_tolerance &&
                (status(solver) <= SOLVE_SUCCEEDED || status(solver) == MAX_ITERATIONS_OUTER)) ||
                opts.force_pn
            solve!(solver.solver_pn)
        end

        # Back-up check
        if status(solver) <= SOLVE_SUCCEEDED
            # TODO: improve this check
            if TO.max_violation(solver) < solver.opts.constraint_tolerance
                solver.stats.status = SOLVE_SUCCEEDED
            end
        end
    end

    terminate!(solver)
    solver
end

# Infeasible methods
function InfeasibleProblem(prob::Problem{RK}, Z0::Traj, R_inf::Real) where RK
    @assert !isnan(sum(sum.(states(Z0))))

    n,m,N = size(prob)  # original sizes

    # Create model with augmented controls
    model_inf = InfeasibleModel(prob.model)

    # Get a trajectory that is dynamically feasible for the augmented problem
    #   and matches the states and controls of the original guess
    Z = infeasible_trajectory(model_inf, Z0)

    # Convert constraints so that they accept new dimensions
    conSet = TO.change_dimension(get_constraints(prob), n, m+n, 1:n, 1:m)

    # Constrain additional controls to be zero
    inf = InfeasibleConstraint(model_inf)
    TO.add_constraint!(conSet, inf, 1:N-1)

    # Infeasible Objective
    obj = infeasible_objective(prob.obj, R_inf)

    # Create new problem
    Problem{RK}(model_inf, obj, conSet, prob.x0, prob.xf, Z, N, prob.t0, prob.tf)
end

function infeasible_objective(obj::Objective, regularizer)
    n,m = TO.state_dim(obj.cost[1]), TO.control_dim(obj.cost[1])
    Rd = [@SVector zeros(m); @SVector fill(regularizer,n)]
    R = Diagonal(Rd)
    cost_inf = TO.DiagonalCost(Diagonal(@SVector zeros(n)), R, checks=false)
    costs = map(obj.cost) do cost
        cost_idx = TO.change_dimension(cost, n, n+m, 1:n, 1:m)
        cost_idx + cost_inf
    end
    TO.Objective(costs)
end


text:
if !isdefined(Main, :InterpolationTestUtils)
    include("InterpolationTestUtils.jl")
    @eval using Main.InterpolationTestUtils
end

using Test, SharedArrays, Random
using StaticArrays, WoodburyMatrices
ambs = detect_ambiguities(StaticArrays, WoodburyMatrices, Base, Core)

using Interpolations
@test isempty(setdiff(detect_ambiguities(Interpolations, Base, Core), ambs))

const isci = get(ENV, "CI", "") in ("true", "True")

@testset "Interpolations" begin
    include("core.jl")

    # b-spline interpolation tests
    include("b-splines/runtests.jl")
    isci && println("finished b-spline")
    include("nointerp.jl")
    # extrapolation tests
    include("extrapolation/runtests.jl")

    # scaling tests
    include("scaling/runtests.jl")

    # monotonic tests
    include("monotonic/runtests.jl")

    # lanczos tests
    include("lanczos/runtests.jl")

    # test gradient evaluation
    include("gradient.jl")
    isci && println("finished gradient")
    # test hessian evaluation
    include("hessian.jl")
    isci && println("finished hessian")

    # gridded interpolation tests
    include("gridded/runtests.jl")

    # test interpolation with specific types
    include("typing.jl")

    include("issues/runtests.jl")

    include("io.jl")
    include("convenience-constructors.jl")
    include("readme-examples.jl")

end


text:
export ComponentReference

import Base.setindex!, Base.getindex

"""
A container for a component, for interacting with it within a model.
"""
type ComponentReference
    m::Model
    component::Symbol
end

"""
Set a component parameter as `setparameter(reference, name, value)`.
"""
function setparameter(c::ComponentReference, name::Symbol, value)
    setparameter(c.m, c.component, name, value)
end

"""
Set a component parameter as `reference[symbol] = value`.
"""
function setindex!(c::ComponentReference, value, name::Symbol)
    setparameter(c.m, c.component, name, value)
end

"""
Connect two components as `connectparameter(reference1, name1, reference2, name2)`.
"""
function connectparameter(target::ComponentReference, target_name::Symbol, source::ComponentReference, source_name::Symbol)
    connectparameter(target.m, target.component, target_name, source.component, source_name)
end

"""
Connect two components as `connectparameter(reference1, reference2, name)`.
"""
function connectparameter(target::ComponentReference, source::ComponentReference, name::Symbol)
    connectparameter(target.m, target.component, name, source.component, name)
end

"""
A container for a name within a component, to improve connectparameter aesthetics.
"""
type VariableReference
    m::Model
    component::Symbol
    name::Symbol
end

"""
Get a variable reference as `reference[name]`.
"""
function getindex(c::ComponentReference, name::Symbol)
    VariableReference(c.m, c.component, name)
end

"""
Connect two components as `reference1[name1] = reference2[name2]`.
"""
function setindex!(target::ComponentReference, source::VariableReference, name::Symbol)
    connectparameter(target.m, target.component, name, source.component, source.name)
end


text:
using LinearAlgebra

mutable struct AndersonAcceleration{R <: Real, C <: Union{R, Complex{R}}, I <: Integer, T <: AbstractArray{C}, M}
    currmem::I
    curridx::I
    s::T
    y::T
    s_M::Vector{T}
    y_M::Vector{T}
end

function AndersonAcceleration(x::T, M::I) where {R <: Real, C <: Union{R, Complex{R}}, I <: Integer, T <: AbstractArray{C}}
    s_M = [zero(x) for i = 1:M]
    y_M = [zero(x) for i = 1:M]
    s = zero(x)
    y = zero(x)
    AndersonAcceleration{R, C, I, T, M}(0, 0, s, y, s_M, y_M)
end

function update!(L::AndersonAcceleration{R, C, I, T, M}, s, y) where {R, C, I, T, M}
    L.s .= s
    L.y .= y
    L.curridx += 1
    if L.curridx > M L.curridx = 1 end
    L.currmem += 1
    if L.currmem > M L.currmem = M end
    copyto!(L.s_M[L.curridx], L.s)
    copyto!(L.y_M[L.curridx], L.y)
    return L
end

function reset!(L::AndersonAcceleration{R, C, I, T, M}) where {R, C, I, T, M}
    L.currmem, L.curridx = zero(I), zero(I)
end

import Base: *

function (*)(L::AndersonAcceleration, v)
    w = similar(v)
    mul!(w, L, v)
end

import LinearAlgebra: mul!

function mul!(d::T, L::AndersonAcceleration{R, C, I, T, M}, v::T) where {R, C, I, T, M}
    if L.currmem == 0
        d .= v
    else
        # TODO: optimize this
        S = hcat(L.s_M[1:L.currmem]...)
        Y = hcat(L.y_M[1:L.currmem]...)
        # H = I + (S - Y)inv(Y'Y)Y'
        d .= v .+ (S - Y)*(pinv(Y'*Y)*(Y'*v))
    end
end


text:
# ---
# title: 576. Out of Boundary Paths
# id: problem576
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming, Depth-first Search
# link: <https://leetcode.com/problems/out-of-boundary-paths/description/>
# hidden: true
# ---
#
# There is an **m** by **n** grid with a ball. Given the start coordinate
# **(i,j)** of the ball, you can move the ball to **adjacent** cell or cross the
# grid boundary in four directions (up, down, left, right). However, you can
# **at most** move **N** times. Find out the number of paths to move the ball
# out of grid boundary. The answer may be very large, return it after mod 109 \+
# 7.
#
#
#
# **Example 1:**
#
#
#
#     Input: m = 2, n = 2, N = 2, i = 0, j = 0
#     Output: 6
#     Explanation:
#     ![](https://assets.leetcode.com/uploads/2018/10/13/out_of_boundary_paths_1.png)
#
#
# **Example 2:**
#
#
#
#     Input: m = 1, n = 3, N = 3, i = 0, j = 1
#     Output: 12
#     Explanation:
#     ![](https://assets.leetcode.com/uploads/2018/10/12/out_of_boundary_paths_2.png)
#
#
#
#
# **Note:**
#
#   1. Once you move the ball out of boundary, you cannot move it back.
#   2. The length and height of the grid is in range [1,50].
#   3. N is in range [0,50].
#
#
## @lc code=start
using LeetCode

function find_paths(m, n, N, i, j)
    dp = zeros(Int, m, n)

    for _ in 1:N
        t = copy(dp)
        inds = CartesianIndices(t)
        for x in inds
            v = 0
            for dir in ((0, -1), (0, 1), (1, 0), (-1, 0))
                p = x + CartesianIndex(dir)
                v += p ∈ inds ? t[p] : 1
            end
            dp[x] = v
        end
    end
    return dp[i, j] % (10^9 + 7)
end
## @lc code=end


text:
# ------------------------------------------------------------------------

# fallback method for coercing `AbstractArrays` to instances of a given MLJ
# scitype `T<:Union{Missing, <:ST.Known}`.
# This fallback fixes issue #39
coerce(X::AbstractArray,
       ::Type{T}; kw...) where {T<:Union{Missing,Nothing,ST.Known}} =
    throw(CoercionError(
        "Coercion of eltype `$(eltype(X))` to element scitype "*
        "`$T` is not supported. "))

# general fallback method for coercing `AbstractArrays` to scitype `T`
coerce(X::AbstractArray, ::Type{T}; kw...) where {T} =
    throw(CoercionError(
        "Coercion of eltype `$(eltype(X))` to element scitype "*
        "`$T` is not supported. Indeed, `$T` does not appear to be a "*
        "scientific type. "))

# ------------------------------------------------------------------------
# FINITE

# Supported types for CategoricalArray{T} under CategoricalArrays 0.9:
const SupportedTypes = Union{Char,AbstractString,Number}

# Arr{T} -> Finite for supported raw label types:
function coerce(v::Arr{T},
                ::Type{T2};
                verbosity::Int=1, tight::Bool=false
                ) where T <: Union{SupportedTypes, Missing} where T2 <: Union{Missing,Finite}
    v1    = _tighten_if_needed(v, T, tight)
    vcat = categorical(v1, ordered=nonmissing(T2)<:OrderedFactor)
    return _finalize_finite_coerce(vcat, verbosity, T2, T)
end

# Arr{T} -> Finite for all other types:
function coerce(v::Arr{T}, ::Type{T2};
                verbosity::Int=1, tight::Bool=true
                ) where T where T2 <: Union{Missing,Finite}
    tight || verbosity < 0 ||
        @warn "Forcing `tight=true`, as `$T` unsupported by "*
        "CategoricalArrays. "
    v1 = _tighten_if_needed(v, T, true)
    if v1 isa Arr{<:Union{SupportedTypes, Missing}}
        vcat = categorical(v1, ordered=nonmissing(T2)<:OrderedFactor)
        return _finalize_finite_coerce(vcat, verbosity, T2, T)
    end
    verbosity < 0 || @warn "Converting array elements to strings before "*
        "wrapping in a `CategoricalArray`, as `$T` unsupported by "*
        "CategoricalArrays. "
#   error("######################## GOT YOU ##########################")
    v_str = string.(v1)
    vcat = categorical(v_str, ordered=nonmissing(T2)<:OrderedFactor)
    return _finalize_finite_coerce(vcat, verbosity, T2, T)
end

# CArr{T} -> Finite
function coerce(v::CArr{T}, ::Type{T2};
                verbosity::Int=1, tight::Bool=false
                ) where T where T2 <: Union{Missing,Finite}
    v = _tighten_if_needed(v, T, tight)
    return _finalize_finite_coerce(v, verbosity, T2, T)
end

# ------------------------------------------------------------------------
# INFINITE

# Arr{Int} -> {Count} is no-op
function coerce(y::Arr{T}, T2::Type{<:Union{Missing,Count}};
                verbosity::Int=1, tight::Bool=false
                ) where T <: Union{Missing,Integer}
    y = _tighten_if_needed(y, T, tight)
    _check_eltype(y, T2, verbosity)
    return y
end

# Arr{Real \ Int} -> {Count} via `_int`, may throw InexactError
function coerce(y::Arr{T}, T2::Type{<:Union{Missing,Count}};
                verbosity::Int=1, tight::Bool=false
                ) where T <: Union{Missing,Real}
    y = _tighten_if_needed(y, T, tight)
    _check_eltype(y, T2, verbosity)
    return _int.(y)
end

# Arr{Float} -> Float is no-op
function coerce(y::Arr{T}, T2::Type{<:Union{Missing,Continuous}};
                verbosity::Int=1, tight::Bool=false
                ) where T <: Union{Missing,AbstractFloat}
    y = _tighten_if_needed(y, T, tight)
    _check_eltype(y, T2, verbosity)
    return y
end

# Arr{Real \ {Float}} -> Float via `float`
function coerce(y::Arr{T}, T2::Type{<:Union{Missing,Continuous}};
                verbosity::Int=1, tight::Bool=false
                ) where T <: Union{Missing,Real}
    y = _tighten_if_needed(y, T, tight)
    _check_eltype(y, T2, verbosity)
    return float(y)
end

# CArr -> Count/Continuous via `_int` or `float(_int)`
function coerce(y::CArr{T}, T2::Type{<:Union{Missing,C}};
                verbosity::Int=1, tight::Bool=false
                ) where T where C <: Infinite
    # here we broadcast and so we don't need to tighten
    iy = _int.(y)
    _check_eltype(iy, T2, verbosity)
    C == Count && return iy
    return float(iy)
end

const MaybeNumber = Union{Missing,AbstractChar,AbstractString}

# Textual => Count / Continuous via parse
function coerce(y::Arr{T}, T2::Type{<:Union{Missing,C}};
                verbosity::Int=1, tight::Bool=false
                ) where T <: MaybeNumber where C <: Infinite
    y = _tighten_if_needed(y, T, tight)
    _check_eltype(y, T2, verbosity)
    C == Count && return _int.(y)
    return _float.(y)
end

## ARRAY OF ANY
# Note: in the categorical case, we don't care, because we broadcast anyway.
# see CArr --> C above.
#
# this is the case where the data may have been badly encoded and resulted
# in an Any[] array a user should proceed with caution here in particular:
#   - if at one point it encounters a type for which there is no AbstractFloat
#     such as a String, it will error.
#   - if at one point it encounters a Char it will **not** error but return a
#     float corresponding to the Char (e.g. 65.0 for 'A') whence the warning
# Also the performances of this should not be expected to be great as we're
# broadcasting an operation on the full vector.
function coerce(y::Arr{Any}, T::Type{<:Union{Missing,C}};
                verbosity=1, tight::Bool=false
                ) where C <: Union{Count,Continuous}
    # to float or to count?
    op, num   = ifelse(C == Count, (_int, "65"), (float, "65.0"))
    has_chars = findfirst(e -> isa(e, Char), y) !== nothing
    if has_chars && verbosity > 0
        @info "Char value encountered, such value will be coerced according to the corresponding numeric value (e.g. 'A' to $num)."
    end
    # broadcast the operation
    c = op.(y)
    # if the container type has  missing but not target, warn
    if (eltype(c) >: Missing) && !(T >: Missing) && verbosity > 0
        @info "Trying to coerce from `Any` to `$T` but encountered missing values.\nCoerced to `Union{Missing,$T}` instead."
    end
    return c
end

# ------------------------------------------------------------------------
# UTILITIES

# If trying to  coerce to a non-union type `T` from a type that >: Missing
# for instance  coerce([missing,1,2], Continuous) will throw a warning
# to avoid that do coerce([missing,1,2], Union{Missing,Continuous})
# Special case with Any which is >: Missing depending on categorical case
function _coerce_missing_warn(::Type{T}, from::Type) where T
    T >: Missing && return
    if from == Any
        @info "Trying to coerce from `Any` to `$T` with categoricals.\n" *
              "Coerced to `Union{Missing,$T}` instead."
    else
        @info "Trying to coerce from `$from` to `$T`.\n" *
              "Coerced to `Union{Missing,$T}` instead."
    end
    return
end

# v is already categorical here, but may need `ordering` changed
function _finalize_finite_coerce(v, verbosity, T, fromT)
    elst = elscitype(v)
    if (elst >: Missing) && !(T >: Missing)
        verbosity > 0 && _coerce_missing_warn(T, fromT)
    end
    if elst <: T
        return v
    end
    return categorical(v, ordered=nonmissing(T)<:OrderedFactor)
end

_int(::Missing)  = missing
_int(x::Integer) = x
_int(x::Cat)     = levelcode(x)
_int(x)          = Int(x)                # NOTE: may throw InexactError

_int(x::AbstractString)   = Int(Meta.parse(x))    # NOTE: may fail
_float(x::AbstractString) = float(Meta.parse(x))
_float(x::Missing) = missing

function _check_eltype(y, T, verb)
    E = eltype(y)
    E >: Missing && verb > 0 && _coerce_missing_warn(T, E)
end

function _tighten_if_needed(v::Arr, T, tight)
    if T >: Missing && tight && findfirst(ismissing, v) === nothing
        v = identity.(v)
    end
    return v
end

function _tighten_if_needed(v::CArr, T, tight)
    if T >: Missing && tight && findfirst(ismissing, v) === nothing
        v = CategoricalArrays.DataAPI.unwrap.(v)
    end
    return v
end

## Images
F = Real
_2d{C} =AbstractArray{C} where C<:Union{ColorTypes.AbstractRGB, ColorTypes.Gray}

# Single Image
function coerce(y::Arr{<:F, 2}, T2::Type{GrayImage})
    return ColorTypes.Gray.(y)
end

function coerce(y::Arr{<:F, 3}, T2::Type{ColorImage})
    return broadcast(ColorTypes.RGB, y[:,:,1], y[:,:,2], y[:,:,3])
end

# Collection

_3Dcollection = AbstractArray{<:F, 3}       #for b/w images
_4Dcollection = AbstractArray{<:F, 4}       #for color images

function coerce(y::_3Dcollection, T2::Type{GrayImage})
    return [ColorTypes.Gray.(y[:,:,idx]) for idx=1:size(y,3)]
end

function coerce(y::_4Dcollection, T2::Type{GrayImage})
    size(y, 3) == 1 || error("Multiple color channels encountered. "*
                      "Perhaps you want to use `coerce(image_collection, ColorImage)`.")
    y = dropdims(y, dims=3)
    return [ColorTypes.Gray.(y[:,:,idx]) for idx=1:size(y,3)]
end

function coerce(y::_4Dcollection, T2::Type{ColorImage})
    return [broadcast(ColorTypes.RGB, y[:,:,1, idx], y[:,:,2,idx], y[:,:,3, idx]) for idx=1:size(y,4)]
end


text:
export constructtriad
export constructtorus
export constructsphere
export constructcylinder
export constructbox
export constructwhirl
export constructframe
export getplane
export constructfiber
export constructtwospinor
export transformg

"""
    constructtriad(q [; length])

Construct a triad with the given configuration `q` and `length`.
"""
function constructtriad(q::Biquaternion; length::Float64 = 1.0)
    map(x -> x + gettranslation(q), rotate(ℝ³.([[0; 0; 0], [length; 0; 0],
                                                [0; 0; 0], [0; length; 0],
                                                [0; 0; 0], [0; 0; length]]),
                                           getrotation(q)))
end

"""
    constructtorus(q, r, R)

Construct a torus of revolution with the given configuration `q`, the smaller radius `r`
and the bigger radius `R`.
"""
function constructtorus(q::Biquaternion,
                        r::Real,
                        R::Real;
                        segments::Int = 36)
    array = Array{ℝ³,2}(undef, segments, segments)
    for i in 1:segments
        for j in 1:segments
            ϕ = i * 2pi / (segments - 1)
            θ = j * 2pi / (segments - 1)
            x₁ = (R + r * cos(ϕ)) * cos(θ)
            x₂ = (R + r * cos(ϕ)) * sin(θ)
            x₃ = r * sin(ϕ)
            array[i, j] = ℝ³(x₁, x₂, x₃)
        end
    end
    map(x -> x + gettranslation(q), rotate(array, getrotation(q)))
end

"""
    constructsphere(q, radius)

Construct a sphere with the given configuration `q` and `radius`.
"""
function constructsphere(q::Biquaternion,
                         radius::Real;
                         segments::Int = 36)
    array = Array{ℝ³,2}(undef, segments, segments)
    lspace = collect(range(float(pi), stop = float(-pi), length = segments))
    for i in 1:segments
        for j in 1:segments
            ϕ = lspace[i]
            θ = lspace[j] / 2
            array[j, i] = ℝ³(Cartesian(Geographic(radius, ϕ, θ)))
        end
    end
    map(x -> x + gettranslation(q), rotate(array, getrotation(q)))
end

"""
    constructcylinder(q, height, radius)

Construct a cylinder with the given configuration `q`, `height` and `radius`.
"""
function constructcylinder(q::Biquaternion,
                           height::Float64,
                           radius::Float64;
                           segments::Int = 36)
    array = Array{ℝ³,2}(undef, segments, segments)
    lspace = collect(range(float(-pi), stop = float(pi), length = segments))
    for i in 1:segments
        for j in 1:segments
            α = lspace[j]
            x₁ = radius * cos(α)
            x₂ = radius * sin(α)
            x₃ = (i - 1) * height / (segments - 1)
            array[i, j] = ℝ³(x₁, x₂, x₃)
        end
    end
    map(x -> x + gettranslation(q), rotate(array, getrotation(q)))
end

"""
    constructbox(p, q, s)

Construct a box with the given translation `p`, rotation `q` and scale `s`.
"""
function constructbox(p::ℝ³, q::Quaternion, s::ℝ³)
    x, y, z = vec(s) ./ 2
    points = [ℝ³(-x, -y, -z) ℝ³(x, -y, -z) ℝ³(-x, -y, -z) ℝ³(x, -y, -z);
              ℝ³(-x, -y, +z) ℝ³(x, -y, +z) ℝ³(-x, -y, +z) ℝ³(x, -y, +z);
              ℝ³(-x, +y, +z) ℝ³(x, +y, +z) ℝ³(-x, +y, +z) ℝ³(x, +y, +z);
              ℝ³(-x, +y, -z) ℝ³(x, +y, -z) ℝ³(-x, +y, -z) ℝ³(x, +y, -z)]
    map(x -> x + p, rotate(points, q))
end

"""
    getactions(point,
               s2tos3map,
               s2tos2map;
               [top, [bottom, [s3rotation, [segments [, scale]]]]])

Take the integral along the Hopf fiber for getting a linear action space
under stereographic projection, with the given 'point' in the base space,
a map from the base space to the total space `s2tos3map`, a conformal map that
takes the base space to itself `s2tos2map`, `top` and `bottom` U(1) actions,
S³ rotation `s3rotation`, the number of `segments` and `scale`.
"""
function getactions(point::S²,
                    s2tos3map,
                    s2tos2map;
                    top::S¹ = U1(-pi),
                    bottom::S¹ = U1(pi),
                    s3rotation::S³ = Quaternion(1, 0, 0, 0),
                    segments::Int = 30,
                    scale::Real = 1.0)
    segments2 = 1000
    array = Array{Float64,2}(undef, segments2, 2)
    minaction = -π
    maxaction = π
    lspace = collect(range(minaction, stop = maxaction, length = segments2))
    array[:, begin] = lspace
    accumulator = 0
    getp(point, α, rotation, scale) = compressedλmap(rotate(S¹action(s2tos3map(s2tos2map(point)), U1(α)), rotation)) * scale
    α₁ = lspace[begin]
    p₁ = getp(point, α₁, s3rotation, scale)
    for i in 1:segments2
        α₂ = lspace[i]
        p₂ = getp(point, α₂, s3rotation, scale)
        accumulator += norm(p₂ - p₁)
        array[i, end] = accumulator
        α₁ = α₂
        p₁ = p₂
    end

    # decide based on cumulative length and then pick the corresponding actions from the array
    totalspan = maxaction - minaction
    normalizedtop = angle(top) / totalspan
    normalizedbottom = angle(bottom) / totalspan
    normalizedsegmentlength = (normalizedbottom - normalizedtop) / (segments - 1)
    actions = Array{Float64,1}(undef, segments)
    for i in 1:segments
        target = (normalizedtop + normalizedsegmentlength * (i - 1)) * accumulator
        for j in 1:segments2
            if array[j, end] ≥ target
                actions[i] = array[j, begin]
                break
            end
        end
    end
    actions
end

"""
    transformg(p, g₁, g₂, segments)

Transform G-acions by finding and applying affine parameters with the given point `p`
in the principal bundle, beginning gauge potential `g₁`, ending gauge potential `g₂`
and the number of segments.
"""
function transformg(p::ComplexPlane, g₁::U1, g₂::U1, segments::Int)
    intervals = 1000
    array = Array{Float64,2}(undef, intervals, 2)
    minaction, maxaction = 0, 2π
    lspace = collect(range(minaction, stop = maxaction, length = intervals))
    array[:, begin] = lspace
    accumulator = 0
    α₁ = lspace[begin]
    z = λ⁻¹map(πmap(Quaternion(p)))
    p₁ = compressedλmap(S¹action(z, U1(α₁)))
    array[begin, end] = 0
    for i in 2:intervals
        α₂ = lspace[i]
        p₂ = compressedλmap(S¹action(p, U1(α₂)))
        accumulator += norm(p₂ - p₁)
        array[i, end] = accumulator
        α₁ = α₂
        p₁ = p₂
    end

    # decide based on cumulative length and then pick the corresponding actions from the array
    totalspan = maxaction - minaction
    n₁ = angle(g₁) / totalspan
    n₂ = angle(g₂) / totalspan
    normalizedsegmentlength = (n₂ - n₁) / (segments - 1)
    actions = Array{U1,1}(undef, segments)
    for i in 1:segments
        target = (n₁ + normalizedsegmentlength * (i - 1)) * accumulator
        for j in 1:intervals
            if array[j, end] ≥ target
                actions[i] = U1(array[j, begin])
                break
            end
        end
    end
    actions
end

"""
    constructwhirl(points, gauge1, gauge2, configuration, segments)

Construct a whirl with the given `points`, `gauge1`, `gauge2`, `configuration` and `segments`.
"""
function constructwhirl(points::Array{ComplexPlane,1},
                        gauge1::Array{U1,1},
                        gauge2::Array{U1,1},
                        configuration::Biquaternion,
                        segments::Int)
    array = Array{ℝ³,2}(undef, segments, length(points))
    for (j, p) in enumerate(points)
        g₁ = gauge1[j]
        g₂ = gauge2[j]
        g = transformg(p, g₁, g₂, segments)
        for (i, α) in enumerate(g)
            array[i, j] = compressedλmap(S¹action(p, α))
        end
    end
    applyconfig(array, configuration)
end

"""
    constructfiber(point,
                   s2tos3map,
                   s2tos2map;
                   [radius, [top, [bottom, [s3rotation, [config, [segments [, scale]]]]]]])

Construct a Hopf fiber with the given `point` in the base space, , a map from the base space to
the total space `s2tos3map`, a conformal map that takes the base space to itself `s2tos2map`,
fiber sectional `radius`, `top` and `bottom` U(1) actions, S³ rotation `s3rotation`,
configuration `config`, the number of `segments` and `scale`.
"""
function constructfiber(point::S²,
                        s2tos3map,
                        s2tos2map;
                        radius::Float64 = 0.02,
                        top::S¹ = U1(-pi),
                        bottom::S¹ = U1(pi),
                        s3rotation::S³ = Quaternion(1, 0, 0, 0),
                        config::Biquaternion = Biquaternion(ℝ³(0, 0, 0)),
                        segments::Int = 30,
                        scale::Real = 1.0)
    segments2 = 10
    α₁ = U1(min(angle(top), angle(bottom)))
    α₂ = U1(max(angle(top), angle(bottom)))
    circle = [ℝ³(cos(α), sin(α), 0) for α in range(0, stop = 2π, length = segments2)] .* radius
    actions = getactions(point, s2tos3map, s2tos2map, top = α₁, bottom = α₂, s3rotation = s3rotation, segments = segments, scale = scale)
    getp(point, α, rotation, scale) = compressedλmap(rotate(S¹action(s2tos3map(s2tos2map(point)), U1(α)), rotation)) * scale
    array = Array{ℝ³,2}(undef, segments, segments2)
    n = ℝ³(0, 0, 1)
    for (i, α) in enumerate(actions)
        p₁ = getp(point, α, s3rotation, scale)
        p₂ = getp(point, α + α * 1e-2, s3rotation, scale)
        w = normalize(p₂ - p₁)
        q = Biquaternion(getrotation(n, w), p₁)
        array[i, :] = applyconfig(circle, q)
    end
    applyconfig(array, config)
end

"""
    constructtwospinor(point, gauge1, gauge2, configuration, radius, segments1, segments2)

Construct a 2-spinor with the given `point` in the principal bundle, `gauge1`, `gauge2`,
`configuration`, `radius`, `segments1` and `segments2`.
"""
function constructtwospinor(point::ComplexPlane, gauge1::U1, gauge2::U1, configuration::Biquaternion,
                            radius::Float64, segments1::Int, segments2::Int)
    angles = angle.([gauge1; gauge2])
    g₁, g₂ = U1(min(angles...)), U1(max(angles...))
    circle = [ℝ³(cos(α), sin(α), 0) for α in range(0, stop = 2π, length = segments2)] .* radius
    actions = transformg(point, g₁, g₂, segments1)
    getp(p::ComplexPlane, g::U1) = compressedλmap(S¹action(p, g))
    array = Array{ℝ³,2}(undef, segments1, segments2)
    n = ℝ³(0, 0, 1)
    ϵ = 1e-2
    for (i, g) in enumerate(actions)
        p₁ = getp(point, g)
        p₂ = getp(point, U1(angle(g) + ϵ))
        w = normalize(p₂ - p₁)
        q = Biquaternion(getrotation(n, w), p₁)
        array[i, :] = applyconfig(circle, q)
    end
    applyconfig(array, configuration)
end

"""
    constructframe(section, configuration, segments)

Construct a frame that is like an S³ cross-section under stereographic projection with the
given `section`, `configuration` and `segments`.
"""
function constructframe(section::Any,
                        configuration::Biquaternion = Biquaternion(ℝ³(0, 0, 0)),
                        segments::Int = 36)
    array = Array{ℝ³,2}(undef, segments, segments)
    lspaceϕ = collect(range(float(-pi), stop = float(pi), length = segments))
    lspaceθ = collect(range(float(pi / 2), stop = float(-pi / 2), length = segments))
    for (i, θ) in enumerate(lspaceθ)
        for (j, ϕ) in enumerate(lspaceϕ)
            p = Geographic(1, ϕ, θ)
            array[i, j] = compressedλmap(section(p))
        end
    end
    applyconfig(array, configuration)
end

getplane(p::Array{Float64}, q::Quaternion, s::Array{Float64}) = begin
    plane = Array{Float64,3}(undef, 2, 2, 3)
    plane[1, 1, :] = vec(rotate(ℝ³([0.0; -1.0; 1.0]), q)) .* s + p
    plane[1, 2, :] = vec(rotate(ℝ³([0.0; 1.0; 1.0]), q)) .* s + p
    plane[2, 1, :] = vec(rotate(ℝ³([0.0; -1.0; -1.0]), q)) .* s + p
    plane[2, 2, :] = vec(rotate(ℝ³([0.0; 1.0; -1.0]), q)) .* s + p
    plane
end

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
import Base.Dates: guess

"""
    guess(start::ZonedDateTime, finish::ZonedDateTime, step) -> Integer

Given a start and end date, indicates how many steps/periods are between them. Defining this
function allows `StepRange`s to be defined for `ZonedDateTime`s.
"""
function guess(start::ZonedDateTime, finish::ZonedDateTime, step)
    guess(start.utc_datetime, finish.utc_datetime, step)
end


text:
module Imports

using PyCall

global const mixture = PyNULL()
global const ensemble = PyNULL()
global const model_selection = PyNULL()
global const gmeans = PyNULL()

export mixture, gmeans, ensemble

function __init__()
  copy!(mixture, pyimport("sklearn.mixture"))
  copy!(ensemble, pyimport("sklearn.ensemble"))
  copy!(model_selection, pyimport("sklearn.model_selection"))
  copy!(gmeans, pyimport("pyclustering.cluster.gmeans"))
end
end


text:
"""
    FixedHomotopy(homotopy, t)

Fix the homotopy at `t`, i.e. consider ``H(x,t)`` as ``H(x)``
"""
struct FixedHomotopy{T<:Number, H<:AbstractHomotopy{T}, S<:Number} <: AbstractPolySystem{T}
    homotopy::H
    t::S
end
function FixedHomotopy(hom::H, t::S) where {T<:Number, H<:AbstractHomotopy{T}, S<:Number}
    FixedHomotopy{T, H, S}(hom, t)
end

@inline evaluate(H::FixedHomotopy, x) = evaluate(H.homotopy, x, H.t)
(H::FixedHomotopy)(x) = evaluate(H,x)

function Base.promote_rule(::Type{FixedHomotopy{T1,H1,S1}}, ::Type{FixedHomotopy{T2, H2, S2}}) where {T1, T2, H1, H2, S1, S2}
    FixedHomotopy{promote_type(T1,T2), promote_type(H1,H2), promote_type(S1,S2)}
end

function Base.promote_rule(::Type{FixedHomotopy{T1,H1,S1}}, ::Type{T2}) where {T1, T2, H1, S1}
    FixedHomotopy{promote_type(T1,T2), promote_type(H1,T2), promote_type(S1,T2)}
end

function Base.convert(::Type{FixedHomotopy{T1,H1,S1}}, hom::FixedHomotopy) where{T1, H1, S1}
    FixedHomotopy(convert(H1, hom.homotopy), convert(S1, hom.t))
end

function Base.show(io::IO, H::FixedHomotopy)
    print(io, typeof(H), ":\n")
    println(io, "* homotopy:")
    println(io, H.homotopy)
    println(io, "* t: ", H.t)
end
function differentiate(H::FixedHomotopy)
    dh = differentiate(H.homotopy)
    (x) -> dh(x, H.t)
end

substitute(H::FixedHomotopy, pair) = FixedHomotopy(substitute(H.homotopy, pair), H.t)
homogenize(H::FixedHomotopy) = FixedHomotopy(homogenize(H.homotopy), H.t)
homogenized(H::FixedHomotopy) = homogenized(H.homotopy)
ishomogenous(H::FixedHomotopy) = ishomogenous(H.homotopy)
nvariables(H::FixedHomotopy) = nvariables(H.homotopy)
degrees(H::FixedHomotopy) = degrees(H.homotopy)

Base.length(H::FixedHomotopy) = length(H.homotopy)
coefftype(::FixedHomotopy{T, H, S}) where {T, H, S} = T