text:
# fast_LP_test.jl
# Interface for running fast_LP methods on metropolis instead of Manhattan
# This will call fast_LP with a modified version of Manhattan

# Default parameters: ignore since everything is loaded from JSON
PROB = 0.7
MODEL = "relaxed_smooth"
LAST_SMOOTH = false
MAX_ROUNDS = 50
MIN_RIDES = 1
MAX_NUM_PATHS_PER_OD = 3
TURN_COST = 2.0
TURN_COST_AS_VARIABLE = false
DELTA_BOUND = [0.6,0.5,0.4]
START_SIMPLE = false
DYNAMIC_CONSTRAINTS = true 			# true for dynamic constraints, false otw
GLOBAL_CONSTRAINT_UPDATE = false 	# true if we remove worst paths globally, false otw
SAMPLE_SIZE = 1000 					# starting number of constraints
NUM_OD_ADDED = 1000 				# number of (O,D) pairs to add
NUM_OD_REMOVED = 1000				# number of (O,D) pairs to remove
UPDATE_EVERY_N_ITERATIONS = 1 		# add number of (O,D) above every $N iterations
ADD_IF_ABOVE_PERCENTILE = 0.7		# only add rides that have an error above this percentile
REMOVE_IF_BELOW_PERCENTILE = 0.3	# only remove rides that have an error below this percentile

function fast_LP_test(;
	model_type::AbstractString=MODEL,
	last_smooth::Bool=LAST_SMOOTH,
	prob::Float64=PROB,
	max_rounds::Int=MAX_ROUNDS,
	min_rides::Int=MIN_RIDES,
	turnCost::Float64=TURN_COST,
	turnCostAsVariable::Bool=TURN_COST_AS_VARIABLE,
	delta_bound::Union{Array{Float64,1},Array{Any,1}}=DELTA_BOUND,
	maxNumPathsPerOD::Int=MAX_NUM_PATHS_PER_OD,
	dynamicConstraints::Bool=DYNAMIC_CONSTRAINTS,
	numPairsToAdd::Int = NUM_OD_ADDED,
	numPairsToRemove::Int = NUM_OD_REMOVED,
	iterationMultiple::Int = UPDATE_EVERY_N_ITERATIONS,
	addOnlyIfAbovePercentile::Float64 = ADD_IF_ABOVE_PERCENTILE,
	removeOnlyIfBelowPercentile::Float64 = REMOVE_IF_BELOW_PERCENTILE,
	globalConstraintUpdate::Bool = GLOBAL_CONSTRAINT_UPDATE,
	sample_size::Int=SAMPLE_SIZE,
	startWithSimpleLP::Bool=START_SIMPLE)

	# This code loads the data
	adjacencyList = load("Inputs/input-graph.jld", "adjList")
	graph = DiGraph(length(adjacencyList))
	for i = 1:nv(graph), j in adjacencyList[i]
		add_edge!(graph, i, j)
	end
	real_times = load("Inputs/input-realTimes.jld", "realTimes")
	minRoadTime = load("Inputs/input-speedLimits.jld", "speedLimits")
	meanRoadTime = load("Inputs/input-meanTimes.jld", "meanTimes")
	distances = load("Inputs/input-distances.jld", "distances")
	coordinates = load("Inputs/input-positions.jld", "coordinates")
	positions = Coordinates[]
	for i = 1:nv(graph)
		push!(positions, Coordinates(coordinates[i,1], coordinates[i,2]))
	end

	travel_times = load("Inputs/input-travelTimes-$(prob).jld", "travelTimes")
	num_rides = load("Inputs/input-numRides-$(prob).jld", "numRides")

	metropolis_manhattan = Manhattan(emptyType = true)
	metropolis_manhattan.network = graph
	metropolis_manhattan.positions = positions
	metropolis_manhattan.distances = distances
	metropolis_manhattan.roadTime = minRoadTime

	status, new_times, directory = fast_LP(metropolis_manhattan, travel_times, num_rides, zeros(nv(graph),nv(graph)), minRoadTime, model_type=model_type, last_smooth = last_smooth, max_rounds=max_rounds, min_rides=min_rides, year=0, sample_size=sample_size, turnCost=turnCost, turnCostAsVariable=turnCostAsVariable, delta_bound=delta_bound, maxNumPathsPerOD=maxNumPathsPerOD, computeFinalSP=true,	startWithSimpleLP=startWithSimpleLP, randomConstraints=false, dynamicConstraints=dynamicConstraints, numPairsToAdd=numPairsToAdd, numPairsToRemove=numPairsToRemove, iterationMultiple = iterationMultiple,	addOnlyIfAbovePercentile = addOnlyIfAbovePercentile, removeOnlyIfBelowPercentile = removeOnlyIfBelowPercentile, globalConstraintUpdate=globalConstraintUpdate, warmstart=false, metropolis=true, real_TOD_metropolis = real_times, real_tij_metropolis = meanRoadTime, prob=prob, last_smooth=last_smooth, warmstart=false, localized_testing=false)

	return status, new_times, directory
end

# Load dependencies
include("metropolis_intro.jl")

t=time()
# Load parameters
function load_parameters()
	if length(ARGS) == 0
		println("Error: shell script failed to pass in arguments")
		return nothing
	else
		JSONfileName = ARGS[1]
		params = JSON.parsefile(JSONfileName)
		return JSONfileName, params
	end
end
JSONfile, params = load_parameters()

# run code
status, new_times, directory = @time fast_LP_test(model_type=params["MODEL"], last_smooth=params["LAST_SMOOTH"], prob=params["PROB"], max_rounds=params["MAX_ROUNDS"], min_rides=params["MIN_RIDES"], turnCost=params["TURN_COST"], turnCostAsVariable=params["TURN_COST_AS_VARIABLE"], delta_bound=params["DELTA_BOUND"], maxNumPathsPerOD=params["MAX_NUM_PATHS_PER_OD"], dynamicConstraints=params["DYNAMIC_CONSTRAINTS"], numPairsToAdd = params["NUM_OD_ADDED"], numPairsToRemove = params["NUM_OD_REMOVED"], iterationMultiple = params["UPDATE_EVERY_N_ITERATIONS"], addOnlyIfAbovePercentile = params["ADD_IF_ABOVE_PERCENTILE"], removeOnlyIfBelowPercentile = params["REMOVE_IF_BELOW_PERCENTILE"], globalConstraintUpdate = params["GLOBAL_CONSTRAINT_UPDATE"], sample_size=params["SAMPLE_SIZE"], startWithSimpleLP=params["START_SIMPLE"])

# Move log files
f=open(string(directory, "timestats.csv"), "a")
write(f, string(time()-t, '\n'))
close(f)
mv(JSONfile, string(directory, "parameters.json"), remove_destination=true)

max_stars_repo_name: sebmart/NYC-timing

max_issues_repo_name: sebmart/NYC-timing

max_forks_repo_name: sebmart/NYC-timing


text:
"""
   HeistEvents is a object that creates a `HiestEvent` and loops over the events.
   It can be used as an event loop. You create it by passing in a `HeistFiles`
   object that represents the files to run over and an integer `nMax` that is the maximum
   number of events to process. Set to zero if you want to process everything.

# Examples

```julia
for anEvent ∈ HeistEvents(hf, 200)
   println( @cxx anEvent.galleryEvent->eventEntry() )
   tracks = getRecord(anEvent, tas)
   ghosts = getRecord(anEvent, gds)
   # ... do something with the objects ...
end
```
"""
struct HeistEvents
   hf::HeistFiles
   nDesired::Int
end

# See https://docs.julialang.org/en/stable/stdlib/collections/#lib-collections-iteration-1 and
# see https://docs.julialang.org/en/stable/manual/interfaces/#man-interface-iteration-1
# for how this works as a for loop. That is,
#   for i = I   # or  "for i in I"
#     ...
#   end
#
#   is translated into
#
#  state = start(I)
#  while !done(I, state)
#      (i, state) = next(I, state)
#      ...
#  end
#
#  So need to define start, next and done. eltype improves efficiency
#
Base.start(::HeistEvents) = nothing  # Don't do anything at the start --
                                     # next will make the HeistEvent -- this allows us to see the first event

Base.next(hes::HeistEvents, state::Void)       = (he = HeistEvent(hes.hf) ; (he, he) )  # Multiple dispatch in action!
Base.next(hes::HeistEvents, state::HeistEvent) = (nextEvent(state), state)

Base.done(hes::HeistEvents, state::Void)       = false
Base.done(hes::HeistEvents, state::HeistEvent) = atEnd(state) ||
                                                 (hes.nDesired > 0 && state.nSeen ≥ hes.nDesired)

Base.eltype(::Type{HeistEvents}) = HeistEvent

ext: jl


text:
using Documenter

makedocs(
   # modules = [NanoDates],
    sitename = "NanoDates.jl",
    authors = "Jeffrey Sarnoff",
    source = "src",
    clean = false,
    strict = !("strict=false" in ARGS),
    doctest = ("doctest=only" in ARGS) ? :only : true,
    format = Documenter.HTML(
        # Use clean URLs, unless built as a "local" build
        prettyurls = !("local" in ARGS),
        highlights = ["yaml"],
        ansicolor = true,
    ),
    pages = [
        "Home" => "index.md",
        "Advantages" => "appropriate/advantages.md",
        "Basic Use" => Any[
            "Construction" => "use/construction.md",
            "Specification" => "use/specify.md",
            "Accessors" => "use/accessors.md",
            "Timestamps" => "use/timestamps.md",
        ],
        "Helpful Extras" => Any[
            "Period Compounding" => "use/compoundperiod.md",
            "Conviences" => "use/convenient.md",
            "Strings" => "use/intostring.md",
        ],
        "tempus fugit" => Any[
             "Outside Time" => "betterinformed/the_representation.md",
             "Inside  Time" =>  "betterinformed/their_dual_natures.md",
        ]
    ]
)

#=
Deploy docs to Github pages.
=#
Documenter.deploydocs(
    branch = "gh-pages",
    target = "build",
    deps = nothing,
    make = nothing,
    repo = "github.com/JeffreySarnoff/NanoDates.jl.git",
)


text:
# # Spectrally accurate Nyström method for the Helmholtz equation in 2D
# ## Carlos Pérez Arancibia (cperez@mat.uc.cl)

using WaveProp, WaveProp.Geometry, WaveProp.Integration, WaveProp.Nystrom, WaveProp.Mesh
using LinearAlgebra
using Plots

pyplot()

#=

## Validation

=#
k   = 20
pde = Helmholtz(dim=2,k=k)
# geo = Kite()
δ = 1e-4
ref = ReferenceLine() # the [0,1] segment
bnd = ParametricEntity(ref) do u
     x = cos(2π*u[1])
     y = (δ+x^2)*sin(2π*u[1])
     SVector(x,y)
end
geo = ParametricBody(boundary=bnd)
G   = SingleLayerKernel(pde)
xₛ  = (0.5,0)
u   = (x) -> G(xₛ,x)
xₒ  = (3,-1)
subfig = plot(geo,label="Γ")
plot!(subfig,[xₛ],m=:o,ms=10,label="source location")
plot!(subfig,[xₒ],m=:o,ms=10,label="test location")
η     = k
Ω,M   = meshgen(geo)
Γ     = boundary(Ω)
niter = 7
er    = zeros(niter)
nn    = 40*collect(1:niter)

for i in 1:niter
     quad_rule = TrapezoidalP(nn[i])
     mesh      = NystromMesh(M,Γ;quad_rule)
     S = SingleLayerOperator(pde,mesh)
     D = DoubleLayerOperator(pde,mesh)
     Smk = Nystrom.assemble_mk(S)
     Dmk = Nystrom.assemble_mk(D)
     L     = I/2 + Dmk - im*η*Smk # create a dense matrix
     rhs = γ₀(u,mesh)
     σ   = L\rhs
     𝒮 = SingleLayerPotential(pde,mesh)
     𝒟 = DoubleLayerPotential(pde,mesh)
     uₙ = (x) -> 𝒟[σ](x) - im*η*𝒮[σ](x)
     xx = (5,5)
     er[i] = abs(u(xx) - uₙ(xx))/abs(u(xx))
end
fig = plot(nn,er,m=:x,yscale=:log10,ylabel="log₁₀(error)",xlabel="n",label="(uₙ - u)/u")

niter = 9
er_dim    = []
er_ldim    = []
nn    = []
quad_rule = GaussLegendre(6)
for i in 1:niter
     Ω,M   = meshgen(geo;gridsize=(1/2)^i)
     Γ     = boundary(Ω)
     mesh      = NystromMesh(M,Γ;quad_rule)
     S = SingleLayerOperator(pde,mesh)
     D = DoubleLayerOperator(pde,mesh)
     Sdim = Nystrom.assemble_dim(S)
     Ddim = Nystrom.assemble_dim(D)
     Sldim = Nystrom.assemble_ldim(S)
     Dldim = Nystrom.assemble_ldim(D)
     Ldim     = I/2 + Ddim - im*η*Sdim # create a dense matrix
     Lldim     = I/2 + Dldim - im*η*Sldim # create a dense matrix
     rhs = γ₀(u,mesh)
     σdim    = Ldim\rhs
     σldim   = Lldim\rhs
     𝒮 = SingleLayerPotential(pde,mesh)
     𝒟 = DoubleLayerPotential(pde,mesh)
     xx = (5,5)
     uₙ = (x) -> 𝒟[σdim](x) - im*η*𝒮[σdim](x)
     push!(er_dim,abs(u(xx) - uₙ(xx))/abs(u(xx)))
     uₙ = (x) -> 𝒟[σldim](x) - im*η*𝒮[σldim](x)
     push!(er_ldim,abs(u(xx) - uₙ(xx))/abs(u(xx)))
     push!(nn,length(mesh))
end
fig = plot(nn,er_dim,m=:x,xscale=:log2,yscale=:log10,ylabel="log₁₀(error)",xlabel="n",label="gdim")
plot!(fig,nn,er_ldim,m=:x,xscale=:log2,yscale=:log10,ylabel="log₁₀(error)",xlabel="n",label="ldim")

plot(fig,subfig,layout=(2,1))

#===============================================
We consider a kite-shaped curve given by:
```math
     \Gamma = \{(\cos(t)+0.65\cos(2t)-0.65, 1.5\sin(t))\in\mathbb R^2,t\in (0,2\pi)\}
```
and we let $\Omega\subset\mathbb{R}^2$ be the (open) domain enclosed by
$\Gamma$, i.e., $\Gamma=\partial\Omega$.

We can construct this **parametric curve** and plot it through the `Geometry` module:
===============================================#

geo = Kite()
plot(geo)

#===============================================
The boundary value problem that we consider corresponds to an incident planewave
$u^{\rm inc}(x_1,x_2) = \exp(ik(x_1\cos\theta+x_2\sin\theta))$ where
$\theta$ is the angle of incidence and $k>0$ is the wavenumber, that impinges on
$\Gamma$ giving rise to a scattered field $u$ which satisfies:
```math
\begin{aligned}
\Delta u+k^2 u &=0\quad\text{in}\quad\mathbb R^2\setminus\overline\Omega,\\
u &=-u^{\rm inc}\quad\text{on}\quad\Gamma,\\
\lim_{|x|\to\infty}|x|^{1/2}\left\{\frac{\partial u}{\partial |x|}-iku\right\} &= 0
\end{aligned}
```

We proceed to define the partial differential equation (**PDE**) that we wish so solve,
as well as the incident wavefield
===============================================#
θ   = 0.
d   = (cos(θ),sin(θ))
uᵢ  = (x) -> exp.(im*k*(x⋅d))
∂uᵢ(x,ν) = im*k*(ν⋅d)*uᵢ(x)

#=
In order to solve the exterior boundary value problem, we resort to the combined field integral equation formulation, whereby we seek the scattered field in the form:

$u(x) = (\mathcal D-i\eta S)[\varphi](x),\quad x\in\mathbb R^2\setminus\Gamma,$

where $\mathcal D$ and $\mathcal S$ are the double- and single-layer potentials defined as:

$(\mathcal D\varphi)(x) = \int_{\Gamma}\frac{\partial G(x,y)}{\partial n(y)}\varphi(y)ds(y)\quad\text{and}\quad (\mathcal S\varphi)(x) = \int_{\Gamma}G(x,y)\varphi(y)ds(y).$

Here $G(x,y) = \frac{i}4 H_0^{(1)}(k|x-y|)$ is the free-space Green's function for Helmholtz equation and $\varphi:\Gamma\to\mathbb C$ is an unknown density function that will be determined by solving an integral equation on curve $\Gamma$.

The constant $\eta>0$ is intruduced so as to ensure uniqueness of solutions of the boundary integral equation. We now define $\eta$ and the tuple _coeff_ that will be used below to construct the combined field potential and operator:

Imposing the Dirichlet boundary condition on $\Gamma$, we arrive at the following second-kind integral equation:

$\left(\frac{I}{2}+K-i\eta S\right)\varphi = -u^{\rm inc}\quad\text{on}\quad\Gamma,$

where operators above are the so-called double- and single-layer integral operators defined as:

$(K\varphi)(x) = \int_{\Gamma}\frac{\partial G(x,y)}{\partial n(y)}\varphi(y)ds(y)\quad\text{and}\quad(S\varphi)(x) = \int_{\Gamma}G(x,y)\varphi(y)ds(y),\qquad x\in\Gamma.$

where we have used the jump relations:

$$\lim_{\epsilon\to 0\pm}(\mathcal D\varphi)(x+\epsilon\nu(x))=\pm\frac{\varphi(x)}{2}+(K\varphi)(x)$$

and

$$\lim_{\epsilon\to 0\pm}(\mathcal S\psi)(x+\epsilon\nu(x))=(S\psi)(x)$$

for $x\in\Gamma.$

We will approximate the integral operators $S$ and $K$ using a Nyström method, for which we must specify (i) a
mesh of $\Gamma$ and (ii) a quadrature rule to be used in each element of the mesh. This is done as follows
=#

n     = 10*k
quad_rule = TrapezoidalP(n)
mesh  = NystromMesh(M,Γ;quad_rule)

# Now that an appropriate mesh is available, we can create the integral operators:
S = SingleLayerOperator(pde,mesh)
D = DoubleLayerOperator(pde,mesh)

#=====
The `S` and `D` variables are `IntegralOperator`s, which are in fact a light wrapper around the
information needed to defined these operators mathematically. They are subtypes of `AbstractMatrix`,
and behave as such. For instance, we can now *solve* for the unknown density $\sigma$ as follows:
=====#

η   = k
L   = I/2 + D - im*η*S # create a dense matrix
rhs = -γ₀(uᵢ,mesh)
σ   = L\rhs

#=

There are two possible issues with the solution method we just employed: **efficiency** and **accuracy**. We will return to those in a bit, but to get to the end of things, let us now visualize the solution obtained:

=#

𝒮 = SingleLayerPotential(pde,mesh)
𝒟 = DoubleLayerPotential(pde,mesh)
uₛ = (x) -> 𝒟[σ](x) - im*η*𝒮[σ](x)
u  = x -> uᵢ(x) + uₛ(x) # the total field
xgrid = ygrid = -5:(2π/k/10):5
U  = [(x,y) ∈ mesh ? NaN+im*NaN : u((x,y)) for y in ygrid, x in xgrid]
fig1 = heatmap(xgrid,ygrid,real.(U),clims=(-2,2))
plot!(fig1,mesh,lw=6,lc=:white)

#=

Because the integral operators we just used are (weakly) singular, special care is needed to compute the
entries of `S` and `D` if high-order convergence is to be obtained. There are several techniques for doing
so, and deciding which to use is highly application dependent. For two-dimensional smooth obstacles admiting a global parametrization, the most accurate method is probably due to Martensen and Kussmaul. To *assemble* the operators `S` and `K` usign this method, we can simply write

=#

Smk = Nystrom.assemble_mk(S)
Dmk = Nystrom.assemble_mk(D)

# The rest follows exactly as before:

L   = I/2 + Dmk - im*η*Smk # create a dense matrix
σ   = L\rhs
poly = copy(qnodes(mesh))
push!(poly,qnodes(mesh)[1])
U    = [ (x,y) in mesh ? NaN+im*NaN : u((x,y)) for y in ygrid, x in xgrid]
fig2 = heatmap(xgrid,ygrid,real.(U),clims=(-2,2))
plot!(fig2,mesh,lw=6,lc=:white)

# The difference between the naive method and the Martensen-Kusmall method is
# evident. In particular the shadow:
plot(fig1,fig2,layout=(1,2),colorbar=:false,ylims=(-5,5))

# Adding more scatterers is simple:
geo   = [Kite(radius=0.1,center=(xc,yc)) for yc in -1:1:1, xc in -1:0.5:1] |> vec
Ω,M   = meshgen(geo)
Γ     = boundary(Ω)
mesh  = NystromMesh(M,Γ;quad_rule)
S     = SingleLayerOperator(pde,mesh)
D     = DoubleLayerOperator(pde,mesh)
Smk   = Nystrom.assemble_mk(S)
Dmk   = Nystrom.assemble_mk(D)
L     = I/2 + Dmk - im*η*Smk # create a dense matrix
rhs   = -γ₀(uᵢ,mesh)
σ     = L\rhs
## now output the solution.
𝒮     = SingleLayerPotential(pde,mesh)
𝒟     = DoubleLayerPotential(pde,mesh)
uₛ = (x) -> 𝒟[σ](x) - im*η*𝒮[σ](x)
u  = x -> uᵢ(x) .+ uₛ(x) # the total field
xgrid = ygrid = -2:(2π/k/10):2
U  = [(x,y) in mesh ? NaN+im*NaN : u((x,y)) for y in ygrid, x in xgrid]
# U  = [u(SVector(x,y)) for y in ygrid, x in xgrid]
heatmap(xgrid,ygrid,real.(U),clims=(-2,2))
plot!(mesh,lw=2,lc=:black)
plot!(size=(2000,2000))

#=

## Adding your own shapes

Not everybody wants to solve scattering by an arrays of kites, so lets create our own shape and solve scattering problem using that.
=#
## define the shape throught its boundary
ref = ReferenceLine() # the [0,1] segment
bnd = ParametricEntity(ref) do u
     r = 1 + 0.25*cos(10π*u[1])
     SVector(r*cos(2π*u[1]),r*sin(2π*u[1]))
end
geo = ParametricBody(boundary=bnd)
## solve the scattering problem
Ω,M   = meshgen(geo)
Γ     = boundary(Ω)
mesh  = NystromMesh(M,Γ;quad_rule=TrapezoidalP(1000))
S     = SingleLayerOperator(pde,mesh)
D     = DoubleLayerOperator(pde,mesh)
Smk   = Nystrom.assemble_mk(S)
Dmk   = Nystrom.assemble_mk(D)
L     = I/2 + Dmk - im*η*Smk # create a dense matrix
rhs   = -γ₀(uᵢ,mesh)
σ     = L\rhs
## now output the solution.
𝒮     = SingleLayerPotential(pde,mesh)
𝒟     = DoubleLayerPotential(pde,mesh)
uₛ = (x) -> 𝒟[σ](x) - im*η*𝒮[σ](x)
u  = x -> uᵢ(x) .+ uₛ(x) # the total field
xgrid = ygrid = -3:(2π/k/10):3
U  = [(x,y) ∈ mesh ? NaN+im*NaN : u(SVector(x,y)) for y in ygrid, x in xgrid]
heatmap(xgrid,ygrid,real.(U),clims=(-2,2))
plot!(mesh,lw=4,lc=:white)
plot!(size=(2000,2000))

# # Composite quadrature rules
k = 40
pde = Helmholtz(k=k,dim=2)
Ω,M   = meshgen(geo;gridsize=0.01)
Γ     = boundary(Ω)
mesh  = NystromMesh(M,Γ;quad_rule=TrapezoidalP(20))
S     = SingleLayerOperator(pde,mesh)
D     = DoubleLayerOperator(pde,mesh)
Smk   = Nystrom.assemble_dim(S)
Dmk   = Nystrom.assemble_dim(D)
L     = I/2 + Dmk - im*k*Smk # create a dense matrix
rhs   = -γ₀(uᵢ,mesh)
σ     = L\rhs
## now output the solution.
𝒮     = SingleLayerPotential(pde,mesh)
𝒟     = DoubleLayerPotential(pde,mesh)
uₛ = (x) -> 𝒟[σ](x) - im*k*𝒮[σ](x)
u  = x -> uᵢ(x) .+ uₛ(x) # the total field
xgrid = ygrid = -3:(2π/k/10):3
U  = [(x,y) ∈ mesh ? NaN+im*NaN : u(SVector(x,y)) for y in ygrid, x in xgrid]
heatmap(xgrid,ygrid,real.(U),clims=(-2,2))
plot!(mesh,lw=1,lc=:white)

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
struct Scale{T}
    value :: T
end

(𝒮::Scale)(a::Sequence) = *(𝒮, a)
Base.:*(𝒮::Scale, a::Sequence) = scale(a, 𝒮.value)

function scale(a::Sequence, γ)
    𝒮 = Scale(γ)
    space_a = space(a)
    new_space = image(𝒮, space_a)
    CoefType = _coeftype(𝒮, space_a, eltype(a))
    c = Sequence(new_space, Vector{CoefType}(undef, dimension(new_space)))
    _apply!(c, 𝒮, a)
    return c
end

function scale!(c::Sequence, a::Sequence, γ)
    𝒮 = Scale(γ)
    space(c) == image(𝒮, space(a)) || return throw(DomainError)
    _apply!(c, 𝒮, a)
    return c
end

function project(𝒮::Scale, domain::VectorSpace, codomain::VectorSpace, ::Type{T}) where {T}
    _iscompatible(domain, codomain) || return throw(DimensionMismatch)
    ind_domain = _findposition_nzind_domain(𝒮, domain, codomain)
    ind_codomain = _findposition_nzind_codomain(𝒮, domain, codomain)
    C = LinearOperator(domain, codomain, sparse(ind_codomain, ind_domain, zeros(T, length(ind_domain)), dimension(codomain), dimension(domain)))
    _project!(C, 𝒮)
    return C
end

function project!(C::LinearOperator, 𝒮::Scale)
    _iscompatible(domain(C), codomain(C)) || return throw(DimensionMismatch)
    coefficients(C) .= zero(eltype(C))
    _project!(C, 𝒮)
    return C
end

_findposition_nzind_domain(𝒮::Scale, domain, codomain) =
    _findposition(_nzind_domain(𝒮, domain, codomain), domain)

_findposition_nzind_codomain(𝒮::Scale, domain, codomain) =
    _findposition(_nzind_codomain(𝒮, domain, codomain), codomain)

# Sequence spaces

image(𝒮::Scale{<:NTuple{N,Any}}, s::TensorSpace{<:NTuple{N,BaseSpace}}) where {N} =
    TensorSpace(map((γᵢ, sᵢ) -> image(Scale(γᵢ), sᵢ), 𝒮.value, spaces(s)))

_coeftype(𝒮::Scale, s::TensorSpace, ::Type{T}) where {T} =
    @inbounds promote_type(_coeftype(Scale(𝒮.value[1]), s[1], T), _coeftype(Scale(Base.tail(𝒮.value)), Base.tail(s), T))
_coeftype(𝒮::Scale, s::TensorSpace{<:Tuple{BaseSpace}}, ::Type{T}) where {T} =
    @inbounds _coeftype(Scale(𝒮.value[1]), s[1], T)

function _apply!(c::Sequence{<:TensorSpace}, 𝒮::Scale, a)
    space_a = space(a)
    A = _no_alloc_reshape(coefficients(a), dimensions(space_a))
    C = _no_alloc_reshape(coefficients(c), dimensions(space(c)))
    _apply!(C, 𝒮, space_a, A)
    return c
end

_apply!(C, 𝒮::Scale, space::TensorSpace{<:NTuple{N₁,BaseSpace}}, A::AbstractArray{T,N₂}) where {N₁,T,N₂} =
    @inbounds _apply!(C, Scale(𝒮.value[1]), space[1], Val(N₂-N₁+1), _apply!(C, Scale(Base.tail(𝒮.value)), Base.tail(space), A))

_apply!(C, 𝒮::Scale, space::TensorSpace{<:Tuple{BaseSpace}}, A::AbstractArray) =
    @inbounds _apply!(C, Scale(𝒮.value[1]), space[1], A)

for (_f, __f) ∈ ((:_nzind_domain, :__nzind_domain), (:_nzind_codomain, :__nzind_codomain))
    @eval begin
        $_f(𝒮::Scale{<:NTuple{N,Any}}, domain::TensorSpace{<:NTuple{N,BaseSpace}}, codomain::TensorSpace{<:NTuple{N,BaseSpace}}) where {N} =
            Base.Iterators.ProductIterator($__f(𝒮, domain, codomain))
        $__f(𝒮::Scale, domain::TensorSpace, codomain) =
            @inbounds ($_f(Scale(𝒮.value[1]), domain[1], codomain[1]), $__f(Scale(Base.tail(𝒮.value)), Base.tail(domain), Base.tail(codomain))...)
        $__f(𝒮::Scale, domain::TensorSpace{<:Tuple{BaseSpace}}, codomain) =
            @inbounds ($_f(Scale(𝒮.value[1]), domain[1], codomain[1]),)
    end
end

function _project!(C::LinearOperator{<:SequenceSpace,<:SequenceSpace}, 𝒮::Scale)
    domain_C = domain(C)
    codomain_C = codomain(C)
    CoefType = eltype(C)
    @inbounds for (α, β) ∈ zip(_nzind_codomain(𝒮, domain_C, codomain_C), _nzind_domain(𝒮, domain_C, codomain_C))
        C[α,β] = _nzval(𝒮, domain_C, codomain_C, CoefType, α, β)
    end
    return C
end

@generated function _nzval(𝒮::Scale{<:NTuple{N,Any}}, domain::TensorSpace{<:NTuple{N,BaseSpace}}, codomain::TensorSpace{<:NTuple{N,BaseSpace}}, ::Type{T}, α, β) where {N,T}
    p = :(_nzval(Scale(𝒮.value[1]), domain[1], codomain[1], T, α[1], β[1]))
    for i ∈ 2:N
        p = :(_nzval(Scale(𝒮.value[$i]), domain[$i], codomain[$i], T, α[$i], β[$i]) * $p)
    end
    return p
end

# Taylor

image(::Scale, s::Taylor) = s

_coeftype(::Scale{T}, ::Taylor, ::Type{S}) where {T,S} = promote_type(T, S)

function _apply!(c::Sequence{Taylor}, 𝒮::Scale, a)
    γ = 𝒮.value
    if isone(γ)
        coefficients(c) .= coefficients(a)
    else
        @inbounds c[0] = a[0]
        γⁱ = one(γ)
        @inbounds for i ∈ 1:order(a)
            γⁱ *= γ
            c[i] = a[i]*γⁱ
        end
    end
    return c
end

function _apply!(C, 𝒮::Scale, space::Taylor, ::Val{D}, A) where {D}
    γ = 𝒮.value
    if !isone(γ)
        γⁱ = one(γ)
        @inbounds for i ∈ 1:order(space)
            γⁱ *= γ
            selectdim(C, D, i+1) .*= γⁱ
        end
    end
    return C
end

function _apply!(C::AbstractArray{T,N}, 𝒮::Scale, space::Taylor, A) where {T,N}
    γ = 𝒮.value
    if isone(γ)
        C .= A
    else
        @inbounds selectdim(C, N, 1) .= selectdim(A, N, 1)
        γⁱ = one(γ)
        @inbounds for i ∈ 1:order(space)
            γⁱ *= γ
            selectdim(C, N, i+1) .= γⁱ .* selectdim(A, N, i+1)
        end
    end
    return C
end

_nzind_domain(::Scale, domain::Taylor, codomain::Taylor) = 0:min(order(domain), order(codomain))
_nzind_codomain(::Scale, domain::Taylor, codomain::Taylor) = 0:min(order(domain), order(codomain))
function _nzval(𝒮::Scale, ::Taylor, ::Taylor, ::Type{T}, i, j) where {T}
    γ = 𝒮.value
    if isone(γ)
        return one(T)
    else
        return convert(T, γ^i)
    end
end

# Fourier

image(𝒮::Scale, s::Fourier) = Fourier(order(s), frequency(s)*𝒮.value)

_coeftype(::Scale, ::Fourier, ::Type{T}) where {T} = T

function _apply!(c::Sequence{<:Fourier}, ::Scale, a)
    coefficients(c) .= coefficients(a)
    return c
end

_apply!(C, ::Scale, ::Fourier, ::Val, A) = C

function _apply!(C, ::Scale, ::Fourier, A)
    C .= A
    return C
end

function _nzind_domain(::Scale, domain::Fourier, codomain::Fourier)
    ord = min(order(domain), order(codomain))
    return -ord:ord
end
function _nzind_codomain(::Scale, domain::Fourier, codomain::Fourier)
    ord = min(order(domain), order(codomain))
    return -ord:ord
end
_nzval(::Scale, ::Fourier, ::Fourier, ::Type{T}, i, j) where {T} = one(T)

# Chebyshev

image(::Scale, s::Chebyshev) = s

_coeftype(::Scale{T}, ::Chebyshev, ::Type{S}) where {T,S} = promote_type(T, S)

function _apply!(c::Sequence{Chebyshev}, 𝒮::Scale, a)
    γ = 𝒮.value
    if isone(γ)
        coefficients(c) .= coefficients(a)
    else # TODO: lift restriction
        return throw(DomainError)
    end
    return c
end

function _apply!(C, 𝒮::Scale, space::Chebyshev, ::Val{D}, A) where {D}
    γ = 𝒮.value
    isone(γ) || return throw(DomainError) # TODO: lift restriction
    return C
end

function _apply!(C::AbstractArray{T,N}, 𝒮::Scale, space::Chebyshev, A) where {T,N}
    γ = 𝒮.value
    if isone(γ)
        C .= A
    else # TODO: lift restriction
        return throw(DomainError)
    end
    return C
end

_nzind_domain(::Scale, domain::Chebyshev, codomain::Chebyshev) = 0:min(order(domain), order(codomain))
_nzind_codomain(::Scale, domain::Chebyshev, codomain::Chebyshev) = 0:min(order(domain), order(codomain))
function _nzval(𝒮::Scale, ::Chebyshev, ::Chebyshev, ::Type{T}, i, j) where {T}
    γ = 𝒮.value
    if isone(γ)
        return one(T)
    else # TODO: lift restriction
        return throw(DomainError)
    end
end

# Cartesian spaces

image(𝒮::Scale, s::CartesianPower) =
    CartesianPower(image(𝒮, space(s)), nb_cartesian_product(s))

image(𝒮::Scale, s::CartesianProduct) =
    CartesianProduct(map(sᵢ -> image(𝒮, sᵢ), spaces(s)))

_coeftype(𝒮::Scale, s::CartesianPower, ::Type{T}) where {T} =
    _coeftype(𝒮, space(s), T)

_coeftype(𝒮::Scale, s::CartesianProduct, ::Type{T}) where {T} =
    @inbounds promote_type(_coeftype(𝒮, s[1], T), _coeftype(𝒮, Base.tail(s), T))
_coeftype(𝒮::Scale, s::CartesianProduct{<:Tuple{VectorSpace}}, ::Type{T}) where {T} =
    @inbounds _coeftype(𝒮, s[1], T)

function _apply!(c::Sequence{<:CartesianPower}, 𝒮::Scale, a)
    @inbounds for i ∈ 1:nb_cartesian_product(space(c))
        _apply!(component(c, i), 𝒮, component(a, i))
    end
    return c
end
function _apply!(c::Sequence{CartesianProduct{T}}, 𝒮::Scale, a) where {N,T<:NTuple{N,VectorSpace}}
    @inbounds _apply!(component(c, 1), 𝒮, component(a, 1))
    @inbounds _apply!(component(c, 2:N), 𝒮, component(a, 2:N))
    return c
end
function _apply!(c::Sequence{CartesianProduct{T}}, 𝒮::Scale, a) where {T<:Tuple{VectorSpace}}
    @inbounds _apply!(component(c, 1), 𝒮, component(a, 1))
    return c
end

function _findposition_nzind_domain(𝒮::Scale, domain::CartesianSpace, codomain::CartesianSpace)
    u = map((dom, codom) -> _findposition_nzind_domain(𝒮, dom, codom), spaces(domain), spaces(codomain))
    len = sum(length, u)
    v = Vector{Int}(undef, len)
    δ = δδ = 0
    @inbounds for (i, uᵢ) in enumerate(u)
        δ_ = δ
        δ += length(uᵢ)
        view(v, 1+δ_:δ) .= δδ .+ uᵢ
        δδ += dimension(domain[i])
    end
    return v
end

function _findposition_nzind_codomain(𝒮::Scale, domain::CartesianSpace, codomain::CartesianSpace)
    u = map((dom, codom) -> _findposition_nzind_codomain(𝒮, dom, codom), spaces(domain), spaces(codomain))
    len = sum(length, u)
    v = Vector{Int}(undef, len)
    δ = δδ = 0
    @inbounds for (i, uᵢ) in enumerate(u)
        δ_ = δ
        δ += length(uᵢ)
        view(v, 1+δ_:δ) .= δδ .+ uᵢ
        δδ += dimension(codomain[i])
    end
    return v
end

function _project!(C::LinearOperator{<:CartesianSpace,<:CartesianSpace}, 𝒮::Scale)
    @inbounds for i ∈ 1:nb_cartesian_product(domain(C))
        _project!(component(C, i, i), 𝒮)
    end
    return C
end


text:
export halfpi, hrmint, hx2dp

function _halfpi()
    ccall((:halfpi_c, libcspice), SpiceDouble, ())
end

"""
    halfpi()

!!! warning "Deprecated"
    Use `π/2` instead.
"""
halfpi

@deprecate halfpi() π/2

"""
    hrmint(xvals, yvals, x)

Evaluate a Hermite interpolating polynomial at a specified abscissa value.

### Arguments ###

- `xvals`: Abscissa values
- `yvals`: Ordinate and derivative values
- `x`: Point at which to interpolate the polynomial

### Output ###

- `f`: Interpolated function value at `x`
- `df`: Interpolated function's derivative at `x`

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/hrmint_c.html)
"""
function hrmint(xvals, yvals, x)
    n = length(xvals)
    @checkdims 2n yvals
    work = Array{SpiceDouble}(undef, 4n)
    f = Ref{SpiceDouble}()
    df = Ref{SpiceDouble}()
    ccall((:hrmint_c, libcspice), Cvoid,
          (SpiceInt, Ref{SpiceDouble}, Ref{SpiceDouble}, SpiceDouble, Ref{SpiceDouble},
           Ref{SpiceDouble}, Ref{SpiceDouble}),
          n, xvals, yvals, x, work, f, df)
    handleerror()
    f[], df[]
end

"""
    hx2dp(str)

Convert a string representing a double precision number in a base 16 "scientific notation"
into its equivalent double precision number.

### Arguments ###

- `str`: Hex form string to convert to double precision

### Output ###

- `dp`: Double precision value to be returned

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/hx2dp_c.html)
"""
function hx2dp(str)
    dp = Ref{SpiceDouble}()
    error = Ref{SpiceBoolean}()
    errmsg = Array{SpiceChar}(undef, 46)
    ccall((:hx2dp_c, libcspice), Cvoid,
        (Cstring, SpiceInt, Ref{SpiceDouble}, Ref{SpiceBoolean}, Ref{SpiceChar}),
        str, 46, dp, error, errmsg)
    if Bool(error[])
        throw(SpiceError(chararray_to_string(errmsg)))
    end
    dp[]
end


text:
using DemeNet: DemeSpec, Deme, Signer, Certificate, save, DemeID, ID
using PeaceVote: KeyChain, members, proposals, attest, voters, BraiderConfig, RecorderConfig, BraidChainConfig, BraidChainServer, Vote, Proposal, BraidChain, load, record, braid!, sync!
using PeaceCypher

### Some clenup first ###
for dir in [homedir() * "/.demenet/"]
    isdir(dir) && rm(dir,recursive=true)
end

### Setting up Deme with PeaceVote ###

demespec = DemeSpec("PeaceDeme",:default,:PeaceCypher,:default,:PeaceCypher,:PeaceVote)
save(demespec) ### Necessary to connect with Mixer
uuid = demespec.uuid
deme = Deme(demespec)

maintainer = Signer(deme,"maintainer")
server = Signer(deme,"server")

MIXER_ID = server.id
SERVER_ID = server.id
MAINTAINER_ID = maintainer.id

MIXER_PORT = 3001 # Self mixing
BRAIDER_PORT = 3000
REGISTRATOR_PORT = 3002
VOTING_PORT = 3003
PROPOSAL_PORT = 3004
SYNC_PORT = 3005

braiderconfig = BraiderConfig(BRAIDER_PORT,MIXER_PORT,UInt8(3),UInt8(64),SERVER_ID,DemeID(uuid,MIXER_ID))
recorderconfig = RecorderConfig([MAINTAINER_ID,SERVER_ID],server.id,REGISTRATOR_PORT,VOTING_PORT,PROPOSAL_PORT)
braidchainconfig = BraidChainConfig(SERVER_ID,MIXER_PORT,SYNC_PORT,braiderconfig,recorderconfig)

braidchain = BraidChain(braidchainconfig,deme)
system = BraidChainServer(braidchain,server)

### Testing the server ###

maintainer = Signer(deme,"maintainer")

for i in 1:3
    account = "account$i"
    keychain = KeyChain(deme,account)
    cert = Certificate(keychain.member.id,maintainer)
    @show record(braidchain,cert)
end

# Now let's test braiding

@sync for i in 1:3
    @async begin
        account = "account$i"
        keychain = KeyChain(deme,account)
        braid!(braidchain,keychain)
    end
end

# Proposing

keychain = KeyChain(deme,"account2")
proposal = Proposal("Found peace for a change?",["yes","no","maybe"])
cert = Certificate(proposal,keychain.member)
record(braidchain,cert)

sleep(1)

loadedledger = load(braidchain)
messages = attest(loadedledger,braidchain.deme.notary)
index = proposals(messages)[1]
proposal = messages[index]

# Now we can vote

for i in 1:3
    account = "account$i"
    keychain = KeyChain(deme,account)

    option = Vote(index,rand(1:length(proposal.document.options)))
    cert = Certificate(option,braidchain,keychain)
    record(braidchain,cert)
end

# Now let's count
sleep(1)

@show tally = count(index,braidchain)

# Let's test synchronization

demesync = Deme(demespec)
bcsync = BraidChain(braidchainconfig,demesync)
sync!(bcsync)
@show tally = count(index,bcsync)


text:
module TestPackage

using SymbolicUtils
include("extra_file.jl")
export(test_f)
export(test_f_2)

end


text:
# Diagonal Approximation for trace of inverse of matrix
# Values from Extrapolation Methods for Estimating Trace of the Matrix Inverse by P. Fika

export diagapp

using LinearAlgebra

struct Diagappspace{T, TM <: AbstractMatrix{T}, TV <: AbstractVector{T}, F<:Function}
    A::TM
    x::TV
    y::TV
    randfunc::F
end

function Diagappspace(A::AbstractMatrix, randfunc::Function)
    x = randfunc(eltype(A) <: Integer ? Float64 : eltype(A), size(A, 1))
    y = similar(x)
    return Diagappspace(A, x, y, randfunc)
end

# Extrapolation for c-1 and calculating value of v0
function v0(w::Diagappspace)

    copyto!(w.x, w.randfunc(-1:2:1, size(w.A, 1)))

    c0 = dot(w.x, w.x)
    mul!(w.y, w.A, w.x)
    c1 = dot(w.y, w.x)
    c2 = dot(w.y, w.y)
    mul!(w.x, w.A, w.y)
    c3 = dot(w.y, w.x)

    p = (c0 * c2)/(c1^2)

    v0 = log10((c1^2)/(c0 * c2)) / log10((c1 * c3)/(c2^2))
end

# Calculating the Approximation for ith value of inverse diagonal
function dfun(A, i, v)
    isum = zero(eltype(A))
    s = size(A, 1)
    for k in 1:s
        isum = isum + A[k, i]^2
    end
    #v = v0(A)
    p = isum / (A[i, i]^2)
    d = 1 / (p^v * A[i, i])
end

# Calculating Diagonal Approximation for Matrix Inverse
# This works a SPD Matrix with low condition number
"""
    diagapp(A::AbstractMatrix; randfunc=Base.rand)

Diagonal Approximation algorithm for inverse of matrix.
# Arguments
 - `A` : Symmetric Positive Definite Matrix with Low Condtion Number (k < 500)
 - `randfunc` : Random function for the particular type of matrix A
                (An example can be found in test/diagapp.jl)
"""
function diagapp(A; randfunc::Function=Base.rand)
    tr = zero(eltype(A))
    s = size(A, 1)
    w = Diagappspace(A, randfunc)
    v = v0(w)
    for i in 1:s
        tr = tr + dfun(w.A, i, v)
    end
    tr
end


text:
using DFTK
using Test

include("testcases.jl")

function test_kernel_unpolarized(termtype)
    Ecut=2
    kgrid = [2, 2, 2]
    testcase = silicon
    ε = 1e-8

    @testset "Kernel $(typeof(termtype))" begin
        spec = ElementPsp(testcase.atnum, psp=load_psp(testcase.psp))
        model = Model(testcase.lattice; atoms=[spec => testcase.positions],
                      terms=[termtype])
        basis = PlaneWaveBasis(model, Ecut; kgrid=kgrid)
        term  = only(basis.terms)

        ρ0 = guess_density(basis)
        dρ = from_real(basis, randn(size(ρ0)))

        ρ_minus = ρ0 - ε * dρ
        ρ_plus  = ρ0 + ε * dρ
        pot_minus = DFTK.ene_ops(term, nothing, nothing; ρ=ρ_minus).ops[1].potential
        pot_plus  = DFTK.ene_ops(term, nothing, nothing; ρ=ρ_plus).ops[1].potential
        dV = (pot_plus - pot_minus) / (2ε)

        dV_apply = DFTK.apply_kernel(term, dρ; ρ=ρ0)[1]
        kernel = DFTK.compute_kernel(term; ρ=ρ0)
        dV_compute = reshape(kernel * vec(dρ.real), size(dρ))
        @test norm(dV - dV_apply.real) < 100ε
        @test norm(dV - dV_compute)    < 100ε
    end
end

function test_kernel_collinear(termtype)
    Ecut=2
    kgrid = [2, 2, 2]
    testcase = silicon
    ε = 1e-8

    @testset "Kernel $(typeof(termtype)) (collinear)" begin
        spec = ElementPsp(testcase.atnum, psp=load_psp(testcase.psp))
        magnetic_moments = [spec => 2rand(2)]
        model = Model(testcase.lattice; atoms=[spec => testcase.positions],
                      terms=[termtype], magnetic_moments=magnetic_moments,
                      spin_polarization=:collinear)
        basis = PlaneWaveBasis(model, Ecut; kgrid=kgrid)
        iup   = DFTK.krange_spin(basis, 1)[1]  # First spin-up k-point
        idown = DFTK.krange_spin(basis, 2)[1]  # First spin-down k-point
        term  = only(basis.terms)

        ρ0      = guess_density(basis)
        dρ      = from_real(basis, randn(size(ρ0)))
        ρspin0  = guess_spin_density(basis, magnetic_moments)
        dρspin  = from_real(basis, randn(size(ρspin0)))

        ρ_minus     = ρ0     - ε * dρ
        ρspin_minus = ρspin0 - ε * dρspin
        ρ_plus      = ρ0     + ε * dρ
        ρspin_plus  = ρspin0 + ε * dρspin

        ops_minus = DFTK.ene_ops(term, nothing, nothing; ρ=ρ_minus, ρspin=ρspin_minus).ops
        ops_plus  = DFTK.ene_ops(term, nothing, nothing; ρ=ρ_plus,  ρspin=ρspin_plus).ops
        dVα = (ops_plus[  iup].potential - ops_minus[  iup].potential) / (2ε)
        dVβ = (ops_plus[idown].potential - ops_minus[idown].potential) / (2ε)

        dVα_apply, dVβ_apply = DFTK.apply_kernel(term, dρ, dρspin; ρ=ρ0, ρspin=ρspin0)

        kernel = DFTK.compute_kernel(term; ρ=ρ0, ρspin=ρspin0)
        dρcat  = vcat(vec(dρ.real), vec(dρspin.real))
        dV_matrix = reshape(kernel * dρcat, size(dρ)..., 2)

        @test norm(dVα - dVα_apply.real) < 100ε
        @test norm(dVβ - dVβ_apply.real) < 100ε
        @test norm(dVα - dV_matrix[:, :, :, 1]) < 100ε
        @test norm(dVβ - dV_matrix[:, :, :, 2]) < 100ε
    end
end

test_kernel_unpolarized(PowerNonlinearity(1.2, 2.0))
test_kernel_unpolarized(Hartree())
test_kernel_unpolarized(Xc(:lda_xc_teter93))

test_kernel_collinear(Hartree())
test_kernel_collinear(PowerNonlinearity(1.2, 2.5))
test_kernel_collinear(Xc(:lda_xc_teter93))


text:
module ReferenceFrameRotations

import Base: +, -, *, /, \, ≈, ==
import Base: conj, convert, copy, display, eltype, firstindex, getindex, imag
import Base: inv, iterate, lastindex, ndims, one, real, setindex!, show
import Base: zero, zeros
import Base: Broadcast.broadcastable
import LinearAlgebra: norm

using Crayons
using LinearAlgebra
using Printf
using StaticArrays

# Re-export `I` from LinearAlgebra.
export I

################################################################################
#                                    Types
################################################################################

include("types.jl")

################################################################################
#                                  Constants
################################################################################

# Pre-defined crayons.
const _reset_crayon = Crayon(reset = true)
const _crayon_bold  = crayon"bold"
const _crayon_g     = crayon"bold green"
const _crayon_u     = crayon"bold blue"
const _crayon_y     = crayon"bold yellow"

# Escape sequences related to the crayons.
const _b = _crayon_bold
const _d = _reset_crayon
const _g = _crayon_g
const _y = _crayon_y
const _u = _crayon_u

################################################################################
#                                   Includes
################################################################################

include("angle.jl")
include("angleaxis.jl")
include("compose_rotations.jl")
include("dcm.jl")
include("inv_rotations.jl")
include("quaternion.jl")

include("deprecations.jl")

const _CONVERSIONS_DIR = joinpath(@__DIR__, "conversions")

include(joinpath(_CONVERSIONS_DIR, "angle_to_angle.jl"))
include(joinpath(_CONVERSIONS_DIR, "angle_to_angleaxis.jl"))
include(joinpath(_CONVERSIONS_DIR, "angle_to_dcm.jl"))
include(joinpath(_CONVERSIONS_DIR, "angle_to_quat.jl"))
include(joinpath(_CONVERSIONS_DIR, "angle_to_rot.jl"))
include(joinpath(_CONVERSIONS_DIR, "angleaxis_to_angle.jl"))
include(joinpath(_CONVERSIONS_DIR, "angleaxis_to_dcm.jl"))
include(joinpath(_CONVERSIONS_DIR, "angleaxis_to_quat.jl"))
include(joinpath(_CONVERSIONS_DIR, "dcm_to_angle.jl"))
include(joinpath(_CONVERSIONS_DIR, "dcm_to_angleaxis.jl"))
include(joinpath(_CONVERSIONS_DIR, "dcm_to_quat.jl"))
include(joinpath(_CONVERSIONS_DIR, "quat_to_angle.jl"))
include(joinpath(_CONVERSIONS_DIR, "quat_to_angleaxis.jl"))
include(joinpath(_CONVERSIONS_DIR, "quat_to_dcm.jl"))
include(joinpath(_CONVERSIONS_DIR, "smallangle_to_dcm.jl"))
include(joinpath(_CONVERSIONS_DIR, "smallangle_to_quat.jl"))
include(joinpath(_CONVERSIONS_DIR, "smallangle_to_rot.jl"))

end # module


text:
using Resample
using Documenter
using PlutoStaticHTML

"""
Run all Pluto notebooks (".jl" files) in `tutorials_dir` and write output to Markdown files.
"""
function build_tutorials()
    println("Building tutorials")
    dir = joinpath(pkgdir(Resample), "docs", "src", "notebooks")
    # Evaluate notebooks in the same process to avoid having to recompile from scratch each time.
    # This is similar to how Documenter and Franklin evaluate code.
    # Note that things like method overrides may leak between notebooks!
    use_distributed = false
    output_format = documenter_output
    bopts = BuildOptions(dir; use_distributed, output_format)
    build_notebooks(bopts)
    return nothing
end

# Build the notebooks; defaults to "true".
if get(ENV, "BUILD_DOCS_NOTEBOOKS", "true") == "true"
    build_tutorials()
end

pages = [
    "Home" => "index.md",
    "Tutorials" => [
        "SMOTE" => "notebooks/smote.md"
    ]
]

makedocs(;
    modules=[Resample],
    sitename="Resample.jl",
    authors = "Rik Huijzer and contributors",
    format=Documenter.HTML(;
        canonical="https://rikhuijzer.github.io/Resample.jl",
        # Using MathJax3 since Pluto uses that engine too.
        mathengine=Documenter.MathJax3(),
        prettyurls=get(ENV, "CI", "false") == "true",
    ),
    pages
)

deploydocs(; repo="github.com/rikhuijzer/Resample.jl", devbranch="main")


text:
# Test of Bellman Ford algorithm for shortest paths

using Graphs
using Base.Test

# g1: the example in CLRS (2nd Ed.)
g1 = simple_inclist(5)

g1_wedges = [
    (1, 2, 10.),
    (1, 3, 5.),
    (2, 3, 2.),
    (3, 2, 3.),
    (2, 4, 1.),
    (3, 5, 2.),
    (4, 5, 4.),
    (5, 4, 6.),
    (5, 1, 7.),
    (3, 4, 9.) ]

ne = length(g1_wedges)
eweights1 = zeros(ne)
for i = 1 : ne
    we = g1_wedges[i]
    add_edge!(g1, we[1], we[2])
    eweights1[i] = we[3]
end

@assert num_vertices(g1) == 5
@assert num_edges(g1) == 10

s1 = bellman_ford_shortest_paths(g1, eweights1, [1])

@test s1.parents == [1, 3, 1, 2, 3]
@test s1.dists == [0., 8., 5., 9., 7.]
@test !has_negative_edge_cycle(g1, eweights1)

immutable MyEdge{V}
    index::Int
    source::V
    target::V
    dist::Float64
end
Graphs.target{V}(e::MyEdge{V}, g::AbstractGraph{V}) = e.target
Graphs.source{V}(e::MyEdge{V}, g::AbstractGraph{V}) = e.source
Graphs.edge_index(e::MyEdge) = e.index
g2 = inclist([i for i=1:10], MyEdge{Int})

for i = 2:10
    for j = 1:(i - 1)
        index = div((i-1)*(i-2), 2) + j
        weight = (i == j + 1) ? 1.0 : (i - j + 1.0)
        edge = MyEdge(index, j, i, weight)
        add_edge!(g2, edge)
    end
end

type MyEdgePropertyInspector{T} <: AbstractEdgePropertyInspector{T} end

Graphs.edge_property{T,V}(inspector::MyEdgePropertyInspector{T}, e::MyEdge, g::AbstractGraph{V}) = e.dist
insp =  MyEdgePropertyInspector{Float64}()
s2 = bellman_ford_shortest_paths(g2, insp, [1])
@test s2.dists == [i * 1.0 for i=0:9]
@test s2.parents == [i ==0 ? 1 : i for i = 0:9]
@test !has_negative_edge_cycle(g2, insp)
add_edge!(g2, MyEdge{Int}(46, 10, 1, -10.0))
@test has_negative_edge_cycle(g2, insp)
@test_throws NegativeCycleError bellman_ford_shortest_paths(g2, insp, [1])

# g1: the example in CLR[~S] (1st Ed. Figure 25.7)
g3 = simple_inclist(5)

g3_wedges = [
    (1, 2,  5.),
    (1, 3, -4.),
    (1, 4,  8.),
    (2, 1, -2.),
    (3, 2,  7.),
    (3, 5,  2.),
    (4, 2, -3.),
    (4, 3,  9.),
    (5, 4,  7.),
    (5, 1,  6.)]

ne3 = length(g3_wedges)
eweights3 = zeros(ne3)
for i = 1 : ne3
    we = g3_wedges[i]
    add_edge!(g3, we[1], we[2])
    eweights3[i] = we[3]
end

@assert num_vertices(g3) == 5
@assert num_edges(g3) == 10

s3 = bellman_ford_shortest_paths(g3, eweights3, [5])

@test s3.parents == [2, 4, 1, 5, 5]
@test s3.dists == [2.0, 4.0, -2.0, 7.0, 0.0]
@test !has_negative_edge_cycle(g3, eweights1)


text:
__precompile__()

module InterpHalo

using JLD, Interpolations

export intH
path = @__DIR__
H1 = load(path*"/H1Data.jld")
H2 = load(path*"/H2Data.jld")

const ORB1 = H1["ORB"]

const ORB2 = H2["ORB"]

intORB1 = interpolate(ORB1, BSpline(Quadratic(Reflect(OnCell()))))
intORB2 = interpolate(ORB2, BSpline(Quadratic(Reflect(OnCell()))))

interpH1(x,y,z) = intORB1(x,y,z)
interpH2(x,y,z) = intORB2(x,y,z)

include("functions_interp.jl")

end # module


text:
module Storage

#
#  The object managing output storage for the system.
#
#  The object has two purposes:
#
#    * gather data that arrives per-frame and package it into
#      readable output files, one per frame.
#
#    * gather data that will be stored in a single file for all
#      time points at the end of the simulation
#
#   Note: the data may overlap between storage types.
#
#

using Calendar
import Calendar.CalendarTime

# the global storage manager object
global sm = nothing

type StorageManager
    output_dir :: String
    log_name :: String
    frame_prefix :: String

    tag_cfg :: Dict
    current_frame :: Dict
    frame_ndx :: Int
    flushed :: Bool

    ts_store :: Dict

    log_io

    StorageManager(out_dir, log_name, frame_prefix) = new(out_dir, log_name, frame_prefix, Dict{String,Any}(), Dict{String,Any}(), 1, false, Dict{String,Any}(), nothing)
end

function sopen(out_dir, log_name, frame_prefix)
    global sm
    sm = StorageManager(out_dir, log_name, frame_prefix)
    log_path = join([out_dir, log_name], "/")
    sm.log_io = open(log_path, "w")
    println(sm.log_io, "Initialized log at [$log_path]")
end

function sclose()
    global sm
    close(sm.log_io)
    sm = nothing
end

function setup_tag(tag :: String, stdout :: Bool, logout :: Bool)
    sm.tag_cfg[tag] = (stdout, logout)
end

function next_frame()
    flush_frame()
    sm.current_frame = Dict{String,Any}()
    sm.flushed = false
end

function flush_frame()

    # if the frame is already flushed, don't call again
    if sm.flushed == false

        f = sm.current_frame
        kk = collect(keys(f))

        # another construction to avoid "," after last item
        if length(kk) > 0

            # render the data in python mode
            ff_name = join([sm.output_dir, string(sm.frame_prefix, string(sm.frame_ndx))], "/")
            io = open(ff_name, "w")
            println(io, "{")

            python_render_item(io, kk[1], f[kk[1]])
            for k in kk[2:]
                println(io, ",")
                python_render_item(io, k, f[k])
            end
            println(io)
            println(io, "}")
            close(io)
        end

        sm.frame_ndx += 1

        # log end of frame to io
        println(sm.log_io, "end-of-frame [$ff_name]")

        # indicate this frame has been flushed
        sm.flushed = true
    end
end

function spush(tag :: String, data)
    c = sm.tag_cfg[tag]
    sm.current_frame[tag] = data

    # if stdout requested print it
    if c[1] println(string(tag, " : ", data)) end

    # if logout is requested
    if c[2] println(sm.log_io, string(tag, " : ", data)) end
end

function python_render_item(io::IO, k::String, v::CalendarTime)
    show(io,k)
    print(io, " : ")
    print(io, "datetime($(year(v)), $(month(v)), $(day(v)), $(hour(v)), $(minute(v)), $(second(v)))")
end

function python_render_item(io::IO, k::String, v::Array)
    show(io, k)
    print(io, " : ")

    # construct to avoid "," after last item in list
    s = size(v)
    print(io, length(s) > 1 ? "np.reshape([" : "[")
    if length(v) > 0
        python_render_item(io, v[1])
        for i in v[2:]
            print(io, ", ")
            python_render_item(io, i)
        end
    end
    println(io, length(s) > 1 ? "], $s, order = 'F')" : "]")
end

function python_render_item(io::IO, k::String, v::Any)
    show(io, k)
    print(io, " : ")
    python_render_item(io, v)
end

# specialize for floats which are otherwise printed with too many digits
python_render_item(io::IO, f::Float64) = isnan(f) ? @printf(io, "np.nan") : @printf(io, "%f", f)
python_render_item(io::IO, f) = show(io, f)

end


text:
#=
build:
- Julia version:
- Author: Dan
- Date: 2021-08-20
=#


text:
mutable struct minHeap
    valArr::Array{Float64,1};
	indArr::Array{Int64,1};
    size::Int64;
    maxSize::Int64;
end

function initHeap(n::Int64)
valArr = zeros(n);
indArr = zeros(Int64,n);
heap = minHeap(valArr,indArr,0,n);
return heap;
end

function initHeap(valArr::Array{Float64,1},indArr::Array{Int64,1})
if length(valArr)!=length(indArr)
	error("initHeap::unevenly lengthed arrays");
end
valArr[:] .= 0.0;
indArr[:] .= 0;
heap = minHeap(valArr,indArr,0,length(valArr));
return heap;
end

function insertToHeap(heap::minHeap,val::Float64,ind::Int64)
if heap.size == heap.maxSize
	println("Allocating more memory in heap.")
	tmpArr = heap.valArr;
	heap.valArr = zeros(2*heap.maxSize);
	heap.valArr[1:heap.maxSize] = tmpArr;
	tmpArr = heap.indArr;
	heap.indArr = zeros(Int64,2*heap.maxSize);
	heap.indArr[1:heap.maxSize] = tmpArr;
	tmpArr = 0; # releasing memory???
	heap.maxSize = 2*heap.maxSize;
end
heap.size += 1;
now = heap.size;
heap.valArr[now] = val;
heap.indArr[now] = ind;

# par = div(now,2);
# while now > 1 && val < heap.valArr[par]
    # # @inbounds if ()
    # @inbounds heap.valArr[now] = heap.valArr[par];
	# @inbounds heap.indArr[now] = heap.indArr[par];
    # @inbounds heap.valArr[par] = val;
	# @inbounds heap.indArr[par] = ind;
    # now = par;
	# par = div(now,2);
    # # else
        # # break;
    # # end
# end

while now > 1
    par = div(now,2);
    @inbounds if (val < heap.valArr[par])
        @inbounds heap.valArr[now] = heap.valArr[par];
		@inbounds heap.indArr[now] = heap.indArr[par];
        @inbounds heap.valArr[par] = val;
		@inbounds heap.indArr[par] = ind;
        now = par;
    else
        break;
    end
end
return;
end

function getMin(heap::minHeap)
if heap.size >= 1
	if (heap.size==0)
		least = NaN;
		ind = 0;
		return least,ind;
	elseif heap.size==1
		least = heap.valArr[1];
		ind = heap.indArr[1];
		heap.valArr[1] = 0.0;
		heap.indArr[1] = 0;
		heap.size = 0;
		return least,ind;
	end
end

# taking the root node as the minimum
@inbounds least = heap.valArr[1];
@inbounds least_ind = heap.indArr[1];
# setting the last element as x and reducing the size of the heap.
@inbounds x = heap.valArr[heap.size];
@inbounds x_ind = heap.indArr[heap.size];
@inbounds heap.valArr[heap.size] = 0.0;
@inbounds heap.indArr[heap.size] = 0;
heap.size = heap.size - 1;

curr = 1;

while true
    leftChild = 2*curr; # 2k is the left child of k
    if leftChild > heap.size # no left child.
		break;
	else
		rightChild = curr*2 + 1;
        @inbounds if rightChild > heap.size || heap.valArr[leftChild] < heap.valArr[rightChild]
            child = leftChild;
        else
            child = rightChild;
        end
        if x > heap.valArr[child]
            @inbounds heap.valArr[curr] = heap.valArr[child];
		    @inbounds heap.indArr[curr] = heap.indArr[child];
            curr = child;
        else
            break;
        end
    end
end

@inbounds heap.valArr[curr] = x;
@inbounds heap.indArr[curr] = x_ind;
return least,least_ind;
end

function getMin(heap::minHeap,Ans::Array{Float64,1})
if heap.size >= 1
	if (heap.size==0)
		least = NaN;
		@inbounds Ans[1] = least;
		return 0;
	elseif heap.size==1
		Ans[1] = heap.valArr[1];
		ind = heap.indArr[1];
		heap.valArr[1] = 0.0;
		heap.indArr[1] = 0;
		heap.size = 0;
		return ind;
	end
end

# taking the root node as the minimum
@inbounds least = heap.valArr[1];
@inbounds least_ind = heap.indArr[1];
# setting the last element as x and reducing the size of the heap.
@inbounds x = heap.valArr[heap.size];
@inbounds x_ind = heap.indArr[heap.size];
@inbounds heap.valArr[heap.size] = 0.0;
@inbounds heap.indArr[heap.size] = 0;
heap.size = heap.size - 1;

curr = 1;
# leftChild = 2;
# while leftChild <= heap.size # no left child.
   # rightChild = curr*2 + 1;
   # child = rightChild + (heap.valArr[leftChild] < heap.valArr[rightChild] || rightChild > heap.size)*(leftChild - rightChild)

   # # @inbounds if heap.valArr[leftChild] < heap.valArr[rightChild] || rightChild > heap.size
      # # child = leftChild;
   # # else
      # # child = rightChild;
   # # end
   # if x > heap.valArr[child]
       # @inbounds heap.valArr[curr] = heap.valArr[child];
	   # @inbounds heap.indArr[curr] = heap.indArr[child];
       # curr = child;
   # else
       # break;
   # end
   # leftChild = 2*curr; # 2k is the left child of k
# end

while true
    leftChild = 2*curr; # 2k is the left child of k
    if leftChild > heap.size # no left child.
		break;
	else
		rightChild = curr*2 + 1;
        @inbounds if rightChild > heap.size || heap.valArr[leftChild] < heap.valArr[rightChild]
            child = leftChild;
        else
            child = rightChild;
        end
        if x > heap.valArr[child]
            @inbounds heap.valArr[curr] = heap.valArr[child];
		    @inbounds heap.indArr[curr] = heap.indArr[child];
            curr = child;
        else
            break;
        end
    end
end

@inbounds heap.valArr[curr] = x;
@inbounds heap.indArr[curr] = x_ind;
@inbounds Ans[1] = least;
return least_ind;
end


text:
"fake registry"
module Registry

using Test
using Comonicon

@cast add(path) = @test path === "abc"
@cast rm(path) = @test path === "abc"

end

@cast Registry


text:
"""
	testhist()
Testing distributions using χ^2
Null hypothesis: no significant difference between sample and parent populaiton
Test is based on a comparison of expected and observed occurrances (histograms)
Following steps are carried out (see Measurements and their Uncertainty by I.G. Hughes & T.P.A. Hase):
	1. create histogram using observations (`x`): O
	2. create histogram using given distributions (`d`): E
	 	* use the same bins as in 1 and distribuiton CDF * number of observations
	3. combine bins with occurance E < 5 (to avoid chi^2 to be skewed by the asymetry of Poisson)
	4. compute χ^2 = sum((O-E)^2/E)
	5. return statistics
		* good fit: χ^2 ≈ 1 and P(χ^2;v) ≈ 0.5
		* poor fit χ^2 << 1 || χ^2 > 3 and P(χ^2;v) → 1 || P(χ^2;v) < 1e-4

**Input**
* `x`: observations
* `distin`: input distribution to be tested agains `x`
* `unkpar`: number of unknown parameters used to define `dising` (e.g. =2 if only estimated mean and deviation define the normal distribution). Will be used to compute degree of freedom
* `bin_method`: is either bin method (e.g. "Scott", see `histbins`), number of bins (e.g. 10), or vector with bin start stop range ([0.1,0.2,0.3,0.4])
* `alpha`: significance level, a probability threshold below which the null hypothesis will be rejected

**Output**
* test summary as TestResult

**Example**
```
# generate sample
d = Distributions.Binomial(100,0.1)
x = rand(d,100) .+ 0.0; # must be float
# compute test
out = testhist(x,d,0,alpha=0.001,bin_method="Scott")
out.estim < out.critical

```
"""
function testhist(x::Vector{Float64},distin,unkpar::Int;
				bin_method="Scott",alpha::Float64=0.05)::TestResult
	xuse = HypoTest.prepdata(x);
	n = length(xuse);
	b = typeof(bin_method)==String || typeof(bin_method)==Int ?
			HypoTest.histbins(xuse,method=bin_method) :
			bin_method;
	# initial E (will be modified to E>5)
	E = HypoTest.theobins(distin,b,n);
	r = findall(x->x.<5.,E);
	i = 1; # count iteration to avoid infinite loop
	while !isempty(r)
		r[1]+1 < length(b) ? deleteat!(b,r[1]+1) : nothing
		r[1]+1 == length(b) ? deleteat!(b,length(b)-1) : nothing
		E = HypoTest.theobins(distin,b,n);
		r = findall(x->x.<5.,E)
		i += 1;
		i > n ? break : nothing
	end
	# Compute observed occurrances using updated bins
	O = HypoTest.countbins(xuse,b);
	# compute χ^2
	chi2 = HypoTest.pearsonschi2(E,O);
	# compare with tabled critical value
	v = (length(b)-1) - unkpar - 1; # degree of  freedom (length(b)-1 == number of bins)
	chi2crit = HypoTest.chi2critical(v,alpha);
	return TestResult("data sampled from given distribution",
					 "",chi2,chi2crit,alpha);
end

"""
	pearsonchi2(E,O)
Calculating the test-statistic using `E` expected and `O` observed values
"""
pearsonschi2(E::Vector{<:Real},O::Vector{<:Real})::Float64 = sum(((O.-E).^2)./E)

"""
	chi2critical(v,alpha)
Compute chritical chi^2 value using `v` degree of freedom and `alpha` significace level
"""
function chi2critical(v::Int,alpha::Float64)::Float64
	c = Distributions.Chisq(v);
	return Distributions.quantile(c,1-alpha);
end

max_stars_repo_name: emenems/HypoTest

max_issues_repo_name: emenems/HypoTest

max_forks_repo_name: emenems/HypoTest


text:
"""Functions for constructing parameter profiles"""
constantParameter(offset::Real=0.0) = x -> offset
export constantParameter
constant = constantParameter
export constant

function heaviside(x::Real, stepOpt::Real=1.0)
    if x < 0
        y = 0
    elseif x > 0
        y = 1
    elseif x == 0
        y = stepOpt
    else
        return nothing
    end
    y = convert(typeof(x), y)
end
heaviside(x::AbstractArray, stepOpt::Real=1.0) = heaviside.(x, stepOpt)

sigmoid(x::Real, width=1) = 1/(1+exp(-x/width))

"""
    unitStep(threshold::Real=0.0, baseline::Real=0.0, stepHeight::Real=1.0, stepOpt::Real=1.0)

Construct a constant function of 'x' that undergoes a step of size 'stepHeight' from a 'baseline' at 'threshold'. 'stepOpt' sets the value at 0. A Discontinuous type is returned, which can be called like a normal function but also contains the x-value for which there is a discontinuity.

# Examples
```julia-repl
julia> D = unitStep();
D([-1, 0, 1])
```
"""
function unitStep(threshold::Real=0.0, baseline::Real=0.0, stepHeight::Real=1.0, stepOpt::Real=1.0)
    Discontinuous(x -> (heaviside(x-threshold, stepOpt).*stepHeight)+baseline, Set([threshold]))
end
unitStep(T::Tuple, args...) = unitStep(T[1], args...) # For compatibility with other parameter profile syntax. If T is tuple, use T[1] as the threshold
export unitStep

# * A unit step that is smooth as butter
function unitsigmoid(threshold::Real=0.0, baseline::Real=0.0, stepHeight::Real=1.0, width::Real=0.1)
    Discontinuous(x -> (sigmoid(x-threshold, width).*stepHeight)+baseline, Set([threshold]))
end
unitsigmoid(T::Tuple, args...) = unitsigmoid(T[1], args...) # For compatibility with other parameter profile syntax. If T is tuple, use T[1] as the threshold
export unitsigmoid

# Same as unitStep, but come back down
function unitBump(T::Tuple, baseline::Real=0.0, bumpHeight::Real=1.0, stepOpt::Real=1.0)
    D = unitStep(T[1], baseline, bumpHeight, stepOpt) - unitStep(T[2], 0.0, bumpHeight, stepOpt)
end
export unitBump

function compositeStep(stepThrs::Vector, baselines::Vector, stepHeights::Vector)
    D = map(unitStep, stepThrs, baselines, stepHeights)
    D = sum(D)
end
compositeStep(T::Tuple, args...) = compositeStep(T[1], args...)
export compositeStep

"""
    stepNoise()

Compose a discontinuous function of many square bumps and troughs, following a white noise process.
    'stepHeight':   the standard deviation of the step height
    'T':            a tuple, (t0, tmax)

# Examples
```julia-repl
julia> D = stepNoise();
D([-1, 0, 1])
```
"""
function stepNoise(T::Tuple, stepWidth::Real=100, stepHeight::Real=1, baseline::Real=0)
    stepIdxs = T[1]:stepWidth:T[2]-stepWidth
    steps = stepHeight.*randn(Float64, (1, length(stepIdxs)))
    ps = sum(map((x, y) -> unitBump((x, x+stepWidth), baseline, y), stepIdxs, steps))
end
export stepNoise

function stepRandomWalk(T::Tuple, stepWidth::Real=100, stepHeight::Real=1, baseline::Real=0)
    stepIdxs = T[1]:stepWidth:T[2]
    steps = stepHeight.*randn(Float64, (1, length(stepIdxs)))
    ps = sum(map((x, y) -> unitStep((x, x+stepWidth), baseline, y), stepIdxs, steps))
end
export stepRandomWalk

function ramp(gradient::Real=1, p0::Real=0, t0::Real=0)
    x -> gradient.*(x.-t0) .+ p0
end
function ramp(p1::Real, p2::Real, t1::Real, t2::Real)
    gradient = (p2 - p1)/(t2 - t1)
    ramp(gradient, p1, t1)
end
export ramp

"""
Define a line over an interval. Saturate before and after this interval
"""
function rampInterval(p1::Real, p2::Real, t1::Real, t2::Real)
    r = ramp(p1, p2, t1, t2)
    b = unitBump((t1, t2), 0, 1)
    s1 = unitStep(t1, p1, -p1)
    s2 = unitStep(t2, 0, p2)
    return r*b + s1 + s2
end
export rampInterval

"""
Define a line over an interval. Before this interval, saturate, but after, extrapolate.
"""
function rampOn(p1::Real, p2::Real, t1::Real, t2::Real)
    r = ramp(p1, p2, t1, t2)
    b = unitStep(t1, 0, 1)
    s1 = unitStep(t1, p1, -p1)
    return r*b + s1
end
export rampOn

"""
Define a line over an interval. After this interval, saturate, but before, extrapolate.
"""
function rampOff(p1::Real, p2::Real, t1::Real, t2::Real)
    r = ramp(p1, p2, t1, t2)
    b = unitStep(t2, 1, -1)
    s2 = unitStep(t2, 0, p2)
    return r*b + s2
end
export rampOff

function sineWave(period::Real=1, amplitude::Real=1, t0::Real=0, baseline::Real=5*amplitude)
    x -> amplitude.*sin.((2π/period).*(x.-t0)) + baseline
end
export sineWave

function triangleWave(period::Real=1, amplitude::Real=1, t0::Real=0, tmax=100, baseline::Real=5*amplitude)
    # This is a discontinuous wave, so we can't (yet) define it over all x
    # Be careful to set the t0 and tmax to match the simulation time
    f = x -> (2amplitude/π).*asin.(sin.((2π/period).*(x.-t0))) + baseline
    d = Set((t0 + π/2):π:tmax) # The set of maxima and minima of the triangle wave
    Discontinuous(f, d)
end
export triangleWave

function lorentzian(A=1.0, Γ=1.0, x₀=0.0, y₀=0.0)
    # Γ is F.W.H.M, x₀ is centre, A is height and y₀ is the baseline
    # (function is not mass normalised, as would be usual)
    γ = Γ/2
    x -> (A.*γ.^2)./((x.-x₀).^2 .+ γ.^2) + y₀
end
export lorentzian


text:
module Integration

import SparseArrays: sparse

import ..Elements: ReferenceElement, quadrule
import ..Functions: callable, CompiledArray, CompiledSparseArray
import ..Topologies: Topology, refelems

export integrate

function integrate(func::CompiledArray{T}, domain::Topology, npts::Int) where T
    quadrules = Dict{ReferenceElement, Tuple{Array{Float64,2}, Vector{Float64}}}(
        elem => quadrule(elem, npts) for elem in refelems(domain)
    )

    result = zeros(T, size(func)...)
    kernel = callable(func)
    for elem in domain
        (pts, wts) = quadrules[elem.reference]
        for i = 1:length(wts)
            result .+= kernel(pts[i,:], elem) .* wts[i]
        end
    end

    result
end

function integrate(func::CompiledSparseArray{T,1}, domain::Topology, npts::Int) where T
    quadrules = Dict{ReferenceElement, Tuple{Array{Float64,2}, Vector{Float64}}}(
        elem => quadrule(elem, npts) for elem in refelems(domain)
    )

    result = zeros(T, length(func))
    ikernel, dkernel = callable(func)
    for elem in domain
        (pts, wts) = quadrules[elem.reference]
        blockI, = ikernel(pts[1,:], elem) :: Tuple{Vector{Int}}
        for i = 1:length(wts)
            data = dkernel(pts[i,:], elem) :: Array{T,1}
            result[blockI] += data .* wts[i]
        end
    end

    result
end

function integrate(func::CompiledSparseArray{T,2}, domain::Topology, npts::Int) where T
    quadrules = Dict{ReferenceElement, Tuple{Array{Float64,2}, Vector{Float64}}}(
        elem => quadrule(elem, npts) for elem in refelems(domain)
    )

    nelems = length(domain)
    I = zeros(Int, func.blockshape..., nelems) :: Array{Int, 3}
    J = zeros(Int, func.blockshape..., nelems) :: Array{Int, 3}
    V = zeros(T, func.blockshape..., nelems) :: Array{T, 3}

    ikernel, dkernel = callable(func)
    for (elemid, elem) in enumerate(domain)
        (pts, wts) = quadrules[elem.reference]

        blockI, blockJ = ikernel(pts[1,:], elem) :: Tuple{Vector{Int}, Vector{Int}}
        I[:, :, elemid] .= blockI
        J[:, :, elemid] .= reshape(blockJ, 1, :)

        for i = 1:length(wts)
            data = dkernel(pts[i,:], elem) :: Array{T,2}
            V[:, :, elemid] += data .* wts[i]
        end
    end

    sparse(I[:], J[:], V[:], size(func)...)
end

end # module

max_stars_repo_name: TheBB/Jutils

max_issues_repo_name: TheBB/Jutils

max_forks_repo_name: TheBB/Jutils


text:
# julia make.jl
# mkdocs build
# mkdocs gh-deploy --clean
using Documenter, Couchzilla

makedocs(modules=[Couchzilla], doctest=true)

deploydocs(
	deps = Deps.pip("mkdocs"),
    repo = "github.com/xpqz/Couchzilla.jl.git",
    julia  = "0.4.6",
    osname = "osx")

max_stars_repo_name: xpqz/couchzilla

max_issues_repo_name: xpqz/couchzilla

max_forks_repo_name: xpqz/couchzilla


text:
using DiffEqJump, DiffEqBase
using Test

rate = (u,p,t) -> u[1]
affect! = function (integrator)
  integrator.u[1] += 1
end
jump = ConstantRateJump(rate,affect!)

rate = (u,p,t) -> 0.5u[1]
affect! = function (integrator)
  integrator.u[1] -= 1
end
jump2 = ConstantRateJump(rate,affect!)

prob = DiscreteProblem([10.0],(0.0,3.0))
jump_prob = JumpProblem(prob,Direct(),jump,jump2)

integrator = init(jump_prob,SSAStepper())
step!(integrator)
integrator.u[1]

# test different saving behaviors

sol = solve(jump_prob,SSAStepper())
@test sol.t[begin] == 0.0
@test sol.t[end] == 3.0

sol = solve(jump_prob, SSAStepper(),save_end=false)
@test sol.t[begin] == 0.0
@test sol.t[end] < 3.0

sol = solve(jump_prob, SSAStepper(),save_start=false)
@test sol.t[begin] > 0.0
@test sol.t[end] == 3.0

jump_prob = JumpProblem(prob,Direct(),jump,jump2,save_positions=(false,false))
sol = solve(jump_prob, SSAStepper(),save_start=false,save_end=false)
@test isempty(sol.t) && isempty(sol.u)

sol = solve(jump_prob,SSAStepper(),saveat=0.0:0.1:2.9)
@test sol.t == collect(0.0:0.1:3.0)


text:
using Zygote
using Zygote: @adjoint
using Yao, Yao.AD

@adjoint function apply!(reg::ArrayReg, block::AbstractBlock)
    out = apply!(reg, block)
    out, function (outδ)
        (in, inδ), paramsδ = apply_back((out, outδ), block)
        return (inδ, paramsδ)
    end
end

@adjoint function Yao.dispatch!(block::AbstractBlock, params)
    out = dispatch!(block, params)
    out, function (outδ)
        (nothing, outδ)
    end
end

@adjoint function Matrix(block::AbstractBlock)
    out = Matrix(block)
    out, function (outδ)
        paramsδ = mat_back(block, outδ)
        return (paramsδ,)
    end
end

@adjoint function ArrayReg{B}(raw::AbstractArray) where B
    ArrayReg{B}(raw), adjy->(reshape(adjy.state, size(raw)),)
end

@adjoint function ArrayReg{B}(raw::ArrayReg) where B
    ArrayReg{B}(raw), adjy->(adjy,)
end

@adjoint function ArrayReg(raw::AbstractArray)
    ArrayReg(raw), adjy->(reshape(adjy.state, size(raw)),)
end

@adjoint function copy(reg::ArrayReg) where B
    copy(reg), adjy->(adjy,)
end

@adjoint state(reg::ArrayReg) = state(reg), adjy->(ArrayReg(adjy),)
@adjoint statevec(reg::ArrayReg) = statevec(reg), adjy->(ArrayReg(adjy),)
@adjoint state(reg::AdjointArrayReg) = state(reg), adjy->(ArrayReg(adjy')',)
@adjoint statevec(reg::AdjointArrayReg) = statevec(reg), adjy->(ArrayReg(adjy')',)
@adjoint parent(reg::AdjointArrayReg) = parent(reg), adjy->(adjy',)
@adjoint Base.adjoint(reg::ArrayReg) = Base.adjoint(reg), adjy->(parent(adjy),)
Zygote.@nograd Yao.nparameters
Zygote.@nograd Yao.zero_state
Zygote.@nograd Yao.rand_state
Zygote.@nograd Yao.uniform_state
Zygote.@nograd Yao.product_state


text:
"""
    LinearCombination(coefficients) <: Transform

Calculates the linear combination of a collection of terms weighted by some `coefficients`.

When applied to an N-dimensional array, `LinearCombination` reduces along the `dim` provided
and returns an (N-1)-dimensional array.

If no `inds` are specified, then the transform is applied to all elements.

!!!note
    The current default is that `dims=1` but this behaviour will be deprecated in a future
    release and the `dims` keyword argument will have to be specified explicitly.
    https://github.com/invenia/FeatureTransforms.jl/issues/82
"""
struct LinearCombination <: Transform
    coefficients::Vector{Real}
end

cardinality(::LinearCombination) = ManyToOne()

function _apply(terms, LC::LinearCombination; kwargs...)
    # Need this check because map will work even if there are more/less terms than coeffs
    if length(terms) != length(LC.coefficients)
        throw(DimensionMismatch(
            "Number of terms $(length(terms)) does not match "*
            "number of coefficients $(length(LC.coefficients))."
        ))
    end

   return sum(map(*, terms, LC.coefficients))
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# Bidiagonal matrices
type Bidiagonal{T} <: AbstractMatrix{T}
    dv::Vector{T} # diagonal
    ev::Vector{T} # sub/super diagonal
    isupper::Bool # is upper bidiagonal (true) or lower (false)
    function Bidiagonal(dv::Vector{T}, ev::Vector{T}, isupper::Bool)
        if length(ev) != length(dv)-1
            throw(DimensionMismatch("length of diagonal vector is $(length(dv)), length of off-diagonal vector is $(length(ev))"))
        end
        new(dv, ev, isupper)
    end
end
"""
    Bidiagonal(dv, ev, isupper::Bool)

Constructs an upper (`isupper=true`) or lower (`isupper=false`) bidiagonal matrix using the
given diagonal (`dv`) and off-diagonal (`ev`) vectors.  The result is of type `Bidiagonal`
and provides efficient specialized linear solvers, but may be converted into a regular
matrix with [`convert(Array, _)`](@ref) (or `Array(_)` for short). `ev`'s length
must be one less than the length of `dv`.

# Example

```jldoctest
julia> dv = [1; 2; 3; 4]
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> ev = [7; 8; 9]
3-element Array{Int64,1}:
 7
 8
 9

julia> Bu = Bidiagonal(dv, ev, true) # ev is on the first superdiagonal
4×4 Bidiagonal{Int64}:
 1  7  ⋅  ⋅
 ⋅  2  8  ⋅
 ⋅  ⋅  3  9
 ⋅  ⋅  ⋅  4

julia> Bl = Bidiagonal(dv, ev, false) # ev is on the first subdiagonal
4×4 Bidiagonal{Int64}:
 1  ⋅  ⋅  ⋅
 7  2  ⋅  ⋅
 ⋅  8  3  ⋅
 ⋅  ⋅  9  4
```
"""
Bidiagonal{T}(dv::AbstractVector{T}, ev::AbstractVector{T}, isupper::Bool) = Bidiagonal{T}(collect(dv), collect(ev), isupper)
Bidiagonal(dv::AbstractVector, ev::AbstractVector) = throw(ArgumentError("did you want an upper or lower Bidiagonal? Try again with an additional true (upper) or false (lower) argument."))

"""
    Bidiagonal(dv, ev, uplo::Char)

Constructs an upper (`uplo='U'`) or lower (`uplo='L'`) bidiagonal matrix using the
given diagonal (`dv`) and off-diagonal (`ev`) vectors.  The result is of type `Bidiagonal`
and provides efficient specialized linear solvers, but may be converted into a regular
matrix with [`convert(Array, _)`](@ref) (or `Array(_)` for short). `ev`'s
length must be one less than the length of `dv`.

# Example

```jldoctest
julia> dv = [1; 2; 3; 4]
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> ev = [7; 8; 9]
3-element Array{Int64,1}:
 7
 8
 9

julia> Bu = Bidiagonal(dv, ev, 'U') #e is on the first superdiagonal
4×4 Bidiagonal{Int64}:
 1  7  ⋅  ⋅
 ⋅  2  8  ⋅
 ⋅  ⋅  3  9
 ⋅  ⋅  ⋅  4

julia> Bl = Bidiagonal(dv, ev, 'L') #e is on the first subdiagonal
4×4 Bidiagonal{Int64}:
 1  ⋅  ⋅  ⋅
 7  2  ⋅  ⋅
 ⋅  8  3  ⋅
 ⋅  ⋅  9  4
```
"""
#Convert from BLAS uplo flag to boolean internal
Bidiagonal(dv::AbstractVector, ev::AbstractVector, uplo::Char) = begin
    if uplo === 'U'
        isupper = true
    elseif uplo === 'L'
        isupper = false
    else
        throw(ArgumentError("Bidiagonal uplo argument must be upper 'U' or lower 'L', got $(repr(uplo))"))
    end
    Bidiagonal(collect(dv), collect(ev), isupper)
end
function Bidiagonal{Td,Te}(dv::AbstractVector{Td}, ev::AbstractVector{Te}, isupper::Bool)
    T = promote_type(Td,Te)
    Bidiagonal(convert(Vector{T}, dv), convert(Vector{T}, ev), isupper)
end

"""
    Bidiagonal(A, isupper::Bool)

Construct a `Bidiagonal` matrix from the main diagonal of `A` and
its first super- (if `isupper=true`) or sub-diagonal (if `isupper=false`).

# Example

```jldoctest
julia> A = [1 1 1 1; 2 2 2 2; 3 3 3 3; 4 4 4 4]
4×4 Array{Int64,2}:
 1  1  1  1
 2  2  2  2
 3  3  3  3
 4  4  4  4

julia> Bidiagonal(A, true) #contains the main diagonal and first superdiagonal of A
4×4 Bidiagonal{Int64}:
 1  1  ⋅  ⋅
 ⋅  2  2  ⋅
 ⋅  ⋅  3  3
 ⋅  ⋅  ⋅  4

julia> Bidiagonal(A, false) #contains the main diagonal and first subdiagonal of A
4×4 Bidiagonal{Int64}:
 1  ⋅  ⋅  ⋅
 2  2  ⋅  ⋅
 ⋅  3  3  ⋅
 ⋅  ⋅  4  4
```
"""
Bidiagonal(A::AbstractMatrix, isupper::Bool)=Bidiagonal(diag(A), diag(A, isupper?1:-1), isupper)

function getindex{T}(A::Bidiagonal{T}, i::Integer, j::Integer)
    if !((1 <= i <= size(A,2)) && (1 <= j <= size(A,2)))
        throw(BoundsError(A,(i,j)))
    end
    if i == j
        return A.dv[i]
    elseif (istriu(A) && (i == j - 1)) || (istril(A) && (i == j + 1))
        return A.ev[min(i,j)]
    else
        return zero(T)
    end
end

function setindex!(A::Bidiagonal, x, i::Integer, j::Integer)
    if i == j
        A.dv[i] = x
    elseif (istriu(A) && (i == j - 1)) || (istril(A) && (i == j + 1))
        return A.ev[min(i,j)] = x
    else
        throw(ArgumentError("cannot set elements outside main and $(istriu(A) ? "super": "sub") diagonals."))
    end
end

## structured matrix methods ##
function Base.replace_in_print_matrix(A::Bidiagonal,i::Integer,j::Integer,s::AbstractString)
    if A.isupper
        i==j || i==j-1 ? s : Base.replace_with_centered_mark(s)
    else
        i==j || i==j+1 ? s : Base.replace_with_centered_mark(s)
    end
end

#Converting from Bidiagonal to dense Matrix
function convert{T}(::Type{Matrix{T}}, A::Bidiagonal)
    n = size(A, 1)
    B = zeros(T, n, n)
    for i = 1:n - 1
        B[i,i] = A.dv[i]
        if A.isupper
            B[i, i + 1] = A.ev[i]
        else
            B[i + 1, i] = A.ev[i]
        end
    end
    B[n,n] = A.dv[n]
    return B
end
convert{T}(::Type{Matrix}, A::Bidiagonal{T}) = convert(Matrix{T}, A)
convert(::Type{Array}, A::Bidiagonal) = convert(Matrix, A)
full(A::Bidiagonal) = convert(Array, A)
promote_rule{T,S}(::Type{Matrix{T}}, ::Type{Bidiagonal{S}})=Matrix{promote_type(T,S)}

#Converting from Bidiagonal to Tridiagonal
Tridiagonal{T}(M::Bidiagonal{T}) = convert(Tridiagonal{T}, M)
function convert{T}(::Type{Tridiagonal{T}}, A::Bidiagonal)
    z = zeros(T, size(A)[1]-1)
    A.isupper ? Tridiagonal(z, convert(Vector{T},A.dv), convert(Vector{T},A.ev)) : Tridiagonal(convert(Vector{T},A.ev), convert(Vector{T},A.dv), z)
end
promote_rule{T,S}(::Type{Tridiagonal{T}}, ::Type{Bidiagonal{S}})=Tridiagonal{promote_type(T,S)}

# No-op for trivial conversion Bidiagonal{T} -> Bidiagonal{T}
convert{T}(::Type{Bidiagonal{T}}, A::Bidiagonal{T}) = A
# Convert Bidiagonal{Told} to Bidiagonal{Tnew} by constructing a new instance with converted elements
convert{Tnew,Told}(::Type{Bidiagonal{Tnew}}, A::Bidiagonal{Told}) = Bidiagonal(convert(Vector{Tnew}, A.dv), convert(Vector{Tnew}, A.ev), A.isupper)
# When asked to convert Bidiagonal{Told} to AbstractMatrix{Tnew}, preserve structure by converting to Bidiagonal{Tnew} <: AbstractMatrix{Tnew}
convert{Tnew,Told}(::Type{AbstractMatrix{Tnew}}, A::Bidiagonal{Told}) = convert(Bidiagonal{Tnew}, A)

broadcast(::typeof(big), B::Bidiagonal) = Bidiagonal(big.(B.dv), big.(B.ev), B.isupper)

similar{T}(B::Bidiagonal, ::Type{T}) = Bidiagonal{T}(similar(B.dv, T), similar(B.ev, T), B.isupper)

###################
# LAPACK routines #
###################

#Singular values
svdvals!{T<:BlasReal}(M::Bidiagonal{T}) = LAPACK.bdsdc!(M.isupper ? 'U' : 'L', 'N', M.dv, M.ev)[1]
function svdfact!{T<:BlasReal}(M::Bidiagonal{T}; thin::Bool=true)
    d, e, U, Vt, Q, iQ = LAPACK.bdsdc!(M.isupper ? 'U' : 'L', 'I', M.dv, M.ev)
    SVD(U, d, Vt)
end
svdfact(M::Bidiagonal; thin::Bool=true) = svdfact!(copy(M),thin=thin)

####################
# Generic routines #
####################

function show(io::IO, M::Bidiagonal)
    # TODO: make this readable and one-line
    println(io, summary(M), ":")
    print(io, " diag:")
    print_matrix(io, (M.dv)')
    print(io, M.isupper?"\n super:":"\n sub:")
    print_matrix(io, (M.ev)')
end

size(M::Bidiagonal) = (length(M.dv), length(M.dv))
function size(M::Bidiagonal, d::Integer)
    if d < 1
        throw(ArgumentError("dimension must be ≥ 1, got $d"))
    elseif d <= 2
        return length(M.dv)
    else
        return 1
    end
end

#Elementary operations
broadcast(::typeof(abs), M::Bidiagonal) = Bidiagonal(abs.(M.dv), abs.(M.ev), abs.(M.isupper))
broadcast(::typeof(round), M::Bidiagonal) = Bidiagonal(round.(M.dv), round.(M.ev), M.isupper)
broadcast(::typeof(trunc), M::Bidiagonal) = Bidiagonal(trunc.(M.dv), trunc.(M.ev), M.isupper)
broadcast(::typeof(floor), M::Bidiagonal) = Bidiagonal(floor.(M.dv), floor.(M.ev), M.isupper)
broadcast(::typeof(ceil), M::Bidiagonal) = Bidiagonal(ceil.(M.dv), ceil.(M.ev), M.isupper)
for func in (:conj, :copy, :real, :imag)
    @eval ($func)(M::Bidiagonal) = Bidiagonal(($func)(M.dv), ($func)(M.ev), M.isupper)
end
broadcast{T<:Integer}(::typeof(round), ::Type{T}, M::Bidiagonal) = Bidiagonal(round.(T, M.dv), round.(T, M.ev), M.isupper)
broadcast{T<:Integer}(::typeof(trunc), ::Type{T}, M::Bidiagonal) = Bidiagonal(trunc.(T, M.dv), trunc.(T, M.ev), M.isupper)
broadcast{T<:Integer}(::typeof(floor), ::Type{T}, M::Bidiagonal) = Bidiagonal(floor.(T, M.dv), floor.(T, M.ev), M.isupper)
broadcast{T<:Integer}(::typeof(ceil), ::Type{T}, M::Bidiagonal) = Bidiagonal(ceil.(T, M.dv), ceil.(T, M.ev), M.isupper)

transpose(M::Bidiagonal) = Bidiagonal(M.dv, M.ev, !M.isupper)
ctranspose(M::Bidiagonal) = Bidiagonal(conj(M.dv), conj(M.ev), !M.isupper)

istriu(M::Bidiagonal) = M.isupper || iszero(M.ev)
istril(M::Bidiagonal) = !M.isupper || iszero(M.ev)

function tril!(M::Bidiagonal, k::Integer=0)
    n = length(M.dv)
    if abs(k) > n
        throw(ArgumentError("requested diagonal, $k, out of bounds in matrix of size ($n,$n)"))
    elseif M.isupper && k < 0
        fill!(M.dv,0)
        fill!(M.ev,0)
    elseif k < -1
        fill!(M.dv,0)
        fill!(M.ev,0)
    elseif M.isupper && k == 0
        fill!(M.ev,0)
    elseif !M.isupper && k == -1
        fill!(M.dv,0)
    end
    return M
end

function triu!(M::Bidiagonal, k::Integer=0)
    n = length(M.dv)
    if abs(k) > n
        throw(ArgumentError("requested diagonal, $k, out of bounds in matrix of size ($n,$n)"))
    elseif !M.isupper && k > 0
        fill!(M.dv,0)
        fill!(M.ev,0)
    elseif k > 1
        fill!(M.dv,0)
        fill!(M.ev,0)
    elseif !M.isupper && k == 0
        fill!(M.ev,0)
    elseif M.isupper && k == 1
        fill!(M.dv,0)
    end
    return M
end

function diag{T}(M::Bidiagonal{T}, n::Integer=0)
    if n == 0
        return M.dv
    elseif n == 1
        return M.isupper ? M.ev : zeros(T, size(M,1)-1)
    elseif n == -1
        return M.isupper ? zeros(T, size(M,1)-1) : M.ev
    elseif -size(M,1) < n < size(M,1)
        return zeros(T, size(M,1)-abs(n))
    else
        throw(ArgumentError("matrix size is $(size(M)), n is $n"))
    end
end

function +(A::Bidiagonal, B::Bidiagonal)
    if A.isupper==B.isupper
        Bidiagonal(A.dv+B.dv, A.ev+B.ev, A.isupper)
    else
        Tridiagonal((A.isupper ? (B.ev,A.dv+B.dv,A.ev) : (A.ev,A.dv+B.dv,B.ev))...)
    end
end

function -(A::Bidiagonal, B::Bidiagonal)
    if A.isupper==B.isupper
        Bidiagonal(A.dv-B.dv, A.ev-B.ev, A.isupper)
    else
        Tridiagonal((A.isupper ? (-B.ev,A.dv-B.dv,A.ev) : (A.ev,A.dv-B.dv,-B.ev))...)
    end
end

-(A::Bidiagonal)=Bidiagonal(-A.dv,-A.ev,A.isupper)
*(A::Bidiagonal, B::Number) = Bidiagonal(A.dv*B, A.ev*B, A.isupper)
*(B::Number, A::Bidiagonal) = A*B
/(A::Bidiagonal, B::Number) = Bidiagonal(A.dv/B, A.ev/B, A.isupper)
==(A::Bidiagonal, B::Bidiagonal) = (A.dv==B.dv) && (A.ev==B.ev) && (A.isupper==B.isupper)

BiTriSym = Union{Bidiagonal, Tridiagonal, SymTridiagonal}
BiTri = Union{Bidiagonal, Tridiagonal}
A_mul_B!(C::AbstractMatrix, A::SymTridiagonal, B::BiTriSym) = A_mul_B_td!(C, A, B)
A_mul_B!(C::AbstractMatrix, A::BiTri, B::BiTriSym) = A_mul_B_td!(C, A, B)
A_mul_B!(C::AbstractMatrix, A::BiTriSym, B::BiTriSym) = A_mul_B_td!(C, A, B)
A_mul_B!(C::AbstractMatrix, A::AbstractTriangular, B::BiTriSym) = A_mul_B_td!(C, A, B)
A_mul_B!(C::AbstractMatrix, A::AbstractMatrix, B::BiTriSym) = A_mul_B_td!(C, A, B)
A_mul_B!(C::AbstractVector, A::BiTri, B::AbstractVector) = A_mul_B_td!(C, A, B)
A_mul_B!(C::AbstractMatrix, A::BiTri, B::AbstractVecOrMat) = A_mul_B_td!(C, A, B)
A_mul_B!(C::AbstractVecOrMat, A::BiTri, B::AbstractVecOrMat) = A_mul_B_td!(C, A, B)

\(::Diagonal, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector"))
\(::Bidiagonal, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector"))
\{TA<:Number,TB<:Number}(::Bidiagonal{TA}, ::RowVector{TB}) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector"))

At_ldiv_B(::Bidiagonal, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector"))
At_ldiv_B{TA<:Number,TB<:Number}(::Bidiagonal{TA}, ::RowVector{TB}) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector"))

Ac_ldiv_B(::Bidiagonal, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector"))
Ac_ldiv_B{TA<:Number,TB<:Number}(::Bidiagonal{TA}, ::RowVector{TB}) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector"))

function check_A_mul_B!_sizes(C, A, B)
    nA, mA = size(A)
    nB, mB = size(B)
    nC, mC = size(C)
    if !(nA == nC)
        throw(DimensionMismatch("sizes size(A)=$(size(A)) and size(C) = $(size(C)) must match at first entry."))
    elseif !(mA == nB)
        throw(DimensionMismatch("second entry of size(A)=$(size(A)) and first entry of size(B) = $(size(B)) must match."))
    elseif !(mB == mC)
        throw(DimensionMismatch("sizes size(B)=$(size(B)) and size(C) = $(size(C)) must match at first second entry."))
    end
end

function A_mul_B_td!(C::AbstractMatrix, A::BiTriSym, B::BiTriSym)
    check_A_mul_B!_sizes(C, A, B)
    n = size(A,1)
    n <= 3 && return A_mul_B!(C, Array(A), Array(B))
    fill!(C, zero(eltype(C)))
    Al = diag(A, -1)
    Ad = diag(A, 0)
    Au = diag(A, 1)
    Bl = diag(B, -1)
    Bd = diag(B, 0)
    Bu = diag(B, 1)
    @inbounds begin
        # first row of C
        C[1,1] = A[1,1]*B[1,1] + A[1, 2]*B[2, 1]
        C[1,2] = A[1,1]*B[1,2] + A[1,2]*B[2,2]
        C[1,3] = A[1,2]*B[2,3]
        # second row of C
        C[2,1] = A[2,1]*B[1,1] + A[2,2]*B[2,1]
        C[2,2] = A[2,1]*B[1,2] + A[2,2]*B[2,2] + A[2,3]*B[3,2]
        C[2,3] = A[2,2]*B[2,3] + A[2,3]*B[3,3]
        C[2,4] = A[2,3]*B[3,4]
        for j in 3:n-2
            Ajj₋1   = Al[j-1]
            Ajj     = Ad[j]
            Ajj₊1   = Au[j]
            Bj₋1j₋2 = Bl[j-2]
            Bj₋1j₋1 = Bd[j-1]
            Bj₋1j   = Bu[j-1]
            Bjj₋1   = Bl[j-1]
            Bjj     = Bd[j]
            Bjj₊1   = Bu[j]
            Bj₊1j   = Bl[j]
            Bj₊1j₊1 = Bd[j+1]
            Bj₊1j₊2 = Bu[j+1]
            C[j,j-2]  = Ajj₋1*Bj₋1j₋2
            C[j, j-1] = Ajj₋1*Bj₋1j₋1 + Ajj*Bjj₋1
            C[j, j  ] = Ajj₋1*Bj₋1j   + Ajj*Bjj       + Ajj₊1*Bj₊1j
            C[j, j+1] = Ajj  *Bjj₊1   + Ajj₊1*Bj₊1j₊1
            C[j, j+2] = Ajj₊1*Bj₊1j₊2
        end
        # row before last of C
        C[n-1,n-3] = A[n-1,n-2]*B[n-2,n-3]
        C[n-1,n-2] = A[n-1,n-1]*B[n-1,n-2] + A[n-1,n-2]*B[n-2,n-2]
        C[n-1,n-1] = A[n-1,n-2]*B[n-2,n-1] + A[n-1,n-1]*B[n-1,n-1] + A[n-1,n]*B[n,n-1]
        C[n-1,n  ] = A[n-1,n-1]*B[n-1,n  ] + A[n-1,  n]*B[n  ,n  ]
        # last row of C
        C[n,n-2] = A[n,n-1]*B[n-1,n-2]
        C[n,n-1] = A[n,n-1]*B[n-1,n-1] + A[n,n]*B[n,n-1]
        C[n,n  ] = A[n,n-1]*B[n-1,n  ] + A[n,n]*B[n,n  ]
    end # inbounds
    C
end

function A_mul_B_td!(C::AbstractVecOrMat, A::BiTriSym, B::AbstractVecOrMat)
    nA = size(A,1)
    nB = size(B,2)
    if !(size(C,1) == size(B,1) == nA)
        throw(DimensionMismatch("A has first dimension $nA, B has $(size(B,1)), C has $(size(C,1)) but all must match"))
    end
    if size(C,2) != nB
        throw(DimensionMismatch("A has second dimension $nA, B has $(size(B,2)), C has $(size(C,2)) but all must match"))
    end
    nA <= 3 && return A_mul_B!(C, Array(A), Array(B))
    l = diag(A, -1)
    d = diag(A, 0)
    u = diag(A, 1)
    @inbounds begin
        for j = 1:nB
            b₀, b₊ = B[1, j], B[2, j]
            C[1, j] = d[1]*b₀ + u[1]*b₊
            for i = 2:nA - 1
                b₋, b₀, b₊ = b₀, b₊, B[i + 1, j]
                C[i, j] = l[i - 1]*b₋ + d[i]*b₀ + u[i]*b₊
            end
            C[nA, j] = l[nA - 1]*b₀ + d[nA]*b₊
        end
    end
    C
end

function A_mul_B_td!(C::AbstractMatrix, A::AbstractMatrix, B::BiTriSym)
    check_A_mul_B!_sizes(C, A, B)
    n = size(A,1)
    n <= 3 && return A_mul_B!(C, Array(A), Array(B))
    m = size(B,2)
    Bl = diag(B, -1)
    Bd = diag(B, 0)
    Bu = diag(B, 1)
    @inbounds begin
        # first and last column of C
        B11 = Bd[1]
        B21 = Bl[1]
        Bmm = Bd[m]
        Bm₋1m = Bu[m-1]
        for i in 1:n
            C[i, 1] = A[i,1] * B11 + A[i, 2] * B21
            C[i, m] = A[i, m-1] * Bm₋1m + A[i, m] * Bmm
        end
        # middle columns of C
        for j = 2:m-1
            Bj₋1j = Bu[j-1]
            Bjj = Bd[j]
            Bj₊1j = Bl[j]
            for i = 1:n
                C[i, j] = A[i, j-1] * Bj₋1j + A[i, j]*Bjj + A[i, j+1] * Bj₊1j
            end
        end
    end # inbounds
    C
end

SpecialMatrix = Union{Bidiagonal, SymTridiagonal, Tridiagonal}
# to avoid ambiguity warning, but shouldn't be necessary
*(A::AbstractTriangular, B::SpecialMatrix) = Array(A) * Array(B)
*(A::SpecialMatrix, B::SpecialMatrix) = Array(A) * Array(B)

#Generic multiplication
for func in (:*, :Ac_mul_B, :A_mul_Bc, :/, :A_rdiv_Bc)
    @eval ($func){T}(A::Bidiagonal{T}, B::AbstractVector{T}) = ($func)(Array(A), B)
end

#Linear solvers
A_ldiv_B!(A::Union{Bidiagonal, AbstractTriangular}, b::AbstractVector) = naivesub!(A, b)
At_ldiv_B!(A::Bidiagonal, b::AbstractVector) = A_ldiv_B!(transpose(A), b)
Ac_ldiv_B!(A::Bidiagonal, b::AbstractVector) = A_ldiv_B!(ctranspose(A), b)
function A_ldiv_B!(A::Union{Bidiagonal, AbstractTriangular}, B::AbstractMatrix)
    nA,mA = size(A)
    tmp = similar(B,size(B,1))
    n = size(B, 1)
    if nA != n
        throw(DimensionMismatch("size of A is ($nA,$mA), corresponding dimension of B is $n"))
    end
    for i = 1:size(B,2)
        copy!(tmp, 1, B, (i - 1)*n + 1, n)
        A_ldiv_B!(A, tmp)
        copy!(B, (i - 1)*n + 1, tmp, 1, n) # Modify this when array view are implemented.
    end
    B
end
for func in (:Ac_ldiv_B!, :At_ldiv_B!)
    @eval function ($func)(A::Union{Bidiagonal, AbstractTriangular}, B::AbstractMatrix)
        nA,mA = size(A)
        tmp = similar(B,size(B,1))
        n = size(B, 1)
        if mA != n
            throw(DimensionMismatch("size of A' is ($mA,$nA), corresponding dimension of B is $n"))
        end
        for i = 1:size(B,2)
            copy!(tmp, 1, B, (i - 1)*n + 1, n)
            ($func)(A, tmp)
            copy!(B, (i - 1)*n + 1, tmp, 1, n) # Modify this when array view are implemented.
        end
        B
    end
end
#Generic solver using naive substitution
function naivesub!{T}(A::Bidiagonal{T}, b::AbstractVector, x::AbstractVector = b)
    N = size(A, 2)
    if N != length(b) || N != length(x)
        throw(DimensionMismatch("second dimension of A, $N, does not match one of the lengths of x, $(length(x)), or b, $(length(b))"))
    end
    if !A.isupper #do forward substitution
        for j = 1:N
            x[j] = b[j]
            j > 1 && (x[j] -= A.ev[j-1] * x[j-1])
            x[j] /= A.dv[j] == zero(T) ? throw(SingularException(j)) : A.dv[j]
        end
    else #do backward substitution
        for j = N:-1:1
            x[j] = b[j]
            j < N && (x[j] -= A.ev[j] * x[j+1])
            x[j] /= A.dv[j] == zero(T) ? throw(SingularException(j)) : A.dv[j]
        end
    end
    x
end

### Generic promotion methods and fallbacks
for (f,g) in ((:\, :A_ldiv_B!), (:At_ldiv_B, :At_ldiv_B!), (:Ac_ldiv_B, :Ac_ldiv_B!))
    @eval begin
        function ($f){TA<:Number,TB<:Number}(A::Bidiagonal{TA}, B::AbstractVecOrMat{TB})
            TAB = typeof((zero(TA)*zero(TB) + zero(TA)*zero(TB))/one(TA))
            ($g)(convert(AbstractArray{TAB}, A), copy_oftype(B, TAB))
        end
        ($f)(A::Bidiagonal, B::AbstractVecOrMat) = ($g)(A, copy(B))
    end
end

factorize(A::Bidiagonal) = A

# Eigensystems
eigvals(M::Bidiagonal) = M.dv
function eigvecs{T}(M::Bidiagonal{T})
    n = length(M.dv)
    Q = Array{T}(n, n)
    blks = [0; find(x -> x == 0, M.ev); n]
    v = zeros(T, n)
    if M.isupper
        for idx_block = 1:length(blks) - 1, i = blks[idx_block] + 1:blks[idx_block + 1] #index of eigenvector
            fill!(v, zero(T))
            v[blks[idx_block] + 1] = one(T)
            for j = blks[idx_block] + 1:i - 1 #Starting from j=i, eigenvector elements will be 0
                v[j+1] = (M.dv[i] - M.dv[j])/M.ev[j] * v[j]
            end
            c = norm(v)
            for j = 1:n
                Q[j, i] = v[j] / c
            end
        end
    else
        for idx_block = 1:length(blks) - 1, i = blks[idx_block + 1]:-1:blks[idx_block] + 1 #index of eigenvector
            fill!(v, zero(T))
            v[blks[idx_block+1]] = one(T)
            for j = (blks[idx_block+1] - 1):-1:max(1, (i - 1)) #Starting from j=i, eigenvector elements will be 0
                v[j] = (M.dv[i] - M.dv[j+1])/M.ev[j] * v[j+1]
            end
            c = norm(v)
            for j = 1:n
                Q[j, i] = v[j] / c
            end
        end
    end
    Q #Actually Triangular
end
eigfact(M::Bidiagonal) = Eigen(eigvals(M), eigvecs(M))

# fill! methods
_valuefields{T <: Diagonal}(S::Type{T}) = [:diag]
_valuefields{T <: Bidiagonal}(S::Type{T}) = [:dv, :ev]
_valuefields{T <: Tridiagonal}(S::Type{T}) = [:dl, :d, :du]
_valuefields{T <: SymTridiagonal}(S::Type{T}) = [:dv, :ev]
_valuefields{T <: AbstractTriangular}(S::Type{T}) = [:data]

SpecialArrays = Union{Diagonal,
    Bidiagonal,
    Tridiagonal,
    SymTridiagonal,
    AbstractTriangular}

@generated function fillslots!(A::SpecialArrays, x)
    ex = :(xT = convert(eltype(A), x))
    for field in _valuefields(A)
        ex = :($ex; fill!(A.$field, xT))
    end
    :($ex;return A)
end

# for historical reasons:
fill!(a::AbstractTriangular, x) = fillslots!(a, x);
fill!(D::Diagonal, x) = fillslots!(D, x);

_small_enough(A::Bidiagonal) = size(A, 1) <= 1
_small_enough(A::Tridiagonal) = size(A, 1) <= 2
_small_enough(A::SymTridiagonal) = size(A, 1) <= 2

function fill!(A::Union{Bidiagonal, Tridiagonal, SymTridiagonal} ,x)
    xT = convert(eltype(A), x)
    (xT == zero(eltype(A)) || _small_enough(A)) && return fillslots!(A, xT)
    throw(ArgumentError("array A of type $(typeof(A)) and size $(size(A)) can
    not be filled with x=$x, since some of its entries are constrained."))
end

max_stars_repo_name: oxinabox/julia

max_issues_repo_name: oxinabox/julia

max_forks_repo_name: oxinabox/julia